Quantopian's community platform is shutting down. Please read this post for more information and download your code.
Back to Community
[Python] Black&Scholes PDE finite difference method

Hi,
I am trying to make again my scholar projet.
One of them was to solve the Black and Scholes PDE with finite different methods.

I tried some codes but didnt get a right result.
See below my last try :

import numpy as np

_vol = 0.40  
_r = 0.10  
_K = 50  
T = 0.5

m = 100 # time  
n = 200 # space

dt = T / m # time step  
dx = 2 * _K / (n+1) # space step

print("dt = ", dt)  
print("dx = ", dx)

l = np.zeros((m, n))


def payoff(S):  
    # Call Euro  
    # strike = 100  
    return max(S - _K, 0)  
# Maturity condition :  
for j in range(n):  
    l[0][j] = payoff(j * dx)

for i in range(m-1):  
    for j in range(n):  
        if (j == 0) :  
            l[i+1][j] = 0  
        elif (j == n-1):  
            l[i+1][j] = j * dx - _K * np.exp(-_r * (T - i * dt))  
        else:  
            # explicit method :  
            l[i+1][j] = (1/2 * (_vol * j)**2 * dt - 1/2 * _r * j * dt) * l[i][j-1] + (1 - (_vol * j)**2 * dt - _r * dt) * l[i][j] + (1/2 * (_vol * j)**2 * dt + 1/2 * _r * j * dt) * l[i][j+1]  

You can see this link : http://planetmath.org/solvingtheblackscholespdebyfinitedifferences

I tried the explicit & implicit and never get the correct result..
For example, the price of a call, t=0, K=50, S=50 gives a, explosive negative result : -5.9.e+117 --"
I waited a price close to zero..
Please someone has already get it with a easy to reade code ?
Do not hesitate if you find a mistake in this code.

Thanks & regards,

5 responses

Cory's python notebooks on github

I have some notebooks on github, link above, that will help. There is an explicit method for european vanilla and barrier options.
Hope they help.

Cory

Hi Cory,

Thank you for you answer.

I looked your code. Even if the result is ok, i didnt understand why the delta and gamma are in a negative sens ?
You wrote a formula like : delta = (f(x-h) - f(x+h) ) / 2*h

When I rewrite your formula for a function u(x,t) the price of an option, I have :
[u(x,t)-u(x,t-1)]/dt + r.x.[u(x-1,t), u(x+1)]/2dx + 0.5.vol².x².[u(x-1,t)-2.u(x,t)+u(x+1,t)]/dx² -r.u(x,t) = 0

This formula seems not match with Black and Scholes PDE ? In particular, the delta and gamma are not in the same sign ?

Thanks a lot.

Srn

I think it just looks confusing the way the y-1 in the code gives you the upper value and y+1 gives you the lower value on the y axis in the grid.
I put in some code to show more explicitly what values it is putting in the calculation and what Delta it is calculating.

Let me know if that clears it up.

Best

Ok thanks a lot. I didnt pay attention to the asset_price decreases.

A last (maybe) question, what is the boundary rule for maximum asset price :

    #Set Ceiling Value  
    value_matrix[0,-x-1] = 2 * value_matrix[1,-x-1] - value_matrix[2,-x-1]  

What is the idea of this formula ? Why not apply directly the discounted rule of payoff ? Like the lower boundary rule ?
Thank you for your spent time.

ps: sorry for my english

For derivatives on stocks, the lower bound is known and trivial as the stock price cannot go below zero.
The upper bound is technically infinity, for practical purposes you should use 2 to 3 times the strike.
Applying the floor rule for the upper bound would only work for simple European calls and puts. If you used finite differences for american options, barrier, binary, etc..., the discounting rule does not work.
The ceiling value is a linear extrapolation based on the two lower vertical points in the grid. In the code: value_matrix[*0*,-x-1] = 2 * value_matrix[*1*,-x-1] - value_matrix[*2*,-x-1] , the 0 is the ceiling value, the 1 is one step down on the y axis and the 2 is two steps down on the y axis.

Hope that helps.
Best.