"""
copulagarch.py -- Copula-GARCH: time-varying volatility + a dependence copula.

Backs the notebooks in  "Copula-GARCH: Time-Varying Volatility Meets Dependence".  It fuses two engines already
in the collection:
  * garch_scratch.py  -- univariate GARCH(1,1)-t for each asset's VOLATILITY,
  * copulas.py        -- a copula for the DEPENDENCE of the standardized residuals.

WHY
---
A static copula fit to raw returns confounds two very different things: the fact
that markets go through calm and turbulent regimes (volatility clustering) and
the fact that assets move together (dependence).  The copula-GARCH decomposition
separates them:
  1. Fit a GARCH marginal to each asset -> the conditional volatility sigma_{j,t}
     and the standardized residual  z_{j,t} = (r_{j,t} - mu_j) / sigma_{j,t}.
     The residuals are close to i.i.d. -- GARCH removes the clustering, so the
     z's are the genuine "invariants" of the horizon project.
  2. Fit a copula to the residuals' dependence.
Re-assembling  r_{j,t} = mu_j + sigma_{j,t} z_{j,t}  with z drawn from the copula
gives a portfolio model whose risk BREATHES with the market: a one-day-ahead VaR
that rises in turbulent regimes and falls in calm ones -- which a static VaR
cannot do, and which is exactly what its clustered backtest violations demanded.
"""

import numpy as np
import garch_scratch as gs
import copulas as cp
from scipy.stats import t as student_t
from scipy.optimize import minimize


# --------------------------------------------------------------------------- #
#  Step 1 -- GARCH marginals                                                    #
# --------------------------------------------------------------------------- #

def fit_marginal(r):
    """Fit a GARCH(1,1) with Student-t innovations to one return series, reusing
    garch_scratch's log-likelihood and variance recursion but a self-contained,
    robust optimiser (avoids the engine's numerical-Hessian step). Returns the
    parameters, the conditional-volatility path, and the standardized residuals z
    (which should be ~ i.i.d.)."""
    r = np.asarray(r, float); mu = r.mean(); rc = r - mu
    r2, s0 = gs.prepare(rc)
    nll = lambda th: -gs.garch_loglik(th, r2, s0, gjr=False)         # 4-vector -> Student-t
    init = np.array([0.05 * np.var(rc), 0.05, 0.90, 7.0])
    best = None
    for start in (init, [0.1 * np.var(rc), 0.10, 0.85, 5.0], [0.02 * np.var(rc), 0.03, 0.95, 10.0]):
        res = minimize(nll, np.asarray(start, float), method="Nelder-Mead",
                       options={"maxiter": 8000, "xatol": 1e-7, "fatol": 1e-7})
        if np.isfinite(res.fun) and (best is None or res.fun < best.fun):
            best = res
    mle = best.x
    s2 = gs.garch_var(mle, r2, s0, gjr=False)
    return dict(theta=mle, mu=mu, sigma=np.sqrt(s2), sigma2=s2,
                z=rc / np.sqrt(s2), nu=float(mle[-1]), se=None,
                persistence=float(mle[1] + mle[2]))


def fit_all(X):
    return [fit_marginal(X[:, j]) for j in range(X.shape[1])]


def std_resid_matrix(marg):
    """Matrix of standardized residuals (T x N)."""
    return np.column_stack([m["z"] for m in marg])


# --------------------------------------------------------------------------- #
#  Step 2 -- copula on the residuals, and re-assembly                           #
# --------------------------------------------------------------------------- #

def residual_pool(marg, R_t, nu_c, n, rng, semiparametric=True):
    """Draw n standardized-residual vectors with the fitted t-copula dependence
    (R_t, nu_c). SEMIPARAMETRIC (default): invert each uniform draw through the
    asset's EMPIRICAL residual distribution -- the standard copula-GARCH choice,
    which calibrates the tails correctly. Set semiparametric=False to use a
    parametric standardized-Student-t marginal instead."""
    N = len(marg)
    U = cp.sim_t_copula(R_t, nu_c, n, rng)
    if semiparametric:
        return np.column_stack([np.quantile(marg[j]["z"], np.clip(U[:, j], 1e-4, 1 - 1e-4))
                                for j in range(N)])
    Z = np.empty((n, N))
    for j in range(N):
        nu = marg[j]["nu"]
        Z[:, j] = student_t.ppf(np.clip(U[:, j], 1e-6, 1 - 1e-6), nu) * np.sqrt((nu - 2) / nu)
    return Z


def dynamic_var(marg, w, Z_pool, alpha=0.01):
    """One-day-ahead copula-GARCH VaR for every day, using each day's conditional
    volatilities to scale a fixed pool of copula residuals. Returns VaR_t (>=0)."""
    w = np.asarray(w, float); N = len(marg)
    sig = np.column_stack([m["sigma"] for m in marg])       # (T, N) conditional vols
    mu = np.array([m["mu"] for m in marg])
    T = sig.shape[0]; var = np.empty(T)
    for t in range(T):
        port = (Z_pool * (w * sig[t])) .sum(1) + float(w @ mu)   # portfolio return sims for day t
        var[t] = -np.quantile(port, alpha)
    return var


def static_var_rolling(r_port, window, alpha=0.01):
    """Rolling historical VaR of the portfolio return series (the static baseline)."""
    var = np.full(len(r_port), np.nan)
    for t in range(window, len(r_port)):
        var[t] = -np.quantile(r_port[t - window:t], alpha)
    return var