Black-Scholes-Merton (BSM) model

bsm_model.py

from math import log, sqrt, exp
import scipy.stats as stats
import numpy as np

# from https://www.quantconnect.com/tutorials/introduction-to-options/options-pricing-black-scholes-merton-model

class BsmModel:
    
    def __init__(self, option_type, price, strike, interest_rate, expiry, volatility=None, target_price=None, dividend_yield=0):
        self.s = price  # Underlying asset price
        self.k = strike  # Option strike K
        self.r = interest_rate  # Continuous risk fee rate
        self.q = dividend_yield  # Dividend continuous rate
        self.T = expiry  # time to expiry (year)
        self.sigma = volatility  # Underlying volatility
        self.type = option_type # option type "p" put option "c" call option
        
        if target_price is not None:
            self.calc_sigma(target_price)
        
    def n(self, d):
        # cumulative probability distribution function of standard normal distribution
        return stats.norm.cdf(d)

    def dn(self, d):
        # the first order derivative of n(d)
        return stats.norm.pdf(d)

    def d1(self):
        d1 = (log(self.s / self.k) + (self.r - self.q + self.sigma ** 2 * 0.5) * self.T) / (self.sigma * sqrt(self.T))
        return d1

    def d2(self):
        d2 = (log(self.s / self.k) + (self.r - self.q - self.sigma ** 2 * 0.5) * self.T) / (self.sigma * sqrt(self.T))
        return d2

    def bsm_price(self):
        d1 = self.d1()
        d2 = d1 - self.sigma * sqrt(self.T)
        if self.type == 'c':
            price = exp(-self.r*self.T) * (self.s * exp((self.r - self.q)*self.T) * self.n(d1) - self.k * self.n(d2))
            return price
        elif self.type == 'p':
            price = exp(-self.r*self.T) * (self.k * self.n(-d2) - (self.s * exp((self.r - self.q)*self.T) * self.n(-d1)))
            return price
        else:
            print("option type can only be c or p")

    def delta(self):
        d1 = self.d1()
        if self.type == "c":
            return exp(-self.q * self.T) * self.n(d1)
        elif self.type == "p":
            return exp(-self.q * self.T) * (self.n(d1)-1)

    def gamma(self, ):
        d1 = self.d1()
        dn1 = self.dn(d1)
        return dn1 * exp(-self.q * self.T) / (self.s * self.sigma * sqrt(self.T))

    def theta(self):
        d1 = self.d1()
        d2 = d1 - self.sigma * sqrt(self.T)
        dn1 = self.dn(d1)

        if self.type == "c":
            theta = -self.s * dn1 * self.sigma * exp(-self.q*self.T) / (2 * sqrt(self.T)) \
                    + self.q * self.s * self.n(d1) * exp(-self.q*self.T) \
                    - self.r * self.k * exp(-self.r*self.T) * self.n(d2)
            return theta
        elif self.type == "p":
            theta = -self.s * dn1 * self.sigma * exp(-self.q * self.T) / (2 * sqrt(self.T)) \
                    - self.q * self.s * self.n(-d1) * exp(-self.q * self.T) \
                    + self.r * self.k * exp(-self.r * self.T) * self.n(-d2)
            return theta

    def vega(self):
        d1 = self.d1()
        dn1 = self.dn(d1)
        return self.s * sqrt(self.T) * dn1 * exp(-self.q * self.T)

    def rho(self):
        d2 = self.d2()
        if self.type == "c":
            rho = self.k * self.T * (exp(-self.r*self.T)) * self.n(d2)
        elif self.type == "p":
            rho = -self.k * self.T * (exp(-self.r*self.T)) * self.n(-d2)

        return rho
    
    def iv(self):
        return self.sigma
    
    def calc_sigma(self, target_price):
        MAX_ITERATIONS = 200
        PRECISION = 1.0e-5
        self.sigma = 0.5
        for i in range(0, MAX_ITERATIONS):
            price = self.bsm_price()
            vega = self.vega()
            diff = target_price - price  # our root
            if (abs(diff) < PRECISION):
                return self.sigma
            self.sigma = self.sigma + diff/vega # f(x) / f'(x)
        return self.sigma # value wasn't found, return best guess so far