An efficient pricing algorithm for American options with double stochastic volatilities and double jumps

Abstract

The purpose of the paper is to provide an efficient pricing algorithm for American options with stochastic volatilities and jumps. This paper extends the double Heston model with double exponential jumps and derives the characteristic function of the model by Feynman–Kac theorem. With the obtained characteristic function, this paper also extends the Fourier-cosine expansion method for pricing Bermudan options to the model. Based on the COS method, this paper approximates American options by using Richardson extrapolation schemes on a series of Bermudan options and provides a pricing algorithm for American put options. Numerical results show that the proposed pricing algorithm is efficient, especially for short-term American put options.

Introduction

It is important for efficiently pricing American options, because a majority of derivative contracts are American style. Single-factor stochastic volatility models, including Stein and Stein, ¹ Heston, ² and Schöbel and Zhu, ³ can explain the volatility smile observed in the real market. Many literatures consider American options pricing under these models by developing some numerical methods, including the finite difference methods of Ikonen and Toivanen,^4,5 Ito and Toivanen, ⁶ Zhu and Chen, ⁷ and some extensions, such as Kunoth et al., ⁸ Rambeerich et al., ⁹ Ballestra and Pacelli, ¹⁰ and Burkovska et al., ¹¹ the Monte Carlo simulation method of Abbas-Turki and Lapeyre ¹² and the tree methods of Beliaeva and Nawalkha ¹³ and Ruckdeschel et al. ¹⁴

However, single-factor stochastic volatility models are not able to fit the implied volatility smile very well. Evidence from Cont and Tankov, ¹⁵ Fonseca et al., ¹⁶ Christoffersen et al., ¹⁷ and Fouque and Lorig ¹⁸ indicate that single-factor models can do a poor job in capturing the term structures of implicit volatilities over time. By introducing another volatility process, Christoffersen et al. ¹⁷ propose the double Heston model which can provide better empirical fit to the market price than the Heston model. However, Yu and Zhang, ¹⁹ Zhou and Zhu, ²⁰ González-Urteaga, ²¹ and Jang et al. ²² provide strong evidence for stochastic volatility and jumps in prices. By adding the log-normal jump to the asset price, the stochastic volatility jump-diffusion models including the Bates model ²³ are proposed and applied in option pricing. Compare to log-normal jump, the double exponential jump-diffusion model proposed by Kou ²⁴ leads to tractable pricing formulas for path dependent options. Motivated by the superior features of the double Heston model and the double exponential jumps, this paper proposes a new model by combining double Heston model and double exponential jumps.

Additional two factors lead to a high-dimensional partial integro-differential equation which makes the aforementioned finite difference methods quite complex and difficult to be extended to double stochastic volatilities and jumps case. In contrast, the tree method and the Monte Carlo method are easier to be extended to this case. However, exponentially rising numbers of nodes and a large amount of simulation times will make these methods time-consuming. The COS method proposed by Fang and Oosterlee ²⁵ is a very efficient method for pricing European option and has been extended to Bermudan options.^26,27 Based on the COS method, this paper takes American option as the combination of several Bermudan options and approximates the option price using a Richardson extrapolation technique. ²⁸

The rest of the paper is organized as follows. The following section develops the underlying pricing model. The subsequent section derives the characteristic function of the model. Then the method for pricing American options is detailed, followed by some numerical experiments. The last section concludes.

The model

Assume that $W_1(t)$ , $W_2(t)$ , $Z_1(t)$ and $Z_2(t)$ are all standard Brownian motions which satisfy $Cov(d{W}_1(t),d{Z}_1(t))={\rho }_{1}dt,$ $Cov(d{W}_2(t),d{Z}_2(t))={\rho }_{2}dt,\text{ and }Cov(d{W}_1(t),d{W}_2(t))=0$ and $Cov(d{Z}_1(t),d{Z}_2(t))=0.$ Suppose that the asset price process $S(t)$ is governed by the following double Heston model with double exponential jumps (DHestonDJ)

$$\{\begin{array}{l}\frac{dS(t)}{S(t)}=(r-\lambda \delta )dt+{\displaystyle \sum _{j=1}^2\sqrt{V_j(t)}d{W}_j(t)}\text{+}d\left({\displaystyle \sum _{j=1}^{N(t)}(J_j-1)}\right),\\ d{V}_1(t)=k_1({\theta }_1-V_1(t))dt+{\sigma }_1\sqrt{V_1(t)}d{Z}_1(t),\\ d{V}_2(t)=k_2({\theta }_2-V_2(t))dt+{\sigma }_2\sqrt{V_2(t)}d{Z}_2(t)\end{array}$$

(1)

where $r$ is the risk-neutral interest rate, $k_j,{\theta }_j,{\sigma }_j(j=1,2)$ is the rates of reversion, long-run mean volatility, and instantaneous volatility of variance process $V_j(t)$ , respectively. $V_j(t)$ remains strictly positive if the Feller conditions $2k_j{\theta }_j\ge {\sigma }_j{\hspace{0.17em}}^2$ is satisfied. $N(t)$ is a Poisson process with constant intensity $\lambda >0$ . $J={(J_j)}_{j\ge 1}$ is a sequence of independent identically distribute nonnegative random variables, such that $Y=\log J$ has an asymmetric double exponential distribution with the density.

$$f_Y(y)=p{\eta }_1{e}^{-{\eta }_{1}y}{\text{1}}_{(y\ge 0)}+q{\eta }_2{e}^{{\eta }_{2}y}{\text{1}}_{(y<0)},{\eta }_1>1,{\eta }_2>0$$

(2)

where $\text{1}$ denotes the indicator function, so ${\text{1}}_{(y\ge 0)}$ equals 1 if $y\ge 0$ , but 0 otherwise. $p,q\ge 0,p+q=1$ are the probability of the up-move jump and down-move jump, respectively. Suppose the processes $W_1(t)$ , $W_2(t)$ , $Z_1(t)$ , $Z_2(t)$ are all independent of $N(t)$ and $J$ , $\delta =E(J-1)$ . The density (2) implies $\delta =\frac{p{\eta }_1}{{\eta }_1-1}+\frac{q{\eta }_2}{{\eta }_2+1}-1$ . Suppose $V_1(0)\text{=}$ $V_1,V_2(0)=V_2,S(0)=S$ _.

Remark The model contains the following known models as special cases.

the Heston model by setting $\lambda =0$ and $k_2=0,{\theta }_2=0,$ ${\sigma }_2=0$ ;

the double Heston model by $\lambda =0$ ;

the double exponential jump-diffusion model by setting $k_1=0,{\theta }_1=0,{\sigma }_1=0,k_2=0,{\theta }_2=0,{\sigma }_2=0$ .

Deriving the characteristic function

Let $K$ denote the strike price. For the asset price $S(t)$ satisfying the DHestonDJ model (1), we define the characteristic function of the asset price $X_t=\log S(\tau )/K$ at time $\tau =T-t$ as follows

$${\varphi }_X(u,x,V_1,V_2,\tau )=E(e^{iu{X}_t})$$

(3)

where $i$ is imaginary unit. By Feynman–Kac theorem, ${\varphi }_X(\cdot )$ satisfies the following partial integro-differential equation

$$\{\begin{array}{l}-{\varphi }_{\tau }+\left[r-\lambda \delta -\frac{1}{2}(V_1+V_2)\right]{\varphi }_x+\frac{1}{2}(V_1+V_2){\varphi }_{xx}\hfill \\ +\frac{1}{2}{\sigma }_1^2{V}_1{\varphi }_{V_1{V}_1}+\frac{1}{2}{\sigma }_2^2{V}_2{\varphi }_{V_2{V}_2}+{\rho }_1{\sigma }_1{V}_1{\varphi }_{x{V}_1}\hfill \\ +{\rho }_2{\sigma }_2{V}_2{\varphi }_{x{V}_2}+k_1({\theta }_1-V_1){\varphi }_{V_1}\hfill \\ +k_2({\theta }_2-V_2){\varphi }_{V_2}+\lambda {\displaystyle {\int }_{-\infty }^{+\infty }(\varphi (x+J)-\varphi (x))\varpi (J)\mathrm{d}J=0,}\hfill \\ \varphi (u,x,V_1,V_2,0)=e^{ux}\hfill \end{array}$$

(4)

By direct calculation, we rewrite the integral term in equation (4) as follows

$$\begin{array}{l}\lambda {\displaystyle {\int }_{-\infty }^{+\infty }(\varphi (x+J)-\varphi (x))\varpi (J)\mathrm{d}J}\hfill \\ =\lambda {\displaystyle {\int }_{-\infty }^{+\infty }(E[e^{u(x+J)}]-E[e^{ux}])\varpi (J)\mathrm{d}J}\\ =\lambda \varphi {\displaystyle {\int }_{-\infty }^{+\infty }(e^{uJ}-1)\varpi (J)\mathrm{d}J}\\ =\lambda \varphi \left(\frac{p}{{\eta }_1-u}+\frac{q}{{\eta }_2+u}-1\right)\end{array}$$

(5)

According to Duffie et al. ²⁹ and Heston, ² ${\varphi }_X(\cdot )$ has the following form

$${\varphi }_X(u,x,V_1,V_2,\tau )=e^{[ux+A(u,\tau )+B(u,\tau )V_1+C(u,\tau )V_2]}$$

(6)

with initial conditions $A(u,0)=0,$ $B(u,0)=0$ , $C(u,0)=0$ .

Substituting equations (5) and (6) into equation (4) yields

$$\begin{array}{l}0=-(A(u,\tau )+B(u,\tau )V_1+C(u,\tau )V_2)+\frac{1}{2}(V_1+V_2)u^2\\ +\frac{1}{2}{\sigma }_1^2{V}_1{B}^2(u,\tau )+\frac{1}{2}{\sigma }_2^2{V}_2{C}^2(u,\tau )+{\rho }_1{\sigma }_1{V}_{1}uB(u,\tau )\\ +{\rho }_2{\sigma }_2{V}_{2}uC(u,\tau )+k_1({\theta }_1-V_1)B(u,\tau )\\ +k_2({\theta }_2-V_2)C(u,\tau )+\lambda \left(\frac{p}{{\eta }_1-u}+\frac{q}{{\eta }_2+u}-1\right)\\ +\left[r-\lambda \delta -\frac{1}{2}(V_1+V_2)\right]u\end{array}$$

The above equation produces the following system of three ordinary differential equations

$$\{\begin{array}{l}\frac{\partial B(u,\tau )}{\partial \tau }=\frac{1}{2}{\sigma }_1^2{B}^2(u,\tau )+(u{\rho }_1{\sigma }_1-k_1)B(u,\tau )-\frac{1}{2}u(i+u),\\ \frac{\partial C(u,\tau )}{\partial \tau }=\frac{1}{2}{\sigma }_2^2{C}^2(u,\tau )+(u{\rho }_2{\sigma }_2-k_2)C(u,\tau )-\frac{1}{2}u(i+u),\\ \frac{\partial A(u,\tau )}{\partial \tau }=ur+B(u,\tau )k_1{\theta }_1+C(u,\tau )k_2{\theta }_2+\lambda \Lambda (u,J)\end{array}$$

where $\Lambda (u,J)=\frac{p}{{\eta }_1-u}+\frac{q}{{\eta }_2+u}-1-u\left(\frac{p}{{\eta }_1-1}+\frac{q}{{\eta }_2+1}-1\right)$ . By solving the above equations, we obtain the characteristic function ${\varphi }_X(u,x,V_1,V_2,\tau )$ as below

$${\varphi }_X(u,x,V_1,V_2,\tau )=e^{\left[ux+A(u,\tau )+{\displaystyle \sum _{j=1}^2{B}_j(u,\tau )V_j}\right]}$$

(7)

where

$$\begin{array}{l}A(u,\tau )=-{\displaystyle \sum _{j=1}^2\frac{k_j{\theta }_j}{{\sigma }_j^2}\left[(D_j+{\beta }_j)\tau +2\ln \left(\frac{1-G_j{e}^{-D_j\tau }}{1-G_j}\right)\right]}\\ +ur\tau +\lambda \Lambda (u,J),\\ B_j(u,\tau )=-\frac{D_j+{\beta }_j}{{\sigma }_j^2}\frac{1-e^{-D_j\tau }}{1-G_j{e}^{-D_j\tau }},\\ D_j=\sqrt{{\beta }_j^2-2{\sigma }_j^2{\gamma }_j},G_j=\frac{{\beta }_j+D_j}{{\beta }_j-D_j},\\ {\beta }_j={\sigma }_j{\rho }_{j}u-k_j,{\gamma }_j=-\frac{1}{2}u(i+u)\end{array}$$

For the forthcoming computation, we also provide the first two cumulants $c_1=\frac{1}{i}\frac{\partial \ln \varphi }{\partial u}|{\hspace{0.17em}}_{u=0}$ , $c_2=\frac{1}{i^2}\frac{{\partial }^2\ln \varphi }{\partial u^2}|{\hspace{0.17em}}_{u=0}$ as follows

$$\begin{array}{l}{c}_1=\log \left(\frac{S}{K}\right)+{\displaystyle \sum _{j=1}^2\frac{{\theta }_j{V}_j}{2k_j}\left(1-e^{-k_j\tau }\right)}\hfill \\ +\left(r-\frac{{\theta }_1-{\theta }_2}{2}+\lambda \left(\frac{p}{{\eta }_1}+\frac{q}{{\eta }_2}-\left(\frac{p}{{\eta }_1-1}+\frac{q}{{\eta }_2+1}-1\right)\right)\right)\tau \hfill \end{array}$$

(8)

$$c_2=\lambda \left(\frac{2p}{{\eta }_1{\hspace{0.17em}}^2}+\frac{2q}{{\eta }_2{\hspace{0.17em}}^2}\right)+{\displaystyle \sum _{j=1}^2\left({\Pi }_{0j}+{\Pi }_{1j}{e}^{-k_j\tau }+{\Pi }_{2j}{e}^{-2k_j\tau }\right)}$$

(9)

where

$$\begin{array}{l}{\Pi }_{0j}={\theta }_j\tau +\frac{1}{k_j}({\theta }_j-V_j+{\theta }_j{\sigma }_j{\rho }_j\tau )\hfill \\ +\frac{1}{k_j{\hspace{0.17em}}^2}\left[{\sigma }_j{\rho }_j(2{\theta }_j-V_j)+\frac{1}{4}{\theta }_j{\sigma }_j{\hspace{0.17em}}^2\tau \right]\hfill \\ -\frac{1}{k_j{\hspace{0.17em}}^3}\frac{{\sigma }_j{\hspace{0.17em}}^2}{8(5{\theta }_j-2V_j)},\hfill \\ {\Pi }_{1j}=\frac{1}{k_1{\hspace{0.17em}}^3{k}_2{\hspace{0.17em}}^3}{\theta }_j{\rho }_j(V_j-2{\theta }_j)+\frac{1}{k_j}(1-{\sigma }_j{\rho }_j\tau )({\theta }_j-V_j)\hfill \\ +\frac{{\theta }_j{\sigma }_j{\hspace{0.17em}}^2}{2k_j{\hspace{0.17em}}^3}+\frac{\tau {\sigma }_j{\hspace{0.17em}}^2({\theta }_j-V_j)}{2k_j{\hspace{0.17em}}^2},\\ {\Pi }_{2j}=\frac{1}{k_j{\hspace{0.17em}}^3}\frac{{\sigma }_j{\hspace{0.17em}}^2}{8({\theta }_j-2V_j)}\hfill \end{array}$$

Pricing method and algorithm for American options

Let $t_0$ be the current time, $T$ be the maturity, and $I=\{t_1,t_2,\dots ,t_M\}$ be the collection of all exercise dates with ${\Delta }_t=t_m-t_{m-1}(m=2,3,\dots ,M)$ . Let $A(t_0,S(t_0))$ denote the price of an American option and $V^{\Delta }(t_0,S(t_0))$ denote the price of a Bermudan option with the exercise times equally spaced, where $\Delta =t_m-t_{m-1}$ and $\Delta \in I$ then

$$\underset{\Delta \to 0}{\lim }{V}^{\Delta }(t_0,S(t_0))=A(t_0,S(t_0))$$

(10)

Assume that $V^{\Delta }(t_0,S(t_0))$ can be expanded to $A(t_0,S(t_0))$ with respect to $\Delta$ as

$$\begin{array}{l}{V}^{\Delta }(t_0,S(t_0))=A(t_0,S(t_0))+a_1\Delta +a_2{\Delta }^2+\dots +a_k{\Delta }^{q_k}\hfill \\ +{\rm O}({\Delta }^{q_{k+1}})\end{array}$$

(11)

with parameters $a_1,a_2,\dots ,a_k.$ To approximate $A(t_0,S(t_0))$ , we compute $V^{\Delta }(t_0,S(t_0))$ a number of times with successively smaller steps, ${\Delta }_1=\Delta ,{\Delta }_2=\Delta /2,{\Delta }_3=\Delta /2^2,$ $\dots$ . In such a way, we obtain the increasing exercise opportunities $N_{{\Delta }_1},N_{{\Delta }_2},N_{{\Delta }_3},\dots$ and a sequence of approximation $V^{{\Delta }_1},V^{{\Delta }_2},V^{{\Delta }_3},\dots$ . Based on polynomial interpolation and an asymptotic $\Delta$ -expansion, we can construct repeated Richardson extrapolation scheme ²⁸ as follows

$$\begin{array}{l}{A}_{j,0}=V^{{\Delta }_j},\\ A_{j,m}=A_{j+1,m-1}+\frac{A_{j+1,m-1}-A_{j,m-1}}{\frac{{\Delta }_j}{{\Delta }_{j+m}}-1}\end{array}$$

(12)

where $j=1,2,\dots$ and $m=1,2,\dots ,j-1$ . Figure 1 indicates the four-point Richardson extrapolation scheme.

Figure 1.

The four-point Richardson extrapolation scheme.

Bermudan options pricing based on the COS method

Assume that $x=\log (S(t)/K)$ , $C(t,x)$ is the continuation value and $g(t,x)$ is the value of the payoff. The value of a Bermudan option with $M$ exercise dates can be expressed by a backward recursion as

$$\{\begin{array}{l}V(t_M,x_M)=g(t_M,x_M),\\ C(t_m,x_m)=e^{-r\Delta t}{\displaystyle {\int }_{R}V(t_{m\text{+}1},y)f(y|x_m)\mathrm{d}y,}\\ V(t_m,x_m)=\max (g(t_m,x),C(t_m,x_m)),\\ V(t_0,x_0)=C(t_0,x_0)\end{array}$$

(13)

By Fourier-cosine series expansion, $f(y|x_m)$ can be approximated by

$$f(y|x_m)\approx \frac{2}{b-a}{\displaystyle \sum _{k=0}^{\infty }\mathrm{Re}}\left\{\varphi \left(\frac{k\pi }{b-a};x_m\right)e^{-ik\pi \frac{a}{b-a}}\right\}\cos \left(k\pi \frac{y-a}{b-a}\right)$$

where $\mathrm{Re}(\cdot )$ denotes taking the real part of the argument. Putting the above equation into equation (13) and interchanging integration and summation gives the COS formula for approximating $C(t_m,x_m)$ by $\widehat{c}(t_m,x_m)$ as

$$\widehat{c}(t_m,x_m)=e^{-r\Delta t}{\displaystyle \sum _{k=0}^{N-1}\mathrm{Re}}\left\{\varphi \left(\frac{k\pi }{b-a};x_m\right)e^{-ik\pi \frac{a}{b-a}}\right\}{V}_k(t_{m\text{+}1})$$

(14)

where

$$V_k(t_{m\text{+}1})\text{=}\frac{2}{b-a}{\displaystyle {\int }_a^b{V}_k(t_{m+1},y)\cos \left(k\pi \frac{y-a}{b-a}\right)}\mathrm{d}y$$

(15)

To obtain Bermudan option price, we need to recover $V_K\left(t_{m-1}\right)$ from $V_K\left(t_m\right)$ . We split $V_K\left(t_m\right)$ into two parts by the early exercise point $x_m^{*}$

$$V_k(t_m)\text{=}\{\begin{array}{l}{C}_k(a,x_m^{*},t_m)+G_k(x_m^{*},b),\text{Call}\\ G_k(a,x_m^{*},t_m)+C_k(x_m^{*},b),\text{Put}\end{array}$$

(16)

where

$$G_k(x_1,x_2)\text{=}\frac{2}{b-a}{\displaystyle {\int }_{x_1}^{x_2}g(t_m,y)\cos \left(k\pi \frac{y-a}{b-a}\right)}\mathrm{d}y$$

(17)

$$C_k(x_1,x_2)\text{=}\frac{2}{b-a}{\displaystyle {\int }_{x_1}^{x_2}\widehat{c}(t_m,y)\cos \left(k\pi \frac{y-a}{b-a}\right)}\mathrm{d}y$$

(18)

$G_k(x_1,x_2)$ has the following analytic expression

$$G_k(x_1,x_2)\text{=}\{\begin{array}{l}\frac{2}{b-a}K[{\chi }_k(x_1,x_2)-{\psi }_k(x_1,x_2)],\text{Call}\\ \frac{2}{b-a}K[{\psi }_k(x_1,x_2)-{\chi }_k(x_1,x_2)],\text{Put}\end{array}$$

(19)

where

$$\begin{array}{l}{\chi }_k(x_1,x_2)=\frac{1}{1+{\left(\frac{k\pi }{b-a}\right)}^2}[\cos \left(\frac{x_2-a}{b-a}\right)e^{x_2}-\cos \left(\frac{x_1-a}{b-a}\right)e^{x_1}\hfill \\ +\frac{k\pi }{b-a}\left(\sin \left(\frac{x_2-a}{b-a}\right)e^{x_2}-\sin \left(\frac{x_1-a}{b-a}\right)e^{x_1}\right)],\hfill \\ {\psi }_k(x_1,x_2)=\{\begin{array}{ll}\frac{b-a}{k\pi }\left[\sin \left(\frac{x_2-a}{b-a}\right)-\sin \left(\frac{x_1-a}{b-a}\right)\right],\hfill & k\ne 0,\hfill \\ x_2-x_1,\hfill & k=0.\hfill \end{array}\end{array}$$

Putting equation (14) into $C_k(x_1,x_2)$ , we have

$$C(t_{m-1},x_{m-1})\approx e^{-r\Delta t}\mathrm{Re}\left[{\displaystyle \sum _{j=0}^{N-1}\varphi \left(\frac{j\pi }{b-a}\right)V_j(t_m)M_{k,j}(x_1,x_2)}\right]$$

(20)

where

$$\begin{array}{l}{V}_j\left(t_m\right)=G_j\left(a,0\right),\hfill \\ M_{k,j}(x_1,x_2)\text{=}\frac{2K}{b-a}{\displaystyle {\int }_{x_1}^{x_2}{e}^{ij\pi \frac{y-a}{b-a}}\cos \left(k\pi \frac{y-a}{b-a}\right)}\text{d}y\end{array}$$

By basic calculation, $M_{j,k}\left(x_1,x_2\right)$ can be formulated as

$$M_{k,j}\left(x_1,x_2\right)=-\frac{i}{k}\left[M_{k,j}^c\left(x_1,x_2\right)+M_{k,j}^c\left(x_1,x_2\right)\right]$$

(21)

where

$$M_{k,j}^c\left(x_1,x_2\right)=\{\begin{array}{cc}\frac{\left(x_2-x_1\right)\pi i}{b-a},& k=j=0,\\ \frac{e^{{\hspace{0.17em}}^{i\left(j+k\right)\frac{\left(x_2-a\right)\pi }{b-a}}}-e^{i\left(j+k\right)\frac{\left(x_{{\hspace{0.17em}}_1}-a\right)\pi }{b-a}}}{j+k},& \mathrm{others}\end{array}$$

(22)

$$M_{k,j}^s\left(x_1,x_2\right)=\{\begin{array}{cc}\frac{\left(x_2-x_1\right)\pi i}{b-a},& k=j,\\ \frac{e^{{\hspace{0.17em}}^{i\left(j-k\right)\frac{\left(x_2-a\right)\pi }{b-a}}}-e^{i\left(j-k\right)\frac{\left(x_{{\hspace{0.17em}}_1}-a\right)\pi }{b-a}}}{j+k},& \mathrm{others}\end{array}$$

(23)

Set

$$u={\left[\varphi \left(\frac{j\pi }{b-a}\right)v_j\left(t_m\right)\right]}_{J=0}^{N-1}$$

(24)

Equation (20) can be written in the following matrix-vector-product form

$$C_k\left(x_1,x_2\right)=\frac{e^{-r\Delta t}}{\pi }\mathbf{I}m\left[\left({\mathbf{M}}_{\mathbf{c}}+{\mathbf{M}}_{\mathbf{s}}\right)\mathbf{u}\right]$$

(25)

where $\mathrm{Im}(\cdot )$ denotes taking the imaginary part of the argument

$${\mathbf{M}}_{\mathbf{c}}=\left(\begin{array}{lllll}{m}_0\hfill & m_1\hfill & m_2\hfill & \cdots \hfill & m_{N-1}\hfill \\ m_1\hfill & m_2\hfill & \cdots \hfill & \cdots \hfill & m_N\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & \ddots \hfill & ⋮\hfill \\ m_{N-2}\hfill & m_{N-1}\hfill & \cdots \hfill & \cdots \hfill & m_{2N-3}\hfill \\ m_{N-1}\hfill & \cdots \hfill & \cdots \hfill & m_{2N-1}\hfill & m_{2N-2}\hfill \end{array}\right)$$

(26)

$${\mathbf{M}}_{\mathrm{s}}=\left(\begin{array}{lllll}{m}_0\hfill & m_1\hfill & m_2\hfill & \cdots \hfill & m_{N-1}\hfill \\ m_{-1}\hfill & m_0\hfill & \cdots \hfill & \cdots \hfill & m_{N-2}\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & \ddots \hfill & ⋮\hfill \\ m_{2-N}\hfill & \cdots \hfill & m_{-1}\hfill & m_0\hfill & m_1\hfill \\ m_{1-N}\hfill & m_{2-N}\hfill & \cdots \hfill & m_{-1}\hfill & m_0\hfill \end{array}\right)$$

(27)

$$m_j=\{\begin{array}{cc}\frac{\left(x_2-x_1\right)\pi i}{b-a},& j=0,\\ \frac{e^{ij\frac{(x_2-a)\pi }{b-a}}-e^{ij\frac{(x_1-a)\pi }{b-a}}}{j},& j\ne 0.\end{array}$$

(28)

It is obvious that ${\mathbf{M}}_{\text{c}}$ is a Hankel matrix and ${\mathbf{M}}_s$ is a Toeplitz matrix. With the help of fast Fourier transform algorithm, ${\mathbf{M}}_{\mathrm{c}}\mathbf{u}$ and ${\mathbf{M}}_s\mathbf{u}$ can be computed efficiently.

The accuracy of the COS method is related to the size of the integrating range $\left[a,b\right]$ . According to Fang and Oosterlee, ²⁵ we define the truncation range as follows

$$[a,b]=\left[c_1\text{-}L\sqrt{\left|c_2\right|},c_1+L\sqrt{\left|c_2\right|}\right]$$

(29)

where $c_j=\frac{1}{i^j}\frac{{\partial }^j\ln \varphi }{\partial u^j}|{\hspace{0.17em}}_{u=0}(j=1,2)$ , L is a proportion constant. For the DHestonDJ model (1), $c_1$ and $c_2$ is provided by equations (8) and (9), respectively.

The hybrid pricing algorithm based on COS and Richardson extrapolation technology

By combining the COS and Richardson extrapolation technology, we provide the following hybrid algorithm for pricing American put options.

Algorithm 1: The algorithm for pricing an American put option

Numerical experiments

We use the four-point Richardson extrapolation schemes and the COS method (Re-COS) to price American put options. The exercise opportunities of each Bermuda option are ${\text{2}}^9,2^8,2^7,2^6$ , respectively. For comparison, we also use four-point Richardson extrapolation schemes and the convolution method proposed by Lord et al. ³⁰ (Re-CONV) to evaluate American options. For COS method, we use grid points $N=2^{10}$ in equation (14) and $L=13$ in equation (29). For CONV method, we use dampening factor 0.5, grid points $N=2^{10}$ and $L=3$ . The model parameter values used in the computation are: $S=100,{\eta }_1={\eta }_2=40,p=0.6,k_1=12,{\theta }_1=0.05,{\sigma }_1=0.9,$ ${\rho }_1=-0.5,V_1=0.05,k_2=16,{\theta }_2=0.03,{\sigma }_2=0.9,{\rho }_2=-0.5,$ $V_2=0.02,r=\mathrm{0.03.}$ We specify three maturities $T=1/4$ , $T=1/2$ , and $T=1$ . Under each maturity, we specify $\lambda =8$ , $\lambda =9$ , and $\lambda =10$ . Table 1 reports the results.

Table 1.

Comparison of the accuracy between the Re-COS method and the Re-CONV method for pricing American put options with different strike, maturity and $\lambda$ under the DHestonDJ model.

Strike	$\lambda =8$		$\lambda =9$		$\lambda =10$
	Re-COS	Re-CONV	Re-COS	Re-CONV	Re-COS	Re-CONV
$T=1/4$
80	0.2717	0.2718	0.2819	0.2819	0.2922	0.2921
90	1.6115	1.6108	1.6407	1.6398	1.6698	1.6689
100	5.2923	5.2908	5.3346	5.3330	5.3767	5.3750
110	11.7460	11.7444	11.7820	11.7804	11.8180	11.8163
120	20.2933	20.2921	20.3110	20.3098	20.3289	20.3277
$T=1/2$
80	1.0166	1.0166	1.0433	1.0433	1.0701	1.0701
90	3.1906	3.1903	3.2386	3.2382	3.2863	3.2860
100	7.2976	7.2968	7.3578	7.3570	7.4176	7.4169
110	13.4182	13.4171	13.4754	13.4743	13.5323	13.5312
120	21.2164	21.2152	21.2582	21.2570	21.3000	21.2989
$T=1$
80	2.5532	2.5542	2.6046	2.6057	2.6560	2.6571
90	5.4997	5.5008	5.5722	5.5733	5.6444	5.6454
100	9.9384	9.9391	10.0229	10.0236	10.1069	10.1076
110	15.8335	15.8338	15.9183	15.9185	16.0025	16.0028
120	23.0097	23.0095	23.0839	23.0837	23.1578	23.1576

DHestonDJ: double Heston model with double exponential jumps; Re-CONV: Richardson extrapolation schemes and the CONV method; Re-COS: Richardson extrapolation schemes and the COS method.

Our numerical experiments show that the two methods have almost the same accuracy. If we take the Re-CONV method as the benchmark, the relative error of Re-COS does not exceed 0.0549%, while the absolute error of Re-COS does not exceed 0.0017. We observe that the effect of $\lambda$ on the error is negligible, while the effect of $T$ on the error is significant. Both the absolute error and relative error of Re-COS increases with the decrease of $T$ .

Furthermore, we measure the convergence of the above two methods by the relative error. We approximate an American put option by the Bermudan put options with daily exercise (direct approximation) and take the obtained price as the benchmark. We compute the relative error of the Re-COS and Re-CONV method under $T=1/12$ (short term), 1/2 (middle term), and 1 (long term), respectively, and compare the convergence of the two methods with grid points. We specify $\lambda =10$ , $K=110$ . Other parameters are the same with the ones used in Table 1. For the direct approximation method, we use exercise opportunities 30, 128, and 256 under $T=1/12$ , 1/2, and 1, respectively and grid points $N=2^{15}$ under each maturity. For Re-COS, we use $L=13$ . For Re-CONV, we use $L=3$ and dampening factor 0.5. For four-point Richardson extrapolation schemes, we use exercise opportunities: $(2^5,2^4,2^3,2^2)$ , $(2^7,2^6,2^5,2^4)$ , and $(2^8,2^7,2^6,2^5)$ under $T=1/12$ , 1/2, and 1, respectively. Figure 2 summarizes the main results. The CPU time wasted by the two methods with grid points is also presented in Figure 2.

Figure 2.

The comparison of convergence and CPU time between the Re-COS method and the Re-CONV method with grid points for pricing American put options. We use (2⁵, 2⁴, 2³, 2²), (2⁷, 2⁶, 2⁵, 2⁴), (2⁸, 2⁷, 2⁶, 2⁵) in four-point Richardson extrapolation schemes under T = 1/12, 1/2, 1, respectively. The benchmark prices are computed by the Bermudan put options with daily exercise, that is, 30, 128, 256 exercise opportunities under T = 1/12, 1/2, 1, respectively. Re-CONV: Richardson extrapolation schemes and the CONV method; Re-COS: Richardson extrapolation schemes and the COS method.

From Figure 2, we see that both the two methods present smooth and stable convergence under each maturity. The Re-COS method converges with fewer grid options compared to the Re-CONV method. It only needs $2^6$ and $2^7$ grid points to obtain stable convergence for the Re-COS method under $T=1/12$ and $T=1/2$ , respectively, while it needs $2^9$ grid points to obtain the same convergence for the Re-CONV method. The CPU time wasted by the Re-COS method is little longer than the Re-CONV method for more than $2^9$ grid options. Since the Re-COS method has converged with $2^7$ grid points, the disadvantage of the Re-COS method in time can be negligible. Table 1 and Figure 2 verify that the Re-COS method is efficient for pricing American put options, especially for short-term American options.

Conclusion

The double Heston model with double exponential jumps incorporates several important features of stock return. We derive the characteristic function of the model. Under the model, we approximate American options by combining Richardson extrapolation technique and the COS method and provide a pricing algorithm for American put options. Numerical results show that the presented algorithm is efficient for pricing American put options, especially for short-term American put options. Although the values of American call options can be obtained by the put-call symmetry, our numerical results (not reported in the paper) show the performance of our algorithm for call options is not as good as for put options.