An Alternating-Direction Implicit Difference Scheme for Pricing Asian Options

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.


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 twodimensional 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 convectiondominated 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 firstorder 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  2 / ≤ 1 or  2 / > 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.

The Continuous Problem
Suppose that the underlying asset price () follows a geometric Brownian motion where  is risk-free interest rate,  is the volatility, and () is a standard Brownian motion under the risk-neutral measure P. Let () denote the underlying asset price running sum given by then the arithmetically averaged price is given by ()/.The stochastic differential equation for evolution of (()) is given by Thus, the Asian call option price V (, , ) for continuous arithmetic average strike satisfies the following PDE [15,18]: with the terminal condition and the left-hand boundary condition 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  ≥ : From this formula the right-hand boundary condition can be obtained: 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 V(, , ) = (/)(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: where Ω 1 = (0, ), Ω 2 = (0, ), Ω 3 = (0, ) and Ω = Ω 1 × Ω 2 × Ω 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 Split the spatial differential operator into the following two operators: Using the partitions, the semidiscrete scheme can be written as follows: (a) This method gives approximations V  (, ) to the solution V(, , ) of ( 9) at the time levels   .The operators ( + Δ  ) and (+Δ  ) satisfy a maximum principle, and consequently This ensures the stability of the temporal semidiscretization ( 12)-( 14) and that each step of the scheme ( 12)-( 14) has a unique solution V  (, ).
The local truncation error is defined as The following lemma gives the local error estimates.
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 and therefore the time semidiscretization method is a first-order convergent scheme.
Proof.It is easy to see that where 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 From we obtain the desired result.

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, ]: where 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 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 ≤  ≤ , 0 ≤  < , (37) where , We can obtain Clearly, Hence, we verify that the matrix associated with ( + Δ   ) is an M-matrix.
The local error of the difference scheme (34) is given by for 1 ≤  < , where the dependence on the parameter   is omitted.Applying the Taylor formula, we get the following estimate for the truncation error: 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 V +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 where  is a constant independent of  and Δ.
In the second half step, we have problem ( 14), whose discretization is In order to find the relation between V  and   , we introduce the auxiliary problem Lemma 5. Let V  be the solution of the problem ( 19)-( 20) and Ṽ the solution of the problem (48); then one has the error estimates 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.

Ṽ𝑛
− we can obtain where we have used Lemma 4. Noting that we can get the following error estimates by (49) and (51).
Theorem 6.Let V  be the solution of the problem (16)- (20) and   the solution of the problem (34)-(38); then one has the error estimates where  is a constant independent of , , and Δ.

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: 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 (  ,   ,   ).

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.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.
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.
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.
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, V(, , ) 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 ,,       ,, ,, − V (  ,   ,   )      , and the convergence rate  ,, = log 2 (  ,,  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.
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.
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 Ω 2 = (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 .

Conclusion
In this paper, a finite difference scheme to solve the twodimensional 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  2 / ≤ 1 or  2 / > 1.We show that the scheme is second-order convergent with respect to the asset price.Numerical results support the theoretical results.

Figure 1 :
Figure 1: Computed Asian option value  for Test 1.

Figure 2 :
Figure 2: Computed Asian option value  for Test 2.

Table 2 :
Numerical results for Test 1.

Table 3 :
Numerical results for Test 2.

Table 4 :
Comparisons of our numerical solutions with the analytical solutions.