# An Alternating-Direction Implicit Difference Scheme for Pricing Asian Options

Authors: Zhongdi Cen; Anbo Le; Aimin Xu
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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 x 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 x=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 areM-matrices, which ensures that the numerical scheme is stable in the maximum norm and free of nonphysical oscillations whether σ2/r≤1 or σ2/r>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 Section3, 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 ∥·∥D, where D is a closed and bounded set.


## 2. The Continuous Problem

Suppose that the underlying asset pricex(t) follows a geometric Brownian motion
(1)dx(t)=rx(t)dt+σx(t)dB(t),
where r  is risk-free interest rate, σ is the volatility, and B(t) is a standard Brownian motion under the risk-neutral measure ℙ. Let A(t) denote the underlying asset price running sum given by
(2)A(t)=∫0t‍x(τ)dτ;
then the arithmetically averaged price is given by A(t)/t. The stochastic differential equation for evolution of (A(t)) is given by
(3)dA(t)=x(t)dt.
Thus, the Asian call option price v(x,A,t) for continuous arithmetic average strike satisfies the following PDE [15, 18]:
(4)Lv(x,A,t)≡-∂v∂t-12σ2x2∂2v∂x2-rx∂v∂x-x∂v∂A+rv=0,
with the terminal condition
(5)v(x,A,T)=max(AT-E,0),
and the left-hand boundary condition
(6)v(0,A,t)=e-r(T-t)max(AT-E,0),
where T is the expiry date and E is the strike price. Note that in the work of Geman and Yor [3], a simple formula is obtained for A≥ET:
(7)v(x,A,t)=xrT(1-e-r(T-t))+(AT-E)e-r(T-t).
From this formula the right-hand boundary condition can be obtained:
(8)v(x,A,t)~xrT(1-e-r(T-t))+e-r(T-t)(AT-E),asx⟶+∞.
See Hugger [19] for the derivation details of boundary value conditions for Asian options.As the exact solution to the problem (4)–(6) for A≥ET is known, we only consider the solution of the PDE for 0<A<ET using (7) for the boundary condition at A=ET. For applying the numerical method, we truncate the domain (0,+∞) of spatial variable x into (0,X) for sufficiently large X. 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 X=4E. The boundary condition at x=X is chosen to be v(X,A,t)=(X/rT)(1-e-r(T-t))+(A/T-E)e-r(T-t). 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:
(9)Lv(x,A,t)=0,(x,A,t)∈Ω,v(x,A,T)=max(AT-E,0),(x,A)∈Ω1×Ω2,v(0,A,t)=e-r(T-t)max(AT-E,0),(A,t)∈Ω2×Ω3,v(X,A,t)=XrT(1-e-r(T-t))+(AT-E)e-r(T-t),(A,t)∈Ω2×Ω3,v(x,ET,t)=xrT(1-e-r(T-t)),(x,t)∈Ω1×Ω3,
where Ω1=(0,X),Ω2=(0,ET),Ω3=(0,T) and Ω=Ω1×Ω2×Ω3.


## 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 thex and A directions.We discrete the above PDE in[0,T] using a uniform mesh, which is defined by
(10)Ω-K={tn=nΔt,0≤n≤K,Δt=T/K}.
Split the spatial differential operator into the following two operators:
(11)Lx≡-12σ2x2∂2∂x2-rx∂∂x+rv,LA=-x∂∂A.
Using the partitions, the semidiscrete scheme can be written as follows:(a)
(12)vK=max(AT-E,0),(b)
(13)(I+ΔtLx)vn+1/2=vn+1,vn+1/2(0,A)=e-r(T-tn)max(AT-E,0),vn+1/2(X,A)=XrT(1-e-r(T-tn))+(AT-E)e-r(T-tn),(c)
(14)(I+ΔtLA)vn=vn+1/2,vn(x,ET)=xrT(1-e-r(T-tn)).This method gives approximations vn(x,A) to the solution v(x,A,t) of (9) at the time levels tn. The operators (I+ΔtLx) and (I+ΔtLA) satisfy a maximum principle, and consequently
(15)∥(I+ΔtLx)-1∥Ω-1×Ω-2≤11+rΔt,∥(I+ΔtLA)-1∥Ω-1×Ω-2≤1.
This ensures the stability of the temporal semidiscretization (12)–(14) and that each step of the scheme (12)–(14) has a unique solution vn(x,A).The local truncation error is defined asen=v(tn)-v-n, where v-n is the solution of
(16)(I+ΔtLx)v-n+1/2=v(x,A,tn+1),(17)v-n+1/2(0,A)=e-r(T-tn)max(AT-E,0),(18)v-n+1/2(X,A)=XrT(1-e-r(T-tn))+(AT-E)e-r(T-tn),(19)(I+ΔtLA)v-n=v-n+1/2,(20)v-n(x,ET)=xrT(1-e-r(T-tn)).
The following lemma gives the local error estimates.Lemma 1.
The local error for the scheme (16)–(20) satisfies
(21)∥en∥Ω-1×Ω-2≤C(Δt)2.Proof.
The functionv-n satisfies
(22)(I+ΔtLx)(I+ΔtLA)v-n(x,A)=v(x,A,tn+1).
On the other hand, since the solution of (4) is smooth enough, we have
(23)v(x,A,tn+1)=v(x,A,tn)+Δt(Lx+LA)v(x,A,tn)+∫tntn+1‍(tn+1-s)∂2v∂s2ds=(I+ΔtLx)(I+ΔtLA)v(x,A,tn)+O((Δt)2).
Hence, en satisfies
(24)(I+ΔtLx)(I+ΔtLA)en=O((Δt)2),en(0,A)=en(X,A)=en(x,ET)=0.
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 En=v(tn)-vn, satisfies
(25)sup0≤n<K∥En∥Ω-1×Ω-2≤CΔt,
and therefore the time semidiscretization method is a first-order convergent scheme.Proof.
It is easy to see that(26)En=en+REn+1,
where
(27)R=(I+ΔtLx)-1(I+ΔtLA)-1,
is the transition operator and REn+1 is the result obtained after one step of scheme (13)-(14) with En+1 as final value. Using this recurrence, we deduce that
(28)En=∑i=nK-1‍Ri-nei.
From
(29)∥R∥Ω-1×Ω-2=∥(I+ΔtLx)-1∥Ω-1×Ω-2·∥(I+ΔtLA)-1∥Ω-1×Ω-2≤11+rΔt,
we obtain the desired result.


## 4. The Spatial Discretization

As discussed in Section1, 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ΩN on the space interval [0,X]:
(30)xi={hi=1,h[1+σ2r(i-1)]i=2,…,N,
where
(31)h=X1+(σ2/r)(N-1).
For the A variable discretization, we use a uniform mesh ΩM on [0,ET] with M mesh elements. It is easy to see that the mesh sizes hx,i=xi-xi-1 and hA,j=Aj-Aj-1 satisfy
(32)hx,i={hfori=1,σ2rhfori=2,…,N,(33)hA,j=ETM,j=1,…,M,
respectively. We denote by Ω-N,M,K=Ω-N×Ω-M×Ω-K the corresponding mesh for x,A, and t.Thus, we apply the central difference scheme on the piecewise uniform mesh and the implicit Euler scheme to approximate problem (16)–(20): (34)(I+ΔtLxN)V-ijn+1/2=vijn+1,1≤i≤N,0≤j<M,(35)V-0,jn+1/2=e-r(T-tn)max(AjT-E,0),1≤i<N,(36)V-N,jn+1/2=XrT(1-e-r(T-tn))+(AjT-E)e-r(T-tn),1≤i<N,(37)(I+ΔtLAM)V-ijn=V-ijn+1/2,1≤i≤N,0≤j<M,(38)V-i,Mn=xirT(1-e-r(T-tn)),1≤i≤N,
where
(39)LxNV-ijn+1/2≡-σ2xi2hx,i+hx,i+1×(V-i+1,jn+1/2-V-ijn+1/2hx,i+1-V-ijn+1/2-V-i-1,jn+1/2hx,i)-rxiV-i+1,jn+1/2-V-i-1,jn+1/2hx,i+hx,i+1+rV-ijn+1/2,(40)LAMV-ijn≡-xiV-i,j+1n-V-ijnhA,j+1.Lemma 3 (discrete maximum principle).
The operator(I+ΔtLxN) defined by (39) on the piecewise uniform mesh ΩN satisfies a discrete maximum principle; that is, if ui and wi are mesh functions that satisfy u0≥w0,uN≥wN, and (I+ΔtLxN)ui≥(I+ΔtLxN)wi(1≤i<N), then ui≥wi for all i.Proof.
Let(41)ai=-σ2xi2Δt(hx,i+hx,i+1)hx,i+rxiΔthx,i+hx,i+1,bi=σ2xi2Δthx,ihx,i+1+rΔt+1,ci=-σ2xi2Δt(hx,i+hx,i+1)hx,i+1-rxiΔthx,i+hx,i+1,1≤i<N.
We can obtain
(42)ai<-σ2x1xiΔt(hx,i+hx,i+1)hx,i+rxiΔthx,i+hx,i+1≤(-αx1+β*hx,i)(hx,i+hx,i+1)hx,ixiΔt=(-αh+β*(α/β*)h)(hx,i+hx,i+1)hx,ixiΔt=0,2≤i<N.
Clearly,
(43)bi>0for1≤i≤N-1,ci<0for1≤i≤N-2,b1+c1>0,ai+bi+ci>0,2≤i≤N-2,aN-1+bN-1>0.
Hence, we verify that the matrix associated with (I+ΔtLxN) is an M-matrix.The local error of the difference scheme (34) is given by
(44)τi=(I+ΔtLxN)(v-n+1/2(xi))-(I+ΔtLx)(v-n+1/2(xi)),
for 1≤i<N, where the dependence on the parameter Aj is omitted. Applying the Taylor formula, we get the following estimate for the truncation error:
(45)|τi|≤CN-2Δt,1≤i<N,
where C is a positive constant independent of the mesh. Hence, applying the discrete maximum principle we can get the following error estimates.Lemma 4.
Letv-n+1/2 be the solution of the problem (16)–(18) and V-n+1/2 be the solution of the problem (34)–(36), then one has the error estimates
(46)|v-n+1/2(xi,Aj)-V-ijn+1/2|≤CN-2Δt,1≤i≤N,0≤j≤M,
where C is a constant independent of N and Δt.In the second half step, we have problem (14), whose discretization is
(47)(I+ΔtLAM)V-ijn=V-ijn+1/2,1≤i≤N,0≤j<M,V-i,Mn=xirT(1-e-r(T-tn)),1≤i≤N.
In order to find the relation between v-n and V-n, we introduce the auxiliary problem
(48)(I+ΔtLAM)V~ijn=v-ijn+1/2,1≤i≤N,0≤j<M,V~i,Mn=xirT(1-e-r(T-tn)),1≤i≤N.Lemma 5.
Letv-n be the solution of the problem (19)-(20) and V~n the solution of the problem (48); then one has the error estimates
(49)|v-n(xi,Aj)-V~ijn|≤CM-1Δt,1≤i≤N,0≤j≤M,
where C is a constant independent of N,M, and Δt.Proof.
It is easy to see that the matrix associated with(I+ΔtLAM) is an M-matrix. Hence, applying the discrete maximum principle we can obtain the desired results.Since(50)V~ijn-V-ijn=(I+ΔtLAM)-1(v-ijn+1/2-V-ijn+1/2),∥(I+ΔtLAM)-1∥Ω-N×Ω-M≤1,
we can obtain
(51)|V~ijn-V-ijn|≤CN-2Δt,
where we have used Lemma 4.Noting that(52)v-n(xi,Aj)-V-ijn=v-n(xi,Aj)-V~ijn+V~ijn-V-ijn,
we can get the following error estimates by (49) and (51).Theorem 6.
Letv-n be the solution of the problem (16)–(20) and V-n the solution of the problem (34)–(38); then one has the error estimates
(53)|v-n(xi,Aj)-V-ijn|≤C(N-2+M-1)Δt,1≤i≤N,0≤j≤M,
where C is a constant independent of N,M, and Δt.


## 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:
(54)VijK=max(AjT-E,0),1≤i≤N,0≤j<M,(I+ΔtLxN)Vijn+1/2=Vijn+1,1≤i≤N,0≤j<M,V0,jn+1/2=e-r(T-tn)max(AjT-E,0),1≤i<N,VN,jn+1/2=XrT(1-e-r(T-tn))+(AjT-E)e-r(T-tn),1≤i<N,(I+ΔtLAM)Vijn=Vijn+1/2,1≤i≤N,0≤j<M,Vi,Mn=xirT(1-e-r(T-tn)),1≤i≤N,forn=K-1,…,1,0,
where the discrete operators LxN,LAM are described in Section 4 and Vijn is the fully discrete approximation to the exact solution of (9) at the mesh point (xi,Aj,tn).Theorem 7 (global error).
Letv(x,A,t) be the exact solution of (9) and V the discrete solution of the fully discrete scheme (54). Then, there exists a positive constant C, independent of N,M, and Δt, such that the global error satisfies
(55)∥v(xi,Aj,tn)-Vijn∥Ω-N,M,K≤C(N-2+M-1+Δt).Proof.
Splitting the global error in the form(56)∥v(xi,Aj,tn)-Vijn∥Ω-N,M,K≤∥v(xi,Aj,tn)-v-ijn∥Ω-N×Ω-M+∥v-ijn-V-ijn∥Ω-N×Ω-M+∥V-ijn-Vijn∥Ω-N×Ω-M.
From Lemma 1 and Theorem 6, we deduce
(57)∥v(xi,Aj,tn)-Vijn∥Ω-N,M,K≤C(N-2+M-1+Δt)Δt+∥V-ijn-Vijn∥Ω-N×Ω-M.
To bound the last term of (57), we take into account that {V-n-Vn} can be written as the solution of one step of (54), taking {v(xi,Aj,tn+1)-Vijn+1} as final value. Applying the stability of the discrete scheme, we have
(58)∥V-ijn-Vijn∥Ω-N×Ω-M≤∥v(xi,Aj,tn+1)-Vijn+1∥Ω-N×Ω-M.
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,r=0.08,T=1, and E=2.Test 2.
Fixed strike Asian call option with parameters:σ=0.4,r=0.06,T=1, and E=2.In order to appropriately selectX, we compute the value of Asian option at one mesh point for different X. From Table 1, we see that the value X=4E is large enough to guarantee the fact that the price does not depend on the position of X. Hence, in our numerical experiments we choose X=4E.Table 1
Comparison of the results at mesh point(E,KT/2,0) for different values of X for fixed parameters: σ=0.2,r=0.08,T=1,E=2,N=1024, and M=K=256.
X
X
=
4
E
X
=
5
E
X
=
6
E
X
=
7
E
X
=
8
E
Price
0.9988279010
0.9988279011
0.9988279013
0.9988279017
0.9988279017For Test1, the computed Asian option value V at t=0.5 with N=M=K=64 is depicted in Figure 1 for the asset price running sum A between 0 and 2, since the exact option value is given by (8) for A≥ET=2.Figure 1
Computed Asian option valueV for Test 1.For Test2, the computed Asian option value V at t=0.5 with N=M=K=64 is depicted in Figure 2 for the asset price running sum A between 0 and 2.Figure 2
Computed Asian option valueV for Test 2.To demonstrate the theoretical rates of convergence numerically, we takeM=K=512 which is a sufficiently large choice so that the error is dominated by the x variable discretization. The exact solutions of our test problems for the asset price running sum A between 0 and ET are not available. We use the approximated solution of N=1024,M=K=512 as the exact solution. We present the error estimates for different N. 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^(x,A,t) denotes “the exact solution” which is a linear interpolation of the approximated solution V1024,512,512. We measure the accuracy in the discrete maximum norm
(59)eN,M,K=maxi,j,n|Vi,j,nN,M,K-V^(xi,Aj,tn)|,
and the convergence rate
(60)RN,M,K=log2(eN,M,Ke2N,M,K).
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 Test1.
M
=
K
N
Error
Rate
512
64
8.0147
e
-
2
—
128
3.9103
e
-
2
1.035
256
1.0132
e
-
2
1.948
512
2.1007
e
-
3
2.270From the figures, it is seen that the numerical solutions by our method are nonoscillatory. From Tables2 and 3, we see that RN,M,K is close to 2 for sufficiently large M and K, which supports the convergence estimate of Theorem 7. They indicate that the theoretical results are fairly sharp.Table 3
Numerical results for Test2.
M
=
K
N
Error
Rate
512
64
3.9971
e
-
2
—
128
1.0620
e
-
2
1.912
256
2.6179
e
-
3
2.020
512
5.3662
e
-
4
2.286Finally, in order to further confirm the accuracy of our method we compare our numerical results with the analytical solutions forET≤A≤2ET. We give the absolute errors of the option values and the analytical solutions for ET≤A≤2ET in Table 4 when σ=0.2,r=0.08,T=1,E=2,X=8, and Ω2=(0,2ET). From Table 4, we see that eN,M,K/e2N,4M,4K 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 A.Table 4
Comparisons of our numerical solutions with the analytical solutions.
M
=
K
N
Error
Rate
64
32
8.9402
e
-
2
—
256
64
2.5838
e
-
2
1.791
1024
128
6.8437
e
-
3
1.917
4048
256
1.7573
e
-
3
1.961


## 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σ2/r≤1 or σ2/r>1. We show that the scheme is second-order convergent with respect to the asset price. Numerical results support the theoretical results.



---

*Source: 605943-2013-06-18.xml*