A Numerical Study of RBFs-DQ Method for Multi-Asset Option Pricing Problems

abstract: In this paper, we propose a numerical scheme to solve multi-dimensional Black-Scholes equation using the global radial basis functions-based differential quadrature (RBFs-DQ) method. Before applying the method, it is needed to remove mixed derivatives from the Black-Scholes equation by making an appropriate change of variables. Then, any spatial derivativeis are approximated by a linear weighted sum of all the function values in the whole physical domain. In the RBFs-DQ method the weighting coefficients are computed by RBFs. The method is very easy to implement and the non-singularity is ensured. The proposed method combines the advantages of the conventional DQ method and the RBFs. It also preserves mesh-free feature of RBFs.


Introduction
In the past several decades, financial derivative products have become increasingly important in the world of finance.Options as a kind of important financial derivatives have a wide range of applications that one of the major concerns in financial markets is determining the value of options.An option contract is an agreement between a buying party (the holder) and a selling party (the underwriter).The holder of the option contract has no obligation to use his option contract, whereas the underwriter is obliged to agree with the holder uses the option contract.Since an investor of options can set the strike price of underlying assets in advance, these options are powerful tools for hedging risk in financial market, therefore, the trading volume of options are increasing all over the world.Options occur in many forms.Examples are vanilla options, barrier options, digital options and multi-asset options.In this paper, in order to implement the method, we consider options where the payoff depends on multi underlying assets.Such option pricing can be modelled by higher dimensional generalization of the original Black-Scholes equation in which Black and Scholes proposed an explicit formula for evaluating European call options without dividends [1] and extended in [2].By assuming that the asset price is risk-neutral, Black and Scholes showed that the European call options value satisfies a log-normal diffusion type partial differential equation (PDE) which is known now as the Black-Scholes equation.If we introduce the notation, C(s 1 , s 2 , • • • , s d , t) to represent the option value at time t which stock prices s 1 , s 2 , • • • , s d , then we get the following linear parabolic PDE known as the Black-Scholes equation for multi-asset option problems.Where s i ,σ i ,ρ ij , and r are respectively, i-th asset price, volatility of the i-th asset price, correlation between the prices of i-th and j-th assets, and risk free interest rate.Since then, methods for option pricing have been discovered and improved by many scholars.Details can be found in reference [3] which is a good review of various models and applications to the option pricing.According to [4], there have been numerous attempts to find the analytic form of solutions of the Black-Scholes equation for various derivative products in financial world.However, it is not easy to get an analytic form of solution for most of financial derivatives because of the complexity of the financial product itself and the system of the financial market.Therefore, in finance, numerical methods such as finite difference method, finite element method and Monte Carlo simulation techniques for pricing derivatives problems, cf.[5,6,7,8,9,10,11,12], have been used.As it has been mentioned in [13], different numerical methods can be applied to price multi-variate derivatives.Higher dimensional generalizations of lattice binomial methods can be used, cf.[14], where European options based on three underlying assets are solved numerically.Stable higher order methods for the Black and Scholes equation have been introduced by [15] and [16], mesh-free methods based on RBFs may also reduce the computational efforts significantly, see [17].In [4] Jo and Kim combined the operator splitting method with parallel computation technique to solve the multi-dimensional Black-Scholes equations.Nielsen, Skavhaug and Tveito, [13] used Penalty methods for the numerical solution of American multi-asset option problems.Authors [18] proposed explicit Runge-Kutta methods for multi-asset American options in 2014.
Numerical Study of RBFs-DQ Method

11
In this work, the global RBFs-DQ method is utilized for the numerical solution of multi-asset option pricing problems.As pointed out in [19], the classical form of the differential quadrature (DQ) method was introduced by [20,21] for the numerical solution of PDEs.The basic idea of the DQ method is that any derivative at a mesh point can be approximated by a weighted linear sum of all the function values along a mesh line.The key procedure in the DQ method is the determination of weighting coefficients for any order derivative discretization.As shown by [22], when the solution of a PDE is approximated by a high order polynomial, the weighting coefficients can be computed by a simple algebraic formulation or by a recurrence relationship.Later, [23] also showed that when the solution of the PDE is approximated by a Fourier series expansion, the weighting coefficients of the first-and second-order derivatives can be computed explicitly by algebraic formulations.The details of the polynomial-based and Fourier series expansionbased DQ methods can be found in the book of [24].The DQ method has been extensively applied in engineering for the rapid and accurate solution of various linear and non-linear differential equations [22,23,24,25,26,27,28,29,30].In general, the polynomial-based and Fourier series expansion-based DQ methods can achieve very accurate results by using a considerable small number of grid points.It is a robust and efficient technique, but the dependency of the method on a mesh leads to complications.To overcome this difficulty, using advantages of the RBF, [19] proposed other set of mesh-less methods which is named the RBFs-DQ method, which combines the mesh-free nature of RBFs with the derivative approximation of DQ method.It seems that the multi-dimensional polynomial approximation as the test function may not be a good choice in the DQ approximation.As it has been shown in [19], RBFs, which have truly mesh-less property and insensitivity to high dimension, could be a good choice in the DQ approximation.The advantages of RBFs-DQ method as combination of DQ approximation and RBFs provide an efficient discretization method, which is a derivative approximation approach and is mesh-free.In this method, the RBFs are taken as the test functions in the DQ approximation to compute the weighting coefficients.Moreover, the method not only inherits the advantages of the DQ method such as high accuracy and efficient computation, but also owns the merits of RBFs such as mesh-free feature and easy extension to high dimension.
The organization of this paper is as follows.In Section 2 we describe RBFs-DQ method for partial derivatives approximation.In Section 3 we apply the method on multi-dimensional option pricing problems by applying Theorem 3.1.Stability of the method to solve multi-asset option pricing problems is studied in Section 4. In the Section 5, the numerical example for two-asset option pricing with this method is given.Finally, Section 6 is dedicated to a brief conclusion.Note that we have computed the numerical results by Matlab programming.

RBFs-based DQ method
In this section, we present the global RBFs-DQ method in detail.In the following, we give the details step by step.

Radial basis functions
In the past decades, the interpolation theory of RBF has undergone intensive research, and nowadays RBFs play an increasingly important role in the field of reconstructing functions from multivariate scattered data.In general, the interpolation theory of RBF can be described as follows: If an unknown function f (X) is only known at a finite set of centers X i , i = 1, ..., N , the approximation of a function f (X) can be written as a linear combination of N , RBFs where N is the number of data points, X = (x 1 , x 2 , ..., x d ), d is the dimension of the problem, λ's are coefficients to be determined and ϕ is the RBF.Also, Eq.
(2.1) can be written without the additional polynomial ψ.The success of RBF interpolant is due to its excellent performance.Based on numerical experiments, [31] gave a comprehensive review on the interpolation methods for scattered data.
From the numerical tests, Franke found that RBFs performed better than other tested methods regarding accuracy, stability, efficiency, memory requirement, and simplicity of implementation.Among the tested RBFs, multi-quadrics (MQs) yields the most accurate results.MQ-RBFs can be written as where r j = |X − X j | is the usual Euclidean distance and c j is a shape parameter (for generality, it is written with a variable shape parameter c j ) and given by the practitioners.It is noted that the shape parameter plays a very important role in the MQ formulation.Its value determines the fundamental shape of the basis function.These effects were initially observed for scattered data interpolation, but we will see that they also occur in the numerical solution of the PDEs.Thus, the problem of how to select a good value for the parameter c is a key question.Several methods for selecting c for the MQ interpolants in the two-dimensional case were suggested in the literature.[32] used c = 0.815d where d = 1/N N i=1 d i and d i is the distance between the i-th node and its neighboring node.The parameter d is replaced by D/ √ N in [31], where D is the diameter of the minimal circle enclosing all supporting points and suggested to use 1.25D/ √ N .Optimization of shape parameter and its distribution are still active research field.If ψ d q denotes the space of d-variate polynomials of order not exceeding q, and letting the polynomials P 1 , • • • , P m be the basis of ψ d q in R d , then the polynomial ψ(X) in Eq. (2.1), is usually written in the following form where m = (q−1+d)!(d!(q−1)!) .To determine the coefficients (2.4)

DQ method for partial derivative approximation in multi-dimension
In this section,we explain the method to approximate the spatial derivative in two-dimension.However, the method can be generalized for multi-dimensional cases.
We note that the basic idea of the DQ method is that any derivative can be Figure 1: A structured mesh for two-dimensional problem approximated by a linear weighted sum of function values at some mesh points.We can keep this idea but release the choice of function values along a mesh line in the conventional DQ approximation.In other words, for a two-dimensional problem shown if Fig. 1, any spatial derivative is approximated by a linear weighted sum of all the function values in the whole two-dimensional domain.In this approximation, a mesh point in the two-dimensional domain is represented by one index, k, while in the conventional DQ approximation, the mesh point is represented by two indexes i, j.If the mesh is structured, it is easy to establish the relationship between i, j and k.For the example shown in Fig. 1, k can be written as k = (i − 1)N 2 + j, i = 1, 2, ..., N 1 ; j = 1, 2, ..., N 2 .Clearly, when i is changed from 1 to N 1 , and j is changed from 1 to N 2 and k is changed from 1 to N = N 1 × N 2 .The DQ approximation for the n-th order derivative with respect to x, ∂ n f ∂x n , and the m-th order derivative with respect to y, ∂ m f ∂y m at (x k , y k ) can be written as In [19], Shu et al. employed this technique and made a novel and effective algorithm for the use of RBFs to solve the PDEs.Instead of using polynomials for determining coefficients, they applied RBFs as test functions.
In the following subsection, we will show that the weighting coefficients in Eq. (2.5) and Eq.(2.6) can be determined by the function approximation of RBFs and the analysis of a linear vector space.

RBFs-DQ approximation
In this subsection, we will use the MQ-RBFs as test functions to determine the weighting coefficients in the DQ approximation of derivatives for a two-dimensional problem.However, the method can be easily extended to the multi-dimensional case as it is dimension-independent, and other RBFs can also be used as test functions.
Suppose that the solution of a PDE is continuous, which can be approximated by MQ-RBFs, and only a constant is included in the polynomial term ψ(x, y).So, for a two-dimensional case, Eq. (2.1) can be reduced to (2.7) To make the problem well-posed, one more additional equation is required.From Eq. (2.4), we have Substituting Eq. (2.8) into Eq.(2.7) gives (2.9) The number of unknowns in Eq. (2.7) is reduced to N .When no confusion rises, µ can be replaced by λ i and Eq.(2.9) can be written as (2.10) can be further written as λ j g j (x, y) + λ i . (2.11) Eq. (2.11) constructs an N -dimensional linear vector space V N .A set of base functions in V N can be taken as From the concept of linear independence, the bases of a vector space can be considered as linearly independent subset that spans the entire space.From the property of a linear vector space, if all the base functions satisfy Eq. (2.5) and Eq.(2.6), so does any function in the space V N represented by Eq. (2.11).There is an interesting feature.From Eq. (2.11), while all the base functions are given, the function f (x, y) is still unknown since the coefficients λ i are unknown.However, when all the base functions satisfy Eq.
(2.6).In other words, we can guarantee that the solution of a PDE approximated by the RBF satisfies Eq. (2.5) and Eq.(2.6).Thus, when the weighting coefficients of DQ approximation are determined by all base functions, they can be used to discretize the derivatives of a PDE.That is the essence of the RBFs-DQ method.
Using the idea of DQ method, the weighting coefficients of the n−th order partial derivatives can be determined by substituting all the base functions q 1 , q 2 , • • • , q N , into Eq.(2.5), as ) For the given i equation system (2.12) and (2.13) has N unknowns with N equations.The matrix form for the weighting coefficients can be written as where Therefore, the coefficient vector W (n)x i can be obtained by Then, the coefficient vector W (n)x i can be used to approximate the n-th-order derivative in the x direction for any unknown smooth function at node i.Clearly, there exists a unique solution only if the collocation matrix [Q] is non-singular.The non-singularity of the collocation matrix [Q] depends on the properties of used RBFs.[33] proved that matrix [Q] is conditionally positive definite for MQ-RBFs.This fact can not guarantee the non-singularity of matrix [Q].[34] showed that cases of singularity are quite rare, and not serious objection to a valuable numerical technique.In a similar manner, all the base functions are substituted into Eq.(2.6) to approximate the m-th-order derivative in the y direction for any unknown smooth function at node i.For MQ-RBFs, q j (x, y) can be written as where second-and higher-order derivatives of q j (x, y) can also be obtained by differentiating Eq. (2.16).One of the most attractive properties in the above method is that the weighting coefficients are only related to the test functions and the position of the collocation points.

Global RBFs-DQ method for multi-dimensional Black-Scholes equation
In this section, first we will introduce the variable transformations s i = e xi and C(s 1 , ..., s d , t) = V (x 1 , ..., x d , t).By this change of variable Eq. (1.1) becomes the parabolic PDE with constant coefficients as as Eq.(3.1) satisfies the conditions of the following theorem, we can remove the crossing terms then we can apply the RBFs-DQ method to option pricing.
assuming that the mixed partial derivatives are equal, we may as well assume that a ij = a ji , by making an appropriate change of variables, we can write the top-order term Therefore, the multi-dimensional Black-Scholes equation transformes into The domain is discretized by taking N knots according to the method used in previous section and we would like to discrete the Eq.(3.2) with respect to time.We introduce a weight θ and apply the θ implicit scheme to the problem as where represents the approximation of the function value , and the time t n .By applying the RBFs-DQ method, Eq. (3.2) can be rewritten in a discrete form as follows: where k = 1, 2, ..., N and w (1) ξ i kj and w (1) ξ i kj represents the computed weighting coefficients in the DQ approximation.

Stability of the method
In this section, we study the stability of the implicit finite difference method described above.Let us assume that U be exact and Ũ is the numerical solution of equation (3.3).Note that the discrete equations in the last section can be rewritten as the following form So we can write Eq. (4.1) as where L. Khodayari and M. Ranjbar is an N × N matrix determined by descretized form of Eq. (3.2) and The scheme for initial value problem is stable if and only if there exists a positive constant C independent of the mesh spacing and initial data such that If E depends on n we get a product of the operators at each time level.Taking a norm Therefore, the numerical scheme is stable if and only if there exists a positive constant C such that Since ρ(E) as the spectral radius of E provides a lower bound to any matrix norm, for the scheme to remain stable, we should have ρ(E) ≤ 1 or equivalently we can say that which holds for η D are located in the right half plane, where η D are eigenvalues of the matrix D. Inequality (4.3) also shows that stability of the scheme, in the case of RBFs with shape parameter like MQ, depends upon shape parameter.

Numerical computation of two asset European option
For two-dimensional option pricing under Geometric Brownian motion framework, consider the PDE (1.1) with d = 2 and the final time condition for European call is where E, T are exercise price and maturity time.The following boundary conditions are imposed where s max 1 , s max 2 are respectively maximum of s 1 , maximum of s 2 .In this case, the variables transformation are s 1 = e x , s 2 = e y and C(s 1 , s 2 , t) = V (x, y, t).Under this change of variables the two-dimensional Eq. (1.1) in which d=2 becomes ∂V ∂y − rV = 0. (5.4) Numerical Study of RBFs-DQ Method

19
Now we can remove the crossing term ρσ 1 σ 2 ∂ 2 V ∂x∂y .In this way, we apply characteristic line to remove mixed derivatives in Eq. (5.4) by considering its characteristic equation for the above PDE [35] as where i = √ −1.Thus we observe the following relation between spatial variables x and y, Further, we have Hence, Thus, the Black-Scholes equation in (5.4) ransformed into where ). where the mixed derivative term vanishes.Eq. (5.5) can be rewritten in a discrete form as follows: where i = 1, 2, ..., N .For illustration of the accuracy of the proposed method, let s 1 ∈ [e −3.5 , e 4 ], s 2 ∈ [e −3.5 , e 4 ], and θ = 0.5.To study the behavior of the method, different structured and random mesh sizes are used for the point distribution in ξ and η direction with shape parameter c 2 = 1.25(D/ √ N ) as selected by Franke (1982), where D is the diameter of the minimal circle enclosing all supporting points.Fig. (2) and Fig. (3), show that our scheme can provide reasonable approximations for Stulz method [36] which provides a closed form solution.The root-mean-square-error (RMSE) defined by the where U Stz. is the solution by Stulz Method, Ũnum. is solution by numerical approximations.As ρ(E), the spectral radius of matrix E, is equivalent or less than unity in each simulation, it guarantees the stability of the method.

Conclusion
In this paper, A mesh-free RBFs-DQ method is presented and used for numerical solution of multi-dimensional Black-Scholes equation.In our approach, any spatial derivative is approximated by a linear weighted sum of all the function values in the whole physical domain.The weighting coefficients in the method are determined by RBF approximation and linear vector space analysis.The proposed method is similar to finite difference schemes in the sense of derivative approximation.Numerical results showed that our RBFs-DQ scheme is an efficient approach for solution of multi-asset option pricing.
extra m equations are required in addition to the N equations resulting from Numerical Study of RBFs-DQ Method 13 the collocating Eq. (2.1) at the N points.This is ensured by the m conditions for Eq.(2.1), N j=1

Theorem 3 . 1 .
Consider the second-order equation n i,j=1 a ij u xixj as n k=1 d k u x k x k , where the coefficients d k are the eigenvalues of the n × n matrix A = (a ij ).