Pricing Asian options using mpmath with automatic precision control (Python recipe) by Dieter Kadelka

ActiveState Code (http://code.activestate.com/recipes/576964/)

# -*- coding: iso-8859-1 -*-

# mp_laplace.py
# laplace.py with mpmath
#   appropriate for high precision

# Talbot suggested that the Bromwich line be deformed into a contour that begins
# and ends in the left half plane, i.e., z \to \infty at both ends.
# Due to the exponential factor the integrand decays rapidly
# on such a contour. In such situations the trapezoidal rule converge
# extraordinarily rapidly.
# For example here we compute the inverse transform of F(s) = 1/(s+1) at t = 1
#
# >>> error = Talbot(1,24)-exp(-1)
# >>> error
#   (3.3306690738754696e-015+0j)
#
# Talbot method is very powerful here we see an error of 3.3e-015
# with only 24 function evaluations
#
# Created by Fernando Damian Nieuwveldt      
# email:fdnieuwveldt@gmail.com
# Date : 25 October 2009
#
# Adapted to mpmath and classes by Dieter Kadelka
# email: Dieter.Kadelka@kit.edu
# Date : 27 October 2009
# Automatic precision control by D. Kadelka 2009-11-26
#
# Reference
# L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.

try:
  import psyco
  psyco.full()
except ImportError:
  print 'Psyco not installed, the program will just run slower'

from mpmath import mp,mpf,mpc,pi,sin,tan,exp,floor,log10

# testfunction: Laplace-transform of exp(-t)
def F(s):
  return 1.0/(s+1.0)

class Talbot(object):

  # parameters from
  # T. Schmelzer, L.N. Trefethen, SIAM J. Numer. Anal. 45 (2007) 558-571
  c1 = mpf('0.5017')
  c2 = mpf('0.6407')
  c3 = mpf('0.6122')
  c4 = mpc('0','0.2645')
  # High precision of these parameters not needed
      
  def __init__(self,F=F,shift=0.0,prec=50):
    self.F = F
    # test = Talbot() or test = Talbot(F) initializes with testfunction F
    # Assumption: F realvalued and analytic

    self.shift = mpf(shift)
    # Shift contour to the right in case there is a pole on the 
    #   positive real axis :
    # Note the contour will not be optimal since it was originally devoloped 
    #   for function with singularities on the negative real axis For example
    #   take F(s) = 1/(s-1), it has a pole at s = 1, the contour needs to be 
    #   shifted with one unit, i.e shift  = 1. 
    # But in the test example no shifting is necessary
 
    self.N = 12
    # with double precision this constant N seems to best for the testfunction
    #   given. For N = 11 or N = 13 the error is larger (for this special
    #   testfunction).

    self.prec = prec
    # calculations with prec more digits

  def __call__(self,t):
  
    with mp.extradps(self.prec):
      t = mpf(t)
      if t == 0:
        print "ERROR:   Inverse transform can not be calculated for t=0"
        return ("Error");
            
      N = 2*self.N
      # Initiate the stepsize (mit aktueller Präsision)
      h = 2*pi/N
   
    # The for loop is evaluating the Laplace inversion at each point theta i
    #   which is based on the trapezoidal rule
      ans =  0.0
      for k in range(self.N):
        theta = -pi + (k+0.5)*h
        z = self.shift + N/t*(Talbot.c1*theta/tan(Talbot.c2*theta) - Talbot.c3 + Talbot.c4*theta)
        dz = N/t * (-Talbot.c1*Talbot.c2*theta/sin(Talbot.c2*theta)**2 + Talbot.c1/tan(Talbot.c2*theta)+Talbot.c4)
        v1 = exp(z*t)*dz
        prec = floor(max(log10(abs(v1)),0))
        with mp.extradps(prec):
          value = self.F(z)
        ans += v1*value
            
      return ((h/pi)*ans).imag

*********************************************************************************
# -*- coding: iso-8859-1 -*-

# asian.py

# Title : Numerical inversion of the Laplace transform for pricing Asian options
#         The Geman and Yor model
# 
# Numerical inversion is done by Asian's method.
#
################################################################################
## Created by Fernando Damian Nieuwveldt                                        
## Date : 26 October 2009                                                      
## email : fdnieuwveldt@gmail.com
## This was part work of my masters thesis (The Asian 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. Asian quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.

# adapted to mp_laplace by D. Kadelka 2009-11-17
# Automatic precision control by D. Kadelka 2009-11-26
# email: Dieter.Kadelka@stoch.uni-karlsruhe.de

# Example:
# from asian import Asian
# f = Asian()
# print f
# Pricing Asian options: The Geman and Yor model with
#   S = 2, K = 2, T = 1, t = 0, sig = 0.1, r = 0.02
# print f()
# 0.0559860415440029
# f.ch_sig('0.05')
# print f
# Pricing Asian options: The Geman and Yor model with
#   S = 2, K = 2, T = 1, t = 0, sig = 0.05, r = 0.02
# print f()
# 0.0345709175410301
# f.N = 100
# print f()
# 0.0339410537085201
# from mpmath import mp
# mp.dps = 50
# f.update()
# 0.033941053708520319031364170122438704213486236188948

try:
  import psyco
  psyco.full()
except ImportError:
  print 'Psyco not installed, the program will just run slower'

from mpmath import mp,mpf,mpc,pi,sin,tan,exp,gamma,hyp1f1,sqrt,log10,floor
from mp_laplace import Talbot

class Asian(object):

  def G(self,s): # Laplace-Transform
    zz = 2*self.v + 2 + s
    mu = sqrt(self.v**2+2*zz)
    a  = mu/2 - self.v/2 - 1
    b  = mu/2 + self.v/2 + 2
    v1 = (2*self.alp)**(-a)*gamma(b)/gamma(mu+1)/(zz*(zz - 2*(1 + self.v)))
    prec = floor(max(log10(abs(v1)),mp.dps))+self.prec
    # additional precision needed for computation of hyp1f1
    with mp.extradps(prec):
      value = hyp1f1(a,mu + 1,self.beta)*v1
    return value
 
  def update(self):
  # Geman and Yor's variable
  # possibly with infinite precision (strings)
    self.S = mpf(self.parameter['S'])
    self.K = mpf(self.parameter['K'])
    self.T = mpf(self.parameter['T'])
    self.t = mpf(self.parameter['t'])
    self.sig = mpf(self.parameter['sig'])
    self.r = mpf(self.parameter['r'])

    self.v = 2*self.r/(self.sig**2) - 1
    self.alp  = self.sig**2/(4*self.S)*self.K*self.T
    self.beta = -1/(2*self.alp)
    
    self.f.shift = self.shift
 
  def __init__(self,S=2,K=2,T=1,t=0,sig='0.1',r='0.02',N=50,shift=0.0,prec=0):
    # Strings allowed for infinite precision
    # prec compensates rounding errors not catched with automatic precision control
    # parameters may be changed later
    # after changing mp.dps or any of these parameters (except prec, N and t),
    #   update (v,alp,beta depend on these parameters)

    self.N = N  
    self.shift = shift
    self.prec = max(prec,0)
    self.parameter = {'S':S,'K':K,'T':T,'t':t,'sig':sig,'r':r}
  # input: possibly strings with infinite precision
    self.f = Talbot(self.G,shift=self.shift,prec=0)
    self.update()

  def __call__(self):
    # Initialize the stepsize (with actual precision)
    self.f.N = self.N
    tau  = ((self.sig**2)/4)*(self.T - self.t)
    # Evaluation of the integral at tau
    return 4*exp(tau*(2*self.v+2))*exp(-self.r*(self.T - self.t))*self.S/(self.T*self.sig**2)*self.f(tau)

  # Update Parameters
  def ch_S(self,S):
    self.parameter['S'] = S
    self.update()

  def ch_K(self,K):
    self.parameter['K'] = K
    self.update()

  def ch_T(self,T):
    self.parameter['T'] = T
    self.update()

  def ch_t(self,t):
    self.parameter['t'] = t
    self.update()

  def ch_r(self,r):
    self.parameter['r'] = r
    self.update()

  def ch_sig(self,sig):
    self.parameter['sig'] = sig
    self.update()
 
  # Actual Parametes
  def __str__(self):
    s = 'Pricing Asian options: The Geman and Yor model with\n'
    s += "  S = %(S)s, K = %(K)s, T = %(T)s, t = %(t)s, sig = %(sig)s, r = %(r)s" % self.parameter
    return s