Monte Carlo Method of Pricing Options

A.C. Jokela

2024-06-19

Pricing American Options with the Monte Carlo Method

Introduction to Monte Carlo Simulation

Monte Carlo simulation is a powerful computational technique that relies on repeated random sampling to solve complex problems. It is widely used in various fields, including finance, physics, and engineering. The core idea behind Monte Carlo simulation is to generate a large number of random scenarios and then analyze the results to gain insights into the problem at hand.

To illustrate the concept, let's consider a simplified example of estimating the value of π (pi) using Monte Carlo simulation. Imagine a square with sides of length 1 and a circle inscribed within the square. The circle's radius is 0.5, and its area is π/4. By randomly generating points within the square and counting the number of points that fall inside the circle, we can estimate the value of π.

Here's a step-by-step breakdown of the process:

Generate a large number of random points (x, y) within the square, where both x and y are between 0 and 1.
For each point, calculate the distance from the point to the center of the square (0.5, 0.5) using the formula:

distance = \sqrt{{(x - 0.5)}^{2} + {(y - 0.5)}^{2}}

If the distance is less than or equal to 0.5 (the radius of the circle), count the point as inside the circle.
After generating all the points, calculate the ratio of points inside the circle to the total number of points.
Multiply the ratio by 4 to estimate the value of π.

As the number of generated points increases, the estimated value of π will converge to its true value. This example demonstrates how Monte Carlo simulation can be used to solve problems by generating random scenarios and analyzing the results.

Pricing American Options with Monte Carlo Simulation

American options are financial contracts that give the holder the right, but not the obligation, to buy (call option) or sell (put option) an underlying asset at a predetermined price (strike price) on or before the expiration date. Unlike European options, which can only be exercised at expiration, American options can be exercised at any time prior to expiration. This early exercise feature makes pricing American options more challenging than pricing European options.

Monte Carlo simulation can be used to price American options by simulating the price paths of the underlying asset and determining the optimal exercise strategy at each time step. Here's a high-level overview of the process:

Generate a large number of price paths for the underlying asset using a suitable stochastic process, such as Geometric Brownian Motion (GBM). Each price path represents a possible future scenario. The price path can be simulated using the following formula:

S_{t} = S_{0} \exp ((r - \frac{σ^{2}}{2}) t + σ \sqrt{t} Z)

where: - $S_{t}$ is the stock price at time t - $S_{0}$ is the initial stock price - $r$ is the risk-free interest rate - $σ$ is the volatility of the stock - $t$ is the time step - $Z$ is a standard normal random variable

At each time step along each price path, compare the intrinsic value of the option (the payoff from immediate exercise) with the continuation value (the expected value of holding the option).
If the intrinsic value is greater than the continuation value, the option should be exercised, and the payoff is recorded.
If the continuation value is greater than the intrinsic value, the option should be held, and the simulation continues to the next time step.
At expiration, the option is exercised if it is in-the-money (i.e., the stock price is above the strike price for a call option or below the strike price for a put option).
The average of the discounted payoffs across all price paths provides an estimate of the option price.

To determine the continuation value at each time step, a common approach is to use the Least Squares Monte Carlo (LSM) method, which involves regressing the discounted future payoffs on a set of basis functions of the current state variables (e.g., stock price and time to expiration).

One of the advantages of using Monte Carlo simulation for pricing American options is its flexibility in handling complex payoff structures and multiple underlying assets. However, it can be computationally intensive, especially for high-dimensional problems or when a large number of price paths is required for accurate results.

In practice, implementing the Monte Carlo method for pricing American options involves several technical considerations, such as choosing an appropriate stochastic process for the underlying asset, selecting suitable basis functions for the LSM algorithm, and applying variance reduction techniques to improve the efficiency of the simulation.

Despite its challenges, Monte Carlo simulation remains a valuable tool for pricing American options and has been widely adopted in the financial industry. Its ability to handle complex scenarios and provide insights into the option's behavior under various market conditions makes it an indispensable technique in the field of quantitative finance.

Implementing the Monte Carlo Method for American Option Pricing

To implement the Monte Carlo method for pricing American options, we need to follow a step-by-step approach. Here's a more detailed breakdown of the process:

Define the parameters of the American option, such as the strike price (K), expiration time (T), risk-free interest rate (r), and volatility (σ).
Generate a large number of price paths for the underlying asset using the Geometric Brownian Motion (GBM) model. Each price path consists of discrete time steps from the current time to the expiration time. The stock price at each time step can be calculated using the following formula:

S_{t + Δ t} = S_{t} \exp ((r - \frac{σ^{2}}{2}) Δ t + σ \sqrt{Δ t} Z_{t})

where: - $S_{t + Δ t}$ is the stock price at time t + Δt - $S_{t}$ is the stock price at time t - $Δ t$ is the size of the time step - $Z_{t}$ is a standard normal random variable at time t

At each time step along each price path, calculate the intrinsic value of the option. For a call option, the intrinsic value is given by:

\max (S_{t} - K, 0)

For a put option, the intrinsic value is given by:

\max (K - S_{t}, 0)

Apply the Least Squares Monte Carlo (LSM) method to estimate the continuation value at each time step. This involves the following substeps:
Select a set of basis functions that depend on the current stock price and time to expiration. Common choices include polynomials, Laguerre polynomials, or exponential functions.
For each in-the-money price path at the current time step, calculate the discounted future payoffs from continuing to hold the option.
Regress the discounted future payoffs on the basis functions to obtain an estimate of the continuation value as a function of the current stock price.
Compare the intrinsic value with the estimated continuation value at each time step. If the intrinsic value is greater, the option should be exercised, and the payoff is recorded. If the continuation value is greater, the option should be held, and the simulation continues to the next time step.
At expiration, the option is exercised if it is in-the-money, and the payoff is recorded.
Discount the payoffs at each exercise time back to the present using the risk-free interest rate. The average of the discounted payoffs across all price paths provides an estimate of the American option price.

The LSM method used in step 4 is a critical component of the Monte Carlo approach for pricing American options. It allows us to estimate the continuation value at each time step without the need for a recursive tree-based method. The choice of basis functions and the number of basis functions used in the regression can impact the accuracy of the estimated continuation value and, consequently, the option price.

Once the Monte Carlo simulation is complete, we can analyze the results to gain insights into the behavior of the American option under various market conditions. This can include examining the distribution of payoffs, the probability of early exercise, and the sensitivity of the option price to changes in the underlying parameters (known as the "Greeks").

It's important to note that the accuracy of the Monte Carlo method for pricing American options depends on several factors, such as the number of price paths simulated, the size of the time steps, the choice of basis functions in the LSM method, and the presence of any variance reduction techniques. Increasing the number of price paths and reducing the size of the time steps can improve the accuracy of the results but also increase the computational burden.

In practice, implementing the Monte Carlo method for American option pricing requires careful consideration of these factors and may involve trade-offs between accuracy and computational efficiency. Nonetheless, the Monte Carlo approach remains a powerful and flexible tool for tackling the complex problem of pricing American options, particularly in cases where analytical solutions are not available or when dealing with multi-dimensional or path-dependent options.

Code

Monte Carlo Options Pricing in Rust

Jupyter Notebook with usage example

use rand::Rng;
use std::f64::consts::E;
use pyo3::prelude::*;

/// Calculates the price of an American call option using the Monte Carlo simulation method.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run.
/// * `num_steps` - The number of steps in each simulation.
///
/// # Returns
///
/// The price of the American call option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::monte_carlo_american_call;
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let price = monte_carlo_american_call(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("American call option price: {:.2}", price);
/// ```
#[pyfunction]
fn monte_carlo_american_call(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let dt = time_to_maturity / num_steps as f64;
    let discount_factor = E.powf(-risk_free_rate * time_to_maturity);

    let mut rng = rand::thread_rng();
    let mut payoffs = Vec::with_capacity(num_simulations);

    for _ in 0..num_simulations {
        let mut price = spot_price;
        let mut max_payoff: f64 = 0.0;

        for _ in 0..num_steps {
            let rand_num = rng.gen_range(-1.0..1.0);
            price *= E.powf(
                (risk_free_rate - 0.5 * volatility.powi(2)) * dt
                    + volatility * rand_num * dt.sqrt(),
            );
            max_payoff = max_payoff.max((price - strike_price).max(0.0));
        }

        payoffs.push(max_payoff);
    }

    let sum_payoffs: f64 = payoffs.iter().sum();
    let avg_payoff = sum_payoffs / num_simulations as f64;

    Ok(avg_payoff * discount_factor)
}

/// Calculates the price of an American put option using the Monte Carlo simulation method.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run.
/// * `num_steps` - The number of steps in each simulation.
///
/// # Returns
///
/// The price of the American put option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::monte_carlo_american_put;
///
/// let spot_price = 100.0;
/// let strike_price = 90.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let price = monte_carlo_american_put(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("American put option price: {:.2}", price);
/// ```
#[pyfunction]
fn monte_carlo_american_put(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let dt = time_to_maturity / num_steps as f64;
    let discount_factor = E.powf(-risk_free_rate * time_to_maturity);

    let mut rng = rand::thread_rng();
    let mut payoffs = Vec::with_capacity(num_simulations);

    for _ in 0..num_simulations {
        let mut price = spot_price;
        let mut max_payoff: f64 = 0.0;

        for _ in 0..num_steps {
            let rand_num = rng.gen_range(-1.0..1.0);
            price *= E.powf(
                (risk_free_rate - 0.5 * volatility.powi(2)) * dt
                    + volatility * rand_num * dt.sqrt(),
            );
            max_payoff = max_payoff.max((strike_price - price).max(0.0));
        }

        payoffs.push(max_payoff);
    }

    let sum_payoffs: f64 = payoffs.iter().sum();
    let avg_payoff = sum_payoffs / num_simulations as f64;

    Ok(avg_payoff * discount_factor)
}

/// Calculates the delta of an American call option using the finite difference method.
///
/// Delta represents the rate of change of the option price with respect to the change in the underlying asset price.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run for the Monte Carlo method.
/// * `num_steps` - The number of steps in each simulation for the Monte Carlo method.
///
/// # Returns
///
/// The delta of the American call option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::{calculate_delta, monte_carlo_american_call};
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let delta = calculate_delta(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("Delta of the American call option: {:.4}", delta);
/// ```
#[pyfunction]
fn calculate_delta(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let h = 0.01;
    let price_up = monte_carlo_american_call(
        spot_price * (1.0 + h),
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    let price_down = monte_carlo_american_call(
        spot_price * (1.0 - h),
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    Ok((price_up - price_down) / (2.0 * h * spot_price))
}

/// Calculates the gamma of an American call option using the finite difference method.
///
/// Gamma represents the rate of change of the option's delta with respect to the change in the underlying asset price.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run for the Monte Carlo method.
/// * `num_steps` - The number of steps in each simulation for the Monte Carlo method.
///
/// # Returns
///
/// The gamma of the American call option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::{calculate_gamma, monte_carlo_american_call};
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let gamma = calculate_gamma(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("Gamma of the American call option: {:.4}", gamma);
/// ```
#[pyfunction]
fn calculate_gamma(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let h = 0.01;
    let price_up = monte_carlo_american_call(
        spot_price * (1.0 + h),
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    let price = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    let price_down = monte_carlo_american_call(
        spot_price * (1.0 - h),
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    Ok((price_up - 2.0 * price + price_down) / (h * spot_price).powi(2))
}

/// Calculates the vega of an American call option using the finite difference method.
///
/// Vega represents the sensitivity of the option price to changes in the volatility of the underlying asset.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run for the Monte Carlo method.
/// * `num_steps` - The number of steps in each simulation for the Monte Carlo method.
///
/// # Returns
///
/// The vega of the American call option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::{calculate_vega, monte_carlo_american_call};
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let vega = calculate_vega(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("Vega of the American call option: {:.4}", vega);
/// ```
#[pyfunction]
fn calculate_vega(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let h = 0.01;
    let price_up = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate,
        volatility * (1.0 + h),
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    let price_down = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate,
        volatility * (1.0 - h),
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    Ok((price_up - price_down) / (2.0 * h * volatility))
}

/// Calculates the theta of an American call option using the finite difference method.
///
/// Theta represents the rate of change of the option price with respect to the passage of time.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run for the Monte Carlo method.
/// * `num_steps` - The number of steps in each simulation for the Monte Carlo method.
///
/// # Returns
///
/// The theta of the American call option (per day).
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::{calculate_theta, monte_carlo_american_call};
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let theta = calculate_theta(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("Theta of the American call option (per day): {:.4}", theta);
/// ```
#[pyfunction]
fn calculate_theta(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let h = 0.01;
    let price_up = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity * (1.0 + h),
        num_simulations,
        num_steps,
    )?;
    let price_down = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate,
        volatility,
        time_to_maturity * (1.0 - h),
        num_simulations,
        num_steps,
    )?;
    Ok(-(price_up - price_down) / (2.0 * h * time_to_maturity) / 365.0)
}

/// Calculates the rho of an American call option using the finite difference method.
///
/// Rho represents the sensitivity of the option price to changes in the risk-free interest rate.
///
/// # Arguments
///
/// * `spot_price` - The current price of the underlying asset.
/// * `strike_price` - The strike price of the option.
/// * `risk_free_rate` - The risk-free interest rate.
/// * `volatility` - The volatility of the underlying asset.
/// * `time_to_maturity` - The time to maturity of the option (in years).
/// * `num_simulations` - The number of simulations to run for the Monte Carlo method.
/// * `num_steps` - The number of steps in each simulation for the Monte Carlo method.
///
/// # Returns
///
/// The rho of the American call option.
///
/// # Example
///
/// ```rust
/// use monte_carlo_options_rs::{calculate_rho, monte_carlo_american_call};
///
/// let spot_price = 100.0;
/// let strike_price = 110.0;
/// let risk_free_rate = 0.05;
/// let volatility = 0.2;
/// let time_to_maturity = 1.0;
/// let num_simulations = 100_000;
/// let num_steps = 252;
///
/// let rho = calculate_rho(
///     spot_price,
///     strike_price,
///     risk_free_rate,
///     volatility,
///     time_to_maturity,
///     num_simulations,
///     num_steps,
/// ).unwrap();
///
/// println!("Rho of the American call option: {:.4}", rho);
/// ```
#[pyfunction]
fn calculate_rho(
    spot_price: f64,
    strike_price: f64,
    risk_free_rate: f64,
    volatility: f64,
    time_to_maturity: f64,
    num_simulations: usize,
    num_steps: usize,
) -> PyResult<f64> {
    let h = 0.01;
    let price_up = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate * (1.0 + h),
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    let price_down = monte_carlo_american_call(
        spot_price,
        strike_price,
        risk_free_rate * (1.0 - h),
        volatility,
        time_to_maturity,
        num_simulations,
        num_steps,
    )?;
    Ok((price_up - price_down) / (2.0 * h))
}

/// Initializes the `libmonte_carlo_options_rs` Python module.
///
/// This function is called when the Python module is imported and adds the Rust functions
/// to the module so they can be accessed from Python.
///
/// # Arguments
///
/// * `_py` - The Python interpreter.
/// * `m` - The Python module to add the functions to.
///
/// # Returns
///
/// `Ok(())` if the module was initialized successfully, or an error if there was a problem.
///
/// # Example
///
/// ```python
/// import libmonte_carlo_options_rs
///
/// # Use the functions from the module
/// price = libmonte_carlo_options_rs.monte_carlo_american_call(...)
/// delta = libmonte_carlo_options_rs.calculate_delta(...)
/// gamma = libmonte_carlo_options_rs.calculate_gamma(...)
/// vega = libmonte_carlo_options_rs.calculate_vega(...)
/// theta = libmonte_carlo_options_rs.calculate_theta(...)
/// rho = libmonte_carlo_options_rs.calculate_rho(...)
/// ```
#[pymodule]
fn libmonte_carlo_options_rs(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
    m.add_function(wrap_pyfunction!(monte_carlo_american_call, m)?)?;
    m.add_function(wrap_pyfunction!(monte_carlo_american_put, m)?)?;
    m.add_function(wrap_pyfunction!(calculate_delta, m)?)?;
    m.add_function(wrap_pyfunction!(calculate_gamma, m)?)?;
    m.add_function(wrap_pyfunction!(calculate_vega, m)?)?;
    m.add_function(wrap_pyfunction!(calculate_theta, m)?)?;
    m.add_function(wrap_pyfunction!(calculate_rho, m)?)?;
    Ok(())
}