I have recently started learning about option pricing and the Black Scholes formula, where stock prices are assumed to be lognormally distributed and returns normally distributed. While trying to do some simulations to learn more about the parameters and their behaviour I noticed that I seem to be misapplying the concept of the stock prices being lognormal. Please see the code below in python where I simulate a 30 day call option with daily price changes, where if I try to use lognormal prices I converge on about 0.4~ EV versus breakeven, which is expected, if I just use regular normally distributed returns, which would allow for negative stock prices.

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

def black_scholes(S, K, r, T, sigma, option="call"):
    d1, d2 = d(S, K, r, T, sigma)

    if option == "call":
        return (S * stats.norm.cdf(d1, 0, 1) - K * np.exp(-r * T) * stats.norm.cdf(d2, 0, 1))
    elif option == "put":
        return (K * np.exp(-r * T) * stats.norm.cdf(-d2, 0, 1) - S * stats.norm.cdf(-d1, 0, 1))

def d(S, K, r, T, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return d1, d2

price = 100
strike = 105
r = 0
T = 30/365
sigma = 0.5
option = "call"

cost = black_scholes(price,strike,r,T,sigma,option) # = 3.69...

avg_gains = []
avg_gains_logn = []
for _ in range(1000):
    runs = {n:None for n in range(1000)}
    runs_logn = {n:None for n in range(1000)}
    for run in runs:
        p = price
        p_logn = price
        for _ in range(30):
            c = np.random.normal(0,sigma/np.sqrt(365))
            p *= 1+c
            p_logn *= np.exp(c)
        runs[run] = p
        runs_logn[run] = p_logn

    gains = []
    gains_logn = []
    for n, result in runs.items():
        gains.append(max(0, result - strike) - cost)
        gains_logn.append(max(0, runs_logn[n] - strike) - cost)


print(sum(avg_gains)/1000, sum(avg_gains_logn)/1000)

enter image description here


Under the BSM model, the terminal stock price is assumed to be lognormally distributed, with expected value equal to $\mathrm{E}_t(S_T)=S_0e^{r(T-t)}$. In order to achieve this in your simulation (of your log-normal stock process), you may want to modify your code:

c = np.random.normal(r-0.5*sigma**2/365, sigma/np.sqrt(365)) # instead of mean=0

This ensures that the drift induced by the exponentiation of the logarithmic return is cancelled, i.e.

$$ \mathrm{E}\left(e^{\left(r-\frac{1}{2}\sigma^2\right)(T-t)+\sigma\sqrt{(T-t)}z_{0,1}}\right)=e^{r(T-t)} $$

In your case, $r=0$, so only the correction for the volatility is required during simulation.

| improve this answer | |
  • $\begingroup$ Thanks Kermit, that should be ? c = np.random.normal(r-0.5*sigma**2/365, sigma/np.sqrt(365)) (sigma**2 instead of sigma**sigma) $\endgroup$ – herbertmm Mar 27 at 17:59
  • $\begingroup$ Yes thanks for pointing that out $\endgroup$ – Kermittfrog Mar 27 at 20:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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