# Unconditional Positive Stable Numerical Solution of Partial Integrodifferential Option Pricing Problems

Authors: M. Fakharany; R. Company; L. Jódar
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

This paper is concerned with the numerical solution of partial integrodifferential equation for
option pricing models under a tempered stable
process known as CGMY model. A double
discretization finite difference scheme is used for the treatment
of the unbounded nonlocal integral term. We also introduce in the
scheme the Patankar-trick to guarantee unconditional nonnegative
numerical solutions. Integration formula of open type is used in
order to improve the accuracy of the approximation of the integral
part. Stability and consistency are also studied. Illustrative examples are included.


---

## Body

## 1. Introduction

In Black-Scholes model, it is assumed that the probability distribution of the stock price is lognormal and the instantaneous log return is a geometric Brownian motion. However the market for the options shows that the geometric Brownian model for the underlying asset leads to underprice or overprice for these options [1]. Several models have been proposed to overcome these shortcomings such as models where the volatility follows a stochastic process like Heston model [2] and models incorporating jumps in the underlying asset following Lévy processes [3] and [4, chap. 14, 15]. There are two main features for the Lévy processes: first, models with finite activity, that is, jump diffusion models [5, 6], and, second, models with infinite intensity measure [7–9]. The numerical solution of these models can be achieved through three main techniques: by partial integrodifferential equation (PIDE) methods, as a solution of numerical integration, and by Monte Carlo simulation. Based on the close relation between the characteristic function of the probability density function and the Fourier-cosine expansion, the so-called COS method has been used to obtain the option pricing for European options [10].In this paper we focus on PIDE approach. The PIDE valuating the option price presents a differential part with reaction, convection, and diffusion terms, while the nonlocal integral part is extended over an infinite or semi-infinite interval. Several finite difference (FD) schemes have been proposed to solve numerically these PIDE problems [11–20]. In order to implement FD methods, there are many challenges to face such as how to treat the unbounded domain for the spatial variable and the possible singularity of the kernel of the integral term. Furthermore, a way of numerical integration has to be chosen and the discretizations of both the integral part and the differential part have to be matched [16, 19].The model proposed by Carr et al., the so-called CGMY model, is one of the most practical and adaptable Lévy models [9]. In the CGMY model, the diffusions and jumps can be of finite or infinite activity. Its Lévy density is given by
(1)
ν
y
=
C
e
-
G
|
y
|
y
1
+
Y
,
y
<
0
,
C
e
-
M
|
y
|
y
1
+
Y
,
y
>
0
,
where C
>
0 measures the overall level of activity, G
≥
0 and M
≥
0 measure the skewness, and Y
<
2 controls 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
)
d
y
<
∞. For 0
≤
Y
≤
1, it is of infinite activity but finite variance; that is, ∫
|
y
|
<
1
‍
y
ν
(
y
)
d
y
<
∞. Finally, for 1
<
Y
<
2, both the activity and variation are infinite. Beside that for Y
=
0 one gets the well known Variance Gamma process proposed by Madan and Milne [8] as a special case.In the innovative paper [16], Cont and Voltchkova provided implicit discretization for the differential part and explicit step for the integral part after truncating its domain. The singularity of the integral kernel and the nonsmoothness of initial conditions are treated using the viscosity solutions and applying this technique to Merton and Variance Gamma models. It is worth mentioning that, in the Variance Gamma model, the integral part is split into singular and nonsingular parts.An implicit FD method for CGMY model has been proposed by Wang et al. in [19] and a semi-Lagrangian discretization has been used for the drift term. On the other hand the integral part is split into singular and nonsingular parts; in the singular part the singularity has been removed using Taylor expansion; after that the trapezoidal rule is implemented. In order to adapt the unknown option price of the integral part to the computational mesh points, it is approximated using the upwind quadratic interpolation.The option pricing for jump diffusion models with finite jump intensity has been treated using ADI finite difference method, accelerated by the fast Fourier transformation [13]. In [20], a three-time-level finite difference method is proposed showing a second order convergence rate in the numerical experiments for infinite activity models. However, the authors in [20] focus the interest on computational issues more than the numerical analysis.An explicit scheme has been used in [15], applying the trapezoidal rule to approximate the integral term after removing the singularity of the kernel and including the unbounded domain using a double discretization technique [14]. Furthermore, the authors provided conditions to guarantee the positivity and stability of the numerical solution.Dealing with option pricing models, ensuring positive solutions is a necessary requirement. Our objective is to construct a stable and conditionally consistent numerical scheme that guarantees positive solutions for the PIDE governing the CGMY model(2)
∂
C
∂
τ
=
σ
2
2
S
2
∂
2
C
∂
S
2
+
r
-
q
S
∂
C
∂
S
-
r
C
+
∫
-
∞
+
∞
‍
ν
y
C
S
e
y
,
τ
-
C
S
,
τ
∂
C
∂
S
-
S
e
y
-
1
∂
C
∂
S
d
y
,
S
∈
0
,
∞
,
τ
∈
0
,
T
,
(3)
C
(
S
,
0
)
=
f
(
S
)
,
S
∈
(
0
,
∞
)
,
with measure ν
(
y
) given by (1). Here C
(
S
,
τ
) is the option price depending on the underlying asset S, the time τ
=
T
-
t, σ is the volatility parameter, 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. Numerical methods that guarantee the positive of solutions for parabolic equations have been studied in [21, 22] following the idea initiated by Patankar, the so-called Patankar-trick [23].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 [19]. Then a variable transformation is developed in order to remove both the convection and reaction terms of the differential part. Following partially the ideas about double discretization developed in [14, 15] an explicit scheme is constructed in Section 3. Unconditionally positivity and stability of the numerical solutions are shown 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. The paper ends with a conclusion section.


## 2. Removing the Reaction and Convection Terms of the PIDE Problem

This section begins with removing the singularity of the kernel of the integral part of PIDE (2). In order to achieve this aim, we split the real line into two regions depending on a parameter ε
>
0; that is, R
1
=
[
-
ε
,
ε
] and R
2
=
(
-
∞
,
ε
)
∪
(
ε
,
∞
); see [16]. By using Taylor expansion of C
(
S
e
y
,
τ
) for z
=
S
e
y in R
1 about z
=
S, one gets a convergency of order O
(
ε
3
-
Y
); see [19].On the other hand, the convection and reaction terms in (2) can be removed, obtaining a simpler PIDE with further numerical advantages, with the following transformation; see [15]:
(5)
x
=
exp
⁡
[
(
r
-
q
-
γ
(
ε
)
)
τ
]
S
,
U
(
x
,
τ
)
=
exp
⁡
[
(
r
+
λ
(
ε
)
)
τ
]
V
(
S
,
τ
)
,
such that
(6)
γ
ε
=
∫
R
2
‍
ν
(
y
)
(
e
y
-
1
)
d
y
,
λ
ε
=
∫
R
2
‍
ν
(
y
)
d
y
.
Under transformation (5) problem (2)-(3) is approximated to the following form:
(7)
∂
U
∂
τ
=
σ
^
2
2
x
2
∂
2
U
∂
x
2
+
J
(
U
)
,
x
∈
(
0
,
+
∞
)
,
τ
∈
(
0
,
T
]
,
(8)
U
(
x
,
0
)
=
f
(
x
)
=
max
⁡
(
x
-
E
,
0
)
,
x
∈
(
0
,
+
∞
)
,
where
(9)
σ
^
2
=
σ
2
+
σ
2
(
ε
)
,
σ
2
(
ε
)
=
∫
-
ε
ε
‍
ν
(
y
)
e
y
-
1
2
d
y
,
(10)
J
U
=
J
(
x
,
τ
,
ε
)
=
∫
R
2
‍
ν
(
y
)
U
(
x
e
y
,
τ
)
d
y
=
∫
-
∞
-
ε
‍
ν
(
y
)
U
(
x
e
y
,
τ
)
d
y
+
∫
ε
∞
‍
ν
(
y
)
U
(
x
e
y
,
τ
)
d
y
.
In order to match both discretizations of the differential and integral part, let us consider the substitution ϕ
=
x
e
y in (10). Then, the resulting expression can be written in the more compact form
(11)
J
(
U
)
=
∫
0
∞
‍
g
(
x
,
ϕ
)
U
(
ϕ
,
τ
)
d
ϕ
,
where the new kernel g
(
x
,
ϕ
) takes the form
(12)
g
(
x
,
ϕ
)
=
ν
(
ln
⁡
(
ϕ
/
x
)
)
ϕ
,
0
<
ϕ
≤
x
e
-
ε
,
0
,
x
e
-
ε
<
ϕ
<
x
e
ε
,
ν
(
ln
⁡
(
ϕ
/
x
)
)
ϕ
,
ϕ
≥
x
e
ε
.
In order to evaluate integral (11) without truncation, following the double discretization technique developed in [14, 15], let us introduce a parameter A
>
0 splitting [
0
,
∞
) into [
0
,
A
]
∪
[
A
,
∞
). The point A can be chosen according to the criteria used by [16, 24, 25] to truncate or split the numerical domain. For instance, in [17] one takes A
=
4
E and in [14] one takes A
=
3
E. The unbounded integral part related to x
>
A is transformed to a finite one by means of the change z
=
A
/
ϕ obtaining an integral of the form
(13)
∫
α
β
‍
g
(
x
,
ϕ
)
U
(
ϕ
,
τ
)
d
ϕ
=
A
∫
ϕ
0
ϕ
1
‍
g
(
x
,
A
z
)
U
(
A
z
,
τ
)
d
z
z
2
,
where ϕ
0
=
A
/
β and ϕ
1
=
A
/
α. In particular if β
→
∞ then ϕ
0
=
0. Since there is a wide class of integrals that can be evaluated using exponential integrals, we recall the definition of the exponential integrals. Let s and z be continuous (real or complex) variables; the exponential integral of order s denoted by E
s
(
z
) is given by [26]
(14)
E
s
(
z
)
=
∫
1
∞
‍
t
-
s
exp
⁡
(
-
z
t
)
d
t
.
From computational point of view, the integrals in (6) and (9) can be evaluated efficiently using exponential integrals [26] and [27, chapter 7], as follows:
(15)
γ
ε
=
C
ε
-
Y
·
∑
k
=
0
1
‍
1
k
-
1
1
-
k
·
E
1
+
Y
ε
G
+
k
+
E
1
+
Y
ε
M
-
k
,
(16)
λ
(
ε
)
=
C
ε
-
Y
E
1
+
Y
G
ε
+
E
1
+
Y
M
ε
,
(17)
σ
2
ε
=
C
∑
k
=
0
2
‍
-
1
k
2
k
·
f
^
1
+
Y
,
G
+
k
,
ε
+
f
^
1
+
Y
,
M
-
k
,
ε
,
M
>
2
,
where
(18)
f
^
(
α
,
x
,
ε
)
=
x
α
-
1
Γ
(
1
-
α
)
-
ε
1
-
α
E
α
(
ε
x
)
,
Γ denotes the gamma function and E
α is the exponential integral (14). Notice that (17) holds for Y
∈
(
0,1
)
∪
(
1,2
), while for particular cases Y
=
0 and Y
=
1 one can find the expression of σ
2
(
ε
) in [15].


## 3. Computing the Numerical Solution

In this section a positivity-preserving explicit difference scheme for problem (7)–(13) is constructed. For the time variable, given τ
∈
(
0
,
T
], let k be the time-step discretization k
=
Δ
τ
=
T
/
L and τ
l
=
l
k, 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 x
i
=
i
h, 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 z
i
=
i
δ, δ
=
1
/
M, 0
<
i
≤
M is mapped into a nonuniform mesh partition of [
A
,
∞
), x
i
=
A
/
z
N
+
M
-
i, N
≤
i
≤
N
+
M
-
1. Hence, we have
(19)
x
i
=
i
h
,
0
≤
i
≤
N
,
A
1
-
(
i
-
N
)
δ
,
N
≤
i
≤
N
+
M
-
1
.
Let us denote U
(
x
i
,
τ
l
)
≈
u
i
l, 0
≤
i
≤
N
+
M
-
1, 0
≤
l
≤
L,
(20)
∂
U
∂
τ
(
x
i
,
τ
l
)
≈
u
i
l
+
1
-
u
i
l
k
∂
2
U
∂
x
2
(
x
i
,
τ
l
)
≈
Δ
i
l
,
where, in the approximation of the second partial derivative, we use Patankar-trick [23] obtaining two-time-level approximations:
(21)
Δ
i
l
≡
u
i
-
1
l
-
2
u
i
l
+
1
+
u
i
+
1
l
h
2
,
1
≤
i
<
N
,
2
u
i
-
1
l
h
i
-
1
h
i
+
h
i
-
1
-
u
i
l
+
1
h
i
h
i
-
1
+
u
i
+
1
l
h
i
h
i
-
1
+
h
i
,
N
≤
i
≤
N
+
M
-
2
,
and h
i
=
x
i
+
1
-
x
i
>
0. Then the difference scheme for (7) has the following norm:
(22)
u
i
l
+
1
=
u
i
l
+
k
σ
^
2
2
x
i
2
Δ
i
l
+
k
J
^
i
l
,
1
≤
i
≤
N
+
M
-
2
.
Note that the first expression of (21) corresponds to spatial zone with uniform discretization, while the second expression of (21) is related to the nonuniform discretization. On the other hand, for the approximation of the integral part of (7), instead of using the trapezoidal rule like in [14, 16, 20], we use a composite four-point integration formula of open type because of the higher order approximation of this rule [28, pp. 92-93]. This higher accuracy comes out because the singularity points of the kernel are not nodes of the integration mesh due to the truncation (see (12)) and the open type nature of the quadrature formula. Thus the approximation of (13) corresponding to the nodes x
=
x
i and τ
=
τ
l is given by
(23)
J
^
i
l
=
5
h
24
∑
j
=
0
N
/
5
‍
11
u
5
j
+
1
l
g
i
,
5
j
+
1
+
u
5
j
+
2
l
g
i
,
5
j
+
2
+
u
5
j
+
3
l
g
i
,
5
j
+
3
+
11
u
5
j
+
4
l
g
i
,
5
j
+
4
+
5
δ
24
A
∑
j
=
N
/
5
(
N
+
M
)
/
5
-
1
‍
11
u
5
j
+
1
l
g
i
,
5
j
+
1
x
5
j
+
1
2
+
u
5
j
+
2
l
g
i
,
5
j
+
2
x
5
j
+
2
2
+
u
5
j
+
3
l
g
i
,
5
j
+
3
x
5
j
+
3
2
+
11
u
5
j
+
4
l
g
i
,
5
j
+
4
x
5
j
+
4
2
,
where g
i
,
j
=
g
(
x
i
,
x
j
). Consequently the corresponding difference scheme for PIDE given by (7) takes the following form:
(24)
β
^
i
u
i
l
+
1
=
α
^
i
u
i
-
1
l
+
u
i
l
+
γ
^
i
u
i
+
1
l
+
k
J
^
i
l
,
where
(25)
α
^
i
=
k
σ
^
2
x
i
2
2
h
2
,
2
≤
i
≤
N
-
1
,
k
σ
^
2
A
2
h
h
+
A
δ
/
(
1
-
δ
)
,
i
=
N
,
k
σ
^
2
x
i
2
h
i
-
1
(
h
i
+
h
i
-
1
)
,
N
+
1
≤
i
≤
N
+
M
-
2
,
0
,
i
=
N
+
M
-
1
,
β
^
i
=
1
+
k
σ
^
2
x
i
2
h
2
,
1
≤
i
≤
N
-
1
,
1
+
k
σ
^
2
A
(
1
-
δ
)
h
δ
,
i
=
N
,
1
+
k
σ
^
2
x
i
2
h
i
h
i
-
1
,
N
+
1
≤
i
≤
N
+
M
-
2
,
1
,
i
=
N
+
M
-
1
,
γ
^
i
=
k
σ
^
2
x
i
2
2
h
2
,
1
≤
i
≤
N
-
1
,
k
σ
^
2
A
(
1
-
δ
)
δ
(
h
+
A
δ
/
(
1
-
δ
)
)
,
i
=
N
,
k
σ
^
2
x
i
2
h
i
(
h
i
+
h
i
-
1
)
,
N
+
1
≤
i
≤
N
+
M
-
2
.
From (24)-(25), one gets
(26)
u
i
l
+
1
=
α
i
u
i
-
1
l
+
β
i
u
i
l
+
γ
i
u
i
+
1
l
+
k
β
^
i
J
^
i
l
,
where
(27)
α
i
=
α
^
i
β
^
i
,
β
i
=
1
β
^
i
,
γ
i
=
γ
^
i
β
^
i
.In order to complete the difference scheme, we include the initial and boundary conditions. From (8), we have
(28)
u
i
0
=
max
⁡
(
x
i
-
E
,
0
)
,
1
≤
i
≤
N
+
M
-
1
.
On the other hand, for a vanilla call option the boundary condition for i
=
0 is
(29)
u
0
l
=
0
,
0
≤
l
≤
L
,
and by assuming the linear behavior of the solution for large values of the spatial variable, we have ∂
2
U
/
∂
x
2
→
0 and thus Δ
N
+
M
-
1
l
=
0 and the null integral term approximation J
N
+
M
-
1
=
0, for all time levels l. Thus from (26) for i
=
N
+
M
-
1, one gets
(30)
u
N
+
M
-
1
l
+
1
=
u
N
+
M
-
1
l
=
u
N
+
M
-
1
0
,
0
≤
l
≤
L
-
1
.
For the sake of convenience in the study of stability, we now introduce the vector formulation of the scheme (26)–(30) for the sake of the stability. Let us denote the vector in R
N
+
M
-
1 as U
l
=
u
1
l
u
2
l
⋯
u
N
+
M
-
1
l
t, and let P
∈
R
(
N
+
M
-
1
)
×
(
N
+
M
-
1
) be the tridiagonal matrix related to the differential part and defined by
(31)
P
=
β
1
γ
1
0
0
⋯
0
α
2
β
2
γ
2
0
⋯
0
0
α
3
β
3
γ
3
⋯
0
⋱
⋱
⋱
⋱
⋱
⋱
0
⋯
α
N
+
M
-
2
β
N
+
M
-
2
γ
N
+
M
-
2
0
⋯
0
α
N
+
M
-
1
β
N
+
M
-
1
.
Let B
=
(
b
i
j
) be the matrix in ∈
R
(
N
+
M
-
1
)
×
(
N
+
M
-
1
) related to the integral part whose entries b
i
j for each fixed i in 1
≤
i
≤
N
+
M
-
2 are defined by b
i
j
=
(
k
/
β
^
i
)
b
^
i
j, where
(32)
b
^
i
j
=
0
,
j
=
5
,
10
,
15,20
,
…
,
N
+
M
-
5
,
55
h
24
g
i
j
,
j
=
1,6
,
11
,
…
,
N
-
4
,
55
δ
24
A
x
j
2
g
i
j
,
j
=
N
+
1
,
N
+
6
,
…
,
M
+
N
-
4
,
5
h
24
g
i
j
,
j
=
2,7
,
12
,
…
,
N
-
3
,
5
δ
24
A
x
j
2
g
i
j
,
j
=
N
+
2
,
N
+
7
,
…
,
M
+
N
-
3
,
5
h
24
g
i
j
,
j
=
3,8
,
13
,
…
,
N
-
2
,
5
δ
24
A
x
j
2
g
i
j
,
j
=
N
+
3
,
N
+
8
,
…
,
M
+
N
-
2
,
55
h
24
g
i
j
,
j
=
4,9
,
14
,
…
,
N
-
1
,
55
δ
24
A
x
j
2
g
i
j
,
j
=
N
+
4
,
N
+
9
,
…
,
M
+
N
-
1
.
Hence scheme (26)–(30) can be written in the form
(33)
U
l
+
1
=
(
P
+
B
)
U
l
=
P
+
B
l
U
0
,
0
≤
l
≤
L
-
1
,
U
0
=
f
x
1
f
x
2
⋯
f
x
N
+
M
-
1
t
.


## 4. Unconditional Positivity and Stability

The numerical solution{
u
i
l
} of scheme (26) is unconditionally nonnegative because all coefficients of (26) and the initial and boundary conditions (28)–(30) are nonnegative.For the sake of clarity, before studying stability, let us recall the definition of the infinite norm for vectors and matrices. For a given vectorv
∈
R
n such that v
=
(
v
1
,
v
2
,
…
,
v
n
)
T, the infinite norm of v is denoted by v
∞ and is defined as v
∞
=
max
⁡
{
v
j
,
1
≤
j
≤
n
}. For a matrix A
=
(
a
i
j
)
∈
R
m
×
n, the infinite norm of A is denoted by A
∞ and is defined to be A
∞
=
max
⁡
{
∑
j
=
1
n
‍
|
a
i
j
|
,
1
≤
i
≤
m
}.In order to study the stability of the scheme given by (26)–(30), we first calculate the norm of the matrices P and B. Since the norm of the matrix P is obtained by
(34)
P
∞
=
max
⁡
i
∑
j
=
1
N
+
M
-
1
‍
P
i
j
=
max
⁡
i
∑
j
=
1
N
+
M
-
1
‍
P
i
j
,
∑
j
=
1
N
+
M
-
1
‍
P
1
j
=
β
1
+
γ
1
<
1 for i
=
1, and ∑
j
=
1
N
+
M
-
1
‍
P
i
j
=
α
i
+
β
i
+
γ
i
=
1 for 2
≤
i
≤
N
+
M
-
1, it follows that P
∞
=
1. On the other hand,
(35)
B
∞
=
max
⁡
i
∑
j
=
1
N
+
M
-
1
‍
b
i
j
=
max
⁡
i
∑
j
=
1
N
+
M
-
1
‍
b
i
j
,
1
≤
i
≤
N
+
M
-
1
.
Let m
~ be the row containing the maximum of (35), as β
^
m
~
≥
1; then we have
(36)
B
∞
=
∑
j
=
1
N
+
M
-
1
‍
b
m
~
j
=
k
β
^
m
~
∑
j
=
1
N
+
M
-
1
‍
b
^
m
~
j
≤
k
∑
j
=
1
N
+
M
-
1
‍
b
^
m
~
j
.
Note that the change of variables y
=
x
m
~
e
ϕ in (6) gives
(37)
λ
(
ε
)
=
∫
0
x
m
~
e
-
ε
‍
g
(
x
m
~
,
ϕ
)
d
ϕ
+
∫
x
m
~
e
ε
∞
‍
g
(
x
m
~
,
ϕ
)
d
ϕ
,
which coincides with (10) when U
(
ϕ
,
τ
)
=
1. Hence from (10), (35), and (37), we conclude that ∑
j
=
1
N
+
M
-
1
‍
b
^
m
~
j is an approximation for λ
(
ε
). Thus, for small enough h and δ, one gets [28]
(38)
∑
j
=
1
N
+
M
-
1
‍
b
^
m
~
j
<
λ
(
ε
)
+
1
.
Hence, from (35) and (37)(39)
B
∞
<
k
(
λ
(
ε
)
+
1
)
,
independently of the size of matrix B.Since there are many definitions of stability, let us point out that the chosen concept of stability is the one used in [15, pag. 7]. Note that scheme (26)–(30) is equivalent to the vector form scheme (33). After taking norms in (33),
(40)
U
l
+
1
∞
≤
P
∞
+
B
∞
U
l
∞
<
1
+
k
λ
ε
+
1
U
l
∞
.
Hence, from (40) and the fact that 0
≤
l
≤
L, k
L
=
T,
(41)
U
l
∞
U
0
∞
<
1
+
k
λ
ε
+
1
l
≤
exp
⁡
⁡
(
l
k
λ
ε
+
1
)
≤
exp
⁡
T
λ
ε
+
1
.
Thus the following result has been established.Theorem 1.
The numerical scheme given by (33) is strongly uniform stable.


## 5. Conditional Consistency of the Scheme

We say that a numerical scheme is consistent with a PIDE if the exact solution of the PIDE approximates well the difference scheme as the stepsizes discretization tend to zero [29, 30].Let us write (22) in the form
(42)
F
i
l
(
u
)
=
u
i
l
+
1
-
u
i
l
k
-
σ
^
2
2
x
i
2
Δ
i
l
-
J
^
i
l
=
0
,
let U
i
l
=
U
(
x
i
,
τ
l
) be the value of the exact price of (7) at (
x
i
,
τ
l
), let A
>
0 such that x
i
<
A
e
-
ε, and let us write the PIDE (7) as
(43)
L
(
U
)
=
J
(
U
)
,
where J
(
U
) is given by (11) and
(44)
L
(
U
)
=
∂
U
∂
τ
-
σ
^
2
2
x
2
∂
2
U
∂
x
2
.
Now we show that the local truncated error T
i
l
(
U
) at (
x
i
,
τ
l
), given by
(45)
T
i
l
(
U
)
=
F
i
l
(
U
)
-
L
(
U
i
l
)
+
J
(
U
i
l
)
,
satisfies
(46)
T
i
l
(
U
)
⟶
0
,
as
h
⟶
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 (
x
i
,
τ
l
), it follows that
(47)
U
i
l
+
1
-
U
i
l
k
=
∂
U
∂
τ
(
x
i
,
τ
l
)
+
k
E
i
l
(
1
)
,
E
i
l
(
1
)
=
1
2
∂
2
U
∂
τ
2
(
x
i
,
ζ
)
,
τ
l
<
ζ
<
τ
l
+
1
,
1
h
2
U
i
-
1
l
-
2
U
i
l
+
1
+
U
i
+
1
l
=
∂
2
U
∂
x
2
x
i
,
τ
l
+
h
2
E
i
l
(
2
)
-
k
h
2
E
i
l
(
3
)
,
where
(48)
E
i
l
(
2
)
=
1
12
∂
4
U
∂
x
4
(
ζ
~
,
τ
l
)
,
x
i
-
h
<
ζ
~
<
x
i
+
h
,
E
i
l
(
3
)
=
2
∂
U
∂
τ
(
x
i
,
ζ
)
,
τ
l
<
ζ
<
τ
l
+
1
.
Let us denote the maximum of the associated errors E
i
l
(
1
), E
i
l
(
2
), and E
i
l
(
3
) by W
i
l
(
1
), W
i
l
(
2
), and W
i
l
(
3
), respectively, where
(49)
E
i
l
1
≤
1
2
W
i
l
1
=
1
2
max
⁡
∂
2
U
∂
τ
2
x
i
,
ζ
;
τ
l
≤
ζ
≤
τ
l
+
1
,
(50)
E
i
l
2
≤
1
12
W
i
l
2
=
1
12
max
⁡
∂
4
U
∂
x
4
ζ
~
,
τ
l
;
x
i
-
h
≤
ζ
~
≤
x
i
+
h
,
(51)
E
i
l
3
≤
2
W
i
l
3
=
2
max
⁡
∂
U
∂
τ
(
x
i
,
ζ
)
;
τ
l
≤
ζ
≤
τ
l
+
1
.
To study the consistency of the integral part, it is convenient to rewrite it in the following form:
(52)
J
=
J
(
U
)
=
J
1
(
U
)
+
J
2
(
U
)
,
where
(53)
J
1
(
U
)
=
∫
0
A
‍
g
(
x
,
ϕ
)
U
(
ϕ
,
τ
l
)
d
ϕ
,
J
2
(
U
)
=
A
∫
0
1
‍
g
(
x
,
A
z
)
U
(
A
z
,
τ
l
)
d
ϕ
.
In accordance with [29] let us denote the local consistency error of J
1
,
i
l:
(54)
C
1
,
i
l
(
h
,
A
)
=
∫
0
A
‍
g
(
x
,
ϕ
)
U
(
ϕ
,
τ
l
)
d
ϕ
-
J
^
1
,
i
l
(
[
0
,
A
]
)
,
where
(55)
J
^
1
,
i
l
0
,
A
=
5
h
24
∑
j
=
0
N
/
5
‍
11
u
5
j
+
1
l
g
i
,
5
j
+
1
+
u
5
j
+
2
l
g
i
,
5
j
+
2
+
u
5
j
+
3
l
g
i
,
5
j
+
3
+
11
u
5
j
+
4
l
g
i
,
5
j
+
4
.
By [28, p. 92] we have
(56)
C
1
,
i
l
h
,
A
≤
95
A
h
4
144
β
^
i
max
⁡
0
≤
x
≤
A
U
x
,
τ
l
g
x
,
x
i
4
=
95
A
h
4
144
W
i
l
(
4
)
,
where U
(
x
,
τ
l
)
g
(
x
,
x
i
)
(
4
) in (56) is the fourth derivative with respect to the variable x.Similarly, the local consistency error for the unbounded region is given by(57)
C
2
,
i
l
(
δ
,
A
)
=
A
∫
0
1
‍
g
x
,
A
z
U
A
z
,
τ
l
d
z
z
2
-
J
^
2
,
i
l
0,1
,
where
(58)
J
^
2
,
i
l
0,1
=
5
δ
24
A
·
∑
j
=
N
/
5
(
N
+
M
)
/
5
-
1
‍
11
u
5
j
+
1
l
g
i
,
5
j
+
1
x
5
j
+
1
2
+
u
5
j
+
2
l
g
i
,
5
j
+
2
x
5
j
+
2
2
+
u
5
j
+
3
l
g
i
,
5
j
+
3
x
5
j
+
3
2
+
11
u
5
j
+
4
l
g
i
,
5
j
+
4
x
5
j
+
4
2
,
C
2
,
i
l
δ
,
A
≤
95
A
δ
4
144
β
^
i
max
⁡
0
≤
z
≤
1
1
z
2
U
A
z
,
τ
l
g
x
i
,
A
z
4
=
95
A
δ
4
144
W
i
l
(
5
)
.
Thus the local truncation error is given by
(59)
T
i
l
U
=
k
E
i
1
+
σ
^
2
x
i
2
2
h
2
E
i
l
2
-
k
h
2
E
i
l
3
+
C
1
,
i
l
(
h
,
A
)
+
C
2
,
i
l
(
δ
,
A
)
,
T
i
l
U
≤
k
W
i
l
1
+
σ
^
2
A
2
h
2
24
W
i
l
2
+
σ
^
2
A
2
k
h
2
W
i
l
(
3
)
+
95
A
144
(
h
4
W
i
l
(
4
)
+
δ
4
W
i
l
(
5
)
)
.
Consequently, the order of the local truncated error is given by
(60)
T
i
l
(
U
)
=
O
(
k
)
+
O
(
h
2
)
+
O
(
k
h
2
)
+
O
(
δ
4
)
,
and the following result has been established.Theorem 2.
Under conditionk
=
O
(
h
2
+
ϵ
), ϵ
>
0, the numerical scheme (42) is consistent with (7).


## 6. Numerical Examples

Based on the double discretization and Patankar-trick, a difference scheme has been established to obtain a numerical solution for the option price. This solution is guaranteed to be nonnegative and stable. In this section, we illustrate with several examples the behavior of the option price obtained by this scheme using MATLAB. All examples are done using CPU with Microprocessor 3.4 GHz Intel Core i7.The following example illustrates that the consistency conditionk
=
O
(
h
2
+
ϵ
) cannot be ignored.Example 1.
Here in this example the parameters have been selected as follows:T
=
1, E
=
10, A
=
3
E, σ
=
0.25, r
=
0.01, q
=
0, C
=
1, G
=
25, M
=
25, Y
=
1.65, ε
=
0.15, h
=
0.25, and δ
=
0.1, for several values of k such that k
=
h
2.5, h
1.5, and h. Figure 1 shows that the consistency condition holds for k
=
h
2.5, while, for the other two values, it is broken and the values of the option price become unreliable.Figure 1
The effect of consistency condition onV.The aim of the following examples is to exhibit the effects of different parameters such ash, k, and ε on the variation of the absolute error in two cases: first, when Y
=
0 (Variance Gamma case) and, second, for CGMY process when Y
=
1.5,1.98. Also the CPU time is given in seconds (sec.).In the next example we calculate the associated error with this numerical scheme for the Variance Gamma model as a particular case (Y
=
0) of CGMY model for which the exact solution is known [31].Example 2.
ForY
=
0 parameters have been selected as follows: T
=
0.5, E
=
80, A
=
3
E, S
=
100, σ
=
0.25, r
=
0.1, q
=
0, C
=
1, G
=
30, M
=
20, and δ
=
0.1; Table 1 shows the variation of the absolute error with h with fixed k
=
0.003 and for two values of ε
=
0.1 and 0.15. From Table 1, it is observed that the associated error exhibits the second order convergence rate α providing that k
/
h
2 is small enough in all the cases.
Table2 reveals the change of the associated error for various values of time stepsize k, while h
=
0.8 for ε
=
0.1 and 0.15. Notice that the associated error due to the change of k satisfies the expected first order convergence rate β.
The aim of Table3 is to show the sensitivity of the associated error of the option price due to the variation of ε, for h
=
0.5 and 0.35, while k
=
0.005.Table 1
Errors and convergence rates due to the change ofh for VG.
h
ε
=
0.1
ε
=
0.15
Absolute error
α
CPU time
Absolute error
α
CPU time
(sec.)
(sec.)
1
.
2
6.835
×
1
0
-
4
—
0.19
6.15
×
1
0
-
4
—
0.19
1
4.821
×
1
0
-
4
1.915
0.26
4.33
×
1
0
-
4
1.925
0.26
0.8
3.138
×
1
0
-
4
1.924
0.38
2.81
×
1
0
-
4
1.938
0.38
0.5
1.266
×
1
0
-
4
1.931
0.44
1.124
×
1
0
-
4
1.95
0.44Table 2
Errors and convergence rates due to the change ofk for VG.
k
ε
=
0.1
ε
=
0.15
Absolute error
β
CPU time
Absolute error
β
CPU time
(sec.)
(sec.)
0
.
1
7.654
×
1
0
-
3
—
0.248
5.321
×
1
0
-
3
—
0.248
0.05
3.962
×
1
0
-
3
0.950
0.256
2.793
×
1
0
-
3
0.930
0.256
0.025
2.041
×
1
0
-
3
0.957
0.263
1.429
×
1
0
-
3
0.967
0.263
0.01
8.367
×
1
0
-
4
0.973
0.271
5.794
×
1
0
-
4
0.985
0.271Table 3
The associated errors for several values ofε for VG.
ε
h
=
0.5
h
=
0.35
Absolute error
CPU time
Absolute error
CPU time
(sec.)
(sec.)
0
.
75
3.495
×
1
0
-
4
0.45
1.473
×
1
0
-
4
0.78
0.5
7.643
×
1
0
-
4
0.45
4.587
×
1
0
-
4
0.78
0.25
5.874
×
1
0
-
4
0.45
3.198
×
1
0
-
4
0.78
0.1
2.382
×
1
0
-
4
0.45
1.258
×
1
0
-
4
0.78Example 3.
Here we compare in Tables4 and 5 in our results with the reference values given in [10, Tables 9, 10] related to accuracy and computational time. We consider the CGMY model for the following parameters: T
=
1, r
=
0.1, q
=
0, C
=
1, G
=
M
=
5, E
=
100, S
=
100, A
=
3
E, δ
=
0.1, and k
=
0.003. Table 4 shows the variation of the associated error for several values of h when Y
=
1.5 and 1.98, while ε
=
0.1.
The variation of the associated error for several values ofε is presented in Table 5 for Y
=
1.5 and 1.98, while h
=
1.Table 4
Comparison of errors and convergence rates due to the change ofh for CGMY model.
h
Y
=
1.5
Y
=
1.98
Absolute error
α
CPU time
Absolute error
α
CPU time
(sec.)
(sec.)
1
.
5
6.1
×
1
0
-
4
—
0.13
6.62
×
1
0
-
4
—
0.13
1.2
3.98
×
1
0
-
4
1.91
0.16
4.32
×
1
0
-
4
1.913
0.16
1
2.8
×
1
0
-
4
1.928
0.2
3.04
×
1
0
-
4
1.927
0.2
0.8
1.8
×
1
0
-
4
1.955
0.25
1.97
×
1
0
-
4
1.944
0.25
0.5
7.18
×
1
0
-
5
1.967
0.47
7.83
×
1
0
-
5
1.962
0.47Table 5
Comparison of errors due to the variation ofε for CGMY model.
ε
Y
=
1.5
Y
=
1.98
Absolute error
CPU time
Absolute error
CPU time
(sec.)
(sec.)
0
.
8
2.52
×
1
0
-
4
0.2
2.63
×
1
0
-
4
0.2
0.4
9.26
×
1
0
-
4
0.2
7.26
×
1
0
-
4
0.2
0.2
5.75
×
1
0
-
4
0.2
4.39
×
1
0
-
4
0.2
0.1
2.8
×
1
0
-
4
0.2
3.04
×
1
0
-
4
0.2The next example reveals that the double discretization strategy reduces the error near the parameterA by changing the stepsize δ.Example 4.
Consider a call option under the Variance Gamma process with parametersC
=
1, G
=
20  M
=
30, T
=
1, r
=
0.01, q
=
0, σ
=
0.2, ε
=
0.12, E
=
10, A
=
3
E, k
=
0.005, and h
=
0.35. Figure 2 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 2
The associated error for several values ofδ.Example 5.
Figures3(a) and 3(b) describe the behavior of the Greek parameters Delta and Gamma for European call option. They exhibit the Greek parameters as functions in the underlying asset S and time t. The parameters have been chosen as follows: T
=
1, E
=
10, A
=
3
E, σ
=
0.25, r
=
0.01, q
=
0, C
=
1, G
=
35, M
=
35, Y
=
1.6, ε
=
0.15, h
=
0.25, δ
=
0.1, and k
=
0.04.The Greek parameters for European call option.
(a)
Δ as a function of S and t
(b)
Γ as a function of S and t


## 7. Conclusion

In this paper we develop a double discretization technique for the numerical solution of option pricing PIDE model under CGMY process, previously used in [15]. This technique allows us to consider the information of the solution of the integral part outside of the bounded numerical domain. In this paper we propose a finite difference scheme that introduces two main advantages with respect to [15] such as the unconditional positivity and stability of the proposed numerical scheme as well as a higher accuracy due to the use of more accurate numerical integration.



---

*Source: 960728-2015-02-04.xml*