Accurate Numerical Method for Pricing Two-Asset American Put Options

Abstract

We develop an accurate finite difference scheme for pricing two-asset American put options. We use the central difference method for space derivatives and the implicit Euler method for the time derivative. Under certain mesh step size limitations, the matrix associated with the discrete operator is an M-matrix, which ensures that the solutions are oscillation-free. We apply the maximum principle to the discrete linear complementarity problem in two mesh sets and derive the error estimates. It is shown that the scheme is second-order convergent with respect to the spatial variables. Numerical results support the theoretical results.

1. Introduction

An option is a financial instrument that gives the holder the right, but not the obligation, to buy (call option) or to sell (put option) an agreed quantity of a specified asset at a fixed price (exercise or strike price) on (European option) or before (American option) a given date (expiry date). It was shown by Black-Scholes [1] that the value of a European option is governed by a second-order parabolic partial differential equation with respect to the time and the underlying asset price. The value of an American option is determined by a linear complementarity problem involving the Black-Scholes operator [2, 3]. Since this complementarity problem is, in general, not analytically solvable, numerical approximation to the solution is normally sought in practice.

Various numerical methods have been proposed for the valuation of single-factor American options. Among them, the lattice method [4], the Monte Carlo method [5], the finite difference method [6–8], the finite element method [9, 10], and the finite volume method [11–13] are the most popular ones in both practice and research.

Finite difference methods applied to the multifactor American option valuation have also been developed. S. O'Sullivan and C. O’Sullivan [14] presented explicit finite difference methods with an acceleration technique for option pricing. Clarke and Parrott [15] and Oosterlee [16] used finite difference schemes along with a projected full approximation scheme (PFAS) multigrid for pricing American options under stochastic volatility. Ikonen and Toivanen [17–19] proposed finite difference methods with componentwise splitting methods on nonuniform grids for pricing American options under stochastic volatility. Hout and Foulon [20] and Zhu and Chen [21] applied finite difference schemes based on the ADI method to price American options under stochastic volatility. Le et al. [22] presented an upwind difference scheme for the valuation of perpetual American put options under stochastic volatility. Yousuf [23] developed an exponential time differencing scheme with a splitting technique for pricing American options under stochastic volatility. Nielsen et al. [24] and Zhang et al. [25] analyzed finite difference schemes with penalty methods for pricing American two-asset options, but their difference methods are first-order convergent.

In part of the domain, the differential operator of the two-asset American option pricing model becomes a convection-dominated operator. The differential operator also contains a second-order mixed derivative term. The classical finite difference methods lead to some off-diagonal elements in the coefficient matrix of the discrete operator due to the dominating first-order derivatives and the mixed derivative. These elements can lead to nonphysical oscillations in the computed solution [17, 18]. In this paper, we present an accurate finite difference scheme for pricing two-asset American options. We use the central difference method for space derivatives and the implicit Euler method for the time derivative. Under certain mesh step size limitations, we obtain a coefficient matrix with an M-matrix property, which ensures that the solutions are oscillation-free. We apply the maximum principle to the discrete linear complementarity problem in two mesh sets and derive the error estimates. We will show that the scheme is second-order convergent with respect to the spatial variables.

The rest of the paper is organized as follows. In the next section, we describe some theoretical results on the continuous complementarity problem for the two-asset American put option pricing model. In Section 3, the discretization method is described. In Section 4, we present a stability and error analysis for the finite difference scheme. In Section 5, numerical experiments are provided to support these theoretical results.

2. The Continuous Problem

We consider the following two-asset American put option pricing model [24, 25]: ℓ𝑃𝑆₁,𝑆₂,𝑡≥0,𝑆₁,𝑆₂[),𝑆>0,𝑡∈0,𝑇ℓ𝑃₁,𝑆₂,𝑡⋅𝑃𝑆₁,𝑆₂,𝑡−𝜙𝑆₁,𝑆₂𝑆=0,₁,𝑆₂[),𝑃𝑆>0,𝑡∈0,𝑇₁,𝑆₂,𝑡−𝜙𝑆₁,𝑆₂≥0,𝑆₁,𝑆₂[],𝑃𝑆≥0,𝑡∈0,𝑇₁,𝑆₂,𝑇=𝜙𝑆₁,𝑆₂,𝑆₁,𝑆₂𝑃≥0,0,𝑆₂,𝑡=𝑔₂𝑆₂,𝑡,𝑆₂[],𝑃𝑆≥0,𝑡∈0,𝑇₁,0,𝑡=𝑔₁𝑆₁,𝑡,𝑆₁≥0,𝑡∈[0,𝑇],lim_𝑆1→∞𝑃𝑆₁,𝑆₂,𝑡=0,𝑆₂≥0,𝑡∈[0,𝑇],lim_𝑆2→∞𝑃𝑆₁,𝑆₂,𝑡=0,𝑆₁≥0,𝑡∈[0,𝑇], (1) where ℓ denotes the two-dimensional Black-Scholes operator defined by ℓ𝑃𝑆₁,𝑆₂,𝑡≡−𝜕𝑃−1𝜕𝑡2𝜎²₁𝑆²₁𝜕²𝑃𝜕𝑆²₁−𝜌𝜎₁𝜎₂𝑆₁𝑆₂𝜕²𝑃𝜕𝑆₁−1𝜕𝑦2𝜎²₂𝑆²₂𝜕²𝑃𝜕𝑆²₂−𝑟𝑆₁𝜕𝑃𝜕𝑆₁−𝑟𝑆₂𝜕𝑃𝜕𝑆₂+𝑟𝑃, (2) and 𝜙(𝑆₁,𝑆₂) is the final (payoff) condition defined by 𝜙𝑆₁,𝑆₂=max𝐸−𝛼₁𝑆₁+𝛼₂𝑆₂,0. (3) Here, 𝑃 is the value of the option, 𝑆_𝑖 is the value of the 𝑖th underlying asset, 𝜌∈[−1,0)∪(0,1] is the correlation of two underlying assets, 𝑟 is the risk-free interest rate, and 𝑔_𝑖(⋅,⋅) is a given function providing suitable boundary conditions. Typically, 𝑔_𝑖(⋅,⋅) is determined by solving the associated one-dimensional American put option problem ℓ_𝑖𝑔_𝑖𝑆_𝑖,𝑡≥0,𝑆_𝑖>0,𝑡∈[0,𝑇),ℓ_𝑖𝑔_𝑖𝑆_𝑖,𝑡⋅𝑔_𝑖𝑆_𝑖,𝑡−max𝐸−𝛼_𝑖𝑆_𝑖𝑆,0=0,_𝑖[),𝑔>0,𝑡∈0,𝑇_𝑖𝑆_𝑖,𝑡−max𝐸−𝛼_𝑖𝑆_𝑖,0≥0,𝑆_𝑖[],𝑔≥0,𝑡∈0,𝑇_𝑖𝑆_𝑖,𝑇=max𝐸−𝛼_𝑖𝑆_𝑖,0,𝑆_𝑖𝑔≥0,_𝑖(0,𝑡)=𝐸,lim_{𝑆𝑖→∞}𝑔_𝑖𝑆_𝑖,𝑡=0,𝑡∈[0,𝑇], (4) where ℓ_𝑖 denotes the one-dimensional Black-Scholes operator defined by ℓ_𝑖𝑔_𝑖𝑆_𝑖,𝑡≡−𝜕𝑔_𝑖−1𝜕𝑡2𝜎²_𝑖𝑆²_𝑖𝜕²𝑔_𝑖𝜕𝑆²_𝑖−𝑟𝑆_𝑖𝜕𝑔_𝑖𝜕𝑆_𝑖+𝑟𝑔_𝑖,𝑖=1,2. (5)

Introducing the logarithmic prices 𝑥=ln𝑆₁ and 𝑦=ln𝑆₂, the linear complementarity problem (1) is transformed as [),⋅𝑢[),𝑢[],𝑢,𝑢𝐿𝑢𝑥,𝑦,𝑡≥0,𝑥,𝑦,𝑡∈ℝ×ℝ×0,𝑇𝐿𝑢𝑥,𝑦,𝑡𝑥,𝑦,𝑡−𝜑𝑥,𝑦=0,𝑥,𝑦,𝑡∈ℝ×ℝ×0,𝑇𝑥,𝑦,𝑡−𝜑𝑥,𝑦≥0,𝑥,𝑦,𝑡∈ℝ×ℝ×0,𝑇𝑥,𝑦,𝑇=𝜑𝑥,𝑦𝑥,𝑦∈ℝ×ℝ,0,𝑦,𝑡=𝑔₂𝑒^𝑦,[],𝑢(),𝑡𝑦,𝑡∈ℝ×0,𝑇𝑥,0,𝑡=𝑔₁𝑒^𝑥,𝑡,(𝑥,𝑡)∈ℝ×[0,𝑇],lim_𝑥→∞𝑢𝑥,𝑦,𝑡=0,𝑦,𝑡∈ℝ×[0,𝑇],lim_𝑦→∞𝑢𝑥,𝑦,𝑡=0,(𝑥,𝑡)∈ℝ×[0,𝑇], (6) where 𝐿𝑢𝑥,𝑦,𝑡≡−𝜕𝑢−1𝜕𝑡2𝜎²₁𝜕²𝑢𝜕𝑥²−𝜌𝜎₁𝜎₂𝜕²𝑢−1𝜕𝑥𝜕𝑦2𝜎²₂𝜕²𝑢𝜕𝑦²−1𝑟−2𝜎²₁𝜕𝑢−1𝜕𝑥𝑟−2𝜎²₂𝜕𝑢𝜑𝛼𝜕𝑦+𝑟𝑢,𝑥,𝑦=max𝐸−₁𝑒^𝑥+𝛼₂𝑒^𝑦,0. (7)

For applying the numerical method, we truncate the infinite domain into Ω≡(𝑥_min,𝑥_max)×(𝑦_min,𝑦_max), where the boundaries 𝑥_min,𝑥_max,𝑦_min, and 𝑦_max are chosen so as not to introduce huge errors in the value of the option [26]. Based on Willmott et al.'s estimate [3] 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𝐸) and 𝑦_max=ln(4𝐸). The artificial boundary conditions at 𝑥=𝑥_min and 𝑥=𝑥_max are chosen to be 𝑢(𝑥_min,𝑦,𝑡)=𝑔₂(𝑒^𝑦,𝑡),𝑢(𝑥_max,𝑦,𝑡)=0. The artificial boundary conditions at 𝑦=𝑦_min and 𝑦=𝑦_max are chosen to be 𝑢(𝑥,𝑦_min,𝑡)=𝑔₁(𝑒^𝑥,𝑡),𝑢(𝑥,𝑦_max,𝑡)=0. Therefore, in the rest of this paper, we will consider the following linear complementary problem: [),⋅𝑢[),𝑢[],𝑢,𝑢𝑥𝐿𝑢𝑥,𝑦,𝑡≥0,𝑥,𝑦,𝑡∈Ω×0,𝑇𝐿𝑢𝑥,𝑦,𝑡𝑥,𝑦,𝑡−𝜑𝑥,𝑦=0,𝑥,𝑦,𝑡∈Ω×0,𝑇𝑥,𝑦,𝑡−𝜑𝑥,𝑦≥0,𝑥,𝑦,𝑡∈Ω×0,𝑇𝑥,𝑦,𝑇=𝜑𝑥,𝑦𝑥,𝑦∈Ω,_min,𝑦,𝑡=𝑔₂𝑒^𝑦,𝑡,𝑦,𝑡∈𝑦_min,𝑦_max×[],𝑢𝑥0,𝑇_max,𝑦,𝑡=0,𝑦,𝑡∈𝑦_min,𝑦_max×[],𝑢0,𝑇𝑥,𝑦_min,𝑡=𝑔₁𝑒^𝑥,𝑡,(𝑥,𝑡)∈𝑥_min,𝑥_max×[],𝑢0,𝑇𝑥,𝑦_max,𝑡=0,(𝑥,𝑡)∈𝑥_min,𝑥_max×[0,𝑇]. (8)

3. Discretization

The operator 𝐿 contains a second-order mixed derivative term. Usual finite difference approximations lead to some positive off-diagonal elements in the matrix associated with the discrete operator due to the mixed derivative, which may lead to nonphysical oscillations in the computed solution. Hence, it is not easy to construct a discretization with good properties and accuracy for problems with mixed derivatives. There are some works dealing with stable difference approximations of mixed derivatives [27, 28]. In this paper, we present an accurate finite difference scheme to discretize the operator 𝐿. We use the technique of [22] to give the mesh step size limitation, which guarantees that the coefficient matrix corresponding to the discrete operator is an 𝑀-matrix.

The discretization is performed using a uniform mesh Ω^{𝑁,𝑀,𝐾} for the computational domain Ω×[0,𝑇]. The mesh steps to the 𝑥 direction, 𝑦 direction, and 𝑡 direction are denoted by Δ𝑥=(𝑥_max−𝑥_min)/𝑁, Δ𝑦=(𝑦_max−𝑦_min)/𝑀, and Δ𝑡=𝑇/𝐾. The mesh point values of the finite difference approximation are denoted by 𝑈^𝑘_𝑖,𝑗≈𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘for𝑖=0,1,…,𝑁;𝑗=0,1,…,𝑀;𝑘=0,1,…,𝐾. (9)

We discretize the differential operator 𝐿 using the central difference scheme on the previous uniform mesh. We set 𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗≡−𝐷⁺_𝑡𝑈^𝑘_𝑖,𝑗−12𝜎²₁𝛿²_𝑥𝑈^𝑘_𝑖,𝑗−𝜎₁𝜎₂𝜌⁺𝛿⁺_𝑥𝑦𝑈^𝑘_𝑖,𝑗+𝜌⁻𝛿⁻_𝑥𝑦𝑈^𝑘_𝑖,𝑗−12𝜎²₂𝛿²_𝑦𝑈^𝑘_𝑖,𝑗−1𝑟−2𝜎²₁𝐷_𝑥𝑈^𝑘_𝑖,𝑗−1𝑟−2𝜎²₂𝐷_𝑦𝑈^𝑘_𝑖,𝑗+𝑟𝑈^𝑘_𝑖,𝑗, (10) where 𝛿²_𝑥𝑈^𝑘_𝑖,𝑗=𝐷⁺_𝑥−𝐷⁻_𝑥𝑈Δ𝑥^{𝑘𝑖,𝑗},𝛿²_𝑦𝑈^𝑘_𝑖,𝑗=𝐷⁺_𝑦−𝐷⁻_𝑦𝑈Δ𝑦^{𝑘𝑖,𝑗},𝛿⁺_𝑥𝑦𝑈^𝑘_𝑖,𝑗=𝐷⁺_𝑥𝐷⁺_𝑦+𝐷⁻_𝑥𝐷⁻_𝑦2𝑈^𝑘_𝑖,𝑗,𝛿⁻_𝑥𝑦𝑈^𝑘_𝑖,𝑗=𝐷⁺_𝑥𝐷⁻_𝑦+𝐷⁻_𝑥𝐷⁺_𝑦2𝑈^𝑘_𝑖,𝑗,𝐷⁺_𝑥𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖+1,𝑗−𝑈^𝑘_𝑖,𝑗Δ𝑥,𝐷⁺_𝑦𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖,𝑗+1−𝑈^𝑘_𝑖,𝑗,𝐷Δ𝑦_−𝑥𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖,𝑗−𝑈^𝑘_{𝑖−1,𝑗}Δ𝑥,𝐷⁻_𝑦𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖,𝑗−𝑈^𝑘_{𝑖,𝑗−1},𝐷Δ𝑦_𝑥𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖+1,𝑗−𝑈^𝑘_{𝑖−1,𝑗}2Δ𝑥,𝐷_𝑦𝑈^𝑘_𝑖,𝑗=𝑈^𝑘_𝑖,𝑗+1−𝑈^𝑘_{𝑖,𝑗−1},𝐷2Δ𝑦_+𝑡𝑈^𝑘_𝑖,𝑗=𝑈^𝑘+1_𝑖,𝑗−𝑈^𝑘_𝑖,𝑗,𝜌Δ𝑡^±=12𝜌±𝜌. (11)

Denote Ω_ℎ=,Ω𝑖,𝑗,𝑘∣0≤𝑖≤𝑁,0≤𝑗≤𝑀,0≤𝑘≤𝐾_ℎ=},𝑖,𝑗,𝑘∣1≤𝑖≤𝑁−1,1≤𝑗≤𝑀−1,1≤𝑘≤𝐾−1𝜕Ω_ℎ=Ω_ℎ⧵Ω_ℎ. (12)

Thus, we apply the central difference scheme on the uniform mesh to approximate the parabolic complementarity problem (8) as follows: 𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗≥0,𝑖,𝑗,𝑘∈Ω_ℎ,𝑈^𝑘_𝑖,𝑗−𝜑_𝑖,𝑗≥0,𝑖,𝑗,𝑘∈Ω_ℎ,𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗⋅𝑈^𝑘_𝑖,𝑗−𝜑_𝑖,𝑗=0,𝑖,𝑗,𝑘∈Ω_ℎ,𝑈^𝐾_𝑖,𝑗=𝜑_𝑖,𝑗𝑈,0≤𝑖≤𝑁,0≤𝑗≤𝑀,_𝑘0,𝑗=𝑔₂^𝑘_𝑗,𝑈^𝑘_𝑁,𝑗𝑈=0,0≤𝑗≤𝑀,0≤𝑘<𝐾,_𝑘𝑖,0=𝑔₁^𝑘_𝑖,𝑈^𝑘_𝑖,𝑀=0,0<𝑖<𝑁,0≤𝑘<𝐾. (13) Here, (𝑔₁)^𝑘_𝑖and(𝑔₂)^𝑘_𝑗 are discrete approximates of 𝑔₁(𝑒^𝑥,𝑡)and𝑔₂(𝑒^𝑦,𝑡), respectively. Hence, (𝑔₁)^𝑘_𝑖 and (𝑔₂)^𝑘_𝑗 can be obtained by solving the corresponding one-dimensional Black-Scholes equations [29]. In the next section, we will prove that the system matrix corresponding to the discrete operator 𝐿^{𝑁,𝑀,𝐾} is an M-matrix. Hence, from the uniqueness theorem of Goeleven [30], we can obtain that there exists a unique solution 𝑈 for the previous linear complementarity problem (13).

4. Analysis of the Method

First, we give the stability analysis for the difference scheme (13).

Lemma 1. If mesh steps satisfy the inequalities 𝜎Δ𝑥≤²¹2𝑟−𝜎²₁𝜎,Δ𝑦≤²²2𝑟−𝜎²₂,2𝜌𝜎₁𝜎₂≤Δ𝑥≤Δ𝑦2𝜎₁𝜌𝜎₂, (14)(15) then the system matrix corresponding to the discrete operator 𝐿^{𝑁,𝑀,𝐾} is an 𝑀-matrix.

Proof. The difference operator 𝐿^{𝑁,𝑀,𝐾} can be written as follows: 𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗𝜌=−⁺𝜎₁𝜎₂𝑈2Δ𝑥Δ𝑦^{𝑘𝑖−1,𝑗−1}+𝜎₁𝜎₂𝜌⁺−𝜌⁻−𝜎2Δ𝑥Δ𝑦2Δ𝑦²²²+𝑟−(1/2)𝜎²₂𝑈2Δ𝑦^{𝑘𝑖,𝑗−1}+𝜌⁻𝜎₁𝜎₂𝑈2Δ𝑥Δ𝑦^{𝑘𝑖+1,𝑗−1}+𝜎₁𝜎₂𝜌⁺−𝜌⁻−𝜎2Δ𝑥Δ𝑦2(Δ𝑥²¹)²+𝑟−(1/2)𝜎²₁𝑈2Δ𝑥^{𝑘𝑖−1,𝑗}+1+𝜎Δ𝑡(Δ𝑥²¹)²−𝜎₁𝜎₂𝜌⁺−𝜌⁻+𝜎Δ𝑥Δ𝑦Δ𝑦²²²+𝑟𝑈^𝑘_𝑖,𝑗+−𝜎²₁2(Δ𝑥)²+𝜎₁𝜎₂𝜌⁺−𝜌⁻−()𝜎2Δ𝑥Δ𝑦𝑟−1/2²₁𝑈2Δ𝑥^{𝑘𝑖+1,𝑗}+𝜌⁻𝜎₁𝜎₂𝑈2Δ𝑥Δ𝑦^{𝑘𝑖−1,𝑗+1}+𝜎₁𝜎₂𝜌⁺−𝜌⁻−𝜎2Δ𝑥Δ𝑦2Δ𝑦²²²−𝑟−(1/2)𝜎²₂𝑈2Δ𝑦^{𝑘𝑖,𝑗+1}−𝜌⁺𝜎₁𝜎₂𝑈2Δ𝑥Δ𝑦^{𝑘𝑖+1,𝑗+1}−1𝑈Δ𝑡^{𝑘+1𝑖,𝑗}. (16) The coefficient of 𝑈_𝑖,𝑗 in the previous expression (which corresponds to the diagonal of the system matrix) is positive since 𝜎²₁(Δ𝑥)²−𝜎₁𝜎₂𝜌⁺−𝜌⁻+𝜎Δ𝑥Δ𝑦Δ𝑦²²²≥0. (17) All the coefficients of the other 𝑈 in the previous expression (which correspond to off-diagonal elements in the system matrix) will be nonpositive once the following inequalities are satisfied: 𝜎²₁4(Δ𝑥)²−𝑟−(1/2)𝜎²₁𝜎2Δ𝑥≥0,²²4Δ𝑦²−𝑟−(1/2)𝜎²₂𝜎2Δ𝑦≥0,₁𝜎₂𝜌⁺−𝜌⁻−𝜎2Δ𝑥Δ𝑦4Δ𝑦²²²≤0,𝜎₁𝜎₂𝜌⁺−𝜌⁻−𝜎2Δ𝑥Δ𝑦4(Δ𝑥²¹)²≤0. (18) Together, they require that the following inequalities hold: 𝜎Δ𝑥≤²¹2𝑟−𝜎²₁𝜎,Δ𝑦≤²²2𝑟−𝜎²₂,2𝜎₁𝜌⁺−𝜌⁻𝜎₂≤Δ𝑥≤Δ𝑦2𝜎₁𝜎₂𝜌⁺−𝜌⁻, (19) which are (14) and (15), respectively. Thus, we have shown that the system matrix, corresponding to the discrete operator 𝐿^{𝑁,𝑀,𝐾} is an 𝑀-matrix and the result follows.

There are only few error estimates for the direct application of finite difference method to linear complementarity problems. Here, we apply the maximum principle to the linear complementarity problem (13) in two mesh sets and derive the error estimates [29, 31].

By using Taylor's formula, we can easily obtain the following truncation error estimate.

Lemma 2. Let 𝑢(𝑥,𝑦,𝑡) be a smooth function defined on Ω^{𝑁,𝑀,𝐾}. Then the truncation error of the difference scheme (10) satisfies 𝐿^{𝑁,𝑀,𝐾}𝑢^𝑘_𝑖,𝑗−𝐿𝑢^𝑘_𝑖,𝑗=𝑂(Δ𝑥)²+Δ𝑦²+Δ𝑥Δ𝑦+Δ𝑡, (20) for all (𝑖,𝑗,𝑘)∈Ω_ℎ.

Now we can derive our main result for the difference scheme.

Theorem 3. Let 𝑢(𝑥,𝑦,𝑡) be the solution of the problem (8) and let 𝑈^𝑘_𝑖,𝑗 be the solution of the problem (13). If mesh steps satisfy conditions (14) and (15), the difference scheme (13) satisfies the following error estimate: max_{(𝑖,𝑗,𝑘)∈Ωℎ}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗≤𝐶(Δ𝑥)²+Δ𝑦²+Δ𝑥Δ𝑦+Δ𝑡, (21) where 𝐶 is a constant independent of Δ𝑥,Δ𝑦, and Δ𝑡.

Proof. Denote Ω⁽¹⁾=𝑖,𝑗,𝑘∈Ω_ℎ∣𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘=𝜑𝑥_𝑖,𝑦_𝑗,Ω⁽²⁾=Ω_ℎ⧵Ω⁽¹⁾. (22) From (8), we have the result 𝐿𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘≥0,𝑖,𝑗,𝑘∈Ω⁽¹⁾,𝐿𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘=0,𝑖,𝑗,𝑘∈Ω⁽²⁾. (23) Denote Ω_ℎ⁽¹⁾=𝑖,𝑗∈Ω_ℎ∣𝑈^𝑘_𝑖,𝑗=𝜑𝑥_𝑖,𝑦_𝑗,Ω_ℎ⁽²⁾=Ω_ℎ⧵Ω_ℎ⁽¹⁾. (24) Obviously, 𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗=0,𝑖,𝑗,𝑘∈Ω_ℎ⁽²⁾. (25)
Define the function on Ω_ℎ by 𝑊^𝑘_𝑖,𝑗=𝐶(Δ𝑥)²+Δ𝑦²+Δ𝑥Δ𝑦+Δ𝑡>0, (26) where 𝐶 is a sufficiently large constant.
For (𝑖,𝑗,𝑘)∈Ω_ℎ⁽²⁾, by the fact that 𝐿𝑢(𝑥_𝑖,𝑦_𝑗,𝑡_𝑘)≥0, (25), (26), and Lemma 2, we obtain 𝐿^{𝑁,𝑀,𝐾}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗+𝑊^𝑘_𝑖,𝑗=𝐿^{𝑁,𝑀,𝐾}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘+𝐿^{𝑁,𝑀,𝐾}𝑊^𝑘_𝑖,𝑗=𝐿^{𝑁,𝑀,𝐾}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝐿𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘+𝐿^{𝑁,𝑀,𝐾}𝑊^𝑘_𝑖,𝑗+𝐿𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘≥0. (27) On the “boundary” of Ω_ℎ⁽²⁾, the nodes (𝑖,𝑗,𝑘)∈Ω_ℎ⁽¹⁾, so 𝑈^𝑘_𝑖,𝑗=𝜑(𝑥_𝑖,𝑦_𝑗), but 𝑢(𝑥_𝑖,𝑦_𝑗,𝑡_𝑘)≥𝜑(𝑥_𝑖,𝑦_𝑗), therefore 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗+𝑊^𝑘_𝑖,𝑗=𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝜑𝑥_𝑖,𝑦_𝑗+𝑊^𝑘_𝑖,𝑗≥0, (28) and the nodes (𝑖,𝑗,𝑘)∈𝜕Ω_ℎ, 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗+𝑊^𝑘_𝑖,𝑗=𝑊^𝑘_𝑖,𝑗≥0. (29) Applying the maximum principle to Ω_ℎ⁽²⁾, we get 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗+𝑊^𝑘_𝑖,𝑗≥0,𝑖,𝑗,𝑘∈Ω_ℎ⁽²⁾. (30) Thus, 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗+𝑊^𝑘_𝑖,𝑗≥0,𝑖,𝑗,𝑘∈Ω_ℎ. (31)
For (𝑖,𝑗,𝑘)∈Ω⁽²⁾,𝐿𝑢(𝑥_𝑖,𝑦_𝑗,𝑡_𝑘)=0, but 𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗≥0, thus, 𝐿^{𝑁,𝑀,𝐾}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗=𝐿^{𝑁,𝑀,𝐾}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝐿𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝐿^{𝑁,𝑀,𝐾}𝑊^𝑘_𝑖,𝑗−𝐿^{𝑁,𝑀,𝐾}𝑈^𝑘_𝑖,𝑗≤0. (32) On the “boundary” of Ω⁽²⁾, the nodes (𝑖,𝑗,𝑘)∈Ω⁽¹⁾, so 𝑢(𝑥_𝑖,𝑦_𝑗,𝑡_𝑘)=𝜑(𝑥_𝑖,𝑦_𝑗), but 𝑈^𝑘_𝑖,𝑗≥𝜑(𝑥_𝑖,𝑦_𝑗), therefore 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗=𝜑𝑥_𝑖,𝑦_𝑗−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗≤0, (33) and the nodes (𝑖,𝑗,𝑘)∈𝜕Ω_ℎ, 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗=−𝑊^𝑘_𝑖,𝑗≤0. (34) Applying the maximum principle to Ω⁽²⁾, we get 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗≤0,𝑖,𝑗,𝑘∈Ω⁽²⁾. (35) Thus, 𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗−𝑊^𝑘_𝑖,𝑗≤0,𝑖,𝑗,𝑘∈Ω_ℎ. (36)
From (31) and (36), we obtain max_{∈𝑖,𝑗,𝑘Ωℎ}𝑢𝑥_𝑖,𝑦_𝑗,𝑡_𝑘−𝑈^𝑘_𝑖,𝑗≤max_{∈𝑖,𝑗,𝑘Ωℎ}𝑊^𝑘_𝑖,𝑗≤𝐶(Δ𝑥)²+Δ𝑦²+Δ𝑥Δ𝑦+Δ𝑡, (37) where 𝐶 is a sufficiently large constant. From this we complete the proof.

5. Numerical Experiments

In this section, we verify experimentally the theoretical results obtained in the preceding section. Errors and convergence rates for the second-order finite difference scheme are presented for two test problems.

Test  1. American put option with parameters: 𝑇=1, 𝑟=0.1, 𝜎₁=0.4, 𝜎₂=0.5, 𝜌=0.5, 𝑥_min=−ln(30), 𝑥_max=ln(30), 𝑦_min=−ln(40), 𝑦_max=ln(40), 𝛼₁=0.3, 𝛼₂=0.7, and 𝐸=10.

Test  2. American put option with parameters: 𝑇=1, 𝑟=0.08, 𝜎₁=0.3, 𝜎₂=0.4, 𝜌=−0.6, 𝑥_min=−ln(30), 𝑥_max=ln(30), 𝑦_min=−ln(40), 𝑦_max=ln(40), 𝛼₁=0.3, 𝛼₂=0.7, and 𝐸=10.

To solve the linear inequality system (13), we use the projection scheme used in [32, page 433]. Since mesh steps need to satisfy conditions (14) and (15), we choose the number of mesh steps in the 𝑦 direction 2𝜌𝜎𝑀=¹/𝜎₂+2𝜎₁/𝜌𝜎₂𝑁2𝑥_max−𝑥_min, (38) where 𝑁 is the number of mesh steps in the 𝑥 direction. The exact solutions of the test problems are not available. Therefore, we use the double mesh principle to estimate the errors and compute the experiment convergence rates in our computed solution. We measure the accuracy in the discrete maximum norm 𝑒^{𝑁,𝑀,𝐾}=max_{𝑖,𝑗,𝑘}𝑈^{𝑁,𝑀,𝐾}_{𝑖,𝑗,𝑘}−𝑈^{2𝑁,2𝑀,𝐾}_{𝑖,𝑗,𝑘}, (39) and the convergence rate 𝑅^{𝑁,𝑀,𝐾}=log₂𝑒^{𝑁,𝑀,𝐾}𝑒^{2𝑁,2𝑀,𝐾}. (40)

The error estimates and convergence rates in our computed solutions of Tests 1 and 2 are listed in Tables 1 and 2, respectively. From Tables 1 and 2, we see that 𝑒^{𝑁,𝑀,𝐾}/𝑒^{2𝑁,2𝑀,𝐾} is close to 4 for sufficiently large 𝐾, which supports the convergence estimate of Theorem 3. However, the numerical results of Nielsen et al. [24] and Zhang et al. [25] verify that their schemes are only first-order convergent. Hence, our scheme is more accurate.

Table 1 

Numerical results for Test 1.

Table 2 

Numerical results for Test 2.

Acknowledgments

The authors would like to thank the anonymous referees for several suggestions for the improvement of this paper. The work was supported by Zhejiang Province Natural Science Foundation of China (Grant no. Y2111160).