A numerical study of the European option by the MLPG method with moving kriging interpolation

In this paper, the meshless local Petrov–Galerkin (MLPG) method is applied for solving a generalized Black–Scholes equation in financial problems. This equation is a PDE governing the price evolution of a European call or a European put under the Black–Scholes model. The θ-weighted method and MLPG are used for discretizing the governing equation in time variable and option pricing, respectively. We show that the spectral radius of amplification matrix with the discrete operator is less than 1. This ensures that this numerical scheme is stable. Numerical experiments are performed with time varying volatility and the results are compared with the analytical and the numerical results of other methods.

numerical solution for valuing options on dividend-paying stocks. The boundary condition is enforced to take into account the fact that the stock price will drop by the amount of the dividend (Schwartz 1977). A numerical upwind scheme for solving the backward time parabolic partial differential equation is proffered by Vazquez (1998). A robust finite difference method for pricing European and American options is presented by Le (2010, 2011). Cen et al. introduced a central difference scheme with moving mesh in the spatial discretization for pricing Asian options (Cen et al. 2013). Lesmana and Wang purposed an upwind finite difference method for a nonlinear Black-Scholes equation (Lesmana and Wang 2013). The uniform cubic B-spline collocation method is implemented to find the numerical solution by a horizontal method of lines, to discretize the temporal variable and the spatial variable (Kadalbajoo et al. 2012). Kumar et al. suggested that the governing equation is discretized by the θ-weighted method and the option price is approximated by the radial basis functions based on the finite difference method (Kumar et al. 2015). Mohammadi presented Quintic B-spline collocation approach for solving generalized Black-Scholes equation governing option pricing, the horizontal method of lines for time integration and θ-method are used for temporal discretization (Mohammadi 2015). The fractional Black-Scholes model are presented by Bjork and Hult (2005), Song and Wang (2013). Kumar et al. presented a numerical computation of fractional Black-Scholes equation arising in financial market (Kumar et al. 2014). The meshless local Petrov-Galerkin (MLPG) method based on moving kriging interpolation for solving fractional Black-Scholes model are purposed by Phaochoo et al. (2016). Kleinert and Korbel presented double-fractional differential equation for the prices of options (Kleinert and Korbel 2016). Zhang et al. introduced the numerical simulation of the tempered fractional Black-Scholes equation for European double barrier option (Zhang et al. 2016).
A truly meshless method for solving boundary value problems based on the local symmetric weak form and the moving least squares (MLS) approximation, called the MLPG method, has been proposed and successfully applied to solve many problems (Atluri and Shen 2002). The MLPG method is one of the most viable methods in which the moving least square (MLS) approach is used to construct the shape functions. Although the MLPG method has been applied to many problems, there exists an inconvenience or disadvantage when using the MLPG because of the difficulty in implementing essential boundary conditions. This is because the MLS shape functions lack the Kronecker delta property. Therefore the moving kriging interpolation (MKI) method (Yimnak and Luadsong 2014) has been proposed to overcome this problem. It uses nodal values in the local support domain to construct shape functions with the Kronecker delta property. The MKI method works well for practical problems.
In this paper, we propose a numerical method based on the MLPG method to solve a generalized 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. MLPG type 2 (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. In addition, we compared numerical solutions among the finite difference method, the cubic spline method and the MLPG method in the numerical experiment.

Problem formulation
The Black-Scholes equation is the outstanding financial equation that solves European option pricing without a transaction cost. Moreover, underlying asset prices distributed on the log indicate normal random walks, risk-free interest rates, no dividends and no arbitrate opportunities as fundamental assumptions. The Black-Scholes equation follows: with terminal and boundary conditions where u(s, τ), T, r and σ are the value of European call options at underlying asset price s at time τ, the expiration date, risk-free interest rate and volatility of underlying asset prices, respectively. From Eq. (1), when s goes to zero then degenerating will occur in approximation. We transform the Black-Scholes equation into a non-degenerate partial differential equation by using a logarithmic transformation, x = ln s, t = T − τ, and define the computational domain for convenience in numerical experiments by (Huang and Cen 2014).

Moving kriging interpolation
The moving kriging interpolation (MKI) is used to construct the shape function. The function u(x) is defined in the domain Ω and the approximate function is u n (x). The subdomain Ω x that encompasses these surrounding nodes is called the interpolation domain of point x. The formulation of the meshless shape function using MKI is given by T is a vector value of the function in the domain Ω. Φ(x) is a 1 × N vector of shape functions, expressed as where matrix A and B are defined as and For matrix P with size N × m, values of the polynomial basis function Eq. (7) at the given set of nodes are collected as follows Matrix R and vector r(x) are defined by the following where γ(x i , x j ) is the correlation between any pair of nodes located at x i and x j , representing the covariance of the field value u(x). A simple and frequently used correlation function is a Gaussian function as where r ij = x i − x j and ǫ > 0 are the correlation and shape parameters, respectively used to fit the model.

Spatial discretization
The MLPG method constructs the local weak form over the local subdomain, which is a small region taken for each node in a global domain. Multiplying test function v i (x) into Eq. (2) and then integrating it over subdomain i s yields the following expression: where v i is a test function that is significant for each node. Rearranging Eq. (12), we have where N is the number of nodes surrounding point x which has an effect on u(x) and û j (t) is value of the option at time t. The shape function, ф j (x), is constructed by the moving kriging interpolation which has the Kronecker delta property, thereby enhancing arrangement of the nodal shape construction accuracy. Rearranging Eq. (14) yields the following results: This research uses MLPG type 2 (MLPG2) (Cen and Le 2011), then the test function v i is chosen by the Kronecker delta function, The test function will define significance for each node in the subdomain. In this case, substituting test function v i (x) to Eq. (15) and then integrating it over subdomain i s yields the following results: Equation (16) can be written in matrix form as follows: Since the shape function that is constructed by the moving kriging interpolation satisfies the Kronecker delta property, A is the identity matrix. Therefore, Eq. (17) can be written as

Temporal discretization
The numerical solution of a European option, using the implicit method, requires the generation of a modified PDE operator through a finite difference approximation of time derivative. We will do this using the θ-weighted method. Consider the following initial-boundary value problem: . By a finite difference approximation made for the time derivative with notation u n (x) that approximates the exact solution u(x, t) at t n and t n = t n−1 + Δt, we obtain For each fixed time level t n , the above equation, Eq. (20), is the system of linear ODEs. Now, using the MLPG2 method with moving kriging interpolation for constructing shape functions is discussed through Eqs. (12)-(18) for spatial discretization of operator Lu leads to where U n = û n 1û n 2û n 3 . . .û n N T , B is the discretization matrix for the space discretization of linear differential operator Lu, and I is the identity matrix.

Stability analysis
In this section, we present an analysis of the stability of the MLPG2 method with moving kriging interpolation using the matrix method. A small fluctuation at the nth time level e n = U n −Ũ n is introduced in Eq. (21), where U n is exact and Ũ n is the numerical solution. The equation of the error e n+1 can be written as e n+1 = Ge n , where the amplifi- The numerical scheme will be stable if as n → ∞, the error e n → 0. This can ensure that ρ(G) < 1 provided, where ρ(G) denotes the spectral radius of G.
It can be seen that stability is assured if all eigenvalues of the matrix [I − θ∆tB] −1 [I + (1 − θ)∆tB] satisfy the following conditions: where λ is the eigen value of matrix B. In the case of the Crank-Nicolson scheme (θ = 1 2 ) the inequality in Eq. (22) is satisfied when Re(λ) max ≤ 0. This shows that the scheme is stable if Re(λ) ≤ 0. The eigenvalues of matrix B highly depends on the mesh spacing parameter h and the shape parameter ǫ, where h is defined to be the minimal distance between any two points in the domain. Since it is not possible to find an explicit relationship among the eigenvalue of matrix B, the number of nodes and the shape parameter ǫ we investigated this dependent numerically, as is given in Fig. 1. Figure 1 shows that the maximum eigenvalue Re(λ) of matrix B varies as a function of shape parameter ǫ, when mesh spacing parameter h is constant. Figure 2 shows the effect of mesh length, h, for eigenvalue of matrix B, when the shape parameter ǫ is constant. It is found that the condition number of the collocation matrix becomes very large and the system leads to ill-conditioning, when ǫ and h become very small. Figure 3 shows that increasing volatility trends to ill-conditioning. In this case, if the shape parameter increases then eigenvalue Re(λ) will decrease. Figure 4 presents that a risk-free interest rate which decreases leads to ill-conditioning and if the shape parameter increases, then Fig. 1 The relation between Re(λ) max and shape parameter (ǫ)

Fig. 2 Relation between Re(λ) max and spatial mesh length (h)
eigenvalue Re(λ) will decrease. Figure 5 shows the relation between mesh length h and the smallest of shape parameter ǫ is Re(λ) max < 0. It shows that ǫ could be large when h becomes smaller in value. It can be seen that we can control the stability of this numerical scheme by choosing the appropriate shape parameter.

Numerical experiments
In this section, we are going to present various numerical results to evaluate proposed meshless approaches. Although the schemes work for all correlation functions, we will use the Gaussian function on different experimental setups. Using the MLPG2 method, the resulting problems for European call options are solved via Crank-Nicolson's method. The computational domain is partitioned with N being equi-spaced spatial nodes with mesh lengths. The temporal domain is divided into K equi-spaced points.
The relation between the smallest of shape parameter ϵ and the mesh length 'h' To illustrate accuracy of the proposed method, numerical simulation was done for the European call option with parameters σ = 0.4, r = 0.8, T = 1, E = 1, x min = − ln(4E), x max = ln(4E). Accuracy is measured in the discrete maximum norm and root mean square error. The discrete maximum norm and maximum of root mean square error are given in Table 1 that is estimated for differences N and K in the MLPG2 methods.
where û i is an approximate solution of option price s and u i is the analytic solution of option price s i .
Tables 1 and 2 show convergence trends of the present method, with N = K where N and K are the number of points in the spatial and temporal domain. From these tabular results one can observe that Crank-Nicolson's method converges to the exact solution and the maximum error and root mean square error decreases while increasing the number of nodes for any ǫ.

The European call option for long maturity time
We consider the Black-Scholes equation for the European call option with parameters r, E, T with a wide range of volatilities σ and strike price E. Figures 6 and 7 shows the relation between maximum error and ǫ for varieties of volatility. Figure 6 presents that ǫ increases and then the error tends to stabilize for any volatility. In Fig. 7, we see that the error will increase when the volatility become small and shape parameter, ǫ becomes large. The other volatility tends to slowly increase the error. Figures 8 and 9 shows the relation between maximum error and ǫ for a variety of riskfree interest rates. These figures are observed at large values of risk-free interest rates due to the convection-dominant nature of the problem.

The European call option with volatility function
We consider the Black-Scholes equation for the European call option with parameters r = 0.06, The volatility function is the same as the given in Kadalbajoo et al. (2012) and Kumar et al. (2015). In this case, the exact solution is not known. Table 3 shows that the comparison of value options from finite differences, the cubic spline and MLPG method. We found that these results are very similar. Table 4 shows that our value options results differ by only 0.54 % as compared to the finite difference method and 0.36 % as compared to the cubic spline method.

Conclusion
In this paper, the MLPG method was proposed for the Black-Scholes equation, which is transformed for a non-degenerate partial differential equation, and used moving kriging shape functions which have Kronecker delta properties. The temporal discretization was chosen by the Crank-Nicolson method. The eigenvalue (λ) of B depends on the mesh spacing parameter (h) and the shape parameter (ǫ  the volatility (σ) is an increasing function for all shape parameters. The relation between Re(λ) max and risk free interest rates (r) tends to decrease functions for all shape parameters. If the mesh length (h) increases, then the smallest of shape parameters ǫ will be decreased. The numerical results have demonstrated the accuracy and efficiency of the present methods. The present method gives the value option in both regular and irregular nodal points. We found that the relation between errors and shape parameters vary by volatility and risk-free interest rates. This method works well for finding the approximate solution of option pricing with a volatility function.