Adding Shocks to a Prospective Mortality Model

1. Introduction

The construction of life expectancy projections has been the subject of much work since the seminal article by Lee and Carter (Lee and Carter [1992]).

With a view to extrapolating trends observed in the past into the future, most of the approaches proposed are based on a “mortality surface”, measuring mortality forces by age and year at the time, which must then be extrapolated over time.

The models inspired by Lee and Carter start by reducing dimension by performing a PCA and then extrapolating one or two time series associated with the projection on the principal axes.

Bongaarts (Bongaarts [2004]) has proposed a different approach, based on parametric adjustments by year of time and extrapolation of the estimated coefficients each year.

In Bongaarts [2004], however, the author uses a rather frustrating parametric representation (Thatcher’s model) which does not allow all ages to be included. Moreover, he limits his extrapolation to 2 out of 3 parameters, treating them independently, which is a questionable approximation.

Models of this type project regular $t\to \mu \left(x,t\right)$ series. However, when we look at annual variations in the mortality rate in France, for example, we see a fairly high degree of volatility:

Figure 1. Annual variation in mortality rate.

Figure 1. Annual variation in mortality rate.

The classic models described above cannot account for these short-term variations. Proposed approaches have been put forward, for example, in Guette [2010] or Currie et al. [2003], but with a slightly different objective, as these works propose to model catastrophes such as wars or epidemics. More recently, an approach using regime-switching models was proposed in Robben and Antonio [2023].

Our aim here is not to model catastrophes, but to incorporate the above volatility into the modeling to provide a more accurate assessment of residual life expectancies when an unbiased estimate of mortality rates is available.

A specific approach is therefore proposed here, with the aim of accounting for this short-term volatility and measuring its impact on the anticipation of prospective residual life expectancies.

2. Proposed Stochastic Mortality Model

We draw inspiration here from the fragility models proposed by Vaupel and his co-authors (Vaupel et al. [1979]), by assigning a regular base hazard function of a shock that depends only on time, under a proportional hazards assumption.

The proposed specification is described below, followed by a method for estimating the parameters within the framework of conditional maximum likelihood.

a. Specification

Consider the following specification of the hazard function for the year of time t:

$$\mu \left(x,t\right)=Z_t\times {\mu }_0\left(x,t\right)$$

with the semi-parametric form of the basic hazard function $\mathbf{ln}{\mu }_0\left(x,t\right)={\alpha }_x+{\beta }_x{k}_t$. The shocks are assumed to be mean-centered, i.e., $E\left(Z_t\right)=1$ and the classical identifiability conditions for the basic hazard function are imposed (cf. Brouhns et al. [2002]):

$${\displaystyle \sum _{x=x_m}^{x_M}{\beta }_x=1}\mathrm{and}{\displaystyle \sum _{t=t_m}^{t_M}{k}_t=0}$$

This involves estimating the parameters of $Z_t$ and the matrix$\left(\alpha ,\beta ,k\right)$, then extrapolating the time series$t\to k\left(t\right)$.

b. Log-Likelihood Determination

For maximum likelihood estimation, we know that everything happens as if the number of deaths observed followed a Poisson distribution,

$$D_{x,t}\sim P\left(E_{x,t}\times \mu \left(x,t\right)\right),$$

which leads to the following expression for the conditional likelihood for an observation, noting $\lambda =E_{x,t}\times {\mu }_0\left(x,t\right)$:

$$P\left(D=d\left|Z\right.\right)=e^{-\lambda Z}\frac{{\lambda }^d}{d!}{Z}^d,$$

Likelihood is easy to deduce:

$$P\left(D=d\right)=E_Z\left[P\left(D=d\left|Z\right.\right)\right]={\displaystyle \int e^{-\lambda z}{z}^d\frac{{\lambda }^d}{d!}d{F}_Z\left(z\right)}$$

We then choose a Gamma distribution of parameters a and b for the distribution of Z, i.e., $f_Z\left(z\right)=z^{a-1}\frac{b^a{e}^{-b\times z}}{\Gamma \left(a\right)}$, which leads to $P\left(D=d\right)=\frac{{\lambda }^d}{d!}\frac{b^a}{\Gamma \left(a\right)}{\displaystyle \underset{0}{\overset{+\infty }{\int }}{e}^{-\left(\lambda +b\right)z}{z}^{d+a-1}dz}$. Using the change of variable $u=\left(\lambda +b\right)z$, we obtain the following expression for the likelihood of an observation

$$P\left(D=d\right)=\frac{{\lambda }^d}{d!}\frac{b^a}{\Gamma \left(a\right)}\frac{1}{{\left(\lambda +b\right)}^{d+a}}{\displaystyle \underset{0}{\overset{+\infty }{\int }}{e}^{-u}{u}^{d+a-1}dz}=\frac{{\lambda }^d}{d!}\frac{b^a}{\Gamma \left(a\right)}\frac{1}{{\left(\lambda +b\right)}^{d+a}}\Gamma \left(d+a\right)$$

which in turn gives log-likelihood

$$\mathbf{ln}P\left(D=d\right)=f\left(a,b\right)+d\mathbf{ln}\left(\lambda \right)-\left(d+a\right)\mathbf{ln}\left(\lambda +b\right)$$

with $f\left(a,b\right)=\mathbf{ln}\left(b^a\frac{\Gamma \left(d+a\right)}{\Gamma \left(d+1\right)\Gamma \left(a\right)}\right)$ .

As a function of the parameters $\left(\alpha ,\beta ,k\right)$ and conditional on $\left(a,b\right)$, the log-likelihood for an observation is of the form $l\left(\alpha ,\beta ,k\right)=\mathbf{ln}P\left(D=d\right)$ with $\lambda =E_{x,t}\times {\mu }_0\left(x,t\right)=E_{x,t}\times e^{{\alpha }_x+{\beta }_x{k}_t}$. Conditional on $\left(a,b\right)$, the log-likelihood has the following form (with one additive constant):

$$\mathbf{ln}L={\displaystyle \sum _{x=x_m}^{x_M}{\displaystyle \sum _{t=t_m}^{t_M}\mathbf{ln}P\left(D_{x,t}=d_{x,t}\left|a,b\right.\right)}}={\displaystyle \sum _{x=x_m}^{x_M}{\displaystyle \sum _{t=t_m}^{t_M}\left[d_{x,t}\mathbf{ln}\left({\lambda }_{x,t}\right)-\left(d_{x,t}+a\right)\mathbf{ln}\left({\lambda }_{x,t}+b\right)\right]}}$$

It can then be maximized by $\left(\alpha ,\beta ,k\right)$ under the constraint that ${\displaystyle \sum _{x=x_m}^{x_M}{\beta }_x=1}$ and ${\displaystyle \sum _{t=t_m}^{t_M}{k}_t=0}$.

c. Parameter Estimation

Parameter estimation can be carried out in two stages: in the first stage, the fragility parameter is estimated, then, in the second stage, the above log-likelihood is maximized at $\left(\alpha ,\beta ,k\right)$.

The condition $E\left(Z_t\right)=1$ implies $a=b$. We also have $V\left(Z_t\right)=\frac{a}{b^2}=\frac{1}{a}$, so the disturbance control parameter $Z_t$ is the inverse of the variance $a={\sigma }_Z^{-2}$. A direct estimate of this parameter can be made as follows, observing that the mean annual output intensities are of the form

$$\overline{\mu }\left(t\right)=Z_t\times {\overline{\mu }}_0\left(t\right)\mathrm{with}{\overline{\mu }}_0\left(t\right)=\frac{{\displaystyle \sum _{x=x_m}^{x_M}{E}_{x,t}{\mu }_0\left(x,t\right)}}{{\displaystyle \sum _{x=x_m}^{x_M}{E}_{x,t}}}$$

from which we derive $E\left(\overline{\mu }\left(t\right)\right)={\overline{\mu }}_0\left(t\right)$, $V\left(\overline{\mu }\left(t\right)\right)=V\left(Z_t\right){\overline{\mu }}_0^2\left(t\right)$ then $V\left(Z_t\right)=\frac{V\left(\overline{\mu }\left(t\right)\right)}{E{\left(\overline{\mu }\left(t\right)\right)}^2}$. The variance of $Z_t$ is therefore equal to the square of the coefficient of variation of $\overline{\mu }\left(t\right)$: $\sigma \left(Z_t\right)=cv\left(\overline{\mu }\left(t\right)\right)$. It is then straightforward to propose as estimator

$${\widehat{\sigma }}_Z^2=\frac{\frac{1}{t_M-t_m+1}{\displaystyle \sum _{t=t_m}^{t_M}{\left(\widehat{\mu }\left(t\right)-\frac{1}{t_M-t_m+1}{\displaystyle \sum _{t=t_m}^{t_M}\widehat{\mu }\left(t\right)}\right)}^2}}{{\left(\frac{1}{t_M-t_m+1}{\displaystyle \sum _{t=t_m}^{t_M}\widehat{\mu }\left(t\right)}\right)}^2}$$

With $\widehat{\mu }\left(t\right)=\frac{{\displaystyle \sum _{x=x_m}^{x_M}{E}_{x,t}\widehat{\mu }\left(x,t\right)}}{{\displaystyle \sum _{x=x_m}^{x_M}{E}_{x,t}}}$ and $\widehat{\mu }\left(x,t\right)=\frac{D_{x,t}}{E_{x,t}}$ as estimators, the Hoem estimator of the hazard function.

Once the fragility parameter has been estimated, the aim is to maximize the previously expressed log-likelihood. The partial derivatives of the log-likelihood for an observation are as follows:

$$\frac{\partial }{\partial p}\mathbf{ln}P\left(D=d\right)=d\frac{1}{\lambda }\frac{\partial \lambda }{\partial p}-\left(d+a\right)\frac{1}{\lambda +b}\frac{\partial \lambda }{\partial p}=\left(\frac{d}{\lambda }-\frac{d+a}{\lambda +b}\right)\frac{\partial \lambda }{\partial p}$$

with p one of the parameters $\left(\alpha ,\beta ,k\right)$. We also have $\frac{\partial \lambda }{\partial \alpha }=\lambda$, $\frac{\partial \lambda }{\partial \beta }=k\lambda$ and$\frac{\partial \lambda }{\partial k}=\beta \lambda$.

Finally, it remains to estimate $\left(\alpha ,\beta ,k\right)$, which is a solution of the first-order conditions:

$$\frac{\partial }{\partial {\alpha }_x}\mathbf{ln}L={\displaystyle \sum _{t=t_m}^{t_M}\left(\frac{d}{{\lambda }_{x,t}}-\frac{d+a}{{\lambda }_{x,t}+b}\right){\lambda }_{x,t}}=0$$

$$\frac{\partial }{\partial {\beta }_x}\mathbf{ln}L={\displaystyle \sum _{t=t_m}^{t_M}\left(\frac{d}{{\lambda }_{x,t}}-\frac{d+a}{{\lambda }_{x,t}+b}\right)k_t{\lambda }_{x,t}}=0$$

$$\frac{\partial }{\partial k_t}\mathbf{ln}L={\displaystyle \sum _{x=x_m}^{x_M}\left(\frac{d}{{\lambda }_{x,t}}-\frac{d+a}{{\lambda }_{x,t}+b}\right){\beta }_x{\lambda }_{x,t}}=0$$

This system is non-linear.

d. Calculating Prospective Residual Life Expectancies

In the proposed model, the calculation of prospective life expectancy

$$e\left(x,t\right)=E\left(e\left(x,t\left|Z\right.\right)\right)={\displaystyle \sum _{i\ge 0}{\displaystyle \prod _{j=0}^{i}E\left[\exp \left(-{\mu }_{x+j,t+j}\right)\right]}}$$

(1)

takes the form

$$e\left(x,t\right)={\displaystyle \sum _{i\ge 0}{\displaystyle \prod _{j=0}^i{\left(\frac{b}{b+{\mu }_0\left(x+j,t+j\right)}\right)}^a}}$$

(2)

because the Laplace transform of a Gamma distribution is $E\left(e^{-x{Z}_t}\right)={\left(\frac{b}{b+x}\right)}^a$. As $E\left(e^{-\mu \left(x,t\right)}\right)=E\left(e^{-Z_t\times {\mu }_0\left(x,t\right)}\right)$, we deduce that:

$$E\left(e^{-\mu \left(x,t\right)}\right)=E\left(e^{-Z_t\times {\mu }_0\left(x,t\right)}\right)={\left(\frac{b}{b+{\mu }_0\left(x,t\right)}\right)}^a$$

(3)

Note that when $b=a\to +\infty$, $E\left(e^{-\mu \left(x,t\right)}\right)\to e^{-{\mu }_0\left(x,t\right)}$ and then we find the classic formula $e\left(x,t\right)={\displaystyle \sum _{i\ge 0}{\displaystyle \prod _{j=0}^i\mathbf{exp}\left(-{\mu }_0\left(x+j,t+j\right)\right)}}$.

3. Numerical Application

We use data for metropolitan France for the period 2000-2020, for ages 0 to 105 inclusive, from the Agalva and Blanpain study [2021]. All calibration was performed in R.

Prospective analyses are then carried out over the entire age range and for the years from 2021 to 2060.

a. Model Adjustment

All these steps are discussed in turn in the following subsections. Throughout the study, all the results obtained are compared with those given by a Lee-Carter model calibrated on the same data.

Estimation of Gamma Distribution Parameters

The estimation of the pair of parameters $\left(a,b\right)$ has been made with the raw data and we find ${\sigma }_Z^2=4,3\%$ , i.e., $a=b=550$.

Estimation of Model Parameters

The calibration of $\left(\alpha ,\beta ,k\right)$ under constraints was then carried out using the Rsolnp package (https://cran.r-project.org/web/packages/Rsolnp/index.html), and more specifically the solnp function (see Ghalanos and Theussl [2015]). In order to carry out this calibration, it was chosen to provide as initial parameters the results of a Lee-Carter model, calibrated using the lca function from the demography package (https://cran.r-project.org/web/packages/demography/index.html) (cf. Hyndman [2023]). All these coefficients are transcribed in the Appendix.

As shown in the following figure, this leads to coefficients $\left(\alpha ,\beta ,k\right)$ very close to those of the reference Lee-Carter model:

Figure 2. Comparison of different coefficients with those of a Lee-Carter model.

Figure 2. Comparison of different coefficients with those of a Lee-Carter model.

It can be seen that the time parameter shows a slower rate of decline from the 18^ème year of the observation period.

Extrapolation of Time Coefficients

Whether in the model studied here or in the Lee-Carter model used as a reference, the projection of coefficients ${\left(k_t\right)}_t$ for $t$ beyond the calibration range has been done by linear regression, by fitting the following equation to the calibrated parameters:

$$k_t=m\times t+p$$

In both cases, we find the coefficients shown in the following table:

Table 1. Results of the extension of the time coefficients k.

	m	p
Model studied	-2,19	4401,98
Lee-Carter reference model	-2,19	4402,33

Table 1. Results of the extension of the time coefficients k.

	m	p
Model studied	-2,19	4401,98
Lee-Carter reference model	-2,19	4402,33

The results are logically very similar in both cases.

b. Projected Mortality Forces

Here, we compare projected mortality forces with and without shocks integrated into the model, as a function of age and year.

First, we look at the evolution of the mortality force as a function of age, for a few fixed years:

Figure 3. Mortality forces by age.

Figure 3. Mortality forces by age.

The mortality forces of the two models are very close, except at the highest ages, where the model with shocks tends to predict higher mortality forces. A closer look at the age range [90; 105] reveals the following:

Figure 4. Mortality rates by age from 90 to 105 years.

Figure 4. Mortality rates by age from 90 to 105 years.

The maximum relative deviation over the entire prospective range is 3.6% at age 80 for the year 2060.

We now compare the mortality forces of the two models over the entire prospective analysis period, for a few selected ages.

Figure 5. Mortality rates by year for selected ages.

Figure 5. Mortality rates by year for selected ages.

Overall, the mortality model studied tends to predict higher mortality forces than the reference Lee-Carter model at older ages, and the gap increases over time.

Before calculating prospective residual life expectancies, it’s worth looking at the overall impact on mortality forces. To this end, we calculate the following ratio:

$$r\left(x,t\right)=\frac{\mu \left(x,t\right)}{{\mu }_{LC}\left(x,t\right)}$$

with ${\mu }_{LC}\left(x,t\right)$ the mortality force derived from the reference Lee-Carter model and $\mu \left(x,t\right)$ model. This ratio is shown in the following figure:

Figure 6. Mortality force ratio over the entire age range and prospective analysis period (2021-2060).

Figure 6. Mortality force ratio over the entire age range and prospective analysis period (2021-2060).

The average value of this ratio for all ages and years combined is 100%, which means that the model studied is equivalent to the Lee-Carter model.

The impact of the introduction of shocks on adjustment is therefore negligible when assessed on a very global basis. However, the difference increases over time, leading the two models to diverge in the medium term and at older ages.

If we restrict ourselves to ages over 65, we obtain a weighted average equal to 99.8%.

In addition, the average mortality rate of the population, calculated on the basis of exposure to risk in 2020, evolves as follows:

Figure 7. Average mortality rate per year, all ages combined, from 2021 to 2060.

Figure 7. Average mortality rate per year, all ages combined, from 2021 to 2060.

It can be seen that the model studied is on average more pessimistic than the Lee-Carter model on mortality improvement in future mortality.

c. Estimating Prospective Residual Life Expectancies

It is then possible to look at the consequences of the mortality model studied on prospective residual life expectancies, first by variable age for a few fixed years, then by variable year for a few fixed ages.

As shown in the following figure, it appears that the mortality model studied does not greatly change prospective residual life expectancies compared to the reference mortality model, except possibly at high ages:

Figure 8. Trends in prospective residual life expectancy by age.

Figure 8. Trends in prospective residual life expectancy by age.

An enlargement of these graphs at older ages is shown below, with the algebraic difference between the two models for each year:

Figure 9. Change in prospective residual life expectancy from age 90 to 105.

Figure 9. Change in prospective residual life expectancy from age 90 to 105.

The maximum absolute difference between the two models is found at age 96 for the year 2060, and is worth 0.18, or around 65 days.

In this analysis, we return to the observation made in the previous section: it is at older ages that the difference with the reference Lee-Carter model is most marked.

Figure 10. Prospective residual life expectancy by year.

Figure 10. Prospective residual life expectancy by year.

Below is an enlargement for age 96 with the algebraic difference between the two models:

Figure 11. Trend in prospective residual life expectancy by year at age 96.

Figure 11. Trend in prospective residual life expectancy by year at age 96.

d. Sensitivity to Fragility Parameter

The volatility of fragility is estimated at 4.3%; however, over a longer period, this parameter can be higher, for example (https://www.ressources-actuarielles.net/C1256F13006585B2/0/39B54166464089AFC12572B0003D88C2/$FILE/20230921_FP.pdf?OpenElement) from 1982 to 2022, it comes out at 5.5%.

With this level of volatility, we note that $P\left(Z_t\ge 1,09\right)\approx 5\%$; 9% being the excess mortality rate (https://actudactuaires.typepad.com/laboratoire/2021/03/mortalit%C3%A9-en-france-en-2020-suite.html) for the year 2020, we can deduce that the probability of observing excess mortality at this level is of the order of 5%. Furthermore, $Va{R}_{99,5\%}\left(Z_t\right)\approx 1,15$, which corresponds to the mortality shock for the “mortality” risk module of the delegated regulation.

On the basis of the central table${\mu }_0\left(x,t\right)$ adjusted above, the prospective residual life expectancies associated with a volatility coefficient of 5.5% are recalculated using

$$e\left(x,t\right)={\displaystyle \sum _{i\ge 0}{\displaystyle \prod _{j=0}^i{\left(\frac{{\sigma }^{-2}}{{\sigma }^{-2}+{\mu }_0\left(x+j,t+j\right)}\right)}^{{\sigma }^{-2}}}}$$

which leads to the following results:

Figure 12. Evolution of prospective residual life expectancies from age 65 to 105 for selected years, with a new volatility coefficient.

Figure 12. Evolution of prospective residual life expectancies from age 65 to 105 for selected years, with a new volatility coefficient.

There is no significant difference between the two models with this new volatility setting.

e. Consequences for the Capital Requirement of an Annuity Plan

The presence of the frailty factor therefore has no material impact on central tendency indicators (mortality forces, residual life expectancies, etc.).

However, the random nature of the mortality distribution in a given year has consequences for the assessment of the capital required to protect against adverse deviations in mortality. In the specific context of a life annuity plan, we are therefore led, following a logic analogous to that of the Solvency 2 standard, to consider the 99.5% quantile of the distribution of residual life expectancies induced by frailty. For each age from 60 to 100, we obtain the following results for the ratio between this quantile and the expectation:

Figure 13. Ratio between 99.5% quantile and expectation (SCR).

Figure 13. Ratio between 99.5% quantile and expectation (SCR).

Weighted by the age structure of the French population, the average ratio is around 101.3%.

For its part, the delegated regulation (EU Delegated Regulation n°2015/35: https://eur-lex.europa.eu/legal-content/FR/TXT/?uri=CELEX:32015R0035) imposes a 20% discount on death rates when calculating the SCR associated with longevity risk (cf. art. 138 of delegated regulation EU n°2015/35), which leads to a capital requirement equal to 10% of the expectation.

This means that the volatility observed in annual death rates explains around 12% of the longevity SCR.

4. Conclusions and Discussion

The use of a Gamma frailty model enables us to correctly account for the annual variations in mortality levels observed throughout France.

Incorporating these variations into the fitting of a forward-looking log-Poisson model poses no major difficulty, and a two-stage parameter estimation process enables us to use conventional likelihood maximization algorithms.

The results obtained show that the impact of this additional volatility is negligible on the central tendency indicators.

On the other hand, there is a material impact on the capital requirement associated with longevity risk, with just under 15% of this requirement being induced by the presence of this volatility. The remainder is associated with uncertainty about the trend in death rates.

Thus, while the main hazard associated with the construction of a prospective mortality table remains the uncertainty attached to the determination of the trend (see Juillard et al. [2008] and Juillard and Planchet [2006] for detailed analyses on this point), taking into account these short-term fluctuations in mortality levels provides a better understanding of the determinants of longevity risk.