3.1. The spatial–temporal autoregressive (STAR) model

Despite its popularity, a general VAR model cannot be directly fitted to a typical mortality data set for two reasons. First, a VAR model requires a stationary response variable, while mortality rates over time are clearly trending and therefore nonstationary. Second, the unknown parameters (p) exceed the observations (T) in a nonconstrained VAR framework. For instance, with an intermediate number of age groups (N), say 50, the $p = 50\times51 >> T$ issue will arise, given that observations are available only for a few dozens of years.

To resolve those issues, Li and Lu (Reference Li and Lu2017) propose the spatial–temporal autoregressive (STAR) model to study and forecast mortality data. On the temporal dimension, it considers the Granger causality and cointegration to resolve the stationarity problem. On the age dimension, the STAR model adopts the sparse spatial information to reduce the dimensionality of p with the follow specification:

(3.1) \begin{equation}\begin{split}y_{1,t} &= \alpha_1 + y_{1,t-1} + \epsilon_{1,t} \\[5pt]y_{2,t} &= \alpha_2 + \beta_{2,1} y_{1,t-1} + (1-\beta_{2,1})y_{2,t-1}+\epsilon_{2,t} \\[5pt]y_{i,t} &= \alpha_i + \beta_{i,i-2} y_{i-2,t-1}+ \beta_{i,i-1} y_{i-1,t-1}+ (1-\beta_{i,i-2}-\beta_{i,i-1})y_{i,t-1}+\epsilon_{i,t},\end{split}\end{equation}

where we let $y_{i,t}=\ln m_{i,t}$ for simplicity, $i=3,4,\ldots,N$ , and $t=1,2,\ldots,T$ . The residual sequence $\epsilon_{i,t}$ is similarly assumed as in the LC model.Footnote 1

Rewriting (3.1) in a VAR(1) form, we have that

(3.2) \begin{equation}\boldsymbol{{y}}_t = \boldsymbol{\alpha} + \boldsymbol{{B}} \boldsymbol{{y}}_{t-1} + \boldsymbol{\epsilon}_t,\end{equation}

where $\boldsymbol{{y}}_t$ = $(y_{1,t},y_{2,t},\ldots,y_{N,t})^\prime$ , $\boldsymbol{\alpha}$ = $(\alpha_{1},\alpha_{2},\ldots,\alpha_{N})^\prime$ , $\boldsymbol{\epsilon}_t$ = $(\epsilon_{1,t},\epsilon_{2,t},\ldots,\epsilon_{N,t})^\prime$ , and

(3.3) \begin{equation}\boldsymbol{{B}} = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}1&0&0&\cdots&\cdots\\[5pt]\beta_{2,1}&1-\beta_{2,1}&0&\cdots&\cdots\\[5pt]\beta_{3,1}&\beta_{3,2}&1-\beta_{3,1}-\beta_{3,2}&0&\cdots\\[5pt]\vdots&\vdots&\ddots&\ddots&0\\[5pt]\cdots&0&\beta_{N,N-2}&\beta_{N,N-1}&1-\beta_{N,N-2}-\beta_{N,N-1}\\[5pt]\end{array}\right]\end{equation}

To ensure invertibility and interpretability, it is required that $0<\beta_{i,i-2}<1$ , $0<\beta_{i,i-1}<1$ , and $0<\beta_{i,i-2}+\beta_{i,i-1}<1$ for all $i>2$ . When $i=2$ , we need that $0<\beta_{2,1}<1$ .

The interpretation of the STAR model parameters is straightforward. For age $i>2$ , $(1-\beta_{i,i-2}-\beta_{i,i-1})$ is the usual temporal effect of the first time lag $y_{i,t-1}$ on $y_{i,t}$ . As for the mortality practice, cohort effects are more important, which are measured by $\beta_{i,i-1}$ (the same cohort) and $\beta_{i,i-2}$ (the younger neighboring cohort). The relevant positive constraints not only ensure the invertibility (and thus stationarity) of $y_{i+1,t}-y_{i,t}$ but also the interpretability of temporal and cohort effects. A negative measure is usually unexplainable, and nor it is revealed in the existing literature.

In addition to the interpretability, the specification described in (3.2) and (3.3) has attractive statistical features. As shown in Li and Lu (Reference Li and Lu2017), with the assumption that all $y_{i,t}$ are I(1) (commonly used in mortality research), all neighboring age pairs $y_{i,t}$ and $y_{i+1,t}$ are cointegrated with order $(\!-\!1,1)$ . This approach successfully resolves the stationarity issue and results in age coherence. More specifically, cointegration suggests that $y_{i+1,t}-y_{i,t}$ is stationary, which leads to a constant long-run mean. Further, the total number of parameters is largely reduced from $N^2+N$ to $3N-3$ , which is feasible to solve under a general penalized least squares framework with closed-form solutions.

The forecasting of STAR is therefore performed in an iterative fashion, where

(3.4) \begin{equation}\begin{split}\widehat{\boldsymbol{{y}}}_{t+1} &= \widehat{\boldsymbol{\alpha}} + \widehat{\boldsymbol{{B}}} \boldsymbol{{y}}_{t}\\[5pt]\widehat{\boldsymbol{{y}}}_{t+h} &= \widehat{\boldsymbol{\alpha}} + \widehat{\boldsymbol{{B}}} \widehat{\boldsymbol{{y}}}_{t+h-1} \end{split}\end{equation}

and $h>1$ . Due to the established co-integration feature, $\hat{y}_{i,t+h}-\hat{y}_{j,t+h}$ is stationary for all $i,j=1,\ldots,N$ and $i\ne j$ , and $h\ge 1$ . This then straightforwardly ensures the age coherence in forecasting.

3.2. Age-coherent extensions of the LC model

Other than working with the VAR framework, a recent study by Gao and Shi (Reference Gao and Shi2021) proposes two effective age-coherent extensions of the LC model. Both are motivated by Li et al. (Reference Li, Lee and Gerland2013) to rotate $\hat{b}_{x,h}$ gradually over time. The rotation is supported by the fact that mortality decline decelerates at younger ages and accelerates at old ages (Li et al., Reference Li, Lee and Gerland2013). Also, to realize age coherence in the long run, Gao and Shi (Reference Gao and Shi2021) require that $\hat{b}_{x,h}$ eventually converge (rotate) to a constant $1/N$ for all ages (a flat line). The rationale of $1/N$ comes from the constraint that $\sum_{x=1}^N\hat{b}_{x,h}=1$ when $h \ge 1$ and $\hat{b}_{i,h}=\hat{b}_{j,h}$ when $h \rightarrow \infty$ .

Specifically, the first extension adopts the autoregressive (AR) framework. Under the AR(1) scenario, the time-varying $\hat{b}_{x,h}$ (represented by $\hat{b}^G_{x,h}$ ) is specified below.

(3.5) \begin{equation}\hat{b}^G_{x,h}={r^l_x}(\hat{b}^G_{x,h-1}-1/N)+1/N=(r^l_x)^h(\hat{b}_x-1/N)+1/N,\end{equation}

where $h\ge 1$ and $\hat{b}^G_{x,0}=\hat{b}_{x}$ . To ensure the stationarity, the coefficient $r^l_x$ strictly falls between 0 and 1. In such a case, $\hat{b}^G_{x,h}$ converges to the long-run mean $1/N$ geometrically. For a larger (smaller) $r^l_x$ , the speed of convergence is slower (faster).

The second extension employs the hyperbolic decay, which may achieve the final convergence slower than the AR specification. This concept is often refer to as long memory in time-series analysis for economic and financial practices (see, for example, Smallwood and Norrbin, Reference Smallwood and Norrbin2006; Ho and Shi, Reference Ho and Shi2020, among others). When assumed to decay hyperbolically, the time-varying $\hat{b}_{x,h}$ (represented by $\hat{b}^H_{x,h}$ ) is described by

(3.6) \begin{equation}\hat{b}^H_{x,h}=\delta_h(d^l_x)(\hat{b}_x-1/N)+1/N,\end{equation}

where

\begin{equation*}\delta_k(d^l_x)=\frac{k-1+d^l_x}{k}\delta _{k-1} (d^l_x)\mbox{ and }\delta _0 (d^l_x)=1.\end{equation*}

The hyperbolic parameter, or $d^l_x$ , needs to fall between 0 and 1 for a stationary rotation. In this case, since $\delta_{h} (d^l_x) \rightarrow 0$ when $h \rightarrow \infty$ , $\hat{b}^H_{x,h}$ will eventually converge to $1/N$ and thus leads to the age-coherent forecasts. The speed of decay is slower (faster) for larger (smaller) $d^l_x$ .

Note that both $r^l_x$ and $d^l_x$ are age dependent, and their estimation is not a trivial issue due to the large sample size on the age dimension. Gao and Shi (Reference Gao and Shi2021) employ the inversed Epanechnikov kernel to reduce the complexity of estimation. More specifically, let $\tau$ defined as a scaled index $x/N$ with $x \in \{1,\ldots,N\}$ , when evaluated at N, the inversed Epanechnikov kernel is determined by $1- K_{b^w}(\tau-1)=1-K(\frac{\tau-1}{{b^w}})$ , where $K(\tau-1)=0.75(1-(\tau-1)^2)_+$ . The parameter ${b^w}$ is the corresponding kernel bandwidth and constrained to fall in the range of (0, 1]. Hence, the parametric structures of $r^l_x$ and $d^l_x$ are simplified to follow the pattern of $1-K_{b^w}(\tau-1)$ across ages.

Once $r^l_x$ and $d^l_x$ become available, the forecast logged mortality rates of the LC extension with $\hat{b}^G_{x,h}$ (LC-G) and that with $\hat{b}^H_{x,h}$ (LC-H) are

(3.7) \begin{equation}\begin{array}{c}\ln \hat{m}^G_{x,T+h} =\hat{a}_x + \hat{b}^G_{x,h} \left( \hat{k}_T + h \hat{d}_k\right)\\[5pt] \;\,\ln \hat{m}^H_{x,T+h} =\hat{a}_x + \hat{b}^H_{x,h} \left( \hat{k}_T + h \hat{d}_k\right),\end{array}\end{equation}

where $\hat{a}_x$ and $\hat{k}_t$ are identical to those of the original LC model. In other words, the in-sample fits of LC-G and LC-H are the same as those described in (2.1). Recall that when $h \rightarrow \infty$ , $\hat{b}^G_{x,h}$ or $\hat{b}^H_{x,h}$ converges to $1/N$ . Thus, following (3.7), using either LC-G or LC-H, forecast rates of different ages will not differ infinitely and are thus age coherent.

3.3. The sparse VAR (SVAR) model and coherent extensions

Despite its attractive features, STAR’s parametric structure as of (3.3) is ad hoc and relatively inflexible. Specifically, the effects of all cohorts other than the same and next younger one are forced to be zero. This may limit a comprehensive analysis of the effects of all possible cohorts.

To address this, Guibert et al. (Reference Guibert, Lopez and Piette2019) recently work on the mortality improvements denoted by $\triangle y_{i,t}=y_{i,t}-y_{i,t-1}$ and propose a sparse VAR (SVAR) model. Since $y_{i,t}$ is assumed I(1), the differentiation resolves the stationarity issue of a VAR system. Moreover, the SVAR model adopts a pure data-driven method to select nonzero coefficients via the ENET penalty estimation, rather than working with the ad hoc structure as in (3.3). The SVAR model is specified below.

(3.8) \begin{equation}\triangle{\boldsymbol{{y}}}_t = \boldsymbol{{a}} + \sum_{v=1}^{V}\boldsymbol{{B}}_v \triangle{\boldsymbol{{y}}}_{t-v} + \boldsymbol{\epsilon}_t,\end{equation}

where $\triangle{\boldsymbol{{y}}}_t$ = $(\triangle y_{1,t},\triangle y_{2,t},\ldots,\triangle y_{N,t})^\prime$ , V is the preselected AR lags and the sparsity of $\boldsymbol{{B}}_v$ is determined by the ENET penalty without any constraints. The forecast is performed similarly to (3.4) as a usual VAR-type model for $\triangle{\boldsymbol{{y}}}_t$ , and the forecast for the original logged mortality rate is computed as $\widehat{\boldsymbol{{y}}}_{t+h} = \boldsymbol{{y}}_{t} + \sum_{l=1}^{h} \triangle{\widehat{\boldsymbol{{y}}}}_{t+l}$ .

However, the desirable feature of age coherence is lost when working with mortality improvements directly, as in the SVAR model. That is, as long as the estimated long-term mean of $\triangle y_{i+1,t}-\triangle y_{i,t}$ (denoted by $\triangle\hat{\mu}_{i+1}$ ) is not $O_p(1/h)$ , the long-run forecast ${\widehat{y}_{i+1,t+h}}-{\widehat{y}_{i,t+h}}={y_{i+1,t}}-{y_{i,t}}+h\triangle\hat{\mu}_{i+1}$ will still grow to reach infinity, when $h \rightarrow \infty$ .

To produce age-coherent forecasts, a recent study of Li and Shi (Reference Li and Shi2021) extends the approach of Gao and Shi (Reference Gao and Shi2021) to the SVAR framework. In particular, the intercept $\hat{\alpha}_{x}$ is allowed to be time-varying for age x, and Li and Shi (Reference Li and Shi2021) adopt the hyperbolic decay as follows:

\begin{equation*}\hat{\alpha}^H_{x,h}=\delta_h(d^s_x)(\hat{\alpha}_x-\hat{\alpha}^*_x)+\hat{\alpha}^*_x,\end{equation*}

where $\hat{\alpha}_x$ is the estimate of a usual SVAR model as of (3.8), and $\hat{\alpha}^*_x$ is the long-term mean of $\hat{\alpha}_x$ . This specification is named SVAR with hyperbolic decay, or SVAR-H model. A similar specification can then be derived following the geometric decay as in (3.5):

\begin{equation*}\hat{\alpha}^G_{x,h}={r^s_x}(\hat{\alpha}^G_{x,h-1}-\hat{\alpha}^*_x)+\hat{\alpha}^*_x,\end{equation*}

which is named SVAR with geometric decay, or SVAR-G model. The age-dependent $r^s_x$ and $d^s_x$ can be specified using the inversed Epanechnikov kernel, as for LC-G and LC-H models.

Forecasts of $\triangle{\widehat{\boldsymbol{{y}}}}_{T+h}$ using the SVAR-G or SVAR-H model are the same as in (3.4). Logged mortality rates can then be forecast by

\begin{equation*}\hat{y}_{x,T+h}=y_{x,T}+\sum_{k=1}^{h} \triangle \hat{y}_{x,T+k}.\end{equation*}

Remark 5. Since $\triangle y_{x,t}$ is stationary, $\hat{y}_{x,T+h}$ will approach $y_{x,T}+ h \hat{\mu}_{x}$ in the long run, where $\hat{\mu}_{x}$ is the estimated mean of $\triangle y_{x,t}$ . In a matrix form, it is easy to see that

(3.9) \begin{equation}\boldsymbol\mu = [I-\boldsymbol{{B}}_1 -\ldots-\boldsymbol{{B}}_V]^{-1} \boldsymbol \alpha.\end{equation}

Further, recall that when $h \rightarrow \infty$ , $\hat{\alpha}^G_{x,T+h}$ or $\hat{\alpha}^H_{x,T+h}$ converges to $\hat{\alpha}^*_x$ . To achieve age coherence, $\hat{\alpha}^*_x$ is obtained such that all elements of $\boldsymbol{\mu}$ are identical. In practice, estimates of $\hat{\alpha}^*_x$ are produced as follows:

\begin{equation*}\hat{\boldsymbol \alpha}^*_x= [I-\hat{\boldsymbol{{B}}}_1 -\ldots-\hat{\boldsymbol{{B}}}_V] \hat{\boldsymbol\mu}^*,\end{equation*}

where $\hat{\boldsymbol{{B}}}_v$ are in-sample estimates of a usual SVAR model, and $\hat{\boldsymbol\mu}^*$ consists of identical elements of $\hat{\mu}^*$ , which can be estimated as the average of $\hat{\mu}_x$ over all ages. $\hat{\mu}_x$ is obtained from the in-sample estimates of a usual SVAR model with (3.9).

Note that there are three tuning parameters in both SVAR-H and SVAR-G models: the sparsity penalty parameter ( $\lambda_s$ ), the hyperbolic or geometric decay parameter of the first age group ( $d^s_1$ or $r^s_1$ ), and the kernel bandwidth ( ${b^w}$ ). See Section 3.4 for how those parameters are selected.

3.4. Tuning parameters selection

All the five age-coherent mortality models described in this section require some preselected tuning parameters, either before the model is fitted (STAR, SVAR-H, and SVAR-G) or at the forecasting stage (LC-H, LC-G, SVAR-H, and SVAR-G). Due to the time-series nature, the usual cross-validation technique to select tuning parameters is not directly applicable.

In practice, a popular strategy for time-series data is to adopt the expanding-window approach explained in Hyndman and Athanasopoulos (Reference Hyndman and Athanasopoulos2018) to collect short-term out-of-sample forecasts. However, since age coherence is a long-term property, we follow Gao and Shi (Reference Gao and Shi2021) and adopt the hold-out-sample approach to select tuning parameters in all cases. Specifically, the selection aims to minimize

(3.10) \begin{equation}\text{RMSFE}=\sqrt{\frac{1}{N(T/4)}\sum_{i=1}^{N}\sum_{h=1}^{T/4}\left(\ln \hat{m}_{i,3T/4+h}-\ln {m}_{i,3T/4+h}\right)^2},\end{equation}

where RMSFE is the root of mean squared forecasting errors, and the evaluation period is given by the last fourth ( $[3T/4+1,T]$ ) of the data in our study. A high-level summary of the selection procedure is listed in Table 1 of Supplement Material for each model.

Table 1. Out-of-sample forecasting performance.

4.1. Background: A solution from global minimum variance portfolio

In his seminal work, Markowitz (Reference Markowitz1952) discusses a solution to achieve global minimum variance for a portfolio, or the GMVP solution. The idea can be straightforwardly applied to forecasting error minimization using an MA model. Suppose that there exist a finite number ( $J>2$ ) of models, the forecasting error of an MA approach is

\begin{equation*}{e}_{M,t} = w_1 {e}_{1,t} + w_2 {e}_{2,t}+\ldots + w_J {e}_{J,t}, \end{equation*}

where $e_{j,t}$ is the forecast error at time t for model j ( $j \in \{1,2,\ldots,J\}$ ), $w_j$ is the assigned weight for the jth model and ${e}_{M,t}$ is then the weighted summation of forecast errors of all J models, or the forecast error of the MA approach. Denote that $\boldsymbol{{w}}=(w_1,\ldots,w_J)^\prime$ , the GMVP solution of the optimal ${\boldsymbol{{w}}}^*$ is

\begin{equation*}{\boldsymbol{{w}}}^*=\boldsymbol{\Sigma}_e^{-1} \boldsymbol\iota (\boldsymbol{\iota}^\prime \boldsymbol{\Sigma}_e^{-1} \boldsymbol{\iota})^{-1},\end{equation*}

where $\boldsymbol{\Sigma}_e$ is a $J \times J$ variance-covariance matrix of forecast errors, and $\boldsymbol\iota$ is a $J \times 1$ vector of ones.

Remark 6. Note that when forecast errors of some models are strongly (and positively) correlated, it is possible to have negative weights in ${\boldsymbol{{w}}}^*$ . Under a financial portfolio optimization scenario, this means a short position should be held for some assets. In the mortality case, although negative weights are sensible from a model combination perspective, their interpretation is not as straightforward. Also, negative weights might cause instability in long-term forecasting. To avoid this, the non-negative constraints should be added to the optimization issue that is usually faced in a GMVP problem:

\begin{equation*}\begin{matrix}\underset{\boldsymbol{{w}}}\min\;\;\;\; & \;\;\;\;\boldsymbol{{w}}^\prime \boldsymbol{\Sigma}_e \boldsymbol{{w}\;\;\;\;} \\[5pt]\mbox{subject to} & \boldsymbol\iota^\prime \boldsymbol{{w}} =1\\[3pt]& \boldsymbol{{w}} \ge \textbf{0}.\end{matrix}\end{equation*}

This is a quadratic optimization issue with equality and inequality constraints, which can be solved via iteration-based algorithms (Gill et al., Reference Gill, Murray and Wright2019).

4.2. An MA approach for mortality forecasting

For an MA model used in mortality forecasting, the most essential issue is the objective function to be employed in the optimization. This is because that for the original GMVP solution, there are two major issues. First, since only the global risk is considered, GMVP may result in large weights in few models. This is also known as the “winner-take-all” problem and may negatively affect diversity in an MA approach. Second, the estimate is highly sensitive to parameter estimation errors of the covariance matrix. Consequently, an inappropriate covariance error matrix may lead to unrobust weights in an MA model.

In this paper, diversity is ensured in two ways. First, we employ both age-coherent and age-incoherent models, which include a total of ten specifications: STAR, LC-H, LC-G, SVAR-H, SVAR-H, LC, SVAR, APC, HU, and RH. Second, different from a single weight for all ages used in existing studies (see, for example, Shang, Reference Shang2012; Kontis et al., Reference Kontis, Bennett, Mathers, Li, Foreman and Ezzati2017; Kessy et al., Reference Kessy, Sherris, Villegas and Ziveyi2021, among others), age-specific weights are employed in all cases. As preliminarily evidenced in Figure 1(b), not all models can perform uniformly well for all age groups. This setting then effectively prevents the existence of one or few dominant models.

To avoid the potential overfitting issue of in-sample errors, we construct $\Sigma_{e}$ using the out-of-sample errors only. Specifically, the objective function will consist of out-of-sample mean squared forecasting errors (MSFE). Since it is well known that MSFE incorporates both process risk and parameter risk, this will effectively address the second problem of GMVP. In addition, we adopt two pieces of prior information in the optimization to further improve reasonableness in mortality forecasting. First, as argued above, the MA approach should (asymptotically) result in age-coherent forecasts. Therefore, the weights assigned to age-incoherent models should be asymptotically approaching zero. A penalty is imposed in the objective function for this aim. Second, mortality rates across neighboring ages tend to demonstrate similar dynamics. However, without further constraints, abrupt changes of forecast rates will be unavoidable, even for neighboring ages, due to potential randomness introduced by small samples. Consequently, inspired by existing studies such as Li and Lu (Reference Li and Lu2017), a smoothness penalty is additionally introduced to the objective function.

The optimization problem to choose optimal weights of an MA mortality model is then stated below. Note that for a mortality data set with N age groups estimated by J models (with the total $J_c$ age-coherent models ordered first and the rest $J_n$ age-incoherent models ordered last), the weight vector to be estimated is ${\boldsymbol{{w}}}=(\boldsymbol{{w}}^\prime_1,\ldots,\boldsymbol{{w}}^\prime_N)^\prime$ (the dimensionality is $(NJ)\times 1$ ), with $\boldsymbol{{w}}_x=(w_{x,1},\ldots,w_{x,J})^\prime$ . To further allow for the age coherence penalty, smoothness penalty and non-negative weights, we have the following optimization problem:

(4.1) \begin{equation} \begin{matrix} \underset{\mathcal{\boldsymbol{{w}}}}\min\;\;\;\; & \boldsymbol{{w}}^\prime \boldsymbol{{\Omega}}_{e} \boldsymbol{{w}}+ + \lambda_1\boldsymbol{{w}}^\prime \boldsymbol{{C}}_1 \boldsymbol{{w}} + \lambda_2\boldsymbol{{w}}^\prime \boldsymbol{{C}}_2 \boldsymbol{{w}}\\[5pt] \mbox{subject to}\;\;\;\; & \boldsymbol{\mathcal{I}} \boldsymbol{{w}} =\textbf{1}\\[5pt] & \boldsymbol{{w}} \ge \textbf{0}, \end{matrix}\end{equation}

where $\boldsymbol{{\Omega}}_{e}$ is an $(NJ)\times(NJ)$ matrix with $\boldsymbol{\Sigma}_{e_x}$ on the diagonal blocks and elsewhere 0. $\boldsymbol{{C}}_1$ and $\boldsymbol{{C}}_2$ are both $(NJ)\times(NJ)$ ancillary matrices to specify the quadratic age coherence penalty and smoothness penalty, respectively. $\lambda_1$ and $\lambda_2$ are the associated tuning parameters. $\boldsymbol{\mathcal{I}}$ is another $(NJ)\times(NJ)$ ancillary matrix to impose the equality constraint for each age. Specifically, $\boldsymbol{{C}}_1$ is a diagonal matrix with the following diagonal items:

\begin{equation*}(\overbrace{\underbrace{0,\ldots,0}_{\text{$J_c$}},\underbrace{1,\ldots,1}_{\text{$J_n$}}\ldots\underbrace{0,\ldots,0}_{\text{$J_c$}},\underbrace{1,\ldots,1}_{\text{$J_n$}}}^{\text{$N\times J$}}).\end{equation*}

The xjth column of $\boldsymbol{{C}}_2$ is

\begin{equation*}\boldsymbol{{c}}_{xj}=(\overbrace{\underbrace{0,..,0}_{\text{$j-1$}},-1,\underbrace{0,..,0}_{\text{$j-1$},-1,\ldots}}^{\text{$xj-1$}},1, \underbrace{0,..,0}_{\text{$j-1$}},\ldots )^\prime.\end{equation*}

That is, there are $(N-1)$ (−1)’s and one 1 in each column, with all other elements being 0. An example is provided below to illustrate the specification when $N=2$ , $J_c=1$ and $J_n=1$ . Also,

\begin{equation*}\boldsymbol{\mathcal{I}}=\begin{pmatrix} \boldsymbol{\iota}^\prime\;\;\;\; & \textbf{0}\;\;\;\; & \cdots\;\;\;\; &\textbf{0}\\[5pt] \textbf{0}\;\;\;\;&\boldsymbol{\iota}^\prime \;\;\;\;& \cdots \;\;\;\;&\textbf{0} \\[5pt] \cdots\;\;\;\; & \cdots\;\;\;\; & \cdots\;\;\;\; & \cdots\\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \cdots\;\;\;\; &\boldsymbol{\iota}^\prime\\[5pt]\end{pmatrix},\end{equation*}

where $\boldsymbol{\iota}$ is a $J \times 1$ vector consisting of 1, as defined previously. The $\boldsymbol{{w}}^*$ can then be obtained using a standard quadratic programming algorithm, with preselected $\lambda_1$ and $\lambda_2$ (Goldfarb and Idnani, Reference Goldfarb and Idnani1983). Note that the coherence penalty term is

\begin{equation*}\lambda_1\boldsymbol{{w}}^\prime \boldsymbol{{C}}_1 \boldsymbol{{w}}=\lambda_1\sum_{x=1}^N\sum_{j=1}^{J_n} w_{x,j}^2,\end{equation*}

which reduces the weights for the $J_n$ age-incoherent models with the influence of penalty $\lambda_1$ , and the smoothness penalty term is

\begin{equation*}\lambda_2\boldsymbol{{w}}^\prime \boldsymbol{{C}}_2 \boldsymbol{{w}}=\lambda_2\sum_{x=2}^N\sum_{j=1}^J (w_{x-1,j}-w_{x,j})^2,\end{equation*}

which increases the smoothness of assigned weights for the same model between adjacent ages for a larger $\lambda_2$ .

Example 1. Consider a simple case with $N=2$ ages, $J_c=1$ and $J_n=1$ models, we have that

\begin{equation*}\boldsymbol{{\Omega}}_{e}=\begin{pmatrix} \boldsymbol{\Sigma}_{e_1}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \boldsymbol{\Sigma}_{e_2} \\[5pt] \end{pmatrix} = \begin{pmatrix} \sigma_{1,1}^2 \;\;\;\;& \sigma_{1,12} \;\;\;\;& 0 \;\;\;\;& 0\\[5pt] \sigma_{1,12} \;\;\;\;& \sigma_{1,2}^2 \;\;\;\;& 0 \;\;\;\;& 0\\[5pt] 0\;\;\;\; & 0\;\;\;\; & \sigma_{2,1}^2 \;\;\;\;& \sigma_{2,12}\\[5pt] 0\;\;\;\; & 0\;\;\;\; & \sigma_{2,12}\;\;\;\; & \sigma_{2,2}^2\\[5pt] \end{pmatrix},\end{equation*}

\begin{equation*}\boldsymbol{{C}}_1= \begin{pmatrix} 0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & 0\\[5pt] 0\;\;\;\; & 1\;\;\;\; & 0\;\;\;\; & 0\\[5pt] 0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & 0\\[5pt] 0\;\;\;\; & 0\;\;\;\; & 0\;\;\;\; & 1\\[5pt] \end{pmatrix} \text{ and } \boldsymbol{{C}}_2= \begin{pmatrix} 1\;\;\;\; & 0\;\;\;\; & -1\;\;\;\; & 0\\[5pt] 0\;\;\;\; & 1\;\;\;\; & 0\;\;\;\; & -1\\[5pt] -1\;\;\;\; & 0\;\;\;\; & 1\;\;\;\; & 0\\[5pt] 0\;\;\;\; & -1\;\;\;\; & 0\;\;\;\; & 1\\[5pt] \end{pmatrix}. \end{equation*}

The optimization problem stated in (4.1) then reduces to

\begin{equation*}\begin{matrix}\underset{\mathcal{\boldsymbol{{w}}}}\min & w_{1,1}^2\sigma_{1,1}^2+w_{1,2}^2\sigma_{1,2}^2+ w_{2,1}^2\sigma_{2,1}^2 + w_{2,2}^2\sigma_{2,2}^2+2w_{1,1}w_{1,2}\sigma_{1,12}+2w_{2,1}w_{2,2}\sigma_{2,12}\\[5pt]&+\lambda_1(w_{1,2}^2+w_{2,2}^2)+\lambda_2 [(w_{1,1}-w_{2,1})^2+(w_{1,2}-w_{2,2})^2]\\[5pt] & w_{1,1}+w_{1,2} =1\\[5pt]{subject \; to}& w_{2,1}+w_{2,2} =1\\[5pt]& w_{1,1}, w_{1,2}, w_{2,1}, w_{2,2} \ge 0.\end{matrix}\end{equation*}

With a larger (smaller) $\lambda_1$ , weights of the second model (assumed age-incoherent) will be smaller (larger). Similarly, with a larger (smaller) $\lambda_2$ , the $\hat{w}_{1,1}$ and $\hat{w}_{1,2}$ will change to $\hat{w}_{2,1}$ and $\hat{w}_{2,2}$ more (less) smoothly.

Denote that the h-step-ahead forecast of the age x’s logged mortality rate by the jth model as $\hat{y}_{x,j,T+h}$ , the forecast of MA is therefore

\begin{equation*}\hat{y}_{x,M,T+h}=\sum_{j=1}^{J}\hat{w}_{x,j}\hat{y}_{x,j,T+h}.\end{equation*}

We now demonstrate that this forecast also achieves the desirable age coherence in an asymptotic fashion.

Finally, the selection of $\lambda_1$ and $\lambda_2$ can still be performed using the hold-out-sample strategy as explained in Section 3.4. The detailed steps to forecast logged mortality rates using the proposed MA approach are listed below.

1. 1. With the training period $[1,\ldots,T/2]$ , fit each of the J models that are included in the ensemble;

2. 2. Collect out-of-sample forecast errors over $[T/2+1,\ldots,3T/4]$ , and calculate the sample estimate $\widehat{\boldsymbol{\Sigma}}_{e}$ ;

3. 3. Perform a grid search of $\lambda_1$ and $\lambda_2$ : for each candidate, obtain estimates of weights as in (4.1) using $\widehat{\boldsymbol{\Sigma}}_{e}$ obtained in Step 2, fit all individual models over the full training sample $[1,\ldots,3T/4]$ , produce the forecasts over the test period $[3T/4+1,\ldots,T]$ , and calculate the corresponding RMSFE of the MA model using the estimated weights and forecasts of individual models as in (3.10);

4. 4. Select the optimal $\lambda_1$ and $\lambda_2$ as that minimize the RMSFE;

5. 5. Repeat steps 1-2 to obtain final estimated weights, with selected tuning parameters over training period $[1,\ldots,3T/4]$ and $\widehat{\boldsymbol{\Sigma}}_{e}$ estimated with out-of-sample forecast errors over $[3T/4+1,\ldots,T]$ ; and

6. 6. Conduct the out-of-sample forecasting ( $[T+1,\ldots,T+h]$ ) with each of the J models, and the forecast rate of MA is their weighted summation using weights estimated at step 5.

To study the uncertainties of the forecast, we may follow the usual bootstrap method as explained in Chang and Shi (Reference Chang and Shi2021). Simply speaking, in-sample errors can be bootstrapped to produce replicates, which are then fitted following the procedure explained above to forecast new rates. The 2.5th and 97.5th percentiles out of those forecasts can then construct the 95th prediction interval (PI).

5.1. Baseline results

To compare the forecasting performance across all models, we follow existing studies such as Li and Lu (Reference Li and Lu2017) and employ the RMSFE examining all age groups and forecasting steps:

\begin{equation*}RMSFE = \sqrt{\dfrac{1}{101 \times h}\sum_{i=1}^{h}\sum_{x=0}^{100}(y_{x,T+i} - \hat{y}_{x,T+i})^2}.\end{equation*}

In the baseline case, we select $h=10$ , such that the training sample of 1950–2006 is fitted, and the out-of-sample forecasts are produced over 2007–2016 for each population. The ten individual models: STAR, LC-H, LC-G, STAR-H, STAR-G, LC, SVAR, APC, HU, and RH are first estimated, following necessary tuning parameter selection procedures described in Section 3.4. Weights of those models are then chosen following steps listed in Section 4.2. All forecasts of each of the eleven models are then collected to compute their corresponding RMSFE. This procedure is executed for all the ten examined populations.

All specific RMSFEs are presented in Table 1. To facilitate the comparison, the ranking of each model is displayed in Figure 2, and a black color indicates a higher place. On average, it can be seen that the five age-coherent models beats the five age-incoherent counterparts. More importantly, contrasting the six age-coherent sets of forecasts, those of MA outperform individual models in eight out of the ten populations and are the second or third best for the rest. This strongly supports the dominant forecasting accuracy using our proposed MA technique. As for the specific metrics, on average, the RMSFE of MA is 0.1707, which is almost 25% smaller than that of LC (0.2360). A summary of those baseline RMSFEs can be found in Table 2, which is presented when the MCS analysis is performed.

Figure 2. Heatmap of model forecasting accuracy rankings under the baseline scenario.

Table 2. Model confidence sets: selected models.

We now compare the individual models included in the ensemble. On average, STAR produces the lowest RMSFE, and results of the other four age-coherent models are relatively close to each other. Specifically, the average ranking of STAR is 3.6, and those of LC-H, SVAR-H, and SVAR-G are close to 5, whereas LC-G ranks somewhat lower at 7. Further, RMSFEs of SVAR-H and SVAR-G are much similar throughout all cases. Those of LC-H and LC-G are also not too far away from each other, although the difference is visually larger than those in SVAR-H and SVAR-G. For the five age-incoherent models, SVAR and HU result in similar overall performance, whereas LC, APC, and RH rank at the end of all considered models.

5.2. Sensitivity analyses

To check the robustness of baseline results, we consider three sets of sensitivity analyses individually:

1. 1. Forecast steps are increased to $h=16$ , as examined in Gao and Shi (Reference Gao and Shi2021);

2. 2. The starting year is truncated to 1960; and

3. 3. The age range is reduced to 0–89.

The ranking of all models examined in each set is plotted in Figure 3(a)–(c), and corresponding RMSFEs are summarized in Table 2 in the Supplement Materials. The ranking plot suggests that our proposed MA approach works consistently well by altering individually the forecasting horizon, weights construction, sample period, and age groups. It almost uniformly ranks among the top 3 models and produces the smallest average RMSFE in all scenarios. Nevertheless, it is worth noting that STAR model ranks the second best in all cases, whereas LC, APC, and RH are the three least preferred models. The overall performance of the five age-coherent models is uniformly better than that of the five age-incoherent models, by averaging the RMSFEs across relevant models and populations.

Figure 3. Heatmaps of model forecasting accuracy rankings under additional scenarios.

In summary, we employ the popular criterion RMSFE and demonstrate that the proposed MA approach can overall improve the forecasting performance of all models included in the ensemble. This conclusion is robust when various sensitivity analyses are conducted with individual adjustment being made.

5.3. Comparison with the model confidence set and Bayesian modeling average approaches

In their seminal work, Hansen et al. (Reference Hansen, Lunde and Nason2011) propose an MCS approach that chooses the subset of superior models out of a range of candidate. Those superior models are identified by testing if they can be assumed to present equal predictive ability at a given confidence level. Such a test may be performed for any loss functions, including the popular squared losses employed in this paper. Briefly speaking, the MCS test works using bootstrap samples sequentially, by eliminating the worst model at each stage. The procedure stops until the null hypothesis of equal predictive ability of the remaining model cannot be rejected.

To employ the MCS for mortality forecasting, we produce squared forecast errors over the last 1/4 of the training period (i.e., 1993–2006) for each individual model. Similar to Kessy et al. (Reference Kessy, Sherris, Villegas and Ziveyi2021), we choose a confidence level of 10%, and the selected models for each population are presented in Table 2. It is interesting to note that LC-H and LC-G are the two mostly selected model in the MCS, followed by STAR. On the other hand, LC, APC, and RH models are not picked in any of the ten populations.

Assuming the equal weights for all ages and selected models, we construct this new MA approach via the MCS method. Together with the baseline results, the RMSFEs are summarized in Table 3, and individual ranks are illustrated in Figure 3(d). Clearly, our proposed MA approach beats the MCS according to all the five metrics presented in Table 3. For RMSFE of individual populations, our approach results in more accurate forecasts in seven out of ten countries. Additionally, given that MCS significantly reduces the number of selected models in the ensemble, we conclude that our approach is preferred for both the achieved diversity and resulted forecasting accuracy.

Table 3. Summary of all investigated models including the MCS.

Note: bold numbers are the smallest quantity for each statistic across the twelve models.

Alternative to the GMVP and MCS, another potentially useful model averaging approach is the Bayesian modeling average (BMA). Specifically, weights of the BMA are essentially posterior model probabilities. Detailed calculations of those weights can be found in Raftery et al. (Reference Raftery, Madigan and Hoeting1997), and a full BMA usually requires a computational intensive approach, such as the Markov chain Monte Carlo. Fortunately, in practice, seminal works including Bates and Granger (Reference Bates and Granger1969) have suggested that weights of a full BMA can be well approximated by employing a much simpler Akaike or Bayesian information criterion. Applications of such an approximation to mortality data can be found in Shang (Reference Shang2012). However, Wagenmakers and Farrell (Reference Wagenmakers and Farrell2004) point out that if one model in the ensemble is dominant according to the used information criterion, the posterior probability of a BMA may tend to be close to one, also known as the “winner-take-all” issue. Consequently, we consider an alternative BMA technique discussed in Kontis et al. (Reference Kontis, Bennett, Mathers, Li, Foreman and Ezzati2017), and age-specific weights are approximated by

\begin{equation*}\hat{w}^b_{x,j}\approx \dfrac{e^{-0.5|PB_{x,j}|}}{\sum_{j=1}^{J}e^{-0.5|PB_{x,j}|}},\end{equation*}

where $PB_{x,j}=\sum_{i=1}^{h}(\hat{y}_{x,j,T+i}-y_{x,j,T+i})/h$ is the projection bias for logged mortality rates of age x produced by model j. Intuitively, due to the closeness of mortality forecasts, no $PB_{x,j}$ could be universally small or close to 0 to cause a dominant weight in a BMA mortality model.

Using the age-specific BMA weights of all models, the RMSFE results are summarized in the last row of Table 3. Overall, the results of BMA outperform all but those of our MA approach. Two conclusions could be drawn. First, BMA forecasts could serve as an additional robustness check to our baseline results. The altered factor is the weighting strategy in the model averaging. Second, the overall outperformance of MA over BMA supports the superior effectiveness of our proposed approach. Despite that no dominant weights may exist in the employed BMA, such an approach is “deterministic” and formula based. Thus, it cannot cope with the coherent and smoothness penalties. Consequently, the desirable asymptotic age coherence will be lost for the resulting forecasts, and abrupt changes in mortality forecasts are inevitable, even for neighboring age groups. This might also contribute to its observed inferior forecasting performance to our MA approach.

6.1. Forecasting results

Instead of examining a single metric, we employ the RMSFE over the age dimension to compare forecasting performance of each model. Consistent with our baseline setting discussed in Section 5.1, the training period of 1950–2006 is employed to fit in-sample estimates, and out-of-sample forecasts are calculated for the test period of 2007–2016.Footnote 3 At each of the 101 one-year age groups (0–100), we calculate the age-specific RMSFE as

\begin{equation*}RMSFE_{x} = \sqrt{\dfrac{1}{10}\sum_{h=1}^{10}(y_{x,T+h} - \hat{y}_{x,T+h})^2}.\end{equation*}

The results are summarized in Table 4.

Table 4. Summary of RMSFE over age groups for the Italian population.

Note: bold numbers are the smallest quantity for each statistic across all models.

In Table 4, it can be seen that the proposed MA approach still leads to the smallest average RMSFE over ages, with the narrowest variation. Specifically, the mean (standard deviation) of MA’s RMSFE is 0.0859 (0.0604), over 50% (80%) smaller than that of LC. Thus, for the Italian population, MA can improve the forecasting accuracy of individuals models in the ensemble.

We explore out-of-sample interval estimates in Figure 4, where the data investigated are the logged mortality rates averaged out of all ages. The results of MA and SVAR-H, the top two best models as evidenced in Table 4 are contrasted, and the 95% PIs of MA and SVAR-H are constructed using bootstrap replicates as outlined in Section 4.2. Clearly, compared to the mean estimates of SVAR-H, those of MA are overall closer to the true values spanning 2007–2016. This is consistent with our observations in Table 4. Also, although both PIs manage to cover the true values at all forecasting steps, the width of MA’s PI is uniformly narrower than that of SVAR-H. This implies that MA is potentially more efficient than the SVAR-H model. Such a finding is consistent with our discussion in Section 3. Specially, MA with weights estimated via the forecast error optimization may improve the forecasting uncertainty over all individual models.

Figure 4. Forecast versus actual Italian log mortality rates: 2001–2016.

Finally, we display the long-term forecasts as of 2051, a 35-step-ahead horizon in Figure 5. The full data of 1950–2016 are fitted to provide the estimates of parameters, and forecasts of MA and LC are generated and contrasted. Compared to the true values of 2016, forecasts of LC in 2051 demonstrate little improvements on both the very old ages and “accidental hump” (ages 20–25). Such results are against the observations in Figure 1(a), which suggests significant improvements even for those ages. In contrast, due to its age-coherence feature, forecasts of MA show more consistent improvements over all ages. Also, even for a considerably long period of 35 years, no abrupt changes are observed at any ages for the MA model. This supports the effectiveness of the smoothness penalty of $\lambda_2$ employed in (4.1).

Figure 5. Actual and Long-term forecast Italian log mortality rates: 2016 vs 2051.

6.2. Validation of the proposed weighting method

In Figure 1(b), we present the $RMSFE_x$ of age-incoherent models over the out-of-sample period 1993–2006. Those corresponding forecast errors are the core part in the loss to be minimized in (4.1). Consistent with this, we present $RMSFE_x$ of age-coherent models over the same period in Figure 6. To facilitate the comparison, RMSFEs of five age-incoherent models are grouped together, and their average value is also plotted (denoted as NC). Summarizing from the two figures, there are two major findings. First, as stated in Section 1, no single model can uniformly outperform the rest across all age groups. This validates the necessity of age-specific weights as employed in our MA approach. Second, the $RMSFE_x$ of age-incoherent model is overall higher than that of age-coherent models. Specifically, age-incoherent $RMSFE_x$ is quite high over young ages 0–25 and is still among the worst for ages beyond 50. This suggests that with the imposed coherent penalty $\lambda_1$ in (4.1), the estimated weights of age-incoherent models may be much lower than those of age-coherent models.

Figure 6. $RMSFE_x$ of fitted models for Italian population: 1993–2006.

In Figure 7(a), we plot the final estimated weights of the MA model. The dynamics are much consistent with our findings in Figure 6. First, the total weights of age-incoherent models are small, due to their unfavorable RMSFE across ages discussed above. Also, fitted weights demonstrate much consistent patterns as their $RMSFE_x$ . For instance, the weights of SVAR-H and SVAR-G models are increasing over the oldest age groups, for their smallest RMSFEs at those ages. Additionally, all weights vary smoothly across ages, due to the impact of smoothness penalty $\lambda_2$ . More importantly, it is verified that no single model dominates in the identified age-specific weights, which supports the diversity realized by our MA approach. Last but not least, since SVAR-H and SVAR-G models produce much similar $RMSFE_x$ at all ages, their corresponding weights are close to each other.

We now contrast the MA results to those of BMA. Recall that in Section 5.3, an issue is identified for BMA, such that abrupt changes are unavoidable, due to the formula-based nature of BMA. To see this, we plot estimated weights of BMA in Figure 7(a). Clearly, those weights are identical to the proportions of $RMSFE_x$ of an individual model taking in the total RMSFEs at the same age, as shown in Figures 1(b) and 6. Consequently, rough dynamics are observed for all fitted weights, and sharp changes are observed in multiple ages. For instance, the weight of LC model in BMA drops from 0.1 to 0.05 from age 0 to 1, and then bounces back to 0.09 for age 2. Further, despite their overall unfavorable performance, age-incoherent models are still assigned considerably large weights. Averaging over all ages, weights of each model are close to 0.1, the value of which would have been assumed for a simple average method. As for the forecasting accuracy, the $RMSFE_x$ estimated by the BMA approach is reported in Table 4, the results are outperformed by the proposed MA counterparts.

In conclusion, to realize the coherent and smoothness penalties, a squared loss on weights will need to be implemented in the objective function. Consequently, the GMVP structure as adopted in (4.1) is a simple but effective approach, compared to other alternatives such as the BMA.Footnote 4 Besides, the incapability to incorporate those penalties will be faced by a novel approach proposed in Kessy et al. (Reference Kessy, Sherris, Villegas and Ziveyi2021). Although this stacked regression employ an Elastic NET (ENET) loss, the inclusion of the smoothness penalty is not straightforward. This is due to the fact that the smoothness penalty is not a standard L2-type loss, and a more complicated computational approach will need to be implemented (Chang and Shi, Reference Chang and Shi2021). The imposed equality and inequality constraints will further increase the computational cost.

Figure 7. Fitted weights of the ensemble average and Bayesian model average approaches.

6.3. A financial application: Fixed-term annuity pricing

To demonstrate the practical usefulness of the MA approach, a financial application is presented in this section. Following Shi (Reference Shi2022b), we consider the pricing of fixed-term annuities using the MA technique. This product has attracted a growing number of policyholders world-widely, especially those planning for their retirements. Comparing to the lifetime annuities, fixed-term products pay a predetermined and guaranteed income of higher level with deferrable options (Shang and Haberman, Reference Shang and Haberman2017). A cohort approach, as adopted in Shi (Reference Shi2022b), is employed here to price annuities. To be consistent with our previous analyses, the maximal survival age is limited to 100.

The specific pricing scheme is explained below. First, the $\tau$ year survival probability of a person aged x at $t = 0$ is

\begin{equation*} {_\tau}p_x = \prod_{j=1}^\tau {_1} p_{x+j-1},\end{equation*}

where $\tau {_1} p_{x+j-1} = e^{-m_{x+j-1}}$ and $m_{x+j-1}$ can be obtained from mortality forecasts. Hence, an annuity with a T-year maturity and written for an x-year-old person including a benefit $1 per year and conditional on the survival is priced as

\begin{equation*} a^{T}_x (m_{x,{1:T}}) = \sum_{\tau=1}^T P_B(0,\tau)E(I(T_x>\tau)|m_{x,{1:T}}),\end{equation*}

where $I(\cdot)$ is the indication function, $T_x$ is the survival time and $P_B(0,\tau)$ is the price of a $\tau$ -year bond with interest rate set to the yield of annuity. Since $E(I(T_x>\tau)|m_{x,{1:T}})={_\tau}p_x(m_{x,{1:T}})$ , the fixed-term annuity price is fully determined by the underlying yield and mortality rates. Thus, to accurately price the product, it is critical to providing precise forecasts of $m_{x,{1:T}}$ , potentially over a long period. That is, the mortality experiences of policyholders need to be accurately projected to minimize the mis-pricing risks.

For illustration purpose, we employ up to 35-step-ahead forecasts (from 2017 to 2051) of the proposed MA model. The estimated prices of fixed-term annuities as of 2016 are presented in Table 5. The calculation is conducted for the Italian population up to the 30-year-maturity, starting from age 65. We follow Shi (Reference Shi2022b) and present results of the four ages: 65, 70, 75 and 80. Both point and 95% interval estimates are reported. For simplicity, a constant interest rate of 3% is assumed in all case. Thus, the only source of uncertainty is limited to the precision of mortality forecasting. Evidenced in Figure 4, the proposed MA technique is more efficient in providing PIs. This supports the small widths observed in Table 5. As argued in Fung et al. (Reference Fung, Peters and Shevchenko2015), underpricing of annuities as small as by 0.1% can lead to dramatic shortfall in reserving with a large portfolio. Consequently, precisely and efficiently determine the uncertainty of premium rate is critical for insurers to optimizing their reserves, such that the associated ruin probability is minimized.

Table 5. Predicted fix-term annuity prices for the Italian population.

Note: this table displays the forecast fix-term annuity price for the average population as of 2016. The forecast mortality rates range from 2017 to 2051. LB and UB stand for the 2.5th and 97.5th percentiles of the prediction interval, respectively. Mean is the point/mean forecast price. T is the maturity term. Value in bracket is the percentage difference compared to the forecast mean annuity price. We only consider contracts with maturity so that age + maturity $\le$ 100.

Specifically, a prospective lifetable constructed using the mortality forecast by the MA approach may largely reduce the longevity risks. As demonstrated in Figure 7, an age-incoherent model like LC may underestimate mortality improvements of old ages. For pension products, a corresponding prospective lifetable will then result in an underestimated reserve and thus increase the ruin probability. Using an age-coherent model like the MA approach to produce such a prospective lifetable will help address this issue.