A machine learning approach to volatility forecasting

A machine learning approach to volatility forecasting

\thispagestyle{empty}

Introduction

\setcounter{page}{1} The unprecedented levels of volatility in financial markets during the last decades has led to increased awareness about the importance of explaining key drivers behind such movements. On the one hand, the structure of financial markets is complex, highly nonlinear and has a low signal-to-noise ratio, which makes it nearly impossible to predict short-term asset returns \citep*[see, e.g.,][]{gu-kelly-xiu:20a, chen-pelger-zhu:24a}. On the other hand, various stylized facts (e.g., the slow decay of autocorrelation in absolute returns) implies that volatility is to a large extent predictable.

With its crucial role in asset pricing, portfolio allocation, and risk management, volatility forecasting has attracted great attention. The origins of this literature go back to \citet*{engle:82a, bollerslev:86a, taylor:82a}, who developed discrete-time GARCH and stochastic volatility models for autoregressive conditional heteroskedasticity. As argued by \citet*{corsi:09a}, however, standard volatility models are not able to reproduce several salient features of financial markets. He proposes a Heterogeneous Auto-Regressive (HAR) model, which combines nonparametric realized variance measured at different frequencies with a parametric autoregressive model, in order to capture heterogeneous types of market participants and approximate long-memory of volatility. This has been an incredibly successful model, which is by now “the benchmark,” to a large extent occupying the former role of the GARCH(1,1) \citep*[e.g.][]{hansen-lunde:05b}. However, because the HAR model is so parsimonious, the complicated structure of the underlying data renders it inadequate in some directions. Therefore, a number of extensions of the baseline HAR have been proposed \citep*[e.g.,][]{andersen-bollerslev-diebold:07a, corsi-reno:12a, patton-sheppard:15a, bollerslev-patton-quaedvlieg:16a}. These are, however, still based on past returns as conditioning information, which is surprising, since a large amount of research has documented an intimate relationship between news and volatility \citep*[e.g.,][]{schwert:89a, engle-ghysels-sohn:13a, bollerslev-li-xue:18a}.

The reason that no or only few additional covariates are typically included in volatility prediction is partly explained by the observation that past volatility is a very powerful predictor of its future self. That is, after all, how this literature evolved. However, it is also motivated by the fact that traditional models, often estimated with linear regression, break down when the explanatory variables are strongly correlated, exhibit low signal-to-noise ratios, or if the underlying structure is nonlinear. As explained above, this is a trademark of financial markets, so it is unsurprising that relatively few studies conduct volatility forecasting in a data-rich environment. The recent access to large datasets therefore highlights a need for investigating more powerful tools that take these problems into account.

As emphasized by \citet*{varian:14a}, machine learning (ML) techniques overcome these deficiencies. Up to now, however, only a few studies combine ML and volatility forecasting using a single approach, e.g., \citet*{audrino-knaus:16a, audrino-sigrist-ballinari:20a, caporin-poli:17a} investigate lasso, \citet*{luong-dokuchaev:18a} use random forests, \citet*{mittnik-robinzonov-spindler:15a} apply component-wise gradient boosting\glsadd{Boosting}, \citet*{bucci:20a, donaldson-kamstra:97a, hillebrand-medeiros:10a, fernandes-medeiros-scharth:14a, rahimikia-poon:20a} explore neural networks. In related work, \citet*{carr-wu-zhang:20a} study the application of ML to low-frequency options data in order to forecast 30-day realized variance of the S&P 500 index.

In contrast, we conduct a comprehensive analysis of the out-of-sample performance of multiple ML techniques for volatility forecasting, i.e. linear regularized models\glsadd{Regularization} (ridge, lasso, elastic net), tree-based algorithms (bagging, random forest, gradient boosting), and neural networks, when a reasonably large number of additional controls appear in the information set. In such a setting, regularization is helpful to deal with the high-dimensionality of the problem and get less noisy parameter estimates. Furthermore, tree-based regression allows for estimation without imposing any distributional or functional form on the model, while neural networks can regularize and capture nonlinearities. We compare these rather different ML approaches to the HAR lineage. We aim not only to do a comparison of HAR models and ML methods but also to provide evidence of how and why ML improves the accuracy of forecasting volatility. We adopt a conservative approach to implement ML in that we do minimal hyperparameter tuning to assess how good the techniques are when implemented “as is”.

The findings are five-fold. Firstly, we show superior out-of-sample forecasting with off-the-shelf implementations of ML compared to HAR. A \citet*{diebold-mariano:95a} test for equal predictive accuracy shows that these gains are often viewed as statistically significant. Secondly, we document that substantial incremental information about future volatility can be extracted with ML from additional volatility predictors. In contrast, the HAR models yield only minor improvements or outright deteriorate in forecast accuracy, unless regularization is applied to account for overfitting. Thirdly, we compute variable importance via accumulated local effects (ALE), which shows a general agreement about the set of the most dominant predictors, but disagreement on their ranking. The ALE also reveal a nonlinear interaction between covariates. Fourthly, we show that the forecast gains are stronger at longer horizons, which we attribute to ML techniques better capturing an underlying long-memory behavior in the realized variance. Fifthly, we conduct a Value-at-Risk application. We show that ML continues to excel at forecasting in this setting, although the differences are smaller and less statistically significant. Moreover, ML produces correct unconditional and conditional coverage.

The paper is organized as follows. Section (ref) sets the theoretical foundation and introduces the various models applied in the paper. Section (ref) describes the high-frequency data used in the empirical analysis. Section (ref) features an out-of-sample one-day-ahead forecast comparison and conducts a series of robustness checks. In Section (ref), we examine the weekly and monthly horizon. In Section (ref), we perform an economic application by forecasting Value-at-Risk. We conclude in Section (ref). In the Appendix, we present supplemental information, including a glossary that can assist readers who are not acquainted with the jargon of ML.

Methodology

The setting

We assume the log-price $X = \left(X_{t} \right)_{t \geq0}$ is supported by a filtered probability space $\big( \Omega, ( \mathcal{F}_{t})_{t \geq 0}, \mathcal{F}, \mathbb{P} \big)$. If the price is determined in an arbitrage-free frictionless market, then $X$ is a semimartingale process, see \citet*{delbaen-schachermayer:94a}. We assume $X_{t}$ can be represented as follows:

where $X_{0}$ is $\mathcal{F}_{0}$-measurable, $\mu = (\mu_{t})_{t \geq0}$ denotes the drift, $\sigma = ( \sigma)_{t \geq 0}$ is the stochastic volatility process, $W = (W_t)_{t \geq 0}$ is a standard Brownian motion, $N = (N_{t})_{t \geq 0}$ is a counting process, which represents the number of jumps in $X$, and $J = (J_{s})_{s=1, \dots,N_{t}}$ is a sequence of nonzero random variables of jump sizes with corresponding jump times $\tau = ( \tau_{s})_{s=1, \dots,N_{t}}$.

Our aim is to forecast the daily quadratic variation, which is $QV_{t} = \int_{t-1}^{t} \sigma_{s}^{2}ds + \sum_{t-1 \leq \tau_{s} \leq t} J_{s}^{2}$, for $t = 1, \ldots, T$, where $t$ is the day and $T$ is the total number of days in the sample. The quadratic variation is not directly observable, since in practice only discretely sampled measurements of $X$ are available. An estimator of QV is given by the realized variance:

where $n$ is the number of intraday logarithmic returns, $\Delta_{t-1,j}^{n}X = X_{t-1+\frac{j}{n}} - X_{t-1 + \frac{j-1}{n}}$. Then, as shown by \citet*{andersen-bollerslev:98a, barndorff-nielsen-shephard:02a}, realized variance is a consistent estimator of the quadratic variation as the sampling frequency increases, i.e. $RV_{t} \overset{ \mathbb{P}}{ \longrightarrow} QV_{t}$ as $n \rightarrow \infty$.

Later on, we also employ a decomposition, where realized variance is broken into its positive and negative semivariance part \citep*[e.g.][]{barndorff-nielsen-kinnebrock-shephard:10a}:

where

The HAR model and its extensions

The parsimonious structure of the HAR model of \citet*{corsi:09a} makes it a default choice in the literature. Therefore, the HAR model is the benchmark throughout this study, ensuring comparability with earlier work.

The HAR is defined as:

where $RV_{t-1 \mid t-h}=\frac{1}{h} \sum_{i=1}^{h}RV_{t-i}$. The explanatory variables are the daily, weekly, and monthly averages of lagged realized variance.

After its inception, a large number of extensions of the basic HAR has been developed. To ensure our results are robust against a broader suite of recent HAR-type models, we include several of the more popular ones as “horses” in the race. In particular, to impose a nonlinear relationship between present and future volatility and attenuate the influence of outliers, \citet*{corsi:09a} also proposed the log-version of the HAR (LogHAR):

To capture a leverage effect, \citet*{corsi-reno:12a} suggested to include past aggregated negative returns into the HAR model, which is dubbed LevHAR:

where

and $r_{t-1 \mid t-h}^{-} = \min(0, r_{t-1 \mid t-h})$.

An alternative way to model an asymmetric impact of past positive and negative returns on future volatility is given by the semivariance HAR (SHAR) of \citet*{patton-sheppard:15a}:

We also add the HARQ by \citet*{bollerslev-patton-quaedvlieg:16a}, which corrects an inherent error-in-variables problem in the HAR that arises because realized variance is a generated regressor. The HARQ is defined as:

where

is the realized quarticity.

As readily seen, the formula for a general HAR model (except the LogHAR) is:

where $Z_{t-1}$ is the information set, $\beta = ( \beta_{1}, \ldots, \beta_{J})'$ are the slope parameters, and $J$ is the number of explanatory variables. The HAR model has $Z_{t-1} = (RV_{t-1}, RV_{t-1 \mid t-5}, RV_{t-1 \mid t-22})$, whereas the HARQ has $Z_{t-1} = (RV_{t-1}, RQ_{t-1}^{1/2}RV_{t-1}, RV_{t-1 \mid t-5}, RV_{t-1 \mid t-22})$.

In general, $Z_{t-1}$ can include covariates that are not functions of past returns. We inspect such a “big data” setting. A complete list of extra explanatory variables included in this paper is available in Table (ref). To ensure comparability with ML, we therefore adopt a distributed lag-type version of each HAR model to allow for a broader selection of variables to be included as conditioning information. We denote by HAR-X the extended version of the basic HAR. However, to keep notation down the other HAR models are still denoted with their original label even in the extended setting.

While the above list is comprehensive, it is necessarily incomplete. We do not include HAR models, where the quadratic variation is split into a continuous and jump part \citep*[e.g.][]{andersen-bollerslev-diebold:07a}.

In order to estimate the coefficients of the HAR models, we employ a least squares approach with objective function:

where $f(Z_{t}; \beta) = \beta'{Z}_{t}$ and in-sample is defined later.

Regularization

When increasing the number of predictors in a setting with a low signal-to-noise ratio, the linear models at some point start to fit noise rather than relevant information, also known as overfitting. A common way to avoid overfitting and increase out-of-sample performance is to shrink\glsadd{Shrinkage} the regression coefficients by imposing a penalty term.

We define the penalized loss function:

where $\phi( \beta; \theta)$ is a penalty term, and $\theta$ is a vector of hyperparameters\glsadd{Hyperparameters}, which are always determined in the validation set\glsadd{Validation Set}.\footnote{A thorough explanation of the validation procedure can be found in Appendix (ref).}

In this paper, we estimate three of the most widely used penalization methods, i.e. ridge regression (RR) by \citet*{hoerl-kennard:70a} and lasso (LA) by \citet*{tibshirani:96a}, where the penalty term is an $L^{2}$ and $L^{1}$ norm, and elastic net (EN) by \citet*{zou-hastie:05a}, which is a convex combination of ridge and lasso.

In ridge regression, the penalty term is given by:

where $\lambda \geq 0$ controls the amount of shrinkage. In the implementation, we choose a range for $\lambda$ that is effectively broad enough to include both the unregularized model and the constant model as potential solutions.

If dealing with a large feature space\glsadd{Feature}, shrinking is not always sufficient, and subset selection\glsadd{Subset selection} is preferable. Lasso addresses this with a penalty term of the form:

Due to the geometry, lasso often forces coefficient estimates to zero and generates sparsity\glsadd{Sparsity} and subset selection. The cost of this is that there is no closed-form solution to the problem, which makes it necessary to conduct numerical optimization.

The final penalization approach is the elastic net:

where $\alpha \in [0,1]$.\footnote{The elastic net nests lasso for $\alpha = 0$ and ridge for $\alpha = 1$.}

To gauge the effect of subset selection as a tool to pre-filter the data, we also include a post lasso HAR (P-LA), in which the penalization in (ref) is applied to the HAR model in a first-stage regression, and then only the variables with a non-zero slope coefficient estimate are included in an unrestricted second-stage HAR regression.

At last, following \citet*{zou:06a} we also estimate an adaptive lasso HAR model (A-LA). Here, the two-stage procedure of the post lasso are basically reversed. Firstly, an unrestricted HAR is estimated to get estimates of the slope parameters, say $\hat{ \beta}^{\text{1st}}$. Secondly, a lasso is run with a modified penalty term linked to the reciprocal of the elements in $\hat{ \beta}^{\text{1st}}$:

where $\hat{w}_{i}= 1/ | \hat{\beta_{i}}^{\text{1st}}|$. Hence, a larger magnitude of the beta coefficient from the first stage leads to it being downweighted in the penalization during the second stage.\footnote{All variables are standardized prior to estimation to make the slope parameters scale invariant.}

Tree-based regression

A concern with linear models is that it is left to the researcher to impose the true association among predictors and the response variable. In contrast, a regression tree is a fully nonparametric approach, which not only allows for nonlinearity, but it also implicitly accounts for interaction effects between explanatory variables.

A tree-based regression is, as the name suggests, based on a decision tree. A tree is grown by partitioning the domain (or feature space in ML parlance)\glsadd{Feature} of the explanatory variables $Z_{t}$, $\text{dom}(Z_{t}) \subseteq \mathbb{R}^{J}$, into smaller and smaller rectangular subspaces as we move through the tree (an illustrative example is provided in Figure (ref)). This process is continued until a stopping criterion is reached. At the end of the tree, a constant prediction is assigned to all observations that fall in a given terminal node (also called a leaf).\glsadd{Node Regression tree}

Suppose that after constructing the tree there are $M$ terminal nodes. They correspond to a sequence of disjoint rectangles $R_{m}$, for $m = 1, \ldots, M$, such that $R_{m} \subseteq \text{dom}(Z_{t})$ and $\bigcup_{m=1}^{M} R_{m} = \text{dom}(Z_{t})$. The tree-based regression then predicts as follows:

where $\widehat{RV}_{m}$ is a constant.

To ensure internal consistency, we adopt the minimum sum of squares criterion for selecting $\widehat{RV}_{m}$, i.e. $\widehat{RV}_{m} = \text{average}(RV_{t+1} \mid Z_{t} \in R_{m})$.

In essence, tree-based regression poses a combinatorial problem. There is no closed-form solution to find the best order in which to sort and split variables, and it may be too time-consuming to search over all possible choices. Hence, the need for a short-sighted, greedy algorithm arises. A popular approach to this problem is to start at the top node (called the root) and independently optimize nearby nodes in a stepwise fashion, known as the Classification And Regression Tree (CART). The CART, however, often yields inferior out-of-sample performance due to high variance from overfitting sample noise.

\citet*{breiman:96a} proposed bagging (bootstrap aggregation, BG) to refine the decision tree. In bagging, we grow several trees by resampling the original data. We subsequently predict as follows:

where $\hat{f}^{b}(Z_{t})$ is the tree-based prediction from the $b$th bootstrap sample.

The serial correlation in volatility suggests a block bootstrap may be appropriate. However, in our experience an i.i.d. bootstrap does not impair the recursive segmentation of the input space compared to more refined resampling.\footnote{As a verification, we applied a block bootstrap with $B$ random samples from a given dataset $\mathcal{M}$, resulting in a sequence $\mathcal{M}^{b}$, for $b = 1, \ldots, B$, where the block length $l$ and the number of blocks in each bootstrap sample $k$ was chosen such that $T = kl$. $l$ was found in the validation set\glsadd{Validation Set}. We estimated a regression tree for each bootstrap and computed the prediction function as $\hat{f}(Z_{t}) = \frac{1}{B} \sum_{b=1}^{B} \hat{f}^{b}(Z_{t})$. There was no discernible change in the results compared to the i.i.d. bootstrap.}

The benefit of bagging compared to CART is related to the reduced correlation between trees. Therefore, an alternative approach building on bagging is random forest (RF), see \citet*{breiman:01a}. Here, attention is restricted to a random subset of input features before finding the best split point. In this way a random forest is able to de-correlate the trees further. In particular, $B$ trees $\big\{ f(Z^{b}_{t}, \theta^{b}) \big\}_{b=1}^{B}$ are grown, where $\theta^{b}$ summarizes the $b$th random forest tree in terms of split variables, split points, and values at the terminal nodes. The prediction is then:

To reduce the complexity of the random forest and bagging, we set all tuning parameters equal to their default values from the original Fortran code by \citet*{breiman-cutler:04a}. This sidesteps validation and facilitates reproduction of our findings.

Another tree-based approach is gradient boosting (GB)\glsadd{Boosting} by \citet*{friedman:01a}.\footnote{He extends the original boosting algorithm from \citet*{schapire:90a} and \citet*{freund:95a}.} Gradient boosting produces a sequential model based on weak learners\glsadd{Weak learner}, where each of $B$ trees is grown on information from the previous one. Firstly, $\hat{f}^{0}(Z_{t})$ is initialized as a constant determined by $\hat{f}^{0}(Z_{t}) = \operatorname{arg\,min}_{ \beta_{0}} \sum_{t \in \text{in-sample}} (RV_{t} - \beta_{0})^{2}$. With the MSE loss function, this constant is the sample average. Secondly, the negative gradient of the loss function with respect to the prediction is calculated. With the MSE loss function, this corresponds to the residuals. A shallow tree is then fitted to the residuals, yielding a set of terminal nodes $R_{jb}$, $j = 1, \ldots, J_{b}$, where $j$ denotes the leaf and $b$ the tree. This is followed by choosing a gradient descent size as $\rho^{jb} = \operatorname{arg\,min}_{ \beta_{0}} \sum_{Z_{t} \in R_{jb}} (RV_{t} - \beta_{0} - \hat{f}^{b-1}(Z_{t}))^{2}$, and in the last step $\hat{f}^b(Z_{t})$ is updated iteratively as follows:

for $b = 1, \ldots, B$, where $\nu$ is a learning rate\glsadd{Learning Rate} and $\hat{f}(Z_{t}) \equiv \hat{f}^{B}(Z_{t})$ the final prediction.

As explained by \citet*{deprado:18a}, gradient boosting attempts to address an underfitting problem, whereas random forest deals with overfitting (generally perceived to be a greater evil in finance). As a consequence, gradient boosting places substantially more weight on misclassification and is therefore susceptible to outliers. With the typical heavy tails in the volatility distribution, we expect this to lead to poor results for gradient boosting in our framework.

Neural network

To ensure a complete investigation of ML space, at last we apply the artificial neural network (NN). With their highly flexible and nonlinear structure, neural networks show astonishing performance in handling complex problems. However, with a broad set of hyperparameters\glsadd{Hyperparameters} the replicability of the neural network is frequently questioned. Therefore, we construct the neural network as streamlined as possible and set the majority of the hyperparameters\glsadd{Hyperparameters} to standard suggestions from the literature. This is presented in Appendix (ref).

We briefly outline the mechanics of a so-called feed-forward network, while a more comprehensive overview is available in, e.g., \citet*{goodfellow-bengio-courville:16a}. A neural network is constructed from an input layer\glsadd{Input layer}, hidden layers\glsadd{Hidden layer}, and an output layer\glsadd{Output layer}. Suppose the total number of such layers is $L$. In the first layer, the neural network receives an input $Z_{t}$. The data is transformed through one or more hidden layers using an activation function\glsadd{Activation function} $g$.

The general equation for the $l$-th layer is denoted as:

where $g_{l}$ is the activation function, $\theta^{(l)}$ is the weight matrix, $b^{(l)}$ is the error serving as an activation threshold for the neurons\glsadd{Nodes Neural network} in the next layer, $N_{l}$ denotes the number of hidden neurons, and $a_{t}^{ \theta_{l+1},b_{l+1}}$ is the prediction.\footnote{By definition, the prediction from the input layer is $a_{t}^{ \theta_{2}, b_{2}} = g_{1} \Big( \sum_{j=1}^{N_{1}} \theta_{j}^{(1)}Z_{t} + b^{(1)} \Big)$.}

The output layer constructs a forecast:

As seen, a neural network applies a series of functional transformations to construct the forecast. In principle, a single hidden layer with a large enough number of neurons (and appropriate activation function) is sufficient to approximate any continuous function. This follows from the Universal Approximation Theorem of \citet*{cybenko:89a}. In practice, however, it is often convenient and computationally more efficient to add extra hidden layers than to arbitrarily increase the number of neurons in a layer. Hence, the optimal architecture of a neural network depends on the problem and is determined through tuning.

To allow inspection of the inner workings of a neural network, we construct four models with an architecture inspired by the geometric pyramid \citep*[see, e.g.,][]{masters:93a, gu-kelly-xiu:20a}. The first, denoted NN$_{1}$, has a shallow structure with a single hidden layer and 2 neurons. NN$_{2}$ is two-layered with 4 and 2 neurons, NN$_{3}$ has three layers with 8, 4, and 2 neurons, and NN$_{4}$ has four layers with 16, 8, 4, and 2 neurons.

We apply the Leaky Rectified Linear Unit (L-ReLU) by \citet*{maas-hannun-ng:13a} as activation function for all layers, as it is able to infer when the activation is zero via gradient-based methods, in contrast to the standard ReLU.\footnote{Other common activation functions include Sigmoid, TanH, or ReLU. Switching to one of these has no material impact on our empirical results.} The L-ReLU is defined as:

where $c \geq 0$.\footnote{In this paper, $x = \sum_{j=1}^{N_{l}} \theta_{j}^{(l)}a_{t}^{ \theta_{l},b_{l}} + b^{(l)}$. We follow the majority of the literature by setting $c = 0.01$.}

The weight parameters are chosen to minimize the squared error loss function.\footnote{Alternatively, a penalized $L^{2}$ objective function can be applied. We found no additional improvement when combining $L^{2}$ and drop-out.} In lack of a closed-form solution, we adopt the Adaptive Moment Estimation (Adam\glsadd{Adam}) by \citet*{kingma-ba:14a}, where the optimizer updates the weights $( \hat{ \theta}, \hat{b})$ as a combination of Momentum\glsadd{Momentum} and RMSprop\glsadd{RMSprop} (Root Mean Square propagation).\footnote{To limit the number of hyperparameters and reduce the complexity of the optimization, we select exponential decay rates in line with \citet*{kingma-ba:14a}.} We further predetermine the learning rate\glsadd{Learning Rate} at 0.001, see \citet*{kingma-ba:14a}.

Regularization of a neural network is crucial to avoid overfitting. There is a large literature proposing different regularization techniques. In this study, we employ four of the most commonly used, namely learning rate shrinkage, drop-out, early stopping, and ensembles. Ensembling is a model averaging technique, in which numerous neural networks are fitted, and the prediction is constructed from a subset of the best performing models from the validation set. The idea is to reduce the variation in the forecast caused by the random initial weights assigned to a neural network by repeating the experiment several times, much as an ensemble average. We investigate ensembles of 1---corresponding to no averaging---and 10 (out of 100), which we denote NN$_{1}^{1}$ and NN$_{1}^{10}$ for the single hidden layer neural network, and so forth. In Appendix (ref) -- (ref), we explain the regularization techniques and choice of hyperparameters in more detail.

\tikzset{ every neuron/.style={ circle, draw, minimum size=0.5cm }, neuron missing/.style={ draw=none, scale=2, text height=0.2cm, execute at begin node=\color{black}$\vdots$ }, }

Forecast comparison

As the mean squared error (MSE) laid the foundation for the optimal structure of the ML algorithms, a logical continuation is to employ it also as an out-of-sample evaluation measure:

where $\hat{f}(Z_{t})$ is the forecast.\footnote{The LogHAR produces a forecast $\hat{f}(Z_{t})$ of log-realized variance, but we require a forecast of realized variance $\exp \big( \hat{f}(Z_{t}) \big)$. This is a nonlinear transformation, so the forecast of realized variance is biased by Jensen's inequality, even if the forecast of log-realized variance is not. We bias-correct as follows: $E\left[ \exp \big( \hat{f}(Z_{t}) \big) \right] = \exp \left(E \big[ \hat{f}(Z_{t}) \big] + 0.5 \text{var} \big[ \hat{f}(Z_{t}) \big] \right)$, where $\text{var} \big( \hat{f}(Z_{t}) \big)$ is the variance of the residuals in the training and validation set. This is suitable when the distribution of log-realized variance is Gaussian, which is approximately true in practice, see, e.g., \citet*{andersen-bollerslev-diebold-ebens:01a, christensen-thyrsgaard-veliyev:19a}.}

If a model predicts volatility to be negative, we replace the forecast with the minimum in-sample realized variance $\min_{t \in \text{in-sample}} RV_{t}$. We also adopt the insanity filter from \citet*{bollerslev-patton-quaedvlieg:16a} for the HARQ.

To measure statistical significance, we compute a pairwise \citet*{diebold-mariano:95a} test for equal predictive accuracy with a one-sided alternative that the comparison model beats the benchmark model, as explained later. We also construct a Model Confidence Set (MCS) of \citet*{hansen-lunde-nason:11a}. It defines a collection of models containing the “best” one with a given level of confidence. Inferior forecasting models are removed via an elimination rule.

Data description

The empirical investigation is based on high-frequency data from 29 of the 30 Dow Jones Industrial Average (DJIA) constituents prior to the recomposition on August 31, 2020. The complete list of ticker symbols is AAPL, AXP, BA, CAT, CSCO, CVX, DIS, DOW, GE, GS, HD, IBM, INTC, JNJ, JPM, KO, MCD, MMM, MRK, MSFT, NKE, PFE, PG, RTX, TRV, UNH, VZ, WMT, XOM.\footnote{We also ran the forecast comparison for the composite DJIA index and the SPY. The latter is a liquid ETF tracking the S&P 500. The outcome of this analysis confirms the findings reported for individual equities, suggesting our conclusions hold more broadly on a market-wide level. Moreover, we attempted to forecast volatility of the DJIA index based on cross-sectional information of the volatilities from the individual member firms. However, we did not find any improvement in forecast accuracy, once we controlled for past DJIA index volatilities.} Visa (V) is excluded due to limited data availability, as its shares were only listed for public trading on NYSE on March 19, 2008. We construct a continuous stock price record for each company on a sample from 29 January 2001 to 31 December 2017, or $T = 4{,}257$ trading days. The sample includes several noteworthy events, such as the financial crisis, the European sovereign debt crisis, several rounds of negotiations about the limit of the U.S. debt ceiling, and the flash crashes of 6 May 2010 and 24 August 2015.

There has been several corporate actions changing the ticker symbol history for some of the current members of DJIA. We include high-frequency information from the main predecessor prior to such events. In particular, we employ 1) Chevron (CHV) prior to its merger with Texaco on October 9, 2001 to forge the company ChevronTexaco (CVX), 2) St. Paul Travelers Companies (STA) prior to its name change to Travelers Inc. (TRV) on February 27, 2007, 3) United Technologies (UTX) prior to its replacement by Raytheon Technologies (RTX) on April 6, 2020 (Raytheon is the fusion of United Technologies and Raytheon Company, which merged as of April 3, 2020), 4) Walgreens (WAG) prior to forming Walgreens Boots Alliance (WBA) by purchasing a majority stake in Alliance Boots on December 31, 2014, and 5) DowDuPont (DWDP) before the spin-off of Dow Chemical Company (DOW) on April 2, 2019. DowDupont itself was the result of a previous merger between Dow Chemical Company and DuPont (DD) on September 1, 2017, and here we take DuPont as predecessor.

The data is extracted from the NYSE Trade and Quote (TAQ) database. Prior to our analysis, the data were pre-processed and filtered for outliers using standard algorithms, see \citet*{barndorff-nielsen-hansen-lunde-shephard:09a, christensen-oomen-podolskij:14a}. Based on the cleaned transaction price record, we construct a five-minute log-return series and compute for each asset a time series of daily realized variance $RV_{t}$ with $n = 78$, for $t = 1, \dots, 4{,}257$. The goal is to forecast realized variance.

We split the data into a training\glsadd{Training Set}, validation\glsadd{Validation Set}, and test set\glsadd{Test Set}. The primary analysis is restricted to a training set of 70% of the dataset (or 2,964 days), a validation set of 10% (or 424 days), and a test set of 20% (or 847 days).\footnote{A month worth of high-frequency data (or 22 days) are “lost” due to the lagging.} The latter are reserved for out-of-sample evaluation. We conduct a robustness check, where the training set is modified to 1,000 and 2,000 days (with a validation set fixed at 200 days). The results, which are broadly in line with our reported findings, are available in Appendix (ref).

In the non-regularized HAR models, we merge the training and validation set and employ time-varying parameters by constructing a rolling window forecast. This approach is also adopted for bagging and random forest to avoid any hyperparameter tuning\glsadd{Hyperparameters}, which results in off-the-shelf implementation of these methods. However, tuning is almost unavoidable for ridge, lasso, elastic net, and gradient boosting. In keeping with the spirit of our paper, a rolling scheme without concatenation of the training and validation set is conducted. An in-depth explanation of our tuning of hyperparameters is available in Appendix (ref). The vast majority of the hyperparameters are set equal to their default values proposed in the original research papers, often based on very different contexts. At last, we employ a fixed window estimation for the neural networks, where the weights are only found once in the initial validation sample and not rolled forward in the out-of-sample window. This implementation is strictly advantageous to the other competing models, which are allowed to adapt to new information. It may be possible to further improve the performance of the neural networks with a rolling window approach, however this requires a very large computer capacity, which is outside our budget.

To study the effect of tagging on additional explanatory variables, we construct two datasets $\mathcal{M}_{ \text{HAR}}$ and $\mathcal{M}_{ \text{ALL}}$. $\mathcal{M}_{ \text{HAR}}$ contains the daily, weekly, and monthly lag of realized variance (RVD, RVW, and RVM).\footnote{$\mathcal{M}_{ \text{HAR}}$ also includes variables that are employed in the various extensions of the HAR model. However, these variables are only passed to the concrete model (i.e., the square-root of realized quarticity appears as an interaction term with realized variance in HARQ, but not in the other models).} This ensures one-to-one comparison between the HAR models and ML algorithms. $\mathcal{M}_{ \text{ALL}}$ extends $\mathcal{M}_{ \text{HAR}}$ by including additional variables examined in previous research. We select a total of nine other key variables that are arguably important predictors of future volatility. In particular, we extract four firm characteristics: Model-free implied volatility (IV),\footnote{The data was downloaded from the OptionMetrics' IvyDB database, which estimates the implied volatility for each traded option contract. We followed the data cleaning suggested by \citet*{lei-wang-yan:20a} to deal with erroneous data points.}\footnote{To account for the correlation between implied volatility and realized variance, we also estimated an auxiliary model: $IV_{t} = \beta_{0} + \beta_{1}RV_{t} + u_{t}$ and employed the residuals from the fitted regression in place of implied volatility. The results were roughly unchanged.} an indicator for earnings announcements (EA), 1-week momentum (M1W), and dollar trading volume (\$VOL). We include five macroeconomic indicators: The CBOE volatility (VIX) index, the Hang Seng stock index daily squared log-return (HSI), the \citet*{aruoba-diebold-scotti:09a} business conditions (ADS) index, the US 3-month T-bill rate (US3M) and the economic policy uncertainty (EPU) index from \citet*{baker-bloom-nicholas:16a}.\footnote{We also constructed an intermediate dataset via subset selection\glsadd{Subset selection}. In particular, we fitted a separate least squares regression for each combination of the feature\glsadd{Feature} space of $\mathcal{M}_{ \text{ALL}}$. The best model was then selected with a BIC criterion. It contains the daily and weekly lag of realized variance, implied volatility, 1-week momentum, and the earnings announcement indicator, i.e. a total of five variables. The conclusions from the forecast comparison on this dataset is consistent with what we report below. Interestingly, if the evaluation criterion is changed from BIC to cross-validation error, the subset selection includes only the standard three HAR variables.} We first-difference US3M to account for possible nonstationarity, and we replace \$VOL by first-differencing its log-transform. Moreover, in the LogHAR we also log-transform VIX and IV.

A full list of explanatory variables and their acronyms is presented in Table (ref), along with some standard descriptive statistics. Prior to estimation, we standardize the input data with the sample mean and sample variance from the training set\glsadd{Training Set} to render the estimation scale invariant and improve the numerical optimization procedure for ridge, lasso, elastic net, and the neural network.

Empirical results

One-day-ahead forecasting

In Table (ref), reading by column we report the out-of-sample MSE of each forecasting model relative to the benchmark in the selected row for the dataset $\mathcal{M}_{ \text{HAR}}$. The ratio is a cross-sectional average of the relative one-day-ahead realized variance forecast MSEs for each stock. Hence, the first row contains the out-of-sample MSE of every model relative to the basic HAR, averaged over stocks. In addition, for each stock and pairwise comparison we compute the Diebold-Mariano test. The number formatting shows whether the null hypothesis of equal predictive accuracy for a fixed pairwise comparison was rejected more than half of the times across stocks at various levels of significance. Note that in $\mathcal{M}_{ \text{HAR}}$ HAR-X is restricted access to lagged realized variance and is therefore identical to HAR.

Several interesting findings emerge. Firstly, the basic HAR is superior to its offspring. Compared to the extant literature, this is surprising. However, the results naturally depend on the split between training, validation, and test set. With a conventional shorter training set of 1,000 days, the extended HAR models outperform HAR by some distance, as reported in Appendix (ref). A possible reconciliation of this is that with a shorter in-sample size, the financial crisis appears in the out-of-sample set, which has a huge impact on the relative performance of the models. For example, the quarticity correction in the HARQ is arguably more crucial in high volatility environments. Secondly, when the covariate list is restricted to $\mathcal{M}_{ \text{HAR}}$ the ML results are also mixed. This is not surprising. The lagged values of realized variance are strong determinants of future volatility, so there is not much room for improvement. In alignment with this interpretation, regularization is at par with HAR sporting an aggregated MSE ratio at or above one. So, as expected, there is no impending need to penalize the HAR model in $\mathcal{M}_{ \text{HAR}}$. The regression trees are also inferior to HAR, especially bagging, whereas the random forest is about equal in performance. In contrast, the cross-sectional relative MSEs for the neural networks are consistently below one by roughly five percent on average, and the Diebold-Mariano test is frequently rejected at the 10% or 5% significance level, when a neural network is evaluated against the HAR as a benchmark. This shows that for short-horizon forecasts of next-day volatility there are minor refinements associated with allowing for nonlinearities in the dynamic. What is perhaps the most striking outcome of Table (ref) is that even a shallow neural network beats the HAR. This can possibly be explained by the fact that a limited complexity of a neural network allows the benefits of early stopping and dropout to be emphasized.

Turning next to Table (ref), we show the corresponding results for $\mathcal{M}_{ \text{ALL}}$. Here, the supplemental exogenous predictors from Table (ref) appear in the information set. The addition of these variables causes a huge impact on the pairwise comparisons. This forcefully demonstrates the complexity and low signal-to-noise ratio in the underlying financial data, as inferred from the MSE ratio of the HAR-X and HAR. As gauged by comparing Table (ref) and Table (ref), there is but a minuscule decrease in MSE of $3.4\%$ moving from $\mathcal{M}_{ \text{HAR}}$ to $\mathcal{M}_{ \text{ALL}}$. Furthermore, the LevHAR, SHAR, and HARQ fare worse than HAR, even though their information set is larger. This supports our argument in Section (ref) that lack of regularization leave the linear HAR models exposed to in-sample overfit. As such, ridge regression is expected to alleviate this problem. This is also confirmed in Table (ref). The $\ell_{2}$ regularization enforced by ridge reduces the average MSE of $8.1\%$ overall. The effect of allowing for subset selection\glsadd{Subset selection} can be uncovered by contrasting lasso with ridge. The latter is marginally better than the former, but combining ridge and lasso via the elastic net is best overall.\footnote{The average value of the weight parameter $\alpha$ is about 0.165 in our sample.} The elastic net ranks on top of the linear predictions with a pronounced $8.4\%$ reduction in MSE relative to the HAR.

An appealing observation that can be made from Table (ref) is that the dimension of the feature space is rather critical to achieve good performance with regression trees. In $\mathcal{M}_{ \text{HAR}}$, their performance were inferior to HAR, but with the expanded database they work better, although bagging and gradient boosting are still worse than a regularized HAR. The best-in-class tree-based approach is the random forest. It delivers a $9.9\%$ reduction in MSE against HAR. This is arguably caused by its ability to decorrelate trees via random selection of the number of explanatory variables, as opposed to bagging and gradient boosting. Comparing random forest to gradient boosting, the latter shows inferior forecast accuracy. Since gradient boosting produces a sequential model based on an ensemble of weak learners\glsadd{Weak learner}, it places more weight on outlying observations in contrast to classical tree-based approaches. Hence, it is unsurprising that random forest does a better job than gradient boosting, given the low signal-noise ratio in the data.\footnote{We also estimated an XGBoost model of \citet*{chen-guestrin:16a} with a refined set of hyperparameters, early stopping and $\ell_{1}$ penalization. However, the results were not much different from gradient boosting.} It should be stressed that even though random forest has the lowest MSE among linear regression and tree-based models, there is only a marginal drop between it and the elastic net, which is moreover typically not statistically significant. This highlights that automatic interactions between variables can be informative about future volatility. Meanwhile, tackling overfitting is a more potent source of forecast improvement in the “high-dimensional” $\mathcal{M}_{ \text{ALL}}$ framework compared to $\mathcal{M}_{ \text{HAR}}$.

At last, we compare the neural networks to HAR. It is evident that the average MSE decreases further, such that the neural networks remain in front. However, the enlargement of the feature space benefits several models, and the edge of the neural networks is roughly as before. What is more, when additional predictors are included, the improvement is statistically significant at the 1% level across more than half of the stocks.

To shed further light on the causes of the decline in MSE for the neural networks, a logical next step is to inspect their complexity. In this aspect the geometric pyramid structure is helpful, as we can readily investigate the effects of changing network architecture. Indeed, our results show that there is no big demand for an overly complex model, because the results are rather stable across network configurations. Initially, adding hidden layers to the network is helpful. Peak performance is attained with a NN$_{2}$ -- NN$_{3}$ architecture, but it starts to deteriorate after that. Moreover, ensemble averaging the forecast over several trained neural networks leads to further declines in MSE, except again for NN$_{4}$.

Finally, we highlight that the best HAR model is the LogHAR. At first glance it remains competitive against the ML techniques even in the $\mathcal{M}_{ \text{ALL}}$ setting. Thus, while a regularized HAR outperforms a non-regularized HAR, it is in turn overtaken by LogHAR. As the latter is nonlinear, this implies that handling nonlinearities is a key source to improvement, also evidenced by the ML techniques. This confirms that regularization is less vital than capturing nonlinearities, which may partially be explained by the fact that $\mathcal{M}_{ \text{ALL}}$ is not too vast in our setting.\footnote{We only detect minor improvements in forecast accuracy when regularizing the LogHAR. For example, in $\mathcal{M}_{ \text{ALL}}$ the relative MSE of the LogHAR regularized with elastic net to the unregularized version is 0.995. We believe this is because the log-transformation reduces the influence of outliers and attenuates the need for penalization.} However, in contrast to ML the true story behind the success of the LogHAR is not to be found in a general leap in forecast accuracy, but rather in a more narrow improvement in high states of realized variance, which is a reflection of the log-transformations capability to bring down outliers. This becomes apparent when we take an in-depth look at the forecast accuracy over the entire realized variance distribution in Figure (ref).

To sum up, the main improvement in forecasting accuracy in $\mathcal{M}_{ \text{ALL}}$ is due to the inclusion of additional explanatory variables, which generally aids the ML techniques and the LogHAR, whereas the unregularized linear HAR models suffer from overfitting. Within the ML models, the ones that capture non-linearities are superior, and regularization turns out to be more potent when growing the feature space.

To examine whether the superior forecast accuracy of the ML methods is driven by a small subset of stocks or it extends to the entire cross-section, we next dissect the aggregate numbers reported in Table (ref) -- (ref). In Figure (ref), we construct a boxplot of the relative MSE for each model versus the HAR across the 29 stocks in our sample. Panel A is for $\mathcal{M}_{ \text{HAR}}$, whereas Panel B is for $\mathcal{M}_{ \text{ALL}}$. The tree-based approaches are slightly inferior to the HAR for $\mathcal{M}_{ \text{HAR}}$, as also noticed above. However, there is a steady and sustained improvement in forecast accuracy of ML as we grow the information set.

In Figure (ref), we compute the MCS at a 90% confidence level. The figure shows the percentage of times a model was retained in the final set across stocks. In line with our earlier results, the neural networks exhibit an inclusion rate of around $80\%$ for the $\mathcal{M}_{ \text{ALL}}$ dataset, compared to about 20% -- 40% for the HAR. The LogHAR is the sole HAR model that is retained at a high rate in the $\mathcal{M}_{ \text{ALL}}$ setting.

To get insights about how good the models forecast across states of the volatility distribution, we split the test set into deciles of the observed values of daily realized variance (i.e., the first subsample holds the $10\%$ smallest values of realized variance in the out-of-sample window and is denoted (0.0,0.1)). For brevity we select a “preferred” model from each main category, i.e. we look at HAR-X, LogHAR, elastic net, random forest, and NN$_{2}^{10}$.

In Figure (ref), we present the outcome of these calculations for selected low, medium, and high deciles. As indicated, the random forest and neural network yield lower MSE compared to the HAR across the entire support of the volatility distribution. Hence, the key driver for the ML results can be found in the majority of the sample. The MSE is only reduced by $7.38\%$ for the highest 10% of the out-of-sample realized variancee, suggesting a closer alignment in forecast accuracy during market turmoil, which is of course crucial in practice. The sustained results are not shared with the HAR-X, LogHAR and elastic net. On the one hand, the HAR-X struggles when volatility is high, but has a much better MSE in tranquil markets. On the other hand, LogHAR yields only a minor improvement compared to the HAR, which appears to contradict the conclusion of Table (ref) and (ref). However, a closer look reveals only a sizable improvement of the LogHAR in the highest decile. As realized variance is highly right-skewed with extreme outliers, it is therefore the ability of the LogHAR to downweigh outliers that brings its average squared error markedly down.

Marginal effects and variable importance

The autonomy of a neural network or tree-based approach severely limits our ability to explain what is going on “under the hood”, which often leads ML to be accused of being a black box approach. With their complex and highly nonlinear structure, a complete revelation of the underlying machinery is of course also rather difficult.

Nonetheless, we can illustrate the inner workings of the ML algorithms and shed at least some initial light on the black box by exploiting a recent contribution of \citet*{apley-zhu:20a} that allows to compute an intuitive measure of marginal effect and variable importance in the ML framework. In particular, we calculate an Accumulated Local Effect (ALE) plot, which is an in-sample statistic based on the fitted model. To explain the concept, suppose for ease of exposition that $f(Z_{t})$ is differentiable and that we are interested in analyzing the marginal effect of the $j$th predictor $Z_{jt}$, $j = 1, \dots, J$, from the fitted model $\hat{f}(Z_{t})$. We divide $Z_{t}$ into $Z_{jt}$ and the complement $Z_{ \complement t} = Z_{t} \setminus Z_{jt}$, i.e. the remaining predictors in $Z_{t}$.

The ALE on $\hat{f}(Z_{t})$ of $Z_{jt}$ is the function:

where $p_{Z_{ \complement t} \mid Z_{jt}}(z_{ \complement t} \mid x)$ is the conditional joint density of $Z_{ \complement t}$ given $Z_{jt} = x$ and the function is defined over the range of observed values of $Z_{jt}$ in the training set, as explained below. The constant is such that $E[f^{ \text{ALE}} (Z_{jt})] = 0$.

The term $\displaystyle \frac{ \partial \hat{f}(x, z_{ \complement t})}{ \partial x}$ is called the local effect of $Z_{jt}$ on $\hat{f}(Z_{t})$. By accumulating the local effect, the ALE can be interpreted as the combined impact of $Z_{jt}$ at the value $z$ compared to the average prediction of the data. The derivative isolates the effect of the target predictor and controls for correlated features. Hence, ALE calculates changes in the prediction and adds them up, rather than computing the prediction directly. As explained in \citet*{apley-zhu:20a}, this circumvents problems with omitted variable biases. Moreover, by exploiting the conditional distribution the ALE is not impaired by unlikely data points, a key weakness of the more common Partial Dependence (PD) function of \citet*{friedman:01a}.

$f^{ \text{ALE}}$ can be estimated via a finite difference approach. The information set used in the estimation is the matrix $Z = (Z_{1}, \ldots, Z_{T_{0}})$, where $T_{0}$ is the number of days in the training sample. The observations of the $j$th predictor correspond to the entries in the $j$th row of $Z$, which we denote $Z_{j} = (z_{j1}, \ldots, z_{jT_{0}})$. $Z_{ \complement}$ is the complement of $Z_{j}$, i.e. the remaining rows in $Z$. We partition the range of $Z_{j}$ such that $\min(Z_{j}) - \epsilon = z_{0} < z_{1} < \dots < z_{K} = \max(Z_{j})$, where $\epsilon$ is a small positive number ensuring that $\min(Z_{j})$ is included in the subsequent calculations.\footnote{In our implementation, we partition $Z_{j}$ into 100 subintervals containing equally many observations.} Next, as a function of $z \in (z_{0}, z_{K}]$ denote by $\text{idx}(z)$ the index of the interval into which $z$ falls, i.e. $z \in ( z_{\text{idx}(z)-1}, z_{\text{idx}(z)}]$. The uncentered ALE is estimated as follows:

where $T_{k}$ is the number of data points of $Z_{j}$ that fall in $(z_{k-1}, z_{k}]$ with $T_{0} = \sum_{k=1}^{K} T_{k}$. As readily seen, the second summation approximates the local effect across small intervals by averaging over complement predictors $Z_{ \complement}$ that are associated with the observations of $z_{jt} \in (z_{k-1}, z_{k}]$. The first sum accumulates these effects.

The centered ALE is then:

An ALE plot is based on the graph of $(z, \hat{f}^{ \text{ALE}} (z))$.

We convert ALE to a variable importance measure following \citet*{greenwell-boehmke-mccarthy:18a} in the PD framework by computing the sample standard deviation of $\hat{f}^{ \text{ALE}}$, i.e.

Intuitively, if our prediction is only weakly related to $Z_{j}$, once we control for the effects of other predictors, $Z_{j}$ has a low importance. Consistent with this, the importance score is then $I(Z_{j}) \simeq 0$. On the other hand, if $\hat{f}^{ \text{ALE}}$ fluctuates a lot over its domain, it means that $Z_{j}$ has a large influence on the response variable, leading to larger values of $I(Z_{j})$.

Proceeding in this way, we construct a sequence $I(Z_{j})$, for $j = 1, \dots, J$. Our measure of variable importance is based on expressing each $I(Z_{j})$ relative to their sum:

Note that $VI(Z_{j}) \geq 0$ and $\sum_{j=1}^{J} VI(Z_{j}) = 1$.

In Figure (ref), we plot the ALE on expected realized variance over the interval [-1,1] standard deviation. We restrict attention to a set of selected characteristics that are the most influential drivers of volatility, namely RVD, RVW, IV and M1W. Moreover, we do not explore the marginal effect averaged over stocks, as it smoothes out the unique relationship at the individual asset level. Instead, we arbitrarily select Apple for the illustration. Hence, this part serves as a showcase and may not be representative for the whole sample, although what we observe here is consistent with the broad pattern across stocks. In the construction of Figure (ref), we again select a single model from each main category.

Firstly, the figure reveals a positive relationship between $RV_{t}$ and RVD, RVW, and IV. In contrast, the marginal association between $RV_{t}$ and M1W is negative. This is corroborated by Table (ref) in Appendix (ref), where parameter estimates for the HAR-X model of Apple are reported. Consistent with Figure (ref), the slope coefficient for M1W is negative and significant. Secondly, the figures illustrate a purely data-driven approach leads to a slightly nonlinear model. Interestingly, the marginal effect for RVD detected by the random forest and NN$_{2}^{10}$, especially the former, are closer to the LogHAR than the HAR. The ML algorithms therefore detect a stronger relationship between $RV_{t}$ and RVD compared to the linear model, whereas the impact of RVW is attenuated. The discrepancy between elastic net and HAR-X is generally caused by the shrinkage of the parameter estimates, such that elastic net ends up suggesting a marginal effect that is more consistent with the random forest.

The cross-sectional average VI measure is reported in Figure (ref). As illustrated, there is agreement about the dominant set of predictors: RVD, RVW, and IV. Intuitively, these covariates relate to the recent history or future expectation of the stock's volatility. In contrast, less weight is assigned to variables that are not directly related to volatility. However, even if there is consensus about the best set of predictors, there are individual differences in their relative importance. It is interesting to note that IV ranks higher in the random forest and neural network compared to the HAR-X and elastic net. Hence, the ML models weigh market expectations of future volatility higher than past realizations of realized variance in forming their forecast, adding an ML perspective to the literature on the informational content of realized variance versus implied volatility \citep*[e.g.,][]{christensen-prabhala:98a,busch-christensen-nielsen:11a}.

With a nonlinear market structure, valuable information can potentially be extracted from a broader set of covariates. In line with this conjecture, the macroeconomic variable ADS contributes more to the prediction of volatility in the random forest relative to RVM, whereas the ADS variable does not weigh much in the HAR-X. This is again consistent with the parameter estimates of the HAR-X model for Apple in Table (ref), where the coefficient in front of ADS is found to be insignificant in the $\mathcal{M}_{ \text{ALL}}$ dataset.

Lasso applies subset selection. Accordingly, it recovers a sparse solution with a median of four beta coefficients set to zero across stocks. Thus, lasso is capable of taming overfitting, but this prevents it from extracting incremental information from several of the explanatory variables, as opposed to elastic net, random forest or neural network.

The EA variable has a low VI score. This highlights a weakness with ALE, since it is mainly a tool for finding the drivers behind the average prediction. Earnings information is released infrequently, but it is arguably very important for making short-run predictions of realized variance in the days surrounding the announcement, see \citet*{siggaard:22a}.

Longer-run forecasting

In this section, we modify the forecast horizon from one-day-ahead to one-week-ahead and one-month-ahead. We replace the dependent variable in all models from next-day realized variance ($RV_{t+1}$) to next-week average realized variance ($RV_{t+5 \mid t+1}$) and next-month average realized variance ($RV_{t+22 \mid t+1}$), as a function of the information set that is available at time $t$, i.e. $Z_{t}$. Table (ref) -- (ref) report the results for the weekly forecasts, whereas the monthly forecast results are available in Table (ref) -- (ref).

The results are striking. On the one hand, changing the forecast horizon does not alter the conclusions for $\mathcal{M}_{ \text{HAR}}$ in any discernible way. The neural networks are still dominating and perhaps marginally better than before. On the other hand, the picture changes for $\mathcal{M}_{ \text{ALL}}$, where regression trees and neural networks now provide outstanding improvements in forecast accuracy. Moreover, they appear to get better as we proceed from the weekly to monthly horizon. Whereas neural networks were ahead at the one-day-ahead forecast horizon, we observe that regression trees start to be competitive at the one-week-ahead, and even surpass neural networks at the one-month-ahead, horizon. This shows that while nonlinearities remain important at longer horizons, functional form and interaction effects inherent in the regression trees are more crucial when expanding the predictor base and forecast horizon. The size of the reduction in relative MSE is large with a 40% or more improvement in forecast accuracy on several occasions, although one should observe that this is partly driven by a deterioration in performance for some HAR models. Even gradient boosting is outperforming most of the HAR models. The differences are often significant at the 1% level for more than 50% of the included stocks.

The LogHAR remains ahead of the HAR pack, but the Diebold-Mariano test of equal forecast accuracy with this model as the benchmark is now routinely rejected in favor of ML models. As indicated in Figure (ref), the main driver behind the LogHAR model is it ability to handle outliers. As longer-run average realized variance is smoother, it is not surprising that LogHAR cannot compete with more sophisticated models, such as random forest and neural networks. This also explains the mediocre results for regularization, since this is less critical here.

Why does random forest suddenly emerge as the lead forecaster of future volatility at the one-month horizon? An understanding of this can help to shed light on where the improved forecast gains are extracted from. In Figure (ref), we plot the in-sample acf of the fitted realized variance series extracted from the initial training set for the HAR, random forecast, and NN$_{2}^{10}$. The acf embodies how each model perceives the persistence of the realized variance series. In Panel A, we plot the acf corresponding to one-day-ahead forecasts, whereas Panel B is for the one-month ahead forecasts. We observe that at the short horizon, the models are roughly equally persistent, albeit a bit lower for random forest. However, at the monthly horizon random forecast and neural networks display a visibly higher acf, suggesting they are better at approximating an underlying long-memory structure in realized variance, which is harder to detect at the one-day-ahead window, but becomes more pronounced at longer horizons.

A Value-at-Risk analysis

We now turn our attention toward a more economically motivated application by looking at volatility forecasting through the lens of a Value-at-Risk (VaR) framework. We follow \citet*{audrino-huang-okhrin:19a} and construct one-day-ahead forecasts $\widehat{ \text{VaR}}_{t+1}^{ \alpha}$ of $\text{VaR}_{t+1}^{ \alpha}$, where $\mathbb{P}(r_{t+1} \leq \text{VaR}_{t+1}^{ \alpha}) = \alpha$, which are obtained by means of filtered historical simulation \citep*{barone-adesi-bourgoin-giannopoulos:98a, barone-adesi-giannopoulos-vosper:99a}. The VaR forecasts are evaluated using a loss function that is commonly applied in quantile regression \citep*[e.g.,][]{koenker-bassett:78a}, namely

where $d_{t+1} = 1_{ \{ r_{t+1} < \widehat{ \text{VaR}}_{t+1}^{ \alpha} \}}$ is the “hit” function. This loss function is asymmetric and associates a weight $1- \alpha$ to observations for which the daily log-return is beneath the forecasted quantile, which is arguably more critical in practice, especially for negative returns. Otherwise, a weight $\alpha$ is assigned. Throughout, we describe the results for $\alpha = 0.05$. We also investigated $\alpha = 0.01$ without any major shift in the outcome.

We also examine the exceedance probabilities of the hit sequence by conducting the likelihood ratio tests of \citet*{kupiec:95a} and \citet*{christoffersen:98a} to inspect the unconditional coverage and independence property of the VaR forecasts. The first inspects whether the empirical rejection rates are statistically speaking equal to $\alpha$, on average, whereas the second test checks this assumption conditionally on the previous outcome.

A graphical representation of the VaR forecasts for Apple are presented in Figure (ref) for a few models. As evident, they capture the temporal clustering of volatility. The forecast comparison and exceedance probability tests are reported in Table (ref) -- (ref). Overall, the main conclusions from the previous sections are still valid, although the VaR loss differentials are smaller and less statistically significant than before. Interestingly, the basic HAR model is remarkably solid and tough to beat in the VaR framework. The ML techiques are at best on par with the HAR in the $\mathcal{M}_{ \text{HAR}}$ setting, but we do again observe some more substantial improvements in the $\mathcal{M}_{ \text{ALL}}$ dataset, where many of the unregularized HAR models continue to struggle. Moreover, from the unconditional coverage test we see that ML is slightly worse than the HAR and LogHAR, whereas conditional on the previous outcome they are equally good. The overall impression of the results is that ML delivers precise measures of the Value-at-Risk.

Conclusion

This paper synthesizes the field of machine learning (ML) with volatility forecasting. We conduct an extensive out-of-sample forecast comparison of current ML algoritms with a broad suite of HAR models for forecasting realized variance. We show that ML is better at handling the highly nonlinear structure governing financial markets, and that they are able to extract valuable incremental information about future volatility in a data-rich environment with a large feature space; a setting in which standard HAR models struggle. The random forest and neural networks are preferred among ML algorithms and significantly improve prediction relative to HAR. The improvements increase at longer horizons and are often statistically significant. We dissect the underlying structure of the data via ALE plots in order to pin down the main covariates leading to the improved prediction. In an economic Value-at-Risk application, we continue to observe reduced loss with ML models, although it is less statistically significant.

To the best of our knowledge, this study is the first to extensively implement a broad selection of the most widely applied ML methods for the problem of volatility prediction and compare them to a broad suite of HAR models. We envision the paper can serve as a starting point for applying ML to forecast stochastic volatility and perhaps even as a benchmark when building more complicated models. Moreover, we look forward to behold how much additional improvement that can be had by doing more extensive hyperparameter tuning than what we do in this paper.