# The Pricing of Asian Options in Uncertain Volatility Model

Authors: Yulian Fan; Huadong Zhang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper studies the pricing of Asian options when the volatility of the underlying
asset is uncertain. We use the nonlinear Feynman-Kac formula in the G-expectation theory
to get the two-dimensional nonlinear PDEs. For the arithmetic average fixed strike Asian
options, the nonlinear PDEs can be transferred to linear PDEs. For the arithmetic average
floating strike Asian options, we use a dimension reduction technique to transfer the two-dimensional
nonlinear PDEs to one-dimensional nonlinear PDEs. Then we introduce the
applicable numerical computation methods for these two classes of PDEs and analyze the
performance of the numerical algorithms.


---

## Body

## 1. Introduction

Model uncertainty aversion has shown to have important consequences in price behavior in capital markets [1–3] and macroeconomics [4, 5]. Denis and Martini [6] prove, in the model uncertainty situation the value of any bounded contingent claim can be expressed as the supremum of its expectation over some set of martingale measures.In finance, there is an important case of model uncertainty called volatility uncertainty in which the uncertainty comes from the volatility coefficient. A major difficulty here is that the probabilities are mutually singular. This type of uncertainty is initially studied by Avellaneda et al. [7] and Lyons [8] for the superhedging of European options with payoffs depending only on the basic assets' terminal values. In this case the values of the derivatives should satisfy the (Black-Scholes-Barenblatt) BSB equations. The corresponding discrete-time case has been studied in Delbaen [9]. But for the path-dependent options, the difficulty is dramatically increased. In this paper, we are going to study the pricing of Asian options in the uncertain volatility model. As far as we know, there is little literature on this problem.In the model without volatility uncertainty, there is a vast literature on pricing of Asian options.In the geometric Brownian motion model, Asian options pricing methods are studied extensively. Geman and Yor [10] develop a dimensional Laplace transform method to price Asian options (see also [11]). Kemna and Vorst [12] use Monte Carlo simulation with a specific variance reduction method to compute the prices of fixed strike average-rate options. Ingersoll [13], and Wilmott et al. [14] prove that the two-dimensional PDE for a floating strike Asian option can be reduced to a one-dimensional PDE. Rogers and Shi [15] formulate a one-dimensional PDE that can model the prices of both floating and fixed strike Asian options and propose a numerical technique to compute tight lower bound on European average-rate options' prices. Barraquand and Pudet [16] describe a forward shooting grid algorithm and prove that it is unconditionally convergent. In Zvan et al. [17] and Zvan et al. [18], a flux limiter is used to retain accuracy while preventing oscillations. Simon et al. [19] propose an easy computable upper bound for the prices of arithmetic Asian options. Kaas et al. [20] and Dhaene et al. [21] compute the sharp lower and upper bounds for the prices of arithmetic Asian options by using the method of “comonotonic approximations.” Večeř [22], and Zhang [23, 24] propose various PDE methods. In Marcozzi [25], the first-order hyperbolic term is discretized by using a first-order upwind type method, resulting in at most first-order accuracy. A related approach based on a combination of the (weighted essentially nonoscillatory) WENO discretization and the grid stretching is used for Asian options in Oosterlee et al. [26]. Linetsky [27] investigates a spectral expansion approach.In the jump diffusion models, methods for pricing Asian options are still under development. Fusai [28] obtains a simple expression for the double transform of the prices of continuously monitored Asian options by means of Fourier and Laplace transform. D’Halluin et al. [29] propose a semi-Lagrangian method to price the continuously observed fixed strike Asian options. Cai and Kou [30] consider the pricing of Asian options for the double exponential jump diffusion model. Bayraktar and Xing [31] obtain a fast numerical approximation scheme and analyze the performance of the numerical algorithm.In this paper, we use the stochastic processes driven byG-Brownian motion to describe the stocks' price processes. G-Brownian motion is introduced in Peng [32–34], and it has the properties that its increments are zero-mean, independent, and stationary and can be proved to be G-normally distributed. The G-normal distribution random variables have no mean uncertainty but can have variance uncertainty. The variance uncertainty of the G-normal distribution is the source of the volatility uncertainty of the stocks' price processes.The paper is organized as follows. In Section2, we use the nonlinear Feynman-Kac formula in the G-expectation theory to get the two-dimensional nonlinear PDEs satisfied by the values of the Asian options. For the arithmetic average fixed strike Asian options, we prove the convexity of the options' values with respect to the stock prices and transfer the nonlinear PDEs to linear PDEs with the uncertain volatilities being replaced by the maximum volatilities. For the arithmetic average floating strike Asian options, we use a dimension reduction technique to transfer the two-dimensional nonlinear PDEs to one-dimensional nonlinear PDEs. In Section 3, we introduce the applicable numerical computation methods for these two classes of PDEs and give the numerical results.


## 2. The Pricing of Asian Options in Uncertain Volatility Model

### 2.1.G-Brownian Motion and G-Expectation

In this subsection, we will introduce theG-Brownian motion and G-expectation defined in Peng [33–36].LetΩ be a given set and let H be a linear space of real valued functions defined on Ω. We suppose that H satisfies c
∈
H for each constant c and |
X
|
∈
H if X
∈
H. The space H can be considered as the space of random variables.A sublinear expectationE on H is a functional E
:
H
→
R satisfying the following properties: for all X
,
Y
∈
H, we have(a)
Monotonicity: IfX
≥
Y then E
[
X
]
≥
E
[
Y
].
(b)
Constant preserving:E
[
c
]
=
c, for c
∈
R.
(c)
Sub-additivity:E
[
X
+
Y
]
≤
E
[
X
]
+
E
[
Y
].
(d)
Positive homogeneity:E
[
λ
X
]
=
λ
E
[
X
]
,
∀
λ
≥
0.The triple(
Ω
,
H
,
E
) is called a sublinear expectation space.Definition 1.
Ad-dimensional process (
B
t
)
t
≥
0 on a sublinear expectation space (
Ω
,
H
,
E
) is called a G-Brownian motion if the following properties are satisfied:(i)
B
0
(
ω
)
=
0;
(ii)
For eacht
,
s
≥
0, the increment B
t
+
s
-
B
t is N
(
{
0
}
×
Σ
)-distributed and is independent from (
B
t
1
,
B
t
2
,
…
,
B
t
n
), for each n
∈
N and 0
≤
t
1
≤
⋯
≤
t
n
≤
t.Remark 2.
Just as in the classical situation, the increments ofG-Brownian motion (
B
t
+
s
-
B
t
)
s
≥
0 are independent from Ω
t.Letπ
t
N
=
{
t
0
N
,
t
1
N
,
…
,
t
N
N
} such that 0
=
t
0
N
<
t
1
N
<
⋯
<
t
N
N
=
t for N
=
1,2
,
…, be a sequence of partitions of [
0
,
t
] with
(1)
μ
(
π
t
N
)
∶
=
max
⁡
⁡
{
|
t
i
+
1
-
t
i
|
:
i
=
0,1
,
…
,
N
-
1
}
.
The quadratic variation process of 1-dimensional G-Brownian motion (
B
t
)
t
≥
0 with B
1
=

=
d
N
(
{
0
}
×
[
σ
_
2
,
σ
¯
2
]
) is defined as
(2)
〈
B
〉
t
∶
=
lim
⁡
μ
(
π
t
N
)
→
0
⁡
∑
j
=
0
N
-
1
‍
(
B
t
j
+
1
N
-
B
t
j
N
)
2
=
B
t
2
-
2
∫
0
t
B
s
d
B
s
.
(
〈
B
〉
t
)
t
≥
0 is an increasing process with 〈
B
〉
0
=
0.LetΩ
=
C
0
d
(
R
+
) be the space of all R
d-valued continuous paths (
ω
t
)
t
∈
R
+, with ω
0
=
0, equipped with the distance
(3)
ρ
(
ω
1
,
ω
2
)
∶
=
∑
i
=
1
∞
‍
2
-
i
[
(
max
⁡
t
∈
[
0
,
i
]
|
ω
t
1
-
ω
t
2
|
)
∧
1
]
.
For each fixed T
∈
[
0
,
∞
), we set Ω
T
:
=
{
ω
·
∧
T
:
ω
∈
Ω
}.From Peng [36], the distribution of the quadratic variation process of the G-Brownian motion has the following properties.Lemma 3.
For each fixeds
,
t
≥
0
,
〈
B
〉
s
+
t
-
〈
B
〉
s is identically distributed with 〈
B
〉
t and independent from Ω
s.Theorem 4.
〈
B
〉
t is N
(
[
σ
_
2
t
,
σ
¯
2
t
]
×
{
0
}
)-distributed; that is, for each φ
∈
C
l
,
L
i
p
(
R
),
(4)
E
^
[
φ
(
〈
B
〉
t
)
]
=
sup
⁡
σ
_
2
≤
v
≤
σ
¯
2
φ
(
v
t
)
.The sublinear expectationE
^ in the above theorem is called G-expectation defined on L
i
p
(
Ω
) via the following procedure. Let
(5)
L
i
p
(
Ω
T
)
∶
=
{
ϕ
(
B
t
1
∧
T
,
…
,
B
t
n
∧
T
)
:
n
∈
N
,
k
t
1
,
…
,
t
n
∈
[
0
,
∞
)
,
ϕ
∈
C
l
,
L
i
p
(
R
d
×
n
)
}
with B
t
=
ω
t
,
t
∈
[
0
,
∞
) for ω
∈
Ω, and
(6)
L
i
p
(
Ω
)
∶
=
⋃
n
=
1
∞
‍
L
i
p
(
Ω
n
)
.Firstly, construct a sequence ofd-dimensional random vectors (
ξ
)
i
=
1
∞ on a sublinear expectation space (
Ω
~
,
H
~
,
E
~
) such that ξ
i is G-normal distributed and ξ
i
+
1 is independent from (
ξ
1
,
…
,
ξ
i
) for each i
=
1,2
,
…. For each X
∈
L
i
p
(
Ω
) with
(7)
X
=
ϕ
(
B
t
1
-
B
t
0
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
,
some ϕ
∈
C
l
,
L
i
p
(
R
d
×
n
) and 0
=
t
0
<
t
1
<
⋯
<
t
n
<
∞, set
(8)
E
^
[
ϕ
(
B
t
1
-
B
t
0
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
]
∶
=
E
~
[
ϕ
(
t
1
-
t
0
ξ
1
,
…
,
t
n
-
t
n
-
1
ξ
n
)
]
.And the related conditional expectation of(9)
X
=
ϕ
(
B
t
1
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
under Ω
t
j is defined by
(10)
E
^
[
X
∣
Ω
t
j
]
∶
=
ψ
(
B
t
1
,
…
,
B
t
j
-
B
t
j
-
1
)
,
where
(11)
ψ
(
x
1
,
…
,
x
j
)
=
E
~
[
ϕ
(
x
1
,
…
,
x
j
,
t
j
+
1
-
t
j
ξ
j
+
1
,
…
,
t
n
-
t
n
-
1
ξ
n
)
]
.
E
^
[
·
] consistently defines a sublinear expectation on L
i
p
(
Ω
) and (
B
t
)
t
≥
0 is a G-Brownian motion.We letL
G
p
(
Ω
t
k
)
,
p
≥
1, denote the completion of
(12)
L
i
p
(
Ω
t
k
)
∶
=
{
φ
(
B
t
1
∧
T
,
…
,
B
t
n
∧
T
)
:
n
∈
N
,
k
t
1
,
…
,
t
n
∈
[
0
,
∞
)
,
φ
∈
C
l
,
L
i
p
(
R
d
×
n
)
}
under the norm ∥
X
∥
p
∶
=
(
E
^
[
|
X
|
p
]
)
1
/
p.


### 2.2. The PDEs

We assume the security market is defined on a sublinear expectation space(
Ω
,
H
,
E
^
) which opens continuously and offers a constant riskless interest rate r to all borrowers and lenders. And we also assume there are no transaction costs and/or taxes. Suppose that the price process (
S
t
)
t
≥
0 of some risky asset is given by
(13)
d
S
ν
=
μ
S
ν
d
ν
+
σ
~
S
ν
d
B
ν
,
S
0
=
S
,
where B is the G-Brownian motion with σ
_
2
=
-
E
^
[
-
B
1
2
] and σ
¯
2
=
E
^
[
B
1
2
], and μ and σ
~ are constants. Hence (
S
ν
)
ν
≥
0 is a geometric G-Brownian motion:
(14)
S
ν
=
S
t
exp
⁡
(
1
2
μ
(
ν
-
t
)
k
k
k
k
k
i
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
)
.Define(15)
I
ν
=
∫
0
ν
‍
S
τ
d
τ
.
Then I
ν satisfies the following stochastic differential equation:
(16)
d
I
ν
=
S
ν
d
ν
,
I
0
=
0
.LetX
ν
=
(
S
ν
,
I
ν
)
T, so
(17)
d
X
ν
=
μ
X
X
ν
d
ν
+
σ
~
X
X
ν
d
B
ν
,
X
0
=
x
=
(
S
,
0
)
,
where μ
X and σ
~
X are 2
×
2 matrixes:
(18)
μ
X
=
(
μ
0
1
0
)
,
σ
~
X
=
(
σ
~
0
0
0
)
.The typical Asian options in this security market are the following four types:the arithmetic average fixed strike call option with final payoff,(19)
max
⁡
(
I
T
T
-
K
,
0
)
;
the corresponding arithmetic average fixed strike put option with payoff,(20)
max
⁡
(
K
-
I
T
T
,
0
)
;
the arithmetic average floating strike call option with final payoff,(21)
max
⁡
(
S
T
-
I
T
T
,
0
)
;
and the corresponding arithmetic average floating strike put option with payoff,(22)
max
⁡
(
I
T
T
-
S
T
,
0
)
.LetV
(
t
,
S
t
,
I
t
) denote the price of one of the options at time t. Take the arithmetic average fixed strike call option for example, we know V
(
T
,
S
T
,
I
T
)
=
max
⁡
(
I
T
/
T
-
K
,
0
).Here we still assume the market operates in a risk neutral world [7].Peng [35] introduces the axiomatic conditions satisfied by the valuation mechanism. In the sublinear expectation space (
Ω
,
H
,
E
^
), E
^
t
[
·
] is a filtration consistent valuation operator. By this valuation mechanism, we have
(23)
V
(
t
,
S
t
,
I
t
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
∣
Ω
t
]
.
Here we assume the riskless interest rate is r.Since we assume the market is risk neutral, we haveμ
=
r. It is easy to verify that the discounted prices of the stock satisfy the valuation mechanism {
E
^
t
[
·
]
}
t
≥
0; that is,
(24)
S
t
=
E
^
t
[
e
-
r
(
ν
-
t
)
S
ν
]
,
0
≤
t
≤
ν
≤
T
.Let us assume(25)
d
X
ν
t
,
x
=
μ
X
X
ν
t
,
x
d
ν
+
σ
~
X
X
ν
t
,
x
d
B
ν
,
X
t
=
x
=
(
S
,
I
)
,
with
(26)
μ
X
=
(
r
0
1
0
)
,
σ
~
X
=
(
σ
~
0
0
0
)
.The following result is a corollary of the nonlinear Feynman-Kac formula of Peng [36].Corollary 5.
Let(27)
Y
s
t
,
x
=
E
^
[
∫
s
T
‍
e
-
r
(
T
-
s
)
Φ
(
X
T
t
,
x
)
k
k
+
∫
s
T
‍
e
-
r
(
ν
-
s
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
k
k
+
∫
s
T
‍
e
-
r
(
ν
-
s
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
∣
Ω
s
]
,
where Φ
:
R
2
→
R is a given Lipschitz function and f
,
g
:
R
2
×
R
→
R are given Lipschitz functions; that is, |
ϕ
(
x
,
y
)
-
ϕ
(
x
′
,
y
′
)
|
≤
K
(
|
x
-
x
′
|
+
|
y
-
y
′
|
), for each x
,
x
′
∈
R
2,  y
,
y
′
∈
R,  ϕ
=
f and g. And denote
(28)
u
(
t
,
x
)
=
Y
t
t
,
x
.
Then
(29)
u
(
t
,
x
)
=
E
^
[
∫
t
t
+
δ
‍
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
]
is a viscosity solution of the following PDE:
(30)
∂
t
u
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
u
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
+
g
(
x
,
u
)
)
σ
2
+
〈
D
u
(
t
,
x
)
,
μ
X
x
〉
-
r
u
+
f
(
x
,
u
)
=
0
,
u
(
T
,
x
)
=
Φ
(
x
)
.Proof.
The proof of expression (29) is the same as that of Peng [36]. By the proof of Peng's theorem of the nonlinear Feynman-Kac formula and the boundedness of e
-
r
(
T
-
t
), we can easily prove that u is a Lipschitz function; that is,
(31)
|
u
(
t
,
x
)
-
u
(
t
,
x
′
)
|
≤
C
|
x
-
x
′
|
,
(32)
E
^
[
|
Y
s
t
,
x
|
]
≤
C
(
1
+
|
x
|
)
.
Let us see the continuity of u with respect to t. Firstly, we have
(33)
|
u
(
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
≤
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
+
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kkk
k
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
|
]
.
Replace e
-
r
δ with its Taylor's expansion in the first term on the right side of (33), we have
(34)
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
=
E
^
[
|
(
1
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
k
k
k
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
=
E
^
[
|
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
k
k
k
+
(
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
|
]
≤
E
^
[
|
X
t
+
δ
t
,
x
-
x
|
]
+
E
^
[
|
(
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
|
]
.
By the following formulas [36]
(35)
E
^
[
|
X
t
+
δ
t
,
x
-
x
|
2
]
≤
C
(
1
+
x
2
)
δ
,
E
^
[
|
X
t
+
δ
t
,
x
|
2
]
≤
C
(
1
+
x
2
)
,
and (32), we have
(36)
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
≤
C
′
(
1
+
|
x
|
2
)
1
/
2
(
δ
1
/
2
+
δ
+
δ
2
2
⋯
)
.
For the second term on the right side of (33), we have
(37)
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
|
]
=
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
kkkik
k
k
k
k
×
(
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
-
f
(
x
,
Y
t
t
,
x
)
)
d
ν
kkkik
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
kkkik
k
k
k
k
k
i
×
(
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
-
g
(
x
,
Y
t
t
,
x
)
)
d
〈
B
〉
ν
∫
t
t
+
δ
‍
|
]
+
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
x
,
Y
t
t
,
x
)
d
ν
kkkik
ki
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
x
,
Y
t
t
,
x
)
d
〈
B
〉
ν
|
]
≤
C
(
1
+
|
x
|
)
δ
.
Therefore u is continuous in t.
Now we are going to prove thatu
(
t
,
x
) is a viscosity solution of PDE (30). For fixed (
t
,
x
)
∈
(
0
,
T
)
×
R
2, let ψ
∈
C
b
2,3
(
(
0
,
T
)
×
R
2
) such that ψ
≥
u and ψ
(
t
,
x
)
=
u
(
t
,
x
). By (29) and Taylor's expansion, for δ
∈
(
0
,
T
-
t
),
(38)
0
≤
E
^
[
e
-
r
δ
ψ
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
ψ
(
t
,
x
)
∫
t
t
+
δ
‍
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
]
≤
E
^
[
1
2
-
r
ψ
(
t
,
x
)
δ
+
∂
t
ψ
(
t
,
x
)
δ
kkk
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
δ
kkk
+
〈
1
2
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
(
〈
B
〉
t
+
δ
-
〈
B
〉
t
)
kkk
1
2
+
f
(
x
,
ψ
)
δ
+
g
(
x
,
ψ
)
(
〈
B
〉
t
+
δ
-
〈
B
〉
t
)
]
+
C
(
1
+
|
x
|
+
|
x
|
2
+
|
x
|
3
)
δ
3
/
2
≤
(
∂
t
ψ
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
kkkkkkkkkkkkkk
+
g
(
x
,
ψ
)
1
2
)
σ
2
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
k
i
-
r
ψ
(
t
,
x
)
+
f
(
x
,
ψ
)
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
)
δ
+
C
(
1
+
|
x
|
+
|
x
|
2
+
|
x
|
3
)
δ
3
/
2
.
It is easy to check that
(39)
∂
t
ψ
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
+
g
(
x
,
ψ
)
)
σ
2
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
-
r
ψ
(
t
,
x
)
+
f
(
x
,
ψ
)
≥
0
.
Thus u is a viscosity subsolution of (30). Similarly we can prove that u is a viscosity supersolution of (30).
For the arithmetic average Asian options, we have(40)
V
(
t
,
x
)
=
V
(
t
,
S
,
I
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
]
.
Then by Corollary 5, V
(
t
,
x
) satisfies the following PDE:
(41)
∂
t
V
(
t
,
x
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
〈
D
2
V
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
σ
2
+
〈
D
V
(
t
,
x
)
,
μ
X
x
〉
-
r
V
(
t
,
x
)
=
0
,
V
(
T
,
x
)
=
φ
(
T
,
x
)
.
By (18), the above PDE is
(42)
∂
t
V
(
t
,
S
,
I
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
S
2
σ
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
φ
(
T
,
S
,
I
)
,
where
(43)
φ
(
T
,
S
,
I
)
=
{
max
⁡
(
I
T
-
K
,
0
)
,
for
the
arithmetic
average
fixed
strike
call
option
,
max
⁡
(
K
-
I
T
,
0
)
,
for
the
arithmetic
average
fixed
strike
put
option
,
max
⁡
(
S
-
I
T
,
0
)
,
for
the
arithmetic
average
floating
strike
call
option
,
max
⁡
(
I
T
-
S
,
0
)
,
for
the
arithmetic
average
floating
strike
put
option
.


### 2.3. The Arithmetic Average Fixed Strike Asian Call Option

Let us see the arithmetic average fixed strike Asian call option. From (14), we have
(44)
S
ν
t
,
S
=
S
exp
⁡
⁡
{
μ
(
ν
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
}
.
Then by the definition of I
ν, we have
(45)
I
ν
t
,
I
=
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
k
k
k
k
k
k
k
k
k
k
+
σ
~
(
B
τ
-
B
t
)
{
μ
(
τ
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
}
d
τ
.
Since
(46)
V
(
t
,
x
)
=
V
(
t
,
S
,
I
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
]
=
E
^
[
e
-
r
(
T
-
t
)
(
I
T
T
-
K
)
+
]
,
V
(
t
,
S
,
I
)
=
E
^
[
×
(
T
)
-
1
(
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
)
-
K
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
)
+
e
-
r
(
T
-
t
)
i
i
i
k
×
(
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
1
2
kkkkkkkkkkkkkkki
i
i
i
k
i
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
kkkkkkkkkkkkkkki
i
i
i
i
k
1
2
+
σ
~
(
B
τ
-
B
t
)
}
d
τ
S
∫
t
ν
‍
)
i
i
k
i
×
T
-
1
(
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
)
-
K
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
)
+
]
.
By the convexities of function f
(
x
)
=
(
x
-
k
)
+ and the sublinear expectation, we can say
(47)
V
(
t
,
λ
S
1
+
(
1
-
λ
)
S
2
,
I
)
≤
λ
V
(
t
,
S
1
,
I
)
+
(
1
-
λ
)
V
(
t
,
S
2
,
I
)
;
that is, V
(
t
,
S
,
I
) is convex about S. So in the numerical computation, the second order difference of V is positive. Then the numerical solutions of (43) are equivalent to those of PDE:
(48)
∂
t
V
(
t
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
(
I
T
-
K
)
+
.Theoretically, the above PDE is defined on domainD
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
∞
)
×
[
0
,
∞
)
}. But generally, closed form solutions cannot be found; we have to find numerical solutions. We have no interest in finding the solution in the entire infinite domain. Therefore we consider the domain D
*
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
S
*
)
×
[
0
,
I
*
)
}.Now we determine the boundary conditions of (48). If I
≥
K
T, that is, I
/
T
≥
K, for the positiveness of S
ν, we have I
T
/
T
≥
K. Therefore the final payoff on the option is certain to be positive. At time T, this payoff can be expressed as
(49)
I
T
T
-
K
=
I
T
-
K
+
1
T
∫
t
T
‍
S
ν
d
ν
.
This payoff can also be obtained by using the following self-financing duplicating portfolio strategy. Assume that an investor commits (
I
/
T
-
K
)
e
-
r
(
T
-
t
) into riskless bonds in order to secure the return promised by the first part of the right of (49), that is, (
I
/
T
-
K
) at time T. If the investor is to be certain of accruing the return promised by the second part of the right of (49), he is obliged to transfer a certain portion of stock (equal to (
1
/
T
)
e
-
r
(
T
-
ν
)
Δ
ν) to a riskless bond in every time interval (
ν
,
ν
+
Δ
ν
). Overall this strategy requires a sum equal to (
I
/
T
-
K
)
e
-
r
(
T
-
t
) plus the portion of a stock which can be expressed as follows:
(50)
∫
t
T
‍
1
T
e
-
r
(
T
-
ν
)
d
ν
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
.
Then the price of the option when I
t
≥
K
T must be equal to
(51)
V
(
t
,
S
,
I
)
=
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.The above discussion is similar to that in the model without volatility uncertainty. Equation (51) obviously satisfies (48). So we only need to find the solution of (48) in the domain D
*
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
S
*
)
×
[
0
,
K
T
)
}.On the boundaryI
=
K
T, the boundary condition from (51) can be written as
(52)
V
(
t
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.We note that ifS
t
=
0, (13) tells us that S
(
ν
) is also equal to zero for ν
∈
[
t
,
T
] and I
(
T
) is thus equal to I
(
t
). By (46), we can write the following boundary condition:
(53)
V
(
t
,
0
,
I
)
=
e
-
r
(
T
-
t
)
(
I
T
-
K
)
+
.We need another boundary condition atS
=
S
max
⁡
=
S
*. Hugger [37] gives the boundary condition at S
=
S
max
⁡ in the model without volatility uncertainty. By the same discussion, we can get the same boundary condition in uncertain volatility model:
(54)
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
S
max
⁡
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
0
}
.So, the values of the arithmetic average fixed strike Asian call option satisfy(55)
∂
t
V
(
t
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
(
I
T
-
K
)
+
,
V
(
t
,
0
,
I
)
=
e
-
r
(
T
-
t
)
(
I
T
-
K
)
+
,
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
S
max
⁡
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
0
}
,
V
(
t
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.For the arithmetic average fixed strike Asian put option with terminal payoffp
(
T
,
S
T
,
I
T
)
=
max
⁡
⁡
(
K
-
I
T
/
T
,
0
), we will value this option by the following put-call parity.It is easily to see that we have(56)
max
⁡
(
I
T
T
-
K
,
0
)
=
max
⁡
(
K
-
I
T
T
,
0
)
+
1
T
∫
0
T
‍
S
ν
d
ν
-
K
.
So
(57)
c
(
t
,
S
t
,
I
t
)
=
e
-
r
(
T
-
t
)
E
^
t
[
max
⁡
(
I
T
T
-
K
,
0
)
]
=
e
-
r
(
T
-
t
)
E
^
t
[
max
⁡
(
K
-
I
T
T
,
0
)
+
1
T
∫
0
T
‍
S
ν
d
ν
]
-
e
-
r
(
T
-
t
)
K
=
p
(
t
,
S
t
,
I
t
)
-
e
-
r
(
T
-
t
)
K
+
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
ν
d
ν
+
e
-
r
(
T
-
t
)
E
t
σ
¯
[
1
T
∫
t
T
‍
S
ν
d
ν
]
.e
-
r
(
T
-
t
)
E
t
σ
¯
[
(
1
/
T
)
∫
t
T
‍
S
ν
d
ν
] is the time t value of time T claim (
1
/
T
)
∫
t
T
‍
S
ν
d
ν, where the stock has certain volatility σ
~
σ
¯. To get claim (
1
/
T
)
∫
t
T
‍
S
ν
d
ν at T, let us see the replicating strategy: in any time interval (
ν
,
ν
+
Δ
ν
), transfer a certain portion (
1
/
T
)
e
-
r
(
T
-
ν
)
Δ
ν of stock S
ν to a riskless bond. This strategy ensures the final payoff of (
1
/
T
)
∫
t
T
‍
S
ν
d
ν. The stocks that should be transferred are
(58)
∫
t
T
‍
1
T
e
-
r
(
T
-
ν
)
d
ν
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
.
Then the time t value of (
1
/
T
)
∫
t
T
‍
S
ν
d
ν equals (
1
/
r
T
)
(
1
-
e
-
r
(
T
-
t
)
)
S
t. So the put-call parity can be written as
(59)
c
(
t
,
S
,
I
)
+
K
e
-
r
(
T
-
t
)
=
p
(
t
,
S
,
I
)
+
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
ν
d
ν
+
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
t
.


### 2.4. The Arithmetic Average Floating Strike Asian Call Option

For the arithmetic average floating strike Asian call (put) option, the PDE (43) still applies with terminal condition
(60)
V
(
T
,
S
,
I
)
=
(
S
-
I
T
)
+
(
V
(
T
,
S
,
I
)
=
(
I
T
-
S
)
+
)
.In Ingersoll [13] where there is no volatility uncertainty, due to the linear homogeneity of the PDE and its terminal condition, their PDE can be reduced to a one-dimensional PDE. In our uncertain volatility model, the PDE (43) is not linear homogeneous anymore. But we can still reduce PDE (43) to a one-dimensional PDE by the following argument.For the arithmetic average floating strike Asian call option, we know(61)
V
(
t
,
S
t
,
I
t
)
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
T
-
I
T
T
)
+
]
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
T
-
I
t
+
∫
t
T
‍
S
ν
d
ν
T
)
+
]
.Denote(62)
M
t
ν
=
exp
⁡
(
μ
(
ν
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
)
.
Then by (14), we have S
ν
=
S
t
M
t
ν and S
t
=
S
ν
M
ν
t. So
(63)
I
t
=
∫
0
t
‍
S
ν
d
ν
=
∫
0
t
‍
S
t
M
t
ν
d
ν
=
S
t
∫
0
t
‍
M
t
ν
d
ν
,
V
(
t
,
S
t
,
I
t
)
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
t
M
t
T
-
S
t
∫
0
t
‍
M
t
ν
d
ν
+
S
t
∫
t
T
‍
M
t
ν
d
ν
T
)
+
]
=
S
t
E
^
t
[
e
-
r
(
T
-
t
)
(
M
t
T
-
∫
0
t
‍
M
t
ν
d
ν
+
∫
t
T
‍
M
t
ν
d
ν
T
)
+
]
∶
=
S
t
H
(
t
,
R
t
)
with R
t
=
I
t
/
S
t.Then the functionH
(
t
,
R
) satisfies the following PDE:
(64)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
which is a one-dimensional PDE.For the European style Asian option, we must impose boundary conditions both atR
=
0 and as R
→
∞. The boundary condition as R
→
∞ is simple. The only way that R can tend to infinity is that S tends to zero. In this case the option will not be exercised, and so
(65)
lim
⁡
R
→
∞
H
(
t
,
R
)
=
0
.
For computational purpose, we truncate the asset region (
0
,
+
∞
) into (
0
,
R
max
⁡
), where R
max
⁡ is sufficiently large to ensure the accuracy of the solution. If R
=
0, the PDE degenerates into the first-order hyperbolic PDE
(66)
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
.SoH
(
t
,
R
) satisfies
(67)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
H
(
T
,
R
)
=
(
1
-
R
T
)
+
,
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
H
(
t
,
R
max
⁡
)
=
0
.By the same argument, for the arithmetic average floating strike Asian put option, its prices could be calculated by(68)
V
(
t
,
S
,
I
)
=
S
H
(
t
,
R
)
,
with H
(
t
,
R
) satisfying the following PDE:
(69)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
H
(
T
,
R
)
=
(
R
T
-
1
)
+
,
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
H
(
t
,
R
max
⁡
)
=
(
R
max
⁡
T
-
exp
⁡
{
-
r
(
T
-
t
)
}
)
.Generally, the put-call parity will not hold any more in the uncertain volatility model. Let us see the following argument. The mature payoff of a portfolio with one European arithmetic average floating strike Asian call holding long and one put holding short is(70)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
-
S
T
max
⁡
(
R
T
T
-
1,0
)
=
S
T
-
R
T
S
T
T
.
The value of this portfolio is equal to one consisting of one stock and a financial product whose payoff is -
R
T
S
T
/
T. The values of the financial product with final payoff -
R
T
S
T
/
T satisfy the following PDE:
(71)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
with terminal condition H
(
T
,
R
T
)
=
-
R
T
S
T
/
T. We seek a solution of (71) of the form
(72)
H
(
t
,
R
)
=
a
t
+
b
t
R
with a
T
=
0 and b
T
=
-
1
/
T. Substituting (72) into (71) and satisfying the boundary conditions, we find that
(73)
a
t
=
-
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
b
t
=
-
1
T
e
-
r
(
T
-
t
)
.
For
(74)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
=
S
T
max
⁡
(
R
T
T
-
1,0
)
+
S
T
-
R
T
S
T
T
,
by the valuation mechanism in space (
Ω
,
H
,
E
^
), we take the discounted G-conditional expectation on both sides and get
(75)
E
^
t
[
e
-
r
(
T
-
t
)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
]
=
e
-
r
(
T
-
t
)
E
^
t
[
S
T
max
⁡
(
R
T
T
-
1,0
)
+
S
T
-
R
T
S
T
T
]
≤
e
-
r
(
T
-
t
)
E
^
t
[
S
T
max
⁡
(
R
T
T
-
1,0
)
]
+
e
-
r
(
T
-
t
)
E
^
t
[
S
T
-
R
T
S
T
T
]
.
That is
(76)
c
(
t
,
S
,
I
)
≤
p
(
t
,
S
,
I
)
+
S
-
S
r
T
(
1
-
e
-
r
(
T
-
t
)
)
-
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
τ
d
τ
,
where c
(
t
,
S
,
I
) is the time t price of the arithmetic average floating strike Asian call option and p
(
t
,
S
,
I
) is the time t price of the arithmetic average floating strike Asian put option.


## 2.1.G-Brownian Motion and G-Expectation

In this subsection, we will introduce theG-Brownian motion and G-expectation defined in Peng [33–36].LetΩ be a given set and let H be a linear space of real valued functions defined on Ω. We suppose that H satisfies c
∈
H for each constant c and |
X
|
∈
H if X
∈
H. The space H can be considered as the space of random variables.A sublinear expectationE on H is a functional E
:
H
→
R satisfying the following properties: for all X
,
Y
∈
H, we have(a)
Monotonicity: IfX
≥
Y then E
[
X
]
≥
E
[
Y
].
(b)
Constant preserving:E
[
c
]
=
c, for c
∈
R.
(c)
Sub-additivity:E
[
X
+
Y
]
≤
E
[
X
]
+
E
[
Y
].
(d)
Positive homogeneity:E
[
λ
X
]
=
λ
E
[
X
]
,
∀
λ
≥
0.The triple(
Ω
,
H
,
E
) is called a sublinear expectation space.Definition 1.
Ad-dimensional process (
B
t
)
t
≥
0 on a sublinear expectation space (
Ω
,
H
,
E
) is called a G-Brownian motion if the following properties are satisfied:(i)
B
0
(
ω
)
=
0;
(ii)
For eacht
,
s
≥
0, the increment B
t
+
s
-
B
t is N
(
{
0
}
×
Σ
)-distributed and is independent from (
B
t
1
,
B
t
2
,
…
,
B
t
n
), for each n
∈
N and 0
≤
t
1
≤
⋯
≤
t
n
≤
t.Remark 2.
Just as in the classical situation, the increments ofG-Brownian motion (
B
t
+
s
-
B
t
)
s
≥
0 are independent from Ω
t.Letπ
t
N
=
{
t
0
N
,
t
1
N
,
…
,
t
N
N
} such that 0
=
t
0
N
<
t
1
N
<
⋯
<
t
N
N
=
t for N
=
1,2
,
…, be a sequence of partitions of [
0
,
t
] with
(1)
μ
(
π
t
N
)
∶
=
max
⁡
⁡
{
|
t
i
+
1
-
t
i
|
:
i
=
0,1
,
…
,
N
-
1
}
.
The quadratic variation process of 1-dimensional G-Brownian motion (
B
t
)
t
≥
0 with B
1
=

=
d
N
(
{
0
}
×
[
σ
_
2
,
σ
¯
2
]
) is defined as
(2)
〈
B
〉
t
∶
=
lim
⁡
μ
(
π
t
N
)
→
0
⁡
∑
j
=
0
N
-
1
‍
(
B
t
j
+
1
N
-
B
t
j
N
)
2
=
B
t
2
-
2
∫
0
t
B
s
d
B
s
.
(
〈
B
〉
t
)
t
≥
0 is an increasing process with 〈
B
〉
0
=
0.LetΩ
=
C
0
d
(
R
+
) be the space of all R
d-valued continuous paths (
ω
t
)
t
∈
R
+, with ω
0
=
0, equipped with the distance
(3)
ρ
(
ω
1
,
ω
2
)
∶
=
∑
i
=
1
∞
‍
2
-
i
[
(
max
⁡
t
∈
[
0
,
i
]
|
ω
t
1
-
ω
t
2
|
)
∧
1
]
.
For each fixed T
∈
[
0
,
∞
), we set Ω
T
:
=
{
ω
·
∧
T
:
ω
∈
Ω
}.From Peng [36], the distribution of the quadratic variation process of the G-Brownian motion has the following properties.Lemma 3.
For each fixeds
,
t
≥
0
,
〈
B
〉
s
+
t
-
〈
B
〉
s is identically distributed with 〈
B
〉
t and independent from Ω
s.Theorem 4.
〈
B
〉
t is N
(
[
σ
_
2
t
,
σ
¯
2
t
]
×
{
0
}
)-distributed; that is, for each φ
∈
C
l
,
L
i
p
(
R
),
(4)
E
^
[
φ
(
〈
B
〉
t
)
]
=
sup
⁡
σ
_
2
≤
v
≤
σ
¯
2
φ
(
v
t
)
.The sublinear expectationE
^ in the above theorem is called G-expectation defined on L
i
p
(
Ω
) via the following procedure. Let
(5)
L
i
p
(
Ω
T
)
∶
=
{
ϕ
(
B
t
1
∧
T
,
…
,
B
t
n
∧
T
)
:
n
∈
N
,
k
t
1
,
…
,
t
n
∈
[
0
,
∞
)
,
ϕ
∈
C
l
,
L
i
p
(
R
d
×
n
)
}
with B
t
=
ω
t
,
t
∈
[
0
,
∞
) for ω
∈
Ω, and
(6)
L
i
p
(
Ω
)
∶
=
⋃
n
=
1
∞
‍
L
i
p
(
Ω
n
)
.Firstly, construct a sequence ofd-dimensional random vectors (
ξ
)
i
=
1
∞ on a sublinear expectation space (
Ω
~
,
H
~
,
E
~
) such that ξ
i is G-normal distributed and ξ
i
+
1 is independent from (
ξ
1
,
…
,
ξ
i
) for each i
=
1,2
,
…. For each X
∈
L
i
p
(
Ω
) with
(7)
X
=
ϕ
(
B
t
1
-
B
t
0
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
,
some ϕ
∈
C
l
,
L
i
p
(
R
d
×
n
) and 0
=
t
0
<
t
1
<
⋯
<
t
n
<
∞, set
(8)
E
^
[
ϕ
(
B
t
1
-
B
t
0
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
]
∶
=
E
~
[
ϕ
(
t
1
-
t
0
ξ
1
,
…
,
t
n
-
t
n
-
1
ξ
n
)
]
.And the related conditional expectation of(9)
X
=
ϕ
(
B
t
1
,
B
t
2
-
B
t
1
,
…
,
B
t
n
-
B
t
n
-
1
)
under Ω
t
j is defined by
(10)
E
^
[
X
∣
Ω
t
j
]
∶
=
ψ
(
B
t
1
,
…
,
B
t
j
-
B
t
j
-
1
)
,
where
(11)
ψ
(
x
1
,
…
,
x
j
)
=
E
~
[
ϕ
(
x
1
,
…
,
x
j
,
t
j
+
1
-
t
j
ξ
j
+
1
,
…
,
t
n
-
t
n
-
1
ξ
n
)
]
.
E
^
[
·
] consistently defines a sublinear expectation on L
i
p
(
Ω
) and (
B
t
)
t
≥
0 is a G-Brownian motion.We letL
G
p
(
Ω
t
k
)
,
p
≥
1, denote the completion of
(12)
L
i
p
(
Ω
t
k
)
∶
=
{
φ
(
B
t
1
∧
T
,
…
,
B
t
n
∧
T
)
:
n
∈
N
,
k
t
1
,
…
,
t
n
∈
[
0
,
∞
)
,
φ
∈
C
l
,
L
i
p
(
R
d
×
n
)
}
under the norm ∥
X
∥
p
∶
=
(
E
^
[
|
X
|
p
]
)
1
/
p.


## 2.2. The PDEs

We assume the security market is defined on a sublinear expectation space(
Ω
,
H
,
E
^
) which opens continuously and offers a constant riskless interest rate r to all borrowers and lenders. And we also assume there are no transaction costs and/or taxes. Suppose that the price process (
S
t
)
t
≥
0 of some risky asset is given by
(13)
d
S
ν
=
μ
S
ν
d
ν
+
σ
~
S
ν
d
B
ν
,
S
0
=
S
,
where B is the G-Brownian motion with σ
_
2
=
-
E
^
[
-
B
1
2
] and σ
¯
2
=
E
^
[
B
1
2
], and μ and σ
~ are constants. Hence (
S
ν
)
ν
≥
0 is a geometric G-Brownian motion:
(14)
S
ν
=
S
t
exp
⁡
(
1
2
μ
(
ν
-
t
)
k
k
k
k
k
i
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
)
.Define(15)
I
ν
=
∫
0
ν
‍
S
τ
d
τ
.
Then I
ν satisfies the following stochastic differential equation:
(16)
d
I
ν
=
S
ν
d
ν
,
I
0
=
0
.LetX
ν
=
(
S
ν
,
I
ν
)
T, so
(17)
d
X
ν
=
μ
X
X
ν
d
ν
+
σ
~
X
X
ν
d
B
ν
,
X
0
=
x
=
(
S
,
0
)
,
where μ
X and σ
~
X are 2
×
2 matrixes:
(18)
μ
X
=
(
μ
0
1
0
)
,
σ
~
X
=
(
σ
~
0
0
0
)
.The typical Asian options in this security market are the following four types:the arithmetic average fixed strike call option with final payoff,(19)
max
⁡
(
I
T
T
-
K
,
0
)
;
the corresponding arithmetic average fixed strike put option with payoff,(20)
max
⁡
(
K
-
I
T
T
,
0
)
;
the arithmetic average floating strike call option with final payoff,(21)
max
⁡
(
S
T
-
I
T
T
,
0
)
;
and the corresponding arithmetic average floating strike put option with payoff,(22)
max
⁡
(
I
T
T
-
S
T
,
0
)
.LetV
(
t
,
S
t
,
I
t
) denote the price of one of the options at time t. Take the arithmetic average fixed strike call option for example, we know V
(
T
,
S
T
,
I
T
)
=
max
⁡
(
I
T
/
T
-
K
,
0
).Here we still assume the market operates in a risk neutral world [7].Peng [35] introduces the axiomatic conditions satisfied by the valuation mechanism. In the sublinear expectation space (
Ω
,
H
,
E
^
), E
^
t
[
·
] is a filtration consistent valuation operator. By this valuation mechanism, we have
(23)
V
(
t
,
S
t
,
I
t
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
∣
Ω
t
]
.
Here we assume the riskless interest rate is r.Since we assume the market is risk neutral, we haveμ
=
r. It is easy to verify that the discounted prices of the stock satisfy the valuation mechanism {
E
^
t
[
·
]
}
t
≥
0; that is,
(24)
S
t
=
E
^
t
[
e
-
r
(
ν
-
t
)
S
ν
]
,
0
≤
t
≤
ν
≤
T
.Let us assume(25)
d
X
ν
t
,
x
=
μ
X
X
ν
t
,
x
d
ν
+
σ
~
X
X
ν
t
,
x
d
B
ν
,
X
t
=
x
=
(
S
,
I
)
,
with
(26)
μ
X
=
(
r
0
1
0
)
,
σ
~
X
=
(
σ
~
0
0
0
)
.The following result is a corollary of the nonlinear Feynman-Kac formula of Peng [36].Corollary 5.
Let(27)
Y
s
t
,
x
=
E
^
[
∫
s
T
‍
e
-
r
(
T
-
s
)
Φ
(
X
T
t
,
x
)
k
k
+
∫
s
T
‍
e
-
r
(
ν
-
s
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
k
k
+
∫
s
T
‍
e
-
r
(
ν
-
s
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
∣
Ω
s
]
,
where Φ
:
R
2
→
R is a given Lipschitz function and f
,
g
:
R
2
×
R
→
R are given Lipschitz functions; that is, |
ϕ
(
x
,
y
)
-
ϕ
(
x
′
,
y
′
)
|
≤
K
(
|
x
-
x
′
|
+
|
y
-
y
′
|
), for each x
,
x
′
∈
R
2,  y
,
y
′
∈
R,  ϕ
=
f and g. And denote
(28)
u
(
t
,
x
)
=
Y
t
t
,
x
.
Then
(29)
u
(
t
,
x
)
=
E
^
[
∫
t
t
+
δ
‍
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
]
is a viscosity solution of the following PDE:
(30)
∂
t
u
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
u
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
+
g
(
x
,
u
)
)
σ
2
+
〈
D
u
(
t
,
x
)
,
μ
X
x
〉
-
r
u
+
f
(
x
,
u
)
=
0
,
u
(
T
,
x
)
=
Φ
(
x
)
.Proof.
The proof of expression (29) is the same as that of Peng [36]. By the proof of Peng's theorem of the nonlinear Feynman-Kac formula and the boundedness of e
-
r
(
T
-
t
), we can easily prove that u is a Lipschitz function; that is,
(31)
|
u
(
t
,
x
)
-
u
(
t
,
x
′
)
|
≤
C
|
x
-
x
′
|
,
(32)
E
^
[
|
Y
s
t
,
x
|
]
≤
C
(
1
+
|
x
|
)
.
Let us see the continuity of u with respect to t. Firstly, we have
(33)
|
u
(
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
≤
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
+
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kkk
k
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
|
]
.
Replace e
-
r
δ with its Taylor's expansion in the first term on the right side of (33), we have
(34)
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
=
E
^
[
|
(
1
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
k
k
k
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
=
E
^
[
|
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
k
k
k
+
(
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
|
]
≤
E
^
[
|
X
t
+
δ
t
,
x
-
x
|
]
+
E
^
[
|
(
-
r
δ
+
(
r
δ
)
2
+
⋯
)
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
|
]
.
By the following formulas [36]
(35)
E
^
[
|
X
t
+
δ
t
,
x
-
x
|
2
]
≤
C
(
1
+
x
2
)
δ
,
E
^
[
|
X
t
+
δ
t
,
x
|
2
]
≤
C
(
1
+
x
2
)
,
and (32), we have
(36)
E
^
[
|
e
-
r
δ
u
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
u
(
t
+
δ
,
x
)
|
]
≤
C
′
(
1
+
|
x
|
2
)
1
/
2
(
δ
1
/
2
+
δ
+
δ
2
2
⋯
)
.
For the second term on the right side of (33), we have
(37)
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
|
]
=
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
kkkik
k
k
k
k
×
(
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
-
f
(
x
,
Y
t
t
,
x
)
)
d
ν
kkkik
k
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
kkkik
k
k
k
k
k
i
×
(
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
-
g
(
x
,
Y
t
t
,
x
)
)
d
〈
B
〉
ν
∫
t
t
+
δ
‍
|
]
+
E
^
[
|
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
x
,
Y
t
t
,
x
)
d
ν
kkkik
ki
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
x
,
Y
t
t
,
x
)
d
〈
B
〉
ν
|
]
≤
C
(
1
+
|
x
|
)
δ
.
Therefore u is continuous in t.
Now we are going to prove thatu
(
t
,
x
) is a viscosity solution of PDE (30). For fixed (
t
,
x
)
∈
(
0
,
T
)
×
R
2, let ψ
∈
C
b
2,3
(
(
0
,
T
)
×
R
2
) such that ψ
≥
u and ψ
(
t
,
x
)
=
u
(
t
,
x
). By (29) and Taylor's expansion, for δ
∈
(
0
,
T
-
t
),
(38)
0
≤
E
^
[
e
-
r
δ
ψ
(
t
+
δ
,
X
t
+
δ
t
,
x
)
-
ψ
(
t
,
x
)
∫
t
t
+
δ
‍
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
f
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
ν
kk
+
∫
t
t
+
δ
‍
e
-
r
(
ν
-
t
)
g
(
X
ν
t
,
x
,
Y
ν
t
,
x
)
d
〈
B
〉
ν
]
≤
E
^
[
1
2
-
r
ψ
(
t
,
x
)
δ
+
∂
t
ψ
(
t
,
x
)
δ
kkk
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
δ
kkk
+
〈
1
2
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
(
〈
B
〉
t
+
δ
-
〈
B
〉
t
)
kkk
1
2
+
f
(
x
,
ψ
)
δ
+
g
(
x
,
ψ
)
(
〈
B
〉
t
+
δ
-
〈
B
〉
t
)
]
+
C
(
1
+
|
x
|
+
|
x
|
2
+
|
x
|
3
)
δ
3
/
2
≤
(
∂
t
ψ
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
kkkkkkkkkkkkkk
+
g
(
x
,
ψ
)
1
2
)
σ
2
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
k
i
-
r
ψ
(
t
,
x
)
+
f
(
x
,
ψ
)
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
)
δ
+
C
(
1
+
|
x
|
+
|
x
|
2
+
|
x
|
3
)
δ
3
/
2
.
It is easy to check that
(39)
∂
t
ψ
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
(
1
2
〈
D
2
ψ
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
+
g
(
x
,
ψ
)
)
σ
2
+
〈
D
ψ
(
t
,
x
)
,
μ
X
x
〉
-
r
ψ
(
t
,
x
)
+
f
(
x
,
ψ
)
≥
0
.
Thus u is a viscosity subsolution of (30). Similarly we can prove that u is a viscosity supersolution of (30).
For the arithmetic average Asian options, we have(40)
V
(
t
,
x
)
=
V
(
t
,
S
,
I
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
]
.
Then by Corollary 5, V
(
t
,
x
) satisfies the following PDE:
(41)
∂
t
V
(
t
,
x
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
〈
D
2
V
(
t
,
x
)
σ
~
X
x
,
σ
~
X
x
〉
σ
2
+
〈
D
V
(
t
,
x
)
,
μ
X
x
〉
-
r
V
(
t
,
x
)
=
0
,
V
(
T
,
x
)
=
φ
(
T
,
x
)
.
By (18), the above PDE is
(42)
∂
t
V
(
t
,
S
,
I
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
S
2
σ
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
φ
(
T
,
S
,
I
)
,
where
(43)
φ
(
T
,
S
,
I
)
=
{
max
⁡
(
I
T
-
K
,
0
)
,
for
the
arithmetic
average
fixed
strike
call
option
,
max
⁡
(
K
-
I
T
,
0
)
,
for
the
arithmetic
average
fixed
strike
put
option
,
max
⁡
(
S
-
I
T
,
0
)
,
for
the
arithmetic
average
floating
strike
call
option
,
max
⁡
(
I
T
-
S
,
0
)
,
for
the
arithmetic
average
floating
strike
put
option
.


## 2.3. The Arithmetic Average Fixed Strike Asian Call Option

Let us see the arithmetic average fixed strike Asian call option. From (14), we have
(44)
S
ν
t
,
S
=
S
exp
⁡
⁡
{
μ
(
ν
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
}
.
Then by the definition of I
ν, we have
(45)
I
ν
t
,
I
=
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
k
k
k
k
k
k
k
k
k
k
+
σ
~
(
B
τ
-
B
t
)
{
μ
(
τ
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
}
d
τ
.
Since
(46)
V
(
t
,
x
)
=
V
(
t
,
S
,
I
)
=
E
^
[
e
-
r
(
T
-
t
)
V
(
T
,
S
T
,
I
T
)
]
=
E
^
[
e
-
r
(
T
-
t
)
(
I
T
T
-
K
)
+
]
,
V
(
t
,
S
,
I
)
=
E
^
[
×
(
T
)
-
1
(
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
)
-
K
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
)
+
e
-
r
(
T
-
t
)
i
i
i
k
×
(
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
1
2
kkkkkkkkkkkkkkki
i
i
i
k
i
-
1
2
σ
~
2
(
〈
B
〉
τ
-
〈
B
〉
t
)
kkkkkkkkkkkkkkki
i
i
i
i
k
1
2
+
σ
~
(
B
τ
-
B
t
)
}
d
τ
S
∫
t
ν
‍
)
i
i
k
i
×
T
-
1
(
I
+
S
∫
t
ν
‍
exp
⁡
⁡
{
μ
(
τ
-
t
)
)
-
K
(
(
(
I
+
S
∫
t
ν
‍
exp
⁡
{
μ
(
τ
-
t
)
)
+
]
.
By the convexities of function f
(
x
)
=
(
x
-
k
)
+ and the sublinear expectation, we can say
(47)
V
(
t
,
λ
S
1
+
(
1
-
λ
)
S
2
,
I
)
≤
λ
V
(
t
,
S
1
,
I
)
+
(
1
-
λ
)
V
(
t
,
S
2
,
I
)
;
that is, V
(
t
,
S
,
I
) is convex about S. So in the numerical computation, the second order difference of V is positive. Then the numerical solutions of (43) are equivalent to those of PDE:
(48)
∂
t
V
(
t
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
(
I
T
-
K
)
+
.Theoretically, the above PDE is defined on domainD
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
∞
)
×
[
0
,
∞
)
}. But generally, closed form solutions cannot be found; we have to find numerical solutions. We have no interest in finding the solution in the entire infinite domain. Therefore we consider the domain D
*
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
S
*
)
×
[
0
,
I
*
)
}.Now we determine the boundary conditions of (48). If I
≥
K
T, that is, I
/
T
≥
K, for the positiveness of S
ν, we have I
T
/
T
≥
K. Therefore the final payoff on the option is certain to be positive. At time T, this payoff can be expressed as
(49)
I
T
T
-
K
=
I
T
-
K
+
1
T
∫
t
T
‍
S
ν
d
ν
.
This payoff can also be obtained by using the following self-financing duplicating portfolio strategy. Assume that an investor commits (
I
/
T
-
K
)
e
-
r
(
T
-
t
) into riskless bonds in order to secure the return promised by the first part of the right of (49), that is, (
I
/
T
-
K
) at time T. If the investor is to be certain of accruing the return promised by the second part of the right of (49), he is obliged to transfer a certain portion of stock (equal to (
1
/
T
)
e
-
r
(
T
-
ν
)
Δ
ν) to a riskless bond in every time interval (
ν
,
ν
+
Δ
ν
). Overall this strategy requires a sum equal to (
I
/
T
-
K
)
e
-
r
(
T
-
t
) plus the portion of a stock which can be expressed as follows:
(50)
∫
t
T
‍
1
T
e
-
r
(
T
-
ν
)
d
ν
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
.
Then the price of the option when I
t
≥
K
T must be equal to
(51)
V
(
t
,
S
,
I
)
=
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.The above discussion is similar to that in the model without volatility uncertainty. Equation (51) obviously satisfies (48). So we only need to find the solution of (48) in the domain D
*
=
{
(
t
,
S
,
I
)
∈
[
0
,
T
]
×
[
0
,
S
*
)
×
[
0
,
K
T
)
}.On the boundaryI
=
K
T, the boundary condition from (51) can be written as
(52)
V
(
t
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.We note that ifS
t
=
0, (13) tells us that S
(
ν
) is also equal to zero for ν
∈
[
t
,
T
] and I
(
T
) is thus equal to I
(
t
). By (46), we can write the following boundary condition:
(53)
V
(
t
,
0
,
I
)
=
e
-
r
(
T
-
t
)
(
I
T
-
K
)
+
.We need another boundary condition atS
=
S
max
⁡
=
S
*. Hugger [37] gives the boundary condition at S
=
S
max
⁡ in the model without volatility uncertainty. By the same discussion, we can get the same boundary condition in uncertain volatility model:
(54)
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
S
max
⁡
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
0
}
.So, the values of the arithmetic average fixed strike Asian call option satisfy(55)
∂
t
V
(
t
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
I
)
+
r
S
V
S
(
t
,
S
,
I
)
+
S
V
I
(
t
,
S
,
I
)
-
r
V
(
t
,
S
,
I
)
=
0
,
V
(
T
,
S
,
I
)
=
(
I
T
-
K
)
+
,
V
(
t
,
0
,
I
)
=
e
-
r
(
T
-
t
)
(
I
T
-
K
)
+
,
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
(
T
-
t
)
+
S
max
⁡
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
0
}
,
V
(
t
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
.For the arithmetic average fixed strike Asian put option with terminal payoffp
(
T
,
S
T
,
I
T
)
=
max
⁡
⁡
(
K
-
I
T
/
T
,
0
), we will value this option by the following put-call parity.It is easily to see that we have(56)
max
⁡
(
I
T
T
-
K
,
0
)
=
max
⁡
(
K
-
I
T
T
,
0
)
+
1
T
∫
0
T
‍
S
ν
d
ν
-
K
.
So
(57)
c
(
t
,
S
t
,
I
t
)
=
e
-
r
(
T
-
t
)
E
^
t
[
max
⁡
(
I
T
T
-
K
,
0
)
]
=
e
-
r
(
T
-
t
)
E
^
t
[
max
⁡
(
K
-
I
T
T
,
0
)
+
1
T
∫
0
T
‍
S
ν
d
ν
]
-
e
-
r
(
T
-
t
)
K
=
p
(
t
,
S
t
,
I
t
)
-
e
-
r
(
T
-
t
)
K
+
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
ν
d
ν
+
e
-
r
(
T
-
t
)
E
t
σ
¯
[
1
T
∫
t
T
‍
S
ν
d
ν
]
.e
-
r
(
T
-
t
)
E
t
σ
¯
[
(
1
/
T
)
∫
t
T
‍
S
ν
d
ν
] is the time t value of time T claim (
1
/
T
)
∫
t
T
‍
S
ν
d
ν, where the stock has certain volatility σ
~
σ
¯. To get claim (
1
/
T
)
∫
t
T
‍
S
ν
d
ν at T, let us see the replicating strategy: in any time interval (
ν
,
ν
+
Δ
ν
), transfer a certain portion (
1
/
T
)
e
-
r
(
T
-
ν
)
Δ
ν of stock S
ν to a riskless bond. This strategy ensures the final payoff of (
1
/
T
)
∫
t
T
‍
S
ν
d
ν. The stocks that should be transferred are
(58)
∫
t
T
‍
1
T
e
-
r
(
T
-
ν
)
d
ν
=
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
.
Then the time t value of (
1
/
T
)
∫
t
T
‍
S
ν
d
ν equals (
1
/
r
T
)
(
1
-
e
-
r
(
T
-
t
)
)
S
t. So the put-call parity can be written as
(59)
c
(
t
,
S
,
I
)
+
K
e
-
r
(
T
-
t
)
=
p
(
t
,
S
,
I
)
+
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
ν
d
ν
+
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
S
t
.


## 2.4. The Arithmetic Average Floating Strike Asian Call Option

For the arithmetic average floating strike Asian call (put) option, the PDE (43) still applies with terminal condition
(60)
V
(
T
,
S
,
I
)
=
(
S
-
I
T
)
+
(
V
(
T
,
S
,
I
)
=
(
I
T
-
S
)
+
)
.In Ingersoll [13] where there is no volatility uncertainty, due to the linear homogeneity of the PDE and its terminal condition, their PDE can be reduced to a one-dimensional PDE. In our uncertain volatility model, the PDE (43) is not linear homogeneous anymore. But we can still reduce PDE (43) to a one-dimensional PDE by the following argument.For the arithmetic average floating strike Asian call option, we know(61)
V
(
t
,
S
t
,
I
t
)
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
T
-
I
T
T
)
+
]
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
T
-
I
t
+
∫
t
T
‍
S
ν
d
ν
T
)
+
]
.Denote(62)
M
t
ν
=
exp
⁡
(
μ
(
ν
-
t
)
-
1
2
σ
~
2
(
〈
B
〉
ν
-
〈
B
〉
t
)
+
σ
~
(
B
ν
-
B
t
)
)
.
Then by (14), we have S
ν
=
S
t
M
t
ν and S
t
=
S
ν
M
ν
t. So
(63)
I
t
=
∫
0
t
‍
S
ν
d
ν
=
∫
0
t
‍
S
t
M
t
ν
d
ν
=
S
t
∫
0
t
‍
M
t
ν
d
ν
,
V
(
t
,
S
t
,
I
t
)
=
E
^
t
[
e
-
r
(
T
-
t
)
(
S
t
M
t
T
-
S
t
∫
0
t
‍
M
t
ν
d
ν
+
S
t
∫
t
T
‍
M
t
ν
d
ν
T
)
+
]
=
S
t
E
^
t
[
e
-
r
(
T
-
t
)
(
M
t
T
-
∫
0
t
‍
M
t
ν
d
ν
+
∫
t
T
‍
M
t
ν
d
ν
T
)
+
]
∶
=
S
t
H
(
t
,
R
t
)
with R
t
=
I
t
/
S
t.Then the functionH
(
t
,
R
) satisfies the following PDE:
(64)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
which is a one-dimensional PDE.For the European style Asian option, we must impose boundary conditions both atR
=
0 and as R
→
∞. The boundary condition as R
→
∞ is simple. The only way that R can tend to infinity is that S tends to zero. In this case the option will not be exercised, and so
(65)
lim
⁡
R
→
∞
H
(
t
,
R
)
=
0
.
For computational purpose, we truncate the asset region (
0
,
+
∞
) into (
0
,
R
max
⁡
), where R
max
⁡ is sufficiently large to ensure the accuracy of the solution. If R
=
0, the PDE degenerates into the first-order hyperbolic PDE
(66)
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
.SoH
(
t
,
R
) satisfies
(67)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
H
(
T
,
R
)
=
(
1
-
R
T
)
+
,
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
H
(
t
,
R
max
⁡
)
=
0
.By the same argument, for the arithmetic average floating strike Asian put option, its prices could be calculated by(68)
V
(
t
,
S
,
I
)
=
S
H
(
t
,
R
)
,
with H
(
t
,
R
) satisfying the following PDE:
(69)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
H
(
T
,
R
)
=
(
R
T
-
1
)
+
,
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
H
(
t
,
R
max
⁡
)
=
(
R
max
⁡
T
-
exp
⁡
{
-
r
(
T
-
t
)
}
)
.Generally, the put-call parity will not hold any more in the uncertain volatility model. Let us see the following argument. The mature payoff of a portfolio with one European arithmetic average floating strike Asian call holding long and one put holding short is(70)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
-
S
T
max
⁡
(
R
T
T
-
1,0
)
=
S
T
-
R
T
S
T
T
.
The value of this portfolio is equal to one consisting of one stock and a financial product whose payoff is -
R
T
S
T
/
T. The values of the financial product with final payoff -
R
T
S
T
/
T satisfy the following PDE:
(71)
∂
t
H
(
t
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
t
,
R
)
+
(
1
-
r
R
)
H
R
(
t
,
R
)
=
0
,
with terminal condition H
(
T
,
R
T
)
=
-
R
T
S
T
/
T. We seek a solution of (71) of the form
(72)
H
(
t
,
R
)
=
a
t
+
b
t
R
with a
T
=
0 and b
T
=
-
1
/
T. Substituting (72) into (71) and satisfying the boundary conditions, we find that
(73)
a
t
=
-
1
r
T
(
1
-
e
-
r
(
T
-
t
)
)
,
b
t
=
-
1
T
e
-
r
(
T
-
t
)
.
For
(74)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
=
S
T
max
⁡
(
R
T
T
-
1,0
)
+
S
T
-
R
T
S
T
T
,
by the valuation mechanism in space (
Ω
,
H
,
E
^
), we take the discounted G-conditional expectation on both sides and get
(75)
E
^
t
[
e
-
r
(
T
-
t
)
S
T
max
⁡
(
1
-
R
T
T
,
0
)
]
=
e
-
r
(
T
-
t
)
E
^
t
[
S
T
max
⁡
(
R
T
T
-
1,0
)
+
S
T
-
R
T
S
T
T
]
≤
e
-
r
(
T
-
t
)
E
^
t
[
S
T
max
⁡
(
R
T
T
-
1,0
)
]
+
e
-
r
(
T
-
t
)
E
^
t
[
S
T
-
R
T
S
T
T
]
.
That is
(76)
c
(
t
,
S
,
I
)
≤
p
(
t
,
S
,
I
)
+
S
-
S
r
T
(
1
-
e
-
r
(
T
-
t
)
)
-
1
T
e
-
r
(
T
-
t
)
∫
0
t
‍
S
τ
d
τ
,
where c
(
t
,
S
,
I
) is the time t price of the arithmetic average floating strike Asian call option and p
(
t
,
S
,
I
) is the time t price of the arithmetic average floating strike Asian put option.


## 3. Numerical Computation

In the model without volatility uncertainty, Asian options are much more difficult to value than regular stock options. Standard techniques tend to be impractical, inaccurate, or slow. For example, traditional binomial lattice methods require such enormous amounts of computer memory for the necessity of keeping track of every possible path throughout the tree that they are effectively unusable. Monte Carlo simulation works well for European style options (see [12]) but is relatively slow. Sun [38] proposed to use the weighted upwind method [39] to price Asian options, but he did not give numerical results. The Asian option pricing in the uncertain volatility model will be even more difficult than in the model without volatility uncertainty.For the arithmetic average fixed strike Asian call options, due to the convexity of the valueV
(
t
,
S
,
I
) with respect to variable S, (43) can be written as (48) which are the equations satisfied by the arithmetic average fixed strike Asian call options with certain volatility σ
~
σ
¯.Zvan et al. [17] point out that the Crank-Nicolson scheme in conjunction with the van Leer flux limiter method can be used to numerically value Asian options with certain volatilities. The van Leer limiter has the property that it is total variation diminishing and thus produces oscillation free second order accurate solutions. But the van Leer limiter is nonlinear, so the solutions should be obtained by using full Newton iteration at each time level. In this paper, we will use the alternating direction implicit (ADI) methods (see Duffy [40]).The arithmetic average fixed strike Asian put options can be valued by the put-call parity.For the arithmetic average floating strike Asian call (put) options, we need to solve (67) (see (69)). These PDEs are one-dimensional in space. And they are not linear equations but HJB equations in some sense. For such PDEs, we will use the fully implicit positive coefficient finite different method [41] to solve them. We will see that this method can ensure the numerical solutions converge to the viscosity solutions, and is accurate, efficient, and quick to be implemented.All the calculations were performed by Matlab 2009, on a computer with 1.83 GHz hard disk and 512 MB memory.


### 3.1. Numerical Computation of Fixed Strike Asian Options

An equivalent formulation of (48) in terms of the average A
t
=
(
1
/
t
)
∫
0
t
‍
S
τ
d
τ is given in Barraquand and Pudet [16] in certain volatility model with volatility σ
=
σ
~
σ
¯:
(77)
∂
t
V
(
t
,
S
,
A
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
A
)
+
r
S
V
S
(
t
,
S
,
A
)
+
S
-
A
t
V
A
(
t
,
S
,
A
)
-
r
V
(
t
,
S
,
A
)
=
0
,
V
(
T
,
S
,
A
)
=
(
A
-
K
)
+
.Equation (77) is convection dominated in the A direction because there is no diffusion effect in this dimension. Barraquand and Pudet [16] find that an explicit centrally weighted scheme for (77) is unstable. In particular, the convective term in the A dimension becomes very large as t
→
0. Barraquand and Pudet also note that implicit centrally weighted schemes will generally produce unsatisfactory results because of the numerical diffusion introduced by this first order accurate in time scheme.Zvan et al. [17] point out that under the condition of convection dominated, solutions generated by using a centrally weighted scheme on the convective term for (77) cannot be ensured to be free of oscillations. And solutions generated by the first order upstream weighting scheme on the convective term are no longer oscillatory, but their profiles are too diffuse. To produce oscillation free solutions without the excessive diffusion of first order upstream weighting, Zvan et al. [17] suggest the nonlinear van Leer flux limiter in conjunction with the Crank-Nicolson scheme. Since the flux limiter is nonlinear, the solutions should be obtained by using full Newton iteration at each time level.The ADI method, pioneered in the United States by Douglas, Rachford, Peaceman, Gunn, and others, has a number of advantages. First, explicit difference methods are rarely used to solve initial boundary value problems because of their poor stability. Implicit methods have superior stability properties but unfortunately they are difficult to solve in two and more dimensions. Consequently, ADI methods become an alternative because they can be programmed by solving a simple trigonal system of equations. For the convection-dominated problems, we could use the exponentially fitted schemes in each leg of the approximate scheme to ensure the stability of the results [40]. Since Barraquand and Pudet [16] have used the forward shooting grid method and Zvan et al. [17] have used the van Leer flux limiter, we will test the ADI method in this paper by comparing our results with theirs.For convenience, let us first convert (55) to be the following forward equation in time by substituting t with τ
=
T
-
t:
(78)
-
∂
τ
V
(
τ
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
τ
,
S
,
I
)
+
r
S
V
S
(
τ
,
S
,
I
)
+
S
V
I
(
τ
,
S
,
I
)
-
r
V
(
τ
,
S
,
I
)
=
0
,
V
(
0
,
S
,
I
)
=
(
I
T
-
K
)
+
,
V
(
τ
,
0
,
I
)
=
e
-
r
τ
(
I
T
-
K
)
+
,
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
τ
+
S
max
⁡
r
T
(
1
-
e
-
r
τ
)
,
0
}
,
V
(
τ
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
τ
)
S
.We letτ
n,  n
=
0,1
,
…
,
N, be a set of portion points in [
0
,
T
] satisfying 0
=
τ
0
<
τ
1
<
⋯
<
τ
N
=
T and use even partition of time; that is, τ
n
=
n
Δ
τ,  n
=
0,1
,
…
,
N. Define the partitions of space variables as 0
=
S
0
<
S
1
<
⋯
<
S
M
=
S
max
⁡ and 0
=
I
0
<
I
1
<
⋯
<
I
J
=
I
max
⁡.With ADI, we march from time leveln to time level n
+
1
/
2 and then from time level n
+
1
/
2 to time level n
+
1. In this case we use exponential fitting in all space variables and implicit Euler in time. The first leg is given by the scheme
(79)
-
V
i
,
j
n
+
1
/
2
-
V
i
,
j
n
(
1
/
2
)
Δ
τ
+
ρ
i
,
j
[
1
Δ
S
i
Δ
S
i
+
1
V
i
-
1
,
j
n
+
1
/
2
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
-
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
+
1
Δ
S
i
+
1
2
V
i
+
1
,
j
n
+
1
/
2
]
+
r
S
i
V
i
+
1
,
j
n
+
1
/
2
-
V
i
,
j
n
+
1
/
2
Δ
S
i
+
1
+
S
i
V
i
,
j
+
1
n
-
V
i
,
j
-
1
n
Δ
I
j
+
1
+
Δ
I
j
-
r
V
i
,
j
n
+
1
/
2
=
0
,
with
(80)
ρ
i
,
j
=
r
S
i
Δ
S
i
+
1
2
coth
⁡
(
r
S
i
Δ
S
i
+
1
σ
¯
2
σ
~
2
S
i
2
)
,
coth
⁡
(
x
)
=
e
2
x
+
1
e
2
x
-
1
.
Let
(81)
α
i
=
-
ρ
i
,
j
Δ
S
i
Δ
S
i
+
1
,
β
i
=
-
ρ
i
,
j
(
Δ
S
i
+
1
)
2
-
r
S
i
Δ
S
i
+
1
,
γ
i
=
-
α
i
-
β
i
+
r
.
Then the scheme (79) can be written as
(82)
α
i
V
i
-
1
,
j
n
+
1
/
2
+
(
γ
i
+
2
Δ
τ
)
V
i
,
j
n
+
1
/
2
+
β
i
V
i
+
1
,
j
n
+
1
/
2
=
2
Δ
τ
V
i
,
j
n
+
S
i
Δ
I
j
+
1
+
Δ
I
j
V
i
,
j
+
1
n
-
S
i
Δ
I
j
+
1
+
Δ
I
j
V
i
,
j
-
1
n
.
We let
(83)
V
j
n
=
[
V
1
,
j
n
,
V
2
,
j
n
,
…
,
V
M
-
1
,
j
n
]
T
,
V
j
n
+
1
/
2
=
[
V
1
,
j
n
+
1
/
2
,
V
2
,
j
n
+
1
/
2
,
…
,
V
M
-
1
,
j
n
+
1
/
2
]
T
,
f
j
n
+
1
/
2
=
[
α
1
V
0
,
j
n
+
1
/
2
,
0
,
…
,
0
,
β
M
-
1
V
M
,
j
n
+
1
/
2
]
T
,
r
j
n
=
[
(
2
Δ
τ
V
1
,
j
n
+
S
1
Δ
I
j
+
1
+
Δ
I
j
V
1
,
j
+
1
n
k
k
-
S
1
Δ
I
j
+
1
+
Δ
I
j
V
1
,
j
-
1
n
)
,
…
,
k
(
2
Δ
τ
V
M
-
1
,
j
n
+
S
M
-
1
Δ
I
j
+
1
+
Δ
I
j
V
M
-
1
,
j
+
1
n
k
.
k
-
S
M
-
1
Δ
I
j
+
1
+
Δ
I
j
V
M
-
1
,
j
-
1
n
)
]
T
,
A
1
=
[
γ
1
β
1
α
2
γ
2
β
2
⋱
⋱
⋱
α
M
-
1
γ
M
-
1
]
(
M
-
1
)
×
(
M
-
1
)
.
Then we can write (82) as the following equivalent matrix form
(84)
A
1
V
j
n
+
1
/
2
+
f
j
n
+
1
/
2
=
r
j
n
.The second leg of the discretized PDE in (78) is given by the scheme
(85)
-
V
i
,
j
n
+
1
-
V
i
,
j
n
+
1
/
2
(
1
/
2
)
Δ
τ
+
ρ
i
,
j
[
1
Δ
S
i
Δ
S
i
+
1
V
i
-
1
,
j
n
+
1
/
2
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
kk
-
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
kk
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
+
1
Δ
S
i
+
1
2
V
i
+
1
,
j
n
+
1
/
2
]
+
r
S
i
V
i
+
1
,
j
n
+
1
/
2
-
V
i
,
j
n
+
1
/
2
Δ
S
i
+
1
+
S
i
V
i
,
j
+
1
n
+
1
-
V
i
,
j
n
+
1
Δ
I
j
+
1
-
r
V
i
,
j
n
+
1
/
2
=
0
.
Let
(86)
δ
j
=
2
Δ
τ
+
S
i
Δ
I
j
+
1
,
ϵ
j
=
-
S
i
Δ
I
j
+
1
.
Then the scheme (85) can be written as
(87)
δ
j
V
i
,
j
n
+
1
-
ϵ
j
V
i
,
j
+
1
n
+
1
=
-
α
i
V
i
-
1
,
j
n
+
1
/
2
-
(
γ
i
-
2
Δ
τ
)
V
i
,
j
n
+
1
/
2
-
β
i
V
i
+
1
,
j
n
+
1
/
2
.We let(88)
V
i
n
+
1
=
[
V
i
,
0
n
+
1
,
V
i
,
1
n
+
1
,
…
,
V
i
,
J
-
1
n
+
1
]
T
,
V
i
n
+
1
/
2
=
[
V
i
,
0
n
+
1
/
2
,
V
i
,
1
n
+
1
/
2
,
…
,
V
i
,
J
-
1
n
+
1
/
2
]
T
,
F
i
=
[
0,0
,
…
,
0
,
ϵ
J
-
1
V
i
,
J
n
+
1
]
J
T
,
R
i
n
+
1
/
2
=
[
(
-
α
i
V
i
-
1,0
n
+
1
/
2

-

(
γ
i
-
2
Δ
τ
)
V
i
,
0
n
+
1
/
2

-

β
i
V
i
+
1,0
n
+
1
/
2
)
,
…
,
(
-
α
i
V
i
-
1
,
J
-
1
n
+
1
/
2
-
(
γ
i
-
2
Δ
τ
)
V
i
,
J
-
1
n
+
1
/
2
-
β
i
V
i
+
1
,
J
-
1
n
+
1
/
2
)
]
T
,
A
2
=
[
δ
0
-
ϵ
0
0
δ
1
-
ϵ
1
⋱
⋱
-
ϵ
J
-
2
δ
J
-
1
]
J
×
J
.
The matrix form of (87) is
(89)
A
2
V
i
n
+
1
+
F
i
n
+
1
=
R
i
n
+
1
/
2
.
The matrix equations (84) and (89) can be solved by using LU decomposition.The goal of the computation is to find the fair price of the Asian option at emission, that is,V
(
S
0
,
I
0
,
0
). We need to make sense of any value I
∈
[
0
,
I
max
⁡
] and stock price S
∈
[
0
,
S
max
⁡
] at time t
=
0. There are no problems with the stock prices, but the integral at time zero can only be 0 for continuous averages; that is, I
0
=
0.Figure1 is the results we obtained by using the above discretization scheme with σ
~
=
0.8,  σ
¯
=
0.5,  r
=
0.10,  T
=
1, and K
=
100. The point marked by * is the price V
(
0
,
S
,
0
) of the Asian option. We choose S
max
⁡
=
3
K to ensure the desirable accuracy.Figure 1
Fixed strike Asian call value whenσ
~
=
0.8,  σ
¯
=
0.5,  r
=
0.10,  T
=
1, and K
=
100.Table1 shows the convergence of the results with the refinement of the grid. From (78), the values of variable I lie in the interval [
0
,
K
T
], where K
T is the boundary and 0 is the point where the continuous averages Asian options are valued, so there is no point that is less important on the interval [
0
,
K
T
]. Hence we use uniform grids in direction I.Table 1
Successive grid refinements demonstrating convergence of the results with the refinements whenr
=
0.10, T
=
1, σ
~
=
0.8, σ
-
=
0.5, and K
=
95. The values are the option prices at S
=
100. Δ
S and Δ
I denote the grid spacing in directions S and I, respectively. Δ
t
=
0.04. Execution times (in seconds) are in parentheses.
Δ
I
=
0.95
Δ
I
=
0.475
Δ
I
=
0.2375
Δ
I
=
0.11875
Δ
I
=
0.059375
Δ
S
=
5.70
14.532
14.19
14.021
13.94
13.89
(1.01)
(2.10)
(7.081)
(27.225)
(90)
Δ
S
=
2.85
14.42
14.08
13.914
13.82
13.781
(3.08)
(7.13)
(23.012)
(76.51)
(292)
Δ
S
=
1.425
14.39
14.053
13.879
13.790
13.74
(10.22)
(24.71)
(76.01)
(238)
(940)In our uncertain volatility model, the pricing PDEs (48) are similar to the PDEs in certain volatility model with volatility σ
=
σ
~
σ
¯. Barraquand and Pudet [16] and Zvan et al. [17] give numerical results of Asian options with certain volatilities. In the uncertain volatility model, when σ
~
=
0.8,  σ
¯
=
0.5, the Asian option PDE is the same as that of the Asian option with certain volatility σ
=
0.4. When the parameters r
=
0.10,  T
=
1,  K
=
95, and grid Δ
I
=
0.11875,  Δ
S
=
2.85, the price of the Asian call option we calculate is 13.82. The results with the same parameters of Barraquand and Pudet [16] and Zvan et al. [17] are 13.825 and 13.721, respectively.Table2 is the comparison of our results with those of Zvan et al. [17]. Z, F, and V in Table 2 refer to the result of Zvan et al. [17]. Our results are denoted by Asian fixed strike, which calculated with Δ
τ
=
0.01,0.02,0.04,  Δ
I
=
0.059375,0.11875,0.2375, and Δ
S
=
2.85 for maturities of one quarter, half a year, and one year, respectively. Since I
max
⁡
=
K
T, the partitions in I direction are related with the maturity T. The execution time is much less than theirs. The grid spacing was chosen to achieve an accuracy of at least 0.10% of S. Our results are close to that of Zvan et al. [17] and the difference is less than 0.10% of S.Table 2
Fixed strike Asian call values whenr
=
0.10 and K
=
95. The values are the option prices at S
=
100. Our results are denoted by Asian fixed strike, which calculated with Δ
τ
=
0.01,0.02,0.04, Δ
I
=
0.059375,0.11875,0.2375, and Δ
S
=
2.85 for maturities of one quarter, half a year, and one year, respectively. Z, F, and V refer to the result of Zvan et al. [17]. Execution times (in seconds) are in parentheses.
σ
~
σ
-
T
Asian fixed strike
Z, F, and V
0.2
0.5
0.25
6.190 (23.0)
6.133
0.5
7.199 (24.0)
7.244
1
9.204 (23.6)
9.316
1
0.25
6.56 (23.0)
6.501
0.5
7.99 (24.1)
7.921
1
10.44 (25.0)
10.309
2
0.25
8.30 (24.1)
8.123
0.5
10.55 (24.4)
10.357
1
13.91 (25.2)
13.721Our calculation procedure is different from Zvan et al. [17], as well as Barraquand and Pudet [16], in the following four aspects. First, our pricing PDE (78) is in terms of the running sum I
t. Zvan et al. [17] and Barraquand and Pudet [16] use the pricing PDE
(90)
-
∂
τ
V
(
τ
,
S
,
A
)
+
1
2
S
2
σ
2
V
S
S
(
τ
,
S
,
A
)
+
r
S
V
S
(
τ
,
S
,
A
)
+
S
-
A
t
V
A
(
τ
,
S
,
A
)
-
r
V
(
τ
,
S
,
A
)
=
0
,
where A
t is defined as A
t
=
I
t
/
t. But PDE (78) and (90) are equivalent for pricing fixed strike European Asian options. Second, since the PDEs are different, PDE (78) and (90) must have different boundary conditions. We do not know the boundary conditions used by Zvan et al. [17]. The forward shooting grid method used by Barraquand and Pudet [16] proceeds without any boundary conditions. Third, our PDE (78) is defined on the domain D
=
{
(
t
,
S
,
I
)
:
0
≤
t
≤
T
,
0
≤
S
≤
S
max
⁡
,
0
≤
I
≤
K
T
}, and I
max
⁡
=
K
T is the boundary and I
=
0 is the point where the options are valued. So there is no point that is less important on the interval [
0
,
K
T
]. Hence our numerical schemes use uniform grid spacing, while Zvan et al. [17] use nonuniform grid spacing of 41
×
45 on domain [
0,300
]
×
[
0,300
]. To achieve an accuracy of at least 0.1% of S, our grids are much finer than Zvan et al. [17].


### 3.2. Numerical Computation of Floating Strike Asian Options

To value the arithmetic average floating strike Asian call (put) options, we need to solve (67) (see (69)) for H
(
t
,
R
), and the values of the options will be V
(
t
,
S
,
I
)
=
S
H
(
t
,
R
) with R
=
I
/
S. Actually, what we really care about is the option values at time t
=
0, while, from the definition of I
t, it is obvious that I
0
=
0 and therefore R
0
=
0, so our goal of the computation is to find the prices of the Asian options at t
=
0; that is, V
(
0
,
S
,
0
)
=
S
H
(
t
,
0
).To solve (67), we use the fully implicit positive coefficient discretization, which will ensure convergence to the viscosity solution.Letτ
=
T
-
t, and then the PDE in (67) can be changed into the following forward equation in time:
(91)
-
∂
τ
H
(
τ
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
τ
,
R
)
+
(
1
-
r
R
)
H
R
(
τ
,
R
)
=
0
.Now, let us consider the discretization of (91). Define a grid 0
=
R
0
<
R
1
<
⋯
<
R
M
=
R
max
⁡ and a set of timesteps 0
=
τ
0
<
τ
1
<
⋯
<
τ
N
=
T. Let the discretization parameters be given by
(92)
Δ
R
i
=
R
i
+
1
-
R
i
,
i
=
0
,
…
,
M
;
we use even discretization in time; that is, τ
n
+
1
-
τ
n
=
Δ
τ, n
=
0
,
…
,
N
-
1.LetH
i
n be the approximate values of the solution at τ
n,  R
i; that is, H
i
n
≃
H
(
τ
n
,
R
i
). If 1
-
r
R
i
≥
0, we use forward difference to term H
R
(
τ
,
R
). Otherwise, we use backward difference to term H
R
(
τ
,
R
).The general form of the discretized PDE (91) is
(93)
α
i
n
+
1
H
i
-
1
n
+
1
+
γ
i
n
+
1
H
i
n
+
1
+
β
i
n
+
1
H
i
+
1
n
+
1
=
1
Δ
τ
H
i
n
,
with
(94)
α
i
n
+
1
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
Δ
R
i
+
1
,
β
i
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
+
1
2
-
1
-
r
R
i
Δ
R
i
+
1
,
kkkkkkkkkk
if
1
-
r
R
i
≥
0
,
α
i
n
+
1
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
Δ
R
i
+
1
-
r
R
i
-
1
Δ
R
i
,
β
i
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
+
1
2
,
ki
if
1
-
r
R
i
<
0
,
γ
i
n
+
1
=
-
α
i
n
+
1
-
β
i
n
+
1
+
1
Δ
τ
,
where
(95)
σ
^
i
n
+
1
=
argmax
σ
_
≤
σ
i
n
+
1
≤
σ
¯
⁡
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
k
k
k
k
k
k
k
k
k
k
i
+
H
i
+
1
n
Δ
R
i
+
1
2
]
.Of course, we haveα
i
n
+
1
≤
0
,
β
i
n
+
1
≤
0 and γ
i
n
+
1
≥
0. So the discretization equation (93) satisfies the positive coefficient condition of Forsyth and Labahn [41].Let(96)
H
n
=
[
H
1
n
,
…
,
H
M
-
1
n
]
T
,
σ
n
+
1
=
[
σ
1
n
+
1
,
…
,
σ
M
-
1
n
+
1
]
T
,
f
n
+
1
=
[
α
1
n
+
1
H
0
n
+
1
,
0
,
…
,
0
,
β
M
-
1
n
+
1
H
M
n
+
1
]
M
-
1
T
,
A
n
+
1
=
A
n
+
1
(
σ
i
n
+
1
)
=
[
γ
1
n
+
1
β
1
n
+
1
α
2
n
+
1
γ
2
n
+
1
β
2
n
+
1
⋱
⋱
⋱
α
M
-
1
n
+
1
γ
M
-
1
n
+
1
]
(
M
-
1
)
×
(
M
-
1
)
.
Then (93) can be written as the following matrix form
(97)
A
n
+
1
H
n
+
1
=
1
Δ
τ
H
n
-
f
n
+
1
,
σ
^
n
+
1
=
argmax
σ
_
≤
σ
n
+
1
≤
σ
¯
{
[
A
n
+
1
H
n
+
1
+
f
n
+
1
]
}
.By the boundary condition(98)
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
we replace H
0
n
+
1 in f
n
+
1 with (
Δ
τ
/
(
Δ
R
1
+
Δ
τ
)
)
H
1
n
+
1
+
(
Δ
R
1
/
(
Δ
R
1
+
Δ
τ
)
)
H
0
n. Then (97) can be written as
(99)
A
¯
n
+
1
H
n
+
1
=
1
Δ
τ
H
n
-
f
¯
n
+
1
,
σ
^
n
+
1
=
argmax
σ
_
≤
σ
n
+
1
≤
σ
¯
{
[
A
¯
n
+
1
H
n
+
1
+
f
¯
n
+
1
]
}
,
with
(100)
f
¯
n
+
1
=
[
α
1
n
+
1
Δ
R
1
Δ
R
1
+
Δ
τ
H
0
n
,
0
,
…
,
0
,
β
M
-
1
n
+
1
H
M
n
+
1
]
M
-
1
T
,
A
¯
n
+
1
=
A
¯
n
+
1
(
σ
n
+
1
)
=
[
γ
1
n
+
1
+
α
1
n
+
1
Δ
τ
Δ
R
1
+
Δ
τ
β
1
n
+
1
α
2
n
+
1
γ
2
n
+
1
β
2
n
+
1
⋱
⋱
⋱
α
M
-
1
n
+
1
γ
M
-
1
n
+
1
]
(
M
-
1
)
×
(
M
-
1
)
.It is clear thatA
¯
n
+
1 is M-matrix. From the property of M-matrix, the solution of (99) exists and is unique.It is important to ensure that we generate a numerical solution which converges to the viscosity solution. It has been shown that (67) satisfies the strong comparison property [42–44]. Then, from Barles and Souganidis [45] and Barles [46], a numerical scheme converges to the viscosity solution if the method is consistent, stable, and monotone. Thus, we will show that our numerical scheme satisfies these conditions.Lemma 6 (Stability).
The discretization (93) is stable.Proof.
It follows from (93) that
(101)
γ
i
n
+
1
|
H
i
n
+
1
|
≤
1
Δ
τ
|
H
i
n
|
-
α
i
n
+
1
|
H
i
-
1
n
+
1
|
-
β
i
n
+
1
|
H
i
+
1
n
+
1
|
≤
1
Δ
τ
∥
H
n
∥
∞
-
(
α
i
n
+
1
+
β
i
n
+
1
)
∥
H
n
+
1
∥
∞
,
since α
i
n
+
1
≤
0
,
β
i
n
+
1
≤
0.
Denote(102)
H
¯
n
+
1
=
[
H
0
n
+
1
,
…
,
H
M
n
+
1
]
T
.
If ∥
H
¯
n
+
1
∥
∞
=
|
H
j
n
+
1
|
,
0
<
j
<
M, then we have
(103)
∥
H
¯
n
+
1
∥
∞
≤
(
1
/
Δ
τ
)
∥
H
n
∥
∞
γ
i
n
+
1
+
α
i
n
+
1
+
β
i
n
+
1
=
∥
H
n
∥
∞
.
Otherwise, ∥
H
¯
n
+
1
∥
∞
=
|
H
0
n
+
1
| or ∥
H
¯
n
+
1
∥
∞
=
|
H
M
n
+
1
|. Therefore,
(104)
∥
H
¯
n
+
1
∥
∞
≤
max
⁡
(
∥
H
0
∥
∞
,
|
H
0
n
|
,
|
H
N
n
|
)
with H
0
n,  H
N
n being the given boundary conditions.Lemma 7 (Monotonicity).
The discretization (93) is monotonic.Proof.
The lemma is trivially true on the boundary, so let us just see the cases0
<
i
<
M and 0
<
n
≤
N. Denote
(105)
G
i
n
+
1
=
H
i
n
+
1
-
H
i
n
Δ
τ
+
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
{
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
H
i
+
1
n
+
1
kkkkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
}
.
For anyϵ
>
0,
(106)
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
+
ϵ
,
H
i
-
1
n
+
1
,
H
i
n
)
-
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
=
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
(
H
i
+
1
n
+
1
+
ϵ
)
kkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
)
-
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
H
i
+
1
n
+
1
kkkkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
)
≤
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
β
i
n
+
1
(
σ
i
n
+
1
)
ϵ
)
≤
0
,
since β
i
n
+
1
(
σ
i
n
+
1
)
≤
0 and sup
⁡
σ
_
≤
σ
≤
σ
¯
G
1
(
σ
)
-
sup
⁡
σ
_
≤
σ
≤
σ
¯
G
2
(
σ
)
≤
sup
⁡
σ
_
≤
σ
≤
σ
¯
(
G
1
(
σ
)
-
G
2
(
σ
)
).
Similarly,(107)
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
+
ϵ
,
H
i
n
)
≤
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
,
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
+
ϵ
)
≤
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
.Lemma 8 (Consistency).
The discretization (93) is consistent.Proof.
Supposeϕ
(
τ
,
R
) is a smooth test function with bounded derivatives of all orders with respect to (
τ
,
R
). Let
(108)
L
H
(
τ
,
R
)
=
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
τ
,
R
)
+
(
1
-
r
R
)
H
R
(
τ
,
R
)
.
Then
(109)
|
G
i
n
+
1
-
H
τ
(
τ
,
R
)
+
L
H
(
τ
,
R
)
|
≤
|
H
i
n
+
1
-
H
i
n
Δ
τ
-
H
τ
(
τ
,
R
)
|
+
|
L
H
(
τ
,
R
)
-
L
i
n
+
1
H
(
τ
,
R
)
|
,
where
(110)
L
i
n
+
1
H
(
τ
,
R
)
=
1
2
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
kkkkkkkikkkkkk
+
H
i
+
1
n
Δ
R
i
+
1
2
]
+
H
i
+
1
n
+
1
-
H
i
n
+
1
Δ
R
(
1
-
r
R
i
)
,
if 1
-
r
R
i
≥
0; otherwise
(111)
L
i
n
+
1
H
(
τ
,
R
)
=
1
2
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
kkkkkkkkkkkkkk
+
H
i
+
1
n
Δ
R
i
+
1
2
]
+
H
i
n
+
1
-
H
i
-
1
n
+
1
Δ
R
(
1
-
r
R
i
)
.
Using Taylor series expansions to (109), we have
(112)
|
G
i
n
+
1
-
H
τ
(
τ
,
R
)
+
L
H
(
τ
,
R
)
|
≤
O
(
Δ
τ
)
+
O
(
Δ
R
)
;
that is, (93) is consistent.Therefore, we have the following.Theorem 9.
The solution of the discretization scheme (93) converges uniformly to the unique viscosity solution of PDE (91).Table3 shows the convergence of the values of the arithmetic average floating strike Asian call and put options with the refinement of the grid. We use uniform grids. To ensure the desirable accuracy, we choose R
max
⁡
=
2. The execution times are less than 40 seconds.Table 3
Successive grid refinements demonstrating convergence of the results with the refinements whenr
=
0.10, σ
_
=
0.5, σ
-
=
1, σ
~
=
0.4, T
=
1, and tolerance
=
1
0
-
4. The values are the option price at S
=
100. Δ
R denotes the grid spacing, and Δ
τ denotes the time stepping. Execution times (in seconds) are in parentheses.
Δ
R
=
0.008
Δ
R
=
0.004
Δ
R
=
0.002
Δ
R
=
0.001
Δ
R
=
0.0005
Δ
τ
=
0.008
Δ
τ
=
0.004
Δ
t
=
0.002
Δ
t
=
0.001
Δ
t
=
0.0005
c
12.698
12.129
11.828
11.673
11.592
(1.00)
(2.19)
(9.08)
(35.66)
(141.38)
p
7.831
7.274
6.981
6.828
6.753
(1.01)
(2.06)
(9.05)
(34.10)
(143.12)Figure2 is the values of function H
(
0
,
R
) for the arithmetic average floating strike Asian call option with the given parameters specified under the figure. The corresponding option values are S
H
(
0,0
).Figure 2
TheH
(
0
,
R
) function of the arithmetic average floating strike Asian call option when r
=
0.10,  T
=
1,  σ
_
=
0.5,  σ
¯
=
1,  σ
~
=
0.4, and tolerance
=
1
0
-
4.


## 3.1. Numerical Computation of Fixed Strike Asian Options

An equivalent formulation of (48) in terms of the average A
t
=
(
1
/
t
)
∫
0
t
‍
S
τ
d
τ is given in Barraquand and Pudet [16] in certain volatility model with volatility σ
=
σ
~
σ
¯:
(77)
∂
t
V
(
t
,
S
,
A
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
t
,
S
,
A
)
+
r
S
V
S
(
t
,
S
,
A
)
+
S
-
A
t
V
A
(
t
,
S
,
A
)
-
r
V
(
t
,
S
,
A
)
=
0
,
V
(
T
,
S
,
A
)
=
(
A
-
K
)
+
.Equation (77) is convection dominated in the A direction because there is no diffusion effect in this dimension. Barraquand and Pudet [16] find that an explicit centrally weighted scheme for (77) is unstable. In particular, the convective term in the A dimension becomes very large as t
→
0. Barraquand and Pudet also note that implicit centrally weighted schemes will generally produce unsatisfactory results because of the numerical diffusion introduced by this first order accurate in time scheme.Zvan et al. [17] point out that under the condition of convection dominated, solutions generated by using a centrally weighted scheme on the convective term for (77) cannot be ensured to be free of oscillations. And solutions generated by the first order upstream weighting scheme on the convective term are no longer oscillatory, but their profiles are too diffuse. To produce oscillation free solutions without the excessive diffusion of first order upstream weighting, Zvan et al. [17] suggest the nonlinear van Leer flux limiter in conjunction with the Crank-Nicolson scheme. Since the flux limiter is nonlinear, the solutions should be obtained by using full Newton iteration at each time level.The ADI method, pioneered in the United States by Douglas, Rachford, Peaceman, Gunn, and others, has a number of advantages. First, explicit difference methods are rarely used to solve initial boundary value problems because of their poor stability. Implicit methods have superior stability properties but unfortunately they are difficult to solve in two and more dimensions. Consequently, ADI methods become an alternative because they can be programmed by solving a simple trigonal system of equations. For the convection-dominated problems, we could use the exponentially fitted schemes in each leg of the approximate scheme to ensure the stability of the results [40]. Since Barraquand and Pudet [16] have used the forward shooting grid method and Zvan et al. [17] have used the van Leer flux limiter, we will test the ADI method in this paper by comparing our results with theirs.For convenience, let us first convert (55) to be the following forward equation in time by substituting t with τ
=
T
-
t:
(78)
-
∂
τ
V
(
τ
,
S
,
I
)
+
1
2
σ
~
2
S
2
σ
¯
2
V
S
S
(
τ
,
S
,
I
)
+
r
S
V
S
(
τ
,
S
,
I
)
+
S
V
I
(
τ
,
S
,
I
)
-
r
V
(
τ
,
S
,
I
)
=
0
,
V
(
0
,
S
,
I
)
=
(
I
T
-
K
)
+
,
V
(
τ
,
0
,
I
)
=
e
-
r
τ
(
I
T
-
K
)
+
,
V
(
t
,
S
max
⁡
,
I
)
=
max
⁡
{
(
I
T
-
K
)
e
-
r
τ
+
S
max
⁡
r
T
(
1
-
e
-
r
τ
)
,
0
}
,
V
(
τ
,
S
,
K
T
)
=
1
r
T
(
1
-
e
-
r
τ
)
S
.We letτ
n,  n
=
0,1
,
…
,
N, be a set of portion points in [
0
,
T
] satisfying 0
=
τ
0
<
τ
1
<
⋯
<
τ
N
=
T and use even partition of time; that is, τ
n
=
n
Δ
τ,  n
=
0,1
,
…
,
N. Define the partitions of space variables as 0
=
S
0
<
S
1
<
⋯
<
S
M
=
S
max
⁡ and 0
=
I
0
<
I
1
<
⋯
<
I
J
=
I
max
⁡.With ADI, we march from time leveln to time level n
+
1
/
2 and then from time level n
+
1
/
2 to time level n
+
1. In this case we use exponential fitting in all space variables and implicit Euler in time. The first leg is given by the scheme
(79)
-
V
i
,
j
n
+
1
/
2
-
V
i
,
j
n
(
1
/
2
)
Δ
τ
+
ρ
i
,
j
[
1
Δ
S
i
Δ
S
i
+
1
V
i
-
1
,
j
n
+
1
/
2
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
-
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
+
1
Δ
S
i
+
1
2
V
i
+
1
,
j
n
+
1
/
2
]
+
r
S
i
V
i
+
1
,
j
n
+
1
/
2
-
V
i
,
j
n
+
1
/
2
Δ
S
i
+
1
+
S
i
V
i
,
j
+
1
n
-
V
i
,
j
-
1
n
Δ
I
j
+
1
+
Δ
I
j
-
r
V
i
,
j
n
+
1
/
2
=
0
,
with
(80)
ρ
i
,
j
=
r
S
i
Δ
S
i
+
1
2
coth
⁡
(
r
S
i
Δ
S
i
+
1
σ
¯
2
σ
~
2
S
i
2
)
,
coth
⁡
(
x
)
=
e
2
x
+
1
e
2
x
-
1
.
Let
(81)
α
i
=
-
ρ
i
,
j
Δ
S
i
Δ
S
i
+
1
,
β
i
=
-
ρ
i
,
j
(
Δ
S
i
+
1
)
2
-
r
S
i
Δ
S
i
+
1
,
γ
i
=
-
α
i
-
β
i
+
r
.
Then the scheme (79) can be written as
(82)
α
i
V
i
-
1
,
j
n
+
1
/
2
+
(
γ
i
+
2
Δ
τ
)
V
i
,
j
n
+
1
/
2
+
β
i
V
i
+
1
,
j
n
+
1
/
2
=
2
Δ
τ
V
i
,
j
n
+
S
i
Δ
I
j
+
1
+
Δ
I
j
V
i
,
j
+
1
n
-
S
i
Δ
I
j
+
1
+
Δ
I
j
V
i
,
j
-
1
n
.
We let
(83)
V
j
n
=
[
V
1
,
j
n
,
V
2
,
j
n
,
…
,
V
M
-
1
,
j
n
]
T
,
V
j
n
+
1
/
2
=
[
V
1
,
j
n
+
1
/
2
,
V
2
,
j
n
+
1
/
2
,
…
,
V
M
-
1
,
j
n
+
1
/
2
]
T
,
f
j
n
+
1
/
2
=
[
α
1
V
0
,
j
n
+
1
/
2
,
0
,
…
,
0
,
β
M
-
1
V
M
,
j
n
+
1
/
2
]
T
,
r
j
n
=
[
(
2
Δ
τ
V
1
,
j
n
+
S
1
Δ
I
j
+
1
+
Δ
I
j
V
1
,
j
+
1
n
k
k
-
S
1
Δ
I
j
+
1
+
Δ
I
j
V
1
,
j
-
1
n
)
,
…
,
k
(
2
Δ
τ
V
M
-
1
,
j
n
+
S
M
-
1
Δ
I
j
+
1
+
Δ
I
j
V
M
-
1
,
j
+
1
n
k
.
k
-
S
M
-
1
Δ
I
j
+
1
+
Δ
I
j
V
M
-
1
,
j
-
1
n
)
]
T
,
A
1
=
[
γ
1
β
1
α
2
γ
2
β
2
⋱
⋱
⋱
α
M
-
1
γ
M
-
1
]
(
M
-
1
)
×
(
M
-
1
)
.
Then we can write (82) as the following equivalent matrix form
(84)
A
1
V
j
n
+
1
/
2
+
f
j
n
+
1
/
2
=
r
j
n
.The second leg of the discretized PDE in (78) is given by the scheme
(85)
-
V
i
,
j
n
+
1
-
V
i
,
j
n
+
1
/
2
(
1
/
2
)
Δ
τ
+
ρ
i
,
j
[
1
Δ
S
i
Δ
S
i
+
1
V
i
-
1
,
j
n
+
1
/
2
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
kk
-
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
kkkkkkkkkkkkkkkk
kki
kk
(
1
Δ
S
i
+
1
2
+
1
Δ
S
i
Δ
S
i
+
1
)
V
i
,
j
n
+
1
/
2
+
1
Δ
S
i
+
1
2
V
i
+
1
,
j
n
+
1
/
2
]
+
r
S
i
V
i
+
1
,
j
n
+
1
/
2
-
V
i
,
j
n
+
1
/
2
Δ
S
i
+
1
+
S
i
V
i
,
j
+
1
n
+
1
-
V
i
,
j
n
+
1
Δ
I
j
+
1
-
r
V
i
,
j
n
+
1
/
2
=
0
.
Let
(86)
δ
j
=
2
Δ
τ
+
S
i
Δ
I
j
+
1
,
ϵ
j
=
-
S
i
Δ
I
j
+
1
.
Then the scheme (85) can be written as
(87)
δ
j
V
i
,
j
n
+
1
-
ϵ
j
V
i
,
j
+
1
n
+
1
=
-
α
i
V
i
-
1
,
j
n
+
1
/
2
-
(
γ
i
-
2
Δ
τ
)
V
i
,
j
n
+
1
/
2
-
β
i
V
i
+
1
,
j
n
+
1
/
2
.We let(88)
V
i
n
+
1
=
[
V
i
,
0
n
+
1
,
V
i
,
1
n
+
1
,
…
,
V
i
,
J
-
1
n
+
1
]
T
,
V
i
n
+
1
/
2
=
[
V
i
,
0
n
+
1
/
2
,
V
i
,
1
n
+
1
/
2
,
…
,
V
i
,
J
-
1
n
+
1
/
2
]
T
,
F
i
=
[
0,0
,
…
,
0
,
ϵ
J
-
1
V
i
,
J
n
+
1
]
J
T
,
R
i
n
+
1
/
2
=
[
(
-
α
i
V
i
-
1,0
n
+
1
/
2

-

(
γ
i
-
2
Δ
τ
)
V
i
,
0
n
+
1
/
2

-

β
i
V
i
+
1,0
n
+
1
/
2
)
,
…
,
(
-
α
i
V
i
-
1
,
J
-
1
n
+
1
/
2
-
(
γ
i
-
2
Δ
τ
)
V
i
,
J
-
1
n
+
1
/
2
-
β
i
V
i
+
1
,
J
-
1
n
+
1
/
2
)
]
T
,
A
2
=
[
δ
0
-
ϵ
0
0
δ
1
-
ϵ
1
⋱
⋱
-
ϵ
J
-
2
δ
J
-
1
]
J
×
J
.
The matrix form of (87) is
(89)
A
2
V
i
n
+
1
+
F
i
n
+
1
=
R
i
n
+
1
/
2
.
The matrix equations (84) and (89) can be solved by using LU decomposition.The goal of the computation is to find the fair price of the Asian option at emission, that is,V
(
S
0
,
I
0
,
0
). We need to make sense of any value I
∈
[
0
,
I
max
⁡
] and stock price S
∈
[
0
,
S
max
⁡
] at time t
=
0. There are no problems with the stock prices, but the integral at time zero can only be 0 for continuous averages; that is, I
0
=
0.Figure1 is the results we obtained by using the above discretization scheme with σ
~
=
0.8,  σ
¯
=
0.5,  r
=
0.10,  T
=
1, and K
=
100. The point marked by * is the price V
(
0
,
S
,
0
) of the Asian option. We choose S
max
⁡
=
3
K to ensure the desirable accuracy.Figure 1
Fixed strike Asian call value whenσ
~
=
0.8,  σ
¯
=
0.5,  r
=
0.10,  T
=
1, and K
=
100.Table1 shows the convergence of the results with the refinement of the grid. From (78), the values of variable I lie in the interval [
0
,
K
T
], where K
T is the boundary and 0 is the point where the continuous averages Asian options are valued, so there is no point that is less important on the interval [
0
,
K
T
]. Hence we use uniform grids in direction I.Table 1
Successive grid refinements demonstrating convergence of the results with the refinements whenr
=
0.10, T
=
1, σ
~
=
0.8, σ
-
=
0.5, and K
=
95. The values are the option prices at S
=
100. Δ
S and Δ
I denote the grid spacing in directions S and I, respectively. Δ
t
=
0.04. Execution times (in seconds) are in parentheses.
Δ
I
=
0.95
Δ
I
=
0.475
Δ
I
=
0.2375
Δ
I
=
0.11875
Δ
I
=
0.059375
Δ
S
=
5.70
14.532
14.19
14.021
13.94
13.89
(1.01)
(2.10)
(7.081)
(27.225)
(90)
Δ
S
=
2.85
14.42
14.08
13.914
13.82
13.781
(3.08)
(7.13)
(23.012)
(76.51)
(292)
Δ
S
=
1.425
14.39
14.053
13.879
13.790
13.74
(10.22)
(24.71)
(76.01)
(238)
(940)In our uncertain volatility model, the pricing PDEs (48) are similar to the PDEs in certain volatility model with volatility σ
=
σ
~
σ
¯. Barraquand and Pudet [16] and Zvan et al. [17] give numerical results of Asian options with certain volatilities. In the uncertain volatility model, when σ
~
=
0.8,  σ
¯
=
0.5, the Asian option PDE is the same as that of the Asian option with certain volatility σ
=
0.4. When the parameters r
=
0.10,  T
=
1,  K
=
95, and grid Δ
I
=
0.11875,  Δ
S
=
2.85, the price of the Asian call option we calculate is 13.82. The results with the same parameters of Barraquand and Pudet [16] and Zvan et al. [17] are 13.825 and 13.721, respectively.Table2 is the comparison of our results with those of Zvan et al. [17]. Z, F, and V in Table 2 refer to the result of Zvan et al. [17]. Our results are denoted by Asian fixed strike, which calculated with Δ
τ
=
0.01,0.02,0.04,  Δ
I
=
0.059375,0.11875,0.2375, and Δ
S
=
2.85 for maturities of one quarter, half a year, and one year, respectively. Since I
max
⁡
=
K
T, the partitions in I direction are related with the maturity T. The execution time is much less than theirs. The grid spacing was chosen to achieve an accuracy of at least 0.10% of S. Our results are close to that of Zvan et al. [17] and the difference is less than 0.10% of S.Table 2
Fixed strike Asian call values whenr
=
0.10 and K
=
95. The values are the option prices at S
=
100. Our results are denoted by Asian fixed strike, which calculated with Δ
τ
=
0.01,0.02,0.04, Δ
I
=
0.059375,0.11875,0.2375, and Δ
S
=
2.85 for maturities of one quarter, half a year, and one year, respectively. Z, F, and V refer to the result of Zvan et al. [17]. Execution times (in seconds) are in parentheses.
σ
~
σ
-
T
Asian fixed strike
Z, F, and V
0.2
0.5
0.25
6.190 (23.0)
6.133
0.5
7.199 (24.0)
7.244
1
9.204 (23.6)
9.316
1
0.25
6.56 (23.0)
6.501
0.5
7.99 (24.1)
7.921
1
10.44 (25.0)
10.309
2
0.25
8.30 (24.1)
8.123
0.5
10.55 (24.4)
10.357
1
13.91 (25.2)
13.721Our calculation procedure is different from Zvan et al. [17], as well as Barraquand and Pudet [16], in the following four aspects. First, our pricing PDE (78) is in terms of the running sum I
t. Zvan et al. [17] and Barraquand and Pudet [16] use the pricing PDE
(90)
-
∂
τ
V
(
τ
,
S
,
A
)
+
1
2
S
2
σ
2
V
S
S
(
τ
,
S
,
A
)
+
r
S
V
S
(
τ
,
S
,
A
)
+
S
-
A
t
V
A
(
τ
,
S
,
A
)
-
r
V
(
τ
,
S
,
A
)
=
0
,
where A
t is defined as A
t
=
I
t
/
t. But PDE (78) and (90) are equivalent for pricing fixed strike European Asian options. Second, since the PDEs are different, PDE (78) and (90) must have different boundary conditions. We do not know the boundary conditions used by Zvan et al. [17]. The forward shooting grid method used by Barraquand and Pudet [16] proceeds without any boundary conditions. Third, our PDE (78) is defined on the domain D
=
{
(
t
,
S
,
I
)
:
0
≤
t
≤
T
,
0
≤
S
≤
S
max
⁡
,
0
≤
I
≤
K
T
}, and I
max
⁡
=
K
T is the boundary and I
=
0 is the point where the options are valued. So there is no point that is less important on the interval [
0
,
K
T
]. Hence our numerical schemes use uniform grid spacing, while Zvan et al. [17] use nonuniform grid spacing of 41
×
45 on domain [
0,300
]
×
[
0,300
]. To achieve an accuracy of at least 0.1% of S, our grids are much finer than Zvan et al. [17].


## 3.2. Numerical Computation of Floating Strike Asian Options

To value the arithmetic average floating strike Asian call (put) options, we need to solve (67) (see (69)) for H
(
t
,
R
), and the values of the options will be V
(
t
,
S
,
I
)
=
S
H
(
t
,
R
) with R
=
I
/
S. Actually, what we really care about is the option values at time t
=
0, while, from the definition of I
t, it is obvious that I
0
=
0 and therefore R
0
=
0, so our goal of the computation is to find the prices of the Asian options at t
=
0; that is, V
(
0
,
S
,
0
)
=
S
H
(
t
,
0
).To solve (67), we use the fully implicit positive coefficient discretization, which will ensure convergence to the viscosity solution.Letτ
=
T
-
t, and then the PDE in (67) can be changed into the following forward equation in time:
(91)
-
∂
τ
H
(
τ
,
R
)
+
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
τ
,
R
)
+
(
1
-
r
R
)
H
R
(
τ
,
R
)
=
0
.Now, let us consider the discretization of (91). Define a grid 0
=
R
0
<
R
1
<
⋯
<
R
M
=
R
max
⁡ and a set of timesteps 0
=
τ
0
<
τ
1
<
⋯
<
τ
N
=
T. Let the discretization parameters be given by
(92)
Δ
R
i
=
R
i
+
1
-
R
i
,
i
=
0
,
…
,
M
;
we use even discretization in time; that is, τ
n
+
1
-
τ
n
=
Δ
τ, n
=
0
,
…
,
N
-
1.LetH
i
n be the approximate values of the solution at τ
n,  R
i; that is, H
i
n
≃
H
(
τ
n
,
R
i
). If 1
-
r
R
i
≥
0, we use forward difference to term H
R
(
τ
,
R
). Otherwise, we use backward difference to term H
R
(
τ
,
R
).The general form of the discretized PDE (91) is
(93)
α
i
n
+
1
H
i
-
1
n
+
1
+
γ
i
n
+
1
H
i
n
+
1
+
β
i
n
+
1
H
i
+
1
n
+
1
=
1
Δ
τ
H
i
n
,
with
(94)
α
i
n
+
1
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
Δ
R
i
+
1
,
β
i
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
+
1
2
-
1
-
r
R
i
Δ
R
i
+
1
,
kkkkkkkkkk
if
1
-
r
R
i
≥
0
,
α
i
n
+
1
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
Δ
R
i
+
1
-
r
R
i
-
1
Δ
R
i
,
β
i
=
-
σ
~
2
R
i
2
(
σ
^
i
n
+
1
)
2
2
Δ
R
i
+
1
2
,
ki
if
1
-
r
R
i
<
0
,
γ
i
n
+
1
=
-
α
i
n
+
1
-
β
i
n
+
1
+
1
Δ
τ
,
where
(95)
σ
^
i
n
+
1
=
argmax
σ
_
≤
σ
i
n
+
1
≤
σ
¯
⁡
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
k
k
k
k
k
k
k
k
k
k
i
+
H
i
+
1
n
Δ
R
i
+
1
2
]
.Of course, we haveα
i
n
+
1
≤
0
,
β
i
n
+
1
≤
0 and γ
i
n
+
1
≥
0. So the discretization equation (93) satisfies the positive coefficient condition of Forsyth and Labahn [41].Let(96)
H
n
=
[
H
1
n
,
…
,
H
M
-
1
n
]
T
,
σ
n
+
1
=
[
σ
1
n
+
1
,
…
,
σ
M
-
1
n
+
1
]
T
,
f
n
+
1
=
[
α
1
n
+
1
H
0
n
+
1
,
0
,
…
,
0
,
β
M
-
1
n
+
1
H
M
n
+
1
]
M
-
1
T
,
A
n
+
1
=
A
n
+
1
(
σ
i
n
+
1
)
=
[
γ
1
n
+
1
β
1
n
+
1
α
2
n
+
1
γ
2
n
+
1
β
2
n
+
1
⋱
⋱
⋱
α
M
-
1
n
+
1
γ
M
-
1
n
+
1
]
(
M
-
1
)
×
(
M
-
1
)
.
Then (93) can be written as the following matrix form
(97)
A
n
+
1
H
n
+
1
=
1
Δ
τ
H
n
-
f
n
+
1
,
σ
^
n
+
1
=
argmax
σ
_
≤
σ
n
+
1
≤
σ
¯
{
[
A
n
+
1
H
n
+
1
+
f
n
+
1
]
}
.By the boundary condition(98)
∂
t
H
(
t
,
R
)
+
∂
R
H
(
t
,
R
)
=
0
,
R
=
0
,
we replace H
0
n
+
1 in f
n
+
1 with (
Δ
τ
/
(
Δ
R
1
+
Δ
τ
)
)
H
1
n
+
1
+
(
Δ
R
1
/
(
Δ
R
1
+
Δ
τ
)
)
H
0
n. Then (97) can be written as
(99)
A
¯
n
+
1
H
n
+
1
=
1
Δ
τ
H
n
-
f
¯
n
+
1
,
σ
^
n
+
1
=
argmax
σ
_
≤
σ
n
+
1
≤
σ
¯
{
[
A
¯
n
+
1
H
n
+
1
+
f
¯
n
+
1
]
}
,
with
(100)
f
¯
n
+
1
=
[
α
1
n
+
1
Δ
R
1
Δ
R
1
+
Δ
τ
H
0
n
,
0
,
…
,
0
,
β
M
-
1
n
+
1
H
M
n
+
1
]
M
-
1
T
,
A
¯
n
+
1
=
A
¯
n
+
1
(
σ
n
+
1
)
=
[
γ
1
n
+
1
+
α
1
n
+
1
Δ
τ
Δ
R
1
+
Δ
τ
β
1
n
+
1
α
2
n
+
1
γ
2
n
+
1
β
2
n
+
1
⋱
⋱
⋱
α
M
-
1
n
+
1
γ
M
-
1
n
+
1
]
(
M
-
1
)
×
(
M
-
1
)
.It is clear thatA
¯
n
+
1 is M-matrix. From the property of M-matrix, the solution of (99) exists and is unique.It is important to ensure that we generate a numerical solution which converges to the viscosity solution. It has been shown that (67) satisfies the strong comparison property [42–44]. Then, from Barles and Souganidis [45] and Barles [46], a numerical scheme converges to the viscosity solution if the method is consistent, stable, and monotone. Thus, we will show that our numerical scheme satisfies these conditions.Lemma 6 (Stability).
The discretization (93) is stable.Proof.
It follows from (93) that
(101)
γ
i
n
+
1
|
H
i
n
+
1
|
≤
1
Δ
τ
|
H
i
n
|
-
α
i
n
+
1
|
H
i
-
1
n
+
1
|
-
β
i
n
+
1
|
H
i
+
1
n
+
1
|
≤
1
Δ
τ
∥
H
n
∥
∞
-
(
α
i
n
+
1
+
β
i
n
+
1
)
∥
H
n
+
1
∥
∞
,
since α
i
n
+
1
≤
0
,
β
i
n
+
1
≤
0.
Denote(102)
H
¯
n
+
1
=
[
H
0
n
+
1
,
…
,
H
M
n
+
1
]
T
.
If ∥
H
¯
n
+
1
∥
∞
=
|
H
j
n
+
1
|
,
0
<
j
<
M, then we have
(103)
∥
H
¯
n
+
1
∥
∞
≤
(
1
/
Δ
τ
)
∥
H
n
∥
∞
γ
i
n
+
1
+
α
i
n
+
1
+
β
i
n
+
1
=
∥
H
n
∥
∞
.
Otherwise, ∥
H
¯
n
+
1
∥
∞
=
|
H
0
n
+
1
| or ∥
H
¯
n
+
1
∥
∞
=
|
H
M
n
+
1
|. Therefore,
(104)
∥
H
¯
n
+
1
∥
∞
≤
max
⁡
(
∥
H
0
∥
∞
,
|
H
0
n
|
,
|
H
N
n
|
)
with H
0
n,  H
N
n being the given boundary conditions.Lemma 7 (Monotonicity).
The discretization (93) is monotonic.Proof.
The lemma is trivially true on the boundary, so let us just see the cases0
<
i
<
M and 0
<
n
≤
N. Denote
(105)
G
i
n
+
1
=
H
i
n
+
1
-
H
i
n
Δ
τ
+
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
{
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
H
i
+
1
n
+
1
kkkkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
}
.
For anyϵ
>
0,
(106)
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
+
ϵ
,
H
i
-
1
n
+
1
,
H
i
n
)
-
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
=
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
(
H
i
+
1
n
+
1
+
ϵ
)
kkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
)
-
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
α
i
n
+
1
(
σ
i
n
+
1
)
H
i
-
1
n
+
1
+
β
i
n
+
1
(
σ
i
n
+
1
)
H
i
+
1
n
+
1
kkkkkkkkk
+
γ
i
n
+
1
(
σ
i
n
+
1
)
H
i
n
+
1
)
≤
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
(
β
i
n
+
1
(
σ
i
n
+
1
)
ϵ
)
≤
0
,
since β
i
n
+
1
(
σ
i
n
+
1
)
≤
0 and sup
⁡
σ
_
≤
σ
≤
σ
¯
G
1
(
σ
)
-
sup
⁡
σ
_
≤
σ
≤
σ
¯
G
2
(
σ
)
≤
sup
⁡
σ
_
≤
σ
≤
σ
¯
(
G
1
(
σ
)
-
G
2
(
σ
)
).
Similarly,(107)
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
+
ϵ
,
H
i
n
)
≤
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
,
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
+
ϵ
)
≤
G
i
n
+
1
(
Δ
R
¯
,
H
i
n
+
1
,
H
i
+
1
n
+
1
,
H
i
-
1
n
+
1
,
H
i
n
)
.Lemma 8 (Consistency).
The discretization (93) is consistent.Proof.
Supposeϕ
(
τ
,
R
) is a smooth test function with bounded derivatives of all orders with respect to (
τ
,
R
). Let
(108)
L
H
(
τ
,
R
)
=
sup
⁡
σ
_
2
≤
σ
2
≤
σ
¯
2
1
2
σ
~
2
R
2
σ
2
H
R
R
(
τ
,
R
)
+
(
1
-
r
R
)
H
R
(
τ
,
R
)
.
Then
(109)
|
G
i
n
+
1
-
H
τ
(
τ
,
R
)
+
L
H
(
τ
,
R
)
|
≤
|
H
i
n
+
1
-
H
i
n
Δ
τ
-
H
τ
(
τ
,
R
)
|
+
|
L
H
(
τ
,
R
)
-
L
i
n
+
1
H
(
τ
,
R
)
|
,
where
(110)
L
i
n
+
1
H
(
τ
,
R
)
=
1
2
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
kkkkkkkikkkkkk
+
H
i
+
1
n
Δ
R
i
+
1
2
]
+
H
i
+
1
n
+
1
-
H
i
n
+
1
Δ
R
(
1
-
r
R
i
)
,
if 1
-
r
R
i
≥
0; otherwise
(111)
L
i
n
+
1
H
(
τ
,
R
)
=
1
2
sup
⁡
σ
_
≤
σ
i
n
+
1
≤
σ
¯
σ
~
2
R
i
2
(
σ
i
n
+
1
)
2
k
k
k
k
k
k
k
k
k
×
[
H
i
-
1
n
+
1
Δ
R
i
Δ
R
i
+
1
-
(
H
i
n
+
1
Δ
R
i
+
1
2
+
H
i
n
+
1
Δ
R
i
Δ
R
i
+
1
)
kkkkkkkkkkkkkk
+
H
i
+
1
n
Δ
R
i
+
1
2
]
+
H
i
n
+
1
-
H
i
-
1
n
+
1
Δ
R
(
1
-
r
R
i
)
.
Using Taylor series expansions to (109), we have
(112)
|
G
i
n
+
1
-
H
τ
(
τ
,
R
)
+
L
H
(
τ
,
R
)
|
≤
O
(
Δ
τ
)
+
O
(
Δ
R
)
;
that is, (93) is consistent.Therefore, we have the following.Theorem 9.
The solution of the discretization scheme (93) converges uniformly to the unique viscosity solution of PDE (91).Table3 shows the convergence of the values of the arithmetic average floating strike Asian call and put options with the refinement of the grid. We use uniform grids. To ensure the desirable accuracy, we choose R
max
⁡
=
2. The execution times are less than 40 seconds.Table 3
Successive grid refinements demonstrating convergence of the results with the refinements whenr
=
0.10, σ
_
=
0.5, σ
-
=
1, σ
~
=
0.4, T
=
1, and tolerance
=
1
0
-
4. The values are the option price at S
=
100. Δ
R denotes the grid spacing, and Δ
τ denotes the time stepping. Execution times (in seconds) are in parentheses.
Δ
R
=
0.008
Δ
R
=
0.004
Δ
R
=
0.002
Δ
R
=
0.001
Δ
R
=
0.0005
Δ
τ
=
0.008
Δ
τ
=
0.004
Δ
t
=
0.002
Δ
t
=
0.001
Δ
t
=
0.0005
c
12.698
12.129
11.828
11.673
11.592
(1.00)
(2.19)
(9.08)
(35.66)
(141.38)
p
7.831
7.274
6.981
6.828
6.753
(1.01)
(2.06)
(9.05)
(34.10)
(143.12)Figure2 is the values of function H
(
0
,
R
) for the arithmetic average floating strike Asian call option with the given parameters specified under the figure. The corresponding option values are S
H
(
0,0
).Figure 2
TheH
(
0
,
R
) function of the arithmetic average floating strike Asian call option when r
=
0.10,  T
=
1,  σ
_
=
0.5,  σ
¯
=
1,  σ
~
=
0.4, and tolerance
=
1
0
-
4.



---

*Source: 786391-2014-06-16.xml*