Diagnosing a Binary Option Finite Difference Pricing Error

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


# 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[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
```

## Answer by StackG (score 3, accepted)

https://quant.stackexchange.com/a/57747

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()))
```