4.1. Hierarchical clustering with the dynamic time warping

Li et al. (Reference Li, Li, Lu and Panagiotelis2019) applied Ward’s hierarchical clustering to normalized age-specific COD mortality rates to identify hierarchical groupings. However, clustering directly on age-specific mortality rates can be inefficient. COD rates are often volatile and noisy, particularly at ages with small case-specific death counts. Hence, clustering such data does not effectively reduce this noise. Moreover, without decomposition, clustering on age–period matrices conflates age and period effects, making interpretation difficult. Finally, Ward’s method employs the Frobenius norm distance, which implicitly weights each dimension by its raw variance. As a result, even after normalization, dimensions with larger fluctuations may dominate, obscuring more interpretable clustering structures.

Figure 4.

LC-corrected parameter estimates for different causes of death.

To resolve these issues, we follow Varga (Reference Varga2025) and perform the dynamic time warping (DTW) method to obtain hierarchical clusters to $b_i(x)$ and $\kappa_{i,t}$ obtained as in Section 2 of the Supplementary Material. DTW is a shape-based distance measure originally developed by Bellman and Kalaba (Reference Bellman and Kalaba1959) for speech recognition applications. It computes an improved alignment between two time-series by allowing nonlinear “warping” along the time axis. In our context, taking $\kappa$ ’s as an example, DTW seeks a warping path $\psi = \{(h_1,v_1), \dots, (h_L,v_L)\}$ through the grid $\{1, 2, \dots,T\}^2$ that aligns each year h of $\kappa_{i,\cdot}$ with a (possibly different) year v of $\kappa_{j,\cdot}$ .

To build this alignment, we follow Varga (Reference Varga2025) and form the local cost matrix (LCM) using the $l_2$ norm:

\begin{equation*}\mathrm{LCM}_{h,v}^\kappa\;=\;\bigl(\kappa_{i,h}^* - \kappa_{j,v}^*\bigr)^2,\end{equation*}

where $\kappa_{i,h}^*$ and $\kappa_{j,v}^*$ are normalized entries from $\kappa_{i,t}$ and $\kappa_{j,t}$ , respectively. Each pair measures the pointwise discrepancy between the year h of one series and the year v of the other.

The DTW distance is then defined as the minimal cost route through this matrix. Let $\psi$ be any admissible path from (1, 1) to (T, T). Assign to each step $(h_\ell,v_\ell)\in\psi$ a weight $\eta_\ell$ (set to unity) and let $\eta = \sum_{\ell=1}^{L} \eta_\ell$ , where L is the length of the warping path $\psi$ . The DTW distance is

\begin{equation*}{d}_{i,j}^\kappa\;=\;\min_{\psi}\Biggl\{\sum_{\ell=1}^L\frac{\eta_\ell \,\mathrm{LCM}_{(h_\ell,v_\ell)}}{\eta}\Biggr\}^{1/2}.\end{equation*}

This construction guarantees that the first and last years match exactly, while intermediate years can “stretch” or “compress” to reveal the strongest common trends, even if the improvement of one cause leads or lags the other.

To equally weight age and temporal dynamics in the cluster analysis, we define a composite DTW distance between two causes i and j as the sum of the univariate DTW distances. That is,

\begin{equation*}{d}_{i,j}={d}_{i,j}^b+{d}_{i,j}^\kappa.\end{equation*}

The resulting DTW distance matrix can be computed via dynamic programming as in Giorgino (Reference Giorgino2009). This matrix is then used as input to the hierarchical clustering as in Sardá-Espinosa (Reference Sardá-Espinosa2019). The final clusters are plotted in Figure 5 and described in Table 2.

Table 2.

Clusters and hierarchical levels by clustering approach. The slope and trend-polarity clusterings are introduced in Section 4.2

Note: CA = cancer, DM = diabetes, EC = external causes, Flu = influenza, MD = mental disorder, NE = nephritis, VA = vascular disease.

The DTW-based clustering analysis reveals two primary hierarchical groups. The first cluster – diabetes, mental disorders, and nephritis – exhibits increasing trends in their $\kappa$ indices (Figure 3). The second cluster – vascular diseases, cancer, and external causes – demonstrates steady improvements over time, as shown in Figure 3. Within the first cluster, mental disorders and nephritis form the closest subcluster, while in the second cluster, cancer and external causes are identified as the closest pair. These groups align with the similarities observed in their estimated $b_i(x)$ and $\kappa_{i,t}$ profiles (Figure 3). This underscores the effectiveness of our DTW-based method in capturing shared mortality dynamics.

Table 3.

Sample mean of $\triangle\kappa_{i,t}$ for each COD.

Figure 5.

Hierarchical clustering using custom DTW distance.

4.2. Alternative hierarchical clustering

In addition to the clusters identified in Figure 5, this section examines alternative clustering strategies. Unlike the DTW approach, which incorporates both age and temporal dimensions, the first alternative clustering method considers only the similarity of temporal dynamics.

Specifically, we computed the sample mean of the differenced $\kappa_{i,t}$ (i.e., $\triangle\kappa_{i,t} = \kappa_{i,t} - \kappa_{i,t-1}$ ) using the corrected LC estimates. These means, shown in Table 3, represent the average annual rate of change in the period index. Consistent with the DTW results, two broad clusters emerge: one with positive mean values (indicating increasing trends) and one with negative mean values (indicating declining trends). Within the former, diabetes (0.23) and mental disorders (0.20) show the closest mean rates. Within the latter, cancer and external causes share similar means around $-1$ , while influenza and vascular diseases form a sub-cluster with means roughly between $-2$ and $-3$ . Based on these similarities, we construct the resulting hierarchical structure, shown in Table 2 (named slope clustering).

Finally, we consider a simplified clustering strategy that leverages only the overall trend. Two groups are identified: the first comprises causes exhibiting exclusively upward trends (diabetes, mental disorders, and nephritis), and the second includes those with downward mortality trends (cancer, influenza, external causes, and vascular diseases). No further hierarchical levels are defined. This two-tier structure (named trend-polarity clustering), together with the more granular cluster configurations explained above, is presented in Table 2 (named trend-polarity clustering).

5.1. Re-estimated LC parameters with cluster-dependent reconciliation

In Section 2 of the Supplementary Materials, corrected LC parameters were introduced to incorporate reconciliation information at the COD level. This correction improved estimation accuracy, which is essential for the clustering analysis in Section 4.1. Building on this, we propose a more flexible cluster-level reconciliation, which relaxes the requirement that fitted death counts match empirical counts for each COD.

For Cluster 2 identified by the DTW method (Figure 5), reconciliation proceeds sequentially:

1. (1) Bottom-level reconciliation. Apply a penalty term to reconcile fitted and empirical death counts for cancer and external causes.

2. (2) Intermediate-level reconciliation. Carry any remaining discrepancy from step (1) to reconcile vascular diseases alongside their own empirical and fitted death counts.

3. (3) Top-level reconciliation. Reconcile the residual difference from step (2) with the empirical and fitted death counts for Influenza.

A similar sequential reconciliation is performed independently for Cluster 1. Within each cluster, reconciliation at a higher aggregation level incorporates discrepancies from all subordinate levels, whereas reconciliation at a lower level does not incorporate information at higher levels.

To improve accuracy, the LC parameters are re-estimated using the sequential and cluster-level reconciliation procedure described above. Taking the clusters identified by the DTW method as an example, we proceed as follows:

1. (1) At the bottom level, LC parameters for mental disorders and nephritis in Cluster 1 are re-estimated by solving

   \begin{equation*}\min_{a,b,\kappa}\;\sum_{i\in\{\mathrm{MD}, \mathrm{NE}\}}\sum_{x}\sum_{t}\bigl(\ln m_{i,t}(x) - [\,a_i(x) + b_i(x)\,\kappa_{i,t}\,]\bigr)^{2}\;+\;\lambda_{1,1}\sum_{i\in\{\mathrm{MD}, \mathrm{NE}\}}\sum_{x}\sum_{t}\bigl(D_{i,t}(x) - \widehat D_{i,t}(x)\bigr)^{2},\end{equation*}
   where $\lambda_{1,1}$ is a cluster- and level-specific penalty that enforces coherence and captures shared structure among causes grouped within Cluster 1 at Level 1 of the hierarchy. Similarly, the parameters for cancer and external causes in Cluster 2 are re-estimated by solving
   \begin{equation*}\min_{a,b,\kappa}\;\sum_{i\in\{\mathrm{CA}, \mathrm{EC}\}}\sum_{x}\sum_{t}\bigl(\ln m_{i,t}(x) - [\,a_i(x) + b_i(x)\,\kappa_{i,t}\,]\bigr)^{2}\;+\;\lambda_{2,1}\sum_{i\in\{\mathrm{CA}, \mathrm{EC}\}}\sum_{x}\sum_{t}\bigl(D_{i,t}(x) - \widehat D_{i,t}(x)\bigr)^{2}.\end{equation*}

2. (2) Define the remaining imbalance for each cluster by

   \begin{align*}R_{1,1,t}(x) &= \sum_{i\in\{M,N\}}\bigl[D_{i,t}(x)-\widehat D_{i,t}(x)\bigr],\\[3pt]R_{2,1,t}(x) &= \sum_{i\in\{C,E\}}\bigl[D_{i,t}(x)-\widehat D_{i,t}(x)\bigr].\end{align*}

   Then, the parameters for diabetes in Cluster 1 are re-estimated by solving

   \begin{equation*}\min_{a,b,\kappa}\;\sum_{t,x}\bigl(\ln m_{\mathrm{DM}, t}(x) - [\,a_\mathrm{DM}(x) + b_\mathrm{DM}(x)\,\kappa_{\mathrm{DM}, t}\,]\bigr)^{2}\;+\;\lambda_{1,2}\sum_{t,x}\bigl(D_{\mathrm{DM}, t}(x) - \widehat D_{\mathrm{DM}, t}(x) - R_{1,1,t}(x)\bigr)^{2},\end{equation*}
   and the parameters for vascular diseases in Cluster 2 are re-estimated by solving
   \begin{equation*}\min_{a,b,\kappa}\;\sum_{t,x}\bigl(\ln m_{\mathrm{VA}, t}(x) - [\,a_\mathrm{VA}(x) + b_\mathrm{VA}(x)\,\kappa_{\mathrm{VA}, t}\,]\bigr)^{2}\;+\;\lambda_{2,2}\sum_{t,x}\bigl(D_{\mathrm{VA}, t}(x) - \widehat D_{\mathrm{VA}, t}(x) - R_{2,1,t}(x)\bigr)^{2}.\end{equation*}

3. (3) Only for cluster 2, define the updated imbalance

   \begin{equation*}R_{2,2,t}(x) = \bigl[D_{\mathrm{VA}, t}(x)-\widehat D_{\mathrm{VA}, t}(x)\bigr] - R_{2,1,t}(x).\end{equation*}
   The parameters for influenza are then re-estimated by solving
   \begin{equation*}\min_{a,b,\kappa}\;\sum_{t,x}\bigl(\ln m_{\mathrm{Flu}, t}(x) - [\,a_\mathrm{Flu}(x) + b_\mathrm{Flu}(x)\,\kappa_{\mathrm{Flu}, t}\,]\bigr)^{2}\;+\;\lambda_{2,3}\sum_{t,x}\bigl(D_{\mathrm{Flu}, t}(x) - \widehat D_{\mathrm{Flu}, t}(x) - R_{2,2,t}(x)\bigr)^{2}.\end{equation*}

All penalty parameters are selected using grid search, as described in Section 2 of the Supplementary Materials. The cluster-dependent LC parameters are also re-estimated using the same sequential reconciliation procedure in both the slope clustering and trend-polarity clustering frameworks. This approach more effectively exploits the hierarchical structure of the COD data than the individual-level correction in Section 2 of the Supplementary Materials. Consequently, the resulting $\kappa_{i,t}$ trajectories exhibit improved fidelity to the underlying mortality dynamics. This may provide more appropriate inputs for the subsequent estimation of cluster-dependent sparse VECM.

5.2. Sparse VECM

As discussed in Section 2, all COD log mortality rates are clearly trending over time. Hence, the estimated period indices $\kappa_{i,t}$ are nonstationary, which is well established in the mortality literature (see, e.g., Lee and Carter Reference Lee and Carter1992; Booth and Tickle Reference Booth and Tickle2008, among others). Moreover, as noted in Remark 1, the trending of the total mortality rates implies the existence of long-run equilibrium relationships, or cointegration, among the $\kappa_{i,t}$ series (Engle and Granger Reference Engle and Granger1987; Johansen Reference Johansen1991). These empirical characteristics motivate the use of a VECM to model the cluster-dependent $\kappa_{i,t}$ estimated in Section 5.1.

Let $\Delta\boldsymbol{\kappa}_t = (\Delta\kappa_{\mathrm{CA}, t}, \Delta\kappa_{\mathrm{DM}, t}, \Delta\kappa_{EC,t}, \Delta\kappa_{\mathrm{MD}, t}, \Delta\kappa_{\mathrm{NE}, t}, \Delta\kappa_{\mathrm{ Flu}, t}, \Delta\kappa_{\mathrm{VA}, t})^{\top}$ . We specify VECM(1) as

(5.1) \begin{equation}\Delta\boldsymbol{\kappa}_t = \boldsymbol{\mu} + \boldsymbol{\alpha}\,\boldsymbol{\beta} ^{\top}\,\boldsymbol{\kappa}_{t-1} + \boldsymbol{\Gamma}_1\,\Delta\boldsymbol{\kappa}_{t-1} + \boldsymbol{\varepsilon}_t, \end{equation}

where $\boldsymbol{\mu}$ denotes the intercept vector; $\boldsymbol{\beta}=(1,\beta_{2},\dots,\beta_{7})^{\top}$ defines the cointegration relationship, and we fix $\beta_{1}=1$ for the unique identification; $\boldsymbol{\alpha}$ is the adjustment (loading) vector; $\boldsymbol{\Gamma}_{1}$ is the matrix of short-run autoregressive coefficients; and $\boldsymbol{\varepsilon}_t$ is a zero-mean innovation term with covariance $\boldsymbol{\Sigma}_\varepsilon$ .

The unrestricted VECM(1) in (5.1) employs a full autoregressive matrix $\Gamma_{1}$ , implying $7\times7=49$ free parameters. With mortality data spanning only a limited period, estimating this many coefficients is often infeasible, risking overfitting, and unstable inference. To mitigate these concerns, we impose a block-structured sparsity pattern on $\Gamma_{1}$ guided by hierarchical clustering. Consistent with the justifications in Remark 3, all cross-cluster coefficients (those linking series with increasing trends to those with declining trends and the other way around) are constrained to zero. Within each cluster, coefficients between series at the same hierarchical level are nonzero. Lastly, we enforce bottom-up propagation of dynamics: lagged changes in lower-level series may affect changes in higher-level series, but not vice versa. The sparsity patterns of $\boldsymbol{\Gamma}_1$ for the three clusters analyzed in Section 4 can be found in Table 4.

When $\Gamma_{1}$ is constrained to the block-structured sparsity pattern described above, we refer to the resulting model as the sparse VECM. This specification is expected to mitigate overfitting while enhancing the estimation stability relative to the unrestricted VECM. More importantly, the nonzero entries reflect hierarchical dependencies among all COD temporal patterns, consistent with the corresponding clustering framework.

Table 4.

Sparsity patterns of $\boldsymbol{\Gamma}_1$ under three clustering structures. The 1’s indicate free parameters in each case.

Further, existing studies indicate that period indices $\kappa_{i,t}$ may exhibit conditional heteroskedasticity (Lee and Miller Reference Lee and Miller2001; Lin et al. Reference Lin, Wang and Tsai2015; Zhou Reference Zhou2019). To accommodate this, we model the residuals $\varepsilon_{i,t}$ in (5.1) using a univariate GARCH(1,1) specification:

(5.2) \begin{align} \varepsilon_{i,t} &= \sigma_{i,t}\,z_{i,t}, \nonumber \\[3pt] \sigma_{i,t}^{2} &= \omega_{i}^{G} + \alpha_{i}^{G}\,\varepsilon_{i,t-1}^{2} + \beta_{i}^{G}\,\sigma_{i,t-1}^{2}, \end{align}

where $z_{i,t}\overset{\mathrm{i.i.d.}}{\sim}N(0,1)$ . For parsimony, the cross-cause terms are omitted from the cross-cause impact in the conditional variance equations.Footnote ⁷

To jointly estimate the parameters of the sparse VECM and the GARCH margins, we leverage MLE, which requires an explicit form for the full joint density of the innovation vector $\boldsymbol{\varepsilon}_t$ . A copula will be used to join all margins and capture the cross-cause dependence. In Section 5.3, we explain the chosen copula and derive the complete log-likelihood function.

5.3. Dependence modeling via hierarchical Archimedean copulae

The rationale for employing a copula lies in its ability to construct a multivariate distribution from specified marginal distributions (Sklar Reference Sklar1959). Let $F_i(\! \cdot \!)$ denote the marginal cumulative density function (CDF) of innovation $\varepsilon_{i,t}$ of $\kappa_{i,t}$ . Then, Sklar’s theorem implies

\begin{equation*}F(\varepsilon_{\mathrm{CA}, t},\dots,\varepsilon_{\mathrm{VA}, t})= C\bigl(F_\mathrm{CA}(\varepsilon_{\mathrm{CA}, t}),\dots,F_\mathrm{VA}(\varepsilon_{\mathrm{VA}, t});\,\boldsymbol{\theta}_C\bigr),\end{equation*}

where $C(\! \cdot \!)$ is the copula function, and $\boldsymbol{\theta}_C$ is the vector of the dependence parameters (Nelsen Reference Nelsen2007). Under the GARCH(1, 1) specification in (5.2), we form standardized residuals $z_{i,t}=\varepsilon_{i,t}/\sigma_{i,t}$ and map them to the unit interval via

\begin{equation*}u_{i,t} = \Phi(z_{i,t}),\end{equation*}

where $\Phi$ is the standard normal CDF. The transformed variables $u_{i,t}$ are approximately uniform in (0,1), making them suitable inputs for the copula.

In high-dimensional applications, hierarchical Archimedean copulae (HAC) offer a tractable and interpretable way to encode nested dependence structures (McNeil and Nešlehová Reference McNeil and Nešlehová2008; Joe Reference Joe1997). The HAC tree structure aligns well with our hierarchical clustering to model the dependence of COD innovations. In particular, we can capture stronger dependencies within the cluster and weaker links between clusters, consistent with the sparsity pattern assumed for the sparse VECM. Let $C_{\text{H}}(u_{1,t},\dots,u_{7,t};\,\boldsymbol{\theta}_C)$ denote a selected type of HAC. Its density is

\begin{equation*}c^{\text{H}}(u_{\mathrm{CA}, t},\dots,u_{\mathrm{VA}, t};\,\boldsymbol{\theta}_C)\;=\;\frac{\partial^7}{\partial u_{\mathrm{CA}, t}\cdots\partial u_{\mathrm{VA}, t}}\;C^{\text{H}}(u_{\mathrm{CA}, t},\dots,u_{\mathrm{VA}, t};\,\boldsymbol{\theta}_C).\end{equation*}

To align the copula with the hierarchical clustering in Section 4, we assign one Archimedean generator to each nesting level within a cluster and a root generator across clusters (Joe Reference Joe1997; McNeil and Nešlehová Reference McNeil and Nešlehová2008). Taking the DTW-based clustering to illustrate, let $\varphi_{1}$ and $\varphi_{2}$ be the bottom- and top-level generators for Cluster 1 (increasing trend group: MD,NE then DM), let $\varphi_{3},\varphi_{4},\varphi_{5}$ be the bottom, intermediate, and top generators for Cluster 2 (declining trend group: CA, EC then VA then Flu), respectively, and let $\varphi_{0}$ denote the root generator.

The bottom-level copula for Cluster 1 is

\begin{equation*}C^{H}_\mathrm{MD,NE}(u_{\mathrm{MD},t},u_{\mathrm{NE},t})=\varphi_{1}^{-1}[\varphi_{1}(u_{\mathrm{MD},t})+\varphi_{1}(u_{\mathrm{NE},t})],\end{equation*}

and the cluster-level copula is

\begin{equation*}C^{H}_{1}(u_{\mathrm{MD},t},u_{\mathrm{NE},t},u_{\mathrm{DM},t})=\varphi_{2}^{-1}\bigl[\varphi_{2}\bigl(C^H_\mathrm{MD,NE}\bigr)+\varphi_{2}(u_{\mathrm{DM},t})\bigr].\end{equation*}

Similarly, the bottom, intermediate, and cluster-level copulae for cluster 2 are

\begin{align*}C^H_\mathrm{CA,EC}(u_{\mathrm{CA},t},u_{\mathrm{EC},t})&= \varphi_{3}^{-1}[\varphi_{3}(u_{\mathrm{CA},t}) + \varphi_{3}(u_{\mathrm{EC},t})],\\[4pt]C^H_\mathrm{CA,EC,Vasc}(u_{\mathrm{CA},t},u_{\mathrm{EC},t},u_{\mathrm{VA},t})&= \varphi_{4}^{-1}\bigl[\varphi_{4}\bigl(C^H_\mathrm{CA,EC}\bigr) + \varphi_{4}(u_{\mathrm{VA},t})\bigr],\\[4pt]C^H_{2}(u_{\mathrm{CA},t},u_{\mathrm{EC},t},u_{\mathrm{VA},t},u_{\mathrm{Flu},t})&= \varphi_{5}^{-1}\bigl[\varphi_{5}\bigl(C^H_\mathrm{CA,EC,VA}\bigr) + \varphi_{5}(u_{\mathrm{Flu},t})\bigr],\end{align*}

respectively. Finally, the root-level copula binds the two clusters together:

\begin{equation*}C^{H}(u_{\mathrm{CA},t},\dots,u_{\mathrm{VA},t})= \varphi_0^{-1} \left\{ \varphi_0 \left[ C^{H}_{1}(u_{\mathrm{MD},t}, u_{\mathrm{NE},t}, u_{\mathrm{DM},t}) \right] + \varphi_0 \left[ C^{H}_{2}(u_{\mathrm{CA},t}, u_{\mathrm{EC},t}, u_{\mathrm{VA},t}, u_{\mathrm{Flu},t}) \right]\right\}.\end{equation*}

This HAC structure concentrates dependence within each nested subgroup before linking the clusters at the root, which matches the dependence structure imposed by the sparse VECM.

In the Archimedean copula, each generator function $\varphi\colon[0,1]\to[0,\infty]$ is a continuous (see, e.g., Okhrin and Ristig Reference Okhrin and Ristig2014), strictly decreasing and convex mapping that satisfies $\varphi(1)=0$ (and, for strict families, $\varphi(0)=\infty$ ) (Joe Reference Joe1997; Nelsen Reference Nelsen2007). The copula itself is constructed by adding the uniforms transformed by the generator and then applying the pseudo-inverse $\varphi^{-1}$ . Thus, $\varphi$ translates marginal ranks into a common scale, and its shape controls the strength and symmetry of dependence.

In our HAC construction, the root generator $\varphi_{0}$ governs the overall dependence between clusters, while nested generators ( $\varphi_{1}, \dots, \varphi_{5}$ ) capture associations within clusters at each level of the hierarchy. In this paper, we assume a common generator family for these $\varphi$ ’s, the parameters differing at each node to reflect the corresponding dependence. This parsimonious specification preserves interpretability while capturing the hierarchical dependence implied by our sparse-VECM structure. Using the HAC package (Okhrin and Ristig Reference Okhrin and Ristig2014) in (R Core Team 2025), we consider four Archimedean families: Frank, Joe, Clayton, and Gumbel. The model selection is then based on the Bayesian information criterion (BIC). Since each candidate copula contains a single dependence parameter, choosing the family with the lowest BIC is equivalent to selecting the copula with the maximized log-likelihood.

Given the selected Archimedean family, the full log-likelihood of the sparse VECM–GARCH model is

(5.3) \begin{equation}Loglik(\boldsymbol{\Theta}) = \sum_{t=1}^T\Biggl[ \ln c^{\text{H}}(u_{\mathrm{CA},t},\dots,u_{\mathrm{VA},t};\,\boldsymbol{\theta}_C) \;+\;\sum_{i \in \text{COD}}\ln\, f_i\bigl(\varepsilon_{i,t};\,\boldsymbol{\theta}_{S})\Biggr], \end{equation}

where $\boldsymbol{\theta}_{S}$ is the parameter vector of the sparse VECM-GARCH model, consisting of $\alpha_i$ , $\beta_i$ , $\boldsymbol{\Gamma}_1$ in the conditional mean equation, and the GARCH parameters $\omega^G_i$ , $\alpha^G_i$ , and $\beta^G_i$ . The marginal density $f_{i}(\! \cdot \!)$ is given by the Gaussian GARCH(1,1) specification. Jointly maximizing $Loglik(\boldsymbol{\Theta})$ with respect to $\boldsymbol{\Theta}=(\boldsymbol{\theta}_C,\boldsymbol{\theta}_S)^{\top}$ delivers the maximum likelihood estimates.

6.1. Forecasting steps in sparse VECM

Given the estimated sparse-VECM parameter vectors and matrices $\widehat{\boldsymbol{\mu}}$ , $\widehat{\boldsymbol{\alpha}}$ , $\widehat{\boldsymbol{\beta}}$ , and $\widehat{\boldsymbol{\Gamma}}_1$ , the h-step ahead forecasts of the period-indices $\kappa_t$ are obtained by iterating the conditional expectation of Equation (5.1), with $\boldsymbol{\varepsilon}_{T+h}$ set to zero. Denote by $\widehat{\kappa}_{T+h}$ and $\widehat{\Delta \kappa}_{T+h}$ the h-step-ahead forecasts of $\kappa_t$ and $\Delta \kappa_t = \kappa_t - \kappa_{t-1}$ , respectively. Then, the one-step forecast is

\begin{equation*}\Delta\widehat{\boldsymbol\kappa}_{T+1}=\widehat{\boldsymbol{\mu}}+\widehat{\boldsymbol{\alpha}}\,\widehat{\boldsymbol{\beta}}'\,\boldsymbol{\kappa}_{T}+\widehat{\boldsymbol{\Gamma}}_{1}\,\Delta\boldsymbol{\kappa}_{T},\quad\widehat{\boldsymbol\kappa}_{T+1}=\boldsymbol{\kappa}_{T}+\Delta\widehat{\boldsymbol\kappa}_{T+1}.\end{equation*}

For $j=2,\dots,h$ , we recursively compute

\begin{equation*}\Delta\widehat{\boldsymbol{\kappa}}_{T+j}= \widehat{\boldsymbol{\mu}}+ \widehat{\boldsymbol{\alpha}}\, \widehat{\boldsymbol{\beta}}'\, \widehat{\boldsymbol{\kappa}}_{T+j-1}+ \widehat{\boldsymbol{\Gamma}}_{1}\, \Delta\widehat{\boldsymbol{\kappa}}_{T+j-1},\quad\widehat{\boldsymbol{\kappa}}_{T+j}= \widehat{\boldsymbol{\kappa}}_{T+j-1}+ \Delta\widehat{\boldsymbol{\kappa}}_{T+j}.\end{equation*}

Given the predicted $\widehat{{\kappa}}_{i,T+h}$ for COD i, the corresponding logarithmic mortality rate is forecasted by

\begin{equation*}\ln \widehat{m}_{i,T+h}(x) = a_{i}(x) + b_{i}(x)\widehat{{\kappa}}_{i,T+h}.\end{equation*}

Table 5.

Cointegration tests.

Note: T.S. is the test statistics. C.V. is the critical value.

$^*$ Means the null hypothesis $H_0$ is rejected at 1% significance level.

Using exposure at time T, the total mortality rate is then forecasted byFootnote ⁹

(6.1) \begin{equation} \widehat{m}_{T+h}(x)=\dfrac{\sum_i \widehat{m}_{i,T+h}(x) E_{i,T}(x)}{\sum_i E_{i,T}(x)}.\end{equation}

6.2. Preliminary analysis

Leveraging the clustering structures from Section 4, we re-estimate the LC parameters using the sequential reconciliation procedure for our training sample. In Figure 6, we plot these reconciled period indices alongside the original LC estimates for comparison. Notably, the three reconciled $\boldsymbol{\kappa}$ ’s for nephritis display a clear upward trend, correcting the spurious decline observed in the original LC fit, as previously highlighted in Figure 3. We also observe that bottom-level indices (e.g., cancer and external causes) exhibit greater variation across clustering schemes, while higher-level indices (e.g., influenza) remain nearly unchanged. Despite these localized differences, all three reconciled clustering structures yield broadly consistent temporal patterns at the COD level.

Before fitting a sparse VECM, we need to validate the cointegration assumption. In addition to the three clustering-based reconciliations, we include the simple approach (no clustering) discussed in Section 2 of the Supplementary Materials for comparison. The re-estimated LC $\kappa$ ’s are then assessed for cointegration. Table 5 reports the results of the Johansen test (Johansen Reference Johansen1991), including the maximum eigenvalue and trace approaches. In each case, the null hypothesis of no cointegration is uniformly rejected, while the null hypothesis of at most one cointegration relationship is never rejected. These findings complement studies, including Arnold and Glushko (Reference Arnold and Glushko2022), and provide strong empirical support for our sparse VECM specification with a single error-correction term.

Next, we examine the best-performing HAC for each clustering approach. As described in Section 5.3, we compare four one-parameter Archimedean families, namely Gumbel, Clayton, Frank, and Joe. Since each family has a single dependence parameter, the BIC ranking is equivalent to the copula log-likelihood ranking. For each clustering, we re-estimate the LC period indices $\kappa_{i,t}$ , fit the sparse VECM in Section 5.2, and extract the residuals $\varepsilon_{i,t}$ . We then fit each candidate copula to the corresponding pseudo-observations and compute the resulting log-likelihoods. Table 6 shows that the Gumbel copula achieves the highest log-likelihood in all three clusters. This result is consistent with previous evidence of a predominant upper-tail dependence on mortality data (Wang et al. Reference Wang, Yang and Huang2015; Zhou Reference Zhou2019; Li et al. Reference Li, Balasooriya and Liu2021). Consequently, we employ the Gumbel-HAC for all subsequent analyses.

Figure 6.

Estimated $\kappa_{i,t}$ (normalized) using clustering information.

Table 6.

Log-likelihood values for HAC types and clustering methods.

Table 7.

Out-of-sample RMSE of total mortality forecasts (2009–2018).

Note: Bold number indicates the smallest RMSE.

The sparse VECM–GARCH–Gumbel-HAC is fitted for each clustering method. The estimated copula parameters are shown in Figure 7. Recall that for the Gumbel family, the dependence parameter $\theta\gt1$ corresponds to upper-tail dependence $\tau=1-1/\theta$ , and the nesting constraint

\begin{equation*}\theta_{\mathrm{bottom}}\;\ge\;\theta_{\mathrm{cluster}}\;\ge\;\theta_{\mathrm{root}}\end{equation*}

ensures the assumed hierarchical relationship, in which immediate siblings exhibit stronger coupling than broader groupings. When comparing the three HAC trees, the DTW and slope clusterings yield similar overall shapes. In DTW, the (CA, EC) node has a higher dependence than the (MD, NE) node, whereas slope clustering produces a slightly lower root level $\theta$ . In contrast, both the slope and trend polarity clusters assign nearly equal bottom-level $\theta$ values to the increasing and declining clusters, indicating symmetric dependence within each cluster.

Figure 7.

Fitted hierarchical structures for HAC modeling.

6.3. Forecasting results

Using the estimated sparse VECM parameters, we generate out-of-sample forecasts of total mortality rates for 2009–2018 following the procedure in Section 6.1. The forecast accuracy is measured by the root of the mean squared errors (RMSE) in all ages $x=1,\dots,85+$ and forecast horizons $h=1,\dots,10$ :

\begin{equation*}\mathrm{RMSE}=\sqrt{\frac{1}{85\times10} \sum_{x=1}^{85}\sum_{h=1}^{10} [\ln m_{T+h}(x)\;-\;\ln \widehat{m}_{T+h\mid T}(x)]^{2}}\,,\end{equation*}

where $\ln m_{T+h}(x)$ is the total mortality rate observed in 2009–2018. The corresponding forecast rate $\ln \widehat{m}_{T+h\mid T}(x)$ is obtained by (6.1).

Table 7 compares the out-of-sample RMSE for four models: the standard LC forecast using total mortality rates only, and the three clustering-specific sparse VECM–GARCH–Gumbel-HAC approaches based on COD rates. The trend-polarity clustering achieves the lowest RMSE (0.1459), followed by the slope clustering (0.1483) and the DTW clustering (0.1500). All three hierarchical models outperform the univariate LC benchmark (0.1601), with an improvement ranging from around 6% (DTW) to 10% (Trend-polarity). These results confirm that incorporating structured information from COD dynamics enhances the accuracy of aggregate mortality forecasts.

To assess the robustness of our findings, we performed four sensitivity tests, with the first three based on RMSE: shifting the training start year from 1970 to 1980; truncating the training sample at 2003 and extending the forecast horizon to 15 years (2004–2018); (iii) aggregating ages into five-year groups; and using the additional metric of mean absolute error (MAE) and mean absolute percentage error (MAPE). Table 8 reports the new out-of-sample RMSE for the four models concerned above. Our findings in the baseline scenario hold across all alternative scenarios.

Table 8.

Sensitivity of RMSE under alternative samples and age groupings and additional accuracy metrics.

In summary, the out-of-sample results reveal three key insights: (1) Models that incorporate reconciled cause-specific information consistently outperform those based solely on aggregated mortality. (2) Among the three clustering schemes, the trend-polarity grouping achieves the greatest reduction in RMSE, as well as in MAE and MAPE. This demonstrates its superior ability to capture short-run sparse dynamics in COD rates while respecting their long-run trajectories. (3) Due to its robust forecast performance, we can adopt the sparse VECM with trend-polarity clustering for in-sample analysis over the entire 1970–2018 period.Footnote ¹⁰ The in-sample analysis is presented in the next section.

7.1. Connectedness measure

To capture system-wide spillovers beyond pairwise interactions, we adopt the connectedness framework of Diebold and Yilmaz (Reference Diebold and Yilmaz2012), which quantifies dynamic spillovers among multiple interdependent variables. By transforming classic forecast-error variance decompositions (FEVD) into a unified weighted directed network, this approach enables the identification of net transmitters and receivers while ensuring invariance to variable ordering. In contrast to traditional VAR matrices or plain FEVD, which are limited to direct pairwise effects or isolated shock magnitudes, the connectedness network maps directional spillovers among all variables simultaneously. This feature helps uncover both direct and indirect transmission channels, thereby identifying systemically important pathways.

We describe how connectedness is computed in the context of COD mortality. Let $\boldsymbol{\varepsilon}_t$ be the N-vector of reduced-form residuals from the sparse VECM in Equation (5.1). Assume that

\begin{equation*}\boldsymbol{\varepsilon}_t = \boldsymbol{B}_0\,\boldsymbol{e}_t,\quad\boldsymbol{e}_t\sim MN(\mathbf{0},\mathbf{I}_N),\quad\boldsymbol{\Sigma} = \mathrm{Var}(\boldsymbol{\varepsilon}_t)=\boldsymbol{B}_0\,\boldsymbol{B}_0',\end{equation*}

where $\boldsymbol{e}_t$ are orthonormal structural shocks.

Rewriting the sparse VECM in vector-moving-average form gives

\begin{equation*}\triangle \boldsymbol{\kappa}_{t+h}=\sum_{k=0}^{h-1}\boldsymbol{\Upsilon}_k\,\boldsymbol{\varepsilon}_{t+h-k}=\sum_{k=0}^{h-1}\boldsymbol{\Upsilon}_k\,\boldsymbol{B}_0\,\boldsymbol{e}_{t+h-k},\end{equation*}

with $\boldsymbol{\Upsilon}_k$ the moving-average coefficient matrices. The orthogonal impulse response function (IRF), $\boldsymbol{e}_i'\,\boldsymbol{o}_h\,\boldsymbol{B}_0\,\boldsymbol{e}_j$ , aims to measure the effect on $\kappa_{i,t+h}$ of an one-unit shock in $e_{j,t}$ . $\boldsymbol{o}_h=\sum_{m=0}^{h}\boldsymbol{\Upsilon}_m$ is the cumulative effect of the moving-average terms. However, because $\boldsymbol{B}_0$ is typically obtained by Cholesky decomposition of $\Sigma_\varepsilon$ , this orthogonalization depends on the variable ordering.

To avoid ordering dependence, we follow Diebold and Yilmaz (Reference Diebold and Yilmaz2012) and use the generalized IRF (GIRF) and generalized FEVD (GFEVD) of Pesaran and Shin (Reference Pesaran and Shin1998). The GIRF is

\begin{equation*}\mathrm{GIRF}_\mathrm{ij}(h)=\frac{\boldsymbol{e}_i'\,\boldsymbol{o}_h\,\boldsymbol{\Sigma}\,\boldsymbol{e}_j}{\sqrt{\boldsymbol{\Sigma}_\mathrm{jj}}},\end{equation*}

and the GFEVD is

\begin{equation*}\tau_\mathrm{ij}(h)=\frac{\displaystyle\sum_{k=0}^{h-1}\bigl(\boldsymbol{e}_i'\,\boldsymbol{o}_k\,\boldsymbol{\Sigma}\,\boldsymbol{e}_j\bigr)^{2}\;/\;\boldsymbol{\Sigma}_\mathrm{jj}}{\displaystyle\sum_{k=0}^{h-1}\boldsymbol{e}_i'\,\boldsymbol{o}_k\,\boldsymbol{\Sigma}\,\boldsymbol{o}_k'\,\boldsymbol{e}_i}.\end{equation*}

Here, $\tau_\mathrm{ij}(h)$ is the fraction of the h-step forecast error variance of $\kappa_{i,t}$ due to shocks in $\kappa_{j,t}$ .

We normalize these contributions as

\begin{equation*}\widetilde{\tau}_\mathrm{ij}(h)=\frac{\tau_\mathrm{ij}(h)}{\sum_{j=1}^{N}\tau_\mathrm{ij}(h)},\end{equation*}

and define the net connectedness of cause j by

\begin{equation*}\rho_{j}=\sum_{h=1}^{H}\sum_{i\neq j}\widetilde{\tau}_\mathrm{ij}(h)-\sum_{h=1}^{H}\sum_{i\neq j}\widetilde{\tau}_\mathrm{ji}(h),\end{equation*}

where H is the limit of the forecast horizon. A positive $\rho_j$ indicates that the cause j is a net transmitter of mortality shocks, while a negative $\rho_j$ indicates that it is a net receiver.

7.2. Empirical results

As shown in Section 6, the sparse-VECM–GARCH–Gumbel-HAC model with trend-polarity clustering delivers the best out-of-sample forecasts. Its strong predictive accuracy implies minimal overfitting, making it well suited for in-sample connectedness analysis over the full 1970–2018 period. Consequently, we re-estimate the LC parameters in the full sample, extract the reconciled $\kappa_{i,t}$ , and fit the sparse-VECM–GARCH–Gumbel-HAC under trend-polarity clustering to compute generalized FEVDs and connectedness measures.

Table 9 presents the GFEVD for Cluster 1 (increasing-trend causes) and Cluster 2 (decreasing-trend causes), respectively. For each horizon $h=1,5,$ and 10, the entries show the percentage of forecast-error variance in each series explained by cause-specific shocks.

Table 9.

GFEVD for Impulse of Cluster 1: Diabetes, Mental Disorders, Nephritis and for Impulse of Cluster 2: Cancer, External, Influenza, Vascular.

Own shocks dominate short-horizon uncertainty across all causes, typically accounting for 70–90% of the variance at $h=1$ . Over longer horizons, the contributions of the own shock decline, while cross-cause spillovers increase. For example, by $h=10$ , diabetes shocks explain more than 18% of nephritis variance, whereas vascular’s own contribution falls below 50%.

Within the increasing-trend cluster, mental disorders and nephritis interact more weakly. Mental shocks account for only about 1–2% of the variance in nephritis, whereas diabetes shocks have a much greater impact as described above. In the decreasing-trend cluster, most causes remain highly self-contained. Cancer, external causes, and influenza each retain their own shock shares close to 90%. Nonetheless, cross-cluster spillovers are present. Nephritis shocks drive over 20% of influenza variance and roughly 10% of vascular variance at $h=10$ . nephritis itself absorbs about 17% of the variance from vascular shocks by that horizon.

Table 10 presents the connectedness measures up to step 10. The table shows that, for each cause i (rows) and shock j (columns), the percentage of cause-i’s forecast-error variance explained by the shocks to cause j, summed over $h=1,\dots,10$ . These percentages are normalized; therefore, each row sums to 100%. The diagonal entries give the own-shock share, and the off-diagonals make the “From Others” contributions (listed in the final column). In contrast, the “To Others” row at the bottom reports, for each cause j, the column sum of off-diagonal entries. This represents the total spillover of cause-j’s shocks transmitted to all other causes. The bold “Net” row shows the difference between “To Others” and “From Others,” indicating whether a cause is a net transmitter (positive) or net receiver (negative) of mortality dynamics.

Table 10.

Connectedness summary for the seven COD categories (reordered).

The “To Others” measures reveal a clear separation among the seven causes for transmitting influence. Diabetes, mental disorders, nephritis, influenza, and vascular diseases each transmit over 30% of their shocks to other causes, indicating a strong systemic influence. In contrast, cancer and external causes transmit only 9.9% and 12.2%, respectively, reflecting their relative isolation.

On the receiving side, we see a similar separation. Mental disorders, nephritis, influenza, and vascular diseases each absorb more than 35% of their variance from other causes. By comparison, diabetes, cancer, and external causes are more self-contained, with less than 25% of their variance explained by shocks from other CODs.

Aggregating the off-diagonal FEVD shares and averaging over the six counterparts yields a total connectedness index of 36%. This signals a moderate level of interdependence among COD trends. This moderate spillover suggests that mortality reduction efforts may benefit from coordinated strategies across multiple causes rather than isolated interventions on only certain COD.

More importantly, the seven causes can be classified into three groups, based on their net connectedness:

\begin{equation*}\begin{aligned}\text{Net transmitters:}&\;\text{Diabetes }({+}10.0\%),\;\text{nephritis }({+}5.4\%),\;\text{vascular disease }({+}6.7\%),\\[3pt]\text{Net receivers:} &\;\text{Cancer }({-}12.3\%),\;\text{external causes }({-}12.8\%),\\[3pt]\text{Near-neutral:} &\;\text{Mental disorders }({+}1.4\%),\;\text{influenza }({+}1.5\%).\end{aligned}\end{equation*}

These findings align well with clinical, epidemiological, and statistical understanding. Chronic diseases such as diabetes and nephritis frequently precipitate cardiovascular complications, generating substantial interrelated mortality risks. A collaborative meta-analysis of 102 prospective studies found that diabetes confers a twofold excess risk of vascular disease, independently of other conventional risk factors (Emerging Risk Factors Collaboration 2010). Similarly, patients with chronic kidney disease face markedly elevated cardiovascular mortality. The corresponding mortality risks increase progressively with declining renal function, reaching two- to threefold higher rates in stages 3–5 of chronic kidney disease (Go et al. Reference Go, Chertow, Fan, McCulloch and Hsu2004). In contrast, population-based cancer registry data show that the vast majority of deaths among cancer patients are attributable to their index malignancy rather than other causes. Improved SEER estimates indicate that net cancer-specific survival accounts for over 80% of the observed deaths in many common cancers, with noncancer mortality comprising only a small fraction (Howlader et al. Reference Howlader, Ries, Mariotto, Reichman, Ruhl and Cronin2010). Consequently, cancer mortality remains largely self-contained, with limited reciprocal impact on other causes. External causes of mortality due to accidents are unlikely to interact with medical causes. In a competing-risk framework, death from one cause precludes death from more isolated causes, such as cancer and external causes. Consequently, these CODs emerge as net receivers of shocks when all causes are modeled jointly.