"""
Basket option pricing via Monte Carlo simulation.
Provides:
- simulate_basket: Multi-asset GBM simulation for basket option pricing
- build_pricing_grid: Construct a multi-dimensional pricing tensor
for Tensor-Train compression
"""
from __future__ import annotations
import numpy as np
[docs]
def simulate_basket(
S0: np.ndarray,
K: float,
T: float,
r: float,
sigma: np.ndarray,
corr: np.ndarray,
weights: np.ndarray,
n_paths: int = 100_000,
q: np.ndarray | None = None,
option_type: str = "call",
seed: int | None = None,
) -> tuple[float, float]:
"""Price a basket option via Monte Carlo simulation.
Uses correlated geometric Brownian motion (GBM) under the
risk-neutral measure to simulate terminal asset prices.
Payoff_call = max(sum(w_i * S_i(T)) - K, 0)
Price = exp(-rT) * E[Payoff]
Args:
S0: Initial spot prices, shape (d,).
K: Strike price.
T: Time to expiry (years).
r: Risk-free rate.
sigma: Volatilities, shape (d,).
corr: Correlation matrix, shape (d, d).
weights: Basket weights, shape (d,).
n_paths: Number of Monte Carlo paths.
q: Dividend yields, shape (d,). Defaults to zero.
option_type: 'call' or 'put'.
seed: Random seed for reproducibility.
Returns:
(price, stderr): Discounted expected payoff and its standard error.
"""
S0 = np.asarray(S0, dtype=float)
sigma = np.asarray(sigma, dtype=float)
corr = np.asarray(corr, dtype=float)
weights = np.asarray(weights, dtype=float)
d = len(S0)
if option_type not in ("call", "put"):
raise ValueError(f"option_type must be 'call' or 'put', got {option_type!r}")
if np.any(S0 <= 0):
raise ValueError("All initial spot prices S0 must be positive")
if K <= 0:
raise ValueError("Strike K must be positive")
if T <= 0:
raise ValueError("Time to expiry T must be positive")
if np.any(sigma <= 0):
raise ValueError("All volatilities must be positive")
if sigma.shape != (d,):
raise ValueError(f"sigma shape {sigma.shape} doesn't match S0 length {d}")
if corr.shape != (d, d):
raise ValueError(f"corr shape {corr.shape} doesn't match (d,d)=({d},{d})")
if weights.shape != (d,):
raise ValueError(f"weights shape {weights.shape} doesn't match S0 length {d}")
if n_paths < 1:
raise ValueError("n_paths must be >= 1")
if q is None:
q = np.zeros(d)
rng = np.random.default_rng(seed)
# Cholesky decomposition of correlation matrix
# Regularize if needed for numerical stability
try:
L = np.linalg.cholesky(corr)
except np.linalg.LinAlgError:
corr_reg = corr + 1e-8 * np.eye(d)
L = np.linalg.cholesky(corr_reg)
# Generate correlated standard normals: (n_paths, d)
Z = rng.standard_normal((n_paths, d))
Z_corr = Z @ L.T
# Simulate terminal prices under risk-neutral measure
# S_i(T) = S_i(0) * exp((r - q_i - sigma_i^2/2)*T + sigma_i*sqrt(T)*Z_i)
drift = (r - q - 0.5 * sigma**2) * T # shape (d,)
diffusion = sigma * np.sqrt(T) * Z_corr # shape (n_paths, d)
S_T = S0 * np.exp(drift + diffusion) # shape (n_paths, d)
# Basket value at maturity
basket = S_T @ weights # shape (n_paths,)
# Payoff
if option_type == "call":
payoff = np.maximum(basket - K, 0.0)
else:
payoff = np.maximum(K - basket, 0.0)
# Discounted price and standard error
discount = np.exp(-r * T)
price = discount * np.mean(payoff)
stderr = discount * np.std(payoff) / np.sqrt(n_paths)
return float(price), float(stderr)
def _price_at_spots(
spots: np.ndarray,
K: float,
T: float,
r: float,
sigma: np.ndarray,
corr: np.ndarray,
weights: np.ndarray,
n_mc_paths: int,
q: np.ndarray | None,
option_type: str,
seed: int,
) -> float:
"""Price a single basket option for given spot values."""
price, _ = simulate_basket(
S0=spots,
K=K,
T=T,
r=r,
sigma=sigma,
corr=corr,
weights=weights,
n_paths=n_mc_paths,
q=q,
option_type=option_type,
seed=seed,
)
return price
[docs]
def build_pricing_grid(
S0_ranges: list[tuple[float, float]],
K: float,
T: float,
r: float,
sigma: np.ndarray,
corr: np.ndarray,
weights: np.ndarray,
n_points: int = 30,
n_mc_paths: int = 1_000,
q: np.ndarray | None = None,
option_type: str = "call",
seed: int = 42,
) -> tuple[np.ndarray, list[np.ndarray]]:
"""Build a multi-dimensional pricing tensor for TT compression.
For each combination of spot prices on the grid, runs a Monte Carlo
simulation to compute the basket option price. Returns a d-dimensional
tensor of shape (n_points, n_points, ..., n_points).
Args:
S0_ranges: List of (min, max) tuples for each asset's spot range.
K: Strike price.
T: Time to expiry.
r: Risk-free rate.
sigma: Volatilities, shape (d,).
corr: Correlation matrix, shape (d, d).
weights: Basket weights, shape (d,).
n_points: Number of grid points per asset.
n_mc_paths: MC paths per grid point (use fewer for speed).
q: Dividend yields.
option_type: 'call' or 'put'.
seed: Base random seed.
Returns:
(grid_tensor, axes): The pricing tensor and list of 1D grid axes.
"""
d = len(S0_ranges)
if q is None:
q = np.zeros(d)
# Create grid axes for each asset
axes = [np.linspace(lo, hi, n_points) for lo, hi in S0_ranges]
# Create meshgrid indices
grid_shape = tuple([n_points] * d)
grid = np.zeros(grid_shape)
# Iterate over all grid points
# Use np.ndindex for clean iteration over d-dimensional grid
n_points**d
for flat_idx, multi_idx in enumerate(np.ndindex(*grid_shape)):
# Construct spot vector for this grid point
spots = np.array([axes[i][multi_idx[i]] for i in range(d)])
# Use a deterministic seed per grid point for reproducibility
point_seed = seed + flat_idx
grid[multi_idx] = _price_at_spots(
spots,
K,
T,
r,
sigma,
corr,
weights,
n_mc_paths,
q,
option_type,
point_seed,
)
return grid, axes
[docs]
def build_pricing_grid_analytic(
S0_ranges: list[tuple[float, float]],
K: float,
T: float,
r: float,
sigma: np.ndarray,
weights: np.ndarray,
n_points: int = 30,
) -> tuple[np.ndarray, list[np.ndarray]]:
"""Build a pricing grid using log-normal moment-matching (Gentle 1993).
Approximates the basket option price by matching the first two moments
of the basket value to a single log-normal distribution, then applying
the Black-Scholes formula to that equivalent asset. This is the standard
"moment-matching" approximation used in practice for European basket options.
The approximation is accurate to within ~1-3% for typical parameters
and produces a smooth, low-rank tensor ideal for TT compression.
For deep ITM/OTM or highly correlated assets (rho > 0.9), consider
using ``build_pricing_grid`` (Monte Carlo) for higher accuracy.
Args:
S0_ranges: List of (min, max) tuples for each asset's spot range.
K: Strike price.
T: Time to expiry.
r: Risk-free rate.
sigma: Volatilities per asset, shape (d,).
weights: Basket weights, shape (d,).
n_points: Grid points per axis.
Returns:
(grid_tensor, axes): Pricing tensor and grid axes.
"""
from scipy.stats import norm as _norm
d = len(S0_ranges)
sigma = np.asarray(sigma, dtype=float)
weights = np.asarray(weights, dtype=float)
axes = [np.linspace(lo, hi, n_points) for lo, hi in S0_ranges]
grids = np.meshgrid(*axes, indexing="ij") # d arrays each shape (n,...,n)
# Forward prices for each asset at each grid point
# F_i = S_i * exp(r * T)
forwards = [grids[i] * np.exp(r * T) for i in range(d)]
# First moment of basket: E[B] = sum_i w_i * F_i
E_B = sum(weights[i] * forwards[i] for i in range(d))
# Second moment via independence assumption (worst case: zero correlation)
# E[B^2] = sum_i sum_j w_i * w_j * F_i * F_j * exp(sigma_i^2 * T if i==j else 0)
# For the diagonal (i==j): E[S_i^2] = F_i^2 * exp(sigma_i^2 * T)
E_B2 = np.zeros_like(E_B)
for i in range(d):
for j in range(d):
if i == j:
E_B2 += (
weights[i] * weights[j] * forwards[i] * forwards[j] * np.exp(sigma[i] ** 2 * T)
)
else:
E_B2 += weights[i] * weights[j] * forwards[i] * forwards[j]
# Equivalent lognormal: match E[B] and E[B^2]
# If X ~ LN(mu_X, sigma_X^2) then E[X] = exp(mu_X + 0.5*sigma_X^2)
# and E[X^2] = exp(2*mu_X + 2*sigma_X^2)
# => sigma_X^2 = log(E[B^2]) - 2*log(E[B])
# => mu_X = log(E[B]) - 0.5 * sigma_X^2
sigma_X2 = np.log(np.maximum(E_B2, 1e-15)) - 2.0 * np.log(np.maximum(E_B, 1e-15))
sigma_X = np.sqrt(np.maximum(sigma_X2, 1e-12))
F_X = E_B # forward price of the equivalent asset equals E[B]
# Black-76 call price: C = exp(-rT) * (F * N(d1) - K * N(d2))
# d1 = (log(F/K) + 0.5*sigma_X^2*T) / (sigma_X * sqrt(T))
# d2 = d1 - sigma_X * sqrt(T)
sqrtT = np.sqrt(T)
log_FK = np.log(np.maximum(F_X, 1e-15) / K)
d1 = (log_FK + 0.5 * sigma_X2 * T) / np.maximum(sigma_X * sqrtT, 1e-12)
d2 = d1 - sigma_X * sqrtT
N_d1 = _norm.cdf(d1)
N_d2 = _norm.cdf(d2)
discount = np.exp(-r * T)
grid = discount * (F_X * N_d1 - K * N_d2)
# Where intrinsic value exceeds BS value (shouldn't happen but clip for safety)
intrinsic = discount * np.maximum(E_B - K, 0.0)
grid = np.maximum(grid, intrinsic)
return grid, axes
