# A Finite Difference Scheme for Pricing American Put Options under Kou's Jump-Diffusion Model

Authors: Jian Huang; Zhongdi Cen; Anbo Le
Journal: Journal of Function Spaces and Applications
 (2013)

---

## Abstract

We present a stable finite difference scheme on a
piecewise uniform mesh along with a penalty method for pricing American
put options under Kou's jump-diffusion model. By adding a penalty term,
the partial integrodifferential complementarity problem arising from pricing
American put options under Kou's jump-diffusion model is transformed into
a nonlinear parabolic integro-differential equation. Then a finite difference
scheme is proposed to solve the penalized integrodifferential equation, which
combines a central difference scheme on a piecewise uniform mesh with respect
to the spatial variable with an implicit-explicit time stepping technique. This
leads to the solution of problems with a tridiagonal M-matrix. It is proved
that the difference scheme satisfies the early exercise constraint. Furthermore,
it is proved that the scheme is oscillation-free and is second-order convergent
with respect to the spatial variable. The numerical results support the theoretical
results.


---

## Body

## 1. Introduction

It is widely recognized that the assumption of log-normal stock diffusion with constant volatility in the standard Black-Scholes model [1] of option pricing is not ideally consistent with that of the market price movement. In particular, the probability distribution of realized asset returns often exhibits features that are not taken into account by the standard Black-Scholes model: heavy tails, volatility clustering, and volatility smile [2]. In order to explain these phenomena, extensions of the Black-Scholes model have been proposed. Generally speaking, two different classes of models have been studied in the finance literature: the stochastic volatility models [3, 4] and the jump-diffusion models [2, 5]. Contrary to the Black-Scholes model, the jump diffusion models allow for a more realistic representation of price dynamics and greater flexibility in modeling. During the last twenty years, research on models with jumps has become very active. Most of such models have been proposed; see [2] and references therein. Here we focus on a jump-diffusion model with finite jump activity proposed by Kou in [6].Unlike the standard Black-Scholes equation, the valuation of options under jump-diffusion models requires solving a partial integrodifferential equation. A fully implicit scheme would lead to full matrices due to the integral term, which makes many methods computationally too expensive. Several numerical methods based on the finite difference method have been proposed for pricing options under jump-diffusion models. Amin [7] gave a multinomial tree method for pricing options under jump-diffusion models, which is actually an explicit type finite difference approach. Zhang [8] and Cont and Voltchkova [9] used implicit-explicit finite difference methods for pricing options under jump-diffusion models. Andersen and Andreasen [10] and Almendral and Oosterlee [11] proposed operator splitting methods coupled with a fast Fourier transformation (FFT) technique for pricing options with jump diffusion processes. d'Halluin et al. [12, 13] developed a second-order accurate numerical method with a fixed-point iteration method and an implicit finite difference scheme along with a penalty method for pricing American options under jump diffusion processes. Toivanen et al. [14–16] introduced a high-order front-fixing finite difference method and an artificial volatility scheme along with an iterative method for pricing American options under jump-diffusion models. Zhang and Wang [17, 18] proposed fitted finite volume schemes coupled with the Crank-Nicolson time stepping method for pricing options under jump diffusion processes.It is well known that the Black-Scholes partial differential operator atx=0 is degenerative. The Black-Scholes partial differential operator becomes a convection-dominated operator when the volatility or the asset price is small. Hence, numerical difficulty can be caused when the standard methods such as the central difference and piecewise linear finite element methods are used to solve those problems. A common and widely used approach by many authors dealing with finite difference/volume/element methods for the Black-Scholes partial differential equaion is to apply an Euler transformation to remove the singularity of the differential operator when the parameters of the Black-Scholes equation are constant or space-independent; see for example, [2, 19]. As a result of the Euler transformation, the transformed interval becomes (-∞,∞). However, the truncation on the left-hand side of the domain to artificially remove the degeneracy may cause computational errors. Furthermore, the uniform mesh on the transformed interval will lead to the originally grid points concentrating around x=0 inappropriately. Moreover, when a problem is space-dependent, this transformation is impossible, and thus the Black-Scholes equation in the original form needs to be solved [20]. The same problem also appears in the partial integral-differential equations resulting from jump-diffusion models [9, 11, 18]. Wang [21] and Angermann and Wang [22] applied a stable fitted finite volume method to deal with the degeneracy and singularity of the Black-Scholes operator. In this paper, we will present a stable finite difference method with a second-order convergency with respect to the spatial variable for solving the partial integrodifferential equation defined on (0,+∞) for arbitrary volatility and arbitrary interest rate.The penalty method was introduced by Zvan et al. [23] for pricing American options with stochastic volatility by adding a source term to the discrete equation. Nielsen et al. [24] presented a refinement of their work by adding a penalty term to the continuous equation and illustrated the performance of various numerical schemes. By adding a penalty term, the linear complementarity problem for pricing the American options can be transformed into a nonlinear parabolic partial differential equation. As the solution approaches the pay-off function at expiry, the penalty term forces the solution to stay above it. When the solution is far from the barrier, the term is small and thus the Black-Scholes equation is approximatively satisfied in this region.In [25] we have presented a robust difference scheme for the penalized Black-Scholes equation governing American put option pricing. In this paper we present a stable finite difference scheme on a piecewise uniform mesh along with a power penalty method for pricing American put options under Kou's jump-diffusion model. By adding a penalty term the partial integrodifferential complementarity problem arising from pricing American put options under Kou's jump-diffusion model is transformed into a nonlinear parabolic integrodifferential equation. Then a finite difference scheme is proposed to solve the penalized integrodifferential equation, which combines a central difference scheme on a piecewise uniform mesh with respect to the spatial variable with an implicit-explicit time stepping technique. This leads to the solution of problems with a tridiagonal M-matrix. It is proved that the difference scheme satisfies the early exercise constraint. Furthermore, it is proved that the scheme is oscillation-free and is second-order convergent with respect to the spatial variable. Numerical results support the theoretical results.The rest of the paper is organized as follows. In the next section, we describe some theoretical results on the continuous problem for pricing American put options under Kou's jump-diffusion model. The discretization method is described in Section3. In Section 4 we prove that the difference scheme satisfies the early exercise constraint. In Section 5, we present a stability and error analysis for the finite difference scheme. In Section 6, numerical experiments are provided to support these theoretical results. Finally, a discussion is indicated in Section 7.


## 2. The Continuous Problem

Letv denote the value of an American put option with strike price E on the underlying asset x and time t. It is known that the price v under a jump-diffusion model satisfies the following partial integrodifferential complementarity problem [14–17]:
(1)Lv(x,t)≥0,x>0,t∈[0,T),(2)v(x,t)-V*(x)≥0,x>0,t∈[0,T],(3)Lv(x,t)·[v(x,t)-V*(x)]=0,x>0,t∈[0,T),(4)v(x,T)=V*(x),x≥0,t=T,(5)v(0,t)=E,x=0,t∈[0,T],(6)v(x,t)→0,x→+∞,t∈[0,T],
where L denotes the partial integrodifferential operator defined by
(7)Lv(x,t)≡-∂v∂t-12σ2x2∂2v∂x2-(r-λζ)x∂v∂x+(r+λ)v-λ∫0∞‍v(xy,t)f(y)dy,V*(x) is the final (payoff) condition defined by
(8)V*(x)=max{E-x,0},σ is the volatility of the underlying asset, r is the risk free interest rate, t is the current time, T is the maturity date, and f(y) is the probability function of the jump amplitude y with the obvious properties that for all y, f(y)≥0 and ∫0∞‍f(y)dy=1, the constant ζ is given by ζ=∫0∞‍(y-1)f(y)dy. In Kou's model, f(y) is the following log-double-exponential density
(9)f(y)={qα2yα2-1,y<1,pα1y-α1-1,y≥1,
where p,q,α1>1 and α2 are positive constants such that p+q=1. It can be shown that, in this case, ζ=pα1/(α1-1)+qα2/(α2+1)-1. When the jump rate λ is zero, the partial integrodifferential operator reduces to the standard Black-Scholes operator [1]. In this paper, we assume that r-λζ≠0.The above linear complementarity problem (1)–(6) can be solved by a penalty approach. Let 0<ε≪1 be a small regularization parameter and consider the following initial-boundary value problem,
(10)Lv(x,t)-Cεv(x,t)+ε-q(x)=0,x>0,t∈[0,T),v(x,T)=V*(x),x≥0,t=T,v(0,t)=E,x=0,t∈[0,T],v(x,t)→0,x→+∞,t∈[0,T],
where C≥rE+λSmax is a positive constant and q(x)=E-x.By adding a penalty term(11)Cεv(x,t)+ε-q(x),
the linear complementarity problem for pricing American options can be transformed into a nonlinear parabolic integrodifferential equation. Essentially, it is of order ε in regions where v(x,t)≫q(x), and hence the partial integrodifferential equation is approximately satisfied. When v(x,t) approaches q(x), this term is approximately equal to C assuring that the early exercise constraint is not violated. For the continuous case, the convergence and the positivity constraint of the penalty method have been proved in [26]. In this paper, we consider a second-order finite difference scheme to discretize the semilinear partial integrodifferential equation (10) and prove that the approximate option values generated by the scheme satisfies a discrete version of (2).For applying the numerical method, we truncate the domain(0,+∞) into (0,Smax). Based on Wilmott et al.’s estimate [19] that the upper bound of the asset price is typically three or four times the strike price, it is reasonable for us to set Smax=4E. The boundary condition at x=Smax is chosen to be v(Smax,t)=0. Normally, this truncation of the domain leads to a negligible error in the value of the option [27].Therefore, in the remaining of this paper, we will consider the following nonlinear parabolic integrodifferential equation:(12)Lv(x,t)-Cεv(x,t)+ε-q(x)=0,(x,t)∈(0,Smax)×(0,T),(13)v(x,T)=V*(x),x∈[0,Smax],(14)v(0,t)=E,v(Smax,t)=0,t∈[0,T].


## 3. Discretization

We now consider the approximation of the solution to the semilinear partial integrodifferential equation (12)–(14) by a central difference scheme for the spatial derivatives.The use of central difference scheme on a uniform mesh may produce nonphysical oscillations in the computed solution. To overcome this oscillation, we use a piecewise uniform meshΩN on the space interval [0,Smax]:
(15)xi={h,i=1,h[1+σ2|r-λζ|(i-1)],i=2,…,N,
where
(16)h=Smax1+(σ2/|r-λζ|)(N-1).
For the time discretization, we use a uniform mesh ΩK on [0,T] with K mesh elements. Then the piecewise uniform mesh ΩN×K on Ω=(0,Smax)×(0,T) is defined to be the tensor product ΩN×K=ΩN×ΩK. It is easy to see that the mesh sizes hi=xi-xi-1 and τj=tj-tj-1 satisfy
(17)hi={h,i=1,σ2|r-λζ|h,i=2,…,N,(18)τ=τj=TK,j=1,…,K,
respectively.The space derivatives of (12) are approximated with central differences on the above piecewise-uniform mesh:
(19)∂v∂x(xi,t)≈vi+1(t)-vi-1(t)hi+hi+1,∂2v∂x2(xi,t)≈2hi+hi+1(vi+1(t)-vi(t)hi+1-vi(t)-vi-1(t)hi).The integral term(20)I=∫0∞‍v(xy,t)f(y)dy
of (12) can be approximated by a fast method as in [14–16]. By making the change of variable y=s/x, we have
(21)I=∫0∞‍v(xy,t)f(y)dy=∫0∞‍v(s,t)f(s/x)xds.
By using the linear interpolation, we can obtain an approximation
(22)Ii≈Ai(t)=∑n=0N-1‍∫xnxn+1‍(xn+1-shn+1v(xn,t)+s-xnhn+1v(xn+1,t))f(s/xi)xids
of I at each mesh point xi for i=1,2,…,N-1.A fully implicit scheme would lead to full matrices due to the integral term, which makes many methods computationally too expensive. Our technique is similar in some respects to Zhang [8], though less constrained in terms of stability restrictions. The integral term is treated explicitly in time, while the differential terms are treated implicitly. This leads to the solution of problems with a tridiagonal M-matrix. We will prove that the resulting time stepping method is unconditionally stable.Our implicit-explicit scheme to discretize the integrodifferential equation (12)–(14) is
(23)LN,KVij-CεVij+ε-qi=0,1≤i<N,0≤j<K,(24)ViK=V*(xi),0≤i≤N,(25)V0j=E,VNj=0,0≤j<K,
where
(26)LN,KVij≡-Vij+1-Vijτ-σ2xi2hi+hi+1(Vi+1j-Vijhi+1-Vij-Vi-1jhi)-(r-λζ)xiVi+1j-Vi-1jhi+hi+1+(r+λ)Vij-λAij+1,Aij+1=∑n=0N-1‍Vnj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍Vn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids.
Then from the above solution Vij, we can obtain the optimal stopping price which is the maximum asset price such that Vij=Vi* for each tj.Remark. Our method can be extended to the more general case of an infinite-activity process, like that of the Lévy-type models such as VG model [28] or CGMY model [29]. For example, the American put option price v under a generalized VG process satisfies the following partial integrodifferential complementarity problem [30, 31]:
(27)L~v(x,t)≥0,x>0,t∈[0,T),v(x,t)-V*(x)≥0,x>0,t∈[0,T],L~v(x,t)·[v(x,t)-V*(x)]=0,x>0,t∈[0,T),v(x,T)=V*(x),x≥0,t=T,v(0,t)=E,x=0,t∈[0,T],v(x,t)→0,x→+∞,t∈[0,T],
where L~ denotes the partial integrodifferential operator defined by
(28)L~v(x,t)≡-∂v∂t-12σ2x2∂2v∂x2-(r+w)x∂v∂x+rv-∫ℝ‍[v(xey,t)-v(x,t)]g(y)dy,w is some “compensation constant” given by
(29)w=∫ℝ‍(1-ey)g(y)dy,g(y) is known as the Lévy density. It is noted in Cont and Tankov [2] and Cont and Voltchkova [9] that for the Lévy densities (in the usual case where the density decays exponentially for large jumps) an infinite-activity process can be arbitrarily well approximated by a finite-activity process and an adjusted volatility. Having done this, our numerical method for the finite-activity process can be used to value options under infinite-activity processes. However, the matrix associated with a discrete operator may not be an M-matrix.


## 4. Positivity Constraint

In this section, we will prove that our scheme satisfies the early exercise constraint.Let(30)ai=σ2xi2(hi+hi+1)hi-(r-λζ)xihi+hi+1,bi=σ2xi2hihi+1,ci=σ2xi2(hi+hi+1)hi+1+(r-λζ)xihi+hi+1,i=1,…,N-1.
We can obtain
(31)ai>σ2x1xi(hi+hi+1)hi-(r-λζ)xihi+hi+1=σ2x1-(r-λζ)hi(hi+hi+1)hixi=σ2h-(r-λζ)(σ2/|r-λζ|)h(hi+hi+1)hixi≥0,1<i≤N-1,ci>σ2x1xi(hi+hi+1)hi+1+(r-λζ)xihi+hi+1=σ2x1+(r-λζ)hi+1(hi+hi+1)hi+1xi=σ2h+(r-λζ)(σ2/|r-λζ|)h(hi+hi+1)hi+1xi≥0,1≤i<N-1.Theorem 1.
The approximate option values generated by the implicit-explicit difference scheme (23)–(25) satisfy
(32)Vij≥V*(xi),(i,j)∈ΩN×K,
provided that C≥rE+λSmax.Proof.
We apply the similar technique of [24] to give the proof.
The schemes (23) can be written as
(33)(1τ+bi+r+λ)Vij=1τVij+1+aiVi-1j+ciVi+1j+λAij+1+CεVij+ε-qi.
The difference
(34)uij=Vij-qi
satisfies the following equation:
(35)(1τ+bi+r+λ)uij=1τuij+1+aiui-1j+ciui+1j+Cεuij+ε+λ(ζ+1)xi-(r+λ)E+λ(∑n=0N-1‍unj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍un+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids)+λ(∑n=0N-1‍qn·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍qn+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids).
Next, by defining
(36)uj=miniuij,
let k be an index such that
(37)ukj=uj.
For i=k, it follows from (35) that
(38)(1τ+bk+r+λ)uj≥1τukj+1+(ak+ck)uj+Cεuj+ε+λ(ζ+1)xk-(r+λ)E+λ(∑n=0N-1‍unj+1·∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=0N-1‍un+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xk)xkds)+λ(∑n=0N-1‍qn·∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=0N-1‍qn+1·∫xnxn+1‍s-xnhn+1·f(s/xk)xkds),
where we have used (31). Since
(39)∑n=0N-1‍qn·∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=0N-1‍qn+1·∫xnxn+1‍s-xnhn+1·f(s/xk)xkds>(E-Smax)·∫0xN‍f(s/xk)xkds=(E-Smax)·∫0xN/xk‍f(y)dy≥E-Smax,ζ+1=∫0∞‍(y-1)f(y)dy+1=∫0∞‍yf(y)dy>0,
we conclude that
(40)(1τ+r+λ)uj-Cεuj+ε+(rE+λSmax)≥1τukj+1+λ(∑n=1N-1‍unj+1·∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=1N-1‍un+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xk)xkds)≥1τuj+1+λuj+1(∑n=1N-1‍∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=1N-1‍∫xnxn+1‍s-xnhn+1·f(s/xk)xkds).
If we assume that
(41)uj+1≥0,
from (40) we can get
(42)F(uj)≥0,
where
(43)F(u)=(1τ+r+λ)u-Cεu+ε+(rE+λSmax).
Obviously,
(44)F′(u)=(1τ+r+λ)+Cε(u+ε)2≥0,F(0)=-C+(rE+λSmax)≤0,
where we have used the assumption C≥(rE+λSmax). Hence, from (42) we can obtain
(45)uj≥0.
Consequently, by induction on j, it follows from (34) that
(46)Vij≥qi,0≤i≤N,
for j=K,…,1,0.
Next we prove thatVij≥0. As above, we define
(47)Vj=miniVij,
and let k be an index such that
(48)Vkj=Vj.
For i=k, from (33) we obtain
(49)(1τ+bk+r+λ)Vkj≥1τVj+1+(ak+ck)Vj+λAkj+1+CεVkj+ε-qk.
Hence, we have
(50)(1τ+r+λ)Vkj≥1τVj+1+λAkj+1+CεVkj+ε-qk.
Since we have proved that Vij≥qi, we get
(51)(1τ+r+λ)Vkj≥1τVj+1+λ(∑n=0N-1‍Vnj+1·∫xnxn+1‍xn+1-shn+1·f(s/xk)xkds+∑n=0N-1‍Vn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xk)xkds).
By the above inequality. Hence, by induction on j, we can obtain
(52)Vij≥0,(i,j)∈ΩN×K,
for j=K,…,1,0.
Combine (46) with (52) to complete the proof.


## 5. Error Estimates

To investigate the convergence of the method, note that the error functionszij=Vij-vij (0≤i≤N,0≤j≤K) are the solutions of the discrete problem
(53)-zij+1-zijτ-σ2xi2hi+hi+1(zi+1j-zijhi+1-zij-zi-1jhi)-(r-λζ)xizi+1j-zi-1jhi+hi+1+(r+λ)zij+Cε(ξij+ε-qi)2zij-λ(∑n=0N-1‍znj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍zn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids)=Rij,1≤i<N,0≤j<K,ziK=0,0≤i≤N,z0j=zNj=0,0≤j<K,
where ξij=vij+κzij, 0<κ<1, and
(54)Rij≡(vij+1-vijτ-∂v∂t(xi,tj))+12σ2xi2[2hi+hi+1(vi+1j-vijhi+1-vij-vi-1jhi)-∂2v∂x2(xi,tj)]+(r-λζ)xi(vi+1j-vi-1jhi+hi+1-∂v∂x(xi,tj))+λ[∑n=0N-1‍v(xn,tj+1)·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍v(xn+1,tj+1)·∫xnxn+1‍s-xnhn+1·f(s/xi)xids-∫0Smax‍v(s,tj)f(s/xi)xids].Let(55)LεN,Kzij≡-zij+1-zijτ-σ2xi2hi+hi+1(zi+1j-zijhi+1-zij-zi-1jhi)-(r-λζ)xizi+1j-zi-1jhi+hi+1+(r+λ)zij+Cε(ξij+ε-qi)2zij-λ(∑n=0N-1‍znj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍zn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids).The operatorLεN,K satisfies the following discrete maximum principle. Hence, the difference scheme is stable and oscillation-free for arbitrary volatility and arbitrary interest rate.Lemma 2 (discrete maximum principle).
The operatorLεN,K defined by (55) on the piecewise uniform mesh ΩN×K satisfies a discrete maximum principle; that is, if uij and wij are mesh functions that satisfy u0j≥w0j, uNj≥wNj (0≤j<K), uiK≥wiK(0≤i≤N), and LεN,Kuij≥LεN,Kwij(1≤i<N,0≤j<K), then uij≥wij for all i,j.Proof.
The discrete operatorLεN,K can be written as follows:
(56)LεN,Kzij=-zij+1-zijτj+1-aizi-1j+[bi+r+λ+Cε(ξij+ε-qi)2]zij-cizi+1j-λ(∑n=0N-1‍znj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍zn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids),
where ai, bi, and ci are defined in Section 4. From (31)  we have
(57)ai>0for1<i≤N-1,ci>0for1≤i<N-1.
Clearly,
(58)bi+r+λ+Cε(ξij+ε-qi)2>0for1≤i≤N-1,b1+c1>0,ai+bi+ci>0,2≤i≤N-2,aN-1+bN-1>0.
Furthermore, we have
(59)λ-λ(∑n=0N-1‍∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍∫xnxn+1‍s-xnhn+1·f(s/xi)xids)=λ(1-∫0xN‍f(s/xi)xids)=λ(1-∫0xN/xi‍f(y)dy)>0.
Hence, we verify that the matrix associated withLεN,K is an M-matrix; see for example [32]. Thus, by the same argument as [33, Lemma 3.1] the result follows.Now we can get the following error estimates.Theorem 3.
Letv be the solution of (12)–(14) and V the solution of the finite difference scheme (23)–(25). Then one has the following error estimates,
(60)|Vij-v(xi,tj)|≤C1(τ+N-2),0≤i≤N,0≤j≤K,
where C1 is a positive constant independent of τ and N.Proof.
Under the regularity assumption onv, the linear interpolation error for any s∈[xi,xi+1] is
(61)v(s,t)-(xi+1-shi+1v(xi,t)+s-xihi+1v(xi+1,t))=12∂2v∂x2(η,t)(s-xi)(s-xi+1)
for some η∈[xi,xi+1]. Thus, we get
(62)|v(s,t)-(xi+1-shi+1v(xi,t)+s-xihi+1v(xi+1,t))|≤C2N-2,
where C2 is a positive constant independent of τ and N. Using this, we have
(63)|∑n=0N-1‍v(xn,tj+1)·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍v(xn+1,tj+1)·∫xnxn+1‍s-xnhn+1·f(s/xi)xids-∫0Smax‍v(s,tj)f(s/xi)xids|≤|∑n=0N-1‍v(xn,tj+1)·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍v(xn+1,tj+1)·∫xnxn+1‍s-xnhn+1·f(s/xi)xids-∫0Smax‍v(s,tj+1)f(s/xi)xids|+∫0Smax‍|v(s,tj+1)-v(s,tj)|f(s/xi)xids≤C2N-2∫0Smax‍f(s/xi)xids+C3τ∫0Smax‍f(s/xi)xids≤C4(τ+N-2)∫0Smax/xi‍f(y)dy≤C4(τ+N-2),
where C3 and C4 are positive constants independent of τ and N. Hence, we use the Taylor expansion to obtain
(64)|Rij|≤C5∫tjtj+1‍|∂2v∂t2(xi,t)|dt+C6h∫xi-1xi+1‍[xi2|∂4v∂x4(x,tj)|+xi|∂3v∂x3(x,tj)|]dx+C4(τ+N-2)≤C6(τ+N-2),
for 0<i<N, 0≤j<K, where C5 and C6 are also positive constants independent of τ and N. Therefore, using the barrier function Wij=C1(τ+N-2) (with the constant C1 sufficiently large), Lemma 2 implies that
(65)|Vij-v(xi,tj)|≤C1(τ+N-2),0≤i≤N,0≤j≤K,
which completes the proof.


## 6. Numerical Experiments

In this section, we verify experimentally the theoretical results obtained in the preceding section. Errors and convergence rates for the finite difference scheme are presented for two test problems.Test 1.
American put option under Kou’s jump-diffusion model with parameters:(66)σ=0.15,r=0.06,T=1,E=25,Smax=100,λ=0.1,α1=2.0465,α2=3.0775,p=0.3445.Test 2.
American put option under Kou’s jump-diffusion model with parameters:(67)σ=0.4,r=0.06,T=1,E=25,Smax=100,λ=0.1,α1=2.0465,α2=3.0775,p=0.3445.To solve the nonlinear problem (23)–(25), we use Newton iterative method. The initial guesses for Tests 1 and 2 are taken as [Vij](0)=Vij+1 for time mesh point j and the stoping criterion is
(68)maxi|[Vij](m)-[Vij](m-1)|≤10-5.For Tests1 and 2 we choose ε=0.0001,C=rE+λSmax. The computed option value V and the constraint V-V* with N=128 and K=64 for Test 1 are depicted in Figures 1 and 2, respectively. The computed option value V and the constraint V-V* with N=128 and K=64 for Test 2 are depicted in Figures 3 and 4, respectively.Figure 1
Computed option valueV for Test 1.Figure 2
The constraintV-V* for Test 1.Figure 3
Computed option valueV for Test 2.Figure 4
The constraintV-V* for Test 2.The exact solutions of our test problems are not available. We use the approximated solution ofN=1024 and K=4096 as the exact solution. We present the error estimates for different N and K at t=0. Because we only know “the exact solution” on mesh points, we use the linear interpolation to get solutions at other points. In this paper V-(x,t) denotes “the exact solution” which is a linear interpolation of the approximated solution V1024,4096. We measure the accuracy in the discrete maximum norm
(69)eN,K=maxi|Vi,0N,K-V-(xi,0)|,
and the convergence rate
(70)RN,K=log2(eN,Ke2N,4K).
The error estimates and the convergence rates in our computed solutions of Tests 1 and 2 are listed in Tables 1 and 2, respectively.Table 1
Numerical results for Test1.
N
K
error
rate
64
16
2.8859
e
-
2
—
128
64
8.3456
e
-
3
1.780
256
256
2.1187
e
-
3
1.978
512
1024
5.5979
e
-
4
1.920Table 2
Numerical results for Test2.
N
K
error
rate
64
16
4.8818
e
-
2
—
128
64
1.3125
e
-
2
1.895
256
256
3.3288
e
-
3
1.979
512
1024
7.2756
e
-
4
2.194From the figures, it is seen that the numerical solutions by our method are nonoscillatory. From Tables1 and 2, we see that eN,K/e2N,4K is close to 4, which supports the convergence estimate of Theorem 3. The numerical results of Zhang and Wang [17, 18] verify that their schemes are also second-order convergent. However, their penalty term is nonsmooth and a smoothing technique is needed for solving the nonlinear discretization system.


## 7. Discussion

### 7.1. The Crank-Nicolson Scheme

For the time discretization, we can use the Crank-Nicolson scheme to improve the accuracy of the scheme. Then the implicit finite difference scheme on the above piecewise-uniform mesh for the integrodifferential equation (12)–(14) is
(71)L-N,KVij-1/2CεVij+ε-qi-1/2CεVij+1+ε-qi=0,1≤i<N,0≤j<K,ViK=V*(xi),0≤i≤N,V0j=E,VNj=0,0≤j<K,
where
(72)L-N,KVij≡-Vij+1-Vijτ-σ2xi2hi+hi+1[12(Vi+1j-Vijhi+1+Vi+1j+1-Vij+1hi+1)-12(Vij-Vi-1jhi+Vij+1-Vi-1j+1hi)]-12(r-λζ)xi(Vi+1j-Vi-1jhi+hi+1+Vi+1j+1-Vi-1j+1hi+hi+1)+(r+λ)Vij+Vij+12-λAij+Aij+12,Aij=∑n=0N-1‍Vnj·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍Vn+1j·∫xnxn+1‍s-xnhn+1·f(s/xi)xids.
The implicit scheme would lead to full matrices due to the integral term, which makes many methods computationally expensive. Furthermore, in order to satisfy the positivity constraint and to avoid spurious oscillations in the Crank-Nicolson method [34–36], the time mesh size should satisfy the following constraint condition [24, 35]
(73)τ≤2hi2σ2Smax2+rSmaxhi+rhi2+(2rE/ε)hi2.


### 7.2. Power Penalty Methods

Some papers [12, 22, 37] have used power penalty methods to the linear complementarity problem arising from pricing American options. Wang et al. [37] prove that the solution to their penalized equation converges to that of the variational inequality problem with an arbitrary order. The power penalty methods [12, 22, 37] can also be applied to valuate American options under Kou's jump-diffusion model. The discretization method can be same as that of Section 3. Therefore, our difference scheme for the power penalty equation is
(74)L^N,KVij+λ[max(V*(xi)-Vij,0)]1/s=0,1≤i<N,0≤j<K,ViK=V*(xi),0≤i≤N,V0j=E,VNj=0,0≤j<K,
where
(75)L^N,KVij≡-Vij+1-Vijτ-σ2xi2hi+hi+1(Vi+1j-Vijhi+1-Vij-Vi-1jhi)-(r-λζ)xiVi+1j-Vi-1jhi+hi+1+(r+λ)Vij-λAij+1,Aij+1=∑n=0N-1‍Vnj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍Vn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids,λ and s are the penalty parameters. Since the power penalty term is nonsmooth, a smoothing technique is needed for solving the nonlinear discretization system.


## 7.1. The Crank-Nicolson Scheme

For the time discretization, we can use the Crank-Nicolson scheme to improve the accuracy of the scheme. Then the implicit finite difference scheme on the above piecewise-uniform mesh for the integrodifferential equation (12)–(14) is
(71)L-N,KVij-1/2CεVij+ε-qi-1/2CεVij+1+ε-qi=0,1≤i<N,0≤j<K,ViK=V*(xi),0≤i≤N,V0j=E,VNj=0,0≤j<K,
where
(72)L-N,KVij≡-Vij+1-Vijτ-σ2xi2hi+hi+1[12(Vi+1j-Vijhi+1+Vi+1j+1-Vij+1hi+1)-12(Vij-Vi-1jhi+Vij+1-Vi-1j+1hi)]-12(r-λζ)xi(Vi+1j-Vi-1jhi+hi+1+Vi+1j+1-Vi-1j+1hi+hi+1)+(r+λ)Vij+Vij+12-λAij+Aij+12,Aij=∑n=0N-1‍Vnj·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍Vn+1j·∫xnxn+1‍s-xnhn+1·f(s/xi)xids.
The implicit scheme would lead to full matrices due to the integral term, which makes many methods computationally expensive. Furthermore, in order to satisfy the positivity constraint and to avoid spurious oscillations in the Crank-Nicolson method [34–36], the time mesh size should satisfy the following constraint condition [24, 35]
(73)τ≤2hi2σ2Smax2+rSmaxhi+rhi2+(2rE/ε)hi2.


## 7.2. Power Penalty Methods

Some papers [12, 22, 37] have used power penalty methods to the linear complementarity problem arising from pricing American options. Wang et al. [37] prove that the solution to their penalized equation converges to that of the variational inequality problem with an arbitrary order. The power penalty methods [12, 22, 37] can also be applied to valuate American options under Kou's jump-diffusion model. The discretization method can be same as that of Section 3. Therefore, our difference scheme for the power penalty equation is
(74)L^N,KVij+λ[max(V*(xi)-Vij,0)]1/s=0,1≤i<N,0≤j<K,ViK=V*(xi),0≤i≤N,V0j=E,VNj=0,0≤j<K,
where
(75)L^N,KVij≡-Vij+1-Vijτ-σ2xi2hi+hi+1(Vi+1j-Vijhi+1-Vij-Vi-1jhi)-(r-λζ)xiVi+1j-Vi-1jhi+hi+1+(r+λ)Vij-λAij+1,Aij+1=∑n=0N-1‍Vnj+1·∫xnxn+1‍xn+1-shn+1·f(s/xi)xids+∑n=0N-1‍Vn+1j+1·∫xnxn+1‍s-xnhn+1·f(s/xi)xids,λ and s are the penalty parameters. Since the power penalty term is nonsmooth, a smoothing technique is needed for solving the nonlinear discretization system.



---

*Source: 651573-2013-02-07.xml*