Version pdf de ce document

Valuation of Asian Options using moments

Gianluca Fusai*
Aldo Tagliani


*Dipartimento SEMEQ, Università del Piemonte Orientale, Via Lanino 1, 28100 Novara, Italia, e-mail: gianluca.fusai@eco.unipmn.it.
DISA, Università di Trento, Via Inama, Trento, Italia.

Contents

 1 Introduction
 2 The general framework
 3 The Laplace transform (Geman-Yor)
 4 The upper and lower bounds for the price (Thompson, Rogers and Shi)
 5 The moment problem
  5.1 Approximations based on the moments of AT _
   5.1.1 The lognormal approximation (Levy)
   5.1.2 The Edgeworth series approximation (Turnbull and Wakeman)
   5.1.3 The reciprocal gamma approximation (Milevsky and Posner)
  5.2 Approximations based on the moments of -1
At.
  5.3 Approximations based on the moments of lnAT _
   5.3.1 The Edgeworth approximation (Fusai and Tagliani)
 6 A comparison between the different approximations
 7 Conclusion

Premia 14

1 Introduction

Asian options, options for which the payoff depends on the arithmetic average value of the asset price over some time period, have had a very large success in the last years, because they reduce the possibility of market manipulation near the expiry date and offer a better hedge to firms having a stream of positions.

Actually there is no closed form solution for the price of this option because in the usual Black-Scholes framework, the arithmetic average is a sum of correlated lognormal distributions for which there is no recognizable density function. Several tentatives have been done:

In the present paper we examine in greater detail the methods based on the moments, i.e. assuming as unique information the moments of the average propose different approximations. The underlying framework is always a Black-Scholes world where the underlying asset evolves following a Geometric Brownian Motion process. We can mainly distinguish the following approaches:

i) using the first two algebraic moments of the average, fit a lognormal density, Levy [21].;

ii) using the first four algebraic moments of the average, approximate the density using an Edgeworth series expansion around a lognormal density, Turnbull and Wakeman [33] and Ritchken and al. [27];

iii) using the first two algebraic moments of the average fit a reciprocal gamma density, Milevsky and Posner [25];

iv) using the first four logarithmic moments of the average, approximate the density using an Edgeworth series expansion around a normal density as in Fusai and Tagliani[14],

v) always using the logarithmic moments of the average, Fusai and Tagliani[14] propose a Maximum Entropy approximation.

vi) using the moments of the reciprocal of the average, Dufresne approximates the density using a Laguerre series expansion.

In using different approximations based on the moments, we would like that:

a) the true density is uniquely determined by the infinite sequence of moments (the so called moments problem);

b) for a given number of moments, the approximated distribution should be positive,

c) for a given number of moments, a bound to the error made in pricing Asian options is provided;

d) increasing the number of moments considered, the sequence of approximants is converging to the true distribution in some norm and the same happens for the Asian option price.

We remark that as proved in Fusai and Tagliani[14] only the Maximum Entropy approximation satisfies all the above requirements. Property a) is satisfied if we use logarithmic moments [14] or the moments of the reciprocal, Dufresne [12], whilst no result is available about arithmetic moments. Property b) is usually not guaranteed when we approximate the true density with an Edgeworth or Laguerre series expansion. Property c) has never been discussed in the literature except in Fusai and Tagliani[14] that give a bound to the error in terms of the number of moments considered. Property d) is a delicate question when we use Edgeworth series, indeed in practical work, the Edgeworth series converges for a small number of terms and then diverges.

Finally, in order to have a reference point for comparing the different approximations we use:

  1. the lower and upper bounds given in Rogers and Shi [28] and in Thompson [31].
  2. the numerical inversion of the Laplace transform given in Geman and Yor, [17],
  3. the numerical solution of the reduced PDE given in Rogers and Shi [28].

We do not consider MonteCarlo simulation because, although there are a number of techniques for variance reduction that allow to limit the number of simulations, it is still quite time consuming requiring in every case several thousands of simulations.

In section 1 we describe the Black-Scholes framework underlying the Asian option problem. In section 2 we discuss Gemany-Yor Laplace transform approach and the related numerical inversion using the Dubner Abate method as discussed in Abate and Whitt, [1]. In section 3 we present the lower and upper bounds given in Rogers and Shi [28] and in Thompson [31]. In section 4 we discuss the general moment problem recalling some basic fact useful when we want approximate a distribution using its moments. Then we describe different approximations based on the moments. We distinguish approaches based on a) algebraic moments of the arithmetic average, b) moments of the reciprocal of arithmetic average, c) moments of the logarithm of the arithmetic average and illustrate the main properties of the different methods. Finally in section 6, we compare the different approximations and then we draw some conclusion. In Appendix A we give some more details about the derivation of the Geman-Yor formula.

2 The general framework

We assume that under the martingale measure the underlying asset price evolves according to a Geometric Brownian Motion:

dS  = rS dt+  σS dW
   t     t      t   t

        (    )
         r- σ22 t+σWt
St = S0e
(1)

where r is the instantaneous risk free rate, σ is the instantaneous volatility and Wt is a Wiener process. If the asset pays a continuous dividend yield q the risk free rate becomes r - q. Assuming that the average is computed continuously in time in the period (0,t), we define:

     ∫ t
At =    e((r-σ2∕2)w+ σWw )dw
      0
(2)

and the payoff at time t of a Asian option is given by max[AtS0∕t - K,0] for a call option and by max[K -AtS0∕t,0] for a put option. For numerical reasons, it is convenient to consider the case of an Asian put option for which the price is computed as:

               ∫ K
E0 [K  - At]+ =    (K - xS0 ∕t)fAt(x;t)dx
                0
(3)

where fAt denotes the density of At. Then the pricing problem is essentially related to the determination of this density.

3 The Laplace transform (Geman-Yor)

In this section we briefly review the analytical result due to Geman and Yor, [17]. Using the fact that a Geometric Brownian Motion is a time-changed squared Bessel process and the stability by additivity of this process, Geman and Yor [17] have been able to provide the Laplace transform of the Asian call option1 . If we define:

      e-rt4S
cgy = -----2-c(h,q)
       t  σ
(4)

the Laplace transform wrto h =σ2t∕4 of the function c(h, q) is given by:

          ∫ +∞                 ∫1∕(2q) -x  1(μ- v)-2         12(μ+v)+1
C (λ,q) =      e-λhc (h, q)dh = -0----e--x-2-------(1(---2qx)-----)--dx-
           0                       λ (λ- 2 - 2v)Γ  12 (μ- v) - 1
(5)

where:

    -2r        σ2t K-     ∘ ------2
v = σ2 - 1;q =  4  S ;μ =   2λ + v

The same result has been obtained recently by Lipton [24] reducing the augmented Black-Scholes PDE through the use of a non-dimensional variable and then solving it using the Laplace transform.

Although there are numerous and well understood techniques for the Laplace Transform inversion, in Finance this technique has received a very limited attention, mainly because it is reputed, against the evidence, to be a method slow and difficult to implement accurately. In the case of Asian option we are aware only of a work by Fu and al. [13], and among the different routines tested we agree with them in using the algorithm due to Dubner and Abate [10] and described in detail in Abate and Whitt [1].

The inversion algorithm is essentially a trapezoidal rule approximation to the Bromwich inversion formula. Since the trapezoidal rule is a quite simple integration procedure, its use can appear surprising. It turns out to be surprisingly effective in this context with periodic and oscillating integrands, because the errors tend to cancel. In particular, it turns out to be better than familiar alternatives such as Simpson’s rule for inversion integrals. Two features of the method are: a) the expression for the maximum error in the computed inverse transform and b) the use of the Euler summation algorithm in order to exploit the alternating character of the series and accelerate its convergence, for details compare Abate and Whitt [1]. The inversion formula is then given by:

             (   )
          ∑m   m   - m
cda (h,q) ≈      k  2   sn+k (h )
          k=0

where:

         A∕2   (   (     ))    A ∕2 ∑n          (  (           ) )
sn(h) = e---Re   C  -A-,q   + e---    (- 1)k Re  C   A-+-2kπi,q
         2h         2h          h  k=1                 2h
(6)

The algorithm requires the specification of three parameters: A, n and m: a good choice results to be A = 18.4, m= 15, n= 11, except when σ = 0.1 where we have used m= 35, n= 15. The algorithm has been implemented in Microsoft Visual C++ Version 5.0 and all the calculations have been done on a Compaq Presario with Pentium 133 processor. The integral in (5) has been computed using Gaussian quadrature (Legendre) with 100 nodes. In table 1 we report the results of the inversion and the lower and upper bound due to Thompson [31] and to Rogers and Shi [28]. More results can be found in tables 5a-5c.













Table1: Laplace Inversion and Bounds
on the Asian Option Price






σ√-
 t rt K Lower AW Upper





0.1 0.09100 4.91508 4.91512 4.91541 0.2 0.09100 6.77700 6.77735 6.77866 0.3 0.09100 8.82755 8.82876 8.83329 0.4 0.0910010.9209010.92377 10.93596 0.5 0.0910013.0225313.02816 13.05680











The inversion did not present problems except for low volatility levels (σ√-
 t < 0.08) when we had numerical problems in the computation of the integral in (5), due to finite precision of computers (double precision provides 16 decimal places). The same problem has been noticed in solving the related PDE with finite difference schemes. Viceversa, in this case Monte Carlo simulation results to be very effective.

Geman and Eydeland [15] suggest to replace e-z∕2q in (5) using a series expansion and the relationship between beta and gamma functions. Then they obtain the following approximate expression for the Laplace transform2 :

                  (    )
-----1-------M∑  1-   -1- 1+δ+i Γ-((μ---v)-∕2--1-+-i)-Γ ((μ-+-v)∕2-+-2)
λ(λ - 2v - 2)   i! - 2q           Γ (μ+ 1 + i)Γ ((μ - v)∕2 - 1)
             i=0

where M is the number of terms used in the series expansion. For the numerical inversion they apply the Fast Fourier Transform. However, we have applied to this expression the AW algorithm described above and the inversion was accurate if σ√ --
  T > 0.3 with M=40. If σ√--
 T < 0.3 the numerical inversion results to be completely inaccurate, independently on the number of terms in the series expansion.

4 The upper and lower bounds for the price (Thompson, Rogers and Shi)

A different approach to the problem of Asian option price is to find a lower and upper bound for the price and practically they result to be surprisingly accurate, i.e. the size of the interval results to be small and in general the Asian option price will be near the lower bound. The idea of finding the bounds is due to Curran [7], [8] and to Rogers and Shi [28]. These bounds have been improved by Thompson [31], [32] that was able to simplify the computation of the lower bound given in Rogers and Shi. Moreover, Thompson gives an upper bound for the option price, much more stronger than the one obtained by Rogers and Shi. Using a lower and upper bound can give us some better feeling in choosing the different methods described above and then understand the limits of the different approaches.

The Asian option price is obtained through the computation of the quantity

  ( S ∫ t (r-σ2∕2)s+σWs       )+
E   t- 0 e            ds - K    .

Making the substitution h = s∕t and using the scaling property of the BM:

  ∫ t                     ∫ 1
1-   e(r-σ2∕2)s+σWsds  =      e(r-σ2∕2)Th+σWThdh
 t 0                       0
                          ∫ 1 (r-σ2∕2)Th+σ√T-Wh
                      =    0 e                dh
we can always consider the case with time to expiry equal to 1 and set the interest rate equal to rt and the volatility equal to σ√ -
  t. So in the following we will assume that the residual life of the option is equal to 1 and the risk-free rate r stands for rt and the volatility σ stands for σ√ -
  t. Moreover we set:
     ∫ 1 (r- σ2∕2)th+ σ√tWh
X  =  0 e               dh - K

and the option price is given by EX+.

Using the iterated conditional expectation and the fact that X+ X and that X+ is a positive quantity, then:

  (    )     [ (       )]    [          ]
E  X+   =  E E   X+ |Z   ≥ E  E (X |Z )+

and the accuracy of the lower bound can be estimated from the quite general bounds:

     (    )    [          ]      [∘ ----------]
0 ≤ E  X+  - E  E (X |Z )+ ≤  1E    V ar(X |Z )
                              2

The method can be very successfully if the conditioning variable Z is chosen to make the conditional variance small. Rogers and Shi realize that a good choice results to be Z = 0tWsds and then they are able to give an expression to E[         ]
 E (X |Z )+ as a double integral:

             ∫ +∞      (    )[∫ 1                                       ]+
c ≥ clow ≡ e-r      √3ϕ  √3z-      S0e3σξ(1-ξ∕2)z+αξ+12σ2(ξ-3ξ2(1-ξ∕2)2)dξ - K   dz
              - ∞              0

where α = (r - σ2∕2), and where ϕ(x) is the normal density.

Thompson gives a simpler version of the above formula, trying to approximate the event that the option makes a positive payout with a simpler event. If A = { 01Stdt > K}, we have for any event A*:

  [∫ 1         ]+     [( ∫ 1         )    ]   ∫ t
E     Sudu - K    ≥ E       Sudu - K   1A* =     E (Su - K; A *)du
    0                     0                    0
(7)

and a natural choice is A* = { 0tWudu > γ} where γ is chosen in an optimal way maximinizing the rhs in (7). The optimal value of γ, γ*, results to be the unique solution to:

∫        (                          (                ))
  1S  exp  3γ*σξ(1 - ξ∕2)+ α ξ + 1-σ2 ξ - 3ξ2(1- ξ∕2)2   dξ = K
 0  0                           2
(8)

and the lower bound is:

         (                     )
    ∫ 1
e-r    E (Su - K; 1(∫1       *)) dt
     0               0 Wudu>γ

i.e.:

          ∫ 1            (   *              )         (    * )
clow = e- r   S0eαξ+12σ2ξΦ   - γ-+-σξ(√1---ξ∕2)- dξ - K Φ  -- γ√--
           0                     1∕  3                  1∕  3
(9)

where Φ is the normal distribution function.

Thompson shows also the lower bound is the same as the one given in Rogers and Shi with the advantage of being expressed as a single integral. The computation of the optimal value of γ can be done very quickly with standard optimization routines, such as bisection or Newton search.

The upper bound given in Thompson is:

          - r∫ 1∫ +∞         [        (a (t,x ))          ( a(t,x)) ]
c ≤ cup ≡ e          2vφ (w)  a(t,x)Φ  -b(t,x)  +  b(t,x) φ  b(t,x)   dwdv
              0  -∞

where v = √ -
  t, w = x∕√ -
  t,

a(t,x)  =  S0 e∘xp(σx-+-αt)---K-(μt + σx )+ K σ(1 - t∕2)x
                 1-           2
b(t,x)  =  K σ   3 - t(1 - t∕2)
            1               √ --
    μt  =  K- (S0exp (αt)+ γ  vt)
            2                            2
    vt  =  ctt+ 2 (Kσ )ctt (1 - t∕2)+ (K σ) ∕3
    ct  =  S0 exp(αt)σ - K σ
                     α          ∫ 1√ --
     γ  =  (K -  S0(e -  1)∕α)∕      vtdt
                                 0

This upper bound results to be tighter than the one given in Rogers and Shi [28].

5 The moment problem

Before describing in detail the different approaches that make use of the moments, we review some questions related to the determination of the distribution of the average using the moments technique.

1. There is some initial reason to think that matching moments will give a good approximation. For example, Lindsay and Roeder [23] show that if two distributions functions have the same first n moments, then they must cross each other at least n times. Moreover, Akhiezer [2], page 66, shows that the possible error can be bounded, and establishes the relevant relationships of the difference between two distributions that share the same 2m moments. However more recently, Lindsay and Basak [22] show that this bound results to be quite strict only in the tails of the distribution.

2. When the random variable assumes only positive values, the distribution recovering from its moments (the so called Stieltjes moment problem) arises both existence and determinacy questions. We consider uniquely the latter one. In Stieltjes moment problem an infinite sequence of moments determines a unique distribution iff this distribution decays asymptotically as an exponential exp(-αxλ), α > 0, λ 12. For example, this does not happen in the case of the lognormal density, that is, there are multiple distributions with exactly the same moments3 . This distribution has indeed tails decaying slower than exp(-αxλ).

3. In the case the given information is represented uniquely by the sequence of moments, there are sufficient conditions for the moment problem to be determinate or indeterminate. The following two sufficient criteria can be useful

a) the popular Carleman’s criterion gives a sufficient condition that can be used to verify if the infinite sequence of moments determines uniquely a distribution. The criterion requires that:

∞∑   1
   --1-=  +∞
n=1μn2n
(10)

b) the Krein condition, compare Akhiezer [2], requires the distribution function to be absolutely continuous with density f(x) > 0, x > 0 (and f(x) = 0, x 0) and the moments of all order exist. If :

∫
  ∞  --ln-f (x)
  0   1 + x2 dx < ∞

then the distribution function is not determined by its moments.

c) In Shohat-Tamarkin, [29] it is possible to find a necessary and sufficient condition for the determinacy problem, but this condition seems too complicated to be very useful.

In next subsections, we review the different moment based approximations that have been proposed in the literature for solving the Asian option pricing problem.

5.1 Approximations based on the moments of A_T

On a practical side, the most used approximations in pricing Asian options are those based on the tentative of recovering the density of the average from its moments. In order to use this method, we need an expression for the moments of the arithmetic average. They can be found in Geman and Yor [17], page 359, and are given by the following expression:

                    ({  n          [(           )  ])}
μ  := E [A_ n] = n!-- ∑  d(v∕λ) exp    λ2j2-+ λjv  t
 n        T      λ2n( j=0 j           2            )
(11)

where:

 (β)      n ∏   [      2         2]-1
dj   =   2       (β + j) - (β + i)
           0i≤⁄=i≤jn
               (        )
  λ  =   σ;v =  r - σ2 ∕2 ∕σ
and in particular for n = 1, we have:
     e2(v+1)t --1
μ1 =  2(v + 1)

For the fact that the moment problem for the lognormal density is undeterminate, in the case of the arithmetic average one might expect that approximations based on a sequence of algebraic moments might either fail to converge, or converge to another distribution with the same moments. A difficult proof would no doubt be required and maybe a similar approximation could be numerically very unstable for more than a few moments. On the other hand, in the case of the arithmetic average we have a convolution of an infinite number of lognormal distributions and this convolution can change the distribution from a heavy tail density (the lognormal) to a light-tail density (the arithmetic average) and so there could a (little) possibility that the moments determine uniquely the distribution. But we need to remark that the effect of averaging can just reduce the value of the constant α, so that at finite values of x the tails seem decay fast, but for very large values of x the term xλ in the product αxλ becomes again important. Geman and Yor [17], page 359, have shown that the Carleman’s criterion is not satisfied in the case of the average of the geometric Brownian motion4 and so nothing can be said about the determinacy problem.

Although the determinacy problem in this case is unsolved, there are several approximations based on trying to recover the density of _
AT using the information content of its algebraic moments. For these approximations there exists no analytical bound for the error term as a function of the number of moments included in the approximation. Then for a fixed number of moments, the relative size of the error needs to be examined using numerical analysis and so nothing can be said about convergence.

5.1.1 The lognormal approximation (Levy)

To use a lognormal density as approximating density of the average is a very common approximation among the practitioners and it has been suggested in Turnbull and Wakeman [33] and in Levy [21].

If we assume that ln  _
AT has a normal distribution with mean m and variance v2, then, for any integer n :

    [ _ n]    nm+0.5n2v2
Elog AT   =  e

and where Elog stands for expected value assuming a (risk-neutral) lognormal distribution for the average. The idea is to choose m and v to fit the true mean and variance of  _
AT∕T. We obtain:

                 1
m   =   2logμ1 - --log μ2 - logT                     (12)
 2               2
v   =   log μ2 - 2 log μ1
and the price of the Asian call option is given by a modified Black-Scholes formula:
         m+v2∕2-rT          -rT
clog = S0e         N (d1) - e   KN  (d2)
(13)

where d1 =(ln S0∕K+m+v2)/v, d2=d1 - v.

5.1.2 The Edgeworth series approximation (Turnbull and Wakeman)

Using a two-parameters density we can capture the mean and the variance of the distribution, but we could have important differences between the fitted and the true skewness and kurtosis as we increase the volatility σ or the time to maturity of the option. For this reason Turnbull and Wakeman [33] and Ritchken and al. [27] propose to use a fourth-order Edgeworth series expansion. If flog(y;m, v2) is the lognormal density with parameters m and v2, then the Edgeworth approximation fedg(y;m,v2) to the true density is given by:

    (       )      (       )   ∑4 ki∂iflog (y;m, v2)
fedg  y;m, v2 = flog y;m, v2  +    ----------i----- + e(y)
                               i=1 i!      ∂y
(14)

where e(y) is the residual error term and ki is the difference between the i-th cumulant of the exact distribution and the approximate distribution:

ki = χi(f)- χi (l)

and:

χ1 (f)  =  μ1
             ( _     )2
χ2 (f)  =  E  AT  - μ1
             ( _     )3
χ3 (f)  =  E  AT  - μ1
             ( _     )4
χ4 (f)  =  E  AT  - μ1  - 3χ2 (f )

The first cumulant is the mean, the second the variance, the third a measure of skewness and the fourth a measure of kurtosis.

No analytical bound for the error term e(y) as a function of the number of terms included in the approximation exists. If all moments exist, it can be shown that e(y) goes uniformly to zero in y as the number of moments included in the approximation increases. For a fixed number of moments, the relative size of the error needs to be examined using numerical simulations. Moreover, if the error between the lognormal model and the true one is small, then the approximating technique on average adds noise, although not large in an absolute sense. For details compare Jarrow and Rudd [19].

Once the parameters m and v2 have been set according to (12) and then k1 = k2 = 0 in the expansion above, the option price formula is given by, compare Jarrow and Rudd [19]5 page 355:

                  ⌊                                                       ⌋
               S     k  ∂f   (y;m, v2)||         k   ∂2f  (y;m, v2)||
cedg = clog+e -rT--0⌈- -3 --log---------||       + --4 ---log--2------||       ⌉+e (K )
                T     6       ∂y      |y=T K∕S0  24       ∂y       |y=T K∕S0
(15)

where clog is given by (13) and where flog(y;m, v2) is the lognormal density with parameters m and v2.

In figure 1 we compare the transition densities of the underlying asset, of the arithmetic average obtained using Laplace inversion technique and the two approximations based on the lognormal density. Note that using the Edgeworth approximation, for high volatility levels ( √ --     )
 σ  T > 0.3the density becomes bimodal and besides it can assume negative values. The density of the underlying asset, of the arithmetic average (obtained with the Laplace inversion) and the lognormal and Edgeworth approximations. Note the negative values assumed by the Edgeworth approximation. Parameters: σ = 0.4, r = 0.09, t = 1.

5.1.3 The reciprocal gamma approximation (Milevsky and Posner)

Milevsky and Posner [25] have proved that the stationary density for the arithmetic average of a geometric Brownian motion is given by a reciprocal gamma density, i.e. the reciprocal of the average has a gamma density. The same result, using different techniques, have been previously obtained by Dufresne [11], De Schepper and al. [9], Yor [36], page 525.

In particular they have shown that _A = lim t→∞ _
Athas a density given by:

              β-αx -1-α
frg(x;α, β) = ----------exp(- 1∕(xβ))
                Γ (α)

where α = 1-2r/σ2and β = σ22. It is easy to verify that 1/_
A has a gamma density. The result is valid provided that r - σ22 < 0. Milevsky and Posner [25] state also that, according to the Kolmogorov-Smirnov criterion, this density is a better fit to the true density than the lognormal, where the true density has been estimated using Monte Carlo simulation.

The idea in Milevsky and Posner [25] is to approximate the density of the average using the stationary density given above, but where the parameters α and β are chosen in order to match the first two moments of the average. The moments of the reciprocal gamma distribution are:

    [(    )n]
Erg   -1--    = ------------1------------
      A ∞       β (α - 1)(α-  2)...(α-  n)

and are positive, provided that α > n. Fitting the first two moments to the corresponding moments of the arithmetic average S0  _
AT∕T, the parameters α and β are given by:

     S02μ2 - μ21      (S0 )2 μ2 - μ21
α =  T--μ----μ2 ,β =   T--  --μ-μ--
         2    1               2 1

and, assuming continuous averaging, the price crg of the Asian option with residual life equal to T is given by:

     1---e-rT    ( 1-        )    -rT    ( 1-     )
crg =    rT   S0G   K ;α - 1,β  - e   KG    K ;α,β
(16)

where G(x;α,β ) is the cumulative density function of the gamma density function and can be expressed in terms of the incomplete Gamma function, Press and al. [26], page 216. However, in the computation above it is preferable to use numerical integration.

5.2 Approximations based on the moments of  A1t.

Recently, Dufresne [12] has proved the following important and surprising result:

Theorem 1. (Dufresne) The distribution of 1∕At is determined by its moments.

Exploiting this result, it could be reasonable to recover the density of 1∕At and then the density of At. Indeed, for the computation of the price of the Asian option Dufresne uses Laguerre expansions of the price in terms of the moments of 1∕At. Unfortunately, the computation of the moments requires to solve numerically a differential-difference equation for the k-th moment of 1∕At and this fact may arise problems in the estimation of higher order moments. Indeed in order to obtain a good approximation, Dufresne has to use at least the first 11 moments for volatility levels around 0.3, and the first 5-6 moments for higher volatility levels. Besides this confirms the fact that a high volatility increases the accuracy of the approximation methods (Laplace transform, numerical solution of the PDE). The main advantage of the Dufresne approach is the convergence of the approximated price to the true one as we increase the number of considered moments. The main disadvantages are: a) no analytical bound on the error introduced by truncating the series after n terms is available, b) it is not guaranteed that the density assumes always positive values (compare figures in the Dufresne paper), c) the approximation works well only for σ√--
 T > 0.28.

5.3 Approximations based on the moments of lnA_T

In this section we show that:

a) the distribution of the logarithm of the average is determined by its moments;

b) we then show how to calculate the logarithmic moments;

c) we propose the Edgeworth expansion around a normal density for recovering the density using the information content of the logarithmic moments.

We have the following result:

Theorem 2. The distribution of the logarithm of the average is determined by its moments.

Proof: The proof is quite simple once we register the result in Dufresne [12] that E[ _ γ]
 AT exists finite for every real value of γ6 .This fact implies that the moment generating function of ln A_T:

  [    _]   ∫ +∞                  ∫ +∞                    [ _ γ]
E  eγlnAT  =      eγxfln _A (x )dx =      eγlnxf _A (x )dx = E  AT
             - ∞        T          0          T

exists finite for every real value of γ. Then the density function to be recovered has an exponentially decaying behavior7 and this is a sufficient condition for moment problem determinacy.

This result appears natural if we recall that the logarithm of a lognormal distribution has a normal distribution that is determined in a unique way from the sequence of its moments. The distribution of the average should be not very different from a lognormal distribution and so it is natural the result that the logarithm of the average has a distribution determined entirely from its moments.

The moments of ln _
AT. In order to calculate the moments of   _
ln AT we use the expression given8 in Geman and Yor [16] for the moments of  _
ATn, n real. If we define:

     ∫ h
Dh =     e2(Ws+vs)ds
      0

we have, by the scaling property of BM:

At =  4-D σ2t
      σ2   4

For every n 0 (n a real number), the Laplace transform of E   n
(D h) with respect to h, compare Geman and Yor[16], is given by:

         n     ∫+∞  -λh     n
Lh (1EΓ ((nD+1h)))Γ ((:=μ+v)0∕2+1e)Γ ((μE-v()D∕2-h)nd)h
=  λ--2nΓ ((μ-v)∕2)Γ (n+1+(μ+v)∕2)
(17)

so we have that:

                (  ( 4    )n)    ( 4 )n
Lh(E (Ant)) = Lh  E   -2Dh     =   --2   Lh(E (Dh )n)
                     σ            σ

If we differentiate m times the above expression with respect to n and then we set n = 0, we obtain:

             |
∂mLh-(E-(Ant))||              m
    ∂mn      |n=0 = Lh (E (ln  At))

i.e. the Laplace transform of the desired moments. Then the moments are obtained by inversion. The procedure is the following:

1) For every fixed λ, we compute the m-order derivative of (17) with respect to n.

2) We operate the numerical inversion choosing the values of T, σ and r and so we fix the corresponding value of h = σ2T∕4.

The differentiation can be done numerically applying the Cauchy-Goursat integral theorem9 , compare Churchill and Brown [5]. The inversion has been done using the Abate-Whitt algorithm with the same parameter settings as for the inversion of the Laplace transform of the Asian option price. Respect to the inversion of the Asian option price, we did not have problems for low volatility levels. Moreover, the computation of the Laplace transform in this case appears much less difficult than respect to the Geman-Yor formula.

5.3.1 The Edgeworth approximation (Fusai and Tagliani)

Once we have calculated the moments of ln _
AT we have to choose how many moments consider for recovering the density of ln  _
AT . In Fusai and Tagliani [14] are discussed three possible choices: a) the normal distribution, b) an Edgeworth expansion around a normal, c) a maximum entropy approximation. In this paper we present only the Edgeworth expansion around a normal density for its greater simplicity, however among the different moment based methods the maximum entropy method is the preferable.

The idea is to use as density for ln  _
AT an Edgeworth series approximation around a normal distribution and not, as usually it is done, to approximate the density of  _
AT using an Edgeworth series around a lognormal distribution.

Then the approximating density for ln _
AT is given by:

                                  4           (      )
f       (y)  =   f    (y;m, v2) + ∑  ki∂ifnorm-y;m,-v2- + e(y)      (18)
 edgnorm          norm             i=1 i!       ∂yi
                  (   _ )
        m   =   E  ln AT                                            (19)
         2        (  2 _ )    2
        v   =   E  ln  AT   - m                                     (20)
where fnorm(y) is the normal density with mean m and variance v2, e(y ) is the residual error term and ki is the difference between the i-th cumulant of the exact distribution and the approximate distribution:
ki = χi(f)- χi (l)

and:

                  (   _ )
χ1 (f ) =   μ1 = E  lnAT
             (   _      )2
χ2 (f ) =   E  lnAT  - μ1
             (   _      )3
χ3 (f ) =   E  lnAT  - μ1
             (   _      )4
χ4 (f ) =   E  lnAT  - μ1  - 3χ2 (f)
The approximating density for A_
  T is given by the following expression:
                (        2)
f _ (y) = fednorm--lny;m,-v--
 AT               y

Using this density we have computed the price of the Asian call option exploiting the fact that:

∫ +∞ ( S     )+               S∫ +∞  (      T)+
 0    yT-- K    fA_T (y)dy  = T- KT ∕S y - K S-  fA_T (y)dy
                           = -S∫ 1(KS-- K  T)+ f _ (KT-) KS-dw
                             T  0  Tw      S    AT  Sw   Tw2

so that the infinite integration range has been reduced to a finite one.

Obviously this approximation has the limit that there is no expression for bounding the error and moreover it is not guaranteed the positivity of the density, but respect to the Edgeworth expansion around a lognormal distribution the current approximation appears much more robust.

6 A comparison between the different approximations

The scheme below compares the different methods that exploit the information content of a sequence of moments and the properties of the approximation used for recovering the density of the average and the option price. Regard to the Edgeworth approximation using the algebraic moments, the convergence to the true density is not known, because the determinacy problem is unsolved. Instead using the Edgeworth approximation with logarithm moments the determinacy problem is positively solved, so that it makes sense to investigate the convergence of the approximation that is uniform as the number of moments considered increase.











MomentsDeterminacyApproximantConvergence Notes





E0(Ant) Not known Edgeworth Not known
very simple;
negative densities;
no error bound.





E0((   )n)
   1At Yes Laguerre Yes
not very fast;
numerical problems;
no error bound.





E0(lnAn )
     t Yes Edgeworth Yes
simple;
good accuracy;
no error bound.





E0     n
(lnA t) Yes ME Yes
not very fast,
very accurate;
error bound.










The different moments based approximations can be compared in several ways:

  1. evaluating with the Kolmogorov-Smirnov criterion the distance between the approximating distribution and the distribution obtained by the Laplace inversion; unfortunately, as we will this measure does not seem to be useful once we consider the pricing problem.
  2. compare the different approximations with the lower and the upper bound to the option price as given in Rogers and Shi [28] and in Thompson [31]. The lower bound appears to be very strict, but from this bound we cannot recover easily information regard to the density function of the average or the Greeks of the contract;
  3. evaluate the difference (mean, absolute, percentage) between the option price coming from the different approximations and the one obtained by the Laplace transform inversion.

As we can see from table (4) where we report the Kolmogorov Smirnov measure between the density obtained with the numerical inversion and the different approximations, we can draw the following conclusions:

a) for low volatility levels (σ√T-- < 0.2), the different approximations perform equally well and the first two moments seems to be sufficient to recover the density of the average. The lognormal approximation based on the logarithmic moments does slightly better than the lognormal approximation based on the power moments. The reciprocal gamma behaves better than the lognormal. However as we increase the volatility level, it is important to include higher order moments.

b) if we restrict our attention to approximations based on the moments of the average the GB2 approximation dominates the Edgeworth approximation. Moreover using the Edgeworth approximation, for high volatility levels we obtain a bimodal density. Besides the density can assume negative values. This fact should discourage from using an Edgeworth approximation when the distribution is not uniquely determined by its moments.

c) In every case the approximations based on the moments of the logarithm of the average perform better and in this sense the Edgeworth approximation around a normal and the ME with four moments are very good.

d) although not implemented in Premia, with regard to the ME approximation to use four moments instead of six seems very reliable.



















Table 4: Kolmogorov-Smirnov distance between
the true distribution and the different approximations









Volatility Fitlog EdgeLogRGamma GB2 FitNormEdgeNorm ME4 ME6








0.1 0.03055 0.03301 0.03441 0.03248 0.03053 0.03282 0.03199 0.03168 0.2 0.01327 0.01821 0.01961 0.01569 0.01246 0.01661 0.01656 0.01647 0.3 0.01071 0.01725 0.01549 0.00990 0.00833 0.01123 0.01142 0.01132 0.4 0.01200 0.02530 0.01397 0.00673 0.00801 0.00851 0.00889 0.00886 0.5 0.01470 0.06956 0.01344 0.00617 0.00859 0.00687 0.00746 0.00740








Maximum0.03055 0.06956 0.03441 0.03248 0.03053 0.03282 0.03199 0.03168








Average 0.01624 0.03266  0.01938 0.01420 0.01358 0.01521  0.01527 0.01514

















If we compare the different approximations when we have to price Asian options, we can see a synthesis of the results in table 5 and in more detail in tables 6 and 6a-6f. In table 5, using as benchmark the result of the Laplace inversion, we report the average of the mean squared error (MSE), the average mean absolute error (MAE) and the absolute maximum error (MAXA), where the average is taken considering the above errors in tables 6a-6f, when we vary the volatility level (σ  --
√ T = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5). In the case σ√T-- = 0.05 the column with the Laplace inversion in table 6 has been substituted by the average of the low and the upper bound, because in this case we had numerical problems in the computation of the Laplace transform. In table 5 we report also the percentage of times the price of a given approximation stays inside the lower and upper bound given in Rogers and Shi [28] and in Thompson [31].

Table 5: Evaluation of different approximations















MSE MAE MAXA %Inside





Abate-Whitt
-
- - 96.67%
ME4
0.00028
0.000910.00380 64.44%
Lower
0.00049
0.001830.00647 -
EdgeNormale
0.00083
0.003040.01279 31.11%
GB2
0.00203
0.006820.03363 31.11%
Upper
0.00235
0.008880.03631 -
RecGamma
0.01362
0.050090.19608 2.22%
Lognormal
0.01661
0.059330.24819 8.89%
Edgeworth
0.01928
0.060850.56647 31.11%










We can draw the following conclusions:

a) Comparing tables 6a-6f, we can see the very good results that we can obtain with the numerical Laplace inversion. Indeed also with a volatility level at 10% when the bounds to the option price are very tight, the price obtained stays inside this small range: e.g. when σ√ -
  t = 0.1, r = 0.05, S0 = 100 and K = 100 the lower bound is 3.64134 and the upper bound is 3.64157 and the Laplace inversion10 gives 3.64139.

b) A second aspect are the accurate results obtained with the maximum entropy estimate based on four moments. The fact that the percentage of times this estimate stays inside the bound is only 64.4% is due to the difficulties in estimating the ME distribution for low volatility levels (σ√-
 t = 5%). However, as we increase the volatility, also the ME estimate stays inside the bounds and the improvement respect to the lognormal approximation is around 65 times. The maximum error committed with a ME approach with four moments is equal at 0.0032 when σ=10% and reduces to 0.0015 as σ increases at 50%.

c) The results of the Edgeworth expansion around the normal distribution compared to the Edgeworth expansion around a lognormal density, confirms the fact that using the (first four) logarithmic moments is a good choice.

d) Among the approximations that use the power moments, the Edgeworth expansion around a lognormal distribution performs very badly, mainly as we increase the volatility level. The use of the Reciprocal Gamma approximation does not seem to improve significantly the results respect to the lognormal approximation. Another interesting point is that also with very low volatility levels the prices obtained with the approximations based on the power moments stay very often outside the bounds. However the maximum absolute difference between the lognormal approximation and the Laplace answer is 0.008 for volatility at the 10% level and increases at 0.248 for a volatility at 50%. This non-enormous absolute error can explain the success of the lognormal approximation among practitioners.

7 Conclusion

In the present paper we have given new insights regard the problem of pricing Asian options using moments. We have showed, both theoretically both empirically, how approximations based on the moments of the logarithm of the arithmetic average are more reliable than approximations based on the algebraic moments of the arithmetic average. An important extension left for future work is to examine if approximations based on the logarithmic moments are possible also in the case of discrete monitoring of the average and in the case of floating strike Asian options. Indeed typically traded options are written on baskets of dividend-paying securities: even the evaluation of moments of the average for which there exists explicit formulae, becomes difficult beyond order three.

References

[1]   Abate, J. and Whitt, W. (1992), The Fourier-series method for inverting transforms of probability distributions, Queueing Systems Theory Appl. 10 5-88.

[2]   Akhiezer N. I. (1965), The Classical Moment Problem and some Related Questions in Analysis, Oliver & Boyd, London.

[3]   Alziary B., J. Decamps and P. Koehl (1997), A P.D.E. Approach to Asian options: Analytical and Numerical Evidence, Journal of Banking and Finance, 21, 613-640.

[4]   Barraquand, J. and Pudet T. (1996), Pricing of American path-dependent contingent claims, Mathematical Finance, Vol. 6, 17-51.

[5]   Churchill R. V. and J. W. Brown (1990), Complex Variables and Applications, Fifth edition, Mc-Graw-Hill International Editions.

[6]   Clelow L. J. and A. P. Carverhill (1994), Quicker on the Curves, Risk, Vol. 7 No. 5, 63-64.

[7]   Curran, M. (1992), Beyond average intelligence, Risk, Vol. 5, No. 10, pag. 60.

[8]   Curran, M. (1992), Valuing Asian and Portfolio Options by Conditioning on the Geometric Mean Price, Management Science, Vol. 40, No. 12, pagg.1705-1711.

[9]   De Schepper A., M. Teunen and M. Goovaerts (1994), An Analytical Inversion of a Laplace transform related to annuities certain, Insurance: Math. Econonom., 14, 33-37.

[10]   Dubner H. and Abate J. (1968), Numerical Inversion of Laplace transforms by relating them to the finite Fourier cosine transform, J. ACM 15, 1 January, 115-123.

[11]   Dufresne D. (1989), Weak Convergence of Random Growth Processes with Applications to Insurance, Insurance: Math. Econonom., 8, 187-201.

[12]   Dufresne D., Laguerre Series for Asian and Other Options, Mathematical Finance v. 10, n. 4, pagg. 407-428., October, 2000.

[13]   Fu M. C., D. Madan and T. Wang (1998), Pricing Continuous Asian Options: a comparison of MonteCarlo and Laplace transform inversion methods, The Journal of Computational Finance, Vol. 2, No. 1, 49-74.

[14]   Fusai G. and A. Tagliani (2000), Accurate valuation of Asian Options using moments, w.p. U. of Novara, Italy, submitted.

[15]   Geman H. and A. Eydeland (1995), Domino Effect, Risk, Vol. 8, n. 4, april, 65-67.

[16]   Geman H. and M. Yor (1992), Quelques relations entre processus de Bessel, options asiatique et fonctions confluentes hypergéométriques, C.R. Acad. Sci. Paris, t. 314, Série I, p. 471-474.

[17]   Geman H. and M. Yor (1993), Bessel Processes, Asian Options and Perpetuities, Mathematical Finance, vol. 3, no. 4, pagg. 349-375.

[18]   Hull, J. and A. White (1993), Efficient Procedures for Valuing European and American Path Dependent Options, Journal of Derivatives, Vol. 1, 21-31.

[19]   Jarrow R. and A. Rudd (1982), Approximate Option Valuation for Arbitrary Stochastic Processes, Journal of Financial Economics, vol. 10, pp. 347-369.

[20]   Kemna A. and A. Vorst (1990), A Pricing method for Options based on Average Values, Journal of Banking and Finance, vol. 14, 113-129.

[21]   Levy E. (1992), Pricing European Average Rate Currency Options, Journal of International Money and Finance, vol. 11, 474-491.

[22]   Lindsay, B. G. and P. Basak (1995), Moments determine the tail of a distribution (but not much else), Technical report 95-7, Center for likelihood Studies, Dept. of Statistics, The Pennsylvania State University, University Park, PA.

[23]   Lindsay, B. G. and K. Roeder (1997), Moment-based oscillation properties of mixture models, Ann. Statist., 25, pagg. 378-386.

[24]   Lipton A. (1999), Similarities via self-similarities, Risk, vol. 12, no. 9, 101-105, September.

[25]   Milevsky M. A. and S. E. Posner (1998), Asian Options, the sum of lognormals and the reciprocal gamma distribution, Journal of Financial and Quantitative Analysis, vol. 33, no. 3, pagg. 409-422, Sept..

[26]   Press, W. H., Teukolsky S. A., Vetterling W. T. and Flannery B. P. (1997). Numerical Recipes in C, Cambridge University Press, version 2.08

[27]   Ritchken P., L. Sankarasubramanian and A. M. Vijh (1990), The Valuation of Path-Dependent contracts on the average, Management Science, vol. 39, no. 10, 1202-1213.

[28]   Rogers L. C. G. and Shi Z. (1992), The Value of an Asian option, J. of Applied Probability, Vol. 32, pagg. 1077-1088.

[29]   Shohat J. A. and J.D. Tamarkin (1943), The problem of moments, AMS Mathematical Survey 1, Providence RI, (1943).627-1628.

[30]   Stoyanov J. (1997), Counterexamples in Probability, 2nd edition, Wiley Series in Probability and Statistics.

[31]   Thompson G. W. P. (1999), Topics in Mathematical Finance, PH.D. Thesis, Queens’ College, January.

[32]   Thompson G. W. P. (1999), Fast narrow bounds on the Value of Asian options, w.p. Judge Institute of Management Studies, U. of Cambridge.

[33]   Turnbull S. and L. Wakeman (1991), A Quick Algorithm for Pricing European Average Options, Journal of Financial and Quantitative Analysis, vol. 26, pagg. 377-389.

[34]   Vorst, T. (1996), Averaging Options, in I. Nelken, The Handbook of Exotic Options, Chpt. 6, pp. 175-199, IRWIN.

[35]   Wilmott, P. Dewynne, J.N. and Howison, S. (1993), Option Pricing: Mathematical Models and Computation, Oxford Financial Press

[36]   Yor M. (1992), On some exponential functionals of Brownian Motion, Adv. Appl. Prob. 24, pagg. 509-531.

[37]   Yor M. (1992), Sur les lois des fonctionnelles exponentielles du mouvement brownien, considérées en certains instants aléatoires, C. R. Acad. Sci Paris Série I, 314, 951-956.

[38]   Yor M. (1992), Some Aspects of Brownian Motion. Part I: Some Special Functionals. Birkhauser Verlag. Basel - Boston - Berlin.

[39]   Zvan R., P. A. Forsyth and K. R. Vetzal (1997), Robust Numerical methods for PDE models of Asian options, The Journal of Computational Finance, Vol. 1, N. 2, 39-77.