EconBase
← Back to paper

GARCHX-NoVaS: A Model-free Approach to Incorporate Exogenous Variables

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

59,939 characters · 16 sections · 42 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

GARCHX-NoVaS: A Model-free Approach to Incorporate Exogenous Variables

abstractIn this work, we explore the forecasting ability of a recently proposed normalizing and variance-stabilizing (NoVaS) transformation with the possible inclusion of exogenous variables. From an applied point-of-view, extra knowledge such as fundamentals- and sentiments-based information could be beneficial to improve the prediction accuracy of market volatility if they are incorporated into the forecasting process. In the classical approach, these models including exogenous variables are typically termed GARCHX-type models. Being a Model-free prediction method, NoVaS has generally shown more accurate, stable and robust (to misspecifications) performance than that compared to classical GARCH-type methods. This motivates us to extend this framework to the GARCHX forecasting as well. We derive the NoVaS transformation needed to include exogenous covariates and then construct the corresponding prediction procedure. We show through extensive simulation studies that bolster our claim that the NoVaS method outperforms traditional ones, especially for long-term time aggregated predictions. We also provide an interesting data analysis to exhibit how our method could possibly shed light on the role of geopolitical risks in forecasting volatility in national stock market indices for three different countries in Europe.

\JELcodes{C32; C53; C63; Q54}

Introduction

In the long history of time series econometrics literature, accurate forecasting has always stood out as a fundamental and important problem. It has a range of applications in various industries, e.g., weather forecasting, climate forecasting, and economic forecasting. Discrete-time series data, e.g., heights of ocean tides, and temperature of a city, is the realization of a stochastic process $\{X_t,t\in \mathbb{Z}\}$. The earliest modern time series analysis could be traced back to the work of yule1927vii where the pattern of the sunspots number was studied. Unlike the prediction of independent data, the prediction of time series gets more complicated due to the inherent data dependence. To get accurate predictions and inferences, it is crucial to model the dependent relationship within the data. Usually, very generally speaking, the time series data is assumed to be generated by some underlying mechanism as follows:

equation[equation omitted — 85 chars of source]

$G(\cdot,\cdot)$ could be any suitable function; $\epsilon_t$ is called innovation and assumed to be $i.i.d.$ with appropriate moments and independent with $X_{t-i}$, $i\geq 1$; $\bm{X}_{t-p}$ represents $\{X_{t-1},\ldots,X_{t-p}\}$ and stands for the historical information. To further simplify the forecasting problem, participators focus on some standard formats of $G(\cdot,\cdot)$, e.g., linear or non-linear. For linear models, such as linear AR, MA and ARMA models, we can apply the Box-Jenkins method of identifying, fitting, checking and predicting models systematically box1976time. However, the prediction of non-linear models is not as trivial as the case of linear models since the innovation must be appropriately included in the prediction process, especially for the multi-step ahead predictions; see wu2023bootstrap for more related discussions.

In this paper, we are exclusively interested in one non-linear type of (ref) which is the so-called Generalized Auto-Regressive Conditional Heteroskedasticity (GARCH) model proposed by bollerslev1986generalized and has a form below:

equation[equation omitted — 152 chars of source]

where, $a \geq 0$, $a_1 > 0$, $b_1 > 0$, and $W_t\sim i.i.d.~N(0,1)$. The GARCH model is a generalization of the famous Autoregressive Conditional Heteroskedasticity (ARCH) model proposed by engle1982autoregressive. Its ability to forecast the absolute magnitude and quantiles or entire density of squared financial log-returns (i.e., equivalent to volatility forecasting to some extent) was shown by articleengle2001 using the Dow Jones Industrial Index. Later, many studies to investigate the performance of different GARCH-type models in predicting volatility of financial series were conducted; see following references peters2001estimating,gonzalez2004forecasting,lim2013comparing, herrera2018forecasting,karmakar2020bayesian. For ARCH/GARCH-type models, it is usual practice to identify and fit models based on quasi-maximum likelihood inference.

Traditionally, economists primarily utilize univariate GARCH-family models to understand dynamics of econometric data such as stock/ index/ price, etc observed for a long time. However, one of the key focuses of financial econometrics is to understand how extra knowledge such as fundamentals- and sentiments-based information could be beneficial to improve the prediction accuracy of market volatility if they are incorporated into the forecasting process; see more discussion from engle2007good and the references therein \footnote{In this regard, there is also a large literature that involves incorporating information of low-frequency variables using the GARCH-Mixed Data Sampling (MIDAS) model, as originally developed by engle2013stock.}. In line with the classical GARCH methods, we can wrap the exogenous covariates into the prediction process to get the GARCH models augmented with additional explanatory variables, GARCHX models in short. The estimation methodology of GARCHX models was discussed thoroughly in the work of francq2019qml; see more details about the GARCHX model in (ref).

Rather than taking the traditional approach (i.e., specifying and fitting a model and then predicting), we consider a recently proposed model-free prediction idea, namely normalizing and variance-stabilizing transformation (NoVaS transformation). The NoVaS method was initially developed by politis2003normalizing and then well discussed under the framework of the Model-free prediction principle in politis2015modelfreepredictionprinciple. In short, the Model-free prediction principle hinges on the idea of applying an inverse transformation function to bridge two equivalent probability spaces. For example, if we observe a univariate time series $\{Y_1,\ldots,Y_T\}$, we can try to find a transformation to map $\{Y_1,\ldots,Y_T\}$ to an $i.i.d.$ series $\{Z_1,\ldots,Z_T\}$. Since the prediction of $i.i.d.$ data is trivial, we can then transform the prediction of $i.i.d.$ data back to the prediction of the original data; see more details about the Model-free prediction principle in (ref).

Following the literature, we call this transformation-based approach a Model-free method in this paper. This necessitates us to emphasize that although the transformation is inspired from a model-assumption, for the prediction step it does not tie with any specific model structure. In the huge literature of applied econometrics, while analyzing data observed over a long time that can show signs of heteroscedasticity, the usual practice is to pick some specific GARCH-type model. Next an estimation of that model is carried out accordingly and subsequently the forecast of future volatility will be made based on the estimated model. Apparently, there is no universal rule for the choice of the specific GARCH model. In other words, a specific GARCH model can not work uniformly well across different datasets compared to other variants. On the other hand, the Model-free prediction could work well for any scenario as long as a suitable transformation function can be found. Therefore, the model-selection stage is shunned and the Model-free prediction approach is more robust against the model misspecification. Moreover, the standard GARCH methods require a relatively large sample size to be estimated well. For the NoVaS method, it tends to work stably even with short data. The existence of such transformation function in the context of predicting with exogenous variables will be analyzed theoretically in (ref).

Empirically speaking, NoVaS methods were mainly applied to forecast volatility in financial econometrics in the past few years. gulay2018comparison showed that the NoVaS method could beat GARCH-type models (GARCH, EGARCH and GJR-GARCH) with generalized error distributions by comparing the pseudo-out of sample (POOS) forecasting of volatility. Here the POOS forecasting analysis means using data up to and including the current time to predict future values. Later, chen2019optimal extended the NoVaS method to do multi-step ahead predictions. wu2021model further substantiated the great performance of NoVaS methods on time-aggregated long-term (30-steps ahead) predictions. wang2022model applied the Model-free idea to provide estimation and prediction inference for a general class of time series. Our present work is motivated by the wu2023model work where the authors recommended a so-called GARCH-NoVaS (GA-NoVaS) transformation structure inspired by the development of GARCH from ARCH. This NoVaS method is significantly robust against different model misspecification. Given this, it was a natural and probably quite an important question to see if such a robust forecasting framework can be built where exogenous covariates can be included and thus improve forecasting accuracy.

In this work, we explore the new methodology of forecasting stock market volatility with additional covariates being available to be included in the volatility dynamics. As far as we know, the NoVaS model-free prediction idea has not been studied when the exogenous variables are featured even in modeling the mean or average let alone the more complicated variance or volatility dynamics of a time-series. Due to the superior performance of the GARCH-NoVaS method in volatility forecasting, for this paper, we stick to the variance part and attempt to further boost the ability of the GA-NoVaS method with the help of exogenous covariate information. Towards this, we propose a so-called GARCHX-NoVaS (abbreviated as GAX-NoVaS henceforth) method which takes the GARCHX model as the starting step to build transformation. To obtain the inference about the future situation at an overall level, we choose the time-aggregated prediction metric. This aggregated metric has been applied to evaluate future predictions of electricity price or financial data fryzlewicz2008normalized,chudy2020long, karmakar2022long; see the formal definition in (ref). We wanted to check if the NoVaS prediction method can incorporate exogenous variables. More importantly, we hope the GAX-NoVaS method can sustain its great performance compared to the GARCHX method.

In addition to comparing our model-free method and classical GARCHX model with several simulated datasets, we also provide an interesting real data analysis. Our goal is to exhibit how our method could possibly shed light on the role of geopolitical risks, which are currently engulfing the global economy with multiple wars taking place, in forecasting the volatility of three stock markets of Europe: Germany and its two neighbors (Austria and Switzerland), based on a daily index of uncertainty associated with the Russia-Ukraine war as perceived by German Twitter activity. We also use newspaper-based metrics of global geopolitical risks due to acts and threats, as developed by caldara2022, to check for the robustness of our result covering a longer data sample. Hence, we add from a methodological perspective to the existing literature on forecasting international stock returns volatility using the information contained in geopolitical events and threats that basically rely on GARCH-type models (see, for example, Salisu2022,zhang2023 for details discussion of this literature). In this regard, note that, caldara2022 pointed out that entrepreneurs, market participants, and central bank officials view geopolitical risks as key determinants of investment decisions and stock market dynamics, with such risks, along with economic and policy uncertainties, forming an “uncertainty trinity” that would adversely impact the economy and the financial sector, as has been traditionally reported in the large existing literature on the impact of terror attacks and threats balcilar2018,bouras2018,bourforthcoming. In our empirical and simulation exercises, we measure the performance of GARCHX and GAX-NoVaS methods by the standard mean square prediction error. Moreover, we apply the forecast comparison tests to compare the two methods in a statistical way.

Our main contributions are summarized as follows:

itemize• We propose a new methodology-- namely GAX-NoVaS--to do the volatility forecasting with exogenous variables. This model-free method depends on a transformation function to connect two equivalent probability spaces instead of relying on any model assumption. The idea behind the GAX-NoVaS method hinges on the Model-free prediction principle. • We show such a transformation function exists under some mild conditions. This serves as the theoretical foundation of our method. • We apply our new method and standard GARCHX model to investigate the role of geopolitical risks in forecasting volatility. It turns out that our new method can be significantly more accurate, especially for long-horizon time aggregated predictions.

We organize the remainder of this article as follows. In (ref), we review the classical forecasting model, namely GARCHX, which is used as the starting point to propose the GAX-NoVaS method. Also, we present more details of the Model-free prediction principle and prove the existence of a transformation function with some exogenous variables existing. In (ref), we delineate the details of proposed GAX-NoVaS method. Then, some simulation studies and model evaluation criteria are collated in (ref). Next, we contrast our methods to existing classical ones on three empirical datasets in (ref). Finally, in (ref), we conclude by discussing the implications of our findings and some future directions.

GARCHX estimation and Model-free prediction principle

Before introducing our GAX-NoVaS method, we first explain the GARCHX model since the transformation of GAX-NoVaS is based on the GARCHX model. In addition, we give more details on the Model-free prediction principle and we prove the existence of a transformation function to achieve the Model-free prediction goal.

GARCHX model

In a seminal work, ARCH was proposed by engle1982autoregressive to model volatility or $\sigma_t^2$ for a time-series in a dynamic way. Following this, many different variants were developed in the econometrics literature. The GARCH model, especially the GARCH(1,1), stands as possibly the most popular one. The classic GARCH(1,1) model as defined by bollerslev1986generalized can be described below:

equation[equation omitted — 158 chars of source]

where $a \geq 0$, $a_1 > 0$, $b_1 > 0$, and $W_t\sim i.i.d.~N(0,1)$. More generally, we can express the GARCH-type models as

equation*[equation* omitted — 43 chars of source]

where a particular distribution of $\eta_t$ is not necessary and we can only assume that $\mathbb{E}(\eta_t^2|\mathcal{F}_{t-1}) = 1$. Usually, $\mathcal{F}_{t-1}$ is taken as the sigma-field generated by previous information $\{Y_{t-1},\ldots\}$. When additional information is available, people would like to utilize this extra knowledge to improve prediction accuracy. Subsequently, the so-called GARCHX model enters the public eye; see francq2019qml for discussions on the quasi-maximum likelihood estimation inference of GARCHX models. To simplify the analysis, we still assume the normality of $W_t$. After taking a vector of exogenous covariates $\bm{X} = (X_{1},\ldots,X_{m})$ into account, we can wrap the exogenous covariates into the prediction process by turning the GARCH(1,1) model into the following GARCHX(1,1,1) model:

equation[equation omitted — 182 chars of source]

where $\bm{X}_{t-1}$ represents $(X_{1,t-1},\ldots,X_{m,t-1})$ and $\bm{c}$ are the coefficients of these exogenous variables to be estimated. To perform a moving-window out-of-sample prediction experiment, we first need to estimate the GARCH(1,1) and GARCHX(1,1,1) models\footnote{For estimation of the GARCH and GARCHX models, we use the fGarch wuertz2013package and garchx packages sucarrat2020garchx in the R language and environment Rlanguage.}, and then we compute predictions iteratively; see (ref) for details. In this process, we assume that we know the true exogenous variables, which is feasible because we generate out-of-sample predictions. For practical applications, if needed, the future exogenous information can be estimated separately.

Model-free prediction principle

The model-free prediction principle was initially well developed by politis2015modelfreepredictionprinciple. Later, chen2019optimal applied this idea to multi-step ahead predictions of financial returns in the context of an ARCH-model structure. In short, the main idea behind the model-free prediction is to apply an invertible transformation function, $H_T$, that can map a non-$i.i.d.$ vector, $\{Y_t~;t = 1,\ldots,T\}$, to a vector, $\{\epsilon_t;~t=1,\ldots,T\}$, with $i.i.d.$ components (chosen as standard normal in this work). Due to the invertibility of the function, $H_T$, it is possible to construct a one-to-one relationship between a future value, $Y_{T+1}$, and $\epsilon_{T+1}$, i.e.,

equation[equation omitted — 77 chars of source]

where $\bm{Y}_{T}$ denotes all historical data $\{Y_t;~t =1,\ldots,T\}$; $\bm{X}_{T+1}$ is the collection of all predictors, and it also contains the value of a future predictor $X_{T+1}$; the form of $f_{T+1}(\cdot)$ depends on $H^{-1}_{T}$. This relationship implies that we can also transform the prediction of $\epsilon_{T+1}$ to the prediction of $Y_{T+1}$. Assume we have $\hat{\epsilon}_{T+1}$ to the be the predictor of $\epsilon_{T+1}$, we can express the predictor of $Y_{T+1}$ as

equation[equation omitted — 94 chars of source]

Because the prediction of $i.i.d.$ data is standard, the $L_1$ (Mean Absolute Deviation) , $L_2$ (Mean Squared Error), or another optimal quantile predictor of $\epsilon_{T+1}$ can easily be found. We, thus, can easily obtain the corresponding optimal predictor of $Y_{T+1}$.

For multi-step ($h$-step) ahead prediction, we simply repeat this prediction process, i.e., we express $Y_{T+h}$ through a function w.r.t. $\bm{Y}_{T}$, $\bm{X}_{T+1}$ and $\{\epsilon_{T+1},\ldots,\epsilon_{T+h}\}$:

equation[equation omitted — 123 chars of source]

In order to compute the prediction of $Y_{T+h}$, we take a distribution-match approach to approximate the distribution of $Y_{T+h}$. Ideally, when we know the exact distribution of the $i.i.d.$ $\epsilon$, we can use a Monte Carlo simulation to approximate the distribution of $Y_{T+h}$ based on (ref). Practically speaking, when we just have the empirical transformation results, i.e., the observed sample $\{\epsilon_t\}_{t=1}^{T}$, bootstrap is an appropriate approach. Moreover, we can even predict $g(Y_{T+h})$, where $g(\cdot)$ is a general continuous function. For example, we can compute the $L_1$ and $L_2$ optimal predictors of $g(Y_{T+h})$ as below:

equation[equation omitted — 330 chars of source]

where $g(Y_{T+h})_{L_2}$ and $g(Y_{T+h})_{L_1}$ represent the optimal $L_2$ and $L_1$ predictor of $g(Y_{T+1})$, the $\{\hat{\epsilon}_{T+1,m}\}_{m=1}^{M}$ are generated by bootstrap or Monte Carlo simulation, and $M$ is some large number (2000 in our empirical analysis). For further discussion, see politis2015modelfreepredictionprinciple.

To the best of our knowledge, the model-free prediction idea has not been studied when the model features exogenous variables. However, it may be beneficial to take into account in the prediction process the additional information embedded in such exogenous variables. To show the NoVaS approach is still applicable, we need a transformation function that maps the targeted variables and exogenous predictors together into some simple $i.i.d. $ random variables. Under some mild conditions, we show the existence of such a transformation function based on the probability integral transform. We assume:

itemize• A1 The joint density of $\{Y_1,\cdots, Y_T \}$ exists for any $T\geq 1$. • A2 For exogenous random vector $\bm{X}: = \{X_1,\cdots, X_{m}\}$, the joint density $\{Y_1,\cdots, Y_T, X_1, \cdots, X_m \}$ exists for any $m\geq 1$.

Then, the feasibility of NoVaS transformation with exogenous variables existing is guaranteed by (ref) shown below:

TheoremUnder A1 and A2, there exists a function $\bm{g}$ such that $\bm{Z} = \bm{g}((\bm{Y},\bm{X}))$ and the corresponding inverse function $\bm{h}$ such that $(\widetilde{\bm{Y}},\widetilde{\bm{X}}) = \bm{h}(\bm{Z})$; $\bm{Z} \sim N(0,\bm{I}_{T+m})$; $\bm{Y} = (Y_1,\cdots,Y_T)$ and $\bm{X} = (X_1,\cdots,X_m)$ are any two random vectors; $(\widetilde{\bm{Y}},\widetilde{\bm{X}})$ have the same joint distribution of $(\bm{Y},\bm{X})$.
proofThe proof of (ref) is based on the probability integral transform; see angus1994probability for a review. Without loss of generality, we start from $Y_1$ to determine the transformation function $\bm{g}$. Let $U_1 := \Tilde{g}_1(Y) = F(Y_1)$; $F(Y_1)$ is the distribution of $Y_1$. According to the probability integral transform, we know $U_1$ has a uniform distribution on $[0,1]$. Then, we make \begin{equation} U_2 := \Tilde{g}_2(Y_1,Y_2) = F(Y_2 | Y_1 ). \end{equation} $F(Y_2 | Y_1 )$ is the conditional distribution of $Y_2$. (ref) implies that $Z_2$ is $\text{Uniform}(0,1)$ conditional on $Y_1$. Thus, the unconditional (marginal) distribution of $U_2$ is still $\text{Uniform}(0,1)$, and $U_2$ and $Y_1$ are independent so that $U_2$ is also independent with $U_1$. This can be seen from the equation below: \begin{equation} p_{U_2, Y_1}(u_2,y_1) = p_{U_2|Y_1}(u_2|y_1) p_{Y_1}(y_1). \end{equation} Integrating both sides w.r.t. $y_1$, we can find $p_{U_2}(u_2) = 1$ on the region $[0,1]$, since $Z_2$ is $\text{Uniform}(0,1)$ conditional on $Y_1 = y_1$ for any $y_1$. We can repeat this process as a Gram–Schmidt-like recursion, i.e., we let $U_3 := \Tilde{g}_3(Y_1,Y_2,Y_3) = F(Y_2 | Y_1,Y_2 )$ and so on. In total, we need $\{\Tilde{g}_1,\cdots,\Tilde{g}_{T+m}\}$ and they are functions of $\bm{Y}$. Thus, there exists a function $\Tilde{\bm{g}}$ which maps $(\bm{Y},\bm{X})$ to $\bm{U}$ which has $i.i.d.$ uniform components $\{U_1,\cdots, U_{m+T}\}$, i.e., $\bm{U} = \Tilde{\bm{g}}( (\bm{Y},\bm{X}))$. Then, $\bm{Z} = \Phi^{-1}\circ\Tilde{\bm{g}}((\bm{Y},\bm{X}))$ has multivariate normal distribution $N(0,\bm{I}_{T+m})$; $\Phi^{-1}$ is the quantile function of $N(0,\bm{I}_{T+m})$. Finally, we can take $\bm{g} = \Phi^{-1}\circ\Tilde{\bm{g}}$. On the other hand, if $F^{-1}_{Y_1}:=\inf \{x: F_{Y_1}(x) \geq y\}, 0\leq y\leq 1$, then $F_{Y_1}^{-1}(U_1)$ has the distribution as the same as $Y_1$. Similarly, we can get the conditional distribution of $Y_2$ on $Y_1$ by taking $F_{Y_2|Y_1}^{-1}(U_2)$. By repeating this process, we can recover the joint distribution of $(\bm{Y},\bm{X})$ by chain rule. In other words, there exists a $\bm{h}$ such that $(\bm{Y},\bm{X}) = \bm{h}(\bm{Z})$.

The direct implication of (ref) is that the Model-free prediction principle is feasible even when exogenous variables are included in the dependence dynamics. Moreover, our theorem is more general than the result from wang2022model where the time series must satisfy some strict conditions. In short, we build a transformation function based on the GARCHX model structure to estimate the oracle functions $\bm{g}$ and $\bm{h}$, so we call our method GAX-NoVaS. We should mention again that the “model-free” in this context means we do not rely on an assumed underlying model to make predictions. Although a transformation function needs to be estimated, it is just a “bridge” that connects original and transformed distributions according to the distribution match idea.

GAX-NoVaS model-free prediction method

We first present the state-of-the-art GARCH-NoVaS method which is based on a so-called NoVaS transformation. Then, we extend the GARCH-NoVaS method to a GAX-NoVaS method, which features the exogenous variables.

NoVaS transformation

For the sake of completeness, we first give a brief introduction to the NoVaS transformation (model) which is a direct application of the Model-free prediction idea explained in (ref). Initially, the NoVaS transformation is developed from the ARCH model:

equation[equation omitted — 77 chars of source]

here, these parameters satisfy $a\geq 0$, $a_i\geq 0$, for all $i = 1,\ldots,p$; $W_t\sim i.i.d.~N(0,1)$. In other words, the structure of the ARCH model gives us a ready-made $H^{-1}_T$. We can express $W_t$ in (ref) using other terms to get a potential $H_{T}$ :

equation[equation omitted — 116 chars of source]

politis2003normalizing further modified (ref) as follows:

equation[equation omitted — 144 chars of source]

here, $\{Y_t;~t=1,\ldots,T\}$ is the sample data; $\{W_{t};~t=p+1,\ldots,T\}$ is the transformed vector; $\alpha$ is a fixed scale invariant constant; $s_{t-1}^2$ is an estimator of the variance of $\{Y_i;~i = 1,\ldots,t-1\}$ and can be calculated by $(t-1)^{-1}\sum_{i=1}^{t-1}(Y_i-\overline{Y})^2$, where $\overline{Y}$ is the sample mean of $\{Y_i;~i = 1,\ldots,t-1\}$. For making (ref) be a qualified function $H_T$, i.e., making $\{W_t\}_{t=p+1}^{T}$ really obey $i.i.d.$ standard normal distribution, we need to impose some restrictions on $\alpha$ and $\beta, a_1,\ldots,a_p$. We first stabilize the variance by requiring:

equation[equation omitted — 128 chars of source]

In application, $\{W_t\}_{t=p+1}^{T}$ transformed from financial log-returns by NoVaS transformation are usually uncorrelated. Therefore, if we make $\{W_t\}_{t=p+1}^{T}$ close to a Gaussian series i.e., normalizing $\{W_t\}_{t=p+1}^{T}$, we can get the desired $i.i.d.$ property. This is why this transformation is called NoVaS.

There are many criteria to measure the normality of a series. Under the observation that the distribution of financial log-returns is usually symmetric, we choose the kurtosis to be a simple distance to measure the departure of a non-skewed dataset from that of the standard normal distribution politis2015modelfreepredictionprinciple. Besides, matching marginal distribution seems sufficient to normalize the joint distribution of $\{W_t\}_{t=p+1}^{T}$ for practical purposes based on empirical results. If we denote the marginal distribution of $\{W_t\}_{t=p+1}^{T}$ and the corresponding kurtosis by $\widehat{F}_w$ and $\text{KURT}(W_t)$, respectively, we then attempt to minimize $|\text{KURT}(W_t)-3|$ to obtain the optimal combination of $\alpha,\beta,a_1,\ldots,a_p$ such that $\widehat{F}_w$ is as close to standard normal distribution as possible. Subsequently, the NoVaS transformation can be determined.

The remaining difficulty is how to finish this optimization step to get optimal coefficients $\alpha,\beta,a_1,\ldots,a_p$, especially when $p$ is large. To simplify this problem, politis2015modelfreepredictionprinciple defined an exponentially decayed form of $\{a_i\}_{i=1}^p$:

equation[equation omitted — 150 chars of source]

The NoVaS transformation based on coefficients defined in (ref) is called Generalized Exponential NoVaS (GE-NoVaS). In other words, we can represent the $p+2$ number of coefficients by two parameters $c$ and $\alpha$, which relieve the optimization burden, but with a sacrifice that the coefficients are fixed in a decayed form. To achieve a balance between the relief of the optimization dilemma and the freedom of coefficients, inspired by the development of GARCH from ARCH, wu2023model built a NoVaS transformation according to the GARCH model, namely GARCH-NoVaS which was shown to be more stable and accurate. Later, we specify the GAX-NoVaS transformation in detail.

GAX-NoVaS transformation method

Starting from (ref), we take similar steps of building GA-NoVaS to find the transformation function of the GAX-NoVaS method. In order to simplify the notation, we consider the case of only one exogenous covariate $X_t$. The case of multiple exogenous covariates can be analyzed in an analogous way. First, we notice that we can rewrite the (ref) as

equation[equation omitted — 117 chars of source]

We also have

equation[equation omitted — 213 chars of source]

so that we can substitute these terms into (ref). We then get

equation[equation omitted — 141 chars of source]

where $p$ is a large constant that is used to truncate the infinite summation (we take $p=q$ in this work), because $a_1$ and $b_1$ need to be less than one to guarantee the stationary property of the GARCHX series. In line with the NoVaS transformation, we finally write the transformation function as follows:

equation[equation omitted — 174 chars of source]

where $s^2_{t-1,Y}$ and $s^2_{t-1,X}$ are the sample variance of $\{Y_1,\ldots,Y_{t-1}\}$ and $\{X_1,\ldots,X_{t-1}\}$, respectively. Thus, we can use (ref) as the transformation function for the GARCHX model, where $\{W_t\}$ is the transformed series. We want to make $\{W_t\}$ $i.i.d.$ normal, so that the one-step (conditional) prediction $\widehat{Y}_{T+1}$ can be expressed as

equation[equation omitted — 192 chars of source]

where $\widehat{W}_{T+1}$ is the optimal point prediction based on $i.i.d.$ $\{W_1,\ldots, W_{T}\}$. Multi-step-ahead predictions can be computed as explained in (ref).

The final question left now is how to find a transformation function that indeed makes $\{W_t\}$ $i.i.d.$ normal. While we have made some brief remarks on this question in (ref), we next provide a full explanation with a focus on the GAX-NoVaS method. Our goal is to determine the coefficients, $\alpha,\beta,a_1,b_1,c_1$, of (ref) to obtain the desired transformation. The most important step is to minimize $|\text{KURT}(W_t) - 3|$, where $3$ is the kurtosis of normal distribution. For this optimization, we use the numerical technique to find the optimal coefficients.\footnote{We use the nloptr package ypma2014nloptr for R.} In operation, one may obtain some extremely large values from the transformed series $\{W_t\}$, and such outliers may spoil the normality of the transformed series and may also influence prediction performance. Thus, before moving on to the prediction step, we truncate the transformed series by the 0.99 and 0.01 quantile values of a normal distribution with mean and standard deviation given by the sample mean and sample standard deviation of $\{W_t\}$.

RemarkIt is not difficult to perceive that the transformed series $\{W_t\}$ may be correlated, such that the decorrelation step is beneficial and necessary. One way to carry out this step is by fitting an AR(p) model on the $\{W_t\}$ series. Then, we record the residuals of the AR fit as $\{\hat{\epsilon}_t\}$. Also, we can approximate the one-step-ahead value $\widehat{W}_{T+1}$ with a fitted AR model. Then, we can create the new series as $\{\epsilon_{t} + \widehat{W}_{T+1} \}$. We use the empirical distribution of this new series to approximate the distribution of $W_{T+1}$. (ref) can be used with the optimal prediction $\widehat{W}_{T+1}$ derived from this empirical distribution. It is still an open question, however, how to extend this decorrelation step to multi-step-ahead predictions.

Simulations

In this section, we deploy several simulations to check the performance of GARCH, GARCHX and GAX-NoVaS methods. Before presenting the data-generating model used to do simulations, we explain the procedure of the moving-window time-aggregated predictions and give the model evaluation metrics to measure the performance of different methods.

Moving-window time-aggregated prediction

If we have sample $\{Y_1,\ldots,Y_{N}\}$ at hand, in order to fully exhaust the dataset, we can focus on moving-window out-of-sample predictions, i.e., we use $\{Y_1,\cdots,Y_{T}\}$ to predict \{$Y_{T+1}^2,\cdots,Y_{T+h}^2\}$, then we use $\{Y_2,\cdots,Y_{T+1}\}$ to predict $\{Y_{T+2}^2,\cdots,Y_{T+h+1}^2\}$, and so on until we reach the end of the sample (that is, until we use $\{Y_{N-T+h+1},\cdots,Y_{N-h}\}$ to predict $\{Y^2_{N-h+1},\ldots, Y^2_{N}\}$). Here, $T$ denotes the moving window size; we fix its size as 250; $h$ is the prediction horizon, i.e., 1, 5, 20 in our setting. Sometimes, we may not have enough data available to perform predictions. Thus, the window size $T = 250$ is designed to see if this method is stable and can still return accurate predictions even with short data. A 250-size moving window is in line with around one year of daily financial data. In this perspective, it is important to keep in mind that a 250-size moving window is practically meaningful since the time series may not be stationary on a wider span.

In addition, we are not only interested in the one-step-ahead prediction $h=1$ but also multi-step-ahead prediction $h>1$. From a practical aspect of forecasting volatility, as mentioned above in our introduction Section (ref), the long-term prediction ($h$ takes a large value) is important and can guide future strategic decisions. Before this evaluation, we start by writing time-aggregated predictions as follows:

equation[equation omitted — 113 chars of source]

where $\overline{\widehat{Y}}_{T,h}^2$ is the $h$-step ahead time-aggregated volatility prediction starting from $Y_T$. For example, if the total number of data $N = 1000$ and we consider the 6-step-ahead moving window time-aggregated predictions with $T = 500,$ we need to find predictions $\overline{\widehat{Y}}_{T,6}^2$ for $T = 500,\ldots,994$.

We hope the time-aggregated volatility prediction is close to the true aggregated value calculated from the realized average squared log-returns $\overline{Y}_{l,h}^2 = \sum_{k=1}^h(Y_{T+k}^2/h)$. To evaluate the accuracy, we can consider the specific mean of squared prediction errors (MSPE) shown below, with this statistic aiming to compare the prediction performance in an absolute way:

equation[equation omitted — 129 chars of source]

where $\overline{\widehat{Y}}_{l,h}^2$ and $\overline{Y}_{l,h}^2$ denote the predicted and true time-aggregated values for each moving-window forecasting, respectively.

Simulation setting

In this part, we present three data-generating models to simulate data and then evaluate the performance of various methods considered in this paper. Due to the true data-generating process being known to us, we can simulate any size of the sample and compare predictions from various methods with oracle values. Besides, we remove the variance term $\beta s^2_{T-1,X}$ of (ref) when we do the transformations since the prediction performance hardly changes with or without $\beta s^2_{T-1,X}$. Three true underlying models are presented below:

itemize• Model 1: Standard GARCH(1,1) with Student-$t$ errors\\ $Y_t = \sigma_t\epsilon_t,$ $~\sigma_t^2 = 0.00001 + 0.73\sigma_{t-1}^2+0.1Y_{t-1}^2 + c|X_{t-1}|,$\\ $X_{t-1} \sim i.i.d.~N(0,1)$;$~\{\epsilon_t\}\sim i.i.d.~t$ $\text{distribution with four degrees of freedom}$\\ $c = 1$. • Model 2: Standard GARCH(1,1) with Student-$t$ errors\\ $Y_t = \sigma_t\epsilon_t,$ $~\sigma_t^2 = 0.00001 + 0.8895\sigma_{t-1}^2+0.1Y_{t-1}^2 + c|X_{t-1}|,$\\ $X_{t-1} \sim i.i.d.~N(0,1)$;$~\{\epsilon_t\}\sim i.i.d.~t$ $\text{distribution with four degrees of freedom}$\\ $c = 1$. • Model 3: Time-varying GARCHX(1,1) with standard normal exogenous variables:\\ $Y_t = \sigma_t\epsilon_t,~\sigma_t^2 = b_{t}\sigma_{t-1}^2+a_{t}Y_{t-1}^2 + c|X_{t-1}|,$\\ $X_{t-1} \sim i.i.d.~N(0,1)$; $\{\epsilon_t\} \sim i.i.d.$ $t$ distribution with five degrees freedom;\\ $c = 1; g_t = t/n$; $a_{t} = 0.1 - 0.05g_t$; $b_{t} = 0.7 + 0.2g_t,~n$ is the total length of the time series.

To check the robustness of methods on model misspecification, we intend to take the innovation distribution of simulation models as the $t$-distribution. Besides this purpose, we argue that the $t$ distribution as the innovation to mimic the real-world cases is more appropriate since real data usually show the heavy tail phenomenon. Models 1 and 2 are from a standard GARCH where in Model 2 we intended to explore a scenario that $\alpha_1 + \beta_1$ is very close to 1 and thus mimics what would happen for the iGARCH situation. In addition, we make the coefficients of the GARCHX model change linearly in Mode-3 so that we can observe the ability of different methods to handle the data generated from a time-varying model which is more coherent to the real-world situation. To sync with the empirical studies later, we simulate a time series with a length $T = 4694$. We take the moving-window size $T = 250$. MSPE of GARCH, GARCHX and GAX-NoVaS methods with three simulation settings are presented below:

table[table omitted — 1,070 chars of source]

From (ref), it is clear that the GAX-NoVaS is much better than GARCH and even the classical GARCHX models according to the MSPE criterion. Interestingly, the GARCHX model is even much worse than the GARCH model for some specific cases, e.g., the 20-step-ahead prediction of Model-3. By taking a deeper analysis, we find the terrible performance GARCHX method is due to some extremely large predictions. On the other hand, the prediction returned by GAX-NoVaS is more stable. The superiority is further shown in (ref) with three real datasets and various exogenous predictors.

Since the focus of this paper is exploring a new approach to incorporate exogenous variables, we take the DM test to evaluate the performance of GARCHX and GAX-NoVaS methods more formally; see diebold2002comparing for the technical details of the DM-test\footnote{We perform the DM-test with the function DM-test in the R package multDM.}. The DM-test results on comparing GARCHX and GAX-NoVaS for forecasting three simulated datasets are tabularized in (ref). These tests further verify the advantage of our methods on forecasting with exogenous variables, especially for a long-prediction horizon.

table[table omitted — 958 chars of source]

\FloatBarrier

Empirical analyses with real data

A summarizing note of our findings in (ref) reads that the GAX-NoVaS method performs better than standard GARCH-type methods, especially for long-term time aggregated predictions. In this section, we deploy an interesting data analysis to exhibit how our method could shed light on the role of geopolitical risks in forecasting volatility with real-world data. We start by describing the data below.

Data description

The ongoing Ukraine-Russia has led many countries in Europe, particularly Germany, to adjust their military and security, as well as energy policies in light of new geopolitical risks, with such adjustments entailing large costs to the macroeconomy and financial markets, as depicted by grebe2024. In this regard, these authors, first, assemble a data set of more than eight million German Twitter posts related to the war in Ukraine to construct a daily index of uncertainty about the war as perceived by German Twitter based on using state-of-the-art methods of textual analysis. grebe2024 show that an increase in uncertainty has strong effects on financial markets, associated with a significant decline in economic activity as well as an increase in expected inflation. We utilize this index (Ukraine)\footnote{The data is available for download from: \url{https://www.uni-giessen.de/de/fbz/fb02/fb/professuren/vwl/tillmann/forschung/ukraine-uncertainty-index}.} in our empirical analysis to forecast stock market volatility of not only Germany, but two of its neighbors namely, Austria and Switzerland, over the daily period of 1st January, 2021 to 28th February, 2023. The national stock market indexes (ATX (Austria), DAX (Germany), SMI (Switzerland)) of these three countries, for which we compute log-returns to feed into our volatility models were derived from the Bloomberg terminal. With the focus being on geopolitical risks, we also utilized the daily newspapers-based geopolitical risks index (GPRD) of caldara2022\footnote{The data can be accessed from: \url{https://www.matteoiacoviello.com/gpr.htm}.}, which, in turn, allowed us to analyze a longer data sample covering 2nd January, 2006 to 10th August, 2023. The starting date of this longer sample, and the choice of these three countries, were also motivated by the availability of Google searches-based daily data on economic activity (Trend) for all three countries, and inflation (Inflation) for Germany and Switzerland,\footnote{The data can be downloaded from: \url{https://www.trendecon.org/}.} which are used as additional predictors to ensure that our results are not only limited to geopolitical risks.

We present three log-return series of Germany, Switzerland and Austria from 2nd January, 2006 to 10th August, 2023 in (ref). The volatility clustering phenomena observed in all plots reveals the heteroskedasticity within these three series. To investigate the property of three long return series, we provide the summary statistics, e.g., the mean, skewness, and kurtosis in (ref). To verify the heteroskedasticity with all series more directly, we split the whole time period into four equal-length sub-periods and denote the sample variance of all four sub-periods by $V_i$, $i = 1,\ldots, 4$. These statistics are also provided in (ref). Towards statistical tests, we also perform modified Ljung-Box (m-LB) and ARCH Lagrange Multiplier (ALM) tests to check the autocorrelation and ARCH effects of squared return series. For the m-LB test, we consider the lag order 20. For the ALM test, we consider the maximum lag order 10. These two tests are performed in R with functions lbtest and Lm.test, respectively. The p-values of tests are presented. Summarizing (ref), the large kurtosis values indicate the heavy-tailed property for all three log-return series. The variance of return series in different time regions changes notably, indicating the heteroskedasticity. The ALM test with a pretty small p-value also confirms the heteroskedasticity for all return series. The m-LB test shows strong evidence of autocorrelation within all squared return series.

table[table omitted — 1,143 chars of source]

To show the fluctuations behind the Ukraine index and GPRD index, we present two plots in (ref) below. As one can see from there, these two indices fluctuate severely around the beginning of 2022 which corresponds with the real-world event. Later, we attempt to use this information to forecast the stock market volatility of three countries.

figure[figure omitted — 210 chars of source]
figure[figure omitted — 222 chars of source]
figure[figure omitted — 210 chars of source]
figure[figure omitted — 251 chars of source]
figure[figure omitted — 181 chars of source]

\FloatBarrier

Empirical results

We first consider the forecasting exercise with the short period (1st January, 2021 to 28th February, 2023) data described in (ref). Then, the analysis of three long returns series is given in (ref).

Short data

To compare the performance of GARCH, GARCHX and GAX-NoVaS methods, we still apply the time aggregated prediction metric described in (ref). We consider $h = 1,5, 20$ and use a 250-size moving window, which is about 1 year of daily data. We start the empirical analysis with the short period of data and we take the Ukraine index as the exogenous predictor. The MSPE results are summarized in (ref).

table[table omitted — 1,070 chars of source]

From (ref), we can see that the GAX-NoVaS method can bring some large improvements, especially for long-horizon aggregated predictions. Meanwhile, it seems that the GARCHX and GARCH models have indistinguishable performance. However, the GAX-NoVaS method is generally better than both GARCH-type methods. The DM-test results on comparing GARCHX and GAX-NoVaS for forecasting short real-world data are tabularized in (ref), which reveals the significant advantage of GAX-NoVaS for long-horizon predictions, especially for the 20-step-ahead predictions of Short Austria data.

table[table omitted — 993 chars of source]

Long data

We continue our real data analysis with long datasets (2nd January, 2006 to 10th August, 2023). We also apply more exogenous variables. Similar to the analysis procedure for short data, we present MSPE ratios and corresponding DM-test results of GARCHX and GAX-NoVaS in (ref). Generally speaking, the GAX-NoVaS method still dominates the other two GARCH-type methods, and this superiority is verified to be significant by the DM-test.

table[table omitted — 1,541 chars of source]
table[table omitted — 1,405 chars of source]

\FloatBarrier

Conclusion

We extend the current NoVaS prediction method to the realm of prediction with exogenous variables. We provide the theoretical foundation to guarantee the feasibility of applying Model-free prediction. Inspired by the GARCHX model, we propose a specifically designed model-free/model-based method namely GAX-NoVaS prediction. The dominance of GAX-NoVaS on the classical GAX-NoVaS method is verified by simulation and empirical datasets. Also, such an advantage is not only exhibited by some MSPE metric but we also show some statistical significance through the parlance of classical DM tests.

We should also mention that going far beyond GAX-NoVaS method might have limitations if the model becomes increasingly complex. Recall that the great performance of the GAX-NoVaS method relies on a successful transformation, the satisfied transformation may not be achievable if the underlying time series is very complicated. However, there is a growing literature on forecasting using more non-parametric neural network based models. It will be an interesting future work that combines the idea of model-free prediction with the state-of-the-art machine learning method, such as Deep neural network (DNN), convolutional neural network (CNN) or LSTM etc. Finally, in the field of binary/categorical/count data INGARCH models have recently garnered significant attention both from theoretical and applied researchers. One could potentially also think of a Model-free INGARCH-X type model and try to integrate the idea of exogenous covariates into it and build a new forecasting framework to challenge the existing ones. In short, our paper remains the first paper to propose this idea of model-free transformation-based forecasting focused on GARCHX type models but the scope of extending this to several directions is ample.