Option pricing using the fast Fourier transform under the double exponential jump model with stochastic volatility and stochastic intensity

Journal of Computational and Applied Mathematics 263 (2014) 152–159

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Option pricing using the fast Fourier transform under the double exponential jump model with stochastic volatility and stochastic intensity✩ Jiexiang Huang ∗ , Wenli Zhu, Xinfeng Ruan School of Economic Mathematics, Southwestern University of Finance and Economics, Chengdu 611130, PR China

article

info

Article history: Received 21 August 2013 Received in revised form 14 November 2013 Keywords: Option pricing Fast Fourier transform Double exponential jump Stochastic volatility Stochastic intensity

abstract This paper is based on the FFT (Fast Fourier Transform) approach for the valuation of options when the underlying asset follows the double exponential jump process with stochastic volatility and stochastic intensity. Our model captures three terms structure of stock prices, the market implied volatility smile, and jump behavior. Via the FFT method, numerical examples using European call options show effectiveness of the proposed model. Meanwhile, numerical results prove that the FFT approach is considerably correct, fast and competent. Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved.

1. Introduction The Black–Scholes model and its extensions comprise one of the major developments in modern finance. The methodology for pricing options is developed by Black and Scholes [1]. Although the pricing in the Black–Scholes model can be done very accurately, these prices do not capture the volatility pattern observed from traded option prices. A number of papers are proposed for volatility dynamics, such as Hull and White [2], and Stein and Stein [3], but the most influential paper in the square root model is developed by Heston [4]. However, Heston’s model has a serious disadvantage: the model does not consider the rare events such as financial crisis. To get more realistic models, Bates [5] proposes the addition of jump in his model. In Bates’ framework, it is hard to disentangle the diffusion and jump risks since they are all driven by a single state variable, the diffusive volatility. Duffie, Pan and Singleton [6] synthesize and significantly extend the literature on affine asset pricing models by deriving a closed-form expression for an ‘‘extended transform’’ of an affine jump diffusion process, and highlight the impact on option ‘smirks’ of the joint distribution of jumps in volatility and jumps in the underlying asset price. Santa-Clara and Yan [7] find that the innovations to the two risks, respectively denoted by diffusion volatility process and jump intensity process in their model, have affected the expected return in the stock market. When they calibrate the model to the Standard & Poor’s 500 index option prices from the beginning of 1996 to the end of 2002, they obtain time series of the implied diffusive volatility and jump intensity. The two components of risk vary substantially over time and show a high degree of persistence. In their model, it is the first model in which the jump intensity follows explicitly its own stochastic process. The empirical results show that the estimated correlation between the increments of the diffusive volatility and jump intensity is quite low. This is an evidence that the two processes are largely uncorrelated and do not support models

✩ This work is supported by the Fundamental Research Funds for the Central Universities (JBK130401).

∗

Corresponding author. Tel.: +86 18215505829. E-mail address: [email protected] (J. Huang).

0377-0427/$ – see front matter Crown Copyright © 2013 Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.009

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

153

Fig. 1. The intensity defined by different levels of extreme returns.

that make jump intensity vary with the level of diffusive volatility. They also obtain the risk-adjusted dynamics of the stock, volatility, and jump intensity processes and use them to price European options. Using stock return data of 10 years, Chang, Fuh and Lin [8] confirm the existence of switches in jump intensity and show that modeling the dynamic nature of jump dependence effectively captures overall changing volatility, supporting the use of this class of models to demonstrate stochastic volatility. Following the same way of Chang, Fuh and Lin [8], we analyze extreme returns of the Standard & Poor’s 500 index from 1983 to 2013. Summary statistics are presented in Fig. 1 which shows the number of days for each year at different levels of returns exceeded, i.e., 2%, 3%, and 5%. We find that the jump intensity is consistently high or low and follows a similar mean-reverting process. Especially, in 2% threshold, jump intensity fluctuates around mean 30. Dates, empirically, demonstrate that the stochastic jump intensity model is better than the constant jump intensity model. Because of several previous studies and the display of the dates, we add the jump intensity which follows explicitly mean-reverting stochastic process to our model. An excellent contribution of our model is developing the model of Santa-Clara and Yan [7], which combines stochastic volatility, jump and mean-reverting jump intensity processes. Thus, our model better corresponds with the real market than the constant jump intensity model. As for the distribution of the jump size, Merton [9] proposed the jump diffusion where the logarithm of jump size is normally distributed. Although Merton’s model makes the asset pricing more realistic and explains the volatility smile in a stationary way, it cannot capture the leptokurtic feature of the return distribution. Another jump diffusion is studied by Kou [10], where the jump size has an asymmetric double exponential distribution. This double exponential jump diffusion model is able to reproduce the leptokurtic feature of the return distribution and the volatility smile observed in option prices. In addition, the empirical tests performed by Kou and Wang in [11] suggest that the double exponential jump diffusion model fits stock data better than the normal jump diffusion model. Because the double exponential distribution has both a high peak and heavy tails, it can be used to model both the overreaction (attributed to the heavy tails) and underreaction (attributed to the high peak) to outside news. Therefore, the double exponential jump diffusion is captured by our model, within stochastic volatility and stochastic intensity. The finite difference method and the Monte Carlo simulation are usually used to value the options. But both ways are difficult to be applied in option pricing because they require substantially more computing time than the FFT approach. A large amount of the current literature on option valuation has fully applied Fourier analysis to determine option prices, such as Scott [12], and Bakshi and Chen [13]. But those papers do not consider the computational power of the FFT approach, which is one of the most fundamental advances in computing. Carr and Madan [14] describe an approach for numerically determining option values, which is designed to use the FFT approach to value options efficiently. Their use of the FFT approach in the inversion stage permits real-time pricing, marking, and hedging using realistic models, even for books with thousands of options. They anticipate that advantages of the FFT approach are generic to the widely known improvements in computation speed attained by this algorithm and is not connected to the particular characteristic function or process they chose to analyze. A number of recent papers use the FFT approach for a solution to the problem of option valuation (see e.g. Wong and Lo [15], Pillay and OHara [16], Zhang and Wang [17]). A study closely related to this paper is that by Pillay and OHara [16], who present a general formula for valuing derivatives with mean reversion and stochastic volatility. Yet we focus on European options pricing by the FFT technique in the double exponential jumps setting with stochastic volatility and stochastic intensity. The model allows the covariance structure of the stocks to the intensity to be unrestricted, which proves to be important in the numerical analysis and empirical estimation. Numerical examples indicate that the formulas are easy to implement, and are accurate. Unfortunately, this approach has not been used for option pricing under double exponential jump process with stochastic volatility and stochastic intensity. Another contribution of our paper is employing the FFT technique to address our model and obtaining the solution of the characteristic function and the expression of the option pricing formula.

154

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

This paper proceeds as follows. In Section 2, we present a model setup under the double exponential jump model with stochastic volatility and stochastic intensity. The characteristic function is derived in Section 3. In Section 4, we discuss approximation solutions for European options pricing by FFT. The proposed algorithm estimation and numerical examples are provided in Section 5. The paper is concluded in Section 6. 2. The model Let (Ω , F , Q) be a probability space where Ft is the filtration generated by the Brownian motion and the jump process at time t , 0 ≤ t ≤ T and Q is a risk-neutral probability. The underlying asset price S (t ) at time t is given by

  V (t )dWs (t ) + (eJ − 1)dN (t ), + dS (t )/S (t ) = (r − d − λ(t )m)dt  dV (t ) = κv (θv − V (t ))dt + εv V (t )dWv (t ), V (0) = V0 > 0,  λ(0) = λ0 > 0, dλ(t ) = κλ (θλ − λ(t ))dt + ελ λ(t )dWλ (t ),

S (0) = S0 > 0, (1)

where d is the dividend rate, N (t ) is a Poisson process with stochastic intensity λ(t ), and jump size J is a random variable. m is the average jump amplitude where m = E Q [eJ − 1]. The stochastic volatility of the stock return is V (t ). The parameters εv and ελ , the mean-reverting rates κv and κλ are positive constants. The constants θv and θλ are the long-term volatility and intensity constants, respectively. Moreover, Ws (t ) and Wv (t ) are a pair of correlated Brownian motions with correlation coefficient ρ which is constant. The Brownian process Wλ (t ) is independent of Ws (t ) and Wv (t ). We suppose that the jump size J has an asymmetric double exponential distribution with the density f (J ), f (J ) = p

1 −1J 1 1J e ηu 1J ≥0 + q e ηd 1J ≤0 ,

ηu

ηd

0 < ηu < 1, ηd > 0,

where p, q ≤ 1, p + q = 1 are probability of up-move jump and probability of down-move jump respectively. η1 , η1 are u d p q mean of positive jumps and mean of negative jumps respectively. Therefore, m = 1−η + 1+η − 1. u

d

Below is a more detailed comparison between our model and alternative models. The Black–Scholes model assumes that under the martingale measure Q the price of the underlying asset follows a standard Brownian motion with constant drift and diffusion parameters. However, in practice, almost all option markets exhibit the volatility smile phenomenon which means that options with different maturities and strikes have different implied Black–Scholes volatilities. The double exponential jump process with stochastic volatility and stochastic intensity model attempts to improve empirical implications of the Black–Scholes model while still retaining its analytical tractability. In Merton’s paper (Merton [9]), variable eJ is log-normally distributed. Both the double exponential and normal jump-diffusion models can lead to the leptokurtic feature. However, it is shown by Kou [10] that the kurtosis from the double exponential jump-diffusion model is significantly more pronounced. The stochastic volatility models (see Heston [4]) are in good agreement with implied volatility surfaces, and implied smiles are quite stable and unchanging over time. Unfortunately, investors have to hedge stochastic volatility to replicate and price the option. Unlike a stock or a currency, volatility is not a traded variable with an observable price. More precisely, in stochastic volatility models, the kurtosis decreases as the sampling frequency increases, while in jump models the instantaneous jumps are independent of the sampling frequency. 3. Deriving the characteristic function Given the log asset price process in (1), it is possible to obtain the characteristic function φ(·) of X (τ ) := ln S (τ ) at time τ := T − t. Under the measure Q, we define an explicit expression for the moment generating function(MGF) M (Φ , X , V , λ, τ ) of X (τ ) at time τ , M (Φ , X , V , λ, τ ) = E Q [eΦ X (τ ) ] = e−r τ E Q [er τ eΦ X (τ ) ],

(2)

and the complex-valued characteristic function is given by φ(u) = M (iu). If no confusion is possible, we write M (Φ ) instead of M (Φ , X , V , λ, τ ). We find that M (Φ ) can be interpreted as a contingent claim that pays off er τ +Φ X (τ ) at time τ . Therefore, it is possible to get M (Φ ) and then obtain the characteristic function φ(u). From (1), the Feynman–Kac formula in (Duffie, Pan and Singleton [6]) gives the following PIDE for the characteristic function,

   1 1   −Mτ + r − d − λ(t )m − 2 V (t ) Mx + 2 VMxx + κv (θv − V )Mv   ∞ 1 2 1 2 + ε VM + ρε VM + κ (θ − λ) M + ε λ M + λ [M (X + J ) − M ]ω(J )dJ = 0,  vv xv λ λ λ λλ v λ   2 2 −∞  M (Φ , X , V , λ, τ ) = eΦ X .

(3)

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

155

We conjecture the solution of PIDE (3) is M (Φ ) = eX Φ +(r −d)τ Φ +A(Φ ,τ )+B(Φ ,τ )V +C (Φ ,τ )+D(Φ ,τ )λ ,

(4)

with initial conditions A(Φ , 0) = 0, B(Φ , 0) = 0, C (Φ , 0) = 0, and D(Φ , 0) = 0. We consider the integral term in (3),

λ

∞



[M (X + J ) − M ]ω(J )dJ = λ



=λ



∞

−∞ ∞

−∞

=λ

−∞  ∞

[E Q [eΦ (X +J ) ] − E Q [eΦ X ]]ω(J )dJ [E Q [eΦ X (eΦ J − 1)]]ω(J )dJ E Q [eΦ X ]E Q [eΦ J − 1]ω(J )dJ

−∞

= λM (Φ )U (Φ ),

(5)

where U (Φ ) := −∞ (eΦ J − 1)ω(J )dJ = 1−Φ η + 1+Φ η − 1. u d Substituting (4) into (3) yields

∞

p

q

− (r − d)Φ − Aτ (Φ , τ ) − Bτ (Φ , τ )V − Cτ (Φ , τ ) − Dτ (Φ , τ )λ + (r − d)Φ 1

1

2

2

1

− V Φ − λmΦ + V Φ 2 + κv θv B(Φ , τ ) − κv VB(Φ , τ ) + εv2 VB2 (Φ , τ ) 2

1

+ ρεV Φ Bτ (Φ , τ ) + κλ θλ D(Φ , τ ) − κλ λD(Φ , τ ) + ελ λD (Φ , τ ) + λU (Φ ) = 0. 2

2

2

If we denote Λ(Φ ) = U (Φ ) − mΦ = ordinary differential equations,

p 1−Φ ηu

+

q 1+Φ ηd

p − 1 − Φ ( 1−η + u

Aτ (Φ , τ ) + Cτ (Φ , τ ) = κv θv B(Φ , τ ) + κλ θλ D(Φ , τ ), 1

q 1+ηd

(6)

− 1), then (6) leads to a system of five (7)

1

Bτ (Φ , τ ) = − (Φ − Φ 2 ) − (κv − ρεv Φ )B(Φ , τ ) + εv2 B2 (Φ , τ ), 2 2 1 2 2 Dτ (Φ , τ ) = −κτ D(Φ , τ ) + ελ D (Φ , τ ) + Λ(Φ ), 2

(8) (9)

where Aτ (Φ , τ ) = κv θv B(Φ , τ ) and Cτ (Φ , τ ) = κλ θλ D(Φ , τ ). We consider the second ODE (8), which the general solution can be derived for the Riccati equation. Making the substitution B(Φ , τ ) = −

2O′ (τ )

εv2 O(τ )

.

(10)

Hence,

ε O′′ (τ ) + (κv − ρεv Φ )O′ (τ ) − v (Φ − Φ 2 )O(τ ) = 0. 4 2

(11)

The ODE (11) has a general solution of the form 1

1

O(τ ) = C1 e− 2 φ− τ + C2 e 2 φ+ τ ,

(12)

where

φ± =∓(κv − ρεv Φ ) + ς , ς = (κv − ρεv Φ )2 + εv2 (Φ − Φ 2 ).

(13)

According to (12) and initial condition C (Φ , 0) = 0, C1 and C2 are two constants to be determined from boundary conditions, O(0) = C1 + C2 , 1 1 O′ (0) = − φ− C1 + φ+ C2 . 2 2



(14)

156

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

We obtain the solution to (14) is C1 = B(Φ , τ ) = −

1 2 − 2 φ−

ε2

φ+ O(0) 2ς

φ+ O(0) − 1 φ− τ e 2 2ς

φ+ O(0) − 1 φ− τ e 2 2ς

= −(Φ − Φ 2 )

Therefore, the general solution for (10) is

+ 12 φ+ φ−2Oς(0) e 2 φ+ τ 1

+

1

= −(Φ − Φ 2 )

φ− O(0) . 2ς

and C2 =

φ− O(0) 1 φ+ τ e2 2ς 1

e − 2 φ− τ + e 2 φ+ τ 1

1

φ+ e− 2 φ− τ + φ− e 2 φ+ τ 1 − e−ςτ

φ− + φ+ e−ςτ

.

(15)

Substituting (15) together with (12) into Aτ (Φ , τ ) = kv θv B(Φ , τ ), we have A(Φ , τ ) = κv θv

τ



B(Φ , s)ds

0

=−

2κv θv



τ

O′ (s)

ds

ε O(s) 0 2κv θv = − 2 (ln O(s) |ss=τ =0 ) ε   2κv θv O(τ ) = − 2 ln ε O(0)   φ+ O(0) − 1 φ− τ φ− O(0) 1 φ+ τ 2 2 + e e 2κv θv 2ς 2ς  = − 2 ln  φ+ O(0) ε + φ− O(0) 2

2ς

κv θv =− 2 ε



2ς

φ− + φ+ e−ςτ φ+ τ + 2 ln 2ς 



.

(16)

The OEDs for C (Φ , τ ) and D(Φ , τ ) are solved by analogy. Hence,

ψ− + ψ+ e−ξ τ ψ+ τ + 2 ln 2ξ 1 − e−ξ τ D(Φ , τ ) = 2Λ(Φ ) , ψ− + ψ+ e−ξ τ ψ± =∓κλ + ξ ,

κλ θλ C (Φ , τ ) = − 2 ε

ξ=







, (17)

κλ2 − 2ελ2 Λ(Φ ).

We conclude that the characteristic function of the mean reverting process with stochastic volatility, jump and stochastic intensity is M (Φ ) = eX Φ +(r −d)τ Φ +A(Φ ,τ )+B(Φ ,τ )V +C (Φ ,τ )+D(Φ ,τ )λ , where

   φ− + φ+ e−ςτ φ+ τ + 2 ln , 2ς 1 − e−ςτ B(Φ , τ ) = −(Φ − Φ 2 ) , φ− + φ+ e−ςτ    κλ θλ ψ− + ψ+ e−ξ τ C (Φ , τ ) = − 2 ψ+ τ + 2 ln , 2ξ ελ A(Φ , τ ) = −

κv θv εv2

D(Φ , τ ) = 2Λ(Φ )

1 − e−ξ τ

, ψ− + ψ+ e−ξ τ φ± = ∓(κv − ρεv Φ ) + ς ,  ς = (κv − ρεv Φ )2 + εv2 (Φ − Φ 2 ),

ψ± = ∓κλ + ξ ,  ξ = κλ2 − 2ελ2 Λ(Φ ).

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

157

Using φ(u) = M (iu), we have

φ(u) = eiuX +iu(r −d)τ +A(u,τ )+B(u,τ )V +C (u,τ )+D(u,τ )λ ,

(18)

where

   φ− + φ+ e−ςτ φ+ τ + 2 ln , 2ς 1 − e−ςτ , B(u, τ ) = −(iu + u2 ) φ− + φ+ e−ςτ    κλ θλ ψ− + ψ+ e−ξ τ C ( u, τ ) = − 2 ψ+ τ + 2 ln , 2ξ ελ A(u, τ ) = −

κv θv εv2

1 − e−ξ τ

, ψ− + ψ+ e−ξ τ φ± = ∓(κv − iuρεv ) + ς ,  ς = (κv − iuρεv )2 + εv2 (iu + u2 ),

D(u, τ ) = 2Λ(u)

ψ± = ∓κλ + ξ ,  ξ = κλ2 − 2ελ2 Λ(u), Λ(u) =

p 1 − iuηu

+

q 1 + iuηd

 − 1 − iu

p 1 − ηu

+

q 1 + ηd



−1 .

4. European option pricing using the FFT Under the risk-neutral measure Q, the price of the European call option C (S , r , V , t ) at time t with strike price K and maturity date T is given by C (S , r , V , t , K ) = E Q [e−r (T −t ) (S − K )+ ]. Then valuation of the option at initial time is c (K ) := C (S , r , V , 0, K ) = E Q [e−rT (S − K )+ ]. Let k denote the log of the strike price K , and let c (k) := c (ek ) = c (K ) be the desired value of an initial European call option with strike ek . The risk-neutral density of the log price x := X (T ) is qT (x). Then the valuation of option c (k) is related to the risk-neutral density qT (x) by ∞



c (k) = e−rT

(ex − ek )qT (x)dx.

(19)

k

Note that c (k) tends to S (0) as k tends to −∞, and hence the call pricing function (19) is not square-integrable. To obtain a square-integrable function, according to Carr and Madan [14], we consider the modified call price function c¯ defined by c¯ (k) ≡ eα k c (k) for α > 0.

(20)

The Fourier transform of c¯ (k) is defined by

ϖT (u) =

∞



eiuk c¯ (k)dk

−∞ −rT

=

e φT (u − (α + 1)i) . α 2 + α − u2 + i(2α + 1)u

Then the inverse transform of ϖT (u) is given by c¯ (k) =

1

∞



π

e−iuk ϖT (u)du. 0

Then c (k) =

e−α k

π

∞



e−iuk ϖT (u)du. 0

(21)

158

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159 Table 1 Call option prices: FFT vs. Monte Carlo. Strike price

FFT

Monte Carlo

90.2830 92.1938 94.1451 96.1377 98.1724 100.2502 102.3720 104.5387 106.7512 107.8750 110.1581

10.4458 9.6731 8.9326 8.2254 7.5522 6.9135 6.3096 5.7407 5.2065 4.9525 4.4700

10.43973 9.67525 8.93118 8.21982 7.54877 6.91015 6.31072 5.74133 5.21243 4.96002 4.47683

%difference 0.058143

−0.02222 0.015899 0.055712 0.045438 0.048479 −0.01775 −0.01097 −0.11377 −0.15161 −0.15256

Using the Trapezoid rule for the integral in (21), the value of CT (k) is approximated as N e−α k  −iuk ϖT (uj )ν, e

c (k) ≈

π

(22)

j =0

where uj = ν(j − 1). The FFT returns N values of k and for a regular spacing size of z where N is a power of 2, the value for k is ku = −b + z (a − 1)

for a = 1, . . . , N ,

which corresponds to log strike prices ranging from −b to b, where b = Substituting (23) into (22) yields c (ku ) ≈

(23) Nz . 2

N e−α ku  −iz ν(j−1)(a−1) ibuj e e ϖT (uj )ν.

π

(24)

j =0

Applying the FFT, we note that ν z = 2Nπ . With Simpson’s rule weightings, the value of European call option is c (ku ) ≈

N ν e−α ku  −i 2π (j−1)(a−1) ibuj e ϖT (uj ) [3 + (−1)j − ℓj−1 ], e N π j =1 3

(25)

where ℓn is the Kronecker delta function that is unity for n = 0 and zero otherwise. Thus ℓn is equivalent to

ℓn =



1, 0,

n = 0; otherwise.

5. Numerical results In this section, we study the numerical difference between using the FFT approach and the Monte Carlo simulation method. Moreover, we compare the performance of the two techniques, such as the accuracy and speed of call option pricing. For our FFT approach, we use N = 4098 and ν = 600 in the FFT scheme, which obtain a log strike spacing of N π z = 300 . The coefficient in the transform of the modified call option prices is α = 1.18. Parameters are set as follows: r = 0.05, d = 0.05, κv = 0.3, θv = 0.6, εv = 0.1, κλ = 5, θλ = 0.6, ελ = 0.3, ρ = −0.25, ηu = 0.03, ηd = 0.13, p = 0.4, q = 0.6, S0 = 100, V0 = 0.15, λ0 = 3, K = 100, T = 0.5. The numerical results are displayed in Table 1. To compare the FFT approach, we also compute numerical results by the Monte Carlo method. Firstly, we get a discrete scheme of system (1) as follows:

     1   ln S (t + △t ) = ln S (t ) + r − d − λ(t )m − V (t ) △t + V (t )ε1 △t ,   2  + (eJ − 1) ln S (t )(N (t + △t )− N (t )),     V (t + △t ) = V (t ) + κv (θv − V (t ))△t + εv V (t )(ρε 1 − ρ 2 ε2 ) △t ,  1 +    λ(t + △t ) = λ(t ) + κλ (θλ − λ(t ))△t + ελ λ(t )dε3 △t ,

(26)

where t is the time interval, and ε1 , ε2 , ε3 are samples from the standard normal distribution with correlation coefficient 0. Then, we use the following pricing formula to get the option valuation, V (S (t ), K , V (t ), λ(t ), t ) = E Q [(S (T ) − K )+ | Ft ].

J. Huang et al. / Journal of Computational and Applied Mathematics 263 (2014) 152–159

159

Here we denote △t = 0.001. Every simulation trial involves calculating two values of the option price. The first value is calculated in the usual way; the second value is calculated by changing the sign of all the random samples from standard normal distributions. Finally, the results are presented in Table 1. From the numerical examples, we find two main results. First, the FFT approach is faster than the Monte Carlo method with different strike prices. The FFT approach takes around approximately 207 s to produce 4098 option prices corresponding to different strike prices. For the same calculation, the Monte Carlo technique takes around 1157 s for each option price. Second, comparing their pricing accuracy, if we consider the Monte Carlo simulation to be the benchmark, the absolute percentage differences of the FFT approach are less than 0.153% for all cases in Table 1. Therefore, numerical results proved that the FFT approach is more correct and efficient. 6. Conclusion In this paper, we derive the characteristic function and obtain a numerical solution of European options pricing under the double exponential jump model with stochastic volatility and stochastic intensity. According to the solution of the characteristic function in the risk-neutral density, we use the FFT approach to get the option price numerically and its time value. Finally, the results of numerical simulations show that the FFT approach is more accurate and efficient than the Monte Carlo method. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Polit. Econ. (1973) 637–654. J. Hull, A. White, The pricing of options on assets with stochastic volatilities, J. Finance 42 (2) (1987) 281–300. E.M. Stein, J.C. Stein, Stock price distributions with stochastic volatility: an analytic approach, Rev. Financ. Stud. 4 (4) (1991) 727–752. S.L. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financ. Stud. 6 (2) (1993) 327–343. D.S. Bates, Jumps and stochastic volatility: exchange rate processes implicit in deutsche mark options, Rev. Financ. Stud. 9 (1) (1996) 69–107. D. Duffie, J. Pan, K. Singleton, Transform analysis and asset pricing for affine jump-diffusions, Econometrica 68 (6) (2000) 1343–1376. P. Santa-Clara, S. Yan, Crashes, volatility, and the equity premium: lessons from S&P 500 options, Rev. Econ. Stat. 92 (2) (2010) 435–451. C. Chang, C.-D. Fuh, S.-K. Lin, A tale of two regimes: theory and empirical evidence for a Markov-modulated jump diffusion model of equity returns and derivative pricing implications, J. Banking Finance (2013). R.C. Merton, Option pricing when underlying stock returns are discontinuous, J. Financ. Econ. 3 (1) (1976) 125–144. S.G. Kou, A jump-diffusion model for option pricing, Manage. Sci. 48 (8) (2002) 1086–1101. S.G. Kou, H. Wang, Option pricing under a double exponential jump diffusion model, Manage. Sci. 50 (9) (2004) 1178–1192. L.O. Scott, Pricing stock options in a jump-diffusion model with stochastic volatility and interest rates: applications of Fourier inversion methods, Math. Finance 7 (4) (1997) 413–426. G. Bakshi, C. Cao, Z. Chen, Empirical performance of alternative option pricing models, J. Finance 52 (5) (1997) 2003–2049. P. Carr, D. Madan, Option valuation using the fast Fourier transform, J. Comput. Finance 2 (4) (1999) 61–73. H.Y. Wong, Y.W. Lo, Option pricing with mean reversion and stochastic volatility, European J. Oper. Res. 197 (1) (2009) 179–187. E. Pillay, J. O’Hara, FFT based option pricing under a mean reverting process with stochastic volatility and jumps, J. Comput. Appl. Math. 235 (12) (2011) 3378–3384. S. Zhang, L. Wang, Fast Fourier transform option pricing with stochastic interest rate, stochastic volatility and double jumps, Appl. Math. Comput. 219 (23) (2013) 10928–10933.