3.1. Forecasting mortality with GAPC models

In the general setting of GAPC mortality models, the time evolution of mortality is determined by the period effects $\kappa_t^{(i)}$ , $i=1,\ldots,N$ , and the cohort effect $\gamma_{t-x}$ . To forecast mortality rates, the standard approach (Villegas et al. Reference Villegas, Kaishev and Millossovich2018) is to assume that the period indices follow a multivariate random walk so that

(3.5) \begin{equation}\boldsymbol{\kappa}_t = \boldsymbol{\delta} + \boldsymbol{\kappa}_{t-1} + \boldsymbol{\xi}_t^\kappa, \qquad \boldsymbol{\kappa}_t =\left(\! \begin{array}{c} \kappa_t^{(1)} \\ \vdots\\ \kappa_t^{(N)} \end{array} \!\right), \qquad \boldsymbol{\xi}_t^\kappa\sim N(\mathbf{0},\Sigma),\end{equation}

where $\boldsymbol{\delta}=\{\delta^{(1)},\ldots, \delta^{(N)}\}'$ is an N-dimension column vector of drift parameters and $\Sigma$ is the $N\times N$ variance–covariance matrix of the multivariate white noise $\boldsymbol{\xi}_t^\kappa$ . For the cohort effect, the standard approach is to assume that it follows a univariate ARIMA model:

(3.6) \begin{equation} {{\Delta^{d} \gamma_{c} = \delta_{0} + \sum_{k=1}^{p} \phi_{k} \Delta^{d} \gamma_{c-k} + \xi_{c} + \sum_{j=1}^{q} \delta_{j} \xi_{c-j}.}}\end{equation}

where $\Delta$ is the difference operator, $\phi_{k}$ are the autoregressive coefficients, $\delta_{j}$ are the moving average coefficients, and $\xi_{c}$ is a Gaussian white noise process with variance $\sigma_{\xi}^2$ .

Let $t_{n_y}$ denote the last year where mortality data are available. To obtain mortality projections for forecast horizon h, one first obtains period index forecasts $\hat{\boldsymbol{\kappa}}_{t_{n_y}+h}$ and, if necessary, a cohort index forecast $\hat{\gamma}_{t_y-x+h}$ using their respective time series models. The mortality projections at horizon h are then given by:

(3.7) \begin{equation}\skew{22}\widehat{\ln(\mu_{x,t_{n_y}+h})} = \hat{\alpha}_{x} + \sum_{i=1}^{N} \hat{\beta}_{x}^{(i)} \hat{\kappa}_{t_{n_y}+h}^{(i)} + \hat{\gamma}_{t_{n_y}-x+h}\end{equation}

or

(3.8) \begin{equation}\skew{30}\widehat{\mathop{\mathrm{logit}}(q_{x,t_{n_y}+h})} = \hat{\alpha}_{x} + \sum_{i=1}^{N} \hat{\beta}_{x}^{(i)} \hat{\kappa}_{t_{n_y}+h}^{(i)} + \hat{\gamma}_{t_{n_y}-x+h}.\end{equation}

4.1. Model framework

Whilst in the GAPC framework the number of deaths is assumed to be drawn from either a Poisson or Binomial distribution, we model the dependent variable, $\ln(D_{x,t}/E^c_{x,t})$ , as a Gaussian generalised linear model. This decision is made for computational efficiency, noting that closed-form solutions for group regularisation estimates exist when using ordinary least squares (Breheny and Huang Reference Breheny and Huang2015).

The main drawback of ordinary least squares estimation is that the errors are assumed to be homoscedastic. As noted by Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002a), this is unrealistic as the logarithm of the observed force of mortality becomes much more variable at older ages due to the significantly smaller absolute number of deaths, which may have important consequences for inference and interval estimates of the model parameters (Brouhns et al. Reference Brouhns, Denuit and Vermunt2002b). Fortunately, for central parameter estimates and central predictions, the differences between the homoscedastic log-Gaussian approach and the heteroscedastic Poisson approach tend to be minimal (Brouhns et al. Reference Brouhns, Denuit and Vermunt2002a; Renshaw and Haberman Reference Renshaw and Haberman2003a).

Therefore, in the interest of computational efficiency, our general base model consists of the following modified version of the GAPC framework described in Equations (3.1)–(3.4):

1. 1. The random component: the logarithm of the force of mortality follows a Normal distribution so that

   (4.1) \begin{equation}\ln(D_{x,t}/E^c_{x,t}) \sim \text{Normal}(\!\ln (\mu_{x,t}), \sigma^2).\end{equation}

2. 2. The link function: the identity link (canonical) is used, such that

   (4.2) \begin{equation}\eta_{x,t} = g(E[\!\ln(D_{x,t}/E^c_{x,t})]) = \ln(\mu_{x,t}).\end{equation}

3. 3. The systematic component: the predictor captures age, period, and cohort effects, using

   (4.3) \begin{equation}\eta_{x,t} = \alpha_x + \sum_{i=1}^B f^{(i)}(x) \kappa_t^{(i)} + \gamma_{t-x},\end{equation}
   where the number of age-period terms, B, is very large. Here,

4. − $\boldsymbol{\alpha} = (\alpha_1,\ldots,\alpha_{n_a})$ is a row vector of coefficients to capture patterns by age;

5. − $\boldsymbol{\kappa^{\mathbf{(i)}}} = (\kappa_1^{(i)},\ldots,\kappa_{n_y}^{(i)})$ are row vectors of coefficients to capture systematic improvements over time;

6. − $f^{(i)}(x)$ are pre-defined functions to modulate trends by age; and

7. − $\boldsymbol{\gamma} = (\gamma_1,\ldots,\gamma_{n_b})$ is a row vector of coefficients to capture fixed effects embedded in years of birth.

Crucially, we include a large number, B, of age-period terms in the initial model specification, which allows the model to select from a wide range of potential trends. ^{Footnote 1} This is done through the specification of the different factor loadings, $f^{(i)}(x)$ , which modulate how the random period factors, $\kappa^{(i)}_t$ , impact different ages to determine how the age-period mortality curve changes through time. As a starting point, we consider five different continuous functional forms for $f^{(i)}(x)$ , which are summarised in Table 2 and illustrated in Figure 1(a). Within each functional form, we also consider various parametrisations, as illustrated in Figure 1(b).

Figure 1.

Examples of age-modulating basis functions, $f^{(i)}(x)$ , included in initial group regularisation model specification (a) by type, and (b) by parametrisation (for call and put option pay-offs).

Table 2.

Basis functions, $f^{(i)}(x)$ , included in initial group regularisation model specification.

The unit period term determines the general level of mortality and is included to capture systematic improvements over time that affect all ages equally, arising from increasing social advances (Cairns et al. Reference Cairns, Blake and Dowd2006). The linear or slope term, however, could capture the recent trend towards an increasing concentration of deaths around the mean age of death, known as “rectangularisation” (Biffis Reference Biffis2005). Further, the inclusion of higher-order polynomials is inspired by the curvature in $\mathop{\mathrm{logit}}(q_{x,t})$ plots, mirroring various extensions to the CBD model (Li and O’Hare Reference Li and O’Hare2017; Cairns et al. Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2010). Similarly, the call and put option pay-offs are motivated by the cPLAT model (Plat Reference Plat2009); however, we consider a greater range of strikes. We note that the framework and its associated implementation are flexible enough to allow the consideration of other functional forms beyond those in Table 1. For example, we could also consider the parametric forms in the “Age Function” toolkit of Hunt and Blake (Reference Hunt and Blake2014), orthogonal piecewise linear splines as in Aro and Pennanen (Reference Aro and Pennanen2011), Wavelets as in Hainaut and Denuit (Reference Hainaut and Denuit2020), or B-Splines (Currie et al. Reference Currie, Durban and Eilers2004). We explore the use of this latter type of basis function in our empirical exercise in Sections 5 and 6.

4.2. Model identifiability

As with most GAPC models, parameter identifiability might be a concern, especially given the diversity of basis functions allowed within the initial model (see Hunt and Blake, Reference Hunt and Blake2020a,b and Currie Reference Currie2020 for a detailed discussion on this matter). For example, we can freely move the general level of mortality across the age, period, and cohort parameters through the transformations:

(4.4) \begin{equation}\left(\alpha_x, \kappa^{(i)}_t\right)\to\left(\alpha_x + b_if^{(i)}(x), \kappa^{(1)}_t - b_i\right), \quad \text{$i =1,\ldots, B$},\end{equation}

and

(4.5) \begin{equation}\left(\alpha_x, \gamma_{t-x}\right)\to\left(\alpha_x + b, \gamma_{t-x} - b\right)\!,\end{equation}

for b and $b_i$ , $i =1,\ldots, B$ , arbitrary constants. Thus, we set $\kappa_{1}^{(i)}=0$ , $i=1,\ldots, B$ so that $\alpha_x$ can be interpreted as the age shape of mortality in the first year of the data. Similarly, we set $\gamma_{1}=0$ .

Besides the shifting of mortality levels, Hunt and Blake (Reference Hunt and Blake2020a) show that if a GAPC model contains basis functions spanning polynomials up to order $M-1$ , then the model possesses invariant transformations which allow the arbitrary movement of polynomial trends of order M across the age, period, and cohort parameters. However, we do not impose any explicit parameter constrains to deal with such invariant transformations, noting that we mitigate the possible impact of this identifiability issue on projections by modelling both the period and cohort indices using time series models with an order of integration $d=1$ . As discussed in Currie (Reference Currie2020), such a choice of the order of integration implies that in most situations, the forecast mortality rate values would be invariant despite the possible non-identifiability of the model parameters. Moreover, the empirical evidence suggests that mortality rates across most ages and many countries have unit roots and are difference stationary, justifying the use of $d=1$ (Njenga and Sherris Reference Njenga and Sherris2011). ^{Footnote 2}

4.3. Generalised linear model representation

Leveraging the result from Currie (Reference Currie2016), we can express the systematic component as a generalised linear model. Let $\boldsymbol{\eta} = \text{vec}(\!\ln(\boldsymbol{\mu}))$ ; here, the vec operator stacks the columns of a matrix in column order on top of each other, so $\boldsymbol{\eta}$ is column vector of dimensions $n\times1$ with $n=n_an_y$ . We can then write Equation (4.3) in matrix form:

(4.6) \begin{equation}\boldsymbol{\eta} = \boldsymbol{X} \boldsymbol{\beta} = \sum_{j=0}^{B+1} \boldsymbol{X}_j \boldsymbol{\beta}_j,\end{equation}

where the coefficients vector, $\boldsymbol{\beta}$ , can be decomposed into $B+2$ non-overlapping groups:

(4.7) \begin{equation}\boldsymbol{\beta} = (\boldsymbol{\beta_0}', \text{ } \boldsymbol{\beta_1}', \text{ } \boldsymbol{\beta_2}', \text{ } \ldots \text{ } , \text{ } \boldsymbol{\beta_B}', \text{ } \boldsymbol{\beta_{B+1}}')' = (\boldsymbol{\alpha}, \text{ } \boldsymbol{\kappa^{(1)}}, \text{ } \boldsymbol{\kappa^{(2)}}, \text{ } \ldots \text{ } , \text{ } \boldsymbol{\kappa^{(B)}}, \text{ } \boldsymbol{\gamma})'.\end{equation}

The first group refers to the vector of static age coefficients, $\boldsymbol{\alpha}$ . The next B groups denote each set period of period coefficients, $\boldsymbol{\kappa^{(i)}}$ , $i=1,\ldots, B$ , while the final group denotes the vector of cohort coefficients, $\boldsymbol{\gamma}$ . The design matrix, $\boldsymbol{X}$ , can also be divided:

(4.8) \begin{equation}\boldsymbol{X} = [\boldsymbol{X_0}:\boldsymbol{X_1}:\boldsymbol{X_2}: \ldots : \boldsymbol{X_B}: \boldsymbol{X_{B+1}}],\end{equation}

where $\boldsymbol{X}_j$ refers to the portion of the design matrix related to the predictors in Group j. For example, the first portion, $\boldsymbol{X_0}$ $(n \times n_a)$ , defines the design matrix corresponding to $\boldsymbol{\alpha}$ . Similarly, $\boldsymbol{X_i}$ $(n \times n_y)$ defines each $\boldsymbol{\kappa^{(i)}}$ for $i = 1,2,...,B$ , and $\boldsymbol{X_{B+1}}$ $(n \times n_b)$ defines $\boldsymbol{\gamma}$ .

4.4. Group regularisation fitting

Typically, to estimate parameters in GAPC-style models, $\boldsymbol{\hat\beta} = (\boldsymbol{\hat\alpha}, \text{ } \boldsymbol{\hat\kappa^{(1)}},$ $ \text{ } \boldsymbol{\hat\kappa^{(2)}}, \text{ } \ldots \text{ }, \text{ } \boldsymbol{\hat\kappa^{(B)}}, \text{ } \boldsymbol{\hat\gamma})'$ , we maximise the log-likelihood or minimise the residual sum of squares. These will estimate all the coefficients in the model, noting that all values will likely be non-zero. Whilst suitable for existing pre-defined GAPC models, these approaches are not appropriate for our framework. If every estimated coefficient from our model were nonzero, we would have a highly over-parametrised and complex model with $n_a + B\times n_y + n_b$ parameters. This would be problematic due to interpretability and predictive accuracy, commonly referred to as the “ $p>n$ ” problem (James et al. Reference James, Witen, Hastie and Tibshirani2013).

Instead, we estimate parameters using group regularisation (Yuan and Lin Reference Yuan and Lin2006; Breheny and Huang Reference Breheny and Huang2015), which minimises the objective function:

(4.9) \begin{equation}Q(\boldsymbol{\beta}|\boldsymbol{X}, \boldsymbol{\eta}) = L(\boldsymbol{\beta}|\boldsymbol{X}, \boldsymbol{\eta}) + \sum_{j=1}^{B+1} p_{\lambda}(\lVert \boldsymbol{\beta_j} \rVert),\end{equation}

where $L(\boldsymbol{\beta}|\boldsymbol{X}, \boldsymbol{\eta})$ denotes the residual sum of squares, and $p_{\lambda}(\cdot)$ is a penalty function applied to the Euclidean ( $\ell_2$ ) norm of each group to prevent the model from becoming too complex.

The key result is that the penalty function carries out group selection, meaning that within a group, coefficients will either all be zero or all be non-zero: if Group j is selected, then $\hat\beta_{jk} \ne 0$ for all k, and otherwise $\hat\beta_{jk} = 0$ for all k. For example, if Group 1 is not selected, $\hat\kappa_t^{(1)}$ will be zero for all t, and the age-period pair $f^{(1)}(x)\kappa_t^{(1)}$ will be eliminated from the model. Therefore, this estimation approach will select basis functions from the initial large model to best capture the period and cohort patterns in each data set.

There are alternative ways of defining the penalty function, $p_{\lambda}(\cdot)$ , including the group lasso penalty (Yuan and Lin Reference Yuan and Lin2006), the smoothly clipped absolute deviation penalty (Fan and Li Reference Fan and Li2001), and the minimax concave penalty (Zhang Reference Zhang2010). We have decided to use the minimax concave penalty (MCP),

(4.10) \begin{equation}p_{\lambda, \gamma}(\theta) =\begin{cases}\lambda \theta - \frac{\theta^2}{2\nu}, & \text{if } \theta \leq \nu \lambda, \\[4pt]\frac{1}{2} \nu \lambda^2, & \text{if } \theta > \nu \lambda,\end{cases}\end{equation}

which is a function of the tuning parameter, $\nu$ , and the penalty term, $\lambda$ . The MCP not only selects which groups of basis function to include but also shrinks estimated coefficients to prevent the model from over-fitting to noise, and thus each unique $(\nu, \lambda)$ defines a unique fitted model. However, for ease of computation, we fix the tuning parameter to be $\nu = 3$ , which aligns with common practice (Breheny and Huang Reference Breheny and Huang2015). Thus, $\lambda$ solely controls the trade-off between goodness of fit and parsimony; the larger the value of $\lambda$ , the heavier the penalisation, resulting in a more parsimonious model (fewer included basis functions and smaller estimated coefficients).

It is worth noting that the objective function in Equation (4.9) depends on the original scale of the design matrix $\boldsymbol{X}$ , which in turn depends on the original scale of the basis functions $f^{(i)}(\cdot)$ , $i=1,\ldots,B$ . To abstract from the original scale of the basis functions, the group-specific design matrices $\boldsymbol{X_j}$ , $j=1,\ldots,B+1$ , are orthonormalised prior to fitting and then transformed back to the original scale after fitting the model (Breheny and Huang Reference Breheny and Huang2015).

4.5. Block cross-validation for $\lambda$ selection

The value of $\lambda$ in the group regularisation defines a regularisation path containing GAPC mortality models with varying levels of complexity among which we need to choose the model considered best for a given application. To determine the optimal degree of penalisation, we propose a block cross-validation-style approach (Racine Reference Racine2000; Bergmeir and Benitez Reference Bergmeir and Bentez2012; Bergmeir et al. Reference Bergmeir, Costantini and Bentez2014) that aims to assess how well a model is likely to perform over a given forecasting horizon. This approach, which is in the spirit of the resampling methods discussed in Atance and DebÓn (Reference Atance and DebÓn2020), involves iteratively dividing the in-sample data into training data to fit the models and test data to assess forecasting abilities.

First, a portion of the data (roughly one-third) at the end is held back for validation, denoted as $\mathcal{V}$ and shown by the green points in Figure 2. This corresponds to time periods $t = n_T+1, n_T+2, \ldots , n_T+n_E$ . The remaining data ( $t = 1, 2, \ldots , n_T$ ) are used for training and testing, where the test sets can take different widths to reflect differing forecasting horizons, as illustrated in Figure 2. For example, to select a model to forecast a 1-year horizon, test sets are defined as 1-year blocks (top row of Figure 2), whereas by contrast, to select a model to forecast a 3-year horizon, test sets are defined as 3-year blocks (bottom row of Figure 2). Therefore, the training sets (blue) for forecasting horizon h are defined as:

(4.11) \begin{equation}\mathcal{T}^k(h) = \{(x,t)\}_{t \notin [k+1, k+h]},\end{equation}

Figure 2.

Iterations of cross-validation for horizon $h=1$ (top row) and $h=3$ (bottom row). for all $k = 1, 2, ..., n_T-h-2,n_T-h-1$ .

and the test sets (red) as:

(4.12) \begin{equation}\mathcal{F}^k(h) = \{(x,t)\}_{t = k+h},\end{equation}

Inspired by Bergmeir et al. (Reference Bergmeir, Hyndman and Koo2018), we use data both to the left and right of each test set to train the model, which contrasts heavily to the more common rolling window out-of-sample evaluation illustrated in Figure 3 (Tashman Reference Tashman2000; Dowd et al. Reference Dowd, Cairns, Blake, Coughlan, Epstein and Khalaf-Allah2009). This allows us to conduct significantly more test sets for a given data set than if we were restricted to fitting models using only data to the left, since this would place a limit on how far to the left test sets can occur. Moreover, we also minimise the portion of data that is left completely unused (grey), thereby extracting as much information as possible and producing an error estimate that is less variable than that of the traditional out-of-sample procedure (Bergmeir et al. Reference Bergmeir, Costantini and Bentez2014).

Figure 3.

Iterations of a rolling window approach for horizon $h=3$ .

Given that every training set will contain at least one observation from every age (row) and cohort (major diagonal) of the data set, the full vector of coefficients will be produced for the age, $\boldsymbol{\hat{\alpha}}$ , and cohort, $\boldsymbol{\hat{\gamma}}$ , terms. In other words, we will be able to estimate a coefficient corresponding to every age group and to every cohort. However, there will not be any observations from the non-training years, and thus all $\boldsymbol{\hat{\kappa}^{(i)}}$ will not be estimated fully. Figure 4 illustrates three examples, where the blue circles correspond to $\hat\kappa_t^{(i)}$ values that can be estimated directly from the training data, and the grey shaded regions to non-training years. Therefore, in order to assess the ability of a model to predict the test sets, we need to impute the missing values. This is achieved by fitting a random walk with drift using all the period indices that are observable (blue circles), both to the left and right of the missing region:

(4.13) \begin{equation}\kappa_t^{(i)} = \kappa_{t-1}^{(i)} + \delta^{(i)} + \epsilon_t^{(i)},\end{equation}

where $\epsilon_t^{(i)}$ is a Gaussian white noise process and ${\delta}^{(i)}$ is the drift. We can then predict the red circles using a forward fill procedure, basing predictions on the value from the previous period:

(4.14) \begin{equation}\hat{\kappa}_t^{(i)} = \hat{\kappa}_{t-1}^{(i)} + \hat{\delta}^{(i)}.\end{equation}

With these estimates, we can then compute the fitted/predicted values:

(4.15) \begin{equation}\widehat{\ln(\mu_{x,t})} = \hat{\alpha}_x + \sum_{i=1}^B f^{(i)}(x) \hat{\kappa}_t^{(i)} + \hat{\gamma}_{t-x},\end{equation}

and the mean squared error (MSE) for the test set:

(4.16) \begin{equation}\text{CV}_k(h) = \frac{1}{n_a} \sum_{(x,t) \in \mathcal{F}^k(h)} \left(\widehat{\ln\left(\mu_{x,t}\right)}-\ln\left(\mu_{x,t}\right)\right)^2.\end{equation}

This procedure is iterated for all test sets, $\mathcal{F}^k(h)$ , $k = 1, 2, ..., n_T-h-2,n_T-h-1$ , to determine the cross-validation error at horizon h:

(4.17) \begin{equation}\text{CV}(h) = \frac{1}{n_T-h-1} \sum_{k=1}^{n_T-h-1} \text{CV}_k(h).\end{equation}

This approach can be used to compare between a number of models and to select the “optimal” model that has the lowest cross-validation error for forecasting horizon h. This includes being able to test a grid of penalties, $\boldsymbol{\Lambda} = \{\lambda_\ell\}_{\ell=1}^{L}$ , in our group regularisation approach to select, for every horizon $h=1,\ldots,H$ , the value which minimises this error:

(4.18) \begin{equation}\hat{\lambda}(h) = \mathop{\mathrm{argmin}}_{\lambda \in \boldsymbol{\Lambda}} \text{CV}(h;\,\lambda).\end{equation}

Figure 4.

Illustration of time series imputation approach. The blue points correspond to $\hat\kappa_t$ estimates that can be estimated directly from the training set, and the grey area to the non-training years that cannot be estimated with the red points indicating the values imputed using a forward fill procedure.

4.6. Rolling window model evaluation

After selection, we must then evaluate our chosen model using completely independent and unseen data to ensure that it performs competitively relative to its peers. We implement out-of-sample evaluations with a rolling window approach, akin to Dowd et al. (Reference Dowd, Cairns, Blake, Coughlan, Epstein and Khalaf-Allah2009), using the validation data, $\mathcal{V}$ .

As shown in Figure 5, we fit the models using all in-sample data (blue) and forecast on the unseen validation set (red), which we assume contains $n_E$ years. We then roll forward one period $n_E-h$ times until all data are exhausted. Similar to the cross-validation procedure, forecasting is enabled by extrapolating the period and cohort parameters over the required time horizon, with a difference being that this is now achieved using conventional time series forecasting as discussed in Section 3.1. The out-of-sample MSEs at horizon h for each validation set can then be computed, using

(4.19) \begin{equation}\text{MSE}_j(h) = \frac{1}{n_a} \sum_{(x,t) \in \mathcal{V}^j(h)} \left(\widehat{\ln(\mu_{x,t})}-\ln(\mu_{x,t})\right)^2,\end{equation}

where $\mathcal{V}^j(h)= \{(x,t)\}_{t = n_T+j+h-1}$ for $j=1,\ldots,n_E-h+1$ . The overall out-of-sample MSE at horizon h is then given by:

(4.20) \begin{equation}\text{MSE}(h) = \frac{1}{n_E-h+1} \sum_{j=1}^{n_E-h+1} \text{MSE}_j(h).\end{equation}

Figure 5.

Iterations of rolling window out-of-sample evaluation for horizon $h=1$ (top row) and $h=3$ (bottom row).

In the cross-validation and evaluation of models, we follow Booth et al. (Reference Booth, Hyndman, Tickle and de Jong2006) and focus on the log-rates when assessing model accuracy. This is to reduce the impact of older ages which would otherwise dominate model selection and comparison due to the exponential increase in mortality rates with age. Moreover, without the log-transformation, older ages can distort model comparison due to their higher mortality uncertainty. However, it is noted that Equations (4.16) and (4.19) can readily be adapted to focus on errors on the untransformed mortality rates or errors in life expectancies.

5.1. Initial basis functions

To implement the GR framework, we consider a base model with $B=37$ basis functions taken from the functional forms in Table 1 as follows:

(5.1) \begin{align}f^{(1)}(x) &= 1; \end{align}

(5.2) \begin{align}f^{(2)}(x) &= x - \bar{x}; \end{align}

(5.3) \begin{align}f^{(2+j)}(x) &= (x - \bar{x})^{j+1}, \quad j=1,\dots,9; \end{align}

(5.4) \begin{align}f^{(11+r)}(x) &= (x - 20 - 5r)^+, \quad r=1,\dots,13; \end{align}

(5.5) \begin{align}f^{(24+p)}(x) &= (20 + 5p - x)^+, \quad p=1,\dots,13. \end{align}

5.2. Regularisation path

To search for the optimal level of regularisation, we consider a grid of 25 values of $\lambda$ equally spaced on the log-scale and in the range $\left[e^{-9},e^{-3.5}\right]$ . This range is wide enough to produce a regularisation path which starts with a very over-parametrised model and ends in a very parsimonious model. The resulting regularisation path is shown in Figure 6, which shows the relative size of each Group j of coefficients $\lVert \boldsymbol{\mathit{{X}_{j}}} \boldsymbol{\beta}_j\rVert/\lVert\boldsymbol{\mathit{X}} \boldsymbol{\beta}\rVert$ as a function of $\ln(\lambda)$ . In this particular exercise, we have decided not to penalise the $\boldsymbol{\alpha}$ parameters so that the model always contains a static age term $\alpha_x$ , independent of the amount of regularisation. This is because we are considering the age range 20 to 89 years, where the static age term is always believed to be important to appropriately capture the accident hump at younger ages (Plat Reference Plat2009; Hunt and Blake 2021 b).

Figure 6.

Regularisation path of mortality models for US males fitted to age range of 20–89 years for years 1960–2000.

Small values of $\lambda$ correspond to low levels of regularisation and more complex models. As such, to the left of Figure 5 when $\lambda=0.0001234098$ , we start with the model:

(5.6) \begin{equation}\begin{split}\ln(\mu_{x,t}) = \alpha_x + \kappa_t^{(1)} + (\bar{x}-x) \kappa_t^{(2)} + (x-70)^+ \kappa_t^{(21)} + (x-85)^+ \kappa_t^{(24)} + \\ (25-x)^+ \kappa_t^{(25)} + (30-x)^+ \kappa_t^{(26)} +(45-x)^+ \kappa_t^{(29)} + \gamma_{c},\end{split}\end{equation}

which, in addition to the default static age term and the cohort effect, includes 7 of the 37 possible age-period terms. Noticeably, this very complex model includes terms that capture the idiosyncrasies of mortality trends and noise at very young and very old ages via the $\kappa_t^{(21)}$ , $\kappa_t^{(24)}$ , $\kappa_t^{(25)}$ , and $\kappa_t^{(26)}$ terms. As we increase the value of $\lambda$ , the additional regularisation induces the elimination of some age-period terms, resulting in more parsimonious models. For example, when we reach a value of $\lambda=0.0024276457$ , we obtain the model:

(5.7) \begin{equation}\ln(\mu_{x,t}) = \alpha_x + \kappa_t^{(1)} + (45-x)^+ \kappa_t^{(29)} + \gamma_{c},\end{equation}

which, as opposed to the seven age-period terms of the initial model, only includes two age-period terms capturing the general level of mortality and the differential improvement of individuals younger than 45 years of age. We note that this model is very similar to the sPLAT model, albeit without the linear term and with the strike of the put function at 45 rather than $\bar{x}=54.5$ .

5.3. Selected values of $\lambda$

Given the regularisation path, the next step is to select a model suitable for each forecasting horizon using the block cross-validation approach described in Subsection 4.5. The selected $\hat\lambda(h)$ for each horizon h is the value that minimises the cross-validation MSE. Figure 7 illustrates this for horizons $h=1, 5,$ and 10, where the selected values correspond to $\hat\lambda(1) = 0.0001234098$ , $\hat\lambda(5) = 0.001220697$ , and $\hat\lambda(10) = 0.0024276457$ , respectively. Note, for $\lambda$ values to the right of each respective $\hat\lambda$ , we have under-fitting, whilst to the left, we have over-fitting. Table 3 presents the selected $\hat\lambda(h)$ for all horizons of interest along with a summary of the corresponding model structure.

Figure 7.

Cross-validation MSE for forecasting horizons $h=1$ , $h=5$ , and $h=10$ for US males data fitted to ages 20–89 years for years 1960–2000.

Table 3.

Summary of optimal model structures for group regurlarisation for ages 20–89 years in the period 1960–2000 in US males.

In the case of 1-year and 2-year forecasting horizons, the selected model has seven age-period terms and corresponds to the specification in Equation (5.6). This contrasts with the 7 to 10-year forecasting horizons for which the selected model corresponds to the simpler structure in Equation (5.7), which has only two age-period terms. In general, as indicated by Table 3, the longer the forecasting horizon, the higher the value of $\hat\lambda(h)$ , the fewer the number of period terms, and the more parsimonious the selected model specification.

5.4. Out-of-sample evaluation

Having selected appropriate models for each horizon, we now evaluate their out-of-sample performance relative to existing GAPC models using the rolling window approach described in Section 4.6, with validation data from 2001 to 2017. Thus, for horizon $h=1$ , there are 17 years for evaluation, which reduces to 8 years for horizon $h=10$ namely, 2010, 2011, $\dots$ , 2017, whose forecasts are derived using the years 2000, 2001, $\dots$ , 2007, as the last training year, respectively.

In the evaluation exercise, we consider two versions of the GR framework obtained by varying the number and form of the initial basis functions. The first version, herein referred to as GR1, corresponds to the implementation described previously in Subsection 5.1 which uses the $B=37$ basis functions in Equations (5.1)–(5.5). The second version, referred to as GR2, explores the use of a larger pool of initial basis functions in order to enrich the types of trends that the model can capture. To do so, in addition to the 37 basis functions in Equations (5.1)–(5.5), we consider additional cubic B-Splines basis functions (Currie et al. Reference Currie, Durban and Eilers2004), which are piecewise polynomials commonly used in the graduation and forecasting of mortality (Macdonald et al. Reference Macdonald, Richards and Currie2018). Figure 8 depicts the basis functions we add, which are comprised of four groups of B-splines, each of which contain 2, 4, 8, and 16 functions, respectively. Thus, the total number of basis in the initial model under the GR2 implementation is $B=37+2+4+8+16=67$ .

Figure 8.

Cubic B-splines basis functions added to the initial model for the GR2 version of the group regularisation implementation.

We consider as benchmarks eight of the ten GAPC models presented in Table 1. We have discarded the CBD and M7 models as they are both designed to fit higher ages (typically 50 years of age and beyond), and thus are not entirely appropriate for the 20–89 years age range considered in our experiments. While we estimate our regularised models using ordinary least squares under a log-Gaussian setting, we estimate all the classical GAPC models using maximum likelihood under a log-Poisson setting. This conforms with the existing literature that favours the use of a generalised (non-)linear framework for estimating GAPC models (Hunt and Blake Reference Hunt and Blake2021b). Appendix A includes the constraints used to ensure the identifiability of these benchmark models.

For the forecasting of the period functions in GAPC and GR models, we follow the standard approach and use a multivariate random walk with drift (Cairns et al. Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2010; Haberman and Renshaw Reference Haberman and Renshaw2011; Villegas et al. Reference Villegas, Kaishev and Millossovich2018). In addition, for models with a cohort effect, we forecast $\gamma_{t-x}$ using an ARIMA model as in Equation (3.6) with level of differentiation $d=1$ , autoregressive order $p=1$ , and moving average order $q=0$ . As discussed before in Subsection 4, choosing $d=1$ tempers the possible impact of model identifiability on projections (Currie Reference Currie2020). Moreover, setting $p=1$ is in line with the level autoregression found by others when modelling cohort effects (Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2008); Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011); Hunt and Blake (2021 a)).

Table 4 presents the out-of-sample MSE by forecasting horizon for each of the eight GAPC models and the two implementations of the GR framework. Among the standard GAPC models, the sPLAT and cPLAT model have the best forecasting performance at every horizon, with the sPLAT model topping the ranking at 7 of the 10 forecasting horizons.

Table 4.

Out-of-Sample Mean Square Error for US males by forecasting horizon with model ranking in brackets. Values are presented as $1000 \times$ MSE.

The GR models perform competitively, outperforming the CBDx, LC, LC2, and M7x models at every horizon and the APC and RH at most horizons. Moreover, the GR approaches show a very consistent performance across all horizons, with MSE values which are very close to the MSE values of the sPLAT and cPLAT models. This is not surprising in the case of GR1, as for the US male data set for horizons $h=7, 8, 9, 10$ , the GR1 approach identified a model which closely resembles the cPLAT model (recall Equation (5.7)).

It is worth noting that for horizons 1 and 2, the GR2 approach outperforms all nine other competing models by a significant margin. By including the 67 additional basis functions depicted in Figure 8, the GR2 model can explore more complex models that are able to better capture the cross-sectional patterns in the mortality data. Hence, the improved performance of the GR2 over the rest of the models highlights the benefits of allowing for a larger variety of age-specific trends and complexity of model when conducting short-term mortality forecasts. However, this extra complexity comes at the cost of ease of interpretability of coefficients relative to the GR1 models.

6.1. Summary of model structures

The results of applying the GR approach reveal that the optimal model varies not only by country, but also within country, depending on the forecasting horizon. To illustrate this, Tables 5 and 6 report, respectively, the average characteristics of the model structures produced by the GR1 and GR2 approaches for different horizons.

Table 5.

Summary of optimal model structures according to GR1 for ages 20–89 years in the period 1960–2000 in 52 data sets.

Table 6.

Summary of optimal model structures according to GR2 for ages 20–89 years in the period 1960–2000 in 52 data sets.

Notably, across all 52 data sets, $\hat\lambda(h)$ tends to increase with the horizon of interest, which means that simpler models are constructed for longer forecasting periods. For example, while at horizon 1 the models have on average 4.77 and 6.98 period terms for GR1 and GR2, respectivley, at horizon 10, the average number of period terms reduces to 2.27 and 3.29. These results conform with the hypothesis of parsimony being preferred for longer horizons (Green and Armstrong 2015). More specifically, for shorter forecasting horizons, the predictive success of a model is determined by its ability to accurately fit the cross section of the mortality curve. By contrast, as the forecasting horizon increases, it becomes more important to filter out the year-on-year noise in the data, and having a closer fit to the initial mortality curve becomes less relevant.

It is also worth noting that the longer the horizon, the less likely the optimal structure will include a cohort effect. While at short horizons the optimal structure will almost certainly include a cohort effect, at horizon 10, only 60% of the 52 data sets will have a cohort effect under the GR1 approach. This proportion is further reduced to 35% under the GR2 approach, indicating that in the GR1 approach, the cohort effect might be acting as a substitute for uncaptured age-period trends.

As an illustration of the diversity of models produced by the GR framework, Table 7 shows the optimal model produced by the GR1 approach at forecasting horizon 10 for all 52 data sets. In some cases, the optimal models are similar in structure to existing GAPC models. For example, for US females, Japanese males, UK (GBR-NP) females, Belarusian males, Austrian males, and Bulgarian males, the structure of the optimal model is the same as the sPLAT model. Similarly, for Swedish males, Czech females, and Irish females, the optimal model structure coincides with the CBDx model. In most cases, though, the optimal models vary quite remarkably from existing models, but with the following commonalities in the type of basis functions used for the different countries:

• 79% of models contain a unit period term, $\kappa_t^{(1)}$ , which captures changes in the general level of mortality that affect all ages.

• 42% of models contain a linear term, $(x-\bar{x})\kappa_t^{(2)}$ which captures changes in the slope of the mortality curve.

• Only 6% of models contain a quadratic term, $(x-\bar{x})^2\kappa_t^{(3)}$ , which captures the curvature of the mortality curve.

• 37% of models contain at least one of the terms, $(x-25)^+\kappa_t^{(23)}$ , $(x-30)^+\kappa_t^{(24)}$ , or $(x-35)^+\kappa_t^{(25)}$ , which specifically target declining mortality rates at ages below 35 years, known as the “flattening of the accident hump.”

Table 7.

Optimal model structures according to GR1 for 10-year forecasting horizon using data for ages 20–89 years in the period 1960–2000.

Existing GAPC models are unable to capture the right combination of the above findings, stressing the merits of our data-driven approach and the pitfalls of using pre-defined structures for all data sets.

6.2. Out-of-sample evaluation

Whilst the previous section implies that a wide range of bespoke models may be more appropriate to explain historical data, it is important to assess whether the GR-constructed models generalise well out-of-sample and if they perform competitively relative to existing GAPC models. For the 52 data sets, we replicate the out-of-sample evaluation described in Section 5.4.

Figure 9 presents heatmaps reporting the ranking of the eight GAPC models and the two GR approaches for forecasting horizons $h=1, 5,$ and 10. In this figure, the GAPC models are ordered in terms of their complexity as indicated by their number of parameters. Similarly, the countries are ordered according to their population size. To facilitate the interpretation of these model rankings, Figure 10 presents the mean ranking of each of the 12 models as a function of the 10 forecasting horizons of interest.

Figure 9.

Ranking of models according to out-of-sample MSE at $h=1$ , $h=5$ , and $h=10$ , for different data sets.

Figure 10.

Average of out-of-sample MSE rank across all 52 data sets for each model.

Over the 1-year horizon, the GR2 approach produces the best result with a median rank of 2 and is in the top three models in 41 of the 52 data sets. Interestingly, the best performance of the GR2 approach is among the countries with the biggest populations. The GR1 approach also has a very strong performance at short horizons, surpassing most GAPC models, but is almost always outperformed by the GR2 approach. This reiterates the benefit of allowing for a larger variety of age-specific trends when focusing on short-term forecasting so that the cross-sectional variability of the mortality curve is better captured. The cPLAT and M7x models are the two best-performing GAPC models, particularly for larger populations, whilst the more parsimonious CBDx, APC, and LC models perform relatively poorly. This again reaffirms our hypothesis that complex models are preferred for shorter horizons and for less noisy data sets.

As we increase the forecasting horizon, the absolute values of the forecasting errors increase for all models, as would be expected. However, the relative rankings of the models also change, which is of far more interest. The strong performance of our GR approaches remain consistent, with the GR1 in particular having the best median rank at horizons 5 and 10, reflecting the ability of our regularisation approach to filter out the year-on-year noise in the mortality trends to produce a parsimonious model which concentrates on the secular longer-term trends. ^{Footnote 3} In contrast, the other models show a significant shift, with more parsimonious models, such as the APC and the CBDx, now performing significantly better than the complex cPLAT, M7x, and LC2 models. This notion of simplicity outperforming complexity over longer horizons is further supported by the improving performance of the GR1 models relative to GR2 as h increases.

Interestingly, the popular cPLAT and LC models rank within the bottom four models almost exclusively for the longer horizons. Notably, one of the only data sets where the cPLAT model performs reasonably well over a 10-year horizon is for the US populations, which is actually the data that Plat (Reference Plat2009) used to construct this model. This lends support to the notion that if we examine data and build a model based on the patterns observed within it, as Plat (Reference Plat2009) did for the USA, we can expect to yield strong forecasting results. However, we cannot just assume that a model that explains US mortality well will fit well to all other countries. Thus, we must repeat this procedure for each data set, which is effectively what our data-driven approach aims to do.

7.1. Effect of generator model

Figure 11 illustrates forecasting errors using data generated by a number of GAPC models, including models defined in the literature (APC, LC, and cPLAT), and others that have been randomly constructed using a suite of basis functions. These models have been created to reflect varying numbers of age-period terms, the inclusion/exclusion of a cohort term, and a range of basis functions with different forms and parametrisations. This includes a basis function not considered in the GR approach in the bottom right panel, namely a unit period term applying to ages 50 years and below only.

Figure 11.

Boxplot of mean squared errors from 100 simulated data sets, involving (a) rates generated from different GAPC models, (b) using a population size of 5 million, and (c) a forecasting horizon of 10 years.

As expected, in the first three panels, the model that generated each data set performs very well with regard to both median errors and the range of errors across the 100 simulations. However, the performance of models on data generated from another model is of far greater interest. Notably, no single GAPC model performs strongly across all scenarios, with the complex models, such as the cPLAT and M7x, tending to perform stronger on the more complex data sets, and more parsimonious models, such as the LC and APC, performing better on the simpler data sets. More strikingly, this study highlights the significance of the cohort term, where models without a cohort term perform extremely poorly where cohort effects are prevalent in the underlying data. This is most pronounced in the bottom left panel of Figure 11, where the LC and CBDx models generate unwieldy errors.

By contrast, the GR1 model performs strongly and consistently, with the lowest median error in six of the eight scenarios, and ranking within the top three models every time. This is hardly surprising, given the adaptability of its framework to select more period and/or cohort terms where more complex patterns exist in the underlying data.

7.2. Effect of population size

When population size decreases (or the level of noise increases), it is obvious that forecasting errors will increase in absolute terms, as illustrated in Figure 12. However, the changes in the relative performance of the models are more insightful. For very small populations, complex models have a tendency to overfit the Poisson noise in the data, resulting in a poor forecasting performance, with higher median forecasting errors and increased tendencies to produce significantly large errors or anomalies. In contrast, the more parsimonious LC and CBDx models perform stably for these smaller populations as their reduced number of parameters allows them to better filter the Poisson noise.

Figure 12.

Boxplot of mean squared errors from 100 simulated data sets, involving (a) rates generated from an APC model, (b) using varying population sizes, and (c) a forecasting horizon of 10 years.

The GR1 model, however, remains competitive under all scenarios, registering the lowest median error on the 500,000 and 5 million populations, and recording extremely similar median errors to the best models for the larger 50 and 500 million populations. Notably, in all scenarios, the GR1 approach tends to select a model structure that is similar to or exactly matches that of the APC model that generated the data, which explains the strikingly similar results of these two models for larger populations. Interestingly, however, the GR1 model significantly outperforms the APC model for the smaller populations, which owes to the ability of the regularisation approach to not only select basis functions but also to shrink parameter estimates to avoid over-fitting to noise in the training data.

7.3. Effect of forecasting horizon

Figure 13 illustrates errors for varying forecasting horizons, where for consistency, the number of years of training data was amended in each scenario such that the GR1 approach had 30 test sets per horizon. For the short-term 1-year horizon, all models perform strongly in absolute terms, with similar median errors observed for all models. However, even for this horizon, it is observed that there is a greater tendency for the complex models such as the cPLAT, M7x, and sPLAT to produce large, unwieldy errors, owing to the potential to overfit to noise given the large number of parameters. This poor and volatile performance of the complex models is even more pronounced as the length of the forecasting horizon increases, where there is a noticeable difference between the two simplest GAPC models – CBDx and LC – and the others. This mirrors the behaviour observed in Section 6. For all horizons, the GR1 approach performs strongly, which again owes to the ability to simplify the underlying structure of the model as the forecasting horizon increases.

Figure 13.

Boxplot of mean squared errors from 100 simulated data sets, involving (a) rates generated from an APC model, (b) using a population size of 5 million, and (c) varying forecasting horizons.