Numerical Methods of Pricing Options

A.C. Jokela

2024-05-26

Stock options are financial derivatives that give the holder the right, but not the obligation, to buy (call options) or sell (put options) an underlying stock at a predetermined price (strike price) on or before a specified date (expiration date). Options are valuable tools for investors and traders to hedge risk, generate income, or speculate on price movements. I may write a post how to use options to generate regular income with only a modicum of risk.

The value of an option is influenced by several factors, including the current stock price, strike price, time to expiration, risk-free interest rate, and the stock's volatility. These factors are quantified using "greeks," which measure the sensitivity of the option's price to changes in these variables. The primary greeks are:

Delta (Δ): The rate of change of the option price with respect to the change in the underlying stock price.
Gamma (Γ): The rate of change of delta with respect to the change in the underlying stock price.
Theta (Θ): The rate of change of the option price with respect to the passage of time.
Vega (ν): The rate of change of the option price with respect to the change in the implied volatility of the underlying stock.
Rho (ρ): The rate of change of the option price with respect to the change in the risk-free interest rate.

To price options and calculate greeks, closed-form models like the Black-Scholes and Bjerksund-Stensland are commonly used. The Black-Scholes model is widely used for European-style options, which can only be exercised at expiration. The model assumes that the underlying stock price follows a geometric Brownian motion with constant volatility and risk-free interest rate. The Black-Scholes formula for a call option is:

Formula
$C = S N (d_{1}) - K e^{- r T} N (d_{2})$
where: $d_{1} = \frac{\ln (\frac{S}{K}) + (r + \frac{σ^{2}}{2}) T}{σ \sqrt{T}}$
$d_{2} = d_{1} - σ \sqrt{T}$

The Bjerksund-Stensland model is an extension of the Black-Scholes model that allows for early exercise, making it suitable for American-style options. The closed-form solution is more complex and involves additional calculations.

The closed-form formulas for the Black-Scholes greeks are:

Greek	Formula
Delta (Δ)	$Δ = N (d_{1})$
Gamma (Γ)	$Γ = \frac{N' (d_{1})}{S σ \sqrt{T}}$
Theta (Θ)	$Θ = - \frac{S N' (d_{1}) σ}{2 \sqrt{T}} - r K e^{- r T} N (d_{2})$
Vega (ν)	$ν = S \sqrt{T} N' (d_{1})$
Rho (ρ)	$ρ = K T e^{- r T} N (d_{2})$

Understanding options and their greeks is crucial for effective risk management and trading strategies. Closed-form models provide a quick and efficient way to price options and calculate greeks, although they rely on certain assumptions that may not always hold in real-world market conditions.

I go over Black-Scholes and Bjerksund-Stensland in more detail in my post on pricing options. In this post, I will focus on numerical methods used to price options, including binomial trees and the Leisen-Reimer method.

Numerical methods and closed-form methods are two different approaches to solving mathematical problems, including the pricing of options and calculation of greeks. Here's how they differ:

Closed-form methods:

Closed-form methods provide an exact solution to a problem using a mathematical formula or equation. These methods give a precise answer that can be calculated directly from the inputs without the need for iterative calculations or approximations. The Black-Scholes model is an example of a closed-form solution for pricing European-style options.

Advantages:

Exact solution
Fast computation
Easy to implement and interpret

Disadvantages:

Rely on assumptions that may not hold in real-world conditions
Limited flexibility to accommodate complex scenarios or non-standard options

Numerical methods:

Numerical methods provide approximate solutions to problems that cannot be solved analytically or have complex closed-form solutions. These methods involve iterative calculations and approximations to arrive at a solution within a specified tolerance level. Examples of numerical methods for option pricing include the Binomial Option Pricing Model (BOPM), Monte Carlo simulations, and Finite Difference Methods%20difference%20equations.) (FDM).

Advantages:

Can handle complex scenarios and non-standard options
Flexible and adaptable to various assumptions and market conditions
Provide insights into option price behavior over time and across various input parameters

Disadvantages:

Approximate solution
Computationally intensive and time-consuming
Requires more programming expertise to implement and interpret

In practice, closed-form methods are preferred when available due to their simplicity and speed. However, numerical methods are essential for dealing with more complex options and market conditions that violate the assumptions of closed-form models.

Numerical methods can also be used to validate closed-form solutions and to provide additional insights into option price behavior. In some cases, a combination of closed-form and numerical methods may be used to balance the trade-offs between accuracy, speed, and flexibility.

We will be focusing on the Binomial Option Pricing Model (BOPM) in this post, which is a popular numerical method for pricing options. The BOPM is based on the concept of a binomial tree, where the option price is calculated at each node of the tree by considering the possible price movements of the underlying stock.

The binomial tree consists of nodes representing the possible stock prices at each time step and branches representing the up and down movements of the stock price. The option price at each node is calculated using the risk-neutral pricing formula, which considers the expected value of the option price at the next time step.

The BOPM can be used to price both European and American-style options and can handle complex scenarios such as dividend payments, early exercise, and multiple time steps. The flexibility and adaptability of the BOPM make it a valuable tool for pricing options in real-world market conditions.

The Leisen-Reimer method is an enhancement of the BOPM that improves the accuracy of option pricing by adjusting the probabilities of up and down movements in the binomial tree. The Leisen-Reimer method uses a modified version of the risk-neutral pricing formula that accounts for skewness and kurtosis in the stock price distribution.

The Leisen-Reimer binomial tree is an improvement over the simple binomial tree method for pricing options. It uses a more accurate approximation for the up and down factors in the tree, leading to faster convergence and more precise results. Here's a simple pseudo-code comparison between the two methods:

Simple Binomial Tree:

function SimpleBinomialTree(S, K, T, r, sigma, N):
    dt = T / N
    u = exp(sigma * sqrt(dt))
    d = 1 / u
    p = (exp(r * dt) - d) / (u - d)

    Initialize asset_prices[N+1] to 0
    asset_prices[0] = S * d^N

    for i from 1 to N:
        asset_prices[i] = u * asset_prices[i-1]

    Initialize option_values[N+1] to 0

    for i from 0 to N:
        option_values[i] = max(0, asset_prices[i] - K)

    for i from N-1 down to 0:
        for j from 0 to i:
            option_values[j] = exp(-r * dt) * (p * option_values[j+1] + (1-p) * option_values[j])

    return option_values[0]

Leisen-Reimer Binomial Tree:

function LeisenReimerBinomialTree(S, K, T, r, sigma, N):
    dt = T / N
    d1 = (log(S / K) + (r + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
    d2 = d1 - sigma * sqrt(T)

    p = 0.5 + copysign(1, d2) * 0.5 * sqrt(1 - exp(-((d2 / (N + 1/3 + 0.1/(N+1)))^2) * (N + 1/6)))
    u = exp(sigma * sqrt(dt / p))
    d = exp(-sigma * sqrt(dt / (1 - p)))

    Initialize asset_prices[N+1] to 0
    asset_prices[0] = S * d^N

    for i from 1 to N:
        asset_prices[i] = u * asset_prices[i-1]

    Initialize option_values[N+1] to 0

    for i from 0 to N:
        option_values[i] = max(0, asset_prices[i] - K)

    for i from N-1 down to 0:
        for j from 0 to i:
            option_values[j] = exp(-r * dt) * (p * option_values[j+1] + (1-p) * option_values[j])

    return option_values[0]

The main differences between the two methods are:

The Leisen-Reimer binomial tree method incorporates the d1 and d2 values from the Black-Scholes model to determine the up and down factors (u and d) in the binomial tree. In the Black-Scholes model, d1 and d2 are used to calculate the probability of the option expiring in-the-money under the assumption of a log-normal distribution of stock prices. By using these values in the binomial tree, the Leisen-Reimer method captures the essence of the Black-Scholes model's assumptions while retaining the flexibility and intuitive appeal of the binomial tree framework. This approach leads to a more accurate approximation of the continuous-time process assumed by the Black-Scholes model, resulting in faster convergence and more precise option pricing compared to the simple binomial tree method.
The probability p is calculated differently in the Leisen-Reimer method, using a more complex formula that includes a correction term (1/3 + 0.1/(N+1)) to improve convergence.

These modifications lead to faster convergence and more accurate results compared to the simple binomial tree method, especially for a smaller number of time steps (N). The Leisen-Reimer method is particularly useful for options with high volatility or near-the-money strike prices, where the standard BOPM may produce inaccurate results.

When comparing the computational complexity of the binomial tree methods (both simple and Leisen-Reimer) to the Black-Scholes model, there is a significant difference in their time and space complexities.

Time Complexity: - Black-Scholes: O(1) The Black-Scholes formula involves a fixed number of arithmetic operations and standard normal distribution function evaluations, resulting in a constant time complexity. Regardless of the input parameters, the computation time remains the same.

Binomial Trees (Simple and Leisen-Reimer): O(N^2) Both the simple binomial tree and the Leisen-Reimer binomial tree have a time complexity of O(N^2), where N is the number of time steps in the tree. The computation time grows quadratically with the number of time steps, making them significantly slower than the Black-Scholes model, especially for large values of N.

Space Complexity: - Black-Scholes: O(1) The Black-Scholes formula only requires a constant amount of memory to store the input parameters and the calculated option price, resulting in a space complexity of O(1).

Binomial Trees (Simple and Leisen-Reimer): O(N) Both the simple binomial tree and the Leisen-Reimer binomial tree have a space complexity of O(N), as they require storing the asset price tree and the option value array, each of size N+1. The memory requirements grow linearly with the number of time steps.

The Black-Scholes model has a significantly lower computational complexity compared to the binomial tree methods. The Black-Scholes model has a constant time complexity of O(1) and a constant space complexity of O(1), while the binomial tree methods have a quadratic time complexity of O(N^2) and a linear space complexity of O(N).

However, it's important to note that the binomial tree methods offer more flexibility in handling complex options and relaxing assumptions, which may justify their use in certain situations despite their higher computational complexity. The choice between the Black-Scholes model and binomial tree methods depends on the specific requirements of the problem and the trade-offs between accuracy, flexibility, and computational efficiency.

// binomial_lr_with_greeks.rs

use crate::binomial_lr_option::BinomialLROption;

/// Represents a binomial LR (Leisen-Reimer) option with Greeks calculation.
///
/// This struct extends the `BinomialLROption` to include the calculation of option Greeks,
/// such as delta, gamma, theta, vega, and rho.
pub struct BinomialLRWithGreeks {
    /// The underlying binomial LR option.
    pub lr_option: BinomialLROption,
}

impl BinomialLRWithGreeks {
    /// Creates a new `BinomialLRWithGreeks` instance with the given `BinomialLROption`.
    ///
    /// # Arguments
    ///
    /// * `lr_option` - The binomial LR option to be used for Greeks calculation.
    pub fn new(lr_option: BinomialLROption) -> Self {
        BinomialLRWithGreeks { lr_option }
    }

    /// Generates a new stock price tree based on the binomial LR option parameters.
    ///
    /// This method calculates the stock prices at each node of the binomial tree using
    /// the up and down factors from the binomial LR option.
    fn new_stock_price_tree(&mut self) {
        let u_over_d = self.lr_option.tree.u / self.lr_option.tree.d;
        let d_over_u = self.lr_option.tree.d / self.lr_option.tree.u;

        self.lr_option.tree.option.sts = vec![vec![
            self.lr_option.tree.option.s0 * u_over_d,
            self.lr_option.tree.option.s0,
            self.lr_option.tree.option.s0 * d_over_u,
        ]];

        for _ in 0..self.lr_option.tree.option.n {
            let prev_branches = &self.lr_option.tree.option.sts[self.lr_option.tree.option.sts.len() - 1];
            let mut st = prev_branches
                .iter()
                .map(|&x| x * self.lr_option.tree.u)
                .collect::<Vec<_>>();
            st.push(prev_branches[prev_branches.len() - 1] * self.lr_option.tree.d);
            self.lr_option.tree.option.sts.push(st);
        }
    }

    /// Calculates the option price and Greeks (delta, gamma, theta, vega, rho).
    ///
    /// This method first sets up the binomial LR option parameters and generates the stock price tree.
    /// It then calculates the option payoffs using the `begin_tree_traversal` method from the binomial LR option.
    /// Finally, it computes the option price and various Greeks based on the calculated payoffs and stock prices.
    ///
    /// # Returns
    ///
    /// A tuple containing the following values:
    /// - `option_value`: The calculated option price.
    /// - `delta`: The option's delta (rate of change of option price with respect to the underlying asset price).
    /// - `gamma`: The option's gamma (rate of change of delta with respect to the underlying asset price).
    /// - `theta`: The option's theta (rate of change of option price with respect to time).
    /// - `vega`: The option's vega (sensitivity of option price to changes in volatility).
    /// - `rho`: The option's rho (sensitivity of option price to changes in the risk-free interest rate).
    pub fn price(&mut self) -> (f64, f64, f64, f64, f64, f64) {
        self.lr_option.setup_parameters();
        self.new_stock_price_tree();

        let payoffs = self.lr_option.tree.begin_tree_traversal();
        let option_value = payoffs[payoffs.len() / 2];
        let payoff_up = payoffs[0];
        let payoff_down = payoffs[payoffs.len() - 1];

        let s_up = self.lr_option.tree.option.sts[0][0];
        let s_down = self.lr_option.tree.option.sts[0][2];

        let ds_up = s_up - self.lr_option.tree.option.s0;
        let ds_down = self.lr_option.tree.option.s0 - s_down;
        let ds = s_up - s_down;
        let dv = payoff_up - payoff_down;

        // Calculate delta as the change in option value divided by the change in stock price
        let delta = dv / ds;

        // Calculate gamma as the change in delta divided by the change in stock price
        let gamma = ((payoff_up - option_value) / ds_up - (option_value - payoff_down) / ds_down)
            / ((self.lr_option.tree.option.s0 + s_up) / 2.0 - (self.lr_option.tree.option.s0 + s_down) / 2.0);

        let dt = 0.0001; // Small perturbation in time
        let original_t = self.lr_option.tree.option.t;
        self.lr_option.tree.option.t -= dt;
        self.lr_option.setup_parameters();
        let payoffs_theta = self.lr_option.tree.begin_tree_traversal();
        let option_value_theta = payoffs_theta[payoffs_theta.len() / 2];

        // Calculate theta as the negative of the change in option value divided by the change in time
        let theta = -(option_value_theta - option_value) / dt;
        self.lr_option.tree.option.t = original_t;

        let dv = 0.01;
        self.lr_option.tree.option.sigma += dv;
        self.lr_option.setup_parameters();
        let payoffs_vega = self.lr_option.tree.begin_tree_traversal();
        let option_value_vega = payoffs_vega[payoffs_vega.len() / 2];

        // Calculate vega as the change in option value divided by the change in volatility
        let vega = (option_value_vega - option_value) / dv;
        self.lr_option.tree.option.sigma -= dv;

        let dr = 0.01;
        self.lr_option.tree.option.r += dr;
        self.lr_option.setup_parameters();
        let payoffs_rho = self.lr_option.tree.begin_tree_traversal();
        let option_value_rho = payoffs_rho[payoffs_rho.len() / 2];

        // Calculate rho as the change in option value divided by the change in interest rate
        let rho = (option_value_rho - option_value) / dr;
        self.lr_option.tree.option.r -= dr;

        (option_value, delta, gamma, theta, vega, rho)
    }
}

Numerical methods like the Binomial Option Pricing Model and the Leisen-Reimer method are essential tools for pricing options and calculating greeks in real-world market conditions. These methods provide flexibility, accuracy, and insights into option price behavior that closed-form models may not capture. By understanding the strengths and limitations of numerical methods, traders and investors can make informed decisions and manage risk effectively in the options market.

The complete source code for this project can be found on Github. Binomial Option Pricing Model with Greeks This code is written in Rust, but is based off of the Python code in the book "Mastering Python for Finance" Second Edition by James Ma Weiming. The book is a great resource for learning how to price options using numerical methods among many finance topic problems all solvable and explorable in Python. The Rust code is a translation of the Python code with the addition of calculating all of the Greeks.

The Rust code compiles to a library that can then be imported in Python. Making this a Python module allows for easier integration with other Python libraries and tools. The Python code can be used to visualize the results and provide a user interface for the Rust code.