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.
102,526 characters · 24 sections · 103 citation commands
Realized Stochastic Volatility Model with Skew-t Distributions for Improved Volatility and Quantile Forecasting
{\it Keywords:} Bayesian estimation, Markov chain Monte Carlo, Realized volatility, Skew-t distribution, Stochastic volatility
Volatility, defined as the conditional standard deviation or variance of asset returns, evolves stochastically over time and plays a central role in financial risk management. Accurate volatility forecasting is essential for evaluating tail risks such as value-at-risk (VaR) and expected shortfall (ES). Traditionally, two major classes of time-series models have been used for modeling time-varying volatility: the generalized autoregressive conditional heteroskedasticity (GARCH) family engle_autoregressive_1982, bollerslev_generalized_1986, and the stochastic volatility (SV) model taylor_modelling_1986. These models successfully capture key features of financial volatility, including volatility clustering and persistence, and have been extended to account for leverage effects---i.e., the negative correlation between current returns and future volatility. A prominent example is the exponential GARCH (EGARCH) model by nelson_conditional_1991.
More recently, realized volatility (RV), calculated as the sum of squared intraday returns, has gained popularity as a nonparametric and more efficient measure of daily volatility. Comprehensive surveys on RV are provided by andersen_realized_2009 and mcaleer_realized_2008. To model the time-series behavior of RV, researchers have proposed long-memory models such as autoregressive fractionally integrated moving average (ARFIMA) models beran_statistics_1994, as well as the heterogeneous autoregressive (HAR) model corsi_simple_2009, which approximates long-memory dynamics using a small number of lags.
Although ARFIMA and HAR models often outperform traditional GARCH and SV models in volatility forecasting, RV estimates are prone to biases caused by market microstructure noise and non-trading hours. Several techniques have been developed to address these issues, including multiscale estimators zhang_tale_2005, zhang_efficient_2006, pre-averaging methods jacod_microstructure_2009, and realized kernel (RK) estimators barndorff-nielsen_designing_2008, barndorff-nielsen_realized_2009. For further refinements and applications, see ait-sahalia_estimating_2009, ubukata_pricing_2014, and liu_does_2015.
To exploit the advantages of both parametric and nonparametric approaches, hybrid models have been proposed. These include the realized stochastic volatility (RSV) model takahashi_estimating_2009, dobrev_information_2010, koopman_analysis_2013, as well as the realized GARCH (RGARCH) and realized EGARCH (REGARCH) models hansen_realized_2012, hansen_exponential_2016. While several studies have compared forecasting performance within each class, comprehensive cross-model evaluations remain limited. Recent studies such as takahashi_stochastic_2023, takahashi_forecasting_2024 have attempted to bridge this gap.
Beyond volatility, accurate prediction of financial tail risks requires appropriate modeling of return distributions, which often exhibit skewness and leptokurtosis. While time-varying volatility captures some of this behavior, the conditional distribution of returns may remain non-Gaussian even after accounting for volatility. To address this, skewed Student’s \(t\)-distributions have been widely applied in volatility models abanto-valle_bayesian_2015, kobayashi_skew_2016. Among them, the generalized hyperbolic (GH) skew-\(t\) distribution has shown strong empirical performance in SV nakajima_stochastic_2012, leao_bayesian_2017 and RSV models trojansebastian_regime_2013, nugroho_realized_2014, nugroho_boxcox_2016, takahashi_volatility_2016.
Score-driven (SD) models, introduced by creal_generalized_2013 and harvey_dynamic_2013, have also gained popularity as a flexible alternative to SV and GARCH.\footnote{An extensive and regularly updated list of SD models proposed in the literature is available at \url{https://www.gasmodel.com/}.} catania_forecasting_2022 apply the GH skew-\(t\) distribution within an SD--GARCH setting, while catania_forecasting_2020 develop a joint SD model for returns and realized volatility using a time-varying bivariate \(t\) specification. These approaches are conceptually related to ours, but their dynamic updating rules and likelihood evaluations are substantially more complex. We therefore focus on SV/RSV-type models, while acknowledging SD-based specifications as an important complementary direction in the broader volatility-modeling literature.
In this paper, we extend the RSV framework by incorporating three types of skew-\(t\) distributions: the Azzalini-type skew-\(t\) distribution azzalini_class_1985, the Fern\'{a}ndez and Steel (FS) skew-\(t\) distribution fernandez_bayesian_1995, and the GH skew-\(t\) distribution aas_generalized_2006. We adopt a Bayesian estimation strategy using Markov chain Monte Carlo (MCMC) simulation and apply the models to daily returns of the Dow Jones Industrial Average (DJIA) and the Nikkei 225 (N225).
Our empirical analysis evaluates the performance of RSV models alongside benchmark SV, EGARCH, and REGARCH models. Volatility forecasts are assessed using the quasi-likelihood loss function, while VaR and ES forecasts are evaluated using the joint loss function proposed by fissler_higher_2016 and specified in patton_dynamic_2019. To test forecast performance, we implement the predictive ability test of giacomini_tests_2006 and the model confidence set procedure by hansen_model_2011. The results show that incorporating realized volatility substantially improves forecasting performance. In particular, RSV models with skew-\(t\) innovations exhibit superior predictive accuracy, confirming and extending findings from prior research.
The remainder of the paper is organized as follows. Section (ref) presents the SV and RSV models, including skewed innovations and a brief overview of realized measures. Section (ref) outlines the Bayesian estimation method and forecast evaluation metrics. Section (ref) reports the empirical results. Section (ref) concludes with a summary of findings and implications.
Let \( y_t \) denote the daily log return of an asset: \[ y_t = \log p_t - \log p_{t-1}, \] where \( p_t \) is the closing price on day \( t \). The conditional return variance is modeled as \[ \mathrm{E}[y_t^2 \mid \mathcal{I}_{t-1}] = \sigma_t^2 = \exp(h_t), \] where $\mathcal{I}_{t-1}$ denotes the information set available at time $t-1$ and \( h_t \) denotes the conditional log volatility, treated as a latent process.
The log volatility \( h_t \) is assumed to follow a stationary AR(1) process:
Here, \( \phi \) is the persistence parameter, and \( \mu \) is the unconditional mean of log volatility. Stationarity is ensured by \( |\phi| < 1 \). The initial state \( h_1 \) is drawn from its unconditional distribution.
The bivariate error term \( (\epsilon_t, \eta_t)' \) follows a jointly normal distribution with correlation \( \rho \). This structure allows for a leverage effect: a negative return shock (\( \epsilon_t < 0 \)) tends to increase future volatility (\( h_{t+1} \)), reflected by \( \rho < 0 \). This feature captures the empirically observed asymmetry in financial markets, where volatility tends to rise following negative returns.
In the SV model, volatility \( \exp(h_t) \) is treated as a latent variable, since the log-volatility \( h_t \) is unobserved. With the increasing availability of high-frequency data, RV has emerged as a prominent nonparametric estimator of true volatility.
Let \( p(s) \) denote the logarithm of the asset price at time \( s \), and suppose it follows the continuous-time diffusion process:
where \( \mu(s) \) and \( \sigma(s)^2 \) are the drift and instantaneous variance, respectively, and \( W(s) \) is a standard Brownian motion. The integrated volatility (IV) over a single day is defined as
Assume that \( m \) intraday returns are observed on day \( t \), denoted by \( \{r_{t-1+1/m}, r_{t-1+2/m}, \ldots, r_{t}\} \). Then, the realized volatility is defined as
Under ideal conditions (e.g., no microstructure noise or missing observations), \( RV_t \) converges to \( IV_t \) as \( m \to \infty \).
However, these ideal conditions are often violated in practice. One source of bias arises from non-trading hours, which are not captured in intraday data. Ignoring these periods leads to an underestimation of daily volatility. To address this, hansen_forecast_2005 proposed a scaling adjustment:
This adjustment ensures that the average of the scaled RV matches the sample variance of daily returns.
Another significant bias stems from microstructure noise. Observed prices deviate from the efficient price due to bid--ask bounce, price discreteness, and asynchronous trading, introducing autocorrelation and distortions in high-frequency returns. The bias tends to worsen as the sampling interval becomes shorter. Optimal sampling frequency and alternative estimators have been studied extensively ait-sahalia_how_2005, bandi_separating_2006, bandi_microstructure_2008, liu_does_2015.
Among various alternatives, the RK estimator barndorff-nielsen_designing_2008 has been widely adopted due to its robustness. It is defined as:
where \( H \) is the bandwidth, \( k(\cdot) \) is a weight function, and
barndorff-nielsen_realized_2009 provide guidelines for choosing the bandwidth \( H \).
To account for these biases in modeling, hybrid frameworks such as the RSV and RGARCH models have been developed. In particular, the RSV model of takahashi_estimating_2009 addresses RV bias by jointly estimating model parameters, as detailed in Section (ref).
takahashi_estimating_2009 proposed enhancing the SV model by incorporating additional information from RV. In their framework, the daily log-return \( y_t \) and the logarithm of realized volatility \( x_t = \log RV_t \) are jointly modeled, with the latent log volatility \( h_t = \log IV_t \) serving as a common driver.
The RSV model is specified as:
The error term \( u_t \) is assumed to be independent of \( (\epsilon_t, \eta_t) \).
Equation (ref) may be generalized as \( x_t = \xi + \psi h_t + u_t \), but empirical evidence suggests that estimating \( \psi \ne 1 \) does not improve forecasting performance. Accordingly, we fix \( \psi = 1 \), following takahashi_estimating_2009.
The parameter \( \xi \) adjusts for systematic bias in \( x_t = \log RV_t \). When \( \xi = 0 \), \( x_t \) is an unbiased estimate of \( h_t \) . In practice, however, \( RV_t \) may be biased downward due to non-trading hours and upward due to microstructure noise. hansen_realized_2006 show that the net bias can be either positive or negative. Therefore, the sign of \( \xi \) reflects the relative magnitude of these two opposing effects.
The distribution of \( \epsilon_t \) captures residual return behavior after controlling for time-varying volatility. While the SV structure already accounts for volatility clustering and excess kurtosis, it is well documented that financial returns often exhibit additional heavy tails and asymmetry. To model these features more effectively, several studies have incorporated non-Gaussian innovations.
For the SV model, nakajima_stochastic_2012 and leao_bayesian_2017 employed the GH skew-\(t\) distribution from aas_generalized_2006. Other approaches include the skew-\(t\) distribution by azzalini_distributions_2003 abanto-valle_bayesian_2015, the skew exponential power distribution kobayashi_skew_2016, and the flexible skew-\(t\) model by fernandez_bayesian_1995, as adopted in steel_bayesian_1998. In the context of RSV models, trojansebastian_regime_2013, nugroho_realized_2014, nugroho_boxcox_2016, and takahashi_volatility_2016 applied the GH skew-\(t\) distribution.
Building on this literature, we extend the RSV model defined in equations (ref)--(ref)---hereafter referred to as the RSV-N model---by incorporating skewed and heavy-tailed distributions for \( \epsilon_t \). This extension allows for more accurate modeling of return asymmetry and tail risk, which is essential for reliable risk forecast applications.
We extend the RSV-N model by allowing the return innovation \( \epsilon_t \) to follow a standardized Student's \( t \) distribution. The resulting model, referred to as the RSV-T model, modifies equations (ref)--(ref) as follows:
where
Here, \( \lambda_t \) is an inverse-gamma mixing variable, and \( \nu \) is the degrees of freedom parameter controlling tail thickness. As \( \nu \to \infty \), the model converges to the Gaussian RSV-N specification.
The use of Student’s \( t \) distribution allows the model to accommodate excess kurtosis beyond what is captured by the stochastic volatility component \( h_t \). For finite \( \nu \), the distribution of \( \epsilon_t \) exhibits heavier tails, improving the model’s ability to reflect extreme return events.
Importantly, the contemporaneous correlation between \( \epsilon_t \) and \( \eta_t \) becomes a function of both \( \rho \) and \( \nu \), given by:
where
Note that the variance of \( \epsilon_t \) exists only when \( \nu > 2 \), and its fourth moment exists when \( \nu > 4 \). In empirical applications, \( \nu \) is typically estimated from the data.
takahashi_volatility_2016 extended the RSV-N model in equations (ref)--(ref) by incorporating the GH skew-\( t \) distribution, resulting in the RSV-GH-ST model. This is achieved by replacing \( \epsilon_t \) and \( \eta_t \) in (ref)--(ref) with the following mixture representation:
where \( \lambda_t \sim \mathcal{IG}(\nu/2, \nu/2) \), \( \mu_\lambda = \nu/(\nu - 2) \), and
A distinctive feature of the GH skew-\( t \) distribution, originally emphasized by aas_generalized_2006, is its asymmetric tail behavior: one tail exhibits polynomial (power-law) decay, while the opposite tail decays exponentially, depending on the sign of the skewness parameter \( \beta \). This feature makes the GH skew-\( t \) distribution particularly flexible in capturing return asymmetry and tail risk in financial applications.
The denominator in (ref) standardizes \( \epsilon_t \) so that its conditional variance equals \( \exp(h_t) \). In the present formulation, the skewness and tail behavior of \( \epsilon_t \) are jointly governed by \( \beta \) and \( \nu \). The parameter \( \beta \) introduces asymmetry and determines which tail becomes heavier, while \( \nu \) controls the degree of polynomial tail thickness. When \( \beta = 0 \), the model reduces to the RSV-T specification in . As \( \nu \to \infty \), or if \( \lambda_t = 1 \) for all \( t \), the distribution converges to the Gaussian case regardless of \( \beta \).
The correlation between \( \epsilon_t \) and \( \eta_t \) depends on \( \rho \), \( \nu \), and \( \beta \), and is given by:
where \( \mathrm{E}[ \sqrt{\lambda_t} ] \) is as defined in equation (ref).
To illustrate the influence of \( \beta \) and \( \nu \), Figure (ref) plots the simulated densities of \( \epsilon_t \) under various parameter values. Panel (i) shows that increasing negative values of \( \beta \) yield more pronounced left-skewness and heavier polynomial tails. In contrast, Panel (ii) indicates that as \( \nu \) increases, the distribution becomes more symmetric and approaches normality, exhibiting thinner tails.
We extend the RSV-N model defined in equations (ref)--(ref) by adopting the skew-\(t\) distribution proposed by azzalini_class_1985.\footnote{See also azzalini_distributions_2003 and sahu_new_2003.} The resulting model, referred to as RSV-AZ-ST, modifies the return innovation \( \epsilon_t \) and its correlation with \( \eta_t \) as follows:
where \( \lambda_t \sim \mathcal{IG}(\nu/2, \nu/2) \), \( \mu_\lambda = \nu / (\nu - 2) \), and \( \mathcal{TN}_{(0,\infty)} \) denotes a standard normal distribution truncated above zero. The constants \[ c = \mathrm{E}[z_{0t}] = \sqrt{ \frac{2}{\pi} }, \quad \mathrm{Var}[z_{0t}] = 1 - c^2 \] are used to ensure that \( \epsilon_t \) is standardized, so that its conditional variance remains \( \exp(h_t) \).
The contemporaneous correlation between \( \epsilon_t \) and \( \eta_t \) is now a function of \( \rho \), \( \nu \), and \( \delta \), and is given by
where \( \mathrm{E}[ \sqrt{ \lambda_t } ] \) is defined in equation (ref). Setting \( \delta = 0 \) recovers the RSV-T model. Moreover, as \( \nu \to \infty \) or \( \lambda_t = 1 \), the model simplifies to the RSV-AZ-SN model with a skew-normal distribution.
The shape of the distribution of \( \epsilon_t \) is governed jointly by \( \delta \) and \( \nu \), which respectively control skewness and tail thickness. From the general moment expression
we derive the third and fourth moments of \( \epsilon_t \) as:
As with the GH skew-\( t \) distribution, the skewness and kurtosis of the return distribution are jointly influenced by the parameters. While \( \nu \) governs tail behavior, \( \delta \) introduces asymmetry.
Figure (ref) illustrates the impact of these parameters on the shape of \( \epsilon_t \). Panel (i) shows that increasing negative values of \( \delta \) lead to stronger left-skewness, with \( \delta = 0 \) yielding a symmetric \( t \)-distribution. Panel (ii) demonstrates that larger values of \( \nu \) reduce tail thickness and approximate the normal distribution. Notably, even when \( \nu \to \infty \), the distribution remains skewed unless \( \delta = 0 \), in which case the skew-normal distribution is recovered.
Following the approach of fernandez_bayesian_1995, we define the Fern\'{a}ndez--Steel (FS) skew-\(t\) density as
where \( f_T(\cdot \mid \nu) \) denotes the probability density function of the standard Student’s \(t\) distribution:
Here, \( \mathbbm{1}\{ \cdot \} \) denotes the indicator function.
The mean and variance of \( w \sim p_T(w \mid \gamma, \nu) \) are:
where
To obtain a standardized version of the FS skew-\(t\) density, we define: \[ w_* = \frac{w - \mu_*}{\sigma_*}, \] so that \( \mathrm{E}[w_*] = 0 \) and \( \mathrm{Var}[w_*] = 1 \). The corresponding standardized density becomes:
Replacing \( f_T \) with the standard normal density \( f_N \) in the above expression yields the standardized FS skew-normal density \( q_N(\cdot \mid \gamma) \).
Following trottier_higher_2016, higher-order moments for the FS skew-$t$ distribution and their standardized forms are derived utilizing the absolute moments of the underlying Student's $t$ distribution. For the absolute moments, we refer to the unified results provided by kirkby_moments_2025. The $k$-th absolute moment, denoted by $M_k = 2 \int_0^\infty \tilde{w}^k f_T(\tilde{w} \mid \nu) d\tilde{w}$, is generally given by:
Specifically, for $k=3$ and $k=4$, we obtain:
Using the general raw-moment formula provided in trottier_higher_2016: \[ \mu_r' = E[w^r] = M_r \frac{\gamma^{r+1} + (-1)^r \gamma^{-(r+1)}}{\gamma + \gamma^{-1}}, \quad r=1,2,3,4, \] the centered third and fourth moments are defined as: \[ \mu_3 = \mu_3' - 3\mu_* \mu_2' + 2\mu_*^3, \qquad \mu_4 = \mu_4' - 4\mu_* \mu_3' + 6\mu_*^2 \mu_2' - 3\mu_*^4. \] Consequently, the skewness and excess kurtosis of the standardized variable $w_* = (w - \mu_*)/\sigma_*$ are given by:
These expressions implicitly depend on $(\gamma, \nu)$ through the moments $M_k$ and the parameter $\gamma$, confirming that $\gamma$ governs asymmetry while $\nu$ determines tail thickness.
We now define the RSV-FS-ST model by assuming:
where \( q_T(\cdot \mid \gamma, \nu) \) is the standardized FS skew-\(t\) density. Replacing \( q_T \) with \( q_N \) yields the RSV-FS-SN model.
The skewness of the distribution is governed by the parameter \( \gamma \), where \( \gamma = 1 \) corresponds to a symmetric \(t\)-distribution. If \( \gamma < 1 \), the distribution becomes left-skewed, while \( \gamma > 1 \) results in right-skewness. The parameter \( \nu \) controls tail heaviness.
Figure (ref) illustrates how the shape of the distribution varies with \( \gamma \) and \( \nu \). Panel (i) shows that decreasing \( \gamma \) induces stronger negative skewness, with \( \gamma = 1 \) yielding symmetry. Panel (ii) demonstrates that increasing \( \nu \) reduces tail thickness. Notably, as \( \nu \to \infty \), the distribution converges to a skew-normal form when \( \gamma \ne 1 \), and to a standard normal when \( \gamma = 1 \).
GARCH models are widely used to describe time-varying volatility. Among numerous variants, we consider an EGARCH-type model. While the original EGARCH specification of nelson_conditional_1991 involves the term $\mathrm{E}[|\epsilon_t|]$, this expectation does not admit a convenient closed-form expression under standardized $t$ innovations. Following hansen_exponential_2016, we therefore estimate the model using the quadratic term $(\epsilon_t^2-1)$ instead, which yields
with the initial condition \(h_1 = \omega\). The parameter \(\tau\) captures the asymmetric response of conditional volatility to return shocks, and when \(\tau < 0\), negative shocks increase volatility more than positive ones, reflecting the leverage effect. For expositional clarity, we present the model above under the Gaussian assumption for $\epsilon_t$. In Section (ref), we relax this assumption and consider a standardized Student-t distribution with degrees-of-freedom parameter $\nu$, allowing for heavy-tailed return behavior.
\paragraph{Remark.} nelson_conditional_1991 noted that estimating EGARCH with Student-t innovations can suffer from stability issues (p. 365), and suggested the generalized error distribution (GED) as a more robust alternative. While nelson_conditional_1991 originally proposed a quasi-maximum likelihood (QML) estimation procedure, in this study we estimate the model parameters using a Bayesian approach. We acknowledge the importance of the EGARCH-GED specification; however, given the additional computational complexity of Bayesian estimation and the main goal of this study to compare realized volatility-based stochastic volatility models, we leave the estimation of EGARCH-GED for future work.
hansen_realized_2012 and hansen_exponential_2016 extend the conventional GARCH and EGARCH models to the RGARCH and REGARCH frameworks, respectively, by incorporating realized measures into the volatility dynamics. Since the REGARCH model allows a more flexible characterization of the joint dynamics between returns and volatility than the RGARCH model, we adopt the REGARCH specification, which can be written as
with the initial condition \(h_1 = \omega\). Similarly to the EGARCH model, in Section (ref) we also consider a standardized Student-t distribution for $\epsilon_t$ in the REGARCH model, with degrees-of-freedom parameter $\nu$, allowing for heavy-tailed return behavior.
The latent log-volatility process \(h_t\) and the realized measure \(x_t\) each include a quadratic form of the leverage function, given by \(\tau_1 \epsilon_t + \tau_2 (\epsilon_t^2 - 1)\) and \(\delta_1 \epsilon_t + \delta_2 (\epsilon_t^2 - 1)\), respectively. Following hansen_exponential_2016, one can introduce an additional scaling parameter \(\psi\) as
but since the estimate of \(\psi\) is typically close to unity, we fix \(\psi=1\) for brevity. The parameter \(\zeta\) captures the bias in \(\log RV_t\); \(\zeta=0\) implies an unbiased realized measure.
Each RSV model introduced in Section (ref) can be formulated as a nonlinear Gaussian state space model with a large number of latent variables. This structure renders the likelihood function intractable in closed form, and consequently makes maximum likelihood estimation computationally difficult. To overcome this, we adopt a Bayesian estimation framework via MCMC simulation, following established approaches in the literature.
We assign prior distributions to the common model parameters \( (\mu, \phi, \sigma_\eta^2, \rho, \xi, \sigma_u^2) \) as follows:
where \( \mathcal{B}(\cdot, \cdot) \) and \( \mathcal{IG}(\cdot, \cdot) \) denote the beta and inverse-gamma distributions, respectively. The hyperparameters \( (m_\mu, s_\mu, a_{\phi 0}, b_{\phi 0}, a_{\rho 0}, b_{\rho 0}, n_\eta, S_\eta, m_\xi, s_\xi, n_u, S_u) \) are chosen to reflect prior beliefs or to be weakly informative when such prior knowledge is unavailable.
For notational convenience, let \( \bm{y} = (y_1, \ldots, y_n)' \), \( \bm{x} = (x_1, \ldots, x_n)' \), and \( \bm{h} = (h_1, \ldots, h_n)' \). Denote the parameter vector by \( \bm{\theta} \), and let \( \pi(\bm{\theta}) \) be its joint prior density. The joint likelihood of \( (\bm{y}, \bm{x}, \bm{h}) \) given \( \bm{\theta} \) is denoted by \( f(\bm{y}, \bm{x}, \bm{h} \mid \bm{\theta}) \).
In the following, we detail the MCMC sampling algorithms specifically for the RSV-AZ-ST and RSV-FS-ST models.\footnote{See takahashi_estimating_2009 and takahashi_volatility_2016 for the algorithms used in the RSV-N and RSV-GH-ST models, respectively.} We also outline the procedure for computing one-step-ahead forecasts of volatility and returns, from which we derive one-step-ahead VaR and ES.
Let \( \bm{\theta} = (\mu, \phi, \rho, \sigma_\eta^2, \delta, \xi, \sigma_u^2, \nu)' \) denote the full vector of model parameters, and define the latent variable vectors \( \bm{\lambda} = (\lambda_1, \ldots, \lambda_n)' \) and \( \bm{z}_0 = (z_{01}, \ldots, z_{0n})' \). For the prior distributions of the model-specific parameters \( \delta \) and \( \nu \), we assume:
where \( \mathcal{G}(\cdot, \cdot) \) denotes the gamma distribution.
The MCMC algorithm proceeds as follows:
Details of each sampling step, including full conditional distributions and implementation strategies (e.g., Gibbs sampling or Metropolis--Hastings updates), are provided in Appendix (ref).
Let \( \bm{\theta} = (\mu, \phi, \rho, \sigma_\eta^2, \gamma, \xi, \sigma_u^2, \nu)' \) denote the full set of model parameters. For the model-specific parameters \( \gamma \) and \( \nu \), we assign the following prior distributions:
The MCMC algorithm proceeds as follows:
The full conditional distributions and implementation details for each step are provided in Appendix (ref).
To obtain one-day-ahead forecasts of financial returns and volatilities, we utilize the predictive distribution within each state-space model. Let $\bm{\theta}^{(i)}$ and $\bm{h}^{(i)}$ represent the $i$th sample of parameters and latent log-volatilities in the MCMC simulation, respectively. The one-step-ahead predictive samples are then generated as follows:
Repetition of the above procedure for $M$ times allows us to obtain the generated samples $\{h_{n+1}^{(i)}\}_{i=1}^{M}$ and $\{y_{n+1}^{(i)}\}_{i=1}^{M}$.
One-day-ahead quantile forecasts, such as VaR and ES, can be derived from the distribution of $\{y_{n+1}^{(i)}\}_{i=1}^{M}$ for each model. The one-day-ahead VaR forecast at time $n+1$ at level $\alpha$, denoted by $\mbox{VaR}_{n+1}(\alpha)$, is defined as
where $\mathcal{I}_{n}$ denotes the information set available at time $n$. The corresponding ES forecast, denoted by $\mbox{ES}_{n+1}(\alpha)$, is given as
The one-day-ahead VaR and ES can be obtained as the $(1-\alpha)$th quantile and conditional average of $\{y_{n+1}^{(i)}\}_{i=1}^{M}$, respectively.
We evaluate the obtained forecasts for volatility, VaR, and ES using appropriate loss or scoring functions. To evaluate the volatility forecasts, we compute the Gaussian quasi-likelihood (QLIKE) loss function, defined as:
where $x$ and $f$ represent a volatility proxy and a volatility forecast, respectively. The QLIKE loss function is robust, as defined by patton_volatility_2011, meaning that it provides consistent rankings regardless of whether the ranking is based on true volatility or a conditionally unbiased volatility proxy. Furthermore, patton_evaluating_2009 demonstrated its superior power compared to the mean squared error (MSE), another robust loss function, in the predictive accuracy test proposed by diebold_comparing_1995.
For evaluating the forecasts of VaR and ES, we employ the joint loss function introduced by fissler_higher_2016 (FZ loss). Following the approach of patton_dynamic_2019, we utilize a specific form of the FZ loss function, referred to as the FZ0 loss function, which is expressed as:
where $y$, $v$, and $e$ represent a return, VaR, and ES, respectively. patton_dynamic_2019 demonstrated that the FZ0 loss function is unique in producing loss differences that are homogeneous of degree zero.
While average losses derived from the aforementioned loss functions offer initial insights into the forecast performance of the models in contention, they do not indicate if the differences in losses are statistically significant. To ascertain this, we utilize the conditional and unconditional predictive ability tests by giacomini_tests_2006, hereinafter referred to as GW tests, on the loss differences. The GW tests are particularly pertinent to this paper's objectives. This is because they accommodate the rolling window methods adopted in Section (ref), enabling a unified assessment of both nested and non-nested models, including the RSV models described in Section (ref).
For the unconditional predictive ability, the GW test statistic conforms to that proposed by diebold_comparing_1995 and is, under the null hypothesis, asymptotically standard normally distributed. For the conditional predictive ability, the GW test defines the null hypothesis as
where $\Delta L_{t}$ is the loss difference between two models at forecast date $t$, and $\mathcal{I}_{t}$ is the information set at time $t$.
Given a test function $\mathbbm{h}_{t}$, a $q \times 1$ vector measurable with respect to $\mathcal{I}_{t}$ stinchcombe_consistent_1998, the null hypothesis translates to
Defining $Z_{t+1} = \mathbbm{h}_{t} \Delta L_{t+1}$ and $\bar{Z}_{n_{f}} = n_{f}^{-1} \sum_{t=n}^{n+n_{f}-1} Z_{t+1}$, the alternative hypothesis becomes
We employ the Wald-type test statistic
where $\hat{\Omega}_{n_{f}} = n_{f}^{-1} \sum_{t=n}^{n+n_{f}-1} Z_{t+1} Z_{t+1}'$ consistently estimates the variance of $Z_{t+1}$. Under the null, $T_{n_{f}}$ follows a chi-square distribution with $q$ degrees of freedom, $\chi_{q}^{2}$. At the significance level $p$, the null hypothesis is rejected if $T_{n_{f}} > \chi_{q,1-p}^{2}$, where $\chi_{q,1-p}^{2}$ is the $(1-p)$ quantile of $\chi_{q}^{2}$.
In practice, the selection of $\mathbbm{h}_{t}$ should discern between the forecast performances of the models. Focusing on one-day-ahead forecasts, we follow the specifications of takahashi_stochastic_2023 and define $\mathbbm{h}_{t} = (1, \Delta L_{t})'$ for $t = n, n+1, \ldots, n+n_{f}-1$. If the null is rejected, it implies that the lagged loss differences, $\Delta L_{t}$, can help predict the subsequent loss differences, $\Delta L_{t+1}$. Using $\hat{b}$ to represent the regression coefficient of $\Delta L_{t+1}$ on $\mathbbm{h}_{t}$, the predicted loss differences $\{\hat{b}' \mathbbm{h}_{t}\}_{t=n}^{n+n_{f}-1}$ serve as indicators for assessing model performances across different time points. In Section (ref), we measure relative performance by the frequency with which one model forecasts larger losses than the other, as expressed by
Following mitsui_bayesian_2003, we use the Metropolis--Hastings (MH) algorithm to sample the parameters of the EGARCH and REGARCH models from the posterior distribution.\footnote{For other Bayesian approaches to RGARCH-type models and their applications to volatility and tail risk forecasting, see chen_bayesian_2023.} First, we obtain the mode of the posterior density, and then set the proposal density to a multivariate normal distribution centered at the posterior mode, with a covariance matrix proportional to the inverse Hessian of the log posterior density evaluated at the mode, scaled by a factor of 1.2. The scaling factor is introduced to ensure sufficient exploration of the tails of the posterior distribution. asai_comparison_2006 shows that this method is efficient for sampling the parameters of GARCH-type models from the posterior distribution.
The prior distributions for the EGARCH model are specified as follows:
For the additional parameters of the REGARCH model, the prior distributions are given by \[ -\tau_1 \sim \mathcal{G}(0.1,0.1), \quad \tau_2 \sim \mathcal{G}(0.1,0.1), \quad \zeta \sim \mathcal{N}(0,10), \quad -\delta_1 \sim \mathcal{G}(0.1,0.1), \quad \delta_2 \sim \mathcal{G}(0.1,0.1), \quad \sigma_\upsilon^2 \sim \mathcal{IG}(2.5,0.5). \] The parameters $\tau$, $\tau_1$, $\tau_2$, $\gamma$, $\delta_1$, and $\delta_2$ are, in principle, unrestricted in sign. However, empirical evidence typically suggests $\tau<0$ for EGARCH, and $\tau_1<0$, $\tau_2>0$, $\gamma>0$, $\delta_1<0$, and $\delta_2>0$ for REGARCH models. Accordingly, we employ Gamma priors with appropriate sign restrictions to confine these parameters to economically meaningful regions.
Based on the posterior draws, one-day-ahead forecasts of volatility, VaR, and ES are obtained. This approach allows us to jointly sample both parameters and the implied conditional volatilities from the posterior distribution, thereby explicitly incorporating parameter uncertainty into volatility and return forecasts.
Importantly, this feature addresses concerns raised in the literature that ignoring parameter uncertainty in GARCH-type models may lead to unfair comparisons with stochastic volatility models; see, for example, ardia_forecasting_2018, who document that accounting for parameter uncertainty improves volatility and risk forecasts in GARCH frameworks. By adopting a fully Bayesian estimation strategy, our framework ensures that parameter uncertainty is treated consistently across models, providing a fair basis for comparison and yielding coherent predictive distributions for one-step-ahead volatility and return forecasts.
We estimate the SV and RSV models described in Section (ref) using daily (close-to-close) returns and RVs for two major stock indices: the DJIA and the N225. Following liu_does_2015, we adopt the 5-minute RV estimator, computed during trading hours only, among several available alternatives. The DJIA data are obtained from the Oxford-Man Institute’s Realized Library,\footnote{The website \url{https://realized.oxford-man.ox.ac.uk/} is no longer accessible.} while the N225 data are constructed from the Nikkei NEEDS-TICK dataset.\footnote{See ubukata_pricing_2014 for details on the construction of the N225 dataset.}
The sample period spans from June 2009 to September 2019 for both indices. Specifically, the DJIA sample covers 2,596 trading days from June 1, 2009, to September 27, 2019, and the N225 sample includes 2,532 trading days from June 1, 2009, to September 30, 2019.
Figure (ref) presents time series plots and histograms of daily returns for the DJIA and N225. Both return series exhibit considerable variation around zero and show frequent large negative returns, indicating negatively skewed distributions.
Table (ref) reports the descriptive statistics for the daily returns. The mean return is significantly different from zero for the DJIA, but not for the N225. However, since the means are negligible in magnitude, we do not demean the series in the subsequent analysis. The $p$-values from the Ljung--Box statistic ljung_measure_1978, adjusted for heteroskedasticity as in diebold_empirical_1988, do not reject the null of no autocorrelation up to 10 lags in either series. This allows us to estimate the models directly using raw returns.
Both return series are characterized by significantly negative skewness and high kurtosis, suggesting leptokurtic distributions---a stylized fact in financial returns. The Jarque--Bera (JB) statistic confirms that normality is strongly rejected for both series. These empirical features motivate our use of skewed-$t$ distributions for the return innovations, as discussed in Section (ref).
Figure (ref) shows the time series and histograms of the 5-minute RVs and their logarithmic transformations for both indices. Both series exhibit strong temporal clustering, high persistence, and occasional sharp spikes.\footnote{For the DJIA, the largest RV spike occurred on August 24, 2015, amid heightened market turbulence (see The New York Times: \url{https://www.nytimes.com/2015/08/25/business/dealbook/daily-stock-market-activity.html}). For the N225, the most notable spike on March 15, 2011, reflects the aftermath of the Great East Japan Earthquake on March 11, 2011.} These spikes lead to positively skewed distributions in the log-RVs.
Table (ref) summarizes the descriptive statistics of the log-RVs. The Ljung--Box test strongly rejects the null of no autocorrelation, confirming the presence of volatility clustering in both markets. The distributions are positively skewed and leptokurtic, and the JB test rejects the normality of the log-RVs. This finding is inconsistent with the normality assumption for the error term \(u_t\) in equation (ref), which we retain for tractability.
\paragraph{Remark.} The positive skewness of the log-RVs suggests that more flexible innovations could be considered. In particular, specifying \(u_t\) as skew-t would allow the model to directly accommodate distributional asymmetry. Incorporating such an extension into the RSV framework would require additional development, so we leave this promising direction for future work.
We evaluate the forecasting performance of the following RSV models, which differ in the distributional assumption for $\epsilon_t$:
The prior distributions for the common parameters are specified as follows:
These priors are chosen to be weakly informative and follow standard specifications commonly adopted in the stochastic volatility literature, allowing the data to dominate posterior inference.
For the additional distributional parameters, we set:
The truncation $\nu>4$ ensures the existence of the fourth moment, which is required for volatility forecasting and risk evaluation.
In addition, to isolate the contribution of RV, we estimate standard SV models with the same innovation distributions as in the RSV specifications. These include SV-N, SV-T, SV-GH-ST, SV-AZ-SN, SV-AZ-ST, SV-FS-SN, and SV-FS-ST, mirroring the distributional assumptions used in the RSV models. This setup allows us to directly assess the added value of incorporating RV in terms of predictive performance.
To further benchmark the RSV models, we also consider the EGARCH and REGARCH models described in Sections (ref) and (ref), respectively. These models are estimated within a Bayesian framework using the MH algorithm under standard normal and Student’s $t$ innovations. One-day-ahead forecasts of volatility, VaR, and ES are then obtained based on the posterior draws. For other Bayesian approaches to RGARCH-type models and their applications to volatility and tail risk forecasting, see chen_bayesian_2023.
Following the approach of takahashi_volatility_2016, we implement a rolling window estimation procedure to evaluate the out-of-sample performance of each model. The window size is kept fixed throughout the forecasting period. For the DJIA, we use a window of 1,993 observations, generating forecasts from May 1, 2017, to September 27, 2019. For the N225, the window consists of 1,942 observations, with forecast dates ranging from May 1, 2017, to September 30, 2019. After each estimation step, we produce one-day-ahead forecasts of volatility, VaR, and ES.
At each forecast point, we draw 15,000 predictive samples from the posterior predictive distribution. For the SV and RSV models, we compute both the posterior means and medians of the volatility forecasts.\footnote{We report the median rather than the mean for volatility forecasts of SV models, as the posterior predictive distribution of the one-step-ahead volatility forecast occasionally exhibits heavy tails, which can disproportionately inflate the posterior mean.} Additionally, we calculate predictive quantiles of the return distribution to obtain VaR and ES estimates. This procedure yields 603 forecasts for the DJIA and 590 for the N225, covering the period from early May 2017 to late September 2019.
To evaluate the accuracy of volatility forecasts, we compute the average loss using the QLIKE loss function, as described in Section (ref). The QLIKE loss is known for its robustness, yielding model rankings that align closely with those based on latent volatility, provided the proxy is conditionally unbiased. Simulation evidence from patton_evaluating_2009 and empirical findings by hansen_forecast_2005 and patton_volatility_2011 suggest that QLIKE has greater statistical power than MSE to discriminate between models.
To mitigate market microstructure noise, we employ multiple volatility proxies: realized kernel (RK) with a flat-top Tukey-Hanning${}_2$ kernel barndorff-nielsen_designing_2008, bipower variation (BV) barndorff-nielsen_power_2004, and median realized volatility (Med) andersen_jump_2012, in addition to the standard 5-minute realized volatility (RV5).
To correct for biases due to non-trading hours, we adopt the adjustment method proposed by hansen_forecast_2005, as given in equation (ref), using the following correction factor:
where $n$ is the window size and $x_s$ denotes the volatility proxy at time $s$.
Figures (ref) and (ref) display the log-scale volatility forecasts alongside the adjusted RV5 series for the DJIA and N225. For the DJIA, SV model forecasts---especially those from the SV-GH-ST specification---exhibit pronounced volatility spikes. Across both indices, forecast patterns are more homogeneous within model families (e.g., SV, RSV, EGARCH, REGARCH) than across different distributional assumptions, indicating that model structure plays a more dominant role than distributional form. Figures with alternative realized volatility proxies are omitted for brevity, as the volatility forecasts are identical.
Table (ref) reports QLIKE scores across all models and volatility proxies. An asterisk ($^*$) denotes inclusion in the 90% model confidence set (MCS) of hansen_model_2011.\footnote{The rolling window scheme satisfies the stationarity assumption required by the MCS bootstrap. We use 1,000 bootstrap replications with a block size of 10. See hansen_model_2011, Section 4.3.} Overall, RSV and REGARCH models attain lower QLIKE scores than SV and EGARCH counterparts. For the DJIA, REGARCH-T achieves the lowest scores under RV5 and BV, while SV-FS-ST achieves the lowest scores under RK and Med. For the N225, RSV models dominate REGARCH models. {On average, the RSV-GH-ST model delivers the best performance across indices, followed by the RSV-AZ-ST and REGARCH-T models, with very small margins.}
To assess time-varying performance, Figures (ref) and (ref) depict cumulative loss differences (CLDs) relative to SV-N, based on RV5 as the volatility proxy.\footnote{Results based on alternative realized volatility proxies are reported in the Supplementary Material.} Positive values indicate superior performance. The RSV and REGARCH models consistently outperform the SV and EGARCH models across the entire forecast horizon for both the DJIA and N225, underscoring the value of incorporating realized volatility in prediction. In particular, the performance advantage of RSV and REGARCH is stable over time, with clear and persistent separations from the SV and EGARCH benchmarks in both markets.
Figures (ref) and (ref) visualize the GW test results, based on RV5 as the volatility proxy, using heatmaps.\footnote{Results based on alternative realized volatility proxies are reported in the Supplementary Material.} The lower triangular part presents unconditional GW test statistics, whereas the upper triangular part reports win proportions as defined in Equation (ref). Lower values (green) favor the row model, while larger values (orange) indicate better performance by the column model. Double and single asterisks ($**$ and $*$) denote significance at the 1% and 5% levels, respectively, based on unconditional and conditional GW test p-values for the lower and upper triangular parts.
For the DJIA, SV-FS-SN and SV-FS-ST outperform other SV variants and EGARCH models, while RSV and REGARCH models generally dominate SV and EGARCH. RSV-AZ-ST is significantly superior to even FS-type SV models, although RSV and REGARCH differ less markedly. Conditional GW tests further confirm the superiority of RSV models, particularly over SV and EGARCH.
For the N225, both unconditional and conditional GW tests support the outperformance of RSV and REGARCH over SV and EGARCH, with RSV models often exceeding REGARCH. FS-type distributions enhance forecast accuracy within the SV and RSV frameworks. Among RSV specifications, RSV-GH-ST lags significantly.
In summary, incorporating RV substantially improves forecast accuracy, with RSV and REGARCH models consistently outperforming SV and EGARCH counterparts. These results corroborate prior findings takahashi_stochastic_2023, takahashi_forecasting_2024. The added flexibility of FS-type distributions further improves forecast accuracy, especially in SV models. The next section examines whether these improvements carry over to tail risk measures, including VaR and ES.
Figure (ref) illustrates the 1% VaR and ES forecasts for the DJIA and N225 indices, as generated by the RSV-AZ-ST and REGARCH-T models. VaR violations---instances where realized returns fall below the forecasted VaR---are observed across both models, occurring during periods of market stress such as early 2018 and late 2019, as well as during more tranquil phases. For clarity, the forecasts from other models, which follow similar patterns with minor differences in scale, are omitted.
Table (ref) presents the empirical violation rates ($\hat{\alpha}$), the average FZ0 loss values (as defined in Equation (ref)), and the corresponding $p$-values from the dynamic quantile (DQ) test of engle_caviar_2004 and the MCS procedure for the DJIA. At the 1% level, all models tend to underestimate tail risk, as indicated by violation rates exceeding the nominal level. In contrast, for $\alpha = 5\%$, the violation rates are generally well aligned with the target.
Regarding predictive accuracy based on the FZ0 loss, the SV-AZ-SN and RSV-AZ-ST models achieve the lowest losses at $\alpha = 1\%$ and 5%, respectively. Most models---except SV-GH-ST at $\alpha = 5\%$---are included in the 75% MCS, indicating comparable performance in tail risk forecasting.
The DQ test, which regresses VaR violations on a constant, the forecasted VaR, and a one-period lagged violation indicator, suggests that the SV-FS-ST and SV-GH-ST models fail the independence test at the 1% level. This implies potential misspecification in their dynamic structures.
Table (ref) summarizes the results for the N225. Across both target levels, most models yield violation rates consistent with the nominal levels. The RSV-AZ-ST model again performs best in terms of FZ0 loss, and both RSV and REGARCH models are included in the 75% MCS. Notably, the SV-GH-ST and EGARCH-N models fall outside the MCS at the 5% level, while the EGARCH-T model remains outside the MCS even at the 1% level, indicating weaker predictive accuracy. Unlike the DJIA case, none of the models are rejected by the DQ test for the N225.
Figures (ref) and (ref) show the CLDs relative to the SV-N benchmark for the DJIA. At the 1% level, REGARCH models perform poorly during the early 2018 volatility spike. From mid-2018 onwards, however, REGARCH and RSV models consistently outperform SV and EGARCH models at the 5% level, highlighting the benefits of incorporating realized volatility into tail risk forecasting. Consistent results are obtained for the N225, reinforcing the benefit of incorporating RV, and are therefore reported in the Supplementary Material.
\paragraph{GW test results}
Figures (ref) and (ref) visualize the GW test results based on the FZ0 loss function for the DJIA using heatmaps. The interpretation of the lower and upper triangular elements follows that of volatility forecasts in Section (ref). The main findings are summarized below.
At the 1% level, the unconditional GW test reveals that SV-AZ-SN significantly outperforms SV-N, while SV-GH-ST is significantly outperformed by SV-AZ-SN. Among RSV models, RSV-N exhibits inferior predictive accuracy relative to all other RSV specifications. For EGARCH models, EGARCH-T significantly outperforms EGARCH-N, while REGARCH models are significantly outperformed by several RSV specifications.
The conditional GW test at the 1% level indicates that SV-T performs significantly worse than SV-AZ-SN, but better than EGARCH-N. RSV-N is significantly outperformed by RSV-T, RSV-AZ-SN, RSV-FS-SN, and RSV-FS-ST. Moreover, RSV models with skew(-t) distributions significantly outperform REGARCH-N and REGARCH-T.
At the 5% level, the unconditional test shows that SV-GH-ST is significantly outperformed by SV-T, SV-AZ-SN, SV-AZ-ST, and all RSV and REGARCH models. In addition, RSV-GH-ST outperforms SV-N, and RSV-AZ-ST significantly outperforms RSV-AZ-SN. Among EGARCH and REGARCH models, the T specifications consistently outperform their N counterparts.
For the conditional GW test at the 5% level, SV-GH-ST is significantly outperformed by RSV-AZ-SN, RSV-AZ-ST, RSV-FS-SN, RSV-FS-ST, and RSV-GH-ST. RSV-T performs significantly worse than RSV-GH-ST. EGARCH-N is significantly outperformed by EGARCH-T, and REGARCH-N by REGARCH-T.
In summary, the results highlight the superior performance of RSV models---especially RSV-AZ-ST, RSV-FS-ST, and RSV-GH-ST---compared to SV, EGARCH, and REGARCH models, particularly under stringent risk levels. The inclusion of RV consistently improves the accuracy of tail risk forecasts. Moreover, as the FZ0 loss function jointly evaluates VaR and ES by emphasizing tail behavior, these results underscore the importance of capturing skewness and heavy tails in return distributions. The results for the N225 are broadly consistent with those for the DJIA at the 5% level and are therefore reported in the Supplementary Material.
\paragraph{Distributional characteristics of daily return forecasts under RSV models}
To explore differences among RSV models, we examine the distributional properties of their one-day-ahead return forecasts. Table (ref) reports the mean and standard deviation of skewness, kurtosis, and several lower percentiles (0.1%, 1%, 5%, and 10%) based on the posterior predictive distributions. Figures (ref) and (ref) present the histograms of these statistics across the forecast horizon for the DJIA and N225, respectively.
For both indices, the RSV-FS-SN and RSV-GH-ST models tend to generate more negatively skewed and heavy-tailed predictive distributions, as reflected in their higher kurtosis and lower percentile values. Among all models, RSV-GH-ST exhibits the most negative skewness, while RSV-FS-SN displays the highest kurtosis.
For the N225 index, the histograms of the 0.1st and 1st percentiles indicate that RSV-FS-SN consistently produces more conservative forecasts in the extreme left tail. Interestingly, its 10th percentile is relatively higher than those of other models. These distributional characteristics may help explain the variation in VaR and ES forecast performance observed across RSV models.
\paragraph{Remark.} Additional results for the COVID-19 period are reported in the Supplementary Material. While RSV models with skew-$t$ specifications improve forecast accuracy for volatility in the main sample, their relative performance for VaR and ES forecasts deteriorates during the COVID-19 period. This finding suggests that the benefits of modeling conditional asymmetry in the return distribution may be less pronounced under extreme market stress, where tail risk dynamics are dominated by abrupt and persistent shocks. A more detailed investigation of this issue is left for future research.
Table (ref) presents the rankings of average losses for volatility, VaR, and ES forecasts. The average ranking in the final column highlights the overall superiority of the RSV models, with RSV-AZ-ST and RSV-GH-ST achieving the best performance.
Our findings demonstrate that RV is a valuable predictor for forecasting not only volatility but also tail risk measures such as VaR and ES. Overall, RSV models tend to deliver superior performance, particularly for the N225, while REGARCH-T remains highly competitive for the DJIA. Incorporating a skewed $t$-distribution within the RSV framework further enhances forecast accuracy, especially for tail risk measures.
This study has evaluated the predictive performance of multiple models for forecasting volatility, VaR, and ES using data from two major financial indices: the DJIA and N225. Our comprehensive assessment, summarized in Table (ref), demonstrates that the RSV models consistently outperform alternative approaches across all forecast objectives.
Among the models analyzed, the RSV-AZ-ST specification achieves the best overall performance. This highlights not only the value of incorporating RV but also the benefit of accounting for skewness and heavy tails in return distributions. Specifically, the use of skewed $t$-distributions within the RSV framework substantially improves forecast accuracy for both volatility and tail risk measures.
A key contribution of this study is the joint validation of these modeling enhancements: integrating RV and adopting flexible return distributions. These improvements offer methodological guidance for constructing robust models for financial risk forecasting.
In sum, our findings suggest that RSV models---particularly those incorporating skewed $t$-distributions---constitute a robust and effective approach for forecasting volatility, VaR, and ES. These results emphasize the critical role of model specification and distributional assumptions in improving the reliability of financial risk forecasts.
An interesting avenue for future research is to extend the proposed framework to allow for alternative or multiple realized measures within the measurement equation. In the present study, the RSV specification relies on log-RV5 as a single realized proxy, while forecast evaluation is conducted using several realized measures, including RK, BV, and Med. Allowing the model to incorporate alternative realized measures, or a multi-measure RSV specification, could help assess the robustness of the reported gains with respect to the choice of realized volatility proxy and may also mitigate the distributional mismatch suggested by the non-normality of log-RV documented in the empirical analysis.
A related line of research is provided by score-driven volatility models with flexible innovation distributions, such as the GH skew-\(t\)-based SD--GARCH specification of catania_forecasting_2022, as well as joint SD models for returns and realized volatility proposed by catania_forecasting_2020. While these approaches require more elaborate dynamic updating rules and likelihood evaluations, they provide informative points of comparison and could offer further insights into the role of realized measures and distributional flexibility in volatility modeling.