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

Authors: Mohammad Mehdizadeh Khalsaraei; Ali Shokri; Zahra Mohammadnia; Hamid Mohammad Sedighi
Journal: Journal of Mathematics
 (2021)

---

## Abstract

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.


---

## Body

## 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 datet≤T. 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:(1)Ut+12σ02S2USS+rSUS−rU=0,where U=US,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. 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:(2)Ut+12σ2t,S,US,USSS2USS+rSUS−rU=0,in which σ is the modified volatility. Various nonlinear volatility models have been presented in [2–4]. We consider the nonlinear Barles and Soner model which is obtained as follows:(3)σ2=σ021+ΨerT−ta2S2USS,where a=μγN, in which γ is the risk-adjusted factor and N is the number of options that are sold. The function Ψ is the solution of the following nonlinear initial value problem [5]:(4)Ψ′w=Ψw+12wΨw−w,w≠0,Ψ0=0.For functionΨ, an exact implicit solution is presented as follows [6]:(5)w=−sinh−1ΨΨ+1+Ψ,forw>0,Ψw>0,−w=sin−1−ΨΨ+1−−Ψ,forw<0,−1<Ψw<0.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:(6)Vτ=12σ2VSSS2VSS+rSVS−rV,for S>0 and 0<τ≤T,(7)σ2VSS=σ021+Ψerτa2S2VSS,as in [5] we consider the initial condition, the first boundary condition, and the second boundary condition for S∈0,Smax and τ∈0,T as follows:(8)L1S=maxS−K,0,for the Call option,maxK−S,0,for the Put option,maxS−K1,0−2maxS−K2,0+maxS−K3,0,for Butterfly Spread,B×HS−K,for Cash or Nothing,(9)L2τ=0,for the Call option,Ke−rτ,for the Put option,0,for Butterfly Spread,0,for Cash or Nothing,(10)L3τ=Smax−Ke−rτ,for the Call option,0,for the Put option,0,for Butterfly Spread,Be−rτ,for Cash or Nothing,where K,K1,K2, and K3 denote the strike prices of the options, B is a positive constant, and H is the Heaviside function. Conditions L10=L20 and L1Smax=L30 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 [8–11]. In [12–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–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. 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 x=lnS/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 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 section2, 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 [3–5, 13, 15, 16, 18, 21, 26–38].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 derivativeVx is replaced by Vx+h−Vx/h, where h is step size. However, in the Mickens schemes, this term is replaced by Vx+h−Vx/ϕh, where ϕh is an increasing continuous function of h, and the function ϕh satisfies the following conditions:(11)ϕh=h+Oh2,0<ϕh<1,h⟶0.Note that in takinglimh⟶0 to obtain the derivative, the use of any of this ϕh will lead to the usual result for the first derivative:(12)dVdx=limh⟶0Vx+ϕ1h−Vxϕ2h=limh⟶0Vx+h−Vxh,where ϕ1h and ϕ2h 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 denominatorh is replaced by a nonnegative function ϕh such that(13)ϕh=h+Oh2,0<ϕh<1,h⟶0,or(14)ϕh2=h2+Oh4,0<ϕh2<1,h⟶0.Examples of functionsϕh that satisfy these conditions are as follows [25]:(15)h,1−e−h,sinh,or4sin2h2.(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(16)V≈3Vj+1n−2Vjn,V≈2Vj+1n+Vjn3,V2≈aVjn2+bVjnVj+1n,a+b=1,a,b∈ℝ.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].


## 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Ω=0,Smax×0,T and discretize it in the following form. We introduce a grid of mesh points S,t=Sj,tn where Sj=jhandtn=nΔt, j=0,1,…,M;n=0,1,…,X, and the spatial step-size is given by h=Smax/M, and the time step size is Δt=T/X. We denote the approximation of VSj,tn by Vjn. The new NSFD is of the following form:(17)−Vjn+1−VjnΔτ+12σ2VSSnjSj2Vj−1n−2Vjn+1+Vj+1nh2+rSjVj+1n−Vjn+1h−rVjn+1=0.The explicit form of (17) is as follows:(18)Vjn+1=1/2σ2VSSnjSj2/h2Vj−1n+1/ΔτVjn+rSj/h+1/2σ2VSSnjSj2/h2Vj+1n1/Δτ+rSj/h+σ2VSSnjSj2/h2+r.The finite difference approximation provides a difference equation as follows:(19)Vn+1=AnVnVn+Bn,with (8)–(10), and we define the initial and boundary conditions as follows:(20)Vj0=L1Sj,V0n=L2τn,VMn=L3τn.The matrix form (18) for all j=1,…,M−1 and n=0,1,…,N−1 is given in the following form:(21)AnVn=y1nz1n0⋯000x2ny2nz2n⋯0000x3ny3n⋯000⋮⋮⋮⋱⋮⋮⋮000⋯yM−3nzM−3n0000⋯xM−2nyM−2nzM−2n000⋯0xM−1nyM−1n,Bn=x1nV0n00⋮00zM−1nVMn,in which AnVn is a tridiagonal matrix with positive components and Bn is known as boundary values. The new scheme can be written in the following form:(22)Vjn+1=xjnVj−1n+yjnVjn+zjnVj+1n,where(23)xjnVn=σ2VSSnjSj2/2h21/Δτ+rSj/h+σ2VSSnjSj2/h2+r,(24)yjnVn=1/Δτ1/Δτ+rSj/h+σ2VSSnjSj2/h2+r,(25)zjnVn=rSj/h+1/2σ2VSSnjSj2/h21/Δτ+rSj/h+σj2VSSnjSj2/h2+r.


## 4. Analysis of the New Method

In this section, we investigate the monotonicity, positivity preserving, stability, consistency, and convergency properties of the new method.Forj=1,…,M−1 and n=0,1,…,N−1, we define the function φjn+1 as follows:(26)φjn+1Vjn+1,Vj+1n,Vj−1n,Vjn≔pjVj+1n+qjVjn+1−1ΔτnVjn−12Sj2σ2VSSnjVSSnj=0,where(27)pj=−rSjhj,qj=1Δτn+rSjhj+r.Lemma 1 (monotonicity).
The new NSFD method is monotone.Proof.
As in [5] which is enough to show that for each ξ>0 and j=1,…,M−1,(28)φjn+1Vjn+1,Vj+1n+ξ,Vj−1n+ξ,Vjn+ξ≤φjn+1Vjn+1,Vj+1n,Vj−1n,Vjn,(29)φjn+1Vjn+1+ξ,Vj+1n,Vj−1n,Vjn≥φjn+1Vjn+1,Vj+1n,Vj−1n,Vjn.
Sincepj≤0, 1/Δτn>0, and qj>0, the first three sentences of the right of the relationship (26), respectively, are nonincreasing in Vj+1n+1, increasing in Vjn+1, and decreasing in Vjn.
SupposeEk=0,0,…,1︸kth,0,…,0T is a column vector M−1×1. We have(30)VSSVn+ξEj−1+ξEj+1j=1h2Vj−1n+ξ−2h2Vjn+1+1h2Vj+1n+ξ=VSSnj+2h2ξ,(31)VSSVn+ξEjj=1h2Vj−1n−2h2Vjn+1+ξ+1h2Vj+1n=VSSnj−2h2ξ.
Assume that we take into account the nonlinear sentence is to the right of the relation (26). From (3), we have(32)12Sj2σ2vjnvjn=12σ02Sj21+ΨKjn+1vjnvjn,where(33)vjn=VSSj,Kjn+1=erτn+1a2Sj2>0,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(34)ddv1+ΨKvv=Ψ′KvKv+1+ΨKv=1+Ψw2wΨw−ww+1+Ψw=1+Ψw2wΨw2wΨw−w=2wΨwΨ′w≥0.
So,1+ΨKvv is an increasing function of v. Using the above symbols and by monotone combination of 1+ΨKvv and properties of the linear sentences of the relationship (26) and (30), we have(35)φjn+1Vjn+1,Vj+1n+ξ,Vj−1n+ξ,Vjn+ξ=pjVj+1n+ξ+qjVjn+1−1ΔτnVjn+ξ−12σ02Sj21+ΨKjn+1vjn+2h2ξvjn+2h2ξ≤φjn+1Vjn+1,Vj+1n,Vj−1n,Vjn.
This is the same as (28). Similarly from monotonicity 1+ΨKvv 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–42]).Proposition 1 (positivity).
SupposeVj−1n, Vjn, and Vj+1n are the real numbers that are nonnegative, in this case, method (18) offers a nonnegative approximation of Vjn+1 for the solution of equation (2).Proof.
Since the parameter(36)S,t,T,r,h,Δτ,σ0≥0,also(37)σ2≥0,then(38)xjn,yjn,zjn≥0.
Therefore, this completes the proof.Lemma 2 (stability).
The new NSFD method is unconditionally stable.Proof.
As in [5], for each n=0,1,…,M−1, we know Vn+1=V0n,VnT,VMn, that Vn+1 is the solution of the NSFD method, and then for the initial and boundary conditions, L1, L2, and L3 are defined, and Vn+1 just shows(39)Vn+1∞≤maxL1∞,L2∞,L3∞.
.∞ represents infinite norm. For each n=0,1,…,N−1 and j=1,…,M−1, we have(40)Vjn+1=xjnVj−1n+yjnVjn+zjnVj+1n.
From (23)–(25), it is clear that xjn, yjn, zjn≥0, so we will have(41)Vjn+1≤xjnVj−1n+yjnVjn+zjnVj+1n≤xjnVn+yjnVn+zjnVn=xjn+yjn+zjnVn.
Now, forj=0,1,…,M, let us consider the following two modes:(i)
Assumej=1,…,M−1, then Vn+1∞=Vkn+1k∈1,2,…,M−1andsince r≥0, (42)xkn+ykn+zkn=1/2σkSk/h2+1/Δτ+rSk/h+1/2σkSk/h21/Δτ+rSk/h+σkSk/h2+r≤1.Therefore,(43)Vn+1∞≤xkn+ykn+zknVn∞≤Vn∞≤Vn−1∞≤⋯≤V0∞≤L1∞.(ii)
Suppose thatj=0,M, then if Vn+1∞=V0n+1 or Vn+1∞=VMn+1, it is seen from the boundary conditions that(44)Vn+1∞≤maxV0n+1,VMn+1≤maxL2∞,L3∞.
By 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 OΔτ,h2.Proof.
The local truncation error for (18) is as follows:(45)Tjn=−VSj,τn+1−VSj,τnΔτ+rSjVSj+1,τn−VSj,τn+1h+12σj2Sj2VSj−1,τn−2VSj,τn+1+VSj+1,τnh2−rVSj,τn+1.
We have Taylor expansions(46)VSj,τn+1=VSj,τn+Δτ∂V∂τjn+Δτ22∂2V∂τ2jn+Δτ36∂3V∂τ3jn+⋯,VSj+1,τn=VSj,τn+h∂V∂Sjn+h22∂2V∂S2jn+h36∂3V∂S3jn+⋯,VSj−1,τn=VSj,τn−h∂V∂Sjn+h22∂2V∂S2jn−h36∂3V∂S3jn+⋯,replacing these expressions in the term Tjn as follows:(47)Tjn=−∂V∂τ+rS∂V∂S+12σ2S2∂2V∂S2−rVjn−rjΔτ+σ2j2Δτ+rΔτ∂V∂τjn+12rjh2∂2V∂S2jn−12Δτ+12rjΔτ2+12σ2j2Δτ2+12rΔτ2∂2V∂τ2jn+⋯.
However,V is the solution of the Black–Scholes equation, so the proposed method is consistent:(48)−∂V∂τ+rS∂V∂S+12σ2S2∂2V∂S2−rVjn=0.
Therefore, the main part of the local truncation error is as follows:(49)−rjΔτ+σ2j2Δτ+rΔτ∂V∂τjn+12rjh2∂2V∂S2jn,so OΔτ,h2. Therefore, this completes the proof.Theorem 1.
The solution of method (18) convergences to the viscosity solution of equation (2), when(50)h=max0≤j≤M−1hj⟶0,Δτ=max0≤n≤N−1Δτn⟶0.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 variablesS,τ 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. The numerical solutions are then transformed back in S,t. 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:(51)r=0.1,a=0.05,σ=0.2,K=40,T=1,Smax=80.Option value with uniform meshesh=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. This is also true in practice. The value of the parameters used in this section has been taken from [5].Figure 1
(a) European call option price charts of the explicit nonstandard method fora=0.05 and (b) European option price charts of the explicit nonstandard method for different transaction costs in t=0.
(a)(b)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:(52)K1=30,K2=40,K3=50,r=0.1,σ0=0.2,T=1,Smax=80.Option value with uniform meshesh=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).Figure 2
(a) European put option price charts of the explicit nonstandard method fora=0.05 and (b) European option price charts of the explicit nonstandard method for different transaction costs in t=0.
(a)(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:(53)r=0.1,a=0.05,σ=0.2,K=40,T=1,Smax=80.Option value with uniform meshesh=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).Figure 3
(a) Butterfly spread option price charts of the explicit nonstandard method fora=0.05 and (b) butterfly spread option price charts of the explicit nonstandard method for different transaction costs in t=0.
(a)(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:(54)r=0.1,σ0=0.2,K=40,T=1,B=1,Smax=80.Option value with uniform meshesh=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 cash-or-nothing option with different transaction costs at t=0 is shown in Figure 4(b).Figure 4
(a) European cash-or-nothing option price charts of the explicit nonstandard method fora=0.05 and (b) cash-or-nothing charts of the explicit nonstandard method for different transaction costs in t=0.
(a)(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.



---

*Source: 6679484-2021-05-11.xml*