2.1. Three shapes of mortality curves

Since a deep-deferred annuity does not start to pay benefits until the annuitant reaches a very advanced age, its cash flow may be greatly affected by the tail of a mortality curve. To examine the impact of mortality tail assumption on the annuity value, we consider the following three models:

(1) Gompertz law (Gompertz, Reference Gompertz1825):

where μ _x is the force of mortality for a life at age x, and a and b are constants. Under Gompertz law, the force of mortality increases exponentially with age. The mortality curve is concave upward without bound. Gompertz law has been used widely because of its simplicity and tractability.

(2) Cubic model with a flat tail (Panjer & Russo, Reference Panjer and Russo1992; Panjer & Tan, Reference Panjer and Tan1995):

where a, b, c, and d are constants. The force of mortality increases as a cubic function of age until it reaches the cap, ln2, and then remains at ln2 for higher ages. This mortality curve models a late-life mortality plateau. The recently released RP-2014 Mortality TablesFootnote ¹ by the Society of Actuaries (SoA) uses a similar mortality curve, and sets the death probability to 0.5 for males over 110-year-old and females over 112-year-old.

(3) Logistic model (Perks, Reference Perks1932)

where a, b, c, and d are constants. The force of mortality accelerates first, and then decelerates at higher ages, gradually approaching an asymptotic maximum. There are two variants of this model: Beard (Reference Beard1963), where e ^d =0; and Kannisto (Reference Kannisto1992), where e ^d =0 and b=c.

2.2. Data and fitting

We use Japanese mortality data in year 2012 with the age range of 75–105Footnote ² to fit the three mortality curves described above.

We recognise that the market for deep-deferred annuities mainly lies in the United States. We choose to use Japanese mortality data because there is no US annuitant mortality data publicly available. Our paper focuses on how the assumptions of mortality tail curve and mortality improvement at high ages affect deep-deferred annuity valuation. High age mortality data are crucial for our study. One alternative is to use the US national population data, but it is recognised that annuitant mortality is much lighter than national population mortality due to adverse selection (Mitchell & McCarthy, Reference Mitchell and McCarthy2002). Using the national population data may result in underestimation of demographic risk in deep-deferred annuities.

Japan is one of the countries that have the world’s lowest mortality rates (Robine et al., Reference Robine, Cheung, Saito, Jeune, Parker and Herrmann2010), and it has a relatively large population. Japanese mortality data are of high quality since Japan implemented a national family registration system in 1872. Researchers (Robine & Saito, Reference Robine and Saito2003) have found that Japanese mortality data offer some unique information on the demographic trends for extremely old population. Japanese data are also reliable for investigating the determinants and limits of human longevity and the important plasticity in longevity with regard to environmental conditions (Arai et al., Reference Arai, Inagaki, Takayama, Abe, Saito, Takebayashi and Hirose2014). These features of Japanese mortality data prompt us to use them in studying demographic risk in deep-deferred annuity valuation.

The number of deaths at each age is assumed to follow a Poisson distribution with mean equal to the expected number of deaths under a mortality curve model:

where D _x,t , E _x,t , and μ _x,t are the number of deaths, the number of exposures-to-risk, and the force of mortality of lives at age x in year t, respectively. Maximum likelihood estimation is used to estimate the parameters in the mortality curve models. The log-likelihood function can be expressed as

To maximise the log-likelihood, we can use standard Newton’s method or optimisation functions in computational softwares. We obtain our estimates using the optimisation function in Matlab. Table 1 summarises the parameter estimates.

We use the Bayesian information criteria (BIC)Footnote ³ to compare the goodness of fit of the three models. Perks’ logistic model has the lowest BIC for females and is thus the best fit for female mortality data. While Cubic model is the best fit for male mortality data, Perks’ model fits male mortality data slightly worse. We have also fitted Beard’s and Kannisto’s logistic models, and found that they do not yield better goodness of fit than Perks’ model. Therefore, we will only consider Perks’ logistic model hereinafter.

The upper panel of Figure 1 plots raw mortality rates and fitted curves for males and females. The three curves have similar fits for ages under 95; however, they show very different patterns for ages over 95. The fitted male mortality curve using Perks’ model is concave upward before age 97 and concave downward after age 97. It gradually approaches an asymptote of 0.686. The fitted female mortality curve using Perks’ model has a reflection point at age 100 and an asymptote of 0.660. The fitted curve using cubic model is concave upward for ages over 75. It reaches the imposed cap of ln2 at age 106 for males and age 109 for females, respectively.

Figure 1 (a) Raw mortality rates and the three fitted mortality curves for a male and female currently at age 75. (b) Probability density function of the future lifetime of a person aged 75 in 2012 using the three mortality curves.

The goodness of fit of Gompertz law is much less accurate than that of the other two models. The fitted curve deviates significantly from the raw mortality rates beyond age 95 for males and age 97 for females. However, Gavrilov & Gavrilova (Reference Gavrilov and Gavrilova2011) found that Gompertz law fits the US Social Security Administration (SSA) Death Master File mortality dataFootnote ⁴ of ages 88–106 better than logistic models. They concluded that improvement of data quality at advanced ages results in a mortality trajectory closer to that of the Gompertz law. Therefore, determining which model fits best appears to be very much data dependent. Practitioners should carefully choose the model for mortality tail curve based on the mortality experience of their own annuity portfolios. The large difference between the three fitted mortality curves suggests that mortality tail curve may be an important source of model risk for deep-deferred annuity valuation.

Denote f ₇₅(t) as the probability density function of the future lifetime of a person aged 75 in 2012, and f ₇₅(t)= _t p ₇₅ μ _75+t,2012, where _t p ₇₅ is the probability that this person survives t years. In the lower panel of Figure 1, we plot f ₇₅(t), t=25, 26, … , 45 for both males and females using the three mortality tail curves. Gompertz law has the lightest tail because it yields the lowest probability of surviving to very high ages among the three mortality tail curves. We also notice that the value of f ₇₅(t) for males is much smaller than that for females, because males have much lower probability of surviving to very high ages. In order to display the difference between the density functions using various mortality tail curves clearly, we use different scales for males and females in these two plots.

We construct the complete mortality curve with raw mortality rates for ages below 75 and fitted mortality rates for ages above 75 using the models described in section 2.1 because the raw mortality rates for younger ages are believed to be reliable. We use raw mortality rates for ages below 75 in the mortality curve to avoid fitting errors. In addition, using the same mortality rates for younger ages ensures that the difference of annuity values comes solely from the mortality tail curve assumption.

3.1. Single-population Lee–Carter model

In the main text, we use the Lee–Carter model to obtain mortality improvement rates, because it is often regarded as the benchmark for mortality modelling. In Appendix A, we introduce the Mitchell et al. (Reference Mitchell, Brockett, Mendoza-Arriaga and Muthuraman2013) model and show that these two models result in similar findings.

Mathematically, the Lee–Carter model can be expressed as follows:

where m _x,t is the central death rate for a life aged x in year t, α _x the age-specific mortality level, κ _t the time-varying mortality index, and β _x the age-specific sensitivity to the mortality index. Following Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002), we assume that D _x,t follows a Poisson distribution with mean equal to E _x,t m _x,t . Note that Poisson assumption is commonly violated by overdispersion, which is often evident in death count data due to population heterogeneity. Li et al. (Reference Li, Chan and Zhou2009) addressed this problem by adding a gamma distribution for the unobserved heterogeneity. For simplicity, we continue to use the classic Poisson distribution and maximum likelihood estimation to obtain parameter estimates. To ensure parameter uniqueness, two identification constraints are imposed: ∑ _x β _x =1 and ∑ _t κ _t =0. κ _t is further modelled by an autoregressive integrated moving average (ARIMA) model.

Mortality improvement rate at age x and time t can be written as follows:

Equation (3.1) indicates that mortality improvement rate for age x is determined by β _x and the change in κ _t . It is important to note that the cohort effect is another important factor driving mortality improvement due to the phenomenon that people born in certain years have experienced more rapid mortality improvement than others. Cohort effect has been observed in both United Kingdom and Japanese populations by Willets (2004). Several researchers (Renshaw & Haberman, Reference Renshaw and Haberman2003; Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) have incorporated cohort effect in mortality modelling, and it is considered an important factor in annuity pricing and valuation. Although the Lee–Carter model cannot capture cohort effect, it has been chosen because the focus of our study is on the uncertainty stemming from the assumption of mortality tail curve and the dependence structure of male and female mortality. Future research may consider a collection of mortality models with cohort effect in the analysis.

The forecast of mortality rates is the multiplication of the base year mortality curve and the projected mortality improvement rate. This two-step approach allows us to separately examine the impact of mortality improvement model and the base mortality curve on deep-deferred annuity valuation. In this paper, we use 2012 as the base year. For a life aged x at the beginning of 2013, his/her death rate in year 2013+s, where s=0, 1, 2, … , can be expressed as follows:

where μ _x+s,2012 is obtained from the base mortality curve in section 2. This approach is analogous to the mortality improvement scale approach used in the SoA’s Mortality Improvement Scale MP – 2015. In SoA’s approach, future mortality rates in year t+s are calculated by

where IS (x; s) is the improvement scale, retrieved from a two-dimensional table with one dimension for age and the other for calendar year. The base year mortality curve is the RP-2014 Mortality Table. However, this tabulated approach can only provide point forecasts of future mortality rates, but not the uncertainty of the forecasts. In contrast, our approach is stochastic and can capture the uncertainty in mortality improvement rates.

If we use the Lee–Carter model directly, the forecasted value of m _x+s,2013+s is

where $$\hat{m}_{{x{\plus}s,2012}} \, {\equals}\, e^{{\alpha _{{x{\plus}s}} {\plus}\beta _{{x{\plus}s}} \kappa _{{2012}} }} $$ and it is the fitted value of m _x+s,2012. Therefore, using the Lee–Carter model to directly forecast mortality rates is equivalent to our two-step approach with the base mortality curve set to the fitted Lee–Carter mortality rates in the base year.

We fit the Lee–Carter model to the Japanese male and female mortality data from the sample period of 1965–2012. Due to the concern of data reliability, we use the sample age range of 60–95. Data before year 1965 were dropped due to greater volatility compared to subsequent years. Figure 2 demonstrates the estimated values of α _x , β _x , and κ _t for males and females, respectively. α _x increases almost linearly with age, indicating that mortality rates increase with age. As shown in equation (3.1), mortality improvement rate is determined by β _x and κ _t −κ _t−1. Since κ _t decreases with t and β _x is positive, mortality rates decrease over time for all ages. β _x attains its maximum value at approximately age 75, thus suggesting that people who were approximately age 75 experienced faster mortality improvement than other age groups. We also observe that female κ _t has a steeper slope than male κ _t , and females over age 77 have higher β _x than males. This observation implies that the mortality improvement of females over age 77 is faster than that of males. Note that other countries may have different experiences. For example, the mortality of UK males has improved more rapidly than females since 1970, according to the 2012-based National Population Projections report by the Office for National Statistics (2013).

Figure 2 Estimates of α _x , β _x , and κ _t .

We selected ARIMA(0, 1, 1) model for both male and female κ _t based on BIC values. The ARIMA(0, 1, 1) model can be expressed as follows:

where Δκ _t =κ _t −κ _t−1 and {z _t } are independently and identically distributed normal random variables with mean 0 and variance σ ². The parameter estimates for the ARIMA model and their standard errors (in parentheses) are shown in Table 2. Depending on the choice of identifiability constraints, the estimates of κ _t may change, thereby changing the estimates of ARIMA parameters as well. However, the structure of the fitted ARIMA model and future morality forecast will remain the same.

Since the sample age range of the data is 60–95, we can only obtain the estimates of α _x and β _x up to age 95. For ages above 95, we extrapolate the patterns of α _x and β _x using deterministic functions. It is reasonable to assume that β _x is non-negative, because mortality is expected to continue to improve, or at least not to deteriorate. With this constraint, we consider two models for β _x : a Gaussian function with an asymptote of 0, and a quadratic function with a minimum value of 0. Using Gaussian function, β _x decreases slowly and approaches 0 as x goes to infinity. Using a quadratic function, β _x decreases quickly and reaches the floor of 0 when x=105. Both cases are displayed in Figure 3. Since α _x does not determine mortality improvement rate, as seen in equation (3.1), the choice of extrapolation function will not affect our conclusion. We can use either a Gaussian function or linear function with a cap of 0 to extrapolate α _x .

Figure 3 Extrapolation of α _x and β _x at high ages.

3.2. Two-population Lee–Carter model

To build a two-population model, we use one Lee–Carter model for each population:

where i=1 for males and 2 for females. Following Zhou et al. (Reference Zhou, Wang, Kaufhold, Li and Tan2014), we further model $$\left[ {\Delta \kappa _{t}^{{(1)}} ,\Delta \kappa _{t}^{{(2)}} } \right]\prime$$ by a pth order vector autoregressive (VAR) process:

where A ₀ is a 2×1 vector, A _j , for j=1, 2, … , p, is a 2×2 matrix, and $$[ {z_{t}^{{(1)}} ,z_{t}^{{(2)}} } ]\prime$$ follows a bivariate normal distribution with 0 mean and variance–covariance matrix of Σ.

We use likelihood ratio tests to determine the lag order p. Using a significance level of 5%, a p-value smaller than 5% implies strong evidence against the null model. Table 3 summarises the test results and shows that the tests favour VAR(1), VAR(2), and VAR(3) over VAR(4), indicating that lag order of <4 should be used. Similarly, the tests favour VAR(1) and VAR(2) over VAR(3), indicating that lag order of <3 should be used. Lastly, the test of VAR(1) versus VAR(2) rejects VAR(1) at the p-value of 0.0309. Therefore, we select VAR(2) to model $$\left[ {\Delta \kappa _{t}^{{(1)}} ,\Delta \kappa _{t}^{{(2)}} } \right]\prime$$ .

The parameter estimates and their standard errors (in parentheses) are displayed below:

This two-population Lee–Carter model can capture the correlation between the mortality rates of the two populations. However, it does not impose any non-divergence conditions. Researchers (Li & Lee, Reference Li and Lee2005; Cairns et al., Reference Cairns, Blake, Dowd, Coughlan and Khalaf-Allah2011; Dowd et al., Reference Dowd, Cairns, Blake, Coughlan and Khalaf-Allah2011; Li et al., 2015 Reference Li, Zhou and Hardyb ) have argued that mortality rates of different populations should not diverge indefinitely due to biological reasonableness. To incorporate the non-divergence assumption and study how it affects deep-deferred annuity valuation, we consider the following constraints:

1. 1. $$\beta _{x}^{{(1)}} \, {\equals}\, \beta _{x}^{{(2)}} $$ ;

2. 2. $$\Delta \kappa _{x}^{{(1)}} $$ and $$\Delta \kappa _{x}^{{(2)}} $$ share the same unconditional mean.

Under these two constraints, the mortality improvement rates of the two populations have the same unconditional mean. The female and male mortality rates will eventually run in parallel. We follow Zhou et al. (Reference Zhou, Wang, Kaufhold, Li and Tan2014) and estimate this model by maximum likelihood estimation.

Figure 4 plots the estimates of $$\alpha _{x}^{{(i)}} $$ , $$\beta _{x}^{{(i)}} $$ , and $$\kappa _{t}^{{(i)}} $$ when the non-divergence constraints are applied. We observe that $$\alpha _{x}^{{(i)}} $$ and the common $$\beta _{x}^{{(i)}} $$ follow similar patterns as those in Figure 2. Therefore, we continue to extrapolate $$\alpha _{x}^{{(i)}} $$ by linear function and $$\beta _{x}^{{(i)}} $$ by Gaussian or quadratic function. The new parameter estimates of the VAR(2) model for $$[ {\Delta \kappa _{t}^{{(1)}} ,\Delta \kappa _{t}^{{(2)}} } ]\prime$$ are $$A_{0} \, {\equals}\, \left[ {\matrix{ {\, {\minus}0.9872} \cr {\, {\minus}1.0340} \cr } } \right],$$ $$A_{1} \, {\equals}\, \left[ {\matrix{ {1.1344} &#x0026; {\, {\minus}1.2043} \cr {0.9852} &#x0026; {\, {\minus}1.0241} \cr } } \right],$$ $$A_{2} \, {\equals}\, \left[ {\matrix{ {\, {\minus}0.7530} &#x0026; {0.6210} \cr {\, {\minus}1.4681} &#x0026; {1.2480} \cr } } \right],$$ and $$\Sigma \, {\equals}\, \left[ {\matrix{ {0.4674} &#x0026; {0.5397} \cr {0.5397} &#x0026; {0.6963} \cr } } \right].$$

Figure 4 Estimates of $$\alpha _{x}^{{(i)}} $$ , $$\beta _{x}^{{(i)}} $$ , and $$\kappa _{t}^{{(i)}} $$ of the two-population Lee–Carter model with non-divergence constraints.

Figure 5 shows the mean mortality forecasts of an 80-year-old Japanese male and an 80-year-old Japanese female over 40 years using the three dependence structures: independent (left), correlated but allowing divergence (centre), and correlated with non-divergence constraints (right). The upper panel uses Gaussian extrapolation for β _x and the lower panel uses quadratic function with a floor of 0. Perks’ model is used in the base mortality curve for illustration. When using Gaussian extrapolation for β _x , there is no obvious difference between the left and centre plots. Incorporating correlation without imposing non-divergence constraints leads to marginal changes in the mean forecasts of both male and female mortality. However, comparing the centre and right plots of the upper panel, we find that the non-divergence constraints shift male mortality significantly downward and female mortality slightly upward. This may have a major impact on the annuity values.

Figure 5 Force of mortality forecasts of an 80-year-old over the next 40 years using the three dependence structures, assuming Gaussian extrapolation (a) and quadratic extrapolation (b) for β _x , and the Perks’ mortality curve.

To understand the reasons behind this difference, we calculate the unconditional mean of $$\Delta \kappa _{t}^{{(1)}} $$ and $$\Delta \kappa _{t}^{{(2)}} $$ under all three structures, as presented in Table 4. When no constraints are imposed, $$E( {\Delta \kappa _{t}^{{(1)}} } )$$ is significantly higher than $$E( {\Delta \kappa _{t}^{{(2)}} } )$$ , indicating slower mortality improvements for males. After imposing the non-divergence constraints, $$E( {\Delta \kappa _{t}^{{(1)}} } )$$ decreases and $$E( {\Delta \kappa _{t}^{{(2)}} } )$$ increases, such that they agree on the same value. With faster mortality improvement, male mortality forecasts become lower. The opposite occurs for female mortality forecasts.

When using quadratic extrapolation for β _x , the three plots in the lower panel of Figure 5 show no difference for ages above 105. In Figure 3, we have seen that β _x decreases quickly and reaches the floor of 0 at x=105. Male and female mortality will not improve further after age 105. Therefore, the choice of correlation structure has no influence on the mortality forecasts for ages above 105.

4.1. Measuring model risk

According to Draper (Reference Draper1995), a model formalises assumptions about how unknown quantities are related to known quantities. A model typically has two parts: structural assumption, such as the form of the link function and the choice of error distribution in a generalised linear model; and parameters, whose meaning is specific to a given choice of structural assumptions. The uncertainty of a model arises from two sources:

1. 1. Uncertainty in the values of parameters given a specific model structure, i.e. parameter risk; and

2. 2. uncertainty in the model structure itself, i.e. model risk.

In this section, we aim to quantify the exposure of annuity value to the second source of uncertainty. To simplify the problem and maintain our focus on model risk, we assume that annuity value is the actuarial present value of all future benefit payments. Consider an n-year deferred annuity issued to a life aged x at the beginning of 2013. This annuity pays $10,000 at the beginning of each year contingent on the survival of the annuitant. Its value at issue can be expressed as

where ω is the limiting age, i the annual effective interest rate, and _k p _x the probability that this life survives k years. When n=0, this annuity becomes an immediate annuity. Assuming constant force of mortality between integral ages, _k p _x can be approximated by

Danielsson et al. (Reference Danielsson, James, Valenzuela and Zer2015) proposed to use risk ratio to quantify model risk. The risk ratio measures the level of disagreement amongst candidate models by calculating the ratio of the maximum to the minimum risk forecasts. The difference between risk ratio and one indicates the extent to which the candidate models disagree. When risk forecasts come from similar models, then risk ratio should be close to one. In this paper, we use 95% value at risk (VaR) of annuity values as the risk forecasts. If we have N candidate models, then the risk ratio can be expressed as

where i denotes the ith model and i=1, 2, … , N. The mortality model we use for annuity valuation is comprised of two parts: mortality tail curve and mortality improvement rate model. We examine the model risk in these two parts separately. Since we do not know the exact distribution of annuity value, we use the 95th percentile of simulated annuity values as the 95% VaR.

4.2. Model risk in mortality tail curves

To start with, we assume there is no mortality improvement and examine how the annuity value changes with the choice of mortality tail curves. With zero mortality improvement rate, we have $_{k} p_{x} \, {\equals}\, \prod\nolimits_{s\, {\equals}\, 0}^{k\, {\minus}\, 1} e^{{\, {\minus}\, m_{{x{\plus}s,2012}} }} $ . Once the mortality curve in 2012 is determined, there is no uncertainty in _k p _x . Therefore, annuity value is a constant and the 95% VaR of annuity value is equal to this constant.

We use immediate annuities and deep-deferred annuities with various deferral periods for illustration. Assume that these annuities are purchased at age 75 and annual effective interest rate is 2%. Table 5 lists risk ratios and average annuity values using various mortality tail curves. Although we observed a large discrepancy of mortality rates at high ages in Figure 1, Table 5 shows that there is insignificant difference in immediate annuity values from using different mortality curves because the three mortality curves only start to diverge at approximately age 95. Immediate annuity values are heavily weighted by benefits paid in earlier years due to less discounting and higher survival probability. The actuarial present value of benefits paid after age 95 is relatively small. Therefore, the choice of mortality tail curve does not have a significant impact on the values of immediate annuities.

As expected, deep-deferred annuities are much less costly than immediate annuities. The value of a 15-year deferred annuity issued to a 75-year-old male is ~9.7% of the corresponding immediate annuity value because a 75-year-old male has only an ~54% chance of surviving to age 90, and the annuity payments at high ages are deeply discounted.

Although deep-deferred annuities have a low-cost advantage, they have a higher model risk than immediate annuities. Figure 6 plots the risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ for n=0, 1, 2, … , 20 under various interest rate assumptions. While the focus of this discussion is on males, the female plot can be analysed in a similar fashion.

Figure 6 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst the three mortality tail curves and zero mortality improvement rate.

The pattern of risk ratio is determined largely by the relative impacts of mortality rates in two age ranges: 82–94 and above 95. As the deferral period increases, the impact of age range 82–94 decreases while that of ages above 95 increases due to increasing weights of payments in later life. Recall that in Figure 1, Gompertz law has the lowest force of mortality for males in the age range of 82–92 and the highest force of mortality for males who are age 95 and above among the three mortality curves. In contrast, Perks’ model has the highest or close to the highest force of mortality for males in the age range of 82–92 and the lowest force of mortality for males who are age 95 and above.

When n≤18, the impact of age range 82–94 outweighs that of ages over 95 and causes Gompertz law to yield the highest annuity value. Because the gap between Gompertz and Perks’ mortality curves increases from age 82 to age 88, the risk ratio increases for n<15. When 15≤n≤18, the risk ratio decreases. The reasons are twofold: (1) the gap between Gompertz mortality curve and Perks’ mortality curve shrinks from age 89 to age 92; (2) the impact of mortality at age 95 and above weighs more in annuity values as n increases. With n continuing to increase, the impact of mortality at age 95 and above eventually outweighs that of mortality in the age range of 82–94. This results in Perks’ model having the highest annuity value for n>18, and hence a kink point at n=18 in the left plot of Figure 6. As the gap between Perks’ and Gompertz mortality curves enlarges quickly for ages above 95, the risk ratio also increases steeply for n>18.

The impact of interest rate also depends on the length of the deferral period. We observe that higher interest rate leads to higher risk ratio for 10≤n≤18, but lower risk ratio for n>18. When interest rate increases, the present value of benefit payment becomes lower. However, the present value of payments made in the age range of 82–94 is less affected than the value of payments made at age 95 and above due to discounting. Therefore, the impact of mortality in the age range of 82–94 on annuity value is amplified for deferred annuities with 10≤n≤18, thus resulting in higher risk ratio. When n>18, annuity value is mostly determined by the mortality at age 95 and above. With a lower interest rate, payments made at higher ages when a larger discrepancy of mortality rate is observed between Gompertz and Perks’ mortality curve gain more weight in annuity value. The risk ratio increases as a result.

Next we incorporate mortality improvement rates and further study the model risk from mortality tail curves. We first use the single-population Lee–Carter model to obtain 1,000 trajectories of future mortality improvement rates. We then multiply the improvement rates by the base year mortality curve to calculate simulated trajectories of mortality rates. Again, we consider the three types of mortality tail curves described in section 2. Simulated annuity values are determined for each combination of mortality improvement trajectory and base mortality curve. Table 6 lists the average simulated annuity value using various mortality tail curves and the risk ratios assuming Gaussian extrapolation for β _x and i=2%.

After incorporating mortality improvement trend, the probability of a 75 years old surviving to advanced ages increases. As a result, annuity prices increase significantly compared to those in Table 5, especially for annuities with a longer deferral period. The risk ratio for 20-year deferred annuity becomes higher. The value of a 20-year deferred annuity is largely driven by mortalities at ages above 95. With higher surviving probabilities, the actuarial present value of payments made at high ages gains more weight in annuity value. The large discrepancy between different mortality curves at age 95 and above exerts more impacts, thereby resulting in higher risk ratios. In contrast, the values of 10-year and 15-year deferred annuities are mostly determined by mortalities at ages below 95, where the discrepancy between mortality curves is small. Hence, their risk ratios are also small.

In Figure 7, we plot the relation between risk ratios amongst mortality tail curves and length of deferral period assuming various interest rates and extrapolation methods for β _x . Due to the reasons explained above, risk ratios increase sharply when the deferral period is >17. We also observe that risk ratios using Gaussian extrapolation for β _x are slightly higher than those using quadratic extrapolation. Gaussian extrapolation admits more mortality improvement for ages above 95 and thus increases the weights of payments at these ages in the annuity value. Since the discrepancy between mortality curves for these ages is particularly large, higher weights of payments at these ages lead to higher risk ratios.

Figure 7 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst the three mortality tail curves, assuming single-population Lee–Carter model for mortality improvement, Gaussian extrapolation for β _x (a) and quadratic extrapolation for β _x (b).

Overall, the model risk from mortality tail curve appears to be low for annuities issued to a 75-year-old male with a deferral period of less than 17 years. When the deferral period is longer than 17 years, the model risk increases sharply. In practice, deferred annuities often have payments starting no later than age 85. Therefore, model risk from mortality tail curve should not be a major concern.

Figure 8 examines the total model risk from the two sources: mortality tail curve and extrapolation method for β _x . The risk ratios are based on six combinations of mortality tail curve and extrapolation method for β _x , since we have three mortality tail curves and two extrapolation methods. Since extrapolation method for β _x is part of the mortality improvement model, the effects of adding mortality improvement are twofold: (1) amplifying the model risk from mortality tail curves as shown in Table 6; and (2) introducing model risk from extrapolation method for β _x . A comparison of Figures 7 and 8 indicates that male annuity values are subject to more model risk from mortality tail curve while female annuity values are subject to more model risk from extrapolation method for β _x .

Figure 8 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst six combinations of mortality tail curve and extrapolation method for β _x , assuming single-population Lee–Carter model for mortality improvement.

4.3. Model risk in two-population Lee–Carter model

In this section, we study model risk in the two-population Lee–Carter model. We consider three different assumptions about the dependence structure between male and female mortality improvement:

• ∙ independent;

• ∙ correlated but allowed to diverge over the long run;

• ∙ correlated and non-divergent over the long run.

Table 7 lists risk ratios and average annuity values using various dependence structures assuming Gaussian extrapolation for β _x , Perks’ mortality curve, and 2% annual effective interest rate. The difference in average annuity values under these three structures can be explained by the mortality forecasts in Figure 5. When imposing the non-divergence constraints with Gaussian extrapolation for β _x , male mortality shifts down significantly and female mortality shifts up slightly. This leads to higher male annuity values and lower female annuity values. A comparison of Tables 5–7 indicates that model risk caused by two-population Lee–Carter model is much more significant than that caused by mortality for male annuities. For female annuities, the model risk from all three sources is low.

Figure 9 plots the risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ , for n=0, 1, … , 20, among the three dependence structures assuming Gaussian extrapolation for β _x and Perks’ mortality curve. Annuities with a longer deferral period have cash flows concentrated at higher ages. Since the difference of mortality forecasts from using the three structures increases with the length of the forecasting period, a longer deferral period results in higher risk ratios. Gaussian extrapolation for β _x allows for more mortality improvement and heavier cash flows at advanced ages than quadratic extrapolation. It amplifies the impact of dependence structures on annuity values, thereby resulting in higher risk ratios than quadratic extrapolation.

Figure 9 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst various correlation structures with Perks’ mortality curve, Gaussian extrapolation for β _x (a) and quadratic extrapolation for β _x (b).

Figure 10 studies how risk ratios change when we gradually add the three sources of model risk. Annual effective interest rate is assumed to be 2%. The left panel shows that the model risk of male annuity values is high, and it is primarily attributable to dependence structure. The model risk from extrapolation method for β _x and mortality curve is almost negligible for annuities with deferral period of less than 10 years; however, this risk becomes more significant when the deferral period increases. The right panel shows that the model risk of female annuity values mainly comes from dependence structure and extrapolation method for β _x , and this risk is much lower than the model risk of male annuity values. Overall, the model risk from the two-population mortality improvement model is much more significant than that from mortality tail curve.

Figure 10 Comparison of risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ in three cases: Case 1: risk ratios amongst dependence structures, assuming Gaussian extrapolation for β _x and Perks’ mortality curve; Case 2: risk ratios amongst six combinations of dependence structure and extrapolation method for β _x , assuming Perks’ mortality curve; and Case 3: risk ratios amongst 18 combinations of dependence structure, extrapolation method for β _x , and mortality tail curve. i=2%.

5.1. Incorporating parameter risk

According to Draper (Reference Draper1995), a major source of uncertainty is parameter risk. The values of parameters given a specific model structure have uncertainty because only a finite set of data are used to estimate these parameters and the estimates cannot be exact. Due to the sparseness of high age mortality data, parameter risk may be more severe for deep-deferred annuities.

To quantify parameter uncertainty in mortality tail curve and mortality improvement model, we adopt the parametric bootstrap method. To apply bootstrap for mortality tail curve, we follow Brouhns et al. (Reference Brouhns, Denuit and Van Keilegom2005) and bootstrap death counts by applying a Poisson noise to the observed numbers of deaths as follows:

1. 1. Generate one bootstrap sample of death counts in year 2012, assuming that the bootstrapped death counts are from Poisson distributions with mean equal to the observed number of deaths.

2. 2. Estimate the three mortality tail curves by using the bootstrap sample.

3. 3. Construct base mortality curves using each estimated mortality tail curve.

Since we use a two-stage method to fit the Lee–Carter model, parameter risk comes from both stages and bootstrapping should be performed for both stages as well. In the first stage, we estimate the Lee–Carter model and obtain estimates of $$\alpha _{x}^{{(i)}} $$ , $$\beta _{x}^{{(i)}} $$ , and $$\kappa _{t}^{{(i)}} $$ . To capture the parameter risk in this stage, we again follow Brouhns et al. (Reference Brouhns, Denuit and Van Keilegom2005) and bootstrap death counts by applying a Poisson noise to the observed numbers of deaths. In the second stage, we estimate the extrapolation function using the estimates of $$\alpha _{x}^{{(i)}} $$ and $$\beta _{x}^{{(i)}} $$ , and estimate the time series model using the estimates of $$\kappa _{t}^{{(i)}} $$ from the first stage. To take into account the parameter risk in this stage, we apply residual bootstrap.

The detailed procedures to apply bootstrap in the Lee–Carter model are as follows:

1. 1. Generate one bootstrap sample of death counts for the sample period of 1965–2012, assuming that the bootstrapped death counts are from Poisson distributions with mean equal to the observed number of deaths.

2. 2. Estimate the mortality improvement rate model using the bootstrap sample and obtain estimates of $$\alpha _{x}^{{(i)}} $$ , $$\beta _{x}^{{(i)}} $$ , and $$\kappa _{t}^{{(i)}} $$ .

3. 3. Estimate the extrapolation functions for $$\alpha _{x}^{{(i)}} $$ and $$\beta _{x}^{{(i)}} $$ .

4. 4. Generate one bootstrap sample of $$\alpha _{x}^{{(i)}} $$ and $$\beta _{x}^{{(i)}} $$ using residual bootstrapFootnote ⁵ and re-estimate the extrapolation functions by using this sample.

5. 5. Extrapolate $$\alpha _{x}^{{(i)}} $$ and $$\beta _{x}^{{(i)}} $$ to desired age using the re-estimated function.

6. 6. Estimate the time series model for $$\kappa _{t}^{{(i)}} $$ .

7. 7. Generate one bootstrap sample of $$\kappa _{t}^{{(i)}} $$ using residual bootstrapFootnote ⁶ and re-estimate the time series model by using this sample.

8. 8. Simulate one path of $$\kappa _{t}^{{(i)}} $$ , for t=2013, 2014, … , using the re-estimated time series model.

9. 9. Calculate the path of future mortality improvement rates based on extrapolated $$\alpha _{x}^{{(i)}} $$ and $$\beta _{x}^{{(i)}} $$ , and simulated path of $$\kappa _{t}^{{(i)}} $$ .

After obtaining the base mortality curve and future mortality improvement rates, we determine the path of future mortality rates and calculate a simulated annuity value. These procedures are repeated 1,000 times to obtain 1,000 simulated annuity values. The new risk ratio is then calculated as follows:

where i=1, 2, … , N and N is the number of models investigated. The extent of parameter risk is indicated by the difference between the new risk ratio and the risk ratio without considering parameter uncertainty as defined in equation (4.1).

5.2. Parameter risk in single-population mortality model

A full single-population mortality model includes a base mortality curve and a single-population mortality improvement model. To examine its parameter risk, we follow the bootstrap procedures described in section 5.1 and obtain simulated annuity values for each combination of mortality tail curve and extrapolation method for β _x . We then determine $\:{\rm VaR}_{{95\,\&#x0025;\,}} \left( {_{{n\,\mid\,}} \ddot{a}_{x} } \right)$ for all six models considered and use the highest one as the numerator in equation (5.1) to calculate the new risk ratio.

Figure 11 plots the risk ratios of $_{{n\,\mid\,}} \ddot{a}_{{75}} $ , for n=0, 1, 2, … , 20, amongst six combinations of mortality tail curve and extrapolation method for β _x . We consider three cases: no bootstrapping, partial bootstrapping, and full bootstrapping. The dark lower part of each bar displays the risk ratio without considering parameter uncertainty. It represents the model risk from mortality tail curve and extrapolation method for β _x used in the single-population Lee–Carter model.

Figure 11 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst six combinations of mortality tail curve and extrapolation method for β _x , with parameter uncertainty incorporated. i=2%.

In partial bootstrapping, we apply bootstrap only to the mortality tail curve, and not to the mortality improvement model. Partial bootstrapping increases risk ratios by adding parameter uncertainty from mortality tail curve. The grey part of each bar depicts the increased risk, i.e. parameter risk from mortality tail curve.

In full bootstrapping, we apply bootstrap to both the mortality tail curve and the mortality improvement model. Full bootstrapping further increases risk ratios by adding parameter uncertainty from the mortality improvement rate model. The upper white part of each bar represents the increased risk, i.e. parameter risk from the mortality improvement model.

For male annuity values, both model risk and parameter risk play important roles. However, parameter risk from mortality tail curve is much smaller than that from mortality improvement rate model. For female annuity values, model risk is of more concern, while parameter risk appears to be low.

5.3. Parameter risk in two-population mortality models

In section 4.3, we examined model risk from the dependence structure in the two-population Lee–Carter model. We found that the model risk is significant, especially for males, and this risk increases steadily with the length of deferral period. In this section, we study the parameter risk in the two-population mortality model and compare it with the model risk.

We consider 18 combinations of mortality tail curve, extrapolation method for β _x , and dependence structure. For each combination, we use the bootstrapping method described in section 5.1 to simulate 1,000 annuity values and calculate ${\rm VaR}\,\left( {_{{n\,\mid\,}} \ddot{a}_{{75}} } \right)$ . The highest ${\rm VaR}\,\left( {_{{n\,\mid\,}} \ddot{a}_{{75}} } \right)$ is used to calculate the new risk ratio.

Figure 12 demonstrates the risk ratios of $_{{n\,\mid\,}} \ddot{a}_{{75}} $ amongst the 18 combinations of mortality tail curve, extrapolation method for β _x , and dependence structure. Three levels of bootstrapping are considered: no bootstrapping, partial bootstrapping, and full bootstrapping. In partial bootstrapping, we again apply bootstrap only to the mortality tail curve, and not to the mortality improvement model. The dark lower part of each bar depicts the risk ratios without parameter uncertainty. It measures the variation of ${\rm VaR}\,\left( {_{{n\,\mid\,}} \ddot{a}_{{75}} } \right)$ stemming from model risk. The grey part is almost invisible indicating that parameter risk from mortality tail curve is negligible. The upper white part represents the increase of risk ratio due to parameter uncertainty in the two-population mortality improvement model, which turns out to be the most significant source of risk in both male and female annuity values. The overall risk in male annuity values is much higher than that in female annuity values.

Figure 12 Risk ratios of $$_{{n\,\mid\,}} \ddot{a}_{{75}} $$ amongst 18 combinations of mortality tail curve, extrapolation method for β _x , and dependence structure, with parameter uncertainty incorporated. i=2%.