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?

  • 1
    $\begingroup$ 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). $\endgroup$
    – nbbo2
    Oct 16, 2020 at 9:21
  • $\begingroup$ @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? $\endgroup$
    – John Smith
    Oct 16, 2020 at 9:37
  • $\begingroup$ 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. $\endgroup$
    – nbbo2
    Oct 16, 2020 at 10:46
  • $\begingroup$ @JohnSmith have you had a chance to try out the code I posted in the answer below? $\endgroup$
    – wgajate
    Oct 31, 2020 at 17:17

1 Answer 1


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):
    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
5.795249085244208  # results vary from run to run (random seed not set)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.