The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
56,461 characters
Dynamic tail risk forecasting: what do realized skewness and kurtosis add?
\begin{frontmatter}
\title{Dynamic tail risk forecasting: what do realized skewness and kurtosis add?}
\author[adrG1]{Giampiero M. Gallo\fnref{fn1}}
\ead{[email removed]}
\author[adrO1,adrO2]{Ostap Okhrin}
\ead{[email removed]}
\author[adrS1]{Giuseppe Storti\corref{cor1}\fnref{fn2}}
\ead{[email removed]}
\fntext[fn2]{Opinions expressed here are personal and do not involve the Corte dei conti.}
\cortext[cor1]{Corresponding author}
\address[adrG1]{Corte dei conti, New York University in Florence, and CRENoS}
\address[adrO1]{Technische Universität Dresden, 01062 Dresden, Germany}
\address[adrO2]{Center for Scalable Data Analytics and Artificial Intelligence (ScaDS.AI) Dresden/Leipzig, Germany}
\address[adrS1]{Universit\`a di Salerno, Department of Economics and Statistics, Fisciano, Italy}
\begin{abstract}
\noindent This paper compares the accuracy of tail risk forecasts with a focus on including realized skewness and kurtosis in "additive" and "multiplicative" models. Utilizing a panel of 960 US stocks, we conduct diagnostic tests, employ scoring functions, and implement rolling window forecasting to evaluate the performance of Value at Risk (VaR) and Expected Shortfall (ES) forecasts. Additionally, we examine the impact of the window length on forecast accuracy. We propose model specifications that incorporate realized skewness and kurtosis for enhanced precision. Our findings provide insights into the importance of considering skewness and kurtosis in tail risk modeling, contributing to the existing literature and offering practical implications for risk practitioners and researchers.
\vspace{.6cm}
\end{abstract}
\begin{keyword}
Value at Risk, CAViaR, Expected Shortfall, Realized Skewness, Realized Kurtosis.
\end{keyword}
\end{frontmatter}
\section{Introduction}
Starting approximately thirty years ago, the issue of capital adequacy has received increased attention, with significant impetus given to supervisory and regulatory functions to closely monitor the impact of volatility and interconnectedness on financial institution portfolios. Modern risk management is based on the principle that increased risks must be adequately covered with sufficient resources to avoid liquidity crises or defaults that could affect other institutions and the financial system as a whole. The consequences of the 2007-2008 financial crisis underscored the need for suitable capital risk measures exhibiting forecastability over relevant time horizons.
The various recommendations of the Basel Committee on Banking Supervision regarding capital risk regulations emphasize that the main parameters of a conditional distribution of returns to be monitored are some position index, specifically the threshold \citep[Value at Risk, $\mathsf{VaR}$, ][]{jorion1997value} corresponding to a certain probability in the tail where losses occur, and the average value of the loss once that threshold has been surpassed \citep[Expected Shortfall, $\mathsf{ES}$, ][]{artzner1999}. In this context, without loss of generality, we assume that the tail in question is the left tail, representing losses in long positions.
Market activity, characterized by price and volume movements in response to news, renders the conditional distribution of returns non-constant over time. Consequently, both Value at Risk ($\mathsf{VaR}$) and Expected Shortfall ($\mathsf{ES}$) become time-varying risk measures. Moreover, observed persistence in market behavior suggests dynamics that leverage valuable past information. From an econometric perspective, it is challenging to determine which features of past market behavior are relevant for predicting $\mathsf{VaR}$\ and $\mathsf{ES}$, as these measures represent conditional quantiles and expectations, respectively, in the tail of the asset return distribution.
Approaches to address this issue can broadly be categorized into three main groups. The first category assumes a known parametric distribution for returns, typically a Student-$t$ distribution, and focuses on the dynamic evolution of the conditional variance of returns. This approach augments the fixed quantile identification with a GARCH process that models the dependence of conditional variance on recent returns and past estimates. Parameters are estimated using (Quasi) Maximum Likelihood (QML) methods.
At the opposite end of the spectrum, parametric assumptions about the return distribution or its dynamics are entirely discarded. So-called historical simulation methods are employed, where future outcomes are simulated by repeating observed past behaviors.
A third stream adopts an intermediate stance, focusing on the dynamics of the risk measure of interest while limiting or avoiding reliance on parametric assumptions about the shape of the conditional return distribution. This semi-parametric approach to financial risk modeling is gaining popularity due to its flexibility and often demonstrates competitive performance compared to more complex parametric models. In what follows, we will position ourselves in this stream of literature, addressing, in particular, the role that higher-order conditional moments, notably skewness and kurtosis, have on the refinement of predictions, hence highlighting the role of the time-varying evolution of asymmetry and tail density of the return distribution in sharpening the projections of $\mathsf{VaR}$\ and $\mathsf{ES}$.
Our synthesis in this field is to identify two main categories of semi-parametric modeling approaches for tail risk measures. The first is called the ``additive'' approach, which utilizes linearized representations of GARCH models, such as in CAViaR models. The second approach, referred to as the ``multiplicative'' approach, involves estimating GARCH-type models via the minimization of a properly defined strictly consistent scoring function. We consider the recent literature \citep[e.g.,][]{neuberger2013,neuberger2020,baelee2020} on the derivation of realized measures of skewness and kurtosis as consistent estimates of the conditional skewness and kurtosis of daily returns. For our purposes, these additional features of the conditional distributions may be relevant when included in the specifications. Given that \cite{amaya2015} provides evidence that realized skewness and kurtosis are useful when forecasting the cross-section distribution of equity returns, our interest here is to assess whether these benefits extend to risk forecasting as well.
Although the additive approach has gained popularity, there is still a lack of extensive forecasting comparison between these two methodologies. Hence, we aim to bridge a gap in the literature by proposing an application that evaluates the accuracy of forecasts generated by additive and multiplicative modeling strategies for a panel of 960 US stocks. To achieve this, we employ various diagnostic tests and scoring functions for both $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasts. Additionally, we investigate the impact of window length on forecasting accuracy, a critical issue for practitioners. Short windows tend to minimize bias but increase variability in risk forecasts, while long windows have the opposite effect: hence, we conduct a rolling window estimation/forecasting exercise and evaluate the performance of three different window lengths, 500, 1000, and 2000 days.
Our novel model specifications using information on realized higher-order moments to forecast tail risk measures are both regression quantile time series models for forecasting $\mathsf{VaR}$, as well as bivariate semi-parametric models for joint $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasting: we are interested in providing specific evidence on the relevance of the realized skewness and kurtosis via Wald-type tests, but also on their contribution in improving the forecast performance, assessed with standard backtesting procedures. Their predictive performances are compared to those of competitors that do not include such information.
In a nutshell, the evidence on the vast panel of stock indicates that multiplicative models are preferred to additive ones, and that the extension to higher moments does not buy a generalized relevant improvement in the outcome. In general, simpler models are to be preferred to more complex ones.
The structure of the paper is as follows. In Section \ref{s:mods}, we propose our models in Subsection \ref{s:setup}, while the related estimation procedures and some properties of the estimators are illustrated in Subsection \ref{s:esti}. In Section \ref{s:revsk}, we present a recent literature review on realized estimators of skewness and kurtosis of financial returns. Section \ref{s:appli} is dedicated to the empirical application, while Section \ref{s:conc} concludes.
\section{Semi-parametric risk modeling}
\label{s:mods}
The literature on semi-parametric risk modeling features a seminal paper by \cite{caviar}, who introduced the Conditional Autoregressive Value-at-Risk (CAViaR) model for forecasting $\mathsf{VaR}$. This model has interesting connections with both quantile regression and GARCH models, in that the CAViaR model can be viewed as a quantile autoregression with a recursive term. By the same token, a linear GARCH model of a given order can be represented as a CAViaR model of the same order. Building on the duality between GARCH and CAViaR, \cite{Xiao_Koenker_2009} present an original approach to estimating parameters of a GARCH model, proposing to minimize the typical quantile loss function used in quantile regression models.
Direct semi-parametric modeling of $\mathsf{ES}$\ is not feasible because, unlike $\mathsf{VaR}$, $\mathsf{ES}$\ is not elicitable relative to a given loss function. However, \cite{fissleretal2015} have derived a class of loss functions that are strictly consistent for the pair ($\mathsf{VaR}$, $\mathsf{ES}$), in the sense that the expected loss is minimized by the true ($\mathsf{VaR}$, $\mathsf{ES}$). Within this framework \cite{tayl2019} proposes a class of semi-parametric models for ($\mathsf{VaR}$, $\mathsf{ES}$), augmenting the standard CAViaR setup with an additional dynamic equation for $\mathsf{ES}$, and replacing the usual quantile loss with a member of the Fissler-Ziegel (FZ) class. In particular, among the available choices, \cite{taylor2017} considers a loss, or scoring, function based on the Asymmetric Laplace quasi-likelihood function, AL for short. \cite{pattonetal2019} extend the work by \cite{tayl2019} in two different directions. First, they consider time-varying semi-parametric ($\mathsf{VaR}$, $\mathsf{ES}$) models based on the Generalized Autoregressive Score (GAS) framework \citep{crealetal2013}. Second, as done by \cite{Xiao_Koenker_2009} for $\mathsf{VaR}$, they consider directly estimating GARCH models minimizing a specific strictly consistent loss function in the FZ class called FZ0 (owing its denomination to the fact that, when using this loss to compare two models, it yields loss differentials that are homogeneous of degree zero). This property can lead to a higher power in Diebold-Mariano tests \citep{diebold1995comparing}.
\subsection{The model setup}
\label{s:setup}
Let $r_t$ be the log-return for the day $t$, for $t=1,\ldots, T$, and $Q_{\alpha,t} = F_{r}^{-1}(\alpha|\mathcal{I}_{t-1})$ indicate the conditional $\alpha$-quantile of $r_t$ (level-$\alpha$ Value-at-Risk --$\mathsf{VaR}$), with $F_r$ being the cdf of $r_t$; correspondingly, $ES_{\alpha,t}= \operatorname{E}(r_{t}|r_{t}<Q_{\alpha,t},\mathcal{I}_{t-1})$ indicates the conditional $\alpha$-tail expectation of $r_t$, given past information $\mathcal{I}_{t-1}$ (level-$\alpha$ Expected Shortfall -- $\mathsf{ES}$).
We let $\operatorname{RV}_t$, $\operatorname{Sk}_t$, and $\operatorname{Ku}_t$ denote, respectively, the conditional variance, skewness, and kurtosis of daily returns $r_t$ as follows
\begin{align*}
\operatorname{RV}_{t} &= \operatorname{E}_0\{(r_t-\mu_{1t})^2|\mathcal{I}_{t-1}\} = \mu_{2t} - \mu_{1t}^2 ,\\
\operatorname{Sk}_{t} &= \operatorname{E}_0 \left\{\left(\frac{r_t-\mu_{1t}}{\operatorname{RV}^{1/2}_t}\right)^3| \mathcal{I}_{t-1}\right\}= \frac{\mu_{3t} - \mu_{1t}\mu_{2t}+2\mu_{1t}^3}{(\mu_{2t} - \mu_{1t}^2)^{3/2}} ,\\
\operatorname{Ku}_{t} &= \operatorname{E}_0 \left\{\left(\frac{r_t-\mu_{1t}}{\operatorname{RV}^{1/2}_t}\right)^4| \mathcal{I}_{t-1}\right\}=\frac{\mu_{4t} - 4\mu_{3t}\mu_{1t}+6\mu_{2t}\mu_{1t}^2-3\mu_{1t}^4}{(\mu_{2t} - \mu_{1t}^2)^{2}},
\end{align*}
\noindent where $\mu_{kt} = \operatorname{E}_0(r_t^k|\mathcal{I}_{t-1})$ indicates the $k$-th conditional noncentered moment under the true measure. Estimates of these quantities can be readily obtained by replacing the involved conditional moments $\mu_{kt}$ with their estimated counterparts, at least using daily observations. In Section \ref{s:revsk}, we will formally address the estimation of $\mu_{kt}$ for $1\leq k \leq 4$.
We can now present the two alternative modeling frameworks under which $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasts are generated, denoted as, for ease of reference, the \emph{additive} and the \emph{multiplicative} models, respectively. We simplify the notation by defining $v_{t}\equiv Q_{\alpha,t}$ and $e_{t}\equiv ES_{\alpha,t}$. Thus, the \emph{additive modeling framework} can be represented as a regression model for the 1-step ahead expected $\alpha$-level of $\mathsf{VaR}$, $v_{t}$:
\begin{align*}
r_t&= v_{t}+\eta_t,
\end{align*}
where the error term $\eta_t$ is controlling the left tail of the conditional distribution of returns, so that, under correct specification of $v_{t}$, the error term $\eta_t$ is such that $F_{\eta}^{-1}(\alpha|\mathcal{I}_{t-1}) = 0$.
This is a general framework since several models can be derived as special cases by varying the dynamic specifications for $v_{t}$. Noting that a \, $\widehat \cdot$ \, is used to indicate an estimate, $\bar r = T^{-1}\sum_{t = 1}^T r_t$ and $\widehat\operatorname{Sk}^{+}_t$ and $\widehat\operatorname{Sk}^{-}_t$, represent negative and positive skewness:
\begin{equation*}
\widehat\operatorname{Sk}_{1t}^+ = \widehat\operatorname{Sk}_{1t} \cdot I\{\widehat\operatorname{Sk}_{1t} > 0\},\qquad\qquad
\widehat\operatorname{Sk}_{1t}^- = -\widehat\operatorname{Sk}_{1t} \cdot I\{\widehat\operatorname{Sk}_{1t} < 0\},
\end{equation*}
in what follows, we investigate three specifications.
The first is the simple additive form of the VaR being driven only by the lagged observation of an estimator of the integrated volatility ($\widehat\operatorname{RV}_{t-1}$)\footnote{In the absence of jumps, the integrated variance coincides with the quadratic variation that, in turn, diverges from the conditional variance $\operatorname{RV}_t$ by a zero mean error, thus motivating our notation \citep{ande01_jasa}. A set of alternative choices for $\widehat\operatorname{RV}_{t-1}$ will be presented and discussed in Section \ref{s:revsk}.}.
\begin{equation}
\mbox{\texttt{add_sim}:}\ \ v_{t}= d_0 + d_1 \widehat\operatorname{RV}^{1/2}_{t-1} + d_2 v_{t-1},\label{eq:add_sim}
\end{equation}
To account for the potential misspecification in \texttt{add_sim}, we can resort to the Cornish-Fisher (CF) expansion \citep{HilDav1968}, which approximates the quantiles of an unknown non-Gaussian distribution using the information on sample skewness and kurtosis to adjust the value of the corresponding Gaussian quantiles. Considering as an illustration a random variable $X\sim(0,1)$, the CF approximation for the $\alpha$-quantile of $X$ reads as
\begin{equation*}
X^{CF}_{\alpha}=z_\alpha+\frac{z^2_\alpha-1}{6}Sk +\frac{z^3_{\alpha}-3 z_{\alpha}}{2}Ku-\frac{2z^3_{\alpha}-5 z_{\alpha}}{36}Sk^2,
\end{equation*}
where $Sk$ and $Ku$ are the usual moment-based sample skewness and kurtosis coefficients of $X$ respectively, and $z_\alpha=\Phi^{-1}(\alpha)$ is the $\alpha$-quantile of a $N(0,1)$ random variable.
Therefore, the second model adds realized negative and positive skewnesses ($\operatorname{Sk}^-_{t-1}$ and $\operatorname{Sk}^+_{t-1}$) and kurtosis ($\operatorname{Ku}_{t-1}$) to the \texttt{add_sim}\footnote{In our approach, we focus on the conditional distribution of returns rather than on their unconditional distribution, as would happen when using the standard CF expansion. Hence, the sample skewness and kurtosis coefficients are replaced by their realized counterparts that provide point estimates of daily conditional skewness and kurtosis.}:
\begin{equation}
\mbox{\texttt{add_skk}:}\ \ v_{t} = d_0 + d_1 \widehat\operatorname{RV}^{1/2}_{t-1} + d_2 v_{t-1}+(a_1 \widehat\operatorname{Sk}^-_{t-1}+a_2 \widehat\operatorname{Sk}^+_{t-1}+a_3 \widehat\operatorname{Ku}_{t-1}).\label{eq:add_skk}
\end{equation}
The inclusion of the skewness and kurtosis terms are thus motivated by a data-driven CF expansion, whose coefficients, as it will be later illustrated, can be estimated in a semi-parametric fashion by minimizing a strictly consistent loss function.
The third model further extends the \texttt{add_skk} with an asymmetric impact of the integrated volatility in correspondence with returns smaller than their average (leverage effect):
\begin{equation}
\mbox{\texttt{add_lev}:}\ \ v_{t} = d_0 + d_1 \widehat\operatorname{RV}^{1/2}_{t-1} + d_2 v_{t-1}+(a_1 \widehat\operatorname{Sk}^-_{t-1}+a_2 \widehat\operatorname{Sk}^+_{t-1}+a_3 \widehat\operatorname{Ku}_{t-1}) + d_3 \widehat\operatorname{RV}^{1/2}_{t-1}I\{r_{t-1}\leq \bar r\}.\label{eq:add_lev}
\end{equation}
By contrast, a \emph{multiplicative modeling framework} can be represented in terms of the following nonlinear regression model
\begin{align*}
r_t&= v_{t}\,\eta_t,
\end{align*}
where, under correct specification of $v_{t}$, $\eta_t$ is such that $F_\eta^{-1}(\alpha|\mathcal{I}_{t-1})=1$.
This framework can also be motivated by a simple location-scale representation of the returns process
\begin{equation*}
r_t=h_t z_t, \qquad\qquad z_t \overset{iid}{\sim} (0,1),
\end{equation*}
where the dynamics of $h_t^2=\mathop{\mbox{\sf Var}}(r_t|\mathcal{I}_{t-1})$ can be modelled by means of GARCH type models. Under the iid assumption for $z_t$ the 1-step ahead $\alpha$-level $\mathsf{VaR}$\ of $r_t$ is given by $v_{t} = h_t z_{\alpha}$ where $z_{\alpha}=F_z^{-1}(\alpha)$. Thus, the first multiplicative model is for the $\alpha$-level 1-step ahead $\mathsf{VaR}$ is
\begin{equation}
\mbox{\texttt{mlt_sim}:}\ v_{t}=h_{t}, \quad
h^2_{t}=d_0+d_1\widehat\operatorname{RV}_{t-1}+d_2h^2_{t|t-1}.\label{eq:mlt_sim}
\end{equation}
When we allow for a time-varying conditional skewness and kurtosis in the returns distribution, this assumption must be generalized to read
\begin{equation*}
v_{t}=h_t z_{\alpha,t},
\end{equation*}
where the time variation in the conditional error quantile $z_{\alpha,t}$ is driven by the time-varying conditional skewness and kurtosis values as in the \texttt{mlt_lev} and \texttt{mlt_skk} specifications introduced below. Thus, within the multiplicative modeling framework, we consider the following alternative specifications for the $\alpha$-level 1-step ahead $\mathsf{VaR}$
\begin{align}
\mbox{\texttt{mlt_skk}:}\ \ & v_{t} = h_{t}(a_1 \operatorname{Sk}^-_{t-1}+a_2 \operatorname{Sk}^+_{t-1}+a_3 \operatorname{Ku}_{t-1}), & &
h^2_{t} = d_0+d_1\widehat\operatorname{RV}_{t-1}+d_2h^2_{t|t-1},\label{eq:mlt_skk}\\
\mbox{\texttt{mlt_lev}:}\ \ &v_{t} = h_{t}(a_1 \operatorname{Sk}^-_{t-1}+a_2 \operatorname{Sk}^+_{t-1}+a_3 \operatorname{Ku}_{t-1}), & & h^2_{t}=d_0+d_1\widehat\operatorname{RV}_{t-1}+d_2h^2_{t|t-1}+d_3 \widehat\operatorname{RV}_{t-1}I\{r_{t-1}\leq \bar r\}\label{eq:mlt_lev}.
\end{align}
Multiplicative models closely mirror additive models (\ref{eq:add_sim}), (\ref{eq:add_skk}) and (\ref{eq:add_lev}). The \texttt{mlt_sim} is similar to (\ref{eq:add_sim}) and is the simplest specification with only the integrated volatility driving the dynamics of the scale. Further \texttt{mlt_skk} assumes similar to (\ref{eq:add_skk}) in the additional information incorporated in realized skewness and kurtosis that drives the dynamics of the scale. The most complex model \texttt{mlt_lev} also controls for the leverage in $h_t$ in the same fashion as in the model (\ref{eq:add_lev}).
For both additive and multiplicative frameworks, the $\mathsf{ES}$\ can be modeled according to two different alternative specifications:
\begin{eqnarray}
\mbox{\texttt{ES_sim}:}\quade_{t}&=&\{1+\exp(b_0)\}v_{t}\label{eq:ES_sim},\\
\mbox{\texttt{ES_skk}:}\quade_{t}&=&\{1+\exp{\left(b_0 + b_1 \operatorname{Sk}_{t-1}+b_2 \operatorname{Ku}_{t-1}\right)}\}v_{t}.\label{eq:ES_skk}
\end{eqnarray}
Here \texttt{ES_sim} is the simple specification assuming that the $\mathsf{ES}$\ is a rescaling of $\mathsf{VaR}$. \cite{tayl2019} shows that this simple specification provides competitive $\mathsf{VaR}$\ forecasts. More recently, \cite{wang_gerl_chen_2023} have extended the framework proposed in \cite{tayl2019} to allow for separate $\mathsf{VaR}$\ and $\mathsf{ES}$\ dynamics as well as for the incorporation of realized measures.
The more complex \texttt{ES_skk} brings the dynamics of the $\mathsf{ES}$\ to be also driven by the skewness and kurtosis, possibly accounting for the misspecification of \texttt{ES_sim}. Differently from $\mathsf{VaR}$, in this case, we did not split the skewness into negative and positive. By construction, both \texttt{ES_skk} and \texttt{ES_sim} specifications avoid the crossing of $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasts.
In the additive case, neglecting to model the $\mathsf{ES}$\ dynamics leads to a pure $\mathsf{VaR}$\ model. This case, labeled as \texttt{ES_no}, corresponds to a model specification that is close in spirit to a \cite{caviar} CAViaR type one where some realized estimator of the integrated variance replaces the volatility measure based on lagged daily returns.
\subsection{Estimation}
\label{s:esti}
Estimation of the vector $\boldsymbol{\theta}$ of unknown parameters describing the models for $v_t$ and $e_t$ both in the additive and multiplicative models (\ref{eq:add_sim})-(\ref{eq:ES_skk}) is done semi-parametrically by minimizing a strictly consistent scoring rule,
\begin{equation}
\hat{\boldsymbol{\theta}}=\underset{\boldsymbol{\theta}}{\arg\min}
\sum_{t=1}^{T}S^{(\alpha)}_t,\label{e:optim}
\end{equation}
where $S^{(\alpha)}_t$ is a member of the general class presented by \cite{Fissler2016}, i.e.
\begin{align*}
S^{(\alpha)}_t \equiv S(v_{t},e_{t}|r_t;\alpha)&=\{I(r_t\leq v_{t})-\alpha\}G_1(v_{t})-I(r_t\leqv_{t})G_1(r_t)+G_2(e_{t})\\
&\left\{e_{t}-v_{t}+I(r_t\leqv_{t})\frac{v_{t}-r_t}{\alpha}\right\}
-\zeta_2(e_{t})+a(r_t).\label{e:scorule}
\end{align*}
In the definition of $S^{(\alpha)}_t$, the functions $G_1$, $\zeta_2$, and $G_2$ satisfy the following conditions: $G_1$ is increasing, $\zeta_2$ is increasing and convex, and $G_2=\zeta^'_2$. In particular, setting $G_1(\cdot)=0$, $G_2(x)=-1/x$, $\zeta_2(x)=-\log(-x)$, $a(r_t)=1-\log(1-\alpha)$, leads to the following scoring rule \citep{tayl2019}
\begin{align*}
AL^{(\alpha)}_t&=\frac{I(r_t\leqv_{t})r_t+v_{t}\{\alpha-I(r_t\leqv_{t})\}}{\alpha e_{t}}+\log(-e_{t})-\log(1-\alpha)\\
&=\frac{I(r_t\leqv_{t})r_t+v_{t}\{\alpha-I(r_t\leqv_{t})\}}{\alpha e_{t}}-\log\left(\frac{1-\alpha}{e_{t}}\right).
\end{align*}
Adding and subtracting $\alpha r_t$ to the numerator of the second term on the right-hand side of the previous equation, we get
\begin{equation*}
AL^{(\alpha)}_t=-\log\left(\frac{\alpha-1}{e_{t}}\right)-\frac{(r_t-v_{t})\{\alpha-I(r_t\leqv_{t})\}}{\alpha e_{t} } + \frac{r_t}{e_{t}}.
\end{equation*}
In the simplified expression of $AL^{(\alpha)}_t$ obtained by \cite{tayl2019}, the last term on the right-hand side is dropped. This simplification arises from the assumption that the conditional mean of returns is zero, as shown in their equation (19). Notably, it can be demonstrated that the negative value of $AL^{(\alpha)}_t$ quantifies the contribution of the $t$-observation to a quasi-likelihood function, which is constructed based on the Asymmetric Laplace distribution \citep{tayl2019}.
With the same choices for $G_1$, $G_2$, and $\zeta_2$ as above, but setting $a(r_t)=0$, leads to the zero-degree homogeneous loss used by \cite{pattonetal2019}
\begin{equation*}
FZ0^{(\alpha)}_t=\frac{I(r_t\leqv_{t})r_t+v_{t}\{\alpha-I(r_t\leqv_{t})\}}{\alpha e_{t}}+\log(-e_{t})-1=
-\frac{I(r_t\leqv_{t})(v_{t}-r_t)}{\alpha e_{t}}+\frac{v_{t}}{e_{t}}+ \log(-e_{t})-1.
\end{equation*}
Therefore, in the case of the model \texttt{EM_no}, where only the $\mathsf{VaR}$\ is estimated, the objective function is given by the quantile loss:
\begin{equation}\label{e:emloss}
EM^{(\alpha)}_t = \{\alpha - I(r_t < v_{t})\} \cdot (r_t - v_{t}).
\end{equation}
In what follows, models estimated relying on loss functions $AL^{(\alpha)}_t$, $FZ0^{(\alpha)}_t$ and $EM^{(\alpha)}_t$ are labeled as \texttt{Loss=ALS}, \texttt{Loss=FZ0} and \texttt{Loss=EM}, respectively. Standard errors are computed using the asymptotic theory developed by \cite{caviar}, for pure $\mathsf{VaR}$\ models, and \cite{pattonetal2019}, for joint $\mathsf{VaR}$-$\mathsf{ES}$\ models. Technical details are provided in Section \ref{s:appA}.
It is worth noting that optimization in the (\ref{e:optim}) is a challenging task irrespective of what function $S_t^{(\alpha)}$ is chosen, be it either $AL^{(\alpha)}_t$, $FZ0^{(\alpha)}_t$ or $EM^{(\alpha)}_t$. In particular, the optimization of these loss functions is typically strongly dependent upon the chosen initial values. For this reason, we implemented an optimization technique similar to \cite{caviar} which, for ease of reference, we call \emph{complete estimation}. Namely, for each model, we evaluated the objective function on $\mathfrak{n} = 5 \cdot 10^4$ uniformly sampled possible parameter constellations, and among them, we selected the $\mathfrak{m} = 10$ parameter vectors that lead to the smallest objective function values. Selecting each of these $\mathfrak{m}$ vectors as a starting point, we re-estimated the model $\mathfrak{m}$-times iterating between a Nelder-Mead and a BFGS optimizer until convergence is achieved, and the final estimates are those delivering the smallest value of the objective function. In a rolling window forecasting exercise, one may be advised to follow a parsimonious estimation strategy, by using the most recent estimates as the starting point for the next estimation round, at regular intervals.
\section{The underlying process and the derived realized measures}
\label{s:revsk}
Having developed a setup where the theoretical estimators of conditional moments such as $\operatorname{RV}$, $\operatorname{Sk}$, and $\operatorname{Ku}$ are considered within suitable models, we are left with the delicate phase to choose which operational counterparts we can count on at daily frequencies, employing rolling windows and sample statistics. The standard framework starts from a true underlying continuous log-price following a diffusion process, disregarding, for example, the presence of structural breaks:
\begin{equation*}
\mathrm{d} X_\mathfrak{t} = \mu(X_\mathfrak{t})\mathrm{d}\mathfrak{t} +\sigma(X_\mathfrak{t})\mathrm{d}W_\mathfrak{t},
\end{equation*}
where, $W_\mathfrak{t}$ represents the standard Brownian motion, $\mu(X_\mathfrak{t})$ is the drift càdlàg finite variation process, and $\sigma(X_\mathfrak{t})$ is the time-varying càdlàg volatility function. It is important to note that $\sigma(X_\mathfrak{t})$ may depend on a separate Brownian motion, which could potentially be correlated with $W_\mathfrak{t}$. This general family encompasses well-known processes such as the Heston or Bates processes (see \cite{heston1993closed, bates1996jumps}). In this context, the parameter $\mathfrak{t}$ represents the continuous temporal component that spans within and across days.
The second moment $\mu_{2t}$ is known as the \emph{integrated variance}, an object of paramount importance to researchers and practitioners. By utilizing the aforementioned process over a one-day interval $[t-1d, t]$, the integrated variance can be computed as $\int_{t-1d}^t \sigma^2(u)\mathrm{d}u$.
The temporal component then needs to be somehow aggregated to get the daily estimators for the relevant moments: in this respect, we ground ourselves in the massive literature that considers the market activity of a day (using the same index $t \in \{1, \ldots, T\}$) between opening and closing to be divided into regularly spaced intervals $i \in \{0, \ldots, N\}$. We then take the high-frequency log-prices $x_{t,i}$ as the elementary information, to be converted into $r_{t, i} = x_{t,i} - x_{t, i-1}$, the corresponding \emph{intraday} log-returns, $i=1,\ldots,N$.
The overwhelming attention of the literature was devoted to the design of consistent estimators of the integrated variance $\mu_{2t}$ of the continuous process over a discrete interval \citep{andersen2010parametric}, with specific care devoted to departures from the standard framework (e.g. jumps) or to the nature of observed prices which are affected by trading mechanisms (so-called market microstructure). There exists a range of options for researchers and practitioners alike seeking to estimate these quantities accurately and efficiently. Starting from the realized variance \citep{AndersenBoll98},
\begin{equation*}
\widehat\mu_{2t}^{RV} = \sum_{i = 1}^{N} r_{t,i}^2,
\end{equation*}
other widely used estimators of the integrated volatility are the, proposed in \cite{BarndorffNielsen_Bipower_2004} and \cite{andersen2012jump}, bipower variation $\widehat\mu_{2t}^{BPV} = \frac{\pi}{2}\frac{N}{N-1}\sum_{i=1}^{N-1}|r_{t,i}||r_{t,i+1}|$, or the upside and downside semivariances $\widehat\mu_{2t}^{SVPOS}= \sum_{i=1}^{N} r_{t,i}^2 \ \cdot \ I \{ r_{t,i} > 0\}$ and $\widehat\mu_{2t}^{SVNEG} = \sum_{i=1}^{N} r_{t,i}^2 \ \cdot \ I \{ r_{t,i} < 0\}$ developed in \cite{Barndorff_Semivariance_2008} and \cite{Bollerslev_SemiCovar_2020}.
Barring a horse race among the many estimators of $\mu_{2t}$ available, we limit ourselves to a single choice, and our preference goes to the median estimator
\begin{equation*}
\widehat\mu_{2t}^{MED}=\frac{\pi}{6-4\sqrt{3}+\pi}\frac{N}{N-2}\sum_{i=2}^{N-1} \mbox{med}(|r_{t,i-1}|,|r_{t,i}|, |r_{t,i+1}|)^2,
\end{equation*}
proposed by \cite{andersen2012jump}, because of its documented robustness properties.
Several new estimators for \emph{higher-order} moments have emerged in recent years. While these estimators do not directly estimate daily skewness or kurtosis, they instead estimate the integrated third or fourth power of intraday returns or the averaged jump component. Empirical evidence suggests that these estimators can be informative in predicting cross-sectional next week's stock returns or in forecasting RV at medium- to long-term horizons, as demonstrated by \cite{MeiLiuMaChen2017}. The simplest estimator of the integrated $k$-order moments was proposed by \cite{amaya2015}, shadowing the relationship between the realized variance and the integrated volatility (case $k = 2$).
\begin{equation*}
\widehat{\mu}_{kt}^{ACJV} = \sum_{i = 1}^{N} r_{t,i}^k.
\end{equation*}
\cite{amaya2015} demonstrate the estimator's consistency, which asymptotically captures only the jump component and the average jump size but does not capture skewness arising from the leverage effect and heavily depends on the sampling frequency. Later \cite{LiuWangLiu2014} derived asymptotic properties of the \cite{amaya2015} estimator and developed their own measures of realized skewness accounting for market microstructure noise. Another extension was provided by \cite{ChoeLee2014}, who showed that the daily third moment is proportional to the quadratic covariation between the squared return and the return process, and the fourth moment is proportional to the quadratic variation of the squared return process with some additional cross‐terms.
Based on some preliminary analysis, our choice for the realized skewness and kurtosis falls on the estimators by
\cite{neuberger2013} and \cite{neuberger2020}:
\begin{align*}
\widehat{\mu}_{3t}^{NP} & = \frac{1}{\tau} \sum_{j = 0}^{\tau-1} \sum_{i = 1}^{N} \left(r_{t-j,i}^3 + 3y_{t-j,i-1}^* r_{t,i}^2\right),
&\widehat{\mu}_{4t}^{NP} & = \frac{1}{\tau} \sum_{j = 0}^{\tau-1} \sum_{i = 1}^{N} \left(r_{t-j,i}^4 + 4 y^*_{t-j,i-1} r_{t-j,i}^3 + 6z_{t-j,i-1}^*r_{t-j,i}^2\right),
\end{align*}
where $y_{t,i-1}^* = \frac{1}{N}\sum_{j=1}^{N}(x_{t,i-1} - x_{t,i-j})$ and $z^*_{t,i-1} = \frac{1}{N}\sum_{j=1}^{N}(x_{t,i-1} - x_{t,i-j})^2$ measure local (daily) trends in simple and squared log-prices. Similar to \cite{ChoeLee2014}, they assume that the conditional mean of the returns is zero. Also, in what follows, we choose $\tau = 5$.
We note that the estimators for realized skewness and kurtosis sometimes produce outliers that can significantly affect the performance of $\mathsf{VaR}$\ and $\mathsf{ES}$\ models. To address this issue, we applied a filter that removes estimated skewness and kurtosis values falling outside the ranges of $(-15;15)$ and $(0;20)$, respectively. Outliers excluded from our analysis are then smoothed out using interpolation techniques accounting for autocorrelation.
\section{Empirical evidence}
\label{s:appli}
\subsection{Data and forecasting design}
In this section, we present the results of our setup to a very large panel of 960 U.S. stocks traded on the New York Stock Exchange (NYSE), included in the S\&P500 index at various times over the considered period. The list of stocks can be found in Web Appendix \textit{List of Tickers}. The original dataset for each stock consists of intra-daily prices adjusted for stock splits and dividends sampled every 5 minutes. We focus only on regular trading hours, from 9:30 am to 4:00 pm, resulting in 78 observations for each trading day. The stocks have different timespans, starting within a range between 1998-01-02 and 2016-10-11, and ending between 1998-01-09 and 2017-02-09. Furthermore, to ensure an adequate sample size, we limit our analysis to assets with a continuous record of at least 500 daily observations. This reduces the cross-sectional size of our sample to 823 assets (marked in the Web Appendix \textit{List of Tickers} in italics).
Our empirical strategy consists of two main steps. In the first, we conduct a full-sample analysis to assess the performance of various models in fitting $\mathsf{VaR}$\ and $\mathsf{ES}$. In the second step, we focus on the out-of-sample forecasting performance using a rolling window approach. We consider three estimation windows: 500, 1000, and 2000 days. For the out-of-sample analysis, we include assets with a continuous record of daily pricing observations from the start date of our sample to its end, 2017-02-09. In this case, we were able to obtain one-step ahead predictions for the dates 2000-01-04 -- 2017-02-09, for $w = 500$, 2002-01-03 -- 2017-02-09 for $w = 1000$, and 2005-12-21 -- 2017-02-09 for $w = 2000$. 406 of the original 960 stocks meet this criterion and are included in the out-of-sample analysis (marked in boldface in the Web Appendix \textit{List of Tickers}).
The model universe considered for both the full-sample and out-of-sample analysis includes all the specifications presented in Section \ref{s:mods}. As discussed, each of these is coupled with three different $\mathsf{ES}$\ specifications, \texttt{ES_sim}, \texttt{ES_sk}, and \texttt{ES_no}, for a total of 18 different models. Implementing the procedures discussed in Section \ref{s:esti}, each model is estimated for three different risk levels, $\alpha \in \{0.01,0.025,0.05\}$.
\subsection{Analysis of the in- and out-of-sample losses and coverage}
Before delving into the assessment of the model performances through the various tests conducted on the extensive dataset, it is first useful to visually assess their in- and out-of-sample coverage.
Figure \ref{fig:insampleCov} presents a comprehensive overview of the aggregated information across all datasets, three coverage levels, and all models for in-sample performance. Each model is estimated for every dataset, and the in-sample empirical coverage ($\hat\alpha$) is calculated. The models are represented by row blocks in the figure, with corresponding names on the y-axis, such as ``\texttt{VaR=m_lev, ES=no, Loss=EM}''. Within each block, three box plots display the coverage for all datasets, with blue indicating $\alpha = 0.01$, green representing $\alpha = 0.025$, and red representing $\alpha = 0.05$. These levels are also depicted by vertical dashed lines. All the models exhibit similar behavior and, on average, achieve the desired coverage level, albeit with slight variations. Simpler models generally exhibit less variability. Some cases encountered convergence difficulties, leading to the inability to estimate certain models. The right panel of Figure \ref{fig:insampleCov} shows the fraction of such problematic cases, consistently below 3\%.
A similar analysis has been conducted for out-of-sample coverage, utilizing three different window sizes of 500, 1000, and 2000 days. Aggregated results are presented in Figure \ref{fig:outofsampleCov}. In addition to the three colors representing coverage levels ($\alpha = 0.01$, $\alpha = 0.025$, and $\alpha = 0.05$), varying color intensities indicate window size (lightest shade = 2000 days; darkest = 500 days).\footnote{Fewer models are considered in the out-of-sample analysis, excluding "\texttt{ES=SkKu}" due to computational complexity and "\texttt{Loss=FZ0}" due to its similar behavior to "\texttt{Loss=ALS}".} The out-of-sample results reveal a less favorable situation than the in-sample, as all models tend to overestimate the coverage on average, less severely so with a wider rolling window. Surprisingly, the variance also increases in this scenario. This can be attributed to longer intervals containing more diverse data from potentially different underlying models, thus imperfectly capturing future behavior. Despite these nuances, all models demonstrate similar behavior based on simple visual inspection. Additionally, Figure \ref{fig:outofsampleLosses} provides aggregated loss information. It is evident that both the values and spreads of the loss function decrease with larger sample sizes.
\subsection{Evaluation metrics}
Following the practice by researchers and risk managers, the in- and out-of-sample performances of the dynamic models\footnote{It should be noted that performing a complete estimation for all rolling windows in the out-of-sample exercise has been highly time-consuming. As a result, we perform a full estimation for the initial window and subsequently at 500-day time intervals, and, instead of repeating the complete estimation process for each subsequent window, we update the parameters at regular intervals. To accomplish this, we perform parameter optimization every 50 observations, starting from the results obtained in the previous step. This parameter updating allows us to refine the estimation without repeating the entire process. In all other rolling windows, we maintain the parameters obtained from the previous window.} for $\mathsf{VaR}$\ and $\mathsf{ES}$\ presented in Section \ref{s:mods} are firstly assessed using some diagnostic tests, whose technical details are summarized in Appendix \ref{s:hit}-\ref{s:patton} for the reader's convenience.
To assess the in-sample $\mathsf{VaR}$\ estimation performance, we consider the in-sample Dynamic Quantile (DQ) test by \cite{caviar}. In particular, we consider the test in its conditional coverage and independence
versions as described in \cite{dumietal2012}. The asymptotic theory for these tests was originally derived by \cite{caviar} for the pure CAViaR models. Therefore, the results presented hereafter refer to the \texttt{ES_no} case only, involving CAViaR models estimated by minimizing the aggregated quantile loss.
While the in-sample DQ test assesses the goodness-of-fit of CAViaR models, its out-of-sample counterpart can be seen as a general test for evaluating the statistical properties of a set of $\mathsf{VaR}$\ forecasts, regardless of the model. This includes testing for unbiasedness, independent hits, and the independence of quantile estimates, as outlined by \cite{caviar}. In our case, the out-of-sample DQ (OOS-DQ) test can effectively evaluate the properties of the $\mathsf{VaR}$\ forecasts generated by joint dynamic $\mathsf{VaR}$-$\mathsf{ES}$\ models.
Further, we jointly assess the statistical accuracy of $\mathsf{VaR}$\ and $\mathsf{ES}$\ estimates, both in and out-of-sample, employing two regression-based testing procedures, i.e., the regression-based calibration tests by \cite{pattonetal2019}, henceforth PZC, and the ESR test by \cite{bay_dim_2020}. The former includes separate calibration tests for $\mathsf{VaR}$\ and $\mathsf{ES}$\ while the latter test is specific for $\mathsf{ES}$\ diagnostics. Moreover, the PZC tests are based on OLS auxiliary regression equations where the standardized generalized residuals \citep[as in][]{pattonetal2019} are regressed on their past values as well as on $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasts, respectively.
The ESR approach, instead, is based on three separate test statistics: the Auxiliary, the Strict and the Strict Intercept Backtest, which can be seen as an extension of Mincer-Zarnowitz regression to a semi-parametric setting, relying on the minimization of the consistent loss functions proposed by \cite{Fissler2016}. The Auxiliary and Strict test statistics are computed regressing returns on the $\mathsf{ES}$\ forecasts and test the $\mathsf{ES}$\ coefficients for joint $(0, 1)$ values. Specifically, the Auxiliary test requires an auxiliary $\mathsf{VaR}$\ forecast, while the Strict Intercept tests whether the expected shortfall of the forecast error $(r_t-ES_t)$ is zero.\footnote{\cite{bay_dim_2020} also consider a one-sided version of the Strict Intercept that is particularly useful for regulatory evaluations. Since our main interest is simply in the assessment of forecasting accuracy, in this paper we only consider the two-sided version of the test.}
Finally, the out-of-sample forecasting accuracy of each of the models is assessed by comparing the average values of the FZ0 loss achieved over the forecasting period.
\subsection{Testing the properties of the risk estimates}
In this section, we analyze the properties of the in-sample risk estimates over the full-sample for the set of 823 assets, postponing to a later subsection the generation of out-of-sample risk forecasts. Also, for the sake of brevity, we only discuss results for $\alpha=0.01$, while the results for $\alpha=0.025$ and $\alpha=0.05$ are contained in the Web Appendix \textit{Tables W.1 and W.2}.
The \emph{In-sample DQ} section of Table \ref{tab:1} reports the results of the in-sample DQ test in its conditional coverage ($CC-DQ_{IS}$) and independence ($ID-DQ_{IS}$) versions, respectively. For each model, the table provides the non-rejection frequency at the $5\%$ significance level. Higher values in the table indicate better performance, as they correspond to models less frequently rejected by the tests.
Although the $CC-DQ_{IS}$ provides a comprehensive evaluation of the risk estimation performance, it has a \textit{portmanteau} nature, which overlooks the clustering features of the hit series and neglects their coverage properties. By contrast, the $ID-DQ_{IS}$ test offers a complementary perspective to the previous test, since it focuses explicitly on clustering. Combining the information from both tests makes it possible to gain deeper insight into the reasons behind any model underperformance.
The test findings for $\alpha=0.01$ (similar results hold for the other risk levels) can be summarized as follows:
\begin{itemize}
\item[] $CC-DQ_{IS}$ : multiplicative models outperform their additive counterparts with the \texttt{mlt_lev} resulting the best model at all risk levels. This model is not rejected at the 5\% level in approximately 70\% of cases, closely followed by the \texttt{mlt_skk}. The \texttt{mlt_sim} yields slightly lower rates than models incorporating information on realized skewness and kurtosis. The non-rejection frequency of additive models is much lower, being on average close to $20\%$, with the highest rate being recorded for the \texttt{add_sim} model.
\item[] $ID-DQ_{IS}$: multiplicative and additive models are characterized by similar performances, suggesting that the high rejection rate of the latter class is mostly due to lack of coverage rather than to hit clustering.
\end{itemize}
Second, to appreciate the contribution of the additional information in the form of realized skewness and kurtosis in the \texttt{skk} models, in the \emph{Wald test} section of Table \ref{tab:1} we assess the significance of the skewness and kurtosis coefficients involved in the $\mathsf{VaR}$\ and $\mathsf{ES}$\ dynamics, respectively. Specifically, we test the null $a_1 = a_2 = a_3=0$, for $\mathsf{VaR}$, and $b_1 = b_2 = 0$, for $\mathsf{ES}$, against a two-sided alternative.
Again, the test results in terms of empirical non-rejection frequencies are summarized over the panel of assets considered. For the plain CAViaR models, the null $a_1 = a_2 = a_3 =0$ is almost always rejected at the usual 5\% significance level, for all risk levels considered. Differently, for joint $\mathsf{VaR}$-$\mathsf{ES}$\ models, the non-rejection frequency increases with the risk level $\alpha$. Namely, the percentage of non-rejections is close to 0 for $\alpha = 0.01$ but it increases to values up to $\approx 40\%$ for $\alpha = 0.05$.
The discrepancy between non-rejection frequencies for pure $\mathsf{VaR}$\ and joint $\mathsf{VaR}$-$\mathsf{ES}$\ models is likely to be due to the fact that, for each class of models, testing is based on a different asymptotic distribution: we rely on the theory derived by \cite{caviar}, for the EM loss, and on \cite{pattonetal2019}, for ALS and FZ0. The test results are only marginally affected by the choice of the joint loss, AL or FZ0, used for estimation.
Moving to the analysis of $\mathsf{ES}$\ dynamics, we find that the non-rejection frequencies of the null $b_1 = b_2 = 0$ are substantially higher than the values observed for $\mathsf{VaR}$\ parameters and are clearly affected by the risk level. Namely, they approximately lie in the range 49\%-61\%, for $\alpha=0.01$, 55\%-72\%, for $\alpha=0.025$, and 62\%-79\%, for $\alpha=0.05$. Results are very close for models based on ALS and FZ0 losses. Overall, we conclude that the inclusion of realized skewness and kurtosis measures in the $\mathsf{ES}$\ equation is less strongly supported than for the $\mathsf{VaR}$.
Next, we focus on the in-sample PZC and ESR calibration tests. First, the \emph{Calibration test} section of Table \ref{tab:1} for $\alpha = 0.01$ (see Web Appendix \textit{Tables W.1 and W.2} for $\alpha= 0.025$ and $\alpha= 0.05$) reports the results of the former tests for $\mathsf{VaR}$\ and $\mathsf{ES}$. The main findings arising from the table can be summarized as follows:
\begin{itemize}
\item for $\alpha \geq 0.025$ (Tables W.1 and W.2), all models yield remarkably good non-rejection frequencies, with values ranging from 72\% to 94\%.
\item For both $\mathsf{VaR}$\ and $\mathsf{ES}$, we record a decay of the non-rejection frequency at the 0.01 risk level (Table \ref{tab:1}). This is particularly relevant for $\mathsf{VaR}$\ since $\alpha=0.01$ is the mandatory level indicated by the Basel Committee.
\item Models based on ALS and FZ0 losses return very close performances.
\item Comparing simpler models (\texttt{*_sim}) with more complicated specifications (\texttt{*_skk} and \texttt{*_lev}), we find that there is no clear winner but the ranking depends on the functional form and risk level.
\end{itemize}
Finally, to assess the "calibration" of $\mathsf{ES}$\ forecasts, the \emph{ES calibration test} section in Table \ref{tab:1} for $\alpha = 0.01$ (see Web Appendix \textit{Tables W.1 and W.2} for $\alpha= 0.025$ and $\alpha= 0.05$) reports the non-rejection frequencies of the three ESR tests proposed by \cite{bay_dim_2020}\footnote{The tests were implemented using the \texttt{esback} \texttt{R}-library provided by the same authors, freely available from CRAN at the URL: https://cran.r-project.org/web/packages/esback/index.html}. It is worth noting that due to numerical problems in the computation of the test statistic, this could not be computed for some of the assets in our panel, in addition to those that had been previously excluded due to convergence issues in the estimation of the reference risk models: the number of valid assets for each configuration, determined by a combination of available models and risk levels, ranges from a minimum of 644 to a maximum of 796 assets out of 823.
Compared to the calibration test by \cite{pattonetal2019}, the ESR reveals a much lower discriminatory power returning non-rejection frequencies very close to unity for all models and risk levels. Again, we do not report any apparent differences in model performances based on the ALS and FZ0 losses.
\subsection{Out-of-sample forecasting comparison}
\label{sec:oos}
This section presents the results of the out-of-sample forecasting analysis. First, the performance of the models under analysis is assessed by computing the following test statistics and diagnostics over the out-of-sample period
\begin{itemize}
\item DQ tests for independence and conditional coverage
\item $\mathsf{VaR}$\ and $\mathsf{ES}$\ calibration tests by \cite{pattonetal2019}
\item ESR tests for $\mathsf{ES}$\ calibration by \cite{bay_dim_2020}.
\end{itemize}
As in the previous section, test results across the whole panel of assets are summarized in terms of empirical non-rejection frequencies. Also, we only discuss results for $\alpha=0.01$ in Table \ref{tab:2} while results for $\alpha=0.025$ and $\alpha=0.05$ have been reported in the Web Appendix \textit{Tables W.3 and W.4}.
Similarly to what was observed in the full sample analysis, in a limited number of cases it has not been possible to calculate the $p$-values of ESR tests due to failures in the estimation of the auxiliary regression model underlying the test. Overall, depending on risk level, specific test of interest, and sample size, the available number of stocks has been found to range between 371 and 406 out of 406 potentially available stocks.
The findings of the analysis can be succinctly summarized as:
\begin{itemize}
\item DQ tests: the non-rejection frequencies are very low for the shortest estimation window $T=500$ but they tend to increase with the sample size although, even for $T=2000$, they barely exceed $40\%$, for independence tests, only in a few isolated cases. Overall some stylized facts arise. Plain $\mathsf{VaR}$\ models on average perform better than joint $\mathsf{VaR}$-$\mathsf{ES}$\ models while the inclusion of information on skewness and kurtosis does not bring any evident advantages.
\item $\mathsf{VaR}$\ calibration tests: the performances are very poor for the shortest sample size $T=500$ but tend to improve as $T$ increases. The model performances also depend on the value of the risk level $\alpha$ with the best results obtained for $\alpha=0.025$. In terms of model specifications, \texttt{add_sim} and \texttt{mlt_sim} yield the highest non-rejection frequencies that exceed $70\%$ for T=2000 and $\alpha=0.025$ when the EM loss is used. When comparing plain $\mathsf{VaR}$\ and joint $\mathsf{VaR}$-$\mathsf{ES}$\ models, there are no clear performance gaps.
\item $\mathsf{ES}$\ calibration tests: the results are qualitatively not different from what was observed for the $\mathsf{VaR}$\ tests. Hence, similar considerations hold.
\item ESR tests: the performance of the ``strict'' and ``auxiliary'' tests improves as the sample size increases although the performance gap across different sample sizes is less evident than for the other regression-based $\mathsf{VaR}$\ and $\mathsf{ES}$\ calibration tests. As above, even in this case, we record the best performances for the \texttt{add_sim} and \texttt{mlt_sim} models reaching, in some cases, non-rejection frequencies close to $80\%$. As far as the ``strict intercept test'' is concerned, the differences across different models and sample sizes are much less evident and the non-rejection frequency is $>80\%$ in all instances. Again, the information on realized skewness and kurtosis does not appear to lead to improvements in terms of forecasting performances.
\end{itemize}
Finally, we assess and compare the forecasting accuracy of the different models based on the out-of-sample values of the following strictly consistent scoring functions: quantile loss (E) for $\mathsf{VaR}$\ and AL-score for joint $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasting (ALS).
In terms of median loss (Table \ref{tab:medianloss}), the multiplicative model without skewness and kurtosis information (\texttt{mlt\_sim}) achieves the minimum loss value in most cases for both quantile and ALS scoring functions. It is only slightly outperformed by its additive counterpart (\texttt{add\_sim}) in one instance for pure $\mathsf{VaR}$\ models and in two instances for joint forecasts of $\mathsf{VaR}$\ and $\mathsf{ES}$. A similar trend is observed when considering average ranks (Table \ref{tab:medianranks}). The \texttt{mlt\_sim} model consistently delivers the minimum average rank, except in the case of joint $\mathsf{VaR}$\ and $\mathsf{ES}$\ forecasts at the 0.05 level and for $T=1000$, where it ranks second, closely following the \texttt{add\_sim} model that also does not use skewness and kurtosis information.
In conclusion, the key insights from the assessment of forecasting performance can be summarized as follows:
\begin{itemize}
\item Incorporating information on realized skewness and kurtosis does not enhance forecasting accuracy;
\item Simpler models are preferable to more complex ones, as the latter are more vulnerable to computational issues;
\item The \emph{multiplicative} specification is generally preferable to the more popular \emph{additive} approach.
\end{itemize}
\section{Concluding remarks}
\label{s:conc}
In this paper, we have presented a forecasting comparison of several semi-parametric risk forecasting models. Our work presents some important elements of novelty and potential interest for practitioners and researchers alike. First, the comparison is based on an unusually large set of 823 stocks: to the best of our knowledge, there are no other contributions relying on such a large dataset in the tail-risk forecasting literature. Also, the availability of such a rich data environment has a positive impact on the reliability of the regularities that emerge from the empirical analysis, giving them a good degree of external validity.
Second, we assess the potential contribution coming from considering information on some recently proposed realized skewness and kurtosis measures. Third, we provide deeper insight into the selection of the functional form of the semi-parametric model used to generate forecasts.
The results of our analysis clearly indicate that, at the forecasting stage, simple models should be preferred to more complicated ones with a preference for multiplicative GARCH-type specifications. Realized skewness and kurtosis measures do not apparently provide valuable information for improving the accuracy of tail risk forecasts even if in most cases, their coefficients turn out to be significant in the full-sample analysis. By the same token, they may prove useful in generating improved density forecasts, a task that we leave for future research.
When we shift the focus to the functional form of the dynamic risk model, an interesting and original finding from our extensive empirical investigation is that the standard \emph{CaViaR-like} additive model specification outperformed by the less commonly used (in a semi-parametric framework) \emph{GARCH-like} multiplicative parameterization.
\bigskip
\bibliographystyle{apalike}
\bibliography{literature}
\newpage