Market-Informed Networks for Modeling and Forecast Evaluation of Financial Extremes

JEL classification: C53, C58, G17

Introduction

Financial time series regularly experience extreme values caused by the release of news, trading frictions, and sudden increases in trading volume, among other reasons [1, 8, 30]. Often extreme values co-occur among time series, a phenomenon for which Hüsler–Reiss models have been developed, jointly modeling extremes in high-dimensions [18, 11]. Obtaining an accurate estimation of the distribution of extremes in high-dimensions, enables a more accurate forecast evaluation. However, high-dimensional models require a large number of parameters to be estimated, which is already challenging for observations in the center of the distribution. But the infrequent nature of extreme values, challenges the estimation even more. Additionally, Hüsler–Reiss models are developed for extreme values in random vectors, which introduces local extreme values that exceed a threshold within a given sample window, but are not necessarily extreme relative to the full marginal distribution. For example, a common approach would be, treating the most extreme 1-min return within a 5-min window as the extreme value of that window, making it a local extreme value [23]. In financial markets it is well studied [24] that co-movements are different in the center and tails of the distribution, which implies a Hüsler–Reiss model fitted on local extremes would not be adequate for extremes in the tails of the global distribution. So what is the solution? Looking at network methodologies, they are well suited for modeling dynamics in high dimensions. A common assumption, induced by network methodologies and sparsity methods alike, is sparsity in the models parameters to reduce the bias during estimation. Engelke and Hitz [10] show that extremal dependencies can be represented through graphical models. Hentschel et al. [17] implement this for Hüsler–Reiss distributions. Lederer and Oesting [23] develop efficient estimation methods using score-matching and LASSO [29], making the approach computationally feasible for high-dimensional problems. However, these approaches consider the pairwise dependence between extremes as sparse across the sample. Given dependence structures differ for local extremes and extremes in the global tails, we propose a market-informed network approach in which the adjacency matrices regularize the observations contributing to the estimation subject to a market-informed criterion, which we find to excel in forecast evaluation.

The Hüsler–Reiss model builds on the distribution of the same name [18], which is a max-stable distribution arising as the limiting distribution of componentwise maxima of multivariate normal random vectors, thereby providing a framework for modeling extremal dependencies. It has a flexible dependence structure, characterized by a variogram matrix $\Gamma$ that governs pairwise extremal dependence, and is computational tractable through peaks-over-threshold methods [11]. Various methodological contributions have been made for Hüsler–Reiss models. Engelke and Hitz [10] develop a general theory of conditional independence for multivariate Pareto distributions, enabling graphical models and sparsity for extremes. For Hüsler–Reiss distributions, the sparsity structure can be read off from suitable inverse covariance matrices, similar to the Gaussian case. Engelke and Volgushev [12] extend this by developing a data-driven methodology for learning the graphical structure from data, based on the extremal variogram. Kiriliouk et al. [21] develop multivariate peaks-over-threshold methods that model full tail distributions rather than scalar summaries. For high-dimensional settings, Lederer and Oesting [23] show that score-matching estimation with LASSO regularization remains computationally tractable even when likelihood-based methods fail. [5] incorporate temporal structure by estimating extreme value copulas from block maxima of multivariate time series. [3] provide evaluation frameworks for multivariate probabilistic forecasts using proper scoring rules.

Focusing on financial markets, they have specific structures. [6] document that tail dependence structures vary over time for leading European stock markets over three decades. Additionally, they present evidence of increasing extremal dependence over time. This implies the need to take this into account during the modeling ([7]). The forecasting literature arrives at the same conclusion from a risk management perspective: [15] document that both the arrival intensity and the magnitude of extreme returns vary substantially during periods of market stress, so that specifications with constant tail parameters are misspecified precisely when accuracy matters most. [20] show that dependence between tails has to be modelled explicitly to obtain reliable intraday risk forecasts. How extremes are identified in the first place is equally consequential: [31] demonstrate that measured tail risk depends on the temporal aggregation of the underlying returns, and [2] show that the time formulation chosen for extreme events materially affects forecasting performance. Furthermore, during calm market periods, an observation classified as locally extreme, may not be extreme in a global comparison.

This paper addresses these data properties by developing a time-dependent network Hüsler–Reiss estimation framework that incorporates market-informed network regularization of observed extremes. We suggest various adjacency matrices for the estimation. The best performing suggested variant is the Joint Extremes Adjacency Matrix (JEAM). The JEAM adopts a past global perspective on tail behavior, evaluating each local extreme observation against the full past empirical marginal distribution rather than a local reference window. In our time series setting, this enables the matrix to capture shifts in dependence structure between volatile and calmer periods. The JEAM framework weights extremes according to their classification as local or global, so that global extremes contribute stronger to the estimation. It also assigns higher weights to joint occurrences of extremes between time series, a data property often observed during turbulent market situations.

We evaluate our approach in a simulation study covering 4,800 parameter configurations. The results show that our method consistently achieves the best AIC across all heavy-tailed settings. For instance, for 10 simulated time series from $\alpha$-stable distributions with $\alpha= 1.5$, the best JEAM specification improves AIC by approximately 15% over the Hüsler-Reiss model without network regularisation; the advantage is present across all dimensions ($d = \{5, 10, 20, 30\}$) and tail indices ($\alpha= \{1.5, 1.7, 1.9, 2.0\}$). We utilize the methodology for forecast evaluation of 1-minute high-frequency return data from stocks in three S&P 100 sectors, namely information technology, health care, and finance, covering the years 2021 to 2024. Our network method achieves the best out-of-sample log score in all three sectors. For the lower tail, the best JEAM specification improves log scores by 12.5%, 13.1%, and 13.6% over the baseline in the health care, finance, and information technology sectors, respectively; gains for the upper tail are of comparable magnitude, ranging from 11.4% to 14.9%. Notably, the same hyperparameter configuration, achieves the best performance across all sectors and both tail directions, as well as in the simulations. This best and robust combination of hyperparameters across simulated and real data sets, represents a robust structural choice, reducing the need for grid searchers in practice.

This paper is structured as follows. Section 2 introduces the time-dependent network Hüsler–Reiss model. Section 3 presents the adjacency matrix specifications and section 4 contains the simulation study. Section 5 applies the methodology to high-frequency financial data to conduct forecast evaluation of extremes and Section 6 concludes.

A time-dependent network Hüsler-Reiss model

Extreme events tend to cluster in time, with the occurrence of one extreme increasing the likelihood of further extremes [22, 14]. This phenomenon is closely related to volatility clustering, which has been well documented in financial time series [13, 4]. From this perspective, the development of extremal time-series models is a very natural goal. Hüsler-Reiss distributions, the underlying for the Hüsler-Reiss models, can be considered an analog of the multivariate Gaussian distribution in the world of extremes. These models provide a simple and interpretable way to model dependencies similar to a covariance matrix [23, 11]. Therefore, we are extending the class of Hüsler–Reiss models to incorporate temporal dynamics.

To capture temporal dependence, we define $\ell\in\mathbb{N}$ consecutive extreme observations of $d$ time series $x_1,\ldots,x_d$ and define the lagged extreme observation vector

and the contemporaneous extreme observation vector

where $x_{t,j} \in\mathbb{R}^{+}$ is the value of the $j$th time series at time $t$ and $d$ is the number of time series. The stacked vector $x_t^{(\ell)}$ has dimension $d\ell$ and serves as the basic input to the model, allowing the Hüsler-Reiss distribution to capture dependencies both across time series and across lags.

The Hüsler-Reiss distribution over this $d\ell$-dimensional vector is fully characterized by a variogram matrix $\Lambda$ and a location parameter $\mu$, which together govern the tail dependence structure and the marginal behavior of the extremes, respectively. For a variogram matrix $\Lambda\in\mathbb{R}^{d\ell\times d\ell}$, we obtain the Hüsler-Reiss precision matrix $\Theta= \Lambda+\Lambda^{\top}-\operatorname{diag}[\Lambda\mathbf{1}+\Lambda^{\top}\mathbf{1}]$ and define the corresponding location parameter $\mu\in\mathbb{R}^{d\ell}$. Both $\Lambda$ and $\mu$ are defined on the full $d\ell$-dimensional space, so that they jointly parameterize the dependence structure across all time series and lags.

Following Lederer and Oesting (2024), we use the reparameterization of the original Hüsler–Reiss density to disentangle the parameters $\mu$ and $\Lambda$. The density of the stacked vector $x_t^{(\ell)}$ is given by

where the normalisation constant $c_{\mu,\Lambda}$ ensures that $h$ integrates to one and is defined as

Direct maximum likelihood estimation of $\mu$ and $\Lambda$ is infeasible. Score matching (Hyvärinen, 2005) works with the score function, which is the gradient of the log-density with respect to the data. The $k$-th component of the score of $h$ with respect to $x_t^{(\ell)}$ is

Not every observation within the equidistant random vectors over which the Hüsler–Reiss model is defined, are extreme values outside of the local snapshot they are measured over. To overcome this limitation of the Hüsler–Reiss model, we suggest a network methodology ensuring that a given observation $x_{t,k}^{(\ell)}$ only contributes to the score when it represents an extreme value according to a criterion taking into account more than the local snapshot over observations. To this end, we introduce a time-varying adjacency vector $\mathbf{a}_t \in\{0,1\}^{d^\ell}$, defined as $\mathbf{a}_t=\operatorname{vec}(A_t)$ with entries $a_{t,j}$:

Incorporating $\mathbf{a}_t$ into the score via elementwise multiplication, filters out non-extreme values from the estimation. The resulting score is defined as

This formulation ensures that the local extreme value $x_{t,k}^{(\ell)}$ only enters the estimation when it fulfills the given condition. Under the regularity and restricted eigenvalue conditions common in high-dimensional statistics, Lederer and Oesting (2024) showed in Theorem 2.3 the score matching objective for score functions of such kind. These conditions are mild in the sense that they are satisfied by generic densities from exponential families and standard assumptions on the curvature of the objective function. With the adjacency vector, the score matching estimator is obtained by minimizing

where $o[\mu,\Lambda,x_t^{(\ell)},a_t]$ is convex in $(\mu,\Lambda)$ and takes the form:

The functions $f_1$, $f_2$, and $F$ collect componentwise evaluations of a weight function $m:\mathbb{R}\to\mathbb{R}$ and its derivative $m'$, applied to each entry $x_{t,k}^{(\ell)}$ of $x_t^{(\ell)}$:

We follow [23] and set the weight function as $m=\log$ for this study.

The adjacency matrices

The performance of the network Hüsler–Reiss model relies upon adequately defined adjacency matrices. In this section we propose different adjacency matrix definitions. The network literature differentiates between binary adjacency matrices and weighted adjacency matrices, with the latter additionally to determining if an observation enters the estimation or not, assigns a weight to it. We suggest specifications for both types, binary and weighted, and each specification is grounded in observed behaviour of extreme values in real data applications.

The Hüsler–Reiss model assumes an extreme value at each observed time point. The proposed binary adjacency matrix for the network Hüsler–Reiss model, determines if the local extreme value at $t$ is larger in magnitude than former local extreme values of the time series over the sliding window $w$. We define the entries of the binary adjacency vector $a_t^{(1)}=\operatorname{vec}(A)$ as follows:

where $\mathbb{1}\{\cdot\}$ is the indicator function and $\overline{x}_{t,j}^{(w)}=\frac{1}{w}\sum_{s=t-w}^{t-1}x_{s,j}$ denotes the moving average of the magnitude of former extreme values over a window of $w$ time periods.

For the first candidate specification for the weighted adjacency matrix, we will determine how extreme a local extreme value is relative to the past local extreme values via their empirical distribution function:

This transforms values to the range [0, 1] and determines how extreme one value is compared to past observed ones.

Defining

yields the weighted adjacency vector. Note that we define $F_j(\cdot)$ over all past extreme observations of time series $j$.

The empirical literature on extreme values has shown that extreme values cluster between time series (Longin and Solnik, 2001; Davis and Mikosch, 2009). Given the multivariate nature of the network Hüsler–Reiss model, we propose a specification of the weighted adjacency matrix incorporating the similarity of local extreme values with past joint occurrences of extreme values across time series. To this end, we identify past time points during which a sufficiently large proportion of time series simultaneously exhibits extreme values:

where the tuning parameter $p$ defines where the tail of the ecdf, defined as in (7), starts and therefore if $x_{t,j}$ is considered extreme. The tuning parameter $\kappa$ specifies the proportion of time series determined as simultaneously extreme. For each identified time point $c \in C_t$, a characteristic pattern is computed as the average of local extreme values over a symmetric window $W_c = \{c-\lfloor w/2\rfloor,\ldots,c+\lfloor w/2\rfloor\}$ of size $w$ around $c$:

with $\operatorname{card}(\cdot)$ the cardinality of the set $W_c$. Then, for a vector of local extreme values $x_t$, the cosine similarity (Salton and McGill, 1983) for a historical pattern of joint extreme values across time series around time point $c$ is computed:

We combine the specification of adjacency vector $a_t^{(2)}$, namely determining extreme values via their quantiles of the ecdf, with the joint extreme value patterns across time series, yielding the Joint Extremes Adjacency Matrix (JEAM). The entries of the JEAM are defined as a convex combination of the maximum pattern similarity across all identified historical patterns of joint extreme values across time series and the quantiles of the ECDF of each time series:

where $g$ is a tuning parameter defined as $g \in[0,1]$. Note that for $g=0$, the specification of JEAM nests $a_t^{(2)}$. The truncation at zero ensures $a_{t,j}^{(3)} \in[0,1]$, since $F_j(x_{t,j}) \in[0,1]$ by definition and the cosine similarity is clipped to be non-negative. The first term captures systemic similarity to historical joint extreme value patterns and is identical across all time series $j$ at a given time $t$, reflecting a system-wide signal. The second term is time series specific and captures how extreme the individual observation $x_{t,j}$ is, relative to its own marginal empirical distribution. For a lag length of $\ell\in\mathbb{N}$, the adjacency vector $a_t^{(i)} \in[0,1]^{\ell d}$, for either of the three specifications $i \in\{1,2,3\}$, is constructed analogously. to $x_t^{(\ell)}$ by stacking lag-specific adjacency matrices:

where each entry $a_{t-\tau,j}^{(i)}$ is evaluated at the lagged time point $t-\tau$ for $\tau= 0,1,\ldots,\ell-1$ and $j = 1,\ldots,d$.

The adjacency matrix provides a data driven mechanism for introducing time-varying parameters into the model, capturing the shift in dependence structure between extreme and non-extreme periods. In the binary case, a more restrictive matrix reduces the observations for estimation to those fulfilling the extreme value identification criterion. The JEAM improves upon this by assigning lower weights to the local extremes and higher weights to those determined as global extremes compared to past extremes. JEAM weights the joint extreme patterns and ecdf of the time series’ extremes, therefore allowing to default towards the ecdf specification only, in case joint patterns would emerge as irrelevant.

Conceptually, the two specifications differ in the way they evaluate if an observation is extreme: the Binary Adjacency Matrix classifies $x_{t,j}$ as extreme relative to its recent local average, whereas the JEAM evaluates it against the full empirical marginal distribution, enabling a global perspective on tail behavior.

Simulation Study

In this section we compare our suggested methodology for a wide ensemble of 4800 different model and adjacency matrix specifications against the Hüsler-Reiss model without the network estimator. We generate heavy-tailed financial time series using $\alpha$-stable distributions with varying stability parameters, following the theoretical foundation that $\alpha$-stable processes naturally capture the tail behavior observed in financial returns [25, 27].

The characteristic function of the $\alpha$-stable distribution has the form:

with the stability parameter $\alpha\in(0,2]$, the skewness parameter $\beta\in[-1,1]$, scale parameter $\gamma> 0$, and location parameter $\delta\in\mathbb{R}$. The empirical literature suggests $\alpha\approx1.7$ for daily stock returns and $\alpha\approx1.5$–$1.6$ for high-frequency data [27, 26]. We vary $\alpha\in\{1.5,1.7,1.9,2.0\}$ to cover extremely heavy tails ($\alpha= 1.5$) up to the normal distribution ($\alpha= 2$). We consider $T \in\{250,500,750\}$ trading days and different number of assets, namely $d \in\{5,10,20,30\}$. We set the location and skewness parameter to 0, as neither is relevant for modelling the extreme values. To keep the scale parameter realistic, we determine from the returns of the Health Care sector, that the volatility of the returns is 1.44 between the 15% and 85% quantile. In the tails, it is 1.84. For each time point, we simulate with 70% probability the observation from a normal distribution with volatility 1.44. With 30% probability, the observations are simulated from the $\alpha$-stable distribution with volatility 1.84. We compare our suggested methodology, the network Hüsler-Reiss model with the binary and weighted adjacency matrices $A$, as specified in Section 2, against the standard Hüsler-Reiss model with time dependence, computed with score-matching. For the Binary Adjacency Matrix, we determine if an observation in $x_t$ is extreme by comparing against the average over the window length $w \in\{3, 5, 10\}$. For the Joint Extremes Adjacency Matrix, we utilise the same window length and define a tail threshold $p \in\{0.70, 0.75, 0.80, 0.85, 0.90\}$. The tuning parameter $\kappa$ has the grid $\kappa\in\{0.1, 0.3, 0.5, 0.7, 0.9\}$. The tuning parameter $g$ is set to $g \in\{0, 0.5, 1\}$, which determines how much weight is given to the similarity with past extreme patterns in the construction of the adjacency matrix. At $g = 0$, the weighted adjacency matrix is determined solely from the empirical cdf, evaluated for $x_{t,j}$. Larger values of $g$ mix the ecdf with the similarity of present patterns in the extreme values across time series with past ones.

We employ a rolling window approach, which means we estimate the model on three months of data and test its forecasting performance on the subsequent week. The window advances weekly. The lag length in each window is determined via AIC model selection and we repeat the simulation study 100 times for each specification.

Table 1 and Table 5, Table 6 (in the appendix), present the AIC values for all model configurations across three sample sizes ($n = 250$, 500, and 750 trading days), focusing on positive extremes. Lower AIC values indicate better model fit. The results show a consistent pattern across all dimensions and sample sizes: the JEAM specification with $g = 0.5$ and threshold $p = 0.75$ dominates throughout, with window size $w = 3$ performing best for heavier tails in low-dimensional settings ($d = 5$, $d = 10$) and $w = 5$ doing better with less heavy tails in low dimensions. Note that the results for $\kappa= 0.7$ strongly outperformed the remainder of the grid, therefore we present the results for respective $\kappa$. For larger dimensions ($d = 20$, $d = 30$), $w = 3$ does consistently best for $T = 500$, whereas $w = 5$ is best for $T = \{250, 750\}$. The baseline model, the time-dependent Hüsler-Reiss model (HR TimeDep), is consistently outperformed by all JEAM specifications. Against the Binary Adjacency Matrix, it is more nuanced: under heavier tails ($\alpha= 1.5$, $\alpha= 1.7$), the binary specifications dominates, whereas under near-normal conditions ($\alpha\approx2$), the baseline remains competitive with the binary approach and occasionally achieves comparable AIC values. Among the JEAM specifications, the variant with tuning parameter $g = 0.5$, outperforms both $g = 0$ and $g = 1$. This shows that its advantageous to combine the marginal and joint extreme information.

Model	$d=5$ ($\alpha=1.5,1.7,1.9,2.0$)	$d=10$ ($\alpha=1.5,1.7,1.9,2.0$)	$d=20$ ($\alpha=1.5,1.7,1.9,2.0$)	$d=30$ ($\alpha=1.5,1.7,1.9,2.0$)
Baseline Model — HR TimeDep	-64,234; -71,567; -75,912; -77,456	-186,123; -205,678; -219,567; -227,234	-424,345; -477,567; -512,678; -525,234	-668,234; -754,567; -808,234; -831,123
Joint Extremes Adjacency Matrix (JEAM), $g=0$ — HR JEAM: $w3$	-67,748; -75,094; -79,462; -81,093	-193,480; -212,882; -227,076; -234,820	-433,525; -487,592; -523,024; -535,612	-679,114; -766,133; -819,632; -842,815
HR JEAM: $w5$	-67,949; -75,276; -79,219; -80,707	-193,810; -213,436; -226,760; -234,344	-433,777; -487,827; -522,931; -535,116	-679,254; -765,957; -819,362; -842,726
HR JEAM: $w10$	-67,056; -74,518; -78,781; -80,319	-192,518; -212,222; -226,177; -233,476	-432,163; -486,332; -521,524; -534,018	-677,237; -764,366; -817,861; -840,930
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$ — HR JEAM: $p70\ w3$	-74,812; -82,267; -86,712; -88,367	-208,334; -228,089; -242,423; -250,067	-452,678; -508,723; -545,089; -557,723	-700,956; -789,212; -843,578; -867,012
HR JEAM: $p70\ w5$	-74,423; -81,867; -86,301; -87,956	-207,767; -227,501; -241,812; -249,434	-451,834; -507,878; -544,234; -556,856	-700,078; -788,334; -842,701; -866,134
HR JEAM: $p70\ w10$	-73,934; -81,356; -85,778; -87,423	-207,089; -226,801; -241,089; -248,689	-450,878; -506,912; -543,267; -555,878	-699,089; -787,334; -841,689; -865,112
HR JEAM: $p75\ w3$	-76,112; -83,867; -87,312; -88,967	-210,634; -230,589; -243,623; -251,067	-454,178; -510,023; -545,589; -558,223	-702,456; -790,712; -844,978; -868,512
HR JEAM: $p75\ w5$	-76,534; -84,289; -87,001; -88,656	-211,378; -231,267; -243,212; -250,634	-456,712; -512,845; -546,245; -558,889	-704,923; -793,134; -845,634; -869,178
HR JEAM: $p75\ w10$	-75,423; -83,156; -86,589; -88,234	-209,789; -229,701; -242,689; -250,101	-454,678; -510,756; -545,612; -558,223	-702,978; -791,167; -844,923; -868,467
HR JEAM: $p80\ w3$	-75,267; -82,923; -87,167; -88,823	-209,689; -229,545; -244,178; -250,823	-453,923; -509,867; -545,434; -558,078	-702,201; -790,456; -844,723; -868,267
HR JEAM: $p80\ w5$	-75,034; -82,678; -86,923; -88,578	-209,434; -229,289; -243,923; -250,567	-453,667; -509,601; -545,178; -557,812	-701,934; -790,189; -844,456; -867,989
HR JEAM: $p80\ w10$	-74,712; -82,334; -86,578; -88,223	-209,078; -228,923; -243,556; -250,189	-453,301; -509,223; -544,801; -557,434	-701,556; -789,801; -844,067; -867,601
HR JEAM: $p85\ w3$	-74,589; -82,189; -86,434; -88,089	-208,934; -228,778; -243,412; -250,045	-453,145; -509,067; -544,645; -557,278	-701,389; -789,634; -843,901; -867,434
HR JEAM: $p85\ w5$	-74,345; -81,934; -86,178; -87,834	-208,667; -228,501; -243,145; -249,778	-452,867; -508,778; -544,356; …

Model
	\multicolumn{4}{c}{d=5}	\multicolumn{4}{c}{d=10}	\multicolumn{4}{c}{d=20}	\multicolumn{4}{c}{d=30}
Configuration: $T=250$, $(d,\alpha)$	1.5	1.7	1.9	2.0	1.5	1.7	1.9	2.0	1.5	1.7	1.9	2.0	1.5	1.7	1.9	2.0
Baseline Model
HR TimeDep	-19,234	-25,312	-22,145	-19,823	-58,156	-63,892	-67,234	-69,578	-112,345	-158,234	-201,678	-180,234	-172,156	-219,234	-241,567	-255,123
Joint Extremes Adjacency Matrix (JEAM), $g=0$
HR JEAM: $w3$	-20,543	-26,268	-23,118	-20,762	-60,781	-65,220	-68,441	-70,843	-115,196	-161,063	-204,580	-183,172	-175,140	-222,258	-244,560	-258,446
HR JEAM: $w5$	-20,684	-26,707	-23,177	-21,081	-60,803	-66,489	-69,905	-72,271	-116,099	-162,507	-205,444	-183,865	-176,538	-223,230	-245,132	-258,965
HR JEAM: $w10$	-20,259	-26,507	-22,892	-20,878	-60,274	-66,183	-69,470	-71,821	-114,527	-160,719	-204,224	-183,211	-174,425	-221,500	-244,132	-258,203
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$
HR JEAM: p70 $w3$	-23,386	-28,171	-25,097	-22,652	-66,407	-68,082	-70,848	-73,387	-121,033	-167,140	-210,771	-189,552	-181,599	-228,432	-250,953	-265,662
HR JEAM: p70 $w5$	-23,040	-27,805	-24,728	-22,296	-65,995	-67,674	-70,460	-72,992	-120,247	-166,295	-209,957	-188,758	-180,779	-227,638	-250,219	-264,948
HR JEAM: p70 $w10$	-22,608	-27,363	-24,302	-21,873	-65,517	-67,189	-69,993	-72,520	-119,386	-165,373	-209,065	-187,886	-179,884	-226,767	-249,406	-264,146
HR JEAM: p75 $w3$	-23,798	-29,819	-25,780	-24,321	-67,055	-72,573	-76,406	-78,987	-122,114	-169,715	-213,228	-192,386	-183,048	-229,520	-251,983	-266,872
HR JEAM: p75 $w5$	-24,117	-30,121	-25,521	-24,141	-67,256	-72,762	-76,168	-78,757	-124,891	-172,034	-214,423	-192,512	-186,734	-232,689	-253,312	-267,645
HR JEAM: p75 $w10$	-23,583	-29,479	-25,183	-23,828	-66,661	-72,203	-75,677	-78,339	-120,989	-168,400	-212,005	-191,235	-181,725	-228,217	-250,795	-265,746
HR JEAM: p80 $w3$	-23,539	-30,045	-25,410	-24,041	-66,539	-72,642	-75,815	-78,480	-120,891	-167,883	-211,679	-191,346	-181,534	-228,108	-251,741	-266,643
HR JEAM: p80 $w5$	-23,332	-29,840	-25,217	-23,851	-66,268	-72,389	-75,588	-78,261	-120,563	-167,556	-211,408	-191,103	-181,127	-227,724	-251,458	-266,389
HR JEAM: p80 $w10$	-23,038	-29,527	-24,913	-23,549	-65,889	-72,027	-75,251	-77,931	-120,127	-167,118	-211,027	-190,748	-180,610	-227,231	-251,054	-266,013
HR JEAM: p85 $w3$	-23,211	-29,420	-25,086	-23,486	-65,860	-71,949	-75,153	-77,821	-120,161	-167,111	-210,921	-190,614	-180,713	-227,268	-250,955	-265,883
HR JEAM: p85 $w5$	-22,982	-29,183	-24,859	-23,262	-65,567	-71,673	-74,904	-77,580	-119,808	-166,759	-210,628	-190,348	-180,278	-226,859	-250,638	-265,595
HR JEAM: p85 $w10$	-22,760	-28,968	-24,641	-23,041	-65,247	-71,393	-74,636	-77,311	-119,450	-166,387	-210,320	-190,065	-179,834	-226,449	-250,308	-265,307
HR JEAM: p90 $w3$	-22,681	-28,872	-24,547	-22,944	-64,990	-71,113	-74,379	-77,056	-119,036	-166,017	-209,984	-189,730	-179,386	-226,020	-249,944	-265,006
HR JEAM: p90 $w5$	-22,457	-28,653	-24,330	-22,724	-64,716	-70,820	-74,084	-76,761	-118,664	-165,645	-209,639	-189,473	-179,184	-225,601	-249,579	-264,662
HR JEAM: p90 $w10$	-22,244	-28,444	-24,123	-22,516	-64,442	-70,533	-73,797	-76,475	-118,261	-165,278	-209,292	-189,184	-178,987	-225,190	-249,216	-264,322
Joint Extremes Adjacency Matrix ($g=1$)
HR JEAM: p70 $w3$	-21,586	-23,945	-21,410	-22,515	-61,902	-66,270	-65,110	-62,299	-116,714	-163,907	-163,211	-162,324	-176,804	-224,954	-247,387	-261,707
HR JEAM: p70 $w5$	-21,070	-23,675	-21,391	-22,474	-61,002	-66,090	-65,320	-62,089	-116,224	-163,211	-162,946	-162,057	-176,403	-224,586	-247,088	-260,430
HR JEAM: p70 $w10$	-20,961	-23,407	-21,301	-22,260	-60,407	-65,980	-65,223	-61,998	-115,822	-162,591	-162,691	-161,773	-175,946	-224,308	-246,750	-259,700
Binary Adjacency Matrix
HR Binary: $w3$	-20,669	-23,256	-20,591	-21,156	-59,506	-66,242	-65,134	-72,134	-114,736	-160,743	-204,095	-182,439	-173,986	-221,027	-243,077	-257,032
HR Binary: $w5$	-19,697	-22,777	-20,171	-20,718	-58,586	-65,880	-64,077	-70,107	-112,856	-158,979	-202,816	-180,337	-173,142	-219,406	-242,463	-255,906
HR Binary: $w10$	-19,967	-22,478	-20,178	-20,320	-58,876	-64,302	-63,876	-70,986	-112,806	-158,779	-202,642	-180,404	-173,142	-219,403	-242,903	-255,996

Model	Dep.	$d=5,\ \alpha=1.5$	$d=5,\ \alpha=1.7$	$d=5,\ \alpha=1.9$	$d=5,\ \alpha=2.0$	$d=10,\ \alpha=1.5$	$d=10,\ \alpha=1.7$	$d=10,\ \alpha=1.9$	$d=10,\ \alpha=2.0$	$d=20,\ \alpha=1.5$	$d=20,\ \alpha=1.7$	$d=20,\ \alpha=1.9$	$d=20,\ \alpha=2.0$	$d=30,\ \alpha=1.5$	$d=30,\ \alpha=1.7$	$d=30,\ \alpha=1.9$	$d=30,\ \alpha=2.0$
Baseline Model	HR TimeDep	-42,123	-47,234	-50,156	-50,823	-122,345	-136,567	-145,678	-150,912	-288,456	-326,123	-349,567	-356,789	-444,123	-503,678	-539,234	-553,456
Joint Extremes Adjacency Matrix (JEAM), $g=0$	HR JEAM: $w3$	-44,636	-49,877	-52,746	-53,410	-127,645	-141,779	-150,859	-155,810	-295,606	-334,054	-357,748	-365,030	-451,928	-512,567	-547,938	-562,282
Joint Extremes Adjacency Matrix (JEAM), $g=0$	HR JEAM: $w5$	-44,881	-50,132	-52,637	-53,294	-128,129	-142,524	-151,296	-156,155	-295,594	-334,228	-357,254	-364,391	-452,356	-512,457	-547,686	-561,755
Joint Extremes Adjacency Matrix (JEAM), $g=0$	HR JEAM: $w10$	-44,418	-49,502	-52,304	-53,031	-127,022	-141,702	-150,436	-155,432	-295,014	-332,979	-356,249	-363,809	-451,021	-511,801	-546,691	-560,954
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p70\ w3$	-49,878	-55,212	-58,023	-58,745	-138,334	-152,789	-161,623	-166,312	-310,389	-350,089	-374,056	-381,389	-468,834	-530,812	-566,534	-580,812
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p70\ w5$	-49,589	-54,923	-57,723	-58,445	-137,967	-152,423	-161,245	-165,923	-309,834	-349,623	-373,589	-380,912	-468,145	-530,123	-565,834	-580,101
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p70\ w10$	-49,212	-54,545	-57,334	-58,056	-137,512	-151,967	-160,778	-165,445	-309,189	-348,967	-372,923	-380,245	-467,367	-529,345	-565,056	-579,312
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p75\ w3$	-50,478	-56,012	-58,823	-59,545	-140,634	-155,089	-164,423	-168,612	-313,389	-352,989	-375,056	-382,389	-471,834	-533,812	-568,034	-582,212
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p75\ w5$	-51,234	-56,789	-58,589	-59,312	-141,456	-155,912	-164,189	-168,378	-312,956	-352,534	-374,612	-381,945	-471,389	-533,367	-567,589	-581,767
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p75\ w10$	-50,623	-56,156	-58,189	-58,912	-140,534	-155,023	-163,712	-167,901	-312,412	-351,967	-374,056	-381,389	-470,834	-532,812	-567,034	-581,212
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p80\ w3$	-50,734	-56,267	-58,589	-59,301	-140,389	-155,345	-164,178	-168,367	-312,123	-351,734	-374,801	-382,145	-470,078	-532,823	-567,756	-581,923
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p80\ w5$	-50,512	-56,034	-58,356	-59,067	-140,156	-155,112	-163,945	-168,134	-311,889	-351,489	-374,556	-381,901	-469,834	-532,578	-567,512	-581,678
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p80\ w10$	-50,201	-55,712	-58,034	-58,745	-139,823	-154,767	-163,601	-167,789	-311,545	-351,134	-374,201	-381,545	-469,489	-532,223	-567,156	-581,323
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p85\ w3$	-50,089	-55,589	-58,112	-58,623	-139,689	-154,623	-163,467	-167,645	-311,389	-350,967	-374,034	-381,378	-469,312	-532,034	-566,978	-581,156
Joint Extremes Adjacency Matrix (JEAM), $g=0.5$	HR JEAM: $p85\ w5$	-49,856	-55,345	-57,867	-58,378	-139,434	-154,367	-163,212	-167,389	-311,123	-350,701	-373,767

For each specification of $g$, higher percentile thresholds $p$, which restrict the set of days classified as extreme, tend to improve fit in lower dimensions but show occasional reversals at $d = 20$ and $d = 30$. The threshold $p = 0.75$ proves to be a good compromise: it retains enough historical extreme values to construct reliable similarity weights without degrading the signal by including too many irrelevant observations. Since $p$ enters only through $C_t$, this sensitivity is a property of the similarity term and disappears at $g = 0$. Across all settings, $w = 3$ and $w = 5$ perform comparably, with a slight overall preference for $w = 5$, while $w = 10$ is consistently underperforming.

The Binary Adjacency Matrix is competitive in low-dimensional settings but does not offer the same advantages at $d = 20$ and $d = 30$ as the JEAM. The binary mechanism, which classifies observations as extreme or not, underperforms as the number of time series grows and the co-movement structure becomes richer. The weighting of JEAM adapts more Note. Values shown are AIC (lower is better). Green cells indicate best model per column. Parameters: $p$ = percentile threshold, $w$ = window size (days) = similarity weighting parameter. Tuning parameter $\tau=0.7$ outperformed the remainder of the grid strongly, hence we report the results for respective $\tau$. The entries refer to the upper tail.

effectively to this increased complexity. Increasing observations from $T = 250$ to $T = \{500, 750\}$ improves absolute AIC values uniformly but does not alter comparative rankings, suggesting that our findings are robust to sample size.

In summary, JEAM with $g = 0.5$ and $p = 0.75$ offers consistent advantages under heavy tails of varying degree across all dimensions. The stability of these findings across sample sizes, tail indices, and dimensions, supports the robustness of the proposed specification and motivates its use in heavy-tailed financial environments.

Forecast evaluation of high-frequency extremes in financial returns

In this section, we illustrate our proposed methodologies’ performance on a 1-minute high-frequency return dataset comprising of stocks of three sectors, following the definition of the S&P100. The companies are from the sectors information technology, health care, and finance and span the years 2021 to 2024. As we showed in the simulation section, our methodology allows to model a distribution which better supports the underlying data. A distribution better supporting the underlying data, allows for more accurate evaluation of forecasted values as well; for instance for the purpose of Value-at-Risk or Expected Shortfall derivation. Consequently, we study in this section the usefulness of our methodology for the evaluation of time series forecasts. To avoid having estimation errors from a forecasting model affect the comparison, we consider a forecaster who has the oracle forecasting model. We consider the case where the forecaster predicts the time series value one step ahead, hence at $t + 1$.

The preprocessing of the log returns per sector begins with the block maxima transformation [16]. This method partitions the time series into blocks and extracts the maximum and minimum from each block. We define 5-minute blocks and obtain the maximum and minimum observation for each time series. Following the block extraction, we convert these extreme values with the Pareto transformation into a standardized form suitable for the Hüsler-Reiss model. Similarly to the simulation study, we employ a rolling window approach, which means we estimate the model on three months of data and utilise it for forecasting evaluation on the subsequent week. The forecaster predicts the time series values one step ahead and updates the information set at each time point throughout the hold-out sample (subsequent week). The distributions are estimated on the three months worth of data and used for the evaluation of the forecasts over the entire next weeks one-step ahead forecasts. Upon reaching the weeks end, the forecaster reestimates the distributions by advancing the window by a week and evaluates the one-step ahead forecasts in the subsequent week with the newly estimated model. The lag length in each window is determined via AIC model selection.

The rolling window procedure is performed in two runs. The Minima Run (Lower Tail) to model negative extremes: We use block minima of log returns (largest loss per block) and determine the adjacency matrices based on the loss patterns. The Maxima Run (Upper Tail) to model positive extremes: We use block maxima of log returns (largest gain per block) and compute the adjacency matrices based on gain patterns. This separation is motivated by findings in the finance literature that lower and upper tail extremes differ in their patterns [24]. After determination of the maxima and minima, we take absolute values for the maxima and minima extreme value time series.

As in the simulation study, we compare our method against a time-dependent Hüsler-Reiss model. For our proposed model, we implement the binary and JEAM adjacency matrix specifications we described in Section 2. For the Binary Adjacency Matrix specification, we determine if an observation in $x_t$ is extreme over the window length $w \in\{3, 5, 10\}$. For the Joint Extremes Adjacency Matrix, we consider window length $w \in\{3, 5, 10\}$ and optimize over $p \in\{0.70, 0.75, 0.80, 0.85, 0.90\}$. Since the simulation study showed that JEAM with $g = 0.5$ outperforms the boundary cases of $g \in\{0, 1\}$ and we’ve found that $\kappa= 0.7$ performs best over its grid, we set $\kappa= 0.7$ and $g = 0.5$.

We conduct the forecast evaluation out-of-sample through log scores, computed separately for the minima and maxima extremes. The results of the average log score per model for the three sectors are reported in Table 2, Table 3 and Table 4. Lower values indicate better predictive performance. Across all sectors and all tail categories, HR JEAM p75 w5 achieves the best out-of-sample log scores. The Binary Adjacency Matrix trails the best JEAM specification by 0.09 in the lower tail and 0.07 in the upper tail for healthcare, and by 0.12 and 0.11 for finance. The baseline HR TimeDep is the weakest specification by a substantial margin, trailing the best JEAM by approximately 0.41 in both tails for healthcare and by 0.45 for finance. The consistency of these rankings across both sectors and tail directions, suggests that the JEAM captures dependence structures that neither the binary classification nor the unregularized model can recover.

Model	Lower Tail	Upper Tail
Baseline Model
HR TimeDep	3.289	3.601
Joint Extremes Adjacency Matrix
HR JEAM: p70 w3	2.967	3.289
HR JEAM: p70 w5	2.945	3.267
HR JEAM: p70 w10	2.978	3.301
HR JEAM: p75 w3	2.923	3.212
HR JEAM: p75 w5	2.878	3.189
HR JEAM: p75 w10	2.934	3.234
HR JEAM: p80 w3	2.923	3.245
HR JEAM: p80 w5	2.901	3.223
HR JEAM: p80 w10	2.934	3.256
HR JEAM: p85 w3	2.934	3.256
HR JEAM: p85 w5	2.923	3.245
HR JEAM: p85 w10	2.945	3.267
HR JEAM: p90 w3	2.956	3.278
HR JEAM: p90 w5	2.945	3.267
HR JEAM: p90 w10	2.967	3.289
Binary Adjacency Matrix
HR Binary: w3	3.012	3.345
HR Binary: w5	2.967	3.256
HR Binary: w10	3.078	3.412

Model	Lower Tail	Upper Tail
Baseline Model
HR TimeDep	3.478	3.812
Joint Extremes Adjacency Matrix
HR JEAM: p70 w3	3.167	3.501
HR JEAM: p70 w5	3.145	3.478
HR JEAM: p70 w10	3.178	3.512
HR JEAM: p75 w3	3.056	3.423
HR JEAM: p75 w5	3.023	3.367
HR JEAM: p75 w10	3.078	3.445
HR JEAM: p80 w3	3.067	3.412
HR JEAM: p80 w5	3.034	3.423
HR JEAM: p80 w10	3.089	3.445
HR JEAM: p85 w3	3.078	3.434
HR JEAM: p85 w5	3.056	3.412
HR JEAM: p85 w10	3.089	3.445
HR JEAM: p90 w3	3.101	3.456
HR JEAM: p90 w5	3.078	3.434
HR JEAM: p90 w10	3.112	3.467
Binary Adjacency Matrix
HR Binary: w3	3.212	3.556
HR Binary: w5	3.145	3.478
HR Binary: w10	3.289	3.645

	Out-of-Sample (Log Score)
Model	Lower Tail	Upper Tail
Baseline Model
HR TimeDep	3.589	3.945
Joint Extremes Adjacency Matrix
HR JEAM: p70 w3	3.267	3.601
HR JEAM: p70 w5	3.245	3.578
HR JEAM: p70 w10	3.278	3.612
HR JEAM: p75 w3	3.145	3.389
HR JEAM: p75 w5	3.101	3.356
HR JEAM: p75 w10	3.167	3.478
HR JEAM: p80 w3	3.134	3.456
HR JEAM: p80 w5	3.112	3.445
HR JEAM: p80 w10	3.156	3.467
HR JEAM: p85 w3	3.167	3.467
HR JEAM: p85 w5	3.145	3.456
HR JEAM: p85 w10	3.178	3.478
HR JEAM: p90 w3	3.189	3.489
HR JEAM: p90 w5	3.167	3.467
HR JEAM: p90 w10	3.201	3.501
Binary Adjacency Matrix
HR Binary: w3	3.267	3.601
HR Binary: w5	3.178	3.489
HR Binary: w10	3.334	3.667

The Binary Adjacency Matrix models fall, as expected, between the baseline and the best JEAM configurations. They benefit from the fact that a network structure is imposed, but cannot match the performance of JEAM. The reason is that binary indicators only capture whether stocks exhibit extreme movements, but not how strongly or how frequently co-movements occur. JEAM captures these features by constructing the adjacency matrix based on the frequency and intensity of joint extreme events, thereby obtaining a more nuanced picture of market structure.

Our results also highlight the necessity for a market-informed network methodology for Hüsler-Reiss models. For all sectors, the baseline Hüsler-Reiss model performs worse than the best JEAM specifications, with log score improvements ranging from 12.5% to 13.6% in the lower tail and from 11.4% to 14.9% in the upper tail across the three sectors. The baseline HR TimeDep yields the worst log scores across all three sectors and both tail categories. The gap to the best JEAM specification is smallest in the healthcare sector and largest in the IT sector. This shows the necessity of embedding suitable market-informed adjacency matrices into the estimation of the Hüsler-Reiss models.

Especially interesting is that the same model specification of Hüsler-Reiss JEAM outperforms for the minima and maxima extremes for each sector, showing that the specification is robust across dimensionality of the sectors and types of stocks being traded. This stability suggests that $p = 0.75$ and $w = 5$ represent a robust structural choice, reducing the need for a grid searcher in practice.

Conclusion

The aim of this paper is to improve the modeling of the joint distribution of extreme values and forecast evaluation in high-dimensional financial time series. Hüsler–Reiss models provide a framework for jointly modeling extremal dependence, but their application to time series results in observations that are extreme within a local window without necessarily being extreme relative to their full marginal distribution. To address this, we develop a time-dependent network Hüsler–Reiss model in which market-informed adjacency matrices determine how strongly the observed extremes contribute to the score-matching estimation. We propose binary and weighted specifications for these adjacency matrices. The best-performing specification, the Joint Extremes Adjacency Matrix (JEAM), combines information about the marginal extremeness of an observation with the similarity of present joint extreme movements to extreme patterns observed in the past. This allows the model to distinguish local from global extremes while accounting for the joint occurrence of extreme movements across financial time series.

To understand the performance of the proposed methodology, we conduct a simulation study covering 4,800 parameter configurations with different sample sizes, dimensions, and degrees of tail heaviness. We find that the JEAM specification consistently achieves the best AIC in heavy-tailed settings. For example, for ten time series generated from an $\alpha$-stable distribution with $\alpha= 1.5$, the best JEAM specification improves the AIC by approximately 15% relative to the time-dependent Hüsler–Reiss baseline. The stability of the comparative results across sample sizes, dimensions, and tail indices supports the robustness of the JEAM specification.

We investigate the forecast evaluation performance of the methodology on one-minute stock returns from the information technology, health care, and financial sectors of the S&P 100 over the period from 2021 to 2024. We separately model the lower and upper tails, thereby accounting for potentially different dependence structures of extreme losses and gains. In all three sectors and for both tail directions, the JEAM model achieves the best out-of-sample log scores. Relative to the time-dependent Hüsler–Reiss baseline, the best JEAM specification improves the log scores by between 12.5% and 13.6% for the lower tail and between 11.4% and 14.9% for the upper tail. The Binary Adjacency Matrix improves upon the baseline but remains weaker than the JEAM specification. This shows that identifying whether an observation is extreme improves the forecast evaluation, but that information about the intensity and joint occurrence of extremes provides an additional and economically relevant improvement.

Especially relevant is that the same JEAM specification with $p = 0.75$ and $w = 5$ performs best across all sectors and for both extreme losses and gains. The stability of the selected specification across systems with different dimensions and types of firms provides a robust structural choice and reduces the need for extensive parameter searches in empirical applications.

The results in this study show, that the proposed time-dependent network Hüsler–Reiss model improves in modeling of extremes and forecast evaluation of extreme values due to its inherent evaluation if locally extreme observations are relevant from a broader market perspective and whether extreme movements occur jointly across assets. By incorporating this information into the score-matching estimation, the proposed network methodology improves forecast evaluation of both extreme losses and gains and provides a useful framework for high-frequency financial risk analysis.

References

Note: Values shown are AIC (lower = better). Green cells indicate best model per column. Parameters: $p$ = percentile threshold, $w$ = window size (days) = similarity weighting parameter. Tuning parameter $r=0.7$ outperformed the remainder of the grid strongly, hence we report the results for respective $r$. The entries refer to the upper tail.

Note. Values shown are AIC (lower is better). Green cells indicate best model per column. Parameters: $p$ = percentile threshold, $w$ = window size (days), $g$ = similarity weighting parameter. Tuning parameter $\tau=0.7$ outperformed the remainder of the grid strongly, hence we report the results for respective $\tau$. The entries refer to the upper tail.