ABSTRACT

The N-body problem is an active research topic in physics for which there are two major algorithms for efficient computation, the fast multipole method and treecode, but these algorithms are not popular in financial engineering. In this article, we apply a fast N-body algorithm called the Cartesian treecode to the computation of the integral operator of integro-partial differential equations to compute option prices under the CGMY model, a generalization of a jump-diffusion model. We present numerical examples to illustrate the accuracy and effectiveness of the method and thereby demonstrate its suitability for application in financial engineering.

Keywords:

Jump-Diffusion Model, CGMY Model, Option Pricing, Cartesian Treecode

1. Introduction

The standard and celebrated model for option pricing in the financial industry has been the Black-Scholes model [1] because of its simplicity. However, it is difficult to capture the realistic behaviors observed in option markets with this model. Recently, jump-diffusion models, which take into account jump behavior of underlying assets, have gained attention in financial engineering (the details are given in [2] , for example). For a European option, researchers have suggested efficient schemes which apply numerical methods such as the fast Fourier transform (FFT) [3] , the Hilbert transform [4] , and the Fourier-cosine (COS) [5] methods. The application of these methods to path-dependent options has also been investigated (for example, [6] ).

One approach to price options under jump-diffusion models is to solve the corresponding partial integro-differential equations (PIDEs) numerically, but the PIDEs include non-local integral operators which make the direct numerical evaluation of PIDEs computationally expensive. An efficient approach for computing the integral term is to use the FFT [7] [8] [9] [10] , but this requires two FFT operations per time-step and direct use of FFT requires either a uniform grid or an additional interpolation scheme [11] . Also, an efficient scheme which transforms the PIDEs to pseudo-differential equations and applies FFT to the equations has been suggested [12] , but this also requires two FFT operations per time-step for path-dependent option pricing.

The present paper introduces an application of the Cartesian treecode [13] which has been used as an efficient algorithm for the N-body problem. The N-body problem, which simulates the interactions between individual particles in a multi-particle system, is an active research topic in physics. However, direct computation for this problem is computationally expensive and there are two major efficient algorithms: the fast multipole method (FMM, [14] [15] [16] [17] ) and treecode [13] [18] . FMM is relatively complicated and the performance is not sufficient for the actual computation. On the other hand, although the com- putational burden of treecode is greater, its implementation is straightforward and easier to use in practice.

Little study has been done on the application of these fast N-body algorithms to financial engineering. One application of FMM called the fast Gauss transform (FGT, [16] ) has been implemented under Merton’s jump-diffusion model [19] to solve the corresponding PIDE [9] . However, its computational efficiency was observed to be poor because the method requires more grid points than the FFT method to achieve similar accuracy. Recently, the improved fast Gauss transform (IFGT, [20] ) has been applied and numerical experiments have shown that this method is more efficient than the FFT method [21] ; however, IFGT can be applied only to Merton’s jump-diffusion models.

In the present paper, we apply another efficient approach, called the Cartesian treecode [13] , to the CGMY model ( [11] [22] ). To our best knowledge, this is the first paper applying treecode to option pricing modeling. Using numerical examples, we examine its computational accuracy and ease of implementation for the computation of the integral operators of the PIDE. The use of the treecode makes the computation of option pricing faster than the original finite difference method without losing accuracy.

Furthermore, we examine whether treecode is also applicable in the field of financial engineering.

2. Cartesian Treecode

In the N-body problem studied in physic, the main goal is to compute the summation of pairwise interactions over the particles in a given system (e.g., gravitational forces of stars in a galaxy or the Coulomb forces of atoms in a molecule). Herein, we consider N particles located at $x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_N$ . Then given coefficients $C_1\mathrm{,}{C}_2\mathrm{,}\cdots \mathrm{,}{C}_N$ and potential $K\left(x\mathrm{,}y\right)$ which describes some relation between particles at positions x and y, we want to compute

$F\left(x_i\right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}}{C}_j\cdot K\left(x_i,z_j\right).$ (1)

However, naive computation of the direct summation results in an $O\left(N^2\right)$ computational cost. The Cartesian treecode [13] is a type of treecode used to evaluate the screened Coulomb (Yukawa) potential efficiently, which involves potential $K\left(x_i\mathrm{,}{z}_j\right)$ given by

$K\left(x_i,z_j\right):=\frac{{\text{e}}^{-\kappa \left|z_j-x_i\right|}}{\left|z_j-x_i\right|}.$ (2)

The intuitive idea behind treecode is that one imagines wanting to compute gravitational force $F\left(x_i\right)$ of a star located at some point $x_i$ . Given clusters of stars which are very far from the target star, instead of summing the gravita- tional force of all the stars in the clusters, if D/r is very small, where D is the size of the cluster and r is the distance of the cluster point from the target point, then we can approximate the force for each cluster as that of a star whose mass is equal to the total sum of the stars in the cluster and is located at the center of the cluster (called the cluster point). Treecode applies the Barnes-Hut algorithm [18] to build a hierarchy of such clusters and employs an efficient recursive compu- tation of (far-field) Taylor coefficients to compute particle-cluster approxima- tions of the screened Coulomb potential, reducing the computational cost from $O\left(N^2\right)$ to $O\left(N\log N\right)$ .

A more concrete explanation of the Cartesian treecode is as follows (see [13] for details). First, particles $\left(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_n\right)$ are divided into a hierarchy of clusters c’s by the Barnes-Hut algorithm,

$F\left(x_i\right)={\displaystyle \underset{c}{\sum }}{F}_c\left(x_i\right),F_c\left(x_i\right)={\displaystyle \underset{z_j\in c}{\sum }}{C}_i\cdot K\left(x_i,z_j\right).$ (3)

Second, truncated expansion of the Taylor series of $K\left(x_i\mathrm{,}{z}_j\right)$ up to order p around cluster point $z_c$ gives

$K\left(x_i,z_j\right)={\displaystyle \underset{k=0}{\overset{p}{\sum }}}\frac{1}{k!}{\frac{{\partial }^{k}K\left(x_i,z\right)}{\partial z^k}|}_{z=z_c}{\left(z_j-z_c\right)}^k.$ (4)

For each cluster c, the Taylor expansion is applied if $\frac{r}{R}\le \theta$ , where r is the

radius of the cluster, R is the distance from $x_i$ to the cluster c, and $\theta$ is the multipole acceptance criterion (MAC) parameter. If the condition is not satisfied, then the children of c are examined, or simple direct summation is applied if the cluster is a leaf of the tree. Finally, substituting (4) into (3) gives

$\begin{array}{c}{F}_c\left(x_i\right)={\displaystyle \underset{z_j\in c}{\sum }}{C}_j\cdot {\displaystyle \underset{k=0}{\overset{p}{\sum }}}\frac{1}{k!}{\frac{{\partial }^{k}K\left(x_i,z\right)}{\partial z^k}|}_{z=z_c}{\left(z_j-z_c\right)}^k\\ ={\displaystyle \underset{k=0}{\overset{p}{\sum }}}\frac{1}{k!}{\frac{{\partial }^{k}K\left(x_i,z\right)}{\partial z^k}|}_{z=z_c}{\displaystyle \underset{z_j\in c}{\sum }}{C}_j{\left(z_j-z_c\right)}^k:={\displaystyle \underset{k=0}{\overset{p}{\sum }}}{L}_k\left(x_i,z_c\right)M_{c,k}.\end{array}$ (5)

Because $L_k\left(x_i\mathrm{,}{z}_c\right)$ is independent of $z_j$ and $M_{c\mathrm{,}k}$ is independent of $x_i$ , $L_k\left(x_i\mathrm{,}{z}_c\right)$ and $M_{c\mathrm{,}k}$ can be calculated independently and computation of $F_c\left(x_i\right)$ is more efficient (notice that in order to compute (1), the summation is required at each $x_i$ ). Furthermore, the upshot of this Cartesian treecode is that the Taylor coefficients $L_k\left(x_i\mathrm{,}{z}_c\right)$ are evaluated recursively.

Generally, FMM is a more sophisticated model for N-body problems in the sense that it also uses near-field expansions, which reduces the computational cost from $O\left(N\log N\right)$ to $O\left(N\right)$ , but the additional use of near-field expansions complicates the algorithm and affects the actual performance, whereas the treecode method is simple and straightforward to implement.

3. Application to Solving a PIDE

We apply the Cartesian treecode to option pricing under the CGMY model with diffusion coefficients ( [11] [22] ). We assume the underlying stock price S to follow a jump process so that the European option price $V\left(S\mathrm{,}t\right)$ satisfies the following PIDE (see [2] for details of the model and the derivation of the PIDE):

$\begin{array}{c}{V}_t=\frac{1}{2}{\sigma }^2{S}^2{V}_{SS}+\left(r-q\right)S{V}_S-rV\\ +{\displaystyle {\int }_R}\left[V\left(S{\text{e}}^z,t\right)-V\left(S,t\right)-S\left({\text{e}}^z-1\right)V_S\right]f\left(z\right)\text{d}z,\end{array}$ (6)

where $\omega ={\displaystyle {\int }_R}\left(1-{\text{e}}^z\right)f\left(z\right)\text{d}z,$ r is the risk-free rate, q is the dividend, t is the time to maturity, and $f\left(z\right)$ is the jump measure. The initial condition is given by

$V\left(S,0\right)=\max \left(S-K,0\right),$ (7)

where K is the strike price. Substituting $x=\ln \left(S/K\right)$ gives

$\begin{array}{c}{V}_t=\frac{1}{2}{\sigma }^2{V}_{xx}+\left(r-q\right)V_x-rV\\ +{\displaystyle {\int }_R}\left[V\left(x+z,t\right)-V\left(x,t\right)-\left({\text{e}}^z-1\right)V_x\left(x,t\right)\right]f\left(z\right)\text{d}z\end{array}$ (8)

with

$V\left(x,0\right)=\max \left(K{\text{e}}^x-K,0\right).$ (9)

The major drawback of this approach is that we need to compute the integral operator given in this PIDE at each time-step; however, use of the Cartesian treecode makes this efficient. Under the CGMY model, $f\left(z\right)$ is expressed as

$f_{CGMY}\left(y\right)=C\frac{{\text{e}}^{-G\left|y\right|}}{{\left|y\right|}^{1+Y}}\cdot 1_{y<0}+C\frac{{\text{e}}^{-My}}{y^{1+Y}}\cdot 1_{y>0}$ (10)

for $C>0$ , $G\ge 0$ , $M\ge 0$ , and $Y<2$ . The CGMY model has both finite and infinite activity: in the case of $Y<0$ , it has finite activity. Also, the CGMY mo- del has both finite and infinite variance: if $0<Y<1$ , the process has finite variance, whereas it has infinite variance in the case of $1<Y<2$ . As mentioned previously, the CGMY model is a generalization of a jump-diffusion model in the sense that $Y=0$ corresponds to what is also referred to as the variance gamma (VG) model [23] and $Y=-1$ corresponds to Kou’s jump-diffusion model [24] .

Furthermore, a more general jump-diffusion model exists known as the generalized tempered stable model [2] ,

$f_{alpha}\left(y\right)=C_1\frac{{\text{e}}^{-G\left|y\right|}}{|y{|}^{1+Y_1}}\cdot 1_{y<0}+C_2\frac{{\text{e}}^{-My}}{y^{1+Y_2}}\cdot 1_{y>0},$ (11)

for $C_1,C_2>0$ and $Y_1,Y_2<2$ , and application of the Cartesian treecode is straightforward.

For numerical evaluation of the PIDE, we employed a simple scheme of the Crank-Nicholson method for the case of a uniform discretization of points, $0\le i\le N_S$ and $0\le j\le N_T$ , where $N_S$ and $N_T$ represent the numbers of dis- cretization points and time steps, respectively. Under this scheme, we decom- pose the integral part as

$\begin{array}{l}{\displaystyle {\int }_R}\left[V\left(x+z,t\right)-V\left(x,t\right)-\left({\text{e}}^z-1\right)V_x\left(x,t\right)\right]f_{CGMY}\left(z\right)\text{d}z\\ \approx {\displaystyle {\int }_{x_{\min }}^{x_i-\Delta x}}\left[V\left(z\right)-V\left(x_i\right)\right]f_{CGMY}\left(z-x_i\right)\text{d}z\\ +{\displaystyle {\int }_{x_i+\Delta x}^{x_{\max }}}\left[V\left(z\right)-V\left(x_i\right)\right]f_{CGMY}\left(z-x_i\right)\text{d}z\\ -V_x\left(x_i\right){\displaystyle {\int }_{x_{\min }}^{-\Delta x}}\left({\text{e}}^z-1\right)f_{CGMY}\left(z\right)\text{d}z-V_x\left(x_i\right){\displaystyle {\int }_{\Delta x}^{x_{\max }}}\left({\text{e}}^z-1\right)f_{CGMY}\left(z\right)\text{d}z\\ +{\displaystyle {\int }_{-\Delta x}^{\Delta x}}\left[V\left(x_i+z\right)-V\left(x_i\right)-\left({\text{e}}^z-1\right)V_x\left(x_i\right)\right]f_{CGMY}\left(z\right)\text{d}z,\end{array}$ (12)

where we employ the following approximation [25] for the last term:

$\begin{array}{l}{\displaystyle {\int }_{-\Delta x}^{\Delta x}}\left[V\left(x_i+z\right)-V\left(x_i\right)-\left({\text{e}}^z-1\right)V_x\left(x_i\right)\right]f_{CGMY}\left(z\right)\text{d}z\\ \cong \frac{{\left(\Delta x\right)}^3}{2}\left[f_{CGMY}\left(-\Delta x\right)+f_{CGMY}\left(\Delta x\right)\right]\cdot V_{xx}\left(x_i\right).\end{array}$ (13)

Then the PIDE is solved numerically via a simple Crank-Nicholson scheme:

$\begin{array}{l}\frac{V_i^{j+1}-V_i^j}{\Delta t}-\frac{r-q-L_2}{2}\left[\frac{V_{i+1}^j-V_{i-1}^j}{2\Delta x}+\frac{V_{i+1}^{j+1}-V_{i-1}^{j+1}}{2\Delta x}\right]+\frac{r}{2}\left[V_i^j+V_i^{j+1}\right]\\ -\frac{L_1}{2}\left[\frac{V_{i-1}^j-2V_i^j+V_{i+1}^j}{\Delta x^2}+\frac{V_{i-1}^{j+1}-2V_i^{j+1}+V_{i+1}^{j+1}}{\Delta x^2}\right]\\ =L_3\left(x_i\right)-L_4\left(x_i\right),\end{array}$ (14)

where

$L_1:=\frac{{\left(\Delta x\right)}^3}{2}\left[f_{CGMY}\left(-\Delta x\right)+f_{CGMY}\left(\Delta x\right)\right],$

$L_2:={\displaystyle {\int }_{x_{\min }}^{-\Delta x}}\left({\text{e}}^z-1\right)f_{CGMY}\left(z\right)\text{d}z+{\displaystyle {\int }_{\Delta x}^{x_{\max }}}\left({\text{e}}^z-1\right)f_{CGMY}\left(z\right)\text{d}z,$

$L_3\left(x_i\right):={\displaystyle {\int }_{x_{\min }}^{x_i-\Delta x}}V\left(z\right)f_{CGMY}\left(z-x_i\right)\text{d}z+{\displaystyle {\int }_{x_i+\Delta x}^{x_{\max }}}V\left(z\right)f_{CGMY}\left(z-x_i\right)\text{d}z,$

$L_4\left(x_i\right):=V\left(x_i\right)\cdot \left[{\displaystyle {\int }_{x_{\min }}^{x_i-\Delta x}}{f}_{CGMY}\left(z-x_i\right)\text{d}z+{\displaystyle {\int }_{x_i+\Delta x}^{x_{\max }}}{f}_{CGMY}\left(z-x_i\right)\text{d}z\right].$

The boundary conditions are given by

$V\left(x\mathrm{,}t\right)(\begin{array}{ll}\to \mathrm{0,}\hfill & x\to -\infty \mathrm{,}\hfill \\ \simeq K{\text{e}}^x-K{\text{e}}^{-r\tau }\mathrm{,}\hfill & x\to +\infty \mathrm{.}\hfill \end{array}$

Use of an iterative scheme [11] is also possible as a means of maintaining second-order accuracy with respect to time, but for simplicity $L_3\left(x_i\right)$ and $L_4\left(x_i\right)$ are computed explicitly. $L_3\left(x_i\right)$ is computed via the Cartesian tree- code and in this case the Taylor coefficients of $f_{CGMY}$ are computed recursively as described in the next section. In the numerical evaluation of $L_2$ , $L_3\left(x_i\right)$ , and $L_4\left(x_i\right)$ , approximation using the trapezoid rule is employed. In addition, the construction of the hierarchy of clusters is executed only at the initial time-step because the grid points are fixed under the Crank-Nicolson scheme.

4. Recursive Computation of $L_3\left(x_i\right)$

The Taylor coefficients can be evaluated as follows. Let $\Phi =\frac{{\text{e}}^{-\kappa \left|z-x\right|}}{{\left|z-x\right|}^{1+Y}}$ and $\Psi =\frac{{\text{e}}^{-\kappa \left|z-x\right|}}{\left|z-x\right|},$ where $\kappa :=G\cdot 1_{z-x<0}+M\cdot 1_{z-x>0}$ . Then ${\left|z-x\right|}^Y\Phi =\Psi$ and differ-

entiating with respect to z k times ( $k\ge 2$ ) gives

$\begin{array}{l}{\left|z-x\right|}^2\frac{{\partial }^k\Phi }{\partial z^k}+\left[2\left(k-1\right)+Y\right]\left(z-x\right)\frac{{\partial }^{k-1}\Phi }{\partial z^{k-1}}\\ +\left(k-1\right)\left[\left(k-2\right)+Y\right]\frac{{\partial }^{k-2}\Phi }{\partial z^{k-2}}=\frac{{\partial }^{k-1}}{\partial z^{k-1}}\left({\left|z-x\right|}^{2-Y}\frac{\partial \Psi }{\partial z}\right)\mathrm{.}\end{array}$ (15)

On the other hand,

${\left|z-x\right|}^{2-Y}\frac{\partial \Psi }{\partial z}=\left[1+\kappa \left|z-x\right|\right]\Phi \mathrm{,}$ (16)

and differentiation with respect to z k times ( $k\ge 2$ ) gives

$\begin{array}{c}\frac{{\partial }^{k-1}}{\partial z^{k-1}}\left({\left|z-x\right|}^{2-Y}\frac{\partial \Psi }{\partial z}\right)=-\left(k-1\right)\left(k-2\right)\kappa {\left|z-x\right|}^{-1}\left(z-x\right)\frac{{\partial }^{k-3}\Phi }{\partial z^{k-3}}\cdot 1_{k\ge 3}\\ -\left(k-1\right)\left(1+2\kappa \left|z-x\right|\right)\frac{{\partial }^{k-2}\Phi }{\partial z^{k-2}}\\ -\left(z-x\right)\left[1+\kappa \left|z-x\right|\right]\frac{{\partial }^{k-1}\Phi }{\partial z^{k-1}}.\end{array}$ (17)

From (15) and (17), we get

$\begin{array}{l}{\left|z-x\right|}^2\frac{{\partial }^k\Phi }{\partial z^k}+\left[2\left(k-1\right)+Y\right]\left(z-x\right)\frac{{\partial }^{k-1}\Phi }{\partial z^{k-1}}+\left(k-1\right)\left[\left(k-2\right)+Y\right]\frac{{\partial }^{k-2}\Phi }{\partial z^{k-2}}\\ =-\left(k-1\right)\left(k-2\right)\kappa {\left|z-x\right|}^{-1}\left(z-x\right)\frac{{\partial }^{k-3}\Phi }{\partial z^{k-3}}\cdot 1_{k\ge 3}\\ -\left(k-1\right)\left(1+2\kappa \left|z-x\right|\right)\frac{{\partial }^{k-2}\Phi }{\partial z^{k-2}}-\left(z-x\right)\left[1+\kappa \left|z-x\right|\right]\frac{{\partial }^{k-1}\Phi }{\partial z^{k-1},}\end{array}$

and

$\begin{array}{c}{\left|z-x\right|}^2\frac{{\partial }^k\Phi }{\partial z^k}=-\left(z-x\right)\left[2k-1+Y+\kappa \left|z-x\right|\right]\frac{{\partial }^{k-1}\Phi }{\partial z^{k-1}}\\ -\left(k-1\right)\left(k-1+Y+2\kappa \left|z-x\right|\right)\frac{{\partial }^{k-2}\Phi }{\partial z^{k-2}}\\ -\left(k-1\right)\left(k-2\right)\kappa {\left|x-z\right|}^{-1}\left(z-x\right)\frac{{\partial }^{k-3}\Phi }{\partial z^{k-3}}\cdot 1_{k\ge 3}.\end{array}$ (18)

Finally, dividing by $k\mathrm{<mo>\; !\; </mo>}$ and substituting $x=x_i$ and $z=z_c$ gives the following recursive equation for $L^k\left(x_i\mathrm{,}{z}_c\right)$ :

$\begin{array}{l}k{\left|z_c-x_i\right|}^2{L}^k\left(x_i,z_c\right)\\ =-\left(z_c-x_i\right)\left[2k-1+Y+\kappa \left|z_c-x_i\right|\right]L^{k-1}\left(x_i,z_c\right)\\ -\left(k-1+Y+2\kappa \left|z_c-x_i\right|\right)L^{k-2}\left(x_i,z_c\right)\\ -\kappa {\left|z_c-x_i\right|}^{-1}\left(z_c-x_i\right)L^{k-3}\left(x_i,z_c\right)\cdot 1_{k\ge 3}.\end{array}$ (19)

5. Numerical Examples

¹The code is available at https://bitbucket.org/cdcooper/yukawatree

We conducted the following numerical experiments, which were executed in Python. The following parameter settings were used [22] : r = 0.06, q = 0.00, C = 0.42, G = 4.37, M = 191.2, Y = 1.0102, and T = 0.25. In the case of the Cartesian treecode, we used publicly available Python code1 using the settings $\theta =0.5$ , P = 15, and MP = 30, where $MP$ is the maximum number of particles in each node of the hierarchy of clusters, and so using a smaller value for MP results in a deeper hierarchy of clusters. The vertical axis in Figure 1 above indicates the price error computed as the difference between the reference price of the Fourier transform method [3] and the price computed using (14), where $L_3\left(x_i\right)$ is computed using the Cartesian treecode (Treecode 1). Treecode 2 indicates the difference between the reference price of the Fourier transform method [3] and the price whose $L_3\left(x_i\right)$ is computed at only one-eighth the number of grid points $x_i$ equally distributed and spline interpolation is employed for the other grid points. The vertical axis in Figure 1 bellow indicates the computation time in case of Treecode 1 and Treecode 2. As shown, there is little difference in Price Error between the cases of FFT and Treecode 1.

The flexibility of treecode allows it to achieve higher efficiency of decreasing computational time by choosing target grid points to compute $L_3\left(x_i\right)$ , as is done in Treecode 2.

For Treecode 2, because $L_3\left(x_i\right)$ is a smooth curve with respect to $x_i$ , it suffices to only compute one-eighth the number of grid points $x_i$ equally distributed and interpolating the other grid points via a spline scheme for evaluation of $L_3\left(x_i\right)$ . As shown in Figure 1, we can achieve similar price errors in the case of Treecode 2, while reducing computational time. This method would be especially useful for path-dependent option pricing and a detailed investigation is left for future research.

Figure 2 shows the price and delta values of a European call under the CGMY model for different values of Y using Treecode 2. Here, Y = 0.75 indicates that the process has finite variance and therefore the prices are lower than the case of Y = 1.0102, whose process has infinite variance. Figure 3 shows the price and delta values of a European call under the CGMY model with Y = 0.0 (therefore this model corresponds to a VG model) using Treecode 2, where the following parameter settings [26] were used: r = 0.00, q = 0.0, C = 5.9311, G = 20.2648, M = 39.784, Y = 0.00, and T = 0.5. Both figures demonstrate the smoothness of price and delta with respect to stock price.

6. Conclusion

In this article, we describe applying one of the fast N-body algorithms used in physics, the Cartesian treecode [13] , to the evaluation of the integral operator in a PIDE under the CGMY model. Numerical examples are presented and the results verify the accuracy and efficiency of the treecode method. The recursive computation of coefficients is efficient and by choosing the grid points to be computed and applying interpolation for the other grid points, we can achieve

high efficiency compared to the original finite difference method without losing accuracy. In addition, the Cartesian treecode can be directly applied to a non- uniform grid, and the use of transformation maps [10] to increase accuracy is also possible. In this paper, we consider a one-dimensional case, but the treecode method is applicable to up to three-dimensional cases, and therefore we can apply the code to a two-dimensional PIDE for pricing basket and spread options. On the other hand, because the treecode method is only applicable to up to three-dimensional cases, the limitation of this approach is that we cannot compute the price of a multi-asset option which involves more than three assets. Application to path-dependent option pricing is also interesting and left for future research.

Acknowledgements

We appreciate Professor Lorena Barba, Simon Layton, and Christopher Cooper for helpful advice and support. I would especially like to note my appreciation of Christopher Cooper for his support in instructing me on how to use his Python code of the Cartesian treecode.

Cite this paper

Sakuma, T. (2017) Application of Fast N-Body Algorithm to Option Pricing under CGMY Model. Journal of Mathematical Finance, 7, 308-318. https://doi.org/10.4236/jmf.2017.72016

References