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