The Black-Scholes model is a cornerstone in financial mathematics, providing a theoretical framework for pricing European call and put options. A clean Python class can effectively encapsulate the model’s parameters and calculations, offering a structured approach to determine the price P(t) of an option. This class typically includes attributes for the strike price K, volatility σ, maturity T, and the underlying stock price S. By leveraging Python’s object-oriented capabilities, the class can efficiently compute the option’s price using the Black-Scholes formula, ensuring accuracy and ease of use for financial analysts and developers.
Python Class for the price P of an equity option with strike K, volatility σ and maturity T of an underlying stock S
Implementing such a class involves defining methods that calculate the necessary components of the Black-Scholes formula, such as the cumulative distribution function and the standard normal distribution. These methods allow for the dynamic computation of option prices as market conditions change. Additionally, the class can be extended to include error handling and input validation, ensuring that users provide valid and realistic parameters. This approach not only enhances the reliability of the model but also makes it accessible to those who may not have a deep understanding of the underlying mathematics.
import math
from dataclasses import dataclass
from typing import Literal
from scipy.stats import norm
@dataclass
class EquityOption:
S: float # underlying price
K: float # strike
sigma: float # volatility
T: float # maturity (in years)
r: float = 0.0 # risk-free rate
option_type: Literal["call", "put"] = "call"
def price(self, t: float = 0.0) -> float:
"""
Returns Black-Scholes price P(t) at time t.
"""
tau = self.T - t # time to maturity
if tau <= 0:
# Option has expired
if self.option_type == "call":
return max(self.S - self.K, 0)
else:
return max(self.K - self.S, 0)
d1 = (
math.log(self.S / self.K)
+ (self.r + 0.5 * self.sigma**2) * tau
) / (self.sigma * math.sqrt(tau))
d2 = d1 - self.sigma * math.sqrt(tau)
if self.option_type == "call":
return self.S * norm.cdf(d1) - self.K * math.exp(-self.r * tau) * norm.cdf(d2)
else:
return self.K * math.exp(-self.r * tau) * norm.cdf(-d2) - self.S * norm.cdf(-d1)
usage:
opt = EquityOption(S=100, K=110, sigma=0.2, T=1, r=0.03, option_type="call")
price_now = opt.price(t=0)
print(price_now)
Advanced Features and Customisation Options
Beyond the basic implementation, a Python class for option pricing can be enhanced with advanced features that cater to more sophisticated financial strategies. For instance, the class can be extended to support the calculation of Greeks, which are derivatives of the option price with respect to various parameters. These metrics, such as delta, gamma, and vega, provide valuable insights into the sensitivity of the option’s price to changes in market conditions. By incorporating these features, the class becomes a powerful tool for risk management and strategic decision-making in trading environments.
Customisation options can further enhance the utility of the class, allowing users to tailor the model to specific market scenarios or personal preferences. This could include the ability to input historical volatility data, adjust for dividends, or simulate different market conditions. By providing a flexible and adaptable framework, the class can accommodate a wide range of use cases, from academic research to real-world trading applications. Such versatility ensures that the model remains relevant and useful in the ever-evolving landscape of financial markets.
Historical Pre-Amble: Market Observations and empirical facts
An examination of market prices of an ATM call compared to OTM (K > ATM) and ITM (K > ATM) calls on reveals interesting facts. In 1974 we claimed that in the Black Scholes model, the volatility does not depend on the strike. Merton claimed that the BS pricing formula will become a self-fulling prophecy.
In 1980 we see the market smiling.
In 1987 we have the market crash (in Oct 1987).
| Year | Author | Model Name(s) | SDE | Martingale | IR |
|---|---|---|---|---|---|
| 1974 | |||||
With the underlying stock S modelled as time-dependent S(t)
Self-contained Python class that models price P(t) of an equity option with strike K, volatility σ, maturity T, and a time-dependent underlying process S(t).
It includes:
- Storage of S(t) as a callable time-dependent function or array
- Black–Scholes price for European calls/puts
- Greeks (Delta, Gamma, Vega, Theta, Rho)
- Path simulation under geometric Brownian motion
- American option pricing using Longstaff–Schwartz Monte Carlo (LSM)
import numpy as np
from scipy.stats import norm
class EquityOption:
def __init__(self, S_func, K, sigma, r, T, option_type="call"):
"""
S_func(t) : function returning underlying price S(t)
(e.g., lambda t: S0*np.exp(...)
or a dict/array interpolator)
K : strike
sigma : volatility
r : risk-free rate
T : maturity
option_type : "call" or "put"
"""
self.S_func = S_func
self.K = K
self.sigma = sigma
self.r = r
self.T = T
self.option_type = option_type.lower()
# ---------------------------
# Black–Scholes price
# ---------------------------
def price(self, t=0):
S = self.S_func(t)
tau = self.T - t
if tau <= 0: return max(0, (S - self.K) if self.option_type=="call"
else (self.K - S)) d1 = (np.log(S/self.K) + (self.r + 0.5*self.sigma**2)*tau) /
(self.sigma*np.sqrt(tau)) d2 = d1 - self.sigma*np.sqrt(tau)
if self.option_type == "call": return S*norm.cdf(d1) - self.K*np.exp(-self.r*tau)*norm.cdf(d2)
else:
return self.K*np.exp(-self.r*tau)*norm.cdf(-d2) - S*norm.cdf(-d1)
# ---------------------------
# Greeks
# ---------------------------
def greeks(self, t=0):
S = self.S_func(t) tau = self.T - t
d1 = (np.log(S/self.K) + (self.r + 0.5*self.sigma**2)*tau) / (self.sigma*np.sqrt(tau))
d2 = d1 - self.sigma*np.sqrt(tau) pdf = norm.pdf(d1) Delta = norm.cdf(d1)
if self.option_type=="call"
else -norm.cdf(-d1) Gamma = pdf / (S*self.sigma*np.sqrt(tau))
Vega = S * pdf * np.sqrt(tau)
Theta = (-S*pdf*self.sigma/(2*np.sqrt(tau)) - self.r*self.K*np.exp(-self.r*tau) * (norm.cdf(d2)
if self.option_type=="call"
else norm.cdf(-d2))) Rho = self.K * tau * np.exp(-self.r*tau) * (norm.cdf(d2)
if self.option_type=="call"
else -norm.cdf(-d2))
return {
"Delta": Delta,
"Gamma": Gamma,
"Vega": Vega,
"Theta": Theta,
"Rho": Rho
}
# ---------------------------
# Simulate underlying paths
# ---------------------------
def simulate_paths(self, S0, n_steps=100, n_paths=10000):
"""
Simulate geometric Brownian motion paths for S(t)
"""
dt = self.T / n_steps S = np.zeros((n_paths, n_steps+1))
S[:,0] = S0 for i in range(1, n_steps+1):
z = np.random.randn(n_paths)
S[:,i] = S[:,i-1] * np.exp( (self.r - 0.5*self.sigma**2)*dt + self.sigma*np.sqrt(dt)*z )
return S
# ---------------------------
# American option via LSM
# ---------------------------
def american_price_LSM(self, S0, n_steps=50, n_paths=20000, poly_order=2):
dt = self.T / n_steps S = self.simulate_paths(S0, n_steps, n_paths)
# Payoff function
if self.option_type == "call":
payoff = lambda s: np.maximum(s - self.K, 0)
else:
payoff = lambda s: np.maximum(self.K - s, 0) cashflows = payoff(S[:, -1])
# Backward induction
for t in range(n_steps-1, 0, -1): itm = payoff(S[:, t]) > 0
X = S[itm, t]
Y = cashflows[itm] * np.exp(-self.r*dt)
# Regression continuation value
A = np.vander(X, poly_order + 1)
coeffs = np.linalg.lstsq(A, Y, rcond=None)[0]
continuation = A @ coeffs
# Exercise decision
exercise = payoff(X) > continuation
cashflows[itm] = np.where(exercise, payoff(X), cashflows[itm] * np.exp(-self.r*dt))
# Initial price
return np.mean(cashflows) * np.exp(-self.r*dt)
usage
# S(t) as a function
S0 = 100
S_t = lambda t: S0 * np.exp(0.05 * t)
opt = EquityOption(S_func=S_t, K=100, sigma=0.2, r=0.05, T=1)
print("European price:", opt.price())
print("Greeks:", opt.greeks())
# Path simulation
paths = opt.simulate_paths(S0=100)
# American option
print("American price:", opt.american_price_LSM(S0=100))🧩 Python Class for Equity Option Pricing (with Stochastic Vol)
Below is a clean, extensible Python class that represents the price of an equity option with
- strike K
- maturity
- stochastic volatility
- time-dependent underlying
- ability to compute Greeks
- ability to simulate paths
- ability to handle American options (via Longstaff–Schwartz)
This is a framework you can extend; full production-grade implementations are long, but this gives you a correct, consistent structure with working components.
import numpy as np
from dataclasses import dataclass
@dataclass
class OptionParams:
K: float # strike
T: float # maturity
r: float # risk-free rate
is_call: bool # True for call, False for put
class StochasticVolOption:
"""
Represents an equity option with:
- time-dependent underlying S(t)
- stochastic volatility (Heston dynamics)
- Greek estimation
- path simulation (Euler or QE)
- American option valuation via Longstaff–Schwartz
"""
def __init__(self, S0, option: OptionParams,
v0=0.04, kappa=2.0, theta=0.04, xi=0.5, rho=-0.5):
"""
Heston parameters:
dS = r*S dt + sqrt(v) S dW1
dv = kappa*(theta - v) dt + xi sqrt(v) dW2
corr(W1, W2) = rho
"""
self.S0 = S0
self.opt = option
# stochastic vol parameters
self.v0 = v0
self.kappa = kappa
self.theta = theta
self.xi = xi
self.rho = rho
# ----------------------------------
# Stochastic Vol Path Simulation
# ----------------------------------
def simulate_paths(self, N=20000, steps=200, seed=None):
if seed is not None:
np.random.seed(seed)
dt = self.opt.T / steps
S = np.zeros((N, steps+1))
v = np.zeros((N, steps+1))
S[:, 0] = self.S0
v[:, 0] = self.v0
# Correlated normals
Z1 = np.random.normal(size=(N, steps))
Z2 = np.random.normal(size=(N, steps))
Z2 = self.rho*Z1 + np.sqrt(1 - self.rho**2)*Z2
for t in range(steps):
vt = np.maximum(v[:, t], 0)
# Heston variance process (Euler)
v[:, t+1] = vt + self.kappa*(self.theta - vt)*dt + \
self.xi*np.sqrt(vt*np.maximum(dt, 1e-10)) * Z2[:, t]
v[:, t+1] = np.maximum(v[:, t+1], 1e-12)
# Stock price
S[:, t+1] = S[:, t] * np.exp(
(self.opt.r - 0.5*vt)*dt + np.sqrt(vt*dt)*Z1[:, t]
)
return S, v
# ----------------------------------
# European Option Price (MC)
# ----------------------------------
def price_european(self, N=30000, steps=200):
S, _ = self.simulate_paths(N, steps)
ST = S[:, -1]
K = self.opt.K
payoff = (ST - K if self.opt.is_call else K - ST)
payoff = np.maximum(payoff, 0)
return np.exp(-self.opt.r*self.opt.T) * payoff.mean()
# ----------------------------------
# Greeks (MC Pathwise)
# ----------------------------------
def delta(self, bump=1e-4, **kwargs):
p_plus = self.shift_price(self.S0 + bump, **kwargs)
p_minus = self.shift_price(self.S0 - bump, **kwargs)
return (p_plus - p_minus) / (2*bump)
def vega(self, bump=1e-4, **kwargs):
p_plus = self.shift_price(self.S0, v0=self.v0 + bump, **kwargs)
p_minus = self.shift_price(self.S0, v0=self.v0 - bump, **kwargs)
return (p_plus - p_minus) / (2*bump)
def shift_price(self, S0_new, v0=None, **kwargs):
tmp = StochasticVolOption(
S0=S0_new,
option=self.opt,
v0=self.v0 if v0 is None else v0,
kappa=self.kappa,
theta=self.theta,
xi=self.xi,
rho=self.rho
)
return tmp.price_european(**kwargs)
# ----------------------------------
# American Option (Longstaff–Schwartz)
# ----------------------------------
def price_american(self, N=30000, steps=200):
S, _ = self.simulate_paths(N, steps)
K = self.opt.K
dt = self.opt.T / steps
r = self.opt.r
if self.opt.is_call:
payoff = lambda x: np.maximum(x - K, 0)
else:
payoff = lambda x: np.maximum(K - x, 0)
cashflow = payoff(S[:, -1])
# Backward induction
for t in reversed(range(1, steps)):
itm = payoff(S[:, t]) > 0
St_itm = S[itm, t]
cf_itm = cashflow[itm] * np.exp(-r*dt)
if len(St_itm) > 0:
# regression basis = [1, S, S^2]
X = np.vstack([np.ones_like(St_itm), St_itm, St_itm**2]).T
coeffs = np.linalg.lstsq(X, cf_itm, rcond=None)[0]
continuation = X @ coeffs
exercise = payoff(St_itm)
# exercise if immediate payoff >= continuation
cashflow[itm] = np.where(exercise >= continuation, exercise, cf_itm)
else:
cashflow = cashflow * np.exp(-r*dt)
return np.exp(-r*dt) * cashflow.mean()