- Email: [email protected]
Finite difference scheme with a moving mesh for pricing Asian options
Applied Mathematics and Computation 219 (2013) 8667–8675
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Finite difference scheme with a moving mesh for pricing Asian options Zhongdi Cen ⇑, Anbo Le, Aimin Xu Institute of Mathematics, Zhejiang Wanli University, Ningbo 315100, Zhejiang, PR China
a r t i c l e
i n f o
Keywords: Asian option Path-dependent option Finite difference Moving mesh method Rannacher time stepping scheme
a b s t r a c t In this paper we propose a stable numerical method for pricing Asian call options, which is based on a central difference scheme with a moving mesh in the spatial discretization and the Rannacher time stepping scheme in the time discretization. At each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We show that the scheme is second-order convergent with respect to both time and spatial variables. Numerical results support the theoretical results. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Asian options are securities whose payoffs depend on the average value of the underlying asset price over some time period. Asian options are used by traders who are interested to hedge against the average price of an asset over a period rather than the price at the end of period. 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. This averaging feature furthermore makes Asian options enjoy lower volatilities than their underlying assets, thus as a hedge instrument they are more economical than standard options on the same underlying assets. Since the payoff of an Asian option depends on the average value of an underlying asset price over some time interval, Asian options have proven to be much more difficult to value than standard options. Standard numerical techniques tend to be impractical, inaccurate, or slow, as discussed in [3,26]. For example, traditional lattice methods [9,15] require such enormous amounts of computer memory (owing to the necessity of keeping track of every possible path throughout the tree) that they are almost impracticable. Monte Carlo simulation [14,16] works well for Asian options, but it is computationally expensive without the enhancement of variance reduction techniques. The classical numerical partial differential equation (PDE) methods are inaccurate since the degeneration of the PDE for pricing Asian options causes numerical diffusion and spurious oscillation [2,26]. Zvan et al. [26] apply a flux-limiting method from computational fluid dynamics to tackle the problem of spurious oscillations that arise in Asian options. Since their method requires that nonlinear equations should be solved, the computational workloads would be heavy and its implementation could be complicated. Hugger [10] proposes an artificial viscosity numerical method for Asian options to avoid the oscillations. Adding artificial viscosity generally smoothens out the numerical solution, but also decreases the precision since a different problem is now considered. Dubois and Lelièvre [6] extend the approach by Rogers and Shi [20] and transform the PDE for an Asian option into the Vecˇerˇ’s PDE [23]. But the new one-dimensional PDE may be convection-dom-
⇑ Corresponding author. E-mail address: [email protected] (Z. Cen). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.02.065
8668
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
inated as the standard Black–Scholes equation when the ratio of diffusion to convection term is small. It is well known that the standard finite difference method for those convection-dominated problems may lead to nonphysical oscillations in the computed solution [24]. Tangman et al. [21] use the exponential time integration scheme in combination with a dimensional splitting technique for pricing Asian options under a variety of pricing models, including the Merton model and the Carr–Geman–Madan–Yor model. However this scheme is expensive for refined meshes. Mudzimbabwe et al. [17] propose an implicit finite difference method which is stable according to a von Neumann stability analysis. Wang [24] and Angermann and Wang [1] present a fitted finite volume method for the Black–Scholes equation, but the scheme is first-order convergent. We have proposed a robust difference scheme for the Black–Scholes equation, which is based on a central difference method for spatial discretization and the implicit Euler method for time discretization [4,5]. In this paper we extend this scheme to solve the Vecˇerˇ’s PDE for Asian options. We propose a stable numerical method for pricing Asian call options, which is based on a central difference scheme with a moving mesh in the spatial discretization and the Rannacher time stepping scheme in the time discretization. Specifically, at each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We will show that the scheme is second-order convergent with respect to both time and spatial variables. The rest of the paper is organised as follows. In the next section we describe some theoretical results on the Asian option pricing model. In Section 3 we prove the convergence of the Rannacher time stepping scheme. In Section 4 we consider the spatial discretization based on a central difference scheme with a moving mesh. It is shown that the difference scheme is stable for arbitrary volatility and arbitrary interest rate, and is second-order convergent with respect to both time and spatial variables. Finally, numerical experiments are provided to support these theoretical results in Section 5. 2. The continuous problem Let S ¼ SðtÞ denote the price of the underlying asset and IðtÞ denote the underlying asset price running sum given by
IðtÞ ¼
Z
t
SðsÞds:
0
The price CðS; I; tÞ of the continuous arithmetic average Asian call option with fixed strike price satisfies the following PDE [11,25]
@C 1 2 2 @ 2 C @C @C þ r S þS rC ¼ 0; þ rS 2 @t 2 @S @I @S
ð2:1Þ
with the final condition
CðS; I; TÞ ¼ max
I K; 0 ; T
ð2:2Þ
where r is the volatility of underlying asset, r is the risk-free interest rate, T is the expiry date, and K is the strike price. Note that the above problem is a two-dimensional PDE which leads to greater computational costs. As shown by Rogers and Shi [20], by making the change of variables
y¼
K I=T ; S
CðS; I; tÞ ¼ Sv ðy; tÞ;
ð2:3Þ
the two-dimensional PDE (2.1) can be reduced to a one-dimensional PDE
@v 1 2 2 @2v 1 @v þ ry þ r y ¼ 0; 2 T @t 2 @y @y
v ðy; TÞ ¼ maxðy; 0Þ:
ð2:4Þ ð2:5Þ
However, the above one-dimensional PDE is convection-dominated when the ratio of diffusion to convection term is small, which may lead to numerical difficulty for constructing a discretization with good properties and accuracy. As shown by Dubois and Lelièvre [6], by making the change of variables
t y¼x ; T
uðx; tÞ ¼ v ðy; tÞ;
ð2:6Þ
the PDE (2.4) and (2.5) can be transformed into the Vecˇerˇ’s PDE [23]
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
uðx; TÞ ¼ maxð1 x; 0Þ:
ð2:7Þ ð2:8Þ
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
8669
The PDE (2.7) and (2.8) is defined on the whole real axis. Geman and Yor [7] have given a formula for the case I P KT as
CðS; I; tÞ ¼
S I ð1 erðTtÞ Þ þ K erðTtÞ : rT T
ð2:9Þ
By making the change of variables as in (2.3) and (2.6), we obtain
uðx; tÞ ¼
1 t rðTtÞ ð1 erðTtÞ Þ x e rT T
for x 6
t : T
ð2:10Þ
Hence we only consider the solution of the PDE (2.7) and (2.8) for x > t=T using (2.10) for the boundary condition at x ¼ t=T. Therefore, we have the following PDE
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
ðx; tÞ 2 ðt=T; 1Þ ð0; TÞ;
uðx; TÞ ¼ maxð1 x; 0Þ; x 2 ½1; 1Þ; 1 uðt=T; tÞ ¼ ð1 erðTtÞ Þ; t 2 ½0; T; rT lim uðx; tÞ ¼ 0; t 2 ½0; T:
x!1
This PDE is also convection-dominated as the standard Black–Scholes equation when the convection term is larger than the diffusion term. We propose a stable finite difference scheme with a moving mesh to solve it. For applying the numerical method we truncate the domain ðt=T; 1Þ into X ðt=T; XÞ with an appropriate choice of X. The boundary condition at x ¼ X is chosen to be uðX; tÞ ¼ 0. Normally, this truncation of the domain leads to a negligible error in the value of the option [12]. Therefore, in the remaining of this paper we will consider the following PDE:
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
uðx; TÞ ¼ maxð1 x; 0Þ; uðt=T; tÞ ¼
x 2 ½1; 1Þ;
1 ð1 erðTtÞ Þ; rT
uðX; tÞ ¼ 0;
ðx; tÞ 2 X ð0; TÞ;
t 2 ½0; T;
t 2 ½0; T:
ð2:11Þ ð2:12Þ ð2:13Þ ð2:14Þ
Once the solution uðx; tÞ is obtained the price of the Asian option is determined by CðS; I; tÞ ¼ Suðx; tÞ. 3. The time semidiscretization As discussed in Giles and Carter [8] and Pooley et al. [18], the Crank–Nicolson method leads to serious degradation in the convergence rates of numerical schemes due to the non-smooth payoff function. The Rannacher time stepping scheme [19] uses a few first time steps with the implicit Euler method, and after that it uses the Crank–Nicolson method. They show that the Rannacher time stepping is simple to implement and greatly improves the stability of the numerical scheme in approximations of the Black–Scholes equations. We also use the Rannacher time stepping scheme to approximate the solution of (2.11)–(2.14). To ensure that the scheme M on the time interval ½0; T: is second-order accurate we use a piecewise uniform mesh x
( tj ¼
T4M 2 j; M4 T Mj ; M2
j ¼ 0; 1; . . . ; M 4; j ¼ M 3; . . . ; M:
ð3:1Þ
It is easy to see that the mesh sizes Mt j ¼ tj tj1 satisfies
( Mt j ¼
T4M2 M4 1 M2
j ¼ 0; 1; . . . ; M 4; j ¼ M 3; . . . ; M:
ð3:2Þ
Then the Rannacher time stepping scheme to approximate the solution of (2.11)–(2.14) reads
8 M u ¼ uðx; TÞ ¼ maxð1 x; 0Þ; > > > > h i > 8 > > < < I þ hj Mtjþ1 Ljx uj ¼ I ð1 hj ÞMt jþ1 Ljþ1 ujþ1 ; x > : j > > u ðt j =TÞ ¼ rT1 ð1 erðTtj Þ Þ; > > > > : for j ¼ M 1; . . . ; 1; 0;
uj ðXÞ ¼ 0;
ð3:3Þ
8670
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
where
Ljx uj ¼ ( hj ¼
2 2 j 1 2 t @ u t j @uj r x j þ r x ; 2 T @x2 T @x
1; j ¼ M 1; . . . ; M 4; 1 ; 2
j ¼ M 5; . . . ; 1; 0;
and uj ðxÞ denotes the approximation of the exact solution uðx; tÞ at the time level t j . Similarly to Kellogg and Tsan [13] we can prove that the differential operator ðI þ hj Mt jþ1 Ljx Þ satisfies a maximum principle, which we will use in the analysis of the convergence of the Rannacher time stepping scheme. Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problem
8 jþ1 ^ ¼ uðx; t jþ1 Þ; u > > > > h i > 8 > > < < I þ hj Mtjþ1 Ljx u ^ j ¼ I ð1 hj ÞMtjþ1 Ljþ1 ^ jþ1 ; u x > : j > > ^ ðt j =TÞ ¼ rT1 1 erðTtj Þ ; u > > > > : for j ¼ M 1; . . . ; 1; 0;
ð3:4Þ
^ j ðXÞ ¼ 0; u
is introduced to define the local error. ^ j ðxÞ satisfies Lemma 1 (Local error estimate). The local error associated to the method (3.1), defined by ej ¼ uðx; t j Þ u
kej kX ¼ OðM 3 Þ;
0 6 j 6 M;
ð3:5Þ
. where kej kX is the maximum norm of ej on the closed set X Proof. Using Taylor expansion we have
uðx; tjþ1 Þ uðx; tj Þ ¼ ut ðx; tj Þ þ OðMt jþ1 Þ ¼ Ljx uðx; tj Þ þ OðMt jþ1 Þ; Mt jþ1
ð3:6Þ
for j ¼ M 1; . . . ; M 4, and
uðx; tjþ1 Þ uðx; tj Þ ut ðx; t j Þ þ ut ðx; tjþ1 Þ ¼ ut ðx; tjþ1=2 Þ þ OððMt jþ1 Þ2 Þ ¼ þ OððMt jþ1 Þ2 Þ Mt jþ1 2 ¼
Ljx uðx; tj Þ þ Ljþ1 x uðx; t jþ1 Þ þ OððMt jþ1 Þ2 Þ 2
for j ¼ M 5; . . . ; 1; 0, where tjþ1=2 ¼ error is the solution of the problem
8 < I þ h Mt j
j jþ1 Lx
t j þt jþ1 . 2
ð3:7Þ
From (3.2)–(3.4) and (3.6) and (3.7), it is straightforward to show that the local
ej ¼ OðM 3 Þ;
: ej ðt =TÞ ¼ ej ðXÞ ¼ 0; j
and therefore the result follows from the maximum principle for the operator ðI þ hj Mtjþ1 Ljx Þ. h Lemma 2 (Global error estimate). The global error associated to the Rannacher time stepping scheme, given by Ej ¼ uðx; tj Þ uj ðxÞ, satisfies
E ¼ supkEj k1 6 CM 2 ;
ð3:8Þ
j6M
and therefore the time semidiscretization method is a second-order convergent scheme. Proof. Using the local error estimate down to the jth time step given by Lemma 1, we get the following global error estimate at the jth time step:
X M j kE k1 ¼ e 6 kej k1 þ kejþ1 k1 þ . . . þ keM k1 6 CM 2 ; k¼j j
1
for 0 6 j 6 M, where C is a positive constant independent of M. h
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
8671
4. The spatial discretization As discussed in introduction, the standard finite difference method to discretise the spatial derivatives of the problem (3.3), numerical difficulty can be caused. The main reason is that when the ratio of diffusion to convection term is small, the differential operator becomes convection-dominated. Since the coefficients of the derivatives in (3.3) are time-dependent, we apply the central difference scheme on a moving mesh to discretise problem (3.3). Specifically, at each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We will show that the scheme is second-order convergent with respect to the spatial variable. At the time mesh point t j , we use a piecewise uniform mesh XNj on the spatial interval ½tj =T; X:
8 tj j þh; i ¼ 1; > > > : 1 þ ði NÞ N1 ; i ¼ N þ 1; . . . ; N þ Q ;
ð4:1Þ
where j
h ¼
1 tj =T ; 2 1 þ rr ðN 1Þ
X ¼1þ
Q : N j
It is easy to see that the mesh sizes hi ¼ xji xji1 satisfy
8 j i ¼ 1; > :1 ; i ¼ N þ 1; . . . ; N þ Q: N
Thus, at each time mesh point tj we apply the central difference scheme on the above piecewise uniform mesh to approximate problem (3.3):
8 h i < I þ hj Mt jþ1 Lj;N U j ¼ I ð1 hj ÞMtjþ1 Ljþ1;N U jþ1 ; x x i i : U j ¼ 1 ð1 erðTtj Þ Þ; 0 rT
U jNþQ
1 6 i < N þ Q;
ð4:2Þ
¼ 0;
where
2 tj r2 j j Lj;N U ¼ x x i i j T h þ hj i
iþ1
U jiþ1 U ji j
hiþ1
U ji U ji1
!
j
hi
j j t j U iþ1 U i1 þ r xji ; j j T h þh i
ð4:3Þ
iþ1
and U jþ1 has been obtained by using the spline interpolation of the numerical results from the time mesh point tjþ1 . i N Lemma 3 (Discrete maximum principle). The operator ðI þ hj Mtjþ1 Lj;N x Þ defined by (4.3) on the piecewise uniform mesh Xj satisfies a discrete maximum principle, i.e. if v i and wi are mesh functions that satisfy v 0 P w0 ; v NþQ P wNþQ , and j;N ðI þ hj Mt jþ1 Lj;N x Þv i P ðI þ hj Mt jþ1 Lx Þwi ð1 6 i < N þ Q Þ, then v i P wi for all i.
Proof. Let
ai ¼
2 tj rhj Mtjþ1 tj j j x x ; i i j j j T T Þh h þh
r2 hj Mtjþ1 j
j
i
iþ1
ðhi þ hiþ1 i i iþ1 2 2 r hj Mtjþ1 j tj xi þ 1; 1 6 i < N þ Q; bi ¼ j j T hi hiþ1 2 r2 hj Mtjþ1 tj rhj Mtjþ1 tj j j : ci ¼ j x þ x i i j j j j T T ðh þ h Þh h þh iþ1
i
iþ1
We can obtain
r2 hj Mtjþ1 tj tj rhj Mt jþ1 tj hj Mt jþ1 tj tj j ci 6 xj1 xji r2 xj1 xji þ j ¼ þ rhiþ1 xji j j j j j j j T T T T T hi þ hiþ1 hi þ hiþ1 hiþ1 hi þ hiþ1 hiþ1 hj Mt jþ1 r2 j tj j ¼ 0; 1 6 i < N þ Q: xji r2 h þ r h ¼ j j j r T ðhi þ hiþ1 Þhiþ1
8672
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
Clearly,
ai < 0 for 2 6 i 6 N þ Q 1;
bi > 0 for 1 6 i 6 N þ Q 1;
and
b1 þ c1 > 0; ai þ bi þ ci > 0; 2 6 i 6 N þ Q 2; aNþQ 1 þ bNþQ1 > 0:
is an M-matrix. Thus by the same argument as Lemma 3.1 in Hence we verify that the matrix associated with I þ hj Mtjþ1 Lj;N x Kellogg and Tsan [13] the result follows. h To prove the convergence of the finite difference scheme we discretise the auxiliary problem (3.4), obtaining
8 b j ¼ I hj Mt jþ1 Lj;N u < I þ hj Mt jþ1 Lj;N U ^ jþ1 ðxji Þ; x x i : bj U 0 ¼ rT1 ð1 erðTtj Þ Þ;
1 6 i < N þ Q;
ð4:4Þ
bj U NþQ ¼ 0:
b j g be the solution of (4.4). Then, we have the following error estimates ^ j ðxÞ be the solution of (3.4) and f U Lemma 4. Let u i
^j j b j 6 CN 2 Mt jþ1 ; u xi U i
0 6 i 6 N þ Q;
where C is a constant independent of N and Mt jþ1 Proof. We remember that the local truncation error associated to the finite difference scheme is given by
j;N b j ¼ hj Mt jþ1 Lj;N Lj u ^ j ðxji Þ U ^ j ðxji Þ 6 CN2 Mt jþ1 ; u I þ hj Mt jþ1 Lx x x i
1 6 i < N þ Q:
Hence, using the discrete maximum principle for the discrete operator ðI þ hj Mtjþ1 Lj;N x Þ we have
b j j 6 CN2 Mt jþ1 ; ^ j xji U ju i
which completes the proof.
0 6 i 6 N þ Q;
h
5. The fully discrete scheme Combining the time semidiscretization scheme (3.3) with the spatial discretization scheme (4.2), the following fully discrete scheme is deduced
8 M U ¼ max 1 xM 1 6 i 6 N þ Q; > i ;0 ; > > 8i h i > > j;N j j;N jþ1 < < I þ h Mt L U ¼ I 1 h Mt L U ; 1 6 i < N þ Q ; j jþ1 x j jþ1 x i i > : j rðTt j Þ 1 > > U 0 ¼ rT 1 e ; U jNþQ ¼ 0; > > : for j ¼ M 1; . . . ; 1; 0;
ð5:1Þ
j where the discrete operators Lj;N x is described in Section 4 and U i is the fully discrete approximation to the exact solution of (2.11)–(2.14) at the mesh point ðxi ; t j Þ.
Theorem 1 (Global error). Let uðx; tÞ be the exact solution of (2.11)–(2.14) and U be the discrete solution of the fully discrete scheme (5.1). Then there exists a positive constant C, independent of N and M, such that the global error satisfies
j j 2 2 u xi ; t j U i 6 CðN þ M Þ;
0 6 i 6 N þ Q ; 0 6 j 6 M:
Proof. First, we split the global error at the time tj in the form
j j j b j U j : b j þ U ^ j ðxji Þ þ u ^ j xji U uðxi ; tj Þ U i 6 uðxi ; t j Þ u i i i Using Lemmas 1 and 4, it follows
bj 1 j j 2 2 j u xi ; t j U i 6 CM ðM þ N Þ þ U i U i :
ð5:2Þ
8673
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
b j U j can To bound the last term of (5.2), we take into account that U as the solution of one step of (5.1), taking as be written source term f ¼ 0 together with zero boundary conditions and u xji ; t jþ1 U jþ1 as final value. Then, from the uniform stai bility of the totally discrete scheme, it follows
b j U j j 6 Cjuðxj ; t jþ1 Þ U jþ1 j: jU i i i i
ð5:3Þ
The required result immediately holds from (5.2) and (5.3) and recurrence relation for the global error. h
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. We also compare our scheme with the existing popular methods in Asian option pricing. Test 1. Fixed strike Asian call Option with parameters: r ¼ 0:05; r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. Test 2. Fixed strike Asian call Option with parameters: r ¼ 0:30; r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. Here we also choose Q ¼ N=2 as in Dubois and Lelièvre [6], where the authors have indicated that the value Q ¼ N=2 is large enough to guarantee the fact that the value of u does not depend on the position of X. The exact solutions of our test problems are not available. We use the double mesh principle to estimate the errors for different N and M at t ¼ 0. We measure the accuracy in the discrete maximum norm 2N;2M eN;M ¼ maxjU N;M j; i;0 U i;0 i
and the convergence rate
RN;M ¼ log2
eN;M : e2N;2M
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 eN;M =e2N;2M is close to 4, which supports the convergence estimate of Theorem 1. The computational time for our scheme are presented in Table 3 on a PC equipped with 2.4 GHz Intel Core2 Duo CPU. The computational speed of our scheme is lower than that of Dubois and Lelièvre [6] since at each time mesh point the spline interpolation is used. But there are two advantages for our scheme. The first one is that our spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. The matrix associated with the discrete operator of the standard finite difference method for convection-dominated problems is not an M-matrix, as discussed in Dubois and
Table 1 Numerical results for Test 1. N
M
error
rate
16 32 64 128
128 256 512 1024
9.3365e-3 2.7742e-3 7.7153e-4 1.8546e-4
– 1.7508 1.8463 2.0566
N
M
error
rate
16 32 64 128
128 256 512 1024
7.1970e-4 1.7409e-4 4.3001e-5 1.0699e-5
– 2.0476 2.0174 2.0069
Table 2 Numerical results for Test 2.
Table 3 Computational times of our method for some N and M. N
M
computational time (s)
100 200 300 400
100 200 300 400
21.47 98.67 254.36 507.44
8674
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
Table 4 Errors of option value and delta value. N
our method
100 200 400 800
Dubois and Lelièvre [6]
error on price
error on delta
error on price
error on delta
0.014408 0.003511 0.000869 0.000216
0.0004863 0.0001236 0.0000323 0.0000008
0.019144 0.004697 0.001168 0.000291
0.0018231 0.0003120 0.0000603 0.0000128
Table 5 Comparisons of the prices obtained with other methods.
r
K
Our method
Zvan et al.
Vecˇerˇ
Dubois and Lelièvre
Thompson (low)
Thompson (up)
0.05
95 100 105 90 100 110
11.094095 (N = 300) 6.794315 (N = 1000) 2.744408 (N = 3000) 16.512726 (N = 300) 10.209776 (N = 1000) 5.729717 (N = 1000)
11.094 6.793 2.748 16.514 10.210 5.729
11.094 6.795 2.744 16.516 10.215 5.736
11.09409 6.7943 2.7444 16.512 10.209 5.7304
11.094094 6.794354 2.744406 16.512024 10.208724 5.728161
11.094096 6.794465 2.744581 16.523720 10.214085 5.735488
0.30
Lelièvre [6], which may lead to nonphysical oscillations in the computed solution. The second one is that at each time mesh point we obtain the spline interpolation function about the option value. We also perform a convergence analysis of the delta values. The delta values at t ¼ 0 are given by (see (2.3) and (2.6))
K K @u K delta ¼ u ;0 ;0 : S0 S0 @x S0 The computed delta values can be obtained by taking derivative of the spline interpolation function for option values. We give the absolute error of the option value and the delta value in Table 4 when r ¼ 0:05; r ¼ 0:02; T ¼ 1; K ¼ 100 and S0 ¼ 100, where the reference values from Dubois and Lelièvre [6] are price ¼ 1:697058; delta ¼ 0:6334899. From Table 4 we see that both the option values and the delta values for our method are second-order convergent and our method is more accurate than the method of Dubois and Lelièvre [6]. Finally, we give a few comparisons of the results obtained with our method and other methods, such as Zvan et al. [26], Vecˇerˇ [23], Dubois and Lelièvre [6]. The lower and upper bounds are obtained with the method of Thompson [22]. The option values for these methods are presented in Table 5 when r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. From Table 5 we see that our method is accurate both for small and large volatilities. Acknowledgement We would like to thank the anonymous referee for several suggestions for the improvement of this paper. The work was supported by National Natural Science Foundation (Grant No. 11201430) of China, Zhejiang Province Natural Science Foundation (Grant No. Y6110310) of China and Ningbo Municipal Natural Science Foundation (Grant Nos. 2012A610035, 2012A610036) of China. References [1] L. Angermann, S. Wang, Convergence of a fitted finite volume method for the penalized Black–Scholes equation governing European and American option pricing, Numerische Mathematik 106 (1) (2007) 1–40. [2] J. Barraquand, T. Pudet, Pricing of American path-dependent contingent claims, Mathematical Finance 6 (1) (1996) 17–51. [3] P. Boyle, A. Potapchik, Price sensitivities of Asian options: a survey, Insurance: Mathematics and Economics 42 (1) (2008) 189–211. [4] Z. Cen, A. Le, A robust finite difference scheme for pricing American put options with singularity-separating method, Numerical Algorithms 53 (4) (2010) 497–510. [5] Z. Cen, A. Le, A robust and accurate finite difference method for a generalized Black–Scholes equation, Journal of Computational and Applied Mathematics 235 (13) (2011) 3728–3733. [6] F. Dubois, T. Lelièvre, Efficient pricing of Asian options by the PDE approach, Journal of Computational Finance 8 (2) (2004–2005) 55–64. [7] H. Geman, M. Yor, Bessel processes Asian options and perpetuities, Mathematical Finance 3 (4) (1993) 349–375. [8] M.B. Giles, R. Carter, Convergence analysis of Crank–Nicolson and Rannacher time-marching, Journal of Computational Finance 9 (4) (2006) 89–112. [9] W.W.Y. Hsu, Y.D. Lyuu, A convergent quadratic time lattice algorithm for pricing European style Asian options, Applied Mathematics and Computation 189 (2) (2007) 1099–1123. [10] J. Hugger, A fixed strike Asian option and comments on its numerical solution, ANZIAM Journal 45 (2004) 215–231. [11] J. Ingersoll, Theory of Financial Decision Making, Roman & Littlefield, Totowa, New Jersey, 1987. [12] R. Kangro, R. Nicolaides, Far field boundary conditions for Black–Scholes equations, SIAM Journal on Numerical Analysis 38 (4) (2000) 1357–1368. [13] R.B. Kellogg, A. Tsan, Analysis of some difference approximations for a singular perturbation problem without turning points, Mathematics of Computation 32 (144) (1978) 1025–1039.
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675 [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
8675
A.G.Z. Kemna, A.C.F. Vorst, A pricing method for options based on average asset values, Journal of Banking and Finance 14 (1) (1990) 113–129. T.R. Klassen, Simple fast and flexible pricing of Asian options, Journal of Computational Finance 4 (3) (2001) 89–124. B. Lapeyre, E. Temam, Competitive Monte Carlo methods for the pricing of Asian options, Journal of Computational Finance 5 (1) (2001) 39–59. W. Mudzimbabwe, K.C. Patidar, P.J. Witbooi, A reliable numerical method to price arithmetic Asian options, Applied Mathematics and Computation 218 (22) (2012) 10934–10942. D.M. Pooley, K. Vetzal, P.A. Forsyth, Remedies for non-smooth payoffs in option pricing, Journal of Computational Finance 6 (4) (2003) 25–40. R. Rannacher, Finite element solution of diffusion problems with irregular data, Numerische Mathematik 43 (2) (1984) 309–327. L.C.G. Rogers, Z. Shi, The value of an Asian option, Journal of Applied Probability 32 (4) (1995) 1077–1088. D.Y. Tangman, A.A.I. Peer, N. Rambeerich, M. Bhuruth, Fast simplified approaches to Asian option pricing, Journal of Computational Finance 14 (4) (2011) 3–36. G.W.P. Thompson, Fast narrow bounds on the value of Asian options, Working paper, Judge Institute of Management, University of Cambridge, Cambridge, UK,1999. J. Vecˇerˇ, A new PDE approach for pricing arithmetic average Asian options, Journal of Computational Finance 4 (4) (2001) 105–113. S. Wang, A novel fitted finite volume method for the Black–Scholes equation governing option pricing, IMA Journal of Numerical Analysis 24 (4) (2004) 699–720. P. Wilmott, J. Dewynne, S. Howison, Option Pricing: Mathematical Models and Computation, Oxford Financial Press, Oxford, UK, 1993. R. Zvan, P.A. Forsyth, K. Vetzal, Robust numerical methods for PDE models of Asian options, Journal of Computational Finance 1 (2) (1998) 39–78.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Finite difference scheme with a moving mesh for pricing Asian options Zhongdi Cen ⇑, Anbo Le, Aimin Xu Institute of Mathematics, Zhejiang Wanli University, Ningbo 315100, Zhejiang, PR China
a r t i c l e
i n f o
Keywords: Asian option Path-dependent option Finite difference Moving mesh method Rannacher time stepping scheme
a b s t r a c t In this paper we propose a stable numerical method for pricing Asian call options, which is based on a central difference scheme with a moving mesh in the spatial discretization and the Rannacher time stepping scheme in the time discretization. At each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We show that the scheme is second-order convergent with respect to both time and spatial variables. Numerical results support the theoretical results. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Asian options are securities whose payoffs depend on the average value of the underlying asset price over some time period. Asian options are used by traders who are interested to hedge against the average price of an asset over a period rather than the price at the end of period. 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. This averaging feature furthermore makes Asian options enjoy lower volatilities than their underlying assets, thus as a hedge instrument they are more economical than standard options on the same underlying assets. Since the payoff of an Asian option depends on the average value of an underlying asset price over some time interval, Asian options have proven to be much more difficult to value than standard options. Standard numerical techniques tend to be impractical, inaccurate, or slow, as discussed in [3,26]. For example, traditional lattice methods [9,15] require such enormous amounts of computer memory (owing to the necessity of keeping track of every possible path throughout the tree) that they are almost impracticable. Monte Carlo simulation [14,16] works well for Asian options, but it is computationally expensive without the enhancement of variance reduction techniques. The classical numerical partial differential equation (PDE) methods are inaccurate since the degeneration of the PDE for pricing Asian options causes numerical diffusion and spurious oscillation [2,26]. Zvan et al. [26] apply a flux-limiting method from computational fluid dynamics to tackle the problem of spurious oscillations that arise in Asian options. Since their method requires that nonlinear equations should be solved, the computational workloads would be heavy and its implementation could be complicated. Hugger [10] proposes an artificial viscosity numerical method for Asian options to avoid the oscillations. Adding artificial viscosity generally smoothens out the numerical solution, but also decreases the precision since a different problem is now considered. Dubois and Lelièvre [6] extend the approach by Rogers and Shi [20] and transform the PDE for an Asian option into the Vecˇerˇ’s PDE [23]. But the new one-dimensional PDE may be convection-dom-
⇑ Corresponding author. E-mail address: [email protected] (Z. Cen). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.02.065
8668
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
inated as the standard Black–Scholes equation when the ratio of diffusion to convection term is small. It is well known that the standard finite difference method for those convection-dominated problems may lead to nonphysical oscillations in the computed solution [24]. Tangman et al. [21] use the exponential time integration scheme in combination with a dimensional splitting technique for pricing Asian options under a variety of pricing models, including the Merton model and the Carr–Geman–Madan–Yor model. However this scheme is expensive for refined meshes. Mudzimbabwe et al. [17] propose an implicit finite difference method which is stable according to a von Neumann stability analysis. Wang [24] and Angermann and Wang [1] present a fitted finite volume method for the Black–Scholes equation, but the scheme is first-order convergent. We have proposed a robust difference scheme for the Black–Scholes equation, which is based on a central difference method for spatial discretization and the implicit Euler method for time discretization [4,5]. In this paper we extend this scheme to solve the Vecˇerˇ’s PDE for Asian options. We propose a stable numerical method for pricing Asian call options, which is based on a central difference scheme with a moving mesh in the spatial discretization and the Rannacher time stepping scheme in the time discretization. Specifically, at each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We will show that the scheme is second-order convergent with respect to both time and spatial variables. The rest of the paper is organised as follows. In the next section we describe some theoretical results on the Asian option pricing model. In Section 3 we prove the convergence of the Rannacher time stepping scheme. In Section 4 we consider the spatial discretization based on a central difference scheme with a moving mesh. It is shown that the difference scheme is stable for arbitrary volatility and arbitrary interest rate, and is second-order convergent with respect to both time and spatial variables. Finally, numerical experiments are provided to support these theoretical results in Section 5. 2. The continuous problem Let S ¼ SðtÞ denote the price of the underlying asset and IðtÞ denote the underlying asset price running sum given by
IðtÞ ¼
Z
t
SðsÞds:
0
The price CðS; I; tÞ of the continuous arithmetic average Asian call option with fixed strike price satisfies the following PDE [11,25]
@C 1 2 2 @ 2 C @C @C þ r S þS rC ¼ 0; þ rS 2 @t 2 @S @I @S
ð2:1Þ
with the final condition
CðS; I; TÞ ¼ max
I K; 0 ; T
ð2:2Þ
where r is the volatility of underlying asset, r is the risk-free interest rate, T is the expiry date, and K is the strike price. Note that the above problem is a two-dimensional PDE which leads to greater computational costs. As shown by Rogers and Shi [20], by making the change of variables
y¼
K I=T ; S
CðS; I; tÞ ¼ Sv ðy; tÞ;
ð2:3Þ
the two-dimensional PDE (2.1) can be reduced to a one-dimensional PDE
@v 1 2 2 @2v 1 @v þ ry þ r y ¼ 0; 2 T @t 2 @y @y
v ðy; TÞ ¼ maxðy; 0Þ:
ð2:4Þ ð2:5Þ
However, the above one-dimensional PDE is convection-dominated when the ratio of diffusion to convection term is small, which may lead to numerical difficulty for constructing a discretization with good properties and accuracy. As shown by Dubois and Lelièvre [6], by making the change of variables
t y¼x ; T
uðx; tÞ ¼ v ðy; tÞ;
ð2:6Þ
the PDE (2.4) and (2.5) can be transformed into the Vecˇerˇ’s PDE [23]
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
uðx; TÞ ¼ maxð1 x; 0Þ:
ð2:7Þ ð2:8Þ
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
8669
The PDE (2.7) and (2.8) is defined on the whole real axis. Geman and Yor [7] have given a formula for the case I P KT as
CðS; I; tÞ ¼
S I ð1 erðTtÞ Þ þ K erðTtÞ : rT T
ð2:9Þ
By making the change of variables as in (2.3) and (2.6), we obtain
uðx; tÞ ¼
1 t rðTtÞ ð1 erðTtÞ Þ x e rT T
for x 6
t : T
ð2:10Þ
Hence we only consider the solution of the PDE (2.7) and (2.8) for x > t=T using (2.10) for the boundary condition at x ¼ t=T. Therefore, we have the following PDE
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
ðx; tÞ 2 ðt=T; 1Þ ð0; TÞ;
uðx; TÞ ¼ maxð1 x; 0Þ; x 2 ½1; 1Þ; 1 uðt=T; tÞ ¼ ð1 erðTtÞ Þ; t 2 ½0; T; rT lim uðx; tÞ ¼ 0; t 2 ½0; T:
x!1
This PDE is also convection-dominated as the standard Black–Scholes equation when the convection term is larger than the diffusion term. We propose a stable finite difference scheme with a moving mesh to solve it. For applying the numerical method we truncate the domain ðt=T; 1Þ into X ðt=T; XÞ with an appropriate choice of X. The boundary condition at x ¼ X is chosen to be uðX; tÞ ¼ 0. Normally, this truncation of the domain leads to a negligible error in the value of the option [12]. Therefore, in the remaining of this paper we will consider the following PDE:
2 @u 1 2 t @2u t @u r x ¼ 0; þ r x @t 2 T @x2 T @x
uðx; TÞ ¼ maxð1 x; 0Þ; uðt=T; tÞ ¼
x 2 ½1; 1Þ;
1 ð1 erðTtÞ Þ; rT
uðX; tÞ ¼ 0;
ðx; tÞ 2 X ð0; TÞ;
t 2 ½0; T;
t 2 ½0; T:
ð2:11Þ ð2:12Þ ð2:13Þ ð2:14Þ
Once the solution uðx; tÞ is obtained the price of the Asian option is determined by CðS; I; tÞ ¼ Suðx; tÞ. 3. The time semidiscretization As discussed in Giles and Carter [8] and Pooley et al. [18], the Crank–Nicolson method leads to serious degradation in the convergence rates of numerical schemes due to the non-smooth payoff function. The Rannacher time stepping scheme [19] uses a few first time steps with the implicit Euler method, and after that it uses the Crank–Nicolson method. They show that the Rannacher time stepping is simple to implement and greatly improves the stability of the numerical scheme in approximations of the Black–Scholes equations. We also use the Rannacher time stepping scheme to approximate the solution of (2.11)–(2.14). To ensure that the scheme M on the time interval ½0; T: is second-order accurate we use a piecewise uniform mesh x
( tj ¼
T4M 2 j; M4 T Mj ; M2
j ¼ 0; 1; . . . ; M 4; j ¼ M 3; . . . ; M:
ð3:1Þ
It is easy to see that the mesh sizes Mt j ¼ tj tj1 satisfies
( Mt j ¼
T4M2 M4 1 M2
j ¼ 0; 1; . . . ; M 4; j ¼ M 3; . . . ; M:
ð3:2Þ
Then the Rannacher time stepping scheme to approximate the solution of (2.11)–(2.14) reads
8 M u ¼ uðx; TÞ ¼ maxð1 x; 0Þ; > > > > h i > 8 > > < < I þ hj Mtjþ1 Ljx uj ¼ I ð1 hj ÞMt jþ1 Ljþ1 ujþ1 ; x > : j > > u ðt j =TÞ ¼ rT1 ð1 erðTtj Þ Þ; > > > > : for j ¼ M 1; . . . ; 1; 0;
uj ðXÞ ¼ 0;
ð3:3Þ
8670
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
where
Ljx uj ¼ ( hj ¼
2 2 j 1 2 t @ u t j @uj r x j þ r x ; 2 T @x2 T @x
1; j ¼ M 1; . . . ; M 4; 1 ; 2
j ¼ M 5; . . . ; 1; 0;
and uj ðxÞ denotes the approximation of the exact solution uðx; tÞ at the time level t j . Similarly to Kellogg and Tsan [13] we can prove that the differential operator ðI þ hj Mt jþ1 Ljx Þ satisfies a maximum principle, which we will use in the analysis of the convergence of the Rannacher time stepping scheme. Estimates for the global error are deduced from appropriate bounds of the local error, where the auxiliary problem
8 jþ1 ^ ¼ uðx; t jþ1 Þ; u > > > > h i > 8 > > < < I þ hj Mtjþ1 Ljx u ^ j ¼ I ð1 hj ÞMtjþ1 Ljþ1 ^ jþ1 ; u x > : j > > ^ ðt j =TÞ ¼ rT1 1 erðTtj Þ ; u > > > > : for j ¼ M 1; . . . ; 1; 0;
ð3:4Þ
^ j ðXÞ ¼ 0; u
is introduced to define the local error. ^ j ðxÞ satisfies Lemma 1 (Local error estimate). The local error associated to the method (3.1), defined by ej ¼ uðx; t j Þ u
kej kX ¼ OðM 3 Þ;
0 6 j 6 M;
ð3:5Þ
. where kej kX is the maximum norm of ej on the closed set X Proof. Using Taylor expansion we have
uðx; tjþ1 Þ uðx; tj Þ ¼ ut ðx; tj Þ þ OðMt jþ1 Þ ¼ Ljx uðx; tj Þ þ OðMt jþ1 Þ; Mt jþ1
ð3:6Þ
for j ¼ M 1; . . . ; M 4, and
uðx; tjþ1 Þ uðx; tj Þ ut ðx; t j Þ þ ut ðx; tjþ1 Þ ¼ ut ðx; tjþ1=2 Þ þ OððMt jþ1 Þ2 Þ ¼ þ OððMt jþ1 Þ2 Þ Mt jþ1 2 ¼
Ljx uðx; tj Þ þ Ljþ1 x uðx; t jþ1 Þ þ OððMt jþ1 Þ2 Þ 2
for j ¼ M 5; . . . ; 1; 0, where tjþ1=2 ¼ error is the solution of the problem
8 < I þ h Mt j
j jþ1 Lx
t j þt jþ1 . 2
ð3:7Þ
From (3.2)–(3.4) and (3.6) and (3.7), it is straightforward to show that the local
ej ¼ OðM 3 Þ;
: ej ðt =TÞ ¼ ej ðXÞ ¼ 0; j
and therefore the result follows from the maximum principle for the operator ðI þ hj Mtjþ1 Ljx Þ. h Lemma 2 (Global error estimate). The global error associated to the Rannacher time stepping scheme, given by Ej ¼ uðx; tj Þ uj ðxÞ, satisfies
E ¼ supkEj k1 6 CM 2 ;
ð3:8Þ
j6M
and therefore the time semidiscretization method is a second-order convergent scheme. Proof. Using the local error estimate down to the jth time step given by Lemma 1, we get the following global error estimate at the jth time step:
X M j kE k1 ¼ e 6 kej k1 þ kejþ1 k1 þ . . . þ keM k1 6 CM 2 ; k¼j j
1
for 0 6 j 6 M, where C is a positive constant independent of M. h
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
8671
4. The spatial discretization As discussed in introduction, the standard finite difference method to discretise the spatial derivatives of the problem (3.3), numerical difficulty can be caused. The main reason is that when the ratio of diffusion to convection term is small, the differential operator becomes convection-dominated. Since the coefficients of the derivatives in (3.3) are time-dependent, we apply the central difference scheme on a moving mesh to discretise problem (3.3). Specifically, at each time mesh point we make a piecewise uniform mesh to discretise the space interval, which ensures that the matrix associated with the discrete operator is an M-matrix. Hence the spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. We will show that the scheme is second-order convergent with respect to the spatial variable. At the time mesh point t j , we use a piecewise uniform mesh XNj on the spatial interval ½tj =T; X:
8 tj j þh; i ¼ 1; > >
ð4:1Þ
where j
h ¼
1 tj =T ; 2 1 þ rr ðN 1Þ
X ¼1þ
Q : N j
It is easy to see that the mesh sizes hi ¼ xji xji1 satisfy
8 j i ¼ 1; >
Thus, at each time mesh point tj we apply the central difference scheme on the above piecewise uniform mesh to approximate problem (3.3):
8 h i < I þ hj Mt jþ1 Lj;N U j ¼ I ð1 hj ÞMtjþ1 Ljþ1;N U jþ1 ; x x i i : U j ¼ 1 ð1 erðTtj Þ Þ; 0 rT
U jNþQ
1 6 i < N þ Q;
ð4:2Þ
¼ 0;
where
2 tj r2 j j Lj;N U ¼ x x i i j T h þ hj i
iþ1
U jiþ1 U ji j
hiþ1
U ji U ji1
!
j
hi
j j t j U iþ1 U i1 þ r xji ; j j T h þh i
ð4:3Þ
iþ1
and U jþ1 has been obtained by using the spline interpolation of the numerical results from the time mesh point tjþ1 . i N Lemma 3 (Discrete maximum principle). The operator ðI þ hj Mtjþ1 Lj;N x Þ defined by (4.3) on the piecewise uniform mesh Xj satisfies a discrete maximum principle, i.e. if v i and wi are mesh functions that satisfy v 0 P w0 ; v NþQ P wNþQ , and j;N ðI þ hj Mt jþ1 Lj;N x Þv i P ðI þ hj Mt jþ1 Lx Þwi ð1 6 i < N þ Q Þ, then v i P wi for all i.
Proof. Let
ai ¼
2 tj rhj Mtjþ1 tj j j x x ; i i j j j T T Þh h þh
r2 hj Mtjþ1 j
j
i
iþ1
ðhi þ hiþ1 i i iþ1 2 2 r hj Mtjþ1 j tj xi þ 1; 1 6 i < N þ Q; bi ¼ j j T hi hiþ1 2 r2 hj Mtjþ1 tj rhj Mtjþ1 tj j j : ci ¼ j x þ x i i j j j j T T ðh þ h Þh h þh iþ1
i
iþ1
We can obtain
r2 hj Mtjþ1 tj tj rhj Mt jþ1 tj hj Mt jþ1 tj tj j ci 6 xj1 xji r2 xj1 xji þ j ¼ þ rhiþ1 xji j j j j j j j T T T T T hi þ hiþ1 hi þ hiþ1 hiþ1 hi þ hiþ1 hiþ1 hj Mt jþ1 r2 j tj j ¼ 0; 1 6 i < N þ Q: xji r2 h þ r h ¼ j j j r T ðhi þ hiþ1 Þhiþ1
8672
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
Clearly,
ai < 0 for 2 6 i 6 N þ Q 1;
bi > 0 for 1 6 i 6 N þ Q 1;
and
b1 þ c1 > 0; ai þ bi þ ci > 0; 2 6 i 6 N þ Q 2; aNþQ 1 þ bNþQ1 > 0:
is an M-matrix. Thus by the same argument as Lemma 3.1 in Hence we verify that the matrix associated with I þ hj Mtjþ1 Lj;N x Kellogg and Tsan [13] the result follows. h To prove the convergence of the finite difference scheme we discretise the auxiliary problem (3.4), obtaining
8 b j ¼ I hj Mt jþ1 Lj;N u < I þ hj Mt jþ1 Lj;N U ^ jþ1 ðxji Þ; x x i : bj U 0 ¼ rT1 ð1 erðTtj Þ Þ;
1 6 i < N þ Q;
ð4:4Þ
bj U NþQ ¼ 0:
b j g be the solution of (4.4). Then, we have the following error estimates ^ j ðxÞ be the solution of (3.4) and f U Lemma 4. Let u i
^j j b j 6 CN 2 Mt jþ1 ; u xi U i
0 6 i 6 N þ Q;
where C is a constant independent of N and Mt jþ1 Proof. We remember that the local truncation error associated to the finite difference scheme is given by
j;N b j ¼ hj Mt jþ1 Lj;N Lj u ^ j ðxji Þ U ^ j ðxji Þ 6 CN2 Mt jþ1 ; u I þ hj Mt jþ1 Lx x x i
1 6 i < N þ Q:
Hence, using the discrete maximum principle for the discrete operator ðI þ hj Mtjþ1 Lj;N x Þ we have
b j j 6 CN2 Mt jþ1 ; ^ j xji U ju i
which completes the proof.
0 6 i 6 N þ Q;
h
5. The fully discrete scheme Combining the time semidiscretization scheme (3.3) with the spatial discretization scheme (4.2), the following fully discrete scheme is deduced
8 M U ¼ max 1 xM 1 6 i 6 N þ Q; > i ;0 ; > > 8i h i > > j;N j j;N jþ1 < < I þ h Mt L U ¼ I 1 h Mt L U ; 1 6 i < N þ Q ; j jþ1 x j jþ1 x i i > : j rðTt j Þ 1 > > U 0 ¼ rT 1 e ; U jNþQ ¼ 0; > > : for j ¼ M 1; . . . ; 1; 0;
ð5:1Þ
j where the discrete operators Lj;N x is described in Section 4 and U i is the fully discrete approximation to the exact solution of (2.11)–(2.14) at the mesh point ðxi ; t j Þ.
Theorem 1 (Global error). Let uðx; tÞ be the exact solution of (2.11)–(2.14) and U be the discrete solution of the fully discrete scheme (5.1). Then there exists a positive constant C, independent of N and M, such that the global error satisfies
j j 2 2 u xi ; t j U i 6 CðN þ M Þ;
0 6 i 6 N þ Q ; 0 6 j 6 M:
Proof. First, we split the global error at the time tj in the form
j j j b j U j : b j þ U ^ j ðxji Þ þ u ^ j xji U uðxi ; tj Þ U i 6 uðxi ; t j Þ u i i i Using Lemmas 1 and 4, it follows
bj 1 j j 2 2 j u xi ; t j U i 6 CM ðM þ N Þ þ U i U i :
ð5:2Þ
8673
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
b j U j can To bound the last term of (5.2), we take into account that U as the solution of one step of (5.1), taking as be written source term f ¼ 0 together with zero boundary conditions and u xji ; t jþ1 U jþ1 as final value. Then, from the uniform stai bility of the totally discrete scheme, it follows
b j U j j 6 Cjuðxj ; t jþ1 Þ U jþ1 j: jU i i i i
ð5:3Þ
The required result immediately holds from (5.2) and (5.3) and recurrence relation for the global error. h
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. We also compare our scheme with the existing popular methods in Asian option pricing. Test 1. Fixed strike Asian call Option with parameters: r ¼ 0:05; r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. Test 2. Fixed strike Asian call Option with parameters: r ¼ 0:30; r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. Here we also choose Q ¼ N=2 as in Dubois and Lelièvre [6], where the authors have indicated that the value Q ¼ N=2 is large enough to guarantee the fact that the value of u does not depend on the position of X. The exact solutions of our test problems are not available. We use the double mesh principle to estimate the errors for different N and M at t ¼ 0. We measure the accuracy in the discrete maximum norm 2N;2M eN;M ¼ maxjU N;M j; i;0 U i;0 i
and the convergence rate
RN;M ¼ log2
eN;M : e2N;2M
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 eN;M =e2N;2M is close to 4, which supports the convergence estimate of Theorem 1. The computational time for our scheme are presented in Table 3 on a PC equipped with 2.4 GHz Intel Core2 Duo CPU. The computational speed of our scheme is lower than that of Dubois and Lelièvre [6] since at each time mesh point the spline interpolation is used. But there are two advantages for our scheme. The first one is that our spatial discretization scheme is maximum-norm stable for arbitrary volatility and arbitrary interest rate. The matrix associated with the discrete operator of the standard finite difference method for convection-dominated problems is not an M-matrix, as discussed in Dubois and
Table 1 Numerical results for Test 1. N
M
error
rate
16 32 64 128
128 256 512 1024
9.3365e-3 2.7742e-3 7.7153e-4 1.8546e-4
– 1.7508 1.8463 2.0566
N
M
error
rate
16 32 64 128
128 256 512 1024
7.1970e-4 1.7409e-4 4.3001e-5 1.0699e-5
– 2.0476 2.0174 2.0069
Table 2 Numerical results for Test 2.
Table 3 Computational times of our method for some N and M. N
M
computational time (s)
100 200 300 400
100 200 300 400
21.47 98.67 254.36 507.44
8674
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675
Table 4 Errors of option value and delta value. N
our method
100 200 400 800
Dubois and Lelièvre [6]
error on price
error on delta
error on price
error on delta
0.014408 0.003511 0.000869 0.000216
0.0004863 0.0001236 0.0000323 0.0000008
0.019144 0.004697 0.001168 0.000291
0.0018231 0.0003120 0.0000603 0.0000128
Table 5 Comparisons of the prices obtained with other methods.
r
K
Our method
Zvan et al.
Vecˇerˇ
Dubois and Lelièvre
Thompson (low)
Thompson (up)
0.05
95 100 105 90 100 110
11.094095 (N = 300) 6.794315 (N = 1000) 2.744408 (N = 3000) 16.512726 (N = 300) 10.209776 (N = 1000) 5.729717 (N = 1000)
11.094 6.793 2.748 16.514 10.210 5.729
11.094 6.795 2.744 16.516 10.215 5.736
11.09409 6.7943 2.7444 16.512 10.209 5.7304
11.094094 6.794354 2.744406 16.512024 10.208724 5.728161
11.094096 6.794465 2.744581 16.523720 10.214085 5.735488
0.30
Lelièvre [6], which may lead to nonphysical oscillations in the computed solution. The second one is that at each time mesh point we obtain the spline interpolation function about the option value. We also perform a convergence analysis of the delta values. The delta values at t ¼ 0 are given by (see (2.3) and (2.6))
K K @u K delta ¼ u ;0 ;0 : S0 S0 @x S0 The computed delta values can be obtained by taking derivative of the spline interpolation function for option values. We give the absolute error of the option value and the delta value in Table 4 when r ¼ 0:05; r ¼ 0:02; T ¼ 1; K ¼ 100 and S0 ¼ 100, where the reference values from Dubois and Lelièvre [6] are price ¼ 1:697058; delta ¼ 0:6334899. From Table 4 we see that both the option values and the delta values for our method are second-order convergent and our method is more accurate than the method of Dubois and Lelièvre [6]. Finally, we give a few comparisons of the results obtained with our method and other methods, such as Zvan et al. [26], Vecˇerˇ [23], Dubois and Lelièvre [6]. The lower and upper bounds are obtained with the method of Thompson [22]. The option values for these methods are presented in Table 5 when r ¼ 0:15; T ¼ 1; K ¼ 100 and S0 ¼ 100. From Table 5 we see that our method is accurate both for small and large volatilities. Acknowledgement We would like to thank the anonymous referee for several suggestions for the improvement of this paper. The work was supported by National Natural Science Foundation (Grant No. 11201430) of China, Zhejiang Province Natural Science Foundation (Grant No. Y6110310) of China and Ningbo Municipal Natural Science Foundation (Grant Nos. 2012A610035, 2012A610036) of China. References [1] L. Angermann, S. Wang, Convergence of a fitted finite volume method for the penalized Black–Scholes equation governing European and American option pricing, Numerische Mathematik 106 (1) (2007) 1–40. [2] J. Barraquand, T. Pudet, Pricing of American path-dependent contingent claims, Mathematical Finance 6 (1) (1996) 17–51. [3] P. Boyle, A. Potapchik, Price sensitivities of Asian options: a survey, Insurance: Mathematics and Economics 42 (1) (2008) 189–211. [4] Z. Cen, A. Le, A robust finite difference scheme for pricing American put options with singularity-separating method, Numerical Algorithms 53 (4) (2010) 497–510. [5] Z. Cen, A. Le, A robust and accurate finite difference method for a generalized Black–Scholes equation, Journal of Computational and Applied Mathematics 235 (13) (2011) 3728–3733. [6] F. Dubois, T. Lelièvre, Efficient pricing of Asian options by the PDE approach, Journal of Computational Finance 8 (2) (2004–2005) 55–64. [7] H. Geman, M. Yor, Bessel processes Asian options and perpetuities, Mathematical Finance 3 (4) (1993) 349–375. [8] M.B. Giles, R. Carter, Convergence analysis of Crank–Nicolson and Rannacher time-marching, Journal of Computational Finance 9 (4) (2006) 89–112. [9] W.W.Y. Hsu, Y.D. Lyuu, A convergent quadratic time lattice algorithm for pricing European style Asian options, Applied Mathematics and Computation 189 (2) (2007) 1099–1123. [10] J. Hugger, A fixed strike Asian option and comments on its numerical solution, ANZIAM Journal 45 (2004) 215–231. [11] J. Ingersoll, Theory of Financial Decision Making, Roman & Littlefield, Totowa, New Jersey, 1987. [12] R. Kangro, R. Nicolaides, Far field boundary conditions for Black–Scholes equations, SIAM Journal on Numerical Analysis 38 (4) (2000) 1357–1368. [13] R.B. Kellogg, A. Tsan, Analysis of some difference approximations for a singular perturbation problem without turning points, Mathematics of Computation 32 (144) (1978) 1025–1039.
Z. Cen et al. / Applied Mathematics and Computation 219 (2013) 8667–8675 [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]
8675
A.G.Z. Kemna, A.C.F. Vorst, A pricing method for options based on average asset values, Journal of Banking and Finance 14 (1) (1990) 113–129. T.R. Klassen, Simple fast and flexible pricing of Asian options, Journal of Computational Finance 4 (3) (2001) 89–124. B. Lapeyre, E. Temam, Competitive Monte Carlo methods for the pricing of Asian options, Journal of Computational Finance 5 (1) (2001) 39–59. W. Mudzimbabwe, K.C. Patidar, P.J. Witbooi, A reliable numerical method to price arithmetic Asian options, Applied Mathematics and Computation 218 (22) (2012) 10934–10942. D.M. Pooley, K. Vetzal, P.A. Forsyth, Remedies for non-smooth payoffs in option pricing, Journal of Computational Finance 6 (4) (2003) 25–40. R. Rannacher, Finite element solution of diffusion problems with irregular data, Numerische Mathematik 43 (2) (1984) 309–327. L.C.G. Rogers, Z. Shi, The value of an Asian option, Journal of Applied Probability 32 (4) (1995) 1077–1088. D.Y. Tangman, A.A.I. Peer, N. Rambeerich, M. Bhuruth, Fast simplified approaches to Asian option pricing, Journal of Computational Finance 14 (4) (2011) 3–36. G.W.P. Thompson, Fast narrow bounds on the value of Asian options, Working paper, Judge Institute of Management, University of Cambridge, Cambridge, UK,1999. J. Vecˇerˇ, A new PDE approach for pricing arithmetic average Asian options, Journal of Computational Finance 4 (4) (2001) 105–113. S. Wang, A novel fitted finite volume method for the Black–Scholes equation governing option pricing, IMA Journal of Numerical Analysis 24 (4) (2004) 699–720. P. Wilmott, J. Dewynne, S. Howison, Option Pricing: Mathematical Models and Computation, Oxford Financial Press, Oxford, UK, 1993. R. Zvan, P.A. Forsyth, K. Vetzal, Robust numerical methods for PDE models of Asian options, Journal of Computational Finance 1 (2) (1998) 39–78.