Coding towards CFA (9) – From Binomial Model to BSM

*More articles can be found from my blog site - https://dataninjago.com

To kick off this blog post, let’s start with a quick experiment where we compare option prices derived from the binomial model with increasing periods to those calculated using the Black-Scholes Model (BSM).

Here, we run two sets of option pricing calculations. In the first set, we use the multi-period binomial model, calculating the option price with an increasing number of steps, ranging from 1 to 100. In the second set, we calculate the price of the same option using the Black-Scholes Model (BSM). * The code for both sets of runs can be found at the bottom of this post.

The diagram below illustrates the results.

From the results, we can observe that as more periods are added to the binomial model, the calculated option price gradually converges to the price derived from the BSM. In this blog post, I aim to explore the connection between the binomial model and the BSM.

In short, the BSM model can be seen as the continuous-time limit of the discrete-time binomial model. In the binomial model, time is divided into discrete intervals, and at each step, the price of the underlying asset can either move up or down. As the number of time steps increases, the up and down factors are adjusted to reflect smaller and more frequent price movements in the asset. When the number of steps becomes sufficiently large, the model’s discrete price movements become finer, and the binomial model begins to approximate the continuous-time process observed in the real world.

In other words, as the number of steps approaches infinity and the time intervals shrink to infinitesimally small values, the binomial model converges to the BSM, providing a closer approximation to real-world price behaviour. Therefore, that would be no problem to view the BSM as an extreme step beyond the binomial model.

In the continuous-time limit, the discrete binomial model transitions to a geometric Brownian motion (GBM). The continuous-time assumption is the foundation of the BSM, where the option price is derived based on the assumption that the underlying asset follows a GBM with constant volatility and drift. Here is the stochastic differential equation (SDE) describing GBM:

In the continuous form, this can be written as:

• S(t) – the price at current step
• S(t+1) – the price at the next step
• u – drift
• sigma – volatility
• t – time step
• Z – random variable from a standard normal distribution

Here is the Python code implementing the formula:

We can kick off the multiple parallel runs of the code to simulate the potential asset prices paths.

We can then draw the distribution of the final prices of those prices paths.

Full Code – Python

Binomial Model vs BSM

import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm
 
def calculate_up_down_factor(time_per_period, sigma):
 
    u = math.exp(sigma * math.sqrt(time_per_period))
    d = math.exp(-sigma * math.sqrt(time_per_period))
 
    return u, d
 
 
def calculate_up_down_probability(u, d, r, time_per_period):
     
    Pi = (math.exp(r*time_per_period)-d)/(u-d)
 
    return Pi, 1-Pi
 
 
def binomial_tree_m_period(s0, k, T, r, sigma, number_of_periods, option_type="call"):
 
    time_per_period = T/number_of_periods
     
    # Calculate the up and down factors
    u, d = calculate_up_down_factor(time_per_period, sigma)
     
    # Calculate the up and down probabilities
    Pi_up, Pi_down = calculate_up_down_probability(u, d, r, time_per_period)
     
    # calculate underlying prices at expiration
    stock_prices = [0] * (number_of_periods+1)
    for i in range(number_of_periods+1):
        stock_prices[i] = s0*(u**i)*(d**(number_of_periods-i))
 
    # calculate the option payoff at expiration
    payoffs = [0] * (number_of_periods+1)
    for i in range(number_of_periods+1):
        if option_type == "call":
            payoffs[i] = max(stock_prices[i] - k, 0)
        else:
            payoffs[i] = max(k - stock_prices[i], 0)
     
    # Backward induction of option value at the initiation
    discount_factor = 1/math.exp(r*time_per_period)
    for i in range(number_of_periods-1, -1, -1):
        for j in range(i+1):
            payoffs[j] = (Pi_up*payoffs[j+1]+Pi_down*payoffs[j])*discount_factor
     
    option_price = payoffs[0]
 
    return option_price
 
 
def plot_tree_vs_bsm(option_prices_binomial_tree, option_price_bsm): 
 
    plt.figure(figsize=(10, 5))
    plt.plot(steps, option_prices_binomial_tree, linestyle='-', color='steelblue', label="Binomial Tree Option Price")
    plt.axhline(y=option_price_bsm, color='orangered', linestyle='--', label="BSM Option Price")
    plt.xlabel("Number of Steps", fontsize=11)
    plt.ylabel("Option Price", fontsize=11)
    plt.legend(fontsize=12)
    plt.show()
 
def value_option_bsm(s, k, T, r, sigma, option_type="call"):
 
    # Calculate d1 and d2
    d1 = (np.log(s / k) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
 
    # Calculate option price
    if option_type == "call":
        option_price = s * norm.cdf(d1) - k * np.exp(-r * T) * norm.cdf(d2)
    else:
        option_price = k * np.exp(-r * T) * norm.cdf(-d2) - s * norm.cdf(-d1)
 
    return option_price
 
 
s = 100
k = 100
T = 1
r = 0.05
sigma = 0.2
steps = np.arange(1, 100) 
 
option_prices_binomial_tree = [binomial_tree_m_period(s, k, T, r, sigma, n) for n in steps]
option_price_bsm = value_option_bsm(s, k, T, r, sigma)
plot_tree_vs_bsm(option_prices_binomial_tree, option_price_bsm)        

GBM

import numpy as np
import matplotlib.pyplot as plt
import time
from concurrent.futures import ThreadPoolExecutor
import seaborn as sns
 
def plot_distribution(data):
     
    kde = sns.kdeplot(data, color='orange', fill=True, alpha=0.5, linewidth=2)
    plt.title('Simulated Prices Distribution', fontsize=14)
    plt.xlabel('Value', fontsize=12)
    plt.ylabel('Density', fontsize=12)
    plt.grid(False)
    plt.show()
 
def plot_simulated_prices_paths(data):
    plt.figure(figsize=(10, 6))
    for path in data:
        plt.plot(np.linspace(0, 1, len(path)), path)
    plt.title('Prices Path Simulation')
    plt.xlabel('Time(t)')
    plt.ylabel('Price')
    plt.grid(False)
    plt.show()
 
def simulate_gbm_path(s0, drift, sigma, step_t, step_n):
    steps = np.random.normal(0, np.sqrt(step_t), size=step_n)
    prices = np.zeros(step_n)
    prices[0] = s0
     
    for i in range(1, step_n):
        prices[i] = prices[i-1] * np.exp((drift-0.5*sigma**2)*step_t+sigma*steps[i-1])
     
    return prices
 
def run_simulations(s0, drift, sigma, T, step_t, step_n, num_simulations):
    with ThreadPoolExecutor() as executor:
        return [executor.submit(simulate_gbm_path, s0, drift, sigma,
                                    step_t, step_n)
                        .result() 
                            for _ in range(num_simulations)
                ]
         
 
s0 = 100 
drift = 0.1  
sigma = 0.2 
T = 1
step_t = 0.01
step_n = int(T/step_t)
num_simulations = 1000
 
simulated_paths =  run_simulations(s0, drift, sigma, T, step_t, step_n, num_simulations)
end_prices = [p[-1] for p in simulated_paths]
mean = np.mean(end_prices)
std_dev = np.std(end_prices)
 
plot_simulated_prices_paths(simulated_paths)
plot_distribution(end_prices)        

Full Code – DolphinDB

Binomial Model vs BSM

def calculate_up_down_factor(time_per_period, sigma){
 
    u = exp(sigma * sqrt(time_per_period))
    d = exp(-sigma * sqrt(time_per_period))
 
    return u, d
}
 
def calculate_up_down_probability(u, d, r, time_per_period){
     
    Pi = (exp(r*time_per_period)-d)/(u-d)
 
    return Pi, 1-Pi
}
 
def binomial_tree_m_period(s0, k, T, r, sigma, number_of_periods, option_type="call"){
 
    time_per_period = T/number_of_periods
     
    // Calculate the up and down factors
    u, d = calculate_up_down_factor(time_per_period, sigma)
     
    // calculate the up and down probabilities
    Pi_up, Pi_down = calculate_up_down_probability(u, d, r, time_per_period)
     
    // calculate underlying prices at expiratoin
    stock_prices = each(def(i): s0*pow(u, i)*pow(d, (number_of_periods-i)),
                0..number_of_periods)
 
    // calculate the option payoff at expiration
    if (option_type == "call"){
        payoffs= each(def(i): max(stock_prices[i] - k, 0), 0..number_of_periods)
    } else {
        payoffs= each(def(i): max(k - stock_prices[i], 0), 0..number_of_periods)
    }
     
    // Backward induction of option value at the initiation
    discount_factor = 1/exp(r*time_per_period)
    for (i in (number_of_periods-1)..0 ){
        for (j in 0..i){
            payoffs[j] = (Pi_up*payoffs[j+1]+Pi_down*payoffs[j])*discount_factor
        }
    }
     
    option_price = payoffs[0]
 
    return option_price
}
 
def value_option_bsm(s, k, T, r, sigma,  option_type="call"){
 
    // Calculate d1 and d2
    d1 = (log(s / k) + (r + 0.5 * pow(sigma, 2)) * T) / (sigma * sqrt(T))
    d2 = d1 - sigma * sqrt(T)
 
    // Calculate option price
    if (option_type == "call")
        option_price = s * cdfNormal(0, 1, d1) - k * exp(-r * T) * cdfNormal(0, 1, d2)
    else
        option_price = k * exp(-r * T) * cdfNormal(0, 1, -d2) - s * cdfNormal(0, 1, -d1)
 
    return option_price
}
 
def price_option_by_steps(func, s, k, T, r, sigma, steps) {
    calculation = partial(binomial_tree_m_period, s, k, T, r, sigma)
    option_prices_binomial_tree = ploop(calculation, 1..steps)
    return option_prices_binomial_tree
}
 
s = 100
k = 100
T = 1.0
r = 0.05
sigma = 0.2
steps = 100
 
results = price_option_by_steps(binomial_tree_m_period, s, k, T, r, sigma, steps)
option_price_bsm = value_option_bsm(s, k, T, r, sigma)
 
print(results)
print("Option price (BSM): "+option_price_bsm)
 
undef all;        

GBM

def simulate_gbm_path(s0, drift, sigma, step_t, step_n, run_no){
 
    steps = randNormal(0, sqrt(step_t), step_n)
    prices = array(double, step_n, step_n, 0)
    prices[0] = s0
     
    for (i in 0:(step_n-1)){
        prices[i+1] = prices[i] * exp((drift-0.5*pow(sigma,2))*step_t+sigma*steps[i])
         
    }
     
    return prices
}
 
def run_simulations(s0, drift, sigma, T, step_t, step_n, num_simulations){
 
    simulate = partial(simulate_gbm_path, s0, drift, sigma, step_t, step_n)
    results = ploop(simulate, 1..num_simulations)
 
    return results
}
         
 
s0 = 100.0
drift = 0.1  
sigma = 0.2
T = 1.0
step_t = 0.25
step_n = int(T/step_t)
num_simulations = 100
 
simulated_paths =  run_simulations(s0, drift, sigma, T, step_t, step_n, num_simulations)
print(simulated_paths)
 
undef all;