Front. Math.

China 2016, 11(1): 237–254

DOI 10.1007/s11464-015-0500-0

Lattice Boltzmann methods for

solving partial differential equations of
exotic option pricing
Zhiqiang ZHOU, Jingtang MA
School of Economic Mathematics, Southwestern University of Finance and Economics,
Chengdu 611130, China


c Higher Education Press and Springer-Verlag Berlin Heidelberg 2015

Abstract This paper establishes a lattice Boltzmann method (LBM) with two
amending functions for solving partial differential equations (PDEs) arising in
Asian and lookback options pricing. The time evolution of stock prices can be
regarded as the movement of randomizing particles in different directions, and
the discrete scheme of LBM can be interpreted as the binomial models. With
the Chapman-Enskog multi-scale expansion, the PDEs are recovered correctly
from the continuous Boltzmann equation and the computational complexity is
O(N ), where N is the number of space nodes. Compared to the traditional
LBM, the coefficients of equilibrium distribution and amending functions are
taken as polynomials instead of constants. The stability of LBM is studied
via numerical examples and numerical comparisons show that the LBM is as
accurate as the existing numerical methods for pricing the exotic options and
takes much less CPU time.
Keywords Exotic option pricing, lattice Boltzmann method, Chapman-Enskog
multi-scale expansion, stability, computational complexity
MSC 91G20, 91G60, 91G80

1 Introduction

An option is a contract which gives the buyer (the owner) the right, but not
the obligation, to buy or sell an underlying asset or instrument at a specified
strike price on or before a specified date. The pricing of options can date back
to well-known Black-Scholes-Merton model ([3,16]). Since then, there have
been numerous papers studying the option pricing methods including binomial
methods, finite difference methods, Monte Carlo methods etc. (see a survey
Received February 4, 2015; accepted August 30, 2015
Corresponding author: Jingtang MA, E-mail: mjt@swufe.edu.cn
238 Zhiqiang ZHOU, Jingtang MA

paper [5]).
Exotic option is a class of options which have features making it more
complex than commonly traded vanilla options. Asian options and lookback
options are two classical examples of exotic options. In risk-neutral sense,
the value of such kinds of exotic options can be written as the solutions of
partial differential equations (PDEs). In general, these PDEs cannot be solved
analytically and usually finite difference methods (FDMs) are used to solve
them (see, e.g., [4,7,9,11,17,21,22]). This paper studies a new method—lattice
Boltzmann method for solving the pricing PDEs of Asian options and lookback
options. The idea of lattice Boltzmann method originates from and has many
applications in computational physics.
Over the past thirty years, the lattice Boltzmann method (LBM) has
achieved much success in studying nonlinear equations and simulating complex
fluid systems (see, e.g., [1,6,12–15,19,20]). The LBM is a new technique based
on the continuous Boltzmann equation: through density distribution functions
and the Chapman-Enskog multi-scale expansion, the macroscopic variables are
derived and the macroscopic equations are recovered. Different from traditional
finite difference methods, LBM regards time evolution behaviors as a special
movement of amount of particles, and the solution process does not need to solve
algebraic system at each time level, which makes the computational complexity
to be only O(N ) for each time level and O(N ) for the memory. Moreover, the
simplicity of programming and intrinsic parallelism are also the notable
advantages for LBM. In addition, the lattice Boltzmann model has its own
financial meaning for the option pricing. To the best of our knowledge, there
has been no reference exploring LBM in exotic option pricing.
The major contributions are as follows: (i) correctly recover the pricing
PDEs from LBM equations by designing special amending functions and
equilibrium distribution functions; (ii) explore the relationship between lattice
Boltzmann model and binomial option pricing model.
The remaining parts of the paper are arranged as follows. In Section 2, we
describe the PDEs for the Asian options and lookback options. In Section 3,
we establish the continuous Boltzmann model and deduce the parameters to
recover macroscopic PDEs of option pricing. In Section 4, we deal with the
boundary conditions and give the implementation of LBM. Numerical results
demonstrating the efficiency, computation cost and stability of LBM are
presented in Section 5. Concluding remarks are given in the final section.

2 PDEs for Asian options and lookback options

Let the stock price S = S(t) follow a geometric Brownian motion in the risk-
neutral sense:
dS(t) = rS(t)dt + σS(t)dW (t),

where r denotes the risk free interest rate, σ the volatility, and dW (t) the
Lattice Boltzmann methods for solving PDEs of exotic option pricing 239

standard Brownian motion.

Denote  t
I(t) = S(τ )dτ.
0
Then the value of the continuous arithmetic average Asian options with payoff
∗)
max( I(T ∗
T ∗ − K, 0), where K is the fixed strike and T is the expiry time of the
option, which is given by

∗ −t)
  I(T ∗ ) 
C(S(t), I(t), t) = e−r(T Et max − K, 0 ,
T∗
where Et denotes the conditional expectation at t. As is well known (see, e.g.,
[11]), the value of the Asian option is formulated to satisfy the following PDE:

∂C 1 ∂2C ∂C ∂C
+ σ 2 S 2 2 + rS +S − rC = 0, (1)
∂t 2 ∂S ∂S ∂I
where C = C(S, I, t) with S and I being dummy variables. The terminal
condition at expiry time T ∗ is given by
 I 
C(S, I, T ∗ ) = max − K, 0 . (2)
T∗
Due to the degeneracy of the PDE (1) in the I direction, it can be verified that
in the region I/T ∗  K, for all t  T ∗ ,

S  I 
−r(T ∗ −t) ∗
C(S, I, t) = ∗
(1 − e ) + ∗
− K e−r(T −t) (3)
rT T
satisfies the PDE (1) (see, e.g., [7,8]). Therefore, we only need to consider the
case I/T ∗  K. Using the transform
I
K− T∗
C(S, I, t) = Su(x, t), x= , (4)
S
Rogers and Shi [17] reduced the PDE (1) with I/T ∗  K for the Asian option
into the following one dimensional PDE:

∂  1 ∂ 1 ∂2
u(x, t) + − rx − ∗ u(x, t) + σ 2 x2 2 u(x, t) = 0,
∂t T ∂x 2 ∂x
x ∈ [0, +∞), 0  t  T ∗ , (5)

with terminal and boundary conditions

⎧
⎪ ∂u
⎨ u(x, T ∗ ) = 0, lim = 0,
x→∞ ∂x
(6)
⎪
⎩ u(0, t) = 1 (1 − e−r(T ∗ −t) ).
rT ∗
240 Zhiqiang ZHOU, Jingtang MA

Now, we describe the PDE for the lookback options. Define

J(t) := max (τ ).
0τ t

Then a lookback option is defined as a European option with payoff max(J(T ∗ )−

S(T ∗ ), 0), where T ∗ is the expiry date. Define a transform

S
C(S, J, t) = Ju(x, t), x= . (7)
J
Then function u satisfies the following terminal-boundary problem of PDE (see,
e.g., [4,9,22]):

∂ ∂ 1 ∂2
u(x, t) − ru + rx u(x, t) + σ 2 x2 2 u(x, t) = 0, x ∈ [0, 1], 0  t  T ∗ ,
∂t ∂x 2 ∂x
(8)
with the terminal and boundary conditions
⎧
⎪ u(x, T ∗ ) = max(1 − x, 0),
⎪
⎨ ∗
u(0, t) = e−r(T −t) , (9)
⎪
⎪ 
⎩ ∂u  = u(1, t).
∂x x=1
In next section, we will develop lattice Boltzmann methods for solving the
terminal-boundary value problem (5) and (6) for pricing the Asian options and
(8) and (9) for pricing the lookback options.

3 Lattice Boltzmann methods

The key procedures of general LBM include four steps (see, e.g., [6,19]): set
an appropriate continuous Boltzmann equation, recover the macroscopic PDEs
by Taylor’s expansion and Chapman-Enskog multi-scale expansion techniques,
deal with the boundary conditions carefully, and then adjust parameters in the
model to ensue computational stability. We follow the above steps to develop
the LBM for Asian option pricing problem (5) and (6) and lookback option
pricing problem (8) and (9).
3.1 Continuous Boltzmann models
We establish the lattice BGK model (see, e.g., Bhatnager, Gross, and Krook
[2]) for Asian option pricing problem (5) and (6) and lookback option pricing
problem (8) and (9). We discretize the one-dimensional space into linear well-
proportioned grids, and discretize the velocity (regarding time variable) into
five directions, which are defined as

[e1 , e2 , e3 , e4 , e5 ] = [0, 1, −1, 2, −2].

Lattice Boltzmann methods for solving PDEs of exotic option pricing 241

The time-space lattice BGK model is defined as

fi (x + cεei , t − ε2 T ) − fi (x, t)
1 (0)
= − (fi (x, t) − fi (x, t)) + εhi (x, t) + ε2 gi (x, t), (10)
τ
(0)
where fi (x, t) and fi (x, t) are defined as density distribution functions (DFs)
and local equilibrium distribution functions (EDFs), respectively, τ is the
dimensionless relaxation time, c is a constant, T is the characteristic time,
and hi (x, t) and gi (x, t) are amending functions. We regard the small
parameter cε as the space step and ε2 T the time step. In the classic kinetic
theory, ε is called the Knudsen number and it is defined as ε = /L,  is the
mean free path and L is the characteristic length.
Similar to the general LBM, the macroscopic variable u satisfies the
following conservation law:
5
 5
 (0)
fi (x, t) = u(x, t), fi (x, t) = u(x, t). (11)
i=1 i=1

(0)
By choosing appropriate local equilibrium distribution function fi and using
the Chapman-Enskog multi-scale expansion techniques, we recover the
macroscopic equation (5) or (8) from the time-space lattice BGK model (10).
To this end, we design the amending functions hi (x, t), gi (x, t) as
⎧
⎪ gi (x, t) = γu, i = 1, 2, 3, 4, 5,
⎪
⎪
⎨ h (x, t) = 0,
1
(12)
⎪
⎪ h2 (x, t) = −h3 (x, t) = 2β(x)u,
⎪
⎩
h4 (x, t) = −h5 (x, t) = β(x)u,

where γ is a constant and β(x) is a polynomial with degree less than three. Let
5
 5
 5

(0) (0) (0)
ei fi = 0, e2i fi = λ(x)u, e3i fi = η(x)u. (13)
i=1 i=1 i=1

Then it is easy to verify that

5
 5
 5

eji gi (x, t) = 0, j = 1, 3, gi (x, t) = 5γu, e2i gi (x, t) = 10γu, (14)
i=1 i=1 i=1
5
 5

eji hi (x, t) = 0, j = 0, 2, ei hi (x, t) = 8β(x)u, (15)
i=1 i=1
5

e3i hi (x, t) = 20β(x)u. (16)
i=1
242 Zhiqiang ZHOU, Jingtang MA

Remark 1 It is noted that different from the standard lattice GBK in [6] and
[2], where λ, β, and η are taken to be constants, the proposed lattice GBK in
this paper uses variable λ(x), β(x), and η(x) in order to correctly recover the
macroscopic equation (5) or (8) from the time-space lattice BGK model (10).
The rest of the work is divided into the following four steps.
Step 1 Applying the Taylor expansion to the left-hand side of equation (10)
and retaining terms up to O(ε3 ), we obtain

∂fi ∂ 1 ∂2
−ε2 T + cε (ei fi ) + c2 ε2 2 (e2i fi )
∂t ∂x 2 ∂x
∂ 2 1 ∂3
−cT ε3 (ei fi ) + c3 ε3 3 (e3i fi )
∂x∂t 6 ∂x
1 (0)
= − (fi − fi ) + εhi + ε gi + O(ε4 ).
2
(17)
τ
Step 2 Substituting the Chapman-Enskog scale expansion in time and space

∂ ∂ ∂
= +ε , (18)
∂t ∂t0 ∂t1
(0) (1) (2) (3)
fi = fi + εfi + ε2 fi + ε3 fi + O(ε4 ), (19)
into (17) and collecting terms for ε, ε2 , and ε3 , we get

∂ (0) 1 (1)
c (ei fi ) = − fi + hi , (20)
∂x τ
(0)
∂fi ∂ (1) 1 ∂2 (0) 1 (2)
−T +c (ei fi ) + c2 2 (e2i fi ) = − fi + gi , (21)
∂t0 ∂x 2 ∂x τ
(1) (0)
∂fi ∂f ∂ (2) 1 2 ∂ 2 2 (1)
−T −T i +c (ei fi ) + c (e f )
∂t0 ∂t1 ∂x 2 ∂x2 i i
∂2 (0) 1 3 ∂ 3 3 (0) 1 (3)
−cT (ei fi ) + c 3
(ei fi ) = − fi . (22)
∂x∂t 6 ∂x τ
(k)
Also, we let the non-equilibrium DFs fi , k = 1, 2, 3, satisfy the solvability
conditions
5
(k)
fi = 0, k = 1, 2, 3. (23)
i=1

Step 3 Using (20), (11)–(16), and (23), we obtain


5  5 5
∂ (0) 1  (1) 
0=c ei fi =− fi + hi , (24)
∂x τ
i=1 i=1 i=1
Lattice Boltzmann methods for solving PDEs of exotic option pricing 243

 5
∂ (1)
−τ c (λ(x)u) + 8τ β(x)u = ei fi , (25)
∂x
i=1
5
∂ (1)
−τ c (η(x)u) = e2i fi . (26)
∂x
i=1
Using (21), (11)–(16), and (23)–(26), we obtain
∂u 1  ∂2 ∂
−T + c2 −τ 2
(λ(x)u) + 8τ c (β(x)u) − 5γu = 0, (27)
∂t0 2 ∂x ∂x
 1  ∂2
5

2 (2)
c τ τ− (η(x)u) = ei fi . (28)
2 ∂x2
i=1
Using (11)–(16) and (22)–(28), we obtain
∂u  1  ∂3
−T + c3 τ 2 − τ + (η(x)u) = 0. (29)
∂t1 6 ∂x3
Multiplying (29) by ε and adding to (27), we obtain
∂u 1  ∂2 ∂
−T + c2 −τ 2
(λ(x)u) + 8τ c (β(x)u)
∂t 2 ∂x ∂x
 1  ∂3
−5γu + εc3 τ 2 − τ + (η(x)u) = 0. (30)
6 ∂x3
Step 4 To derive the macroscopic equations, we select
λ(x) = λx2 , β(x) = β1 x + β2 , η(x) = η, (31)
where parameters λ, β1 , β2 , and η are constants. Substituting (31) into (30),
we obtain
∂u 1  ∂u ∂2u 
−T + c2 − τ 2λu + 4λx + λx2 2
∂t 2 ∂x ∂x
 ∂u ∂u   1  ∂3u
3 2
+ 8τ c β1 u + β1 x + β2 − 5γu + εc η τ − τ + = 0. (32)
∂x ∂x 6 ∂x3
Observing the left-hand side of equation (32), the terms u, ∂u ∂u ∂u
∂t , ∂x , x ∂x , and
x2 ∂u
∂x are just the parts we needed. So we are able to recover the macroscopic
PDEs (5) and (8) arising from the Asian options and lookback options,
respectively (see next section for the details.).
3.2 Recovery of macroscopic PDEs
First, we recover the Asian option pricing PDE (5). To make that (32) and
(5) are equivalent, the parameters c, τ, T ∗ , λ, β1 , β2 , and γ must satisfy the
following identities:
T 1  T σ2
2
η = 0, 8τ cβ2 = ∗ , c λ −τ =− ,
T 2 2
244 Zhiqiang ZHOU, Jingtang MA

1  1 
4c2 λ − τ + 8τ cβ1 = rT, 2c2 λ − τ + 8τ cβ1 − 5γ = 0.
2 2
Denote by Δx and Δt the meshsize for space and time and let Δx = cε and
Δt = ε2 T. Then the parameters in time-space lattice BGK model (10) can be
calculated as follows:
σ 2 Δt Δt
λ= , ε2 γ = (r + σ 2 ), (33)
(Δx)2 (2τ − 1) 5
(r + 2σ 2 )Δt Δt
εβ1 = , εβ2 = , η = 0, (34)
8τ Δx 8τ T ∗ Δx
where τ is the relaxation time to be selected. It is used for adjusting the
local equilibrium distribution functions (EDFs) and amending functions. If we
choose a proper value of τ, then the numerical accuracy and stability may be
improved (see, e.g., [12]).
Similarly, recovering (8) for the lookback options from (32) gives the
following identities:
σ 2 Δt Δt 2
λ= , ε2 γ = (σ − 2r), (35)
Δx2 (2τ − 1) 5
(2σ 2 − r)Δt
εβ1 = , εβ2 = 0, η = 0. (36)
8τ Δx
(0)
Now, it is ready to derive the formulas for local EDFs fi . Combining (11),
(13), and (31), and letting
(0) (0)
f4 = f5 = δ(x)u(x, t),
we can derive the local EDFs as follows:
⎧ (0)
⎪
⎪ f = (1 − λx2 + 6δ(x))u(x, t)u =: ρ1 (x)u,
⎪ 1
⎪
⎪
⎪ 1
⎪
⎪ (0)
f2 = (λx2 − 8δ(x))u(x, t) =: ρ2 (x)u,
⎪
⎪
⎨ 2
(0) 1 (37)
⎪
⎪ f3 = (λx2 − 8δ(x))u(x, t) =: ρ3 (x)u,
⎪
⎪ 2
⎪
⎪ (0)
⎪
⎪ f = δ(x)u(x, t) =: ρ4 (x)u,
⎪
⎪
4
⎩ (0)
f5 = δ(x)u(x, t) =: ρ5 (x)u.
To ensue the computation stability of lattice BGK, it is suggested to select δ(x)
such that 0  ρi (x)  1. For example, set
1
δ(x) = λx2 .
16
Then
 5 2 1 2 1
(0) (0) (0) (0) (0)
f1 = 1 − λx u, f2 = f3 = λx u, f4 = f5 = λx2 u. (38)
8 4 16
Lattice Boltzmann methods for solving PDEs of exotic option pricing 245

Finally, the evolution equation (10) can be rewritten as

 1 1
fi (x + Δxei , t − Δt) = 1 − fi (x, t) + ρi (x)u(x, t)
τ τ
+ εe∗i β(x)u(x, t) + ε2 γu(x, t), i = 1, 2, 3, 4, 5, (39)

where e∗ := [0, 2, −2, 1, −1], β(x) and ρi (x) are defined by (31) and (37),
respectively, and the values of involved parameters are given by (33)-(34) for
Asian options and (35)-(36) for lookback options, respectively.
The desired evolution equation (39) is the approximation of PDEs (5) and
(8), respectively, for Asian options and lookback options, with truncation errors
O(ε4 ).

4 Implementation and financial explanation of LBM

4.1 Implementations of the algorithm
In this subsection, we deal with the terminal and boundary conditions and the
implement algorithm. Let Xmax be a sufficiently large value such that
∂
u(Xmax , t) ≈ 0
∂x
for the Asian options or let Xmax = 1 for the lookback options. Let
T∗ Xmax
Δt = , Δx = ,
M N
where M and N are the numbers of nodes for time and space, respectively.
Define uniform time and space meshes:

tk = kΔt, k = M, M − 1, . . . , 0;
xj = jΔx, j = 0, 1, . . . , N,

and denote
5

k
fi,j = fi (xj , tk ), ukj = k
fi,j ≈ u(xj , tk ).
i=1
Now, we describe the methods for dealing with terminal and boundary
conditions (6) and (9) for the Asian options and lookback options, respectively.
The values uM ∗
j of macroscopic variables u(x, T ) in equations (5) and (8) are
set by the terminal conditions (6) or (9) and set
(0)
fi (x, T ∗ ) := fi (x, T ∗ ).

For the left boundary conditions, we take

1
uk0 = u(0, tk ), uk1 = (u(0, tk ) + uk2 ). (40)
2
246 Zhiqiang ZHOU, Jingtang MA

On the right boundary x = Xmax , we let

ukN = ukN −1 = ukN −2 (41)

for the Asian options and take

ukN −1 1 − Δx k
ukN = , ukN −1 = u , (42)
1 − Δx 1 − 2Δx N −2
for the lookback options. The unknown DFs, on boundary nodes xj , are
computed by using the non-equilibrium extrapolation scheme (see [10]), i.e.,
(0) (0)
fi,j = fi,j + fi,j − f , (43)
i,
j

where j represents the adjacent node with j.

Before the ending this subsection, we give the implementation of LBM for
pricing the Asian options and the lookback options. Since the iteration process
does not require to solve algebra system at each time level, the computational
complexity is only O(N ) for each time and O(N ) for memory.
Algorithm of LBM
Step 1 Initialization of model parameters and DFs:

• using formulas (33)-(34) or (35)-(36) to calculate parameters of LBM;

(0)
• calculating DFs: fi (x, T ) by setting to equal fi (x, t) for all nodes at
t=T ; ∗

• computing local EDFs at t = T ∗ from (37).

Step 2 Evolution computation backward in time: for k = M, . . . , 0,

• calculating DFs fi (xj , tk−1 ) for all nodes xj by (39);


• computing solutions uk−1j = 5i=1 fi,j k−1
and deal with boundary by (40)–
(43);
(0)
• computing local EDFs: fi (xj , tk−1 ) at tk−1 by (37).

Step 3 Output solutions u0j and option prices:

• using C(S(0), 0, 0) = S(0)u(K/S(0), 0) for Asian options;

• using C(S(0), 0, 0) = S(0)u(1, 0) for lookback options.

4.2 Financial meaning of lattice Boltzmann models

In this subsection, with Asian options for example, we give the financial
explanation for the lattice BGK models. Similar arguments can be made for
lookback options.
Lattice Boltzmann methods for solving PDEs of exotic option pricing 247

For the Asian options, denote

C(S(t), I(t), t)
R(t) =
S(t)
as the option-stock relative return from time t to expiry time T ∗ . From (4), we
have
C(S(t), I(t), t)  K − I(t) 
T∗
R(t) = = u(x, t) = u ,t . (44)
S(t) S(t)
Define
 1  fi,j
k
1
k
ωi,j = 1−  k + ρi (xj ) + εe∗i β(xj ) + ε2 γ, (45)
τ i fi ,j τ
and it is verified that
   1  fi,j
k
1

k
ωi,j = 1−  k + ρi (xj ) + O(ε) → 1, Δt, Δx → 0.
τ i fi ,j τ
i i

Applying (45), the evolution equation (39) can be rewritten as

k k+1
fi,j = ωi,j−e uk+1 ,
i j−ei
(46)

and thus,  
k+1
ukj = k
fi,j = ωi,j−e uk+1 .
i j−ei
(47)
i i
Denoting stochastic process

K − I(t)
T∗
X(t) =
S(t)
and noting (44), we have

k+1
R(tk )|X=xj = ωi,j−e i
R(tk+1 )|X=xj−ei , (48)
i

which means that the stock return

R(tk ) = u(X(tk ), tk )|X=xj

is an average value of

R(tk+1 ) = u(X(tk+1 ), tk+1 )|X=xj−ei

k+1
with weights ωi,j−e i
. Therefore, we have the following remark.
Remark 2 If we treat ωi,j k as an analog of risk-neutral probabilities of stock

price S(t), then the lattice GBK model (10) can be interpreted by a special
form of binomial option pricing model ([18]).
248 Zhiqiang ZHOU, Jingtang MA

5 Numerical experiments

In this section, we test the lattice Boltzmann model proposed in this paper. In
comparison with the results in existing literature, the stability, efficiency, and
computation cost of the lattice Boltzmann models are verified. In the following
examples, we take Xmax = 2 for Asian options, T ∗ = 1.0, and the EDFs are
taken as (38).
Example 1 This example is used to verify the stability of LBM proposed by
this paper. The values of the parameters are set

σ = 0.15, r = 0.05, Δt = 0.0002, Δx = 0.005, τ = 1.05.

As shown in Fig. 1, the computational solutions do not have vibration or

blow-up and agree with the results listed in [4,9]. However, the algorithm
of LBM is not stable for Asian options if Δt = 1/1400. Hence, the numerical
scheme is conditionally stable, and experiments show that the factors
influencing stability are the time-space mesh ratio Rm = Δt/Δx2 and the

Fig. 1 Numerical solutions of PDE problem (5)-(6) Asian options (a) (simulated time t is
taken as 0.1667, 0.2500, 0.3333, 0.4167, 0.5, 0.5833, 0.6667, 0.75, 0.8333, 0.9167, and 1,
respectively) and (8)-(9) lookback options by LBM (b) (simulated time t is taken as 0.0333,
0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0, respectively)
Lattice Boltzmann methods for solving PDEs of exotic option pricing 249

relaxation time τ. The stability of the equation requires that τ > 0.5 (see [2]),
and to get higher convergence rate, it is suggested that τ < 1.5.
Table 1 lists the threshold values of time-space mesh ratio Rm for fixed
τ = 1.05. Figs. 2 and 3 plot the stable and unstable areas for Asian options and
Table 1 Threshold value of stable time-space mesh ratio Rm with τ = 1.05
Asian options lookback options
N M Rm M Rm
100 280 8.92 248 40.32
140 670 7.31 503 38.96
180 1160 6.98 846 38.29
220 1780 6.79 1280 37.81
260 2540 6.65 1803 37.49
300 3420 6.57 2417 37.23
340 4450 6.49 3122 37.02
380 5575 6.47 3915 36.88

Fig. 2 Stable areas of LBM for Asian options (a)

and lookback options (b) with fixed τ = 1.05
250 Zhiqiang ZHOU, Jingtang MA

Fig. 3 Stable areas of LBM for Asian options (a)

and lookback options (b) with fixed N = 100

lookback options with fixed τ = 1.05 and fixed N = 100, respectively. Fig. 2
shows that the Δt decreases greatly as the Δx decreases, and Table 1 tells us
that the value of threshold time-space mesh ratio Rm is close to a constant.
From Fig. 3, we see that there is a positive linear correlation between Rm and
τ. Moreover, the LBM has better stability for lookback options than that for
Asian options, which is observed from that the time-space mesh ratio Rm ≈ 7
for computing lookback options is smaller than Rm ≈ 37 for Asian options.
The main reason is that the convection term in Asian option pricing PDEs
dominates the diffusion term, compared to lookback option pricing PDEs.
Example 2 In comparison with the option prices derived by Rogers and Shi
[17] and Goldman et al. [9], we list the computational results by LBM in Tables
2 and 3.
Table 2 presents numerical results of prices for Asian options at time t = 0
with S(0) = 100. The numbers in column 1 give the strike prices, columns 2
Lattice Boltzmann methods for solving PDEs of exotic option pricing 251

Table 2 Prices for Asian options at t = 0 for Example 2 (r = 0.05, τ = 0.9)

K LB UB PDE PDE LBM
($) (Rogers) (Rogers) (Rogers) (Boyle) (this paper)
σ = 0.30, M = 5000, N = 200
90 13.952 14.1610 13.951 13.964 13.9753
100 7.9440 8.1530 7.9440 7.9595 7.96833
110 4.0700 4.2790 4.0740 4.0856 4.09796
σ = 0.20, M = 5000, N = 200
90 12.595 12.687 12.589 12.607 12.6014
100 5.7620 5.8540 5.7600 5.7835 5.77767
110 1.9890 2.0800 1.9960 2.0087 2.03079
σ = 0.10, M = 9000, N = 500
90 11.951 11.973 11.942 11.955 11.962
100 3.6410 3.6630 3.6240 3.6566 3.6191
110 0.3310 0.3530 0.3590 0.3455 0.3426

Table 3 Prices for lookback options at t = 0 with r = 0.05 for Example 2

(τ = 1.25, M = 5000, N = 500)
σ = 0.15 σ = 0.25
S(0) Goldman Boyle LBM Goldman Boyle LBM
80 8.0085 9.5903 8.0604 14.978 13.8243 15.0330
90 9.0096 10.789 9.0679 16.851 15.5524 16.9121
100 10.010 11.987 10.076 18.723 17.2804 18.7912
110 11.011 13.186 11.083 20.595 19.0085 20.6704
120 12.012 14.385 12.091 22.467 20.7365 22.5495

and 3, respectively, give the lower and upper bounds of the option prices
obtained by Rogers and Shi [17], and columns 4 to 6, respectively, give the
results from Rogers-Shi PDE approach [17], Boyle’s PDE approach [4], and
the LBM approach proposed in this paper. As observed, the option values
obtained by the LBM approach are within the lower and upper bounds
provided by Rogers and Shi [17]. The relative change in option prices increases
as the volatility σ increases, and this is more so for smaller strike prices.
Table 3 lists the values of the lookback options obtained by Goldman et al.
[9], Boyle et al. ([4], prices under regime switching) and by the LBM approach
in this paper. All the numerics in the last column agree with those of Goldman
et al. [9], which verify the accuracy of LBM. Similar to Asian option, the value
of the lookback options increases as the volatility σ increases.
From Table 4, we see that our results (without brackets) under different
expired time T ∗ are close to the prices (with brackets) given by Boyle and
Draviam [4]. The relative change in option prices increases as the maturity
period increases.
Example 3 In this example, we compare the CPU time between LBM and
FDM. The CPU time growth rate is defined by
252 Zhiqiang ZHOU, Jingtang MA

log(CPU time on finer mesh/CPU time on coarser mesh)

RCPU = .
log(number of finer mesh points/number of coarser mesh points)

Table 4 Prices for exotic options with different expired time T ∗ for Example 2
Asian option (K = 90) lookback option (S(0) = 80)
T∗ σ1 = 0.15 σ2 = 0.25 σ1 = 0.15 σ2 = 0.25
0.25 10.515 (10.512) 10.639 (10.699) 4.422 (4.395) 7.6626 (7.7613)
0.50 11.041 (11.051) 11.551 (11.569) 5.961 (5.985) 10.732 (10.825)
0.75 11.597 (11.607) 12.425 (12.423) 6.920 (7.115) 12.981 (13.106)
1.00 12.152 (12.161) 13.239 (13.228) 8.064 (8.008) 15.033 (14.978)
1.25 12.697 (12.705) 14.001 (13.984) 8.800 (8.750) 16.282 (16.587)

All the numerical experiments are run in MATLAB 7.9.0 (R2009b) on a

PC with the configuration: Intel(R) Core(TM) i5-4200 CPU @ 2.5 GHz and
4.00 GB RAM. The FDM is designed as the algorithm proposed by Boyle et al.
[4] and the conjugate gradient normal residual (CGNR) method is applied to
solve the linear algebraic system. The FDM of Boyle et al. [4] is constructed
by the following idea: integrate equation (5) over a cell (xj− 1 , xj+ 1 ), where
2 2

xj−1 + xj xj + xj+1
xj− 1 := , xj+ 1 := ,
2 2 2 2
for all available index j, then use quadrature rules to approximate the integrals,
and finally obtain the numerical scheme using the following exponential
interpolation:
e−P xj+1 − e−P x e−P x − e−P xj
u(x, t) ≈ u(xj , t) + u(xj+1 , t) ,
e−P xj+1 − e−P xj e−P xj+1 − e−P xj
where
1
−rxj+ 1 − T
P = 1
2
.
2 σ 2 x2j+ 1
2

From Table 5, we see that the time growth rate RCPU is larger than 2 for
FDM and only 1 for LBM. So, the CPU time needed for LBM increases much
more slowly than that of the FDM.

Table 5 Comparison of CPU time between LBM and FDM for Asian options
(σ = 0.15, r = 0.05, τ = 1.2, M = 8000)
LBM FDM
N CPU time(s) RCPU CPU time(s) RCPU
100 11.9028 0.9886 125.6588 1.5256
200 23.6185 0.9446 361.7819 2.0495
300 35.3498 1.0180 830.5181 2.0718
400 47.3775 1.0127 1507.397 2.1040
500 59.3895 — 2410.449 —
Lattice Boltzmann methods for solving PDEs of exotic option pricing 253

6 Conclusions

In this paper, we apply lattice BGK model with two amending functions to
solve the PDEs arising from Asian option and lookback option pricing. By
the Taylor approximation and the Chapman-Enskog multi-scale expansion,
both PDEs of Asian options and lookback options are precisely recovered from
the lattice BGK model. The most important reasons to introduce lattice
Boltzmann models into financial area are the lower computation workloads
and the good financial meanings. In the future, we will consider the nonlinear
PDEs of option pricing such as the options with transaction costs, stochastic
volatility, and stochastic interest rates. We believe that more and more
advantages and challenges will be found for lattice Boltzmann models for
financial problems.

Acknowledgements The authors were grateful to the anonymous referees for their
valuable suggestions that led to a greatly improved paper. This work was supported by the
National Natural Science Foundation of China (Grant No. 11171274) and the Program for
New Century Excellent Talents in University (Grant No. NCET-12-0922).

References
1. Al-Zoubi A, Brenner G. Comparative study of thermal flows with different finite
volume and lattice Boltzmann schemes. Int J Mod Phys C, 2004, 15: 307–319
2. Bhantnagar J, Gross E P, Krook M K. A model for collision processes in gas I: Small
amplitude processes in charged and neutral one-componet systems. Phys Rev, 1954,
94: 511–525
3. Black F, Scholes M. The pricing of options and corporate liabilities. J Polit Econom,
1973, 81: 637–654
4. Boyle P, Draviam T. Pricing exotic options under regime switching. Insurance: Math
Econom, 2007, 40: 267–282
5. Broadie M, Detemple J. Option pricing: valuation models and applications.
Management Sci, 2004, 50: 1145–1177
6. Chen L, Ma C F. A lattice Boltzmann model with an amending function for simulating
nonlinear partial differential equations. Chin Phys B, 2010, 19: 1–8
7. Dubois F, Lelièvre T. Efficient pricing of Asian options by the PDE approach.
J Comput Finan, 2005, 8: 55–64
8. Geman H, Yor M. Bessel processes Asian options and perpetuities. Math Finan, 1993,
3: 349–375
9. Goldman M B, Sosin H B, Gatto M A. Path-dependent options buy at the low, sell at
the high. J Finan, 1979, 34: 1111–1127
10. Guo Z L, Zheng C G, Shi B C. Non-equilibrium extrapolation method for velocity and
pressure boundary conditions in the lattice Boltzmann method. Chin Phys, 2002, 11:
366–374
11. Ingersoll J. Theory of Financial Decision Making. Totowa/New Jersey: Roman &
Littlefield, 1987
12. Lai H L, Ma C F. A higher order lattice BGK model for simulating some nonlinear
partial differential equations. Sci China Ser G, 2009, 52: 1053–1061
13. Lai H L, Ma C F. Lattice Boltzmann method for the generalized Kuramoto-Sivashinsky
equation. Phys A, 2009, 388: 1405–1412
254 Zhiqiang ZHOU, Jingtang MA

14. Lai H L, Ma C F. Lattice Boltzmann model for generalized nonlinear wave equations.
Phys Rev E, 2011, 84: 046708
15. Lai H L, Ma C F. A new lattice Boltzmann model for solving the coupled viscous
Burgers’ equation. Phys A, 2014, 395: 445–457
16. Merton R. Theory of rational option pricing. Bell J Econom Management Sci, 1973, 4:
141–183
17. Rogers L C G, Shi Z. The value of an Asian option. J Appl Probab, 1995, 32: 1077–1088
18. Shreve S E. Stochastic Calculus for Finance: the Binomial Asset Pricing Model. New
York: Springer, 2003
19. Succi S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Oxford:
Oxford University Press, 2001
20. Sukop M C, Thorne D T. Lattice Boltzmann Modeling: An Introduction for
Geoscientists and Engineers. Heidelberg/Berlin/New York: Springer, 2006
21. Vecer J. A new PDE approach for pricing arithmetic average Asian options. J Comput
Finan, 2001, 4: 105–113
22. Wilmott P, Dewynne J, Howison S. Option Pricing: Mathematical Models and
Computation. Oxford: Oxford Financial Press, 1997