In this paper, we present a detailed comparison of the performance of the numerical solvers such as the biconjugate gradient stabilized, operator splitting, and multigrid methods for solving the twodimensional BlackScholes equation. The equation is discretized by the finite difference method. The computational results demonstrate that the operator splitting method is fastest among these solvers with the same level of accuracy.
1. INTRODUCTION
Black and Scholes
[2]
derived the BlackScholes (BS) partial differential equation for the valuation of a European option under the noarbitrage assumption. Various types of exotic options are popular in the market. Finding the analytic closedform solution of the BS equation is not easy. Therefore, it is necessary to apply numerical methods to obtain the values of exotic options. The finite difference methods (FDM), which converts the differential equations into a system of difference equations, are very popular to approximate the solution of the BS equations
[5]
. There have been many numerical methods and among them, we focus on biconjugate gradient stabilized (BiCGSTAB) [21], operator splitting (OS)
[11]
and multigrid (MG)
[16]
methods in this paper.
BiCGSTAB method was introduced by H.A. van der Vorst
[21]
, which is similar to conjugate gradient stabilized (CGS) method with favorable stability properties. As a iterative type method, BiCGSTAB method is appropriate to solve the problem when the coe±cient matrix of problem is large and sparse. MG method introduced by R.P. Fedorenko
[6
,
7]
is numerical algorithm using a hierarchy of discretizations. By employing different mesh size, a multigrid algorithms are combined by smoothers and coarsegrid correction procedures. For this reason, this method provides rapid convergence rates than the standard iterative techniques such as the Jacobi and Gauss{Seidel schemes. There have been applications in option pricing by many researchers
[12
,
15
,
16
,
17]
. On option pricing, OS method was proposed by S. Ikonen and J. Toivanen
[11]
. This method is by decoupling a complex equation in various simpler equations and solving the simpler equation with discretization. Since then, many researchers
[4
,
5
,
12]
have applied OS method to the BS equation.
For different types of problems, different system solvers gain advantages over the other methods, see
[19]
. To show the performance of the finite difference schemes for the twodimensional problems, we compare the wellknown solvers, BiCGSTAB, OSM, and MG methods, for the twodimensional BS equations. There also have been other system solvers, such as alternating direction method (ADI)
[3]
and generalized minimal residual algorithm (GMRES)
[13
,
18]
, however we omit the comparison in this work since GMRES and ADI methods are similar to BiCGSTAB and OS methods, respectively. The outline of the paper is as follows. In Section 2, we first set up the problem to price stock options. In Section 3, we describe the general setting of numerical strategies and explain different solvers of linear system. In Section 4, we show the comparison of the numerical experiments between the solvers. The conclusions are drawn in Section 5.
2. BLACKSCHOLES EQUATIONS
Let
s_{i}
(
t
), i = 1, 2, ... ,
n
, be the price of ith asset at time
t
and be the unique solution to a geometric Brownian motion with a constant volatility
σ_{i}
> 0,
i
= 1, 2, ... ,
n
. Let
S
= (
s
_{1}
,
s
_{2}
, ... ,
s
_{n}
) be the vector of asset prices and
ρ_{ij}
,
i
,
j
= 1, ... ,
n
, be the correlation coeffcients between Brownian motions. We assume that the interest rate is constant,
V
(
S
,
t
) is the value of a European option that underlies assets 1, ... ,
n
,
T
is the expiration date, and
A
(
S
) is the payoff function. By following the `noarbitrage' argument for the BS equation, a partial differential equation for
V
is derived to be
In this paper, we use the original BS model with two underlying assets to keep this presentation simple. However, we can easily extend the current method for more than two underlying assets
[1]
. Let us consider the computational domain Ω = (0,
L
) × (0,
M
) and
x
=
s
_{1}
and
y
=
s
_{2}
. Let us first convert the given backward equation (2.1) to the following forward equation by a change of variable 𝜏 =
T

t
,
u
(
x
,
y
, 𝜏) =
V
(
s
_{1}
,
s
_{2}
,
T
 𝜏):
We use the following linear boundary conditions on all boundaries,
3. NUMERICAL METHOD
In this paper, we discretize the partial derivatives in Eq. (2.2) using finite difference methods that have been used in option pricing.
3.1. Discretization
Let us first discretize the given computational domain Ω = (0,
L
) × (0,
M
) as a uniform grid with a space step
h
=
L
/
N_{x}
=
M
/
N^{y}
and a time step ∆τ =
T/N_{τ}
. Here,
N_{x}
and
N_{y}
are the number of grid points, and
N_{τ}
is the total number of time steps. Let the numerical approximations of the solution be
≈
u
((i0.5)
h
, (
j
 0.5)
h
,
n
∆
τ
), where
i
= 1, ... ,
N_{x}
,
j
= 1, ... ,
N_{y}
and
n
= 0, 1, ... ,
N_{τ}
. We use 𝜕
u
/𝜕
x
≈ (
u
_{i+1,j}

u
_{i,j}
)/
h
, 𝜕
^{2}
u
/𝜕
x
^{2}
≈ (
u
_{i1, j}
 2
u_{ij}
+
u
_{i+1, j}
)/
h
^{2}
, 𝜕
^{2}
u
/𝜕
x
𝜕
y
≈ (
u
_{i+1, j+1}
+
u_{ij}

u
_{i,j+1}

u
_{i+1, j}
)/
h
^{2}
, and 𝜕
u
/𝜕
τ
≈ (
u
^{n+1}

u^{n}
)/∆
τ
.
3.2. BiCGSTAB
The biconjugate gradient stabilized method (BiCGSTAB) was developed to solve nonsymmetric linear systems
[21]
. We solve Eq. (2.2) by BiCGSTAB method. We write Eq. (2.2) in a discretizad form:
where
Next, to renumber the multiindexed data
u_{ij}
as the singleindexed data
U_{l}
, we denote by
U_{l}
=
U
_{Nx(j1)+i}
=
u_{ij}
, where
l
= 1, ... ,
N_{x}
×
N_{y}
,
i
= 1, ... ,
N_{x}
, and
j
= 1, ... ,
N_{y}
. Consequently, we get the following system
where
U
^{n+1}
=
,
b
^{n}
=
, and matrix A is composed of coe±cients of
U
. To solve the linear system (3.1), BiCGSTAB starts with an intial guess
U
^{0}
and proceeds as follows:
BiCGSTAB cycle
Define the maximum number of iteration
ITER
and the error tolerance
TOL
Set
r
^{0}
=
b

A
U
^{0}
,
=
r
^{0}
,
𝜌
^{0}
= α = 𝜔
^{0}
= 1,
v
^{0}
=
p
^{0}
= 0,
k
= 1
While (
k
≤
ITER
& 
r
^{k}

_{2}
>
TOL
)
End While
3.3. Operator splitting method
The basic idea of operator splitting method is to split the spatial operator into onedimensional operators and then fractional time steps are performed with these simpler operators. The operator splitting method computes the solutions in two time steps:
where the discrete difference operators
and
are defined by
In OS method, we first solve (

)/∆τ =
, and then we solve (

)/∆τ =
.
3.4. Multigrid method
Multigrid methods belong to the class of fastest iterations, because their convergence rate is independent of the step size
h
, see
[8]
. We define a discrete domain by Ω
_{k}
= {(
h
(
i
 0.5),
h
(
j
 0.5))1≤
i
,
j
≤ 2
^{k+1}
}. Ω
_{k1}
is coarser than Ω
_{k}
by factor 2. The multigrid solution of the discrete BS equation
makes use of a hierarchy of meshes created by successively coarsening the original mesh, see
Fig. 1
.
A sequence of coarse grids starting with h.
We use a multigrid cycle to solve the discrete system at the implicit time level. A pointwise GaussSeidel relaxation scheme is used as the smoother in the multigrid method. We first rewrite the above equation (3.2) by
L
(
) =
for each (
i, j
) ∈ Ω
_{k}
, where
L
(
) =
 ∆τ𝓛
_{BS}
Given the number 𝜈
_{1}
and 𝜈
_{2}
of pre and post smoothing relaxation sweeps, an iteration step for the multigrid method using the Vcycle is formally written as follows
[20]
. We use a notation
as a numerical solution on the discrete domain Ω
_{k}
at time
t
=
n
∆τ. Given
, we want to find
solution which satisfies equation (3.2). At the very beginning of the multigrid cycle the solution from the previous time step is used to provide an initial guess for the multigrid procedure. First, let
The algorithm of the multigrid method for solving the discrete BS equation (3.2) is following:
Multigrid cycle
Step 1) Presmoothing:
perform 𝜈
_{1}
GaussSeidel relaxation steps.
Step 2) Coarse grid correction

⋅ Compute the residual on Ωk:

⋅ Restriction to Ωk1:,

⋅ Compute an approximation soultion on Ωk1:

⋅ Solve the equation (3.4):


⋅ Interpolate the correction:

⋅ Compute the corrected approximation on Ωk:
Step 3) Postsmoothing:
=
SMOOTH^{ν2}
4. Computational Results
In this section, we compare the performance of the numerical methods (BiCGSTAB, OS, and MG) using CPU times. Each method is implemented using MATLAB
[14]
. We consider three types of twoasset cashornothing options. The cashornothing options are useful building blocks for constructing more complex exotic option products and they are widely traded in the real world financial market.
Case 1: A two asset cashornothing call pays out a fixed cash amount K if asset one, x, is above the strike X_{1} and asset two, y, is above strike X_{2} at expiration. The payoff is given by
Case 2. and Case 3.:
Figures 2(a)
,
(b)
, and
(c)
show the payoff function A(
x, y
) for
Case 1
,
Case 2
, and
Case 3
, respectively. The closedform solutions
[9]
are
Case 1
:
u
(
x
,
y
,
T
) =
K_{e}^{rT}M
(
α
,
β
; 𝜌),
Case 2
:
u
(
x
,
y
,
T
) =
K_{e}^{rT}M
(
α
,
β
; 𝜌),
Case 3
:
u
(
x
,
y
,
T
) =
K_{e}^{rT}M
(
α
,
β
), 𝜌), where
[10]
. Let 𝜌 be the correlation between the two variables, then
Payoff functions of (a)Case 1, (b)Case 2, and (c)Case 3, respectively.
We computed the numerical solution on uniform grids,
h
= 300/2
^{n}
for
n
= 5, 6, 7, and 8 on the computational domain Ω = [0, 300] × [0, 300]. For each case, we ran the calculation to time
T
= 1 with a uniform time step ∆τ = 0.01 with a given strike price of
X
_{1}
= 100,
X
_{2}
= 100 and cash amount
K
= 1. The volatilities are 𝜎
_{1}
= 0.3, 𝜎
_{2}
= 0.3 with a correlation 𝜌 = 0.5, and the riskless interest rate
r
= 0.03.
Figure 3
shows the numerical solution at
T
= 1 case by case. We let e be the matrix with components
e_{ij}
=
u
(
x_{i}
,
y_{j}
)
U_{ij}
and compute its discrete
l
_{2}
norm of the error, 
e

_{2}
:
Numerical solutions at time T = 1 of (a) Case 1, (b) Case 2, and (c) Case 3, respectively.
We test the numerical experiments of different case with three solvers, BiCGSTAB, OSM and MG. To make a fair comparison of these solvers, we match the accuracy of these solvers by changing iteration parameters.
In this figures, the solid line with triangles, the dashdot line with squares, and the dashed line with stars express OSM, BICGSTAB, and MG, respectively. Next, let us check the CPU times to compare efficiency of these solvers.
Table 1
also shows the CPU times and
l
_{2}
error with each method. We can confirm that OS method has a linear CPU time cost as the spatial domain is doubled in each direction.
Table 2
and
Table 3
also show the CPU times and
l
_{2}
error with
Case 2
and
Case 3
. And the corresponding results are plotted in
Figs. 4(b)
and
(c)
, respectively. From all these results, we can confirm that OS method is faster than other methods under the same accuracy.
(Case 1)Comparison of l2error and CPU time.
(Case 1) Comparison of l_{2} error and CPU time.
(Case 2)Comparison of l2error and CPU time.
(Case 2) Comparison of l_{2} error and CPU time.
(Case 3)Comparison of l2error and CPU time.
(Case 3) Comparison of l_{2} error and CPU time.
CPU times of (a) Case 1, (b) Case 2, and (c) Case 3, respectively.
5. CONCLUSION
The main purpose of this paper is to present the performance comparison of finite difference schemes of the BS equation for stock option pricing. The large linear system, derived from the discrete BS equation, was solved by biconjugate gradient stabilized, operator splitting, and multigrid methods. The performance of these methods was compared for two asset option problems based on twodimensional BS equations. The numerical results indicated that although BiCGSTAB and multigrid solvers are accurate, they need a lot of computational times. On the other hand, operator splitting is faster than the other two methods under the same accuracy.
Acknowledgements
The corresponding author (J.S. Kim) was supported by the National Institute for Mathematical Sciences(NIMS) grant funded by the Korea government(No. A21301). And this work is based on the first author’s Ph. D. thesis[12].
Amstera P.
,
Averbuj C.
,
de Napoli P.
,
Mariani M.
2010
A parabolic problem arising in financial mathematics
Nonlinear Anal. Real World Appl.
11
759 
763
DOI : 10.1016/j.nonrwa.2009.01.019
Black F.
,
Scholes M.
1973
The pricing of options and corporate liabilities
J. Polit. Econ.
81
637 
659
DOI : 10.1086/260062
Chin R.
,
Manteuffel T.
,
de Pillis J.
1984
ADI as a preconditioning for solving the convectiondiffusion equation
SIAM J. Sci. Comput.
5
(2)
281 
299
DOI : 10.1137/0905020
Daoud Y.
,
Öziş T.
2011
The operator splitting method for BlackScholes equation
Appl. Math.
2
(6)
771 
778
DOI : 10.4236/am.2011.26103
Duffy D.
2006
Finite Difference Methods in Financial Engineering, A Partial Differential Equation Approach
John Wiley and Sons
New York, USA
Fedorenko R.P.
1962
A relaxation method for solving elliptic difference equations
USSR Computational Math. and Math. Phys.
1
(4)
1092 
1096
DOI : 10.1016/00415553(62)900319
Fedorenko R.P.
1964
The speed of convergence of one iterative process
USSR Computational Math. and Math. Phys.
4
(3)
227 
235
Hackbusch W.
1994
Iterative Solution of Large Linear Systems of Equations
Springer
New York, USA
Heynen R.
,
Kat H.
1996
Brick by Brick
Risk Magazine
9
(6)
28 
31
Haug E.
2007
The Complete Guide to Option Pricing Formulas
McGrawHill
New York, USA
Jeong D.
2012
Mathematical model and numerical simulation in computational finance. Ph.D. Thesis
Department of Mathematics, Korea University
Korea
Lotstedt P.
,
Persson J.
,
von Sydow L.
,
Tysk J.
2007
Spacetime adaptive finite difference method for European multiasset options
Comput. Math. Appl.
53
(8)
1159 
1180
DOI : 10.1016/j.camwa.2006.09.014
1990
Users Guide: Natick
Mathworks Inc
MA, USA
Oosterlee C.W.
2003
On multigrid for linear complementarity problems with application to Americanstyle options
Electron. Trans. Numer. Anal.
15
(27)
165 
185
Ramage A.
,
von Sydow L.
2011
A multigrid preconditioner for an adaptive BlackScholes solver
BIT Numer. Math.
51
(1)
217 
233
DOI : 10.1007/s1054301103166
Reisinger C.
,
Wittum G.
2004
On multigrid for anisotropic equations and variational inequalities “pricing multidimensional european and american options”
Comput. Vis. Sci.
7
(34)
189 
197
Saad Y.
,
Schultz M.
1986
GMRES, a generalized minimal residual algorithm for solving nonsymmetric linear systems
SIAM J. Sci. & Stat. Comput.
7
(3)
856 
869
DOI : 10.1137/0907058
Trottenberg U.
,
Oosterlee C.
,
Schüller A.
2001
Multigrid
Academic press
London, England
van der Vorst H.
1992
BiCGSTAB, A fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems
SIAM J. Sci. & Stat. Comput.
13
(2)
631 
644
DOI : 10.1137/0913035