\documentclass[twoside]{article} \usepackage{graphics} \usepackage{graphicx} \usepackage{latexsym} \usepackage{amssymb} \usepackage{amsmath} \usepackage{BSPMsample} \begin{document} \title[BSPM \LaTeX{} Style] {A Numerical Study of RBFs-DQ Method for Multi-Asset Option Pricing Problems } %between [] goes the title that will appear on the top of every %odd page, between {} goes the real title. %\thanks{...} put your grant \author[L. Khodayari and M. Ranjbar]{Leila Khodayari and Mojtaba Ranjbar} %between [] goes the authors that will appear on the top of every %even page, between {} goes the real author. \address{Mojtaba Ranjbar\\Department of Mathematics, \\ Azarbaijan Shahid Madani University \\ Tabriz - Iran\\ \email: m\_ranjbar@azaruniv.edu} \maketitle \begin{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 derivatives 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 remains mesh-free feature of RBFs. \end{abstract} %%put here Key words \keywords Radial basis functions, multi-dimensional Black-Scholes equation, differential quadrature, European option \tableofcontents %% put here the subject class \2000mathclass{91G20, 65M06, 65M12} \section{Introduction} \label{sec1} 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 which Black and Scholes proposed an explicit formula for evaluating European call options without dividends \cite{black1973pricing} and extended by \cite{merton1973theory}. By assuming that the asset price is risk-neutral, Black and Scholes showed that the European call options value satisfies a lognormal 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,\cdots, s_d, t) $ to represent the option value at time $t$ and stock prices $s_1, s_2,\cdots, s_d$, then we get the following linear parabolic PDE \begin{equation}\label{eq1} \frac{\partial C}{\partial t}+ \dfrac{1}{2} \sum_{i,j=1}^{d} \rho_{ij} \sigma_i\sigma_j s_i s_j\frac{\partial^2 C}{\partial s_i \partial s_j}+ \sum_{i=1}^{d} r s_i \frac{\partial C}{\partial s_i}-rC=0, \end{equation} known as the Black-Scholes equation for multi-asset option problems. Where $s_i$,$\sigma_i$,$\rho_{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 \cite{broadie2004anniversary} which is a good review of valuation models and applications to the option pricing. According to \cite{jo2013comparison}, there have been numerous attempts to find the analytic form solutions of the Black-Scholes equation for various derivative products in financial world. However, it is not easy to get an analytic form 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 (FDM), Finite Element Method (FEM) and Monte Carlo (MC) simulation techniques for pricing derivatives problems, cf. \cite{glasserman2004monte,ibanez2004monte,choi2001numerical,choi2003valuation,clarke1999multigrid,forsyth1999finite,wilmott1993option,zhu1999singularity}, have been used. As has been mentioned in \cite{nielsen2008penalty}, different numerical methods can be applied to price multi-variate derivatives. Higher dimensional generalizations of lattice binomial methods can be used, cf. \cite{boyle1989numerical}, 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 \cite{voss2003fourth} and \cite{khaliq2006linearly}, mesh-free methods based on RBFs may also reduce the computational efforts significantly, see \cite{fasshauer2004using}. In \cite{jo2013comparison} Jo and Kim combined the operator splitting method with parallel computation technique to solve the multi-dimensional Black-Scholes equations. \cite{nielsen2008penalty} used Penalty methods for the numerical solution of American multi-asset option problems. \cite{martin2014stabilized} proposed explicit Runge-Kutta methods for multi-asset American options in 2014. In this work, the global RBFs-DQ method is utilized for the numerical solution of multi-asset option pricing problems. As printed out in \cite{shu2004solution}, the classic form of the differential quadrature (DQ) method was introduced by \cite{bellman1971differential,bellman1972differential} 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 \cite{shu1992application}, 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, \cite{shu1997fourier} 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 expansion-based DQ methods can be found in the book of \cite{shu2000differential}. The DQ method has been extensively applied in engineering for the rapid and accurate solution of various linear and non-linear differential equations \cite{shu1992application,shu1997fourier,shu2000differential,bert1988two,sherbourne1991differential,gutierrez1994method,bert1996differential,han1997eight,shu2000free}. 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 difficult, using advantages of the RBF, \cite{shu2004solution} proposed other set of meshless methods which is named the RBFs-DQ method, which combines the mesh-free nature of RBFs with the derivative approximation of differential quadrature (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 has been shown in \cite{shu2004solution}, RBFs, which have truly meshless property and insensitivity to high dimension, could be a good choice in the DQ approximation. The advantages of RBF-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 \ref{teo1}. Stability of the method to solve multi-asset option pricing problems is studied in Section 4. In the following 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. %%%%%%%%%%%%%%%%%%5 \section{RBFs-based DQ method } \label{sec2} In this section, we will show in detail the global RBFs-DQ method. In the following, we will show the details step by step. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{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 centres ${X_i}, i=1, . . . , N$, the approximation of a function $f(X) $ can be written as a linear combination of $N$, RBFs \begin{equation}\label{eq3} f(X)\cong \sum_{j=1}^{N}\lambda_j\varphi(X,X_j)+\mathbb{P}(X),\quad X\in\Omega\subset \mathbb{R}^d, \end{equation} where $N$ is the number of data points, $X=(x_1, x_2,..., x_d)$, $d$ is the dimension of the problem, $\lambda$'s are coefficients to be determined and $\varphi$ is the RBF. Also, Eq. (\ref{eq3}) can be written without the additional polynomial $\mathbb{P}$. The success of RBF interpolant is due to its excellent performance. Based on numerical experiments, \cite{franke1982scattered} 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 henceforth) yields the most accurate results. MQ-RBFs can be written as \begin{equation}\label{eq4} \varphi(X,X_j)=\varphi(r_j)=\sqrt{r_j^2+c_j^2}, \end{equation} where $r_j=\vert X-X_j\vert$ 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$ appears in front of us. Several methods for selecting $c$ for the MQ interpolants in the two-dimensional case were suggested in the literature. \cite{hardy1971multiquadric} used $c=0.815d$ where $d=1/N\sum_{i=1}^{N} d_i$ and $d_i$ is the distance between the $i$ th node and its neighboring node. \cite{franke1982scattered} replaced $d$ by $D/\sqrt{N}$ where $D$ is the diameter of the minimal circle enclosing all supporting points and suggested to use $1.25D/\sqrt{N}$. Up to now, optimization of shape parameter and its distribution are still under research. If $\mathbb{P}_q^d$ denotes the space of $d$-variate polynomials of order not exceeding $q$, and letting the polynomials $P_1, \cdots, P_m$ be the basis of $\mathbb{P}_q^d$ in $\mathbb{R}^d$, then the polynomial $\mathbb{P}(X)$ in Eq. (\ref{eq3}), is usually written in the following form \begin{equation}\label{eq5} \mathbb{P}(X)=\sum_{i=1}^{m}\xi_iP_i(X), \end{equation} where $m=\frac{(q-1+d)!}{(d!(q-1)!)} $. To determine the coefficients $(\lambda_1,\lambda_2,\cdots, \lambda_N)$ and $(\xi_1, \xi_2, \cdots, \xi_m)$, extra $m$ equations are required in addition to the $N$ equations resulting from the collocating Eq. (\ref{eq3}) at the $N$ points. This is ensured by the $m$ conditions for Eq. (\ref{eq3}), \begin{equation}\label{eq6} \sum_{j=1}^{N}\lambda_j P_i(X_j)=0, i=1,\cdots, m. \end{equation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{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. It is noted that the basic idea of the DQ method is that any derivative can be 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, 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, $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\times N_2$. The DQ approximation for the $n$-th order derivative with respect to $x,\frac{\partial^n f}{\partial x^n}$, and the $m$-th order derivative with respect to $y, \frac{\partial^m f}{\partial y^m}$ at $(x_k,y_k)$ can be written as \begin{equation}\label{eq8} \frac{\partial^n f}{\partial x^n}(x_k,y_k)=\sum_{s=1}^{N}w_{k,s}^{(n)_x}f(x_{s},y_{s}), \end{equation} \begin{equation}\label{eq9} \frac{\partial^m f}{\partial y^m}(x_k,y_k)=\sum_{s=1}^{N}w_{k,s}^{(m)_y}f(x_{s},y_{s}). \end{equation} \cite{shu2004solution} 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. (\ref{eq8}) and Eq. (\ref{eq9}) and can be determined by the function approximation of RBFs and the analysis of a linear vector space. %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{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 in Eq. (\ref{eq4}), and only a constant is included in the polynomial term $\mathbb{P}(x,y)$. So, for a two-dimensional case, the function approximation can be written as \begin{equation}\label{eq10} f(x,y)= \sum_{j=1}^{N} \lambda_j \varphi(\vert (x,y)-(x_j,y_j)\vert)+ \mu. \end{equation} Eq. (\ref{eq10}) has $(N +1)$ unknown coefficients and can only be applied at $N$ nodes. So, we need an additional equation to close the system. To make the problem be well-posed, this additional equation can be made by Eq. (\ref{eq6}) ( the sum of the expansion coefficients to be zero). As a result, we have \begin{equation}\label{eq11} \sum_{j=1}^{N}\lambda_j =0\Longrightarrow \lambda_i=-\sum_{j=1, j\neq i}^{N}\lambda_j. \end{equation} Substituting Eq. (\ref{eq11}) into Eq. (\ref{eq10}) gives \begin{equation}\label{eq12} f(x,y)=\sum_{j\neq i, j,i=1}^{N}\lambda_j(\varphi(\vert (x,y)-(x_j,y_j)\vert)-\varphi(\vert (x,y)-(x_i,y_i)\vert))+\mu. \end{equation} The number of unknowns in Eq. (\ref{eq10}) is reduced to $N$. As no confusion rises, $ \mu $ can be replaced by $\lambda_i$ and Eq. (\ref{eq12}) can be written as \begin{equation}\label{eq13} f(x,y)=\sum_{j\neq i, j,i=1}^{N}\lambda_j(\varphi(\vert (x,y)-(x_j,y_j)\vert)-\varphi(\vert (x,y)-(x_i,y_i)\vert))+\lambda_i \end{equation} By setting $ g_j(x,y)=\varphi(\vert (x,y)- (x_j,y_j)\vert)-\varphi(\vert (x,y)- (x_i,y_i)\vert), j=1, \cdots,i-1,i+1,\cdots , N$, Eq. (\ref{eq13}) can be further written as \begin{equation}\label{eq14} f(x,y)=\sum_{{j\neq i=1}}^{N}\lambda_j g_j(x,y)+\lambda_i. \end{equation} The form of Eq. (\ref{eq14}) constructs an $N$-dimensional linear vector space $\textbf{V}^N$. A set of base functions in $\textbf{V}^N$ can be taken as $q_i=1, q_j= g_j(x,y)=\varphi(\vert (x,y)- (x_j,y_j)\vert)-\varphi(\vert (x,y)- (x_i,y_i)\vert), j=1, \cdots,i-1,i+1,\cdots , N.$ 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 the linear Eq. (\ref{eq8}) and Eq. (\ref{eq9}), so does any function in the space $\textbf{V}^N$ represented by Eq. (\ref{eq14}). There is an interesting feature. From Eq. (\ref{eq14}), while all the base functions are given, the function $f(x,y)$ is still unknown since the coefficients $\lambda_i$ are unknown. However, when all the base functions satisfy Eq. (\ref{eq8}) and Eq. (\ref{eq9}), we can guarantee that $f(x,y)$ also satisfies Eq. (\ref{eq8}) and Eq. (\ref{eq9}). In other words, we can guarantee that the solution of a PDE approximated by the RBF satisfies Eq. (\ref{eq8}) and Eq. (\ref{eq9}). Thus, when the weighting coefficients of DQ approximation are determined by all the base functions, they can be used to discretize the derivatives in 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,\cdots,q_N,$ into Eq. (\ref{eq8}), as \begin{equation}\label{eq15} \frac{\partial^nq_i(x_i,y_i)}{\partial x^n}=\sum_{k=1}^{N}w_{i,k}^{(n)_x}=0, \end{equation} \begin{equation}\label{eq16} \frac{\partial^nq_j(x_i,y_i)}{\partial x^n}=\sum_{k=1}^{N} w_{i,k}^{(n)_x}q_j(x_k,y_k), ~~~~j=1,2, \cdots, N, j\neq i. \end{equation} For the given $i$ equation system (\ref{eq15}) and (\ref{eq16}) has $N$ unknowns with $N$ equations. The matrix form for the weighting coefficients can be written as \begin{equation}\label{eq17} {\textbf{q}_i}^{(n)_x}=[\textbf{Q}]\mathbf{W}_i^{(n)_x}, \end{equation} where $$ \textbf{q}_i^{(n)_x}=[0,\frac{\partial^{n} q_1(x_i,y_i)}{\partial x^n}, \cdots,\frac{\partial^{n} q_{i-1}(x_i,y_i)}{\partial x^n}, \frac{\partial^{n} q_{i+1}(x_i,y_i)}{\partial x^n}, \cdots \frac{\partial^{n} q_N(x_i,y_i)}{\partial x^n}]^T, $$ $$ [\textbf{Q}]=\left[ \begin{array}{ccc} 1 &\cdots &1 \\ q_1(x_1,y_1) &\cdots & q_1(x_N,y_N)\\ \vdots & \ddots &\vdots \\ q_{i-1}(x_1,y_1) &\cdots & q_{i-1}(x_N,y_N) \\ q_{i+1}(x_1,y_1) &\cdots & q_{i+1}(x_N,y_N) \\ \vdots & \ddots &\vdots \\ q_N(x_1,y_1) &\cdots & q_N(x_N,y_N) \end{array} \right], $$ and $$ \mathbf{W}_i^{(n)_x}=[w_{i,1}^{(n)},w_{i,2}^{(n)},...,w_{i,N}^{(n)}]^T. $$ Clearly, there exists a unique solution only if the collocation matrix $[\textbf{Q}]$ is non-singular. The non-singularity of the collocation matrix $[\textbf{Q}]$ depends on the properties of used RBFs. \cite{micchelli1984interpolation} proved that matrix $[\textbf{Q}]$ is conditionally positive definite for MQ-RBFs. This fact can not guarantee the non-singularity of matrix $[\textbf{Q}]$. \cite{hon2001unsymmetric} showed that cases of singularity are quite rare, and not serious objection to a valuable numerical technique. Therefore, the coefficient vector $\mathbf{W}_i^{(n)_x}$ can be obtained by \begin{equation}\label{eq18} \mathbf{W}_i^{(n)_x}=[\textbf{Q}]^{-1}\textbf{q}_i^{(n)_x}. \end{equation} Then, the coefficient vector $\mathbf{W}_i^{(n)_x}$ can be used to approximate the $n$-th-order derivative in the $x$ direction for any unknown smooth function at node $i$. In a similar manner, all the base functions are substituted into Eq. (\ref{eq9}) 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 \begin{equation}\label{eq20} q_j(x,y)=\sqrt{(x-x_j)^2+(y-y_j)^2+c_j^2}-\sqrt{(x-x_i)^2+(y-y_i)^2+c_i^2}, \end{equation} where second- and higher-order derivatives of $q_j(x,y)$ can also be obtained by differentiating Eq. (\ref{eq20}). 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. %%%%%%%%%%%%%%%%%%% \section{Global RBFs-DQ method for multi-dimensional Black-Scholes equation } \label{sec2} In this section, first we will introduce the variable transformation $S_i=e^{x_i }$ and $C(s_1,...,s_d,t)=V(x_1,...,x_d,t)$. Under this change of variables Eq. (\ref{eq1}) becomes the parabolic PDE with constant coefficients as \begin{equation}\label{eq201} \frac{\partial V}{\partial t}+ \dfrac{1}{2} \sum_{i,j=1}^{d} \rho_{ij} \sigma_i\sigma_j \frac{\partial^2 V}{\partial x_i \partial x_j}+ \sum_{i=1}^{d}(r- \dfrac{1}{2}\sigma_i^2) \frac{\partial V}{\partial x_i}-rV=0, \end{equation} as Eq. (\ref{eq201}) satisfies the conditions of the following theorem, we can remove the crossing terms $$\rho_{ij} \sigma_i\sigma_j \frac{\partial^2 V}{\partial x_i \partial x_j}, i\neq j,$$ then we can apply the RBF-DQ method to option pricing. \begin{thm} \label{teo1} Consider the second-order equation $$ \sum_{i,j=1}^{n}a_{ij}u_{x_ix_j}+\sum_{i=1}^{n}a_{i}u_{x_i}+a_0u=0, $$ 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 $\sum_{i,j=1}^{n}a_{ij}u_{x_ix_j}$ as $\sum_{k=1}^{n}d_ku_{x_kx_k}$, where the coefficients $d_k$ are the eigenvalues of the $n\times n$ matrix $A = (a_{ij})$. \end{thm} Therefore, the multi-dimensional Black-Scholes equation transformed into \begin{equation}\label{eq202} -\dfrac{\partial U}{\partial t}=\sum_{i=1}^{d}\lbrace a_i\dfrac{\partial U}{\partial \xi_i}+b_i\dfrac{\partial^2 U}{\partial \xi_i^2}\rbrace-rU\equiv \mathcal{D}U, \end{equation} The domain is discretized by taking $N$ knots according to the method used in previous section and we would like to discrete the Eq. (\ref{eq202}) with respect to time. We introduce a weight $\theta$ and apply the $\theta$ implicit scheme to the problem as \begin{equation}\label{eq203} -\dfrac{\partial U}{\partial t}\vert_{(t_n, \chi_k)}\cong \dfrac{\tilde{U}_k^n-\tilde{U}_k^{n-1}}{\Delta t}=(1-\theta)\mathcal{D}\tilde{U}\vert_{(t_n, \chi_k)}+\theta\mathcal{D}\tilde{U}\vert_{(t_{n-1}, \chi_k)} \end{equation} Denote $\tilde{U}\vert_{(t_n, \chi_k)}=\tilde{U}_k^n$, where represents the approximation of the function value $U$ at knot $k$, $\chi_k=(\xi_1^k,...,\xi_d^k )$, and the time $t_n$. By applying the RBF-DQ method, Eq. (\ref{eq202}) can be rewritten in a discrete form as follows: \begin{align*} (1+r(1-\theta) \triangle t)\tilde{U}^{n}_k + (1-\theta){\triangle t}\lbrace-\sum_{i=1}^{d}\lbrace a_i\sum_{j=1}^{N}w_{kj}^{(1)_{\xi_i}}\tilde{U}^n_j+ b_i\sum_{j=1}^{N}w_{kj}^{(2)_{\xi_i}}\tilde{U}^n_j\rbrace\rbrace = \\ (1-r\theta \triangle t)\tilde{U}^{n-1}_k -\theta{\triangle t}\lbrace-\sum_{i=1}^{d}\lbrace a_i\sum_{j=1}^{N}w_{kj}^{(1)_{\xi_i}}\tilde{U}^{n-1}_j+ b_i\sum_{j=1}^{N}w_{kj}^{(2)_{\xi_i}}\tilde{U}^{n-1}_j\rbrace\rbrace , \end{align*}\\ where $k=1,2,...,N$ and $w_{kj}^{(1)_{\xi_i}}$ and $w_{kj}^{(1)_{\xi_i}}$ represents the computed weighting coefficients in the DQ approximation. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Stability of the method} \label{sec3} In this section, we study the stability of the implicit finite difference method described above. Let us assume that $\textbf{U}$ be exact and $\tilde{\textbf{U}}$ is the numerical solution of equation (\ref{eq203}). Note that the discrete equations in the last section can be rewritten as the following form \begin{equation}\label{eq27} {[\textbf{I}+(1-\theta) \Delta t \textbf{D}]}\tilde{\textbf{U}}_n={[\textbf{I}-\theta \Delta t \textbf{D}]}\tilde{\textbf{U}}_{n-1}. \end{equation} So we can write Eq. (\ref{eq27}) as \begin{equation}\label{eq28} \tilde{\textbf{U}}_n=[\textbf{I}+(1-\theta) \Delta t \textbf{D}]^{-1}[\textbf{I}-\theta \Delta t \textbf{D}]\tilde{\textbf{U}}_{n-1}=\textbf{E}\tilde{\textbf{U}}_{n-1}, \end{equation} where $$\textbf{D}=r-\sum_{i=1}^{d}\lbrace a_i\lbrace [\textbf{Q}]^{-1}[\textbf{Q}]_{\xi_i^{(1)}}\rbrace^T+b_i\lbrace [\textbf{Q}]^{-1}[\textbf{Q}]_{\xi_i^{(2)}}\rbrace^T\rbrace$$ is $N\times N$ matrix determined by descretized form of Eq. (\ref{eq202}) and $$ [\textbf{Q}]_{\xi_i^{(j)}}=[\textbf{q}^{(j)_{\xi_i}}_1,\textbf{q}^{(j)_{\xi_i}}_2, \cdots, \textbf{q}^{(j)_{\xi_i}}_N], j=1,2.$$ 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 $$ \Vert \tilde{\textbf{U}}_n\Vert\leq C\Vert \tilde{\textbf{U}}_0\Vert,n\rightarrow\infty, \Delta t, \Delta \chi\rightarrow 0, $$ If $\textbf{E}$ depended on $n$ we would get a product of the operators at each time level. Taking a norm $$ \Vert \tilde{\textbf{U}}_n\Vert=\Vert\textbf{E}^n \tilde{\textbf{U}}_0\Vert\leq\Vert\textbf{E}^n\Vert\vert \tilde{\textbf{U}}_0\Vert,n\rightarrow\infty, \Delta t, \Delta \chi\rightarrow 0, $$ Therefore, the numerical scheme is stable if and only if there exists positive constant $C$ such that $$ \Vert\textbf{E}^n\Vert\leq C,n\rightarrow\infty, \Delta t, \Delta \chi\rightarrow 0, $$ Since $\rho(\textbf{E})$ as the spectral radius of $\textbf{E}$ provides a lower bound to any matrix norm, for the scheme to remain stable, we should have $\rho(\textbf{E})\leq 1$ or equivalently we can say that \begin{equation}\label{eq29} \vert\frac{1-\theta\Delta t \eta_{\textbf{D}}}{1+(1-\theta)\Delta t \eta_{\textbf{D}}}\vert\leq1, \end{equation} which holds for $ \eta_{ \textbf{D}}$, eigenvalues of the matrix $\textbf{D}$, located in the right half plane. Inequality (\ref{eq29}) also shows that stability of the scheme, in the case of RBFs with shape parameter like MQ, depends upon shape parameter. %%%%%%%%%%%%%%%%%%%% \section{Numerical computation of two asset European option} \label{sec3} For two-dimensional option pricing under Geometric Brownian motion framework, consider the PDE (\ref{eq1}) with $d=2$ and the final time condition for European call is \begin{equation} C(s_1,s_2,T)=max(max(s_1,s_2)-E,0), \end{equation} where $E$, $T$ are exercise price and maturity time. The following boundary conditions are imposed\ \begin{equation} C(s_{1\max},s_2,T)=s_{1\max}-Ee^{-r(T-t)}, \end{equation} \begin{equation} C(s_1,s_{2\max},T)=s_{2\max}-Ee^{-r(T-t)}, \end{equation} where $s_{1\max}$,\ $s_{2\max}$ are respectively maximum of $s_1$, maximum of $s_2$. In this case, the variable transformation $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. (\ref{eq1}) in which d=2 becomes \begin{equation}\label{eq22} \frac{\partial V}{\partial t}+ \rho \sigma_1 \sigma_2 \frac{\partial^2 V}{\partial x \partial y}+ \dfrac{1}{2} \sigma_1 ^2 \frac{\partial^2 V}{\partial x^2}+ \dfrac{1}{2} \sigma_2^2 \frac{\partial^2 V}{ \partial y^2}+ (r-\dfrac{1}{2}\sigma_1^2) \frac{\partial V}{\partial x}+(r-\dfrac{1}{2}\sigma_2^2) \frac{\partial V}{\partial y}-r V=0.\ \end{equation} Now we can remove the crossing term $\rho \sigma_1 \sigma_2 \frac{\partial^2 V}{\partial x \partial y}$ . In this way, we apply characteristic line to remove mixed derivatives in Eq. (\ref{eq22}) by considering its characteristic equation for the above PDE \cite{Jinghui} as $$\dfrac{dy}{dx}=\dfrac{\rho \sigma_1 \sigma_2 +\sqrt{ \sigma_1^2 \sigma_2^2(\rho^2-1) }}{\sigma_1^2}=\dfrac{\sigma_2}{\sigma_1}(\rho+\sqrt{1-\rho^2}i),$$\\ where $i=\sqrt{-1}.$ Thus we observe the following relation between spatial variables $x$ and $y$,\\ $$\sigma_1y-\sigma_2\rho x-i\sigma_2\sqrt{1-\rho^2}x=c^{\prime} .$$\\ Let $\xi=\sigma_1y-(\sigma_2\rho) x, \eta=-\sigma_2\sqrt{1-\rho^2}x,$ and $V(x,y,t)=U(\xi,\eta,t)$ then \begin{equation} \begin{split} \left\{\!\! \begin{array}{cc} \dfrac{\partial V}{\partial x}=\dfrac{\partial U}{\partial \xi}\dfrac{\partial \xi}{\partial x}+\dfrac{\partial U}{\partial \eta}\dfrac{\partial \eta}{\partial x}=-\sigma_2\rho \dfrac{\partial U}{\partial \xi}-\sigma_2\sqrt{1-\rho^2}\dfrac{\partial U}{\partial \eta} \\\\ \dfrac{\partial V}{\partial y}=\dfrac{\partial U}{\partial \xi}\dfrac{\partial \xi}{\partial y}+\dfrac{\partial U}{\partial \eta}\dfrac{\partial \eta}{\partial y}=\sigma_1\dfrac{\partial U}{\partial \xi}. \end{array}\right.\quad \nonumber \end{split} \end{equation}\\ Further, we have \begin{equation} \begin{split} ~~~ \left\{\!\! \begin{array}{ccc} \dfrac{\partial^2 V}{\partial x^2}=\sigma_2^2\rho^2\dfrac{\partial^2 U}{\partial \xi^2}+2\sigma_2^2\rho\sqrt{1-\rho^2}\dfrac{\partial^2 U}{\partial \xi \partial \eta}+\sigma_2^2(1-\rho^2) \dfrac{\partial^2 U}{\partial \eta^2} \\\\ \dfrac{\partial^2 V}{\partial x\partial y}=-\sigma_1\sigma_2\rho\dfrac{\partial^2 U}{\partial \xi^2}-\sigma_1\sigma_2\sqrt{1-\rho^2}\dfrac{\partial^2 U}{\partial \xi \partial \eta} \\\\ \dfrac{\partial^2 V}{\partial y^2}=\sigma_1^2\dfrac{\partial^2 U}{\partial \xi^2}. \end{array}\right.\quad \nonumber \end{split} \end{equation}\\ Hence, $$ \dfrac{1}{2} \sigma_1 ^2 \frac{\partial^2 V}{\partial x^2}+\rho \sigma_1 \sigma_2 \frac{\partial^2 V}{\partial x \partial y}+ \dfrac{1}{2} \sigma_2^2 \frac{\partial^2 V}{ \partial y^2}=\dfrac{1}{2} \sigma_1 ^2 \sigma_2 ^2(1-\rho^2)(\frac{\partial^2 U}{\partial \xi^2}+\frac{\partial^2 U}{\partial \eta^2}),$$\\ and $$(r-\dfrac{1}{2}\sigma_1^2) \frac{\partial V}{\partial x}+(r-\dfrac{1}{2}\sigma_2^2) \frac{\partial V}{\partial y} =-(r\sigma_2\rho-\dfrac{1}{2}\sigma_1^2\sigma_2\rho-r\sigma_1+\dfrac{1}{2}\sigma_2^2\sigma_1)\dfrac{\partial U}{\partial \xi}-(r\sigma_2\sqrt{1-\rho^2}-\dfrac{1}{2}\sigma_1^2\sigma_2\sqrt{1-\rho^2})\dfrac{\partial U}{\partial \eta}.$$\\ Thus, the Black-Scholes equation in (\ref{eq22}) ransformed into \\ \begin{equation}\label{eq23} \dfrac{\partial U}{\partial t}+a\dfrac{\partial U}{\partial \xi}+b\dfrac{\partial U}{\partial \eta}+c(\dfrac{\partial^2 U}{\partial \xi^2}+\dfrac{\partial^2 U}{\partial \eta^2})-rU=0, \end{equation} where\\ $$a=-(r\sigma_2\rho-\dfrac{1}{2}\sigma_1^2\sigma_2\rho-r\sigma_1+\dfrac{1}{2}\sigma_2^2\sigma_1), $$ $$b=-(r\sigma_2\sqrt{1-\rho^2}-\dfrac{1}{2}\sigma_1^2\sigma_2\sqrt{1-\rho^2}),$$ $$c=\dfrac{1}{2}\sigma_1^2\sigma_2^2(1-\rho^2).$$ where the mixed derivative term vanishes. Eq. (\ref{eq23}) can be rewritten in a discrete form as follows: \begin{align*} (1+r(1-\theta) \triangle t)\tilde{U}^{n}_k +(1-\theta) \triangle t\{-a\sum_{i=1}^{N}w_{ki}^{(1)_{\xi}}\tilde{U}^n_i-b\sum_{i=1}^{N}\tilde{w}_{ki}^{(1)_{\eta}}\tilde{U}^n_i -c(\sum_{i=1}^{N}w_{ik}^{(2)_{\xi}}\tilde{U}^n_i+\sum_{i=1}^{N}\tilde{w}_{ki}^{(2)_{\eta}}\tilde{U}^n_i)\} = \\ (1-r\theta \triangle t)\tilde{U}^{n-1}_k -\theta \triangle t\{-a\sum_{i=1}^{N}w_{ki}^{(1)_{\xi}}\tilde{U}^{n-1}_i-b\sum_{i=1}^{N}\tilde{w}_{ki}^{(1)_{\eta}}\tilde{U}^{n-1}_i-c(\sum_{i=1}^{N}w_{ik}^{(2)_{\xi}}\tilde{U}^{n-1}_i+\sum_{i=1}^{N}\tilde{w}_{ki}^{(2)_{\eta}}\tilde{U}^{n-1}_i)\} \end{align*}\\ where $i=1,2,...,N$. For illustration of the accuracy of the proposed method, we consider a European call option with two underlying assets with $E = 10, r = 0.05, \sigma_1 = 0.22,\sigma_2 = 0.14,\rho=0.5$ and $T = 0.5$ (year) with 100 time steps . Let $s_1 \in [e^{-3.5},e^{4}]$, $s_2 \in [e^{-3.5},e^{4}]$, and $\theta=0.5$. To study the behavior of the method, different structured mesh sizes are used for the point distribution in $\xi$ and $\eta$ direction with shape parameter $c^2=1.25(D/\sqrt{N}) $ as selected by Franke (1982), where $D$ is the diameter of the minimal circle enclosing all supporting points. Fig. 1 shows that our scheme can provide reasonable approximations for Stulz method \cite{stulz1982options} which provides a closed form solution and is available in Matlab. The root-mean-square-error (RMSE) defined by the $$ RMSE=\dfrac{1}{N}\sqrt{\sum_{i=1}^{N}(U_{Stz.}^i-\tilde{U}_{num.}^i)^2}, $$ where $U_{Stz.}$ is the solution by Stulz Method, $\tilde{U}_{num.}$ is solution by numerical approximations. As $\rho(\textbf{E})$, the spectral radius of matrix $\textbf{E}$, is equivalent or less than unity in each simulation, it guarantees the stability of the method. \begin{figure} \begin{center} \includegraphics[width=5 in,height=2 in]{fig1.eps} \caption{Comparison of two-asset option pricing results by RBF-DQ method and Stulz method}\label{fig1} \end{center} \end{figure} \section{Conclusion} \label{sec4} A mesh-free RBFs-DQ method is presented and used for numerical solution of multi-dimensional Black-Scholes equation in this paper. 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. %\bibliographystyle{plain} %\bibliography{mybib} \begin{thebibliography}{99} \bibitem{black1973pricing} F. Black, M. Scholes, {\it The pricing of options and corporate liabilities}. J. of political economy, 637-654, (1973). \bibitem{merton1973theory} R. C. Merton, {\it Theory of rational option pricing}. The Bell J. economics and management science, 141-183, (1973). \bibitem{broadie2004anniversary} M. Broadie , J. B. Detemple, {\it Anniversary article: Option pricing: Valuation models and applications}. Management Science 50(9), 1145-1177, (2004). \bibitem{jo2013comparison} J. Jo , Y. Kim, {\it Comparison of numerical schems on multi-dimensional black-scholes equation}. Bulletin of the Korean mathematical society 50(6), 2035-2051, (2013). \bibitem{glasserman2004monte}P. Glasserman, {\it Monte Carlo methods in financial engineering}. Springer vol 53, (2004). \bibitem{ibanez2004monte}A. Ibanez , F. Zapatero, {\it Monte carlo valuation of american options through computation of the optimal exercise frontier}. J. Financial and Quantitative Analysis 39(02), 253-275, (2004). \bibitem{choi2001numerical} S. Choi , M. D. Marcozzi, {\it A numerical approach to american currency option valuation}. J. Derivatives 9(2), 19-29, (2001). \bibitem{choi2003valuation}S. Choi , M. Marcozzi, {\it The valuation of foreign currency options under stochastic interest rates}. Computers \& Mathematics with Applications 46(5), 741-749, (2003). \bibitem{clarke1999multigrid} N. Clarke , K. Parrott, {\it Multigrid for american option pricing with stochastic volatility}. Applied Mathematical Finance 6(3), 177-195, (1999). \bibitem{forsyth1999finite} P. Forsyth , K. Vetzal , R. Zvan, {\it A finite element approach to the pricing of discrete lookbacks with stochastic volatility}. Applied Mathematical Finance 6(2), 87-106, (1999). \bibitem{wilmott1993option} P. Wilmott , J. Dewynne , S. Howison, {\it Option pricing: mathematical models and computation}. Oxford financial press, (1993) \bibitem{zhu1999singularity} Y. I. Zhu , Y. Sun, {\it The singularity-separating method for two-factor convertible bonds}. J. Computational Finance 3(1), 91-110, (1999). \bibitem{nielsen2008penalty} B. F. Nielsen , O. Skavhaug , A. Tveito, {\it Penalty methods for the numerical solution of american multi-asset option problems}. J. Computational and Applied Mathematics 222(1), 3-16, (2008). \bibitem{boyle1989numerical}P. P. Boyle , J. Evnine , S. Gibbs, {\it Numerical evaluation of multivariate contingent claims}. Review of Financial Studies 2(2), 241-250, (1989). \bibitem{voss2003fourth}D. A Voss , A. Q. M. Khaliq , S. Kazmi , H. He, {\it A fourth order l-stable method for the black-scholes model with barrier options}. In Computational Science and Its ApplicationsICCSA 2003, Springer, 199-207, (2003). \bibitem{khaliq2006linearly} A. Khaliq , D. Voss , S. Kazmi, {\it A linearly implicit predictor-corrector scheme for pricing american options using a penalty method approach}. J. Banking \& Finance 30(2), 489-502, (2006). \bibitem{fasshauer2004using} G. E. Fasshauer , A. Q. M. Khaliq , D. A. Voss, {\it Using meshfree approximation for multiasset american options}. J. the Chinese Institute of Engineers 27(4), 563-571, (2004). \bibitem{martin2014stabilized}J. Martin-Vaquero , A. Khaliq , B. Kleefeld, {\it Stabilized explicit runge-kutta methods for multi-asset american options}. Computers \& Mathematics with Applications 67(6), 1293-1308, (2014). \bibitem{shu2004solution} C. Shu , H. Ding , K. Yeo, {\it Solution of partial differential equations by a global radial basis function-based differential quadrature method}. Engineering analysis with boundary elements 28(10), 1217-1226, (2004). \bibitem{bellman1971differential} R. Bellman, J. Casti, {\it Differential quadrature and long-term integration}. J. Mathematical Analysis and Applications 34(2), 235-238, (1971). \bibitem{bellman1972differential}R. Bellman , B. Kashef , J. Casti, {\it Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations}. J. Computational Physics 10(1), 40-52, (1972). \bibitem{shu1992application}C. Shu , B. E. Richards, {\it Application of generalized differential quadrature to solve two-dimensional incompressible navier-stokes equations}. International J. Numerical Methods in Fluids 15(7), 791-798, (1992). \bibitem{shu1997fourier} C. Shu , Y. Chew, {\it Fourier expansion-based differential quadrature and its application to helmholtz eigenvalue problems}. Communications in Numerical Methods in Engineering 13(8), 643?653, (1997). \bibitem{shu2000differential} C. Shu, {\it Differential quadrature and its application in engineering}. Springer (2000). \bibitem{bert1988two} C. W. Bert , S. K. Jang , A. G. Striz, {\it Two new approximate methods for analyzing free vibration of structural components}. AIAA J. 26(5), 612-618, (1988). \bibitem{sherbourne1991differential} A. Sherbourne ,M. Pandey , {\it Differential quadrature method in the buckling analysis of beams and composite plates}. Computers \& Structures 40(4), 903-913, (1991). \bibitem{gutierrez1994method} R. Gutierrez, P. Laura , R. Rossi, {\it The method of differential quadrature and its application to the approximate solution of ocean engineering problems}. Ocean Engineering 21(1), 57-66, (1994). \bibitem{bert1996differential}C.W. Bert , M. Malik, {\it The differential quadrature method for irregular domains and application to plate vibration}. International J. Mechanical Sciences 38(6), 589-606, (1996). \bibitem{han1997eight}J. B. Han , K. Liew, {\it An eight-node curvilinear differential quadrature formulation for reissner/mindlin plates}. Computer Methods in Applied Mechanics and Engineering 141(3), 265-280, (1997). \bibitem{shu2000free} C. Shu, W. Chen , H. Du, {\it Free vibration analysis of curvilinear quadrilateral plates by the differential quadrature method}. J. Computational Physics 163(2), 452-466, (2000). \bibitem{franke1982scattered} R. Franke, {\it Scattered data interpolation: Tests of some methods}. Mathematics of computation 38(157), 181-200, (1982). \bibitem{hardy1971multiquadric} R. L. Hardy, {\it Multiquadric equations of topography and other irregular surfaces}. J. geophysical research 76(8), 1905-1915, (1971). \bibitem{micchelli1984interpolation} C. A. Micchelli, {\it Interpolation of scattered data: distance matrices and conditionally positive definite functions}. Springer (1984). \bibitem{hon2001unsymmetric} Y. Hon , R. Schaback, {\it On unsymmetric collocation by radial basis functions}. Applied Mathematics and Computation 119(2), 177-186, (2001). \bibitem{Jinghui} J. Zhou , {\it Multi-asset option pricing with Levy process}. Ph. D. Thesis, National University of Singapore, (2009). \bibitem{stulz1982options} R. Stulz, {\it Options on the minimum or the maximum of two risky assets: analysis and applications}. J. Financial Economics 10(2), 161-185, (1982). \end{thebibliography} \end{document}