"""M08 教学图：零溢价下的最大夏普，以及同一因子两段窗口的 t 值。

左图是独立正态夏普的模拟，不是某次回测。右图复用 M07 的月频解析，
只读已缓存的 French 文件。t 值按月收益独立同分布，没有做 Newey–West。
"""

from __future__ import annotations

import importlib.util
import math
import zipfile
from io import BytesIO
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["font.sans-serif"] = ["Microsoft YaHei", "SimHei", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False

OUT = Path(__file__).resolve().parent / "multiple_testing.png"
M07 = Path(__file__).resolve().parents[1] / "M07" / "plot_ff5_windows.py"


def load_m07():
    spec = importlib.util.spec_from_file_location("m07plot", M07)
    module = importlib.util.module_from_spec(spec)
    spec.loader.exec_module(module)
    return module


def norm_cdf(z: float) -> float:
    return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))


def norm_ppf(p: float) -> float:
    # 默认环境没有 scipy。二分法够用来核对 Bonferroni 的正态临界值。
    lo, hi = -8.0, 8.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if norm_cdf(mid) < p:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


def bonferroni_z(m: int, alpha: float = 0.05) -> float:
    # 双侧检验：m 次里至少一次误拒的概率控制在 alpha。
    return norm_ppf(1.0 - alpha / (2.0 * m))


def holm_bh(p_values: list[float], alpha: float = 0.05) -> None:
    m = len(p_values)
    order = sorted(range(m), key=lambda i: p_values[i])
    print("bonferroni_cut", f"{alpha / m:.6f}")
    holm_stop = None
    for rank, idx in enumerate(order, start=1):
        cut = alpha / (m - rank + 1)
        ok = p_values[idx] <= cut
        print(
            "holm", rank, f"{p_values[idx]:.4f}",
            "cut", f"{cut:.6f}", "pass", int(ok),
        )
        if not ok and holm_stop is None:
            holm_stop = rank
            break
    # BH 看排序后最大的 k，使 p_(k) <= alpha * k / m。后面更大的 p 不会回头。
    best = 0
    for rank, idx in enumerate(order, start=1):
        cut = alpha * rank / m
        if p_values[idx] <= cut:
            best = rank
        print("bh", rank, f"{p_values[idx]:.4f}", "cut", f"{cut:.6f}")
    print("bh_reject_first", best)


def factor_t(m07) -> dict[str, tuple[float, float]]:
    raw = zipfile.ZipFile(BytesIO(m07.fetch_ff5())).read(
        "F-F_Research_Data_5_Factors_2x3.csv"
    )
    ff5 = m07.monthly_block(raw.decode("latin1"), ",Mkt-RF")
    mom_raw = zipfile.ZipFile(m07.MOM_ZIP).read("F-F_Momentum_Factor.csv")
    mom = m07.monthly_block(mom_raw.decode("latin1"), ",Mom")
    data = ff5.join(mom.rename(columns={mom.columns[0]: "Mom"}), how="inner")
    out: dict[str, tuple[float, float]] = {}
    late = data.loc["2016-01-01":"2025-12-01"]
    for name in m07.NAMES:
        full = m07.stats(data[name])["t"]
        recent = m07.stats(late[name])["t"]
        out[name] = (full, recent)
        print(name, "t_full", f"{full:.2f}", "t_2016_2025", f"{recent:.2f}")
    return out


def deflated_sharpe(
    sr: float,
    sr_std: float,
    n_trials: int,
    n_obs: int,
    skew: float,
    kurtosis: float,
) -> dict[str, float]:
    """Bailey 与 Lopez de Prado (2014) 的 DSR。夏普必须和收益同一频率。

    年化夏普直接代进来，会把门槛和标准误乘到另一套尺度上。Euler 常数用正文给出的
    0.5772156649，不再加更长的尾数，这样打印才能和手算对上。
    """
    gamma = 0.5772156649
    # 独立试验次数 M 下，标准正态最大值的期望。乘的是尝试集合的截面标准差，不是这一条的标准误。
    z_m = norm_ppf(1.0 - 1.0 / n_trials)
    z_me = norm_ppf(1.0 - 1.0 / (n_trials * math.e))
    z0 = (1.0 - gamma) * z_m + gamma * z_me
    sr0 = sr_std * z0
    # Pearson 峰度。正态是 3，所以 (kurtosis-1)/4 在正态下回到 SR^2/2。超额峰度不能直接代。
    scale = 1.0 - skew * sr + (kurtosis - 1.0) / 4.0 * sr * sr
    se = math.sqrt(scale / (n_obs - 1))
    z_stat = (sr - sr0) / se
    return {
        "z_m": z_m,
        "z_me": z_me,
        "z0": z0,
        "sr0": sr0,
        "scale": scale,
        "se": se,
        "z_stat": z_stat,
        "dsr": norm_cdf(z_stat),
        "psr": norm_cdf(sr / se),
    }


def combinations(indices: list[int], k: int) -> list[tuple[int, ...]]:
    """从 indices 里取 k 个的全部组合。CSCV 用 k=S/2，互补的两组都要各算一次。"""
    if k == 0:
        return [()]
    picked: list[tuple[int, ...]] = []
    for pos, item in enumerate(indices):
        for tail in combinations(indices[pos + 1 :], k - 1):
            picked.append((item, *tail))
    return picked


def probability_of_backtest_overfitting(block_returns: np.ndarray) -> dict[str, float]:
    """样本内冠军落到样本外后半区的比例。得分是选中分块之和，等长时与均值排序相同。

    名次从 1 起，1 最好。后半区指名次大于列数的一半。并列会让名次跟着排序稳定性走，
    教学矩阵里四列得分都不要并列。
    """
    n_blocks, n_strat = block_returns.shape
    if n_blocks % 2 != 0:
        raise ValueError("CSCV 的分块数要是偶数，才能切成同样长的两半")
    splits = combinations(list(range(n_blocks)), n_blocks // 2)
    n_bottom = 0
    for split in splits:
        holdout = tuple(i for i in range(n_blocks) if i not in split)
        in_sample = block_returns[list(split)].sum(axis=0)
        out_sample = block_returns[list(holdout)].sum(axis=0)
        champ = int(np.argmax(in_sample))
        order = np.argsort(-out_sample, kind="mergesort")
        rank = int(np.where(order == champ)[0][0]) + 1
        bottom = int(rank > n_strat / 2)
        n_bottom += bottom
        print(
            "split",
            ",".join(str(i + 1) for i in split),
            "is",
            [int(v) for v in in_sample],
            "champ",
            champ,
            "oos",
            [int(v) for v in out_sample],
            "rank",
            rank,
            "bottom",
            bottom,
        )
    return {
        "n_bottom": float(n_bottom),
        "n_splits": float(len(splits)),
        "pbo": n_bottom / len(splits),
    }


def selection_bias_demo() -> None:
    """月频手算。截面标准差 0.12 是本例给定的尝试集合离散程度，不是上一段模拟的年化标准误。"""
    monthly = 1.20 / math.sqrt(12.0)
    result = deflated_sharpe(
        sr=monthly,
        sr_std=0.12,
        n_trials=100,
        n_obs=120,
        skew=-0.5,
        kurtosis=5.0,
    )
    print("sr_monthly", f"{monthly:.6f}")
    print("z_m", f"{result['z_m']:.6f}")
    print("z_me", f"{result['z_me']:.6f}")
    print("Z0", f"{result['z0']:.6f}")
    print("SR0", f"{result['sr0']:.6f}")
    print("se_num", f"{result['scale']:.6f}")
    print("SE", f"{result['se']:.6f}")
    print("z_stat", f"{result['z_stat']:.6f}")
    print("PSR", f"{result['psr']:.6f}")
    print("DSR", f"{result['dsr']:.6f}")
    try:
        from scipy.stats import norm

        # 对照二分法。未安装时不改上面的六位。
        z_ref = float(norm.ppf(1.0 - 1.0 / 100.0))
        print("scipy_ppf", f"{z_ref:.6f}")
    except ImportError:
        print("scipy_norm", "skipped")
    # 四块百分比收益。列是四条策略，行是四段等长时间。数字是构造的，不是行情。
    blocks = np.array(
        [
            [1.0, -3.0, 1.0, -1.0],
            [7.0, 7.0, -2.0, 2.0],
            [0.0, 2.0, 8.0, 6.0],
            [-6.0, 4.0, 6.0, -2.0],
        ],
        dtype=float,
    )
    summary = probability_of_backtest_overfitting(blocks)
    print("n_bottom", int(summary["n_bottom"]))
    print("n_splits", int(summary["n_splits"]))
    print("PBO", f"{summary['pbo']:.6f}")


def snooping() -> tuple[np.ndarray, float]:
    rng = np.random.default_rng(7)
    years = 10.0
    se = 1.0 / math.sqrt(years)
    # 零溢价、独立、年化夏普的标准误约 1/√年数。这里模拟的是估计量，不是日收益路径。
    draws = rng.normal(0.0, se, size=(20_000, 200))
    maximum = draws.max(axis=1)
    print("se", f"{se:.4f}")
    print("mean_max", f"{float(maximum.mean()):.3f}")
    print("p95_max", f"{float(np.quantile(maximum, 0.95)):.3f}")
    print("share_any_above_1", f"{float((maximum > 1.0).mean()):.3f}")
    print("share_slots_above_0.5", f"{float((draws > 0.5).mean()):.3f}")
    print("expected_count_above_0.5", f"{float((draws > 0.5).sum(axis=1).mean()):.2f}")
    lead = se * math.sqrt(2.0 * math.log(200))
    print("sqrt_2ln_rule", f"{lead:.3f}")
    # 一次有限样本：日收益均值真为 0。夏普用样本标准差，rf 已经是 0。
    daily = rng.normal(0.0, 0.01, size=(200, 10 * 252))
    mean = daily.mean(axis=1)
    std = daily.std(axis=1, ddof=1)
    sharpe = mean / std * math.sqrt(252.0)
    print(
        "one_sample_max", f"{float(sharpe.max()):.3f}",
        "count_above_1", int((sharpe > 1.0).sum()),
        "count_above_0.5", int((sharpe > 0.5).sum()),
    )
    return maximum, se


def main() -> None:
    selection_bias_demo()
    ps = [0.001, 0.004, 0.012, 0.02, 0.03] + [0.11] * 15
    holm_bh(ps)
    for m in (1, 20, 100, 316):
        z = bonferroni_z(m)
        extreme = math.sqrt(2.0 * math.log(m)) if m > 1 else 0.0
        print("m", m, "bonf_z", f"{z:.3f}", "sqrt2ln", f"{extreme:.3f}")
    for years in (5, 10, 20, 50):
        print(
            "years", years,
            "sr_at_t2", f"{1.96 / math.sqrt(years):.3f}",
            "sr_at_t3", f"{3.0 / math.sqrt(years):.3f}",
        )
    m07 = load_m07()
    ts = factor_t(m07)
    maximum, _se = snooping()

    labels = ["市场", "规模", "价值", "盈利", "投资", "动量"]
    full = [ts[name][0] for name in m07.NAMES]
    late = [ts[name][1] for name in m07.NAMES]
    fig, axes = plt.subplots(1, 2, figsize=(11.2, 4.8))
    ax = axes[0]
    ax.hist(maximum, bins=40, color="#1f4e79", alpha=0.85)
    ax.axvline(1.0, color="#c1121f", lw=1.2, label="夏普 = 1")
    ax.axvline(float(maximum.mean()), color="#2a9d8f", lw=1.2, label="最大夏普的均值")
    ax.set_xlabel("200 条零溢价策略里的最大年化夏普")
    ax.set_ylabel("模拟次数")
    ax.legend(frameon=False, fontsize=8)
    ax.grid(True, axis="y", alpha=0.3)
    ax = axes[1]
    xpos = np.arange(len(labels))
    ax.bar(xpos - 0.18, full, width=0.36, color="#1f4e79", label="1963–2026.08")
    ax.bar(xpos + 0.18, late, width=0.36, color="#c1121f", label="2016–2025")
    ax.axhline(1.96, color="#666666", lw=0.8, ls="--")
    ax.axhline(3.0, color="#333333", lw=0.8, ls=":")
    ax.axhline(0.0, color="#333333", lw=0.6)
    ax.set_xticks(xpos)
    ax.set_xticklabels(labels, fontsize=8)
    ax.set_ylabel("iid t 值")
    ax.legend(frameon=False, fontsize=8)
    ax.grid(True, axis="y", alpha=0.3)
    fig.suptitle("M08 多次尝试下的最大夏普，以及同一因子的两段 t 值", fontsize=12)
    fig.tight_layout()
    fig.savefig(OUT, dpi=140, bbox_inches="tight")
    print("wrote", OUT)


if __name__ == "__main__":
    main()