Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
56,480 characters · 12 sections · 60 citation commands
Dynamic tail risk forecasting: what do realized skewness and kurtosis add?
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 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 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 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 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), we propose our models in Subsection (ref), while the related estimation procedures and some properties of the estimators are illustrated in Subsection (ref). In Section (ref), we present a recent literature review on realized estimators of skewness and kurtosis of financial returns. Section (ref) is dedicated to the empirical application, while Section (ref) concludes.
The literature on semi-parametric risk modeling features a seminal paper by 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, 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, 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 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, taylor2017 considers a loss, or scoring, function based on the Asymmetric Laplace quasi-likelihood function, AL for short. pattonetal2019 extend the work by 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 crealetal2013. Second, as done by 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 diebold1995comparing.
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
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), 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 additive and the multiplicative models, respectively. We simplify the notation by defining $v_{t}\equiv Q_{\alpha,t}$ and $e_{t}\equiv ES_{\alpha,t}$. Thus, the additive modeling framework can be represented as a regression model for the 1-step ahead expected $\alpha$-level of $\mathsf{VaR}$, $v_{t}$:
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:
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 ande01_jasa. A set of alternative choices for $\widehat\operatorname{RV}_{t-1}$ will be presented and discussed in Section (ref).}.
To account for the potential misspecification in add_sim, we can resort to the Cornish-Fisher (CF) expansion 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
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 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.}:
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 add_skk with an asymmetric impact of the integrated volatility in correspondence with returns smaller than their average (leverage effect):
By contrast, a multiplicative modeling framework can be represented in terms of the following nonlinear regression model
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
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
When we allow for a time-varying conditional skewness and kurtosis in the returns distribution, this assumption must be generalized to read
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 mlt_lev and 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}$
Multiplicative models closely mirror additive models ((ref)), ((ref)) and ((ref)). The mlt_sim is similar to ((ref)) and is the simplest specification with only the integrated volatility driving the dynamics of the scale. Further mlt_skk assumes similar to ((ref)) in the additional information incorporated in realized skewness and kurtosis that drives the dynamics of the scale. The most complex model mlt_lev also controls for the leverage in $h_t$ in the same fashion as in the model ((ref)).
For both additive and multiplicative frameworks, the $\mathsf{ES}$\ can be modeled according to two different alternative specifications:
Here ES_sim is the simple specification assuming that the $\mathsf{ES}$\ is a rescaling of $\mathsf{VaR}$. tayl2019 shows that this simple specification provides competitive $\mathsf{VaR}$\ forecasts. More recently, wang_gerl_chen_2023 have extended the framework proposed in tayl2019 to allow for separate $\mathsf{VaR}$\ and $\mathsf{ES}$\ dynamics as well as for the incorporation of realized measures.
The more complex ES_skk brings the dynamics of the $\mathsf{ES}$\ to be also driven by the skewness and kurtosis, possibly accounting for the misspecification of ES_sim. Differently from $\mathsf{VaR}$, in this case, we did not split the skewness into negative and positive. By construction, both ES_skk and 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 ES_no, corresponds to a model specification that is close in spirit to a caviar CAViaR type one where some realized estimator of the integrated variance replaces the volatility measure based on lagged daily returns.
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))-((ref)) is done semi-parametrically by minimizing a strictly consistent scoring rule,
where $S^{(\alpha)}_t$ is a member of the general class presented by Fissler2016, i.e.
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 tayl2019
Adding and subtracting $\alpha r_t$ to the numerator of the second term on the right-hand side of the previous equation, we get
In the simplified expression of $AL^{(\alpha)}_t$ obtained by 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 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 pattonetal2019
Therefore, in the case of the model EM_no, where only the $\mathsf{VaR}$\ is estimated, the objective function is given by the quantile loss:
In what follows, models estimated relying on loss functions $AL^{(\alpha)}_t$, $FZ0^{(\alpha)}_t$ and $EM^{(\alpha)}_t$ are labeled as Loss=ALS, Loss=FZ0 and Loss=EM, respectively. Standard errors are computed using the asymptotic theory developed by caviar, for pure $\mathsf{VaR}$\ models, and pattonetal2019, for joint $\mathsf{VaR}$-$\mathsf{ES}$\ models. Technical details are provided in Section (ref).
It is worth noting that optimization in the ((ref)) 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 caviar which, for ease of reference, we call 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.
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:
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 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 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 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 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 AndersenBoll98,
other widely used estimators of the integrated volatility are the, proposed in BarndorffNielsen_Bipower_2004 and 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 Barndorff_Semivariance_2008 and 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
proposed by andersen2012jump, because of its documented robustness properties.
Several new estimators for 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 MeiLiuMaChen2017. The simplest estimator of the integrated $k$-order moments was proposed by amaya2015, shadowing the relationship between the realized variance and the integrated volatility (case $k = 2$).
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 LiuWangLiu2014 derived asymptotic properties of the amaya2015 estimator and developed their own measures of realized skewness accounting for market microstructure noise. Another extension was provided by 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 neuberger2013 and neuberger2020:
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 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.
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 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 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 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). As discussed, each of these is coupled with three different $\mathsf{ES}$\ specifications, ES_sim, ES_sk, and ES_no, for a total of 18 different models. Implementing the procedures discussed in Section (ref), each model is estimated for three different risk levels, $\alpha \in \{0.01,0.025,0.05\}$.
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) 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 “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) 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). 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 "ES=SkKu" due to computational complexity and "Loss=FZ0" due to its similar behavior to "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) provides aggregated loss information. It is evident that both the values and spreads of the loss function decrease with larger sample sizes.
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) are firstly assessed using some diagnostic tests, whose technical details are summarized in Appendix (ref)-(ref) for the reader's convenience.
To assess the in-sample $\mathsf{VaR}$\ estimation performance, we consider the in-sample Dynamic Quantile (DQ) test by caviar. In particular, we consider the test in its conditional coverage and independence versions as described in dumietal2012. The asymptotic theory for these tests was originally derived by caviar for the pure CAViaR models. Therefore, the results presented hereafter refer to the 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 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 pattonetal2019, henceforth PZC, and the ESR test by 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 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 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{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.
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 Tables W.1 and W.2.
The In-sample DQ section of Table (ref) 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 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:
Second, to appreciate the contribution of the additional information in the form of realized skewness and kurtosis in the skk models, in the Wald test section of Table (ref) 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 caviar, for the EM loss, and on 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 Calibration test section of Table (ref) for $\alpha = 0.01$ (see Web Appendix 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:
Finally, to assess the "calibration" of $\mathsf{ES}$\ forecasts, the ES calibration test section in Table (ref) for $\alpha = 0.01$ (see Web Appendix 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 bay_dim_2020\footnote{The tests were implemented using the esback 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 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.
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
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) while results for $\alpha=0.025$ and $\alpha=0.05$ have been reported in the Web Appendix 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:
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)), the multiplicative model without skewness and kurtosis information (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 (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)). The 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 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:
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 CaViaR-like additive model specification outperformed by the less commonly used (in a semi-parametric framework) GARCH-like multiplicative parameterization.