Qualitatively Stable Nonstandard Finite Difference Scheme for Numerical Solution of the Nonlinear Black–Scholes Equation

In this paper, we use a numerical method for solving the nonlinear Black–Scholes partial diﬀerential equation of the European option under transaction costs, which is based on the nonstandard discretization of the spatial derivatives. The proposed scheme, in addition to the unconditional positivity, is stable, consistent, and monotone. In order to illustrate the eﬃciency of the new method, numerical results have been performed by four models.


Introduction
Financial mathematics is a branch of applied mathematics that deals with financial markets. Option contracts are the most traditionally financial instruments. Hence, with the increasing demand for this financial instrument, the pricing of contracts has very important role to play in the economy. An option on an underlying asset traded on financial markets is an agreement that gives its holder the option to buy or sell the asset mentioned in the contract on a specified date at a certain and predetermined price. e date specified in the contract also calls to the maturity date or expiration date, and the price specified in the contract is called the strike price or exercise price. e biggest advantage of an option contract is limited loss and unlimited profits. Option contracts are divided into two types such as European and American types according to the time of exercise. A European option is only applicable at the time specified in the contract, but the American option is applicable at any time until maturity date t ≤ T. One of the models in the discussion of financial derivatives pricing is the Black-Scholes model. is model was introduced by Fisher Black and Miran Scholes and developed by Robert Merton [1]. It is represented by a partial differential equation in the following form: where U � U(S, t) is the value of the option, S is the current underlying asset price, t t � t is the time, T > 0 is the expiration date of the option, r ≥ 0 is the risk-free interest rate, and σ 0 is the volatility (also called implicit volatility). r and σ 0 are fixed. Partial differential equation paid attention in the field of economics with the Black-Scholes model. erefore, in order to obtain the price of these contracts, we need a simple and accurate method to solve this differential equation. e standard Black-Scholes equation is considered in a complete market under certain simple assumptions such as the cash market, without transaction cost. Also, the most important feature in the classical Black-Scholes model, the volatility of an underlying asset is assumed to be constant. In 1987, all over the world, stock the value decreased observably in the short course of time. After this decrease, analysts have concluded that the options with profits and losses are traded with higher implicit volatility versus the option that do not have profit and loss. So, this fact is not compatible with the assumption of constant volatility in the Black-Scholes model. In addition to despite the transaction cost, the Black-Scholes equation changes to an equation with a nonlinear volatility function, which depends on the time, current underlying asset price, and the price derivatives of the option relative to the underlying asset price. en, the classic Black-Scholes equation changes to the following nonlinear Black-Scholes equation: in which σ is the modified volatility. Various nonlinear volatility models have been presented in [2][3][4]. We consider the nonlinear Barles and Soner model which is obtained as follows: where a � μ ��� cN, in which c is the risk-adjusted factor and N is the number of options that are sold. e function Ψ is the solution of the following nonlinear initial value problem [5]: For function Ψ, an exact implicit solution is presented as follows [6]: Let us consider equation (2) with the nonlinear modified volatility of Barles and Soner (3) and (4). By applying the variable τ � T − t, equation (1) is converted to the following equation: for S > 0 and 0 < τ ≤ T, as in [5] we consider the initial condition, the first boundary condition, and the second boundary condition for S ∈ (0, S max ) and τ ∈ (0, T] as follows: for the Call option, max(K − S, 0), for the Put option, where K, K 1 , K 2 , and K 3 denote the strike prices of the options, B is a positive constant, and H is the Heaviside function. Conditions L 1 (0) � L 2 (0) and L 1 (S max ) � L 3 (0) are established. So far, no classical solution has been given for equation (2), but the viscosity solution is used as the solution of equation (2). e numerical solution of the nonlinear Black-Scholes equation has been considered in some works. For example, in [7], an explicit method has been presented that preserves positivity property and is unconditional stable but not consistent under any conditions. A variational principle is established for the generalized KdV-Burgers equation by the semiinverse method in [8][9][10][11]. In [12][13][14], the homotopy perturbation method is designed to obtain a solution to the Black-Scholes equation with boundary conditions for a European option pricing problem and nonlinear problems. Methods based on the variational iteration method (VIM) in [15][16][17] are proposed for the Black-Scholes model. In [18], an explicit nonstandard finite difference method has been proposed that it is stable and conditionally acceptable. Also, a finite difference method has been proposed which is the stability condition for it and is established by limiting the length of the step time and the spatial step size. In [5], a completely implicit time step method has been presented for (6) and represented convergence approximate solution to the viscosity solution and used Newton's algorithm to solve a nonlinear system. ey caused numerical instability in the solution that has been resolved by an iterative method. is process leads to very expensive computing costs. Also in [6], the MOL method has been used to solve (6). In [19,20], equation (6) is solved by the variable x � ln(S/K) which in addition to having very expensive computational cost, convergence, and other results of this work is under relatively unrealistic conditions or limits the step-size time and space.
In [21], considering a class of nonlinear option pricing models, Newton and Picard iterative procedures are proposed for solving the nonlinear systems of algebraic equations. Also, to improve the computational fast, two-grid algorithms are developed. In [22], combination of a sixthorder finite difference scheme in space and a third-order strong stability preserving Runge-Kutta in time has been used to obtain numerical solutions of the linear and nonlinear European put option models. In [23], a new fourthorder finite difference approximation is contributed and equipped with the fourth-order Runge-Kutta scheme for the nonlinear Black-Scholes. ey proved that under several criteria, the procedure is time stable. Also, the proposed technique reduces the computational cost when more accurate results are requested. In [24], nonlinear generalization of the Black-Scholes equation for pricing and American-style call options is investigated. ey proposed a numerical method that involves transforming the free boundary problem for a nonlinear Black-Scholes equation into the so-called Gamma variational inequality with a new variable depending on the Gamma of the option. Also, a modified projected successive overrelaxation method in order to construct a numerical scheme for discretization of the Gamma variational inequality is presented.
In this paper, we solve equation (6) by using an explicit method that has very low computational cost, and the method is free of any spurious oscillations. It is unconditionally positive, stable, consistent, and monotone. We show that the solution of the proposed method converges to the solution of the viscosity of the equation. e rest of the paper is organized as follows. In section 2, we give a brief summary of the NSFD methods. In section 3, we examine the application of nonstandard methods through spatial nonlocalized discretization for the nonlinear Black-Scholes under transaction costs. In section 4, we analyze the new method. e numerical results are given in section 5.

Nonstandard Finite Difference Scheme
In this section, we give some preliminaries including nonstandard finite difference methods. e initial foundation of NSFD schemes came from the exact finite difference schemes. Numerical methods based on the standard finite difference approach are consistent with the original differential equation and guarantee convergence of the discrete solution to the exact one, but they impose severe restrictions on the time step, and essential qualitative properties of the solution are not transferred to the numerical solution. One response to this situation was the initiation by Mickens [25] of a research program for the investigation of new methods for constructing finite difference schemes that are convergent for any step size. Nonstandard finite difference methods (NSFDs) in addition to the usual properties of consistency, stability, and hence convergence produce numerical solutions that also exhibit essential properties of solutions [3-5, 13, 15, 16, 18, 21, 26-38].
is class of schemes and their formulations centers on two issues: first, how should discrete representations for derivatives be determined and second what are the proper forms to be used for nonlinear and reaction terms. e forward finite difference method is one of the simplest discretization schemes. In this method, the derivative V x is replaced by However, in the Mickens schemes, this term is replaced by is an increasing continuous function of h, and the function ϕ(h) satisfies the following conditions: Note that in taking lim h ⟶ 0 to obtain the derivative, the use of any of this ϕ(h) will lead to the usual result for the first derivative: where ϕ 1 (h) and ϕ 2 (h) are continuous functions of the step size h. A scheme is called a nonstandard finite difference method if at least one of the following conditions is met: (i) In the discrete derivatives, the traditional denominator h is replaced by a nonnegative function ϕ(h) such that or Examples of functions ϕ(h) that satisfy these conditions are as follows [25]: (ii) Nonlinear terms are approximated in a nonlocal way, i.e., by a suitable function of several points of the mesh. If there are reaction terms of the differential equation [3-5, 13, 15, 16, 18, 21, 26-28], these are replaced by

Journal of Mathematics
One can say that there is no appropriate general method to choose the function ϕ(h) or to choose which nonlinear terms are to be replaced, and some special techniques may be found in [3-5, 13, 15, 16, 18, 21, 26-28].

Construction of the New NSFD
In this section, by using the mentioned rules for NSFD schemes, we propose our new NSFD scheme. Let us consider the computational domain Ω � [0, S max ] × [0, T] and discretize it in the following form. We introduce a grid of mesh points (S, t) � (S j , t n ) where S j � jh and t n � nΔt, j � 0, 1, . . . , M; n � 0, 1, . . . , X, and the spatial step-size is given by h � S max /M, and the time step size is Δt � T/X. We denote the approximation of V(S j , t n ) by V n j . e new NSFD is of the following form: e explicit form of (17) is as follows: e finite difference approximation provides a difference equation as follows: with (8)-(10), and we define the initial and boundary conditions as follows: e matrix form (18) in which A n (V n ) is a tridiagonal matrix with positive components and B n is known as boundary values. e new scheme can be written in the following form: where

Analysis of the New Method
In this section, we investigate the monotonicity, positivity preserving, stability, consistency, and convergency properties of the new method.

Journal of Mathematics
where v n j � V SS (j), K n+1 j � e rτ n+1 a 2 S 2 j > 0, (33) for each j � 0, 1, . . . , M − 1 and K > 0; independent of v, by derivative phrase (1 + Ψ(Kv)) relative to the variable v and using (5) and (24) and w � Kv, we have So, (1 + Ψ(Kv))v is an increasing function of v. Using the above symbols and by monotone combination of (1 + Ψ(Kv))v and properties of the linear sentences of the relationship (26) and (30), we have is is the same as (28). Similarly from monotonicity (1 + Ψ(Kv))v and (31), it can be shown that (29) is confirmed, and this completes the proof. Now, we investigate the positivity property of the constructed methods in the previous section. With positivity, we mean that the component-wise non-negativity of the initial vector is preserved in time for the approximated solution. Several NSFD schemes have been constructed in the literature (see, for example, [1,25,29,[39][40][41][42]). □ Proposition 1 (positivity). Suppose V n j−1 , V n j , and V n j+1 are the real numbers that are nonnegative, in this case, method (18) offers a nonnegative approximation of V n+1 j for the solution of equation (2).
Proof. Since the parameter also then x n j , y n j , z n j ≥ 0.
erefore, this completes the proof.
Proof. As in [5], for each n � 0, 1, . . . , M − 1, we know V n+1 � (V n 0 , (V n ) T , V n M ), that V n+1 is the solution of the NSFD method, and then for the initial and boundary conditions, L 1 , L 2 , and L 3 are defined, and V n+1 just shows ‖.‖ ∞ represents infinite norm. For each n � 0, 1, . . . , N − 1 and j � 1, . . . , M − 1, we have From (23)- (25), it is clear that x n j , y n j , z n j ≥ 0, so we will have Now, for j � 0, 1, . . . , M, let us consider the following two modes: Journal of Mathematics erefore, M |, it is seen from the boundary conditions that By combining two states (43) and (44), stability is achieved and this completes the proof.
e local truncation error for (18) is as follows: the approximate solution converges to the viscosity solution of the equation [43].

Numerical Results
To illustrate the effectiveness of the new NSFD, the problem in the independent variables (S, τ) is solved for various values of parameter transaction cost a on a number of uniform meshes by the numerical method presented in the previous sections. e numerical solutions are then transformed back in (S, t). We have programmed these methods in MATLAB.

Example 1.
e first test problem is which we consider European call option contracts for solving equation (6) with initial (8) and boundary conditions (9) and (10) by the following parameters: Option value with uniform meshes h � 4 and Δt � 0.1 is drawn up in Figure 1(a). It can be seen that the solution of the proposed method is positive and numerically stable. Also, the comparison of price options with different transaction costs at t � 0 is shown in Figure 1(b). It can be seen that with the increase in the transaction cost, the price of the option will also increase. is is also true in practice. e value of the parameters used in this section has been taken from [5].  Example 2. As our the second numerical experiment, consider put option contracts for solving equation (6) with initial (8) and boundary conditions (9) and (10) by the following parameters: Option value with uniform meshes h � 2 and Δt � 0.05 is drawn up in Figure 2(a). As can be seen from Figure 2, the function of the value of the option is an increasing function of the transaction cost parameter. Also, the price of the put option with different transaction costs at t � 0 is shown in Figure 2(b).
Example 3. Let us consider European butterfly spread option contracts for solving equation (6) with initial (8) and boundary conditions (9)   Option value with uniform meshes h � 4 and Δt � 0.1 is drawn up in Figure 3(a). As can be seen from Figure 3, the function of the value of the option is an increasing function and is the transaction cost parameter. Also, the price of the European put option with different transaction costs at t � 0 is shown in Figure 3(b).
Example 4. In our last example, we consider cash-ornothing option contracts for solving equation (6) with initial (8) and boundary conditions (9) and (10) by the following parameters: Option value with uniform meshes h � 2 and Δt � 0.05 is drawn up in Figure 4(a). As can be seen from Figure 4, the function of the value of the option is an increasing function of the transaction cost parameter. Also, the price of a cashor-nothing option with different transaction costs at t � 0 is shown in Figure 4(b).

Discussion and Conclusion
In this paper, we proposed a new nonstandard finite difference method for the Black-Scholes nonlinear equation under transaction costs, which, in addition to being free from spurious oscillation and preserving positivity property for the desired step size, is unconditionally stable. Also, we showed the convergence of the solution of the new method to the viscosity solution of the equation by proving stability, consistency, and monotonicity and observed that the price of a European option contract is an increasing function of the parameter transaction cost. Our results indicate that a properly implemented version of our scheme is useful for the numerical integration of the considered Black-Scholes equation. In our future research, we intend to solve time fractional Black-Scholes European option pricing equation using the two-scale fractal derivatives.

Data Availability
e data used to support this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.