"""Dependence-aware inference: Diebold–Mariano with HAC (overlap-aware lag), cluster bootstrap by trading day
(pooling all assets of a day into one cluster), Holm and Benjamini–Hochberg controls."""
from __future__ import annotations

import numpy as np
from scipy import stats


def hac_variance(d: np.ndarray, lag: int) -> float:
    d = d - d.mean()
    n = len(d)
    v = (d * d).sum() / n
    for k in range(1, lag + 1):
        w = 1 - k / (lag + 1)
        v += 2 * w * (d[k:] * d[:-k]).sum() / n
    return float(max(v, 1e-30))


def diebold_mariano(loss_a: np.ndarray, loss_b: np.ndarray, order: np.ndarray, lag: int) -> dict:
    """DM on the loss differential ordered chronologically (`order` sorts the observations), HAC lag = overlap."""
    d = (loss_a - loss_b)[np.argsort(order, kind="stable")]
    n = len(d)
    v = hac_variance(d, lag)
    t = d.mean() / np.sqrt(v / n)
    return {"mean_diff": float(d.mean()), "t": float(t), "p": float(2 * (1 - stats.norm.cdf(abs(t)))), "n": int(n),
            "lag": int(lag)}


def cluster_bootstrap_mean(x: np.ndarray, clusters: np.ndarray, n_boot: int = 2000, seed: int = 0,
                           alpha: float = 0.05) -> dict:
    """Percentile CI of mean(x) resampling whole clusters (trading days) with replacement."""
    rng = np.random.default_rng(seed)
    uc, inv = np.unique(clusters, return_inverse=True)
    sums = np.bincount(inv, weights=x, minlength=len(uc)); cnt = np.bincount(inv, minlength=len(uc)).astype(float)
    C = len(uc)
    boots = np.empty(n_boot)
    for b in range(n_boot):
        idx = rng.integers(0, C, C)
        boots[b] = sums[idx].sum() / cnt[idx].sum()
    return {"mean": float(x.mean()), "lo": float(np.quantile(boots, alpha / 2)), "hi": float(np.quantile(boots, 1 - alpha / 2)),
            "se": float(boots.std(ddof=1)), "n_clusters": int(C), "n": int(len(x))}


def holm(pvals: list[float]) -> list[float]:
    p = np.asarray(pvals, dtype=float); m = len(p)
    order = np.argsort(p); adj = np.empty(m)
    running = 0.0
    for rank, i in enumerate(order):
        running = max(running, (m - rank) * p[i])
        adj[i] = min(1.0, running)
    return adj.tolist()


def benjamini_hochberg(pvals: list[float]) -> list[float]:
    p = np.asarray(pvals, dtype=float); m = len(p)
    order = np.argsort(p)[::-1]; adj = np.empty(m); running = 1.0
    for rank, i in zip(range(m, 0, -1), order):
        running = min(running, m / rank * p[i]); adj[i] = running
    return adj.tolist()