Iterative Speedup by Utilizing Symmetric Data in Pricing Options with Two Risky Assets

The Crank–Nicolson method can be used to solve the Black–Scholes partial differential equation in one-dimension when both accuracy and stability is of concern. In multi-dimensions, however, discretizing the computational grid with a Crank–Nicolson scheme requires significantly large storage compared to the widely adopted Operator Splitting Method (OSM). We found that symmetrizing the system of equations resulting from the Crank–Nicolson discretization help us to use the standard pre-conditioner for the iterative matrix solver and reduces the number of iterations to get an accurate option values. In addition, the number of iterations that is required to solve the preconditioned system, resulting from the proposed iterative Crank–Nicolson scheme, does not grow with the size of the system. Thus, we can effectively reduce the order of complexity in multidimensional option pricing. The numerical results are compared to the one with implicit Operator Splitting Method (OSM) to show the effectiveness.


Introduction
The multidimensional Black-Scholes equation is often used to model options written on multiple assets.One of the traditional methods both in practice and research for discretizing multidimensional Black-Scholes equations is the Operator Splitting Method (OSM) [1][2][3].To solve multidimensional Black-Scholes equations, the OSM solves one-dimensional Black-Scholes equations in turn.Thus, it is possible to use a highly efficient tridiagonal matrix solver as in one dimension [2].The OSM converges at first order in time and second order in space if we discretize multidimensional Black-Scholes equations with an implicit central difference method.In one-dimensional cases, in which the option is written on a single asset, the order of convergence in time can be improved to the second order without too many difficulties if one replaces a time integration scheme with Crank-Nicolson.
In multi-dimensions, however, replacing the time integration scheme is not straightforward as in one-dimension.In general, multiple assets are correlated, and thus the multidimensional Black-Scholes equation has corresponding cross partial derivative terms, which do not appear in one-dimensional Black-Scholes equations.Sometimes, these partial derivatives are not calculated and are assumed to be known, which is easy to implement but leads to inaccuracy under high correlation and large volatilities.If an implicit scheme along with OSM is applied to discretize the mixed partial derivatives appearing in the equation, the resulting system becomes no longer tridiagonal and the Thomas algorithm is not applicable.Therefore, another matrix solver or more advanced multidimensional modeling with radial basis functions [4,5] have to be used at the cost of computation time.In practice, however, the Thomas algorithm is indispensable because of its highly efficient nature.Therefore, one avoids fully implicit discretization of multidimensional Black-Scholes equations by replacing the mixed partial derivative terms with known values so that other partial derivative terms can be discretized implicitly.In this way, practitioners could use the OSM with the most efficient matrix solver, the Thomas algorithm, but end up with only first order convergence in time.
Having second order convergence both in space and time is important if highly accurate option values are of concern.A second order convergence in time is helpful for the practitioners if an option has complex payoff structures or parameter adjustment is needed before the maturity.In addition, Greeks (∆, Γ, etc.) are usually calculated in the post-processing stage using the computed option values, and more significant digits on the option price will improve the accuracy of Greeks.However, in practice, the time to obtain one more significant digit on the option price grows exponentially with the first order convergence speed in time.Therefore, it would be natural to try second order schemes such as Crank-Nicolson, BDF-2 [6,7], etc.Each of the second order schemes has its own advantages and disadvantages.We use the Crank-Nicolson scheme, which has second-order convergence in space and time.
A straightforward Crank-Nicolson discretization of the multidimensional Black-Scholes equation produces a system that makes the direct solver unattractive in terms of the computational effort.However, a simple modification to symmetrize the system helps to solve the system more efficiently with an iterative solver.We found that a standard preconditioner for the iterative solver significantly reduces the number of iterations for those problems that we tested after symmetrization.Finally, we compare the computational complexity of the iterative Crank-Nicolson method and OSM.

Nomenclature
We use bold uppercase letters to represent matrix M, bold lowercase letters to denote vector v, superscript' to denote transpose, and superscript (n) to indicate n-th time period.All vectors that are used in this paper are assumed to be a column vector.
Let the price of the derivative V(t, x, y), where x and y are two different asset prices, and V(t, x, y) is the solution of the following two-dimensional Black-Scholes partial differential equation [8,9]: where the domain is defined by {(t, x, y) | t ∈ (0, T ] , x ≥ 0 , y ≥ 0} .In Equation (1), σ i is the volatility of the i-th asset , ρ is the correlation coefficient, and r is the risk-free interest rate.The maturity of the option V is denoted by subscript T. We express the final payoff of the option V in the following form: Big O notation is used to measure the growth rate of algorithm in terms of input size.For example, a function f (N) = O(g(N)) means that there exist positive numbers C and N 0 such that f (N) ≤ C × g(N) for all N > N 0 .

Implicit OSM
Operator Splitting Method (OSM) finds the solution of multi-dimensional version of Equation ( 1) by splitting the differential operator so that the multi-dimensional problem becomes several one-dimensional problems [2].For brevity, let us consider two-dimensional version of OSM.Equation (1) can be rewritten as follows: where Given the final condition V T , we find the solution at previous time, T − ∆t, in two steps: Step (1) Find the solution at T − 1 2 ∆t, by discretizing the equation ∂ t V + L x V = 0 with the given V T : ∆t .
Step (2) Find the solution at T − ∆t, by discretizing the equation ∂ t V + L y V = 0 with the solution found in Step (1) as the given condition: Depending on the size of the time step ∆t, we need to repeat the above Steps (1) and (2) to find the option price at t = 0.In the current presentation, we have only two steps to obtain solution at T − ∆t from time T because we have two operators in Equation (3).For a general m-dimensional problem, we need m steps to obtain a solution at T − ∆t.In each Step, we have to solve the following linear system where M (n) and v (n+1) are known: The above discretization given in Equations ( 5) and ( 6) are not a full implicit discretization because of the partial derivative terms.Note that the mixed partial derivative terms are assumed to be known so that other partial derivative terms can be discretized implicitly.Thus, we should call it semi implicit discretization to be more precise, but, for the rest of this study, we will call the discretization given in Equation (7) an implicit discretization.This implicit discretization is unconditionally stable and has truncation error, O(∆x 2 , ∆t).In addition, the matrix M (n) in Equation ( 7) is tridiagonal (depending on the boundary condition, the matrix M (n) may not be exactly tridigonal; however, it could be converted to tridiagonal [10]).Thus, the system can be effectively solved by Thomas Algorithm [11], which is the main feature of implicit OSM.

Iterative Crank-Nicolson Method
We propose an iterative Crank-Nicolson finite difference discretization of Equation (1) on a general non-uniform grid.Instead of solving a smaller size one-dimensional problem repeatedly, as in the Operator Splitting Method (OSM), we propose to fully discretize the two asset Black-Scholes equation with a Crank-Nicolson scheme.The discretization can be solved efficiently by a GMRES (Generalized Minimal Residual Method) [12,13] solver with preconditioning after symmetrization.
In the following, we use x i and y i for the price of the first and second asset on a (i, j)-th finite difference stencil.We define h i = x i+1 − x i and k j = y j+1 − y j .We approximate the differential operator in the Equation (1) as follows: Applying Equations ( 8)−( 14) to Equation (1), we obtain the following Crank-Nicolson discretization of two-dimensional Black-Scholes equation: where and The symmetrized iterative Crank-Nicolson method for the two asset Black-Scholes equation is described as follows.
Step (1) Get a linear system by discretizing Equation ( 1) with the Crank-Nicolson scheme.The following equation is matrix-vector form of the Equation ( 15) : where Note that the data, L (n) and R (n) , grows quadratically in terms of the grid points.In addition, the system in Equation ( 16) is non-symmetric.
Step (2) Apply appropriate boundary conditions to the L (n) and R (n+1) .Depending on the option type, the boundary condition is either given as a linear boundary condition on the truncated interface or an essential boundary condition where the price of the option is zero.We denote the boundary condition imposed system as follows: Step (3) Symmetrize the system given in Equation ( 19) as follows: where L bc is the transpose of L bc (n) .
Step (4) Create preconditioned matrix P with L bc (n) L bc (n) using incomplete LU factorization (where LU stands for lower and upper triangular matrix).The choice of the preconditioner is more important than the choice of the Krylov iterative method such as GMRES [14,15].The effectiveness of preconditioner P created by incomplete LU factorization is measured by how Step (5) Solve Equation (20) repetitively using GMRES with the preconditioner P and the previous solution vector v (n+1) as an initial guess until we find the option price v (0) .Use the final condition v (N) = V T = V(T, x i , y j ) to start the iteration.We use the following split preconditioning with the incomplete LU factors, P −1 L and P −1 R : The previous solution vector, v (n+1) , is used as an initial guess of v (n) to solve the Equations ( 21) and ( 22).
Thus far, we have explained the general procedure of the iterative Crank-Nicolson method for two asset Black-Scholes equations.The idea can be extended to the three asset Black-Scholes equation: where {(t, x, y, z) | t ∈ (0, T ] , x ≥ 0 , y ≥ 0 , z ≥ 0} .The difference with two asset case is that we have four additional terms to discretize.The discretization is essentially the same with different indices.Thus, we obtain equations that are similar to Equation ( 15) but have 27 terms on each side instead of nine terms.In vector notation, the procedure given above for Steps (1) to ( 5) is the same for the three asset case.
The oscillations in the solution due to non-smooth initial data are a well-known drawback of the Crank-Nicolson method.Thus, if the final condition of the given option has non-smoothness, the computational grid can be prepared so that the option strike price agrees to one of the midpoints in the grid or the v (n+1) in Equation ( 16) can be replaced with ṽ(n+1) by a simple moving average-for example, on a uniform grid, we can use the following equation: (24)

Computational Perspective of the Iterative Crank-Nicolson Method
We will demonstrate the iterative Crank-Nicolson method with some European style options with two assets.In the following, we study computational cost and order of convergence of the iterative Crank-Nicolson method.

Some Numerical Examples
All numerical computations in this section are performed on the finite domain [0, 300] × [0, 300].The relative error (%) in the maximum norm is calculated by the following equation: where u app. is the computed numerical solution and u ref.
is the reference solution.We present numerical results with the following three different final payoffs V T and reference solution V ref [16,17].M is the bivariate normal distribution and K is the strike price.The payoffs in Equations ( 26), ( 28) and (30) are carefully chosen so that (1) V T in Equation ( 26) is symmetric and discontinuous; (2) V T in Equation ( 28) is non-symmetric and discontinuous; and (3) V T in Equation ( 30) is symmetric and continuous: (1) Cash or Nothing with parameters: ρ = 0.5, r = 0.02, K = 75, σ 1 = 0.15, σ 2 = 0.2, T = 1.0.Figure 1a,b shows V T and V ref , respectively.
(3) Basket V ref is approximated with the formula given in [18] with parameters: ρ = 0.5, r = 0.02, K = 150, We use sparse storage for all matrices to hold the data throughout the test.In addition, the drop tolerance is set to 10 −7 in the incomplete LU factorization to create preconditioner and 10 −8 is used for GMRES stopping tolerance.We observed that the solution is reached within two iterations for each time step in all examples that we considered.These tolerances show dependency on the grid size for the problems that we considered and could be optimized, but we did not investigate the effect, as these parameters gave us an accurate numerical solution.Before we proceed to test the iterative Crank-Nicolson method for different payoffs, we show the symmetrization effect, Equation (20).The comparison between non-symmetrized system, given in Equation (19), and the symmetrized system is shown in Figure 2. The cash or nothing payoff is used with those parameters given in Equation ( 27).The slope for both symmetric and non-symmetric case is approximately constant, which suggests that the number of iteration per time steps and the number of iteration per spatical discretization does not grow as grid refinement.Computed solution with the iterative Crank-Nicolson method and its errors for three different options are shown in Figure 3.We can observe there is no oscillation in the computed solution with non-smooth payoffs.The result of numerical tests is summarized in Table 1.The ratio column shows the dropping rate of relative error for both methods.Note that the number of time steps required for OSM is significantly larger to maintain the same level of accuracy compared to the iterative Crank-Nicolson method.For a coarse grid, note that the time required to obtain the same level of accuracy for OSM is shorter than the iterative Crank-Nicolson.In the beginning, when the grid size is still coarse, the iterative Crank-Nicolson scheme takes more time to precondition the system than to actually solve it.However, the iterative Crank-Nicolson method reached the accuracy level faster than OSM as the grid becomes finer.Table 1 suggests that the iterative Crank-Nicolson method becomes more favorable in time if the accuracy is of concern.
We further compare both methods while keeping the ∆x/∆t ratio constant.The results are summarized in Table 2.We see that the iterative Crank-Nicolson method needs about half of the time that is required for OSM to reach the same level of relative error.However, the memory consumption of the iterative Crank-Nicolson is significantly larger than OSM.

Computational Cost
We compare the computational complexity of the iterative Crank-Nicolson method and implicit OSM.Throughout this section, we use N and M for the number of discretization steps for space and time, respectively.Both methods have an order of complexity O (MN 2 ) and storage requirement O (N 2 ), while the former has a second order convergence rate both in space and time and the latter has a second order convergence rate in space and a first order convergence rate in time.Theorem 1.The implicit OSM for the two-dimensional Black-Scholes equation has an order of computational complexity O (MN 2 ).
Proof.Suppose we partition the space and time domain with (N − 1) and (M − 1) intervals.For each time slice, we have to solve 2N (N times for each xand y-direction) tridiagonal systems with the Thomas algorithm, which requires O (N) computational cost.In other words, 2N 2 operations are needed to solve the system of equations for a fixed time period, and we have a total of M time steps.As a consequence, the total computational cost to solve Equation (3) becomes O (MN 2 ).
With O (MN 2 ) computational complexity, the implicit OSM achieves second order convergence in space but first order convergence in time.The computational demand for implicit OSM increases to O (M 2 N 2 ), if we increase the number of time steps to M 2 to obtain an overall second order convergence rate in error defined in Equation (25). Figure 4a-c shows the growth rate of computational complexity measured in time for implicit OSM, where M was chosen to be the same as N.The slope in Figure 4 supports the fact that implicit OSM, which has a second order convergence rate in terms of maximum norm, has computational complexity O (M 2 N 2 ) = O (N 4 ).Proof.Let us assume that we partition the space and time domain with (N − 1) and (M − 1) intervals.Then, we have to solve the system given in Equation ( 16) of size N 2 , for which the preconditioning and GMRES solver in Equations ( 21) and (22) require O (N 2 ) computation for each time slice.Therefore, total computational cost to solve Equation (1), using the implicit Crank-Nicolson method is O (MN 2 ).
The growth rate shown in Figure 4d-f supports the iterative Crank-Nicolson method having computational complexity O (MN 2 ) = O (N 3 ) when M was chosen to be the same as N.
Thus far, we have seen how the computational complexity grows for the iterative Crank-Nicolson method.Memory consumption is another factor that one should consider for the actual computation.Theorem 3. The memory requirement of the fully discretized Crank-Nicolson method is O(N 4 ) for the two-dimensional Black-Scholes equation.
Proof.The fully discretized Crank-Nicolson method, Equation ( 16), requires holding P L , P R , L (n) , R (n+1) in Equations ( 21) and ( 22) for a fixed (n).The exact dimension of these matrices is N 2 × N 2 , requiring approximately N 4 amount of space when fully populated.However, P, L (n) , R (n+1) are all sparse matrices so that we can relax the storage requirement to N 2 by storing only the non-zeros.The actual number of non-zeros are on the order of N 2 , not N 4 .The ratio of non-zeros in the matrix, P L , P R , L (n) , R (n+1) , are similar.To see the growth of non-zeros in these matrices in terms of N, see, for example, Figure 5.The non-zero terms in sparse matrix (a) L (n) , (b) R (n+1) , and (c) P grow quadratically in number of grid points.The dotted line shown in Figure 5 has slope 2 and shows us that the actual storage requirement can be reduced to the order of N 2 from N 4 .Proof.To solve the tridiagonal system given in Equation ( 7), we only need 4N spaces in memory.However, to do this, we have to hold the data v (n+1) .Since the entire computational domain is a grid of size N × N, N 2 amount of storage is needed in memory to hold v (n+1) throughout the time period, ∆t, to obtain v (n) in Equation (7).Therefore, the amortized cost is on the order of N 2 .
Before closing this section, we summarize the computational complexity and memory requirement for the implicit OSM and iterative Crank-Nicolson method.The former is first order, and the latter is second order method in time.Both have the computational complexity of O (MN 2 ).However, implicit OSM has O (N 2 ) memory requirement and iterative Crank-Nicolson has O (N 4 ) memory need for the two-dimensional Black-Scholes equation.Therefore, the growth of the data can not be compared for large N but, after symmetrizing and using a sparse storage structure, we can reduce the storage need to O (N 2 ) for the iterative Crank-Nicolson method, which is worth the effort.

Order of Convergence in Space and Time
The numerical solution obtained by implicit OSM has a second order convergence rate in space and a first order convergence rate in time because the truncation error of the finite difference discretization given in Section 3 are all in second order and first order in time (see, for example, Figure 6a-c).The relative error in maximum norm for implicit OSM drops with slope one, which means it is first order in time.On the other hand, the iterative Crank-Nicolson method has a second order convergence rate both in space and time.The slope of the relative error curve, shown in Figure 6a-c, support the proposed method that has a second order convergence rate in time for all three options that we have considered while the implicit OSM converges first order in time.The second order of convergence can be proved by calculating the truncation error.The truncation error of iterative Crank-Nicolson discretization with non-uniform grid size is indeed second order in space and time (see Appendix in [19]).

Conclusions
A second order method is essential for pricing options when highly accurate option price is needed in time.The second order convergence even in the simplest case, a European style, is of great importance in practice when the option has complex payoff structures.For example, a multidimensional Black-Scholes equation has to be solved many times in a row to price and hedge an equity-linked security.A three-dimensional extension of the current study would be an interesting and important future work.A straightforward implementation of the multi-dimensional Crank-Nicolson scheme could be thought inefficient.However, with a simple symmetrization and a standard preconditioner, we have found that the order of complexity to solve the system is about the same for the fully discretized iterative Crank-Nicolson scheme and (semi-)implicit Operator Splitting Method.In other words, the order of complexity for the first and second order method turned out to be the same.However, note that the second order method, the iterative Crank-Nicolson method, needs more storage compared to the first order method, OSM, but reaches the same level of error in significantly less time.It is interesting to observe the trade-off between storage and computational time in the context of option pricing.

Figure 1 .
Figure 1.The surface of option payoff V T and V ref : (a) cash or nothing V T ; (b) cash or nothing V ref ; (c) two asset call V T ; (d) two asset call V ref ; (e) basket call V T ; and (f) basket call V ref .

Figure 2 .
Figure 2. The effect of symmetrization: (a) number of time steps (M) versus total iterations; and (b) number of spatial discretization versus total iterations.

Theorem 2 .
The iterative Crank-Nicolson method for the two-dimensional Black-Scholes equation has the computational complexity O (MN 2 ).

Figure 4 .
Figure 4.The order of complexity measured in time for implicit OSM with a different payoff: (a) cash or nothing; (b) two asset call; (c) basket call.The order of complexity measured in time for the iterative Crank-Nicolson method with a different payoff: (d) cash or nothing; (e) two asset call; and (f) basket call.The dotted line shows the growth rate.

Theorem 4 .
The memory requirement of the implicit OSM is O (N 2 ) for the two-dimensional Black-Scholes equation.

Figure 5 .
Figure 5. Number of non-zeros in the matrices, given in Equation (16), versus N, the number of grid points: (a) L (n) ; (b) R (n+1) ; and (c) P. The quadratic growth rate is shown with the dotted line.

Figure 6 .
Figure 6.The relative error measured in maximum norm versus the time step size for: (a) cash or nothing option; (b) two asset call option; and (c) basket call option.

Table 1 .
Computational time comparison between iterative Crank-Nicolson and Operator Splitting Method (OSM) to reach targeted accuracy while maintaining convergence rate.

Table 2 .
CPU time, Memory, and Relative Error comparison with the use of identical ∆x and ∆t for both methods.The memory and CPU time are accumulated for the entire simulation.