quant-trading/research/lab/report.py

"""Performance statistics, factor IC analysis and overfitting controls."""

from __future__ import annotations

import math

import numpy as np
import pandas as pd
from scipy import stats as sps

import metrics


def stats(equity: pd.Series, turnover: pd.Series | None = None) -> dict:
    eq = equity.dropna()
    r = eq.pct_change().dropna()
    out = {
        "cagr": metrics.annualized_return(eq),
        "vol": metrics.annualized_volatility(r),
        "sharpe": metrics.sharpe_ratio(r),
        "sortino": metrics.sortino_ratio(r),
        "mdd": metrics.max_drawdown(eq),
    }
    out["calmar"] = out["cagr"] / abs(out["mdd"]) if out["mdd"] < 0 else np.inf
    if turnover is not None:
        t = turnover.loc[eq.index[0]:eq.index[-1]]
        years = max((eq.index[-1] - eq.index[0]).days / 365.25, 1e-9)
        out["turnover"] = t.sum() / years          # one-way, x per year
    return out


def fmt(s: dict) -> str:
    base = (f"CAGR {s['cagr']*100:6.1f}%  vol {s['vol']*100:5.1f}%  Sharpe {s['sharpe']:5.2f}  "
            f"Sortino {s['sortino']:5.2f}  MDD {s['mdd']*100:6.1f}%  Calmar {s['calmar']:5.2f}")
    if "turnover" in s:
        base += f"  turn {s['turnover']:5.1f}x"
    return base


def yearly(equity: pd.Series) -> pd.Series:
    e = equity.dropna()
    ye = e.resample("YE").last()
    first = e.iloc[0]
    return ye.pct_change().fillna(ye.iloc[0] / first - 1)


# ---------------------------------------------------------------------------
# Factor analytics
# ---------------------------------------------------------------------------

def forward_returns(panel, horizon: int = 5) -> pd.DataFrame:
    """Return from the open of t+1 to the open of t+1+horizon — the return a
    signal formed at the close of t can actually earn."""
    o = panel.open.where(panel.open > 0).ffill()
    return o.shift(-(1 + horizon)) / o.shift(-1) - 1


def ic_series(score: pd.DataFrame, fwd: pd.DataFrame, member: pd.DataFrame,
              every: int = 5) -> pd.Series:
    """Cross-sectional Spearman rank IC on every ``every``-th date."""
    idx = score.index[::every]
    out = {}
    for d in idx:
        s = score.loc[d].where(member.loc[d])
        f = fwd.loc[d]
        ok = s.notna() & f.notna()
        if ok.sum() < 30:
            continue
        out[d] = sps.spearmanr(s[ok], f[ok]).statistic
    return pd.Series(out)


def top_minus_universe(score: pd.DataFrame, fwd: pd.DataFrame, member: pd.DataFrame,
                       top_n: int = 20, every: int = 5) -> pd.Series:
    """Mean forward return of the top-N names minus the member average."""
    idx = score.index[::every]
    out = {}
    for d in idx:
        s = score.loc[d].where(member.loc[d]).dropna()
        f = fwd.loc[d]
        s = s[f.reindex(s.index).notna()]
        if len(s) < top_n * 3:
            continue
        top = s.nlargest(top_n).index
        out[d] = f[top].mean() - f[s.index].mean()
    return pd.Series(out)


def ic_summary(ic: pd.Series, periods_per_year: float = 52) -> dict:
    ic = ic.dropna()
    if len(ic) < 5:
        return {"ic": np.nan, "icir": np.nan, "t": np.nan, "hit": np.nan, "n": len(ic)}
    return {"ic": ic.mean(), "icir": ic.mean() / ic.std() * math.sqrt(periods_per_year),
            "t": ic.mean() / ic.std() * math.sqrt(len(ic)), "hit": (ic > 0).mean(), "n": len(ic)}


# ---------------------------------------------------------------------------
# Overfitting controls (Bailey & Lopez de Prado 2014)
# ---------------------------------------------------------------------------

def deflated_sharpe(returns: pd.Series, n_trials: int, sr_trials_var: float | None = None,
                    periods: int = 252) -> dict:
    """Probability that the true Sharpe exceeds the best Sharpe expected from
    ``n_trials`` unskilled strategies. Uses per-period Sharpe internally."""
    r = returns.dropna()
    T = len(r)
    sr = r.mean() / r.std()
    g3 = sps.skew(r)
    g4 = sps.kurtosis(r, fisher=False)
    if sr_trials_var is None:
        sr_trials_var = 1.0 / T          # variance of SR estimates under H0
    emc = 0.5772156649
    z1 = sps.norm.ppf(1 - 1.0 / n_trials) if n_trials > 1 else 0.0
    z2 = sps.norm.ppf(1 - 1.0 / (n_trials * math.e)) if n_trials > 1 else 0.0
    sr0 = math.sqrt(sr_trials_var) * ((1 - emc) * z1 + emc * z2)
    denom = math.sqrt(max(1 - g3 * sr + (g4 - 1) / 4 * sr ** 2, 1e-12))
    dsr = sps.norm.cdf((sr - sr0) * math.sqrt(T - 1) / denom)
    return {"sr_annual": sr * math.sqrt(periods), "sr0_annual": sr0 * math.sqrt(periods), "dsr": dsr}