# A Shannon Wavelet Method for Pricing American Options under Two-Factor Stochastic Volatilities and Stochastic Interest Rate

Authors: Huang Shoude; Xunxiang Guo
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In the paper, the pricing of the American put options under the double Heston model with Cox–Ingersoll–Ross (CIR) interest rate process is studied. The characteristic function of the log asset price is derived, and thereby Bermuda options are well evaluated by means of a state-of-the-art Shannon wavelet inverse Fourier technique (SWIFT), which is a robust and highly efficient pricing method. Based on the SWIFT method, the price of American option can be approximated by using Richardson extrapolation schemes on a series of Bermudan options. Numerical experiments show that the proposed pricing method is efficient, especially for short-term American put options.


---

## Body

## 1. Introduction

It is well known that both the academic researchers and financial engineers have been concerned with the study of option pricing theory which is one of the most important and useful research fields in financial engineering. The well-known Black–Scholes–Merton (BS) model has been widely applied in the financial industry since it was proposed in the 1970s [1]. Unfortunately, the model has many shortcomings, such as logarithm yield of the stock price are assumed to be normal distribution and volatility of stock price is assumed to be constant which is different from the real data from the financial market. This is an important reason that brings a lot of economic losses to many financial institutions and even leads to the bankruptcy of some financial institutions. And, to make matters worse, the stability and sustainability of the whole financial system are influenced. Therefore, in order to maintain the stability and sustainability of the country’s financial system, the government implemented tighter financial regulation and many scholars continue to study the option pricing theory.Considering the abovementioned shortcomings, many researchers generalized the BS model by introducing stochastic volatility among which the popular models are the models raised by Stein and Stein [2], Heston [3], and Schöbel and Zhu [4]. However, single-factor stochastic volatility models are unable to fit the implied volatility smile very well. Many empirical research studies, such as Cont and Tankov [5], Fonseca et al. [6], Christoffersen et al. [7], and Fouque and Lorig [8] show that single-factor models can do a bad job in capturing the term structures of implicit volatilities. Christoffersen et al. [7] proposed the double Heston model by adding another uncorrelated volatility process which can provide better empirical fit to the market price and more flexibility to model the volatility term structure than the Heston model. In addition, spot interest rate is a fundamental economic variable in the financial market and it cannot be treated as a constant. Some researchers introduce stochastic interest rate into option pricing models in recent years. Grzelak and Oosterlee [9] examined the correlation issues of European product pricing with the Heston–Hull–White and Heston-CIR hybrid models [9]. Recently, Kim et al. [10] showed that incorporation of stochastic interest rates into a stochastic volatility model gave better results compared with the constant interest rate case in any maturity. It has been verified that considering stochastic interest rate into the Heston model improves the pricing performance of the model for European options in Bakshi et al. [11]. Consequently, the model which combines the two-factor stochastic volatilities and stochastic interest rate is proposed.Because many exchange-traded derivative contracts are of American style, accurate and fast pricing of American options is important, particularly, for high frequency traders. But the pricing of American options is one of the most challenging problems. The essential difficulty is that the pricing of American options involves a difficulty-free boundary problem in mathematics. There is no exact formula for evaluating most of the American options, and some numerical pricing techniques can be developed. The existing numerical methods under stochastic volatility models can be classified into three major groups: the Monte Carlo simulation, partial integro-differential equation (PIDE) methods, and numerical integration methods. One important type of numerical integration is the so-called Fourier transform method, which has attracted more and more attention since it is efficient and often used for calibration purposes. One of the best known examples of this class is the Carr–Madan method [12] based on fast Fourier transform. Zhylyevskyy [13] also evaluated American options under stochastic volatility by the fast Fourier transform technique. However, fast Fourier transform is usually unstable affected by the dampening factor. A more efficient pricing method based on Fourier cosine series expansion is proposed by Fang and Oosterlee [14]. The main achievement of this method is that it maintains an exponential convergence rate. But, it brings forth disadvantages in the vicinity of integration boundaries because cosines exhibit periodicity behavior. And, round-off errors may accumulate near domain boundaries for long-maturity options. Gracia and Oosterlee [15] further apply the Shannon wavelet inverse Fourier technique (SWIFT) method to price European option under exponential Lévy dynamics, and SWIFT was proved to be a robust, accurate, and highly efficient financial option valuation method. In this paper, we firstly utilize the SWIFT method for pricing Bermudan options, and then the price of American option can be approximated by using Richardson extrapolation schemes on a series of Bermudan options.An outline of the paper is as follows. In Section2, we introduce the basic model specification. The characteristic function of the underlying log-asset price is also derived. In Section 3, we briefly introduce the SWIFT method and the procedures to price the American put option. We present a set of numerical results and analysis in Section 4 to evaluate the accuracy and efficiency of the SWIFT method benchmarked to the Monte Carlo method. We conclude the paper in Section 5.


## 2. Model Specification and Derivation of the Characteristic Function

Assume thatΩ,ℱ,ℱt0≤t≤T,P denotes a complete probability space with a filtration continuous on the right, where P is a risk-neutral probability. Suppose that no dividends are paid in the market and the underlying asset can be traded in continuous time and without transaction cost. Let St denote the underlying asset price process. Let W1,ts,W2,ts,W1,tv,W2,tv, and Wtr be Brownian motions. Moreover, CovdW1,ts,dW1,tv=ρ1dt,CovdW2,ts,dW2,tv=ρ2dt; any other two processes are uncorrelated with each other. The double Heston with CIR interest rate process model assumes the following dynamics:(1)dStSt=rtdt+V1,tdW1,ts+V2,tdW2,ts,(2)dV1,t=k1θ1−V1,tdt+σ1V1,tdW1v,t,(3)dV2,t=k2θ2−V2,tdt+σ2V2,tdW2v,t(4)drt=krθr−rtdt+σrrtdWtrwhere kj,θj,σjj=1,2 mean the speeds of the mean reversion, long-run volatility, and instantaneous volatility of process Vj,tj=1,2, respectively. The stochastic volatility processes Vj,tj=1,2 satisfy the Feller conditions [16], i.e., 2kjθj≥σj2, to make Vj,tj=1,2 to be positive. In addition, kr,θr,σr are the mean reversion speed, long-term interest rate, and volatility of the interest rate process rt, respectively. Let us suppose V1,0=V1,V2,0=V2,r0=r,S0=S.Remark 1.
(1) The BS model is a special case of the above model by settingk1=θ1=σ1=0, k2=θ2=σ2=V2=ρ2=0, and kr=θr=σr=0. (2) The Heston model is also a special case of the above model by setting kr=θr=σr=0 and k2=θ2=σ2=V2=ρ2=0. (3) The double Heston is another special case of the above model by setting kr=θr=σr=0. According to Brigo and Mercurio [17], when the short-term interest rate rt evolves according to equation (4), the price of a risk-free zero-coupon bond Pt,T with maturity T at time 0 is given by(5)P0,T=eAT+BTrwhere(6)AT=kr−ξkrθrTσr2+2krθrσr2ln2ξ2ξ+kr−ξ1−e−ξT,BT=2e−ξT−12ξ+kr−ξ1−e−ξT,ξ=kr2+2σr2.
The conditional characteristic function of the log asset price process is the key component in SWIFT to price the options. We can simplify calculation and change from the numeraire measureP to the T-forward measure Q by applying the following Radon–Nikodym derivative:(7)dQdP=exp−∫0TrtdtP0,T.
Then, we derive the conditional characteristic function of the log asset price process in the following result.Theorem 1.
Suppose the price of the underlying asset follows the dynamics (1)–(4). Let τ=T−t,xT=lnST. Under the T-forward measure Q, the conditional characteristic function ϕ⋅ of xt is given by(8)ϕx,V1,V2,r,τ;u=eiux+Aτ,u+B1τ,uV1+B2τ,uV2+Cτ,ur,where(9)Aτ,u=−k1θ1α2ζ+α12τ+lnζ−α1+ζ−α1e−ζτ2ζ−k2θ2β2η+β12τ+lnη−β1+η+β1e−ητ2η−krθrγ2Λ+γ12τ+lnΛ−γ1+Λ−γ1e−Λτ2Λ,B1τ,u=2α01−e−ζτζ−α1+ζ−α1e−ζτ,B2τ,u=2β01−e−ητη−β1+η−β1e−ητ,Cτ,u=2γ01−e−ΛτΛ−γ1+Λ−γ1e−Λτ,α2=12σ12,α1=iuρ1σ1−k1,α0=−12ui+u,ζ=α12−4α0α2,β2=12σ22,β1=iuρ2σ2−k2,β0=−12ui+u,η=β12−4β0β2,γ2=12σr2,γ1=−kr,γ0=iu−1,Λ=γ12−4γ0γ2.Proof.
By applying Feynman-Kac theorem,ϕ⋅ satisfies the following PDE:(10)−rϕ+∂ϕ∂t+r−V1+V22∂ϕ∂x+V1+V22∂2ϕ∂x2+∑j=1212σj2Vj∂2ϕ∂Vj2+ρjσjVj∂2ϕ∂xVj+kjθj−Vj∂ϕ∂Vj+krθr−rt∂ϕ∂r+12σr2r∂2ϕ∂r2=0ϕx,V1,V2,r,0;u=eiux
According to Duffie et al. [18], we conjecture ϕ⋅ has the following form:(11)ϕx,V1,V2,r,τ;u=eiux+Aτ,u+B1τ,uV1+B2τ,uV2+Cτ,ur,with initial conditions A0,u=B10,u=B20,u=C0,u=0. Substituting (22) into (21) yields(12)−r+∂A∂t+∂B1∂tV1+∂B1∂tV2+∂C∂tr−12V1+V2u2+12σ12V1∂B1∂t2+12σ22V2∂B2∂t2+12σr2rC2+ρ1σ1V1iuB1+ρ2σ2V2iuB2+k1θ1−V1B1+k2θ2−V2B2+krθr−rC+r−12V1+V2iu=0.
From the above equation, we have a system of four ordinary differential equations as follows:(13)−∂B1∂τ+12σ12B12+iuρ1σ1−k1B1−12ui+u=0,−∂B2∂τ+12σ22B22+iuρ2σ2−k2B2−12ui+u=0,−∂C∂τ+iu+12σr2C2−Ckr−1=0,−∂A∂τ+B1k1θ˜1+B2k2θ2+Ckrθr=0.
Then, we can obtain(14)Aτ,u=−k1θ1α2ζ+α12τ+lnζ−α1+ζ+α1e−ζτ2ζ−k2θ2β2η+β12τ+lnη−β1+η+β1e−ητ2η−krθrγ2Λ+γ12τ+lnΛ+γ1+Λ+γ1e−Λτ2ΛB1τ,u=2α01−e−ζτζ+α1+ζ+α1e−ζτ,B2τ,u=2β01−e−ητη−β1+η+β1e−ητ,Cτ,u=2γ01−e−ΛτΛ−γ1+Λ+γ1e−Λτ,α2=12σ12,α1=iuρ1σ1−k1,α0=−12ui+u,ζ=α12−4α0α2,β2=12σ22,β1=iuρ2σ2−k2,β0=−12ui+u,η=β12−4β0β2,γ2=12σr2,γ1=−kr,γ0=iu−1,Λ=γ12−4γ0γ2.


## 3. American Option Pricing with Shannon Wavelet Expansions

### 3.1. Shannon Wavelet Expansions for European Options

Before we discuss the American options, we firstly introduce briefly the SWIFT method for European pricing problems as in Gracia and Oosterlee [15]. The price of a European option can be found by means of the Feynman–Kac formula as the discounted expectation of the option value at final time T. According to the risk-neutral option valuation formula,(15)vx,t=Ee∫tTrsdsvy,Tx=Pt,T∫ℝvy,Tfyxdy,where(16)vy,T=gy=αKey−1+,α=1,for a call,−1,for a put,where v denotes the option value, T is the maturity, t is the initial date, Δt≔T−t is the time to maturity, E is the risk-neutral expectation, x and y are state variables at time t and T respectively, and fyx is the probability density of y given x. As can be seen from the above expression, the key to price European option is to calculate the conditional density function. The SWIFT method can approximate the density function by recovering it from its characteristic function f^ω;x in three steps.In the first step,f is approximated by Shannon wavelet, that is,(17)fyx≈f1yx≔Pmfyx=∑k∈ZDm,kxDm,ky,where Dm,kx≔f⋅x,ϕm,k⋅, ϕm,ky is the Shannon wavelet.In the second step, the infinite summation in equation (25) has to be truncated for some κ=2mc+1, thus we obtain(18)f1yx≈f2yx≔∑k=1−κκDm,kxϕm,ky.The final step is to compute the density coefficientsDm,kx, i.e.,(19)Dm,kx≈Dm,k∗x≔∫ℝfyxϕm,ky¯dy≈2m/2J∑j=1JRe∫ℝfyxe−iωj2my−kdy=2m/2J∑j=1JRef^ωj2m;xeiωjk,where J≔⌈log2πN⌉. How to compute the vector of coefficients Dm,kx≔f⋅x,ϕm,k⋅ efficiently using the fast Fourier transform has been showed in by Maree [19]. To sum up, by means of the SWIFT method, the resulting option value at time t can be written as follows:(20)vx,t≈v0x,t≔Pt,T∫y≤cfyxvy,Tdy,(21)≈Pt,T∫y≤cf2yxvy,Tdy,(22)=Pt,T∑k=1−κκDm,k∗x∫y≤cvy,Tϕm,kydy.Let(23)Vm,kT≔∫y≤cvy,Tϕm,kydy.It is closely related to the wavelet coefficients of the value functionvy,T. For put options, the following approximate solutions can be obtained:(24)Vm,kTVm,k∗T≔−K2m/22J−1∑j=12J−1I1k12m,k¯22m−I2k12m,k¯22m,ifk1<0,0,ifk1≥0,where detailed expression of I1k1/2m,k¯2/2m−I2k1/2m,k¯2/2m can be found in [15]. Thus, we can obtain the following result:(25)vx,t≈Pt,T∑k=1−κκDm,k∗xVm,k∗T.


### 3.2. Shannon Wavelet Expansions for Bermudan Options

According to Chang and Chung [20], the prices of American options can be obtained by applying a Richardson extrapolation on the prices of a small number of Bermudan options. So, we need obtain the price of the Bermudan option firstly.A Bermudan option is a financial contract that allows the holder to exercise the option for a predetermined period of times before the option expires. Without loss of generality, we can assume a Bermudan option with strike priceK and a set of N exercise moments t1,…,tN, and strictly ordered, 0=t0<t1<t2<⋯<tN=T, where T is the option’s maturity. In the case of Bermudan option pricing, the dynamic programming principle reveals that the price of the Bermudan option at any exercise moment is the maximum of the spot payoff and the so-called continuation value. The valuation process can be regarded as a European option pricing problem between the two exercise times and can be priced by means of risk-neutral option pricing formula. For simplicity, we can divide time grid Δtn=tn−tn−1=Δt and we define x≔Xtn−1=logStn−1/K, and y≔Xtn=logStn−1/K, where St is the price process of the underlying asset price process. The payoff of the option is given by gy, and for the European option, the value of the option at maturity is given by(26)vy,T=gy=αKey−1+,α=1,for a call,−1,for a put.It is worth noting that we consider the stochastic interest rate in this article. By the dynamic programming method, the option value ahead of maturity can be expressed recursively forn=N,N−1,…,2,1, by(27)vx,tn−1=maxgx,cx,tn−1,cx,tn−1=P0,Δt∫ℝvy,tnfy|xdy,and the option value at t0 is finally followed by(28)vx,t0=P0,Δt∫ℝvy,t1fy|xdywhere cx,tn−1 is the continuation value and the probability function of y given x is written as fyx. The European option value is given by the integral representation in (24). According to the SWIFT pricing formula (26), we obtain(29)vx,t0≈v∗x,t0=P0,Δt∑k=1−κκDm,k∗xVm,kt1.So, the option value can be known once the value coefficientsVm,kt1 are determined. According to the characterization of Bermudan options, a backward recursion will be applied to recover these coefficients Vm,ktjj=1N−1. We know that the option value at maturity tN=T equals the payoff of the option, i.e., vy,T=gy. Thus, the value coefficients for a put option is given by(30)Vm,ktN=Vm,kT≔∫y≤cvy,Tϕm,kydy=∫−c0gyϕm,kydy=Gm,k−c,0,where Gm,ky1,y2 means the payoff coefficients, and it can be efficiently computed using the fast Fourier transform as Oosterlee mentioned in [19]. Once a formula for the value coefficients at maturity is obtained, we can calculate the coefficients at any time tn, for n=N−1,…,1, ahead of maturity recursively. When the early-exercise point xn∗ is known, which is the point where the continuation value equals the payoff, i.e., cxn∗,tn=gxn∗, where(31)cx,tn−1≈c∗x,tn−1≔e−rΔt∑j=1JRef^ωj2mU˜jtne−iωj2mx,U˜jtn≔2m/2J∑p≤κVm,ptneiωjp.Oncexn∗ is obtained, the integral for the value coefficients Vm,k⋅ can be splitted into two parts, i.e.,(32)Vm,ktn−1≔∫x≤cvx,tn−1ϕm,kxdx=∫x≤cmaxgx,cx,tn−1ϕm,kxdx=∫−cxn∗gxϕm,kxdx+∫xn∗ccx,tn−1ϕm,kxdx≕Gm,k−c,xn∗+Cm,kxn∗,c,tn−1,where Cm,kx1,x2,tn−1 denotes the continuation coefficients at time tn−1 over the interval x1,x2. It also can be efficiently approximated when the value coefficients at the next step Vm,ktn are known, i.e.,(33)Cm,k∗x1,x2,tn−1≔e−rΔtRe∑q=1−JJJqx1,x2eiωqk,where(34)Jqx1,x2≔∑i=1Jf^ωj2mU˜jtnℳj,qx1,x2,with(35)ℳj,qx1,x2=x2−x1,forq=1−j,ie−iωj+ωq2mx2−e−iωj+ωq2mx1ωj+ωq2m,else.We can see thatℳj,q is a Hankel matrix and Jq represents a matrix-vector product. We can compute this matrix-vector product by means of FFT since the matrix-vector product ℳx is equal to the first N elements of the circular convolution mh⊛xh, with the 2N-vectors(36)mh≔m0,m−1,m−2,…,m1−N,0,mN−1,mN−2,…,m1T,xh≔x0,x1,x2,…,xN1,0,…,0T.At the same time, a circular convolution of two vectors is equal to the inverse discrete Fourier transformDFT−1 of the product of the forward DFTs, i.e.,(37)x⊛y=DFT−1DFTx⋅DFTy,Vm,ktn−1=Gm,k∗−c˜,xn∗+Cm,k∗xn∗,c˜,tn.Thus, according to the abovementioned formulas,Vm,kt1 can be obtained by the backward recursion. Between t0 and t1, the valuation of the Bermudan option process can be regarded as a European option pricing problem and can be priced with the help of the risk-neutral option valuation formula (15). Thus, the price of Bermudan option at t0 can be obtained as follows:(38)vx,t0≈v∗x,t0=P0,Δt∑k=1−κκDm,k∗xVm,kt1.The prices of American options can be obtained by applying a Richardson extrapolation on the prices of a few Bermudan options [18], which is also demonstrated in [19]. Let vM denote the value of a Bermudan option with M early exercise dates. We will use the following 4-point Richardson extrapolation scheme:(39)vAMd=12164v2d+3−56v2d+2+14v2d+1−v2d,where d is a positive integer and vAMd denotes the approximated value of the American option. So, we can approximate the price of an American option by four Bermudan options with certain number of exercise dates.


## 3.1. Shannon Wavelet Expansions for European Options

Before we discuss the American options, we firstly introduce briefly the SWIFT method for European pricing problems as in Gracia and Oosterlee [15]. The price of a European option can be found by means of the Feynman–Kac formula as the discounted expectation of the option value at final time T. According to the risk-neutral option valuation formula,(15)vx,t=Ee∫tTrsdsvy,Tx=Pt,T∫ℝvy,Tfyxdy,where(16)vy,T=gy=αKey−1+,α=1,for a call,−1,for a put,where v denotes the option value, T is the maturity, t is the initial date, Δt≔T−t is the time to maturity, E is the risk-neutral expectation, x and y are state variables at time t and T respectively, and fyx is the probability density of y given x. As can be seen from the above expression, the key to price European option is to calculate the conditional density function. The SWIFT method can approximate the density function by recovering it from its characteristic function f^ω;x in three steps.In the first step,f is approximated by Shannon wavelet, that is,(17)fyx≈f1yx≔Pmfyx=∑k∈ZDm,kxDm,ky,where Dm,kx≔f⋅x,ϕm,k⋅, ϕm,ky is the Shannon wavelet.In the second step, the infinite summation in equation (25) has to be truncated for some κ=2mc+1, thus we obtain(18)f1yx≈f2yx≔∑k=1−κκDm,kxϕm,ky.The final step is to compute the density coefficientsDm,kx, i.e.,(19)Dm,kx≈Dm,k∗x≔∫ℝfyxϕm,ky¯dy≈2m/2J∑j=1JRe∫ℝfyxe−iωj2my−kdy=2m/2J∑j=1JRef^ωj2m;xeiωjk,where J≔⌈log2πN⌉. How to compute the vector of coefficients Dm,kx≔f⋅x,ϕm,k⋅ efficiently using the fast Fourier transform has been showed in by Maree [19]. To sum up, by means of the SWIFT method, the resulting option value at time t can be written as follows:(20)vx,t≈v0x,t≔Pt,T∫y≤cfyxvy,Tdy,(21)≈Pt,T∫y≤cf2yxvy,Tdy,(22)=Pt,T∑k=1−κκDm,k∗x∫y≤cvy,Tϕm,kydy.Let(23)Vm,kT≔∫y≤cvy,Tϕm,kydy.It is closely related to the wavelet coefficients of the value functionvy,T. For put options, the following approximate solutions can be obtained:(24)Vm,kTVm,k∗T≔−K2m/22J−1∑j=12J−1I1k12m,k¯22m−I2k12m,k¯22m,ifk1<0,0,ifk1≥0,where detailed expression of I1k1/2m,k¯2/2m−I2k1/2m,k¯2/2m can be found in [15]. Thus, we can obtain the following result:(25)vx,t≈Pt,T∑k=1−κκDm,k∗xVm,k∗T.


## 3.2. Shannon Wavelet Expansions for Bermudan Options

According to Chang and Chung [20], the prices of American options can be obtained by applying a Richardson extrapolation on the prices of a small number of Bermudan options. So, we need obtain the price of the Bermudan option firstly.A Bermudan option is a financial contract that allows the holder to exercise the option for a predetermined period of times before the option expires. Without loss of generality, we can assume a Bermudan option with strike priceK and a set of N exercise moments t1,…,tN, and strictly ordered, 0=t0<t1<t2<⋯<tN=T, where T is the option’s maturity. In the case of Bermudan option pricing, the dynamic programming principle reveals that the price of the Bermudan option at any exercise moment is the maximum of the spot payoff and the so-called continuation value. The valuation process can be regarded as a European option pricing problem between the two exercise times and can be priced by means of risk-neutral option pricing formula. For simplicity, we can divide time grid Δtn=tn−tn−1=Δt and we define x≔Xtn−1=logStn−1/K, and y≔Xtn=logStn−1/K, where St is the price process of the underlying asset price process. The payoff of the option is given by gy, and for the European option, the value of the option at maturity is given by(26)vy,T=gy=αKey−1+,α=1,for a call,−1,for a put.It is worth noting that we consider the stochastic interest rate in this article. By the dynamic programming method, the option value ahead of maturity can be expressed recursively forn=N,N−1,…,2,1, by(27)vx,tn−1=maxgx,cx,tn−1,cx,tn−1=P0,Δt∫ℝvy,tnfy|xdy,and the option value at t0 is finally followed by(28)vx,t0=P0,Δt∫ℝvy,t1fy|xdywhere cx,tn−1 is the continuation value and the probability function of y given x is written as fyx. The European option value is given by the integral representation in (24). According to the SWIFT pricing formula (26), we obtain(29)vx,t0≈v∗x,t0=P0,Δt∑k=1−κκDm,k∗xVm,kt1.So, the option value can be known once the value coefficientsVm,kt1 are determined. According to the characterization of Bermudan options, a backward recursion will be applied to recover these coefficients Vm,ktjj=1N−1. We know that the option value at maturity tN=T equals the payoff of the option, i.e., vy,T=gy. Thus, the value coefficients for a put option is given by(30)Vm,ktN=Vm,kT≔∫y≤cvy,Tϕm,kydy=∫−c0gyϕm,kydy=Gm,k−c,0,where Gm,ky1,y2 means the payoff coefficients, and it can be efficiently computed using the fast Fourier transform as Oosterlee mentioned in [19]. Once a formula for the value coefficients at maturity is obtained, we can calculate the coefficients at any time tn, for n=N−1,…,1, ahead of maturity recursively. When the early-exercise point xn∗ is known, which is the point where the continuation value equals the payoff, i.e., cxn∗,tn=gxn∗, where(31)cx,tn−1≈c∗x,tn−1≔e−rΔt∑j=1JRef^ωj2mU˜jtne−iωj2mx,U˜jtn≔2m/2J∑p≤κVm,ptneiωjp.Oncexn∗ is obtained, the integral for the value coefficients Vm,k⋅ can be splitted into two parts, i.e.,(32)Vm,ktn−1≔∫x≤cvx,tn−1ϕm,kxdx=∫x≤cmaxgx,cx,tn−1ϕm,kxdx=∫−cxn∗gxϕm,kxdx+∫xn∗ccx,tn−1ϕm,kxdx≕Gm,k−c,xn∗+Cm,kxn∗,c,tn−1,where Cm,kx1,x2,tn−1 denotes the continuation coefficients at time tn−1 over the interval x1,x2. It also can be efficiently approximated when the value coefficients at the next step Vm,ktn are known, i.e.,(33)Cm,k∗x1,x2,tn−1≔e−rΔtRe∑q=1−JJJqx1,x2eiωqk,where(34)Jqx1,x2≔∑i=1Jf^ωj2mU˜jtnℳj,qx1,x2,with(35)ℳj,qx1,x2=x2−x1,forq=1−j,ie−iωj+ωq2mx2−e−iωj+ωq2mx1ωj+ωq2m,else.We can see thatℳj,q is a Hankel matrix and Jq represents a matrix-vector product. We can compute this matrix-vector product by means of FFT since the matrix-vector product ℳx is equal to the first N elements of the circular convolution mh⊛xh, with the 2N-vectors(36)mh≔m0,m−1,m−2,…,m1−N,0,mN−1,mN−2,…,m1T,xh≔x0,x1,x2,…,xN1,0,…,0T.At the same time, a circular convolution of two vectors is equal to the inverse discrete Fourier transformDFT−1 of the product of the forward DFTs, i.e.,(37)x⊛y=DFT−1DFTx⋅DFTy,Vm,ktn−1=Gm,k∗−c˜,xn∗+Cm,k∗xn∗,c˜,tn.Thus, according to the abovementioned formulas,Vm,kt1 can be obtained by the backward recursion. Between t0 and t1, the valuation of the Bermudan option process can be regarded as a European option pricing problem and can be priced with the help of the risk-neutral option valuation formula (15). Thus, the price of Bermudan option at t0 can be obtained as follows:(38)vx,t0≈v∗x,t0=P0,Δt∑k=1−κκDm,k∗xVm,kt1.The prices of American options can be obtained by applying a Richardson extrapolation on the prices of a few Bermudan options [18], which is also demonstrated in [19]. Let vM denote the value of a Bermudan option with M early exercise dates. We will use the following 4-point Richardson extrapolation scheme:(39)vAMd=12164v2d+3−56v2d+2+14v2d+1−v2d,where d is a positive integer and vAMd denotes the approximated value of the American option. So, we can approximate the price of an American option by four Bermudan options with certain number of exercise dates.


## 4. Numerical Examples

In this section, some numerical examples are performed to test the performance of the SWIFT method for pricing American puts against its alternative competitor, the LSM Monte Carlo method [21], as a benchmark. In particular, we make use of the same constant parameters as specified by Zhang and Zhang [22]. These parameter values are listed in Table 1 for all our numerical examples. The computer used in the experiments equips an Intel Core i7 CPU with a 2.2 GHz processor. And, all of these numerical examples were performed in Matlab 2016a. For the SWIFT method, unless otherwise mentioned, we use L = 10, where L is the parameter used to define the truncation range [a, b] as follows:(40)c≔c1+Lc2+c4,where cnn=1,…,4 is the nth cumulant of log-price process logS. It can be computed by(41)cn=1in∂nlnϕu∂unu=0.Table 1
Parameter values for the numerical analyses.
Parameter
k
1
θ
1
σ
1
v
1
ρ
1
k
2
θ
2
σ
2
Value
0.9
0.1
0.1
0.6
−0.5
1.2
0.15
0.2
Parameter
v
2
ρ
2
k
r
θ
r
σ
r
s
r
—
Value
0.7
−0.5
0.26
0.04
0.08
10
0.03
—The relative error and CPU time information, comparing the SWIFT and the LSM method, are presented in Table2 for the pricing of American put options. Table 2 compares relative errors between the SWIFT method and LSM for pricing American put options when d=5 in equation (27), and it shows that the SWIFT method is significantly faster than LSM. Furthermore, Table 2 compares the pricing accuracy among the two methods across a range of strike prices, and the relative percentage pricing differences of SWIFT are all less than 0.13%. Our numerical results show that the SWIFT method is fast and accurate for pricing American put options.Table 2
Comparisons of the CPU time and accuracy for of SWIFT and LS methods.
K
T
SWIFT
LSM (s)
Abs relative error (%)
8
1
/
12
0.4143
0.4146
±
0.0007
0.07
8
1
/
4
1.0533
1.0523
±
0.0014
0.10
8
1
/
2
1.6220
1.6223
±
0.0018
0.02
10
1
/
12
1.2711
1.2725
±
0.0013
0.11
10
1
/
4
2.0856
2.0844
±
0.0020
0.06
10
1
/
2
2.7570
2.7601
±
0.0024
0.11
12
1
/
12
2.6051
2.6072
±
0.0018
0.08
12
1
/
4
3.4101
3.4056
±
0.0025
0.13
12
1
/
2
4.1064
4.1044
±
0.0030
0.05
CPU time
—
0.1936 s
2011.6351
—


## 5. Conclusion

This article proposed an efficient method for pricing American put options under double stochastic volatility with stochastic interest rate. The conditional characteristic function of the log-asset price is derived, and an analytical approximation method for evaluating American puts is proposed. Numerical results of this article confirm that the SWIFT method is efficient and accurate for pricing the American put options. The future research is mainly about parameter calibration and hedge strategy.



---

*Source: 8531959-2020-04-09.xml*