#### 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 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 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 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]: where denotes the partial integrodifferential operator defined by is the final (payoff) condition defined by 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 , and , the constant is given by . In Kou's model, is the following log-double-exponential density where and are positive constants such that . It can be shown that, in this case, . When the jump rate is zero, the partial integrodifferential operator reduces to the standard Black-Scholes operator [1]. In this paper, we assume that .

The above linear complementarity problem (1)–(6) can be solved by a penalty approach. Let be a small regularization parameter and consider the following initial-boundary value problem, where is a positive constant and .

By adding a penalty term 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 into . 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 . The boundary condition at is chosen to be . 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:

#### 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 : where For the time discretization, we use a uniform mesh on with mesh elements. Then the piecewise uniform mesh on is defined to be the tensor product . It is easy to see that the mesh sizes and satisfy respectively.

The space derivatives of (12) are approximated with central differences on the above piecewise-uniform mesh:

The integral term of (12) can be approximated by a fast method as in [14–16]. By making the change of variable , we have By using the linear interpolation, we can obtain an approximation of at each mesh point for .

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 where 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]:
where denotes the partial integrodifferential operator defined by
is some “compensation constant” given by
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 We can obtain

Theorem 1. *The approximate option values generated by the implicit-explicit difference scheme (23)–(25) satisfy
**
provided that .*

*Proof. *We apply the similar technique of [24] to give the proof.

The schemes (23) can be written as
The difference
satisfies the following equation:
Next, by defining
let be an index such that
For , it follows from (35) that
where we have used (31). Since
we conclude that
If we assume that
from (40) we can get
where
Obviously,
where we have used the assumption . Hence, from (42) we can obtain
Consequently, by induction on , it follows from (34) that
for .

Next we prove that . As above, we define
and let be an index such that
For , from (33) we obtain
Hence, we have
Since we have proved that , we get
By the above inequality. Hence, by induction on , we can obtain
for .

Combine (46) with (52) to complete the proof.

#### 5. Error Estimates

To investigate the convergence of the method, note that the error functions () are the solutions of the discrete problem where , , and

Let

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 , (), , and , then for all .*

* Proof. *The discrete operator can be written as follows:
where , , and are defined in Section 4. From (31) we have
Clearly,
Furthermore, we have

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,
**
where is a positive constant independent of and .*

*Proof. *Under the regularity assumption on , the linear interpolation error for any is
for some . Thus, we get
where is a positive constant independent of and . Using this, we have
where and are positive constants independent of and . Hence, we use the Taylor expansion to obtain
for , , where and are also positive constants independent of and . Therefore, using the barrier function (with the constant sufficiently large), Lemma 2 implies that
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:

*Test 2. * American put option under Kou’s jump-diffusion model with parameters:

To solve the nonlinear problem (23)–(25), we use Newton iterative method. The initial guesses for Tests 1 and 2 are taken as for time mesh point and the stoping criterion is

For Tests 1 and 2 we choose . The computed option value and the constraint with and for Test 1 are depicted in Figures 1 and 2, respectively. The computed option value and the constraint with and for Test 2 are depicted in Figures 3 and 4, respectively.

The exact solutions of our test problems are not available. We use the approximated solution of and as the exact solution. We present the error estimates for different and at . 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 . We measure the accuracy in the discrete maximum norm and the convergence rate The error estimates and the convergence rates in our computed solutions of Tests 1 and 2 are listed in Tables 1 and 2, respectively.

From the figures, it is seen that the numerical solutions by our method are nonoscillatory. From Tables 1 and 2, we see that 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 where 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]

##### 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 where 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).