Forecasting the CBOE VIX and SKEW Indices Using Heterogeneous Autoregressive Models

Abstract

We analyze the predictability of daily data on the CBOE $VIX$ and $SKEW$ indices, used to capture the average level of risk-neutral risk and downside risk, respectively, as implied by S&P 500 index options. In particular, we use forecast models from the Heterogeneous Autoregressive ($HAR$) class to test whether and how lagged values of the $VIX$ and of the $SKEW$ may increase the forecasting power of $HAR$ for the $SKEW$ and the $VIX$. We find that a simple $HAR$ is very hard to beat in out-of-sample experiments aimed at forecasting the $VIX$. In the case of the $SKEW$, the benchmarks (the random walk and an $AR\left(1\right)$) are clearly outperformed by $HAR$ models at all the forecast horizons considered and there is evidence that special definitions of the $SKEW$ index based on put options data only yield superior forecasts at all horizons.

1. Introduction

Volatility modeling and forecasting is an important issue of research in financial markets, and it has received much attention of researchers and practitioners over the last three decades. A range of different models have been used, starting from simple one-factor Autoregressive Conditional Heteroskedasticity ($ARCH$) type models, extending to multiple-factor stochastic volatility models, and eventually–owing to the increasing availability of high-frequency intraday data–to the use of realized variance for measuring and modeling volatility (see the review of the literature standard financial econometrics textbooks, e.g., Guidolin and Pedio [1], Hautsch [2]). In spite of such a rich drive to develop and test progressively more complex models, the literature has increasingly emphasized how the forecasts of future volatility of the underlying returns that (whenever possible) can be inferred from the price of traded options can and do help to sharpen the predictive accuracy of many models (see, e.g., Han and Park [3] and references therein). Therefore one portion of such an expanding body of research has investigated the dynamic properties and predictability of the main indices that have been proposed to measure implied volatility and its main features, such as the Chicago Board Options Exchange (CBOE) Volatility Index (in what follows, referred to simply as the $VIX$) and the CBOE Skew Index (the $SKEW$) in the case of S&P 500 index options (see Fassas and Siriopoulos [4] for a review of implied volatility indices).

One of the most important characteristics of both realized and option-implied variance is its long memory feature, i.e., its strong persistence of their sample autocorrelations. As argued by Bandi and Perron [5], this may be an indication of the presence of fractional integration in both series. Fractional integration implies that the autocorrelations of a time series decay at a hyperbolic rate rather than at an exponential rate, as seen in non-fractionally integrated series. This slower decay indicates that even distant past values have a significant impact on the current and future states of the series. As it is well known, in the case of physical realized variance, to account for this stylized fact, Corsi [6] developed a simple and parsimonious Heterogeneous Autoregressive ($HAR$) model that builds on the Heterogeneous Market Hypothesis. With reference to spot, physical realized variance, Corsi showed that the model provided more accurate out-of-sample forecasts of realized variance than the alternatives. Even though $HAR$ has turned out to be considerably successful, both in its original structure and in its subsequent extensions (see, e.g., the review in Corsi et al. [7]), much less is known about the performance of $HAR$ models when applied to the classical, key implicit volatility indices, such as the $VIX$ and especially the $SKEW$. The goal of our paper is exactly to provide a thorough, state-of-the-art assessment of the performance of the baseline $HAR$ model and a few of its variations in the prediction of the CBOE Volatility Index and the Skew Index. The choice to stick to rather simple forms of $HAR$ modelling derives from the fact that it is generally challenging to provide more complex specifications of the pure $HAR$ model that clearly outperform it in terms of forecasting accuracy. The model presented by Psaradellis and Sermpinis [8], the $HAR-GASVR$ model, incorporates a Support Vector Regression ($SVR$) algorithm optimized using a Genetic Algorithm ($GA$) to capture nonlinear relationships in the data. Among all tested hybrid models, the $HAR$-$GASVR$ stands out as the only one achieving significantly superior forecasting performance compared to the pure $HAR$ model. In particular, and quite distinctively, we rely on the $HAR$ model not only for modeling and forecasting the daily behaviour of the $VIX$ but also for analysing the time series of the $SKEW$ index. Since both indices display similar properties in their time series, we believe the same rationale that led to the development of the $HAR$ model for realized volatilities (and its extension to implied volatilities) can be extended to risk-neutral skewness, and thus to the $SKEW$ index.

The Chicago Board Options Exchange Volatility Index is a real-time market index representing the market’s expectations for volatility over the coming 30 days. Often referred to as the “fear gauge”, the $VIX$ measures anticipated market volatility derived from the prices of S&P 500 index options. The focus is primarily on Out-of-The-Money (OTM) options, which are considered more sensitive to market volatility expectations. The $VIX$ was first introduced in 1993 by the CBOE and it was initially based on the implied volatility of eight S&P 100 (OEX) options. However, in 2003, the CBOE updated the $VIX$’s calculation methodology to better reflect market conditions and expectations by shifting the underlying index to the S&P 500. This change allowed for a broader and more comprehensive measure of market volatility, as the S&P 500 is a more widely followed index with a larger number of component stocks. The $VIX$ calculation is based on the prices of a wide range of S&P 500 index options, particularly focusing on options (both puts and calls) with 23 to 37 days until expiration. It uses a model-free approach that does not assume any specific distribution for the returns of the underlying assets, see, e.g., Whaley [9] for a general discussion. The $VIX$ serves as a crucial tool for investors and analysts, providing insights into market sentiment and expectations of volatility. It is widely used for hedging purposes and as a benchmark for volatility-based derivatives (see, e.g., Zhang et al. [10]).

The CBOE Skew Index was designed (and officially launched in 1990) to measure the market’s perception of tail risk, which refers to the risk of extreme events that are unlikely but potentially very impactful on aggregate stock prices. The concept of “skew” in this context refers to the tendency of OTM put options to be more expensive than OTM call options, reflecting the market’s greater concern for significant downside risks than upside potential. Specifically, the Skew Index is calculated using the prices of out-of-the-money options on the S&P 500 index and it quantifies the market’s expectations of a significant deviation from the average, particularly focusing on the risk of extreme negative returns. (The index values typically hover around 100 in stable market conditions. Values much higher than 100 suggest that market participants are expecting higher chances of large drops in the index than large rises, indicating greater perceived tail risk). The $SKEW$ is considered a “fear index” similar to the $VIX$, but it specifically measures the perceived risk of significant negative market moves. A higher Skew indicates that investors are willing to pay more for downside protection, suggesting concerns about potential large market downturns. Therefore the $SKEW$ can provide valuable insights into market sentiment and the likelihood of tail events, which are important for risk management and investment decision-making: it serves as a signal to investors about the prevailing market fears and can indicate periods of higher financial market vulnerability, see also the discussion in Cao et al. [11].

The literature regarding the use of the $VIX$ to forecast a range of key financial market variables is massive and we shall make no effort to summarize it here because of space constraints. There is however also a literature on forecasting (see also the discussion in Majmudar and Banerjee [12]), the $VIX$ itself that was started with Blair et al. [13], who had emphasized the $VIX$’s predictability, comparing it with historical volatility measures, and hence its superior performance as a predictor of future volatility. Degiannakis [14] was the first known paper to apply Autoregressive Integrated Moving Average ($ARIMA$) models to predict the $VIX$, demonstrating that an $ARIMA(1,1,1)$ model could effectively forecast the $VIX$’s direction with over 58 percent accuracy, while adding other financial or macroeconomic variables would not significantly improve the model’s performance. Ahoniemi [15] generalized this effort by examining various time-series models for predicting the $VIX$, including $ARIMA$ and a Generalized Autoregressive Conditional Heteroskedasticity ($GARCH$) specifications. The study found that these models, particularly when combined with financial indicators like S&P 500 returns, can predict the direction of the $VIX$ with notable accuracy. Fernandes et al. [16] have evaluated different $HAR$-type models, showing that $HAR$ consistently outperforms simpler models like the autoregressive model and the random walk in terms of predictive power for the $VIX$, especially in short-term forecasts. Recently, Taylor [17] has proposed a model that includes various economic indicators and financial variables and applied multivariate regressions and machine learning methods, to analyze and predict $VIX$ futures returns. He reports that his model significantly outperforms traditional linear models in both statistical accuracy and economic value terms.

There is also a small literature that has focused on predicting the CBOE Skew Index itself. Mora-Valencia et al. [18] have used $ARIMA$ $GARCH$ models to forecast the CBOE Skew Index and found that a Moving Average ($MA$) model and an Exponential Generalized Autoregressive Conditional Heteroskedasticity ($EGARCH$) model, specifically a $MA\left(4\right)-EGARCH(1,3)$, with Student’s t-distributed innovations was particularly effective in capturing the characteristics of the Skew Index. They found that the $VIX$, alongside other variables, could help forecast the $SKEW$, indicating that changes in market volatility, as measured by the $VIX$, are related to expectations of tail risk reflected in the Skew Index. They also report that accurate forecasting of the Skew Index could provide warnings about potential market downturns. Bevilacqua and Tunaru [19] have examined the predictive power of decomposed components of the $SKEW$, the Positive Skew Index ($SKE{W}^{+}$) and the Negative Skew Index ($SKE{W}^{-}$), for forecasting the index itself. They found that including these components in a First-Order Autoregressive ($AR\left(1\right)$) model improved out-of-sample forecasting performance, suggesting that the decomposition provides more nuanced predictive information about market sentiments and tail risks. Moreover, including the lagged value of the $VIX$ in their models improved the out-of-sample forecasting accuracy for the $SKEW$ Index. This suggests that the $VIX$, which captures market volatility expectations, provides significant information that can be used to predict changes in the $SKEW$ (Elyasiani et al. [20] have explored the $SKEW$ as a measure of market greed and fear, indicating its usefulness in predicting market sentiment shifts).

We use daily data on the CBOE $VIX$, the $SKEW$, and the panel of all put and call options that were tradable on the S&P 500 index for the sample period from 4 January 1996, to 31 December 2019, to perform a classical recursive Out-Of-Sample (OOS) forecasting exercise. We limit our data to 2019 to avoid dealing with the disruptions in the role and informational effectiveness of the US index options market, from which the $VIX$ and $SKEW$ are derived, see, e.g., John and Li [21]. Of course, it would be interesting to extend our analysis and formally contrast the pre- and post-2020 samples, when enough data will become available. Where appropriate, besides data from the CBOE website, to decompose the $SKEW$ index and extract only its negative part or distill a version of the index only based on put options, we compute afresh the entire daily history of the index. In the recursive, prediction exercise, we employ a 1000-day rolling window recursive forecasting scheme. On each day starting with observation 1001, we compute the 1-, 5-, 10-, and 21-day ahead forecasts. To produce predictions beyond the 1-day horizon, we employ a classical, indirect dynamic forecast scheme.

We obtain three key results. First, within an otherwise standard $HAR$ framework, the $SKEW$ index does not contain useful information to improve the in-sample prediction of the $VIX$ index. However, when it comes to perform in-sample prediction of the $SKEW$ (especially of its logarithm), both the lagged value of $SKE{W}^{-}$ and the lagged value of the $VIX$ improve the predictive power of the $HAR$ model. Second, a simple, Corsi [6]-style, simple $HAR$ is very hard to beat in OOS experiments aimed at forecasting the $VIX$. Specifically, using Hansen’s test for Superior Predictive Ability (SPA) test for the squared error loss function, neither the $HAR$-$SK$ nor the $HAR$-$SK$-$PUTS$ models turn out to be clearly superior in any of the cases investigated. The same conclusion holds under an absolute error loss function. Furthermore, when forecasting the levels of the $VIX$, the $HAR$ model is clearly superior to its two extensions, except for the 1-day ahead forecast horizon. In fact, when the MCS methodology is deployed, the naive $HAR$ model turns out to be the only model included in the confidence set 85% of time. Third, when it comes to forecast the $SKEW$, the benchmarks (here, the random walk and an $AR\left(1\right)$) are clearly outperformed by $HAR$ models in predicting all the transformations of the $SKEW$ used as a target variable and at all the forecast horizons considered. The benchmarks generally fail to be included in the MCS, in fact. Fourth, again in the case of the $SKEW$, when we compare the performance of the $HAR$ model with the $HAR$-$SK$-$PUTS$ model, the latter produces superior forecasts at all horizons and this tends to be confirmed by MCS-based evidence.

The rest of the paper is organized as follows. Section 2 adds to the paper already surveyed here and performs an additional, short review of the literature. Section 3 introduces the models tested and the OOS research design adopted in our paper. Section 4 describes the data and presents preliminary summary statistics. Section 5 shows our main empirical results, especially those concerning the OOS performance of the different models. Section 7 concludes.

2. Brief Literature Review

As already mentioned, there exists a relatively large body of literature that has employed econometric methods to forecast the $VIX$ index. For instance, Ahoniemi [15] estimated various time series models for the $VIX$ using daily data from 1990 to 2007. The best model in terms of in-sample goodness-of-fit and OOS correct sign prediction turned out to be a simple $ARIMA(1,1,1)-GARCH(1,1)$. However, the best performance in terms of the Mean Squared Forecast Error ($MSFE$) criterion was given by homoskedastic $ARIMA(1,1,1)$ and a Autoregressive Integrated Moving Average with Exogenous variables $ARIMAX(1,1,1)$. Konstantinidi et al. [22] expanded this effort by studying the predictability of U.S. and European daily implied volatility indices using six models, i.e., of $VIX$ indices extended to a range of markets. They found that while point and interval forecasts offered a good predictive accuracy under statistical loss functions, particularly when vector autoregressive models employing principal components were used, these same predictions turned out to be unable to support economically significant OOS results. For instance, they failed to yield abnormal profits when used to support simple volatility trading strategies.

In the body of research that has used the same type of models explored in this paper, Fernandes et al. [16] examined the daily predictability of the US $VIX$ index using a range of $HAR$-type models. They reported that it was challenging to outperform a pure $HAR$ model in terms of forecasting performance, especially at the shortest horizons. Psaradellis and Sermpinis [8] confirmed these findings and reported that a $HAR$-$GASVR$ (genetic algorithm-support vector regression) hybrid model and a $GASVR$ regression applied to the residuals of the pure $HAR$ model produced the most accurate forecasts and the largest economic value. However, these findings did open the door to the idea that by extending $HAR$ models to include sufficient levels of complexity, forecasting accuracy might have been gained. Our paper returns on this idea but works by including lags of the target variables within the $HAR$ framework to yield a stronger predictive power.

Byun and Kim [23] opened to the idea that combining information from both (conditional) variance and skeweness may offer optimal predictions of the former. Therefore they investigated the impact of risk-neutral skewness on forecasting future, spot realized volatility. They compared three models and found that the coefficient of risk-neutral skewness is positive and significant for daily volatility forecasts, with its significance decreasing with the forecasting horizon. They adopted the methodology to compute the risk-neutral skewness proposed by Bakshi et al. [24], which is the same methodology that underlies the construction of the $SKEW$ index we follow. Their OOS analysis showed that risk-neutral skewness significantly improves forecasting performance for daily, weekly, and monthly horizons. In our empirical applications, we also pursue this idea, even though our effort is entirely devote to the prediction of risk-neutral moments, and hence of the $VIX$ instead of spot, realized variance.

Liu and Faff [25] have developed SIX, an index forecasting 30-day realized skewness of the spot returns on the S&P 500. They compared SIX with the CBOE $SKEW$ index and a non-parametric measure, studying their relationship with the $VIX$ index. They argued that daily changes in $VIX$ and SIX are negatively related as opposed to changes in $VIX$ and $SKEW$, which are instead positively related. They also examined the predictive power of SIX and $SKEW$ for one-day-ahead S&P 500 returns, concluding that a lower risk-neutral skewness implies (more) negative returns. Lastly, they explored whether risk-neutral skewness may forecast the future physical skewness of the S&P 500 return distribution and found that SIX efficiently predicts future realized (physical) skewness, whereas $SKEW$ did not. Even though our paper only focuses on the prediction of the $VIX$ and of the $SKEW$ and hence of risk-neutral, option-implied moments and refrains from embarking in the prediction of realized quantities, a few of the ideas in Liu and Faff [25] are pursued and expanded further.

Recently, Bevilacqua and Tunaru [19], also extending the work by Liu and Faff [25], have decomposed the CBOE $SKEW$ index into positive ($SKE{W}^{+}$) and negative ($SKE{W}^{-}$) skewness. They find $SKE{W}^{-}$ to be a strong predictor of market downturns and recessions, while $SKE{W}^{+}$ would capture market sentiment (see also Han [26]). Additionally, $SKE{W}^{+}$ dampens the actual level of risk conveyed by $SKE{W}^{-}$ alone, thereby resulting in an underestimation of actual tail-risk. In studying the OOS one step-ahead forecasting performance of the resulting models, Bevilacqua and Tunaru [19] compare an $AR\left(1\right)$ benchmark model to three competing models. These models include the $AR\left(1\right)$ model augmented with the lagged values of $SKEW$, $SKE{W}^{-}$, and $SKE{W}^{+}$ as additional regressors. They report that all three competing models outperformed the benchmark $AR\left(1\right)$ model in terms of forecasting accuracy. Although we work within the tradition of using $HAR$ models as a predictive framework, we also pursue the idea that $SKE{W}^{-}$ may contain different and additional forecasting power vs. the overall $SKEW$ index. Byun and Kim [23] have focused on the predictive power of risk-neutral skewness for future volatility. They assess whether incorporating this skewness information into forecasting models can improve predictions of future realized volatility.

While the CBOE $SKEW$ index is a specific measure of perceived tail risk and skewness in the S&P 500 options market, the paper by Byun and Kim broadly examines risk-neutral skewness across different contexts and its role in forecasting volatility, rather than focusing exclusively on the $SKEW$ index itself. Therefore, the paper does not primarily aim to forecast the CBOE $SKEW$ index but rather uses skewness as a factor to enhance volatility forecasting models.

Finally, in this paper we test the separate predictive power of $SKEW$-$PUTS$, the risk-neutral measure of skewness that can be derived from using put options only. Doran et al. [27] had already investigated the existence of predictive informational content of the implied volatility skew using options on the S&P 100 index, evaluating whether the shape of the volatility skew could be used in forecasting future market movements. Using differences in implied volatilities between OTM and In-The-Money (ITM) options for both call and put options, they found that the put volatility skew exhibited strong predictive power for short-term market crashes. In contrast, the call volatility skew displayed weaker predictive power for short-term upward spikes in volatility. We extend this evidence on the separate market value of risk-neutral skewness extracted from put vs. call options to the problem of forecasting the $VIX$, hence risk-neutral volatility, instead of aggregate stock market returns.

3. Methodology

3.1. The Heterogeneous Autoregressive Models

The heterogeneous autoregressive model was proposed by Corsi [6] with the goal to provide accurate forecasts of (spot) realized volatilities, which in this paper we extend to the modelling and forecasting of The $HAR$ model is well known to be able to capture the highly implied persistent nature of volatility and it is logically grounded on the heterogeneous market hypothesis of Müller et al. [28], by which there would exist (at least) three types of market participants: short-term (characterized by high trading frequencies, typically daily), medium-term (who have longer trading horizons, typically weekly, including some types of speculative funds), and long-term (who generally have monthly trading frequencies, for instance institutional investors) traders and investors. The heterogeneity derives from the different time horizons of market participants and it leads to three distinct types of volatility: daily, weekly, and monthly. This creates an overall volatility cascade from low to high frequencies. Corsi [6] shows that this structure for the volatility cascade leads to a simple, restricted linear autoregressive model that employs volatilities realized over different time horizons. Typically, the $HAR$ model features daily, weekly, and monthly moving averages of realized volatilities to replicate the three components of volatility associated with different types of trading activity. In this paper, we follow the approach of Fernandes et al. [16], expanding the original $HAR$ model by also including biweekly and quarterly components of implied volatility. Hence, the $HAR$ model specification we employ is:

$$Y_t={\beta }_0+{\beta }_1{Y}_{t-1}+{\beta }_{5d}{Y}_{5d}+{\beta }_{10d}{Y}_{10d}+{\beta }_{21d}{Y}_{21d}+{\beta }_{63d}{Y}_{63d}+{ϵ}_t$$

(1)

where $Y_t$ is daily implied volatility recorded at time t and:

$$Y_{5d}={\displaystyle \frac{Y_{t-1}+Y_{t-2}+Y_{t-3}+Y_{t-4}+Y_{t-5}}{5}}$$

(2)

$$Y_{10d}={\displaystyle \frac{Y_{t-1}+Y_{t-2}+\cdots +Y_{t-9}+Y_{t-10}}{10}}$$

(3)

$$Y_{21d}={\displaystyle \frac{Y_{t-1}+Y_{t-2}+\cdots +Y_{t-20}+Y_{t-21}}{21}}$$

(4)

$$Y_{63d}={\displaystyle \frac{Y_{t-1}+Y_{t-2}+\cdots +Y_{t-62}+Y_{t-63}}{63}}$$

(5)

In words, the implied option volatilities associated with the weekly, biweekly, monthly, and quarterly components are represented by the 5-day, 10-day, 21-day, and 63-day moving averages, respectively, assuming 5 trading days per week and 21 trading days per month.

The main benefit of the $HAR$ model is its ability to capture complex long-range dependence in time series of volatility estimates while preserving a manageable number of parameters. In contrast, because of their long memory and the corresponding slow decay of the SAC and PAC functions, when modeling time series of both realized and implied volatilities using $ARMA$-type models, the number of parameters in both the autoregressive and moving average components tends to be very large. Additionally, Fernandes et al. [16] found that Autoregressive Fractionally Integrated Moving Average ($ARFIMA$) models perform very poorly both in-sample and out-of-sample when analyzing the time series of asset returns volatility because these models impose a linear form of long memory that depends exclusively on a single parameter.

The objective of this paper is however to extend the $HAR$ modelling strategy to involve the $SKEW$ index, both as a potential predictor of the $VIX$ and as a target variable to be forecast on its own. As discussed in the Introduction, the results in Byun and Kim [23] showed that risk-neutral skewness significantly enhances the forecasting accuracy of $HAR$ models for realized volatility across daily, weekly, and monthly horizons. Building on this finding, we employ a similar framework but adapted to predict implied volatility rather than realized volatility. This model, that we shall call $HAR$-$SK$, integrates daily, weekly, and monthly moving averages of realized volatilities, supplemented by one additional predictor: the lagged value of risk-neutral skewness, i.e.,

$$\begin{array}{c}VI{X}_t={\beta }_0+{\beta }_{1}VI{X}_{t-1}+{\beta }_{5d}VI{X}_{5d}+{\beta }_{10d}VI{X}_{10d}\\ +{\beta }_{21d}VI{X}_{21d}+{\beta }_{63d}VI{X}_{63d}+{\beta }_{\mathrm{SKEW}}SKE{W}_{t-1}+{ϵ}_t\end{array}$$

(6)

Furthermore, also following Bevilacqua and Tunaru [19], to explore the distinct informational content of the $SKEW$ index and of its negative component calculated using put options only ($SKE{W}^{-}$), we introduce a similar model named $HAR$-$SK$-$PUTS$. This model replaces the lagged value of $SKE{W}^{-}$ for the lagged value of the $SKEW$ index in the $HAR$-$SK$ model:

$$\begin{array}{cc}VI{X}_t\hfill & ={\beta }_0+{\beta }_{1}VI{X}_{t-1}+{\beta }_{5d}VI{X}_{5d}+{\beta }_{10d}VI{X}_{10d}+{\beta }_{21d}VI{X}_{21d}\hfill \\ \hfill & +{\beta }_{63d}VI{X}_{63d}+{\beta }_{\mathrm{SKEW}-\mathrm{PUTS}}SKE{W}_{t-1}^{-}+{ϵ}_t,\hfill \end{array}$$

(7)

where $SKE{W}_{t-1}^{-}$ indicates the lagged value of $SKEW$ calculated using put options only.

All the models specified above are used in modelling and forecasting the daily behaviour of the $VIX$ index. We shall also apply the $HAR$ model to forecast the time series of the $SKEW$ index. To investigate whether the negative part of the $SKEW$ is a significant driver of its evolution, we also adopt the hybrid $HAR$-$SK$-$PUTS$ model mentioned earlier:

$$\begin{array}{cc}SKE{W}_t\hfill & ={\beta }_0+{\beta }_{1}SKE{W}_{t-1}+{\beta }_{5d}SKE{W}_{5d}+{\beta }_{10d}SKE{W}_{10d}\hfill \\ \hfill & +{\beta }_{21d}SKE{W}_{21d}+{\beta }_{63d}SKE{W}_{63d}+{\beta }_{\mathrm{SKEW}-\mathrm{PUTS}}SKE{W}_{t-1}^{-}+{ϵ}_t\hfill \end{array}$$

(8)

Finally, to explore the relationship between the $VIX$ and $SKEW$ indices, we estimate and produce recursive forecasts using one additional model, the $HAR$-$IV$, as in Byun and Kim [23]. This model extends the pure $HAR$ model applied to forecast the $SKEW$ by including the lagged value of the $VIX$ index as an additional predictor:

$$\begin{array}{cc}SKE{W}_t\hfill & ={\beta }_0+{\beta }_{1}SKE{W}_{t-1}+{\beta }_{5d}SKE{W}_{5d}+{\beta }_{10d}SKE{W}_{10d}\hfill \\ \hfill & +{\beta }_{21d}SKE{W}_{21d}+{\beta }_{63d}SKE{W}_{63d}+{\beta }_{\mathrm{VIX}}VI{X}_{t-1}+{ϵ}_t.\hfill \end{array}$$

(9)

3.2. Research Design

The core of our empirical study consists of comparing the forecasting power of the models described above. Following the set up in Ahoniemi [15], we employ a 1000-day rolling window recursive forecasting scheme. On each day starting with observation 1001, We compute the 1-, 5-, 10-, and 21-day ahead forecasts. To produce predictions beyond the 1-day horizon, we employ a classical, indirect dynamic forecast scheme, producing the forecast for the $t+1$ period and treating this as actual data to compute the $t+2$ prediction, and so on.

To evaluate the forecasting accuracy of the various models introduced earlier, three main indicators shall be considered: the $MSFE$, the Mean Absolute Forecast Error ($MAFE$), and the percentage of correct predicted change of direction. The latter indicator is particularly useful from a trader’s point of view because many trading strategies rely on the correct direction prediction rather than on the magnitude of such changes (see, e.g., Costantini et al. [29] and Leitch and Tanner [30]).

We define the point forecast error for every observation in our (pseudo) OOS period, variable j, forecast horizon h, and model $M_i$ at time t as:

$$\begin{array}{c}{\widehat{e}}_{t,t+h}^{j,M_i}=y_{t+h}^j-{\widehat{y}}_{t,t+h}^{j,M_i},\end{array}$$

(10)

where $y_{t+h}^j$ is the actual value of variable j observed at time $t+h$ and ${\widehat{y}}_{t,t+h}^{j,M_i}$ is the predicted value as of time t for the variable j at time $t+h$.

Next, we define $MSFE$ as:

$$\begin{array}{c}MSF{E}^{j,h,M_i}={\displaystyle \frac{1}{T-h}}\sum _{t=1}^{T-h}{\left({\widehat{e}}_{t,t+h}^{j,M_i}\right)}^2,\end{array}$$

(11)

and $MAFE$ as

$$\begin{array}{c}MAF{E}^{j,h,M_i}={\displaystyle \frac{1}{T-h}}\sum _{t=1}^{T-h}\left|{\widehat{e}}_{t,t+h}^{j,M_i}\right|,\end{array}$$

(12)

where T is the size of the out-of-sample period and t the first observation of such out-of-sample period. Lastly, the percentage of correct change of direction is defined as:

$$\begin{array}{c}\mathrm{Correct}\mathrm{change}\mathrm{of}\mathrm{direction}={\displaystyle \frac{1}{T-h}}\sum _{t=1}^{T-h}{x}_{t,t+h}\end{array}$$

(13)

where

$$\begin{array}{c}{x}_{t,t+h}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}(y_{t+h}^j-y_t^j)({\widehat{y}}_{t,t+h}^{j,M_i}-y_t^j)>0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(14)

The indicators introduced above are not particularly explanatory by themselves, specifically in assessing whether differences between different models are significant. In order to assess the different OOS forecasting performances of the various models, we rely on two different tests.

The first type of testing framework employed is Hansen’s test for Superior Predictive Ability Hansen [31]. In this framework, all alternative forecasts are compared to a benchmark prediction model. Forecasting power is defined by expected loss. The framework is the same as White’s Reality Check (RC) test White [32] but Hansen’s uses a different (studentized) test statistic and a sample-dependent distribution under the null hypothesis. Compared to White’s test, Hansen’s SPA test is more powerful and less sensitive to irrelevant alternative models, see Hansen [31]. In this framework, our rolling estimation scheme is acceptable when comparison of forecasts is interpreted as being conditional on the estimated parameters. The composite null hypothesis tested using Hansen’s SPA test is that the benchmark model is not inferior to any alternative models, namely:

$$\begin{array}{c}{H}_0:\mathbb{E}\left(d_t\right)\le 0,\end{array}$$

(15)

where $\mathbb{E}\left(d_t\right)$ is the vector of expected values of the distances in performances $d_t={(d_{1,t},d_{2,t},\dots ,d_{m,t})}^{\prime }$ defined as:

$$\begin{array}{c}{d}_{k,t}=L(Y_t,{\delta }_{0,t-h})-L(Y_t,{\delta }_{k,t-h}),\end{array}$$

(16)

where $L(Y_t,{\delta }_{0,t-h})$ is the loss function of a selected benchmark model ${\delta }_{0,t-h}$ for the prediction made h-steps in advance for the target variable $Y_t$ and $L(Y_t,{\delta }_{k,t-h})$ is the loss function for the k-th alternative model with $k=1,2,\dots ,m$. The distribution of the test statistic is approximated by implementing the stationary bootstrapping method of Politis and Romano [33].

Our second framework to test superior predictive accuracy is Giacomini and White [34]’s who proposed an alternative approach after emphasizing that other frameworks for evaluating alternative OOS performances may not be necessarily useful in real-time forecast selection. Giacomini and White’s (GW) approach is based on inference about conditional expectations of forecasts and forecast errors rather than unconditional expectations. The authors argue that their framework assesses forecasting methods rather than forecasting models, thereby including evaluation choices about estimation procedures. Another advantage of GW test is that it can be applied to both nested and non-nested models, for multi-step point, interval, and density forecasting without a need for major adjustments. In fact, a fundamental assumption of GW is the heterogeneity of the process of the time series to be predicted instead of its stationarity. (The source of heterogeneity considered by GW is that induced by a distribution that changes over time. Furthermore, even if the underlying economic processes were stationary, heterogeneity in the observed time series can arise from changes in the measurement process Giacomini and White [34]. In this context, they suggest a rolling estimation scheme rather than an expanding window one).

The null hypothesis tested in GW’s approach is specified as:

$$\begin{array}{c}{H}_0:\mathbb{E}\left[L_{t+\tau }(Y_{t+\tau },{\widehat{f}}_{m,t})-L_{t+\tau }\left(Y_{t+\tau },{\widehat{g}}_{m,t}\right)|F_t\right]=0\end{array}$$

(17)

where $L_{t+\tau }(Y_{t+\tau },{\widehat{f}}_{m,t})$ and $L_{t+\tau }\left(Y_{t+\tau },{\widehat{g}}_{m,t}\right)$ are the loss functions relative to the predicted variables from forecasting methods f and g. $F_t$ is defined as $F_t\equiv [1L_{t+\tau -1}({\widehat{Y}}_{t+\tau -1},{\widehat{f}}_{m,t-1})$$-L_{t+\tau -1}({\widehat{Y}}_{t+\tau -1},{\widehat{g}}_{m,t-1})]$. If the null hypothesis fails to be rejected, this implies that the two models evaluated are, on average, equally accurate. A positive sign of the test statistics indicates that model f forecasts losses are, on average, larger than those produced by model g, while a negative sign of the test statistics indicates the opposite.

In the final part of the paper, we also resort to Hansen et al. [35]’s model confidence set methodology to identify a set of models that are statistically indistinguishable in terms of predictive accuracy. The MCS process starts with a broad set of candidate models and iteratively eliminates them based on performance, until reaching a set of models that are equivalent in terms of forecasting precision. It is interesting to note how this approach differs from the SPA method described above, which focuses on detecting the superiority of a specific model compared to a benchmark, evaluating whether a single model has significantly better predictive ability. Essentially, while the MCS seeks to find a group of equally valid models, the SPA test aims to identify the best model, comparing it with a benchmark and rejecting the others. Even though the MCS method is particularly useful when there is no clear ”winner” among the models, as it provides a robust set of optimal models, if it yields a set composed of only one model (i.e., a singleton), this would clearly be a very strong indication in favor of that model.

4. Data

4.1. Original Data and Filtering Procedures

The data consist of daily observations for $VIX$, $SKEW$, and $SKE{W}^{-}$ indices, over the period 4 January 1996–31 December 2019, for a total of 6005 daily observations. These data come from different sources: the time series of the $VIX$ and $SKEW$ indices are collected from the CBOE website. To decompose the $SKEW$ index and extract only its negative part, we re-compute the entire daily history of the index. This involved running the filtering algorithm described in Appendix A and calculating afresh the $SKEW$ index using the prices of the power payoffs detailed in the appendix. The only difference vs. the official procedure adopted by the CBOE is that we apply this procedure considering put options only. Because we have S&P 500 options available only with reference to a 4 January 1996–31 December 2019 sample, we limit our analysis to this range of dates even though the full sample of $VIX$ and $SKEW$ data from the CBOE are in principle available for a longer, 2 January 1990–30 November 2020 sample. In any event, this sample provides more than six years of additional data vs. the sample covered by Fernandes et al. [16]. The options price data are sourced from the Optionmetrics database available through WRDS (Wharton Research Data Services).

To calculate the $SKE{W}^{-}$ index, besides the prices of S&P 500 put options, the only other additional of information needed are the rates implicit in the U.S. Treasury yield curve, taken as a proxy for the risk-free rates at a range of relevant maturities. Such yield data are collected from the Federal Reserve Bank Reports, also available in WRDS. In practice, the relevant Constant Maturity Treasury rates are interpolated using Nelson and Siegel [36]’s model, which consists of fitting a parametric model to yield data to obtain fitted (interpolated) continuously compounded rates for all relevant option maturities. The model parameters are estimated by converting observed market Treasury rates into discount factors and minimizing the sum of squared errors derived as a difference between the model-implied and market discount factors.

Once the novel $SKE{W}^{-}$ index has been calculated, we obtain the complete time series for all the three variables of interest. In principle, the series should contain a total of 6042 observations. However, the datasets must also be adjusted for missing data points by stepwise deletion, considering only those dates for which data is available for all series. The resulting, final dataset contains a total of 6005 observations.

4.2. Descriptive Statistics

Figure 1 and Figure 2 show the evolution of both the $VIX$ and $SKEW$ indices between 4 January 1996 and 31 December 2019. Because in what follows we shall forecast both the level of the indices and their natural logarithm, Figure 1 concerns levels, while Figure 2 covers the case of the log of the indices. Importantly, the relatively long sample period covered in our paper covers a series of financial crises and turmoil periods in the financial markets that were characterized by high levels of the $VIX$ index: the Asian crisis in 1997, the Russian crisis in 1998, the dot-com bubble burst in 2000, the 9/11 terrorist attack in 2001, the Lehman Brothers bankruptcy in September 2008, and the subsequent global financial crisis. From the plots, we observe that at least until the 2008 financial crisis, it seems that the $SKEW$ index tended to be increase contemporaneously with the events marked by spikes in the $VIX$. Starting in 2009, however, we note that an upward trend in the $SKEW$ sets in in spite of the fact that $VIX$ remains rather stable or even declines between 2011 and 2018. Hence Figure 1 fails to provide any clear insights concerning the overall relationship between the two implied volatility-based CBOE indices.

Table 1. Summary Statistics: The table reports the key sample statistics for the levels and the logarithms of the $VIX$ and $SKEW$ indices calculated on the full sample, from 4 January 1996 to 31 December 2019. The table also shows the p-values for the Jarque-Bera test of normality. Under the null hypothesis of zero skewness and of kurtosis equal to three, the Jarque-Bera test statistics is distributed as a chi-squared with 2 degrees of freedom.

	VIX	SKEW	LOG(VIX)	LOG(SKEW)
Mean	19.9434	120.5813	2.9255	4.7902
Median	18.3200	119.3200	2.9080	4.7818
Maximum	80.8600	159.0300	4.3927	5.0691
Minimum	9.1400	104.0900	2.2127	4.6453
Std. Dev.	8.1434	8.0220	0.3548	0.0648
Skewness	2.0343	0.9538	0.5673	0.7590
Kurtosis	10.1969	4.0143	3.2928	3.5171
Jarque-Bera	0.0000	0.0000	0.0000	0.0000
Observations	6005	6005	6005	6005

reports standard summary statistics for the $VIX$ and $SKEW$ indices and their logarithms. The distribution of the $VIX$ index turns out to be very skewed to the right and leptokurtic as the Jarque-Bera test strongly rejects the null hypothesis of normality. Applying the logarithm to the series of the $VIX$ removes as substantial portion of the non-normal features, but the Jarque-Bera test still rejects the null hypothesis of a normal distribution. Regarding the time series of the $SKEW$ index, it also displays positive skewness and a positive excess of kurtosis. However, in this case applying a log transformation has very weak effects on the measurable sample deviations from normality, at least in comparison to the $VIX$ index.

The (unreported, to save space) Sample Autocorrelation Function (SACF) and Sample Partial Autocorrelation Function (SPACF) of the series all display high levels of persistence, which may be due to the presence of a unit root. To address this, we conduct unit root tests, specifically the Augmented Dickey–Fuller (ADF) and Phillips–Perron (PP) tests for unit roots, and the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test for stationarity. The ADF and PP tests strongly reject the null hypothesis of a unit root, whereas the KPSS rejects the null of stationarity.

As it is well known, because structural change and unit roots are closely related phenomena, we should be mindful that conventional unit root tests are biased toward a false unit root null when the data are trend stationary but contain a structural break. However, unit root tests exist that remain valid in the presence of a break. In particular, we apply a modified augmented Dickey-Fuller tests implemented in two steps, which allow either a change in level or a change in both level and trend when the breaks occur instantaneously (additive breaks), which we consider more typical of financial markets. In the first-step ADF test, the number of lags is selected by minimization of standard information criteria (we cannot find any difference between BIC and Hannan-Quinn results), given a maximum of 33 lags. The break date is selected as the one providing the strongest evidence of a break when all the break dates after day 34 in the sample are tested. In all cases, we reject the null of a unit root in the series considered. For instance, in the case of the raw of $VIX$ index, the baseline ADF test stistic (with 7 lags selected by the BIC) is −5.142 (p-value of 0.000); when the occurrence of breaks the ADF test statistics declines to −7.818 (p-value of 0.000). The corresponding stastistics in the case of the $SKEW$ index are −4.503 and −14.543 (p-values are both 0.000). Detailed results are available upon request. This contradiction holds for both indices and persists even when the logarithmic transformation is applied. The tests are also applied to the first differences of the series. In this case, the KPSS test fails to reject the null hypothesis of stationarity for all series under analysis, the differenced series and the differenced logarithm series. All three tests lead to the same conclusion. To account for this, we conduct our out-of-sample analysis using both the original and differenced series.

5. Empirical Results

5.1. In-Sample Analysis

We start our analysis by inspecting the in-sample predictability properties of (a few transformations of) the time series of the $VIX$ and $SKEW$ indices.

5.1.1. In-Sample Forecasting of the $VIX$ Index

Because the literature has shown that $ARIMA$ and $ARFIMA$ fail to outperform the $HAR$ model and its extensions (see Fernandes et al. [16]), to save space in the following we shall focus exclusively on the $HAR$ model in both in-sample and out-of-sample analyses. This choice aligns with our objective to explore whether the $SKEW$ index offers any additional information for modeling and forecasting the $VIX$ index beyond what the $HAR$ model already incorporates. The other predictive models employed are the $HAR$-$SK$ and the $HAR$-$SK$-$PUTS$. These models are applied to both the time series of the $VIX$ index as well as to its logarithmic transformation. (For consistency, when modelling the logarithm of the $VIX$, the $HAR$-$SK$ and $HAR$-$SK$-$PUTS$ models include as additional regressors the lagged value of the logarithm of the $SKEW$ index and of the $SKE{W}^{-}$ respectively).

Table 2. In-sample modelling of $VIX$ and $LOG\left(VIX\right)$ including all predictors in Equations (1), (6), and (7). p-values are reported in parenthesis. In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	VIX			LOG(VIX)
	HAR	HAR-SK	HAR-SK-PUTS	HAR	HAR-SK	HAR-SK-PUTS
${\beta }_0$	0.2454	0.6436	0.3832	0.0283	0.1083	0.0709
	(0.000) ***	(0.083) *	(0.178)	(0.000) ***	(0.159)	(0.258)
${\beta }_1$	0.8623	0.8614	0.8616	0.8985	0.8975	0.8973
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{5d}$	0.0189	0.0190	0.0190	−0.0050	−0.0048	−0.0047
	(0.479)	(0.476)	(0.476)	(0.847)	(0.853)	(0.855)
${\beta }_{10d}$	0.1175	0.1176	0.1175	0.0829	0.0831	0.0830
	(0.000) ***	(0.000) ***	(0.000) ***	(0.004) ***	(0.004) ***	(0.004) ***
${\beta }_{21d}$	−0.0191	−0.0187	−0.01889	−0.0049	−0.0048	−0.0047
	(0.335)	(0.345)	(0.341)	(0.796)	(0.801)	(0.805)
${\beta }_{63d}$	0.0081	0.0071	0.0078	0.0188	0.0178	0.0184
	(0.311)	(0.379)	(0.332)	(0.022) **	(0.031) **	(0.025) **
${\beta }_{SKEW}$	−0.0031			−0.0158
	(0.277)			(0.296)
${\beta }_{SKEW-PUTS}$			−0.0009			−0.0079
			(0.620)			(0.493)
$R_{adjusted}^2$	0.96118	0.96115	0.96115	0.96495	0.96495	0.96495
	$\mathbf{\Delta }$ VIX			$\mathbf{\Delta }$ LOG(VIX)
	HAR	HAR-SK	HAR-SK-PUTS	HAR	HAR-SK	HAR-SK-PUTS
${\beta }_0$	−0.0009	−0.0009	−0.0009	−0.0001	−0.0001	−0.0001
	(0.965)	(0.964)	(0.965)	(0.942)	(0.941)	(0.942)
${\beta }_1$	−0.0516	−0.0485	−0.0463	−0.0358	−0.0332	−0.0321
	(0.000) ***	(0.001) ***	(0.001) ***	(0.011) **	(0.020) **	(0.025) **
${\beta }_{5d}$	−0.2327	−0.2328	−0.2333	−0.1301	−0.1304	−0.1306
	(0.000) ***	(0.000) ***	(0.000) ***	(0.004) ***	(0.003) ***	(0.003) ***
${\beta }_{10d}$	−0.1049	−0.1069	−0.1072	−0.1943	−0.1959	−0.1957
	(0.195)	(0.187)	(0.185)	(0.012) **	(0.011) **	(0.011) **
${\beta }_{21d}$	−0.0663	−0.0668	−0.0655	−0.0981	−0.0977	−0.0977
	(0.589)	(0.586)	(0.593)	(0.408)	(0.410)	(0.410)
${\beta }_{63d}$	−0.7251	−0.7244	−0.7261	−0.8042	−0.8039	−0.8048
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{SKEW}$		0.0109			0.0443
		(0.128)			(0.236)
${\beta }_{SKEW-PUTS}$			0.0059			0.0229
			(0.032) **			(0.180)
$R_{adjusted}^2$	0.0263	0.0267	0.0270	0.0213	0.0215	0.0216

reports the estimated coefficients of the three models obtained over the whole sample. The first insight we can draw is that neither $SKEW$ nor $SKE{W}^{-}$ are characterized by statistically significant coefficients, when these are used as predictors. Hence, their addition to the pure $HAR$ model does not add any explanatory power in terms of adjusted $R^2$. Furthermore, among all the specified regressors, only the coefficients of the lagged value of the $VIX$ and the 10-day average are significant. Only when modeling the logarithm of the $VIX$, the coefficient of the 63-day average is also significant at the 5% test size level. However, this result may be influenced by high levels of linear dependence among the regressors (i.e., multicollinearity) that delivers a good adjusted $R^2$ (in fact, always in excess of 96 percent in Table 2) but modest significance of the individual coefficients. The lower panels of Table 2) devoted to the estimation of models for changes in $VIX$ and log-$VIX$, apart from a structurally lower R-squared vs. the case of levels and log levels (2–3% vs. 96–97%, as one would expect), disclose the same rates of statistical significance as the upper panels, but now tilted away from the first lag of the series and 10-day volatility, towards 5- and 63-days.

To address these issues in our OOS analysis, we shall consider two specifications of the $HAR$ model (and hence of the $HAR$-$SK$ and $HAR$-$SK$-$PUTS$ models): one which includes all the regressors, and one which includes only those predictors that display significant coefficients in the in-sample specification. Even though this is clearly sub-optimal and may imply some degree of hindsight bias, it provides an upper bound to the forecasting performance of the models when a substantial part of their tendency to over-fit is removed from the problem. (The first-best solution would consist of re-estimating and re-testing for significance the coefficients of all models in a recursive fashion throughout the pseudo OOS exercise). By re-estimating the models using only the significant predictors, we aim at obtaining additional evidence of the robustness of our baseline findings that use predictions from the full models that include all predictors. However, in unreported tables, we observe that even in the parsimonious specifications, the inclusion of the lagged values of the $SKEW$ index and of the $SKE{W}^{-}$ does not add any explanatory power to the pure $HAR$ model and that their coefficients are not significantly different from zero. (Complete and tabulated results concerning the case of the models trimmed of the insignificant coefficients obtained over the full sample are available from the authors upon request).

5.1.2. In-Sample Forecasting of the $SKEW$ Index

As discussed in Section 4 earlier, the time series of the $SKEW$ index displays, as it happens for the $VIX$ index, high levels of persistence. Furthermore, from the estimation of $ARMA$ type models, we learn that model selection based on standard information criteria typically leads to highly parameterized models, including dozens of autoregressive and moving average terms. Because of this, we also employ the $HAR$ model directly for the in-sample analysis of the $SKEW$ index. The predictors adopted are the same employed for the $VIX$ index, as specified in Equation (1). Moreover, by using a decomposition of the $SKEW$ index, Bevilacqua and Tunaru [19] have argued that $SKE{W}^{-}$ plays a major role in driving the evolution of the overall $SKEW$ index, we also implement the $HAR$-$SK$-$PUTS$ specification in Equation (8). Finally, as there is mixed evidence in the literature about the usefulness of the $SKEW$ in modelling the time series of the $VIX$ index, we also test whether the opposite relationship is somewhat clearer and we recursively estimate the $HAR$-$IV$ model in Equation (9).

Table 3. In-sample modelling of $SKEW$ and $LOG\left(SKEW\right)$ including all the regressors specified in Equations (1), (8), and (9). p-values are reported in parenthesis. In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	SKEW			LOG(SKEW)
	HAR	HAR-SK-PUTS	HAR-IV	HAR	HAR-SK-PUTS	HAR-IV
${\beta }_0$	2.3106	2.6646	2.6718	0.0878	0.0951	0.1152
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_1$	0.6231	0.6115	0.6226	0.6180	0.6082	0.6168
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{5d}$	0.2666	0.2643	0.2667	0.2786	0.2768	0.2787
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{10d}$	−0.0900	−0.0913	−0.0911	−0.091	−0.0922	−0.0930
	(0.012) **	(0.011) **	(0.011) **	(0.011) **	(0.010) **	(0.010) ***
${\beta }_{21d}$	0.1574	0.1590	0.1584	0.1490	0.1506	0.1505
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{63d}$	0.0239	0.0204	0.0221	0.0271	0.0244	0.0240
	(0.107)	(0.170)	(0.140)	(0.065) *	(0.097) *	(0.105)
${\beta }_{SKEW-PUTS}$		0.0119			0.0118
		(0.010) ***			(0.035) **
${\beta }_{VIX}$			−0.0046			−0.0016
			(0.337)			(0.089) *
$R_{adjusted}^2$	0.88068	0.88079	0.88068	0.88447	0.88454	0.88451
	$\Delta$ SKEW			$\Delta$ LOG(SKEW)
	HAR	HAR-SK-PUTS	HAR-IV	HAR	HAR-SK-PUTS	HAR-IV
${\beta }_0$	0.0146	0.0146	0.0146	0.0001	0.0001	0.0001
	(0.685)	(0.684)	(0.686)	(0.692)	(0.684)	(0.685)
${\beta }_1$	−0.2022	−0.2017	−0.2016	−0.2081	−0.2070	−0.2078
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{5d}$	−0.1832	−0.1831	−0.1830	−0.1847	−0.1845	−0.1846
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{10d}$	−0.5442	−0.5442	−0.5440	−0.5342	−0.5342	−0.5340
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{21d}$	−0.7632	−0.7632	−0.7636	−0.7093	−0.7092	−0.7093
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
${\beta }_{63d}$	−1.2132	−1.2129	−1.2127	−1.2832	−1.2827	−1.2832
	(0.001) **	(0.001) **	(0.001) **	(0.001) ***	(0.001) ***	(0.001) ***
${\beta }_{SKEW-PUTS}$		−0.0006			−0.0015
		(0.909)			(0.794)
${\beta }_{VIX}$			0.0076			0.0007
			(0.733)			(0.876)
$R_{adjusted}^2$	0.1204	0.1204	0.1204	0.1220	0.1220	0.1220

shows the results of parameters estimation over the full sample of the three models cited above, both for the $SKEW$ index and for its logarithm. Consistently to what was done in the case of the $VIX$, when the $SKEW$ is taken in logs, the additional predictors ($VIX$ and $SKE{W}^{-}$) are taken in logs too. The coefficient estimates show that the only non-significant coefficient in the $HAR$ model is the one associated to the quarterly moving average. As for the overall goodness of fit, also in Table 3, the $HAR$ model produces satisfactory results in terms of adjusted $R^2$, always in excess of 88 percent. In the case of the $HAR$-$SK$-$PUTS$ and the $HAR$-$IV$ specifications, we observe that the first predictor adds some explanatory power in terms of adjusted $R^2$, while the second specification does not. Hence, it seems that the $VIX$ index does not contain useful information for the in-sample fit of the $SKEW$ index. Indeed, the coefficient related to the $SKE{W}^{-}$ is significantly different from zero while the coefficient associated to the $VIX$ index is not. However, by looking at the results produced by the in-sample analysis of the logarithm of the $SKEW$ index, the results change: in the pure $HAR$ model, all the predictors are characterized by significant coefficients and both the lagged value of the $SKE{W}^{-}$ and the lagged value of the $VIX$ improve the explanatory power of the $HAR$ model, although marginally (at least in terms of the in-sample, adjusted $R^2$). The two lower panels of Table 3, devoted to in-sample modelling of changes in the $SKEW$, apart from a comparable drop in R-squared coefficients (here from 88% to 12%), show qualitatively similar results in terms of widespread coefficient significance (apart from the loss of estimation accuracy of the first lag, as one would expect when the series are first-differenced). Yet, the point estimates tend to turn negative, as one would expect as a result of over-differencing the SKEW in levels.

In unreported results based on employing only those predictors that turned out to be significant in the $HAR$ model (thus excluding the 63-day moving average), the conclusions regarding the in-sample predictive power of the $VIX$ and $SKE{W}^{-}$ reported above are unchanged: the coefficient associated with $SKE{W}^{-}$ in the model for $HAR$-$SK$-$PUTS$ is significantly different from zero, while the coefficient associated with the $VIX$ index in the $HAR$-$IV$ model is not.

5.2. Out-of-Sample Analysis

After evaluating the models in terms of their in-sample forecasting accuracy, as measured by their $R^2$, we focus our attention on a pseudo out-of-sample analysis. The pseudo nature of the exercise stems from the fact that the predictive $HAR$ models have been selected and fine tuned in Section 5.1.2 using the full sample of data. Nonetheless the recursive implementation of the OOS forecasting scheme described in Section 3.2 guarantees that no hindsight bias occurs in the calculation of the point forecasts used to assess the accuracy of alternative models. As planned, in what follows, we present the results obtained in terms of $MSFE$, $MAFE$, and the percentage of correctly predicted changes in direction. The OOS forecasting performance of the various models are evaluated using Hansen’s Hansen [31] test for superior predictive ability and GW’s Giacomini and White [34] test for conditional predictive ability. As a robustness check of our results, and given their persistence, we also employ the models to predict the first differences of the variables.

All the results that follow refer to the prediction of the level of the $VIX$ and $SKEW$ indices. Hence, when we employ transformations of such variables, we first transform back the outputs of the models in levels of the $VIX$ and $SKEW$ before calculating the indicators of forecasting accuracy. This choice is justified by the fact that the results obtained in this manner have a more solid practical meaning when compared to those derived from transformed variables.

5.2.1. Forecasting the $VIX$ Index

Table 4. VIX and $LOG\left(VIX\right)$ out-of-sample forecasting performance; models specified using all the regressors detailed in Equations (1), (6), and (7). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: VIX	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	2.83505	0.0010 ***	1.02704	0.0208 **	48.071%
AR(1)	2.79258	0.0354 **	1.02883	0.0778 *	48.071%
HAR	2.77097	1.0000	1.02599	1.0000	48.151%
HAR-SK	2.77736	0.1628	1.02861	0.0348 **	48.091%
HAR-SK-PUTS	2.77991	0.1616	1.02884	0.0020 ***	48.231%
5 Days Ahead
RW	9.73185	0.0000 ***	2.01384	0.0088 ***	49.770%
AR(1)	9.43394	0.0000 ***	2.03524	0.0081 ***	49.510%
HAR	8.45600	1.0000	1.92459	1.0000	49.750%
HAR-SK	8.55837	0.0512 *	1.93905	0.0548 *	49.790%
HAR-SK-PUTS	8.58210	0.0058 ***	1.94787	0.0002 ***	49.750%
10 Days Ahead
RW	14.90454	0.0000 ***	2.52102	0.0007 ***	50.850%
AR(1)	14.22340	0.0000 ***	2.53821	0.0005 ***	51.150%
HAR	10.9335	1.0000	2.25484	1.0000	50.100%
HAR-SK	11.1682	0.0216 **	2.28500	0.0196 **	50.160%
HAR-SK-PUTS	11.2505	0.001 ***	2.30479	0.0000 ***	50.040%
21 Days Ahead
RW	26.81030	0.0000 ***	3.34224	0.0000 ***	49.789%
AR(1)	24.58905	0.0000 ***	3.38506	0.0000 ***	49.970%
HAR	12.7081	1.0000	2.50222	1.0000	48.164%
HAR-SK	13.20000	0.006 ***	2.55884	0.0032 ***	48.545%
HAR-SK-PUTS	13.4585	0.0000 ***	2.59602	0.0004 ***	48.364%
Panel B: LOG(VIX)	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	2.80133	0.0017 ***	1.02091	0.0243 **	48.071%
AR(1)	2.76920	0.0184 **	1.01743	0.0394 **	48.471%
HAR	2.72317	1.0000	1.01171	1.0000	48.171%
HAR-SK	2.72818	0.3314	1.01243	0.3734	48.151%
HAR-SK-PUTS	2.73619	0.1304	1.01227	0.4496	48.411%
5 Days Ahead
RW	9.73185	0.0000 ***	2.01384	0.0075 **	49.230%
AR(1)	9.36455	0.0000 ***	1.99364	0.0114 **	49.490%
HAR	8.50527	1.0000	1.90265	1.0000	49.430%
HAR-SK	8.55933	0.3136	1.90857	0.2724	49.610%
HAR-SK-PUTS	8.62174	0.1098	1.91658	0.0240 **	49.710%
10 Days Ahead
RW	14.90454	0.0000 ***	2.52102	0.0007 **	50.901%
AR(1)	14.07639	0.0000 ***	2.46299	0.0006 **	51.181%
HAR	11.2748	1.0000	2.22273	1.0000	50.120%
HAR-SK	11.4111	0.2030	2.24202	0.0370 **	49.940%
HAR-SK-PUTS	11.5662	0.0520 *	2.25992	0.0020 ***	49.900%
21 Days Ahead
RW	26.81030	0.0000 ***	3.34224	0.0000 **	49.789%
AR(1)	24.14271	0.0000 ***	3.22153	0.0000 **	49.970%
HAR	13.3846	1.0000	2.44966	1.0000	48.304%
HAR-SK	13.7287	0.0382 **	2.49199	0.0070 ***	48.164%
HAR-SK-PUTS	14.0684	0.0110 **	2.53122	0.0004 ***	48.385%

and

Table 5. D(VIX) and $D\left(LOG\right(VIX\left)\right)$ out-of-sample forecasting performance; models specified using all the regressors detailed in Equations (1), (6), and (7). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: D(VIX)	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	2.80620	0.04420 **	1.02230	0.4780	48.321%
AR(1)	2.79258	0.03545 **	1.02883	0.0778 **	48.071%
HAR	2.80063	1.00000	1.02206	1.0000	48.221%
HAR-SK	2.80472	0.17340	1.02409	0.0016 ***	48.001%
HAR-SK-PUTS	2.80485	0.31820	1.02572	0.0108 ***	48.361%
5 Days Ahead
RW	3.38879	0.00580 ***	1.13876	0.0097 **	47.939%
AR(1)	9.43394	0.00000 ***	2.03524	0.0081 **	49.510%
HAR	2.80291	1.00000	1.02163	0.4430	48.579%
HAR-SK	2.80512	0.20240	1.02154	1.0000	48.419%
HAR-SK-PUTS	2.81609	0.00500 ***	1.02708	0.0002 ***	48.299%
10 Days Ahead
RW	3.11284	0.00000 ***	1.08347	0.0009 ***	47.927%
AR(1)	14.22340	0.00000 ***	2.53821	0.0005 ***	51.150%
HAR	2.79450	1.00000	1.02080	0.4532	48.508%
HAR-SK	2.79720	0.15540	1.02074	1.0000	48.448%
HAR-SK-PUTS	2.80761	0.00540 ***	1.02626	0.0008 ***	48.248%
21 Days Ahead
RW	3.05253	0.00000 ***	1.07647	0.0000 ***	47.927%
AR(1)	24.58905	0.00000 ***	3.38506	0.0000 ***	49.970%
HAR	2.78851	1.00000	1.02142	0.2554	48.535%
HAR-SK	2.79037	0.24460	1.02105	1.0000	48.234%
HAR-SK-PUTS	2.80247	0.00280 ***	1.02700	0.0002 ***	48.173%
Panel B: D(LOG(VIX))	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	7.13495	0.00000 ***	1.53507	0.0000 **	49.100%
AR(1)	2.76920	0.01840 **	1.01743	0.0411 **	48.471%
HAR	2.74492	0.48740	1.01354	1.0000	48.221%
HAR-SK	2.74757	0.34780	1.01454	0.0326 **	48.141%
HAR-SK-PUTS	2.74307	1.00000	1.01418	0.34420	48.421%
5 Days Ahead
RW	3.65120	0.00000 ***	1.14897	0.0064 ***	48.459%
AR(1)	9.36455	0.00000 ***	1.99364	0.01140 **	49.490%
HAR	2.74398	1.00000	1.01123	1.0000	48.279%
HAR-SK	2.74630	0.41880	1.01230	0.0284 **	48.139%
HAR-SK-PUTS	2.74534	0.53680	1.01188	0.3506	48.499%
10 Days Ahead
RW	3.20295	0.00000 ***	1.08643	0.0088 ***	48.087%
AR(1)	14.07639	0.00000 ***	2.46299	0.0000 ***	51.181%
HAR	2.74120	1.00000	1.01075	1.0000	48.288%
HAR-SK	2.74296	0.49580	1.01186	0.0184 **	48.207%
HAR-SK-PUTS	2.74309	0.47800	1.01151	0.3136	48.468%
21 Days Ahead
RW	3.00392	0.00060 ***	1.02554	0.0380 ***	48.314%
AR(1)	24.14271	0.00000 ***	3.22153	0.0000 ***	49.970%
HAR	2.74088	1.00000	1.01013	1.0000	48.294%
HAR-SK	2.74323	0.34840	1.01131	0.0116 **	48.254%
HAR-SK-PUTS	2.74226	0.51900	1.01076	0.3466	48.454%

report the results in terms of MSFE, MAFE, and the percentage of correct direction changes when forecasting the $VIX$ index using four different specifications of the target variable, i.e., levels, changes, logs, and change in logs. (The predictors in these tables include all the variables used in the in-sample analysis presented above, without excluding those found to be insignificant. Nonetheless we shall comment, where appropriate, on the main differences revealed by supplementary analyses based on models that have been trimmed, as described in Section 5.1.1 and Section 5.1.2, to exclude any insignificant predictors in the in-sample analyses). The tables also report the p-values of Hansen’s SPA test using both the squared error (SPA(SE)) and the absolute error (SPA(AE)) loss functions. A rejection of the null hypothesis in Hansen’s test as signalled by a low p-value indicates that the model taken as a benchmark is inferior to at least one of the alternative models treated in the set. Forecasts are recursively computed at the one-, five-, 10-, and 21-day ahead horizons.

Table 4 and Table 5 show that all in all, the results obtained in the literature about forecasting the $VIX$ index are confirmed by our findings: the pure $HAR$ model is very hard to beat. Specifically, in evaluating the results using Hansen’s SPA test for the squared error loss function, neither the $HAR$-$SK$ nor the $HAR$-$SK$-$PUTS$ models turn out to be clearly superior in any of the cases investigated. The same conclusion holds when employing Hansen’s test under the absolute error loss function. In fact, the results for SPA(AE) reveal that at several horizons (always in the case of the $VIX$ taken in levels), the $HAR$-$SK$ and $HAR$-$SK$-$PUTS$ models are rejected with small p-values from the model confidence set.

Furthermore, in forecasting the levels of the $VIX$ index (Panel A in Table 4), the $HAR$ model is clearly superior to its two extensions, except for the 1-day ahead forecast horizon. In predicting the log transformation of the $VIX$ (Panel B, Table 4), the superior forecasting performance of the simple $HAR$ model becomes more evident at the longest horizons. Specifically, the $HAR$-$SK$-$PUTS$ model performs the worst, being outperformed in terms of absolute errors at all forecast horizons except in the 1-day ahead set of exercises. Also in terms of squared errors, the $HAR$-$SK$-$PUTS$ model is outperformed at all forecast horizons except for the 1-day ahead horizon. Interestingly, forecasting the level of the $VIX$ at a 5-day horizon and longer, yields more accuracy under a squared loss function (the $MSFE$) vs. forecasting the log-$VIX$ and then transforming such predictions back to levels; yet, the opposite occurs in the case of absolute loss (the $MAFE$ criterion). This is probably caused by the fact that a squared loss is more influenced by extreme outlier prediction errors that tend to appear more frequently as a result of recursive model estimations based on the logarithm of the $VIX$ that get then compounded in the exponential transformation.

In both panels A and B of Table 4, the average correct change of direction statistics fails to give unequivocal indications, apart from showing that when the target variable is the level of the $VIX$, applying transformations appears to generally hurt performance. In particular, in panel A for horizons of five days and higher, $HAR$-$SK$ slightly outperforms. In panel B, every horizon implies a different ranking, even though the differences are small. As pointed out by one referee, our results on the correct change of direction tend to be below 50%, which is relatively disappointing in the light of previous literature. Casual inspection of the results when forecasting the $VIX$ and $SKEW$ in levels using the plain vanilla $HAR$ model, we note that it is the inclusion of the 2008–2010 data that erodes the correct change of direction performance of the $HAR$ model. For instance, in the case of the $VIX$, the fraction of correctly predicted signs exceeds 52% when these three years are excluded; the corresponding statistic is 51% in the case of $SKEW$ forecasts. Very similar conclusions are achieved with reference to the prediction of the first differences of the $VIX$ in Panel A of Table 5). The $HAR$-$SK$-$PUTS$ model performs the worst, being outperformed (i.e., excluded from the model confidence set) both in terms of squared and absolute errors at all the forecast horizons except for the 1 day-ahead (but just in terms of squared errors).

Panel B in Table 5 shows instead the results for predicting the first difference of the logarithm of the $VIX$ (i.e., its continuous growth rate). In this case, there is no clearly superior model in terms of squared errors, while the $HAR$-$SK$ model is outperformed at all forecast horizons in terms of absolute errors, i.e., when the models are ranked in terms of their $MAFE$. However, in panel B forecasting the first-difference of the logarithm of the $VIX$ always yields, despite the need of an exponential transformation to convert the predictions back to levels, more accurate predictions vs. the case in which one directly predicts the first-difference of the $VIX$. Hence, whether the exponentiation required by the prediction of variables in logs may eventually penalized realized predictive performance or not, remains unclear as probably data-dependent.

Also with reference to Table 5, when it comes to the average percentage of correct predicted changes in direction, none of the models tested yields particularly precise prediction results, as this percentage ranges between 48% and 49%. Hence, none of the models appears to be specifically suitable for setting up a purely directional trading strategy as they would fail to outperform a naive coin-tossing device. This is in some contrast with earlier literature and appears to depend on poor performance of all models during the Great Financial Crisis.

We have also examined the OOS realized accuracy of forecasting the $VIX$ index by including only those predictors that were found to be significant in the in-sample analysis in Section 5.1.1. In unreported tables, we find results that are qualitatively similar to the ones commented above. According to Hansen’s test applied to both the squared and absolute loss functions, the $HAR$ model is not outperformed by the other two at any of the forecast horizons. Moreover, in the case of the logarithm of the $VIX$, in terms of the absolute error loss function, the $HAR$-$SK$-$PUTS$ model is outperformed at all horizons except the 1-day. Nonetheless, also in this case, we find no structurally superior model when their performance is assessed under a squared loss function. Additionally, the forecasting results obtained from the first difference of both the $VIX$ and its logarithm show that the inclusion of non-significant predictors does not affect the OOS results: the $HAR$-$SK$-$PUTS$ model is outperformed in predicting the $VIX$ via its first difference (according to both the SPA (SE) and SPA (AE) tests), whereas the $HAR$-$SK$ model is inferior in predicting the logarithm of the $VIX$ (yet, according to the SPA (AE) test only). Once again, no particularly useful information emerges from the percentage of correct predicted changes in direction indicator.

The picture of the rankings of OOS performances that has emerged so far appears to be clear: the lagged value of the $SKEW$ index does not add any useful information to the OOS forecasting of the $VIX$ index beyond what is already accounted for by the pure $HAR$ model. This conclusion is robust across different specifications of the target variables predicted and across different specifications. These results extend one of the findings in Bevilacqua and Tunaru [19], who found that the inclusion of the $SKEW$ index, its “negative” part ($SKE{W}^{-}$), and its positive part ($SKE{W}^{+}$) produce superior out-of-sample forecasting performances compared to a simple $AR\left(1\right)$ process. Yet, we also found that when employing a more sophisticated model, the effect of the $SKEW$ index in predicting the $VIX$ vanishes. Also when including the the $SKEW$ index obtained from put options only, the predictive power of the $SKEW$ index in predicting the $VIX$ fails to improve and, in fact, occasionally worsens the OOS accuracy of the otherwise simpler, standard $HAR$.

Table 6. Giacomini-White’s type test for conditional predictive ability comparing $VIX$ and $LOG\left(VIX\right)$ models including all the regressors in Equations (1), (6), and (7). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: VIX			Panel B: LOG(VIX)
1 Day Ahead	HAR-SK	HAR-SK-PUTS	1 Day Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.231	0.274	HAR	0.302	0.169
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.387	HAR-SK		0.516
		(negative)			(negative)
5 Days Ahead	HAR-SK	HAR-SK-PUTS	5 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.086 *	0.016 **	HAR	0.475	0.150
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.584	HAR-SK		0.243
		(negative)			(negative)
10 Days Ahead	HAR-SK	HAR-SK-PUTS	10 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.067 *	0.001 ***	HAR	0.357	0.144
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.429	HAR-SK		0.152
		(negative)			(negative)
21 Days Ahead	HAR-SK	HAR-SK-PUTS	21 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.021 **	0.000 ***	HAR	0.0076 ***	0.003 ***
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.007 ***	HAR-SK		0.078 *
		(negative)			(negative)

and

Table 7. Giacomini-White’s type test for conditional predictive ability comparing $VIX$ and $LOG\left(VIX\right)$ models including all the regressors in Equations (1), (6), and (7). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: D(VIX)			Panel B: D(LOG(VIX))
1 Day Ahead	HAR-SK	HAR-SK-PUTS	1 Day Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.298	0.479	HAR	0.581	0.671
	(negative)	(negative)		(negative)	(positive)
HAR-SK		0.865	HAR-SK		0.651
		(negative)			(positive)
5 Days Ahead	HAR-SK	HAR-SK-PUTS	5 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.695	0.013 **	HAR	0.686	0.040 **
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.044 **	HAR-SK		0.279
		(negative)			(positive)
10 Days Ahead	HAR-SK	HAR-SK-PUTS	10 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.312	0.008 ***	HAR	0.121	0.482
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.061 *	HAR-SK		0.822
		(negative)			(negative)
21 Days Ahead	HAR-SK	HAR-SK-PUTS	21 Days Ahead	HAR-SK	HAR-SK-PUTS
HAR	0.702	0.007 ***	HAR	0.225	0.379
	(negative)	(negative)		(negative)	(negative)
HAR-SK		0.022 **	HAR-SK		0.696
		(negative)			(positive)

show the results of the Giacomini and White’s Giacomini and White [34] test of conditional predictive ability. This test differs from Hansen’s SPA test because it evaluates only models taken in pairs and because it is a conditional test rather than an unconditional one, i.e., the exploits the cross-serial correlation structure linking a set of instruments with the forecast errors generated by a predictive model. The loss function employed for this test is the squared error function but the results we have obtainted were essentially similar for an absolute value loss. In a GW test, a rejection of the null hypothesis due to a low p-value of the test statistic indicates that a model (say A) produces a larger average loss than another model (say B) if the test statistic is positive, and vice versa, model B produces a larger average loss than model A if the test statistic is negative. In the tables, models in the rows correspond to model A, and models in the columns correspond to model B when performing the test. If both null hypotheses fails to be rejected, the two models produce equally accurate forecasts.

The results of the test reported in Table 6 and Table 7 confirm the key conclusions already commented: when forecasting both levels and first differences of $VIX$ and log-$VIX$, respectively, the $HAR$-$SK$ and $HAR$-$SK$-$PUTS$ models never produce lower losses than the HAR model, and the $HAR$-$SK$ model is never significantly inferior to the $HAR$-$SK$-$PUTS$ model. The test statistic is non-zero (generally negative, thus favoring $HAR$-$SK$ over $HAR$-$SK$-$PUTS$) only when comparing these last two models at longer horizons and predicting the differenced series, but the null hypothesis is never rejected, at least in the case differencing is applied before forecasting.

Also in this case, When considering the results obtained by including only the significant predictors in the models, there are no notable discrepancies vs. the results detailed in Table 6 and Table 7.

5.2.2. Out-of-Sample Forecasting of the $SKEW$ Index

We now focus on evaluating the OOS forecasting power of the different models in predicting the $SKEW$ index. Since no model has been established in the literature that promises a particularly good predictive performance for the $SKEW$ index, we compare the predictions made using the three models in Equations (1), (8), and (9) with those of two additional, simple benchmarks: a Random Walk ($RW$) model and a $AR\left(1\right)$ model. (Even though Taylor [17] is applied to VIX futures returns (i.e., differences in the log-VIX) only, his fractionally integrated white noise model may represent a further, important benchmark, especially because Taylor shows that this model provides more accurate predictions of price changes and proportional returns, even when incorporating economic predictors like VIX futures price changes, S&P 500 index returns, interest rates, and the futures basis).

Table 8. SKEW and $LOG\left(SKEW\right)$ out-of-sample out-of-sample forecasting performance; models specified using all the regressors detailed in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: SKEW	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	9.41630	0.0000 ***	2.06954	0.0006 ***	42.674%
AR(1)	9.03937	0.0000 ***	2.07538	0.0000 ***	42.854%
HAR	8.32678	0.2972	2.00196	1.0000	42.395%
HAR-SK-PUTS	8.25872	1.0000	2.00742	0.3126	42.395%
HAR-IV	8.34485	0.1448	2.00652	0.0488	42.455%
5 Days Ahead
RW	22.4008	0.0000 ***	3.41638	0.0000 ***	48.870%
AR(1)	23.0523	0.0000 ***	3.54032	0.0000 ***	49.030%
HAR	13.6714	0.0174 **	2.68243	0.201	46.349%
HAR-SK-PUTS	13.0664	1.0000	2.66209	1.0000	47.590%
HAR-IV	13.7928	0.0072 ***	2.70624	0.0044 ***	46.509%
10 Days Ahead
RW	30.0114	0.0000 ***	4.02928	0.0000 ***	49.840%
AR(1)	33.0707	0.0000 ***	4.32587	0.0000 ***	49.780%
HAR	13.9164	0.0144 **	2.71174	0.2158	44.814%
HAR-SK-PUTS	13.2893	1.0000	2.69068	1.0000	46.656%
HAR-IV	14.0907	0.0066 ***	2.74054	0.0054 ***	45.334%
21 Days Ahead
RW	36.1040	0.0000 ***	4.52303	0.0000 ***	49.910%
AR(1)	41.4957	0.0000 ***	4.93136	0.0000 ***	50.110%
HAR	13.9770	0.0368 **	2.72108	0.2914	44.714%
HAR-SK-PUTS	13.4169	1.0000	2.70722	1.0000	46.596%
HAR-IV	14.2065	0.0148 **	2.75435	0.0022 ***	44.854%
Panel B: LOG(SKEW)	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	9.41630	0.0000 ***	2.06954	0.0002 ***	42.674%
AR(1)	9.04601	0.0000 ***	2.07234	0.0000 ***	42.854%
HAR	8.34038	0.1890	2.00136	1.0000	42.415%
HAR-SK-PUTS	8.25255	1.0000	2.00376	0.4378	42.375%
HAR-IV	8.34977	0.1504	2.00395	0.2728	42.534%
5 Days Ahead
RW	22.4008	0.0000 ***	3.41638	0.0000 ***	48.870%
AR(1)	23.0969	0.0000 ***	3.53040	0.0000 ***	49.030%
HAR	13.6335	0.0142 **	2.67487	0.2510	46.349%
HAR-SK-PUTS	12.9894	1.0000	2.65051	1.0000	47.389%
HAR-IV	13.6889	0.0040 ***	2.69241	0.0410 **	46.789%
10 Days Ahead
RW	30.0114	0.0000 ***	4.02928	0.0000 ***	49.840%
AR(1)	33.1226	0.0000 ***	4.31466	0.0000 ***	49.760%
HAR	13.8701	0.0202 **	2.70257	0.1684	44.794%
HAR-SK-PUTS	13.2060	1.0000	2.67817	1.0000	46.676%
HAR-IV	13.9511	0.009 ***	2.72385	0.0286 **	45.294%
21 Days Ahead
RW	36.1040	0.0000 ***	4.52303	0.0000 ***	49.910%
AR(1)	41.6049	0.0000 ***	4.92647	0.0000 ***	50.211%
HAR	13.9312	0.0434 **	2.71198	0.2452	45.053%
HAR-SK-PUTS	13.3310	1.0000	2.69433	1.0000	46.899%
HAR-IV	14.0464	0.0148 **	2.73713	0.0206 **	44.873%

and

Table 9. D(SKEW) and $D\left(LOG\right(SKEW\left)\right)$ out-of-sample out-of-sample forecasting performance; models specified using all the regressors detailed in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

Panel A: D(SKEW)	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	24.28220	0.000 ***	3.30213	0.0000 ***	43.802%
AR(1)	8.65950	0.000 ***	2.06680	0.1280	42.383%
HAR	8.34900	1.0000	1.99530	0.7776	42.603%
HAR-SK-PUTS	8.35399	0.4744	1.99649	0.4824	42.223%
HAR-IV	8.35871	0.3104	1.99495	1.0000	42.803%
5 Days Ahead
RW	19.43840	0.000 ***	3.06528	0.0000 ***	45.258%
AR(1)	9.41923	0.000 ***	2.06980	0.0000 ***	42.817%
HAR	8.60468	1.0000	2.00717	1.0000	42.717%
HAR-SK-PUTS	8.63325	0.1226	2.00964	0.2256	42.597%
HAR-IV	8.60680	0.6044	2.00926	0.1702	42.637%
10 Days Ahead
RW	18.68600	0.000 ***	3.01757	0.0000 ***	45.564%
AR(1)	9.42973	0.000 ***	2.07105	0.0000 ***	42.740%
HAR	8.60865	1.0000	2.00746	1.0000	42.620%
HAR-SK-PUTS	8.63791	1.0000	2.01005	0.2098	42.580%
HAR-IV	8.61016	0.6278	2.00892	0.2946	42.660%
21 Days Ahead
RW	18.67410	0.000 ***	3.01037	0.0000 ***	46.126%
AR(1)	9.44570	0.000 ***	2.07259	0.0000 ***	42.854%
HAR	8.62356	1.0000	2.00960	1.0000	42.714%
HAR-SK-PUTS	8.64947	0.1528	2.01178	0.2594	42.794%
HAR-IV	8.62377	0.6734	2.01055	0.4300	42.774%
Panel B: D(LOG(SKEW))	MSFE	SPA(SE)	MAFE	SPA(AE)	Correct Change of Direction
1 Day Ahead
RW	24.4788	0.000 ***	3.30586	0.0000 ***	43.862%
AR(1)	8.66328	0.000 ***	2.00726	0.0906 *	42.363%
HAR	8.34660	1.0000	1.99398	0.8326	42.463%
HAR-SK-PUTS	8.35007	0.5300	1.99496	0.5646	42.143%
HAR-IV	8.35830	0.2868	1.99373	0.0000 ***	42.743%
5 Days Ahead
RW	24.4932	0.000 ***	3.30680	0.0000 ***	43.838%
AR(1)	9.42282	0.000 ***	2.07033	0.0000 ***	42.857%
HAR	8.60552	1.0000	2.00607	1.0000	42.557%
HAR-SK-PUTS	8.62885	0.1538	2.00790	0.2780	42.757%
HAR-IV	8.60913	0.5668	2.00763	0.3088	42.577%
10 Days Ahead
RW	24.5129	0.000 ***	3.30849	0.0000 ***	43.861%
AR(1)	9.42973	0.000 ***	2.07105	0.0000 ***	42.740%
HAR	8.60958	1.0000	2.00633	1.0000	42.600%
HAR-SK-PUTS	8.63365	0.1522	2.00821	0.2742	42.740%
HAR-IV	8.61467	0.5172	2.00789	0.3068	42.339%
21 Days Ahead
RW	24.55210	0.000 ***	3.31047	0.0000 ***	43.938%
AR(1)	9.44570	0.000 ***	2.07259	0.0000 ***	42.854%
HAR	8.62369	1.0000	2.00836	1.0000	42.633%
HAR-SK-PUTS	8.64449	0.1862	2.00991	0.3358	42.834%
HAR-IV	8.62796	0.5448	2.00976	0.3370	42.493%

document the results obtained using all the predictors of our earlier in-sample analysis, in Section 5.1.2. From the tables, we can draw one conclusion: the benchmark $RW$ and $AR\left(1\right)$ are clearly outperformed by all other models in predicting all the specifications of the $SKEW$ at all the forecast horizons evaluated. This is shown by Hansen’s SPA test performed using both the squared and the absolute loss functions. The only exception is the $AR\left(1\right)$ model, which is not clearly inferior to other models in predicting the first difference of the $SKEW$ and the first difference of the logarithm of the $SKEW$ at the 1-day horizon under the absolute loss.

Panel A in Table 8 displays the results obtained in predicting the levels of the $SKEW$. According to the SPA (SE) test, the $HAR$-$SK$-$PUTS$ model outperforms the benchmarks at predicting the $SKEW$ index at all horizons except at the 1-day according to the SPA (AE) test. Yet, there is no clear superior model between the $HAR$ and the $HAR$-$SK$-$PUTS$. The same results can be observed from Panel B in Table 8 when forecasting the $SKEW$ via its logarithm.

Again, as we did when evaluating the out-of-sample forecasts of the $VIX$ index, we also evaluate our predictions using the first differences of both the $SKEW$ and of its logarithm. Panels A and B of Table 9 report the corresponding results. In this case, there is no clear superior model among the three $HAR$ frameworks, neither in terms of squared errors nor in terms of absolute errors, even though they all outperform the benchmarks.

Unlike the case with the $VIX$ index in Section 5.2.1, the percentage of correct predicted direction of change Table 8 and Table 9, offers some useful insights. At the 1-day horizon, this percentage is approximately 42% for all the models employed, including the benchmark models. Paradoxically, this suggests that establishing a trading strategy based on the opposite of the direction predicted by the models might yield superior returns. However, at longer horizons, the percentage obtained by the benchmark models hovers around 50%, indicating that no clear indication about the direction of the $SKEW$ index can be extracted using either a $RW$ or an $AR\left(1\right)$ model. On the other hand, employing the $HAR$ model and its other two specifications yields a percentage ranging from 44% to 46%, which reveals that while the benchmarks may face problems when forecasting the $SKEW$ as whole, matters are more complicated when it comes to predict the sign of the change in the $SKEW$. These conclusions apply to both the levels of the $SKEW$ index (Panel A of Table 8) and to its logarithm (Panel B of Table 8). When the forecasts are instead to be based on the first difference of the $SKEW$ series, the percentage of correct sign predictions remains around 42–43% at all horizons for all the alternative models and for the $AR\left(1\right)$ model. The only exception appears to be the RW model, which ranges from 43% to 46%.

Considering only those predictors that were found to be significant when analyzing the in-sample properties of the time series of the $SKEW$ index in Section 5.1.2, we obtain additional evidence of the robustness of our results. The $HAR$-$SK$-$PUTS$ model is superior to the others according to Hansen’s test using the squared errors loss function at all forecast horizons except the 1-day horizon. Nonetheless, when we assess predictive performances under the absolute error loss function, there is no clear superior model between the $HAR$ and the $HAR$-$SK$-$PUTS$. So results remain mixed when it comes to compare different strands of $HAR$ models and in terms of realized signed change prediction performance, even though there is never any doubt that $HAR$-type models are superior to standard benchmarks.

Table 10. Giacomini-White’s type test for conditional predictive ability comparing $SKEW$ including all the regressors in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	1 Day Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.003 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)		0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR			0.000 *** (positive)	0.165 (negative)
HAR-SK-PUTS				0.004 *** (negative)
	5 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.439 (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)		0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR			0.016 ** (positive)	0.036 ** (negative)
HAR-SK-PUTS				0.004 *** (negative)
	10 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.114 (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)		0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR			0.002 *** (positive)	0.022 ** (negative)
HAR-SK-PUTS				0.001 *** (negative)
	21 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.088 * (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)		0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR			0.003 *** (positive)	0.051 * (negative)
HAR-SK-PUTS				0.003 *** (negative)

,

Table 11. Giacomini-White’s type test for conditional predictive ability comparing $LOG\left(SKEW\right)$ including all the regressors in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	1 Day Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.006 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.000 *** (positive)	0.000 *** (positive)	0.543 (negative)
HAR-SK-PUTS				0.015 ** (negative)
	5 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.388 (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.008 *** (positive)	0.000 *** (positive)	0.154 (negative)
HAR-SK-PUTS				0.004 *** (negative)
	10 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.095 * (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.001 *** (positive)	0.000 *** (positive)	0.132 (negative)
HAR-SK-PUTS				0.002 *** (negative)
	21 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.074 * (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.003 *** (positive)	0.000 *** (positive)	0.083 * (negative)
HAR-SK-PUTS				0.003 *** (negative)

,

Table 12. Giacomini-White’s type test for conditional predictive ability comparing $D\left(SKEW\right)$ including all the regressors in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	1 Day Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.137 (negative)	0.000 *** (positive)	0.600 (negative)
HAR-SK-PUTS				0.264 (negative)
	5 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.351 (negative)	0.000 *** (positive)	0.011 ** (negative)
HAR-SK-PUTS				0.462 (positive)
	10 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.086 (negative)	0.000 *** (positive)	0.402 (negative)
HAR-SK-PUTS				0.224 (positive)
	21 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.211 (negative)	0.000 *** (positive)	0.790 (negative)
HAR-SK-PUTS				0.272 (positive)

and

Table 13. Giacomini-White’s type test for conditional predictive ability comparing $D\left(LOG\right(SKEW\left)\right)$ including all the regressors in Equations (1), (8), and (9). In the table, * means statistical significance at 10%, ** and *** indicate statistical significance at 5 and 1%, respectively.

	1 Day Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.155 (negative)	0.000 *** (positive)	0.460 (negative)
HAR-SK-PUTS				0.239 (negative)
	5 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.439 (negative)	0.000 *** (positive)	0.163 (negative)
HAR-SK-PUTS				0.575 (positive)
	10 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.083 (negative)	0.000 *** (positive)	0.431 (negative)
HAR-SK-PUTS				0.331 (positive)
	21 Days Ahead
	AR(1)	HAR	HAR-SK-PUTS	HAR-IV
RW	0.000 *** (negative)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
AR(1)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)	0.000 *** (positive)
HAR		0.306 (negative)	0.000 *** (positive)	0.727 (negative)
HAR-SK-PUTS				0.440 (positive)

report the results of the GW test Giacomini and White [34] for conditional forecasting ability in the case of $SKEW$ predictions. As in the case of the $VIX$, the models in the rows correspond to model A and the models in columns correspond to model B when performing the test. A clear finding is that also in the perspective of a GW test, all the $HAR$ models consistently produce superior forecasts compared to the two benchmarks, as indicated by the rejection of the null hypothesis and by the positive test statistics when comparing the benchmarks vs. $HAR$ models. Furthermore, when comparing the performances of the $HAR$ model with the $HAR$-$SK$-$PUTS$ model, the latter produces superior forecasts at all horizons in predicting the $SKEW$ index and its logarithm (Table 10 and Table 11). However, this conclusion cannot be confirmed when we forecast the first difference of the $SKEW$ and the logarithm of the $SKEW$ (see Table 12 and Table 13). Another interesting result can be inferred from Table 10 and Table 11 is that the $HAR$-$IV$ model is found to be always inferior to the $HAR$-$SK$-$PUTS$ model in predicting the $SKEW$ index and of its logarithm. Once more, this conclusion does not hold when assessing the two models in terms of their ability to predict the first difference of the series.

When we investigate the models that include only the predictors found to be statistically significant over the in-sample estimation of the $HAR$ model (the $RW$ and $AR\left(1\right)$ benchmark models specifications and realized predictive performaces do not change because these models contain all predictors by construction) in Section 5.1.2, we cannot find any relevant difference in the OOS result with respect to what concluded above.

Finally, by comparing the conclusions derived from Hansen [31]’s SPA test with those from the GW test Giacomini and White [34], we find that the two types of predictive accuracy tests yield consistent results. The only difference arises at the 1-day horizon, where the Hansen’s test cannot detect any superior model among the three different specifications of the $HAR$ model, while GW’s test does.

6. Model Confidence Set Results

As a final, conclusive step, we apply to our problem the Model Confidence Set (henceforth, MCS) methodology by Hansen et al. [35]. Our question becomes then not really testing whether a plain vanilla $HAR$ may represent the best performing model in terms of OOS accuracy, but whether it is possible to identify a subset of models that are indistinguishable in terms of predictive accuracy. Therefore, instead of selecting a single “best” model, MCS recognizes that multiple models can perform similarly. The method begins with a set of candidate models and iteratively removes those that are statistically worse than others. It uses a test statistic, essentially a t-statistic in our implementation, to compare predictive performance differences. Models with higher test statistics are eliminated, and the process continues until only models that are not significantly different remain. To ensure robustness and handle sample variability, block bootstrap techniques are used to calculate p-values and assess the significance of the differences. MCS provides a set of reliable models, avoiding premature selection of a model that may appear better due to sample noise.

In particular, in our implementation, we use 20,000 bootstraps and select heterogeneous block lengths based on the minimization of the BIC information criterion. The threshold to determine which models are removed from the MCS is set at 0.05. The test statistic is the Selection Quantile (SQ) one, the first one proposed in Hansen et al. [35].

Table 14. Model Confidence Set Results for Alternative Models and Horizon—Squared and Absolute Loss Functions for VIX and SKEW. Boldfaced p-values indicate models included in the MCS.

	Squared Loss Function					Absolute Loss Function
	HAR	HAR-SK	HAR-SK-PUT	RW	AR(1)	HAR	HAR-IV	HAR-SK-PUT	RW	AR(1)
Panel A: VIX
Horizon = 1-day
VIX in levels	1.000	0.000	0.000	0.000	0.000	1.000	0.125	0.000	0.125	0.125
$\Delta$VIX	0.143	0.000	0.000	0.000	1.000	1.000	0.000	0.000	0.143	0.000
log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
$\Delta$log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
Horizon = 5-day
VIX in levels	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000	1.000
$\Delta$VIX	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
$\Delta$log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
Horizon = 10-day
VIX in levels	1.000	0.000	0.000	0.000	0.000	1.000	0.125	0.000	0.000	0.000
$\Delta$VIX	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
log(VIX)	1.000	0.125	0.125	0.000	0.000	1.000	0.125	0.125	0.000	0.000
$\Delta$log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
Horizon = 21-day
VIX in levels	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
$\Delta$VIX	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
$\Delta$log(VIX)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
	HAR	HAR-SK	HAR-SK-PUT	RW	AR(1)	HAR	HAR-IV	HAR-SK-PUT	RW	AR(1)
Panel B: SKEW
Horizon = 1-day
SKEW in levels	0.100	0.000	1.000	0.000	0.000	1.000	0.000	0.100	0.000	0.000
$\Delta$SKEW	1.000	0.000	0.000	0.000	0.000	0.111	1.000	0.000	0.000	0.000
log(SKEW)	0.100	0.000	1.000	0.000	0.000	1.000	0.000	0.100	0.000	0.000
$\Delta$log(SKEW)	1.000	0.111	0.111	0.000	0.000	0.111	1.000	0.000	0.000	0.000
Horizon = 5-day
SKEW in levels	0.000	0.000	1.000	0.000	0.000	0.100	0.000	1.000	0.000	0.000
$\Delta$SKEW	1.000	0.111	0.111	0.000	0.000	1.000	0.000	0.111	0.000	0.000
log(SKEW)	0.000	0.000	1.000	0.000	0.000	0.100	0.000	1.000	0.000	0.000
$\Delta$log(SKEW)	1.000	0.000	0.000	0.000	0.000	1.000	0.111	0.111	0.000	0.000
Horizon = 10-day
SKEW in levels	0.000	0.000	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000
$\Delta$SKEW	1.000	0.111	0.000	0.000	0.000	1.000	0.000	0.000	0.000	0.000
log(SKEW)	0.000	0.000	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.000
$\Delta$log(SKEW)	1.000	0.000	0.000	0.000	0.000	1.000	0.000	0.111	0.000	0.000
Horizon = 21-day
SKEW in levels	0.000	0.000	1.000	0.000	0.000	0.100	0.000	1.000	0.000	0.000
$\Delta$SKEW	1.000	0.111	0.000	0.000	0.000	1.000	0.111	0.000	0.000	0.000
log(SKEW)	0.000	0.000	1.000	0.000	0.000	0.100	0.000	1.000	0.000	0.000
$\Delta$log(SKEW)	1.000	0.111	0.000	0.000	0.000	1.000	0.111	0.000	0.000	0.000

presents the empirical results for both the $VIX$ and the $SKEW$ in their four alternative definitions (levels, differences, logs, and differences in logs). The table lists empirical findings for both the squared and the absolute value loss functions and in correspondence to each model we report the SQ statistic p-value. The p-values that are boldfaced refer only to the model(s) that end up belonging to the MCS.

Table 14 documents the main results of the MCS tests and shows that out of a total of 64 instances in which the MCS was applied to the forecast errors derived from the three $HAR$-style models ($HAR$, $HAR$-$SK$/$IV$, and $HAR$-$SKEW$-$PUT$) and the two benchmarks ($RW$ and $AR\left(1\right)$), in 54 cases the plain vanilla $HAR$ model turns out to belong to the confidence set when a 5% p-value threshold is imposed. This is all the more impressive because out of 64 sets estimated, only in 38 cases the MCS includes only a model and of such 38 instances, in 29 of them it is the plain vanilla $HAR$ model that represents the singleton. This means that a simple $HAR$ framework is indeed very hard to dismiss from the set of forecasts a researcher ought to entertain. For instance, the first row of the table shows that when forecasting the $VIX$ level OOS, under a squared loss, the MCS contains only the plain vanilla $HAR$; under an absolute value loss, the MCS contains instead not only the simple $HAR$, but also $HAR$-$SK$ and the two benchmarks, albeit the latter carry rather modest p-values of 0.125 (yet exceeding 0.05). This conclusion does not depend on the assumed loss function and applies with particular strength to the $VIX$ in its various transformations (30 cases out of 32).

Interestingly, and this is especially the case when we forecast the $SKEW$ (in 21 cases out of 32), also the $HAR$ model that includes put-based only skewness measures tends to belong to the MCS, in total in 23 cases out of 64. The two benchmarks are featured in the MCS very rarely: even when we combine both $RW$ and the $AR\left(1\right)$ model, this occurs only in 5 cases out of 64, which may almost be attributed to pure chance (even though such an application of the concept of size of a test may not be a straightforward in the case of the MCS).

7. Conclusions

We have assessed the OOS forecasting power of simple $HAR$ models to forecast (a number of transformations of) the CBOE $VIX$ and $SKEW$ indices. Besides ranking the realized, OOS performance of the standard $HAR$ and of variations of $HAR$ that allow lagged $SKEW$ to forecast the $VIX$ and viceversa, we use Hansen [31]’s Superior Predictive Accuracy, Giacomini and White [34]’s, and Hansen et al. [35]’s model confidence set methods to test whether the (unconditional and conditional, respectively) differential in predictive accuracy are statistically significant or whether any model can be removed from the set that is worthwhile to consider in forecasting applications. We devote special attention to the idea that $HAR$ models for the $VIX$ ($SKEW$) that also include lagged (transformations of) the $SKEW$ ($VIX$) may provide superior predictive performance vs. standard $HAR$ models and that special definitions of the $SKEW$ index (based only on put options or just measuring downside risk) may carry special information content that is valuable in OOS applications. Nonetheless, we conclude that—especially in $VIX$ forecasting applications—a plan vanilla $HAR$ model delivers the best predictive performance and that only in specific cases involving specific transformations of the $SKEW$ may provide a role to $HAR$ models that are augmented to include additional option-implied information. Because all $HAR$ models can outperform simple but common benchmarks in the literature (such as the random walk and AR(1) models), we lean towards the idea that it is the robust predictive power of $HAR$ that drives our results. This finding is confirmed by the MCS results in which in almost 85% of the combinations of target forecast variables, loss functions, and horizons, a plain vanilla $HAR$ turns out to be impossible to remove from the set to be taken into account. Therefore applications to $VIX$ and $SKEW$ forecasting confirm once more a classical results of applications of $HAR$ modelling to realized, physical volatility, i.e., that systematically outperforming the standard framework popularized by Corsi [6] remains difficult in most cases.

Our objective was not to claim superiority of the $HAR$ approach but to stress that a considerable predictive mileage can also be obtained from a very simple approach. Obviously, while $HAR$ is justified for forecasting the $VIX$ and $SKEW$ indices due to its ability to capture long-memory properties and the impact of volatility at different time scales but remains intuitive and simple, wavelet modeling is more complex, focusing on formally decomposing time series into different frequency components. While wavelets offer insights into localized variations at multiple resolutions, they can be computationally intensive and require careful tuning of parameters and remain an exciting avenue for future research.