I transformed Blacks Scholes equation to a Heat equation. I try to use explicit finite difference method to solve this PDE and get the price of a call option. I also solve for this by using black schols equation "analytically".
The problem is that I cannot get more accurate in the numerical result. Here is my Python code.
Here is the note of my algorithm: https://drive.google.com/file/d/0B5h3oewtgjFgdVFpNFJRNTB5LXM/view?usp=sharing
import math
import numpy as np
from scipy.stats import norm
s0 = 15
sigma = 0.2
r = 0.01
t = 1
Xmax = 10
'''B-S price'''
def C(s,k,t):
d1 = (math.log(s/k)+(r+sigma*sigma/2)*t)/(sigma*math.sqrt(t))
d2 = (math.log(s/k)+(r-sigma*sigma/2)*t)/(sigma*math.sqrt(t))
return s*norm.cdf(d1)-math.exp(-r*t)*k*norm.cdf(d2)
print('B-S', C(s0,10,t))
'''Explicit_finite_difference'''
EFD_n_x = 500
EFD_n_t = 100
EFD_k = Xmax/EFD_n_x
EFD_h = t/EFD_n_t
EFD_xx = np.linspace(Xmax,-Xmax, 2 * EFD_n_x + 1)
EFD_xx = EFD_xx[1:2 * EFD_n_x]
def EFD_T0_Bound(x):
return max(math.exp(x)-10*math.exp(-r*t),0)
def EFD_U_Bound(tao):
return math.exp(Xmax)-10*math.exp(-r*(t-tao))
def EFD_L_Bound(tao):
return 0
EFD_T0bound = np.vectorize(EFD_T0_Bound)
EFD_lambda = EFD_h*sigma*sigma/2/EFD_k/EFD_k
EFD_A = (np.eye(2 * EFD_n_x - 1) * (1-2*EFD_lambda)
+ np.eye(2 * EFD_n_x - 1, k=1)*EFD_lambda
+ np.eye(2 * EFD_n_x - 1, k=-1)*EFD_lambda)
EFD_Y = np.zeros(2 * EFD_n_x - 1)
EFD_U = EFD_T0bound(EFD_xx)
for i in range(EFD_n_t):
EFD_Y[0] = EFD_lambda*EFD_U_Bound(EFD_h*i)
EFD_Y[2 * EFD_n_x - 2] = EFD_lambda*EFD_L_Bound(EFD_h*i)
EFD_U = np.dot(EFD_A,EFD_U) + EFD_Y #U_t_i+1 = A * U_t_i + Y
print('Explicit_finite_difference',EFD_U[EFD_n_x - 1 - round(math.log(s0)/EFD_k)])
in my opinion, since you are using an explicit scheme (that is not unconditionnaly stable), you cant set your nb of asset step and your nb of time step indepently. usually people link them together with volatility in order to keep the scheme stable