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

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 datetime import date, timedelta
#import datetime
import yfinance as yf
import scipy.stats as si

# Get stock price data
# Calculate log annual returns (mu) and log volatility (sigma)
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

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


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

• Wow...! What an absolutely awesome answer. This is a much better, clearer and more helpful answer than I expected. Thanks so much. N.B. (1) I use the 13 Week Treasury Bill (^IRX) rate for the risk-free rate, for which I used 0.135, But guess this must be annualised. (2) For implied volatility, can you use the actually quoted implied volatility for the option under consideration, from Yahoo for example? Or should I calculate it myself. – twhale Jun 26 at 16:33
• You write: "when using the GBM model for pricing options: you should switch from the real world probability measure to the risk-neutral measure." Since I am using the risk-free rate for r (=0.135, which I guess I forgot to annualise to 0.54), this is already done in the code, right? – twhale Jun 26 at 16:42
• If your option maturity is 1-year, then using a 13-month T-bill yield is a reasonably good proxy. For the implied vol, using Yahoo finance implied vols is a good idea: as long as the implied vol corresponds to the maturity of the option you're trying to price. PS: yes, you use the risk-free rate which is the equivalent of switching to the risk-neutral measure (strictly speaking under the risk-neutral measure, you need to use the risk-free rate as well as the implied volatility, i.e. not historical vol and not historical mu). PPs: the T-bill yield should already be quoted as annualized. – Jan Stuller Jun 26 at 16:51
• Got it. And then instead of using the sigma calculated based on sigma = np.std(apple.log_price), I should switch to sigma = implied volatility` – twhale Jun 26 at 17:15
• Exactly, you got it. – Jan Stuller Jun 26 at 17:16