Cubic Spline Method for a Generalized Black-Scholes Equation

Abstract

We develop a numerical method based on cubic polynomial spline approximations to solve a a generalized Black-Scholes equation. We apply the implicit Euler method for the time discretization and a cubic polynomial spline method for the spatial discretization. We show that the matrix associated with the discrete operator is an M-matrix, which ensures that the scheme is maximum-norm stable. It is proved that the scheme is second-order convergent with respect to the spatial variable. Numerical examples demonstrate the stability, convergence, and robustness of the scheme.

1. Introduction

An option is a tradable financial contract whose value depends on the value of an underlying asset. The buyer of the contract obtains the right, but not the obligation, to buy or to sell an asset at a specified price on or before a maturity date. A call option provides the right to buy the underlying asset for a certain price, whereas a put option confers the right to sell the underlying asset for a certain price. A European option can only be exercised at the maturity date, while an American option can be exercised at any time prior to its maturity date. Black and Scholes [1] showed that the value of a European option is governed by a second-order parabolic differential equation with respect to the underlying asset price and time, which is known as the Black-Scholes equation. The value of an American option is determined by a linear complementarity problem or as a free boundary value problem involving the Black-Scholes differential operator [2]. It is often necessary to use numerical methods to solve these partial differential equations, as analytic solutions are not generally available.

There are several numerical methods for the valuation of European and American options. The first numerical method for option pricing was the lattice method proposed in Cox et al. [3] and was improved in Hull and White [4], which is equivalent to an explicit time-stepping scheme. Since the Black-Scholes equation with constant or space-independent parameters can be transformed into a diffusion equation, the finite difference methods applied to constant-coefficient heat equations have also been developed (see, e.g., Schwartz [5], Courtadon [6], Wilmott et al. [2], and Rogers and Talay [7]) for pricing options. Vázquez [8] presented an upwind numerical approach for the Black-Scholes equation. Cen and Le [9–11] presented stable finite difference schemes on a piecewise uniform mesh for pricing European and American options. Wang [12] and Angermann and Wang [13] proposed a fitted finite volume method for the discretization of the Black-Scholes equation. Rambeerich et al. [14] applied the exponential time integration scheme to price options. Other methods, such as meshless approach [15, 16], element-free kp-Ritz method [17, 18], and element-free Galerkin method [19], also can be used to solve the generalized Black-Scholes equation.

The spline collocation methods are useful methods for solving partial differential equations. Spline solutions have their own advantages. For example, once the solution has been computed, the information between mesh points is available. Numerical methods based on spline collocation methods also have been used to solve option pricing problems. Christara et al. [20] proposed a quadratic spline collocation method to the American option pricing problems. Holtz and Kunoth [21] developed a B-spline-based monotone multigrid method for the valuation of American options. Khabir and Patidar [22] applied a B-spline collocation method to solve the heat equation which is obtained from the Black-Scholes equation by an Euler transformation. Kadalbajoo et al. [23, 24] used cubic B-spline collocation methods for the Black-Scholes equation.

In this paper, we present a numerical method based on cubic polynomial spline to solve a generalized Black-Scholes equation. We combine the implicit Euler method for discretizing the time variable with the cubic polynomial spline scheme for discretizing the spatial variable. The matrix associated with the discrete operator is an M-matrix, which ensures that the scheme is maximum-norm stable. We will show that the scheme is second-order convergent with respect to the spatial variable. Numerical results support the theoretical results.

The rest of the paper is organized as follows. In the next section, a generalized Black-Scholes equation is introduced. The temporal semidiscretization is described in Section 3. The spatial discretization is constructed in Section 4. The fully discrete scheme is presented in Section 5. Finally, numerical experiments are provided to support these theoretical results in Section 6.

2. The Continuous Problem

We consider the following generalized Black-Scholes equation: −𝜕𝑣−1𝜕𝑡2𝜎²(𝑆,𝑡)𝑆²𝜕²𝑣𝜕𝑆²−𝑟(𝑡)𝑆𝜕𝑣(𝑡)()𝜕𝑆+𝑟𝑣=0,𝑆,𝑡∈ℝ⁺×(0,𝑇), (1) equipped with the terminal and boundary conditions 𝑣(𝑆,𝑇)=max(𝑆−𝐸,0),𝑥∈ℝ⁺,𝑣(0,𝑡)=0,𝑡∈[0,𝑇]. (2) Here, 𝑣(𝑆,𝑡) is the value of European call option at asset price 𝑆 and at time 𝑡, 𝐸 is the exercise price, 𝑇 is the maturity, 𝑟(𝑡) is the risk-free interest rate, and 𝜎(𝑆,𝑡) is the volatility function of underlying asset. When 𝜎 and 𝑟 are constant functions, it becomes the classical Black-Scholes model. The existence and uniqueness of a classical solution of (1)-(2) are well known (see [25, 26]).

Note that (1) degenerates when 𝑆 goes to zero. We transform the Black-Scholes equations (1)-(2) into a nondegenerate partial differential equation by using a log transformation 𝑥=ln𝑆−𝜕𝑢−1𝜕𝑡2𝜎²()𝜕𝑥,𝑡²𝑢𝜕𝑥²−𝑟(𝑡)−12𝜎²(𝑥,𝑡)𝜕𝑢(𝑡)()(),𝑢()𝑒𝜕𝑥+𝑟𝑢=0,𝑥,𝑡∈ℝ×0,𝑇𝑥,𝑇=max^𝑥𝑢()[],−𝐸,0,𝑥∈ℝ,𝑥,𝑡=0,as𝑥⟶−∞,𝑡∈0,𝑇 (3) where 𝑢(𝑥,𝑡)=𝑣(𝑆,𝑡)=𝑣(𝑒^𝑥,𝑡)and𝜎(𝑥,𝑡)=𝜎(𝑒^𝑥,𝑡).

For applying the numerical method, we truncate the infinite domain ℝ×(0,𝑇) into a truncated domain Ω=(𝑥_min,𝑥_max)×(0,𝑇), where 𝑥_min and 𝑥_max are chosen properly so that, for practical purposes, they do not affect the option price. Based on Willmott et al.'s estimate [2] that the upper bound of the asset price is typically three or four times the strike price, it is reasonable for us to set 𝑥_max=ln(4𝐸), while, for the lower bound of the asset price, since −ln(4𝐸) is negative enough, we take 𝑥_min=−ln(4𝐸) for convenience in numerical experiments. Therefore, in the remaining of this paper, we will consider the following problem: −𝜕𝑢−1𝜕𝑡2𝜎²()𝜕𝑥,𝑡²𝑢𝜕𝑥²−𝑟(𝑡)−12𝜎²(𝑥,𝑡)𝜕𝑢(𝑡)()𝑢()𝑒𝜕𝑥+𝑟𝑢=0,𝑥,𝑡∈Ω,𝑥,𝑇=max^𝑥−𝐸,0,𝑥∈𝑥_min,𝑥_max,𝑢𝑥_min[],𝑢𝑥,𝑡=0,𝑡∈0,𝑇_max,𝑡=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑟(𝑠)d𝑠},𝑡∈[0,𝑇]. (4) Here, the right boundary condition is chosen according to Vázquez [8]. Normally, this truncation of the domain leads to a negligible error in the value of the option [25].

3. The Temporal Semidiscretization

To approximate the solution (4), first, we apply the implicit Euler method to discretize the temporal variable. This scheme, on a uniform mesh Ω^𝐾=𝑡_𝑗𝑇=𝑗Δ𝑡,0≤𝑗≤𝐾,Δ𝑡=𝐾, (5) reads 𝑢^𝐾=𝑢(𝑥,𝑇)=max𝑒^𝑥−𝐸,0,𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗(𝑥)=𝑢^𝑗+1(𝑥),𝑢^𝑗𝑥_min𝑢=0,_𝑗𝑥_max=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑗𝑟(𝑠)d𝑠},for𝑗=𝐾−1,…,1,0, (6) where 𝐿_𝑥𝑢^𝑗(𝑥)1=−2𝜎^𝑗(𝑥)²𝑑²𝑢^𝑗𝑑𝑥²−𝑟^𝑗−12𝜎^𝑗(𝑥)²𝑑𝑢^𝑗𝑑𝑥+𝑟^𝑗𝑢^𝑗(𝑥), (7) and 𝑢^𝑗(𝑥) denotes the approximation of the exact solution 𝑢(𝑥,𝑡) at the time level 𝑡_𝑗.

Similarly to Kellogg and Tsan [27], we can prove that the differential operator (𝐼+Δ𝑡𝐿_𝑥) satisfies a maximum principle, and, consequently, 𝐼+Δ𝑡𝐿_𝑥⁻¹_∞≤1.1+𝑟Δ𝑡 (8) Hence, we can obtain the following result.

Lemma 1. The temporal semidiscretization scheme (6) is unconditionally stable.

Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problem 𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗=𝑢𝑥,𝑡_𝑗+1,𝑢^𝑗𝑥_min𝑢=0,_𝑗𝑥_max=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑗𝑟(𝑠)d𝑠}, (9) is introduced to define the local error.

Lemma 2 (local error estimate). The local error associated with the method (9), defined by 𝑒^𝑗=𝑢(𝑥,𝑡_𝑗)−𝑢^𝑗(𝑥), satisfies 𝑒^𝑗_∞=𝑂(Δ𝑡)². (10)

Proof. Using Taylor expansion, we have −𝑢𝑥,𝑡_𝑗+1−𝑢𝑥,𝑡_𝑗Δ𝑡=−𝑢_𝑡𝑥,𝑡_𝑗+𝑂(Δ𝑡)=−𝐿_𝑥𝑢𝑥,𝑡_𝑗+𝑂(Δ𝑡), (11) for 1≤𝑗<𝐾. From (9) and (11), it is straightforward to show that the local error is the solution of the problem 𝐼+Δ𝑡𝐿_𝑥𝑒^𝑗=𝑂(Δ𝑡)²,𝑒^𝑗𝑥_min=𝑒^𝑗𝑥_max=0, (12) and, therefore, the result follows from the maximum principle for the operator (𝐼+Δ𝑡𝐿_𝑥).

Lemma 3 (global error estimate). The global error associated with the implicit Euler method (6), given by 𝐸^𝑗=𝑢(𝑥,𝑡_𝑗)−𝑢^𝑗(𝑥), satisfies 𝐸=sup_𝑗≤𝐾𝐸^𝑗_∞=𝑂(Δ𝑡), (13) and, therefore, the temporal semidiscretization scheme is a first-order convergent scheme.

Proof. The global truncation error at time 𝑡_𝑗 can be decomposed as 𝐸^𝑗=𝑢𝑥,𝑡_𝑗−𝑢^𝑗(𝑥)=𝑢𝑥,𝑡_𝑗−𝑢^𝑗(𝑥)+𝑢^𝑗(𝑥)−𝑢^𝑗(𝑥). (14) By relations (6) and (9), we have 𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗(𝑥)−𝑢^𝑗(𝑥)=𝑢𝑥,𝑡_𝑗+1−𝑢^𝑗+1(𝑥). (15) Applying the maximum principle for the operator (𝐼+Δ𝑡𝐿_𝑥), we can obtain 𝑢^𝑗(𝑥)−𝑢^𝑗(𝑥)_∞≤11+𝑟Δ𝑡𝑢(𝑥,𝑡_𝑗+1)−𝑢^𝑗+1(𝑥)_∞. (16) Thus, from (14)–(16), we have 𝐸^𝑗_∞≤𝐶𝑒^𝑗_∞+𝑒^𝑗+1_∞+⋯+𝑒^𝐾_∞≤𝐶Δ𝑡, (17) for 0≤𝑗≤𝐾, where 𝐶 is a positive constant independent from Δ𝑡.

4. The Spatial Discretization

For the approximate solution of the semidiscretization problem (6), the spatial discretization is performed on a uniform mesh Ω^𝑁=𝑥_𝑖𝑥=𝑖ℎ,0≤𝑖≤𝑁,ℎ=^max−𝑥_min𝑁, (18) for the computational domain [𝑥_min,𝑥_max]. Thus, at each time point 𝑡_𝑗, we apply a cubic spline scheme on the above uniform mesh Ω^𝑁 to approximate problem (6).

Let 𝑆^𝑗_Δ(𝑥) be the approximate solution of the exact solution 𝑢^𝑗(𝑥) of the boundary value problem (6) at the 𝑗th time level. At each subinterval [𝑥_𝑖,𝑥_𝑖+1], the cubic spline function 𝑆^𝑗_Δ(𝑥) has the following form: 𝑆^𝑗_Δ(𝑥)=𝑎^𝑗_𝑖+𝑏^𝑗_𝑖𝑥−𝑥_𝑖+𝑐^𝑗_𝑖𝑥−𝑥_𝑖²+𝑑^𝑗_𝑖𝑥−𝑥_𝑖³,𝑖=0,1,…,𝑁−1, (19) where 𝑎^𝑗_𝑖,𝑏^𝑗_𝑖,𝑐^𝑗_𝑖, and 𝑑^𝑗_𝑖 are constants. Using the notation 𝑈^𝑗_𝑖 for approximation of 𝑢^𝑗(𝑥) at mesh points 𝑥_𝑖 and 𝑆^𝑗_Δ(𝑥_𝑖)=𝑈^𝑗_𝑖,𝑆^𝑗_Δ(𝑥_𝑖+1)=𝑈^𝑗_𝑖+1 as interpolatory constraints. From algebraic manipulation, we can obtain 𝑎^𝑗_𝑖=𝑈^𝑗_𝑖,𝑏^𝑗_𝑖=𝑈^𝑗_𝑖+1−𝑈^𝑗_𝑖ℎ−ℎ𝑀^𝑗_𝑖+1−𝑀^𝑗_𝑖6,𝑐^𝑗_𝑖=𝑀^𝑗_𝑖2,𝑑^𝑗_𝑖=𝑀^𝑗_𝑖+1−𝑀^𝑗_𝑖,6ℎ (20) where 𝑀^𝑗_𝑖=(𝑆^𝑗_Δ)^′′(𝑥_𝑖). Using the continuity of the first derivative at mesh point 𝑥_𝑖, we get the following equation: 𝑀^𝑗_𝑖+1+4𝑀^𝑗_𝑖+𝑀^𝑗_𝑖−1=6ℎ²𝑈^𝑗_𝑖+1−2𝑈^𝑗_𝑖+𝑈^𝑗_𝑖−1. (21) Substituting 𝑀^𝑗_𝑖=−2𝑢_𝑖^𝑗+1−Δ𝑡2𝑟^𝑗−𝜎^𝑗_𝑖²𝑢^𝑗_𝑥,𝑖+21+Δ𝑡𝑟^𝑗𝑢^𝑗_𝑖Δ𝑡𝜎^𝑗_𝑖² (22) in (21) and using the following approximations for first-order derivative of 𝑢^𝑗𝑢^𝑗_𝑥,𝑖≃𝑈^𝑗_𝑖+1−𝑈^𝑗_𝑖−1,𝑢2ℎ_{𝑗𝑥,𝑖+1}≃3𝑈^𝑗_𝑖+1−4𝑈^𝑗_𝑖+𝑈^𝑗_𝑖−1,𝑢2ℎ_{𝑗𝑥,𝑖−1}≃−𝑈^𝑗_𝑖+1+4𝑈^𝑗_𝑖−3𝑈^𝑗_𝑖−1,2ℎ (23) we get the following spline difference scheme: 𝐿^𝑁𝑈^𝑗_𝑖=2ℎ²𝜎^𝑗_𝑖−1²𝑈^𝑗+1_𝑖−1+8ℎ²𝜎^𝑗_𝑖²𝑈_𝑖^𝑗+1+2ℎ²𝜎^𝑗_𝑖+1²𝑈^𝑗+1_𝑖+1𝑈,𝑖=1,2,…,𝑁−1,_𝑗0=0,𝑈^𝑗_𝑁=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑗𝑟(𝑠)d𝑠}, (24) where 𝐿^𝑁𝑈^𝑗_𝑖=1−6Δ𝑡−2ℎΔ𝑡𝑝^𝑗_𝑖+1+2ℎΔ𝑡𝑝^𝑗_𝑖+32ℎΔ𝑡𝑝^𝑗_𝑖−1+2ℎ²1+Δ𝑡𝑟^𝑗𝜎^𝑗_𝑖−1²𝑈^𝑗_𝑖−1+12Δ𝑡+2ℎΔ𝑡𝑝^𝑗_𝑖+1−2ℎΔ𝑡𝑝^𝑗_𝑖−1+8ℎ²1+Δ𝑡𝑟^𝑗𝜎^𝑗_𝑖²𝑈^𝑗_𝑖+3−6Δ𝑡−2ℎΔ𝑡𝑝^𝑗_𝑖+1−2ℎΔ𝑡𝑝^𝑗_𝑖+12ℎΔ𝑡𝑝^𝑗_𝑖−1+2ℎ²1+Δ𝑡𝑟^𝑗𝜎^𝑗_𝑖+1²𝑈^𝑗_𝑖+1 (25) and 𝑝^𝑗_𝑖=2𝑟^𝑗/(𝜎^𝑗_𝑖)²−1.

It is easy to see that the matrix associated with the discrete operator 𝐿^𝑁 is an M-matrix for sufficiently small ℎ. Hence, the following discrete maximum principle holds true.

Lemma 4 (discrete maximum principle). For sufficiently small ℎ, the operator 𝐿^𝑁 defined by (25) on the uniform mesh Ω^𝑁 satisfies a discrete maximum principle; that is, if 𝑤 is a mesh function that satisfies 𝑤₀≥0,𝑤_𝑁≥0, and 𝐿^𝑁𝑤_𝑖≥0(1≤𝑖<𝑁), then 𝑤_𝑖≥0, for all 𝑖.

From the above lemma, we can conclude that the spatial discretization scheme (24) is maximum-norm stable.

To prove the convergence of the spline difference scheme, we discretize the auxiliary problem (9) and obtain 𝐿^𝑁𝑈^𝑗_𝑖=2ℎ²𝜎^𝑗_𝑖−1²𝑢𝑥_𝑖−1,𝑡_𝑗+1+8ℎ²𝜎^𝑗_𝑖²𝑢𝑥_𝑖,𝑡_𝑗+1+2ℎ²𝜎^𝑗_𝑖+1²𝑢𝑥_𝑖+1,𝑡_𝑗+1𝑈,𝑖=1,2,…,𝑁−1,_𝑗0=0,𝑈^𝑗_𝑁=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑗𝑟(𝑠)d𝑠}. (26)

Lemma 5. Let 𝑢^𝑗(𝑥) be the solution of (9) and {𝑈^𝑗_𝑖} be the solution of (26). Then, we have the following error estimate: 𝑢^𝑗𝑥_𝑖−𝑈^𝑗_𝑖≤𝐶ℎ²Δ𝑡,0≤𝑖≤𝑁. (27)

Proof. We use a Taylor expansion at 𝑥=𝑥_𝑖 to obtain the following local truncation error estimate: 𝐿^𝑁𝑢^𝑗_𝑖−𝑈^𝑗_𝑖=𝐿^𝑁𝑢^𝑗_𝑖−2ℎ²𝜎^𝑗_𝑖−1²𝑢𝑥_𝑖−1,𝑡_𝑗+1−8ℎ²𝜎^𝑗_𝑖²𝑢𝑥_𝑖,𝑡_𝑗+1−2ℎ²𝜎^𝑗_𝑖+1²𝑢𝑥_𝑖+1,𝑡_𝑗+1=𝐿^𝑁𝑢^𝑗_𝑖−2ℎ²𝜎^𝑗_𝑖−1²𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗_𝑖−1−8ℎ²𝜎^𝑗_𝑖²𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗_𝑖−2ℎ²𝜎^𝑗_𝑖+1²𝐼+Δ𝑡𝐿_𝑥𝑢^𝑗_𝑖+1≤13ℎ⁴Δ𝑡𝑝^𝑗_𝑖−1+2𝑝^𝑗_𝑖+𝑝^𝑗_𝑖+1×𝑑³𝑢^𝑗𝑑𝑥³𝜉_𝑖≤𝐶ℎ⁴Δ𝑡. (28) Hence, using the discrete maximum principle (Lemma 4) for the discrete operator 𝐿^𝑁, we have 𝑢^𝑗𝑥_𝑖−𝑈^𝑗_𝑖≤𝐶ℎ²Δ𝑡,0≤𝑖≤𝑁, (29) which completes the proof.

5. The Fully Discrete Scheme

Combining the time semidiscretization scheme (6) with the spatial discretization scheme (24), we can obtain the following fully discretization scheme: 𝑈^𝐾_𝑖=max𝑒^𝑥,𝐿−𝐸,0_𝑁𝑈^𝑗_𝑖=2ℎ²𝜎^𝑗_𝑖−1²𝑈^𝑗+1_𝑖−1+8ℎ²𝜎^𝑗_𝑖²𝑈_𝑖^𝑗+1+2ℎ²𝜎^𝑗_𝑖+1²𝑈^𝑗+1_𝑖+1𝑈,1≤𝑖<𝑁,_𝑗0=0,𝑈^𝑗_𝑁=𝑒^𝑥max−𝐸𝑒^{−∫𝑇𝑡𝑗𝑟(𝑠)d𝑠},for𝑗=𝐾−1,…,1,0, (30) where the discrete operators 𝐿^𝑁 are described in Section 4 and 𝑈^𝑗_𝑖 is the fully discrete approximation to the exact solution of (4) at the mesh point (𝑥_𝑖,𝑡_𝑗).

Now, we can get the main result for our difference scheme.

Theorem 6. Let 𝑢(𝑥,𝑡) be the exact solution of (4) and let 𝑈 be the discrete solution of the fully discrete scheme (30). Then, the global error of our difference scheme satisfies 𝑢𝑥_𝑖,𝑡_𝑗−𝑈^𝑗_𝑖≤𝐶ℎ²+Δ𝑡,0≤𝑖≤𝑁,0≤𝑗≤𝐾, (31) where 𝐶 is a positive constant independent of ℎ and Δ𝑡.

Proof. The global error at the time 𝑡_𝑗 can be decomposed in the form 𝑢𝑥_𝑖,𝑡_𝑗−𝑈^𝑗_𝑖≤𝑢𝑥_𝑖,𝑡_𝑗−𝑢^𝑗𝑥_𝑖+𝑢^𝑗𝑥_𝑖−𝑈^𝑗_𝑖+𝑈^𝑗_𝑖−𝑈^𝑗_𝑖. (32) From Lemmas 2 and 5, we can obtain 𝑢𝑥_𝑖,𝑡_𝑗−𝑈^𝑗_𝑖≤𝐶Δ𝑡ℎ²+Δ𝑡+𝑈^𝑗_𝑖−𝑈^𝑗_𝑖. (33) Further, it is easy to see that 𝑈^𝑗−𝑈^𝑗 can be written as the solution of one step of (30) with zero boundary conditions and 𝑢(𝑥,𝑡_𝑗+1)−𝑈^𝑗+1 as the final value. Applying the stability of the discrete scheme (Lemma 4), we have 𝑈^𝑗_𝑖−𝑈^𝑗_𝑖≤𝐶𝑢𝑥,𝑡_𝑗+1−𝑈^𝑗+1_∞. (34) Then, from (33) and (34), a recurrence relation for the global errors follows, and, from it, the result of Theorem 6 can be obtained immediately.

6. Numerical Experiments

In this section, we present some numerical results to examine the performance and convergence of the cubic spline method. Errors and convergence rates for the numerical scheme are presented for three examples.

Example 1. European call option with parameters: 𝜎=0.4, 𝑟=0.08, 𝑇=1, 𝐸=1, 𝑥_min=−ln(4𝐸), and 𝑥_max=ln(4𝐸): in this case, the analytical solution is 𝑣(𝑆,𝜏)=𝑆𝑁𝑑₁−𝐸𝑒^{−𝑟(𝑇−𝜏)}𝑁𝑑₂, (35) where 𝑁(𝑥)=1√2𝜋^𝑥−∞𝑒^−𝑦2/2𝑑d𝑦,₁()=()+()𝜎𝑆,𝜏ln𝑆/𝐸𝑟+1/2²(𝑇−𝜏)𝜎√,𝑑𝑇−𝜏₂(𝑆,𝜏)=𝑑₁(𝑆,𝜏)−𝜎√𝑇−𝜏. (36)
The maximum error is given in Table 1. The analytical and numerical solution profiles are given in Figure 1.

Table 1 

Numerical results for Example 1.

Figure 1 

Option value at 𝑡=0 for Example 1.

Example 2. European call option with parameters: 𝜎=0.1, 𝑟=0.06, 𝑇=1, 𝐸=1, 𝑥_min=−ln(4𝐸) and 𝑥_max=ln(4𝐸): as in Example 1, in this case, the analytical solution is known.
The maximum error is given in Table 1. The analytical and numerical solution profiles are given in Figure 2.

Figure 2 

Option value at 𝑡=0 for Example 2.

Example 3. European call option with parameters: 𝜎(𝑆,𝜏)=0.15(0.5+2𝜏)((𝑆/100−1.2)²/((𝑆/100)²+1.44)), 𝑟=0.06, 𝑇=1, 𝐸=1, 𝑥_min=−ln(4𝐸), and 𝑥_max=ln(4𝐸): here, the volatility function 𝜎(𝑆,𝜏) is the same as the one given in Toivanen [28] and Kadalbajoo et al. [23, 24].
In this case, the exact solution is not known. We use the approximated solution of 𝑁=2048 and 𝐾=4096 as the exact solution. We present the error estimates for different 𝑁 and 𝐾. Let 𝑈^2048,4096 denote “the exact solution." We measure the accuracy in the discrete maximum norm 𝑒^𝑁,𝐾=max_𝑖,𝑗𝑈^𝑁,𝐾_𝑖𝑗−𝑈^2048,4096_𝑖,𝑗 (37) and the convergence rate 𝑅^𝑁,𝐾=log₂𝑒^𝑁,𝐾𝑒^2𝑁,4𝐾. (38) The error estimates and convergence rates are listed in Table 3. The analytical and numerical solution profiles are given in Figure 3.
From Figures (1)–(3), it is seen that the numerical solutions by our method are nonoscillatory. From Tables 1, 2, and 3, we see that 𝑒^𝑁,𝐾/𝑒^2𝑁,4𝐾 is close to 4, which supports the convergence estimate of Theorem 6. They indicate that the theoretical results are fairly sharp.

Table 2 

Numerical results for Example 2.

Table 3 

Numerical results for Example 3.

Figure 3 

Option value at 𝑡=0 for Example 3.

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 referee for several suggestions for the improvement of this paper. The work was supported by the National Natural Science Foundation (Grant no. 11201430) of China, Ningbo Municipal Natural Science Foundation, Projects in Science and Technique of Ningbo Municipal (Grant no. 2012B82003) of China, and Key Research Center of Philosophy and Social Science of Zhejiang Province—Modern Port Service Industry and Creative Culture Research Center.