Welcome, guest | Sign In | My Account | Store | Cart

Solving the Black-Scholes PDE with laplace inversion:Revised (Python recipe) by Fernando Nieuwveldt
ActiveState Code (http://code.activestate.com/recipes/577142/)

I originally posted this code in Recipe 577132 and this is a repost of that recipe with corrections since there was an error in the original recipe. Added here is an error analysis to show the effectiveness of the Laplace inversion method for pricing European options. One can test the accuracy of this method to the finite difference schemes. The laplace transform of Black-Scholes PDE was taken and the result was inverted using the Talbot method for numerical inversion. For a derivation of the laplace transform of the Black-Scholes PDE, see for instance www.wilmott.com/pdfs/020310_skachkov.pdf.

Python, 49 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
## Laplace inversion of a European option
## Fernando Nieuwveldt
## email: fdnieuwveldt@gmail.com
## Example :
## >>> BSLaplace(8.0,10.0,1,0,0.03,0.3,50,1)
## mpf('0.4123752414202130413263844663926875')
## Referenece:
## A Survey of computational methods for pricing Asian options,masters thesis
## Fernando Nieuwveldt, University of Stellenbosch, South Africa

from mpmath import mp,mpf,mpc,pi,sin,tan,sqrt,exp,ln
mp.dps = 32
  
def bs(x,b1,b2,eps1,eps2,s,r,phi):
    """phi = 1 for a call -1 for a put"""    
    if x <= 0:
       return b1*exp(eps1*x) + (phi-1)/2*(exp(x)/s - 1/(s + r))
    else:
       return b2*exp(eps2*x) + (phi+1)/2*(exp(x)/s - 1/(s + r))

def BSLaplace(S,K,T,t,r,sig,N,phi): 
        """Solving the Black Scholes PDE in the Laplace domain"""
        x   = ln(S/K)     
        r = mpf(r);sig = mpf(sig);T = mpf(T);t=mpf(t)
        S = mpf(S);K = mpf(K);x=mpf(x)
        mu  = r - 0.5*(sig**2)
       
        tau = T - t   
        c1 = mpf('0.5017')
        c2 = mpf('0.6407')
        c3 = mpf('0.6122')
        c4 = mpc('0','0.2645')        
        
        ans = 0.0
        h = 2*pi/N
        h = mpf(h)
        for k in range(N/2): # Use symmetry
            theta = -pi + (k+0.5)*h
            z     =  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)
            eps1  =  (-mu + sqrt(mu**2 + 2*(sig**2)*(z+r)))/(sig**2)
            eps2  =  (-mu - sqrt(mu**2 + 2*(sig**2)*(z+r)))/(sig**2)
            b1    =  1/(eps1-eps2)*(eps2/(z+r) + (1 - eps2)/z)
            b2    =  1/(eps1-eps2)*(eps1/(z+r) + (1 - eps1)/z)
            ans  +=  exp(z*tau)*bs(x,b1,b2,eps1,eps2,z,r,phi)*dz
            val = (K*(h/(2j*pi)*ans)).real
           
            
        return 2*val

For the error analysis I used the call option for numerical purposes. We have two cases one for S > K or x > 0, and S <= K or x<=0. For the x<=0, the relative error is 

>>> (BSLaplace(10.0,10.0,1,0,0.03,0.3,50,1)-call(10.0,10.0,0.03,0.3,1.0))/call(10.0,10.0,0.03,0.3,1.0)
mpf('-9.2423932491484155320148874543823695e-29')

For the case x > 0, relative error obtained is 

>>> (BSLaplace(12.0,10.0,1,0,0.03,0.3,50,1)-call(12.0,10.0,0.03,0.3,1.0))/call(12.0,10.0,0.03,0.3,1.0)
mpf('-1.2995251591940985425649439119791792e-16')

For the exact solution I used the following 

>>> def call(S,K,r,sig,tau):
         d1 = (ln(S/K) + (r + 0.5*sig**2)*(tau))/(sig*sqrt(tau))
         d2 = d1 - sig*sqrt(tau)
         N1 = 0.5*(1+erf(d1/sqrt(2)))
         N2 = 0.5*(1+erf(d2/sqrt(2)))
         C = S*N1-K*exp(-r*(tau))*N2
         return C

Since an exact solution exists, this method should be used as a comparison to existing methods such as finite difference methods for instance.

Tags: black, laplace, scholes, talbot