Gompertz-Makeham Survival Model - Infinity is Really Big

Studying for Exam LTAM, Part 1.7

You susceptibility to death generally grows higher as you age. This is an unfortunate fact of life. However, life insurance can help protect your family if you should die. It is a good idea to purchase life insurance, especially if you have dependents.

So far we have explored modeling human lifetimes with the following types of continuous distributions: uniform (De Moivre’s Law), exponential (constant force of mortality), and triangular distributions. The first two of these are simpler and commonly used over short periods of time. The last is more complicated but has the potential to be more accurate over longer periods of time.

However, despite this potential, the triangular distribution is not as accurate over long periods as other models that are more commonly used. In this article we will consider a better and more commonly used model defined by both Gompertz and Makeham. It is called the Gompertz-Makeham survival distribution. This model is both logical and reasonably accurate.

Force of Mortality for Gompertz-Makeham

The Gompertz-Makeham survival distribution starts with the assumption that “instantaneous risk of death” has two components: 1) a constant term that everyone is susceptible to, and 2) a term that increases exponentially over time.

The constant term can be thought of as arising from a kind of “general average” risk of death that affects everyone basically evenly, perhaps from, for example, accidents. Obviously such a term is a simplification. However, such simplifications are often necessary in mathematical modeling.

The exponential term can be thought of as arising from the fact that aging affects the body in detrimental ways. The most obvious way this occurs is an increased risk of terminal disease as a person ages.

Therefore, in the Gompert-Makeham model, we assume that the force of mortality is $\mu_{t}=A+Bc^{t}$ for some constants $A,B>0$ and $c>1$ . The values of these parameters would have to be determined from real-life data. Plausible values of these parameters for human lifetimes might be $A=0.0001, B=0.0003$ , and $c=1.07$ .

Obviously, we can also write this model as $\mu_{t}=A+Be^{kt}$ , where $k=\ln(c)>0$ . However, it is sometimes better to keep it in the form $\mu_{t}=A+Bc^{t}$ to emphasize the percent that the term $Bc^{t}$ grows every year. That percent, as a decimal, would be the value of $c-1$ .

Historically, Gompertz was the first to study this model, but only in the case where $A=0$ (from the year 1825, see “On the Nature of the Function Expressive of the Law of Human Mortality”). Makeham modified it later to include the case $A>0$ (from the year 1860, see “On the Law of Mortality, and the Construction of Annuity Tables”).

Survival Function for Gompertz-Makeham

Given survival to age $x>0$ and corresponding continuous remaining lifetime random variable $T_{x}$ , the survival function is the complement of the cumulative distribution function $\,_{t}p_{x}=P[T_{x}>t]=e^{-\int_{0}^{t}\mu_{x+s}\, ds}$ . The last equality follows from the definition of the force of mortality as $\mu_{x+t}=-\frac{\frac{d}{dt}\left(\,_{t}p_{x}\right)}{\,_{t}p_{x}}=-\frac{d}{dt}\left(\ln\left(\,_{t}p_{x}\right)\right)$ .

When we assume that $\mu_{t}=A+Bc^{t}$ , the integral to evaluate is $\displaystyle\int_{0}^{t}\left(A+Bc^{x}c^{s}\right)\, ds=\left(As+\frac{Bc^{x}}{\ln(c)}c^{s}\right)\biggr\rvert_{s=0}^{s=t}=At+\frac{B}{\ln(c)}c^{x+t}-\frac{Bc^{x}}{\ln(c)}$ . We can then conclude that $\,_{t}p_{x}=e^{-At-\frac{Bc^{x}}{\ln(c)}(c^{t}-1)}=\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)$ .

Fixing $B=.0003, c=1.07$ , and $x=0$ while letting $A$ increase from .0001 to .1 results in the following animation of the graph of $\,_{t}p_{x}$ for $0\leq t\leq 120$ . Certainly the resulting survival function is more accurate for small values of $A$ . Also note that, when $A$ is large and $t$ is small, the value of $A$ is so dominant in $\mu_{t}=A+Bc^{t}$ that the survival function $\,_{t}p_{0}$ appears exponential (constant force).

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $A$ increases from .0001 to .1, for $x=0, B=.0003$ , and $c=1.07$ . Note that, when $A$ is large and $t$ is small, the value of $A$ is so dominant in $\mu_{t}=A+Bc^{t}$ that the survival function $\,_{t}p_{0}$ appears exponential (constant force).

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $A$ increases from .0001 to .1, for $x=0, B=.0003$ , and $c=1.07$ . Note that, when $A$ is large and $t$ is small, the value of $A$ is so dominant in $\mu_{t}=A+Bc^{t}$ that the survival function $\,_{t}p_{0}$ appears exponential (constant force).

Fixing $A=.0001, c=1.07$ , and $x=0$ while letting $B$ increase from .0003 to .1 results in the following animation of the graph of $\,_{t}p_{x}$ for $0\leq t\leq 120$ .

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $B$ increases from .0003 to .1, for $x=0, A=.0001$ , and $c=1.07$ . The visual effect of increasing $B$ here seems more like a horizontal leftward shift than the previous animation, which might be described as more ‘wave-like’. Also note that when $B$ is close to .1, the decay in this graph is greater than in the first graph. That’s because of the factor of $c^{t}.$

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $B$ increases from .0003 to .1, for $x=0, A=.0001$ , and $c=1.07$ . The visual effect of increasing $B$ here seems more like a horizontal leftward shift than the previous animation, which might be described as more ‘wave-like’. Also note that when $B$ is close to .1, the decay in this graph is greater than in the first graph. That’s because of the factor of $c^{t}.$

The effect of the increase in $B$ has a somewhat similar overall qualitative effect of the increases in $A$ . If I were to describe any subtle differences in these effects, I would say that the effect of $B$ increasing “seems more like” a leftward horizontal shift of the graph than the effect of $A$ increasing, which seems more like a slanted “oscillating wave”. Also note that when $B$ is close to .1, the decay in this graph is greater than in the first graph. That’s because of the factor of $c^{t}$ .

Increasing the value of $c$ produces an effect that seems to be even more similar to a leftward horizontal shift compared to increasing $B$ . The animation below shows this as $c$ increases from 1.07 to 1.20.

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $c$ increases from 1.07 to 1.07, for $x=0, A=.0001$ , and $B=.0003$ . The visual effect of increasing $c$ here seems even more like a horizontal leftward shift than the previous animation.

The graph of $\,_{t}p_{x}$ over $0\leq t\leq 120$ as $c$ increases from 1.07 to 1.07, for $x=0, A=.0001$ , and $B=.0003$ . The visual effect of increasing $c$ here seems even more like a horizontal leftward shift than the previous animation.

Probability Density Function for Gompertz-Makeham

Since $\,_{t}p_{x}$ is the complement of the cumulative distribution function (CDF) of $T_{x}$ , we can say that the probability density function (PDF) of $T_{x}$ is $f_{x}(t)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)$ .

But did you know that the PDF of $T_{x}$ also equals $\,_{t}p_{x}\mu_{x+t}$ ? In fact, this is a very common way that actuaries write the PDF.

Let us see if this really is true for our example. From above, we know that $\mu_{x+t}=A+Bc^{x+t}=A+Bc^{x}c^{t}$ and $\,_{t}p_{x}=\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)$ .

Now compute:

$-\frac{d}{dx}\left(\,_{t}p_{x}\right)=-\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)\cdot (-A-Bc^{x}c^{t})$

$=(A+Bc^{x+t})\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)$ .

And clearly, the other one equals the same thing:

$\,_{t}p_{x}\mu_{x+t}=(A+Bc^{x+t})\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)$ .

But why are these equal in general? It just goes back to the definition of $\mu_{x+t}$ . We know that $\mu_{x+t}=-\frac{\frac{d}{dt}\left(\,_{t}p_{x}\right)}{\,_{t}p_{x}}$ . Hence, $\,_{t}p_{x}\mu_{x+t}=\,_{t}p_{x}\cdot \left(- \frac{\frac{d}{dt}\left(\,_{t}p_{x}\right)}{\,_{t}p_{x}}\right)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)=f_{x}(t)$ .

The corresponding animations of the graphs of the PDFs are even more interesting than the ones above. Take the time to think about these in relation to the animations of the graphs of the survival functions above.

The graph of $f_{x}(t)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)=\,_{t}p_{x}\mu_{x+t}$ over $0\leq t\leq 120$ as $A$ increases from .0001 to .1, for $x=0, B=.0003$ , and $c=1.07$ . Note that, when $A$ is large and $t$ is small, the value of $A$ is so dominant in $\mu_{t}=A+Bc^{t}$ that the PDF $f_{0}(t)$ appears exponential (constant force).

The graph of $f_{x}(t)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)=\,_{t}p_{x}\mu_{x+t}$ over $0\leq t\leq 120$ as $A$ increases from .0001 to .1, for $x=0, B=.0003$ , and $c=1.07$ . Note that, when $A$ is large and $t$ is small, the value of $A$ is so dominant in $\mu_{t}=A+Bc^{t}$ that the PDF $f_{0}(t)$ appears exponential (constant force).

The graph of $f_{x}(t)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)=\,_{t}p_{x}\mu_{x+t}$ over $0\leq t\leq 120$ as $B$ increases from .0003 to .1, for $x=0, A=.0001$ , and $c=1.07$ . This looks similar to an exponential decay when $B$ is large, except for the initial concave down part of the graph near $t=0$ .

The graph of $f_{x}(t)=-\frac{d}{dt}\left(\,_{t}p_{x}\right)=\,_{t}p_{x}\mu_{x+t}$ over $0\leq t\leq 120$ as $c$ increases from 1.07 to 1.20, for $x=0, A=.0001$ , and $B=.0003$ . This has a similar “shape” no matter how big $c$ gets, though it gets more “compact”.

Complete Expectation of Life for Gompertz-Makeham

The complete expectation of life is the mean (expected value) of $T_{x}$ . It is denoted by $\stackrel{\circ}e_{x}$ and equals $\stackrel{\circ}e_{x}=\displaystyle\int_{0}^{\infty}\,_{t}p_{x}\, dt$ .

For the Gompertz-Makeham model, this becomes $\stackrel{\circ}e_{x}=\displaystyle\int_{0}^{\infty}\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)\, dt$ . Suffice it to say, this is basically an impossible integral to do (in fact, Mathematica “gives up”). It seems that the function $\exp\left(\frac{Bc^{x}}{\ln(c)}(1-c^{t})-At\right)$ has no elementary antiderivative in terms of $t$ .

However, numerical integration can be used to approximate this integral for various values of $x$ . When this is done, and the resulting values are plotted, we get the following graph when $A=.0001, B=.0003$ , and $c=1.07$ . Note that the expected remaining lifetime at birth is about 72 and decreases down toward 0 as $x$ increases toward 120. The graph is always concave up and always has a slope between $-1$ and $0$ . You might want to see if you can figure out why this is true.

The graph of the complete expectation of life for the Gompertz-Makeham model when $A=.0001, B=.0003$ , and $c=1.07$ . Note that the expected remaining lifetime at birth is about 72 and decreases down toward 0 as $x$ increases toward 120. The graph is always concave up and always has a slope between $-1$ and $0.$ You might want to see if you can figure out why this is true.

The graph of the complete expectation of life for the Gompertz-Makeham model when $A=.0001, B=.0003$ , and $c=1.07$ . Note that the expected remaining lifetime at birth is about 72 and decreases down toward 0 as $x$ increases toward 120. The graph is always concave up and always has a slope between $-1$ and $0.$ You might want to see if you can figure out why this is true.

Finally, below we graph the complete expectation of life $\stackrel{\circ}e_{x}$ along with the two functions $\stackrel{\circ}e_{x}\pm 2\sigma_{x}$ , where $\sigma_{x}$ is the standard deviation of $T_{x}$ . Much numerical integration was necessary here as well.

The graphs of $\stackrel{\circ}e_{x}$ and $\stackrel{\circ}e_{x}\pm 2\sigma_{x}$ for the Gompertz-Makeham model when $A=.0001, B=.0003$ , and $c=1.07$ . For any fixed value of $x$ , `almost all’ the probability is within 2 standard deviations of the mean.

This is meant to illustrate the “rule of thumb” that “almost all” the probability is within 2 standard deviations of the mean. You should compare this with the graphs of the PDFs $f_{x}(t)$ for various (fixed) values of $x$ . The animation of these graphs as $x$ increases from 0 to 120 is shown below (when $A=.0001, B=.0003$ , and $c=1.07$ ).

The graph of the PDF $f_{x}(t)$ as $x$ increases from 0 to 120. The values of the parameters are fixed at $A=.0001, B=.0003$ , and $c=1.07$ . Compare this with the graphs of $\stackrel{\circ}e_{x}$ and *$\stackrel{\circ}e_{x}\pm 2\sigma_{x}$* above.

The graph of the PDF $f_{x}(t)$ as $x$ increases from 0 to 120. The values of the parameters are fixed at $A=.0001, B=.0003$ , and $c=1.07$ . Compare this with the graphs of $\stackrel{\circ}e_{x}$ and *$\stackrel{\circ}e_{x}\pm 2\sigma_{x}$* above.