Fast Fourier Transform Based Power Option Pricing with Stochastic Interest Rate, Volatility, and Jump Intensity

Abstract

Firstly, we present a more general and realistic double-exponential jump model with stochastic volatility, interest rate, and jump intensity. Using Feynman-Kac formula, we obtain a partial integrodifferential equation (PIDE), with respect to the moment generating function of log underlying asset price, which exists an affine solution. Then, we employ the fast Fourier Transform (FFT) method to obtain the approximate numerical solution of a power option which is conveniently designed with different risks or prices. Finally, we find the FFT method to compute that our option price has better stability, higher accuracy, and faster speed, compared to Monte Carlo approach.

1. Introduction

In financial markets, option pricing problem becomes important for market practitioners ranging from individuals to financial institutes. The classic methodology for option pricing is developed by Black and Scholes [1]. In their model, the stock price follows a geometric Brownian motion and the option pricing formula in no arbitrage. However, a large number of empirical examples show that the initial assumption of their model generates high pricing deviations between theory and practice, for example the volatility clustering in [2] and the incompleteness of markets in [3]. Due to those drawbacks, some market practitioners begin to propose new underlying asset processes to fit empirical facts and develop new option types to control and reduce risks.

Double-exponential jump model with stochastic volatility, interest rate, and jump intensity could overcome these drawbacks of the Black-Scholes model in [2–8]. For example, Santa-Clara and Yan [4, 5] find that innovations to two risks, respectively, denoted by diffusion volatility process and jump intensity process in their model, have affected the expected return in the stock market. The empirical results show that those two processes are largely uncorrelated and do not support models that make jump intensity vary with the level of diffusive volatility. Based on the results of their theory, we believe that two components of risk are more conducive to the optimization of our model. An excellent contribution of our model is developing the model of Santa-Clara and Yan in [4, 5], combing stochastic volatility, interest rate, and mean-reverting jump intensity processes. Thus, our model better corresponds with the real market than other constant jump intensity models.

As for new option types, compared with other exotic options, power option with considerable operability is simpler and more convenient in practical financial market [9]. Unfortunately, research on power option pricing is literally nil according to our best knowledge. Therefore, in this paper, we define a power option and creatively study and explain how the power coefficient affects the valuation of option.

However, it is difficult to get an exact solution expression of option price. Several numerical methods have been proposed to solve it, including finite difference methods by Zhang and Wang in [10, 11], Monte Carlo simulation by Green in [12] and Glasserman in [13] and the FFT method by Carr and Madan in [14]. However, the finite difference method and Monte Carlo simulation are difficult to be applied in option pricing because they require substantially more computing time [14–17]. Compared to those approaches, the FFT is more effective and simpler to compute in [14–17]. Therefore, in this paper, we employ the FFT to compute the valuation of power potion.

In conclusion, the prominent contribution in our paper is modeling a more general and realistic model with stochastic interest rate, volatility, and jump intensity. Defining the power option, discussing the effects of power coefficient, and using the FFT to obtain its numerical results are another three contributions.

This paper proceeds as follows. In Section 2, we build our model which links double-exponential jump, stochastic interest rate, volatility, and jump intensity. The characteristic function is derived in Section 3. In Section 4, we discuss the approximate solutions for power option pricing by the FFT. Numerical examples and discussions of the power coefficient are shown in Section 5. Conclusions are in Section 6.

2. Our Model

Let (Ω,ℱ_𝑡,ℚ) be a probability space where ℱ_𝑡 is the filtration generated by the Brownian motion and jump process at time 𝑡, 0≤𝑡≤𝑇 and ℚ is a risk-neutral probability. The underlying asset price 𝑆(𝑡) at time 𝑡 is given by 𝑑𝑆(𝑡)/𝑆(𝑡)=(𝑟(𝑡)−𝜆(𝑡)𝑚)𝑑𝑡+𝑉(𝑡)𝑑𝑊^𝑠(𝑡)+𝑒^𝑦−1𝑑𝑁^𝑠(𝑡),𝑆(0)=𝑆₀(𝑡)>0,𝑑𝑟=𝜅_𝑟𝜃_𝑟−𝑟(𝑡)𝑑𝑡+𝜀_𝑟𝑟(𝑡)𝑑𝑊^𝑟(𝑡),𝑟(0)=𝑟₀(𝑡)>0,𝑑𝑉=𝜅_𝑣𝜃_𝑣−𝑉(𝑡)𝑑𝑡+𝜀_𝑣𝑉(𝑡)𝑑𝑊^𝑣(𝑡),𝑉(0)=𝑉₀(𝑡)>0,𝑑𝜆=𝜅_𝜆𝜃_𝜆−𝜆(𝑡)𝑑𝑡+𝜀_𝜆𝜆(𝑡)𝑑𝑊^𝜆(𝑡),𝜆(0)=𝜆₀>0, (1) where 𝑊^𝑠(𝑡), 𝑊^𝑟(𝑡), 𝑊^𝑣(𝑡), and 𝑊^𝜆(𝑡) are Brownian motions; the correlation coefficient of 𝑊^𝑠(𝑡) and 𝑊^𝑉(𝑡) is 𝜌>0, and other motions are independent on each other; 𝑁^𝑠(𝑡) is a Poisson process with stochastic intensity 𝜆(𝑡), and jump size 𝑦 is a random variable; 𝑚=𝐸^ℚ[𝑒^𝑦−1] is average jump amplitude; 𝑟(𝑡) is the stochastic interest rate and 𝑉(𝑡) is stochastic volatility of the underlying asset return; the parameters 𝜀_𝑟, 𝜀_𝑣, and 𝜀_𝜆 and the mean-reverting rates 𝜅_𝑟, 𝜅_𝑣, and 𝜅_𝜆 are positive constants; the constants 𝜃_𝑟, 𝜃_𝑣, and 𝜃_𝜆 are the long-term interest rate, volatility, and jump intensity, respectively.

We suppose that jump size 𝑦 has an asymmetric double exponential distribution 𝑓(𝑦), 𝑓𝑦1=𝑝𝜂_𝑢𝑒^{(−1/𝜂𝑢)𝑦}1_𝑦≥01+𝑞𝜂_𝑢𝑒^{(1/𝜂𝑑)𝑦}1_𝑦<0,0<𝜂_𝑢<1,𝜂_𝑑>0, (2) where 𝑝,𝑞≥0, 𝑝+𝑞=1 are probability of up-move jump and down-move jump, respectively. 1/𝜂_𝑢and 1/𝜂_𝑑 are mean of positive jumps and negative jumps, respectively. Therefore, 𝑚=𝑝/(1−𝜂_𝑢)+𝑞/(𝜂_𝑑+1)−1.

3. The Characteristic Function

Given the log asset price process in (1), it is possible to obtain the characteristic function 𝜙(⋅) of 𝑋(𝜏) at time 𝜏=𝑇−𝑡. Under the risk-neutral measure ℚ, firstly, we define an explicit expression of the moment generating function (MGF) 𝑀(Φ,𝑋,𝑟,𝑉,𝜆,𝜏) of 𝑋(𝜏) at time 𝜏, 𝑀(Φ,𝑥,𝑟,𝑉,𝜆,𝜏)=Ε^ℚ𝑒^Φ𝑋(𝜏)=𝑒^−𝑟𝜏Ε^ℚ𝑒^𝑟𝜏𝑒^Φ𝑋(𝜏), (3) and the complex-valued characteristic function is given by 𝜙(𝑢)=𝑀(𝑖𝑢).

From (3), we find that 𝑀(Φ,𝑋,𝑟,𝑉,𝜆,𝜏) can be interpreted as a contingent claim that pays off 𝑒^{𝑟𝜏+Φ𝑋(𝑡)} at time 𝜏. Therefore, it is possible to get the MGF 𝑀(Φ,𝑋,𝑟,𝑉,𝜆,𝜏) and then obtain the characteristic function 𝜙(𝑢) in following theorem.

Theorem 1. The characteristic function 𝜙(𝑢) of 𝑋(𝜏) under the measure ℚ is given by 𝜙(𝑢)=𝑒^{𝑖𝑢𝑋+𝑖𝑢(𝑟−𝑑)𝜏+𝐴(𝑢,𝜏)+𝐵(𝑢,𝜏)𝑉+𝐶(𝑢,𝜏)+𝐷(𝑢,𝜏)𝜆}, (4) where 𝐴()𝜅𝑢,𝜏=−^𝑣𝜃_𝑣𝜀²_𝑣𝜑₊𝜑𝜏+2ln⁻+𝜑₊𝑒^−𝜍𝜏,𝐵()2𝜍𝑢,𝜏=−𝑖𝑢+𝑢²1−𝑒^−𝜍𝜏𝜑₋+𝜑₊𝑒^−𝜍𝜏,𝐶()𝜅𝑢,𝜏=−^𝜆𝜃_𝜆𝜀²_𝜆𝜓₊𝜓𝜏+2ln⁻+𝜓₊𝑒^−𝜉𝜏,𝐷()(𝑢)2𝜉𝑢,𝜏=2Λ1−𝑒^−𝜉𝜏𝜓₋+𝜓₊𝑒^−𝜉𝜏,𝜑_±=∓𝜅_𝑣−𝑖𝑢𝜌𝜀_𝑣+𝜍,𝜍=𝜅_𝑣−𝑖𝑢𝜌𝜀_𝑣²+𝜀_𝑣²𝑖𝑢+𝑢²,𝜓_±=∓𝜅_𝜆+𝜉,𝜉=𝜅²_𝜆−2𝜀²_𝜆Λ(𝑢),Λ(𝑢)=𝑝1−𝑖𝑢𝜂_𝑢+𝑞1+𝑖𝑢𝜂_𝑑𝑝−1−𝑖𝑢1−𝜂_𝑢+𝑞1+𝜂_𝑑−1. (5)

Proof. From (1), Feynman-Kac formula gives the following PIDE for the characteristic function, −𝑀_𝜏+1𝑟−2𝑉−𝜆𝑚𝑀_𝑥+12𝑉𝑀_𝑥𝑥+𝜅_𝑣𝜃_𝑣−𝑉𝑀_𝑣+12𝜀²_𝑣𝑉𝑀_𝑣𝑣+𝜌𝜀𝑉𝑀_𝑥𝑣+𝜅_𝜆𝜃_𝜆−𝜆𝑀_𝜆+12𝜀²_𝜆𝜆𝑀_𝜆𝜆+𝜅_𝑟𝜃_𝑟−𝑟𝑀_𝑟+12𝜀²_𝑟𝑟𝑀_𝑟𝑟+𝜆^∞_−∞𝑀𝜔𝑦𝑀()𝑋+𝑦−𝑀𝑑𝑦=0,Φ,𝑋,𝑟,𝑉,𝜆,0=𝑒^Φ𝑋. (6) We conjecture that the solution of PDE (6) is 𝑀(Φ,𝑋,𝑟,𝑉,𝜆,𝜏)=𝑒^{𝑋Φ+𝑟𝜏Φ+𝐴(Φ,𝜏)+𝐵(Φ,𝜏)𝑉+𝐶(Φ,𝜏)+𝐷(Φ,𝜏)𝜆+𝐸(Φ,𝜏)+𝐹(Φ,𝜏)𝑟}, (7) with initial conditions 𝐴(Φ,0)=0, 𝐵(Φ,0)=0, 𝐶(Φ,0)=0, 𝐷(Φ,0)=0, 𝐸(Φ,0)=0, and 𝐹(Φ,0)=0.
We consider the integral term in (6): &𝜆^∞_−∞𝑀𝑋+𝑦−𝑀𝜔𝑦𝑑𝑦=𝜆𝑀(Φ,𝑋,𝑟,𝑉,𝜆,𝜏)𝑈(Φ), (8) where 𝑈(Φ)∶=∫^∞_−∞(𝑒^Φ𝑦−1)𝜔(𝑦)𝑑𝑦=𝑝/(1−Φ𝜂_𝑢)+𝑞/(1+Φ𝜂_𝑑)−1.
Substituting (7) and (8) into (6) yields −𝐴_𝜏(Φ,𝜏)−𝐵_𝜏(Φ,𝜏)𝑉−𝐶_𝜏(Φ,𝜏)−𝐷_𝜏(Φ,𝜏)𝜆−𝐸_𝜏(Φ,𝜏)−𝐹_𝜏()1Φ,𝜏𝑟−2+1𝑉Φ−𝜆𝑚Φ2𝑉Φ²+𝜅_𝑣𝜃_𝑉𝐵(Φ,𝜏)−𝜅_𝑣()+1𝑉𝐵Φ,𝜏2𝜀²_𝑣𝑉𝐵²(Φ,𝜏)+𝜌𝜀𝑉Φ𝐵(Φ,𝜏)+𝜅_𝑟𝜃_𝑟𝐹(Φ,𝜏)−𝜅_𝑟()+1𝑟𝐹Φ,𝜏2𝜀²_𝑟𝑟𝐹²(Φ,𝜏)+𝜅_𝜆𝜃_𝜆𝐷(Φ,𝜏)−𝜅_𝜆()+1𝜆𝐷Φ,𝜏2𝜀²_𝜆𝜆𝐷²()(Φ)Φ,𝜏+𝜆𝑈=0. (9) If we denote Λ(Φ)∶=𝑈(Φ)−𝑚Φ=𝑝/(1−Φ𝜂_𝑢)+𝑞/(1+Φ𝜂_𝑑)−1−Φ(𝑝/(1−𝜂_𝑢)+𝑞/(1+𝜂_𝑑)−1), then (9) leads to a system of seven ordinary differential equations: 𝐴_𝜏(Φ,𝜏)+𝐶_𝜏(Φ,𝜏)+𝐸_𝜏(Φ,𝜏)=𝜅_𝑣𝜃_𝑣𝐵(Φ,𝜏)+𝜅_𝜆𝜃_𝜆𝐷(Φ,𝜏)+𝜅_𝑟𝜃_𝑟𝐹(),𝐵Φ,𝜏_𝜏()1Φ,𝜏=−2Φ−Φ²−𝜅_𝑣()+1−𝜌𝜀Φ×𝐵Φ,𝜏2𝜀²𝐵²(),𝐷Φ,𝜏_𝜏(Φ,𝜏)=−𝜅_𝜆𝐷()+1Φ,𝜏2𝜀²_𝜆𝐷²()(Φ),𝐹Φ,𝜏+Λ_𝜏(Φ,𝜏)=−𝑘_𝑟𝐹()+1Φ,𝜏2𝜀²_𝑟𝐹²(),𝐴Φ,𝜏_𝜏(Φ,𝜏)=𝑘_𝑣𝜃_𝑣𝐵(),𝐶Φ,𝜏_𝜏(Φ,𝜏)=𝑘_𝜆𝜃_𝜆𝐷(),𝐸Φ,𝜏_𝜏(Φ,𝜏)=𝑘_𝑟𝜃_𝑟𝐹(Φ,𝜏). (10)(11)(12)(13)(14)(15)(16)
We consider the second ODE (11), in which the general solution can be derived for the Riccati equation. Making the substitution, 𝐵(Φ,𝜏)=−2𝑂^′(𝜏)𝜀²𝑂(𝜏). (17) Hence, 𝑂^′′(𝜏)+𝑘_𝑣−𝜌𝜀Φ𝑂^′(𝜏)−𝜀²4Φ−Φ²𝑂(𝜏)=0. (18) The ODE (11) has a general solution of the form 𝑂(𝜏)=𝐶₁𝑒^{(−1/2)𝜑−𝜏}+𝐶₂𝑒^{(1/2)𝜑+𝜏}, (19) where 𝜑_±=∓𝑘_𝑣−𝜌𝜀Φ+𝜍,𝜍=𝑘_𝑣−𝜌𝜀Φ²+𝜀²Φ−Φ². (20)
According to (15) and initial condition 𝐶(Φ,0)=0, 𝐶₁ and 𝐶₂ are two constants to be determined from boundary conditions, 𝑂(0)=𝐶₁+𝐶₂,𝑂^′(0)1=−2𝜑₋𝐶₁+12𝜑₊𝐶₂. (21) We obtain that the solutions of (21) is 𝐶₁=𝜑₊𝑂(0),𝐶2𝜍₂=𝜑₋𝑂(0).2𝜍 (22) Therefore, the exact solution of (11) is 𝐵(Φ,𝜏)=−Φ−Φ²1−𝑒^−𝜍𝜏𝜑₋+𝜑₊𝑒^−𝜍𝜏. (23) Substituting (23) into (14) and using initial conditions 𝐴(Φ,0)=0, we have 𝐴()𝜅Φ,𝜏=−^𝑣𝜃_𝑣𝜀²𝜑₊𝜑𝜏+2ln⁻+𝜑₊𝑒^−𝜍𝜏.2𝜍 (24) The OEDs for 𝐶(Φ,𝜏), 𝐷(Φ,𝜏), 𝐸(Φ,𝜏), and 𝐹(Φ,𝜏) are solved by analogy. Hence 𝐶()𝜅Φ,𝜏=−^𝜆𝜃_𝜆𝜀²𝜓₊𝜓𝜏+2ln⁻+𝜓₊𝑒^−𝜉𝜏,𝐷()(Φ)2𝜉Φ,𝜏=2Λ1−𝑒^−𝜉𝜏𝜓₋+𝜓₊𝑒^−𝜉𝜏,𝐸()𝐹()𝜓Φ,𝜏=0,Φ,𝜏=0,_±=∓𝜅_𝜆+𝜉,𝜉=𝜅²_𝜆−2𝜀²_𝜆Λ(Φ). (25)
Using 𝜙(𝑢)=𝑀(𝑖𝑢), it is easy to get Theorem 1.

4. Power Option Pricing Using the Fast Fourier Transform

4.1. Power Option

Under the risk-neutral measure ℚ, the price of the power call option 𝐶(𝑆,𝑟,𝑉,𝜆,𝑡) at time 𝑡 with strike price 𝐾 and maturity date 𝑇 is given by 𝐶(𝑆,𝑟,𝑉,𝜆,𝑡)=𝐸^ℚ𝑒^{−∫𝑇𝑡𝑟𝑢𝑑𝑢}𝑆^𝛽−𝐾⁺. (26)

As a European option, the valuation of power at initial time is 𝑐𝐾,𝑇,𝛽∶=𝐶(𝑆,𝑟,𝑉,𝜆,0)=𝐸^ℚ𝑒^{−∫𝑇0𝑟𝑢𝑑𝑢}𝑆^𝛽−𝐾⁺. (27)

Remark 2. In particularly, when 𝛽=1, the power option degenerates into the vanilla option.

4.2. Fast Fourier Transform

Let 𝑘 denote the log of the strike price 𝐾, and let 𝑐(𝑘) be the desired value of an initial power call option with strike price 𝑒^𝑘. The risk-neutral density of the price 𝑋(𝑇) at time 𝑇 is 𝑞_𝑇(𝑋). From (27), the initial value of power call option is related to the risk-neutral density 𝑞_𝑇(𝑋) by 𝑐(𝑘)=𝑒^{−∫𝑇0𝑟𝑢𝑑𝑢∞}_𝑘𝑒^{𝛽𝑋(𝑇)}−𝑒^𝑘𝑞_𝑇(𝑋)𝑑𝑋. (28)

Here we employ the price of a zero coupon bond 𝑄(𝑇) at time 𝑇 to express the stochastic discounted process 𝑒^{−∫𝑇0𝑟𝑢𝑑𝑢}𝑄(𝑇)=𝑃(𝑇)𝑒^{−𝐻(𝑇)𝑟0}, (29) where 𝑃(𝑇)=2𝑚𝑒^{(𝑇/2)(𝜅𝑟+𝑚)}𝜅_𝑟+𝑚𝑒^𝑚𝑇−1+2𝑚^{2𝜅𝑟𝜃𝑟/𝜀2𝑟},𝐻(𝑇)=2𝑒^𝑚𝑇−1𝜅_𝑟+𝑚𝑒^𝑚𝑇,−1+2𝑚𝑚=𝜅²_𝑟+2𝜀²_𝑟. (30) Then (28) can be written as 𝑐(𝑘)=𝑄(𝑇)^∞_𝑘𝑒^{𝛽𝑋(𝑇)}−𝑒^𝑘𝑞_𝑇(𝑋)𝑑𝑋. (31)

Note that 𝑐(𝑘) tends to 𝑆^𝛽₀ as 𝑘 tends to −∞, and hence the call pricing function (31) is not square integrable. According to Wong and Lo [15], we now consider the modified call price function 𝑐(𝑘) defined by 𝑐(𝑘)≡𝑒^𝛼𝑘𝑐(𝑘),for𝛼>0. (32) The Fourier transform of 𝑐(𝑘) is defined by 𝜛_𝑇(𝑢)=^∞_−∞𝑒^𝑖𝑢𝑘𝑐(𝑘)=𝑑𝑘_∞−∞𝑒^𝑖𝑢𝑘𝑒^𝛼𝑘𝑐(𝑘)=𝑑𝑘_∞−∞𝑒^𝑖𝑢𝑘𝑒^𝛼𝑘𝑄(𝑇)^∞_𝑘𝑒^𝛽𝑋−𝑒^𝑘𝑞_𝑇(𝑋)(𝑇)𝑑𝑋𝑑𝑘=𝑄^∞_−∞𝑞_𝑇(𝑋)^𝑋_−∞𝑒^{𝛽𝑋+𝛼𝑘}−𝑒^(1+𝛼)𝑘𝑒^𝑖𝑢𝑘(𝑇)𝑑𝑘𝑑𝑋=𝑄^∞_−∞𝑒^−𝑟𝑇𝑞_𝑇(𝑋𝑒)^{(𝛽+𝛼+𝑖𝑢)𝑋}−𝑒𝛼+𝑖𝑢^{(1+𝛼+𝑖𝑢)𝑋}(𝑇𝜙1+𝛼+𝑖𝑢𝑑𝑋=𝑄)^𝑇𝑢−𝛽+𝛼𝑖−𝜙𝛼+𝑖𝑢^𝑇(𝑢−(1+𝛼)𝑖).1+𝛼+𝑖𝑢 (33) Then the inverse transform of 𝜛_𝑇(𝑢) is given by 𝑐(𝑘)=1𝜋^∞₀𝑒^{−𝑖𝑢𝑘}𝜛_𝑇(𝑢)𝑑𝑢. (34) Then 𝑐(𝑘)=𝑒^−𝛼𝑘𝜋^∞₀𝑒^{−𝑖𝑢𝑘}𝜛_𝑇(𝑢)𝑑𝑢. (35) Using the Trapezoid rule for the integral in (35), the value of 𝑐(𝑘) is approximated as 𝑐(𝑘)≈𝑒^−𝛼𝑘𝜋^𝑁_𝑗=0𝑒^{−𝑖𝑢𝑘}𝜛_𝑇𝑢_𝑗𝜈, (36) where 𝑢_𝑗=𝜈(𝑗−1)

The FFT returns 𝑁 values of 𝑘 and for a regular spacing size of 𝑧 where 𝑁 is a power of 2. The value for 𝑘 is 𝑘_𝑢=−𝑏+𝑧(𝑎−1),for𝑎=1,…,𝑁, (37) which corresponds to log strike prices ranging from −𝑏 to 𝑏, where 𝑏=𝑁𝑧/2.

Substituting (37) into (36) yields 𝑐𝑘_𝑢≈𝑒^{−𝛼𝑘𝑢}𝜋^𝑁_𝑗=0𝑒^{−𝑖𝑧𝜈(𝑗−1)(𝑎−1)}𝑒^{𝑖𝑏𝑢𝑗}𝜛_𝑇𝑢_𝑗𝜈; (38) here we denote 𝜈𝑧=2𝜋/𝑁.

With Simpson's rule weightings, the value of European call option is 𝑐𝑘_𝑢≈𝑒^{−𝛼𝑘𝑢}𝜋^𝑁_𝑗=1𝑒^{−𝑖(2𝜋/𝑁)(𝑗−1)(𝑎−1)}𝑒^{𝑖𝑏𝑢𝑗}𝜛_𝑇×𝑢_𝑗𝜈33+(−1)^𝑗−ℓ_𝑗−1, (39) where ℓ_𝑛=1,𝑛=0,0,otherwise. (40)

5. Numerical Examples and Discussion

In this section, we present several numerical experiments to compare numerical differences between the FFT and Monte Carlo simulation and discuss the effects of power coefficient. Now, we perform some Monte Carlo simulations for the endowment dynamics by using the discrete scheme of system (1) as follows: ()(𝑡)+(𝑡)1ln𝑆𝑡+Δ𝑡=ln𝑆𝑟−𝑑−𝜆𝑚−2𝑉(𝑡+)Δ𝑡𝑉(𝑡)𝑑𝑊^𝑠(𝑡)+𝑒^𝑦(𝑡𝑁()(𝑡,𝑉()(𝑡)−1ln𝑆)(𝑡+Δ𝑡−𝑁))𝑡+Δ𝑡=𝑉+𝜅_𝑣𝜃_𝑣−𝑉(𝑡)Δ𝑡+𝜀_𝑣𝑉(𝑡)𝑑𝑊^𝑣(𝑡),𝜆(𝑡+Δ𝑡)=𝜆(𝑡)+𝜅_𝜆𝜃_𝜆−𝜆(𝑡)𝑑𝑡+𝜀_𝜆𝜆(𝑡)𝑑𝑊^𝜆(𝑡),𝑟(𝑡+Δ𝑡)=𝑟(𝑡)+𝜅_𝑟𝜃_𝑟−𝑉(𝑡)Δ𝑡+𝜀_𝑟𝑟(𝑡)𝑑𝑊^𝑟(𝑡), (41) where 𝑑𝑊^𝑠(𝑡)=𝜀₁√Δ𝑡, 𝑑𝑊^𝑣(𝑡)=𝜌𝜀₁√Δ𝑡, 𝑑𝑊^𝜆(𝑡)=𝜀₂√Δ𝑡, and 𝑑𝑊^𝑟(𝑡)=𝜀₃√Δ𝑡; Δ𝑡 is the time interval; 𝜀₁, 𝜀₂, and 𝜀₃ are samples from the standard normal distribution with correlation coefficient 0; 𝑁(𝑡)=𝑗, if 𝑆_𝑗≤𝑡≤𝑆_𝑗+1, 𝑆_𝑗=∑^𝑗_𝑘=1𝜏_𝑘, 𝜏_𝑘 is the 𝑘th jump interval which obeys an exponential distribution with parameter 𝜆.

The option prices are calculated by Monte Carlo simulation using 𝑉(𝑆(𝑡),𝐾,𝑉(𝑡),𝜆(𝑡),𝑟(𝑡),𝑡)=𝐸^ℚ𝑆^𝜀_𝑇+𝑆_𝑇^𝜀−2−𝐾⁺∣𝐹_𝑡. (42)

In our numerical examples by the FFT, we take 𝑁=4098 and 𝜈=600/𝑁 which lead to the log strike spacing of 𝑧=𝜋/300 and the damping coefficient 𝛼=1.18. Other parameters are denoted as 𝑟₀=0.05, 𝜅_𝑣=0.3, 𝜃_𝑣=0.15, 𝜀_𝑣=0.1, 𝜅_𝜆=5, 𝜃_𝜆=0.6, 𝜀_𝜆=0.3, 𝜅_𝑟=0.5, 𝜃_𝑟=0.05, 𝜀_𝑟=0.1, 𝜌=−0.25, 𝜂_𝑢=0.03, 𝜂_𝑑=0.13, 𝑝=0.4, 𝑞=0.6, 𝑆₀=100, 𝑉₀=0.15, 𝜆₀=3, 𝑇=0.5, 𝜅=100, and Δ𝑡=0.01. The results are presented in Tables 1, 2, and 3.

Table 1 

Power option prices in 𝛽=0.98 case: FFT versus Monte Carlo.

Table 2 

Power option prices in 𝛽=1 case: FFT versus Monte Carlo.

Table 3 

Power option prices in 𝛽=1.02 case: FFT versus Monte Carlo.

Numerical examples presented that the FFT is faster than Monte Carlo method with different parameters. To different strike prices, the FFT takes around approximate 340 seconds to produce 4098 option prices in three numerical experiments. However, Monte Carlo technique takes around 1260 seconds. On the level of the accuracy, comparing with Monte Carlo, we show that the percentage differences of the FFT are less than 0.1% in Tables 1, 2, and 3. Therefore, the FFT is more accurate and faster than Monte Carlo.

To examine effects of the power coefficient 𝛽, we calculate valuations with 𝛽=0.98, 𝛽=1, and 𝛽=1.02. Firstly, we can see that with the 𝛽 rising, the gaps between 𝑆^𝛽 and 𝐾 become greater in Figures 1, 2, and 3. Secondly, accordoing to different 𝛽, we find that valuations of option are increasing with the growth of 𝛽. It is because the fluctuations and risks of the underlying asset price are enhancing with the growth of 𝛽, and in order to compensate risk, the option seller need gets higher price of option. It means that the seller can design different power options with different risks or prices by defining different the power coefficients. Therefore, power option is one of convenient options which are designed easily.

Figure 1 

Power option prices in 𝛽=0.98 case.

Figure 2 

Power option prices in 𝛽=1 case.

Figure 3 

Power option prices in 𝛽=1.02 case.

6. Conclusion

In this paper, we study a more general and realistic double exponential-jump model which links the stochastic interest rate, volatility, and jump intensity. If we decide the constant interest rate, volatility, or jump intensity, our model can be degenerated into many previous models, for example, Black and Scholes’ model in [1], Santa-Clara and Yan’ model in [4, 5], Carr and Madan’s model in [14], Zhang and Wang’s model in [10, 11], and Zhang and Wang’s model in [17].

Moreover, we obtain the characteristic function of the log asset price, which is important to not only take the FFT in [14–17] but also derive the option pricing formula, for example, in [18]. Then, in our paper, a power option is defined and its valuations are computed with the FFT by characteristic function. Through numerical examples, comparing the results produced from the FFT and the Monte Carlo simulations of power option valuation, we find that the FFT is approximate with Monte Carlo approach but has better stability, higher accuracy, and faster speed in pricing power option. In addition, we show that the seller can design different power options with different risks or prices by defining different the power coefficients.

Acknowledgment

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