# Option Pricing under Two-Factor Stochastic Volatility Jump-Diffusion Model

Authors: Guohe Deng
Journal: Complexity
 (2020)

---

## Abstract

Empirical evidence shows that single-factor stochastic volatility models are not flexible enough to account for the stochastic behavior of the skew, and certain financial assets may exhibit jumps in returns and volatility. This paper introduces a two-factor stochastic volatility jump-diffusion model in which two variance processes with jumps drive the underlying stock price and then considers the valuation on European style option. We derive a semianalytical formula for European vanilla option and develop a fast and accurate numerical algorithm for the computation of the option prices using the fast Fourier transform (FFT) technique. We compare the volatility smile and probability density of the proposed model with those of alternative models, including the normal jump diffusion model and single-factor stochastic volatility model with jumps, respectively. Finally, we provide some sensitivity analysis of the model parameters to the options and several calibration tests using option market data. Numerical examples show that the proposed model has more flexibility to capture the implied volatility term structure and is suitable for empirical work in practice.


---

## Body

## 1. Introduction

Options’ pricing has played an important issue in the general theory of asset pricing since the celebrated work of Black and Scholes [1] (called by BS model), who assume the asset returns satisfying a Geometric Brownian motion with constant parameters. However, option pricing formulas in the BS model do not perform well empirically in practice because of the stylized facts, such as the asymmetric leptokurtic features and the volatility smile. Due to the limitations of the BS model, there has been a wide range of research carried out to improve upon this model. The BS model has been extended to incorporate either constant volatility plus jumps or stochastic volatility with/without jumps to capture these stylized facts. Merton [2] was the first to introduce a normal jump-diffusion model which can be consistent with the extreme movements sometimes observed in financial time series and lead to the leptokurtic feature, volatility smile, and closed-form solutions for options. Heston [3] employed the square-root process for modeling the variance of the underlying asset upon the BS model and showed how to use the Fourier transform to obtain the analytical formula of vanilla options. Bates [4], Scott [5], Zhang and Wang [6, 7], Deng [8], and Pillay and O’Hara [9] combined Heston’s stochastic volatility model with Merton’s normal jump-diffusion model to model the asset price dynamics and obtained the explicit formulas for European options. Eraker et al. [10] and Chernov [11] investigated Heston’s stochastic volatility model incorporating jumps both in returns and volatility. Duffie et al. [12] suggested a generalized affine jump-diffusion model (AJD), in which the drift vector, covariance matrix, and jump intensities all have affine dependence on the state vector. Consequently, the AJD model has become a main topic for research studies on equity returns and on currencies, among other applications (see Glasserman and Kim [13] and and Kaeck and Alexandex [14], for example). In the case of single-factor stochastic volatility model, one of the parameters of the model is the correlation coefficient between the underlying asset and its instantaneous variance. However, the assumption of a single factor in this model does not provide sufficient flexibility to account for the stochastic behavior of the volatility skew observed in option markets and leads to incorrect pricing formulas for deep in-and-out-of-the money European options.In a parallel development, two-factor stochastic volatility models (TSV) have generated attention from the option-pricing literature. Christonffersen et al. [15] and Fouque et al. [16] extend Heston’s model to generate a two-factor stochastic volatility model, which is able to account for the stochastic correlation between the asset return and its instantaneous variance. Another extension to Heston’s model is proposed by [17] who consider a matrix Wishart diffusion process to introduce a correlation structure between the asset returns and its volatility factors (denoted as the MSV model). Gourieroux [18] and Gourieroux and Sufana [19] consider the Wishart Autoregressive (WAR) diffusion processes to model volatility process of the underlying asset. Recently, Muhle-Karbe et al. [20] develop a multivariate Ornstein–Uhlenbeck stochastic volatility model based on Lévy process. As demonstrated in [17], the main advantage of the MSV model is that it can calibrate short-term and long-term volatility levels better than a single-factor one. Recently, the MSV model has been applied to value exotic options (see Wong and Chan [21], Romo [22], and Zhong and Deng [23], for example).Although the MSV model may be useful for option pricing, it is unsuccessful in the modeling of certain financial asset returns. There is empirical evidence suggesting that volatility jumps cannot be recaptured by a diffusion model only (see Jacod and Todorov [24]). Bates [25] and Chernov et al. [26] find that neither one volatility factor model with jumps nor affine two-factor specifications are well equipped to handle the very short-lived but erratic extreme events typically found during financial crisis. For these reasons, this paper aims to extend the MSV model incorporating jump components in asset returns and its volatilities, which can improve the pricing performance of the AJD and MSV models, and also offer a potentially convenient framework for modeling term structure.This paper contributes to the literature in the following ways. We derive a semianalytic formulas for the prices of European vanilla options in a two-factor stochastic volatility jump-diffusion model by using the martingale method and Fourier transform technique. This enables us to understand how our proposed model can be adjusted to fit volatility smiles and the real data. We also obtain a fast and accurate numerical solution for the European option pricing by the Fast Fourier transform (FFT) technique. This allows us to assess the impacts of the jump components in the two-factor volatility on option price.The outline of this paper is as follows. Section2 introduces the model. The characteristic function of the log-stock price and the semianalytical formula for the European call option are derived. An approximation solution for this option pricing by FFT is also provided. Computational results and empirical analysis are reported in Section 3. Section 4 concludes the paper. In this section, we proposed the following model for the asset price dynamics which is driven by two-factor stochastic volatility jump-diffusion models. Under this proposed model, we derive the semianalytic form for option price and then provide a fast and accurate numerical method using FFT technique.


## 2. Asset Price Dynamics

Consider an arbitrage-free, frictionless financial economy consisting of two assets, namely, a bank accountB and a risky stock S, which are traded continuously over a finite time interval 0,T. Let Ω,ℱ,ℱtt≥0,Q be a complete probability space with filtrated σ-field ℱt, on which are defined a four-dimensional standard Wiener process Wt=W1t,W2t,W3t,W4t, three-dimensional Poisson processes Nt=N1,t,N2,t,N3,t, and three-dimensional random jump size sequences J=J1,J2,J3. Assume that all uncertainties on this probability space are mutually independent and ℱt is generated by these uncertainties.We assume that the joint dynamics of the stock priceSt and its variance processes are specified under some risk-neutral probability measure Q (see [12]) by the following system of stochastic differential equation:(1)dStSt=r−δ−λμdt+v1,tρ1dW1,t+1−ρ12dW3,t+v2,tρ2dW2,t+1−ρ22dW4,t+eJ1−1dN1,t,dv1,t=κ1θ1−v1,tdt+σ1v1,tdW1,t+J2dN2,t,dv2,t=κ2θ2−v2,tdt+σ2v2,tdW2,t+J3dN3,t,where r and δ are constant risk-free interest and dividend yields, respectively, v1,t and v2,t are stochastic variances with the mean reversion speeds κ1 and κ2, the long-term variances θ1 and θ2, and the volatilities of the variances σ1 and σ2, respectively, and ρ1 and ρ2 denote the correlation coefficients between the Wiener processes driving the log-stock price and the two variance processes, respectively. The Poisson processes Ni,t are assumed to have constant intensities λi for i=1,2,3. The jump size distribution has an associated transform of Θc1,c2,c3, such that μ=Θ1,0,0−1, where the jump transform Θ is defined by(2)Θc1,c2,c3=λ1Θ1c1+λ2Θ2c2+λ3Θ3c3λ,with λ=λ1+λ2+λ3 and(3)Θ1c=expμsc+12σs2c2,Θ2c=11−η1c,Θ3c=11−η2c.All parameters in this specification are assumed to be constant. To ensure positivity ofvj,t the parameters must be chosen such that 2κjθj≥σj2 for j=1,2. We will refer this specification as the TSVIJ model.The TSVIJ model nests some popular models in option pricing. Without considering jumps in both returns and volatilities, the TSVIJ model can be reduced to the TSV models (called also the double Heston model) considered by Chernov et al. [26] and LeBaron [27]. If we consider a single-factor volatility process (denoted the SVIJ model for short) only in (1), the SVIJ model discussed by Chernov [11] can be regarded as a particular case of the TSVIJ model. When we consider the single-factor volatility process without jumps, the TSVIJ model can be reduced to the SVJ model discussed by Bates [4]. Without jumps in both returns and variance process and with considering a single-factor volatility process only, the TSVIJ model is degenerated to the SV model considered by Heston [3]. We have the normal jump-diffusion model (denoted as Merton model) considered by Merton [2] without considering the stochastic volatility model in (1).The TSVIJ model specified in (1), unlike the SVIJ and SV model, accounts for a richer variance-covariance structure. In particular, the total conditional variance of the stock return process is the sum of the two variance factors plus jumps:(4)Vt≔1dtVardlnStℱt=^1dtVartdlnSt=v1,t+v2,t+λ1μs2+σs2,whereas, as shown by Eraker et al. [10], the conditional variance of the variance processes and the correlation between the stock return and the variance process are given by(5)1dtVartdVt=∑j=12σj2vj,t+2λj+1ηj2,(6)ρSt,Vt=CorrtdlnSt,dVt=ρ1σ1v1,t+ρ2σ2v2,tσ12v1,t+σ22v2,t+2λ2η12+λ3η22×1v1,t+v2,t+λ1μs2+σs2,respectively. Note that the correlation between the stock return and variance is determined by ρ1 and ρ2 and depends on the current levels of both the factors and jumps. It is an important scale for modeling the leverage effect which implies a stochastic correlation.There are two important advantages of the TSVIJ model with respect to the SVIJ/or SV model. First, the TSVIJ model displays not only stochastic variance but also stochastic correlation between stock return and variance, and this feature potentially enables the model to capture fluctuations in option skewness. Second, the TSVIJ model provides more flexibility to modeling volatility term structure; one of the factors determines the correlation between short-term stock return and variance process, while the other factor determines long-term stock return and variance process.


### 2.1. The Characteristic Function

Given the asset price dynamics specified in (1), it is possible to obtain the characteristic function of the log-stock price. Let Xt=lnSt. Applying the generalized Itô formula to the first equation in (1), the evolution of the log-stock price is described by(7)dXt=r−δ−v1,t+v2,t2−λμdt+v1,tρ1dW1,t+1−ρ12dW3,t+v2,tρ2dW2,t+1−ρ22dW4,t+J1dN1,t.In the following, we use the martingale method and generalized Itô formula to derive the characteristic function whose solution is present in the lemma below.Lemma 1.
Suppose that the dynamics of processesXt,v1,t and v2,t are given by (1). Define the conditional characteristic function of Xt as(8)ϕt≔ϕt,T,x,v1,v2;u=EeiuXTℱt=EeiuXTXt=x,v1,t=v1,v2,t=v2,where u is a real-valued argument. Then,(9)ϕt,T,x,v1,v2;u=expiux+At,T,u+B1t,T,uv1+B2t,T,uv2,where At,T,u,B1t,T,u, and B2t,T,u are given by(10)Bjt,T,u=1σj2aju+γju−2γju1−djuexp−γjuT−t,(11)At,T,u=κ1θ1σ12a1u−γ1uT−t−2κ1θ1σ12ln1−d1uexp−γ1uT−t1−d1u+κ2θ2σ22a2u−γ2uT−t−2κ2θ2σ22ln1−d2uexp−γ2uT−t1−d2u+λ1Θ1iu+iur−δ−λμ−λ+λ2g1u+λ3g2uT−t+λ2γ1u1g1u+1e1u×lng1u−e1ud1ue−γ1uT−tg1u−e1ud1u+λ3γ2u1g2u+1e2u×lng2u−e2ud2ue−γ2uT−tg2u−e2ud2u,with aju=κj−iuρjσj,bu=1/2iuiu−1, γju=aj2u−2σj2bu, dju=aju−γju/aju+γju,eju=1−ηj/σj2aju+γju, and gju=eju+2ηjγju/σj2 for j=1,2.Proof.
See Appendix A.
With the characteristic function available, European options can be valued using the Fourier transform method (FT). There are increasing research efforts in option pricing using the FT method. For example, Carr and Madan [28] illustrated the fundamental idea of using the Fast Fourier Transform (FFT) for computing European vanilla options based on the characteristic function of the log-asset value. Duffie et al. [12] provide a comprehensive survey on the FT approach to option pricing in a wide range of stochastic processes. Dempster and Hoog [29] developed an approximation method for pricing European options on spreads using the two-dimensional FFT. Pillay and O’Hara [9] and Zhang and Wang [6, 7] also use the FFT technique to pricing European options under the stochastic volatility model with/without jumps. Wong and Zhao [30] use a modified FFT method (FRFT), namely, the fractional FFT, to value currency options. Černý [31] presented a detailed discussion on the implementation of FFT to option pricing. In Section 2.2, we examine how these methods can be used to our case.


### 2.2. Semianalytical Formula

Now, we consider a European vanilla call option pricing. LetK be the strike price and T the maturity of a European call option. The price at time 0 of the European call option, denoted by CS, is computed as the discounted risk-neutral conditional expectation of the terminal payoff maxST−K,0=ST−K+ as follows:(12)CS=Ee−rTST−K+.Letk=lnK. Then, the call option function (12) can be expressed as a function of k in integral form as follows:(13)Ck=e−rT∫k+∞exT−ekqxTdxT,where qxT is the probability density function of the process XT. We rewrite (13) as two expectations of indicator functions so that(14)Ck=e−δTEe−r−δTeXT1XT≥k−Ke−rTE1XT≥k.For the first term in (14), we choose the spot price St as numeraire and switch from probability measure Q to Q1 on Ω,ℱ. The Radon–Nikodym derivative is then given by(15)dQ1dQℱT=e−r−δTSTS0=e−r−δTeXTEQe−r−δTeXT.Under this new equivalent martingale measureQ1, the call option price (14) can be restated as(16)Ck=S0e−δTQ1XT≥k−Ke−rTQXT≥k.The two probabilities in (16) can be calculated using the Fourier inversion transform:(17)Q1XT≥k=12+1π∫0+∞ℜe−iukf1uiudu,(18)QXT≥k=12+1π∫0+∞ℜe−iukf0uiudu,where ℜ⋅ denotes the real part, and f1u and f0u are defined by(19)f1u=EQ1eiuXT=Eeiu+1XTEeXT=ϕ0,T,x,v1,v2;u−iϕ0,T,x,v1,v2;−i,f0u=EeiuXT=ϕ0,T,x,v1,v2;u,which are the corresponding characteristic functions of XT under the measures Q1 and Q, respectively. Substituting (17) and (18) into (16), we summarize the result as follows.Theorem 1.
Assuming the dynamics of the asset price specified by (1), the semianalytic formula for the European call option can be given by(20)Ck=S0e−δT12+1π∫0+∞ℜe−iukf1uiudu−Ke−rT12+1π∫0+∞ℜe−iukf0uiudu.From the Fourier transform formula, the probability density function of the log-asset price is given by(21)qxT=1π∫0+∞ℜe−iuxTf0udu.


### 2.3. FFT Method

Note that the option priceCk in (20) does not decay to 0 as k⟶−∞. Therefore, the FFT method cannot be used directly to compute the integrals in (20). As in Carr and Madan [28] and Dempster and Hoog [29], Ck is multiplied by an exponentially damping factor eαk for α>0 so that it is square-integrable in k over −∞,+∞. Define a modified call option price function ck=eαkCk. Then, the Fourier transform of ck is(22)ψu=∫−∞∞eiukckdk=e−rT∫−∞∞eiukeαk∫k∞eXT−ekqxTdxTdk=e−rT∫−∞∞exTqxTdxT∫−∞xTeiu+αkdk−e−rT∫−∞∞qxTdxT∫−∞xTeiu+α+1kdk=e−rT∫−∞∞eiu+α+1xTqxTdxTiu+αiu+α+1=e−rTϕ0,T,x,v1,v2;u−iα−iiu+αiu+α+1.Given the Fourier transform functionψu, the integral Ck can be calculated by the inverse Fourier transform:(23)Ck=e−αk2π∫−∞∞e−iukψudu=e−αkπ∫0∞e−iukψudu.Setuj=jΔ,j=0,1,…,N−1, where Δ denotes the integration step width and N is the total number of grid points and is assumed to be an arbitrary power of 2. Then, invoking the trapezoid rule, the Fourier integral in (23) can be approximated by the summation:(24)Ck≈e−αkπ∑j=0N−1e−iujkψujΔ.Recall that a one-dimensional Fast Fourier transform (FFT) computes for any complex (input) array,χj∈C:j=0,1,…,N−1, and the following (output) array of identical structure is obtained:(25)Yl=∑j=0N−1e−2πi/Njlχj,for all l=0,1,…,N−1. In order to use the FFT algorithm to evaluate the sum in (14), we define a grid of size N×1 along the modified log-strike price k by Λ=kv:0≤v≤N−1, where kv=v−N/2ω for ω>0. Choosing ωΔ=2π/N gives the following values of the sum Ykv on the grid Λ as(26)Ykv=∑j=0N−1e−iujkvψuj.Substituting the representations ofuj and kv in (25) and using the choice of ωΔ, we obtain(27)Ykv=∑j=0N−1e−ijΔv−N/2ωψuj=∑j=0N−1e−2πi/Njv−1jψuj.Using Simpson’s rule for numerical integration, define a sequence of weighting factors by(28)wj=13,j=0,43,jis odd,23,jis even.Then, (26) can be rewritten as(29)Ykv=∑j=0N−1e−2πi/Njv−1jψjΔwj.This can be computed from the expression of the FFT (25) by taking the input array as(30)χj=−1jψjΔwj,∀j=0,1,…,N−1,such that(31)Ykv=∑j=0N−1e−2πi/Njvχj.The result is an approximation for the European call option price atN×1 different (log) strikes given by(32)Ckv≈e−αkvπYkvΔ.All that remains is to choose the value ofα. Carr and Madan [28] suggest that α is chosen such that ESTα+1<∞, that is, ψ0<∞. Therefore, ck is well behaved.Note that the FFT algorithm above depends on the intrinsic value of the option. However, the call options with very short maturity or out-of-the money options have no intrinsic value. This causes the integrand in the Fourier inversion to be high oscillate. As shown in Carr and Madan [28] and Zhang and Wang [6, 7], an FFT based on the time value of the option has to be obtained.LetC˜k be the time value of an out-of-the-money option; then,(33)C˜k=e−rT∫−∞∞ek−exT1xT<k,k<0+exT−ek1xT>k,k>0qxTdxT.Similar to (22), the Fourier transform of the modified time value c˜k=sinhαkC˜k is given by(34)ψ˜u=∫−∞∞eiukc˜kdk=e−rT2ωu−iα−ωu+iα,where ωu=1/1+iu−erT/iu−ϕ0,T,x,v1,v2;u−i/u2−iu. So, the time value, C˜k, can be calculated by the inverse Fourier transform:(35)C˜k=1πsinhαk∫0∞e−iukψ˜udu.Similar to (31), the FFT algorithm for the out-of-the-money options is(36)C˜kv≈Δπsinhαk∑j=0N−1e−2πi/Njvχ˜j,where χ˜j=−1jψ˜jΔwj.Remark 1.
In practice, one may determine an upper bound on the damping parameterα from the analytical expression for the characteristic function and the condition ESTα+1<∞. Carr and Madan [28] found that one-fourth of this upper bound serves as a good choice for α (>1 for calls, >0 for puts) and the impact of the damping factor on the option price is not significant.


## 2.1. The Characteristic Function

Given the asset price dynamics specified in (1), it is possible to obtain the characteristic function of the log-stock price. Let Xt=lnSt. Applying the generalized Itô formula to the first equation in (1), the evolution of the log-stock price is described by(7)dXt=r−δ−v1,t+v2,t2−λμdt+v1,tρ1dW1,t+1−ρ12dW3,t+v2,tρ2dW2,t+1−ρ22dW4,t+J1dN1,t.In the following, we use the martingale method and generalized Itô formula to derive the characteristic function whose solution is present in the lemma below.Lemma 1.
Suppose that the dynamics of processesXt,v1,t and v2,t are given by (1). Define the conditional characteristic function of Xt as(8)ϕt≔ϕt,T,x,v1,v2;u=EeiuXTℱt=EeiuXTXt=x,v1,t=v1,v2,t=v2,where u is a real-valued argument. Then,(9)ϕt,T,x,v1,v2;u=expiux+At,T,u+B1t,T,uv1+B2t,T,uv2,where At,T,u,B1t,T,u, and B2t,T,u are given by(10)Bjt,T,u=1σj2aju+γju−2γju1−djuexp−γjuT−t,(11)At,T,u=κ1θ1σ12a1u−γ1uT−t−2κ1θ1σ12ln1−d1uexp−γ1uT−t1−d1u+κ2θ2σ22a2u−γ2uT−t−2κ2θ2σ22ln1−d2uexp−γ2uT−t1−d2u+λ1Θ1iu+iur−δ−λμ−λ+λ2g1u+λ3g2uT−t+λ2γ1u1g1u+1e1u×lng1u−e1ud1ue−γ1uT−tg1u−e1ud1u+λ3γ2u1g2u+1e2u×lng2u−e2ud2ue−γ2uT−tg2u−e2ud2u,with aju=κj−iuρjσj,bu=1/2iuiu−1, γju=aj2u−2σj2bu, dju=aju−γju/aju+γju,eju=1−ηj/σj2aju+γju, and gju=eju+2ηjγju/σj2 for j=1,2.Proof.
See Appendix A.
With the characteristic function available, European options can be valued using the Fourier transform method (FT). There are increasing research efforts in option pricing using the FT method. For example, Carr and Madan [28] illustrated the fundamental idea of using the Fast Fourier Transform (FFT) for computing European vanilla options based on the characteristic function of the log-asset value. Duffie et al. [12] provide a comprehensive survey on the FT approach to option pricing in a wide range of stochastic processes. Dempster and Hoog [29] developed an approximation method for pricing European options on spreads using the two-dimensional FFT. Pillay and O’Hara [9] and Zhang and Wang [6, 7] also use the FFT technique to pricing European options under the stochastic volatility model with/without jumps. Wong and Zhao [30] use a modified FFT method (FRFT), namely, the fractional FFT, to value currency options. Černý [31] presented a detailed discussion on the implementation of FFT to option pricing. In Section 2.2, we examine how these methods can be used to our case.


## 2.2. Semianalytical Formula

Now, we consider a European vanilla call option pricing. LetK be the strike price and T the maturity of a European call option. The price at time 0 of the European call option, denoted by CS, is computed as the discounted risk-neutral conditional expectation of the terminal payoff maxST−K,0=ST−K+ as follows:(12)CS=Ee−rTST−K+.Letk=lnK. Then, the call option function (12) can be expressed as a function of k in integral form as follows:(13)Ck=e−rT∫k+∞exT−ekqxTdxT,where qxT is the probability density function of the process XT. We rewrite (13) as two expectations of indicator functions so that(14)Ck=e−δTEe−r−δTeXT1XT≥k−Ke−rTE1XT≥k.For the first term in (14), we choose the spot price St as numeraire and switch from probability measure Q to Q1 on Ω,ℱ. The Radon–Nikodym derivative is then given by(15)dQ1dQℱT=e−r−δTSTS0=e−r−δTeXTEQe−r−δTeXT.Under this new equivalent martingale measureQ1, the call option price (14) can be restated as(16)Ck=S0e−δTQ1XT≥k−Ke−rTQXT≥k.The two probabilities in (16) can be calculated using the Fourier inversion transform:(17)Q1XT≥k=12+1π∫0+∞ℜe−iukf1uiudu,(18)QXT≥k=12+1π∫0+∞ℜe−iukf0uiudu,where ℜ⋅ denotes the real part, and f1u and f0u are defined by(19)f1u=EQ1eiuXT=Eeiu+1XTEeXT=ϕ0,T,x,v1,v2;u−iϕ0,T,x,v1,v2;−i,f0u=EeiuXT=ϕ0,T,x,v1,v2;u,which are the corresponding characteristic functions of XT under the measures Q1 and Q, respectively. Substituting (17) and (18) into (16), we summarize the result as follows.Theorem 1.
Assuming the dynamics of the asset price specified by (1), the semianalytic formula for the European call option can be given by(20)Ck=S0e−δT12+1π∫0+∞ℜe−iukf1uiudu−Ke−rT12+1π∫0+∞ℜe−iukf0uiudu.From the Fourier transform formula, the probability density function of the log-asset price is given by(21)qxT=1π∫0+∞ℜe−iuxTf0udu.


## 2.3. FFT Method

Note that the option priceCk in (20) does not decay to 0 as k⟶−∞. Therefore, the FFT method cannot be used directly to compute the integrals in (20). As in Carr and Madan [28] and Dempster and Hoog [29], Ck is multiplied by an exponentially damping factor eαk for α>0 so that it is square-integrable in k over −∞,+∞. Define a modified call option price function ck=eαkCk. Then, the Fourier transform of ck is(22)ψu=∫−∞∞eiukckdk=e−rT∫−∞∞eiukeαk∫k∞eXT−ekqxTdxTdk=e−rT∫−∞∞exTqxTdxT∫−∞xTeiu+αkdk−e−rT∫−∞∞qxTdxT∫−∞xTeiu+α+1kdk=e−rT∫−∞∞eiu+α+1xTqxTdxTiu+αiu+α+1=e−rTϕ0,T,x,v1,v2;u−iα−iiu+αiu+α+1.Given the Fourier transform functionψu, the integral Ck can be calculated by the inverse Fourier transform:(23)Ck=e−αk2π∫−∞∞e−iukψudu=e−αkπ∫0∞e−iukψudu.Setuj=jΔ,j=0,1,…,N−1, where Δ denotes the integration step width and N is the total number of grid points and is assumed to be an arbitrary power of 2. Then, invoking the trapezoid rule, the Fourier integral in (23) can be approximated by the summation:(24)Ck≈e−αkπ∑j=0N−1e−iujkψujΔ.Recall that a one-dimensional Fast Fourier transform (FFT) computes for any complex (input) array,χj∈C:j=0,1,…,N−1, and the following (output) array of identical structure is obtained:(25)Yl=∑j=0N−1e−2πi/Njlχj,for all l=0,1,…,N−1. In order to use the FFT algorithm to evaluate the sum in (14), we define a grid of size N×1 along the modified log-strike price k by Λ=kv:0≤v≤N−1, where kv=v−N/2ω for ω>0. Choosing ωΔ=2π/N gives the following values of the sum Ykv on the grid Λ as(26)Ykv=∑j=0N−1e−iujkvψuj.Substituting the representations ofuj and kv in (25) and using the choice of ωΔ, we obtain(27)Ykv=∑j=0N−1e−ijΔv−N/2ωψuj=∑j=0N−1e−2πi/Njv−1jψuj.Using Simpson’s rule for numerical integration, define a sequence of weighting factors by(28)wj=13,j=0,43,jis odd,23,jis even.Then, (26) can be rewritten as(29)Ykv=∑j=0N−1e−2πi/Njv−1jψjΔwj.This can be computed from the expression of the FFT (25) by taking the input array as(30)χj=−1jψjΔwj,∀j=0,1,…,N−1,such that(31)Ykv=∑j=0N−1e−2πi/Njvχj.The result is an approximation for the European call option price atN×1 different (log) strikes given by(32)Ckv≈e−αkvπYkvΔ.All that remains is to choose the value ofα. Carr and Madan [28] suggest that α is chosen such that ESTα+1<∞, that is, ψ0<∞. Therefore, ck is well behaved.Note that the FFT algorithm above depends on the intrinsic value of the option. However, the call options with very short maturity or out-of-the money options have no intrinsic value. This causes the integrand in the Fourier inversion to be high oscillate. As shown in Carr and Madan [28] and Zhang and Wang [6, 7], an FFT based on the time value of the option has to be obtained.LetC˜k be the time value of an out-of-the-money option; then,(33)C˜k=e−rT∫−∞∞ek−exT1xT<k,k<0+exT−ek1xT>k,k>0qxTdxT.Similar to (22), the Fourier transform of the modified time value c˜k=sinhαkC˜k is given by(34)ψ˜u=∫−∞∞eiukc˜kdk=e−rT2ωu−iα−ωu+iα,where ωu=1/1+iu−erT/iu−ϕ0,T,x,v1,v2;u−i/u2−iu. So, the time value, C˜k, can be calculated by the inverse Fourier transform:(35)C˜k=1πsinhαk∫0∞e−iukψ˜udu.Similar to (31), the FFT algorithm for the out-of-the-money options is(36)C˜kv≈Δπsinhαk∑j=0N−1e−2πi/Njvχ˜j,where χ˜j=−1jψ˜jΔwj.Remark 1.
In practice, one may determine an upper bound on the damping parameterα from the analytical expression for the characteristic function and the condition ESTα+1<∞. Carr and Madan [28] found that one-fourth of this upper bound serves as a good choice for α (>1 for calls, >0 for puts) and the impact of the damping factor on the option price is not significant.


## 3. Numerical Experiments

In this section, we present some numerical results. First, we report the numerical experiments that were performed on pricing options using two techniques: the FFT and the Monte Carlo (MC) simulation method. Next, we use the FFT to further investigate the volatility smile and the distributional properties implied by our proposed model (TSVIJ). We analyze further the effects of volatility jump components and correlation coefficients on option prices under the TSVIJ model. Finally, we describe the market data used for empirical study and estimate the risk-neutral model parameters from option prices on a set of days in various models including the TSVIJ, the SVIJ, SV, Merton, and the BS models as benchmarks.


### 3.1. Performance of FFT

Having presented the FFT algorithm as outlined above, it would be interesting to investigate the performance of FFT implemented to price option under our proposed model. We have also implemented the MC method for comparison purposes. For the FFT method, according to Remark1 and by means of the given setting of model parameters in Table 1, we used the damping parameter α=1.9, Δ=0.5, and N=2048. For the MC simulation, we apply the Euler scheme for the numerical simulation of the log-asset price and variance processes. In order to achieve high accuracy, the simulation is repeated 10,000 times and the time step is kept at 0.02. Other parameter values used in the numerical computation are listed in Table 1. Two methods were implemented in MATLAB R2015b and on the PC with Intel(R) Pentium(R) D CPU 2.80 GHz. Results of the numerical experiment are presented in Table 2.Table 1
Parameter settings for numerical examples.
Model parametersκ1θ1σ1η1ρ1λ2v1,0κ2θ2σ2η2ρ2λ3v2,0μsσsλ130.060.10.2-0.40.30.0240.080.20.40.60.90.04−0.10.152Market dataS0rδ1000.030.02Table 2
Option prices using FFT and MC methods.
Strike priceKT=0.5T=1FFTMC% differenceFFTMC% difference8024.453524.44360.04129.359929.3706−0.0378521.052821.04930.01726.591526.5940−0.0089018.022817.98400.21624.043524.0488−0.0219515.276015.25620.13121.725821.72590.00010012.803812.8628−0.45919.604319.6138−0.04610510.847910.79020.53817.712217.69960.0731109.02309.01640.07815.971015.96950.0131157.43967.5140−0.98514.421214.40940.0821206.27566.25260.36813.015413.00530.077CPU time (s)0.086906.250.080889.76Table2 reports the results for 9 call option prices with different strike prices (from deep-in-the-money (ITM) to at-the-money (ATM) and deep-out-of-money (OTM)) obtained using the FFT method and the Monte Carlo simulations. Column one lists the strike price K for the options. Columns two to five list the FFT and the MC simulation prices for both T=0.5 and T=1. The numbers in parentheses are difference between the two methods. If we use the Monte Carlo simulation as the benchmark, we can see from Table 2 that the absolute percentage differences of the FFT method are less than 1% for all cases.The FFT method takes less than 2 seconds to produce 150 option prices corresponding to different strike prices. The Monte Carlo simulation takes more than 90 seconds for each option price. These numerical experiments show that our analytic solution is correct, and the FFT method is accurate and efficient.


### 3.2. Implied Volatility

The clear advantage of the FFT method enables us to further investigate the volatility smile implied by our proposed model. For the implied volatility calculations, we use two terms to maturity date: 0.5 year (short-term) and one year (long-term). We use the model parameters taken from Table1 and vary the strike price from 50 to 150, to generate an implied volatility graph. Figure 1 displays the implied volatility curves for three different models, including the Merton, SVIJ, and TSVIJ models.Figure 1
Implied volatility for three models with short and long maturities.
(a)(b)As can be seen from Figure1, for both maturities, the curvature of the implied volatility curves for the TSVIJ model are significantly decreased, and the difference between the Merton, SVIJ, and TSVIJ implied volatility curves is quite large. This should not be a surprise, as jumps in two variance processes increase the conditional variance of the underlying asset. Furthermore, the TSVIJ model generates very greater the left-tail skewness, especially for short maturity, which should be valid for the real-world option market and the best framework describing asset price dynamics. This indicates the advantage of the two-factor stochastic volatility model with jumps, with respect to the jumps only in returns (the Merton model) and the single-factor stochastic volatility model with/without jumps (the SVIJ/SV model), is that the two-factor stochastic volatility model with jumps in volatility provides more flexibility to capture the volatility smile since one of the factors determines the correlation between short-term asset’s returns and variance, whereas the other factor determines the correlation between long-term asset’s returns and variance.Next, we examine how the volatility smile varies with the model parameter values of the variance process. We focus on investigating the effects of three parameters, two volatility jump components,λ2 and η1, and correlation coefficient, ρ1, on the volatility smile. We compute the prices for call option with T=1 year maturity by FFT under the TSVIJ model; then, we take those prices as market prices of option and compute the implied volatilities using the BS model. The results are plotted in Figures 2–4.Figure 2
Implied volatility for different values ofλ2.Figure 3
Implied volatility for different values ofη1.Figure 4
Implied volatility for different values ofρ1.First, we can seen from Figure2 that an increase in λ2, the arrive intensity of the volatility jumps, increases the implied volatility. This implies that the BS model systematically undervalues all options. A similar argument holds for η1, the mean of the volatility jumps (see Figure 3). Second, the levels of the implied volatilities depend on the sign and magnitudes of correlation coefficient ρ1 between the volatility and asset (see Figure 4), whereas the levels of the implied volatilities at the lower strikes drop with increasing correlation ρ1 and ones at higher strikes go up. This feature is consistent with the models with stochastic volatility without jumps in volatility. Thus, the proposed model maintains some properties of the models being of stochastic volatility. We stress, however, that this model enables option prices to fit the forward price term structure exactly.


### 3.3. Probability Densities

We now go one step further and examine the distributional properties of the TSVIJ model using the FFT method. It is also instructive to compare the short-term and long-term probability density functions of TSVIJ, SVIJ, and Merton’s models. From (21), the FFT algorithm of qxT can be approximated by(37)qxT=1π∑j=0N−1e−2πi/Njvf0uj,for v=0,1,…,N−1. The results are plotted in Figures 5–8, which examine the values of the probability density function for a range of values of two parameters when all the others are kept constant and taken from Table 1.Figure 5
Probability density function with different correlationsρ1 and ρ2.Figure 6
Probability density function with different jump intensitiesλ2 and λ3.Figure 7
Probability density function with different jump meansη1 and η2.Figure 8
Comparison of probability density function in three different models.
(a)(b)Figure5 shows how the positive correlations of variance process with the stock return (which is calculated in equation (3)) create a fat-right tail and thin-left tail distribution of the log-return. Positive correlations lead to a high variance when the stock return increases, which fattens the right tail of the density function. This is consistent with those in Wong and Zhao [30].Figures6 and 7 provide the impacts of the volatility jump components, λ2,λ3,η1, and η2 in the TSVIJ model on probability densities of the log-return. Figure 6 shows that increasing the arrive intensities of volatility jumps, λ2 and λ3, will decrease the kurtosis only, but create a fat-right-tail and thin-left-tail. A similar argument holds for the means of the volatility jumps, η1 and η2 (see Figure 7).Figure8 provides comparisons of the short-term (T=0.5 year) and long-term (T=1 year) probability densities of the TSVIJ, SVIJ, and Merton’s models. Figure 8 shows that the probability density of the TSVIJ model is distinctly different; it has very high peak and skewness features compared to those of the SVIJ and Merton models for both maturities. Moreover, for both maturity density curves in Figure 8 still show that the density of the TSVIJ model attenuates heavy tails on both sides. Therefore, the TSVIJ model can correct the pricing bias in the BS model and has a potential application in practice.


### 3.4. Sensitivity Analysis

As said previously, the TSVIJ model provides more flexibility than the single-factor specification to model the volatility term structure and probability density of asset returns. This section considers the effects of the model parameters, including the correlation coefficientsρ1 and ρ2, the arrive intensities λ2 and λ3 of volatility jumps, and the means η1 and η2 of volatility jumps, on the option prices. We analyze the prices of T = 1 year call option only. Model parameter values are taken from Table 1 apart from ρ1,ρ2,λ2,λ3,η1, and η2. The numerical results are displayed in Table 3.Table 3
The impacts ofλ2,λ3, η1,η2, and ρ1,ρ2 on option prices.
Strike priceKρ1=−0.4,ρ2=−0.6ρ1=0.4,ρ2=0.6λ2=0.3λ2=1.3λ2=0.3λ2=1.3λ3=0.9λ3=1.9λ3=0.9λ3=1.9η1=0.02,η2=0.048027.409727.828427.275927.69828524.428724.893524.328924.79989021.707222.207821.649422.15779519.215019.748919.208019.741010016.922117.486016.974417.547910514.958715.530615.059715.641511013.109113.693613.269313.863811511.423212.014811.643712.245612010.057310.638710.317410.9092η1=0.2,η2=0.48029.442632.529329.341432.46208526.642729.978026.580029.95459024.059227.582924.040027.60439521.731525.479921.700025.408810019.531323.417019.620323.542210517.601521.535417.736321.702211015.814019.843616.004920.061611514.175618.285414.424018.555912012.779616.855913.065417.1603From Table3, we can see that (i) option prices are very sensible to η1 and η2. Since η1 and η2 indicate the means of the volatility jumps, this sensitivity can be considered as the sensitivity of option prices to their volatilities, which are verified by equation (2). A similar argument holds for the arrive intensities λ2 and λ3 of the volatility jumps. (ii) Option prices with different strike prices have different sensitivities to the correlations ρ1 and ρ2, whereas in-the-money options (K=80,85,90,95) decrease in values with increasing correlations ρ1 and ρ2, at-the-money options (K=100), and out-the-money options (K=105,110,115,120) go up. This finding is consistent with those analyzed by Zhang and Wang [6] and among others.


### 3.5. Calibration of the TSVIJ Model to Market Data

With a successful implementation of the option pricing formula (15) under the TSVIJ model, we can apply the above calculation method to calibrate the TSVIJ model and estimate the model parameters on a set of market data of standard European-style call options by minimizing the difference between market prices and model prices in a least squared error fit.The dataset used in this study consists of AAPL plain vanilla call option prices traded on the NASDAQ. These option prices are based on the Apple stock without paying dividends, having a value of $ 364.11 on July 5, 2020. We employ the market quotes over the period from July 10, 2020, to July 24, 2020, as the in-sample data to calibrate the model parameters, and use those of July 31, 2020, for the out-of-sample test. The riskless dollar interest rater=0.0175, in this period. To ensure sufficient liquidity in option trading, we disregard option trading volume lower than 100 in the sample data. Finally, the option sample contains 128 call options, with 4 trading days and maturies from 5 days to 26 days. The description of AAPL options data are reported in Table 4.Table 4
Numbers of AAPL option data.
MaturityMoneynessS/KITMATMOTM>1.010.99∼1.01<0.995 days1161312 days1422519 days433526 days2310The TSVIJ model has 17 parameters that need estimation, viz.,Θ=v1,0,κ1,θ1,σ1,η1,ρ1,λ2,v2,0,κ2,θ2,σ2,η2,ρ2,λ3,μs,σs,λ1. The estimation procedure is to minimize the least squared error over the parameter space Θ, i.e., we evaluate(38)minΘSqΘ=1M∑m=1MCS0,r,δ,Km,Tm;Θ−CMarketKm,Tm2+PenaltyΘ,Θ0,where M is the number of the sample data and CMarketKm,Tm and CS0,r,δ,Km,Tm;Θ denote the market price and the model price at time t=0 of a call option, respectively, with maturity Tm and strike price Km. The model price CS0,r,δ,Km,Tm;Θ depending on the parameter space Θ is taken as input and calculated by the FFT method above. As the penalty function, we choose the distance to the initial parameter vector Θ0 (see Mikhailov and Nögel [32]):(39)PenaltyΘ,Θ0=Θ−Θ02.We employ the Lsqnonlin function (least-squares nonlinear) in MATLAB to implement (38). To ensure the success of the local optimization scheme (38), we need effective choices of the initial parameters Θ0. Based on reasonable results from the literature, we set the initial parameters Θ0 used the empirical results studied by Duffie et al. [12]. In addition, we also define lower and upper bounds for the optimal parameters which avoid possible solutions that, while mathematically feasible, are not acceptable in finance sense in the calibration. The resulting calibration results are reported in Table 5.Table 5
Estimated risk-neutral parameters for five given models.
Model parametersTSVIJSVIJSVMertonBSv1,00.03510.04920.07340.06260.1011κ124.372621.91522.9990θ10.06360.06900.5005σ10.31040.15620.0613η10.26000.3073ρ1−0.5315−0.4821−0.5012λ20.51640.9979v2,00.1123κ282.3024θ20.0310σ20.3308η20.6327ρ20.5036λ31.0304μs−0.2839−0.1089−0.1138σs0.37950.20920.1978λ10.89631.96951.9943RMSE0.13060.14210.13090.31160.3625ARPE3.25303.54073.25993.78063.8520CPU time (s)732.5932268.204346.892131.52715.5643Note. RMSE=∑m=1MModel price−Market price2/M and ARPE=1/M∑m=1MModel  price−Market price/Market price.Table5 gives the estimated risk-neutral parameters and pricing performances for the TSVIJ, the SVIJ, SV, Merton, and the BS models using the in-sample data. For comparative purposes, we provide the root mean square error (RMSE), the average relative percentage error (ARPE), and the CPU time. Considering the RMSE and ARPE values, the TSVIJ model best fits the market prices in the in-sample data, followed by the SV, SVIJ, Merton, and the BS models. Taking a closer look at the estimates, the jump-related parameter λ1,λ2,λ3,μs,σs,η1,η2 in the TSVIJ model is significant, suggesting that modeling the asset dynamics and its variance processes cannot ignore the jump components.Using these results in Table5, the model predicted values and its comparison with market prices of the out-of sample data are shown in Table 6. As the table shows, the TSVIJ model provides a good fit for most traded options in a lower RMSE value.Table 6
Out-of-sample pricing results.
Option id.KOption priceTSVIJSVIJSVMertonBS132048.8048.625447.637645.231747.490345.2890236015.6515.658315.641915.664815.455714.6615336512.6512.664112.845113.615812.666612.11574367.511.4011.484611.583911.901011.408710.9635537010.2010.292510.415010.320710.24239.89046372.59.459.48739.338110.21559.16668.894673758.508.60778.35118.76548.18007.97428377.57.207.23167.45137.47697.28007.290593806.406.37676.63677.38386.46346.349410382.56.656.32995.90276.24425.72655.6390113854.955.09675.24516.16745.06494.9924123903.653.69714.13914.29293.94853.8765133953.003.24253.27673.48903.07422.9724144002.342.39552.61482.80882.40082.2510154051.801.84192.11272.23931.88931.6838RMSE0.01670.18391.22290.19471.0379


## 3.1. Performance of FFT

Having presented the FFT algorithm as outlined above, it would be interesting to investigate the performance of FFT implemented to price option under our proposed model. We have also implemented the MC method for comparison purposes. For the FFT method, according to Remark1 and by means of the given setting of model parameters in Table 1, we used the damping parameter α=1.9, Δ=0.5, and N=2048. For the MC simulation, we apply the Euler scheme for the numerical simulation of the log-asset price and variance processes. In order to achieve high accuracy, the simulation is repeated 10,000 times and the time step is kept at 0.02. Other parameter values used in the numerical computation are listed in Table 1. Two methods were implemented in MATLAB R2015b and on the PC with Intel(R) Pentium(R) D CPU 2.80 GHz. Results of the numerical experiment are presented in Table 2.Table 1
Parameter settings for numerical examples.
Model parametersκ1θ1σ1η1ρ1λ2v1,0κ2θ2σ2η2ρ2λ3v2,0μsσsλ130.060.10.2-0.40.30.0240.080.20.40.60.90.04−0.10.152Market dataS0rδ1000.030.02Table 2
Option prices using FFT and MC methods.
Strike priceKT=0.5T=1FFTMC% differenceFFTMC% difference8024.453524.44360.04129.359929.3706−0.0378521.052821.04930.01726.591526.5940−0.0089018.022817.98400.21624.043524.0488−0.0219515.276015.25620.13121.725821.72590.00010012.803812.8628−0.45919.604319.6138−0.04610510.847910.79020.53817.712217.69960.0731109.02309.01640.07815.971015.96950.0131157.43967.5140−0.98514.421214.40940.0821206.27566.25260.36813.015413.00530.077CPU time (s)0.086906.250.080889.76Table2 reports the results for 9 call option prices with different strike prices (from deep-in-the-money (ITM) to at-the-money (ATM) and deep-out-of-money (OTM)) obtained using the FFT method and the Monte Carlo simulations. Column one lists the strike price K for the options. Columns two to five list the FFT and the MC simulation prices for both T=0.5 and T=1. The numbers in parentheses are difference between the two methods. If we use the Monte Carlo simulation as the benchmark, we can see from Table 2 that the absolute percentage differences of the FFT method are less than 1% for all cases.The FFT method takes less than 2 seconds to produce 150 option prices corresponding to different strike prices. The Monte Carlo simulation takes more than 90 seconds for each option price. These numerical experiments show that our analytic solution is correct, and the FFT method is accurate and efficient.


## 3.2. Implied Volatility

The clear advantage of the FFT method enables us to further investigate the volatility smile implied by our proposed model. For the implied volatility calculations, we use two terms to maturity date: 0.5 year (short-term) and one year (long-term). We use the model parameters taken from Table1 and vary the strike price from 50 to 150, to generate an implied volatility graph. Figure 1 displays the implied volatility curves for three different models, including the Merton, SVIJ, and TSVIJ models.Figure 1
Implied volatility for three models with short and long maturities.
(a)(b)As can be seen from Figure1, for both maturities, the curvature of the implied volatility curves for the TSVIJ model are significantly decreased, and the difference between the Merton, SVIJ, and TSVIJ implied volatility curves is quite large. This should not be a surprise, as jumps in two variance processes increase the conditional variance of the underlying asset. Furthermore, the TSVIJ model generates very greater the left-tail skewness, especially for short maturity, which should be valid for the real-world option market and the best framework describing asset price dynamics. This indicates the advantage of the two-factor stochastic volatility model with jumps, with respect to the jumps only in returns (the Merton model) and the single-factor stochastic volatility model with/without jumps (the SVIJ/SV model), is that the two-factor stochastic volatility model with jumps in volatility provides more flexibility to capture the volatility smile since one of the factors determines the correlation between short-term asset’s returns and variance, whereas the other factor determines the correlation between long-term asset’s returns and variance.Next, we examine how the volatility smile varies with the model parameter values of the variance process. We focus on investigating the effects of three parameters, two volatility jump components,λ2 and η1, and correlation coefficient, ρ1, on the volatility smile. We compute the prices for call option with T=1 year maturity by FFT under the TSVIJ model; then, we take those prices as market prices of option and compute the implied volatilities using the BS model. The results are plotted in Figures 2–4.Figure 2
Implied volatility for different values ofλ2.Figure 3
Implied volatility for different values ofη1.Figure 4
Implied volatility for different values ofρ1.First, we can seen from Figure2 that an increase in λ2, the arrive intensity of the volatility jumps, increases the implied volatility. This implies that the BS model systematically undervalues all options. A similar argument holds for η1, the mean of the volatility jumps (see Figure 3). Second, the levels of the implied volatilities depend on the sign and magnitudes of correlation coefficient ρ1 between the volatility and asset (see Figure 4), whereas the levels of the implied volatilities at the lower strikes drop with increasing correlation ρ1 and ones at higher strikes go up. This feature is consistent with the models with stochastic volatility without jumps in volatility. Thus, the proposed model maintains some properties of the models being of stochastic volatility. We stress, however, that this model enables option prices to fit the forward price term structure exactly.


## 3.3. Probability Densities

We now go one step further and examine the distributional properties of the TSVIJ model using the FFT method. It is also instructive to compare the short-term and long-term probability density functions of TSVIJ, SVIJ, and Merton’s models. From (21), the FFT algorithm of qxT can be approximated by(37)qxT=1π∑j=0N−1e−2πi/Njvf0uj,for v=0,1,…,N−1. The results are plotted in Figures 5–8, which examine the values of the probability density function for a range of values of two parameters when all the others are kept constant and taken from Table 1.Figure 5
Probability density function with different correlationsρ1 and ρ2.Figure 6
Probability density function with different jump intensitiesλ2 and λ3.Figure 7
Probability density function with different jump meansη1 and η2.Figure 8
Comparison of probability density function in three different models.
(a)(b)Figure5 shows how the positive correlations of variance process with the stock return (which is calculated in equation (3)) create a fat-right tail and thin-left tail distribution of the log-return. Positive correlations lead to a high variance when the stock return increases, which fattens the right tail of the density function. This is consistent with those in Wong and Zhao [30].Figures6 and 7 provide the impacts of the volatility jump components, λ2,λ3,η1, and η2 in the TSVIJ model on probability densities of the log-return. Figure 6 shows that increasing the arrive intensities of volatility jumps, λ2 and λ3, will decrease the kurtosis only, but create a fat-right-tail and thin-left-tail. A similar argument holds for the means of the volatility jumps, η1 and η2 (see Figure 7).Figure8 provides comparisons of the short-term (T=0.5 year) and long-term (T=1 year) probability densities of the TSVIJ, SVIJ, and Merton’s models. Figure 8 shows that the probability density of the TSVIJ model is distinctly different; it has very high peak and skewness features compared to those of the SVIJ and Merton models for both maturities. Moreover, for both maturity density curves in Figure 8 still show that the density of the TSVIJ model attenuates heavy tails on both sides. Therefore, the TSVIJ model can correct the pricing bias in the BS model and has a potential application in practice.


## 3.4. Sensitivity Analysis

As said previously, the TSVIJ model provides more flexibility than the single-factor specification to model the volatility term structure and probability density of asset returns. This section considers the effects of the model parameters, including the correlation coefficientsρ1 and ρ2, the arrive intensities λ2 and λ3 of volatility jumps, and the means η1 and η2 of volatility jumps, on the option prices. We analyze the prices of T = 1 year call option only. Model parameter values are taken from Table 1 apart from ρ1,ρ2,λ2,λ3,η1, and η2. The numerical results are displayed in Table 3.Table 3
The impacts ofλ2,λ3, η1,η2, and ρ1,ρ2 on option prices.
Strike priceKρ1=−0.4,ρ2=−0.6ρ1=0.4,ρ2=0.6λ2=0.3λ2=1.3λ2=0.3λ2=1.3λ3=0.9λ3=1.9λ3=0.9λ3=1.9η1=0.02,η2=0.048027.409727.828427.275927.69828524.428724.893524.328924.79989021.707222.207821.649422.15779519.215019.748919.208019.741010016.922117.486016.974417.547910514.958715.530615.059715.641511013.109113.693613.269313.863811511.423212.014811.643712.245612010.057310.638710.317410.9092η1=0.2,η2=0.48029.442632.529329.341432.46208526.642729.978026.580029.95459024.059227.582924.040027.60439521.731525.479921.700025.408810019.531323.417019.620323.542210517.601521.535417.736321.702211015.814019.843616.004920.061611514.175618.285414.424018.555912012.779616.855913.065417.1603From Table3, we can see that (i) option prices are very sensible to η1 and η2. Since η1 and η2 indicate the means of the volatility jumps, this sensitivity can be considered as the sensitivity of option prices to their volatilities, which are verified by equation (2). A similar argument holds for the arrive intensities λ2 and λ3 of the volatility jumps. (ii) Option prices with different strike prices have different sensitivities to the correlations ρ1 and ρ2, whereas in-the-money options (K=80,85,90,95) decrease in values with increasing correlations ρ1 and ρ2, at-the-money options (K=100), and out-the-money options (K=105,110,115,120) go up. This finding is consistent with those analyzed by Zhang and Wang [6] and among others.


## 3.5. Calibration of the TSVIJ Model to Market Data

With a successful implementation of the option pricing formula (15) under the TSVIJ model, we can apply the above calculation method to calibrate the TSVIJ model and estimate the model parameters on a set of market data of standard European-style call options by minimizing the difference between market prices and model prices in a least squared error fit.The dataset used in this study consists of AAPL plain vanilla call option prices traded on the NASDAQ. These option prices are based on the Apple stock without paying dividends, having a value of $ 364.11 on July 5, 2020. We employ the market quotes over the period from July 10, 2020, to July 24, 2020, as the in-sample data to calibrate the model parameters, and use those of July 31, 2020, for the out-of-sample test. The riskless dollar interest rater=0.0175, in this period. To ensure sufficient liquidity in option trading, we disregard option trading volume lower than 100 in the sample data. Finally, the option sample contains 128 call options, with 4 trading days and maturies from 5 days to 26 days. The description of AAPL options data are reported in Table 4.Table 4
Numbers of AAPL option data.
MaturityMoneynessS/KITMATMOTM>1.010.99∼1.01<0.995 days1161312 days1422519 days433526 days2310The TSVIJ model has 17 parameters that need estimation, viz.,Θ=v1,0,κ1,θ1,σ1,η1,ρ1,λ2,v2,0,κ2,θ2,σ2,η2,ρ2,λ3,μs,σs,λ1. The estimation procedure is to minimize the least squared error over the parameter space Θ, i.e., we evaluate(38)minΘSqΘ=1M∑m=1MCS0,r,δ,Km,Tm;Θ−CMarketKm,Tm2+PenaltyΘ,Θ0,where M is the number of the sample data and CMarketKm,Tm and CS0,r,δ,Km,Tm;Θ denote the market price and the model price at time t=0 of a call option, respectively, with maturity Tm and strike price Km. The model price CS0,r,δ,Km,Tm;Θ depending on the parameter space Θ is taken as input and calculated by the FFT method above. As the penalty function, we choose the distance to the initial parameter vector Θ0 (see Mikhailov and Nögel [32]):(39)PenaltyΘ,Θ0=Θ−Θ02.We employ the Lsqnonlin function (least-squares nonlinear) in MATLAB to implement (38). To ensure the success of the local optimization scheme (38), we need effective choices of the initial parameters Θ0. Based on reasonable results from the literature, we set the initial parameters Θ0 used the empirical results studied by Duffie et al. [12]. In addition, we also define lower and upper bounds for the optimal parameters which avoid possible solutions that, while mathematically feasible, are not acceptable in finance sense in the calibration. The resulting calibration results are reported in Table 5.Table 5
Estimated risk-neutral parameters for five given models.
Model parametersTSVIJSVIJSVMertonBSv1,00.03510.04920.07340.06260.1011κ124.372621.91522.9990θ10.06360.06900.5005σ10.31040.15620.0613η10.26000.3073ρ1−0.5315−0.4821−0.5012λ20.51640.9979v2,00.1123κ282.3024θ20.0310σ20.3308η20.6327ρ20.5036λ31.0304μs−0.2839−0.1089−0.1138σs0.37950.20920.1978λ10.89631.96951.9943RMSE0.13060.14210.13090.31160.3625ARPE3.25303.54073.25993.78063.8520CPU time (s)732.5932268.204346.892131.52715.5643Note. RMSE=∑m=1MModel price−Market price2/M and ARPE=1/M∑m=1MModel  price−Market price/Market price.Table5 gives the estimated risk-neutral parameters and pricing performances for the TSVIJ, the SVIJ, SV, Merton, and the BS models using the in-sample data. For comparative purposes, we provide the root mean square error (RMSE), the average relative percentage error (ARPE), and the CPU time. Considering the RMSE and ARPE values, the TSVIJ model best fits the market prices in the in-sample data, followed by the SV, SVIJ, Merton, and the BS models. Taking a closer look at the estimates, the jump-related parameter λ1,λ2,λ3,μs,σs,η1,η2 in the TSVIJ model is significant, suggesting that modeling the asset dynamics and its variance processes cannot ignore the jump components.Using these results in Table5, the model predicted values and its comparison with market prices of the out-of sample data are shown in Table 6. As the table shows, the TSVIJ model provides a good fit for most traded options in a lower RMSE value.Table 6
Out-of-sample pricing results.
Option id.KOption priceTSVIJSVIJSVMertonBS132048.8048.625447.637645.231747.490345.2890236015.6515.658315.641915.664815.455714.6615336512.6512.664112.845113.615812.666612.11574367.511.4011.484611.583911.901011.408710.9635537010.2010.292510.415010.320710.24239.89046372.59.459.48739.338110.21559.16668.894673758.508.60778.35118.76548.18007.97428377.57.207.23167.45137.47697.28007.290593806.406.37676.63677.38386.46346.349410382.56.656.32995.90276.24425.72655.6390113854.955.09675.24516.16745.06494.9924123903.653.69714.13914.29293.94853.8765133953.003.24253.27673.48903.07422.9724144002.342.39552.61482.80882.40082.2510154051.801.84192.11272.23931.88931.6838RMSE0.01670.18391.22290.19471.0379


## 4. Conclusion

This paper proposes a multifactor variance processes model (TSVIJ model) for European option pricing, with jumps in both asset’s returns and variances, which is a reasonable and natural extension of the two-factor stochastic volatility model, Merton’s jump-diffusion model, and the affine jump-diffusion model. Based on the martingale method and Fourier transform technique, we also obtain closed-form solutions for the characteristic function and European call option. Using FFT, we obtain accurate and efficient numerical algorithm for pricing option within this proposed model. Comparisons of both the implied volatility curves and probability densities under alternative models show that the proposed model has better pricing performance. The analysis of the proposed model empirically shows its ability to handle pricing problem much better than the SVIJ, SV, Merton, and the BS models. Once the model parameters have been calibrated to fit market prices, the proposed model can be used to price other products that are not actively traded in the market. An analysis of this model also reveals that the volatility jump components have a strong impact on option prices. Although the independent arriving jumps in volatilities and jumps in returns are considered in the proposed model, we recognize that it could also be extended to jumps occurring simultaneously in both returns and variances by slightly modifying the framework of this model.



---

*Source: 1960121-2020-09-01.xml*