# Noob Question - Monte Carlo vs BAPM European Option Pricing Discrepancy

It's winter break (happy new year!), and I'm trying implement a few options pricing models (bapm, tapm, monte carlo, Fast Fourier etc.) for practice.

The issue: My BAPM CRR model converges to 8.45544504853379 for "ec", which is consistent with online results. My Monte Carlo implementation converges to 4.748831589623564 :(. I can't find the error in my code...I was hoping someone could tell me if such a large error is possible when using ~100k GBM simulations.

type = ["ec", "ep", "ac", "ap"] (Euro call, euro put, etc.)

Binomial Asset Pricing Model

# Initial parameters
S0 = 100      # initial stock price
K = 100       # strike price
T = 1         # time to maturity in years
r = 0.06      # annual risk-free rate
N = 100         # number of time steps
sig = .13      # sigma - std of stock annualized

import numpy as np
import matplotlib.pyplot as plt

class BAPM():
def __init__(self, S0, K, T, N, r, c, sig) -> None:
self.S0 = S0
self.K = K
self.T = T
self.N = N
self.r = r
self.c = c
self.sig = sig
self.dt = self.T/self.N

def JarrowRudd(self, type):
u = np.exp((self.r - self.sig**2 / 2)*self.dt + self.sig*np.sqrt(self.dt))
d = np.exp((self.r - self.sig**2 / 2)*self.dt - self.sig*np.sqrt(self.dt))
p = (np.exp(self.r*self.dt) - d)/(u - d)
q = 1-p
disc = np.exp(-self.r*self.dt)
return self.optionPrice(type, p, q, u, d, disc)

def CRR(self, type, drift = 0):
u = np.exp(drift*self.dt + self.sig * np.sqrt(self.dt))
d = np.exp(drift*self.dt - self.sig * np.sqrt(self.dt))
p = (np.exp(self.r*self.dt) - d)/(u - d)
q = 1-p
disc = np.exp(-self.r*self.dt)
return self.optionPrice(type, p, q, u, d, disc)

def Tian(self, type):
v = np.exp(self.sig **2 * self.dt)
u = .5 * np.exp(self.r * self.dt) * v * (v + 1 + np.sqrt(v**2 + 2*v -3))
d = .5 * np.exp(self.r * self.dt) * v * (v + 1 - np.sqrt(v**2 + 2*v -3))
p = (np.exp(self.r * self.dt) - d)/(u - d)
q = 1-p
disc = np.exp(-self.r * self.dt)
return self.optionPrice(type, p, q, u, d, disc)

def optionPrice(self, type, p, q, u, d, disc):
# European Option
if type == "ec" or type == "ep":
# terminal value of stock
S = self.S0 * (d ** (np.arange(self.N,-1,-1))) * (u ** (np.arange(0,self.N+1,1)))
V = 0
# value of option at terminal time
if type == "ec":
V = np.maximum(S - self.K, np.zeros(self.N+1))
elif type == "ep":
V = np.maximum(self.K - S, np.zeros(self.N+1))
# take 2 at a time and get weighted-discounted value
for i in np.arange(self.N,0,-1):
V = disc * ( p * V[1:i+1] + q * V[0:i])
# final result is option price
return V[0]

elif type == "ac" or type == "ap":
print("test successful")
# terminal value of stock
S = self.S0 * (d**(np.arange(self.N,-1,-1))) * (u ** (np.arange(0,self.N+1,1)))

V = 0
# value of payoff at terminal time
if type == "ac":
V = np.maximum(0, S - self.K)
elif type == "ap":
V = np.maximum(0, self.K - S)

# backtrack through options tree
for i in np.arange(self.N-1, -1, -1):
# recalculate prices at current level
S = self.S0 * d**(np.arange(i,-1,-1)) * u**(np.arange(0,i+1,1))

V[:i+1] = disc * ( p*V[1:i+2] + q*V[0:i+1] )
V = V[:-1]

if type == 'ap':
V = np.maximum(V, self.K - S)
else:
V = np.maximum(V, S - self.K)

return V[0]

bapm = BAPM(S0, K, T, N, r, c, sig)
print(bapm.CRR(type = "ec"))


8.45544504853379

Monte Carlo Pricing (vanilla only...American not implemented yet)

import numpy as np
import matplotlib.pyplot as plt

class MCPricing():
def __init__(self, S0, drift, sig, T, tIncrement, r) -> None:
self.S0 = S0
self.mu = drift
self.sig = sig
self.T = T
self.dt = tIncrement
self.r = r

self.N = round(self.T / self.dt)
self.t = np.linspace(0, self.T, self.N)

def simulateGBM(self):
W = np.random.standard_normal(size = self.N)
W = np.cumsum(W)*np.sqrt(self.dt)
X = (self.mu - .5*self.sig**2)*self.t + self.sig*W
S = self.S0*np.exp(X)
return S.tolist()

def calculateOptionPrice(self, K, nsim, type):
if type == "ec" or type == "ep":
sims = [self.simulateGBM() for i in range(nsim)]
payoffs = []
for path in sims:
if type == "ec":
payoffs.append(np.maximum(path[-1] - K, 0))
elif type == "ep":
payoffs.append(np.maximum(K - path[-1], 0))
payoffs = np.array(payoffs)
return np.average(np.exp(-1*self.r * self.T) * payoffs)

mcpm = MCPricing(S0 = 100, drift=0, sig =.13, T=1, tIncrement=.01, r=.06)
path1 = mcpm.simulateGBM()
plt.plot(mcpm.t, path1)

mcpm.calculateOptionPrice(100, 100000, "ec")


4.748831589623564

I'd appreciate any help! Thanks

• Your mc drift is scaled by self.t, whereas I think dt would be more appropriate, given the sqrt dt scaling for the Brownian Dec 31, 2021 at 21:43
• @JamesSpencer-Lavan I made that change, and the MC is now consistent at \$5.30. This is much better but is still far off from the BAPM. I tried decreasing the dt but get the same result. Dec 31, 2021 at 22:23
• Your drift term is wrong. We're pricing the option in a risk-neutral world meaning that the drift needs to be equal to r, not 0. Try setting it to .06 and you will probably get the correct answer. Jan 1 at 21:19

Really wish I knew why my previous version didn't work...but I found some modifications to GBM that got me the answer.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

class MCPricing():
def __init__(self, S0, drift, sig, T, tIncrement, r) -> None:
self.S0 = S0
self.mu = drift
self.sig = sig
self.T = T
self.dt = tIncrement
self.r = r

self.N = round(self.T / self.dt)
self.t = np.linspace(0, self.T, self.N)

def simulateGBM(self, nsims):
Xt = np.log(self.S0) +\
np.cumsum(((self.mu - self.sig**2/2)*self.dt + \
self.sig*np.sqrt(self.dt) * \
np.random.normal(size=(self.N,nsims))),axis=0)

return np.exp(Xt)

def calculateOptionPrice(self, K, nsims, type):
if type == "ec" or type == "ep":
paths = self.simulateGBM(nsims)
if type == "ec":
payoffs = np.maximum(paths[-1]-K, 0)
elif type == "ep":
payoffs = np.maximum(K - paths[-1], 0)

option_price = np.mean(payoffs)*np.exp(-self.r * self.T) #discounting back to present value

return option_price

mcpm = MCPricing(S0 = 100, drift=0.06, sig =.13, T=1, tIncrement=.001, r=.06)
paths = mcpm.simulateGBM(1000)
plt.plot(paths)
mcpm.calculateOptionPrice(100, 10000, "ec")


8.373535799654391

• It works now because you set drift and r to be equal, which is exactly what I told you to do in a comment a few days ago. Jan 5 at 21:31
• @Oscar yeah I made that change in the previous implementation but was still getting an incorrect answer. But yeah, the mismatched rates was a stupid mistake. Thank you for the help. Jan 7 at 7:41