# 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