In this paper, we use a numerical method for solving the nonlinear Black–Scholes partial differential 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 efficiency of the new method, numerical results have been performed by four models.

1. 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. The 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. The 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 . One of the models in the discussion of financial derivatives pricing is the Black–Scholes model. This 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 is the value of the option, is the current underlying asset price, t = t is the time, is the expiration date of the option, is the risk-free interest rate, and is the volatility (also called implicit volatility). and are fixed. Partial differential equation paid attention in the field of economics with the Black–Scholes model. Therefore, in order to obtain the price of these contracts, we need a simple and accurate method to solve this differential equation. The 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. Then, 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 [24]. We consider the nonlinear Barles and Soner model which is obtained as follows:where , in which is the risk-adjusted factor and is the number of options that are sold. The 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 , equation (1) is converted to the following equation:for and ,as in [5] we consider the initial condition, the first boundary condition, and the second boundary condition for and as follows:where , and denote the strike prices of the options, is a positive constant, and is the Heaviside function. Conditions and 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).

The 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 [811]. In [1214], 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 [1517] 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. They caused numerical instability in the solution that has been resolved by an iterative method. This 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 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 sixth-order 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 fourth-order finite difference approximation is contributed and equipped with the fourth-order Runge–Kutta scheme for the nonlinear Black–Scholes. They 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. They 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.

The 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. The numerical results are given in section 5.

2. Nonstandard Finite Difference Scheme

In this section, we give some preliminaries including nonstandard finite difference methods. The 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 [35, 13, 15, 16, 18, 21, 2638].

This 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. The forward finite difference method is one of the simplest discretization schemes. In this method, the derivative is replaced by , where is step size. However, in the Mickens schemes, this term is replaced by , where is an increasing continuous function of , and the function satisfies the following conditions:

Note that in taking to obtain the derivative, the use of any of this will lead to the usual result for the first derivative:where and are continuous functions of the step size . 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 is replaced by a nonnegative function such thatorExamples of functions 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 [35, 13, 15, 16, 18, 21, 2628], these are replaced by

One can say that there is no appropriate general method to choose the function or to choose which nonlinear terms are to be replaced, and some special techniques may be found in [35, 13, 15, 16, 18, 21, 2628].

3. 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 and discretize it in the following form. We introduce a grid of mesh points where , , and the spatial step-size is given by , and the time step size is . We denote the approximation of by . The new NSFD is of the following form:

The explicit form of (17) is as follows:

The finite difference approximation provides a difference equation as follows:with (8)–(10), and we define the initial and boundary conditions as follows:

The matrix form (18) for all and is given in the following form:in which is a tridiagonal matrix with positive components and is known as boundary values. The new scheme can be written in the following form:where

4. Analysis of the New Method

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

For and , we define the function as follows:where

Lemma 1 (monotonicity). The new NSFD method is monotone.

Proof. As in [5] which is enough to show that for each and ,Since , , and , the first three sentences of the right of the relationship (26), respectively, are nonincreasing in , increasing in , and decreasing in .
Suppose is a column vector . We haveAssume that we take into account the nonlinear sentence is to the right of the relation (26). From (3), we havewherefor each and ; independent of , by derivative phrase relative to the variable and using (5) and (24) and , we haveSo, is an increasing function of . Using the above symbols and by monotone combination of and properties of the linear sentences of the relationship (26) and (30), we haveThis is the same as (28). Similarly from monotonicity 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, 3942]).

Proposition 1 (positivity). Suppose , , and are the real numbers that are nonnegative, in this case, method (18) offers a nonnegative approximation of for the solution of equation (2).

Proof. Since the parameteralsothenTherefore, this completes the proof.

Lemma 2 (stability). The new NSFD method is unconditionally stable.

Proof. As in [5], for each , we know , that is the solution of the NSFD method, and then for the initial and boundary conditions, , , and are defined, and just shows represents infinite norm. For each and , we haveFrom (23)–(25), it is clear that , , , so we will haveNow, for , let us consider the following two modes:(i)Assume , then since , Therefore,(ii)Suppose that , then if or , it is seen from the boundary conditions thatBy combining two states (43) and (44), stability is achieved and this completes the proof.

Lemma 3 (consistency). The proposed method (18) is consistent, and its local truncation error is .

Proof. The local truncation error for (18) is as follows:We have Taylor expansionsreplacing these expressions in the term as follows:However, is the solution of the Black–Scholes equation, so the proposed method is consistent:Therefore, the main part of the local truncation error is as follows:so . Therefore, this completes the proof.

Theorem 1. The solution of method (18) convergences to the viscosity solution of equation (2), when

Proof. If the discretized method for nonlinear partial differential equations be stable, consistent, and monotone, then the approximate solution converges to the viscosity solution of the equation [43].

5. Numerical Results

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

Example 1. The 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 and 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 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. This is also true in practice. The 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 and 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 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) and (10) by the following parameters:

Option value with uniform meshes and 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 is shown in Figure 3(b).

Example 4. In our last example, we consider cash-or-nothing option contracts for solving equation (6) with initial (8) and boundary conditions (9) and (10) by the following parameters:

Option value with uniform meshes and 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 cash-or-nothing option with different transaction costs at is shown in Figure 4(b).

6. 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

The data used to support this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.