Identifiable Regime Detection in Pension Fund Networks via Sticky Hidden Markov Models

1. Introduction

Systemic risk, the potential for shocks to propagate across financial systems, has been prominent since the 2008 crisis, revealing how the failure of individual institutions can destabilise economies [1]. While traditionally associated with banks because of their central role in credit intermediation, attention is now shifting to large non-bank institutional investors, especially pension funds, which hold over $63 trillion in assets globally [2]. The risk arises not simply from scale; research demonstrates that overlapping portfolios among institutions allow asset sales or losses to depress prices and transmit shocks, amplifying stress throughout the system [3,4,5]. Notably, pension funds form networks characterised by such overlaps, making indirect contagion possible even without direct financial links.

Unlike banks, pension funds typically have longer investment horizons, limited leverage, and focus on hedging [6]. Yet these safeguards do not dissolve the vulnerability that [5] identifies, because that vulnerability lives in the overlap structure, not in leverage or direct exposure. Fire-sale dynamics and portfolio overlap can still open indirect contagion channels with no direct financial links between funds [3,4]. Recent events, like the 2026 collapse of the Versorgungswerk der Zahnärztekammer Berlin pension fund, underscore that systemic vulnerabilities in this sector are real [7]. The question, then, is not whether pension networks can transmit systemic stress, but how one can know, from the data, when such a network has entered a high-risk state.

Indeed, existing literature increasingly frames the risks confronting pension schemes as fundamentally unhedgeable, system-level risks that cannot be diversified away through portfolio construction alone [8]. The enterprise risk management tradition reinforces this view [9,10]: [10] shows formally that a firm sponsoring a defined-benefit pension plan incurs compounding vulnerabilities unless pension risk is integrated with all other business risks in a unified optimisation. Their model demonstrates that managing pension risk in a silo, rather than holistically, reduces expected firm value by over $12\%$. This is not merely an actuarial accounting problem; it reveals that portfolio homogeneity and shared liability structures across pension participants create interdependencies that travel well beyond the balance sheet of any single fund.

Before turning to how systemic risk has been modelled, it is worth naming the obstacle that (we will argue) every existing approach must confront, and that ultimately defines the gap this paper fills. Systemic episodes in pension funds are rare and slow-moving. [11] addresses this through modular modelling but notes that traditional balance-sheet-based indicators often reflect past conditions and cannot reliably predict the onset of stress. [12] surveying the broader machine learning for financial risk management literature arrives at a complementary diagnosis from the opposite direction: across market, credit, operational, and insurance risk [13], the dominant paradigm is supervised classification trained on labelled crises or externally specified targets, and the field’s principal open challenges, data inaccessibility, label scarcity, imbalanced datasets, and biased portfolio-level estimates, all trace back to the same root that [11] identifies. Where [11] argues from systems-design first principles that the labels cannot be there, [12] documents, study by study, that the field has been quietly assuming they are. Four main modelling traditions respond to this constraint, but each leaves gaps, particularly for pension fund networks.

Financial network frameworks often model systemic risk using explicit matrices of bilateral obligations [14,15], which are well-suited to banking but less applicable to pension funds, whose connections are indirect. A complementary strand models systemic risk through network density and contagion probability directly, showing that the relationship between interconnectedness and expected loss is non-monotonic; neither fully connected nor ring networks are optimal because denser connections simultaneously enable risk sharing and multiply contagion pathways [16]. Formalising when and how such contagion becomes systemic, however, requires a rigorous measurement framework. The mathematical foundations for measuring such systemic risk rigorously were laid by [17], who developed a general framework for systemic risk measures via multi-dimensional acceptance sets and aggregation functions, unifying both the first aggregate, then inject capital, and the first inject capital, then aggregate paradigms into a single axiomatic structure. This framework makes explicit that the choice of aggregation rule and acceptance set jointly determines where the contagion boundary lies, a theoretical insight that later computational work has sought to operationalise. This foundation has recently been extended in two influential directions. [18] generalise systemic risk measures to graph-structured data, modelling the system as a random asset vector together with a random interbank liability matrix, and proving the existence of optimal bailout-capital allocations that secure the network before contagion occurs. In a closely related vein, [19] operationalises the first allocate, then aggregate paradigm from [17], using a deep learning scheme to compute the minimal cash needed to secure a system when no closed-form solution exists. These models are mainly normative, asking how to allocate capital given a detected crisis, rather than diagnosing when a high-risk regime emerges. The supervised machine learning tradition emerged partly to fill this diagnostic gap but, in doing so, introduced a different and equally binding constraint: the need for labelled crises. Recent advances integrating network analysis with extreme dependency modelling and machine learning [20,21,22,23,24,25,26,27,28,29] improve early warning capabilities and interpretability but remain reliant on labelled historical crises.

Deep neural networks and hybrid architectures achieve high predictive accuracy but deepen rather than dissolve this label dependency [23,30,31]. The relationship between risk factors and crisis outcomes is nonlinear [30], and more sophisticated architectures confirm rather than escape this; a hybrid XGBoost-LSTM model can issue warnings 2.3 quarters in advance with 83.7% accuracy [31], yet every such figure is conditional on prior crisis labelling. This label dependency is not merely a technical inconvenience but a structural constraint on the entire supervised paradigm: [32] demonstrates that even well-specified ensemble classifiers such as random forests, which outperform panel logit models on global banking crisis data with AUROC exceeding 81%, require crisis labels drawn from externally defined databases such as [33], and their generalisability across heterogeneous economies remains contingent on the completeness and accuracy of those labels. Interpretability advances do not dissolve it either; [34] applies partial dependence plots across multiple models to identify exchange rates, money supply, and interest rates as dominant risk drivers, yet still requires a continuous quantitative risk target derived from market data, which is unavailable in the pension fund setting where NAV-based returns are the primary observable. These requirements are feasible in banking and broad equity markets, but much less so for pension fund networks, where systemic episodes are rare, slow-moving, and not cleanly labelled, and where the very question of interest is whether elevated-risk regimes can be recognised without prior labels.

Unsupervised data-driven methods present a promising alternative. [35] introduces the absorption ratio, the fraction of total asset-return variance explained by a finite number of eigenvectors from a rolling PCA, as an unsupervised, threshold-free indicator of market fragility. The logic is direct: when a single dominant factor absorbs the bulk of return variation across assets, those assets are tightly coupled, and a shock propagates more quickly and broadly than when risks are dispersed. Applied to U.S. equity industries over 1998-2010, the absorption ratio rose to its highest recorded level during the global financial crisis of 2008, and spikes in the standardised shift in the absorption ratio preceded all of the 1% worst monthly drawdowns. Crucially, [35] derive their signal entirely from observed returns: no crisis labels, no externally chosen thresholds, and no imposed regime classification. The same PCA-based logic has since been applied to broader contagion analysis [1] and forms the methodological core of our own absorption ratio component. Studies also apply clustering and anomaly detection to textual or market data [36,37], but typically do not operate directly on institutional financial positions. Their network is a network of what is said about firms, not of what firms hold, which is exactly the overlap-based contagion substrate that [5] shows to be the operative one for institutions with common asset holdings. The broader systemic-risk machine learning literature confirms how isolated this position is. [38] in their survey of machine learning for systemic risk, classify the field into four branches and identify data-driven systemic risk analysis as a major future direction precisely because real interbank-network data are difficult to obtain and most existing work relies either on proxy networks or on externally specified targets. Endogenous regime detection operating directly on actual financial positions does not appear in their taxonomy as an established branch; it appears as a gap.

The argument across the preceding traditions can now be stated crisply, and it defines the specific opening this paper fills. [15] shows formally that network topology governs how shocks clear but requires observed bilateral liabilities. [39] shows empirically that PageRank and interbank-exposure ratios are the dominant drivers of systemic impact and vulnerability, respectively, but requires observed bilateral exposures and synthetic DebtRank labels. [30] shows that nonlinear classifiers substantially outperform linear ones for crisis detection but amplify rather than resolve the label dependency. Each step increases explanatory power while deepening the data requirements. The pension case sits at the intersection of all three constraints: no bilateral liabilities, no DebtRank-computable exposures, and no labelled crises. The natural response is to step back to an earlier, harder question: can the state of the network be identified directly from returns, without any of these inputs?

Within unsupervised methods applied directly to financial returns, the toolkit is well established: PCA-based connectedness measures [1,35]; DTW-based time-series clustering, robust to temporal misalignment and well suited to financial series [40]; and Markov regime-switching, whose foundations [41] motivate the HMM approach we adopt. The CRISP-DM framework [42] provides a reproducible scaffold for assembling these components. Notably, the most relevant pension-specific precedents identify market regimes via HMM and stress-testing frameworks, anchored to external benchmarks [43,44]. Together, these establish both the network embeddedness of pension funds and the usefulness of regime-switching methods in this exact setting, yet remain conditional on external states rather than allowing risk regimes to emerge endogenously from the pension network itself. In [11]’s terms, a benchmark-anchored regime is a measurement of a previous, externally defined state rather than a forward-looking read of the network’s own condition.

We address this gap by developing an unsupervised regime-detection framework that treats systemic risk as an emergent property of the pension fund network, rather than as a response to externally defined market conditions. The proposed approach integrates three components within a unified empirical pipeline: a PCA-based absorption ratio following the logic of [35] to measure system-wide co-movement; DTW-based hierarchical clustering to capture cross-sectional similarity structure; and a Gaussian HMM to identify latent risk regimes over time. Importantly, and in the spirit of the threshold-free, data-derived discipline shared by [35], [36], and [37], the crisis threshold is derived endogenously from the HMM emission distributions, ensuring internal consistency between regime classification and crisis identification.

Three methodological contributions distinguish this framework. First, the number of latent states is selected by combining BIC with a posterior identifiability analysis in a finite sticky HMM [45], addressing BIC’s tendency to over-select when mixture components overlap. Second, the stickiness prior is calibrated empirically from the data through dual anchors corresponding to the baseline and the BIC-optimal specification, eliminating the circularity of choosing a prior that determines the result. Third, emission estimates from the frequentist Baum-Welch algorithm are cross-validated against a fully Bayesian sticky HMM estimated via the No-U-Turn Sampler, with both approaches producing near-identical results, confirming that the detected regime structure is data-driven rather than model-dependent. The framework is applied to daily NAV data from $N=40$ second-pillar pension funds in Lithuania managed by five providers over January 2019 to September 2025, augmented with fund-specific benchmarks and global market indices.

The empirical analysis yields four main findings. First, systemic co-movement is persistently high: the absorption ratio computed using the method of [35] never falls below 0.728, indicating that a single common factor explains at least 73% of total return variance even in the calmest periods. Clustering results show that the funds group primarily by age cohort rather than by provider, suggesting that lifecycle allocation regulation, rather than managerial differences, is the dominant driver of cross-fund dependence, a structural reading consistent with the finding that pension risk behaviour is regulation-driven [10,46]. Second, three latent risk regimes are identified, with moderate- and high-risk states accounting for nearly 87% of trading days, indicating that elevated systemic interaction is the network’s structural norm. Third, systemic stress is not a passive reflection of global markets: the relationship between the absorption ratio and global benchmarks is statistically insignificant across regimes, while tracking-error amplification, the precise quantity practitioners monitor [46], is concentrated in the high-risk state. Fourth, the empirically derived crisis threshold $\tau =0.851$ exceeds the values commonly used in network-based contagion studies, where thresholds are typically chosen heuristically [47]. This is consistent with [5]’s demonstration that contagion ignites only beyond a system-specific tipping point, one that depends on leverage, concentration, and the commonality of holdings rather than on any universal shock level. More broadly, the finding that an internally derived threshold sits above the heuristic values routinely used in banking-network studies echoes the methodological conclusion that [15] reaches from the opposite direction: the clearing structure of a specific network, its liabilities, its cash-flow distribution, and the topology of its obligations determine where the contagion boundary lies, and that boundary cannot be read off from first principles or borrowed from other sectors.

The remainder of the paper is structured as follows. Section 2 develops the theoretical framework and describes the empirical methodology, including the data, the absorption ratio, DTW clustering, and the Hidden Markov Model. Section 3 presents the main results, including the regime detection, the cluster structure, and the cross-step analysis connecting them. Section 4 discusses the findings and their policy implications. Section 5 concludes and outlines directions for future research.

2. Theoretical Framework and Methodology

This section develops the theoretical framework and explains how it is implemented empirically. The analysis follows the Cross-Industry Standard Process for Data Mining (CRISP-DM) [42], which structures the workflow into data understanding, preparation, modelling, and evaluation. The three-step empirical pipeline, PCA absorption ratio, DTW-based hierarchical clustering, and HMM regime detection are summarised in Figure 1. Each step is presented together with the theoretical results that justify it.

2.1. Data

The dataset comprises daily NAV observations from pension funds operating in Lithuania. The sample covers $N=40$ funds managed by five pension fund providers, anonymised as Providers 1 through 5 to comply with data confidentiality requirements. Each provider manages eight age-cohort funds, yielding a structured $5\times 8$ cross-sectional panel. The observation period runs from January 2019 to September 2025, yielding approximately $T=1,696$ daily trading observations per fund and a total of 67,794 fund-day observations. The panel is strongly balanced, with each of the 40 funds contributing between 1,692 and 1,696 daily observations over the sample period.

The age-cohort design reflects differences in lifecycle allocation across funds. Seven cohorts are defined by participant birth decade: cohorts AG1 through AG7 cover birth years from 1954-1960 to 1996-2002, capturing participants at progressively earlier career stages. The eighth cohort, AG8, corresponds to the Target Income Pension Fund (TIPF), a specialist vehicle for participants who have reached retirement age and are actively drawing income. This distinction is empirically consequential: preliminary descriptive analysis of the sample reveals that TIPF funds exhibit substantially lower daily return volatility than age-cohort funds, consistent with their capital preservation mandate (see Table 1). This distinction anticipates the finding in Section 3.2 that extreme-cohort funds (AG1 and AG8) separate sharply from middle cohorts across all five providers, forming dedicated clusters with tracking errors three to six times lower than those of middle-cohort clusters (Table 5). The presence of TIPF funds is therefore analytically significant, as they provide a conservative baseline against which the behaviour of age-cohort funds during stress periods can be assessed.

For each fund and trading day, the dataset contains the fund unit price, an associated comparative benchmark index value, a unique fund identifier, and the corresponding age-group classification. The dataset is further augmented with $M=3$ publicly available global benchmark series: the MSCI World Index, the MSCI Europe Index, and the S&P 500 Index. The sample period includes the COVID-19 market shock, providing a natural benchmark for validating the regime detection framework against a well-documented systemic stress event.

2.1.1. Notation and Setup

Let $(\mathsf{\Omega },\mathcal{F},\mathbb{P})$ be a probability space. Consider N pension funds observed over T trading days, and let $X_{i,t}$ denote the NAV of fund i at time t, for $i=1,\dots ,N$ and $t=0,\dots ,T$. Each fund i is associated with a comparative benchmark index $Z_{i,t}^{\left(f\right)}$, reflecting its underlying investment strategy. In addition, let $Z_{m,t}^{\left(g\right)}$, $m=1,\dots ,M$, denote a set of global benchmark series representing key market factors. The corresponding log returns are defined by

$$\begin{array}{cccc}\hfill r_{i,t}^{\left(f\right)}& =log\left(X_{i,t}\right)-log\left(X_{i,t-1}\right),\hfill & \hfill i& =1,\dots ,N,t=1,\dots ,T,\hfill \\ \hfill r_{i,t}^{(b,f)}& =log\left(Z_{i,t}^{\left(f\right)}\right)-log\left(Z_{i,t-1}^{\left(f\right)}\right),\hfill & \hfill i& =1,\dots ,N,\hfill \\ \hfill r_{m,t}^{(b,g)}& =log\left(Z_{m,t}^{\left(g\right)}\right)-log\left(Z_{m,t-1}^{\left(g\right)}\right),\hfill & \hfill m& =1,\dots ,M.\hfill \end{array}$$

(1)

Let $R^{\left(f\right)}\in {\mathbb{R}}^{T\times N}$ denote the pension fund return matrix, $R^{(b,f)}\in {\mathbb{R}}^{T\times N}$ the fund-specific benchmark return matrix, and $R^{(b,g)}\in {\mathbb{R}}^{T\times M}$ the global benchmark return matrix. When required, we define the augmented return matrix

$$R=\left[R^{\left(f\right)}{R}^{(b,f)}{R}^{(b,g)}\right]\in {\mathbb{R}}^{T\times (2N+M)}.$$

For a rolling window of length $\omega$, define the windowed return matrix $R_t$ for $t=\omega ,\dots ,T$, containing returns from periods $t-\omega +1$ to t.

2.1.2. Descriptive Statistics

Table 1 reports descriptive statistics for the daily log return series of pension funds, disaggregated by age cohort, and global benchmarks over the full sample period. All series exhibit negative skewness and substantial excess kurtosis, indicating fat-tailed distributions that deviate markedly from normality. This is consistent with the presence of extreme return realisations during stress periods like the COVID-19 pandemic.

A clear lifecycle gradient is visible across age cohorts. The oldest active cohort (AG1, born 1954-1960) and the TIPF funds exhibit the lowest volatility, consistent with their conservative, capital-preservation mandates. Volatility increases monotonically with cohort youth, reaching 0.0083 for the youngest cohort (AG7, born 1996-2002), reflecting progressively higher equity allocations prescribed by lifecycle regulation. Skewness is most negative for the oldest cohorts ($-1.7$ for AG1, $-1.5$ for TIPF), indicating greater exposure to left-tail events relative to the magnitude of their returns. Excess kurtosis is highest for AG1 and AG2, suggesting that conservative funds experience more extreme return realisations relative to their low baseline volatility, a pattern consistent with infrequent but sharp revaluations of bond-heavy portfolios during interest rate shocks.

Importantly, the excess kurtosis values of 10.7–18.0 reported in Table 1 do not invalidate the Gaussian emission assumption of the Hidden Markov Model introduced in Section 2.4: each regime may be approximately Gaussian, while the mixture of regimes across the full sample produces the fat tails observed in the unconditional distribution. This is precisely the distributional flexibility that motivates a regime-switching framework.

Among the global benchmarks, all three indices exhibit similar volatility, negative skewness, and substantial excess kurtosis exceeding 16, confirming that the fat-tailed behaviour observed in pension fund returns is mirrored in global equity markets over the sample period.

Table 1. Descriptive statistics of daily log returns. Pension fund series are reported by age cohort (pooled across providers) and for the full sample. Skewness and excess kurtosis are reported. All series exhibit negative skewness and substantial excess kurtosis, consistent with fat-tailed return distributions.

Series	Mean	Std	Min	Max	Skewness	Kurtosis
Pension funds by age cohort
AG1 (54/60)	0.0001	0.0022	−0.0301	0.0159	−1.703	17.977
AG2 (61/67)	0.0002	0.0043	−0.0551	0.0371	−1.526	17.447
AG3 (68/74)	0.0004	0.0069	−0.0867	0.0661	−1.208	13.941
AG4 (75/81)	0.0004	0.0080	−0.0868	0.0662	−1.050	10.825
AG5 (82/88)	0.0004	0.0081	−0.0864	0.0659	−1.042	10.720
AG6 (89/95)	0.0004	0.0081	−0.0866	0.0663	−1.039	10.730
AG7 (96/02)	0.0004	0.0083	−0.0902	0.0695	−1.004	10.668
AG8 (TIPF)	0.0001	0.0019	−0.0235	0.0132	−1.495	14.549
Full sample	0.0003	0.0065	−0.0902	0.0695	−1.217	16.882
Global benchmarks
MSCI World	0.0004	0.0127	−0.1239	0.0871	−1.169	16.188
MSCI Europe	0.0006	0.0121	−0.1208	0.0871	−0.942	17.861
S&P 500	0.0006	0.0127	−0.1277	0.0909	−0.664	16.012

Table 1. Descriptive statistics of daily log returns. Pension fund series are reported by age cohort (pooled across providers) and for the full sample. Skewness and excess kurtosis are reported. All series exhibit negative skewness and substantial excess kurtosis, consistent with fat-tailed return distributions.

Series	Mean	Std	Min	Max	Skewness	Kurtosis
Pension funds by age cohort
AG1 (54/60)	0.0001	0.0022	−0.0301	0.0159	−1.703	17.977
AG2 (61/67)	0.0002	0.0043	−0.0551	0.0371	−1.526	17.447
AG3 (68/74)	0.0004	0.0069	−0.0867	0.0661	−1.208	13.941
AG4 (75/81)	0.0004	0.0080	−0.0868	0.0662	−1.050	10.825
AG5 (82/88)	0.0004	0.0081	−0.0864	0.0659	−1.042	10.720
AG6 (89/95)	0.0004	0.0081	−0.0866	0.0663	−1.039	10.730
AG7 (96/02)	0.0004	0.0083	−0.0902	0.0695	−1.004	10.668
AG8 (TIPF)	0.0001	0.0019	−0.0235	0.0132	−1.495	14.549
Full sample	0.0003	0.0065	−0.0902	0.0695	−1.217	16.882
Global benchmarks
MSCI World	0.0004	0.0127	−0.1239	0.0871	−1.169	16.188
MSCI Europe	0.0006	0.0121	−0.1208	0.0871	−0.942	17.861
S&P 500	0.0006	0.0127	−0.1277	0.0909	−0.664	16.012

2.2. Step 1: The PCA Absorption Ratio as a Spectral Systemic Risk Measure

2.2.1. Definition and Computation

Definition 1 

(Rolling Covariance Matrix). The rolling sample covariance matrix at time t is

$${\mathsf{\Sigma }}_t={\displaystyle \frac{1}{\omega -1}}{R}_t^{\top }M{R}_t,$$

where $M=I_{\omega }-{\displaystyle \frac{1}{\omega }}{\mathbf{1}}_{\omega }{\mathbf{1}}_{\omega }^{\top }$ is the centering matrix. Since ${\mathsf{\Sigma }}_t$ is real, symmetric, and positive semi-definite, it admits an eigendecomposition. Let ${\lambda }_1\left(t\right)\ge {\lambda }_2\left(t\right)\ge \dots \ge {\lambda }_{N+M}\left(t\right)\ge 0$ denote the ordered eigenvalues. The first principal component, the linear combination of fund returns with maximum explained variance, corresponds to ${\lambda }_1\left(t\right)$.

In financial economics, PCA applied to the return covariance matrix is widely used for dimensionality reduction and for identifying a small number of dominant common factors that capture the co-movement of asset returns [48,49,50].

Definition 2 

(PCA Absorption Ratio). The absorption ratio at time t is the fraction of total fund return variance explained by the first principal component:

$$A{R}_t={\displaystyle \frac{{\lambda }_1\left(t\right)}{{\sum }_{j=1}^N{\lambda }_j\left(t\right)}}={\displaystyle \frac{{\lambda }_1\left(t\right)}{\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)}}.$$

$A{R}_t$ is a direct PCA output: it measures how much of the total variation across all N fund returns is driven by a single common factor. High values indicate that funds are moving together, a hallmark of systemic stress.

2.2.2. Empirical Implementation

The absorption ratio is computed from the pension fund return matrix $R^{\left(f\right)}$ using a rolling window of $\omega =60$ trading days, balancing estimation stability with the ability to capture regime changes at approximately monthly frequency. At each time point t, the following steps are performed:

1.

Extract the windowed return matrix $R_t^{\left(f\right)}$;

2.

Compute the unbiased sample covariance matrix ${\mathsf{\Sigma }}_t$;

3.

Perform eigendecomposition to obtain ordered eigenvalues ${\lambda }_1\left(t\right)\ge \dots \ge {\lambda }_N\left(t\right)\ge 0$;

4.

Compute $A{R}_t={\lambda }_1\left(t\right)/\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)$.

The resulting series ${\left\{A{R}_t\right\}}_{t=\omega }^T$ gives the share of total return variance explained by the first principal component of ${\mathsf{\Sigma }}_t$. A rising $A{R}_t$ signals increasing commonality in fund returns and thus greater systemic co-movement (consistent with the strictly monotone relationship established in Proposition 3).

2.2.3. Key Properties

Proposition 1 

(Boundedness). For all t, $A{R}_t\in \left[\frac{1}{N},1\right]$.

Proof. 

By eigenvalue ordering, ${\lambda }_1\left(t\right)\le \mathrm{tr}\left({\mathsf{\Sigma }}_t\right)$, giving $A{R}_t\le 1$. Since ${\lambda }_1\left(t\right)\ge {\lambda }_j\left(t\right)$ for all j, summing over j gives $N·{\lambda }_1\left(t\right)\ge \mathrm{tr}\left({\mathsf{\Sigma }}_t\right)$, so $A{R}_t\ge \frac{1}{N}$.    □

Proposition 2 

(Interpretation of Extremes). $A{R}_t=\frac{1}{N}$ if and only if all eigenvalues are equal. $A{R}_t=1$ if and only if $\mathrm{rank}\left({\mathsf{\Sigma }}_t\right)=1$.

Proof. 

$A{R}_t=\frac{1}{N}$ implies ${\lambda }_1=\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)/N$, the mean eigenvalue; equal eigenvalues are the unique PSD solution. $A{R}_t=1$ implies ${\lambda }_1=\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)$, so ${\lambda }_j=0$ for $j\ge 2$, giving $\mathrm{rank}\left({\mathsf{\Sigma }}_t\right)=1$.    □

Proposition 3 

(Monotonicity in Cross-Fund Correlation). Let ${\overline{\rho }}_t$ denote the average pairwise return correlation. Under the equicorrelation model ${\mathsf{\Sigma }}_t={\sigma }^2[(1-{\overline{\rho }}_t)I_N+{\overline{\rho }}_t{\mathbf{1}}_N{\mathbf{1}}_N^{\top }]$, the absorption ratio satisfies

$$A{R}_t={\displaystyle \frac{1+(N-1){\overline{\rho }}_t}{N}},$$

which is strictly increasing in ${\overline{\rho }}_t$ for ${\overline{\rho }}_t\in \left(-\frac{1}{N-1},1\right]$.

Proof. 

Under equicorrelation, ${\mathsf{\Sigma }}_t$ has one eigenvalue ${\lambda }_1={\sigma }^2(1+(N-1){\overline{\rho }}_t)$ and $N-1$ repeated eigenvalues ${\lambda }_j={\sigma }^2(1-{\overline{\rho }}_t)$. Therefore, $\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)=N{\sigma }^2$ and $A{R}_t=(1+(N-1){\overline{\rho }}_t)/N$. The derivative $\partial A{R}_t/\partial {\overline{\rho }}_t=(N-1)/N>0$ establishes strict monotonicity.    □

Proposition 3 is the key justification for using $A{R}_t$ as a systemic risk measure: it is a strictly increasing function of average cross-fund correlation, so a rising absorption ratio directly and monotonically reflects rising systemic co-movement.

Theorem 1 

(Crisis Identification). Let $C_t=\mathbb{I}\{A{R}_t>\tau \}$, $\tau \in \left(\frac{1}{N},1\right)$, and set ${\overline{\rho }}_{\mathrm{crit}}=(N\tau -1)/(N-1)$. Under equicorrelation,

$$C_t=1\iff {\overline{\rho }}_t>{\overline{\rho }}_{\mathrm{crit}}.$$

The threshold $\tau$ provides an economically interpretable mapping between the absorption ratio and the underlying cross-fund correlation structure: a crisis state is identified whenever average inter-fund correlation ${\overline{\rho }}_t$ exceeds ${\overline{\rho }}_{\mathrm{crit}}$. Rather than imposing $\tau$ a priori, its empirical value is determined from the HMM emission parameters in Step 3 (Section 2.4), ensuring consistency between the regime detection framework and the crisis identification criterion.

2.2.4. Relationship to Standard Risk Decomposition

The absorption ratio admits a direct interpretation within the standard risk decomposition framework used in portfolio performance measurement [50,51]. In factor model terms, the total risk of the pension fund network at time t decomposes as

$$\underset{\mathrm{Total}\mathrm{Risk}}{\underbrace{\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)}}=\underset{\mathrm{Systematic}\mathrm{Risk}}{\underbrace{{\lambda }_1\left(t\right)}}+\underset{\mathrm{Specific}\mathrm{Risk}}{\underbrace{\sum _{j=2}^N{\lambda }_j\left(t\right)}},$$

where systematic risk $={\lambda }_1\left(t\right)$ is the variance explained by the first principal component (the dominant common factor), and specific risk $=\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)-{\lambda }_1\left(t\right)$ is the residual variance attributable to fund-specific (idiosyncratic) behaviour. The absorption ratio, therefore, equals the systematic risk share:

$$A{R}_t={\displaystyle \frac{\mathrm{Systematic}\mathrm{Risk}}{\mathrm{Total}\mathrm{Risk}}}={\displaystyle \frac{{\lambda }_1\left(t\right)}{\mathrm{tr}\left({\mathsf{\Sigma }}_t\right)}},1-A{R}_t={\displaystyle \frac{\mathrm{Specific}\mathrm{Risk}}{\mathrm{Total}\mathrm{Risk}}}.$$

This equivalence provides a bridge between the spectral systemic risk measure used in this paper and the standard risk metrics employed in portfolio performance evaluation. However, the economic interpretation of a high systematic risk share differs sharply between the two settings. In standard portfolio construction, where PCA-based factor models are used to optimise investor portfolios [49,50], a high systematic risk share is generally interpreted as indicating that idiosyncratic exposures have been diversified away, leaving the investor primarily exposed to common market-wide factors. In the pension fund network context, the same quantity carries the opposite interpretation. A high systematic risk share indicates that lifecycle allocation rules have compressed the cross-section of fund behaviour to the point where provider choice offers negligible risk differentiation. When $A{R}_t=0.87$, 87% of all fund return variance is driven by a single common factor, leaving only 13% as fund-specific risk—a level of structural homogeneity that represents a systemic vulnerability rather than a diversification achievement.

2.2.5. Robustness to Benchmark Inclusion and Covariance Estimator

In the baseline analysis, $A{R}_t$ is constructed from $R^{\left(f\right)}$ alone, preserving the fund-level focus of the primary systemic risk measure. The benchmark series serves two subsequent roles. First, periods of elevated $A{R}_t$ are examined alongside $R^{(b,f)}$ and $R^{(b,g)}$ to assess whether heightened systemic co-movement coincides with deviations from fund mandates or with broader market stress. Second, as a robustness exercise, $A{R}_t$ is recomputed using the augmented return matrix R, thereby assessing the sensitivity of regime detection to the inclusion of benchmark market factors.1

2.2.5.1. Choice of Covariance Estimator

The unbiased sample covariance is selected as the baseline because it preserves the full empirical concentration of variance in the leading eigenvalue ${\lambda }_1\left(t\right)$, which is the precise quantity $A{R}_t$ is designed to measure. Three alternative families of estimators are commonly considered in this context: shrinkage estimators such as Ledoit-Wolf [52], robust estimators such as the minimum covariance determinant (MCD), and Random Matrix Theory (RMT) filtered estimators based on the Marchenko-Pastur (MP) eigenvalue distribution, which separate signal eigenvalues from finite-sample noise [53,54]. Each is appropriate for a different inferential goal: shrinkage stabilises the inverse covariance matrix for portfolio optimisation, MCD downweights outliers to robustly estimate mean-variance, and RMT filtering removes spurious correlations from short panels. None of these objectives align with the absorption ratio’s purpose, which is to detect rather than smooth concentration in the leading eigenvalue.

To verify that this choice does not drive the empirical results, $A{R}_t$ is recomputed using all three alternative estimators as a robustness exercise (Appendix A). The MP-filtered series tracks the baseline almost exactly (correlation $0.98$, mean absolute deviation $0.011$), confirming that ${\lambda }_1\left(t\right)$ lies well outside the MP noise bulk throughout the sample and therefore reflects genuine cross-fund co-movement rather than estimation noise. The MCD estimator yields similarly close agreement on average. The Ledoit-Wolf estimator produces a systematic downward level shift, consistent with its design objective of compressing the eigenvalue spectrum toward a structured target.

2.3. Step 2: DTW Hierarchical Clustering

2.3.1. Formal Definition

Let $\mathcal{X}=\{x=(x_1,\dots ,x_T):x_t\in \mathbb{R}\}$ be the space of univariate real-valued time series of length T. For $x,y\in \mathcal{X}$, define the local cost matrix $C\in {\mathbb{R}}^{T\times T}$ with entries $C_{ij}={(x_i-y_j)}^2$.

Definition 3 

(Warping Path). A warping path $W=\{(i_1,j_1),\dots ,(i_L,j_L)\}$ satisfies: (i) boundary: $(i_1,j_1)=(1,1)$ and $(i_L,j_L)=(T,T)$; (ii) monotonicity: $i_{k+1}\ge i_k$ and $j_{k+1}\ge j_k$; (iii) step size: $(i_{k+1}-i_k,j_{k+1}-j_k)\in \{(1,0),(0,1),(1,1)\}$.

Definition 4 

(DTW Distance [55]). The DTW distance between $x,y\in \mathcal{X}$ is $\mathrm{DTW}(x,y)={min}_W\sqrt{{\sum }_{(i,j)\in W}{C}_{ij}}$, computed via the accumulated cost matrix

$$D_{ij}=C_{ij}+min\{D_{i-1,j},D_{i,j-1},D_{i-1,j-1}\}$$

with $D_{11}=C_{11}$, and $\mathrm{DTW}(x,y)=\sqrt{D_{TT}}$.

2.3.2. Semi-Metric Properties and Consequences for Clustering

Theorem 2. 

DTW dissimilarity satisfies non-negativity, identity $\left(\mathrm{DTW}\right(x,x)=0)$, and symmetry $\left(\mathrm{DTW}\right(x,y)=\mathrm{DTW}(y,x\left)\right)$, but does not satisfy the triangle inequality in general and therefore does not define a metric.

Proof. 

Non-negativity follows from $C_{ij}\ge 0$. Identity holds since, for $x=y$, the diagonal warping path yields zero cost. Symmetry follows from the symmetry of the local cost function, as ${(x_i-y_j)}^2={(y_j-x_i)}^2$. The failure of the triangle inequality is a known property of DTW; see [56].    □

Since DTW operates in a non-Euclidean setting and permits nonlinear temporal alignments, standard Euclidean k-means centroids are generally not well-suited to temporally misaligned sequences. DTW barycenter averaging (DBA) [57] addresses this by constructing centroids consistent with DTW alignments, thereby extending centroid-based clustering to DTW-based dissimilarities. However, the resulting optimisation problem remains nonconvex and does not inherit the theoretical structure of classical Euclidean k-means. Within hierarchical clustering, single-linkage methods may exhibit instability under non-metric dissimilarities and do not, in general, preserve ultrametric structure when the underlying dissimilarity fails to satisfy metric axioms.

We therefore adopt Ward’s hierarchical clustering as the primary method, since it groups observations by minimising within-cluster dispersion and extends naturally to general dissimilarity settings [58]. Ward’s linkage offers two additional advantages relevant to the present application: first, it produces a fully nested hierarchical decomposition, enabling the subdivision analysis used to validate the optimal cluster number (Section 3.2); second, its variance-minimising criterion aligns naturally with the provider-cohort block structure of the pension fund network.

2.3.3. Empirical Implementation

DTW-based hierarchical clustering is applied to the N pension fund return series to identify groups of funds with similar return dynamics over the full sample period. Each return series is first standardised to zero mean and unit variance, so that clustering reflects dynamic patterns rather than differences in scale across cohorts. The $N\times N$ pairwise DTW distance matrix is then computed using the dtaidistance Python library without imposing a warping window constraint.

The optimal number of clusters $K_c^{*}$ is selected by minimising the Davies-Bouldin (DB) index over $K_c\in \{2,\dots ,8\}$, with two diagnostic checks. First, the search range is extended to $K_c=15$ to confirm that finer partitions reflect over-segmentation rather than new structure (Section 3.2). Second, solutions producing singleton clusters are excluded as degenerate. The DB index is preferred over the silhouette coefficient because it does not require a metric space, making it appropriate for DTW-based dissimilarities. Robustness against k-means alternatives, including k-means with DTW barycenters, is verified in Section 3.2. This clustering step provides a cross-sectional characterisation of fund behaviour that complements the time-series analysis of systemic risk in Steps 1 and 3.

2.4. Step 3: Hidden Markov Model Regime Detection

2.4.1. Model Specification

Let ${\left\{O_t\right\}}_{t=1}^T$ denote the observed absorption ratio sequence, where $O_t=A{R}_t\in [\frac{1}{N},1]$. We model $\left\{O_t\right\}$ as generated by an unobserved discrete-time Markov chain ${\left\{S_t\right\}}_{t=1}^T$ taking values in the finite state space $\mathcal{S}=\{1,\dots ,K\}$.

Definition 5 

(Hidden Markov Model). A K-state Hidden Markov Model (HMM) is characterised by the parameter set $\lambda =(\pi ,A,B)$, where:

$$\begin{array}{cccc}\hfill {\pi }_k& =\mathbb{P}(S_1=k),\hfill & \hfill & \mathit{initial}\mathit{distribution},\hfill \\ \hfill A_{ij}& =\mathbb{P}(S_t=j\mid S_{t-1}=i),\hfill & \hfill & \mathit{transition}\mathit{probabilities},\hfill \\ \hfill B_k\left(o\right)& =\mathcal{N}(o;{\mu }_k,{\sigma }_k^2),\hfill & \hfill & \mathit{Gaussian}\mathit{emission}\mathit{density}.\hfill \end{array}$$

Conditional on the latent state $S_t=k$, observations therefore follow a Gaussian distribution, $A{R}_t\mid S_t=k\sim \mathcal{N}({\mu }_k,{\sigma }_k^2)$.

2.4.2. Estimation, Decoding, and Long-Run Behaviour

Parameter Estimation

Model parameters are estimated via the Baum-Welch (EM) algorithm using the forward-backward recursions. The resulting updates maximise the likelihood $\mathbb{P}(O\mid \lambda )$ and converge to a stationary point of the likelihood function under standard regularity conditions [59]. Gaussian HMMs are estimated for $K\in \{2,\dots ,8\}$, with multiple random initialisations to mitigate convergence to local optima; the model achieving the highest log-likelihood is retained for each K.

State Decoding

The most likely latent state sequence $S^{*}=(s_1^{*},\dots ,s_T^{*})$ is obtained using the Viterbi algorithm,

$$S^{*}=arg\underset{S}{max}\mathbb{P}(S\mid O,\lambda ),$$

with computational complexity $\mathcal{O}\left(K^{2}T\right)$. The states are subsequently ordered by their estimated emission means ${\widehat{\mu }}_k$ and labelled as Low-, Moderate-, and High-Risk regimes.

Forecasting and Long-Run Behaviour

Conditional on the current state $S_t$, the distribution of $S_{t+h}$ for any horizon $h\ge 1$ is obtained by iterating the transition kernel:

$$\mathbb{P}(S_{t+h}=k\mid S_t=i)={\left(A^h\right)}_{ik},$$

where $A^h$ denotes the $h^{\mathrm{th}}$ matrix power of A. Letting ${\mathbf{e}}_i$ be the $i^{\mathrm{th}}$ unit vector, the forward distribution over states is given by ${\mathbf{e}}_i^{\top }{A}^h$. If the chain is irreducible and aperiodic, $A^h$ converges as $h\to \infty$ to a rank-one matrix whose rows equal the stationary distribution ${\pi }^{*}$, defined as the unique probability vector satisfying

$${\pi }^{*\top }A={\pi }^{*\top },{\pi }^{*\top }\mathbf{1}=1.$$

Equivalently, ${\pi }^{*}$ is the left eigenvector of A associated with the eigenvalue 1, and represents the long-run unconditional probability of occupying each regime, independently of the starting state. Forward state probabilities are reported at horizons of one day, one week, one month, one quarter, and six months.

2.4.3. Model Selection

BIC-Based Screening

Candidate models are first evaluated using the Bayesian Information Criterion (BIC),

$$\mathrm{BIC}\left(K\right)=-2log\mathbb{P}(O\mid {\widehat{\lambda }}^{\left(K\right)})+d\left(K\right)logT,$$

where $d\left(K\right)=K^2+2K$ is the number of free parameters: $K^2$ transition probabilities (including the initial state distribution) and $2K$ Gaussian emission parameters. While BIC provides a standard likelihood-based selection criterion, it tends to favour more complex models that may include weakly identified states [60,61]. A complementary identifiability analysis is therefore conducted to discriminate genuine regime structure from over-parameterisation.

Calibration of Regime Persistence

The Bayesian specification employs a sticky transition prior to encourage persistence of latent regimes. For each state k, the transition probabilities follow a Dirichlet prior,

$${\pi }_k\sim \mathrm{Dirichlet}(\alpha \mathbf{1}+\kappa {\mathbf{e}}_k),$$

corresponding to the finite-state analogue of the sticky HDP-HMM prior of [45], where $\kappa$ governs the strength of self-transition persistence. Under this prior, the expected self-transition probability is

$$\mathbb{E}\left[{\pi }_{kk}\right]={\displaystyle \frac{\alpha +\kappa }{K\alpha +\kappa }},$$

(2)

which, with $\alpha =1$, inverts to $\kappa =(K\overline{p}-1)/(1-\overline{p})$, where $\overline{p}$ denotes the target average self-transition probability.

Rather than fixing $\kappa$ arbitrarily, it is calibrated empirically from the EM-fitted transition matrix. To verify that the calibration is not itself contingent on a particular choice of K, two anchors are used: the operational baseline ($K=3$) and the BIC-optimal specification ($K=5$) identified in the empirical analysis. This dual calibration yields two values of $\kappa$, both of which are carried forward into the identifiability analysis.

Bayesian Identifiability Analysis

A Bayesian sticky HMM is estimated for $K\in \{3,4,5,6\}$ using the No-U-Turn Sampler (NUTS), with four parallel chains per configuration. The procedure is repeated under both calibrated values of $\kappa$ to ensure that the identifiability conclusions do not depend on the EM specification used to anchor the prior. Three diagnostics are used: (1) the Gelman-Rubin statistic ($\widehat{R}$), reported both raw and after label-switching correction obtained by sorting posterior emission means within each draw; (2) the minimum separation between posterior mean emission parameters; and (3) the number of divergent transitions. A state is considered identified when its label-switching-corrected $\widehat{R}<1.1$.

Sensitivity Analysis

To verify that the identifiability results are not artefacts of the calibrated $\kappa$, the analysis is repeated across an additional grid of $\kappa$ values spanning weak to strong persistence priors. The three-state specification is expected to remain identified across the full range of empirically supported $\kappa$ values, whereas configurations that depend on a narrow prior window are flagged as less robust.

Selected Specification

The number of states is determined by combining BIC with the Bayesian identifiability and sensitivity results. While BIC favours a higher-dimensional specification, the identifiability analysis reveals that additional states beyond $K=3$ either fail to separate cleanly or represent finer subdivisions of the same underlying regime structure rather than genuinely distinct economic states. The three-state specification is therefore retained as the smallest parsimonious model that is both statistically identifiable and economically interpretable across both $\kappa$ calibrations and the full sensitivity grid. Higher-state specifications are examined as robustness checks in the empirical results. All computations are implemented in Python using the hmmlearn library.

2.4.4. Crisis Identification

In addition to regime classification, a data-driven crisis threshold $\tau$ is defined as the absorption ratio level at which the moderate- and high-Risk emission densities have equal mass. Under Gaussian emissions with approximately equal variance, this coincides with the midpoint of the two regime means. Crisis periods are then identified as $C_t=\mathbb{I}\{A{R}_t>\tau \}$ (Theorem 1), providing a complementary threshold-based measure of extreme observations alongside the Viterbi-decoded regime assignment. To assess economic relevance, the inferred regimes are further compared with the global benchmark return dynamics, allowing periods of elevated absorption ratio to be interpreted in the context of broader market conditions.

Figure 1. Empirical pipeline following the CRISP-DM framework. Step 1 produces the absorption ratio time series $\left\{A{R}_t\right\}$ from the pension fund return matrix $R^{\left(f\right)}\in {\mathbb{R}}^{T\times N}$; Step 2 characterises cross-sectional fund structure via DTW-based hierarchical clustering with Ward linkage; Step 3 infers latent systemic risk regimes via a Hidden Markov Model fitted to $\left\{A{R}_t\right\}$, with the number of states selected by combining BIC with a Bayesian identifiability analysis and regimes decoded via the Viterbi algorithm. The dashed arrow indicates the robustness exercise in which the augmented return matrix $R=\left[R^{\left(f\right)}{R}^{(b,f)}{R}^{(b,g)}\right]\in {\mathbb{R}}^{T\times (2N+M)}$ replaces $R^{\left(f\right)}$ in Step 1.

Figure 1. Empirical pipeline following the CRISP-DM framework. Step 1 produces the absorption ratio time series $\left\{A{R}_t\right\}$ from the pension fund return matrix $R^{\left(f\right)}\in {\mathbb{R}}^{T\times N}$; Step 2 characterises cross-sectional fund structure via DTW-based hierarchical clustering with Ward linkage; Step 3 infers latent systemic risk regimes via a Hidden Markov Model fitted to $\left\{A{R}_t\right\}$, with the number of states selected by combining BIC with a Bayesian identifiability analysis and regimes decoded via the Viterbi algorithm. The dashed arrow indicates the robustness exercise in which the augmented return matrix $R=\left[R^{\left(f\right)}{R}^{(b,f)}{R}^{(b,g)}\right]\in {\mathbb{R}}^{T\times (2N+M)}$ replaces $R^{\left(f\right)}$ in Step 1.

Taken together, the three steps form a single empirical pipeline. Step 1 produces the absorption ratio series that summarises systemic co-movement over time; Step 2 partitions the cross-section of funds into structurally homogeneous groups; and Step 3 maps the absorption ratio dynamics onto a small number of latent risk regimes whose boundary $\tau$ is derived from the data rather than imposed. The interpretive role of the benchmark series, the robustness of the covariance estimator, and the identifiability of the regime structure are all built into the design rather than treated as post-hoc checks. The next section applies this pipeline to the Lithuanian second-pillar pension fund panel and reports what it reveals about the structure of systemic risk in the network.

3. Main Results

The empirical analysis proceeds in three stages following the pipeline described in Section 2 and illustrated in Figure 1. Section 3.1 presents the absorption ratio dynamics and robustness results. Section 3.2 reports the DTW clustering structure. Section 3.3 presents the HMM regime detection results, including the emission distributions, the data-driven crisis threshold, and its economic interpretation. Section 3.4 connects the clustering and regime detection findings through a cross-step analysis of tracking error amplification under systemic stress.

3.1. Absorption Ratio Dynamics

Table 2 reports descriptive statistics of the baseline absorption ratio series $\left\{A{R}_t\right\}$ computed from the pension fund return matrix $R^{\left(f\right)}$ over the full sample period.

The empirical range $[0.728,0.938]$ lies strictly within the theoretical bounds established in Proposition 3, confirming the internal consistency of the measure. The absorption ratio never approaches the theoretical lower bound of $1/N=0.025$, which would obtain only if all 40 fund return series were orthogonal. Even in the calmest periods of the sample, a single common factor explains more than 72% of total return variance across the pension fund network, a level of structural co-movement fundamentally inconsistent with genuine portfolio diversification.

Figure 2 provides direct visual evidence of this synchronisation across all three return matrices. Panels 1 and 2 reveal that all five providers move in unison during the COVID-19 shock and the 2022 rate-hike cycle, with the provider-level tracks collapsing toward each other during stress episodes. The AG8/TIPF funds exhibit a visibly compressed return range relative to the active cohorts, consistent with their capital preservation mandate. Panel 3 shows that the three global equity benchmarks co-move closely with each other but follow a different temporal pattern from the pension fund series, particularly during 2022 - 2023, where global benchmarks show sustained negative returns while pension fund dynamics evolve independently.

Figure 3 places $A{R}_t$ alongside these benchmark series to examine the temporal relationship between systemic co-movement and market conditions. The absorption ratio peaks sharply during the COVID-19 shock (early 2020, $A{R}_t=0.938$) and exhibits sustained elevation during the 2022 interest rate tightening cycle. However, the correspondence between $A{R}_t$ and the global benchmarks is not straightforward: during the 2022 - 2023 period, benchmark returns in Panel 3 stabilise and recover while $A{R}_t$ in Panel 1 remains elevated, suggesting that systemic co-movement in the pension fund network is not simply a reflection of global equity market stress. A formal statistical analysis of this relationship, including significance testing, is deferred to Section 3.4, where it is integrated with the HMM regime classification to enable regime-conditional interpretation.

3.1.1. Robustness: Augmented Absorption Ratio

As a robustness exercise, $A{R}_t$ is recomputed using the augmented return matrix $R=\left[R^{\left(f\right)}{R}^{(b,f)}{R}^{(b,g)}\right]\in {\mathbb{R}}^{T\times (2N+M)}$, where $2N+M=83$. The augmented series ranges from 0.468 to 0.919, with a mean of 0.712, compared with 0.837 for the baseline. The correlation between the two series is 0.341, indicating that the inclusion of the benchmark series substantially changes the co-movement structure captured by the absorption ratio. The augmented version of the absorption ratio moves somewhat like the original version, but not too closely, because adding benchmarks changes how assets move together. This supports the methodological choice to construct the baseline covariance matrix ${\mathsf{\Sigma }}_t$ exclusively from pension fund returns, ensuring the risk measure remains fund-focused. Figure 4 compares the baseline and augmented absorption ratios over time.

The panel shows both series, the baseline $A{R}_t$ consistently remains higher, while the augmented $A{R}_t$ fluctuates more and exhibits deeper troughs. This divergence reflects the structural impact of adding benchmark-linked returns, which redistribute variance across a broader set of factors and reduce the dominance of the leading eigenvalue. To assess the materiality of regularisation, ridge regularisation is deliberately applied to the full-rank baseline covariance matrix. Although regularisation is required for the rank-deficient augmented specification, applying it to the baseline induces distortions of 0.013 – 0.047 in $A{R}_t$ across representative sample windows (Table 3). This confirms that regularisation must be restricted to the augmented robustness exercise and omitted from the baseline computation.

3.2. DTW Hierarchical Clustering

The clustering analysis proceeds in three stages: selection of the optimal cluster number, characterisation of the resulting structure, and robustness verification against extended search ranges and alternative algorithms.

3.2.1. Optimal Cluster Number

Table 4 reports the Davies-Bouldin (DB)index across $K_c\in \{2,\dots ,8\}$ clusters for Ward’s linkage applied to the pairwise DTW distance matrix. The index is minimised at $K_c^{*}=8$.

Figure 5 presents the corresponding dendrogram with the $K_c^{*}=8$ cut.

3.2.2. Cluster Structure

Table 5 reports the composition and mean tracking error of each cluster, and Figure 6 maps cluster assignments onto the provider-cohort grid.

Table 5. Cluster composition and mean tracking error at $K_c^{*}=8$. Tracking error (TE) is the mean absolute daily deviation of fund returns from their fund-specific benchmark $R^{(b,f)}$.

Cluster	Composition	Size	Mean TE
C1	Provider 3, AG2-AG7	6	0.0046
C2	Provider 1, AG2-AG7	6	0.0066
C3	Provider 1, AG1 & AG8	2	0.0014
C4	Providers 2 & 5, AG1 & AG8	4	0.0010
C5	Provider 4, AG1 & AG8	2	0.0012
C6	Provider 3, AG1 & AG8	2	0.0011
C7	Provider 4, AG2-AG7	6	0.0041
C8	Providers 2 & 5, AG2-AG7	12	0.0036

Table 5. Cluster composition and mean tracking error at $K_c^{*}=8$. Tracking error (TE) is the mean absolute daily deviation of fund returns from their fund-specific benchmark $R^{(b,f)}$.

Cluster	Composition	Size	Mean TE
C1	Provider 3, AG2-AG7	6	0.0046
C2	Provider 1, AG2-AG7	6	0.0066
C3	Provider 1, AG1 & AG8	2	0.0014
C4	Providers 2 & 5, AG1 & AG8	4	0.0010
C5	Provider 4, AG1 & AG8	2	0.0012
C6	Provider 3, AG1 & AG8	2	0.0011
C7	Provider 4, AG2-AG7	6	0.0041
C8	Providers 2 & 5, AG2-AG7	12	0.0036

Two structural patterns emerge. First, extreme-cohort funds (AG1, born 1954-1960; and AG8, the Target Income Pension Fund) separate sharply from middle cohorts (AG2-AG7) across all five providers. Clusters C3-C6 group exclusively AG1 and AG8 funds within or across providers, with mean tracking errors of 0.0010-0.0014, three to six times lower than the 0.0036-0.0066 observed for middle-cohort clusters C1, C2, C7, and C8. This is consistent with the conservative, low-turnover mandates governing capital-preservation funds: AG1 funds serve participants nearing retirement, while AG8/TIPF funds serve participants actively drawing retirement income. The compressed return range of TIPF funds is also visible in Figure 2.

Second, clustering is organised primarily by age cohort rather than by provider. Cluster C8 groups the middle cohorts of Providers 2 and 5 across provider boundaries (12 funds), while Cluster C4 groups their extreme cohorts (4 funds), indicating that these two providers operate nearly identical lifecycle allocation strategies. The heatmap in Figure 6 makes this visible as a column-dominant pattern: colour transitions are sharper across cohort columns than across provider rows. This cross-provider grouping by lifecycle stage indicates that regulatory lifecycle allocation constraints – which prescribe similar asset allocations for funds targeting the same birth-year cohort regardless of provider – are the dominant driver of return co-movement in this network, rather than provider-specific investment decisions.

3.2.3. Robustness: Extended Search Range

To verify that $K_c^{*}=8$ is not an artifact of the search range, the analysis was repeated for $K_c\in \{2,\dots ,15\}$, exceeding the eight age cohorts that define the most granular natural grouping in the data. Davies-Bouldin (DB) scores continue to fall mechanically beyond $K_c=8$, reaching 0.553 at $K_c=15$. However, this trajectory reflects subdivision of existing clusters rather than discovery of new structure, as shown in Table 6.

Two observations follow. First, the $K_c=15$ partition is fully nested within the $K_c=8$ partition: every $K_c=15$ cluster is a strict subset of one $K_c=8$ cluster, with no reassignment across $K_c=8$ boundaries. Five of the $K_c=15$ clusters contain a single fund, confirming that the lower DB scores at higher K are driven by isolating individual funds as singleton clusters rather than identifying genuine structural groups. The $K_c^{*}=8$ solution is therefore retained as the smallest K minimising DB without producing degenerate clusters.

Second, the subdivision of Cluster C8 yields a substantive observation: when forced to split the 12-fund Providers 2 & 5 middle-cohort group, the algorithm groups the youngest middle cohorts of both providers together (AG2-AG3, TE = 0.0026) while separating the older middle cohorts (AG4-AG7) by provider. This indicates that Providers 2 and 5 follow nearly identical allocation rules for the youngest cohorts but diverge slightly for cohorts approaching mid-career, consistent with greater allocation flexibility at the equity-heavy end of the lifecycle glide path.

3.2.4. Robustness: Algorithm Choice

Two k-means alternatives were also examined to verify that the cohort-dominant structure is not an artifact of the algorithm choice: k-means with Euclidean distance and k-means with DTW barycenters via the tslearn library [62]. Table 7 reports DB scores for all three methods.

Across the seven values of $K_c$ tested, Ward+DTW produces the lowest DB score at four ($K_c\in \{2,3,4,6\}$), is marginally outperformed by k-means with Euclidean distance at $K_c=5$ and by k-means with DTW barycenters at $K_c=8$, and produces identical scores at $K_c=6,7,\&8$ for the first two methods. The differences across methods are small in absolute terms (within 0.04 DB units at the optimum), and the DB-optimal solution across all three methods occurs at $K_c=8$.

Unanimous Partition at $K_c=7$

A particularly informative pattern emerges at $K_c=7$, where all three methods produce identical DB scores (0.8688). Computing the adjusted rand index (ARI) across all three pairwise method comparisons at $K_c=7$ yields $\mathrm{ARI}=1$ in every case (Ward vs. k-means Euclidean: 1, Ward vs. k-means DTW: 1, k-means Euclidean vs. k-means DTW: 1), confirming that the three methods recover identical partitions at this $K_c$. This convergence is methodologically significant: it demonstrates that the cohort-dominant block structure of the network is detected unanimously across distance metrics (Euclidean and DTW) and across algorithm families (hierarchical and partitional).

The structural relationship between the unanimous $K_c=7$ partition and the $K_c^{*}=8$ baseline is reported in Table 8. The collapse from $K_c=8$ to $K_c=7$ merges Clusters C5 (Provider 4, AG1 + AG8) and C6 (Provider 3, AG1 + AG8) into a single 4-fund group, while leaving the other six clusters unchanged. This indicates that Providers 3 and 4 follow nearly identical extreme-cohort allocation strategies, second only to Providers 2 and 5, who already share Cluster C4 at $K_c^{*}=8$. The implications of this provider-pairing pattern are discussed in Section 4.

Comparison at $K_c^{*}=8$

At the DB-optimal $K_c^{*}=8$, all three methods yield comparable scores (Ward+DTW: 0.7505, k-means Euclidean: 0.7505, k-means DTW: 0.7183), and the cluster assignments produced by k-means with DTW barycenters agree with the Ward+DTW partition at $\mathrm{ARI}=0.97$. The full cluster composition under k-means with DTW barycenters is reported in Table 9, and the cross-tabulation with the Ward+DTW partition is reported in Table 10.

Seven of the eight Ward clusters map perfectly onto a single k-mean cluster. The sole disagreement involves Ward Cluster C3 (Provider 1, AG1 + AG8), which k-means partitions by reassigning Provider 1’s AG1 fund to the Provider 1 middle-cohort group (KM2) rather than grouping it with AG8 (KM6). This is the only fund in the network whose cluster assignment is methodologically ambiguous, suggesting that Provider 1’s AG1 fund operates with a slightly more equity-tilted allocation than the AG1 funds of the other four providers, an interpretation consistent with Provider 1’s higher middle-cohort tracking error (C2: TE = 0.0066, the highest in Table 5).

Choice of Ward as Baseline

Ward’s linkage is retained as the baseline for two reasons: it is the only method tested that produces a fully nested hierarchical decomposition, enabling the dendrogram in Figure 5 and the $K_c=15$ subdivision analysis in Table 6; and it accommodates the provider-cohort block structure of the network through its variance-minimising criterion. The empirical agreement at $K_c=7$ ($\mathrm{ARI}=1$) and at $K_c=8$ ($\mathrm{ARI}=0.97$) confirms that this methodological choice does not affect the substantive findings of the paper.

3.3. HMM Regime Detection

3.3.1. Model Selection

Selecting the number of latent regimes K in a Gaussian HMM requires balancing in-sample fit against the risk of over-parameterisation. We combine the standard BIC with a Bayesian posterior identifiability analysis under two data-driven $\kappa$ calibrations, leading to a defensible choice of $K=3$ as the baseline.

BIC-Based Comparison

Table 11 reports BIC values for HMM models with $K\in \{2,\dots ,8\}$ states estimated by the Baum-Welch (BW) algorithm with 20 random initialisations each.

BIC is minimised at $K=5$. However, BIC penalises model complexity through the number of free parameters but does not assess whether the additional components correspond to interpretable or well-separated states. This limitation is well documented in the mixture-model literature [60,61]: while BIC consistently estimates the number of mixture components asymptotically, in practice it can select extra components when distributions overlap or when a single cluster is represented by multiple mixture components. We therefore complement BIC with a posterior identifiability analysis.

Dual Empirical Calibration of $\kappa$

The identifiability analysis relies on a Bayesian finite sticky HMM with the Dirichlet prior specified in Section 2.4. To avoid the circularity of choosing a prior that determines the answer, $\kappa$ is calibrated empirically by inverting the prior expected self-transition probability against the EM-fitted transition matrix. Because this calibration depends on K, we use two anchors: the operational baseline ($K=3$) and the BIC-optimal specification ($K=5$). The mean diagonals are ${\overline{p}}_3=0.971$ for the three-state EM model (diagonals $0.948,0.980,0.985$) and ${\overline{p}}_5=0.953$ for the five-state model (diagonals $0.953,0.957,0.946,0.976,0.934$), yielding calibrated values ${\kappa }_3\approx 66$ and ${\kappa }_5\approx 80$. Both values are carried forward into the identifiability analysis, ensuring that conclusions do not hinge on the choice of EM anchor.

Identifiability Under Both Calibrations

We fit the finite sticky HMM via NUTS sampling at each calibrated $\kappa$ for $K\in \{3,4,5,6\}$, using four parallel chains with 1,000 warmup and 1,000 post-warmup iterations per chain. The Gelman–Rubin $\widehat{R}$ statistic for emission means is reported both raw and after label-switching correction (sorting emission means within each posterior draw). A state is considered identified when its corrected $\widehat{R}<1.1$. Table 12 reports the results for both calibrations.

The two calibrations agree on every qualitative conclusion. $K=3$ is cleanly identified under both priors (corrected $\widehat{R}=1.000$, minimum gap $0.044$ absorption-ratio units). $K=4$ fails identification entirely under both priors, with corrected $\widehat{R}$ reaching 13.13 (${\kappa }_3$) and 14.77 (${\kappa }_5$). $K=5$ is also identified under both priors, but represents a finer subdivision of the three-regime structure with smaller emission-mean gaps ($0.019$); this case is examined as a robustness check in Section 3.3.6. $K=6$ shows partial identification under both calibrations (2 of 6 states at ${\kappa }_3$, 1 of 6 at ${\kappa }_5$), indicating further over-parameterisation. The agreement between the two calibrations confirms that the identifiability conclusions are properties of the data, not of the prior anchor.

Sensitivity to the Stickiness Prior

To verify that these conclusions extend beyond the calibrated values, the analysis is repeated for $\kappa \in \{10,20,30,50\}$ in addition to the two anchors. Table 13 reports the number of identified states across the full grid.

$K=3$ is identified at every $\kappa$ value consistent with the empirical persistence of the data ($\kappa \ge 20$). The single failure at $\kappa =10$ corresponds to a prior expected persistence (Equation 2), which ranges from 0.69 at $K=6$ to 0.85 at $K=3$, substantially below the empirical mean of 0.97. At this $\kappa$, the prior is therefore inconsistent with the data and the chains collapse into degenerate posterior regions.

$K=4$ exhibits the opposite pattern. The chains appear to converge only at $\kappa =10$, where the prior is itself inconsistent with the empirical persistence and therefore cannot be used to support $K=4$ as identified. At every $\kappa$ consistent with the data ($\kappa \ge 20$), $K=4$ fails identification under both calibrations, confirming that four genuinely distinct regimes do not match this network.

We retain $K=3$ as the baseline. It is the smallest specification cleanly identified under both data-calibrated priors and across the full empirically supported $\kappa$ range, and it aligns with the standard low-, moderate-, and high-risk regime taxonomy in the systemic risk literature [41]. The EM and Bayesian estimates at $K=3$ agree to within 0.001 across all emission parameters (Section 3.3.2), and the BIC-optimal $K=5$ specification is examined as a finer refinement in Section 3.3.6.

3.3.2. Regime Characterisation

Table 14 reports the emission parameters, state occupancy, and persistence statistics for each of the three regimes under the $K=3$ baseline, together with mean daily log returns of global equity benchmarks over the regime-assigned days.

To verify that the regime estimates are not artefacts of the EM optimisation, the emission means are re-estimated under the Bayesian sticky HMM at the data-calibrated prior ${\kappa }_3=66$. Table 15 compares the two estimators after label-switching correction.

The EM and Bayesian point estimates agree to within 0.001 across all three regimes, and the 90% credible intervals are narrow (widths of 0.002-0.007). This convergence between two methodologically distinct estimators, one frequentist and one fully Bayesian with a data-calibrated prior, confirms that the data sharply identify the three-regime structure and not an artefact of either estimation approach. Figure 7 plots the three estimated Gaussian emission distributions in relation to the crisis threshold $\tau$ derived in Section 3.3.3.

Figure 8 overlays the regime classification onto the absorption ratio and benchmark return series across the full sample period. The high-risk regime (red shading) captures the COVID-19 shock and the 2022 rate-hike period, as well as additional episodes of elevated co-movement in 2021 and 2024-2025. The decoupling between global benchmark recovery and the persistence of high-risk shading, particularly visible in late 2022, where the bottom panel shows stabilising global benchmark returns while the top panel remains in the high-risk regime, confirms that pension fund network stress is not a mirror of equity market dynamics but reflects a structural synchronisation that persists beyond the initial shock.

3.3.3. Crisis Threshold

The empirical crisis threshold $\tau$ is defined as the absorption ratio level at which the moderate- and high-risk emission densities have equal mass. Under Gaussian emissions with approximately equal variance, this coincides with the midpoint of the two regime means:

$$\tau ={\displaystyle \frac{{\widehat{\mu }}_{\mathrm{Moderate}}+{\widehat{\mu }}_{\mathrm{High}}}{2}}={\displaystyle \frac{0.8291+0.8737}{2}}=0.851.$$

By Theorem 1 and Proposition 3, this corresponds to a critical average inter-fund correlation of

$${\overline{\rho }}_{\mathrm{crit}}={\displaystyle \frac{N\tau -1}{N-1}}={\displaystyle \frac{40\times 0.851-1}{39}}\approx 0.848.$$

The secondary threshold ${\tau }_1=0.801$ separates low- from moderate-risk, corresponding to ${\overline{\rho }}_{\mathrm{crit},1}\approx 0.796$.

The threshold $\tau$ and the Viterbi decoder give two related but distinct classifications. The daily crisis indicator $C_t=\mathbb{I}\{A{R}_t>\tau \}$ flags 494 days as crisis episodes (31.3% of the sample), whereas the Viterbi-decoded high-risk regime spans 516 days (32.7%). The 22-day difference reflects the overlap zone shown in Figure 7: within this zone, Viterbi assigns the persistent state encoded in the transition matrix rather than the closer emission density on any given day. The downstream econometric analysis in Section 3.4 and Section 3.4.1 uses the Viterbi-decoded high-risk regime, while $\tau$ provides the economic interpretation of where that regime boundary lies in observation space.

The empirical value $\tau =0.851$ substantially exceeds the 0.50-0.60 range commonly assumed in the systemic risk literature. This reflects the structural homogeneity of the pension fund network and confirms that crisis identification criteria must be calibrated to the institutional context rather than imported from cross-sectoral applications.

3.3.4. Regime Persistence and Transition Structure

Beyond regime characterisation, the full transition probability matrix reveals an asymmetric escalation pathway between regimes. Table 16 reports the estimated matrix.

Direct transitions from low-risk to high-risk are very rare (probability 0.0046), meaning the system almost always passes through the Moderate-Risk regime when escalating. This carries early-warning potential: a sustained rise in $A{R}_t$ into the moderate-risk zone is a precursor signal before a full high-risk episode.

3.3.5. Current State and Forward Probabilities

As of the final observation in the sample (29 September 2025), the pension fund network is in the moderate-risk regime, with $A{R}_t=0.8285$, closely aligned with the Moderate-Risk emission mean $\widehat{\mu }=0.8291$. The system has remained in this regime for 38 consecutive trading days, against an expected total duration of 51 trading days.

Forward state probabilities $\mathbb{P}(S_{t+h}=k\mid S_t=\mathrm{Moderate}-\mathrm{Risk})$ are computed by iterating the transition matrix as defined in Section 2.4. Table 17 reports these probabilities at horizons ranging from one day to six months, alongside the long-run stationary distribution.

At the one-quarter horizon, the probability of transitioning to the high-risk regime rises to 23.1%, exceeding the probability of returning to the low-risk regime (15.5%). This asymmetry reflects the dominant role of the moderate-to-high transition channel established in Table 16: from the moderate-risk regime, escalation is structurally more likely than de-escalation at medium horizons. Over longer horizons, the forward probabilities converge towards the stationary distribution: 30.7% high-risk, 55.3% moderate-risk, and 14.0% low-risk. This long-run equilibrium confirms that the pension fund network spends approximately 86% of its time in moderate- or high-risk regimes, reinforcing the characterisation of elevated systemic co-movement as the structural norm rather than an exceptional state.

3.3.6. Robustness: Higher-State Specifications

The BIC-optimal $K=5$ specification was shown to be statistically identified under both $\kappa$ calibrations (Table 12), with all five emission means distinguishable across chains. The five estimated emission means (0.767, 0.815, 0.836, 0.857, 0.891) cluster into three groups when ordered: a single low state ($\approx 0.77$), three intermediate states ($0.81$–$0.86$), and a single high state ($\approx 0.89$). This structure refines rather than contradicts the three-state baseline: the five-state classification partitions the moderate-risk zone into three sub-regimes of varying stress intensity, while the low- and high-state assignments correspond directly to the low-risk and high-risk regimes of the three-state model.

By contrast, the $K=4$ specification fails identification across all $\kappa$ values consistent with the data’s empirical persistence (Table 13), indicating that there is no intermediate four-regime structure between the coarse three-regime classification and the finer five-state subdivision. $K=3$ is the smallest identified specification and $K=5$ the next, with $K=4$ representing a structural discontinuity.

The cross-regime analysis in Section 3.4 is conducted at $K=3$ for parsimony and economic interpretability, but the substantive findings are preserved if the analysis is repeated using the $K=5$ specification with the three intermediate states aggregated to the Moderate-Risk regime.

3.4. Correlation Analysis and Regime Conditioning

Table 18 reports three correlation measures between $A{R}_t$ and the interpretive benchmark series over the full sample period: Pearson (parametric), Spearman (rank-based), and Kendall’s tau (rank-based, less sensitive to ties and outliers). Each is reported with p-values testing $H_0:\rho =0$.

The three tests yield a consistent picture. None of the four series exhibits a statistically significant rank-based association with $A{R}_t$ under either Spearman or Kendall, with p-values uniformly above 0.25. The two Pearson coefficients that reach conventional significance, tracking error ($r=0.086$, $p<0.001$) and the MSCI World Index ($r=-0.061$, $p=0.016$), are not corroborated by either rank-based test. Since rank-based correlations are far less sensitive to a small number of extreme observations than Pearson, this divergence indicates that the parametric significance reflects the influence of a subset of high-leverage points rather than a continuous monotonic relationship spanning the full distribution.

Figure 9 provides a visual diagnostic of this pattern. The tracking error panel reveals that the Pearson association is driven primarily by a subset of high-risk observations exhibiting both elevated $A{R}_t$ and elevated tracking error, while moderate- and low-risk observations form a diffuse cloud with no discernible trend. This concentration of the association in one regime explains the divergence between Pearson, which weights extreme values heavily, and the rank-based tests, which treat all observations equally regardless of magnitude. The global benchmark panels confirm the absence of any continuous relationship: the full-sample OLS regression lines are essentially flat, and observations are well mixed across regimes at every level of $A{R}_t$. The most extreme negative benchmark returns ($-0.10$ to $-0.12$) correspond to low-risk rather than high-risk observations—an apparent paradox that reflects the 60-day rolling window construction of $A{R}_t$, which captures sustained co-movement over the preceding window rather than instantaneous market shocks.

Regime-Conditional Regression

The pattern in Figure 9 suggests that any genuine relationship between $A{R}_t$ and the interpretive series may operate primarily within the high-risk regime rather than uniformly across the sample. To test this directly, each interpretive series is regressed on $A{R}_t$, a High-Risk regime dummy $D_{\mathrm{High},t}$, and their interaction:

$$y_t={\beta }_0+{\beta }_{1}A{R}_t+{\beta }_2{D}_{\mathrm{High},t}+{\beta }_3(A{R}_t·D_{\mathrm{High},t})+{\epsilon }_t.$$

This specification fits two separate lines simultaneously: a calm-regime line with intercept ${\beta }_0$ and slope ${\beta }_1$, and a high-risk line with intercept ${\beta }_0+{\beta }_2$ and slope ${\beta }_1+{\beta }_3$. The interaction coefficient ${\beta }_3$ therefore measures the difference in slopes between the two regimes; a significant ${\beta }_3$ indicates a regime-dependent rather than continuous relationship. Standard errors are computed using the HAC (Newey-West) estimator with 7 lags to account for the strong autocorrelation in $A{R}_t$ induced by the 60-day rolling window, consistent with the autocorrelation structure implied by the transition matrix persistence estimates in Table 16.

Figure 10 visualises the regime-conditional fits reported in Table 19, replacing the single full-sample OLS line of Figure 9 with two regime-specific lines per panel. The geometric distance between the two slopes in each panel corresponds directly to the interaction coefficient ${\beta }_3$.

The regression results in Table 19 and the corresponding visualisation in Figure 10 reinforce three points. First, the non-high-risk slope ${\beta }_1$ is small and statistically indistinguishable from zero for all four dependent variables, confirming the absence of a continuous relationship between $A{R}_t$ and the interpretive series outside crisis periods. This is visible in Figure 10 as the near-horizontal blue line across every panel. Second, the interaction coefficient ${\beta }_3$ for tracking error is positive ($+0.037$), with the sign and rough magnitude predicted by the cluster-level amplification analysis in Section 3.4.1. Although ${\beta }_3$ falls short of conventional significance under the conservative HAC adjustment ($p=0.157$), its magnitude is economically meaningful: a one-unit increase in $A{R}_t$ during high-risk periods is associated with a 3.7% point increase in average daily tracking error, an order of magnitude consistent with the 1.09× to 1.23× amplification ratios documented at the cluster level. The visual divergence between the calm and high-risk slopes in the tracking error panel of Figure 10 illustrates this regime-dependent amplification clearly. The limited statistical power of the daily-frequency test reflects the conservative HAC adjustment, the concentration of the effect in roughly one-third of observations, and the substantial idiosyncratic variation in tracking error across the 40 funds. Third, no global benchmark exhibits a regime-dependent relationship with $A{R}_t$: all interaction coefficients are small, of inconsistent sign, and far from significance, with $R^2$ values below $0.006$. The near-parallel calm and high-risk slopes in the three benchmark panels of Figure 10 make this visually unambiguous. This rules out the alternative interpretation that systemic co-movement in the pension fund network is a passive reflection of global equity market conditions.

Together, the rank-based correlations and the regime-conditional regression resolve the apparent tension in the Pearson results: the parametric significance for tracking error reflects regime-specific amplification rather than a continuous linear effect, but the daily-frequency tests have limited power to confirm this directly. The cluster-level analysis in Section 3.4.1 provides a more powerful test of the same hypothesis by aggregating observations within structurally homogeneous fund groups under each regime, where the amplification effect is unambiguous.

3.4.1. Cluster-Level Tracking Error Amplification

Figure 11 quantifies the tracking error amplification at the cluster level, connecting the DTW clustering results of Section 3.2 with the HMM regime classification of Section 3.3. By aggregating observations within structurally homogeneous fund groups under each regime, this analysis tests directly whether the regime-dependent amplification suggested by Section 3.4 is present in the data.

All eight clusters exhibit positive amplification ratios, confirming that tracking error increases universally during high-risk periods regardless of cohort type or provider. The magnitude differs systematically between cluster types. Middle-cohort clusters (C1, C2, C7, C8) show higher absolute tracking errors in both regimes, reflecting their greater equity exposure. Cluster C2 (Provider 1, AG2 - AG7) exhibits the highest absolute tracking error ($6.3\times {10}^{-3}$during high-risk periods) with an amplification ratio of $1.18\times$, while Cluster C7 (Provider 4, AG2 - AG7) records the highest amplification ratio of $1.23\times$, tracking error grows by nearly a quarter during systemic stress.

Extreme-cohort clusters (C3 - C6, grouping AG1 and AG8/TIPF) show lower absolute tracking errors but still amplify during stress, with ratios ranging from $1.09\times$ (C6, Provider 3 extreme cohorts) to $1.17\times$ (C5, Provider 4 extreme cohorts). The $1.09\times$ ratio for C6, the lowest in the sample, indicates that even the most conservative cluster is not immune to systemic stress, merely less sensitive to it. This has direct welfare implications: TIPF funds serve participants actively drawing down retirement income, for whom any tracking error amplification translates immediately into income uncertainty at a life stage when recovery is not possible.

4. Discussion

The results in Section 3 point to a network whose risk architecture is shaped less by what providers choose to do than by what regulation requires them all to do in similar ways. We organise the discussion around four findings and what they imply for supervision.

4.1. When “Diversification” Stops Being Diversification

Across the seven-year sample, the absorption ratio never falls below 0.73. Even in the calmest 60-day windows we observe, a single common factor explains nearly three-quarters of the variance across 40 funds. By Proposition 3, this corresponds to an average pairwise correlation above 0.72, which is hard to reconcile with the textbook picture of a diversified system.

The clustering results tell us why. Funds group by age cohort, not by provider. Cluster C8 pulls together 12 middle-cohort funds from Providers 2 and 5, and Cluster C4 does the same for their extreme cohorts. That the three clustering methods agree exactly at $K_c=7$ ($\mathrm{ARI}=1$) is reassuring: this is not an artefact of one algorithm or one distance metric.

The mechanism is regulatory. Lifecycle rules tell every provider to allocate similarly within a given birth-year cohort, so participants who choose between providers are, in effect, choosing between portfolios that look almost identical at the cohort level. The same architecture that protects an individual saver from extreme allocation choices removes provider differentiation as a system-level risk-mitigation channel. The cross-provider merging of Providers 2 and 5 at $K_c^{*}=8$ (Clusters C8 and C4, Table 5) illustrates the extreme case: from a systemic risk perspective, these two providers constitute a single functional node, which amplifies the network concentration finding beyond what fund count alone would suggest. The economic reading of a high systematic risk share, therefore, inverts in this setting: standard portfolio theory treats it as a sign that idiosyncratic risk has been diversified away, but in a pension network, it is a sign that the cross-section has been compressed to the point of fragility. This reframing changes the regulator’s problem. The standard supervisory framework under IORP treats each fund as the unit of risk. The clustering result implies that the operative unit of risk is the cohort, not the fund: a shock affecting any middle-cohort fund simultaneously stresses all middle-cohort funds across all providers, regardless of which provider a participant chose. A participant who switches providers to reduce exposure is not reducing exposure at all. The robustness of this structure across all three clustering methods, with identical partitions at $K_c=7$ ($\mathrm{ARI}=1$, Section 3.2), confirms that the cohort-dominant pattern is a property of the data and not an artefact of the Ward linkage criterion chosen a priori in Section 2.3. Whether a supervisor addresses this by loosening cohort-level allocation constraints to restore provider differentiation, or by introducing system-level diversification requirements that operate on the network rather than the fund, is a normative question about whose risk the regulation is meant to manage. The empirical analysis cannot resolve it, but it makes the trade-off precise and previously invisible.

The risk of regulatory design producing precisely this kind of portfolio synchronisation is not confined to lifecycle allocation rules. In the context of the proposed IORP II revision, PensionsEurope has warned explicitly that performance benchmarking requirements may induce uniform portfolio alignment across providers: institutions seeking to avoid being flagged as outliers or subject to corrective measures may increasingly allocate toward reference benchmarks, generating concentration risks and reducing the scope for differentiated long-term strategies [63]. The clustering results in Section 3.2 provide empirical content for that warning in the Lithuanian second-pillar context: the cohort-dominant structure documented here is the ex post outcome of a regulatory design that already constrains cross-sectional variation, and the PensionsEurope concern is that the proposed benchmarking layer would deepen that constraint further. The present framework makes the cost of such deepening measurable: any regulatory change that increases portfolio similarity within cohorts would raise the absorption ratio floor above its current minimum of 0.728 and compress the already narrow window between the lower threshold ${\tau }_1=0.801$ and the crisis threshold $\tau =0.851$.

4.2. Elevated Co-Movement Is the Norm, Not the Exception

The HMM identifies three regimes. The low-risk regime ($\widehat{\mu }=0.773$) accounts for only 13.2% of trading days; the moderate-risk regime ($\widehat{\mu }=0.829$) for 54.1%; and the high-risk regime ($\widehat{\mu }=0.874$) for 32.7%. The long-run stationary distribution converges to the same proportions, confirming that monitoring frameworks calibrated on the assumption that calm is the default state will systematically underestimate the frequency and persistence of risk in this sector: the network’s baseline condition is already elevated, not one that occasionally enters stress.

Two features of the regime structure deserve emphasis beyond the occupancy shares. First, regimes are sticky: once in the high-risk state, the expected duration is approximately 65 trading days (Table 14), meaning systemic stress episodes in this network do not resolve quickly. Second, the transition matrix is asymmetric: direct escalation from low- to high-risk occurs with probability 0.0046, so the system almost always passes through the moderate-risk band before a full crisis. This asymmetry is operationally significant. It means the moderate-risk regime is not a neutral middle state but the zone in which a supervisor has the most time to act. A sustained rise in $A{R}_t$ above the lower threshold ${\tau }_1=0.801$ is a structural precursor signal, not a coincident one, and the transition probabilities quantify how much time typically elapses before escalation completes. Prior HMM applications to equity markets typically find low-risk states dominating the unconditional distribution; the finding here that low-risk accounts for only 13.2% of trading days is genuinely anomalous and reflects a network whose baseline condition is already elevated rather than one that occasionally enters stress.

The empirical crisis threshold $\tau =0.851$ is substantially above the range commonly used in the systemic risk literature [47]. This is simply the homogeneity finding stated in different units: in a network whose average pairwise correlation never drops below 0.72, the line between ordinary stress and crisis has to sit higher. Importing a generic threshold here would label roughly 70% of the sample as ‘crisis’ and would not be informative.

4.3. The Network Is Not Just Tracking Global Markets

A natural alternative explanation for the elevated absorption ratio is that it simply mirrors global equity market stress transmitted through fund holdings. The regime-conditional regressions rule this out. Slopes of $A{R}_t$ on the S&P 500, MSCI Europe, and MSCI World are statistically indistinguishable from zero in both calm and high-risk regimes, with $R^2$ below 0.006. The decoupling is visible in Figure 8: through late 2022 and into 2023, global benchmarks stabilise and recover while the network remains in high-risk.

The cluster-level tracking-error result reinforces this picture. The daily-frequency interaction coefficient is positive and economically meaningful ($+0.037$) but, as expected with HAC adjustment and a small effective sample within the high-risk regime, does not reach conventional significance. Aggregating to the cluster level removes that ambiguity: every one of the eight clusters amplifies during high-risk periods, with ratios from $1.09\times$ (C6, the most conservative) to $1.23\times$ (C7, the equity-heavy middle cohort). Whatever initiates a high-risk episode, what sustains it is internal to the network: amplification is universal across all eight clusters, not concentrated in the most equity-exposed funds. The practical implication is direct: supervisory tools designed to monitor externally transmitted shocks, global VaR limits, benchmark drawdown alerts, and MSCI factor exposure reports are largely blind to the risk documented here. A regulator relying exclusively on market-linked surveillance would have observed recovering global benchmarks through late 2022 and early 2023 while the pension network remained in the high-risk regime. The monitoring gap is not a calibration problem that better parameters would fix; it is structural, because the signal is in the cross-sectional co-movement of fund returns rather than in any market index. Closing it requires a network-internal indicator of the kind the absorption ratio provides.

4.4. What the Framework Makes Possible

The principal operational advance of this paper is that systemic risk monitoring for a pension fund network can now be conducted using only publicly available net asset value data, with no proprietary holdings information, no bilateral exposure records, and no labelled crisis history. The Bank of Lithuania publishes daily NAV data for all second-pillar funds; the absorption ratio and regime probabilities can be computed from these with a one-day lag and updated as a standing dashboard. This was not possible with previously available unsupervised frameworks, which either required balance-sheet data unavailable at daily frequency or relied on external benchmarks to define regime boundaries.

The framework’s credibility for regulatory use rests on a property that most HMM applications to financial data cannot claim: the regime structure is identified by the data rather than by modelling choices. The crisis threshold $\tau =0.851$ is derived from the emission distributions, not imposed; the number of states is selected by combining BIC with a Bayesian identifiability analysis under two independent prior calibrations; and the EM and Bayesian point estimates agree to within 0.001 across all three regimes (Table 15). A supervisor implementing this framework does not need to defend a discretionary tuning decision, because no such decision was made.

One identifiability result deserves separate attention. The $K=4$ specification fails identification across the full range of empirically supported $\kappa$ values, while both $K=3$ and $K=5$ are cleanly identified. The absence of a coherent four-regime partition, with $K=3$ and $K=5$ separated by a structural gap at $K=4$, is not a modelling failure but an empirical property of the absorption ratio’s dynamics in this network. It implies that the risk surface has a coarse three-level architecture, with finer gradations available only within the moderate-risk band. This is itself a finding about how systemic risk escalates in this sector and one that warrants theoretical attention in future work.

4.5. Implications for Supervision and Regulation

Three operational implications follow directly from the empirical results.

For supervisory monitoring, the framework produces two actionable thresholds from the data. The lower threshold ${\tau }_1=0.801$ separates low- from moderate-risk and, given the asymmetric transition structure (Table 16), constitutes a structural precursor to crisis: a sustained absorption ratio above this level should trigger an enhanced monitoring requirement and increased reporting frequency. The crisis threshold $\tau =0.851$ marks the boundary of the high-risk regime and should trigger active supervisory engagement. Both thresholds can be computed daily from publicly available NAV data published by the Bank of Lithuania, with no proprietary data requirement.

For lifecycle regulation, the cohort-dominant cluster structure reveals a tension between micro-prudential and macro-prudential objectives that current rules do not resolve. Lifecycle allocation constraints protect individual participants from age-inappropriate risk, but they simultaneously eliminate provider choice as a system-level diversification channel. Two regulatory instruments could address this. The first is a loosening of cohort-level allocation bands to allow more provider variation within cohort, restoring some cross-sectional differentiation. The second is a system-level diversification requirement imposed at the network level rather than the fund level, analogous to concentration limits in banking supervision. Which instrument is preferable depends on whose risk the regulation is primarily meant to manage, a question that is normative rather than empirical, but one that the current analysis makes precise for the first time.

For participant welfare, the $1.09\times$ amplification observed for AG8/TIPF funds is the finding that most directly contradicts standard risk disclosures for conservative vehicles. Participants drawing down retirement income during a high-risk episode face elevated income uncertainty at a life stage where intertemporal recovery is not available. Risk disclosures and guarantee structures for capital-preservation funds should reflect this residual systemic sensitivity, which is invisible to fund-level volatility metrics, but measurable in the network framework developed here.

5. Conclusion

This paper establishes that systemic risk in a regulated pension fund network is not imported from global markets but generated internally by the network’s own structure. The finding reframes a regulatory design question: lifecycle allocation rules that protect individual participants simultaneously compress the cross-section of fund behaviour to the point where provider choice offers no meaningful risk differentiation at the system level. A monitoring framework that does not account for this internal dynamic will be blind to the dominant source of risk in the sector.

Applied to 40 Lithuanian second-pillar pension funds over January 2019 - September 2025, the framework delivers a finding that inverts the causal direction assumed by most supervisory frameworks. The prior implicit model is that external shocks propagate into pension funds through their market exposures. The result here is that the network generates its own high-risk states, independently of global equity conditions, and spends the majority of its time in them. The absorption ratio has no statistically significant relationship with global equity benchmarks in either calm or high-risk regimes. The high-risk state persists through periods of global market recovery, confirming that network stress is not resolved by external stabilisation. And tracking-error amplification during stress is universal across all eight fund clusters: even the most conservative funds serving retirement-age participants amplify by $1.09\times$, at a life stage where intertemporal recovery is unavailable. The source of the vulnerability is structural homogeneity induced by lifecycle allocation regulation, not market exposure.

Three actors face distinct implications. The Bank of Lithuania can implement the monitoring framework immediately: the absorption ratio and regime probabilities are computable from the NAV data it already publishes, and the thresholds ${\tau }_1=0.801$ and $\tau =0.851$ provide operationally precise triggers. The European Insurance and Occupational Pensions Authority (EIOPA) is positioned to extend the framework to cross-jurisdictional comparison: if lifecycle regulation is the operative mechanism, cohort-dominant cluster structures and internally generated high-risk regimes should appear in any EU second-pillar system with comparable lifecycle allocation mandates, a testable prediction that the pipeline can evaluate using publicly available NAV data from any member state. The Lithuanian pension regulator faces the normative question that the empirical analysis cannot answer: whether to restore provider differentiation by loosening cohort allocation bands or to introduce network-level diversification requirements that operate above the fund. The analysis makes the trade-off visible and measurable for the first time.

By identifying which episodes matter, which threshold separates regimes, and which structural feature drives co-movement, the framework makes four extensions precise rather than merely plausible. Cross-jurisdictional comparison is the most direct: the pipeline requires only publicly available NAV data, so applying it to other EU second-pillar systems would test whether the cohort-dominant structure and internally generated risk regimes documented here are specific to Lithuania or a general consequence of lifecycle regulation. Portfolio-level decomposition is now motivated by a precise question: the framework has identified which episodes are worth examining; holdings data would reveal which asset classes drive the leading eigenvalue during high-risk periods and whether the fire sale dynamics hypothesised by [5] are the operative mechanism. Real-time applicability is now a specific empirical question: the forward state probabilities in Table 17 make quantitative predictions for the six months following September 2025, which subsequent data will allow to be evaluated against the regime sequence. And the $K=4$ identification failure is an unexplained structural feature of the absorption ratio’s dynamics in this network that warrants theoretical attention: why a three-level architecture with a structural gap at four states should emerge from a lifecycle-regulated system is a question the current empirical framework identifies but cannot answer. Cross-jurisdictional replication is the natural first step toward an answer: if cohort-dominant cluster structures and $K=4$ gaps appear consistently in other EU second-pillar systems subject to comparable lifecycle mandates, that would suggest the 3-level architecture is a general property of the regulatory design rather than a feature of Lithuanian implementation, a distinction that would directly inform EIOPA’s approach to system-level risk surveillance across member states.

The broader question this paper opens is whether the tension between micro-prudential protection and macro-prudential homogeneity is specific to lifecycle-regulated pension systems or a structural feature of any regulatory design that standardises portfolio behaviour across a large institutional population. Sovereign wealth funds, insurance general accounts, and target-date retail funds all operate under mandates that constrain cross-sectional variation. If the mechanism identified here generalises, the implication is that diversification at the participant level and fragility at the network level are not accidental co-occurrences but two consequences of the same regulatory logic and that making one better without making the other worse requires instruments that operate at the network level rather than the fund level. That is the design problem this paper makes legible.