Interest Rate Curve Builder

# curve_builder.py — Multi-curve bootstrapping for EUR OIS and EURIBOR
from __future__ import annotations

from dataclasses import dataclass, field
from enum import Enum, auto
from typing import Optional

import numpy as np
from scipy.interpolate import CubicSpline


class InstrumentType(Enum):
    DEPOSIT  = auto()   # Money market deposit: ACT/360, simple compounding
    OIS_SWAP = auto()   # OIS: fixed vs daily-compounded ESTR (ACT/360 float leg)
    FRA      = auto()   # Forward Rate Agreement: 3M EURIBOR, ACT/360
    IRS      = auto()   # Interest Rate Swap: fixed (ann, 30/360) vs 3M EURIBOR


@dataclass(frozen=True)
class Instrument:
    """
    Single market instrument input for yield curve bootstrapping.

    Conventions:
        rate     — decimal (0.035 = 3.5%)
        maturity — end maturity in years from today (ACT/365 approximation)
        start    — start maturity in years (used for FRAs; 0 for all others)
    """
    itype:    InstrumentType
    rate:     float
    maturity: float
    start:    float = 0.0
    label:    str   = ""


def _act360(t: float) -> float:
    """Convert a maturity in years (ACT/365 basis) to an ACT/360 day count fraction."""
    return t * 365.0 / 360.0


@dataclass
class DiscountCurve:
    """
    Piecewise-cubic interpolated discount curve built from (t_i, P_i) pillars.

    Interpolation: natural cubic spline on ln P(0,t).

    Rationale for log-space interpolation:
        The forward rate between t₁ and t₂ is:
            f(t₁,t₂) = −[ln P(t₂) − ln P(t₁)] / (t₂ − t₁)
        Interpolating ln P directly (rather than P or zero rates) ensures that
        the forward rate is controlled by the slope of a smooth cubic — it stays
        finite and bounded even between sparse pillars.

    Extrapolation: flat instantaneous forward rate beyond the last pillar.
        This is the industry-standard choice: it avoids the explosive growth
        that polynomial extrapolation would produce.

    Assumptions:
        - P(0,0) = 1 by definition.
        - All rates: continuously compounded, annualised.
        - No day-count adjustment beyond the ACT/360 conversion in the bootstrapper.
    """
    _pillars: list[tuple[float, float]] = field(default_factory=list)
    _spline: Optional[CubicSpline]      = field(default=None, init=False, repr=False)

    def add_pillar(self, t: float, P: float) -> None:
        """Append a (t, P(0,t)) pillar and rebuild the cubic spline."""
        if not (0 < P <= 1.0 + 1e-9):
            raise ValueError(
                f"Discount factor must satisfy 0 < P ≤ 1, got P({t:.4f}) = {P:.8f}. "
                "Check market data for implied negative forward rates."
            )
        self._pillars.append((t, max(P, 1e-12)))
        self._pillars.sort(key=lambda x: x[0])
        self._rebuild()

    def _rebuild(self) -> None:
        if len(self._pillars) < 2:
            return
        ts  = np.array([0.0] + [p[0] for p in self._pillars])
        Ps  = np.array([1.0] + [p[1] for p in self._pillars])
        self._spline = CubicSpline(ts, np.log(Ps), bc_type="not-a-knot")

    def df(self, t: float) -> float:
        """Discount factor P(0,t)."""
        if t <= 0.0:
            return 1.0
        if self._spline is None:
            raise RuntimeError("Curve has fewer than 2 pillars")
        t_max = self._pillars[-1][0]
        if t <= t_max:
            return float(np.exp(self._spline(t)))
        # Flat forward extrapolation
        f_last   = float(-self._spline(t_max, 1))   # −d(ln P)/dt at last pillar
        lnP_last = float(self._spline(t_max))
        return float(np.exp(lnP_last - f_last * (t - t_max)))

    def zero_rate(self, t: float) -> float:
        """Continuously compounded zero-coupon rate r(t) = −ln P(0,t) / t."""
        if t < 1e-9:
            return float(-self._spline(1e-6, 1)) if self._spline else 0.0
        return -np.log(self.df(t)) / t

    def inst_forward(self, t: float) -> float:
        """Instantaneous forward rate f(t) = −d(ln P)/dt."""
        if self._spline is None:
            raise RuntimeError("Curve not built")
        return float(-self._spline(max(t, 1e-9), 1))

    def forward_rate(self, t1: float, t2: float) -> float:
        """Simply-compounded forward rate F(0; t₁, t₂) = [P(t₁)/P(t₂) − 1]/(t₂−t₁)."""
        if t2 <= t1:
            raise ValueError(f"t2={t2} must exceed t1={t1}")
        return (self.df(t1) / self.df(t2) - 1.0) / (t2 - t1)

    @property
    def pillars(self) -> list[tuple[float, float]]:
        return list(self._pillars)

def bootstrap_ois(instruments: list[Instrument]) -> DiscountCurve:
    """
    Bootstrap the EUR OIS discount curve from ESTR-linked market instruments.

    Instrument types handled:
        DEPOSIT  — direct inversion: P(0,T) = 1 / (1 + r·δ(0,T))
        OIS_SWAP — iterative: solve for P(0,T_N) given all prior pillars

    OIS swap pricing equation (fixed annual vs daily-compounded ESTR):

        Fixed leg PV  =  K · Σᵢ δᵢ · P(0,Tᵢ)
        Float leg PV  =  P(0,T₀) − P(0,T_N)    [OIS ≈ compounding of overnight rates]

    Setting Fixed = Float and isolating the last unknown P(0,T_N):

        P(0,T_N) = [1 − K · Σᵢ₌₁^{N−1} δᵢ · P(0,Tᵢ)] / (1 + K · δ_N)

    This is the closed-form sequential bootstrap — no root-finding required.

    Limitations:
        - Fixed payment frequency assumed annual (standard EUR OIS convention).
        - Day count approximated as ACT/360 without exact business day adjustment.
          Use QuantLib for production schedules with holiday calendars (TARGET2).
        - For OIS tenors > 1Y with semi-annual fixed frequency, modify T_i grid.
    """
    curve = DiscountCurve()
    insts = sorted(instruments, key=lambda i: i.maturity)

    for inst in insts:
        T, K = inst.maturity, inst.rate

        if inst.itype == InstrumentType.DEPOSIT:
            # Direct inversion: P(0,T) = 1 / (1 + K·δ)
            P = 1.0 / (1.0 + K * _act360(T))
            curve.add_pillar(T, P)

        elif inst.itype == InstrumentType.OIS_SWAP:
            n = max(1, int(round(T)))
            T_i = [min(k, T) for k in range(1, n + 1)]
            δ_i = [_act360(T_i[0])] + [_act360(T_i[k] - T_i[k-1]) for k in range(1, n)]

            # Annuity from already-bootstrapped pillars (indices 0 … N-2)
            annuity = sum(δ_i[k] * curve.df(T_i[k]) for k in range(n - 1))

            P_N = (1.0 - K * annuity) / (1.0 + K * δ_i[-1])
            curve.add_pillar(T_i[-1], P_N)

    return curve

def bootstrap_euribor(
    instruments: list[Instrument],
    ois_curve: DiscountCurve,
) -> DiscountCurve:
    """
    Bootstrap the EURIBOR 3M projection curve using the OIS discount curve.

    Post-2008 multi-curve framework (Mercurio 2010):
        All cash flows are discounted on the OIS (risk-free) curve.
        EURIBOR forward rates are projected from a separate EURIBOR tenor curve.

    Instrument types:
        DEPOSIT — defines the 3M pillar directly: P_E(0,T) = 1/(1+K·δ)

        FRA — 3M EURIBOR forward, p×q months:
            F(0; T₁, T₂) = K  ⟹  P_E(0,T₁)/P_E(0,T₂) = 1 + K·δ(T₁,T₂)
            Given P_E(0,T₁): P_E(0,T₂) = P_E(0,T₁) / (1 + K·δ(T₁,T₂))

        IRS — par swap rate K (annual fixed vs quarterly 3M EURIBOR):

            Fixed PV  = K · Σᵢ δᵢ_fix · P_OIS(0,Tᵢ)
            Float PV  = Σⱼ [P_E(0,Tⱼ₋₁)/P_E(0,Tⱼ) − 1] · P_OIS(0,Tⱼ)

            Setting equal and solving for the last unknown P_E(0,T_N):

                Let  A = Σᵢ δᵢ_fix · P_OIS(0,Tᵢ)           [fixed-leg annuity]
                     B = Σⱼ₌₁^{N−1} [P_OIS(0,Tⱼ)·P_E(0,Tⱼ₋₁)/P_E(0,Tⱼ) − P_OIS(0,Tⱼ)]
                         [float PV from known pillars]

                P_E(0,T_N) = P_OIS(0,T_N) · P_E(0,T_{N−1}) / (K·A − B + P_OIS(0,T_N))

    Limitations:
        - OIS discounting assumes no CVA/DVA (collateralised trades, CSA assumed).
        - Tenor mismatch: if hedging with 6M instruments, need a separate 6M curve.
        - Basis between EURIBOR 3M and 6M is real (≈5–15bp pre-2022) and must be
          captured if pricing cross-tenor products (EURIBOR 3M vs 6M basis swaps).
    """
    curve = DiscountCurve()
    insts = sorted(instruments, key=lambda i: i.maturity)

    for inst in insts:
        T, T0, K = inst.maturity, inst.start, inst.rate

        if inst.itype == InstrumentType.DEPOSIT:
            P = 1.0 / (1.0 + K * _act360(T - T0))
            curve.add_pillar(T, P)

        elif inst.itype == InstrumentType.FRA:
            δ    = _act360(T - T0)
            P_T0 = curve.df(T0) if T0 > 1e-9 else 1.0
            curve.add_pillar(T, P_T0 / (1.0 + K * δ))

        elif inst.itype == InstrumentType.IRS:
            n_fix = max(1, int(round(T)))
            T_fix = [min(k, T) for k in range(1, n_fix + 1)]
            δ_fix = [_act360(T_fix[0])] + [_act360(T_fix[k] - T_fix[k-1]) for k in range(1, n_fix)]
            A = sum(d * ois_curve.df(t) for d, t in zip(δ_fix, T_fix))

            n_flt = max(1, int(round(T * 4)))
            T_flt = [min(k * 0.25, T) for k in range(1, n_flt + 1)]

            B, T_prev = 0.0, 0.0
            for j in range(len(T_flt) - 1):
                Tj     = T_flt[j]
                P_prev = curve.df(T_prev) if T_prev > 1e-9 else 1.0
                P_j    = curve.df(Tj)
                if P_j <= 0:
                    break
                B     += ois_curve.df(Tj) * (P_prev / P_j - 1.0)
                T_prev = Tj

            T_prev_N = T_flt[-2] if len(T_flt) > 1 else 0.0
            P_E_prev = curve.df(T_prev_N) if T_prev_N > 1e-9 else 1.0
            denom    = K * A - B + ois_curve.df(T)
            if abs(denom) < 1e-12:
                raise ValueError(f"Degenerate IRS bootstrap at T={T:.2f}Y (denom={denom:.2e})")
            curve.add_pillar(T, ois_curve.df(T) * P_E_prev / denom)

    return curve

# nelson_siegel.py — Nelson-Siegel (1987) and Svensson (1994) parametrisations
from __future__ import annotations

import numpy as np
from scipy.optimize import least_squares
from typing import NamedTuple


class NSParams(NamedTuple):
    """
    Nelson-Siegel (1987) yield curve parametrisation.

    r(τ) = β₀ + β₁·g₁(τ/λ) + β₂·g₂(τ/λ)

    where:
        g₁(x) = (1 − e^{−x}) / x                  [Level loading]
        g₂(x) = (1 − e^{−x}) / x − e^{−x}         [Curvature loading]

    Factor interpretation:
        β₀ > 0  : Long-run level    (r → β₀ as τ → ∞)
        β₁      : Short/long spread (r(0⁺) − r(∞) = β₁; β₁ < 0 → normal curve)
        β₂      : Hump/trough       (β₂ < 0 → mid-curve trough; β₂ > 0 → hump)
        λ > 0   : Decay speed       (g₂ peaks at τ* = λ·ln 2 ≈ 0.693λ years)

    Limitations:
        - One hump only. Cannot fit curves with two local extrema (use Svensson).
        - β₀ must be positive for positive long-run forward rates.
        - Extrapolation to τ → ∞ collapses to β₀ regardless of market structure.
    """
    beta0:   float
    beta1:   float
    beta2:   float
    lambda_: float


class SvenssonParams(NamedTuple):
    """
    Svensson (1994) extension: adds a second hump factor β₃ with decay λ₂.

    r(τ) = β₀ + β₁·g₁(τ/λ₁) + β₂·g₂(τ/λ₁) + β₃·g₂(τ/λ₂)

    Used by the ECB and most central banks for publishing official yield curves.
    """
    beta0:    float
    beta1:    float
    beta2:    float
    beta3:    float
    lambda1:  float
    lambda2:  float


def _g1(x: np.ndarray) -> np.ndarray:
    return np.where(x < 1e-9, 1.0, (1 - np.exp(-x)) / x)


def _g2(x: np.ndarray) -> np.ndarray:
    ex = np.exp(-np.clip(x, 0, 500))
    return np.where(x < 1e-9, 0.0, (1 - ex) / x - ex)


def ns_zero(tau: np.ndarray | float, p: NSParams) -> np.ndarray:
    """Nelson-Siegel zero-coupon yield. Continuously compounded, annualised."""
    tau = np.asarray(tau, dtype=float)
    x   = tau / p.lambda_
    return p.beta0 + p.beta1 * _g1(x) + p.beta2 * _g2(x)


def ns_forward(tau: np.ndarray | float, p: NSParams) -> np.ndarray:
    """Instantaneous forward rate f(τ) = β₀ + β₁·e^{−τ/λ} + β₂·(τ/λ)·e^{−τ/λ}."""
    tau = np.asarray(tau, dtype=float)
    x   = tau / p.lambda_
    ex  = np.exp(-np.clip(x, 0, 500))
    return p.beta0 + p.beta1 * ex + p.beta2 * x * ex


def svensson_zero(tau: np.ndarray | float, p: SvenssonParams) -> np.ndarray:
    """Svensson zero-coupon yield."""
    tau = np.asarray(tau, dtype=float)
    return (
        p.beta0
        + p.beta1 * _g1(tau / p.lambda1)
        + p.beta2 * _g2(tau / p.lambda1)
        + p.beta3 * _g2(tau / p.lambda2)
    )


def fit_nelson_siegel(
    maturities: np.ndarray,
    zero_rates:  np.ndarray,
    weights:     np.ndarray | None = None,
    n_restarts:  int = 8,
) -> tuple[NSParams, float]:
    """
    Fit Nelson-Siegel parameters to a bootstrapped zero-rate curve.

    Method: Levenberg-Marquardt on weighted residuals.
    Minimises: Σᵢ wᵢ · (r_NS(τᵢ; θ) − rᵢ)²

    Args:
        maturities: Pillar maturities (years).
        zero_rates:  Continuously compounded zero rates (decimal).
        weights:    Per-pillar weights (default: uniform).
        n_restarts: Random restarts to avoid local minima.

    Returns:
        (params, rmse_in_bps)
    """
    if weights is None:
        weights = np.ones_like(maturities)

    def residuals(x: np.ndarray) -> np.ndarray:
        p = NSParams(beta0=x[0], beta1=x[1], beta2=x[2], lambda_=max(x[3], 0.05))
        return weights * (ns_zero(maturities, p) - zero_rates)

    bounds = ([0.0, -0.20, -0.20, 0.1], [0.20, 0.20, 0.20, 10.0])

    rng     = np.random.default_rng(42)
    r_short = float(zero_rates[0])
    r_long  = float(np.mean(zero_rates[-3:]) if len(zero_rates) >= 3 else zero_rates[-1])
    best, best_cost = None, np.inf

    for i in range(n_restarts):
        x0 = (
            [max(r_long, 0.001), r_short - r_long, 0.0, 2.0]
            if i == 0 else
            [rng.uniform(0.001, 0.12), rng.uniform(-0.15, 0.15),
             rng.uniform(-0.15, 0.15), rng.uniform(0.5, 8.0)]
        )
        try:
            res = least_squares(residuals, x0, bounds=bounds, method="trf",
                                ftol=1e-12, xtol=1e-12, gtol=1e-12)
            if res.cost < best_cost:
                best, best_cost = res, res.cost
        except Exception:
            continue

    if best is None:
        raise RuntimeError("Nelson-Siegel fitting failed on all restarts")

    p = NSParams(*best.x[:3], lambda_=best.x[3])
    rmse_bps = float(np.sqrt(np.mean(residuals(best.x) ** 2))) * 10_000
    return p, rmse_bps

# sensitivities.py — DV01 per pillar via Jacobian bump-reval
from __future__ import annotations

from copy import deepcopy
from dataclasses import dataclass

import numpy as np

from curve_builder import Instrument, InstrumentType, DiscountCurve, bootstrap_ois


@dataclass
class PillarSensitivity:
    label:         str
    maturity:      float          # years
    zero_rate_pct: float          # base zero rate at this pillar, in %
    df:            float          # base discount factor P(0, T)
    dv01_eur:      float          # DV01 in EUR (10M notional zero-coupon bond)
    dv01_curve:    dict[str, float]  # ΔP(0,Tⱼ) per 1bp bump of this instrument


def dv01_ladder(
    ois_instruments: list[Instrument],
    bump_bp: float = 1.0,
    notional: float = 10_000_000,
) -> list[PillarSensitivity]:
    """
    DV01 ladder via Jacobian bump-reval.

    For each instrument i:
        1. Bump instrument i's rate by +bump_bp basis points.
        2. Re-bootstrap the OIS curve.
        3. DV01(i, j) = ΔP(0,Tⱼ) × notional / bump_bp  for all pillars j.

    This gives the Jacobian ∂P(0,Tⱼ) / ∂Kᵢ (per 1bp) — the full sensitivity
    matrix of bootstrapped discount factors to market instrument quotes.

    Analytic approximation (for a ZCB at pillar T):
        DV01(T) = notional · T · P(0,T) · 1bp
    This is the modified duration formula and is accurate for parallel shifts
    of the zero curve; the Jacobian method is exact for instrument-specific bumps.

    Args:
        ois_instruments: List of OIS market instruments.
        bump_bp:         Bump size in basis points.
        notional:        EUR notional per pillar (default 10M).

    Returns:
        List of PillarSensitivity, one per instrument.
    """
    bump    = bump_bp * 1e-4
    base    = bootstrap_ois(ois_instruments)
    result  = []

    for i, inst in enumerate(ois_instruments):
        bumped_insts       = deepcopy(ois_instruments)
        bumped_insts[i]    = Instrument(
            itype=inst.itype, rate=inst.rate + bump,
            maturity=inst.maturity, start=inst.start, label=inst.label,
        )
        bumped = bootstrap_ois(bumped_insts)

        # Jacobian row: change in each pillar DF per 1bp bump of instrument i
        dv01_curve: dict[str, float] = {}
        for t, P_base in base.pillars:
            delta_P = bumped.df(t) - P_base
            dv01_curve[f"{t:.3f}Y"] = delta_P * notional / bump_bp

        T = inst.maturity
        result.append(PillarSensitivity(
            label         = inst.label or f"OIS_{T:.2f}Y",
            maturity      = T,
            zero_rate_pct = base.zero_rate(T) * 100,
            df            = base.df(T),
            dv01_eur      = notional * T * base.df(T) * 1e-4,
            dv01_curve    = dv01_curve,
        ))

    return result

# market_data.py — Live EUR rate data from ECB SDMX REST API
from __future__ import annotations

from typing import Optional
import requests
import pandas as pd


ECB_BASE = "https://data-api.ecb.europa.eu/service/data"


def _fetch_ecb(flow: str, key: str, start: str, end: Optional[str] = None) -> pd.Series:
    """
    Fetch a single time series from the ECB Data Portal (SDMX-JSON format).

    Args:
        flow:  Dataset code (e.g. "FM" = financial markets, "EST" = €STR).
        key:   Series key within the dataset.
        start: ISO date string "YYYY-MM-DD" for start period.
        end:   ISO date string for end period (default: today).

    Returns:
        pd.Series with DatetimeIndex and float values in decimal (not %).

    Raises:
        requests.HTTPError on non-200 response.
        ValueError if the API returns an empty or malformed payload.

    Note on rate units:
        The ECB publishes EURIBOR and €STR in percent (e.g. 3.90 for 3.90%).
        This function divides by 100 before returning. Always check units
        when mixing sources.
    """
    params: dict[str, str] = {"format": "jsondata", "startPeriod": start}
    if end:
        params["endPeriod"] = end

    r = requests.get(f"{ECB_BASE}/{flow}/{key}", params=params, timeout=30)
    r.raise_for_status()
    data = r.json()

    try:
        # The series key varies by dataset; take the first (and only) series.
        series_key = next(iter(data["dataSets"][0]["series"]))
        obs        = data["dataSets"][0]["series"][series_key]["observations"]
        date_vals  = data["structure"]["dimensions"]["observation"][0]["values"]
    except (KeyError, IndexError, StopIteration) as e:
        raise ValueError(f"Unexpected ECB API response for {flow}/{key}: {e}") from e

    index  = pd.to_datetime([d["id"] for d in date_vals])
    values = [obs[str(i)][0] if str(i) in obs else float("nan") for i in range(len(date_vals))]
    s = pd.Series(values, index=index, dtype=float).dropna()

    if s.empty:
        raise ValueError(f"No data returned for {flow}/{key} from {start}")
    return s / 100.0   # percent → decimal


# ---------------------------------------------------------------------------
# Public API
# ---------------------------------------------------------------------------

EURIBOR_SERIES = {
    "1W":  "FM/B.U2.EUR.RT0.MM.EURIBOR1WD_.HSTA",
    "1M":  "FM/B.U2.EUR.RT0.MM.EURIBOR1MD_.HSTA",
    "3M":  "FM/B.U2.EUR.RT0.MM.EURIBOR3MD_.HSTA",
    "6M":  "FM/B.U2.EUR.RT0.MM.EURIBOR6MD_.HSTA",
    "12M": "FM/B.U2.EUR.RT0.MM.EURIBOR12MD_.HSTA",
}

ECB_YC_SERIES = {
    "3M":  "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_3M",
    "6M":  "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_6M",
    "1Y":  "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_1Y",
    "2Y":  "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_2Y",
    "5Y":  "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_5Y",
    "10Y": "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_10Y",
    "20Y": "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_20Y",
    "30Y": "YC/B.U2.EUR.4F.G_N_A.SV_C_YM.SR_30Y",
}


def fetch_euribor(tenor: str = "3M", start: str = "2020-01-01",
                  end: Optional[str] = None) -> pd.Series:
    """
    EURIBOR fixing from the ECB Statistical Data Warehouse.

    Convention: ACT/360, simple compounding.
    Tenors available: "1W", "1M", "3M", "6M", "12M".

    Example:
        euribor3m = fetch_euribor("3M", start="2023-01-01")
        print(f"Latest 3M EURIBOR: {euribor3m.iloc[-1]:.4%}")
    """
    if tenor not in EURIBOR_SERIES:
        raise ValueError(f"Tenor {tenor!r} not in {list(EURIBOR_SERIES)}")
    flow, key = EURIBOR_SERIES[tenor].split("/", 1)
    return _fetch_ecb(flow, key, start, end)


def fetch_estr(start: str = "2020-01-01", end: Optional[str] = None) -> pd.Series:
    """
    Euro Short-Term Rate (€STR) from the ECB.

    €STR replaced EONIA effective 2 January 2022.
    Convention: ACT/360, simple compounding.
    Published by ECB by 08:00 CET each TARGET2 business day.

    Historical note: EONIA = €STR + 8.5bp was maintained transitionally until
    3 January 2022. For historical back-fill before 2 October 2019 (€STR launch),
    use EONIA − 8.5bp as an approximation.
    """
    return _fetch_ecb("EST", "B.EU000A2X2A25.WT", start, end)


def fetch_ecb_yield_curve(start: str = "2020-01-01",
                           end: Optional[str] = None) -> pd.DataFrame:
    """
    ECB AAA-rated government bond spot yield curve at standard tenors.

    The ECB fits a Svensson (1994) model daily to a basket of euro-area
    government bonds rated AAA. This function extracts spot rates at
    standard tenors from the published curve.

    Returns:
        DataFrame with DatetimeIndex and columns = tenor strings
        (e.g. "3M", "1Y", "10Y"). Rates in decimal.

    Usage:
        df = fetch_ecb_yield_curve(start="2023-01-01")
        latest = df.iloc[-1]   # most recent yield curve
        print(latest * 100)    # in percent

    Reference:
        https://www.ecb.europa.eu/stats/financial_markets_and_interest_rates/
        euro_area_yield_curves/html/index.en.html
    """
    frames: list[pd.Series] = []
    for label, series_id in ECB_YC_SERIES.items():
        flow, key = series_id.split("/", 1)
        try:
            s = _fetch_ecb(flow, key, start, end)
            s.name = label
            frames.append(s)
        except Exception as e:
            print(f"Warning: could not fetch {label}: {e}")
            continue

    if not frames:
        raise ValueError("No ECB yield curve data returned")
    return pd.DataFrame(frames).T.sort_index()


def build_ois_instruments_from_estr(
    estr_rate: float,
    swap_quotes: dict[str, float],
) -> list["Instrument"]:
    """
    Construct a list of OIS Instrument objects from live market data.

    Args:
        estr_rate:   Latest €STR fixing (decimal).
        swap_quotes: Dict of {tenor_label: par_rate_decimal}, e.g.:
                     {"1M": 0.039, "3M": 0.038, "6M": 0.037, "1Y": 0.035, ...}

    Returns:
        List[Instrument] ready for bootstrap_ois().

    Tenor to maturity mapping (approximate; use QuantLib for exact dates):
        "1M" → 1/12, "3M" → 0.25, "6M" → 0.5, "1Y" → 1.0, etc.
    """
    from curve_builder import Instrument, InstrumentType

    tenor_to_years = {
        "ON": 1/365, "TN": 2/365, "1W": 1/52,
        "1M": 1/12,  "3M": 0.25,  "6M": 0.5,
        "9M": 0.75,  "1Y": 1.0,   "2Y": 2.0,
        "3Y": 3.0,   "5Y": 5.0,   "7Y": 7.0,
        "10Y": 10.0, "15Y": 15.0, "20Y": 20.0, "30Y": 30.0,
    }

    insts: list[Instrument] = [
        Instrument(InstrumentType.DEPOSIT, rate=estr_rate,
                   maturity=1/365, label="ESTR_ON"),
    ]
    for tenor, rate in swap_quotes.items():
        t = tenor_to_years.get(tenor)
        if t is None:
            continue
        itype = InstrumentType.DEPOSIT if t <= 0.25 else InstrumentType.OIS_SWAP
        insts.append(Instrument(itype=itype, rate=rate, maturity=t, label=f"OIS_{tenor}"))

    return sorted(insts, key=lambda i: i.maturity)

# visualisation.py — Plotly interactive yield curve charts
from __future__ import annotations

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from curve_builder import DiscountCurve
from nelson_siegel import NSParams, ns_zero, ns_forward, SvenssonParams, svensson_zero


TAU_GRID = np.linspace(0.25, 30, 300)


def plot_curves(
    ois_curve: DiscountCurve,
    euribor_curve: DiscountCurve | None = None,
    ns_params: NSParams | None = None,
) -> go.Figure:
    """
    Interactive yield curve: zero rates + instantaneous forward rates.

    Panels:
        Left  — Zero coupon rates (OIS, EURIBOR 3M, Nelson-Siegel fit)
        Right — Instantaneous forward rates (same three curves)

    All rates in % on y-axis. Continuously compounded, ACT/365 convention.
    """
    fig = make_subplots(
        rows=1, cols=2,
        subplot_titles=["Zero Coupon Rates (%)", "Instantaneous Forward Rates (%)"],
        column_widths=[0.5, 0.5],
    )

    ois_z = [ois_curve.zero_rate(t) * 100 for t in TAU_GRID]
    ois_f = [ois_curve.inst_forward(t) * 100 for t in TAU_GRID]

    for col, ys in enumerate([ois_z, ois_f], start=1):
        fig.add_trace(go.Scatter(
            x=TAU_GRID, y=ys, mode="lines", name="OIS (ESTR)",
            line=dict(color="#6366f1", width=2.5), showlegend=(col == 1),
        ), row=1, col=col)

    # Pillar markers on zero curve
    for t, P in ois_curve.pillars:
        fig.add_trace(go.Scatter(
            x=[t], y=[-np.log(P) / t * 100], mode="markers",
            marker=dict(color="#6366f1", size=7, symbol="circle-open", line=dict(width=2)),
            name="OIS pillars", showlegend=False,
        ), row=1, col=1)

    if euribor_curve is not None and len(euribor_curve.pillars) > 1:
        eib_z = [euribor_curve.zero_rate(t) * 100 for t in TAU_GRID]
        eib_f = [euribor_curve.inst_forward(t) * 100 for t in TAU_GRID]
        for col, ys in enumerate([eib_z, eib_f], start=1):
            fig.add_trace(go.Scatter(
                x=TAU_GRID, y=ys, mode="lines", name="EURIBOR 3M",
                line=dict(color="#10b981", width=2.5), showlegend=(col == 1),
            ), row=1, col=col)

        # EURIBOR-OIS spread (basis) overlay on left panel
        spread_bps = [(e - o) * 100 for e, o in zip(eib_z, ois_z)]
        fig.add_trace(go.Scatter(
            x=TAU_GRID, y=spread_bps, mode="lines",
            name="EURIBOR-OIS spread (%)",
            line=dict(color="#f59e0b", width=1.5, dash="dot"),
            showlegend=True,
        ), row=1, col=1)

    if ns_params is not None:
        ns_z = ns_zero(TAU_GRID, ns_params) * 100
        ns_f = ns_forward(TAU_GRID, ns_params) * 100
        for col, ys in enumerate([ns_z, ns_f], start=1):
            fig.add_trace(go.Scatter(
                x=TAU_GRID, y=ys, mode="lines", name="Nelson-Siegel fit",
                line=dict(color="#f97316", width=2, dash="dot"), showlegend=(col == 1),
            ), row=1, col=col)

    fig.update_xaxes(title_text="Maturity τ (years)")
    fig.update_yaxes(title_text="Rate (%)", row=1, col=1)
    fig.update_yaxes(title_text="Forward Rate (%)", row=1, col=2)
    fig.update_layout(
        height=430, template="plotly_dark",
        title="EUR Yield Curves — OIS & EURIBOR 3M",
        hovermode="x unified",
        legend=dict(orientation="h", yanchor="bottom", y=1.02),
    )
    return fig


def plot_dv01_ladder(
    curve: DiscountCurve,
    pillars: list[float] | None = None,
    notional: float = 10_000_000,
) -> go.Figure:
    """
    Bar chart of DV01 per pillar (EUR per 1bp, ZCB approximation).

    DV01(τ) = Notional · τ · P(0,τ) · 1bp

    This is the modified duration formula for a zero-coupon bond. For coupon
    bonds or par swaps, use the full Jacobian from sensitivities.dv01_ladder().

    Args:
        curve:    Bootstrapped OIS (or EURIBOR) discount curve.
        pillars:  Maturities in years (default: standard EUR swap tenors).
        notional: EUR notional per pillar.
    """
    if pillars is None:
        pillars = [0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30]

    labels = [f"{int(t*12)}M" if t < 1 else f"{int(t)}Y" for t in pillars]
    dv01s  = [notional * t * curve.df(t) * 1e-4 for t in pillars]
    zeros  = [curve.zero_rate(t) * 100 for t in pillars]
    n      = len(pillars)
    colors = [f"hsl({int(240 - i * 200 / (n-1))}, 65%, 55%)" for i in range(n)]

    fig = make_subplots(
        rows=2, cols=1,
        subplot_titles=["DV01 per Pillar (EUR, 10M notional)", "Zero Coupon Rate (%)"],
        row_heights=[0.65, 0.35], vertical_spacing=0.12,
    )
    fig.add_trace(go.Bar(
        x=labels, y=dv01s, marker=dict(color=colors),
        text=[f"€{d:,.0f}" for d in dv01s], textposition="outside",
        name="DV01",
    ), row=1, col=1)
    fig.add_trace(go.Scatter(
        x=labels, y=zeros, mode="lines+markers",
        line=dict(color="#6366f1", width=2), marker=dict(size=6),
        name="Zero Rate",
    ), row=2, col=1)

    fig.update_yaxes(title_text="DV01 (EUR)", row=1, col=1)
    fig.update_yaxes(title_text="Rate (%)", row=2, col=1)
    fig.update_layout(
        height=520, template="plotly_dark",
        title="OIS Curve — DV01 Ladder", showlegend=False,
    )
    return fig


def plot_historical_surface(df: pd.DataFrame) -> go.Figure:
    """
    3-D yield surface: rate = f(maturity, date).

    Args:
        df: DataFrame from fetch_ecb_yield_curve(); index = dates,
            columns = tenor strings, values = rates in decimal.
    """
    def to_years(s: str) -> float:
        return int(s[:-1]) / 12 if s.endswith("M") else float(s[:-1])

    x = [to_years(c) for c in df.columns]
    z = df.values * 100
    date_labels = [d.strftime("%Y-%m-%d") for d in df.index]

    fig = go.Figure(data=[go.Surface(
        x=x, y=list(range(len(df))), z=z,
        colorscale="RdYlBu_r",
        colorbar=dict(title="Rate (%)"),
        hovertemplate="Maturity: %{x:.1f}Y<br>Rate: %{z:.3f}%<extra></extra>",
    )])
    tick_step = max(1, len(df) // 10)
    fig.update_layout(
        title="EUR Yield Surface — ECB AAA Government Bonds",
        scene=dict(
            xaxis_title="Maturity (years)",
            yaxis=dict(
                title="Date",
                tickvals=list(range(0, len(df), tick_step)),
                ticktext=date_labels[::tick_step],
            ),
            zaxis_title="Rate (%)",
            camera=dict(eye=dict(x=1.8, y=-1.5, z=0.8)),
        ),
        height=560, template="plotly_dark",
    )
    return fig


def plot_animated_curve(df: pd.DataFrame) -> go.Figure:
    """
    Animated yield curve with a date slider. Useful for
    visualising ECB rate-hiking / cutting cycles.

    Args:
        df: Same format as fetch_ecb_yield_curve() output.
    """
    def to_years(s: str) -> float:
        return int(s[:-1]) / 12 if s.endswith("M") else float(s[:-1])

    x = [to_years(c) for c in df.columns]
    frames = [
        go.Frame(
            data=[go.Scatter(x=x, y=df.iloc[i].values * 100,
                             mode="lines+markers",
                             line=dict(color="#6366f1", width=2.5),
                             marker=dict(size=6))],
            name=df.index[i].strftime("%Y-%m-%d"),
            layout=go.Layout(title_text=f"EUR Yield Curve — {df.index[i].strftime('%Y-%m-%d')}"),
        )
        for i in range(len(df))
    ]
    fig = go.Figure(
        data=frames[0].data,
        frames=frames,
        layout=go.Layout(
            title=f"EUR Yield Curve — {df.index[-1].strftime('%Y-%m-%d')}",
            xaxis_title="Maturity (years)",
            yaxis_title="Spot Rate (%)",
            template="plotly_dark", height=430,
            updatemenus=[dict(type="buttons", buttons=[
                dict(label="▶ Play",  method="animate",
                     args=[None, dict(frame=dict(duration=180, redraw=True), fromcurrent=True)]),
                dict(label="⏸ Pause", method="animate",
                     args=[[None], dict(frame=dict(duration=0), mode="immediate")]),
            ])],
            sliders=[dict(
                steps=[dict(args=[[f.name], dict(frame=dict(duration=180), mode="immediate")],
                            label=f.name, method="animate") for f in frames],
                currentvalue=dict(prefix="Date: "), len=0.9, x=0.05,
            )],
        ),
    )
    return fig

# tests/test_curve_builder.py — pytest test suite
# Run: pytest tests/ -v
import numpy as np
import pytest

from curve_builder import (
    Instrument, InstrumentType, DiscountCurve,
    bootstrap_ois, bootstrap_euribor, _act360,
)
from nelson_siegel import (
    NSParams, ns_zero, ns_forward, fit_nelson_siegel,
)


# ── Fixtures ──────────────────────────────────────────────────────────────────

@pytest.fixture
def ois_instruments() -> list[Instrument]:
    """EUR OIS market snapshot (approximate 2024 inverted curve)."""
    return [
        Instrument(InstrumentType.DEPOSIT,  rate=0.0390, maturity=1/12,  label="DEP_1M"),
        Instrument(InstrumentType.DEPOSIT,  rate=0.0388, maturity=0.25,  label="DEP_3M"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0370, maturity=0.5,   label="OIS_6M"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0350, maturity=1.0,   label="OIS_1Y"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0310, maturity=2.0,   label="OIS_2Y"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0280, maturity=5.0,   label="OIS_5Y"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0265, maturity=10.0,  label="OIS_10Y"),
        Instrument(InstrumentType.OIS_SWAP, rate=0.0255, maturity=30.0,  label="OIS_30Y"),
    ]

@pytest.fixture
def euribor_instruments() -> list[Instrument]:
    """EUR EURIBOR 3M curve instruments (approximate 2024 levels)."""
    return [
        Instrument(InstrumentType.DEPOSIT, rate=0.0393, maturity=0.25, label="EURIBOR_3M"),
        Instrument(InstrumentType.FRA,     rate=0.0375, maturity=0.50, start=0.25, label="3x6"),
        Instrument(InstrumentType.FRA,     rate=0.0360, maturity=0.75, start=0.50, label="6x9"),
        Instrument(InstrumentType.FRA,     rate=0.0345, maturity=1.00, start=0.75, label="9x12"),
        Instrument(InstrumentType.IRS,     rate=0.0310, maturity=2.0,  label="IRS_2Y"),
        Instrument(InstrumentType.IRS,     rate=0.0285, maturity=5.0,  label="IRS_5Y"),
        Instrument(InstrumentType.IRS,     rate=0.0268, maturity=10.0, label="IRS_10Y"),
    ]


# ── Discount curve fundamentals ────────────────────────────────────────────────

class TestDiscountCurve:
    def test_df_at_zero(self, ois_instruments):
        curve = bootstrap_ois(ois_instruments)
        assert abs(curve.df(0.0) - 1.0) < 1e-10, "P(0,0) must equal 1"

    def test_dfs_monotone_decreasing(self, ois_instruments):
        curve = bootstrap_ois(ois_instruments)
        ts    = [0.25, 0.5, 1, 2, 5, 10, 20]
        dfs   = [curve.df(t) for t in ts]
        for i in range(len(dfs) - 1):
            assert dfs[i] > dfs[i+1], (
                f"Non-monotone: P({ts[i]}) = {dfs[i]:.6f} ≤ P({ts[i+1]}) = {dfs[i+1]:.6f}"
            )

    def test_zero_rates_positive(self, ois_instruments):
        curve = bootstrap_ois(ois_instruments)
        for t in [0.5, 1, 5, 10]:
            assert curve.zero_rate(t) > 0, f"Negative zero rate at t={t}"

    def test_forward_rates_positive(self, ois_instruments):
        """No negative simply-compounded forward rates (necessary for no-arb)."""
        curve = bootstrap_ois(ois_instruments)
        for t in [0.25, 0.5, 1, 2, 5, 10]:
            f = curve.forward_rate(t, t + 0.25)
            assert f > 0, f"Negative forward F(0;{t},{t+0.25}) = {f:.6f}"


# ── OIS bootstrap self-consistency ────────────────────────────────────────────

class TestOISBootstrap:
    def test_deposit_reprices_at_par(self, ois_instruments):
        """
        A bootstrapped deposit must satisfy: (1 + K·δ) · P(0,T) = 1.
        This is the definition of par repricing.
        """
        curve = bootstrap_ois(ois_instruments)
        for inst in ois_instruments:
            if inst.itype == InstrumentType.DEPOSIT:
                T   = inst.maturity
                pv  = (1 + inst.rate * _act360(T)) * curve.df(T)
                assert abs(pv - 1.0) < 1e-8, (
                    f"Deposit {T:.3f}Y PV = {pv:.10f} ≠ 1"
                )

    def test_ois_swap_reprices_at_par(self, ois_instruments):
        """Par OIS swap: fixed leg PV = float leg PV = P(0,0) − P(0,T) = 1 − P(0,T)."""
        curve = bootstrap_ois(ois_instruments)
        for inst in ois_instruments:
            if inst.itype == InstrumentType.OIS_SWAP:
                T      = inst.maturity
                K      = inst.rate
                n      = max(1, int(round(T)))
                T_i    = [min(k, T) for k in range(1, n + 1)]
                d_i    = [_act360(T_i[0])] + [_act360(T_i[k] - T_i[k-1]) for k in range(1, n)]
                fixed  = K * sum(d * curve.df(t) for d, t in zip(d_i, T_i))
                flt    = 1.0 - curve.df(T)
                assert abs(fixed - flt) < 1e-6, (
                    f"OIS {T}Y NPV = {fixed - flt:.2e} (should be 0)"
                )


# ── Multi-curve bootstrap ─────────────────────────────────────────────────────

class TestMultiCurve:
    def test_euribor_above_ois(self, ois_instruments, euribor_instruments):
        """EURIBOR zero rates must exceed OIS zero rates (positive credit/liquidity basis)."""
        ois   = bootstrap_ois(ois_instruments)
        eibor = bootstrap_euribor(euribor_instruments, ois)
        for t in [0.25, 0.5, 1.0, 2.0]:
            assert eibor.zero_rate(t) >= ois.zero_rate(t) - 1e-6, (
                f"EURIBOR < OIS at t={t}: {eibor.zero_rate(t):.4%} < {ois.zero_rate(t):.4%}"
            )

    def test_euribor_deposit_reprices(self, ois_instruments, euribor_instruments):
        ois   = bootstrap_ois(ois_instruments)
        eibor = bootstrap_euribor(euribor_instruments, ois)
        dep   = next(i for i in euribor_instruments if i.itype == InstrumentType.DEPOSIT)
        pv    = (1 + dep.rate * _act360(dep.maturity)) * eibor.df(dep.maturity)
        assert abs(pv - 1.0) < 1e-7, f"EURIBOR deposit PV = {pv:.10f}"


# ── Nelson-Siegel ─────────────────────────────────────────────────────────────

class TestNelsonSiegel:
    def test_short_rate_limit(self):
        """r(0⁺) = β₀ + β₁ (g₁ → 1, g₂ → 0 as x → 0)."""
        p = NSParams(beta0=0.04, beta1=-0.015, beta2=0.008, lambda_=2.5)
        r0 = float(ns_zero(1e-9, p))
        assert abs(r0 - (p.beta0 + p.beta1)) < 1e-6

    def test_long_rate_limit(self):
        """r(∞) = β₀ (g₁ → 0, g₂ → 0 as x → ∞)."""
        p = NSParams(beta0=0.04, beta1=-0.015, beta2=0.008, lambda_=2.5)
        assert abs(float(ns_zero(1000.0, p)) - p.beta0) < 1e-6

    def test_hump_location(self):
        """The curvature factor g₂ peaks at τ* = λ · ln 2."""
        p    = NSParams(beta0=0.04, beta1=0.0, beta2=1.0, lambda_=2.0)
        tau* = p.lambda_ * np.log(2)
        taus = np.linspace(0.1, 10, 500)
        # The NS forward rate with β₀=β₁=0 is just β₂·(τ/λ)·e^{−τ/λ},
        # which peaks at τ = λ. The NS zero rate (g₂) peaks near τ = λ·ln2.
        g2_vals = np.where(
            taus / p.lambda_ < 1e-9, 0.0,
            (1 - np.exp(-taus / p.lambda_)) / (taus / p.lambda_) - np.exp(-taus / p.lambda_)
        )
        peak_tau = taus[np.argmax(g2_vals)]
        assert abs(peak_tau - tau*) < 0.15, f"g₂ peak at {peak_tau:.3f}, expected {tau*:.3f}"

    def test_roundtrip_fitting(self):
        """Generate zero rates from known NS params, fit, verify recovery."""
        true = NSParams(beta0=0.04, beta1=-0.015, beta2=0.008, lambda_=2.5)
        taus = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30])
        rates = ns_zero(taus, true)
        fitted, rmse_bps = fit_nelson_siegel(taus, rates)
        assert rmse_bps < 0.1, f"Fit RMSE {rmse_bps:.3f}bps > 0.1bps"
        np.testing.assert_allclose(ns_zero(taus, fitted), rates, atol=5e-6)

    @pytest.mark.parametrize("t", [0.25, 1.0, 5.0, 10.0, 30.0])
    def test_discount_factors_in_unit_interval(self, t):
        p   = NSParams(beta0=0.04, beta1=-0.015, beta2=0.008, lambda_=2.5)
        P   = float(np.exp(-float(ns_zero(t, p)) * t))
        assert 0 < P <= 1, f"P(0,{t}) = {P:.6f} out of (0,1]"