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

Authors: M.-C. Casabán; R. Company; L. Jódar; J.-V. Romero
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

A new discretization strategy is introduced for the numerical solution of partial integrodifferential equations appearing in option pricing jump diffusion models. In order to consider the unknown behaviour of the solution in the unbounded part of the spatial domain, a double discretization is proposed. Stability, consistency, and positivity of the resulting explicit scheme are analyzed. Advantages of the method are illustrated with several examples.


---

## Body

## 1. 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 jump-diffusion 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(1.1)dSS=μdt+σdz+(η-1)dq,
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]:
(1.2)∂V∂t+12σ2S2∂2V∂S2+(r-λK)S∂V∂S-(r+λ)V+λ∫0∞‍V(Sη,t)g(η)dη=0,0<S<∞,0≤t<T,(1.3)V(S,T)=f(S),0<S<∞,
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,
(1.4)g(η)=exp(-1/2((ln(η)-μJ)/σJ)2)σJη2π.In this paper we consider Merton’s jump-diffusion model for a vanilla call option with payoff function(1.5)f(S)=max(S-E,0).
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. Section2 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 J1 and an infinite one J2, 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 Xi,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 Xi,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 ∂2U/∂X2 is centered in [0,A] of the unique parameter h, and in [A,A/δ] the approximation of ∂2U/∂X2 involves a nonuniform stepsize discretization hj 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 Section4 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.Ifv=(v1,v2,…,vn)T is a vector in ℝn, we denote its infinite norm ∥v∥∞=max{|vj|;i≤j≤n}. Vector v is said to be nonnegative if vj≥0 for all 1≤j≤n. Then we denote v≥0. For a matrix A=(aij)m×n in ℝm×n, we denote by ∥A∥∞=max1≤i≤m{∑j=1n‍|aij|}. Matrix A is said to be nonnegative if aij≥0 for all 1≤i≤m,1≤j≤n, and we denote A≥0.


## 2. 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(2.1)X=exp((r-λK)(T-t))S,τ=T-t,U(X,τ)=exp((r+λ)(T-t))V(S,t),
and note that problem (1.2)-(1.3) is transformed into the problem
(2.2)∂U∂τ=12σ2X2∂2U∂X2+λ∫0∞‍U(Xη,τ)g(η)dη,0<X<∞,0<τ≤T,(2.3)U(X,0)=f(X),0<X<∞.In order to approximate the integral appearing in (2.2) and further discretization it is convenient the change of the variable
(2.4)ϕ=Xη,∫0∞‍U(Xη,τ)g(η)dη=1X∫0∞‍U(ϕ,τ)g(ϕX)dϕ.Let us denote(2.5)J=J(X,τ)=∫0∞‍U(ϕ,τ)g(ϕX)dϕ.TakingA>0, let us decompose
(2.6)J=J1+J2,J1=∫0A‍U(ϕ,τ)g(ϕX)dϕ,J2=∫A∞‍U(ϕ,τ)g(ϕX)dϕ,J1=J1(X,τ,A),J2=J2(X,τ,A).
Following [33, page 201] let us consider the substitution
(2.7)z=Aϕ,
into J2, obtaining the expression
(2.8)J2=A∫01‍U(Az,τ)g(AXz)1z2dz.
Taking into account (2.4)–(2.8), the problem (2.2)-(2.3) can be written in the form
(2.9)∂U∂τ=σ2X22∂2U∂X2+λX(J1+J2),0<X<∞,0<τ≤T,(2.10)U(X,0)=f(X),0<X<∞.


## 3. 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 Xj=jh,0≤j≤N, being N integer.Note that the integralJ2(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 zj=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
(3.1)Xj=AzN+M-j,N≤j≤N+M-1,
and since zj=jδ,
(3.2)Xj=A1-(j-N)δ,N≤j≤N+M-1.
Thus the spatial domain [0,A]∪[A,A/δ]=[0,A/δ] is split into N+M points
(3.3)Xj={jh,0≤j≤N,A1-(j-N)δ,N≤j≤N+M-1,
from those only the first N+1 are equidistributed.Let us denote(3.4)ujn≈U(Xj,τn),0≤j≤N+M-1,0≤n≤L,∂U∂τ(Xj,τn)≈ujn+1-ujnk,1≤j≤N+M-1,0≤n≤L.
For the approximation of ∂2U/∂X2 we consider two types of finite differences:
(3.5)∂2U∂X2(Xj,τn)≈uj+1n-2ujn+uj-1nh2=Δjn,1≤j≤N-1,
for the internal points of [0,A], and denoting hj=Xj+1-Xj>0,
(3.6)∂2U∂X2(Xj,τn)≈2(uj+1nhj(hj+hj-1)+uj-1nhj-1(hj+hj-1)-ujnhjhj-1)=Δjn,N≤j≤N+M-2,
for the points Xj lying in [A,A/2δ]. Note that the discrete operator Δjn 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(3.7)uin+1=uin+k2σ2Xi2Δin+kλXi(J1,in+J2,in),1≤i≤N+M-2,
where J1,in and J2,in are approximations of the composite trapezoidal type of integrals appearing in (2.6), (2.8):
(3.8)J1,in≈∫0A‍U(ϕ,τn)g(ϕXi)dϕ,(3.9)J2,in≈∫A∞‍U(ϕ,τn)g(ϕXi)dϕ=A∫01‍U(Az,τn)g(AXiz)1z2dz.
Let us denote gi,j=g(Xj/Xi). The approximation J1,in takes the form
(3.10)J1,in=h(∑j=1N-1‍ujngi,j+12uNngi,N),1≤i≤N+M-2,
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
(3.11)U(Az,τn)g(AXiz)1z2→0,asz→0.
Consequently, the last term of J2,in related to j=N+M vanishes, and one gets
(3.12)J2,in=δA(12uNngi,NXN2+∑j=N+1N+M-1‍ujngi,jXj2),1≤i≤N+M-2.
The numerical scheme (3.7)–(3.12) needs to incorporate the transformed initial condition
(3.13)ui0=f(Xi)=max(Xi-E,0),1≤i≤N+M-1,
and the boundary conditions for i=0:
(3.14)u0n=0,0≤n≤L,
and assuming linear behaviour of the solution for large values of the spatial variable at any time, we have ΔN+M-1n=0 and null integral term approximation J1,N+M-1n=J2,N+M-1n=0. Hence, considering (3.7) for i=N+M-1, one gets
(3.15)uN+M-1n+1=uN+M-1n=uN+M-10,0≤n≤L-1.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 ℝN+M-1 as
(3.16)Un=[u1nu2n…uN+M-1n]t,
and let P∈ℝ(N+M-1)×(N+M-1) be the tridiagonal matrix related to the differential part and defined by
(3.17)P=[β1γ100⋯0α2β2γ20⋯00α3β3γ3⋯0⋱⋱⋱⋱⋱⋱0⋯αN+M-2βN+M-2γN+M-20⋯0αN+M-1βN+M-1],
where
(3.18)γ1=k2σ2,αN+M-1=0,βN+M-1=1,γi=αi=k2h2σ2Xi2,2≤i≤N-1,βi=1-kh2σ2Xi2,1≤i≤N-1,αN=kσ2A2h(h+Aδ/(1-δ)),βN=1-kσ2A(1-δ)δh,γN=kσ2A(1-δ)δ(h+Aδ/(1-δ)),αi=kσ2Xi2hi-1(hi+hi-1),βi=1-kσ2Xi2hihi-1,γi=kσ2Xi2hi(hi+hi-1),N+1≤i≤N+M-2.LetB=(bij) be the matrix in ∈ℝ(N+M-1)×(N+M-1) related to the integral part whose entries bij for each fixed i in 1≤i≤N+M-2 are defined by
(3.19)bij={kλXihgi,j,1≤j≤N-1,kλ2Xi(h+δA)gi,N,j=N,kλXiδAgi,jXj2,N+1≤j≤N+M-1,bN+M-1,j=0,1≤j≤N+M-1.
From the previous notation the scheme (3.7)–(3.15) can be written in the form
(3.20)Un+1=(P+B)Un=(P+B)n+1U0,0≤n≤L-1,U0=[f(X1)f(X2)⋯f(XN+M-1)]t.


## 4. 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.Lemma 4.1.
With previous notation, assume that stepsizesk=Δτ, h=ΔX in [0,A] and 0<δ≤1/3, and δ=Δz in ]0,1] satisfy (C1)
k/h2≤1/σ2A2,(C2)
k≤min{δ2/σ2(1-2δ),δh/σ2A(1-δ)}. Then matrix P given by (3.17) is nonnegative.Proof.
From (3.18), for 1≤i≤N-1 one has γi>0 and αi>0 for 2≤i≤N-1. On the other hand, for 1≤i≤N-1, we have that
(4.1)βi≥0,iffkh2≤1σ2Xi2.
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
(4.2)βN≥0,iffk≤δhσ2A(1-δ).
In order to guarantee the nonnegativeness of the remaining entries of matrix P, let us introduce the function
(4.3)H(i)=hihi-1Xi2=δ2[1-(i-(N+1))δ][1-(i-(N-1))δ]>0,
for N+1≤i≤N+M-2. With this notation, we have that βi appearing in (3.18) satisfies
(4.4)βi≥0,iffk≤H(i)σ2.
Note that
(4.5)dH(i)di=2(1-(i-N)δ)δ3[1-(i-(N+1))δ]2[1-(i-(N-1))δ]2.
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
(4.6)dH(i)di>0,N+1≤i≤N+M-2.
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
(4.7)H(N+1)=δ21-2δ,
and condition (4.4) holds true under the condition
(4.8)k≤δ2σ2(1-2δ).
From condition (C2), properties (4.2) and (4.8) are satisfied and matrix P is nonnegative.Note that as matrixB defined by (3.19) is always nonnegative, from Lemma 4.1 and (3.20) starting from nonnegative initial vector U0, the following result is established.Theorem 4.2.
With the hypotheses and notation of Lemma4.1, the solution {uin} of the scheme (3.7)–(3.15) is nonnegative if the initial values ui0≥0,1≤i≤N+M-1.The next result will be used below to guarantee stability.Lemma 4.3.
MatricesP and B defined by (3.17), (3.18), and (3.19) satisfy the following. (1)
Under conditions (C1) and (C2) of Lemma4.1, ∥P∥∞=1.(2)
∥B∥∞≤kλC, where C=1+2/π(exp(-μJ+σJ2/2)/σJ)(1+exp(2μJ)).Proof.
By Lemma4.1, under hypotheses (C1) and (C2), all the entries of matrix P are nonnegative. Thus,
(4.9)|β1|+|γ1|=β1+γ1=1-k2h2σ2X12≤1,|αN+M-1|+|βN+M-1|=αN+M-1+βN+M-1=1.
We also have
(4.10)|αi|+|βi|+|γi|=αi+βi+γi=1,for2≤i≤N+M-2.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
(4.11)∑j=1N+M-1‍|bij|=∑j=1N+M-1‍bij=kλXi[TN(g(XXi),[0,A])+TM(g(AXi1z)Az2,[0,1])],
where
(4.12)TN(g(XXi),[0,A])=h(∑j=1N-1‍gi,j+12gi,N),(4.13)TM(g(AXi1z)Az2,[0,1])=δA(∑j=N+1N+M-1‍gi,jXj2+gi,NXN22)=δ(∑j=1M-1‍g(AXi1zj)Azj2+12g(AXi)A)
are the trapezoidal approximation rules with N and M points, approximating the integrals
(4.14)∫0A‍g(XXi)dX,∫A∞‍g(XXi)dX=∫01‍g(AXi1z)Az2dz,
respectively.LetXq(i) be the maximum of g(X/Xi) in ℝ+, given by
(4.15)Xq(i)=Xiexp(μJ-σJ2);g(Xq(i)Xi)=12πσJexp(-μJ+σJ22).
Note that the integrand g(X/Xi) is increasing in the interval ]0,Xq(i)[ and decreasing for ]Xq(i),+∞[. In order to upper bound (4.12) we consider two cases. Firstly, let us assume that 2h<Xq(i)<A-h and let us denote j0 to be the first integer with 1≤j0≤N-3, such that
(4.16)Xj0+h=Xj0+1≤Xq(i)<Xj0+2=Xj0+2h.
From the properties of the lower Riemann sums, it follows that
(4.17)h∑j=1j0‍gi,j≤∫0Xq(i)‍g(XXi)dX,h∑j=j0+3N‍gi,j≤∫Xq(i)A‍g(XXi)dX.
Taking into account the values gi,j0+1 and gi,j0+2 located just before and after Xq(i), together with (4.17), from (4.12) it follows that
(4.18)TN(g(XXi),[0,A])=h∑j=1j0‍gi,j+h(gi,j0+1+gi,j0+2)+h∑j=j0+3N‍gi,j-h2gi,N≤∫0A‍g(XXi)dX+2hg(Xq(i)Xi).
In an analogous way, for the second situation where Xq(i)≤2h, one gets again (4.18). Finally, if Xq(i)≥A-h is also true that
(4.19)h∑j=1N‍gi,j≤h(gi,N-1+gi,N)+∫0A‍g(XXi)dX≤2hg(Xq(i)Xi)+∫0A‍g(XXi)dX.
Now we will upper bound (4.13). Let us denote
(4.20)Gi(z)=g(AXi1z)Az2,0<z≤1,
with limz→0+Gi(z)=0. Let zq(i) be the maximum of Gi(z):
(4.21)zq(i)=AXiexp(-μJ-σJ2),Gi(zq(i))=Xi22πAσJexp(μJ+σJ22).
In this case we could distinguish the three possible situations 2δ<zq(i)<1-δ; zq(i)≤2δ;zq(i)≥1-δ, and upper bounding the lower Riemann sums relative to (4.13) one gets
(4.22)TM(Gi(z),[0,1])≤2δGi(zq(i))+∫01‍Gi(z)dz=2δGi(zq(i))+∫A∞‍g(XXi)dX.
From (4.11), (4.18), and (4.22) together with the fact that ∫0∞‍g(X/Xi)dX=Xi, one gets for each value of i,
(4.23)∑j=1N+M-1‍|bij|≤kλ[1+2πexp(-μJ+σJ2/2)σJ(hXi+δAXiexp(2μJ))].
Taking into account that h=X1≤Xi≤XN+M-1=A/δ, for 1≤i≤N+M-1, it follows that h/Xi≤1 and (δ/A)Xi≤1. Hence, from (4.23) one gets
(4.24)∥B∥∞≤kλC,
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{uin} 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 {Un} be the corresponding vector form, that is, Un=[u1nu2n⋯uN+M-1n]t of (3.20). We say that {uin} is strongly uniformly ∥·∥∞ stable, if
(4.25)∥Un∥∞≤W∥U0∥∞,0≤n≤L,
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{uin} of the scheme (3.7)–(3.15) is strongly uniformly ∥·∥∞ stable if one satisfies the condition 0<δ≤1/3 together with
(4.26)kh2≤1σ2A2,k≤min{δ2σ2(1-2δ),δhσ2A(1-δ)}.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),
(4.27)∥Un+1∥∞≤(∥P∥∞+∥B∥∞)∥Un∥∞≤(1+kλC)∥Un∥∞.
Hence, from (4.27), Bernouilli’s inequality, and 0≤n≤L,kL=τ≤T,  one gets
(4.28)∥Un∥∞∥U0∥∞≤(1+kλC)n≤1+nkλC≤exp(nkλC)≤exp(λCT).
Thus the conditional strong uniform stability is established.


## 5. 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(uin)=0, where
(5.1)F(uin)=uin+1-uink-12σ2Xi2Δin-λXi(J1,in+J2,in),
and let us write the PIDE (2.9) in the form
(5.2)L(U)=I(U),
where
(5.3)L(U)=∂U∂τ-12σ2X2∂2U∂X2,I(U)=λX(J1+J2),
where J1 and J2 are given by (2.6)–(2.8).Let us denoteUin=U(Xi,τn) to be the value of the theoretical solution of PIDE (5.2). Let A>0 such that Xi<A. We denote by the following expression the local truncation error Tin(U):
(5.4)Tin(U)=F(Uin)-L(Uin)+I(Uin).
In order to prove the consistency, we must show that
(5.5)Tin(U)→0,ash→0,δ→0,k→0.
Assuming that U is twice continuously partially differentiable with respect to τ and using Taylor’s expansions about (Xi,τn) one gets
(5.6)Uin+1-Uink=∂U∂τ(Xi,τn)+kEin(1),
where
(5.7)Ein(1)=12∂2U∂τ2(Xi,η),τn<η<τn+1,|Ein(1)|≤12Win(1)=12max{|∂2U∂τ2(Xi,η)|;τn≤η≤τn+1}.Let us assume thatU admits four times continuous partial derivatives with respect to X, and let us denote
(5.8)Δin(U)=∂2U∂X2(Xi,τn)+h2Ein(2),Ein(2)=112∂4U∂X4(ν,τn),Xi-h<ν<Xi+h,|Ein(2)|≤112Wn(2)=112max{|∂4U∂X4(X,τn)|;0≤X≤A}.
In accordance with [34, page 101] let us explain the local consistency error of J1, see (2.6), by
(5.9)C1,in(h,A)=∫0A‍U(ϕ,τn)g(ϕXi)dϕ-T1,in([0,A]),(5.10)T1,in([0,A])=h∑j=1N-1‍Ujngi,j+h2UNngi,N,(5.11)|C1,in(h,A)|≤h212Amax0≤X≤A{|(U(X,τn)g(XXi))(2)|}=h2A12Win(3),
where (U(X,τn)g(X/Xi))(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 integralJ2, see (2.8), by
(5.12)C2,in(δ,A)=A∫01‍1z2U(Az,τn)g(AXi1z)dz-T2,in([0,1]),T2,in([0,1])=δA(12UNngi,NXN2+∑j=N+1N+M-1‍Ujngi,jXj2),|C2,in(δ,A)|≤Aδ212max0≤z≤1{|(1z2U(Az,τn)g(AXi1z))(2)|}=Aδ212Win(4).
Summarizing, one gets
(5.13)T1n(U)=kEin(1)-h22σ2Xi2Ein(2)+λXi(C1,in(h,A)+C2,in(δ,A)),|Ti(n)(U)|≤h224Wn(2)σ2XN2+h212λAXiWin(3)+δ212λAXiWin(4)+kWin(1).Thus(5.14)Tin(U)=O(h2)+O(δ2)+O(k),
which proves the consistency of the scheme with the PIDE.


## 6. Numerical Results

In the following examples the code was run onMatlab. 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.Figure 1
Satisfying and breaking stability conditions.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 Example6.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.Figure 2
Variation of the jump intensityλ.In the next example, the error is the difference between the numerical solutionV(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 separation point A=30 remains stationary when h decreases. This fact agrees with facts illustrated in [28, pages 15-16].Example 6.3.
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, K=0.00501, and λ=0.1. For A=30, k=0.001, and δ=0.0625, the Figure 3 shows the variation of the absolute error of the solution under changes of the stepsize h.Figure 3
Absolute errors with several values ofh and a fixed δ.The next Example6.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 Example6.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.Figure 4
Absolute errors with several values ofδ and a fixed h.


## 7. 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 Example6.4.



---

*Source: 120358-2012-12-04.xml*