An Alternating-Direction Implicit Difference Scheme for Pricing Asian Options

Abstract

We propose a fast and stable numerical method to evaluate two-dimensional partial differential equation (PDE) for pricing arithmetic average Asian options. The numerical method is deduced by combining an alternating-direction technique and the central difference scheme on a piecewise uniform mesh. The numerical scheme is stable in the maximum norm, which is true for arbitrary volatility and arbitrary interest rate. It is proved that the scheme is second-order convergent with respect to the asset price. Numerical results support the theoretical results.

1. Introduction

An Asian option is a derivative product whose payoff depends on the average price of the underlying asset, which can be stocks, commodities, or financial indices. A company that has to purchase huge amount of an asset such as oil at a fixed date, but has to sell it by small amount during some period, could think of Asian options as an effective way to hedge the risk that comes from mismatch of cash flows. The price of the Asian option is less subject to price manipulation. Hence, the averaging feature is popular in many thinly traded markets and embedded in complex derivatives such as “refix” clauses in convertible bonds.

Since no general analytical solution for the price of arithmetic Asian option is known, a variety of techniques have been proposed. Binomial lattice methods require such enormous amounts of computer memory (owing to the necessity of keeping track of every possible path throughout the tree); thus, they are effectively unusable. Simulation methods such as Monte Carlo methods have difficulty in achieving high accuracy for Asian options. Other methods include sharp bounds [1], closed-form approximations [2], and Laplace transforms of Asian option prices [3]. The usual numerical partial differential equation (PDE) methods are inaccurate since the degeneration of the PDE for pricing Asian options causes numerical diffusion and spurious oscillation [4, 5]. Zvan et al. [5] apply a flux-limiting method from computational fluid dynamics to tackle the problem of spurious oscillations that arise in Asian options. Hugger [6] proposes an artificial viscosity numerical method for Asian options to avoid the oscillations. Večeř [7] presents that the pricing techniques of an option on a traded account could be applied to price the Asian option, and that the price of the Asian option is characterized by a simple one-dimensional PDE. Dubois and Lelièvre [8] use a characteristic method to solve the Večeř PDE for the Asian option. Cen et al. [9] present a robust finite difference scheme with a moving mesh for the Večeř PDE. Mudzimbabwe et al. [10] propose an implicit finite difference method for the one-dimensional PDE which is obtained by the exponential transformed for pricing Asian options. Tangman et al. [11] use the exponential time integration scheme in combination with a dimensional splitting technique for pricing Asian options under a variety of pricing models.

In general, the Asian options pricing model is a two-dimensional PDE. European-style Asian options and American-style floating strike options may be valued using one-dimensional PDEs by making the change of variables. But, for the American-style fixed strike options, a two-dimensional PDE needs to be solved. For the sake of generality, in this paper we consider the two-dimensional PDE for pricing the Asian options.

Note that the two-dimensional PDE leads to greater computational costs. It is very natural to consider dimension splitting: at each time-step, one alternately solves independent one-dimensional problems. Such alternating-direction implicit (ADI) scheme for classical problems is presented in detail by, for example, Marchuk [12], Samarskiĭ [13], and Strikwerda [14]. Here, we take advantage of the computational cost reduction yielded by the use of alternating directions and of the robust difference schemes for both the one-dimensional Black-Scholes equation and the advection equation.

When one uses the standard finite difference method to solve the Black-Scholes equation, numerical difficulty rises, especially when the volatility 𝜎 is small. The main reason is that when the volatility 𝜎 or spatial variable 𝑥 is small, the partial differential operator becomes a convection-dominated operator. Hence, the implicit Euler scheme with central spatial difference method may lead to nonphysical oscillations in the computed solution. The implicit Euler scheme with upwind spatial difference method does not have this disadvantage, but this difference scheme is only first-order convergent. Applying an Euler transformation [15] to remove the singularity of the differential operator when the parameters of the Black-Scholes equation are constant or space independent, the truncation on the left-hand side of the transformed domain to artificially remove the degeneracy may cause computational errors. Furthermore, the uniform mesh on the transformed interval will lead to the originally grid points concentrating around 𝑥=0 inappropriately. We have proposed robust difference schemes for the Black-Scholes equation, which is based on a central difference spatial discretization on a piecewise uniform mesh and an implicit time stepping technique [16, 17].

The numerical method is deduced by combining an alternating-direction technique and the central difference scheme on a piecewise uniform mesh. The scheme is stable for arbitrary volatility and arbitrary interest rate. It is proved that the scheme is second-order convergent with respect to the asset price. Numerical results support the theoretical results.

In this paper, we propose a fast and stable numerical scheme to evaluate two-dimensional partial differential equation arising in pricing arithmetic average Asian options. To avoid solving a large discrete linear system, we apply the ADI scheme to alternately solve a one-dimensional Black-Scholes equation and an advection equation. We apply the central difference scheme on a piecewise uniform mesh to discrete the one-dimensional Black-Scholes equation and apply the implicit Euler scheme to discrete the advection equation. We will show that the matrices associated with discrete operators are M-matrices, which ensures that the numerical scheme is stable in the maximum norm and free of nonphysical oscillations whether 𝜎²/𝑟≤1 or 𝜎²/𝑟>1. It is proved that the scheme is second-order convergent with respect to the asset price. Numerical results support the theoretical results.

The rest of the paper is organized as follows. In the next section, we describe some theoretical results on the Asian option pricing model. In Section 3, we give the corresponding stability and convergence property of the time semidiscretization. The robust difference schemes for the spatial discretization are presented 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.

Notation. We always use the (pointwise) maximum norm ‖⋅‖_𝐷, where 𝐷 is a closed and bounded set.

2. The Continuous Problem

Suppose that the underlying asset price 𝑥(𝑡) follows a geometric Brownian motion 𝑑𝑥(𝑡)=𝑟𝑥(𝑡)𝑑𝑡+𝜎𝑥(𝑡)𝑑𝐵(𝑡), (1) where 𝑟 is risk-free interest rate, 𝜎 is the volatility, and 𝐵(𝑡) is a standard Brownian motion under the risk-neutral measure ℙ. Let 𝐴(𝑡) denote the underlying asset price running sum given by 𝐴(𝑡)=^𝑡₀𝑥(𝜏)d𝜏; (2) then the arithmetically averaged price is given by 𝐴(𝑡)/𝑡. The stochastic differential equation for evolution of (𝐴(𝑡)) is given by 𝑑𝐴(𝑡)=𝑥(𝑡)𝑑𝑡. (3) Thus, the Asian call option price 𝑣(𝑥,𝐴,𝑡) for continuous arithmetic average strike satisfies the following PDE [15, 18]: 𝐿𝑣(𝑥,𝐴,𝑡)≡−𝜕𝑣−1𝜕𝑡2𝜎²𝑥²𝜕²𝑣𝜕𝑥²−𝑟𝑥𝜕𝑣𝜕𝑥−𝑥𝜕𝑣𝜕𝐴+𝑟𝑣=0, (4) with the terminal condition 𝑣()𝐴𝑥,𝐴,𝑇=max𝑇−𝐸,0, (5) and the left-hand boundary condition 𝑣(0,𝐴,𝑡)=𝑒^{−𝑟(𝑇−𝑡)}𝐴max𝑇−𝐸,0, (6) where 𝑇 is the expiry date and 𝐸 is the strike price. Note that in the work of Geman and Yor [3], a simple formula is obtained for 𝐴≥𝐸𝑇: 𝑣()=𝑥𝑥,𝐴,𝑡𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡)}+𝐴𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡)}. (7) From this formula the right-hand boundary condition can be obtained: 𝑣()∼𝑥𝑥,𝐴,𝑡𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡)}+𝑒^{−𝑟(𝑇−𝑡)}𝐴𝑇−𝐸,a𝑠𝑥⟶+∞. (8) See Hugger [19] for the derivation details of boundary value conditions for Asian options.

As the exact solution to the problem (4)–(6) for 𝐴≥𝐸𝑇 is known, we only consider the solution of the PDE for 0<𝐴<𝐸𝑇 using (7) for the boundary condition at 𝐴=𝐸𝑇. For applying the numerical method, we truncate the domain (0,+∞) of spatial variable 𝑥 into (0,𝑋) for sufficiently large 𝑋. Based on Willmott et al.’s estimate [15] that the upper bound of the asset price is typically three or four times the strike price, it is reasonable for us to set 𝑋=4𝐸. The boundary condition at 𝑥=𝑋 is chosen to be 𝑣(𝑋,𝐴,𝑡)=(𝑋/𝑟𝑇)(1−𝑒^{−𝑟(𝑇−𝑡)})+(𝐴/𝑇−𝐸)𝑒^{−𝑟(𝑇−𝑡)}. Normally, this truncation of the domain leads to a negligible error in the value of the option [20]. Therefore, in the remaining of this paper we will consider the following PDE: ()()𝑣()𝐴𝐿𝑣𝑥,𝐴,𝑡=0,𝑥,𝐴,𝑡∈Ω,𝑥,𝐴,𝑇=max𝑇−𝐸,0,(𝑥,𝐴)∈Ω₁×Ω₂,𝑣(0,𝐴,𝑡)=𝑒^{−𝑟(𝑇−𝑡)}𝐴max𝑇−𝐸,0,(𝐴,𝑡)∈Ω₂×Ω₃,𝑣()=𝑋𝑋,𝐴,𝑡𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡)}+𝐴𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡)},(𝐴,𝑡)∈Ω₂×Ω₃,𝑣()=𝑥𝑥,𝐸𝑇,𝑡𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡)},(𝑥,𝑡)∈Ω₁×Ω₃, (9) where Ω₁=(0,𝑋),Ω₂=(0,𝐸𝑇),Ω₃=(0,𝑇) and Ω=Ω₁×Ω₂×Ω₃.

3. The Time Semidiscretization

Note that the above problem is a two-dimensional PDE which leads to greater computational costs. It is very natural to consider dimension splitting: at each time-step, one alternately solves independent one-dimensional problems in the 𝑥 and 𝐴 directions.

We discrete the above PDE in [0,𝑇] using a uniform mesh, which is defined by Ω^𝐾=𝑡_𝑛=𝑛Δ𝑡,0≤𝑛≤𝐾,Δ𝑡=𝑇/𝐾. (10) Split the spatial differential operator into the following two operators: 𝐿_𝑥1≡−2𝜎²𝑥²𝜕²𝜕𝑥²𝜕−𝑟𝑥𝜕𝑥+𝑟𝑣,𝐿_𝐴𝜕=−𝑥.𝜕𝐴 (11) Using the partitions, the semidiscrete scheme can be written as follows:(a)𝑣^𝐾𝐴=max𝑇−𝐸,0, (12)(b)𝐼+Δ𝑡𝐿_𝑥𝑣^𝑛+1/2=𝑣^𝑛+1,𝑣^𝑛+1/2(0,𝐴)=𝑒^{−𝑟(𝑇−𝑡𝑛)}𝐴max𝑇,𝑣−𝐸,0_𝑛+1/2()=𝑋𝑋,𝐴𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}+𝐴𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡𝑛)}, (13)(c)𝐼+Δ𝑡𝐿_𝐴𝑣^𝑛=𝑣^𝑛+1/2,𝑣^𝑛()=𝑥𝑥,𝐸𝑇𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}. (14)This method gives approximations 𝑣^𝑛(𝑥,𝐴) to the solution 𝑣(𝑥,𝐴,𝑡) of (9) at the time levels 𝑡_𝑛. The operators (𝐼+Δ𝑡𝐿_𝑥) and (𝐼+Δ𝑡𝐿_𝐴) satisfy a maximum principle, and consequently 𝐼+Δ𝑡𝐿_𝑥⁻¹_Ω1×Ω2≤1,1+𝑟Δ𝑡𝐼+Δ𝑡𝐿_𝐴⁻¹_Ω1×Ω2≤1. (15) This ensures the stability of the temporal semidiscretization (12)–(14) and that each step of the scheme (12)–(14) has a unique solution 𝑣^𝑛(𝑥,𝐴).

The local truncation error is defined as 𝑒_𝑛=𝑣(𝑡_𝑛)−𝑣^𝑛, where 𝑣^𝑛 is the solution of 𝐼+Δ𝑡𝐿_𝑥𝑣^𝑛+1/2=𝑣𝑥,𝐴,𝑡_𝑛+1,𝑣^𝑛+1/2(0,𝐴)=𝑒^{−𝑟(𝑇−𝑡𝑛)}𝐴max𝑇−𝐸,0,𝑣^𝑛+1/2()=𝑋𝑋,𝐴𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}+𝐴𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡𝑛)},𝐼+Δ𝑡𝐿_𝐴𝑣^𝑛=𝑣^𝑛+1/2,𝑣^𝑛()=𝑥𝑥,𝐸𝑇𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}. (16)(17)(18)(19)(20) The following lemma gives the local error estimates.

Lemma 1. The local error for the scheme (16)–(20) satisfies 𝑒_𝑛Ω1×Ω2≤𝐶(Δ𝑡)². (21)

Proof. The function 𝑣^𝑛 satisfies 𝐼+Δ𝑡𝐿_𝑥𝐼+Δ𝑡𝐿_𝐴𝑣^𝑛(𝑥,𝐴)=𝑣𝑥,𝐴,𝑡_𝑛+1. (22) On the other hand, since the solution of (4) is smooth enough, we have 𝑣𝑥,𝐴,𝑡_𝑛+1=𝑣𝑥,𝐴,𝑡_𝑛+Δ𝑡𝐿_𝑥+𝐿_𝐴𝑣𝑥,𝐴,𝑡_𝑛+^𝑡𝑛+1_𝑡𝑛𝑡_𝑛+1𝜕−𝑠²𝑣𝜕𝑠²=d𝑠𝐼+Δ𝑡𝐿_𝑥𝐼+Δ𝑡𝐿_𝐴𝑣𝑥,𝐴,𝑡_𝑛+𝑂(Δ𝑡)². (23) Hence, 𝑒_𝑛 satisfies 𝐼+Δ𝑡𝐿_𝑥𝐼+Δ𝑡𝐿_𝐴𝑒_𝑛=𝑂(Δ𝑡)²,𝑒_𝑛(0,𝐴)=𝑒_𝑛(𝑋,𝐴)=𝑒_𝑛(𝑥,𝐸𝑇)=0. (24) The application of the stability results (15) proves (21).

The following theorem gives the convergence result of the method (12)–(14).

Theorem 2. The global error associated to the method (12)–(14), defined by 𝐸_𝑛=𝑣(𝑡_𝑛)−𝑣^𝑛, satisfies sup_{0≤𝑛<𝐾}𝐸_𝑛Ω1×Ω2≤𝐶Δ𝑡, (25) and therefore the time semidiscretization method is a first-order convergent scheme.

Proof. It is easy to see that 𝐸_𝑛=𝑒_𝑛+𝑅𝐸_𝑛+1, (26) where 𝑅=𝐼+Δ𝑡𝐿_𝑥⁻¹𝐼+Δ𝑡𝐿_𝐴⁻¹, (27) is the transition operator and 𝑅𝐸_𝑛+1 is the result obtained after one step of scheme (13)-(14) with 𝐸_𝑛+1 as final value. Using this recurrence, we deduce that 𝐸_𝑛=^𝐾−1_𝑖=𝑛𝑅^𝑖−𝑛𝑒_𝑖. (28) From ‖𝑅‖_Ω1×Ω2=𝐼+Δ𝑡𝐿_𝑥⁻¹_Ω1×Ω2⋅𝐼+Δ𝑡𝐿_𝐴⁻¹_Ω1×Ω2≤1,1+𝑟Δ𝑡 (29) we obtain the desired result.

4. The Spatial Discretization

As discussed in Section 1, the standard finite difference method to solve the Black-Scholes equation may cause numerical difficulty. Here, we apply the central difference scheme on a piecewise uniform mesh [16, 17] to discrete problem (13) and apply the implicit Euler scheme to discrete problem (14). The matrices associated with discrete operators are M-matrices, which ensure that the scheme is stable for arbitrary volatility and arbitrary interest rate without any extra conditions.

The use of central difference scheme on a uniform mesh may produce nonphysical oscillations in the computed solution. To overcome this oscillation, we use a piecewise uniform mesh Ω^𝑁 on the space interval [0,𝑋]: 𝑥_𝑖=ℎ𝜎ℎ𝑖=1,1+²𝑟(𝑖−1)𝑖=2,…,𝑁, (30) where 𝑋ℎ=1+𝜎²()./𝑟𝑁−1 (31) For the 𝐴 variable discretization, we use a uniform mesh Ω^𝑀 on [0,𝐸𝑇] with 𝑀 mesh elements. It is easy to see that the mesh sizes ℎ_𝑥,𝑖=𝑥_𝑖−𝑥_𝑖−1 and ℎ_𝐴,𝑗=𝐴_𝑗−𝐴_𝑗−1 satisfy ℎ_𝑥,𝑖=𝜎ℎfor𝑖=1,₂𝑟ℎℎfor𝑖=2,…,𝑁,_𝐴,𝑗=𝐸𝑇𝑀,𝑗=1,…,𝑀, (32)(33) respectively. We denote by Ω^{𝑁,𝑀,𝐾}=Ω^𝑁×Ω^𝑀×Ω^𝐾 the corresponding mesh for 𝑥,𝐴, and 𝑡.

Thus, we apply the central difference scheme on the piecewise uniform mesh and the implicit Euler scheme to approximate problem (16)–(20): 𝐼+Δ𝑡𝐿^𝑁_𝑥𝑉^𝑛+1/2_𝑖𝑗=𝑣^𝑛+1_𝑖𝑗,1≤𝑖≤𝑁,0≤𝑗<𝑀,𝑉^𝑛+1/2_0,𝑗=𝑒^{−𝑟(𝑇−𝑡𝑛)}𝐴max^𝑗𝑇−𝐸,0,1≤𝑖<𝑁,𝑉^𝑛+1/2_𝑁,𝑗=𝑋𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}+𝐴_𝑗𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖<𝑁,𝐼+Δ𝑡𝐿^𝑀_𝐴𝑉^𝑛_𝑖𝑗=𝑉^𝑛+1/2_𝑖𝑗,1≤𝑖≤𝑁,0≤𝑗<𝑀,𝑉^𝑛_𝑖,𝑀=𝑥_𝑖𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖≤𝑁, (34)(35)(36)(37)(38) where 𝐿^𝑁_𝑥𝑉^𝑛+1/2_𝑖𝑗𝜎≡−²𝑥²_𝑖ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1×𝑉^𝑛+1/2_𝑖+1,𝑗−𝑉^𝑛+1/2_𝑖𝑗ℎ_𝑥,𝑖+1−𝑉^𝑛+1/2_𝑖𝑗−𝑉^𝑛+1/2_{𝑖−1,𝑗}ℎ_𝑥,𝑖−𝑟𝑥_𝑖𝑉^𝑛+1/2_𝑖+1,𝑗−𝑉^𝑛+1/2_{𝑖−1,𝑗}ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1+𝑟𝑉^𝑛+1/2_𝑖𝑗,𝐿^𝑀_𝐴𝑉^𝑛_𝑖𝑗≡−𝑥_𝑖𝑉^𝑛_𝑖,𝑗+1−𝑉^𝑛_𝑖𝑗ℎ_𝐴,𝑗+1. (39)(40)

Lemma 3 (discrete maximum principle). The operator (𝐼+Δ𝑡𝐿^𝑁_𝑥) defined by (39) on the piecewise uniform mesh Ω^𝑁 satisfies a discrete maximum principle; that is, if 𝑢_𝑖 and 𝑤_𝑖 are mesh functions that satisfy 𝑢₀≥𝑤₀,𝑢_𝑁≥𝑤_𝑁, and (𝐼+Δ𝑡𝐿^𝑁_𝑥)𝑢_𝑖≥(𝐼+Δ𝑡𝐿^𝑁_𝑥)𝑤_𝑖(1≤𝑖<𝑁), then 𝑢_𝑖≥𝑤_𝑖 for all 𝑖.

Proof. Let 𝑎_𝑖𝜎=−²𝑥²_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1ℎ_𝑥,𝑖+𝑟𝑥_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1,𝑏_𝑖=𝜎²𝑥²_𝑖Δ𝑡ℎ_𝑥,𝑖ℎ_𝑥,𝑖+1𝑐+𝑟Δ𝑡+1,_𝑖𝜎=−²𝑥²_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1ℎ_𝑥,𝑖+1−𝑟𝑥_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1,1≤𝑖<𝑁. (41) We can obtain 𝑎_𝑖𝜎<−²𝑥₁𝑥_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1ℎ_𝑥,𝑖+𝑟𝑥_𝑖Δ𝑡ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1≤−𝛼𝑥₁+𝛽^∗ℎ_𝑥,𝑖ℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1ℎ_𝑥,𝑖𝑥_𝑖=Δ𝑡−𝛼ℎ+𝛽^∗𝛼/𝛽^∗ℎℎ_𝑥,𝑖+ℎ_𝑥,𝑖+1ℎ_𝑥,𝑖𝑥_𝑖Δ𝑡=0,2≤𝑖<𝑁. (42) Clearly, 𝑏_𝑖𝑐>0for1≤𝑖≤𝑁−1,_𝑖𝑏<0for1≤𝑖≤𝑁−2,₁+𝑐₁𝑎>0,_𝑖+𝑏_𝑖+𝑐_𝑖𝑎>0,2≤𝑖≤𝑁−2,_𝑁−1+𝑏_𝑁−1>0. (43) Hence, we verify that the matrix associated with (𝐼+Δ𝑡𝐿^𝑁_𝑥) is an M-matrix.

The local error of the difference scheme (34) is given by 𝜏_𝑖=𝐼+Δ𝑡𝐿^𝑁_𝑥𝑣^𝑛+1/2𝑥_𝑖−𝐼+Δ𝑡𝐿_𝑥𝑣^𝑛+1/2𝑥_𝑖, (44) for 1≤𝑖<𝑁, where the dependence on the parameter 𝐴_𝑗 is omitted. Applying the Taylor formula, we get the following estimate for the truncation error: 𝜏_𝑖≤𝐶𝑁⁻²Δ𝑡,1≤𝑖<𝑁, (45) where 𝐶 is a positive constant independent of the mesh. Hence, applying the discrete maximum principle we can get the following error estimates.

Lemma 4. Let 𝑣^𝑛+1/2 be the solution of the problem (16)–(18) and 𝑉^𝑛+1/2 be the solution of the problem (34)–(36), then one has the error estimates 𝑣^𝑛+1/2𝑥_𝑖,𝐴_𝑗−𝑉^𝑛+1/2_𝑖𝑗≤𝐶𝑁⁻²Δ𝑡,1≤𝑖≤𝑁,0≤𝑗≤𝑀, (46) where 𝐶 is a constant independent of 𝑁 and Δ𝑡.

In the second half step, we have problem (14), whose discretization is 𝐼+Δ𝑡𝐿^𝑀_𝐴𝑉^𝑛_𝑖𝑗=𝑉^𝑛+1/2_𝑖𝑗,1≤𝑖≤𝑁,0≤𝑗<𝑀,𝑉^𝑛_𝑖,𝑀=𝑥_𝑖𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖≤𝑁. (47) In order to find the relation between 𝑣^𝑛 and 𝑉^𝑛, we introduce the auxiliary problem 𝐼+Δ𝑡𝐿^𝑀_𝐴𝑉^𝑛_𝑖𝑗=𝑣^𝑛+1/2_𝑖𝑗𝑉,1≤𝑖≤𝑁,0≤𝑗<𝑀,_{𝑛𝑖,𝑀}=𝑥_𝑖𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖≤𝑁. (48)

Lemma 5. Let 𝑣^𝑛 be the solution of the problem (19)-(20) and 𝑉^𝑛 the solution of the problem (48); then one has the error estimates 𝑣^𝑛𝑥_𝑖,𝐴_𝑗−𝑉^𝑛_𝑖𝑗≤𝐶𝑀⁻¹Δ𝑡,1≤𝑖≤𝑁,0≤𝑗≤𝑀, (49) where 𝐶 is a constant independent of 𝑁,𝑀, and Δ𝑡.

Proof. It is easy to see that the matrix associated with (𝐼+Δ𝑡𝐿^𝑀_𝐴) is an M-matrix. Hence, applying the discrete maximum principle we can obtain the desired results.

Since 𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗=𝐼+Δ𝑡𝐿^𝑀_𝐴⁻¹𝑣^𝑛+1/2_𝑖𝑗−𝑉^𝑛+1/2_𝑖𝑗,𝐼+Δ𝑡𝐿^𝑀_𝐴⁻¹_{Ω𝑁×Ω𝑀}≤1, (50) we can obtain 𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗≤𝐶𝑁⁻²Δ𝑡, (51) where we have used Lemma 4.

Noting that 𝑣^𝑛𝑥_𝑖,𝐴_𝑗−𝑉^𝑛_𝑖𝑗=𝑣^𝑛𝑥_𝑖,𝐴_𝑗−𝑉^𝑛_𝑖𝑗+𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗, (52) we can get the following error estimates by (49) and (51).

Theorem 6. Let 𝑣^𝑛 be the solution of the problem (16)–(20) and 𝑉^𝑛 the solution of the problem (34)–(38); then one has the error estimates 𝑣^𝑛𝑥_𝑖,𝐴_𝑗−𝑉^𝑛_𝑖𝑗≤𝐶𝑁⁻²+𝑀⁻¹Δ𝑡,1≤𝑖≤𝑁,0≤𝑗≤𝑀, (53) where 𝐶 is a constant independent of 𝑁,𝑀, and Δ𝑡.

5. The Fully Discrete Scheme

Combining the time semidiscretization scheme (12)–(14) and the spatial discretization scheme (34)–(38), the following fully discrete scheme is deduced: 𝑉^𝐾_𝑖𝑗𝐴=max^𝑗𝑇−𝐸,0,1≤𝑖≤𝑁,0≤𝑗<𝑀,𝐼+Δ𝑡𝐿^𝑁_𝑥𝑉^𝑛+1/2_𝑖𝑗=𝑉^𝑛+1_𝑖𝑗𝑉,1≤𝑖≤𝑁,0≤𝑗<𝑀,_{𝑛+1/20,𝑗}=𝑒^{−𝑟(𝑇−𝑡𝑛)}𝐴max^𝑗𝑇𝑉−𝐸,0,1≤𝑖<𝑁,_{𝑛+1/2𝑁,𝑗}=𝑋𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)}+𝐴_𝑗𝑇−𝐸𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖<𝑁,𝐼+Δ𝑡𝐿^𝑀_𝐴𝑉^𝑛_𝑖𝑗=𝑉^𝑛+1/2_𝑖𝑗𝑉,1≤𝑖≤𝑁,0≤𝑗<𝑀,_{𝑛𝑖,𝑀}=𝑥_𝑖𝑟𝑇1−𝑒^{−𝑟(𝑇−𝑡𝑛)},1≤𝑖≤𝑁,for𝑛=𝐾−1,…,1,0, (54) where the discrete operators 𝐿^𝑁_𝑥,𝐿^𝑀_𝐴 are described in Section 4 and 𝑉^𝑛_𝑖𝑗 is the fully discrete approximation to the exact solution of (9) at the mesh point (𝑥_𝑖,𝐴_𝑗,𝑡_𝑛).

Theorem 7 (global error). Let 𝑣(𝑥,𝐴,𝑡) be the exact solution of (9) and 𝑉 the discrete solution of the fully discrete scheme (54). Then, there exists a positive constant 𝐶, independent of 𝑁,𝑀, and Δ𝑡, such that the global error satisfies 𝑣𝑥_𝑖,𝐴_𝑗,𝑡_𝑛−𝑉^𝑛_{𝑖𝑗Ω𝑁,𝑀,𝐾}≤𝐶𝑁⁻²+𝑀⁻¹+Δ𝑡. (55)

Proof. Splitting the global error in the form 𝑣𝑥_𝑖,𝐴_𝑗,𝑡_𝑛−𝑉^𝑛_{𝑖𝑗Ω𝑁,𝑀,𝐾}≤𝑣𝑥_𝑖,𝐴_𝑗,𝑡_𝑛−𝑣^𝑛_{𝑖𝑗Ω𝑁×Ω𝑀}+𝑣^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗Ω^𝑁_×Ω^𝑀+𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗Ω^𝑁_×Ω^𝑀_. (56) From Lemma 1 and Theorem 6, we deduce 𝑣𝑥_𝑖,𝐴_𝑗,𝑡_𝑛−𝑉^𝑛_{𝑖𝑗Ω𝑁,𝑀,𝐾}≤𝐶𝑁⁻²+𝑀⁻¹+Δ𝑡Δ𝑡+𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗Ω^𝑁_×Ω^𝑀_. (57) To bound the last term of (57), we take into account that {𝑉^𝑛−𝑉^𝑛} can be written as the solution of one step of (54), taking {𝑣(𝑥_𝑖,𝐴_𝑗,𝑡_𝑛+1)−𝑉^𝑛+1_𝑖𝑗} as final value. Applying the stability of the discrete scheme, we have 𝑉^𝑛_𝑖𝑗−𝑉^𝑛_𝑖𝑗Ω^𝑁_×Ω^𝑀≤𝑣𝑥_𝑖,𝐴_𝑗,𝑡_𝑛+1−𝑉^𝑛+1_{𝑖𝑗Ω𝑁×Ω𝑀}. (58) Then, from (57) and (58) a recurrence relation for the global errors follows, and from it the result of Theorem 7 follows immediately.

6. Numerical Experiments

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

Test 1. Fixed strike Asian call option with parameters: 𝜎=0.2,𝑟=0.08,𝑇=1, and 𝐸=2.

Test 2. Fixed strike Asian call option with parameters: 𝜎=0.4,𝑟=0.06,𝑇=1, and 𝐸=2.

In order to appropriately select 𝑋, we compute the value of Asian option at one mesh point for different 𝑋. From Table 1, we see that the value 𝑋=4𝐸 is large enough to guarantee the fact that the price does not depend on the position of 𝑋. Hence, in our numerical experiments we choose 𝑋=4𝐸.

Table 1 

Comparison of the results at mesh point (𝐸,𝐾𝑇/2,0) for different values of 𝑋 for fixed parameters: 𝜎=0.2,𝑟=0.08,𝑇=1,𝐸=2,𝑁=1024, and 𝑀=𝐾=256.

For Test 1, the computed Asian option value 𝑉 at 𝑡=0.5 with 𝑁=𝑀=𝐾=64 is depicted in Figure 1 for the asset price running sum 𝐴 between 0 and 2, since the exact option value is given by (8) for 𝐴≥𝐸𝑇=2.

Figure 1 

Computed Asian option value 𝑉 for Test 1.

For Test 2, the computed Asian option value 𝑉 at 𝑡=0.5 with 𝑁=𝑀=𝐾=64 is depicted in Figure 2 for the asset price running sum 𝐴 between 0 and 2.

Figure 2 

Computed Asian option value 𝑉 for Test 2.

To demonstrate the theoretical rates of convergence numerically, we take 𝑀=𝐾=512 which is a sufficiently large choice so that the error is dominated by the 𝑥 variable discretization. The exact solutions of our test problems for the asset price running sum 𝐴 between 0 and 𝐸𝑇 are not available. We use the approximated solution of 𝑁=1024,𝑀=𝐾=512 as the exact solution. We present the error estimates for different 𝑁. Because we only know “the exact solution” on mesh points, we use the linear interpolation to get solutions at other points. In this paper, 𝑉(𝑥,𝐴,𝑡) denotes “the exact solution” which is a linear interpolation of the approximated solution 𝑉^1024,512,512. We measure the accuracy in the discrete maximum norm 𝑒^{𝑁,𝑀,𝐾}=max_{𝑖,𝑗,𝑛}𝑉^{𝑁,𝑀,𝐾}_{𝑖,𝑗,𝑛}−𝑉𝑥_𝑖,𝐴_𝑗,𝑡_𝑛, (59) and the convergence rate 𝑅^{𝑁,𝑀,𝐾}=log₂𝑒^{𝑁,𝑀,𝐾}𝑒^{2𝑁,𝑀,𝐾}. (60) The error estimates and convergence rates in our computed solutions of Tests 1 and 2 are listed in Tables 1 and 2, respectively.

Table 2 

Numerical results for Test 1.

From the figures, it is seen that the numerical solutions by our method are nonoscillatory. From Tables 2 and 3, we see that 𝑅^{𝑁,𝑀,𝐾} is close to 2 for sufficiently large 𝑀 and 𝐾, which supports the convergence estimate of Theorem 7. They indicate that the theoretical results are fairly sharp.

Table 3 

Numerical results for Test 2.

Finally, in order to further confirm the accuracy of our method we compare our numerical results with the analytical solutions for 𝐸𝑇≤𝐴≤2𝐸𝑇. We give the absolute errors of the option values and the analytical solutions for 𝐸𝑇≤𝐴≤2𝐸𝑇 in Table 4 when 𝜎=0.2,𝑟=0.08,𝑇=1,𝐸=2,𝑋=8, and Ω₂=(0,2𝐸𝑇). From Table 4, we see that 𝑒^{𝑁,𝑀,𝐾}/𝑒^{2𝑁,4𝑀,4𝐾} is close to 4, which indicates that our method is second-order convergent with respect to the asset price and is first-order convergent with respect to both time variable and spatial variable 𝐴.

Table 4 

Comparisons of our numerical solutions with the analytical solutions.

7. Conclusion

In this paper, a finite difference scheme to solve the two-dimensional PDE arising in pricing arithmetic Asian options is proposed. To avoid solving a large discrete linear system, the ADI scheme is applied to Asian option pricing; that is, at each time-step, a one-dimensional Black-Scholes equation and an advection equation are alternately solved. Since the standard finite difference method to solve the Black-Scholes equation leads to numerical difficulty, we apply the central difference scheme on a piecewise uniform mesh to discrete the one-dimensional Black-Scholes equation and apply the implicit Euler scheme to discrete the advection equation. It is proved that the matrices associated with discrete operators are M-matrices, which ensures that the numerical scheme is stable in the maximum norm and free of nonphysical oscillations whether 𝜎²/𝑟≤1 or 𝜎²/𝑟>1. We show that the scheme is second-order convergent with respect to the asset price. Numerical results support the theoretical results.

Acknowledgments

The authors would like to thank the anonymous referees 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 (Grants nos. 2012A610035, 2012A610036), 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.