################################################################################
## Created by Fernando Damian Nieuwveldt
## Date : 26 October 2009
## email : fdnieuwveldt@gmail.com
## This was part work of my masters thesis 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 :
# >>>Hypergeometric1F1(5,2,100-1000j,50)
# (7.0028644420385242e+050+8.9737757674587575e+050j)
#
# using mpmath
# Hypergeometric1F1(-1000,1,1000,2000)
# mpc(real='-2.5938207833620713892854e+215', imag='-1.1411010130233474563877e+202')
#
# Other references :
# A. Erdélyi, W. Magnus, F. Oberhettinger, and F. Tricomi. Higher
# transcendental functions. Vol. I. Robert E. Krieger Publishing Co.
# Inc., Melbourne, Fla., 1981, pg 273
#
# 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
mp.dps = 100
def F(a,b,t,s):
return s**(-b)*(1+1.0/s)**(a-b)
def Hypergeometric1F1(a,b,t,N):
# Initiate the stepsize
h = 2*pi/N;
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):
theta = -pi + (k+0.5)*h
z = N/t*(c1*theta/tan(c2*theta) - c3 + c4*theta)
dz = N/t*(-c1*c2*theta/sin(c2*theta)**2 + c1/tan(c2*theta)+c4)
ans += exp(z*t)*F(a,b,t,z)*dz
return gamma(b)*t**(1-b)*exp(t)*((h/(2j*pi))*ans)
Diff to Previous Revision
--- revision 2 2010-03-22 04:48:19
+++ revision 3 2010-03-22 05:00:20
@@ -35,10 +35,8 @@
def Hypergeometric1F1(a,b,t,N):
# Initiate the stepsize
- h = 2*pi/N;
+ h = 2*pi/N;
-
-
c1 = mpf('0.5017')
c2 = mpf('0.6407')
c3 = mpf('0.6122')
