Abstract

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.

1. 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 and 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 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 with the terminal condition 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., [59]). 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, 1315]) 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.

2. Transformation and Spatial Semidiscretization

2.1. Preliminaries

We first consider the properties of the function appearing in (2), which will play an important role in the following numerical analysis.

Lemma 1 (see [14, 15]). The solution of ODE (3) exists and is unique, and it satisfies the following: (1) is an increasing function mapping the real line onto the interval ;(2) is implicitly defined by (3)if , then the function is bounded and where

Lemma 2 (see [14, 15]). The function is continuously differentiable at and satisfies where and are given by (8) and (9), respectively, and

Combining Lemma 2 and the results obtained in [3] leads to the following lemma.

Lemma 3. The function is a nondecreasing function of and continuously differentiable at . This implies that

2.2. The Transformed Problem

After considering the change of variable we transform the problem (4)-(5) into the following initial value problem The boundary conditions take the form of 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

2.3. The Semidiscrete Nonlinear System

To numerically approximate the solution of (14)–(16), we should consider a bounded numerical domain . Then, the problem (14)-(15) is equipped with the boundary conditions 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 , , so that . After performing the second-order central finite difference approximation of the partial derivative as we obtained the corresponding ODEs system for the semidiscrete solution as with where 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: where is the solution of the perturbation problem with the initial data .

Proof. To obtain the stability of the semidiscrete difference scheme (23), we first calculate Now, let nonlinear operator be defined by . For system (23), we have Since by Lemma 3, we obtain Then, inequality (26) follows from [18], where is the Jacobian matrix of and is the logarithmic maximum norm.

3. 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 and 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, 1315]. 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]).

3.1. Linearization System

Let us set , , 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: with where For this linearized system, we have and therefore

3.2. Numerical Method Construction

In this subsection, we construct a splitting time-stepping method based on (37). In this paper, for matrices , we define Then, the following is not true in general: Let us consider the following splitting: with where is the zero vector. Obviously, other choices for are possible; for example, , .

With the splitting (40), we solve the subproblems: starting from , and we take to complete the splitting integration step.

Then, if the Backward Euler method is used to solve every subproblem, we get 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 : 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]).

4. Properties of the Numerical Scheme

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

4.1. 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 , corresponding to the eigenvalues of , satisfies Thus, it is easy to obtain and , where denotes the spectral radius of 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]).

4.2. Positivity

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

Theorem 6. The numerical scheme (45) is unconditionally positivity preserving; that is, the solution of (45) is positive on each time level , , if is positive.

Proof. Since all entries of the matrices , , are nonnegative, a nonnegative solution of (45) on each time level can be obtained in view of the nonnegative property of and . This completes the proof.

4.3. 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).

Theorem 8. The numerical scheme (45) is unconditionally monotone.

Proof. From the proof of Theorem 6, it is known that all entries of the matrices , , are nonnegative. Then, is nondecreasing in each argument, and therefore the scheme is monotone.

4.4. Maximum Principle

From (43) and (46), it is not difficult to obtain Then, we have the maximum principle which also implies maximum norm estimate; that is, .

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

Lemma 9 (see [22]). For splitting , one has with being the commutator of and , and

Theorem 10. With the previous notation, the local discretization error of the numerical scheme (45) is given by

Proof. It follows from (37) that Because of the expansion formula we have We also have the following expansion formula similar to (61): It is easy to obtain the expansion formula of , : Comparing (63) with (64) yields, for , which implies that From (63) and (62), one gets For splitting (41), since , , we have Substituting (66) and (68) into (60) obtains Then, (59) follows from the above equation and (56).

We can now immediately conclude that the splitting scheme (45) is consistent.

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

Theorem 11. The splitting Backward Euler scheme (45) unconditionally converges to the viscosity solution of (14)-(15).

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

5. 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 , the volatility , the interest rate , the maturity date year, and the artificial boundary location . We first consider a small time stepsize and a spatial stepsize as done in [14]. The numerical results computed by the numerical scheme (45) are plotted in Figure 1 in coordinate. The case of is linear model, and the cases of and 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 and , 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 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 (or at at time to maturity being 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 , which confirm our theoretical analysis results. The numerical results for the cases of and are similar and we omit these here.

From the data presented in Figures 13 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 , and a short position in two calls at strike . The payoff function is given by Let , , the volatility , the interest rate , the maturity date , and the artificial boundary location . We first compute the numerical solution by scheme (45) with the time stepsize and the spatial stepsize 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 where and denote the numerical solutions computed by our splitting scheme (45) and the nonstandard scheme proposed in [16] at the maturity date , 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 and . 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 .

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.

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

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the anonymous referees and the editor for the important comments that led to a greatly improved paper. This work was supported by the National Natural Science Foundation of China (Grants nos. 11371074 and 11001033), the Hunan Provincial Natural Science Foundation of China (Grant no. 13JJ1020), the Research Foundation of Education Bureau of Hunan Province, China (Grant no. 13A108), and the Open Fund Project of Key Research Institute of Philosophies and Social Sciences in Hunan Universities.