ブラックショールズ方程式のライブラリーです。自動取引の支援に使う予定です。基本的には一番下にある出典に全部載っています。

# -*- coding: utf-8 -*-
"""
Created on Sat Feb 10 13:13:31 2024

@author: Lune 2121

IB/option_utils.py

"""
import numpy as np
from scipy.stats import norm

#https://www.codearmo.com/python-tutorial/options-trading-black-scholes-model
#https://www.macroption.com/black-scholes-formula/


"""
right = 'C' or 'P'
s = Spot Price
k =  Strike Price
t = Days to expiration
rfr = Risk-free Rate
sigma = volatility
div = Annual dividend rate. Defaulted to zero.
price = Known option price. Necessary for implied_volatility function
"""


def d_one(s, k, t, rfr, sigma, div=0):
    """d1 calculation"""
    d_1 = (np.log(s / k) +
       (rfr - div + sigma ** 2 / 2) * t) / (sigma * np.sqrt(t))
    return d_1


def d_two(s, k, t, rfr, sigma, div=0):
    """d2 calculation"""
    d_2 = d_one(s, k, t, rfr, sigma, div) - sigma * np.sqrt(t)
    return d_2


def nd_one(right, s, k, t, rfr, sigma, div=0):
    """nd1 calculation"""
    if right == 'C':
        nd_1 = norm.cdf(d_one(s, k, t, rfr, sigma, div), 0, 1)
    elif right == 'P':
        nd_1 = norm.cdf(-d_one(s, k, t, rfr, sigma, div), 0, 1)
    return nd_1


def nd_two(right, s, k, t, rfr, sigma, div=0):
    """nd2 calculation"""
    if right == 'C':
        nd_2 = norm.cdf(d_two(s, k, t, rfr, sigma, div), 0, 1)
    elif right == 'P':
        nd_2 = norm.cdf(-d_two(s, k, t, rfr, sigma, div), 0, 1)
    return nd_2


def option_price(right, s, k, t, rfr, sigma, div=0):
    """option price"""
    right = right.upper()
    t /= 365.0
    s = abs(s)
    k = abs(k)
    if right == 'C':
        if t <= 0:
            if s > k:
                price = s - k
            else:
                price = 0
        else:
            price = (s * np.exp(-div * t) *
                 nd_one(right, s, k, t, rfr, sigma, div)
                 - k * np.exp(-rfr * t) *
                 nd_two(right, s, k, t, rfr, sigma, div))
    elif right == 'P':
        if t <= 0:
            if k > s:
                price = k - s
            else:
                price = 0
        else:            
            price = (k * np.exp(-rfr * t) *
                 nd_two(right, s, k, t, rfr, sigma, div)
                 - s * np.exp(-div * t) *
                 nd_one(right, s, k, t, rfr, sigma, div))
    return price

def option_delta(right, s, k, t, rfr, sigma, div=0):
    """option delta"""
    t /= 365
    if right == 'C':
        delta = np.exp(-div * t) * norm.cdf(d_one(s,k,t,rfr,sigma,div))
    else:
        delta = - np.exp(-div * t) * norm.cdf(-d_one(s,k,t,rfr,sigma,div))
    return delta         


def option_gamma(right, s, k, t, rfr, sigma, div=0):
    """option gamma"""
    t /= 365
    gamma = np.exp(-div*t)/s/sigma/np.sqrt(t) * norm.pdf(d_one(s,k,t,rfr,sigma,div))
    return gamma


def option_theta(right, s, k, t, rfr, sigma, div=0): #r, S, K, T, sigma, type="c"):
    """option theta"""
    t /= 365
    if right == "C":
        thetat = -s*np.exp(-div*t)*norm.pdf(d_one(s,k,t,rfr,sigma,div), 0, 1)*sigma/(2*np.sqrt(t)) \
            - rfr*k*np.exp(-rfr*t)*norm.cdf(d_two(s,k,t,rfr,sigma), 0, 1) \
            + div*s*np.exp(-div*t)*norm.cdf(d_one(s,k,t,rfr,sigma,div),0,1)
    else:
        thetat = -s*np.exp(-div*t)*norm.pdf(d_one(s,k,t,rfr,sigma,div), 0, 1)*sigma/(2*np.sqrt(t)) \
            + rfr*k*np.exp(-rfr*t)*norm.cdf(-d_two(s,k,t,rfr,sigma,div), 0, 1) \
            - div*s*np.exp(-div*t)*norm.cdf(-d_one(s,k,t,rfr,sigma,div),0,1)
    theta = thetat/365
    return theta
        
def option_vega(s, k, t, rfr, sigma, div=0):
    """option vega"""
    t /= 365
    vega = (.01 * s * np.exp(-div * t) * np.sqrt(t)
            * norm.pdf(d_one(s, k, t, rfr, sigma, div)))
    return vega


def implied_vol(right, s, k, t, rfr, market_price, div=0, tol=0.00001):
    """Compute the implied volatility of a European Option
        s: stock price
        k:  strike
        t:  maturity
        rfr:  risk-free rate
        market_price: market observed price
        tol: tolerance
    """
    max_iter = 200 #max number of iterations
    vol_old = 0.30 #initial guess
    for i in range(max_iter):
        bs_price = option_price(right, s, k, t, rfr, vol_old, div)
        Cprime = option_vega(s,k,t,rfr,vol_old, div)*100
        C = bs_price - market_price
        vol_new = vol_old - C/Cprime
        bs_new = option_price(right, s, k, t, rfr, vol_new, div)
        if (abs(vol_old - vol_new) < tol or abs(bs_new - market_price) < tol):
            break
        vol_old = vol_new
    implied_vol = vol_old
    return implied_vol

基本下のサイトを参考に作りました。