Using Monte Carlo Simulation to Estimate Option Probabilities

Introduction

When selling options, one of the most crucial pieces of information is the probability of the option expiring out-of-the-money (OTM). While we commonly use the delta of an option as a quick approximation for this probability, Monte Carlo simulation offers a more detailed look at potential outcomes.

What is Monte Carlo Simulation?

Monte Carlo simulation is a computational technique that uses repeated random sampling to produce numerical results. In options trading, it simulates thousands of possible price paths that a stock might take until expiration, giving us a distribution of potential outcomes.

The Math Behind the Simulation

The simulation uses geometric Brownian motion (GBM), which is a mathematical model commonly used to describe stock price movements. The key components are:

1. Drift: The general direction of price movement (based on risk-free rate)
2. Volatility: The measure of price variation
3. Random Walk: Daily random price changes following a normal distribution

For each simulated path, we calculate daily prices using:

next_price = current_price * exp((r - 0.5 * σ²) * dt + σ * √dt * Z)

Where:

• r = risk-free rate
• σ = volatility
• dt = time step (1/252 for daily trading)
• Z = random number from standard normal distribution

Python Implementation

Let's dive deep into each function that makes up our Monte Carlo simulation for options probability analysis. We'll explore how each function works and why each component is essential for accurate probability estimation.

Time to Expiry

The calculate_time_to_expiry() function converts the time between today and option expiration into years. This conversion is crucial because option pricing models require time expressed in years rather than days.

def calculate_time_to_expiry():
    expiry = datetime(2025, 1, 17)
    today = datetime(2024, 12, 5)
    return (expiry - today).days / 365

Time plays a vital role in option pricing through two primary mechanisms. First, it affects the rate of time decay (theta) of the option's value. Second, the impact of volatility on the option's price scales with the square root of time. This time component is essential for both our Monte Carlo simulation and the Black-Scholes model comparison.

Daily Price Path Generation

The heart of our simulation lies in the generate_daily_paths() function. This function creates multiple possible price paths that the underlying asset might take from today until expiration.

def generate_daily_paths(S0, r, sigma, T, num_simulations):
    days_to_expiry = int(T * 252)  # Trading days
    dt = T/252  # Time step size
    
    times = np.linspace(0, T, days_to_expiry + 1)
    paths = np.zeros((num_simulations, len(times)))
    paths[:, 0] = S0  # Set starting price
    
    Z = np.random.standard_normal((num_simulations, len(times)-1))
    drift = (r - 0.5 * sigma**2) * dt
    diffusion = sigma * np.sqrt(dt)
    
    for t in range(1, len(times)):
        paths[:, t] = paths[:, t-1] * np.exp(drift + diffusion * Z[:, t-1])
    
    return paths, times

The function implements geometric Brownian motion, the standard model for stock price movements in financial mathematics. It takes several key parameters: the current stock price (S0), risk-free rate (r), volatility (sigma), time to expiry (T), and the number of simulations to run.

We use 252 trading days per year rather than 365 calendar days because financial markets are only open on business days. This adjustment provides more accurate modeling of how prices actually evolve in the market.

The price path calculation combines two components: drift and diffusion. The drift term, (r - 0.5 * sigma**2) * dt, represents the expected movement in the stock price based on the risk-free rate, adjusted for volatility. The diffusion term, sigma * sqrt(dt), models the random price movements based on the stock's volatility.

Option Expiry Simulation

The simulate_option_expiry() function analyzes the generated price paths to determine the probabilities of different outcomes at expiration.

def simulate_option_expiry(S0, K, r, sigma, T, num_simulations, option_type='put', confidence_level=0.95):
    paths, times = generate_daily_paths(S0, r, sigma, T, num_simulations)
    final_prices = paths[:, -1]  # Get all final prices
    
    if option_type.lower() == 'call':
        itm_count = np.sum(final_prices > K)
    else:  # put
        itm_count = np.sum(final_prices < K)
    
    probability_itm = itm_count / num_simulations
    probability_otm = 1 - probability_itm
    
    z_score = norm.ppf((1 + confidence_level) / 2)
    margin_of_error = z_score * np.sqrt((probability_itm * (1 - probability_itm)) / num_simulations)
    
    confidence_interval = (
        max(0.0, probability_itm - margin_of_error),
        min(1.0, probability_itm + margin_of_error)
    )
    
    return probability_itm, probability_otm, confidence_interval, final_prices, paths, times

This function examines the final prices from each simulated path and compares them to the option's strike price. For a call option, it counts how many paths end above the strike price; for a put option, it counts paths ending below the strike price. These counts are converted to probabilities by dividing by the total number of simulations.

The function also calculates confidence intervals using the normal approximation method. This gives us a range within which we can be confident the true probability lies, accounting for the random nature of our simulation.

Black-Scholes Delta Calculation

Finally, the calculate_black_scholes_delta() function provides a theoretical benchmark for our Monte Carlo results.

def calculate_black_scholes_delta(S0, K, r, sigma, T, option_type='put'):
    d1 = (np.log(S0/K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    
    if option_type.lower() == 'call':
        delta = norm.cdf(d1)
    else:  # put
        delta = norm.cdf(d1) - 1
    
    return delta

The Black-Scholes delta represents both the option's price sensitivity to changes in the underlying asset's price and approximately the probability of the option expiring in-the-money. This provides a quick sanity check for our Monte Carlo probability estimates. The delta calculation uses the same inputs as our simulation but arrives at the probability through analytical means rather than simulation.

Analyzing GOOG Put Option Probabilities

Let's apply our Monte Carlo simulation to analyze a real-world example: a Google (GOOG) put option with a $185 strike price expiring in January 2025.

Monte Carlo simulation showing potential GOOG price paths until January 2025 expiration, with strike price (red) and current price (green).

We can see from the price paths chart above that while most simulated paths stay above our strike price, the model also captures some extreme scenarios where the price could drop significantly. The red dashed line at $185 represents our put strike price, while the green dashed line shows the current price at $199.

Price distribution of GOOG's simulated prices at January 2025 expiration.

The histogram above provides a different perspective on our simulation results. The bulk of the distribution lies well above our strike price, with the peak centered around the current price level. This aligns with our Monte Carlo probability estimate of 98.71% for the put expiring worthless, though we should remember that real-world market movements can deviate significantly from these theoretical projections.

Results Analysis

Using 100,000 simulations for statistical significance, our analysis reveals some interesting insights about this put option trade.

1. Monte Carlo Probability Estimate:

   • 98.71% probability of the put expiring worthless (above $185)
   • Only 1.29% probability of the put finishing in-the-money
   • 95% Confidence Interval: 1.22% to 1.36% for loss probability

2. Black-Scholes Delta Comparison:

   • Delta: -0.1976
   • Delta-implied success probability: 80.24%

Interesting Divergence

There's a notable divergence between our Monte Carlo results and the Black-Scholes delta implication. The Monte Carlo simulation suggests a much higher probability of success (98.71%) compared to what the delta implies (80.24%). This difference warrants further investigation.

Potential Explanations:

1. The difference between spot price ($199) and strike price ($185) is significant relative to the short time to expiration (0.118 years)
2. The implied volatility of 28.33% might not fully capture the market's view of potential price movements
3. Monte Carlo simulation assumes log-normal returns, which might not perfectly match market dynamics

Limitations of Monte Carlo Simulation

The model assumes geometric Brownian motion with constant volatility, but real markets move in jumps, especially during earnings and major market events. Returns aren't truly log-normal - extreme events happen more frequently than the model suggests. This leads to underestimation of tail risks that are particularly important in options trading.

Moreover, Monte Carlo simulations model securities in isolation, ignoring the conditions of the broader market. The simulation's static assumptions miss these shifting market dynamics, leading to potential risk underestimation when market conditions change.

Conclusion

In all, Monte Carlo simulations are a comprehensive tool for analyzing option expiration probabilities. The simulation approach allows us to visualize potential price paths and understand the distribution of possible outcomes, while the Black-Scholes delta provides a theoretical benchmark.

While no probability estimation method is perfect, Monte Carlo simulation provides valuable insights beyond simple delta approximation. It's particularly useful for visualizing potential outcomes and understanding the range of possibilities when selling options.