The Fractional Step Method versus the Radial Basis Functions for Option Pricing with Correlated Stochastic Processes

: In option pricing models with correlated stochastic processes, an option premium is commonly a solution to a partial di ﬀ erential equation (PDE) with mixed derivatives in more than two space dimensions. Alternating direction implicit (ADI) ﬁnite di ﬀ erence methods are popular for solving a PDE with more than two space dimensions; however, it is not straightforward to employ the ADI method for solving a PDE with mixed derivatives. The aim of this study is to ﬁnd out which numerical method would be appropriate to solve PDEs with mixed derivatives based on the accuracy of the solutions and the computation time. This study applies the fractional step method and the radial basis functions to solve a PDE with a mixed derivative, and investigates the e ﬃ ciency of these numerical methods. Numerical experiments are conducted by applying these methods to exchange option pricing; exchange options are selected because the exchange option premium has an analytical form. The numerical results show that the both methods calculate premiums with high accuracy in the presence of mixed derivatives. The fractional step method calculates the option premium more accurately and much faster than the radial basis functions. Therefore, from the numerical experiments, this study concludes that the fractional step method is more appropriate than the radial basis functions for solving a PDE with a mixed derivative.


Introduction
Many financial option models are formulated using two or more stochastic processes. Margrabe (1978) derived a valuation formula for exchange option, which is a function of two stochastic asset prices. Heston (1993) developed an equity option pricing model in which the logarithm of an equity price follows a Wiener process with stochastic volatility. Amin and Jarrow (1992) formulated an equity option pricing model with a stochastic equity price and the stochastic risk-free interest rate. Heath et al. (1992) proposed a multi-factor risk-free interest rate model. Lando (1998) evaluated credit derivatives with a stochastic default intensity and the stochastic risk-free interest rate. In these models, option premiums are given as solutions to partial differential equations (PDEs) with mixed derivatives in two or more space dimensions.
PDEs seldom have analytical solutions and hence, numerical procedures are essential to calculate option premiums. Various numerical methods have been applied to option pricing, including tree models (Shreve 2004), finite difference methods (Duffy 2006), and Monte Carlo simulations (Jäckel 2002), among others. When option pricing models are formulated using two or more stochastic processes, 2 of 13 correlations between the stochastic processes affect option premiums. However, numerical procedures for correlated stochastic processes are not straightforward. The construction of a correlated tree is not a simple task. Correlated stochastic processes introduce mixed derivatives to PDEs. Among the finite difference methods, the alternating direction implicit (ADI) finite difference method is popular in finance studies; however, the original ADI method is unable to solve PDEs with mixed derivatives (Duffy 2006, p. 224). Monte Carlo simulation can generate correlated random variables; however, it is a time-consuming way to obtain accurate option premiums. Thus, it is essential to establish efficient methods to solve PDEs with mixed derivatives.
The purpose of this study is to examine the accuracy of the solutions and the computation time for solving PDEs with mixed derivatives. In this study, the fractional step method and the radial basis function are used as the basis of the numerical methods. The fractional step method is unfamiliar in the field of finance research, which usually employs the ADI finite difference method to solve PDEs with two or more space dimensions. As noted above, it is not a straightforward process to use the ADI finite difference method to solve PDEs with mixed derivatives. 1 The radial basis functions, also known as the meshfree method, were not used until 2000; since then, however, they have become very popular among researchers in many fields (see, e.g., Fasshauer 2007;Liu 2009;Liu and Gu 2005). The radial basis functions are easy to implement, even in higher space dimensions. Thus, the fractional step method and the radial basis functions are applied to solve PEDs with mixed derivatives. Both numerical methods are conducted to evaluate the exchange option premium; the exchange option premium is a solution to a PDE with mixed derivatives, and an analytical solution for the exchange option is known. Therefore, exchange options are suitable as a benchmark to evaluate the efficiency of the numerical methods.
The rest of the paper is organized as follows. Section 2 explains payoffs of exchange options and presents the analytical pricing formula for exchange option. The fractional step method and the radial basis functions are reviewed to solve a PDE with mixed derivatives in Sections 3 and 4, respectively. The drawback of the ADI method is also explained in Section 3. Section 5 presents the results from numerical experiments, which examine the accuracy and computation time of the fractional step method and the radial basis functions. Both numerical methods are applied to pricing of exchange options. Section 5.1 describes the general setting for the numerical experiments. Section 5.2 investigates the sufficient number of partitions of the space dimensions to obtain the accurate numerical results. The calculation times are examined in this subsection. Section 5.3 investigates the accuracy of both numerical methods for various correlations. Section 6 concludes.

Exchange Option
An exchange option allows its holder to exchange one asset for another. Consider an exchange option to exchange stock 2 for stock 1 at maturity, T. Payoff at maturity, f T , is expressed as The exchange option premium is calculated using the Black-Scholes model (Margrabe 1978). In the Black-Scholes framework, the time evolution of a stock price follows a log-normal distribution, 1 Several studies employ the ADI finite difference method to evaluate PDEs with mixed derivatives (Hendricks et al. 2016;In't Hout and Foulon 2010;Jeong and Kim 2013). where r is the risk-free rate, σ i is its volatility, and W i,t is a Wiener process under a martingale measure (i = 1, 2). The correlation between Wiener processes is ρ 12 , where , stands for a quadratic variation. Stochastic calculus and the arbitrage condition show that the time-t exchange option premium, f , satisfies the following PDE, is a logarithm transformed price. A pricing formula is developed by Margrabe (1978) as, Sections 3 and 4 explain the fractional step method and the radial basis functions, respectively. In these sections, the PDE in two space dimensions is written as which includes the PDE, as can be seen in Equation (4). The terminal time condition is given as intrinsic values of the exchange option. The boundary conditions for space variables are set following Jeong and Kim (2013), where x i,L and x i,U are the lower and upper bounds of space variables, respectively (i = 1, 2).

The Fractional Step Method
The fractional step method is a variation of the finite difference methods. The fractional step method was developed by Yanenko and his collaborators (Yanenko 1971). As far as we are concerned, few previous studies applied the fractional step method to option pricing: American option pricing (Ikonen and Toivane 2007) and barrier option pricing (Itkin and Carr 2008;Toivanen 2010).
First, the finite difference method is reviewed to explain the fractional step method. Again, consider a PDE, x ≤ x U , and y L ≤ y ≤ y U , with an initial condition, and boundary conditions, lim y→y L ∂ 2 u ∂y 2 (t, x, y) = 0, lim y→y U ∂ 2 u ∂y 2 (t, x, y) = 0 To find a solution to Equation (8) using a numerical procedure, the time axis is discretized as t n = n(∆t), n = 0, 1, . . . , N, ∆t = T N .
Likewise, the two-dimensional space is discretized to construct a mesh, The numbers of partitions, N and M, are called the number of mesh points. In the following, u(x i , y j , t n ) is denoted as u n i,j . Here, the ADI method is discussed. When the mixed derivative is absent, the PDE reduces to where The ADI method inserts an intermediate time t n+ 1 2 between t n and t n+1 . Then, from t n to t n+ 1 2 , the implicit method is applied to x and the explicit method is applied to y. Next, from t n+ 1 2 to t n+1 , the implicit method is applied to y and the explicit method is applied to x. These procedures are rewritten as To apply the ADI method, the derivative terms of x and of y should be decoupled. Therefore, the ADI method is not capable of solving a PDE with mixed derivatives. The fractional step method was developed by Yanenko (1971). This method is also known as the splitting method or the locally one-dimensional method. The fractional step method is applied to solve PDEs with mixed derivatives. By using the fractional step method with a θ scheme, the PDE of Equation (8) is written as The fully implicit scheme corresponds to θ = 1 and the Crank-Nicolson scheme is obtained by setting θ = 1 2 . Expressing the differential terms in Equation (19) by finite differences, it can be rewritten as where Similarly, Equation (20) is rewritten as where Equations (21) and (29) are successively solved from time t 1 to t N , given an initial condition.

The Radial Basis Functions
The radial basis functions, also known as the meshfree method, is a relatively new numerical method (Fasshauer 2007;Liu 2009;Liu and Gu 2005), and the first application of the radial basis functions to option pricing was performed by Hon and Mao (1998). The radial basis functions have been applied to option pricing with the Levy process model (Brummelhuis and Chan 2014;Chan 2016;Golbabai et al. 2012;Kadalbajoo et al. 2016;Lung and Hubbert 2014;Saib et al. 2012), the Heston model (Ballestra and Pacelli 2013;Mollapouras et al. 2019), and the Cox-Ingersoll-Ross model (Jebreen 2019). The radial basis functions can be applied to various types of options other than plain vanilla options. For example, the radial basis functions can be applied to the barrier option (Goto et al. 2007), the basket option (Company et al. 2018;Shcherbakov and Larsson 2016), and the Asian option (Goto et al. 2007).
Given a function Here, the basis functions, φ j (x) , are dependent on the distances between x and ξ j (ξ j ∈ R d ); thus, they are called the radial basis functions. The set of points {x 1 , . . . , x n } and {ξ 1 , . . . , ξ n } are called the collocation points and centers, respectively. Formally, the collocation points may be different from the centers; however, they coincide with each other in many cases. There are many candidates for the radial basis functions. In this study, following Hon and Mao (1998) and Brummelhuis and Chan (2014), we chose Hardy's multiquadratics, The radial basis functions were applied to PDEs in two space dimensions. In the following, it is convenient to denote (x, y) in place of x(∈ R 2 ) and redefine t as a term to maturity (t := T − t). Consider a PDE in two-space dimensions, Substituting Equation (37) into Equation (39), the PDE is written as Defining the following vectors and matrices, k, l ∈ x, y .
Equation (37) is written as y = Φu (47) and the matrix representation of the PDE of Equation (40) is expressed as, where Then the next step is to discretize the time differential using the θ scheme, Arranging Equation (50), we obtain the recursive relation Similar to the case of the fractional step method, the fully implicit scheme corresponds to θ = 1 and the Crank-Nicolson scheme is obtained by setting θ = 1 2 . Successive iterations of solving u(t i+1 ) using u(t i ) from i = 1 to i = N give the option premium at the initial time. Note that the boundary conditions for space dimensions are not necessary to the radial basis functions.
The radial basis functions do not ensure the positivity of the time-t i option premium at x k , To remedy the negative premium, the following procedure is conducted at each time step: (i) once u j (t i ) is obtained, f = Φu is calculated, (ii) negative values of f (t i , x k ) are replaced by zero, and (iii) redefine u j (t i ) as u = Φ −1 f .

Numerical Experiments
This section examines the accuracy and computation time of the fractional step method and the radial basis function. Both numerical methods are applied to pricing of exchange options. Section 5.1 describes the general setting for the numerical experiments. Section 5.2 investigates the sufficient number of partitions of the space dimensions to obtain the accurate numerical results, and examines the calculation time. Section 5.3 investigates the accuracy of both numerical methods for various correlations.

General Setting
Numerical experiments examine pricing errors and computation times using the fractional step method and the radial basis functions to calculate exchange option premiums. Both methods are applied with the fully implicit scheme (θ = 1) and the Crank-Nicolson scheme (θ = 1 2 ). The domain is set common to the fractional step method and the radial basis functions. In the application of the radial basis functions, the collocation points and the centers are the same, and are set to the mesh points for the fractional step method. Thus, the number and location of mesh points are the same in the fractional step methods and the radial basis functions, and, therefore, the accuracy and computation time of the fractional step method and the radial basis functions are comparable.
The specification of the exchange option is given as follows unless otherwise specified. The initial prices of stock 1 and 2 are 10 dollars, S 1 = S 2 = 10, and their volatilities are 0.3, σ 1 = σ 2 = 0.3. The risk-free interest rate is 0.05, r = 0.05. To analyze the effect of positive and negative correlations on the exchange option premium, two cases of correlation coefficient, ρ = −0.75 and ρ = 0.75, are studied.
The maturity of the exchange option is 1 year, T = 1. The domain of x is set to [x L , x U ] = [log S L , log S U ] so that the time-T stock prices evolve within ±5σ bounds: The domain of y is taken as the same as that of x, [y L , y U ] = [log S L , log S U ]. The number of partitions of time is set to N = 30. The number of partitions of space is set to M = 80, which produces percentage pricing errors for the exchange option less than 1% both for the radial basis functions and the fractional step methods with the fully implicit and the Crank-Nicolson schemes. The radial basis functions have some parameters to be fixed. Following Hon and Mao (1998) and Brummelhuis and Chan (2014), we set Hardy's multiquadratics to c = 4d min where d min is the minimum distance between any two collocation points x j . As can be shown in the following subsections, this choice gives accurate numerical results. The analytical expression of the exchange option premium, Equation (5), gives theoretical values of 2.2100 and 0.8447 for ρ = −0.75 and ρ = 0.75, respectively. The premium for the negative correlation is greater than that for the positive correlation because the negative correlation produces a wider price difference between stocks 1 and 2. Pricing errors are calculated by subtracting the theoretical value from numerical results. Percentage pricing errors are calculated by dividing the pricing errors by the theoretical value.
The fractional step method and the radial basis functions are implemented using Matlab to solve the PDE with mixed derivatives, as can be seen in Equation (4). Exchange option premiums are calculated under various parameter settings. Numerical calculations are performed on a laptop PC with Intel Xero E-2176M CPU @2.70 GHz and 64 GB RAM.

Number of Partitions of Space Dimension
The minimum number of partitions of the space dimensions is investigated to obtain accurate numerical results. Pricing errors, percentage pricing errors, and computation times are examined by changing the number of partitions of the space dimension from M = 10 to M = 100 with 1 increment. The decimal and percentage pricing errors are shown in Figures 1 and 2, respectively. The computation times are shown in Figure 3. applying the radial basis functions with the fully implicit and the Crank-Nicolson scheme, respectively. The minimum numbers of partition size to achieve percentage pricing errors less than 1% are 66 and 78 for 0.75   applying the radial basis functions with the fully implicit scheme and the Crank-Nicolson scheme, respectively. In conclusion, the fractional step method is more accurate than the radial basis functions when the same mesh is used.   The results in this subsection confirm that both the fractional step method and the radial basis functions can appropriately solve PDEs with mixed derivative terms. Comparing their results, the fractional step method is more accurate than the radial basis functions with the same mesh and the fractional step method is much faster than the radial basis functions.

Errors with Respect to Correlation
The main purpose of this study is to examine whether the fractional step method and the radial basis functions can solve PDEs with mixed derivatives. Numerical results of exchange option premiums are obtained by changing  from 1  to 0.99 with 0.01 increments, and the results are shown in Figure 4. The exchange option premium is a decreasing function with respect to the correlation. The analytical premium varies from 2.3582 for 1   to 0.1692 for 0.99   . Figure 4. Theoretical exchange option premium with respect to  . First, the results on decimal pricing errors are analyzed. The decimal pricing errors are shown in Figure 1. The pricing errors of the fractional step methods monotonically converge to zero as the number of partitions increases. The pricing errors of the radial basis functions converge to zero with fluctuations as the number of partitions increases; when the number of partitions is smaller than 50, the exchange option premium is smaller than the theoretical value for the even number of partitions, while the exchange option premium is bigger than the theoretical value for the odd number of partitions. Thus, an average of numerical value for the partition size m and m + 1 reduces pricing error. Comparing the decimal pricing errors, the fractional step method provides smaller decimal errors than the radial basis functions for a common mesh. The pricing errors for ρ = 0.75 are greater than those for ρ = −0.75 when there is the same number of partitions. Option premiums are smaller for ρ = 0.75 than for ρ = −0.75, which results in larger errors for ρ = 0.75 than for ρ = −0.75.
Next, the results on percentage pricing errors are analyzed. The percentage pricing errors, shown in Figure 2, display similar behavior to the decimal pricing errors. The percentage pricing errors of the fractional step methods monotonically converge to zero. The percentage pricing errors of the radial basis functions converge to zero with fluctuations as the number of partitions increases. The minimum numbers of partition to achieve percentage pricing errors less than 1% are 22 and 21 for ρ = −0.75 applying the fractional step methods with the fully implicit scheme and the Crank-Nicolson scheme, respectively. The minimum numbers of partition size to achieve percentage pricing errors less than 1% are 45 and 54 for ρ = 0.75 applying the fractional step methods with the fully implicit scheme and the Crank-Nicolson scheme, respectively. The minimum numbers of partition size to achieve percentage pricing errors less than 1% are 28 and 34 for ρ = −0.75 applying the radial basis functions with the fully implicit and the Crank-Nicolson scheme, respectively. The minimum numbers of partition size to achieve percentage pricing errors less than 1% are 66 and 78 for ρ = 0.75 applying the radial basis functions with the fully implicit scheme and the Crank-Nicolson scheme, respectively. In conclusion, the fractional step method is more accurate than the radial basis functions when the same mesh is used.
Computation times of the fractional step method and the radial basis functions are shown in Figure 3. Here, the computation times are calculated as an average of computation times for ρ = −0.75 and ρ = 0.75. The computation time of the fractional step method increases linearly as the number of partitions increases, while that of the radial basis functions grows exponentially as the number of partitions increases. Thus, when the number of partitions is 10, the computation times of the fractional step method are comparable with those of the radial basis functions; the fractional step method is much faster than the radial basis functions as the number of partitions increases. As mentioned at the end of Section 4, when the radial basis functions are conducted to evaluate exchange options, the non-negativity of option premium f = Φu is not ensured. To recover the non-negativity of option premium, negative values of option premium at the intermediate times are replaced by zero and re-defined u as u = Φ −1 f . This procedure takes time and increases the computation time of the radial basis functions as the number of partitions increases.
The results in this subsection confirm that both the fractional step method and the radial basis functions can appropriately solve PDEs with mixed derivative terms. Comparing their results, the fractional step method is more accurate than the radial basis functions with the same mesh and the fractional step method is much faster than the radial basis functions.

Errors with Respect to Correlation
The main purpose of this study is to examine whether the fractional step method and the radial basis functions can solve PDEs with mixed derivatives. Numerical results of exchange option premiums are obtained by changing ρ from −1 to 0.99 with 0.01 increments, and the results are shown in Figure 4. The exchange option premium is a decreasing function with respect to the correlation. The analytical premium varies from 2.3582 for ρ = −1 to 0.1692 for ρ = 0.99. exchange options, the non-negativity of option premium fu  is not ensured. To recover the non-negativity of option premium, negative values of option premium at the intermediate times are  The results in this subsection confirm that both the fractional step method and the radial basis functions can appropriately solve PDEs with mixed derivative terms. Comparing their results, the fractional step method is more accurate than the radial basis functions with the same mesh and the fractional step method is much faster than the radial basis functions.

Errors with Respect to Correlation
The main purpose of this study is to examine whether the fractional step method and the radial basis functions can solve PDEs with mixed derivatives. Numerical results of exchange option premiums are obtained by changing  from 1  to 0.99 with 0.01 increments, and the results are shown in Figure 4. The exchange option premium is a decreasing function with respect to the correlation.   First, the decimal pricing errors are examined with respect to the correlations, as depicted in Figure 5. The decimal pricing errors of the fractional step method using the fully implicit scheme are less than 0.01 for −1 ≤ ρ ≤ 0.88, and those using the Crank-Nicolson scheme are less than 0.01 for −1 ≤ ρ ≤ 0.85. The decimal pricing errors of the radial basis functions using the fully implicit scheme are less than 0.01 for −1 ≤ ρ ≤ 0.98, and those using the Crank-Nicolson scheme are less than 0.01 for −1 ≤ ρ ≤ 0.97. Thus, the differences between the fully implicit scheme and the Crank-Nicolson scheme are negligible for the fractional step method and for the radial basis functions. For strong correlations, the decimal and percentage pricing errors of the fractional step method are greater than those of the radial basis functions. The decimal pricing errors of the fractional step method at ρ = 0.99 are 0.138251 and 0.142960 for the fully implicit scheme and the Crank-Nicolson scheme, respectively. The decimal pricing errors of the radial basis functions at ρ = 0.99 are 0.01183 and 0.013793 for the fully implicit scheme and the Crank-Nicolson scheme, respectively. Thus, the radial basis functions give accurate numbers, even for an extremely strong correlation.
Crank-Nicolson scheme are negligible for the fractional step method and for the radial basis functions. For strong correlations, the decimal and percentage pricing errors of the fractional step method are greater than those of the radial basis functions. The decimal pricing errors of the fractional step method at 0.99   are 0.138251 and 0.142960 for the fully implicit scheme and the Crank-Nicolson scheme, respectively. The decimal pricing errors of the radial basis functions at 0.99   are 0.01183 and 0.013793 for the fully implicit scheme and the Crank-Nicolson scheme, respectively. Thus, the radial basis functions give accurate numbers, even for an extremely strong correlation. Next, the percentage pricing errors are examined with respect to the correlations, shown in Figure 5. The behavior of the percentage pricing errors is similar to that of the decimal pricing errors. The percentage pricing errors of the fractional step method using the fully implicit scheme are less   Next, the percentage pricing errors are examined with respect to the correlations, shown in Figure 5. The behavior of the percentage pricing errors is similar to that of the decimal pricing errors. The percentage pricing errors of the fractional step method using the fully implicit scheme are less than 1% for −1 ≤ ρ ≤ 0.86, and those using the Crank-Nicolson scheme are less than 1% for −1 ≤ ρ ≤ 0.82. The percentage pricing errors of the radial basis functions using the fully implicit scheme are less than 1% for −1 ≤ ρ ≤ 0.95, and those using the Crank-Nicolson scheme are less than 1% for −1 ≤ ρ ≤ 0.91. Thus, the fractional step method gives sufficiently accurate results except for extremely strong correlations. The radial basis functions give more accurate results than the fractional step method with the same mesh.
The pricing errors of the fractional step method are reduced by increasing the number of partitions for time and space dimensions to N = 120 and M = 250, respectively. The decimal and percentage pricing errors are depicted in Figure 6. Using this mesh, the percentage errors of the fractional step method using the fully implicit scheme are less than 1% for −1 ≤ ρ ≤ 0.94, and those using the Crank-Nicolson scheme are less than 1% for −1 ≤ ρ ≤ 0.93. Despite the large number of partitions, the computation time is less than the radial basis functions with N = 30 and M = 80.

Conclusions
This study conducts numerical experiments to solve PDEs with mixed derivatives. The experiments show that both the fractional step method and the radial basis functions are efficient methods to solve PDEs with mixed derivatives. The computation time of the radial basis functions is longer than that of the fractional step method when the numbers of partitions are the same. The computation time of the radial basis functions increases exponentially as the number of partitions increases. The radial basis functions do not ensure the non-negativity of the option premium, and it

Conclusions
This study conducts numerical experiments to solve PDEs with mixed derivatives. The experiments show that both the fractional step method and the radial basis functions are efficient methods to solve PDEs with mixed derivatives. The computation time of the radial basis functions is longer than that of the fractional step method when the numbers of partitions are the same. The computation time of the radial basis functions increases exponentially as the number of partitions increases. The radial basis functions do not ensure the non-negativity of the option premium, and it takes time to recover the non-negativity of the option premium at the intermediate times. The computation time of the fractional step method increases linearly as the number of partitions increases. Moreover, the fractional step method gives accurate results even with relatively small numbers of partitions. Both the fractional step method and the radial basis functions give accurate numbers except for the extremely strong correlations. Even under the extremely strong correlations, the fractional step method gives accurate results by increasing the number of partitions with computation time less than the radial basis functions. Therefore, this study concludes that the fractional step method is superior to the radial basis function.
Both numerical methods are applicable to various option pricings in two space dimensions such as equity options under stochastic interest rates, multi-factor interest rate models, credit derivatives under stochastic interest rates and default intensity, and so forth. Moreover, these methods can be applied to a PDE with more than two dimensions. We leave these problems for further research.
Funding: This work was supported by JSPS KAKENHI Grant Numbers 18K01707.