pythonfinancemontecarloquantitative-finance

My Monte Carlo simulation keeps giving me the incorrect value for Theta. The other estimated option greeks are correct, but not Theta


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.


Solution

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