Welcome, guest | Sign In | My Account | Store | Cart
```# 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
```