Joint calibration to SPX and VIX Derivative Markets with Composite Change of Time Models
Liexin Cheng Xue Cheng Xianhua Peng


3D Guidance Law for Maximal Coverage and Target Enclosing with Inherent Safety
    Praveen Kumar Ranjan, Abhinav Sinha, Member, IEEE, Yongcan Cao, Senior Member, IEEE
P. K. Ranjan and Y. Cao are with the Unmanned Systems Lab, Department of Electrical and Computer Engineering, The University of Texas at San Antonio, TX- 78249 USA. Emails: praveen.ranjan@my.utsa.edu, yongcan.cao@utsa.edu.A. Sinha is with the GALACxIS Lab, Department of Aerospace Engineering and Engineering Mechanics, University of Cincinnati, OH- 45221 USA. Email: abhinav.sinha@uc.edu.

    
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


   The Chicago Board Options Exchange Volatility Index (VIX) is calculated from SPX options and derivatives of VIX are also traded in market, which leads to the so-called “consistent modeling" problem. This paper proposes a time-changed Lévy model for log price with a composite change of time structure to capture both features of the implied SPX volatility and the implied volatility of volatility. Consistent modeling is achieved naturally via flexible choices of  jumps and leverage effects, as well as the composition of time changes. Many celebrated models are covered as special cases. From this model, we derive an explicit form of the characteristic function for the asset price (SPX) and the pricing formula for European options as well as VIX options. The empirical results indicate great competence of the proposed model in the problem of joint calibration of the SPX/VIX Markets.
    
     
Keywords: Time change; Lévy process; Option pricing; Consistent Modeling


















§ MOTIVATION AND FORMULATION



 §.§ Consistent Modeling Problem

By the definition from CBOE, Volatility Index (VIX), as an indicator of implied volatility in the following 30 days, is given as 

    VIX_t^2  = - 2/τ E^ℚ[ lne^-r τS_t+τ/S_t | ℱ_t ],

where τ = 30/365 and ℚ is the risk neutral measure of the equity market. If continuity of the asset price process is assumed, then equivalently

    VIX_t^2     = 1/τ E^ℚ[[ln S]_t+τ - [ln S]_t  | ℱ_t]
        =: 1/τ E^ℚ[∫_t^t+τ v^VIX_sds | ℱ_t],

where v^VIX is the squared volatility of S. Even though the asset price is not continuous, formula (<ref>) only leads to a third-order error O((d S_t/S_t-)^3), as shown in <cit.>. Hence the analysis below is still effective for general jump models to a large degree.

Likewise, we have the formula of VVIX, the volatility of volatility computed from VIX market:

    VVIX_t^2     =  - 2/τ E^ℚ[lne^-r τVIX_t+τ/VIX_t | ℱ_t]
       = 1/τ E^ℚ[[lnVIX]_t+τ - [lnVIX]_t  | ℱ_t]
        =: 1/4τ E^ℚ[∫_t^t+τ v^VVIX_sds | ℱ_t],

where v^VVIX is the variance of lnVIX^2. It is important to note that the first equality holds if we believe that measure ℚ is risk-neutral in both SPX option market and the VIX market. That is, the two markets can be consistently modeled.

The problem of joint calibration for SPX market and VIX market is equivalent to (or at least incorporates) the calibration of current VIX and VVIX, which can be approximated by √(v^VIX) and √(v^VVIX)/2, volatility and volatility of volatility respectively, if we consider a Markov setup and that τ is small. More generally, we have VIX_t^2 = g^M_1(v^VIX_t) and VVIX_t^2 = g^M_2(v^VVIX_t) by equation (<ref>) and (<ref>), where g_1^M(·), g_2^M(·) are functions determined by model parameters. Therefore, it is crucial to study the relationship of v^VIX and v^VVIX in the problem of consistent modeling. 

In the past literature on consistent modeling, there have been a lot of research on such volatility relationship. One line of work is aimed at reconstructing the widely used stochastic volatility models (SVMs) by allowing for more realistic and flexible vol-of-vol functionals. One such example is the 3/2 model in <cit.> and in <cit.>. <cit.> proposed a Heston vol-of-vol model, where the vol-of-vol consists of additional stochastic factors. Additionally, following the work of <cit.>, authors including <cit.>, <cit.> and <cit.> proposed rough volatility models in the joint calibration of stock and volatility smiles. And according to <cit.> and <cit.>, the role of jumps were studied and highlighted in the consistent modeling. Moreover, <cit.> studied the consistency condition of recovering SVMs from market models of the VIX futures
term structure. Another category of models suggest a multi-factor specification of volatility. The first attempt was made by <cit.>
 and <cit.>, who adopted a continuous diffusion model with double mean reverting structure. Multifactor affine specification was considered in <cit.>, <cit.> and <cit.>. And <cit.> considered a regime-switching Heston model, where sharp volatility regime shifts captures both volatility skews. Finally, there is some recent research that characterizes the volatility relationship using a non-parametric framework. <cit.> introduced a time-continuous formulation of the joint calibration problem, followed by <cit.>, who built a non-parametric discrete-time model that achieved exact joint calibration.

While many models above achieves satisfactory consistent modeling results, our approach, apart from having good joint calibration performance, is capable of theoretically decoupling smiles from the two markets. Such nice interpretation of our model is achieved via time change technique by representing the vol-of-vol v^VVIX in a linear mixture form of v^VIX and another free factor. We present the following examples to show how restrictive or implicit the volatility relationship is in certain SVMs.
[Heston Model]
    
    {[       dS_t/ S_t =r dt + √(v)_t dW_t,                                     ; d v_t = κ (θ - v_t)dt + η√(v)_t dZ_t;                                      ].

with E[W_tZ_t] = ρ t. We may compute the volatility of volatility under the assumption that v_t approximates VIX_t:

    v^VVIX_t ≈d[ln v]_t/dt = η^2/v_t = η^2/v_t^VIX.

Such an inverse relationship is unrealistic for VIX and VVIX, and therefore explains the unfavorable calibration result of Heston model. In fact, empirical results show that Heston models generates downward-sloping volatility smiles in VIX market as opposed to the upward-sloping observed smile, as shown in <cit.> and <cit.>. The model is poorly fitted under consistent modeling even if jump structures are incorporated, see <cit.>.

[3/2 Model]
        
    {[            dS_t/ S_t =r dt + √(v)_t dW_t,                                          ; d v_t = κ v_t(θ - v_t)dt + η v_t^3/2 dZ_t;                                           ].

with E[W_tZ_t] = ρ t. By the same reasoning, we obtain

    v^VVIX_t ≈d[ln v]_t/dt = η^2 v_t = η^2 v_t^VIX,

which is more close to empirical data, see figure <ref>. In fact, the correlation between VIX and VVIX is high in the past ten years (2013-2023), usually around 0.7. 

Moreover, 3/2 model generates upward sloping volatility smiles and captures the behavior of VIX better than a variety of stochastic volatility models in <cit.>. 

Despite its success in volatility market, the 3/2 model is restrictive in the situation of consistent modeling in the sense that v^VVIX is directly determined by v^VIX.

 

[Multi-factor Heston Model]
In the multi-factor Heston specification for volatility process, we have, by a reasoning similar to the Heston model,
    
    v_t^VVIX≈∑_i=1^nη_i^2v_t^(i)/(∑_i=1^nv_t^(i))^2,

where v^(i),   i=1,…, n are volatility factors and η_i,  i=1, …, n the corresponding coefficients. However, it is implicit how each factor v^(i) acts upon the volatility of volatility.


[Heston vol-of-vol; <cit.>]
    
    {[ d S_t / S_t=r  d t+√(v_t) S_t  d W_t,;  d v_t=κ(θ-v_t) d t+η_t √(v_t) d B_t,;                    d W_t  d B_t=ρ d t ].

where η_t=η(Y_t^ε, Z_t^δ) is a stochastic factor correlated with W and B. Then

    v_t^VVIX≈d[ln v]_t/dt = η_t^2/v_t^VIX.

Although the randomness of η improves the calibration of v^VVIX, η directly depends on v^VIX and hence is not flexible enough.







 §.§ Decoupling Smiles via Composite Time Change Approach


To study the volatility relationship in composite time change (CTC) models . We first introduce some background knowledge on time change approaches. Time changes can be interpreted as the intensity of business activities that drives the variation in volatility of an asset. The original clock {t, t ≥ 0} is referred to as calendar time and time change process T is called business time. The log-price process of an asset is originally thought to be stationary and of independent increment, i.e. a Lévy process L. Every time a market event happens and drives the variation in volatility, the change is reflected in the time change, either via accelerating or slowing the business clock. And the real market price of the asset is updated under the business clock, namely L_T. 

Originally, <cit.> and <cit.> proposed a subordinated Brownian motion model for log price. The time change introduces jumps in volatility. <cit.> and <cit.> introduced time-changed Lévy models, where the time change is absolutely continuous, through which
stochastic volatility is introduced. <cit.>, <cit.> and <cit.> modeled the dependence of multi-assets via correlation of
subordinated Brownian motions. <cit.> extented the time-chanegd Lévy model by considering the
combination of time change and a certain type of composite time change. And recently, <cit.> established a unified TCLM structure that allows leverage via diffusion as well as jumps.

But like SVMs, TCMs with specifications shown in <cit.> are restricted by the relationship v^VVIX_t =f(v^VIX_t) for a model-determined function f. But when extended to composite time change, the model may generate the form of vol-of-vol as

    v^VVIX_t = av^VIX_t + bv^I_t,

where v^I is the idiosyncratic component of v^VVIX and a, b are determined by model parameters, which are typically considered steady in short periods. The linear combination satisfies both the need for the dependence of VVIX on VIX and the flexibility of VVIX. v^VIX is calibrated to the SPX market while a free factor v^I is calibrated to the VIX market. 








    A composite time change has the form T^co = T_V,
    where T = {T_t, t≥ 0} and V = {V_t, t≥ 0} are time changes, i.e., {T_t}, {V_t} are non-decreasing and right-continuous, and T_t, V_t are stopping times for every t ≥ 0.



Here we assume time changes to be absolutely continuous and the base Lévy process to be a standard Brownian motion. Specifically, dT_t = u_tdt and dV_t = v_tdt, where u and v are two independent Itô processes. Then the instantaneous variance of the model is v_t^VIX = d T_V_t/dt = u_V_tv_t. Then the product rule

    dv_t^VIX = d(u_V_tv_t) = u_V_tdv_t + v_tdu_V_t +d[u_V_·, v_·]_t

gives

    v_t^VVIX≈d[ln v^VIX]_t/dt =1/v_t^2d[v]_t/dt + 1/(u_V_t)^2d[u_V]_t/dt.

Equation (<ref>) results in different linear mixture forms according to the specification of activity rate process. And the ideal form of equation (<ref>) can be achieved by the composite 3/2 model proposed in section 3.2. In addition, the model parameters in time change T and V are naturally separated and linearly combined in form.



 §.§ Organization of the Paper

The article is organized as follows. In section2, we summarize the theory of time-changed Lévy models developed by <cit.> and <cit.> and we show how the technique of leverage neutral measure change helps form the characteristic function; section 3 develops the theory of composite time change models, where a general form is considered and characteristic fuctions derived. We also discusses some useful specifications of CTC models. In Section 4, we introduce the application of the model in derivative pricing, including European options and VIX options. We show that the European option pricing in CTC models can be conducted quite efficiently. In section 5, we perform real-market joint calibration and discusses the results. And the last section concludes.



§ PRELIMINARIES: TIME-CHANGED LÉVY PROCESSES



 §.§ Lévy processes


    A Lévy process, L(t), on a filtered probability space (Ω, 𝔽,{𝔽_t}_t ≥ 0, ℙ) is a continuous-time process with independent and stationary increments with a characteristic function ϕ_L(m ; t)=e^t Ψ_L(m), m ∈ℝ with characteristic exponent

    Ψ_L(m)=i α m-m^2/2σ^2+∫_ℝ(e^i m x-1-i m x 1_|x| ≤ 1) ν(d x),

where α∈ℝ, σ∈ℝ^+ and ν is a positive measure on ℝ such that v({0})=0,   ∫_ℝ(|x|^2 ∧ 1) v(d x)<∞.
The triplet (α, σ^2, v(d x)) determines the Lévy process and is referred to as differential characteristics.


A subordinator T = {T_t, t≥ 0} is a Lévy process such that t ↦ T_t is non-decreasing.



 §.§ Time-changed Lévy Processes

Leverage neutral measure, first introduced in <cit.>, is a complex-valued measure change technique that enables the explicit computation of the characteristic function (Ch.f.) of time-changed Lévy process L_T, especially when L and T are not independent.

 The Lévy processes L considered in this paper are sufficiently integrable. That is, ∫_|x|>1 e^m xν(d x)<∞ for all m∈ℝ.




There are three versions of assumptions in the order of restrictiveness.

    
  * (Type 1) T satisfies sufficient regularity conditions. That is, the Laplace transforms considered in this paper always exist for every t > 0.


  * (Type 2) condition of type 1 and T is in synchronization with some adapted process X, i.e. X is constant a.s. on the interval [T_t-, T_t] for each t > 0.


  * (Type 3) condition of type 1 and the time change is absolutely continuous with T_t = ∫_0^t v_s ds,   v_t > 0 a.s. or an independent Lévy subordinator.




    Let the Lévy process L satisfy the assumption and time change T be of type 1. The log-price is X_t = L_T_t. 

    1. If T is independent of L, then the Ch.f. of X_t is
    
    ϕ_X(m;t) := E[e^i m X_t] = ϕ_T(-iΨ_L(m); t)

     
2. Otherwise,
    
    ϕ_X(m; t) := E[e^i m X_t]=E^ℚ[e^T_t Ψ_L(m)] = ϕ_T^ℚ(-iΨ_L(m); t),

where E[·] and E^ℚ[·] denote expectations under measures ℙ and ℚ, respectively. The new class of complex-valued measures ℚ(m) is absolutely continuous with respect to ℙ and is defined by

    .dℚ(m)/dℙ|_ℱ_t = M_t(m),

with

    M_t(m) := exp(i m L_t+t Ψ_L(m)),    m ∈ℝ.

Moreover, the Ch.e. of L under the leverage neutral measure is given by

    Ψ_L^ℚ(z)=Ψ_L(z+m)-Ψ_L(m),    m, z ∈ℝ.



    Since the result above is a bit different from those proposed in <cit.> and <cit.>, we leave the proof in Appendix A.


Extra conditions for T is required if we want to characterise the time-changed process X. Specifically, if T is of type 2, then semimartingale X has local characteristics 
    (α dT_t-, σ^2 dT_t-, ν(dx)dT_t-),

where (α, σ^2, ν(dx)) is the Lévy characteristic of L (see <cit.>). Such an extra assumption makes the model identifiable, i.e., L and time change T are unique up to a constant multiplication.



§ COMPOSITE TIME CHANGE MODELS

In this section, we assume a composite time change model X_t = L_T_V_t, where both time changes T and V are absolutely continuous with activity rates u and v, respectively. Furthermore, we assume ℙ to be the risk-neutral measure and ℱ the original filtration under which the base Lévy process L is adapted and T, V are time changes. We denote by ℱ^X with ℱ^X_t := σ(ℱ_T_V_t) ∨σ(ℱ_V_t) the filtration generated by X, T_V and V. Denote by E_t[·] the expectation taken under filtration ℱ^X_t.



 §.§ Model Theory

Under such assumptions, the model is specified as follows

    {[                         d L(t) =-Ψ(-i) dt + σ dW(t) + η dJ(t),                                                               ; d u(t) = α_T(u(t))dt + β_T(u(t))dZ(T_t) + γ_T(u(t)) dJ_T(T_t),; d v(t) = α_V(v(t))dt + β_V(v(t))dB(V_t) + γ_V(v(t)) dJ_V(V_t),                                                                ].

where W,Z, B are (ℙ, ℱ)-Brownian motions, J, J_T, J_V are (ℙ, ℱ)-subordinators and Ψ(m)=-m^2/2σ^2+Ψ_J(η m). In addition, α_i(·), β_i(·), γ_i(·), i = T, V are unspecified functions. Finally, γ_i(·) ≥ 0,   i = T,V to guarantee the positivity of u and v.

To introduce leverage effect, we assume that [W, Z]_t = ρ_Tt, [W, B]_t = ρ_Vt but Z and B are independent. We also assume that J and J_T have joint distribution F_T, J and J_V have joint distribution F_V. J_T and J_V are independent.

Note that if we let α_i(·) be affine and β_i(·), γ_i(·) be constant, then u, v both have an affine structure.






Next, we introduce the leverage neutral measure that will be useful for the computation with composite time change. The leverage neutral measure ℚ(m) is defined by

    .dℚ(m)/dℙ|_ℱ_t = M^m_t,

where

    M^m_t = exp(imL_t + tΨ_L(m)).

is a (ℙ, ℱ)-exponential martingale.

Then we have the following result of the Ch.f. of X.

The Ch.f. of the composite time-changed process X is given by

    ϕ_X(m; t) : = E[e^imX_t] =E^ℚ[ϕ_T^ℚ(-i Ψ_L(m); V_t)]

 where, under the leverage neutral measure ℚ(m), u^ℚ and v^ℚ are given by

    
    du^ℚ(t)    = (α_T(u^ℚ(t)) + imρ_T σ v^ℚ(t)β_T(u^ℚ(t)))dt + β_T(u^ℚ(t))dZ^ℚ(T_t)+ γ_T(u^ℚ(t))  dJ_T^ℚ(V_t)
    
                dv^ℚ(t)    = (α_V(v^ℚ(t)) + imρ_V σ v^ℚ(t)β_V(v^ℚ(t)))dt + β_V(v^ℚ(t))dB^ℚ(V_t)+ γ_V(v^ℚ(t))dJ^ℚ_V(V_t),

where B^ℚ, Z^ℚ are (Q,ℱ)-Brownian motions. The Ch.e. of J_i^ℚ,  i=T,V are given by

    Ψ_J_i^ℚ(z) = Ψ_J, J_i(η m, z) - Ψ_J(η m),

where Ψ_J, J_i(m,z) is the Ch.e. of mJ + zJ_i.



    Note that the time changes T and V are still real-valued and positive under the new measure. The dynamics of activity rates above is only effective when Ch.f. is computed.



 By the independence of T and V,

    ϕ_X(m; t)    = E^ℚ[e^T_V_tΨ_L(m)]
       = E^ℚ[E^ℚ[e^T_V_tΨ_L(m)| V_t]]
        = E^ℚ[ϕ_T^ℚ(-i Ψ_L(m); V_t)].


Next, we derive the dynamics of T and V under the risk-neutral measure ℚ. By theorem 1, the Ch.e. of L under measure ℚ is given by 
    
    Ψ_L^ℚ(z) = Ψ_L(z+m) - Ψ_L(m).

    By the independence of Brownian motions and jump processes, the Ch.e. of W under measure ℚ is 
    
    Ψ_W^ℚ(z) = Ψ_W(z+σ m) - Ψ_W(σ m) = -1/2z^2 -σ mz.

Therefore W^ℚ(t) := W(t) - imσ t is a (ℚ, ℱ)-Brownian motion. And by the correlation [W, Z]_t = ρ_T t, Z^ℚ(t) := W(t) - imρ_Tσ t is a (ℚ, ℱ)-Brownian motion. B^ℚ is defined likewise. 

Next we show the Ch.f. of J_T and J_V. We note that the single joint distribution determines the joint distribution of J(t) and J_i(t),   i = T, V at all time t > 0. Under the leverage neutral measure, the Ch.f. of J_V becomes

    ϕ^ℚ_J_V(z;t)     = Eexp(izJ_V(t) + imη J(t)-tΨ_J(η m))
        = exp(t(Ψ_J, J_V(η m,z) - Ψ_J(η m)),

where Ψ_J, J^V(m,z) is the Ch.e. of mJ + zJ_V. It follows that 

    Ψ_J_V^ℚ(z) = Ψ_J, J_V(η m, z) - Ψ_J(η m).




When affine structure is imposed, that is, α_T(u(t)) = κ_T(θ_T - u(t)), α_V(v(t)) = κ_V(θ_V - v(t)) and β_i(·) ≡σ_i, γ_i(·) ≡η_i, i = T,V, then Ch.f. of X is explicit as follows:


    ϕ_X(m; t) =   1/π∫_0^∞∫_0^∞ f(s;m)Re[e^-z sϕ_V^ℚ(-iz;t)] d z_Id s,

where z = z_R + iz_I with z_R > 0 and 

    ϕ_V^ℚ(m;t) = e^b_V(t)v(0) +c_V(t),

with the affine exponents b_V(t), c_V(t) solutions to the system of Riccati-type ODEs

    b_V(t)^'   = im - κ^ℚ_V b_V(t) + σ_V^2/2b_V(t)^2 + Ψ_J_V^ℚ(iη_V b_V(t))
    
        c_V(t)^'    = κ_V θ_V b_V(t).

with b_V(0) = c_V(0) = 0, κ^ℚ_V = κ_V - imρ_Vσ_Vσ. The function f is given by

    f(t;m) = e^b_T(t)u(0) + c_T(t)

with coefficients satisfying 

    b_T(t)^'   = Ψ_L(m) - κ^ℚ_Tb_T(t) + σ_T^2/2b_T(t)^2 + Ψ_J_T^ℚ(iη_T b_T(t))
    
        c_T(t)^'    = κ_T θ_T b_T(t).

with b_T(0) = c_T(0) = 0, κ^ℚ_T = κ_T - imρ_Tσ_Tσ. The Ch.e. of J_i, i=T,V are given in equation (<ref>).



By the inverse Laplace transform,

    ϕ_X(m; t)      = E^ℚ[ϕ_T^ℚ(-i Ψ_L(m); V_t)]
       = 1/π∫_0^∞∫_0^∞ f(s;m)Re[e^-z sϕ_V^ℚ(-iz;t)] d z_Id s,

where z = z_R + iz_I with z_R > 0 and f(t) = ϕ_T^ℚ(-iΨ_L(m);t) (Here we denote f(t;m) by f(t) for simplicity). If an affine structure is imposed for T and V, then as is given in <cit.>, the Laplace transform of T_t is given by

    f(t) = e^b_T(t)u(0)+ c_T(t),

where coefficients b_T and c_T are given by the equation (<ref>). Likewise, the Ch.f. of ϕ_V has coefficients given by the equation (<ref>).







 §.§ Specifications

[Composite 3/2 Model]
        
    {[                        L(t) =-1/2 dt + dW(t),; d u_t = κ_T u_t(θ_T - u_t)dt + σ_T u_tdZ_T_t,;  d v_t = κ_V v_t(θ_T - v_t)dt + σ_V v_tdB_V_t ].

with E[W_tZ_t] = ρ_T t, E[W_t B_t] =ρ_Vt and E[Z_tB_t] = 0. 


In this specification, we have by equation (<ref>) that 

    v_t^VVIX≈d[ln v^VIX]_t/dt   = (σ_Vu_V_t)^2v_t^3 + (σ_Tv_t)^2 (u_V_t^3 v_t)/(v_t^VIX)^2
       = σ_V^2 v_t + σ_T^2 v_t^VIX,

a linear combination of factors as in equation (<ref>). 

We clearly see a separation of effects, including the effect of VIX and the idiosyncratic component. The SPX market calibrates v_t^VIX and implies a general relationship of model parameters σ_V and σ_T, while the VIX market calibrates v and determines the model parameters. In other specifications likewise, we also obtain a linear combination of factors.

To derive the Ch.f. under composite 3/2, we note that, according to <cit.>, the Laplace transform of V has

    E^ℚ(m)(e^-λ∫_t^τ v_s  d s| v_t) 
         = Γ(γ_V-α_V)/Γ(γ_V)(2/σ_V^2 y(t, v_t))^α_V M(α_V, γ_V, -2/σ_V^2 y(t, v_t)),

where

    y(t, v_t)     =v_t  e^κ_V θ_V(τ-t) - 1/κ_V θ_V, 
    α_V     =-(1/2-p_V/σ_V^2)+√((1/2-p_V/σ_V^2)^2+2 λ/σ_V^2)
    γ_V     =2(α+1-p_V/σ_V^2), 
    
    p_V     =-κ^ℚ_V :=-κ_V +i σ_V ρ_V m,

and M(α, γ, z) is the confluent hypergeometric function, defined as

    M(α, γ, z)=∑_n=0^∞(α)_n/(γ)_nz^n/n !,

and

    (x)_n=x(x+1)(x+2) ⋯(x+n-1) .

And 

    f(t;m)    = ϕ_T^ℚ(-iΨ_L(m);t)
        =Γ(γ_T-α_T)/Γ(γ_T)(2/σ_T^2 y_T(0, u_0))^α_T M(α_T, γ_T, -2/σ_T^2 y_T(0, u_0)))

where

    y_T(0, u_0)     =u_0 e^κ_T θ_Tt - 1/κ_T θ_T, 
    α_T    =-(1/2-p_T/σ_T^2)+√((1/2-p_T/σ_T^2)^2+2 q_T/σ_T^2)
    γ_T     =2(α_T+1-p_T/σ_T^2), 
    
    p_T     =-κ^ℚ_T :=-κ_T +i σ_Tρ_T m,
    
    q_T     = i m/2+m^2/2.


Despite its nice interpretation, the composite 3/2 can be time-consuming in pricing, particularly in VIX pricing. A modified version of composite 3/2 model is constructed:
    
    {[                        L(t) =-1/2 dt + dW(t),; d u_t = κ_T u_t(θ_T - u_t)dt + σ_T u_tdZ_T_t,;       d v_t = κ_V (θ_T - v_t)dt + σ_V dB_V_t. ].

That is, the second time change V is substituted by a CIR process. It's shown in the section of VIX pricing that such formation is efficient in the method of exact simulation.
[Composite 3/2 + Jump]
            
    {[            L(t) =-Ψ(-i) dt + dW(t) + η dJ(t),; d u_t = κ_T u_t(θ_T - u_t)dt + σ_T u_tdZ_T_t,;     d v_t = κ_V v_t(θ_T - v_t)dt + σ_V dB_V_t ].

with E[W_tZ_t] = ρ_T t, E[W_t B_t] =ρ_Vt, E[Z_tB_t] = 0 and J is a Lévy process. Under the specification, we only change q_T as

    q_T = -Ψ_L(m) = im + m^2/2 + im Ψ_J(-η i) - Ψ(η m).



[Composite Heston Model]
    
    {[                     L(t) =-1/2 dt + dW(t),                                           ; d u(t) = κ_T (θ_T - u(t))dt + σ_T dZ(T_t),; d v(t) = κ_V (θ_V - v(t)dt + σ_V d B(V_t). ].

where E[W_t Z_t] = ρ_T t, E[W_t B_t] =  ρ_V t and E[Z_t B_t] = 0.


Likewise, we have

    v_t^VVIX   ≈d[ln v^VIX]_t/dt
       = (σ_Vu_V_t)^2v_t + (σ_Tv_t)^2 (u_V_t v_t)/(v_t^VIX)^2
       = σ_V^2/v_t + σ_T^2/v_t^VIX

The composite version of Heston inherits the inverse relationship between VIX and VVIX, but a linear mixture effect is incorporated.


Brought up in <cit.>, the composite time change T^co = T_V is a pure-jump process. T and V are independent of L as well as of each other.

    𝔼[e^iu X_t]=𝔼[𝔼[exp(T_V_tΦ_L(u)) | V_t]]=𝔼[ϕ_T (-i Ψ_L(u); V_t)].

The model exhibits stochastic jump intensity even if L has no jumps.

However, the model is theoretically redundant in the sense that L_T is in fact a Lévy process under the model assumption. Nevertheless, the composite version of Lévy process is convenient in exhibiting flexible moments of distribution, which is usually meaningful in practice.    


[Leverage via Reflected Jumps]
            
    {[              L_t =-Ψ(-i) dt + η dJ(t),                                       ; d v_t = κ (θ - v_t)dt - η_J dJ^T(T_t), ].

where J^T = J^- is the negative part of the CGMY process J. 

The model has shown superior performance in <cit.>.


[Compound Poisson With General Leverage]
    
    {[                        L(t) =-Ψ(-i)t + η∑_i = 1^N_t J_i,                                                         ; d u(t) = κ_T(θ_T - u(t))dt + η_T d(∑_i = 1^N_T_t J^T_i),;                d v(t) = κ_V (θ_V - v(t))dt + σ_V dB(V_t) ].

where the jump sizes J and J^T are correlated with the joint distribution F(x,y).

Under the leverage neutral measure, the Ch.e. of Y_t:=∑_i=1^N_tJ^T_i becomes

    Ψ^ℚ_Y(z) = λ (E exp(izJ_1 + iuJ^T_1) - ϕ_J(u)) = λϕ_J(u)(ϕ^u_J^T(z) - 1),

where 
    ϕ_J^T^u(z) = ∫ e^izy(e^iux-Ψ_J(u)dF(x,y)).

It can be interpreted as a new (complex-valued) compound Poisson process with jump intensity λΨ_J(u) and jump size with a tilted distribution ϕ^u_J^T(z).

[Composite Jump Heston]
      
    {[                     L(t) =-Ψ(-i)t + η J(t),                                            ; d u(t) = κ_T(θ_T - u(t))dt - η_T dJ(T_t)^-,;   d v(t) = κ_V(θ_V - v(t))dt + η_V dZ(V_t), ].

where J(t)^- is the negative component of the CGMY processes J(t), and Z is a Brownian motion. 

In this specification, leverage is introduced purely by simultaneous jumps of return and volatility. <cit.> demonstrate a superior performance of a single time change JH model over other classic models, e.g. Heston, BNS.





[Multifactor]
    It has been shown in empirical study that risks reflected in diffusion and jumps are of different sources. It is therefore natural to consider a TCL of the form 
    
    X = B_T + L_V.



    




§ APPLICATION IN DERIVATIVES PRICING



 §.§ CTC-COS Method

Since we've obtained the expression for the Ch.f. of log price process X, the pricing of European options follows directly. For example, readers may consider acceleration methods such as FFT ( <cit.>) or COS method (<cit.>), to efficiently compute option prices. 

In our CTC model, we show how the Ch.f. and option prices can be efficiently computed with a CTC-COS method as follows. 

Given current time t and expiry date s, the price of a European call option with strike K is numerically approximated by

    C(K, τ) ≈ e^-r τ2/c.∫_0^c  (∑_k=0^N-1^'Re{f(y; kπ/b-a) A_k} V_k)(  ∑_l=0^M - 1^'Re{ϕ_V^ℚ(lπ/c; s)}cos (lπ y/c))dy.

where a,b,c are integration range and f is given in (<ref>).


























According to the cosine expansion method,

    C(K, τ)    ≈ e^-r τ∑_k=0^N-1^'Re{ϕ_X(k π/b-a ; s) A_k}V_k
        =  e^-r τ∑_k=0^N-1^'Re{E^ℚ[f(V_s; k π/b-a))]A_k} V_k
        = e^-r τ∑_k=0^N-1^'∑_l=0^M - 1^'Re{ϕ_V^ℚ(lπ/c; s)}Re{U_kl A_k}V_k
       ≈ e^-r τ2/c.∫_0^c  (∑_k=0^N-1^'Re{f(y; kπ/b-a) A_k} V_k)(  ∑_l=0^M - 1^'Re{ϕ_V^ℚ(lπ/c; s)}cos (lπ y/c))dy.


    U_kl = 2/c∫_0^c f(y; kπ/b-a)cos(lπ y/c)dy,



    A_k = exp{-i k πa/b-a+ikπ(ln S_t/K + rτ)/b-a}

and

    V_k = 2/b-a∫_0^b K(e^y - 1)cos(k πy-a/b-a)dy.



What we do in the transformation above is performing double cosine series expansions and then re-ordering the summation and integration. Since the double summations are separate, the overall complexity is O(ND) under affine specifications, where D is the discretization degree of numerical integration. Such computation cost is only slightly higher than those of one-factor models with affine structure, e.g. Heston model.

In practice, when pricing the whole volatility surface, the number of strike prices does not add computational complexity since fourier-based methods like COS allows for a separation of strikes and the underlying asset. Meanwhile, the temporal discretization in the solving process of ODEs in function f and g, enabling the computation on all maturities at once. Therefore, we may price the whole volatility surface as efficently as in the case of a single option.

It's surprising that the CTC models have the same order of complexity as the single TCMs. Table <ref> shows the speed of pricing the volatility surface (containing 4316 option quotes on Janurary, 9, 2023), where the composite Jump Heston and composition Heston models finish the pricing within 15s and 10s respectively.















 §.§ VIX Options





















To price VIX derivatives, it's common practice to first derive the relationship between the VIX and spot variance at maturity since the distribution of VIX is not directly available. It's known that VIX^2 is linearly dependent on the spot variance for the Heston model. However, for most other models, such relationship is implicit and requires simulation or a numerical procedure to derive.


[Rough Bergomi]
    As shown in <cit.> for rough Bergomi models, the VIX_T^2 is expressed as an integral form and a computationally costly simulation is needed to obtain samples of VIX_T.


[3/2 Model]
    The VIX-Spot relationship is also implicit for a 3/2 model. A numerical differentiation is needed and leads to additional computational cost in numerical pricing.


Next, we will show that both single time change models and CTC models have simple VIX-Spot forms if the time changes have affine activity rates. Non-affine VIX-Spot relationship are also discussed.



  §.§.§ Single Time Change Models




When a single time change T = ∫_0^· v_sds is considered with an affine process v,

    VIX_t^2 = av_t + b
 with

    a = 2 (Ψ(-i) - η E J(1))ϕ(m)


    b =  2 (Ψ(-i) - η E J(1))κθ/m (ϕ(m) - 1)

and ϕ(x;τ) = e^xτ - 1/x τ and we denote ϕ(x) for short if there is no confusion.




    VIX_t^2     = - 2/τ E_t[lne^-r τS_t+τ/S_t]
        = - 2/τ E_t[X_t+ τ - X_t ] 
        =2/τ( Ψ(-i) E_t[ T_t+τ - T_t] - E_t[η(J(T_t+τ) - J(T_t))])
        = 2 (Ψ(-i) - η E J(1))g(v_t ; τ),

where τ = 30 / 365 and g(v_t;τ) =  1/τE_t[ T_t+τ - T_t]. Since

    Ev_t = v_0 + κθ t + m ∫_0^t Ev_sds,

where m = η_J EJ^T_1 - κ < 0. Differential form gives

    g^'_t = κθ + m g_t,

with initial condition g_0 = 0. The solution is

    E[v_t|v_s] = v_s e^m (t- s) + κθ/m (e^m (t-s) - 1).

Then we have

    g(v_t;τ) = v_t ϕ(m) + κθ/m (ϕ(m) - 1).



Since the Ch.f. of v is easily obtained by solving an ODE system, the pricing formula of VIX options is given as:


    C^V(K, τ)=    e^-rτ/2 a √(π)
       ×∫_0^∞Re[e^φ  b / aϕ_v(φ; s) 1- erf(K √(φ / a))/(√(φ  / a))^3] dφ_I

where φ=φ_R+φ_I i is a complex variable with ϕ_R>0. ϕ_v(φ; s) is the characteristic function of v_s and a, b are given in Proposition <ref>.




  §.§.§ CTC Models

For the case of affine CTC models, the VIX-spot relationship is still explicit and easy to obtain.

For a CTC model with affine activity rates,

    VIX_t^2 = A(v_t) u_V_t + B v_t + C(v_t),

where coefficients 

    A(v_t) = Mϕ_V(-im_T;t,t+τ)-1/m_Tτ,


    B = M κ_Tθ_Tϕ(m_V)/m_T

and 

    C(v_t) = M (- κ_Tκ_Vθ_Tθ_V(ϕ(m_V) - 1)/m_Tm_V+ κ_Tθ_T/m_T^2τ(ϕ_V(-im_T;t,t+τ) - 1))

with a common multiplier

    M = 2 (Ψ(-i) - η E J(1)).




    
    VIX_t^2     = 2 (Ψ(-i) - η E J(1))/τE_t [T_V_t+τ - T_V_t] 
       =  2 (Ψ(-i) - η E J(1))/τE_t [E[T_V_t+τ - T_V_t|ℱ_t, V_t+τ]]
        =2 (Ψ(-i) - η E J(1))/τE_t[ u_V_tΔ_V(τ)ϕ(m_T; Δ_V(τ)) + κ_T θ_T/m_TΔ_V(τ)(ϕ(m_T; Δ_V(τ)) - 1)]
        =2 (Ψ(-i) - η E J(1)){(m_T u_V_t + κ_T θ_T)(E_t e^m_TΔ_V(τ) -1)/(m_T)^2τ.
       .       - κ_Tθ_T/m_T(v_tϕ(m_V) + κ_Vθ_V/m_V (ϕ(m_V) - 1))}
        = 2 (Ψ(-i) - η E J(1)){ϕ_V(-im_T;t,t+τ)-1/m_Tτu_V_t - κ_Tθ_Tϕ(m_V)/m_T v_t + Const},

where Δ_V(τ) = V_t+τ - V_t and constant

    Const = - κ_Tκ_Vθ_Tθ_V(ϕ(m_V) - 1)/m_Tm_V+ κ_Tθ_T/m_T^2τ(ϕ_V(-im_T;t,t+τ) - 1).



For other specification of time changes, there generally does not exist explicit expression for VIX_t. Still, we could recover it from the Laplace transform. That is,

    VIX_t = 2 (Ψ(-i) - η E J(1))/τg(v_t, u_V_t, τ),

where 
    g(v_t, u_V_t, τ) := -.∂/∂ lE[exp{-l(T_V_t+τ - T_V_t)}| u_V_t, v_t]|_l=0


The Laplace transform can be explicitly computed via

    E[exp{-l(T_V_t+τ - T_V_t)}| u_V_t, v_t] 
       ≈∑_k=0^M-1Re{ϕ_V(kπ/c; τ)}(2/c∫_0^c ℒ^T(l; y)cos(kπ y/c)dy)
        = 2/c∫_0^c ℒ^T(l; y) ∑_k=0^M-1(Re{ϕ_V(kπ/c; τ)}cos(kπ y/c)) dy

where ϕ_V and ℒ^T are computed conditional on u_V_t, v_t. 

In practice, the formula above is computed based on a large amount of simuated v_t and u_V_t and the whole time cost is high. But in the special case where an affine v is considered, but not necessarily u, the Ch.f. ϕ_V is exponentially-affine with respect to v_t and, based on such explicit linear dependence, the computation is significantly faster.



  §.§.§ Exact Simulation of Spot Variances

Despite its explicit VIX-Spot relationship, a simulation procedure for u_V_s and v_s is still needed. While Monte Carlo method is commonly slow in practice, we argue that there exists exact simulation methods for certain models. That is, we only need to simulation the distributions at maturity as opposed to the whole trajectories. As a result, there is no discretization error and the simulation procedure can be quite fast.

We assume that the Ch.f. of u is available and v follows a CIR process, our main concern is to simulate v_T and u_V_T. Then the sample of VIX is obtained from formula (<ref>). 

    
  Step 1 Simulate v_T from a non-central chi-square distribution.
    
  Step 2 Simulate the conditional distribution (V_T | v_0, v_T) by the method of <cit.>.
    
    By the independence of u and v, we may simulate V_T first and then simulate u_V_T with a non-central chi-square distribution. Therefore, the conditional distribution of V_T, in the form of (∫_0^T v_s ds | v_0, v_T), needs to be derived. This is the exactly same problem faced in the Heston simulation, as brought up in <cit.>. We then apply the method of gamma expansion in <cit.>, where the conditional distribution is efficiently simulated as a sum of independent variables.
    
  Step 3 Simulate u_V_T by inverting the Ch.f. of u at V_T. 

    The Ch.f. of u is explicit for some models, e.g. Heston  and 3/2 models given in <cit.>. For composite JH model (<ref>), we can ease the computation by precomputing the Ch.f. of u_t for a range of t ∈ [0, maxV_T] along the numerical discretization of the corresponding ODE. The sample of u_V_T is then obtained by matching the precomputed values with V_T samples.
    
  Step 4 Obtain VIX_T by formula (<ref>) and compute C^V(K, T) by taking the mean of the payoff.










To sum up, we price European options by






§ JOINT CALIBRATION



 §.§ Data

The data we consider spans two years, containing the implied volatility surfaces of the S&P 500 and the VIX from January 3, 2023, to February 28, 2023. The implied volatility is calculated using futures prices that are inferred from the put-call parity with ATM options.

Following the procedure of <cit.>, we clean the data by removing options with time to maturity less than 7 days or more than one year. Options quotes with negative bid-ask spreads, zero bid price, zero traded volume / open interest are also removed. For SPX options, we further exclude data with moneyness below 0.5 or above 1.4. We also remove all the ITM options and recover the corresponding futures prices using the ATM put-call parity. For VIX options, we exclude data with moneyness below 0.4 or above 3.3 and all the put options. If a VIX ITM call is illiquid, we use the put-call parity to infer the liquid price of the call from a more liquid VIX OTM put (<cit.>).

The final sample is made of 162,438 SPX options (daily average 4,165) and 10,254 VIX options (daily average 263).


We calibrate on a sample by solving the optimization problem below.

    Θ=  Θmin1/N_S∑_i=1^N_S(σ_S^i(Θ)-σ_S^i)^2+1/N_V∑_i=1^N_V(σ_V^i(Θ)-σ_V^i)^2,

where N_S and N_V are the number of SPX and VIX options in the sample.



 §.§ Calibration Procedure

Calibrate the model based on the data on some specific dates. For example, we may choose the dates with remarkably high or low VIX / VVIX to test the robustness of the model in different market scenarios.

Since the calibration of VIX market is more challenging, we typically calibrate to VIX option data first, and then take the initial set of parameters in joint calibration based on the result in VIX calibration.


























































 §.§ Calibration Results

The next table is the calibration result of SPX market. The RMSE in the table below is the root mean squared error of implied volatilities.












 §.§ Joint Calibration

First, we test the performance of joint calibration on date 2023.01.13 (with VIX = 18.35) and date 2020.04.24 (with VIX = 35.93). We specifically consider the JH model, the composite JH model and the 2-factor JH model.



To make comparison with other benchmark models, we compare the performances with the results in <cit.> on 2008.10.22 and 2012.05.16.




It's shown in the results that, compared with 2-factor JH, composite JH models generally performs better. This demonstrates the effectiveness of time composition v.s. time combination. In particular, when the market volatility is large, composite JH models become even better.









§ CONCLUSION

Specifically, the contributions of our work are three-folded.

Firstly, we develop a generalized form of composite time-changed Lévy models. These models have the advantage of exhibiting various variance, skewness and kurtosis. Moreover, as we show in detail in section 4, CTC models typically have good tractability. It only requires a complexity of O(N) if the affine structure of time changes is imposed.

Secondly, we theoretically demonstrate the effectiveness of our model in consistent modeling. Unlike historical works of consistent modeling, the composite time change models provide explicit interpretation in its decoupling mechanism of volatility and volatility of volatility. 

And finally, we validates the superiority of our proposed model in the consistent modeling problem. We develop its option pricing theory and test its performance in the joint calibration problem of SPX and VIX option market. As shown in section 5, the composite time change models successfully calibrate the joint smiles in real market.




§ PROOF OF THEOREM 1

    
    
        Consider a filtered space (Ω, ℱ, ℙ, ℚ) with probability measure ℙ and a complex-valued measure ℚ. Assume that ℚ is locally dominated by ℙ with Radon-Nikodym derivative M, i.e.,
        
    M_t=.d ℚ/d ℙ|_ℱ_t,    t ≥ 0

Then for any finite stopping time T, we have ℚ_T ≪ℙ_T and 

    M_T = .d ℚ/d ℙ|_ℱ_T

    

    The lemma is an extension of <cit.> (Theorem III.3.4.(ii)).
    
Since ℱ_T ∧ n⊆ℱ_n, we see that ℚ_T ∧ n≪ℙ_T ∧ n, and as τ∧ n is a bounded stopping time, it follows by the optional stopping theorem that
  
    .d ℚ/d ℙ|_ℱ_T ∧ n = E^ℙ[M_n|ℱ_T ∧ n] = M_T ∧ n.


To prove the theorem, it is enough to show ℚ (A) = ∫_A M_T dℙ for A ∈ℱ_T. Then choose A ∈ℱ_T, we have

    E^ℙ(1_A M_T)     =∑_n ≥ 1 E^ℙ(1_A 1_{n-1 ≤ T<n} E^ℙ[M_n |ℱ_T])=∑_n ≥ 1 E^ℙ(1_A 1_{n-1 ≤ T<n} M_n) 
        =∑_n ≥ 1ℚ(A ∩{n-1 ≤ T<n})=ℚ(A) .

And the result follows.
    

    Now we begin our proof of Theorem 1.
    

    When T is an independent time change, use iterated conditioning:
        
    E[e^i u X_t] =E[E[e^i u L_T_t| T_t]] = E [e^T_tΨ_L(u)]=  ϕ_T(-iΨ_L(u); t).


    When X is not independent of T. First we show that M is a (ℙ, ℱ_T)-martingale. Denote
    
    N_t = exp(iuL_t - tΨ_L(u)).

    which is a complex-valued martingale. 
 As a result of lemma <ref>, M = N_T is a (ℙ, ℱ_T)-martingale. And the Ch.f. of X_t follows from direct computation. 

 Since L is a Lévy process under filtration ℱ, we have that 
 
    ϕ^ℚ_L(z;t)    = E exp(iz L_t)N_t
       = Eexp(i(u+z)L_t) / e^tΨ_L(u)
        = exp(t(Ψ_L(u+z) - Ψ_L(u)).


 Thus, we have Ψ^ℚ_L(z) = Ψ_L(u+z) - Ψ_L(u).




§ COS METHOD FOR VIX OPTIONS

According to COS method in <cit.>, the price at time t_0 of a European style option is


    v(x, t_0) ≈ e^-r Δ t∑_k=0^N-1Re{ϕ(k π/b-a ; x) e^-i k πa/b-a} V_k,

where a,b are integration range and 

    V_k:=2/b-a∫_a^b v(y, T) cos(k πy-a/b-a) d y,

where v(y, T) is the payoff function w.r.t. the underlying asset with value y.
The analytically formula of V_k for call options is known if the log-asset price has analytical Ch.f.

    V_k^call = 2/b-a∫_0^b K(e^y - 1)cos (k πy-a/b-a)dy.


Since S_t = S_0e^rt + L_V_t - Ψ_L(-i)V_t,
        
    VIX_t     = - 2/τ E[lne^-r τS_t+τ/S_t|ℱ_t]
        = 2/τ(Ψ_L(-i)E[ V_t+τ - V_t|ℱ_t] - E[L_V_t+ τ - L_V_t|ℱ_t])
        = 2/τ E[ V_t+τ - V_t|ℱ_t](Ψ_L(-i) - EL_1)
        = 2/τ(Ψ_L(-i) - EL_1) f(v_t).

In affine models, we can show that VIX_t = mv_t + n for some constants m and m. 

Let a = max(0, K^2 - n/m), for VIX options, we have the analytical expression for the Ch.f. of VIX^2, then 

    V_k^VIX = 2/b∫_a^b (√(my+n) - K)cos (k π y/b)dy = 2/b∫_a^b √(my+n)cos (k π y/b)dy - 2K/bψ_k(a,b),

the integral can be done numerically. Then the pricing formula becomes

    C^VIX(v_t_0, t_0) ≈ e^-r Δ t∑_k=0^N-1Re{ϕ_v(k π/b ;v_t_0)} V_k^VIX.



§ CALIBRATION DETAILS

Spot index prices like SPX and VIX are not used in the calibration because they are not directly traded in the market, see also <cit.>. They are recovered from the market according to the put-call parity as discounted futures price. Meanwhile, such implied spot prices contain the term structure of future dividend expectations.

    
  * SPX options with AM settlement are specially treated with 1 day less maturity
    
  * The risk-free rate is quoted from daily U.S. treasury bond rates with Spline interpolation
    
  * Implied volatility is a function of moneyness k_t = K/F_t (and no additional rates) because
    
    e^-rτ E_t[S_T - K]_+      = C_t(K, T) ≡ C^B S(k_t, K, τ, I V_t) 
        =K e^-r τ[e^k_tΦ(d_1)-Φ(d_2)]

yields

    e^k_tΦ(d_1(k_t, I V_t, τ))-Φ(d_2(k_t, I V_t, τ))=E_t[e^k_T-1]_+.


To price and compute IV, it's enough to obtain the moneyness. By put-call parity,

    C - P = e^-rτ(F - K),

if r is known, then 

    k_t = K/K + (C_t-P_t)e^rτ.






















40


urlstyle

[Baldeaux and Badran(2014)]baldeaux2014consistent
Jan Baldeaux and Alexander Badran.
Consistent modelling of vix and equity derivatives using a 3/2 plus
  jumps model.
Applied Mathematical Finance, 210 (4):0
  299–312, 2014.

[Ballotta and Rayée(2022)]ballotta2022smiles
Laura Ballotta and Grégory Rayée.
Smiles & smirks: Volatility and leverage by jumps.
European Journal of Operational Research, 2980
  (3):0 1145–1161, 2022.

[Bardgett et al.(2019)Bardgett, Gourier, and
  Leippold]bardgett2019inferring
Chris Bardgett, Elise Gourier, and Markus Leippold.
Inferring volatility dynamics and risk premia from the s&p 500 and
  vix markets.
Journal of Financial Economics, 1310 (3):0
  593–618, 2019.

[Bayer et al.(2013)Bayer, Gatheral, and Karlsmark]bayer2013fast
Christian Bayer, Jim Gatheral, and Morten Karlsmark.
Fast ninomiya–victoir calibration of the double-mean-reverting
  model.
Quantitative Finance, 130 (11):0 1813–1829,
  2013.

[Bayer et al.(2016)Bayer, Friz, and Gatheral]bayer2016pricing
Christian Bayer, Peter Friz, and Jim Gatheral.
Pricing under rough volatility.
Quantitative Finance, 160 (6):0 887–904,
  2016.

[Broadie and Kaya(2006)]broadie2006exact
Mark Broadie and Özgür Kaya.
Exact simulation of stochastic volatility and other affine jump
  diffusion processes.
Operations research, 540 (2):0 217–231, 2006.

[Carr and Madan(1999)]carr1999option
Peter Carr and Dilip Madan.
Option valuation using the fast fourier transform.
Journal of computational finance, 20 (4):0
  61–73, 1999.

[Carr and Sun(2007)]carr2007new
Peter Carr and Jian Sun.
A new approach for option pricing under stochastic volatility.
Review of Derivatives Research, 10:0 87–150, 2007.

[Carr and Wu(2004)]carr2004time
Peter Carr and Liuren Wu.
Time-changed lévy processes and option pricing.
Journal of Financial economics, 710 (1):0
  113–141, 2004.

[Carr and Wu(2009)]carr2009variance
Peter Carr and Liuren Wu.
Variance risk premiums.
The Review of Financial Studies, 220 (3):0
  1311–1341, 2009.

[Carr et al.(2003)Carr, Geman, Madan, and Yor]carr2003stochastic
Peter Carr, Hélyette Geman, Dilip B Madan, and Marc Yor.
Stochastic volatility for lévy processes.
Mathematical finance, 130 (3):0 345–382,
  2003.

[Cheng(2019)]cheng2019vix
Ing-Haw Cheng.
The vix premium.
The Review of Financial Studies, 320 (1):0
  180–227, 2019.

[Clark(1973)]clark1973subordinated
Peter K Clark.
A subordinated stochastic process model with finite variance for
  speculative prices.
Econometrica: journal of the Econometric Society, pages
  135–155, 1973.

[Cont and Kokholm(2013)]cont2013consistent
Rama Cont and Thomas Kokholm.
A consistent pricing model for index options and volatility
  derivatives.
Mathematical Finance: An International Journal of Mathematics,
  Statistics and Financial Economics, 230 (2):0 248–274,
  2013.

[Drimus(2012)]drimus2012options
Gabriel G Drimus.
Options on realized variance by transform methods: a non-affine
  stochastic volatility model.
Quantitative Finance, 120 (11):0 1679–1694,
  2012.

[Eberlein and Madan(2009)]eberlein2009correlating
Ernst Eberlein and Dilip B Madan.
On correlating lévy processes.
Robert H. Smith School Research Paper No. RHS, pages 06–118,
  2009.

[Fang and Oosterlee(2009)]fang2009novel
Fang Fang and Cornelis W Oosterlee.
A novel pricing method for european options based on fourier-cosine
  series expansions.
SIAM Journal on Scientific Computing, 310
  (2):0 826–848, 2009.

[Filipović(2001)]filipovic2001general
Damir Filipović.
A general characterization of one factor affine term structure
  models.
Finance and Stochastics, 5:0 389–412, 2001.

[Fouque and Saporito(2018)]fouque2018heston
J-P Fouque and Yuri F Saporito.
Heston stochastic vol-of-vol model for joint calibration of vix and
  s&p 500 options.
Quantitative Finance, 180 (6):0 1003–1016,
  2018.

[Gatheral(2008)]gatheral2008consistent
Jim Gatheral.
Consistent modeling of spx and vix options.
In Bachelier congress, volume 37, pages 39–51, 2008.

[Gatheral et al.(2018)Gatheral, Jaisson, and
  Rosenbaum]gatheral2018volatility
Jim Gatheral, Thibault Jaisson, and Mathieu Rosenbaum.
Volatility is rough.
Quantitative finance, 180 (6):0 933–949,
  2018.

[Gatheral et al.(2020)Gatheral, Jusselin, and
  Rosenbaum]gatheral2020quadratic
Jim Gatheral, Paul Jusselin, and Mathieu Rosenbaum.
The quadratic rough heston model and the joint s&p 500/vix smile
  calibration problem.
arXiv preprint arXiv:2001.01789, 2020.

[Geman et al.(2001)Geman, Madan, and Yor]geman2001time
Hélyette Geman, Dilip B Madan, and Marc Yor.
Time changes for lévy processes.
Mathematical Finance, 110 (1):0 79–96, 2001.

[Glasserman and Kim(2011)]glasserman2011gamma
Paul Glasserman and Kyoung-Kuk Kim.
Gamma expansion of the heston stochastic volatility model.
Finance and Stochastics, 15:0 267–296, 2011.

[Goard and Mazur(2013)]goard2013stochastic
Joanna Goard and Mathew Mazur.
Stochastic volatility models and the pricing of vix options.
Mathematical Finance: An International Journal of Mathematics,
  Statistics and Financial Economics, 230 (3):0 439–458,
  2013.

[Guo et al.(2022)Guo, Loeper, Obłój, and Wang]guo2022joint
Ivan Guo, Grégoire Loeper, Jan Obłój, and Shiyi Wang.
Joint modeling and calibration of spx and vix by optimal transport.
SIAM Journal on Financial Mathematics, 130
  (1):0 1–31, 2022.

[Guyon(2020)]guyon2020joint
Julien Guyon.
The joint s&p 500/vix smile calibration puzzle solved.
Risk, April, 2020.

[Jacod and Shiryaev(2013)]jacod2013limit
Jean Jacod and Albert Shiryaev.
Limit theorems for stochastic processes, volume 288.
Springer Science & Business Media, 2013.

[Jacquier et al.(2018)Jacquier, Martini, and Muguruza]jacquier2018vix
Antoine Jacquier, Claude Martini, and Aitor Muguruza.
On vix futures in the rough bergomi model.
Quantitative Finance, 180 (1):0 45–61, 2018.

[Kokholm and Stisen(2015)]kokholm2015joint
Thomas Kokholm and Martin Stisen.
Joint pricing of vix and spx options with stochastic volatility and
  jump models.
The Journal of Risk Finance, 160 (1):0 27–48,
  2015.

[Küchler and Sorensen(2006)]kuchler2006exponential
Uwe Küchler and Michael Sorensen.
Exponential families of stochastic processes.
Springer Science & Business Media, 2006.

[Lian and Zhu(2013)]lian2013pricing
Guang-Hua Lian and Song-Ping Zhu.
Pricing vix options with stochastic volatility and random jumps.
Decisions in Economics and Finance, 36:0 71–88, 2013.

[Lian and Zhu(2012)]lian2012pricing
Guanghua Lian and Song-Ping Zhu.
Pricing vix options with stochastic volatility and random jumps.
Available at SSRN 2065065, 2012.

[Lin and Chang(2010)]lin2010consistent
Yueh-Neng Lin and Chien-Hung Chang.
Consistent modeling of s&p 500 and vix derivatives.
Journal of Economic Dynamics and Control, 340
  (11):0 2302–2319, 2010.

[Luciano and Schoutens(2006)]luciano2006multivariate
Elisa Luciano and Wim Schoutens.
A multivariate jump-driven financial asset model.
Quantitative finance, 60 (5):0 385–402, 2006.

[Luciano and Semeraro(2007)]luciano2007extending
Elisa Luciano and Patrizia Semeraro.
Extending time-changed lévy asset models through multivariate
  subordinators.
2007.

[Mendoza-Arriaga et al.(2010)Mendoza-Arriaga, Carr, and
  Linetsky]mendoza2010time
Rafael Mendoza-Arriaga, Peter Carr, and Vadim Linetsky.
Time-changed markov processes in unified credit-equity modeling.
Mathematical Finance: An International Journal of Mathematics,
  Statistics and Financial Economics, 200 (4):0 527–569,
  2010.

[Pacati et al.(2018)Pacati, Pompa, and Renò]pacati2018smiling
Claudio Pacati, Gabriele Pompa, and Roberto Renò.
Smiling twice: the heston++ model.
Journal of Banking & Finance, 96:0 185–206, 2018.

[Papanicolaou(2022)]papanicolaou2022consistent
Andrew Papanicolaou.
Consistent time-homogeneous modeling of spx and vix derivatives.
Mathematical Finance, 320 (3):0 907–940,
  2022.

[Papanicolaou and Sircar(2014)]papanicolaou2014regime
Andrew Papanicolaou and Ronnie Sircar.
A regime-switching heston model for vix and s&p 500 implied
  volatilities.
Quantitative Finance, 140 (10):0 1811–1827,
  2014.






