cppi.py

from __future__ import annotations
import numpy as np
import pandas as pd
from dataclasses import dataclass

from scipy.stats import norm


@dataclass(frozen=True)
class Params:
    V0: float = 100.0
    G: float = 100.0
    T: float = 5.0
    r: float = 0.04
    mu: float = 0.08
    sigma: float = 0.20
    m: float = 4.0
    N_mc: int = 20_000
    steps_per_year: int = 252
    tc_bps: float = 10.0

    @property
    def M(self) -> int:
        return int(self.T * self.steps_per_year)

    @property
    def dt(self) -> float:
        return self.T / self.M

    @property
    def tc(self) -> float:
        return self.tc_bps * 1e-4


def metrics(VT: np.ndarray, G: float, label: str, materiality: float = 1.0) -> dict:
    """Compute risk metrics on terminal portfolio values.

    Parameters
    ----------
    VT : np.ndarray of shape (N_mc,)
        Simulated terminal portfolio values.
    G : float
        Floor (guaranteed amount).
    label : str
        Identifier of the model/scenario, included in the output dict.
    materiality : float, default 1.0
        Threshold below the floor to define a material breach.

    Returns
    -------
    dict
        Mapping of metric names to values: mean, standard deviation,
        quantiles (0.5%, 5%, 95%), raw and material breach probabilities,
        conditional expected shortfall on breaches, and tail-5% expected shortfall.
    """

    VT = np.asarray(VT)
    breach = VT < G
    breach_mat = VT < (G - materiality)
    loss = np.maximum(G - VT, 0.0)
    p_breach = breach.mean()
    p_breach_mat = breach_mat.mean()
    es_cond = loss[breach].mean() if breach.any() else 0.0
    q05 = np.quantile(VT, 0.05)
    tail_losses = np.maximum(G - VT[VT <= q05], 0.0)
    es_tail5 = tail_losses.mean() if len(tail_losses) else 0.0
    return {
        "Model": label,
        "E[V_T]": VT.mean(),
        "Std(V_T)": VT.std(),
        "q_0.5%": np.quantile(VT, 0.005),
        "q_5%": q05,
        "q_95%": np.quantile(VT, 0.95),
        "P(V_T<G) % raw": 100 * p_breach,
        "P(V_T<G) % material": 100 * p_breach_mat,
        "ES_cond": es_cond,
        "ES_tail5%": es_tail5,
    }


def bs_closed_form(p: Params, n_sim: int = 20_000, seed: int = 42) -> np.ndarray:
    """Simulate terminal CPPI values under continuous-time Black-Scholes.

    Uses the closed-form solution for the CPPI value at maturity, which
    expresses V_T as a power function of the terminal underlying S_T.
    No transaction costs, no discretization error.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    n_sim : int, default 20_000
        Number of Monte Carlo paths.
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (n_sim,)
        Terminal portfolio values.
    """
    rng = np.random.default_rng(seed)
    Z = rng.standard_normal(n_sim)
    ST = p.V0 * np.exp((p.mu - 0.5 * p.sigma**2) * p.T + p.sigma * np.sqrt(p.T) * Z)
    Cm = np.exp((1 - p.m) * (p.r + p.m * p.sigma**2 / 2) * p.T)
    VT = p.G + (p.V0 - p.G * np.exp(-p.r * p.T)) * Cm * (ST / p.V0) ** p.m
    return VT


def p_lsf_bs(p: Params, rebal_every: int) -> tuple[float, float]:
    """Compute local and global breach probabilities under discrete BS.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    rebal_every : int
        Rebalancing interval in days.

    Returns
    -------
    p_lsf : float
        Local single-step failure probability.
    p_bf : float
        Global breach probability over [0, T], assuming independence
        between rebalancing periods.
    """
    dt_r = rebal_every / p.steps_per_year
    N = int(p.T / dt_r)
    d2 = (np.log(p.m / (p.m - 1)) + (p.mu - p.r) * dt_r - 0.5 * p.sigma**2 * dt_r) / (
        p.sigma * np.sqrt(dt_r)
    )
    p_lsf = norm.cdf(-d2)
    p_bf = 1 - (1 - p_lsf) ** N
    return p_lsf, p_bf


def p_bf_kou(p: Params, lam: float, p_neg: float, eta_minus: float) -> float:
    """Theoretical upper bound on global breach probability under Kou.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    lam : float
        Jump intensity (jumps per year).
    p_neg : float
        Probability of a downward jump.
    eta_minus : float
        Rate parameter of the exponential distribution of downward jumps.

    Returns
    -------
    float
        Upper bound on P(V_T < G) over [0, T].
    """

    return 1 - np.exp(-p.T * lam * p_neg * (1 - 1 / p.m) ** eta_minus)


def cppi_bs_discrete(p: Params, rebal_every: int = 1, seed: int = 42) -> np.ndarray:
    """Simulate CPPI under discrete-time Black-Scholes with transaction costs.

    Parameters
    ----------
    p : Params
        Portfolio parameters (V0, G, T, m, sigma, tc_bps, ...).
    rebal_every : int, default 1
        Rebalancing frequency in time steps (1 = daily).
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (N_mc,)
        Terminal portfolio values.
    """

    rng = np.random.default_rng(seed)
    V = np.full(p.N_mc, p.V0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)

    for k in range(p.M):
        t = k * p.dt
        if k % rebal_every == 0:
            F = p.G * np.exp(-p.r * (p.T - t))
            C = np.maximum(V - F, 0.0)
            e_new = np.minimum(p.m * C, V)
            tc_cost = p.tc * np.abs(e_new - e)
            V = V - tc_cost
            e = e_new
            cash = V - e

        Z = rng.standard_normal(p.N_mc)
        S_ret = np.exp((p.mu - 0.5 * p.sigma**2) * p.dt + p.sigma * np.sqrt(p.dt) * Z)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)
    return V


def cppi_heston(
    p: Params,
    rebal_every: int = 1,
    v0: float = 0.04,
    kappa: float = 0.5,
    theta: float = 0.09,
    xi: float = 1.0,
    rho: float = -0.9,
    seed: int = 42,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate CPPI under the Heston stochastic volatility model.

    Variance follows a mean-reverting square-root process correlated with
    the asset Brownian motion. Variance is truncated at zero (full
    truncation scheme) to avoid negative values from the Euler discretization.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    rebal_every : int, default 1
        Rebalancing frequency in time steps.
    v0 : float, default 0.04
        Initial variance.
    kappa : float, default 0.5
        Mean reversion speed.
    theta : float, default 0.09
        Long-term variance.
    xi : float, default 1.0
        Volatility of variance.
    rho : float, default -0.9
        Correlation between asset and variance Brownian motions.
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (N_mc,)
        Terminal portfolio values.
    realized_vol : np.ndarray of shape (N_mc,)
        Path-wise realized volatility over [0, T].
    """

    rng = np.random.default_rng(seed)
    V = np.full(p.N_mc, p.V0)
    v = np.full(p.N_mc, v0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)
    sum_v = np.zeros(p.N_mc)

    for k in range(p.M):
        t = k * p.dt
        if k % rebal_every == 0:
            F = p.G * np.exp(-p.r * (p.T - t))
            C = np.maximum(V - F, 0.0)
            e_new = np.minimum(p.m * C, V)
            tc_cost = p.tc * np.abs(e_new - e)
            V = V - tc_cost
            e = e_new
            cash = V - e

        Z1 = rng.standard_normal(p.N_mc)
        Z2 = rng.standard_normal(p.N_mc)
        dW_s = np.sqrt(p.dt) * Z1
        dW_v = np.sqrt(p.dt) * (rho * Z1 + np.sqrt(1 - rho**2) * Z2)
        v_pos = np.maximum(v, 0.0)
        sum_v += v_pos * p.dt
        v = v + kappa * (theta - v_pos) * p.dt + xi * np.sqrt(v_pos) * dW_v
        v = np.maximum(v, 0.0)

        S_ret = np.exp((p.mu - 0.5 * v_pos) * p.dt + np.sqrt(v_pos) * dW_s)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)

    realized_vol = np.sqrt(sum_v / p.T)
    return V, realized_vol


def _kou_compensator(lam: float, p_neg: float, eta_up: float, eta_dn: float) -> float:
    """Martingale compensator of the Kou jump process."""

    p_up = 1.0 - p_neg
    return lam * (p_up * eta_up / (eta_up - 1) + p_neg * eta_dn / (eta_dn + 1) - 1)


def cppi_kou(
    p: Params,
    rebal_every: int = 1,
    lam: float = 1.0,
    p_neg: float = 0.6,
    eta_up: float = 25.0,
    eta_dn: float = 10.0,
    seed: int = 42,
) -> np.ndarray:
    """Simulate CPPI under the Kou jump-diffusion model.

    Log-returns combine a diffusion part with a compound Poisson process
    whose jump sizes follow a double-exponential distribution. The drift
    is adjusted by the Kou compensator to ensure the discounted asset
    is a martingale under the risk-neutral measure.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    rebal_every : int, default 1
        Rebalancing frequency in time steps.
    lam : float, default 1.0
        Jump intensity (jumps per year).
    p_neg : float, default 0.6
        Probability of a downward jump.
    eta_up : float, default 25.0
        Rate of the exponential distribution for upward jump sizes.
    eta_dn : float, default 10.0
        Rate of the exponential distribution for downward jump sizes.
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (N_mc,)
        Terminal portfolio values.
    """

    rng = np.random.default_rng(seed)
    comp = _kou_compensator(lam, p_neg, eta_up, eta_dn)
    V = np.full(p.N_mc, p.V0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)

    for k in range(p.M):
        t = k * p.dt
        if k % rebal_every == 0:
            F = p.G * np.exp(-p.r * (p.T - t))
            C = np.maximum(V - F, 0.0)
            e_new = np.minimum(p.m * C, V)
            tc_cost = p.tc * np.abs(e_new - e)
            V = V - tc_cost
            e = e_new
            cash = V - e

        Z = rng.standard_normal(p.N_mc)
        logret = (p.mu - 0.5 * p.sigma**2 - comp) * p.dt + p.sigma * np.sqrt(p.dt) * Z

        Nj = rng.poisson(lam * p.dt, size=p.N_mc)
        mask = Nj > 0
        if mask.any():
            nj = int(mask.sum())
            p_up = 1.0 - p_neg
            up = rng.random(nj) < p_up
            J = np.where(
                up, rng.exponential(1 / eta_up, nj), -rng.exponential(1 / eta_dn, nj)
            )
            logret[mask] += J

        S_ret = np.exp(logret)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)
    return V


def cppi_kou_var_ewma(
    p: Params,
    alpha: float = 0.995,
    lambda_ewma: float = 0.94,
    m_max: float | None = None,
    m_min: float = 1.0,
    rebal_every: int = 1,
    lam: float = 1.0,
    p_neg: float = 0.6,
    eta_up: float = 25.0,
    eta_dn: float = 10.0,
    seed: int = 42,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate CPPI under Kou with a dynamic VaR+EWMA multiplier.

    The multiplier m_t is recomputed at each rebalancing date using a
    VaR-based formula, with volatility estimated by an EWMA on observed
    log-returns. The underlying jumps follow a Kou double-exponential process.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    alpha : float, default 0.995
        VaR confidence level.
    lambda_ewma : float, default 0.94
        EWMA smoothing factor (RiskMetrics convention).
    m_max : float or None, default None
        Upper cap on the multiplier. None means no cap.
    m_min : float, default 1.0
        Lower bound on the multiplier.
    rebal_every : int, default 1
        Rebalancing frequency in time steps.
    lam : float, default 1.0
        Jump intensity (jumps per year).
    p_neg : float, default 0.6
        Probability of a downward jump.
    eta_up, eta_dn : float, default 25.0, 10.0
        Rates of the exponential distributions for upward/downward jump sizes.
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (N_mc,)
        Terminal portfolio values.
    m_history : np.ndarray of shape (M,)
        Multiplier trajectory m_t for the first simulation
        (visualization only, not for statistical analysis).
    """
    rng = np.random.default_rng(seed)
    comp = _kou_compensator(lam, p_neg, eta_up, eta_dn)
    z_alpha = norm.ppf(alpha)

    V = np.full(p.N_mc, p.V0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)

    sigma_ewma = np.full(p.N_mc, p.sigma**2 * p.dt)
    prev_logret = np.zeros(p.N_mc)
    m_history = []

    for k in range(p.M):
        t = k * p.dt

        if k > 0:
            sigma_ewma = lambda_ewma * sigma_ewma + (1 - lambda_ewma) * prev_logret**2

        sigma_t = np.sqrt(sigma_ewma / p.dt)

        tau = max(p.T - t, 1.0 / 12.0)

        exponent = (p.mu - p.r - 0.5 * sigma_t**2) * tau - sigma_t * z_alpha * np.sqrt(
            tau
        )
        denom = 1 - np.exp(exponent)
        denom = np.where(denom > 1e-8, denom, 1e-8)
        m_t = 1.0 / denom
        m_t = np.clip(m_t, m_min, m_max if m_max is not None else 1e9)
        m_history.append(float(m_t[0]))

        if k % rebal_every == 0:
            F = p.G * np.exp(-p.r * (p.T - t))
            C = np.maximum(V - F, 0.0)
            e_new = np.minimum(m_t * C, V)
            tc_cost = p.tc * np.abs(e_new - e)
            V = V - tc_cost
            e = e_new
            cash = V - e

        Z = rng.standard_normal(p.N_mc)
        logret = (p.mu - 0.5 * p.sigma**2 - comp) * p.dt + p.sigma * np.sqrt(p.dt) * Z
        Nj = rng.poisson(lam * p.dt, size=p.N_mc)
        mask = Nj > 0
        if mask.any():
            nj = int(mask.sum())
            p_up = 1.0 - p_neg
            up = rng.random(nj) < p_up
            J = np.where(
                up, rng.exponential(1 / eta_up, nj), -rng.exponential(1 / eta_dn, nj)
            )
            logret[mask] += J

        S_ret = np.exp(logret)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)
        prev_logret = logret

    return V, np.array(m_history)


def _bs_put_price(S, K, T, r, sigma):
    """Black-Scholes price of a European put option."""
    if T <= 0:
        return np.maximum(K - S, 0)
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)


def cppi_kou_hedged(
    p: Params,
    sigma_hedge: float = 0.25,
    rebal_every: int = 1,
    lam: float = 1.0,
    p_neg: float = 0.6,
    eta_up: float = 25.0,
    eta_dn: float = 10.0,
    seed: int = 42,
) -> np.ndarray:
    """Simulate CPPI under Kou with short-dated OTM put hedging.

    At each rebalancing date, OTM puts are bought to hedge the gap risk
    between rebalancings. Strikes are set so that a jump down to the
    payoff barrier triggers the puts. Put prices are computed under
    Black-Scholes with implied volatility sigma_hedge.

    Parameters
    ----------
    p : Params
        Portfolio parameters.
    sigma_hedge : float, default 0.25
        Implied volatility used to price the hedging puts.
    rebal_every : int, default 1
        Rebalancing frequency in time steps. Also defines the put maturity.
    lam, p_neg, eta_up, eta_dn : float
        Kou jump parameters (see cppi_kou).
    seed : int, default 42
        Random seed.

    Returns
    -------
    V_T : np.ndarray of shape (N_mc,)
        Terminal portfolio values, including final put payoffs.
    """

    rng = np.random.default_rng(seed)
    comp = _kou_compensator(lam, p_neg, eta_up, eta_dn)
    V = np.full(p.N_mc, p.V0)
    e = np.zeros(p.N_mc)
    cash = np.zeros(p.N_mc)
    S = np.full(p.N_mc, p.V0)
    q_puts = np.zeros(p.N_mc)
    K_put = np.zeros(p.N_mc)
    T_put = rebal_every / p.steps_per_year

    for k in range(p.M):
        t = k * p.dt
        if k % rebal_every == 0:
            if k > 0:
                payoff = q_puts * np.maximum(K_put - S, 0.0)
                V = V + payoff
            F = p.G * np.exp(-p.r * (p.T - t))
            C = np.maximum(V - F, 0.0)
            e_new = np.minimum(p.m * C, V)

            K_new = (1 - 1 / p.m) * np.exp(p.r * T_put) * S
            q_new = e_new / S
            put_price = _bs_put_price(S, K_new, T_put, p.r, sigma_hedge)
            put_cost = q_new * put_price

            tc_cost = p.tc * np.abs(e_new - e)
            V = V - tc_cost - put_cost
            e = e_new
            cash = V - e
            q_puts = q_new
            K_put = K_new

        Z = rng.standard_normal(p.N_mc)
        logret = (p.mu - 0.5 * p.sigma**2 - comp) * p.dt + p.sigma * np.sqrt(p.dt) * Z
        Nj = rng.poisson(lam * p.dt, size=p.N_mc)
        mask = Nj > 0
        if mask.any():
            nj = int(mask.sum())
            p_up = 1.0 - p_neg
            up = rng.random(nj) < p_up
            J = np.where(
                up, rng.exponential(1 / eta_up, nj), -rng.exponential(1 / eta_dn, nj)
            )
            logret[mask] += J

        S_ret = np.exp(logret)
        S = S * S_ret
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        e = e * S_ret
        cash = cash * np.exp(p.r * p.dt)

    payoff = q_puts * np.maximum(K_put - S, 0.0)
    V = V + payoff
    return V


def one_path_bs(p: Params, seed: int = 2):
    """Simulate a single CPPI trajectory under BS for visualization.

    Returns
    -------
    t_grid : np.ndarray
        Time grid of length M+1.
    V, F, E, S : np.ndarray
        Portfolio value, floor, risky exposure, underlying price along the path.
    """
    rng = np.random.default_rng(seed)
    M = p.M
    t_grid = np.linspace(0, p.T, M + 1)
    V = p.V0
    e = 0.0
    cash = 0.0
    S = p.V0
    Vs = [V]
    Fs = [p.G * np.exp(-p.r * p.T)]
    Es = [0.0]
    Ss = [p.V0]

    for k in range(M):
        t = k * p.dt
        F = p.G * np.exp(-p.r * (p.T - t))
        C = max(V - F, 0.0)
        e = min(p.m * C, V)
        cash = V - e
        Z = rng.standard_normal()
        S_ret = np.exp((p.mu - 0.5 * p.sigma**2) * p.dt + p.sigma * np.sqrt(p.dt) * Z)
        V = e * S_ret + cash * np.exp(p.r * p.dt)
        S = S * S_ret
        Vs.append(V)
        Fs.append(p.G * np.exp(-p.r * (p.T - (t + p.dt))))
        Es.append(e * S_ret)
        Ss.append(S)
    return t_grid, np.array(Vs), np.array(Fs), np.array(Es), np.array(Ss)