Price of Binary Option using Explicit Finite Difference Method not matching with closed Form solution

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
coeffs[0,1] -= alpha
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
• Does your pde get the forward correct?
– will
Sep 2 '20 at 13:58

The price close to 0.93 is correct, here is a reimplementation of both FD and analytic using QuantLib:

import QuantLib as ql

# World State for Vanilla Pricing
spot = 50
vol = 0.2
rate = 0.01
dividend = 0.0

today = ql.Date(1, 9, 2020)

day_count = ql.Actual365Fixed()
calendar = ql.NullCalendar()

# Set up the vol and risk-free curves
volatility = ql.BlackConstantVol(today, calendar, vol, day_count)
riskFreeCurve = ql.FlatForward(today, rate, day_count)

flat_ts = ql.YieldTermStructureHandle(riskFreeCurve)
dividend_ts = ql.YieldTermStructureHandle(riskFreeCurve)
flat_vol = ql.BlackVolTermStructureHandle(volatility)

process = ql.BlackScholesMertonProcess(ql.QuoteHandle(ql.SimpleQuote(spot)), dividend_ts, flat_ts, flat_vol)

# And define the option
expiry_date = ql.Date(1, 3, 2021)
strike = 40
payoff = ql.CashOrNothingPayoff(ql.Option.Call, strike, 1)
european_exercise = ql.EuropeanExercise(expiry_date)

binary_option = ql.VanillaOption(payoff, european_exercise)

# Run with Analytic Engine
engine = ql.AnalyticEuropeanEngine(process)
binary_option.setPricingEngine(engine)
print("Analytic Price: {}".format(binary_option.NPV()))

# Run with FD Engine
tGrid, xGrid = 2000, 200
engine = ql.FdBlackScholesVanillaEngine(process, tGrid, xGrid)
binary_option.setPricingEngine(engine)
print("Finite Differences Price: {}".format(binary_option.NPV())) • Perfect! This is what I was looking for. Just 1 quick query. The FdBlackScholesVanillaEngine in QuantLib is based on explicit or implicit Finite Difference method? Sep 3 '20 at 4:30