Quant Companion MCP

/** * Local volatility Monte Carlo simulation * * Uses a simplified Dupire-style local vol approach where volatility * varies with both spot level and time based on the implied vol surface. */ import type { VolSurface } from "./types"; import { randomNormal, mean, stdDev, percentile } from "./utils"; export interface LocalVolSimParams { /** Current spot price */ spot: number; /** Risk-free rate (annualized) */ rate: number; /** Dividend yield (annualized) */ dividendYield?: number; /** Time horizon in years */ timeHorizonYears: number; /** Vol surface for local vol lookup */ surface: VolSurface; /** Number of simulation paths (max 200000) */ paths?: number; /** Number of time steps */ steps?: number; /** Optional target price for probability calculation */ targetPrice?: number; } export interface LocalVolSimResult { /** Terminal price distribution statistics */ distribution: { mean: number; median: number; stdDev: number; percentile5: number; percentile25: number; percentile75: number; percentile95: number; min: number; max: number; }; /** Probability of reaching target (if specified) */ targetProbability?: number; /** VaR metrics */ riskMetrics: { var95: number; var99: number; expectedShortfall95: number; }; /** Sample of terminal prices (max 1000) */ samplePrices: number[]; /** Average local vol used across simulation */ avgLocalVol: number; /** Simulation parameters used */ params: { paths: number; steps: number; timeHorizonYears: number; }; } const MAX_PATHS = 200_000; const DEFAULT_PATHS = 50_000; const DEFAULT_STEPS = 50; /** * Simulate price paths using local volatility from the vol surface */ export function simulateWithLocalVol(params: LocalVolSimParams): LocalVolSimResult { const { spot, rate, dividendYield = 0, timeHorizonYears, surface, paths = DEFAULT_PATHS, steps = DEFAULT_STEPS, targetPrice, } = params; if (paths > MAX_PATHS) { throw new Error(`paths must be <= ${MAX_PATHS}`); } const dt = timeHorizonYears / steps; const sqrtDt = Math.sqrt(dt); const baseDrift = rate - dividendYield; const terminalPrices: number[] = []; let totalLocalVol = 0; let volSamples = 0; for (let i = 0; i < paths; i++) { let S = spot; let t = 0; for (let j = 0; j < steps; j++) { // Get local vol from surface at current (S, t) const localVol = getLocalVol(surface, S, t); totalLocalVol += localVol; volSamples++; // GBM Euler-Maruyama step with local vol // dS = S * ((r - q - 0.5*σ²)*dt + σ*dW) for risk-neutral drift const drift = baseDrift - 0.5 * localVol * localVol; const dW = randomNormal() * sqrtDt; const dS = S * (drift * dt + localVol * dW); S = Math.max(S + dS, 0.01); // Floor to prevent negative prices t += dt; } terminalPrices.push(S); } // Sort for percentile calculations const sorted = [...terminalPrices].sort((a, b) => a - b); // Calculate returns for VaR const returns = terminalPrices.map((p) => (p - spot) / spot); const sortedReturns = [...returns].sort((a, b) => a - b); // VaR (as loss, so negative of left tail) const var95 = -percentile(sortedReturns, 0.05); const var99 = -percentile(sortedReturns, 0.01); // Expected Shortfall (average of worst 5%) const worstN = Math.floor(paths * 0.05); const worstReturns = sortedReturns.slice(0, worstN); const expectedShortfall95 = -mean(worstReturns); // Target probability let targetProbability: number | undefined; if (targetPrice !== undefined) { const countAbove = terminalPrices.filter((p) => p >= targetPrice).length; targetProbability = countAbove / paths; } // Sample prices (max 1000) const sampleSize = Math.min(1000, paths); const samplePrices = terminalPrices.slice(0, sampleSize); return { distribution: { mean: mean(terminalPrices), median: percentile(sorted, 0.5), stdDev: stdDev(terminalPrices), percentile5: percentile(sorted, 0.05), percentile25: percentile(sorted, 0.25), percentile75: percentile(sorted, 0.75), percentile95: percentile(sorted, 0.95), min: sorted[0], max: sorted[sorted.length - 1], }, targetProbability, riskMetrics: { var95, var99, expectedShortfall95, }, samplePrices, avgLocalVol: totalLocalVol / volSamples, params: { paths, steps, timeHorizonYears, }, }; } /** Minimum allowed local vol (1% annualized) */ const MIN_LOCAL_VOL = 0.01; /** Maximum allowed local vol (300% annualized) */ const MAX_LOCAL_VOL = 3.0; /** * Clamp volatility to reasonable bounds */ function clampVol(vol: number): number { return Math.max(MIN_LOCAL_VOL, Math.min(MAX_LOCAL_VOL, vol)); } /** * Get local volatility at a specific (spot, time) point * * Uses bilinear interpolation on the vol surface. * Falls back to ATM vol if point is outside surface bounds. * Clamps result to [1%, 300%] to prevent numerical instability. */ function getLocalVol(surface: VolSurface, S: number, t: number): number { const { smiles, spot, stats } = surface; if (smiles.length === 0) { return stats.avgAtmIv; } // Sort smiles by maturity const sortedSmiles = [...smiles].sort( (a, b) => a.timeToMaturityYears - b.timeToMaturityYears ); // Find bracketing maturities let smileBelow = sortedSmiles[0]; let smileAbove = sortedSmiles[sortedSmiles.length - 1]; for (const smile of sortedSmiles) { if (smile.timeToMaturityYears <= t) { smileBelow = smile; } if (smile.timeToMaturityYears >= t) { smileAbove = smile; break; } } // Get IV at current spot level from each smile const moneyness = S / spot; const ivBelow = interpolateSmileAtMoneyness(smileBelow, moneyness); const ivAbove = interpolateSmileAtMoneyness(smileAbove, moneyness); // Interpolate in time if (smileBelow === smileAbove || smileBelow.timeToMaturityYears === smileAbove.timeToMaturityYears) { return ivBelow; } const w = (t - smileBelow.timeToMaturityYears) / (smileAbove.timeToMaturityYears - smileBelow.timeToMaturityYears); // Clamp weight to [0, 1] for extrapolation cases const clampedW = Math.max(0, Math.min(1, w)); const interpolatedVol = ivBelow + clampedW * (ivAbove - ivBelow); return clampVol(interpolatedVol); } /** * Interpolate IV at a specific moneyness on a smile */ function interpolateSmileAtMoneyness( smile: { points: Array<{ moneyness: number; iv: number }>; atmIv: number }, targetMoneyness: number ): number { const { points, atmIv } = smile; if (points.length === 0) { return atmIv; } // Find bracketing points let below: typeof points[0] | null = null; let above: typeof points[0] | null = null; for (const p of points) { if (p.moneyness <= targetMoneyness) { if (!below || p.moneyness > below.moneyness) below = p; } if (p.moneyness >= targetMoneyness) { if (!above || p.moneyness < above.moneyness) above = p; } } if (below && above && below !== above) { const w = (targetMoneyness - below.moneyness) / (above.moneyness - below.moneyness); return below.iv + w * (above.iv - below.iv); } // Extrapolation: use closest point with some dampening toward ATM if (below) return below.iv; if (above) return above.iv; return atmIv; }