# An Unconditionally Stable, Positivity-Preserving Splitting Scheme for Nonlinear Black-Scholes Equation with Transaction Costs

Authors: Jianqiang Guo; Wansheng Wang
Journal: The Scientific World Journal
 (2014)

---

## Abstract

This paper deals with the numerical analysis of nonlinear Black-Scholes equation with transaction costs. An unconditionally stable and monotone splitting method, ensuring positive numerical solution and avoiding unstable oscillations, is proposed. This numerical method is based on the LOD-Backward Euler method which allows us to solve the discrete equation explicitly. The numerical results for vanilla call option and for European butterfly spread are provided. It turns out that the proposed scheme is efficient and reliable.


---

## Body

## 1. Introduction

One of the modern financial theory’s biggest successes in terms of both approach and applicability has been the Black-Scholes option pricing model developed by Black and Scholes in 1973 [1] and previously by Merton [2]:
(1)
∂
V
∂
τ
+
1
2
σ
0
2
s
2
∂
2
V
∂
s
2
+
r
s
∂
V
∂
s
-
r
V
=
0
,
hhhhhhh
(
s
,
τ
)
∈
Ω
T
∶
=
Ω
×
(
0
,
T
]
,
where r
>
0 and σ
0 are given real constants that represent the interest rate and the volatility, respectively, s is the price of the underlying asset, and T is the maturity date. The celebrated Black-Scholes model is based on several restrictive assumptions such as liquid, frictionless, and complete markets. In recent years, nonlinear Black-Scholes models have been used to build transaction costs, market liquidity, or volatility uncertainty into the celebrated Black-Scholes concept.In this paper, we are interested in the option pricing model with transaction costs proposed by Barles and Soner [3] that are motivated by Hodges and Neuberger [4]. In practice, transaction costs arise when trading securities. Recent studies of their influence reveal that they result in a nonnegligible increase in the option price, although they are generally small for institutional investors. To show this increase, in Barles and Soner’s model, the constant volatility σ
0 in the linear model is replaced by the nonlinear volatility σ that reads
(2)
σ
2
=
σ
0
2
(
1
+
Ψ
(
e
r
(
T
-
τ
)
a
2
s
2
∂
2
V
∂
s
2
)
)
,
where a
=
μ
γ
P with the proportional transaction cost being μ, the risk aversion factor being γ, and the number of options to be sold being P, and Ψ
(
x
) is the solution of the following ordinary differential equation:
(3)
Ψ
′
(
x
)
=
Ψ
(
x
)
+
1
2
x
Ψ
(
x
)
-
x
,
x
≠
0
,
Ψ
(
0
)
=
0
.
Barles and Soner’s option pricing model now reads
(4)
∂
V
∂
τ
+
1
2
σ
2
s
2
∂
2
V
∂
s
2
+
r
s
∂
V
∂
s
-
r
V
=
0
,
(
s
,
τ
)
∈
Ω
T
∶
=
Ω
×
(
0
,
T
]
,
with the terminal condition
(5)
V
(
s
,
T
)
=
f
(
s
)
,
s
∈
Ω
∶
=
(
0
,
+
∞
)
.
The payoff function f
(
s
) is assumed to be a continuous piecewise linear function.Many results have been reported for the numerical solutions of linear Black-Scholes equations (see, e.g., [5–9]). On the other hand, because of the nonlinear nature of this model, numerical methods are mandatory to price derivatives and portfolios. The strong nonlinearity of problem (4) makes it difficult to obtain the reliable numerical solutions. Implicit numerical schemes have been used for numerically solving nonlinear option pricing PDEs [10], and an iterative approach is required to solve the nonlinear algebraic equation resulting from the discretization, which results in more computational cost. Consequently, some high-order compact semi-implicit difference schemes are proposed for numerically solving the nonlinear option pricing model [11, 12]. Some researchers (see, e.g., [3, 13–15]) have also constructed explicit finite difference schemes for (4)-(5) and investigated their consistency and stability. However, these explicit schemes have the disadvantage that strictly restrictive conditions on the discretization parameters are needed to guarantee stability and positivity. To relax the restrictive conditions, in [16], Zhou et al. proposed an unconditionally stable explicit finite difference scheme based on a nonstandard approximation of the second partial derivative. However, this scheme is conditionally consistent, and the truncation error depends on the ratio of the time stepsize and the square of the space stepsize.In this paper, we will consider a splitting method with inhomogeneous boundary conditions. This method proposed here is unconditionally stable, monotone, and positivity preserving. It is also consistent and essentially a “limit” version of the LOD-Backward Euler method (see Chapter IV on splitting methods in [17]) and therefore allows us to solve the discrete equation explicitly.The paper is organized as follows: we begin by transforming the original equations into nonlinear heat equations and considering the spatial semidiscretization. The splitting scheme will be discussed in Section3 after linearization of semidiscrete system. The stability, the monotonicity, the positivity-preserving property, and the convergence of this scheme are analysed in Section 4. To illustrate our method, we present some numerical experiments in Section 5. Finally, we give a summary.


## 2. Transformation and Spatial Semidiscretization

### 2.1. Preliminaries

We first consider the properties of the functionΨ appearing in (2), which will play an important role in the following numerical analysis.Lemma 1 (see [14, 15]).
The solutionΨ of ODE (3) exists and is unique, and it satisfies the following: (1)
Ψ is an increasing function mapping the real line onto the interval (
-
1
,
+
∞
);
(2)
Ψ
=
Ψ
(
x
) is implicitly defined by
(6)
x
=
(
-
arcsinh
Ψ
Ψ
+
1
+
Ψ
)
2
,
i
f
Ψ
>
0
,
x
=
-
(
arcsinh
-
Ψ
Ψ
+
1
-
-
Ψ
)
2
,
i
f
-
1
<
Ψ
<
0
;
(3)
ifx
>
0, then the function Ψ
(
x
) is bounded and
(7)
0
<
Ψ
(
x
)
<
Ψ
′
(
x
2
)
x
+
d
2
,
where
(8)
x
2
=
(
sinh
⁡
2
-
2
(
sinh
⁡
2
)
2
+
1
)
2
≈
9.58
,
Ψ
(
x
2
)
=
(
sinh
⁡
2
)
2
,
(9)
Ψ
′
(
x
2
)
=
(
e
8
+
2
e
4
+
1
)
2
e
16
-
66
e
8
+
1
≈
1.10
,
d
2
=
Ψ
(
x
2
)
-
Ψ
′
(
x
2
)
≈
2.62
.Lemma 2 (see [14, 15]).
The functiong
(
x
)
=
x
Ψ
(
x
) is continuously differentiable at x
=
0 and satisfies
(10)
|
g
′
(
x
)
|
≤
max
⁡
{
G
,
2
|
x
|
Ψ
′
(
x
2
)
+
d
2
}
,
where x
2 and d
2 are given by (8) and (9), respectively, and
(11)
G
=
max
⁡
{
|
g
′
(
x
)
|
;
x
1
≤
x
≤
x
2
}
,
x
1
=
-
(
4
π
-
3
3
)
2
36
.Combining Lemma2 and the results obtained in [3] leads to the following lemma.Lemma 3.
The functionw
(
x
)
=
x
+
x
Ψ
(
x
) is a nondecreasing function of x and continuously differentiable at x
=
0. This implies that
(12)
1
+
Ψ
(
x
)
+
x
Ψ
′
(
x
)
=
d
d
x
[
x
+
x
Ψ
(
x
)
]
≥
0
.


### 2.2.  The Transformed Problem

After considering the change of variable(13)
S
=
e
r
(
T
-
τ
)
s
,
t
=
T
-
τ
,
U
(
S
,
t
)
=
e
r
(
T
-
τ
)
V
(
s
,
τ
)
,
we transform the problem (4)-(5) into the following initial value problem
(14)
∂
U
∂
t
-
σ
0
2
S
2
2
(
1
+
Ψ
(
a
2
S
2
∂
2
U
∂
S
2
)
)
∂
2
U
∂
S
2
=
0
,
mmmmmfmfmm
(
S
,
t
)
∈
Ω
×
(
0
,
T
]
,
(15)
U
(
S
,
0
)
=
f
(
S
)
,
S
∈
Ω
.
The boundary conditions take the form of
(16)
U
(
0
,
t
)
=
0
,
U
(
S
,
t
)
=
α
(
+
∞
,
t
)
,
t
∈
[
0
,
T
]
,
where the second boundary condition in (16) is derived by asymptotic considerations in (
s
,
τ
) coordinates for extreme value s
→
∞ and subsequent transformation. For example, the boundary conditions for a vanilla call are given by
(17)
U
(
0
,
t
)
=
0
,
lim
⁡
S
→
+
∞
α
(
S
,
t
)
S
=
1
,
t
∈
[
0
,
T
]
,
and the boundary conditions for a butterfly spread are given by
(18)
U
(
0
,
t
)
=
0
,
lim
⁡
S
→
+
∞
α
(
S
,
t
)
=
0
,
t
∈
[
0
,
T
]
.


### 2.3. The Semidiscrete Nonlinear System

To numerically approximate the solution of (14)–(16), we should consider a bounded numerical domain (
S
,
t
)
∈
[
0
,
b
]
×
[
0
,
T
]. Then, the problem (14)-(15) is equipped with the boundary conditions
(19)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
,
t
)
,
t
∈
[
0
,
T
]
.
Note that the second equality in (19) derives from the second equality in (16). Taking these into consideration, the boundary conditions for a vanilla call are given by
(20)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
)
=
max
⁡
{
0
,
b
-
E
}
,
m
f
m
m
t
∈
[
0
,
T
]
,
where E is the strike price, and the boundary conditions for a butterfly spread are given by
(21)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
)
=
f
(
b
)
,
t
∈
[
0
,
T
]
.
We introduce the spatial grid Ω
h with step h by the nodes S
i
=
i
h, i
=
0,1
,
…
,
M, so that M
h
=
b. After performing the second-order central finite difference approximation of the partial derivative (
∂
2
U
/
∂
S
2
)
(
S
,
t
) as
(22)
∂
2
U
∂
S
2
=
U
(
S
i
+
1
,
t
)
-
2
U
(
S
i
,
t
)
+
U
(
S
i
-
1
,
t
)
h
2
+
O
(
h
2
)
=
Δ
i
U
(
t
)
+
O
(
h
2
)
,
we obtained the corresponding ODEs system for the semidiscrete solution u
(
t
)
=
[
U
1
(
t
)
,
U
2
(
t
)
,
…
,
U
M
-
1
(
t
)
]
T as
(23)
u
′
(
t
)
=
A
(
u
(
t
)
)
u
(
t
)
+
g
(
u
(
t
)
)
,
t
∈
[
0
,
T
]
,
with
(24)
A
(
u
)
=
1
h
2
tridiag
(
β
i
(
u
)
,
-
2
β
i
(
u
)
,
β
i
(
u
)
)
,
β
i
(
u
)
=
σ
i
2
(
u
)
S
i
2
,
σ
i
2
(
u
)
=
1
2
σ
0
2
(
1
+
Ψ
i
)
,
Ψ
i
=
Ψ
(
a
2
S
i
2
Δ
i
u
(
t
)
)
,
where g
∈
ℝ
M
-
1 is a vector, generated by the boundary conditions,
(25)
g
(
u
(
t
)
)
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
α
(
b
)
]
T
.Obviously, the semidiscrete difference scheme (23) is consistent. For the stability, we have the following theorem.Theorem 4.
For system (23), the following maximum norm contractivity holds:
(26)
∥
u
(
t
)
-
u
~
(
t
)
∥
∞
≤
∥
u
(
0
)
-
u
~
(
0
)
∥
∞
,
∀
t
∈
[
0
,
T
]
,
where u
~ is the solution of the perturbation problem
(27)
u
~
′
(
t
)
=
A
(
u
~
(
t
)
)
u
~
(
t
)
+
g
(
u
~
(
t
)
)
,
t
∈
[
0
,
T
]
,
with the initial data u
~
(
0
).Proof.
To obtain the stability of the semidiscrete difference scheme (23), we first calculate
(28)
∂
β
i
∂
U
i
=
1
2
σ
0
2
S
i
2
Ψ
i
′
·
(
-
2
a
2
S
i
2
h
2
)
=
-
1
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
;
∂
β
i
∂
U
i
-
1
=
∂
β
i
∂
U
i
+
1
=
1
2
σ
0
2
S
i
2
Ψ
i
′
·
(
a
2
S
i
2
h
2
)
=
1
2
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
.
Now, let nonlinear operator F
(
u
) be defined by F
(
u
)
=
A
(
u
)
u
+
g
(
u
). For system (23), we have
(29)
∂
F
i
∂
U
i
=
1
h
2
(
-
2
β
i
-
2
∂
β
i
∂
U
i
U
i
+
∂
β
i
∂
U
i
(
U
i
-
1
+
U
i
+
1
)
)
=
1
h
2
(
-
2
β
i
-
U
i
-
1
-
2
U
i
+
U
i
+
1
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
)
=
-
1
h
2
(
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
)
,
∂
F
i
∂
U
i
-
1
=
∂
F
i
∂
U
i
+
1
=
1
h
2
(
β
i
+
U
i
-
1
-
2
U
i
+
U
i
+
1
h
2
1
2
σ
0
2
a
2
S
i
4
Ψ
i
′
)
=
1
2
h
2
(
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
)
.
Since
(30)
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
=
1
+
Ψ
i
(
a
2
S
i
2
Δ
i
u
(
t
)
)
+
a
2
S
i
2
Δ
i
u
(
t
)
Ψ
i
′
(
a
2
S
i
2
Δ
i
u
(
t
)
)
,
by Lemma 3, we obtain
(31)
∂
F
i
∂
U
i
≤
0
,
∂
F
i
∂
U
i
-
1
=
∂
F
i
∂
U
i
+
1
≥
0
.
Then, inequality (26) follows from μ
∞
[
J
F
(
u
)
]
=
0 [18], where J
F
(
u
) is the Jacobian matrix of F and μ
∞
[
·
] is the logarithmic maximum norm.


## 2.1. Preliminaries

We first consider the properties of the functionΨ appearing in (2), which will play an important role in the following numerical analysis.Lemma 1 (see [14, 15]).
The solutionΨ of ODE (3) exists and is unique, and it satisfies the following: (1)
Ψ is an increasing function mapping the real line onto the interval (
-
1
,
+
∞
);
(2)
Ψ
=
Ψ
(
x
) is implicitly defined by
(6)
x
=
(
-
arcsinh
Ψ
Ψ
+
1
+
Ψ
)
2
,
i
f
Ψ
>
0
,
x
=
-
(
arcsinh
-
Ψ
Ψ
+
1
-
-
Ψ
)
2
,
i
f
-
1
<
Ψ
<
0
;
(3)
ifx
>
0, then the function Ψ
(
x
) is bounded and
(7)
0
<
Ψ
(
x
)
<
Ψ
′
(
x
2
)
x
+
d
2
,
where
(8)
x
2
=
(
sinh
⁡
2
-
2
(
sinh
⁡
2
)
2
+
1
)
2
≈
9.58
,
Ψ
(
x
2
)
=
(
sinh
⁡
2
)
2
,
(9)
Ψ
′
(
x
2
)
=
(
e
8
+
2
e
4
+
1
)
2
e
16
-
66
e
8
+
1
≈
1.10
,
d
2
=
Ψ
(
x
2
)
-
Ψ
′
(
x
2
)
≈
2.62
.Lemma 2 (see [14, 15]).
The functiong
(
x
)
=
x
Ψ
(
x
) is continuously differentiable at x
=
0 and satisfies
(10)
|
g
′
(
x
)
|
≤
max
⁡
{
G
,
2
|
x
|
Ψ
′
(
x
2
)
+
d
2
}
,
where x
2 and d
2 are given by (8) and (9), respectively, and
(11)
G
=
max
⁡
{
|
g
′
(
x
)
|
;
x
1
≤
x
≤
x
2
}
,
x
1
=
-
(
4
π
-
3
3
)
2
36
.Combining Lemma2 and the results obtained in [3] leads to the following lemma.Lemma 3.
The functionw
(
x
)
=
x
+
x
Ψ
(
x
) is a nondecreasing function of x and continuously differentiable at x
=
0. This implies that
(12)
1
+
Ψ
(
x
)
+
x
Ψ
′
(
x
)
=
d
d
x
[
x
+
x
Ψ
(
x
)
]
≥
0
.


## 2.2.  The Transformed Problem

After considering the change of variable(13)
S
=
e
r
(
T
-
τ
)
s
,
t
=
T
-
τ
,
U
(
S
,
t
)
=
e
r
(
T
-
τ
)
V
(
s
,
τ
)
,
we transform the problem (4)-(5) into the following initial value problem
(14)
∂
U
∂
t
-
σ
0
2
S
2
2
(
1
+
Ψ
(
a
2
S
2
∂
2
U
∂
S
2
)
)
∂
2
U
∂
S
2
=
0
,
mmmmmfmfmm
(
S
,
t
)
∈
Ω
×
(
0
,
T
]
,
(15)
U
(
S
,
0
)
=
f
(
S
)
,
S
∈
Ω
.
The boundary conditions take the form of
(16)
U
(
0
,
t
)
=
0
,
U
(
S
,
t
)
=
α
(
+
∞
,
t
)
,
t
∈
[
0
,
T
]
,
where the second boundary condition in (16) is derived by asymptotic considerations in (
s
,
τ
) coordinates for extreme value s
→
∞ and subsequent transformation. For example, the boundary conditions for a vanilla call are given by
(17)
U
(
0
,
t
)
=
0
,
lim
⁡
S
→
+
∞
α
(
S
,
t
)
S
=
1
,
t
∈
[
0
,
T
]
,
and the boundary conditions for a butterfly spread are given by
(18)
U
(
0
,
t
)
=
0
,
lim
⁡
S
→
+
∞
α
(
S
,
t
)
=
0
,
t
∈
[
0
,
T
]
.


## 2.3. The Semidiscrete Nonlinear System

To numerically approximate the solution of (14)–(16), we should consider a bounded numerical domain (
S
,
t
)
∈
[
0
,
b
]
×
[
0
,
T
]. Then, the problem (14)-(15) is equipped with the boundary conditions
(19)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
,
t
)
,
t
∈
[
0
,
T
]
.
Note that the second equality in (19) derives from the second equality in (16). Taking these into consideration, the boundary conditions for a vanilla call are given by
(20)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
)
=
max
⁡
{
0
,
b
-
E
}
,
m
f
m
m
t
∈
[
0
,
T
]
,
where E is the strike price, and the boundary conditions for a butterfly spread are given by
(21)
U
(
0
,
t
)
=
0
,
U
(
b
,
t
)
=
α
(
b
)
=
f
(
b
)
,
t
∈
[
0
,
T
]
.
We introduce the spatial grid Ω
h with step h by the nodes S
i
=
i
h, i
=
0,1
,
…
,
M, so that M
h
=
b. After performing the second-order central finite difference approximation of the partial derivative (
∂
2
U
/
∂
S
2
)
(
S
,
t
) as
(22)
∂
2
U
∂
S
2
=
U
(
S
i
+
1
,
t
)
-
2
U
(
S
i
,
t
)
+
U
(
S
i
-
1
,
t
)
h
2
+
O
(
h
2
)
=
Δ
i
U
(
t
)
+
O
(
h
2
)
,
we obtained the corresponding ODEs system for the semidiscrete solution u
(
t
)
=
[
U
1
(
t
)
,
U
2
(
t
)
,
…
,
U
M
-
1
(
t
)
]
T as
(23)
u
′
(
t
)
=
A
(
u
(
t
)
)
u
(
t
)
+
g
(
u
(
t
)
)
,
t
∈
[
0
,
T
]
,
with
(24)
A
(
u
)
=
1
h
2
tridiag
(
β
i
(
u
)
,
-
2
β
i
(
u
)
,
β
i
(
u
)
)
,
β
i
(
u
)
=
σ
i
2
(
u
)
S
i
2
,
σ
i
2
(
u
)
=
1
2
σ
0
2
(
1
+
Ψ
i
)
,
Ψ
i
=
Ψ
(
a
2
S
i
2
Δ
i
u
(
t
)
)
,
where g
∈
ℝ
M
-
1 is a vector, generated by the boundary conditions,
(25)
g
(
u
(
t
)
)
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
α
(
b
)
]
T
.Obviously, the semidiscrete difference scheme (23) is consistent. For the stability, we have the following theorem.Theorem 4.
For system (23), the following maximum norm contractivity holds:
(26)
∥
u
(
t
)
-
u
~
(
t
)
∥
∞
≤
∥
u
(
0
)
-
u
~
(
0
)
∥
∞
,
∀
t
∈
[
0
,
T
]
,
where u
~ is the solution of the perturbation problem
(27)
u
~
′
(
t
)
=
A
(
u
~
(
t
)
)
u
~
(
t
)
+
g
(
u
~
(
t
)
)
,
t
∈
[
0
,
T
]
,
with the initial data u
~
(
0
).Proof.
To obtain the stability of the semidiscrete difference scheme (23), we first calculate
(28)
∂
β
i
∂
U
i
=
1
2
σ
0
2
S
i
2
Ψ
i
′
·
(
-
2
a
2
S
i
2
h
2
)
=
-
1
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
;
∂
β
i
∂
U
i
-
1
=
∂
β
i
∂
U
i
+
1
=
1
2
σ
0
2
S
i
2
Ψ
i
′
·
(
a
2
S
i
2
h
2
)
=
1
2
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
.
Now, let nonlinear operator F
(
u
) be defined by F
(
u
)
=
A
(
u
)
u
+
g
(
u
). For system (23), we have
(29)
∂
F
i
∂
U
i
=
1
h
2
(
-
2
β
i
-
2
∂
β
i
∂
U
i
U
i
+
∂
β
i
∂
U
i
(
U
i
-
1
+
U
i
+
1
)
)
=
1
h
2
(
-
2
β
i
-
U
i
-
1
-
2
U
i
+
U
i
+
1
h
2
σ
0
2
a
2
S
i
4
Ψ
i
′
)
=
-
1
h
2
(
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
)
,
∂
F
i
∂
U
i
-
1
=
∂
F
i
∂
U
i
+
1
=
1
h
2
(
β
i
+
U
i
-
1
-
2
U
i
+
U
i
+
1
h
2
1
2
σ
0
2
a
2
S
i
4
Ψ
i
′
)
=
1
2
h
2
(
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
)
.
Since
(30)
2
β
i
+
Δ
i
u
(
t
)
σ
0
2
a
2
S
i
4
Ψ
i
′
=
1
+
Ψ
i
(
a
2
S
i
2
Δ
i
u
(
t
)
)
+
a
2
S
i
2
Δ
i
u
(
t
)
Ψ
i
′
(
a
2
S
i
2
Δ
i
u
(
t
)
)
,
by Lemma 3, we obtain
(31)
∂
F
i
∂
U
i
≤
0
,
∂
F
i
∂
U
i
-
1
=
∂
F
i
∂
U
i
+
1
≥
0
.
Then, inequality (26) follows from μ
∞
[
J
F
(
u
)
]
=
0 [18], where J
F
(
u
) is the Jacobian matrix of F and μ
∞
[
·
] is the logarithmic maximum norm.


## 3. Splitting Time-Stepping Method

In this section, we will focus on the time integration methods of system (23). The implicit numerical schemes are generally viewed as stable methods, but [10] showed that the stability of the BDF schemes with orders 1 and 2 and the Crank-Nicolson scheme is still restricted by a condition. Additionally, for the implicit schemes, the nonlinear iteration will be required and will produce the additional computational cost in each time step. Fully explicit finite difference schemes for the PDEs (14)-(15) are also constructed in [3, 13–15]. However, these schemes are stable only for the severe restriction on the time stepsize Δ
t. For example, in [14], the condition reads
(32)
Δ
t
h
3
[
(
1
+
d
2
)
h
+
4
a
2
E
2
Ψ
′
(
x
2
)
]
≤
1
8
E
2
.
To relax the condition, in this paper, we will propose an unconditionally stable method, which allows us to solve the discrete equation explicitly, based on the LOD-Backward Euler method (see Chapter IV on splitting methods in [17]).


### 3.1. Linearization System

Let us sett
n
=
n
Δ
t, n
=
1,2
,
…
,
N, for the temporal stepsize Δ
t
=
T
/
N. To linearize the nonlinear system (23), we allow the nonlinearities in (23) to lag one step behind and obtain the following linear system:
(33)
u
′
(
t
)
=
A
(
u
(
t
n
)
)
u
(
t
)
+
g
(
u
(
t
n
)
)
=
A
n
u
(
t
)
+
g
n
,
m
t
∈
[
t
n
,
t
n
+
1
]
,
with
(34)
A
n
=
1
h
2
tridiag
(
β
i
n
,
-
2
β
i
n
,
β
i
n
)
,
g
n
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
n
α
(
b
)
]
T
,
where
(35)
β
i
n
=
σ
i
,
n
2
S
i
2
,
σ
i
,
n
2
=
1
2
σ
0
2
(
1
+
Ψ
i
n
)
,
Ψ
i
n
=
Ψ
(
a
2
S
i
2
Δ
i
U
n
)
,
Δ
i
U
n
=
U
i
-
1
n
-
2
U
i
n
+
U
i
+
1
n
h
2
,
U
n
=
[
U
1
n
,
U
2
n
,
…
,
U
M
-
1
n
]
T
.
For this linearized system, we have
(36)
u
(
t
)
=
e
(
t
-
t
n
)
A
n
u
(
t
n
)
+
∫
t
n
t
‍
e
(
t
-
s
)
A
n
g
n
d
s
,
and therefore
(37)
u
(
t
)
=
e
(
t
-
t
n
)
A
n
u
(
t
n
)
+
[
e
(
t
-
t
n
)
A
n
-
I
]
A
n
-
1
g
n
.


### 3.2. Numerical Method Construction

In this subsection, we construct a splitting time-stepping method based on (37). In this paper, for matrices C
i, we define
(38)
∏
i
=
1
M
-
1
‍
C
i
∶
=
C
M
-
1
C
M
-
2
⋯
C
2
C
1
.
Then, the following is not true in general:
(39)
∏
i
=
1
M
-
1
‍
C
i
=
C
1
C
2
⋯
C
M
-
2
C
M
-
1
.
Let us consider the following splitting:
(40)
A
n
=
∑
i
=
1
M
-
1
‍
A
n
,
i
,
g
n
=
∑
i
=
1
M
-
1
‍
g
n
,
i
,
with
(41)
A
n
,
1
=
1
h
2
[
-
2
β
1
n
β
1
n
⋯
0
0
0
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
0
]
,
A
n
,
M
-
1
=
1
h
2
[
0
0
⋯
0
0
0
0
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
0
0
0
0
⋯
β
M
-
1
n
-
2
β
M
-
1
n
]
,
A
n
,
i
=
1
h
2
[
0
⋯
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯
⋯
⋯
0
⋮
⋮
⋱
⋮
⋮
⋮
⋮
0
⋯
β
i
n
-
2
β
i
n
β
i
n
⋯
0
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
⋯
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯
⋯
⋯
0
]
,
m
m
m
m
m
m
m
m
m
i
=
2,3
,
…
,
M
-
2
,
g
n
,
i
=
0
,
i
=
1,2
,
…
,
M
-
2
,
g
n
,
M
-
1
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
n
α
(
b
)
]
T
,
where 0
∈
ℝ
M
-
1 is the zero vector. Obviously, other choices for g
n
,
i are possible; for example, g
n
,
1
=
g
n, g
n
,
i
=
0
,
i
=
2,3
,
…
,
M
-
1.With the splitting (40), we solve the M
-
1 subproblems:
(42)
d
u
i
(
t
)
d
t
=
A
n
,
i
u
i
(
t
)
+
g
n
,
i
,
u
i
(
t
n
)
=
u
i
-
1
(
t
n
)
,
mmmm
i
=
1,2
,
…
,
M
-
1
,
starting from u
0
(
t
n
)
=
u
(
t
n
), and we take u
(
t
n
+
1
)
=
u
M
-
1
(
t
n
+
1
) to complete the splitting integration step.Then, if the Backward Euler method is used to solve every subproblem, we get(43)
u
n
i
=
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
i
-
1
+
[
I
-
Δ
t
A
n
,
i
]
-
1
Δ
t
g
n
,
i
,
where u
n
i is an approximation of u
i
(
t
n
). It is easily verified that the matrix I
-
Δ
t
A
n
,
i is M-matrix and therefore is nonsingular. By induction, we have
(44)
u
n
+
1
=
u
n
+
1
0
=
u
n
M
-
1
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
+
∑
j
=
1
M
-
1
‍
∏
i
=
j
+
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
[
I
-
Δ
t
A
n
,
j
]
-
1
Δ
t
g
n
,
j
,
which implies that for the splitting (41)
(45)
u
n
+
1
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
+
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
Δ
t
g
n
,
M
-
1
.
For other choices of g
n
,
i, we can similarly obtain corresponding splitting schemes.By brief calculation, one can obtain an explicit expression of matrix[
I
-
Δ
t
A
n
,
i
]
-
1:
(46)
[
I
-
Δ
t
A
n
,
1
]
-
1
=
[
h
2
h
2
+
2
Δ
t
β
1
n
Δ
t
β
1
n
h
2
+
2
Δ
t
β
1
n
⋯
0
0
1
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
1
]
,
[
I
-
Δ
t
A
n
,
i
]
-
1



=

[


1
0
⋯
0
0
0
⋯
0
0
0
1
⋯
0
0
0
⋯
0
0
⋮
⋮
⋱
⋮
⋮
⋮
⋮
⋮
⋮
0
0
⋯
1
0
0
⋯
0
0
0
0
⋯
Δ
t
β
i
n
h
2


+


2
Δ
t
β
i
n
h
2
h
2


+


2
Δ
t
β
i
n
Δ
t
β
i
n
h
2


+


2
Δ
t
β
i
n
⋯
0
0
0
0
⋯
0
0
1
⋯
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
0
⋯
0
0
0
⋯
1
0
0
0
⋯
0
0
0
⋯
0
1


]
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
2,3
,
⋯
,
M
-
2
,
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
=
[
1
0
⋯
0
0
0
1
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
1
0
0
0
⋯
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
h
2
h
2
+
2
Δ
t
β
M
-
1
n
]
.
Then, we get the following explicit formula for computing u
n
+
1
=
[
U
1
n
+
1
,
U
2
n
+
1
,
…
,
U
M
-
1
n
+
1
]
T, where U
i
n is an approximation of U
(
S
i
,
t
n
):
(47)
U
1
n
+
1
=
h
2
h
2
+
2
Δ
t
β
1
n
U
1
n
+
Δ
t
β
1
n
h
2
+
2
Δ
t
β
1
n
U
2
n
,
U
i
n
+
1

=

Δ
t
β
i
n
h
2
+
2
Δ
t
β
i
n
U
i
-
1
n
+
1
+
h
2
h
2
+
2
Δ
t
β
i
n
U
i
n
+
Δ
t
β
i
n
h
2
+
2
Δ
t
β
i
n
U
i
+
1
n
,
m
m
m
m
m
m
f
m
m
i
=
2,3
,
…
,
M
-
2
,
U
M
-
1
n
+
1
=
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
U
M
-
2
n
+
1
+
h
2
h
2
+
2
Δ
t
β
M
-
1
n
U
M
-
1
n
+
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
α
(
b
)
.For the finite-dimensional discrete system (33), the splitting in the numerical scheme (45) is such that all computations become effectively “truly” one-dimensional. The numerical scheme (45) can be viewed as a “limit” version of LOD-Backward Euler method (see Chapter IV on splitting methods in [17]).


## 3.1. Linearization System

Let us sett
n
=
n
Δ
t, n
=
1,2
,
…
,
N, for the temporal stepsize Δ
t
=
T
/
N. To linearize the nonlinear system (23), we allow the nonlinearities in (23) to lag one step behind and obtain the following linear system:
(33)
u
′
(
t
)
=
A
(
u
(
t
n
)
)
u
(
t
)
+
g
(
u
(
t
n
)
)
=
A
n
u
(
t
)
+
g
n
,
m
t
∈
[
t
n
,
t
n
+
1
]
,
with
(34)
A
n
=
1
h
2
tridiag
(
β
i
n
,
-
2
β
i
n
,
β
i
n
)
,
g
n
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
n
α
(
b
)
]
T
,
where
(35)
β
i
n
=
σ
i
,
n
2
S
i
2
,
σ
i
,
n
2
=
1
2
σ
0
2
(
1
+
Ψ
i
n
)
,
Ψ
i
n
=
Ψ
(
a
2
S
i
2
Δ
i
U
n
)
,
Δ
i
U
n
=
U
i
-
1
n
-
2
U
i
n
+
U
i
+
1
n
h
2
,
U
n
=
[
U
1
n
,
U
2
n
,
…
,
U
M
-
1
n
]
T
.
For this linearized system, we have
(36)
u
(
t
)
=
e
(
t
-
t
n
)
A
n
u
(
t
n
)
+
∫
t
n
t
‍
e
(
t
-
s
)
A
n
g
n
d
s
,
and therefore
(37)
u
(
t
)
=
e
(
t
-
t
n
)
A
n
u
(
t
n
)
+
[
e
(
t
-
t
n
)
A
n
-
I
]
A
n
-
1
g
n
.


## 3.2. Numerical Method Construction

In this subsection, we construct a splitting time-stepping method based on (37). In this paper, for matrices C
i, we define
(38)
∏
i
=
1
M
-
1
‍
C
i
∶
=
C
M
-
1
C
M
-
2
⋯
C
2
C
1
.
Then, the following is not true in general:
(39)
∏
i
=
1
M
-
1
‍
C
i
=
C
1
C
2
⋯
C
M
-
2
C
M
-
1
.
Let us consider the following splitting:
(40)
A
n
=
∑
i
=
1
M
-
1
‍
A
n
,
i
,
g
n
=
∑
i
=
1
M
-
1
‍
g
n
,
i
,
with
(41)
A
n
,
1
=
1
h
2
[
-
2
β
1
n
β
1
n
⋯
0
0
0
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
0
]
,
A
n
,
M
-
1
=
1
h
2
[
0
0
⋯
0
0
0
0
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
0
0
0
0
⋯
β
M
-
1
n
-
2
β
M
-
1
n
]
,
A
n
,
i
=
1
h
2
[
0
⋯
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯
⋯
⋯
0
⋮
⋮
⋱
⋮
⋮
⋮
⋮
0
⋯
β
i
n
-
2
β
i
n
β
i
n
⋯
0
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
⋯
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯
⋯
⋯
0
]
,
m
m
m
m
m
m
m
m
m
i
=
2,3
,
…
,
M
-
2
,
g
n
,
i
=
0
,
i
=
1,2
,
…
,
M
-
2
,
g
n
,
M
-
1
=
1
h
2
[
0,0
,
…
,
0
,
β
M
-
1
n
α
(
b
)
]
T
,
where 0
∈
ℝ
M
-
1 is the zero vector. Obviously, other choices for g
n
,
i are possible; for example, g
n
,
1
=
g
n, g
n
,
i
=
0
,
i
=
2,3
,
…
,
M
-
1.With the splitting (40), we solve the M
-
1 subproblems:
(42)
d
u
i
(
t
)
d
t
=
A
n
,
i
u
i
(
t
)
+
g
n
,
i
,
u
i
(
t
n
)
=
u
i
-
1
(
t
n
)
,
mmmm
i
=
1,2
,
…
,
M
-
1
,
starting from u
0
(
t
n
)
=
u
(
t
n
), and we take u
(
t
n
+
1
)
=
u
M
-
1
(
t
n
+
1
) to complete the splitting integration step.Then, if the Backward Euler method is used to solve every subproblem, we get(43)
u
n
i
=
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
i
-
1
+
[
I
-
Δ
t
A
n
,
i
]
-
1
Δ
t
g
n
,
i
,
where u
n
i is an approximation of u
i
(
t
n
). It is easily verified that the matrix I
-
Δ
t
A
n
,
i is M-matrix and therefore is nonsingular. By induction, we have
(44)
u
n
+
1
=
u
n
+
1
0
=
u
n
M
-
1
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
+
∑
j
=
1
M
-
1
‍
∏
i
=
j
+
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
[
I
-
Δ
t
A
n
,
j
]
-
1
Δ
t
g
n
,
j
,
which implies that for the splitting (41)
(45)
u
n
+
1
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
n
+
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
Δ
t
g
n
,
M
-
1
.
For other choices of g
n
,
i, we can similarly obtain corresponding splitting schemes.By brief calculation, one can obtain an explicit expression of matrix[
I
-
Δ
t
A
n
,
i
]
-
1:
(46)
[
I
-
Δ
t
A
n
,
1
]
-
1
=
[
h
2
h
2
+
2
Δ
t
β
1
n
Δ
t
β
1
n
h
2
+
2
Δ
t
β
1
n
⋯
0
0
1
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
1
]
,
[
I
-
Δ
t
A
n
,
i
]
-
1



=

[


1
0
⋯
0
0
0
⋯
0
0
0
1
⋯
0
0
0
⋯
0
0
⋮
⋮
⋱
⋮
⋮
⋮
⋮
⋮
⋮
0
0
⋯
1
0
0
⋯
0
0
0
0
⋯
Δ
t
β
i
n
h
2


+


2
Δ
t
β
i
n
h
2
h
2


+


2
Δ
t
β
i
n
Δ
t
β
i
n
h
2


+


2
Δ
t
β
i
n
⋯
0
0
0
0
⋯
0
0
1
⋯
0
0
⋮
⋮
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
0
⋯
0
0
0
⋯
1
0
0
0
⋯
0
0
0
⋯
0
1


]
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
2,3
,
⋯
,
M
-
2
,
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
=
[
1
0
⋯
0
0
0
1
⋯
0
0
⋮
⋮
⋱
⋮
⋮
0
0
⋯
1
0
0
0
⋯
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
h
2
h
2
+
2
Δ
t
β
M
-
1
n
]
.
Then, we get the following explicit formula for computing u
n
+
1
=
[
U
1
n
+
1
,
U
2
n
+
1
,
…
,
U
M
-
1
n
+
1
]
T, where U
i
n is an approximation of U
(
S
i
,
t
n
):
(47)
U
1
n
+
1
=
h
2
h
2
+
2
Δ
t
β
1
n
U
1
n
+
Δ
t
β
1
n
h
2
+
2
Δ
t
β
1
n
U
2
n
,
U
i
n
+
1

=

Δ
t
β
i
n
h
2
+
2
Δ
t
β
i
n
U
i
-
1
n
+
1
+
h
2
h
2
+
2
Δ
t
β
i
n
U
i
n
+
Δ
t
β
i
n
h
2
+
2
Δ
t
β
i
n
U
i
+
1
n
,
m
m
m
m
m
m
f
m
m
i
=
2,3
,
…
,
M
-
2
,
U
M
-
1
n
+
1
=
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
U
M
-
2
n
+
1
+
h
2
h
2
+
2
Δ
t
β
M
-
1
n
U
M
-
1
n
+
Δ
t
β
M
-
1
n
h
2
+
2
Δ
t
β
M
-
1
n
α
(
b
)
.For the finite-dimensional discrete system (33), the splitting in the numerical scheme (45) is such that all computations become effectively “truly” one-dimensional. The numerical scheme (45) can be viewed as a “limit” version of LOD-Backward Euler method (see Chapter IV on splitting methods in [17]).


## 4. Properties of the Numerical Scheme

In this section, we investigate some properties of the numerical scheme proposed here.


### 4.1. Stability

In [17], the stability of LOD-Backward Euler method is provided. The following theorem shows the stability of the numerical scheme (45).Theorem 5.
The numerical scheme (45) is unconditionally stable in both the spectral norm ∥
·
∥ and the maximum norm ∥
·
∥
∞; that is, one has the following stability inequalities:
(48)
∥
u
n
-
u
~
n
∥
≤
∥
u
0
-
u
~
0
∥
,
∀
n
≥
0
,
(49)
∥
u
n
-
u
~
n
∥
∞
≤
∥
u
0
-
u
~
0
∥
∞
,
∀
n
≥
0
.
It should be pointed out that (49) can be regarded as a discrete version of maximum norm contractivity (26).Proof.
To obtain the spectral norm stability inequality (48), we apply the Gerschgorin theorem to the matrix A
n
,
i and get that the nonzero matrix eigenvalues line in the disc
(50)
|
z
+
2
σ
i
,
n
2
S
i
2
|
≤
2
σ
i
,
n
2
S
i
2
,
and therefore they are nonpositive. Then, any of the eigenvalues ζ
i of the matrix [
I
-
Δ
t
A
n
,
i
]
-
1, corresponding to the eigenvalues η
i of A
n
,
i, satisfies
(51)
|
ζ
i
|
≤
1
|
1
-
Δ
t
η
i
|
≤
1
.
Thus, it is easy to obtain ∏
i
=
1
M
-
1
‍
|
ζ
i
|
≤
1 and ρ
(
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
)
≤
1, where ρ
(
·
) denotes the spectral radius of the matrix. The stability inequality (48) follows from this estimate.
We note that(52)
μ
∞
[
A
n
,
i
]
≤
0
,
i
=
1,2
,
…
,
M
-
1
,
which directly leads to the maximum norm stability inequality (49) (see [17]).


### 4.2.  Positivity

A nice property of the numerical scheme for the pricing equation is positivity preserving, since the value of option is nonnegative.Theorem 6.
The numerical scheme (45) is unconditionally positivity preserving; that is, the solution of (45) is positive on each time level t
n
+
1, n
=
0,1
,
…
,
N
-
1, if u
0 is positive.Proof.
Since all entries of the matrices[
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1, are nonnegative, a nonnegative solution u
n
+
1 of (45) on each time level can be obtained in view of the nonnegative property of g
n
,
i and u
0. This completes the proof.


### 4.3.  Monotonicity

Let us now consider the monotonicity of the numerical scheme (45). To do this, we note that scheme (45) can be written as
(53)
U
i
n
+
1
=
H
(
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
)
,
i
=
1,2
,
…
,
M
-
1
.
For scheme (53), we have the following definition of monotonicity.Definition 7 (see [19]; see also [20, 21]).
Scheme (53) is said to be monotone if and only if H is nondecreasing in each argument.The following theorem shows the monotonicity of the numerical scheme (45).Theorem 8.
The numerical scheme (45) is unconditionally monotone.Proof.
From the proof of Theorem6, it is known that all entries of the matrices [
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1, are nonnegative. Then, H is nondecreasing in each argument, and therefore the scheme is monotone.


### 4.4.  Maximum Principle

From (43) and (46), it is not difficult to obtain
(54)
min
⁡
{
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
}
≤
U
i
n
+
1
≤
max
⁡
{
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
}
,
m
m
i
=
1,2
,
…
,
M
-
1
.
Then, we have the maximum principle
(55)
min
⁡
j
{
U
j
0
}
≤
U
i
n
≤
max
⁡
j
{
U
j
0
}
,
i
=
1,2
,
…
,
M
-
1
,
∀
n
≥
0
,
which also implies maximum norm estimate; that is, ∥
u
n
∥
∞
≤
∥
u
0
∥
∞.


### 4.5.  Local Error Analysis and Consistency

Now, we consider the error. Letϵ
n
=
u
(
t
n
)
-
u
~
n denote the local discretization error, that is, the error introduced in one single step of the method. To bound the local discretization errors ϵ
n, we need the following truncation error estimation.Lemma 9 (see [22]).
For splittingA
n
=
∑
i
=
1
M
-
1
A
n
,
i, one has
(56)
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
=
Δ
t
2
2
∑
i
=
1
M
-
1
‍
∑
j
=
i
+
1
M
-
1
‍
[
A
n
,
i
,
A
n
,
j
]
u
(
t
n
)
+
O
(
Δ
t
3
)
,
with
(57)
[
A
n
,
i
,
A
n
,
j
]
=
A
n
,
i
A
n
,
j
-
A
n
,
j
A
n
,
i
being the commutator of A
n
,
i and A
n
,
j, and
(58)
[
e
Δ
t
A
n
-
1
2
(
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
+
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
M
-
i
)
]
u
(
t
n
)
=
O
(
Δ
t
3
)
.Theorem 10.
With the previous notation, the local discretization errorϵ
n of the numerical scheme (45) is given by
(59)
ϵ
n
=
O
(
Δ
t
2
)
.Proof.
It follows from (37) that
(60)
u
(
t
n
+
1
)
=
e
Δ
t
A
n
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
u
(
t
n
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
.
Because of the expansion formula
(61)
e
Δ
t
A
n
=
∑
k
=
0
∞
‍
1
k
!
(
Δ
t
A
n
)
k
=
I
+
Δ
t
A
n
+
1
2
(
Δ
t
)
2
A
n
2
+
⋯
,
we have
(62)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
Δ
t
∑
k
=
0
∞
‍
1
(
k
+
2
)
!
(
Δ
t
A
n
)
k
g
n
=
Δ
t
∑
i
=
1
M
-
1
‍
(
I
+
1
2
Δ
t
A
n
+
1
6
(
Δ
t
)
2
A
n
2
+
⋯
)
g
n
,
i
.
We also have the following expansion formula similar to (61):
(63)
e
Δ
t
A
n
,
i
=
∑
k
=
0
∞
‍
1
k
!
(
Δ
t
A
n
,
i
)
k
=
I
+
Δ
t
A
n
,
i
+
1
2
(
Δ
t
)
2
A
n
,
i
2
+
1
6
(
Δ
t
)
3
A
n
,
i
3
+
⋯
,
m
m
m
m
m
m
m
m
m
m
m
m
m
m
i
=
1,2
,
…
,
M
-
1
.
It is easy to obtain the expansion formula of [
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1:
(64)
[
I
-
Δ
t
A
n
,
i
]
-
1
=
I
+
∑
k
=
1
∞
‍
(
Δ
t
A
n
,
i
)
k
=
I
+
Δ
t
A
n
,
i
+
(
Δ
t
)
2
A
n
,
i
2
+
(
Δ
t
)
3
A
n
,
i
3
+
⋯
.
Comparing (63) with (64) yields, for i
=
1,2
,
…
,
M
-
1,
(65)
e
Δ
t
A
n
,
i
=
[
I
-
Δ
t
A
n
,
i
]
-
1
-
1
2
(
Δ
t
)
2
A
n
,
i
2
+
⋯
,
which implies that
(66)
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
-
1
2
∑
i
=
1
M
-
1
‍
(
Δ
t
)
2
A
n
,
i
2
+
⋯
.
From (63) and (62), one gets
(67)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
Δ
t
∑
i
=
1
M
-
1
(
∏
j
=
i
M
-
1
‍
[
I
-
Δ
t
A
n
,
j
]
-
1
+
∑
j
=
1
i
-
1
‍
1
2
Δ
t
A
n
,
j
‍
mmmmmmmm
-
∑
j
=
i
M
-
1
‍
1
2
Δ
t
A
n
,
j
+
O
(
Δ
t
2
A
n
2
)
)
g
n
,
i
.
For splitting (41), since g
n
,
i
=
0, i
=
1,2
,
…
,
M
-
2, we have
(68)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
(
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
+
∑
j
=
1
M
-
2
‍
1
2
Δ
t
A
n
,
j
m
m
m
f
m
m
-
1
2
Δ
t
A
n
,
M
-
1
+
O
(
Δ
t
2
A
n
2
)
∑
j
=
1
M
-
2
‍
)
Δ
t
g
n
,
M
-
1
.
Substituting (66) and (68) into (60) obtains
(69)
u
(
t
n
+
1
)
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
(
t
n
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
+
O
(
Δ
t
2
)
=
u
~
n
+
1
-
1
2
Δ
t
2
A
n
,
M
-
1
g
n
,
M
-
1
+
∑
j
=
1
M
-
2
‍
1
2
Δ
t
2
A
n
,
j
g
n
,
M
-
1
+
O
(
Δ
t
2
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
.
Then, (59) follows from the above equation and (56).We can now immediately conclude that the splitting scheme (45) is consistent.


### 4.6.  Viscosity Solutions and Convergence

The so-called viscosity solutions are the meaningful solutions in financial applications. This has been shown in [23, 24]. Putting all results together, we can formulate the following convergence result.Theorem 11.
The splitting Backward Euler scheme (45) unconditionally converges to the viscosity solution of (14)-(15).Proof.
The convergence of the fully-discrete scheme (45) is a direct result of the consistency, the stability, and the monotonicity of this scheme [25].


## 4.1. Stability

In [17], the stability of LOD-Backward Euler method is provided. The following theorem shows the stability of the numerical scheme (45).Theorem 5.
The numerical scheme (45) is unconditionally stable in both the spectral norm ∥
·
∥ and the maximum norm ∥
·
∥
∞; that is, one has the following stability inequalities:
(48)
∥
u
n
-
u
~
n
∥
≤
∥
u
0
-
u
~
0
∥
,
∀
n
≥
0
,
(49)
∥
u
n
-
u
~
n
∥
∞
≤
∥
u
0
-
u
~
0
∥
∞
,
∀
n
≥
0
.
It should be pointed out that (49) can be regarded as a discrete version of maximum norm contractivity (26).Proof.
To obtain the spectral norm stability inequality (48), we apply the Gerschgorin theorem to the matrix A
n
,
i and get that the nonzero matrix eigenvalues line in the disc
(50)
|
z
+
2
σ
i
,
n
2
S
i
2
|
≤
2
σ
i
,
n
2
S
i
2
,
and therefore they are nonpositive. Then, any of the eigenvalues ζ
i of the matrix [
I
-
Δ
t
A
n
,
i
]
-
1, corresponding to the eigenvalues η
i of A
n
,
i, satisfies
(51)
|
ζ
i
|
≤
1
|
1
-
Δ
t
η
i
|
≤
1
.
Thus, it is easy to obtain ∏
i
=
1
M
-
1
‍
|
ζ
i
|
≤
1 and ρ
(
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
)
≤
1, where ρ
(
·
) denotes the spectral radius of the matrix. The stability inequality (48) follows from this estimate.
We note that(52)
μ
∞
[
A
n
,
i
]
≤
0
,
i
=
1,2
,
…
,
M
-
1
,
which directly leads to the maximum norm stability inequality (49) (see [17]).


## 4.2.  Positivity

A nice property of the numerical scheme for the pricing equation is positivity preserving, since the value of option is nonnegative.Theorem 6.
The numerical scheme (45) is unconditionally positivity preserving; that is, the solution of (45) is positive on each time level t
n
+
1, n
=
0,1
,
…
,
N
-
1, if u
0 is positive.Proof.
Since all entries of the matrices[
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1, are nonnegative, a nonnegative solution u
n
+
1 of (45) on each time level can be obtained in view of the nonnegative property of g
n
,
i and u
0. This completes the proof.


## 4.3.  Monotonicity

Let us now consider the monotonicity of the numerical scheme (45). To do this, we note that scheme (45) can be written as
(53)
U
i
n
+
1
=
H
(
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
)
,
i
=
1,2
,
…
,
M
-
1
.
For scheme (53), we have the following definition of monotonicity.Definition 7 (see [19]; see also [20, 21]).
Scheme (53) is said to be monotone if and only if H is nondecreasing in each argument.The following theorem shows the monotonicity of the numerical scheme (45).Theorem 8.
The numerical scheme (45) is unconditionally monotone.Proof.
From the proof of Theorem6, it is known that all entries of the matrices [
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1, are nonnegative. Then, H is nondecreasing in each argument, and therefore the scheme is monotone.


## 4.4.  Maximum Principle

From (43) and (46), it is not difficult to obtain
(54)
min
⁡
{
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
}
≤
U
i
n
+
1
≤
max
⁡
{
U
i
-
1
n
,
U
i
n
,
U
i
+
1
n
}
,
m
m
i
=
1,2
,
…
,
M
-
1
.
Then, we have the maximum principle
(55)
min
⁡
j
{
U
j
0
}
≤
U
i
n
≤
max
⁡
j
{
U
j
0
}
,
i
=
1,2
,
…
,
M
-
1
,
∀
n
≥
0
,
which also implies maximum norm estimate; that is, ∥
u
n
∥
∞
≤
∥
u
0
∥
∞.


## 4.5.  Local Error Analysis and Consistency

Now, we consider the error. Letϵ
n
=
u
(
t
n
)
-
u
~
n denote the local discretization error, that is, the error introduced in one single step of the method. To bound the local discretization errors ϵ
n, we need the following truncation error estimation.Lemma 9 (see [22]).
For splittingA
n
=
∑
i
=
1
M
-
1
A
n
,
i, one has
(56)
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
=
Δ
t
2
2
∑
i
=
1
M
-
1
‍
∑
j
=
i
+
1
M
-
1
‍
[
A
n
,
i
,
A
n
,
j
]
u
(
t
n
)
+
O
(
Δ
t
3
)
,
with
(57)
[
A
n
,
i
,
A
n
,
j
]
=
A
n
,
i
A
n
,
j
-
A
n
,
j
A
n
,
i
being the commutator of A
n
,
i and A
n
,
j, and
(58)
[
e
Δ
t
A
n
-
1
2
(
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
+
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
M
-
i
)
]
u
(
t
n
)
=
O
(
Δ
t
3
)
.Theorem 10.
With the previous notation, the local discretization errorϵ
n of the numerical scheme (45) is given by
(59)
ϵ
n
=
O
(
Δ
t
2
)
.Proof.
It follows from (37) that
(60)
u
(
t
n
+
1
)
=
e
Δ
t
A
n
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
u
(
t
n
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
.
Because of the expansion formula
(61)
e
Δ
t
A
n
=
∑
k
=
0
∞
‍
1
k
!
(
Δ
t
A
n
)
k
=
I
+
Δ
t
A
n
+
1
2
(
Δ
t
)
2
A
n
2
+
⋯
,
we have
(62)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
Δ
t
∑
k
=
0
∞
‍
1
(
k
+
2
)
!
(
Δ
t
A
n
)
k
g
n
=
Δ
t
∑
i
=
1
M
-
1
‍
(
I
+
1
2
Δ
t
A
n
+
1
6
(
Δ
t
)
2
A
n
2
+
⋯
)
g
n
,
i
.
We also have the following expansion formula similar to (61):
(63)
e
Δ
t
A
n
,
i
=
∑
k
=
0
∞
‍
1
k
!
(
Δ
t
A
n
,
i
)
k
=
I
+
Δ
t
A
n
,
i
+
1
2
(
Δ
t
)
2
A
n
,
i
2
+
1
6
(
Δ
t
)
3
A
n
,
i
3
+
⋯
,
m
m
m
m
m
m
m
m
m
m
m
m
m
m
i
=
1,2
,
…
,
M
-
1
.
It is easy to obtain the expansion formula of [
I
-
Δ
t
A
n
,
i
]
-
1, i
=
1,2
,
…
,
M
-
1:
(64)
[
I
-
Δ
t
A
n
,
i
]
-
1
=
I
+
∑
k
=
1
∞
‍
(
Δ
t
A
n
,
i
)
k
=
I
+
Δ
t
A
n
,
i
+
(
Δ
t
)
2
A
n
,
i
2
+
(
Δ
t
)
3
A
n
,
i
3
+
⋯
.
Comparing (63) with (64) yields, for i
=
1,2
,
…
,
M
-
1,
(65)
e
Δ
t
A
n
,
i
=
[
I
-
Δ
t
A
n
,
i
]
-
1
-
1
2
(
Δ
t
)
2
A
n
,
i
2
+
⋯
,
which implies that
(66)
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
-
1
2
∑
i
=
1
M
-
1
‍
(
Δ
t
)
2
A
n
,
i
2
+
⋯
.
From (63) and (62), one gets
(67)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
Δ
t
∑
i
=
1
M
-
1
(
∏
j
=
i
M
-
1
‍
[
I
-
Δ
t
A
n
,
j
]
-
1
+
∑
j
=
1
i
-
1
‍
1
2
Δ
t
A
n
,
j
‍
mmmmmmmm
-
∑
j
=
i
M
-
1
‍
1
2
Δ
t
A
n
,
j
+
O
(
Δ
t
2
A
n
2
)
)
g
n
,
i
.
For splitting (41), since g
n
,
i
=
0, i
=
1,2
,
…
,
M
-
2, we have
(68)
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
=
(
[
I
-
Δ
t
A
n
,
M
-
1
]
-
1
+
∑
j
=
1
M
-
2
‍
1
2
Δ
t
A
n
,
j
m
m
m
f
m
m
-
1
2
Δ
t
A
n
,
M
-
1
+
O
(
Δ
t
2
A
n
2
)
∑
j
=
1
M
-
2
‍
)
Δ
t
g
n
,
M
-
1
.
Substituting (66) and (68) into (60) obtains
(69)
u
(
t
n
+
1
)
=
∏
i
=
1
M
-
1
‍
[
I
-
Δ
t
A
n
,
i
]
-
1
u
(
t
n
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
+
[
e
Δ
t
A
n
-
I
]
A
n
-
1
g
n
+
O
(
Δ
t
2
)
=
u
~
n
+
1
-
1
2
Δ
t
2
A
n
,
M
-
1
g
n
,
M
-
1
+
∑
j
=
1
M
-
2
‍
1
2
Δ
t
2
A
n
,
j
g
n
,
M
-
1
+
O
(
Δ
t
2
)
+
[
e
Δ
t
A
n
-
∏
i
=
1
M
-
1
‍
e
Δ
t
A
n
,
i
]
u
(
t
n
)
.
Then, (59) follows from the above equation and (56).We can now immediately conclude that the splitting scheme (45) is consistent.


## 4.6.  Viscosity Solutions and Convergence

The so-called viscosity solutions are the meaningful solutions in financial applications. This has been shown in [23, 24]. Putting all results together, we can formulate the following convergence result.Theorem 11.
The splitting Backward Euler scheme (45) unconditionally converges to the viscosity solution of (14)-(15).Proof.
The convergence of the fully-discrete scheme (45) is a direct result of the consistency, the stability, and the monotonicity of this scheme [25].


## 5. Numerical Experiments

In order to illustrate the stability and convergence properties of our proposed scheme, in this section, we present several numerical experiments in which the vanilla call option and the European butterfly spread are considered. To obtain the numerical solution of nonlinear Black-Scholes equations (14)–(16), we should first solve the ODE (3). In this paper, we use the symmetric midpoint scheme to solve numerically the ODE (3). The utility function Ψ
(
x
) can be done by linear interpolation.Example 1.
Let us consider a vanilla European call option with(70)
f
(
s
)
=
max
⁡
(
s
-
E
,
0
)
.
This example comes from [14]. The parameters are the following: the strike price E
=
100, the volatility σ
0
=
0.2, the interest rate r
=
0.02, the maturity date T
=
1 year, and the artificial boundary location b
=
200. We first consider a small time stepsize Δ
t
=
1
/
N
=
0.0002 and a spatial stepsize h
=
200
/
M
=
4 as done in [14]. The numerical results V
(
s
,
τ
) computed by the numerical scheme (45) are plotted in Figure 1 in (
s
,
τ
) coordinate. The case of a
=
0 is linear model, and the cases of a
=
0.015 and a
=
0.1 are nonlinear model. The numerical results are the same as those showed in [14].Figure 1
Option pricing of a vanilla European call option for several values of parametera at maturity time T
=
1 year, where h
=
4 and Δ
t
=
0.0002.To further confirm that this scheme is unconditionally stable, monotone, and positivity preserving, we also calculate the numerical solutions, which are presented in Figures2 and 3, by scheme (45) with larger time stepsizes Δ
t
=
0.002 and Δ
t
=
0.02, respectively. From Figures 2 and 3, we observe that there are no stability issues for scheme (45). We also note that the scheme proposed in [14] produces the wrong numerical solution when the stepsize Δ
t
=
0.00027027 such that the stability condition (32) is not satisfied (see [14]). This reveals that the numerical scheme (45) has better stability properties than the scheme proposed in [14].Figure 2
Option pricing of a vanilla European call option for several values of parametera at maturity time T
=
1 year, where h
=
4 and Δ
t
=
0.002.Figure 3
Option pricing of a vanilla European call option for several values of parametera at maturity time T
=
1 year, where h
=
4 and Δ
t
=
0.02.To illustrate the convergence of scheme (45), we show the option value at S
=
E
=
100 (or at s
=
98.01987 at time to maturity being 1 year), the differences between successive approximation processes, and the ratios between the differences in Table 1. Since scheme (45) is of order two in space and order one in time, the number of spatial mesh points is double, but the number of temporal grid points is quadruple. Table 1 only shows the numerical results for the case a
=
0.015, which confirm our theoretical analysis results. The numerical results for the cases of a
=
0 and a
=
0.1 are similar and we omit these here.Table 1
Convergence results for the scheme (45) for the transaction costs model together with a vanilla call option, where a
=
0.015.
M
N
Value
Difference
Ratio
5
0
50
8.391857
—
—
100
200
8.416984
0.025127
—
200
800
8.423089
0.006105
4.115806
400
3200
8.424567
0.001487
4.105581From the data presented in Figures1–3 and Table 1, the conclusion can be drawn: scheme (45) is unconditionally stable and convergent.Example 2 (european butterfly spread).
This example comes from [16]. A butterfly spread consists of two long positions in two calls with E
1, E
2 and a short position in two calls at strike S
=
(
E
1
+
E
2
)
/
2. The payoff function f
(
S
) is given by
(71)
f
(
S
)
=
max
⁡
(
S
-
E
1
,
0
)
-
2
max
⁡
(
S
-
E
,
0
)
+
max
⁡
(
S
-
E
2
,
0
)
.
Let E
1
=
0.8, E
2
=
1.2, the volatility σ
0
=
0.5, the interest rate r
=
0.04, the maturity date T
=
0.5, and the artificial boundary location b
=
10. We first compute the numerical solution by scheme (45) with the time stepsize Δ
t
=
0.5
/
N
=
0.0005 and the spatial stepsize h
=
10
/
M
=
0.1 similar to those in [16]. The numerical results are presented in Figure 4. These numerical data together with those in Figures 5 and 6 show that this scheme proposed here is unconditionally stable.Figure 4
Option pricing of a European butterfly spread for several values of parametera at maturity time T
=
0.5, where h
=
0.1 and Δ
t
=
0.0005.Figure 5
Option pricing of a European butterfly spread for several values of parametera at maturity time T
=
0.5, where h
=
0.1 and Δ
t
=
0.005.Figure 6
Option pricing of a European butterfly spread for several values of parametera at maturity time T
=
0.5, where h
=
0.1 and Δ
t
=
0.05.In [16], the authors compared the numerical results obtained by their nonstandard scheme and the schemes proposed in [14, 15] and showed that their scheme produces better numerical solutions than the schemes in [14, 15] with the same stepsizes. To compare our scheme with the nonstandard scheme proposed in [16], we first observe that both of the numerical schemes are unconditionally stable. Then, we will investigate their convergence and compare their convergence order. To this end, we calculate the discrete maximum norms of the errors
(72)
E
∞
Sp
=
max
⁡
i
∥
U
i
N
,
Sp
-
U
~
i
N
∥
,
E
∞
Ns
=
max
⁡
i
∥
U
i
N
,
Ns
-
U
~
i
N
∥
,
where U
i
N
,
Sp and U
i
N
,
Ns denote the numerical solutions computed by our splitting scheme (45) and the nonstandard scheme proposed in [16] at the maturity date T
=
0.5, respectively, and U
~
i
N denotes the numerical reference solution computed by the Backward Euler method with the standard second-order finite differences on a finer grid with N
=
2560 and M
=
320. Numerical results of both schemes (45) and the nonstandard scheme proposed in [16] are listed in Table 2. This table shows the maximum norm of the errors and the ratios between these errors. It is evident from the numerical data obtained by the splitting scheme (45) that this scheme is convergent with the temporal order one and the spatial order two. This is consistent with the theoretical results presented in this paper. For the nonstandard scheme proposed in [16], however, a reduction in the error is observed. This is not surprising since, as Zhou et al. pointed out in [16], the nonstandard scheme is conditionally consistent and the truncation error really depends on the ratio Δ
t
/
h
2.Table 2
Convergence results for the scheme (45) and the nonstandard scheme proposed in [16] for the transaction costs model together with a butterfly spread, where a
=
0.05.
M
N
E
∞
Sp
Ratio
E
∞
Ns
Ratio
2
0
10
8.205076
e
-
002
—
8.209482
e
-
002
—
40
40
1.753266
e
-
002
4.679881
1.776302
e
-
002
4.621670
80
160
4.409681
e
-
003
3.975947
4.857916
e
-
003
3.656510
160
640
1.127720
e
-
003
3.910262
2.229942
e
-
003
2.178494From the theoretical analysis given in this paper and the numerical results shown in this section, we come to the following remark: the proposed scheme is efficient and reliable.


## 6. Concluding Remarks

In this paper, an unconditionally stable splitting scheme has been proposed to solve the nonlinear option pricing model with transaction costs. This method can be viewed as a “limit” version of LOD-Backward Euler method. This “limit” property that all subproblems are one-dimensional allows us to solve the discrete equation explicitly. As a consequence, this method is computationally efficient. The theoretical analysis carried out in this paper shows that this method is unconditionally stable, monotone, and positivity preserving. We also present several numerical experiments in which the vanilla call option and the European butterfly spread are considered. The theoretical analysis presented and the numerical results shown in this paper confirm that the proposed scheme here is efficient and reliable.



---

*Source: 525207-2014-05-07.xml*