Quant Companion MCP

/** * SABR Stochastic Volatility Model * * Implements the SABR implied volatility approximation (Hagan et al. 2002) * and smile calibration. */ import { normalCDF } from "./utils"; export interface SabrParams { /** Forward price */ forward: number; /** Strike price */ strike: number; /** Time to maturity in years */ timeToMaturityYears: number; /** Alpha (ATM vol level) */ alpha: number; /** Beta (CEV exponent, typically 0 to 1) */ beta: number; /** Rho (correlation between spot and vol) */ rho: number; /** Nu (vol of vol) */ nu: number; } export interface SabrCalibrationInput { /** Forward price */ forward: number; /** Time to maturity in years */ timeToMaturityYears: number; /** Array of strikes */ strikes: number[]; /** Array of market implied volatilities (same length as strikes) */ marketIvs: number[]; /** Fixed beta (optional, default 1.0 for lognormal) */ beta?: number; /** Initial guess for alpha */ alphaGuess?: number; /** Initial guess for rho */ rhoGuess?: number; /** Initial guess for nu */ nuGuess?: number; } export interface SabrCalibrationResult { /** Calibrated alpha */ alpha: number; /** Beta used (fixed or calibrated) */ beta: number; /** Calibrated rho */ rho: number; /** Calibrated nu */ nu: number; /** Root mean square error of fit */ error: number; /** Model IVs at input strikes */ modelIvs: number[]; } /** * Compute SABR implied volatility using Hagan's approximation * * Reference: Hagan et al. (2002) "Managing Smile Risk" * * @param params SABR parameters * @returns Black implied volatility */ export function sabrImpliedVol(params: SabrParams): number { const { forward: F, strike: K, timeToMaturityYears: T, alpha, beta, rho, nu } = params; // Handle ATM case if (Math.abs(F - K) < 1e-10) { return sabrAtmVol(F, T, alpha, beta, rho, nu); } // Handle zero time if (T <= 0) { return alpha; } // Validate parameters if (alpha <= 0 || nu < 0 || beta < 0 || beta > 1 || rho < -1 || rho > 1) { return NaN; } const logFK = Math.log(F / K); const FK_mid = Math.sqrt(F * K); const FK_beta = Math.pow(FK_mid, 1 - beta); // z parameter const z = (nu / alpha) * FK_beta * logFK; // x(z) function let xz: number; if (Math.abs(z) < 1e-10) { xz = 1; } else { const sqrtTerm = Math.sqrt(1 - 2 * rho * z + z * z); const numerator = sqrtTerm + z - rho; const denominator = 1 - rho; if (numerator <= 0 || denominator <= 0) { xz = 1; } else { xz = z / Math.log(numerator / denominator); } } // Correction terms const oneBeta = 1 - beta; const logFK2 = logFK * logFK; const FK_2beta = Math.pow(FK_mid, 2 * oneBeta); // Denominator term const denom1 = 1 + (oneBeta * oneBeta / 24) * logFK2 + (Math.pow(oneBeta, 4) / 1920) * logFK2 * logFK2; // Numerator correction const term1 = (oneBeta * oneBeta / 24) * (alpha * alpha / FK_2beta); const term2 = 0.25 * rho * beta * nu * alpha / FK_beta; const term3 = (2 - 3 * rho * rho) * nu * nu / 24; const correction = 1 + (term1 + term2 + term3) * T; // Final SABR vol const sigmaB = (alpha / (FK_beta * denom1)) * xz * correction; return Math.max(sigmaB, 0.001); // Floor at 0.1% } /** * SABR ATM volatility approximation */ function sabrAtmVol( F: number, T: number, alpha: number, beta: number, rho: number, nu: number ): number { const oneBeta = 1 - beta; const F_beta = Math.pow(F, oneBeta); const term1 = (oneBeta * oneBeta / 24) * (alpha * alpha / (F_beta * F_beta)); const term2 = 0.25 * rho * beta * nu * alpha / F_beta; const term3 = (2 - 3 * rho * rho) * nu * nu / 24; const correction = 1 + (term1 + term2 + term3) * T; return (alpha / F_beta) * correction; } /** * Build a SABR smile (IV curve) for a range of strikes * * @param params Base SABR parameters (strike will be varied) * @param strikes Array of strikes to compute IV for * @returns Array of implied volatilities */ export function buildSabrSmile( params: Omit<SabrParams, "strike">, strikes: number[] ): number[] { return strikes.map(strike => sabrImpliedVol({ ...params, strike })); } /** * Calibrate SABR parameters to market smile * * Uses Levenberg-Marquardt style optimization with simple gradient descent fallback. * Beta is typically fixed (common choices: 0, 0.5, or 1). * * @param input Calibration input data * @returns Calibrated parameters and fit error */ export function calibrateSabrSmile(input: SabrCalibrationInput): SabrCalibrationResult { const { forward, timeToMaturityYears, strikes, marketIvs, beta = 1.0, alphaGuess, rhoGuess, nuGuess, } = input; if (strikes.length !== marketIvs.length) { throw new Error("strikes and marketIvs must have same length"); } if (strikes.length < 3) { throw new Error("Need at least 3 points for SABR calibration"); } // Find ATM IV for initial alpha guess const atmIdx = findClosestIndex(strikes, forward); const atmIv = marketIvs[atmIdx]; // Initial parameter guesses let alpha = alphaGuess ?? atmIv * Math.pow(forward, 1 - beta); let rho = rhoGuess ?? -0.3; let nu = nuGuess ?? 0.4; // Bounds const bounds = { alpha: { min: 0.001, max: 5.0 }, rho: { min: -0.999, max: 0.999 }, nu: { min: 0.001, max: 3.0 }, }; // Objective function: sum of squared errors const objective = (a: number, r: number, n: number): number => { let sse = 0; for (let i = 0; i < strikes.length; i++) { const modelIv = sabrImpliedVol({ forward, strike: strikes[i], timeToMaturityYears, alpha: a, beta, rho: r, nu: n, }); if (isNaN(modelIv)) return 1e10; const diff = modelIv - marketIvs[i]; sse += diff * diff; } return sse; }; // Simple coordinate descent optimization const maxIter = 200; const tol = 1e-8; let prevError = objective(alpha, rho, nu); for (let iter = 0; iter < maxIter; iter++) { // Optimize alpha alpha = goldenSectionSearch( (a) => objective(a, rho, nu), bounds.alpha.min, bounds.alpha.max, tol ); // Optimize rho rho = goldenSectionSearch( (r) => objective(alpha, r, nu), bounds.rho.min, bounds.rho.max, tol ); // Optimize nu nu = goldenSectionSearch( (n) => objective(alpha, rho, n), bounds.nu.min, bounds.nu.max, tol ); const error = objective(alpha, rho, nu); if (Math.abs(prevError - error) < tol) { break; } prevError = error; } // Compute final model IVs and RMSE const modelIvs = strikes.map(strike => sabrImpliedVol({ forward, strike, timeToMaturityYears, alpha, beta, rho, nu, }) ); const rmse = Math.sqrt( modelIvs.reduce((sum, iv, i) => sum + Math.pow(iv - marketIvs[i], 2), 0) / strikes.length ); return { alpha, beta, rho, nu, error: rmse, modelIvs, }; } /** * Golden section search for 1D minimization */ function goldenSectionSearch( f: (x: number) => number, a: number, b: number, tol: number ): number { const phi = (1 + Math.sqrt(5)) / 2; const resphi = 2 - phi; let x1 = a + resphi * (b - a); let x2 = b - resphi * (b - a); let f1 = f(x1); let f2 = f(x2); while (Math.abs(b - a) > tol) { if (f1 < f2) { b = x2; x2 = x1; f2 = f1; x1 = a + resphi * (b - a); f1 = f(x1); } else { a = x1; x1 = x2; f1 = f2; x2 = b - resphi * (b - a); f2 = f(x2); } } return (a + b) / 2; } /** * 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; } /** * Generate SABR-based terminal price distribution via sampling * * Uses Black-Scholes with SABR-implied vol at each simulated strike * to generate a distribution consistent with the SABR smile. * * @param params SABR parameters (strike not used, forward is spot) * @param paths Number of samples to generate * @returns Array of terminal prices */ export function sampleSabrDistribution( params: Omit<SabrParams, "strike">, rate: number, paths: number ): number[] { const { forward, timeToMaturityYears: T, alpha, beta, rho, nu } = params; // For SABR, we sample by: // 1. Generate uniform random percentiles // 2. Map to strikes via inverse CDF approximation // 3. Use those as terminal prices // Simpler approach: MC with local vol from SABR smile const sqrtT = Math.sqrt(T); const dt = T / 50; // 50 steps const sqrtDt = Math.sqrt(dt); const prices: number[] = []; for (let i = 0; i < paths; i++) { let S = forward; let t = 0; for (let j = 0; j < 50; j++) { // Get local vol from SABR at current spot const localVol = sabrImpliedVol({ forward, strike: S, timeToMaturityYears: Math.max(T - t, 0.001), alpha, beta, rho, nu, }); // GBM step with SABR local vol const z = randomStdNormal(); const drift = rate - 0.5 * localVol * localVol; S = S * Math.exp(drift * dt + localVol * sqrtDt * z); S = Math.max(S, 0.01); t += dt; } prices.push(S); } return prices; } /** * 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); }