1 Introduction
Analytic expressions for the value of an annuity in terms of the relatively small number of parameters that define a parametric mortality “law” can be convenient for several practical reasons. For example, if large numbers of evaluations are required very quickly for interactive purposes or in simulations, analytic approaches may offer a speed advantage over methodologies such as numerical integration or other approximations.
Comparable cases where expressions for annuities based on parametric mortality laws have been used include Milevsky et al. (Reference Milevsky, Salisbury and Chigodaev2016), who derived an expression for single life annuities based on the Gompertz–Makeham family of mortality laws, and Seal (Reference Seal1962), who derived expressions for joint life annuities based on the Makeham law.
In this note, analytic expressions are derived for the values of single life and joint life annuities based on a more general class of parametric mortality laws than those above, namely, the Makeham–Beard family as defined by Richards (Reference Richards2008, Reference Richards2012).
The High Age Working Party of the CMI (CMI Working Paper 100, 2017) and Richards (Reference Richards2008) show that laws that incorporate a logistic function at higher ages, i.e., featuring a deceleration in the rate of increase in mortality for the higher ages, provide a better fit to UK pensioner data than the Gompertz–Makeham model that implies exponential increases in mortality rates. Whelan (Reference Whelan2009) reports similar findings for Irish data.
Various models that include the logistic function and hence the feature of decelerating increases in mortality rates at higher ages have been proposed over the years starting with Perks (Reference Perks1932); the Makeham–Beard family as defined by Richards (Reference Richards2008) includes many of these as special cases.
The expressions for single and joint life annuities that are derived in Section 4 turn out to involve the Gauss hypergeometric function and a two-variable generalisation, the Appell function of the first kind, respectively. These functions belong to a well-studied class of functions that may not be very familiar to some practitioners even though they have appeared in the actuarial literature in similar contexts.
For example, the expression derived by Seal (Reference Seal1962) involves the confluent hypergeometric function, which can be interpreted as a special case of the Gauss hypergeometric function. An analytic expression for life expectancy using a gamma-Gompertz multiplicative frailty model investigated by Missov (Reference Missov2013) also contains the Gauss hypergeometric function.
In Section 2, a summary of the Makeham–Beard family of mortality laws as defined by Richards (Reference Richards2008) is set out along with integral expressions for single life and joint lives annuities. In Section 3, the definitions of the Gauss hypergeometric function and Appell function of the first kind are given along with a proposition and corollary that are based on standard results relating to them. In Section 4, the propositions from Section 3 are applied to derive the main results of the note, namely, the analytic expressions for the values of single and joint annuities. A numerical comparison of the times required to evaluate the analytic expressions as well as other common approximate and numerical methods is reported in Section 5. Concluding remarks are given in Section 6.
2 Makeham–Beard Mortality Laws
Richards (Reference Richards2008) defines the Makeham–Beard family of laws to be those that can be written as:
\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)
\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
\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
\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
\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
\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
\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
\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:
\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
\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 1, 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
\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:
\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:
\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:
\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
\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
\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):
\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}:$
\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
\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
\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
\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.
Acknowledgements
The author is grateful for very helpful and detailed suggestions from an anonymous reviewer and the associate editor.
Appendix A
Proof of Proposition 3.1 The proof depends on three sets of standard results relating to Appell functions (see Schlosser, Reference Schlosser2013), including a special case, the derivatives of the Appell function and a so-called contiguous relationship:
Standard result 1 (special case when a = c):
\begin{align} {F}_{1}\left(\text{c};\,b,{b}^{\prime };\,\text{c};\,\text{z},\text{w}\right)&=\sum\limits_{m=0}^{\infty }{\left(b\right)}_{m}\frac{{z}^{m}}{m!}\sum\limits_{j=0}^{\infty \text{ }}{\left({b}^{\prime }\right)}_{j}\frac{{w}^{j}}{j!}={}_{2}{}^{}F{}_{1}\left(1,b;\,1;\,z\right).{}_{2}{}^{}F{}_{1}\left(1,{b}^{\prime };\,1;\,w\right)\nonumber\\[3pt] &={\left(1-z\right)}^{-b}{\left(1-w\right)}^{-{b}^{\prime}}\end{align}
Standard result 2 (derivatives of the Appell function with respect to z and w):
\begin{align} \frac{\partial {F}_{1}\left(\text{a};\,b,b^{\prime };\,{c};\,{z},{w}\right)}{\partial z}=\frac{ab}{c}{F}_{1}\left(a+1;\,b+1,{b}^{\prime };\,c+1;\,z,w\right)\end{align}
\begin{align} \frac{\partial {F}_{1}\left(a;\,b,b^{\prime };\,{c};\,{z},{w}\right)}{\partial w}=\frac{a{b}^{\prime }}{c}{F}_{1}\left(a+1;\,b,{b}^{\prime }+1;\,c+1;\,z,w\right)\end{align}
Standard result 3 (so-called “contiguous” relationship):
\begin{align} {F}_{1}\!\left(a+1;\,b,{b}^{\prime };\,{c};\,{z},{w}\right)&={F}_{1}\left(a;\,b,{b}^{\prime };\,{c};\,{z},{w}\right)+\frac{bz}{c}{F}_{1}\left(a+1;\,b+1,{b}^{\prime };\,c+1;\,z,w\right)\nonumber\\[4pt] &\quad +\frac{{b}^{\prime }w}{c}{F}_{1}\left(a+1;\,b,{b}^{\prime }+1;\,c+1;\,z,w\right)\end{align}
The proof follows by differentiating
$ {F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\!.g(t)$
with respect to t and then making use of the standard results above to simplify, i.e.,
\begin{align} \frac{d}{dt}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\!.g(t)&={g}^{\prime }(t).{F}_{1}\left(\text{a};\,b,{b}^{\prime };\,c;\,z,w\right)+g(t).\frac{\partial }{\partial z}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right).{z}^{\prime }(t)\nonumber\\[3pt] &\quad +g(t).\frac{\partial }{\partial w}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right).{w}^{\prime }(t)\end{align}
Then substitute the conditions
$ {z}^{\prime }(t)=z(t).h,$
$ {w}^{\prime }(t)=w(t).h\text{ }$
and
$ {g}^{\prime }(t)=g(t).{{h}}_{g}$
into equation A5 to obtain
\begin{align} \frac{d}{dt}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\!.g(t)=g(t)\left[\begin{array}{c}{{h}}_{g}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\\[3pt] \!\!\!+h.z(t).\frac{ab}{c}{F}_{1}\left(a+1;\,b+1,{b}^{\prime };\,c+1;\,z,w\right)\\[3pt] +h.w(t).\frac{a{b}^{\prime }}{c}{F}_{1}\left(a+1;\,b,{b}^{\prime }+1;\,c+1;\,z,w\right)\end{array}\right]\end{align}
Using the standard results in equations A2, A3 and A4 implies
\begin{align} \frac{d}{dt}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\!.g(t)&=g(t)\left[\begin{array}{c}{\text{h}}_{g}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\\ +ah\left\{\begin{array}{c}{F}_{1}\left(a+1;\,b,{b}^{\prime };\,c;\,z,w\right)\\ -{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\end{array}\right\}\end{array}\right]\nonumber\\[3pt] &={h}_{g}.g(t)\left[\begin{array}{c}{F}_{1}\left(a;\,b,{b}^{\prime };\,c;\,z,w\right)\left(1-\frac{ah}{{h}_{g}}\right)\\ +\frac{ah}{{h}_{g}}{F}_{1}\left(a+1;\,b,{b}^{\prime };\,c;\,z,w\right)\end{array}\right]\end{align}
Putting
$ a=\left(c-1\right)=\frac{{h}_{g}}{h}$
and applying the special case in equation A1 completes the proof:
\begin{align} \frac{d}{dt}{F}_{1}\left(c-1;\,b,b^{\prime };\,c;\,z,w\right)\!.g(t)={g}^{\prime }(t).{\left(1-z\right)}^{-b}{\left(1-w\right)}^{-{b}^{\text{'}}}\end{align}
Appendix B
This appendix contains some further information about the functions used in the numerical comparison, including the R and Mathematica code for some key functions.
Implementation of Gauss hypergeometric and Appell functions in R
Functions for evaluating general Gauss hypergeometric and Appell functions in R can be obtained from several different “packages” including, e.g., the appell package (Bove, Reference Bove2013). These and similar packages are designed to deal with many of the complications and continuations required for the most general set of parameters, including dealing with the fact that they may be complex. However, for the application in this note where the parameters are all real and converge, it is quicker to use a straightforward coding of the series and to use limiting cases when appropriate. Even in this simpler context, there are some special cases where care is needed.
We set out below two functions that deal with the Gauss hypergeometric and Appell functions for cases where z and w (when appropriate) are much less than one, and two functions that deal with the limiting case where z is very close to one.
The implementation of the Gauss hypergeometric function based on equation (9) is . Please note that the implementation shown here also does not check for some special cases, such as c being a negative integer. The number of terms used to check for convergence in this implementation depends on only z; values of c that are close to (but not exactly) a negative integer are also likely to mean that this implementation is not accurate.
If we note the definition of c in Proposition 4.1, we see that c could be close to a negative integer if
$ \delta \approx n\beta -{e}^{\epsilon}$
where n is a positive integer. However, for the choices of the Makeham–Beard parameters (see later) and discount rate in Section 5, this is not a concern for this application.
If the terminal age (denoted
$ \omega $
and chosen to be 250 in the comparison in Section 5) is very high and hence
$ {\mu }_{\omega }^{\ast}$
in equation (15) is very close to unity, evaluation is quicker if we use an approximation derived from a series of transformations (see Schlosser, Reference Schlosser2013). When the parameters are such that c − a − b < 0, the relevant approximation is given in f21limit.
An implementation of the Appell function (for real z and w that are both less than one, with a = 1 and b = b from equation (10)) is set out in .
When z is close to one, the Appell function is approximated by the limiting value that is obtained through another series of standard transformations and reductions (see Schlosser, Reference Schlosser2013) and is set out in . The function calls the Gauss hypergeometric function (f21hyperR) and another function (f21limitP) that evaluates the hypergeometric function when z = 1 and the parameters satisfy the condition that c – a – b > 0.
Numerical integration
The numerical integration routines used are the integral function from the pracma package (Borchers, Reference Borchers2019) and the integrate function from the stats package (R Core Team, 2018). Both these functions use Gauss–Kronrod quadrature. The function integral is implemented in R, whereas the function is a R wrapper that calls the quadrature routines compiled in C and is therefore much faster.
Makeham–Beard parameters
The parameters used in the numerical analysis are
$ \epsilon=-5.07444,$
$ \alpha =-13.2666$
,
$ \rho =0.0931624$
and
$ \beta = 0.131902.$
These values are taken from Richards (Reference Richards2019) and represent the maximum likelihood estimates based on 125,535 male pensioners in England and Wales (31,207 deaths).
Processing environment
R version 3.5.2 (R Core Team, 2018) was used and run on a laptop with an Intel(R) Core™ i5-5300U CPU @ 2.30 GHz processor with 8.00 GB RAM installed and running Windows 10 Enterprise N.
Mathematica code
The code used for evaluating the single life annuities (shown here for age 55) using the analytic approach is
The code used for evaluating the single life annuities (shown here for age 55) using the numerical integration approach is
The code for the analytic evaluation of the joint life annuity (shown for the pair of ages 85 and 82) is:
Appendix C
The following Mathematica code yields a (slightly differently parameterized) version of equation (15) in terms of the Mathematica Hypergeometric2F1 function.
The result of this run is:
\begin{align*} &\frac{1}{\delta+e^{\epsilon}}\left(e^{\alpha+\rho+\beta{\textbf{x}}}+1\right)^{\frac{e^{-\rho}-e^{\epsilon}}{\beta}}\,\bigg(\texttt{Hypergeometric2F1}\,\bigg[\frac{e^{\epsilon}+e^{-\rho}}{\beta}, -\frac{\delta+e^{\epsilon}}{\beta}, -\frac{-\beta + e^{\epsilon}+\delta}{\beta},\\
&\quad -e^{\alpha+\textbf{x}\beta+\rho}\bigg]\, -e^{\left(\delta+e^{\epsilon}\right)(\textbf{x}-\omega)} \texttt{Hypergeometric2F1}\,\bigg[\frac{-e^{\epsilon}+e^{-\rho}}{\beta}, -\frac{\delta+e^{\epsilon}}{\beta}, -\frac{-\beta + e^{\epsilon}+ \delta}{\beta},\\
&\quad -e^{\alpha+\rho+\beta\omega}\bigg]\,\bigg)\, \texttt{if}\ \ e^{\alpha+\rho+\beta\textbf{x}} \geq -\texttt{1\&\&} e^{\alpha+\beta\omega+\rho} \geq -\texttt{1} \end{align*}
