I am trying to price the Binary option using Explicit Finite Difference Method. However, the output is not matching with the closed form solution formula.
Here is the code for the same:
import numpy as np
import math
import scipy.stats as si
# set up parameters
S0 = 50
K = 40
r = 0.01
T = 0.5
sigma = 0.2
Smax = 100
M = 100 # S
N = 1000 # t
is_call = True
def ExplicitFiniteDifferences(S0, K, r, T, sigma, Smax, M, N, is_call):
""" Shared attributes and functions of FD """
M, N = int(M), int(N) # Ensure M&N are integers
dS = Smax / float(M)
dt = T / float(N)
iValues = np.arange(1, M)
jValues = np.arange(N)
grid = np.zeros(shape=(M+1, N+1)) # grid is M+1 by N+1
SValues = np.linspace(0, Smax, M+1)
alpha = 0.5*dt * (sigma**2 * iValues**2 - r*iValues)
beta = -dt * (sigma**2 *iValues**2 + r)
gamma = 0.5*dt * (sigma**2 *iValues**2 + r*iValues)
coeffs = np.diag(alpha[1:], -1) + np.diag(1 + beta) + np.diag(gamma[:-1], 1)
# terminal condition
if (S0 > K):
grid[:, -1] = 1
else:
grid[:, -1] = 0
# side boundary conditions
coeffs[0,0] += 2*alpha[0]
coeffs[0,1] -= alpha[0]
coeffs[-1,-1] += 2*gamma[-1]
coeffs[-1,-2] -= gamma[-1]
for j in reversed(jValues):
grid[1:-1, j] = np.dot(coeffs, grid[1:-1, j+1])
grid[0, j] = 2 * grid[1, j] - grid[2, j]
grid[-1, j] = 2 * grid[-2, j] - grid[-3, j]
return np.interp(S0,SValues,grid[:, 0])
x = ExplicitFiniteDifferences(S0, K, r, T, sigma, Smax, M, N, is_call)
print("The price of Binary Option using Explicit Finite Difference Method is ", x*math.exp(-r*T))
d1 = (math.log(S0/K) + (r + (sigma**2)*0.5)*T)/(sigma*math.sqrt(T))
d2 = d1 - sigma*math.sqrt(T)
Nd2 = si.norm.cdf(d2)
print("The price of Binary Option using closed form solution is ",math.exp(-r * T)*Nd2)
Output:
The price of Binary Option using Explicit Finite Difference Method is 0.990049821373499
The price of Binary Option using closed form solution is 0.9338439709795566