Black-Scholes Model: Options Pricing for Quant Traders
The Black-Scholes (see our options calculator) model is the foundational framework for options pricing and risk management. Published in 1973 by Fischer Black, Myron Scholes, and Robert Merton, it provides closed-form solutions for European option prices under specific assumptions about the underlying asset's behavior. Despite its well-known limitations, Black-Scholes remains the lingua franca of options markets: traders quote in Black-Scholes implied volatility, risk is managed using Black-Scholes Greeks, and more sophisticated models are often expressed as extensions of the Black-Scholes framework. This guide implements the complete Black-Scholes toolkit, from pricing and Greeks through implied volatility calculation and the model's practical limitations.Key Takeaways
- Black-Scholes provides closed-form solutions for European calls and puts, enabling instant pricing and Greek computation.
- Implied volatility is the market's consensus about future uncertainty, extracted by inverting the Black-Scholes formula.
- The Greeks (Delta, Gamma, Theta, Vega, Rho) quantify an option's sensitivity to each input parameter.
- The model's assumptions are violated in practice, but it remains the standard framework by convention.
The Black-Scholes Formula
``python
import numpy as np
from scipy.stats import norm
from scipy.optimize import brentq
from dataclasses import dataclass
@dataclass
class BlackScholes:
"""
Black-Scholes European option pricing model.
Parameters:
S: Current stock price
K: Strike price
T: Time to expiration (years)
r: Risk-free interest rate (annual, continuous)
sigma: Volatility (annual)
q: Continuous dividend yield
"""
S: float
K: float
T: float
r: float
sigma: float
q: float = 0.0
def _d1(self) -> float:
"""Calculate d1 parameter."""
return (
np.log(self.S / self.K)
+ (self.r - self.q + 0.5 self.sigma2) self.T
) / (self.sigma * np.sqrt(self.T))
def _d2(self) -> float:
"""Calculate d2 parameter."""
return self._d1() - self.sigma * np.sqrt(self.T)
def call_price(self) -> float:
"""European call option price."""
if self.T <= 0:
return max(self.S - self.K, 0)
d1, d2 = self._d1(), self._d2()
return (
self.S np.exp(-self.q self.T) * norm.cdf(d1)
- self.K np.exp(-self.r self.T) * norm.cdf(d2)
)
def put_price(self) -> float:
"""European put option price."""
if self.T <= 0:
return max(self.K - self.S, 0)
d1, d2 = self._d1(), self._d2()
return (
self.K np.exp(-self.r self.T) * norm.cdf(-d2)
- self.S np.exp(-self.q self.T) * norm.cdf(-d1)
)
def price(self, option_type: str = "call") -> float:
"""Price call or put."""
if option_type.lower() == "call":
return self.call_price()
elif option_type.lower() == "put":
return self.put_price()
else:
raise ValueError(f"Unknown option type: {option_type}")
Example: AAPL call option
bs = BlackScholes(S=185.0, K=190.0, T=30/365, r=0.05, sigma=0.25)
print(f"Call Price: ${bs.call_price():.2f}")
print(f"Put Price: ${bs.put_price():.2f}")
Verify put-call parity: C - P = Sexp(-qT) - Kexp(-rT)
parity_lhs = bs.call_price() - bs.put_price()
parity_rhs = bs.S np.exp(-bs.q bs.T) - bs.K np.exp(-bs.r bs.T)
print(f"Put-Call Parity check: {parity_lhs:.4f} = {parity_rhs:.4f}")
`
Greeks: Sensitivity Measures
The Greeks quantify how the option price changes when each input parameter changes by a small amount.
`python
class BlackScholesGreeks(BlackScholes):
"""Extended Black-Scholes with all Greeks."""
def delta(self, option_type: str = "call") -> float:
"""
Delta: dC/dS (call) or dP/dS (put)
Measures price sensitivity to underlying price change.
"""
d1 = self._d1()
if option_type == "call":
return np.exp(-self.q self.T) norm.cdf(d1)
else:
return np.exp(-self.q self.T) (norm.cdf(d1) - 1)
def gamma(self) -> float:
"""
Gamma: d2C/dS2
Rate of change of delta. Same for calls and puts.
"""
d1 = self._d1()
return (
np.exp(-self.q self.T) norm.pdf(d1)
/ (self.S self.sigma np.sqrt(self.T))
)
def theta(self, option_type: str = "call") -> float:
"""
Theta: dC/dT (daily)
Time decay per calendar day.
"""
d1, d2 = self._d1(), self._d2()
common = (
-self.S np.exp(-self.q self.T) norm.pdf(d1) self.sigma
/ (2 * np.sqrt(self.T))
)
if option_type == "call":
theta = (
common
+ self.q self.S np.exp(-self.q self.T) norm.cdf(d1)
- self.r self.K np.exp(-self.r self.T) norm.cdf(d2)
)
else:
theta = (
common
- self.q self.S np.exp(-self.q self.T) norm.cdf(-d1)
+ self.r self.K np.exp(-self.r self.T) norm.cdf(-d2)
)
return theta / 365 # Per calendar day
def vega(self) -> float:
"""
Vega: dC/dsigma (per 1% vol change)
Same for calls and puts.
"""
d1 = self._d1()
return (
self.S np.exp(-self.q self.T) * norm.pdf(d1)
- np.sqrt(self.T) / 100 # Per 1% vol change
)
def rho(self, option_type: str = "call") -> float:
"""
Rho: dC/dr (per 1% rate change)
"""
d2 = self._d2()
if option_type == "call":
return self.K self.T np.exp(-self.r self.T) norm.cdf(d2) / 100
else:
return -self.K self.T np.exp(-self.r self.T) norm.cdf(-d2) / 100
def all_greeks(self, option_type: str = "call") -> dict:
"""Return all Greeks as a dictionary."""
return {
"price": self.price(option_type),
"delta": self.delta(option_type),
"gamma": self.gamma(),
"theta": self.theta(option_type),
"vega": self.vega(),
"rho": self.rho(option_type),
}
Calculate all Greeks
greeks = BlackScholesGreeks(S=185.0, K=190.0, T=30/365, r=0.05, sigma=0.25)
result = greeks.all_greeks("call")
print("Option Greeks:")
for name, value in result.items():
print(f" {name:>8s}: {value:+.4f}")
`
Implied Volatility
Implied volatility is found by inverting the Black-Scholes formula: given the market price, solve for sigma.
`python
class ImpliedVolatility:
"""
Compute implied volatility from market prices.
"""
@staticmethod
def solve(
market_price: float,
S: float,
K: float,
T: float,
r: float,
option_type: str = "call",
q: float = 0.0,
tol: float = 1e-8,
) -> float:
"""
Find implied volatility using Brent's method.
"""
def objective(sigma):
bs = BlackScholes(S=S, K=K, T=T, r=r, sigma=sigma, q=q)
return bs.price(option_type) - market_price
try:
iv = brentq(objective, 0.001, 5.0, xtol=tol)
return iv
except ValueError:
return np.nan
@staticmethod
def newton_raphson(
market_price: float,
S: float,
K: float,
T: float,
r: float,
option_type: str = "call",
q: float = 0.0,
max_iter: int = 100,
tol: float = 1e-8,
) -> float:
"""
Newton-Raphson method for implied volatility.
Faster than Brent for well-behaved inputs.
"""
# Initial guess using Brenner-Subrahmanyam approximation
sigma = np.sqrt(2 np.pi / T) market_price / S
for _ in range(max_iter):
bs = BlackScholesGreeks(S=S, K=K, T=T, r=r, sigma=sigma, q=q)
price = bs.price(option_type)
vega = bs.vega() * 100 # Convert back from per-1%
diff = price - market_price
if abs(diff) < tol:
return sigma
if abs(vega) < 1e-12:
break
sigma -= diff / vega
sigma = max(0.001, min(sigma, 5.0))
return sigma
Example: compute IV from market price
iv_solver = ImpliedVolatility()
market_call_price = 3.50
iv = iv_solver.solve(
market_price=market_call_price,
S=185.0, K=190.0, T=30/365, r=0.05,
option_type="call",
)
print(f"Implied Volatility: {iv:.2%}")
Verify: price at IV should match market
bs_check = BlackScholes(S=185.0, K=190.0, T=30/365, r=0.05, sigma=iv)
print(f"Model Price at IV: ${bs_check.call_price():.2f}")
print(f"Market Price: ${market_call_price:.2f}")
`
Option Chain Analysis
Analyze a full option chain to extract the implied volatility surface.
`python
def analyze_option_chain(
chain: pd.DataFrame,
spot: float,
r: float = 0.05,
q: float = 0.0,
) -> pd.DataFrame:
"""
Compute implied volatilities and Greeks for an option chain.
Input columns: strike, expiry, mid_price, option_type
"""
iv_solver = ImpliedVolatility()
results = []
for _, row in chain.iterrows():
T = (pd.to_datetime(row["expiry"]) - pd.Timestamp.now()).days / 365
if T <= 0:
continue
iv = iv_solver.solve(
market_price=row["mid_price"],
S=spot, K=row["strike"], T=T, r=r,
option_type=row["option_type"], q=q,
)
if np.isnan(iv):
continue
greeks = BlackScholesGreeks(
S=spot, K=row["strike"], T=T, r=r, sigma=iv, q=q
)
greek_values = greeks.all_greeks(row["option_type"])
moneyness = np.log(row["strike"] / spot) / (iv * np.sqrt(T))
results.append({
"strike": row["strike"],
"expiry": row["expiry"],
"type": row["option_type"],
"mid_price": row["mid_price"],
"implied_vol": iv,
"moneyness": moneyness,
"T": T,
**greek_values,
})
return pd.DataFrame(results)
`
Portfolio Greeks Management
`python
class OptionsPortfolio:
"""Manage Greeks at the portfolio level."""
def __init__(self):
self.positions = []
def add_position(
self,
S: float, K: float, T: float, r: float, sigma: float,
option_type: str, quantity: int, q: float = 0.0,
):
"""Add an option position to the portfolio."""
self.positions.append({
"S": S, "K": K, "T": T, "r": r, "sigma": sigma,
"type": option_type, "quantity": quantity, "q": q,
})
def portfolio_greeks(self) -> dict:
"""Compute aggregate portfolio Greeks."""
total = {"delta": 0, "gamma": 0, "theta": 0, "vega": 0, "rho": 0, "value": 0}
for pos in self.positions:
greeks = BlackScholesGreeks(
S=pos["S"], K=pos["K"], T=pos["T"],
r=pos["r"], sigma=pos["sigma"], q=pos["q"],
)
g = greeks.all_greeks(pos["type"])
total["delta"] += g["delta"] * pos["quantity"]
total["gamma"] += g["gamma"] * pos["quantity"]
total["theta"] += g["theta"] * pos["quantity"]
total["vega"] += g["vega"] * pos["quantity"]
total["rho"] += g["rho"] * pos["quantity"]
total["value"] += g["price"] * pos["quantity"]
return total
def delta_hedge_shares(self) -> float:
"""Number of shares needed to delta-hedge the portfolio."""
greeks = self.portfolio_greeks()
return -greeks["delta"] * 100 # Options are on 100 shares
Example: construct a portfolio
portfolio = OptionsPortfolio()
portfolio.add_position(S=185, K=190, T=30/365, r=0.05, sigma=0.25,
option_type="call", quantity=10)
portfolio.add_position(S=185, K=180, T=30/365, r=0.05, sigma=0.25,
option_type="put", quantity=-5)
greeks = portfolio.portfolio_greeks()
print("Portfolio Greeks:")
for k, v in greeks.items():
print(f" {k}: {v:+.4f}")
print(f"\nShares to delta-hedge: {portfolio.delta_hedge_shares():.0f}")
`
Black-Scholes Limitations
Understanding where Black-Scholes fails is as important as knowing how to use it.
`python
def demonstrate_bs_limitations():
"""
Key assumptions violated in practice:
1. Constant volatility (reality: volatility smile/skew)
2. Log-normal returns (reality: fat tails, negative skew)
3. Continuous trading (reality: discrete, gaps)
4. No transaction costs (reality: bid-ask, impact)
5. Known constant interest rate (reality: stochastic)
"""
# Demonstrate volatility smile
spot = 100
strikes = np.arange(80, 121, 2)
# Simulated market IVs (showing skew)
atm_vol = 0.20
skew = -0.001 # Negative skew: OTM puts have higher IV
market_ivs = atm_vol + skew * (strikes - spot)
print("Volatility Skew (Black-Scholes assumes flat line):")
print(f"{'Strike':>8s} {'IV':>8s} {'Moneyness':>12s}")
for k, iv in zip(strikes, market_ivs):
moneyness = "ITM" if k < spot else "ATM" if k == spot else "OTM"
print(f"{k:>8.0f} {iv:>8.1%} {moneyness:>12s}")
``