A Fast Fourier Transform Technique for Pricing European Options with Stochastic Volatility and Jump Risk

Abstract

We consider European options pricing with double jumps and stochastic volatility. We derived closed-form solutions for European call options in a double exponential jump-diffusion model with stochastic volatility (SVDEJD). We developed fast and accurate numerical solutions by using fast Fourier transform (FFT) technique. We compared the density of our model with those of other models, including the Black-Scholes model and the double exponential jump-diffusion model. At last, we analyzed several effects on option prices under the proposed model. Simulations show that the SVDEJD model is suitable for modelling the long-time real-market changes and stock returns are negatively correlated with volatility. The model and the proposed option pricing method are useful for empirical analysis of asset returns and managing the corporate credit risks.

1. Introduction

The classical Black-Scholes (BS) model [1] has long been known to result in systematically biased option valuation. By adding jumps to the archetypal price process with Gaussian innovations Merton [2] is able to partly explain the observed deviations from the benchmark model which are characterized by fat tail and excess kurtosis in the returns distribution. For an overview of “stylized facts” on asset returns see Cont [3]. Statistical properties of implied volatilities are summarized in Cont et al. [4]. In the sequel also other authors develop more realistic models, for example, the pure jump models of Eberlein and Keller [5], Madan et al. [6], and Duffie et al. [7], stochastic volatility models of Steven [8], and stochastic volatility model with normal jumps of Bates [9] and Keppo et al. [10]. The double exponential jump-diffusion (DEJD) model, recently proposed by Kou [11], generates a highly skewed and leptokurtic distribution and is capable of matching key features of stock and index returns. Moreover, the DEJD model leads to tractable pricing formulas for exotic and path dependent options [12]. Accordingly, the DEJD model has gained wide acceptance. However, the DEJD model cannot capture the volatility clustering effects, which can be captured by stochastic volatility models [13]. Jump-diffusion models and the stochastic volatility model complement each other: the stochastic volatility model can incorporate dependent structures better, while the DEJD model has better analytical tractability, especially for path-dependent options. Since allowing interest rates to be stochastic does not improve pricing performance any further [14], the model that combines stochastic volatility and double exponential jump-diffusion (SVDEJD) may be more reasonable.

In the BS setting, the probability measure has a well-known analytic form [15], but, under stochastic volatility, it can only be obtained numerically. Monte Carlo simulation and the finite difference method are usually used to value the options. But, the two techniques require substantially more computing time and thus are difficult to be applied in real option pricing. Recently, being fast, accurate, and easy to implement, Fourier transforms have been widely used in valuing financial derivatives, for example, Carr and Madan [16] propose Fourier transforms with respect to log-strike price; Duffie et al. [7] offer a comprehensive survey that the Fourier methods are applicable to a wide range of stochastic processes; Carr and Wu [17] apply the transforms to time-changed Lévy processes and the class of generalized affine models. Hurd and Zhou [18] express the spread option payoff in terms of the gamma function and FFT technique. For an overview of option pricing using Fourier transforms, see Schmelzle [19].

The current paper extends the study of option pricing under the DEJD model in three ways. First, we propose a model which combines the double jumps and stochastic volatility. Second, using the martingale method, Fourier transform formula, and Feynman-Kac theorem, we obtain a closed-form solution for European call options pricing under the proposed model. Third, we obtain fast and accurate numerical solutions for European call options pricing by FFT technique.

The rest of the paper is organized as follows. Section 2 develops the underlying pricing model. Section 3 derives a closed-form solution for European call options pricing under the proposed model. Section 4 provides approximation solutions for European call options pricing by FFT technique. Section 5 numerically compares the density of the solutions to the alternative models and analyzes several effects on potion prices. Section 6 concludes. Applied program codes in Matlab package are presented in the appendix.

2. The Model

We consider an arbitrage-free, frictionless financial market where only riskless asset 𝐵 and risky asset 𝑆 are traded continuously up to a fixed horizon date 𝑇. Let {Ω,ℱ,{ℱ_𝑡}_{0≤𝑡≤𝑇},𝑃} be a complete probability space with a filtration satisfying the usual conditions, that is, the filtration is continuous on the right and ℱ₀ contains all 𝑃-null sets. Suppose 𝑊(𝑡), 𝑊_𝑣(𝑡) are both standard Brownian motion, which is ℱ_𝑡 adapted, and 𝑊(𝑡) has correlation 𝜌 with 𝑊_𝑣(𝑡).

Let 𝑆(𝑡) represent the price for a stock or a stock portfolio. Generally, instantaneous variance of asset returns in financial markets shows randomness; thus, the continuous part of the price process, defined as 𝑆^𝑐(𝑡), is𝑑𝑆^𝑐(𝑡)=𝑟𝑆^𝑐(𝑡)𝑑𝑡+𝜎𝑉(𝑡)𝑆^𝑐(𝑡)𝑑𝑊(𝑡), (2.1) where 𝑟 is risk-free rate and 𝜎 is nonnegative constant, and suppose 𝑆^𝑐(0)=𝑠, which can be set equal to 1 without any loss of generality. The size of the diffusion component is determined by 𝑉(𝑡), which represents, absent of any jump occurring, the level of (stochastic) return variance attributable to diffusion variations. For tractability, let 𝑉(𝑡) follow a square-root process:𝑑𝑉(𝑡)=𝜃_𝑣−𝛼_𝑣𝑉(𝑡)𝑑𝑡+𝜎_𝑣𝑉(𝑡)𝑑𝑊_𝑣(𝑡), (2.2) where nonnegative constants 𝜃_𝑣,𝜃_𝑣/𝛼_𝑣, and 𝜎_𝑣, respectively, reflect the speed of adjustment, the long-run mean, and the variation coefficient of 𝑉(𝑡), and suppose 𝑉(0)=𝑉₀.

It has been suggested from extensive empirical studies that markets tend to have both overreaction and underreaction to various good or bad news. One may interpret the jump part of the model as the market response to outside news. Good or bad news arrives according to a Poisson process, and the asset price changes in response according to the jump size distribution. According to Kou [11], the jumps in the log-price are modeled as a sequence of i.i.d. nonnegative random variables that occur at times determined by an independent Poisson process 𝑁(𝑡) with constant intensity 𝜆>0 such that 𝑌=ln𝑈 has an asymmetric double exponential distribution with the density𝑓_𝑌𝑦=𝑝𝜂₁𝑒^−𝜂1𝑦𝟏_𝑦≥0+𝑞𝜂₂𝑒^𝜂2𝑦𝟏_y<0,𝜂₁>1,𝜂₂>0, (2.3) where 𝟏 denotes the indicator function, so 𝟏_𝑦≥0 equals 1 if 𝑦≥0, but 0 otherwise. 𝑝,𝑞≥0,𝑝+𝑞=1 are up-move jump and down-move jump, respectively. Except for 𝑊(𝑡) which has correlation with 𝑊_𝑣(𝑡), all sources of randomness, 𝑊(𝑡),𝑊_𝑣(𝑡),𝑁(𝑡),𝑌_𝑗, and 𝑁(𝑡), are assumed to be independent.

Because of jumps and stochastic volatility, the risk-neutral probability measure is not unique. Following Naik and Lee [20] and Kou [11], by using the rational expectations argument with a HARA-type utility function for the representative agent, one can choose a particular risk-neutral measure 𝑃^∗ so that the equilibrium price of an option is given by the expectation under this risk-neutral measure of the discounted option payoff. Throughout this paper, we assume that there exists a martingale probability measure 𝑃^∗ being equivalent to 𝑃. Let 𝑋(𝑡) be the sum of all the jumps which occur up to and including time 𝑡, 𝐽(𝑡)=exp[𝑋(𝑡)]−𝐸[exp(𝑋(𝑡))], we have𝐽(𝑡)=exp𝑋(𝑡)−𝜆𝑡𝑝𝜂₁𝜂₁+−1𝑞𝜂₂𝜂₂+1−1. (2.4) Obviously, 𝐽(𝑡) is a 𝑃-martingale. Finally, the price process 𝑆(𝑡) is defined as𝑆(𝑡)=𝑆^𝑐(𝑡)𝐽(𝑡). (2.5)

Remark 2.1. The model contains most existing models as special cases. For example, we obtain (1) the BS model by setting 𝜆=0 and 𝜃_𝑣=𝛼_𝑣=𝜎_𝑣=0; (2) the SV model by setting 𝜆=0; (3) the DEJD model by setting 𝜃_𝑣=𝛼_𝑣=𝜎_𝑣=0.

Let 𝑑𝑊(𝑡)=𝜌𝑑𝑊_𝑣(𝑡)+1−𝜌²𝑑𝑍(𝑡), where 𝑍(𝑡) is standard Brownian motion that is ℱ_𝑡 adapted, and independent of 𝑊_𝑣(𝑡), 𝑁(𝑡), and random variables 𝑈_𝑗. From Itŏ’s formula, we haveln𝑆(𝑡)=ln𝐽(𝑡)+ln𝑆^𝑐(𝑡)=𝑋(𝑡)−𝜆𝑡𝑝𝜂₁𝜂₁+−1𝑞𝜂₂𝜂₂++1−1+𝑟𝑡𝜌𝜎^𝑡₀𝑉(𝑡)𝑑𝑊_𝑣(𝑡)−12𝜎²𝜌^2𝑡₀𝑉(𝑡)+𝑑𝑡1−𝜌²𝜎^𝑡₀𝑉(𝑡)(𝑡)−1𝑑𝑍2𝜎²1−𝜌^2𝑡₀𝑉(𝑡)(𝑡)𝑑𝑡=𝑋−𝜆𝑡𝑝𝜂₁𝜂₁+−1𝑞𝜂₂𝜂₂+1−1+𝑟𝑡+𝜉_𝑡+𝜍_𝑡. (2.6)

3. A Closed-Form Solution of European Option Pricing

In this section, we derive closed-form solution of a European call option pricing under the SVDEJD model. For a European put option, we can obtain easily corresponding result by the put-call parity [1]. For this purpose, we need the following results.

Lemma 3.1. Supposing the variance process 𝑉(𝑡) follows (2.2) and 𝑠₁,𝑠₂ are any complex, one has 𝐸exp−𝑠₁^𝑇₀𝑉(𝑡)𝑑𝑡−𝑠₂𝑉(𝑡)=exp𝐴(𝑇)−𝐵(𝑇)𝑉₀, (3.1) where 𝐴(𝑇)=2𝜃_𝑣𝜎²_𝑣ln2𝛾𝑒^{(1/2)(𝛼𝑣−𝛾)𝑇}2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇,𝐵(𝑇)=1−𝑒^−𝛾𝑇2𝑠₁−𝛼_𝑣𝑠₂+𝛾𝑠₂1+𝑒^−𝛾𝑇2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇,𝛾=𝛼²_𝑣+2𝜎²_𝑣𝑠₁. (3.2)

Proof. Let 𝐹(𝑉,0,𝑇)=𝐸{exp[−𝑠₁∫^𝑇₀𝑉(𝑡)𝑑𝑡−𝑠₂𝑉(𝑡)]}. Because of the affine structure of the variance process (2.2), we obtain that 𝐹(𝑉,0,𝑇) has a solution of the following form: 𝐹(𝑉,0,𝑇)=exp𝐴(𝑇)−𝐵(𝑇)𝑉₀. (3.3) From the Feynman-Kac formula, 𝐹(𝑉,0,𝑇) is the solution of the following backward Parabolic partial differential equation with the Cauchy problem: 𝜕𝐹+𝜃𝜕𝑡^𝑣−𝛼_𝑣𝑉𝜕𝐹+1𝜕𝑉2𝜎²_𝑣𝑉𝜕²𝐹𝜕𝑉²−𝑠₁𝐹()𝑉𝐹=0,𝑉,0,0=exp−𝑠₂𝑉₀. (3.4) Putting (3.3) in (3.4), we have 𝐴_𝑡(𝑇)−𝜃_𝑣𝐵(𝑇)=0,𝐴(0)=0,−𝐵_𝑡(𝑇)+12𝜎²_𝑣𝐵²(𝑇)+𝛼_𝑣𝐵(𝑇)−𝑠₁=0,𝐵(0)=𝑠₂. (3.5) Solving (3.5), we can obtain the result of Lemma 3.1.

Lemma 3.2. Supposing the asset price 𝑆(𝑇) follows (2.6) and 𝑧 is any complex, one has 𝐸exp[−𝑟𝑇+𝑧ln𝑆(𝑇)]=exp(𝑧−1)𝑟𝑇+𝜆𝑇𝑝𝜂₁𝜂₁+−𝑧𝑞𝜂₂𝜂₂+𝑧−1−𝑧𝜆𝑇𝑝𝜂₁𝜂₁+−1𝑞𝜂₂𝜂₂+1−1−𝑧𝜌𝜎𝜎_𝑣𝑉(𝑇)+𝜃_𝑣𝑇+2𝜃_𝑣𝜎²_𝑣ln2𝛾𝑒^{(1/2)(𝛼𝑣−𝛾)𝑇}2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇−1−𝑒^−𝛾𝑇2𝑠₁−𝛼_𝑣𝑠₂+𝛾𝑠₂1+𝑒^−𝛾𝑇2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇𝑉₀, (3.6) where𝑠₁()𝑧1=−𝑧−12𝜎²1−𝜌²−𝑧𝜌𝜎𝜎_𝑣𝛼_𝑣−12𝜎²𝜌²,𝑠₂=−𝑧𝜌𝜎𝜎_𝑣. (3.7)

Proof. Let 𝜙(𝑧)=𝐸{exp[−𝑟𝑇+𝑧ln𝑆(𝑇)]}. Because 𝑁(𝑡) is independent of 𝑊(𝑡),𝑊_𝑣(𝑡), and𝑍(𝑡), we have 𝜙(𝑇)=𝑒^{(𝑧−1)𝑟𝑇}𝐸𝑒^{𝑧ln𝐽(𝑇)}𝐸𝑒^{𝑧(𝜉𝑇+𝜍𝑇)}=𝑒^{(𝑧−1)𝑟𝑇}𝐶(𝑇)𝐷(𝑇). (3.8) From (2.2) and (2.3), we have 𝐶(𝑇)=exp𝜆𝑇𝑝𝜂₁𝜂₁+−𝑧𝑞𝜂₂𝜂₂+𝑧−1−𝑧𝜆𝑇𝑝𝜂₁𝜂₁+−1q𝜂₂𝜂₂.+1−1 (3.9) Because 𝑊_𝑣(𝑡) is a standard Brownian motion, we have 𝐸𝜍_𝑇1=−2𝜎²1−𝜌^2𝑇₀𝑉(𝑇)𝑑𝑡,Var𝜍_𝑇=𝜎²1−𝜌^2𝑇₀𝑉(𝑇)𝑑𝑡. (3.10) Then, 𝐷(𝑇))𝑧1=𝐸exp(𝑧−12𝜎²1−𝜌^2𝑇₀𝑉(𝑡)𝑑𝑡+𝑧𝜉_𝑇=exp−𝑧𝜌𝜎𝜎_𝑣𝑉(𝑇)+𝜃_𝑣𝑇)𝑧1×𝐸exp(𝑧−12𝜎²1−𝜌²+𝑧𝜌𝜎𝜎_𝑣𝛼_𝑣−12𝜎²𝜌^2𝑇₀𝑉(𝑡)𝑑𝑡−𝑧𝜌𝜎𝜎_𝑣𝑉₀. (3.11) Let 𝑠₁=−(𝑧−1)𝑧(1/2)𝜎²(1−𝜌²)−𝑧((𝜌𝜎/𝜎_𝑣)𝛼_𝑣−(1/2)𝜎²𝜌²), and 𝑠₂=−𝑧(𝜌𝜎/𝜎_𝑣). From Lemma 3.1, we have 𝐷(𝑇)=exp−𝑧𝜌𝜎𝜎_𝑣V(𝑇)+𝜃_𝑣𝑇+2𝜃_𝑣𝜎²_𝑣ln2𝛾𝑒^{(1/2)(𝛼𝑣−𝛾)𝑇}2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇−1−𝑒^−𝛾𝑇2𝑠₁−𝛼_𝑣𝑠₂+𝛾𝑠₂1+𝑒^−𝛾𝑇2𝛾𝑒^−𝛾𝑇+𝛼_𝑣+𝛾+𝜎²_𝑣𝑠₂1−𝑒^−𝛾𝑇𝑉₀. (3.12) From (3.8), (3.9), (3.10) and (3.12), we can obtain the required Lemma 3.2.

Lemma 3.3. Suppose 𝜑(𝑢)=𝐸[exp(𝑖𝑢ln𝑆(𝑇)] is the characteristic function of ln𝑆(𝑇); then 𝜑(𝑢)=2𝛿2𝛿+(𝛼_𝑣−𝛿−𝑖𝑢𝜌𝜎𝜎_𝑣)(1−𝑒^−𝛿𝑇)^{2𝜃𝑣/𝜎2𝑣}(𝑡)+𝜃×exp𝑖𝑢ln𝑆^𝑣𝛼_𝑣−𝛿𝑇𝜎²_𝑣−𝑖𝑢𝜃_𝑣𝜎𝜌𝑇𝜎_𝑣+𝜆𝑇𝑝𝜂₁𝜂₁+−𝑖𝑢𝑞𝜂₂𝜂₂+𝑖𝑢−1−𝑖𝑢𝑝𝜂₁𝜂₁+−1𝑞𝜂₂𝜂₂+1−1+𝑖𝑢𝑟𝑇+𝜖𝑉₀, (3.13) where 𝛿=𝛼_𝑣−𝑖𝑢𝜌𝜎𝜎_𝑣²+𝑖𝑢(1−𝑖𝑢)𝜎²𝜎²_𝑣,()𝜎𝜖=𝑖𝑢𝑖𝑢−1²1−𝑒^−𝛿𝑇2𝛿+𝛼_𝑣−𝛿−𝑖𝑢𝜌𝜎𝜎_𝑣1−𝑒^−𝛿𝑇. (3.14)

Proof. Let 𝜙(𝑧)=𝐸{exp[−𝑟𝑇+𝑧ln𝑆(𝑇)]}. Because 𝜑(𝑢)((𝑇=𝐸[(𝑇=𝐸exp𝑖𝑢ln𝑆))exp−𝑟𝑇+𝑖𝑢ln𝑆)]𝐸(=𝜙()exp−𝑟𝑇)𝑖𝑢𝜙(0), (3.15) from Lemma 3.2, we can obtain the required Lemma 3.3.

Theorem 3.4. Let 𝑘 denote the log of the strike price 𝐾, 𝑥_𝑇=ln(𝑆(𝑇)), and 𝐶_𝑇(𝑘) the desired value of a 𝑇-maturity call option with strike exp(𝑘). Assume that, under 𝑃^∗, the underlying nondividend-paying stock price 𝑆(𝑡) and its components are given by (2.1)–(2.5), 𝜑(𝑢) is the characteristic function of 𝑥_𝑇, 𝑞(𝑥) is the density of 𝑥_𝑇; then the initial call value 𝐶_𝑇(𝑘) is written as 𝐶_𝑇(𝑘)=12𝑆(𝑡)−𝑒^−𝑟𝑇𝐾+1𝜋^∞₀𝑆(𝑡)ℜ𝑒^𝑖𝑢𝑘𝜑_𝑇(𝑢−𝑖)𝑖𝑢−𝑒^−𝑟𝑇𝑒𝐾ℜ^𝑖𝑢𝑘𝜑_𝑇(𝑢)𝑖𝑢𝑑𝑢, (3.16) where ℜ[⋅] represents real part.

Proof. From the risk-neutral theory, we have 𝐶_𝑇(𝑘)=𝐸𝑒^−𝑟𝑇(𝑆(𝑇)−𝐾)⁺=𝑒^−𝑟𝑇₀^+∞(𝑆(𝑇)−𝐾)⁺𝑞(𝑆(𝑇))𝑑𝑆(𝑇)=𝑒^−𝑟𝑇_𝑘^+∞𝑒^𝑥𝑇𝑞(𝑥)𝑑𝑥−𝑒^−𝑟𝑇𝐾_𝑘^+∞𝑞(𝑥)𝑑𝑥=𝑆Π₁−𝑒^−𝑟𝑇Π₂. (3.17) Introducing a change of measure from 𝑃^∗ to 𝑄^∗ by a Radon-Nikodym derivative, we get 𝑑𝑄^∗𝑑𝑃^∗=𝑒^𝑥𝑇𝐸[𝑒^𝑥𝑇]. (3.18) With this new measure 𝑄^∗, the Fourier transform of Π₁ is defined as 𝐸^𝑄∗𝑒^{𝑖𝑢𝑥𝑇}=𝜑(𝑢−𝑖)𝜑().−𝑖 (3.19) Because of the no-arbitrage condition, we can obtain Π₁=12+1𝜋^∞₀ℜ𝑒^{−𝑖𝑢𝑘}𝜑_𝑇(𝑢−𝑖)𝑖𝑢𝜑_𝑇(−𝑖)𝑑𝑢. (3.20) From the Fourier transform formula, the probability density for our model is given by 𝑞(𝑥)=1𝜋₀^+∞𝑒^{−𝑖𝑢𝑘}𝜑(𝑢)𝑑𝑢. (3.21) Hence, Π₂=_𝑘^+∞1𝜋₀^+∞𝑒^{−𝑖𝑢𝑘}𝜑(𝑢)𝑑𝑢𝑑𝑥. (3.22) Changing the order of integration, we have Π₂=12+1𝜋^∞₀ℜ𝑒^{−𝑖𝑢𝑘}𝜑_T(𝑢)𝑖𝑢𝑑𝑢. (3.23) From (3.17), (3.20), and (3.23), we can obtain the required Theorem 3.4.

Remark 3.5. In (3.16), 𝐶_𝑇(𝑘) tends to 𝑆₀ not zero as 𝑘 goes to −∞. Hence, 𝐶_𝑇(𝑘) is not 𝐿¹ (absolutely integrable) and a Fourier transform does not exist.

4. Fast Fourier Transform for European Option Pricing

Since the integrand in (3.16) is singular at the required evaluation point 𝑢=0, the FFT cannot be applied directly to evaluate the integrals we mentioned above. Therefore, instead of solving for the risk-neutral exercise probabilities of finishing in-the-Money (ITM), Carr and Madan [16] introduce a new technique with the key idea to calculate the Fourier transform of a modified call option price with respect to the logarithmic strike price. With this specification and a FFT routine, a whole range of option prices can be obtained within a single Fourier inversion. In this section, we develop the numerical solutions of the prices by using the idea of Carr and Madan [16].

4.1. Fourier Transform of ITM and at-the-Money (ATM) Option Prices

By introducing an exponential damping factor 𝑒^𝛼𝑘 with 𝛼>0, it is possible to make the integrand in (3.16) be square integrable. We modified the pricing function (3.16) by𝐶_𝑇(𝑘)=exp(−𝛼𝑘)𝜋^∞₀𝑒^{−𝑖𝑣𝑘}𝜓_𝑇(𝑣)𝑑𝑣, (4.1) where 𝜓_𝑇(𝑣)=𝑒^−𝑟𝑇𝜑_𝑇(𝑣−(𝛼+1)𝑖)/(𝛼²+𝛼−𝑣²+𝑖(2𝛼+1)𝑣).

This method is viable when 𝛼 is chosen in a way that the damped option price is well behaved. Damping the option price with 𝑒^𝛼𝑘 makes it integrable for the negative axis 𝑘<0. On the other hand, for 𝑘>0 the option prices increase by the exponential 𝑒^𝛼𝑘, which influences the integrability for the positive axis. A sufficient condition of 𝑐_𝑇(𝑘) to be integrable for both sides (square integrability) is given by 𝜓(0) being finite, that is,𝜓(0)=𝑒^−𝑟𝑇𝜑_𝑇(−(𝛼+1)𝑖)𝛼²+𝛼<∞. (4.2) Thus we need 𝜑_𝑇(−(𝛼+1)𝑖)<∞, which is equivalent to𝐸^𝑄𝑆(𝑇)^1+𝛼<∞. (4.3) Therefore, 𝑐_𝑇(𝑘) is well behaved when the moments of order 1+𝛼 of the underlying asset exist and are finite. If not all moments of 𝑆(𝑇) exist, this will impose an upper bound on 𝛼. We find that one quarter of this upper bound serves as a good choice for 𝛼.

Using the trapezoid rule for the integral on the right-hand side of (4.1) and setting 𝑣_𝑗=𝜂(𝑗−1), an approximation for 𝐶_𝑇(𝑘) is𝐶_𝑇(𝑘)≈exp(−𝛼𝑘)𝜋^𝑁_𝑗=1𝑒^{−𝑖𝑣𝑗𝑘}𝜓_𝑇𝑣_𝑗𝜂. (4.4) The FFT returns 𝑁 values of 𝑘, and we employ a regular spacing of size ℎ so that our values for 𝑘 are𝑘_𝑢=−𝑏+ℎ(𝑢−1)for𝑢=1,…,𝑁. (4.5) This gives us log-strike levels ranging from −𝑏 to 𝑏, where1𝑏=2𝑁ℎ. (4.6)

In order to apply FFT we define𝜂ℎ=2𝜋𝑁. (4.7)

To obtain an accurate integration with larger values of 𝜂, we incorporate Simpson’s rule weightings into our summation. From (4.1)–(4.7)and Simpson’s rule weightings, we obtain ATM and ITM call value as𝐶𝑘_𝑢=exp−𝛼𝑘_𝑢𝜋^𝑁_𝑗=1𝑒^{−𝑖(2𝜋/𝑁)(𝑗−1)(𝑖−1)}𝑒^{𝑖𝑏𝑣𝑗}𝜓𝑣_𝑗𝜂33+(−1)^𝑗−𝜔_𝑗−1, (4.8) where 𝜔_𝑛 is the Kronecker delta function that is unity for 𝑛=0 and zero otherwise. The summation in (4.8) is an exact application of the FFT.

4.2. Fourier Transform of out-of-the-Money (OTM) Option Prices

In the previous section call values are calculated by an exponential function to obtain square integrable function whose Fourier transform is an analytic function of the characteristic function of the log-price. But, for very short maturities, the call value approaches its non analytic intrinsic value causing the integrand in the Fourier inversion to be high oscillate, and therefore difficult to integrate numerically. We introduce an alternative approach that works with time values only, which is quite similar to the previous approach. But in this case the call price is obtained via the Fourier transform of a modified time value, where the modification involves a hyperbolic sine function instead of an exponential function.

Let 𝑧_𝑇(𝑘) denote the time value of an OTM option, that is, for 𝑘<𝑥_𝑇 we have the put price for 𝑧_𝑇(𝑘) and for 𝑘<𝑥_𝑇 we have the call price. Scaling 𝑆₀=1 for simplicity, 𝑧_𝑇(𝑘) is defined by𝑧_𝑇(𝑘)=𝑒^{−𝑟𝑇∞}_−∞𝑒^𝑘−𝑒^𝑥𝑇𝟏_{𝑥𝑇<𝑘,𝑘<0}+𝑒^𝑥𝑇−𝑒^𝑘𝟏_{𝑥𝑇>𝑘,𝑘>0}𝑞(𝑥)𝑑𝑥, (4.9) where 𝑞(𝑥) is the risk-neutral density of the log-price 𝑥_𝑇. Let 𝜁_𝑇(𝑢) be the Fourier transform of 𝑧_𝑇(𝑘):𝜁_𝑇(𝑢)=^∞_−∞𝑒^𝑖𝑢𝑘𝑧_𝑇(𝑘)𝑑𝑘. (4.10) By considering a damping function sinh(𝛼𝑘), the time value of an option follows a Fourier inversion:𝑧_𝑇(𝑘)=1()1sinh𝛼𝑘𝜋^∞₀𝑒^{−𝑖𝑢𝑘}Υ_𝑇(𝑢)𝑑𝑢, (4.11) where Υ_𝑇(𝑢)=(𝜁_𝑇(𝑢−𝑖𝛼)−𝜁_𝑇(𝑢+𝑖𝛼))/2.

The use of the FFT for calculating OTM option prices is similar to (4.8). The only differences are that they replace the multiplication by exp(−𝛼𝑘_𝑢) with a division by sinh(𝛼𝑘) and the function call to 𝜓(𝑣) is replaced by a function call to Υ_𝑇(𝑢).

5. Simulation Studies

In this section, to compare across the BS, DEJD, and SVDEJD models, we analyze the probability densities of these models. Then, we analyze mainly the impact of 𝜌 and volatility of volatility 𝜎_𝑣 on option pricing under the SVDEJD model. For our FFT methods, we used 𝑁=4096 points in our quadrature, implying a log-strike spacing of ℎ=𝜋/300=0.01047, which is adequate for practice. For the choice of the dampening coefficient in the transform of the modified call price, we used a value of 𝛼=2.55. For the modified time value, we used 𝛼=1.55. Other parameter values used in the computation are listed in Table 1. (We have used analytic moments to set plausible parameter values for the model. For a formal econometric estimator, one could use these moments to develop a generalized method of moments estimator within the framework of Hall and Inoue [21].)

Table 1 

Default parameters for simulation of option prices.

5.1. Probability Densities under Alternative Models

We compare the probability densities of the SVDEJD model, the BS model, and the DEJD model to verify the rationality of our model. Suppose 𝜑(𝑢) is the characteristic function of 𝑥_𝑇 and 𝑞(𝑥) the probability density of our model. From FFT algorithm, 𝑞(𝑥) can be approximated by𝑞(𝑥)≈1𝜋^𝑁_𝑗=1𝑒^{−𝑖(2𝜋/𝑁)(𝑗−1)(𝑘−1)}𝜑(𝑢)𝑘=1,…,𝑁. (5.1) The density 𝑞 has the mean and variance given by𝐸_𝑞(𝑄)=𝜑^′(0)𝑖,Var_𝑞(𝑄)=−𝜑^′′(0)+𝜑^′(0)². (5.2)

Figure 1 shows the figures of the probability density 𝑞(𝑥), compared with the normal density with the same mean and variance given by (5.2). The first figure compares the overall shapes of the densities of the SVDEJD model and the BS model, the second one details the shapes around the peak areas, and the last one shows the right tail. From Figure 1, we can see that the leptokurtic and skewness feature of the density of our model is quite evident. Moreover, additional numerical plots suggest that the feature of skewness becomes more significant if |𝜌| increases, which is impossible for the DEJD model.

Figure 1 

Comparison of the probability densities of the SVDEJD model and the BS model. Except for the maturity time 𝑇=2 years, the parameters used here are shown in Table 1.

We also compare the short-term and long-term densities of the SVDEJD model, the BS model, and the DEJD one. Figure 2 shows their densities under 𝑇=3 months and 𝑇=2 years. From Figure 2, we can see that the SVDEJD model and the DEJD model generate virtually identical densities for short-term options, with a slight departure occurring between the two densities in the upper tail. This means that differential pricing performance between the SVDEJD model and the DEJD model is unlikely to occur when they are applied to price short-term OTM puts and that only when they are applied to deep ITM puts (and deep OTM calls) can differences be observed between these models. Yet, compared to the BS model density, the densities of the two models are distinctly different: they all have leptokurtic and skewness feature. Therefore, the two models can potentially correct the BS model’s tendency to underprice deep OTM puts and overprice deep OTM calls. The long-term density curves in Figure 2 still show significantly different pricing structures between the BS and its two alternatives. But, more importantly, the densities of the SVDEJD model and the DEJD model also exhibit different shapes now. The SVDEJD density has higher peak and assigns more weight to both the entire lower tail and the far upper tail, but less weight to those payoffs than the DEJD.

Figure 2 

Comparison of the short-term and long-term probability densities of the SVDEJD model, the BS model, and the DEJD model. Except for the maturity time 𝑇=3 months and 𝑇=2 years, the parameters used here are shown in Table 1.

Our simulation studies have demonstrated that the SVDEJD model has better performance than the DEJD one on pricing long-term options, while both the DEJD model and the SVDEJD model have better performance than the BS model.

5.2. Effects of the Main Parameter on Option Values

In Table 2, we use the SVDEJD model to examine the effects of volatility of volatility 𝜎_𝑣, the correlation coefficient 𝜌, exercise price 𝐾, and maturity time 𝑇 on option values. We analyze the prices of three-month call options and two-year call options. To examine the effect of the negative correlation coefficient, we have calculated the model with 𝜌=−0.8, 𝜌=0, and 𝜌=0.8. The prices for three-month call options associated with volatility of volatility 𝜎_𝑣=0.15 and 𝜎_𝑣=0.25 are relatively close. With 𝜌=−0.8, the largest price difference is 0.0149; with 𝜌=0, the largest price difference is only 0.0008; with 𝜌=0.8, the largest price difference is 0.0213. However, the difference is significantly larger when longer-time horizons such as two-year call options are valued. With 𝜌=−0.8, the largest price difference is an increase on 0.0149 to 0.0545 and the effect of volatility of volatility is an increase for ITM calls and a decrease for OTM calls. With 𝜌=0, the largest price difference is increase of 0.0008 to 0.01 and the effect of volatility of volatility is a small decreas that is negligible for option values. With 𝜌=0.8, the largest price difference is an increase of 0.0213 to 0.0633 and the effect of volatility of volatility is a decrease for ITM and ATM calls and an increase for OTM calls. The correlation parameter 𝜌 has several effects depending on the relation between the strike price and the current stock price. A negative 𝜌 tends to produce higher values for ITM calls and lower values for OTM money calls.

Table 2 

The effects of volatility of volatility, exercise price 𝐾, maturity time 𝑇, and correlation 𝜌 on option values.

We have also compared the model with the BS model, which can be interpreted as a first-order approximation with no jumps, and 𝜌=0. A common practice is to set the implied volatility in the BS model so that the model matches the price for the option with a strike price closest to the current stock price. For some comparisons not reported here, we have set the implied volatility in the BS model so that it matches the price generated by the stochastic volatility model for an ATM option. The BS implied volatility is very close to the expected volatility under the risk-neutral distribution when short-term options are valued. When longer-term options are used, there is a significant difference between the BS implied volatility and the expected volatility. As an approximation, the BS model tends to undervalue ITM calls and overvalue OTM calls.

6. Conclusion

The SVDEJD model incorporates several important features of stock returns. We derive a closed-form solution for European call options in the model by using the martingale method, Fourier inversion transform formula, and Feynman-Kac theorem. Using FFT, we obtain fast and accurate numerical solution to European option under the model. The comparison of densities of the alternative models shows that the SVDEJD model has better pricing performance on long-time options. An analysis of the model reveals that volatility of volatility 𝜎_𝑣 and the correlation coefficient 𝜌 have significant impact on option values, especially long-time option, stock returns are negatively correlated with volatility, and these negative correlations are important for option valuation.

Appendix

A.1. Matlab Codes for ITM and ATM Options Pricing by FFT
function CV =inSVDexpJ(ata1, ata2, lamta, sigma, thetav, alphav,rho, sigmav, r, p, s0, v0, strike, T)x0 = log(s0)alpha = 2.55N = 4096c = 600eta = c/Nb =pi/etau = [0:N-1]*etalamda = 2*b/Nposition = (log(strike) + b)/lamda + 1v = u - (alpha+1)*ik=p*ata1/(ata1-1)+(1-p)*ata2/(ata2+1)-1l=p*ata1./(ata1-i*v)+(1-p)*ata2./(ata2+i*v)m=sqrt((alphav-i*v*rho*sigma*sigmav).^∧2+i*v.*(1-i*v)*(sigma * sigmav)^∧2)n=2*m+(alphav-m-i*v*rho*sigma*sigmav).*(1-exp(-m*T))A=(2*m./n).^∧(2*thetav/sigmav^∧2)B=i*v*x0+(thetav*(alphav-m)*T)/sigmav^∧2-(i*v*rho*sigma*thetav*T)/sigmav…+lamta*T*(l-i*v*k-1)+i*v*r*TC=(i*v.*(i*v-1)*sigma^∧2.*(1-exp(-m*T)))./ncharFunc=A.*exp(C*v0+B)ModifiedCharFunc = charFunc*exp(-r*T)./(alpha^∧2+alpha - u.^∧2 + i*(2*alpha +1)*u)SimpsonW = 1/3*(3 + (−1).^∧[1:N] - [1, zeros(1,N-1)])FftFunc = exp(i*b*u).*ModifiedCharFunc*eta.*SimpsonWpayoff = real(fft(FftFunc))CallValueM = (exp(-log(strike)*alpha))’*payoff/piformat shortCV= CallValueM(round(position)).

A.2. Matlab Codes for OTM Options Pricing by FFT
function CV=outSVDexpJ(ata1,ata2,lamta,sigma,thetav,alphav, rho,sigmav,r,p,s0,v0,strike,T)x0 = log(s0)alpha = 1.55N=4096c = 600eta = c/Nb =pi/etau = [0:N-1]*etalamda = 2*b/Nposition = (log(strike) + b)/lamda + 1w1 = u-i*alphaw2 = u+i*alphav1 = u-i*alpha -iv2 = u+i*alpha -ik=p*ata1/(ata1-1)+(1-p)*ata2/(ata2+1)-1l1=p*ata1./(ata1-i*v1)+(1-p)*ata2./(ata2+i*v1)m1=sqrt((alphav-i*v1*rho*sigma*sigmav).^∧2+i*v1.*(1-i*v1)*(sigma*sigmav)^∧2)n1=2*m1+(alphav-m1-i*v1*rho*sigma*sigmav).*(1-exp(-m1*T))A1=(2*m1./n1).^∧(2*thetav/sigmav^∧2)B1=i*v1*x0+(thetav*(alphav-m1)*T)/sigmav^∧2-(i*v1*rho*sigma*thetav*T)/sigmav…+lamta*T*(l1-i*v1*k-1)+i*v1*r*TC1=(i*v1.*(i*v1-1)*sigma^∧2.*(1-exp(-m1*T)))./n1charFunc1=A1.*exp(C1*v0+B1)ModifiedCharFunc1 = exp(-r*T)*(1./(1+i*w1)…-exp(r*T)./(i*w1) - charFunc1./(w1.^∧2 - i*w1))l2=p*ata1./(ata1-i*v2)+(1-p)*ata2./(ata2+i*v2)m2=sqrt((alphav-i*v2*rho*sigma*sigmav).^∧2+i*v2.*(1-i*v2)*(sigma*sigmav)^∧2)n2=2*m2+(alphav-m2-i*v2*rho*sigma*sigmav).*(1-exp(-m2*T))A2=(2*m2./n2).^∧(2*thetav/sigmav^∧2)B2=i*v2*x0+(thetav*(alphav-m2)*T)/sigmav^∧2-(i*v2*rho*sigma*thetav*T)/sigmav…+lamta*T*(l2-i*v2*k-1)+i*v2*r*TC2=(i*v2.*(i*v2-1)*sigma^∧2.*(1-exp(-m2*T)))./n2charFunc2=A2.*exp(C2*v0+B2)ModifiedCharFunc2 = exp(-r*T)*(1./(1+i*w2)- exp(r*T)./(i*w2)…- charFunc2./(w2.^∧2 - i*w2))ModifiedCharFuncCombo = (ModifiedCharFunc1 - ModifiedCharFunc2)/2SimpsonW = 1/3*(3 + (−1).^∧[1:N] - [1, zeros(1,N-1)])FftFunc = exp(i*b*u).*ModifiedCharFuncCombo*eta.*SimpsonWpayoff = real(fft(FftFunc))CallValueM = payoff/pi/sinh(alpha*log(strike))format shortCV= CallValueM(round(position)).

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 11171266) and the Science Plan Foundation of the Education Bureau of Shanxi Province (nos. 11JK0491 and 11JK0493). The authors are grateful to the reviewers as well as the editor for their valuable comments and suggestions, which led to a substantial improvement of the paper.