A numerical scheme for pricing American options with transaction costs under a jump diffusion process

In this paper we develop a numerical method for a nonlinear partial integro-differential complementarity problem arising from pricing American options with transaction costs when the underlying assets follow a jump diffusion process. We first approximate the complementarity problem by a nonlinear partial integro-differential equation (PIDE) using a penalty approach. The PIDE is then discretized by a combination of a spatial upwind finite differencing and a fully implicit time stepping scheme. We prove that the coefficient matrix of the system from this scheme is an M -matrix and that the approximate solution converges to the viscosity solution to the PIDE by showing that the scheme is consistent, monotone, and unconditionally stable. We also propose a Newton's iterative method coupled with a Fast Fourier Transform for the computation of the discretized integral term for solving the fully discretized system. Numerical results will be presented to demonstrate the convergence rates and usefulness of this method.


1.
Introduction. In a complete market without transaction costs, Black-Scholes [6] proposed a partial differential equation pricing model for European options, based on the assumptions that the underlying stock price follows a geometric Brownian motion with a log-normal diffusion and the volatility is constant. However, there is evidence that these assumptions are not consistent with the market price movements in practice, which is often called the volatility skew or smile [20]. A number of models have been proposed to remedy this problem. These improved models can be categorized into three classes: stochastic volatility model [16,21], the deterministic volatility function model [12], and the jump diffusion model [10,27,42]. Among them, the jump diffusion model is very popular.
It has been shown that the price of a European option under a jump diffusion process without transaction costs satisfies a linear partial integro-differential equation (PIDE) and the price of its American counterpart satisfies a linear complementarity problem involving the partial integro-differential operator. Since exact solutions to such problems can hardly be found, numerical approximations are usually sought in practice in order to price such options. This is challenging as the PIDE and differential complementarity problems involve non-local integral terms. In [11,42], the authors treat the integral term in the conventional jump-diffusion model explicitly and the remaining terms implicitly in their numerical approaches. However, the method is only conditionally stable. An operator splitting method coupled with a Fast Fourier Transform (FFT) for the evaluation of the integral term has been used in [2]. The method is unconditionally stable and 2nd-order accurate. There are various other methods for the linear PIDE such as [14,39,40].
In [28,4] the authors showed that, in the presence of transaction costs, the price of a European option whose underlying asset price satisfies a jump diffusion process is governed by a nonlinear PIDE. Likewise, the price of an American option of this type is governed by a nonlinear complementarity problem. Although some theoretical properties of the non-linear model are established in [28,4], to our best knowledge, no numerical methods for nonlinear partial integro-differential complementarity problem can be found in the open literature. In this work, we propose a numerical method for pricing American options under transaction costs and jump diffusion based on a penalty approach and some finite difference schemes for the time and spatial derivatives as well as the integral term. We show that the system matrix of the fully discretized equation is an M -matrix and establish a convergence theory for the numerical method. To evaluate the discretized integral term, we use the FFT technique. A Newton iterative method coupled with a regular operator splitting is used to avoid the inversion of the dense matrix arising from the non-local term. The rest of this work is organized as follows.
In the next section, we will give a brief account of the American option pricing model. In Section 3, we will present some discretization schemes for the PIDE. A convergence analysis is given in Section 4 in which we show that the numerical scheme is consistent, unconditionally stable and monotone. In Section 5, we propose a Newton iterative algorithm for solving the nonlinear algebraic system and show it is convergent. In Section 6, we present some numerical experimental results to demonstrate the rates of convergence and usefulness of the numerical method for solving a number of model problems.
2. The continuous model. Consider a European option with strike price K and expiry time T on an asset whose price S follows the following stochastic differential where dZ is an increment of the standard Gauss-Wiener process and dq is the independent Poisson process with deterministic jump intensity λ. Here, σ 0 is the volatility, ν is the drift rate, η − 1 is an impulse function producing a jump from S to Sη, and κ = E(η − 1) with E the expectation operator. Following [37], it is easy to show using (2.1) that the value V (S, t) of the European option at any t ∈ [0, T ) and S > 0 when transactions do not incur any costs satisfies the following linear PIDE: where r is the risk-free interest rate, τ = T − t, and g(η) is the probability density function of the jump amplitude η, satisfying ∞ 0 g(η) dη = 1. Note that (2.2) is defined on an infinite domain. In practice, we are only interested in the computational domain 0 < S < S max , where S max is a positive number much greater than the strike price of the option. Also, various choices of g are available, but for brevity, we only consider, in this work, the following lognormal density function in Merton's model: where µ and σ J determine the mean and variance of the jumps. In this case, κ = E(η − 1) = exp(µ + σ 2 J /2) − 1. When transactions on selling and buying the asset involve costs, the value of a European option V (S, τ ) is governed by the following nonlinear PIDE: where σ is the modified volatility as a function of τ, S, V S and V SS . Assuming that the transaction cost is proportional to the value of the transaction and the portfolio is revised every ∆t units of time, where ∆t denotes a small, but non-infinitesimal, fixed time interval, Mocioalca [28] showed that the option price is the solution to (2.4) with the modified volatility models for, respectively, long and short positions in the option, where ρ is the transaction cost parameter, measured as a fraction of the volume of transactions. It is shown in [22] that these can be written as where Le is the Leland number given by Clearly, Le has to be such that σ 2 (V SS ) > 0, that is 0 ≤ Le < 1.
In [4], the author assumes that the transaction costs behave like a non-increasing positive linear function, h(x) = a − bx, where a, b ≥ 0 are two transaction cost parameters. Using this cost function, the author derived the following modified volatility We assume that k 1 , k 2 and V SS satisfy the following inequalities Under these conditions, One can easily see from (2.6) that σ 2 (S, V SS ) > 0. Note that V SS ≥ 0 is usually satisfied by European and American vanilla options as used in [4] and thus the 2nd condition in (2.7) is usually not needed in practice unless V SS is negative. If we define then the American option pricing problem can be stated as the following nonlinear partial integro-differential complementarity problem: where V * is a given non-negative function (usually the payoff function) defining a lower bound for the value of the option. Note that (2.9) contains (2.4) as the special case when V * = 0.
To solve (2.4) or (2.9), boundary and payoff conditions need to be defined. There are various types of boundary and initial conditions and payoff functions depending on the type of contingent claims. For European vanilla call and put options, they are given by Penalty methods have long been used in solving constrained optimization problems. In recent year, penalty approaches have been developed and used for solving nonlinear complementarity problems arising in financial engineering [34,41,7,8,24]. We apply the linear penalty method to (2.9) in which the complementarity problem is approximated by the following nonlinear PIDE: satisfying the initial and boundary conditions (2.10)-(2.12), where ϑ > 1 is the penalty constant and [z] + = max(z, 0) for any z and L is defined in (2.8).
The unique solvability and convergence of the solution of (2.16) to that of (2.9) when ϑ → +∞ in both infinite and finite dimensions have been discussed extensively in the open literature. In [34,43], the authors show the solution to a problem similar to (2.16) is unique and V ϑ converges to the solution V of (2.9) at the rate of order O(ϑ 1/2 ) in a Sobolev norm. In [32,17,18,30,33,31], the authors have developed convergence theories of the penalty approach for various finite-dimensional complementarity problems similar to the discretized form of (2.16) and the rates of convergence established in these works are of the order O(ϑ). Extensive numerical results have also been reported in the aforementioned works to demonstrate the theoretical convergence rates of V ϑ to V as ϑ → +∞. For brevity, we will omit these discussion, but refer the reader to these papers.
We comment that this penalty method, which can be regarded as a special case of the power penalty methods proposed in the references mentioned in the above paragraph, is different from the regularization technique appearing in the variational adjoint method such as that in [36]. The former is a technique for solving complementarity problems while the latter is usually for inverse problems.
3. Discretization of the nonlinear PIDE. We divide this discussion into two parts: the discretization of the integral term and that of the other terms in (2.16). For notational simplicity, we omit the subscript ϑ.

Discretization of the integral.
We follow the idea in [39] to approximate the integral in (2.16) or (2.4) numerically. First, we use a logarithmic transformation to transform the integral into a correlation integral so that it can be computed by the FFT algorithm. For clarity, we will suppress the time variable τ for most of the discussion in this subsection.
Let x = ln S and y = ln η. Then the integral in (2.4) can be written as . . and ∆y = ∆x. Note that the density function defined in (2.3) satisfies g(η) → 0 exponentially as |η| → ∞. So does f (y). Therefore, we approximate the correlation integral in (3.1) at x i by the following finite sumQ for all feasible i, where J 0 is a positive even integer,V k denotes a nodal approximation ofV (k∆x, ·) for any k and Eq.(3.2) defines nodal approximations toQ(x, ·) at the mesh nodes of the transformed region using nodal values ofV . However, as can be seen in the next subsection, we will discretize (2.16) in the original solution domain (0, S max ). Therefore, it is necessary to define nodal approximations to V (S) using those ofV (x) defined in (3.2) which requires interpolatingV using nodal values of V and Q using nodal values ofQ. Given a set of nodes {S 1 , S 2 , . . . , S n , . . .} on the S-axis satisfying 0 = S 1 < S 2 < . . . < S n < . . ., let V i be a nodal approximation to V (S i ) for i = 1, 2, .... We use the following steps to approximate Q(S i ).

Approximation of {V j } by the linear interpolation of {V i }.
For any integer j, let p(j) be the index such that S p(j) ≤ e xj ≤ S p(j)+1 . Then, we define the following approximation toV j From the definition, it is clear that 0 ≤ φ p(j) ≤ 1 and 0 ≤ ψ q(i) ≤ 1. Using (3.2), (3.3) and (3.4), it is easy to verify that Q i is defined by The operator Π j i (V ) is linear in V . This approximation will be used in the spatial discretization of the problem as given below.

Full discretization.
To discretize (2.16) with the nonlinear volatility (2.5) or (2.6) and boundary and initial conditions (2.10)-(2.12), we first define a mesh for (0, S max ) × (0, T ). For a given positive even integer M , let (0, S max ) be divided uniformly into M sub-intervals with mesh nodes .., M + 1 and n = 1, ..., N + 1, we define the following finite difference operators on the mesh defined above: Using these operators and Q i defined in (3.5) for a sufficiently large even integer J, we approximate (2.4) by the following finite difference equation: for i = 2, ..., M and n = 1, ..., N , wherē and V n i denotes an approximation to V (S i , τ n ) to be determined for any feasible index pair (i, n). The discretization of V S is based on the so-called upwind finite difference scheme which has been used for solving various problems in classic and financial engineering [25]. For a detailed discussion of the above discretization scheme and its convergence for the nonlinear Black-Scholes equation without any jump, we refer the reader to [23]. We comment that there are various other spatial discretization schemes for solving Black-Scholes equations such as [3,9,13,15,19,35] and any of these methods can be used for the discretization of the differential operators in (2.16). Using (2.10)-(2.12), we define the initial and boundary conditions for the above system as follows The system (3.7) can be rewritten as where J is a large even positive number as used in (3.5), (3.14) The nonlinear finite-difference system (3.11) along with (3.8)-(3.10) form a recursive system that determines the unknown vectorV n+1 = (V n+1 2 , . . . , V n+1 M ) for n = 1, ..., N . From (3.5) and (3.6) we see that the sum on the LHS of (3.11) is linear in Then, (3.11) can be written in a matrix form as follows where I denotes the identity matrix,V * = (V * 2 , ..., V * M ) and A n+1 and B n+1 are given by Clearly, B n+1 contains contributions from the boundary conditions in (3.9)-(3.10).
For the system matrix of (3.16) we have the following results.
since 1 τn > 0 and the other terms are nonnegative. Therefore, (3.18) is also satisfied. Also, it is obvious that A n is irreducible and thus it is irreducibly diagonally dominant with positive diagonal and non-positive off-diagonal entries. By [29], A n is an M -matrix for any given V n . Now we investigate d ij defined in by (3.15). From Since f is a probability density function, we have Using these inequalities we can conclude that all the off-diagonal entries of A n (V n )− λD are non-positive and This implies Therefore, all the diagonal entries of A n (V n ) − λD are positive, off-diagonal entries are non-positive, and the matrix is strictly diagonally dominant. Hence, it is an M -matrix by [29].
Using Theorem 3.1 and the fact that the initial and boundary conditions in (3.8)-(3.10) are all non-negative, it is trivial to show the following corollary. 4. Convergence of the numerical scheme. In [4], the author shows the existence and uniqueness of the viscosity solution to (2.4). In this section we will prove that the solution to (3.16) converges to the viscosity solution to (2.4) as the mesh sizes approach zero. It has been shown in [5] that the convergence of the fully discretized system (3.16) to the viscosity solution of a full nonlinear 2nd-order PDE is guaranteed if the discretization is consistent, stable and monotone. Thus, in the rest of this section we will prove the convergence of our numerical scheme by showing that it satisfies these properties. For i = 2, 3, ..., M and n = 1, ..., N , introduce a functional H n+1 i defined by where

DONNY CITRA LESMANA AND SONG WANG
Then, it is easy to see that (3.11) (or (3.7)) becomes for all feasible i and n. For this scheme, we have the following lemma.
Proof. To prove this lemma, we show that, for any ε > 0 and i = 2, 3, . . . , M, Let us consider the first nonlinear term on the RHS of (4.1). We discuss this term in the following two cases corresponding to the choices of σ defined in (2.5) and (2.6).
Combining the monotonicity of σ 2 (S, z)z and the properties of the linear terms on the RHS of (4.1) and (4.4) we have Similarly, using the monotonicity of σ 2 (S, z)z and (4.5) and noticing that [ i ] + for any ε ≥ 0, it is easy to show that (4.3) also holds true. Hence, the discretization scheme is monotone.
The stability of the method is given in the following lemma.
where ∆τ = max n ∆τ n , g 1 , g 2 and g 3 are the initial and boundary conditions defined in (3.8)-(3.10) and · ∞ denotes either the discrete or continuous L ∞ -norm.

DONNY CITRA LESMANA AND SONG WANG
Using (3.21) we have from this inequality It is easy to verify that the above system can be written as where p n = ϑ∆τ n or 0. Repeatedly using this relationship, we obtain where p = ϑ∆τ . In the above we used the fact that z/(1 + z) is increasing for z ≥ 0.
In this case, from (3.8), (3.9) and (3.10) it is easy to see that Combining the above two cases we have (4.6).
Note that for American Vanilla put options, V * = g 1 as defined in (2.13). The consistency of the numerical scheme is given in the following lemma: The proof is standard since both of the time and spatial discretization schemes are standard and have been used extensively in the literature for 2nd-order partial differential equations. Therefore, we omit the proof of this lemma. Combining the above three lemmas we have the following convergence result. Proof. In [4] the authors show that if a discretization scheme for a fully nonlinear 2nd order PDE is monotone, stable and consistent, then the solution to the fully discretized system converges to the viscosity solution to the PDE. Therefore, this theorem is just a consequence of Lemmas 4.1, 4.2 and 4.3.

5.
Solution of the nonlinear system. 5.1. The European case. In this case ϑ = 0 in (3.11). Note that D arising from discretization of correlation product (integral) term is a dense matrix. Therefore, the solution of (3.16) is computationally expensive because of the inversion of D.
To solve this, we use a Newton iterative method coupled with a regular splitting technique used in [38]. First, we write (3.16) in the following form and V n+1 M +1 are defined in (3.9) and (3.10). The Jacobian matrix of j for all feasible i and j and D is as defined before.
For the Jacobian E n+1 (V n+1 ) − λD, we have the following results.
Proof. For simplicity of notation, we omit the superscript n+1. To show that E is an M -matrix, from [29] we see that it suffices to prove that | for at least one index i. For volatility given by (2.5), it is easy to check that E n+1 (V n+1 ) = A n+1 (V n+1 ). Thus, from Theorem 3.1, both E and E − λD are M -matrices. Now we will prove that E and E − λD are M -matrices for the volatility given in (2.6).
Let us first consider J i,i−1 . Using the definition of α n+1 i in (3.12) and (2.6) we have Similarly it can be shown that From these expressions we see that for any i = 2, 3, ..., M with the convention that E 2,1 = 0 = J M,M +1 . Therefore, E is an M -matrix by [29]. Now, let us consider E − λD. From the definition of D, it is clear that all the off-diagonal entries of E − λD are non-positive and This implies Therefore, all the diagonal entries of E n+1 (V n+1 ) − λD are positive and the matrix is strictly diagonally dominant. Hence, it is an M -matrix.
Using the Jacobian of H n+1 , we propose the following Newton algorithm for (3.16):
Note that D is a dense matrix. Therefore, the solution of (5.1) is computationally expensive because of the inversion of D. To solve this, we use a regular splitting technique. Let Now, we introduce this following definition From Theorem 5.1 we have E n+1 (V n+1 ) is an M -matrix for any given V n+1 and hence P −1 ≥ 0. Also from (3.11) and (3.16) we have D > 0 and thus R ≥ 0. Therefore, we have a regular splitting. Using this splitting, we define an iterative scheme for (5.1) as follows. 3) The following lemma establishes the convergence of the iterative method (5.3).
The product R(δW ) k is computed in the following two FFT operations. Let F F T (u) denote the FFT of u and define the vector which generates the row vector of C by permutation. First, we compute F F T (F ) and F F T ((δU ) k ). Then, we compute the inverse FFT of the product of F F T (F ) and F F T ((δU ) k ).
The numerical implementation for the system (3.16) can be summarized in this following algorithm.

Algorithm N2
1. Choose tolerances ε 1 , ε 2 > 0. Let n = 1 and evaluate the discrete initial conditionV 1 = (V 1 2 , ..., V 1 M ) using (3.8). 6. Compute F F T ((δU ) k ) and the inverse FFT on the product of F F T (F ) and [1,(δW ) k+1 ∞ ] < ε 1 then stop. Otherwise, set k := k + 1, go to Step 6. 9. Set (δW ) l = (δW ) k+1 . If (δW ) l ∞ < ε 2 then stop. Otherwise, compute V l+1 = V l + (δW ) l and set l = l + 1, go to Step 4. 10. SetV n+1 =V l+1 . If n < N , let n := n + 1 and go to Step 3. Otherwise, stop. 5.2. The American case. From (3.7) we see that the discretized nonlinear system from American options differs from that from the European options by a penalty term −ϑ[V * − V ] + . Since this term is non-smooth in V n+1 for each n = 1, ..., N , a smoothing technique such as those proposed in [34,17] is usually used to locally smooth out the penalty term near zero. It has been shown in [34,17] that the smoothed penalty term is an increasing function of V n+1 and thus its derivative is non-negative. Therefore, all the analysis presented for European options above hold true in this case.

Compute
6. Numerical results. In this section we demonstrate the rates of convergence and numerical performance of our method using European call and put options as well as an American put option.   [30,50] in Figure 6.2 for the European call option and put, respectively. From these plots we see that the prices of the options increase as the transaction cost parameter b increases, as expected in practice.
We now investigate the rate of convergence numerically, which requires an exact or reference solution. Since the exact solution to this problem is unknown, we use the numerical solution to (3.16) on the uniform mesh with M = 5121 (h = 1 5120 ) and N = 2560 (∆τ = 1 2560 ) as an "exact" or reference solution. To determine these rates, we choose a sequence of meshes generated by successively halving the mesh sizes of the previous ones, starting from a given coarse mesh. Using the reference solution we then calculate the following ratios of the numerical solutions from two consecutive meshes:  where V β α denotes the computed solution on the mesh with spatial mesh size α and time mesh size β and || · || h,2 denotes the discrete L 2 -norm defined by The computed errors in || · || h,2 and the ratios two consecutive errors are listed in Tables 6.1 and 6.2 for the European call and put respectively. From the tables we see that the rates of convergence of our method is about 2, showing that our numerical method is 2nd order accurate in the L 2 -norm.
Test Problem 2. American put option with the system parameters: r = 0.1, σ 0 = 0.2, K = 40, T = 1, S max = 80, µ = 0.01, σ J = 0.5, λ = 0.1, and ϑ = 1000. The payoff and boundary conditions are given in (2.13)-(2.15). To see the difference between the European and American options, we plot the computed values of both of the options at t = 0 (or τ = T ) for S ∈ [20, 80] in Figure  6.3. From the figure it is clear that the American put option is more valuable than the European put option as expected.
Furthermore, we investigate the effect of the transaction cost parameters on the value of the option. To see this, we plot the values of the option at t = 0 (or τ = T ) for three different values of b on the interval [30,50] in Figure 6.4 in which the curve for a = 0 and b = 0 is the price of the American put option without transaction cost. From this figure we see that the price of the option increases as the transaction parameter b increases as expected in practice.    7. Conclusion. In this paper we have proposed and analyzed an upwind finite difference scheme for a nonlinear Partial Integro-Differential Equation (PIDE) and a nonlinear partial integro-differential complementarity problem arising from European and American option pricing with transaction costs under a jump diffusion process. The complementarity problem is first approximated by a nonlinear PIDE with a penalty term which contains the nonlinear PIDE for pricing European options as a special case. A finite difference scheme has then been developed for the penalized PIDE and the convergence of the solution to discretized system to the viscosity solution of the penalty equation has been proven. To solve the discretized system with a dense coefficient matrix, we have developed a fast iterative method coupled with a regular splitting technique. To speed up the computation of the discretized integral term, we have proposed to use the FFT algorithm. Numerical experiments have been performed to confirm the theoretical results and demonstrate the usefulness of the method. The numerical results show that the prices of both European and American options are increasing functions of the transaction cost parameters.