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:
- Can both models be used with mu and sigma based on log calculations?
- Can the time to maturity T be entered in the same unit for both models?
- 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)```