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, riskfree 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 riskfree interest rate.
To price options and calculate greeks, closedform models like the BlackScholes and BjerksundStensland are commonly used. The BlackScholes model is widely used for Europeanstyle 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 riskfree interest rate. The BlackScholes formula for a call option is:
Formula 

$C=SN\left({d}_{1}\right)K{e}^{rT}N\left({d}_{2}\right)$ 
where: ${d}_{1}=\frac{\mathrm{ln}\left(\frac{S}{K}\right)+(r+\frac{{\sigma}^{2}}{2})T}{\sigma \sqrt{T}}$ 
${d}_{2}={d}_{1}\sigma \sqrt{T}$ 
The BjerksundStensland model is an extension of the BlackScholes model that allows for early exercise, making it suitable for Americanstyle options. The closedform solution is more complex and involves additional calculations.
The closedform formulas for the BlackScholes greeks are:
Greek  Formula 

Delta (Δ)  $\Delta =N\left({d}_{1}\right)$ 
Gamma (Γ)  $\Gamma =\frac{N\text{'}\left({d}_{1}\right)}{S\sigma \sqrt{T}}$ 
Theta (Θ)  $\Theta =\frac{SN\text{'}\left({d}_{1}\right)\sigma}{2\sqrt{T}}rK{e}^{rT}N\left({d}_{2}\right)$ 
Vega (ν)  $\nu =S\sqrt{T}N\text{'}\left({d}_{1}\right)$ 
Rho (ρ)  $\rho =KT{e}^{rT}N\left({d}_{2}\right)$ 
Understanding options and their greeks is crucial for effective risk management and trading strategies. Closedform 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 realworld market conditions.
I go over BlackScholes and BjerksundStensland 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 LeisenReimer method.
Numerical methods and closedform methods are two different approaches to solving mathematical problems, including the pricing of options and calculation of greeks. Here's how they differ:
Closedform methods:
Closedform 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 BlackScholes model is an example of a closedform solution for pricing Europeanstyle options.
Advantages:
 Exact solution
 Fast computation
 Easy to implement and interpret
Disadvantages:
 Rely on assumptions that may not hold in realworld conditions
 Limited flexibility to accommodate complex scenarios or nonstandard options
Numerical methods:
Numerical methods provide approximate solutions to problems that cannot be solved analytically or have complex closedform 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 nonstandard 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 timeconsuming
 Requires more programming expertise to implement and interpret
In practice, closedform 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 closedform models.
Numerical methods can also be used to validate closedform solutions and to provide additional insights into option price behavior. In some cases, a combination of closedform and numerical methods may be used to balance the tradeoffs 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 riskneutral 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 Americanstyle 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 realworld market conditions.
The LeisenReimer 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 LeisenReimer method uses a modified version of the riskneutral pricing formula that accounts for skewness and kurtosis in the stock price distribution.
The LeisenReimer 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 pseudocode 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[i1] 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 N1 down to 0: for j from 0 to i: option_values[j] = exp(r * dt) * (p * option_values[j+1] + (1p) * option_values[j]) return option_values[0]
LeisenReimer 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[i1] 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 N1 down to 0: for j from 0 to i: option_values[j] = exp(r * dt) * (p * option_values[j+1] + (1p) * option_values[j]) return option_values[0]
The main differences between the two methods are:

The LeisenReimer binomial tree method incorporates the d1 and d2 values from the BlackScholes model to determine the up and down factors (u and d) in the binomial tree. In the BlackScholes model, d1 and d2 are used to calculate the probability of the option expiring inthemoney under the assumption of a lognormal distribution of stock prices. By using these values in the binomial tree, the LeisenReimer method captures the essence of the BlackScholes 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 continuoustime process assumed by the BlackScholes 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 LeisenReimer 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 LeisenReimer method is particularly useful for options with high volatility or nearthemoney strike prices, where the standard BOPM may produce inaccurate results.
When comparing the computational complexity of the binomial tree methods (both simple and LeisenReimer) to the BlackScholes model, there is a significant difference in their time and space complexities.
Time Complexity:  BlackScholes: O(1) The BlackScholes 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 LeisenReimer): O(N^2) Both the simple binomial tree and the LeisenReimer 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 BlackScholes model, especially for large values of N.
Space Complexity:  BlackScholes: O(1) The BlackScholes 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 LeisenReimer): O(N) Both the simple binomial tree and the LeisenReimer 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 BlackScholes model has a significantly lower computational complexity compared to the binomial tree methods. The BlackScholes 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 BlackScholes model and binomial tree methods depends on the specific requirements of the problem and the tradeoffs between accuracy, flexibility, and computational efficiency.
// binomial_lr_with_greeks.rs use crate::binomial_lr_option::BinomialLROption; /// Represents a binomial LR (LeisenReimer) 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 riskfree 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 LeisenReimer method are essential tools for pricing options and calculating greeks in realworld market conditions. These methods provide flexibility, accuracy, and insights into option price behavior that closedform 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.