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. 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. Let s_i ( 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 ₁ , s ₂ , ... , 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 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 ₁ and y = s ₂ . 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 ₁ , s ₂ , T - 𝜏): We use the following linear boundary conditions on all boundaries, 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 ((i-0.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 , 𝜕 ² u /𝜕 x ² ≈ ( u _{i-1, j} - 2 u_ij + u _{i+1, j} )/ h ² , 𝜕 ² u /𝜕 x 𝜕 y ≈ ( u _{i+1, j+1} + u_ij - u _i,j+1 - u _{i+1, j} )/ h ² , and 𝜕 u /𝜕 τ ≈ ( u ⁿ⁺¹ - uⁿ )/∆ τ . 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: where Next, to renumber the multi-indexed data u_ij as the single-indexed data U_l , we denote by U_l = U _Nx(j-1)+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 ⁿ⁺¹ = , b ⁿ = , 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 ⁰ and proceeds as follows: Bi-CGSTAB cycle Define the maximum number of iteration ITER and the error tolerance TOL Set r ⁰ = b - A U ⁰ , = r ⁰ , 𝜌 ⁰ = α = 𝜔 ⁰ = 1, v ⁰ = p ⁰ = 0, k = 1 While ( k ≤ ITER & || r ^k || ₂ > TOL ) 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: 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 }. ­Ω _k-1 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 . 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 ( ) = for each ( i, j ) ∈ Ω­ _k , where L ( ) = - ∆τ𝓛 _BS Given the number 𝜈 ₁ and 𝜈 ₂ 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 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 𝜈 ₁ Gauss-Seidel relaxation steps. Step 2) Coarse grid correction Step 3) Postsmoothing: = SMOOTH^ν2 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 X₁ and asset two, y, is above strike X₂ 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 closed-form solutions [9] are Case 1 : u ( x , y , T ) = K_e^-rTM ( α , β ; 𝜌), Case 2 : u ( x , y , T ) = K_e^-rTM ( -α , -β ; 𝜌), Case 3 : u ( x , y , T ) = K_e^-rTM ( -α , β ), -𝜌), where [10] . Let 𝜌 be the correlation between the two variables, then We computed the numerical solution on uniform grids, h = 300/2 ⁿ 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 ₁ = 100, X ₂ = 100 and cash amount K = 1. The volatilities are 𝜎 ₁ = 0.3, 𝜎 ₂ = 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 ₂ -norm of the error, || e || ₂ : 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 ₂ 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 ₂ 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. 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.