2 Makeham–Beard Mortality Laws

Richards (Reference Richards2008) defines the Makeham–Beard family of laws to be those that can be written as:

(1) \begin{align} {\mu }_{x}=\frac{{\text{e}}^{\epsilon}+{\text{e}}^{\alpha +\beta x}}{1+{\text{e}}^{\alpha +\rho +\beta x}}\end{align}

dwhere $ {\mu }_{x}$ is the force of mortality for a life aged x, with parameters $ \epsilon\lt0,$ $ \alpha \lt0$ , $ \rho \gt0$ and $ \beta \gt 0.$

Richards (Reference Richards2008, Reference Richards2012) finds that the Makeham–Beard model provides the best statistical fit to the mortality rates of life office pensioners and defined benefit pension scheme members in the UK, when compared with other well-known parametric models.

The $ {e}^{\epsilon}$ term in equation (1) can be interpreted as representing non-age-related hazards (Makeham, Reference Makeham1860). The $ \rho $ term is the “heterogeneity” term introduced by Beard (Reference Beard1959) and can be interpreted as allowing for a mix of initial frailties in the population being considered. The $ \alpha $ and $ \beta $ terms are from the Gompertz (Reference Gompertz1825) model of mortality and set the initial “baseline” rate of mortality and the rate at which it increases with age, respectively. The overall logistic form that accommodates the decelerating increases in mortality at higher ages is due to Perks (Reference Perks1932).

The related integrated hazard function is given by (see Richards, Reference Richards2008)

(2) \begin{align} {H}_{x}^{MB}(t)=\underset{0}{\overset{t}{{\displaystyle \int }}}{\mu }_{x+s} ds=t {e}^{\epsilon}+\frac{{e}^{-\rho }-{e}^{\epsilon}}{\beta }\mathrm{log}\left(\frac{1+{e}^{\alpha +\rho +\beta \left(x+t\right)}}{1+{e}^{\alpha +\rho +\beta x}}\right)\end{align}

and hence the survival function for a life aged x is given by

(3) \begin{align} {S}^{MB}\left(x,t\right)=\mathrm{exp}\left(-{H}_{x}^{MB}(t)\right)={e}^{-t{e}^{\epsilon}}.{\left(\frac{1+{e}^{\alpha +\rho +\beta \left(x+t\right)}}{1+{e}^{\alpha +\rho +\beta x}}\right)}^{-\frac{{e}^{-\rho }-{e}^{\epsilon}}{\beta }}\end{align}

A continuously payable life annuity issued to a life aged x valued at a continuously compounded rate of interest δ (assumed known and constant) and subject to the Makeham–Beard mortality law with a terminal age denoted $ \omega $ is therefore obtained from

(4) \begin{align} {\bar{a}}_{x}^{MB}=\underset{0}{\overset{\omega -x}{{\displaystyle \int }}}{S}^{MB}\left(x,t\right){e}^{-\delta t}dt=\underset{0}{\overset{\omega -x}{{\displaystyle \int }}}{e}^{-t\left({e}^{\epsilon}+\delta \right)}.{\left(\frac{1+{e}^{\alpha +\rho +\beta \left(x+t\right)}}{1+{e}^{\alpha +\rho +\beta x}}\right)}^{-\frac{{e}^{-\rho }-{e}^{\epsilon}}{\beta }}dt\end{align}

The value of the annuity therefore depends on only a small number of parameters, which have some interpretation as mentioned earlier. Provided the integration in equation (4) can be performed quickly, the values of annuities with different combinations of parameter values can be easily and quickly obtained.

In a similar fashion, the value of an annuity written on joint lives x and y (assuming without loss of generality that y < x) subject to the same Makeham–Beard law is

(5) \begin{align} {\bar{a}}_{x,y}^{MB}&=\underset{0}{\overset{{T}_{x}}{{\displaystyle \int }}}{S}^{MB}\left(x,t\right){S}^{MB}\left(y,t\right){e}^{-\delta t}dt\nonumber\\[3pt] &=\underset{0}{\overset{{T}_{x}}{{\displaystyle \int }}}{e}^{-t\left(2{e}^{\epsilon}+\delta \right)}.{\left[\left(\frac{1+{e}^{\alpha +\rho +\beta \left(x+t\right)}}{1+{e}^{\alpha +\rho +\beta x}}\right)\left(\frac{1+{e}^{\alpha +\rho +\beta \left(y+t\right)}}{1+{e}^{\alpha +\rho +\beta y}}\right)\right]}^{-\frac{{e}^{-\rho }-{e}^{\epsilon}}{\beta }}dt\end{align}

The evaluation of the annuities in equations (4) and (5) can be done either by finding analytic expressions for the integrals (see later in Section 4) or by alternative methods such as discrete-time approximations or numerical integration.

The discrete-time approximations follow from observing that equation (4) can be written as

(6) \begin{align} {\bar{a}}_{x}^{MB}=\underset{m\to \infty }{\mathrm{lim}}\sum\limits_{t=0}^{m{T}_{x}-1}\frac{1}{m}.{S}^{MB}\left(x,t/m\right).{e}^{-\frac{\delta t}{m}}\end{align}

and hence a discrete-time approximation (for a given m) is

(7) \begin{align} {\bar{a}}_{x}^{MB}\approx \sum\limits_{t=0}^{m{T}_{x}-1}\frac{1}{m}.{S}^{MB}\left(x,(t+0.5)/m\right).{e}^{-\frac{\delta \left(t+0.5\right)}{m}}\end{align}

In principle, the approximation can be made arbitrarily good by increasing the size of m. In practice, the choice of m is often made by reference to a conventional time interval, e.g., a choice of m = 1 represents a mid-year approximation based on annual payments or m = 12 for monthly or m = 52 for weekly.

As an aside, most annuities are in practice paid at discrete intervals, so it would be perhaps more accurate to refer to the continuous formulation as the approximation. However, since we use the continuous-time value as the benchmark in this comparison, it is convenient to continue to refer to the discrete-time calculations as the approximations.

A similar approximation for the joint life annuity is

(8) \begin{align} {\bar{a}}_{x,y}^{MB}\approx \sum\limits_{t=0}^{m{T}_{x}-1}\frac{1}{m}.{S}^{MB}\left(x,(t+0.5)/m\right).{S}^{MB}\left(y,(t+0.5)/m\right).{e}^{-\frac{\delta \left(t+0.5\right)}{m}}\end{align}

The relevant analytic expressions are derived in Section 4 and depend on functions from a broad class known as the hypergeometric functions. A brief introduction to the two specific functions needed, namely, the Gauss hypergeometric function and Appell function of the first kind, is included in the next Section, along with a basic proposition and corollary based on standard properties of the functions that are used in the later derivations.

3 Gauss Hypergeometric Function and Appell Function of the First Kind

The Gauss hypergeometric function (see, e.g., Abramowitz & Stegun, Reference Abramowitz and Stegun1972) is defined as:

(9) \begin{align} {}_{2}F{}_{1}\left(a,b;\,c;\,z\right)=\sum\limits_{j=0}^{\infty }\frac{{\left(a\right)}_{j}{\left(b\right)}_{j}}{{\left(c\right)}_{j}}\frac{{z}^{j}}{j!}\end{align}

where $ {\left(n\right)}_{j}=n\left(n+1\right)\dots \left(n+j-1\right),$ c is not zero or a negative integer and the function converges for arbitrary a, b and c if $ \left|z\right|\lt1,$ or at $ \left|z\right|=1$ if the real part of $[c-a-b]>1$ .

An extension of the Gauss hypergeometric function for two variables is known as the Appell function (Appell, Reference Appell1925). There are four series of such extensions, the first of which is

(10) \begin{align} {F}_{1}\left(a;\,b,b^{\prime };\,c;\,z,w\right)=\sum\limits_{m=0}^{\infty }\sum\limits_{j=0}^{\infty }\frac{{\left(a\right)}_{m+j}{\left(b\right)}_{m}{\left(b^{\prime }\right)}_{j}}{{\left(c\right)}_{m+j}}\frac{{z}^{m}}{m!}\frac{{w}^{j}}{j!}\end{align}

The relationship between the Gauss hypergeometric function and the Appell function, F ₁, is given by (Schlosser, Reference Schlosser2013)

\begin{align*} {F}_{1}\left(a;\,b,0;\,c;\,z,w\right)={}_{2}{}^{}F{}_{1}\left(a,b;\,c;\,z\right)\end{align*}

In general, z and w can be complex, and the functions will converge if the arguments lie within the unit circle. For the application relevant to this note, the arguments are real functions of an underlying parameter, t, and are positive but less than one.

A further well-known result (see Schlosser, Reference Schlosser2013) relating to the Appell and Gauss hypergeometric functions is also used in obtaining the analytic results that follow

(11) \begin{align} {F}_{1}\left(a;\,b,b^{\prime };\,c;\,z,w\right)={\left(1-z\right)}^{-b}{\left(1-w\right)}^{-{b}^{\text{'}}}{F}_{1}\left(c-a;\,b,{b}^{\prime };\,c;\,\frac{z}{z-1},\frac{w}{w-1}\right)\end{align}

and in the special case $ {b}^{\prime }=0,$ the transformation is sometimes known as the Pfaff–Kummer transformation:

(12) \begin{align} {}_{2}F{}_{1}\left(a,b;\,c;\,z\right)={\left(1-z\right)}^{-b}{}_{2}F{}_{1}\left(c-a,b;\,c;\,\frac{z}{z-1}\right)\end{align}

Using the standard properties of Appell functions, the following proposition that we use later in Section 4 can be obtained.

Proposition 3.1 Let $ g, z$ and $ w$ be functions of t such that $ {g}^{\prime }(t)=g(t).{h}_{g},$ $ {z}^{\prime }(t)=z(t).h$ and $ {w}^{\prime }(t)=w(t).h$ and if $ \left(1-c\right)={h}_{g}/h$ then:

(13) \begin{align} \frac{d}{dt}{F}_{1}\left(c-1;\,b,b^{\prime };\,c;\,z(t),w(t)\right)\!.g(t)={g}^{\prime }(t).{\left(1-z(t)\right)}^{-b}{\left(1-w(t)\right)}^{{b}^{\prime }}\end{align}

Proof. See Appendix A.

A special case of Proposition 3.1 that is used in deriving the expression for a single life annuity in Section 4 follows as a corollary to Proposition 3.1.

Corollary 3.1 Let $ g$ and z be functions of t such that $ {g}^{\prime }(t)=g(t).{\text{h}}_{\text{g}}$ and $ {z}^{\prime }(t)=z(t).h$ and if $ \left(1-c\right)={h}_{g}/h$ then:

(14) \begin{align} \frac{d}{dt}{}_{2}{}^{}F{}_{1}\left(c-1,b;\,c;\,z(t)\right)\!.g(t)=g^{\prime}(t){\left(1-z(t)\right)}^{-b}\end{align}

Proof. Set $ {b}^{\prime }=0$ in equation (13) and use the fact noted earlier that the Gauss hypergeometric function is a special case of the Appell function.

Remark: The results in Proposition 3.1 and Corollary 3.1, and hence of Propositions 4.1 and 4.2 below, can also be obtained from symbolic mathematics software. For instance, the Mathematica code (Wolfram, 2019) for obtaining equation (15) is given in Appendix C. However, it is instructive to give the explicit proofs in Section 4.

4 Analytic Expressions for Single Life and Joint Life Annuities

The proposition and corollary obtained in Section 3 enable us to derive analytic expressions for evaluating single and joint life annuities where the underlying lives are subject to a Makeham–Beard mortality law.

Proposition 4.1 The value of a single life annuity issued to a life aged x can be evaluated using

(15) \begin{align} {\bar{a}}_{x}^{MB}={\left({\text{e}}^{\epsilon}+\delta \right)}^{-1}\left({}_{2}{}^{}F{}_{1}\left(1,b;\,c;\,{\mu }_{x}^{\ast}\right)-{e}^{-\delta {T}_{x}}.{S}^{MB}\left(x,{T}_{x}\right).{}_{2}{}^{}F{}_{1}\left(1,b;\,c;\,{\mu }_{\omega }^{\ast}\right)\right)\end{align}

where $ {\mu }_{x}^{\ast}=\frac{{\text{e}}^{\alpha +\rho +\beta x}}{1+{\text{e}}^{\alpha +\rho +\beta x}},$ $ {T}_{x}=\omega -x$ , $ b=\frac{{\text{e}}^{-{\rho}}-{\text{e}}^{\epsilon}}{\beta}$ and $ \left(1-c\right)=\frac{{\delta}+{\text{e}}^{\epsilon}}{\beta }.$

Proof. Substitute $ g(t)=-\frac{{\text{e}}^{-t\left({\text{e}}^{\epsilon}+{\delta}\right)}}{{\text{e}}^{\epsilon}+{\delta}}$ and $ z(t)=-{\text{e}}^{\alpha +\rho +\beta \left(x+t\right)}$ into equation (4) and rearrange to obtain

(16) \begin{align} {\bar{a}}_{x}^{MB}={\left(1-z\left(0\right)\right)}^{b}\underset{0}{\overset{{T}_{x}}{{\displaystyle \int }}}g^{\prime}(t).{\left(1-z(t)\right)}^{-b}dt\end{align}

Note that g and z satisfy the conditions of Corollary 3.1 and so, using equation (14):

(17) \begin{align} {\bar{a}}_{x}^{MB}={\left(1-z\left(0\right)\right)}^{b}\left({}_{2}{}^{}F{}_{1}\left(\text{c}-1,b;\,c;\,z\left(0\right)\right)\!.g\left(0\right)-{}_{2}{}^{}F{}_{1}\left(\text{c}-1,b;\,c;\,z\left({T}_{x}\right)\right)\!.g\left({T}_{x}\right)\right)\end{align}

Apply the transformation in equation (12) to equation (17) and simplify, noting that $ {\mu }_{x+t}^{\ast}=\frac{z(t)}{z(t)-1}:$

(18) \begin{align} {\bar{a}}_{x}^{MB}=g\left(0\right).\left({}_{2}{}^{}F{}_{1}\left(1,b;\,c;\,{\mu }_{x}^{\ast}\right)-{}_{2}{}^{}F{}_{1}\left(1,b;\,c;\,{\mu }_{\omega }^{\ast}\right).\frac{g\left({T}_{x}\right)}{g\left(0\right)}.\frac{{\left(1-z\left(0\right)\right)}^{b}}{{\left(1-z\left({T}_{x}\right)\right)}^{b}}\right)\end{align}

Further noting that $ {S}^{MB}\left(x,t\right)$ defined in Section 2 can be written as $ {\text{ e}}^{-t{\text{e}}^{\epsilon}}.{\left(\frac{1-\text{z}\left(\text{t}\right)}{1-\text{z}\left(0\right)}\right)}^{-b},$ equation (18) simplifies to the required result.

Incidentally, the expression for $ {\mu }_{x}^{\ast}$ can be recognised as being the force of mortality at age x suggested by Perks (Reference Perks1932) with parameters $ \alpha +\rho $ and $ \beta.$ The argument, $ {\mu }_{x}^{\ast}$ , is less than one for all x and so the hypergeometric functions will converge (subject to c not being a negative integer).

Proposition 4.2 The value of a joint life annuity issued to lives aged x and y with x > y can be evaluated using

(19) \begin{align} {\bar{a}}_{x,y}^{MB}=\frac{1}{\left(2{\text{e}}^{\epsilon}+\delta \right)}\left(F\left({\mu }_{x}^{\ast},{\mu }_{y}^{\ast}\right)-{e}^{-\delta {T}_{x}}.{S}^{MB}\left(x,{T}_{x}\right).{S}^{MB}\left(y,{T}_{x}\right).F\left({\mu }_{\omega }^{\ast},{\mu }_{y+{T}_{x}}^{\ast}\right)\right)\end{align}

where $ {\mu }_{x}^{\ast}=\frac{{\text{e}}^{\alpha +\rho +\beta x}}{1+{\text{e}}^{\alpha +\rho +\beta x}},$ $ {T}_{x}=\omega -x$ , $ b=b^{\prime }=\frac{{\text{e}}^{-{\rho}}-{\text{e}}^{\epsilon}}{\beta }$ , $ \left(1-c\right)=\frac{{\delta}+2{\text{e}}^{\epsilon}}{\beta }$ and defining $ F\left(z,w\right)={F}_{1}\left(1;\,b,{b}^{\prime };\,\text{c};\,z,w\right)$

Proof. Substitute $ g(t)=-\frac{{\text{e}}^{-t\left(2{\text{e}}^{\epsilon}+{\delta}\right)}}{2{\text{e}}^{\epsilon}+{\delta}},$ $ z(t)=-{\text{e}}^{\alpha +\rho +\beta \left(x+t\right)}$ , $ w(t)=-{\text{e}}^{\alpha +\rho +\beta \left(y+t\right)}$ into equation (5) and rearrange to obtain

(20) \begin{align} {\bar{a}}_{x,y}^{MB}={\left(1-z\left(0\right)\right)}^{b}{\left(1-w\left(0\right)\right)}^{{b}^{\prime }}.\underset{0}{\overset{{T}_{x}}{{\displaystyle \int }}}{g}^{\prime }(t).{\left(1-z(t)\right)}^{-b}.{\left(1-w(t)\right)}^{-{b}^{\prime }}dt\end{align}

The proof then follows along the same lines as Proposition 4.1, but using Proposition 3.1 instead of Corollary 3.1 and transformation (11) instead of (12).

5 Numerical Comparison of Evaluation Times

The formulae for the annuities in Propositions 4.1 and 4.2 are analytic in that they are expressed in terms of convergent series, i.e., the hypergeometric and Appell functions. However, because those functions involve an explicit summation of multiple terms, it is relevant to understand if their evaluation is quicker in practice than alternative numerical or approximate methodologies.

The first pair of alternatives is the discrete-time approximations described in equations (7) and (8) in Section 2. The second alternative method is numerical integration to evaluate the integrals in equations (4) and (5) directly.

To test the relative speeds of the different evaluation methodologies, we set up a simple simulation wherein we generate randomly a large number of parameter combinations and then evaluate the annuities associated with them. The total time required to perform the evaluations is used to compare the effective speeds of the methodologies.

In fact the elapsed time is a measure of both the methodology and the way in which it is implemented in code, including, e.g., the programming language and the hardware used. It is therefore important to try to keep the implementation details as consistent as possible across the methodologies. The comparisons shown in Table 1 are all implemented in R, except for a much faster implementation of numerical integration based on code compiled in C. More details of how the key functions are coded and of the computing environment are set out in Appendix B.

Table 1.

Comparison of the times taken to evaluate annuities using different approaches.

Note: Times are shown in seconds for the evaluation of 1,000,000 annuities with different discount rates randomly drawn between 0 and 0.1. The mean absolute deviations (MAD, calculated relative to the relevant analytic values) are also shown. All the functions are coded in R (denoted (R) in the column headings) apart from the second version of the numerical integration approach that uses compiled C code. The values of m for the discrete-time approximations refer to the parameter m in equations (7) and (8). Total elapsed time is measured using the system.time function in R.

Note: Times are shown in seconds for the evaluation of 1,000,000 annuities with different discount rates randomly drawn between 0 and 0.1. The mean absolute deviations (MAD, calculated relative to the relevant analytic values) are also shown. All the functions are coded in R (denoted (R) in the column headings) apart from the second version of the numerical integration approach that uses compiled C code. The values of m for the discrete-time approximations refer to the parameter m in equations (7) and (8). Total elapsed time is measured using the system.time function in R.

5.1 Results of numerical comparison

Table 1 shows the elapsed times (in seconds) for 1,000,000 evaluations of annuities at a selection of different ages using fixed mortality parameters (see Appendix B for the values chosen) but randomly drawn discount rates (δ) between 0 and 0.1. In principle, all the parameters could be varied jointly but that does not affect the timings.

As a measure of consistency in the values obtained across the methodologies, the mean absolute differences (MAD) expressed relative to the analytic approach in the evaluations are also shown. Small values of MAD indicate that the values of the annuities from the two approaches are very similar; larger values indicate that there are some differences. The calculation of the MAD for approach j is given by

(21) \begin{align} {\text{MAD}}_{\text{j}}=\frac{1}{\text{N}}\sum\limits_{i=1}^{N}\left|{A}_{j}(i)-{A}_{0}(i)\right|\end{align}

where j denotes one of the discrete-time approximations or one of the numerical integration approaches; N is the number of evaluations (1,000,000 in this case); $ {A}_{j}(i)$ is the value obtained for the annuity in evaluation i for approach j and $ {A}_{0}(i)$ is the value from the analytic approach.

In all cases, the analytic expressions execute most quickly. The annual discrete approximation (m = 1) is the next quickest, taking around three times as long for the single life case, but broadly the same time for the older joint life annuities. The monthly versions of the discrete-time approximation take around ten times longer than the annual version because there are around 12 times as many terms in the summation. For a similar reason, the times for the younger life (55) take longer than those for the older life (85).

For the single life annuities, the compiled implementation of numerical integration is the next fastest and takes around five times as long, but the version coded in R is much slower, taking ten times as long again as the compiled version.

The MAD values indicate that there is little or no practical difference in the values of the annuities calculated using the different approaches, with the possible exception of the annual discrete-time approximation which is less accurate.

For the joint life annuities, the analytic approach is also the fastest. The compiled numerical integration method takes about twice as long depending on the age of the annuitant. The analytic approach is relatively slower for higher ages as convergence of the Appell function using the implementation in Appendix B takes longer. The use of a compiled language for the summations is likely to speed up the evaluation of the Appell functions and discrete-time approximations considerably.

Table 2.

Comparison of analytic and numerical integration using Mathematica.

Note: Times are shown in seconds for the evaluation of 10,000 annuities with different discount rates randomly drawn between 0 and 0.1. All the functions are as standard in Mathematica.

A further check on the robustness of the timing comparisons was obtained by undertaking a similar experiment using Mathematica (Wolfram, 2019). Evaluations were implemented for the single life annuities for ages 55 and 85, and for the joint life annuity for a pair of lives aged 85 and 82 years; the joint life annuity for a pair of ages 65 and 55 could not be evaluated using the same code. The results are summarised in Table 2, and the code used can be found in Appendix B.

Broadly in keeping with the results using R, the numerical integration approach is over 80 times as long as the analytic approach for the single life annuities when using the standard implementation of these functions in Mathematica. For the joint life annuity that was able to be evaluated using the same code, the times for the two approaches are very similar.

6 Concluding Remarks

In this note, expressions involving the Gauss hypergeometric function and Appell function of the first kind have been derived to evaluate analytically annuities based on a general family of mortality models (the so-called Makeham–Beard family) that incorporate decelerating increases in mortality at the higher ages.

Hypergeometric functions are readily available in many mathematical packages used by actuaries or can be coded relatively simply, making the expressions straightforward to use in many cases. However, it should be noted that some implementations of the functions in mathematical packages do not accommodate all the possible combinations of inputs or limiting results as standard, so some care may be needed.

A comparison of evaluation times based on the analytic expressions coded in R show that numerical integration, even when compiled in C, takes twice (for joint life) to five times (for single life) as long as the analytic calculations. Numerical integration coded in R takes 15 to 60 times as long. Discrete-time approximations, coded in R with comparable levels of accuracy to the analytic approach, take 10 to 30 times longer than the analytic calculations.

A more limited comparison using Mathematica confirms that evaluation using numerical integration is around 80 times as long as the analytic approach for single life annuities. In Mathematica, there is little difference in evaluation times for the joint life annuity although it is difficult to investigate these results in more detail given the proprietary nature of that software.

For applications that require only a small number of, or only occasional, evaluations, the absolute lengths of time may not be an issue as the average times per calculation for all methods are considerably less than a second. However, the analytic approaches offer material improvements if a very large number of calculations are required, or if speed is of the essence.