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

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.

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 at 𝑥=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 𝑥=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 𝑀-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 Section 3. 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

Let 𝑣 denote the value of an American put option with strike price 𝐸 on the underlying asset 𝑥 and time 𝑡. It is known that the price 𝑣 under a jump-diffusion model satisfies the following partial integrodifferential complementarity problem [14–17]: ()[),𝑣()𝐿𝑣x,𝑡≥0,𝑥>0,𝑡∈0,𝑇𝑥,𝑡−𝑉^∗(𝑥)[],()⋅𝑣()≥0,𝑥>0,𝑡∈0,𝑇𝐿𝑣𝑥,𝑡𝑥,𝑡−𝑉^∗(𝑥[),𝑣())=0,𝑥>0,𝑡∈0,𝑇𝑥,𝑇=𝑉^∗(𝑥)𝑣()[],𝑣()[],,𝑥≥0,𝑡=𝑇,0,𝑡=𝐸,𝑥=0,𝑡∈0,𝑇𝑥,𝑡→0,𝑥→+∞,𝑡∈0,𝑇 (1)(2)(3)(4)(5)(6) where 𝐿 denotes the partial integrodifferential operator defined by 𝐿𝑣(𝑥,𝑡)≡−𝜕𝑣−1𝜕𝑡2𝜎²𝑥²𝜕²𝑣𝜕𝑥²−(𝑟−𝜆𝜁)𝑥𝜕𝑣+()𝜕𝑥𝑟+𝜆𝑣−𝜆^∞₀𝑣𝑥𝑦,𝑡𝑓𝑦d𝑦, (7)𝑉^∗(𝑥) is the final (payoff) condition defined by 𝑉^∗(𝑥)=max{𝐸−𝑥,0}, (8)𝜎 is the volatility of the underlying asset, 𝑟 is the risk free interest rate, 𝑡 is the current time, 𝑇 is the maturity date, and 𝑓(𝑦) is the probability function of the jump amplitude 𝑦 with the obvious properties that for all 𝑦, 𝑓(𝑦)≥0 and ∫^∞₀𝑓(𝑦)d𝑦=1, the constant 𝜁 is given by 𝜁=∫^∞₀(𝑦−1)𝑓(𝑦)d𝑦. In Kou's model, 𝑓(𝑦) is the following log-double-exponential density 𝑓𝑦=𝑞𝛼₂𝑦^𝛼2−1,𝑦<1,𝑝𝛼₁𝑦^−𝛼1−1,𝑦≥1, (9) where 𝑝,𝑞,𝛼₁>1 and 𝛼₂ are positive constants such that 𝑝+𝑞=1. It can be shown that, in this case, 𝜁=𝑝𝛼₁/(𝛼₁−1)+𝑞𝛼₂/(𝛼₂+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 𝑟−𝜆𝜁≠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, 𝐿𝑣(𝑥,𝑡)−𝐶𝜀𝑣()(𝑥)[),𝑣()𝑥,𝑡+𝜀−𝑞=0,𝑥>0,𝑡∈0,𝑇𝑥,𝑇=𝑉^∗(𝑥)𝑣()[],𝑣()[],,𝑥≥0,𝑡=𝑇,0,𝑡=𝐸,𝑥=0,𝑡∈0,𝑇𝑥,𝑡→0,𝑥→+∞,𝑡∈0,𝑇 (10) where 𝐶≥𝑟𝐸+𝜆𝑆_max is a positive constant and 𝑞(𝑥)=𝐸−𝑥.

By adding a penalty term 𝐶𝜀𝑣()(𝑥),𝑥,𝑡+𝜀−𝑞 (11) 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 𝑣(𝑥,𝑡)≫𝑞(𝑥), and hence the partial integrodifferential equation is approximately satisfied. When 𝑣(𝑥,𝑡) approaches 𝑞(𝑥), this term is approximately equal to 𝐶 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,𝑆_max). 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 𝑆_max=4𝐸. The boundary condition at 𝑥=𝑆_max is chosen to be 𝑣(𝑆_max,𝑡)=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: 𝐿𝑣(𝑥,𝑡)−𝐶𝜀𝑣()(𝑥)()∈𝑥,𝑡+𝜀−𝑞=0,𝑥,𝑡0,𝑆_max×(),𝑣()0,𝑇𝑥,𝑇=𝑉^∗(𝑥),𝑥∈0,𝑆_max,𝑣(0,𝑡)=𝐸,𝑣𝑆_max,𝑡=0,𝑡∈[0,𝑇]. (12)(13)(14)

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 Ω^𝑁 on the space interval [0,𝑆_max]: 𝑥_𝑖=ℎ𝜎ℎ,𝑖=1,1+²(𝑟−𝜆𝜁𝑖−1),𝑖=2,…,𝑁, (15) where 𝑆ℎ=^max1+𝜎²/().𝑟−𝜆𝜁𝑁−1 (16) For the time discretization, we use a uniform mesh Ω^𝐾 on [0,𝑇] with 𝐾 mesh elements. Then the piecewise uniform mesh Ω^𝑁×𝐾 on Ω=(0,𝑆_max)×(0,𝑇) is defined to be the tensor product Ω^𝑁×𝐾=Ω^𝑁×Ω^𝐾. It is easy to see that the mesh sizes ℎ_𝑖=𝑥_𝑖−𝑥_𝑖−1 and 𝜏_𝑗=𝑡_𝑗−𝑡_𝑗−1 satisfy ℎ_𝑖=𝜎ℎ,𝑖=1,₂𝑟−𝜆𝜁ℎ,𝑖=2,…,𝑁,𝜏=𝜏_𝑗=𝑇𝐾,𝑗=1,…,𝐾, (17)(18) respectively.

The space derivatives of (12) are approximated with central differences on the above piecewise-uniform mesh: 𝜕𝑣𝑥𝜕𝑥^𝑖≈𝑣,𝑡^𝑖+1(𝑡)−𝑣_𝑖−1(𝑡)ℎ_𝑖+ℎ_𝑖+1,𝜕²𝑣𝜕𝑥²𝑥_𝑖≈2,𝑡ℎ_𝑖+ℎ_𝑖+1𝑣_𝑖+1(𝑡)−𝑣_𝑖(𝑡)ℎ_𝑖+1−𝑣_𝑖(𝑡)−𝑣_𝑖−1(𝑡)ℎ_𝑖. (19)

The integral term 𝐼=^∞₀𝑣𝑥𝑦,𝑡𝑓𝑦d𝑦 (20) of (12) can be approximated by a fast method as in [14–16]. By making the change of variable 𝑦=𝑠/𝑥, we have 𝐼=^∞₀𝑣𝑥𝑦,𝑡𝑓𝑦d𝑦=^∞₀𝑣(𝑠,𝑡)𝑓(𝑠/𝑥)𝑥d𝑠. (21) By using the linear interpolation, we can obtain an approximation 𝐼_𝑖≈𝐴_𝑖(𝑡)=^𝑁−1_𝑛=0^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1𝑣𝑥_𝑛+,𝑡𝑠−𝑥_𝑛ℎ_𝑛+1𝑣𝑥_𝑛+1𝑓,𝑡𝑠/𝑥_𝑖𝑥_𝑖d𝑠 (22) of 𝐼 at each mesh point 𝑥_𝑖 for 𝑖=1,2,…,𝑁−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 𝑀-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 𝐿^𝑁,𝐾𝑉^𝑗_𝑖−𝐶𝜀𝑉^𝑗_𝑖+𝜀−𝑞_𝑖𝑉=0,1≤𝑖<𝑁,0≤𝑗<𝐾,_𝐾𝑖=𝑉^∗𝑥_𝑖𝑉,0≤𝑖≤𝑁,_𝑗0=𝐸,𝑉^𝑗_𝑁=0,0≤𝑗<𝐾, (23)(24)(25) where 𝐿^𝑁,𝐾𝑉^𝑗_𝑖𝑉≡−^𝑖𝑗+1−𝑉^𝑗_𝑖𝜏−𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖ℎ_𝑖+1−𝑉^𝑗_𝑖−𝑉^𝑗_𝑖−1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1+(𝑟+𝜆)𝑉^𝑗_𝑖−𝜆𝐴_𝑖^𝑗+1,𝐴_𝑖^𝑗+1=^𝑁−1_𝑛=0𝑉_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑉^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠. (26) Then from the above solution 𝑉^𝑗_𝑖, we can obtain the optimal stopping price which is the maximum asset price such that 𝑉^𝑗_𝑖=𝑉^∗_𝑖 for each 𝑡_𝑗.

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 𝑣 under a generalized VG process satisfies the following partial integrodifferential complementarity problem [30, 31]: ()[),𝑣()𝐿𝑣𝑥,𝑡≥0,𝑥>0,𝑡∈0,𝑇𝑥,𝑡−𝑉^∗(𝑥)[],()⋅𝑣()≥0,𝑥>0,𝑡∈0,𝑇𝐿𝑣𝑥,𝑡𝑥,𝑡−𝑉^∗(𝑥[),𝑣())=0,𝑥>0,𝑡∈0,𝑇𝑥,𝑇=𝑉^∗(𝑥)𝑣()[],𝑣()[],,𝑥≥0,𝑡=𝑇,0,𝑡=𝐸,𝑥=0,𝑡∈0,𝑇𝑥,𝑡→0,𝑥→+∞,𝑡∈0,𝑇 (27) where 𝐿 denotes the partial integrodifferential operator defined by 𝐿𝑣(𝑥,𝑡)≡−𝜕𝑣−1𝜕𝑡2𝜎²𝑥²𝜕²𝑣𝜕𝑥²−(𝑟+𝑤)𝑥𝜕𝑣𝜕𝑥+𝑟𝑣−_ℝ𝑣𝑥𝑒^𝑦,𝑡−𝑣(𝑥,𝑡)𝑔𝑦d𝑦, (28)𝑤 is some “compensation constant” given by 𝑤=_ℝ1−𝑒^𝑦𝑔𝑦d𝑦, (29)𝑔(𝑦) 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 𝑀-matrix.

4. Positivity Constraint

In this section, we will prove that our scheme satisfies the early exercise constraint.

Let 𝑎_𝑖=𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1,𝑏_𝑖=𝜎²𝑥²_𝑖ℎ_𝑖ℎ_𝑖+1,𝑐_𝑖=𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖+1+(𝑟−𝜆𝜁)𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1,𝑖=1,…,𝑁−1. (30) We can obtain 𝑎_𝑖>𝜎²𝑥₁𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1=𝜎²𝑥₁−(𝑟−𝜆𝜁)ℎ_𝑖ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖𝑥_𝑖=𝜎²ℎ−(𝑟−𝜆𝜁)𝜎²/𝑟−𝜆𝜁ℎℎ_𝑖+ℎ_𝑖+1ℎ_𝑖𝑥_𝑖𝑐≥0,1<𝑖≤𝑁−1,_𝑖>𝜎²𝑥₁𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖+1+(𝑟−𝜆𝜁)𝑥_𝑖ℎ_𝑖+ℎ_𝑖+1=𝜎²𝑥₁+(𝑟−𝜆𝜁)ℎ_𝑖+1ℎ_𝑖+ℎ_𝑖+1ℎ_𝑖+1𝑥_𝑖=𝜎²ℎ+(𝑟−𝜆𝜁)𝜎²/𝑟−𝜆𝜁ℎℎ_𝑖+ℎ_𝑖+1ℎ_𝑖+1𝑥_𝑖≥0,1≤𝑖<𝑁−1. (31)

Theorem 1. The approximate option values generated by the implicit-explicit difference scheme (23)–(25) satisfy 𝑉^𝑗_𝑖≥𝑉^∗𝑥_𝑖,𝑖,𝑗∈Ω^𝑁×𝐾, (32) provided that 𝐶≥𝑟𝐸+𝜆𝑆_max.

Proof. We apply the similar technique of [24] to give the proof.
The schemes (23) can be written as 1𝜏+𝑏_𝑖+𝑟+𝜆𝑉^𝑗_𝑖=1𝜏𝑉_𝑖^𝑗+1+𝑎_𝑖𝑉^𝑗_𝑖−1+𝑐_𝑖𝑉^𝑗_𝑖+1+𝜆𝐴_𝑖^𝑗+1+𝐶𝜀𝑉^𝑗_𝑖+𝜀−𝑞_𝑖. (33) The difference 𝑢^𝑗_𝑖=𝑉^𝑗_𝑖−𝑞_𝑖 (34) satisfies the following equation: 1𝜏+𝑏_𝑖+𝑟+𝜆𝑢^𝑗_𝑖=1𝜏𝑢_𝑖^𝑗+1+𝑎_𝑖𝑢^𝑗_𝑖−1+𝑐_𝑖𝑢^𝑗_𝑖+1+𝐶𝜀𝑢^𝑗_𝑖()𝑥+𝜀+𝜆𝜁+1_𝑖−(𝑟+𝜆)𝐸+𝜆^𝑁−1_𝑛=0𝑢_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑢^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠+𝜆^𝑁−1_𝑛=0𝑞_𝑛⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓s/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑞_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠. (35) Next, by defining 𝑢^𝑗=min_𝑖𝑢^𝑗_𝑖, (36) let 𝑘 be an index such that 𝑢^𝑗_𝑘=𝑢^𝑗. (37) For 𝑖=𝑘, it follows from (35) that 1𝜏+𝑏_𝑘+𝑟+𝜆𝑢^𝑗≥1𝜏𝑢_𝑘^𝑗+1+𝑎_𝑘+𝑐_𝑘𝑢^𝑗+𝐶𝜀𝑢^𝑗()𝑥+𝜀+𝜆𝜁+1_𝑘−(𝑟+𝜆)𝐸+𝜆^𝑁−1_𝑛=0𝑢_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=0}𝑢^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘d𝑠+𝜆^𝑁−1_𝑛=0𝑞_𝑛⋅^𝑥n+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=0}𝑞_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘d𝑠, (38) where we have used (31). Since ^𝑁−1_𝑛=0𝑞_𝑛⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=0}𝑞_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘>d𝑠𝐸−𝑆_max⋅^𝑥𝑁₀𝑓𝑠/𝑥_𝑘𝑥_𝑘=d𝑠𝐸−𝑆_max⋅^{𝑥𝑁/𝑥𝑘}₀𝑓𝑦d𝑦≥𝐸−𝑆_max,𝜁+1=^∞₀𝑓𝑦=𝑦−1d𝑦+1_∞0𝑦𝑓𝑦d𝑦>0, (39) we conclude that 1𝜏+𝑟+𝜆𝑢^𝑗−𝐶𝜀𝑢^𝑗++𝜀𝑟𝐸+𝜆𝑆_max≥1𝜏𝑢_𝑘^𝑗+1+𝜆^𝑁−1_𝑛=1𝑢_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=1}𝑢^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘≥1d𝑠𝜏𝑢^𝑗+1+𝜆𝑢^{𝑗+1𝑁−1}_𝑛=1^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=1𝑥𝑛+1𝑥𝑛}𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘d𝑠. (40) If we assume that 𝑢^𝑗+1≥0, (41) from (40) we can get 𝐹𝑢^𝑗≥0, (42) where 𝐹(𝑢)=1𝜏+𝑟+𝜆𝑢−𝐶𝜀+𝑢+𝜀𝑟𝐸+𝜆𝑆_max. (43) Obviously, 𝐹^′(𝑢)=1𝜏+𝑟+𝜆+𝐶𝜀(𝑢+𝜀)²𝐹(0)≥0,=−𝐶+𝑟𝐸+𝜆𝑆_max≤0, (44) where we have used the assumption 𝐶≥(𝑟𝐸+𝜆𝑆_max). Hence, from (42) we can obtain 𝑢^𝑗≥0. (45) Consequently, by induction on 𝑗, it follows from (34) that 𝑉^𝑗_𝑖≥𝑞_𝑖,0≤𝑖≤𝑁, (46) for 𝑗=𝐾,…,1,0.
Next we prove that 𝑉^𝑗_𝑖≥0. As above, we define 𝑉^𝑗=min_𝑖𝑉^𝑗_𝑖, (47) and let 𝑘 be an index such that 𝑉^𝑗_𝑘=𝑉^𝑗. (48) For 𝑖=𝑘, from (33) we obtain 1𝜏+𝑏_𝑘+𝑟+𝜆𝑉^𝑗_𝑘≥1𝜏𝑉^𝑗+1+𝑎_𝑘+𝑐_𝑘𝑉^𝑗+𝜆𝐴_𝑘^𝑗+1+𝐶𝜀𝑉^𝑗_𝑘+𝜀−𝑞_𝑘. (49) Hence, we have 1𝜏+𝑟+𝜆𝑉^𝑗_𝑘≥1𝜏𝑉^𝑗+1+𝜆𝐴_𝑘^𝑗+1+𝐶𝜀𝑉^𝑗_𝑘+𝜀−𝑞_𝑘. (50) Since we have proved that 𝑉^𝑗_𝑖≥𝑞_𝑖, we get 1𝜏+𝑟+𝜆𝑉^𝑗_𝑘≥1𝜏𝑉^𝑗+1+𝜆^𝑁−1_𝑛=0𝑉_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘+d𝑠_{𝑁−1𝑛=0}𝑉^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑘𝑥_𝑘d𝑠. (51) By the above inequality. Hence, by induction on 𝑗, we can obtain 𝑉^𝑗_𝑖≥0,𝑖,𝑗∈Ω^𝑁×𝐾, (52) for 𝑗=𝐾,…,1,0.
Combine (46) with (52) to complete the proof.

5. Error Estimates

To investigate the convergence of the method, note that the error functions 𝑧^𝑗_𝑖=𝑉^𝑗_𝑖−𝑣^𝑗_𝑖 (0≤𝑖≤𝑁,0≤𝑗≤𝐾) are the solutions of the discrete problem −𝑧_𝑖^𝑗+1−𝑧^𝑗_𝑖𝜏−𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1𝑧^𝑗_𝑖+1−𝑧^𝑗_𝑖ℎ_𝑖+1−𝑧^𝑗_𝑖−𝑧^𝑗_𝑖−1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖𝑧^𝑗_𝑖+1−𝑧^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1+(𝑟+𝜆)𝑧^𝑗_𝑖+𝐶𝜀𝜉^𝑗_𝑖+𝜀−𝑞_𝑖²𝑧^𝑗_𝑖−𝜆^𝑁−1_𝑛=0𝑧_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑧^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠=𝑅^𝑗_𝑖𝑧,1≤𝑖<𝑁,0≤𝑗<𝐾,_𝐾𝑖𝑧=0,0≤𝑖≤𝑁,_𝑗0=𝑧^𝑗_𝑁=0,0≤𝑗<𝐾, (53) where 𝜉^𝑗_𝑖=𝑣^𝑗_𝑖+𝜅𝑧^𝑗_𝑖, 0<𝜅<1, and 𝑅^𝑗_𝑖≡𝑣_𝑖^𝑗+1−𝑣^𝑗_𝑖𝜏−𝜕𝑣𝑥𝜕𝑡^𝑖,𝑡_𝑗+12𝜎²𝑥²_𝑖2ℎ_𝑖+ℎ_𝑖+1𝑣^𝑗_𝑖+1−𝑣^𝑗_𝑖ℎ_𝑖+1−𝑣^𝑗_𝑖−𝑣^𝑗_𝑖−1ℎ_𝑖−𝜕²𝑣𝜕𝑥²𝑥_𝑖,𝑡_𝑗+(𝑟−𝜆𝜁)𝑥_𝑖𝑣^𝑗_𝑖+1−𝑣^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1−𝜕𝑣𝑥𝜕𝑥^𝑖,𝑡_𝑗+𝜆^𝑁−1_𝑛=0𝑣𝑥_𝑛,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑣𝑥_𝑛+1,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠−^𝑆max₀𝑣𝑠,𝑡_𝑗𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠. (54)

Let 𝐿_𝜀^𝑁,𝐾𝑧^𝑗_𝑖𝑧≡−^𝑖𝑗+1−𝑧^𝑗_𝑖𝜏−𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1𝑧^𝑗_𝑖+1−𝑧^𝑗_𝑖ℎ_𝑖+1−𝑧^𝑗_𝑖−𝑧^𝑗_𝑖−1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖𝑧^𝑗_𝑖+1−𝑧^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1+(𝑟+𝜆)𝑧^𝑗_𝑖+𝐶𝜀𝜉^𝑗_𝑖+𝜀−𝑞_𝑖²𝑧^𝑗_𝑖−𝜆^𝑁−1_𝑛=0𝑧_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑧^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠. (55)

The operator 𝐿_𝜀^𝑁,𝐾 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 operator 𝐿_𝜀^𝑁,𝐾 defined by (55) on the piecewise uniform mesh Ω^𝑁×𝐾 satisfies a discrete maximum principle; that is, if 𝑢^𝑗_𝑖 and 𝑤^𝑗_𝑖 are mesh functions that satisfy 𝑢^𝑗₀≥𝑤^𝑗₀, 𝑢^𝑗_𝑁≥𝑤^𝑗_𝑁 (0≤𝑗<𝐾), 𝑢^𝐾_𝑖≥𝑤^𝐾_𝑖(0≤𝑖≤𝑁), and 𝐿_𝜀^𝑁,𝐾𝑢^𝑗_𝑖≥𝐿_𝜀^𝑁,𝐾𝑤^𝑗_𝑖(1≤𝑖<𝑁,0≤𝑗<𝐾), then 𝑢^𝑗_𝑖≥𝑤^𝑗_𝑖 for all 𝑖,𝑗.

Proof. The discrete operator 𝐿_𝜀^𝑁,𝐾 can be written as follows: 𝐿_𝜀^𝑁,𝐾𝑧^𝑗_𝑖𝑧=−^𝑖𝑗+1−𝑧^𝑗_𝑖𝜏_𝑗+1−𝑎_𝑖𝑧^𝑗_𝑖−1+𝑏_𝑖+𝑟+𝜆+𝐶𝜀𝜉^𝑗_𝑖+𝜀−𝑞_𝑖²𝑧^𝑗_𝑖−𝑐_𝑖𝑧^𝑗_𝑖+1−𝜆^𝑁−1_𝑛=0𝑧_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑧^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠, (56) where 𝑎_𝑖, 𝑏_𝑖, and 𝑐_𝑖 are defined in Section 4. From (31) we have 𝑎_𝑖𝑐>0for1<𝑖≤𝑁−1,_𝑖>0for1≤𝑖<𝑁−1. (57) Clearly, 𝑏_𝑖+𝑟+𝜆+𝐶𝜀𝜉^𝑗_𝑖+𝜀−𝑞_𝑖²𝑏>0for1≤𝑖≤𝑁−1,₁+𝑐₁𝑎>0,_𝑖+𝑏_𝑖+𝑐_𝑖𝑎>0,2≤𝑖≤𝑁−2,_𝑁−1+𝑏_𝑁−1>0. (58) Furthermore, we have 𝜆−𝜆^𝑁−1_𝑛=0^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0𝑥𝑛+1𝑥𝑛}𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠=𝜆1−^𝑥𝑁₀𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠=𝜆1−^{𝑥𝑁/𝑥𝑖}₀𝑓𝑦d𝑦>0. (59)
Hence, we verify that the matrix associated with 𝐿_𝜀^𝑁,𝐾 is an 𝑀-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. Let 𝑣 be the solution of (12)–(14) and 𝑉 the solution of the finite difference scheme (23)–(25). Then one has the following error estimates, 𝑉^𝑗_𝑖−𝑣𝑥_𝑖,𝑡_𝑗≤𝐶₁𝜏+𝑁⁻²,0≤𝑖≤𝑁,0≤𝑗≤𝐾, (60) where 𝐶₁ is a positive constant independent of 𝜏 and 𝑁.

Proof. Under the regularity assumption on 𝑣, the linear interpolation error for any 𝑠∈[𝑥_𝑖,𝑥_𝑖+1] is 𝑣()−𝑥𝑠,𝑡^𝑖+1−𝑠ℎ_𝑖+1𝑣𝑥_𝑖,𝑡+𝑠−𝑥_𝑖ℎ_𝑖+1𝑣𝑥_𝑖+1=1,𝑡2𝜕²𝑣𝜕𝑥²𝜂,𝑡𝑠−𝑥_𝑖𝑠−𝑥_𝑖+1 (61) for some 𝜂∈[𝑥_𝑖,𝑥_𝑖+1]. Thus, we get 𝑣()−𝑥𝑠,𝑡^𝑖+1−𝑠ℎ_𝑖+1𝑣𝑥_𝑖,𝑡+𝑠−𝑥_𝑖ℎ_𝑖+1𝑣𝑥_𝑖+1,𝑡≤𝐶₂𝑁⁻², (62) where 𝐶₂ is a positive constant independent of 𝜏 and 𝑁. Using this, we have ^𝑁−1_𝑛=0𝑣𝑥_𝑛,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑣𝑥_𝑛+1,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠−^𝑆max₀𝑣𝑠,𝑡_𝑗𝑓𝑠/𝑥_𝑖𝑥_𝑖≤d𝑠_{𝑁−1𝑛=0}𝑣𝑥_𝑛,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑣𝑥_𝑛+1,𝑡_𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠−^𝑆max₀𝑣𝑠,𝑡_𝑗+1𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_𝑆max0𝑣𝑠,𝑡_𝑗+1−𝑣𝑠,𝑡_𝑗𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠≤𝐶₂𝑁^−2𝑆max₀𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠+𝐶₃𝜏^𝑆max₀𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠≤𝐶₄𝜏+𝑁^{−2𝑆max/𝑥𝑖}₀𝑓𝑦d𝑦≤𝐶₄𝜏+𝑁⁻², (63) where 𝐶₃ and 𝐶₄ are positive constants independent of 𝜏 and 𝑁. Hence, we use the Taylor expansion to obtain 𝑅^𝑗_𝑖≤𝐶₅^𝑡𝑗+1_𝑡𝑗𝜕²𝑣𝜕𝑡²𝑥_𝑖,𝑡d𝑡+𝐶₆ℎ^𝑥𝑖+1_𝑥𝑖−1𝑥²_𝑖𝜕⁴𝑣𝜕𝑥⁴𝑥,𝑡_𝑗+𝑥_𝑖𝜕³𝑣𝜕𝑥³𝑥,𝑡_𝑗d𝑥+𝐶₄𝜏+𝑁⁻²≤𝐶₆𝜏+𝑁⁻², (64) for 0<𝑖<𝑁, 0≤𝑗<𝐾, where 𝐶₅ and 𝐶₆ are also positive constants independent of 𝜏 and 𝑁. Therefore, using the barrier function 𝑊^𝑗_𝑖=𝐶₁(𝜏+𝑁⁻²) (with the constant 𝐶₁ sufficiently large), Lemma 2 implies that 𝑉^𝑗_𝑖−𝑣𝑥_𝑖,𝑡_𝑗≤𝐶₁𝜏+𝑁⁻²,0≤𝑖≤𝑁,0≤𝑗≤𝐾, (65) 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: 𝜎=0.15,𝑟=0.06,𝑇=1,𝐸=25,𝑆_max𝛼=100,𝜆=0.1,₁=2.0465,𝛼₂=3.0775,𝑝=0.3445. (66)

Test 2. American put option under Kou’s jump-diffusion model with parameters: 𝜎=0.4,𝑟=0.06,𝑇=1,𝐸=25,𝑆_max𝛼=100,𝜆=0.1,₁=2.0465,𝛼₂=3.0775,𝑝=0.3445. (67)

To solve the nonlinear problem (23)–(25), we use Newton iterative method. The initial guesses for Tests 1 and 2 are taken as [𝑉^𝑗_𝑖]⁽⁰⁾=𝑉_𝑖^𝑗+1 for time mesh point 𝑗 and the stoping criterion is max_𝑖𝑉^𝑗_𝑖^(𝑚)−𝑉^𝑗_𝑖^(𝑚−1)≤10⁻⁵. (68)

For Tests 1 and 2 we choose 𝜀=0.0001,C=𝑟𝐸+𝜆S_max. The computed option value 𝑉 and the constraint 𝑉−𝑉^∗ with 𝑁=128 and 𝐾=64 for Test 1 are depicted in Figures 1 and 2, respectively. The computed option value 𝑉 and the constraint 𝑉−𝑉^∗ with 𝑁=128 and 𝐾=64 for Test 2 are depicted in Figures 3 and 4, respectively.

Figure 1 

Computed option value 𝑉 for Test 1.

Figure 2 

The constraint 𝑉−𝑉∗ for Test 1.

Figure 3 

Computed option value 𝑉 for Test 2.

Figure 4 

The constraint 𝑉−𝑉∗ for Test 2.

The exact solutions of our test problems are not available. We use the approximated solution of 𝑁=1024 and 𝐾=4096 as the exact solution. We present the error estimates for different 𝑁 and 𝐾 at 𝑡=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 𝑉(𝑥,𝑡) denotes “the exact solution” which is a linear interpolation of the approximated solution 𝑉^1024,4096. We measure the accuracy in the discrete maximum norm 𝑒^𝑁,𝐾=max_𝑖𝑉^𝑁,𝐾_𝑖,0−𝑉𝑥_𝑖,0, (69) and the convergence rate 𝑅^𝑁,𝐾=log₂𝑒^𝑁,𝐾𝑒^2𝑁,4𝐾. (70) 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 Test 1.

Table 2 

Numerical results for Test 2.

From the figures, it is seen that the numerical solutions by our method are nonoscillatory. From Tables 1 and 2, we see that 𝑒^𝑁,𝐾/𝑒^2𝑁,4𝐾 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 𝐿^𝑁,𝐾𝑉^𝑗_𝑖−1/2𝐶𝜀𝑉^𝑗_𝑖+𝜀−𝑞_𝑖−1/2𝐶𝜀𝑉_𝑖^𝑗+1+𝜀−𝑞_𝑖𝑉=0,1≤𝑖<𝑁,0≤𝑗<𝐾,_𝐾𝑖=𝑉^∗𝑥_𝑖𝑉,0≤𝑖≤𝑁,_𝑗0=𝐸,𝑉^𝑗_𝑁=0,0≤𝑗<𝐾, (71) where 𝐿^𝑁,𝐾𝑉^𝑗_𝑖𝑉≡−^𝑖𝑗+1−𝑉^𝑗_𝑖𝜏−𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+112𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖ℎ_𝑖+1+𝑉^𝑗+1_𝑖+1−𝑉_𝑖^𝑗+1ℎ_𝑖+1−12𝑉^𝑗_𝑖−𝑉^𝑗_𝑖−1ℎ_𝑖+𝑉_𝑖^𝑗+1−𝑉^𝑗+1_𝑖−1ℎ_𝑖−12(𝑟−𝜆𝜁)𝑥_𝑖𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1+𝑉^𝑗+1_𝑖+1−𝑉^𝑗+1_𝑖−1ℎ_𝑖+ℎ_𝑖+1+()𝑉𝑟+𝜆^𝑗𝑖+𝑉_𝑖^𝑗+12𝐴−𝜆^𝑗𝑖+𝐴_𝑖^𝑗+12,𝐴^𝑗_𝑖=^𝑁−1_𝑛=0𝑉^𝑗_𝑛⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑉^𝑗_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠. (72) 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] 𝜏≤2ℎ²_𝑖𝜎²𝑆²_max+𝑟𝑆_maxℎ_𝑖+𝑟ℎ²_𝑖+(2𝑟𝐸/𝜀)ℎ²_𝑖. (73)

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 𝐿^𝑁,𝐾𝑉^𝑗_𝑖+𝜆max𝑉^∗𝑥_𝑖−𝑉^𝑗_𝑖,0^1/𝑠𝑉=0,1≤𝑖<𝑁,0≤𝑗<𝐾,_𝐾𝑖=𝑉^∗𝑥_𝑖𝑉,0≤𝑖≤𝑁,_𝑗0=𝐸,𝑉^𝑗_𝑁=0,0≤𝑗<𝐾, (74) where 𝐿^𝑁,𝐾𝑉^𝑗_𝑖𝑉≡−^𝑖𝑗+1−V^𝑗_𝑖𝜏−𝜎²𝑥²_𝑖ℎ_𝑖+ℎ_𝑖+1𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖ℎ_𝑖+1−𝑉^𝑗_𝑖−𝑉^𝑗_𝑖−1ℎ_𝑖−(𝑟−𝜆𝜁)𝑥_𝑖𝑉^𝑗_𝑖+1−𝑉^𝑗_𝑖−1ℎ_𝑖+ℎ_𝑖+1+(𝑟+𝜆)𝑉^𝑗_𝑖−𝜆𝐴_𝑖^𝑗+1,𝐴_𝑖^𝑗+1=^𝑁−1_𝑛=0𝑉_𝑛^𝑗+1⋅^𝑥𝑛+1_𝑥𝑛𝑥_𝑛+1−𝑠ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖+d𝑠_{𝑁−1𝑛=0}𝑉^𝑗+1_𝑛+1⋅^𝑥𝑛+1_𝑥𝑛𝑠−𝑥_𝑛ℎ_𝑛+1⋅𝑓𝑠/𝑥_𝑖𝑥_𝑖d𝑠, (75)𝜆 and 𝑠 are the penalty parameters. Since the power penalty term is nonsmooth, a smoothing technique is needed for solving the nonlinear discretization system.

Acknowledgments

The authors would like to thank the anonymous referees for their helpful comments and suggestions for the improvement of this paper. The work was supported by the Zhejiang Province Natural Science Foundation (Grant nos. Y6100021, Y6110310), Research Project (Grant no. Y201016720) of the Zhejiang Provincial Education Department and the Ningbo Municipal Natural Science Foundation (Grant nos. 2012A610035, 2012A610036).