A Positivity-Preserving Numerical Scheme for Nonlinear Option Pricing Models

Abstract

A positivity-preserving numerical method for nonlinear Black-Scholes models is developed in this paper. The numerical method is based on a nonstandard approximation of the second partial derivative. The scheme is not only unconditionally stable and positive, but also allows us to solve the discrete equation explicitly. Monotone properties are studied in order to avoid unwanted oscillations of the numerical solution. The numerical results for European put option and European butterfly spread are compared to the standard finite difference scheme. It turns out that the proposed scheme is efficient and reliable.

1. Introduction

It is widely recognized that the value of a European option can be obtained by solving the linear Black-Scholes equation under quite restrictive assumptions (such as liquid, frictionless, and complete markets) [1, 2]. However, these restrictive assumptions are never fulfiled in reality. In order to conform the actual situation, many modified Black-Scholes models have been proposed in recent years, such as transaction costs (Leland [3], Palmer [4], Hoggard et al. [5], Barles and Soner [6], and Jandačka and Ševčovič [7]), illiquid market (Frey and Patie [8], Sircar and Papanicolaou [9], Liu and Yong [10], etc.), and volatility uncertainty (Avellaneda et al. [11]). These models result in quasilinear or fully nonlinear Black-Scholes equations.

In this paper, we are interested in the option pricing model with transaction costs proposed by Barles and Soner [6] that are motivated by Hodges and Neuberger [12]. In their model, the value 𝑉(𝑆,𝑡) of the option satisfies the following partial differential equation: 𝜕𝑉+1𝜕𝑡2𝜎𝑉_𝑆𝑆²𝑆²𝜕²𝑉𝜕𝑆²+𝑟𝑆𝜕𝑉𝜕𝑆−𝑟𝑉=0 (1) with the nonlinear volatility 𝜎(𝑉_𝑆𝑆) that reads 𝜎𝑉_𝑆𝑆²=𝜎²₀1+Ψe^{𝑟(𝑇−𝑡)}𝑎²𝑆²𝜕²𝑉𝜕𝑆² (2) and the terminal condition 𝑉(𝑆,𝑇)=𝑓(𝑆), (3) where 𝑆 is the price of the underlying asset, 𝑇 is the maturity date, 𝑟 is the risk-free interest rate, 𝜎₀ is the asset volatility, and 𝑎=𝜇𝛾𝑁, with the proportional transaction cost 𝜇, the risk aversion factor 𝛾, and the number 𝑁 of options to be sold. Ψ(𝑥) is the solution of the following ordinary differential equation: Ψ^′(𝑥)=Ψ(𝑥)+12(𝑥)(0)𝑥Ψ−𝑥,𝑥≠0,Ψ=0. (4)

The existence and uniqueness of the solution of (1.1)–(1.3) have been shown by the theory of stochastic optimal control in [6]. However, analytical solutions cannot be found because of fully nonlinear properties of (1.1); thus, we need to compute the option values numerically.

There have been rich achievements for the numerical method of linear Black-Scholes equations (e.g., see [13–17]). As for the nonlinear situation, only a few results can be found. There is a stable numerical scheme developed in [7] (see also [18] for application to a general class of nonlinear Black-Scholes equations) for the so-called gamma equation (a quasilinear parabolic equation for the 𝑆𝑉_𝑆𝑆). Recently Kútik and Mikula [19] did some progress in showing its stability and accuracy for nonlinear Black-Scholes equation. Company et al. [20–22] construct explicit finite difference schemes for (1.1)–(1.3), and consistency and stability are studied. However, they have the disadvantage that strictly restrictive conditions on the discretization parameters are needed to guarantee stability and positivity. The implicit schemes do not have this disadvantage, but they are quite time-consuming. Yousuf et al. [23] develop a new second-order exponential time differencing (ETD) scheme to avoid unwanted oscillations near the non-smooth nodes for the Hoggard-Whalley-Wilmott (HWW) model [5] based on the Cox and Matthews approach [24] and partial fraction version of the matrix exponential, but the theoretical analysis of stability and convergence are not studied. Some authors (for instance, Düring et al. [25] for European options, Dremkova and Ehrhardt [26] for American options) construct high-order compact difference schemes with frozen values of the nonlinear coefficient of the nonlinear Black-Scholes equation to make the scheme linear and show that the resulted linearized problem is stable.

On the other hand, since the value of option is nonnegative, it is very important to make numerical schemes preserve the positivity of solution. Several authors have developed some schemes that guarantee the positivity of solutions for ordinary differential equations [27, 28] and parabolic equations [29]. In [30], Chen-Charpentier and Kojouharov propose an unconditionally positivity-preserving scheme for linear advection-diffusion reaction equations. They construct a nontraditional discretization of the advection and diffusion terms by the approximation of the spatial derivatives using values at different time levels. Motivated by this work, we will develop the method to a nonlinear Black-Scholes equation, and some properties (such as stability, monotonicity, and consistency,) of numerical scheme are studied in this paper. The new numerical method unconditionally preserves the positivity of the solutions, the stability, and monotonicity of the scheme. In addition, the designed numerical approximations allow us to solve the discrete equation explicitly, which reduces the time of calculation and increases the efficiency of the methods.

The rest of the paper is organized as follows. In the next section the original problem (1.1)–(1.3) is transformed into a nonlinear diffusion problem by an appropriate change of variables and some properties of the function Ψ(𝑥) are given. In Section 3, the discretization method is constructed. In Section 4, we prove the boundedness of coefficients, positivity, and monotonicity of the numerical scheme. Stability and consistency are studied in Section 5. In Section 6, numerical experiments for European put option and a European butterfly spread are presented to support these theoretical results. Finally, some conclusions are drawn in Section 7.

2. The Transformed Problem

For the convenience in the numerical processing and the study of the numerical analysis, we are going to transform the problem (1.1)–(1.3) into a nonlinear diffusion equation. Taking the variable transformation 𝑥=exp(𝑟(𝑇−𝑡))𝑆,𝜏=𝑇−𝑡,𝑢=exp(𝑟(𝑇−𝑡))𝑉. (5) the original problem (1.1)–(1.3) is transformed into ℒ(𝑢)=𝜕𝑢−1𝜕𝜏2𝜎²₀1+Ψ𝑎²𝑥²𝜕²𝑢𝜕𝑥²𝑥²𝜕²𝑢𝜕𝑥²=0,𝑥∈[0,+∞),𝜏∈[0,𝑇], (6) with the initial condition 𝑢(𝑥,0)=𝑓(𝑥) (7) and the boundary condition Calloption∶𝑢(0,𝑡)=0,lim_𝑥→∞𝑢(𝑥,𝜏)𝑥()=1,Putoption∶𝑢0,𝑡=𝐾,lim_𝑥→∞𝑢(𝑥,𝜏)=0. (8)

The following two lemmas give the properties of the function Ψ appearing in (1.2), which will play an important role in the numerical analysis and numerical calculation.

Lemma 2.1 (see [20]). The solution Ψ of ordinary differential equation (1.4) exists and is unique, and it satisfies,(i)Ψ is an increasing function mapping the real line onto the interval (−1,+∞).(ii)Ψ=Ψ(𝐴) is implicitly defined by −√𝐴=arcsinhΨ√+√Ψ+1Ψ²√𝑖𝑓Ψ>0,𝐴=−arcsin−Ψ√−√Ψ+1−Ψ²𝑖𝑓−1<Ψ<0, (9)(iii)if 𝐴>0, then the function Ψ(𝐴) is bounded and satisfies 0<Ψ(𝐴)<Ψ^′𝐴₂𝐴+𝑑₂, (10)where 𝐴₂=2sinh2−(sinh2)²+1²≃9.58,Ψ𝐴₂=(sinh2)²,Ψ^′𝐴₂=e⁸+2e⁴+1²e¹⁶−66e⁸+1≃1.10,𝑑₂=Ψ𝐴₂−Ψ^′𝐴₂𝐴₂≃2.62. (11)

Lemma 2.2 (see [20]). Let 𝑔(𝐴)=𝐴Ψ(𝐴), then 𝑔(𝐴) is continuously differentiable at 𝐴=0 and satisfies 𝑔^′(𝐴)≤max𝐺,2|𝐴|Ψ^′𝐴₂+𝑑₂,𝐴∈ℝ, (12) where 𝐴₂ and 𝑑₂ are given by (2.7), and 𝐴₁√=−4𝜋−33²𝑔36,𝐺=max^′(𝐴);𝐴₁≤𝐴≤𝐴₂. (13)

3. The Unconditionally Positivity-Preserving Scheme

Since the value of an option is nonnegative, it is important that numerical scheme is positivity preserving.

We see that the problem (2.2) is described in an infinite domain ℝ⁺×[0,𝑇], which makes it difficult to construct scheme effectively. Let us consider the truncated numerical domain Ω=[0,𝐵]×[0,𝑇] and discretize it in the following form. We introduce a grid of mesh points (𝑥,𝜏)=(𝑥_𝑖,𝜏_𝑛), where 𝑥_𝑖=𝑖ℎ,𝜏_𝑛=𝑛𝑘,𝑖=0,1,…,𝑁,𝑛=0,1,…,𝐿, and the spatial step size given by ℎ=𝐵/𝑁, and the time step size is 𝑘=𝑇/𝐿. Let us denote the approximation of 𝑢(𝑥_𝑖,𝜏_𝑛) by 𝑢^𝑛_𝑖 and define the following finite difference approximations of derivatives: 𝜕𝑢𝑥𝜕𝜏^𝑖,𝜏_𝑛≈𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖𝑘,𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛≈𝑢^𝑛_𝑖+1−2𝑢_𝑖^𝑛+1+𝑢^𝑛_𝑖−1ℎ²=𝛿^𝑛_𝑖. (14)

Clearly, the approximation used to calculate the second partial derivative of 𝑢 with respect to 𝑥 is nonstandard. From (2.2), we can obtain the positivity-preserving finite-difference numerical scheme as follows: 𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖𝑘−12𝜎²₀1+Ψ𝑎²𝑥²_𝑖𝛿^𝑛_𝑖𝑥²_𝑖𝑢^𝑛_𝑖+1−2𝑢_𝑖^𝑛+1+𝑢^𝑛_𝑖−1ℎ²=0. (15)

Scheme (3.2) is equivalent to 𝑢_𝑖^𝑛+1=𝜌𝛼^𝑛_𝑖𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖−1+𝑢^𝑛_𝑖1+2𝜌𝛼^𝑛_𝑖, (16) where 𝑘𝜌=ℎ²,𝛼^𝑛_𝑖=12𝜎²₀1+Ψ𝑎²𝑥²_𝑖𝛿^𝑛_𝑖𝑥²_𝑖. (17)

Remark 3.1. From property (i) of Lemma 2.1, Ψ takes values in the interval (−1,+∞), so the coefficients are nonnegative, that is, 𝛼^𝑛_𝑖≥0 for any 𝑖,𝑛.
Obviously, numerical scheme (3.3) is unconditionally positive for a nonnegative payoff 𝑢⁰_𝑖. However, we cannot obtain the numerical solution explicitly since the nonlinear term 𝛼^𝑛_𝑖 involves the value 𝑢 at the time level 𝑛+1, which makes it quite difficult to prove the stability of the scheme presented previously. In fact, the numerical scheme (3.3) can only be solved by a nonlinear iteration in each time step which is quite time-consuming.
In order to obtain an efficient scheme, we correct the approximation of the nonlinear coefficients in (2.2) by using the standard second-order central difference. Thus the corrected numerical scheme is as follows: 𝑢_𝑖^𝑛+1=𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖−1+𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖, (18) where 𝛽^𝑛_𝑖=12𝜎²₀1+Ψ𝑎²𝑥²_𝑖Δ^𝑛_𝑖𝑥²_𝑖,Δ^𝑛_𝑖=𝑢^𝑛_𝑖+1−2𝑢^𝑛_𝑖+𝑢^𝑛_𝑖−1ℎ². (19)
Since the calculation of (3.5) for 𝑖=0 and 𝑖=𝑁 requires us to know the fictitious values 𝑢^𝑛₋₁ and 𝑢^𝑛_𝑁+1, we obtain them using the linear extrapolation as follows: 𝑢^𝑛₋₁=2𝑢^𝑛₀−𝑢^𝑛₁,𝑢_𝑁+1=2𝑢^𝑛_𝑁−𝑢^𝑛_𝑁−1. (20)
Thus the values of boundary points turn out 𝑢₀^𝑛+1=𝑢^𝑛₀=⋯=𝑢⁰₀=𝑓(0),𝑢_𝑁^𝑛+1=𝑢^𝑛_𝑁=⋯=𝑢⁰_𝑁=𝑓(𝐵). (21)
The left boundary condition is not needed and in fact must not be prescribed in the case of a parabolic equation with degenerating diffusion term at 𝑥=0. This is known in the literature as the Fichera condition [31] (see [18] for application in Black-Scholes equations, Chapter 8, (8.25)). The Fichera condition is just a comparison of the speed of degeneration versus. speed of advection at the boundary 𝑥=0. Fortunately, the left boundary condition in (3.8) is a consequence of (3.5) with 𝑥₀=0,𝑖=0,𝑢⁰₀=𝑓(0). The only boundary condition that can be prescribed is the right boundary condition at 𝑥=𝐵. For the convenience in the numerical calculation, let us denote the vectors 𝑢^𝑛=[𝑢^𝑛₀,𝑢^𝑛₁,…,𝑢^𝑛_𝑁]^𝑇, then numerical scheme (3.5), (3.8) can be written in matrix form 𝑢^𝑛+1=𝐶(𝑛)𝑢^𝑛,𝑢⁰=𝑓(0),𝑓𝑥₁,…,𝑓(𝐵)^𝑇, (22) where 𝐶(𝑛)=𝑎1000⋯0₁𝑏₁𝑐₁0⋯00𝑎₂𝑏₂𝑐₂⋯0⋱⋱⋱0⋯0𝑎_𝑁−1𝑏_𝑁−1𝑐_𝑁−1𝑎0⋯0001_𝑖,=𝑐_𝑖=𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖,𝑏_𝑖=11+2𝜌𝛽^𝑛_𝑖,𝑖=1,2,…,𝑁−1. (23)(24)
In order to study the related properties of numerical scheme (such as monotonicity and stability), we need to know the behaviour of Δ^𝑛_𝑖 appearing in the nonlinear team 𝛽^𝑛_𝑖.

Lemma 3.2. With the previous notation, Δ^𝑛_𝑖 appearing in the nonlinear team 𝛽^𝑛_𝑖 satisfies the scheme Δ_𝑖^𝑛+1=𝜌𝛽^𝑛_𝑖+11+2𝜌𝛽^𝑛_𝑖+1Δ^𝑛_𝑖+1+𝜌𝛽^𝑛_𝑖−11+2𝜌𝛽^𝑛_𝑖−1Δ^𝑛_𝑖−1+11+2𝜌𝛽^𝑛_𝑖Δ^𝑛_𝑖Δ,1≤𝑖≤𝑁−1,_𝑛0=Δ^𝑛_𝑁=0,𝑛=0,1,2,…,𝐿. (25)

Proof. From (3.6) and (3.7) we can know for 𝑖=0 and 𝑖=𝑁 that Δ^𝑛₀=Δ^𝑛_𝑁=0,𝑛=0,1,2,…,𝐿. (26)
Putting (3.5) into the expression (3.6) of Δ^𝑛_𝑖, and after a simple calculation we can obtain (3.12).
On the other hand, the most common option price sensitivities are the first and second derivatives with respect to the price of the underlying asset, that is, “delta” and “gamma”, respectively. These are important features in risk management and are challenging to compute numerically. From the transformation (2.1) we can obtain the approximations of Greeks of the option as follows: 𝜕Gamma=²𝑉𝜕𝑆²𝑆_𝑖,𝑡_𝑛=exp(𝑟(𝑇−𝑡))Δ^𝑛_𝑖(𝑢),Delta=𝜕𝑉𝑆𝜕𝑆^𝑖,𝑡_𝑛(𝑟(𝑢=exp𝑇−𝑡))^𝑛𝑖+1−𝑢^𝑛_𝑖−1.2ℎ (27)

4. Properties of the Numerical Scheme

4.1. Boundedness of the Coefficients

For the sake of convenience, we introduce the following definition.

Definition 4.1 (see [32]). If 𝑥=[𝑥₁,𝑥₂,…,𝑥_𝑁] is a vector in ℝ^𝑁, then its 1-norm is denoted by ‖𝑥‖₁=∑^𝑁_𝑖=10𝑥0200𝑑|𝑥_𝑖|, and maximum-norm is denoted by ‖𝑥‖_∞=max_{1≤𝑖≤𝑁}|𝑥_𝑖|.
The following theorem shows that the nonlinear team 𝛽^𝑛_𝑖 appearing in (3.5) is bounded.

Theorem 4.2. Let Δ^𝑛=[Δ^𝑛₀,Δ^𝑛₁,…,Δ^𝑛_𝑁]^𝑇, then the following properties hold true.(i)‖Δ^𝑛‖₁ is nonincreasing.(ii)The nonlinear team 𝛽^𝑛_𝑖 appearing in (3.5) satisfies 0≤𝛽^𝑛_𝑖≤𝐿(ℎ),𝑖=0,1,…,𝑁,𝑛=0,1,…,𝐿, (28)where 𝐿(ℎ)=𝐵²𝜎²₀21+𝑑₂+𝑎²𝐵²Δ⁰₁Ψ^′𝐴₂, (29) with 𝐴₂, 𝑑₂, and Ψ^′(𝐴₂) given by (2.7).

Proof. Property (i) is proved using the induction principle over index 𝑛.
For 𝑛=0, from Lemma 3.2 and 𝛽^𝑛_𝑖≥0 by Remark 4.7, it follows that Δ¹_𝑖≤𝜌𝛽⁰_𝑖+11+2𝜌𝛽⁰_𝑖+1Δ⁰_𝑖+1+𝜌𝛽⁰_𝑖−11+2𝜌𝛽⁰_𝑖−1Δ⁰_𝑖−1+11+2𝜌𝛽⁰_𝑖Δ⁰_𝑖,1≤𝑖≤𝑁−1. (30) Taking into account (3.13) and (4.3), we have Δ¹₁≤^𝑁−1_𝑖=10𝑥0200𝑑Δ⁰_𝑖−𝜌𝛽⁰₁1+2𝜌𝛽⁰₁Δ⁰₁−𝜌𝛽⁰_𝑁−11+2𝜌𝛽⁰_𝑁−1Δ⁰_𝑁−1≤Δ⁰₁. (31) Thus property (i) is proved for 𝑛=0.
Now, let us assume that property (i) holds true up 𝑛, that is, Δ^𝑛₁≤Δ^𝑛−1₁≤Δ^𝑛−2₁≤⋯≤Δ⁰₁. (32) For 1≤𝑖≤𝑁−1, from Lemma 3.2, it follows that Δ_𝑖^𝑛+1≤𝜌𝛽^𝑛_𝑖+11+2𝜌𝛽^𝑛_𝑖+1Δ^𝑛_𝑖+1+𝜌𝛽^𝑛_𝑖−11+2𝜌𝛽^𝑛_𝑖−1Δ^𝑛_𝑖−1+1−2𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖Δ^𝑛_𝑖. (33) Taking into account Δ₀^𝑛+1=Δ_𝑁^𝑛+1=0, we have Δ^𝑛+1₁≤^𝑁_𝑖=20𝑥0200𝑑𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖Δ^𝑛_𝑖+^𝑁−2_𝑖=00𝑥0200𝑑𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖Δ^𝑛_𝑖+^𝑁−1_𝑖=10𝑥0200𝑑1−2𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖Δ^𝑛_𝑖=^𝑁−1_𝑖=10𝑥0200𝑑Δ^𝑛_𝑖−𝜌𝛽^𝑛₁1+2𝜌𝛽^𝑛₁Δ^𝑛₁−𝜌𝛽^𝑛_𝑁−11+2𝜌𝛽^𝑛_𝑁−1Δ^𝑛_𝑁−1≤Δ^𝑛₁. (34) Hence property (i) is proved completely.
On the other hand, from (3.6) and the monotonic increasing property of Ψ (see Lemma 2.1), it follows that 𝛽^𝑛_𝑖=12𝜎²₀1+Ψ𝑎²𝑥²_𝑖Δ^𝑛_𝑖𝑥²_𝑖≤12𝜎²₀1+Ψ𝑎²𝑥²_𝑖Δ^𝑛₁𝑥²_𝑖≤12𝜎²₀1+Ψ𝑎²𝑥²_𝑖Δ⁰₁𝑥²_𝑖. (35)
Since 𝑥_𝑖∈[0,𝐵] and from the property (iii) and property (i) of Lemma 2.2, we have 0≤𝛽^𝑛_𝑖≤𝐿(ℎ),𝑖=0,1,…,𝑁,𝑛=0,1,…,𝐿. (36) Thus the proof of theorem is complete.

4.2. Positivity

Since the value of option is nonnegative, a nice property of the numerical scheme for the pricing equation is positivity-preserving.

Clearly, all the coefficients of the numerical scheme (3.5) are nonnegative (see Theorem 4.2). Hence, for a nonnegative payoff 𝑢⁰_𝑖, the following result has been established.

Proposition 4.3. The numerical scheme (3.5), (3.8) is unconditionally positive.

Remark 4.4. Δ^𝑛_𝑖 is also unconditionally positive, that is, if Δ⁰_𝑖≥0 for 0≤𝑖≤𝑁, then Δ^𝑛_𝑖≥0 for 0≤𝑖≤𝑁,0≤𝑛≤𝐿.

4.3. Monotonicity

For the sake of convenience in the presentation, we introduce the following definition of a monotonicity preserving numerical scheme.

Definition 4.5 (see [31]). In numerical scheme 𝐹(𝑢^𝑛_𝑖)=0, we say that it is monotonicity-preserving. If each time that 𝑢^𝑛_𝑖≤𝑢^𝑛_𝑖+1 or 𝑢^𝑛_𝑖≥𝑢^𝑛_𝑖+1 for all 𝑖, then it occurs that 𝑢_𝑖^𝑛+1≤𝑢^𝑛+1_𝑖+1 or 𝑢_𝑖^𝑛+1≥𝑢^𝑛+1_𝑖+1 for all 𝑖.
The next result shows the monotonicity of the numerical scheme.

Theorem 4.6. The numerical scheme (3.5), (3.8) is unconditionally monotonicity-preserving with 0≤𝑖≤𝑁,0≤𝑛≤𝐿.

Proof. Let us write 𝑢^𝑛+1_𝑖+1−𝑢_𝑖^𝑛+1=𝑢^𝑛+1_𝑖+1−𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖−𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖. (37) Assuming that 𝑢^𝑛_𝑖+1≥𝑢^𝑛_𝑖 for 0≤𝑖≤𝑁−1,0≤𝑛≤𝐿−1, then from (3.5), it follows that 𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖=𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖−1−2𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖≤𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖𝑢,1≤𝑖≤𝑁−1,_𝑖𝑛+1−𝑢^𝑛_𝑖≥−𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖−𝑢^𝑛_𝑖−11+2𝜌𝛽^𝑛_𝑖,1≤𝑖≤𝑁−1, (38)(39) that substituting 𝑖 by 𝑖+1 one gets 𝑢^𝑛+1_𝑖+1−𝑢^𝑛_𝑖+1≥−𝜌𝛽^𝑛_𝑖+1𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖+1,0≤𝑖≤𝑁−2. (40) From (4.10), (4.11), and (4.13), it follows that 𝑢^𝑛+1_𝑖+1−𝑢_𝑖^𝑛+1≥−𝜌𝛽^𝑛_𝑖+1𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖+1+𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖−𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖=1−𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖−𝜌𝛽^𝑛_𝑖+11+2𝜌𝛽^𝑛_𝑖+1𝑢^𝑛_𝑖+1−𝑢^𝑛_𝑖. (41) Taking into account (4.1) and (4.14) for 𝑖 and 𝑖+1, we have 𝑢^𝑛+1_𝑖+1−𝑢_𝑖^𝑛+1≥(ℎ)1−2𝜌𝐿(ℎ)𝑢1+2𝜌𝐿^𝑛𝑖+1−𝑢^𝑛_𝑖≥0,1≤𝑖≤𝑁−2. (42)
The monotonicity of the scheme in the internal mesh points has been proved.
In an analogous way, we can verify that 𝑢₁^𝑛+1≥𝑢₀^𝑛+1,𝑢_𝑁^𝑛+1≥𝑢^𝑛+1_𝑁−1. (43)
Similarly, we can prove that if 𝑢^𝑛_𝑖≥𝑢^𝑛_𝑖+1, then 𝑢_𝑖^𝑛+1≥𝑢^𝑛+1_𝑖+1 for 0≤𝑖≤𝑁−1,0≤𝑛≤𝐿−1.
Thus the proof of theorem is complete.

Remark 4.7. If the payoff 𝑓(𝑥) is nondecreasing with 𝑓(0)=0 (e.g., 𝑓(𝑥)=max{𝑥−𝐾,0}), then 0=𝑢^𝑛₀≤𝑢^𝑛₁≤⋯≤𝑢^𝑛_𝑖≤𝑢^𝑛_𝑖+1≤⋯≤𝑢^𝑛_𝑁 for a fixed 𝑛 with 0≤𝑛≤𝐿.

5. Stability and Consistency

5.1. Stability

Theorem 5.1. The difference scheme (3.5), (3.8) is unconditionally ‖⋅‖_∞-stable.

Proof. From (3.5), let us write 𝑢_𝑖^𝑛+1=𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖−1+1−2𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖. (44) Since all the coefficients of (5.1) are nonnegative, then using triangle inequality, it follows that 𝑢_𝑖^𝑛+1≤𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1+𝑢^𝑛_𝑖−1+1−2𝜌𝛽^𝑛_𝑖1+2𝜌𝛽^𝑛_𝑖𝑢^𝑛_𝑖≤sup_{0≤𝑖≤𝑁}𝑢^𝑛_𝑖, (45) so that sup_{0≤𝑖≤𝑁}𝑢_𝑖^𝑛+1≤sup_{0≤𝑖≤𝑁}𝑢^𝑛_𝑖. (46) According to the definition of the maximum norm, it follows that 𝑢^𝑛+1_∞≤𝑢^𝑛_∞. (47) Thus the result is established.

5.2. Consistency

The proposed new difference scheme (3.5) is explicit, unconditionally positive, and unconditionally stable; however, it is not unconditionally consistent. There are extra truncation error terms since the approximations to second derivatives with respect to 𝑥 are at different time levels.

The following theorem gives the consistency condition of the difference scheme (3.5).

Theorem 5.2. With the previous notation, suppose that the exact solution 𝑢 of (1.1)–(1.3) satisfies 𝑢∈𝐶^4,2(Ω). Then the local truncation error is given by 𝑇_ℎ,𝑘𝑢^𝑛_𝑖=𝑂ℎ²(𝑘)𝑘+𝑂+𝑂ℎ². (48)

Proof. Let us write the scheme (3.5) in the form 𝐹_ℎ,𝑘𝑢^𝑛_𝑖=𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖𝑘−𝛽^𝑛_𝑖𝛿^𝑛_𝑖=0. (49) Using Taylor's expansion about (𝑥_𝑖,𝜏_𝑛) and 𝑢∈𝐶^4,2(Ω), it follows that 𝛿^𝑛_𝑖=𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛+ℎ²𝜕12⁴𝑢𝜕𝑥⁴𝜂,𝜏_𝑛−2𝑘ℎ²𝜕𝑢𝑥𝜕𝜏^𝑖=𝜕,𝜉₂𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛+ℎ²𝑀^𝑛_𝑖(1)−𝑘ℎ²𝑀^𝑛_𝑖(2),Δ^𝑛_𝑖=𝑢^𝑛_𝑖+1−2𝑢_𝑖^𝑛+1+𝑢^𝑛_𝑖−1ℎ²=𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛+ℎ²𝜕12⁴𝑢𝜕𝑥⁴𝜂,𝜏_𝑛,𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖𝑘=𝜕𝑢𝑥𝜕𝜏^𝑖,𝜏_𝑛+𝑘2𝜕²𝑢𝜕𝜏²𝑥_𝑖,𝜉=𝜕𝑢𝑥𝜕𝜏^𝑖,𝜏_𝑛+𝑘𝑀^𝑛_𝑖(3), (50)(51)(52) where 𝑥_𝑖−ℎ<𝜂<𝑥_𝑖+ℎ,𝜏_𝑛<𝜉<𝜏_𝑛𝑀+𝑘,_𝑛𝑖(1≤1)𝜕12max⁴𝑢𝜕𝑥⁴𝑥,𝜏_𝑛;0≤𝑥≤𝐵=𝑊^𝑛(1)_max,𝑀^𝑛_𝑖(2)≤2max𝜕𝑢𝑥𝜕𝜏^𝑖,𝜏;𝜏_𝑛≤𝜏≤𝜏_𝑛+1=𝑊^𝑛_𝑖(2)_max,𝑀^𝑛_𝑖(3≤1)2𝜕max²𝑢𝜕𝜏²𝑥_𝑖,𝜏;𝜏_𝑛≤𝜏≤𝜏_𝑛+1=𝑊^𝑛_𝑖(3)_max. (53) From (5.5)–(5.8), it follows that the local truncation error takes the form 𝑇_ℎ,𝑘𝑢^𝑛_𝑖=𝐹_ℎ,𝑘𝑢^𝑛_𝑖−ℒ𝑢^𝑛_𝑖1=−2𝜎²₀𝑥²_𝑖1+Ψ𝑎²𝑥²_𝑖Δ^𝑛_𝑖𝛿^𝑛_𝑖−1+Ψ𝑎²𝑥²_𝑖𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛+𝑘𝑀^𝑛_𝑖(3). (54) Let us introduce the notation 𝐴^𝑛_𝑖=𝑎²𝑥²_𝑖𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛,𝛿₁𝐴^𝑛_𝑖=𝑎²𝑥²_𝑖ℎ²𝑀^𝑛_𝑖(1),𝛿₂𝐴^𝑛_𝑖=−𝑎²𝑥²_𝑖𝑘ℎ²𝑀^𝑛_𝑖(2). (55) Using function 𝑔(𝐴) introduced in Lemma 2.2, it follows that 1+Ψ𝑎²𝑥²_𝑖Δ^𝑛_𝑖𝛿^𝑛_𝑖−1+Ψ𝑎²𝑥²_𝑖𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛𝜕²𝑢𝜕𝑥²𝑥_𝑖,𝜏_𝑛=𝑎⁻²𝑥_𝑖⁻²𝑔𝐴^𝑛_𝑖+𝛿₁𝐴^𝑛_𝑖−𝑔𝐴^𝑛_𝑖+𝛿₁𝐴^𝑛_𝑖+1+Ψ𝐴^𝑛_𝑖+𝛿₁𝐴^𝑛_𝑖𝛿₂𝐴^𝑛_𝑖=𝑔^′𝐴^𝑛_𝑖+𝜃𝛿₁𝐴^𝑛_𝑖ℎ²𝑀^𝑛_𝑖(1)+ℎ²𝑀^𝑛_𝑖(1)−𝑘ℎ²𝛽^𝑛_𝑖,0<𝜃<1. (56) From Lemma 2.2, (5.11), and 𝛽^𝑛_𝑖≤𝐿(ℎ), it is easy to know that 𝑇_ℎ,𝑘𝑢^𝑛_𝑖≤𝑂ℎ²(𝑘)𝑘+𝑂+𝑂ℎ². (57) Thus the result has been proved.

From Theorem 5.2, we can see that the meshes should satisfy 𝑘/ℎ²→0 as 𝑘 and ℎ go to zero in order to ensure the consistency. Therefore, the key to the convergence of the scheme is the consistency rather than the stability. In actual calculation, we can choose the time step depending on the spatial size so that inconsistent terms go to zero.

6. Numerical Experiments

In this section, we implement the positivity-preserving scheme (3.5), (3.8) on European put option and European butterfly spread. We analyse the effectiveness of the method and compare it with the numerical scheme (SFD1) given in [20] and a standard forward Euler finite difference scheme (SFD2) 𝑢_𝑖^𝑛+1−𝑢^𝑛_𝑖𝑘−𝛽^𝑛_𝑖𝑢^𝑛_𝑖+1−2𝑢^𝑛_𝑖+𝑢^𝑛_𝑖−1ℎ²=0, (58) which is analysed in [21]. The function Ψ is calculated with (2.5) using the “fsolve” function in Matlab. Both the experiments are performed for large values of 𝜎₀ to visualize the sensitivity of the methods towards high volatility.

6.1. European Put Option

A European put option is a contract where the owner of the option has the right to sell an underlying asset 𝑆(𝑡) for a fixed amount, known as the strike price 𝐾, at the expiry date 𝑇. The payoff function 𝑓(𝑆) is given by 𝑓(𝑆)=max{𝐾−𝑆,0}. (59) We choose the parameters as 𝐾=2,𝑇=0.5,𝜎₀=0.5,𝐵=10,𝑟=0.04. (60)

Equations (1.1)–(1.3) give analytical solution for only 𝑎=0 (see [6]). Figure 1 (the left one: 𝑎=0, the right one: 𝑎=0.02) gives the option value using scheme (3.5) and schemes (SFD1 and SFD2) in [20, 21], which shows that our scheme stable, monotonous, and is able to produce solution that is close to the exact solution, but numerical solution of scheme (6.1) appears as spurious oscillation for 𝑘=0.0005 and ℎ=0.1.

(a) 𝑎=0

(b) 𝑎=0.02

Figure 1 

The values of European put option under different schemes with the computational parameters ℎ=0.1 and 𝑘=0.0005.

Moreover, Figure 2 presents the related hedging parameters of European put option using the two schemes. We can see that our proposed scheme produces smoother solutions than standard scheme for the delta and gamma. In addition, the gamma is positivity preserving and is maximal as it closes the strike price 𝐾=2.

(a) Standard scheme for delta  

(b) New scheme for delta

(c) Standard scheme for gamma

(d) New scheme for gamma

Figure 2 

Time evolution profiles of the Greeks of European put option with the computational parameters ℎ=0.1, 𝑘=0.0005, and 𝑎=0.02.

Next we consider the influence of transaction costs (parameter 𝑎) on the value of European put option (Figure 3). The left one shows the change of the option value with the parameter 𝑎, and the right one presents an evolution profile of the difference 𝑉_nonlinear(𝑆,𝑡)−𝑉_linear(𝑆,𝑡) between our proposed scheme with transaction costs and without transaction costs. We can see that the difference is not symmetric, but decreases towards the expiry date, and is maximal close to the strike price 𝐾=2, where the nonlinear value is significantly higher than the linear value.

(a) Linear and nonlinear cases  

(b) Nonlinear-linear (𝑎=0.02)

Figure 3 

Influence of the transaction costs (parameter 𝑎) on the value of European put option with the computational parameters ℎ=0.1 and 𝑘=0.0005.

6.2. European Butterfly Spread

A butterfly spread is a combination of three-call options with three-strike prices, in which one contract is purchased with two outside strike prices and two contracts are sold at the middle strike price. The payoff function 𝑓(𝑆) is given by 𝑓(𝑆)=max𝑆−𝐾₁,0−2max𝑆−𝐾₂,0+max𝑆−𝐾₃,0, (61) where 𝐾₁,𝐾₂, and 𝐾₃ are the strike prices that satisfy 𝐾₁<𝐾₂<𝐾₃ and 𝐾₂=(𝐾₁+𝐾₃)/2.

We choose the following parameters: 𝑟=0.04,𝐾₁=0.8,𝐾₂=1,𝐾₃𝜎=1.2,𝑇=0.5,₀=0.5,𝐵=10, (62) the time step 𝑘=0.000455, and the spatial step ℎ=0.1.

Figure 4 shows the different solutions of two schemes at 𝑡=0 for 𝑎=0.02. We can see that our proposed scheme is smooth and stable at the same step sizes. Moreover, Figure 5 shows the Greeks of the numerical solution calculated with scheme (3.14), which is different from the vanilla option. The standard scheme (the left one) is unstable and produces unwanted oscillations which are not present while using our proposed scheme (the right one).

Figure 4 

The value of European butterfly spread under different schemes with the computational parameters ℎ=0.1, 𝑘=0.000455, and 𝑎=0.02.

(a)   Standard scheme for delta

(b) New scheme for delta

(c) Standard scheme for gamma

(d) New scheme for gamma

Figure 5 

Time evolution profiles of the Greeks of European butterfly spread with the computational parameters ℎ=0.1, 𝑘=0.000455, and 𝑎=0.02.

The last two figures show the value of European butterfly spread with different transaction costs (parameter 𝑎). It is evident from Figure 6 that butterfly spread becomes more expensive in the presence of transaction cost and the difference is maximal as it closes the strike price 𝐾₂=1, which is similar with the vanilla option.

(a)   Linear and nonlinear cases

(b) Nonlinear-linear (𝑎=0.02)

Figure 6 

Influence of the transaction costs (parameter 𝑎) on the value of European butterfly spread with the computational parameters ℎ=0.1 and 𝑘=0.000455.

7. Discussions and Conclusions

In this paper, we have extended the numerical method in [30] to nonlinear situation and presented the numerical scheme for a nonlinear Black-Scholes equation in the presence of transaction costs. The numerical method is based on a nonstandard approximation of the second partial derivative 𝜕²𝑢(𝑥_𝑖,𝜏_𝑛)/𝜕𝑥² by (𝑢^𝑛_𝑖+1−2𝑢_𝑖^𝑛+1+𝑢^𝑛_𝑖−1)/ℎ², and the nonlinear team is treated explicitly, which guarantees to solve the original problem without iteration. The scheme is unconditionally positive and stable, but it is conditionally consistent. In fact, as 𝑢_𝑖^𝑛+1≈𝑢^𝑛_𝑖+𝑘(𝜕𝑢/𝜕𝜏)(𝑥_𝑖,𝜏_𝑛), the scheme effectively solves the nonlinear parabolic equation 𝑘1+ℎ²𝛽(𝑥)𝜕𝑢−1𝜕𝜏2𝛽(𝑥)𝜕²𝑢𝜕𝑥²=0, (63) where 𝛽(𝑥)=𝜎²₀1+Ψ𝑎²𝑥²𝜕²𝑢𝜕𝑥²𝑥². (64) It means that the new scheme converges to a solution of (2.2) if 𝑘/ℎ²→0. Otherwise, (if 𝑘/ℎ² is fixed) it converges to a solution of (2.2) in a different time scale. In fact, it can be seen from Theorem 5.2, where the truncation error really depends on the ratio 𝑘/ℎ², which is also the reason that we consider the smaller ratio 𝑘/ℎ² in the experiment. The numerical results show that our method produces better numerical solutions than the schemes in [20, 21] with the same step sizes. In the future work, we will extend the method to the problem of American option pricing.

Acknowledgments

The authors would like to thank the anonymous referees for several suggestions for the improvement of this paper. This work is supported by the Fundamental Research Funds for the Central Universities of China (Grant no. JGK101677).