import numpy as np from scipy.stats import norm def d1(F, K, sigma, T): return (np.log(F / K) + sigma**2 * T / 2) / (sigma * np.sqrt(T)) def d2(F, K, sigma, T): return d1(F, K, sigma, T) - sigma * np.sqrt(T) def black(F, K, T, sigma, option_type = "payer"): chi = 1 if option_type == "payer" else -1 if option_type == "payer": return F * norm.cdf(d1(F, K, sigma, T)) - K * norm.cdf(d2(F, K, sigma, T)) else: return K * norm.cdf(- d2(F, K, sigma, T)) - F * norm.cdf(- d1(F, K, sigma, T))