2.1 HS models

Spline fitting is an interpolation method under which segments of low-degree polynomials connected at chosen control points are fitted to describe the data shape. One widely used category known as the cubic HS obtains the fitted curve using cubic polynomials defined in Hermite form (Hermite, Reference Hermite1864). In more detail, one specifies the values and first derivatives of the start and end points for each segment, and the shape of the fitted curve is controlled by the Hermite basis functions (Marschner & Shirley, Reference Marschner and Shirley2018). In matrix form, the interpolation polynomial $p\!\left( u \right)$ over $u \in \left[ {0,1} \right]$ can be expressed as

(1) \begin{align} p\!\left( u \right) = \left( {{u^3},{u^2},u,1} \right)\left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 2 & -2 & 1 & 1\\[4pt]-3 & 3 & -2 & -1\\[4pt]0 & 0 & 1 & 0\\[4pt]1 & 0 & 0 & 0\\[4pt]\end{array}\right)\left( \begin{array}{c} {{p_0}} \\[4pt] {{p_1}} \\[4pt] {{\delta _0}} \\[4pt] {{\delta _1}} \\[4pt] \end{array} \right) = \left( \begin{array}{c} {2{u^3} - 3{u^2} + 1} \\[4pt] { - 2{u^3} + 3{u^2}} \\[4pt] {{u^3} - 2{u^2} + u} \\[4pt] {{u^3} - {u^2}} \end{array} \right)^{\prime} \left( \begin{array}{c} {{p_0}} \\[4pt] {{p_1}} \\[4pt] {{\delta _0}} \\[4pt] {{\delta _1}} \end{array} \right) = \boldsymbol{f}\!\left( \boldsymbol{u} \right)\!\boldsymbol{P}, \end{align}

where $\boldsymbol{P}$ is the $\left( {4 \times 1} \right)$ control matrix, ${p_0}$ and ${p_1}$ represent the start and end points of the fitted function, ${\delta _0}$ and ${\delta _1}$ are the corresponding tangents at $u = 0$ and $u = 1$ and $\boldsymbol{f}\!\left( \boldsymbol{u} \right)$ contains basis functions of the cubic HS. The resulting HS $p\!\left( u \right)$ is a linear combination of the four basis functions.

Richards (Reference Richards2020) proposed the HS model that smooths post-retirement mortality ${m_x}$ in a given year using four third-order polynomials in a Hermite form,

(2) \begin{align}{ln}\, {m_x} = \alpha {h_{00}}\!\left( {{x_k}} \right) + \omega {h_{01}}\!\left( {{x_k}} \right) + {s_0}{h_{10}}\!\left( {{x_k}} \right) + {s_1}{h_{11}}\!\left( {{x_k}} \right),\end{align}

where

(3) \begin{align}\left\{\begin{array}{l}{{h_{00}}\!\left( {{x_k}} \right) = \left( {1 + 2{x_k}} \right){{\left( {1 - {x_k}} \right)}^2}} \\[4pt] {{h_{01}}\!\left( {{x_k}} \right) = x_k^2\!\left( {3 - 2{x_k}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;} \\[4pt] {{h_{10}}\!\left( {{x_k}} \right) = {x_k}{{\left( {1 - {x_k}} \right)}^2}\;\;\;\;\;\;\;\;\;\;\;\;\;\;}, \\[4pt] {{h_{11}}\!\left( {{x_k}} \right) = x_k^2\!\left( {{x_k} - 1} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \end{array} \right.\end{align}

${x_k} = \left( {x - {x_0}} \right)/\left( {{x_1} - {x_0}} \right)$ for a predetermined age interval $\left[ {{x_0},{x_1}} \right]$ and ${x_k} \in \left[ {0,1} \right]$ . There are four coefficients $\alpha $ , $\omega $ , ${s_0}$ and ${s_1}$ for the Hermite basis functionsFootnote 2 ${h_{ij}}\!\left( {{x_k}} \right),\;i = 0,1,\;j = 0,1$ to be estimated. Thereinto, $\alpha $ and $\omega $ correspond to ${p_0}$ and ${p_1}$ in (1) and reflect the level of $\ln {m_x}$ at the two boundary ages ${x_0}$ and ${x_1}$ , respectivelyFootnote 3.

We denote respectively $\boldsymbol{D}$ , $\boldsymbol{E}$ and $\boldsymbol{M}$ as the $\left( {{n_a} \times {n_y}} \right)$ matrices of the death counts, central exposures and central death rates sorted by age and year, where the total number of ages (years) equals ${n_a}$ ( ${n_y}$ ). It is assumed that the number of deaths is a random variable following a Poisson distribution. Suppose $\boldsymbol{d}$ , $\boldsymbol{e}$ and $\boldsymbol{m}$ are the vectors with a length of ${n_a}$ in the above matrices corresponding to a given year. In a matrix form, one may write the HS model as

(4) \begin{align}ln\ \boldsymbol{m} = {\boldsymbol{x}_{\textbf{0}}}{\boldsymbol\beta} ,\end{align}

where $\boldsymbol{m}$ is a vector of $\ln {m_x}$ with a length of ${n_a} = {x_1} - {x_0} + 1$ , ${\boldsymbol{x}_{\textbf{0}}}$ refers to a ( ${n_a} \times 4$ ) matrix containing the four Hermite basis functions at each age, $\boldsymbol{\beta} = \left( {\alpha ,\omega ,{s_0},{s_1}} \right)^{\prime}$ is the coefficient vector.

A smooth path between the youngest and oldest death rates in the interval is specified by ${s_0}$ and ${s_1}$ which give the gradients at the two boundary ages. As discussed in the original paper of the HS method, the fitted mortality curve may be too flexible and may deviate from the typical shape when both slope parameters are employed. Accordingly, four versions of the HS model can be constructed by setting one or both gradient coefficients to zero. The four models, namely HS1 ( ${s_0} = {s_1} = 0$ ), HS2 ( ${s_1} = 0$ ), HS3 ( ${s_0} = 0$ ) and HS4 can be defined by adjusting the coefficient vector $\boldsymbol\beta $ :

\begin{align*}\boldsymbol\beta^{\prime} = \left\{ {\begin{array}{*{20}{c}}{HS1:\left( {\alpha ,\omega ,0,\;0} \right)}\\[3pt]{HS2:\left( {\alpha ,\omega ,{s_0},0} \right)}\\[3pt]{HS3:\left( {\alpha ,\omega ,0,{s_1}} \right)}\\[3pt]{HS4:\left( {\alpha ,\omega ,{s_0},{s_1}} \right)}\end{array}} \right..\end{align*}

As explained in section 1, the classical Gompertz law requires two parameters to model the decreasing mortality differentials over age and can potentially cause an unjustified crossover at advanced ages. On the other hand, the four HS models can ensure the convergence using only the parameter $\alpha $ . In the work of Richards (Reference Richards2020), the initial mortality level ${\alpha ^{\left( i \right)}}$ for population $i$ can be expressed as:

(5) \begin{align}{\alpha ^{\left( i \right)}} = {\alpha ^{\left( 0 \right)}} + \mathop \sum \nolimits_{j = 1}^J {\alpha ^{\left( {{r_j}} \right)}}{z^{\left( {i,j} \right)}},\end{align}

where ${z^{\left( {i,j} \right)}}$ is an indicator variable which equals 1 if population $i$ has risk factor $j$ ( $j = 1,2, \ldots ,J$ ) and 0 otherwise. When a population does not exhibit any of the $J$ risk characteristics, the intercept of the fitted mortality curve is equivalent to the baseline value ${\alpha ^{\left( 0 \right)}}$ . The estimated mortality level at the starting age increases by ${\alpha ^{\left( {{r_j}} \right)}}$ for each additional risk factor j. By way of illustration, suppose that mortality data with the same gender of two states in the same country only differ in one risk factor – state. Under the simplest model HS1, their mortality differentials can be incorporated using the individual starting values ${\alpha ^{(1)}} = {\alpha ^{\left( 0 \right)}}$ and ${\alpha ^{(2)}} = {\alpha ^{\left( 0 \right)}} + {\alpha ^{\left( {{r_1}} \right)}}$ and a common ending value $\omega $ . Under the compensation law of mortality, the two sets of state-level mortality rates (in a specific year) should converge gradually as age increases. As shown in Figure 1 that plots the basis functions of HSs, ${h_{00}}\!\left( {{x_k}} \right)$ is monotonically decreasing over the defined age range. Consequently, the gap between the fitted mortality curves of the two populations vanishes over age by default, without the need of adopting a “slope” parameter.

Figure 1

Hermite basis splines ${h_{ij}}\!\left( {{x_k}} \right)$ for ${x_k}$ ranging from 0 to 1, $i = 0,1$ and $j = 0,1$ .

In addition to the parsimonious modelling of mortality differentials, the potential crossover problem under the Gompertz model can be avoided due to the monotonicity of the Hermite function ${h_{00}}\!\left( {{x_k}} \right)$ . Specifically, the difference between the log central death rates of the two populations at age $x$ under the HS1 model is calculated as

(6) \begin{align}ln\ m_x^{(1)} - ln\ m_x^{(2)} & = \left[ {{\alpha ^{(1)}}{h_{00}}\!\left( {{x_k}} \right) + \omega {h_{01}}\!\left( {{x_k}} \right)} \right] - \left[ {{\alpha ^{(2)}}{h_{00}}\!\left( {{x_k}} \right) + \omega {h_{01}}\!\left( {{x_k}} \right)} \right]\nonumber\\[3pt]& = \left( {{\alpha ^{(1)}} - {\alpha ^{(2)}}} \right){h_{00}}\!\left( {{x_k}} \right), \end{align}

where ${x_k} = \left( {x - {x_0}} \right)/\left( {{x_1} - {x_0}} \right)$ . Note that the two populations are assumed to share a common death rate $\ln {m_{{x_1}}} = \omega $ at an advanced age ${x_1}$ Footnote 4. Since ${h_{00}}\!\left( {{x_k}} \right)$ is a decreasing function ranging between 0 and 1, the difference in equation (6) narrows down to 0 as $x$ reaches ${x_1}$ but never switches its sign over the given range $\left[ {{x_0},{x_1}} \right]$ . Nonetheless, one may argue that mortality differentials resulted from certain risk factors such as gender reduce at higher ages but do not vanish entirely (Tickle, Reference Tickle1997). That is, one population always has lower death rates than the other for the entire age range considered $\big($ assuming ${\alpha ^{(1)}} \gt {\alpha ^{(2)}}$ and $\left.{\omega ^{(1)}} \gt {\omega ^{(2)}}\right)$ . To arrive at this pattern, a population-specific ending value ${\omega ^{\left( i \right)}}$ may be employed, producing mortality differentials as follows:

(7) \begin{align}ln\ m_x^{(1)} - ln\ m_x^{(2)} & = \left[ {{\alpha ^{(1)}}{h_{00}}\!\left( {{x_k}} \right) + {\omega ^{(1)}}{h_{01}}\!\left( {{x_k}} \right)} \right] - \left[ {{\alpha ^{(2)}}{h_{00}}\!\left( {{x_k}} \right) + {\omega ^{(2)}}{h_{01}}\!\left( {{x_k}} \right)} \right]\nonumber\\[3pt]& = \left( {{\alpha ^{(1)}} - {\alpha ^{(2)}}} \right){h_{00}}\!\left( {{x_k}} \right) + \left( {{\omega ^{(1)}} - {\omega ^{(2)}}} \right){h_{01}}\!\left( {{x_k}} \right).\end{align}

As illustrated in Figure 1, ${h_{00}}\!\left( {{x_k}} \right)$ is monotonically decreasing, while ${h_{01}}\!\left( {{x_k}} \right)$ is an increasing function. Besides, they are symmetric about the line $h\!\left( {{x_k}} \right) = 0.5$ $\left({h_{00}}\!\left( {{x_k}} \right) + {h_{01}}\!\left( {{x_k}} \right) = 1\right)$ . It implies that when $\;\left( {{\alpha ^{(1)}} - {\alpha ^{(2)}}} \right) = \left( {{\omega ^{(1)}} - {\omega ^{(2)}}} \right)$ , the sum of the two weighted functions in (7) becomes a horizontal line. Since $\left( {{\alpha ^{(1)}} - {\alpha ^{(2)}}} \right)$ and $\left( {{\omega ^{(1)}} - {\omega ^{(2)}}} \right)$ refer to the differences in the log mortality rates between the two populations at the youngest and oldest ages, respectively, one would expect the former to be more pronounced than the latter due to the compensation law of mortality (Gavrilov & Gavrilova, Reference Gavrilov and Gavrilova2001). Accordingly, ${h_{00}}\!\left( {{x_k}} \right)$ contributes more to the weighted sum in equation (7), leading to declining but non-zero differentials. Again, the no-crossover property can be guaranteed, as the weighted sum of the two Hermite basis functions does not change in sign over the domain interval. The other three HS models also possess this property if one sets common gradient parameters ${s_0}$ and ${s_1}$ between related populationsFootnote 5.

The above HS models were initially proposed to model mortality over the age range rather than over time. Next, we show that the HS approach can be extended to model two-dimensional population mortality data. Similar to the widely adopted CBD model, the HS structure is fitted to the mortality data of each year. By doing so, the estimated values of each model parameter over different years form a time series, based on which the future death rates can be projected or simulated. The modified versions of the HS models are given below:

(8) \begin{align}ln\ \boldsymbol{M} = \boldsymbol{XB} = \left( {{\boldsymbol{I}_{{\boldsymbol{n}_{\boldsymbol{y}}}}} \otimes {\boldsymbol{X}_{\textbf{0}}}} \right)\boldsymbol{B},\end{align}

where $\boldsymbol{I}_{{\boldsymbol{n}_{\boldsymbol{y}}}}$ is an identity matrix of size ${n_y}$ , ${\boldsymbol{x}_{\textbf{0}}}$ has a dimension of $\left( {{n_a} \times 4} \right)$ and comprises the Hermite basis functions at each age, the design matrix $\boldsymbol{x}$ is the Kronecker product of $\boldsymbol{I}_{{\boldsymbol{n}_{\boldsymbol{y}}}}$ and ${\boldsymbol{x}_{\textbf{0}}}$ , the coefficient vector ${\textbf{B}} = \left( {\boldsymbol\beta^{\prime}_{\boldsymbol{t}_{\textbf{0}}},\boldsymbol\beta^{\prime}_{{\boldsymbol{t}_{\textbf{0}}} + \textbf{1}}, \ldots , \boldsymbol\beta^{\prime}_{{\boldsymbol{t}_{\textbf{1}}}}} \right)^{\prime}$ concatenates the year-specific HS coefficients ${\boldsymbol\beta _\boldsymbol{t}}$ . Again, the four HS versions can be obtained by adjusting ${\boldsymbol\beta _\boldsymbol{t}}$ . The modified versions of the HS models are given below:

\begin{align*}\left\{ {\begin{array}{l}{HS1:\;\;\;ln\ {m_{x,t}} = {\alpha _t}{h_{00}}\!\left( {{x_k}} \right) + {\omega _t}{h_{01}}\!\left( {{x_k}} \right)}\\[4pt]{HS2:\;\;\;ln\ {m_{x,t}} = {\alpha _t}{h_{00}}\!\left( {{x_k}} \right) + {\omega _t}{h_{01}}\!\left( {{x_k}} \right) + {s_{0,t}}{h_{10}}\!\left( {{x_k}} \right)}\\[4pt]{HS3:\;\;\;ln\ {m_{x,t}} = {\alpha _t}{h_{00}}\!\left( {{x_k}} \right) + {\omega _t}{h_{01}}\!\left( {{x_k}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; + {s_{1,t}}{h_{11}}\!\left( {{x_k}} \right)}\\[4pt]{HS4:\;\;\;ln\ {m_{x,t}} = {\alpha _t}{h_{00}}\!\left( {{x_k}} \right) + {\omega _t}{h_{01}}\!\left( {{x_k}} \right) + {s_{0,t}}{h_{10}}\!\left( {{x_k}} \right) + {s_{1,t}}{h_{11}}\!\left( {{x_k}} \right)}\end{array}} \right..\end{align*}

Similar to the one-dimensional case, ${\alpha _t}$ and ${\omega _t}$ can be interpreted as the fitted values of $\ln {m_{{x_0},t}}$ and $\ln {m_{{x_1},t}}$ . The other two parameters ${s_{0,t}}$ and ${s_{1,t}}$ specify the smoothed pathway between the two mortality rates at the two boundary ages in year $t$ .

When the time dimension is included in modelling mortality data, avoiding any mortality crossover between related populations is not straightforward. Particularly, when the sample period covers more distant years, the observed mortality rates at the oldest ages are more volatile, and the data themselves may exhibit some crossover. Such features would not be helpful in determining the fitted mortality surface. We propose one possible solution to cope with this issue by imposing constraints on the second-level parameter ${\omega _t}$ . Consider males (population 1) and females (population 2) in the same country, when the two populations are assumed to share (only) the common gradient parameters ${s_{0,t}}$ and ${s_{1,t}}$ , for year t, the initial mortality differential $\left( {\alpha _t^{(1)} - \alpha _t^{(2)}} \right)$ should decrease to $\left( {\omega _t^{(1)} - \omega _t^{(2)}} \right)$ without changing its sign over the age range considered. This property is true if $\left( {\alpha _t^{(1)} - \alpha _t^{(2)}} \right)\left( {\omega _t^{(1)} - \omega _t^{(2)}} \right) \gt 0$ . Given $\alpha _t^{(1)} \gt \alpha _t^{(2)}$ for all $t$ , $\omega _t^{(1)} \gt \omega _t^{(2)}$ is required to avoid the potential crossover. In other words, the female population with lower mortality rates at the starting age is expected to maintain this “advantage” at the ending age. As we will discuss in more detail later, the maximum likelihood estimation with the Newton–Raphson iterative scheme is used to calibrate the model parameters. Nevertheless, the inequality constraints are not guaranteed to be satisfied during the estimation procedure. Accordingly, we adopt the barrier method (Wright & Nocedal, Reference Wright and Nocedal1999) under which the inequality constraints are incorporated using the so-called barrier functions that serve as a penalty term. It has been applied in mortality modelling by Li & Liu (Reference Li and Liu2020, Reference Li and Liu2021) to constrain their heat wave model parameters. As the parameter estimates approach the inequality bounds, the magnitude of the barrier functions rises to infinity, and so the log-likelihood is penalised more by such estimates. Then the optimised parameters are “forced” to fall within the boundaries.

For ease of exposition, consider an optimisation problem searching for parameter values that maximise a function $f\!\left( \theta \right)$ , subject to the constraint $\theta \gt a$ . Under the barrier method, the optimisation problem is redesigned as maximising the penalised function $f\!\left( \theta \right) - r\!\left( {\frac{1}{{\theta - a}}} \right)$ , where $B\!\left( \theta \right) = r\!\left( {\frac{1}{{\theta - a}}} \right)$ refers to the barrier function with $r \gt 0$ . Given a feasible initial value that over-satisfies the constraint, as $\theta $ decreases to its lower bound $a$ , the barrier function tends to infinity, which then contradicts the aim of maximisation. Accordingly, this formulation would avoid the violation of the parameter restriction. Similarly, the barrier functions for the constraints $\theta \lt b$ and $a \lt \theta \lt b$ can be constructed as $r\!\left( {\frac{1}{{b - \theta }}} \right)$ and $r\!\left( {\frac{1}{{b - \theta }} + \frac{1}{{\theta - a}}} \right)$ , respectively.

Setting the objective function to the Poisson log-likelihood function, multiple barrier functions can be embedded to manage the inequality constraints imposed on the population-specific $\omega _t^{\left( i \right)}$ . For example, when the male and female populations in the same country are considered, the following constraints are imposed:

\begin{align*}\left\{ {\begin{array}{*{20}{c}}{\omega _t^{(1)} \gt a_t^{\left( {1,2} \right)}}\\[6pt]{\omega _t^{(2)} \lt a_t^{\left( {1,2} \right)}}\end{array}} \right.\!,\end{align*}

where $\omega _t^{(1)}$ and $\omega _t^{(2)}$ are the male and female level parameters at the ending age, respectively, $a_t^{\left( {1,2} \right)}$ is calculated as the average of the initial values of $\omega _t^{(1)}$ and $\omega _t^{(2)}$ . Under the careful design of the barrier method, the no-crossover condition $\omega _t^{(1)} \gt \omega _t^{(2)}$ can be met. The specifications of the barrier functions, initial values, penalised log-likelihood function and the updating equations for the HS4 model are provided in Appendix A. As we will show in section 3.2, the above barrier constraints can be applied similarly to risk factors other than gender.

2.2 Gompertz model

Gompertz (Reference Gompertz1825) discovered that the age-specific hazard rate of adult mortality increases exponentially with age and proposed a two-parameter model

\begin{align*}{\mu _x} = B{c^x}\end{align*}

or equivalently

(9) \begin{align}ln\ {\mu _x} = {k_1} + {k_2}x.\end{align}

Since this discovery, actuaries and researchers have developed various modifications to the original designFootnote 6. For instance, Makeham (Reference Makeham1867) added an age-independent term to the above model for the age-specific hazard to account for deaths that do not result from senescence. Besides including more parameters, adjustments were made to the format of the original function. The so-called Perks model (Perks, Reference Perks1932) with four parameters adopted a logistic function, which is described as follows:

\begin{align*}{\mu _x} = \frac{{A + B{c^x}}}{{1 + D{c^x}}}.\end{align*}

The well-known CBD model (Cairns et al., Reference Cairns, Blake and Dowd2006) can be obtained by replacing the hazard rate with the death rate and setting $A = 0$ and $B = D$ under the Perks law and adding a time dimension. Specifically, the probability of death ${q_{x,t}}$ is modelled as:

\begin{align*}{q_{x,t}} = \frac{{{e^{{k_{1,t}} + {k_{2,t}}x}}}}{{1 + {e^{{k_{1,t}} + {k_{2,t}}x}}}}\end{align*}

or equivalently

\begin{align*}logit\;{q_{x,t}} = ln \left( {\frac{{{q_{x,t}}}}{{1 - {q_{x,t}}}}} \right) = {k_{1,t}} + {k_{2,t}}x.\end{align*}

It is evident that the CBD model can be regarded as an application of the Gompertz law (equation (9)) applied to a transformed measure, in modelling population mortality data with both the age and time dimensions. This model has been widely applied and extended in the literature and has been found particularly suitable for explaining mortality rates at higher ages (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009).

For comparison purposes, we utilise the logarithm of the central death rate rather than the logit probability of death for the two-dimensional Gompertz law in this paper. The log central death rate $\ln {m_{x,t}}$ of a life aged $x$ in year $t$ is expressed as

(10) \begin{align}ln\ {m_{x,t}} = {k_{1,t}} + {k_{2,t}}x,\end{align}

where ${k_{1,t}}$ refers to the level of the mortality curve in year $t$ and ${k_{2,t}}$ describes the speed at which mortality increases for each extra year of age in year t. The fitted log central death rates in a given year are assumed to grow linearly with age.

2.3 PS model

The two mortality models introduced above make assumptions about the functional form of mortality rates and only impose smoothness across adjacent ages. However, some academics are sceptical about specifying the mortality surface using defined functions. To address this concern, one branch of mortality models originated from the spline fitting technique that smooths mortality rates only based on features of data and avoids the identification of the model structure.

One of the most commonly applied methods in the actuarial literature is the PS model (Currie et al., Reference Currie, Durban and Eilers2004). The concept of PS was initially proposed by Eilers & Marx (Reference Eilers and Marx1996). They applied difference penalties on the coefficients of B-splines to strike a balance between the goodness of fit and smoothness of the fitted curve. This one-dimensional spline fitting approach was extended to a bivariate regression model in Eilers & Marx (Reference Eilers and Marx2003) and Eilers et al. (Reference Eilers, Currie and Durbán2006). Using the same notation as in section 2.1, $\boldsymbol{D}$ , $\boldsymbol{E}$ and $\boldsymbol{M}$ represent the matrices of the death counts, central exposures and central death rates, and $\boldsymbol{d}$ , $\boldsymbol{e}$ and $\boldsymbol{m}$ are the vector components corresponding to a given year $t$ .

The log death counts in year $t$ can be modelled as:

(11) \begin{align}ln\ \boldsymbol{d} = ln\ \boldsymbol{e} + ln\ \boldsymbol{m} = ln\ \boldsymbol{e} + \boldsymbol{Ba},\end{align}

where $\boldsymbol{B}$ is the $\left( {{n_a} \times k} \right)$ B-spline regression matrix containing polynomials connected at equally spaced knots, and $\boldsymbol{a}$ is the $\left( {k \times 1} \right)$ coefficient vector. The column dimension $k$ of $\boldsymbol{B}$ is equivalent to the sum of the polynomial degree $q$ and the number of (internal) knots ${n_{{d_x}}}$ (i.e., $k = q + {n_{{d_x}}}$ ). The more knots the regression employs, the more flexible the fitted curve is. One may adopt a sufficient number of knots for B-splines to provide more flexibility than needed and avoid over-fitting by constructing a difference penalty term that imposes smoothness (Eilers & Marx, Reference Eilers and Marx2002). Accordingly, the coefficient vector $a$ can be estimated by maximising the penalised log-likelihood

\begin{align*}{l_p} = l\!\left( {\boldsymbol{a};\, \boldsymbol{d}} \right) - \frac{1}{2}\boldsymbol{a}^{\prime}\boldsymbol{Pa}, \end{align*}

where $l\!\left( {\boldsymbol{a};\,\boldsymbol{d}} \right)$ refers to the log-likelihood of Poisson distribution and the penalty matrix $\boldsymbol{P} = \lambda {\boldsymbol{\Delta }}^{\prime}{\boldsymbol{\Delta }}$ is composed of the smoothing parameter $\lambda $ and difference matrix ${\boldsymbol{\Delta }}$ with order 2.

The PS approach can also be applied to fitting the mortality surface across ages and years. We denote the $\left( {{n_a} \times {k_a}} \right)$ and $\left( {{n_y} \times {k_y}} \right)$ B-spline basis matrices associated with age $x$ and year $t$ by ${\boldsymbol{B}_a}$ and ${\boldsymbol{B}_y}$ , respectively. The number of equally spaced knots can be determined separately for the age and year horizons based on the data range. Then the B-spline matrix in the two-dimensional setting is constructed as the Kronecker product

\begin{align*}{\boldsymbol{B}^{(\textbf{2})}} = {\boldsymbol{B}_y} \otimes {\boldsymbol{B}_a}.\end{align*}

Given mortality data across both ages and years, the PS model is expressed as:

(12) \begin{align}ln\ \boldsymbol{D} = \ln \boldsymbol{E} + \ln \boldsymbol{M} = \ln \boldsymbol{E} + \boldsymbol{B}^{(2)} \boldsymbol{a} = \ln \boldsymbol{E} + {\boldsymbol{B}_a} \boldsymbol{AB}^{\prime}_y,\end{align}

where $\boldsymbol{A}$ is the $\left( {{k_a} \times {k_y}} \right)$ matrix obtained by rearranging the coefficient vector $\boldsymbol{a}$ with a length of ${k_a}{k_y}$ . The parameters in the PS model are estimated by maximising the same Poisson log-likelihood function using a smoothing matrix that penalises the two dimensions of $\boldsymbol{A}$ separately. The penalty matrix is given as

\begin{align*}\boldsymbol{P} = {\lambda _a}\!\left( {{\boldsymbol{I}_{{k_y}}} \otimes \boldsymbol\Delta^{\prime}_a{\boldsymbol\Delta _a}} \right) + {\lambda _y}\!\left( {\boldsymbol\Delta^{\prime}_y{\boldsymbol\Delta _y} \otimes {\boldsymbol{I}_{{k_a}}}} \right),\end{align*}

where ${\boldsymbol{I}_{{k_y}}}$ ( ${\boldsymbol{I}_{{k_a}}}$ ), ${\lambda _a}$ ( ${\lambda _y}$ ) and ${\boldsymbol\Delta _a}$ ( ${\boldsymbol\Delta _y}$ ) are the identity matrix, smoothing parameter and difference matrix contributing to the columns (rows) of $\boldsymbol{A}$ .

We have briefly reviewed the theoretical background of the PS method, more technical details and improvements on the fitting algorithm can be found in Currie et al. (Reference Currie, Durban and Eilers2004, Reference Currie, Durban and Eilers2006) and Eilers et al. (Reference Eilers, Currie and Durbán2006).

3.1 Single-population fitting

To fit the above three models, we collect population mortality data of Australia, England and Wales, France, Japan and the United States from the Human Mortality Database (HMD, 2021). Specifically, the death counts and exposure to risk of both genders from 1950 to 2018 over ages 56–95 are employed for this analysis. Following Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002), we assume that the number of deaths is a random variable following a Poisson distribution with the mean equal to the expected number of deaths. This Poisson assumption enables a rigorous framework to estimate the model parameters by maximising the log-likelihood function

(13) \begin{align}l = ln\ L = \mathop \sum \limits_{x,t} \left( {{d_{x,t}}\ ln\ {{\hat d}_{x,t}} - {{\hat d}_{x,t}} - ln\! \left( {{d_{x,t}}!} \right)} \right),\;\end{align}

where ${d_{x,t}}$ represents the observed death counts. The fitted number of deaths ${{\hat d}_{x,t}}$ is calculated as the product of the observed exposures and fitted mortality rates $\left({\hat d_{x,t}} = {E_{x,t}}{\hat m_{x,t}}\right)$ . As an illustration, the fitted death count for lives aged $x$ in year $t$ under the HS4 model is ${\hat d_{x,t}} = {E_{x,t}}{\hat m_{x,t}} = {E_{x,t}}\;{\rm{exp}}\!\left( {{\alpha _t}{h_{00}}\!\left( {{x_k}} \right) + {\omega _t}{h_{01}}\!\left( {{x_k}} \right) + {s_{0,t}}{h_{10}}\!\left( {{x_k}} \right) + {s_{1,t}}{h_{11}}\!\left( {{x_k}} \right)} \right)$ . For a parameter $\theta $ , the Newton–Raphson method is adopted to update its estimate iteratively as ${\theta ^*} = \theta - \dfrac{{\partial l}}{{\partial \theta }}/\dfrac{{{\partial ^2}l}}{{\partial {\theta ^2}}}$ . We provide below a summary of estimation steps under the HS4 model.

1. (1) Assign initial values to the parameters under the HS4 model and compute using the initial estimates. Specifically, we set ${\hat s_{0,t}} = {\hat s_{1,t}} = 0$ . The initial values of the other two parameters ${\alpha _t}$ and ${\omega _t}$ are set to the observed mortality levels at the two boundary ages ${x_0}$ and ${x_1}$ of the sample range in year $t$ (i.e., ${\hat \alpha _t} = \ln {m_{{x_0},t}}$ and ${\hat \omega _t} = \ln {m_{{x_1},t}}$ ). For our sample, ${x_k} = \dfrac{{x - {x_0}}}{{{x_1} - {x_0}}} = \dfrac{{x - 56}}{{95 - 56}}$ .

2. (2) Update ${\alpha _t}$ as $\hat \alpha _t^* = {\hat \alpha _t} + {\dfrac{\mathop \sum \nolimits_x \!\left( {{d_{x,t}} - {{\hat d}_{x,t}}} \right){h_{00}}\!\left( {{x_k}} \right)}{\mathop \sum \nolimits_x {{\hat d}_{x,t}}{{\left[ {{h_{00}}\!\left( {{x_k}} \right)} \right]}^2}}}$ .

3. (3) Update ${\omega _t}$ as $\hat \alpha _t^* = {\hat \omega _t} + {\dfrac{\mathop \sum \nolimits_x \!\left( {{d_{x,t}} - {{\hat d}_{x,t}}} \right){h_{01}}\!\left( {{x_k}} \right)}{\mathop \sum \nolimits_x {{\hat d}_{x,t}}{{\left[ {{h_{01}}\!\left( {{x_k}} \right)} \right]}^2}}}$ .

4. (4) Update ${s_{0,t}}$ as $\hat s_{0,t}^* = {\hat s_{0,t}} + {\dfrac{\mathop \sum \nolimits_x \!\left( {{d_{x,t}} - {{\hat d}_{x,t}}} \right){h_{10}}\!\left( {{x_k}} \right)}{\mathop \sum \nolimits_x {{\hat d}_{x,t}}{{\left[ {{h_{10}}\!\left( {{x_k}} \right)} \right]}^2}}}$ .

5. (5) Update ${s_{1,t}}$ as $\hat s_{1,t}^* = {\hat s_{1,t}} + {\dfrac{\mathop \sum \nolimits_x \!\left( {{d_{x,t}} - {{\hat d}_{x,t}}} \right){h_{11}}\!\left( {{x_k}} \right)}{\mathop \sum \nolimits_x {{\hat d}_{x,t}}{{\left[ {{h_{11}}\!\left( {{x_k}} \right)} \right]}^2}}}$ .

6. (6) Calculate the log-likelihood function $l = \mathop \sum \nolimits_{x,t} \!\left( {_{x,t}\ln {{\hat d}_{x,t}} - {{\hat d}_{x,t}} - \ln \left( {{d_{x,t}}!} \right)} \right)$ using the parameters estimated in (2)–(5).

7. (7) Repeat steps (2)–(6) until the improvement in $l$ is small enough, e.g., the change in $l$ between two iterations is less than ${10^{ - 7}}$ .

This study calibrates the HS and Gompertz models using the Poisson updating approach described above, while the PS model is fitted and forecasted using an R package “MortalitySmoothFootnote 7”. Both fitting methods are based on the Poisson maximum likelihood estimation, while the latter uses a so-called generalised linear array modelling that acts on matrices (Camarda, Reference Camarda2012). To examine the goodness of fit of the candidates, we calculate two measures – Akaike information criterion (AIC) and Bayesian information criterion (BIC) – as below:

\begin{align*}AIC = - 2l + {n_p} \times 2\\[3pt]BIC = - 2l + {n_p}\ ln\ {n_d},\end{align*}

where $l$ refers to the Poisson log-likelihood function, and ${n_p}$ is the number of effective parameters. These two criteria penalise the likelihood by the number of parameters involved in the model with a penalty factor of 2 and the log of data counts ( ${n_d}$ ), respectively. The information criteria and residuals plots are provided in Table 1 and Figure 2, respectively.

Table 1.

Effective number of parameters ( ${n_p}$ ), log-likelihood values ( $l$ ), AIC and BIC values under the HS1-4, Gompertz, PS10 and PS3 models. The mortality models are ranked from 1 (the lowest AIC/BIC) to 6 (the highest AIC/BIC) and the ranks are given in brackets.

Some interesting remarks can be made from Table 1 that lists the model comparison criteria values for both genders of the ten datasets. Firstly, HS1-3 models always produce greater AICs and BICs than the full HS4 model. Note that ${\alpha _t}$ and ${\omega _t}$ are interpreted as the mortality rates at the two boundaries of the age interval. Essentially, the fitted mortality curves in each year under the HS1 model are assumed to have a shape determined by the two Hermite functions linking the starting and ending points. Such a specification may be too rigid for modelling population mortality data, which is vividly demonstrated by the worst performance of the HS1 model. Employing additional gradient factors can introduce more flexibility to the model and so can improve the goodness-of-fit, which can be seen from the better ranks of the HS2-4 models.

Comparing between the HS2 and HS3 models with one gradient factor ${s_{0,t}}$ and ${s_{1,t}}$ , respectively, the former only outperforms the latter for males in Australia and England and Wales. Since ${s_{0,t}}$ ( ${s_{1,t}}$ ) controls the pattern of the fitted mortality curves leaving from (approaching) the starting (ending) age, one may suspect that there have been more variations in the mortality development of older age groups. Under the full HS model, employing two gradient factors yields significantly better fitting performance than the other HS candidates. As a benchmark, the Gompertz model is specified by two time-varying components and gives smaller AIC and BIC values than the HS1-2 models in all cases. However, it consistently underperforms the HS4 model (and HS3 model in some cases).

Figure 2

Standardised deviance residuals plots under the HS4, Gompertz, PS10 and PS3 models for Australian males (top panel) and females (bottom panel).

The above five models all propose a functional form for mortality rates in the age dimension, while no smoothness assumption across periods is imposed. On the other hand, the PS approach produces a smoothed mortality surface using splines connected by knots. We consider the PS10 and PS3 models with a different number of knots. Specifically, the PS10 (PS3) refers to the PS model with one internal knot per 10 (3) data points in each dimension. The PS3 model tends to be the optimal choice among the seven candidates for males in four out of the five countries, while the HS4 model appears to be preferred for female populations. When fewer internal knots are adopted, the PS10 model tends to underperform the HS4 model but still performs better than the Gompertz model in all cases.

Figure 2 depicts the deviance residuals against age, period and cohort year for Australia. Only the HS4, Gompertz, PS10 and PS3 models, which have similar AIC/BIC values, are presented. It is apparent that the residuals are generally quite random under the HS4 and two PS models. However, those produced by the Gompertz method do exhibit obvious patterns, which may suggest the need for more parameters. Moreover, the residuals plotted against age under the four models (and the residuals against year under the PS models) show some “bars” at certain ages (in certain years). Since the HS4 and Gompertz models do not assume smoothed mortality rates in the period dimension, the time-varying parameters would capture the pattern in each year separately and lead to more scattered residuals against year. This difference is a direct consequence of different smoothness assumptions between the mortality models.

Based on these observations, one can see that imposing smoothness on both the age and time dimensions can lead to missing out some important patterns, though having the advantage of a smaller number of effective parameters. For those mortality models that only impose the smoothness assumption on the age dimension (HS and Gompertz models), more parameters are used to describe the temporal developments more adequately. It is noteworthy that an outstanding fitting performance does not guarantee accurate forecasting. As we will see later, the forecasting performance under the PS approach is subject to great uncertainty, though the fitting results are superior. Backtesting will be conducted and presented in section 5 to investigate the forecast accuracy of the mortality models.

3.2 No-crossover property of the HS models

As discussed in section 2.1, the key strength of the HS approach is that mortality crossover is avoided via the properties of Hermite basis functions. Mortality differentials due to risk factors other than gender decrease as age increases, portrayed by the level parameter ${\alpha _t}$ specifically for each population. When the gender effect is considered, mortality differentials diminish with increasing age but never vanish entirely. By imposing constraints on the level parameter specifically for each population $\left\{ {\begin{array}{*{20}{c}}{\omega _t^{(1)} \gt a_t^{\left( {1,2} \right)}}\\[6pt]{\omega _t^{(2)} \lt a_t^{\left( {1,2} \right)}}\end{array}} \right.\!,$ the HS models can capture such a feature even if there is occasional mortality crossover in the data. In contrast, the Gompertz and PS models are incapable of modelling this pattern.

As an illustration of this crossover issue, we choose the mortality data for ages between 56 and 108Footnote 8 and years from 1950 to 2018. Figure 3 presents the observed and fitted log central death rates for the United States in 2018 under the six models. The PS approach here adopts one internal knot per 10 data points in each dimension to alleviate the potential crossover issue caused by overfitting. It can be seen that the HS1-3 models tend to underestimate the mortality rates at younger ages. The constructed shape of the HS1 and HS2 curves fits poorly at the oldest ages. On the other hand, the Gompertz model proposes a linear mortality curve of age and fails to capture the curvature observed at the highest ages. By comparison, the fitted curves under the HS4 and PS models follow the historical patterns more closely. However, the curves fitted under the PS model adhere to the data very closely, potentially indicating overfitting and lack of smoothness, even after applying the penalty structure. Another concern from Figure 3 is the intersection of the fitted female and male central death rates under the Gompertz and PS models at advanced ages. Under the compensation law of mortality, one would expect the gap between two genders to decline gradually over increasing age. But there is no biological basis to justify any crossover. Overall, the four HS models can generate a smoothed mortality curve that captures the relationships between males and females in the same country. Given appropriate constraints imposed on $\omega _t^{\left( i \right)}$ , the no-crossover property can be ensured.

Figure 3

Historical (dots) and fitted (solid lines) log central death rates in 2018 for the United States.

In addition to gender, the no-crossover property of the HS method holds when more risk factors are considered. As explained in equation (7) and footnote 6, when common gradient parameters are assumed for associated populations, the fitted mortality curves from the HS models would demonstrate narrowed gaps without an unreasonable crossover. This assumption should only be imposed if one believes that the populations under study exhibit similar paths of reaching the boundary mortality rates over the age interval. We have illustrated that the HS4 model tends to produce a good fit for mortality rates of males and females in the United States. Next, the fitting performances of the HS4, Gompertz and PS models are investigated when two risk factors – gender and country – are employed. Since the five countries considered in this study may not be highly related, the common-gradient assumption can be inappropriate. Accordingly, mortality data between 2008 and 2018 of females and males in Scotland and England and Wales are used for this particular illustrationFootnote 9. Based on the observations on $\alpha _t^{\left( i \right)}$ of the four populations, the following restriction is required for the HS4 model to avoid the potential crossover:

\begin{align*}\omega _t^{\left( {m2} \right)} \gt \omega _t^{\left( {m1} \right)} \gt \omega _t^{\left( {f2} \right)} \gt \omega _t^{\left( {f1} \right)},\end{align*}

where $\omega _t^{\left( {m1} \right)}$ and $\omega _t^{\left( {f1} \right)}$ $\left(\omega _t^{\left( {m2} \right)}\right.$ and $\left.\omega _t^{\left( {f2} \right)}\right)$ represent the level parameter of the female and male populations in England and Wales (Scotland). To conform with the constraint format of the barrier method, the above inequality can be converted into:

\begin{align*}\left\{ {\begin{array}{*{20}{c}}{\omega _t^{\left( {m2} \right)} \gt a_t^{\left( {1,2} \right)}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}\\[5pt]{a_t^{\left( {2,3} \right)} \lt \omega _t^{\left( {m1} \right)} \lt a_t^{\left( {1,2} \right)}}\\[5pt]{a_t^{\left( {3,4} \right)} \lt \omega _t^{\left( {f2} \right)} \lt a_t^{\left( {2,3} \right)}}\\[5pt]{\omega _t^{\left( {f1} \right)} \lt a_t^{\left( {3,4} \right)}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}\end{array}} \right.\!,\end{align*}

where $a_t^{\left( {j,j + 1} \right)}$ is computed as the average of the ${j^{th}}$ and ${\left( {j + 1} \right)^{th}}\;\left( {j = 1,2,3} \right)$ highest initial values of the level parameter at the ending age. Note that other relationships among $\omega _t^{\left( i \right)}$ ’s can be assumed. For example, if one believes that populations with the same gender will converge to the same level at the ending age (regardless of the country), a common $\omega _t^{\left( m \right)}$ $\left(\omega _t^{\left( f \right)}\right)$ could be employed.

Figure 4 displays the observed and fitted mortality rates in 2018 under the three mortality models. The PS curves show unwarranted bends due to the lack of smoothness. Although the Gompertz model captures the decreasing mortality differentials over age, evident crossovers are observed at around age 100. Under the HS4 model, the fitted mortality differentials for populations with the same gender in the two countries (red and green curves, blue and purple curves) tend to reduce towards a small gap at advanced ages. For males and females in the same country, the differences in fitted mortality rates diminish with increasing age but never vanish entirely. To conclude, the HS4 model seems to strike a good balance between the goodness-of-fit and smoothness. More importantly, the HS4 model can describe the law of mortality properly and avoid the unjustified mortality crossover that presents under the Gompertz and PS models.

Figure 4

Fitted log central death rates in 2018 of females and males in Scotland and England and Wales.

4.1 Time series models

This study now proceeds to examine the long-term performance of the six mortality models. The PS model extrapolates mortality rates not covered by the sample range during the fitting process by treating them as missing values. The R package “MortalitySmooth” offers the facility of both model estimation and projection, given a predetermined penalty order. In this study, we follow the suggestion by Currie et al. (Reference Currie, Durban and Eilers2004) and use the default penalty order of 2 for both the age and period dimensions, which implies that (log) mortality rates are predicted via (approximately) linear extrapolation. Under the HS and Gompertz models, the time-varying parameters are modelled by time series processes and projected to a future period to predict mortality levels. We first conduct a unit root test for each of the period factors ${\alpha _t}$ , ${\omega _t}$ , ${s_{0,t}}$ , ${s_{1,t}}$ , ${k_{1,t}}$ and ${k_{2,t}}$ .

The results suggest that all the period factors are not stationary at 5% significance level under the five models except for the female ${\omega _t}$ of Australia under the HS3 model, female ${\omega _t}$ of France under the HS4 model, female ${s_{1,t}}$ of Japan under the HS3 and HS4 models and male ${s_{0,t}}$ of Japan under the HS4 model. When using a 1% significance level, the unit root is detected for all the time series processes. We decide to fit the time-varying factors in the same mortality model by a multivariate (bivariate) random walk with drift (RWD):

(14) \begin{align}{\boldsymbol\theta _{\boldsymbol{t}}} = {\boldsymbol\theta _{\boldsymbol{t} - \textbf{1}}} + \boldsymbol\mu + {\boldsymbol{Z}_{\boldsymbol{t}}},\end{align}

where ${\boldsymbol\theta _{\boldsymbol{t}}}$ is a $\left( {n \times 1} \right)$ vector of the $n$ time-varying factors in a mortality model in year $t$ , $\boldsymbol\mu $ refers to the corresponding drift vector, ${\boldsymbol{Z}_{\boldsymbol{t}}}$ contains $n$ multivariate Gaussian error terms with mean zero. Although other time series choices can be considered, the RWD model has been widely adopted in the literature to extrapolate the decreasing trend of mortality (e.g., Cairns et al., Reference Cairns, Blake and Dowd2006; Dowd et al., Reference Dowd, Cairns, Blake, Coughlan, Epstein and Khalaf-Allah2010c). In addition, the assumed linear extrapolation under the multi-dimensional RWD process allows for a fair comparison between the PS approach with the second-order penalty and the other five models.

Figure 5

Estimated (solid lines) and projected (dashed lines) period factors in the HS1-4 and Gompertz models for males in England and Wales.

Figure 5 displays the estimated and projected parameters in the HS1-4 and GompertzFootnote 10 models for males in England and Wales. As illustrated, the two level parameters ${\alpha _t}$ and ${\omega _t}$ show a clear declining trend under the HS models, though the downward tendency in ${\omega _t}$ becomes less prominent under the HS3 model. This observation aligns with our expectations because these parameters are interpreted as the fitted mortality rates at the youngest and oldest ages of the sample which have been improving over time. Furthermore, the behaviour of ${s_{0,t}}$ and ${s_{1,t}}$ varies to some extent under different models. For example, the estimated ${s_{0,t}}$ values under the HS2 model (second row) demonstrate a generally decreasing trend over time. However, under the full model (fourth row), the declining tendency becomes less clear and more fluctuations can be found in the estimated ${s_{0,t}}$ values. Similarly, although the estimated ${s_{1,t}}$ values have been increasing since around 1980 under both the HS3 (third row) and HS4 (fourth row) models, the latter tends to show more fluctuating patterns. Besides, one may notice that the intercept and slope parameters ${k_{1,t}}$ and ${k_{2,t}}$ in the Gompertz model (first row) show an opposite trend to each other. It indicates that the average mortality level has been improving, and this reduction has been larger at the younger ages. Though not presented here, the parameter estimates exhibit similar trends for females and males in the five countries. Then we model the time-varying parameters as the multi-dimensional RWD. The estimated drift terms, standard deviations and correlation matrices are presented in Table B.1 in Appendix B.

4.2 Long-term projections

After obtaining the time series forecasts, future mortality rates can be calculated by updating the predicted time-varying components. It is of interest to investigate the impact of employing different model structures on mortality forecasts. We calculate life expectancies at key ages over the projection period. Since our data only cover ages up to 95, the Coale–Kisker approach (Coale & Kisker, Reference Coale and Kisker1990) is applied to extend the mortality rates to an “ultimate” age 110 with a presumed log central death rate of 0.7 (Gampe, Reference Gampe2010)Footnote 11. Figure 6 shows the Australia and England and Wales observed life expectancy at age 80 ( ${e_{80}}$ ), with projections up to the year 2050.

Figure 6

Observed and projected life expectancy at age 80 for Australia and England and Wales from 2010 to 2050. For demonstration purposes, the projections under the HS4 and Gompertz models are displayed in dashed lines.

It is apparent that the PS10 model (yellow curves) predicts significantly longer life expectancy than the other model candidates for Australia. We remark that the projections produced under the PS10 model become particularly unreasonable when mortality data of England and Wales are employed – declining life expectancies over time. As demonstrated in Figure B.1 in Appendix B that plots the fitted and extrapolated mortality rates at ages 60, 70 and 80, the projected mortality rates for England and Wales exhibit an increasing tendency, which is not the case for Australia. The potential reason could be the slowing down of mortality improvement over the most recent few years in England and Wales. Since mortality rates are not predicted by selecting appropriate time series processes to ensure biological reasonableness, it is possible to observe unrealistic patterns in the forecasts under the PS approach.

The potentially unstable projections of the PS model were highlighted in previous studies. For example, based on Italian populations, Cocevar (Reference Cocevar2007) obtained insensible mortality forecasts for females, while the performance is reasonable for males. Richards et al. (Reference Richards, Ellam, Hubbard, Lu, Makin and Miller2007) examined the prediction performance of the two-dimensional PS method using mortality data of seven countries and found erratic forecasts in some cases. They argued that the PS method tends to balance irregular effects in data with smoothing and those features are carried through to projections.

We also investigate the long-term prediction performance under the PS model using different number of knots, while there is no optimal knot choice that can produce sensible forecasts for all datasets. In addition, though excessive knots are penalised to generate a reasonable effective dimension, they tend to yield even more volatile forecasts. Therefore, only the prediction performance of the PS10 model is presented in the remainder of the paper.

Besides this noticeable (and unreasonable) deviation under the PS method, there exist some differences among the other models. Firstly, the HS1 model that obtains future mortality rates by projecting the mortality levels at the initial (56) and ending (95) ages results in the longest life expectancy of the four HS models. As depicted in Figure 6, adding one or two gradient parameters tends to reduce the predicted values. For example, the green curves under the HS3 model almost always lie below the other four sets of predictions (excluding the PS method) for the two ages. This can be attributed primarily to the less linear trend of the second level parameter ${\omega _t}$ . For example, it can be seen from Figure 5 that the projected values of ${\hat{\omega} _t}$ (the projected $ln\ {\widehat{m}}_{x_{1},t}$ ) from the HS3 model are generally higher than those from the other HS candidates, while the HS3 $ {\widehat{\alpha}}_{t}$ (the projected $ln\ {\widehat{m}}_{x_{0},t}$ ) values are similar to those from the HS2 model. One would expect the former to predict longer life expectancy than the latter. Although not presented here, for all the datasets the projections under the HS2 (blue) and HS4 (purple) models tend to lie within the two boundary curves formed by the other two HS candidates for all the datasets. Gompertz projections tend to be similar to those for HS1, with differences of less than 0.1–0.2 years in 2050.

Overall, the HS1-4 and Gompertz models give sensible forecasts of life expectancies. The PS model predicts unreasonably high or decreasing life expectancies, possibly caused by the automatic extrapolation process that may have overfitted the data and misused the random noises therein.