1 Introduction↩︎

Value-at-Risk (VaR) and Expected Shortfall (ES) are central measures of financial tail risk in both academic research and regulatory practice. Considering the loss of a portfolio or investment, VaR characterizes only a quantile of the distribution of this loss, while ES summarizes expected losses beyond that quantile and therefore captures the severity of extreme downside outcomes. Moreover, ES is a coherent risk measure, satisfying properties like subadditivity, monotonicity, positive homogeneity, and translation invariance [1]. In contrast, VaR is generally not subadditive, implying that portfolio risk measured by VaR may exceed the sum of individual risks. These advantages make ES more consistent with the diversification principles of portfolio theory and capture the tail risk in a more comprehensive manner, explaining its growing role in modern risk management [2]¹.

Although conceptually attractive, ES remains challenging to forecast in practice. A fundamental obstacle is that ES on its own is not elicitable, whereas VaR is [9]. In other words, there exists no loss function for which ES is the unique minimizer of expected loss. This lack of direct elicitability partly explains why the development of ES forecasting methods has lagged behind that of VaR. A major breakthrough is provided by [10], which introduce a class of consistent scoring rules, known as Fissler-Ziegel (FZ) loss functions, to utilize the joint elicitability of VaR and ES. This breakthrough leads to a growing literature on joint VaR-ES forecasting. For example, [11] employ FZ loss functions in semi-parametric settings based on the generalized autoregressive score (GAS) model and the generalized autoregressive conditional heteroskedasticity (GARCH) model; [12] exploits the link between the asymmetric Laplace likelihood and the FZ loss for joint VaR-ES estimation; [13] further extend the Laplace-based approach to multivariate regression settings. Until now, most existing work on joint VaR-ES forecasting remains confined to parametric or semi-parametric models. In such frameworks, the number of parameters is kept modest to ensure interpretability, and each asset is often modeled separately with its own set of parameters. While tractable, their performance depends critically on correct model specification, where the true data-generating process for tail risks may be far more complex than a GARCH or GAS model can capture.

Fortunately, the era of big data brings new opportunities for tail risk modeling. Risk forecasters now have access to decades of financial data, covering tens of thousands of assets, along with hundreds of firm characteristics. This abundance of information naturally raises a key question: Can we develop a large tail risk model, analogous to large language models (e.g., ChatGPT, Gemini, and Claude) in natural language processing, that learns from rich and diverse data to predict VaR and ES more accurately? In other words, instead of relying on asset-specific (semi-)parametric models, could a single high-capacity universal model trained on massive data serve as a general engine for joint VaR-ES forecasting?

To shed light on large tail risk models, we first develop a unified learning framework for joint VaR-ES forecasting. This framework assumes a universal function indexed by unknown parameters, mapping from predictors (e.g., firm characteristics) to two latent but unconstrained risk scores. Then, it transforms these two risk scores into valid VaR and ES forecasts through a fixed one-to-one mapping. This data-driven framework enables linear regression models, machine learning methods, and deep learning architectures to be treated in the same manner, with parameters estimated by minimizing the FZ loss.

Under this framework, we propose a new large tail risk model, the retrieval-enhanced self-grouping autoencoder (ReSGA), designed to fully exploit complex financial information for joint VaR-ES forecasting. ReSGA first learns latent group structures across assets, allowing tail-risk information to be shared among assets with similar financial behavior. It then incorporates a retrieval mechanism that enables each asset to draw on relevant long-term historical “experiences” from both its own past and other assets within the same group, thereby enriching usable temporal information. By leveraging cross-sectional dependence and long-term temporal dynamics in a flexible data-driven manner, ReSGA exploits spatial-temporal information more effectively than existing approaches.

Empirically, we study the monthly VaR and ES of United States (US) equity returns using a comprehensive dataset provided by [14], covering more than 40,000 stocks from January 1926 to December 2023. Each stock-month observation is associated with 153 firm characteristics that are widely used in the asset pricing literature. We evaluate the out-of-sample VaR and ES forecasting performance of ReSGA and its 12 competing models through the average loss. From this analysis, we find that ReSGA consistently attains the lowest out-of-sample average loss across all stock universes considered. Besides the loss analysis, we further assess the statistical performance of all models via several statistical hypothesis tests, including the Diebold–Mariano (DM; [15]) test, the model confidence set (MCS; [16]) test, the conditional coverage (CC; [17]) test, and the auxiliary expected shortfall regression (AESR; [18]) test. The out-of-sample testing results also provide consistent evidence in favor of ReSGA.

Beyond statistical assessments, we are interested in whether the tail risk forecasts generated by ReSGA can lead to economically meaningful gains. Following [19], we construct a left-side momentum strategy by sorting stocks into deciles based on the predicted VaR or ES. The resulting VaR- and ES-sorted portfolios exhibit sharp nonlinearities concentrated in the extreme deciles: Stocks with the highest predicted tail risk (i.e., the lowest VaR or ES) deliver substantially lower returns and markedly worse downside risk, while a long-short strategy that buys the lowest-risk decile and sells the highest-risk decile, yields economically sizable Sharpe ratios. These patterns indicate that ReSGA’s tail risk forecasts capture cross-sectional variation in downside risk, which is also priced in the cross-section. Building on this insight, we propose a new size-enhanced left-side momentum signal that combines predicted ES with firm size. This signal is motivated by the “too big to fail” phenomenon in the US market: Tail risk in large firms is more likely to reflect compensated systematic risk, whereas similar tail risk in small firms often stems from idiosyncratic fragility. Based on this signal, we propose long-short portfolios for each considered model and find that all proposed non-econometric portfolios exhibit significant alphas relative to the Fama–French five-factor model. Hence, it indicates that our size-enhanced left-side momentum signal could broaden the mean-variance frontier. More importantly, our portfolio performances confirm that ReSGA delivers the best economic performance among all competing models, with the highest Sharpe ratio for the long-short decile portfolio and a clear monotonic return pattern across deciles. All of the above findings show that the forecasting advantage from ReSGA is not only statistically relevant but also economically meaningful.

Taken together, our statistical and economic findings above highlight that models can deliver better tail risk forecasts by exploiting richer information in the input. We view the informativeness of model input as one aspect of data complexity. This phenomenon connects to a broader debate on the virtue of complexity. In artificial intelligence, a growing body of evidence shows that the out-of-sample predictive performance of models often improves systematically with increases in either data complexity (with respect to in-sample data availability) or model complexity (with respect to model parameter size); see, for example, [20]. This empirical regularity, known as the scaling law, underpins the success of large language models in the last five years. In finance, there has been an increasing focus on the investigation of model complexity, while the examination of data complexity remains very rare. [21] provide early theoretical and empirical support for the virtue of model complexity in asset pricing, demonstrating that simple models with few parameters can severely understate return predictability relative to more complex alternatives. See more empirical evidence for this view in [22] and [23]. However, at the same time, other studies also question whether such gains from model complexity are economically meaningful or robust in asset pricing applications [24]–[27]. In terms of risk forecasting, evidence on the role of model complexity is far more limited; [28] show that machine learning models outperform linear models in realized volatility forecasting by capturing nonlinear temporal dynamics. However, to date, the role of model complexity and data complexity in tail risk forecasting remains largely unexplored.

To fill the gap, we further investigate the scaling performance of ReSGA and its competitors by varying model parameter size or in-sample data size to assess the role of model complexity or data complexity, respectively. First, our experiment results show that, holding the in-sample data fixed, simply increasing parameter size alone does not yield consistent performance improvements: Larger models do not reliably outperform smaller ones, and the best results in terms of out-of-sample average loss and portfolio Sharpe ratio typically occur at intermediate scales (around or below the in-sample sample size). In particular, models that exploit spatial–temporal inputs, such as ReSGA, benefit more from additional parameters and consistently outperform simpler architectures once sufficient parameters are available. Second, our complementary experiments vary the available in-sample data size while keeping the model architectures fixed. The corresponding results demonstrate that, in terms of out-of-sample average loss and portfolio Sharpe ratio, the performance of each model exhibits a clear improvement trend with the size of the in-sample data used. Together, the aforementioned scaling results indicate that our out-of-sample performance gains are primarily driven by richer data rather than by parameter growth alone. Hence, they provide limited support for a general virtue of model complexity, instead emphasizing the central role of data complexity in tail risk forecasting.

Lastly, we illustrate the usefulness of ReSGA through a group-importance analysis and a transfer-learning analysis. In the group-importance analysis, we explore economic model interpretability by investigating which types of firm characteristics drive tail risk predictability in the ReSGA model. It turns out that the predictive power for VaR and ES is concentrated in a small set of economic themes (or groups): Value, Low Risk, Momentum, Quality, and Short-Term Reversal. Additionally, the importance of these groups varies systematically across firm size and over time, reflecting shifts in the underlying economic sources of tail risk. In the transfer-learning analysis, we deploy the ReSGA model trained only on US equity data, without any re-estimation or fine-tuning, to five major international equity markets: China, Japan, the United Kingdom (UK), Australia and Canada. Despite substantial cross-market heterogeneity, ReSGA consistently presents reliable out-of-sample VaR-ES forecasting performance in every market, indicating that its predictive advantage is not US-specific and remains robust under market shifts. Hence, ReSGA, as a universal large tail risk model, has great generalizability in learning VaR and ES. At the portfolio level, however, return patterns differ substantially across countries. In China and Japan, higher size-enhanced left-side momentum is associated with higher returns, while in the UK, Australia, and Canada, large return dispersion concentrates in the extreme deciles. These findings suggest that tail-risk pricing varies across markets, likely reflecting differences in market structure, investor behavior, and risk premia.

The remaining paper proceeds as follows. 2 introduces a unified VaR-ES learning framework and the ReSGA model. 3 presents the empirical results for the US equity market, while 4 examines the transferability of the ReSGA model trained on US data to other markets. 5 concludes. Technical details are provided in Appendix 6 and other appendices in the supplementary materials.

2 Methodology↩︎

Let \(r_{i,t}\) denote the return of asset \(i\) at time \(t\), where \(t = 1, \ldots, T\) and \(i = 1, \ldots, N_t\), with \(N_t\) representing the number of available assets at time \(t\). The variables of interest in this paper are the conditional \(\tau\)-th VaR and ES of \(r_{i,t}\), defined respectively as \[\begin{align} \mathrm{VaR}_{i,t} = F_{i,t}^{-1}(\tau) \,\,\, \text{and} \,\,\, \mathrm{ES}_{i,t} = E \left( r_{i,t} \mid r_{i,t} \le \mathrm{VaR}_{i,t}, \mathcal{F}_{t-1} \right), \end{align}\] where \(\tau\in(0,1)\) is the quantile level, \(F_{i,t}(\cdot)\) is the conditional distribution of \(r_{i,t}\) given \(\mathcal{F}_{t-1}\), and \(\mathcal{F}_{t-1}\) is the information set containing all available information up to time \(t-1\).

In the following subsections, 2.1 introduces a unified learning framework for \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\). 2.2 presents our new ReSGA model. 2.3 describes a set of representative benchmark tail risk models that are used for comparison in our empirical study.

2.1 A Unified Learning Framework↩︎

Following the seminal work of [29], a common way to learn VaR is based on the check loss function. However, as shown in [9], ES alone is not elicitable, meaning that it cannot be the unique minimizer of an expected loss function. [10] solve this issue by showing that VaR and ES are jointly elicitable: \((\mathrm{VaR}, \mathrm{ES})\) is the unique minimizer of the general FZ loss function. Specifically, to jointly learn the VaR and ES of a univariate random variable \(Y\), [10] define the general FZ loss function (also termed as the scoring rule) as follows: \[\begin{align} \label{general95FZ} \begin{aligned} \ell(Y, (v, e)) &= \big(\boldsymbol{1}\{ Y \le v\} - \tau\big)g_1(v) - \boldsymbol{1}\{Y \le v\} g_1(Y) \\ &\quad + g_2(e)\left[e-v+\frac{1}{\tau}\boldsymbol{1}\{ Y \le v\} (v - Y)\right] - \tilde{g}_2(e), \end{aligned} \end{align}\tag{1}\] where \(\tau\in(0, 1)\) is the quantile level, \(\boldsymbol{1}\{\cdot\}\) is the indicator function, and \(g_1,g_2:\mathbb{R} \to \mathbb{R}\) are two pre-determined functions, with \(g_2\) being strictly increasing and convex, and \(\tilde{g}_2\) being its primitive. In particular, when \(g_1(x) = 0\) and \(g_2(x) = -1/x\), the general FZ loss function in (1 ) becomes the commonly used degree-0 FZ loss function: \[\begin{align} \label{eq:FZ0} \ell_{\mathrm{FZ0}}(Y, (v, e)) = -\frac{1}{\tau e}\,\boldsymbol{1}\{Y \le v\}(v - Y) + \frac{v}{e} + \log(-e) - 1, \end{align}\tag{2}\] where \((v, e)\in \Gamma\) for an admissible region \(\Gamma=\{(v, e): e < v < 0\}\). Using the degree-0 FZ loss function, a joint learning of VaR and ES can be achieved by noting \[\begin{align} \label{es95optim} (\mathrm{VaR} of Y,\, \mathrm{ES} ofY)=\mathop{\mathrm{argmin}}_{(v,e)\in\Gamma} \, \mathbb{E} \left[\ell_{\mathrm{FZ0}}(Y, (v, e)) \right]. \end{align}\tag{3}\]

As shown in 3 , the learning of \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\) under \(\ell_{\mathrm{FZ0}}\) needs to meet a key admissible condition: \(\mathrm{ES}_{i,t} < \mathrm{VaR}_{i,t} < 0\). To achieve this goal, we propose a new unified learning framework for \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\). Specifically, let \(\boldsymbol{y}_{i,t} = (y_{i,t,1}, y_{i,t,2})'\in\mathbb{R}^2\) represent a vector of two latent risk scores \(y_{i,t,1}\) and \(y_{i,t,2}\). Then, we define a Softplus-based one-to-one mapping \(p\): \(\boldsymbol{y}_{i,t} \rightarrow (\mathrm{VaR}_{i,t}, \mathrm{ES}_{i,t})\), satisfying \[\begin{align} \label{activate} \begin{aligned} \mathrm{VaR}_{i, t} & = -\mathrm{Softplus}(y_{i,t,1}), \\ \mathrm{ES}_{i, t} & = - \big[\mathrm{Softplus}(y_{i,t,1}) + \mathrm{Softplus}(y_{i,t,2})\big], \end{aligned} \end{align}\tag{4}\] where \(\mathrm{Softplus}(\cdot)\) is a function [30] defined as \[\begin{align} \mathrm{Softplus}(x) = \log (1 + e^x) \in \mathbb{R}^{+}, \,\,\, \text{for} \,\,\, x \in \mathbb{R}. \end{align}\] Since both \(\mathrm{Softplus}(y_{i,t,1})\) and \(\mathrm{Softplus}(y_{i,t,2})\) are strictly positive, the mapping in 4 guarantees that \(\mathrm{ES}_{i,t} < \mathrm{VaR}_{i,t} < 0\), regardless of the values of \(y_{i,t,1}\) and \(y_{i,t,2}\). Consequently, we are able to learn \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\) through a model for \(\boldsymbol{y}_{i,t}\), the form of which can be specified without imposing restrictions on \(y_{i,t,1}\) and \(y_{i,t,2}\).

Following the above idea, our learning framework only needs to specify a model for \(\boldsymbol{y}_{i,t}\) indexed by unknown parameters \(\boldsymbol{\theta}\). Therefore, unless otherwise stated, all tail risk models in this paper are designed for \(\boldsymbol{y}_{i,t}\). By writing \(\boldsymbol{y}_{i,t} \equiv \boldsymbol{y}_{i,t}(\boldsymbol{\theta})\), we have \((\mathrm{VaR}_{i,t}, \mathrm{ES}_{i,t})=p(\boldsymbol{y}_{i,t}(\boldsymbol{\theta}))\) according to (4 ), so we can make use of result (3 ) to estimate \(\boldsymbol{\theta}\) by minimizing the empirical degree-0 FZ loss: \[\begin{align} \label{eq:train95objective} \widehat{\boldsymbol{\theta}} & = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}} \sum_{t=1}^{T} \sum_{i=1}^{N_t} \ell_{\mathrm{FZ0}}\!\big(r_{i,t}, \big( \mathrm{VaR}_{i,t}, \mathrm{ES}_{i,t} \big) \big) = \mathop{\mathrm{argmin}}_{\boldsymbol{\theta}} \sum_{t=1}^{T} \sum_{i=1}^{N_t} \ell_{\mathrm{FZ0}}\!\big(r_{i,t}, p(\boldsymbol{y}_{i,t}(\boldsymbol{\theta}))\big). \end{align}\tag{5}\] In practice, due to the massive data volume, we use the adaptive moment estimation (Adam) algorithm in [31] to solve optimization problem 5 . Given \(\widehat{\boldsymbol{\theta}}\), \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\) are then learned by \(\widehat{\mathrm{VaR}}_{i,t}\) and \(\widehat{\mathrm{ES}}_{i,t}\), respectively, where \((\widehat{\mathrm{VaR}}_{i,t}, \widehat{\mathrm{ES}}_{i,t})=p(\boldsymbol{y}_{i,t}(\widehat{\boldsymbol{\theta}}))\).

2.2 The ReSGA Model↩︎

Let \(\boldsymbol{Y}_{t} \in \mathbb{R}^{N_t \times 2}\) denote the matrix of cross-sectional risk scores at time \(t\), where its \(i\)-th row is \(\boldsymbol{y}_{i,t}'\). Under the learning framework in Section 2.1, we design a new ReSGA model for \(\boldsymbol{Y}_{t}\), which aims to capture the nonlinear relationships between asset characteristics and risk scores, while simultaneously accounting for spatial and temporal dependencies of assets. To facilitate the construction of ReSGA, we let \(\boldsymbol{X}_{i,t-1} \in \mathbb{R}^{S \times P}\) be the feature matrix of asset \(i\) up to time \(t-1\), with its \(s\)-th row \(\boldsymbol{x}_{i,t-1-S+s}'\), where \(\boldsymbol{x}_{i,t} = (x_{i,t,1}, \ldots, x_{i,t,P})'\in\mathbb{R}^{P}\) contains \(P\) characteristics of this asset at time \(t\). Here, \(P\) is the number of characteristics and \(S\) is the number of time lags. Collecting all individual feature matrices, we obtain the tensor of asset characteristics \(\mathcal{X}_{t-1} = [\boldsymbol{X}_{1,t-1}, \dots, \boldsymbol{X}_{N_t,t-1}] \in \mathbb{R}^{N_t \times S \times P}\).

Using \(\mathcal{X}_{t-1}\) as the input, ReSGA, parameterized by \(\boldsymbol{\theta}=(\boldsymbol{\phi}',\boldsymbol{\psi}')'\), learns \(\boldsymbol{Y}_{t}\) through three core modules, encoder, retriever, and decoder, which are defined as follows: \[\begin{align} \label{resga95structure} \begin{aligned} (\mathcal{H}_t, \mathbb{M}_t) &= \mathrm{Encoder}(\mathcal{X}_{t-1}; \boldsymbol{\phi}), \\ \boldsymbol{Z}_t &= \mathrm{Retriever}(\mathcal{H}_{t}; \mathbb{M}_t), \\ \boldsymbol{Y}_t &= \mathrm{Decoder}\left(\mathcal{H}_t, \boldsymbol{Z}_t; \mathbb{M}_t, \boldsymbol{\psi}\right), \end{aligned} \end{align}\tag{6}\] where \(\mathrm{Encoder}(\cdot; \boldsymbol{\phi})\), \(\mathrm{Retriever}(\cdot; \mathbb{M}_t)\), and \(\mathrm{Decoder}(\cdot, \cdot; \mathbb{M}_t, \boldsymbol{\psi})\) denote the functional forms of encoder, retriever, and decoder, respectively (see their detailed architectures in Sections 6.1–6.3). Generally speaking,

• \(\mathrm{Encoder}(\cdot; \boldsymbol{\phi})\), parameterized by \(\boldsymbol{\phi}\), extracts a temporal feature \(\mathcal{H}_t\in \mathbb{R}^{N_t \times S \times D}\) for all assets across lags in an autoregressive manner, along with a set \(\mathbb{M}_t\) that captures latent group structures among assets according to their similarities in characteristics, where the dimension \(D\) serves as a user-specific hyperparameter controlling model complexity;

• \(\mathrm{Retriever}(\cdot;\mathbb{M}_t)\) leverages the learned group structure \(\mathbb{M}_t\) to retrieve long-term historical information not only from each asset’s own temporal history but also from other assets within the same group, yielding the retrieval feature \(\boldsymbol{Z}_t \in \mathbb{R}^{N_t \times D}\);

• \(\mathrm{Decoder}(\cdot, \cdot; \mathbb{M}_t, \boldsymbol{\psi})\), parameterized by \(\boldsymbol{\psi}\), hierarchically aggregates the group-, retrieval-, and asset-level information to generate the risk score \(\boldsymbol{Y}_t\).

It is worth noting that \(\mathbb{M}_t\) depends on the specific choice of \(\mathcal{X}_{t-1}\) and loss function. In this study, under the degree-0 FZ loss in (5 ), \(\mathbb{M}_t\) reveals the dynamic topological relationships among assets for tail risk analysis, guided by the similarities in asset characteristics. When ReSGA is applied to other studies, different inputs and loss functions can lead to different interpretations for \(\mathbb{M}_t\).

2.3 Competing Models↩︎

To assess the performance of ReSGA, we consider its 12 competing models. These competitors are classified into four categories: Point-wise, Temporal, Spatial–temporal and Econometric models, depending on whether they exploit temporal or cross-sectional information and whether they employ firm characteristics. Models in the first three categories are based on our learning framework in Section 2.1, while models in the last category learn VaR and ES using existing econometric approaches. For more details about these 12 competing models, one can refer to Appendix 7 in the supplementary materials.

2.3.0.1 Point-wise Models

The point-wise models aim to capture the relationship between asset characteristics and tail risks, but ignore both temporal and cross-sectional dependencies. They directly map the most recent asset characteristics to risk scores as follows: \[\begin{align} \label{point-wise} \boldsymbol{y}_{i,t} = f(\boldsymbol{x}_{i,t-1}; \boldsymbol{\theta}), \end{align}\tag{7}\] where \(\boldsymbol{x}_{i,t-1} \in\mathbb{R}^{P}\) is the vector of \(P\) characteristics for asset \(i\) at time \(t-1\), and \(f(\cdot;\boldsymbol{\theta})\) is a parametric function indexed by \(\boldsymbol{\theta}\). Particularly, we focus on two models in (7 ): (i) the classical linear regression model, labeled as “Linear”; (ii) a three-hidden-layer neural network model proposed by [32], labeled as “NN”.

2.3.0.2 Temporal Models

Unlike point-wise models, the temporal models explicitly account for the dynamic evolution of asset characteristics over time by assuming \[\begin{align} \label{temporal95model} \boldsymbol{y}_{i,t} = g(\boldsymbol{X}_{i,t-1}; \boldsymbol{\theta}), \end{align}\tag{8}\] where \(\boldsymbol{X}_{i, t - 1}\in\mathbb{R}^{S\times P}\) is the feature matrix of asset \(i\) at time \(t-1\), and \(g(\cdot;\boldsymbol{\theta})\) is a parametric function indexed by \(\boldsymbol{\theta}\). Note that the hyperparameter \(S\) specifies the look-back time period of characteristics. When \(S=1\), \(\boldsymbol{X}_{i,t-1}'\) reduces to \(\boldsymbol{x}_{i,t-1}\), which is the vector of the most recent characteristics (i.e., the input to point-wise models in (7 )). Compared with point-wise models, the temporal models in (8 ) exploit not only the predictive power of the latest characteristics, but also the historical evolution of characteristics, allowing them to capture dynamic patterns that may enhance predictive performance.

The temporal models require specific designs to effectively capture temporal dynamics. In this paper, we employ several representative temporal deep learning models in 8 , including:

1. two extended models: a lag-augmented neural network model that flattens \(\boldsymbol{X}_{i,t-1}\) into a vector as input (labeled as “LANN”), and a decomposition-based neural network model proposed by [33] that separately builds trend and seasonal components (labeled as “DLinear”);

2. two recurrent neural network models: long short-term memory network model from [34] (labeled as “LSTM”) and gated recurrent unit network model from [35] (labeled as “GRU”);

3. three transformer-type neural network models: the “Informer” model proposed by [36], which leverages an attention-based encoder and decoder to capture temporal dependencies, along with its encoder-only and decoder-only variants (labeled as “EInformer” and “DInformer”, respectively).

2.3.0.3 Spatial-temporal Models

For the spatial-temporal models, they simultaneously capture both the temporal evolution of asset characteristics and the cross-sectional interactions among assets using advanced deep learning architectures. Specifically, they are defined as \[\begin{align} \label{spatial95temporal95model} \boldsymbol{Y}_{t} = m(\mathcal{X}_{t-1}; \boldsymbol{\theta}), \end{align}\tag{9}\] where \(\mathcal{X}_{t-1} \in \mathbb{R}^{N_t \times S \times P}\) is the tensor of asset characteristics, and \(m(\cdot;\boldsymbol{\theta})\) is a parametric function indexed by \(\boldsymbol{\theta}\). Compared to the temporal models, which handle assets separately, the spatial–temporal models allow each asset’s risk scores to depend not only on its own historical information but also on that of other assets. By explicitly modeling these inter-asset dependencies, the spatial–temporal models are capable of capturing common shocks, contagion effects, and network spillovers, which are particularly relevant in financial markets. This spatial information sharing is expected to lead to improved robustness and predictive accuracy, especially during periods of heightened market co-movement or systemic stress.

Given the rich information contained in \(\mathcal{X}_{t-1}\), the spatial–temporal models must be carefully designed to efficiently exploit both temporal and cross-sectional dependencies. Besides the proposed ReSGA model, we consider another spatial–temporal model in 9 : the self-grouping autoencoder model from [37] (labeled as “SGA”). In short, the SGA model can be viewed as a simplified version of ReSGA. Both models share a common encoder–decoder architecture that learns time-varying group structures across assets. However, SGA removes all retrieval-related components and focuses solely on contemporaneous spatial dependencies. As a result, it could perform well when the temporal length \(S\) is relatively short. In contrast, ReSGA extends this framework by incorporating a retrieval mechanism to capture long-term spatial–temporal interactions, making it suitable for longer sequences and more complex systemic risk patterns.

2.3.0.4 Econometric Models

In the econometrics literature, [11] apply two widely used semi-parametric models, GAS and GARCH, to study the dynamics of VaR and ES. These models do not rely on machine learning or deep learning techniques to utilize the information of asset characteristics, yet they need to impose the constraint of \(\mathrm{ES}_{i,t}<\mathrm{VaR}_{i,t}<0\) when implementing optimization in model estimation. Including them allows us to assess whether more complex models deliver meaningful improvements over well-established econometric approaches. Notably, unlike the aforementioned non-econometric models, which employ a unified learning framework with a single parameter set shared across assets, the GAS and GARCH models are estimated separately for each asset to learn VaR and ES directly.

1 summarizes the features of all considered tail risk models, such as nonlinearity from deep learning, temporal dynamics from lagged characteristics, and spatial information sharing within cross-sectional input. In addition, it also distinguishes models by whether they support long-history inputs (i.e., capturing long memory effect). By clarifying the incremental modeling capabilities across point-wise, temporal, spatial–temporal, and econometric models, this table clearly illustrates how ReSGA integrates all these features within a single model to exploit available information efficiently.

3 Empirical Study of US Equity↩︎

3.1 Data↩︎

This section analyzes a comprehensive US equity dataset compiled by [14]. This dataset covers more than 40,000 stocks, providing a broad and representative cross-section of the US equity market for large-scale risk forecasting. Meanwhile, it spans the period from January 1926 to December 2023 at a monthly frequency, offering nearly a century of observations across different market regimes, including economic expansions and major crisis episodes. In this dataset, each stock in every month is associated with 153 firm characteristics². As in [32], we rank-normalize all characteristics cross-sectionally to mitigate the impact of outliers, and then impute missing characteristic values using contemporaneous cross-sectional medians. After this preprocessing treatment, we let \(r_{i,t}\) denote the excess return of stock \(i\) in month \(t\), and construct the characteristics tensor \(\mathcal{X}_{t-1}\) from all 153 characteristics.

Throughout the analysis, we apply the ReSGA model to learn \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\) at the quantile level \(\tau=0.05\), with a comparison to the 12 competing models outlined in Section 2.3. To evaluate out-of-sample forecasting performance for all models, we follow [32] to adopt an expanding-window process, with the last ten years of the sample (2014–2023) reserved for out-of-sample evaluation. Specifically, each model is trained on data from 1926 to 1995, validated on data from 1996 to 2013, and tested on data from 2014; this process then shifts forward by one year (trained on 1926–1996, validated on 1997–2014, and tested on 2015), and it is repeated until reaching the end of the out-of-sample period³. For the detailed implementation of the above training, validation, and testing process, one can refer to 9.

3.2 Statistical Performance Evaluation↩︎

We first evaluate the out-of-sample forecasting performance of ReSGA and its 12 competitors from a statistical perspective. As in [11], we consider the following out-of-sample average loss as a natural assessment criterion: \[\begin{align} \label{oos95loss} \ell_{\text{oos}} = \frac{1}{\sum_{t \in \text{oos}} |\mathbb{S}_t|}\sum_{t \in \text{oos}} \sum_{i \in \mathbb{S}_t} \ell_{\mathrm{FZ0}}\left(r_{i,t}, \left(\widehat{\mathrm{VaR}}_{i,t}, \widehat{\mathrm{ES}}_{i,t}\right) \right), \end{align}\tag{10}\] where “oos” denotes the out-of-sample period, \(\mathbb{S}_t\) denotes a group of stocks at time \(t\), with \(|\mathbb{S}_t|\) being its cardinality, \(\ell_{\mathrm{FZ0}}\) is the degree-0 FZ loss function defined in 2 , and \(\widehat{\mathrm{VaR}}_{i,t}\) and \(\widehat{\mathrm{ES}}_{i,t}\) are forecasts of \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\) from each learned model.

2 reports the results of \(\ell_{\text{oos}}\) across models and different stock groups. From this table, we can draw the following interesting findings:

• ReSGA achieves the lowest loss across the full stock group, as well as within five size-based stock groups: mega, large, small, micro, and nano. This consistent dominance indicates that ReSGA provides a broadly effective engine for VaR and ES forecasting, with its combination of learned cross-sectional group structure and retrieval-based temporal augmentation.

• Two spatial-temporal models tend to outperform seven temporal models, which in turn generally outperform two point-wise models. A notable exception is LANN, which performs worse than NN in most cases despite having richer inputs. By flattening the \(S \times P\) characteristics matrix into a single vector, LANN discards temporal ordering and therefore fails to model sequential dependence directly. As a result, it fails to capture dynamic patterns such as persistence or volatility clustering, illustrating that exploiting temporal information requires appropriate model architecture. Nevertheless, the overall ranking among all eleven universal models provides strong empirical support for the virtue of data complexity: Improvements from point-wise to temporal models reflect the value of incorporating temporal information, while further gains from temporal to spatial-temporal models highlight the importance of cross-sectional information.

• All eleven universal models substantially outperform two econometric models: GAS and GARCH. A natural explanation is that GAS and GARCH are primarily designed to capture return dynamics and largely overlook rich firm-level information, which is crucial for VaR and ES forecasting. Note that GAS and GARCH are estimated separately for each asset, whereas all universal models are trained on the pooled samples across assets. Hence, the above advantage of universal models over econometric models also points to a clear data scaling effect, highlighting the advantage of larger effective sample sizes and richer information content.

• Within the universal model family, moving beyond linear specifications yields clear gains, since nonlinear models (e.g., NN, LSTM, SGA, and ReSGA) deliver lower losses than the linear model in most cases. This finding suggests that nonlinearities and interaction effects among firm characteristics play an important role in tail risk forecasting.

• In terms of size-based stock groups, a clear monotonic pattern emerges: Larger firms systematically exhibit lower losses. This pattern is expected, as larger firms have greater liquidity, more standard published financial statements, and more stable trading environments, all of which contribute to their better predictability of tail risk.

In addition to the out-of-sample average loss, we further evaluate the out-of-sample performance of all models using four statistical hypothesis tests: the DM in [15], MCS in [16], CC in [17], and AESR in [18]. See their detailed descriptions in 10 of the supplementary materials. In terms of the values of \(\ell_{\text{oos}}\) with \(\mathbb{S}_t\) being the full stock group, the DM test aims to compare the prediction performance of any two models, and the MCS test searches for a batch of models that perform significantly better than the others. The CC and AESR tests provide statistical evidence of the validity of VaR and ES forecasts, respectively. Notably, we do not adopt the e-backtesting approach of [38] to examine the validity of ES forecasts, since that framework is mainly designed for sequential monitoring and non-fixed-sample sizes, whereas our empirical exercise is based on fixed-sample out-of-sample evaluation.

3 reports the pairwise DM test statistics for all considered models. From this table, we have two clear findings. First, except for SGA, ReSGA dominates all other models with significantly lower values of \(\ell_{\text{oos}}\). Second, all universal models deliver statistically significant improvements over econometric models.

Next, starting from the full set of models, our MCS test iteratively eliminates inferior models at the 90% confidence level, until only those statistically indistinguishable models remain to form the final model set. It turns out that only SGA and ReSGA are retained in the final model set. This outcome indicates that no other forecasting models perform as well as SGA and ReSGA in a statistical sense, underscoring the importance of spatial-temporal modeling for tail risk forecasting.

Moreover, we apply the CC and AESR tests to check the validity of VaR and ES forecasts for each stock at the significance level \(\alpha \in \{0.01, 0.05, 0.10\}\), respectively, and then denote the proportion of stocks having valid VaR and ES forecasts as the pass rate. 4 reports pass rates with respect to VaR and ES forecasts across all models. From this table, we find that, except for ES at the level \(\alpha=0.01\), the two spatial-temporal models, SGA and ReSGA have the largest values of pass rates. In particular, ReSGA shows a substantially larger pass rate than SGA for the prediction of ES at the level \(\alpha=0.05\) or \(0.10\).

Overall, the above statistical performance evaluations consistently demonstrate the state-of-the-art performance of ReSGA in tail risk forecasting, which in turn implies the virtue of data complexity.

3.3 Economic Performance Evaluation↩︎

To explore the economic implications of the VaR and ES, it is important to evaluate all tail risk models from an economic perspective. Following [19], we first assess economic gains using decile portfolios sorted on the predicted VaR or ES. This trading strategy is referred to as “left-side momentum”. It guides us to monthly rebalance ten decile portfolios (P1–P10), after sorting all stocks in descending order according to their predicted VaR or ES at the end of each month. Note that our predicted VaR and ES are negative in this paper, so P1 corresponds to the lowest-risk portfolio and P10 represents the highest-risk portfolio.

To keep the presentation concise, 5 reports portfolio results only for ReSGA, as the performance for other models is qualitatively similar. As shown in this table, under VaR-based sorting, average portfolio returns remain relatively stable up to P8 (\(1.090\%\)), but drop sharply to \(0.619\%\) in P9 and turn negative in P10 (\(-0.284\%\)), while the Sharpe ratio falls to \(-0.085\) in P10. This decline in returns is accompanied by a marked increase in downside risk: Maximum drawdown worsens from \(0.499\) in P8 to \(0.730\) in P9 and \(0.807\) in P10. A similar pattern appears under ES-based sorting, with average returns decreasing from \(1.196\%\) in P8 to \(0.512\%\) in P9 and \(-0.024\%\) in P10, and maximum drawdown deepening to \(0.710\) in P9 and \(0.810\) in P10.

The above findings indicate P9 and P10, which include stocks with high predicted tail risk, exhibit a substantially worse performance in terms of return, Sharpe ratio, and drawdown. [19] attribute this pattern to behavioral underreaction to bad news and the persistence of left-tail risks, which concentrate losses in the most risk-exposed portfolios. Our findings mirror this stylized fact under both VaR- and ES-based sorting, supporting the view that tail-risk forecasts from the ReSGA can capture economically meaningful variation in downside risk across the cross-section.

Beyond the above left-side momentum, we further design a new size-enhanced left-side momentum strategy, which utilizes the trading signal motivated by the “too big to fail” phenomenon in the US stock market. Our key trading idea is to exploit the different economic implications of tail risk across firms of different sizes. Specifically, we assign strong buying signal to large-cap stocks with small values of ES, reflecting that downside risk in large firms often attracts market attention, policy support, or investor demand, and it thus can limit extreme losses and help preserve long-term value. In contrast, we assign strong selling signal to small-cap stocks with significantly negative ES, which often reflects idiosyncratic fragility, limited liquidity, or elevated default risk, rather than priced systematic risk.

Motivated by this idea, we construct the following trading signal: \[\begin{align} \label{Cap95ES} \alpha_{i,t} = \bigl(\mathrm{Cap}_{i,t-1} - \overline{\mathrm{Cap}}_{t-1}\bigr) \times \bigl[1 - \exp(\widehat{\mathrm{ES}}_{i,t})\bigr], \end{align}\tag{11}\] where \(\mathrm{Cap}_{i,t-1}\) denotes the log market capitalization of asset \(i\) at time \(t-1\), and \(\overline{\mathrm{Cap}}_{t-1}\) is the corresponding cross-sectional mean. Here, the term \([1 - \exp(\widehat{\mathrm{ES}}_{i,t})]\) lies between \(0\) and \(1\), with its value decreasing as \(\widehat{\mathrm{ES}}_{i,t}\) increases, serving as a proxy for tail-loss severity; the term \((\mathrm{Cap}_{i,t-1} - \overline{\mathrm{Cap}}_{t-1})\) tends to be positive for large-cap stocks and negative for small-cap stocks, acting as a classifier for firm size. Combining these two terms, the signal \(\alpha_{i,t}\) in (11 ) assigns high values to large-cap stocks with severe predicted tail risk, and low values to small-cap stocks with similar tail risk.

To assess the economic performance of our size-enhanced left-side momentum strategy, we adopt the signal \(\alpha_{i,t}\) to propose monthly rebalanced decile portfolios as before, with \(\widehat{\mathrm{ES}}_{i,t}\) computed from different models. Moreover, we regress the monthly returns of the resulting H–L portfolios on the Fama–French five factors ([39]) to examine whether there is an abnormal return (alpha) unexplained by standard factors. 6 reports the results of this analysis, and it delivers the following findings:

• The ReSGA-based strategy delivers the highest value of monthly average return, annualized Sharpe ratio and alpha for the H–L portfolio, outperforming all competing strategies. This result indicates that the best ES forecasting accuracy from ReSGA can translate into the strongest economic gains when combined with firm size information.

• Regardless of the choice of models, decile portfolios exhibit a clear and robust monotonic pattern: Moving from P1 to P10, average portfolio returns and annualized Sharpe ratios decline steadily and become strongly negative in the lowest deciles (P9 and P10). This pattern confirms that the signal \(\alpha_{i,t}\) can effectively rank stocks that have economically meaningful size-enhanced tail risk exposure.

• For all strategies based on the universal models, the \(p\)-values of the estimated alphas are below \(0.05\), whereas the corresponding \(p\)-values for the GAS- and GARCH-based strategies are \(0.051\) and \(0.128\), respectively. This contrast indicates that the econometric models exhibit relatively weak predictive power, highlighting the advantage of our proposed unified learning framework.

Taken together, 5 6 show that the superior advantage of ReSGA in learning tail risks can translate into economically meaningful gains.

4pt

3.4 Scaling Performance Evaluation↩︎

To examine whether there is virtue of model complexity in learning VaR and ES, we begin by evaluating scaling performance with respect to parameter size for all considered machine-learning-based models (i.e., all universal models excluding Linear). Specifically, we adopt a within-architecture scaling strategy: For each model, we vary the number of model parameters by adjusting width/depth-related hyperparameters, while keeping all other mechanisms unchanged⁴. This approach produces multiple variants of each model with substantially different numbers of learnable parameters. The connection between hyperparameter choices and parameter size for each model is reported in 15 of the supplementary materials.

7 details the out-of-sample average loss of different models when their parameter size varies over several orders of magnitude. It includes exact loss values, while its accompanying figure visualizes the corresponding scaling patterns⁵. Combining both pieces of evidence, we obtain three main interesting findings:

• ReSGA consistently outperforms all competitors once the parameter size reaches \(10^5\) or above, with SGA ranking second, whereas Informer achieves the lowest loss at smaller scales. This pattern is evident from the bold entries in the table and the crossing behavior in the figure: ReSGA and SGA improve sharply from small to intermediate scales, while Informer deteriorates once its parameter size becomes large. These findings indicate that models with richer spatial-temporal inputs, such as SGA and ReSGA, need sufficient capacity in parameters to exploit useful information from inputs, while temporal models can perform well with relatively few parameters. In contrast, the point-wise model of NN benefits little from expanded capacity, likely because its input information is inherently limited.

• There is no systematic evidence for a general virtue of model complexity. Except for LSTM, the parameter-richest specification (around \(10^7\) parameters) does not deliver the best performance within each model family. The figure makes this non-monotonicity especially visible: Most loss curves flatten or turn upward after intermediate scales. Even for LSTM, where loss performance improves with scale, the gains are modest. Meanwhile, LSTM with around \(10^7\) parameters still underperforms ReSGA with \(10^5\) parameters. This absence of monotonic and substantial improvement suggests that simply increasing model size alone does not guarantee better tail-risk forecasts.

• As shown by [20], in natural language processing, optimal performance is often achieved when model size is of the same order as the in-sample data size. This principle is only partially supported in our financial task. To be specific, ReSGA and SGA achieve their lowest out-of-sample losses at around \(10^6\) parameters, which is comparable to the effective training sample size (approximately \(2\times 10^6\)); however, this alignment is not universal, since other models show no clear match between optimal parameter size and sample size. These observations point to the heterogeneity in the virtue of data complexity: Spatial-temporal models could exploit richer inputs, allowing for more capacity in parameters to be used effectively, whereas temporal and point-wise models rely on less information and consequently obtain fewer benefits from increased parameterization.

Meanwhile, 8 reports the out-of-sample Sharpe ratios of H–L portfolios as parameter size varies, with an accompanying figure to visualize the trend of Sharpe ratio across parameter scales. From 8, we find that Sharpe ratio fails to exhibit a clear monotonic pattern as the parameter size increases. First, the best-performing parameter size differs across models: Most models peak at intermediate scales, while LSTM and Informer attain their highest Sharpe ratios at the largest scale. Second, the leading model also changes across parameter scales. These mixed patterns indicate that increasing model complexity alone does not reliably improve portfolio-level gains.

To further examine the virtue of data complexity in learning VaR and ES, we next evaluate scaling performance with respect to in-sample data size. Specifically, we hold each model fixed and vary the fraction of available in-sample data used for training and validation, from \(1\%\) to \(100\%\) of the full in-sample dataset. For each fraction, we draw the in-sample data uniformly at random from the full in-sample dataset, and keep the hyperparameter configuration fixed at the values selected using the full in-sample dataset. This design holds model complexity fixed, and it isolates the role of in-sample data size that is one key dimension of data complexity.

9 reports the out-of-sample average loss across different in-sample data sizes, while its accompanying figure visualizes how these losses evolve as more in-sample data become available. From this table and its accompanying figure, we can reach the following findings:

• All models benefit from more in-sample data (from \(1\%\) to \(100\%\)). With the exception of EInformer, each model attains its best performance when the full in-sample dataset is available. Meanwhile, the loss curves in the figure mostly present monotonic pattern and decline as the sample fraction increases.

• The best-performing model changes with the amount of in-sample data. When in-sample data are scarce (only \(1\%\), \(5\%\), or \(10\%\)), temporal models such as GRU and Informer outperform the other models. Once the fraction of used in-sample data reaches \(25\%\) or above, the spatial-temporal models (SGA and especially ReSGA) become dominant, indicating that these models require sufficient data to fully exploit cross-sectional sharing and long-term temporal dependencies.

• Point-wise model, NN, exhibits the largest improvement as the amount of in-sample data grows, consistent with its reliance on pooled cross-sectional information to compensate for weaker inductive structure. Nevertheless, even with the full in-sample dataset, NN remains clearly outperformed by the best temporal and spatial-temporal models. Thus, more data help broadly, with the greatest statistical gains when the model can use the richest in-sample information effectively.

Beyond evaluation from loss, 10 reports the out-of-sample Sharpe ratios of H–L portfolios when the in-sample data size varies, while its accompanying figure visualizes the scaling patterns of Sharpe ratio. Being broadly consistent with the loss results in 9, this table delivers the following notable results from the economic perspective:

• ReSGA delivers the highest Sharpe ratio across different models once the in-sample data size reaches at least \(25\%\) of the full dataset.

• The fitted curves show that portfolio performance generally improves as more in-sample data become available, although the pattern is not strictly monotonic for every model.

Taken together, the above results show that a larger in-sample dataset can lead to stronger economic performance, especially for models that can effectively exploit rich information.

Overall, our analysis results provide little support for the virtue of model complexity. Instead, they highlight the importance of data complexity: Performance gains arise mainly from models that can exploit richer information with the use of a larger in-sample data size, rather than a larger parameter size alone. While the available sample size is always limited in finance, this motivates the next scaling dimension: expanding the effective information content of the input with appropriate model specification and enough parameters, which appears to be a more effective path to improving tail risk forecasting accuracy than simply increasing model size.

3.5 Group Importance↩︎

A further empirical question concerns which firm characteristics are most important for forecasting tail risk. Due to the superior performance of ReSGA from both statistical and economic perspectives, we utilize this model to tackle the above question. Assessing importance at the level of individual characteristics, however, is challenging due to the well-known dilution problem: When many characteristics are highly correlated, their marginal contributions become spread across related variables, obscuring the true economic drivers of predictive performance [28]. To address this issue, we assess importance at the group level rather than at the level of individual characteristics. Specifically, we adopt the 13-category taxonomy of [14], which organizes the 153 firm characteristics into economically interpretable groups, including Low Risk, Value, Quality, Low Leverage, Momentum, Size, Profit Growth, Short-Term Reversal, Seasonality, Investment, Profitability, Debt Issuance, and Accruals. As argued by [28], grouping characteristics reduces redundancy among highly correlated signals and provides a clearer decomposition of predictive content. This approach allows us to evaluate the incremental value of entire economic themes, rather than attributing importance to individual characteristics that may proxy for similar information.

We quantify group importance via a drop-group procedure. To be specific, for each group (denoted as \(\mathbb{G}\)), we follow [32] to construct perturbed testing samples in which all characteristics in group \(\mathbb{G}\) are set to zero, while all remaining characteristics are left unchanged. By passing these perturbed samples through the trained model to generate new forecasts, we recompute the out-of-sample average loss in 10 on these new forecasts, with the set \(\mathbb{S}_t\) containing all examined stocks. Then, we measure the importance of group \(\mathbb{G}\) by the increase in loss of perturbed samples relative to that obtained using the original testing samples. After computing importance for all 13 groups, we set negative values to zero to avoid spurious attribution and normalize the remaining values to sum to one.

1 exhibits the ten most important groups from ReSGA over the out-of-sample period (2014.01–2023.12). From this figure, two clear and economically intuitive patterns emerge.

• The top five groups (Value, Low Risk, Momentum, Quality, and Short-Term Reversal) together account for more than 90% of importance for the full dataset. This concentration indicates that tail risk predictability is driven by a small set of groups rather than being evenly distributed across all characteristics. Similar findings can be found in the asset pricing literature ([32], [41]; [42]). Note that the concentration is even stronger for mega- and large-cap stocks, with top five groups receiving over \(99\%\) importance, whereas the top five groups in small-, micro-, and nano-stocks contribute \(91\%\), \(89\%\), and \(94\%\), respectively.

• The ranking of the leading groups shows both stability and systematic variation across different levels of market capitalization. Generally speaking, Low Risk and Momentum are consistently important regardless of the level of market capitalization. This is intuitive, as volatility- and beta-related characteristics in Low Risk are closely linked to tail risk, while Momentum captures persistent return dynamics that affect tail outcomes. Apart from these two common drivers, for mega- and large-cap stocks, Quality and Low Leverage are more influential, suggesting that tail risk variation among large firms is closely tied to balance-sheet strength and financial stability; in contrast, for small-, micro-, and nano-cap stocks, Value and Short-Term Reversal become more crucial, where valuation performance and short-horizon price reversals are expected to closely relate to distress risk, illiquidity, and left-tail outcomes of those stocks.

Moreover, 2 traces the month-by-month evolution of group importance for ReSGA, where we compute group importance using the increase of loss from each month. Its heatmap reinforces that Low Risk remains important almost throughout the entire out-of-sample period, consistent with its role documented from 1. As one may expect, the importance of Low Risk was especially pronounced from April to September 2020, since the COVID-19 shock had severely disrupted global financial markets. In contrast, Value and Quality display a visible “complementary” pattern over time, where the elevated period of Quality importance often coincides with muted Value, and vice versa. This finding suggests that the dominant predictors of tail risk shift between valuation-driven repricing forces and balance-sheet resilience/financial strength. Additionally, Quality and Low Leverage exhibit similar dynamics in certain months, a similarity that is expected due to their overlapping economic content related to profitability, safety, and leverage constraints.

4 Transfer Learning↩︎

Fundamentally, our universal models are driven by a principle that tail risk is governed by a mapping structure from asset characteristics to risk scores. As long as relevant asset characteristics are available, the model, by construction, can be applied to different assets without re-estimation. This contrasts sharply with econometric models, such as GARCH and GAS, which are inherently asset-specific. Hence, it raises a natural question: Can a universal model, trained on the largest and most information-rich US equity market, generalize to other equity markets without re-estimation?

To address this question, we conduct a transfer learning exercise. To be specific, we train each universal model exclusively on US equity dataset as in Section 3 and then apply it, without any re-training or fine-tuning, to five major international equity markets: China, Japan, UK, Australia, and Canada. This exercise allows us to assess whether the universal models capture fundamental relationships between firm characteristics and future tail risk across different institutional, regulatory, and liquidity environments.

11 reports the out-of-sample average loss across universal models and international equity markets. From this table, we have the following noteworthy findings:

• Despite substantial cross-market heterogeneity, ReSGA delivers the lowest out-of-sample loss in Japan, UK, and Australia, and remains highly competitive in China (ranked fourth) and Canada (ranked second). This strong and stable performance indicates that the predictive structure learned by ReSGA from US data generalizes well across international equity markets, even without re-estimation.

• In contrast, SGA performs poorly in most non-US markets, except for Canada. In particular, it is often dominated by simpler temporal models. This pattern suggests that cross-sectional dependence learned in the source market does not necessarily transfer well to other markets, potentially diluting asset-level tail-risk signals and introducing noise. The relatively strong performance of SGA in Canada is plausibly attributable to the close economic integration and high similarity between the US and Canadian equity markets, which makes cross-sectional structures learned from US data more transferable in this case. Note that ReSGA mitigates the above limitation of SGA by augmenting cross-sectional grouping with a retrieval mechanism that extends the effective temporal memory. The ability to exploit long historical information appears to compensate for potential misspecification in the cross-sectional structure, resulting in more robust transfer performance across markets.

• Temporal models, such as GRU, LSTM, and DLinear, also exhibit relatively stable out-of-sample performance across markets and often rank among the top predictors. This robustness further supports the view that temporal dynamics constitute a more transferable source of predictive power than cross-sectional dependence when models are deployed across markets without re-estimation.

Besides the loss analysis, 12 further displays the models in the final model set from the MCS test for the transfer learning across markets. For reference, we also include the earlier results for the US market, where only SGA and ReSGA remain in the final model set. The cross-market evidence from 12 delivers a clear message: ReSGA is the only model selected in all final model sets. This result indicates that the advantage of ReSGA is not confined to the US market, but generalizes robustly across international equity universes. Moreover, the CC and AESR testing results in 13 provide further evidence of the effectiveness of ReSGA in transfer learning: Across the five target markets, ReSGA delivers consistently strong pass rates for both VaR and ES forecasts and ranks among the top performers in nearly all cases.

Lastly, we examine the economic gains of ReSGA in the transfer learning exercise. 14 reports cross-market portfolio performance based on the ES-driven trading signal defined in 11 . Its results show that portfolio return patterns differ substantially between US and non-US markets. In China and Japan, the decile portfolios exhibit a reversed pattern, with higher predicted signal associated with higher returns, leading to negative H–L portfolio performance. In the UK, Australia, and Canada, the results are more heterogeneous. Although monotonicity is weaker, the extreme deciles display substantial return dispersion, especially in the lowest-signal portfolios, which drive large negative H–L returns. These findings suggest that the pricing of tail risk varies across markets, likely reflecting differences in market structure, investor behavior, and risk premia.

Taken together, the transfer learning results reinforce our central conclusions: ReSGA functions as a general engine for tail risk forecasting.

5 Conclusion↩︎

This paper studies tail risk forecasting, with a particular focus on VaR and ES, in the era of big data. Methodologically, we develop a unified joint VaR-ES learning framework supervised by the FZ loss. Based on this framework, we propose a large tail risk model, ReSGA, which captures the nonlinear relationships between asset characteristics and VaR-ES by accounting for spatial and temporal dependencies among assets. Empirically, using nearly a century of US equity data covering over 40,000 stocks and 153 firm characteristics, we find that ReSGA consistently delivers the best out-of-sample forecasting performance among a broad set of competing models. These forecasting gains are economically meaningful: Trading strategies based on ReSGA forecasts exhibit pronounced left-tail momentum and deliver superior long-short decile portfolio performances using a newly proposed size-enhanced left-side momentum signal.

Moreover, we provide systematic evidence on the scaling behavior of ReSGA and other machine-learning-based universal models. Our results offer little support for a general virtue of model complexity measured by parameter size alone. Instead, they favor the virtue of data complexity in learning tail risk, as improvements in forecasting VaR and ES can be achieved through larger in-sample datasets or more informative model inputs. In other words, the success of ReSGA is not only due to its access to rich available information but also to its ability to effectively extract, share, and utilize this information. Finally, we conduct group-importance analysis and transfer-learning experiments to illustrate the interpretability and cross-market generalizability of ReSGA.

6 The Architecture of ReSGA↩︎

6.1 Encoder↩︎

The encoder of ReSGA first applies the widely used recurrent neural network, GRU [35], to extract the temporal feature \(\mathcal{H}_{t}\). The GRU is designed to capture temporal dynamics among asset characteristics through its recurrent hidden states. Let \(\underline{\boldsymbol{x}}_{i, t-1}^s = \boldsymbol{x}_{i, t-1-S+s}\) denote the input characteristic vector of asset \(i\) at time lag \(s\). Then, the updates within the GRU are formulated as \[\begin{align} \label{gru} \begin{aligned} \text{(update gate)} \quad \boldsymbol{z}_{i,t}^s &= \mathrm{sigmoid}\left(\boldsymbol{W}_{z} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{z} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_{z}\right), \\ \text{(reset gate)} \quad \boldsymbol{r}_{i,t}^s &= \mathrm{sigmoid}\left(\boldsymbol{W}_{r} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{r} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_{r}\right), \\ \text{(candidate hidden state)} \quad \tilde{\boldsymbol{h}}_{i,t}^s &= \tanh \left(\boldsymbol{W}_{h} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{h}\left(\boldsymbol{r}_{i,t}^s \odot \boldsymbol{h}_{i,t}^{s-1}\right) + \boldsymbol{b}_{h} \right), \\ \text{(hidden state)} \quad \boldsymbol{h}_{i,t}^s &= \boldsymbol{z}_{i,t}^s \odot \boldsymbol{h}_{i,t}^{\,s-1} + \left(1 - \boldsymbol{z}_{i,t}^s\right) \odot \tilde{\boldsymbol{h}}_{i,t}^s, \end{aligned} \end{align}\tag{12}\] where \(\boldsymbol{W}_* \in \mathbb{R}^{D \times P}\), \(\boldsymbol{U}_* \in \mathbb{R}^{D \times D}\), and \(\boldsymbol{b}_* \in \mathbb{R}^D\) for \(* \in \{z,r,h\}\) are weight matrices and bias vectors shared across time lags. The hidden-state dimension \(D\) is a tunable hyperparameter that controls model complexity, and \(\odot\) denotes element-wise multiplication. The element-wise activation functions \(\mathrm{sigmoid}\) and \(\tanh\) are defined as \(\mathrm{sigmoid}(x) = \frac{1}{1 + e^{-x}}\) and \(\tanh(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}\) for \(x\in \mathbb{R}\).

In 12 , the hidden state \(\boldsymbol{h}_{i,t}^{s} \in \mathbb{R}^D\) serves as the key feature, summarizing the accumulated temporal information up to lag \(s\), while \(\tilde{\boldsymbol{h}}_{i,t}^s\) represents the candidate hidden state, used to store newly inputted information \(\boldsymbol{x}_{i, t-1}^s\) at lag \(s\). The update gate \(\boldsymbol{z}_{i,t}^{s} \in (0,1)^D\) controls how much of the past information from \(\boldsymbol{h}_{i,t}^{s - 1}\) should be kept, whereas the reset gate \(\boldsymbol{r}_{i,t}^{s} \in (0,1)^D\) determines the degree to which the previous hidden state \(\boldsymbol{h}_{i,t}^{s - 1}\) contributes to the candidate state \(\tilde{\boldsymbol{h}}_{i,t}^{s} \in (-1,1)^D\). Together, these gating mechanisms enable the GRU to adaptively balance memory retention and new information integration, effectively capturing complex temporal dependencies in an autoregressive manner.

For each asset \(i = 1, \dots, N_t\) and lag \(s = 1, \dots, S\), we obtain the hidden states \(\{\boldsymbol{h}_{i,t}^1, \dots, \boldsymbol{h}_{i,t}^S\}\) from the GRU, and then we stack them into the temporal feature across all assets and lags to construct \[\begin{align} \label{temporal95feature} \mathcal{H}_t = [\boldsymbol{H}_{1,t}, \dots, \boldsymbol{H}_{N_t,t}] \in \mathbb{R}^{N_t \times S \times D} \,\,\, \text{with} \,\,\, \boldsymbol{H}_{i,t} = [\boldsymbol{h}_{i,t}^1, \dots, \boldsymbol{h}_{i,t}^S] \in \mathbb{R}^{S \times D}. \end{align}\tag{13}\]

In addition to \(\mathcal{H}_t\), the encoder of ReSGA uses the most recent hidden states \(\{\boldsymbol{h}_{1,t}^S, \dots, \boldsymbol{h}_{N_t,t}^S\}\) to infer the latent group structure among assets. Specifically, it constructs a set \(\mathbb{M}_t\) that captures inter-asset relationships based on the similarity of temporal features. To fulfill this, we first treat each asset as a potential group center, with its node feature defined as \[\begin{align} \boldsymbol{u}_{k,t} = \boldsymbol{h}_{k,t}^S, \,\,\, \text{for} \,\,\, k = 1, \dots, N_t. \end{align}\] Then, we measure the similarity between asset \(i\) and center \(k\) by the cosine similarity: \[\begin{align} \gamma_{i,k,t} = \text{Cosine} \left( \boldsymbol{u}_{k,t}, \boldsymbol{h}_{i,t}^S\right) = \frac{\boldsymbol{u}_{k,t}' \boldsymbol{h}_{i,t}^S}{\|\boldsymbol{u}_{k,t}\| \, \|\boldsymbol{h}_{i,t}^S\|}, \,\,\, \text{for} \,\,\, i, k = 1, \dots, N_t, \end{align}\] where \(\|\cdot\|\) denotes the \(L_2\) norm. To mitigate the impact of weak or noisy similarities, we retain only the \(K^\ast\) strongest connections for each asset \(i\): \[\begin{align} \label{topK} \gamma^\ast_{i,k,t} = \begin{cases} 1, & \text{if } k = i, \\ \gamma_{i,k,t}, & \text{if } |\gamma_{i,k,t}| > \text{Top}(\mathbb{A}_{i,t}, K^\ast), \\ 0, & \text{otherwise}, \end{cases} \end{align}\tag{14}\] where \(K^\ast\) is a hyperparameter controlling the number of retained connections, \(\mathbb{A}_{i,t} = \{|\gamma_{i,k,t}| : k \neq i\}\), and \(\text{Top}(\mathbb{A}, b)\) denotes the \(b\)-th largest value in the set \(\mathbb{A}\).

After thresholding, we form each group \(\mathbb{M}_{k,t}= \{ i : \gamma^\ast_{i,k,t} \neq 0 \}\) by gathering all assets connected to center \(k\). Meanwhile, we discard those groups containing only their own center, as single-member groups carry no sharing information. Consequently, we reach the group structure: \[\begin{align} \label{group95structure} \mathbb{M}_t = \{\mathbb{M}_{k,t} : k = k_1, \dots, k_{K_t} \}, \end{align}\tag{15}\] where \(k_1, \dots, k_{K_t}\) denote the indices of all identified group centers.

Our data-driven way to discover group structure \(\mathbb{M}_{t}\) matches many empirical practices. In financial markets, assets such as stocks naturally form groups based on their shared characteristics and economic linkages. Typically, stocks within the same group \(\mathbb{M}_{k,t}\) tend to co-move due to common factors such as industrial background, financial statement, or ownership background. For instance, technology companies like NVIDIA, Apple, Microsoft, and AMD tend to cluster into the same group, as their prices are jointly affected by factors such as semiconductor supply chains, innovation cycles, and market sentiment toward the tech sector. These co-movements underscore the importance of modeling group-specific relationships to capture heterogeneity across different groups. For more empirical elaborations on the group structures in stock return forecasting, one can refer to [43] and references therein.

Combining 13 and 15 , we can equivalently write the workflow of encoder in a brief form \[\begin{align} \label{resga95encoder} (\mathcal{H}_t, \mathbb{M}_t) = \mathrm{Encoder}(\mathcal{X}_{t-1}; \boldsymbol{\phi}), \end{align}\tag{16}\] where \(\boldsymbol{\phi}\) includes all trainable parameters in 12 .

6.2 Retriever↩︎

Although the temporal feature \(\boldsymbol{h}_{i,t}^s\) in 12 summarizes the accumulated temporal information up to lag \(s\) (i.e. time \(t-S+s\)), it inevitably “forgets” long-term dependencies due to the effects of the update and reset gates. Moreover, our encoder only extracts the temporal feature for each asset independently, without incorporating useful temporal information from other related assets. For instance, during the 2022 technology sell-off, the price movement of NVIDIA exhibited a pattern strikingly similar to that of Apple during the 2008 financial crisis: Both experienced rapid growth followed by sharp corrections driven by macroeconomic tightening. Retrieving such long-ago but structurally similar patterns across assets can provide valuable insights that purely temporal models may overlook.

To address these limitations, we introduce the retriever of ReSGA, which captures structurally similar patterns from the historical trajectories of related assets. Specifically, we partition the \(S\) lags into \(L^\ast\) non-overlapping segments, forming \[\begin{align} \boldsymbol{H}_{i,t} = [\boldsymbol{h}_{i,t}^1, \dots, \boldsymbol{h}_{i,t}^{S/L^\ast}, \dots, \boldsymbol{h}_{i,t}^{2S/L^\ast}, \dots, \boldsymbol{h}_{i,t}^{(L^\ast-1)S/L^\ast}, \dots, \boldsymbol{h}_{i,t}^{S}], \end{align}\] where \(\boldsymbol{h}_{i,t}^{lS/L^\ast}\) for \(l = 1, \dots, L^\ast\) represents the historical context up to time \(lS/L^\ast\), and \(L^\ast\) is a user-specific hyperparameter satisfying that \(S\) is divisible by \(L^\ast\). The temporal features from the \(L^\ast\) segments thus act as reference checkpoints, each summarizing distinct stages of past information for subsequent retrieval and comparison.

Next, we retrieve information from its related assets defined by the group structure \(\mathbb{M}_t\). For each asset \(i\), the related assets are collected into a set \[\begin{align} \label{center95set} \mathbb{N}_{i,t} = \{k : i \in \mathbb{M}_{k,t}\}, \end{align}\tag{17}\] where \(\mathbb{M}_{k,t}\) is the group centered at asset \(k\). By construction, \(\mathbb{N}_{i,t}\) contains the \(K^{\ast}\) most similar assets to \(i\) identified by the encoder, indicating which assets’ historical trajectories should be searched in the retrieval stage.

We now define the retrieval feature for asset \(i\) at time \(t\) by aggregating information from a set of related assets. For each related asset \(j \in \mathbb{N}_{i,t}\), its contribution to the retrieval feature of asset \(i\) is given by \[\begin{align} \label{related95R} \boldsymbol{z}_{i,j,t} = \sum_{l = 1}^{L^\ast - 1} \alpha_{i,j,l,t}\, \boldsymbol{h}_{j,t}^{(lS/L^\ast) + 1} \in \mathbb{R}^{D}, \end{align}\tag{18}\] where \(\alpha_{i,j,l,t}\) is the normalized attention weight defined as \[\begin{align} \label{attn95weight} \alpha_{i,j,l,t} = \frac{\exp(\beta_{i,j,l,t})}{ \sum_{\ell = 1}^{L^\ast - 1} \exp(\beta_{i,j,\ell,t})}, \,\,\, \text{with} \,\,\, \beta_{i,j,l,t} = \mathrm{Cosine}\!\left(\boldsymbol{h}_{i,t}^S,\, \boldsymbol{h}_{j,t}^{lS/L^\ast}\right). \end{align}\tag{19}\]

Here, the attention weight \(\alpha_{i,j,l,t}\) quantifies the similarity between the current temporal feature of asset \(i\), \(\boldsymbol{h}_{i,t}^S\), and the checkpoint feature of asset \(j\), \(\boldsymbol{h}_{j,t}^{lS/L^\ast}\). This similarity measures how much the present feature of asset \(i\) resembles the past feature of asset \(j\) at segment \(l\). Importantly, while the attention weights are based on these checkpoint features, \(\boldsymbol{z}_{i,j,t}\) uses the subsequent feature \(\boldsymbol{h}_{j,t}^{(lS/L^\ast) + 1}\) that immediately follows each checkpoint⁶. The rationale is that if asset \(i\)’s current feature mirrors asset \(j\)’s historical feature at time \(lS/L^\ast\), then asset \(j\)’s next state after that point provides valuable predictive information for asset \(i\)’s future evolution. In this way, the retriever aligns historical features across assets and integrates their forward states, allowing the model to infer potential future behavior from similar past behaviors.

Subsequently, we aggregate the contributions from all related assets to form the retrieval feature. To fulfill this, we compute the average similarity between each related asset \(j\) and the target asset \(i\) across all checkpoints: \[\begin{align} \bar{\beta}_{i,j,t} = \frac{1}{L^\ast-1}\sum_{l=1}^{L^\ast - 1} \beta_{i,j,l,t}, \end{align}\] which captures the overall relevance of asset \(j\) to asset \(i\). Then, following 18 and 19 , we calculate the relative contribution of each related asset and aggregate it into the retrieval feature for each asset \(i\): \[\begin{align} \boldsymbol{z}_{i,t} = \sum_{j\in\mathbb{N}_{i,t}} \omega_{i,j,t} \boldsymbol{z}_{i,j,t} \in \mathbb{R}^{D}, \,\,\, \text{with} \,\,\, \omega_{i,j,t} = \frac{\exp(\bar{\beta}_{i,j,t})}{ \sum_{m \in \mathbb{N}_{i,t}}\exp(\bar{\beta}_{i,m,t})}. \end{align}\]

Finally, we stack all \(\boldsymbol{z}_{i,t}\) across assets to obtain the retrieval feature: \[\begin{align} \label{resga95retriever} \boldsymbol{Z}_{t} = \mathrm{Retriever}(\mathcal{H}_t; \mathbb{M}_t) = [\boldsymbol{z}_{1,t}, \dots, \boldsymbol{z}_{N_t,t}] \in \mathbb{R}^{N_t \times D}. \end{align}\tag{20}\] By design, \(\boldsymbol{Z}_{t}\) aggregates information through a two-stage attention mechanism. In the first stage, the model adaptively selects the most informative temporal segments from the history of each related asset, focusing on those that resemble the current behavior of asset \(i\). In the second stage, it re-weights the related assets based on their overall relevance, emphasizing those whose past dynamics are most predictive for asset \(i\). This retrieval structure allows the retriever to recall and integrate inter-asset temporal patterns that are similar to the present behavior of the target asset, thereby improving its ability to capture and forecast risk dynamics at time \(t\). See 3 for the entire process of retriever.

6.3 Decoder↩︎

The decoder integrates the temporal features \(\mathcal{H}_t\), group structure \(\mathbb{M}_t\), and retrieval features \(\boldsymbol{Z}_t\) to forecast the risk scores \(\boldsymbol{y}_t\). Specifically, it adopts a residual connecting framework ([44]), including four layers: (i) a group layer capturing shared group-level co-behavior, (ii) a retrieval layer utilizing the retrieval features, (iii) an individual layer modeling stock-specific behavior, and (iv) a fully connected layer producing the final risk score forecasts.

6.3.0.1 Group Layer

As the retriever captures spatial–temporal dependencies for each asset by linking its current temporal feature to the historical features of related assets, the group layer instead focuses on modeling contemporaneous interactions among assets within the same group. In other words, while the retriever aligns the present state of asset \(i\) with the past trajectories of its peers to recover long-term temporal relations, the group layer exploits the most recent temporal features across assets to capture shared short-term co-movements.

Specifically, the group layer aims to capture common influences within each group based on the most recent temporal feature \(\{\boldsymbol{h}_{1, t}^S, \dots, \boldsymbol{h}_{N_t,t}^S\}\) and group structure \(\mathbb{M}_t\). For each group \(\mathbb{M}_{k,t}\) in \(\mathbb{M}_t\), we first design a group-level feature by aggregating member features: \[\begin{align} \boldsymbol{G}_{k,t} = \mathrm{LeakyReLU}\!\left[\boldsymbol{W}_1 \left(\sum_{i \in \mathbb{M}_{k,t}} \gamma_{i,k,t}^{*} \boldsymbol{h}^S_{i,t} \right) + \boldsymbol{b}_1 \right] \in \mathbb{R}^D, \end{align}\] where \(\boldsymbol{h}^S_{i,t}\) is the temporal feature of asset \(i\) from 12 , \(\boldsymbol{W}_1 \in \mathbb{R}^{D \times D}\) and \(\boldsymbol{b}_1 \in \mathbb{R}^D\) are learnable parameters, and \(\mathrm{LeakyReLU}(\cdot)\) is an element-wise activation defined as \[\begin{align} \mathrm{LeakyReLU}(x) = \begin{cases} x, & x > 0, \\ 0.01 x, & x \leq 0. \end{cases} \end{align}\] This aggregation produces \(\boldsymbol{G}_{k,t}\) as a concise summary of the short-term shared dynamics within group \(\mathbb{M}_{k,t}\).

Next, we propagate the group-level effects to individual assets through a graph attention mechanism [45]: \[\begin{align} \boldsymbol{g}_{i,t} = \sum_{k \in \mathbb{N}_{i,t}} \alpha^\ast_{i,k,t}\boldsymbol{G}_{k,t}, \,\,\, \text{with} \,\,\, \alpha^\ast_{i,k,t} = \frac{\exp(\beta^\ast_{i,k,t})}{ \sum_{m \in \mathbb{N}_{i,t}}\exp(\beta^\ast_{i,m,t})}, \,\,\, \beta^\ast_{i,k,t} = \text{Cosine}(\boldsymbol{h}^S_{i,t}, \boldsymbol{G}_{k,t}), \end{align}\] where \(\mathbb{N}_{i,t}\) in 17 denotes the set of group centers containing asset \(i\). Here, \(\beta^\ast_{i,k,t}\) measures the similarity between the current temporal feature of asset \(i\) and the group-level feature of group \(k\), while \(\alpha^\ast_{i,k,t}\) provides normalized attention weights to indicate the relative importance of each group to asset \(i\).

In the final step of the graph layer, we feed the group feature \(\boldsymbol{g}_{i,t}\) into a two-headed fully connected network to produce the forward and backward components: \[\begin{align} \label{graph95forecast} \begin{aligned} \widetilde{\boldsymbol{g}}_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_2\boldsymbol{g}_{i,t} + \boldsymbol{b}_2), \\ \boldsymbol{h}^G_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_3\widetilde{\boldsymbol{g}}_{i,t} + \boldsymbol{b}_3), \\ \boldsymbol{o}^G_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_4\widetilde{\boldsymbol{g}}_{i,t} + \boldsymbol{b}_4), \end{aligned} \end{align}\tag{21}\] where \(\boldsymbol{W}_2, \boldsymbol{W}_3, \boldsymbol{W}_4 \in \mathbb{R}^{D \times D}\) and \(\boldsymbol{b}_2, \boldsymbol{b}_3, \boldsymbol{b}_4 \in \mathbb{R}^D\) are learnable parameters. Here, \(\boldsymbol{o}^G_{i,t}\) captures the group-level information directly used for forecasting \(\boldsymbol{y}_{i,t}\), while \(\boldsymbol{h}^G_{i,t}\) represents the portion of \(\boldsymbol{h}^S_{i,t}\) already explained by group co-behavior. This residual design separates the group effect from subsequent retrieval and individual layers, enhancing both model interpretability and learning efficiency [44].

6.3.0.2 Retrieval Layer

After the graph layer, the retrieval layer employs the retrieval feature \(\boldsymbol{Z}_{t}\) defined in 20 , which captures long-term spatial–temporal dependencies through a two-stage attention mechanism. To be specific, this retrieval layer passes \(\boldsymbol{Z}_{t}\) through a two-headed fully connected network, mirroring the structure used in the group layer: \[\begin{align} \label{ret95forecast} \begin{aligned} \widetilde{\boldsymbol{z}}_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_5\boldsymbol{z}_{i,t} + \boldsymbol{b}_5), \\ \boldsymbol{h}^R_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_6\widetilde{\boldsymbol{z}}_{i,t} + \boldsymbol{b}_6), \\ \boldsymbol{o}^R_{i,t} &= \mathrm{LeakyReLU}(\boldsymbol{W}_7\widetilde{\boldsymbol{z}}_{i,t} + \boldsymbol{b}_7), \end{aligned} \end{align}\tag{22}\] where \(\boldsymbol{z}_{i,t}\) is the retrieval feature of asset \(i\) (i.e., \(i\)-th row of \(\boldsymbol{Z}_t\)), and \(\boldsymbol{W}_5, \boldsymbol{W}_6, \boldsymbol{W}_7 \in \mathbb{R}^{D \times D}\) and \(\boldsymbol{b}_5, \boldsymbol{b}_6, \boldsymbol{b}_7 \in \mathbb{R}^D\) are learnable parameters. Here, \(\boldsymbol{h}^R_{i,t}\) serves as the backward component of retrieval layer, representing the portion of \(\boldsymbol{h}_{i,t}\) explained by retrieved information, while \(\boldsymbol{o}^R_{i,t}\) denotes the forward component of retrieval layer that contributes directly to forecasting \(\boldsymbol{y}_{i,t}\). The parameter sets \((\boldsymbol{W}_5,\boldsymbol{W}_6,\boldsymbol{W}_7)\) and \((\boldsymbol{b}_5,\boldsymbol{b}_6,\boldsymbol{b}_7)\) play analogous roles to \((\boldsymbol{W}_2,\boldsymbol{W}_3,\boldsymbol{W}_4)\) and \((\boldsymbol{b}_2,\boldsymbol{b}_3,\boldsymbol{b}_4)\) in the group layer, respectively.

6.3.0.3 Individual Layer

The individual layer models the idiosyncratic part of each asset after accounting for the effects already captured by the group and retrieval layers. Specifically, we define: \[\begin{align} \label{ind95forecast} \boldsymbol{o}^{I}_{i,t} = \mathrm{LeakyReLU}(\boldsymbol{W}_8 \boldsymbol{h}^I_{i,t} + \boldsymbol{b}_8), \,\,\,\text{with} \,\,\, \boldsymbol{h}^I_{i,t} = \boldsymbol{h}^S_{i,t} - \boldsymbol{h}^{G}_{i,t} - \boldsymbol{h}^{R}_{i,t}, \end{align}\tag{23}\] where \(\boldsymbol{W}_8 \in \mathbb{R}^{D \times D}\) and \(\boldsymbol{b}_8 \in \mathbb{R}^D\) are learnable parameters. The output (or the forward component of individual layer) \(\boldsymbol{o}^{I}_{i,t}\) captures asset-specific behavior unique to asset \(i\), ensuring that the residual temporal information unexplained by group or retrieval effects is explicitly modeled.

6.3.0.4 Fully Connected Layer

Finally, the decoder aggregates the outputs from the group, retrieval, and individual layers through additive composition to generate the final forecasts: \[\begin{align} \label{resga95decoder} \boldsymbol{Y}_t = [\boldsymbol{y}_{1,t}, \dots, \boldsymbol{y}_{N_t,t}] = \mathrm{Decoder}(\mathcal{H}_{t}, \boldsymbol{Z}_t; \mathbb{M}_t, \boldsymbol{\psi}) \in \mathbb{R}^{N_t \times 2}, \end{align}\tag{24}\] where each individual forecast is given by \(\boldsymbol{y}_{i,t} = \boldsymbol{W}_9 \left(\boldsymbol{o}^{G}_{i,t} + \boldsymbol{o}^{R}_{i,t} + \boldsymbol{o}^{I}_{i,t}\right) + \boldsymbol{b}_9\), with \(\boldsymbol{W}_9 \in \mathbb{R}^{2 \times D}\) and \(\boldsymbol{b}_9 \in \mathbb{R}^2\) as learnable parameters. The collection \(\boldsymbol{\psi}= \{(\boldsymbol{W}_j, \boldsymbol{b}_j)\}_{j=1}^{9}\) represents all trainable parameters within the decoder.

Intuitively, this hierarchical design in (24 ) allows the decoder to jointly capture group-level co-behavior, long-term spatial-temporal dependencies, and idiosyncratic behavior within a unified learning framework.

6.4 Summary↩︎

Accompanied with \(\mathrm{Encoder}(\cdot; \boldsymbol{\phi})\) in 16 , \(\mathrm{Retriever}(\cdot; \mathbb{M}_t)\) in 20 , and \(\mathrm{Decoder}(\cdot, \cdot; \mathbb{M}_t, \boldsymbol{\psi})\) in 24 , we get the specification of ReSGA in 6 . See Fig 4 for its entire network architecture.

7 Model Architectures↩︎

In this appendix, we present the details of all competing models considered in Section 2.3.

7.1 Point-Wise Models↩︎

We consider two point-wise models in (7 ): Linear and NN.

7.1.1 Linear↩︎

The Linear model assumes \(f\) in (7 ) satisfies \[\begin{align} \label{linear} \boldsymbol{y}_{i,t} = f(\boldsymbol{x}_{i,t - 1}; \boldsymbol{\theta}) = \boldsymbol{\alpha}+ \boldsymbol{\beta}' \boldsymbol{x}_{i,t - 1} , \end{align}\tag{25}\] where \(\boldsymbol{\alpha}\in \mathbb{R}^{2}\) and \(\boldsymbol{\beta}\in \mathbb{R}^{P \times 2}\) are the vectors of intercepts and regression coefficients, respectively. Here, \(\boldsymbol{\theta}\) includes all parameters in \(\boldsymbol{\alpha}\) and \(\boldsymbol{\beta}\). To mitigate the potential high-dimensional issue, we apply \(\ell_1\)-penalization on \(\boldsymbol{\theta}\) during training, promoting sparsity in the learned parameters.

7.1.2 NN↩︎

Following [32], the NN model applies a multi-layer fully connected neural network as the function \(f\) in (7 ): \[\begin{align} \tag{26} \boldsymbol{x}^{(l)}_{i, t - 1} &= \mathrm{ReLU} \left (\boldsymbol{W}^{(l)} \boldsymbol{x}^{(l - 1)}_{i, t-1} + \boldsymbol{b}^{(l)} \right) \in \mathbb{R}^{D^{(l)}}, \,\,\, \text{for}\,\,\, l = 1, \dots, L, \\ \tag{27} \boldsymbol{y}_{i,t} & = f(\boldsymbol{x}_{i,t-1}; \boldsymbol{\theta}) = \boldsymbol{W}^{(L + 1)} \boldsymbol{x}^{(L)}_{i, t-1} + \boldsymbol{b}^{(L + 1)} \in \mathbb{R}^2, \end{align}\] where \(\boldsymbol{x}_{i,t-1}^{(l)}\) denotes the output from the \(l\)-th hidden layer, with \(\boldsymbol{x}^{(0)}_{i,t-1}=\boldsymbol{x}_{i,t-1}\), and \(D^{(l)}\) is its dimension. Here, \(L\) is the number of hidden layers, the activation function ReLU is defined as \[\begin{align} \mathrm{ReLU}(z) = \max(0, z), \,\,\, \text{for}\,\,\, z \in \mathbb{R}, \end{align}\] the weight matrix \(\boldsymbol{W}^{(l)} \in \mathbb{R}^{D^{(l)} \times D^{(l-1)}}\) and bias vector \(\boldsymbol{b}^{(l)} \in \mathbb{R}^{D^{(l)}}\), for \(l=1, \dots, L+1\), are network parameters, with \(D^{(0)} = P\) and \(D^{(L+1)} = 2\), and \(\boldsymbol{\theta}\) contains all the parameters in \(\{\boldsymbol{W}^{(l)}, \boldsymbol{b}^{(l)}: l=1, \dots, L+1\}\).

From 26 , we can see that the hidden layers recursively apply linear transformations and nonlinear activation to extract predictive features from the raw characteristics, while the final linear projection in 27 is almost the same as the linear model in 25 . In particular, the dimensionalities \(\{D^{(l)}:l = 1, \dots, L\}\) control the “width” of NN, while the number of hidden layers \(L\) controls the “depth” of NN. Needless to say, they are crucial for the complexity of NN (i.e., the dimensionality of \(\boldsymbol{\theta}\)).

According to the geometric pyramid rule proposed by [46], an exponentially decaying sequence of network widths helps balance computational efficiency and nonlinear approximation capability. Specifically, this rule suggests \[\begin{align} D^{(l)} = \left\lfloor \gamma^{l - 1} D \right\rfloor, \,\,\, \text{for}\,\,\, l = 1, \dots, L, \end{align}\] where \(\lfloor \cdot \rfloor\) is the floor operator, and \(\gamma\) is a user-defined hyperparameter. [32] employ \(L = 3\), \(D = 32\), and \(\gamma = 0.5\) for \(P = 94\). Following this setup, we also fix \(L = 3\) and \(\gamma = 0.5\), but vary \(D\) to investigate the virtue of model complexity in tail risk forecasting. Additionally, as in [32], we also apply the batch normalization technique [47] to control the variability of network inputs across mini-batches. For more technical details on fully connected neural networks, one can refer to [48].

7.2 Temporal Models↩︎

We consider seven temporal models in (8 ): LANN, DLinear, LSTM, GRU, Informer, EInformer, and DInformer.

7.2.1 LANN↩︎

The LANN model extends the NN model in 7.1.2 by incorporating temporal information through lagged firm characteristics. Architecturally, LANN is identical to the NN, and is defined as \[\begin{align} \boldsymbol{x}^{(l)}_{i, t - 1} &= \mathrm{ReLU} \left (\boldsymbol{W}^{(l)} \boldsymbol{x}^{(l - 1)}_{i, t-1} + \boldsymbol{b}^{(l)} \right) \in \mathbb{R}^{D^{(l)}}, \quad l = 1, \dots, L, \\ \boldsymbol{y}_{i,t} & = f(\boldsymbol{x}_{i,t-1}; \boldsymbol{\theta}) = \boldsymbol{W}^{(L + 1)} \boldsymbol{x}^{(L)}_{i, t-1} + \boldsymbol{b}^{(L + 1)} \in \mathbb{R}^2, \end{align}\] where the input \(\boldsymbol{x}^{(0)}_{i,t-1}=\mathrm{vec}(\boldsymbol{X}_{i,t-1})\in\mathbb{R}^{SP}\) flattens the matrix \(\boldsymbol{X}_{i,t-1}\) into a vector. As a result, the input dimension of the first layer increases from \(D^{(0)} = P\) to \(D^{(0)} = SP\), and the corresponding weight matrices and bias vectors adjust accordingly. All remaining network architecture and notation follow the NN model in 7.1.2.

Intuitively, the LANN model treats all lagged characteristics as “static” covariates. It does not model temporal dynamics explicitly and therefore ignores potential temporal dependence structures across lags.

7.2.2 DLinear↩︎

The DLinear model first decomposes the matrix \(\boldsymbol{X}_{i,t-1}\) into a trend component \(\boldsymbol{T}_{i,t-1} \in \mathbb{R}^{S \times P}\) and a seasonal component \(\boldsymbol{S}_{i,t-1} \in \mathbb{R}^{S \times P}\). For the trend component, the \(s\)-th row of \(\boldsymbol{T}_{i,t-1}\) is defined as \[\begin{align} \label{trend95comp} \boldsymbol{T}_{i,t-1}^s = \frac{1}{|\mathcal{N}_K(s)|}\sum_{j\in \mathcal{N}_K(s)} \underline{\boldsymbol{x}}_{i,t-1}^j \in \mathbb{R}^P,\,\,\, \text{for}\,\,\, s = 1, \dots, S, \end{align}\tag{28}\] where \(\mathcal{N}_K(s) = \left\{ j \in \{1,\dots,S\} : |j-s| \le \lfloor K/2 \rfloor \right\}\) is the local index set for index \(s\), and \(\underline{\boldsymbol{x}}_{i,t-1}^j \equiv \boldsymbol{x}_{i,t-1-S+j} \in \mathbb{R}^P\) denotes the \(j\)-th row of \(\boldsymbol{X}_{i,t-1}\), which also corresponds to the \(j\)-th lagged characteristics in the sequence input. Then, for the seasonal component, we have \[\begin{align} \label{season95temp} \boldsymbol{S}_{i,t - 1} = \boldsymbol{X}_{i, t-1} - \boldsymbol{T}_{i,t-1} \in \mathbb{R}^{S \times P}. \end{align}\tag{29}\] In 28 and 29 , \(\boldsymbol{T}_{i,t-1}\) aims to capture the slow-moving and persistent temporal patterns, while \(\boldsymbol{S}_{i,t-1}\), as the residual, retains short-term fluctuations and high-frequency variations.

After decomposing \(\boldsymbol{X}_{i,t-1}\) into interpretable trend and seasonal components, the DLinear model aggregates these components and projects them into a vector of risk scores: \[\begin{align} \label{dlinear95output} \begin{aligned} \boldsymbol{z}_{i,t} &= \boldsymbol{W}_{S} \boldsymbol{S}_{i,t-1} + \boldsymbol{W}_{T} \boldsymbol{T}_{i,t-1} \in \mathbb{R}^P, \\ \widetilde{\boldsymbol{z}}_{i,t} & = \mathrm{ReLU} \left( \boldsymbol{W}_{\mathrm{hid}} \boldsymbol{z}_{i,t} + \boldsymbol{b}_{\mathrm{hid}} \right) \in \mathbb{R}^D, \\ \boldsymbol{y}_{i,t} &= g(\boldsymbol{X}_{i,t-1};\boldsymbol{\theta}) = \boldsymbol{W}_{\mathrm{out}} \widetilde{\boldsymbol{z}}_{i,t} + \boldsymbol{b}_{\mathrm{out}}, \end{aligned} \end{align}\tag{30}\] where \(\boldsymbol{W}_T\in\mathbb{R}^{1\times S}\) and \(\boldsymbol{W}_S\in\mathbb{R}^{1\times S}\) are the weight vectors that aggregate the trend and seasonal components along the temporal dimension, respectively, the weight matrices \(\boldsymbol{W}_{\mathrm{hid}} \in \mathbb{R}^{D \times P}\) and \(\boldsymbol{W}_{\mathrm{out}} \in \mathbb{R}^{2 \times D}\), together with the bias vectors \(\boldsymbol{b}_{\mathrm{hid}} \in \mathbb{R}^D\) and \(\boldsymbol{b}_{\mathrm{out}} \in \mathbb{R}^2\), define a one-hidden-layer neural network from the aggregated components to the vector of risk scores.

From 28 –30 , we can see that the DLinear model captures the dynamics of characteristics through trend–seasonal decomposition and then applies a one-hidden-layer neural network to exploit these components in forecasting.

7.2.3 LSTM↩︎

In the LSTM model, the temporal dynamics of the features are captured through a hidden state vector \(\boldsymbol{h}^s_{i,t} \in \mathbb{R}^D\) in an autoregressive manner: \[\begin{align} \begin{aligned}\label{lstm} \boldsymbol{z}_{i,t}^s &= \mathrm{tanh} (\boldsymbol{W}_{z} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{z} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_z),\\ \boldsymbol{i}_{i,t}^s &= \mathrm{sigmoid} (\boldsymbol{W}_{i} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{i} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_i),\\ \boldsymbol{f}_{i,t}^s &= \mathrm{sigmoid} (\boldsymbol{W}_{f} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{f} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_f),\\ \boldsymbol{c}_{i,t}^s &= \boldsymbol{f}_{i,t}^s \odot \boldsymbol{c}_{i,t}^{s-1} + \boldsymbol{i}_{i,t}^s \odot \boldsymbol{z}_{i,t}^s,\\ \boldsymbol{o}_{i,t}^s &= \mathrm{sigmoid} (\boldsymbol{W}_{o} \underline{\boldsymbol{x}}_{i,t-1}^s + \boldsymbol{U}_{o} \boldsymbol{h}_{i,t}^{s-1} + \boldsymbol{b}_o),\\ \boldsymbol{h}_{i,t}^s &= \boldsymbol{o}_{i,t}^s \odot \mathrm{tanh}(\boldsymbol{c}_{i,t}^s), \end{aligned} \end{align}\tag{31}\] for \(s = 1, \dots, S\), where \(\underline{\boldsymbol{x}}_{i,t-1}^s = \boldsymbol{x}_{i,\,t-1-S+s} \in \mathbb{R}^P\) is the \(s\)-th row of \(\boldsymbol{X}_{i,t-1}\), \(\boldsymbol{W}_* \in \mathbb{R}^{D \times P}\), \(\boldsymbol{U}_* \in \mathbb{R}^{D \times D}\), and \(\boldsymbol{b}_* \in \mathbb{R}^{D}\) for \(* \in \{z,i,f,o\}\) are the weight matrices and bias vectors, respectively. Here, the operator \(\odot\) denotes element-wise multiplication, while \(\tanh(\cdot)\) and \(\mathrm{sigmoid}(\cdot)\) are element-wise activation functions. Following [34], the initial states \(\boldsymbol{h}_{i,\,t-1}^{\,0}\) and \(\boldsymbol{c}_{i,\,t-1}^{\,0}\) are set to \(D\)-dimensional zero vectors.

In 31 , the hidden state \(\boldsymbol{h}_{i,t}^{s}\) serves as the key representation, summarizing the accumulated information up to step \(s\). The cell state \(\boldsymbol{c}_{i,t}^{s}\) models the internal long-term memory by carrying information across steps. Through the \(\mathrm{sigmoid}\) function, the forget gate \(\boldsymbol{f}_{i,t}^s \in (0,1)^D\), input gate \(\boldsymbol{i}_{i,t}^s \in (0,1)^D\), and output gate \(\boldsymbol{o}_{i,t}^s \in (0,1)^D\) regulate how much past information is discarded, how much new information is stored, and how much information is revealed to the hidden state \(\boldsymbol{h}_{i,t}^{s}\), respectively.

After the recursive extraction in 31 , the final hidden state \(\boldsymbol{h}_{i,t}^S\) is passed through a linear transformation to obtain the risk scores: \[\begin{align} \label{lstm95fc} \boldsymbol{y}_{i,t} = g(\boldsymbol{X}_{i,t-1}; \boldsymbol{\theta}) = \boldsymbol{W}_L \boldsymbol{h}_{i,t}^S + \boldsymbol{b}_L, \end{align}\tag{32}\] where \(\boldsymbol{W}_L \in \mathbb{R}^{2 \times D}\) and \(\boldsymbol{b}_L \in \mathbb{R}^2\) are weight matrix and bias vector, and \(\boldsymbol{\theta}\) collects all parameters within the weight matrices and bias vectors in 31 and 32 .

7.2.4 GRU↩︎

As another widely used recurrent neural network, the GRU model simplifies the LSTM model by decreasing the number of gates. Its autoregressive update for the hidden state is given in 12 . Although the ReSGA model uses the full cross-sectional input \(\mathcal{X}_{t-1}\), the GRU component itself operates independently on each asset’s temporal information \(\boldsymbol{X}_{i, t-1}\). This makes the GRU, when used alone, a purely temporal model.

Then, as in the LSTM model, the GRU model passes the final hidden state \(\boldsymbol{h}_{i,t}^S\) through a linear transformation to compute the risk scores: \[\begin{align} \label{gru95fc} \boldsymbol{y}_{i,t} = g(\boldsymbol{X}_{i,t-1}; \boldsymbol{\theta}) = \boldsymbol{W}_G \boldsymbol{h}_{i,t}^S + \boldsymbol{b}_G, \end{align}\tag{33}\] where \(\boldsymbol{W}_G \in \mathbb{R}^{2 \times D}\) and \(\boldsymbol{b}_G \in \mathbb{R}^2\) are the output weight matrix and bias vector, respectively, and \(\boldsymbol{\theta}\) includes all parameters within weight matrices and bias vectors in 12 and 33 .

7.2.5 Informer↩︎

While the vanilla Transformer model [49] has achieved remarkable success in natural language processing through the attention mechanism, several variants have been developed for different tasks, such as the Vision Transformer [50] for computer vision and the Speech Transformer [51] for audio recognition. In time-series forecasting, a milestone Transformer-type architecture is the Informer model [36], which follows the Transformer’s encoder–decoder framework but replaces the standard attention mechanism with a probabilistic sparse self-attention mechanism to handle long sequences more efficiently, making it particularly suitable for time-series modeling.

Let \(\boldsymbol{X}_{i,t - 1}^\ast = [\boldsymbol{x}_{i, t-S^\ast}, \boldsymbol{x}_{i,t-2}, \dots, \boldsymbol{x}_{i,t-1}] \in \mathbb{R}^{S^\ast \times P}\) be a truncated version of \(\boldsymbol{X}_{i,t-1}\) with \(S^\ast\) being a constant smaller than \(S\). The general architecture of Informer model is defined as: \[\begin{align} \begin{aligned} \boldsymbol{y}_{i, t} & = g(\boldsymbol{X}_{i, t-1}; \boldsymbol{\theta}) = \mathrm{InfDecoder}(\boldsymbol{X}_{i,t - 1}^\ast, \boldsymbol{H}_{i, t}; \boldsymbol{\psi}) \in \mathbb{R}^2 \\ \,\,\, \text{with}\,\,\, \boldsymbol{H}_{i, t} & = \mathrm{InfEncoder}(\boldsymbol{X}_{i, t-1}; \boldsymbol{\phi}) \in \mathbb{R}^{S \times D}, \end{aligned} \end{align}\] where \(\mathrm{InfEncoder}(\cdot; \boldsymbol{\phi}): \mathbb{R}^{S \times P} \to \mathbb{R}^{S \times D}\) represents the encoder parameterized by \(\boldsymbol{\phi}\), and \(\mathrm{InfDecoder}(\cdot, \cdot; \boldsymbol{\psi}): \mathbb{R}^{S^\ast \times P} \times \mathbb{R}^{S \times D} \to \mathbb{R}^2\) denotes the decoder that depends on the encoder output \(\boldsymbol{H}_{i,t}\) and its parameters \(\boldsymbol{\psi}\). Here, the encoder extracts the past-aware feature \(\boldsymbol{H}_{i,t}\) from the full historical input \(\boldsymbol{X}_{i,t-1}\), summarizing the underlying dynamics from the global viewpoint; afterwards, the decoder then leverages both the learned feature \(\boldsymbol{H}_{i,t}\) and the shorter input \(\boldsymbol{X}_{i,t-1}^\ast\) to predict the target vector \(\boldsymbol{y}_{i,t}\).

Specifically, the encoder first applies an embedding layer to extract embeddings from \(\boldsymbol{X}_{i,t-1}\), which are then combined as the input for subsequent layers: \[\begin{align} \label{encoder95embedding} \begin{aligned} \boldsymbol{H}_{i, t}^{(0)} & = \boldsymbol{E}^{\mathrm{val}}_{i,t-1} + \boldsymbol{E}^{\mathrm{pos}}_{i,t-1} + \boldsymbol{E}^{\mathrm{temp}}_{i,t-1} \\ & = \mathrm{VE}(\boldsymbol{X}_{i,t-1}; \boldsymbol{\zeta}) + \mathrm{PE} (\boldsymbol{X}_{i,t-1}) + \mathrm{TE}(\boldsymbol{X}_{i,t-1}) \in \mathbb{R}^{S\times D}, \end{aligned} \end{align}\tag{34}\] where \(\boldsymbol{E}^{\mathrm{val}}_{i,t-1}\), \(\boldsymbol{E}^{\mathrm{pos}}_{i,t-1}\) and \(\boldsymbol{E}^{\mathrm{temp}}_{i,t-1}\) are three \(S \times D\) matrices, referred to as the value, positional and temporal embeddings, respectively. Meanwhile, \(\mathrm{VE}(\cdot; \boldsymbol{\zeta}): \mathbb{R}^{S \times P} \rightarrow \mathbb{R}^ {S \times D}\), indexed by parameters \(\boldsymbol{\zeta}\), applies a one-dimensional temporal convolution to produce the value embedding, while \(\mathrm{PE}(\cdot): \mathbb{R}^{S \times P} \rightarrow \mathbb{R}^ {S \times D}\) and \(\mathrm{TE}(\cdot): \mathbb{R}^{S \times P} \rightarrow \mathbb{R}^ {S \times D}\) are two parameter-free mappings to provide the positional and temporal embeddings⁷.

Given \(\boldsymbol{H}_{i, t}^{(0)}\), a multi-layer self-attention block is employed to further capture complex temporal patterns: \[\begin{align} \label{encoder95attn} \boldsymbol{H}_{i,t}^{(l)} = \mathrm{FFN}\big( \mathrm{SelfAttn}\big( \boldsymbol{H}_{i,t}^{(l-1)}; \boldsymbol{\omega}^{(l)} \big) ; \boldsymbol{\delta}^{(l)} \big) \in \mathbb{R}^{S \times D}, \,\,\, \text{for}\,\,\, l = 1, \dots, L, \end{align}\tag{35}\] where \(L\) is the number of self-attention layers, and \(\boldsymbol{H}_{i,t}^{(l)}\) denotes the output from the \(l\)-th layer. Here, two operators are involved: \(\mathrm{SelfAttn}\) and \(\mathrm{FFN}\). Particularly, \(\mathrm{SelfAttn}(\cdot; \boldsymbol{\omega}): \mathbb{R}^{S \times D} \to \mathbb{R}^{S \times D}\) is the self-attention operator based on the probabilistic sparse technique⁸, indexed by parameters \(\boldsymbol{\omega}^{(l)}\); meanwhile, \(\mathrm{FFN}(\cdot;\boldsymbol{\delta}): \mathbb{R}^{S \times D} \to \mathbb{R}^{S \times D}\) denotes a row-wise feed-forward network, that is, a one-hidden-layer fully connected network applied independently to each row of the input⁹, parameterized by \(\boldsymbol{\delta}\).

After the recursive extraction, the output of encoder is defined as \[\begin{align} \boldsymbol{H}_{i,t} \equiv \boldsymbol{H}_{i,t}^{(L)} = \mathrm{InfEncoder}(\boldsymbol{X}_{i, t-1}; \boldsymbol{\phi}), \end{align}\] where \(\boldsymbol{\phi}\) encompass all parameters in 34 –35 . The matrix \(\boldsymbol{H}_{i,t}\) is referred to as the past-aware feature in our paper, as it summarizes the general information contained in \(\boldsymbol{X}_{i,t-1}\).

So far, we have illustrated the architecture of encoder. Compared to the encoder, the decoder also first applies an embedding layer to map \(\boldsymbol{X}_{i,t-1}^\ast\) into an aggregated feature: \[\begin{align} \label{decoder95embedding} \boldsymbol{F}_{i, t}^{(0)} = \mathrm{VE}(\boldsymbol{X}_{i,t-1}^\ast; \boldsymbol{\zeta}^\ast) + \mathrm{PE} (\boldsymbol{X}_{i,t-1}^\ast) + \mathrm{TE}(\boldsymbol{X}_{i,t-1}^\ast) \in \mathbb{R}^{S^\ast \times D}, \end{align}\tag{36}\] where \(\mathrm{VE}\), \(\mathrm{PE}\), and \(\mathrm{TE}\) operators are consistent with those in 34 , \(\boldsymbol{\zeta}^\ast\) denotes the parameters of the decoder’s embedding layer, and \(\boldsymbol{F}_{i,t}^{(0)}\) is the resulting aggregated feature.

Next, a multi-layer cross-attention block is applied to compute the hidden states within the decoder: \[\begin{align} \label{decoder95attn} \begin{aligned} \boldsymbol{F}_{i,t}^{(l)} &= \mathrm{FFN} \big( \mathrm{CrossAttn} \big( \overline{\boldsymbol{F}}_{i,t}^{(l)}, \boldsymbol{H}_{i,t} ; \boldsymbol{\gamma}^{(l)} \big) ; \boldsymbol{\delta}^{(l + L)} \big) \in \mathbb{R}^{S^\ast \times D} \\ \text{with}\,\,\, \overline{\boldsymbol{F}}_{i,t}^{(l)} &= \mathrm{CausalAttn} \big( \boldsymbol{F}_{i,t}^{(l-1)}; \boldsymbol{\omega}^{(l + L)} \big) \in \mathbb{R}^{S^\ast \times D}, \,\,\, l = 1, \dots, L, \end{aligned} \end{align}\tag{37}\] where \(\boldsymbol{F}_{i,t}^{(l)}\) denotes the output from the \(l\)-th layer. Compared to the self-attention block in 35 for the encoder, the cross-attention block has the same depth and applies the same \(\mathrm{FFN}\) operator, but introduces two additional operators: \(\mathrm{CausalAttn}\) and \(\mathrm{CrossAttn}\).

Specifically, \(\mathrm{CausalAttn}(\cdot; \boldsymbol{\omega}): \mathbb{R}^{S^\ast \times D} \to \mathbb{R}^{S^\ast \times D}\) is the causal attention operator, which is identical to \(\mathrm{SelfAttn}\) except that it imposes a no-forward-looking constraint, forcing the decoder to follow a strictly autoregressive structure¹⁰. Meanwhile, \(\mathrm{CrossAttn}(\cdot, \cdot; \boldsymbol{\gamma}): \mathbb{R}^{S^\ast \times D} \times \mathbb{R}^{S \times D} \to \mathbb{R}^{S^\ast \times D}\) denotes the cross-attention operator, which functions similarly to self-attention but also attends to the past-aware feature \(\boldsymbol{H}_{i,t}\) (i.e., the encoder’s output) rather than only its own input. This mechanism allows the decoder to selectively integrate information from the full historical context, aligning its input with the temporal dependencies captured by the encoder.

Lastly, to obtain the risk scores, the most up-to-date information extracted by the decoder is fed into a linear transformation: \[\begin{align} \label{informer95fc} \boldsymbol{y}_{i,t} = \mathrm{InfDecoder}(\boldsymbol{X}_{i,t - 1}^\ast, \boldsymbol{H}_{i, t}; \boldsymbol{\psi}) = \boldsymbol{W}_{\mathrm{in}} \boldsymbol{f}_{i,t} + \boldsymbol{b}_{\mathrm{in}}, \end{align}\tag{38}\] where \(\boldsymbol{f}_{i,t} \in \mathbb{R}^{D}\) denotes the last row of \(\boldsymbol{F}_{i,t}^{(L)}\), representing the latest information available for time \(t\), \(\boldsymbol{W}_{in} \in \mathbb{R}^{2 \times D}\) and \(\boldsymbol{b}_{in} \in \mathbb{R}^2\) are the weight matrix and bias vector, respectively, and \(\boldsymbol{\psi}\) includes all parameters involved in 36 –37 . Clearly, the Informer model can be regarded as a special case of the general temporal model, with \(\boldsymbol{\theta}\) encompassing all parameters within \(\boldsymbol{\phi}\) and \(\boldsymbol{\psi}\).

7.2.6 EInformer and DInformer↩︎

As the vanilla Informer follows the encoder–decoder framework, recent studies have highlighted the potential of encoder-only and decoder-only variants in time-series forecasting [52]–[54]. Motivated by these works, we implement two variants, the encoder-only Informer (EInformer) and the decoder-only Informer (DInformer), for VaR and ES forecasting.

Specifically, the EInformer model consists solely of the encoder but modifies the original input \(\boldsymbol{X}_{i,t-1}\) by appending an additional row of learnable parameters at the end: \[\begin{align} \overline{\boldsymbol{X}}_{i,t-1} = [\boldsymbol{X}_{i,t-1}, \boldsymbol{\nu}] \in \mathbb{R}^{(S + 1) \times P}, \end{align}\] where \(\overline{\boldsymbol{X}}_{i,t-1}\) is the augmented input, and \(\boldsymbol{\nu} \in \mathbb{R}^P\) is a vector of learnable parameters. This additional vector \(\boldsymbol{\nu}\), referred to as the predictive token, serves as a placeholder representing the characteristics at next (yet unseen) time step. It allows the model to generate a forecast without explicitly relying on a decoder. During training, the predictive token interacts with the historical information through the self-attention mechanism, enabling the encoder to learn how the future prediction depends on the past characteristic. In this sense, the predictive token acts as a query that asks, “given the past, what comes next”?

Next, the augmented input \(\overline{\boldsymbol{X}}_{i,t-1}\) is passed through the encoder as in 34 –35 , and the last row of the encoder output, corresponding to the predictive token, is used to compute the risk scores: \[\begin{align} \label{EInformer} \boldsymbol{y}_{i,t} = g(\boldsymbol{X}_{i,t-1}; \boldsymbol{\theta}) = \boldsymbol{W}_e \overline{\boldsymbol{h}}_{i,t} + \boldsymbol{b}_e \in \mathbb{R}^2, \end{align}\tag{39}\] where \(\boldsymbol{W}_e \in \mathbb{R}^{2 \times D}\) and \(\boldsymbol{b}_e \in \mathbb{R}^2\) are the output weight matrix and bias vector, respectively. Here, \(\overline{\boldsymbol{h}}_{i,t} \in \mathbb{R}^{D}\) denotes the last row of the encoder output \[\begin{align} \overline{\boldsymbol{H}}_{i,t} = [\overline{\boldsymbol{h}}_{i,t-S}, \dots, \overline{\boldsymbol{h}}_{i,t}] = \mathrm{InfEncoder}(\overline{\boldsymbol{X}}_{i,t-1}; \boldsymbol{\phi}) \in \mathbb{R}^{(S + 1) \times D}, \end{align}\] and \(\boldsymbol{\theta}\) includes all parameters in \(\boldsymbol{\nu}\), \(\boldsymbol{W}_e\), \(\boldsymbol{b}_e\), and \(\boldsymbol{\phi}\).

In contrast, the DInformer model relies solely on the decoder and is defined as \[\begin{align} \label{DInformer} \boldsymbol{y}_{i,t} = g(\boldsymbol{X}_{i,t-1}; \boldsymbol{\theta}) = \mathrm{InfDecoder}(\boldsymbol{X}_{i,t-1}; \boldsymbol{\psi}) \in \mathbb{R}^2. \end{align}\tag{40}\] Compared with the decoder in the vanilla Informer model, the DInformer uses \(\boldsymbol{X}_{i,t-1}\) rather than \(\boldsymbol{X}_{i,t-1}^\ast\) as input. Moreover, since the past-aware representation \(\boldsymbol{H}_{i,t}\) from the encoder is absent, the cross-attention operator in 37 is omitted.

7.3 Spatial-Temporal Models↩︎

Besides our proposed ReSGA, we consider another spatial-temporal model in (9 ): SGA. The SGA model can be viewed as a simplified version of the ReSGA model. It has exactly the same encoder as that in the ReSGA model, but omits the retriever component and hence removes the retrieval layer in the decoder of ReSGA. All remaining components, especially the mechanism for learning and utilizing dynamic group structures, are retained. Therefore, comparing ReSGA with SGA isolates the role of the retrieval mechanism and highlights the additional value of incorporating long-term cross-asset information beyond only considering contemporaneous group effects.

7.4 Econometric Models↩︎

As econometric benchmarks, we consider two univariate models that are estimated separately for each asset using only its own historical returns. Let \(r_{i,t}\) denote the excess return of asset \(i\) at time \(t\), and let \(v_{i,t}\) and \(e_{i,t}\) denote the corresponding conditional \(\mathrm{VaR}\) and \(\mathrm{ES}\) at level \(\tau\).

7.4.1 GAS↩︎

Following [11], we adopt the one-factor GAS model by assuming that both \(\mathrm{VaR}\) and \(\mathrm{ES}\) are driven by a single latent factor \(\kappa_{i,t}\): \[\begin{align} v_{i,t} &= a_i \exp(\kappa_{i,t}), \qquad e_{i,t} = b_i \exp(\kappa_{i,t}), \qquad b_i < a_i < 0, \end{align}\] with dynamics \[\begin{align} \kappa_{i,t} = \omega_i + \beta_i \kappa_{i,t-1} + \frac{\gamma_i}{b_i \exp(\kappa_{i,t-1})} \left[ \frac{1}{\tau}\mathbf{1}\!\left\{r_{i,t-1} \le a_i \exp(\kappa_{i,t-1})\right\} r_{i,t-1} - b_i \exp(\kappa_{i,t-1}) \right]. \end{align}\] where \((\omega_i, a_i, b_i, \beta_i, \gamma_i)\) is a parameter vector. We set \(\omega_i = 0\) as in [11] for identification, and estimate the remaining parameters by minimizing the FZ loss.

The one-factor GAS specification assumes that the time variation in \(\mathrm{VaR}\) and \(\mathrm{ES}\) is governed by a common latent factor \(\kappa_{i,t}\), so that the two risk measures move proportionally over time. The recursion for \(\kappa_{i,t}\) is score-driven: Under the one-factor parameterization, the updating term is obtained from the derivative of the FZ loss with respect to \(\kappa_{i,t}\). The parameter \(\beta_i\) controls the persistence of the latent factor, while the parameter \(\gamma_i\) determines the sensitivity of tail risk forecasts to new return realizations. The parameters \(a_i\) and \(b_i\) scale the latent factor into the corresponding \(\mathrm{VaR}\) and \(\mathrm{ES}\) forecasts.

7.4.2 GARCH↩︎

[11] also consider a GARCH(\(1, 1\)) specification estimated by FZ loss minimization. This model is defined as \[\begin{align} r_{i,t} &= \sigma_{i,t}\eta_{i,t}, \qquad \eta_{i,t} \sim \text{i.i.d. } F_{\eta,i}(0,1), \\ \sigma_{i,t}^2 &= \omega_i + \beta_i \sigma_{i,t-1}^2 + \gamma_i r_{i,t-1}^2. \end{align}\] Under this location-scale structure, \[\begin{align} v_{i,t} &= a_i \sigma_{i,t}, \qquad a_i = F_{\eta,i}^{-1}(\tau), \\ e_{i,t} &= b_i \sigma_{i,t}, \qquad b_i = E(\eta_{i,t} \mid \eta_{i,t} \le a_i), \end{align}\] where \((\omega_i, a_i, b_i, \beta_i, \gamma_i)\) is a parameter vector, \(\sigma_{i,t}\) is the conditional volatility of asset \(i\) at time \(t\), and \(\eta_{i,t}\) is an i.i.d. innovation with distribution \(F_{\eta,i}(0,1)\). As in [11], we set \(\omega_i = 1\).

Under this GARCH specification, both \(\mathrm{VaR}\) and \(\mathrm{ES}\) are proportional to the conditional volatility \(\sigma_{i,t}\). Therefore, their time variation is fully determined by the GARCH recursion for \(\sigma_{i,t}\). The parameter \(\beta_i\) captures the persistence of conditional volatility, whereas the parameter \(\gamma_i\) measures the impact of lagged squared returns on current risk. The parameters \(a_i\) and \(b_i\) correspond to the \(\tau\)-quantile and the associated lower-tail conditional mean of the standardized innovation \(\eta_{i,t}\), respectively. In [11], this model is estimated by minimizing the FZ loss, so that the GARCH dynamics are fitted to the joint forecasting problem for \(\mathrm{VaR}\) and \(\mathrm{ES}\) rather than to the full conditional distribution of returns.

8 Details of Operators in the Informer Model↩︎

8.1 Embedding Layer↩︎

As described in 34 and 36 , the inputs to the encoder and decoder, \(\boldsymbol{X}_{i,t-1}\) and \(\boldsymbol{X}_{i,t-1}^\ast\), are first processed by the embedding layer that consists of three embedding operators, \(\mathrm{VE}\), \(\mathrm{PE}\), and \(\mathrm{TE}\), which generate the value, positional, and temporal embeddings, respectively. For notational simplicity, we assume a generic input matrix \(\boldsymbol{X}\in \mathbb{R}^{S \times P}\), which contains \(P\) different characteristics observed over \(S\) consecutive time steps. Based on this formulation, the details of each embedding operator are described below.

8.1.1 Value Embedding Operator↩︎

The value embedding operator \(\mathrm{VE}\) applies a simple one-dimensional convolution to the input matrix \(\boldsymbol{X}\), mapping it into an output matrix \(\boldsymbol{E}^{\mathrm{val}} \in \mathbb{R}^{S \times D}\). The convolution employs \(D\) learnable filters (also referred to as kernels), where \(D\) corresponds to the dimensionality of the hidden state in the Informer model. Each filter is a smaller matrix of size \(W \times P\) that scans the input along the row axis. Specifically, each filter slides over the input, computes the dot product between its weights and the corresponding input sub-vectors, and aggregates the results to produce the convolution output.

The resulting matrix \(\boldsymbol{E}^{\mathrm{val}}\) effectively forms a weighted average of the input among nearby \(W\) time points with data-driven weights and bias. The learnable parameters of these filters are collected in \(\boldsymbol{\zeta}\). Figure 5 provides a visual illustration of this one-dimensional convolution process.

8.1.2 Positional Embedding Operator↩︎

The positional embedding operator \(\mathrm{PE}\) is an essential component in Transformer-type models, as the attention mechanism itself treats the inputs as an unordered set of numbers rather than an ordered matrix. In other words, the attention mechanism does not inherently recognize that one element corresponds to a specific position or another position, unless positional information is explicitly injected¹¹.

Following [49], a fixed sinusoidal mapping is adopted to construct the positional embedding representation, which encodes the positional information of each element of input \(\boldsymbol{X}\): \[\boldsymbol{E}^{\mathrm{pos}} = \mathrm{PE}(\boldsymbol{X}) \in \mathbb{R}^{S \times D},\] where the \((i,j)\)-th element of \(\boldsymbol{E}^{\mathrm{pos}}\) is defined as \[\begin{align} e_{i,j}^{\mathrm{pos}} = \begin{cases} \sin \!\left(\dfrac{i - 1}{10000^{\frac{j - 1}{2D}}} \right), & \text{if } j \text{ is odd}, \\[6pt] \cos \!\left(\dfrac{i - 1}{10000^{\frac{j - 2}{2D}}} \right), & \text{if } j \text{ is even}. \end{cases} \end{align}\] By definition, the operator \(\mathrm{PE}\) contains no learnable parameters; it depends solely on the shape of \(\boldsymbol{X}\) and provides deterministic positional representation that helps the attention mechanism preserve the sequential structure of the input.

8.1.3 Temporal Embedding Operator↩︎

Motivated by the positional embedding operator, [36] introduce the temporal embedding operator \(\mathrm{TE}\), which transforms absolute time stamps into representations that capture potential seasonality and tendency. To construct the temporal embedding, we associate each row of \(\boldsymbol{X}\) with a time stamp consisting of its month, day, and corresponding weekday. Each of these calendar components can be represented by an integer: For example, months are encoded from \(0\) (January) to \(11\) (December), days from \(0\) (1st) to \(30\) (31st), and weekdays from \(0\) (Monday) to \(6\) (Sunday). Consequently, each time stamp is represented as a three-dimensional integer vector encoding its month, day, and weekday information.

Let \(\boldsymbol{M}\in \mathbb{R}^{S \times 3}\) denote the matrix containing these integer-encoded time stamps for all \(S\) observations. Given \(\boldsymbol{M}\), the temporal embedding representation is defined as \[\begin{align} \boldsymbol{E}^{\mathrm{temp}} = \mathrm{TE}(\boldsymbol{X}) \in \mathbb{R}^{S \times D}, \end{align}\] where the \((i,j)\)-th element of \(\boldsymbol{E}^{\mathrm{temp}}\) is given by \[\begin{align} e_{i,j}^{\mathrm{temp}} = \begin{cases} \displaystyle \sum_{k = 1}^3 \sin\!\left(10m_{i,k} \exp\!\left(\frac{4(j - 1)}{D}\right) \right), & \text{if } j \text{ is odd}, \\[6pt] \displaystyle \sum_{k = 1}^3 \cos\!\left(10m_{i,k} \exp\!\left(\frac{4(j - 2)}{D} \right) \right), & \text{if } j \text{ is even}, \end{cases} \end{align}\] where \(m_{i,k}\) denotes the \((i,k)\)-th element of \(\boldsymbol{M}\).

Similar to \(\mathrm{PE}\), \(\mathrm{TE}\) contains no learnable parameters. It deterministically maps absolute time stamps into continuous temporal representations, allowing the attention mechanism to better capture seasonal patterns and potential trends inherent in the data.

8.2 Probabilistic Sparse Attention↩︎

As discussed in Section 7.2.5, the Informer model employs three types of attention operators: \(\mathrm{SelfAttn}\), \(\mathrm{CausalAttn}\), and \(\mathrm{CrossAttn}\). To efficiently handle long time input, these operators are further equipped with a probabilistic sparse (ProbSparse) technique, resulting in the so-called probabilistic sparse attention [36]. The details of each attention operator are introduced below, while we assume a generic matrix input \(\boldsymbol{X}\in \mathbb{R}^{S \times D}\) for these attention operators for notational simplicity.

8.2.1 Self-Attention Operator↩︎

First, the input \(\boldsymbol{X}\) is linearly projected into three matrices: the query \(\boldsymbol{Q}\), the key \(\boldsymbol{K}\), and the value \(\boldsymbol{V}\), defined as \[\begin{align} \label{QKV} \boldsymbol{Q}= \boldsymbol{X}\boldsymbol{W}_Q \in \mathbb{R}^{S \times D}, \quad \boldsymbol{K}= \boldsymbol{X}\boldsymbol{W}_K \in \mathbb{R}^{S \times D}, \quad \boldsymbol{V}= \boldsymbol{X}\boldsymbol{W}_V \in \mathbb{R}^{S \times D}, \end{align}\tag{41}\] where \(\boldsymbol{W}_Q, \boldsymbol{W}_K, \boldsymbol{W}_V \in \mathbb{R}^{D \times D}\) are learnable weight matrices. In the multi-head attention setting with \(M\) heads, these matrices are further projected into head-specific low-dimensional subspaces: \[\begin{align} \label{head95QKV} \boldsymbol{Q}^{(m)} = \boldsymbol{Q}\boldsymbol{W}_Q^{(m)} \in \mathbb{R}^{S \times \underline{D}}, \quad \boldsymbol{K}^{(m)} = \boldsymbol{K}\boldsymbol{W}_K^{(m)} \in \mathbb{R}^{S \times \underline{D}}, \quad \boldsymbol{V}^{(m)} = \boldsymbol{V}\boldsymbol{W}_V^{(m)} \in \mathbb{R}^{S \times \underline{D}}, \end{align}\tag{42}\] for \(m = 1, \dots, M\), where \(\boldsymbol{W}_Q^{(m)}, \boldsymbol{W}_K^{(m)}, \boldsymbol{W}_V^{(m)} \in \mathbb{R}^{D \times \underline{D}}\) are the head-specific weight matrices, and \(\underline{D} = D / M\) denotes the dimensionality of each head-specific subspace.

Next, the ProbSparse technique is then applied independently within each head. Different from the standard Transformer attention where each query attends to all keys, the ProbSparse technique reduces complexity by sampling a subset of keys for each query and computing a sparsity score to identify the most informative queries.

Specifically, for \(i = 1, \dots, S\), let \(\boldsymbol{q}_i^{(m)}\), \(\boldsymbol{k}_i^{(m)}\), and \(\boldsymbol{v}_i^{(m)}\) denote the \(i\)-th rows of \(\boldsymbol{Q}^{(m)}\), \(\boldsymbol{K}^{(m)}\), and \(\boldsymbol{V}^{(m)}\), respectively. Moreover, let \(\mathcal{S}_i\) be a uniformly random sampled subset of \(\{1, \dots, S\}\) with \(|\mathcal{S}_i| = \lfloor c \log S \rfloor\), where \(c\) is a hyperparameter that balances computational efficiency and approximation accuracy. The sparsity score is then defined as \[\begin{align} \label{attn95sparse95score} \gamma_i^{(m)} = \max_{j \in \mathcal{S}_i} \beta^{(m)}_{i,j} - \frac{1}{|\mathcal{S}_i|} \sum_{j \in \mathcal{S}_i} \beta^{(m)}_{i,j}, \end{align}\tag{43}\] where \(\beta^{(m)}_{i,j}\), referred to as the attention coefficient, measures the similarity between the query vector \(\boldsymbol{q}_i^{(m)}\) and the key vector \(\boldsymbol{k}_j^{(m)}\) through their scaled inner product: \[\begin{align} \beta^{(m)}_{i,j} = \frac{\boldsymbol{q}_i^{(m)\prime} \boldsymbol{k}_j^{(m)}}{\sqrt{\underline{D}}}. \end{align}\]

Intuitively, the sparsity score \(\gamma_i^{(m)}\) quantifies how much the maximum of similarity between the query vector \(\boldsymbol{q}_i^{(m)}\) and its sampled key vectors \(\{\boldsymbol{k}_j^{(m)} : j \in \mathcal{S}_i\}\) deviates from its mean. A larger sparsity score indicates that the query vector has a few highly similar key vectors, implying that these key–value pairs are more informative and thus require full attention computation. Conversely, a smaller sparsity score suggests that the query and key vectors exhibit more uniform similarity, allowing the corresponding attention computation to be pruned to reduce complexity.

To implement this idea, let \(\mathcal{I}^{(m)}\) denote the set of indices corresponding to the largest \(\lfloor c \log S \rfloor\) values among \(\{\gamma_i^{(m)} : i = 1, \dots, S\}\). The output of the \(m\)-th attention head is then defined as \(\boldsymbol{C}^{(m)} \in \mathbb{R}^{S \times \underline{D}}\) with its \(i\)-th row satisfying \[\begin{align} \label{attn95head} \boldsymbol{c}_{i}^{(m)} = \begin{cases} \displaystyle \sum_{j=1}^S \alpha_{i,j}^{(m)} \boldsymbol{v}_j^{(m)}, & \text{if } i \in \mathcal{I}^{(m)}, \\[8pt] \displaystyle \frac{1}{S} \sum_{j=1}^S \boldsymbol{v}_j^{(m)}, & \text{otherwise}, \end{cases} \end{align}\tag{44}\] where \(\alpha_{i,j}^{(m)}\) denotes the normalized attention weight, defined as \[\begin{align} \alpha_{i,j}^{(m)} = \mathrm{SoftMax}\left( \beta_{i,j}^{(m)} \right) = \frac{\exp(\beta_{i,j}^{(m)})}{\sum_{k = 1}^S \exp\left( \beta_{i,k}^{(m)} \right)}. \end{align}\]

Notably, by employing the sparsity score defined in 43 , only \(\lfloor c \log S \rfloor\) head-specific outputs \(\boldsymbol{c}_{i}^{(m)}\) require full attention computation involving \(\{\alpha_{i,j}^{(m)} : j = 1, \dots, S\}\), while the remaining ones are approximated by the average of all value vectors. Consequently, the computational complexity is reduced from \(\mathcal{O}(S^2)\) in the standard self-attention mechanism to \(\mathcal{O}(S \log S)\) in the ProbSparse self-attention, providing a significant efficiency advantage for handling long input sequences (i.e., large value of \(S\)) in the Informer model.

Finally, the outputs from all attention heads are concatenated and linearly projected to obtain the final self-attention output \(\boldsymbol{Y}\): \[\begin{align} \label{attn95fc} \boldsymbol{Y}= \mathrm{SelfAttn}(\boldsymbol{X}; \boldsymbol{\omega}) = \boldsymbol{C}\boldsymbol{W}_O \in \mathbb{R}^{S \times D}, \,\,\, \text{with} \,\,\, \boldsymbol{C}= \left[ \boldsymbol{C}^{(1)}, \boldsymbol{C}^{(2)}, \dots, \boldsymbol{C}^{(M)} \right] \in \mathbb{R}^{S \times D}, \end{align}\tag{45}\] where \(\boldsymbol{W}_O \in \mathbb{R}^{D \times D}\) is a learnable output projection matrix. The parameter set \(\boldsymbol{\omega}\) includes all learnable weights involved in 41 , 42 , and 45 .

8.2.2 Causal Attention Operator↩︎

The causal-attention operator \(\mathrm{CausalAttn}\) extends the self-attention operator by imposing a no-forward-looking constraint. This restriction enforces a strictly autoregressive structure within the attention computation. Specifically, following the same procedure in 41 –43 , we obtain \(\big\{ \big( \boldsymbol{q}_i^{(m)}, \boldsymbol{k}_i^{(m)}, \boldsymbol{v}_i^{(m)} \big) \big\}\) and compute the sparsity score \(\overline{\gamma}_i^{(m)}\) in \(\mathrm{CausalAttn}\) as \[\begin{align} \label{causal95sparse95score} \overline{\gamma}_i^{(m)} = \max_{j \in \mathcal{S}_i} \beta^{(m)}_{i,j} - \frac{1}{|\mathcal{S}_i|} \sum_{j \in \mathcal{S}_i} \beta^{(m)}_{i,j}, \end{align}\tag{46}\] where the attention coefficient \(\beta_{i,j}\) is defined as before, but the sampling set \(\mathcal{S}_i\) is drawn from \(\{1, \dots, i\}\) instead of \(\{1, \dots, S\}\), with size \(\lfloor c \log S \rfloor\).

Similarly, let \(\overline{\mathcal{I}}^{(m)}\) denote the indices of the largest \(\lfloor c \log S \rfloor\) values among \(\{\overline{\gamma}_i^{(m)} : i = 1, \dots, S\}\). The head-specific output \(\overline{\boldsymbol{C}}^{(m)}\) is then obtained row-wisely as \[\begin{align} \label{causal95head} \overline{\boldsymbol{c}}_{i}^{(m)} = \begin{cases} \displaystyle \sum_{j=1}^i \overline{\alpha}_{i,j}^{(m)} \boldsymbol{v}_j^{(m)}, & \text{if } i \in \overline{\mathcal{I}}^{(m)}, \\[8pt] \displaystyle \frac{1}{i} \sum_{j=1}^i \boldsymbol{v}_j^{(m)}, & \text{otherwise}, \end{cases} \end{align}\tag{47}\] where the normalized attention weight \(\overline{\alpha}_{i,j}^{(m)}\) is given by \[\begin{align} \overline{\alpha}_{i,j}^{(m)} = \mathrm{SoftMax}\left(\beta_{i,j}^{(m)} \right) = \frac{\exp(\beta_{i,j}^{(m)})}{\sum_{k = 1}^i \exp(\beta_{i,k}^{(m)})}. \end{align}\]

The outputs from all attention heads are concatenated and linearly projected to produce the causal-attention output \(\boldsymbol{Y}\): \[\begin{align} \label{causal95fc} \begin{aligned} \boldsymbol{Y}&= \mathrm{CausalAttn}(\boldsymbol{X}; \boldsymbol{\omega}) = \overline{\boldsymbol{C}} \boldsymbol{W}_O \in \mathbb{R}^{S \times D}, \\[3pt] \text{with} \quad \overline{\boldsymbol{C}} &= \left[ \overline{\boldsymbol{C}}^{(1)}, \overline{\boldsymbol{C}}^{(2)}, \dots, \overline{\boldsymbol{C}}^{(M)} \right] \in \mathbb{R}^{S \times D}, \end{aligned} \end{align}\tag{48}\] where \(\boldsymbol{W}_O\) is defined in 45 , and \(\boldsymbol{\omega}\) includes all learnable weights involved in 41 , 42 , and 48 . The causal constraint guarantees that the model respects temporal order, making it suitable for autoregressive sequence modeling.

8.2.3 Cross-Attention Operator↩︎

The cross-attention operator \(\mathrm{CrossAttn}\) is similar to the self-attention operator, but it requires two generic inputs, \(\boldsymbol{X}\in \mathbb{R}^{S \times D}\) and \(\boldsymbol{Z}\in \mathbb{R}^{S^\ast \times D}\), where \(S^\ast\) is another constant and can be different from \(S\). The key distinction lies in that the key and value matrices, \(\boldsymbol{K}\) and \(\boldsymbol{V}\), are computed from the external input \(\boldsymbol{Z}\) rather than from \(\boldsymbol{X}\). This design enables the model to align and integrate information between two different representations (e.g., between the encoder’s output \(\boldsymbol{H}_{i,t}\) and the decoder’s hidden state \(\boldsymbol{F}_{i,t}\) in the Informer architecture).

Specifically, the query \(\boldsymbol{Q}\), key \(\boldsymbol{K}\), and value \(\boldsymbol{V}\) for \(\mathrm{CrossAttn}\) are defined as \[\begin{align} \boldsymbol{Q}= \boldsymbol{X}\boldsymbol{W}_Q \in \mathbb{R}^{S \times D}, \,\,\, \boldsymbol{K}= \boldsymbol{Z}\boldsymbol{W}_K \in \mathbb{R}^{S^\ast \times D}, \,\,\, \text{and} \,\,\, \boldsymbol{V}= \boldsymbol{Z}\boldsymbol{W}_V \in \mathbb{R}^{S^\ast \times D}. \end{align}\] Given \(\boldsymbol{Q}\), \(\boldsymbol{K}\), and \(\boldsymbol{V}\), following the same computation steps as in 42 –45 , the output of the cross-attention operator is obtained as \[\begin{align} \boldsymbol{Y}= \mathrm{CrossAttn}(\boldsymbol{X}, \boldsymbol{Z}) = \boldsymbol{C}\boldsymbol{W}_O \in \mathbb{R}^{S \times D}. \end{align}\] Through this mechanism, the cross-attention operator allows the model to selectively attend to the most relevant parts of \(\boldsymbol{Z}\) while processing \(\boldsymbol{X}\), effectively establishing an information bridge between different representation spaces.

9 Implementation Details↩︎

9.1 Additional Missing-Value Preprocessing↩︎

As described in 3.1, firm characteristics are rank-normalized cross-sectionally each month to reduce the influence of outliers, and missing characteristic values are imputed using contemporaneous cross-sectional medians. For the models under the unified learning framework, we drop stock-month observations when their related \(r_{i,t}\) are missing. For the econometric models, GAS and GARCH, which are estimated separately for each asset using return histories, we remove missing \(r_{i,t}\) observations and concatenate the remaining non-missing observations into a valid return series.

9.2 Regularization Learning↩︎

As discussed in 2.3, we apply a range of deep learning models to VaR and ES forecasting. While the deep learning models offer strong capacity and flexibility in extracting nonlinear relationship from data, they are also prone to overfitting. In addition, their complex architectures and high-dimensional parameter spaces may lead to unstable estimates due to sensitivity to initialization and convergence to local optima. To address the above issues, we adopt two regularization techniques for these models.

First, a standard approach to mitigate overfitting is to partition the data into three disjoint parts: training, validation, and testing samples, while preserving the temporal ordering. The training sample is used for parameter estimation, the validation sample for hyperparameter tuning and early stopping, and the testing sample exclusively for evaluating out-of-sample performance. In this study, we apply this split within each expanding window in chronological order. Details of the window construction are provided in 3.1.

Second, following [32], we adopt an ensemble learning strategy for regularization. Model parameters \(\boldsymbol{\theta}\) are estimated by solving the optimization problem in 5 using the Adam algorithm [31]. Due to the nonlinearity of the models and the sensitivity of Adam to initial values, the resulting parameter estimates may vary across different initializations. To mitigate the impact of initialization, we repeat the optimization from multiple random initializations. For each run, we estimate the parameters and obtain forecasts \(\widehat{\mathrm{VaR}}_{i,t}\) and \(\widehat{\mathrm{ES}}_{i,t}\). The final forecasts are then computed as the average across these runs (five in our implementation). This ensemble procedure reduces stochasticity caused by random initialization and improves the stability of the predictions.

9.3 Hyperparameters↩︎

As with any deep learning model, performance depends on the choice of hyperparameters, for which optimal selection is generally infeasible in theory. We therefore adopt standard settings commonly used in the literature for the learning rate, batch size, and other non-structural hyperparameters. Our main focus is on structural hyperparameters that govern model capacity, primarily the hidden-state dimension \(D\), the depth \(L\) (number of layers), and, for the Informer family, the number of attention heads \(M\). To reduce computational costs, hyperparameters are tuned only in the first rolling window, and based on performance on the validation sample, ensuring that no look-ahead bias is introduced. Then, they are held fixed across subsequent windows, reflecting a realistic forecasting setting in which frequent re-tuning is impractical and time-consuming.

Moreover, these hyperparameters are also varied to study the virtue of model complexity discussed in 3.4. See 15 for a summary of the corresponding hyperparameter settings for each parameter-scale tier.

10 Statistical Testing↩︎

In this appendix, we introduce the details of four tests used in 3.

10.1 DM Test↩︎

The DM test evaluates whether two competing models exhibit statistically different predictive accuracy under a chosen loss function. Suppose we compare two forecasting models, denoted \(a\) and \(b\), for their ability to predict \(\mathrm{VaR}_{i,t}\) and \(\mathrm{ES}_{i,t}\). Let \(\left (\widehat{\mathrm{VaR}}^{(a)}_{i,t}, \widehat{\mathrm{ES}}^{(a)}_{i,t} \right)\) and \(\left(\widehat{\mathrm{VaR}}^{(b)}_{i,t}, \widehat{\mathrm{ES}}^{(b)}_{i,t} \right)\) be the forecasts of \(\left (\mathrm{VaR}_{i,t}, \mathrm{ES}_{i,t} \right)\) given by models \(a\) and \(b\), respectively. The loss differential between models \(a\) and \(b\) at time \(t\) is defined as \[\begin{align} & d_{t}(a, b) = \frac{1}{N_t} \sum_{i = 1}^{N_t} d_{i,t}(a, b), \\ \,\,\, \text{with} \,\,\, & d_{i,t}(a, b)= \ell_{\mathrm{FZ0}} \left(r_{i,t}, \left( \widehat{\mathrm{VaR}}^{(a)}_{i,t}, \widehat{\mathrm{ES}}^{(a)}_{i,t} \right) \right) - \ell_{\mathrm{FZ0}} \left(r_{i,t}, \left( \widehat{\mathrm{VaR}}^{(b)}_{i,t}, \widehat{\mathrm{ES}}^{(b)}_{i,t} \right) \right), \end{align}\] where \(\ell_{\mathrm{FZ0}}(\cdot)\) is the FZ loss function defined in 2 .

Since the loss differential measures the difference of predictive powers of models \(a\) and \(b\), the DM test can examine the following null hypothesis \(\mathcal{H}_{0}\) to check whether two competing models exhibit statistically different predictive accuracy: \[\begin{align} \mathcal{H}_{0}: \mathbb{E}\left[d_t(a, b) \right] = 0. \end{align}\] To detect \(\mathcal{H}_{0}\), the DM test statistic is defined as \[\begin{align} \mathrm{DM} = \frac{\bar{d}(a,b)}{\widehat{\sigma}(a,b)}, \end{align}\] where \(\bar{d}(a,b)\) and \(\widehat{\sigma}(a,b)\) are the mean and Newey–West standard error [55] of \(d_t(a, b)\) over the testing sample.

Notably, a significantly positive (or negative) DM statistic indicates that model \(a\) performs worse (or better) than model \(b\) under the FZ loss. Under some mild conditions, the critical value for the DM test statistic is the \((1 - \alpha/2)\)-th quantile of standard normal distribution, where \(\alpha\) is the confidence level.

10.2 MCS Test↩︎

The DM test compares models pairwise, whereas the Model Confidence Set (MCS) test of [16] evaluates a collection of competing models jointly. Starting from the full model set, the MCS procedure sequentially removes models that are significantly dominated under the chosen loss function and returns a set of models that cannot be statistically distinguished from the best-performing model(s).

Let \(\mathbb{M}^0\) denote the initial set of competing models. We use the FZ loss differential \(d_{i,t}(a,b)\) defined in the DM test section to compare any two models \(a,b\in\mathbb{M}^0\). We aggregate over the evaluation panel and define \[\begin{align} \bar d(a,b) = \frac{1}{|\mathcal{P}|} \sum_{(i,t)\in\mathcal{P}} d_{i,t}(a,b), \end{align}\] where \(\mathcal{P}\) is the set of stock-month observations in the test sample. Let \(D(a,b)=\mathbb{E}[d_{i,t}(a,b)]\). Since lower loss indicates better predictive performance, the population superior set is \[\begin{align} \mathbb{M}^{\ast} = \left\{ a\in\mathbb{M}^0: D(a,b)\le 0 \;\text{for all }\; b\in\mathbb{M}^0 \right\}. \end{align}\]

For any candidate set \(\mathbb{M}\subseteq\mathbb{M}^0\), the MCS test considers the null hypothesis \[\mathcal{H}_{0,\mathbb{M}}: D(a,b)=0, \quad \text{for all } a,b\in\mathbb{M}.\] We use the range statistic \[\begin{align} T_{\mathbb{M}} = \max_{a,b\in\mathbb{M}} \left|t(a,b)\right|, \qquad t(a,b) = \frac{\bar d(a,b)}{\widehat{\sigma}(a,b)}, \end{align}\] where \(\widehat{\sigma}(a,b)\) denotes the bootstrap standard error of \(\bar d(a,b)\).

The critical value is obtained from a panel block bootstrap as in [42]. In each bootstrap replication, we resample stocks with replacement and resample months using contiguous time blocks, so that both cross-sectional heterogeneity and serial dependence in the loss differentials are accounted for. Denote by \(\bar d^{(s)}(a,b)\) the average loss differential computed from the \(s\)-th bootstrap sample, \(s=1,\ldots,B\). We estimate the bootstrap standard error by \[\widehat{\sigma}(a,b) = \left[ \frac{1}{B} \sum_{s=1}^{B} \left\{ \bar d^{(s)}(a,b)-\bar d(a,b) \right\}^{2} \right]^{1/2}.\] The bootstrap distribution of \[T_{\mathbb{M}}^{(s)} = \max_{a,b\in\mathbb{M}} \left| \frac{\bar d^{(s)}(a,b)-\bar d(a,b)}{\widehat{\sigma}(a,b)} \right|\] is used to obtain the critical value \(CV_{\alpha,\mathbb{M}}\).

If \(T_{\mathbb{M}}\le CV_{\alpha,\mathbb{M}}\), the null is not rejected and the current set \(\mathbb{M}\) is retained. Otherwise, the model with the largest standardized average loss relative to the remaining models is removed. Specifically, the eliminated model is selected according to \[e_{\mathbb{M}} = \mathop{\mathrm{argmax}}_{a\in\mathbb{M}} \frac{\bar d_{a,\cdot}}{\widehat{\sigma}_{a,\cdot}}, \qquad \bar d_{a,\cdot} = \frac{1}{|\mathbb{M}|-1} \sum_{\substack{b\in\mathbb{M}\\ b\ne a}} \bar d(a,b),\] where \(\widehat{\sigma}_{a,\cdot}\) is the bootstrap standard error of \(\bar d_{a,\cdot}\). The procedure is repeated on \(\mathbb{M}\setminus\{e_{\mathbb{M}}\}\) until the null of equal predictive ability is no longer rejected. The final set is reported as the model confidence set.

10.3 CC Test↩︎

The CC test aims to evaluate the validity of VaR forecasts. Let \(I_{i,t} = \mathbf{1}\big\{ r_{i,t} \leq \widehat{\mathrm{VaR}}_{i,t} \big\}\), where \(\mathbf{1}(\cdot)\) denotes the indicator function. Under correct specification, the violation indicators \(\{I_{i,t}\}\) should form an i.i.d. Bernoulli sequence with success probability \(\tau\). Based on this fact, we examine the validity of VaR forecasts by checking the null hypothesis \[\begin{align} \mathcal{H}^{(v)}_{0, i}: \{I_{i,t}\} \text{ is an i.i.d. Bernoulli sequence.} \end{align}\]

Following [17], let \(T_{jk}\) denote the number of observations for which \(I_{i,t-1}=j\) and \(I_{i,t}=k\), with \(j,k\in\{0,1\}\). Define \[T_1 = T_{01} + T_{11}, \quad T_0 = T_{00} + T_{10}, \quad \hat{\pi}_0 = \frac{T_{01}}{T_{00}+T_{01}}, \quad \hat{\pi}_1 = \frac{T_{11}}{T_{10}+T_{11}}, \quad \hat{\pi} = \frac{T_{01}+T_{11}}{T^\ast},\] where \(T^\ast\) is the total number of observations in the evaluation period. Then, the CC test jointly assesses unconditional coverage and independence by combining the corresponding likelihood ratio statistics: \[\begin{align} \mathrm{LR}_{\mathrm{cc}} = \mathrm{LR}_{\mathrm{uc}} + \mathrm{LR}_{\mathrm{ind}}, \end{align}\] where \[\mathrm{LR}_{\mathrm{ind}} = -2 \log \left( \frac{(1-\hat{\pi})^{T_{00}+T_{10}} \hat{\pi}^{T_{01}+T_{11}}}{(1-\hat{\pi}_0)^{T_{00}} \hat{\pi}_0^{T_{01}} (1-\hat{\pi}_1)^{T_{10}} \hat{\pi}_1^{T_{11}}} \right)\] and \[\mathrm{LR}_{\mathrm{uc}} = -2 \log \left( \frac{(1-\tau)^{T_0}\tau^{T_1}}{(1-\hat{\pi})^{T_0}\hat{\pi}^{T_1}} \right),\] Under the null hypothesis \(\mathcal{H}^{(v)}_{0, i}\), \(\mathrm{LR}_{\mathrm{cc}}\) is asymptotically distributed as \(\chi^2(2)\). Therefore, we reject \(\mathcal{H}^{(v)}_{0, i}\) when the value of \(\mathrm{LR}_{\mathrm{cc}}\) is larger than the \(\alpha\)-th upper quantile of \(\chi^2(2)\).

10.4 AESR test↩︎

Following [18], we assess the validity of ES forecasts via the AESR test, which is based on two regressions, linking observed returns to the corresponding VaR and ES forecasts. Specifically, these regressions are defined as \[\begin{align} r_{i,t} = \beta_{i,1} + \beta_{i,2}\widehat{\mathrm{VaR}}_{i,t} + u^{(q)}_{i,t}, \qquad r_{i,t} = \gamma_{i,1} + \gamma_{i,2}\widehat{\mathrm{ES}}_{i,t} + u^{(e)}_{i,t}, \end{align}\] where \(\beta_{i,1}, \beta_{i,2}, \gamma_{i,1}\), and \(\gamma_{i,2}\) are unknown coefficients, and \(u^{(q)}_{i,t}\) and \(u^{(e)}_{i,t}\) are model errors. The AESR test makes use of the fact that, conditional on the information set \(\mathcal{F}_{t-1}\), the \(\tau\)-th conditional VaR of \(u^{(q)}_{i,t}\) and the \(\tau\)-th conditional ES of \(u^{(e)}_{i,t}\) are equal to zero almost surely.

Consequently, for each asset \(i\), the AESR test assesses the adequacy of \(\widehat{\mathrm{ES}}_{i,t}\) by testing the null hypotheses: \[\begin{align} \mathcal{H}^{(e)}_{i,0}: (\gamma_{i,1}, \gamma_{i,2})' = (0, 1)'. \end{align}\] Intuitively, under \(\mathcal{H}^{(e)}_{i,0}\), the ES forecasts are correctly specified in the sense that \(\widehat{\mathrm{ES}}_{i,t}\) behaves close to the conditional ES of \(r_{i,t}\) given the information set. Hence, for each model, rejection of \(\mathcal{H}^{(e)}_{i,0}\) at the significance level \(\alpha\) indicates that the corresponding ES forecasts fail to capture the tail distribution of returns, suggesting deficiencies in the model’s calibration.

Notably, because ES is not elicitable on its own, \(\gamma_{i,1}\) and \(\gamma_{i,2}\) cannot be estimated separately. Instead, the coefficients \((\beta_{i,1}, \beta_{i,2}, \gamma_{i,1}, \gamma_{i,2})\) are estimated jointly by minimizing the FZ loss. Let \(\widehat{\boldsymbol{\gamma}}_i = (\widehat{\gamma}_{i,1}, \widehat{\gamma}_{i,2})'\) denote the resulting estimates, and let \(\widehat{\Sigma}_{\gamma,i}\) be a consistent estimator of the asymptotic covariance matrix of \(\widehat{\boldsymbol{\gamma}}_i\). Then, we can define a Wald-type test statistic: \[\begin{align} W^{(e)}_i = \big(\widehat{\boldsymbol{\gamma}}_i - (0,1)'\big)^{\!\top} \widehat{\Sigma}_{\gamma,i}^{-1} \big(\widehat{\boldsymbol{\gamma}}_i - (0,1)'\big), \end{align}\] which asymptotically follows a \(\chi^2(2)\) distribution under \(\mathcal{H}^{(e)}_{i,0}\). The null hypothesis is rejected when \(W^{(e)}_i\) exceeds the corresponding critical value. Like CC test, the AESR test is applied to each asset \(i\) separately.

References↩︎

---

1. For more studies on VaR and ES, one can refer to [3], [4], [5], [6], [7], and [8].↩︎

2. See https://jkpfactors.s3.amazonaws.com/documents/Documentation.pdf for more details.↩︎

3. As an exception, two econometric models, GAS and GARCH, have no need for validation, so they are estimated using all training and validation data.↩︎

4. For example, in the NN model, parameter size is controlled by the hidden-state dimension and the number of hidden layers. See 7.1.2 for more details.↩︎

5. For each model, we fit a polynomial function linking loss to log parameter size, with the polynomial order selected by the Akaike information criterion [40].↩︎

6. This is why only the first \(L^\ast-1\) checkpoints are used, as \(\boldsymbol{h}_{j,t}^{S + 1}\) is not available at time \(t\).↩︎

7. For clarity, their precise definitions are deferred to Section 8.1.↩︎

8. See 8.2 for its details.↩︎

9. See 7.1.2 for its details.↩︎

10. The autoregressive structure means that the feature for time \(s\) relies only on information available up to time \(s\), for any \(s = 1, \dots, S\). In contrast, the encoder operates in a non-autoregressive manner, where the feature at time \(s\) can also depend on information from future steps \(s+1, \dots, S\).↩︎

11. See Section 8.2 for a detailed discussion on this property of the attention mechanism.↩︎