多重检验夏普折扣：用 Bonferroni/BHY 给数据挖掘出来的夏普打折

import numpy as np
from scipy import stats

def sharpe_to_tstat(sharpe_annual, years):
    """年化夏普率 -> t 统计量"""
    return sharpe_annual * np.sqrt(years)

def tstat_to_pvalue(t, two_sided=True):
    """t 统计量 -> p 值（用正态近似，大样本足够）"""
    p = stats.norm.sf(abs(t))
    return 2 * p if two_sided else p

sr, years = 1.0, 10
t = sharpe_to_tstat(sr, years)
p = tstat_to_pvalue(t)
print(f"夏普 {sr} / {years}年 => t = {t:.2f}, p = {p:.4f}")

def haircut_sharpe(sharpe, years, M, method='holm', alpha=0.05):
    """
    对数据挖掘出来的夏普率打折。
    返回：原始 p、校正后 p、是否仍显著、折扣后的夏普率
    """
    t = sharpe * np.sqrt(years)
    p_raw = 2 * stats.norm.sf(abs(t))

    if method == 'bonferroni':
        p_adj = min(p_raw * M, 1.0)
    elif method == 'holm':
        # 单个策略的 Holm 上界近似等于 Bonferroni（最悲观位置）
        p_adj = min(p_raw * M, 1.0)
    elif method == 'bhy':
        # BHY 带相关性修正项 c(M) = sum(1/i)
        cM = np.sum(1.0 / np.arange(1, M + 1))
        p_adj = min(p_raw * M * cM / 1.0, 1.0)
    else:
        raise ValueError(method)

    # 由校正后的 p 反推"折扣后"的 t 和夏普
    t_adj = stats.norm.isf(p_adj / 2)
    sr_adj = t_adj / np.sqrt(years)
    haircut_pct = (1 - sr_adj / sharpe) * 100 if sharpe else 0

    return {
        'p_raw': p_raw,
        'p_adj': p_adj,
        'significant': p_adj < alpha,
        'sharpe_raw': sharpe,
        'sharpe_adj': max(sr_adj, 0),
        'haircut_pct': haircut_pct,
    }

# 那个"闪闪发光"的夏普 2.3，其实是 500 次搜索的产物
res = haircut_sharpe(sharpe=2.3, years=5, M=500, method='bhy')
for k, v in res.items():
    print(f"{k:>12}: {v}")

np.random.seed(7)
years = 8

# 45 个纯噪声策略（真夏普=0）+ 5 个真 alpha（真夏普≈1.2）
true_sr = np.r_[np.zeros(45), np.full(5, 1.2)]
# 观测夏普 = 真夏普 + 抽样噪声（噪声标准差≈1/√years）
obs_sr = true_sr + np.random.normal(0, 1/np.sqrt(years), size=50)
M = len(obs_sr)

t_all = obs_sr * np.sqrt(years)
p_all = 2 * stats.norm.sf(np.abs(t_all))

# BHY / FDR 校正
order = np.argsort(p_all)
cM = np.sum(1.0 / np.arange(1, M + 1))
p_sorted = p_all[order]
# BHY 阈值：p_(k) <= (k / (M*cM)) * alpha
alpha = 0.10
thresh = (np.arange(1, M+1) / (M * cM)) * alpha
passed = p_sorted <= thresh
k_max = np.max(np.where(passed)[0]) + 1 if passed.any() else 0

selected = order[:k_max]
print(f"校正前 p<0.05 的策略数: {(p_all < 0.05).sum()}")
print(f"BHY 校正后入选策略数:   {len(selected)}")
print(f"其中真 alpha (index>=45) 命中: "
      f"{(selected >= 45).sum()} / 5")
print(f"其中假阳性 (index<45):        "
      f"{(selected < 45).sum()}")