I present a numerical method for pricing Asian options. The method is based on the numerical inversion of the Laplace transform. The inversion method that is used is based on Talbot contours. It is known that Geman and Yor's formula is computational expensive for low volatility cases. By using Talbots method we can reduce the timing for the low volatility cases, at least to \sigma ~ 0.05. Afterwards the method start to converge slowly. In the literature for \sigma = 0.1 the Geman and Yor formula converges slowly.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | ```
# Title : Numerical inversion of the Laplace transform for pricing Asian options
# The Geman and Yor model
#
# Numerical inversion is done by Talbot's method.
#
################################################################################
## Created by Fernando Damian Nieuwveldt
## Date : 26 October 2009
## email : fdnieuwveldt@gmail.com
## This was part work of my masters thesis (The Talbot method not mpmath part)
## in Applied Mathematics at the University of Stellenbosch, South Africa
## Thesis title : A Survey of Computational Methods for Pricing Asian Options
## For reference details contact me via email.
################################################################################
# Example :
# Asian(2,2,1,0,0.1,0.02,100)
# 0.0559860415440030213974642963090994900722---mp.dps = 100
# Asian(2,2,1,0,0.05,0.02,250)
# 0.03394203103227322980773---mp.dps = 150
#
# NB : Computational time increases as the volatility becomes small, because of
# the argument for the hypergeometric function becomes large
#
# H. Geman and M. Yor. Bessel processes, Asian options and perpetuities.
# Mathematical Finance, 3:349 375, 1993.
# L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.
from mpmath import mp,mpf,mpc,pi,sin,tan,exp,gamma,hyp1f1,sqrt
def Asian(S,K,T,t,sig,r,N):
# Assigning multi precision
S = mpf(S); K = mpf(K)
sig = mpf(sig)
T = mpf(T); t = mpf(t)
r = mpf(r)
# Geman and Yor's variable
tau = mpf(((sig**2)/4)*(T - t))
v = mpf(2*r/(sig**2) - 1)
alp = mpf(sig**2/(4*S)*K*T)
beta = mpf(-1/(2*alp))
tau = mpf(tau)
N = mpf(N)
# Initiate the stepsize
h = 2*pi/N;
mp.dps = 100
c1 = mpf('0.5017')
c2 = mpf('0.6407')
c3 = mpf('0.6122')
c4 = mpc('0','0.2645')
# The for loop is evaluating the Laplace inversion at each point theta which is based on the trapezoidal
# rule
ans = 0.0
for k in range(N/2): # N/2 : symmetry
theta = -pi + (k+0.5)*h
z = 2*v+2 + N/tau*(c1*theta/tan(c2*theta) - c3 + c4*theta)
dz = N/tau*(-c1*c2*theta/sin(c2*theta)**2 + c1/tan(c2*theta)+c4)
zz = N/tau*(c1*theta/tan(c2*theta) - c3 + c4*theta)
mu = sqrt(v**2+2*z)
a = mu/2 - v/2 - 1
b = mu/2 + v/2 + 2
G = (2*alp)**(-a)*gamma(b)/gamma(mu+1)*hyp1f1(a,mu + 1,beta)/(z*(z - 2*(1 + v)))
ans += exp(zz*tau)*G*dz
return 2*exp(tau*(2*v+2))*exp(-r*(T - t))*4*S/(T*sig**2)*h/(2j*pi)*ans
``` |

It is known that Geman and Yor's formula is computational expensive for low volatility cases. By using Talbots method we can reduce the timing for the low volatility cases, at least to \sigma ~ 0.05. Afterwards the method start to converge slowly. In the literature for \sigma = 0.1 the Geman and Yor formula converges slowly.

Also see Recipe 576964: Pricing Asian options using mpmath with automatic precision control, which based on this algorithm but with automatic precision.