3.1 Compositional data analysis

In the field of demography, Oeppen (Reference Oeppen2008) and Bergeron-Boucher et al. (Reference Bergeron-Boucher, Canudas-Romo, Oeppen and Vaupel2017) have laid the foundations by presenting a modelling and forecasting framework. Following Oeppen (Reference Oeppen2008) and Bergeron-Boucher et al. (Reference Bergeron-Boucher, Canudas-Romo, Oeppen and Vaupel2017), a CoDa method can be summarised as follows:

1. (1) We begin from a data matrix $\textbf{\textit{D}}$ of size $n\times K$ of the life-table death counts $(d_{t,x})$ with n rows representing the number of years and K columns representing the age x. The sum of each row adds up to the life-table radix, such as 100,000. Since we are working with life-table death counts, it is not necessary to have population-at-risk estimates.

2. (2) We compute the geometric mean at each age, given by

   (1) \begin{equation} \alpha_x = \exp^{\frac{1}{n}\sum^n_{t=1}\ln (d_{t,x})}, \qquad x = 1,\dots, K,\quad t=1,\dots,n \label{eqn1} \end{equation}
   For a given year t, we divide $(d_{t,1},\dots,d_{t,K})$ by the corresponding geometric means $(\alpha_1,\dots,\alpha_K)$ , and then standardise all elements so as to sum up to unity. As with compositional data, it is more important to know a relative proportion of each component (i.e. age) than the sum of all components. Although the life-table radix is 100,000 customarily, we model and forecast relative proportions, which are then multiplied by 100,000 to obtain forecast life-table death counts. Via standardisation, this is expressed as
   \begin{equation*} f_{t,x} = \frac{d_{t,x}/\alpha_x}{\sum^K_{x=1}d_{t,x}/\alpha_x} \end{equation*}
   where $f_{t,x}$ denotes de-centred data.

3. (3) Log-ratio transformation. Aitchison (Reference Aitchison1982, Reference Aitchison1986) showed that compositional data are represented in a restricted space where the components can only vary between 0 and a positive constant. Therefore, Aitchison (Reference Aitchison1982, Reference Aitchison1986) proposed a log-ratio transformation to transform the data into a real-valued space. We apply the centred log-ratio transformation, given by

   \begin{equation*} z_{t,x} = \ln\left(\frac{f_{t,x}}{g_t}\right) \end{equation*}
   where $g_t$ denotes the geometric mean over the age at time t, given by
   \begin{equation*} g_t = \exp^{\frac{1}{K}\sum^K_{x=1}\ln (f_{t,x})} \end{equation*}
   The transformed data matrix is denoted as $\textbf{\textit{Z}}$ with elements $z_{t,x}\in R$ , where R denotes real-valued space.

4. (4) Principal component analysis. Principal component analysis is then applied to the matrix $\textbf{\textit{Z}}$ , to obtain the estimated principal components and their scores,

   (2) \begin{equation} z_{t,x} = \sum_{\ell=1}^{\min(n,K)}\widehat{\beta}_{t,\ell}\widehat{\phi}_{\ell,x} = \sum^L_{\ell=1}\widehat{\beta}_{t,\ell}\widehat{\phi}_{\ell,x}+\widehat{\varpi}_{t,x} \label{eqn2} \end{equation}
   where $\widehat{\varpi}_{t,x}$ denotes model residual term for age x in year t, $\{\widehat{\phi}_{1,x},\dots,\widehat{\phi}_{L,x}\}$ represents the first L sets of estimated principal components, $\{\widehat{\beta}_{t,1},\dots,\widehat{\beta}_{t,L}\}$ represents the first L sets of estimated principal component scores for time t and L denotes the number of retained principal components.

5. (5) Forecast of principal component scores. Via a univariate time series forecasting method, such as ETS, we obtain the h-step-ahead forecast of the $\ell^{\text{th}}$ principal component score $\widehat{\beta}_{n+h, \ell}$ , where h denotes forecast horizon (see also Hyndman & Ullah, Reference Hyndman and Ullah2007). We utilise an automatic algorithm developed by Hyndman & Khandakar (Reference Hyndman and Khandakar2008) to determine the optimal ETS model based on the corrected Akaike information criterion (Hurvich & Tsai, Reference Hurvich and Tsai1993). Conditioning on the estimated principal components and observed data, the forecast of $z_{n+h,x}$ can be obtained by

   (3) \begin{equation} \skew2\widehat{z}_{n+h|n, x} = \sum^{L}_{\ell=1}\widehat{\beta}_{n+h|n,\ell}\widehat{\phi}_{\ell,x} \label{eqn3} \end{equation}
   In (3), the principal component scores can be modelled and forecasted by an ETS method. To select the optimal ETS parameter, we use the automatic search algorithm of Hyndman & Khandakar (Reference Hyndman and Khandakar2008) by the smallest corrected Akaike information criterion. In Table 1, we also consider three other forecasting methods, namely autoregressive integrated moving average (ARIMA), random walk with drift (RWD) and random walk without drift (RW).

6. (6) Transform back to the compositional data. We take the inverse centred log-ratio transformation, given by

   \begin{equation*} \skew4\widehat{\textbf{\textit{f}}}_{n+h|n} = \left[\frac{\exp^{\skew2\widehat{z}_{n+h|n, 1}}}{\sum^K_{x=1}\exp^{\skew2\widehat{z}_{n+h|n, x}}}, \frac{\exp^{\skew2\widehat{z}_{n+h|n, 2}}}{\sum^K_{x=1}\exp^{\skew2\widehat{z}_{n+h|n, x}}}, \dots,\frac{\exp^{\skew2\widehat{z}_{n+h|n, K}}}{\sum^K_{x=1}\exp^{\skew2\widehat{z}_{n+h|n, x}}}\right] \end{equation*}
   where $\skew2\widehat{z}_{n+h|n, x}$ denotes the forecasts in (3).

7. (7) Finally, we add back the geometric means, to obtain the forecasts of the life-table death matrix $d_{n+h}$ ,

   \begin{align*} \skew9\widehat{\textbf{\textit{d}}}_{n+h|n} &= \left[\frac{\skew4\widehat{f}_{n+h|n, 1}\times \alpha_1}{\sum^K_{x=1}\skew4\widehat{f}_{n+h|n, x}\times \alpha_x}, \frac{\skew4\widehat{f}_{n+h|n,2}\times \alpha_2}{\sum^K_{x=1}\skew4\widehat{f}_{n+h|n,x}\times \alpha_x},\dots,\frac{\skew4\widehat{f}_{n+h|n,K}\times \alpha_K}{\sum^K_{x=1}\skew4\widehat{f}_{n+h|n,x}\times \alpha_x}\right] \end{align*}
   where $\alpha_x$ is the age-specific geometric mean given in (1).

Table 1. A comparison of the point and interval forecast accuracy, as measured by the overall MAPE and mean interval score, among the CoDa, HU, LC and two naïve RW methods (in italic) using the holdout sample of the Australian female and male data.

Note: Further, we consider four univariate time series forecasting methods for the CoDa and HU methods. In terms of interval forecast accuracy, we compare the finite-sample performance between the proposed bootstrap approach and an existing bootstrap approach of Bergeron-Boucher et al. (Reference Bergeron-Boucher, Canudas-Romo, Oeppen and Vaupel2017). The smallest errors are highlighted in bold.

To determine the number of components L in (2) and (3), we consider a criterion known as cumulative percentage of variance (CPV) – that is, we define the value of L as the minimum number of components that reaches a certain level of the proportion of total variance explained by the L leading components such that

$$L = {\rm{ }}\mathop {{\rm{argmin}}}\limits_{_{L:L \ge 1}} \left\{ {\sum\limits_{\ell = \mathds{1}}^L {{{\hat \lambda }_\ell }} /\sum\limits_{\ell = 1}^n {{{\hat \lambda }_\ell }} {1_{\left\{ {{{\hat \lambda }_\ell } > 0} \right\}}} \ge \delta } \right\}$$

where $\delta = 85\%$ (see also Horváth & Kokoszka, Reference Horváth and Kokoszka2012, 41) and $\mathds{1}_{\{\cdot\}}$ denotes the binary indicator function that excludes possible zero eigenvalues. For the Australian female and male data, the chosen number of components is $L=1$ and $L=2$ , respectively.

We highlight some similarities and differences between Oeppen’s (Reference Oeppen2008) approach and Lee & Carter’s (Reference Lee and Carter1992) approach. From Steps 3–5, Oeppen’s (Reference Oeppen2008) approach uses the Lee & Carter’s (Reference Lee and Carter1992) approach to model and forecast the log-ratio of life-table death counts. The difference is that the Oeppen’s (Reference Oeppen2008) approach works with life-table death counts in a constrained space, whereas the Lee & Carter’s (Reference Lee and Carter1992) approach works with real-valued log mortality rates.

When the number of components $L=1$ , our proposed CoDa method corresponds to the one presented by Oeppen (Reference Oeppen2008) and Bergeron-Boucher et al. (Reference Bergeron-Boucher, Canudas-Romo, Oeppen and Vaupel2017). However, we allow the possibility of using more than one pair of principal component and principal component scores in Steps 4 and 5. As a sensitivity analysis, we also consider setting the number of components to be $L=6$ (see also Hyndman & Booth, Reference Hyndman and Booth2008). From this aspect, our proposal shares some similarity with the Hyndman & Ullah’s (Reference Hyndman and Ullah2007) approach. The difference is that our proposal works with life-table death counts in a constrained space, whereas the Hyndman & Ullah’s (Reference Hyndman and Ullah2007) approach works with real-valued log mortality rates.

3.2 Construction of prediction interval for the CoDa

Prediction intervals are a valuable tool for assessing the probabilistic uncertainty associated with point forecasts. The forecast uncertainty stems from both systematic deviations (e.g. due to parameter or model uncertainty) and random fluctuations (e.g. due to model error term). As was emphasised by Chatfield (Reference Chatfield1993, Reference Chatfield2000), it is essential to provide interval forecasts as well as point forecasts to

1. (1) assess future uncertainty levels;

2. (2) enable different strategies to be planned for a range of possible outcomes indicated by the interval forecasts;

3. (3) compare forecasts from different methods more thoroughly; and

4. (4) explore different scenarios based on various assumptions.

3.2.1 Proposed bootstrap method

We consider two sources of uncertainty: truncation errors in the principal component decomposition and forecast errors in the projected principal component scores. Since principal component scores are regarded as surrogates of the original functional time series, these principal component scores capture the temporal dependence structure inherited in the original functional time series (see also Salish & Gleim, Reference Salish and Gleim2015; Paparoditis, Reference Paparoditis2018; Shang, Reference Shang2018). By adequately bootstrapping the forecast principal component scores, we can generate a set of bootstrapped $\textbf{\textit{Z}}^*$ , conditional on the estimated mean function and estimated functional principal components from the observed Z in (3).

Using a univariate time series model, we can obtain multi-step-ahead forecasts for the principal component scores, $\{\widehat{\beta}_{1,\ell},\dots,\widehat{\beta}_{n,\ell}\}$ for $l=1,\dots,L$ . Let the h-step-ahead forecast errors be given by $\widehat{\vartheta}_{t,h,\ell} = \widehat{\beta}_{t,\ell}-\widehat{\beta}_{t|t-h,\ell}$ for $t=h+1,\dots,n$ . These can then be sampled with replacement to give a bootstrap sample of $\beta_{n+h,\ell}$ :

\begin{equation*} \widehat{\beta\,}_{n+h|n,\ell}^{(b)} = \widehat{\beta}_{n+h|n,\ell}+\widehat{\vartheta}_{\ast, h, \ell}^{(b)}, \qquad b=1,\dots,B \end{equation*}

where $B=1,000$ symbolises the number of bootstrap replications and $\widehat{\vartheta}_{\ast, h, \ell}^{(b)}$ are sampled with replacement from $\{\widehat{\vartheta}_{t, h, \ell}\}$ .

Assuming the first L principal components approximate the data $\textbf{\textit{Z}}$ relatively well, the model residual should contribute nothing but random noise. Consequently, we can bootstrap the model fit errors in (2) by sampling with replacement from the model residual term $\{\widehat{\varpi}_{1,x},\dots,\widehat{\varpi}_{n,x}\}$ for a given age x.

Adding all two components of variability, we obtain B variants for $z_{n+h,x}$ ,

\begin{equation*} \skew2\widehat{z\,\,}_{n+h|n,x}^{(b)} = \sum^L_{\ell=1}\widehat{\beta\,}^{(b)}_{n+h|n,\ell}\widehat{\phi}_{\ell,x}+\widehat{\varpi}_{n+h,x}^{(b)} \end{equation*}

where $\widehat{\beta\,}^{(b)}_{n+h,\ell}$ denotes the forecast of the bootstrapped principal component scores. The construction of the bootstrap samples $\{\skew1\widehat{z\ }^{(1)}_{n+h|n, x}, \dots, \skew1\widehat{z\ }^{(B)}_{n+h|n,x}\}$ is in the same spirit as the construction of the bootstrap samples for the LC method in (5) and for the HU method in Hyndman & Shang (Reference Hyndman and Shang2010).

With the bootstrapped $\skew1\widehat{z}_{n+h|n,x}^{\ (b)}$ , we follow Steps 6 and 7 in section 3.1, in order to obtain the bootstrap forecasts of $d_{n+h,x}$ . At the $100(1-\gamma)\%$ nominal coverage probability, the pointwise prediction intervals are obtained by taking $\gamma/2$ and $1-\gamma/2$ quantiles based on $\{\skew1\widehat{d\,}_{n+h|n,x}^{(1)},\dots,\skew1\widehat{d\,}_{n+h|n,x}^{(B)}\}$ .

3.3 Other forecasting methods

3.3.1 LC method

As a comparison, we revisit Lee & Carter’s (Reference Lee and Carter1992) method. To stabilise the higher variance associated with mortality at advanced old ages, it is necessary to transform the raw data by taking the natural logarithm. For those missing death counts at advanced old ages, we use a simple linear interpolation method to approximate the missing values. We denote by $m_{t,x}$ the observed mortality rate at year t at age x calculated as the number of deaths in the calendar year t at age x, divided by the corresponding mid-year population aged x. The model structure is given by

(4) \begin{equation} \ln (m_{t,x}) = a_x + b_x\kappa_t + \epsilon_{t,x} \label{eqn4} \end{equation}

where $a_x$ denotes the age pattern of the log mortality rates averaged across years; $b_x$ denotes the first principal component reflecting relative change in the log mortality rate at each age; $\kappa_t$ is the first set of principal component scores by year t and measures the general level of the log mortality rates; and $\epsilon_{t,x}$ denotes the residual at year t and age x.

The LC model in (4) is over-parametrised in the sense that the model structure is invariant under the following transformations:

\begin{align*} \{a_x, b_x, \kappa_t\} &\mapsto \{a_x, b_x/c, c\kappa_t\} \\ \{a_x, b_x, \kappa_t\} &\mapsto \{a_x-cb_x, b_x, \kappa_t+c\} \end{align*}

To ensure the model’s identifiability, Lee & Carter (Reference Lee and Carter1992) imposed two constraints, given as

\begin{equation*} \sum^n_{t=1}\kappa_t=0, \qquad \sum^{x_p}_{x=x_1}b_x=1 \end{equation*}

where n is the number of years, and p is the number of ages in the observed data set.

Also, the LC method adjusts $\kappa_t$ by refitting to the total number of deaths. This adjustment gives more weight to high rates, thus roughly counterbalancing the effect of using a log transformation of the mortality rates. The adjusted $k_{t}$ is then extrapolated using RW models. Lee & Carter (Reference Lee and Carter1992) used an RWD model, which can be expressed as

\begin{equation*} \kappa_{t}=\kappa_{t-1}+d+e_{t} \end{equation*}

where d is known as the drift parameter and measures the average annual change in the series, and $e_{t}$ is an uncorrelated error. It is notable that the RWD provides satisfactory results in many cases (see e.g. Tulkapurkar et al., Reference Tulkapurkar, Li and Boe2000; Lee & Miller, Reference Lee and Miller2001). From this forecast of the principal component scores, the forecast age-specific log mortality rates are obtained using the estimated age effects $a_x$ and $b_x$ in (4).

Note that the LC model can also be formulated within a Generalised Linear Model framework with a generalised error distribution (see e.g. Tabeau, Reference Tabeau, Tabeau, Van Den Berg Jeths and Heathcote2001; Brouhns et al., Reference Brouhns, Denuit and Vermunt2002; Renshaw & Haberman, Reference Renshaw and Haberman2003). In this setting, the LC model parameters can be estimated by maximum likelihood methods based on the choice of the error distribution. Thus, in line with traditional actuarial practice, this approach assumes that the age- and period-specific number of deaths are independent realisations from a Poisson distribution. In the case of Gaussian error, the estimation method based on either the singular value decomposition or maximum likelihood method leads to the same parameter estimates.

Two sources of uncertainty are considered in the LC model: errors in the parameter estimation of the LC model and forecast errors in the forecasted principal component scores. Using a univariate time series model, we can obtain multi-step-ahead forecasts for the principal component scores, $\{\widehat{\kappa}_1,\dots,\widehat{\kappa}_n\}$ . Let the h-step-ahead forecast errors be given by $\nu_{t,h}=\widehat{\kappa}_t-\widehat{\kappa}_{t|t-h}$ for $t=h+1,\dots,n$ . These can then be sampled with replacement to give a bootstrap sample of $\kappa_{n+h}$ :

\begin{equation*} \widehat{\kappa}_{n+h|n}^{(b)} = \widehat{\kappa}_{n+h|n}+\widehat{\nu}_{\ast,h}^{(b)},\qquad b=1,\dots,B \end{equation*}

where $\widehat{\nu}_{\ast,h}^{(b)}$ are sampled with replacement from $\{\widehat{\nu}_{t,h}\}$ .

Assuming the first principal component approximates the data $\ln \textbf{\textit{m}}$ relatively well, the model residual should contribute nothing but random noise. Consequently, we can bootstrap the model fit error in (4) by sampling with replacement from the model residual term $\{\,\widehat{e}_{1,x},\dots,\widehat{e}_{n,x}\}$ for a given age x.

Adding the two components of variability, we obtain B variants of $\ln m_{n+h,x}$ :

(5) \begin{equation} \ln \widehat{m}_{n+h|n,x}^{(b)} = \widehat{\kappa}_{n+h|n}^{\,(b)}\widehat{b}_x + \widehat{e}_{n+h,x}^{\,(b)} \label{eqn5} \end{equation}

where $\widehat{\kappa}_{n+h|n}^{(b)}$ denotes the forecast of the bootstrapped principal component scores.

Since the LC model forecasts age-specific central mortality rate, we convert forecast mortality rate to the probability of dying via a simple approximation. The formula is given as

\begin{equation*} \widehat{q}_{n+h|n,x} = 1 - \exp^{-\widehat{m}_{n+h|n,x}} \end{equation*}

where $\widehat{q}_{n+h|n,x}$ denotes the probability of dying in age x and year $n+h$ , and $\widehat{m}_{n+h|n,x}$ denotes the forecast of central mortality rate in age x and year $n+h$ . With an initial life-table death count of 100,000, the forecast of $\skew1\widehat{d}_{n+h|n,x}$ can be obtained from $\widehat{q}_{n+h|n,x}$ .

3.3.3 Random walk with and without drift

Given the findings of linear life expectancy (White, Reference White2002; Oeppen & Vaupel, Reference Oeppen and Vaupel2002) and the debate about its continuation (Bengtsson, Reference Bengtsson, Bengtsson and Keilman2003; Lee, Reference Lee2003), it is pertinent to compare the forecast accuracy of the CoDa method with a linear extrapolation method (see Alho & Spencer, Reference Alho and Spencer2005, 274–276 for an introduction). Using the centred log-ratio transformation, the linear extrapolation of $z_{t, x}$ in (2) is achieved by applying the RW and RWD models for each age:

(6) \begin{equation} z_{t+1,x} =\zeta + z_{t,x} + e_{t+1,x} \label{eqn6} \end{equation}

where $z_{t,x}$ represents the centred log ratio transformed data at age $x=0, 1, \dots, p$ in year $t=1,2,\dots,n-1$ , $\zeta$ denotes the drift term and $e_{t+1,x}$ denotes model error term. The h-step-ahead point and interval forecasts are given by

(7) \begin{align} \skew2\widehat{z}_{n+h|n,x} &= \text{E}\left[z_{n+h,x}|z_{1,x},\dots,z_{n,x}\right] = \zeta h + z_{n,x} \notag\\ \text{var}(\skew2\widehat{z}_{n+h|n,x}) &= \text{var}[z_{n+h,x}|z_{1,x},\dots,z_{n,x}] = \text{var}(z_{n,x}) + \text{var}(e_{n+h,x}) \label{eqn7}\end{align}

Based on (6) and (7), the RW can be obtained by omitting the drift term. Computationally, the forecasts of the RW and RWD are obtained by the rwf function in the forecast package (Hyndman et al., Reference Hyndman, Athanasopoulos, Bergmeir, Caceres, Chhay, O’Hara-Wild, Petropoulos, Razbash, Wang, Yasmeen, Team, Ihaka, Reid, Shaub, Tang and Zhou2019) in R (R Core Team, 2019). Then, by way of the back-transformation of the centred log ratio, naïve forecasts of age-specific life-table death counts can be obtained.

5.1 Forecast error criteria

An expanding window analysis of a time-series model is commonly used to assess model and parameter stability over time, and prediction accuracy. The expanding window analysis determines the constancy of a model’s parameter by computing parameter estimates and their resultant forecasts over an expanding window of a fixed size through the sample (for details, Zivot & Wang, Reference Zivot and Wang2006: 313–314). Using the first 74 observations from 1921 to 1994 in the Australian female and male age-specific life-table death counts, we produce 1- to 20-step-ahead forecasts. Through an expanding-window approach, we re-estimate the parameters in the time series forecasting models using the first 75 observations from 1921 to 1995. Forecasts from the estimated models are then produced for 1- to 19-step-ahead points and intervals. We iterate this process by increasing the sample size by 1 year until reaching the end of the data period in 2014. This process produces 20 one-step-ahead forecasts, 19 two-step-ahead forecasts, …, and one 20-step-ahead forecast. We compare these forecasts with the holdout samples to determine the out-of-sample point forecast accuracy.

To evaluate the point forecast accuracy, we consider the MAPE that measures how close the forecasts are in comparison to the actual values of the variable being forecast, regardless of the direction of forecast errors. The error measure can be written as

\begin{equation*} \text{MAPE}(h) = \frac{1}{111\times (21-h)}\sum^{20}_{\varsigma=h}\sum^{111}_{x=1}\left|\frac{d_{n+\varsigma, x}-\skew4\widehat{d}_{n+\varsigma,x}}{d_{n+\varsigma, x}}\right|\times 100 \end{equation*}

where $d_{n+\varsigma,x}$ denotes the actual holdout sample for the $x^{\text{th}}$ age and $\varsigma^{\text{th}}$ forecasting year, while $\skew4\widehat{d}_{n+\varsigma,x}$ denotes the point forecasts for the holdout sample.

To evaluate and compare the interval forecast accuracy, we consider the interval score of Gneiting & Raftery (Reference Gneiting and Raftery2007). For each year in the forecasting period, the h-step-ahead prediction intervals are calculated at the $100(1-\gamma)\%$ nominal coverage probability. We consider the common case of the symmetric $100(1-\gamma)\%$ prediction intervals, with lower and upper bounds that are predictive quantiles at $\gamma/2$ and $1-\gamma/2$ , denoted by $\skew4\widehat{d}^l_{n+\varsigma, x}$ and $\skew1\widehat{d}^u_{n+\varsigma, x}$ . As defined by Gneiting & Raftery (Reference Gneiting and Raftery2007), a scoring rule for the interval forecasts at time point $d_{n+\varsigma,x}$ is

\begin{array}{*{20}{c}} {{S_{\gamma ,\varsigma }}\left[ {\hat d_{n + \varsigma ,x}^l,\hat d_{n + \varsigma ,x}^u,{d_{n + \varsigma ,x}}} \right] = \left( {\hat d_{n + \varsigma ,x}^u - \hat d_{n + \varsigma ,x}^l} \right)} & { + \frac{2}{\gamma }\left( {\hat d_{n + \varsigma ,x}^l - {d_{n + \varsigma ,x}}} \right)\mathds{1}\left\{ {{d_{n + \varsigma ,x}} < \hat d_{n + \varsigma ,x}^l} \right\}} \\ {} & { + \frac{2}{\gamma }\left( {{d_{n + \varsigma ,x}} - \hat d_{n + \varsigma ,x}^u} \right)\mathds{1}\left\{ {{d_{n + \varsigma ,x}} > \hat d_{n + \varsigma ,x}^u} \right\}} \\ \end{array}

where $\mathds{1}\{\cdot\}$ represents the binary indicator function and $\gamma$ denotes the level of significance, customarily $\gamma=0.2$ or $\gamma=0.05$ . The interval score rewards a narrow prediction interval, if and only if the true observation lies within the prediction interval. The optimal interval score is achieved when $d_{n+\varsigma, x}$ lies between $\widehat{d}_{n+\varsigma,x}^l$ and $\widehat{d}_{n+\varsigma, x}^u$ , and the distance between $\widehat{d}_{n+\varsigma, x}^l$ and $\widehat{d}_{n+\varsigma, x}^u$ is minimal.

For different ages and years in the forecasting period, the mean interval score is defined by

\begin{equation*} \overline{S}_{\gamma}(h) = \frac{1}{111\times (21-h)} \sum_{\varsigma = h}^{20} \sum_{x=1}^{111} S_{\gamma,\varsigma}\left[\widehat{d}_{n+\varsigma,x}^{l}, \widehat{d}_{n+\varsigma,x}^{u}; d_{n+\varsigma,x}\right] \end{equation*}

where $S_{\gamma,\varsigma}\left[\widehat{d}_{n+\varsigma,x}^{l}, \widehat{d}_{n+\varsigma,x}^{u}; d_{n+\varsigma,x}\right]$ denotes the interval score at the $\varsigma^{\text{th}}$ curve in the forecasting year.

5.2 Forecast results

Using the expanding window approach, we compare the point and interval forecast accuracies among the CoDa, HU, LC and two naïve RW methods based on the MAPE and mean interval score. Also, we consider forecasting each set of the estimated principal component scores by the ARIMA, ETS, RWD and RW. The overall point and interval forecast error results are presented in Table 1, where we average over the 20 forecast horizons. From Table 1, we find that the CoDa method with the ETS forecasting method generally performs the best among a range of methods considered. Among the four univariate time series forecasting methods, the ETS method produces the smallest errors and thus is recommended to be used in practice. Also, setting $L=6$ , including additional principal component decomposition pairs, generally provides more accurate point and interval forecast accuracies than setting $L=1$ or $L=2$ for the Australian female and male data.

The CoDa method tends to provide less bias forecasts than the LC method as it allows the rate of mortality improvements to change over time (see Bergeron-Boucher et al., Reference Bergeron-Boucher, Canudas-Romo, Oeppen and Vaupel2017). With respect to the positivity and summability constraints, the CoDa method can adopt temporal changes of the age distribution of the life-table death counts over the years. The unsatisfactory performance of the HU and LC methods may be because of the approximation of converting the forecasts of central mortality rate to the probability of dying at higher ages. The inferior performance of the two naïve RW methods may be because of their slow responses to adapt to the change of age distribution of death counts. Between the two naïve RW methods, we found that the RW is preferable for producing point forecasts, while the RWD is preferable for producing interval forecasts.

While Table 1 presents the average over 20 forecast horizons, we show the 1-step-ahead to 20-step-ahead point and interval forecast errors in Figure 3. Since it is advantageous to set $L=6$ , we report the CoDa method with the ETS forecasting method and $L=6$ in Figure 3. The difference in forecast accuracy between the CoDa and the other methods is widening over the forecast horizon. We suspect that in the relatively longer forecast horizon, the errors associated with all the methods become larger, but they are relatively smaller for the CoDa method with the ETS forecasting method and $L=6$ .

Figure 3. A comparison of the point and interval forecast accuracy, as measured by the MAPE and mean interval score, among the CoDa, HU, LC and two naïve RW methods using the holdout sample of the Australian female and male data. In the CoDa method, we use the ETS forecasting method with the number of principal components $L=6$ .