Diagnosing Monte Carlo European Option Pricing Errors
Article Quant Q&A · Author: NotYoAvgJoe
Summary
The document compares a binomial tree price for a European call with a much lower estimate from a Monte Carlo implementation. The key diagnostic is the simulated stock process: the original code builds paths with time points and Brownian increments that do not align correctly, and it uses a drift inconsistent with risk-neutral pricing. For option valuation under the stated risk-free-rate framework, the simulation should use the risk-free rate as the drift and apply time-scaled increments.
A revised implementation simulates log prices using cumulative normally distributed increments, computes terminal call or put payoffs, and discounts their average at the risk-free rate. Its reported call estimate is close to the binomial result. This example illustrates how simulation setup can dominate sampling noise; it does not establish a general convergence guarantee or provide a formal error analysis. The code also does not address American exercise, and model assumptions and parameter choices remain relevant to any comparison.
Key ideas
- Monte Carlo option valuation should simulate under the risk-neutral measure when pricing from the risk-free rate.
- Brownian increments must be scaled by the square root of the time step.
- European option prices can be estimated by discounting the average simulated terminal payoff.
- A discrepancy between a tree and simulation estimate can indicate a path construction or drift error, not just sampling noise.
- The example covers vanilla European options and does not implement early exercise.
Tags
Full text
# Noob Question - Monte Carlo vs BAPM European Option Pricing Discrepancy
# 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
## Answer by NotYoAvgJoe (score 0, accepted)
https://quant.stackexchange.com/a/69354
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.373535799654391Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.