Gompertz, annuities, and special functions

This is part 2 of a 2 part series on Gompertz’s law. Part 1 is here.

The first part of this series discussed the Gompertz distribution and gave a formula for “How likely am I to live until at least age \(t_0 + t\) conditional on having lived until age \(t_0\)?”. In this post, we will use this to compute the present value of an annuities. This is a riff on chapter 3 of Moshe Milevsky’s book The 7 Most Important Equations for Your Retirement with a more mathy twist.

Present value

A standard idea in finance is that getting \(\$1000\) today is worth more than getting \(\$1000\) in a year from now. This should be true even if these are inflation-adjusted dollars: If we assume some bank will give us the expected rate of inflation as their interest rate when we deposit the \(\$1000\), we could just do that and create the “\(\$1000\) a year from now” situation. Ignoring the dilemma of choice, the first situation gives us more options and should therefore be worth more to us.

How much more should it be worth? Or, equivalently, what’s the value for \(y\) at which point we become ambivalent about \(\$x\) now and \(\$y\) in one year?

The answer will be different for different people (you may really need that money now), but generally this is modeled with an interest/discount rate as well: \(\$x\) today will be worth \(\$x\cdot(1+r)\) a year from now (and therefore \(\$x\cdot (1+r)^2\) two years from now) for some discount rate \(r\). Why the temporal structure is exponential as opposed to some other increasing function is a topic all for itself (look up hyperbolic discounting), but it’s a common choice in various fields including finance, economics, psychology, neuroscience, and reinforcement learning.

Conversely, \(\$y\) in one year should be worth \(\frac{\$y}{1+r}\) today.

Getting money forever

A nice property of this discounting rule is that the option of getting a neverending stream of money (say, \(a = \$1000\) every month) is worth a finite amount today by the magic of the geometric series:

for a monthly discount rate \(r\ne 0\). Just to be extra clear: The \(n\)th term of this sum is the present value of the payment we receive in the \(n\)th period.

Such an arragement is called an annuity. We talked about annuity loans before on this blog and just as in that case, the present value of a neverending stream of \(a\) per month is \(a\) times an annuity factor, in this case \(\af(r)=1 + \frac{1}{r}\).

The same is true if the annuity runs out after a fixed number of periods. However, as a component for retirement planning, annuities that pay for the rest of ones natural life (or the life of ones spouse) tend to be a better option, as they help covering longevity risk. In fact, annuities are (in theory, although maybe not psychologically in practice) a great component of a retirement plan as they allow pensioners to take more risk with the rest of their funds, e.g., stay invested in the stock market where expected returns are higher (but so is the variance of returns).

Getting money until you die

This raises a question: What’s the present value of an annuity that pays for the rest of ones life? It ought to be less than \(1+\frac{1}{r}\) times the annual amount (unless one expects to live forever), but how much less exactly?

Our answer, of course, is the survival function of the Gompertz distribution, as per the last blog post. Of course this depends on the current age of the pensioner: The survival function of the Gompertz distribution, conditioned on having lived until year \(t_0\), is

where \(b \approx 1/9.5\) is the inverse dispersion coefficient and \(t_m \approx 87.25\) is the modal value of human life.

The value of a pension that starts now (with a pensioner at age \(t_0\)) and pays \(\$1\) per year for the rest of the life is therefore

The \(n\)th term in this series is the probabilty of being alive at year \(t_0 + n\) given having been alive at year \(t_0\) times the discount rate that turns it into a present value.

This is equation #3 in Milevsky’s book The 7 Most Important Equations for Your Retirement.

Some insurance contracts allow the pensioner to decide between a lump sum payment or an annuity; this formula can help to compare these two options. Some contracts provide payments until death, but with a minimum of (say) 10 years. Replacing \(p_{t_0}(0), \ldots, p_{t_0}(9)\) with \(1\) in \(\eqref{discraf}\) would model this situation.

Milevsky suggests using a spreadsheet program to sum this series, along the lines of

In this example a pensioner of age 65 buys an annuity that gives them \(\$1\) every 5 years (starting at age 70) with a discount rate of \(r=7.25\%\). The present value of that annuity is \(\$1.51\). If the payment every 5 years is \(\$50.000\) instead, the present value of the annuity at age 65 is \(1.507\cdot \$50.000 = \$ 75369.90\).

The example uses 5 year intervals in order to not get unduly long. But typically, annuities will pay monthly. Computing the sum in a spreadsheet then becomes somewhat annoying: Here’s an monthly example in Google sheets. (Notice that the difference to the 5 year interval version actually matters: In the latter case, one gets the money for five years in advance; in the former case one gets it only for the current month. A lot can happen in 5 years.)

I wasn’t quite satisfied with this.

The hunt for a closed-form solution

Equation \(\eqref{discraf}\) is nice and all, but it would be much better to have a closed-form expression for it. Since \(p_{t_0}\) is a doubly-exponential function, the solution isn’t immediately obvious. In fact, I don’t know a closed-form solution for \(\eqref{discraf}\)¹.

But let’s imagine a world with continuous banking, where instead of once a month we get a small portion of our annuity every moment. The continuous annuity factor would then be

where \(\eta = e^{b(t_0 - t_m)}\) as in the last blog post.²

This, too, doesn’t look super easy. In fact, it’s not solvable with “elemantary” functions. But somewhere within the large zoo of special functions there is the right one for us.

In this case, Wikipedia already tells us that the moment-generating function (aka Laplace transform) of the Gompertz distribution is

with the generalized exponential integral

The Gompertz distribution has the nice property that “left-truncated” versions of itself are still Gompertz distributed (see the last blog post for details). This is a consequence of the fact that its survival function \(S(x) = \exp(-\eta(e^{bx}-1))\) shows up in its pdf \(f(x) = b\eta S(x)e^{bx}\).

Hence,

So the moment-generating function at \(t = r + b\) gives us a closed-form for \(\eqref{contaf}\):

The annuity factor \(\af_c(r, t_0)\) for different interest rates \(r\) and starting ages \(t_0\). Retiring early is expensive. Gnuplot source

Trying to use this

Formula \(\eqref{contaf-gef}\) is in fact a closed-form solution to our problem. So let’s try to use it in a computation with Python. SciPy has the generalized exponential integral function, so this should be easy:

import numpy as np
from scipy import special

b = 1 / 9.5
t_m = 87.25


def af(r, t_0):
    eta = np.exp(b * (t_0 - t_m))
    return np.exp(eta) / b * special.expn(1 + r / b, eta)


print(af(0.025, t_0=65))

and the output is:

gompertz_discount.py:10: RuntimeWarning: floating point number truncated to an integer
  return np.exp(eta) / b * special.expn(1 + r / b, eta)
19.439804660538815

Ah, dang! Even though the generalized exponential function is part of SciPy, the implementation there only supports integer orders.

We’ll have to visit the special functions zoo for a bit longer.

More zoo animals: The incomplete gamma function

The classic reference for special functions is Abramowitz and Stegun. Wikipedia has a quote about it from the (American) National Institute of Standards and Technology:

More than 1,000 pages long, the Handbook of Mathematical Functions was first published in 1964 and reprinted many times […] [W]hen New Scientist magazine recently asked some of the world’s leading scientists what single book they would want if stranded on a desert island, one distinguished British physicist said he would take the Handbook. […] During the mid-1990s, the book was cited every 1.5 hours of each working day.

In the internet age, the Digital Library of Mathematical Functions (DLMF) hosts an updated version of this classic work and it is very useful for situations like the one we are in now.

In fact, looking at equation 8.19.1 in DLMF, we see that the generalized exponential integral is nothing but the (upper) incomplete gamma function:

and thus

Alas, this won’t do either. SciPy’s implementation of the incomplete gamma function doesn’t allow for a negative first argument.

Why might this be the case? Looking at the definition of the incomplete gamma functions in DLMF (8.2.2 there), we read

Note how this is an “incomplete” version of the standard gamma function \(\Gamma(a) = \Gamma(a, 0)\). Also note how this integral won’t work for negative \(a\): There’s a non-integrable singularity when \(z=0\).

So did our previous formulae not make any sense? They did. We just have to understand these functions a little bit better.

Analytic continuation

Complex analysis, the theory of differentiable complex functions, is one of the most beautiful areas of mathematics. In fact, I believe I once saw a video of Donald Knuth saying it’s such a great topic that he asked his daughter to attend a complex analysis course – although she didn’t otherwise study at any university.³

One of the results of complex analysis is that if two nice (aka complex-differentiable, aka holomorphic, aka analytic) functions coincide on a set that’s not totally discrete, they are the same! This means any analytic function like \(a\mapsto\Gamma(a, z)\) has at most one analytic continuation. This is an important principle which also applies to the famous Riemann zeta function \(\zeta\). The zeros of its analytic continuation to the left complex half-plane are what the most important open problem in mathematics is about (you’ll get $1M if you solve this one, which should help with retirement planning). In fact, \(\zeta\) and \(\Gamma\) are closely linked.

In the case of \(\zeta(s) := \sum_{n=1}^\infty n^{-s}\) its unique analytic continuation to the left half-plane yields for instance \(\zeta(-1) = -\frac{1}{12}\) which gives some sense to Ramanujan’s mysterious equation \(1 + 2 + 3 + \cdots = -\frac{1}{12}.\)

How is such an analytic continuation found? In the case of both the zeta and the gamma functions, it’s via functional equations. For instance, \(\Gamma(n) = (n-1)!\) for integer \(n\), and more generally \(\Gamma(z+1) = z\Gamma(z)\) for complex \(z\). If we start with some \(z\) in the left half-plane (i.e., \(\Im z \le 0\)), repeatedly applying this formula eventually yields only terms of \(\Gamma\) at “known” arguments, which establishes the value of the analytic continuation of \(\Gamma\) to that point \(z\).

For the incomplete Gamma function, a similar recurrence relation exists, namely (8.8.2 in DLMF)

for, say, \(a\) and \(z\) in the right half-plane. This equation can be used to extend the incomplete gamma function to negative values of \(a\). It also gives rise to a corresponding recurrence relation for the generalized exponential integral (8.19.12 in DLMF)

Applying the functional equation

Since the implementations of \(\mathrm{E}_{p}(z)\) and \(\Gamma(a,z)\) that we want to use don’t implement the continued versions of these functions, we can use the functional equations that define the continuations ourselves. The recurrence relation for \(\mathrm{E}_{p+1}\) yields

For moderate values of \(r\) (and \(b\)), this is enough. For very large \(r\) we’d have to apply the recurrence relation again.

This allows us to finally run this in Python:

import numpy as np
from scipy import special

b = 1 / 9.5
t_m = 87.25


def af(r, t_0):
    eta = np.exp(b * (t_0 - t_m))
    return (
        1
        - np.exp(eta)
        * eta ** (r / b)
        * special.gamma(1 - r / b)
        * special.gammaincc(1 - r / b, eta)  # A normalized version of \Gamma(a, z).
    ) / r


print(af(0.025, t_0=65))

Which prints 14.79901377449508.

What about spreadsheets?

The above works for SciPy, but Excel (or Google Sheets) is a different matter as it doesn’t directly have support for the incomplete gamma function.⁴ However, it does support the Gamma distribution which has a normalized version of the lower incomplete gamma function as its cdf. As the helpful authors of Wikipedia note, \(\Gamma(s, x)\) can be computed as GAMMA(s)*(1-GAMMA.DIST(x,s,1,TRUE)).

This finally gives us closed-form we can use in spreadsheets. Here’s an example in Google Sheets. Note that the relative difference of \(\af(r,t_0)\) and \(\af_c(r,t_0)\) for \(r = 0.25\%\) and \(t_0 = 65\) is only \(0.64\%\).

What does it all mean?

If you (or your parents) bought an annuity in your name and the issuer offers the option for a lump sum payment instead, you can multiply the yearly (or monthly) payments with the annuity factor as determined above, treat that as the annuity’s value, and compare it to the lump sum. (The computations here mostly assume the annuity will start right away, extending this to an annuity that starts in a number of years is left as an exercise for the reader.)

That raises the question of which interest rate to pick? And why that should be treated as constant over time?

Also, the annuity provider will have done some math as well and will be aware of adverse selection: People who self-assess as being in good health might be right about this and choose the annuity more often. Milevsky quotes from Jane Austen’s Sense and Sensibility,

People always live forever when there is an annuity to be paid them.

Milevsky also ends his discussion of his equation #3 with a warning: This is a model. If market prices are different, they are likely to be right and the model wrong.

That said, if you got an annuity a while back, for instance pre-2008 when the effective Fed funds rate was 5%+, you likely got a good deal by today’s standards. The future might be different. But at least some people argue that interest rates will stay permanentely low, as they did in Japan for the last decades (demographic factors might play a role here). Others disagree, and of course people some always predict stronger inflation.

That’s it for today, thanks for reading. Reach out if you have comments.