2.1. Background and description

In response to the COVID-19 pandemic, the Human Mortality Database launched the STMF series (STMF, 2021) in May 2020, attracting growing interest in the analysis of excess mortality. Jdanov et al. (Reference Jdanov, Galarza, Shkolnikov, Jasilionis, Németh, Leon, Boe and Barbieri2021) develop an online tool to visualize weekly excess deaths across countries, using 2010–2019 as the reference period. Based on STMF data from 2016 to 2020, Islam et al. (Reference Islam, Shkolnikov, Acosta, Klimkin, Kawachi, Irizarry, Alicandro, Khunti, Yates and Jdanov2021) find that excess deaths in 29 high-income countries substantially exceed officially reported COVID-19 deaths, particularly among children under 15 years of age. This result is supported by Karlinsky and Kobak (Reference Karlinsky and Kobak2021), who underscore the need for timely and reliable mortality data to improve forecast accuracy. Kontis et al. (Reference Kontis, Bennett, Rashid, Parks, Pearson-Stuttard, Guillot, Asaria, Zhou, Battaglini and Corsetti2020) focus on countries with populations over four million and incorporate seasonal patterns, temperature effects, long-term trends, and weekly dependencies into their model. Their findings reveal that individuals aged 65 years and older account for 94% of total excess deaths, both in absolute and relative terms. Vanella et al. (Reference Vanella, Basellini and Lange2021) apply the principal component analysis to mortality data from 19 European countries. They demonstrate that country-specific models often overlook shared common trends, reducing forecast accuracy. Their study underscores the importance of accounting for mortality correlations across countries, particularly given the transnational spread of pathogens. They also identify substantial cross-country variation in excess mortality, driven by differences in public health responses and demographic structures. Similarly, Arik et al. (Reference Arik, Klein and Li2026) present a cross-country analysis of excess mortality during 2020–2023, showing broadly similar patterns across measures while highlighting that methodological choices can affect the magnitude and variability of the estimates.

While much of the literature has focused on explaining pandemic-related excess mortality, this paper aims to address the broader need for long-term mortality forecasts using high-frequency data. In this paper, we collect mortality data of 30 countries from the STMF series, excluding those with populations under two million due to the frequent occurrence of zero weekly deaths. Such data sparsity introduces substantial uncertainty and can result in unstable model estimates. We extract national-level mortality rates spanning from the second week of 2015 to the final week of 2019, comprising exactly 52 weeks per year. To address the issue of zero mortality values which result in undefined logarithmic transformations, we choose to exclude the youngest age group (0–14 years old) from our analysis. Consequently, we focus on the remaining four age groups: 15–64, 65–74, 75–84, and 85+ years.

2.2. Visualization

The weekly mortality rates across four age groups for all 30 countries over the period 2015–2019 are presented in Appendix A. However, the seasonal patterns are not easily discernible due to the high dimensions of the data. To facilitate a clearer examination of seasonal patterns, we present weekly mortality trends for four representative countries, Australia, England and Wales, Italy, and the USA, in Figure 1. The selection of these countries reflects both their geographic diversity, covering major regions such as Oceania, North America, and Europe, and the opportunity to compare contrasting seasonal patterns between the Northern and Southern Hemispheres. A recurring seasonal pattern is observed across all countries, with each cycle spanning approximately one year. Mortality rates tend to peak during the winter and decline in the summer, reflecting the influence of climate factors. As age increases, the seasonal pattern becomes smoother and more pronounced, suggesting a stronger and more apparent relationship between weather conditions and mortality among elderly individuals. This observation aligns with findings from prior research on weather-related mortality (Huang et al., Reference Huang, Chu, Wang and Barnett2015). Historical evidence indicates that events such as influenza outbreaks, heatwaves, and extreme cold spells have contributed to short-term fluctuations in mortality rates (Jdanov et al., Reference Jdanov, Galarza, Shkolnikov, Jasilionis, Németh, Leon, Boe and Barbieri2021), and the impacts are more pronounced on individuals aged 55 years and older (Li and Tang, Reference Li and Tang2022).

Figure 1.

2015–2019 mortality rates for the four representative countries by age groups.

Figure 1. Long description

The image contains four line graphs showing mortality rates for different age groups in four countries from 2015 to 2020. The graphs are divided into four panels, each representing a different age group: 15-64 years old, 65-74 years old, 75-84 years old, and 85+ years old. Each graph plots the mortality rate on the y-axis against the date on the x-axis. The countries represented are Australia, England and Wales, Italy, and the USA, each depicted by a different colored line: blue for Australia, gray for England and Wales, pink for Italy, and green for the USA. The mortality rates vary across the age groups and countries, with noticeable peaks and troughs over the years. The graphs illustrate how mortality rates fluctuate over time and differ by age group and country. All values are approximated.

We then take a closer look at the seasonal timing differences between the Southern and Northern Hemispheres. In Australia, peak mortality rates occur around the middle of the year, corresponding to the winter in the Southern Hemisphere. In contrast, Italy, England and Wales, as well as the USA, exhibit peak mortality at the start of the year, reflecting their winter seasons in the Northern Hemisphere. Notably, beyond the typical U-shaped annual mortality pattern, a W-shaped seasonal trajectory appears in certain countries during specific years. This pattern reflects two distinct mortality peaks within a year: a primary and more pronounced peak in winter, commonly associated with cold waves, and a secondary, smaller summer peak likely driven by heatwaves. The W-shaped seasonal pattern is particularly evident among older age groups (75–84 and 85+ years) and is especially pronounced in Italy and England and Wales.

3.1. LL model

Li and Lee (Reference Li and Lee2005) propose the multi-population version of the Lee–Carter model to simultaneously estimate and forecast the mortality rates of a group of populations with similar economic conditions and close-knit connections. The Lee–Carter model proposed by Lee and Carter (Reference Lee and Carter1992) serves as the foundation for the estimation of the LL model, and its definition and estimation algorithm are summarized in Appendix B. Li and Lee (Reference Li and Lee2005) develop an augmented common factor model: for age group $x\in\{1,...,N\}$ , time $t\in\{1,...,T\}$ , and country $j\in\{1,...,J\}$ ,

\begin{equation*}\log\!(m^{j}_{x,t})=A^{\,j}_x + b^{\,p}_x \kappa^{\,p}_t + b^{\,j}_x \kappa^{\,j}_t + \varepsilon^{\,j}_{x,t},\end{equation*}

where $m_{x,t}^{\,j}$ is the mortality rate for age group x and country j at time t, $A^{\,j}_x$ is the age-specific mean of log mortality rate over time for country j, and $\varepsilon^{\,j}_{x,t}$ is the corresponding error term for age group x, time t, and country j. The common trend is captured by $b^{\,p}_x$ and $\kappa^{\,p}_t$ , obtained by applying the Lee–Carter method for the combined group including all the populations to avoid long-term divergence in mean mortality forecasts. In addition, $b^{\,j}_x$ and $\kappa^{\,j}_t$ represent the shorter-term excessive rate of mortality change in country j over the common trend.

The original estimation method proposed by Li and Lee (Reference Li and Lee2005) first extracts the age-specific mean of log mortality rates $A_x^{\,j}$ . Then the common trend $b^{\,p}_x \kappa^{\,p}_t$ across all subpopulations is computed using the singular value decomposition (SVD) on the pooled, mean-adjusted log mortality rates $\log\!(m^{j}_{x,t})-A_x^{\,j}$ . Population-specific residuals are then analyzed via the SVD to obtain additional components $b^{\,j}_x \kappa^{\,j}_t$ to capture the deviations from the common trend. However, to ensure more coherent and simultaneous modeling across populations with other benchmark models, we follow the product-ratio functional method proposed by Hyndman et al. (Reference Hyndman, Booth and Yasmeen2013) with the order chosen to be 1, to estimate the LL model. The detailed estimation is outlined as follows. First, let ${\boldsymbol{y}}^{\,j}$ denote the $T \times N$ (time $\times$ age group) matrix of the observed log mortality rates for country j, whose (t, x)-th element is given by $y^{\,j}_{x,t} = \log\!(\tilde{m}^{\,j}_{x,t})$ . We then re-express the observed mortality rate as $\tilde{m}^{\,j}_{x,t} = \tilde{p}_{x,t}\tilde{r}^{\,j}_{x,t}$ where $\tilde{p}_{x,t}=({\tilde{m}^{1}_{x,t}\times...\times \tilde{m}^{\,j}_{x,t}})^{1/{J}}$ is the product term, calculated as the geometric mean of the observed mortality rate across all countries, and $\tilde{r}^{\,j}_{x,t}=\tilde{m}^{\,j}_{x,t}/\tilde{p}_{x,t}$ represents the ratio term for country j. We use ${\boldsymbol{y}}^{\,j}$ to fit separate Lee–Carter models to the product and ratio terms, to get the common factors $\boldsymbol{\hat{\kappa}}^{\,p}_t, {\hat{\boldsymbol{a}}}^{\,p}_x, {\hat{\boldsymbol{b}}}^{\,p}_x$ and country-specific factors $\boldsymbol{\hat{\kappa}}^{\,j}_t, {\hat{\boldsymbol{a}}}^{\,j}_x, {\hat{\boldsymbol{b}}}^{\,j}_x$ , respectively; see Appendix B for the model specification $a_x$ , $b_x$ , and $\kappa_t$ of the Lee–Carter model. We then compute the overall age-specific mean as ${\hat{\boldsymbol{a}}}^{\,j}_x={\hat{\boldsymbol{a}}}^{\,p}_x + {\hat{\boldsymbol{a}}}^{\,j}_x$ and the final predicted log mortality rate is given by

$$\hat{y}_{x,t}^{\,j}=\log\!(\hat{m}^{\,j}_{x,t})=\hat{A}^{\,j}_x+\hat{b}^{\,p}_x \hat{\kappa}^{\,p}_t+\hat{b}^{\,j}_x \hat{\kappa}^{\,j}_t,$$

for age group x, time t, and country j. The complete estimation procedure is summarized in Algorithm 1.

Algorithm 1:

Estimation of the LL model (Product-Ratio Method with order 1)

3.2. GBLL model

Given concerns about the capacity of existing stochastic mortality models to accurately capture seasonality of mortality rates, we adopt a gradient boosting framework to iteratively model different layers of mortality trends. Gradient boosting is an additive modeling technique that iteratively improves predictive accuracy by correcting residual errors from previous fits (Friedman, Reference Friedman2001). Due to its flexibility and predictive power, the gradient boosting method can capture complex, nonlinear relationships and enhance both model fitting and forecast accuracy. Consequently, this data-driven approach has become increasingly popular in mortality rate modeling (see e.g., Blanes-Selva et al., Reference Blanes-Selva, Ruiz-Garca, Tortajada, Bened, Valdivieso and Garca-Gómez2021; Asgari et al., Reference Asgari, Miri and Asgari2022; Qiu et al., Reference Qiu, Chen, Dincer, Lundberg, Kaeberlein and Lee2022; Neshat et al., Reference Neshat, Phipps, Jha, Khojasteh, Tong and Gandomi2026).

There are two main components in the gradient boosting framework, namely the weak learner and the associated learning rate. Traditionally, a weak learner refers to a simple model with limited predictive power that performs slightly better than random guessing, such as a shallow decision tree. More recently, boosting has been used with other weak learners, including generalized linear models (Bühlmann and Hothorn, Reference Bühlmann and Hothorn2007), generalized additive models (Tutz and Binder, Reference Tutz and Binder2006), and copulas (Hans et al., Reference Hans, Klein, Faschingbauer, Schneider and Mayr2023). By sequentially aggregating each fitted weak learner to the final model with a learning rate that controls the contribution of each weak learner to the final model, it results in a stronger and more comprehensive model with enhanced forecast accuracy.

In this paper, we propose to employ the LL model as the weak learner to leverage domain expertise in mortality modeling, enhancing the interpretability of the overall model. Gradient boosting seeks to construct an ensemble of weak learners defined as $\sum_{g=1}^{l}\gamma_g{\hat{\boldsymbol{y}}}_{x,g}^{\,j}$ , where $\hat{{\boldsymbol{y}}}_{x,g}^{\,j}$ denotes the g-th predicted log mortality rate from the LL model for country j and age group x, $\gamma_g$ represents the associated coefficient, and l is the number of iterations. This procedure involves fitting each weak learner to the residuals and solving a one-dimensional optimization problem to determine its contribution.

To estimate the GBLL model, we first let ${\boldsymbol{y}}^{\,j} = ({\boldsymbol{y}}_1^{\,j},..., {\boldsymbol{y}}_N^{\,j})$ denote the input mortality matrix with dimension $T \times N$ (time $\times$ age group) for country j. In the first iteration, the parameter set $\{\boldsymbol{\hat{\kappa}}_{t,1}^p, {\hat{\boldsymbol{a}}}_{x,1}^p, {\hat{\boldsymbol{b}}}_{x,1}^p, \boldsymbol{\hat{\kappa}}_{t,1}^{\,j}, {\hat{\boldsymbol{a}}}_{x,1}^{\,j}, {\hat{\boldsymbol{b}}}_{x,1}^{\,j}\}$ is obtained by fitting the LL model to ${\boldsymbol{y}}^{\,j}$ . Let ${\hat{\boldsymbol{y}}}_1^{\,j}=({\hat{\boldsymbol{y}}}_{1,1}^{\,j},...,{\hat{\boldsymbol{y}}}_{N,1}^{\,j})$ denote the fitted value of the first iteration. The learning rate $\gamma_1$ is then determined by minimizing the following quadratic loss function:

\begin{equation*}L_1(\gamma) =\sum_{j=1}^{\,J} \left\| {\boldsymbol{y}}^{\,j} - \gamma {\hat{\boldsymbol{y}}}_1^{\,j} \right\|_F^2,\end{equation*}

where $\left\| \cdot \right\|_F$ represents the Frobenius norm that measures the overall magnitude of a matrix, that is, $L_1$ can be equivalently written as:

\begin{equation*}L_1(\gamma)=\sum_{j=1}^{\,J}\sum_{x=1}^{N}({\boldsymbol{y}}_{x}^{\,j}-\gamma {\hat{\boldsymbol{y}}}_{x,1}^{\,j})'({\boldsymbol{y}}_{x}^{\,j}-\gamma {\hat{\boldsymbol{y}}}_{x,1}^{\,j}).\end{equation*}

In the $(g+1)$ -th iteration for $g\ge 1$ , the parameter set $\{\boldsymbol{\hat{\kappa}}_{t,{g+1}}^p, {\hat{\boldsymbol{a}}}_{x,{g+1}}^p, {\hat{\boldsymbol{b}}}_{x,{g+1}}^p, \boldsymbol{\hat{\kappa}}_{t,{g+1}}^{\,j}, {\hat{\boldsymbol{a}}}_{x,{g+1}}^{\,j}, {\hat{\boldsymbol{b}}}_{x,{g+1}}^{\,j}\}$ and the fitted value ${\hat{\boldsymbol{y}}}^{\,j}_{g+1} = ({\hat{\boldsymbol{y}}}_{1,g+1}^{\,j}, \ldots, {\hat{\boldsymbol{y}}}_{N,g+1}^{\,j})$ are obtained by fitting the LL model to residual ${\boldsymbol{E}}^{\,j}_g = ({\boldsymbol{E}}_{1,g}^{\,j}, ..., {\boldsymbol{E}}_{N,g}^{\,j})$ from the g-th iteration, defined by:

\begin{equation*}{\boldsymbol{E}}^{\,j}_{g} =\begin{cases}{\boldsymbol{y}}^{\,j} - \gamma_1 {\hat{\boldsymbol{y}}}^{\,j}_{1} = ({\boldsymbol{E}}_{1,1}^{\,j},..., {\boldsymbol{E}}_{N,1}^{\,j})=({\boldsymbol{y}}_{1}^{\,j} - \gamma_1 {\hat{\boldsymbol{y}}}_{1,1}^{\,j}, ..., {\boldsymbol{y}}_{N}^{\,j} - \gamma_1 {\hat{\boldsymbol{y}}}_{N,1}^{\,j}), & \text{if } g = 1, \\[3pt]{\boldsymbol{E}}^{\,j}_{g-1} - \gamma_g {\hat{\boldsymbol{y}}}^{\,j}_{g} = ({\boldsymbol{E}}_{1,g}^{\,j},..., {\boldsymbol{E}}_{N,g}^{\,j})=({\boldsymbol{E}}_{1,g-1}^{\,j} - \gamma_g {\hat{\boldsymbol{y}}}_{1,g}^{\,j}, ..., {\boldsymbol{E}}_{N,g-1}^{\,j} - \gamma_g {\hat{\boldsymbol{y}}}_{N,g}^{\,j}), & \text{if } g \geq 2.\end{cases}\end{equation*}

The learning rate $\gamma_{g}$ for $g\ge 2$ is determined by minimizing the following quadratic loss function:

\begin{equation*}L_g(\gamma) =\sum_{j=1}^{\,J} \left\| {\boldsymbol{E}}_{g-1}^{\,j} - \gamma {\hat{\boldsymbol{y}}}_g^{\,j} \right\|_F^2.\end{equation*}

The algorithm terminates when (a) all 120 residual series (30 countries, 4 age groups each) have Ljung–Box test p-values of at least 0.05 (Ljung and Box, Reference Ljung and Box1978), or (b) the maximum in-sample fitting MAPE across all age groups, countries and training periods falls below a small threshold. For condition (a), the p-value $p_{x,g}^{\,j}$ for each age group x, country j and iteration g is given by $1 - F_{\chi^2_{52}}\left(T(T+2) \sum_{k=1}^{52} (\rho_{k,x,g}^{(j)})^2/(T-k)\right)$ , where $F_{\chi^2_{52}}(\cdot)$ represents the cumulative distribution function (CDF) of a chi-squared distribution with 52 degrees of freedom, T is the whole training period, lag is chosen at 52, and $\rho_{k,x,g}^{(j)}$ is the sample autocorrelation at lag k for the residual ${\boldsymbol{E}}_{x,g}^{\,j}$ for country j, age group x, and the g-th iteration. For condition (b), the maximum in-sample fitting MAPE across all age groups $x\in\{1,...,N\}$ , countries $j\in\{1,...,J\}$ , and training periods $t\in\{1,...,T\}$ is defined as (see Section 4.2 for more details on the definition of MAPE):

\begin{equation*} \underset{x,j,t}{\max}\left(\mathrm{MAPE}_{x,j,t}\right)=\underset{x,j,t}{\max}\left(\frac{|\hat{m}^{\,j}_{x,t}-\tilde{m}^{\,j}_{x,t}|}{\tilde{m}^{\,j}_{x,t}}\right).\end{equation*}

The algorithm stops if the maximum in-sample fitting MAPE is less than $0.001$ . This constraint prevents unnecessary iterations once a sufficiently accurate fit is reached and avoids reducing the fitting error to values that may reflect overfitting rather than meaningful improvement. The 0.001 threshold aligns with common optimization tolerances, generally ranging from $10^{-6}$ to $10^{-3}$ (Nocedal and Wright, Reference Nocedal and Wright2006). It suggests that the residuals contain no remaining systematic patterns and the model adequately captures the underlying data structure. We set the maximum number of boosting iterations to $G = 50$ . A detailed summary of this procedure is provided in Algorithm 2. The fitted log mortality rates are obtained as follows: for age group $x\in\{1,...,N\}$ , time $t\in\{1,...,T\}$ , and country $j\in\{1,...,J\}$ ,

\begin{equation*}\log\!(\hat{m}^{\,j}_{x,t}) = \sum_{g=1}^{l} \gamma_g\left(\hat{a}^{\,p}_{x,g}+\hat{b}^{\,p}_{x,g}\hat{\kappa}^{\,p}_{t,g}+\hat{a}^{\,j}_{x,g}+\hat{b}^{\,j}_{x,g}\hat{\kappa}^{\,j}_{t,g}\right).\end{equation*}

Algorithm 2:

Estimation of the GBLL Model

To ensure consistency in seasonal mortality patterns across countries, particularly between the hemispheres, we address the phase reversal in seasonal cycles observed in weekly mortality data. Specifically, winter-related mortality peaks occur in opposite calendar months across hemispheres. During model estimation, we therefore apply a temporary monotonic reciprocal transformation (i.e., $\tilde{m}_{x,t}^{\,j*}=1/\tilde{m}_{x,t}^{\,j}$ where $\tilde{m}_{x,t}^{\,j*}$ and $\tilde{m}_{x,t}^{\,j}$ are the transformed and raw observed mortality rates, respectively) to the Southern Hemisphere mortality rates (Australia and New Zealand) to align the phase of their seasonal patterns with those of Northern Hemisphere countries. This transformation is used solely for estimation purposes to facilitate the identification of a common seasonal component within the multi-population framework. After forecasting, all series are transformed back to their original scale (i.e., $\hat{m}_{x,t}^{\,j}=1/\hat{m}_{x,t}^{\,j*}$ where $\hat{m}_{x,t}^{\,j}$ is the fitted value of the mortality rate and $\hat{m}_{x,t}^{\,j*}$ is the fitted value of the transformed mortality rate), ensuring that the resulting mortality forecasts remain fully interpretable.

Forecasting mortality rates primarily relies on projecting both the common and country-specific time trend component $\kappa_t$ . In classical mortality modeling, $\kappa_t$ is often assumed to follow a random walk with drift. For high-frequency weekly mortality data, a natural extension is to employ a seasonal-ARIMA (SARIMA) model to account for annual seasonality. However, SARIMA models require seasonal differencing at lag 52, which effectively removes a full year of observations from the already short five-year sample (260 weeks), potentially leading to unstable parameter estimates and unreliable forecasts. Moreover, SARIMA with a long seasonal period is computationally intensive, especially in a multi-population gradient boosting framework. To overcome this limitation, we adopt a hybrid forecasting approach in which seasonal patterns are modeled deterministically using Fourier regressors (Serfling, Reference Serfling1963), while the residual stochastic dynamics are captured by a non-seasonal ARIMA(p,d,q) process. This strategy allows us to retain the full sample for estimation while flexibly capturing both annual and biannual seasonal cycles. Specifically, $\kappa_t$ is modeled as

\begin{equation*}\kappa_t = \beta_1 \sin\!\big(2\pi w(t)\big) + \beta_2 \cos\!\big(2\pi w(t)\big) + \beta_3 \sin\!\big(4\pi w(t)\big) + \beta_4 \cos\!\big(4\pi w(t)\big) + \text{ARIMA}(p,d,q),\end{equation*}

where $w(t)\in[0,1)$ denotes the fractional part of the calendar year corresponding to time t. There are four Fourier regressors, including both sine and cosine functions to capture both the annual and biannual cycles. The statistical significance of the Fourier components is assessed using a Wald test on the corresponding regression coefficients (Wald, Reference Wald1943). Specifically, for each Fourier term, the null hypothesis that its coefficient is equal to zero is tested using the associated z-statistic $z_i = \hat{\beta}_i / \mathrm{SE}(\hat{\beta}_i)$ , for $i = 1,2,3,4$ , along with the corresponding p-value obtained from the coeftest function in R. Under the large-sample approximation, the test statistic follows $Z \sim \mathcal{N}(0,1)$ . A Fourier regressor is deemed statistically significant at the 5% level if its p-value is less than 0.05. If none of the Fourier components satisfy this criterion, the forecasting procedure relies solely on the ARIMA(p,d,q) model. Therefore, standing on the final observed time T, the final h-step-ahead forecasts of the log mortality rate is given by: for age group $x\in\{1,...,N\}$ , $j\in\{1,...,J\}$ ,

\begin{equation*}\log\!(\hat{m}^{\,j}_{x,T+h}) = \sum_{g=1}^{l} \gamma_g\left(\hat{a}^{\,p}_{x,g}+\hat{b}^{\,p}_{x,g}\hat{\kappa}^{\,p}_{T+h|T,g}+\hat{a}^{\,j}_{x,g}+\hat{b}^{\,j}_{x,g}\hat{\kappa}^{\,j}_{T+h|T,g}\right).\end{equation*}

There are alternative choices for each component of the GBLL model, including the loss function, stopping criterion, models used to fit and forecast the time trend $\kappa_t$ , and seasonal adjustments for countries in the Southern Hemisphere. Detailed descriptions of these alternatives are provided in Section A of the Supplementary Material. Nevertheless, our main conclusions remain unchanged, and the proposed GBLL model is shown to be insensitive to the choice of individual components.

4.1. Fitting performance

Empirically, the GBLL algorithm terminates due to condition (b), typically after around eight iterations. The maximum in-sample fitting error across all horizons, countries, and age groups is below 0.001, indicating that the model has captured all underlying significant patterns and meaningful variations in the data. This adequacy can be further assessed by examining the number of residual series exhibiting white-noise behavior across age groups and countries using the Ljung–Box test at a 5% significance level (Ljung and Box, Reference Ljung and Box1978).

For each expanding window, we count the number of residual series identified as white noise (out of 120, corresponding to 30 countries and 4 age groups). The average number of white-noise residuals across the 10 expanding windows show that the LL and HBY models produce 54 and 15 white-noise residuals, respectively, while the proposed GBLL model (with eight iterations) produces 73. This indicates that the GBLL model captures nearly all significant and meaningful patterns in the mortality data, as also reflected by the small maximum in-sample fitting error (below 0.001). In contrast, the HBY model fails to capture some structures, particularly for the 65–74 and 75–84 years of age groups, while the LL model performs better but still leaves more residual structure unexplained. The higher number of white-noise residuals in the GBLL model, together with the relatively small fitting error, demonstrates its superior ability to extract the underlying trends and seasonality across countries and age groups. It is worth noting that our comparison focuses on the HBY model with order 6; alternative specifications may perform differently, but exploring the full range of configurations is beyond the scope of this paper.

Figure 2.

The residual plot of Canada under the LL model.

Figure 2. Long description

The image contains four separate line graphs, each representing residual plots for different age groups in Canada under the LL model. The graphs are arranged in a 2x2 grid. The top-left graph shows residuals for the age group 15-64 years old, with the x-axis labeled ’Date’ ranging from 2015 to 2019 and the y-axis labeled ’Residual’ ranging from -2e-04 to 3e-04. The top-right graph shows residuals for the age group 65-74 years old, with the x-axis labeled ’Date’ ranging from 2015 to 2019 and the y-axis labeled ’Residual’ ranging from -5e-04 to 5e-04. The bottom-left graph shows residuals for the age group 75-84 years old, with the x-axis labeled ’Date’ ranging from 2015 to 2019 and the y-axis labeled ’Residual’ ranging from -0.003 to 0.002. The bottom-right graph shows residuals for the age group 85+ years old, with the x-axis labeled ’Date’ ranging from 2015 to 2019 and the y-axis labeled ’Residual’ ranging from -0.010 to 0.010. Each graph displays fluctuations in residuals over time for the respective age groups.

Figure 3.

The residual plot of Canada under the HBY model.

Figure 3. Long description

The residual plot of Canada under the HBY model consists of four separate line graphs, each representing different age groups: 15-64 years old, 65-74 years old, 75-84 years old, and 85+ years old. The x-axis for all graphs indicates the date ranging from 2015 to 2019, while the y-axis represents the residual values. The graph for the 15-64 years old group shows fluctuations in residuals with notable peaks and troughs throughout the period. The 65-74 years old group exhibits similar patterns with some distinct peaks. The 75-84 years old group’s residuals show a significant peak around 2017 and subsequent fluctuations. The 85+ years old group’s residuals demonstrate a gradual increase with periodic peaks. Each graph highlights the short-term mortality fluctuations and seasonal dynamics within the specified age groups. All values are approximated.

Looking closely at those residuals that are not white noise, there remain either seasonal cyclical patterns or clear downward/upward sloping trends or both. This means that, given the complexity and high dimensions of the dataset, it is not suitable to use the existing methods, as they fail to capture all the patterns. For instance, in the case of Canada, residual plots from both the LL and HBY models reveal remaining seasonal cycles and trends. Figure 2 shows clear seasonal cycles for the 65–74 years old and a downward sloping trend for the 75–84 years old, under the LL model. Similarly, Figure 3 gives the residual plot of Canada under the HBY model, and there are clear seasonal cyclical patterns for all the four age groups and a slightly decreasing trend for the 65–74 years old. Notably, the age group of 15–64 years exhibits an inverse seasonal pattern. This may be due to the model imposing seasonal structure to accommodate for other age groups, even though the raw data show little seasonal variation. Additionally, the residuals for individuals aged 85 years and over are approximately around $\pm 0.01$ , clearly indicating that both the LL and HBY models struggle to capture the underlying patterns for this age group. In contrast, the residuals across all four age groups under the GBLL model exhibit white-noise behaviors, indicating improved model adequacy. Additional visualizations of in-sample fits for the GBLL model and the benchmark models for selected countries are provided in Section C.1 of the Supplementary Material.

Finally, Figure 4 illustrates the contribution of each component at different stages of the model fitting process. We present the decomposition of the GBLL model for the first two iterations here, focusing on the common trends. The country-specific components breakdown are provided in Appendix E. The first column (left) displays, from top to bottom, the estimated age-specific intercepts ( $\hat{a}_x$ ), age loadings ( $\hat{b}_x$ ), and the time-varying index ( $\hat{\kappa}_t$ ) from the first iteration. The second column (right) presents the corresponding components from the second iteration.

Figure 4.

Components of the common trends under the GBLL model: first iteration (left) and second iteration (right).

Figure 4. Long description

The image contains six line graphs arranged in two columns. The left column represents the first iteration, and the right column represents the second iteration. Each column contains three graphs. The top graphs show the parameter alpha (a) for different age groups: 15-64, 65-74, 75-84, and 85+. The middle graphs display the parameter beta (b) for the same age groups. The bottom graphs illustrate the parameter kappa (k) over time from 2015 to 2019. In the first iteration, alpha and beta increase with age, while kappa fluctuates over time. In the second iteration, alpha and beta decrease with age, and kappa shows smaller fluctuations over time. All values are approximated.

• • For the age-specific intercepts $\hat{a}_x$ , the first iteration captures a generally increasing trend. It highlights the vulnerability of older people, resulting in higher average mortality rates. In the second iteration, a predominantly decreasing pattern is observed. Importantly, the magnitude of these adjustments is on the order of $10^{-6}$ , suggesting that the primary age-related mortality pattern is effectively captured in the first iteration.

• • For the age loadings $\hat{b}_x$ , the first iteration exhibits an upward trend, indicating the greater exposure experienced by elderly individuals to higher mortality rates. In the second iteration, a declining trend is observed. This may be related to the original smoother and more pronounced seasonal patterns typically exhibited by older age groups. As these dominant effects are captured in the first iteration, only a small portion of the residual variation remains to be explained. Additionally, some spikes exist in the 15–64 and 85+ age groups, indicating that this iteration primarily captures residual effects concentrated in the extreme age cohorts.

• • For the common time trend component $\hat{\kappa}_t$ , the first iteration effectively captures the dominant seasonal cycles. The second iteration provides a refinement by identifying smaller-scale cyclical patterns with short-term fluctuations and a mild downward trend, indicative of gradual improvements in mortality rates over time. A prominent W-shaped seasonal pattern is evident in the first iteration across all years, arising from the aggregation of mortality dynamics across all countries. This aggregation enhances structural features that may be weakly present but commonly shared across populations, thereby making the W-shape (with dual seasonal peaks) more pronounced. Additionally, an M-shaped pattern is discernible, particularly around early 2015 and 2018. This may be attributed to the reciprocal transformation applied to the Southern Hemisphere countries, where a secondary mortality peak, likely associated with heatwaves, typically occurs at the beginning of the calendar year.

4.2. Forecast performance with expanding windows

To evaluate the forecast performance, we adopt the Mean Absolute Percentage Error (MAPE) as the error measure to quantify the average magnitude of the percentage deviation between forecasted and actual values. MAPE is scale-independent, allowing direct comparison of forecast accuracy across countries with different mortality levels. It has been widely used in the mortality forecasting literature as a standard metric for evaluating predictive performance (see e.g., McNown and Rogers, Reference McNown and Rogers1992; Neves et al., Reference Neves, Fernandes and Hoeltgebaum2017; Li et al., Reference Li, Li, Lu and Panagiotelis2019; Tsai and Cheng, Reference Tsai and Cheng2021; Qiao et al., Reference Qiao, Wang and Zhu2024). We first define the set $H= \{4, 9, 13, 17, 22, 26, 30, 35, 39, 43, 48, 52\}$ , which corresponds to the cumulative number of weeks until the end of each month, thereby facilitating the conversion from weekly to monthly forecast horizons. For month $h\in\{1,...,12\}$ , expanding window $r\in\{0,...,9\}$ , the MAPE for h-step-ahead forecasts across all age groups, countries, and expanding windows is given by:

\begin{equation*}\mathrm{MAPE}_h=\frac{1}{4\times H_{(h)} \times 30 \times 10}\sum_{x=1}^4\sum_{u=1}^{H_{(h)}}\sum_{j=1}^{30}\sum_{r=0}^{9}\frac{|\hat{m}^{\,j}_{x,169+u+H_{(r)}}-\tilde{m}^{\,j}_{x,169+u+H_{(r)}}|}{\tilde{m}^{\,j}_{x,169+u+H_{(r)}}},\end{equation*}

where $H_{(i)}$ represents the i-th order statistic of the set H with $H_{(0)}=0$ and 169 is the number of weeks in the first training sample. An additional scale-invariant error metric, the Mean Absolute Scaled Error (MASE), is also employed and leads to comparable conclusions to those obtained using MAPE. Full results and discussions are provided in Section B of the Supplementary Material.

Table 1.

Mean MAPE of out-of-sample forecasts across 30 countries ( $\times 100$ scale).

Table 1. Long description

The table presents the mean Mean Absolute Percentage Error (MAPE) of out-of-sample forecasts across 30 countries for various models over a forecast horizon ranging from 1 to 12 months. The models compared include LL, HBY, GBLL, VAR, and GVAR. The table has 12 rows and 5 columns. The columns are labeled LL, HBY, GBLL, VAR, and GVAR, representing different models. Each row corresponds to a different forecast horizon (h) from 1 to 12 months. The values in the table are the mean MAPE scaled by a factor of 100. The lowest mean MAPE for each horizon, indicating the best-performing model, is highlighted in bold. The GBLL model consistently shows the lowest mean MAPE across all forecast horizons, indicating superior performance compared to the LL, HBY, VAR, and GVAR models. This suggests that the GBLL model, leveraging the iterative structure of the gradient boosting framework, is more effective in capturing complex patterns in the data. In contrast, the VAR and GVAR models, which are more generic and require differencing to ensure stationarity, exhibit comparatively weaker performance. The results highlight the strengths of factor-based models like GBLL in forecasting long-run mortality trends, while VAR-type models are better suited for analyzing short-run dynamics.

In addition to the factor-based LL and HBY models, we include further comparisons of forecast performance using VAR models. The specifications of both the single-population VAR and the multi-population GVAR models are provided in Appendix D. Table 1 presents the mean MAPE of out-of-sample forecasts over the forecast horizon h from 1 to 12 months (scaled by a factor of 100). For each horizon, the lowest mean MAPE, indicating the best-performing model, is highlighted in bold. The results indicate that, based on the STMF dataset, the proposed GBLL model outperforms the two benchmark factor models and the two VAR models across all forecast horizons. This superior performance is achieved by the iterative structure of the gradient boosting framework, which enables the model to sequentially identify and incorporate multiple layers of complex patterns. The comparatively weaker performance of VAR and GVAR models can be attributed to their generic time series nature and the differencing required to ensure stationarity. In contrast to factor models (LL, HBY, and GBLL), which are specifically designed to capture and forecast long-run mortality trends, VAR-type models are typically estimated on differenced series and are therefore better suited to analyzing short-run dynamics, such as spillover effects and shock transmission as characterized by impulse response functions. Moreover, the relatively short time span of mortality data exacerbates the curse of dimensionality in VAR and GVAR models, leading to increased estimation uncertainty. Accordingly, subsequent detailed comparisons of model performance will focus only on the class of factor-based models.

To gain a better understanding of the enhanced forecast accuracy achieved by the proposed GBLL model, Figure 5 presents the improvement in MAPE across countries and forecast horizons. Compared to the standard LL model, the GBLL model enhances forecast accuracy. Substantial gains are observed in countries such as South Korea and Taiwan, while Lithuania is slightly better under the LL model. Relative to the HBY model, GBLL again shows superior performance in nearly all cases, except for Hungary in the first two forecast months. Overall, the GBLL model achieves the highest accuracy across the three models evaluated. Additional visualizations of out-of-sample forecasts for the GBLL model and the benchmark models for selected countries are provided in Section C.2 of the Supplementary Material. Interestingly, the LL model outperforms the HBY model despite its simpler structure. A likely explanation is that the use of only four age groups, as opposed to single-year-of-age data for which the HBY model is designed, reduces the number of principal components needed to explain the variation in data. Consequently, the additional principal components may contribute little beyond the first. Therefore, given the complexity and heterogeneity of mortality data across 30 countries, simpler or progressive approaches such as the gradient boosting framework, may be better suited to capturing shared and country-specific patterns than more complex one-step models. We have conducted an additional analysis to evaluate forecast performance across age groups, with a detailed discussion provided in Section D of the Supplementary Material. Overall, the proposed GBLL model demonstrates superior forecast accuracy for most age groups, particularly among younger cohorts, while the simpler LL model performs comparably for the oldest age group, likely due to its smoother mortality patterns.

Figure 5.

Improvement in MAPE from LL to GBLL (left) and from HBY to GBLL (right).

Figure 5. Long description

The heat map consists of two panels, each representing a different comparison of improvement in mean absolute percentage error (MAPE). The left panel shows the improvement from LL to GBLL, while the right panel shows the improvement from HBY to GBLL. Each panel is a grid with countries listed along the x-axis and horizons along the y-axis. The color scale ranges from light to dark purple, indicating increasing levels of improvement. Darker colors represent higher improvements. The left panel shows a general trend of improvement across most countries and horizons, with some countries showing higher improvements at specific horizons. The right panel also shows a similar trend but with different patterns of improvement. Notable improvements are observed in specific countries at certain horizons, indicating the effectiveness of the GBLL model in reducing MAPE compared to LL and HBY models.

In addition to evaluating point forecast accuracy, we assess each model’s ability to quantify predictive uncertainty via marginal calibration of 95% prediction intervals (Czado et al., Reference Czado, Gneiting and Held2009). For each model, 100 forecast paths for $\kappa_t$ are simulated, and prediction intervals for $\hat{m}^{\,j}_{x,t}$ are constructed from the empirical 2.5th and 97.5th percentiles, denoted $\hat{p}^{\,j}_{2.5,x,t}$ and $\hat{p}^{\,j}_{97.5,x,t}$ . Marginal calibration evaluates whether the nominal probability of the interval aligns with the empirical coverage, that is, the proportion of observed outcomes that fall within the interval. Formally, for week $h_w \in \{1,...,52\}$ , expanding window $r \in \{0,...,9\}$ , age group $x \in \{1,...,4\}$ , and country $j \in \{1,...,30\}$ , the marginal calibration $f_{[\hat{p}^{\,j}_{2.5,x,t},\hat{p}^{\,j}_{97.5,x,t}]}$ is given by:

\begin{align*}f_{[\hat{p}^{\,j}_{2.5,x,t},\hat{p}^{\,j}_{97.5,x,t}]} = \frac{1}{52 \times 10 \times 4 \times 30} \sum_{h_w=1}^{52} \sum_{r=0}^{9} \sum_{x=1}^{4} \sum_{j=1}^{30} \unicode{x1D7D9}\Big(\hat{p}^{\,j}_{2.5,x,169+h_w+H_{(r)}} \le \tilde{m}^{\,j}_{x,169+h_w+H_{(r)}} \le \hat{p}^{\,j}_{97.5,x,169+h_w+H_{(r)}} \Big), \\[-10pt] \end{align*}

where $\unicode{x1D7D9}$ denotes the indicator function. Mean empirical marginal calibration rates show that the LL, HBY, and GBLL models achieve 90%, 94%, and 96% coverage, respectively. These results indicate that the GBLL model is the most properly calibrated, producing prediction intervals that reliably capture realized mortality rates and providing improved uncertainty quantification alongside superior point forecast performance.

5.1. Setup

Clustering typically involves four key components: the data input, clustering algorithm, distance measure, and the optimal number of clusters (Yin et al., Reference Yin, Aryani, Petrie, Nambissan, Astudillo and Cao2024). In this paper, we first extract the time trend component $\kappa_t$ for each country using the Lee–Carter model. We then consider three different variations of the raw time series as the data input to investigate population similarity, while keeping the remaining three clustering components fixed.

Clustering algorithms can be broadly classified into several categories, including partitional, hierarchical, density-based, grid-based, and model-based approaches (Ghosal et al., Reference Ghosal, Nandy, Das, Goswami and Panday2020). In this study, we adopt the K-means algorithm, a representative method within the partitional clustering framework, which minimizes within-cluster sum of squares through an iterative process of centroid updating and point re-assignment (Lloyd, Reference Lloyd1982). The accompanying method to find the optimal number of clusters is the elbow method, which plots the inertia (within-cluster sum of squares) against a range of values for the number of clusters (Thorndike, Reference Thorndike1953). Inertia is calculated as:

\begin{equation*}\text{Inertia} = \sum_{q=1}^{Q} \sum_{\mathbf{x} \in C_q} (\mathbf{x}-\boldsymbol{\mu}_q)'(\mathbf{x}-\boldsymbol{\mu}_q),\end{equation*}

where Q is the total number of clusters, $\mathbf{x}$ is an observation for each country assigned to cluster $C_q$ , and $\boldsymbol{\mu}_q$ is the centroid of cluster $C_q$ . The final optimal number is identified at the “elbow” point of the curve, where the rate of decrease in inertia slows significantly, indicating diminishing returns in clustering quality with additional clusters. Finally, the distance measure used to define similarity between data points is the Euclidean distance metric. It calculates the straight-line distance between two points in an n-dimensional space.

5.1.1. Method 1: Clustering based on raw time trend series

We begin by analyzing the raw time series of the country-specific $\kappa_t$ as the data input. This approach identifies three distinct clusters, as detailed in Table 2. Figure 6 presents the mean value of the time trends for each cluster. Notably, Cluster 3, which comprises only two countries from the Southern Hemisphere, exhibits seasonal patterns that are opposite to those of the other clusters. Cluster 2 consists exclusively of European countries and displays more pronounced seasonal peaks relative to Cluster 1, highlighting its higher mortality rates in winter.

Table 2.

List of countries in each cluster under Method 1.

Table 2. Long description

The table categorizes countries into three clusters based on a specific method. It has three rows and three columns. The columns are labeled Cluster, Number of countries, and Country name. Cluster 1 includes 11 countries: Canada, Denmark, England and Wales, Finland, Netherlands, Norway, Scotland, South Korea, Sweden, Taiwan, and the USA. Cluster 2 includes 17 countries: Austria, Belgium, Bulgaria, Croatia, Czech Republic, France, Germany, Greece, Hungary, Italy, Lithuania, Poland, Portugal, Slovakia, Slovenia, Spain, and Switzerland. Cluster 3 includes 2 countries: Australia and New Zealand.

Figure 6.

Mean value of the time trends in each cluster.

Figure 6. Long description

The line graph illustrates the mean value of time trends for three clusters labeled C 1 mean, C 2 mean, and C 3 mean. The x axis represents the date ranging from 2015 to 2020, while the y axis represents the mean value ranging from negative 0.5 to 1.0. The graph shows three distinct lines, each representing a different cluster. The C 1 mean line is red, the C 2 mean line is green, and the C 3 mean line is blue. The lines fluctuate over time, indicating varying trends for each cluster. The C 2 mean line shows the highest peak around 2017, while the C 3 mean line exhibits more frequent fluctuations. All values are approximated.

It is noteworthy that Cluster 2 exhibits a pronounced spike during the winter of 2017, which can be attributed to an extreme cold-weather event, affecting most of Europe in January. Temperatures remain significantly below seasonal averages for several days, particularly impacting southern and southeastern Europe. Greece, for instance, records its highest January mortality rate in 2017 compared to any other January between 2013 and 2018, with individuals aged over 80 years being most affected (Kostopoulou, Reference Kostopoulou2023). Italy also experiences abnormal snowfall due to cold air moving across the relatively warmer Adriatic Sea, picking up moisture and generating intense precipitation (Murtas and Russo, Reference Murtas and Russo2019). Demirtaş (Reference Demirtas2022) further investigates the broader effects of severe snowstorms and cold temperatures across southeastern Europe. A similar pattern is observed in 2019, with another cold wave, triggered by a Sudden Stratospheric Warming (SSW) event. It results in widespread snowfall and freezing temperatures across central and southern Europe, contributing to increased mortality and significant societal disruptions (Knight et al., Reference Knight, Scaife, Bett, Collier, Dunstone, Gordon, Hardiman, Hermanson, Ineson and Kay2021).

5.1.2. Method 2: Clustering based on estimated trend slopes

Building on the previous analysis of the raw time series, we now focus on the trend reflecting mortality improvement over the years. We first apply Seasonal and Trend decomposition using LOESS (STL) to the country-specific $\kappa_t$ series (Cleveland et al., Reference Cleveland, Cleveland, McRae and Terpenning1990). This method uses locally estimated scatterplot smoothing (LOESS) to decompose each time series into three distinct components: a long-term trend, a recurring seasonal pattern, and a residual representing irregular fluctuations. To classify countries based on their trend dynamics, we compute the slope of each decomposed long-term trend using linear regression. This procedure identifies three distinct clusters, as detailed in Table 3. Figure 7 displays the distribution of trend slopes by clusters. Countries in Cluster 1 exhibit the most substantial mortality improvement, whereas Cluster 2 countries have moderate improvement and Cluster 3 countries are more stable over the five-year period.

Table 3.

List of countries in each cluster under Method 2.

Table 3. Long description

The table presents a list of countries categorized into three clusters based on their mortality improvement trends over a five-year period. The table has three columns: Cluster, Number of countries, and Country name. Cluster 1 includes five countries: Croatia, Lithuania, Slovakia, South Korea, and Sweden. Cluster 2 consists of twelve countries: Australia, Austria, Belgium, Czech Republic, Hungary, Italy, Norway, Scotland, Slovenia, Spain, Switzerland, and Taiwan. Cluster 3 comprises thirteen countries: Bulgaria, Canada, Denmark, England and Wales, Finland, France, Germany, Greece, Netherlands, New Zealand, Poland, Portugal, and the USA. The clusters are determined by the slope of the long-term trend in mortality improvement, with Cluster 1 showing the most substantial improvement, Cluster 2 moderate improvement, and Cluster 3 stability.

Figure 7.

Trend slopes by clusters.

Figure 7. Long description

A scatter plot titled ’Trend slopes by clusters’ displays data points representing trend slopes for various countries. The x-axis represents the trend slope, ranging from approximately -0.0020 to 0. The y-axis does not have a specific label but appears to categorize the data points by country. The data points are color-coded into three clusters: red, blue, and green. The red cluster includes countries such as Korea, Sweden, and Lithuania, with trend slopes around -0.0020. The blue cluster includes countries like Spain, Austria, and Belgium, with trend slopes around -0.0010. The green cluster includes countries such as the United States, Germany, and Poland, with trend slopes closer to 0. Each data point is labeled with the country’s abbreviation. The plot shows that countries in the red cluster have the steepest negative trend slopes, while those in the green cluster have trend slopes closest to zero. All values are approximated.

All countries have negative trend slopes, indicating the mortality improvement across five years. The rapid mortality improvement observed in Cluster 1 may be attributed to factors such as healthcare system reforms and advances in medical technology and treatment (see e.g., OECD, 2021; Pekarcikova, Reference Pekarcikova2024). Given that countries like Croatia, Lithuania, and Slovakia have historically exhibited relatively higher mortality rates, these advancements suggest a particularly substantial improvement in health outcomes for these populations, resulting in higher trend slopes in absolute values. In Clusters 2 and 3, some geographical patterns are evident, including the grouping of Canada and the USA, as well as several neighboring European countries. These similarities likely reflect shared characteristics such as public health infrastructure, climatic conditions, and socioeconomic contexts, which may contribute to comparable rates of mortality improvement over the period from 2015 to 2019.

Table 4.

List of countries in each cluster under Method 3.

Figure 8.

Clusters based on min-max scaled trend slopes and seasonal strength.

Long description

A scatter plot illustrates clusters based on min-max scaled trend slopes and seasonal strength. The x-axis represents trend slope values ranging from 0 to 1, while the y-axis represents seasonal strength values also ranging from 0 to 1. The plot contains several data points, each representing a different country. The data points are color-coded into three clusters: red, green, and blue. The red cluster is positioned towards the lower left, indicating lower trend slopes and seasonal strength. The green cluster is spread across the middle, showing moderate trend slopes and seasonal strength. The blue cluster is located towards the upper right, indicating higher trend slopes and seasonal strength. Notable countries in the blue cluster include USA, Canada, and Australia, which have high seasonal strength and trend slopes. Countries in the red cluster, such as South Korea and Slovakia, exhibit low values in both metrics. The green cluster includes countries like Germany, Poland, and Greece, with varying trend slopes and seasonal strengths. All values are approximated.

5.2. Forecast performance with expanding windows

To mitigate randomness in the sampling procedure, we again employ a 10-fold expanding window approach. Table 5 summarizes the forecast performance across the clustering methods described above in Section 5.1. Among the three clustering approaches, Method 1, which uses raw time trend series, yields the poorest performance, likely due to its higher dimensionality and associated noise. In contrast, Method 2, which relies solely on the slopes of country-specific time trends, achieves the best forecast accuracy. This result suggests that utilizing a more distinct and parsimonious clustering criterion may facilitate the grouping of structurally similar countries, thereby enhancing predictive performance. In this experiment, regardless of the clustering strategy employed, based on the available data, our proposed GBLL model outperforms both the LL and HBY models, with the latter exhibiting the weakest forecast accuracy.

Table 5.

Mean MAPE of out-of-sample forecasts across 30 countries with clustering ( $\times 100$ scale).

Table 5. Long description

The table presents the mean mean absolute percentage error (MAPE) of out-of-sample forecasts across 30 countries using three different clustering methods and a whole method. The table has 12 rows and 13 columns, including headers. The columns are labeled as Method 1, Method 2, Method 3, and Whole, each with sub-columns for LL, HBY, and GBLL. The rows are labeled with the variable h ranging from 1 to 12. Notable trends include Method 2 consistently showing lower MAPE values compared to Method 1 and Method 3, indicating better forecast accuracy. The GBLL model generally outperforms the LL and HBY models across all methods. The data suggests that utilizing distinct and parsimonious clustering criteria enhances predictive performance.

To further assess the impact of clustering, we compare the results with the case without clustering, as presented in Table 1. The overall performance remains comparable; however, Method 2 is the only clustering strategy that improves forecast accuracy across all models and forecast horizons relative to the baseline case without clustering. The varying performance may be attributed to the heterogeneity in the rate of mortality improvement and the degree of seasonality across individual countries, particularly given that weekly mortality data tend to exhibit higher-frequency fluctuations and varying levels of noise. These findings suggest that, within this dataset of 30 countries, trend slopes serve as the most effective factor for grouping countries with similar mortality experience together. Additionally, we include a reference column reporting the forecast MAPE of the GBLL model without any clustering. Notably, this unclustered version of the GBLL model still outperforms all variants of the LL and HBY models under each of the three clustering approaches and across all forecast horizons. Based on the available data, the empirical evidence highlights that the GBLL framework achieves superior forecast performance without the need for clustering, surpassing baseline models that rely on clustering to enhance coherence and forecast accuracy.

Focusing on the HBY model, which emphasizes the need for grouping countries with similar characteristics to enhance forecast performance, the results lend support to this premise. All three clustering methods improve forecast accuracy compared to the unclustered baseline, with the exception of the first month under Method 1. However, identifying an optimal clustering strategy that yields substantial improvement on all models over the baseline case without clustering remains challenging. This may be attributed to the limited span of data with only five years, during which there is minimal divergence in mortality trends, resulting in relatively homogeneous raw time trend series across countries. These observations underscore the importance of a modeling framework that is both adaptive and insensitive, regardless of clustering choices. The LL model demonstrates relatively similar performance across various clustering methods and may be regarded as a stable, though less flexible, approach. However, this stability comes at the cost of lower predictive accuracy. Conversely, the HBY model exhibits substantial sensitivity to country groupings and continues to deliver the weakest forecast performance overall. In contrast, the proposed GBLL model distinguishes itself by maintaining high forecast accuracy across all clustering strategies. Its performance is largely unaffected by the social, economic, or environmental heterogeneity of the countries, highlighting its insensitivity and superior forecasting capability.