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.
60,386 characters · 13 sections · 100 citation commands
Generalized Autoregressive Score Trees and Forests
\doublespacing
Keywords: Forecasting, machine learning, random forest, regression tree, volatility, copula, durations.
J.E.L. Codes: C22, C32, C53.
Models for economic time series data that can capture time variation in features of the predictive density are widely used for policy making, investment decisions, risk management, and in many other applications. Such models include the autoregressive-moving average model of box1970time, the ARCH/GARCH models of engle1982autoregressive and bollerslev1986generalized, and many others. The family of “generalized autoregressive score” (GAS) models, proposed by creal2013generalized and harvey2013dynamic, nests these time series models and others, and has been applied to a wide range of problems. artemova2022scoreTheory,artemova2022scoreApplications and harvey2022score provide recent surveys of this large and growing literature.
Despite their success, score-driven models are inevitably only approximations to the true data generating process. We propose to use data mining methods from the machine learning literature to improve the performance of these models. Specifically, we propose a “GAS tree,” that combines the parsimonious structure of the GAS model with the flexible, data-driven learning of decision trees breiman1984classification,breiman2017classification. A GAS tree allows the parameters of the model to vary across “branches” of the tree, which are formed using a possibly large collection of state variables. This leads to a model that can incorporate information from outside the GAS model, and that allows for potentially complicated nonlinearities and interactions. We further propose “GAS forests,” analogous to the “random forests” of breiman2001random for linear regression, where we create many GAS trees using bootstrap samples of the original data and then average the forecasts from these trees. In many applications random forests have been found to improve upon regression trees due to the reduction in variance obtained via averaging, see e.g. hastie2009elements.
The estimation of GAS trees and GAS forests is computationally demanding. It involves finding the optimal state variables and thresholds from the set of candidate variables, as well as estimating the parameters of the GAS model. We use cluster computing and a “greedy” estimation algorithm related to that of breiman1984classification for regression and audrino2001tree for GARCH trees. This algorithm finds a near-optimal solution and converges quickly. A key hyper-parameter in tree and forest models is the maximum depth of the tree (essentially, how many subsamples of the data will be considered) and we tune this parameter using a validation sample, separate from our forecast evaluation sample.
We apply the proposed GAS tree and GAS forest models in four empirically relevant problems: forecasting stock return volatility, the distribution of stock returns, the joint distribution of stock and bond returns, and high-frequency trade durations. As baseline models for these applications we use the GARCH model of bollerslev1986generalized, the t-GAS model of creal2011dynamic, a joint distribution model with Student's $t$ margins and a Student's $t$ copula, as in janus2014long, and the ACD model of engle1998autoregressive. We then consider tree and forest extensions of these models, and in all four cases we find that the baseline model is significantly out-performed. For the two stock return applications, we find that the GAS tree provides the best out-of-sample forecasts. The estimated tree structures provide significantly better forecasts, and turn out to be relatively simple: we find evidence of a leverage effect, where the GAS model parameters differ depending on whether the lagged stock return was positive or negative, and a variance risk premium effect, where the model parameters differ depending on whether the difference between option-implied and historical volatilities is large or small.
In our study of the joint predictive distribution of stock and bond returns, we find that the GAS forest produces the best out-of-sample forecasts. Variable importance analyses indicate that the most important variables for the GAS forest are the lagged stock and bond returns themselves, indicating omitted nonlinearity in the baseline GAS model. We find evidence of a flight-to-quality effect, where higher bond returns or lower stock returns are associated with even more negative long-run correlations between the stock and bond markets. In our analysis of trade durations, defined as the time taken for 10,000 shares of the S&P 500 exchange traded fund, SPY, to be transacted, we again find the forest-based extension to be the preferred model. In this application the most important state variables are both measures of volatility, consistent with the well-known volume-volatility relationship (see, e.g., Karpoff87).
This paper is part of the fast-growing literature using tools from machine learning in econometrics, see varian2014big and athey2019machine for recent surveys. Various studies have found that machine learning techniques bring significant gains over traditional econometric methods for forecasting applications. For example, medeiros2021forecasting, goulet2020macroeconomy and huber2020nowcasting show that tree-based methods, including random forests, can produce more accurate forecasts of important macroeconomic variables like unemployment and inflation. gu2020empirical and bianchi2021bond show how machine learning methods can improve forecasts of stock and bond returns.
In addition to macroeconomic and financial forecasting, some recent papers have found success applying machine learning methods to volatility models such as the GARCH model of bollerslev1986generalized and the HAR model of corsi2009simple. For instance, christensen2021machine shows that neural networks and random forests significantly improve over HAR model, and nguyen2022statistical,nguyen2022recurrent create hybrid stochastic volatility and GARCH models with recurrent neural networks. reisenhofer2022harnet and tetereva2022forest use convolutional neural networks and random forests, respectively, combined with the HAR model to obtain improved out-of-sample forecasts.
This study also relates to a broadly defined “local estimation” literature. tibshirani1987local, fan1998local and fan2009local use kernel-based methods to localize (quasi-) maximum likelihood models. A more recent strand of this literature includes breiman2001random, schlosser2019distributional and athey2019generalized, who use decision trees and random forests to localize regressions, parametric distributions, and GMM models respectively. Our paper is related to oh2021better, which is part of the first strand of this literature. That paper's approach suffers from the curse of dimensionality, due to its use of kernel-based methods, and it additionally requires that all (or none) parameters of the baseline model are localized. In the next section, we show that our proposed approach can deal with a large number of state variables and permits a subset of parameters to be localized, allowing the researcher to impose more or less structure on the model as needed.
The remainder of the paper is structured as follows. In Section (ref) we review the class of generalized autoregressive score models of creal2013generalized and harvey2013dynamic and introduce our new GAS tree and GAS forest models. Section (ref) also includes computational details on the implementation of these models. Section (ref) presents four empirical analyses, applying the new methods to forecasting volatility, correlation, and univariate and bivariate distributions. Section (ref) concludes, and the appendix presents details on the derivations for the third application. A supplemental appendix contains additional results.
The class of generalized autoregressive score (GAS) models of creal2013generalized and harvey2013dynamic provide a parsimonious and powerful way to capture time variation in the parameter(s) of a given probability density function. We describe this model below, and in Sections (ref) and (ref) we introduce tree- and forest-based extensions of this class of models.
Let the dependent variable be denoted $\mathbf{y}_t \in \mathbb{R}^K$. Conditional on the information set $\mathcal{F}_t$, this variable is assumed to have a parametric predictive density $p$, with $d$-dimensional time-varying parameter $f_t$, and potentially a static parameter $\nu$. The GAS($p,q$) model specifies the evolution of $f_t$ as:
It is the appearance of the score, $\nabla_t$, in the evolution equation for $f_t$ that gives this class of models its name.\footnote{We follow creal2011dynamic and use the inverse information matrix to scale the score in all of our applications, though other choices for this matrix are possible, such as the square root of this matrix, or simply the identity matrix.} Similar to the well-known Newton-Raphson algorithm for numerical optimization, at each date $t$, $f_t$ moves in the direction that most improves the model fit.
Let $\theta =(\omega, \mathrm{vec}(B_1), ..., \mathrm{vec}(B_q), \mathrm{vec}(A_1), ..., \mathrm{vec}(A_p))$ denote the vector of all GAS parameters of this model, making $(\theta,\nu)$ the full set of unknown parameters. Since GAS models are “observation driven,” as opposed to “parameter-driven”, the likelihood function is available in closed form, and $(\theta,\nu)$ can be estimated by maximum likelihood with low computational cost. This feature makes it feasible to consider tree- and forest-based extensions of this class of models, which we introduce below.
Regression trees (breiman1984classification; breiman2017classification) are a type of nonparametric regression based on sequentially splitting the available data into partitions. The partitions are formed using one or more state variables, $Z_t$, and estimated threshold value(s), $c$. Figure (ref) illustrates a simple tree structure. The left panel shows a tree with two state variables and specific thresholds, and the right panel shows the corresponding partition of the support of state variables. This hypothetical tree has three “terminal nodes” and implies a specific partition of the data, denoted $\boldsymbol{\mathcal{P}}=\{\mathcal{P}_1,\mathcal{P}_2,\mathcal{P}_3\}$. Given a tree structure, a “regression tree” is obtained by estimating a linear regression separately for each of the terminal nodes in the tree. In so doing, regression trees allow for nonlinearities and multi-way interactions, greatly generalizing the baseline regression model. Naturally this flexibility makes trees prone to overfit the training data, and therefore, trees must be regularized, or “pruned.” We describe the estimation and regularization methods we use for GAS trees and forests in Section (ref).
We adapt the idea of regression trees for application to generalized autoregressive score (GAS) models. For a given tree structure with $J$ terminal nodes, $\boldsymbol{\mathcal{P}}=\{\mathcal{P}_1,..., \mathcal{P}_J\}$, the GAS(1,1) tree is based on the evolution equation:
where $\theta_j \equiv [\omega_j,\beta_j,\alpha_j]$ are the GAS parameters for partition $j$, and $s_t$ is as in equation ((ref)). By allowing the parameters of the GAS model to vary across partitions we greatly increase the flexibility of this class of models to fit the data. Furthermore, by retaining the GAS structure for each partition, we can more easily interpret how the tree structure improves the fit of the model, in contrast with more “black box” machine learning algorithms.
The parameters of the predictive density that are assumed constant in the baseline GAS model, denoted $\nu$ above, can either be held constant across partitions or can be allowed to vary.\footnote{In the kernel-based local M-estimation approach of oh2021better it is not possible to allow only a subset of parameters to vary with the state variable(s); the framework adopted in that paper requires an “all-or-nothing” assumption on parameter variation.} In our description below we impose they are fixed across partitions.
“Random forests” (breiman2001random) are an extension of regression trees designed to reduce the estimation error in predictions, while retaining the information contained in the tree-based forecast, see hastie2009elements for example. Similar to bootstrap aggregation, or “bagging,” a random forest is populated by trees that are each estimated on a bootstrap sample of the original data. In addition, each tree uses only a randomly-selected subset of the original state variables. The predictions from each of these trees are then averaged to obtain the random forest forecast.
If we denote trees in the random forest as $\boldsymbol{\mathcal{P}_b}$ for $b=1,...,B$, and the forecast from each tree for a given value of the vector of state variables, $\mathbf{Z_t}$ as $f_t^{(b)}(\mathbf{Z_t})$, obtained using equation ((ref)), then the GAS forest forecast is obtained simply as
We next turn to the estimation of the tree structure used in GAS trees and forests.
The estimation of a GAS tree requires finding the optimal state variables and thresholds from the set of candidate variables, as well as estimating the parameters of the GAS model. Finding the global optimum of this optimization problem is computationally infeasible in even moderately-sized regression tree applications, and to reduce the computational burden breiman1984classification proposed a greedy estimation algorithm that finds a near-optimal solution and converges quickly. The algorithm finds a state variable and a threshold to locally minimize the prediction error at each splitting step, continuing until a stopping criteria is satisfied.
Standard regression tree estimation involves estimating a regression separately for each terminal node in the tree, but given the autoregressive nature of GAS models, this is not possible for our application. We propose a modified estimation algorithm similar to audrino2001tree that uses a tree structure for the GAS model parameters and retains the autoregressive structure for $f_t$.
In the case that the number of terminal nodes, $J$, is one, there is no tree structure and the original GAS model is estimated via maximum likelihood:
For $J\geq 2$ we use the following estimation algorithm to estimate the tree structure or, equivalently, to find the optimal partition $\boldsymbol{\mathcal{P}}$. Estimation of the GAS tree involves Steps 1--5 below, and the GAS forest additionally uses Step 6.
Step 1: Denote the entire sample as the trivial partition $\boldsymbol{\mathcal{P}}^{(0)}$. Estimate the parameters of the model as in equation ((ref)), and denote these as $(\hat{\theta}^0,\hat{\nu}^0)$.
Step 2: Define a new partition: $\boldsymbol{\mathcal{P}}_{j,k}^{(m+1)}=\boldsymbol{\mathcal{P}}^{(m)}_{-j} \cup \{\mathcal{P}_{j,k,L}^{(m)},\mathcal{P}_{j,k,R}^{(m)}\}$ where $\boldsymbol{\mathcal{P}}^{(m)}_{-j}=\boldsymbol{\mathcal{P}}^{(m)} / \mathcal{P}_j$ contains all the partitions of $\boldsymbol{\mathcal{P}}^{(m)}$ except for the $j^{th}$, and the $j^{th}$ partition is split into “left” and “right” subpartitions based on the $k^{th}$ state variable and a threshold $c$
Step 3: Estimate the parameters for new subpartitions, taking the parameters of the other partitions, $\hat{\theta}^{(m)}_{-j}$, as fixed:\footnote{Optimizing the split is the most demanding step in the entire algorithm. This assumption significantly reduces the computational burden without hurting the results. Similar ideas are also implemented in the literature, for example athey2019generalized uses a gradient approximation in the split selection step.}
The compute the log-likelihood value at estimated parameter values.
Step 4: Maximize equation (ref) over the partition $j$, state variable $k$, and threshold $c$. Denote the optimized new partition as $\boldsymbol{\mathcal{P}}^{(m+1)}$ and estimate all the model parameters using equation ((ref)), denote these as $\hat{\theta}^{(m+1)}$.\footnote{In steps 3 and 4, parameter estimations are done by nonlinear solvers available in Scipy package of Python language. Following the warm start idea in the optimization literature, we use the estimates from the previous iteration, $\hat{\theta}^{(m)}$, as starting values for optimization procedure to accelerate the algorithm.}
Step 5: Repeat steps 2-4 until the depth of the tree, $m$, reaches a prespecified maximum value, $M$. The depth of the tree controls the model complexity, and we consider values of $M$ between one and six. We tune this parameter using a validation sample.
Step 6: For the GAS forest, repeat steps 2-5 for $B=200$ trees.\footnote{In preliminary analyses we obtained very similar results for $B=500$.} Each tree in the forest uses bootstrap data obtained from a circular block bootstrap (see, e.g., politis1999subsampling), with block length of 100 observations, and a random selection of one-third of the total state variables. One-third is a common choice in the machine learning literature, see hastie2009elements for example. The forecasts from each of the bootstrap trees are then averaged to obtain the GAS forest forecast.
All computations are done using the Duke Computing Cluster exploiting multiple computing nodes. We parallelize the split optimization steps, and use the Numba package to speed up the code. Estimating a single tree takes around five minutes with forty CPUs. We apply fixed-window estimation all models: we estimate the model parameters using the estimation sample and use those parameters to compute all out-of-sample forecasts.
We apply our new GAS tree and forest models in four out-of-sample forecasting analyses. We firstly consider forecasting S&P 500 return volatility using the GARCH model of bollerslev1986generalized, followed by predicting the entire conditional density of S&P 500 returns using the “t-GAS” model of creal2013generalized. In our third application we consider a flexible model for the joint distribution of S&P 500 returns and 10-year U.S. government bond returns, motivated in part by work on the switching sign of this correlation, see guidolin2006econometric and baele2010determinants, using the t-GAS copula to link t-GAS models for the marginal distributions, as in janus2014long. Finally, we consider the “autoregressive conditional duration” model of engle1998autoregressive, using high frequency transaction data on the exchange traded fund tracking the S&P 500 index, the SPY. These four applications represent a range of predictive environments, and we provide evidence of the merits of GAS trees and forests in each of them.
In addition to the baseline GAS model for each application, we consider two other benchmark models. The first is a low-dimensional GAS tree (“small GAS tree”), in which we only consider the lag of the dependent variable(s) as a state variable. This model is similar to that of audrino2001tree, and comparing the GAS tree and forest forecasts with this benchmark reveals the benefits of using a larger set of state variables. The third benchmark model is the “distributional random forest” of schlosser2019distributional. This model has no time series dynamics, but can provide a flexible distributional fit through the forest structure.
We compare all models in terms of one-step-ahead predictive performance. For the volatility forecasting application, we use the QLIKE loss function with realized volatility as the volatility proxy, see Patton2011 for details. For the remaining applications we use the negative log-likelihood, which is a consistent scoring rule for density forecasts, see Gneiting2007. We conduct dieboldMariano1995 tests of equal predictive accuracy, using newey1987simple standard errors based on 10 lags.
We consider GAS trees and forests in four empirical applications. The first three of these use daily data on S&P 500 index returns and 10-year U.S. government bond returns from January 2000 to December 2021, a total of 5447 observations. Our fourth application uses high frequency trade durations for the S&P 500 index tracker fund, SPY, during the calendar year 2021, and has 5,100 observations. In all applications we split the sample into three sub-periods: an estimation sample (first 30% of observations), a validation sample for optimizing hyperparameters (next 30%), and a test sample for out-of-sample forecast comparisons (remaining 40%).
We consider ten state variables for use in the applications based on daily data, and we add three high frequency state variables in the fourth application. We firstly include the (lagged) return on the S&P 500 index and the 10-year U.S. government bond, to capture any nonlinearities omitted by the GAS models. We next consider three measures of volatility: 5-min subsampled realized volatility (RVOL) on the S&P 500 index, a one-month (backward-looking) rolling average of RVOL, motivated by prominent HAR model of corsi2009simple, and the VIX index, a measure of S&P 500 volatility implied by options prices. We then consider three measures from the fixed income market: the federal funds rate, the difference between 10-year and 3-month bond yields, representing the level and slope of the yield curve, and the “default spread” defined as the difference between BAA and AAA rate corporate bond yields. For our ninth state variable we include the economic policy uncertainty index proposed by baker2016measuring, based on newspaper coverage. This index tracks important policy related events like the failure of Lehman Brothers or presidential elections. We take a rolling monthly average of policy uncertainty index to eliminate the noise in the data. Our tenth state variable is time, to capture potential structural breaks, see for example coulombe2020machine and goulet2020macroeconomy. In our fourth application, we additionally consider three high-frequency state variables: the first lag of duration, which can capture nonlinearities missed by the benchmark model, the return on SPY over the last trade event period, which can capture leverage-type effects, and the market liquidity of amihud2002illiquidity, which can gauge whether the ACD model parameters differ during periods of high versus low liquidity.
We use the augmented Dickey-Fuller test (DickeyFuller1979) for each state variable to test for the presence of a unit root. We fail to reject the null of a unit root for the federal funds rate and the difference between 10 year and 3 month yield, and we take the first difference of these two variables. To avoid look-ahead bias, we use a one-period lag of the state variables when forming the tree and forest forecasts.
All of our data comes from the FRED database at the Federal Reserve Bank of St. Louis, with the following exceptions: the realized volatility data comes from the Oxford Realized Library; the high frequency data comes from the New York Stock Exchange's TAQ database; the 10 year bond return series is from liu2021reconstructing; and the policy uncertainty series is from baker2016measuring.\footnote{The data for the latter two variables is available at \url{https://sites.google.com/view/jingcynthiawu/yield-data} and \url{www.policyuncertainty.com}.}
The GARCH model of bollerslev1986generalized is widely used for forecasting asset return volatility, and has been shown to be difficult to beat in a range of applications, see hansen2005forecast. Assuming a zero conditional mean, the model is:
creal2013generalized show that this model can be interpreted as a GAS model for the scale parameter of the Normal distribution. Given this equivalence, the GAS tree and forest models for this case can also be labeled GARCH tree and forest models. The “distributional random forest” (DRF) of schlosser2019distributional in this application sets $\beta =\alpha =0$ and allows the intercept, $\omega$ to vary with the forest structure, while the “small GAS tree” model of audrino2001tree uses a decision tree with only $y_{t-1}$ as a state variable. We compare forecasts from these three models with those from the new GAS tree and GAS forest models introduced in Sections (ref) and (ref) in Table (ref).
We observe that each variant of tree-based GARCH model (small GAS tree, GAS tree and GAS forest) significantly outperforms the benchmark GARCH model with $t$ statistics all less than $-2.5$. Moreover, the GARCH tree and forest models significantly beat the DRF specification. Thus, in this first application, we find that tree structured models improve the out-of-sample forecast accuracy over simple GARCH, a conventional econometrics model, and DRF, a machine learning tool. Table (ref) also shows that the “small GAS tree” outperformed by the GAS tree and forest models, with $t$-statistics below $-2.7$, revealing that external variables carry important information about future volatility.
Interestingly, and in contrast with both the econometrics and the machine learning literatures which generally find ensemble methods tend to outperform forecasts from individual models, we find that the GAS tree outperforms the GAS forest, with a $t$-statistic of nearly five. We interpret this result by noting that random forests have the potential to improve forecast accuracy through variance reduction at the cost of increasing bias, see for example hastie2009elements. In our case, the variance reduction attained by the GAS forest cannot compensate the associated increased bias, leading to less accurate forecasts.
To understand the source of forecast gains from the GAS tree model, Figure (ref) presents the estimated tree structure. The optimal tree depth was found to be three, with three different splitting variables. The algorithm first chooses the S&P 500 return with a threshold value $-0.16$ (its $40^{th}$ percentile) which approximately splits the sample using positive and negative market returns, consistent with an asymmetric reaction of future volatility to past returns, also known as a “leverage effect” (black1976). The second split in the tree is for positive returns and uses realized volatility with a threshold of 0.47 (10.9% in annualized standard deviation form), corresponding to the $45^{th}$ percentile of RVOL (conditional on the first split), thus approximately splitting positive return days into “high” and “low” volatility days. The third and final split is for low volatility days and uses VIX with a threshold of 14.01, corresponding to its $45^{th}$ conditional percentile. Recalling that the “variance risk premium” (carrWu2008; BollerslevVRP2009) can be approximated as the difference between VIX$^2$ and RVOL, the four terminal nodes of the tree in Figure (ref) can be interpreted, approximately, as those associated with (1) negative returns, (2) positive returns and high realized volatility, (3) positive returns, low realized volatility and low variance risk premium, (4) positive returns, low realized volatility and high variance risk premium.\footnote{Employing a semi-structural regime switching model, baele2010determinants finds that variance risk premium is an important economic factor in explaining stock return volatility.}
We next consider the problem of forecasting the entire distribution of daily returns on the S&P 500 index. Our baseline model is the t-GAS model introduced by creal2013generalized, which captures both excess kurtosis, through the use of the Student's $t$ distribution for the standardized residuals, and time-varying volatility, through the GAS structure for the scale parameter. Assuming a zero conditional mean, the t-GAS model is:
where $\nu$ is the degrees of freedom parameter for the $t$ distribution. As in creal2013generalized, $\nu$ is assumed constant, while $\sigma_t^2$ varies over time. The dynamics of $\sigma_t^2$ differs from the familiar GARCH structure when $\nu < \infty$, and simplifies to the GARCH model when $\nu \rightarrow \infty$. The $\{ \cdot \}$ term in equation ((ref)) implies a more moderate reaction to a large past return than in the GARCH model, as large returns are more common under the $t$ distribution than the Normal distribution.
Table (ref) presents comparisons of the t-GAS model with the distributional random forest (DRF), the “small GAS tree” (which only uses the lagged return as a state variable), and the GAS tree and forest models, both of which use all ten state variables described in Section (ref). The first column shows that the t-GAS model is significantly out-performed by all four competing models, with Diebold-Mariano $t$-statistics less than -3.5 in all cases. We also observe that the GAS tree significantly outperforms the “small GAS tree” and also the GAS forest, both with $p$-values less than 0.01. The GAS tree also outperforms the DRF forecast, but the difference is not statistically significant at the 5% level.
Figure (ref) shows the estimated t-GAS tree. Unlike the structure for the GARCH tree in the previous section, this tree only has depth of two, but those first two levels are identical to those of the GARCH tree.\footnote{We use a grid of $19$ values, corresponding to the 0.05, 0.10, ..., 0.95 quantiles of the state variable, for the threshold for each state variable, making our finding of identical threshold values less surprising.} The three terminal nodes of the tree have roughly equal numbers of observations, and can be interpreted as (1) negative returns, (2) positive returns and low realized volatility, (3) positive returns and high realized volatility.
We now focus on forecasting the joint distribution of stock and bond returns, using the S&P 500 index and 10-year Treasury bond for this purpose. We construct this model by combining t-GAS models for the marginal distributions, as utilized in the previous section, see equation (ref), with a Student's $t$ copula, as in janus2014long. Copulas are convenient tools for capturing the dependence between variables separately from the marginal distributions of each of the variables; see PattonHBOOK2013 for a review of these methods for economic time series. The use of a $t$ copula allows for the possibility of tail dependence, or “joint crashes and joint booms.” Although the marginal distributions and the copula are all from the Student's $t$ family, this joint distribution is not bivariate Student's $t$ distribution unless all three degrees-of-freedom parameters are identical; instead, this joint distribution allows for varying degrees of fat tails in each marginal distribution and the joint tails. The Student's $t$ copula with GAS dynamics for the correlation parameter is:
where $x_{i,t} \equiv F_{\text{Student}}^{-1}(u_{i,t};\nu_i)$ uses the inverse $t$ CDF with degrees-of-freedom parameter $\nu_i$, and $g$ and $w_t$ are scalars defined in Appendix (ref). As in the univariate $t$-GAS model, we impose that the copula degrees of freedom parameter, $\nu$ is constant over time.
Table (ref) shows the out-of-sample performance of competing models. We see that the benchmark GAS model significantly beats the distributional random forest (DRF), unlike in the two univariate applications, but it is beaten by both the “small tree” and the GAS tree models. The best-performing model is the GAS forest, which significantly outperforms the GAS model, with a $t$-statistic of -3.7. The GAS forest also significantly beats the DRF, but does not significantly beat either of the tree models. Interestingly, in this application the GAS tree reduces to the “small GAS tree” model, in that the only state variables selected for use in the tree structure are lags of the stock and bond returns.
In contrast with the univariate applications, the best-performing model is the GAS forest, not the GAS tree, and so we cannot present a tree diagram to better understand the structure of the best model. In its place, we consider two methods for interpreting the optimal model. Firstly, we conduct a leave-one-out analysis to measure the importance of each state variable. Specifically, we drop each state variable from the analysis, one at a time, and re-compute the optimal GAS forest forecasts. We then compare the average out-of-sample average loss from the original GAS forest and the GAS forest using one fewer state variable. If the difference is small, then the omitted state variable is unimportant, while if the difference is positive, then the omitted variable is important for forecast performance.\footnote{As this is an out-of-sample comparison of models, it is possible that the difference is negative, meaning that the smaller GAS forest is preferred to the original GAS forest. With a large enough sample size, including irrelevant state variables leads to no change, positive or negative, in out-of-sample loss, as such state variables will never be selected. In finite samples, however, irrelevant state variables may be mistakenly included.} We can use Diebold-Mariano tests to determine whether the change in out-of-sample loss is statistically significant. Figure (ref) presents the results of this analysis, and shows that the most important state variable for the GAS forest is the lagged bond return, T10Y, followed by the lagged stock market return, SPX. Omitting either of these significantly (at the $5\%$ level) deteriorates the GAS forest forecasts. The slope of the term structure (T10Y3M) and the volatility index (VIX) are also found to be important for the quality of GAS forest forecasts. Interestingly, we observe a statistically significant improvement in forecast performance by omitting the Federal funds rate (FFR) as a state variable, indicating that this variable is unhelpful for out-of-sample forecasting, but is selected for inclusion in the forest often enough to deteriorate the forecast.
We next analyze the impact of the most important state variable on the GAS forest model by plotting the parameters of the GAS model (recall equation (ref)) as a function of the state variable, see Figure (ref).\footnote{Similar plots for the other variables found to be significant in Figure (ref) are presented in the supplemental appendix.} The parameters $\beta$ and $\alpha$ are interpretable as the persistence and reaction-to-news of the model. The intercept, $\omega$ is not directly interpretable, and we instead plot $\omega/(1-\beta)$ which is interpretable as the long-run level of the GAS process. As the GAS forest involves averaging 200 bootstrap samples, each based on a random subset of state variables, we construct this plot by averaging the GAS parameters within each 1% quantile of the state variable.
Figure (ref) presents the results for the stock return as a state variable. The upper-left panel shows that the long-run correlation peaks when the stock market return is around zero, at about $0.35$. It declines to about $-0.45$ as the stock return increases to $2\%$, while it declines markedly to nearly $-0.6$ when the stock return is $-2\%$. This is interpretable as a “flight-to-quality” effect, with low stock market returns leading to more negative comovements between the stock and bond markets. Figure (ref) in the supplemental appendix plots the corresponding results with the bond return as the state variable, and that figure is also consistent with a flight-to-quality effect.
The upper-right panel of Figure (ref) shows that the persistence of the GAS model is roughly unrelated to the stock market return. The lower-left panel shows that the GAS model reacts about $20\%$ more strongly to news when the stock market is down versus up: the $\alpha$ parameter is $0.046$ when stocks are down, while it is $0.038$ when stocks are up. This is consistent with investors paying closer attention to bad news than good news, a finding similar to that of PattonSheppard2015 in a different context.
The lower-right panel of Figure (ref) shows the predicted correlation from the GAS forest as a function of the stock market return, and reveals an inverted U-shaped pattern, though with substantial noise. Without an underlying model to guide interpretation, one might be hesitant to draw too much from this panel. With the benefit of the GAS structure underlying our forest forecast, we know that this shape is primarily coming from the long-run level, $\omega/(1-\beta)$, in the upper-left panel, and that that relationship is strong. This reveals an important benefit of combining machine learning tools with economically motivated, and/or empirically successful, econometric models.
Finally, we consider the problem of forecasting the time between consecutive trade events, known as a “trade duration.”\footnote{A “trade event” could be a single transaction occurring, or a total of $x$ transactions occuring, or a total of $\$y$ value of transactions occuring, or a total of $z$ shares being transacted, or some other event defined as a function of characteristics of transactions.} Trade durations are a measure of market activity, and are important for high-frequency risk management and transaction cost minimization. We take as our benchmark the “autoregressive conditional duration” (ACD) model of engle1998autoregressive. Denoting $y_t$ as the time (in minutes) between consecutive trade events, the ACD model assumes an exponential distribution for $y_t$ with a time-varying conditional mean, $\mu_{t}$:
creal2013generalized show that the ACD model is also a special case of a GAS model, allowing us to consider it in our study of tree- and forest-based extensions of GAS models. See BauwensHautsch09 for a review of ACD and related models.
For our empirical analysis in this section, we use high-frequency data on SPY, an exchange traded fund tracking the S&P 500 index, between January 1$^{st}$ 2021 and December 31$^{st}$ 2021. We study the time taken for 10,000 shares of SPY to be transacted, leading to 5,100 durations during this sample period, and corresponding to an average duration of 14.9 minutes.
Table (ref) presents the out-of-sample forecast performance of the baseline ACD model as well as the competing models considered in previous sections: the distributional random forest (DRF), the ACD tree using only lagged durations as a state variable (Small ACD tree), the ACD using all 13 state variables, and the ACD forest model. We firstly observe that the DRF model for durations, which is pure machine learning tool, is beaten by all other models including the benchmark with $t$-statistics all less than $-2.4$.
The baseline ACD model has higher loss in the out-of-sample compared with to tree- and forest-based extensions. Interestingly, the “small ACD tree” model, which only uses lagged duration as a state variable, significantly beats the baseline ACD model, while the “ACD tree” model, which considers 13 state variables including lagged duration, does not significantly beat the baseline model. This reveals the value in imposing some structure (namely, reducing the number of potential state variables) on the tree-based extension in this application. The “ACD forest” model significantly beats the baseline ACD model, revealing the forecast gains available from averaging forecasts from randomly formed trees, consistent with breiman2001random in a linear regression setting.
To understand how forecast accuracy is improved in the ACD forest model, we calculate the variable importance measure for each of the state variables, introduced in the previous section, and present the results in Figure (ref). The horizontal bars show the increase in average loss function from omitting a state variable, and a positive value indicates that that state variable is important for forecasting. The lines refer to 95% confidence intervals computed from Diebold-Mariano tests. We find that the volatility variables RVOL and VIX are two most important state variables, despite the fact that these are measured only daily, and so are constant within a trade day. The next two most important state variables are both high-frequency variables: amihud2002illiquidity liquidity, and duration. Interestingly, we find that omitting default spread and time from the set of potential state variables actually improves forecast accuracy, indicating these are harmful when used in a forest-based ACD model.\footnote{Figure (ref) in the supplemental appendix presents the optimal tree structure for the ACD tree. We find three terminal nodes in this structure, all reflecting the direction of the stock and bond markets. We find one state when the stock market is up (representing 55% of the sample), another when the stock market is down and the bond market is not strongly up (representing 92% of the remaining sample), and a small third state when the stock market is down and the bond market is strongly up (representing just 3% of the total sample).}
In Figure (ref) we plot the parameters of the forest ACD model as a function of the most important state variable, realized volatility (RVOL). We see that the long-run average duration implied by the model is highest when volatility is low, at around 23 minutes, and it steeply declines as volatility increases to about 10%, at around 19 minutes. Persistence, measured by $\alpha+\beta$ in the ACD model, is lowest when volatility is low, and it increases sharply with volatility to around 10% and is approximately flat beyond that. The parameter governing the reaction of the model to news, $\alpha$, is essentially flat as a function of volatility. The pattern for the forecasts from the ACD forest model, in the lower-right panel, shows that predicted durations are around 12 minutes when volatility is above 10%, while they are around double that when volatility is low. This is consistent with the positive volume-volatility relationship (see, e.g., Karpoff87): volume and durations are negatively correlated (longer durations correspond to lower volumes, and vice versa) and so in periods of lower volatility average trade durations tend to be longer.
Since its publication a decade ago, the class of generalized autoregressive score (GAS) models of creal2013generalized and harvey2013dynamic has proven to be a popular, parsimonious way to capture time variation in the parameter(s) of a given model. Its parsimonious nature, however, means that some important exogenous information or nonlinearities may be neglected, and how to best incorporate such additional features is difficult to determine ex ante. We propose adapting methods from machine learning to search across a wide range of exogeneous variables and capture various forms of nonlinearity: the “GAS tree” combines the parsimonious structure of the GAS model with the flexibility of decision trees (breiman1984classification,breiman2017classification), and the “GAS forest,” analogous to the random forests of breiman2001random, averages the forecasts from many GAS trees each produced on a bootstrap sample of the original data. Our GAS tree and GAS forest models can be applied whenever a GAS model is considered, and require from the researcher only a set of exogenous variables that are thought to be possibly useful.
We apply the proposed GAS tree and GAS forest models in four diverse applications: forecasting stock return volatility, the distribution of stock returns, the joint distribution of stock and bond returns, and high-frequency trade durations. We find that the proposed extensions lead to significantly improved forecasts in all four applications. We moreover uncover economic explanations for the sources of these forecast gains. Through inspections of the optimal GAS tree structures, and variable importance and parameter sensitivity analyses for GAS forest forecasts, we find that the best-performing GAS tree and forest models are those that incorporate well-known empirical regularities, such as the leverage effect in volatility, the flight-to-quality effect in stock-bond correlations, and the volume-volatility relationship in trade durations.
Faced with evolving data generating processes and the resulting “small” data sets, the success of machine learning methods in economics and finance relies on good, parsimonious benchmark models as reference points, an observation made nicely in IsraelKellyMoskowitz2020. We used GAS models for this purpose; future work may consider augmenting a different class of forecasting models with machine learning methods.