Quant Companion MCP

/** * Heston Stochastic Volatility Model * * Implements Heston Monte Carlo simulation and basic calibration. * * Model dynamics: * dS = r*S*dt + sqrt(V)*S*dW1 * dV = kappa*(theta - V)*dt + xi*sqrt(V)*dW2 * dW1*dW2 = rho*dt */ import { mean, stdDev, percentile } from "./utils"; export interface HestonParams { /** Current spot price */ spot: number; /** Risk-free rate */ rate: number; /** Mean reversion speed */ kappa: number; /** Long-term variance */ theta: number; /** Volatility of volatility */ xi: number; /** Correlation between spot and variance */ rho: number; /** Initial variance */ v0: number; /** Time to maturity in years */ timeToMaturityYears: number; /** Number of time steps */ steps: number; /** Number of simulation paths */ paths: number; /** Dividend yield (optional) */ dividendYield?: number; } export interface HestonSimResult { /** Terminal prices from simulation */ terminalPrices: number[]; /** Terminal variances from simulation */ terminalVariances: number[]; /** Average variance along paths */ avgVariance: number; } export interface HestonPriceResult { /** Option price */ price: number; /** Standard error of estimate */ standardError: number; /** 95% confidence interval */ confidenceInterval: { lower: number; upper: number }; } export interface HestonCalibrationInput { /** Forward prices for each maturity */ forwards: number[]; /** Maturities in years */ maturitiesYears: number[]; /** Strikes (same for all maturities, or 2D array) */ strikes: number[] | number[][]; /** Market implied vols [maturity][strike] */ marketIvs: number[][]; /** Initial guesses (optional) */ initialGuess?: Partial<{ kappa: number; theta: number; xi: number; rho: number; v0: number; }>; } export interface HestonCalibrationResult { kappa: number; theta: number; xi: number; rho: number; v0: number; /** Root mean square error */ error: number; /** Whether calibration converged */ converged: boolean; } /** * Simulate Heston model paths using Euler-Maruyama with full truncation * * Uses the QE (Quadratic Exponential) scheme for variance to ensure positivity */ export function simulateHestonPaths(params: HestonParams): HestonSimResult { const { spot, rate, kappa, theta, xi, rho, v0, timeToMaturityYears: T, steps, paths, dividendYield = 0, } = params; if (paths > 200000) { throw new Error("paths must be <= 200000"); } const dt = T / steps; const sqrtDt = Math.sqrt(dt); const drift = rate - dividendYield; // Correlation structure const sqrtOneMinusRhoSq = Math.sqrt(1 - rho * rho); const terminalPrices: number[] = []; const terminalVariances: number[] = []; let totalVariance = 0; for (let i = 0; i < paths; i++) { let S = spot; let V = v0; for (let j = 0; j < steps; j++) { // Generate correlated Brownian increments const z1 = randomStdNormal(); const z2 = rho * z1 + sqrtOneMinusRhoSq * randomStdNormal(); // Ensure variance stays positive (full truncation) const sqrtV = Math.sqrt(Math.max(V, 0)); // Update spot S = S * Math.exp((drift - 0.5 * Math.max(V, 0)) * dt + sqrtV * sqrtDt * z1); // Update variance (Euler with reflection) const dV = kappa * (theta - V) * dt + xi * sqrtV * sqrtDt * z2; V = Math.abs(V + dV); // Reflection to keep positive totalVariance += V; } terminalPrices.push(S); terminalVariances.push(V); } return { terminalPrices, terminalVariances, avgVariance: totalVariance / (paths * steps), }; } /** * Price a European option using Heston Monte Carlo */ export function priceOptionHestonMonteCarlo( params: HestonParams & { strike: number; optionType: "call" | "put" } ): HestonPriceResult { const { strike, optionType, rate, timeToMaturityYears: T } = params; const simResult = simulateHestonPaths(params); const { terminalPrices } = simResult; // Calculate payoffs const payoffs = terminalPrices.map((S) => optionType === "call" ? Math.max(S - strike, 0) : Math.max(strike - S, 0) ); // Discount factor const df = Math.exp(-rate * T); // Discounted payoffs const discountedPayoffs = payoffs.map((p) => p * df); // Price and standard error const price = mean(discountedPayoffs); const se = stdDev(discountedPayoffs) / Math.sqrt(params.paths); return { price, standardError: se, confidenceInterval: { lower: price - 1.96 * se, upper: price + 1.96 * se, }, }; } /** * Calibrate Heston model to an IV surface * * Uses a simplified grid search + coordinate descent approach. * For production, consider using Levenberg-Marquardt or differential evolution. */ export function calibrateHestonSurface( input: HestonCalibrationInput ): HestonCalibrationResult { const { forwards, maturitiesYears, strikes, marketIvs, initialGuess } = input; // Validate inputs if (forwards.length !== maturitiesYears.length) { throw new Error("forwards and maturitiesYears must have same length"); } if (marketIvs.length !== maturitiesYears.length) { throw new Error("marketIvs must have one row per maturity"); } // Normalize strikes to 2D if needed const strikes2D: number[][] = Array.isArray(strikes[0]) ? (strikes as number[][]) : maturitiesYears.map(() => strikes as number[]); // Get ATM vol for initial guess const atmVols: number[] = []; for (let t = 0; t < maturitiesYears.length; t++) { const atmIdx = findClosestIndex(strikes2D[t], forwards[t]); atmVols.push(marketIvs[t][atmIdx]); } const avgAtmVol = mean(atmVols); // Initial parameters let kappa = initialGuess?.kappa ?? 2.0; let theta = initialGuess?.theta ?? avgAtmVol * avgAtmVol; let xi = initialGuess?.xi ?? 0.5; let rho = initialGuess?.rho ?? -0.7; let v0 = initialGuess?.v0 ?? avgAtmVol * avgAtmVol; // Parameter bounds const bounds = { kappa: { min: 0.1, max: 10 }, theta: { min: 0.001, max: 1 }, xi: { min: 0.01, max: 2 }, rho: { min: -0.99, max: 0.99 }, v0: { min: 0.001, max: 1 }, }; // Feller condition check: 2*kappa*theta > xi^2 const checkFeller = (k: number, t: number, x: number): boolean => { return 2 * k * t > x * x; }; // Objective: sum of squared IV errors (simplified - uses limited MC) const objective = ( k: number, th: number, x: number, r: number, vz: number ): number => { // For speed, use a small number of paths during calibration const calibPaths = 5000; const calibSteps = 20; let totalError = 0; let count = 0; for (let t = 0; t < maturitiesYears.length; t++) { const T = maturitiesYears[t]; const F = forwards[t]; // Sample a few strikes per maturity (for speed) const strikeIndices = [0, Math.floor(strikes2D[t].length / 2), strikes2D[t].length - 1]; for (const si of strikeIndices) { if (si >= strikes2D[t].length || si >= marketIvs[t].length) continue; const K = strikes2D[t][si]; const marketIv = marketIvs[t][si]; try { // Price option with Heston const result = priceOptionHestonMonteCarlo({ spot: F, strike: K, rate: 0, // Already forward optionType: K > F ? "call" : "put", kappa: k, theta: th, xi: x, rho: r, v0: vz, timeToMaturityYears: T, steps: calibSteps, paths: calibPaths, }); // Convert price back to IV (approximate) const modelIv = approxIvFromPrice(result.price, F, K, T, K > F ? "call" : "put"); if (!isNaN(modelIv) && isFinite(modelIv)) { const diff = modelIv - marketIv; totalError += diff * diff; count++; } } catch { // Skip failed calibration points } } } return count > 0 ? totalError / count : 1e10; }; // Coordinate descent const maxIter = 30; const tol = 1e-6; let prevError = objective(kappa, theta, xi, rho, v0); for (let iter = 0; iter < maxIter; iter++) { // Optimize each parameter in turn kappa = gridSearch((k) => objective(k, theta, xi, rho, v0), bounds.kappa.min, bounds.kappa.max, 10); theta = gridSearch((t) => objective(kappa, t, xi, rho, v0), bounds.theta.min, bounds.theta.max, 10); xi = gridSearch((x) => objective(kappa, theta, x, rho, v0), bounds.xi.min, bounds.xi.max, 10); rho = gridSearch((r) => objective(kappa, theta, xi, r, v0), bounds.rho.min, bounds.rho.max, 10); v0 = gridSearch((v) => objective(kappa, theta, xi, rho, v), bounds.v0.min, bounds.v0.max, 10); const error = objective(kappa, theta, xi, rho, v0); if (Math.abs(prevError - error) < tol) { return { kappa, theta, xi, rho, v0, error: Math.sqrt(error), converged: true }; } prevError = error; } return { kappa, theta, xi, rho, v0, error: Math.sqrt(prevError), converged: false, }; } /** * Simple grid search for 1D optimization */ function gridSearch( f: (x: number) => number, min: number, max: number, points: number ): number { let bestX = min; let bestF = f(min); for (let i = 1; i <= points; i++) { const x = min + (i / points) * (max - min); const fx = f(x); if (fx < bestF) { bestF = fx; bestX = x; } } return bestX; } /** * Approximate IV from option price using Brenner-Subrahmanyam approximation */ function approxIvFromPrice( price: number, F: number, K: number, T: number, optionType: "call" | "put" ): number { if (price <= 0 || T <= 0) return NaN; // For ATM, sigma ≈ price / (0.4 * F * sqrt(T)) const sqrtT = Math.sqrt(T); // Brenner-Subrahmanyam approximation (works best near ATM) const approxIv = price / (0.4 * F * sqrtT); // Clamp to reasonable range return Math.max(0.01, Math.min(3, approxIv)); } /** * Find index of element closest to target */ function findClosestIndex(arr: number[], target: number): number { let minDist = Infinity; let minIdx = 0; for (let i = 0; i < arr.length; i++) { const dist = Math.abs(arr[i] - target); if (dist < minDist) { minDist = dist; minIdx = i; } } return minIdx; } /** * Standard normal random number (Box-Muller) */ function randomStdNormal(): number { let u1: number, u2: number; do { u1 = Math.random(); u2 = Math.random(); } while (u1 === 0); return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2); } /** * Get Heston-implied distribution statistics */ export function hestonDistributionStats( params: Omit<HestonParams, "strike"> ): { mean: number; median: number; std: number; percentile5: number; percentile25: number; percentile75: number; percentile95: number; } { const simResult = simulateHestonPaths(params as HestonParams); const prices = simResult.terminalPrices; const sorted = [...prices].sort((a, b) => a - b); return { mean: mean(prices), median: percentile(sorted, 0.5), std: stdDev(prices), percentile5: percentile(sorted, 0.05), percentile25: percentile(sorted, 0.25), percentile75: percentile(sorted, 0.75), percentile95: percentile(sorted, 0.95), }; }