A Nonstandard Finite Difference Method for a Generalized Black–Scholes Equation

: An implicit ﬁnite difference scheme for the numerical solution of a generalized Black– Scholes equation is presented. The method is based on the nonstandard ﬁnite difference technique. The positivity property is discussed and it is shown that the proposed method is consistent, stable and also the order of the scheme respect to the space variable is two. As the Black–Scholes model relies on symmetry of distribution and ignores the skewness of the distribution of the asset, the proposed method will be more appropriate for solving such symmetric models. In order to illustrate the efﬁciency of the new method, we applied it on some test examples. The obtained results conﬁrm the theoretical behavior regarding the order of convergence. Furthermore, the numerical results are in good agreement with the exact solution and are more accurate than other existing results in the literature.


Introduction
An option is the right, but not the obligation, to buy or sell an underlying asset at a specific price on or before a certain date at a fixed strike price. European call (put) option and American call (put) option are two types of known options. European option can be exercised only on expiration date T. For American option, exercise is permitted at any time t ≤ T. It was shown in [1] that option pricing can be modeled using a partial differential equation of second order. The resulting equation is called the Black-Scholes equation. If the coefficients are constant or independent of the spatial variable, it can be transformed into the standard heat equation and can be solved exactly. However, in the financial market, the coefficients can depend on the time and the asset price. In such situation the analytical solution of the generalized Black-Scholes model is not available. Therefore, numerical simulations are essential to obtain information about the behavior of the solutions. It is important to develop numerical methods that reproduce the qualitative properties of the solution.
Equation (1) is a partial differential equation used for pricing options. The function C(x, t) is the European call option price with asset price x and time t. T is the maturity, K is the price of strike, r(t) is the risk-free interest rate, and σ(x, t) is the volatility function of the underlying asset and time. When the coefficients are constant functions, we get the classical Black-Scholes model. Following [2] we consider the above model on a truncated domain R = (0, S max ) × (0, T), where S max > K is a suitably chosen positive number so that The result about existence and uniqueness of a solution of (4)- (7) when σ(x, t) and r(t) are constant functions can be found in [3,4]. It is proved in [5] that if C and C * are solutions of (1)-(3) and (4)-(7), respectively, then The terminal condition (5) is not a smooth one. Therefore, to guarantee the convergence of the numerical solution, we substitute the terminal condition by a four order smooth function (see in [2]) and obtain the following model: where the function ψ ε is defined as follows: for a sufficiently small ε > 0, ψ ε (y) satisfying By using these conditions one can easily find that The proof of the existence and uniqueness of a solution of (9)- (12) can be found in [3]. A proof of the following result can also be found there: for (x, t) ∈R.
From (8) and (16) we have that we can take the solution of the modified models (9)-(12) closer to that of the original models (1)-(3) by choosing a sufficiently large value for S max and a sufficiently small value for ε.
Traditionally, important requirements for constructing a numerical method are the investigation of the consistency of the discrete scheme with the original differential equation and linear stability analysis with smooth solutions [6][7][8][9][10]. These requirements are important, because they are necessary to ensure the convergence of the discrete solution to the exact one, but for the explicit numerical methods based on standard finite difference approaches. However, most of the time, the qualitative properties of the true solution are not reproduced by the numerical solution. Therefore, this might result in a calamitous erroneous outcome. One way to overcome this disadvantage is to use nonstandard finite difference methods (NSFDs). The origin of the NSFD methods began with Mickens in 1989 [11]. A summary of the known results up to 1994 are given in [12]. The simplicity in construction and the possiblity of their extensions allowed researchers to apply NSFD schemes to solve several complex differential equation models that arise at different disciplines. A scheme is called nonstandard if at least one of the following conditions is satisfied: The usual denominator h or h 2 of the discrete first or second order derivatives is replaced by a non-negative function Φ such that Some examples of Φ(h) satisfying these conditions are 2.
The nonlinear terms that appear in the differential equation are approximated in a non-local way, for instance, the nonlinear terms u 2 and u 3 can be replaced by (see [13]): In addition to the usual properties of consistency, stability and convergence, the NSDF methods give rise to numerical solutions that also maintain essential qualitative properties of the solutions [14][15][16][17][18][19][20][21].
In this paper, we introduce a new scheme for the numerical solution of the generalized Black-Scholes equation. We propose a non-standard finite difference method, where we use a non-local approximation of the reaction term of the generalized Black-Scholes equation, combined with an implicit time step technique. The method is conditionally stable, positivity preserving and of second-order of convergence with respect to the space variable. Our proposed method has some advantages: A uniform mesh is used; thus, the analysis is simple and it has better accuracy and convergence rate.
The rest of the paper is organized as follows. In Section 2, we review the works done on the generalized Black-Scholes equation and the real options. In Section 3, we propose a NSFD scheme for such a model and we present the analysis of the method concerning the positivity preserving property, Stability, and truncation error. Furthermore, the numerical results obtained from the new method are presented and compared with the results in [2]. The last section describes the features of the new method, the results obtained, and ends with the presentation of ideas for future works.

Literature Review
In this section, we review the works done on the generalized Black-Scholes equation.
In addition, at the end of the section, we will mention some of the work done related to real options. There are several methods for solving the generalized Black-Scholes equation in the literature. Cen and Le [2] presented an accurate finite difference method, based on a central difference discretization in the space variable using a piecewise uniform mesh and an implicit method in the temporal variable. Cen et al. [22] proposed an exponential scheme in the temporal variable combined with a central difference approach on a spatial piecewise uniform mesh. Kadalbajoo et al. [23] suggested a cubic B-spline collocation method for the generalized model with second order accuracy in both space and time. A fitted finite volume method was used in [24] for the generalized model, showing a second-order convergence. A Cubic spline method was proposed in [25] for a generalized Black-Scholes equation. In that work, the implicit Euler method was used in the temporal discretization, together with a cubic polynomial spline approximation for the spatial discretization. In [26], Mohammadi achieved third order of accuracy in space considering a quintic B-spline collocation method and first-and second-order of accuracy in time using the backward Euler method and the Crank-Nicolson method, respectively. In [27], a highorder numerical method was obtained. In this method a two-step backward differentiation formula is used in the temporal discretization while a high-order difference approximation was used in the spatial discretization. This gives a second order of accuracy in time and third order of accuracy in space.
In [28], Clevenhaus et al. added a penalty term to the American option problem, by doing this the problem was reduced to a PDE on a fixed domain. Their presented penalty term is bounded by an initial guess of the free boundary value and is thus independent of the solution at the last time step. The novelty and advantage of their approach consist in introducing a time-dependent penalty function which led to the construction of an efficient adaptive numerical scheme. In [29] a new method for pricing American options by Boire et al. is presented. The numerical results show that their method successfully reduces the bias plaguing the standard importance sampling method across a wide range of moneyness and maturities, with negligible change to estimator variance.
In [30], considering several uncertainty factors: a renewable energy production amount, an R and D investment amount, a unit price, a risk-free interest rate, and the R and D investment in the six types of renewable energy sources: biomass energy, waste energy, photovoltaic energy, solar thermal energy, marine energy, and wind power energy are evaluated. They utilize the system dynamics approach using the Black-Scholes model and conducts sensitivity analyses of the uncertainty factors. Larissa et al. [31] investigate the contribution of adjusted net savings to sustainable economic growth for 10 Central and Eastern European and Baltic nations using data corresponding to the period 2005-2016. Their results indicated that adjusted net savings impacted on the GDP across the 10 countries analyzed. A bibliometric analysis to examine the features of Carbon capture and storage (CCS) literature including the research focus and trends as well the real options uncertainty and models, types of options, and valuation techniques in [32] are presented. Their results provide energy and environmental policy-makers and CCS project planners with valuable insights into various aspects of CCS policy and project design. For real options in the field of agriculture, we can refer to the article [33] in which the authors investigated a hydroponic farm considering uncertainty, sustainability, and the system's utility in the dominant desert geography.

Construction of the New Scheme
In this subsection, by using the strategy of nonstandard discretization methods for generalized Black-Scholes models (9)-(12), we consider the forward difference for ∂u ∂t , the centered differences for ∂u ∂x , ∂ 2 u ∂x 2 , and ξ(U to approximate the u in the reaction term. Thus, we propose a new discretization scheme for (9) given by where U j i , σ j i and r j are approximations of u(x, t), σ(x, t) and r(t), respectively at the grid The above discretization may be written for each j in a matrix form as where U j = (U T are vectors containing approximations to u at time levels t j and t j+1 , respectively. Furthermore, d j and d j+1 are column vectors that contain the known boundary-values and zeros, and A j and B j are tridiagonal matrices of the form

Characteristics of the New Scheme
This subsection is concerned with the positivity preservation, stability, consistency and convergence properties of the proposed method. Firstly, we present some definitions and results that will be needed in the proof of Theorem 3.

Definition 2.
A matrix A = (a ij ) is called an L-matrix if a ii > 0, ∀i ∈ N and a ij ≤ 0, i = j.

Theorem 2 ([34]
). Let A be an L-matrix which is strongly row or column diagonally dominant, i.e., Ae > 0 or e T A > 0 T . Then, A is a nonsingular M-matrix.
then, the method given in (21) is positivity preserving.
Proof. From (21), we have In matrix form, we have where As the final condition is positive, if (A j ) −1 > 0, B j ≥ 0, and(d j+1 − d j ) > 0, then using backward recursion we get that U j is positive for every j. The condition (A j ) To be A j an M-matrix it should happen that From (28), we should have and from (29), we should have On the other hand, for the non-negativity of B j , it is sufficient to have Then, from (30)-(32), if we choose ξ so that The first boundary is zero, therefore it is to show that f  Proof. Under condition (24), A j = [A j ir ] is similar to a symmetric tridiagonal matrix, so the eigenvalues of A j , λ i (A j ), i = 1, . . . , N, are real [34]. Furthermore, A j is row diagonally dominant with [34], Proposition 1.3). Therefore, we have and by combining with B j ∞ = 1 ∆t + 2ξr j − r j , we obtain where ρ((A j ) −1 B j ) is the spectral radius of (A j ) −1 B j . Therefore, the method is stable (see [35], p. 54) and then, via the Lax-theorem, convergent. The local truncation error is given by By Taylor's expansion, we have and by substituting these expressions in (38), we obtain As u is the solution of Equation (9), then Therefore, the principle term of the local truncation error is which completes the proof.

Numerical Results
In this subsection, we give some numerical examples to confirm the theoretical results. In the last two examples, the results of our method are compared with the results obtained in [2]. The values of volatility and interest rate parameters have been taken from the works in [2,22,23,25,27]. These values are given in Table 1 with references. The other parameters for calculating the option price of the (9)-(12) model to all examples are equal T = 1, K = 25, S max = 100 and ε = 10 −4 . As in [2], we have used the approximate solution obtained for N = 2 11 and M = 2 10 as a reference solution on the mesh points. We compute the maximum error as follows: whereŪ(x, t) is a linear interpolation of the approximate solution U 2 11 ,2 10 to get solutions at other points. Furthermore, we have used the following formula to compute the convergence rate:R The error estimates along with the convergence rates, obtained with the proposed NSFD method, are listed in Tables 2-8. From Tables 2-8, we see that the convergence order in space is almost 2, which confirms the estimate given in Theorem 4.  It can be seen from Tables 6 and 7

Conclusions
In this paper, the generalized Black-Scholes equation for European option pricing is studied. To numerically solve this equation, a non-standard finite difference method has been used. In the proposed method, a non-local expression to approximate the reaction sentence of the generalized Black-Scholes and an implicit time step technique is used.
The theory results show that the proposed scheme is positivity preserving, stable and the order of the scheme with respect to the space variable is two. To confirm the theory results and good performance of the new method the option price for different values of volatility and interest rates is calculated. The numerical results confirm the theoretical behavior regarding the order of convergence for Theorem 4 (see Tables 2-8). For the last two examples, the option price values are plotted for all time levels, which shows that the new method is positive and non-oscillation (see Figures 1 and 2). Furthermore, the obtained results are more accurate than other existing results in literature.
The results obtained indicate that non-standard difference schemes may be promising for solving problems which may have an impact on the stock price such as nonlinear Black-Scholes equations. In our subsequent investigation, we intend to address these problems using non-standard methods.