# Algorithms of Finite Difference for Pricing American Options under Fractional Diffusion Models

Authors: Jun Xi; Yanqing Chen; Jianwen Cao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

It is well known that linear complementarity problem (LCP) involving partial integro differential equation (PIDE) arises from pricing American options under Lévy Models. In the case of infinite activity process, the integral part of the PIDE has a singularity, which is generally approximated by a small Brownian component plus a compound Poisson process, in the neighborhood of origin. The PIDE can be reformulated as a fractional partial differential equation (FPDE) under fractional diffusion models, including FMLS (finite moment log stable), CGMY (Carr-Madan-Geman-Yor), and KoBol (Koponen-Boyarchenko-Levendorskii). In this paper, we first present a stable iterative algorithm, which is based on the fractional difference approach and penalty method, to avoid the singularity problem and obtain numerical approximations of first-order accuracy. Then, on the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Numerical tests are performed to demonstrate the effectiveness of the algorithm and the extrapolation method. We believe that
this can be used as necessary tools by the engineers in research.


---

## Body

## 1. Introduction

The holder of American options has the right to exercise at any date prior to maturity and not only at the expiry date, while the holder of European options can exercise his right only at the expiry date. Thus, American options give the holder more freedom; hence, the price of American put is always greater or equal to the price of its European counterpart. There are two types of American options: put option and call option. The put option is a contract that the owner has the right but not the obligation to sell, while the holder of a call option has the right but not the obligation to buy.Black and Scholes [1] and Merton [2] proposed the classic option pricing model called Black-Scholes model in 1973. In order to fit observed empirical market data better, many modifications to the Black-Scholes model have been proposed to build new models like jump-diffusion models proposed by Kou in [3] and Merton in [4] and infinite jump activity models such as KoBol [5, 6] and CGMY [7]. More general jump-diffusion models with stochastic volatility are considered in [8].Compared with the jump-diffusion models, infinite activity models like KoBol and CGMY [7, 9] are more flexible to characterize the price process and give a more realistic description of the price process at various time scales. Both time-series and option data indicate that market indices lack a diffusion component and confirm the value of the infinite activity Lévy model [7].Due to the optimal stopping problem, pricing American options lead to a a linear complementarity problem (LCP) or a variational inequality, which involves a partial integrodifferential equation (PIDE) in the case of models with jumps. Merton et al. [10] proposed a finite difference method for this type of formulation in 1977 under Black-Scholes model. More solutions to the LCP, such as the operator splitting method (OS) [11] under the Black-Scholes model, the PSOR method under the Black-Scholes model, the penalty method [12] under the jump-diffusion models, and the Lagrange multiplier method [13] under the Black-Scholes model, have been studied as well. Generally second-order accurate finite differences are used to discretize the PDE of the LCP for option pricing. In [14–16], fourth-order accurate finite difference methods have been considered for American options.As we know, pricing American options are a free boundary problem [6]; hence, the free boundary and the price of an option must be solved simultaneously. If the models satisfy the smooth pasting principle, we can use two types of methods to find the moving location of free boundary: front tracking [17–19] and front fixing [20–22]. However smooth pasting principle does not always hold under models with jumps unless the models satisfy certain conditions.In the case of pricing American options under the infinite activity jump models, how to deal with the integral part of the PIDE is important because it has a singularity in the neighborhood of origin. Almendral [23] approximates the pure-jump process with a compound Poisson process plus a small Brownian component and accelerates the computation by FFT. An integration by parts technique for rewriting the integrodifferential operator in terms of Volterra operators with a weakly singular kernel is proposed. These methods achieve second-order accuracy.Fractional derivatives provide a more delicate instrument to capture the characteristics of processes or materials in comparison with the integer-order derivatives, which proves to be very useful in many fields [24]. Under certain conditions, the Grünwald-Letnikov definition of fractional derivative is equivalent to the Riemann-Liouville definition of fractional derivative; hence, the Grünwald-Letnikov definition, which is convenient for numerical valuation, is often used for discretization of the fractional derivatives. Due to the instability caused by standard Grünwald-Letnikov estimates, a shifted Grünwald-Letnikov definition, which is unconditionally stable, is proposed [25].Cartea and del-Castillo-Negrete [26] presented the fractional partial differential equation (FPDE) under some infinite activity Lévy models (FMLS, KoBol, and CGMY) as a generalization to PIDE and gave a finite difference method using the numerical valuation of fractional derivatives. The numerical solutions of three fractional partial differential equations have been compared by Marom and Momoniat [27]. Their methods are of first-order accuracy and are only for evaluation of European options. Second-order accurate approximation schemes for fractional derivatives have been proposed in [28, 29]. The approach in [28] is based on the classical Crank-Nicholson method combined with spatial extrapolation to obtain temporally and spatially second-order accurate numerical estimates for fractional diffusion equation. By using proper shifting parameter, the generalized shifted Grünwald-Letnikov definition has been proven to be second-order accurate [29].In this paper, our study focuses on the discretization of fractional derivatives, spatial extrapolation, and design of stable iterative algorithms for pricing American options in the scheme of FPDE. We first present a stable iterative algorithm, which is based on the fractional difference approach and penalty method, to avoid the singularity problem and obtain numerical approximations of first-order accuracy. On the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Discretization of Fractional derivatives results in dense matrices which make the computation expensive, so the fast Fourier transform (FFT) method can be employed to decrease the computation cost to gain an almost linear complexityO
(
N
log
⁡
N
) from O
(
N
2
) using direct solution.The outline of this paper is as follows. In Section2, linear complementarity problem and FPDE arising from pricing American options are introduced. In Section 3, we present the finite difference time and space discretization method. Section 4 describes our method used to solve the systems of linear equations and linear complementarity problems for American options. Numerical experiments are given in Section 5, and finally Section 6 contains conclusions. The main contributions of this paper are the discretization method in Section 3, the algorithm proposed in Section 4, and the numerical experiments in Section 5.


## 2. Mathematical Model for American Options

Because the value of an American call on an asset paying no dividends is equal to the value of the European call with same strike and maturity [2], only American put options are considered here.To introduce the mathematical model, let us first recall the Riemann-Liouville definition [26] which is the most widely known definition of the fractional derivatives:
(1)
a
D
x
γ
f
(
x
)
=
1
Γ
(
n
-
γ
)
(
∂
∂
x
)
n
∫
a
x
‍
(
x
-
y
)
n
-
γ
-
1
f
(
y
)
d
y
,
n
-
1
≤
γ
≤
n
,
x
D
a
γ
f
(
x
)
=
1
Γ
(
n
-
γ
)
(
-
∂
∂
x
)
n
∫
x
a
‍
(
y
-
x
)
n
-
γ
-
1
f
(
y
)
d
y
,
n
-
1
≤
γ
≤
n
,
where γ is an arbitrary order, n is a positive integer, and a is the boundary.Under certain conditions, Grünwald-Letnikov definition is equivalent to the Riemann-Liouville definition; hence, we can use Grünwald-Letnikov definition to evaluate the Riemann-Liouville derivatives of the LCP for pricing American options. Grünwald-Letnikov definition [24] is more convenient for numerical valuation than the Riemann-Liouville definition.Using the known definition of the binomial coefficients(2)
(
γ
r
)
=
γ
(
γ
-
1
)
(
γ
-
2
)
⋯
(
γ
-
r
+
1
)
r
!
,
we have the Grünwald-Letnikov definitions of fractional derivatives
(3)
a
D
x
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
x
-
a
h
-
γ
∑
r
=
0
n
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
)
,
x
D
a
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
a
-
x
h
-
γ
∑
r
=
0
n
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
)
,
where γ is an arbitrary order and a is the boundary.It has been proved that the standard Grünwald-Letnikov definition has stability problem if being used to approximate fractional derivate in the space-fractional advection-dispersion equation. Consider the following:(4)
∂
c
(
x
,
t
)
∂
t
=
-
v
(
r
)
∂
c
(
r
,
t
)
∂
r
+
d
(
r
)
∂
γ
c
(
r
,
t
)
∂
r
γ
+
f
(
r
,
t
)
,
where
1
<
γ
⩽
2
.
So the shifted Grünwald-Letnikov definition was proposed to approximate Liouville fractional derivative A
f
(
x
)
=
d
γ
f
(
x
)
/
d
x
γ:
(5)
A
h
f
(
x
)
=
1
Γ
(
-
γ
)
1
h
γ
∑
k
=
0
∞
Γ
(
k
-
γ
)
Γ
(
k
+
1
)
f
(
x
-
(
k
-
p
)
h
)
,
where p is a nonnegative integer. A
h
f
(
x
) is of first-order accuracy.Generally,p is set to 1. We have shifted Grünwald-Letnikov definitions for the left and right fractional derivatives:
(6)
a
D
x
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
x
-
a
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
+
h
)
,
x
D
a
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
a
-
x
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
-
h
)
,
where γ is an arbitrary order and a is the boundary.In the following, a LCP involving FPDE for pricing American put options under fractional diffusion model is given. We denote the value of American put option byV, the exercise price by K, the expiry time by T, the risk free interest rate by r, and the underlying asset price by e
x. The solution of the following LCP with additional constraints gives the price V of the option.Using the fractional partial differential equation (FPDE) given in [26] for pricing European options
(7)
∂
V
(
x
,
t
)
∂
t
+
A
∂
V
(
x
,
t
)
∂
x
+
B
(
x
)
x
D
∞
γ
(
b
(
x
)
V
(
x
,
t
)
)
+
C
(
x
)
-
∞
D
x
γ
(
c
(
x
)
V
(
x
,
t
)
)
+
D
V
(
x
,
t
)
=
0
,
we have the following LCP for pricing American options:
(8)
g
(
x
)
=
max
⁡
(
K
-
e
x
,
0
)
,
x
∈
R
,
V
-
g
(
x
)
≥
0
,
x
∈
R
,
L
V
×
[
V
-
g
(
x
)
]
=
0
,
x
∈
R
,
t
∈
[
0
,
T
]
,
L
V
=
0
,
e
x
>
f
(
t
)
,
t
∈
[
0
,
T
]
,
V
-
g
(
x
)
=
0
,
e
x
≤
f
(
t
)
,
t
∈
[
0
,
T
]
,
L
V
=
0
,
e
x
>
f
(
t
)
,
t
∈
[
0
,
T
]
,
V
⟶
K
,
x
⟶
-
∞
,
t
∈
[
0
,
T
]
,
V
⟶
0
,
x
⟶
∞
,
t
∈
[
0
,
T
]
,
V
(
x
,
T
)
=
g
(
x
)
,
x
∈
R
,
where
(9)
L
V
=
∂
V
(
x
,
t
)
∂
t
+
A
∂
V
(
x
,
t
)
∂
x
+
B
(
x
)
x
D
∞
γ
(
b
(
x
)
V
(
x
,
t
)
)
+
C
(
x
)
-
∞
D
x
γ
(
c
(
x
)
V
(
x
,
t
)
)
+
D
V
(
x
,
t
)
and γ is the fractional order.For KoBol model, we have(10)
L
V
=
∂
V
(
x
,
t
)
∂
t
+
(
r
-
v
-
λ
α
-
1
(
q
-
p
)
)
∂
V
(
x
,
t
)
∂
x
+
1
2
σ
α
[
p
e
λ
x
x
D
∞
α
(
e
-
λ
x
V
(
x
,
t
)
)
000000000
+
q
e
-
λ
x
-
∞
D
x
α
(
e
λ
x
V
(
x
,
t
)
)
]
-
(
r
+
1
2
σ
α
λ
α
)
V
(
x
,
t
)
,
where
(11)
v
=
1
2
σ
α
{
p
(
λ
-
1
)
α
+
q
(
λ
+
1
)
α
-
λ
α
}
,
0
<
α
<
1
,
v
=
1
2
σ
α
{
p
(
λ
-
1
)
α
+
q
(
λ
+
1
)
α
-
λ
α
-
α
λ
α
-
1
(
q
-
p
)
}
,
1
<
α
<
2
.For CGMY model, we have(12)
L
V
=
∂
V
(
x
,
t
)
∂
t
+
(
r
-
d
)
∂
V
(
x
,
t
)
∂
x
+
C
Γ
(
-
Y
)
e
M
x
x
D
∞
Y
(
e
-
M
x
V
(
x
,
t
)
)
+
C
Γ
(
-
Y
)
e
-
G
x
-
∞
D
x
Y
(
e
G
x
V
(
x
,
t
)
)
-
(
r
+
C
Γ
(
-
Y
)
(
M
Y
+
G
Y
)
)
V
(
x
,
t
)
,
where
(13)
d
=
C
Γ
(
Y
)
(
M
-
1
)
Y
-
M
Y
+
(
G
+
1
)
Y
-
G
Y
.The final valueV is given by V
(
x
,
T
)
=
g
(
x
), where g
(
x
) is the payoff function for the option contract. For a put option, it is g
(
x
)
=
max
⁡
(
K
-
e
x
,
0
). We should note that the function f
(
t
) is the free boundary which divides R into two parts. If x
≤
f
(
t
), then we have the equation V
(
x
,
t
)
-
g
(
x
)
=
0, which means that the price V of the option is the exercise price. If x
>
f
(
t
), then equation L
V
=
0 means that the price V is the solution of the above FPDE.


## 3. Space and Time Discretization

The above mathematical model involving the fractional derivatives assumes an infinite domain,x
∈
(
-
∞
,
+
∞
). In order to apply finite difference methods, we will need to discretize and truncate the domain of V
(
x
,
t
). For L
≤
x
≤
R and finite time interval [
0
,
T
], the terminal condition is given by V
(
x
,
T
)
=
g
(
x
). Further, the boundary conditions are V
(
L
,
t
)
=
K and V
(
R
,
t
)
=
0. By previous work in the field [11, 30], we take e
R at four times the exercise price of the option to ensure that the errors generated by the truncation are small enough so as to be negligible. In fact, the choice of L is simpler than that of R because it is easy to choose proper L so as to keep e
L far from the moving free boundary, hence, the errors caused by the left truncation are small enough too.For the space discretization, we use a uniform grid on interval[
L
,
R
] with n
+
1 grid points. The grid step size is denoted by h. Let i
=
0,1
,
2
,
…
,
n. Then we denote space grid point by x
i
=
L
+
i
h. Analogous to the space discretization, a uniform time step Δ
t is used to get m
+
1 gird points on interval [
0
,
T
]. Let j
=
0,1
,
2
,
…
,
m. We denote time grid point by t
j
=
j
Δ
t. In order to perform a finite difference approximation of the partial derivatives, we introduce a simplifying piece of notation along with the discrete approximation
(14)
V
(
x
i
,
t
j
)
=
V
(
L
+
i
h
,
j
Δ
t
)
≈
V
i
,
j
.


### 3.1. Singularity of the Fractional Derivatives

Using the above truncation, we actually assume that-
∞
D
x
γ
f
(
x
)
≈
L
D
x
γ
f
(
x
) and x
D
+
∞
γ
f
(
x
)
≈
x
D
R
γ
f
(
x
). As we know, L
D
x
γ
f
(
x
) is singular at the lower, x
=
L, boundary and x
D
R
γ
f
(
x
) is singular at the upper, x
=
R, boundary. To understand the nature of the singularity, we take the left fractional derivative a
D
x
γ
f
(
x
) as an example.Expanding in Taylor series aroundx
=
a, we have
(15)
a
D
x
γ
f
(
x
)
=
∑
k
=
0
m
f
k
(
a
)
(
x
-
a
)
-
γ
+
k
Γ
(
-
γ
+
k
+
1
)
+
1
Γ
(
-
γ
+
m
+
1
)
∫
a
x
‍
(
x
-
τ
)
m
-
p
f
(
m
+
1
)
(
τ
)
d
τ
.
For γ
>
0, at least, the first term on the right-hand side is singular. Truncation of domain is a must for discretization; so singularity is an important problem that we will clarify in the following discussion on left and right fractional derivatives seperately.The LCP (8) has a function f
(
t
) which is the free boundary at time t. Generally, the free boundary f
(
t
) is far from the lower, x
=
L, boundary. The price of American option V
(
x
,
t
) is equal to g
(
x
), when x
≤
f
(
t
), so we only need to deal with those fractional derivatives which satisfy the condition x
≥
f
(
t
)
>
L; hence, no singularity problems exist for left fractional derivatives.As to the right fractional derivatives, the singularity exists at the upper,x
=
R, boundary if we try to discretize the derivative. Taking into account the boundary condition V
(
R
,
t
)
=
0 and the assumption V
(
x
>
R
,
t
)
=
0, the equality x
D
+
∞
γ
f
(
x
)
=
x
D
R
γ
f
(
x
) holds in such a case.


### 3.2. First-Order Approximation for the Fractional Derivatives

Using the shifted Grünwald-Letnikov definition of fractional derivatives(16)
a
D
x
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
x
-
a
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
+
h
)
,
x
D
a
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
a
-
x
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
-
h
)
,
we have for x
=
n
h
(17)
a
D
x
γ
f
(
x
)
≈
a
Δ
x
γ
f
(
x
)
=
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
+
h
)
,
(18)
x
D
a
γ
f
(
x
)
≈
x
Δ
a
γ
f
(
x
)
=
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
-
h
)
,
which give the first-order accurate approximations.


### 3.3. First-Order Approximation for the Integer-Order Derivatives

There is a first-order space derivative∂
V
(
x
,
t
)
/
∂
x in the LCP (8) to be approximated and we have many choices to use to obtain approximation of order 1 or the higher.For first-order approximation, we have the forward difference(19)
∂
V
i
,
j
∂
x
=
V
i
+
1
,
j
-
V
i
,
j
h
+
O
(
h
)
or the backward difference
(20)
∂
V
i
,
j
∂
x
=
V
i
,
j
-
V
i
-
1
,
j
h
+
O
(
h
)
.In this paper, we choose the backward difference method.


### 3.4. Crank-Nicolson Discretization

AssumingR is the matrix arising from the finite difference scheme, a general time stepping scheme for pricing American options is
(21)
V
i
,
j
-
V
i
,
j
+
1
Δ
t
+
Θ
R
V
i
,
j
+
(
1
-
Θ
)
R
V
i
,
j
+
1
=
0
,
x
i
>
f
(
t
j
)
,
V
i
,
j
=
K
-
e
x
i
,
x
i
≤
f
(
t
j
)
,
where j
=
0,1
,
…
,
m
-
1 and f
(
t
j
) is the free boundary at time step t
j. For above Θ
=
1 we recover the explicit forward Euler scheme. Θ
=
0 gives the implicit backward Euler method. We have the Crank-Nicolson method when Θ
=
1
/
2.


### 3.5. Matrix Formulation

Using (18), (20), and (21), we can now have the discretization of the full FPDE (7):
(22)
V
i
,
j
-
V
i
,
j
+
1
Δ
t
+
Θ
[
A
V
i
,
j
-
V
i
-
1
,
j
h
+
B
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
b
i
+
r
-
1
V
i
+
r
-
1
,
j
000000000
+
C
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
c
i
-
r
+
1
V
i
-
r
+
1
,
j
+
D
V
i
,
j
]
+
(
1
-
Θ
)
[
A
V
i
,
j
+
1
-
V
i
-
1
,
j
+
1
h
∑
r
=
0
i
+
1
‍
00000.00000000
+
B
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
b
i
+
r
-
1
V
i
+
r
-
1
,
j
+
1
000000.0000000
+
C
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
c
i
-
r
+
1
V
i
-
r
+
1
,
j
+
1
+
D
V
i
,
j
+
1
]
,
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
.
v
v
v
i
=
1,2
,
3
,
…
,
n
-
1
,
where
(23)
w
r
=
1
2
π
∫
0
2
π
‍
(
1
-
e
-
i
ψ
)
γ
e
i
r
ψ
d
ψ
,
which can be efficiently computed by FFT.Define matricesG, E, and F such that
(24)
[
G
V
j
]
i
=
(
A
*
Δ
t
h
+
D
*
Δ
t
)
V
i
,
j
-
A
*
Δ
t
h
V
i
-
1
,
j
,
[
E
V
j
]
i
=
Δ
t
B
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
b
i
-
r
+
1
V
i
-
r
+
1
,
j
+
1
,
[
F
V
j
]
i
=
Δ
t
C
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
c
i
+
r
-
1
V
i
+
r
-
1
,
j
+
1
.
Equation (22) then can be rewritten in matrix form as follows:
(25)
[
I
-
(
1
-
Θ
)
G
]
V
j
+
1
=
[
I
+
Θ
G
]
V
j
+
(
1
-
Θ
)
[
E
+
F
]
V
j
+
1
+
Θ
[
E
+
F
]
V
j
.IfΘ is set to 1
/
2, this fractional Crank-Nicholson discretization approach will produce a method with a local truncation error, that is, O
(
Δ
t
2
)
+
O
(
h
).


### 3.6. Spatially Second-Order Approximation Obtained by Extrapolation

In numerical analysis, extrapolation is used to improve the rate of convergence of a sequence. Based on the spatially first-order accurate numerical approximation in 3.5, extrapolation can be employed to improve the low order of spatial convergence. First, we state below the result obtained in [28].Proposition 1.
Let1
≤
α
≤
2, f
∈
C
n
+
3
(
R
) such that all derivatives of f up to order n
+
3 belong to L
1
(
R
). For any integer p
≥
0, define the shifted Grünwald operator by
(26)
Δ
h
,
p
α
f
(
x
)
=
∑
j
=
0
∞
‍
(
-
1
)
j
(
α
j
)
f
(
x
-
(
j
-
p
)
h
)
.
Then, one has for some constants a
l independent of h, f, and x the fact that
(27)
h
-
α
Δ
h
,
p
α
f
(
x
)
=
d
α
d
x
α
f
(
x
)
+
∑
l
=
1
n
-
1
‍
(
a
l
d
α
+
l
d
x
α
+
l
f
(
x
)
)
h
l
+
O
(
h
n
)
uniformly in x
∈
R.
Writing (27) as
(28)
h
-
α
Δ
h
,
p
α
f
(
x
)
=
d
α
d
x
α
f
(
x
)
+
c
1
h
+
c
2
h
2
+
O
(
h
3
)
,
where c
1 and c
2 are independent of h, and (28) shows that the truncation error in the shifted Grünwald finite difference approximation is c
1
h
+
c
2
h
2
+
O
(
h
3
). Now, one can use the Richardson extrapolation method in the following way. First, the Crank-Nicholson method is applied for grid size h and again for grid size h
/
2. This gives two solutions of first-order accuracy for the spatial dimension. We denote by V
h the solution of grid size h and by V
h
/
2 the solution of gird size h
/
2. Then, one has
(29)
V
h
=
V
+
b
1
h
+
b
2
h
2
+
O
(
h
3
)
,
V
h
/
2
=
V
+
b
1
h
2
+
b
2
h
2
4
+
O
(
h
3
)
,
where b
1 and b
2 are real numbers. Taking V
=
2
*
V
h
/
2
-
V
h, it is obvious that V is spatially 2nd-order accurate.
Based on the classical Crank-Nicholson method combined with spatial extrapolation, we obtain a solution with local truncation errorO
(
Δ
t
2
)
+
O
(
h
2
).The extrapolation method can be used at each time step before the maturity of options; because we can obtain two solutions of first-order estimates for diffenent spatial grid size at each time step, then second-order approximation can be computed at each time step.


## 3.1. Singularity of the Fractional Derivatives

Using the above truncation, we actually assume that-
∞
D
x
γ
f
(
x
)
≈
L
D
x
γ
f
(
x
) and x
D
+
∞
γ
f
(
x
)
≈
x
D
R
γ
f
(
x
). As we know, L
D
x
γ
f
(
x
) is singular at the lower, x
=
L, boundary and x
D
R
γ
f
(
x
) is singular at the upper, x
=
R, boundary. To understand the nature of the singularity, we take the left fractional derivative a
D
x
γ
f
(
x
) as an example.Expanding in Taylor series aroundx
=
a, we have
(15)
a
D
x
γ
f
(
x
)
=
∑
k
=
0
m
f
k
(
a
)
(
x
-
a
)
-
γ
+
k
Γ
(
-
γ
+
k
+
1
)
+
1
Γ
(
-
γ
+
m
+
1
)
∫
a
x
‍
(
x
-
τ
)
m
-
p
f
(
m
+
1
)
(
τ
)
d
τ
.
For γ
>
0, at least, the first term on the right-hand side is singular. Truncation of domain is a must for discretization; so singularity is an important problem that we will clarify in the following discussion on left and right fractional derivatives seperately.The LCP (8) has a function f
(
t
) which is the free boundary at time t. Generally, the free boundary f
(
t
) is far from the lower, x
=
L, boundary. The price of American option V
(
x
,
t
) is equal to g
(
x
), when x
≤
f
(
t
), so we only need to deal with those fractional derivatives which satisfy the condition x
≥
f
(
t
)
>
L; hence, no singularity problems exist for left fractional derivatives.As to the right fractional derivatives, the singularity exists at the upper,x
=
R, boundary if we try to discretize the derivative. Taking into account the boundary condition V
(
R
,
t
)
=
0 and the assumption V
(
x
>
R
,
t
)
=
0, the equality x
D
+
∞
γ
f
(
x
)
=
x
D
R
γ
f
(
x
) holds in such a case.


## 3.2. First-Order Approximation for the Fractional Derivatives

Using the shifted Grünwald-Letnikov definition of fractional derivatives(16)
a
D
x
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
x
-
a
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
+
h
)
,
x
D
a
γ
f
(
x
)
=
lim
⁡
h
→
0
n
h
=
a
-
x
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
-
h
)
,
we have for x
=
n
h
(17)
a
D
x
γ
f
(
x
)
≈
a
Δ
x
γ
f
(
x
)
=
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
-
r
h
+
h
)
,
(18)
x
D
a
γ
f
(
x
)
≈
x
Δ
a
γ
f
(
x
)
=
h
-
γ
∑
r
=
0
n
‍
(
-
1
)
r
(
γ
r
)
f
(
x
+
r
h
-
h
)
,
which give the first-order accurate approximations.


## 3.3. First-Order Approximation for the Integer-Order Derivatives

There is a first-order space derivative∂
V
(
x
,
t
)
/
∂
x in the LCP (8) to be approximated and we have many choices to use to obtain approximation of order 1 or the higher.For first-order approximation, we have the forward difference(19)
∂
V
i
,
j
∂
x
=
V
i
+
1
,
j
-
V
i
,
j
h
+
O
(
h
)
or the backward difference
(20)
∂
V
i
,
j
∂
x
=
V
i
,
j
-
V
i
-
1
,
j
h
+
O
(
h
)
.In this paper, we choose the backward difference method.


## 3.4. Crank-Nicolson Discretization

AssumingR is the matrix arising from the finite difference scheme, a general time stepping scheme for pricing American options is
(21)
V
i
,
j
-
V
i
,
j
+
1
Δ
t
+
Θ
R
V
i
,
j
+
(
1
-
Θ
)
R
V
i
,
j
+
1
=
0
,
x
i
>
f
(
t
j
)
,
V
i
,
j
=
K
-
e
x
i
,
x
i
≤
f
(
t
j
)
,
where j
=
0,1
,
…
,
m
-
1 and f
(
t
j
) is the free boundary at time step t
j. For above Θ
=
1 we recover the explicit forward Euler scheme. Θ
=
0 gives the implicit backward Euler method. We have the Crank-Nicolson method when Θ
=
1
/
2.


## 3.5. Matrix Formulation

Using (18), (20), and (21), we can now have the discretization of the full FPDE (7):
(22)
V
i
,
j
-
V
i
,
j
+
1
Δ
t
+
Θ
[
A
V
i
,
j
-
V
i
-
1
,
j
h
+
B
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
b
i
+
r
-
1
V
i
+
r
-
1
,
j
000000000
+
C
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
c
i
-
r
+
1
V
i
-
r
+
1
,
j
+
D
V
i
,
j
]
+
(
1
-
Θ
)
[
A
V
i
,
j
+
1
-
V
i
-
1
,
j
+
1
h
∑
r
=
0
i
+
1
‍
00000.00000000
+
B
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
b
i
+
r
-
1
V
i
+
r
-
1
,
j
+
1
000000.0000000
+
C
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
c
i
-
r
+
1
V
i
-
r
+
1
,
j
+
1
+
D
V
i
,
j
+
1
]
,
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
v
.
v
v
v
i
=
1,2
,
3
,
…
,
n
-
1
,
where
(23)
w
r
=
1
2
π
∫
0
2
π
‍
(
1
-
e
-
i
ψ
)
γ
e
i
r
ψ
d
ψ
,
which can be efficiently computed by FFT.Define matricesG, E, and F such that
(24)
[
G
V
j
]
i
=
(
A
*
Δ
t
h
+
D
*
Δ
t
)
V
i
,
j
-
A
*
Δ
t
h
V
i
-
1
,
j
,
[
E
V
j
]
i
=
Δ
t
B
i
h
-
γ
∑
r
=
0
i
+
1
‍
w
r
b
i
-
r
+
1
V
i
-
r
+
1
,
j
+
1
,
[
F
V
j
]
i
=
Δ
t
C
i
h
-
γ
∑
r
=
0
n
-
i
+
1
‍
w
r
c
i
+
r
-
1
V
i
+
r
-
1
,
j
+
1
.
Equation (22) then can be rewritten in matrix form as follows:
(25)
[
I
-
(
1
-
Θ
)
G
]
V
j
+
1
=
[
I
+
Θ
G
]
V
j
+
(
1
-
Θ
)
[
E
+
F
]
V
j
+
1
+
Θ
[
E
+
F
]
V
j
.IfΘ is set to 1
/
2, this fractional Crank-Nicholson discretization approach will produce a method with a local truncation error, that is, O
(
Δ
t
2
)
+
O
(
h
).


## 3.6. Spatially Second-Order Approximation Obtained by Extrapolation

In numerical analysis, extrapolation is used to improve the rate of convergence of a sequence. Based on the spatially first-order accurate numerical approximation in 3.5, extrapolation can be employed to improve the low order of spatial convergence. First, we state below the result obtained in [28].Proposition 1.
Let1
≤
α
≤
2, f
∈
C
n
+
3
(
R
) such that all derivatives of f up to order n
+
3 belong to L
1
(
R
). For any integer p
≥
0, define the shifted Grünwald operator by
(26)
Δ
h
,
p
α
f
(
x
)
=
∑
j
=
0
∞
‍
(
-
1
)
j
(
α
j
)
f
(
x
-
(
j
-
p
)
h
)
.
Then, one has for some constants a
l independent of h, f, and x the fact that
(27)
h
-
α
Δ
h
,
p
α
f
(
x
)
=
d
α
d
x
α
f
(
x
)
+
∑
l
=
1
n
-
1
‍
(
a
l
d
α
+
l
d
x
α
+
l
f
(
x
)
)
h
l
+
O
(
h
n
)
uniformly in x
∈
R.
Writing (27) as
(28)
h
-
α
Δ
h
,
p
α
f
(
x
)
=
d
α
d
x
α
f
(
x
)
+
c
1
h
+
c
2
h
2
+
O
(
h
3
)
,
where c
1 and c
2 are independent of h, and (28) shows that the truncation error in the shifted Grünwald finite difference approximation is c
1
h
+
c
2
h
2
+
O
(
h
3
). Now, one can use the Richardson extrapolation method in the following way. First, the Crank-Nicholson method is applied for grid size h and again for grid size h
/
2. This gives two solutions of first-order accuracy for the spatial dimension. We denote by V
h the solution of grid size h and by V
h
/
2 the solution of gird size h
/
2. Then, one has
(29)
V
h
=
V
+
b
1
h
+
b
2
h
2
+
O
(
h
3
)
,
V
h
/
2
=
V
+
b
1
h
2
+
b
2
h
2
4
+
O
(
h
3
)
,
where b
1 and b
2 are real numbers. Taking V
=
2
*
V
h
/
2
-
V
h, it is obvious that V is spatially 2nd-order accurate.
Based on the classical Crank-Nicholson method combined with spatial extrapolation, we obtain a solution with local truncation errorO
(
Δ
t
2
)
+
O
(
h
2
).The extrapolation method can be used at each time step before the maturity of options; because we can obtain two solutions of first-order estimates for diffenent spatial grid size at each time step, then second-order approximation can be computed at each time step.


## 4. The Penalty Method

In order to apply penalty method, we need to replace the following equation in LCP:(30)
L
V
×
[
V
-
g
(
x
)
]
=
0
with
(31)
L
V
=
∂
V
(
x
,
t
)
∂
t
+
A
∂
V
(
x
,
t
)
∂
x
+
B
(
x
)
x
D
∞
γ
(
b
(
x
)
V
(
x
,
t
)
)
+
C
(
x
)
-
∞
D
x
γ
(
c
(
x
)
V
(
x
,
t
)
)
+
D
V
(
x
,
t
)
=
ρ
max
⁡
(
V
(
x
)
*
-
V
(
x
,
t
)
,
0
)
,
where V
(
x
)
* is the payoff upon the exercise and in the limit as the penalty parameter ρ
→
∞, the solution satisfies V
(
x
,
t
)
≥
V
(
x
)
*. Let P be the matrix given by the above penalty item. Consider the following:
(32)
P
(
V
j
+
1
)
i
i
=
{
large
number
when
V
i
,
j
+
1
<
V
i
*
0
other
case
.
Then the matrix form for the penalty method is obtained using ((25), (31), (32)). Consider the following:
(33)
[
I
-
(
1
-
Θ
)
G
+
P
(
V
j
+
1
)
]
V
j
+
1
=
[
I
+
Θ
G
]
V
j
+
(
1
-
Θ
)
[
E
+
F
]
V
j
+
1
+
Θ
[
E
+
F
]
V
j
+
P
(
V
j
+
1
)
V
*
.


### 4.1. Iterative Algorithm

At each time step, we need to solve (33) to obtain the valuation of American options at grid points of x
i
>
f
(
t
j
). Option values at grid points of x
i
≤
f
(
t
j
) are given by payoff function g
(
x
i
). Iterative methods are introduced in Algorithm 1.Algorithm 1:Fixed point iteration at each timestep.
(1)Let (
V
j
+
1
)
0
=
V
j
(2)Let V
^
k
=
(
V
j
+
1
)
k
(3)Let P
^
k
=
P
(
(
V
j
+
1
)
k
)
(4)For k
=
0,1
,
2
,
…
,
u
n
t
i
l
c
o
n
v
e
r
g
e
n
c
e
d
o
(5) Solve [
I
-
(
1
/
2
)
G
+
P
^
k
]
V
^
k
+
1
=
(
1
/
2
)
[
G
+
E
+
F
]
V
j
+
(
1
/
2
)
[
E
+
F
]
V
^
k
+
P
^
k
V
*
(6) if  ma
x
i
(
|
V
^
k
+
1
-
V
^
k
|
/
max
⁡
(
1
,
|
V
^
k
+
1
|
)
)
<
t
o
l
e
r
a
n
c
e
t
h
e
n
(7)  Quit
(8) end if
(9)end for


## 4.1. Iterative Algorithm

At each time step, we need to solve (33) to obtain the valuation of American options at grid points of x
i
>
f
(
t
j
). Option values at grid points of x
i
≤
f
(
t
j
) are given by payoff function g
(
x
i
). Iterative methods are introduced in Algorithm 1.Algorithm 1:Fixed point iteration at each timestep.
(1)Let (
V
j
+
1
)
0
=
V
j
(2)Let V
^
k
=
(
V
j
+
1
)
k
(3)Let P
^
k
=
P
(
(
V
j
+
1
)
k
)
(4)For k
=
0,1
,
2
,
…
,
u
n
t
i
l
c
o
n
v
e
r
g
e
n
c
e
d
o
(5) Solve [
I
-
(
1
/
2
)
G
+
P
^
k
]
V
^
k
+
1
=
(
1
/
2
)
[
G
+
E
+
F
]
V
j
+
(
1
/
2
)
[
E
+
F
]
V
^
k
+
P
^
k
V
*
(6) if  ma
x
i
(
|
V
^
k
+
1
-
V
^
k
|
/
max
⁡
(
1
,
|
V
^
k
+
1
|
)
)
<
t
o
l
e
r
a
n
c
e
t
h
e
n
(7)  Quit
(8) end if
(9)end for


## 5. Numerical Results

In this section, we consider an American put option. A series of tests were carried out to study the effectiveness of our method under KoBol model.In Table1, column “change” shows the convergence of our approximations along with the temporal and spatial grid getting finer. Convergence rates of second-order estimates obtained by extrapolation V
=
2
V
h
/
2
-
V
h are significantly improved, compared with those of first-order approximations. Hence, the test results confirm the effectiveness of our method.Table 1
First-order numerical estimates computed by iterative algorithm and second-order estimates obtained by extrapolation.
m
First-order approximationV
h
First-order approximationV
h
/
2
Second-orderV
=
2
V
h
/
2
-
V
h
n
itns
Value
Change
n
itns
Value
Change
Value
Change
25
128
5.92
0.5386
256
8.00
0.5267
0.5148
50
256
6.04
0.5275
0.0111
512
8.94
0.5199
0.0068
0.5123
0.0025
100
512
6.18
0.5203
0.0072
1024
9.90
0.5167
0.0032
0.5131
0.0008
200
1024
6.63
0.5169
0.0034
2048
11.03
0.5149
0.0018
0.5129
0.0002
400
2048
7.20
0.5150
0.0019
4096
12.24
0.5139
0.0010
0.5128
0.0001
“itns” is the average number of iterations per time step.m is the number of time steps. n is the number of spatial points. “value” is a sample of the put option values.Figure1 is an example of first- and second-order approximations under same model conditions.Figure 1
Option values for first-order approximations computed by iterative algorithm and second-order estimates obtained by extrapolation withT
=
0.5, λ
=
2, α
=
1.5, σ
=
0.32, K
=
e
2, and r
=
0.08.Figure2 gives the evolution of pricing process for second-order approximation at all time steps.Figure 2
Option values for second-order approximations at all timesteps withm
=
25, n
=
128, T
=
0.5, λ
=
2, α
=
1.5, σ
=
0.32, K
=
e
2, and r
=
0.08.


## 6. Conclusions

In this paper, starting from the LCP involving FPDE, we first present the algorithms of first-order approximations for American options under fractional diffusion models. Then, on the basis of the first-order accurate algorithm, spatial extrapolation is employed to obtain second-order accurate numerical estimates. Numerical results in Section6 prove that the algorithm is feasible and the extrapolation method is effective.Using the Richardson extrapolation method, we even can obtain numerical approximations of higher order. Hence, the algorithm and extrapolation method form a useful tool box for the engineer. The use of fractional derivatives discretization conveniently avoids the singularity problem in the integral part of the PIDE. With the development of the fractional derivatives, we could see more possibilities of its application in finance.Finally, as for American options, how to extend our method to other models and get combined with other existing methods might be interesting topics. In order to improve computation performance, new technology like parallel computing or new iterative algorithm having quicker convergence rate are important in further study.



---

*Source: 364868-2014-09-29.xml*