###### Abstract

In financial risk management, Value at Risk (VaR) is widely used to estimate
potential portfolio losses. VaR’s limitation is its inability to account for
the magnitude of losses beyond a certain threshold. Expected Shortfall (ES)
addresses this by providing the conditional expectation of such exceedances,
offering a more comprehensive measure of tail risk. Despite its benefits, ES
is not elicitable on its own, complicating its direct estimation. However,
joint elicitability with VaR allows for their combined estimation. Building on
this, we propose a new methodology named Conditional Autoregressive Expected
Shortfall (CAESar), inspired by the CAViaR model. CAESar handles dynamic
patterns flexibly and includes heteroskedastic effects for both VaR and ES,
with no distributional assumption on price returns. CAESar involves a three-
step process: estimating VaR via CAViaR regression, formulating ES in an
autoregressive manner, and jointly estimating VaR and ES while ensuring a
monotonicity constraint to avoid crossing quantiles. By employing various
backtesting procedures, we show the effectiveness of CAESar through extensive
simulations and empirical testing on daily financial data. Our results
demonstrate that CAESar outperforms existing regression methods in terms of
forecasting performance, making it a robust tool for financial risk
management.

Keywords: Autoregressive Models, Conditional Value at Risk, Expected
Shortfall, Risk Management

##  1 Introduction

Risk management is a crucial aspect of the financial literature and industry,
with Value at Risk (VaR) historically serving as the predominant measure. VaR,
a tail quantile of the conditional distribution of portfolio returns, has been
extensively studied and widely used by banks and other financial institutions
for setting regulatory capital requirements over the past two decades. Despite
its intuitive appeal, VaR lacks information regarding exceedances beyond the
quantile. Furthermore, it is not a coherent risk measure as it is not sub-
additive. To address this limitation, Expected Shortfall (ES) has emerged as a
complementary measure, providing the conditional expectation of these
exceedances. ES offers several theoretical advantages [2] as it is a coherent
risk measure [3]; in particular, in contrast to VaR, ES is subadditive,
meaning that ES for a portfolio cannot be larger than the sum of ES’s for the
single assets in the portfolio. Coherence makes ES a more comprehensive
measure of tail risk compared to VaR and, consequently, regulatory frameworks
like Basel IV are increasingly emphasizing it [19].

ES poses unique challenges, particularly because it is not elicitable on its
own, meaning it cannot be directly estimated as the minimum of a functional
[44, 24]. However, recent research has shown that VaR and ES are jointly
elicitable, enabling their combined estimation using suitable scoring
functions [20]. This joint elicitability facilitates the evaluation and
comparison of forecasting methods, a significant advancement in risk
management. Building upon this research, some recent works have proposed a
number of methods that can be grouped into static estimation, relying on Monte
Carlo simulations [40] or high-frequency data [15], and regressions, which can
be seen as general extensions of the Conditional Autoregressive Value-at-risk
by Regression quantiles (CAViaR) [18], in turn based on the general idea of
quantile/expectile regression [41, 12, 6]. Everything stems from a simple
intuition, which originates from the optimization of a jointly elicitable
function permitting the inference of the pair (VaR, ES). In a regression
context, price returns are observed as a time series on a daily time scale,
and the two risk measures are regressed upon such observations (or nonlinear
functions of these) and other covariates. Directly modeling VaR and ES avoids
the need for a distributional assumption on price returns, achieving a more
flexible description. This is key for financial time series, displaying fat
tails and skewness, indicating that extreme returns are more common than
predicted by a normal distribution and that the distribution may be
asymmetrical. Additionally, while returns exhibit zero autocorrelation,
volatility shows long-term memory and clustering; see [9] for a review.
Moreover, higher-order moments, like skewness and kurtosis, can be time-
varying as well, see e.g. [28, 45]. These facts may have a crucial impact on
the description of the tail of the distribution of returns, thus posing unique
challenges for the estimation of VaR and ES [4, 38]. For example, when
considering VaR, a Bayesian approach including nonlinear effects has been
proposed for the general quantile regression problem [23]. Regression methods
for conditional autoregressive Value at Risk (VaR) and Expected Shortfall
(ES), similar to the approach of Engle and Manganelli [18], provide a well-
suited and flexible framework for capturing time-varying patterns associated
with tail risk, thanks to the inclusion of heteroskedastic effects. In this
context, it is worth mentioning the recent work of Patton et al. [39], which
proposes a GAS framework based on the model by Creal et al. [10]. This
framework imposes a parametric structure for the dynamics similar to the
CAViaR model and defines a scoring rule that drives the evolution based on a
jointly elicitable loss function for VaR and ES. Two model specifications are
considered: one with a single factor and another with two factors, where the
factors are score-driven parameters determining the joint dynamics of VaR and
ES. The proposed Monte Carlo and empirical studies demonstrate the superior
performance of the one-factor model, suggesting co-movement dynamics for the
two tail risk measures. Within a similar context, Barrera et al. [5] propose a
general convergence analysis of a two-step approach to learn VaR and ES via
neural networks, minimizing a proper jointly elicitable function. Leveraging
machine learning permits one to naturally account for general nonlinear
patterns between regressors and tail risk. Finally, a natural approach for ES
estimation is based on its definition as a tail mean, suggesting the adoption
of a quantile regression approach for VaR estimation at different tail levels
and then obtaining a consistent estimation of the ES as the mean of the VaR
estimates, see [17, 33]. See the section below for a detailed review of these
works, which are benchmarks for our proposed method.

In this paper, we contribute to the literature on regression methods for the
joint estimation of the Value-at-Risk and the Expected Shortfall by proposing
a methodology, named Conditional Autoregressive Expected Shortfall (CAESar),
inspired by the CAViaR model that can handle general dynamic patterns in a
flexible way, without distributional assumptions on price returns. We include
heteroskedastic effects for VaR and ES, whose dynamics are driven by
regressors defined as nonlinear responses to the observed process. The
methodology is in line with the aim of Patton et al. [39]; in particular, we
consider the same jointly elicitable loss function in the estimation process.
However, our formulation is more general and less limiting in the description
of the tail risk by relaxing the constraint of the constant VaR-to-ES ratio of
the one-factor specification. Furthermore, our VaR and ES forecasts are
updated by always taking into account information about the price movement, in
contrast to the Patton et al. two factors model that only uses extreme events
-i.e., those breaking the quantile. CAESar parameters estimation is a three-
step approach. First, the conditional quantile VaR estimator for a given
probability level is obtained using the well-known CAViaR regression [18].
Second, we propose an autoregressive formulation for the ES based on quantile
information. Third, VaR and ES are jointly estimated by exploiting joint
elicitability. We impose the inequality $ES_{t}\leq VaR_{t}$ for all times
$t$, for a fixed probability level (the “monotonicity” constraint), using a
soft constraint on the loss function. In other words, a penalized estimation
process for CAESar is adopted to ensure “monotonicity” constraint, thus
avoiding the so-called ’crossing quantiles’ problem, a typical issue in
financial time series models[25]. This results in more realistic dynamics for
the VaR-ES pair.

Backtesting expected shortfall is a long-standing problem. In fact, non-
elicitability has long been considered a limitation for backtesting [8], only
recently surpassed, as pointed out in [1]. A number of testing procedures have
been proposed in recent years, see e.g. [13, 16], even if no general agreement
exists on the best approach. We corroborate our methodology on simulations and
daily financial data by using a plethora of testing procedures (reviewed in
detail below) recently proposed for backtesting ES, each one validating tail
risk estimates from different perspectives. Interestingly, after an extensive
evaluation and comparison with the state-of-the-art, we find a general
agreement in terms of forecasting performance ranking, with CAESar
outperforming any other regression method.

This paper is organized as follows. Section 2 reviews the related literature
and formalizes the problem. Section 3 introduces the CAESar methodology with a
discussion of the theoretical advantages with respect to previous works.
Section 4 reviews the testing procedures for ES backtesting adopted in the
experimental stage. Section 5 shows the Monte-Carlo and empirical results on
the stock indexes dataset. Section 6 concludes. The Supplementary Material is
organized as follows. Section A discusses the specific functional form of our
estimator. In section B, the parameters optimization scheme is explained.
Section C provides further insights into the results obtained by the
statistical tests used to validate our proposal. Finally, in section D we show
the CAESar performance on the stock dataset.

##  2 Problem Formulation and Related Literature

The target problem involves the dynamic prediction of the ES for univariate
time series. Consider a time series $\bm{y}=\\{y_{t}\\}_{t=1}^{T}$ and a
probability level $\theta\in[0,1]$. If $\bm{y}$ represents the asset returns
and $F_{t}$ the distribution function of $y_{t}$ conditioned on the
information set up to time $t-1$ ($\mathcal{F}_{t-1}$), then the VaR is
defined as the conditional quantile of $\bm{y}$, that is,
$VaR_{t}(\theta):=F_{t}^{-1}(\theta)$. Instead, the ES (historically known as
Conditional Value-at-Risk) is defined as the tail mean, that is,
$ES_{t}(\theta):=\mathbb{E}_{t-1}[y_{t}|y_{t}\leq VaR_{t}(\theta)]$, where
$\mathbb{E}_{t-1}[\cdot]$ represents the mean conditioned on
$\mathcal{F}_{t-1}$. In the following, when no confusion can arise, we omit
the explicit dependence of VaR and ES from $\theta$.
In the study of risk measures, a key role is played by the concept of
elicitability [7, 27]. In a nutshell, a functional is said to be elicitable if
it can be obtained as the minimum of a loss function. To be more precise, we
consider a set $Y\subset\mathbb{R}^{n}$ with the Borel $\sigma$-algebra. Let
$\mathcal{F}$ be a family of $Y$-valued random variables, and a functional
$R:\mathcal{F}\rightarrow\mathbb{R}^{m}$. The intuition is that $Y$ represents
the observable quantity, $F\in\mathcal{F}$ is a distribution for the elements
in $Y$, and $R$ is the risk measure we want to estimate. Next, consider a
function $\mathcal{L}:Y\times\mathbb{R}^{m}\rightarrow\mathbb{R}$.
$\mathcal{L}$ is said to be a scoring function if $\forall
z\in\mathbb{R}^{m},\ \forall F\in\mathcal{F}$, we have that
$\mathbb{E}\left[|\mathcal{L}(F,z)|\right]<\infty$. From a Machine Learning
perspective, $\mathcal{L}$ is a loss function. Specifically, we are interested
in a loss function that measures how well our guess $z$ fits the target
quantity $R$. So, we call a scoring function $\mathcal{L}$ strictly consistent
for the functional $R$ if $\forall z\in\mathbb{R}^{m},\ \forall
F\in\mathcal{F}$ with $z\neq R(F)$, we have that
$\mathbb{E}\left[\mathcal{L}(F,z)\right]>\mathbb{E}\left[\mathcal{L}(F,R(F))\right]$.
Finally, the functional $R$ is said to be elicitable if it has a strictly
consistent scoring function, that is if there exists a loss function that is
minimized by $R$.
The estimation of VaR in a time series context has been widely discussed in
the literature, starting from the seminal work on regression quantiles [31].
In fact, VaR is elicitable thanks to the well-known Pinball (Quantile) loss
function [31], defined as:

$\mathcal{L}^{\theta}_{q}(\bm{\hat{q}},\bm{y}):=\frac{1}{T}\sum_{t=1}^{T}(y_{t}-\hat{q}_{t})\left(\theta-\mathbbm{1}_{\\{y_{t}<\hat{q}_{t}\\}}\right)$
(1)

where $\bm{\hat{q}}=\\{\hat{q}_{t}\\}_{t=1}^{T}$ is the estimator of the
Value-at-Risk time series $\bm{VaR}=\\{VaR_{t}\\}_{t=1}^{T}$. Among the
different models for $\bm{\hat{q}}$, it is noteworthy to mention the
Conditional Autoregressive Value at Risk (CAViaR) [18]. It is an
autoregressive model, where the target variable $\hat{q}_{t}$ is regressed
against a (possibly multivariate) non-linear transformation of the lagged
observable prices $\bm{f}(y_{t-i})\in\mathbb{R}^{d},\ i=1,\cdots,p$ and the
estimates at previous time steps $\hat{q}_{t-j},\ j=1,\cdots,u$ that is:

$\hat{q}_{t}=\beta_{0}+\sum_{i=1}^{p}\bm{\beta}_{(i-1)d+1:id}\bm{\cdot}\bm{f}(y_{t-i})+\sum_{j=1}^{u}\beta_{j+pd}\hat{q}_{t-j}$
(2)

where
$\bm{\beta}_{(i-1)d+1:id}=\left(\beta_{(i-1)d+1},\cdots,\beta_{id}\right)$,
and $\bm{\cdot}$ represents the usual scalar product. Despite the great
flexibility of CAViaR, researchers usually focus on the case $p=u=1$, as it is
empirically shown to provide fast and reliable (that is, not prone to
overfitting) dynamics estimation. Furthermore, $\bm{f}$ is usually chosen as
the Asymmetric Slope (AS), which accounts for the asymmetry in asset returns
by mapping $y_{t-1}$ into the pair (hence $d=2$)
$\left((y_{t-1})^{+},(y_{t-1})^{-}\right)$, where
$(y_{t-1})^{+}:=\max(0,y_{t-1})$ and $(y_{t-1})^{-}:=\max(0,-y_{t-1})$,
respectively. Thus, the CAViaR specification (also used below as a baseline
model) is:

$\hat{q}_{t}=\beta_{0}+\beta_{1}(y_{t-1})^{+}+\beta_{2}(y_{t-1})^{-}+\beta_{3}\hat{q}_{t-1}$
(3)

As for the ES, the problem is far more complex. Indeed, estimating ES in a
dynamic context is an open problem in the literature, as ES is not elicitable
[24]. However, it has been proved the joint elicitability of VaR and ES pair
[20] (that is, the pair (VaR, ES) is an elicitable functional with values in
$\mathbb{R}^{2}$), giving rise to plenty of works in this field. Among others,
it is worth citing the approach proposed by Barrera et al. in [5], which is
named as BCGNS in the following. Given a probability space, the conditioning
event $X$ and the asset return $Y$ are seen as random variables with values in
a proper Polish (separable and completely metrizable) space $\mathcal{S}$ and
in $\mathbb{R}$. Then, the conditional VaR and ES are seen as Borel measurable
functions $\mathcal{S}\rightarrow\mathbb{R}$, $VaR(S)$ and $ES(S)$
respectively. The authors in [5] prove that the difference
$r(S):=ES(S)-VaR(S)$ is in the argmin, among all the Borel measurable non-
positive functions, of the following loss:

$\mathbb{E}\left[\left(r(S)+\frac{1}{\theta}\left(VaR(S)-Y\right)^{+}\right)^{2}\right]$
(4)

Thus, as the estimator $\bm{\hat{r}}$ of the time series difference
$\bm{r}:=\bm{ES}-\bm{VaR}$, it is used the argmin, within a fixed hypothesis
class $\mathcal{R}$, of the following loss function:

$\mathcal{L}^{\theta}_{B}(\bm{\hat{r}},\bm{y};\bm{VaR}):=\frac{1}{T}\sum_{t=1}^{T}\left(\hat{r}_{t}+\frac{1}{\theta}\left(VaR_{t}-y_{t}\right)^{+}\right)^{2}$
(5)

The loss can be interpreted as the mean square error between the
$\bm{\hat{r}}$ prediction and the ”excess loss” weighted by the confidence
level $\theta$. The excess loss is the difference between the asset value and
the quantile when there is a tail event, or 0 otherwise. As for the hypothesis
class, the authors in [5] use feedforward neural networks fed with previous
$y_{t}$ values. It is important to observe that this loss requires separate
training for VaR and ES estimators. The quantile needs to be estimated in the
first step. This task is achieved using a neural network and minimizing the
Pinball loss. Specifically, they use a network with 3 hidden layers, each one
with the double of the input layer neurons. Then, $\bm{\hat{q}}$ is fixed to
obtain the $\bm{ES}=\\{ES_{t}\\}_{t=1}^{T}$ estimates as
$\bm{\hat{e}}=\\{\hat{e}_{t}\\}_{t=1}^{T}$. It is noteworthy that the
computation cost in terms of training time is reduced by initializing the ES
network weights as equal to the ones of the VaR network, (provided that the
two networks share the same architecture). Despite the possibility of a joint
approach for the estimation of VaR and ES (as discussed below), such a two-
step approach represents a limitation, as claimed in [5], but the choice is
driven by technical reasons. Finally, it is worth underlining that the number
of lags used as input for the neural networks is a hyperparameter set by the
user that can strongly impact the overall result. Below, the model by [5] is
used as a comparative benchmark.
Another approach to ES estimation with a different choice for the loss
function is proposed by Patton et al. [39]. In this case, both ES and VaR are
jointly estimated by minimizing the following loss function:

$\mathcal{L}^{\theta}_{P}(\bm{\hat{e}},\bm{\hat{q}},\bm{y}):=\frac{1}{T}\sum_{t=1}^{T}\left(\frac{\hat{q}_{t}}{\hat{e}_{t}}-\frac{\hat{q}_{t}-y_{t}}{\theta\hat{e}_{t}}\mathbbm{1}_{\\{y_{t}\leq\hat{q}_{t}\\}}+\log\left(-\hat{e}_{t}\right)\right)$
(6)

The $\mathcal{L}^{\theta}_{P}$ loss can be viewed as a transformation of the
$\mathcal{L}^{\theta}_{B}$ one. Indeed, $\mathcal{L}^{\theta}_{P}$ can be
obtained by replacing $\bm{\hat{e}}-\bm{\hat{q}}$ and $\bm{\hat{q}}$ instead
of $\bm{\hat{r}}$ and $\bm{VaR}$ in the square root of (5) and normalizing by
the $\bm{\hat{e}}$ value. Then, a barrier term
($\log\left(-\hat{e}_{t}\right)$) is added to prevent the ES estimator from
assuming arbitrary big (in absolute) values. The hypothesis class used for
$\bm{\hat{q}}$ and $\bm{\hat{e}}$ is the Generalized Autoregressive Score
(GAS) model [10]. Unlike the GAS standard approach, here the forcing variable
is derived from the loss in Equation 6 rather than from the log-likelihood.
That is, the innovation term is obtained from the product between the inverse
of the Hessian and the gradient of the loss function with respect to the model
factors, net of algebraic simplifications. Specifically, the models proposed
in [39] include the one-factor and the two-factor models, which we indicate
with GAS1 and GAS2, respectively. The former uses a latent factor $k_{t}$ and
computes the target variables $\hat{q}_{t}$ and $\hat{e}_{t}$ as $k_{t}$
functions:

$\hat{q}_{t}=ae^{k_{t}};\quad\hat{e}_{t}=be^{k_{t}}\quad{\text{w}ith}\quad\
k_{t}=\beta
k_{t-1}+\frac{\gamma}{\hat{e}_{t-1}}\left(\frac{1}{\theta}\mathbbm{1}_{\\{y_{t-1}\leq\hat{q}_{t-1}\\}}y_{t-1}-\hat{e}_{t-1}\right)$
(7)

where the innovation term in the last equation is the so-called score, namely
a term proportional to the derivative of the loss function in Eq. (6) as a
function of $k_{t}$, see [10] for more details. In particular, notice that the
ratio $\frac{\hat{q}_{t}}{\hat{e}_{t}}=\frac{a}{b}$ is constant for the one-
factor model.
The two-factor model describes the factors as the target variables VaR and ES
directly, that is:

$\left[\begin{array}[]{c}\hat{q}_{t}\\\
\hat{e}_{t}\end{array}\right]=\left[\begin{array}[]{c}w_{1}\\\
w_{2}\end{array}\right]+\left[\begin{array}[]{c}b_{1}\hat{q}_{t-1}\\\
b_{2}\hat{e}_{t-1}\end{array}\right]+\left[\begin{array}[]{cc}a_{1,1}&a_{1,2}\\\
a_{2,1}&a_{2,2}\end{array}\right]\left[\begin{array}[]{c}\hat{q}_{t-1}\left(\theta-\mathbbm{1}_{\\{y_{t-1}\leq\hat{q}_{t-1}\\}}\right)\\\
\frac{y_{t-1}}{\theta}\mathbbm{1}_{\\{y_{t-1}\leq\hat{q}_{t-1}\\}}-\hat{e}_{t-1}\end{array}\right]$
(8)

where the innovation term is the score computed in the bivariate framework.
The parameters in Eqs. (7) and (8) are estimated by minimizing the loss
function in Eq. (6), under proper constraints to ensure the monotonicity
condition $\hat{e}_{t}<\hat{q}_{t}$. Below, we use the model [39] as a second
comparative benchmark.
One last approach for estimating the ES relates to the definition itself of
the ES as the tail mean: $ES_{t}(\theta):=\mathbb{E}_{t-1}[y_{t}|y_{t}\leq
VaR_{t}(\theta)]$. One might consider approximating the mean with the sample
average, which involves averaging the quantile estimates for the $n$ quantiles
$VaR(\theta_{j})$, $j=1,\cdots,n$. The confidence levels $\theta_{j}$ form an
equi-spaced partition of the tail $(0,\theta]$, that is:
$0<\theta_{1}<\cdots<\theta_{n}=\theta$. This method is suggested in [17] and
[33]. In the following, we have implemented it by estimating the quantiles
using both the CAViaR (we name this approach K-CAViaR) or the Quantile
Regression Neural Network (giving rise to the K-QRNN approach) -that is, a
feedforward neural network trained to minimize the Pinball loss in (1).
Observe that no monotonicity constraints are explicitly imposed, so,
theoretically, this approach is subject to the above-mentioned crossing
quantiles problem. However, the key intuition is that quantile estimation
errors for different $\theta_{j}$ are averaged out. Lastly, it is worth
mentioning a couple of recent approaches to ES estimation that, even if not
directly linked to the context of this paper, answer interesting questions. In
particular, a branch of literature focuses on ES estimation from the
perspective of extreme value theory. For example, the authors of [11] propose
an estimator for the tail index, namely a statistical index characterizing
heavy-tailed distributions. The same index is then used for the estimation of
the expectiles of a distribution, connected to the ES inference. The estimator
is then exploited to study the SPX index under different market regimes,
showing that the fitted parameters strongly change when looking at the periods
before, during, and after the 2008 crisis. A similar perspective is followed
in [37]. The peculiarity of this work is considering covariates in the
estimation of the tail index. Other works focus on the joint estimation of VaR
and ES via the maximization of the Asymmetric Laplace log-likelihood.
Noteworthy examples are [42] and [35]. The former introduces the approach,
while the latter generalizes it to the multivariate case. Both of them assess
the goodness of their model by using a real-world dataset made up of stock
indexes.
The authors in [21] compare univariate and multivariate ES predictions.
Specifically, they attempt to predict the weekly ES of a stock portfolio. The
univariate approach directly works with the portfolio time series. Instead,
the multivariate approach proposes a regression analysis for the portfolio
against the Fama–French and momentum factors. Then, a copula model is used to
learn the portfolio ES based on these factors. A noteworthy conclusion
suggests that the univariate approach performs the best. This is mainly due to
the high variance of estimation errors in the multivariate approach that lead
to wrong forecasts.
The study in [22] attempts to model the arrival intensity and magnitude of
tail events by using a score-driven approach. The proposed model is then
tested on a real-world dataset, showing promising results as the event
dynamics seem to capture the major financial crisis. That is, there are spikes
in the arrivals and magnitude forecasts corresponding to well-recognized
turbulent financial periods, such as the 2008 crisis.

##  3 The CAESar Model

In this Section, we present the CAESar model. First, we describe the general
setting, then we highlight the differences with existing competitors, and
finally, we discuss the general idea behind the optimization scheme. More
details can be found in the Supplementary Material, sections A and B, which
provide a deeper insight into the modeling choices.
The main idea behind CAESar is to generalize the CAViaR regression to convert
it into a joint (VaR, ES) estimator. Specifically, both the targets
$\hat{q}_{t}$ and $\hat{e}_{t}$ are regressed against multi-dimensional non-
linear transformations of the asset value at previous time steps $t-i,\
i=1,\cdots,p$. We indicate these transformations with $\bm{f}$ for the
quantile and $\bm{g}$ for the ES. Furthermore, $\hat{q}_{t}$ and $\hat{e}_{t}$
are also regressed on previous estimates $\hat{q}_{t-j},\ j=1,\cdots,u$ and
$\hat{e}_{t-l},\ l=1,\cdots,v$ to account for heteroskedastic effects. The
general form of the proposed estimator is

$\displaystyle\hat{q}_{t}=\beta_{0}+\sum_{i=1}^{p}\bm{\beta}_{(i-1)d+1:id}\bm{\cdot}\bm{f}(y_{t-i})+\sum_{j=1}^{u}\beta_{j+pd}\hat{q}_{t-j}+\sum_{l=1}^{v}\beta_{l+pd+u}\hat{e}_{t-l}$
(9)
$\displaystyle\hat{e}_{t}=\gamma_{0}+\sum_{i=1}^{p}\bm{\gamma}_{(i-1)d+1:id}\bm{\cdot}\bm{g}(y_{t-i})+\sum_{j=1}^{u}\gamma_{j+pd}\hat{q}_{t-j}+\sum_{l=1}^{v}\gamma_{l+pd+u}\hat{e}_{t-l}$
(10)

where two distinct sets of parameters
$\bm{\beta}=(\beta_{0},\cdots,\beta_{pd+u+v})$ and
$\bm{\gamma}=(\gamma_{0},\cdots,\gamma_{pd+u+v})$ describe the dynamics of
$\bm{\hat{q}}$ and $\bm{\hat{e}}$, respectively. As before,
$\bm{\beta}_{(i-1)d+1:id}$ represents the d-dimensional vector
$(\beta_{(i-1)d+1},\cdots,\beta_{id})$ and $\bm{\gamma}_{(i-1)d+1:id}$ is for
$(\gamma_{(i-1)d+1},\cdots,\gamma_{id})$. For simplicity, it is convenient to
choose $\bm{f}\equiv\bm{g}$ and to set $u=v$. Thus, the joint VaR-ES estimator
can be read as:

$\left[\begin{array}[]{c}\hat{q}_{t}\\\
\hat{e}_{t}\end{array}\right]=\left[\begin{array}[]{c}\beta_{0}\\\
\gamma_{0}\end{array}\right]+\sum_{i=1}^{p}\left[\begin{array}[]{c}\bm{\beta}_{(i-1)d+1:id}\\\
\bm{\gamma}_{(i-1)d+1:id}\end{array}\right]\bm{f}(y_{t-i})+\sum_{j=1}^{u}\left[\begin{array}[]{cc}\beta_{pd+j}&\beta_{pd+u+j}\\\
\gamma_{pd+j}&\gamma_{pd+u+j}\end{array}\right]\left[\begin{array}[]{c}\hat{q}_{t-j}\\\
\hat{e}_{t-j}\end{array}\right]$ (11)

It is worth noting that the specification we propose for the ES is similar to
the GAS models in [39] (Eqs. (7) and (8)) as they are both autoregressive
models. However, our formulation is more general and less limiting. As such,
it should be able to capture better the tail behavior. In fact, the one-factor
model in (7) assumes that the ratio $\frac{Var_{t}}{ES_{t}}$ is constant.
However, in principle, this ratio can be nonconstant in time, depending on the
evolution of the shape of the returns’ conditional distribution, while the
previous assumption might lead to model misspecification errors. For instance,
the ratio for a non-centered normal distribution changes as a function of the
standard volatility. Moreover, the ratio for a centered generalized Student’s
t distribution varies with the degrees of freedom. These considerations are
illustrated in Figure 1.

Figure 1: Evolution of the ratio $\frac{Var_{t}}{ES_{t}}$ for $\theta=0.05$ as
a function of the distribution parameters. In the left plot, the ratio for a
standard normal distribution is plotted as a function of the volatility
$\sigma$. The right plot shows the ratio for a centred generalized Student’s t
distribution against the degrees of freedom $\nu$.

The GAS two-factor model in (8) addresses this issue but restricts the
information set to the magnitude of $y_{t}$ in the case of a tail event, i.e.
when the realization is below the quantile. Furthermore, as stated also in the
original paper [39], when tested on empirical data, the two-factor model is
outperformed by the one-factor model. Below, we show using both simulations
and empirical data that our specification improves significantly the GAS
models by providing a more robust and flexible framework.
On the other side, the model proposed by Barrera et al. in [5] exploits neural
networks to catch non-linear dependencies in the data. However, it uses dense
layers, so it exploits information only from the more recent time steps -that
is, the lags used as input of the network. In contrast, CAESar exploits an
exponentially weighted moving average that keeps memory of the entire time
series’ past history. Moreover, the model estimation (see below) is far
simpler and less prone to hyperparameter misspecification. A possible drawback
of our approach refers to the high sensibility to the choice of the initial
parameters -similar to the case of GAS models in [39]. To mitigate this issue,
we propose a three-step approach for the minimization of the loss in (6).
First the $\bm{VaR}$ prediction, $\bm{\hat{q}}$, is obtained using the CAViaR
regression. Then, in the second step, quantile information is used to
construct the $\bm{ES}$ estimator $\bm{\hat{e}}$. This is done following the
idea in [5], that is by modelling the residual time series $\bm{\hat{r}}$:

$\hat{r}_{t}=\hat{e}_{t}-\hat{q}_{t}\quad{\text{w}ith}\quad\
\hat{r}_{t}=\tilde{\gamma_{0}}+\sum_{i=1}^{p}\bm{\tilde{\gamma}}_{(i-1)d+1:id}\bm{\cdot}\bm{f}(y_{t-i})+\sum_{j=1}^{u}\tilde{\gamma}_{j+pd}\hat{q}_{t-j}+\sum_{l=1}^{u}\tilde{\gamma}_{l+pd+u}\hat{r}_{t-l}$
(12)

As for the loss function, it is derived from the joint elicitability of $VaR$
and $ES$. Specifically, Equation (5) is used with an added soft constraint
term to enforce the monotonicity condition $\hat{e}_{t}<\hat{q}_{t}$:

$\mathcal{L}^{\theta}_{r}(\bm{\hat{r}},\bm{y};\bm{VaR}):=\mathcal{L}^{\theta}_{B}(\bm{\hat{r}},\bm{y};\bm{VaR})+\lambda_{r}\sum_{t=1}^{T}(\hat{r}_{t})^{+}$
(13)

In the last step, the previously estimated coefficients are used as a suitable
starting point for minimizing the loss in (6). Indeed, by rearranging
Equations (2) and (12), the formulation for $\hat{e}_{t}$ can be recovered:

$\hat{e}_{t}=(\tilde{\gamma_{0}}+\beta_{0})+\sum_{i=1}^{p}\left(\left(\bm{\tilde{\gamma}}+\bm{\beta}\right)_{(i-1)d+1:id}\bm{\cdot}\bm{f}(y_{t-i})\right)+\sum_{j=1}^{u}\left(\tilde{\gamma}_{pd+j}+\beta_{pd+j}-\tilde{\gamma}_{pd+u+j}\right)\hat{q}_{t-j}+\sum_{j=1}^{u}\gamma_{pd+u+j}\hat{e}_{t-j}$

By comparing the previous expression with (11), we can recover the
$\gamma_{j}$ coefficients as:

$\gamma_{j}=\begin{cases}\tilde{\gamma}_{j}+\beta_{j}&j=0,\cdots,pd\\\
\gamma_{j}=\tilde{\gamma}_{j}+\beta_{j}-\tilde{\gamma}_{j+u}&j=pd+1,\cdots,pd+u\\\
\gamma_{j}=\tilde{\gamma}_{j}&j=pd+u+1,\cdots,pd+2u\end{cases}$

The optimal coefficients for $\bm{\hat{q}}$ and $\bm{\hat{e}}$ obtained in the
previous steps serve as the starting point for the final step. We jointly
optimize the two estimators coefficients in Equation (11). The additional ones
$\beta_{pd+u+1},\cdots,\beta_{pd+2u}$ (i.e. those absent in the original
CAViaR formulation in Eq. (2)) are initialized to zero. The loss function used
is the one in Eq. (6), with the addition of soft constraints ensuring (i) the
monotonicity condition and (ii) the non-positivity of the quantile, which is a
straightforward assumption for financial returns:

$\mathcal{L}^{\theta}_{q,e}(\bm{\hat{e}},\bm{\hat{q}},\bm{y}):=\mathcal{L}^{\theta}_{P}(\bm{\hat{e}},\bm{\hat{q}},\bm{y})+\lambda_{e}\sum_{t=1}^{T}(\hat{e}_{t}-\hat{q}_{t})^{+}+\lambda_{q}\sum_{t=1}^{T}(\hat{q}_{t})^{+}$
(14)

At this point, one may wonder whether the updating in the quantile predictor
leads to a degradation in the VaR estimate. We show that this is not the case,
as discussed in Section 5.2, but there is a slight improvement in the
forecasting performance.
Thus far, the CAESar regression has been discussed in the most general
setting. However, for practical implementations, we suggest focusing on a
particular specification: the AS-(1,1) model -i.e. AS specification for
$\bm{f}$ and $p=u=1$. Indeed, the same considerations previously made for the
CAViaR hold here, as the AS specification is able to capture asymmetries in
the return distribution. Finally, a Markovian specification of the model,
namely using the information on the previous time lag only, reduces both the
training time and the overfitting risk; at the same time, empirical
applications corroborate this modeling choice, displaying very good
forecasting performances. The estimator adopted in the applications below is:

$\left[\begin{array}[]{c}\hat{q}_{t}\\\
\hat{e}_{t}\end{array}\right]=\left[\begin{array}[]{c}\beta_{0}\\\
\gamma_{0}\end{array}\right]+\left[\begin{array}[]{cc}\beta_{1}&\beta_{2}\\\
\gamma_{1}&\gamma_{2}\end{array}\right]\left[\begin{array}[]{c}(y_{t-1})^{+}\\\
(y_{t-1})^{-}\end{array}\right]+\left[\begin{array}[]{cc}\beta_{3}&\beta_{4}\\\
\gamma_{3}&\gamma_{4}\end{array}\right]\left[\begin{array}[]{c}\hat{q}_{t-1}\\\
\hat{e}_{t-1}\end{array}\right]$ (15)

A discussion on the adopted specifications for both the CAESar and CAViaR
models can be found in the Supplementary Material A. And empirical
justification for using the cross terms in Equation (15) is also in
Supplementary Material A. The detailed report of the optimization scheme is
presented in Supplementary Material B. The Python implementation of the
proposed methodology is available online; see the “Data and Code Availability”
Section below.

##  4 Statistical Tests

As discussed in the Introduction, assessing the effectiveness of an ES
estimation model is a challenging task due to the lack of a universally
accepted test. Furthermore, comparing different models to determine whether
one statistically outperforms the other is not straightforward. An attempt to
clarify this is made in [13], where existing tests are surveyed and
classified. Following their taxonomy, we evaluate our proposal using two
separate classes of tests. Some tests directly evaluate the ES approximation,
aiming to assess whether the model’s estimation is coherent with the
definition of the ES itself. Within this context, there is no need for a loss
function as an evaluation metric. However, it is still not clear how to study
the limiting behavior of the considered statistics. Hence, a bootstrap
approach is used to approximate the test p-value (see [43] for further
details). On the other side, some testing procedures have been introduced for
model comparison. Such a comparison is typically carried out using the
previously discussed loss functions. In other words, these tests are based on
the concept of joint elicitability of VaR and ES.

###  4.1 Statistical Tests - Direct Approximation

The first category of tests is those that directly evaluate the ES
approximation. Specifically, three statistical tests are conducted. We
consider first the McNeil and Frey test [34], which is a one-sided test for
the null hypothesis of no risk underestimation. The intuition is testing the
definition of ES as a tail mean. In particular, from the index set
$\mathcal{T}=\\{1,\cdots,T\\}$, we extract the subset of quantile violations
$\tilde{\mathcal{T}}=\\{t\in\mathcal{T}\ s.t.\ y_{t}<\hat{q}_{t}\\}$. Then, on
this subset, we use a bootstrapping approach to test the null hypothesis that
$y_{t}-\hat{e}_{t}$ has a non-negative mean, i.e.
$H_{0}:\mathbb{E}\left[y_{t}-\hat{e}_{t}|t\in\tilde{\mathcal{T}}\right]\geq
0$. Specifically, we generate 10,000 samples (as detailed in [13]) to estimate
the “empirical” p-value for the bootstrapped distribution of the statistic.
Acerbi and Szekely [1] proposed two similar tests using the $Z_{1}$ and
$Z_{2}$ statistics . In both cases, we consider a two-sided test for the null
hypothesis of correct ES predictions. Specifically, both statistics are
derived by rearranging the ES definition. For $Z_{1}$:

$ES_{t}=\mathbb{E}\left[y_{t}|y_{t}\leq
VaR_{t}\right]\implies\mathbb{E}\left[\frac{y_{t}}{ES_{t}}\bigg{|}y_{t}\leq
VaR_{t}\right]=1$ (16)

For $Z_{2}$, we use the conditional probability formula and the definition of
VaR to obtain:

$ES_{t}=\mathbb{E}\left[y_{t}|y_{t}\leq
VaR_{t}\right]=\mathbb{E}\left[\frac{y_{t}}{\theta}\mathbbm{1}_{\\{y_{t}\leq
VaR_{t}\\}}\right]\implies\mathbb{E}\left[\frac{y_{t}}{ES_{t}\theta}\mathbbm{1}_{\\{y_{t}\leq
VaR_{t}\\}}\right]=1$ (17)

Based on the previous equations, the following statistics are defined:

$Z_{1}:=\frac{1}{|\tilde{\mathcal{T}}|}\sum_{t\in\tilde{\mathcal{T}}}\frac{y_{t}}{\hat{e}_{t}};\quad\quad
Z_{2}:=\frac{1}{T}\sum_{t\in\tilde{\mathcal{T}}}\frac{y_{t}}{\hat{e}_{t}\theta}.$
(18)

In both cases, the null hypothesis assumes the statistic to be equal to 1. We
employ a standard bootstrap approach to translate this into a testing
procedure. Again, 10000 bootstrap samples are generated to evaluate the
“empirical” p-value for the bootstrapped distribution of the statistic.
Finally, it is worth noting that, according to the empirical analysis carried
out in [13], in our setting (six years of training/validation and one year of
testing for each fold), the $Z_{1}$ and McNeil and Frey tests are expected to
exhibit the best size and power.

###  4.2 Statistical Tests - Model Comparison

We discuss here a number of testing procedures for comparing ES estimators.
The underlying null hypothesis is about equal predictive accuracy in comparing
two ES predictors. A rejection of the null signals one model outperforming the
other (alternative hypothesis). It is important to underline that this kind of
test tends to be ’conservative’, resulting in a low number of rejections [13],
a phenomenon also observed in the applications below, see Section 5.
Here, we consider four different testing procedures. First, we apply the
Diebold-Mariano test [14], which is a well-known two-sided statistical test
used to compare two time series forecasting estimators, A and B. The null
hypothesis is that A and B are equivalent. The focus lies on the time series
defined by the difference between the losses associated with the two
predictors, under the assumption that such a difference is a stationary
process. The statistic is then defined as the mean value of the difference,
normalized by the sum of the autocorrelations up to a certain lag.
Furthermore, we adjust the statistic with the small-size sample correction
discussed in [26], a method also known as the Harvey adjustment.
The second test we consider is the corrected resampled Student’s t-test,
introduced in [36]. This statistical test is specifically designed for
frameworks where cross-validation is used, as in our case (further details on
the experimental setup are in the next section). It is the only test where we
do not perform an iteration for each fold. Rather, we conduct the test only
once. The population used is made up of the losses among the folds. As for the
test statistic, it is obtained from the Student’s t-test by properly rescaling
the variance. This scaling, which is smaller than 1, embodies the fact that,
when performing cross-validation, the variance in the predictions (and thus,
in the losses) is also affected by the variance in the train set. As the test
is designed to take into account only the volatility in the out-of-sample
data, we have to filter the total variance to remove its train component. This
task is achieved by the scaling factor. The third considered test is the Loss
Difference test. It aims to determine if the difference between the losses
obtained by two predictors (A and B) is statistically smaller than 0. In other
words, the null hypothesis is that the loss difference
$\mathcal{L}\left(\bm{\hat{e}}^{A},\bm{\hat{q}}^{A},\bm{y}\right)-\mathcal{L}\left(\bm{\hat{e}}^{B},\bm{\hat{q}}^{B},\bm{y}\right)$
is greater or equal to 0. If the null is rejected, we can say that the
predictor $A$ outperforms $B$. This test is performed by considering the
population of the loss at time $t\ \forall t\in\mathcal{T}$, and studying its
mean with a one-sided bootstrap test, following the same rationale as in the
McNeil and Frey test.
The last approach is based on the Encompassing test, which generally indicates
a broad class of testing procedures aimed to determine if predictor A
outperforms predictor B by studying a linear combination of the two predictors
[30]. Specifically, in our case, the idea is to introduce a new predictor, C,
defined as a linear combination of A and B: $C=\alpha A+\beta B$. The non-
negative weights $\alpha$ and $\beta$ are determined through an Ordinary Least
Squares regression applied to the first half of the data in the test set. The
forecast provided by C is then evaluated in the second half of the data. As
long as the information in B is already incorporated in A, adding B should
only add noise. Thus, the ultimate prediction provided by C should be inferior
to that of A. Thus, we perform a bootstrapping exercise similar to the Loss
Difference test, but now we compare predictors A and C. If the null is
rejected, then the addition of predictor B only adds noise, indicating that A
outperforms B.

##  5 Results

In this section, we assess the validity of the proposed framework by comparing
the CAESar model with different benchmarks both on synthetic and real-world
datasets. We consider five competitor models of the proposed CAESar model:

  * •

K-CAViaR Based on the definition of the $ES(\theta)$ as the mean value of the
quantiles below the target confidence $\theta$, we construct an $ES(\theta)$
estimator by averaging the CAViaR predictions for $n=10$ equispaced quantiles
$VaR(\theta_{j})\ j=1,\cdots,n$ with $0<\theta_{1}<\cdots<\theta_{n}=\theta$.
This is the same approach proposed in [33]. The implementation is a simple
extension of the CAViaR methodology; the CAViaR code is available online at
http://www.simonemanganelli.org/Simone/Research.html.

  * •

BCGNS We use the BCGNS as presented in [5] by using the code available online
at https://github.com/BouazzaSE/Learning-VaR-and-ES. The model is fed by the
previous 25 lags of each timestep, as done in the original paper.

  * •

K-QRNN This method mirrors the K-CAViaR approach. One single network is
trained, with $n$ outputs, one for each quantile $VaR(\theta_{j})\
j=1,\cdots,n$. To preserve the coherence with the other neural network model,
i.e. the BCGNS, we use as input 25 lagged values.

  * •

GAS1 and GAS2 We employ both the GAS one-factor and two-factor models proposed
in [39]. The methodology has been implemented by ourselves and is available
online; see Data and Code Availability Section below.

The details of the comparison are in the following subsections.

###  5.1 Simulation Study

First, we carry out a simulation experiment to assess whether the predictions
provided by the ES models are coherent with the ground truth. Specifically, we
use the GARCH as the data-generating process, with both Gaussian and Student’s
t innovations. Specifically:

$GARCH(1,1):\quad y_{t}=\sigma_{t}\varepsilon_{t}\quad
with\quad\sigma_{t}=\omega+\beta\sigma_{t-1}^{2}+\gamma\varepsilon_{t-1}^{2}$
(19)

The coefficients of the process with Gaussian innovation are
$(\omega,\beta,\gamma)$. When considering the t innovations, also the degrees
of freedom $\nu$ have to be considered. These parameters are fitted from the
real data (see next section for details). Specifically, we fit the parameters
on the SPX (I set of coefficients), FTSE (II set of coefficients), and DAX
(III set of coefficients). Thus, there are two different innovations and, for
each of them, three sets of coefficients, i.e. we consider six data-generating
processes. For each of them, we generate 20 series made up of 1750 points,
which roughly correspond to seven years, divided into 1500 observations in the
train set and 250 in the test set.
The algorithms are tested on three different confidence levels $\theta$,
namely $\theta=0.05,0.025,0.01$. Specifically, the ground truth for the
Gaussian innovations reads:

$e_{t}=-\frac{\sigma_{t}}{\theta}f_{N}\left(F^{-1}_{N}(1-\theta)\right)$ (20)

where $F_{N}$ and $f_{N}$ are, respectively, the cumulative and probability
density functions of a standard normal distribution. Instead, as for the
Student’s t innovations:

$e_{t}=-\sigma_{t}\frac{\nu+\left(F^{-1}_{T}(\theta)\right)^{2}}{(\nu-1)\theta}\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu}}\left(1+\frac{\left(F^{-1}_{T}(\theta)\right)^{2}}{\nu}\right)^{-\frac{\nu+1}{2}}$
(21)

where $F_{T}$ is the cumulative distribution function of a standard
t-distribution. The ES forecasts are then compared with the ground truth using
both the Mean Absolute Error (MAE), the Root Mean Squared Error (RMSE), the
loss functions $\mathcal{L}^{\theta}_{B}$ (Equation (5)), and the loss
$\mathcal{L}^{\theta}_{P}$ (Equation (6)).

$\bm{\theta=0.05}$ |  | N | t |  | N | t
---|---|---|---|---|---|---
| I | II | III | I | II | III |  | I | II | III | I | II | III
CAESar | MAE | 0.002 | 0.0022 | 0.0017 | 0.3423 | 0.0282 | 0.0081 | $\mathcal{L}_{B}^{\theta}$ | 0.0007 | 0.001 | 0.0005 | 25.9028 | 0.9505 | 0.0267
RMSE | 0.003 | 0.003 | 0.002 | 0.433 | 0.038 | 0.011 | $\mathcal{L}_{P}^{\theta}$ | 0.781 | 0.955 | 0.657 | 3.369 | 2.362 | 1.633
K-CAViaR | MAE | 0.0015 | 0.0017 | 0.0012 | 0.1317 | 0.0345 | 0.0096 | $\mathcal{L}_{B}^{\theta}$ | 0.0007 | 0.0009 | 0.0005 | 28.2971 | 1.2366 | 0.0285
RMSE | 0.002 | 0.002 | 0.002 | 0.19 | 0.044 | 0.012 | $\mathcal{L}_{P}^{\theta}$ | 0.779 | 0.948 | 0.658 | 3.247 | 2.372 | 1.634
BCGNS | MAE | 0.0297 | 0.0332 | 0.0356 | 0.829 | 0.1453 | 0.0542 | $\mathcal{L}_{B}^{\theta}$ | 0.0047 | 0.0066 | 0.0062 | 97.4735 | 3.3069 | 0.071
RMSE | 0.047 | 0.053 | 0.059 | 1.101 | 0.189 | 0.09 | $\mathcal{L}_{P}^{\theta}$ | 2.553 | 2.987 | 2.365 | 5.578 | 5.779 | 3.643
K-QRNN | MAE | 0.0083 | 0.0085 | 0.0065 | 0.9368 | 0.1767 | 0.0371 | $\mathcal{L}_{B}^{\theta}$ | 0.0008 | 0.0009 | 0.0006 | 198.3513 | 10.0393 | 0.065
RMSE | 0.01 | 0.01 | 0.008 | 1.163 | 0.2 | 0.044 | $\mathcal{L}_{P}^{\theta}$ | 0.94 | 1.076 | 0.795 | 9.001 | 5.187 | 2.693
GAS1 | MAE | 0.0892 | 0.0291 | 0.081 | 0.3503 | 0.0684 | 0.0222 | $\mathcal{L}_{B}^{\theta}$ | 0.0013 | 0.0376 | 0.4944 | 46.5417 | 0.7035 | 0.0375
RMSE | 0.011 | 0.068 | 0.01 | 0.167 | 0.1 | 0.028 | $\mathcal{L}_{P}^{\theta}$ | 1.36 | 1.54 | 1.514 | 4.028 | 2.46 | 1.846
GAS2 | MAE | 0.0199 | 0.2465 | 0.0044 | 0.1288 | 0.0724 | 0.0385 | $\mathcal{L}_{B}^{\theta}$ | 0.0037 | 0.1283 | 0.0007 | 92.1 | 0.9577 | 0.6581
RMSE | 0.021 | 0.248 | 0.005 | 0.15 | 0.09 | 0.047 | $\mathcal{L}_{P}^{\theta}$ | 1.543 | 2.529 | 1.463 | 4.182 | 3.795 | 2.225

$\bm{\theta=0.025}$ |  | N | t |  | N | t
---|---|---|---|---|---|---
| I | II | III | I | II | III |  | I | II | III | I | II | III
CAESar | MAE | 0.0022 | 0.0025 | 0.0018 | 0.2286 | 0.0385 | 0.011 | $\mathcal{L}_{B}^{\theta}$ | 0.0011 | 0.0016 | 0.0009 | 57.4856 | 3.1664 | 0.083
RMSE | 0.003 | 0.003 | 0.002 | 0.329 | 0.052 | 0.015 | $\mathcal{L}_{P}^{\theta}$ | 0.916 | 1.081 | 0.795 | 3.458 | 2.738 | 1.899
K-CAViaR | MAE | 0.0017 | 0.0018 | 0.0013 | 0.0911 | 0.0432 | 0.0125 | $\mathcal{L}_{B}^{\theta}$ | 0.0012 | 0.0014 | 0.0008 | 56.2634 | 3.0354 | 0.0856
RMSE | 0.002 | 0.002 | 0.002 | 0.129 | 0.059 | 0.016 | $\mathcal{L}_{P}^{\theta}$ | 0.919 | 1.097 | 0.798 | 3.356 | 2.748 | 1.904
BCGNS | MAE | 0.0258 | 0.0331 | 0.0365 | 0.8211 | 0.1742 | 0.0655 | $\mathcal{L}_{B}^{\theta}$ | 0.0049 | 0.0094 | 0.0094 | 71.0887 | 4.0117 | 0.1256
RMSE | 0.043 | 0.051 | 0.058 | 1.085 | 0.223 | 0.105 | $\mathcal{L}_{P}^{\theta}$ | 2.651 | 2.985 | 2.358 | 4.881 | 5.104 | 4.098
K-QRNN | MAE | 0.0062 | 0.007 | 0.0071 | 0.9112 | 0.1423 | 0.0462 | $\mathcal{L}_{B}^{\theta}$ | 0.0014 | 0.0016 | 0.0011 | 87.2144 | 7.1658 | 0.1064
RMSE | 0.007 | 0.008 | 0.008 | 1.157 | 0.166 | 0.054 | $\mathcal{L}_{P}^{\theta}$ | 0.916 | 1.062 | 0.801 | 7.215 | 4.668 | 2.986
GAS1 | MAE | 0.0534 | 0.0188 | 0.0792 | 0.3301 | 0.0885 | 0.0264 | $\mathcal{L}_{B}^{\theta}$ | 0.0058 | 0.0521 | 0.6869 | 25.0527 | 1.5622 | 0.0956
RMSE | 0.085 | 0.057 | 0.077 | 0.432 | 0.11 | 0.038 | $\mathcal{L}_{P}^{\theta}$ | 1.585 | 2.471 | 1.516 | 3.148 | 2.724 | 2.415
GAS2 | MAE | 0.0177 | 0.2184 | 0.0036 | 0.0795 | 0.0742 | 0.0453 | $\mathcal{L}_{B}^{\theta}$ | 0.0039 | 0.1532 | 0.0009 | 58.2252 | 2.4199 | 0.9466
RMSE | 0.02 | 0.219 | 0.004 | 0.092 | 0.089 | 0.053 | $\mathcal{L}_{P}^{\theta}$ | 1.704 | 2.628 | 1.602 | 3.215 | 3.846 | 2.874

$\bm{\theta=0.01}$ |  | N | t |  | N | t
---|---|---|---|---|---|---
| I | II | III | I | II | III |  | I | II | III | I | II | III
CAESar | MAE | 0.0028 | 0.0036 | 0.003 | 0.2857 | 0.0625 | 0.0172 | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.0031 | 0.0015 | 174.7012 | 15.5411 | 0.2973
RMSE | 0.004 | 0.005 | 0.004 | 0.424 | 0.082 | 0.023 | $\mathcal{L}_{P}^{\theta}$ | 1.056 | 1.223 | 0.937 | 3.729 | 2.829 | 2.171
K-CAViaR | MAE | 0.0032 | 0.0042 | 0.0028 | 0.3341 | 0.0672 | 0.0242 | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.003 | 0.0014 | 288.1228 | 22.8529 | 0.426
RMSE | 0.004 | 0.005 | 0.003 | 0.436 | 0.085 | 0.03 | $\mathcal{L}_{P}^{\theta}$ | 1.083 | 1.239 | 0.937 | 3.681 | 2.839 | 2.181
BCGNS | MAE | 0.0119 | 0.0133 | 0.0096 | 1.0183 | 0.1483 | 0.04 | $\mathcal{L}_{B}^{\theta}$ | 0.0147 | 0.0145 | 0.0102 | 672.7932 | 20.6229 | 0.5798
RMSE | 0.017 | 0.018 | 0.013 | 1.336 | 0.183 | 0.051 | $\mathcal{L}_{P}^{\theta}$ | 2.001 | 2.778 | 1.8 | 7.743 | 4.573 | 3.316
K-QRNN | MAE | 0.0129 | 0.0129 | 0.01 | 1.3303 | 0.257 | 0.0552 | $\mathcal{L}_{B}^{\theta}$ | 0.0013 | 0.0013 | 0.001 | 2320.0497 | 219.8131 | 0.6772
RMSE | 0.014 | 0.014 | 0.011 | 1.664 | 0.29 | 0.065 | $\mathcal{L}_{P}^{\theta}$ | 1.27 | 1.391 | 1.111 | 18.16 | 11.147 | 3.494
GAS1 | MAE | 0.0263 | 0.0584 | 0.0201 | 0.1313 | 0.1911 | 0.0581 | $\mathcal{L}_{B}^{\theta}$ | 0.0684 | 0.6865 | 0.1794 | 67.8674 | 0.5957 | 0.2118
RMSE | 0.031 | 0.023 | 0.021 | 0.171 | 0.284 | 0.081 | $\mathcal{L}_{P}^{\theta}$ | 1.435 | 6.035 | 1.767 | 4.474 | 3.964 | 2.416
GAS2 | MAE | 0.4575 | 0.3327 | 0.1983 | 0.4976 | 0.3689 | 0.2655 | $\mathcal{L}_{B}^{\theta}$ | 0.1853 | 0.3492 | 0.0177 | 0.9937 | 1.3381 | 0.3063
RMSE | 0.458 | 0.334 | 0.2 | 0.508 | 0.337 | 0.271 | $\mathcal{L}_{P}^{\theta}$ | 2.23 | 2.204 | 1.661 | 4.902 | 4.783 | 3.573

Table 1: Simulation results. The algorithms are compared over six different
generating processes. All the generating processes are based on GARCH, both
with Gaussian innovations (N) and t innovations (t). The coefficients of the
process are fitted on the SPX (I set of coefficients), FTSE (II set), or DAX
(III set). The value in each cell represents the mean loss of the row
predictor across the twenty time series generated according to the column
process. The best result is in bold, the second to best is underlined.

The results are displayed in Table 1. The MAE and RMSE values indicate that
the best accuracy results are achieved by CAESar and K-CAViaR. Specifically,
K-CAViaR outperforms for $\theta=0.05$, CAESar outperforms for $\theta=0.01$,
and there is a balanced situation for middle $\theta$. It is worth noting
that, as stated in [39], the GAS models exhibit strong instability in the
parameters optimization stage. Specifically, we have found that for some
series, the dynamic given by these models is strongly unrealistic. Thus, in
the comparison, we have not kept into account these series when computing the
GAS losses.
If we look at the $\mathcal{L}^{\theta}_{B}$ and $\mathcal{L}^{\theta}_{P}$
values, we recover a predominance of CAESar and K-CAViaR. We also find that
often, one of them outperforms the other according to a loss, but the
situation is reverted according to the other loss. Nonetheless, comparing the
best algorithms according to the two classes of metrics (MAE and RMSE on one
side, $\mathcal{L}^{\theta}_{B}$ and $\mathcal{L}^{\theta}_{P}$ on the other
side), we find that the $\mathcal{L}^{\theta}_{P}$ used by Patton et al. in
[39] provides results which are more similar to those of MAE and RMSE. Such an
incoherence between the $\mathcal{L}^{\theta}_{B}$ and
$\mathcal{L}^{\theta}_{P}$ losses can be explained by taking in mind that they
are designed for different purposes. Indeed, $\mathcal{L}^{\theta}_{P}$ is
constructed to jointly estimate the VaR and ES. Instead,
$\mathcal{L}^{\theta}_{B}$ is defined to estimate the ES given the knowledge
of the VaR. However, we do not have the true VaR, but just an estimate for it.
This estimate could sometimes not be correct or even show a large error. Thus,
in such a context, the value expressed by the $\mathcal{L}^{\theta}_{B}$ loss
may not be indicative of the true predictor’s capability of approximating the
ES. To verify this, we evaluate the $\mathcal{L}^{\theta}_{B}$ loss by using
the actual $VaR_{t}$ instead of its estimate $\hat{q}_{t}$. Then, we count how
many times the Barrera and Patton losses agree on which algorithm is the best
and the second to best. This happens 19 times over a sample of 20 time series.
On the contrary, when using the quantile estimate $\hat{q}_{t}$, the agreement
number is reduced to 9, thus confirming the above intuition.

###  5.2 Empirical application

####  5.2.1 Investigated datasets and estimation procedure

As for the real-world dataset, we use daily data from 10 market indexes.
Specifically, we consider: Standard & Poor 500 (SPX); IPC Mexico (MXX); BOlsa
de Valores do Estado de São PAulo (BOVESPA); Cotation Assistée en Continu 40
(CAC); Korea Composite Stock Price Index (KOSPI); Iberia Index 35 (IBEX);
Nihon Keizai Shinbun 225 (NIKKEI); Hang Seng Index (HSI); Deutscher
Aktienindex 40 (DAX); Financial Times Stock Exchange 100 (FTSE). The dataset
spans from 01-07-1993 to 30-06-2023, and it is divided into 24 block folds,
each one seven-year long, obtained by rolling a seven-year-long window over
the whole period covered by the dataset and shifting it by one year time by
time. Then, each fold is split into training and test sets. We define the
training set as the first-part six-year-long time series, which is further
divided into five and one-year segments when a validation set is needed. The
test set covers the last year. It is worth highlighting that each fold is
considered an individual experiment. However, we aggregate the out-of-sample
results when discussing the performance in the next subsections.

####  5.2.2 Main results

$\bm{\theta=0.05}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0022 | 0.0045 | 0.0025 | 0.0033 | 0.0032 | 0.0031 | 0.0025 | 0.0024 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.847 | 0.901 | 1.25 | 0.985 | 1.01 | 1.074 | 1.061 | 1.05 | 0.998 | 0.797
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.0023 | 0.006 | 0.0024 | 0.0035 | 0.0032 | 0.0031 | 0.0028 | 0.0025 | 0.0017
$\mathcal{L}_{P}^{\theta}$ | 0.844 | 0.908 | 1.286 | 0.978 | 1.054 | 1.039 | 1.092 | 1.061 | 1.005 | 0.792
BCGNS | $\mathcal{L}_{B}^{\theta}$ |  0.0111∗ |  0.0091∗ | 0.0117 |  0.0224∗ | 0.0071 |  0.0101∗ |  0.0159∗ |  0.0136∗ |  0.0158∗ |  0.018∗
$\mathcal{L}_{P}^{\theta}$ |  2.44∗ |  2.571∗ |  8.451∗ |  5.984∗ |  2.204∗ |  2.168∗ |  4.249∗ |  3.126∗ |  15.583∗ |  2.516∗
K-QRNN | $\mathcal{L}_{B}^{\theta}$ | 0.0049 | 0.0026 | 0.0064 | 0.0048 | 0.0041 | 0.0048 | 0.0048 | 0.0038 | 0.0045 | 0.0037
$\mathcal{L}_{P}^{\theta}$ | 1.152 | 1.071 | 1.427 | 1.303 | 1.218 | 1.26 | 1.23 | 1.282 |  1.339∗ | 1.104
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.0057 | 0.0024 | 0.0036 | 0.0043 | 0.0028 | 0.0033 | 0.0038 | 0.0035 | 0.0031 | 0.0022
$\mathcal{L}_{P}^{\theta}$ | 1.216 | 0.985 | 1.31 |  1.72∗ | 1.154 | 1.296 | 1.166 | 1.108 | 1.071 | 0.944
GAS2 | $\mathcal{L}_{B}^{\theta}$ | 0.0057 |  0.0183∗ | 0.005 |  0.0268∗ |  0.05∗ | 0.0075 | 0.0033 |  0.0544∗ |  0.0099∗ |  0.0086∗
$\mathcal{L}_{P}^{\theta}$ |  1.313∗ |  1.718∗ |  1.656∗ |  1.672∗ |  1.775∗ |  2.065∗ |  1.568∗ |  2.484∗ |  1.957∗ |  1.484∗

$\bm{\theta=0.025}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0047 | 0.0191 | 0.0058 | 0.0091 | 0.0075 | 0.0072 | 0.0048 | 0.0056 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.064 | 1.096 | 1.475 | 1.148 | 1.275 | 1.191 | 1.243 | 1.255 | 1.169 | 0.965
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.0054 | 0.0057 | 0.0188 | 0.0057 | 0.0093 | 0.0072 | 0.0068 | 0.0103 | 0.0057 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.06 | 1.119 | 1.487 | 1.151 | 1.27 | 1.191 | 1.268 | 1.261 | 1.172 | 0.96
BCGNS | $\mathcal{L}_{B}^{\theta}$ |  0.0149∗ | 0.0106 | 0.0228 |  0.0164∗ | 0.017 | 0.0171 |  0.0187∗ |  0.0158∗ |  0.0159∗ |  0.0128∗
$\mathcal{L}_{P}^{\theta}$ |  5.263∗ |  2.218∗ |  2.536∗ |  2.023∗ |  2.523∗ |  1.97∗ |  4.342∗ |  1.84∗ |  3.0∗ |  2.807∗
K-QRNN | $\mathcal{L}_{B}^{\theta}$ |  0.016∗ | 0.0061 | 0.0212 | 0.0133 | 0.0116 | 0.014 | 0.0133 | 0.0103 | 0.0112 |  0.0114∗
$\mathcal{L}_{P}^{\theta}$ |  1.466∗ | 1.294 | 1.693 |  1.518∗ | 1.492 | 1.492 | 1.467 | 1.517 |  1.507∗ |  1.379∗
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.0068 | 0.0043 | 0.018 | 0.0054 | 0.0117 | 0.0204 | 0.0129 | 0.0085 |  0.0467∗ | 0.0055
$\mathcal{L}_{P}^{\theta}$ | 1.25 | 1.328 |  1.886∗ | 1.415 | 1.401 | 1.376 | 1.505 | 1.443 |  1.636∗ |  1.386∗
GAS2 | $\mathcal{L}_{B}^{\theta}$ |  0.0232∗ |  0.0242∗ | 0.0382 |  0.0417∗ |  0.0637∗ | 0.0114 |  0.0744∗ |  0.0238∗ |  0.0225∗ |  0.0113∗
$\mathcal{L}_{P}^{\theta}$ |  2.65∗ |  2.288∗ | 1.696 |  1.975∗ |  2.451∗ |  1.551∗ |  1.887∗ |  2.303∗ |  1.834∗ |  2.187∗

$\bm{\theta=0.01}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0224 | 0.0188 | 0.0451 | 0.0273 | 0.0353 | 0.026 | 0.0236 | 0.0212 | 0.0333 | 0.011
$\mathcal{L}_{P}^{\theta}$ | 1.427 | 1.376 | 1.652 | 1.394 | 1.493 | 1.389 | 1.482 | 1.476 | 1.461 | 1.191
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.0286 | 0.0191 | 0.0441 | 0.0225 | 0.0364 | 0.0249 | 0.0214 | 0.0153 | 0.0331 | 0.0095
$\mathcal{L}_{P}^{\theta}$ | 1.554 | 1.412 | 1.756 | 1.387 | 1.541 | 1.394 | 1.499 | 1.433 | 1.441 | 1.172
BCGNS | $\mathcal{L}_{B}^{\theta}$ |  0.0743∗ | 0.0341 | 0.1062 | 0.044 | 0.0777 | 0.0688 |  0.067∗ |  0.0628∗ | 0.0487 |  0.0382∗
$\mathcal{L}_{P}^{\theta}$ |  3.48∗ |  2.262∗ |  2.335∗ |  1.993∗ |  2.539∗ |  2.161∗ |  2.603∗ |  2.438∗ |  2.222∗ |  2.131∗
K-QRNN | $\mathcal{L}_{B}^{\theta}$ |  0.0631∗ | 0.0131 | 0.0941 | 0.0453 | 0.0396 | 0.0504 | 0.0426 | 0.037 | 0.0355 |  0.0378∗
$\mathcal{L}_{P}^{\theta}$ | 1.815 | 1.508 |  2.049∗ | 1.77 | 1.782 | 1.735 | 1.752 | 1.811 | 1.717 |  1.659∗
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.0248 | 0.0131 | 0.0496 | 0.0704 | 0.0921 | 0.0304 | 0.0325 | 0.049 | 0.1029 |  0.0675∗
$\mathcal{L}_{P}^{\theta}$ | 1.443 | 1.78 |  2.046∗ |  2.103∗ | 1.774 |  1.82∗ |  2.203∗ |  2.246∗ | 1.612 | 1.555
GAS2 | $\mathcal{L}_{B}^{\theta}$ | 0.0479 | 0.034 | 0.0131 | 0.0211 | 0.0602 | 0.0533 | 0.0383 | 0.0319 | 0.0387 |  0.0385∗
$\mathcal{L}_{P}^{\theta}$ |  2.312∗ | 1.84 |  2.026∗ |  1.903∗ |  2.297∗ |  1.766∗ |  2.235∗ |  1.962∗ | 1.9 |  2.153∗

Table 2: Out-of-sample loss. It is shown the mean value across the 24 folds.
$\mathcal{L}^{\theta}_{B}$ represents the loss used in [5] (Equation 5), while
$\mathcal{L}^{\theta}_{P}$ denotes the loss in [39] (Equation (6)). Different
confidence levels are studied: $\theta=0.05,0.025,0.01$. The best result is
highlighted in bold, and the second best is underlined. The ∗ is used for
those predictors whose mean loss is outside the CAESar confidence interval
with one standard deviation.

The results for the indexes dataset at different levels of $\theta$ are shown
in Table 2. The table reports the mean value of the out-of-sample loss across
the 24 folds, considering both losses in Eqs. (5) and (6)) used in [5] and
[39]. For the latter, we have used the percentage returns to obtain values
coherent with [39]. As the table indicates, better results are typically
achieved using the CAESar and K-CAViaR approaches. Specifically, CAESar
overcomes the competitor more frequently, it is outside the top two algorithms
fewer times, and it displays better performances when using the Patton loss
$\mathcal{L}_{P}^{\theta}$ for evaluation. However, for some time series, the
GAS1 approach also performs well. Widely, GAS1 seems to very often outperforms
GAS2. This is in line with the findings in [39]. Generally, all models show a
decline in performance when $\theta$ is low. This trend is coherent with the
ES definition. Indeed, it averages values in the tail, and we may assume that
the deep tail is caused by very extreme and rare (hence unpredictable) events.
However, they are softened by the less severe events near the central body of
the distribution. As the confidence value $\theta$ decreases, the contribution
of these mitigating values vanishes, making the prediction problem even
tougher. This leads to worsened performance. Similar results are also observed
when working with stock data, as shown in the Supplementary Material D.
As previously discussed in Section 3, it is interesting to understand what
happens to the quantile estimate when we update the $\bm{\hat{q}}$ parameters
in the third step of CAESar optimization procedure. Specifically, a
degradation in the VaR forecast capability could undermine CAESar’s
reliability. Actually, we find out this is not the case. Specifically, we
study the Pinball loss (Equation (1)) of all the competing methods. For each
confidence level $\theta$, asset, and fold, we compute the mean loss
difference between the competitor and CAESar. Figure 2 shows the histogram of
these differences. Specifically, for every competitor $Comp$, the histogram
shows the distribution of the vector
$\left[\mathcal{L}_{q}^{\theta}(\bm{\hat{q}}^{Comp}_{\theta,a,f},\bm{y}_{a,f})-\mathcal{L}_{q}^{\theta}(\bm{\hat{q}}^{CAESar}_{\theta,a,f},\bm{y}_{a,f})\right]_{\theta\in\Theta,\
a\in A,\ f\in F}$, where $\theta\in\Theta=\\{0.05,0.025,0.01\\}$ represents
the confidence level, $a\in A=\\{SPX,MXX,\cdots\\}$ represents the asset, and
$f\in F=\\{1,\cdots,24\\}$ is for the fold. Finally,
$\bm{\hat{q}}^{CAESar}_{\theta,a,f}$ is for the CAESar prediction, while
$\bm{\hat{q}}^{Comp}_{\theta,a,f}$ is for the competitor forecast.

Figure 2: Histogram of the Pinball loss difference. CAESar is compared with
all the competitors by means of the quantile estimate. Specifically, the
difference between the competitor and CAESar’s loss is plotted. The blue bars
in the histogram correspond to negative values, which represent the case when
the competitor outperforms CAESar. On the other hand, the orange bars are
associated with positive values, which means that CAESar overcomes the
competitor. Finally, the x-axis is scaled according to the magnitude of the
loss difference, and the y-axis is on a logarithmic scale to improve the
readability.

Three qualitative conclusions can be drawn. The first one is that conditional
autoregressive models overcome the competitors even when looking at the
quantile estimate. Then, neural network models (BCGNS and K-QRNN) outperform
GAS models. Finally, keep in mind that, by construction, the CAESar quantile
estimate before the joint updating step is exactly the K-CAViaR quantile
estimate (as both of them are simply the CAViaR estimate at the confidence
level $\theta$). Thus, from the upper-left plot, we can observe how the VaR
forecast by CAESar is not degraded by the update in the $\bm{\hat{q}}$
parameters. On the contrary, it is even slightly improved.
Finally, Figure 3 shows the average computational times of the different
methods. The times are computed on a Ubuntu notebook equipped with a 13th Gen
Intel Core™ i7-13620H x 16 processor and a NVIDIA GeForce RTX 4060 video card.
The CAESar model is parallelized over 3 processes, one for each random vector
of coefficients used as the starting point in the optimization step (more
details in Supplementary Material B). The K-CAViaR approach exploits 10
processes, one for each CAViaR model corresponding to a quantile $\theta_{j}$.
As for the GAS1 and GAS2 models, we initialize their parameters to those
provided by the authors, as their performance is highly dependent on the
choice of initial coefficients (as pointed out in [39]). Therefore,
multiprocessing is not necessary. Finally, the K-QRNN and BCGNS models are
trained on the GPU.

Figure 3: Comparison of the computational time required for each algorithm.
The plot is in log scale. The training time is measured in seconds. Every box
summarizes a 720-dimensional vector containing 10 assets, 24 folds, and 3
$\theta$ values.

###  5.3 Statistical Tests - Direct Approximation

We use the statistical tests previously described to further evaluate the
proposed approach’s performance. First, we utilize those that evaluate the
goodness of the ES fit, namely the McNeil and Frey test, the $Z_{1}$ statistic
from [1], and the $Z_{2}$ statistic. It is worth to underline that these tests
are performed asset by asset, fold by fold. So, for every asset, the test is
performed 24 times, for a total of 240 iterations. Table 3 shows the ratio of
rejection of the null hypothesis. The lower the ratio, the better the model.
The confidence level of each test is set to 0.05.

| $\bm{\theta=0.05}$ | $\bm{\theta=0.025}$ | $\bm{\theta=0.01}$
---|---|---|---
| $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$ | $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$ | $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$
CAESar | 0.09 | 0.14 | 0.18 | 0.13 | 0.16 | 0.11 | 0.31 | 0.32 | 0.19
K-CAViaR | 0.2 | 0.16 | 0.21 | 0.19 | 0.19 | 0.14 | 0.4 | 0.4 | 0.22
BCGNS | 0.25 | 0.1 | 0.06 | 0.3 | 0.19 | 0.13 | 0.55 | 0.4 | 0.32
K-QRNN | 0.43 | 0.37 | 0.48 | 0.53 | 0.52 | 0.44 | 0.71 | 0.66 | 0.59
GAS1 | 0.1 | 0.34 | 0.16 | 0.2 | 0.36 | 0.15 | 0.49 | 0.5 | 0.27
GAS2 | 0.26 | 0.49 | 0.57 | 0.42 | 0.52 | 0.41 | 0.62 | 0.62 | 0.47
Table 3: Direct Approximation Tests - Results from the McNeil and Frey test
[34] (MNF) and the Acerbi and Szekely [1], both with the Z1 (AS1) and Z2
statistics (AS2). A separate test is conducted for each asset and each fold,
focusing only on the out-of-sample predictions. Each cell value represents the
ratio of null hypothesis rejections (with a threshold of 0.05). The best
result (i.e. the lowest number of rejections) is in bold, and the second to
best is underlined.

For large values of $\theta$, several models perform well. Specifically, for
$\theta=0.05$, the three top-performing algorithms are CAESar, BCGNS, and
GAS1. The former shows superior performances in the McNeil and Frey (MNF) and
the $Z_{1}$ (AS1) tests. On the other side, BCGNS obtain good performances in
both the Acerbi and Szekely statistics. Instead, the GAS1 model performs well
in the McNeil and Frey and $Z_{2}$ tests (AS2). For $\theta=0.025$, the
situation for BCGNS is almost the same. Instead, CAESar achieves good results
in all the tests, while GAS1 shows a degradation of its performance and
performs well only in the MNF test. Instead, for $\theta=0.01$, the best
results on both the tests are obtained almost always by CAESar and the second
to best by the K-CAViaR. As pointed out in Section 4.1, the discrepancy
between the tests’ results can be explained by taking in mind the empirical
study carried out in [13]. They compare several direct approximation tests
applied to different time series lengths in order to have a picture of their
finite-sample properties. The ultimate result, when using the length the most
similar to that of our folds, is that the AS1 test is the most stable both
regarding size and power, followed by the MNF test. Finally, observe that the
number of rejections follows the same previously discussed pattern: model
performance tends to degrade when $\theta$ is small. However, as the relative
performances of CAESar improve (i.e. its results look better than the ones of
the competitors), we can conclude that it experiences less degradation in
performance. This is coherent also with the results in Table 2.

###  5.4 Statistical Tests - Diebold-Mariano for Model Comparison

In this subsection, we evaluate the proposed ES estimators using the
previously discussed Diebold-Mariano test, which is based on the loss
function. As we use two losses in this work, the test is repeated twice.
Furthermore, it is applied once for each asset and fold. The null hypothesis
is the statistical equivalence of the two ES estimators. Table 4 shows the
number of times this hypothesis is rejected in the Diebold-Mariano test (with
a confidence level equal to 0.05). The left value in each cell indicates the
number of times CAESar has statistically outperformed the competitor, while
the right value is the count of the competitor wins. For a quick overview, we
have only reported the sum of the number of rejections for each asset. The
results clearly show the supremacy of CAESar, which outperforms all the
competitors, by considering both losses and for all $\theta$ levels. Finally,
the Loss Difference, the Encompassing, and the corrected resampled t-test
previously described are shown in Supplementary Material C. The results are in
line with the Diebold-Mariano test.

|  | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---
$\bm{\theta=0.05}$ | $\mathcal{L}_{B}^{\theta}$ | 20 / 0 | 176 / 0 | 27 / 0 | 52 / 0 | 78 / 0
$\mathcal{L}_{P}^{\theta}$ | 28 / 15 | 163 / 0 | 134 / 3 | 74 / 2 | 128 / 5
$\bm{\theta=0.025}$ | $\mathcal{L}_{B}^{\theta}$ | 26 / 0 | 99 / 0 | 22 / 0 | 41 / 0 | 82 / 0
$\mathcal{L}_{P}^{\theta}$ | 29 / 24 | 133 / 2 | 134 / 3 | 102 / 5 | 126 / 0
$\bm{\theta=0.01}$ | $\mathcal{L}_{B}^{\theta}$ | 38 / 0 | 40 / 0 | 19 / 0 | 57 / 0 | 65 / 0
$\mathcal{L}_{P}^{\theta}$ | 54 / 15 | 70 / 8 | 122 / 2 | 114 / 6 | 119 / 4
Table 4: Model Comparison Tests - Diebold-Mariano test with confidence 0.05.
The results obtained across various assets, for a given probability level
$\theta$, have been aggregated by summing them. Each cell value is in the form
$g/b$, where $g$ represents the number of folds where CAESar statistically
outperforms the competitor, and $b$ is the contrary.

##  6 Conclusions

Forecasting the ES in a dynamic context is known to be a hard task for several
reasons. One of them is the intrinsic complexity of financial time series,
which are characterized by heteroskedasticity and fat tails, which complicates
the modelling and forecasting process. Moreover, the fact that ES is not
elicitable further adds to the challenge. This lack has historically slowed
down the research in this direction, at least until the finding that ES and
VaR are jointly elicitable. This breakthrough has opened new doors in the
field. Despite this progress, a noticeable gap still exists between the
theoretical works and the practical needs of the industry for reliable models.
Adding to these challenges are the issues associated with ES backtesting,
mainly driven by the lack of elicitability. Furthermore, the fact that ES is a
latent variable -meaning it is not directly observed from the data- adds
another layer of complexity to its backtesting. Thus, backtesting is still an
open point, and a universally recognized statistical procedure to assess
models’ performance in forecasting the ES of a conditional distribution is
missing.
We contribute to this discussion by proposing CAESar, a three-step approach
for predicting the ES at a given time, conditional on the information
available at the previous time. CAESar draws significant inspiration from the
CAViaR approach to quantile regression, a statistical technique used to
estimate the conditional quantile. Indeed, the specification for the
$\bm{\hat{e}}$ estimator is designed as the asymmetric slope with two
autoregressive addends: one for the previous step quantile and the other for
the previous step ES. This allows for a robust and dynamic estimation of risk.
The first step of CAESar involves obtaining the VaR forecast, $\bm{\hat{q}}$,
using CAViaR. Then, the $\bm{\hat{e}}$ coefficients are trained while the VaR
model remains fixed. Finally, the two estimators are jointly trained. In both
steps, a soft constraint is added to the loss function to ensure the
monotonicity condition - that is, the ES estimate has to be less than the
quantile forecast.
To assess the effectiveness of our proposal, we have conducted extensive
testing. Most of the effort has been put into the daily index dataset. We
collected a dataset comprising daily observations from 10 stock indexes over
the last 30 years. This dataset was divided into block folds, and several
models were compared. Specifically, three types of evaluation were performed.
Firstly, we have compared the model loss functions. Next, we have conducted
statistical tests to assess whether one model statistically outperforms
another based on out-of-sample losses. Then, we performed loss-agnostic tests
to directly evaluate the ES estimator (i.e., the coherence between the
approximated values and the ES definition). The main finding is that CAESar
outperforms the other models throughout the comparative analysis we present.
In particular, looking at the direct approximation tests’ rejection frequency,
CAESar almost always displays the lowest value. Moreover, the gap is far more
evident in the very tail of the returns’ distribution (low values of
$\theta$), when all other ES models start to fail. It is worth noting that
both BCGNS (that is, a neural network-based approach) and K-CAViaR (which
approximates the expectation in the ES definition with the sample mean of
quantile forecasting at different confidence levels) models sometimes display
good forecasting performance, but the results for these models are not stable
and strongly depend on the value of $\theta$. Additionally, we conducted
experiments on a dataset composed of stocks in the U.S. banking sector. Once
again, CAESar has on average superior performance.
There are two primary directions for further research. On the algorithmic
side, the next step would be to predict ES for multiple correlated assets.
This would involve updating CAESar to handle multivariate time series,
providing valuable insights into tail risk networks. Furthermore, moving from
the ES of single assets to that of a portfolio is key. In this context, CAESar
could serve as a tool to assist investors and risk managers in adjusting
portfolio weights to minimize overall tail risk. On the theoretical side, it
would be useful to gain a deeper understanding of the proposed estimator,
particularly its asymptotic properties. This could enhance our understanding
of the dynamics captured by CAESar.

## Data and Code Availability

The data utilized in this paper is publicly accessible. Specifically, index
data have been downloaded from three sources: Investing
https://www.investing.com/, Yahoo Finance https://finance.yahoo.com, and
MarketWatch https://www.marketwatch.com/. This data has been combined to
complete any gaps in the time series. As for the banking sector stocks, the
data has been obtained from Investing and Yahoo Finance. The CAESar code,
along with the codes for competing strategies (unless otherwise stated in the
paper) and model comparisons, are publicly accessible at the GitHub repository
https://github.com/fgt996/CAESar.

## Acknowledgments

We would like to thank the Quantitative Finance Workshop 2024. FL acknowledges
”SoBigData.it receives funding from European Union – NextGenerationEU –
National Recovery and Resilience Plan (Piano Nazionale di Ripresa e
Resilienza, PNRR) – Project: “SoBigData.it – Strengthening the Italian RI for
Social Mining and Big Data Analytics” – Prot. IR0000013 – Avviso n. 3264 del
28/12/2021.”

SUPPLEMENTARY MATERIAL

Supplementary Material:

Discussion about the specific functional form of our estimator; examination of
the parameters optimization scheme; further insights into the results obtained
by the statistical tests used to validate our proposal; analysis of the CAESar
performance on the stock dataset. (.pdf file)

## References

  * [1] C. Acerbi and B. Szekely.  Back-testing expected shortfall.  Risk, 27(11):76–81, 2014.
  * [2] C. Acerbi and D. Tasche.  On the coherence of expected shortfall.  Journal of Banking & Finance, 26(7):1487–1503, 2002.
  * [3] P. Artzner, F. Delbaen, J.-M. Eber, and D. Heath.  Coherent measures of risk.  Mathematical Finance, 9(3):203–228, 1999.
  * [4] D. Bams, G. Blanchard, and T. Lehnert.  Volatility measures and value-at-risk.  International Journal of Forecasting, 33(4):848–863, 2017.
  * [5] D. Barrera, S. Crépey, E. Gobet, H.-D. Nguyen, and B. Saadeddine.  Learning value-at-risk and expected shortfall.  arXiv preprint arXiv:2209.06476, 2022.
  * [6] F. Bellini and E. Di Bernardino.  Risk management with expectiles.  The European Journal of Finance, 23(6):487–506, 2017.
  * [7] J. Brehmer.  Elicitability and its application in risk management.  arXiv preprint arXiv:1707.09604, 2017.
  * [8] L. Carver.  Mooted var substitute cannot be back-tested, says top quant.  Risk, 8:184–201, 2013.
  * [9] R. Cont.  Empirical properties of asset returns: stylized facts and statistical issues.  Quantitative Finance, 1(2):223–236, 2001.
  * [10] D. Creal, S. J. Koopman, and A. Lucas.  Generalized autoregressive score models with applications.  Journal of Applied Econometrics, 28(5):777–795, 2013.
  * [11] A. Daouia, S. Girard, and G. Stupfler.  Expecthill estimation, extreme risk and heavy tails.  Journal of Econometrics, 221(1):97–117, 2021.
  * [12] G. De Rossi and A. Harvey.  Quantiles, expectiles and splines.  Journal of Econometrics, 152(2):179–185, 2009.
  * [13] K. Deng and J. Qiu.  Backtesting expected shortfall and beyond.  Quantitative Finance, 21(7):1109–1125, 2021.
  * [14] F. X. Diebold and R. S. Mariano.  Comparing predictive accuracy.  Journal of Business & Economic Statistics, 20(1):134–144, 2002\.
  * [15] T. Dimitriadis and R. Halbleib.  Realized quantiles.  Journal of Business & Economic Statistics, 40(3):1346–1361, 2022\.
  * [16] Z. Du, P. Pei, X. Wang, and T. Yang.  Powerful backtests for historical simulation expected shortfall models.  Journal of Business & Economic Statistics, 42(3):864–874, 2023\.
  * [17] S. Emmer, M. Kratz, and D. Tasche.  What is the best risk measure in practice? a comparison of standard measures.  Journal of Risk, 18(2):31–60, 2015.
  * [18] R. F. Engle and S. Manganelli.  Caviar: Conditional autoregressive value at risk by regression quantiles.  Journal of Business & Economic Statistics, 22(4):367–381, 2004\.
  * [19] M. Feridun and A. Özün.  Basel iv implementation: a review of the case of the european union.  Journal of Capital Markets Studies, 4(1):7–24, 2020.
  * [20] T. Fissler and J. F. Ziegel.  Higher order elicitability and osband’s principle.  The Annals of Statistics, 44(4):1680–1707, 2016.
  * [21] A.-P. Fortin, J.-G. Simonato, and G. Dionne.  Forecasting expected shortfall: Should we use a multivariate model for stock market factors?  International Journal of Forecasting, 39(1):314–331, 2023.
  * [22] F. Fuentes, R. Herrera, and A. Clements.  Forecasting extreme financial risk: A score-driven approach.  International Journal of Forecasting, 39(2):720–735, 2023.
  * [23] R. H. Gerlach, C. W. Chen, and N. Y. Chan.  Bayesian time-varying quantile forecasting for value-at-risk in financial markets.  Journal of Business & Economic Statistics, 29(4):481–492, 2011\.
  * [24] T. Gneiting.  Making and evaluating point forecasts.  Journal of the American Statistical Association, 106(494):746–762, 2011.
  * [25] H. Han, O. Linton, T. Oka, and Y.-J. Whang.  The cross-quantilogram: Measuring quantile dependence and testing directional predictability between time series.  Journal of Econometrics, 193(1):251–270, 2016.
  * [26] D. Harvey, S. Leybourne, and P. Newbold.  Testing the equality of prediction mean squared errors.  International Journal of Forecasting, 13(2):281–291, 1997.
  * [27] S. Jaimungal, S. M. Pesenti, Y. F. Saporito, and R. S. Targino.  Risk budgeting allocation for dynamic risk measures.  arXiv preprint arXiv:2305.11319, 2023.
  * [28] T. C. Johnson.  Volatility, momentum, and time-varying skewness in foreign exchange returns.  Journal of Business & Economic Statistics, 20(3):390–411, 2002\.
  * [29] G. Keilbar and W. Wang.  Modelling systemic risk using neural network quantile regression.  Empirical Economics, 62(1):93–118, 2022.
  * [30] T. Kışınbay.  The use of encompassing tests for forecast combinations.  Journal of Forecasting, 29(8):715–727, 2010.
  * [31] R. Koenker and G. Bassett Jr.  Regression quantiles.  Econometrica, 46(1):33–50, 1978.
  * [32] D. Kraft.  A software package for sequential quadratic programming.  Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt, 1988.
  * [33] M. Kratz, Y. H. Lok, and A. J. McNeil.  Multinomial var backtests: A simple implicit approach to backtesting expected shortfall.  Journal of Banking & Finance, 88:393–407, 2018.
  * [34] A. J. McNeil and R. Frey.  Estimation of tail-related risk measures for heteroscedastic financial time series: an extreme value approach.  Journal of Empirical Finance, 7(3-4):271–300, 2000.
  * [35] L. Merlo, L. Petrella, and V. Raponi.  Forecasting var and es using a joint quantile regression and its implications in portfolio allocation.  Journal of Banking & Finance, 133:106248, 2021.
  * [36] C. Nadeau and Y. Bengio.  Inference for the generalization error.  Machine Learning, 52(3):239–281, september 2003.
  * [37] J. Nicolau, P. M. Rodrigues, and M. Z. Stoykov.  Tail index estimation in the presence of covariates: Stock returns’ tail risk dynamics.  Journal of Econometrics, 235(2):2266–2284, 2023.
  * [38] N. D. Pankratova and N. G. Zrazhevska.  Method of dynamic var and cvar risk measures forecasting for long range dependent time series on the base of the heteroscedastic model.  Intelligent Control and Automation, 8(2):126–138, 2017.
  * [39] A. J. Patton, J. F. Ziegel, and R. Chen.  Dynamic semiparametric models for expected shortfall (and value-at-risk).  Journal of Econometrics, 211(2):388–413, 2019.
  * [40] R. T. Rockafellar, S. Uryasev, et al.  Optimization of conditional value-at-risk.  Journal of Risk, 2:21–42, 2000.
  * [41] J. W. Taylor.  Estimating value at risk and expected shortfall using expectiles.  Journal of Financial Econometrics, 6(2):231–252, 2008.
  * [42] J. W. Taylor.  Forecasting value at risk and expected shortfall using a semiparametric approach based on the asymmetric laplace distribution.  Journal of Business & Economic Statistics, 37(1):121–133, 2019\.
  * [43] R. J. Tibshirani and B. Efron.  An introduction to the bootstrap.  Monographs on Statistics and Applied Probability, 57(1):1–436, 1993\.
  * [44] S. Weber.  Distribution-invariant risk measures, information, and dynamic consistency.  Mathematical Finance, 16(2):419–441, 2006.
  * [45] A. Wilhelmsson.  Value at risk with time varying variance, skewness and kurtosis—the nig-acd model.  The Econometrics Journal, 12(1):82–104, 2009.

CAESar:
Conditional Autoregressive Expected Shortfall
Supplementary Material
Federico Gatta
Classe di Scienze, Scuola Normale Superiore
and
Fabrizio Lillo
Dipartimento di Matematica, Università di Bologna & Scuola Normale Superiore
and
Piero Mazzarisi
Dipartimento di Economia Politica e Statistica, Università di Siena
July 4, 2024 0

##  Appendix A Supplementary Material - Model Specifications

In Section 3, we have discussed a general formulation for the CAESar
regressor, as we have done for the CAViaR in Section 2. Nonetheless, we
specifically propose to use the (1,1) model with Asymmetric Slope (AS) for
$y_{t}$ transformation. Indeed, as it is a popular choice for CAViaR in the
literature, we expect it to work very well also for CAESar. The reasons behind
such a success are twofold. As for the $\bm{f}$ and $\bm{g}$ functions
specification, Asymmetric Slope is the only one (among those proposed in [18])
that can efficiently capture the asymmetry in market movements. Indeed, the
other popular alternatives, namely the Symmetric Absolute Value (SAV) and the
indirect GARCH specifications, are not able to do this, typically resulting in
poorer performance, or at least not overcoming AS. Specifically, the SAV (1,1)
specification assumes symmetry in the information contained in positive and
negative market returns, thus resulting in:

$\hat{q}_{t}=\beta_{0}+\beta_{1}|y_{t-1}|+\beta_{2}\hat{q}_{t-1}\quad
and\quad\hat{e}_{t}=\gamma_{0}+\gamma_{1}|y_{t-1}|+\gamma_{2}\hat{q}_{t-1}+\gamma_{3}\hat{e}_{t-1}$
(22)

On the other hand, there is the GARCH (1,1):

$\hat{q}_{t}=-\sqrt{\beta_{0}+\beta_{1}y_{t-1}^{2}+\beta_{2}\hat{q}_{t-1}}\quad
and\quad\hat{e}_{t}=-\sqrt{\gamma_{0}+\gamma_{1}y_{t-1}^{2}+\gamma_{2}\hat{q}_{t-1}+\gamma_{3}e_{t-1}^{2}}$
(23)

Observe that according to the formulation in Equation (2), with a slight abuse
of notation, we are applying this specification to compute the squared target.
To corroborate our choice, we repeat the experiments using both the AS, SAV,
and GARCH (1,1) specifications. The results are in Table 5 and show that the
other specifications do not significantly outperform AS.

$\bm{\theta=0.05}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (AS) | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0022 | 0.0045 | 0.0025 | 0.0033 | 0.0032 | 0.0031 | 0.0025 | 0.0024 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.847 | 0.901 | 1.25 | 0.985 | 1.01 | 1.074 | 1.061 | 1.05 | 0.998 | 0.797
CAESar (SAV) | $\mathcal{L}_{B}^{\theta}$ | 0.0028 | 0.0025 | 0.0049 | 0.0032 | 0.0037 | 0.0035 | 0.0032 | 0.003 | 0.0028 | 0.0021
$\mathcal{L}_{P}^{\theta}$ | 0.994 | 0.923 | 1.266 | 1.026 | 1.071 | 1.057 | 1.097 | 1.055 | 1.036 | 0.831
CAESar (G) | $\mathcal{L}_{B}^{\theta}$ | 0.0055 | 0.0022 | 0.0058 | 0.0023 | 0.005 | 0.0025 | 0.0026 | 0.0037 | 0.0035 | 0.0024
$\mathcal{L}_{P}^{\theta}$ | 1.025 | 1.003 | 1.356 | 1.261 | 1.396 | 1.073 | 1.113 | 1.393 | 1.152 | 1.15

$\bm{\theta=0.025}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (AS) | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0047 | 0.0191 | 0.0058 | 0.0091 | 0.0075 | 0.0072 | 0.0048 | 0.0056 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.064 | 1.096 | 1.475 | 1.148 | 1.275 | 1.191 | 1.243 | 1.255 | 1.169 | 0.965
CAESar (SAV) | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0049 | 0.0194 | 0.0071 | 0.011 | 0.0095 | 0.0068 | 0.0056 | 0.0065 | 0.004
$\mathcal{L}_{P}^{\theta}$ | 1.061 | 1.12 | 1.549 | 1.19 | 1.278 | 1.248 | 1.277 | 1.222 | 1.201 | 1.0
CAESar (G) | $\mathcal{L}_{B}^{\theta}$ | 0.014 | 0.005 | 0.0464 | 0.0148 | 0.0111 | 0.0065 | 0.0059 | 0.0052 | 0.0054 | 0.0031
$\mathcal{L}_{P}^{\theta}$ | 1.376 | 1.165 | 2.113 | 1.468 | 1.592 | 1.196 | 1.289 | 1.2 | 1.258 | 1.137

$\bm{\theta=0.01}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (AS) | $\mathcal{L}_{B}^{\theta}$ | 0.0224 | 0.0188 | 0.0451 | 0.0273 | 0.0353 | 0.026 | 0.0236 | 0.0212 | 0.0333 | 0.011
$\mathcal{L}_{P}^{\theta}$ | 1.427 | 1.376 | 1.652 | 1.394 | 1.493 | 1.389 | 1.482 | 1.476 | 1.461 | 1.191
CAESar (SAV) | $\mathcal{L}_{B}^{\theta}$ | 0.018 | 0.0135 | 0.0882 | 0.0306 | 0.0378 | 0.0328 | 0.0202 | 0.0178 | 0.0251 | 0.0123
$\mathcal{L}_{P}^{\theta}$ | 1.334 | 1.391 | 1.862 | 1.432 | 2.078 | 1.429 | 1.474 | 1.44 | 1.509 | 1.195
CAESar (G) | $\mathcal{L}_{B}^{\theta}$ | 0.0563 | 0.0172 | 0.061 | 0.1043 | 0.0251 | 0.0494 | 0.0207 | 0.0518 | 0.0206 | 0.0135
$\mathcal{L}_{P}^{\theta}$ | 1.613 | 1.501 | 1.706 | 1.984 | 1.62 | 1.609 | 1.5 | 2.018 | 1.478 | 1.335

Table 5: Comparison of the three specifications (SAV, AS, and GARCH) with one
lag for both the asset value and previous quantile and expected shortfall on
the index dataset. The same specification is used both for the CAViaR
regressor in the first step and for the $\bm{\hat{r}}$ and $\bm{\hat{e}}$
models in the last steps. The mean value across the 24 folds is displayed. The
best result is in bold, the second to best is underlined.

As for the SAV and GARCH specifications, the former seems to perform slightly
better for $\theta=0.01$, while the latter works better with higher $\theta$.
However, none of them can consistently overcome the AS specification, which
shows superior performances for all $\theta$ and almost all the considered
assets, being always able to give performances close to the best one.
Another point when choosing the CAESar model is related to the lags to
consider. The choice is related to the observations lag $p$ and the previous
estimates lag $u$. We suggest to use $p=u=1$. Indeed, it is the cheapest.
Moreover, it avoids the overfitting risk, and its estimates are often the most
accurate. We empirically show this in Table 6, where the AS-(1,1) is compared
with the AS-(1,2), AS-(2,1), and AS(2,2) specifications.

$\bm{\theta=0.05}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (1,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0022 | 0.0045 | 0.0025 | 0.0033 | 0.0032 | 0.0031 | 0.0025 | 0.0024 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.847 | 0.901 | 1.25 | 0.985 | 1.01 | 1.074 | 1.061 | 1.05 | 0.998 | 0.797
$C.T.$ | 1.3 | 2.5 | 2.0 | 1.5 | 1.3 | 2.0 | 1.9 | 1.4 | 1.3 | 1.4
CAESar (1,2) | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0024 | 0.0062 | 0.0026 | 0.0035 | 0.0034 | 0.0032 | 0.0027 | 0.0027 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.867 | 1.009 | 1.29 | 1.073 | 1.115 | 1.114 | 1.099 | 1.405 | 1.013 | 0.825
$C.T.$ | 5.0 | 5.1 | 4.4 | 4.4 | 2.8 | 5.6 | 3.1 | 4.9 | 4.2 | 4.3
CAESar (2,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.002 | 0.0047 | 0.0026 | 0.0032 | 0.0035 | 0.0032 | 0.0025 | 0.0048 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.846 | 0.91 | 1.247 | 0.984 | 1.046 | 1.037 | 1.086 | 1.035 | 1.488 | 0.805
$C.T.$ | 4.8 | 7.1 | 7.3 | 5.3 | 6.3 | 7.8 | 4.5 | 6.2 | 5.4 | 3.3
CAESar (2,2) | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0023 | 0.0062 | 0.0027 | 0.0037 | 0.0039 | 0.0033 | 0.0028 | 0.0025 | 0.002
$\mathcal{L}_{P}^{\theta}$ | 0.882 | 0.925 | 1.324 | 1.011 | 1.056 | 1.32 | 1.105 | 1.059 | 1.097 | 0.898
$C.T.$ | 14.2 | 16.8 | 16.0 | 11.3 | 10.7 | 13.9 | 10.5 | 14.6 | 11.9 | 8.7

$\bm{\theta=0.025}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (1,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0047 | 0.0191 | 0.0058 | 0.0091 | 0.0075 | 0.0072 | 0.0048 | 0.0056 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.064 | 1.096 | 1.475 | 1.148 | 1.275 | 1.191 | 1.243 | 1.255 | 1.169 | 0.965
$C.T.$ | 1.2 | 1.6 | 1.3 | 1.2 | 1.3 | 1.6 | 1.4 | 1.5 | 1.0 | 1.2
CAESar (1,2) | $\mathcal{L}_{B}^{\theta}$ | 0.007 | 0.0051 | 0.0127 | 0.0067 | 0.0082 | 0.0331 | 0.0079 | 0.0054 | 0.0081 | 0.008
$\mathcal{L}_{P}^{\theta}$ | 1.193 | 1.133 | 1.442 | 1.17 | 1.211 | 1.611 | 1.562 | 1.227 | 6.456 | 2.296
$C.T.$ | 3.7 | 3.5 | 3.5 | 3.2 | 4.8 | 3.9 | 3.5 | 3.0 | 3.6 | 5.1
CAESar (2,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0053 | 0.0042 | 0.0156 | 0.0063 | 0.0076 | 0.0073 | 0.0068 | 0.0044 | 0.0061 | 0.0035
$\mathcal{L}_{P}^{\theta}$ | 1.079 | 1.061 | 1.451 | 1.419 | 1.221 | 1.194 | 1.272 | 1.195 | 1.168 | 0.988
$C.T.$ | 6.8 | 3.3 | 7.4 | 5.4 | 4.7 | 8.9 | 4.8 | 5.7 | 4.6 | 3.8
CAESar (2,2) | $\mathcal{L}_{B}^{\theta}$ | 0.0086 | 0.0049 | 0.0131 | 0.0055 | 0.012 | 0.0083 | 0.0069 | 0.0057 | 0.0063 | 0.004
$\mathcal{L}_{P}^{\theta}$ | 3.891 | 1.083 | 1.459 | 1.206 | 1.413 | 1.256 | 1.283 | 1.321 | 1.224 | 1.027
$C.T.$ | 8.5 | 11.9 | 11.7 | 11.0 | 9.1 | 15.2 | 8.8 | 11.3 | 8.1 | 8.6

$\bm{\theta=0.01}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar (1,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0224 | 0.0188 | 0.0451 | 0.0273 | 0.0353 | 0.026 | 0.0236 | 0.0212 | 0.0333 | 0.011
$\mathcal{L}_{P}^{\theta}$ | 1.427 | 1.376 | 1.652 | 1.394 | 1.493 | 1.389 | 1.482 | 1.476 | 1.461 | 1.191
$C.T.$ | 1.8 | 1.3 | 1.3 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 0.9 | 1.2
CAESar (1,2) | $\mathcal{L}_{B}^{\theta}$ | 0.0239 | 0.0116 | 0.1044 | 0.0212 | 0.0278 | 0.0541 | 0.0267 | 0.0179 | 00256 | 0.0166
$\mathcal{L}_{P}^{\theta}$ | 1.473 | 1.352 | 1.893 | 1.373 | 1.619 | 1.564 | 1.521 | 1.551 | 1.394 | 1.374
$C.T.$ | 3.4 | 3.2 | 3.0 | 5.0 | 3.8 | 4.1 | 3.1 | 3.1 | 2.6 | 3.5
CAESar (2,1) | $\mathcal{L}_{B}^{\theta}$ | 0.0657 | 0.0129 | 0.043 | 0.0403 | 0.0449 | 0.0231 | 0.0251 | 0.0154 | 0.0306 | 0.0128
$\mathcal{L}_{P}^{\theta}$ | 1.74 | 1.338 | 1.674 | 1.56 | 1.658 | 1.411 | 1.495 | 1.434 | 1.414 | 1.233
$C.T.$ | 4.7 | 4.5 | 6.1 | 3.1 | 6.0 | 6.8 | 4.2 | 4.3 | 3.3 | 3.5
CAESar (2,2) | $\mathcal{L}_{B}^{\theta}$ | 0.0336 | 0.0158 | 0.0944 | 0.0238 | 0.0293 | 0.029 | 0.0314 | 0.0251 | 0.031 | 0.016
$\mathcal{L}_{P}^{\theta}$ | 1.555 | 1.415 | 1.881 | 1.48 | 1.57 | 1.46 | 1.63 | 20.276 | 2.127 | 1.379
$C.T.$ | 7.0 | 8.0 | 9.2 | 10.8 | 10.0 | 10.5 | 8.3 | 5.8 | 6.4 | 5.9

Table 6: Comparison between different numbers of lags for the AS
specification. The specifications are compared by means of the
$\mathcal{L}_{B}^{\theta}$ loss, the $\mathcal{L}_{B}^{\theta}$ loss, and the
computational time ($C.T.$). The best result is in bold, and the second to
best is underlined.

As shown in the table, the (1,1) model is by far the fastest. Furthermore, its
performances are often at the top. Indeed, for $\theta=0.05$, even when it is
not the best specification, its performances are close to the best one.
Instead, for lower $\theta$, sometimes it is clearly overcome by the (1,2) or,
more frequently, the (2,1) specifications. However, if we look at the mean
behaviour, it looks pointless to consider over-complicated models.
Interestingly, the (2,2) model exhibits the worst performance. This can be due
to both misspecification or issues in the parameters optimization process due
to the increased complexity of the model.
Finally, it is worth observe we have taken for granted the cross terms in the
$\hat{q}_{t}$ and $\hat{e}_{t}$ estimators in Equation (11). That is, we have
assumed that $\hat{q}_{t}$ depends on $\hat{e}_{t-1}$ and $\hat{e}_{t}$ on
$\hat{q}_{t-1}$. To verify this assumption, we compare CAESar with its No
Cross version. The parameters of the latter are optimized by only looking at
the $\mathcal{L}_{P}^{\theta}$ loss, as the optimization of
$\mathcal{L}_{B}^{\theta}$ would necessarily require $\hat{e}_{t}$ dependent
on $\hat{q}_{t-1}$. The result is shown in Table 7.

$\bm{\theta=0.05}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0022 | 0.0045 | 0.0025 | 0.0033 | 0.0032 | 0.0031 | 0.0025 | 0.0024 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.847 | 0.901 | 1.25 | 0.985 | 1.01 | 1.074 | 1.061 | 1.05 | 0.998 | 0.797
No Cross | $\mathcal{L}_{B}^{\theta}$ | 0.0021 | 0.0023 | 0.006 | 0.0024 | 0.0034 | 0.0033 | 0.0032 | 0.0026 | 0.0023 | 0.0018
$\mathcal{L}_{P}^{\theta}$ | 0.862 | 0.892 | 1.266 | 0.988 | 1.049 | 1.026 | 1.085 | 1.023 | 1.003 | 0.793

$\bm{\theta=0.025}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0047 | 0.0191 | 0.0058 | 0.0091 | 0.0075 | 0.0072 | 0.0048 | 0.0056 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.064 | 1.096 | 1.475 | 1.148 | 1.275 | 1.191 | 1.243 | 1.255 | 1.169 | 0.965
No Cross | $\mathcal{L}_{B}^{\theta}$ | 0.0062 | 0.0051 | 0.0187 | 0.0057 | 0.0093 | 0.0071 | 0.0069 | 0.0126 | 0.0057 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.174 | 1.133 | 1.489 | 1.148 | 3.054 | 1.196 | 1.258 | 1.256 | 1.173 | 0.985

$\bm{\theta=0.01}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0224 | 0.0188 | 0.0451 | 0.0273 | 0.0353 | 0.026 | 0.0236 | 0.0212 | 0.0333 | 0.011
$\mathcal{L}_{P}^{\theta}$ | 1.427 | 1.376 | 1.652 | 1.394 | 1.493 | 1.389 | 1.482 | 1.476 | 1.461 | 1.191
No Cross | $\mathcal{L}_{B}^{\theta}$ | 0.0289 | 0.0128 | 0.0592 | 0.0225 | 0.036 | 0.0246 | 0.0245 | 0.0175 | 0.0246 | 0.0096
$\mathcal{L}_{P}^{\theta}$ | 1.682 | 1.328 | 1.847 | 1.404 | 1.672 | 1.402 | 1.504 | 1.471 | 1.416 | 1.228

Table 7: Comparison of CAESar with AS-(1,1) specification and the No Cross
specification, i.e., that obtained from Equation 15 by forcing
$\beta_{4}=\gamma_{3}=0$, that means no cross dependency between the VaR and
ES estimators. The best result is in bold.

Even if the results are not completely unidirectional, the table still shows
that the specification that keeps into account cross terms roughly outperforms
the other. This suggests that the intrinsic relationship between VaR and ES is
not merely related to elicitability, as information on the previous estimate
of one of them could improve the prediction for the other, suggesting a kind
of Granger causality between these variables.

##  Appendix B Supplementary Material - Optimization Scheme

In this Supplementary Material section, we discuss the optimization scheme
developed to learn the CAESar parameters. First, we outline the optimization
routine. Then, the choice of the loss functions is discussed. CAESar’s first
step involves obtaining the CAViaR quantile estimator. This is accomplished
using the code available on the personal website of its
author111simonemanganelli.org. As for the second step, we follow this
guideline:

  1. 1.

The $\hat{e}_{0}$ estimate, which will serve as the starting point in the
CAESar loop, is obtained using the empirical ES computed over the first 10% of
points in the training set.

  2. 2.

As for the initial guess of the coefficients vector
$\bm{\gamma}=(\gamma_{0},\cdots,\gamma_{4})$, 100 random initializations are
drawn from uniform (between -1 and 1) and standard normal distributions. Then,
they are sorted according to the
$\mathcal{L}^{\theta}_{r}(\bm{\hat{r}},\bm{y};\bm{VaR}(\theta))$ loss, and
only the top three are kept.

  3. 3.

Each of the three $\bm{\gamma}$ initial guesses is used as the starting point
for optimization. Specifically, a series of six sequential optimizations are
performed. Sequential here means the output of one optimization is used as the
starting point for the next. The optimizer used is the Sequential Least
SQuares Programming [32]. Note that three processes are simultaneously run,
one for each starting $\bm{\gamma}$ choice.

After the last point, we have three optimal $\bm{\gamma}$, each one obtained
from a different starting point. They are sorted by means of the
$\mathcal{L}^{\theta}_{r}(\bm{\hat{r}},\bm{y};\bm{VaR}(\theta))$ loss, and the
most promising one is selected as the output for the entire step. Lastly, in
the third step, we repeat the previous routine with using the output from the
previous step as the only starting and the loss
$\mathcal{L}^{\theta}_{q,e}(\bm{\hat{e}},\bm{\hat{q}},\bm{y})$ instead of
$\mathcal{L}^{\theta}_{r}(\bm{\hat{r}},\bm{y};\bm{VaR}(\theta))$.
Three alternatives have been considered for the loss function used in the
training process. These include training only with the loss in Equation (13),
training with the loss in Equation (14), and training with both losses. The
loss function in Equation (13) involves only the optimization of the ES
coefficients, so it is much faster. On the other side, we expect the loss
function in Equation (14) to be more precise as it refines the VaR prediction
by incorporating information on the expected tail mean from the previous step.
So, as standard in the optimization literature, the faster method is initially
utilized to achieve a coarse approximation of the solution. This approximation
is then used as the starting point for the more precise method.
We empirically compare three different versions of CAESar. These include the
original version presented in Section 3, the version optimized based on the
loss in Equation (13) (B-CAESar), and the version derived solely from the loss
in Equation (14) (P-CAESar). The comparison is shown in Table 8.

$\bm{\theta=0.05}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0023 | 0.0022 | 0.0045 | 0.0025 | 0.0033 | 0.0032 | 0.0031 | 0.0025 | 0.0024 | 0.0019
$\mathcal{L}_{P}^{\theta}$ | 0.847 | 0.901 | 1.25 | 0.985 | 1.01 | 1.074 | 1.061 | 1.05 | 0.998 | 0.797
B-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.0023 | 0.006 | 0.0024 | 0.0035 | 0.0033 | 0.0032 | 0.0028 | 0.0025 | 0.0017
$\mathcal{L}_{P}^{\theta}$ | 0.834 | 0.908 | 1.297 | 0.983 | 1.053 | 1.041 | 1.095 | 1.069 | 1.007 | 0.798
P-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0022 | 0.0024 | 0.0045 | 0.0024 | 0.0034 | 0.0035 | 0.0032 | 0.0022 | 0.0024 | 0.0017
$\mathcal{L}_{P}^{\theta}$ | 0.844 | 0.902 | 1.247 | 0.983 | 1.027 | 1.034 | 1.102 | 1.018 | 1.003 | 0.796

$\bm{\theta=0.025}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0051 | 0.0047 | 0.0191 | 0.0058 | 0.0091 | 0.0075 | 0.0072 | 0.0048 | 0.0056 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.064 | 1.096 | 1.475 | 1.148 | 1.275 | 1.191 | 1.243 | 1.255 | 1.169 | 0.965
B-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0054 | 0.0057 | 0.0188 | 0.0057 | 0.0093 | 0.0072 | 0.0068 | 0.0105 | 0.0057 | 0.0033
$\mathcal{L}_{P}^{\theta}$ | 1.048 | 1.133 | 1.473 | 1.149 | 1.239 | 1.183 | 1.276 | 1.295 | 1.169 | 0.962
P-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0055 | 0.0044 | 0.0146 | 0.006 | 0.0081 | 0.0077 | 0.0072 | 0.005 | 0.0055 | 0.0032
$\mathcal{L}_{P}^{\theta}$ | 1.074 | 1.071 | 1.438 | 1.148 | 1.179 | 1.17 | 1.281 | 1.213 | 1.161 | 0.962

$\bm{\theta=0.01}$ |  | SPX | MXX | BOVESPA | CAC | KOSPI | IBEX | NIKKEI | HSI | DAX | FTSE
---|---|---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0224 | 0.0188 | 0.0451 | 0.0273 | 0.0353 | 0.026 | 0.0236 | 0.0212 | 0.0333 | 0.011
$\mathcal{L}_{P}^{\theta}$ | 1.427 | 1.376 | 1.652 | 1.394 | 1.493 | 1.389 | 1.482 | 1.476 | 1.461 | 1.191
B-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0285 | 0.0192 | 0.0439 | 0.0223 | 0.0363 | 0.0248 | 0.0214 | 0.0153 | 0.0332 | 0.0094
$\mathcal{L}_{P}^{\theta}$ | 1.787 | 1.455 | 1.743 | 1.374 | 1.498 | 1.378 | 1.493 | 1.438 | 1.461 | 1.166
P-CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0253 | 0.0105 | 0.0506 | 0.0275 | 0.0398 | 0.0225 | 0.023 | 0.0157 | 0.0252 | 0.0103
$\mathcal{L}_{P}^{\theta}$ | 1.867 | 1.261 | 1.698 | 1.397 | 1.591 | 1.407 | 1.473 | 1.45 | 1.401 | 1.184

Table 8: Comparison between different versions of CAESar. CAESar refers to the
version discussed in the main text, Section 3. B-CAESar is the version trained
only with the Barrera loss
$\mathcal{L}^{\theta}_{r}(\bm{\hat{r}},\bm{y};\bm{VaR}(\theta))$. P-CAESar is
the version that solely uses the Patton loss
$\mathcal{L}^{\theta}_{q,e}(\bm{\hat{e}},\bm{\hat{q}},\bm{y})$. The best
result is in bold, while the second to best is underlined.

This comparison shows a balanced situation, with all competing models closely
aligned. However, a deeper inspection of the results shows that the current
version appears to be the most stable, that is, less susceptible to negative
outliers. This is the rationale for our choice.

##  Appendix C Supplementary Material - Model Comparison Tests

This Supplementary Material section shows the results obtained by the Loss
Difference, the Encompassing, and the corrected resampled Student’s t-test.
Both are model comparison tests. The formers are applied asset by asset and
fold by fold, while the latter is applied once for each asset. These tests
look at the out-of-sample losses to asses whether the compared models are
statistically equivalent (null hypothesis). As for the Loss Difference test,
the exceptions’ number (p-value ¡ 0.05) is reported in Table 9.

|  | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---
$\bm{\theta=0.05}$ | $\mathcal{L}_{B}^{\theta}$ | 37 / 81 | 185 / 12 | 50 / 112 | 97 / 47 | 97 / 50
$\mathcal{L}_{P}^{\theta}$ | 44 / 26 | 228 / 0 | 132 / 5 | 88 / 3 | 127 / 8
$\bm{\theta=0.025}$ | $\mathcal{L}_{B}^{\theta}$ | 48 / 70 | 145 / 32 | 54 / 87 | 84 / 53 | 105 / 37
$\mathcal{L}_{P}^{\theta}$ | 40 / 40 | 165 / 2 | 113 / 5 | 115 / 9 | 126 / 5
$\bm{\theta=0.01}$ | $\mathcal{L}_{B}^{\theta}$ | 39 / 67 | 107 / 38 | 44 / 83 | 90 / 36 | 71 / 31
$\mathcal{L}_{P}^{\theta}$ | 50 / 31 | 106 / 9 | 99 / 9 | 117 / 13 | 109 / 8
Table 9: Model Comparison Tests - Loss Difference test with the p-value
threshold equal to 0.05. The results obtained across various assets, for a
given probability level $\theta$, have been aggregated by summing them. The
first value in each cell represents the number of folds where CAESar
outperforms the competitor. The second value is the number of folds where the
competitor outperforms CAESar.

From the table, we can recover a certain discrepancy between losses. Indeed,
although CAESar generally performs better than K-CAViaR and K-QRNN when
looking at $\mathcal{L}_{P}^{\theta}$, the opposite is true for the
$\mathcal{L}_{B}^{\theta}$ loss. The intuition behind this is the same as
before: as discussed in Section 5.1 of the main test, the Barrera loss
requires the correctness of the quantile estimate, which sometimes may be an
incorrect assumption. Finally, the results for the encompassing test are shown
in Table 10, which follows the same structure and content as the previous
tables.

|  | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---
$\bm{\theta=0.05}$ | $\mathcal{L}_{B}^{\theta}$ | 3 / 21 | 55 / 0 | 4 / 50 | 13 / 5 | 26 / 17
$\mathcal{L}_{P}^{\theta}$ | 29 / 18 | 127 / 0 | 29 / 18 | 40 / 9 | 26 / 4
$\bm{\theta=0.025}$ | $\mathcal{L}_{B}^{\theta}$ | 20 / 21 | 33 / 7 | 18 / 37 | 34 / 19 | 54 / 23
$\mathcal{L}_{P}^{\theta}$ | 23 / 24 | 50 / 5 | 41 / 20 | 51 / 24 | 38 / 20
$\bm{\theta=0.01}$ | $\mathcal{L}_{B}^{\theta}$ | 41 / 43 | 34 / 44 | 30 / 73 | 63 / 36 | 57 / 23
$\mathcal{L}_{P}^{\theta}$ | 46 / 41 | 46 / 27 | 57 / 32 | 73 / 35 | 57 / 23
Table 10: Model Comparison Tests - Encompassing test. The testing procedure
and the results presentation are the same as in the Loss Difference test.

The test results qualitatively are the same as those obtained by the Loss
Difference test. The only difference is in the rejection numbers. Indeed, the
Encompassing test is the most prudent test, while the Loss Difference test is
the most permissive.
Another perspective in comparing models is provided by the corrected resampled
Student’s t-test. We have applied this test once for each asset. This is
because it takes as input a vector where each component represents the mean
loss in a fold (as its population is made up of the losses among the folds).
Its results are shown in Table 11. The table shows the number of rejections of
the null hypothesis about models’ equivalence. Each cell represents the number
of indexes for which the row predictor overcomes the column one. Therefore,
the better models are those with low values in the corresponding column and
high values in their row. CAESar and K-CAViaR are confirmed to be the best-
performing models, with CAESar slightly outperforming. However, the test often
fails to reject the null hypothesis. This aligns with the findings in [13],
indicating that a testing procedure based on the joint eligibility of VaR and
ES requires very strong evidence for rejection.

$\bm{\theta=0.05}$ |  | CAESar | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---|---
CAESar |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | \- / - | 0 / 0 | 1 / 5 | 0 / 8 | 0 / 0 | 0 / 1
K-CAViaR |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | \- / - | 1 / 4 | 0 / 8 | 0 / 0 | 0 / 1
BCGNS |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | \- / - | 0 / 0 | 0 / 0 | 0 / 0
K-QRNN |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 4 | \- / - | 0 / 0 | 0 / 0
GAS1 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / - | 0 / 0
GAS2 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / -

$\bm{\theta=0.025}$ |  | CAESar | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---|---
CAESar |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | \- / - | 0 / 0 | 0 / 1 | 0 / 6 | 0 / 0 | 0 / 0
K-CAViaR |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | \- / - | 0 / 2 | 0 / 5 | 0 / 0 | 0 / 0
BCGNS |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | \- / - | 0 / 0 | 0 / 0 | 0 / 0
K-QRNN |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | \- / - | 0 / 0 | 0 / 0
GAS1 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / - | 0 / 0
GAS2 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / -

$\bm{\theta=0.01}$ |  | CAESar | K-CAViaR | BCGNS | K-QRNN | GAS1 | GAS2
---|---|---|---|---|---|---|---
CAESar |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | \- / - | 0 / 0 | 0 / 0 | 0 / 4 | 0 / 0 | 0 / 0
K-CAViaR |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | \- / - | 0 / 0 | 0 / 3 | 0 / 0 | 0 / 0
BCGNS |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | \- / - | 0 / 0 | 0 / 0 | 0 / 0
K-QRNN |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | \- / - | 0 / 0 | 0 / 0
GAS1 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / - | 0 / 0
GAS2 |  $\mathcal{L}_{B}^{\theta}$ / $\mathcal{L}_{P}^{\theta}$ | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | 0 / 0 | \- / -

Table 11: Model Comparison Tests - Corrected resampled t-test with confidence
0.05. Each cell value represents the number of assets where the row algorithm
performs statistically better than column one.

##  Appendix D Supplementary Material - Stock Dataset

To further validate our proposal, we carry out an additional experiment using
a stock dataset. Specifically, we select a dataset of stocks from the US
banking sector, as done in [29]. This new dataset consists of 8 stocks and
covers the period from 2000-01-01 to 2023-12-31. Again, we use block cross-
validation to evaluate performance, with 6 years for training/validation and 1
year for the test, resulting in 18 folds. The results from ES forecasting are
shown in Table 12. As before, both $\mathcal{L}_{B}^{\theta}$ and
$\mathcal{L}_{P}^{\theta}$ are shown.

$\bm{\theta=0.05}$ |  | BAC | BK | C | GS | JPM | MS | STT | WFC
---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0126 | 0.0138 | 0.0133 | 0.0087 | 0.0065 | 0.0262 | 0.0338 | 0.0088
$\mathcal{L}_{P}^{\theta}$ | 1.467 | 1.758 | 1.468 | 1.364 | 1.253 | 2.59 | 1.563 | 1.327
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.0128 | 0.0141 | 0.0134 | 0.0084 | 0.0067 | 0.0258 | 0.0347 | 0.0079
$\mathcal{L}_{P}^{\theta}$ | 1.419 | 1.627 | 1.483 | 1.408 | 1.26 | 3.345 | 1.547 | 1.337
BCGNS | $\mathcal{L}_{B}^{\theta}$ | 0.0532 | 0.0386 | 0.0948 | 0.02 | 0.041 | 0.0431 | 0.0716 | 0.0365
$\mathcal{L}_{P}^{\theta}$ | 3.424 | 3.201 | 3.33 | 2.298 | 2.497 | 6.646 | 3.0 | 4.818
K-QRNN | $\mathcal{L}_{B}^{\theta}$ | 0.0639 | 0.0257 | 0.075 | 0.0222 | 0.0265 | 0.0433 | 0.0565 | 0.0369
$\mathcal{L}_{P}^{\theta}$ | 2.153 | 1.797 | 2.211 | 1.74 | 1.906 | 1.987 | 1.984 | 2.146
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.0268 | 0.0111 | 0.0391 | 0.0143 | 0.0183 | 0.0201 | 0.0334 | 0.0074
$\mathcal{L}_{P}^{\theta}$ | 1.793 | 2.023 | 1.989 | 1.444 | 1.575 | 1.532 | 1.853 | 1.499
GAS2 | $\mathcal{L}_{B}^{\theta}$ | 0.5711 | 0.0255 | 0.1433 | 0.0238 | 0.0415 | 0.0236 | 0.2587 | 0.028
$\mathcal{L}_{P}^{\theta}$ | 2.0 | 1.97 | 2.123 | 3.412 | 1.871 | 2.489 | 2.213 | 2.123

$\bm{\theta=0.025}$ |  | BAC | BK | C | GS | JPM | MS | STT | WFC
---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.0425 | 0.0766 | 0.0292 | 0.0221 | 0.0159 | 0.0323 | 0.1072 | 0.0182
$\mathcal{L}_{P}^{\theta}$ | 2.352 | 2.034 | 1.679 | 1.597 | 1.476 | 1.714 | 1.962 | 1.543
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.0315 | 0.0773 | 0.0623 | 0.0236 | 0.0155 | 0.077 | 0.1075 | 0.0179
$\mathcal{L}_{P}^{\theta}$ | 1.639 | 1.759 | 1.763 | 1.587 | 1.476 | 2.229 | 1.769 | 1.544
BCGNS | $\mathcal{L}_{B}^{\theta}$ | 0.1456 | 0.0979 | 0.1475 | 0.0459 | 0.0609 | 0.1322 | 0.1908 | 0.0725
$\mathcal{L}_{P}^{\theta}$ | 2.587 | 3.35 | 2.849 | 2.757 | 2.051 | 2.832 | 3.311 | 3.084
K-QRNN | $\mathcal{L}_{B}^{\theta}$ | 0.2679 | 0.1076 | 0.3197 | 0.0798 | 0.1072 | 0.1852 | 0.2295 | 0.1709
$\mathcal{L}_{P}^{\theta}$ | 3.691 | 2.492 | 3.676 | 2.273 | 2.592 | 2.685 | 2.68 | 3.431
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.0413 | 0.0259 | 0.0402 | 0.0295 | 0.0116 | 0.117 | 0.1228 | 0.0243
$\mathcal{L}_{P}^{\theta}$ | 1.837 | 1.949 | 1.865 | 1.791 | 1.521 | 2.186 | 2.217 | 1.773
GAS2 | $\mathcal{L}_{B}^{\theta}$ | 0.0293 | 0.3891 | 0.2874 | 0.0367 | 0.057 | 0.0257 | 0.2466 | 0.0418
$\mathcal{L}_{P}^{\theta}$ | 3.188 | 2.456 | 2.829 | 2.516 | 3.164 | 2.168 | 2.068 | 2.255

$\bm{\theta=0.01}$ |  | BAC | BK | C | GS | JPM | MS | STT | WFC
---|---|---|---|---|---|---|---|---|---
CAESar | $\mathcal{L}_{B}^{\theta}$ | 0.2374 | 0.4018 | 0.1063 | 0.1102 | 0.0596 | 0.3967 | 0.4823 | 0.0352
$\mathcal{L}_{P}^{\theta}$ | 2.44 | 2.519 | 1.985 | 1.886 | 1.977 | 2.224 | 2.2 | 1.748
K-CAViaR | $\mathcal{L}_{B}^{\theta}$ | 0.4017 | 0.4149 | 0.1034 | 0.1473 | 0.0722 | 0.3981 | 0.6182 | 0.062
$\mathcal{L}_{P}^{\theta}$ | 2.828 | 3.104 | 1.969 | 1.951 | 1.82 | 2.211 | 2.23 | 1.817
BCGNS | $\mathcal{L}_{B}^{\theta}$ | 0.7912 | 0.3912 | 0.7121 | 0.1945 | 0.2417 | 0.5786 | 0.8724 | 0.2629
$\mathcal{L}_{P}^{\theta}$ | 5.339 | 3.434 | 3.387 | 2.471 | 2.318 | 2.877 | 3.309 | 2.565
K-QRNN | $\mathcal{L}_{B}^{\theta}$ | 0.564 | 0.4053 | 0.7435 | 0.2382 | 0.3498 | 0.6308 | 0.8121 | 0.3415
$\mathcal{L}_{P}^{\theta}$ | 2.54 | 2.828 | 2.695 | 2.267 | 2.898 | 2.504 | 2.77 | 2.537
GAS1 | $\mathcal{L}_{B}^{\theta}$ | 0.4972 | 0.3312 | 0.098 | 0.0469 | 0.3558 | 0.4366 | 0.8826 | 0.1189
$\mathcal{L}_{P}^{\theta}$ | 2.426 | 2.14 | 3.583 | 4.111 | 2.006 | 3.073 | 2.77 | 3.538
GAS2 | $\mathcal{L}_{B}^{\theta}$ | 0.4444 | 0.0755 | 0.3093 | 0.199 | 0.1359 | 0.0841 | 0.1363 | 0.2065
$\mathcal{L}_{P}^{\theta}$ | 2.557 | 2.241 | 2.545 | 3.019 | 2.649 | 2.342 | 2.675 | 2.71

Table 12: Out-of-sample loss for stock data. The best result is highlighted in
bold, the second to best is underlined. The $*$ is for those values that are
outside the CAESar mean + one standard deviation confidence interval.

Also in this case, the CAESar and K-CAViaR models outperform the competitors
(even if the gap between them is narrower). Also the GAS models perform well,
even if often they are overcome. Specifically, this rise in their performance
can be explained by taking in mind that the tails for indexes are usually
thinner than stocks ones [13]. Furthermore, as seen in Section 5, we recover
the inverse dependence of the loss function value by the confidence level
$\theta$. Finally, we report the direct approximation test result in Table 13.

| $\bm{\theta=0.05}$ | $\bm{\theta=0.025}$ | $\bm{\theta=0.01}$
---|---|---|---
| $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$ | $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$ | $\mathbf{MNF}$ | $\mathbf{AS1}$ | $\mathbf{AS2}$
CAESar | 0.13 | 0.12 | 0.17 | 0.1 | 0.9 | 0.14 | 0.32 | 0.26 | 0.2
K-CAViaR | 0.19 | 0.12 | 0.16 | 0.14 | 0.14 | 0.12 | 0.38 | 0.35 | 0.24
BCGNS | 0.28 | 0.13 | 0.12 | 0.41 | 0.25 | 0.25 | 0.61 | 0.39 | 0.4
K-QRNN | 0.43 | 0.43 | 0.53 | 0.48 | 0.48 | 0.37 | 0.62 | 0.58 | 0.41
GAS1 | 0.12 | 0.31 | 0.13 | 0.2 | 0.31 | 0.14 | 0.5 | 0.49 | 0.18
GAS2 | 0.34 | 0.49 | 0.42 | 0.49 | 0.52 | 0.4 | 0.69 | 0.6 | 0.46
Table 13: Direct approximation test for stock data. The table shows the
rejection of the null hypothesis in the McNeil and Frey test (MNF) and in the
tests based on the Acerbi and Szekely statistics: Z1 (AS1) and Z2 (AS2). As
the null hypothesis of all the tests is related to the correct specification
of the ES estimator, the lower the rejection rate, the better the model.
Finally, the best result is in bold, and the second to best is underlined.

When using the rejection rates to compare the algorithms, we find an
improvement in CAESar performances (with respect to the results in Table 12).
Indeed, while the competition with K-CAViaR is still close for $\theta=0.05$,
CAESar certainly defeats the competitors for lower confidence levels. Finally,
it is worth observing that BCGNS achieves a good result in the AS2 test with
$\theta=0.05$.