2.1 Long memory for time series of counts

In time series literature, the modelling of time series persistence over long lags is termed long memory. There are two typical types of long memory model structure one often considers, the Fractional autoregressive integrated moving average (ARIMA) model known as the autoregressive fractionally integrated moving average (ARFIMA) model class and the more general Gegenbauer ARIMA class of models. This latter class of models simply gets its name from the fact that the representation of long memory time series structure is expressed naturally in terms of Gegenbauer polynomials as will be detailed further mathematically below. Granger & Joyeux (Reference Granger and Joyeux1980) and Hosking (Reference Hosking1981) extended the classical ARIMA model to the ARFIMA model. Hosking (Reference Hosking1981) further extended to a generalised ARFIMA model called Gegenbauer ARMA (GARMA) model, which can model data with an oscillatory damping autocorrelation function (ACF).

An important stylised fact of mortality data by age groups is the presence of a strong persistence in observed death counts. This is easily discernible by simply evaluating the sample autocorrelation structures of the time series of death counts for each age group at a national level. This feature is almost universally observed in all age groups for the majority of countries that have reporting records in the Human Mortality Database (HMD) (Shkolnikov, Reference Shkolnikov2017) and is present to differing degrees in each gender and across age groups. Based on this observation, we therefore seek to incorporate such a structure into a GLM regression model based around extensions of the LC model framework. This will involve modelling the strong autocorrelation persistence in the form of long memory time series structures. If such features are statistically meaningful in helping to more accurately explain mortality projection in-sample and for out of sample forecasting, one would expect improved model selection for such model features and more accurate forecasts. We indeed show this is the case. We will begin by explaining the GLM extension to incorporate this structure.

2.1.1 Modelling persistence via latent linear predictors in regression

In time series settings, extended temporal dependence or persistence in a time series is termed long memory. Given a stationary time series process ${\textbf{\textit{Y}}}_{1:T} \equiv (Y_1, Y_2,\dots, Y_T)$, with ${\textbf{\textit{Y}}}_{1:T} \in (\mathbb{N} \cup \{0\})^T$, Beran (Reference Beran1994) defined a condition for a long memory stationary process in terms of the divergence of the ACF for $Y_t$ and $Y_{t+j}$ at lag j, such that

\begin{equation*}\underset{n\rightarrow \infty }{\mbox{lim}}\sum_{j=-n}^{n}|\rho (j) |\rightarrow \infty \ \ \mbox{where} \ \ \rho(j)=\frac{\mbox{Cov}\left(Y_t,Y_{t+j}\right)}{\sqrt{\mathbb{V}\mbox{ar}(Y_t)\mathbb{V}\mbox{ar}\left(Y_{t+j}\right)}}\end{equation*}

The particular form of trend regression we consider involves an extension of the classical ARMA model structures typically used for period and cohort effects to accommodate long memory. This naturally brings one to consider the family of ARFIMA (p, d, q) models with mean $\mu$ as proposed by Granger & Joyeux (Reference Granger and Joyeux1980) and Hosking (Reference Hosking1981), which describe a long memory stationary process with integrated order $d \in (0,1/2)$.

Definition 2.1 (ARFIMA). Consider a stationary time series process with constant $c \in \mathbb{R}$, an ARFIMA model with order (p, d, q) is defined by

\begin{equation*}\Phi (B)(Y_t-c )=\Theta (B)(1-B)^{-2d}\varepsilon _t, \ \ \varepsilon _t\overset{\mbox{i.i.d}}{\sim }N\left(0,\sigma _\varepsilon ^2\right)\end{equation*}

where B is the backshift operator, such that $B Y_t=Y_{t-1}$ and

(1) \begin{equation}\Phi (B) = 1-\phi_{1}B- \cdots - \phi_{p}B^{p} \ \ \mbox{and} \ \ \Theta (B) = 1+\theta_{1}B+ \cdots +\theta_{q}B^{q} \end{equation}

are the autoregressive and moving average characteristic polynomials, respectively, with no common roots.

For $d \in (0,1/2)$, the ARFIMA process is said to exhibit long memory. We note that ARFIMA model will become a short memory ARMA model for $d=0$ with the long memory operator $(1-B)^{-2d}=1$. A useful GARMA (p, d, q) model was proposed by Hosking (Reference Hosking1981) to describe data showing a slowly damping ACF with an oscillatory pattern.

Definition 2.2 (GARMA) Consider a stationary time series process with constant $c \in \mathbb{R}$, a GARMA model with order (p, d, q) is defined by

(2) \begin{equation}\Phi (B)(Y_t-c )=\Theta (B)\left(1-2uB+B^2\right)^{-d}\varepsilon _t\equiv \Theta (B) \left ( \sum_{j=0 }^{\infty }\psi _{j}\varepsilon_{t-j} \right ), \ \ \varepsilon _t\overset{\mbox{i.i.d}}{\sim }N\left(0,\sigma _\varepsilon ^2\right) \end{equation}

where $\Phi (B)$ and $\Theta (B)$ are defined in equation (1). $\psi _{j}$ denote the coefficients of the generating function for the Gegenbauer polynomials $(1-2uB+B^{2})^{-d}$ (Stein & Weiss 1971). These coefficients are formulated as

(3) \begin{equation}\psi _{j}=\sum_{q=0}^{[j/2]}\frac{({-}1)^{q}(2u)^{j-2q}\Gamma (d-q+j)}{q!(j-2q)!\Gamma (d)} \end{equation}

where $[j/2]$ represents the integral part of $j/2$.

The coefficients $\psi_j$ in equation (3) are functionally dependent on d that controls the strength of long memory and the Gegenbauer parameter u that controls the oscillation of ACF (Rainville, Reference Rainville1960). The ARFIMA model is a special case of the GARMA model when $u = 1$, such that the operator $(1-2uB+B^2)^{-d}=(1-2B+B^2)^{-d}=(1-B)^{-2d}$. For a long memory process, ACF declines at a hyperbolic rate to zero, which is a much slower rate of decay than the exponential decay of the ARMA process. The long memory strength is controlled by the value of d and the value of u determines the cycle of the long memory process. Investigating the behaviours of ACF with different values of d and u helps actuaries to understand how long memory improves forecasting performance. In particular, when $u=-1$ and $0<d<1/4$, the ACF can be written as

\begin{equation*}\rho (j) = ({-}1)^{\,j}\frac{\Gamma (1-2d)\Gamma (j+2d)}{\Gamma (2d)\Gamma (j-2d+1)} \sim \mbox{constant}\cdot ({-}1)^{\,j} j^{4d-1} \mbox{, as} \ \ j\rightarrow \infty \nonumber\end{equation*}

For ARFIMA(0, d, 0) with $u=1$ and $0<d<1/4$, it is given asymptotically by

\begin{equation*}\rho (j) = \frac{\Gamma (1-2d)\Gamma (j+2d)}{\Gamma (2d)\Gamma (j-2d+1)} \sim \mbox{constant}\cdot j^{4d-1} \mbox{, as} \ \ j\rightarrow \infty \nonumber\end{equation*}

In addition, we note that the closed form ACF for the GARMA(0, d, 0) model with the constraint given by $|u|<1$, $0<d<1/2$ is unavailable (Woodward et al., Reference Woodward, Cheng and Gray1998), but is asymptotically given by

\begin{equation*}\rho (j) \sim \mbox{constant}\cdot j^{2d-1}{\sin}(\pi d-j\lambda _{0}) \mbox{, as} \ \ j\rightarrow \infty\end{equation*}

where $\lambda_0=\cos^{-1}(u)$. For different values of d and u, the ACF plots of long memory time series show different patterns. Figures 1 and 2 show that the long memory strength is controlled by the value of d and the value of u determines the oscillation pattern of long memory process as displayed in the ACF.

Figure 1 The ACF plot for simulated data with different values of d and $u=0.9$ for GARMA model.

Figure 2 The ACF plot for simulated data with different values of u and $d=0.49$ for GARMA model.

2.2 Incorporating long memory in extended GLM mortality models

The extensions of GLM regression models, in which the linear predictor incorporates a time series structure, are often distinguished from standard GLM regression structures by the terminology generalised linear GARMA (GLGARMA) to denote that they contain a linear predictor with an Autoregressive Moving Average (ARMA) structure, see Davis et al., (Reference Davis, Dunsmuir and Wang1999). Examples of such structures in mortality modelling include latent process dynamics in the trend, such as in the LC model with a period effect. One can also think of the GLARMA models as a subclass of generalised state-space models for non-Gaussian time series, see discussions in Fung et al., (Reference Fung, Peters and Shevchenko2017).

In this paper, we are interested in the class of long memory GLMs with time series structure in the linear predictor which are denoted by the terminology GLGARMA models. The “GARMA” component refers to a type of long memory ARMA model based on Gegenbauer polynomials that act as coefficients in the GARMA model to express the persistence through a fractional difference operator.

Hence, the GLGARMA model is basically a generalised state-space model for a time series ${\textbf{\textit{Y}}}_{1:T} \equiv (Y_1, Y_2,\dots, Y_T)$ consisting of an observation variable (death counts) and state variables (linear predictors with latent process structures for period effects and cohort effects). The canonical log link is adopted to link this function to the mean of a count distribution, to ensure the intensity is strictly positive.

In order to model discrete time series with long memory features, Yan et al., (Reference Yan, Chan and Peters2017) extended the GLM to GLGARMA model which combines the GLM and the GARMA time series model.

Definition 2.3 (GLGARMA model). Given a discrete stationary time series process and a natural $\sigma$-algebra for the observed data filtration $\mathcal{F}_{1:t-1}=\sigma( Y_{1}, Y_{2},\cdots, Y_{t-1})$, a GLGARMA model with order (p, d, q) is defined by

\begin{eqnarray*}Y_{t} |\mathcal{F}_{1:t-1} & \sim & \mbox{GP}(\mu_t, \nu), \nonumber \\\Phi (B)(\ln(\mu _{t})-c) &=& (1-2uB+B^{2})^{-d} \Theta (B)\varepsilon _t \equiv \Theta (B) \sum_{j=0 }^{\infty }\psi _{j}\varepsilon _{t-j} \nonumber\end{eqnarray*}

Sigma algebra on a set X is a collection of subsets of X that includes X itself and is closed under complement and countable unions. Data filtration is the process of choosing a subset of the data set for analysis. So natural sigma algebra $\mathcal{F}_{1:t-1}$ for the observed data filtration means the analysis requires any subset of a given data $(Y_1, Y_2,\cdots , Y_{t-1})$ ($Y_t \in (\mathbb{N} \cup \{0\})$) including empty set. We have included this information in P.8, Definition 2.3. In this paper, the data filtration $\mathcal{F}_{1:t-1}$ corresponds to observed realisations of time series up to time $t-1$.

The GP distribution, GP$(\mu,\nu)$ has the probability mass function (pmf), mean and variance given by Definition 2.4.

Definition 2.4 (GP). A random variable Y follows a GP distribution with support on $(\mathbb{N} \cup \{0\})$ if it has a pmf, mean and variance given by, respectively

\begin{eqnarray*}f(y;\ \mu ,\nu )& = & \mu(1-\nu) [\mu(1-\nu) +\nu y]^{y-1}e^{-\mu(1-\nu)-\nu y}/y!,\quad \mu>0, \ -1 \leq \nu<1, \\ \mathbb{E}(Y)& = & \mu \ \ \textrm{and} \ \ {\mathbb{V}}ar(Y)=\mu(1-\nu)^{-2}\end{eqnarray*}

GP distribution is over-, under- and equi-dispersed when the dispersion parameter $\nu \in ({-}1,1)$ is greater than, less than and equal to 0, respectively.

We develop new mortality models for death counts by combining GLGARMA model in Definition 2.3 and LC model. Let the random vector ${\textbf{\textit{Y}}}_{x,1:T}=(Y_{x,1},Y_{x,2},\cdots,Y_{x,T})$ denote the death count series for age group $x \in \left\{ x_1,\cdots , x_g \right\}$ with $Y_{x,t} \in (\mathbb{N} \cup \{0\})$ and year $t \in \left\{1,\cdots , T \right\}$. The natural sigma algebra for the observed data filtration is defined as

(4) \begin{equation}\mathcal{F}_{x,1:t-1}=\sigma\left( Y_{x,1}, Y_{x,2},\ldots, Y_{x,t-1}\right) \end{equation}

The vector of central exposures to death risk in the middle of age interval x is ${\textbf{\textit{E}}}_{x,1:T}= \left(E_{x,1}, E_{x,2}, \ldots, E_{x,T}\right)$.

2.2.1 Single age group long memory Lee–Carter model

The Gegenbauer long memory LC model for single age group is given by the following definition.

Definition 2.5 (Proposed model of order (p, d, q)). Condition on the central death rate $ \mu_{x,t}(\kappa_{x,t})$ and central exposure $ E_{x,t}$ where $\kappa_{x,t}$ denotes the period effect, we assume that $Y_{x,t} |{\mathcal{F}_{x,1:t-1}},\break E_{x,t},\mu_{x,t}(\kappa_{x,t})$ forms a Markov process,

(5) \begin{eqnarray}Y_{x,t} |{\mathcal{F}_{x,1:t-1}},E_{x,t},\mu_{x,t}(\kappa_t) &\sim& {\textrm{GP}}( E_{x,t}\mu_{x,t}(\kappa_t),\nu_{x} ), \nonumber \\ \Phi_x(B){\ln \mu_{x,t}} &=& \mathbb{B}_x(t;\ r,l, {\boldsymbol{\tau}})+\kappa_{x,t}+ (1-2u_xB+B^{2})^{-d_x} \Theta_x (B)( \varepsilon _{x,t}),\\ \kappa_{x,t} & =& \eta_x \kappa_{x,t-1} +\varepsilon_{x,t}^{\kappa}, \nonumber \\ \mathbb{B}_x(t;\ r,l, {\boldsymbol{\tau}}) & = & \sum_{i=1}^{l} \left ( a_{x,i}^{(0)}+a_{x,i}^{(1)}(t-\tau_i)\right.\nonumber\\ && \qquad + \left. a_{x,i}^{(2)}(t-\tau_i)^2+a_{x,i}^{(3)}(t-\tau_i)^3 +\cdots a_{x,i}^{(r)}(t-\tau_i)^r \right )\cdot \mathbf{1}_{t \in [\tau_{i},\tau_{i+1}]} \nonumber \end{eqnarray}

where $0<d_x<1/2$, $|u_x|<1$, $\Phi_x(B)$ and $\Theta_x(B)$ are defined in equation (1) with age-specific coefficients for each age group x. The error terms are defined as $\varepsilon _{x,t} \overset{\mbox{i.i.d}}{\sim } {N}(0,\sigma_{x,\varepsilon}^2)$ and $\varepsilon_{x,t}^{\kappa} \overset{\mbox{i.i.d}}{\sim } {N}(0,\sigma_{x,\varepsilon^{\kappa}}^2\!)$. In order to present a graduation (time) effect, the constant c is replaced by the B-spline $\mathbb{B}_x(t;\ r,l, {\boldsymbol{\tau}})$. This temporal spline is of order r with l knot points ${\boldsymbol{\tau}}=(\tau_1,\ldots,\tau_l)$ and it can be reduced to an intercept $\alpha_x$ as a special case. Furthermore, the period effect $\kappa_{x,t}$ for each age group is a latent AR(1) process in which we enforce stationarity through $|\eta_x|<1$.

Remark 2.1. We remark that higher order variants of the latent intensity such as AR(p) models for the period effect can easily be accommodated. In practice, we find it suitable to focus our attention here on the classical LC-type model which has an AR(1) structure, except we extend the classical LC model by restricting ourselves to strictly second-order stationary processes.

Considering some special cases of equation (5), the long memory LC model can be further divided into eight different sub-models with their mean functions defined in Table 2. In Appendix A, we consider these eight sub-models and perform simulation studies to assess the accuracy of the parameter estimates for these models.

Table 2. Eight different sub-models where models 1,4,8 are graduation models, models 3,4,7,8 are LC-type models, models 1–4 are short memory models and models 5–8 are long memory models

2.3 Relationship between LC period effect models and long memory extensions

In this section, we explain the connection between LC model with a non-stationary period effect and the long memory model. There are several disturbance terms in our proposed models, which come from the latent state period effect $\kappa_{x,t}$ specification. These disturbance terms in a linear univariate structural time series model can be combined in an appropriate manner to obtain a reduced form canonical time series, accommodating the fractional difference operator that produces the persistence in the long memory component of the model (Harvey, Reference Harvey1990). For instance, Hyndman et al., (Reference Hyndman, Koehler, Ord and Snyder2008) showed that an ARIMA form is obtained by eliminating the state variables from a linear innovation state-space model. Such a reduced form when available can be beneficial for interpreting aspects of the state-space model such as the ACF and the autocovariance function, which in some cases may be obtained in closed form.

To see why this is possible, we note that such a form can be easily represented by Wold representation (Wold, Reference Wold1938) and hence variance, autocorrelation structure and the latent process can be obtained. The Wold theorem plays a central role in time series analysis, and it indicates that the dynamic of any purely non-deterministic covariance-stationary process can be arbitrarily well approximated by an ARMA process. This simplified form helps us to understand better the impact of long memory features on the period effect. We study the reduced form of the mean function in model 7 as an example since other models are just special cases. Model 7 has the following structure:

(6) \begin{eqnarray}\Phi_x(B){\ln \mu_{x,t}} & = & \alpha_x+\kappa_{x,t}+ (1-2u_xB+B^{2})^{-d_x} \Theta_x (B) \varepsilon _{x,t},\\ \kappa_{x,t} & = & \eta_x \kappa_{x,t-1} +\varepsilon_{x,t}^{\kappa} \nonumber \end{eqnarray}

We have

(7) \begin{align}\kappa_{x,t}-\eta_x \kappa_{x,t-1} & = (1- \eta_x B)\kappa_{x,t} = \varepsilon_{x,t}^{\kappa} \Rightarrow \kappa_{x,t} = \frac{\varepsilon_{x,t}^{\kappa}}{(1- \eta_x B)},\nonumber\\[4pt] \kappa_{x,t} -\kappa_{x,t-1} & = \frac{(1-B)}{(1- \eta_x B)}\varepsilon_{x,t}^{\kappa}\end{align}

By taking differences of equation (6) with respect to time t and $t-1$, we obtain

\begin{equation*}\Phi_x(B)(\ln \mu_{x,t}-\ln \mu_{x,t-1}) = (\kappa_{x,t}- \kappa_{x,t-1})+ (1-2u_xB+B^{2})^{-d_x} \Theta_x (B)( \varepsilon _{x,t}- \varepsilon _{x,t-1}) \nonumber \\\end{equation*}

By substituting equation (7) and eliminating $(1-B)$, we obtain the following reduced form:

\begin{eqnarray}\Phi_x(B)\ln \mu_{x,t} & = & \frac{1}{(1- \eta_x B)} \varepsilon_{x,t}^{\kappa}+ (1-2u_xB+B^{2})^{-d_x} \Theta_x (B) \varepsilon _{x,t} \nonumber\end{eqnarray}

For the case that $\varepsilon_{x,t}^{\kappa}$ and $\varepsilon _{x,t}$ are both i.i.d white noise, we have a transformed GARMA model represented as

\begin{eqnarray}\Phi_x(B)\ln \mu_{x,t} & = & \left [ \frac{1}{(1- \eta_x B)}+ (1-2u_xB+B^{2})^{-d_x} \Theta_x (B) \right ] \varepsilon _{x,t} \nonumber\end{eqnarray}

Based on the generating function $\displaystyle \frac{1}{(1-cx)}=\sum_{i=0}^{\infty}(cx)^i$ for constant c and equation (2), we have

\begin{equation*}\frac{1}{(1- \eta_x B)} = \sum_{i=0}^{\infty}(\eta_x^i)B^i \ \ \mbox{and} \ \ (1-2u_xB+B^{2})^{-d_x} = \sum_{i=0 }^{\infty }\psi_{x,i}B^i\nonumber\end{equation*}

For a transformed GARMA (0, d, 0) model, we have

(8) \begin{eqnarray}\ln \mu_{x,t} & = & \sum_{i=0 }^{\infty }(\eta_x^i+\psi_{x,i}) \varepsilon _{x,t-i} \end{eqnarray}

Figure 3 shows the impact of the period effect on the ACF and spectral density function of GARMA model. The spectral density represents an autocovariance function in the frequency domain, which is obtained via Fourier transform. This frequency domain approach involves projection of the signal onto a sinusoidal basis and is then able to capture and describe the oscillatory features of data in an interpretable manner. For the GARMA model, the peak indicates that a substantial source of variability is contributed by this frequency band around the peaks. Furthermore, the value of u can be easily determined by the location of the peak. For the transformed GARMA model in equation (8), the strength of the long memory is reduced and the spectral density plots show that the long memory pattern becomes less symmetric around the dominant frequency.

Figure 3 The ACF plot and spectral density plot for GARMA model with $d_x=0.45$ and $u_x=0.9$ and transformed GARMA model with $d_x=0.45$, $u_x=0.9$ and $\eta_x=0.9$.

3.1 Posterior predictive distribution

One clear advantage of a Bayesian approach is that it replaces the high-dimensional integrals in evaluating the marginal likelihood functions in the ML approach with posterior sampling. This advantage particularly applies to our proposed model because the model involves latent variables in both mean $\mu_{x,t}$ and period effect $\kappa_{x,t}$ functions. In this section, we introduce the Bayesian approach to estimate the proposed model in Definition 2.5. The basic idea of the Bayesian approach can be described by the following definition and equations.

Definition 3.1 (Posterior distribution). Consider a set of observed values ${\textbf{\textit{y}}}_{x,1:T} = (y_{x,1},\break y_{x,2}, \ldots, y_{x,T})$ with each $y_{x,t} \in (\mathbb{N} \cup \{0\})$, central exposure ${\textbf{\textit{E}}}_{x,1:T}$ and the vector of unknown parameters

\begin{equation*}{\boldsymbol{\vartheta}}_x =\left(\alpha_x, a_{x,i}^{(0)}, \cdots, a_{x,i}^{(r)}, \eta_x, \sigma^2_{x,\varepsilon}, \sigma^2_{x,\varepsilon^{\kappa}}, d_x, u_x,\nu_x \right)\end{equation*}

and denote the state parameters ${\textbf{\textit{z}}}_{x,t}=[\mu_{x,t}(\kappa_{x,t}), \kappa_{x,t}]$ and ${\textbf{\textit{z}}}_x =[ {\textbf{\textit{z}}}_{x,1},\cdots, {\textbf{\textit{z}}}_{x,T}]$ and

\begin{equation*}{\boldsymbol{\vartheta}}_{x}^*=[{\boldsymbol{\vartheta}}_x, {\textbf{\textit{z}}}_x]\end{equation*}

the posterior distribution for $\boldsymbol{\vartheta}_x^*$ conditional on ${\textbf{\textit{y}}}_{x,1:T} $ is given by

\begin{equation*}\pi\!\left(\boldsymbol{\vartheta}_x^* | \mathcal{F}_{x,1:T} \right)=\frac{f\!\left({\textbf{\textit{y}}}_{x,1:T} | {\textbf{\textit{E}}}_{x,1:T},\boldsymbol{\vartheta}_x^*\right)\ \pi\!\left(\boldsymbol{\vartheta}_x^*\right)}{\int f\!\left({\textbf{\textit{y}}}_{x,1:T} | {\textbf{\textit{E}}}_{x,1:T},\boldsymbol{\vartheta}_x^*\right)\ \pi\!\left({\boldsymbol{\vartheta}}_x^*\right) \ d\boldsymbol{\vartheta}_x^*}\propto f\!\left({\textbf{\textit{y}}}_{x,1:T}| {\textbf{\textit{E}}}_{x,1:T},\boldsymbol{\vartheta}_x^*\right)\pi\!\left(\boldsymbol{\vartheta}_x^*\right)\end{equation*}

where the prior densities $\pi(\boldsymbol{\vartheta}_x)$ can be chosen based on available information or past data.

Even if there is no prior information, non-informative or reference priors can be also applied. Furthermore, the credible intervals for all parameters of interest can be constructed from posterior distributions. A credible interval is the Bayesian equivalent of the confidence interval in frequentist statistics. Credible intervals capture our current uncertainty in the location of the parameter values and thus can be interpreted as a probabilistic statement about the parameter. In contrast, confidence intervals capture the uncertainty about the interval we have obtained (i.e. whether it contains the true value or not). Thus, they cannot be interpreted as a probabilistic statement about the true parameter values. The credible interval is defined in the domain of a posterior probability distribution or a predictive distribution. Under the model structure in equation (5), the complete likelihood function for model 8, as an example, is

\begin{eqnarray*} f\!\left({\textbf{\textit{y}}}_{x,1:T}| {\textbf{\textit{E}}}_{x,1:T}, \boldsymbol{\vartheta}_x^*\right)& = & \prod_{t=1}^{T}f\!\left(y_{x,t} | E_{x,t},{\boldsymbol{\vartheta}}_x,\mu_{x,t}\left(\kappa_{x,t}\right)\right)f\!\left(\mu_{x,t}|{\boldsymbol{\vartheta}}_x,\kappa_{x,t}\right)f\!\left(\kappa_{x,t}|{\boldsymbol{\vartheta}}_x\right) \\[3pt] & =& \prod_{t=1}^{T}\frac{E_{x,t}\exp\!\left( \mathbb{B}_x\left(t;\ 0,l,{\boldsymbol{\tau}}\right)+\kappa_{x,t}+ \left(1-2u_xB+B^{2}\right)^{-d_x} \varepsilon _{x,t}\right)\left(1-\nu_x\right) }{\Gamma\! \left(y_{x,t}+1\right)\exp\! \left(y_{x,t}\nu_x\right) } \\[3pt] & \times& \frac{\left[E_{x,t}{\exp}\!\left( \mathbb{B}_x\left(t;\ 0,l,{\boldsymbol{\tau}}\!\right){+}\kappa_{x,t}{+} \left(1{-}2u_xB+B^{2}\right)^{-d_x}\! \varepsilon _{x,t}\right)\!\left(1{-}\nu_x\right) +y_{x,t}\nu_x\right]^{y_{x,t}-1}}{\exp\! \left\{E_{x,t}\exp\!\left[ \mathbb{B}_x\left(t;\ 0,l,{\boldsymbol{\tau}}\right)+\kappa_{x,t}+ \left(1-2u_xB+B^{2}\right)^{-d_x} \varepsilon _{x,t}\right]\left(1-\nu_x\right) \right\}} \\[3pt] & \times& \frac{1}{\sqrt{2 \pi}\sigma^2_{x,\varepsilon}}\exp\!\left ( \frac{\varepsilon _{x,t}^{2}}{2\sigma^2_{x,\varepsilon}} \right )\frac{1}{\sqrt{2 \pi}\sigma^2_{x,\varepsilon^{\kappa}}}\exp\!\left ( \frac{\left( \kappa_{x,t}-\eta_x\kappa_{x,t-1} \right)^{2}}{2\sigma^2_{x,\varepsilon^{\kappa}}} \right )\end{eqnarray*}

where ${\boldsymbol{\mu}}_{x,1:t}({\boldsymbol{\kappa}}_{x,1:t})=(\mu_{x,1}(\kappa_{x,1}), \mu_{x,2}(\kappa_{x,2}), \cdots, \mu_{x,t}(\kappa_{x,t}))$. The following priors $\pi(\boldsymbol{\vartheta}_x)$ are adopted in this paper:

\begin{eqnarray*} \nu_x \sim {\textrm{U}}({-}1,1), u_x\sim {\textrm{U}}({-}1,1), d_x\sim {\textrm{U}}(0,1/2), \alpha_x, a_{x,i}^{(0)}, \cdots, a_{x,i}^{(r)} \sim {\textrm{N}}({-}5,1), \\ \eta_x \sim {\textrm{U}}({-}1,1), \sigma^2_{x,\varepsilon} \sim \mbox{Gamma}(1,1)\ \mbox{and}\ \sigma^2_{x,\varepsilon^{\kappa}} \sim \mbox{Gamma}(1,1)\end{eqnarray*}

where U$(a_{\theta},b_{\theta})$ denotes the uniform priors on the range $(a_{\theta}, b_{\theta})$ for parameter $\theta$, which represents the shape parameter $\nu_x$, the oscillatory parameter $u_x$, the long memory parameter $d_x$ and the AR parameter $\eta_x$ in period effect because these parameters are uninformative and the interval for these uniform priors are the bounds for these parameters. In time series and regression settings, it is fairly common to use a normal distribution for coefficients. The model coefficients $\alpha_x, a_{x,i}^{(0)}, \cdots, a_{x,i}^{(r)}$ follow a normal distribution with mean equal to $-5$ and variance equal to 1. The setting for the mean was informed by the empirical study of the log intersect of real-life mortality data sets. Moreover, the variance setting of 1 allows for great flexibility with regard to this mean specification and makes the prior relatively uninformative when viewed on the log scale. Gamma(a, b) denoting the gamma prior with shape and scale parameters a and b is an adequate choice for positive variances $\sigma^2_{x,\varepsilon}$ and $\sigma^2_{x,\varepsilon^{\kappa}}$ (Gelman, Reference Gelman2006). To implement the proposed models efficiently, we choose the Bayesian R package Rstan, which utilises the stan program within R developed in the C++ language. The HMC sampler (Duane et al., Reference Duane, Kennedy, Pendleton and Roweth1987; Neal, Reference Neal1994) as a class of MCMC sampling methods (Metropolis et al., Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953) is adopted in Rstan. For some complicated models with many parameters, the HMC sampler converges faster than the conventional samplers such as random-walk Metropolis (Metropolis et al., Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953) and Gibbs sampler (Geman & Geman, Reference Geman and Geman1984). In order to closely monitor the dependence, precision and convergence of posterior samples, three measures are reported in Rstan. The first measure is the number of effective samples which indicates the effective posterior sample size after allowing for the dependence within a Monte Carlo sample. The second measure is the Monte Carlo standard error (MCSE)

\begin{equation*}\mbox{MCSE}=\frac{\mbox{posterior standard deviation}}{\sqrt {\mbox{number of effective samples}}}\end{equation*}

which reports the error of estimation for the posterior mean. To monitor the convergence for $k > 2$ chains of length 2n each, Gelman & Rubin (Reference Gelman and Rubin1992) proposed $\widehat{R}$. If $\widehat{R}$ is close to 1, the parameter ${\boldsymbol{\vartheta}}_x$ has converged.

Bayesian forecasting is another application of Bayesian inference. The predictive distribution can be derived from the Bayesian filtering equations $f(y_{x,t}|E_{x,t},\mathcal{F}_{x,1:t-1})$, which compute the marginal posterior distribution or filtering distribution of the state $y_{x,t}$ at each time $t-1$ given the history of the measurements up to the time t. The Bayesian filter for computing the predicted distribution at the time T is given by the following Bayesian filtering equation (9). Suppose we want m-step forecasts ${\textbf{\textit{y}}}_{x,T+1:T+m}$ with the observed data filtration $\mathcal{F}_{x,1:T}$ for age group x defined in equation (4) and observed exposures ${\textbf{\textit{E}}}_{x,1:T+m}$. We note that ${\textbf{\textit{E}}}_{x,T+1:T+m}$ can not be directly observed in general. However, there are many methods to estimate ${\textbf{\textit{E}}}_{x,T+1:T+m}$. For example, the rational expectation model defined as a piecewise constant extension $E_{x,T+s}=\mathbb{E}({\textbf{\textit{E}}}_{x,1:T+s-1})$ for $s=1,2,\cdots,m$ can be applied to estimate future exposures ${\textbf{\textit{E}}}_{x,T+1:T+m}$. Alternatively, we can assume a stochastic model $f({\textbf{\textit{E}}}_{x,T+1:T+m}|{\textbf{\textit{E}}}_{x,1:T},{\boldsymbol{\vartheta}}_x,\mathcal{F}_{x,1:T})$ for future exposures. The posterior predictive distribution for ${\textbf{\textit{y}}}_{x,T+1:T+m}$ is defined as

(9) \begin{eqnarray}& & f({\textbf{\textit{y}}}_{x,T+1:T+m}|{\textbf{\textit{E}}}_{x,1:T+m},\mathcal{F}_{x,1:T}) \nonumber \\[3pt] &=& \int \! {\dots}\! \int \prod_{s=1}^{m} \left[f( y_{x,T+s}|{\textbf{\textit{E}}}_{x,1:T+s},{\boldsymbol{\mu}}_{x,1:T+s}({\boldsymbol{\kappa}}_{x,1:T+s}), {\boldsymbol{\vartheta}}_x, \mathcal{F}_{x,1:T})\right.\nonumber\\[3pt] &&\times\ f( \mu_{x,T+s}( \kappa_{x,T+s})|\mu_{x,T+s-1}( \kappa_{x,T+s-1}), {\boldsymbol{\vartheta}}_x, \mathcal{F}_{x,1:T}) \nonumber \\[3pt] & &\left. \times\ f( \kappa_{x,T+s}| \kappa_{x,T+s-1}, {\boldsymbol{\vartheta}}_x , \mathcal{F}_{x,1:T})\right] f({\boldsymbol{\mu}}_{x,1:T}({\boldsymbol{\kappa}}_{x,1:T})| {\boldsymbol{\vartheta}}_x , \mathcal{F}_{x,1:T}) f({\boldsymbol{\kappa}}_{x,1:T}| {\boldsymbol{\vartheta}}_x , \mathcal{F}_{x,1:T}) \nonumber \\[3pt] & &\times\ \pi({\boldsymbol{\vartheta}}_x|\mathcal{F}_{x,1:T}) d {\boldsymbol{\mu}}_{x,1:T+s}({\boldsymbol{\kappa}}_{x,1:T+s})d {\boldsymbol{\kappa}}_{x,1:T+s}d{\boldsymbol{\vartheta}}_x \end{eqnarray}

and this integral can be approximated by the Monte Carlo estimator, constructed from posterior samples according to

(10) \begin{eqnarray}& & \hat{f}\left( y_{x,T+s}|{\textbf{\textit{E}}}_{x,1:T+s},\mathcal{F}_{x,1:T}\right) = \frac{1}{L}\sum_{l=1}^{L}f\left(y_{x,T+s}|{\textbf{\textit{E}}}_{x,1:T+s},{\boldsymbol{\mu}}_{x,1:T+s}^{(l)}\left({\boldsymbol{\kappa}}_{x,1:T+s}^{(l)}\right), {\boldsymbol{\vartheta}}_{x}^{(l)},\mathcal{F}_{x,1:T}\right) \end{eqnarray}

where L is the number of iterations after burn-in in each MCMC sampler run given the current window of information and ${\boldsymbol{\mu}}_{x,1:T+s}^{(l)}({\boldsymbol{\kappa}}_{x,1:T+s}^{(l)})$, ${\boldsymbol{\kappa}}_{x,1:T+s}^{(l)}$ and ${\boldsymbol{\vartheta}}^{(l)}_{x}$ are the lth draw in the posterior sample of ${\boldsymbol{\mu}}_{x,1:T+s}({\boldsymbol{\kappa}}_{x,1:T+s})$, ${\boldsymbol{\kappa}}_{x,1:T+s}$ and ${\boldsymbol{\vartheta}}_x$, respectively.

3.2 Model selection and forecast performance

The performance of each model is evaluated through a popular Bayesian model selection criteria called Deviance Information Criterion (DIC) (Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and Van Der Linde2014). As a generalisation of Akaike’s Information Criterion (AIC), DIC can deal with models containing informative priors, such as hierarchical models. As the priors can effectively restrict the freedom of model parameters, the number of parameters as required in the calculation of AIC is generally unclear. DIC overcomes such problems by providing an estimate for the effective number of parameters. The DIC can be calculated using the following equation:

\begin{equation*}\text{DIC}=\bar{D}+p_D=2\bar{D}-D(\bar{{\boldsymbol{\vartheta}}}_x) \end{equation*}

where the deviance is defined as $D(\boldsymbol{\vartheta_x} )=-2\ln(f({\textbf{\textit{y}}}_x|{\boldsymbol{\vartheta}}_x))$, $\bar{D}=\mathbb E_{\vartheta_x|y_x}[{-}2\ln(f({\textbf{\textit{y}}}_x|{\boldsymbol{\vartheta}}_x))]$ measures the model fit and the estimated number of parameters $p_{D}=\bar{D}-D(\bar{\boldsymbol{\vartheta}}_x)$ measures model complexity (Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and Van Der Linde2002).

Consider the m-step ahead forecasts $\hat{y}_{x,t}$ given by the posterior mean or median and the observations $y_{x,t}$ with T time point and g age groups, the forecast performance can be evaluated by adopting three types of measures, namely residuals $r_{x,t}=y_{x,t}-\hat{y}_{x,t}$, percentage errors $p_{x,t}=\frac{r_{x,t}}{y_{x,t}}\times 100$ and scaled errors $\epsilon_{x,t}$ defined in equation (13). Based on $r_{x,t}$ and $p_{x,t}$, three popular criteria, namely mean absolute error (MAE), root mean squared error (RMSE) and mean absolute percentage error (MAPE), are defined, respectively, below

(11) \begin{align}\mbox{MAE} & =\frac{1}{g}\sum_{x=1}^{g}\left[\frac{1}{m}\sum_{t=1}^{m}|r_{x,T+t}|\right], \mbox{RMSE}=\sqrt{\frac{1}{g}\sum_{x=1}^{g}\left[\frac{1}{m}\sum_{t=1}^{m}r_{x,T+t}^2\right]} \nonumber \\[5pt] \mbox{and}\quad \mbox{MAPE} & =\frac{1}{g}\sum_{x=1}^{g}\left[\frac{1}{m}\sum_{t=1}^{m}|p_{x,T+t}|\right] \end{align}

However $r_{x,t}$ are scale-dependent making comparison difficult. Although $p_{x,t}$ are scale-free, they are sensitive to observations close to zero. Hence, the fourth criterion we adopt is mean absolute scaled error (MASE), which is defined as

(12) \begin{equation}\mbox{MASE}=\frac{1}{g}\sum_{x=1}^{g}\left[\frac{1}{m}\sum_{t=1}^{m}|\epsilon_{x,T+t}|\right] \end{equation}

making use of the scaled errors

(13) \begin{equation}\epsilon_{x,T+t}=\frac{r_{x,T+t}}{\frac{1}{m-1}\sum\limits_{t=2}^{m}|y_{x,T+t}-y_{x,T+t-1}|} \end{equation}

proposed by Hyndman & Koehler (Reference Hyndman and Koehler2006). Furthermore, similar approach can also be applied to evaluate estimated results $\hat{\mu}_{x,t}$ calculated by the posterior mean or median. Hence, the residuals $r_{x,t}^s=\mu_{x,t}-\hat{\mu}_{x,t}$, percentage errors $p_{x,t}^s=\frac{r_{x,t}^s}{\mu_{x,t}}\times 100$ and scaled errors

\begin{equation*}\epsilon_{x,t}^s=\frac{r_{x,t}^s}{\frac{1}{m-1}\sum\limits_{t=2}^{m}|\mu_{x,t}-\mu_{x,t-1}|}\end{equation*}

can be used to construct similar criteria for the $\mu$ estimator, namely the MAE, RMSE, MAPE and MASE by using the same formulas in equations (11)–(12).

4.1 Life expectancy calculation

4.1.1 Definition of life expectancy

Using the posterior mean of the forecast mortality rates $ \mu_{x,t}$, we can evaluate the period life expectancy at different ages and construct life tables (Koissi et al., Reference Koissi, Shapiro and Högnôs2006). We consider the same age groups $x \in \{0, 1-4, 5-9, \cdots , 95-99\}$ and let $\tilde{x}$ to represent the initial age of each age group, that is $\tilde x \in \mathcal{S}=\{0, 1, 5, \cdots , 90, 95\}$. Defining further $n_{\tilde x}$ as the length of the interval of age group x (corresponds to $\tilde x$), we have $n_0 = 1, n_1 = 4, n_5 = 5, \cdots , n_{95} = 5$. We then calculate the probability that a person aged $\tilde x$ in year t will die in the next $n_{\tilde x}$ years as

\begin{equation*}q_{\tilde x,t}^{n_{\tilde x}}=\frac{n_{\tilde x}\mu_{x,t}}{1+n_{\tilde x}(1-a(\tilde x,n_{\tilde x}))\mu_{x,t}}\end{equation*}

where $a(\tilde x,n_{\tilde x})$ is the average fraction of the $n_{\tilde x}$ years lived by the people who aged initially at $\tilde x$ will die during that age interval. Using the assumption that deaths are distributed uniformly in the interval, we set $a(\tilde x,n_{\tilde x})=0.5$ for all $\tilde x$. The hypothetical number of people $l_{\tilde x +n_{\tilde x},t}$ alive at age $\tilde x +n_{\tilde x}$ is determined by

\begin{eqnarray*}l_{\tilde x +n_{\tilde x},t}&=&l_{\tilde x ,t}\left(1-q_{\tilde x,t}^{n_{\tilde x}}\right)=l_{0,t}\prod_{ \left\{k \le \tilde x : k \in \mathcal{S} \right\} }\left(1-q_{k,t}^{n_{k}}\right) \nonumber\end{eqnarray*}

where $l_{0,t}=100,000$. We can then calculate the number of deaths

\begin{equation*}d_{\tilde x,t}^{n_{\tilde x}}=l_{\tilde x ,t}-l_{\tilde x +n_{\tilde x},t}=l_{\tilde x ,t}q_{\tilde x,t}^{n_{\tilde x}} \nonumber\end{equation*}

and the person-years lived

\begin{eqnarray*}L_{\tilde x,t}^{n_{\tilde x}}&=&n_{\tilde x}\left(l_{\tilde x +n_{\tilde x},t}+a(\tilde x,n_{\tilde x})\times d_{\tilde x,t}^{n_{\tilde x}}\right)=n_{\tilde x}l_{\tilde x ,t}\left(1-0.5q_{\tilde x,t}^{n_{\tilde x}}\right) \nonumber\end{eqnarray*}

The total future lifetime of the $l_{\tilde x ,t}$ people who attain age $\tilde x$ is $T_{\tilde x ,t}=\sum_{ \left\{j > \tilde x : j \in \mathcal{S} \right\} }L_{j,t}^{n_{\tilde x}}$. Hence, a period life expectancy at age $\tilde x$ can be derived by

(14) \begin{equation}e_{\tilde x ,t}=\frac{T_{\tilde x ,t}}{l_{\tilde x ,t}} \end{equation}

4.1.2 Bayesian approach to life expectancy calculation

According to the definition of life expectancy $e_{\tilde x ,t}$ in equation (14), we can now obtain a relationship directly between $e_{\tilde x ,t}$ and mortality rates $\mu_{x,t}$ as

(15) \begin{equation}e_{\tilde x ,t}=\frac{\sum_{ \left\{j > \tilde x : j \in \mathcal{S} \right\} }n_{j}\left(1-0.5\frac{n_{j}\mu_{j,t}}{1+0.5n_{j}\mu_{j,t}}\right)l_{0,t}\prod_{ \left\{k \le j-n_j : k\in \mathcal{S} \right\} }\left(1-\frac{n_{k}\mu_{k,t}}{1+0.5n_{k}\mu_{k,t}}\right)}{l_{0,t}\prod_{ \left\{ k \le \tilde x-n_{\tilde x } : k \in \mathcal{S} \right\} }\left(1-\frac{n_{k}\mu_{k,t}}{1+0.5n_{k}\mu_{k,t}}\right)} \end{equation}

In order to calculate life expectancy, mortality rates $\mu_{x,t} $ for all age groups are needed. For calculation convenience, we denote ${\boldsymbol{\mu}}_t \in \mathbb{R}^g$ to be a vector of all mortality rates at year t, ${\textbf{\textit{e}}}_t \in \mathbb{R}^g$ the corresponding vector for life expectancy and $\mathcal{F}_{1:T} $ to be the filtration for all age groups. Therefore, in equation (15), this shows that we can represent the random vector, ${\textbf{\textit{e}}}_{t} \in \mathbb{R}^g$, by a vector value transformation of the random mortality rates ${\boldsymbol{\mu}}_{t}$, which we denote by ${\textbf{\textit{e}}}_{t} =\varrho\!\left( {\boldsymbol{\mu}}_{t} \right)$ with $\varrho\,:\, \mathbb{R}^g \rightarrow \mathbb{R}^g $. Each element of the vector function $\varrho$ corresponds to a mapping for age group x given by equation (15). We are interested in obtaining the posterior predictive distribution of ${\textbf{\textit{e}}}_{T+1\,:\,T+m}$ given $\mathcal{F}_{1\,:\,T}$. We first need to obtain the posterior predictive distribution of ${\boldsymbol{\mu}}_{T+1\,:\,T+m}$ given $\mathcal{F}_{1\,:\,T}$, which is defined as follows:

(16) \begin{eqnarray} && \hspace*{-25pt}f({\boldsymbol{\mu}}_{T+1\,:\,T+m}|\mathcal{F}_{1\,:\,T}) \nonumber \\ & =& \int \!{\dots}\! \int \prod_{s=1}^{m} \left[f({\boldsymbol{\mu}}_{T+s}({\boldsymbol{\kappa}}_{T+s})|{\boldsymbol{\mu}}_{T+s-1}({\boldsymbol{\kappa}}_{T+s-1}),{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T}) f({\boldsymbol{\kappa}}_{T+s}|{\boldsymbol{\kappa}}_{T+s-1},{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})\right] \nonumber \\ & & \times\ f({\boldsymbol{\mu}}_{1\,:\,T}({\boldsymbol{\kappa}}_{1\,:\,T})| {\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})f({\boldsymbol{\kappa}}_{1\,:\,T}| {\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})\pi({\boldsymbol{\vartheta}}|\mathcal{F}_{1\,:\,T}) d {\boldsymbol{\mu}}_{1\,:\,T+m}({\boldsymbol{\kappa}}_{1\,:\,T+m})d {\boldsymbol{\kappa}}_{1\,:\,T+m}d{\boldsymbol{\vartheta}} \end{eqnarray}

Remark 4.1. Equation (16) is an example for autoregressive order 1, which means that ${\boldsymbol{\mu}}_{T+s}$ only depends on ${\boldsymbol{\mu}}_{T+s-1}$.

Our interest is in obtaining the posterior predictive distribution of the life expectancy ${\textbf{\textit{e}}}_{T+1\,:\,T+m}$ given $\mathcal{F}_{1\,:\,T}$. We note the distribution of ${\textbf{\textit{e}}}_t$ conditional on $\mathcal{F}_{1\,:\,T}$ is then obtained via the following generic expression:

\begin{eqnarray*} && f_{{\textbf{\textit{e}}}_t}({\textbf{\textit{e}}}_{t}|\mathcal{F}_{1\,:\,t}) = f_{{\boldsymbol{\mu}}_t}(\varrho ^{-1}({\textbf{\textit{e}}}_{t}) |\mathcal{F}_{1\,:\,t}) |J(\varrho({\boldsymbol{\mu}}_{t}))| \nonumber\end{eqnarray*}

where $f_{{\textbf{\textit{e}}}_t}$ represents the multivariate probability density function (pdf) of life expectancies ${\textbf{\textit{e}}}_t $, and $f_{{\boldsymbol{\mu}}_t}$ is the multivariate pdf of mortality rates ${\boldsymbol{\mu}}_t$. In addition, J is the Jacobian of the transformation $\varrho ({\cdot})$ which maps mortality rates to life expectancy according to equation (15). The distribution summaries for the predictive distribution of life expectancy for each age group, such as the marginal mean and variance, can be estimated in Algorithm 1 step 4. In practice, this procedure can be approximated by the Monte Carlo estimator, constructed from posterior samples according to Algorithm 1. Furthermore, unlike the Bayesian filtering which is based on filtration $\mathcal{F}_{1\,:\,T}$, the Bayesian smoothing uses also the future measurements for computation. The estimator of states $y_{x,t}$ at each time $t<T$ conditional on all the measurements up to time T is $f(y_{x,t}|\mathcal{F}_{x,1\,:\,T})$ (Särkkä, Reference Särkkä2013). The theorem proposed by Kitagawa (Reference Kitagawa1987) provides the marginal posterior distributions for Bayesian smoothing.

Algorithm 1 Compute posterior predictive distribution for life expectancy $f({\textbf{\textit{e}}}_{T+1\,:\,T+m}|\mathcal{F}_{1\,:\,T})$

Theorem 4.1 (The Bayesian smoothing). The Bayesian smoother for computing the smoothed mortality rates distribution $f({\boldsymbol{\mu}}_{t}|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})$ for any $t<T$ is given by the following Bayesian smoothing equation:

\begin{eqnarray*}f({\boldsymbol{\mu}}_{t}|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})&=&f({\boldsymbol{\mu}}_{t}|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,t})\int \left [ \frac{f({\boldsymbol{\mu}}_{t+1\,:\,T}({\boldsymbol{\kappa}}_{t+1\,:\,T})|{\boldsymbol{\mu}}_{t}({\boldsymbol{\kappa}}_{t}),{\boldsymbol{\vartheta}} )f({\boldsymbol{\mu}}_{t+1\,:\,T}({\boldsymbol{\kappa}}_{t+1\,:\,T})|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})}{f({\boldsymbol{\mu}}_{t+1\,:\,T}({\boldsymbol{\kappa}}_{t+1\,:\,T})|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,t})} \right ] \nonumber \\[4pt] & & \times f({\boldsymbol{\kappa}}_{1\,:\,T}|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,T})d{\boldsymbol{\mu}}_{t+1\,:\,T}({\boldsymbol{\kappa}}_{t+1\,:\,T}) d {\boldsymbol{\kappa}}_{1\,:\,T}\end{eqnarray*}

where $f({\boldsymbol{\mu}}_{t+1\,:\,T}({\boldsymbol{\kappa}}_{t+1\,:\,T})|{\boldsymbol{\vartheta}}, \mathcal{F}_{1\,:\,t})$ is the predictive distribution defined in equation (16).

Since the Bayesian smoothing of life expectancy $f(e_{\tilde x,t}|\mathcal{F}_{x,1\,:\,T})$ can be calculated by mortality rates $\mu_{x,t}$, the smoothed posterior distribution of life expectancy $e_{\tilde x ,t}$ is given by Algorithm 1.

Table 3. Data length and abbreviation for country names

4.2 Mortality data

The mortality data sets of 16 countries are downloaded from the HMD, which provides detailed mortality and population data to the public. Table 3 reports the data size for each country together with its abbreviated name.

Each data set contains the number of deaths $y_{x,t}^*$ and the exposure $E_{x,t}^*$ cross-classified by genders and years and they are aggregated into 21 age groups, namely age 0 group, age 1–4 groups, age 5–10 groups, age 11–15 groups, etc. Values of $y_{x,t}^*$ and $E_{x,t}^*$ are scaled according to the following equation (17):

(17) \begin{equation}y_{x,t}=[y_{x,t}^*/c_x] \ \ \mbox{and} \ \ E_{x,t}=E_{x,t}^*/c_x \end{equation}

where $[{\cdot}]$ represents rounding to the nearest integer, and the scaling factor $c_x$ is determined by

\begin{equation*}c_x=10^k \ \ \mbox{if} \ \ \underset{t}{\min} \ y_{x,t}^* \in \left(10^{k+1},10^{k+2}\right)\end{equation*}

such that the range of $y_{x,t}$ is approximately within (0,100) and k is an integer selected to achieve this. This re-scaling of data only changes the unit and preserves the data features. In other words, if the number of deaths $y_{x,t}^*$ and the exposure $E_{x,t}^*$ are scaled by the same amount $c_x$, the features of data will maintain. This can be explained by following equation (18):

(18) \begin{equation}\mathbb E\left(y_{x,t}^*/c_x|E_{x,t}^*\right)\approx \mathbb E\left(\left[y_{x,t}^*/c_x\right]|E_{x,t}^*\right)=E_{x,t}^*\mu_{x,t}/c_x=E_{x,t}\mu_{x,t} \end{equation}

Furthermore, we find comprehensively statistically significant evidence for the presence of long memory features, which are present irrespective of how the data sets are aggregated or disaggregated. To illustrate this, we take Australia female death count as an example. There are different ways to aggregate age groups. For example, taking five age groups for each way, one can form single age $\{0, 1, 2, 3, 4\}$ for type 1 or alternatively $\{0, 1-4, 5-9, 10-14, 15-19\}$ for type 2 or $\{0-10, 11-20, 21-30, 31-40, 41-50\}$ for type 3. However, since some features like weak correlation in a time series may be diminished or removed via aggregation, we need to test if the persistence in the form of long memory that exists in the data will also be reduced or eliminated by aggregation. To do this, we apply Hurst exponent H (Hurst, Reference Hurst1951), a non-parametric measure of long memory, to measure the strength of long memory for the three types of aggregation. Davidson & De Jong (Reference Davidson and De Jong2000) demonstrated the simple relationship $d = H - 0.5$ between d and H, confirming the role of H to measure the strength of long memory in a given time series. Table 4 shows that the estimated H are similar across the three types of aggregation, testifying the prevalence of long memory as a structural feature in aggregation typically.

Table 4. Estimated Hurst exponent for three types of aggregation. Age group differs across different types of aggregation, for example, age group 2 for Type 1 is 1, for Type 2 is 1–4 and for Type 3 is 11–20

Figure 4 Estimated $d_x$ (first row) and $u_x$ (second row) using female (first two rows) and male (next two rows) death counts fitted by model 6 for the five selected countries. The estimated value of $u_x$ for Iceland female in the last age group is 0.08313 and for Iceland male in the last two age groups are $-0.03799$ and $-0.07389$. For US male, the last reading is 0.9875.

5.1 Out-of-sample forecasting

In order to measure the forecast performance of the eight models, the number of death counts ${\textbf{\textit{Y}}}_{x,1\,:\,T}$ by 21 age groups and both genders of the 5 selected countries are each divided into 2 parts, the training ${\textbf{\textit{Y}}}_{x,1\,:\,(T-20)}$ and forecast ${\textbf{\textit{Y}}}_{x,(T-19)\,:\,T}$ for years 1994 to 2014. Similar to Table 5 , Table 6 reports the four criteria that compare the mortality rate forecasts $\hat{{\boldsymbol{\mu}}}_{x,(T-19)\,:\,T}$ with the observed from HMD for the eight models. Model 8 again provides the best forecasts in all except seven cases coming from Iceland in which five cases favour model 3 as Iceland data has a short memory.

To view the improved forecast performance of model 8 compared with model 3, Figure 10 plots the time series of observed ${\boldsymbol{\mu}}_{x,T-19\,:\,T}$ from HMD and the forecasted $\hat{{\boldsymbol{\mu}}}_{x,T-19\,:\,T}$ using model 3 and model 8 with 95% credible interval for the young, middle and senior age groups, respectively. The forecast performance of model 8 with a long memory component is much better than the traditional LC model (model 3) because the mortality rate forecasts are closer to the observed and have sharper credible intervals, which enclose the observed mortality rates. Moreover, the traditional LC model consistently overestimates the mortality rate.

Figure 11 displays the 20-step ahead forecasts of Australia life expectancy using model 8 for both genders and 4 age groups (0, 20–24, 40–44 and 60–64).

The performance of out-of-sample life expectancy forecasts is also reasonable as the forecasted values are close to the life expectancies calculated using observed mortality rates from HMD and these “observed” life expectancies are contained within the credible interval.

A range of additional forecast results for mortality rates and life expectancies for Iceland, Japan, UK and the US is provided in Yan et al., (Reference Yan, Peters and Chan2018b). This includes further analysis of the four model performance criteria for the forecasted number of deaths $\widehat{{\textbf{\textit{Y}}}}_{x,(T-19)\,:\,T}$ using all the eight models.

Table 6. Four selection criteria to measure the accuracy of forecast $\hat{{\boldsymbol{\mu}}}_{x,(T-19)\,:\,T}$ for the five selected countries

Note: the best models amongst models 1-8 are shaded by light grey colour.

Figure 10 Observed mortality rates ${\boldsymbol{\mu}}_{x,T-19\,:\,T}$ (black line) from HMD and forecasted $\hat{{\boldsymbol{\mu}}}_{x,T-19\,:\,T}$ by model 3 (female in purple and male in skyblue) and model 8 (female in red and male in blue) for both female (first row) and male (second row) and three age groups (10–14, 35–39 and 70–74) in Australia.

Figure 11 Observed life expectancy $e_{\tilde x ,t}$ from HMD (black line) and life expectancy calculated using the forecasted mortality rates from model 8 (red in female and blue in male) for female (first row) and male (second row) and four age groups (0, 20–24, 40–44 and 60–64) in Australia.

Table 7. The forecasted Australia life expectancy by model 3 and model 8

5.2 Life expectancy of 16 countries forecasted 20 years ahead

In this analysis, we calculate life expectancy based on the 20 years ahead of forecasts of mortality rates for 16 countries using model 8, the best performing model in the previous 2 analyses and model 3, the LC model for comparison. Table 7 reports the life table forecast for selected age groups 0, 15 – –19, 60 – –69 in Australia using model 3 and model 8. The life expectancy for all 16 countries and all age groups are given in Yan et al., (Reference Yan, Peters and Chan2018b).

The trend of life expectancy forecast is slowly increasing over the forecast period for all age groups, gender and models. Compared with model 3, the life expectancy forecast using model 8 is higher consistently for all age groups, gender and years.