An Unconditionally Stable, Positivity-Preserving Splitting Scheme for Nonlinear Black-Scholes Equation with Transaction Costs

This paper deals with the numerical analysis of nonlinear Black-Scholes equation with transaction costs. An unconditionally stable and monotone splitting method, ensuring positive numerical solution and avoiding unstable oscillations, is proposed. This numerical method is based on the LOD-Backward Euler method which allows us to solve the discrete equation explicitly. The numerical results for vanilla call option and for European butterfly spread are provided. It turns out that the proposed scheme is efficient and reliable.


Introduction
One of the modern financial theory's biggest successes in terms of both approach and applicability has been the Black-Scholes option pricing model developed by Black and Scholes in 1973 [1] and previously by Merton [2]: where > 0 and 0 are given real constants that represent the interest rate and the volatility, respectively, is the price of the underlying asset, and is the maturity date. The celebrated Black-Scholes model is based on several restrictive assumptions such as liquid, frictionless, and complete markets. In recent years, nonlinear Black-Scholes models have been used to build transaction costs, market liquidity, or volatility uncertainty into the celebrated Black-Scholes concept.
In this paper, we are interested in the option pricing model with transaction costs proposed by Barles and Soner [3] that are motivated by Hodges and Neuberger [4]. In practice, transaction costs arise when trading securities.
Recent studies of their influence reveal that they result in a nonnegligible increase in the option price, although they are generally small for institutional investors. To show this increase, in Barles and Soner's model, the constant volatility 0 in the linear model is replaced by the nonlinear volatility that reads where = √ with the proportional transaction cost being , the risk aversion factor being , and the number of options to be sold being , and Ψ( ) is the solution of the following ordinary differential equation: Barles and Soner's option pricing model now reads The Scientific World Journal with the terminal condition ( , ) = ( ) , ∈ Ω := (0, +∞) .
The payoff function ( ) is assumed to be a continuous piecewise linear function. Many results have been reported for the numerical solutions of linear Black-Scholes equations (see, e.g., [5][6][7][8][9]). On the other hand, because of the nonlinear nature of this model, numerical methods are mandatory to price derivatives and portfolios. The strong nonlinearity of problem (4) makes it difficult to obtain the reliable numerical solutions. Implicit numerical schemes have been used for numerically solving nonlinear option pricing PDEs [10], and an iterative approach is required to solve the nonlinear algebraic equation resulting from the discretization, which results in more computational cost. Consequently, some high-order compact semi-implicit difference schemes are proposed for numerically solving the nonlinear option pricing model [11,12]. Some researchers (see, e.g., [3,[13][14][15]) have also constructed explicit finite difference schemes for (4)-(5) and investigated their consistency and stability. However, these explicit schemes have the disadvantage that strictly restrictive conditions on the discretization parameters are needed to guarantee stability and positivity. To relax the restrictive conditions, in [16], Zhou et al. proposed an unconditionally stable explicit finite difference scheme based on a nonstandard approximation of the second partial derivative. However, this scheme is conditionally consistent, and the truncation error depends on the ratio of the time stepsize and the square of the space stepsize.
In this paper, we will consider a splitting method with inhomogeneous boundary conditions. This method proposed here is unconditionally stable, monotone, and positivity preserving. It is also consistent and essentially a "limit" version of the LOD-Backward Euler method (see Chapter IV on splitting methods in [17]) and therefore allows us to solve the discrete equation explicitly.
The paper is organized as follows: we begin by transforming the original equations into nonlinear heat equations and considering the spatial semidiscretization. The splitting scheme will be discussed in Section 3 after linearization of semidiscrete system. The stability, the monotonicity, the positivity-preserving property, and the convergence of this scheme are analysed in Section 4. To illustrate our method, we present some numerical experiments in Section 5. Finally, we give a summary.

Preliminaries.
We first consider the properties of the function Ψ appearing in (2), which will play an important role in the following numerical analysis.
we transform the problem (4)-(5) into the following initial value problem The Scientific World Journal 3 The boundary conditions take the form of (0, ) = 0, ( , ) = (+∞, ) , where the second boundary condition in (16) is derived by asymptotic considerations in ( , ) coordinates for extreme value → ∞ and subsequent transformation. For example, the boundary conditions for a vanilla call are given by and the boundary conditions for a butterfly spread are given by Note that the second equality in (19) derives from the second equality in (16). Taking these into consideration, the boundary conditions for a vanilla call are given by where is the strike price, and the boundary conditions for a butterfly spread are given by We introduce the spatial grid Ω ℎ with step ℎ by the nodes = ℎ, = 0, 1, . . . , , so that ℎ = . After performing the second-order central finite difference approximation of the partial derivative ( 2 / 2 )( , ) as we obtained the corresponding ODEs system for the semidis- with ( ) = 1 ℎ 2 tridiag ( ( ) , −2 ( ) , ( )) , where ∈ R −1 is a vector, generated by the boundary conditions, Obviously, the semidiscrete difference scheme (23) is consistent. For the stability, we have the following theorem.
Theorem 4. For system (23), the following maximum norm contractivity holds: wherẽis the solution of the perturbation problem with the initial datã(0).

Splitting Time-Stepping Method
In this section, we will focus on the time integration methods of system (23). The implicit numerical schemes are generally viewed as stable methods, but [10] showed that the stability of the BDF schemes with orders 1 and 2 and the Crank-Nicolson scheme is still restricted by a condition. Additionally, for the implicit schemes, the nonlinear iteration will be required and will produce the additional computational cost in each time step. Fully explicit finite difference schemes for the PDEs (14)- (15) are also constructed in [3,[13][14][15]. However, these schemes are stable only for the severe restriction on the time stepsize Δ . For example, in [14], the condition reads To relax the condition, in this paper, we will propose an unconditionally stable method, which allows us to solve the discrete equation explicitly, based on the LOD-Backward Euler method (see Chapter IV on splitting methods in [17]).

Linearization
System. Let us set = Δ , = 1, 2, . . . , , for the temporal stepsize Δ = / . To linearize the nonlinear system (23), we allow the nonlinearities in (23) to lag one step behind and obtain the following linear system: For this linearized system, we have and therefore
Then, if the Backward Euler method is used to solve every subproblem, we get The Scientific World Journal 5 where is an approximation of ( ). It is easily verified that the matrix −Δ , is -matrix and therefore is nonsingular. By induction, we have which implies that for the splitting (41) For other choices of , , we can similarly obtain corresponding splitting schemes. By brief calculation, one can obtain an explicit expression of matrix [ − Δ , ] −1 : Then, we get the following explicit formula for computing , where is an approximation of ( , ): For the finite-dimensional discrete system (33), the splitting in the numerical scheme (45) is such that all computations become effectively "truly" one-dimensional. The numerical scheme (45) can be viewed as a "limit" version of LOD-Backward Euler method (see Chapter IV on splitting methods in [17]).

Properties of the Numerical Scheme
In this section, we investigate some properties of the numerical scheme proposed here.

Stability.
In [17], the stability of LOD-Backward Euler method is provided. The following theorem shows the stability of the numerical scheme (45).

Theorem 5.
The numerical scheme (45) is unconditionally stable in both the spectral norm ‖ ⋅ ‖ and the maximum norm ‖ ⋅ ‖ ∞ ; that is, one has the following stability inequalities: It should be pointed out that (49) can be regarded as a discrete version of maximum norm contractivity (26).
Proof. To obtain the spectral norm stability inequality (48), we apply the Gerschgorin theorem to the matrix , and get that the nonzero matrix eigenvalues line in the disc and therefore they are nonpositive. Then, any of the eigenvalues of the matrix [ − Δ , ] −1 , corresponding to the eigenvalues of , , satisfies Thus, it is easy to obtain ∏ −1 =1 | | ≤ 1 and The Scientific World Journal the matrix. The stability inequality (48) follows from this estimate.
We note that which directly leads to the maximum norm stability inequality (49) (see [17]).

Positivity.
A nice property of the numerical scheme for the pricing equation is positivity preserving, since the value of option is nonnegative.

Monotonicity.
Let us now consider the monotonicity of the numerical scheme (45). To do this, we note that scheme (45) can be written as For scheme (53), we have the following definition of monotonicity.
Definition 7 (see [19]; see also [20,21]). Scheme (53) is said to be monotone if and only if is nondecreasing in each argument.
The following theorem shows the monotonicity of the numerical scheme (45).
which also implies maximum norm estimate; that is,

Local Error Analysis and
Consistency. Now, we consider the error. Let = ( ) −̃denote the local discretization error, that is, the error introduced in one single step of the method. To bound the local discretization errors , we need the following truncation error estimation.
We can now immediately conclude that the splitting scheme (45) is consistent.

Viscosity Solutions and Convergence.
The so-called viscosity solutions are the meaningful solutions in financial applications. This has been shown in [23,24]. Putting all results together, we can formulate the following convergence result.
Proof. The convergence of the fully-discrete scheme (45) is a direct result of the consistency, the stability, and the monotonicity of this scheme [25].

Numerical Experiments
In order to illustrate the stability and convergence properties of our proposed scheme, in this section, we present several numerical experiments in which the vanilla call option and the European butterfly spread are considered. To obtain the numerical solution of nonlinear Black-Scholes equations (14)-(16), we should first solve the ODE (3). In this paper, we use the symmetric midpoint scheme to solve numerically the ODE (3). The utility function Ψ( ) can be done by linear interpolation.

Example 1. Let us consider a vanilla European call option with
This example comes from [14]. The parameters are the following: the strike price = 100, the volatility 0 = 0.2, the interest rate = 0.02, the maturity date = 1 year, and the artificial boundary location = 200. We first consider a small time stepsize Δ = 1/ = 0.0002 and a spatial stepsize ℎ = 200/ = 4 as done in [14]. The numerical results ( , ) computed by the numerical scheme (45) are plotted in Figure 1 in ( , ) coordinate. The case of = 0 is linear model, and the cases of = 0.015 and = 0.1 are nonlinear model. The numerical results are the same as those showed in [14].
To further confirm that this scheme is unconditionally stable, monotone, and positivity preserving, we also calculate the numerical solutions, which are presented in Figures 2 and  3, by scheme (45) with larger time stepsizes Δ = 0.002 and Δ = 0.02, respectively. From Figures 2 and 3, we observe that there are no stability issues for scheme (45). We also note that the scheme proposed in [14] produces the wrong numerical solution when the stepsize Δ = 0.00027027 such that the stability condition (32) is not satisfied (see [14]). This reveals that the numerical scheme (45) has better stability properties than the scheme proposed in [14].
To illustrate the convergence of scheme (45), we show the option value at = = 100 (or at = 98.01987 at time to maturity being 1 year), the differences between successive approximation processes, and the ratios between the differences in Table 1. Since scheme (45) is of order two in space and order one in time, the number of spatial mesh points is double, but the number of temporal grid points is quadruple. Table 1 only shows the numerical results for the case = 0.015, which confirm our theoretical analysis results. The numerical results for the cases of = 0 and = 0.1 are similar and we omit these here. From the data presented in Figures 1-3 and Table 1, the conclusion can be drawn: scheme (45) is unconditionally stable and convergent.
Example 2 (european butterfly spread). This example comes from [16]. A butterfly spread consists of two long positions in two calls with 1 , 2 and a short position in two calls at strike = ( 1 + 2 )/2. The payoff function ( ) is given by Let 1 = 0.8, 2 = 1.2, the volatility 0 = 0.5, the interest rate = 0.04, the maturity date = 0.5, and the artificial boundary location = 10. We first compute the numerical solution by scheme (45) with the time stepsize Δ = 0.5/ = 0.0005 and the spatial stepsize ℎ = 10/ = 0.1 similar to those in [16]. The numerical results are presented in Figure 4. These numerical data together with those in Figures 5 and  6 show that this scheme proposed here is unconditionally stable.
In [16], the authors compared the numerical results obtained by their nonstandard scheme and the schemes proposed in [14,15] and showed that their scheme produces better numerical solutions than the schemes in [14,15] with the same stepsizes. To compare our scheme with the nonstandard scheme proposed in [16], we first observe that both of the numerical schemes are unconditionally stable. Then, we will investigate their convergence and compare their convergence order. To this end, we calculate the discrete maximum norms of the errors The Scientific World Journal 9 where ,Sp and ,Ns denote the numerical solutions computed by our splitting scheme (45) and the nonstandard scheme proposed in [16] at the maturity date = 0.5, respectively, and̃denotes the numerical reference solution computed by the Backward Euler method with the standard second-order finite differences on a finer grid with = 2560 and = 320. Numerical results of both schemes (45) and the nonstandard scheme proposed in [16] are listed in Table 2. This table shows the maximum norm of the errors and the ratios between these errors. It is evident from the numerical data obtained by the splitting scheme (45) that this scheme is convergent with the temporal order one and the spatial order two. This is consistent with the theoretical results presented in this paper. For the nonstandard scheme proposed in [16], however, a reduction in the error is observed. This is not surprising since, as Zhou et al. pointed out in [16], the nonstandard scheme is conditionally consistent and the truncation error really depends on the ratio Δ /ℎ 2 .
From the theoretical analysis given in this paper and the numerical results shown in this section, we come to the following remark: the proposed scheme is efficient and reliable.

Concluding Remarks
In this paper, an unconditionally stable splitting scheme has been proposed to solve the nonlinear option pricing model with transaction costs. This method can be viewed as a "limit" version of LOD-Backward Euler method. This "limit" property that all subproblems are one-dimensional allows us to solve the discrete equation explicitly. As a consequence, this method is computationally efficient. The theoretical analysis carried out in this paper shows that this method is unconditionally stable, monotone, and positivity preserving. We also present several numerical experiments in which the vanilla call option and the European butterfly spread are considered. The theoretical analysis presented and the numerical results shown in this paper confirm that the proposed scheme here is efficient and reliable.

10
The Scientific World Journal