I am trying to run a Monte Carlo simulation to help me price option greeks in python. Most of the values I am getting are correct, but I am not sure why I never get the right result for Theta. Plus, the estimated values for Theta keep changing but it's never around the correct answer.
Correct answer: {'Delta': -0.3611, 'Gamma': 0.0177, 'Theta': -0.745, 'Vega': 26.483, 'Rho': -49.635}
What my code gives me: {'Delta': -0.36114866036456306, 'Gamma': 0.01765540297071766, 'Theta': 4.969082186799323, 'Vega': 26.482584628353578, 'Rho': -49.649267741554404}
import numpy as np
from scipy.stats import norm
def monte_carlo_option_greeks(option_type, S, K, T, r, sigma, steps, N, dS=0.01, dSigma=0.01, dT=0.01):
"""
Calculate option Greeks using Monte Carlo simulation.
Inputs:
- option_type: 'call' or 'put'
- S: Current stock price
- K: Strike price
- T: Time to maturity (in years)
- r: Risk-free interest rate
- sigma: Volatility
- steps: Time steps for the simulation
- N: Number of simulation trials
- dS: Increment for calculating Delta (optional, default is 0.01)
- dSigma: Increment for calculating Vega (optional, default is 0.01)
- dT: Increment for calculating Theta (optional, default is 0.01)
Returns a dictionary containing the estimated Greeks: {'Delta', 'Gamma', 'Theta', 'Vega', 'Rho'}
"""
def geo_paths(S, T, r, sigma, steps, N):
dt = T / steps
ST = np.cumprod(np.exp((r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) *
np.random.normal(size=(steps, N))), axis=0) * S
return ST
def black_scholes_price(S, K, T, r, sigma):
d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type == 'call':
price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else: # option_type == 'put'
price = K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
return price
# Simulate underlying asset price paths
paths = geo_paths(S, T, r, sigma, steps, N)
# Calculate option payoffs
if option_type == 'call':
payoffs = np.maximum(paths[-1] - K, 0)
else: # option_type == 'put'
payoffs = np.maximum(K - paths[-1], 0)
# Discount to present value
option_price = np.mean(payoffs) * np.exp(-r * T)
# Calculate Greeks
delta = (black_scholes_price(S + dS, K, T, r, sigma) - black_scholes_price(S - dS, K, T, r, sigma)) / (2 * dS)
# Corrected gamma calculation
gamma = (black_scholes_price(S + dS, K, T, r, sigma) - 2 * black_scholes_price(S, K, T, r, sigma) + black_scholes_price(S - dS, K, T, r, sigma)) / (dS**2)
# Corrected theta calculation with negative sign for a European put option
dT_theta = dT / 365 # Convert to a daily increment for theta calculation
theta = -(black_scholes_price(S, K, T - dT_theta, r, sigma) - option_price) / dT
vega = (black_scholes_price(S, K, T, r, sigma + dSigma) - black_scholes_price(S, K, T, r, sigma - dSigma)) / (2 * dSigma)
rho = (black_scholes_price(S, K, T, r + dT, sigma) - black_scholes_price(S, K, T, r - dT, sigma)) / (2 * dT)
greeks = {
'Delta': delta,
'Gamma': gamma,
'Theta': theta,
'Vega': vega,
'Rho': rho
}
return greeks
# Example usage:
option_type = 'put' # 'call' or 'put'
S = 50 # stock price S_{0}
K = 52 # strike
T = 2 # time to maturity (in years)
r = 0.05 # risk-free interest rate
sigma = 0.3 # annual volatility in decimal form
steps = 100 # time steps
N = 10000 # number of simulation trials
greeks = monte_carlo_option_greeks(option_type, S, K, T, r, sigma, steps, N)
print(greeks)
I was initially getting very large values for Theta, so I used a larger increment (dT = 0.001 to dT = 0.01) to get more sensible numbers. I am still unsure of the actual formula though.
If for the calculation of theta as an offset to T the amount dT_theta = 0.01 / 365 is used, i.e. only 2.739726027e-05 years (as that is the unit which black_scholes_price takes for its argument T, which is varied for the theta calculation), i.e. a relative difference in T of only 1.3698630135e-05, then my guess is that the precision (floating point arithmetic) starts to become an issue, as the price differences are too small. This hypothesis is also supported by the huge variance of theta one gets when repeatedly running monte_carlo_option_greeks (where the other greeks have no to very small variance):
print(np.var([monte_carlo_option_greeks(
option_type, S, K, T, r, sigma, steps, N
)['Theta'] for k in range(100)]))
prints
65.54872356722292
I also recommend sticking to the central method (i.e. dividing by double the delta of the varying parameter and offsetting in both directions) for calculating the derivative, as you did for the other greeks.
If those adjustments are done, i.e. replacing the theta calculation with
theta = -(black_scholes_price(S, K, T + dT, r, sigma) - black_scholes_price(S, K, T - dT, r, sigma)) / (dT * 2)
I get a theta value of -0.7453608249772259, i.e. the expected result. Also, repeating the calculation shows that now theta has a variance of zero.
(Also, the calculation of the greeks if I am not mistaken should also factor in the option type (put or call), something that still seems missing in the provided code (i.e. my suggestion for the theta calculation works for the put case). However, note that I am not sure on this, as it has been a while since I last dabbled with this kind of stuff, and couldn't be bothered to double check. I merely got this sensation from here.)