Capturing Heterogeneity: Machine Learning Approaches to Implied Volatility Forecasting*

Abstract:

1 Introduction

Accurate forecasting of the implied volatility (IV) surface is crucial for practitioners, as it directly informs option pricing, risk management decisions, and trading strategies. Recent work by Dufays et al., 2025 proposes the Surface Heterogeneous Autoregressive (SHAR) model, which effectively and efficiently captures volatility surface dynamics by exploiting persistence over daily, weekly, and monthly horizons. However, SHAR and other existing models typically apply one set of forecasting model parameters uniformly across all strikes and maturities. This homogeneous specification may be restrictive given that option markets exhibit substantial heterogeneity: short-dated, long-dated, out-of-the-money (OTM), and near-the-money (NTM) options behave very differently in levels and in how they respond to news Pan, 2002; Christoffersen et al., 2009. Indeed, due to violations of Black and Scholes, 1973 assumptions, implied volatilities vary systematically across strikes and maturities, forming the so-called implied volatility surface. In the time series, S&P 500 OTM put implied volatilities tend to be much higher and more volatile than at-the-money (ATM) or OTM call implied volatilities. Short-dated option implied volatilities are, on average, significantly more volatile and less persistent than long-dated option implied volatilities. When dynamics systematically vary across moneyness and maturity dimensions, imposing homogeneous model parameters across the entire surface may limit the potential for performance gains.

In this paper, we propose a machine learning approach to optimally partition the cross section of S&P 500 index options in implied volatility forecasting by extending the SHAR model. Specifically, we identify regions of the moneyness-maturity space where the SHAR model requires different parameter sets, acknowledging that distinct option segments exhibit distinct SHAR dynamics. We allow SHAR parameters to vary across the moneyness-maturity space, using regression trees to identify data-driven partitions where distinct forecasting models perform best. This approach substantially improves forecasting accuracy by adapting the model to regional heterogeneity across the surface and leveraging these localized dynamics.

Market segmentation and variation in investor clienteles across the option surface are well-documented in the literature. Studies have shown that option markets are segmented across strikes and maturities, with different investor clienteles focusing on specific regions of the option surface Bollen and Whaley, 2004; G\^arleanu et al., 2009; Bakshi et al., 2010. These clienteles—including volatility sellers focusing on short-dated options, directional speculators trading out-of-the-money options, and institutional hedgers using near-the-money contracts—create distinct supply and demand patterns across the surface. Buraschi and Jiltsov, 2006 show that considering different beliefs/information from heterogeneous agents improves explaining option market dynamics. Such documented market segmentation provides compelling motivation for our approach: allowing for localized dynamics substantially improves the performance of implied volatility forecasting across the surface.

Our empirical approach employs the partitioning logic of regression trees of Breiman et al., 1984; Breiman et al., 2017 and lets the data decide where the option surface should be split, using option characteristics as splitting variables. Concretely, we use moneyness and maturity as splitting variables and build a tree that recursively partitions the option panel into subregions that are more homogeneous in terms of forecasting performance. At each terminal node, we thereby obtain a separate model estimated using only the options that fall into that region. This tree-based local estimation scheme preserves the structure and economic interpretation of the underlying models, while allowing their parameters to adapt flexibly across the surface. We particularly assess whether the tree structure can improve out-of-sample forecasts of the implied volatility surface.

We apply this tree-based partitioning approach to the SHAR framework for implied volatility surfaces, and we call this approach a tree-augmented Surface HAR model (Tree-SHAR).¹ Dufays et al., 2025 propose a two-step approach in which each daily surface is first fitted by a chosen option model or nonparametric method, and then the SHAR model linearly combines past fitted surfaces to forecast future implied volatilities on the unbalanced option panel, exploiting strong persistence over daily, weekly, and monthly horizons. We extend their setup by allowing the SHAR coefficients to vary across the tree’s moneyness–maturity nodes. In the resulting Tree-SHAR model, each node has its own set of HAR-type parameters, so that short-dated, long-dated, OTM call, or OTM put options can have different implied volatility dynamics among themselves. We repeat the analysis by further adapting advanced machine learning approaches, a bootstrap-aggregated ensemble of tree models following Breiman, 1996 and a boosted tree model following Friedman, 2001.

Our empirical analysis using a large cross section of daily S&P 500 index option data shows that tree-based SHAR models deliver statistically significant improvements in out-of-sample forecast accuracy relative to the Random Walk and the baseline global SHAR model, while the Boosted Tree-SHAR model achieves the lowest forecast errors overall. The improvements are particularly pronounced in regions with greater cross-sectional heterogeneity: OTM calls, deep OTM puts, and short-maturity options. The estimated tree structures reveal systematic partitioning of the implied volatility surface by moneyness and maturity, capturing distinct dynamics across different option segments. For example, region-specific HAR coefficient estimates show that short-term option implied volatilities exhibit very high daily persistence, while long-term option implied volatilities display relatively higher weekly and monthly persistence.

The tree structure remains relatively stable under similar market conditions but adjusts substantially following market shocks. For example, (deep) OTM call IVs depend more on the weekly moving averages and are negatively associated with the monthly moving averages during bad times, whereas they depend more on the previous-day information during good times. These different dynamics of implied volatilities may be due to different market expectations toward upside potentials during good and bad periods. For long-term OTM puts, implied volatilities show sizable dependence on both past weekly and monthly information during good times, but the dependence on monthly information disappears during bad times. Moreover, the tree algorithm tends to divide the volatility surface into more heterogeneous dynamics in OTM calls and NTM options but more homogeneous dynamics in a wide range of short-term OTM puts during bad times. These types of variations are not captured by the global SHAR specification.

Our Tree-SHAR results convey the following messages: appropriately incorporating cross-sectional option characteristics—here, via tree partitions in moneyness and maturity—helps identify the segmentation of option markets and their heterogeneous features and systematically improves time-series implied-volatility forecasting. Our approach is deliberately conservative: it builds on a simple but well-understood benchmark model, preserves their economic structure, and uses the tree only as a data-driven device to choose where separate local models should be estimated. At the same time, the method is flexible and modular; it can be combined with alternative surface fitters, forecasting models, or other splitting variables, and it can be implemented in real time using only modest computational resources.
A large body of empirical literature documents heterogeneity in option markets. Dumas et al., 1998, Constantinides et al., 2009, and Fan and Mancini, 2009 show that standard option pricing models produce unstable parameters over time and systematic pricing errors that vary with moneyness and maturity, underscoring the limitations of a single functional form. Additional examples of inferring heterogeneous implied model parameters include Bakshi et al., 1997, Hurn et al., 2015, and Dufays et al., 2023. Our paper differs from these studies in that we focus on forecasting the time series evolution of the volatility surface rather than fitting option prices.

Our paper is most closely related to the literature on forecasting implied volatility surfaces, particularly Dufays et al., 2025, who propose the SHAR forecasting model. Recent works by Almeida et al., 2023 and Chen et al., 2025 apply neural networks to forecast implied volatility surfaces. The key distinction of our approach is that we allow forecasting model parameters to vary across the moneyness-maturity dimension through data-driven surface partitions. This not only improves forecasting accuracy compared to benchmark models but also reveals economically meaningful surface segmentation. In contrast to purely machine learning methods, our approach leverages machine learning to capture local heterogeneity while preserving the structural tractability of standard models.

Methodologically, our framework aligns with the literature on “tree-structured” parametric models. Audrino and B\"uhlmann, 2001 pioneered the use of regression trees to estimate local GARCH parameters, and Patton and Simsek, 2026 extend this approach to GAS models. Medeiros et al., 2021 and Goulet Coulombe, 2024 demonstrate that tree-based methods successfully identify regime-dependent dynamics in macroeconomic forecasting. While these studies primarily exploit tree-based partitioning in the time-series dimension, Brownlees and Souza, 2025 and Oh and Patton, 2026 apply similar methods to the cross-section of stocks and mutual funds, respectively. We extend this logic to the cross section of options, using trees to partition the moneyness-maturity space and estimate local SHAR parameters. To the best of our knowledge, this is the first application of tree-structured parameter heterogeneity to option-implied volatilities. Our work also relates to the broader “local estimation” literature, which allows model parameters to vary with state variables. This literature has developed along two paths: kernel-based localization for (quasi-) maximum likelihood estimation Tibshirani and Hastie, 1987; Fan et al., 1998; Fan et al., 2009; Oh and Patton, 2024; Kim and Oh, 2025, and tree-based localization for regression, parametric distributions, and GMM models Breiman, 2001; Athey et al., 2019; Schlosser et al., 2019.

The remainder of the paper is organized as follows. Section 2 reviews the SHAR framework for dynamic implied volatility surfaces and outlines the Tree-SHAR extension, together with its bagged and boosted enhancements, emphasizing improvements in time-series forecasts. Section 3 describes the option data and filtering rules. Section 4 presents the empirical results for forecast performance, including out-of-sample forecasts of implied volatility for benchmarks, Tree-SHAR and its extensions. Section 5 concludes.

2 Tree-Based Surface HAR Models

This section develops tree-based extensions of the Surface HAR (SHAR) model for implied-volatility (IV) surface forecasting. Section 2.1 introduces the baseline SHAR model, which extends the HAR structure of Corsi, 2009 to the option panel by modeling IV as a function of lagged fitted-surface components. This specification provides an interpretable decomposition of IV dynamics into short-, medium-, and long-run components. Section 2.2 provides empirical evidence of divergent IV dynamics across the moneyness and maturity dimensions. This empirical observation motivates considering heterogeneous dynamics over the volatility surface for forecasting. Section 2.3 generalizes the baseline model by allowing the HAR dynamics to vary across the moneyness–maturity space. Specifically, we partition the implied-volatility surface using a regression tree and estimate leaf-specific HAR coefficients, thereby capturing cross-sectional heterogeneity while preserving the HAR interpretation within each region. Section 2.4 introduces a bagged Tree-SHAR model that improves stability by averaging forecasts across bootstrap trees. Finally, Section 2.5 develops a boosted Tree-SHAR model that constructs an additive expansion of local HAR corrections, sequentially reducing remaining forecast errors and capturing finer nonlinearities and interactions over the option surface.

2.1 Baseline Surface HAR model

Our benchmark time-series model for the implied volatility surface is the Surface Heterogeneous Autoregressive (SHAR) framework of Dufays et al., 2025. Their forecasting setup is explicitly two-step.

Step 1: daily cross-sectional surface fit. From an initial day onward, we observe a daily stream of implied volatility (IV) surfaces. On each day \(\ell\), we observe \(N_{\ell}\) option-implied volatilities with option \(i\)’s characteristics \(\mathbf{z}_{i,\ell}=(m_{i,\ell},\tau_{i,\ell})'\), where \(m_{i,\ell}\) denotes moneyness and \(\tau_{i,\ell}\) denotes time-to-maturity. Let \(\mathrm{IV}(\mathbf{z}_{i,\ell})\) denote the Black–Scholes implied volatility of option \(i\) on day \(\ell\).

In the first step, each day’s IV surface is fitted cross-sectionally using an arbitrary model \(M\) for implied volatilities: \[\begin{equation} \mathrm{IV}(\mathbf{z}_{i,\ell}) = \mathrm{IV}^{M}\bigl(\mathbf{z}_{i,\ell};\Theta_{\ell}\bigr) + \varepsilon_{i,\ell}, \qquad i = 1,\dots,N_{\ell}, \quad \ell = 1,\dots,T,\tag{1} \end{equation}\] where \(\mathrm{IV}^{M}\bigl(\cdot;\Theta_{\ell}\bigr)\) is the model-implied volatility surface generated by parameter vector \(\Theta_{\ell}\) (estimated on day \(\ell\)), and \(\varepsilon_{i,\ell}\) is an error. The model \(M\) can be a structural option pricing model (for example, Heston, 1993’s (Heston, 1993) stochastic volatility model), a parametric IV surface, or a flexible nonparametric or machine-learning method. Equation (1) therefore delivers, for each day \(\ell\), a fitted surface \(\mathrm{IV}^{M}\bigl(\mathbf{z}_{i,\ell};\widehat{\Theta}_{\ell}\bigr)\) evaluated on the unbalanced set of observed \(\mathbf{z}_{i,\ell}=(m_{i,\ell},\tau_{i,\ell})'\) points. In this paper, we employ the ad-hoc Black-Scholes (AHBS) specification, which provides an easy-to-implement but sufficiently flexible approach used by Dumas et al., 1998 and Dufays et al., 2025, among others. The AHBS model approximates the IV surface using a second-degree polynomial in moneyness and maturity: \[\begin{equation*} \mathrm{IV} (\mathbf{z}_{i,\ell}) =\beta_0^{\text{M}} + \beta_1^{\text{M}} m_{i,\ell} + \beta_2^{\text{M}} m_{i,\ell}^2 + \beta_3^{\text{M}} \tau_{i,\ell} + \beta_4^{\text{M}} \tau_{i,\ell}^2 + \beta_5^{\text{M}} m_{i,\ell}\tau_{i,\ell} + \varepsilon_{i,\ell}. \end{equation*}\]

A widely used benchmark forecast treats the fitted surface as a random walk in shape: on day \(t\), estimate \(\widehat{\Theta}_{t}\) using only the day-\(t\) cross-section, and then forecast the IV surface at day \(t+h\) by evaluating the same fitted model at the future option characteristics: \[\begin{equation} \widehat{\mathrm{IV}}^{\mathrm{RW}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) = \mathrm{IV}^{M}\bigl(\mathbf{z}_{i,t+h};\widehat{\Theta}_{t}\bigr), \qquad i=1,\dots,N_{t+h}.\tag{2} \end{equation}\] That is, \(\mathrm{IV}^{M}(\mathbf{z}_{i,t+h};\widehat{\Theta}_{t})\) denotes the model-implied volatility obtained by evaluating the day-\(t\) fitted volatility surface, parameterized by \(\widehat{\Theta}_{t}\), at the target moneyness and time-to-maturity point \(\mathbf{z}_{i,t+h}=(m_{i,t+h},\tau_{i,t+h})'\) for date \(t+h\). The random-walk benchmark therefore assumes that the shape of the fitted volatility surface remains unchanged between dates \(t\) and \(t+h\). This “RW” model is computationally convenient because it only requires re-estimating the cross-sectional surface model on a single day’s option panel.

As an alternative approach, we also consider generating the volatility surface using the parsimonious Stochastic Volatility Inspired (SVI) method of Gatheral, 2004. This method uses a different functional form of IV than the AHBS method to span the implied volatility curve given maturity. We specifically follow Beason and Schreindorfer, 2022 to extend it to the volatility surface, \(i.e.\), as a function of both moneyness and maturity. We discuss more details of the method in Appendix A. Overall, we obtain consistent forecasting performance regardless of whether we generate the surface following the AHBS or the SVI method.

The surface-fitting procedure in Step 1 is common to all forecasting specifications considered in this paper. The baseline SHAR model, Tree-SHAR, Bagged Tree-SHAR, and Boosted Tree-SHAR differ only in Step 2, which maps the history of fitted IV surfaces into forecasts of future implied volatility. Thus, the improvements from the tree-based methods come from modifying the second-step forecasting rule—in particular, by allowing the SHAR coefficients to vary across the moneyness–maturity surface—while holding the first-step surface construction fixed. We now describe this second step.

Step 2: dynamic modeling of fitted surfaces. The SHAR approach replaces the RW benchmark (2) with an autoregressive forecasting rule that exploits persistence in the sequence of fitted IV surfaces. In our empirical implementation, we follow Dufays et al., 2025 and adopt a heterogeneous autoregressive (HAR) specification in the spirit of Corsi, 2009. For a given forecast horizon \(h\), define the HAR regressor vector \[\begin{equation} \mathbf{x}_{h,i,\ell} = \Biggl( 1,\ \mathrm{IV}^{M}\bigl(\mathbf{z}_{i,\ell};\widehat{\Theta}_{\ell-h}\bigr),\ \frac{1}{5}\sum_{r=0}^{4}\mathrm{IV}^{M}\bigl(\mathbf{z}_{i,\ell};\widehat{\Theta}_{\ell-h-r}\bigr),\ \frac{1}{22}\sum_{r=0}^{21}\mathrm{IV}^{M}\bigl(\mathbf{z}_{i,\ell};\widehat{\Theta}_{\ell-h-r}\bigr) \Biggr)'.\tag{3} \end{equation}\] The first element is an intercept, while the remaining three elements are the fitted IV surface from the previous day, the average of fitted IV surfaces over the previous week, and the average of fitted IV surfaces over the previous month, all evaluated at the same option characteristics \(\mathbf{z}_{i,\ell}\).

Using this notation, the baseline SHAR model is \[\begin{equation} \mathrm{IV}(\mathbf{z}_{i,\ell}) = \mathbf{x}_{h,i,\ell}'\boldsymbol{\beta}_{h} + \varepsilon_{h,i,\ell}, \qquad i=1,\ldots,N_{\ell},\quad \ell=h+22,\ldots,T,\tag{4} \end{equation}\] where \(\boldsymbol{\beta}_{h}=(\beta^{(0)}_{h},\beta^{(d)}_{h},\beta^{(w)}_{h},\beta^{(m)}_{h})'\) depends on the forecast horizon \(h\) but is constant across options and days.

To estimate the SHAR coefficients at horizon \(h\) using information available at date \(t\), we pool observations across days \(\ell=h+22,\ldots,t\) and across all options within each day. The pooled ordinary least squares estimator is \[\begin{equation} \widehat{\boldsymbol{\beta}}_{h|t} = \arg\min_{\boldsymbol{\beta}\in\mathbf{R}^{4}} \sum_{\ell=h+22}^{t} \sum_{i=1}^{N_{\ell}} \left( \mathrm{IV}(\mathbf{z}_{i,\ell}) - \mathbf{x}_{h,i,\ell}'\boldsymbol{\beta} \right)^2.\tag{5} \end{equation}\]

Given \(\widehat{\boldsymbol{\beta}}_{h|t}\) and the sequence of fitted surfaces up to date \(t\), the \(h\)-step-ahead SHAR forecast for an option with characteristics \(\mathbf{z}_{i,t+h}\) is²

\[\begin{equation} \widehat{\mathrm{IV}}^{\mathrm{SHAR}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) = \mathbf{x}_{h,i,t+h}'\widehat{\boldsymbol{\beta}}_{h|t},\tag{6} \end{equation}\] where \[\mathbf{x}_{h,i,t+h} = \Biggl( 1,\ \mathrm{IV}^{M}\bigl(\mathbf{z}_{i,t+h};\widehat{\Theta}_{t}\bigr),\ \frac{1}{5}\sum_{r=0}^{4}\mathrm{IV}^{M}\bigl(\mathbf{z}_{i,t+h};\widehat{\Theta}_{t-r}\bigr),\ \frac{1}{22}\sum_{r=0}^{21}\mathrm{IV}^{M}\bigl(\mathbf{z}_{i,t+h};\widehat{\Theta}_{t-r}\bigr) \Biggr)'.\] Thus, the SHAR forecast is a linear combination of the daily, weekly, and monthly fitted-surface components evaluated at the future option characteristics \(\mathbf{z}_{i,t+h}\).

For later use in the tree-based extensions, it is convenient to collect the observations used in the second-step SHAR estimation into a single notation. Fix the forecast horizon \(h\). At forecast origin \(t\), define the available SHAR sample as \[\begin{equation} \mathcal{D}_{t} = \left\{ \left( \mathrm{IV}(\mathbf{z}_{i,\ell}), \mathbf{z}_{i,\ell}, \mathbf{x}_{h,i,\ell} \right): i=1,\ldots,N_{\ell},\; \ell=h+22,\ldots,t \right\}.\tag{7} \end{equation}\] Each element of \(\mathcal{D}_{t}\) contains the implied volatility \(\mathrm{IV}(\mathbf{z}_{i,\ell})\), the option characteristics \(\mathbf{z}_{i,\ell}=(m_{i,\ell},\tau_{i,\ell})'\), and the SHAR regressor vector \(\mathbf{x}_{h,i,\ell}\) defined in (3). This notation will be used below when constructing Tree-SHAR, Bagged Tree-SHAR, and Boosted Tree-SHAR.

In our empirical analysis, this SHAR specification serves as the baseline time-series model for the implied volatility surface, alongside the RW benchmark in (2). The tree-augmented version introduced in Section 2.3 retains the SHAR structure in (4)–(6) but allows the coefficients to differ across regions of the moneyness–maturity space.

2.2 Heterogeneous characteristics in data

Before we describe details of the tree-augmented SHAR model, we briefly discuss heterogeneous characteristics of the implied volatility surface observed in data.

Figure 1 shows average implied volatilities by moneyness or maturity groups over time. We use three groups where the moneyness or maturity is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile, on each date. Panel A shows the IV by moneyness, and Panel B shows the IV by maturity. In Panel A, the blue solid line intrinsically represents the average IV for OTM calls, and the black dotted line represents the average IV for deep OTM puts. We observe that the average IV increases and is more volatile as moneyness increases (or moves toward deeper OTM puts). In Panel B, the level of IVs across maturity dimension remains similar, but the volatility of IVs differs significantly. Shorter-maturity IVs (blue solid line) show a much more volatile pattern than relatively longer-maturity IVs (red dashed and black dotted lines).

Figure 2 depicts autocorrelation functions (ACF) of implied volatilities and their first differences up to 22 lags. Panels A and B plot the ACF of IVs by moneyness and maturity groups, where the groups are identically defined as in Figure 1. Although the ACFs of IVs look similar across different moneyness levels, they increase as maturity increases. This pattern is more pronounced as the number of lags increases. On the other hand, daily changes in IVs in Panels C and D show a different pattern. Although their ACFs appear similar across different maturity levels, the ACFs vary across different moneyness levels. Specifically, OTM call IVs (blue bars) show relatively more negative first-order and positive second-order autocorrelation coefficients, respectively, implying that they are more likely to depend on lagged IVs.³

All of the above observations clearly indicate that the time series dynamics of IVs are heterogeneous across different moneyness and maturity. We now discuss how we incorporate such heterogeneity of the volatility surface into the forecasting framework through a regression tree.

2.3 Tree-augmented SHAR model

We next extend the benchmark SHAR specification in Section 2.1 by allowing the HAR dynamics to vary across regions of the moneyness–maturity space. The key idea is as follows: we use a regression tree \(\mathcal{T}\) to partition the option observations into terminal nodes (leaves), and we estimate a separate SHAR model within each leaf. This yields a localized forecasting rule that can differ across economically meaningful regions of the option surface, such as short-dated, long-dated, OTM, and near-the-money options.

Let \(m_{i,\ell}\) and \(\tau_{i,\ell}\) denote the moneyness (underlying price divided by strike price) and time to maturity of option \(i\) on day \(\ell\). We collect these option characteristics in⁴ \[\mathbf{z}_{i,\ell} \equiv (z^{(1)}_{i,\ell},z^{(2)}_{i,\ell})' \;=\; (m_{i,\ell},\tau_{i,\ell})' \in \mathcal{Z}.\] A tree \(\mathcal{T}\) on \(\mathcal{Z}\) induces a recursive partition \[\mathcal{Z} \;=\; R_1 \cup \cdots \cup R_{J(\mathcal{T})},\] where terminal nodes (leaves) \(R_j\) are mutually exclusive, and \(J(\mathcal{T})\) denotes the number of leaves. Each region \(R_j\) corresponds to a terminal node (leaf) and is characterized by a sequence of binary splits of the form \(m \le c\) or \(\tau \le c\) applied to the coordinates of \(\mathbf{z}_{i,\ell}\).

Using the HAR regressor vector \(\mathbf{x}_{h,i,\ell}\) defined in (3), the Tree-SHAR model replaces the global coefficient vector in (4) with a leaf-specific coefficient vector: \[\begin{equation} \mathrm{IV}(\mathbf{z}_{i,\ell}) \;=\; \mathbf{x}_{h,i,\ell}'\,\boldsymbol{\beta}_{h,j(i,\ell)} \;+\; \varepsilon_{h,i,\ell},\tag{8} \end{equation}\] where \(\boldsymbol{\beta}_{h,j}\) collects the intercept and HAR coefficients for leaf \(j\), and \(j(i,\ell)\) denotes the index of the terminal node (leaf) to which option \((i,\ell)\) is assigned under the tree partition, i.e., the unique \(j\) such that \(\mathbf{z}_{i,\ell}\in R_{j(i,\ell)}\).

For a fixed tree \(\mathcal{T}\), the leaf-specific coefficients are estimated by ordinary least squares within each leaf. That is, for each leaf \(j\), \[\begin{equation} \widehat{\boldsymbol{\beta}}_{h,j}(\mathcal{T}) \;=\; \arg\min_{\boldsymbol{\beta}\in\mathbf{R}^{4}} \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R_j} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\boldsymbol{\beta}\bigr)^2,\tag{9} \end{equation}\] and the Tree-SHAR objective is the sum of these squared residuals across leaves.

Directly solving equation (9) over all possible trees is infeasible, so we construct the tree greedily in the spirit of Classification and Regression Trees (CART; Breiman et al., 1984; Breiman et al., 2017). We start from the stump tree with a single leaf, denoted by \(R^{(0)}_1\), which contains all option-day observations in the sample \(\mathcal{D}_{t}\) defined in (7). The coefficient vector associated with this initial leaf coincides with the SHAR coefficient vector from Section 2.1. At a generic splitting step \(s\geq 1\), suppose we have a tree \(\mathcal{T}^{(s-1)}\) with \(J^{(s-1)}\) leaves and associated parameter vectors \(\{\widehat{\boldsymbol{\beta}}_{h,j}^{(s-1)}\}_{j=1}^{J^{(s-1)}}\).⁵ For each current leaf \(R_j^{(s-1)}\) and each splitting variable \(k\in\{1,2\}\) (moneyness or maturity), we consider splitting \(R_j^{(s-1)}\) into left and right child nodes at a candidate threshold \(c\) taken from a grid of \(z^{(k)}_{i,\ell}\) within that leaf.⁶ For each candidate split \((j,k,c)\) we:

1. Partition the observations in \(R_j^{(s-1)}\) into left and right subsets according to \(z^{(k)}_{i,\ell} \le c\) or \(> c\), discarding the split if either child contains fewer than a pre-specified minimum number of options: \[\begin{align*} R^{(s)}_{j,k,c,L} &= \{\, \mathbf{z}_{i,\ell} \in R^{(s-1)}_j : z^{(k)}_{i,\ell} \le c \,\}, \\ R^{(s)}_{j,k,c,R} &= \{\, \mathbf{z}_{i,\ell} \in R^{(s-1)}_j : z^{(k)}_{i,\ell} > c \,\}. \end{align*}\]

2. Re-estimate the SHAR parameters separately in the left and right child nodes by minimizing the leaf-specific sum of squared errors (SSE) within each child node, while leaving the parameter estimates in all other leaves unchanged at their previous values \(\widehat{\boldsymbol{\beta}}_{h,j'}^{(s-1)}\): \[\begin{align*} \widehat{\boldsymbol{\beta}}^{(s)}_{h,j,k,c,L} &= \arg\min_{\boldsymbol{\beta}\in\mathbf{R}^{4}} \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(s)}_{j,k,c,L}} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\boldsymbol{\beta}\bigr)^2, \\ \widehat{\boldsymbol{\beta}}^{(s)}_{h,j,k,c,R} &= \arg\min_{\boldsymbol{\beta}\in\mathbf{R}^{4}} \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(s)}_{j,k,c,R}} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\boldsymbol{\beta}\bigr)^2. \end{align*}\]

3. Compute the resulting total loss \(\mathcal{L}(\mathcal{T}^{(s)}_{j,k,c})\) implied by the candidate tree \(\mathcal{T}^{(s)}_{j,k,c}\), obtained by replacing \(R_j^{(s-1)}\) with its two children and updating the corresponding parameter vectors: \[\begin{align*} \mathcal{L}(\mathcal{T}^{(s)}_{j,k,c}) &= \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(s)}_{j,k,c,L}} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\widehat{\boldsymbol{\beta}}^{(s)}_{h,j,k,c,L}\bigr)^2 + \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(s)}_{j,k,c,R}} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\widehat{\boldsymbol{\beta}}^{(s)}_{h,j,k,c,R}\bigr)^2 \\ &\quad + \sum_{j' \neq j}\;\sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(s-1)}_{j'}} \bigl(\mathrm{IV}(\mathbf{z}_{i,\ell})-\mathbf{x}_{h,i,\ell}'\widehat{\boldsymbol{\beta}}^{(s-1)}_{h,j'}\bigr)^2. \end{align*}\]

Among all candidate splits \((j,k,c)\) we select the one that produces the largest reduction in total loss, i.e. \[(j_s,k_s,c_s) \;=\; \arg\min_{(j,k,c)} \mathcal{L}\bigl(\mathcal{T}^{(s)}_{j,k,c}\bigr),\] and set \(\mathcal{T}^{(s)} = \mathcal{T}^{(s)}_{j_s,k_s,c_s}\), where the two child leaves use \(\widehat{\boldsymbol{\beta}}^{(s)}_{h,j_s,k_s,c_s,L}\) and \(\widehat{\boldsymbol{\beta}}^{(s)}_{h,j_s,k_s,c_s,R}\), and all other leaves keep \(\widehat{\boldsymbol{\beta}}^{(s-1)}_{h,j'}\). This procedure is repeated until the number of leaves reaches a pre-specified optimal tree size.

To select the optimal tree size, denoted by \(J^{\ast}_h\), we use the first half of the full sample and further split it into training and validation subsamples. We consider candidate trees with at most \(15\) terminal leaves and impose a minimum leaf size of \(10{,}000\) option observations. Among these candidates, we choose \(J^{\ast}_h\in\{2,\ldots,15\}\) using the CART cost-complexity criterion based on IV forecasting errors in the validation subsample, as described in Appendix B.

Given the selected tree \(\widehat{\mathcal{T}}\) and coefficients \(\{\widehat{\boldsymbol{\beta}}_{h,j}(\widehat{\mathcal{T}})\}\), the \(h\)-step-ahead Tree-SHAR IV forecast for any option with characteristics \(\mathbf{z}_{i,t+h}\) is obtained by assigning it to its leaf under \(\widehat{\mathcal{T}}\) and applying the corresponding local HAR coefficients, \[\begin{equation} \widehat{\mathrm{IV}}^{\mathrm{Tree}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) \;=\; \mathbf{x}_{h,i,t+h}'\,\widehat{\boldsymbol{\beta}}_{h,j(i,t+h)}(\widehat{\mathcal{T}}),\tag{10} \end{equation}\] which is the direct analogue of the SHAR forecast in (6), but with coefficients that adapt across the moneyness–maturity surface.

2.4 Bagged Tree-SHAR model

The Tree-SHAR model in Section 2.3 introduces cross-sectional heterogeneity by estimating a single tree partition and a set of leaf-specific SHAR coefficients. While this approach improves forecasting performance and remains interpretable, a single estimated tree can be sensitive to sampling variation in the data. To reduce this estimation noise, we construct a bootstrap-aggregated (“bagged”) ensemble of Tree-SHAR models following Breiman, 1996. The procedure is closely related to random forests Breiman, 2001, with one important distinction: we do not randomly subsample the set of splitting variables at each split. Instead, every tree uses the same two splitting variables, moneyness and maturity.

We fix the target number of terminal leaves, denoted by \(J^{\ast}_h\), at the value selected for Tree-SHAR using the CART cost-complexity criterion described in Appendix B. We then construct an ensemble of \(B=100\) trees as follows.

Step 1: Bootstrap resampling. For each \(b=1,\ldots,B\), we draw a bootstrap sample \(\mathcal{D}^{[b]}_{t}\) from \(\mathcal{D}_{t}\) in (7) by sampling observations independently with replacement.⁷ Each bootstrap observation therefore retains the full triplet consisting of the observed implied volatility, its associated option characteristics, and its SHAR regressor vector.

Step 2: Estimate a fixed-size Tree-SHAR on each bootstrap sample. For each bootstrap sample \(\mathcal{D}^{[b]}_{t}\), we estimate a fixed-size tree \(\mathcal{T}^{[b]}\) with \(J^{\ast}_h\) terminal leaves and the corresponding leaf-specific SHAR coefficients \(\{\widehat{\boldsymbol{\beta}}^{[b]}_{h,j}\}_{j=1}^{J^{\ast}_h}\) using the same greedy splitting algorithm as in Section 2.3. For all bootstrap trees, we restrict the splitting variables to moneyness and maturity and continue splitting until the tree reaches \(J^{\ast}_h\) leaves, subject to the same minimum leaf-size constraint imposed in Tree-SHAR: each terminal leaf must contain at least \(10{,}000\) option observations. When evaluating a candidate split, we re-estimate the SHAR coefficients only in the two candidate child nodes, while holding the coefficient estimates in all unaffected leaves fixed at their values in the current tree before the split.

Step 3: Bagged forecasts by averaging across trees. To forecast the implied volatility at date \(t+h\) for option \(i\), conditional on its target moneyness–maturity point \(\mathbf{z}_{i,t+h}\), each bootstrap tree \(\mathcal{T}^{[b]}\) assigns the option to a terminal leaf \(j^{[b]}(i,t+h)\). Given this leaf assignment, the corresponding leaf-specific HAR coefficients yield the \(h\)-step-ahead forecast \[\widehat{\mathrm{IV}}^{[b]}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) = \mathbf{x}_{h,i,t+h}'\,\widehat{\boldsymbol{\beta}}^{[b]}_{h,j^{[b]}(i,t+h)}.\] The Bagged Tree-SHAR forecast is obtained by averaging these tree-specific forecasts: \[\begin{equation} \widehat{\mathrm{IV}}^{\mathrm{Bag}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) = \frac{1}{B}\sum_{b=1}^{B} \widehat{\mathrm{IV}}^{[b]}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}).\tag{11} \end{equation}\]

This aggregation preserves the interpretability of Tree-SHAR: each component tree partitions the moneyness–maturity surface and estimates a local HAR specification within each terminal leaf. At the same time, averaging across bootstrap trees reduces the sensitivity of the forecast to any single sample or tree partition.

2.5 Boosted Tree-SHAR model

This subsection introduces a boosted version of Tree-SHAR. The construction follows the gradient boosting framework of Friedman, 2001 (see also Hastie et al., 2009), which builds an additive predictor by repeatedly fitting a tree-based model to the current fitting error and updating the forecast in a shrunken step in that direction. In our setting, each boosting step fits a local SHAR regression to the remaining errors, allowing the HAR coefficients to be refined sequentially across the moneyness–maturity surface.

Boosting is useful in this setting because a single Tree-SHAR model must capture all heterogeneity using one partition of the option surface. The IV surface, however, may contain several distinct nonlinear patterns across moneyness and maturity. Rather than relying on one partition to capture all such patterns, the Boosted Tree-SHAR model adds trees sequentially, with each new tree targeting the errors left by the previous additive fit. Under squared-error loss, this procedure corresponds to a regularized gradient-descent improvement of the forecasting function (see, e.g., Chapter 10 in Hastie et al., 2009).

Fix the forecast horizon \(h\). Let \(\mathcal{D}_{t}\) be the sample available at forecast origin \(t\), as defined in (7). The initial forecast function is the baseline SHAR fit,\[F_{0}\!\left(\mathbf{z}_{i,\ell}\right) = \mathbf{x}_{h,i,\ell}'\,\widehat{\boldsymbol{\beta}}_{h|t},\] where \(\widehat{\boldsymbol{\beta}}_{h|t}\) is the pooled OLS estimate from (4) using \(\mathcal{D}_{t}\).

For boosting iteration \(m=1,\ldots,M\), we first compute the current fitting error, \[e^{(m)}_{i,\ell} = \mathrm{IV}(\mathbf{z}_{i,\ell})- F_{m-1}\!\left(\mathbf{z}_{i,\ell}\right), \qquad (i,\ell)\in\mathcal{D}_{t}.\] We then fit a Tree-SHAR model to these errors using the same splitting variables \(\mathbf{z}_{i,\ell}\) as in Section 2.3. Let \(\mathcal{T}^{(m)}\) denote the tree fitted at iteration \(m\), and let \(\{R^{(m)}_j\}_{j=1}^{J^{\ast}_h}\) denote its terminal leaves. The leaf assignment for observation \((i,\ell)\) is denoted by \(j^{(m)}(i,\ell)\), so that \(\mathbf{z}_{i,\ell}\in R^{(m)}_{j^{(m)}(i,\ell)}\).

Within each terminal leaf \(j\), we estimate a local SHAR correction by regressing the current fitting errors on the SHAR regressor vector: \[\begin{equation} \widehat{\boldsymbol{\delta}}^{(m)}_{h,j} = \arg\min_{\boldsymbol{\delta}} \sum_{(i,\ell):\,\mathbf{z}_{i,\ell}\in R^{(m)}_j} \left( e^{(m)}_{i,\ell} - \mathbf{x}_{h,i,\ell}'\boldsymbol{\delta} \right)^2 .\tag{12} \end{equation}\] This produces the tree-based correction \[g_m\!\left(\mathbf{z}_{i,\ell}\right) = \mathbf{x}_{h,i,\ell}'\, \widehat{\boldsymbol{\delta}}^{(m)}_{h,j^{(m)}(i,\ell)}.\] The boosted forecast function is updated using a learning rate \(\nu\in(0,1]\): \[\begin{equation} F_m\!\left(\mathbf{z}_{i,\ell}\right) = F_{m-1}\!\left(\mathbf{z}_{i,\ell}\right) + \nu\,g_m\!\left(\mathbf{z}_{i,\ell}\right).\tag{13} \end{equation}\]

After \(M\) boosting iterations, the Boosted Tree-SHAR forecast of implied volatility at date \(t+h\) for option \(i\), conditional on its target moneyness–maturity point \(\mathbf{z}_{i,t+h}\), is given by \[\begin{equation} \widehat{\mathrm{IV}}^{\mathrm{Boost}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h}) = F_M\!\left(\mathbf{z}_{i,t+h}\right) = \mathbf{x}_{h,i,t+h}'\,\widehat{\boldsymbol{\beta}}_{h|t} + \nu\sum_{m=1}^{M} \mathbf{x}_{h,i,t+h}'\, \widehat{\boldsymbol{\delta}}^{(m)}_{h,j^{(m)}(i,t+h)} .\tag{14} \end{equation}\]

Tree fitting within each boosting iteration follows the same greedy splitting logic as in Tree-SHAR, except that candidate splits are evaluated using the residual-regression criterion implied by (12). In the empirical implementation, we keep each tree component deliberately simple. We set the number of terminal leaves equal to \(J^{\ast}_h\), the value selected for Tree-SHAR using the CART cost-complexity criterion described in Appendix B, and impose the same minimum leaf-size constraint as in Tree-SHAR: each terminal leaf must contain at least \(10{,}000\) option observations. For the boosting regularization, we set the learning rate to \(\nu=0.1\) and the number of boosting iterations to \(M=100\), which are standard benchmark choices under squared-error loss and provide a stable bias–variance trade-off in practice.

This boosted specification differs from a standard boosted regression tree in one important respect. In a usual boosted tree, each terminal leaf contributes a constant fitted value. Here, each terminal leaf contributes a local HAR regression, summarized by a vector of coefficient adjustments \(\widehat{\boldsymbol{\delta}}^{(m)}_{h,j}\). As a result, the contribution of each boosting tree depends not only on the moneyness–maturity region assigned by \(\mathbf{z}_{i,\ell}\), but also on the SHAR regressor vector \(\mathbf{x}_{h,i,\ell}\). The model therefore preserves the HAR interpretation within each region while allowing the HAR coefficients to vary flexibly over the implied-volatility surface through repeated, region-specific coefficient adjustments.

3 Data

We use daily S&P 500 out-of-the-money (OTM) call and put options between January 2011 and August 2023 obtained from OptionMetrics. Due to the large cross-sectional size of the options data, which include less reliable data points, we use a filtered dataset, mainly following Dufays et al., 2025. We keep (i) options with maturities between 20 and 240 days, (ii) options with implied volatility above 5% and below 150%, (iii) options with moneyness (defined as underlying price/strike) between 0.8 and 1.6, (iv) options with mid price and bid price at least $0.50 and $0.375, respectively, and (v) options not subject to fundamental data errors, where the bid price exceeds the offer price or where a negative price is implied by put-call parity. We further exclude options whose volume and open interest are less than 100 contracts to ensure we include only sufficiently liquid contracts. The final filtered dataset comprises 786,867 option contracts.

We obtain the time series of S&P 500 index prices and returns from CRSP. The risk-free rate is proxied by the one-month Treasury bill rate obtained from CRSP as well. We download the VIX data from Bloomberg. Following the standard implementation in the literature, options are valued using a maturity-specific risk-free rate obtained from OptionMetrics.

4 Empirical Results

In this section, we discuss the forecasting performance of tree-based SHAR models. Section 4.1 reports the overall forecasting results and compare different models, especially between models with global homogeneous specifications and tree-based localized models. Section 4.2 examines cross-sectional performance variations across moneyness and maturity dimensions. Lastly, Section 4.3 illustrates the estimated tree structures on two dates representing good and bad periods. We further discuss the leaf-specific model coefficient estimates and their economic interpretations.

4.1 IV forecasting performance

This subsection evaluates the out-of-sample implied-volatility (IV) forecasting performance of the baseline Surface HAR (SHAR) model and its tree-based extensions proposed in Section 2. We assess the models’ ability to predict the future implied volatilities over horizons of \(h \in \{1, 5, 22\}\) forecasting days.

We evaluate forecasting performance over the OOS period 2018–2023 using an expanding-window scheme. For each forecast origin \(t\) in the OOS period, every model is estimated using only data available up to date \(t\), and the resulting estimates are used to forecast implied volatility at date \(t+h\). For the SHAR benchmark, this daily re-estimation updates the global HAR coefficients. For the tree-based specifications, this daily re-estimation updates both the tree partition structure—including splitting thresholds and leaf assignments—and the leaf-specific HAR coefficients. Thus, the Tree-SHAR, Bagged-SHAR, and Boosted-SHAR forecasts are generated from tree structures and coefficient estimates that are re-optimized at each OOS forecast origin.

We compare the proposed machine-learning specifications against three benchmarks:

1. Random Walk (RW): Assumes the future surface equals the current fitted surface.

2. SHAR: The baseline model of Dufays et al., 2025 with a single set of HAR coefficients for the entire surface.

3. SHAR with \(3{\times}3\) Grids: An ad-hoc local model that partitions the moneyness–maturity space into nine fixed rectangular bins and estimates separate HAR coefficients within each bin. The bins are constructed by dividing each dimension into three empirical terciles: the moneyness cutoffs are the \(1/3\) and \(2/3\) sample quantiles of all observed moneyness values, and the maturity cutoffs are the \(1/3\) and \(2/3\) sample quantiles of all observed maturities. This produces \(3{\times}3=9\) fixed regions and serves as a control to distinguish the benefits of any localization from the benefits of data-driven partitioning.

The complexity of the tree models is determined via the cost-complexity pruning described in Appendix B. The algorithm selects parsimonious structures, \(J^{\ast}_{1}=6\) leaves for the one-day horizon, \(J^{\ast}_{5}=8\) leaves for the five-day horizon, and increasing to \(J^{\ast}_{22}=11\) leaves for the one-month horizon (\(h=22\)).

Table 1 reports the OOS Root Mean Squared Error (RMSE) in percentage points, revealing a strict hierarchy of model performance that holds across all forecast horizons. First, incorporating local heterogeneity is uniformly superior to global modeling. At the one-day horizon (\(h=1\)), SHAR and Random Walk models produce aggregate RMSEs of roughly 2.17. Introducing the fixed \(3{\times}3\) partition reduces this to 2.119. However, the data-driven Tree-SHAR significantly outperforms the fixed grid, lowering the RMSE to 1.985. This suggests that the optimal boundaries for volatility dynamics do not align with arbitrary moneyness–maturity cutoffs but are better identified algorithmically. Second, ensemble learning delivers substantial incremental gains. Bagging the Tree-SHAR (Bagged-SHAR) reduces the RMSE to 1.965, while gradient boosting (Boosted-SHAR) achieves the lowest error of 1.852. The superiority of boosting is consistent with its ability to sequentially correct biases in the surface fit, rather than simply reducing variance through averaging. The economic magnitude of these improvements increases with the forecast horizon. At the one-month horizon (\(h=22\)), the SHAR model yields an RMSE of 5.910. The Tree-SHAR and Boosted-SHAR specifications reduce this value to 5.323 and 5.124, respectively. The Boosted-SHAR result represents a reduction in forecast error of approximately 13% relative to the global SHAR benchmark and over 18% relative to the Random Walk (RMSE 6.248).

The year-by-year decomposition in Table 1 further highlights the robustness of the tree-based methods during the extreme volatility period of 2020 (the onset of the COVID-19 pandemic). In 2020, all models experience elevated errors. However, the adaptability of the tree partitions proves crucial. For the 22-day horizon, the SHAR model suffers an RMSE of 12.177 in 2020. The Tree-SHAR and Boosted-SHAR models mitigate this error to 11.091 and 10.745, respectively. These outcomes indicate that when the implied volatility surface becomes steep and distorted during crises, the ability to isolate specific surface regions (such as deep OTM or short-term options) allows for more accurate dynamic modeling than a rigid global specification. The same argument holds for shorter-horizon forecasts (\(h=1\) or \(5\)).⁸

Table 2 formalizes these comparisons using Diebold-Mariano Diebold and Mariano, 1995 tests and the Model Confidence Set Hansen et al., 2011 confirming that the observed reductions in error are statistically significant. The DM statistics (comparing row vs. column) are positive and large across the board, rejecting the null hypothesis of equal predictive accuracy in favor of the more sophisticated models. For example, at \(h=22\), the test statistic comparing SHAR against Boosted-SHAR is 92.48, providing overwhelming evidence against the global model. Furthermore, the tests confirm that the gains from machine learning are statistically distinct: Boosting outperforms Bagging (statistic 45.98), and Bagging outperforms the single Tree (statistic 30.13). Finally, the Model Confidence Set results (last column of Table 2) are definitive. For every forecast horizon (\(h=1, 5, 22\)), the Boosted-SHAR model is the only model retained in the confidence set at the 95% level. All other specifications—including the single Tree-SHAR and the Bagged-SHAR—are eliminated. This confirms that the stagewise additive expansion provided by gradient boosting captures complex nonlinearities in the evolution of the implied volatility surface that other methods miss.

In summary, the empirical evidence demonstrates that allowing HAR dynamics to vary across algorithmically selected regions of the moneyness–maturity space systematically improves forecasting accuracy.⁹ The data-driven partitions outperform ad-hoc grids, and gradient boosting provides the most effective method for aggregating these local dynamics.

Panel A of Figure 3 provides a time series perspective on these forecast gains by plotting the daily out-of-sample RMSEs of the SHAR and Tree-SHAR models for the one-day horizon. The figure shows that the Tree-SHAR error series is almost uniformly below the SHAR error series throughout the evaluation period. This pattern indicates that the improvement from tree-based localization is not driven by a small number of isolated dates, but instead reflects a persistent reduction in forecast errors over time. The gains are visible both during relatively calm periods and during episodes of market stress, suggesting that the tree partition improves forecast accuracy systematically rather than only in crisis periods.

The advantage of Tree-SHAR is especially pronounced on days with large forecasting errors. Whenever daily RMSEs spike, the Tree-SHAR model typically produces a smaller error than the baseline SHAR model. This suggests that the tree partition is particularly valuable when the implied-volatility surface is harder to forecast, because it allows different regions of the surface to respond through region-specific HAR coefficients. In other words, during periods when a global SHAR specification generates large errors, the localized SHAR coefficients help mitigate those errors by adapting the forecasting dynamics across the moneyness–maturity surface.

Panels B and C of Figure 3 show that the same pattern holds across longer forecast horizons. For \(h=5\) and \(h=22\), the Tree-SHAR RMSE series remains generally below the SHAR series, and the gap is again most visible around dates with large forecast errors. Thus, the benefits of tree-based localization are not confined to one-day-ahead forecasting. They also appear at weekly and monthly horizons, where the implied-volatility surface is more difficult to predict and where regional differences in HAR dynamics become especially important.

Although Figure 3 reports only SHAR and Tree-SHAR for visual clarity, the aggregate results in Table 1 show that the ensemble extensions further improve on the single-tree specification. Bagged Tree-SHAR reduces the sensitivity of forecasts to any single tree partition, while Boosted Tree-SHAR sequentially corrects remaining forecast errors. Therefore, the figures should be interpreted as conservative illustrations of the gains from machine-learning-based localization: even the single Tree-SHAR model already dominates the global SHAR benchmark over time, and the bagged and boosted versions deliver additional reductions in forecast error.

4.2 Cross-sectional variations in IV forecasting performance

To further investigate the sources of forecasting improvements, we examine how the tree-based models perform across different regions of the moneyness-maturity space. Tables 3 and 4 decompose out-of-sample forecasting errors for the one-day horizon across moneyness and maturity buckets, respectively. At each point in time, we partition the sample into three groups based on the 20th and 80th percentiles of the respective dimension. This cross-sectional analysis reveals economically meaningful patterns that help explain why localized SHAR dynamics deliver superior forecasting performance.

In the moneyness dimension (Table 3), the benefits of tree-based localization are most pronounced at the extremes of the moneyness distribution. For OTM calls (Moneyness \(\leq q_{20}\)) and deep OTM puts (Moneyness \(\geq q_{80}\)), Boosted-SHAR consistently outperforms all competing specifications across all years. Tree-SHAR also dominates the Random Walk, SHAR, and SHAR with \(3\times 3\) grids benchmarks in these regions. These OTM options are particularly challenging to forecast using globally homogeneous models, as they exhibit distinct volatility dynamics driven by tail-risk concerns, crash premia, or jackpot expectations that differ fundamentally from NTM options. The localized HAR coefficients estimated by tree-based methods successfully capture these region-specific patterns, leading to substantial forecast improvements.

For NTM options and slightly OTM puts (\(q_{20} < \text{Moneyness} < q_{80}\)), Boosted-SHAR continues to perform best in most years, with exceptions in 2020 and 2022. During the 2020 market stress period, the Random Walk assumption actually produces better short-horizon forecasting accuracy for NTM IVs compared to dynamic forecasting models. This pattern suggests that when the volatility surface undergoes rapid and dramatic shifts, high persistence dominates over dynamic patterns for NTM options. Importantly, the differences among model specifications in this moneyness bucket are substantially smaller than in the extreme moneyness regions, indicating that the Random Walk already performs reasonably well for NTM options. Nevertheless, the overall performance still improves when we apply tree-based SHAR approaches, particularly at longer horizons (see Tables A.4 and A.6 in Appendix). The outperformance of the Random Walk specification in the NTM bucket during 2020 remains confined to 5-day horizon forecasting, while it underperforms tree-based models in 22-day horizon forecasting. This result implies that during highly volatile periods, the very high persistence of NTM IVs is well-captured by the Random Walk over short horizons; the economic implication embedded in the HAR structure becomes valuable as the forecast horizon extends.

The largest improvements occur at the tails of the moneyness distribution when we employ Boosted-SHAR. Under the Random Walk or global SHAR specifications, OTM call and deep OTM put IVs exhibit substantially larger forecasting errors than NTM options. However, this pattern reverses when we use Boosted-SHAR. Specifically, the improvements of Boosted-SHAR relative to the global SHAR benchmark are approximately 30% and 26% for OTM calls and deep OTM puts, respectively, while the improvement in the NTM region is only about 2%. This stark contrast confirms that allowing for heterogeneous volatility dynamics across the moneyness dimension is particularly valuable in regions where global models struggle and where the tree algorithm identifies the need for separate local specifications.

In the maturity dimension (Table 4), we observe a different pattern. Overall, IV forecast errors are larger for shorter-maturity options regardless of the forecasting approach. This baseline observation reflects the higher volatility of volatility for near-term contracts. Nevertheless, Boosted-SHAR performs best across all maturity buckets and all years for one-day horizon forecasts. The same pattern persists at longer horizons (\(h=5\) and \(22\)), with only a few isolated cases where Bagged-SHAR achieves marginally better performance than Boosted-SHAR (see Tables A.5 and A.7).

Taken together, the cross-sectional results demonstrate that the tree-based partitioning successfully identifies economically meaningful regions of the implied volatility surface where distinct forecasting dynamics prevail. The localized SHAR parameters adapt to the unique persistence and mean-reversion characteristics of different option segments, thereby improving forecast accuracy most substantially in deep OTM and short-maturity regions.

4.3 Heterogeneity of the volatility surface

To understand how the tree-based approaches achieve their forecasting improvements, we examine Tree-SHAR’s estimated surface partitions and the corresponding local HAR dynamics by focusing on one-day forecast horizon.¹⁰ Figures 4 and 5 illustrate representative Tree-SHAR partitions estimated for two contrasting dates: January 4, 2018 (a calm market) and March 17, 2020 (a stressed market). Table 5 reports the estimated HAR coefficients for both the global SHAR model and the leaf-specific parameters from Tree-SHAR on these dates. These data-driven partitions, by and large, align with the regions where we observed the noticeable forecast improvements in Tables 3 and 4. This alignment confirms that the tree algorithm identifies economically meaningful surface segments where a localized model outperforms global specifications.

In Figure 4, Panel A shows the tree structure containing six terminal leaves, and Panel B maps these leaves onto the moneyness-maturity space, during a stable market state. The tree algorithm first isolates deep OTM calls (Leaf 1: Moneyness \(< 0.96\)), which constitute 10% of observations. This region experiences distinct volatility dynamics, as evidenced by the substantially different HAR coefficients in Table 5. While the global SHAR model estimates a daily persistence coefficient of 0.909, Leaf 1 exhibits a notably lower daily coefficient of 0.744 and a negative monthly component of -0.096. In other words, deep OTM call IVs respond relatively less to yesterday’s level and exhibit faster mean reversion than the surface average implied in the SHAR model. Deep OTM call options during good times are likely to be held by speculators seeking leveraged directional exposure. Although our tree structure does not particularly divide between short- and long-term OTM calls, Eaton et al., 2026 show that retail investors exhibit strong lottery-like preferences (e.g., OTM calls) and concentrate on trading short-term options. If IVs in Leaf 1 are primarily driven by such trading activities, we conjecture that option prices rapidly incorporate new information that is not captured by historical IVs and do not simply inherit yesterday’s level, hence relatively lower daily coefficient, \(\hat{\beta}^{(d)}_{h=1}\).

The tree then partitions the remaining surface along maturity and moneyness dimensions, creating economically interpretable regions that correspond to different investor clienteles and hedging demands. Leaf 6 (Moneyness \(\geq 0.96\), Maturity \(\geq\) 140 days) captures the broad region of long-maturity options spanning NTM to deep OTM puts, representing the segment where institutional investors typically seek longer-term portfolio protection and hedging exposure. This region displays distinctive dynamics in Table 5. The daily persistence coefficient is only 0.412, but weekly and monthly coefficients are 0.442 and 0.279. In comparison, other leaves containing shorter-maturity options show much higher daily persistence and lower weekly or monthly persistence. This structure indicates that while these long-dated options do not simply inherit yesterday’s volatility level, they exhibit substantial persistence at weekly and monthly horizons, consistent with the observation in Figure 2. A potential interpretation is that long-term hedging demand responds to persistent shifts in risks and macroeconomic uncertainty that evolve gradually over time, rather than reacting primarily to daily market fluctuations.

The shorter-maturity segment (Maturity \(<\) 39 days) is further partitioned by moneyness into Leaves 2, 3, and 4. Leaves 2 and 3, which cover NTM and slightly OTM puts, exhibit very high daily persistence with coefficients of 0.907 and 0.996, respectively. Leaf 4 captures short-dated deep OTM puts and shows an even higher daily persistence coefficient of 1.010. These short-maturity regions account for substantial trading volume from both hedgers and market makers, and their IVs track relatively more recent market conditions closely. The uniformly high daily persistence of IV across short-term options suggests that these IVs are strongly anchored to yesterday’s level, consistent with short-horizon market expectations that evolve smoothly day-to-day.¹¹

Figure 5 shows how the tree partition adapts dramatically during the 2020 COVID-19 crisis period. While the tree still produces six leaves, the splitting thresholds and structure differ markedly from the good period. The tree algorithm now isolates deep OTM calls at a higher moneyness threshold (Moneyness \(<\) 0.92 for Leaf 1, compared to 0.96 in Figure 4), creating a more extreme partition. The moneyness cutoffs for other leaves shift leftward overall, implying more heterogeneous dynamics in OTM calls and NTM options but more homogeneous dynamics in a wide range of short-term OTM puts. In fact, the finer partitions in OTM put regions during good times align with empirical observations that investors are more fearful of negative market shocks during good periods, differentiating among various tail-risk scenarios and embedding more heterogeneous information in OTM puts Kim, 2025; Schreindorfer and Sichert, 2025. Conversely, during bad times, the coarser short-term OTM put partitions suggest that investors respond more uniformly to crash-risk concerns.

The estimated HAR coefficients in Table 5 for March 17, 2020 present even more pronounced heterogeneity than the good period. The global SHAR model now estimates substantially lower daily persistence (0.843 versus 0.909 in 2018) and a sizable positive weekly coefficient (0.191 versus -0.012 in 2018), reflecting the general increase in the importance of weekly dynamics during the crisis. However, this global pattern masks extreme variation across the surface.

Leaf 1, which captures deep OTM calls, now displays radically different dynamics compared to the deep OTM calls segment during good times (Leaf 1 in Figure 4). The daily persistence coefficient is only 0.277, but a weekly coefficient is 0.892, and a monthly coefficient is largely negative, -0.831. This structure implies that during extreme stress, such speculative call options’ IVs depend far more on the weekly moving average of recent surfaces than on yesterday’s observation alone. The highly negative monthly coefficient suggests sharp mean reversion from any sustained elevation of volatilities. We conjecture that this pattern aligns with a market where speculative demand for lottery-like upside exposure persists for some days but not a long period of time, along with market expectations that a crisis may persist for some period.

Similarly, Leaf 6 (long-maturity OTM puts and NTM options) shows a daily coefficient of only 0.532 but a weekly coefficient of 0.635 during the crisis, again emphasizing the shift away from daily toward weekly persistence. However, its monthly persistence disappears in contrast to the sizable monthly persistence observed for long-term OTM put IVs during good times (-0.008 versus 0.279 in 2018). On the other hand, Leaf 4 (short-maturity OTM puts) maintains very high daily persistence of 0.939 even during the crisis, implying high persistence of short-horizon expectation of crash risk. All in all, this heterogeneity explains why the tree-based approach substantially outperforms global models during stressed periods: no single set of HAR coefficients can adequately describe the divergent dynamics across the surface when market conditions fracture the usual relationships.

Across our OOS evaluation period, we find that the tree structure remains relatively stable under similar market conditions but adjusts noticeably following major market shocks or volatility regime shifts. This combination of persistence in tree structure (reducing estimation noise over time) and adaptability when needed (capturing evolution of market segmentation) contributes to the tree-based models’ robust OOS forecasting performance. The localized HAR specifications effectively allow short-term options, long-term options, and deep OTM calls and puts to exhibit distinct dynamics reflecting their specialized clienteles. A single global set of coefficients cannot accommodate such adjustments, but the data consistently identify these adjustments as improving forecast accuracy across the surface.

5 Conclusion

Standard implied volatility forecasting models typically impose homogeneous dynamics across the entire option surface. While this assumption is convenient, it overlooks the substantial heterogeneity in volatility dynamics across moneyness and maturity. In this paper, we propose a machine-learning approach that uses regression trees to identify data-driven regions of the moneyness–maturity space and estimate local forecasting dynamics within each region. Building on the SHAR framework, our Tree-SHAR model preserves the interpretability of HAR-type implied volatility dynamics while allowing the coefficients to vary across economically meaningful regions of the option surface.

Using daily S&P 500 index options, we show that incorporating cross-sectional heterogeneity substantially improves out-of-sample implied volatility forecasts. Tree-SHAR outperforms the global SHAR benchmark and the Random Walk benchmark across forecast horizons. It also improves on an ad-hoc \(3\times 3\) grid specification, indicating that the gains are not simply due to local estimation, but rather to the data-driven selection of moneyness–maturity partitions. The ensemble extensions further strengthen these results. Bagged Tree-SHAR reduces the sensitivity of forecasts to any single estimated partition, while Boosted Tree-SHAR delivers the lowest forecast errors overall by sequentially correcting remaining prediction errors. These gains are statistically significant and are observed under a daily expanding-window forecasting exercise in which both the tree structures and the local coefficients are re-estimated at each out-of-sample forecast origin.

The forecasting improvements are economically informative. The gains from tree-based localization are especially large in regions where global models perform poorly, such as OTM calls, deep OTM puts, and short-maturity options. They are also pronounced during periods with large forecast errors, including the COVID-19 stress period, but remain visible during calmer market conditions. The estimated tree partitions and leaf-specific coefficient estimates show that different regions of the option surface display distinct persistence patterns. Short-dated options tend to place more weight on recent fitted surfaces, whereas longer-dated options exhibit stronger weekly and monthly components. These findings suggest that the implied-volatility surface is governed by heterogeneous dynamic behavior that cannot be adequately captured by a single global set of coefficients.

Overall, the results highlight the value of combining interpretable financial forecasting models with machine-learning-based localization. Regression trees provide a flexible but transparent way to determine where local models should be estimated, while the SHAR structure maintains a clear economic interpretation of daily, weekly, and monthly volatility persistence. One natural extension is to enrich the cross-sectional splitting variables used to partition the option surface. Beyond moneyness and maturity, additional option-level characteristics such as option Greeks, trading volume, open interest, or bid–ask spreads may help identify more refined forms of heterogeneity across the option panel. A second direction is to extend the tree-based approach from the cross-sectional dimension to the time-series dimension. For example, one could use market-wide state variables such as the VIX, past returns, realized volatility, or other measures of financial conditions to determine when forecasting dynamics change over time. This would connect our framework to the broader literature on tree-structured time-series models, which uses observed state variables to identify regimes with distinct dynamic parameters Audrino and B\"uhlmann, 2001; Patton and Simsek, 2026. Combining cross-sectional option characteristics with time-series state variables may provide a richer framework for modeling how implied-volatility dynamics vary both across the option surface and across market conditions.

6 Figures and Tables

Figure 1: Implied Volatilities by Moneyness and Maturity in Data
Panel A: IV by Moneyness
Panel B: IV by Maturity
Notes: This figure depicts the average implied volatilities (IV) by moneyness and maturity groups over time from January 2018 to August 2023. On each date, three groups are classified such that the moneyness or maturity is less than or equal to its 20th percentile (blue solid line), between the 20th and 80th percentiles (red dashed line), and greater than or equal to the 80th percentile (black dotted line). Panel A shows IV by moneyness and Panel B shows IV by maturity. All IVs are expressed in percentage points.

Accessible version

Figure 2: Autocorrelation Function of Implied Volatilities
Panel A: IV Level by Moneyness
Panel B: IV Level by Maturity
Panel C: ΔIV by Moneyness
Panel D: ΔIV by Maturity
Notes: This figure depicts the autocorrelation functions (ACF) of implied volatilities (IVs) and their first differences up to 22 lags. Sample ACFs are computed using the average IVs within three groups classified such that the moneyness or maturity is less than or equal to its 20th percentile (blue bar), between the 20th and 80th percentiles (green bar), and greater than or equal to the 80th percentile (orange bar). Panels A and B plot the ACF of IVs by moneyness and maturity groups. Panels C and D plot the ACF of daily changes (first differences) in IVs. The red dashed lines represent the 95% confidence intervals.

Accessible version

Figure 3: Daily Out-of-Sample RMSEs for SHAR and Tree-SHAR Across Forecast Horizons
Panel A: h = 1
Panel B: h = 5
Panel C: h = 22
Notes: This figure plots the daily out-of-sample RMSEs of implied-volatility forecasts for the SHAR and Tree-SHAR models across forecast horizons \(h=1\), \(5\), and \(22\). Panel A reports the one-day-ahead horizon, Panel B reports the five-day-ahead horizon, and Panel C reports the 22-day-ahead horizon. The blue solid line reports the daily RMSE of the baseline SHAR model, while the red dashed line reports the daily RMSE of the Tree-SHAR model. Daily RMSEs are computed across all options observed on each forecast target date and are expressed in percentage points of implied volatility.

Accessible version

Figure 4: Estimated Tree-SHAR Partition on January 4, 2018: h = 1
Panel A: Estimated Tree Structure
Panel B: Induced Moneyness–Maturity Partition
Notes: This figure illustrates the estimated Tree-SHAR partition for the one-day-ahead forecast horizon on January 4, 2018, using AHBS-fitted implied-volatility surfaces. Panel A reports the estimated regression tree. Each internal node shows the selected splitting variable and threshold, while each terminal leaf reports the number of option observations assigned to that leaf and the corresponding sum of squared errors (SSE). Panel B maps the same terminal leaves into the moneyness–maturity space. The figure shows how the Tree-SHAR algorithm partitions the option surface into economically interpretable regions and estimates separate local HAR coefficients within each region.

Accessible version

Figure 5: Estimated Tree-SHAR Partition on March 17, 2020: h = 1
Panel A: Estimated Tree Structure
Panel B: Induced Moneyness–Maturity Partition
Notes: This figure illustrates the estimated Tree-SHAR partition for the one-day-ahead forecast horizon on March 17, 2020, using AHBS-fitted implied-volatility surfaces. Panel A reports the estimated regression tree. Each internal node shows the selected splitting variable and threshold, while each terminal leaf reports the number of option observations assigned to that leaf and the corresponding sum of squared errors (SSE). Panel B maps the same terminal leaves into the moneyness–maturity space. Compared with Figure 4, the partition highlights how the tree structure adapts during the COVID-19 stress period by isolating more extreme regions (Leaf 1: deep OTM calls) of the moneyness–maturity surface.

Accessible version

Table 1: Out-of-Sample Implied Volatility Forecasting Performance (RMSE in %)

\(h=1\)

\(h=5\)

\(h=22\)

Table 2: Pairwise Diebold-Mariano Test Statistics and Model Confidence Set (MCS) for SHAR

Notes: This table reports the pairwise Diebold-Mariano (DM) test statistics comparing the out-of-sample forecasting errors of the competing models across 1, 5, and 22-day horizons. A negative t-statistic indicates that the model in the row outperforms (has lower error than) the model in the column, while a positive value indicates that the column model outperforms the row model. Statistical significance at the 1% level corresponds to an absolute t-statistic greater than 2.33 (\(|t| > 2.33\)). The table also indicates with a checkmark the models that are included in the Model Confidence Set (MCS) at the 95% confidence level.

\(h = 1\)

\(h = 5\)

\(h = 22\)

Table 3: Out-of-Sample Implied Volatility Forecasting Performance by Moneyness (RMSE in %, \(h=1\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of one-day horizon (\(h=1\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three moneyness groups. On each date, three groups are classified such that the moneyness is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Moneyness \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Moneyness \(<\) \(q_{80}\)

Moneyness \(\geq\) \(q_{80}\)

Table 4: Out-of-Sample Implied Volatility Forecasting Performance by Maturity (RMSE in %, \(h=1\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of one-day horizon (\(h=1\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three maturity groups. On each date, three groups are classified such that the maturity is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Maturity \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Maturity \(<\) \(q_{80}\)

Maturity \(\geq\) \(q_{80}\)

Table 5: SHAR and Tree-SHAR Coefficient Estimates on Representative Forecast Dates

Notes: This table reports coefficient estimates from the baseline SHAR model and the Tree-SHAR model for two representative one-day-ahead forecasting dates. The first panel reports estimates for January 4, 2018, a relatively normal market day, while the second panel reports estimates for March 17, 2020, during the COVID-19 stress period. The SHAR column reports the global coefficient estimates from the baseline SHAR model. The Tree-SHAR columns report leaf-specific coefficient estimates for the six terminal leaves selected by the Tree-SHAR partition. These leaves correspond to the moneyness–maturity partitions displayed in Figures 4 and 5. Newey–West \(t\)-statistics computed with 10 lags are reported in parentheses. The results illustrate how Tree-SHAR allows the HAR dynamics to vary across regions of the moneyness–maturity surface, while the baseline SHAR model imposes a single set of coefficients across all options.

January 4th, 2018

March 17th, 2020

Appendix

A SVI Method

Although we use the AHBS method to interpolate and extrapolate IVs to generate the entire volatility surface in our main analyses, we consider the SVI method of Gatheral, 2004 as an alternative approach. We specifically use its extended version suggested by Beason and Schreindorfer, 2022 to span IVs along both the moneyness and maturity dimensions.

The SVI method differs from the AHBS method in two main ways. First, the functional form of IVs is given as \[\begin{equation} \mathrm{IV}^2 (\mathbf{z}) = a + b\left(\rho(\kappa- \mu) + \sqrt{(\kappa-\mu)^2 + \sigma^2} \right),\tag{15} \end{equation}\] where \(\mathbf{z} = (m,\tau)\) represents the option characteristics (moneyness and time to maturity), \(\kappa= \log(K/F_\tau)/(VIX/100 \times \sqrt{\tau})\) is the standardized moneyness, which is a function of the strike \(K\), forward price \(F_\tau\), and \(VIX\), and the remaining five coefficients (\(a,\ b,\ \rho,\ \mu,\ \sigma)\) are linear functions of \(\tau\). For example, \(a = a_0 + a_1 \tau\). Second, the SVI method incorporates no-arbitrage constraints, which are absent in the AHBS method. Specifically, we impose the following restrictions: IVs should be nonnegative, monotonic in \(K\), and convex in \(K\); additionally, \(\tau\)-maturity model-implied volatility \(\mathrm{IV}^M (\mathbf{z}) \times \sqrt{\tau}\) should be monotonic in \(\tau\).

For implementation, on each day, we find the optimal set of parameters \(\{a_0,\ a_1,\ b_0,\ b_1,\ \rho_0,\ \rho_1,\ \mu_0,\) \(\mu_1,\ \sigma_0,\ \sigma_1\}\) by minimizing the IV errors calculated through equation (15). Then, we replace the model-implied volatility \(\mathrm{IV}^M (\mathbf{z})\) with the fitted value from equation (15) for any option characteristics \(\mathbf{z}\) of interest.

B Selecting the Number of Leaves for Tree–SHAR

We select the number of terminal leaves in Tree–SHAR by combining a validation-risk measure in implied-volatility space with the cost-complexity principle of CART Breiman et al., 1984; Breiman et al., 2017. The key object is a validation-risk curve indexed by tree size. Let \(d\) denote the number of accepted splits, so that a tree with size \(d\) has \(L(d)=d+1\) terminal leaves. For each candidate size \(d\), we compute a validation forecasting loss and then add a linear penalty in the number of splits. The selected tree size is the one that minimizes this penalized validation criterion.

Training–validation split with expanding-window evaluation. We choose the Tree–SHAR tree size based on validation forecasting performance. We take the first half of the full sample as the in-sample period and split this in-sample period chronologically at its midpoint. The first half of the in-sample period is used as the initial training subsample, while the second half of the in-sample period is used as the validation subsample. For each forecast horizon \(h\) and each candidate tree size corresponding to \(J\in\{2,\ldots,15\}\) terminal leaves, equivalently \(d\in\{1,\ldots,14\}\) accepted splits, we impose a minimum leaf size of \(10{,}000\) option observations.

For each candidate size \(d\), we conduct an expanding-window forecasting exercise within the in-sample period. Let \(\mathcal{I}_{\mathrm{val},h}\) denote the set of validation forecast origins \(t\) such that \(t\) lies in the validation subsample and the target date \(t+h\) also lies in the validation subsample. At each forecast origin \(t\in\mathcal{I}_{\mathrm{val},h}\), we estimate the Tree–SHAR partition and the leaf-specific SHAR coefficients using only observations available up to date \(t\). This expanding-window estimation sample starts with the initial training subsample and then expands through the validation period by adding all validation observations dated no later than \(t\). We then form \(h\)-step-ahead forecasts of implied volatility at date \(t+h\), conditional on information available at date \(t\). Thus, the forecast is produced at time \(t\), whi le its forecast error is evaluated ex post using the implied volatilities at the validation target date \(t+h\).

Let \(\widehat{\mathrm{IV}}^{\mathrm{Tree}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h};d)\) denote the Tree–SHAR forecast for option \(i\) at date \(t+h\), produced at forecast origin \(t\) using a tree with \(L(d)=d+1\) terminal leaves. Let \(n_{\mathrm{val},h} = \sum_{t\in\mathcal{I}_{\mathrm{val},h}} N_{t+h}\) be the total number of options in the validation samples. We summarize validation performance by the average of squared forecast errors in implied-volatility space: \[\begin{equation*} \mu^{\mathrm{Tree}}_h(d) = \frac{1}{n_{\mathrm{val},h}} \sum_{t\in\mathcal{I}_{\mathrm{val},h}} \sum_{i=1}^{N_{t+h}} \left( \widehat{\mathrm{IV}}^{\mathrm{Tree}}_{t+h\,|\,t}(\mathbf{z}_{i,t+h};d) - \mathrm{IV}(\mathbf{z}_{i,t+h}) \right)^2. \end{equation*}\] The curve \(\mu^{\mathrm{Tree}}_h(d)\) therefore measures the validation forecasting loss associated with each candidate tree size at horizon \(h\). Because the validation exercise is conducted entirely within the first half of the sample, the second-half out-of-sample evaluation period is not used when selecting the tree size.

Cost-complexity criterion. Following the cost-complexity logic of CART Breiman et al., 1984; Breiman et al., 2017, we regularize tree size by penalizing additional splits. For each forecast horizon \(h\), we define the penalized validation criterion \[\begin{equation*} K_h(d) = \mu^{\mathrm{Tree}}_h(d) + \lambda\cdot\frac{L(d)-1}{n_{\mathrm{val},h}}, \end{equation*}\] where \(\lambda\geq 0\) governs the strength of the complexity penalty. Since \(L(d)-1=d\), the penalty increases linearly with the number of accepted splits. Intuitively, an additional split is selected only if the reduction in validation forecasting loss is large enough to offset the added complexity. In practice, we set \(\lambda=0.25\) as a rule-of-thumb complexity penalty. Our results are not sensitive to small variations around this value.

We then select the optimal number of splits as \[\begin{equation*} \widehat{d}_h = \arg\min_{d\in\{1,\ldots,14\}} K_h(d), \end{equation*}\] and set the selected number of terminal leaves to \[\begin{equation*} J_h^{\ast} = L(\widehat{d}_h) = \widehat{d}_h+1. \end{equation*}\] We use this selected value \(J_h^{\ast}\) for the Tree–SHAR model at horizon \(h\), and we use the same value when constructing the bagged and boosted Tree–SHAR specifications.

Implementation and selected \(J^{\ast}_h\). Table A.1 reports the penalized criterion values \(K_h(d)\), equivalently as a function of the number of leaves \(J=L(d)\), for each candidate Tree–SHAR size and each forecast horizon \(h\) when we use the AHBS method to generate the volatility surface. We set \(J_h^{\ast}\) equal to the number of leaves associated with the smallest reported criterion value. Accordingly, Table A.1 implies \(J_1^{\ast}=6\), \(J_5^{\ast}=8\), and \(J_{22}^{\ast}=11\) for forecast horizons \(h=1\), \(5\), and \(22\), respectively. Using the SVI method to generate the volatility surface, we obtain \(J_1^{\ast}=6\), \(J_5^{\ast}=4\), and \(J_{22}^{\ast}=9\) for forecast horizons \(h=1\), \(5\), and \(22\), respectively.

Table A.1: CART \(K_h(d)\) \((10^{-4})\) with \(\lambda=0.25\) to select optimal number of leaves

C Additional Tables

Table A.2: Out-of-Sample Implied Volatility Forecasting Performance (RMSE in %) (SVI)

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models. The daily volatility surface is fitted using the SVI method as the first step of the forecasting procedure described in Section 2.1. The forecasting performance is evaluated across three different horizons (1 day, 5 days, and 22 days) and decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

\(h=1\)

\(h=5\)

\(h=22\)

Table A.3: Pairwise Diebold-Mariano Test Statistics and Model Confidence Set (MCS) for SHAR (SVI)

Notes: This table reports the pairwise Diebold-Mariano (DM) test statistics comparing the out-of-sample forecasting errors of the competing models across 1, 5, and 22-day horizons. The daily volatility surface is fitted using the SVI method as the first step of the forecasting procedure described in Section 2.1. A negative t-statistic indicates that the model in the row outperforms (has lower error than) the model in the column, while a positive value indicates that the column model outperforms the row model. Statistical significance at the 1% level corresponds to an absolute t-statistic greater than 2.33 (\(|t| > 2.33\)). The table also indicates with a checkmark the models that are included in the Model Confidence Set (MCS) at the 95% confidence level.

\(h = 1\)

\(h = 5\)

\(h = 22\)

Table A.4: Out-of-Sample Implied Volatility Forecasting Performance by Moneyness (RMSE in %, \(h=5\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of five-day horizon (\(h=5\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three moneyness groups. On each date, three groups are classified such that the moneyness is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Moneyness \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Moneyness \(<\) \(q_{80}\)

Moneyness \(\geq\) \(q_{80}\)

Table A.5: Out-of-Sample Implied Volatility Forecasting Performance by Maturity (RMSE in %, \(h=5\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of five-day horizon (\(h=5\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three maturity groups. On each date, three groups are classified such that the maturity is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Maturity \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Maturity \(<\) \(q_{80}\)

Maturity \(\geq\) \(q_{80}\)

Table A.6: Out-of-Sample Implied Volatility Forecasting Performance by Moneyness (RMSE in %, \(h=22\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of 22-day horizon (\(h=22\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three moneyness groups. On each date, three groups are classified such that the moneyness is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Moneyness \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Moneyness \(<\) \(q_{80}\)

Moneyness \(\geq\) \(q_{80}\)

Table A.7: Out-of-Sample Implied Volatility Forecasting Performance by Maturity (RMSE in %, \(h=22\))

Notes: This table reports the out-of-sample Root Mean Squared Error (RMSE) in % of 22-day horizon (\(h=22\)) implied volatility forecasts for the Random Walk, SHAR, SHAR with \(3\times3\) grids, Tree–SHAR, Bagged–SHAR, and Boosted–SHAR models by three maturity groups. On each date, three groups are classified such that the maturity is less than or equal to its 20th percentile, between the 20th and 80th percentiles, and greater than or equal to the 80th percentile. The forecasting performance is also decomposed by year to identify the specific periods where improvements occurred. Boldface values represent the lowest RMSEs among the six forecasting models.

Maturity \(\leq\) \(q_{20}\)

\(q_{20}\) \(<\) Maturity \(<\) \(q_{80}\)

Maturity \(\geq\) \(q_{80}\)

Footnotes