Finding an Efficient Computational Solution for the Bates Partial Integro-Differential Equation Utilizing the RBF-FD Scheme

Abstract

This paper proposes a computational solver via the localized radial basis function finite difference (RBF-FD) scheme and the use of graded meshes for solving the time-dependent Bates partial integro-differential equation (PIDE) arising in computational finance. In order to avoid facing a large system of discretization systems, we employ graded meshes along both of the spatial variables, which results in constructing a set of ordinary differential equations (ODEs) of lower sizes. Moreover, an explicit time integrator is used because it can bypass the need to solve the large discretized linear systems in each time level. The stability of the numerical method is discussed in detail based on the eigenvalues of the system matrix. Finally, numerical tests revealed the accuracy and reliability of the presented solver.

1. Introductory Notes

The Bates model for option pricing considers that the underlying asset $S_t$, the volatility $V_t$, the riskless constant r and $N_t$ as the Poisson process satisfy the following system of stochastic differential equations (SDEs) [1]:

$$\begin{array}{c}d{S}_t=\sqrt{V_t}{S}_{t}d{W}_t^1+(-\lambda \xi -q+r)S_{t}dt+(\varrho -1)S_{t}d{N}_t,\hfill \\ \\ d{V}_t=\sigma \sqrt{V_t}d{W}_t^2+\kappa (-V_t+\theta )dt,\hfill \end{array}$$

(1)

wherein $W_t^2$ and $W_t^1$ are standard Brownian motions having $d{W}_t^{1}d{W}_t^2=\rho dt.$ Here $\kappa$ is the reversion’s rate of the variance $V_t$, $\lambda$ is the Poisson process intensity, $\xi$ is the mean jump, q is the dividend, $\varrho$ is the jump size, while $\theta$ is the mean level and $\sigma$ stands for the volatility fixed value.

Financial derivatives such as European call or put options play pioneer roles in the risk management of some portfolios and their pricing as efficiently as possible is of importance. On the other hand for the financial derivative price, since analytical relations are available only in limited settings, one is in need for the construction and the application of fast and stable numerical solvers. More concretely, starting from the initial time zero, we must numerically solve a second-order high-dimensional time-dependent partial integro-differential equation (PIDE) or a partial differential equation (PDE) and then compute the present value of the financial derivative [2,3,4].

The Heston model, which could be considered as a generalization of the Black–Scholes model [5], can be extended further if one follows the consideration of Bates [1,6] by imposing the jump component into the modeling. In fact, in the stochastic volatility jump (SVJ) model, the price of an option is computed by solving a time-dependent 2D PIDE [7,8]. It is requisite to recall some related models [9,10] discussing stochastic volatility for PDEs in control theory and AI.

The Bates PIDE based on the price function $u(x,y,\tau )$ for European options is expressed by the following [11]:

$$\begin{array}{cc}\hfill \frac{\partial u(x,y,\tau )}{\partial \tau }& =\frac{1}{2}y{x}^2\frac{{\partial }^{2}u(x,y,\tau )}{\partial x^2}+\frac{1}{2}{\sigma }^{2}y\frac{{\partial }^{2}u(x,y,\tau )}{\partial y^2}\hfill \\ \hfill & +\rho \sigma yx\frac{{\partial }^{2}u(x,y,\tau )}{\partial x\partial y}\hfill \\ \hfill & +(-\lambda \xi -q+r)x\frac{\partial u(x,y,\tau )}{\partial x}+\kappa (\theta -y)\frac{\partial u(x,y,\tau )}{\partial y}\hfill \\ \hfill & -(\lambda +r)u(x,y,\tau )+\lambda {\int }_0^{\infty }u(x\varrho ,y,\tau )b\left(\varrho \right)d\varrho =\mathcal{A}u(x,y,\tau ),\hfill \end{array}$$

(2)

wherein T is the time to maturity and $\tau =T-t$ is a time transformation to have forward in time PIDE formulation, unlike the original Bates PIDE, which is backward in time. Besides, both the differential and integral parts of the PIDE problem have been encapsulated in the operator $\mathcal{A}$, that is to say, we also can write

$$\mathcal{A}u(x,y,\tau )={\mathcal{A}}_{D}u(x,y,\tau )+\lambda {\mathcal{A}}_{I}u(x,y,\tau ),$$

(3)

in which ${\mathcal{A}}_D$ and ${\mathcal{A}}_I$ stand for the differential and integral portions of the PIDE problem. The probability density function is $b\left(\varrho \right)=\frac{1}{\sqrt{2\pi }\widehat{\sigma }\varrho }exp\left[-\frac{{(ln\left(\varrho \right)-\gamma )}^2}{2{\widehat{\sigma }}^2}\right],$ where it reads ${\int }_0^{\infty }b\left(\varrho \right)d\varrho =1.$ Here, $\widehat{\sigma }$ and $\gamma$ are the standard deviation and the mean, respectively, which are positive constants. Additionally, we have $\xi =exp\left(\gamma +\frac{1}{2}{\widehat{\sigma }}^2\right)-1$.

The so-called payoff which is the the initial condition for the PIDE problem in call-type option pricing can be expressed as [12]:

$$u(x,y,0)={(0,x-K)}^{+},$$

(4)

wherein K is the strike price. The payoff for a put option could be written similarly. The point is that the initial condition is written only on x and does not rely on the second independent variable of the PIDE, i.e., y.

The side conditions for x and y could be given as follows [12]:

$$\begin{array}{cc}\hfill & u(x,y,\tau )\simeq 0,x\to 0,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & u(x,y,\tau )\simeq x_{\max }exp(-q\tau )-Kexp(-r\tau ),x\to x_{max},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & \frac{\partial u(x,y,\tau )}{\partial y}\simeq 0,y\to +y_{max}.\hfill \end{array}$$

(7)

Note that for the case when $y=0$, the PIDE (2) is degenerate and no boundaries should be incorporated while $x_{max}$ and $y_{max}$ are large constants. Similarly for the put option, the boundary conditions are described by the following:

$$\begin{array}{cc}\hfill & u(x,y,\tau )\simeq Kexp(-r\tau )-x_{\max }exp(-q\tau ),x\to 0,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & u(x,y,\tau )\simeq 0,x\to x_{max},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & \frac{\partial u(x,y,\tau )}{\partial y}\simeq 0,y\to y_{max}.\hfill \end{array}$$

(10)

The Bates PIDE (2) is given on $(x,y,\tau )$$\in [0,+\infty )\times [0,+\infty )$$\times (0,T]$. To solve our high-dimensional linear PDE, we must truncate the unbounded domain while quite delicately ignoring the error caused by imposing the boundary conditions. This can be pursued as follows:

$$\mathsf{\Omega }=[0,x_{max}]\times [0,y_{max}],$$

(11)

wherein $x_{max},y_{max}$ are fixed values. The values for $x_{max}$, and $y_{max}$ should be considered large enough to be able to neglect the effect of imposing artificial boundary conditions or imposing the boundary conditions for truncated domains. Some choices are $\mathsf{\Omega }=[0,4K]\times [0,1]$ or $\mathsf{\Omega }=[0,3K]\times [0,1]$.

Assume that ${\left\{x_i\right\}}_{i=1}^m$ is a mesh of nodes for x. The hyperbolic stretching of nodes [13] can be expressed as follows ($1\le i\le m$):

$$x_i=csinh\left({\beta }_i\right)+K,$$

(12)

wherein $c>0$ stands for a fixed value that controls the density around $x=K$ and $m\gg 3$. In implementations, one can employ c as in [14], i.e., $c=K/5$. This puts a focus around the strike price, in which the initial condition of the PIDE has nonsmoothness. Moreover, ${\left\{{\beta }_i\right\}}_{i=0}^m$ stands for the uniform points given by the following:

$${\beta }_i=(i-1)\Delta \beta +{sinh}^{-1}\left(-\frac{K}{c}\right),1\le i\le m,$$

(13)

wherein $\Delta \beta ={(m-1)}^{-1}\left({sinh}^{-1}\left(\frac{S-K}{c}\right)-{sinh}^{-1}\left(\frac{-K}{c}\right)\right)$.

Also, if ${\left\{y_j\right\}}_{j=1}^n$ is a partition for y, then this stretching strategy can be expressed by the following:

$$y_j=sinh\left({\varsigma }_j\right)\nu ,1\le j\le n,$$

(14)

wherein $\nu >0$ is a fixed value that controls the density around $y=0$ and $n\gg 3$. Basically, we use $\nu =K/500$ [14]. Additionally, the ${\varsigma }_j$ are equidistant nodes provided by ${\varsigma }_j=(\Delta \varsigma )(j-1)$, and for any $1\le j\le n$ we have $\Delta \varsigma =\frac{1}{n-1}{sinh}^{-1}\left(\frac{K}{\nu }\right)$.

Numerical solution methods generally utilize the discretization means to realize the approximate calculation. When the computational domain/interval is partitioned more finely, the calculated result is closer to the theoretical solution. Indeed, the time required for the calculation increases. For high-dimensional PIDE problems with kink behavior at the initial conditions, sometimes special solvers such as high-order sparse numerical methods are necessary, see, e.g., [15]. Noting finite difference (FD) methods are discussed in [16,17].

In this paper, the main aim is to propose a novel computational method for resolving (2) via the radial basis function generated finite difference (RBF-FD) methodology [18]. This is mostly because (2) is a (1+2)D problem with variable coefficients, in which there is one cross derivatives. Hence, the computational solvers should be constructed for this aim with much attention. In fact, the motivation of this work lies in the fact that literature lacks the application of efficient RBF-FD methodology which result in fast and sparse procedures for solving the Bates PIDE model. Hence, such an application and investigation on the theoretical stability issues will help price option under stochastic volatility in equity markets.

The RBF-FD formulations in this paper, see, e.g., [19], are written so they can be applied to graded meshes in which there is a clear concentration on the hot zone. The procedure taken here is to employ tensor grids and then time discretize the semi-discretized constructed problem. We note that the present work is related to the pioneering works in [20,21,22]. Meanwhile, these works motivate us to propose a new variant of the RBF-FD scheme for the Bates PIDE problem that competes with these efficient works.

In this paper, after reviewing the well-resulted maps for generating graded meshes along spatial variables with a clear focus around the hot area, the rest of this article is unfolded as follows. The RBF-FD formulas associated with the GMQ RBF are given in Section 2. Then, the semi-discretization of the two-dimensional PDE (2) is described in Section 3. Then in Section 4, an explicitly quadratically convergent method is taken into account. In fact, an explicit time integrator is used because it can avoid to the need to solve the large discretized linear systems in each time level. It is shown that the proposed solver is fast and conditionally stable. The numerical pieces of evidence are given in Section 5, which overwhelmingly uphold the theoretical discussions of the paper. Concluding notes are provided in Section 6.

2. RBF-FD: The Weights

Generally speaking, for computing the weights ${\alpha }_i$ in the methodology of the RBF-FD, one must consider L as a linear operator and at $\underline{x}={\underline{x}}_p$, for the node locations ${\underline{x}}_i$, the following is written down [20]:

$$\left[\begin{array}{c}{\mathsf{\Lambda }}_1\left({\underline{x}}_1\right){\mathsf{\Lambda }}_1\left({\underline{x}}_2\right)\cdots {\mathsf{\Lambda }}_1\left({\underline{x}}_m\right)\\ {\mathsf{\Lambda }}_2\left({\underline{x}}_1\right){\mathsf{\Lambda }}_2\left({\underline{x}}_2\right)\cdots {\mathsf{\Lambda }}_2\left({\underline{x}}_m\right)\\ ⋮⋮⋮\\ {\mathsf{\Lambda }}_m\left({\underline{x}}_1\right){\mathsf{\Lambda }}_m\left({\underline{x}}_2\right)\cdots {\mathsf{\Lambda }}_m\left({\underline{x}}_m\right)\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_m\end{array}\right]=\left[\begin{array}{c}L{\mathsf{\Lambda }}_1\left(\underline{x}\right){|}_{\underline{x}={\underline{x}}_p}\\ L{\mathsf{\Lambda }}_2\left(\underline{x}\right){|}_{\underline{x}={\underline{x}}_p}\\ ⋮\\ L{\mathsf{\Lambda }}_m\left(\underline{x}\right){|}_{\underline{x}={\underline{x}}_p}\end{array}\right],$$

(15)

where the underlined $\underline{x}$ shows a vector quantity in the dimension d and $1\le k\le m$ for some set of test functions $\mathsf{\Lambda }\left(\underline{x}\right)$. It is noted that the extension of RBF-FD methodology for solving computational finance models was revived by the works of Soleymani and co-authors, see for instance [23,24].

Now, we consider the famous generalized multiquadric RBF (GMQ RBF) as follows ([25] Chapter 4):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathsf{\Lambda }\left(r_i\right)={(p^2+r_i^2)}^l,i=1,2,\dots ,\mathfrak{m},\end{array}\end{array}$$

(16)

where l is a suitable parameter, the parameter of shape is p and $r_i=\parallel \mathbf{y}-{\mathbf{y}}_i\parallel$ shows the Euclidean distance.

It is now focused on computing the weights for the GMQ RBF (in the 1D case without loss of generality). So, we consider a graded mesh including three points along the first spatial variable. For finding the weights of the RBF-FD methodology, by taking into account L as an operator, we could write down [26]:

$$L\left[\mathsf{\Lambda }\left(y_j\right)\right]\simeq \sum _{i=1}^{\psi }{\alpha }_i\mathsf{\Lambda }\left(y_i\right),j=1,2,\dots ,\psi .$$

(17)

This gives us $\psi$ unknowns for $\psi$ equations while the solutions will be ${\alpha }_i$. For computing the 1st derivative, three graded nodes are considered ($\psi =3$) as comes next: $\{y_i-h,y_i,y_i+wh\},w>0,h>0,$ and find (17) as follows:

$$g^{\prime }\left(y_i\right)\simeq {\alpha }_{i-1}g\left(y_{i-1}\right)+{\alpha }_{i}g\left(y_i\right)+{\alpha }_{i+1}g\left(y_{i+1}\right).$$

(18)

Noting that we assume that the function g is smooth sufficiently. In estimating the 1st derivative of a function, the analytical weighting coefficients associated to this RBF can be given as follows [22]:

$$\begin{array}{cc}\hfill & {\alpha }_{i-1}={\displaystyle \frac{\omega \left(p^2(9-6l)-h^2(l-1)(4(l-5)\omega -10l+29)\right)}{3p^{2}h(2l-3)(\omega +1)}},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & {\alpha }_i={\displaystyle \frac{(\omega -1)\left(p^2(6l-9)+4h^2(l-5)(l-1)\omega \right)}{3p^{2}h(2l-3)\omega }},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & {\alpha }_{i+1}={\displaystyle \frac{p^2(6l-9)-h^2(l-1)\omega (2l(5\omega -2)-29\omega +20)}{3p^{2}h(2l-3)\omega (\omega +1)}}.\hfill \end{array}$$

(21)

Similarly, in estimating the function’s second derivative, we can obtain

$$g^{″}\left(y_i\right)\simeq \sum _{j=i-1}^{i+1}{\mathsf{\Theta }}_{j}g\left(y_j\right),$$

(22)

along with the following weighting coefficients:

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}_{i-1}=\frac{2\left(p^2(6l-9)-h^2(l-1)\left(4(l-5){\omega }^2+(34-8l)\omega +10l-29\right)\right)}{3p^2{h}^2(2l-3)(\omega +1)},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}_i=\frac{2\left(p^2(9-6l)+h^2(l-1)\left(4(l-5){\omega }^2+(25-2l)\omega +4(l-5)\right)\right)}{3p^2{h}^2(2l-3)\omega },\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}_{i+1}=\frac{2\left(p^2(6l-9)-h^2(l-1)(2l(\omega (5\omega -4)+2)+\omega (34-29\omega )-20)\right)}{3p^2{h}^2(2l-3)\omega (\omega +1)}.\hfill \end{array}$$

(25)

Also noting that the given RBF-FD formulations are valid for the interior nodes and at boundary points, similar formulations must be constructed. We give the derivation for the independent variable y and it would be similar for the other cases. The formulations (19)–(21) and (23)–(24) are useful for the rows two to the row before the last one, while for the 1st and the last rows of the derivative matrices (30) and (31), the weighting coefficients could not be valid on boundaries and sided estimations should be incorporated. Hence, by the work [21] on the stencil $\{y_1,y_2,y_3\}$, we have:

$$g^{\prime }\left(y_1\right)=g[y_2,y_1]-g[y_3,y_2]+g[y_3,y_1]+\mathcal{O}\left({(y_2-y_1)}^2\right),$$

(26)

and

$$g^{\prime }\left(y_m\right)=-g[y_{m-1},y_{m-2}]+g[y_{m-2},y_m]+g[y_m,y_{m-1}]+\mathcal{O}\left({(y_{m-1}-y_m)}^2\right),$$

(27)

wherein $g[l,p]=\left(g\right(l)-g(p\left)\right)/(l-p)$. In a similar manner, for the four nodes $\{\{y_1,g\left(y_1\right)\}$, $\{y_2,g\left(y_2\right)\}$, $\{y_3,g\left(y_3\right)\}$, $\{y_4,g\left(y_4\right)\}\}$, we can obtain

$$\begin{array}{cc}\hfill g^{″}\left(y_1\right)=& \frac{2(\delta y_{1,2}+\delta y_{1,3}+\delta y_{1,4})}{\delta y_{1,2}\delta y_{1,3}\delta y_{1,4}}g\left(y_1\right)+\frac{2(\delta y_{3,1}+\delta y_{4,1})}{\delta y_{1,2}\delta y_{2,3}\delta y_{2,4}}g\left(y_2\right)\hfill \\ \hfill & +\frac{2(\delta y_{2,1}+\delta y_{4,1})}{\delta y_{1,3}\delta y_{3,2}\delta y_{3,4}}g\left(y_3\right)+\frac{2(\delta y_{2,1}+\delta y_{3,1})}{\delta y_{1,4}\delta y_{4,2}\delta y_{4,3}}g\left(y_4\right)+\mathcal{O}\left(h^2\right),\hfill \end{array}$$

(28)

where $\delta y_{l,q}=y_l-y_q$, h is the maximum space width for the considered stencil nodes. Similarly, we have:

$$\begin{array}{cc}\hfill g^{″}\left(y_m\right)=& \frac{2\left(\delta y_{m-3,m}+\delta y_{m-2,m}+\delta y_{m-1,m}\right)}{\delta y_{m-3,m}\delta y_{m,m-2}\delta y_{m,m-1}}g\left(y_m\right)+\frac{2\left(\delta y_{m-3,m}+\delta y_{m-2,m}\right)}{\delta y_{m-3,m-1}\delta y_{m-1,m-2}\delta y_{m-1,m}}g\left(y_{m-1}\right)\hfill \\ \hfill & +\frac{2\left(\delta y_{m-3,m}+\delta y_{m-1,m}\right)}{\delta y_{m-3,m-2}\delta y_{m-2,m-1}\delta y_{m-2,m}}g\left(y_{m-2}\right)+\frac{2\left(\delta y_{m-2,m}+\delta y_{m-1,m}\right)}{\delta y_{m-2,m-3}\delta y_{m-1,m-3}\delta y_{m,m-3}}g\left(y_{m-3}\right)\hfill \\ \hfill & +\mathcal{O}\left(h^2\right).\hfill \end{array}$$

(29)

3. A New Solution Method

Let us use the well-known procedure of method of lines (MOL) for semi discretization of the time-dependent problem [27,28] and convert the PIDE problem into a set of linear ordinary differential equations (ODEs). Hence, the following derivative matrices for the 1st and 2nd derivatives of the function in approximating the PDE problem (2) via semi-discretization are considered on non-uniform stencils given in Section 2 as comes next:

$$M_x=\left\{\begin{array}{ccc}{\alpha }_{i,j}\mathrm{using}\left(19\right)\hfill & \hfill & i-j=1,\hfill \\ {\alpha }_{i,j}\mathrm{using}\left(20\right)\hfill & \hfill & i-j=0,\hfill \\ {\alpha }_{i,j}\mathrm{using}\left(21\right)\hfill & \hfill & j-i=1,\hfill \\ 0\hfill & \hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(30)

and

$$M_{xx}=\left\{\begin{array}{ccc}{\mathsf{\Theta }}_{i,j}\mathrm{using}\left(23\right)\hfill & \hfill & i-j=1,\hfill \\ {\mathsf{\Theta }}_{i,j}\mathrm{using}\left(24\right)\hfill & \hfill & i-j=0,\hfill \\ {\mathsf{\Theta }}_{i,j}\mathrm{using}\left(25\right)\hfill & \hfill & j-i=1,\hfill \\ 0\hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(31)

Consider the $N\times N$ unit matrix $I=I_x\otimes I_y,$ while $N=m\times n$, $I_x$ and $I_y$ are unit matrices of appropriate sizes. The MOL can be resulted in the following coefficient matrix for the 1 + 2 dimensional PIDE:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B=& \frac{1}{2}\mathcal{Y}{\mathcal{X}}^2(M_{xx}\otimes I_n)+\frac{1}{2}{\sigma }^2\mathcal{Y}(I_m\otimes M_{yy})+\rho \sigma \mathcal{Y}\mathcal{X}{M}_{x,y}\hfill \\ \hfill & +(-\lambda \xi -q+r)\mathcal{X}(M_x\otimes I_n)+\kappa (\theta I_N-\mathcal{Y})(I_m\otimes M_y)-(-\lambda +r)I_N,\hfill \end{array}\end{array}$$

(32)

where ⊗ stands for the Kronecker product. The square matrices $M_x$, $M_y$, $M_z$, $M_{xx}$, and $M_{yy}$, are constructed by the associated weights similarly. Additionally, the sparse diagonal matrices $\mathcal{Y}$ and $\mathcal{X}$ are written as:

$$\mathcal{Y}=I_x\otimes \mathrm{diag}\left(y_1,y_2,\cdots ,y_n\right),$$

(33)

$$\mathcal{X}=\mathrm{diag}\left(x_1,x_2,\cdots ,x_m\right)\otimes I_y.$$

(34)

Here the weights corresponding the cross derivative term in the structure of the PIDE (2) can be obtained by employing the Kronecker product as follows:

$$M_{x,y}=M_x\otimes M_y.$$

(35)

Now it is possible to find the following system of ODEs for pricing (2):

$$u^{\prime }\left(\tau \right)=Bu\left(\tau \right).$$

(36)

Now, note that we can use the work of [22,29] to discretize the integral part as follows. By a linear interpolation for $u(x\varrho ,y,\tau )$ among the adaptive numerical grid nodes, the nonlocal integral given in (2) can be solved using

$${\mathcal{A}}_I\left(u\right)={\int }_0^{\infty }u(x\varrho ,y,\tau )b\left(\varrho \right)d\varrho .$$

(37)

Employing $z=x\varrho ,$ one can transform (37) into the integral below:

$${\mathcal{A}}_I\left(u\right)={\int }_0^{\infty }u(z,y,\tau )b\left(\frac{z}{x}\right)\left(\frac{1}{x}\right)dz.$$

(38)

Using a linear interpolation for (38), we can find the following:

$${\mathcal{A}}_i\left(u\right)\simeq \sum _{l=1}^{m-1}{Q}_{i,l},$$

(39)

for every node $x_i$, $i=2,\dots ,m-1$, wherein

$$Q_{i,l}={\int }_{x_l}^{x_{l+1}}\left(\frac{x_{l+1}-z}{\Delta x_l}u(x_l,y,\tau )+\frac{z-x_l}{\Delta x_l}u(x_{l+1},y,\tau )\right)b\left(\frac{z}{x_i}\right)\left(\frac{1}{x_i}\right)dz,$$

(40)

wherein $\Delta x_l=x_{l+1}-x_l$ is the graded step size. Hence, we have

$$\begin{array}{cc}\hfill Q_{i,l}=& \frac{1}{2\Delta x_l}\left(exp\left(\gamma +\frac{{\widehat{\sigma }}^2}{2}\right)\left(-\mathrm{erf}\left(\frac{-ln\left(\frac{x_l}{x_i}\right)+\gamma +{\widehat{\sigma }}^2}{\sqrt{2}\widehat{\sigma }}\right)\right.\right.\hfill \\ \hfill & \left.+\mathrm{erf}\left(\frac{-ln\left(\frac{x_{l+1}}{x_i}\right)+\gamma +{\widehat{\sigma }}^2}{\sqrt{2}\widehat{\sigma }}\right)\right)x_i\left(u\left(x_l,y,\tau \right)-u\left(x_{l+1},y,\tau \right)\right)+\left(\mathrm{erf}\left(\frac{\gamma -ln\left(\frac{x_l}{x_i}\right)}{\sqrt{2}\widehat{\sigma }}\right)\right.\hfill \\ \hfill & \left.\left.\left.-\mathrm{erf}\left(\frac{\gamma -ln\left(\frac{x_{l+1}}{x_i}\right)}{\sqrt{2}\widehat{\sigma }}\right)\right)\left(x_{l+1}u\left(x_l,y,\tau \right)-x_{l}u\left(x_{l+1},y,\tau \right)\right)\right)\right),\hfill \end{array}$$

(41)

wherein $\mathrm{erf}(·)$ stands for the Gaussian distribution.

So, (36) turns into

$$u^{\prime }\left(\tau \right)=\overline{B}u\left(\tau \right),$$

(42)

where $\overline{B}$ is the system matrix after imposing the integral part. Finally, after considering the boundary conditions, a set of ODEs can be attained as follows:

$$\begin{array}{c}{u}^{\prime }\left(\tau \right)=F(\tau ,u\left(\tau \right))=\overline{B}u\left(\tau \right)+b,\hfill \end{array}$$

(43)

wherein b consists of the boundary conditions.

4. The Time-Stepping Solver

Time stepping schemes must be used to solve (43). Although very recently some optimal time stepping solvers have been proposed in literature [30,31,32] for solving system of ODEs, here we focus on a basic but efficient one. Now it is considered that ${\mathbf{u}}^{\iota }$ as an approximation to $u\left({\tau }_{\iota }\right)$, then we could derive our final (explicit) time-integrator method. Select $k+1$ uniform temporal nodes and $0\le \iota \le k$, ${\tau }_{\iota +1}={\tau }_{\iota }+\zeta$, $\zeta =\frac{T}{k}>0$ with ${\mathbf{u}}^0=\left(4\right)$, then the second-order RK solver (RK2) also known as the mid-point explicit method is given by [33] (p. 95):

$${\mathbf{u}}^{\iota +1}={\mathbf{u}}^{\iota }+{\psi }_2+\mathcal{O}\left({\zeta }^3\right),$$

(44)

where

$${\psi }_2=\zeta F\left({\tau }_{\iota }+\frac{1}{2}\zeta ,{\mathbf{u}}^{\iota }+\frac{1}{2}{\psi }_1\right),{\psi }_1=\zeta F({\tau }_{\iota },{\mathbf{u}}^{\iota }).$$

(45)

The approach (44) is useful because its explicit procedure helps programming in lower computational time than many of its competitors from the RK methods. This is a motivation of choosing (44) and not other higher order members of the RK family since their computational cost per time level increases. Anyhow, the investigation for finding the best time-stepping solver from the RK family of an optimal order for our specific PIDE problem remains an open question which could be focused on forthcoming works. Now the most important thing is to investigate that under what conditions this stability can be kept.

Theorem 1.

Let us assume that (43) satisfies the Lipschitz condition, then we have a conditional time-stable iteration process using (44) for solving (43).

Proof. 

Considering the time-stepping method (44) on the set of ODEs (43) gives:

$$\begin{array}{c}\hfill {\mathbf{u}}^{\iota +1}=\left(\frac{{\left(\zeta \overline{B}\right)}^0}{1}+\frac{{\left(\zeta \overline{B}\right)}^1}{1}+\frac{{\left(\zeta \overline{B}\right)}^2}{2}\right){\mathbf{u}}^{\iota }.\end{array}$$

(46)

The explicit method (46) is clearly time-stable if the matrix eigenvalues of

$$\left(I+\zeta \overline{B}+\frac{{\left(\zeta \overline{B}\right)}^2}{2}\right)$$

(47)

have modulus less than or equal to one. Viewing (46) as an iterative map, it would be clear that the eigenvalues of this matrix are $1+\zeta {\overline{B}}_{\mathfrak{i}}+\frac{{\left(\zeta {\overline{B}}_{\mathfrak{i}}\right)}^2}{2},$ where ${\overline{B}}_{\mathfrak{i}}$ are the eigenvalues of matrix $\overline{B}$. Thus, for $\mathfrak{i}=1,2,\dots ,m$, the A-stability is simplified to

$$\left|1+\zeta {\overline{B}}_{\mathfrak{i}}+\frac{{\left(\zeta {\overline{B}}_{\mathfrak{i}}\right)}^2}{2}\right|\le 1.$$

(48)

Therefore, our proposed method is time-stable if the time step size $\zeta$ reads as (48). The stability function in (48) shows a conditional stable behavior for (44). Using (48) along with $\zeta >0$ we have the following:

$$0<\zeta \le \frac{2}{\mathrm{Re}\left({\lambda }_{max}\left(\overline{B}\right)\right)},$$

(49)

where $\mathrm{Re}(·)$ is the real part and ${\lambda }_{max}(·)$ is the largest eigenvalue (in the absolute value sense). Note that we also obtain

$$-{\xi }_{\mathfrak{i}}\le \mathrm{Im}\left({\overline{B}}_{\mathfrak{i}}\right)\le {\xi }_{\mathfrak{i}},$$

(50)

while

$${\xi }_{\mathfrak{i}}={\left(2{\left(-\frac{\mathrm{Re}\left({\overline{B}}_{\mathfrak{i}}\right)(\zeta \mathrm{Re}\left({\overline{B}}_{\mathfrak{i}}\right)+2)}{{\zeta }^3}\right)}^{1/2}-\frac{\mathrm{Re}\left({\overline{B}}_{\mathfrak{i}}\right)(\zeta \mathrm{Re}\left({\overline{B}}_{\mathfrak{i}}\right)+2)}{\zeta }\right)}^{1/2}.$$

(51)

These inequalities on the real and the imaginary parts of the eigenvalues will determine the conditional time stability bounds of the proposed solver when pricing (2). This ends the proof. □

To discuss about the advantages of the proposed approach, we briefly express that our solver has now been expressed all in matrix notations as in (43) which is a system of linear ODEs. When it couples by the ODE solver (44) with the stability condition (50), it solves (2) and the stability relied only on the largest eigenvalue of the system matrix.

5. Numerical Aspects

The goal here is to resolve (2) for at-the-money options, i.e., the value of u at $x_0=K$ and also $y_0=0.04$ and $K=100\$$. The comparing methods are given below:

• The 2nd-order FD scheme with equidistant stencils for space and the explicit 1st order Euler’s scheme denoted by FDM,
• The method of scalable algebraic multigrid discussed in [34] and shown by AFFT.
• The scheme recently proposed by Soleymani et al. in [21] based on efficient non-uniform procedure denoted by SM.
• The presented solver in Section 2, Section 3 and Section 4 shown via RBF-FD-PM in this section.

Noting that all the programs have been written carefully under similar conditions in Mathematica 13 [35,36]. Here, the whole CPU time (for constructing the meshes, the derivative matrices, the set of ODEs and employing the time-stepping method) is in second. We use more number of nodes along x than y, since its computational domain is larger. The criterion given below is used for computing the errors

$${\epsilon }_{i,j,\iota }=\left|u_{\mathrm{approx}}(x_i,y_j,{\tau }_{\iota })-u_{\mathrm{ref}}(x,y,\tau )\right|,$$

(52)

where $u_{\mathrm{approx}}$ and $u_{\mathrm{ref}}$ are the approximate and exact solutions. $u_{\mathrm{ref}}$ is selected from the already well-known literature [14,34].

It is remarked that one efficient way to compute the shape parameter is to calculate it adaptively via the number of discretization points, the numerical domain as well as the structure of the PIDE problem. Hence, here we use ($1\le i\le m-1$):

$$p=4max\{\Delta x_i\},$$

(53)

where $\Delta x_i$ are the increments along the variable mesh. We can write and use (53) similarly for the other variable. Throughout the tables of this paper, a E-b stands for the scientific notation $a\times {10}^{-b}$.

Example 1

([14]). Let us investigate the computational results for the call option of (2) using the following settings:

$$\rho =-0.9,r=0.025,\lambda =0,\sigma =0.3,\kappa =1.5,\theta =0.04,q=0,T=1.$$

(54)

The reference price, which is obtained by the FFT approach [14], is 8.894869 at the point $(x_0,y_0)=(100,0.04)$. The numerical truncated domain is $\mathsf{\Omega }=[0,3K]\times [0,1]$ and $\psi =1.5$. Economically speaking, the values for the variance (for domain truncating) that are larger than one are not significant. The results in this case are provided in

Table 1. Numerical results for Example 1.

Solver	m	n	N	$\mathit{k}+1$	u	$\mathit{\epsilon }$	Time
FDM
	20	20	400	401	8.700	1.94E-1	0.87
	40	25	1000	2001	8.597	2.97E-1	2.33
	40	40	1600	2001	8.673	2.20E-1	4.92
	65	45	2925	4001	8.860	3.39E-2	14.07
	80	50	4000	10,001	8.874	2.03E-2	31.09
AFFT
	10	10	100	251	8.346	5.48E-1	0.41
	15	15	225	251	8.698	1.96E-1	0.54
	25	20	500	401	8.860	3.47E-2	0.56
	30	30	900	601	8.870	2.43E-2	1.49
	50	30	1500	2001	8.885	9.62E-3	4.46
	80	30	2400	5001	8.890	4.32E-3	11.49
SM
	10	10	100	251	8.388	5.06E-1	0.55
	15	15	225	251	8.746	1.47E-1	0.81
	25	20	500	401	8.877	1.71E-2	1.83
	30	30	900	601	8.888	6.09E-3	3.64
	80	30	2400	2501	8.894	8.19E-4	27.19
RBF-FD-PM
	10	10	100	251	8.389	5.05E-1	0.52
	15	15	225	251	5.753	1.41E-1	0.80
	25	20	500	401	8.876	1.88E-2	1.69
	30	30	900	601	8.889	5.56E-3	3.29
	80	30	2400	2501	8.894	6.69E-4	25.74

, which shows the superiority of the proposed solver RBF-FD-PM.

Example 2

([34]). Let us investigate the computational results of a European put option for (2) using the following settings:

$$\begin{array}{c}\gamma =-0.5,\widehat{\sigma }=0.4,\rho =-0.5,\lambda =0.2,\sigma =0.25,\\ r=0.03,T=0.5,\kappa =2.0,\theta =0.04,q=0.\end{array}$$

(55)

The reference prices for specific locations of the domain are 11.302917 at (90, 0.04, 0.5), 6.589881 at (100, 0.04, 0.5) and 4.191455 at (110, 0.04, 0.5) using [34]. The convergence results are provided in

Table 2. Numerical results of the different solvers in Example 2.

Solver	m	n	N	$\mathit{k}+1$	$\mathit{\epsilon }$ at $\mathit{x}=90$	$\mathit{\epsilon }$ at $\mathit{x}=100$	$\mathit{\epsilon }$ at $\mathit{x}=110$	Time
FDM
	8	5	40	26	4.54E0	3.10E0	3.56E-1	0.37
	16	8	128	51	1.23E0	1.31E-1	1.02E0	0.52
	32	16	512	101	1.34E-1	2.77E-1	2.26E-1	1.49
	32	32	1024	201	1.25E-1	2.02E-1	2.23E-1	2.93
	64	32	2048	501	5.35E-3	5.85E-2	1.87E-2	9.47
SM
	8	5	32	51	2.57E-1	7.48E-1	7.13E-1	0.55
	16	8	128	101	7.31E-2	2.51E-1	3.89E-2	0.81
	32	16	512	501	2.18E-3	6.02E-3	5.64E-3	2.40
	64	32	2048	1001	1.79E-3	6.24E-4	1.39E-3	16.71
RBF-FD-PM
	8	5	32	51	2.49E-1	6.89E-1	7.11E-1	0.50
	16	8	128	101	6.25E-2	2.39E-1	3.84E-2	0.76
	32	16	512	501	2.07E-3	5.44E-3	5.04E-3	2.11
	64	32	2048	1001	1.08E-3	5.81E-4	1.04E-3	15.34

and

Table 3. Numerical results of the AFFT solver in Example 2.

Solver	m	n	$\mathit{\epsilon }$ at $\mathit{x}=90$	$\mathit{\epsilon }$ at $\mathit{x}=100$	$\mathit{\epsilon }$ at $\mathit{x}=110$
AFFT	17	9	1.08E0	1.57E0	1.96E-1
	33	17	2.80E-2	5.20E-1	1.38E-1
	65	33	4.78E-3	1.25E-1	2.84E-2
	129	65	7.38E-3	3.09E-2	5.25E-3

and confirm the superiority of the proposed solver with $\psi =2$ in this paper.

The FDM solver is back-of-the-envelope accounting because it is clear that the uniform grids for the PIDE problem are not a fast calculation to obtain highly accurate prices. To check the stability and positivity of the numerical solution for RBF-FD-PM, the numerical solution for Example 2 is plotted in Figure 1, which shows the stable behavior of RBF-FD-PM using $m=16$, $n=8$ and $k=1001$.

The reason for providing Figure 1 is twofold. We must first reveal that the numerical solution obtained by RBF-FD-PM using some m and n is stable and does not have oscillations. This is important since the PIDE model has a mixed derivative term, which can lead to oscillations in the numerical solution as long as a careless numerical method is employed. Second, we must reveal how the graded meshes (the green points in Figure 1) located on the numerical solution are obtained by employing an automatic interpolation on the obtained solutions.

An inquiry might arise by analyzing the results in Table 1 and Table 2. tt is not easy to find out the advantages of the proposed approach since the numerical values are given for different values of the parameters m, n, N and $k+1$. In fact, larger time step sizes (lower k) are taken for SM and RBF-FD-PM since their ODE solver, i.e., (44), has a larger stability region, and the overall solvers must be compared by fixing an accuracy for the errors and then checking the computational times.

To also show how the instability may ruin the numerical pricing using the stability bound (49), we provide the numerical results of solving (2) by the RBF-FD-PM using $m=16$ and $n=8$, but with $k=25$ uniform discretization nodes along time in Figure 2. This shows that all the involved solvers have some limitations, but the proposed solver sounds more efficient than others.

However, due to nonsmoothness at the strike price in the initial condition, it might be useful to employ a time-stepping solver that works on graded meshes over time with more focus at the beginning of the starting time, i.e., zero (the solution near the initial time point has a weak singularity). One such method is the Rannacher time-marching method [37]. Although such an application will help our solver a lot, we will try to focus on this in forthcoming related works.

6. Concluding Remarks

PIDEs arise in the mathematical modeling of many processes in different fields of engineering and finance. This paper has presented an approximate solution of the linear Bates PIDE with clear application in financial option pricing using a local integral term. The solution method was considered on graded meshes at which there is a clear concentration of the discretization nodes on the financially important are of the problem. Then, an RBF-FD solver using semi-discretization via sparse arrays have been constructed for solving the Bates PIDE. The numerical results were furnished and supported the theoretical discussions. These results have been provided in Table 1 and Table 2 which implicitly state that the proposed approach can compete the most efficient solver (SM) for the same purpose. Additionally, the prospects for future research can be focused on how to obtain RBF-FD weights on stencils having five/six adjacent nodes on graded meshes or employing the Rannacher time-marching method in order to obtain higher accuracies for solving the PIDE problem (2).