# Positive Solutions of European Option Pricing with CGMY Process Models Using Double Discretization Difference Schemes

Authors: R. Company; L. Jódar; M. Fakharany
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper deals with the numerical analysis of PIDE option pricing models with CGMY process using double discretization schemes. This approach assumes weaker hypotheses of the solution on the numerical boundary domain than other relevant papers. Positivity, stability, and consistency are studied. An explicit scheme is proposed after a suitable change of variables. Advantages of the proposed schemes are illustrated with appropriate examples.


---

## Body

## 1. Introduction

The hypothesis that asset prices behave according to the geometric Brownian motion when one derives the option prices is inconsistent with market prices [1]. This drawback has been overcome using Lévy process models [2–9] allowing the calibration of the model to the option market price and the reproduction of a wide variety of implied volatility skews/smiles; see [10] and [11, chapters 14, 15]. Among the Lévy process models, it is remarkable to distinguish these with finite activity, that is, jump diffusion models [2, 3], and those where the intensity of the jumps is not a finite measure [4–9]. These models are characterized by the fact that option price is given by the solution of a partial integrodifferential equation (PIDE) involving a second-order differential operator part and a nonlocal integral term that presents additional difficulties. In [12] wavelet methods are applied to infinite Lévy models. Monte Carlo approaches are developed by [13, 14]. Interesting analytic-numerical treatments have been introduced in [15–17]. The so-called COS method for pricing European options is presented in [15]. This is based on the knowledge of the characteristic function and its relation with the coefficients of the Fourier-cosine expansion of the density function. In [16], an expansion of the characteristic function of local volatility models with Lévy jumps is developed. The authors in [17] 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.Many authors used FD schemes for solving these PIDE problems [18–28]. Dealing with FD methods for such PIDEs, the following challenges should be addressed, for instance, how to approximate the integral term and how to localize a bounded computational domain in order to consider relevant information like large jumps. In addition, the possible singularities of the integral kernel should be carefully treated [18, 19].The nonlocal character of the integral part involves a dense discretization matrix. In the outstanding paper [18], Cont and Voltchkova presented an explicit-implicit method (explicit into the integral part and implicit into the differential one) to obtain the numerical approximation of viscosity solutions for European and barrier options. An improvable issue of [18] is that, in order to approximate the truncated integral term, they assume a particular behavior of the solution outside of the bounded numerical domain. This last drawback is experienced by most of the authors; see [22, 26, 27].Implicit FD methods for the numerical solution of the CGMY model have been used by Wang et al. [19] who proposed an implicit time stepping method avoiding dense linear systems but involving the iteration methods drawbacks of the implicit methods such as ungranted positivity. They also assume that, for large enough values of S, the solution behaves like Black-Scholes.In [21], the authors use an unconditional ADI FD method and accelerate it using fast Fourier transform (FFT) for jump diffusion models with finite jump intensity. Tavella and Randall in [24] use an implicit time discretization and propose a stationary rapid convergent iterative method to solve the full matrix problem quoted above, but with poor numerical analysis. A generalization of their iterative method to price American options is proposed in [25].One of the most relevant and versatile Lévy models is the one proposed by Carr et al., the so-called CGMY model [8], that belongs to the family of KoBoL models [9]. It is considered a prototype of the general class of models with jumps and enjoys widespread applicability. The CGMY model allows diffusions and jumps of both finite and infinite activity. The CGMY Lévy density is given by
(1)ν(y)={Ce-G|y||y|1+Y,y<0,Ce-M|y||y|1+Y,y>0,
where C>0, G≥0, M≥0, and Y<2. The parameter Y allows controling the fine structure of asset return distribution. For Y<0, the Lévy process is of finite activity; that is, the measure is finite; ∫‍ν(y)dy<∞. For 0≤Y≤1, it is of infinity activity but finite variance; that is, ∫|y|<1‍yν(y)dy<∞. Finally, for 1<Y<2, both the activity and variation are infinite. Note that, for Y=0, one gets the well known Variance Gamma process proposed by Madan and Seneta [29] as a particular case. So CGMY model is an extension of the Variance Gamma model [7].The authors in [22] use FD methods discretizing the equation in space by the collocation method and using explicit difference backward schemes focused on the case of infinite activity and finite variation.Very recently [28] proposed an efficient three-time-level finite difference method for the infinite activity Lévy model. Second-order convergence rate is shown in numerical experiments although the numerical analysis of the method is not developed.The aim of this paper is the construction and numerical analysis of an explicit finite difference scheme of the PIDE for the CGMY model with measureν(y) given by (1):
(2)∂V∂τ=σ22S2∂2V∂S2+(r-q)S∂V∂S-rV+∫-∞+∞‍ν(y)[V(Sey,τ)-V(S,τ)-S(ey-1)∂V∂S]dy,S∈(0,∞),τ∈(0,T],(3)V(S,0)=f(S),S∈(0,∞).HereV(S,τ) is the option price depending on the underlying asset S, the time τ=T-t,σ is the volatility parameter, and r and q are the risk-free interest and the continuous dividend paid by the asset, respectively. The payoff function f(S) for a vanilla call option is given by
(4)f(S)=max(S-E,0),
where E is the strike price.Like [20, 23] for jump diffusion models we transform the original PIDE problem in order to remove the convection term to avoid possible numerical oscillations. With respect to the singularity of the integral kernel quoted above, the jump component in the neighborhood of log jump size zero is approximated by using a Taylor expansion, like [18, 19].The selection of the boundary conditions of the numerical domain, the discretization of the infinite domain of the integration, and matching the discretization of both the differential and the integral part are important challenges. Some authors, like those of [18], assume a particular behavior of the solution outside of the bounded numerical domain. In order to weaken these hypotheses we do not truncate the infinite integral and we use a nonuniform partition of the complete unbounded domain, allowing a proper matching of the discretizations of the differential and integral parts by assuming asymptotic linear behavior of the solution. This strategy involves a double discretization with two spatial stepsize parameters that will allow a better flexibility to improve the approximation in different zones of the domain.This paper is organized as follows. In Section2, the integral part of (2) is approximated in a neighborhood of y=0 to obtain a new PIDE integral part extended outside a neighborhood of y=0. Then a variable transformation is developed in order to remove both the convection and reaction terms of the differential part. Following the idea developed in [20] for the finite activity case, a double discretization explicit scheme for solving PIDE problem (2) is constructed in Section 3 for the infinite activity case affording the new challenges. Positivity and stability of the numerical solutions given by the scheme are studied in Section 4. Consistency of the scheme is treated in Section 5. In Section 6, some illustrative numerical examples show the advantages of the new discretization approach showing how the double discretization allows flexible improvement of the accuracy in different zones of the domain.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 ∥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.The exponential integrals have a major role in evaluating important class of integrals. Lets and z be continuous (real or complex) variables; the exponential integral of order s denoted by Es(z) is given by [30]
(5)Es(z)=∫1∞‍t-sexp(-zt)dt.


## 2. Transformation of the PIDE Problem

We begin this section by removing the singularity of the kernel of the integral term of PIDE (2). Let ε>0 and let us split the real line into two regions R1=[-ε,ε] and R2=(-∞,ε)∪(ε,∞). For the term V(Sey,τ) in R1, taking Taylor expansion for z=Sey about z=S one gets
(6)V(Sey,τ)=V(S,τ)+S(ey-1)∂V∂S+(S(ey-1))22∂2V∂S2+𝒪(y3),|y|<ε.
Taking into account (1) the integral part of (2) can be written as
(7)I(V)=σ2(ε)2S2∂2V∂S2-γ(ε)S∂V∂S-λ(ε)V(S,τ)+∫R2‍ν(y)V(Sey,τ)dy+𝒪(ε3-Y),
where the integrals
(8)σ2(ε)=∫-εε‍ν(y)(ey-1)2dy,γ(ε)=∫R2‍ν(y)(ey-1)dy,λ(ε)=∫R2‍ν(y)dy,
can be evaluated with high accuracy using the exponential integrals ([30, 31], chapter 7). Let us denote
(9)f^(α,M,ε)=Mα-1Γ(1-α)-ε1-αEα(εM),
where Γ denotes the gamma function and Eα is the exponential integral. Considering the expression (1) and the exponential integral (5) for M>2, one gets
(10)σ2(ε)=C∑k=02‍(-1)k(2k)×[f^(1+Y,G+k,ε)+f^(1+Y,M-k,ε)].
Notice that (10) holds for Y∈(0,1)∪(1,2). For the particular case where Y=0, one gets
(11)σ2(ε)=C[ln(G+1)2(M-1)2G(G+2)M(M-2)-∑k=02‍(2k)(-1)k×(E1[ε(G+k)]+E1[ε(M-k)])(G+1)2(M-1)2G(G+2)M(M-2)],
while for Y=1, we have
(12)σ2(ε)=2C[ln(G+2)(M-1)(G+1)(M-2)+E1[ε(G+2)]-E1[ε(G+1)]-E1[ε(M-2)]+E1[ε(M-1)](G+2)(M-1)(G+1)(M-2)]+CG[∑k=02‍(2k)lnG(G+2)(G+1)2+∑k=02‍(2k)(-1)kE1[ε(G+k)]]-Cεe-Gε(e-ε-1)2+CM[∑k=02‍(2k)lnM(M-2)(M-1)2+∑k=02‍(2k)(-1)kE1[ε(M-k)]]-Cεe-Mε(eε-1)2.
For the remaining integrals in (8), we have
(13)γ(ε)=Cε-Y∑k=01‍(1k)(-1)1-k×[E1+Y(ε(G+k))+E1+Y(ε(M-k))],(14)λ(ε)=Cε-Y(E1+Y(Gε)+E1+Y(Mε)).
Hence, the problem (2) takes the following form:
(15)∂V∂τ=σ^22S2∂2V∂S2+(r-q-γ(ε))S∂V∂S-(r+λ(ε))V+∫R2‍ν(y)V(Sey,τ)dy+𝒪(ε3-Y),
where σ^2=σ^2(ε)=σ2+σ2(ε).In order to remove the convection and reaction terms from (15), let us introduce the following transformation of variables:
(16)x=exp[(r-q-γ(ε))τ]S,U(x,τ)=exp[(r+λ(ε))τ]V(S,τ).
Hence the problem (15) is approximated by the following form:
(17)∂U∂τ=σ^22x2∂2U∂x2+J,x∈(0,+∞),τ∈(0,T],U(x,0)=f(x),x∈(0,+∞),
where
(18)J=J(x,τ,ε)=∫R2‍ν(y)U(xey,τ)dy=∫-∞-ε‍ν(y)U(xey,τ)dy+∫ε∞‍ν(y)U(xey,τ)dy.
Finally in order to combine both discretizations of the differential and integral part, we use ϕ=xey to change the integrand J as follows:
(19)J=J1+J2=∫0xe-ε‍g(x,ϕ)U(ϕ,τ)dϕ+∫xeε∞‍g(x,ϕ)U(ϕ,τ)dϕ,
where g(x,ϕ)=ν(ln(ϕ/x))/ϕ. For evaluating the integrals in all the positive real line, let us introduce a parameter A>0 that separates [0,∞) into [0,A]∪[A,∞). The point A can be chosen according to the criteria used by [18, 32, 33] to truncate the numerical domain. For instance, in [26] one takes A=4E and in [20] one takes A=3E. To evaluate the integrals related to x>A, they are transformed to finite integrals by using the substitution z=A/ϕ consequently, obtaining integrals of the form
(20)∫αβ‍g(x,ϕ)U(ϕ,τ)dϕ=A∫ϕ0ϕ1‍g(x,Az)U(Az,τ)dzz2,
where ϕ0=A/β, ϕ1=A/α. In particular if β→∞ then ϕ0=0. Hence, the problem (17) takes the form
(21)∂U∂τ=σ^22x2∂2U∂x2+J,x∈(0,∞),τ∈(0,T],(22)U(x,0)=f(x),x∈(0,∞).


## 3. Numerical Scheme Construction

In this section a difference scheme for the problem (21)-(22) is designed. For the time variable, given τ∈(0,T], let k be the time-step discretization k=Δτ=τ/L and τl=lk, 0≤l≤L, with L integer. With respect to the spatial variable x and for an arbitrary fixed A>0, we divide the interval [0,A] into N equal intervals with a spatial step h=Δx=A/N, with xi=ih, 0≤i≤N. Note that the unbounded domain [A,∞) is transformed into (0,1] by the above quoted change z=A/x. Thus a uniform distributed mesh partition of the interval (0,1] of the form zi=iδ, δ=1/M, and 0<i≤M is mapped into a nonuniform mesh partition of [A,∞), xi=A/zN+M-i, N≤i≤N+M-1. Hence, we have
(23)xi={ih,0≤i≤N,A1-(i-N)δ,N≤i≤N+M-1.
Let us denote U(xi,τl)≈uil, 0≤i≤N+M-1, 0≤l≤L,
(24)∂U∂τ(xi,τl)≈uil+1-uilk,∂2U∂x2(xi,τl)≈Δil,
where(25)Δil≡{ui-1l-2uil+ui+1lh2,0≤i<N,2[ui-1lhi-1(hi+hi-1)-uilhihi-1+ui+1lhi(hi-1+hi)],N≤i≤N+M-2,and hi=xi+1-xi>0. With respect to the approximation of (19), note that for each xi we need to evaluate two integrals corresponding to [0,xie-ε] and [xieε,∞), denoted by Ji,rl=Jr(xi,τl,ε), r=1,2. Let i1(i) be the biggest j with 0≤j≤N+M-1 such that xj≤xie-ε and let i2(i) be the first j such that xieε≤xj. Then the expression (19) for the point (xi,τl) has the following form:
(26)Jil=Ji,1l+Ji,2l,Ji,1l=∫0xi1‍g(xi,ϕ)U(ϕ,τl)dϕ+∫xi1xie-ε‍g(xi,ϕ)U(ϕ,τl)dϕ,Ji,2l=∫xi2∞‍g(xi,ϕ)U(ϕ,τl)dϕ+∫xieεxi2‍g(xi,ϕ)U(ϕ,τl)dϕ.
Then we apply the trapezoidal rule for the integrals over (0,xi1] and [xi2,∞) because of (20) and using the first mean value theorem for integrals [34, page 1063], the two remaining integrals are approximated by
(27)∫xi1xie-ε‍g(xi,ϕ)U(ϕ,τl)dϕ≈Iil=(∫xi1xie-ε‍g(xi,ϕ)dϕ)ui1l=g~i,i1ui1l,∫xieεxi2‍g(xi,ϕ)U(ϕ,τl)dϕ≈Hil=(∫xieεxi2‍g(xi,ϕ)dϕ)ui2l=g~i,i2ui2l.
Let us denote
(28)gi,j=ν(lnxj/xi)xj.
Depending on the location of xi for each i with 1<i≤N+M-2, we approximate Jil given by (26) in the following form.Case 1 (xi<Ae-ε).
Note that in this casexi1<A, and thus Ji,1l is approximated by h∑j=1i1-1‍gijujl+((h/2)gi,i1+g~i,i1)ui1l. Also one has xi2<A in the domain of the integral Ji,2l and is approximated by ((h/2)gi,i2+g~i,i2)ui2l+h∑j=i2+1N-1‍gijujn+(giN/2)(h+Aδ)uNl+(δ/A)∑j=N+1N+M-1‍gijxj2ujl, taking into account (20) for xj>A. Thus
(29)J^il=h∑j=1i1-1‍gijujl+(h2gi,i1+g~i,i1)ui1l+(h2gi,i2+g~i,i2)ui2l+h∑j=i2+1N-1‍gijujl+giN2(h+Aδ)uNl+δA∑j=N+1N+M-1‍gijxj2ujl.Case 2 (Ae-ε≤xi<Aeε).
Asxi1<A and xi2≥A, the approximation of Jil becomes
(30)J^il=h∑j=1i1-1‍gijujl+(h2gi,i1+g~i,i1)ui1l+(xi22δ2Agi,i2+g~i,i2)ui2l+δA∑j=i2+1N+M-1‍gijxj2ujl.Case 3 (xi≥Aeε).
Herexi2>xi1≥A and the approximation of Jil is given by
(31)J^il=h∑j=1N-1‍gijujl+giN2(h+Aδ)uNl+δA∑j=N+1i1-1‍gijxj2ujl+(xi12δ2Agi,i1+g~i,i2)ui1l+(xi22δ2Agi,i2+g~i,i2)ui2l+δA∑j=i2+1N+M-1‍gijxj2ujl.
Assuming that U(ϕ,τ) tends to zero at least linearly as ϕ tends to zero one has g(x,ϕ)u(ϕ,τ)→0 by (1) and (28). On the other hand, assuming linear behavior of the solution for large values of ϕ, the integrand of (20) g(x,A/z)U(A/z,τ)(1/z2)→0, as z→0. Thus, both the terms involving u0l and uN+Ml do not appear in the expressions of (29)–(31). Taking into account (25)–(31) the resulting difference scheme for the PIDE problem (21) takes the form
(32)uil+1=uil+kσ^22xi2Δil+kJ^il,1≤i≤N+M-2.
In order to obtain a complete difference scheme, we include the initial and boundary conditions. From (22), we have
(33)ui0=max(xi-E,0),1≤i≤N+M-1.
On the other hand, for a vanilla call option the boundary condition for i=0 is
(34)u0l=0,0≤l≤L,
and by assuming the linear behavior of the solution for large values of the spatial variable, we have ∂2U/∂x2→0 and thus ΔN+M-1l=0 and the null integral term approximation JN+M-1=0, for all time level l. Thus from (32) for i=N+M-1, one gets
(35)uN+M-1l+1=uN+M-1l=uN+M-10,0≤l≤L-1.
For the sake of convenience to study the stability, we now introduce the vector formulation of scheme (32)–(35). Let us denote the vector in ℝN+M-1 as
(36)Ul=[u1lu2l⋯uN+M-1l]t,
and let P=(pij) be a tridiagonal matrix in ℝ(N+M-1)×(N+M-1) related to the differential part, defined by
(37)pij={αi,2≤i≤N+M-1,j=i-1,βi,1≤i≤N+M-1,j=i,γi,1≤i≤N+M-2,j=i+1,0,otherwise,
where
(38)γ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.
Let B=(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
(39)bij={kbij(1),1≤i≤i1(N)-1,kbij(2),i1(N)≤i≤i2(N)-1,kbij(3),i2(N)≤i≤N+M-2,0,i=M+N-1,
where
(40)bij(1)={hgij,1≤j≤i1-1,h2gi,j+g~i,j,j=i1,i2,0,i1+1≤j≤i2-1,hgij,i2+1≤j≤N-1,12(h+Aδ)giN,j=N,δAxj2gij,N+1≤j≤N+M-1,(41)bij(2)={hgij,1≤j≤i1-1,h2gi,j+g~i,j,j=i1,0,i1+1≤j≤i2-1,xi22δ2Agi,j+g~i,j,j=i2,δAxj2gij,i2+1≤j≤N+M-1,(42)bij(3)={hgij,1≤j≤N-1,12(h+Aδ)giN,j=N,δAxj2gij,N+1≤j≤i1-1,δxj22Agi,j+g~i,j,j=i1,i2,0,i1+1≤j≤i2-1,δAxj2gij,i2+1≤j≤N+M-1.
Hence scheme (32)–(35) can be written in the form
(43)Ul+1=(P+B)Ul=(P+B)lU0,0≤l≤L-1,U0=[f(x1)f(x2)⋯f(xN+M-1)]t.


## 4. Positivity and Stability of the Numerical Solution

The price of contracts modelled by PIDE must be a nonnegative value. Our objective here is to demonstrate that the solution of scheme (32)–(35) is conditionally nonnegative and stable.First we study the positivity of the matrixP. The following lemma has been proved in [20].Lemma 1.
With previous notation, assume that stepsizesk=Δτ, h=Δx in [0,A] and 0<δ≤1/3, δ=Δz in (0,1], satisfy (C1)
k/h2≤(1/σ^2A2);(C2)
k≤min{δ2/σ^2(1-2δ),δh/σ^2A(1-δ)}.
Then matrixP given by (37) is nonnegative.Note that as the matrixB defined by (40)-(39) is always nonnegative, from Lemma 1 and (43) starting from nonnegative initial vector U0, the following result is established.Theorem 2.
With the hypotheses and notation of Lemma1, the solution {uil} of scheme (32)–(35) is nonnegative if the initial values ui0≥0, 1≤i≤N+M-1.The next result will be used below to guarantee stability.Lemma 3.
Let matricesP and B be defined by (37)–(39), and let ε>0; then the following results hold. (1)
Under conditions (C1) and (C2) of Lemma1, ∥P∥∞=1.(2)
∥B∥∞≤k(λ(ε)+1), where λ(ε) is defined by (14).By [20, lemma 2], part 1 is proved. Since the norm of B is given by
(44)∥B∥∞=maxi∑j=1N+M-1‍|bij|=maxi∑j=1N+M-1‍bij,1≤i≤N+M-1,
if m~ denotes the row containing the maximum of (44), one gets
(45)∥B∥∞=∑j=1N+M-1‍bm~j=k∑j=1N+M-1‍bm~j(r),r=1,2,or3,
the elements of the summation in (45) are given by (40)-(39). To upper bound (45), we apply the change of variables y=xm~eϕ in (8), resulting in
(46)λ(ε)=∫0xm~e-ε‍g(xm~,ϕ)dϕ+∫xm~eε∞‍g(xm~,ϕ)dϕ,
which coincides with (19) when U(ϕ,τ)=1. Hence from (19), (45), and (46), we conclude that ∑j=1N+M-1‍bm~j(r) is an approximation for λ(ε). Thus, for small enough h and δ, one gets [35]
(47)∑j=1N+M-1‍bm~j(r)<λ(ε)+1.
Hence
(48)∥B∥∞<k(λ(ε)+1),
independently of the value of the size of matrix B.There are many definitions of stability in the literature; we recall our concept of stability in the next definition.Definition 4.
Let{uil} be a numerical solution of the PIDE (21), (22) computed from scheme (32)–(35) with stepsizes h=Δx in [0,A], δ=Δz in ]0,1], and k=Δτ in [0,τ]. Let {Ul} be the corresponding vector form; that is, Ul=[u1lu2l⋯uN+M-1l]t; of (43). One says that {uil} is strongly uniformly ∥·∥∞ stable if
(49)∥Ul∥∞≤W∥U0∥∞,0≤l≤L,
where W>0 is independent of l, h, δ, and k.If the property (49) is satisfied for appropriate relationships between the stepsizes h, δ and k, then one says that the strong uniform stability is conditional.Theorem 5.
With the previous notation, the numerical solution{uil} of scheme (32)–(35) is strongly uniformly ∥·∥∞ stable if one satisfies the condition 0<δ≤1/3 together with
(50)kh2≤1σ^2A2,k≤min{δ2σ^2(1-2δ),δhσ^2A(1-δ)}.Proof.
Note that scheme (32)–(35) is equivalent to the vector form scheme (43). Under condition (50), by Lemma 3 one gets, after taking norms in (43),
(51)∥Ul+1∥∞≤(∥P∥∞+∥B∥∞)∥Ul∥∞≤(1+k(λ(ε)+1))∥Ul∥∞.
Hence, from (51) and that 0≤l≤L, kL=τ≤T,
(52)∥Ul∥∞∥U0∥∞≤(1+k(λ(ε)+1))l≤exp(lk(λ(ε)+1))≤exp(T(λ(ε)+1)).
Thus the conditional strong uniform stability is established.


## 5. Consistency

A numerical scheme is consistent with a PIDE if an exact theoretical solution of the PIDE approximates well the difference scheme as the stepsizes discretization tend to zero [36, 37].After recalling the definition of consistency, let us write (32) in the form
(53)Fil(u)=uil+1-uilk-σ^22xi2Δil-J^il=0.
Let us denote Uil=U(xi,τl) as the value of the theoretical solution of (21) and let A>0 such that xi<Ae-ε, and let us write the PIDE (21) as
(54)L(U)=J(U),
where
(55)L(U)=∂U∂τ-σ^22x2∂2U∂x2,J(U)=J=J1+J2.
The local truncated error Til(U) at (xi,τl) is defined by
(56)Til(U)=Fil(U)-L(Uil)+J(Uil).
In order to prove the consistency, we must show that
(57)Til(U)⟶0,ash⟶0,δ⟶0,k⟶0.
Assuming that U is twice continuously partially differentiable with respect to τ and four times partially differentiable with respect to x, and using Taylor’s expansion about (xi,τl), it follows that
(58)Uil+1-Uilk=∂U∂τ(xi,τl)+kEil(1),Eil(1)=12∂2U∂τ2(xi,ζ),τl<ζ<τl+1,|Eil(1)|≤12Wil(1)=12max{|∂2U∂τ2(xi,ζ)|;τl≤ζ≤τl+1},Δil=∂2U∂x2(xi,τl)+h2Eil(2),Eil(2)=112∂4U∂x4(ζ~,τl),xi-h<ζ~<xi+h,|Eil(2)|≤112Wil(2)=112max{|∂4U∂x4(ζ~,τl)|;xi-h≤ζ~≤xi+h}.
In accordance with [36, page 101] let us denote the local consistency error of Ji,1l (see (29)) by
(59)Ci,1l=Ji,1l-Ti,1l([0,xie-ε]),Ti,1([0,xie-ε])=h∑j=1i1-1‍gijUjl+(h2gi,i1+g~i,i1)Ui1l.
By (26) and (29), the local consistency error for Ji,2l is given by
(60)Ci,2l=Ji,2l-(Ti,2l([xieε,A])+Ti,3l([0,1])),
where
(61)Ti,2l([xieε,A])=(h2gi,i2+g~i,i2)Ui2l+h∑j=i2+1N-1‍gijUjl+h2giNUNl,Ti,3l([0,1])=δA(12giNxN2UNl+∑j=N+1N+M-1‍gijxj2Ujl).
From the first mean value theorem for integrals [34, page 1063], one gets
(62)I(xi,ε)=∫xi1xie-ε‍g(xi,ϕ)U(ϕ,τl)dϕ=(∫xi1xie-ε‍g(xi,ϕ)dϕ)U(c,τl)=g~i,i1U(c,τl),xi1<c<xie-ε,
and since
(63)U(c,τl)=U(xi1,τl)+(c-xi1)∂U∂x(ξ,τl),xi1<ξ<c,
it follows that
(64)|I(xi,ε)-g~i,i1Ui1l|≤g~i,i1hΛil(1)≤h2Wil(3),
where
(65)Λil(1)=max{|∂U∂x(x,τl)|;xi1≤x≤xie-ε},Wil(3)=Λil(1)max{g(xi,x);xi1≤x≤xie-ε}.
Analogously,
(66)|H(xi,ε)-g~i,i2Ui2l|≤g~i,i1hΛil(1)≤h2Wil(4),H(xi,ε)=∫xieεxi2‍g(xi,ϕ)U(ϕ,τl)dϕ,Wil(4)=(max{|∂U∂x(x,τl)|;xie-ε≤x≤xi2})×(max{g(xi,x);xie-ε≤x≤xi2}).
Let Wil(5), Wil(6), and Wil(7) be defined as
(67)Wil(5)=sup{|(g(xi,x)U(x,τl))(2)|;0<x≤xie-ε},Wil(6)=sup{|(g(xi,x)U(x,τl))(2)|;xieε≤x≤A},Wil(7)=sup{|(g(xi,Az)U(Az,τl)1z2)(2)|;0<x≤1},
where the second derivatives appearing in (67) are taken with respect to the variable x for Wil(5) and Wil(6), and with respect to the variable z for Wil(7). From the expression of the error of the trapezoidal rule, [35, page 54], (59)–(67), one gets
(68)|Ci,1l|≤h2(Wil(3)+xie-ε12Wil(5)),|Ci,2l|≤h2(Wil(4)+112(A-xie-ε)Wil(6))+Aδ212Wil(7).
Summarizing, one gets
(69)Til(U)=kEil(1)-h22σ^2xi2Eil(2)+Ci,1l+Ci,2l,|Til(U)|≤h2(Wil(2)24+Wil(3)+Wil(4)+xie-ε12Wil(5)+112(A-xie-ε)Wil(6)Wil(2)24+Wil(3)+Wil(4)+xie-ε12)+Aδ212Wil(7)+kWil(1).
Thus,
(70)Til(U)=𝒪(h2)+𝒪(δ2)+𝒪(k),
showing the consistency of the scheme with PIDE.


## 6. Numerical Results

In the following examples the Matlab program has been used. The first example reveals the effect of Yor parameter on the option price.Example 6.
Consider the vanilla call option problem (2)–(4) under CGMY process with parameters T=1, E=80, A=3E, σ=0.2, r=0.01, q=0, C=0.08, G=M=25.04, ε=0.05, N=100, δ=0.15, and k=0.002. Figure 1 exhibits the variation of the option price V versus the underlying asset at various values of Yor parameter.Figure 1
The variation ofV with S for several values of Y.The next example illustrates the importance of positivity conditions given by Lemma1.Example 7.
Here in this example the parameters have been selected as follows:T=1, E=80, A=3E, σ=0.2, r=0.01, q=0, C=1, G=20, M=30, Y=1.5, ε=0.1, N=100, and δ=0.15. Positivity conditions hold for k=0.002, while for k=0.01, the positivity conditions are broken and the values of the option price become unreliable as shown in Figure 2.Figure 2
The effect of positivity conditions onV.Next example shows that the variation of the absolute and relative error of the solution in light of the stability and positivity conditions hold at the strike for two cases: first, for several values of the stepsize discretizationh and second, for different values of the parameter ε.Example 8.
Consider the European call option for CGMY process with the following values:C=1, G=M=5, E=100, T=1, r=0.1, q=0, k=0.001, δ=0.1, and A=3E, for several values of Yor parameter Y=0.5,1.5, and 1.98.
We consider the evaluation of the price option at the strike andτ=T. Table 1 reveals the deviation between our numerical solutions and the reference values used in [15, tables 8–10] for different stepsizes h and fixed ε=0.12. Notice that the numerical solution exhibits the expected second-order convergence rate.
Table2 shows the deviation for several values of ε, while h=0.5.Table 1
Errors and convergence rates.
(a)
h
Y
=
0.5
Y
=
1.5
Absolute error
Relative error
α
Absolute error
Relative error
α
0.8
4.38
×
1
0
-
4
2.2
×
1
0
-
5
—
7.35
×
1
0
-
5
1.48
×
1
0
-
6
—
0.4
1.16
×
1
0
-
4
5.85
×
1
0
-
6
1.92
1.9
×
1
0
-
5
3.92
×
1
0
-
7
1.952
0.2
2.95
×
1
0
-
5
1.49
×
1
0
-
6
1.98
4.79
×
1
0
-
6
9.62
×
1
0
-
8
1.988
(b)
h
Y
=
1.98
Absolute error
Relative error
α
0.8
3.87
×
1
0
-
5
3.87
×
1
0
-
7
—
0.4
9.76
×
1
0
-
6
9.76
×
1
0
-
8
1.9873
0.2
2.46
×
1
0
-
6
2.46
×
1
0
-
8
1.9882Table 2
Errors due to the variation ofε.
(a)
ε
Y
=
0.5
Y
=
1.5
Absolute error
Relative error
Absolute error
Relative error
0.8
3.91
×
1
0
-
3
1.97
×
1
0
-
4
6.37
×
1
0
-
4
1.28
×
1
0
-
5
0.4
7.18
×
1
0
-
4
3.62
×
1
0
-
5
8.54
×
1
0
-
5
1.72
×
1
0
-
6
0.2
9.32
×
1
0
-
6
4.7
×
1
0
-
7
7.16
×
1
0
-
6
1.44
×
1
0
-
7
(b)
ε
Y
=
1.98
Absolute error
Relative error
0.8
4.19
×
1
0
-
4
4.19
×
1
0
-
6
0.4
5.76
×
1
0
-
5
5.76
×
1
0
-
7
0.2
5.92
×
1
0
-
6
5.92
×
1
0
-
8In the next example, we consider the Variance Gamma model as a particular case (Y=0) of CGMY model for which the exact solution is known [38].Example 9.
Consider a call option under the Variance Gamma process with parametersC=1, G=M=25, T=1, r=0.01, q=0, σ=0.2, ε=0.12, E=10, A=3E, k=0.01, and δ=0.15. Figure 3 displays the associated error of the numerical solution for several values of the stepsize h.
The next example shows that the double discretization strategy reduces the error near the parameterA by changing the stepsize δ.Figure 3
The associated error for various values ofh.Example 10.
Consider Example9 with fixed h=0.5; Figure 4 shows the variation of the error of the numerical solution for various values of δ. Notice that the error decreases near the right boundary A of the numerical domain by decreasing the stepsize δ, while the error near the strike E remains stationary.Figure 4
The associated error for several values ofδ.



---

*Source: 517480-2013-11-21.xml*