Consistent and Efficient Pricing of SPX and VIX Options under Multiscale Stochastic Volatility
CConsistent and Efficient Pricing of SPX and VIX Options under MultiscaleStochastic Volatility
Jaegi Jeon a , Geonwoo Kim b , Jeonggyu Huh c, ∗ a Department of Mathematics, Yonsei University, Seoul 03722, Republic of Korea b School of Liberal Arts, Seoul National University of Science and Technology, Seoul 01811, Republic of Korea c School of Computational Sciences, Korea Institute for Advanced Study, Seoul 02455, Republic of Korea
Abstract
This study provides a consistent and efficient pricing method for both Standard & Poor’s 500 Index (SPX)options and the Chicago Board Options Exchange’s Volatility Index (VIX) options under a multiscalestochastic volatility model. To capture the multiscale volatility of the financial market, our model adds afast scale factor to the well-known Heston volatility and we derive approximate analytic pricing formulasfor the options under the model. The analytic tractability can greatly improve the efficiency of calibrationcompared to fitting procedures with the finite difference method or Monte Carlo simulation. Our experimentusing options data from 2016 to 2018 shows that the model reduces the errors on the training sets of theSPX and VIX options by 9.9% and 13.2%, respectively, and decreases the errors on the test sets of theSPX and VIX options by 13.0% and 16.5%, respectively, compared to the single-scale model of Heston. Theerror reduction is possible because the additional factor reflects short-term impacts on the market, which isdifficult to achieve with only one factor. It highlights the necessity of modeling multiscale volatility.
Keywords:
SPX option; VIX option; multiscale volatility; Asymptotic method
1. Introduction
Volatility in the financial market is one of the most important measures for investors. The VolatilityIndex (VIX) of the Chicago Board Options Exchange (CBOE) was introduced as a measure to capture thevolatility of the Standard & Poor’s 500 Index (SPX) in 1993. The VIX was revised to be independent of anymodel in 2003. After the revision, it is calculated using SPX option prices and the SPX over the next 30calendar-day period. Since the introduction of the VIX, the index has become a standard gauge of marketfear, because the VIX and the SPX are negatively correlated. In addition, VIX options also have receiveda lot of attention as the popularity of the VIX has grown, and this interest in the derivatives has led to thedevelopment of an appropriate pricing method for them.It is noteworthy that many researchers have invented pricing methods for VIX options under the modelsoriginally proposed to price SPX options. This may be because an arbitrage relationship exists betweenthe SPX option market and the VIX option market. In the early research period for VIX options, therewere the studies to derive pricing formulas for VIX options under several one-factor affine jump diffusion(AJD) models; Lin and Chang (2009), Lian and Zhu (2013), Lin and Chang (2010), Cheng et al. (2012) andCont and Kokholm (2013). As it is widely known, the one-factor AJD model has a special structure froma mathematical standpoint, so they often give a pricing formula for a derivative for SPX options as wellas VIX options (refer Duffie et al., 2000). It is highly desirable in the sense that they make it possible toevaluate the SPX options and the VIX options without arbitrage. However, it is recognized that the one-factor AJD models have a variety of fundamental limitations. Above all, the models cannot accommodate ∗ corresponding author. Tel: 82-2-958-3750E-mail address: aifi[email protected] Preprint submitted to arXiv September 24, 2019 a r X i v : . [ q -f i n . M F ] S e p he multiscale property of volatility, which highlights the necessity of multi-factor affine diffusion models.Jumps are usually excluded in the models, but it is frequently stated that the models are superior to theone-factor AJD models in literature such as Kaeck and Alexander (2012).In this light, Gatheral (2008) proposed the double-mean-reverting (DMR) Heston model to capturethe multiscales of volatility. But it does not provide any pricing formula and depends on the Monte-Carlosimulation totally for a pricing and a calibration. Although the simulation method has been greatly improvedin Bayer et al. (2013), the model still fails to suggest a viable solution to a calibration using a large numberof option prices. On the other hand, Fouque and Saporito (2018) proposed the Heston stochastic vol-of-vol model (Heston SVV model) with fast and slow time scales to reflect skews in both the SPX data andVIX data, and derived approximation solutions for the prices of options on the SPX and VIX, respectively.Moreover, it showed that the Heston SVV model calibrated jointly to real market data for the SPX andVIX fits to the data very well. Their solutions, however, cannot be also computed within a short period oftime. Interestingly, Huh et al. (2018) obtained an approximate quasi-closed formula for VIX options underthe DMR Heston model, and furthermore, it confirmed with real market data that the calibration resultutilizing the approximate solution is analogous to the case using simulation methods. But the approach canbe applied to pricing VIX options only.Our novel model is the first multi-factor affine diffusion model to give practically feasible pricing formulasfor both SPX options and VIX options. To be concrete, we propose a new two-factor affine diffusion model,and induce approximate quasi-closed formulas for both the options based upon the perturbation theory. Theanalytic tractability greatly improve the efficiency of calibration. As a result, we could achieve satisfactorycalibration results with long-term data for the SPX and VIX from 2016 to 2018, consisting of hundreds ofthousands of options.The remainder of the paper is organized as follows. We introduce a new two-factor model with multiscalevolatility to jointly price VIX options and SPX options in Section 2. We derive the analytic pricing formulasof SPX options and VIX options under the proposed model using asymptotic methods in Section 3. Wepresent calibration with real options data for the SPX and VIX in Section 4. We provide concluding remarksin Section 5.
2. A new two-factor multiscale stochastic volatility model
Before introducing our model, we briefly review the double Heston model of Christoffersen et al. (2009).The model was originally devised to reflect stochastic correlation as well as stochastic volatility. The im-provement is achievable because the model has two factors. By contrast, the Heston model cannot capturestochastic correlation because it has only one factor. The double Heston model, however, cannot capture themultiscale characteristics of a market, because it was not designed for this purpose (see Gallant et al., 1999;Chernov et al., 2003). Thus, the model needs to be extended to depict well-separated time scales. Moreimportantly, the model provides an analytic formula for VIX options, which requires heavy computation fora double integral. In particular, when practitioners want to fit the model to prices for both SPX optionsand VIX options consistently within a limited time, the time cost raises a severe problem.To overcome these limitations of the double Heston model, we propose a two-factor multiscale stochasticvolatility model, which is given by the following stochastic differential equations: dX t = rX t dt + (cid:112) Y t X t dW X, t + (cid:112) Z t X t dW X, t ,dY t = 1 (cid:15) ( Z t − Y t ) dt + √ ν √ (cid:15) (cid:112) Y t dW Yt , (1) dZ t = κ ( θ − Z t ) dt + σ (cid:112) Z t dW Zt for 0 < (cid:15) (cid:28) Q , where r is the risk-free rate, and κ (cid:28) /(cid:15) is assumed. Thecorrelation structure of our model is dW X, t dW Yt = ηdt, dW X, t dW Zt = ρdt,dW X, t dW Zt = dW X, t dW Yt = dW X, t dW X, t = dW Yt dW Zt = 0 .
2t should be noted that Y t rapidly reverts to Z t under our model, which is a major difference with thedouble Heston model. Because the time scales are clearly separated due to condition κ (cid:28) /(cid:15) , the short-term variance ourY t + Z t rapidly reverts to mid-term variance 2 Z t , and 2 Z t regresses to long-term variance2 θ at a relatively slow rate. On the other hand, one can consider that our model is obtained by adding a fastmean-reverting factor to Heston’s volatility or otherwise by modifying the double Heston model to expressthe well-separated time scales.
3. Pricing SPX and VIX options
By the risk-neutral pricing principle, the price P (cid:15)s for an SPX call option with strike K and maturity T at time t is P (cid:15)s ( t, x, y, z ) = e − r ( T − t ) E Q x,y,z (cid:104) ( X T − K ) + (cid:105) , (2)where ( X t , Y t , Z t ) = ( x, y, z ). Applying the Feynman–Kac theorem to our model (1), it can be shown that P (cid:15)s should satisfy the following partial differential equation (PDE): (cid:18) (cid:15) L s, + 1 √ (cid:15) L s, + L s, (cid:19) P (cid:15)s ( t, x, y, z ) = 0 , where L s, = ( z − y ) ∂ y + ν y∂ yy , L s, = √ ηνxy∂ xy , (3) L s, = ∂ t + rx∂ x + κ ( θ − z ) ∂ z + 12 x ( y + z ) ∂ xx + 12 σ z∂ zz + ρσxz∂ xz − r · with the final condition P (cid:15)s ( T, x, y, z ) = ( x − K ) + . Deriving an analytic solution for this PDE would bedesirable; however, we confirm that the solution cannot be easily induced.Thus, we attempt to obtain an approximate solution for P (cid:15)s based on asymptotic methods (refer toFouque et al. (2011)). Let us expand P (cid:15)s with respect to √ (cid:15) , that is, P (cid:15)s = P s, + √ (cid:15)P s, + (cid:15)P s, + · · · . Wewish to approximate P (cid:15)s to the sum of the leading term P s, and the first non-zero correction term P s, . Todo this, we show that P s, and P s, follow the following PDEs independent of y : • leading term P s, (cid:104)L s, (cid:105) P s, = 0 ,P s, ( T, x, z ) = ( x − K ) + , (4) • first non-zero correction term P s, (cid:104)L s, (cid:105) P s, = W (cid:15) zx∂ x (cid:0) x ∂ xx (cid:1) P s, ,P s, ( T, x, z ) = 0 , (5)where (cid:104)L s, (cid:105) = ∂ t + rx∂ x + κ ( θ − z ) ∂ z + x z∂ xx + 12 σ z∂ zz + ρσxz∂ xz − r,W (cid:15) = − √ ην √ (cid:15).
3f setting 2 z = ξ , the operator (cid:104)L s, (cid:105) changes to L H = ∂ t + rx∂ x + κ (2 θ − ξ ) ∂ ξ + 12 x ξ∂ xx + 12 (cid:16) √ σ (cid:17) ξ∂ ξξ + (cid:18) √ ρ (cid:19) (cid:16) √ σ (cid:17) xξ∂ xξ − r. The resulting operator L H can be then regarded as that generated from the Heston model whose parametersare the mean-reverting rate κ , the long-term variance level 2 θ , the volatility of volatility √ σ , and theleverage correlation ρ/ √ P s, and P s, aresummarized in the following theorem. Theorem 1.
The price P (cid:15)s in (2) can be approximated with an accuracy of order O ( (cid:15) ) , that is, | P (cid:15)s − ( P s, + P s, ) | < C(cid:15) for some C > . Furthermore, P s, and P s, have the following forms: P s, ( t, x, z ) = e − rτ π (cid:90) e − ikq ˆ G ( τ, k, z ) ˆ h ( k ) dk,P s, ( t, x, z ) = e − rτ π (cid:90) e − ikq b ( k ) (cid:16) κθ ˆ f ( τ, k ) + 2 z ˆ f ( τ, k ) (cid:17) ˆ G ( τ, k, z ) ˆ h ( k ) dk, where τ = T − t , q = rτ + log x , ˆ G ( τ, k, y ) = e C ( τ,k )+ yD ( τ,k ) , ˆ h ( k ) = K ik ik − k , ˆ f ( τ, k ) = 2 τ d ( k ) g ( k ) + g ( k ) − d ( k ) g ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) + τ d ( k ) g ( k ) − g ( k ) − d ( k ) g ( k ) , ˆ f ( τ, k ) = e τd ( k ) (cid:16) g ( k ) (cid:0) e τd ( k ) − (cid:1) − τ d ( k ) g ( k ) + 1 (cid:17) − d ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) , and C ( τ, k ) = κθσ (cid:18) ( κ + ρikσ − d ( k )) τ − (cid:18) e − τd ( k ) − g ( k )1 − e − τd ( k ) g ( k ) (cid:19)(cid:19) ,D ( τ, k ) = κ + ρikσ + d ( k )2 σ (cid:18) e − τd ( k ) − e − τd ( k ) − g ( k ) (cid:19) ,b ( k ) = − W (cid:15) (cid:0) ik + k (cid:1) ,d ( k ) = (cid:113) σ ( k − ik ) + ( κ + ρikσ ) ,g ( k ) = κ + ρikσ + d ( k ) κ + ρikσ − d ( k ) . Proof.
The detailed proof is in Appendix A. 4 .2. Relationship between the VIX and our model
Before proceeding, we show the relationship between the VIX and our model. Because the VIX at time t , VIX (cid:15)t , is the square root of an integration of spot variance from t to t + τ ( τ := 30 / (cid:15)t ) = 100 × E Q (cid:20) τ (cid:90) t + τ t ( Y s + Z s ) ds (cid:21) = 100 × τ (cid:90) t + τ t (cid:0) E Q [ Y s ] + E Q [ Z s ] (cid:1) ds (6)= 100 × ( a (cid:15) Y t + a (cid:15) Z t + ( a (cid:15) + a (cid:15) ) θ ) , where a (cid:15) = (cid:15)τ (cid:16) − e − τ /(cid:15) (cid:17) ,a (cid:15) = 1 κτ (cid:0) − e − κτ (cid:1) + 11 − κ(cid:15) (cid:18) κτ (cid:0) − e − κτ (cid:1) − (cid:15)τ (cid:16) − e − τ /(cid:15) (cid:17)(cid:19) ,a (cid:15) = 1 − a (cid:15) ,a (cid:15) = 1 − a (cid:15) . The following facts induced by Ito’s lemma are used for the derivation of the above formula (6).E [ Y s ] = e − (cid:15) ( s − t ) Y t + 11 − κ(cid:15) (cid:16) e − κ ( s − t ) − e − (cid:15) ( s − t ) (cid:17) Z t + (cid:20)(cid:16) − e − (cid:15) ( s − t ) (cid:17) − − κ(cid:15) (cid:16) e − κ ( s − t ) − e − (cid:15) ( s − t ) (cid:17)(cid:21) θ, E [ Z s ] = e − κ ( s − t ) Z t + θ (cid:16) − e − κ ( s − t ) (cid:17) . It is noteworthy that (VIX (cid:15)t ) converges to(VIX ∗ t ) = 100 × ( a ∗ Z t + (1 + a ∗ ) θ ) , as (cid:15) →
0, where a ∗ = κτ (1 − e − κτ ) and a ∗ = 1 − a ∗ . In addition, VIX t and VIX ∗ t have the followingrelationship:(VIX (cid:15)t ) − (VIX ∗ t ) = 100 × (cid:18) (cid:15)τ (cid:16) − e − τ /(cid:15) (cid:17) ( Y t − Z t ) + (cid:15)τ (cid:0) − e − κτ (cid:1) ( Z t − θ ) (cid:19) + O (cid:0) (cid:15) (cid:1) := ∆VIX t . Denoting the price for a VIX call option with strike K and maturity T as P (cid:15)v at time t , P (cid:15)v is expressedby P (cid:15)v ( t, y, z ) = e − r ( T − t ) E Q y,z [ h (cid:15) ( Y T , Z T )] , (7)where ( Y t , Z t ) = ( y, z ), and h (cid:15) ( u, v ) = (VIX (cid:15)T ( u, v ) − K ) + . If H ( x ) := ( √ x − K ) + , h (cid:15) ( u, v ) = H ((VIX (cid:15)T ) ).Then, by expanding H at (VIX ∗ t ) , h (cid:15) can be transformed to h (cid:15) ( u, v ) = H ((VIX ∗ T ) ) + H (cid:48) ((VIX ∗ T ) )∆VIX t + 12 H (cid:48)(cid:48) ((VIX ∗ T ) ) (cid:0) ∆VIX t (cid:1) + · · · = ( (cid:113) a ∗ v + (1 + a ∗ ) θ × − K ) { v ≥ ( K − ( a ∗ ) θ ) /a ∗ } + (cid:15) (cid:32) (cid:0) − e − τ /(cid:15) (cid:1) ( u − v ) + (1 − e − κτ ) ( v − θ )2 τ (cid:112) a ∗ v + (1 + a ∗ ) θ (cid:33) { v ≥ ( K − ( a ∗ ) θ ) /a ∗ } ×
100 + O (cid:0) (cid:15) (cid:1) = h ( v ) + (cid:15)h ( u, v ) + O (cid:0) (cid:15) (cid:1) , h and h are the first and second terms of the second line, respectively. Meanwhile, utilizing Ito’slemma, we can show Y s = Z s + e − (cid:15) ( s − t ) ( Y t − Z t ) + √ ν √ (cid:15) (cid:90) st e − (cid:15) ( s − u ) (cid:112) Y u dW Yu + ∞ (cid:88) n =1 ( κ(cid:15) ) n (cid:104) ( Z s − θ ) − e − (cid:15) ( s − t ) ( Z t − θ ) (cid:105) − σ ∞ (cid:88) n =0 ( κ(cid:15) ) n (cid:90) st e − (cid:15) ( s − u ) (cid:112) Z u dW Zu . Using the tower property of conditional expectation and the linearity of h with respect to u , we obtainE Q y,z [ h ( Y T , Z T )] = E Q y,z (cid:2) E Q y,z [ h ( Y T , Z T ) | Z T ] (cid:3) = E Q y,z (cid:2) E Q y,z [ h ( Z T ) + (cid:15)h ( Y T , Z T ) | Z T ] (cid:3) + O (cid:0) (cid:15) (cid:1) = E Q y,z (cid:2) h ( Z T ) + (cid:15)h (cid:0) E Q y,z [ Y T | Z T ] , Z T (cid:1)(cid:3) + O (cid:0) (cid:15) (cid:1) = E Q y,z (cid:104) h ( Z T ) + (cid:15)h (cid:16) Z T + e − ( T − t ) /(cid:15) ( y − z ) + O ( (cid:15) ) , Z T (cid:17)(cid:105) + O (cid:0) (cid:15) (cid:1) = E Q y,z [ h ( Z T ) + h ∗ ( Z T )] + O (cid:0) (cid:15) (cid:1) , where h ∗ ( v ) = (cid:15)h (cid:0) v + e − ( T − t ) /(cid:15) ( y − z ) , v (cid:1) . Then, the facts presented so far lead to the following theorem. Theorem 2.
The price P (cid:15)v in (7) can be approximated with accuracy of order O ( (cid:15) ) , that is, | P (cid:15)v − ( P v, + P v, ) | < C(cid:15) . for some C > . Furthermore, P v, and P v, are given by P v, ( t, z ) = e − rT (cid:90) ∞ h ( δζ ) f ( ζ ; k, λ ) dζ,P v, ( t, y, z ) = e − rT (cid:90) ∞ h ∗ ( δζ ) f ( ζ ; k, λ ) dζ, where f ( ζ ; k, λ ) is the density function of a non-central chi-square distribution with k = κθσ degrees offreedom and non-centrality parameter λ = ze − κτ δ with τ = T − t , where δ = ( − e − κτ ) σ κ , and h ( v ) = ( (cid:113) a ∗ v + (1 + a ∗ ) θ × − K ) + { v ≥ ( K − ( a ∗ ) θ ) /a ∗ } ,h ∗ ( v ) = (cid:15) (cid:32) e − τ/(cid:15) a (cid:15) ( y − z ) + κ(cid:15)a ∗ ( v − θ )4 (cid:112) a ∗ v + (1 + a ∗ ) θ (cid:33) { v ≥ ( K − ( a ∗ ) θ ) /a ∗ } × . Proof.
Because it is known that future values of the Cox–Ingersoll–Ross (CIR) process Z t follow a non-centralchi-squared distribution, the values of P v, and P v, can be expressed as the expectations with respect tothe distribution (c.f. Brigo and Mercurio, 2007).Note that the above analytic formulas in Theorem 2 are involved in single Integrals, as in the case of theformulas for SPX options. This means that it is possible to compute the values without difficulty. Wealso strongly stress that our model has only six parameters κ , θ , σ , ρ , (cid:15) , and W (cid:15) , which is only two morethan the Heston model. Some of the parameters ( W (cid:15) , ρ ) and (cid:15) control only of SPX and VIX option prices,respectively.
4. Empirical test
In this section, we verify the outperformance of our model through an empirical test using real marketdata. To show the benefits of two-factor models clearly, we choose the Heston model, which is the most6 S P X SPXVIX 101520253035 V I X Figure 1: Series of the SPX and the VIX for 2016–2018. The vertical red line indicates the dividing date, the left and the rightof which correspond to the training period (2016–2017) and the test period (2018), respectively.
SPX VIXmean deviation skewness kurtosis mean deviation skewness kurtosistraining 0.001 0.007 -0.530 4.379 -0.001 0.073 0.816 6.230test 0.000 0.011 -0.435 3.319 0.004 0.100 2.170 13.289all 0.000 0.008 -0.581 5.295 0.000 0.083 1.669 12.385
Table 1: Moments for the log-return series of the SPX and the VIX for the data period. well-known one-factor model, as a benchmark. While it may be considered more reasonable to compareour model with another two-factor model, there is no existing two-factor model that gives analytic pricingformulas that are available and efficient for both SPX and VIX options. For example, the double Hestonmodel (Christoffersen et al., 2009) provides a double integral form for the pricing of VIX options but cannotbe considered in practice because it takes a huge amount of time to compute, as its parameters are foundbased on a lot of option prices.We obtain options data on the SPX and the VIX from 2016 to 2018 from the CBOE data shop . Wefilter them out if their daily trades are below 50, if their prices do not reach 0.5 point, or if they expirewithin 3 days. We then split the data into a training set for 2016–2017 (235,227 SPX options and 46,516VIX options) and a test set for 2018 (176,866 SPX options and 23,347 VIX options). Figure 1 shows theseries of the SPX and the VIX for the data period. The vertical red line indicates the dividing date, the leftand the right of which correspond to the training period and the test period, respectively. Many financialcrises arose during the period: Chinese stock market turbulence (early 2016), Brexit (late 2016), and USstock market downturn (2018). It is clear that the SPX falls sharply and the VIX increases dramaticallyduring each crisis. Table 1 summarizes the moments for the log-returns of both series for the data period.The VIX returns move in a much more volatile way than the SPX returns move. In addition, the SPXreturns show negative skewness and weak kurtosis ( > (cid:29) κ , θ , σ , and (cid:15) to control VIX option prices, and the other parameters ρ , W (cid:15) are soughtwith SPX option prices. Note that the VIX option prices may be more sensitive to the parameters κ , θ , and σ than the SPX option prices. Because a highly precise optimization needs highly sensitive Greeks, it seems https://datashop.cboe.com/ κ , θ , and σ . We firstdefine the method for the Heston model with the four parameters κ, θ, σ , and ρ as follows.1. First, perform calibration with VIX optionsargmin κ,θ,σ (cid:88) i ∈ I (cid:88) j ∈ J v,i (cid:18) P mdlv ( z i , κ, θ, σ ; t i, K j , T j ) − P mktv ( t i, , K j , T j )0 . P mktv ( t i, , K j , T j ) (cid:19) such that VIX mkti = (cid:112) b ∗ z i + b ∗ θ. (8)where b ∗ = a ∗ / b ∗ = (1 + a ∗ ) / ρ (cid:88) i ∈ I (cid:88) j ∈ J s,i (cid:18) P mdls ( ρ ; t i, K j , T j ; z i , κ, θ, σ ) − P mkts ( t i, , K j , T j )0 . P mkts ( t i, , K j , T j ) (cid:19) . The superscripts mkt and mdl are used to indicate the market price and the model price for an option,respectively. I is the index set of the dates for the training set. J v,i and J s,i are the index sets of the VIXoptions and the SPX options, respectively, on the i th date. We explain the two-step method for our modelwith the six parameters κ, θ, σ, ρ, (cid:15) , and W (cid:15) as follows.1. First, perform calibration with VIX optionsargmin κ,θ,σ,(cid:15) (cid:88) i ∈ I argmin y i ,z i (cid:88) j ∈ J v,i (cid:18) P mdlv ( y i , z i , κ, θ, σ, (cid:15) ; t i, K j , T j ) − P mktv ( t i, , K j , T j )0 . P mktv ( t i, , K j , T j ) (cid:19) such that VIX mkti = (cid:113) a (cid:15) y i + a (cid:15) z i + ( a (cid:15) + a (cid:15) ) θ. (9)2. Second, perform calibration with SPX optionsargmin ρ,W (cid:15) (cid:88) i ∈ I (cid:88) j ∈ J s,i (cid:18) P mdls ( ρ, W (cid:15) ; t i, K j , T j ; y i , z i , κ, θ, σ ) − P mkts ( t i, , K j , T j )0 . P mkts ( t i, , K j , T j ) (cid:19) . Note that, under the Heston model, z i is uniquely determined by the relationship (8) between the modeland the VIX. By contrast, under our model, y i and z i cannot be uniquely determined by the relationship(9), which is why we also minimize the daily objective with respect to y i and z i for each i th date. Becausethis process requires large computational resources, the calibrations for both models are implemented basedupon C++ and OpenMP, and executed in parallel on 24 CPU cores of two Intel Xeon Gold 5118. As aresult, we obtain the following parameters for each model: κ = 3 . θ = 0 . σ = 0 . ρ = − κ = 3 . θ = 0 . σ = 0 . ρ = − (cid:15) = 0 . W (cid:15) = 0 . κ = 3 .
58 and 1 /(cid:15) = 104 . √ z i for the8 (a) errors for SPX options (b) errors for VIX optionsFigure 2: Comparison of daily errors of both models by product types for 2016–2018. The vertical red line indicates the dividingdate, the left and the right of which are the training period (2016–2017) and the test period (2018), respectively. (a) √ z i for the Heston model (b) √ y i + z i for our modelFigure 3: Spot volatilities of both models for 2016–2018. The vertical red line indicates the dividing date. Heston model is strongly correlated with the VIX but √ y i + z i for our model is not. Our model capturesshort-term volatility that is difficult to detect and produces a more volatile process. If our model is assumedto be correct, we conclude that the spot volatilities for the Heston model are fairly often biased, especiallywhen sudden strong shocks impact the market. The bias eventually results in large fitting or predictionerrors of the Heston model for short-term products, as shown in Table Appendix C. The table sums up theerrors separately by time to maturity. The table results confirm that our model performs better than theHeston model as the time to maturity for an option becomes shorter, provided that the time to expirationis not too short. Note that the asymptotic method is valid when the time to maturity is not too short.Furthermore, implied volatility surfaces for SPX options and VIX options clearly show how appropriateour model is for fitting to short-term products. To define the implied volatility for VIX options consistently, we first consider the following simple model: d VIX t = σdW t . In other words, VIX T follows the normal distribution N (cid:0) VIX t , σ √ T − t (cid:1) under the model. Then, it is easyto show that the price c of a VIX call option with maturity T and strike K is given by c (VIX t ) = (VIX t − K ) N (cid:18) VIX t − Kσ √ T − t (cid:19) + σ √ T − t n (cid:18) VIX t − Kσ √ T − t (cid:19) where N and n are the cumulative distribution function and the probability density function of the unitGaussian distribution, respectively. As a result, for the market price c mkt , the corresponding impliedvolatility for the VIX option can be defined as σ imp to match the model price with c mkt .9 (a) y = 0 . z = 0 . T 0.00.51.0 K 9095100105 0.120.100.080.060.040.020.000.020.04 (b) y = 0 . z = 0 . z = 0 . T 0.00.51.0 K 17.019.522.025.0 0.050.040.030.020.010.00 (a) y = 0 . z = 0 . T 0.00.51.0 K 17.019.522.025.0 0.010.020.030.040.05 (b) y = 0 . z = 0 . We now generate implied volatility surfaces σ imp , from the corrected prices P s, ( z )+ P s, ( z ) and P v, ( z )+ P v, ( y, z ) with the estimated parameters κ = 3 . θ = 0 . σ = 0 . ρ = − (cid:15) = 0 . W (cid:15) = 0 . y, z ) = (0 . , . y, z ) = (0 . , . σ imp from the uncorrected prices P s, ( z ) and P v, ( z ) for hiddenstate z = 0 . (cid:15) = 0 and W (cid:15) = 0. All the hidden states y , z are chosen so that themodel value of the VIX is 20, that is, VIX t = 20 in the relationship (6).Figures 4 and 5 present the differences σ imp , − σ imp for the SPX cases and the VIX cases, respectively.The figures confirm once again that the corrected terms P s, and P v, allow excellent short-term flexibility,and thereby capture the fluctuations in the slope and the level of the volatility smirk. Based on Figure 4, P s, may make short-term in-the-money (SITM) options more expensive than short-term out-of-the-money(SOTM) options in a somewhat robust way. However, based on Figure 5, a subtle change of z in our modelseems to lead to a significantly different short-term VIX market. When the short-term state y is smallerthan the mid-term state z , as in the right figure, P v, makes the SITM options more expensive than the10OTM options, as in the SPX cases, while P v, works in the opposite way when y is larger than z , as inthe left figure. These phenomena for VIX volatility surfaces may occur because Y t quickly converges Z t inour model. If y > z , the spot volatility might decline in the short term, which means that the VIX SITMcalls would probably not be exercised. If y < z , the spot volatility might increase for a short while, whichwould increase the value of the VIX SOTM call. This is simply the mechanism of our model showing diverseexpression for the short-term VIX market in spite of the constraint (6).
5. Concluding remarks
In this paper, we study consistent and efficient pricing of SPX and VIX options under a new two-factorstochastic volatility model. Specifically, this two-factor model is proposed by adding a fast mean-revertingfactor into the Heston model. Doing so facilitates joint pricing of the SPX option and the VIX option.In practice, joint modeling of both options is important, because an arbitrage relationship exists betweenthe SPX option market and the VIX option market. Moreover, joint modeling leads to a calibration basedon extensive market data, including SPX data and VIX data. Since our analytic solutions are derived asone-dimensional integrals, it is obvious that the pricing solutions are computationally very efficient. Ourexperiment using hundreds of thousands of options shows that the model reduces the errors by 9.9-16.5%,compared to the single-scale model of Heston. The error reduction is possible because the additional factorreflects short-term impacts on the market, which is difficult to achieve with only one factor.In fact, non-affine models have been less studied for explicit pricing formulas because the involvedproblems are hard to solve. But the models are more suitable to express actual dynamics (see Mencia andSentana (2013) and Kaeck and Alexander (2013)). In this context, our models should be extended to havea non-affine form. We leave the topic as future work.
Acknowledgment
Jaegi Jeon received financial support from the National Research Foundation of Korea (NRF) of theKorean government (Grant No. NRF-2019R1I1A1A01062911). Geonwoo Kim received financial supportfrom the NRF (Grant No. NRF-2017R1E1A1A03070886). Jeonggyu Huh received financial support fromthe NRF (Grant No. NRF-2019R1F1A1058352).
Declarations of interest : none.ReferencesReferences
Bayer, C., Gatheral, J., Karlsmark, M., 2013. Fast ninomiya–victoir calibration of the double-mean-reverting model. Quanti-tative Finance 13 (11), 1813–1829.Brigo, D., Mercurio, F., 2007. Interest rate models-theory and practice: with smile, inflation and credit. Springer Science &Business Media.Carr, P., Madan, D., 1998. Towards a theory of volatility trading. Volatility: New estimation techniques for pricing derivatives29, 417–427.Cheng, J., Ibraimi, M., Leippold, M., Zhang, J. E., 2012. A remark on lin and chang’s paper ‘consistent modeling of s&p 500and vix derivatives’. Journal of Economic Dynamics and Control 36 (5), 708–715.Chernov, M., Gallant, A. R., Ghysels, E., Tauchen, G., 2003. Alternative models for stock price dynamics. Journal of Econo-metrics 116 (1-2), 225–257.Christoffersen, P., Heston, S., Jacobs, K., 2009. The shape and term structure of the index option smirk: Why multifactorstochastic volatility models work so well. Management Science 55 (12), 1914–1932.Cont, R., Kokholm, T., 2013. A consistent pricing model for index options and volatility derivatives. Mathematical Finance:An International Journal of Mathematics, Statistics and Financial Economics 23 (2), 248–274.Duffie, D., Pan, J., Singleton, K., 2000. Transform analysis and asset pricing for affine jump-diffusions. Econometrica 68 (6),1343–1376.Fouque, J.-P., Lorig, M. J., 2011. A fast mean-reverting correction to heston’s stochastic volatility model. SIAM Journal onFinancial Mathematics 2 (1), 221–254. ouque, J.-P., Papanicolaou, G., Sircar, R., Sølna, K., 2011. Multiscale stochastic volatility for equity, interest rate, and creditderivatives. Cambridge University Press.Fouque, J.-P., Saporito, Y. F., 2018. Heston stochastic vol-of-vol model for joint calibration of vix and s&p 500 options.Quantitative Finance 18 (6), 1003–1016.Gallant, A. R., Hsu, C.-T., Tauchen, G., 1999. Using daily range data to calibrate volatility diffusions and extract the forwardintegrated variance. Review of Economics and Statistics 81 (4), 617–631.Gatheral, J., 2008. Consistent modeling of spx and vix options. In: Bachelier Congress. p. 3.Goutte, S., Ismail, A., Pham, H., 2017. Regime-switching stochastic volatility model: estimation and calibration to vix options.Applied Mathematical Finance 24 (1), 38–75.Huh, J., Jeon, J., Kim, J.-H., 2018. A scaled version of the double-mean-reverting model for vix derivatives. Mathematics andFinancial Economics 12 (4), 495–515.Kaeck, A., Alexander, C., 2012. Volatility dynamics for the s&p 500: Further evidence from non-affine, multi-factor jumpdiffusions. Journal of Banking & Finance 36 (11), 3110–3121.Kaeck, A., Alexander, C., 2013. Continuous-time vix dynamics: On the role of stochastic volatility of volatility. InternationalReview of Financial Analysis 28, 46–56.Lian, G.-H., Zhu, S.-P., 2013. Pricing vix options with stochastic volatility and random jumps. Decisions in Economics andFinance 36 (1), 71–88.Lin, Y.-N., Chang, C.-H., 2009. Vix option pricing. Journal of Futures Markets: Futures, Options, and Other DerivativeProducts 29 (6), 523–543.Lin, Y.-N., Chang, C.-H., 2010. Consistent modeling of s&p 500 and vix derivatives. Journal of Economic Dynamics and Control34 (11), 2302–2319.Mencia, J., Sentana, E., 2013. Valuation of vix derivatives. Journal of Financial Economics 108 (2), 367–391. Appendix A. Derivation of analytic formula for SPX options
Putting P (cid:15)s = P s, + √ (cid:15)P s, + (cid:15)P s, + (cid:15) √ (cid:15)P s, + · · · into the PDE (3),1 (cid:15) L s, P s, + 1 √ (cid:15) ( L s, P s, + L s, P s, ) + ( L s, P s, + L s, P s, + L s, P s, )+ √ (cid:15) ( L s, P s, + L s, P s, + L s, P s, ) + · · · = 0 . As mentioned in Theorem 1, we intend for the sum of the leading term P s, and the first non-zero correctionterm P s, to approximate P (cid:15)s within accuracy O ( (cid:15) ), that is, | P (cid:15)s − ( P s, + √ (cid:15)P s, ) | < C(cid:15) for some positive C . For this purpose, the following equations should be satisfied: L s, P s, = 0 , (A.1) L s, P s, + L s, P s, = 0 , (A.2) L s, P s, + L s, P s, + L s, P s, = 0 , (A.3) L s, P s, + L s, P s, + L s, P s, = 0 . (A.4) P s, should not depend on y so that the Poisson equation (A.1) has a solution. This means that P s, alsodoes not depend on y owing to (A.2). From (A.3), we then obtain L s, P s, + L s, P s, = 0 . (A.5)Thus, the centering condition for L s, for the Poisson equation (A.5) results in the following PDE (4) for P s, : (cid:104)L s, P s, (cid:105) = (cid:104)L s, (cid:105) P s, = 0 , (A.6) P s, ( T, x, z ) = ( x − K ) + , where (cid:104)·(cid:105) is the expectation with respect to the invariant distribution Φ for Y t , that is, (cid:104) f (cid:105) = (cid:82) f ( y ) Φ ( y ) dy .Moreover, because Y t is a CIR process, Y ∞ ∼ Γ (cid:0) γ, ν (cid:1) , which means that the density function for Y ∞ isΦ ( w ) = 1 ν γ Γ ( γ ) w γ − e − ( w/ν ) (0 , ∞ ) ( w ) , γ = zν . The centering condition indicates that (cid:104) g (cid:105) = 0 should hold for Poisson equation L f + g = 0.The condition is necessary for a Poisson equation to have a solution.In addition, (A.5) yields the following equation: P s, = −L − s, ( L s, − (cid:104)L s, (cid:105) ) P s, + c ( t, x, z )= − x φ ( y, z ) ∂ xx P s, + c ( t, x, z ) , (A.7)where L s, φ ( y, z ) = y − z. (A.8)In fact, we show φ ( y, z ) = z − y , the proof of which is given in Appendix B. On one hand, if we put (A.7)into (A.4) and use the centering condition for L s, in (A.4), we achieve the PDE (5) for P s, in the followingways: (cid:104)L s, (cid:105) P s, = − (cid:104)L s, P s, (cid:105) = 1 √ ην (cid:104) yφ y ( y, z ) (cid:105) x∂ x (cid:0) x ∂ xx P (cid:1) = W zx∂ x (cid:0) x ∂ xx (cid:1) P s, . (A.9)with the corresponding final condition P s, ( T, x, z ) = 0 , where W = − √ ην (recall that W (cid:15) = − √ ην √ (cid:15) , i.e., W = W (cid:15) / √ (cid:15) ).On the other hand, if ξ := 2 z , equations (A.6) and (A.9) are transformed as follows, and they areassociated with the PDE operator L H for the Heston model, whose parameters are the mean reversion rateof κ , the long-run variance 2 θ , the volatility of variance √ σ , and the correlation √ ρ between stock priceand its variance. L H ˜ P s, ( t, x, ξ ) = 0 , L H ˜ P s, ( t, x, ξ ) = 12 W ξx∂ x (cid:0) x ∂ xx (cid:1) ˜ P s, ( t, x, ξ ) , (A.10)where L H = ∂ t + rx∂ x + κ (2 θ − ξ ) ∂ ξ + 12 x ξ∂ xx + 12 (cid:16) √ σ (cid:17) ξ∂ ξξ + (cid:18) √ ρ (cid:19) (cid:16) √ σ (cid:17) xξ∂ xξ − r. Similar to Fouque and Lorig (2011), by utilizing the feasibility of the Heston model, we can achieve thefollowing solutions of the abovementioned PDEs:˜ P s, ( t, x, ξ ) = e − rτ π (cid:90) e − ikq ˆ G ( τ, k, ξ ) ˆ h ( k ) dk, ˜ P s, ( t, x, ξ ) = e − rτ π (cid:90) e − ikq b ( k ) (cid:16) κθ ˆ f ( τ, k ) + 2 z ˆ f ( τ, k ) (cid:17) ˆ G ( τ, k, ξ ) ˆ h ( k ) dk, where τ = T − t , q = rτ + log x , b ( k ) = − W (cid:0) ik + k (cid:1) , ˆ f ( τ, k ) = (cid:90) τ ˆ f ( t, k ) dt, ˆ f ( τ, k ) = (cid:90) τ (cid:18) g ( k ) e sd ( k ) − g ( k ) e τd ( k ) − (cid:19) e d ( k )( τ − s ) ds, G ( τ, k, z ), ˆ h ( k ), d ( k ), and g ( k ) are defined in Theorem 1. To compute ˜ P s, and ˜ P s, , numericalintegrations need to be associated with a single integration and a triple integration. However, as in theforegoing discussion in 3.1, numerical methods for the triple integration are too computationally intensive.Fortunately, we can further simplify ˆ f and ˆ f , because the right-hand side of (A.10) is linear with respectto ξ . The induction process is given in detail as follows:ˆ f ( τ, k ) = (cid:90) τ (cid:0) g ( k ) e sd ( k ) − (cid:1) (cid:0) g ( k ) e τd ( k ) − (cid:1) e d ( k )( τ − s ) ds = e τd ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) (cid:90) τ (cid:0) g ( k ) e sd ( k ) − (cid:1) e sd ( k ) ds = (cid:32) e τd ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) (cid:33) (cid:18) d ( k ) (cid:16) g ( k ) e τd ( k ) − τ d ( k ) g ( k ) − e − τd ( k ) + 1 − g ( k ) (cid:17)(cid:19) + 1= e τd ( k ) (cid:16) g ( k ) (cid:0) e τd ( k ) − (cid:1) − τ d ( k ) g ( k ) + 1 (cid:17) − d ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) . If I ( τ, k ) := e τd ( k ) ( g ( k ) e τd ( k ) − ) , ˆ f ( τ, k ) = I ( t, k ) (cid:82) t I ( s,k ) ds from the second line of the above equations. Thus,we can obtainˆ f ( τ, k ) = (cid:90) τ I ( t, k ) (cid:90) t I ( s, k ) dsdt = (cid:90) τ I ( s, k ) (cid:90) τs I ( t, k ) dtds = (cid:90) τ (cid:0) g ( k ) e sd ( k ) − (cid:1) e sd ( k ) (cid:32) d ( k ) g ( k ) − d ( k ) g ( k ) e τd ( k ) − d ( k ) g ( k ) − d ( k ) g ( k ) e sd ( k ) (cid:33) ds = 2 τ d ( k ) g ( k ) + g ( k ) − d ( k ) g ( k ) (cid:0) g ( k ) e τd ( k ) − (cid:1) + τ d ( k ) g ( k ) − g ( k ) − d ( k ) g ( k ) . Appendix B. Solution for Poisson equation (A.8)
By the spectral theory, solution φ for Poisson equation (A.8) can be found as follows: φ ( y, z ) = z − y. Now, we briefly explain the way to find φ . It is known that the operator L s, has the eigenvalues λ n = − n ( n ∈ N ) and the family of eigenfunctions ψ n associated with eigenvalue λ n (i.e., L s, ψ n = λ n ψ n ) is ψ n ( y ) = (cid:115) n !Γ ( γ )Γ ( n + γ ) L n (cid:16) yν (cid:17) , where γ = zν and L n is an n -order Legendre polynomial, that is, L n ( w ) = w − γ e w n ! d n dw n (cid:0) w n + γ − e − w (cid:1) . Itis also known that ψ n form a complete orthogonal basis of the Hilbert space L (Φ). Φ is the invariantdistribution of Y t , which is defined in the foregoing section. Thus, y − z ∈ L (Φ) can be expressed as y − z = (cid:88) n ≥ c n ψ n ( y ) , (B.1)14here c n = (cid:104) ( y − z ) ψ n (cid:105) . For any n , c n are explicitly calculated as follows: c n = (cid:40) − ν √ z if n = 10 otherwise (B.2)We provide the induction process for the above formula (B.2). First, c can be easily obtained as follows: c = (cid:104) ( y − z ) ψ (cid:105) = (cid:104) y − z (cid:105) = 0 . In addition, c is computed as follows: c = (cid:104) ( y − z ) ψ (cid:105) = (cid:115) Γ ( γ )Γ (1 + γ ) (cid:90) ∞ ( y − z ) (cid:16) γ − yν (cid:17) Φ ( y ) dy = (cid:114) γ (cid:90) ∞ (cid:18) − ν y + (cid:16) zν + γ (cid:17) y − γz (cid:19) Φ ( y ) dy = − ν √ z. It is also proved that c n for n ≥ (cid:90) ∞ y (cid:20)(cid:16) yν (cid:17) − γ e yν d n dy n (cid:18)(cid:16) yν (cid:17) n + γ − e − yν (cid:19)(cid:21) (cid:104) y γ − e − ( y/ν ) (cid:105) dy = ν γ − (cid:90) ∞ y d n dy n (cid:18)(cid:16) yν (cid:17) n + γ − e − yν (cid:19) dy = − ν γ − (cid:90) ∞ d n − dy n − (cid:18)(cid:16) yν (cid:17) n + γ − e − yν (cid:19) dy = − ν γ − (cid:20) d n − dy n − (cid:18)(cid:16) yν (cid:17) n + γ − e − yν (cid:19)(cid:21) ∞ = ν γ − (cid:104) ( n + γ −
1) th degree polynomial without constant term × e − yν (cid:105) ∞ = 0 . Therefore, we eventually obtain the following φ by (B.1) and (B.2): φ ( y, z ) = c ψ ( y ) λ = zγ (cid:16) γ − yν (cid:17) = z − y. pp e nd i x C . E rr o r s i n t h ee m p i r i c a l t e s t τ < . . ≤ τ < . . ≤ τ < . τ ≥ . t o t a l H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ hS P X m e a n . . % . . . % . . . % . . % . . . % s t d . . . . . . . . . . H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h V I X m e a n . . . % . . . % . . . % . . . % . . . % s t d . . . . . . . . . . ( a ) I n - s a m p l ee rr o r s f o r p t i o n s d a t a τ < . . ≤ τ < . . ≤ τ < . τ ≥ . t o t a l H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ hS P X m e a n . . . % . . . % . . . % . . . % . . . % s t d . . . . . . . . . . H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h H e s t o n o u r s o/ h V I X m e a n . . . % . . . % . . . % . . . % . . . % s t d . . . . . . . . . . ( b ) O u t - o f - s a m p l ee rr o r s f o r p t i o n s d a t a T a b l e C . : E rr o r s i n e m p i r i c a l t e s t b a s e d o n p t i o n s d a t a . T h e r e s u l t s a r e s h o w n s e p a r a t e l y b y t i m e t o m a t u r i t y . T h e a bb r e v i a t i o n “o/ h ” m e a n s t h e v a l u e o b t a i n e db y d i v i d i n g t h ee rr o r f o r o u r m o d e l b y t h ec o rr e s p o nd i n g e rr o r f o r t h e H e s t o n m o d e l.l.