Double Discretization Difference Schemes for Partial Integrodifferential Option Pricing Jump Diffusion Models

and Applied Analysis 3 For the sake of clarity in the presentation we recall that in a jump-diffusion model, the modified stochastic differential equation SDE for the underlying asset is dS S μdt σdz ( η − 1dq, 1.1 where S is the underlying stock price, μ is the drift rate, σ is the volatility, dz is the increment of Gauss-Wiener process, and dq is the Poisson process. The random variable representing the jump amplitude is denoted by η, and the expected relative jump size is denoted by K E η − 1 . The jump intensity of the Poisson process is denoted by λ. Based on the SDE 1.1 the resulting PIDE for a contingent claim V S, t is given by 7, 14, 29 : ∂V ∂t 1 2 σ2S2 ∂2V ∂S2 r − λK S ∂S − r λ V


Introduction
Since empirical studies revealed that the normality of the log returns, as assumed by Black and Scholes, could not capture features like heavy tails and asymmetries observed in market-data log-returns densities 1 , a number of models try to explain these empirical observations: stochastic volatility 2, 3 , deterministic local volatility 4, 5 , jump diffusion 6, 7 , and infinite activity Lévy models 8-11 .The two last types of models, discussed in 12 and 13, chapters 14, 15 allow to calibrate the model to market price of options and reproduce a wide variety of implied volatility skews/smiles.These models are characterized by partial integrodifferential equations PIDEs that involve a second-order differential operator and a nonlocal integral term that requires specific treatment and presents additional difficulties.
In order to solve the PIDE problem numerically, Andersen and Andreasen 14 use an unconditionally stable ADI finite difference method and accelerate it using fast Fourier transform FFT .In 15-17 wavelet methods are applied to infinite activity jump-diffusion models.Interesting analytic-numerical treatments for Lévy models have been introduced in 18-20 .The so-called COS method for pricing European options is presented in 18 .This is based on the knowledge of the characteristic function of the jump operator and the close relation of the characteristic function with the series coefficients of the Fourier-cosine expansion of the density function.In 19 , an expansion of the characteristic function of local volatility models with Lévy jumps is developed.The authors in 20 derive an analytical formula for the price of European options for any model including local volatility and Poisson jump process by using Malliavin calculus techniques.Various authors apart of 14 used finite difference schemes for PIDEs in 21-27 .Discretization of the integral term leads to full matrices due to its nonlocal character.Dealing with the integral term several challenges arise, for instance, how to approximate the integral term and how to localizate a bounded computational domain, also the selection of the boundary conditions of the numerical domain and the problem of the double discretization of the differential and integral part of the PIDE.
Tavella and Randall in 26 used an implicit time discretization and propose a stationary fairly rapid convergent iterative method to solve the full matrix problem quoted above but without a careful numerical analysis.A generalization of this iterative method to price American options is proposed in 25 .
In the outstanding paper 22 the authors propose an explicit-implicit finite difference scheme for solving parabolic PIDEs with possibly singular kernels when the random evolution of the underlying asset is driven by a time-inhomogeneous jump-diffusion process.The authors study stability and convergence of the proposed scheme as well as rates of convergence.However, they use backward or forward difference quotients of only first order depending on the sign of the coefficient of the convection term in order to avoid oscillations.An improvable issue of 22 is that in order to approximate the truncated integral term, they assume a particular behavior of the solution outside of the bounded numerical domain.
An efficient solution of PIDEs for the jump-diffusion Merton model is proposed in 24 with a very efficient treatment of the resulting dense linear system by using a circulant preconditioned conjugate gradient method.However, in 24 , they only consider the particular case where the jump sizes have zero mean, μ J 0. They also assume a particular behavior of the solution outside the bounded numerical domain.
Almendral and Oosterlee 28 present an implicit discretization of the PIDE jumpdiffusion model on an uniform grid using finite differences, where a splitting technique combined with FFT is used to accelerate the dense matrix-vector product.The authors also assume a particular behaviour of the solution outside of the bounded numerical domain, in a similar way to 24 .
In 21 a finite difference method for PIDE associated with the CGMY infinite activity Lévy model is treated.The equations are discretized in space by the collocation method and in time by an explicit backward differentiation formula.The integral part is transformed into a Volterra equation.After integration by parts and taking advantage of the vanishing derivative behaviour of the payoff function for large asset values, the authors are able to truncate properly the integral for the case of put and butterfly options.
In 27 the price of European and American options under PIDE Kou's jump-diffusion model is solved using finite differences on nonuniform grids, and time stepping is performed using the implicit Rannacher scheme.The evaluation of the integral term is efficient from the computational cost point of view, assuming that the behaviour of the solution for large values of the underlying asset follows the asymptotic behaviour.
For the sake of clarity in the presentation we recall that in a jump-diffusion model, the modified stochastic differential equation SDE for the underlying asset is where S is the underlying stock price, μ is the drift rate, σ is the volatility, dz is the increment of Gauss-Wiener process, and dq is the Poisson process.The random variable representing the jump amplitude is denoted by η, and the expected relative jump size is denoted by K E η − 1 .The jump intensity of the Poisson process is denoted by λ.Based on the SDE 1.1 the resulting PIDE for a contingent claim V S, t is given by 7, 14, 29 : ∂V ∂t where r is the risk-free interest rate, the probability density of the jump amplitude is given by g η , and f S is the payoff function.Merton's jump-diffusion model assumes that jump sizes are log-normally distributed with mean μ J and standard deviation σ J , that is, In this paper we consider Merton's jump-diffusion model for a vanilla call option with payoff function The aim of the paper is the construction and numerical analysis of an explicit finite difference numerical scheme of the PIDE 1.2 -1.3 , with a different treatment of the integral part from previously quoted authors.Instead of assuming long-term information about the solution we perform a full discretization of the integral part, involving the unknown function values in the numerical scheme, discriminating the finite truncation domain and the infinite remaining one.As a consequence, this strategy involves a double discretization with respect to the spatial variable.With respect to the time variable an explicit forward approximation is used.An account of the advantages of this explicit approach has been explained and applied in 30 .This paper is organized as follows.Section 2 deals with a transformation of variables in order to eliminate both the advection and the reaction terms of the PIDE 1.2 .Then the integral part is split in two parts: a finite integral J 1 and an infinite one J 2 , and the last is again transformed into a finite integral.The separation point A of the two split integrals is becoming a parameter that could be chosen according to the criteria used by 16, 22, 31 .A suitable choice of parameter A is the one used by other authors when they truncate the numerical domain.For instance in 27 , one takes A 4E; in Section 6 we take A 3E. Section 3 deals with the construction of the numerical scheme and the selection of the numerical domain that always involves the difficulty of the consideration of the boundary conditions.For the case of a PIDE this issue is even more relevant because the values throughout all the unbounded integral domain are unknown.The spatial numerical domain is divided in two parts by the parameter A: 0, A and A, ∞ .In the 0, A domain, the stepsize discretization is h Δx, consisting in N equidistributed mesh points X i , 1 ≤ i ≤ N. The domain A, ∞ is transformed into the 0, 1 by transformation z A/X.In the transformed domain 0, 1 we consider a stepsize discretization δ Δz and M mesh points with Mδ 1.When the interval 0, 1 is reversed to the domain A, ∞ , the reversed mesh points X i , N ≤ i ≤ N M − 1 become nonuniformly distributed.Hence, the numerical scheme for problem 2.9 -2.10 is forward in time with time-step discretization k Δτ.The approximation of ∂ 2 U/∂X 2 is centered in 0, A of the unique parameter h, and in A, A/δ the approximation of ∂ 2 U/∂X 2 involves a nonuniform stepsize discretization h j depending on δ and the value of j.The numerical approximation of the integrals is evaluated using trapezoidal quadrature rules with stepsizes h and δ, respectively.The boundary conditions at the boundary of our numerical domain are as follows.At X 0 we assume that the solution is zero according to the vanilla call option problem.At our largest finite value considered X A/δ we assume a linear behavior of the solution.This hypothesis has been previously used in 32 .
In Section 4 sufficient conditions for stability and positivity of the numerical solution are given in terms of the three parameters h ΔX, k Δτ, and δ Δz as well as of the parameter A. Consistency of the scheme is treated in Section 5. Section 6 includes illustrative examples showing the possible advantages of our new discretization approach.Finally conclusions are shown in Section 7. If Vector v is said to be nonnegative if v j ≥ 0 for all 1 ≤ j ≤ n.Then we denote v ≥ 0. For a matrix A a ij m×n in R m×n , we denote by A ∞ max 1≤i≤m { n j 1 |a ij |}.Matrix A is said to be nonnegative if a ij ≥ 0 for all 1 ≤ i ≤ m, 1 ≤ j ≤ n, and we denote A ≥ 0.

Transformation of the Integrodifferential Problem
For the sake of convenience we introduce a transformation of variables to remove both the advection and the reaction terms of the PIDE problem 1.2 -1.3 .
Let us consider the transformation In order to approximate the integral appearing in 2.2 and further discretization it is convenient the change of the variable φ Xη,

2.4
Let us denote

2.6
Following 33, page 201 let us consider the substitution into J 2 , obtaining the expression Taking into account 2.4 -2.8 , the problem 2.2 -2.3 can be written in the form

Numerical Scheme Construction
In this section a difference scheme for problem 2.9 -2.10 is constructed.With respect to the time variable, given τ with 0 < τ ≤ T , let k be the time-step discretization k Δτ τ/L, and τ n nk, 0 ≤ n ≤ L, with L integer.With respect to the spatial variable X, given an arbitrary positive fixed A > 0, we construct a uniform grid in 0, A , with the spatial step discretization h ΔX A/N, with X j jh, 0 ≤ j ≤ N, being N integer.Note that the integral J 2 X, τ, A given by 2.6 requires the evaluation of the unknown U φ, τ at points φ ∈ A, ∞ .As this integral has been transformed into 2.8 over the interval 0, 1 for the variable z, see 2.7 , we consider a uniform mesh of 0, 1 into M points, of the form z j jδ, 1 ≤ j ≤ M, where M is integer, Mδ 1, with M ≥ 3. Taking into account 2.7 , for the original variable X in A, ∞ , one has and since z j jδ, Thus the spatial domain 0,

3.4
For the approximation of ∂ 2 U/∂X 2 we consider two types of finite differences: for the internal points of 0, A , and denoting h j X j 1 − X j > 0, for the points X j lying in A, A/2δ .Note that the discrete operator Δ n j has different expressions depending on the ubication of j, see 3.5 and 3.6 .
From the previous approximations, for the internal points we have where J n 1,i and J n 2,i are approximations of the composite trapezoidal type of integrals appearing in 2.6 , 2.8 : Let us denote g i,j g X j /X i .The approximation J n 1,i takes the form where the first term for j 0 does not appear due to the null value of the limit of function g η given by 1.4 as η tends to zero.On the other hand, considering the assumption that U φ, τ n has asymptotic linear behaviour as φ → ∞ z → 0 and using 1.4 it follows that the integrand of 3.9 verifies Consequently, the last term of J n 2,i related to j N M vanishes, and one gets

3.12
The numerical scheme 3.7 -3.12 needs to incorporate the transformed initial condition and the boundary conditions for i 0: and assuming linear behaviour of the solution for large values of the spatial variable at any time, we have Δ n N M−1 0 and null integral term approximation J n

3.15
For the sake of convenience in the study of stability we introduce a vector formulation of the scheme 3.7 -3.15 .Let us consider the vector in R N M−1 as and let P ∈ R N M−1 × N M−1 be the tridiagonal matrix related to the differential part and defined by where

3.18
Let B b ij be the matrix in

3.19
From the previous notation the scheme 3.7 -3.15 can be written in the form 3.20

Positive and Stability of the Numerical Solution
Dealing with prices of contracts modeled by PIDE, the solution must be nonnegative.In this section we show that numerical solution provided by scheme 3.7 -3.15 is conditionally positive and stable.We begin with the following result.

4.1
Thus, under condition C1 , condition 4.1 holds true.With respect to the nonuniform grid, note that for N ≤ i ≤ N M − 2, from 3.18 one gets that α i > 0, γ i > 0, and α N M−1 0. From 3.18 we also have that In order to guarantee the nonnegativeness of the remaining entries of matrix P , let us introduce the function With this notation, we have that β i appearing in 3.18 satisfies Taking into account that M δ 1, with M ≥ 3, then for N 1 ≤ i ≤ N M − 2, both numerator and denominator of 4.5 are positive.Thus Thus H i is strictly increasing for N 1 ≤ i ≤ N M − 2, and its minimum is achieved at i N 1 with the value and condition 4.4 holds true under the condition From condition C2 , properties 4.2 and 4.8 are satisfied and matrix P is nonnegative.
Note that as matrix B defined by 3.19 is always nonnegative, from Lemma 4.1 and 3.20 starting from nonnegative initial vector U 0 , the following result is established.

Theorem 4.2. With the hypotheses and notation of Lemma 4.1, the solution {u
The next result will be used below to guarantee stability.
Lemma 4.3.Matrices P and B defined by 3.17 , 3.18 , and 3.19 satisfy the following.
1 Under conditions (C1) and (C2) of Lemma 4.1, Proof.By Lemma 4.1, under hypotheses C1 and C2 , all the entries of matrix P are nonnegative.Thus,

4.9
We also have Hence and from the definition of • ∞ , it follows that P ∞ 1.This proves part 1.From 3.19 , for a fixed i, with 1 ≤ i ≤ N M − 2 one gets where are the trapezoidal approximation rules with N and M points, approximating the integrals Let X q i be the maximum of g X/X i in R , given by Note that the integrand g X/X i is increasing in the interval 0, X q i and decreasing for X q i , ∞ .In order to upper bound 4.12 we consider two cases.Firstly, let us assume that 2h < X q i < A − h and let us denote j 0 to be the first integer with 1 ≤ j 0 ≤ N − 3, such that X j 0 h X j 0 1 ≤ X q i < X j 0 2 X j 0 2h.

4.16
From the properties of the lower Riemann sums, it follows that 4.17 Taking into account the values g i,j 0 1 and g i,j 0 2 located just before and after X q i , together with 4.17 , from 4.12 it follows that

4.18
In an analogous way, for the second situation where X q i ≤ 2h, one gets again 4.18 .Finally, if X q i ≥ A − h is also true that

4.19
Now we will upper bound 4.13 .Let us denote with lim z → 0 G i z 0. Let z q i be the maximum of G i z :

4.21
In this case we could distinguish the three possible situations 2δ < z q i < 1 − δ; z q i ≤ 2δ; z q i ≥ 1 − δ, and upper bounding the lower Riemann sums relative to 4.13 one gets

4.22
Abstract and Applied Analysis 13 From 4.11 , 4.18 , and 4.22 together with the fact that ∞ 0 g X/X i dX X i , one gets for each value of i, Hence, from 4.23 one gets independently of the value of the size of matrix B.
For the sake of clarity and as there are many definitions of stability in the literature we recall our concept of stability in the next definition.Definition 4.4.Let {u n i } be a numerical solution of the PIDE 2.9 , 2.10 computed from the scheme 3.7 -3.15 with stepsizes h ΔX in 0, A , δ Δz in 0, 1 , and k Δτ in 0, τ .Let {U n } be the corresponding vector form, that is, where W > 0 is independent of n, h, δ, and k.
If the property 4.25 is satisfied for appropriate relationships between the stepsizes h, δ, and k, then one says that the strong uniform stability is conditional.Theorem 4.5.With the previous notation, the numerical solution {u n i } of the scheme 3.7 -3.15 is strongly uniformly • ∞ stable if one satisfies the condition 0 < δ ≤ 1/3 together with Proof.Note that scheme 3.7 -3.15 is equivalent to the vector form scheme 3.20 .Under condition 4.26 , by Lemma 4.3 one gets, after taking norms in 3.20 , Hence, from 4.27 , Bernouilli's inequality, and 0 ≤ n ≤ L, kL τ ≤ T , one gets Thus the conditional strong uniform stability is established.
14 Abstract and Applied Analysis

Consistency
We say that a numerical difference scheme is consistent with a PIDE, if the exact theoretical solution of the PIDE approximates well to the exact solution of the difference scheme as the stepsize discretization tends to zero, 34, 35 .
Let us write the scheme 3.7 -3.12 in the form F u n i 0, where and let us write the PIDE 2.9 in the form where where J 1 and J 2 are given by 2.6 -2.8 .Let us denote U n i U X i , τ n to be the value of the theoretical solution of PIDE 5.2 .Let A > 0 such that X i < A. We denote by the following expression the local truncation error T n i U :

5.4
In order to prove the consistency, we must show that Assuming that U is twice continuously partially differentiable with respect to τ and using Taylor's expansions about X i , τ n one gets where

5.7
Let us assume that U admits four times continuous partial derivatives with respect to X, and let us denote

5.8
In accordance with 34, page 101 let us explain the local consistency error of J 1 , see 2.6 , by where U X, τ n g X/X i 2 in 5.11 means the second derivative with respect to the variable X, see 33 .
In an analogous way, let us explain the local consistency error of the unbounded integral J 2 , see 2.8 , by W n i 4 .

16
Abstract and Applied Analysis Thus which proves the consistency of the scheme with the PIDE.

Numerical Results
In the following examples the code was run on Matlab.The first example illustrates that stability conditions of Theorem 4.5 cannot be removed.
Example 6.1.Consider the vanilla call option problem 1.2 -1.5 under Merton jump diffusion model with parameters T 1, r 0.05, E 10, σ 0.1, μ J 0.5, K 0.7, and λ 0.1.Taking A 3E, δ 0.1, and h 0.3, Figure 1 shows that when the stability conditions 4.26 are satisfied results are good k 0.01 , while if the stability conditions are broken k 0.015625 results are unreliables.
The next example shows the robustness of our numerical scheme under changes of the jump intensity λ of the model.Example 6.2.Taking the same parameters of Example 6.1 apart from λ and the stepsize discretizations h 0.3, δ 0.1, and k 0.01, Figure 2 shows the variation of the solution with parameter λ, where λ 0 corresponds to the Black-Scholes case.
In the next example, the error is the difference between the numerical solution V S, 0 computed by 3.7 -3.15 and 2.1 and the exact solution given by Merton's formula 7 .Example 6.3 shows that the error of the numerical solution with fixed δ decreases with the uniform stepsize h ΔX about the strike E 10, while the error close to the truncation   The next Example 6.4 shows that the errors in the right boundary of the numerical domain when one uses finite difference schemes, quoted by 28 , can be reduced with our double spatial discretization by decreasing the stepsize δ.Example 6.4.Taking the problem of Example 6.3 with fixed h 0.3, Figure 4 shows the error reduction of the numerical solution about the right boundary of the numerical domain when parameter δ decreases, while the error about the strike remains stationary.

Conclusions
This work introduces a new discretization strategy for solving partial integrodifferential equations which involves the discretization of the unknown in the unbounded part of the integral.This fact increases the accuracy of the numerical solution in the boundary of the numerical domain as it is shown in Example 6.4.

Figure 1 :
Figure 1: Satisfying and breaking stability conditions.

Figure 4 :
Figure 4: Absolute errors with several values of δ and a fixed h.