# Black-Scholes explicit Euler implementation python

I've written some code for the explicit finite difference method to solve the BS equation.

For certain sets of parameters (time-steps and asset-steps) I get a stable but wrong solution. For others, I get everything diverging to nonsense. I can't spot any errors in my implementation. Can you?

import numpy as np
import matplotlib.pyplot as plt

def A(vol, asset_step, risk_free_rate, dividend_rate, time_step):
return 0.5*time_step*(vol**2*asset_step**2 - (risk_free_rate - dividend_rate)*asset_step)

def B(vol, asset_step, risk_free_rate, time_step):
return 1 - time_step*(vol**2*asset_step**2 + risk_free_rate)

def C(vol, asset_step, risk_free_rate, dividend_rate, time_step):
return 0.5*time_step*(vol**2*asset_step**2 + (risk_free_rate - dividend_rate)*asset_step)

def call_payoff(asset_price, strike_price):
return max(asset_price - strike_price, 0)

S_max          = 140.
strike         = 100.
expiry         = 2.
asset_step     = 1
risk_free_rate = 0.05
dividend_rate  = 0.0
vol            = 0.2
t_step         = 0.001

no_asset_steps = int(S_max/asset_step)
print "Number of asset steps: ", no_asset_steps

no_time_steps = int(expiry/t_step)
print "Number of time steps: ", no_time_steps

S = np.zeros(no_asset_steps+1)
V = np.zeros((no_asset_steps+1, no_time_steps+1))

for i in xrange(no_asset_steps+1):
S[i] = i*asset_step
V[i,no_time_steps] = call_payoff(S[i], strike)

for k in range(no_time_steps, 0, -1):
for i in xrange(1,no_asset_steps):
A_co = A(vol, i, risk_free_rate, dividend_rate, t_step)
B_co = B(vol, i, risk_free_rate, t_step)
C_co = C(vol, i, risk_free_rate, dividend_rate, t_step)
V[i, k-1] = A_co*V[i-1, k] + B_co*V[i,k] + C_co*V[i+1,k]

V[0, k-1] = 0
V[no_asset_steps, k-1] = 2*V[no_asset_steps-1, k-1] - V[no_asset_steps-2, k-1]

• For convergence of a parabolic equation the step sizes in x and t cannot be arbitrary but must be chosen to satisfy the Friederichs Lewy condition Oct 16, 2015 at 19:30
• And what is the condition in this case? I tested against the stability condition given in Wilmott dt < 1/(i^2 vol^2) which is correct correct in this case. Oct 16, 2015 at 19:46
• FWIW you can see one form of the condition on Page 14 here people.math.gatech.edu/~meyer/MA6635/chap4.pdf I am not sure which condition applies to you. Oct 16, 2015 at 20:04
• Turns out my code is actually correct. I was comparing against the wrong analytical values. The stability issue is certainly still there, but this seems to be a feature of the explicit method rather than my code. Thanks for the links! Oct 17, 2015 at 13:24