How to best predict option prices using Brownian motion and compare it to the Black and Scholes model?

Question

I am trying to use Brownian motion to predict option prices and compare the outcomes to Black and Scholes. For this purpose, I would like to calculate the average returns (mu) and volatility (sigma) of the underlying asset based on continuous compounding - hence I use log functions.

However, I think that there could be several mistakes in my approach which I am unable to confirm. Unfortunately, I cannot find clear answers to my questions on the web, on sites such as quantconnect, as well as this forum. These are my doubts:

1. Can both models be used with mu and sigma based on log calculations?
2. Can the time to maturity T be entered in the same unit for both models?
3. When using Geometric Brownian motion for simulating stock prices, we loop say 10,000 times and take the average outcome. Most GBM models for option prices I find on the web don't seem to loop. Why is that?

My Python code is below. Any suggestions?

from math import log, e
from pandas_datareader import data
from datetime import date, timedelta
#import datetime
import yfinance as yf
import scipy.stats as si

# Get stock price data
apple = data.DataReader('AAPL', 'yahoo', '2018/1/1')
spot = apple["Adj Close"][-1]
# Calculate log annual returns (mu) and log volatility (sigma)
apple['log_price'] = np.log(apple['Adj Close'])
apple['log_return'] = apple['log_price'].diff()

mu = apple.log_return.sum()/apple.log_return.count()
mu = mu*365 + 0.5*apple.log_return.var()*np.sqrt(365)
sigma = np.std(apple.log_price)   #sigma: volatility of underlying

spot = 463.94    #spot: spot price
K = 460    # strike price
T = 1   # time to maturity
r = 0.135   # risk-free interest rate

# Black and Scholes calculation
s0 = spot
def euro_vanilla_call(S, K, T, r, sigma):    

    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    
    call = (S * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))
    
    return call

BSM1 = euro_vanilla_call(spot, K, T, r, sigma)
print(BSM1)

# Monte Carlo simulation
def mc_euro_options(option_type,s0,strike,maturity,r,sigma,num_reps):
    payoff_sum = 0
    for j in range(num_reps):
        st = s0
        st = st*e**((r-0.5*sigma**2)*maturity + sigma*np.sqrt(maturity)*np.random.normal(0, 1))
        if option_type == 'c':
            payoff = max(0,st-strike)
        elif option_type == 'p':
            payoff = max(0,strike-st)
        payoff_sum += payoff
    premium = (payoff_sum/float(num_reps))*e**(-r*maturity)
    return premium

MCP1 = mc_euro_options('c', spot, K, T, r, sigma, 100)
print(MCP1)```

Answer 1

The GBM model can be written as:

$$ \delta S_t= \mu S_t \delta t+\sigma S_t\delta t $$

The above is short-hand for the following SDE:

$$ S(t)=S(0)+\int^{t}_{0}\mu S(h)dh+\int^{t}_{0}\sigma S(h)dW(h) $$

Solving the above SDE yields an expression that you implemented in your code:

$$ S(t)=S_0exp\left((\mu-0.5 \sigma^2)t+\sigma \sqrt{t} Z\right) $$

The Black-Scholes formula can be derived directly by applying the option pay-off to the above solution of the SDE (below I use the real-world measure for simplicity*, see asterix note further below in the text for more details):

$$ Call(t_0)=e^{-rt}\mathbb{E}\left[ (S_t-K)I_{ \left( S_t>K \right) } \right] = \\ = e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }-KI_{ \left( S_t>K \right) } \right]=\\=e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }\right]-e^{-rt}K\mathbb{E}\left[ I_{ \left( S_t>K \right) }\right]$$

Focusing on the second term:

$$ e^{-rt}K\mathbb{E}\left[ I_{ \left( S_t>K \right) }\right] = e^{-rt}K\mathbb{P}\left( S_t>K \right) = \\ = e^{-rt}K\mathbb{P}\left( S_0 exp\left((\mu-0.5 \sigma^2)t+\sigma \sqrt{t} Z\right)>K \right) = \\ = e^{-rt}K\mathbb{P}\left( (\mu-0.5 \sigma^2)t+\sigma \sqrt{t} Z>ln \left(\frac{K}{S_0} \right) \right) = \\ = e^{-rt}K\mathbb{P}\left( Z>\frac{ln \left(\frac{K}{S_0} \right) -\mu t + 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}K\mathbb{P}\left( Z> (-1)\frac{ln \left(\frac{S_0}{K} \right) +\mu t - 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}K\mathbb{P}\left( Z \leq \frac{ln \left(\frac{S_0}{K} \right) +\mu t - 0.5\sigma^2 t}{\sigma \sqrt{t} } \right) = \\ = e^{-rt}KN(d_2) $$

The first term $e^{-rt}\mathbb{E}\left[ S_tI_{ \left( S_t>K \right) }\right]$ requires a tiny little bit more work to evaluate, but using a similar technique this term comes out as $S_0N(d_1)$.

So what this tedious usage of formulas was meant to demonstrate is that the Black-Scholes formula can be shown to be a direct consequence of the GBM model for the underlying stock price: therefore this answers your first and second questions:

(i) Yes, the mu and sigma in both models are identical, because the BS formula is based on the GBM model

(ii) Yes, both models need to be consistent with one another in terms of units of time.

*Word of warning: there is one additional step that needs to be performed when using the GBM model for pricing options: you should switch from the real world probability measure to the risk-neutral measure. In practical terms, it means that your drift $\mu$ needs to be replaced with drift $r$, where $r$ should be the "risk-free" rate corresponding to the option maturity. If you don't have access to the entire OIS curve for USD, then I would just take the FED funds rate as a proxy for $r$ (right now, the FED funds rate is 0.25%).

You should also use implied volatility to price the option. But using historical volatility (as you do in your code) as a proxy is ok if you just want to experiment.

Your third question: if you want to price the option by Monte-Carlo (i.e. simulating stock price first, then taking expectation of the option pay-off at maturity), you need to run "n" simulations (i.e. loops). But because you know the analytical solution to the GBM model as shown above and you can plug this directly into the option pay-off and analytically compute the option price that way, you don't actually need to run an MC simulation. You can just price the option directly via the B-S formula.

It's basically up to you if you want to evaluate the expectation in the Option pay-off formula via Monte-Carlo or analytically (which leads to the BS formula directly). Obviously, analytical evaluation is more accurate than numerical approximation. Running a numerical simulation on a problem which you know how to solve analytically is a bit like hiding your own Easter eggs and then searching for them.

PS: last but not least, you should not use 365 days, but rather 260 days per year (because there are only roughly 260 trading days in a calendar year).