2.1. The Hyndman-Ullah (HU) model

The HU model is a generalized extension of the LC model, employing the functional data paradigm (Ramsay and Silverman, Reference Ramsay and Silverman2005) and nonparametric smoothing to improve the forecasting of mortality rates. This model is distinguished by its use of multiple principal components, capturing a broader variation within the data, and applying more advanced time series forecasting methods than the random walk with drift model typically used in the LC model (Hyndman and Ullah, Reference Hyndman and Ullah2007).

At its core, the HU model assumes that log mortality rates $y_t(x)$ observed at discrete ages x and years t have an underlying smooth function $ f_t(x) $ , observed with error:

\begin{equation*} y_t(x) = f_t(x) + \sigma_t(x)\epsilon_t, \end{equation*}

where $\epsilon_t$ is the standard normal error term and $\sigma_t(x)$ denotes the observational standard deviation which is allowed to vary with age. This lets the HU model handle data heterogeneity and ensures demographic consistency by applying monocity constraints when using penalized regression splines to smooth the mortality curves (Hyndman and Ullah, Reference Hyndman and Ullah2007).

FPCA is used to decompose the smoothed mortality rates into principal components. By applying FPCA to the centered smooth mortality would yield:

\begin{equation*} f_t(x) = \mu(x) + \sum_{k=1}^K \beta_{t,k} \phi_k(x) + e_t(x). \end{equation*}

Here, $\mu(x)$ represents the mean function, $\phi_k(x)$ are orthonormal basis functions, $\beta_{t,k}$ are time-dependent coefficients, and $e_t(x)$ is the error function. Forecasts of mortality rates are done through the time-dependent coefficients, $\beta_{t,k}$ which are fitted into univariate ARIMA models and extrapolated h-steps-ahead.

A critical consideration in applying FPCA is its sensitivity to outliers. Extreme values can disproportionately influence the principal components derived from the data. This sensitivity may lead to biased or misleading results, especially when modeling phenomena such as mortality rates where spikes might occur due to extraordinary events like epidemics or wars. To address these issues, various adaptations of the HU model have been put forth. These include the robust HU (HUrob) model with robust FPCA (Hyndman and Ullah, Reference Hyndman and Ullah2007), which seeks to minimize the impact of outliers by using the $L_1$ -median as the location measure, paired with the reflection-based algorithm for principal components analysis (Hubert et al., Reference Hubert, Rousseeuw and Verboven2002) for a robust set of principal components; the weighted HU (wHU) model with weighted FPCA (Hyndman and Shang, Reference Hyndman and Shang2009), which adjusts the influence of data points by placing more weight on recent data with the measure of location being a weighted average using geometrically decreasing weights. These modifications aim to enhance the stability and reliability of the model under conditions where data anomalies are present. Having the same objective in mind to make the model less sensitive to outliers, we extend the literature by exploring the use of signature regression on the signature of age-specific mortality rates as an alternative to enhancing the robustness of mortality forecasts.

2.2. The signature of a path

2.2.1. The concept of signatures

In statistical analysis, particularly when dealing with time series or sequential data, understanding the concept of a “path” is crucial. A path $X\,:\, [0,1] \rightarrow \mathbb{R}^d$ is defined as a continuous mapping from an interval to a topological space, which is assumed to be piecewise differentiable and of bounded variation. Unlike time series or stochastic processes, which are explicitly governed by time, a path emphasizes its geometrical movements through space, making the term “path” more suitable. From a computational perspective, think of a path as input data with d dimensions, analogous to multidimensional time series data.

Signatures are mathematical tools used to encode the information contained in a path. Using Chen’s identity (Chen, Reference Chen1957), the signature of a piecewise linear path is computed by taking the tensor product of exponentials of the linear segments of the path. Thus, the signature of a path, or simply signatures, is a collection of all iterated integrals between coordinates of a path. This collection is often described as an abstract summary of the input signal’s variability.

Definition 2.1 (Signature of a path (Fermanian, Reference Fermanian2022)). Let $X\,:\,[0,1]\rightarrow \mathbb{R}^d$ be a path of bounded variation, and let $I=(i_1,\ldots,i_k)\in\{1,\ldots,d\}^k$ , where $k\in\mathbb{N}^*$ , be a multi-index of length k. The signature coefficient of X corresponding to the index I is defined by:

\begin{align*} S^I(X) & = \idotsint\limits_{0 \leq u_1 \lt \ldots \lt u_k \leq 1} dX_{u_1}^{i_1} \ldots dX_{u_k}^{i_k} \\ & = \int_{0}^{1}\int_{u_1}^{1} \int_{u_2}^{1} \ldots \int_{u_{k-1}}^{1} dX_{u_k}^{i_k} \ldots dX_{u_2}^{i_2} dX_{u_1}^{i_1}\text{,} \notag \end{align*}

where $S^I(X)$ is called the signature coefficient of order k. The signature of X is then the sequence containing all the signature coefficients:

\begin{equation*} S(X)=(1, S^{(1)}(X), \ldots, S^{(d)}(X), S^{(1,1)}(X), \ldots, S^{(i_1,\ldots,i_k)}(X), \ldots). \end{equation*}

However, from a practical standpoint, it is not feasible to use all the signature coefficients in the infinite sequence of the signature of X due to computational limits and also due to the factorial decay of the coefficients as the order increases (Bonnier et al., Reference Bonnier, Kidger, Arribas, Salvi and Lyons2019). Therefore, it is often sufficient to consider the signature of a path up to a truncated order of m, denoted as:

\begin{equation*} S^m(X)=(1, S^{(1)}(X), \ldots, S^{(d)}(X), S^{(1,1)}(X), \ldots, S^{(i_1,\ldots,i_m)}(X)), \end{equation*}

which contains $s_d(m)=\frac{d^{m+1}-1}{d-1}$ number of signature coefficients. Note that the number of signature coefficients grows exponentially with the truncation order.

A fundamental attribute of the signature is its ability to encode the geometric properties of a path. The first-order terms, $S^{(i)}(X)$ , are straightforward, representing the total change along each dimension of the path or the increments observed in each variable. The second-order terms, $S^{(i_1, i_2)}(X)$ , provide insights into the interaction between pairs of dimensions, similar to capturing the area enclosed by the path as it travels through these dimensions, and higher-order interactions are particularly valuable for understanding complex dependencies and dynamics within the data.

2.3. Embedding with signatures

Embedding paths before computing its signature is a standard practice used to maximize the potential of signatures by addressing inherent invariances such as translation and time reparametrization that can result in significant information loss. These embedding transformation ensure that all relevant information are preserved and encoded into the signature and in turn enrich the signature’s capability to reflect the complex dynamics within the path.

There are various embedding techniques that can be used to augment paths, each suited to specific applications based on the aspects of the data they capture. Different embedding choices highlight different features of the data, making them appropriate for particular tasks (Fermanian, Reference Fermanian2021). In this paper, we focus mainly on three types of embeddings: basepoint augmentation, time augmentation, and the lead-lag transformation.

The basepoint augmentation (Kidger and Lyons, Reference Kidger, Lyons, Abernethy and Agarwal2020) involves inserting a fixed initial point to the path data. By anchoring the sequence, basepoint augmentation allows the signature to overcome the translation invariance of the signature.

Time augmentation counteracts the signature’s invariance to time reparametrization, by introducing a monotone coordinate that serves as a virtual timestamp to each point in the path. This embeds precise temporal information into the path, making the signature sensitive to the sequence’s timing. Besides, this also guarantees the uniqueness of the signatures computed (Hambly and Lyons, Reference Hambly and Lyons2010).

To further aid the sequential learning capabilities of the signature, the lead-lag transformation extends the path’s dimensionality by interleaving the original path with its lagged version, which enriches the data representation, allowing the signature to capture not only the path’s immediate direction but also its past states. Fermanian (Reference Fermanian2021) noted the lead-lag transformation improves a model’s predictive capabilities and is beneficial for tasks involving sequential data.

In the proceeding sections, we will use these embedding transformation to prepare our data paths before computing its signature so as to enhance the signatures’ ability to capture and reflect the complex dynamics of the sequences.

2.4 Signature regression

One of the many practical applications of signatures are its usage as feature maps in a regression framework. This rich representation allows for the modeling of complex relationships within functional data. The ability of signatures to compactly encode information about paths makes them particularly powerful for handling high-dimensional and complex datasets, enabling the application of regression techniques in scenarios where conventional methods might be inadequate due to the complexity and dimensionality of the data. These are motivated and guaranteed by many results as pointed out by the shuffle identity, the signature approximation theorem (Levin et al. Reference Levin, Lyons and Ni2016), and the universal nonlinearity of paths (Bonnier et al., Reference Bonnier, Kidger, Arribas, Salvi and Lyons2019).

Definition 2.2 (Universal nonlinearity (Bonnier et al., Reference Bonnier, Kidger, Arribas, Salvi and Lyons2019)). Let f be a real-valued continuous function on continuous piecewise smooth paths in $\mathbb{R}^d$ and let $\kappa$ be a compact set of such paths. Assume that $X_0 = 0, \forall X \in \kappa$ , let $\epsilon \gt 0$ , there exists a linear functional L such that,

\begin{equation*} |f(X) - L(S(X))| \lt \epsilon . \end{equation*}

Kiraly and Oberhauser (Reference Kiraly and Oberhauser2019) then gathered that a linear combination of signature features are capable of approximating continuous functions arbitrarily well, meaning signatures are able to linearize functions of a path, further emphasizing the universality (Levin et al., Reference Levin, Lyons and Ni2016) of signatures. To further clarify, consider a classical regression setting with a path X that maps some space into $\mathbb{R}^d$ and Y a response, we have

\begin{equation*}Y=f(X)+\epsilon, \end{equation*}

where f is some functional that characterizes the relationship between the response and the path. Assuming f is nonlinear, it is then generally challenging to determine the function f accurately due to the increased complexity in parameter estimation, model interpretation, and computational resources required for nonlinear regression or possibly neural networks, compared to linear regression. However, if the signature of path X is considered as the predictors instead, we can redefine the relationship with a linear functional g such that

\begin{equation*}Y=g(S(X))+\epsilon.\end{equation*}

Although in theory linear regression is sufficient for estimating g, in practice, regularization techniques are used to handle the multicollinearity of signature features (Levin et al., Reference Levin, Lyons and Ni2016) as linear regression models can become unstable when variables are not independent. This instability often results in large variances in the estimated coefficients, making the model unreliable.

2.5. The Hyndman–Ullah model with truncated signatures

We are now prepared to integrate signatures into the HU model, which we will refer to as the HU with truncated signatures (HUts) model. The integration of signatures with functional data models have been documented by Fermanian (Reference Fermanian2022) and Frévent (Reference Frévent2023). However, these studies primarily applied signature regression in a functional regression framework with scalar responses and functional covariates, where signatures of the functional covariates were used directly. In contrast, the HU model operates more like a functional regression model with functional responses, making the direct application of signature regression more complex. We plan to consider using the signature of cross-sectional data from these functional covariates to apply signature regression at each point throughout the functional domain and then reconstruct a functional response. The HUts model is presented as follows:

\begin{equation*} y_t(x)=f_t(x)+\sigma_t(x)\epsilon_t, \end{equation*}

\begin{equation*} f_t(x)=\mu(x)+\sum_{k=1}^K \beta_{t,k} Z_k(x) +e_t(x), \end{equation*}

where $Z_k(x)$ is the $k\mathrm{th}$ principal component of the signature of the transformed path of age-specific rates, and its respective scores $\beta_{t,k}$ .

The smoothing approach of the mortality rates are remained as the original HU model, which is with the constrained weighted penalized regression splines to obtain $\hat{f}_t(x)$ . This ensures the model retains the capacity to account for handling heteroskedasticity in mortality data as variability is different across different ages or time periods. Although more advanced smoothing techniques such as bivariate smoothing (Dokumentov et al., Reference Dokumentov, Hyndman and Tickle2018) have been put forth, they also introduce additional complexity and computational demands.

Additionally, the estimation of the mean function is also kept the same as in the HU model, by averaging all the smoothed curves:

\begin{equation*} \hat{\mu}(x)=\frac{1}{n}\sum_{t=1}^{n}\hat{f}_t(x). \end{equation*}

This method provides a straightforward and effective way to capture the central tendency of the functional responses. Now, we will decompose the centered curves, $f^*_t(x)=\hat{f}_t(x)-\hat{\mu}(x)$ with signature regression over a fine grid of q equally spaced values $x^*_{1}, \ldots, x^*_{q} \in[x_1,x_p]$ .

Consider the time series of the centered smoothed mortality rates for age $ x_i $ , $ \{f^*_t(x_i)\}_{t=1}^{n} $ . Combining the basepoint augmentation, time augmentation, and the lead-lag transformation, the augmented path can be transformed as follows:

\begin{align*} \tilde{X_i} = \Bigl(&(t_0, 0, 0), (t_1, f^*_1(x_i), f^*_1(x_i)), (t_2, f^*_2(x_i), f^*_1(x_i)), (t_3, f^*_2(x_i), f^*_2(x_i)), \\ & (t_4, f^*_3(x_i), f^*_2(x_i)), \ldots, (t_n, f^*_n(x_i), f^*_n(x_i)) \Bigr), \end{align*}

where t is an ordered sequence within the interval [0,1] and $t_0 = 0 \lt t_1 \lt t_2 \lt \ldots \lt t_n = 1$ , marking the times of observations. This structured augmentation ensures that each data point $x_i$ is not only represented by its current value but also linked to its immediate past, enhancing the path’s dimensionality and capturing both the direction and history within the data, thereby enriching the data representation for signature computation. The transformed time series $\tilde{X}_i$ is then used to compute the truncated signature $S^m(\tilde{X}_i)$ .

To estimate $\beta_{t,k}$ , we then set up our signature regression framework:

\begin{equation*}f^*_t(x_{i})=g_t(S^m(\tilde{X}_i))+\hat{e}_t(x_i)\end{equation*}

for some linear function $g_t$ , which we will use principal component regression to solve. Applying principal component analysis (PCA) helps by reducing the dimensionality of the signature coefficient feature space to a smaller set of k uncorrelated principal components denoted as $Z_k(x_i)$ . This reduction not only simplifies the model but also enhances its stability and efficiency by focusing on the most informative aspects of the dependency within the signatures. Additionally, using PCA aids in orthogonalizing the terms, which simplifies the calculation of forecast variance later and resembles the approaches used in both the original HU and LC models.

For computational simplicity, we define the matrices:

(2.1) \begin{align}S = \begin{bmatrix} S^m(\tilde{X}_1) \\ \vdots \\ S^m(\tilde{X}_q)\end{bmatrix}_{q \times s_d(m)} \;, \hspace{1cm}Y_t = \begin{bmatrix} f^*_t(x_1^{{*}}) \\ \vdots \\ f^*_t(x_q^{{*}})\end{bmatrix}_{q \times 1}. \end{align}

We then normalize S to obtain $S^*$ before applying PCA to obtain the orthogonal principal component scores by singular value decomposition, $S^* = U \Sigma V^T$ where the first K principal component scores are the first K columns of $Z = U \Sigma$ . Then, to obtain $\hat{\beta}_t$ is equivalent to solving

\begin{equation*} \hat{\beta}_t = (Z_K^\top Z_K )^{-1}Z_K^\top Y_t, \end{equation*}

where $Z_K$ is a $q \times K$ matrix where the $i\mathrm{th}$ row contains the principal component of $S^*$ of order m of age $x_i$ .

By combining everything together, we have

\begin{equation*} \hat{f}_t(x)=\hat{\mu}(x)+\sum_{k=1}^K \hat{\beta}_{t,k} \hat{Z}_k(x) +\hat{e}_t(x), \end{equation*}

where $\hat{Z}_k(x)$ is obtained via interpolation.

2.6. Forecasting with the HUts model

2.6.1. Point forecast

Similar to that of the HU model, the time-dependent coefficients are extrapolated by fitting them into univariate time series $\{\hat{\beta}_{t,k}\}_{t=1}^n$ . We preserve the choice of using univariate ARIMA models as in the HU model to forecast each coefficient h-steps-ahead, $\hat{\beta}_{n+h|n,k}$ , for each $k=1,\ldots,K$ . Subsequently, the forecasted log mortality rates can be calculated with:

\begin{equation*} \hat{y}_{n+h|n}(x)=\hat{\mu}(x)+\sum_{k=1}^K \hat{\beta}_{n+h|n,k} \hat{Z}_k(x). \end{equation*}

Since each term is constructed to be orthogonal to one another, the forecast variance can also be approximated, assuming all sources of error are normally distributed, by summing the component variances similar to that of the forecast variance of the HU model (Hyndman and Ullah, Reference Hyndman and Ullah2007):

(2.2) \begin{equation} \text{Var}(y_{n+h|n}(x)) \approx \hat{\sigma}^2_{\mu}(x) + \sum_{k=1}^K \hat{Z}_k^2(x) u_{n+h|n,k} + v(x) + \sigma^2_{n+h|n}(x),\end{equation}

where $\hat{\sigma}^2_{\mu}(x)$ is the variance obtained from estimating the mean function, $ u_{n+h|n,k}$ is the variance from the time series used to forecast $\hat{\beta}_{n+h|n,k}$ , v(x) are model fit errors estimated by averaging $\hat{e}_t(x)$ , and finally $\sigma^2_{n+h|n}(x)$ , the observational variance.

2.6.3. Bootstrap prediction interval

Since the HUts model, like the wHU model, do not satisfy the normality assumption (see Appendix A in the Supplementary Material), we employ the same method used to construct bootstrap prediction intervals for the wHU model in Hyndman and Shang (Reference Hyndman and Shang2009) to construct prediction intervals for the HUts model.

The procedure begins by utilizing univariate time series models to generate multi-step-ahead forecasts for the scores of principal components, $\{\beta_{t,k}\}_{t=1}^n$ . Subsequently, for $t = h+1, \ldots, n$ , we calculate the forecast errors for h-step-ahead forecasts as:

\begin{equation*} \hat{\xi}_{t,h,k} = \hat{\beta}_{t,k} - \hat{\beta}_{t|t-h,k}. \nonumber\end{equation*}

To create a bootstrap sample of $\beta_{n+h,k}$ , we randomly resample these forecast errors with replacement, designated as $\hat{\xi}_{*,h,k}^{(l)}$ , and add them to the h-step-ahead forecast, $\hat{\beta}_{n+h|n,k}$ , yielding:

\begin{equation*} \hat{\beta}_{n+h|n,k}^{(l)} = \hat{\beta}_{n+h|n,k} + \hat{\xi}_{*,h,k}^{(l)}, \nonumber\end{equation*}

for $l=1, \ldots, L$ .

The model’s residual, $\hat{e}_{n+h|n}^{(l)}(x)$ is bootstrapped from the residuals $\{\hat{e}_1(x), \ldots, \hat{e}_n(x)\}$ by sampling with replacement. Similarly, smoothing errors $\hat{\epsilon}_{n+h,i}^{(l)}$ can be bootstrapped by resampling with replacement from $\{\hat{\epsilon}_{1,i}, \ldots, \hat{\epsilon}_{n,i}\}$ .

When constructing the L variants for $\hat{y}_{n+h|n}^{(l)}(x)$ , all possible components of variability are combined, assuming they are uncorrelated to each other by using:

\begin{equation*}\hat{y}_{n+h|n}^{(l)}(x_i) = \hat{\mu}(x_i) + \sum_{k=1}^K \hat{\beta}_{n+h|n,k}^{(l)} \hat{Z}_k(x_i) + \hat{e}_{n+h|n}^{(l)}(x_i) + \hat{\sigma}_{n+h}(x_i)\hat{\epsilon}_{n+h,i}^{(l)}. \nonumber\end{equation*}

Subsequently, the $(1-\alpha)100\%$ prediction intervals, $[l_h(x),u_h(x)]$ can then be obtained from the $(1-\alpha/2)$ and $(\alpha/2)$ quantiles of these bootstrap variants. We also consider the adjustment made to the bootstrap prediction intervals in Hyndman and Shang (Reference Hyndman and Shang2009):

\begin{equation*} [0.5\{l_h(x)+u_h(x)\}-\{u_h(x)-l_h(x)\}p(x),0.5\{l_h(x)+u_h(x)\}+\{u_h(x)-l_h(x)\}p(x)] \nonumber\end{equation*}

to account for the uncertainty in our forecast by calibrating the intervals using one-step-ahead forecast errors. Here, $p(x)=d(x)/[u_1(x)-l_1(x)]$ , where d(x) is the empirical differences between the $(1-\alpha/2)$ and $(\alpha/2)$ quantiles of the in-sample one-step-ahead forecast errors $\{\hat{f}_{t+1}(x)-\hat{f}_{t+1|t}(x); t=r,\cdots,n-1\}$ and r is the smallest number of observations used to fit the model.

3.1. Data description

This study utilized datasets from 12 different countries, all sourced from the Human Mortality Database (2024). Of the 12 countries considered, data spanning all available years up to 2015 were extracted, with age groups ranging from zero to a carefully determined maximum age to avoid zero or missing values in elderly data. For ages beyond this established maximum, they were combined into an upper age group. In addressing intermediate missing values within the data for Denmark, Finland, and Norway, a uniform approach was applied. A constant, represented by the smallest death rate value in the dataset, was added to each death rate prior to taking its logarithm. This adjustment facilitated a consistent treatment of missing values across the dataset. The selected countries are tabulated in Table 1, along with the commencing years and age groups.

Table 1. Mortality data by country.

3.2. Point forecast

The log mortality rates of 12 countries were forecasted with the HUts model using a truncation order $m=2$ (see Appendix B in the Supplementary Material). It is therefore intuitive to compare the forecast accuracy of the proposed method with the HU model, HUrob model, and the wHU model, which are commonly used on mortality rates of a single population. The selection of the number of principal components K can impact the forecasting performance of the HUts model and the HU model variants. They can be chosen based off a threshold function that measures the proportion of variance explained by each principal component, or by minimizing the integrated squared forecast error (Hyndman and Ullah, Reference Hyndman and Ullah2007). However, Hyndman and Booth (Reference Hyndman and Booth2008) noted that the HU models are relatively insensitive to the choice of K as long as it is sufficiently large, and that setting $K=6$ is more than adequate to capture the essential features of the data. For consistency across all models, we adopt $K=6$ for both the HU model variants and the proposed HUts model. The models were fitted using an expanding window approach, incorporating data from the respective commencing years of each country as outlined in Table 1. The initial forecasting period covers the most recent two decades which ends at 2015. We fit this data using both models and produce forecasts up to 10 years ahead. The accuracy of these forecasts are then evaluated. Subsequently, the fitting period is then extended by a year and the forecasts for the same horizons are recalculated. This procedure of expanding the training window continues until it includes data up to the year 2014.

In order to gauge the accuracy of the point forecast of the HUts model, the mean squared errors (MSE) are calculated and compared with those of the forecasts of the HU model variants. MSE quantifies the average squared difference between the predicted and observed values, providing a measure of the overall accuracy and precision of predictions with lower values indicating better forecasting performance. The MSE is computed using:

(3.1) \begin{equation} MSE(h) = \frac{1}{pq}\sum^q_{t=1}\sum^p_{i=1}(y_t(x_i)-\hat{y}_{t|t-h}(x_i))^{2},\end{equation}

where h is the forecasting horizon, p is the number of ages, q is the number of years considered in the forecasting period, $y_t(x_i)$ is the observed log mortality, and $\hat{y}_{t|t-h}(x_i)$ is the predicted value.

Figure 1 presents the plot of the root MSE of $h=1,5,10$ forecasts across the different models. The square roots of the MSE values are plotted instead of the MSE themselves to ensure the plot scale accommodates all data. Readers are directed to Appendix D in the Supplementary Material for the tabulated MSE values of all countries. The results indicate that while the wHU model often achieves the lowest MSE values for one-step-ahead forecasts, the HUts model still performs competitively, frequently achieving values close to the best-performing models. In comparison, the HU model shows higher error rates, followed by the HUrob model.

Figure 1. Root MSE of (a) one, (b) five, and (c) ten-step-ahead point forecasts of log mortality rates by model and country.

For five-step-ahead forecasts, the HUts model demonstrates more competitive performance. Figure 1(b) shows that overall the HUts model achieves lower error metrics compared to the other models, suggesting an improvement in the HUts model’s forecasting ability as the forecast horizon increases.

The results for ten-step-ahead forecasts, shown in Figure 1(c), further confirms the strong performance of the HUts model. It consistently achieves lower or comparable MSE, reinforcing its effectiveness in long-term mortality rate forecasting. In most countries, the HUts model has lower MSE values, with exceptions to a few cases such as in the United States, where the HU model occasionally performs better. This consistent performance across a majority of countries highlights the robustness of the HUts model in capturing long-term trends in mortality rates.

Across different forecasting horizons, the best model is generally a choice between the HUts and wHU models, as these models tend to have the lowest MSE values. The wHU model performs better at shorter forecast horizons, as seen in Figure 1(a). However, as the forecasting horizon increases, the HUts model demonstrates improved performance and eventually surpasses the wHU model. This shift indicates that the HUts model’s accuracy and reliability become more pronounced with longer-term forecasts.

Notably, in the cases of France and Ireland, the HUts model performs consistently well across all forecasting horizons, maintaining lower MSE values compared to other models. Conversely, for Australia, the HUts model does not appear to be the most favorable at any point across the forecast horizons, with other models, such as the wHU model, often yielding better results. The following section will take a closer look at the mortality rate forecasts for France and Australia, exploring the nuances and potential reasons behind the differing performances of the HUts model against the HU model variants.

3.2.1. A closer look at the French and Australian mortality data

To assess the robustness of the HUts model, we examine the mortality rates of France and Australia more closely. Figures 2(a) and 3(a) include the observed log mortality rates, and Figures 2(b) and 3(b) show the smoothed log mortality rates with red representing the earlier years. Figures 2(c) and 3(c) the smoothed log mortality rates viewed as a time series with red representing the earlier ages. France’s mortality data, notably influenced by periods of war and disease outbreaks, presents distinct outliers which are evident in the volatility observed in their mortality trends, as shown in Figure 2 which are most visible around the 1920s and 1940s. This provides an oppurtunity to test for the HUts model’s ability to handle anomalies in data. In contrast, Australia’s mortality data, starting from 1921 appears to have a smoother trend, without any significant outliers as depicted in Figure 3.

Figure 2. France (1899–2015); (a) observed log mortality rates; (b) smoothed log mortality rates; (c) times series for ages.

Figure 3. Australia (1921–2015); (a) observed log mortality rates; (b) smoothed log mortality rates; (c) times series for ages.

The ability of the HUts model to effectively handle anomalies in French mortality is illustrated in Figure 4(a), where the model demonstrates a substantial reduction in MSE across various forecast horizons when compared to the HU model. This significant improvement shows the HUts model’s robustness against data irregularities, making it a reliable choice for forecasting in scenarios with potential outliers. Additionally, the HUts model showcases lower MSE than both the HUrob and wHU models, further affirming its superior performance in handling French mortality. Conversely, Australia presents a more stable dataset with mortality rates following a consistent trend and lacking significant outliers. The MSE plot for Australia in Figure 4(b) shows less pronounced differences compared to France. While the HUts model demonstrates superior performance over the HU and HUrob models, affirming its effectiveness even in more uniform datasets, the wHU model performs better than the HUts model. This suggests that the wHU model might be more suited for datasets with consistent trends and fewer irregularities. This comparative analysis highlights the importance of selecting appropriate models based on the characteristics of the dataset, with the HUts model excelling in more variable datasets and the wHU model in more stable ones. The MSE plots of the remaining countries can be found in Appendix E in the Supplementary Material.

Figure 4. MSE plots for (a) French and (b) Australian mortality of the HUts, wHU, HUrob, and HU models.

To gain deeper insights into the model performance, the MSE values are considered across the age spectrum. Instead of using Equation (3.1), the MSE values are averaged only over the years of the forecasting period. This approach will allow us to observe how the models perform at different age groups, providing a more granular understanding of the forecast accuracy. The MSE across ages are computed with the following:

\begin{equation*} MSE(h,x_i) = \frac{1}{q}\sum_{t=1}^{q}(y_t(x_i)-\hat{y}_{t|t-h}(x_i))^2. \nonumber\end{equation*}

Upon a closer inspection of the MSE across different age groups of the French mortality data in Figures 5(a), 6(a), and 7(a), the HUts model generally surpasses every other model in terms of forecasting performance. For one-step-ahead forecasts (Figure 5(a)), the HUts and wHU models show similar performance in terms of MSE, with both models outperforming the HU and HUrob models substantially. For five-step-ahead forecasts (Figure 6(a)), the HUts model continues to perform well, maintaining lower MSE values across most age groups compared to the other models. The HUts model shows improved performance at middle ages (30–50 years) over the wHU model, while the HU and HUrob models, although better at these ages, struggle during early and older ages, contributing to their higher overall error metric. For ten-step-ahead forecasts (Figure 7(a)), the HUts model’s performance remains robust, achieving lower MSE values across most age groups. However, it still faces competition from the wHU model, which performs better in certain age ranges. Overall, these observations suggest that while the HUts model performs exceptionally well, particularly as the forecast horizon increases, the wHU model also shows strong performance, especially at shorter forecast horizons. The choice between these models often depends on the specific age group and forecasting horizon, with the HUts model becoming more favorable as the horizon increases.

Figure 5. One-step-ahead forecast MSE averaged over years in forecasting period of (a) French and (b) Australian mortality for the HUts, wHU, HUrob, and HU models.

Figure 6. Five-step-ahead forecast MSE averaged over years in forecasting period of (a) French and (b) Australian mortality for the HUts, wHU, HUrob, and HU models.

Figure 7. Ten-step-ahead forecast MSE averaged over years in forecasting period of (a) French and (b) Australian mortality for the HUts, wHU, HUrob, and HU models.

Similarly, MSE for the Australian mortality forecasts in Figures 5(b), 6(b), and 7(b) highlight the competitive nature of the HUts model, particularly at shorter forecast horizons. For one-step-ahead forecasts (Figure 5(b)), the HUts model performs comparably to the wHU model, with both showing nearly identical MSE values. This indicates that the HUts model is effective at short-term forecasting. At the five-step-ahead horizon (Figure 6(b)), the HUts model continues to perform well but begins to show higher MSE values in the 20–40 age range compared to the wHU model. This underperformance in mid-age groups becomes more evident, reflecting some potential biases that start to emerge at this forecast horizon. By the ten-step-ahead forecasts (Figure 7(b)), the HUts model, although still competitive, shows further increased MSE values around ages 20–40 years. This trend highlights the model’s underestimation (see Appendix F of the Supplementary Material) in this specific age group, significantly affecting the overall error metric. Despite this, the HUts model maintains reasonable performance across other age groups, yet the wHU model generally achieves better results at these longer horizons, suggesting that the wHU model might be better suited for datasets with more consistent trends and fewer irregularities.

The findings of French and Australian mortality data present a comparative analysis of the HUts model’s performance under differing demographic conditions and historical events influencing mortality data. While the HUts model is able to perform well in stable demographic conditions, its advantages are particularly noticeable when confronting data with potential outliers, making it a robust tool for forecasting when compared with the HU model variants.

3.3. Interval forecast

To evaluate the interval forecasts, the empirical coverage probability is calculated using the following expression (Hyndman and Shang, Reference Hyndman and Shang2009):

\begin{equation*} \frac{1}{qph} \sum_{t=1}^q \sum_{j=1}^h \sum_{i=1}^p \mathbf{1}\left(\hat{y}_{t+j|t}^{(\alpha/2)}(x_i) \lt y_{t+j}(x_i) \lt \hat{y}_{t+j|t}^{(1-\alpha/2)}(x_i)\right), \nonumber\end{equation*}

where $\mathbf{1}(\cdot)$ is the indicator function and $\hat{y}_{t+j|t}^{(\alpha/2)}(x_i)$ is the $\alpha/2$ -quantile from the bootstrapped samples. It quantifies whether the intervals are adequately wide to encompass the actual observations at a specified nominal coverage of $(1-\alpha)100\%$ . 95% prediction intervals of the HUts model and the wHU model are constructed for all the 12 countries and tabulated in Table 2. The minimum observations required to obtain the one-step-ahead forecast errors used for the adjustment, r is set to 30.

Table 2. Empirical coverage probabilities of bootstrapped prediction intervals.

The bootstrapping procedure generally performs well in obtaining prediction intervals that are close to the nominal 95% coverage level. Across various countries and forecasting methods, the unadjusted 95% prediction intervals yield empirical coverage probabilities that align reasonably well with the expected nominal level. These results demonstrate that the bootstrap approach reliably captures the underlying variability in mortality data, ensuring that the prediction intervals accurately achieve the desired 95% coverage.

However, when applying the Hyndman and Shang (Reference Hyndman and Shang2009) adjustment to these prediction intervals, the effectiveness of the adjustment varies across different countries and methods. This adjustment aims to refine the intervals based on one-step-ahead forecast errors, potentially bringing the empirical coverage probabilities closer to the nominal level. Yet, its impact is not uniformly beneficial and depends on the specific characteristics of the mortality data and the forecasting models used.

In the case of Australia, both the prediction intervals of the HUts and wHU models show improvements in coverage probabilities after the adjustment. For the HUts model, the empirical coverage increased from 92.70% to 94.77%, moving closer to the nominal 95% level. Similarly, for the wHU model, the coverage improves from 88.30% to 89.22%. This indicates that the adjustment effectively enhances the reliability of the prediction intervals in these instances. A similar positive effect is observed for Denmark with the wHU model, where the coverage probability significantly increases from 78.06% to 91.17% after adjustment, as well as from 92.99% to 95.84% with the HUts model.

On the other hand, the adjustment does not always lead to better coverage probabilities. For example, in Belgium using the wHU model, the coverage probability decreases from 91.05% to 81.13% after adjustment, moving further away from the nominal level. Similarly, in Bulgaria with the HUts model, the coverage dropped from 91.89% to 85.44% post-adjustment. These instances indicate that the adjustment may not universally improve the prediction intervals possibly due to the nature of the forecast errors of the fitted models.

A notable case is France, where the empirical coverage probabilities reaches 100% for both unadjusted and adjusted intervals with the HUts model. This phenomenon is caused by the prediction intervals absorbing a significant amount of uncertainty from the forecasts, resulting in wider intervals as seen in Figure 8. For the unadjusted prediction intervals, the uncertainty arises from bootstrapping the model errors. Due to the presence of outliers in the French mortality data, these model errors are larger. When the prediction intervals are adjusted, they incorporate additional variability from the in-sample one-step-ahead forecast errors. These forecast errors are also inflated by outliers, resulting in even wider intervals. This further widening makes the predictions more conservative to ensure coverage despite the high variability in the data. The same effect is observed with the wHU model, where the prediction intervals are similarly wide which suggests that outliers can impact the bootstrapping approach, potentially reducing the reliability of the intervals.

Figure 8. Bootstrapped (a) 95% (b) adjusted 95% prediction intervals of one-step-ahead forecast for French mortality (1899–2014) for the HUts model and the wHU model.

In contrast, when examining the prediction intervals for Australia in Figure 9, they are significantly narrower than those for France. Australia’s mortality data is more stable, with fewer outliers, leading to less uncertainty in the forecasts. Consequently, the prediction intervals are narrower. This comparison highlights the point that uncertainty within mortality data leads to wider bootstrapped prediction intervals, as seen in France, while more stable data results in narrower intervals, as observed in Australia.

Figure 9. Bootstrapped (a) 95% (b) adjusted 95% prediction intervals of one-step-ahead forecast for Australian mortality (1921–2014) for the HUts model and the wHU model.