key: cord-0060990-08ilid30
authors: Milevsky, Moshe Arye
title: The Lifetime Ruin Probability (LRP)
date: 2020-09-27
journal: Retirement Income Recipes in R
DOI: 10.1007/978-3-030-51434-1_9
sha: 0c5fdc73ba0f163109c83da36455b6dd1d0b0c6d
doc_id: 60990
cord_uid: 08ilid30

This chapter returns to the realm of portfolio longevity and focuses on computational algorithms for success and failure rates associated with various retirement income strategies, but accounting for longevity risk. The chapter begins by defining the so-called lifetime ruin probability (LRP), which is the simplest retirement risk metric, widely used by practitioners. After reviewing the underlying probability concepts, the chapter provides a number of analytic expressions and simulation-based recipes for computing, interpreting, and understanding the limitations of the LRP.

• G(a,c) computes the incomplete Gamma function. • a(nu,x,m,b) computes a Gompertz random lifetime present value.

• LRPG(v,xi,x,m,b) computes the lifetime ruin probability (LRP) at age x, under a valuation rate v, initial withdrawal rate ξ (xi), and Gompertz (m, b) lifetime.

• VARPHI.SM(N,x,m,b,xi,nu,sigma), using simulation, approximates the lifetime ruin probability when investment returns are random with parameters (ν, σ ), under a Gompertz (m, b) model. • VARPHI.MM(x,m,b,xi,nu,sigma), same as VARPHI.SM, using moment-matching techniques.

Recall from the discussion in Chap. 5, that if you have an investment portfolio (i.e. financial capital) worth F , earning a fixed interest rate v, and are withdrawing a constant c per year, the longevity of your portfolio satisfies the following equation:

, as long as: v < ξ, (9.1) where ln [.] denotes the natural logarithm, the (new) symbol ξ = c/F is the (previously introduced) initial withdrawal rate, and: L ξ := ∞, when v ≥ ξ . This familiar equation was coded-up and stored in R as the: PL(v,c,F) function. So, for example, if the initial withdrawal rate is (the famous) ξ = 4%, and the real rate of interest is v = 1%, the portfolio longevity is PL(0.01,0.04,1)=28.77 years. Notice that I used F = 1 and c = 0.04, because portfolio longevity scales in c/F . So, is 28.77 years enough time? In Chap. 5 I alluded to a convention of 30 years, but didn't provide much justification, which was rather ad hoc and arbitrary. By this point, though, after spending (pun intended) two chapters discussing models of human longevity, we are in a better position to address lifetimes head-on. The pertinent question is: If indeed the longevity of your portfolio is 28.77 years, will you still be alive at that age? The answer obviously depends on the exact age at which you begin those withdrawals. If you start depleting your nest egg at the age of x = 70, then 28.77 years of longevity is likely enough. But if you "retire" early, at the age of x = 55, for example, it's unlikely the money will last for the rest of your life. This is because there is a high probability you will live to age y = 83.8, which is when the money runs out. This is the question motivating the entire chapter.

Will you be alive when the money runs out?

The answer isn't a binary yes or no. Rather, a proper response should be expressed as a probability of living to the longevity of your portfolio. So, if the portfolio's longevity is 28.77 years, I can compute the probability of living to that point, using explicit expressions for Pr[T x ≥ t], which was the focus of the prior chapter. I'll start with a (very) simple example and then move on to more general statements and formal definitions. Assume for the moment that your remaining lifetime random variable T x is exponentially distributed with a constant mortality hazard rate λ = 1/20 = 5%. And, although I dismissed constant mortality rates as being inconsistent with real-world data (unless you are a lobster), I'll leverage the analytic simplicity of this assumption to extract some financial insights. In a later section I'll return to more realistic models of human longevity, such as the Gompertz assumption. Anyway, assuming exponential lifetimes, the probability of surviving to any given t can be written explicitly as Pr[T x ≥ t] = e −λ t , where the subscript x is redundant but retained for consistency. Combining the survival probability with Eq. (9.1), we arrive at the formal definition of the lifetime ruin probability, for exponential lifetimes, denoted by the Greek letter varphi: (9.2) with the same restrictions imposed on portfolio longevity in Eq. (9.1), namely that v < ξ. Of course, in the event v ≥ ξ , which means that L ξ is defined to be L ξ = ∞, then (obviously) ϕ = 0. You can never die before the money runs out. Here is the interpretation. When the mortality rate λ = 1%, which is quite low (and means you will live longer), the probability of reaching 28.7 years (the longevity of the portfolio) is e −0.2877 , which works out to 75%. This is (very) high and probably unacceptable to most retirees. But, when the mortality rate is increased (in the script) to λ = 5% and then λ = 10%, the probability of living to the end of the life of the portfolio drops to 23.7% and a mere 5.6%, which is much lower (and better) than 75%. Now, few people control their mortality rate (at least to that extent) so this obviously isn't a strategy for reducing failure rates. Rather, this illustrates how mortality rates-and the associated human life expectancy-affect the sustainability of the retirement income plan. How's that? Well, when lifetimes T x are exponentially distributed with a mortality rate parameter of λ, the value of E[T x ] = 1/λ, which is the mean (human) lifetime, and SD[T x ] = 1/λ as well. So, another way to interpret the number in the above script box is as follows: If you are at an age with a human life expectancy of 100 = 1/(0.01) years, there is a 75% probability you will outlive your money. But, if your life expectancy is a mere 10 = 1/(0.1) years, with a standard deviation of 10 years, there is only a 5.6% probability of that (undesirable) event occurring. Of course, all of the above (simple) math and formulas assume that lifetimes are exponentially distributed (they are not) and that your financial capital is allocated to an asset earning a constant risk-free rate of v (it rarely is). Thus, to increase the level of realism (in the next section) I'll graduate to a Gompertz law of mortality and show how it affects results. Then, (in the section after that) I'll move on to random portfolio returns, per the methodology that was described in Chap. 5.

The lifetime ruin probability, ϕ, is defined in the same way as before, namely ϕ = Pr[T x ≥ L ξ ]. However, under the Gompertz law the survival probability involves both parameters (m, b), which are the mode and dispersion coefficients. An explicit expression for Pr[T x ≥ t] was derived in Sect. 8.4, and was coded-up in R as the TPXG(.) function. So, applying the exact same logic that led to Eq. (9.2), namely to compute the probability of living to the time at which the portfolio is exhausted L ξ , we arrive at the equivalent expression:

which is a bit messier than Eq. (9.2), but does share similarities in structure.

Recall that the Gompertz survival probability to any time t, can be expressed as:

The extra step in Eq. (9.3) is to replace e t with the exponential of portfolio longevity e L ξ . Replacing the L ξ with 1 v ln ξ ξ −v per Eq. (9.1) leads to (ξ/(ξ − v)) 1/v . Bear in mind that the same warnings about ξ > v apply with this version of ϕ as well. If t = L ξ := ∞, then e t → ∞, then (1 − e t ) → −∞, and the lifetime ruin probability equation (9. 3) goes to zero. This is exactly what you would expect (and want) from a risk metric when the portfolio lasts forever. Now, on a technical note, the probability purest might wonder why the lifetime ruin probability ϕ also includes the (yes, measure zero) event in which T x = L ξ . After all, if you go broke the same instant in which you die, isn't that good enough? The answer here is (1) it doesn't make any numerical difference, and more importantly (2) a successful retirement plan is usually defined as dying with positive wealth. Either way this is semantics. For the Gompertz case I'll create a user-defined function for Eq. (9.3).

Note that I forced LRPG(.) = 0 when ξ ≤ v, because the portfolio is never depleted. In fact, this function is relatively simple, so I could have used R to evaluate ϕ for a variety of parameter values (by plugging in the portfolio longevity function PL(.)

as the t value in the Gompertz survival probabilityTPXG(.)) and stopped there. But, the lifetime ruin probability "concept," as well as its converse the retirement success rate, is quite central to many of the recipes in this book, which is why it garners a stand-alone formula. Notice how the early retiree (at age x = 55) is subjecting himself (or herself) to a lifetime ruin probability of ϕ = 32%, whereas the late retiree (at age x = 70) is faced with a ϕ = 1% probability, which is a 99% retirement success rate. Ponder this question again. Is the 4% rule safe or is it risky? Well, first of all it depends on the age at which you begin those withdrawals. Start at the age of x = 55 and it looks extremely risky (to me), but delay until age x = 70 and it seems "safer," with the proviso the safety is in the eye (and tail) of the beholder. More on that in a later section. Now, you might have noticed that in the above script I forced R to only report two significant decimals, which might appear like excessive rounding to some readers. I did this quite deliberately to remind everyone about all of the assumptions that have gone into this simple formula (a.k.a. how I made the sausage). Not only did I assume constant withdrawals, which I critiqued in Chap. 5 (and will do so again), I also made assumptions about Gompertz (m, b) parameters, which one never really knows with certainty. And, while most academic studies (and commercial software) that measure retirement success rates do disclose these assumptions, it's hard to over-emphasize the importance of these implicit parameters in the final results. In fact, to further illustrate how mortality ( The lifetime ruin probability, increases from a comforting ϕ = 4% to a rather scary ϕ = 25%, which is a one-in-four chance of living longer than your money. And remember that all I have changed in the above script is the modal value m, from 88 to 96 years. I have not changed the initial withdrawal rate ξ , nor have I changed the assumed (real) interest rate (a.k.a. investment return), v. You can iterate and change those numbers yourself, just to get a sense of how they impact the results. Finally, for the sake of completeness, let's see how the Gompertz dispersion parameter b affects the lifetime ruin probability ϕ, without changing anything else. Increasing the dispersion coefficient b, all else being equal, increases the value of ϕ. This should make intuitive sense, since we are increasing the uncertainty (think standard deviation) associated with your remaining lifetime, so it's no surprise that greater uncertainty makes things worse. Of course, to prove this rigorously (for all values of ξ, v, x, m), one would have to take partial derivative of ϕ, with respect to b, in Eq. (9.3), and show that all "signs" are positive. I leave this (calculus) problem as a technical exercise for ambitious students.

To extract some final insights from the topic of lifetime ruin probabilities under deterministic returns, I offer a summary plot in Fig. 9 .1. For the sake of space and brevity, I haven't displayed the usual R commands that generated this particular plot, but the syntax should be (more than) familiar by now. In fact, going forward I will refrain from displaying any more plot(.) scripts unless there is something new and important within the R code itself. Now, as far as Fig. 9 .1 is concerned, I have assumed that (m = 88, b = 10), as well as v = 2%, which is deterministic interest, and have varied the initial withdrawal rate ξ , as well as the retirement (withdrawals begin at) age x. Once again there should be no surprises, or any counter-intuitive results. At higher values of ξ , the lifetime ruin probability is higher. Likewise, if you start the withdrawals at more advanced (later) ages, the LRP values drop, and quite dramatically. The main takeaway is the sheer difference a few extra years of work, or withdrawal percentage points, can make on the lifetime ruin probabilitieseven without varying the asset allocation or changing the investment returns. Now, I will certainly get to (what I called) stochasticṽ ∼ N(ν, σ ), in a few pages, but the objective here is develop full intuition for the deterministic case, first.

I can replicate (closely) the above analytic results via simulation by actually randomizing human longevity and then counting the number of people (i.e. lives) who exceeded portfolio longevity. Of course, there is no need to estimate probabilities via simulation when you have access to a closed-form analytic expression. However, since I plan to eventually migrate to simulations when analytic expressions aren't available, this gives us an opportunity to measure some errors, but when the "truth" is known. So, let's go back to the case when (the interest rate) v = 1% and (the withdrawal rate) ξ = 4%, starting at the age of x = 70. We already know that portfolio longevity is PL(0.01,0.04,1)=28.77 years, and the Gompertz survival probability to age y = 98.76 is TPXG(70,28.76,88,10)=0.06281, which is also the value of LRPG(0.01,0.04,70,88,10)=0.0626. But now, I'll simulate N = 1000 random Gompertz lifetimes, conditional on age x = 70, and count the number of those people who actually lived beyond age 98.76 = (70 + L 0.04 ). The script is simply:

zdat<-GRAN(1000,70,88,10) sum(zdat>=PL(0.01,0.04,1))/length(zdat) [1] 0.051

The estimated value (in this one simulation run) of ϕ = 5.1% is obviously lower than the "true" (analytic) value of ϕ = 6.26%. This is to be expected in a simulation with only N = 1000 lifetimes. Indeed, if you take the time to execute this small script, your numbers might be closer to (or farther from) the 6.26%, which is why it's always preferable to use an analytic expression if you have (or know) one. But, what I would like to investigate now is exactly how many random lifetimes I have to simulate before the estimated LRP values are "close enough" to the true value of ϕ. I'm driven by more than simple curiosity. This is a rather vexing issue in practice, when researchers and practitioners report simulated (a.k.a. Monte Carlo) values for retirement success and failure rates. Obviously the bigger the N , the better, but simulations take time and slow down the algorithms. So, how many runs are enough?

Well, Fig. 9 .2 offers some perspective on this question, in the context of our (simple) ϕ. It displays the results of a repeated experiment which simulated samples of N random lifetimes, from N = 100 all the way to N = 10,000. In each experiment (that is for each tested value of N ), I computed all the estimated ϕ values, and then plotted the sample mean plus/minus two standard deviations. In Fig. 9 .2, I also traced out the horizontal line corresponding to the (true) value of ϕ = 6.26%. And, for the sake of replication, here is the core portion (recipe) of the R script generating Fig. 9 .2.

for (i in 1:100){ for (j in 1:50){ zdat<-GRAN(100 * i,70,88,10) y[j]<-sum(zdat>=PL(0.01,0.04,1))/length(zdat) } y.max<-mean(y)+2 * sd(y); y.min<-mean(y)-2 * sd(y) segments(i,y.min,i,y.max,col="black",lwd=2) }

If you want to reproduce (a similar graph), don't forget to add the usual plot(.) commands at the beginning. For now, look carefully at the two loops within the script. The inner one: for (j in 1:50), generates 50 simulations of N = i × 100 random lifetimes, and the outer loop is: for (i in 1:100).

In each run of the outer loop (that is for each value of N ) the script computes the mean(.) and sd(.) of the 50 samples, and then plots the relevant segment from mean(.)+2 * sd(.) to mean(.)-2 * sd. Notice that when N < 2000, the upper and lower bounds are well above and below the true value of 6.26%, by quite a large margin. But, as N increases the variation dissipates, although never reaches zero for the range of N displayed. There are a number of practical takeaways here, some obvious and some less so. First, simulation-based estimates should always be displayed with ranges (obviously). Second, and perhaps more controversially, if we accept the fact that retirement income simulation probability values are rather meaningless after the 2nd or 3rd digit, it might not be necessary to generate more than a few thousand lifetimes to get a reasonable estimate (plus/minus ranges) for the value of ϕ. Thus, while I concede that Monte Carlo simulations (in statistics, physics, engineering, etc.) should be conducted with values of N in the hundreds of thousands, and possibly millions, it's overkill when it comes to retirement income lifetimes, and will create a false sense of precision. After all, how much do you really trust your inputs?

Up to this point I have made a number of assumptions worth repeating. First, in addition to the constant withdrawal rate ξ , I assumed a known, fixed real investment return v. In practice, retirees own pensions and annuities, which guarantee an income for life-a topic I'll return to in Chaps. 10 and 11. But even in the absence of these unique longevity-hedging products, they also hold diversified portfolios which are subject to unknown rates of return. Indeed, that was the focus of Chaps. 5, and 6.

So, I will now combine the technology (scripts) for simulating human lifetimes with portfolio lifetimes to compute more realistic values of ϕ = Pr[T x ≥ L ξ ]. You will see that the algorithm is quite similar to what I did in Sect. 8 The variable zdat.pl contains N = 5000 simulated values of portfolio longevity, assuming an initial withdrawal rate of c/F = ξ = 4%, and ν = 2.5% (which, recall, is the expected value of the log return) with a standard deviation of σ = 15%. The summary statistics for portfolio longevity displays the usual spread between very low values (11 years) and the (artificial, forced) maximum of 100 years. In parallel, the variable zdat.hl contains N = 5000 simulated values of human longevity, conditional on age x = 65 under a Gompertz law of mortality, with m = 88 years and b = 10 years. A summary of these 5000 numbers offers no surprises. The earliest death occurred in 1 week, and the longest (of the 5000) survivors made it to 45.34 years, which is a few months after age 110. Think of two large urns from which pieces of paper are drawn at retirement. The first urn informs you how long you will live and the second urn informs you how long your money will last. Bottom line: If the former number is greater than the later number, you are ruined. Figure 9 .3 displays (a sample of) those 5000 pairs of values, with a diagonal line representing the (rare) cases when the two numbers (in the two urns) match exactly. The first thing to notice is all the dots clustered at the portfolio longevity value of 100, which remember implies that the money lasted for at least 100 years of retirement; more than enough. The dots under the (red) diagonal line represent "urn pairs" for which the human longevity (x-axis) was greater than portfolio longevity (y-axis). Falling under the diagonal is obviously not a desired outcome, but there aren't that many pairs in that region. The majority of dots are above the diagonal, which means that the money lasted longer than your retirement, and in many cases far longer. Zoom in and look at the number of dots (a.k.a. pairs) that are clustered just above the diagonal (phew, that was a close call!) and those clustered under the diagonal (ugh, if only the money had lasted another year or two!). Indeed, those dots that are just under the diagonal make a strong case for portfolio insurance (a.k.a. collars, options, and derivatives), topics which were discussed towards the end of Chap. 6. Finally, I count the number of dots within the triangle under the diagonal, or more precisely the overall fraction for which portfolio longevity L ξ was lower than human longevity T x . This is computed as:

# Crude Simulation-based Estimate of Varphi. > sum(zdat.hl>=zdat.pl)/length(zdat.hl) [1] 0.1542

In this simulation run, the value of ϕ = 15.4%, under the above-noted parameter values. Or, more optimistically, the retirement success rate was (1 − ϕ) = 84.6%. Notice that I am quite cautious about declaring whether this (15.4%) is a good or a bad result, nor am I opining on whether a proper and safe ϕ metric should be higher or lower. I'll get to that (thorny, controversial) topic at the very end of this chapter. For now, I'm just "cooking" numbers. To conclude this particular discussion, Fig. 9 .4 combines 9 individual figures into one, using the par(mfrow=c(3,3)) command in R. After you type that, the next 9 figures created will be stacked into 3 × 3 boxes, similar to Fig. 9 .4. It might not be easy to read the labels or axis, but at least you get all the pictures into one figure. After you are done generating the (nice) 3 × 3, make sure to run par(mfrow=c(1,1)), returning to normal settings for your next figure. Figure 9 .4 varies the initial withdrawal rate ξ , from a low value of 1% to a (very high) value of ξ = 9%. Obviously I'm not suggesting (or even hinting) that 9% would ever be suitable, but rather I'm interested in the pattern of the resulting ϕ estimates. Notice how a larger fraction of the dots fall under the diagonal, as ξ increases. In fact, the cluster of dots at the horizontal value of 100 years declines Fig. 9.4 Visualizing LRP: initial withdrawals from ξ = 1% to ξ = 9% (i.e. drops from the ceiling) as the initial withdrawal rate increases. Likewise, the lifetime ruin probability ranges from ϕ ≈ 0%, all the way to ϕ = 63%, when ξ = 9%. Remember, this represents an x = 65 year-old, subject to Gompertz parameters m = 88 years and b = 10 years, withdrawing ξ from a portfolio earning ν = 2.5% (real) with a standard deviation of σ = 15%. If you don't like any of these assumptions (a.k.a. ingredients), you have the algorithm (a.k.a. recipe) to generate (a.k.a. cook) your own results and figures.

Although there aren't any known closed-form expressions for ϕ, when investments (and L ξ ) are stochastic, there are a number of decent approximations that have been developed over the years (some that are quite elegant). One of those approximations, based on the concept of Statistical Moment Matching, is the focus of the next two sections. In fact, ϕ can also be expressed as the solution to a partial differential equation (PDE), which takes us far beyond the scope of this (recipe) book. See the reference cited as [3] . So, for those readers who might not relish the stochastic calculus, and for whom simulation estimates are more than sufficient, I suggest skipping to the final section (9.8) of this chapter, where I take an axe to ϕ. But for those who are brave enough, I move on to some definitions and formalities.

Recall from Chap. 6, that the vector or stochastic process Z t , t ≥ 0 denoted the real (inflation adjusted) investment process driving portfolio longevity. For any given path of Z t , I defined the annualized real investment growth (ARIG) as: (Z T /Z 0 ) 1/T − 1, for any period T . In this section I'll begin by delving (more) deeply into the randomness of Z t . I'll start by defining the investment process in continuous time, via the following stochastic process: + σ B t ) . (9.4) The two parameters in Eq. (9.4) are the familiar ν, σ , which are the (continuously compounded) expected rate of return and the volatility. The new (foreign) B t is the continuous-time analog of the normal distribution which generates the investment returns, which is called a Brownian Motion. And, while I am casually introducing this symbol (and process), one could fill an entire book on the proper definition, construction, and properties of Brownian Motion. I refer you to the reference [8] and [4] for (much) more information about B t , which is famous in its own right. If you have never heard of B t before, you will have to take it on faith that the process Z t can be written and expressed in differential form as: (9.5) where (ν + 1 2 σ 2 ) is known as the drift coefficient, which in the mathematical finance literature is often denoted by the letter μ, and σ is known as the diffusion coefficient of the stochastic process. So, for example, if the expected continuously compounded rate of return is ν = 2% (for example) and the volatility (a.k.a. standard deviation of the continuously compounded return) is σ = 20%, then the drift coefficient in Eq. (9.5) is 4%. This might seem odd at first (Does Z t grow by 2% or by 4%?), but relates to the difference between the geometric (2%) and arithmetic mean (4%), and on a more deeper level (something called) Ito calculus. Nevertheless, bearing in mind that my goal here is recipes and algorithms, versus theory, if the return generating process is governed by Eq. (9.5), then the trajectory of the retirement portfolio will satisfy the following stochastic differential equation (SDE):

where F 0 is the initial portfolio value at retirement and c is the constant consumption rate (a.k.a. ξF 0 ). Compare equation (9.6), which is the stochastic trajectory of wealth, to the deterministic trajectory dF t = (vF t − c)dt, derived in Sect. 5.7, and coded-up as DTRJ(.) in R. Notice the similarities between the two and the fact they match when σ = 0. Now, for the final bit of stochastic calculus, the solution to the Stochastic Differential Equation (SDE) in Eq. (9.6) can be written explicitly as:F

In words: The value of the retirement portfolio at any future time t, is equal to the value of the investment process Z t > 0 multiplied by the difference between the original portfolio value F 0 , minus the consumption c, times the integral of the inverse of the investment return process Z s , until time t. What this means, among other things, is that you can simulate portfolio values ofF τ , at any time τ > 0, by simulating a single path for Z s ; 0 ≤ s ≤ τ , then integrating (i.e. adding up) the reciprocal of Z s from s = 0 to s = τ , and plugging them all into Eq. (9.7). Without any loss of generality, τ can also be the random lifetime based on a Gompertz model, andF T x becomes the value of the portfolio at the time of death. So, ifF T x ≤ 0, this retiree was ruined. But here is the main "trick" in this section. Look closely (again) at Eq. (9.7). The only condition under whichF T x ≤ 0 is if the item in the brackets is negative, because Z t > 0. So, ϕ can be expressed as the probability that the (yes, random) integral expression in Eq. (9.7) is greater than F 0 /c, which is also 1/ξ (and measured in units of initial retirement wealth). Voila! Simulate many (many) values of the random integral, by simulating many paths for Z s , and you have (yet another) estimate of ϕ. Formally, here is yet another expression for the lifetime ruin probability:

Now, strictly speaking you can't really define an integral with an upper bound of integration that is random, which is the time of death. So, the first expression to the right of the equal sign in the above equation isn't a proper integral (to a proper mathematician), which is why I added the second and final expression integrated to infinity, but with an indicator function 1 {T x >s} , to denote the person is still alive. To bring some clarity to the integral in Eq. In the first run (which you will notice, is a bit slower than our prior simulations), the retiree spends ξ = 5%, but only earns an expected continuously compounded real return of 2% with a standard deviation of σ = 20%. In that scenario the lifetime ruin probability is ϕ = 35.25%, which is quite high. In the second simulation the withdrawal rate is reduced to ξ = 4%, resulting in a (better) ϕ = 24.19%. Finally, if the volatility can be reduced from 20% to σ = 15% (somehow), without affecting the expected (continuously compounded) return ν, the lifetime ruin probability can be reduced to 20.79%.

The approach I just described, simulating the integral in Eq. (9.8), suggests a third and final way of approximating the lifetime ruin probability ϕ, when lifetimes are Gompertz distributed and (the logarithm) of investment returns are normally distributed. Without getting lost in the details, there are some good (theoretical) reasons for why the reciprocal of the integral, that is Y = 1/ Z −1 s ds can be approximated by a Gamma distribution. The cumulative distribution function (CDF) can be expressed as:

ds, (9.9) where α > 0 is known as the shape parameter and β > 0 is known as the scale parameter. The mean of the Gamma distribution is αβ, and the variance is αβ 2 So, if I make the (bold, and to this point completely unjustified) assumption that 1/ Z −1 s ds has a Gamma distribution, then the probability Pr[ Z −1 s ds ≥ 1/ξ ], per Eq. (9.8) can be expressed as Pr[1/ Z −1 s ds ≤ ξ ], which is the expression in Eq. (9.9). To get the parameters (α, β), I will match the moments of Z −1 s ds to the Reciprocal Gamma distribution, which is explained in (much) more detail in a technical research paper cited as [5] as well as [3] . Once again, my interest here isn't proofs and theorems, but recipes. To implement this in R, there are a number of functions I have to create first, before I can define a new VARPHI.MM. First, I need a robust function G(a,c) for the incomplete Gamma function. 

To confirm this is working properly, here are two values: a(0.01,65,88,10) which should yield a value of 17.888 as well as the value for a higher discount rate, a(0.03,65,88,10) leading to 14.408, both of which should confirm your code is working properly. I will actually have (much) more to say about this little a(.) function in the next Chap. 10, but for now it's introduced as a stepping stone for the moment-matching approximation. Finally, the user-defined function that computes that approximation can be constructed as follows:

VARPHI.MM<-function(x,m,b,xi,nu,sigma){ mu<-nu+(0.5) * sigma^2 M1<-a(mu-sigma^2,x,m,b) M2<-(a(mu-sigma^2,x,m,b) -a(2 * mu-3 * sigma^2,x,m,b))/(mu/2-sigma^2) alpha<-(2 * M2-M1^2)/(M2-M1^2) beta<-(M2-M1^2)/(M2 * M1) #Shape is Alpha, Scale is Beta pgamma(xi,shape=alpha,scale=beta,lower.tail=TRUE) }

The source, proof, and overall justification for these expressions are in the article cited as [5] , and in particular equations (18) The moment-matching estimates for ϕ are slightly lower (by 2-3% points) compared to the simulation-based results. Note that VARPHI.SM>VARPHI.MM isn't always the case, and is actually reversed at higher withdrawal rates ξ . For example, when ξ = 10%, ν = 5%, and σ = 25%, the LRP under the MM technique is higher by a few percentage points compared to SM. I report results and conclude the section. 

Starting with the Swedish mathematicians Filip Lundberg (b. 1876, d. 1965) and Harald Cramer (b. 1893, d. 1985), insurance actuaries have used ruin probabilities to measure liabilities and insurance company's risk for over a century. In the insurance context it represents the probability a stochastic process subjected to random additions (i.e. premiums) and subtractions (i.e. claims) hits zero (bankruptcy) before some terminal date. Bottom line, the credit for using ϕ in a business context (as opposed to gambling problems) goes to actuaries. But, in the last 2-3 decades these ruin probabilities (and their converse, success rates) have transitioned from the insurance world to quite popular metrics in personal finance, wealth management, and retirement income planning. Indeed, in addition to my own work cited as reference [6] , the first articles in this area were [1] and [2] as well as [9] in the portfolio management literature. Although many didn't use the phrase ruin, nor did they account for mortality risk, the motivation was similar. So, while this entire chapter (and most of this book) has focused on the benefits of thinking about retirement income plans thru the lens of probabilities, I do have some concerns about how ϕ might be abused if taken too far. I have expressed these concerns elsewhere, for example, the article cited as [7] , but in this concluding section I will only highlight caveats about this (fickle, temperamental) metric.

To start, the underlying fix-ξ -until-you-die assumption, namely that people will adhere to a deterministic spending schedule and then wake up one morning, check their account online, and find the money has reached zero, i.e. the portfolio has failed, is silly and naive. Nobody waits until that point to take restorative action to repair their finances. I mentioned this earlier, in Chap. 5, when I introduced portfolio longevity, and emphasize it here again. It's a big assumption, but most commercial software packages and algorithms in the retirement income business suffer from the same weakness. They assume static and pre-determined consumption behavior that is at odds with economic theory. I'll return to the economics in later chapters.

But, here are some additional things to keep in mind, especially now that you understand all the assumptions that go into producing (simulating) ϕ. To be blunt, nobody has any idea what lifetime ruin probability is acceptable in the context of retirement income planning. Is ϕ = 10% too high? Is a healthy ϕ = 1%, or is the one-in-four ϕ = 25% acceptable? I have actually heard practitioners argue that a 51% chance of success (which is a 49% LRP) is good enough for a retirement plan. But, then think about it this way. Would you get on an airplane with a 0.1% probability of crashing? Why should a retirement income plan be any different?

Another issue arises when two different retirement income "strategies" are compared side by side. Assume that Plan A is estimated to have a ϕ = 5%, but for Plan B the number is ϕ = 4%. Ergo, the reader concludes that Plan B is preferred to Plan A. But is a 1% reduction in a lifetime ruin probability really that much better? Think back to the scattered dots in Fig. 9 .3. The 1% reduction in ϕ only implies that a few dots have moved from under diagonal to over the diagonal. Will this make a real difference?

In fact, even if the industry (or the financial engineers) could agree on acceptable thresholds and bounds for ϕ, there are many different statistical distributions or financial plans that can generate the same failure rate. Should a retiree be indifferent between all plans of equal success or failure rate? To echo George Orwell, all 10% shortfall probabilities are equal, but some are more equal than others. Or, in less Orwellian but slightly more technical terms, different statistical distributions can share the same tail probability but have distinct risk and return profiles. Figure 9 .5 illustrates three such curves. It models the distribution of financial legacy in discounted or present value terms, representing (very abstractly) the amount of money left at death (think:F T x ) if you adhere to a fixed ξ forever. All three curves are Gaussian and very well behaved. For Case A the legacy has an average of $100,000 and a standard deviation of $78,000. For Case B, the legacy curve has an average of $200,000 and a standard deviation of $156,000 and for Case C, the average is $300,000 with a standard deviation of $234,000. If you want something more concrete to imagine, think of a very conservative retirement portfolio (Case A); a more balanced asset allocation mix (Case B) and a very aggressive asset allocation (Case C). All three portfolios are subjected to the same withdrawal rate ξ . The bulk Fig. 9 .5 Know your tail magnitudes of the probability distribution is to the right of zero, which means that you can expect to leave a positive financial legacy. Notice that all three curves have the same area to the left of zero (i.e. the dashed line). The lifetime ruin probability is the same 10% under all three, but the severity of the potential rescue operation is quite different. Would retirees be indifferent between the plans? Indeed, it is not uncommon to read simulation-based studies in which counter-intuitive Strategy A is shown to have a higher success rate than more intuitively reasonable Strategy B, because the authors assume a 9% lifetime ruin probability (the area to the left of zero) is better than an 11% LRP, when in fact the entire distribution of outcomes is quite different. In particular, riskier stocks have a nasty habit of reducing lifetime ruin probabilities but wreaking havoc on higher statistical moments. Oddly enough, lower ϕ values are not necessarily better. In plain English, the lifetime ruin probability doesn't seem to care about what happens to all the dots above the diagonal in Fig. 9 .4, but most humans do.

In sum, whether your Monte Carlo is using log-normal returns, bootstrapping historical data, or using a forward-looking equilibrium model for investment returns, you are implicitly making assumptions about the world in the very distant future. Effectively, your black box is forecasting how interest rates, stock prices, inflation as well as mortality will evolve over the next 50 years and how they will co-vary with each other. If all these assumptions were questionable before COVID-19, they are quite dubious afterwards. Please use ruin with caution!

Determining withdrawal rates using historical data

Sustainable withdrawal rates from your retirement portfolio

Ruined moments in your life: How good are the approximations

Numerical solutions of SDE through computer experiments

Self-annuitization and ruin in retirement

A sustainable spending rate without simulation

It's time to retire ruin (probabilities)

Stochastic differential equations: An introduction with applications

Sustainable investment withdrawals