# 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
