# Why am I struggling to replicate the Black-Scholes price of an option stochastically?

I am currently trying to replicate the Black-Scholes price of a call option using stochastic simulations of the price moves of the underlying. My code is as follows:

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import stats

class MarketReturnModel(object):
def __init__(self, current_price, drift, volatility, steps_py):
self.current_price = current_price
self.drift = drift
self.volatility = volatility
self.steps_py = steps_py

def simulate(self, steps):
rand_norm = np.random.normal(0, 1 / self.steps_py, steps)
s = [self.current_price]
for t in range(0, steps):
ds_t = self.drift * s[t] * (1 / self.steps_py) + self.volatility * s[t] * rand_norm[t]
s.append(s[t] + ds_t)
return s

if __name__ == '__main__':

current_price = 120
strike = 135
r = 0.05
volatility = 0.3
T = 0.5

d1 = (math.log(current_price / strike) + (r + 0.5 * volatility ** 2) * T) / (volatility * math.sqrt(T))
d2 = d1 - volatility * math.sqrt(T)
bs_call_price = current_price * stats.norm(0, 1).cdf(d1) - strike * math.exp(-r * T) * stats.norm(0, 1).cdf(d2)

print(f'BS Price: {bs_call_price}')

payoffs = []

for i in range(0, 100000):
mr = MarketReturnModel(current_price, r, volatility, 365)
sim = mr.simulate(365)
payoffs.append(max(sim[-1] - strike, 0))
print(f'Stochastic Price: {math.exp(-r * T) * sum(payoffs) / len(payoffs)}')


This returns the output

BS Price: 5.590757390745406
Stochastic Price: 0.0


I've taken a look at my stock price simulations and the never reach 135 (i.e. the option strike), which is why the stochastic option price is 0. I'm therefore guessing that there is something wrong with my stock price simulations. Can anyone see what this error is?

• Your simulation of GBM is not correct. There should be a $\sigma \varepsilon \sqrt{\Delta t}$ term where $\varepsilon$ is the random normal term, but I don't see the $\sqrt{\Delta t}$. investopedia.com/articles/07/montecarlo.asp (Also the code seems more complex than necessary). – noob2 Oct 16 '20 at 9:21
• @noob2 Thanks for the reply. I used the following SDE to construct ds_t: $dS_t = \mu S_t dt + \sigma S_t dW_t$ Is my calculation not consistent with this? – John Smith Oct 16 '20 at 9:37
• Right. $dW_t$ expands to $\varepsilon \sqrt{\Delta t}$ in discrete time, where $\varepsilon$ is a normal N(0,1). Look at examples of code for simulating GBM, there are many examples online, for example en.wikipedia.org/wiki/Geometric_Brownian_motion Note the np.sqrt(dt) I don't see that in your code. – noob2 Oct 16 '20 at 10:46
• @JohnSmith have you had a chance to try out the code I posted in the answer below? – wgajate Oct 31 '20 at 17:17

Try the code below. It consumes less memory and takes about 6 minutes to run.

Please note that when simulating the stock price for purposes of pricing derivatives, we use the risk-neutral stock price process. The drift of the risk-neutral process becomes the risk free rate minus half the square of annualized volatility as indicated in the code.

import numpy as np
import pandas as pd

T = .5 # half a year
number_time_steps = 365  #evenly spaced time steps
dt = T/number_time_steps
sigma = .3 # percent per year
S_0 = 120 #intial stock price
r = .05 # determinstic growth per year
K = 135  #call option strike price

def price_call_with_stock_sim(num_trials):
avg_call_payoff=0
for idx in range(1,num_trials+1):
df = pd.DataFrame(data={"dW":np.random.normal(0, 1, number_time_steps)*np.sqrt(dt), "dt":np.ones(number_time_steps)*dt})
df["t"] = df["dt"].cumsum()
df["W"] = df["dW"].cumsum()
df["exp((r-.5s^2)t+sW)"] = np.exp((r-.5*(sigma**2))*df["t"]+sigma*df["W"])
df["S(t)"] = S_0*df["exp((r-.5s^2)t+sW)"]
S_T = df["S(t)"].values[-1]
call_payoff = max(S_T-K,0)
avg_call_payoff = (avg_call_payoff*(idx-1) + call_payoff)/idx
return avg_call_payoff*np.exp(-r*T) #discount payoff back to t=0

num_trials = 100000
price_call_with_stock_sim(num_trials)
5.795249085244208  # results vary from run to run (random seed not set)