1. Introduction and Motivation
Asset prices are typically modeled with the Geometric Brownian motion (GBM) of the form
(1)
where St is the asset price, μ is the drift, σ is the volatility, and dBt is a standard Brownian motion.
Rapid developments in vanilla and exotic options markets, over the past several decades, have fundamentally challenged the static assumptions of μ and σ parameters in GBM. Merton [1] introduces jumps and shows that if the logarithm of the percentage jump is normally distributed, a closed form solution for European style options exists. Cox and Ross [2] create the constant elasticity of variance (CEV) model, where an exponential parameter α is added to the asset price. The value of α determines the dependence between asset price and volatility. In a pure jump extension, Madan et al. [3] follow a variance-gamma approach, which creates heavier tails and provides semi-analytic expressions for European style options.
Heston [4] in a seminal model correlates the asset return process with stochastic variance. Many Heston model extensions have since developed with Zhou [5], Hagan et al. [6], Brigo and Pallacinini [7], and Langnau [8] counted among many prominent examples.
The availability of market data of implied volatility surface1 has ushered in a new era of stochastic volatility modeling to reproduce a multitude of finer properties observed in the real world. Cont et al. [9] prove that 1-factor model is insufficient to represent the true dynamics seen in SPX volatility surfaces. Furthermore, the eigenfactors of the vol-surface time series are found not to be perfectly correlated to the underlying asset price movements—hence concluding—“Vega” risk cannot be reduced to “Delta” risk.
In a series of ground-breaking papers, Bergomi [10] and others set out to capture the term-structure of skew and vol-of-vol in a Forward Variance Swap framework, an analog of the Libor Market Model. A class of models has since multiplied aimed to price the “smiling” volatility derivatives consistently to vanilla options, therefore solving the “Joint S&P/VIX Smile Calibration Puzzle”. By adding simultaneous jumps to the Ornstein-Uhlenbeck process of the forward variance swap and the underlier GBM process, the flexible Lévy specification by Cont et al. [11] offers a greater analytical tractability and more efficient control of the vol-skew and the shifting correlation between spot and implied volatility.
Joint modeling of asset volatility and asset correlation should attract more research interest, as the need arises not only for asset allocation but also for derivatives pricing and risk management. By comparison to stochastic volatility, less progress has been made in the literature regarding asset correlation as a state-dependent market factor.
The summary statistics of realized stock correlations in Table 1 reveals that the level and variability of cross-sectional asset correlation both move higher conspicuously when market is in distress, from a distinctly lower regime observed during expansionary economic times. Meissner [12] further asserts that elevated correlations persist without necessarily accompanied by rising volatility during recessionary periods.
Table 1. Dow stock correlation level and Std. Dev. from Jan. 1972 to Mar. 20202.
Analogous to the irreducible “Vega” risk in the Implied Volatility (IV), the empirical evidence clearly supports the role of the Implied Correlation (IC) as a unique source of randomness to the financial system.
Hull et al. [13] model the asset process in a GBM setting and then sample the asset correlation from a beta distribution, with the distribution parameters exogenously derived and not time-varying. Ma [14] models stochastic volatility and stochastic correlation processes to price exchange rate options, and Ma [15] applies stochastic correlation to price and hedge multi-asset options, the volatility and correlation are however assumed independent from each other.
Emmerich [16] highlights the unique mathematical properties of a stochastic correlation process. Dϋllman et al. [17] model stochastic correlation with a Vasicek process. However, both papers do not combine the correlation process with the asset process.
Buraschi et al. [18] and Fonseca et al. [19] extend the Heston [4] model, by correlating a n-dimensional stochastic correlation process with a n-dimensional stochastic asset process.
The term of correlation is ubiquitous. It should be emphasized that the asset correlation in our study is restricted to the equity price returns. In contrast to Burtschell et al. [20], the local stochastic correlation model targets the default time correlation in a portfolio of Credit Default Swaps (CDS) instead.
The remaining paper is structured as follows: In Section 2 we introduce our unified stochastic asset volatility—stochastic asset correlation model. In Section 3 we show the real world fit of the model. Section 4 discusses the calibration of the eight parameters. Section 5 conducts the significance tests. Section 6 applies the model in an enhanced version of GBM.
2. The USVSC Model
Our proposed Unified Stochastic Volatility—Stochastic Correlation model (USVSC henceforth) consists of three stochastic differential equations:
where stochastic instantaneous volatility
and stochastic instantaneous correlation
are set within a local volatility framework, mean-reverting with a rate of
to a long-term mean of
. The diffusion coefficient is denoted as
, and the positive skew in
is captured by a power parameter
. The modeled stochasticity is generated by three 1-dimensional standard Brownian motions
, where
and
orthogonal,
and
correlated with a coefficient
.
Equation (2), referred to as CIRCEV hereafter, extends the Cox-Ingersoll-Ross (CIR) model with a β exponent parameter, otherwise known as Chan-Karolyi-Longstaff-Sanders (CKLS) model proposed by Chan et al. [21], or a variant of the Constant Elasticity of Variance (CEV) model. Evoking CIRCEV model specification is motivated by its non-negativity and mean-reversion properties fitting the implied volatility, and the presence of a positive skew found empirically and prominently displayed in Figure 1(c).
Equation (3) specifies a Jacobi process, otherwise known as Wright-Fisher diffusion process. Jacobi is a natural model choice for implied correlation, in which not only the mean-reverting empirical feature is preserved, the limiting behavior for a correlation is also obeyed naturally at −1 or +1. As discussed by Ahdida et al. [22], Emmerich [16] and Zetocha [23], Jacobi specification permits analytical solutions of the first two moments, as well as the average over an arbitrary length of time. Although, like Heston [4], a Feller-like inequality constraint should apply over model parameters to ensure the desired numerical stability.
Equation (4) gives a dependence structure between the implied volatility and the implied correlation, analogous to Heston [4] in treating the joint distribution between spot and local variance.
In summary, our USVSC model is formulated as a correlated CIRCEV and Jacobi processes in an attempt to capture the equity derivatives dynamics. The rational for disallowing the joint-dynamics between derivatives and underlying spot market is two-fold: 1) to keep the model at the lowest possible dimension in support of historical calibration; and 2) to expose the risk factors exclusively within derivatives, as opposed to a joint model traditionally done for efficient hedging purpose.
Our model consists of a total of eight parameters, namely,
, forming the set of unknowns to be estimated empirically, with a goal to better understand the historically realized marginal and joint distributions of the risk neural option implied volatility and implied correlation through time.
3. Real World Fit
To measure the efficacy of our model, we let VIX represent IV (the option-implied volatility
), and the CBOE Implied Correlation as IC (the option-implied correlation
)3.
Figure 1. (a) Empirical relationship between IV and IC. (b) Time series plot of the empirical IV and IC. (c) Density of empirical IV. (d) Density of empirical IC.
Our dataset is comprised of the SPX, the VIX and the ICJ/JCJ indexes observed from 01/03/2007 to 08/01/2019, daily sampled and sourced from http://www.cboe.com/. CBOE publishes both historical data series of the Implied Volatility (known as the VIX—the “fear gauge”), and the Implied Correlation, derived from standard 1-month options and dispersion trading strategies respectively, under the methodologies given by CBOE [24] and CBOE [25].
We first display the excellent real world fit of our model, which is followed by a detailed description of a 2-staged calibration and the correspondent statistical test results.
Figure 1 describes historical properties of the Implied Volatility (IV) and the Implied Correlation (IC).
Figure 2 displays the simulated properties of the Implied Volatility (IV) and the Implied Correlation (IC) using our calibrated model results.
With respect to Figure 2(a), the USVSC model is designed to produce a positive-“triangular” relationship between IV and IC as shown in Figure 1(a). A positivity in
is expected to enable synchronized co-movements between
and
. Furthermore, as the correlation approaches upper or lower bounds, i.e. +1 or −1, its stochasticity decreases due to the diffusion term
specified in Equation (3).
After outlining our calibration procedure in Section 3, we perform a series of statistical tests associated with Figures 2(a)-(d) in Section 4.
4. Model Calibration
The estimation scheme is described in this section. We approach the calibration in three Stages using the Maximum Likelihood Estimation (MLE) framework. We believe MLE is simplest to implement, numerically stable while allowing for the generation of the desired statistical inferences.
Figure 2. (a) Modeled relationship between IV and IC. (b) Time series plot of the modeled IV and IC. (c) Density of modeled IV. (d) Density of modeled IC.
4.1. Stage I—Marginal Distributions of IV and IC
By design of the proposed USVSC model, except for the correlated Brownian motions anchored by Equation (4), IV and IC have no interactions in the drift as well as the diffusion terms, divergent from that in Heston [4]. Under this specification, therefore, a decoupled and 1-dimensional historical calibration, separated for IV or IC process, is permissible.
To illustrate, we first suppose a generalized univariate stochastic differential equation (SDE) with drift and diffusion terms fully parameterized as follows:
(5)
where
represents the time series of empirical observations, and
the collection of model parameters. We assume the increments of observations
conditionally Gaussian
. This means, based on the Euler discretization scheme, we can derive explicitly the first two moments of the increments, similar to the approach taken by Ait-Sahalia [26]:
(6)
The transition density, under our assumption, can therefore be expressed in a closed-form:
(7)
The optimization objective function immediately follows, resulting in a maximization problem of the log-likelihood function to solve
:
(8)
With an unconstrained optimization routine4, we first report in Table 2 the results of USVSC model parameter estimation of
, along with the corresponding estimation errors.
The calibrated parameters of IV, per Equation (2)—CIRCEV specification, including the long-term mean, the reversion rate (half-life), the magnitude in vol-of-vol, together with the CEV skew parameter
, are all of expected values with a good match to existing literature, such as Ait-Sahalia and Kimmel [27], and Bu et al. [28].
The calibrated parameters of IC, per Equation (3) Jacobi process, are shown to have met the inequality condition of the mean-reverting rate,
, consistent to relevant discussions in Emmerich [16], Ma [14] and [15], as well as Zetocha [23].
To test the robustness, we produce 5k simulations to visualize in Figure 3. We confirm the fulfillment of correlation boundedness achieved by the calibrated parameters. In addition, the mean is fairly stable around the randomly assigned initial value of 0.2.
The red line in Figure 3 represents the time series mean of all simulated IC paths, which exhibits a relative stability as IC converges to a long-term mean prescribed by a Jacobi process as in Equation (3). Measured by each simulated single path shown in Figure 2(b), and relative to IV, however, a greater magnitude of variability in IC empirically observed in Figure 1(b) is largely captured, helped in part by a faster mean-reverting rate
.
Figure 3. Time series of simulated implied correlation.
Table 2. Univariate calibration of IV and IC.
4.2. Stage II—Correlated Brownian Motions
Having estimated the seven parameter values in
,
, required by the model, we now turn our focus to the eighth parameter
and the correlated Brownian motions in Equation (4).
We reconstruct the realized increments of Brownian motions,
and
, making use of the estimated parameters from Stage I. Under the normality assumption, it follows that the asymptotic correlation coefficient estimation as specified in Equation (5) is just the expectation of the product of
and
scaled by dt:
(9)
The estimated correlation coefficient parameter
, together with the t-statistic and p-value are presented in Table 3. The zero-value null hypothesis is strongly rejected at 95+% confidence level. This is an important finding from our Stage II calibration.
The histogram of the product of the white noises is displayed in Figure 4 for illustration purpose.
5. Goodness-of-Fit Tests
In this section, we test if our calibrated model is statistically accurate in representing the empirical data. The marginal distributions of IV and IC resulted from calibration Stage I are first examined, which are followed by normality tests concerning the correlation parameter
estimated in calibration Stage II.
Figure 4. Histogram of
.
Table 3. Correlation calibration result of pw of Equation (4).
.
5.1. Testing for Stage I
To test the equality in aggregation distribution between the historical IV and IC densities in Figure 1(c) and Figure 1(d), and our model-simulated densities in Figure 2(c) and Figure 2(d), we proceed as follows.
Firstly, in Table 4, we perform non-parametric Kolmogorov-Smirnov (KS) test for the difference in univariate distributions. With the cumulative distribution functions (CDF) denoted as F(IV) and F(IC), respectively, between the empirical distributions and the modeled 1-path time series simulation, denoted as
and
, started with a pair of arbitrary points (IV0, IC0). Below 5% are the max distances between two CDFs measured by the D-Scores associated with IV and IC. Nearing 1% p-values found by the tests, however, suggest the null hypothesis should be rejected with 95% confidence level.
Secondly, in Table 5, we utilize a form of randomization test—the Permutation Test of Equality. With paired samples of univariate distributions, the high p-values registered for both IV and IC suggest a high degree of equality in distribution between the model generated IV or IC process and the correspondent empirical counterpart.
Graphical depiction rendered in Figure 5 shows the density functions of IV and IC, respectively, between the simulated series and the empirical observations. The blue color identifies the largest distances in density, along the x-axis continuous values of likely (IV Î [0, +1]) and feasible (IC Î [−1, +1]) ranges. This should shed additional light on the locality of the discrepancies in the interest of future studies.
Lastly, we study the sensitivity of inference test results to the numerical noises embedded in Monte Carlo simulations. When the number of simulation paths increases, the replicability of the distribution by our parameterized model is expected to significantly improve. This effect can be demonstrated via an unpaired 2-sample Wilcoxon test, also known as Mann-Whitney test.
Table 4. 1-Path Kolomogorov-Smirnov (KS) test for equivalence between empirical and the simulated 1-path time series with an arbitrary starting point of IV0 and IC0.
.
.
Figure 5. 1-Path IV and IC permutation test.
Judging by the large test statistics W and the high p-values (>>10%) reported in Table 6, resulted in by comparing 5K paths simulation against the single-path historical realizations, one can’t reject the null hypothesis of identical distribution for both IV and IC at 95% confidence level.
5.2. Testing for Stage II
Our SDE Equation (4) assumes a bivariate normal distribution generated by
and
, two Wiener processes underlying IV and IC observables with an instantaneous and constant correlation parameter
. We now examine the normality of the model derived Brownian increments.
First we present the density function plots with a normal curve overlay, labeled as Normality Test (1) in Figure 6. We observe a very mild case of fat-tail (leptokurtic) in distributions of
and
.
Normality Test (2) is consisted of a pair of QQ-plots and ACF autocorrelation graphs. Figure 7 displays that a wide-ranged linearity and low autocorrelation have been tested in the time series of estimated residual, in support of the Gaussian distribution assumptions for
and
.
With a reasonably high degree of confidence in normality, we conclude that our proposed Stage II calibration method for
is consistent with empirical evidence.
5.3. Testing the Joint Distribution of IV and IC
Clearly Figure 2(a) makes a strong case in support of our USVSC model, as it successfully captures the very key feature of the real word implied volatility-correlation relationship as depicted by Figure 1(a).
To statistically measure the fit, we apply Fisher’s Z-transformation to test the Pearson correlation difference of two un-paired independent samples: one is the historical realization sample of IV and IC, and the other is estimated from the 5K paths of 2-dimensional simulation of IV and IC.
.
Figure 6.
and
normality tests (1).
Figure 7.
and
normality tests (2).
With a zero correlation difference as the null hypothesis, we obtain a Z-Score of 1.59 reported in Table 7. This test statistics indicates that the Pearson correlation estimation, between simulated IV and IC, is indistinguishable in distribution from that of the empirical observed IV and IC, over the sampled time period, at 95% confidence level.
Table 7. IV and IC dependency test.
.
We, therefore, conclude that the empirical joint distribution of IV and IC is replicated by our proposed USVSC model with sufficient accuracy.
6. Application of the Model—An Enhanced Geometric Brownian Motion
With the stochastic IV and IC fully parameterized and tested, we can now re-visit GBM model log-normal specification in Equation (1), the most foundational claim of the asset price returns in the finance literature.
Using GBM as the benchmark, Equation (1) parameters {μ, σ} are first calibrated through MLE. Table 8 gives the estimation results along with the corresponding standard errors.
Although the calibrated spot price volatility σ under GBM is robust and statistically significant, the same conclusion may not be drawn with respect to the parameter μ. The large estimation error suggests that the null hypothesis (μ = 0) is rather unrejectable with confidence, a subject to which we will return later.
Integrating our stochastic local volatility and correlation model, we proceed to replacing the constant diffusion coefficient σ with our time-varying and stochastic IV denoted by σt, as defined in Equation (2). The GBM Equation (1) thereby evolves into Equation (1a)—enhanced-GBM model, allowing the cash spot market infused with additional information sourced from the derivatives market:
(1a)
where
denotes a new Brownian motion for price returns, distinct from the original dBt from Equation (1).
To incorporate the proposed USVSC model, we further assume that the risk factors inherent to IV and IC markets are also integral to the dynamics of index trading, therefore can be thought of as priced state variables for asset pricing. This is not unreasonable, as the derivatives IV and IC have become standalone asset classes, albeit with a shorter history and a narrower investor base relative to the underlying cash market. VIX, for example, has evolved from being a latent risk parameter to a set of publicly investable trading product group largely helped by the creation of VIX Index and VIX futures commenced in 2004. Not unexpectedly, the stochasticity of IC is argued by Emmerich [16] as a fundamental source of risk, which is furthered by Driessen et al. [29] in a hunt for risk premium in the context of asset pricing. The roles of derivatives markets affecting index price discovering and the broader market have become more evident recently. Details and in-depth discussions on the market structure of VIX fixing, dispersion and worst-of option pricing, can be found in Osterrieder et al. [30], Li [31] and Zetocha [23].
Table 8. Equation (1) GBM model calibration.
Across three stylized state variables, i.e. the Brownian motion increments
together with dWt and dZt, a deterministic dependence structure is superimposed. This is manifested by Equation (10) and Equation (11) sequentially, with ρB and ρZ denoting the instantaneous correlation coefficients between spot-IV and spot-IC, respectively.
Since dWt and dZt are interconnected through ρw, as previously established by Equation (4), one can write:
(10)
with Wt □
as usual. Helped with the introduction of a third orthogonal factor
, we can now decompose the Weiner process
in a 3-independent-factor representation:
(11)
So that the recoveries of
,
and the identity property of
are satisfied. This particular setup should allow us to exploit the connections between index price dynamics with IV and IC observables5.
It should be noted that by now dWt and dZt have been attained from calibration Stage II. Therefore
and
are the two added sources of independent Gaussian white noise under the newly enhanced-GBM model Equation (1a).
Our task is therefore reduced to estimating the drift μ and correlation coefficients ρB and ρZ, and testing if the enhanced-GBM Equation (1a) is valid given our model setting. We proceed as follows.
Let xt = ln(St) to simplify the notation for asset price returns, after applying Ito’s lemma, Equation (1a) is expanded to a 3-factor model Equation (1b), admitting the information transformation from IV and IC markets rather explicitly:
(1b)
Taking the expectations on both sides, aided by the assumed Gaussianity of
, we arrive at the drift rate estimator
(12)
where all RHS variables are known, we obtain a calibration result and a statistical inference test in Table 9. At 95% confidence level, the null is once again unrejectable-same conclusion under the original GBM model. We thereby argue that the constant drift assumption under GBM (1) or (1a) is unsupported by the sampled historical data. That is, the linear drift term under GBM is not robust and likely misspecified, irrespective of the diffusion term being constant or locally stochastic.
It is clear now that the calibration of
and
is equivalent to locating the maximization of the log-likelihood (LL) function of a standard normal density over the full historical data sample:
(13)
Figure 8 pictorializes the calibration approach, once again by utilizing the Maximum Likelihood Estimation approach. The RHS expression of Equation (13) containing 2 correlation parameters is assumed to follow a standard Gaussian distribution. Hence the inverse of a max likelihood function, i.e. LL = −Log(L), reaches the minimum while the likelihood function L is maximized.
An implementation for a non-linear optimization algorithm yields a solution pair of
and
, displayed in a matrix form in Table 10 together with the estimated
for completeness:
The estimated negative correlation
, between spot and IV, corroborates the leverage effect puzzle termed by Ait-Sahalia et al. [33], and converges with the results by Zetocha [23]. This is also a widely known phenomenon by retail investors, as described in CBOE [34].
The negative correlation
, between spot and IC, comes as no big surprises as stock prices tend to correlate more as the market is down than the reverse, consistent with a unique and well-known IC market feature—correlation skew by Reghai [35], Zetocha [23] and Delanoe [36]. In the academic literature, however, the interactivities of spot-IC have not been extensively studied relative to spot-IV.
Several normality tests undertaken for
estimated under GBM (1), as well as
under GBM (1a): Table 11 gives the numerical improvements in Jargue Bera statistics; Figure 9 charts the normalized density plots benchmarked against the standard normal; Figure 10 supplements QQ-plots. The substantial gains in performance, of our joint stochastic Implied Volatility and Implied Correlation model over the original GBM model, are demonstratively significant.
Our enhanced GBM model of Equation (1a) may also underlie a portfolio selection process. To avoid penalizing upside volatility potential, a semivariance volatility approach could be applied, see for example Estrada [37] or Chen et al [38].
7. Concluding Summary
The dynamics of asset prices is typically modeled by stochastic processes; however their correlations are exogenously modeled and ad hoc assigned to the price process. This is mathematically and conceptually unsatisfying.
Table 9. GBM (1a) calibration and inference test of drift μ.
.
.
Figure 9. Density distribution comparison.
We propose a simple and unified framework to model local stochastic volatility and local stochastic correlation, similar to the approach by Heston [4]. Powered by coupled CIRCEV and Jacobi processes, our structural USVSC model is parsimonious and capable of re-producing the non-trivial volatility and correlation skews observed in the real world. The calibration scheme developed is intuitive and tractable, with results proven significant. The goodness of fit is backed by robust inference results after undertaking multiple stringent statistical tests in both marginal and joint distributions of the implied volatility and implied correlation.
We show our well-calibrated model performs markedly better over the classical GBM approach, thus a suitable and practical challenger for asset pricing, especially under current market structure.
Future research of our USVSC model may entail analyzing the sensitivity and stability to historical data sampling periods, and testing its capabilities in fitting the volatility smiles, term structure and correlation skews, for the purpose of pricing basket options and risk management. Additionally, new trading strategies may be developed if the equity option-implied correlation skews have awareness of the CDO tranche-implied correlation skews, across firm capital structure and asset classes. Lastly, placing the implied correlation and joint dynamics of volatility-correlation in a wider context of risk premium study should add useful insights in the field of financial economics.
NOTES
1The 3-D surface
is implied from vanillaoption pricesvia Black-Scholes-Merton formula, with strike K and maturity T.
2Correlations estimated are the average of Pearson 30 × 30 correlation matrix of Dow stock price returns.
3Strictly speaking, the VIX only corresponds to the estimated implied volatility of 30-day SPX call/put contracts through a convergence framework between variance swap pricing and log-Price option contract. In order to model VIX Options and Futures over tenor T, consistently with SPX options, it is most common to express VIX as the result of the integration of an instantaneous forward variance rate
, as in Cont et al. [11] and Bergormi [10], rather than the instantaneous volatility stochastic factor
as we proposed in Equation (2):
where T = 30 days. A similar expression exists for IC with respect to the quantity of
. To simplify the calibration, without evoking the instantaneous forward variance curve nor the required instantaneous forward correlation curve, we implicitly assume
and
for illustration purpose. Essentially a 1-month forward term-structure is embedded in both stochastic variables. This slight abuse of notation is acceptable, since 1) T is fixed at 30 days everywhere, and 2) the emphasis of our study is structured around the systemic interactions between IC and IV, impacts from term-structure and optimal risk hedging strategies are left out for future studies.
4We thank Prof. Boukhetala for R package Sim. Diff. Proc. and the useful discussions.
5By this design, the risk transmission from IV into spot market is channeled through correlated Brownian motions, as the stochastic
in Equation (1a) is set to historical observations directly. For additional insights of CKLS, readers are referred to Chan et al. [21] and Hu et al. [32].