Advanced
COMPARISON OF NUMERICAL METHODS (BI-CGSTAB, OS, MG) FOR THE 2D BLACK-SCHOLES EQUATION
COMPARISON OF NUMERICAL METHODS (BI-CGSTAB, OS, MG) FOR THE 2D BLACK-SCHOLES EQUATION
Journal of the Korean Society of Mathematical Education Series B The Pure and Applied Mathematics. 2014. Apr, 21(2): 129-139
Copyright © 2014, Korean Society of Mathematical Education
  • Received : March 01, 2014
  • Accepted : April 25, 2014
  • Published : April 27, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
DARAE JEONG
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:tinayoyo@korea.ac.kr
SUNGKI KIM
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:sgkim82@korea.ac.kr
YONGHO CHOI
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:poohyongho@korea.ac.kr
HYEONGSEOK HWANG
DEPARTMENT OF FINANCIAL ENGINEERING, KOREA UNIVERSITY, SEOUL 136-071, REPUBLIC OF KOREAEmail address:hhs288@naver.com
JUNSEOK KIM
DEPARTMENT OF MATHEMATICS, KOREA UNIVERSITY, SEOUL 136-713, REPUBLIC OF KOREAEmail address:cfdkim@korea.ac.kr

Abstract
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 two-dimensional Black-Scholes 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.
Keywords
1. INTRODUCTION
Black and Scholes [2] derived the Black-Scholes (BS) partial differential equation for the valuation of a European option under the no-arbitrage assumption. Various types of exotic options are popular in the market. Finding the analytic closed-form 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 (Bi-CGSTAB) [21], operator splitting (OS) [11] and multigrid (MG) [16] methods in this paper.
Bi-CGSTAB 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, Bi-CGSTAB 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 coarse-grid 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 two-dimensional problems, we compare the well-known solvers, Bi-CGSTAB, OSM, and MG methods, for the two-dimensional 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 Bi-CGSTAB 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. BLACK-SCHOLES EQUATIONS
Let si ( t ), i = 1, 2, ... , n , be the price of i-th 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 `no-arbitrage' argument for the BS equation, a partial differential equation for V is derived to be
PPT Slide
Lager Image
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 - 𝜏):
PPT Slide
Lager Image
We use the following linear boundary conditions on all boundaries,
PPT Slide
Lager Image
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 / Nx = M / Ny and a time step ∆τ = T/Nτ . Here, Nx and Ny are the number of grid points, and Nτ is the total number of time steps. Let the numerical approximations of the solution be
PPT Slide
Lager Image
u ((i-0.5) h , ( j - 0.5) h , n τ ), where i = 1, ... , Nx , j = 1, ... , Ny and n = 0, 1, ... , Nτ . We use 𝜕 u /𝜕 x ≈ ( u i+1,j - u i,j )/ h , 𝜕 2 u /𝜕 x 2 ≈ ( u i-1, j - 2 uij + u i+1, j )/ h 2 , 𝜕 2 u /𝜕 x 𝜕 y ≈ ( u i+1, j+1 + uij - u i,j+1 - u i+1, j )/ h 2 , and 𝜕 u /𝜕 τ ≈ ( u n+1 - un )/∆ τ .
3.2. Bi-CGSTAB The bi-conjugate gradient stabilized method (Bi-CGSTAB) 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:
PPT Slide
Lager Image
where
PPT Slide
Lager Image
Next, to renumber the multi-indexed data uij as the single-indexed data Ul , we denote by Ul = U Nx(j-1)+i = uij , where l = 1, ... , Nx × Ny , i = 1, ... , Nx , and j = 1, ... , Ny . Consequently, we get the following system
PPT Slide
Lager Image
where U n+1 =
PPT Slide
Lager Image
, b n =
PPT Slide
Lager Image
, and matrix A is composed of coe±cients of U . To solve the linear system (3.1), Bi-CGSTAB starts with an intial guess U 0 and proceeds as follows:
Bi-CGSTAB cycle
Define the maximum number of iteration ITER and the error tolerance TOL Set r 0 = b - A U 0 ,
PPT Slide
Lager Image
= r 0 , 𝜌 0 = α = 𝜔 0 = 1, v 0 = p 0 = 0, k = 1
While ( k ITER & || r k || 2 > TOL )
PPT Slide
Lager Image
End While
3.3. Operator splitting method The basic idea of operator splitting method is to split the spatial operator into one-dimensional operators and then fractional time steps are performed with these simpler operators. The operator splitting method computes the solutions in two time steps:
PPT Slide
Lager Image
where the discrete difference operators
PPT Slide
Lager Image
and
PPT Slide
Lager Image
are defined by
PPT Slide
Lager Image
In OS method, we first solve (
PPT Slide
Lager Image
-
PPT Slide
Lager Image
)/∆τ =
PPT Slide
Lager Image
PPT Slide
Lager Image
, and then we solve (
PPT Slide
Lager Image
-
PPT Slide
Lager Image
)/∆τ =
PPT Slide
Lager Image
PPT Slide
Lager Image
.
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 }. ­Ω k-1 is coarser than ­Ω k by factor 2. The multigrid solution of the discrete BS equation
PPT Slide
Lager Image
makes use of a hierarchy of meshes created by successively coarsening the original mesh, see Fig. 1 .
PPT Slide
Lager Image
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 Gauss-Seidel relaxation scheme is used as the smoother in the multigrid method. We first rewrite the above equation (3.2) by L (
PPT Slide
Lager Image
) =
PPT Slide
Lager Image
for each ( i, j ) ∈ Ω­ k , where L (
PPT Slide
Lager Image
) =
PPT Slide
Lager Image
- ∆τ𝓛 BS
PPT Slide
Lager Image
Given the number 𝜈 1 and 𝜈 2 of pre- and post- smoothing relaxation sweeps, an iteration step for the multigrid method using the V-cycle is formally written as follows [20] . We use a notation
PPT Slide
Lager Image
as a numerical solution on the discrete domain Ω­ k at time t = n ∆τ. Given
PPT Slide
Lager Image
, we want to find
PPT Slide
Lager Image
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
PPT Slide
Lager Image
The algorithm of the multigrid method for solving the discrete BS equation (3.2) is following:
Multigrid cycle
PPT Slide
Lager Image
Step 1) Presmoothing: perform 𝜈 1 Gauss-Seidel relaxation steps.
PPT Slide
Lager Image
Step 2) Coarse grid correction
  • ⋅ Compute the residual on Ωk:
  • ⋅ Restriction to Ωk-1:,
  • ⋅ Compute an approximation soultion on Ωk-1:
PPT Slide
Lager Image
  • ⋅ Solve the equation (3.4):
  • ⋅ Interpolate the correction:
  • ⋅ Compute the corrected approximation on Ωk:
Step 3) Postsmoothing:
PPT Slide
Lager Image
= SMOOTHν2
PPT Slide
Lager Image
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 two-asset cash-or-nothing options. The cash-or-nothing 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 cash-or-nothing call pays out a fixed cash amount K if asset one, x, is above the strike X1 and asset two, y, is above strike X2 at expiration. The payoff is given by
PPT Slide
Lager Image
Case 2. and Case 3.:
PPT Slide
Lager Image
Figures 2(a) , (b) , and (c) show the payoff function A( x, y ) for Case 1 , Case 2 , and Case 3 , respectively. The closed-form solutions [9] are Case 1 : u ( x , y , T ) = Ke-rTM ( α , β ; 𝜌), Case 2 : u ( x , y , T ) = Ke-rTM ( , ; 𝜌), Case 3 : u ( x , y , T ) = Ke-rTM ( , β ), -𝜌), where
PPT Slide
Lager Image
[10] . Let 𝜌 be the correlation between the two variables, then
PPT Slide
Lager Image
PPT Slide
Lager Image
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 eij = u ( xi , yj )- Uij and compute its discrete l 2 -norm of the error, || e || 2 :
PPT Slide
Lager Image
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, Bi-CGSTAB, 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 dash-dot line with squares, and the dashed line with stars express OSM, BI-CGSTAB, 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.
PPT Slide
Lager Image
(Case 1) Comparison of l2 error and CPU time.
(Case 2)Comparison of l2error and CPU time.
PPT Slide
Lager Image
(Case 2) Comparison of l2 error and CPU time.
(Case 3)Comparison of l2error and CPU time.
PPT Slide
Lager Image
(Case 3) Comparison of l2 error and CPU time.
PPT Slide
Lager Image
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 two-dimensional BS equations. The numerical results indicated that although Bi-CGSTAB 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].
References
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 convection-diffusion equation SIAM J. Sci. Comput. 5 (2) 281 - 299    DOI : 10.1137/0905020
Daoud Y. , Öziş T. 2011 The operator splitting method for Black-Scholes 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/0041-5553(62)90031-9
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 McGraw-Hill New York, USA
Ikonen S. , Toivanen J. 2004 Operator splitting methods for American option pricing Appl. Math. Lett. 17 (7) 809 - 814    DOI : 10.1016/j.aml.2004.06.010
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 Space-time adaptive finite difference method for European multi-asset 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 American-style options Electron. Trans. Numer. Anal. 15 (2-7) 165 - 185
Ramage A. , von Sydow L. 2011 A multigrid preconditioner for an adaptive Black-Scholes solver BIT Numer. Math. 51 (1) 217 - 233    DOI : 10.1007/s10543-011-0316-6
Reisinger C. , Wittum G. 2004 On multigrid for anisotropic equations and variational inequalities “pricing multi-dimensional european and american options” Comput. Vis. Sci. 7 (3-4) 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
Saad Y. , van der Vorst H. 2000 Iterative solution of linear systems in the 20th century J. Comput. Appl. Math. 123 (1) 1 - 33    DOI : 10.1016/S0377-0427(00)00412-X
Trottenberg U. , Oosterlee C. , Schüller A. 2001 Multigrid Academic press London, England
van der Vorst H. 1992 Bi-CGSTAB, A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems SIAM J. Sci. & Stat. Comput. 13 (2) 631 - 644    DOI : 10.1137/0913035