import numpy as np  # 导入numpy用于数值计算
import matplotlib.pyplot as plt  # 导入绑图库
from scipy import stats  # 导入统计模块

plt.rcParams['font.sans-serif'] = ['Source Han Serif SC']  # 设置中文字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

np.random.seed(42)  # 设置随机种子保证可重现
num_tests = 1000  # 假设检验次数（1000个因子）
null_test_statistics = np.random.standard_normal(num_tests)  # 从标准正态分布生成检验统计量（所有H0为真）
null_p_values = 2 * (1 - stats.norm.cdf(np.abs(null_test_statistics)))  # 计算双侧p值

num_false_positives = np.sum(null_p_values < 0.05)  # 统计p<0.05的假阳性个数

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))  # 创建1行2列子图

# 左图：p值直方图
axes[0].hist(null_p_values, bins=50, color='#008080', alpha=0.7, edgecolor='white')  # 绘制p值频率直方图
axes[0].axvline(x=0.05, color='#E3120B', linestyle='--', linewidth=2, label='α = 0.05')  # 标注显著性阈值
axes[0].set_xlabel('p 值')  # x轴标签
axes[0].set_ylabel('频次')  # y轴标签
axes[0].set_title(f'p值分布（假阳性: {num_false_positives} 个）')  # 标题（含假阳性数）
axes[0].legend()  # 显示图例

# 右图：-log10(p)散点图
axes[1].scatter(range(num_tests), -np.log10(null_p_values), s=5, alpha=0.5, color='gray')  # 绘制-log10(p)散点
axes[1].axhline(y=-np.log10(0.05), color='#E3120B', linestyle='--', linewidth=2, label='α = 0.05')  # 显著性阈值线
axes[1].set_xlabel('因子编号')  # x轴标签
axes[1].set_ylabel('$-\\log_{10}(p)$')  # y轴标签
axes[1].set_title('所有因子都无效，但假阳性大量出现')  # 标题
axes[1].legend()  # 显示图例

plt.tight_layout()  # 自动调整间距
plt.show()  # 显示图表

from scipy import stats  # 导入统计库
from statsmodels.stats.multitest import multipletests  # 导入多重检验校正工具

# 对每只股票执行单侧 t 检验 (H1: μ > 0)
p_value_list = []  # 存储原始p值
t_stat_list = []  # 存储t统计量
for stock_code in valid_stock_list[:500]:  # 取前500只股票（控制计算时间）
    stock_returns = stock_price_2020[stock_price_2020['order_book_id'] == stock_code]['daily_return'].dropna()  # 提取收益率
    t_statistic, two_sided_p = stats.ttest_1samp(stock_returns, 0)  # 单样本t检验
    one_sided_p = two_sided_p / 2 if t_statistic > 0 else 1 - two_sided_p / 2  # 转为单侧p值
    p_value_list.append(one_sided_p)  # 记录p值
    t_stat_list.append(t_statistic)  # 记录t统计量

p_values_array = np.array(p_value_list)  # 转为numpy数组

# Bonferroni 校正
bonferroni_reject, bonferroni_p_adjusted, _, _ = multipletests(p_values_array, alpha=0.05, method='bonferroni')  # Bonferroni校正

fig, axes = plt.subplots(1, 2, figsize=(13, 5))  # 创建1行2列子图

# 左图：未校正
is_significant_raw = p_values_array < 0.05  # 未校正下的显著性判断
axes[0].scatter(range(len(p_values_array)), -np.log10(p_values_array),  # 绘制-log10(p)散点
    c=['#E3120B' if s else 'gray' for s in is_significant_raw], s=8, alpha=0.6)  # 显著=红色
axes[0].axhline(y=-np.log10(0.05), color='blue', linestyle='--', label='α=0.05')  # 原始阈值线
axes[0].set_title(f'未校正（显著: {sum(is_significant_raw)} 只）')  # 标题
axes[0].set_xlabel('股票编号')  # x轴
axes[0].set_ylabel('$-\\log_{10}(p)$')  # y轴
axes[0].legend()  # 图例

# 右图：Bonferroni 校正后
axes[1].scatter(range(len(p_values_array)), -np.log10(p_values_array),  # 同样的-log10(p)
    c=['#E3120B' if s else 'gray' for s in bonferroni_reject], s=8, alpha=0.6)  # Bonferroni显著=红色
bonferroni_threshold = 0.05 / len(p_values_array)  # 计算Bonferroni校正阈值
axes[1].axhline(y=-np.log10(bonferroni_threshold), color='blue', linestyle='--', label=f'α/m={bonferroni_threshold:.2e}')  # 校正阈值线
axes[1].set_title(f'Bonferroni 校正（显著: {sum(bonferroni_reject)} 只）')  # 标题
axes[1].set_xlabel('股票编号')  # x轴
axes[1].set_ylabel('$-\\log_{10}(p)$')  # y轴
axes[1].legend()  # 图例

plt.tight_layout()  # 自动调整间距
plt.show()  # 显示图表

fig, axes = plt.subplots(1, 2, figsize=(13, 5))  # 创建1行2列子图

# 左图：BH 自适应阈值线
sorted_indices = np.argsort(p_arr)  # 按p值从小到大排序的索引
sorted_p = p_arr[sorted_indices]  # 排序后的p值
m_total = len(p_arr)  # 总检验次数
bh_threshold_line = np.arange(1, m_total + 1) / m_total * 0.05  # BH阈值线: k/(m*q)，线性递增

axes[0].scatter(range(m_total), sorted_p, s=3, c='gray', alpha=0.5, label='排序后的p值')  # 排序p值散点图
axes[0].plot(range(m_total), bh_threshold_line, color='#008080', linewidth=2, label='BH 阈值线')  # BH阈值线
axes[0].set_xlabel('排序位置 k')  # x轴
axes[0].set_ylabel('p 值')  # y轴
axes[0].set_title('BH 自适应阈值线')  # 标题
axes[0].set_ylim(0, 0.1)  # 限制y轴范围便于观察
axes[0].legend()  # 图例

# 右图：火山图
for i in range(len(t_arr)):  # 遍历每个检验
    color = 'gray'  # 默认灰色（不显著）
    if bh_reject[i] and t_arr[i] > 0:  # BH显著且t为正
        color = '#E3120B'  # 红色（显著正效应）
    elif bh_reject[i] and t_arr[i] < 0:  # BH显著且t为负
        color = '#008080'  # 青色（显著负效应）
    axes[1].scatter(t_arr[i], -np.log10(p_arr[i]), s=6, c=color, alpha=0.5)  # 绘制散点

axes[1].set_xlabel('t 统计量')  # x轴
axes[1].set_ylabel('$-\\log_{10}(p)$')  # y轴
axes[1].set_title('火山图（BH 校正）')  # 标题

plt.tight_layout()  # 自动调整间距
plt.show()  # 显示图表

# 获取海康威视(002415.XSHE)的股价数据
hikvision_data = stock_price_data[  # 筛选海康威视
    (stock_price_data['order_book_id'] == '002415.XSHE') &  # 股票代码
    (stock_price_data['date'] >= '2020-01-01')  # 2020年以来
].copy().sort_values('date')  # 按日期排序

hikvision_data['daily_return'] = hikvision_data['close'].pct_change()  # 计算日收益率

# 测试 50 条均线策略 (MA5 到 MA54)
ma_periods = range(5, 55)  # MA周期从5到54
strategy_p_values = []  # 存储各策略的p值

for ma_period in ma_periods:  # 遍历每个MA周期
    hikvision_data[f'MA_{ma_period}'] = hikvision_data['close'].rolling(ma_period).mean()  # 计算移动平均
    signal_mask = hikvision_data['close'] > hikvision_data[f'MA_{ma_period}']  # 价格在均线上方时持仓
    strategy_returns = hikvision_data.loc[signal_mask, 'daily_return'].dropna()  # 持仓期间的收益率
    if len(strategy_returns) > 30:  # 样本量充足时
        t_stat, p_val = stats.ttest_1samp(strategy_returns, 0)  # t检验：均值是否>0
        strategy_p_values.append(p_val / 2 if t_stat > 0 else 1)  # 转为单侧p值
    else:  # 样本量不足
        strategy_p_values.append(1.0)  # 设为不显著

strategy_p_arr = np.array(strategy_p_values)  # 转numpy数组
bonf_reject_ma, _, _, _ = multipletests(strategy_p_arr, alpha=0.05, method='bonferroni')  # Bonferroni校正
bh_reject_ma, _, _, _ = multipletests(strategy_p_arr, alpha=0.05, method='fdr_bh')  # BH校正

fig, ax = plt.subplots(figsize=(11, 5))  # 创建图表
x_pos = list(ma_periods)  # x轴位置
ax.bar(x_pos, -np.log10(strategy_p_arr), color=[  # 绘制-log10(p)柱状图
    '#E3120B' if strategy_p_arr[i] < 0.05 else 'gray' for i in range(len(strategy_p_arr))  # 未校正显著=红色
], alpha=0.7)  # 透明度0.7
ax.axhline(y=-np.log10(0.05), color='blue', linestyle='--', label='未校正 α=0.05')  # 未校正阈值
ax.axhline(y=-np.log10(0.05/len(strategy_p_arr)), color='#E3120B', linestyle='-', label='Bonferroni')  # Bonferroni阈值
ax.set_xlabel('移动平均周期 (MA)')  # x轴
ax.set_ylabel('$-\\log_{10}(p)$')  # y轴
ax.set_title(f'海康威视50条MA策略: 未校正显著 {sum(strategy_p_arr < 0.05)} 条 → 校正后 {sum(bonf_reject_ma)} 条')  # 标题
ax.legend()  # 图例
plt.tight_layout()  # 自动调整间距
plt.show()  # 显示图表

# 三种校正
industry_raw_sig = industry_p_arr < 0.05  # 未校正显著判断
industry_bonf_rej, _, _, _ = multipletests(industry_p_arr, alpha=0.05, method='bonferroni')  # Bonferroni
industry_bh_rej, _, _, _ = multipletests(industry_p_arr, alpha=0.05, method='fdr_bh')  # BH

fig, axes = plt.subplots(1, 3, figsize=(15, 5))  # 1行3列子图
method_configs = [  # 三种方法的配置
    ('未校正 (α=0.05)', industry_raw_sig, int(sum(industry_raw_sig))),  # 未校正配置
    ('Bonferroni', industry_bonf_rej, int(sum(industry_bonf_rej))),  # Bonferroni配置
    ('BH (q=0.05)', industry_bh_rej, int(sum(industry_bh_rej)))  # BH配置
]  # 三种方法

for ax_idx, (method_name, reject_mask, count) in enumerate(method_configs):  # 遍历三种方法
    colors = ['#E3120B' if r else 'gray' for r in reject_mask]  # 显著=红色，不显著=灰色
    axes[ax_idx].barh(range(len(industry_p_arr)), -np.log10(industry_p_arr), color=colors, alpha=0.7)  # 水平柱状图
    axes[ax_idx].set_title(f'{method_name}（显著: {count}）')  # 标题含显著数量
    axes[ax_idx].set_xlabel('$-\\log_{10}(p)$')  # x轴标签

plt.tight_layout()  # 自动调整间距
plt.show()  # 显示图表