Simulating Risk-Neutral GBM to Estimate a Black–Scholes Call Price

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


# 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?

## Answer by wgajate (score 1)

https://quant.stackexchange.com/a/58824

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)
```