NBIn the academic literature there are two common approaches for the evaluation of financial options. These are stochastic calculus and partial differential equations. The former is the method of choice for statisticians and theoreticians, while the latter is the principal tool of physicists and computer scientists because it lends itself to practical implementation schemes. Occasionally small modifications such as linear regression and binomial trees are used, but these are usually treated within either of the two previously mentioned fields. Rarely do the practitioners of these fields compare and contrast methodologies, let alone admit completely different approaches. While Radial Basis Function (RBF) methodology has previously been applied to solving some differential equations, there are very few papers considering its applicability to financial mathematics. The purpose of this article is to show not only that RBF can solve many of the evaluation problems for financial options, but that with Mathematica it can do so with accuracy and speed.
Introduction
A radial basis function (RBF) is a function 
that depends only on the distance 
between 
and a fixed point 
 |
(1) |
Each function 
is radially symmetric about the center 
. Since their discovery in the early 1970s by Hardy [1, 2], RBFs have become a primary tool for the interpolation of multidimensional scattered data. In the 1990s, Kansa [3, 4] showed that RBF methods are applicable to solving elliptic parabolic and some hyperbolic partial differential equations (PDEs). While this approach has some similarities with finite difference (FD) formulas, there are significant differences in that FD stencils typically extend over a subset of the data points at which derivative approximations are sought, and FD formulas are obtained by differentiating polynomial interpolants rather than RBF interpolants. But the fact that RBF methods allow the solution of parabolic PDEs means that they are applicable to the Black-Scholes PDE and hence to the evaluation of financial options.
The purpose of this article is to apply global RBFs within Mathematica as a spatial collocation scheme for solving European and American option pricing models by extending and implementing the work in this area that has recently been done by Hon and Mao [5] and Fasshauer, Khaliq, and Voss [6]. While some results have been presented, none of the papers mentioned actually exhibit any programs or discuss the programming difficulties that are inherent in RBF models of parabolic PDEs. A number of Runge-Kutta time integration schemes are adopted for the time derivatives of the option model. We show that these schemes result in highly accurate approximations when compared with existing numerical techniques and are inherently more stable than the more commonly used finite element methods.
Radial Basis Function Methodology for the Black-Scholes Partial Differential Equation
It has been shown by Black and Scholes [7] that, assuming the underlying asset price is risk-neutral, all European options satisfy a lognormal parabolic PDE, now called the Black-Scholes PDE. If we let 
be the value of the option with underlying asset 
and elapsed time 
, risk-free interest rate 
, dividend yield 
, and volatility or annualized standard deviation of the asset price 
, then the equation governing all European options is:
 |
(2) |
Like all other PDEs, the specification of the boundary conditions determines the type of option studied. The initial condition for the backward PDE has the following maturity payoff function, with 
being the strike and 
the time of maturity:
 |
(3) |
It is shown in [7] for European options that with 
as the cumulative distribution function of the standard normal distribution, the explicit solution is given by
 |
(4) |
By taking the mathematical derivatives 
of the European options in equation (4), we readily obtain the hedging greek deltas
 |
(5) |
Let us make a transformation into the log space so that 
, where 
. Then observing that 
and furthermore that 
, equations (2) and (3) now become
 |
(6) |
with initial condition
 |
(7) |
The RBF methodology is to represent the option function 
as a linear combination of RBFs 
at the collocation points 
 |
(8) |
There are a number of different choices for the RBFs; the most common include Hardy’s multiquadratic (MQ) 
and the Gaussian 
, where 
is the shape parameter. Extensive research in [3, 4] and Goldberg et al. [8] has demonstrated that the MQ interpolation is superior when solving inhomogeneous PDEs, as is the case with the Black-Scholes equation. The specification of 
is the basis of ongoing research, but here we adopt the recommendation in [1] that 
for collocation points 
. Substituting equation (8) into (6) we now have a system of 
linear equations, 
:
 |
(9) |
which can be further reduced by observing that
 |
(10) |
where the partial derivatives of the MQ RBFs, initially described in equation (1) and in the preceding paragraph with the shape parameter 
, can be easily determined in Mathematica.
Now define the following elements of the 
matrices 
, 
, and 
with the 
-vector 
 |
(11) |
Substituting equations (10), (11), and the previous differentiation output into equation (9), while observing that the elements of equation (10) actually define matrix products, results in the following matrix equation:
 |
(12) |
Solving equation (12) for 
allows the estimation of 
and hence of the original option price 
. It has been shown by Powell [9] that the matrix 
is invertible. This allows us to rewrite equation (12) as a matrix differential equation for 
 |
(13) |
is the following 
matrix and since all of its components, such as 
and the identity matrix 
, are known it is readily computed:
 |
(14) |
There are two approaches to solving equation (13): analytical and a backward time integration scheme. Using the work in [9], it follows that, if 
is invertible and the eigenvectors of 
are independent, then equation (12) has a solution for 
in terms of its eigenvectors 
and eigenvalues 
. Using Mathematica‘s DSolve operator we can explicitly analyze possible solutions to equation (13); however, the expression is large.
Choosing larger values for 
, the dimension of the matrix 
suggests the following general solution:
 |
(15) |
Then taking the time derivative of both sides and recalling that 
because of the definition of 
as the eigensystem of 
, we get:
 |
(16) |
This result satisfies equation (13) and shows that (15) yields an analytical solution for 
. Further, define the vectors 
, 
, and 
that consist of the components 
, 
, and 
as well as the matrix 
that has as columns the eigenvectors 
, 
. Now observe that collocating equation (8) about the central points 
converts it into the explicit matrix equation
 |
(17) |
It can be seen from equation (17) that 
, and hence 
, can be calculated once 
is determined. Since the result holds for all values of time 
, then it also holds for the initial condition expressed by equation (7) at maturity time 
. Putting 
in the above result yields
 |
(18) |
While this result is sufficient for an explicit solution that will cover European options, it is not capable of calculating path-dependent options that have to be updated in a time-dependent manner. In fact, if we choose to again use the DSolve operator for larger values of the dimension of the matrix 
, then we will get very complicated polynomial expressions that have to be solved by the RootSum function. Here is an example; again, the expression is large.
For the purpose of calculating American-style options, we will consider times 
and define 
, 
. Now we use backward-difference time integration schemes that can be any of the explicit first-order (BD1):
 |
(19) |
the explicit second-order backward-difference time integration scheme (BD2):
 |
(20) |
the explicit fourth-order backward-difference time integration scheme (BD4):
 |
(21) |
or the implicit backward-difference time integration scheme (ID
):
 |
(22) |
European Options Numerical Specification
Analytical Results
In order to evaluate European options, we need to specify their boundary conditions. Remember that all options will satisfy the Black-Scholes PDE of equation (2) with the initial condition of equation (3), but that for further differentiation between options we also need to specify the boundary conditions as 
and 
:
 |
(23) |
By using these equations we can now implement the numerical solution of European options. Even though we have given equations for both call- and put-style options, for the sake of brevity and to avoid what would essentially be very similar programs, we restrict ourselves to put options only. The time domain is subdivided into 
subintervals, 
, so that for any prior time 
, 
lies in the range 
, with 
corresponding to the maturity date 
. The spatial discretization will be for 
collocation points 
, 
. The Mathematica program for the Black-Scholes European options and their deltas according to equations (4) and (5) can be readily implemented.
For the sake of reproducing actual prices, we specify the following market parameters. The risk-free interest rate is 1%, the dividend yield is 0, the volatility is 30%, the strike is $100, 
of 1 gives a 
of 0, and 
. Let there be 
collocation points (Npts) and subdivide the time to maturity 
year into 
subintervals. We chose rather large values for 
and 
that can be reused later for the American options, but are more than sufficient for evaluating European options. We later consider the implicit time recursion model with 
, which corresponds to the Crank-Nicolson FD scheme. An upper limit on recursion is set at 1000. Further, observe that we calculate not just one value for the analytical and approximation results but all of the possible values for the option over its spatial and time domains.
Our purpose is to compare the exact results for European options from the given Black-Scholes formulas with the analytical results of equation (17), explicit formulas of equations (19) to (21), and the implicit time integration scheme of equation (22). Note that 
is the RBF as described earlier in equations (1), (10), and (11). With 
and 
the smallest and largest values of the stock price, then the spatial increment in the 
domain with 
(or Npts) collocation points is 
and choosing 
in 
, we can now define 
, 
, and 
according to (11) as:
The definition of the matrix 
is then determined from (14) and used to calculate the eigensystem of 
as 
. The initial condition of (7) is used to determine 
(or UT), which is substituted into (18) to calculate 
(or kvec), and this is further substituted into (15) to yield 
.
Initial and boundary conditions of (7) and (23) are used to obtain 
and 
. The general formulas for 
, 
, 
can now be specified by (17). Finally, recalling that 
, 
is the elapsed time, and 
is the time remaining, we can get the analytical formula for 
and its hedging delta in terms of 
by simply taking the derivatives of the RBFs 
.
Explicit and Implicit Results
The explicit and implicit backward-difference time integration schemes given by equations (19) to (22) can now be implemented using straightforward Do loops. We start at time 
with 
set to UT as determined by (7) and step backward through time, considering earlier times 
, while defining 
, 
. There is one further modification that needs to be undertaken in each step of the Do loop; namely, that the boundary condition (23) must be satisfied. Note that 
in the exponent. This modification occurs at 
and hence when 
. This necessitates the rewriting of 
(or U[[1]]).
Now the initial and boundary conditions according to equations (7) and (23) are specified.
Again the general formulas for 
, 
; 
can now be specified according to the left-hand side of equation (17). We can also get the explicit and implicit approximation formulas for 
by using the updated values for 
. The hedging delta formulas are given in terms of 
by simply taking the derivatives of the RBFs 
.
American Options Numerical Specification
Analytical Results
American options will satisfy the same boundary conditions of equation (23) as European options when 
and 
, but must also satisfy the early exercise condition. This means that the current American value 
is the maximum of the conditional expected future discounted value 
, which is the same as the European option price and the exercise price 
. That is,
 |
(24) |
While we can implement this formula directly into the RBF formalism simply by taking the European option code in the previous section and comparing it with the exercise prices, this will not yield satisfactory results. The American options are path-dependent and can be exercised at any time, so they do not have one final boundary at maturity, but an infinite number of boundaries that describe the moving free boundary. There is no explicit formula for the American put, and it is likely that there will only ever be numerical approximations.
Explicit and Implicit Results
The explicit and implicit backward-difference time integration schemes given by equations (19) to (22) can again be implemented using Do loops. As before, start at time 
with 
set to UT as determined by (7) and step backward through time considering earlier times 
while defining 
, 
. Again observe that 
in the exponent. As with the European options, there are further modifications that need to be undertaken in each step of the Do loop; namely, that the exercise condition (23) must be satisfied for all points 
. At each time step 
and each collocation point 
we determine the 
European option value of 
and then compare it with the exercise value 
.
Here the initial and boundary conditions according to equations (7) and (23) are specified.
Again, the general formulas for 
, 
, 
can now be specified according to the left-hand side of equation (17). We can also get the explicit and implicit approximation formulas for 
by using the updated values for 
. The hedging delta formulas are given in terms of 
by simply taking the derivatives of the RBFs 
.
For the sake of calculating actual prices, we use the same market parameters as before. As with the European options, observe that we calculate not just one value for the analytical and approximation results but all of the possible values for the American put option over its spatial and time domains. The time taken for these calculations is only incrementally longer than for the European case. In order to compare the RBF values for the American put with other schemes, we use the equal jumps additive binomial model described in Clewlow and Strickland [10], Chapter 2 and implemented here in Mathematica.
We use this function with a large number 
steps to calculate highly accurate American put prices at the same stock values as for the European options. By varying the stock values by the incremental amount of 
, we can additionally get values for the American put deltas.
In Tables 4 and 5 we can see that while the RBF time integration schemes work very well to approximate the American put, the analytical result is not a sufficiently sensitive measure of the influence of the moving boundary.
A final check on the accuracy of the results for American options can be obtained by loading the powerful financial software UnRisk (www.unriskderivatives.com). It uses compiled backward-recursive path integration over a large grid to obtain highly accurate estimates of path-dependent and American option values. The values returned by the calculation are {Put value, Delta, Gamma, Theta, Vega, VolConvexity, DeltaVega}. In Tables 6, 7, and 8 we can see that while the RBF backward-integration methodology works very well to approximate the American put, the analytical result is still not sufficiently sensitive to the behavior of the moving boundary.
Optimal Exercise Boundary for the American Put
Furthermore, these programs can now be used to efficiently determine the difficult numerical problem of the moving free boundary for the American option. The concept of the free boundary is that it is the function 
that describes the point at time 
where for values of 
the stock price is not so small as to warrant early exercise, and it remains better to hold onto the option, so that it behaves like the European option and satisfies the Black-Scholes PDE of equation (2). However, if 
, then the stock has become sufficiently small as to warrant immediate exercise. On the boundary the value of the option is the exercise price resulting in the equation
 |
(25) |
Consider the time 
; then at that time 
and 
. We now define these 
functions.
To check that the zero function 
is working properly, we use FindRoot, starting at 
.
At time 
when 
and hence 
. For successive values of 
, we move further back in time to the present and 
determines those values 
such that 
. This procedure is implemented using NestList, in which successive solutions are used as the starting point for the next evaluation of the FindRoot function.
Figure 1. Moving free boundary determined by RBF.
3D Graphs of American Options
The given code evaluates all American option values for every point in time 
and at every collocation point 
. This yields a complete picture of the American put as it evolves in time and space, as shown by the following.
Figure 2. 3D plot of the American put determined by the implicit RBF.
Conclusion
In this article we investigated the use of radial basis functions (RBFs) as a means of solving parabolic partial differential equations (PDEs) associated with the solution of European- and American-style financial options. Unlike the usual approach, which typically involves finite difference (FD) schemes, we have used a spatial collocation approach that has resulted in analytical, explicit, and implicit evaluation schemes. We have shown that the RBF methodology can be easily represented using Mathematica‘s powerful programming idiom. Furthermore, there are a number of advantages of this approach over that of the FD method. It is quicker and yields more accurate results in the same time. It is comprehensive, resulting in a complete description of the path of the option for a complete range of values in space and time, as was displayed in Figures 1 and 2. And because the RBFs are themselves readily differentiated, the hedging parameters are also immediately calculable. Lastly, there is a range of parameters such as the shape parameter 
, the number of time steps 
, and the spatial dimension 
that can be adjusted so as to further improve results.
References
| [1] | R. L. Hardy, “Multiquadratic Equations of Topography and Other Irregular Surfaces,” Journal of Geophysical Research, 76(8), 1971 pp. 1905-1915. |
| [2] | R. L. Hardy, “Theory and Applications of the Multiquadric-Biharmonic Method,” Computers and Mathematics with Applications, 19(8-9), 1990 pp. 163-208. |
| [3] | E. J. Kansa, “Multiquadrics—A Scattered Data Approximation Scheme with Applications to Computational Fluid Dynamics I: Surface Approximations and Partial Derivative Estimates,” Computers and Mathematics with Applications, 19(8-9), 1990 pp. 127-145. |
| [4] | E. J. Kansa, “Multiquadrics—A Scattered Data Approximation Scheme with Applications to Computational Fluid Dynamics II: Solutions to Parabolic, Hyperbolic and Elliptic Partial Differential Equations,” Computers and Mathematics with Applications, 19(8-9), 1990 pp. 147-161. |
| [5] | Y-C. Hon and X-Z. Mao, “A Radial Basis Function Method for Solving Options Pricing Models,” Financial Engineering, 8, 1999 pp. 31-49. |
| [6] | G. Fasshauer, A. Q. M. Khaliq, and D. A. Voss, “Using Meshfree Approximations for Multi-Asset American Option Problems,” Journal of Chinese Institute of Engineers (JCIE), 27(4), 2004 pp. 563-571. |
| [7] | F. Black and M. Scholes, “The Pricing of Options and Corporate Liabilities,” Journal of Political Economy, 81(3), 1973 pp. 637-654. |
| [8] | M. A. Goldberg and C. S. Chen, “On a Method of Atkinson for Evaluating Domain Integrals in the Boundary Element Method,” Applied Mathematics and Computation, 80(2-3), 1994 pp. 125-138. |
| [9] | M. J. D. Powell, “The Theory of Radial Basis Function Approximations in 1990,” Advances in Numerical Analysis: Wavelets, Subdivision Algorithms, and Radial Basis Functions, Vol. II (W. A. Light, ed.), Oxford: Oxford University Press, 1992 pp. 105-210. |
| [10] | L. Clewlow and C. Strickland, Implementing Derivatives Models (Wiley Series in Financial Engineering), New York: John Wiley & Sons, 1998. |
| Michael Kelly, “Evaluation of Financial Options Using Radial Basis Functions in Mathematica,” The Mathematica Journal, 2010. dx.doi.org/doi:10.3888/tmj.11.3-3. | |
About the Author
Michael Kelly is currently the Research Director for Unetich Trading LLC, which operates at the Chicago Board of Trade. He designs automated trading programs using Mathematica that implement trading strategies based upon statistical rules. Previously he was a professor of mathematical and computational finance in the Stuart School of Business, Illinois Institute of Technology and a derivatives consultant. Prior to that he was senior lecturer in both mathematics and finance in what is now the School of Computing and Mathematics and the School of Economics and Finance, respectively, at the University of Western Sydney, Australia. He has also been a visiting research fellow at the University of Sydney, the Australian National University and the University of Warwick (U.K.). For at least 15 years his lectures and research have been based upon computational mathematics and maths programs such as Mathematica and Matlab.
Michael Kelly
Research Director
Unetich Trading LLC, Suite 2020, Chicago Board of Trade
141 W. Jackson Boulevard, Chicago, Illinois 60604-2994
mk1444@rcn.com