fx-risk-terminal/risk_engine.py

"""
FX risk engine — Value-at-Risk, Expected Shortfall, Monte Carlo, stress
testing and minimum-variance hedging. Generalised to N currencies.

Everything here is pure Python (stdlib only): the covariance pipeline, the
Cholesky decomposition used to correlate Monte-Carlo shocks, and the three
independent VaR methodologies. No numpy, no pandas — the maths is on show.

Conventions
-----------
* The currency order is fixed by the caller (a sorted list) and every vector /
  matrix is built in that order.
* "Exposure" is the signed USD value of each currency sleeve. A +1% move in
  CCY/USD changes that sleeve's USD value by exposure * 0.01.
* Returns are daily log returns. VaR / ES are positive loss numbers.
* Horizon defaults to 1 trading day; multi-day scales by sqrt(h).
"""

import math
import random
from typing import Sequence

# Standard-normal quantiles for the confidence levels we report.
Z = {0.95: 1.6448536269514722, 0.99: 2.3263478740408408}


# ── Linear algebra ────────────────────────────────────────────────────────────

def _mean(xs: Sequence[float]) -> float:
    return sum(xs) / len(xs) if xs else 0.0


def covariance_matrix(returns: dict[str, list[float]], currencies: list[str]) -> list[list[float]]:
    series = [returns[c] for c in currencies]
    if not series:
        return []
    n = min(len(s) for s in series)
    series = [s[-n:] for s in series]
    means = [_mean(s) for s in series]
    k = len(series)
    cov = [[0.0] * k for _ in range(k)]
    if n < 2:
        return cov
    for i in range(k):
        for j in range(k):
            cov[i][j] = sum((series[i][t] - means[i]) * (series[j][t] - means[j])
                            for t in range(n)) / (n - 1)
    return cov


def correlation_matrix(cov: list[list[float]]) -> list[list[float]]:
    k = len(cov)
    corr = [[0.0] * k for _ in range(k)]
    for i in range(k):
        for j in range(k):
            d = math.sqrt(cov[i][i] * cov[j][j])
            corr[i][j] = cov[i][j] / d if d > 0 else 0.0
    return corr


def cholesky(matrix: list[list[float]]) -> list[list[float]]:
    """Lower-triangular L with L·Lᵀ = matrix. Jitters non-PSD diagonals."""
    n = len(matrix)
    L = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = sum(L[i][k] * L[j][k] for k in range(j))
            if i == j:
                diag = matrix[i][i] - s
                L[i][j] = math.sqrt(diag if diag > 0 else 1e-18)
            else:
                L[i][j] = (matrix[i][j] - s) / L[j][j] if L[j][j] else 0.0
    return L


def _matvec(M, v):
    return [sum(M[i][k] * v[k] for k in range(len(v))) for i in range(len(M))]


def _percentile(sorted_xs: list[float], q: float) -> float:
    if not sorted_xs:
        return 0.0
    if len(sorted_xs) == 1:
        return sorted_xs[0]
    pos = q * (len(sorted_xs) - 1)
    lo, hi = math.floor(pos), math.ceil(pos)
    if lo == hi:
        return sorted_xs[int(pos)]
    frac = pos - lo
    return sorted_xs[lo] * (1 - frac) + sorted_xs[hi] * frac


# ── Exposure ──────────────────────────────────────────────────────────────────

def exposure_vector(portfolio: dict, currencies: list[str]) -> list[float]:
    exp = portfolio.get("exposure_usd", {})
    return [float(exp.get(c, 0.0)) for c in currencies]


# ── Parametric VaR + component VaR ────────────────────────────────────────────

def parametric_var(exposure, cov, currencies, horizon=1) -> dict:
    sig = _matvec(cov, exposure)
    var_usd = sum(exposure[i] * sig[i] for i in range(len(exposure)))
    sigma = math.sqrt(max(var_usd, 0.0)) * math.sqrt(horizon)

    out = {"sigma_usd": round(sigma, 4), "levels": {}, "components": {}}
    for cl, z in Z.items():
        out["levels"][str(cl)] = round(z * sigma, 4)

    z95 = Z[0.95]
    if sigma > 0 and var_usd > 0:
        denom = math.sqrt(var_usd)
        for i, c in enumerate(currencies):
            out["components"][c] = round(exposure[i] * sig[i] / denom * z95 * math.sqrt(horizon), 4)
    else:
        for c in currencies:
            out["components"][c] = 0.0
    return out


# ── Historical-simulation VaR ─────────────────────────────────────────────────

def historical_var(exposure, returns, currencies, horizon=1) -> dict:
    series = [returns[c] for c in currencies]
    n = min((len(s) for s in series), default=0)
    series = [s[-n:] for s in series]

    pnl = []
    for t in range(n):
        p = sum(exposure[i] * (math.exp(series[i][t]) - 1) for i in range(len(exposure)))
        pnl.append(p * math.sqrt(horizon))
    pnl.sort()

    out = {"levels": {}, "es": {}}
    for cl in Z:
        vq = -_percentile(pnl, 1 - cl)
        tail = [x for x in pnl if x <= -vq]
        out["levels"][str(cl)] = round(vq, 4)
        out["es"][str(cl)] = round(-_mean(tail) if tail else vq, 4)
    return out


# ── Monte Carlo VaR ───────────────────────────────────────────────────────────

def monte_carlo_var(exposure, cov, horizon=1, n_sims=50_000, seed=42, bins=41) -> dict:
    rng = random.Random(seed)
    L = cholesky(cov)
    k = len(exposure)
    h = math.sqrt(horizon)

    pnl = []
    for _ in range(n_sims):
        z = [rng.gauss(0, 1) for _ in range(k)]
        shock = _matvec(L, z)
        pnl.append(sum(exposure[i] * (math.exp(shock[i] * h) - 1) for i in range(k)))
    pnl.sort()

    out = {"levels": {}, "es": {}, "n_sims": n_sims}
    for cl in Z:
        vq = -_percentile(pnl, 1 - cl)
        tail = [x for x in pnl if x <= -vq]
        out["levels"][str(cl)] = round(vq, 4)
        out["es"][str(cl)] = round(-_mean(tail) if tail else vq, 4)

    lo, hi = pnl[0], pnl[-1]
    width = (hi - lo) / bins if hi > lo else 1.0
    counts = [0] * bins
    for p in pnl:
        counts[min(int((p - lo) / width), bins - 1)] += 1
    out["histogram"] = {
        "centers": [round(lo + (i + 0.5) * width, 4) for i in range(bins)],
        "counts": counts,
        "var95": out["levels"]["0.95"], "var99": out["levels"]["0.99"],
        "mean": round(_mean(pnl), 4),
    }
    return out


# ── Hedging: best single-instrument min-variance hedge ────────────────────────

def best_hedge(exposure, cov, currencies) -> dict:
    """
    Find the single FX forward that, sold against the book, minimises residual
    variance — and report its effectiveness (R² = variance reduction).

    For instrument j the optimal short notional is
        H* = Cov(book, r_j) / Var(r_j) = (Σe)_j / Σ_jj
    and the variance reduction equals Corr(book, r_j)².
    """
    sig = _matvec(cov, exposure)                       # Σe
    unhedged_var = sum(exposure[i] * sig[i] for i in range(len(exposure)))

    best = None
    for j, c in enumerate(currencies):
        var_j = cov[j][j]
        if var_j <= 0:
            continue
        cov_bj = sig[j]                                # Cov(book, r_j)
        h_star = cov_bj / var_j
        residual_var = unhedged_var - cov_bj ** 2 / var_j
        eff = (1 - residual_var / unhedged_var) if unhedged_var > 0 else 0.0
        cand = {
            "instrument": c,
            "notional_short_usd": round(h_star, 2),
            "effectiveness_pct": round(eff * 100, 2),
            "residual_sigma": round(math.sqrt(max(residual_var, 0)), 4),
        }
        if best is None or cand["effectiveness_pct"] > best["effectiveness_pct"]:
            best = cand

    # Direct full hedge: short each sleeve in its own pair → residual ≈ 0.
    direct = {c: round(exposure[i], 2) for i, c in enumerate(currencies)}

    return {
        "unhedged_sigma": round(math.sqrt(max(unhedged_var, 0)), 4),
        "best_single": best or {},
        "direct_hedge": direct,
    }


# ── Stress testing ────────────────────────────────────────────────────────────

# Stylised one-shot FX shocks (% move in each CCY/USD) drawn from notable
# historical regimes. Note the safe-haven behaviour of JPY/CHF (they rally as
# the dollar bid hits EM/high-beta currencies) — that is the educational point.
STRESS_SCENARIOS = [
    {"name": "GFC — Oct 2008", "note": "Lehman aftermath, USD funding squeeze",
     "shocks": {"EUR": -8.5, "SEK": -13, "GBP": -10, "CHF": -2, "BRL": -20, "MXN": -15, "JPY": 12, "KRW": -25}},
    {"name": "Euro debt crisis 2011", "note": "Peripheral spreads blow out; SNB caps CHF",
     "shocks": {"EUR": -4.5, "SEK": -3, "GBP": -2, "CHF": 8, "BRL": -10, "MXN": -8, "JPY": 5, "KRW": -4}},
    {"name": "USD bull 2014–15", "note": "Fed exit vs. global easing; broad dollar rally",
     "shocks": {"EUR": -10, "SEK": -10, "GBP": -5, "CHF": -8, "BRL": -15, "MXN": -12, "JPY": -12, "KRW": -7}},
    {"name": "COVID crash — Mar 2020", "note": "Dash for dollars; EM liquidity drain",
     "shocks": {"EUR": -3, "SEK": -7.5, "GBP": -10, "CHF": -2, "BRL": -15, "MXN": -20, "JPY": 2, "KRW": -6}},
    {"name": "USD surge — Sep 2022", "note": "Fed hiking into energy shock; gilt crisis; super-peso",
     "shocks": {"EUR": -6.5, "SEK": -9, "GBP": -8, "CHF": -3, "BRL": 2, "MXN": -1, "JPY": -10, "KRW": -8}},
    {"name": "Risk-on rally", "note": "Soft landing / broad USD weakness",
     "shocks": {"EUR": 4, "SEK": 6, "GBP": 4, "CHF": -1, "BRL": 8, "MXN": 6, "JPY": -3, "KRW": 5}},
]


def stress_test(exposure, currencies) -> list[dict]:
    e = dict(zip(currencies, exposure))
    total = sum(e.values())
    out = []
    for sc in STRESS_SCENARIOS:
        impact = sum(e[c] * (sc["shocks"].get(c, 0) / 100) for c in currencies)
        out.append({
            "name": sc["name"], "note": sc["note"],
            "shocks": {c: sc["shocks"].get(c, 0) for c in currencies},
            "impact_usd": round(impact, 2),
            "impact_pct": round(impact / total * 100, 3) if total else None,
        })
    return out


# ── Orchestration ─────────────────────────────────────────────────────────────

def full_risk_report(portfolio: dict, returns: dict[str, list[float]], horizon: int = 1) -> dict:
    currencies = portfolio.get("currencies") or sorted(portfolio.get("exposure_usd", {}))
    exposure = exposure_vector(portfolio, currencies)
    cov = covariance_matrix(returns, currencies)
    corr = correlation_matrix(cov) if cov else []

    total = sum(exposure)
    param = parametric_var(exposure, cov, currencies, horizon)
    hist  = historical_var(exposure, returns, currencies, horizon)
    mc    = monte_carlo_var(exposure, cov, horizon)
    hedge = best_hedge(exposure, cov, currencies)

    port_vol_ann = param["sigma_usd"] * math.sqrt(252) / total * 100 if total else None

    # Diversification ratio: weighted-avg sleeve vol / portfolio vol. >1 means
    # correlations are doing work; ~1 means the book moves as one bet.
    div_ratio = None
    if total and cov:
        w_sigma = sum(abs(exposure[i]) * math.sqrt(cov[i][i]) for i in range(len(exposure)))
        p_sigma = param["sigma_usd"]
        div_ratio = round(w_sigma / p_sigma, 3) if p_sigma > 0 else None

    def pct(v): return round(v / total * 100, 3) if total else None

    return {
        "horizon_days": horizon,
        "currencies": currencies,
        "total_exposure_usd": round(total, 2),
        "portfolio_vol_annual_pct": round(port_vol_ann, 2) if port_vol_ann else None,
        "diversification_ratio": div_ratio,
        "covariance": cov,
        "correlation": corr,
        "var": {
            "parametric": {
                "var95": param["levels"]["0.95"], "var99": param["levels"]["0.99"],
                "var95_pct": pct(param["levels"]["0.95"]), "var99_pct": pct(param["levels"]["0.99"]),
            },
            "historical": {
                "var95": hist["levels"]["0.95"], "var99": hist["levels"]["0.99"],
                "es95": hist["es"]["0.95"], "es99": hist["es"]["0.99"],
                "var95_pct": pct(hist["levels"]["0.95"]), "var99_pct": pct(hist["levels"]["0.99"]),
            },
            "monte_carlo": {
                "var95": mc["levels"]["0.95"], "var99": mc["levels"]["0.99"],
                "es95": mc["es"]["0.95"], "es99": mc["es"]["0.99"],
                "var95_pct": pct(mc["levels"]["0.95"]), "var99_pct": pct(mc["levels"]["0.99"]),
                "n_sims": mc["n_sims"],
            },
        },
        "component_var95": param["components"],
        "es95": mc["es"]["0.95"],
        "distribution": mc["histogram"],
        "hedge": hedge,
    }