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.
88,151 characters · 23 sections · 91 citation commands
A GARCH model with two volatility components and two driving factors
Empirical evidence indicates that the conditional volatility of financial assets is influenced by several stochastic components. Researchers have demonstrated the superiority of employing two factors over a single factor in capturing and modeling the rich and multifaceted dynamics of volatility. For example, alizadeh02, bollerslev02, and magistris15 support stochastic volatility specifications that incorporate two sources of uncertainty, highlighting the inadequacy of a simple one-factor volatility model to fully account for the dynamic dependencies observed in the daily volatility of exchange rates. Additionally, gallant99 and chernhov03 find that two-component volatility specifications provide higher accuracy than single-factor models when fitting equity market data and conclude that at least two components are necessary to capture the dynamics of volatility.
In addition, bates96, taylor99, fouque04, schwartz09, and fouque11 demonstrate that the introduction of a second volatility factor is crucial for capturing the behavior of the volatility surface implied by options. Moreover, bates00 and christoffersen09 propose two-factor stochastic volatility models that significantly outperform closely related one-factor models both in-sample and out-of-sample when applied to equity options data. Furthermore, christoffersen12 highlight the potential of multi-factor models for capturing time and cross-sectional variations in the implied volatility structure.
The aforementioned research documents the significantly superior performance achieved by multi-factor volatility models in continuous time. However, theoretical and empirical studies have also explored the multi-factor nature of volatility in discrete time. The generalized autoregressive conditional heteroskedasticity (GARCH) model, originally proposed by bollerslev86, is one of the most popular and effective single-component models for filtering volatility in discrete time, as it inherently allows for direct estimation of volatility from historical returns. Consequently, given the complex and multifaceted nature of volatility observed empirically in financial markets, researchers have generalized the original GARCH model to incorporate two volatility components.
ding96 propose a novel approach to modeling volatility persistence by employing two GARCH components, demonstrating a significant likelihood improvement over the standard GARCH model when applied to S&P500 daily returns data. Another two-component GARCH model is developed by engle99, who provide evidence in support of volatility decomposition into two components by investigating the US and Japanese stock equity markets. adrian08 further advance this area by proposing a GARCH model with short-run and long-run components that significantly outperforms one-component GARCH-type models. Leveraging this approach and building on the work of engle99, christoffersen08 introduce a Heston-Nandi GARCH model with two volatility components, demonstrating excellent option pricing performance compared to a one-component GARCH model. The model of christoffersen08 was generalized by christoffersen14 and bormetti15 by also including the realized variance. Additional evidence supporting the effectiveness of two-component models on equity returns data is presented in christoffersen10 and conrad20.
All the two-component models mentioned so far have considered a single innovation term in the return equation. However, there is no inherent rationale to assume that both volatility components are driven by the same stochastic factor. In contrast, fouque00 and fouque03 extensively document the presence of two distinct volatility scales in returns driven by different random processes. Employing two different innovation terms seems appropriate, as the volatility components may relate to diverse sources of uncertainty. For example, adrian08, engle08, and engle13, who propose two-component models with single innovation factors, indicate that one volatility component captures market skewness risk, which may be interpreted as a measure of the tightness of financial constraints, while the other volatility component models business cycle risk.
This underscores the advantage of employing a different stochastic factor for each volatility component. Notably, zhang24 recently developed a non-affine model utilizing two variance components with two different stochastic factors by incorporating the realized variance. However, their approach implicitly assumes that both variance components are predictable, as both the returns and the realized variance are observable. Moreover, ghabani24 proposed a two-factor GARCH model that aims to capture the volatility of returns through two independent stochastic components. Nevertheless, the equations describing the two volatility components are independent, so the model does not account for spillovers between the volatility components, which may limit its ability to fully capture the complexities of the interaction mechanisms driving volatility dynamics.
Therefore, in the present work, filling a gap in the GARCH literature, we propose a GARCH model that incorporates two distinct components, each driven by an independent and unobservable innovation factor. In particular, by including two different stochastic factors, consistent with engle08, adrian08, and engle13, we allow the model to effectively adapt to diverse market conditions and discern different market trends. For instance, as evidenced by the empirical analysis conducted in this paper, we can accurately identify two volatility components with different reactions to external shocks. Furthermore, our approach provides flexible modeling of volatility spillovers, incorporating the influence of each volatility component's lagged value on the other, thereby capturing dynamic interactions between the two volatility components.
From a mathematical standpoint, the model we propose is a bivariate affine specification that extends the one-factor GARCH model of heston00, allowing for a closed-form expression of the moment-generating function, which allows for quasi closed-form derivatives pricing. Notably, our proposed framework nests the recent model by ghabani24, as it allows for a richer interaction structure between the volatility components. Additionally, the continuous-time limit of our volatility process recovers the popular stochastic volatility specification by christoffersen09. Importantly, this study is the first in the literature to derive conditions for geometric ergodicity and strict stationarity in a GARCH model driven by two stochastic factors and two volatility components with coupled dynamics.
We test the empirical performance of the proposed two-factor GARCH model by comparing it against closely related models such as the one-factor model of heston00, the component GARCH model of christoffersen08, and the two-factor GARCH model of ghabani24. We conduct extensive in-sample and out-of-sample exercises, focusing on the S&P500 total return time series. The results show that our model exhibits superior performance compared to the benchmark models in explaining the cross-section of equity returns. Specifically, the introduction of a second innovation term significantly enhances the empirical fit compared to single-innovation benchmarks. Furthermore, the advantage of incorporating volatility spillovers leads to a marked improvement in predictive accuracy, both in-sample and out-of-sample, when applied to returns data.
Finally, after deriving the risk-neutralized version of the proposed model, we utilize it for pricing options written on the S&P500 index. The empirical findings show that the inclusion of a second factor generally enhances the option pricing performance compared to single-innovation benchmark models across various moneyness and maturity levels.
The remainder of the paper is organized as follows. In Section (ref), we introduce the GARCH model with two innovations, establish sufficient conditions ensuring ensuring strict stationarity and geometric ergodicity, and derive a continuous-time limit. Section (ref) discusses the parameter estimation procedure in detail. In Section (ref), we briefly review some popular (affine) GARCH models that we use as benchmarks. Section (ref) presents the empirical performance of the proposed model and the benchmarks when applied to the S&P500 total return time series. In Section (ref), we perform the risk-neutralization, derive the expression for the characteristic function, and in Section (ref), we test the models in pricing options written on the S&P500. Finally, Section (ref) concludes. All mathematical proofs are gathered in the appendix.
In this section, we introduce a novel bivariate GARCH model with two innovation factors to capture the different components of the volatility of log-returns. This model will be labeled {GARCH with two factors} (GARCH-2F).
Let $S_t$ denote the price of a risky asset at time $t \in \mathbb{Z}$ and let us consider the total log-return $R_t = \ln (\frac{S_t + D_t}{S_{t-1}})$, including the dividend $D_t$. We model the return process as follows:
where $v_{1,t}$ and $v_{2,t}$ denote the time-varying variance components. Moreover, we specify the conditional expected return as $\mu_t = r + \lambda (v_{1,t} + v_{2,t})$, where $r$ represents the (constant) risk free rate and $\lambda$ is the (constant) risk premia parameter. Unlike traditional GARCH models and similarly to continuous time stochastic volatility models, the specification in (ref) contains two sources of risk, $Z_{1,t}$ and $Z_{2,t}$, and we assume $Z_{1,t} \overset{\mathrm{i.i.d.}}{\sim} N(0,1)$, $Z_{2,t} \overset{\mathrm{i.i.d.}}{\sim} N(0,1)$ with $Z_{1,t} \perp Z_{2,t}$ for all $t \in \mathbb{Z}$. Note that the (conditional) total variance of returns is equal to $v_{1,t} + v_{2,t}$.
We model the time-varying variance components via the following bivariate system of equations:
Equation (ref) specifies two volatility processes governed by Heston-Nandi type GARCH(1,1) dynamics, see heston00, allowing spillovers between each variance component and innovation component. The specification includes only one lag in order to keep a parsimonious structure, but further lags can be added. Note that there is no spillover if $\beta_{12} = \beta_{21} = \alpha_{12} = \alpha_{21} = 0$. Moreover, if $\omega_2 = \beta_{12} = \beta_{22} = \beta_{21} = \alpha_{12} = \alpha_{22} = \alpha_{21} = 0$ we obtain a model with a single volatility component driven by the standard Heston-Nandi dynamic, see heston00.
In this section, we derive sufficient conditions under which the model given by (ref) and (ref) generates stationary and ergodic trajectories.
We prove the geometric ergodicity of $ \{R_t, \boldsymbol{v}_t\}$, where $\boldsymbol{v}_t = (v_{1,t}, v_{2,t})^T$, based on the Markov chain stability theory of nummelin84 and tweedie93. To this aim, let us define $\mathcal{D} = [\omega_1,\infty) \times [\omega_2,\infty)$ and let us consider the following (non-linear) state representation of model (ref)-(ref):
where $\boldsymbol{W}_{t} = \boldsymbol{Z}_{t-1}$ and ${G}: \mathcal{D} \times \mathbb{R}^2 \rightarrow \mathbb{R}$, and $\boldsymbol{F}: \mathcal{D} \times \mathbb{R}^2 \rightarrow \mathcal{D}$, are defined as follows:
where $\boldsymbol{\omega} =
^T$, $\boldsymbol{\beta} =
$ and $ \boldsymbol{\alpha} =
$.
This process is a homogeneous Markov chain with state space $(\mathcal{D}, \mathcal{B})$, where $\mathcal{B}$ is the Borel $\sigma$-algebra on $\mathcal{D}$. Following Chapter 7 of tweedie93, we define inductively a sequence of functions $\boldsymbol{F}_t$ for $t=1,2,3,\ldots$, by $\boldsymbol{F}_{t+1}(\boldsymbol{x},\boldsymbol{W}_1,\ldots,\boldsymbol{W}_{t+1}) = \boldsymbol{F}(\boldsymbol{F}_t(\boldsymbol{x},\boldsymbol{W}_1,\ldots,\boldsymbol{W}_{t}),\boldsymbol{W}_{t+1})$, so that, for any initial condition $\boldsymbol{v}_0 = \boldsymbol{x}$ and $t \geq 1$, we can use equation (ref) recursively to obtain $\boldsymbol{v}_t = \boldsymbol{F}_{t}(\boldsymbol{x},\boldsymbol{W}_1,...,\boldsymbol{W}_{t})$. If we replace the random disturbances $\{ \boldsymbol{W}_t\}$ with a deterministic control sequence, say $\{ \boldsymbol{w}_t\}$, we obtain the so-called deterministic control model associated to the non-linear state space model (ref), see Chapter 7 of tweedie93 for more details.
We use $\mu_{\text{Leb}}(\cdot)$ to denote the Lebesgue measure and $P^n\left(\boldsymbol{x}, A \right) = P \left( \boldsymbol{v}_n \in A \middle| \boldsymbol{v}_0 = \boldsymbol{x} \right)$, for $\boldsymbol{x} \in \mathcal{D}$ and $A \in \mathcal{B}$, to denote the $n$–step transition probability measure of the Markov chain $\boldsymbol{v}_t$. For $n=0$ we have $P^0(\boldsymbol{x},A) = \mathrm{1}_A(\boldsymbol{x})$, that is the indicator of the set A: $1_A(\boldsymbol{x})=1$ if $\boldsymbol{x} \in A$, $1_A(\boldsymbol{x})=0$ if $\boldsymbol{x} \not\in A$. For $n = 1$ the notation $P (\boldsymbol{x}, A)$ is used. We denote with $\left\lVert\cdot\right\rVert$ the $L_2$-norm of any vector and matrix and by $\rho(\boldsymbol{A})$ the spectral radius of any square matrix $\boldsymbol{A}$, i.e., $\rho(\boldsymbol{A}) = \max \{ |\varphi_i|: \varphi_i$ is an eigenvalue of $\boldsymbol{A}$\}.
For later purposes, we introduce the autoregressive matrix of the model (ref):
We state the definitions of the main concepts related to stationarity and ergodicity.
If $\{ \boldsymbol{v}_t \}$ is $\psi$-irreducible there exists a maximal irreducibility measure $M$ on $(\mathcal{D}, \mathcal{B})$, i.e., an irreducibility measure such that all other irreducibility measures are absolutely continuous with respect to $M$ (tweedie93). Moreover, we set $\mathcal{B}^+ = \{ A \in \mathcal{B} : M(A) > 0 \}$.
The following lemma appropriately defines the measure $\psi$ so that we achieve $\psi$-irreducibility and aperiodicity under mild-assumptions.
The proof of the Lemma, together with the other proofs, is reported in the appendix of the article. Here we note that in the proof of Lemma (ref) we also derive an explicit form for the density of the transition probability associated to (ref), see equation (ref).
Studying the small sets of an irreducible chain allows us to derive its long-run probabilistic behavior. Specifically, to prove geometric ergodicity, we first establish the following proposition and preliminary lemma.
Then, given that $ \{ \boldsymbol{v}_t \} $ is irreducible and aperiodic, an appropriate small set exists and we can provide sufficient conditions for the geometric ergodicity of the process $\{R_t, \boldsymbol{v}_t\}$.
We note that a sufficient ergodicity condition analogous to $\varphi < 1$ is obtained also in hafner09.
In this section, we derive the continuous time limit of the model in equations (ref)-(ref). We adopt the convergence scheme of nelson90 and also use the same notation.
We note that the set of equations (ref) is analogous to the bivariate square-root stochastic volatility specification considered, for example, in christoffersen09.
In this section, we present three parsimonious GARCH specifications nested within the GARCH-2F model. This enables a comparison between the full GARCH-2F model and its nested versions, where spillovers are either restricted or entirely absent.
In this model we partially remove the spillover between the two volatilities in Equation (ref), by imposing:
In this model we remove the spillover between each volatility and the innovation driving the other. This amounts to imposing:
In this model we remove all the spillover components given by the $\alpha$ and $\beta$ coefficients. This amounts to setting:
By implementing these restrictions, we can completely eliminate the spillovers between the two volatility components, resulting in a more parsimonious variance process. Notably, under the specification in (ref), we recover the interesting two-factor GARCH model introduced by ghabani24.
In this section, we illustrate a feasible and computationally efficient approach to estimate model (ref)-(ref). Equation (ref) contains two sources of unobservable uncertainty, $Z_{1,t}$ and $Z_{2,t}$, which implies that we cannot estimate the model as done in traditional GARCH models, rather we need a filtering method. Following christoffersen12, we note that a key property of a filter is that the filtered states are equal to their expected values conditional on the relevant information set. Hence, we define the filtered estimates of ${Z}_{1,t}$ and ${Z}_{2,t}$ conditional to the filtration $\mathcal{F}_t = \sigma(\{R_t, R_{t-1}, \ldots \})$ at time $t$ as
Since $Z_{i,t}$ is normally distributed, then also $R_t$ is normally distributed so that the conditional expectation in equation (ref) can be computed using a well-known result on Normal distributions conditioning, see Chapter 2 of bda16, that is:
where $\widetilde{\mu}_t = r + \lambda (\widetilde{v}_{1,t} + \widetilde{v}_{2,t})$ and $\widetilde{v}_{1,t}$ and $\widetilde{v}_{2,t}$ denoting the filtered estimates of the conditional variances ${v}_{1,t}$ and ${v}_{2,t}$ respectively, which can be computed using equation (ref)
The log-likelihood function can now be constructed as the product of the conditional distributions across the sample. Specifically, conditional on $\widetilde{v}_{1,0}$ and $\widetilde{v}_{2,0}$, the log-likelihood function of returns is given by
where $T$ denote the length of the daily log-returns series.
Given the affine and two-component structure of model (ref)-(ref), we will compare it with the following models: the GARCH model developed by heston00, hereafter GARCH-HN, the two-component GARCH model introduced by christoffersen08, hereafter GARCH-CJOW model, and the GARCH-2F$\alpha\beta$ of ghabani24. Both the GARCH-HN and GARCH-CJOW models employ a single innovation factor, and, for the reader's convenience, they are briefly recalled below.
The return and the volatility processes are defined as follows
where $Z_t \overset{\mathrm{iid}}{\sim} N(0,1)$.
The component GARCH model proposed by christoffersen08 comprises three equations, one for the return process, one for the long-term variance component $q_{t}$, and one for the short-term variance component $s_{t}$:
where $Z_t \overset{\mathrm{iid}}{\sim} N(0,1)$.
In this section, we present empirical results for the total S&P500 daily adjusted log-returns time series. Our data span from January 5, 1988 to December 29, 2023 (9069 daily observations), as shown in Figure (ref).
The estimation procedure described in Section (ref) requires, besides the time series of log-returns, also the risk-free interest rate, which we proxy using the 3-month US Treasury Bill rate. The data for the daily levels of the Adjusted Close Price for the S&P500 total Returns series is retrieved from Refinitiv Datastream, whereas the data for the 3-month Treasury Bill rate are gathered from the Federal Funds Effective Rate (FRED) dataset.
The in-sample results, based on the maximum likelihood estimation described in Section (ref), are illustrated in Table (ref). To assess the model performance we consider the Akaike (AIC) and the Bayesian information criteria (BIC). To further assess the goodness-of-fit of our model specification we performed a likelihood ratio (LR) test for the $\textit{GARCH-2F}$ model and its three nested models versus the one factor GARCH-HN model. The LR statistics are reported at the bottom of Table (ref).
As we may see, the $\textit{GARCH-2F}$ model and its nested versions provide superior goodness-of-fit to the returns data compared to the one-factor benchmark models. The LR statistic is significantly higher, and the AIC and BIC values are lower by at least 60 points. Moreover, the most parsimonious model, the $\textit{GARCH-2F}\alpha \beta$ proposed by ghabani24, achieves 30 likelihood extra points compared to $\textit{GARCH-CJOW}$ and 130 points compared to $\textit{GARCH-HN}$. However, in comparison to $\textit{GARCH-2F}\alpha \beta$, results also highlight the significance of including volatility spillovers, either through $\alpha$ for the GARCH-2F$\beta$ model or the $\beta$ coefficient for the GARCH-2F$\alpha$ model specification.
We further assessed the volatility persistence for all the considered models. The persistences of the two variance components $v_{1,t}$ and $v_{2,t}$ in model (ref)-(ref) are defined as the two eigenvalues of the autoregressive matrix $\boldsymbol{B}$ in (ref). Whereas, the persistences, of the GARCH-CJOW model, as discussed in christoffersen08, are equal to $\beta_{11} + \alpha_{11} \gamma^2_1$ for the one volatility component and $\beta_{22} + \alpha_{22} \gamma^2_2$ for the other volatility component. At the bottom of Table (ref) we present the persistences of the volatilities for all the models. A feature that characterize all the component models is that the component displaying the highest $\beta$ coefficient is also the one possessing the lowest reaction to innovation shocks, as indicated by the $\alpha$ coefficient. However, the persistence values do not differ significantly among the nested model specifications and are aligned to the peristences of the GARCH-CJOW model.
When analyzing the model parameters, we observe that for all considered models, the first volatility component returns a $\beta_{11}$ estimate around 0.9, which is commonly encountered in the GARCH literature (see heston00, christoffersen08). A similar pattern holds for $\alpha_{11}$ and $\gamma_1$. The parameter $\omega_1$ is largely non-significant across most models, except for the GARCH-CJOW. In the component models, the second volatility component exhibits higher sensitivity to innovations, as reflected by $\alpha_{22}$ and $\gamma_2$. The autoregressive parameter $\beta_{22}$ for the second component is generally lower than $\beta_{11}$ and is not always significant in the GARCH-2F, $\textit{GARCH-2F}\beta$ and $\textit{GARCH-2F}\alpha$ models.
If we consider the spillover effects between the two volatility components in the $\textit{GARCH-2F}$ and the $\textit{GARCH-2F}\beta$ models, the coefficients $\alpha_{12}$ and $\alpha_{21}$ are highly significant. In the $\textit{GARCH-2F}\alpha$ model, the spillovers through the $\beta_{12}$ and $\beta_{21}$ coefficients are also significant. Consequently, the more parsimonious GARCH-2F$\alpha \beta$ specification proposed by ghabani24, which eliminates spillovers between the two volatility components, is not fully capable of capturing the dynamics of return data.
In Figure (ref), we plot the filtered conditional volatility of the $\textit{GARCH-HN}$ model and the two volatility components of the $\textit{GARCH-CJOW}$ model and $\textit{GARCH-2F}$ model. One of the volatility components of the $\textit{GARCH-2F}$ model mimics the dynamics of the one-factor $\textit{GARCH-HN}$, and the other volatility component is a slower moving process with much lower variation over the considered time span.
In Figure (ref), we compare the conditional variances of the one-factor $\textit{GARCH-HN}$ model and of the two-factor $\textit{GARCH-2F}$ model, that is $v_{1,t} + v_{2,t}$. The variance of both models displays a similar pattern, however during highly volatile periods, as in 2008 and 2020, the higher spikes of the $\textit{GARCH-2F}$ might yield a better descprition of the high volatility of such periods.
To test the forecasting performances of the rival models, we forecast the Value-at-Risk ($\mathrm{VaR}$) at the 1% and 5% significance levels from one to five days ahead. We test the models during different periods of crisis by splitting the dataset into two subperiods (4534 observations each) and running separate out-of-sample analyses: from January 5, 1988 to December 21, 2005 (to include the September 11, 2001 crisis) and from December 22, 2005 to December 29, 2023 (to include the 2008 financial crisis and the COVID-19 period). For each of the considered subsamples, we estimate the model parameters using a rolling window equal to 60% of the days in the period. That is, we start by using the first 60% of the daily S&P500 total returns to calibrate the models, leaving the remaining data for the out-of-sample analysis. Then, for each day in the out-of-sample, we re-estimate the model parameters using the past 60% daily observations and compute forecasts from one to five days ahead.
In particular, for all the competing models, at each forecasting time, we simulate $100,000$ sample paths of $R_{t+\ell}$, with $\ell = 1,2,3,4,5$, according to the return equation (ref). This is done by randomly simulating the $Z_{1,t}$ and $Z_{2,t}$ standard normal random variables. Then, we estimate the $\mathrm{VaR}$ at the 1% and 5% significance levels by computing the empirical quantiles of the simulated sample paths of $R_{t+\ell}$.
We then measure the statistical accuracy of the VaR forecasts by performing the proportion of failures test of kupiec95 and the correct conditional coverage test of christoffersen98. The results obtained for the proposed $\textit{GARCH-2F}$, $\textit{GARCH-2F}\alpha$ and $\textit{GARCH-2F}\beta$, together with the competing $\textit{GARCH-2F}\alpha \beta$, $\textit{GARCH-HN}$ and $\textit{GARCH-CJOW}$ models, are reported in Tables (ref) and (ref).
We note that, overall, the novel two-factor specifications result in a significantly lower number of failures compared to the one-factor models. For the period from January 1988 to December 2005, both the Kupiec and Christoffersen tests confirm the superior performance of the GARCH-2F models, particularly the GARCH-2F$\beta$, which consistently performs well at both the 95% and 99% confidence levels, outperforming the GARCH-2F$\alpha \beta$ model with no spillovers between the volatility components. For the period from December 2005 to December 2023, the GARCH-2F and GARCH-2F$\beta$ models exhibit similar performances.
In addition, we further compare the accuracy of the VaR forecast using the Model Confidence Set (MCS) introduced by hansen11. Following gonzales04, we measure performance by means of the quantile loss function as in koenker78, defined as
where $\mathrm{VaR}_{t+\ell}(\vartheta)$ denotes the $\mathrm{VaR}$ on day $t+\ell$ at the significance level $\vartheta$. Moreover, to check if the forecasts of the models $\textit{GARCH-CJOW}$, $\textit{GARCH-2F}$, $\textit{GARCH-2F}\alpha$, $\textit{GARCH-2F}\beta$ and $\textit{GARCH-2F}\alpha \beta$ significantly differ from the ones of the $\textit{GARCH-HN}$ benchmark we performed the test of diebold02. The quantile losses and the p-values are reported in Table (ref).
Overall, the $\textit{GARCH-2F}$ model, its three nested specifications, and the $\textit{GARCH-CJOW}$ model return significantly different forecasts compared to the ones of the $\textit{GARCH-HN}$ as indicated by the p-values for both periods and significance levels. For the period from January 5, 1988 to December 21, 2005, all the four $\textit{GARCH-2F}$ models return lower quantile lossess, thus indicating a superior VaR forecasting accuracy and are included in the MCS, contrary to the $\textit{GARCH-CJOW}$ model. In the period from December 22, 2005 to December 29, 2023, only the $\textit{GARCH-2F}\beta$ is always included in the MCS and outperform the benchmark models by returning lower quantile losses.
In this section, we derive results that allow us to value option contracts using the four proposed two-factor models. First, we perform the risk-neutralization of these models. Second, we estimate the GARCH-2F model (and its nested models) by jointly using returns and options data.
Equations (ref)-(ref) define the return and the volatility dynamics under the physical measure. In this section, we risk-neutralize the model, so that we can compute option prices. Following christoffersen12, we consider a pricing kernel with affine dynamic
where $\mathbb{E}$ denotes the expectation operator under the physical measure. Accordingly, we guess the conditional Radon–Nikodym derivative in the following form
and we impose
so that using equations (ref)-(ref) we obtain $\Lambda = - \xi$. Therefore, we have the following result.
Finally, the following proposition yields the dynamic of the risk-neutralized processes $ \{R_t\}$, $ \{v_{1,t}\}$ and $ \{v_{2,t}\}$.
Equations (ref)-(ref) allow us to derive the moment generating function using a standard procedure as in heston00.
Let us consider a European Call option on the underlying asset $S_t$ with strike price $K$ and maturity $T$. Using the inversion theorem in gilpelaez51, the option price at a generic time $t < T$, which we denote by $C_t$, is computed as follows:
where, as in heston00, $f^*_t$ denotes the (conditional) moment generating function under the risk-neutral measure, that is function (ref) evaluated using risk-neutral parameters. Put option values can be calculated using the well-known put-call parity as in heston00.
We consider European options, both Put and Call, written on the S&P500 index, with data retrieved from Thomson Reuters Eikon Datastream. Specifically, we consider options with maturities ranging from 2021 to 2023 and the time series of their daily prices from February 10, 2021 to December 29, 2023.
As a common practice, see for example christoffersen12, enzo23 and enzo24, we apply several exclusion filters to obtain the final panel of option contracts. We keep only the options with time-to-maturity between 14 and 365 days and we select only out-of-the-money Put and Call options (we compute the moneyness as $K/S_t$, where $K$ is the strike price and $S_t$ is the underlying index level), and we filter out illiquid quotes by selecting only the six most liquid strikes at each maturity, and we consider option quotes only on Wednesday. Finally, we remove price quotes lower than $3.8\$$.
In Table (ref), we report the resulting number of option prices on the S&P500 index, sorted by moneyness and days to maturity, for a total of $N=9372$ options prices.
A common practice for estimating parameters using daily time series of log-returns and a panel of option contracts is using a suitably weighted log-likelihood function, see christoffersen12 and enzo23. The log-likelihood for the returns series, $\ell_{\text{returns}}(\boldsymbol{\theta})$, is already available in equation (ref), whereas, for the $N$ option prices, we proceed as follows. First, we define the error on the volatility implied by the $i$-th option price as
where $\textit{IV}^{MKT}_i$ and $\textit{IV}^{MOD}_i$ denote the market and the model volatilities implied by the $i$-th option price, which we compute via the Black-Scholes model. Then, by assuming a Gaussian distribution for $e_i$, we can derive the log-likelihood associated to the option pricing errors as
where $\sigma^2_e$ represents the variance of the option pricing errors.
The number $N$ of data points available from the option panel could be significantly higher than the number $T$ of daily returns, see Section 4. Therefore, following christoffersen12, we assign an equal weight to returns and option prices by considering the following weighted joint (total) log-likelihood:
The joint estimation results are presented in Table (ref). We calculated standard information metrics AIC and BIC, and we also report the persistence values, as we did in Table (ref). We observe different persistence values around 0.98 and 0.97 for the two volatility components, compared to the returns data-only estimation, which were 0.99 and 0.92. Additionally, we report the variance of the option pricing errors. The proposed GARCH-2F$\beta$ model return the highest joint likelihood values and lower AIC and BIC compared to the benchmark models, indicating improved performance in fitting both option and returns data. The GARCH-2F models allowing for spillovers have highly significant parameters and provide a superior fit compared to the GARCH-2F$\alpha \beta$ model.
To assess the performance of the six competing models (including the models nested by the GARCH-2F) in option pricing, we follow christoffersen12 and employ the (percentage) implied volatility root mean square error:
The results obtained, shown in Table (ref), reveal that the novel two-factor specifications yield better performance compared to the benchmark approaches when pricing options. In particular, the GARCH-2F model and its nested models return lower (from 2% to 4%) overall implied volatility pricing errors compared to the one factor GARCH-HN benchmark model. By sorting the options for moneyness and days to maturity, the GARCH-2F and its nested models provide consistent results, that is lower IVRMSE than the one-factor benchmark models, and reveal more robust performance in pricing out-of-the money options ($0.8 \leq K/S_t \leq 0.9$). Notably, the parsimonious GARCH-2F$\alpha\beta$ model demonstrates slightly better performance across both moneyness and days to maturity.
Accurately capturing the the complex and multi-faceted nature of volatility is crucial for effective financial analysis and risk management. In this paper, we introduce a novel GARCH specification with two volatility components each driven by an independent stochastic factor, an approach commonly used in the continuous time (stochastic) volatility models, see, e.g., christoffersen08 and fouque11, yet still overlooked in discrete-time GARCH models.
We conducted a theoretical investigation that establishes, for the first time, sufficient conditions for strict stationarity and geometric ergodicity in a two-component GARCH model driven by independent stochastic factors. Additionally, we derived a continuous-time limit for the proposed GARCH-2F model, which recovers a well-known bivariate square-root stochastic volatility specification. Finally, we considered three more parsimonious nested models: the GARCH-2F$\alpha$, GARCH-2F$\beta$, and GARCH-2F$\alpha \beta$, the latter introduced by ghabani24.
We empirically tested the two-factor models on S&P500 total log-returns data and compare their in-sample and out-sample performance against one-factor affine benchmark models. Overall, the two-factor models demonstrate superior in-sample and out-sample $\mathrm{VaR}$ prediction accuracy. Then, we also considered the pricing of options on the S&P500 index. The two-factors model provide lower implied volatility pricing errors, particularly for out-of-the money options, when estimated jointly using returns and options data.
The approach proposed in this article can be extended in many directions. Potential extensions include include incorporating fat-tailed distributions for the innovations while preserving model affinity, exploring different specifications for the drift of the return process (e.g., an autoregressive model), and considering additional lags.