Multiple testing and backtest-overfitting diagnostics

"""M04 exercise: multiple testing and backtest-overfitting diagnostics (notes §4.6).

Implement every function, then run from the repo root:

    cd exercises && .venv/bin/python -m pytest modules/M04-systematic-research-swing/tests/test_m04_multiple_testing.py

p-values are two-sided unless stated. Every adjust function returns
(reject: bool array, p_adj: float array) in the ORIGINAL order of the input.
The tests compare your adjusted p-values with statsmodels.stats.multitest.multipletests,
but do not call statsmodels inside your functions.
"""

from __future__ import annotations

from itertools import combinations

import numpy as np
from scipy.stats import norm


def bonferroni(pvals, alpha: float = 0.05):
    """Reject H_i iff m p_i <= alpha. Adjusted p = min(1, m p_i). Controls FWER (Prop. 4.14)."""
    raise NotImplementedError


def holm(pvals, alpha: float = 0.05):
    """Holm (1979) step-down (Prop. 4.15). Sort p_(1) <= ... <= p_(m).

    Adjusted p_(i) = max_{j <= i} min(1, (m - j + 1) p_(j)). Reject where adjusted <= alpha.
    """
    raise NotImplementedError


def benjamini_hochberg(pvals, alpha: float = 0.05):
    """Benjamini-Hochberg (1995) step-up (Prop. 4.16).

    Adjusted p_(i) = min_{j >= i} min(1, m p_(j) / j). Reject where adjusted <= alpha.
    """
    raise NotImplementedError


def benjamini_yekutieli(pvals, alpha: float = 0.05):
    """Benjamini-Yekutieli (2001): BH with m replaced by m c(m), c(m) = sum_{i=1}^m 1/i. Adjusted p capped at 1."""
    raise NotImplementedError


def t_to_p(t) -> np.ndarray:
    """Two-sided normal p-value 2 (1 - Phi(|t|)), vectorised."""
    raise NotImplementedError


def bonferroni_t_threshold(m: int, alpha: float = 0.05) -> float:
    """Normal t cut-off for a two-sided Bonferroni test of m hypotheses at FWER alpha: Phi^{-1}(1 - alpha/(2m))."""
    raise NotImplementedError


def simulate_error_rates(m: int = 100, m0: int = 80, effect: float = 3.0, alpha: float = 0.05, n_sims: int = 2000, seed: int = 0) -> dict:
    """Monte Carlo of FWER, FDR and power for methods 'none', 'bonferroni', 'holm', 'bh'.

    Use rng = np.random.default_rng(seed). In each simulation draw z = standard normal (m,)
    and add `effect` to hypotheses m0..m-1 (the false nulls); hypotheses 0..m0-1 are true nulls.
    p = t_to_p(z). 'none' rejects p <= alpha.
    For each method: fwer = P(V >= 1), fdr = E[V / max(R, 1)], power = E[(R - V) / (m - m0)],
    V = number of true nulls rejected, R = number rejected.
    Return {method: {"fwer": ..., "fdr": ..., "power": ...}}.
    """
    raise NotImplementedError


def pbo_cscv(perf: np.ndarray, n_blocks: int = 16) -> float:
    """Probability of backtest overfitting via CSCV (notes §4.6.6, Algorithm 4.3).

    perf: T x N matrix of per-period returns of N strategy variants (rows = time).
    1. Raise ValueError if n_blocks is odd. Drop the last T mod S rows; split rows into
       S = n_blocks contiguous equal blocks.
    2. For each of the C(S, S/2) choices of S/2 blocks as in-sample (IS), rest out-of-sample (OOS):
       Sharpe (mean/std, ddof=1) of each variant on IS rows and on OOS rows;
       n* = argmax IS Sharpe; omega = (rank of n* among OOS Sharpes, 1 = worst) / (N + 1);
       lambda = ln(omega / (1 - omega)).
    3. Return the fraction of splits with lambda <= 0.
    Hint: precompute per-block sums and sums of squares, then use an indicator matrix.
    """
    raise NotImplementedError


def stationary_bootstrap_indices(T: int, mean_block: float, rng: np.random.Generator) -> np.ndarray:
    """Politis-Romano (1994) stationary bootstrap index path of length T.

    idx_0 uniform on {0..T-1}; for t >= 1, with probability 1/mean_block jump to a new uniform
    index, otherwise idx_t = (idx_{t-1} + 1) mod T.
    """
    raise NotImplementedError


def reality_check(excess: np.ndarray, n_boot: int = 500, mean_block: float = 10.0, seed: int = 0) -> float:
    """White's (2000) Reality Check p-value for H0: max_k E[f_k] <= 0 (notes §4.6.5).

    excess: T x K, f_{t,k} = rule k return minus benchmark return.
    V = sqrt(T) max_k fbar_k; V*_b = sqrt(T) max_k (fbar*_{b,k} - fbar_k) with rows resampled by
    stationary_bootstrap_indices (rng = np.random.default_rng(seed)). Return mean(V*_b >= V).
    """
    raise NotImplementedError