The meshless local Petrov–Galerkin based on moving kriging interpolation for solving fractional Black–Scholes model

1 Introduction

The first idea of fractional calculus is considered to be the Leibniz’s letter to L’Hospital in 1695. Fractional calculus is a name for the theory of derivatives and integrals of arbitrary order. The famous definitions of a fractional calculus are the Riemann–Liouville and Grunwald–Letnikov definition (Ghandehari and Ranjbar, 2014). Caputo reformulated the definition of the Riemann–Liouville in order to use integer order initial conditions to solve fractional differential equation (Ishteva, 2005). The definitions of Riemann–Liouville and the first Caputo version have the weakness for singular kernel. Caputo and Fabrizio proposed a new fractional order derivative without a singular kernel (Atangana and Alkahtani Badr Saad, 2015a,b). Fractional differential equations have attracted much attention during the past few decade. This is the fact that fractional calculus supplies an competent and excellent tool for the description of many important phenomena such as electromagnetic, physics, chemistry, biology, economy and many more.

Black–Scholes equation, which is proposed by Black and Scholes (1973), is the financial model that concern with option. An option is a contract between the seller and the buyer. It consists of a call option and a put option. Option valuation depends on the underlying asset price and time. The European option can only be exercised at the expiration date, but the American option can be exercised at any time before expiration date. The solution of Black–Scholes equation provides an option pricing formula for European option. The analytic solution is used in general case with basic assumption but it is not satisfied in some conditions. Some restrictions were appeared in the classical Black–Scholes equation that is the weakness of this model. Original assumptions were relieved by other models such as models with transaction cost (Barles and Soner, 1998; Davis et al., 1993), Jump-diffusion model (Merton, 1976), Stochastic volatility model (Hull and White, 1987) and fractional Black–Scholes model (Bjork and Hult, 2005; Wang, 2010).

Fractional Black–Scholes model is derived by many researchers. Some restrictions were appeared in the classical Black–Scholes equation that is the weakness of this model (Song and Wang, 2013). The fractional Black–Scholes models are derived by substituting the standard Brownian motion with fractional Brownian motion.

In this paper, we propose a numerical method based on meshless local Petrov–Galerkin (MLPG) method to solve a fractional Black–Scholes equation. The MLPG is a truly meshless method, which involves not only a meshless interpolation for the trial functions, but also a meshless integration of the weak-form, (Atluri and Shen, 2002). MLPG2 is chosen for this research so the Kronecker delta is the test function. This method will avoid the domain integral in the weak-form.

2 Problem formulation

The Black–Scholes equation is the outstanding financial equation that solves the European option pricing without a transaction cost. Moreover, underlying asset price distributed on the lognormal random walk, risk-free interest rate, no dividend and no arbitrate opportunity are fundamental assumption. Fractional calculus is used in financial market for description the probability of log-price, which is a benefit to specify the variability in prices. Fractional Black–Scholes models are derived by substitution standard Brownian motion with fractional Brownian motion. Some restrictions of original Black–Scholes equation are improved by many methods. The fractional Black–Scholes equation is a choice for reinforcement of original equation. The fractional Black–Scholes equation is following

From Eq. (2.1), when s goes to zero then degenerating will occur in approximation. We transform the Black–Scholes equation into a nondegenerate partial differential equation by using a logarithmic transformation x = ln s, t = T − τ, and define the computational domain for convenient in numerical experiments by Ω = [x_min, x_max] × [0, T], where x_min = −ln (4E), x_max = ln (4E) , Haung and Cen, 2014.

The fundamental definition of fractional calculus as following

3 Spatial discretization

The MLPG method are used for spatial discretization. We create the local weak form over local subdomain, which is a small region taken for each node in global domain. Multiplying test function v_i into Eq. (2.2) and then integrate over subdomain ${\mathrm{\Omega }}_s^i$ yields the following expression

The test function will define significance for each node in subdomain. In this case, substituting test function v_i(x) to Eq. (3.4)and then integrate over subdomain ${\mathrm{\Omega }}_s^i$ yields the following result

Eq. (3.5) can be written in the matrix form as following

4 Temporal discretization

The numerical solution of European option use the implicit finite difference method. By a finite approximation made for the time fractional derivative with notation $\frac{{\partial }^{\alpha }u(x_i\text{,}{t}_n)}{\partial t^{\alpha }}$ that approximates the exact solution u(x_i, t_n) at time level n, we restrict attention to the finite space domain x_min < x < x_max with 0 < α < 1. The time fractional derivative use the implicit finite difference (Murio, 2008), defined by

Hence, $\frac{d^{\alpha }U}{{\mathit{dt}}^{\alpha }}=D_t^{(\alpha )}{U}_i^n+O(\mathrm{\Delta }t)$ .

The first-order approximation method for the computation of Caputo’s fractional derivative is given by

Consider the Eq. (4.2) and substitute time fractional derivative that following $${\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=1}^n}{\omega }_j^{(\alpha )}(U^{n-j+1}-U^{n-j})+{\mathit{BU}}^n=0\text{,}$$ $${\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=1}^n}{\omega }_j^{(\alpha )}{U}^{n-j+1}-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=1}^n}{\omega }_j^{(\alpha )}{U}^{n-j}+{\mathit{BU}}^n=0\text{,}$$ $${\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}(U^n-U^{n-1})=-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}(U^{n-j+1}-U^{n-j})-{\mathit{BU}}^n\text{,}$$ $${\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}(U^n-U^{n-1})=-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}(U^{n-j+1}-U^{n-j})-{\mathit{BU}}^n\text{,}$$ We consider the first case for n = 1,

Case n = 1 $${\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{U}^1+{\mathit{BU}}^1={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{U}^0\text{,}$$ $$({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)U^1={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{U}^0\text{.}$$ Case $n⩾2$ $$\left({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B\right)U^n={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{U}^{n-1}-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}(U^{n-j+1}-U^{n-j})\text{.}$$

5 Stability analysis

In this section, we are propose an analysis of the stability of implicit finite difference method and MLPG2 by using the matrix method. A small perturbation at nth time level is $e^n=U^n-{\tilde{U}}^n$ , where Uⁿ is exact and ${\tilde{U}}^n$ is the numerical solution. The equation of the error eⁿ can be written as

If $n=1$ then $({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)e^1={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^0\text{,}$ $$e^1={({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)}^{-1}{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^0\text{.}$$ where $G={({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)}^{-1}{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}$ .

It can be seen that the stability is assured if all eigenvalues of the matrix $G={({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)}^{-1}{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}\text{.}$ satisfy the following condition:

For case n > 1, $$({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B)e^n={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^{n-1}-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^{n-1}}{\omega }_j^{(\alpha )}(e^{n-j+1}-e^{n-j})-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_n^{(\alpha )}(e^1-e^0)\text{.}$$ $$={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^{n-1}-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^{n-1}}\left({\omega }_{j+1}^{(\alpha )}-{\omega }_j^{(\alpha )}\right)e^{n-j}+{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_n^{(\alpha )}{e}^0\text{.}$$ $$⩽\left({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}-{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\displaystyle \sum _{j=2}^{n-1}}\left({\omega }_{j+1}^{(\alpha )}-{\omega }_j^{(\alpha )}\right)+{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_n^{(\alpha )}\right)e^0\text{.}$$ $$={\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^0\text{.}$$ $$e^n={\left({\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}I+B\right)}^{-1}{\sigma }_{\alpha \text{,}\mathrm{\Delta }t}{\omega }_1^{(\alpha )}{e}^0\text{.}$$ $$e^n={\mathit{Ge}}^0\text{.}$$ where ‖G‖ < 1. $$e^n⩽e^0\text{.}$$ From case n = 1 and n > 1, we conclude that $e^n⩽e^0$ for all n where $\text{Re}{(\lambda )}_{\text{max}}>0$ .

The eigenvalues of matric B highly depend on the mesh spacing parameter ‘h’ (‘h’ is defined to be the minimal distance between any two points in the domain) and the shape parameter ε. The present local approximation is free from these complexities. Since it is not possible to find and explicit relationship among the eigenvalues of matrix B, the number of node and the shape parameter ε we investigate this dependent numerically and is given in Fig. 1.

Close

Figure 1

The relation between $\text{Re}{(\lambda )}_{\max }$ and shape parameter $(\in )$ .

Export to PPT

Fig. 1 shows maximum eigenvalue $\text{Re}(\lambda )$ of matrix B varies as a function of shape parameter ε, when mesh spacing parameter h is constant. Fig. 2 shows the effect of mesh length ‘h’ for eigenvalue of matrix B, when the shape parameter $\in$ is constant. Fig. 3 shows that the increasing of volatility trends to decreasing of maximum eigenvalue at the end. In this case, if the shape parameter increase then eigenvalue $\text{Re}(\lambda )$ will increase. Fig. 4 presents that the risk free interest rate is stable for all case of shape parameter and the shape parameter increase then eigenvalue $\text{Re}(\lambda )$ will increase (see Figs. 5–14).

Close

Figure 2

The relation between $\text{Re}{(\lambda )}_{\max }$ and spatial mesh length (h).

Export to PPT

Close

Figure 3

The relation between $\text{Re}{(\lambda )}_{\max }$ and the volatility (σ).

Export to PPT

Close

Figure 4

The relation between $\text{Re}{(\lambda )}_{\max }$ and risk free interest rate (r).

Export to PPT

Close

Figure 5

The approximate solution compare with the exact solution for $\sigma =0.2\text{,}r=0.04\text{,}\alpha =0.5\text{,}t=T$ .

Export to PPT

Close

Figure 6

The approximate solution compare with the exact solution for $\sigma =0.2\text{,}r=0.04\text{,}\alpha =0.5\text{,}0⩽t⩽T$ .

Export to PPT

Close

Figure 7

The approximate solution compare with the exact solution for σ = 0.2, r = 0.04, α = 0.99, t = T.

Export to PPT

Close

Figure 8

The approximate solution compare with the exact solution for $\sigma =0.2\text{,}r=0.04\text{,}\alpha =0.99\text{,}0⩽t⩽T$ .

Export to PPT

Close

Figure 9

The approximate solution compare with the exact solution for $\sigma =0.2\text{,}r=0.1\text{,}\alpha =0.99\text{,}t=T$ .

Export to PPT

Close

Figure 10

The approximate solution compare with the exact solution for $\sigma =0.2\text{,}r=0.1\text{,}\alpha =0.99\text{,}0⩽t⩽T$ .

Export to PPT

Close

Figure 11

The comparison of the approximate solutions of the fractional and standard Black–Scholes equation for $r=0.06\text{,}\alpha =0.99\text{,}t=T$ .

Export to PPT

Close

Figure 12

The comparison of the approximate solutions of the fractional and standard Black–Scholes equation for $r=0.06\text{,}\alpha =0.99\text{,}0⩽t⩽T$ .

Export to PPT

Close

Figure 13

The comparison of the approximate solutions of the fractional and standard Black–Scholes equation for $r=0.06\text{,}\alpha =0.99\text{,}t=T$ .

Export to PPT

Close

Figure 14

The comparison of the approximate solutions of the fractional and standard Black–Scholes equation for $r=0.06\text{,}\alpha =0.99\text{,}0⩽t⩽T$ .

Export to PPT

6 Numerical examples

In this section, we are going to present various numerical results to evaluate proposed meshless approaches. Using the MLPG2 method, the resulting problems for European call options are solved via implicit finite difference method.

The European call option can be modeled by Black–Scholes PDE

To illustrate accuracy of proposed method numerical simulation was done for European call option with parameters x_min = −ln (4E), x_max = ln (4E).

7 Conclusion

In this paper, the fractional Black–Scholes equation are solved by the implicit finite difference method and MLPG2 for discretizing in time variable and option price, respectively. The stability analysis present relation between maximum eigenvalues of matrix and variety of parameters.

The numerical results are presented in two examples. Example 1 presents numerical results for varieties of parameters in four cases. For case 1 and case 2 present comparison of numerical result for different volatility σ = 0.1, 0.2 while variables are fixed, α = 0.5, 0.99 show for case 2 and case 3. Difference of risk free interest rate present for case 3 and case 4, respectively. Example 2 shows the value of option for volatility function in both cases. Volatility function is σ = 0.15(0.5 + 2τ)((s/100 − 1.2)²/((s/100)² + 1.44)) and σ = 0.4(2 + sin x), respectively. Moreover, we found that the MLPG give the value option in both volatility constant and volatility function.