Quantile-based modeling of scale dynamics in financial returns for Value-at-Risk and Expected Shortfall forecasting

\thispagestyle{empty}

Abstract: We introduce a semiparametric approach for forecasting {Value-at-Risk} (VaR) and {Expected Shortfall} (ES) by modeling the conditional scale of financial returns, defined as the difference between two specified quantiles, via restricted quantile regression. Focusing on downside risk, VaR is derived from the left-tail quantile of rescaled returns, and ES is approximated by averaging quantiles below the VaR level. The method delivers robust, distribution-free estimates of extreme losses and captures skewness, heavy tails, and leverage effects. Simulation experiments and empirical analysis show that it often outperforms established models, including GARCH and joint VaR-ES conditional-quantile approaches. An application to daily returns on major international stock indices, spanning the COVID-19 period, highlights its effectiveness in capturing risk dynamics.

{\bf JEL classification:} C14, C22, G17, G32

{\bf Keywords:} Conditional scale dynamics, {CAViaR} models, Multiple quantiles, Robust risk estimation, Financial risk forecasting

\thispagestyle{empty}

\pagenumbering{arabic}

\setcounter{footnote}{0}

Introduction

Forecasting risk measures such as {Value-at-Risk} (VaR) and {Expected Shortfall} (ES) is crucial for financial institutions, risk managers, and regulators. VaR, which estimates a conditional quantile in the lower tail of the return distribution, has long been a standard tool for assessing potential losses. Specifically, the VaR for period $t$ at probability level $\alpha$ is defined as

where $r_{t}$ denotes the asset return at time $t$ and $Q_{r_{t}}(\alpha \,\vert \, \mathcal{I}_{t-1})$ is the $\alpha$th conditional quantile given the information set $\mathcal{I}_{t-1}$. However, VaR has notable limitations, particularly its inability to account for extreme losses beyond the quantile threshold defined by $\alpha$, and its potential failure to capture diversification benefits Embrechts-McNeil-Straumann:2002. This is because VaR is not a coherent risk measure; it does not satisfy the property of subadditivity, meaning that the VaR of a combined portfolio can exceed the sum of the VaRs of its individual components. As a result, VaR may inadequately reflect the benefits of diversification.

In contrast, ES is a coherent risk measure, addressing these shortcomings by estimating the expected loss beyond the VaR threshold and, in particular, satisfying subadditivity Artzner-Delbaen-Eber-Heath:1999, Acerbi-Tasche:2002. Conditional on information $\mathcal{I}_{t-1}$, ES for period $t$ at probability level $\alpha$ is defined as

providing the average loss when $r_t$ falls below the VaR threshold in ((ref)). This makes ES a more theoretically sound measure that incentivizes diversification and better accounts for extreme risks. Recognizing these advantages, the {Basel Committee}, through its {Basel III} framework and the {Fundamental Review of the Trading Book (FRTB)}, has replaced VaR with ES as the primary metric for market risk capital requirements Basel:2019. This regulatory shift underscores the need for robust ES forecasting models that can accurately quantify tail risk under changing market conditions.

ES forecasts can often be derived as a natural extension of many VaR forecasting methods. Nonparametric approaches, including historical simulation and kernel density estimation, generate density forecasts from which both VaR and ES predictions can be obtained. Similarly, parametric methods, typically based on models for the conditional variance, such as GARCH specifications combined with a precise assumption about the conditional return distribution, provide concurrent forecasts for both risk measures.

A prominent semiparametric method for VaR forecasting is quantile regression, specifically the {Conditional Autoregressive VaR ({CAViaR})} models of Engle-Manganelli:2004, which model the conditional quantile directly without relying on specific distributional assumptions. Let $y_t = r_t - \mu_t$, where $\mu_t$ represents the conditional location of returns, often assumed to be zero, a constant, or a dynamic process such as AR(1), as in Kuester-Mittnik-Paolella:2006. A first-order symmetric absolute value (SAV) {CAViaR} specification is expressed as

with the VaR in ((ref)) given by $Q_{r_t}(\alpha \,\vert \, \mathcal{I}_{t-1}) = \mu_t + Q_{y_t}(\alpha \,\vert \, \mathcal{I}_{t-1})$. This model treats positive and negative deviations of $y_{t-1}$ symmetrically, relying solely on $\lvert y_{t-1} \rvert$ and thereby implicitly responding to volatility symmetrically. Alternatively, an asymmetric slope (AS) {CAViaR} specification is given by

where $\mathds{1}\{\cdot\}$ denotes the indicator function. This formulation allows the model to capture asymmetries in the conditional quantile dynamics by adjusting the slope based on whether the previous centered return was positive or negative, thus making it sensitive to directional changes in past returns. The {CAViaR} models in ((ref)) and ((ref)) allow the dynamics of $Q_{y_{t}}(\alpha \,\vert\,\mathcal{I}_{t-1})$ to vary across different probability levels $\alpha $ and have shown strong performance in empirical studies of VaR forecast accuracy Sener:2012. However, while {CAViaR} models excel in predicting VaR, they do not inherently provide a mechanism for generating ES forecasts, as they focus solely on a particular quantile of the return distribution.

Recent advances in risk forecasting have focused on methods that jointly estimate VaR and ES, leveraging the fact that these two risk measures are jointly elicitable Fissler-Ziegel:2016. A risk measure is said to be elicitable if it can be uniquely identified as the minimizer of a well-defined expected loss function, also known as a scoring function. In simpler terms, an elicitable risk measure can be evaluated and compared using a loss function that ranks different forecasting methods based on their accuracy. While ES is not elicitable on its own, the joint elicitability of VaR and ES allows for their simultaneous estimation and evaluation using a common loss function. This property has been pivotal in recent developments in financial risk forecasting.

Taylor:2019 proposes a joint VaR-ES model based on the asymmetric {Laplace} (AL) distribution, where the VaR component is modeled using a {CAViaR} specification. This approach uses the AL density to construct a likelihood function, whose maximization is equivalent to minimizing a strictly consistent Fissler-Ziegel:2016 scoring function. Similarly, Patton-Ziegel-Chen:2019 introduce a semiparametric approach within the {Generalized Autoregressive Score} (GAS) framework of Creal-Koopman-Lucas:2013, which also exploits the Fissler-Ziegel:2016 class of loss functions. Their GAS model dynamically adjusts both VaR and ES estimates based on past observations, specifically reacting to VaR violations while reverting to a long-run mean in the absence of breaches.

In this paper, we introduce a novel semiparametric {CAViaR}-based framework for forecasting VaR and ES, which we term the {Quantile-based Scale Dynamics} (QbSD) approach. The QbSD framework extends the widely used GARCH class by replacing the conditional standard deviation with a conditional scale---a measure analogous to the interquartile range (IQR), but more flexible as it is based on the difference between two specified quantiles of the return distribution, similar to Taylor:2005. To model the dynamics of the conditional scale, QbSD employs global {CAViaR} specifications in which persistence parameters are shared across quantiles to ensure parsimony, while local parameters account for level-specific adjustments. Parameter constraints enforce non-crossing quantile ordering and strictly positive scale values, ensuring internal consistency and robust tail risk estimation.

Beyond its theoretical advantages, the QbSD approach aligns with practical needs in financial risk management. Financial institutions rely on multiple model validations in modern risk management frameworks, particularly for stress-testing and regulatory reporting under Basel III guidelines. Given its flexibility and lack of reliance on specific distributional assumptions, QbSD can serve as a complementary tool for financial institutions in benchmarking VaR and ES estimates against conventional models. This is particularly relevant in periods of financial distress, where traditional parametric models may struggle to adapt to evolving risk dynamics. Furthermore, by providing more accurate tail risk estimates, QbSD can aid institutions in regulatory capital calculations and internal risk monitoring, supporting more resilient decision-making.

To achieve this, QbSD models the dynamics of conditional quantiles, enabling a structured estimation of both VaR and ES. Specifically, VaR is estimated from the left-tail quantile of the rescaled returns, while ES is approximated by averaging quantiles across probability levels below the VaR threshold. This approach leverages the representation of ES as the integral of conditional quantiles across the left tail of the distribution up to the VaR probability level Acerbi-Tasche:2002. As a result, the method delivers valid, well-ordered risk estimates across all probability levels, including those in the extreme tails of the distribution.

This approach shares similarities with the Filtered Historical Simulation (FHS) method described by Christoffersen:2012, which first applies a volatility model---typically a GARCH specification with an assumed innovation distribution---to filter historical returns. FHS then standardizes the returns by dividing them by the conditional volatility estimates, and subsequently computes the quantiles of these residuals to estimate VaR. The average residual loss beyond the VaR estimate is used to calculate ES. In contrast, the semiparametric QbSD approach does not rely on a specific distributional assumption. Instead, it models the conditional scale dynamics using restricted {CAViaR} specifications, offering a robust and distribution-free framework for forecasting VaR and ES.

The remainder of the paper is organized as follows. Section 2 develops the QbSD approach for forecasting VaR and ES. Section 3 presents simulation experiments designed to evaluate the performance of the proposed method across a range of market conditions. In Section 4, we assess the QbSD method's effectiveness through an empirical application, using daily returns from major international stock indices, including the volatile period surrounding the COVID-19 pandemic. Finally, Section 5 concludes the paper. Additional numerical and empirical results are provided in the Supplementary material.

Quantile-based modeling

Background

Consider first the unconditional distribution of asset returns and suppose it belongs to the location-scale family, so that

where $\mu$ is the location parameter, $s>0$ is the scale parameter, and $\varepsilon_{t}$ is an i.i.d. innovation with distribution function $F_{\varepsilon}$ and zero location.

There are several ways to define the scale parameter. A commonly used one is the standard deviation, which is based on moments. Despite its widespread usage, the standard deviation is sensitive to the presence of outliers and may even be undefined or infinite for distributions with heavy tails. This sensitivity arises from the fact that the standard deviation has a zero breakdown point, meaning that even a single outlier can significantly distort the estimate Huber-Ronchetti:2009. In contrast, scale estimators based on interquantile ranges have non-zero breakdown points, making them more robust to outliers and extreme values, which are common in financial return data.\footnote{Kim-White:2004 also advocate using quantiles to estimate characteristics of return distributions, focusing on robust measures for capturing skewness and kurtosis. Through simulations, they demonstrate that empirical moment-based estimators are significantly more sensitive to outliers. See also White-Kim-Manganelli:2010 and Ghysels:2014. }

With the location-scale representation in ((ref)), the $p$th quantile of the unconditional return distribution is $Q_{r_t}(p)=\mu+sF_{\varepsilon}^{^{-1}}(p)$. Let $0<p<0.5$ be a chosen quantile level. Following Taylor:2005, consider the quantile at the $(1-p)$th level as a natural and compelling counterpart. Even when $\mu$ is unknown, we can obtain $s$ as

where $c_p=F_{\varepsilon}^{-1}(1-p)-F_{\varepsilon}^{-1}(p)$ is a scaling term. Note that when $r_{t}\sim N(\mu,\sigma^{2})$, then ((ref)) yields $s=\sigma$. The quantity $\text{IQR}=Q_{r_t}(0.75)-Q_{r_t}(0.25)$ is the well-known interquartile range, which takes the form in ((ref)) by setting $c_p=1$. Another measure that fits this format is the PearsonTukey1965 scale measure that uses $p=0.05$ and $c_p=q_{0}(1-p)-q_{0}(p)$, where $q_{0}(p)$ is the $p$th quantile of the standard normal distribution. See Taylor:2005 and KotzvanDorp2005 for other examples.

Conditional dynamics

We extend the basic specification in ((ref)) to allow the location and scale parameters to be conditionally dynamic so that

where $\mu_{t}$ and $s_{t}>0$ are functions of past information, $\mathcal{I}_{t-1}$. This structure mirrors the typical GARCH framework, where both the conditional mean and the conditional standard deviation evolve over time based on past information. In the following, we focus primarily on first-order specifications when modeling the dynamics in $\mu_{t}$ and $s_{t}$. While higher lag orders can be taken into account if required, such specifications are widely favored and commonly used in practice.

In some applications, it may be reasonable to set $\mu_{t}$ to zero or to assume that $\mu_{t} = \mu$, as done in Taylor:2019. However, autocorrelation in financial returns is often non-negligible Kuester-Mittnik-Paolella:2006. This property can be incorporated into our quantile-based modeling framework by allowing the returns to have a time-varying conditional median of the form $\mu_{t}=Q_{r_{t}}(0.5\,\vert\,\mathcal{I}_{t-1})$, which may be captured by a quantile autoregression (QAR) model, as proposed by KoenkerXiao2006. For example, a first-order QAR model is formulated as

where $\mu$ represents a baseline median, and $\phi$ captures the impact of $r_{t-1}$ on the current period's median.

Our analysis concentrates on downside market risk. Accordingly, we define $0 < \tau \leq 0.05$, where $\tau$ represents the typical left-tail region of interest in financial risk management, capturing extreme losses. We then set $\tau \leq p < 0.5$, where $p$ specifies the quantile levels used to construct the scale measure. Consider the conditional version of ((ref)) written here as

where $c_p=F_{\varepsilon}^{-1}(1-p)-F_{\varepsilon}^{-1}(p)$ is a time-invariant scaling term, as before.\footnote{As Taylor:2005 observes, estimating scale (or volatility) using the interval between tail quantiles shares similarities with range-based volatility estimation Parkinson:1980, Garman-Klass:1980, Alizadeh-Brandt-Diebold:2002. This class of methods relies on the difference between the highest and lowest log prices, which essentially correspond to the quantiles at $1-p = 1$ and $p = 0$, respectively. } It is interesting to note here again that $s_{t}$ is invariant to the specific process governing $\mu_{t}$ in ((ref)). We can then express the model specification as

where $y_t = r_t - \mu_t$ and $\varepsilon_{t}^{\ast}= y_t/s_{t}^{\ast}=\varepsilon_{t}/c_p$ is simply a rescaled version of $\varepsilon_{t}$. The $\tau$th conditional quantile of $y_{t}= s_{t}^{\ast}\varepsilon_{t}^{\ast}$ is given by $Q_{y_{t}}(\tau \,|\, \mathcal{I}_{t-1}) = s_{t}^{\ast} F^{-1}_{\varepsilon^{\ast}}(\tau) $, where $F^{-1}_{\varepsilon^{\ast}}(\tau)$ is a constant. Similar to the conditional standard deviation in a GARCH model, $s_{t}^{\ast}$ serves as a common factor across quantiles, ensuring coherence and building dependence between them in a parsimonious manner.

To model the dynamics of $s^{\ast}_t$, we draw a parallel with the absolute value GARCH model of Taylor:1986 and Schwert:1990, which captures the conditional standard deviation using past absolute deviations. We assume

so that the scale of returns adjusts dynamically in a manner similar to the absolute value GARCH framework. Observing that $ s_{t-j}^{\ast} F^{-1}_{\varepsilon^{\ast}}(k) = Q_{y_{t-j}}(k \,|\, \mathcal{I}_{t-j-1}) $, we adopt the approach of Xiao-Koenker:2009, who use an absolute value GARCH structure, to express the two key quantiles in ((ref)) with the following global SAV {CAViaR} structure:

for both $k=p$ and $k=1-p$, subject to $\omega(p) < \omega(1-p)$, $\beta \ge 0$, and $\gamma \ge 0$. These restrictions ensure non-crossing quantiles and guarantee $s^{\ast}_t > 0 $ in ((ref)), provided the recursion begins with $Q_{y_{1}}(p \,|\, \mathcal I_{0}) < Q_{y_{1}}(1-p \,|\, \mathcal I_{0})$. Here, the persistence parameters $\beta$ and $\gamma$ are global, meaning they do not vary with $k$, while $\omega(k)= \omega F^{-1}_{\varepsilon^{\ast}}(k)$ is local and dependent on the quantile level $k$.

As an alternative to ((ref)), we also consider a threshold-based specification, similar to the threshold GARCH model of Zakoian:1994: \[ s^{\ast}_t = \omega + \beta s^{\ast}_{t-1} + \big( \gamma_{+} \mathds{1}\{y_{t-1} > 0\} + \gamma_{-} \mathds{1}\{y_{t-1} \leq 0\} \big) |y_{t-1}|, \] which allows $s^{\ast}_t$ to respond differently to positive and negative values of $y_{t-1}$. This leads to the following global AS {CAViaR} structure:

for both $k=p$ and $k=1-p$, subject to $\omega(p) < \omega(1-p)$, $\beta \ge 0$, $\gamma_+ \ge 0$, and $\gamma_- \ge 0$. Here, $\beta$, $\gamma_{+}$, and $\gamma_{-}$ are global parameters. This specification allows an asymmetric response in the conditional quantile process, capturing the differing impacts of positive and negative past shocks. Empirical studies by Nelson:1991, Glosten-Jagannathan-Runkle:1993, and Engle-Ng:1992, among others, provide strong evidence that asymmetry is crucial for accurately modeling volatility dynamics, as negative stock returns are typically associated with greater increases in return volatility. This justifies the inclusion of an asymmetric term in financial time series models such as the AS {CAViaR}. For simplicity, the notation in ((ref)), ((ref)), ((ref)), ((ref)), and ((ref)) omits the dependence on the involved parameters.

For estimation purposes, we recommend the following sequential strategy: (i) estimate the parameters of ((ref)) and obtain the fitted values $\hat{\mu}_{t}=\hat{Q}_{r_t}(0.5\,\vert\,\mathcal{I}_{t-1})$, and then (ii) estimate the parameters of the {CAViaR} specifications with $\mu_{t}$ replaced by $\hat{\mu}_{t}$. While this approach may not be as statistically efficient as joint estimation, it offers the advantage of being simpler and more numerically robust. We begin with the estimation of the QAR model. Let ${\boldsymbol{\theta}}_{k}$ represent the parameter vector for ((ref)) with $k=0.5$. The estimate $\hat{\boldsymbol{\theta}}_{k}$ can be obtained by solving the minimization problem

where $\rho_{k}[u]=u(k-\mathds{1}\{u \leq 0\})$ is the quantile regression “check" loss function KoenkerBassett1978.

The estimation of ((ref)) and ((ref)) proceeds with restricted quantile regression. When fitting these {CAViaR} specifications to data, a choice needs to be made for the initial value $Q_{y_{1}}(k \,|\, \mathcal I_{0}) $. A simple and reasonable approach is to set $Q_{y_{1}}(k \,|\, \mathcal I_{0}) = \hat Q_{y}(k)$, the corresponding sample quantile. This choice ensures that all subsequent conditional quantiles will satisfy the non-crossing condition, both in and out of sample. For instance, in the case of ((ref)), the parameter estimates for $\boldsymbol{\theta} = \big( \omega(p), \omega(1-p), \beta, \gamma \big)'$ are found by solving:

with $Q_{y_{1}}(k \,|\, \mathcal I_{0}) = \hat Q_{y}(k )$, for $k \in \{p, 1-p \}$. Note that the constraints apply only to the model parameters, without involving the data. The estimation of ((ref)) follows a similar procedure, with the constraint on $\gamma$ replaced by $\gamma_+ \ge 0$ and $\gamma_- \ge 0$.

With the estimated parameters, the fitted value of ${s}_{t}^{\ast}= \hat{s}_{t}^{\ast}(p)$ in ((ref)) is calculated as \[ \hat{s}_{t}^{\ast}(p) = \hat{Q}_{y_{t}}(1-p\,|\,\mathcal{I}_{t-1}) - \hat{Q}_{y_{t}}(p\,|\,\mathcal{I}_{t-1}). \] Then, using the fitted location $\hat{\mu}_{t}$ and scale numerator $\hat{s}_{t}^{\ast}(p)$, it is natural to use the unconditional quantiles of $\hat{\varepsilon}_{t}^{\ast}(p) = (r_{t} - \hat{\mu}_{t}) / \hat{s}_{t}^{\ast}(p)$ as estimates for $F_{\varepsilon^{\ast}(p)}^{-1}(\tau)$, since we are not assuming any particular distribution for $\varepsilon_{t}$ in ((ref)). The estimated left-tail $\tau$-quantile of the conditional return distribution is then given by

for $0 < \tau \leq 0.05$.

Estimating VaR

We estimate the VaR in ((ref)) based on the left-tail quantile of the conditional return distribution. In our quantile-based framework, the VaR estimate for a given probability level $\alpha$ is defined as

where $\hat{Q}_{r_{t}}(\alpha, p \,|\, \mathcal{I}_{t-1})$ is derived from the dynamic quantile estimate in ((ref)), evaluated at $\tau=\alpha$.

While the VaR estimate in ((ref)) can be computed for a given value of $p$, relying on a single scale-defining quantile level may lead to inefficiencies. Since the seminal work of PearsonTukey1965, it has been well understood that using the interval between symmetric tail quantiles to estimate standard deviation is influenced by the skewness and kurtosis of the distribution. Furthermore, as early as Mosteller:1946, it was suggested that interquantile range-based estimates of standard deviation could be improved by considering multiple quantile pairs. Building on these insights, we extend this principle to VaR estimation by computing estimates across multiple quantile levels and aggregating the results using an averaging scheme.

Specifically, we define a set of scale-defining quantile levels $\mathcal{P}$ and compute VaR for each $p \in \mathcal{P}$. The final VaR estimate is then given by

where $|\mathcal{P}|$ denotes the number of quantile levels used. For our implementation, we use the set $\mathcal{P} = \{ 0.05, 0.10, 0.15, 0.20, 0.25 \}$. While this choice is somewhat ad hoc, it corresponds to the relatively more robust values reported in Huber-Ronchetti:2009 for estimating scale using interquantile ranges.

We also considered alternative methods, such as using a single quantile level or aggregating via the median across $\mathcal{P}$. However, simulations indicate that averaging yields greater stability and lower overall estimation error. This approach mitigates fluctuations in tail behavior and enhances out-of-sample performance, particularly for heavy-tailed distributions. A detailed comparison of these methods is provided in Section A of the Supplementary material.

Approximating ES

The ES in ((ref)) represents the conditional expectation of loss, given that the loss has exceeded the VaR threshold at a specified probability level $\alpha$. As shown by Acerbi-Tasche:2002, among others, ES can be expressed as

where ${Q}_{r_{t}}(\tau \,\vert\, \mathcal{I}_{t-1})$ represents the conditional quantile for $\tau \leq \alpha$.

To compute this integral in practice, we approximate it using a sum over $N$ equally spaced subdivisions between 0 and $\alpha$, similar to a Riemann sum approximation: \[ \widetilde{\text{ES}}_{r_t}(\alpha, p, N) = \frac{1}{N} \sum_{i=1}^{N} \hat{Q}_{r_{t}}(\tau_{i}, p \,\vert\, \mathcal{I}_{t-1}), \] where $\tau_{i} = i \alpha / N$ for $i = 1, \ldots, N$, and each $\hat{Q}_{r_{t}}(\tau_{i}, p \,\vert\, \mathcal{I}_{t-1})$ is the corresponding quantile estimate from ((ref)).

As with VaR, we compute ES for several values of $p$ to obtain a more robust estimate. Specifically, for each fixed $N$, we take the average ES across different quantile levels $p$: \[ \widehat{\text{ES}}_{r_t}(\alpha, N) = \frac{1}{|\mathcal{P}|} \sum_{p \in \mathcal{P}} \widetilde{\text{ES}}_{r_t}(\alpha, p, N). \] This averaging process reduces estimation errors associated with any single $p$, leading to a more efficient and reliable ES estimate. The order of operations is important: by first approximating the integral over $\tau$ before averaging across $p$, the quantiles used in each $\widetilde{\text{ES}}_{r_t}(\alpha, p, N)$ share a consistent scale structure. Averaging across $p$ only afterward preserves the internal coherence of each quantile specification and avoids mixing quantiles constructed under differing scale definitions.

Rather than fixing $N$ a priori, we employ an adaptive selection method. We initialize with $N=4$ and iteratively increase it, computing ES at each step. The iteration stops when successive averaged ES estimates stabilize within a predefined tolerance level $\epsilon$: \[ \big\vert \widehat{\text{ES}}_{r_t}(\alpha, N) - \widehat{\text{ES}}_{r_t}(\alpha, N-1) \big\vert < \epsilon, \] where we set $\epsilon = 0.0001$ in our implementation. This ensures sufficient numerical accuracy while avoiding unnecessary computations. The stabilized ES estimate is denoted $\widehat{\text{ES}}_{r_t}(\alpha)$.

An important advantage of the QbSD approach is its ability to prevent inconsistencies common in tail risk estimation. Specifically, it avoids: (i) quantile crossings, (ii) crossings among ES values at different probability levels, and (iii) crossings between quantiles and their corresponding ES values. This method ensures valid, well-ordered risk estimates across all probability levels $\alpha$, even in the extreme tails of the distribution.

Simulation experiments

In this section, we investigate the accuracy of the proposed QbSD approach with respect to the one-step-ahead VaR and ES forecasts at probability level $\alpha$. The performance of the QbSD models is evaluated using mean absolute error (MAE) and root mean squared error (RMSE), and compared to several benchmark models, including GARCH and joint VaR-ES models. The analysis also assesses the models' ability to capture skewness, heavy tails, and leverage effects in financial returns.

Simulation design

The data-generating process (DGP) is the asymmetric power ARCH (APARCH) model of Ding-Granger-Engle:1993, specified as \begingroup {4pt} {7pt}

\endgroup where $\omega, \beta, \gamma, \delta > 0$ and $-1 < \theta < 1$. The innovations $\varepsilon_t$ follow the skewed $t$-distribution of Hansen1994. This model nests several well-known specifications as special cases. In particular, when $\delta = 2$ and $\theta = 0$, the model reduces to a standard GARCH model. When $\delta = 1$, it becomes the absolute value GARCH model of Taylor:1986 and Schwert:1990, which implies direct {CAViaR} representations for the conditional quantiles Xiao-Koenker:2009. Furthermore, when $\theta \neq 0$, the model implies an asymmetric response of volatility to past returns. When $\theta > 0$, this asymmetry takes the form of a leverage effect, whereby negative returns have a larger impact on future volatility than positive returns of the same magnitude Black:1976, Christie:1982, French-Schwert-Stambaugh:1987.

The density function of the innovation term $\varepsilon_t$ appearing in ((ref)) is defined as

where the constants $a$, $b$, and $c$ are themselves defined as \[ a = 4 \lambda c \frac{v-2}{v-1}, \quad b = \sqrt{1 + 3 \lambda^2 - a^2}, \quad c=\frac{\Gamma(\frac{v+1}{2})}{\sqrt{\pi (v-2)} \Gamma(\frac{v}{2})}, \] and the parameters $v > 2$ and $-1 < \lambda < 1$ represent the degrees of freedom and the asymmetry of the distribution, respectively. The skewed $t$-distribution has a mean of zero and unit variance. When $ \lambda = 0 $, it reduces to a standardized version of the traditional Student-$t$ distribution. A positive $ \lambda $ implies right skewness. Financial returns typically have a higher probability of large negative returns than large positive ones, corresponding to left skewness ($\lambda < 0$).

The distribution in ((ref)) is piecewise-defined around its mode at $x=-a/b$, which separates the left and right parts of the density and plays an important role in the derivation of VaR and ES. Jondeau-Rockinger2003 show that the associated cumulative distribution function (CDF) is defined by \[ F_{\varepsilon}(x;v,\lambda)=\left\{

\right. \] where $F_{\mathcal T}(x;v)$ is the CDF of the Student-$t$ distribution with $v$ degrees of freedom, and the quantile function of the skewed $t$-distribution is given by \[ F_{\varepsilon}^{-1}(u;v,\lambda)=\left\{

\right. \]

Under this DGP, both the VaR and ES of $r_t$ are proportional to $\sigma_{t}$. The true one-step-ahead VaR forecast is $ \text{VaR}_{r_{t+1}}(\alpha; v, \lambda) = \sigma_{t+1} F_{\varepsilon}^{-1}(\alpha;v,\lambda) $ and the true one-step-ahead ES forecast is $ \text{ES}_{r_{t+1}}(\alpha; v, \lambda) = \sigma_{t+1} \text{ES}_{\varepsilon}(\alpha; v, \lambda), $ where the ES for Hansen1994's (Hansen1994) skewed $t$-distribution is given by Patton-Ziegel-Chen:2019 as\footnote{The expression in Appendix B.4 of Patton-Ziegel-Chen:2019 contains some typos. It can be verified, either analytically or by simulation, that (ref) is the correct ES expression for Hansen’s skewed $t$-distribution. In the second branch of (ref), the parameters are flipped to $(v,-\lambda)$ (so the constants change to $a_2,b_2$ accordingly) when evaluating $\text{ES}^{\ast}_{\varepsilon}$. The prefactor $(1-\alpha)/\alpha$ converts that right-tail ES at level $1-\alpha$ back to the desired left-tail level $\alpha$.}

with \[ \text{ES}^{\ast}_{\varepsilon}(\alpha; v, \lambda) = \frac{\tilde \alpha}{\alpha}(1-\lambda) \Big(-\frac{a}{b} + \frac{1-\lambda}{b} \text{ES}_{\mathcal T}(\tilde \alpha; v) \Big) \] and \[ \tilde \alpha = F_{\varepsilon} \bigg( \frac{b}{1-\lambda}\Big( F_{\varepsilon}^{-1}(\alpha;v,\lambda) + \frac{a}{b} \Big) ;v,0 \bigg). \] Here $\text{ES}_{\mathcal T}(\alpha; v)$ is the ES for a Student-$t$ distribution, standardized to have unit variance, given by \[ \text{ES}_{\mathcal T}(\alpha; v) = -\frac{1}{\alpha} f_{\mathcal T}( q_{\alpha}; v) \Big( \frac{v + q_{\alpha}^2}{v-1} \Big) \sqrt{\frac{v-2}{v} }, \] where $f_{\mathcal T}( x; v)$ and $q_{\alpha} = F_{\mathcal T}^{-1}(\alpha; v)$ are, respectively, the density and $\alpha$-quantile of the Student-$t$ distribution with $v$ degrees of freedom BrodaPaolella2011, Dobrev-Nesmith-Oh:2017.

In the DGP outlined in ((ref)), we set the parameters as follows: $\omega = 0.05$, $\beta = 0.85$, $\gamma = 0.10$, and $\delta = 1.5$, while allowing $\theta$ to take values in $\{0, 0.5\}$. For the skewed $t$-distribution in ((ref)), we select $v \in \{5, 20 \}$ to represent different tail thicknesses, with $v = 5$ corresponding to heavy tails and $v = 20$ representing thinner tails. We use $\lambda \in \{0, -0.5\}$ to capture different levels of skewness, with $\lambda = 0$ indicating no skewness and $\lambda = -0.5$ implying left skewness. The sample size is set to $T = 1250$, which corresponds to roughly five years of daily returns and matches the primary rolling-window length used in our empirical application. The results presented in the main text focus on this baseline case; additional simulations for $T \in \{250, 2500\}$ are reported in Section B of the Supplementary material. We consider three probability levels for the forecasts: $\alpha \in \{0.01, 0.025, 0.05\}$. The choice of $\delta = 1.5$ allows us to explore a setting where the power parameter lies between the standard GARCH model ($\delta = 2$) and the absolute value GARCH model ($\delta = 1$).\footnote{This value of $\delta$ is close to the estimate of $1.43$ reported by Ding-Granger-Engle:1993 in their study of S&P 500 returns.} We set $\theta \in \{0, 0.5 \}$ to explore the role of asymmetry in volatility dynamics. As noted earlier, $\theta = 0$ yields a symmetric response, while $\theta = 0.5$ introduces a leverage effect.

For each DGP configuration, we evaluate the forecasting accuracy by computing the MAE and the RMSE across all 1,000 generated return series. These metrics allow us to compare the performance of the models in capturing the true one-step-ahead VaR and ES forecasts. In our simulations, we consider the QbSD model with the two global {CAViaR} specifications in ((ref)) and ((ref)). Consistent with the simulation DGP in ((ref)), all QbSD and benchmark models, presented next, assume a zero location for returns. This assumption simplifies the modeling process, allowing us to focus on volatility and tail risk dynamics.

Benchmark models

To assess the performance of our QbSD models, we compare them with several established models commonly used in financial risk forecasting. The benchmarks include the VaR and ES forecasts from the standard GARCH model, the asymmetric GJR-GARCH model Glosten-Jagannathan-Runkle:1993, and the EGARCH model Nelson:1991, each with normal innovations. Additionally, for the GARCH and GJR-GARCH models, we consider alternative distributional assumptions for the innovations, including the Student-$t$ distribution and Hansen1994's (Hansen1994) more general skewed $t$-distribution.\footnote{We consider the EGARCH model only with normal innovations, as it tends to exhibit instability when applied with $t$-distributed errors Nelson:1991.}

We also include the joint VaR-ES model proposed by Taylor:2019, which employs a pseudo-maximum likelihood estimation approach based on the AL density function for the joint estimation of VaR and ES. The conditional AL density function is given by \[ f(r_{t}) = \frac{(\alpha-1)}{\text{ES}_{t}}\exp\left(\frac{\left(r_{t}-\text{VaR}_{t}\right)\left(\alpha- \mathds{1}\{ r_{t} \leq \text{VaR}_{t} \} \right)}{\alpha \text{ES}_{t}}\right), \] where $\text{VaR}_{t}$ is modeled using either the SAV or AS {CAViaR} specifications from expressions ((ref)) and ((ref)), with $ y_{t} $ replaced by $ r_{t} $ under the mean zero assumption. The complete model includes two alternative specifications for the dynamics of the ES:

or

where $ x_{t} $ follows a dynamic updating process: \[ x_{t} =

\] This formulation provides a smooth adjustment to the magnitude of exceedances beyond the VaR, dynamically updating based on past returns and exceedances.

To differentiate between the two AL-based specifications, we use the notation “$\text{AL}_{\text{Mult.}}$” for the multiplicative version in ((ref)), where ES is modeled as a multiple of VaR, and “$\text{AL}_{\text{AR}}$” for the autoregressive version in ((ref)). For these AL-based models, we further indicate the choice of VaR specification using the suffixes “-SAV” for the SAV {CAViaR} model in ((ref)) and “-AS” for the AS {CAViaR} model in ((ref)). For example, $\text{AL}_{\text{Mult.}}$-SAV refers to the multiplicative AL-based model with the SAV {CAViaR} specification. The QbSD models, in contrast, rely on global {CAViaR} specifications, where the scale dynamics are shared across quantiles. Here, the suffix “-gSAV” denotes the global SAV {CAViaR} specification ((ref)), while “-gAS” refers to the global AS {CAViaR} specification ((ref)). For example, QbSD-gSAV refers to the QbSD model with the global SAV {CAViaR} specification.

In addition to the quantile-based models, we also consider the one-factor GAS model proposed by Patton-Ziegel-Chen:2019, which provides an alternative framework for jointly modeling VaR and ES. This model is based on the following member of the Fissler-Ziegel:2016 class of loss functions: \[ L_{\text{FZ0}} (r_{t}, \text{VaR}_{t}, \text{ES}_{t} ) = -\frac{1}{\alpha \text{ES}_{t}} \mathds{1}\{ r_{t} \leq \text{VaR}_{t}\} (\text{VaR}_{t} - r_{t}) + \frac{\text{VaR}_{t}}{\text{ES}_{t}} + \log(- \text{ES}_{t}) - 1. \] The $\text{FZ}_0$ loss function is strictly consistent for the pair (VaR, ES). Although the loss itself is not homogeneous of degree zero---due to the $\log(-\text{ES}_t)$ term---the differences in loss values between competing forecasts are. This ensures that forecast rankings based on average loss are invariant to positive scaling of both returns and risk measures.

The joint specification for VaR and ES is \[ \text{VaR}_{t} = \zeta \exp(\kappa_{t}), \quad \text{ES}_{t} = \xi \exp(\kappa_{t}), \quad \text{with} \ \xi < \zeta < 0, \] where the common factor $ \kappa_{t} $ evolves dynamically as \[ \kappa_{t} = \omega + \beta \kappa_{t-1} + \gamma g_{t-1}, \] and the score function $ g_{t} $, derived from the FZ0 loss function, is given by \[ g_{t} = -\frac{1}{\text{ES}_{t}} \left( \frac{1}{\alpha} \mathds{1}\{r_{t} \leq \text{VaR}_t \} r_{t} - \text{ES}_{t} \right). \] This approach dynamically updates both VaR and ES based on the information in past losses.\footnote{Since the Hessian of the FZ0 loss is constant, the scaling matrix in the GAS framework of Creal-Koopman-Lucas:2013 can be set to one without loss of generality. This simplifies estimation while preserving the model's essential dynamics.} Following Patton-Ziegel-Chen:2019, we set $\omega = 0$ for identification and estimate the model by minimizing the average FZ0 loss.

These competing models offer a range of benchmarks for evaluating the forecasting performance of the QbSD models. In the next section, we present the simulation results and assess how well each model captures key features of financial returns across the various DGP configurations. It is interesting to note that none of the models exactly matches the DGP in ((ref)), which makes the forecasting task more challenging and offers a stricter test of adaptability.

Numerical results

This section presents the forecasting results for the QbSD and benchmark models across the simulation scenarios. Tables 1--4 report out-of-sample one-step-ahead VaR and ES accuracy (MAE and RMSE) under symmetric and asymmetric volatility dynamics.

Table 1 reports VaR forecast accuracy with no leverage (symmetric volatility dynamics). Parametric GARCH-family models with correctly specified innovations (Student-$t$ under symmetry; skew-$t$ under left skew) lead across configurations. With symmetric innovations ($v=20,\lambda=0$), GARCH/GJR with Student-$t$ is best across $\alpha$ (e.g., MAE $0.071,0.054,0.044$; RMSE $0.095,0.073,0.059$), while QbSD-gSAV is competitive but not top. With left skewness ($v=20,\lambda=-0.5$), skew-$t$ dominates (e.g., GJR-skew-$t$ MAE $0.098,0.068,0.049$; RMSE $0.132,0.093,0.069$), with only minor differences between GARCH and GJR. Under heavier tails ($v=5$), the pattern persists: when $\lambda=0$, GARCH/GJR-$t$ remains best (MAE $0.103,0.068,0.050$; RMSE $0.143,0.097,0.072$); when $\lambda=-0.5$, GJR-GARCH with skew-$t$ leads (MAE $0.145,0.084,0.054$; RMSE $0.213,0.134,0.093$). Across all blocks, QbSD-gSAV improves on AL-based and GAS alternatives but trails the best-specified GARCH variants.

Turning to ES under symmetric volatility (Table 2), the same ordering emerges: matching the innovation distribution to the DGP is pivotal. In the symmetric case ($v=20,\lambda=0$), GARCH/GJR-$t$ attains the lowest errors (MAE $0.100,0.075,0.060$; RMSE $0.129,0.099,0.081$). Introducing left skewness ($v=20,\lambda=-0.5$) shifts the frontier to skew-$t$ (GJR-skew-$t$ MAE $0.146,0.105,0.080$; RMSE $0.191,0.140,0.108$). Under heavier tails ($v=5$), Student-$t$ remains best when $\lambda=0$ (MAE $0.180,0.118,0.086$; RMSE $0.241,0.162,0.120$), while skew-$t$ variants lead for $\lambda=-0.5$ (e.g., GJR-skew-$t$ MAE $0.117,0.173,0.117$; RMSE $0.174,0.247,0.383$). AL-based and GAS models are consistently less accurate.

Table 3 reports VaR accuracy with leverage ($\theta=0.5$). QbSD-gAS delivers the lowest errors in most configurations featuring skewness and/or heavy tails. For example, at $v=20,\lambda=-0.5$, it is best across $\alpha$ (MAE $0.172,0.125,0.096$; RMSE $0.238,0.183,0.138$). In the heavy-tailed but symmetric case ($v=5,\lambda=0$), QbSD-gAS is generally strongest, with EGARCH edging it at $\alpha=2.5\%$ (MAE $0.105$ vs.\ $0.110$). Under heavy tails with left skewness ($v=5,\lambda=-0.5$), GJR-skew-$t$ attains the lowest MAE ($0.233,0.152,0.108$), with QbSD-gAS close behind. A notable exception is the thin-tailed, symmetric DGP ($v=20,\lambda=0$), where EGARCH dominates (MAE $0.105,0.075,0.062$; RMSE $0.148,0.107,0.086$).

Table 4 presents ES accuracy with leverage. QbSD-gAS leads in most configurations. It is best across $\alpha$ at $v=20,\lambda=-0.5$ (MAE $0.235,0.174,0.140$; RMSE $0.314,0.238,0.195$) and at $v=5,\lambda=0$ (MAE $0.292,0.186,0.136$). In the heavy-tailed, left-skewed case $v=5,\lambda=-0.5$, GJR-skew-$t$ has the edge (MAE $0.388,0.263,0.191$; RMSE $0.526,0.352,0.256$). For the thin-tailed symmetric DGP $v=20,\lambda=0$, EGARCH is competitive and surpasses QbSD-gAS at higher $\alpha$ (e.g., MAE $0.113,0.086$ and RMSE $0.157,0.122$ at $2.5\%$ and $5\%$, vs.\ QbSD-gAS $0.125,0.103$ and $0.165,0.137$); QbSD-gAS is slightly better at $1\%$ (MAE $0.165$ vs.\ $0.167$).

Overall, the evidence suggests a clear division: when volatility responds symmetrically and the innovation distribution is correctly specified (e.g., appropriately heavy-tailed or skewed), parametric GARCH variants set the benchmark for VaR and ES; once leverage/asymmetry matters, allowing asymmetric-slope (gAS) scale dynamics is pivotal---the proposed QbSD-gAS specification is typically strongest, particularly for ES, with EGARCH competitive and occasionally best in thin-tailed symmetric settings. Across these scenarios, AL-based and GAS models remain less accurate. We next examine whether these patterns carry over to international stock index returns.

Empirical Application

We evaluate one-day-ahead VaR and ES forecasts for the daily log returns of the S&P 500, DJIA, NASDAQ, EURO STOXX 50, FTSE 100, DAX, CAC 40, and TSX stock indices. The daily adjusted closing prices for these indices were obtained from Yahoo Finance, spanning October 4, 2002, to February 2, 2024, across indices. Due to differences in trading days and data availability, the number of observations initially varies across the series.

To ensure comparability across stock indices, we standardized the sample size to match the index with the shortest available data. Specifically, the EURO STOXX 50 has the shortest sample, starting on October 21, 2002, and containing 5{,}357 daily observations, ending on February 2, 2024. For consistency, we truncated the data for indices with longer samples by counting backward from February 2, 2024, to retain exactly 5{,}357 observations for each index. This adjustment implies the following starting dates: October 22, 2002, for the S&P 500, DJIA, and NASDAQ; October 21, 2002, for the EURO STOXX 50; November 18, 2002, for the FTSE 100; January 3, 2003, for the DAX; February 27, 2003, for the CAC 40; and October 4, 2002, for the TSX. As returns are calculated as first differences of log prices, this sample size yields 5{,}356 daily return observations for all indices.

Figure 1 illustrates the time series of daily log returns (in percentages) for the S&P 500, DJIA, NASDAQ, and EURO STOXX 50 indices, each exhibiting notable spikes in volatility during significant economic events such as the 2008 financial crisis and the COVID-19 pandemic. Figure 2 shows the daily log returns for the FTSE 100, DAX, CAC 40, and TSX indices, which similarly exhibit heightened volatility during periods of financial stress.

We use a rolling window of $R=1250$ estimation observations (approximately five trading years) to produce one-day-ahead forecasts for VaR and ES. This choice balances sampling variability against responsiveness to structural change. Model parameters are estimated using the most recent 1{,}250 observations, and forecasts are updated daily as the window moves through the sample period. Given the common sample of 5{,}356 daily returns, this setup leaves an out-of-sample evaluation period of 4{,}106 observations. We forecast VaR and ES at probability levels of 1%, 2.5%, and 5%. The 1% and 5% levels are widely used in the literature, whereas the 2.5% level corresponds to the regulatory standard for ES adopted under Basel III. We also examined alternative rolling-window sizes of $R=250$ and $R=2500$. To maintain focus in the main text, we present results for the 1{,}250-observation window, while these alternative specifications appear in Section C of the Supplementary material.

We evaluate a total of twenty-eight competing risk models for forecasting VaR and ES, including: (i) fourteen GARCH-type models (GARCH and GJR-GARCH with normal, Student-$t$, and Hansen1994 skew-$t$ innovations under zero or AR(1) conditional means, plus EGARCH with normal innovations under zero or AR(1) means), (ii) two GAS variants with zero and AR(1) means Patton-Ziegel-Chen:2019, (iii) eight AL-based models following Taylor:2019 (SAV and AS CAViaR for VaR dynamics, multiplicative and autoregressive specifications for ES dynamics, each under zero or AR(1) means), and (iv) four QbSD models (gSAV and gAS scale dynamics under zero or QAR(1) conditional medians).

In the following tables, we use the prefixes “AR-” and “QAR-” to denote an AR(1) conditional mean and a QAR(1) conditional median, respectively; models without a prefix assume a zero conditional location.

Evaluation of VaR forecasts

We evaluate the accuracy of VaR forecasts using the model confidence set (MCS) procedure of Hansen-Lunde-Nason:2011, which identifies a subset of models that are not significantly outperformed at a given confidence level. The MCS framework is particularly useful in a multi-model comparison setting, as it does not require selecting a single best model but instead acknowledges model uncertainty by retaining all models that cannot be statistically eliminated.

Following Hansen-Lunde-Nason:2011's (Hansen-Lunde-Nason:2011) suggestion, we apply the range-based procedure using the maximum absolute $t$-statistic, which iteratively removes the model with the largest loss difference relative to the others until only models that cannot be statistically distinguished from the best-performing ones remain. The final MCS consists of models that survive this elimination process at the chosen confidence level.\footnote{Specifically, we base the MCS on the $T_{R, \mathcal{M}}$ test statistic developed by Hansen-Lunde-Nason:2011 and compute it with 1{,}000 bootstrap iterations using the R package `MCS' Bernardi-Catania:2018.}

The test statistic is computed using the quantile score loss, given by

where $r_{t+1}$ is the realized return and $\text{VaR}_{t+1}$ is the predicted VaR at a given probability level $\alpha$. This loss function is asymmetric and penalizes underpredictions more heavily than overpredictions, making it well suited to the evaluation of downside risk.

Tables 5--7 present the MCS results for 1%, 2.5%, and 5% VaR forecasts across the eight international stock indices. Models excluded from the 90% MCS ($\widehat{\mathcal{M}}_{90\%}^{*}$) for a given index are left blank. The remaining models are ranked by their average quantile score $t$-statistics; the “\#” column counts the number of indices for which a model remains in the MCS, while the “Avg. rank” and “Final rank” summarize overall performance.

For the 1% tail (Table 5), QbSD-gAS attains the best overall rank and appears in the MCS for all eight indices; its closest competitors are GARCH variants with skew-$t$ innovations (GJR-GARCH-skew-$t$, GARCH-skew-$t$), followed by GJR-GARCH-$t$, while AL-based specifications generally trail. At 2.5% (Table 6), QbSD-gAS again ranks first overall (in the MCS for all eight indices), with AL-based models featuring AS dynamics (AL$_{\text{Mult.}}$-AS, AR-AL$_{\text{Mult.}}$-AS) close behind; among GARCH-type models, skew-$t$ remains the most competitive and EGARCH sits in the upper middle of the table. By 5% (Table 7), QbSD-gAS still leads, EGARCH moves into second place, AL$_{\text{Mult.}}$-AS and the QAR-driven QbSD-gAS (QAR-QbSD-gAS) remain prominent, and skew-$t$ GARCH variants continue to perform well.

Across probability levels, the symmetric QbSD specification (QbSD-gSAV) underperforms its asymmetric counterpart, underscoring the value of allowing an asymmetric scale response in practice. Within the GARCH family, skew-$t$ innovations generally improve performance at tighter tails (1% and 2.5%), whereas the strong showing of EGARCH at 5% suggests that volatility asymmetry can partly compensate when attention shifts away from the far tail.

While these results provide important insights into tail-risk quantile prediction, a fuller assessment requires joint evaluation of VaR and ES, which we consider next.

Joint evaluation of VaR and ES forecasts

To jointly evaluate VaR and ES forecasts, we use the AL log score proposed by Taylor:2019. This scoring rule is specifically designed to assess both VaR and ES, reflecting their joint elicitability. The AL log score is given by

where $\mu_{t+1}$ is the model's conditional mean and $\text{ES}_{t+1} < \mu_{t+1}$ is required. This proper scoring rule penalizes both the incidence and severity of forecast errors in the left tail, encouraging accurate joint predictions of VaR and ES.

Tables 8--10 summarize the MCS rankings under AL log scores. At the 1% tail (Table 8), the best overall performer is {AR-GJR-GARCH-skew-$t$}, with {QbSD-gAS} a close second; both appear in the MCS for all eight indices. Other skew-$t$ GARCH variants ({GARCH-skew-$t$}, {GJR-GARCH-skew-$t$}) are also highly competitive, whereas AL-based specifications sit further down the table. At 2.5% (Table 9), leadership shifts to {QbSD-gAS}, which ranks first overall and is included in the MCS for all eight indices. AL-based models with {asymmetric VaR dynamics} (AS CAViaR), namely {AL$_{\text{Mult.}}$-AS} and {AR-AL$_{\text{Mult.}}$-AS}, are close followers, with skew-$t$ GARCH entries forming the next tier. By 5% (Table 10), {QbSD-gAS} again takes the top spot (in the MCS for all indices), followed by {AL$_{\text{AR}}$-AS} and {AL$_{\text{Mult.}}$-AS}; the QAR-driven {QAR-QbSD-gAS} also ranks near the front. Skew-$t$ GARCH specifications remain competitive, while EGARCH occupies an upper-mid-table position.

Taken together, these results indicate that joint VaR-ES forecasting favors models that accommodate {asymmetric tail dynamics}. {QbSD-gAS} is consistently among the best across tail levels and leads at 2.5% and 5%, while skew-$t$ GARCH variants are particularly strong at the most extreme tail (1%). AL-based models with AS CAViaR VaR dynamics provide robust competitors at 2.5% and 5%, reinforcing the importance of allowing asymmetry in the VaR/ES system.

Conclusion

This paper introduced QbSD for forecasting VaR and ES: a semiparametric, distribution-free approach that models conditional scale via restricted quantile regression and accommodates skewness, heavy tails, and leverage---features that intensify during periods of market stress.

Across simulations and international equity indices, a consistent picture emerges. In simulation designs without leverage, parametric GARCH models with heavy-tailed or skewed innovation distributions set a high benchmark for both VaR and ES; Student-$t$ and skew-$t$ specifications are especially effective. With leverage, allowing asymmetric-slope quantile dynamics is pivotal: QbSD-gAS is frequently among the best, particularly for ES, while for VaR, EGARCH and skew-$t$ GJR-GARCH can also be highly competitive depending on the tail level and distributional setting.

In the empirical evaluation with a five-year rolling window ($R=1250$), QbSD-gAS is the clearest overall performer for VaR under quantile-score MCS: it attains the best average rank at 1%, 2.5%, and 5% across the eight indices, with skew-$t$ GARCH variants forming the next tier and EGARCH rising near the front at 5%. When VaR and ES are assessed jointly via AL log scores, performance at the most extreme tail (1%) tilts toward skew-$t$ GARCH (AR-GJR-GARCH-skew-$t$), with QbSD-gAS a close second; at 2.5% and 5%, QbSD-gAS leads, and AL-based specifications with AS CAViaR VaR dynamics (both multiplicative and autoregressive ES linkages) are robust followers.

Two practical lessons follow. First, asymmetry helps---but not all asymmetries are equal: QbSD-gAS delivers the most consistent improvements; skew-$t$ innovations help when empirical skewness is material; and EGARCH’s volatility asymmetry tends to help when return innovations are near-symmetric with relatively light tails, whereas GJR’s benefits appear mainly when combined with skew-$t$ under skewed, heavy-tailed conditions. Second, the preferred specification varies with the tail probability and the evaluation target: skew-$t$ GARCH may perform well at the far tail (1%) under joint VaR-ES scoring, whereas QbSD-gAS provides the most consistent gains across tail levels for VaR and remains among the leaders at 2.5% and 5% in the joint evaluation.

Overall, QbSD---especially with asymmetric-slope (gAS) quantile dynamics---offers a robust, flexible addition to the risk-management toolkit. Its distribution-free construction, strong showing under leverage and during turbulent periods, and straightforward extensibility to log-location-scale settings (e.g., strictly positive losses such as credit losses and insurance claims) make it well suited for practical deployment.

Acknowledgments

We thank the editor, {Dick van Dijk}, the associate editor, and an anonymous reviewer for their constructive feedback and insightful suggestions. {Richard Luger} is supported in part by funding from the {Social Sciences and Humanities Research Council of Canada}.

Data and code availability

The data and code needed to reproduce the results presented in this paper are available at {https://github.com/richardluger/QbSD}. The code is written in {R} with {C++} functions for computational efficiency.

\setcounter{page}{1}

Section A

Tables A1--A12 present comparisons of the mean absolute error (MAE) and root mean squared error (RMSE) for the QbSD forecasting approach, examining the impact of the scale-defining quantile level $p$ on the computation of VaR and ES. Specifically, we consider forecasts based on: (i) individual values of $p \in \{ 0.05, 0.10, 0.15, 0.20, 0.25 \}$; (ii) the mean over these $p$ values; and (iii) the median over them.

The results are obtained from 1,000 replications for each configuration of the asymmetric power ARCH model with innovations drawn from the skewed $t$-distribution, as specified in the main text. The lowest MAE and RMSE values in each table are highlighted in bold.

Overall, the results in Tables A1--A12 suggest that combining forecasts over multiple quantile levels $p$---either by taking the mean or, typically to a slightly lesser degree, the median---tends to yield lower forecast errors than relying on a single $p$ in isolation. This advantage of averaging is most noticeable when $T$ is relatively small (e.g., $T=250$; see especially Tables A1--A4), where individual-$p$ estimates exhibit higher variability. As $T$ increases to 1250 (Tables A5--A8) or 2500 (Tables A9--A12), the performance gap narrows, yet the multi-$p$ strategies remain competitive or superior under most DGP configurations.

Although occasional scenarios arise where an individual $p$ may match or slightly outperform the combined forecasts (e.g., $p=0.10$ in Table A7), no single $p$ exhibits consistent dominance across the wide range of parameter settings (leverage $\theta$, tail thickness $v$, skewness $\lambda$) or sample sizes $T$. Consequently, our findings recommend the “mean over $p$” method as a robust, practically simple strategy for balancing estimation accuracy and stability in quantile-based VaR and ES forecasting.

Section B

Tables B1--B8 report one-step-ahead VaR/ES forecast errors (MAE, RMSE) for the proposed QbSD models (gSAV, gAS) and benchmarks (skew-$t$ GARCH/GJR/EGARCH, AL-based joint VaR-ES, and GAS) across data-generating designs varying by $T\in\{250,2500\}$, $v\in\{5,20\}$, $\lambda\in\{0,-0.5\}$, $\theta\in\{0,0.5\}$, and $\alpha\in\{0.01,0.025,0.05\}$. Minima are bolded.

For $T=250$ without leverage (Tables B1--B2), conventional GARCH-type models dominate. With symmetric innovations, Normal/Student-$t$ GARCH/GJR achieve most minima; under left-skew, their skew-$t$ counterparts lead. QbSD models are competitive but rarely best, while AL-based estimators and GAS trail.

Introducing leverage at $T=250$ (Tables B3--B4) shifts the lead within the GARCH family: EGARCH frequently attains the lowest errors in symmetric cases, and \sloppy skew-$t$ GARCH/GJR often prevail under skewed or heavy-tailed DGPs. QbSD-gAS improves relative to QbSD-gSAV but seldom displaces the top GARCH variants at this sample size.

With larger samples and no leverage (Tables B5--B6, $T=2500$, $\theta=0$), GARCH/GJR remain the primary winners. Student-$t$ specifications account for many minima under symmetry and heavy tails, whereas skew-$t$ versions excel when $\lambda=-0.5$. QbSD narrows the gap but generally does not surpass the best GARCH-type models here.

With larger samples under leverage (Tables B7--B8, $T=2500$, $\theta=0.5$), QbSD-gAS delivers many of the lowest ES MAE/RMSE values across $\alpha$ and $(v,\lambda)$ settings, outperforming AL and GAS and often surpassing standard GARCH models. Two caveats remain: (i) under symmetric innovations, EGARCH can match or exceed QbSD-gAS for ES RMSE at some $\alpha$ (e.g., $v=20$, $\lambda=0$ at $2.5\%$ and $5\%$); and (ii) for VaR the best model can shift with tail level---EGARCH under symmetry and skew-$t$ GJR-GARCH/GARCH under skewed, heavy-tailed designs---while QbSD-gAS remains the most robust ES performer under leverage.

Overall, while GARCH-type models (especially skew-$t$ and EGARCH) are very strong at $T=250$ and in several $T=2500$ designs without leverage, QbSD-gAS becomes particularly effective for ES once leverage is present and the sample is large, delivering broad gains except in the most extreme skew/heavy-tail scenario. AL-based estimators and GAS underperform across configurations.

Section C

This appendix reports additional rolling-window evaluations with $R=250$ and $R=2500$. Tables C1--C12 present 90% Model Confidence Set (MCS) rankings for VaR and joint VaR-ES at the 1%, 2.5%, and 5% levels.

Across windows, indices, and tails, the quantile-based scale-dynamics class performs very strongly overall. QbSD-gAS is the front-runner with $R=250$: it ranks first in most VaR and joint VaR-ES exercises and is almost always in the top two. Its QAR extension (QAR-QbSD-gAS) typically follows closely-often finishing in the top three---but its long-window 1% VaR results slip on a few indices (see Table C7).

With the short window ($R=250$), patterns are clear. For VaR-only (quantile scores; Tables C1--C3), QbSD-gAS is the overall top performer at 1%, 2.5%, and 5%. For joint VaR-ES (AL log scores; Tables C4--C6), skewed-$t$ GARCH/GJR lead at 1% (Table C4) with QbSD-gAS a close second, while QbSD-gAS returns to the top at 2.5% and 5% (Tables C5--C6).

With the long window ($R=2500$), leadership varies by tail and score. For VaR (Tables C7--C9), AL-based specifications dominate at 1% (Table C7), QbSD-gAS ranks first at 2.5% (Table C8), and at 5% EGARCH narrowly takes the lead, with QbSD-gAS essentially tied on average but second in the final ranking (Table C9); note that both have the same average rank (1.9), with QbSD-gAS recording more \#1 finishes (4 of 8) while EGARCH is slightly more uniform across indices. For joint VaR-ES (Tables C10--C12), AL-based specifications lead at 1% and 2.5% (Tables C10--C11), whereas QbSD-gAS regains first place at 5% (Table C12).

Among benchmarks, skewed-$t$ GARCH/GJR are the strongest non-QbSD competitors and are frequently retained in the MCS, especially for $R=250$ and in VaR-only comparisons. Normal-based GARCH variants generally underperform on average, with two caveats: (i) EGARCH is a strong contender at 5% VaR for $R=2500$ (Table C9), and (ii) GJR-GARCH-normal is competitive at $R=250$ for 2.5% and 5% VaR (Tables C2--C3). Within the AL-based family, AS specifications consistently outperform their SAV counterparts. The gSAV variant of QbSD (QbSD-gSAV) typically trails QbSD-gAS, while QAR-QbSD-gSAV most often attains mid-table results.

Overall, these results align with the paper’s main finding: QbSD-gAS is frequently the strongest---especially with the short window---and remains consistently competitive across indices and tail levels under the long window; QAR-QbSD-gAS is often close behind. At the extreme tail and with long windows, AL-based variants tend to lead, and EGARCH is a notable runner-up at 5% VaR.