This simulation explores delta hedging costs under the assumption of constant volatility, drawing from a financial engineering report on OTC options hedging strategies.

The total cost of hedging is comprised of three main components:

1. Dynamic Delta Hedging P&L: This arises from managing futures positions. When delta increases, futures are bought; when delta decreases, futures are sold. Losses from buying high and selling low, plus any terminal payment (K * sign(S-K)), approximate the option premium.
2. Transaction Costs: Frequent adjustments to the futures position, driven by dynamic delta changes with underlying price movements, incur trading fees.
3. Interest Expense on Margin: Financing the margin required for futures positions generates interest costs.

Simulation Parameters:

• Underlying Price Process: Geometric Brownian Motion
• Monte Carlo Simulations: 1000 paths
• Initial Stock Price (S0): 1
• Maturity (T): 0.25 years
• Strike Price (K): 1
• Constant Volatility (Sigma): 0.2
• Risk-Free Rate (r): 0.05
• Transaction Fee Rate: 0.0006

The simulation generated 1000 distinct price paths for the underlying asset. For each path, the following were calculated:

• Hedged Portfolio P&L
• Transaction Costs
• Capital Usage Costs (Interest)
• Un-hedged P&L from Gamma and Theta exposures.

Code Implementation:

First, the Black-Scholes-Merton (BSM) model is implemented to calculate option Greeks, specifically Gamma and Theta.

from math import log, sqrt, exp
from scipy import stats
import numpy as np

class BSMOptionValuation:
    def __init__(self, S0, K, T, r, sigma, div=0.0):
        self.S0 = float(S0)
        self.K = float(K)
        self.T = float(T)
        self.r = float(r)
        self.sigma = float(sigma)
        self.div = float(div)
        
        if T <= 0 or sigma <= 0:
            self.d1 = float('inf')
            self.d2 = float('inf')
        else:
            self.d1 = ((log(self.S0 / self.K) + (self.r - self.div + 0.5 * self.sigma ** 2) * self.T) / (self.sigma * sqrt(self.T)))
            self.d2 = self.d1 - self.sigma * sqrt(self.T)
    
    def delta(self):
        if self.T <= 0 or self.sigma <= 0:
            return 0.0, 0.0
        delta_call = stats.norm.cdf(self.d1) * exp(-self.div * self.T)
        delta_put = delta_call - 1.0
        return delta_call, delta_put

    def gamma(self):
        if self.T <= 0 or self.sigma <= 0:
            return 0.0
        gamma_val = exp(-self.div * self.T) * stats.norm.pdf(self.d1) / (self.S0 * self.sigma * sqrt(self.T))
        return gamma_val

    def theta(self):
        if self.T <= 0 or self.sigma <= 0:
            return 0.0, 0.0
        
        theta_call_part1 = self.div * self.S0 * exp(-self.div * self.T) * stats.norm.cdf(self.d1)
        theta_call_part2 = self.r * self.K * stats.norm.cdf(self.d2) * exp(-self.r * self.T)
        theta_call_part3 = (self.S0 * exp(-self.div * self.T) * stats.norm.pdf(self.d1) * self.sigma) / (2 * sqrt(self.T))
        theta_call = theta_call_part1 - theta_call_part2 - theta_call_part3
        
        theta_put = theta_call + self.r * self.K * exp(-self.r * self.T) - self.div * self.S0 * exp(-self.div * self.T)
        return theta_call, theta_put

Next, the Monte Carlo simmulation is preformed to generate asset price paths.

import math
import datetime
import numpy as np
import pandas as pd
from scipy.stats import norm
import matplotlib.pyplot as plt

# Assume BSMOptionValuation class is defined above

if __name__ == '__main__':
    
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # Simulation parameters
    S0 = 1.0                                     # Initial asset price
    K = 1.0                                      # Strike price
    T = 0.25                                     # Option maturity (years)
    r = 0.05                                     # Risk-free rate
    sigma = 0.2                                  # Volatility
    fee_rate = 6.0 / 10000.0                     # Transaction fee rate
    num_simulations = 1000                       # Number of simulation paths
    simulation_steps = 63                        # Steps per path (e.g., daily for 3 months)

    dt = T / simulation_steps
    random_shocks = np.random.standard_normal((simulation_steps, num_simulations))

    # Generate asset price paths using Geometric Brownian Motion
    log_returns = (r - 0.5 * sigma**2) * dt + sigma * sqrt(dt) * random_shocks
    price_paths_log = np.cumsum(log_returns, axis=0)
    montecarlo_path = S0 * np.exp(price_paths_log)
    
    # Prepend initial price and format for analysis
    initial_row = np.ones((1, num_simulations))
    stock_price_paths = np.vstack((initial_row, montecarlo_path))
    
    df_prices = pd.DataFrame(stock_price_paths)
    df_prices.plot(legend=False)
    plt.title('Simulated Asset Price Paths')
    plt.xlabel('Time Step')
    plt.ylabel('Asset Price')
    plt.show()

    # Hedging simulation and cost calculation
    times_to_expiry = np.array([(simulation_steps - i) * dt for i in range(simulation_steps + 1)])
    times_to_expiry[simulation_steps] = 0 # Ensure last step is expiry

    # Initialize accounts for various costs
    hedge_portfolio_pnl = np.zeros((simulation_steps + 1, num_simulations))
    transaction_costs = np.zeros((simulation_steps, num_simulations))
    interest_costs = np.zeros((simulation_steps + 1, num_simulations))
    gamma_pnl = np.zeros((simulation_steps, num_simulations))
    theta_pnl = np.zeros((simulation_steps + 1, num_simulations))
    delta_greeks = np.zeros((simulation_steps + 1, num_simulations))
    gamma_greeks = np.zeros((simulation_steps + 1, num_simulations))
    theta_greeks = np.zeros((simulation_steps + 1, num_simulations))

    # Calculate initial Greeks
    initial_option = BSMOptionValuation(S0, K, T, r, sigma)
    initial_delta, _ = initial_option.delta()
    initial_gamma = initial_option.gamma()
    initial_theta, _ = initial_option.theta()

    delta_greeks[0, :] = initial_delta
    gamma_greeks[0, :] = initial_gamma
    theta_greeks[0, :] = initial_theta
    hedge_portfolio_pnl[0, :] = 0 # Start with zero P&L

    for sim_idx in range(num_simulations):
        current_prices = stock_price_paths[:, sim_idx]
        current_times_to_expiry = times_to_expiry

        for step in range(1, simulation_steps + 1):
            prev_price = current_prices[step - 1]
            current_price = current_prices[step]
            prev_time_to_expiry = current_times_to_expiry[step - 1]
            current_time_to_expiry = current_times_to_expiry[step]

            # Update Greeks at curent step
            current_option_state = BSMOptionValuation(current_price, K, current_time_to_expiry, r, sigma)
            current_delta, _ = current_option_state.delta()
            current_gamma = current_option_state.gamma()
            current_theta, _ = current_option_state.theta()

            delta_greeks[step, sim_idx] = current_delta
            gamma_greeks[step, sim_idx] = current_gamma
            theta_greeks[step, sim_idx] = current_theta

            # Calculate P&L components
            # Hedge P&L (futures position adjustment)
            delta_change = delta_greeks[step, sim_idx] - delta_greeks[step - 1, sim_idx]
            hedge_portfolio_pnl[step, sim_idx] = hedge_portfolio_pnl[step - 1, sim_idx] - delta_greeks[step - 1, sim_idx] * (current_price - prev_price)
            
            # Transaction Costs
            transaction_costs[step - 1, sim_idx] = fee_rate * abs(delta_change * current_price)

            # Interest Costs (on margin, approximated by delta * price * r * dt)
            # Interest is accrued based on the position at the start of the period
            interest_costs[step, sim_idx] = interest_costs[step - 1, sim_idx] + delta_greeks[step, sim_idx] * current_price * r * dt

            # Gamma P&L approximation
            gamma_pnl[step - 1, sim_idx] = 0.5 * gamma_greeks[step, sim_idx] * (current_price - prev_price)**2
            
            # Theta P&L approximation
            theta_pnl[step, sim_idx] = theta_pnl[step - 1, sim_idx] + theta_greeks[step - 1, sim_idx] * dt

        # Final settlement at expiry
        if current_prices[simulation_steps] > K:
            hedge_portfolio_pnl[simulation_steps, sim_idx] = hedge_portfolio_pnl[simulation_steps - 1, sim_idx] - delta_greeks[simulation_steps - 1, sim_idx] * (current_prices[simulation_steps] - current_prices[simulation_steps - 1]) + (current_prices[simulation_steps] - K)
        else:
            hedge_portfolio_pnl[simulation_steps, sim_idx] = hedge_portfolio_pnl[simulation_steps - 1, sim_idx] - delta_greeks[simulation_steps - 1, sim_idx] * (current_prices[simulation_steps] - current_prices[simulation_steps - 1])
        
        # Ensure theta is accounted for until the very end
        theta_pnl[simulation_steps, sim_idx] = theta_pnl[simulation_steps - 1, sim_idx] + theta_greeks[simulation_steps - 1, sim_idx] * dt

    # Calculate average costs across all simulations
    avg_hedge_pnl = np.mean(hedge_portfolio_pnl[-1, :])
    avg_transaction_costs = np.mean(np.sum(transaction_costs, axis=0))
    avg_interest_costs = np.mean(interest_costs[-1, :])
    avg_gamma_pnl = np.mean(np.sum(gamma_pnl, axis=0))
    avg_theta_pnl = np.mean(theta_pnl[-1, :])

    print(f