Quantifying Extreme Risks in High-Frequency Financial, Energy, and Commodity Markets

2.1 Multivariate Skewed Generalised Hyperbolic Distributions

The skewed GH distribution can be constructed as a mean-variance mixture based on a normal distribution, where the mixing component follows a Generalised Inverse Gaussian (GIG) distribution. Specifically, a random vector X = ( X 1 , … , X d ) T follows a multivariate normal mean-variance mixture distribution if

where Z ∼ N _k (0, I _k ), A ∈ R^d×k; μ and γ are vectors of parameters in R ^d representing location and skewness, respectively. When γ is zero, the distribution is symmetric around μ (see Ignatieva and Landsman 2015). W ≥ 0 is a random variable (independent of Z) with a GIG distribution, GIG(λ, χ, ψ). The probability density function (pdf) of W is given by

where w > 0, χ ≥ 0, ψ ≥ 0, λ ∈ R, and the normalisation constant c _λ is defined as

Here, K _λ (⋅) denotes the modified Bessel function of the third kind with order λ (Abramowitz and Stegun 1972). The covariance structure is captured by the matrix Σ = AA ^T , which is positive definite. Following McNeil and Embrechts (2005), the pdf of a GH distribution can be represented as

where Q(x) = (x − μ ) ^T Σ⁻¹(x − μ ). We write X ∼ GH _d (λ, χ, ψ, μ , Σ, γ ). The parameters λ, χ, and ψ are derived from the GIG mixing distribution, while μ , Σ, and γ control location, scale, and skewness, respectively. The tail behavior is influenced by λ, while χ and ψ impact skewness and kurtosis.

Another useful property is that a linear transformation of a GH random variable also yields a GH distribution (McNeil and Embrechts 2005). This feature facilitates the calculation of portfolio allocation, VaR, and TCE. For example, setting B = w ^T = (w₁, …, w _d ) and b = 0, the portfolio y = w ^T X follows a one-dimensional GH distribution:

Once the multivariate vector X ∼ GH _d (λ, χ, ψ, μ , Σ, γ ) is fitted to the data, the univariate marginals can be derived as X _i ∼ GH₁(λ, χ, ψ, μ _i , Σ _ii , γ _i ), where μ _i and γ _i denote the ith elements of the mean and skewness vectors, respectively, and Σ _ii is the ith diagonal element of the covariance matrix.

We note that the skewed GH family captures asymmetric distributional features in energy and financial markets better than Gaussian or non-skewed Student-t distributions, making it suitable for computing portfolio risk measures such as standard deviation, VaR, TCE, and other tail-related metrics analytically.

2.2 The Univariate Case and Limiting Distributions

The pdf of a univariate random variable X ∼ GH₁(λ, χ, ψ, μ, σ, γ) is derived from the general form of the multivariate GH pdf given in Eq. (2.4), expressed as:

where here | Σ | 1 2 = σ , and the term Q represents the quadratic expression Q = (x − μ)²σ⁻².

Several well-known distributions emerge as special cases of the GH distribution under specific parameter choices. For instance, the (skewed) Student-t distribution is obtained by setting λ < 0, ψ = 0, χ = ν − 2, and ν = −2λ (Praetz 1972). Similarly, the Variance Gamma (VG) distribution emerges when λ > 0, χ = 0, and ψ = 2λ (Madan and Seneta 1990). Other notable cases include the hyperbolic (HYP) distribution, derived by setting λ = (d + 1)/2 (Eberlein and Keller 1995), and the normal inverse Gaussian (NIG) distribution, obtained by choosing λ = −0.5 (Barndorff-Nielsen and Stelzer 2005). The density functions for these special cases are outlined below.

In the case of the non-symmetric Student-t distribution (with γ ≠ 0), the density is derived by substituting ψ = 0, χ = ν − 2, and ν = −2λ into Eq. (2.4), assuming λ < 0:

In cases where γ = 0, an asymptotic expansion of the Bessel function K _λ (x) transforms the density in Eq. (2.6) into the classical form of a symmetric Student-t distribution. For more details, refer to Ignatieva and Landsman (2015).

The Variance Gamma (VG) distribution emerges by setting χ = 0, ψ = 2λ, and assuming λ > 0 in Eq. (2.4):

3.1 TCE in the Univariate Case

For a portfolio loss variable X, we denote its cumulative distribution function (cdf) as F _X (x) and the tail distribution function as F ̄ X ( x ) = 1 − F X ( x ) . The tail conditional expectation, representing the average loss beyond a certain quantile, is defined by

where x _q is the quantile at level q of the loss distribution, so F ̄ X ( x q ) = 1 − q . Unlike VaR, which identifies only the minimum loss at a given confidence level, TCE captures the expected severity of losses exceeding this threshold.

For a univariate GH distribution, an explicit form for TCE can be found in Ignatieva and Landsman (2019). Let X ∼ GH₁(λ, χ, ψ, μ, σ², γ), where F GH 1 ( x ; λ , χ , ψ , μ , σ 2 , γ ) denotes the cumulative distribution function, and F ̄ GH 1 ( x ; λ , χ , ψ , μ , σ 2 , γ ) = 1 − F GH 1 ( x ; λ , χ , ψ , μ , σ 2 , γ ) is the survival function. Given x _q , the quantile for VaR at level q, we solve

We define a constant:

where expressions of k _λ for limiting cases such as the Student-t (ψ = 0) and Variance Gamma (χ = 0) distributions are provided in Ignatieva and Landsman (2019).

The TCE for the univariate GH distribution X ∼ GH₁(λ, χ, ψ, μ, σ², γ) is given by:

where f GH 1 ( ⋅ ) and F GH 1 ( ⋅ ) are the pdf and cdf of the GH random variable X ̃ ∼ GH 1 ( λ + 1 , χ , ψ , μ , σ 2 , γ ) .

3.2 Decomposition of Portfolio Risk with TCE

Consider a portfolio where the risk measure TCE is applied to a multivariate GH vector X = ( X 1 , … , X d ) T , with each component X _i representing an asset (such as an equity or energy market variable). We define the portfolio S = ∑ i = 1 d X i as the aggregate of the individual assets.^[2] If X ∼ GH _d (λ, χ, ψ, μ , Σ, γ ), then S also follows a univariate GH distribution with parameters:

where μ _i are the means of each X _i , σ _ij are elements of the covariance matrix Σ, and γ _i are components of the skewness vector γ .^[3] We denote the variance of S by σ S S = ∑ i = 1 d ∑ j = 1 d σ i j .

Using the linearity of expectation, the TCE can be distributed across individual components:

where s _q is the quantile of S at level q, and E(X _i ∣S > s _q ) represents the individual risk contributions. The contribution E(X _i ∣S > s _q ) is key in portfolio allocation. According to Ignatieva and Landsman (2019), this conditional expectation can be computed as

for i = 1, …, d. Note that ( X i , S ) ∼ GH 2 ( λ , χ , ψ , μ = ( μ i , μ S ) T , Σ i S , γ = ( γ i , γ S ) T ) , where

with σ i S = ∑ j = 1 d σ i j and γ S = ∑ j = 1 d γ j .

Summing K _i over i = 1, …, d recovers the TCE for the total portfolio:

This approach enables the TCE to decompose portfolio risk into individual asset contributions. We note that the decomposition of the multivariate TCE into individual asset contributions incurs negligible computational cost, and we clarify its empirical implementation and relevance in Section 4.3. Such risk decomposition is a powerful tool in portfolio management, allowing for more informed capital allocation based on the specific risk profiles of each asset. By identifying the assets that contribute most significantly to overall portfolio risk, this method supports the development of targeted and effective risk mitigation strategies. The flexibility of the GH distribution in modelling both asset-level and portfolio-level risks makes it a critical framework for addressing the complexities of modern risk management.

4.1 Data

The dataset includes high-frequency data from January 2, 2015, to January 28, 2022, for two commodity ETFs – United States Oil (USO) and gold (GLD) – and one equity ETF, SPY. All ETFs are traded on the NYSE Arca. USO tracks West Texas Intermediate (WTI) light sweet crude oil by primarily investing in lead-month WTI crude oil futures on the NYMEX. GLD tracks the price of gold bullion, while SPY represents the S&P 500 Index, which comprises 500 large-cap U.S. stocks. To reduce the effects of microstructure noise, we use 5-min interval quotes from 9:30 am to 4:00 pm, yielding 78 observations per trading day.

The time frame covers significant global events that heavily influenced energy, commodity, and financial markets. The 2014–2016 Oil Price Decline, driven by a supply surplus and booming U.S. shale oil production, saw crude oil prices dropped by 70 % – the largest decline since the 1980s, causing instability in oil-exporting and importing economies. In 2016, the Brexit Referendum introduced substantial uncertainty into global markets, increasing volatility across asset classes. The COVID-19 Pandemic of 2020 triggered unprecedented market disruptions, including sharp equity market declines, extreme fluctuations in commodity prices, and a brief instance of negative crude oil prices. These events profoundly affected the distributional properties of ETF returns, resulting in heightened volatility, skewness, and kurtosis. Table 1 summarises the key statistics for USO, GLD, and SPY, alongside aggregate returns. Figure 1 visualises these returns, illustrating the pronounced fluctuations and extreme events during the sample period.

Figure 1:

Returns of USO (crude oil), GLD (gold), SPY (S&P 500), and the aggregate portfolio, covering January 2, 2015, to January 28, 2022.

From Table 1, we observe that USO exhibits the highest volatility, with a standard deviation of 11.703, reflecting its pronounced price fluctuations compared to GLD (3.436) and SPY (4.259). The aggregate portfolio shows a moderate standard deviation of 4.902, indicating diversification benefits. The kurtosis values highlight the presence of extreme events, particularly for USO, which has a remarkably high kurtosis of 134.67, signaling significant outliers and extreme fluctuations. GLD and SPY exhibit kurtosis values of 20.999 and 40.785, respectively, while the aggregate portfolio’s kurtosis of 69.756 suggests residual tail risk despite diversification. Skewness further reflects the asymmetry in returns: USO (−0.403) and GLD (−0.218) have negative skewness, indicating a tendency for larger negative returns, while SPY (0.508) shows positive skewness, suggesting more frequent positive returns. The aggregate portfolio exhibits a slight negative skewness (−0.102), demonstrating that diversification reduces asymmetry.

The analysis of these ETFs over the period from January 2, 2015, to January 28, 2022, underscores the influence of extreme market events such as the COVID-19 pandemic and the 2020 oil price crash. These events contributed to heightened kurtosis and volatility, particularly for USO. The observed skewness and kurtosis emphasise the asymmetric and heavy-tailed nature of returns during this period. These findings highlight the critical importance of adopting a modelling framework capable of capturing extreme risks and asymmetric behaviours to mitigate potential losses during periods of market instability and to enhance portfolio resilience against adverse conditions.

4.2 Univariate Case

We apply maximum likelihood estimation (MLE) to fit GH distributions to univariate high-frequency loss data (i.e., negative log-returns) for USO, GLD, and SPY.^[4] For various special cases within the GH family, we fix specific parameter values corresponding to the Student-t, VG, NIG, and HYP distributions. The log-likelihood function is then maximized with respect to the remaining parameters. This approach aims to determine the GH distribution variant that best represents the data. Table 2 reports the estimated parameters for USO, GLD, SPY, and the aggregate portfolio in Panels A through D, respectively. We observe that the skewness parameter γ is negative for all ETFs and the composite portfolio, indicating left-skewed returns. In the Student-t distribution, degrees of freedom ν are expressed as ν = −2λ, which results in values ranging from 2.2 for SPY to 2.6 for GLD.

To evaluate the accuracy of our model fit, we apply three goodness-of-fit tests as proposed by D’Agostino and Stephens (1986) and Stephens (1974). These tests assess the alignment of data with a specified distribution within the GH family. Specifically, we examine the null hypothesis that the sample data corresponds to a hypothesised distribution, represented as H 0 : F ̂ n ( x ) = F n ( x ) , where F ̂ n ( x ) denotes the empirical cumulative distribution function (cdf) and F _n (x) the theoretical cdf for one of the GH family distributions, such as GH, Student-t, VG, NIG, or HYP. Test statistics are derived from z-values z _i = F _n (x _i ),^[5] where parameters are estimated as follows.

The Anderson-Darling statistic A² is calculated with the expression:

The Kolmogorov statistic D measures fit by taking the maximum of D + = max 1 ≤ i ≤ n i n − z i and D − = max 1 ≤ i ≤ n z i − i − 1 n , where a sample size adjustment factor, proposed by Stephens (1974), modifies the statistic as D ( n + 0.12 + 0.11 / n ) .

The Cramér-von Mises statistic W² is computed by

with an adjustment, reported in Stephens (1974), of W 2 − 0.4 n + 0.6 n 2 1 + 1 n .

Alongside these tests, we evaluate model quality with the Akaike information criterion (AIC) as introduced by Akaike (1974). This criterion is given by:

where l ̂ represents the maximised log-likelihood, and k is the number of estimated parameters. AIC weighs the trade-off between model fit (reflected in log-likelihood) and complexity (penalized by the parameter count), where lower AIC values suggest a better fit.

Additionally, the Bayesian information criterion (BIC) is employed, which penalizes model complexity more heavily. For BIC, the penalty term is ln(n)k, where n is the sample size:

Table 3 summarises the results for USO, GLD, SPY, and an aggregate portfolio, including AIC and BIC values. The smallest test statistics, highlighted in bold, identify the best-fitting distributions. The results indicate that none of the GH family distributions are rejected, with the GH distribution consistently achieving the lowest test statistics. This superior performance is likely due to its greater flexibility, enabled by additional parameters compared to its special cases. Furthermore, the GH distribution outperforms its special cases in terms of lower AIC and BIC values, reinforcing its suitability for our data.

Figure 2 illustrates the histogram of losses (in black) overlaid with the fitted probability density function (pdf) of the GH distribution (in red). The strong alignment between the empirical pdf and the fitted GH density highlights the GH distribution’s effectiveness in modelling the observed data. To complement this, Figure 3 presents the quantile-to-quantile (QQ) plot, comparing the empirical quantiles to those predicted by the GH distribution. The QQ plot further confirms the close agreement, demonstrating the GH distribution’s ability to capture the distribution of losses, including the tails, for both the aggregate portfolio and individual ETFs.

Figure 2:

Histogram representing the empirical distribution (black) versus GH pdf (red) fitted to USO, GLD, SPY and aggregate portfolio returns.

Figure 3:

QQ plots comparing empirical quantiles with quantiles from the GH distribution for USO, GLD, SPY, and the aggregate portfolio returns.

Using the parameter estimates from Table 2, we calculate TCEs for several distributions within the GH family. The results are summarised in Table 4 across a range of quantile levels q from 0.9 to 0.995, with Figure 4 illustrating the univariate TCEs for the ETFs under different GH family distributions. For quantile levels below 0.98, the Student-t distribution generally produces lower TCE values compared to the other distributions. However, for extreme quantiles (q ≥ 0.99), the general GH distribution yields consistently higher TCE estimates, demonstrating a more conservative approach to tail risk assessment. This pattern indicates that at the most extreme quantiles, the GH distribution offers a prudent estimation of tail risk, capturing potential losses in the distributional tail more effectively and reflecting heightened sensitivity to extreme market conditions.

Figure 4:

Comparison of TCE values calculated for USO (top-left panel), GLD (top-right panel), SPY (bottom-left panel), and the aggregated portfolio returns (bottom-right panel) across different distributions: GH (solid red), Student-t (solid yellow), VG (dotted blue), NIG (dashed purple), and HYP (dash-dotted green).

These findings underscore the critical role of the GH distribution in capturing tail risks, offering a conservative and reliable measure under extreme market conditions. This is particularly valuable for portfolios with high-volatility assets like USO, which are prone to extreme fluctuations. The elevated TCE values at higher quantiles highlight the GH distribution’s ability to account for severe market stress events, such as the 2020 oil price crash or the COVID-19 pandemic, providing enhanced protection against potential losses. By effectively addressing tail risks, the GH distribution supports the development of more resilient risk management strategies, equipping portfolios to navigate periods of pronounced market instability.

4.3 Multivariate Analysis

In this section, we analyse a portfolio comprised of USO, GLD, and SPY losses by fitting a multivariate GH distribution, along with its specific variants, to this three-dimensional portfolio. Unlike the univariate analysis in Section 4.2, where each asset’s losses were modelled independently, this approach treats the entire portfolio as a single multivariate loss distribution. Using the parameters estimated from the multivariate model, we compute the TCE for the combined three-dimensional portfolio, capturing the joint tail risks and interdependencies among the assets.

Table 5 presents the estimated parameters from fitting the multivariate GH model to the data. The model includes key parameters λ, χ, and ψ, along with a mean vector μ = ( μ 1 , μ 2 , μ 3 ) T , a skewness vector γ = ( γ 1 , γ 2 , γ 3 ) T , and a variance-covariance matrix Σ. The aggregate portfolio’s mean (μ _S ) and skewness (γ _S ) are computed as the sums of the elements in μ and γ , respectively, while the portfolio variance σ S 2 is derived from the sum of all entries in Σ. The shape parameter λ = −1.602 suggests a heavy-tailed distribution, while χ = 0.527 and ψ = 0.029 indicate significant skewness and kurtosis, capturing the pronounced tail risks inherent in the data. The components of the mean vector suggest small positive average returns, whereas the negative skewness points to a higher likelihood of large negative returns. The variance-covariance matrix reveals individual variances of σ 11 2 = 1.660 , σ 22 2 = 0.500 , and σ 33 2 = 1.280 , alongside moderate covariances such as σ₁₂ = 0.100 and σ₁₃ = 0.070, indicating stronger correlations between USO and SPY than between GLD and the other assets. Compared to univariate estimates, the multivariate model captures the interdependencies among assets more effectively, resulting in an aggregate variance σ S 2 = 3.314 and skewness γ _S = −0.283, emphasizing that while diversification lowers risk, it does not eliminate exposure to extreme market events.

To evaluate the fit of each multivariate distribution, the log-likelihood values are presented in the final row of Table 5. Among the distributions, the GH distribution achieves the highest log-likelihood, signifying the best fit to the data. These findings, consistent with the univariate analysis, reinforce the GH distribution’s ability to effectively model the complex dynamics of energy, commodity, and financial markets.

Using the parameter estimates from Table 5, we calculate TCEs for the aggregate portfolio comprising USO, GLD, and SPY losses. The results are presented in Table 6 and visualised in Figure 5 for various GH family distributions across quantile levels ranging from 0.9 to 0.999. Consistent with the log-likelihood results indicating a superior multivariate fit, the GH distribution produces larger TCE values at extreme quantiles, offering a more conservative risk estimation. This pattern aligns with the findings from the univariate TCE analysis (see Section 4.2).

Figure 5:

TCE values for the aggregate portfolio loss S, comprised of USO, GLD, and SPY losses, evaluated across various quantile levels using distributions within the GH family.

4.4 Multivariate Time-Varying Analysis

Building on the superior fitting performance of the general GH distribution for the combined loss data from the USO, GLD, and SPY portfolio, we calculate the time-varying TCE for this multivariate portfolio. We employ a rolling window approach, with a window size of 780 observations (equivalent to 10 trading days, given 78 observations per day). This approach dynamically updates the parameter estimates over time, enabling the computation of a time series of TCE values. The TCEs are calculated for quantile levels of 0.9, 0.95, 0.975, 0.99, and 0.995 (represented by red, blue, green, purple, and orange lines, respectively) and are depicted alongside the portfolio losses (grey circles) in Figure 6.

Figure 6:

Time-varying TCE for a multivariate portfolio (USO, GLD, SPY) estimated using multivariate GH distribution.

While the visual representation of our findings is compelling, it is crucial to validate the statistical robustness of our methodology. We rely on the TCE backtesting framework outlined in McNeil and Frey (2000). Our primary objective is to compare the realised loss with the estimated TCE at each point in time across various quantile levels q ∈ {0.9, 0.95, 0.975, 0.99, 0.995}. A violation occurs whenever the loss exceeds the estimated TCE. These violations are visually represented in Figure 6 as grey circles above the solid lines indicating TCEs at different quantile levels. The results of the backtesting are summarised in Table 7, which presents the actual violation rate, the expected TCE, the actual TCE, and the bootstrap p-values for each quantile level. The actual violation rate reflects the proportion of times the loss exceeded the estimated TCE during the backtesting period. For instance, at the 0.9 quantile level, the actual violation rate is 0.0667, meaning that 6.67 % of the time, the losses were greater than the TCE estimate. The expected TCE corresponds to the theoretically predicted TCE at each quantile level, while the actual TCE represents the average of the observed losses that exceeded the TCE threshold during the backtesting period. For example, at the 0.9 quantile level, the actual TCE is 9.0178, compared to an expected TCE of 10.8203. The comparison between expected and actual TCE provides insight into the model’s accuracy in predicting extreme losses. As shown in the table, the actual TCE is consistently lower than the expected TCE across all quantile levels, suggesting that the model may be somewhat conservative in its risk predictions. This pattern is consistent across all quantile levels, reinforcing the observation that the model tends to predict that extreme losses are more likely than they actually are.

To provide a more comprehensive validation, we apply a bootstrap procedure. The bootstrap test evaluates whether the mean of the exceedance residuals – defined as the difference between the realised losses and the TCE for the violations – is zero. Under the alternative hypothesis, a positive mean would indicate that the model systematically underestimates the TCE. The bootstrap test is particularly valuable as it makes no assumptions about the underlying distribution of residuals, ensuring robustness. The p-values from the bootstrap test, reported in the last row of Table 7, are all greater than 0.05, indicating that we fail to reject the null hypothesis of zero-mean of the exceedance residuals. Our results confirm that losses exceeded the TCE less frequently than anticipated, implying that the TCE estimates provided by the model are generally higher than the actually realised values. While this suggests a conservative approach, which might lead to higher capital reserves or more cautious risk management strategies, it is preferable to the alternative of underestimating risk. A conservative model ensures that extreme losses are adequately covered, thus offering better protection in risk management.