14 缺失数据处理

# 为“可视化警告：均值填充如何扭曲数据分布”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“可视化警告：均值填充如何扭曲数据分布”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“可视化警告：均值填充如何扭曲数据分布”导入 `seaborn` 并绑定 `sns`，用于编码当前任务的统计分布或分组关系。
import seaborn as sns

# 启用白底网格主题，使原始分布与均值填充分布的形状差异清晰可比。
sns.set_style("whitegrid")
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号

# 固定随机种子，使正态分布与后续 MCAR 遮蔽位置可重复。
np.random.seed(42)
# 抽取均值 100、标准差 20 的 1000 个完整值，作为插补前真值分布。
original_data = np.random.normal(loc=100, scale=20, size=1000)

# 复制完整样本后注入 MCAR 缺失，保留原始数组作为真实值参照。
data_with_missing = original_data.copy()
# 无放回抽取 300 个位置制造 MCAR 缺失，保留 70% 观测值。
missing_indices = np.random.choice(1000, 300, replace=False)
# 在指定索引写入缺失值，构造可与真值对照的插补恢复实验。
data_with_missing[missing_indices] = np.nan

# 仅用已观测值计算均值，作为单次填补的固定替代值。
mean_val = np.nanmean(data_with_missing)
# 提取 `imputed_data` 的插补结果，用于与原始非缺失分布和离散程度比较。
imputed_data = np.nan_to_num(data_with_missing, nan=mean_val)

# 创建共享纵轴的双面板，直接比较原始分布与均值填补后的集中堆积。
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), sharey=True)

# 绘制原始数据
sns.histplot(original_data, bins=30, kde=True, color='skyblue', ax=ax1)
# 固定 `ax1.axvline(np.mean(original_data), color` 的业务取值与顺序，作为“可视化警告：均值填充如何扭曲数据分布”的可复算输入。
ax1.axvline(np.mean(original_data), color='red', linestyle='--', label=f'真实均值: {np.mean(original_data):.2f}')
# 用图题概括“可视化警告：均值填充如何扭曲数据分布”的比较对象与当前计算结果。
ax1.set_title(f'真实数据分布\n(标准差 = {np.std(original_data):.2f})', fontsize=14)
# 为“可视化警告：均值填充如何扭曲数据分布”，将纵轴标为“频数”，明确纵向编码的变量。
ax1.set_ylabel('频数')
# 显示“可视化警告：均值填充如何扭曲数据分布”图例，使颜色或线型与比较对象一一对应。
ax1.legend()

# 绘制填充后数据
sns.histplot(imputed_data, bins=30, kde=True, color='salmon', ax=ax2)
# 固定 `ax2.axvline(mean_val, color` 的业务取值与顺序，作为“可视化警告：均值填充如何扭曲数据分布”的可复算输入。
ax2.axvline(mean_val, color='red', linestyle='--', label=f'填充均值: {mean_val:.2f}')
# 用图题概括“可视化警告：均值填充如何扭曲数据分布”的比较对象与当前计算结果。
ax2.set_title(f'均值填充后的数据分布\n(标准差 = {np.std(imputed_data):.2f})', fontsize=14)
# 显示“可视化警告：均值填充如何扭曲数据分布”图例，使颜色或线型与比较对象一一对应。
ax2.legend()

# 将图题限定到当前 MCAR 生成实验，避免外推为一般偏差定理。
plt.suptitle('当前 MCAR 生成实验：均值填充压缩离散度并制造中心尖峰', fontsize=16)
# 为“可视化警告：均值填充如何扭曲数据分布”，压缩子图留白，避免标题和坐标标签相互遮挡。
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
# 显示原始、缺失与均值填充后的分布，检查填充值在均值处堆积及方差被压缩的失真。
plt.show()

# 为“KNN填充的可视化解释”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np
# 为“KNN填充的可视化解释”导入 `matplotlib.pyplot` 并绑定 `plt`，用于构建当前任务的坐标轴并呈现比较结果。
import matplotlib.pyplot as plt
# 为“KNN填充的可视化解释”，从 `sklearn.neighbors` 导入 `NearestNeighbors`，查找 KNN 插补所依赖的邻近观测。
from sklearn.neighbors import NearestNeighbors
# 为“KNN填充的可视化解释”导入 `seaborn` 并绑定 `sns`，用于编码当前任务的统计分布或分组关系。
import seaborn as sns

# 启用白底网格主题，突出 KNN 距离连线与候选邻居位置。
sns.set_style("whitegrid")
plt.rcParams['font.sans-serif'] = ['Source Han Serif SC'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号

# 固定随机种子，生成两维相关点云作为 KNN 填补的几何机制输入。
np.random.seed(0)
# 在 0—10 方形区域均匀抽取 50 个二维点，作为最近邻几何底图。
X = np.random.rand(50, 2) * 10
# 给第二维加入与第一维相关的扰动，形成 KNN 邻域可见的二维结构。
X[:, 1] = X[:, 1] + 0.5 * X[:, 0] + np.random.randn(50) * 0.5

# 固定第一维已知、第二维缺失的目标点，明确最近邻检索只使用可观测坐标。
point_with_missing = np.array([[5, np.nan]])
# 保留全部其他点作为邻居候选集，待距离搜索选出最近三点。
all_other_points = X

# 固定五个邻居，准备在第一维空间检索局部参照样本。
k = 5
# 建立最近邻搜索器 `nn`，固定邻居数量以演示 KNN 插补机制。
nn = NearestNeighbors(n_neighbors=k)
# 用其余样本的第一特征建立最近邻索引，随后为缺失点寻找距离最小的参照观测。
nn.fit(all_other_points[:, 0].reshape(-1, 1))
# 取得 `distances, indices` 的最近邻结果，用相邻观测解释 KNN 插补值来源。
distances, indices = nn.kneighbors(point_with_missing[:, 0].reshape(-1, 1))

# 取得 `neighbors` 的最近邻结果，用相邻观测解释 KNN 插补值来源。
neighbors = all_other_points[indices.flatten()]
# 对五个最近邻的完整第二维取均值，得到缺失坐标的 KNN 插补值。
imputed_value = np.mean(neighbors[:, 1])
# 固定 `imputed_point` 的数值分量，作为“KNN填充的可视化解释”的可手算输入。
imputed_point = np.array([[5, imputed_value]])

# 创建 KNN 几何图画布，展示待填补点、邻居和局部连线。
plt.figure(figsize=(10, 7))

# 绘制所有点
plt.scatter(all_other_points[:, 0], all_other_points[:, 1], c='lightblue', label='完整数据点', s=60, alpha=0.8, edgecolors='black')
# 突出邻居
plt.scatter(neighbors[:, 0], neighbors[:, 1], c='orange', s=120, label=f'K={k}个最近邻', marker='*')
# 绘制填充后的点
plt.scatter(imputed_point[:, 0], imputed_point[:, 1], c='red', s=200, label='填充后的数据点', marker='X', edgecolors='black', linewidth=2)

# 连接到邻居
for neighbor in neighbors:
    # 从待填补点向当前最近邻画虚线，逐条展示局部参照集合。
    plt.plot([imputed_point[0,0], neighbor[0]], [imputed_point[0,1], neighbor[1]], 'k--', alpha=0.5)

# 将图题设为“K-近邻 (KNN) 填充的直观解释”，直接说明当前图形的比较目的。
plt.title('K-近邻 (KNN) 填充的直观解释', fontsize=16)
# 将横轴标为“账面市值比（已知）”，明确横向编码的变量。
plt.xlabel('账面市值比（已知）', fontsize=12)
# 将纵轴标为“资产回报率 (部分缺失)”，明确纵向编码的变量。
plt.ylabel('资产回报率 (部分缺失)', fontsize=12)
# 显示“KNN填充的可视化解释”图例，使颜色或线型与比较对象一一对应。
plt.legend(fontsize=12)
# 为“KNN填充的可视化解释”添加辅助网格，便于比较位置、斜率或组间差异。
plt.grid(True)
# 显示缺失点、候选样本与最近邻连线，检查 KNN 填充值来自局部相似观测而非总体均值。
plt.show()

# 为“步骤1：创建并观察我们的数据集”导入 `pandas` 并绑定 `pd`，用于整理当前任务的表格、字段与时间索引。
import pandas as pd
# 为“步骤1：创建并观察我们的数据集”导入 `numpy` 并绑定 `np`，用于执行当前任务的数组、数值或随机机制计算。
import numpy as np

# 固定随机种子，使 MAR 与 MCAR 机制演示的完整基准数据可复算。
np.random.seed(42)
# 为 `data` 建立包含 3 项的键值记录。
data = {
    'pps': np.random.normal(20, 5, 1000),      # 每股股价
    'book_to_market_ratio': np.random.uniform(0.1, 1.5, 1000),   # 账面市值比
    'roa': np.random.normal(0.05, 0.02, 1000)  # 资产回报率
}
# 将三项完整特征整理为插补前的真值参照表。
df_orig = pd.DataFrame(data)
# 创造变量间的关联：bm不仅是随机的，还与pps和roa相关
df_orig['book_to_market_ratio'] = df_orig['book_to_market_ratio'] + df_orig['pps'] * 0.01 + df_orig['roa'] * 2

# 让 B/M 缺失概率依赖已观测每股股价 `pps`，构造受控 MAR 近似机制。
# 股价(pps)越高的公司，其账面市值比(book_to_market_ratio)越可能缺失
# 我们使用sigmoid函数来创建一个平滑的概率
missing_prob = 1 / (1 + np.exp(-(df_orig['pps'] - 25)))
# 按给定概率抽取 `is_missing` 的 0/1 标签，固定“步骤1：创建并观察我们的数据集”的分类结果。
is_missing = np.random.binomial(1, missing_prob.values, size=1000).astype(bool)
# 深拷贝原始机制数据，隔离后续缺失注入与填补操作，保留未修改分布作为比较基准。
df = df_orig.copy()
# 按与每股股价 `pps` 相关的缺失概率遮蔽账面市值比，构造 MAR 情景。
df.loc[is_missing, 'book_to_market_ratio'] = np.nan

# 断言 MAR 驱动变量在最终观测表中始终可见，避免联合遮蔽把机制改成混合机制。
assert df['pps'].notna().all(), 'MAR 驱动变量 pps 必须保持完全观测'

# 标出注入受控 MAR 缺失后的前五行。
print('数据集的前5行 (包含缺失值NaN):')
# 输出含缺失值数据前五行，定位账面市值比的空缺位置。
print(df.head())
# 报告账面市值比的缺失占比，核对 MAR 遮蔽强度。
print(f'\n账面市值比(book_to_market_ratio)的缺失比例: {df["book_to_market_ratio"].isnull().sum() / len(df):.2%}')
# 报告驱动变量完整性，证明当前 B/M 遮蔽只依赖观测量。
print(f'MAR驱动变量每股股价(pps)的缺失比例: {df["pps"].isnull().sum() / len(df):.2%}')

# 保存原始的bm列，用于后续比较
original_book_to_market_for_plot = df_orig.loc[df['book_to_market_ratio'].isnull(), 'book_to_market_ratio']

# 为“步骤2：使用简单均值填充”，从 `sklearn.impute` 导入`SimpleImputer` 用于按指定统计量填补缺失值。
from sklearn.impute import SimpleImputer

# 建立逐列均值填充器，作为不利用变量间关系的简单基准。
mean_imputer = SimpleImputer(strategy='mean')

# 用完整机制表的已观测值逐列估计均值，并生成单次均值填补结果。
#    Transform: 学习到的均值被用来填充NaN
df_mean_imputed = pd.DataFrame(mean_imputer.fit_transform(df), columns=df.columns)

# 标出均值填充后的前五行。
print('使用均值填充后的数据集的前5行:')
# 输出均值填充后的前五行，确认缺失单元已被训练样本均值替代。
print(df_mean_imputed.head())
# 检查均值填充后账面市值比是否已无缺失值。
print(f'\n填充后账面市值比(book_to_market_ratio)是否还有缺失值? {df_mean_imputed["book_to_market_ratio"].isnull().sum() == 0}')

# 为“步骤3：单次迭代插补演示”，启用 scikit-learn 的实验性迭代插补器接口。
from sklearn.experimental import enable_iterative_imputer
# 为“步骤3：单次迭代插补演示”，从 `sklearn.impute` 导入 `IterativeImputer`，拟合缺失值填补器。
from sklearn.impute import IterativeImputer

# 建立单次确定性 IterativeImputer，明确其不生成多重插补数据集与 Rubin 合并推断。
#    max_iter=10: 迭代10轮
#    random_state=0: 保证结果可以复现
iterative_imputer = IterativeImputer(max_iter=10, random_state=0)

# 在同一机制表上执行单次确定性迭代填补，并保留字段名用于分布比较。
df_iterative_imputed = pd.DataFrame(iterative_imputer.fit_transform(df), columns=df.columns)

# 标出单次确定性迭代插补后的前五行。
print('使用单次确定性迭代插补后的数据集的前5行:')
# 输出单次迭代插补后的前五行，检查预测式填充值与其他变量共同变化。
print(df_iterative_imputed.head())
# 检查单次迭代插补后账面市值比是否已无缺失值。
print(f'\n填充后账面市值比(book_to_market_ratio)是否还有缺失值? {df_iterative_imputed["book_to_market_ratio"].isnull().sum() == 0}')

# 提取被填充的值
mean_imputed_values = df_mean_imputed.loc[is_missing, 'book_to_market_ratio']
# 提取 `iterative_imputed_values` 的插补结果，用于与原始非缺失分布和离散程度比较。
iterative_imputed_values = df_iterative_imputed.loc[is_missing, 'book_to_market_ratio']
# 从未破坏的真值表提取被遮蔽位置，作为插补恢复误差的对照。
original_values = df_orig.loc[is_missing, 'book_to_market_ratio']

# 为“条件性比较：检查填充结果”建立投影画布，给坐标、图例与标注预留空间。
plt.figure(figsize=(12, 7))

# 累计分布不会因均值填充产生大量相同值而出现不可读的密度尖峰。
sns.ecdfplot(original_values, label='原始真实值', color='#2C7BB6', linewidth=3)
sns.ecdfplot(mean_imputed_values, label='均值填充值', color='#D7191C', linewidth=3, linestyle='--')
sns.ecdfplot(iterative_imputed_values, label='单次迭代插补值', color='#1B9E77', linewidth=3, linestyle='-.')

# 少量极端插补值会把主体分布压成一条线；图中明确只放大 1%—99% 区间。
all_comparison_values = pd.concat([original_values, mean_imputed_values, iterative_imputed_values])
display_lower, display_upper = all_comparison_values.quantile([0.01, 0.99])
plt.xlim(display_lower, display_upper)
outside_count = int(((all_comparison_values < display_lower) | (all_comparison_values > display_upper)).sum())
plt.text(0.99, 0.96, f'放大显示 1%—99% 区间；区间外 {outside_count} 个值',
         transform=plt.gca().transAxes, ha='right', va='top', fontsize=11)

# 将图题设为“不同方法填充值的累计分布”，直接说明当前图形的比较目的。
plt.title('不同方法填充值的累计分布', fontsize=16)
# 将横轴标为“账面市值比 (book_to_market_ratio) 值”，明确横向编码的变量。
plt.xlabel('账面市值比 (book_to_market_ratio) 值', fontsize=12)
# 纵轴表示小于等于当前横轴取值的样本比例。
plt.ylabel('累计样本比例', fontsize=12)
# 显示“条件性比较：检查填充结果”图例，使颜色或线型与比较对象一一对应。
plt.legend(fontsize=12)
# 为“条件性比较：检查填充结果”，压缩子图留白，避免标题和坐标标签相互遮挡。
plt.tight_layout()
# 呈现多种插补值的累计分布，检查方法间的位置与离散程度差异。
plt.show()

# 提取原始完整 B/M 列，作为方差保持程度的未缺失参照。
book_to_market_original_full = df_orig['book_to_market_ratio']
# 提取 `book_to_market_mean_imputed_full` 的插补结果，用于与原始非缺失分布和离散程度比较。
book_to_market_mean_imputed_full = df_mean_imputed['book_to_market_ratio']
# 提取 `book_to_market_iterative_imputed_full` 的插补结果，用于与原始非缺失分布和离散程度比较。
book_to_market_iterative_imputed_full = df_iterative_imputed['book_to_market_ratio']

# 计算原始 B/M 标准差，作为比较各填补方法分布压缩程度的基准。
std_original = book_to_market_original_full.std()
# 计算 `std_mean`（标准差），概括“量化分析：标准差的对比”中的样本特征。
std_mean = book_to_market_mean_imputed_full.std()
# 计算 `std_iterative`（标准差），概括“量化分析：标准差的对比”中的样本特征。
std_iterative = book_to_market_iterative_imputed_full.std()

# 打印原始 B/M 标准差，随后与均值、KNN 和迭代填补结果并列核对。
print(f"原始数据的标准差: {std_original:.4f}")
# 报告均值填充后的标准差及其相对收缩幅度。
print(f"均值填充后标准差: {std_mean:.4f} (被人为降低了 { (1 - std_mean/std_original):.2%} )")
# 报告单次迭代插补后的标准差，只用于描述性比较。
print(f"单次迭代插补后标准差: {std_iterative:.4f} （仅报告差异，不作优劣宣称）")

from pathlib import Path  # 管理固定公开估值数据版本路径
from urllib.request import urlretrieve  # 复用本章隐藏设置单元安装的浏览器标识下载器
import time  # 记录三种插补方法的运行时间
import numpy as np  # 处理对数变换与确定性遮蔽
import pandas as pd  # 构造公司年度面板与结果表
# 按 Linux 共享数据、Windows 共享数据、项目缓存的顺序选择估值因子文件。
valuation_path = next((candidate_path for candidate_path in [Path('/home/ubuntu/r2_data_mount/data/stock/valuation_factors_quarterly_15_years.h5'), Path('C:/qiufei/data/stock/valuation_factors_quarterly_15_years.h5'), Path('data/course/valuation_factors_quarterly_15_years.h5')] if candidate_path.exists()), Path('data/course/valuation_factors_quarterly_15_years.h5'))
if not valuation_path.exists():
    valuation_path.parent.mkdir(parents=True, exist_ok=True)
    urlretrieve('https://assets.qiufei.site/data/stock/valuation_factors_quarterly_15_years.h5', valuation_path)
valuation_key = 'valuation_factors'  # 固定 HDF5 数据键
company_pool = ['600276.XSHG', '600009.XSHG', '600585.XSHG', '603259.XSHG', '002415.XSHE', '600104.XSHG', '600019.XSHG', '600406.XSHG']  # 确定八家公司
raw_valuations = pd.read_hdf(valuation_path, key=valuation_key)  # Fixed 键必须物理读取后再立即筛选
valuation_columns = ['book_to_market_ratio_lf', 'market_cap', 'dividend_yield_ttm']  # 确定建模字段
valuation_panel = raw_valuations.loc[raw_valuations.index.get_level_values('order_book_id').isin(company_pool), valuation_columns].reset_index()  # 立即缩减为固定公司与字段
valuation_panel['year'] = valuation_panel['date'].dt.year  # 从观测日提取时间外切分年份
valuation_panel = valuation_panel.query('2019 <= year <= 2024').sort_values(['order_book_id', 'year', 'date'])  # 限定数据版本实际覆盖期并排序
annual_valuations = valuation_panel.groupby(['order_book_id', 'year'], as_index=False).tail(1).copy()  # 每个公司年度保留最后观测日
annual_valuations['log_market_cap'] = np.log(annual_valuations['market_cap'])  # 将市值转为稳定的对数尺度
annual_valuations['next_bm'] = annual_valuations.groupby('order_book_id')['book_to_market_ratio_lf'].shift(-1)  # 对齐下一年 B/M 目标
imputation_panel = annual_valuations.dropna(subset=['next_bm']).copy()  # 删除没有下一年目标的 2024 特征行
assert set(company_pool) == set(imputation_panel['order_book_id'].unique()), '固定公司池不完整'  # 缺公司时明确停止
print({'文件': valuation_path.name, '键': valuation_key, '格式': 'Fixed', '公司数': imputation_panel['order_book_id'].nunique(), '样本期': (imputation_panel['year'].min(), imputation_panel['year'].max())})  # 输出输入说明
print(imputation_panel[valuation_columns].isna().mean().rename('缺失率'))  # 输出三字段原始缺失率
print(imputation_panel.groupby('year').size().rename('公司年度样本数'))  # 输出逐年样本数量变化

from sklearn.experimental import enable_iterative_imputer  # 启用 sklearn 的实验性迭代插补接口
from sklearn.impute import SimpleImputer, KNNImputer, IterativeImputer  # 比较中位数、邻居与单次链式插补
from sklearn.linear_model import Ridge  # 建立稳定的小样本预测模型
from sklearn.metrics import mean_absolute_error  # 导入平均绝对误差，评价连续目标的预测偏差
from sklearn.pipeline import make_pipeline  # 约束缩放与回归仅在训练折拟合
from sklearn.preprocessing import StandardScaler  # 在每个训练折逐特征估计均值与标准差
imputation_features = ['book_to_market_ratio_lf', 'log_market_cap', 'dividend_yield_ttm']  # 固定三项特征顺序

class FoldScaledKNNImputer:  # 把折内尺度拟合、邻居搜索与原量纲还原封装为单一候选
    def fit(self, feature_matrix):  # 只从当前训练折学习尺度、邻居库与指示器列
        feature_array = np.asarray(feature_matrix, dtype=float)  # 固定数值矩阵与列顺序
        self.scaler_ = StandardScaler().fit(feature_array)  # 仅用训练折逐特征估计均值和标准差
        scaled_features = self.scaler_.transform(feature_array)  # 在无量纲空间建立邻居距离
        self.imputer_ = KNNImputer(n_neighbors=3, add_indicator=False).fit(scaled_features)  # 确定三邻居规则
        self.indicator_indices_ = np.flatnonzero(np.isnan(feature_array).any(axis=0))  # 确定训练折出现缺失的指示器列
        return self  # 返回已确定候选以复用验证窗转换
    def transform(self, feature_matrix):  # 用训练折参数转换训练或验证特征
        feature_array = np.asarray(feature_matrix, dtype=float)  # 保持输入列与训练说明一致
        scaled_features = self.scaler_.transform(feature_array)  # 不从验证窗重估尺度
        scaled_imputed = self.imputer_.transform(scaled_features)  # 在确定无量纲空间搜索邻居并插补
        original_scale_imputed = self.scaler_.inverse_transform(scaled_imputed)  # 将三个业务特征还原到原量纲
        missing_indicators = np.isnan(feature_array[:, self.indicator_indices_]).astype(float)  # 按确定列追加缺失指示器
        return np.column_stack([original_scale_imputed, missing_indicators])  # 完成内容原量纲特征与稳定指示器

imputer_factories = {  # 确定三种候选插补器及其超参数
    # 注册带缺失指示器的中位数插补器工厂。
    '中位数': lambda: SimpleImputer(strategy='median', add_indicator=True),
    # 注册三邻居 KNN 插补器工厂。
    'KNN（折内Z-score距离）': lambda: FoldScaledKNNImputer(),
    # 注册固定种子的十轮单次迭代插补器工厂。
    '单次迭代': lambda: IterativeImputer(max_iter=10, random_state=42, sample_posterior=False, add_indicator=True)
}
validation_years = [2021, 2022, 2023]  # 确定三个连续时间外验证年

def evaluate_imputation_fold(validation_year, method_name, imputer_factory):  # 在一个扩展窗口内完成遮蔽、插补与预测
    training_rows = imputation_panel.query('year < @validation_year').copy()  # 仅使用验证年之前的公司年度
    validation_rows = imputation_panel.query('year == @validation_year').copy()  # 固定下一年时间外验证窗
    assert len(training_rows) > 0 and len(validation_rows) > 0, '扩展窗口为空'  # 样本数量变化异常时明确停止
    training_features = training_rows[imputation_features].copy()  # 复制训练特征以实施人工遮蔽
    observed_positions = np.flatnonzero(training_features['book_to_market_ratio_lf'].notna().to_numpy())  # 定位原本可见 B/M
    random_generator = np.random.default_rng(42 + validation_year)  # 为每折固定且可复算的遮蔽序列
    masked_positions = random_generator.choice(observed_positions, size=max(1, round(len(observed_positions) * 0.2)), replace=False)  # 遮蔽训练期 20% 可见 B/M
    masked_truth = training_features.iloc[masked_positions, 0].to_numpy()  # 保存遮蔽真值供还原评估
    training_features.iloc[masked_positions, 0] = np.nan  # 仅在训练折制造 MCAR 遮蔽
    start_time = time.perf_counter()  # 从拟合前开始记录运行时间
    fitted_imputer = imputer_factory().fit(training_features)  # 插补器只观察当前训练折
    imputed_training = fitted_imputer.transform(training_features)  # 转换训练特征并附缺失指示变量
    imputed_validation = fitted_imputer.transform(validation_rows[imputation_features])  # 用训练折参数转换验证特征
    prediction_model = make_pipeline(StandardScaler(), Ridge(alpha=1.0)).fit(imputed_training, training_rows['next_bm'])  # 只由当前训练折逐特征估计缩放统计量并拟合 Ridge
    validation_predictions = prediction_model.predict(imputed_validation)  # 生成下一年时间外 B/M 预测
    restored_bm = imputed_training[masked_positions, 0]  # 取回人工遮蔽位置的插补值
    elapsed_seconds = time.perf_counter() - start_time  # 记录完整插补与预测耗时
    return {'fold': validation_year, 'method': method_name, 'train_n': len(training_rows), 'valid_n': len(validation_rows), 'masked_mae': mean_absolute_error(masked_truth, restored_bm), 'oos_mae': mean_absolute_error(validation_rows['next_bm'], validation_predictions), 'bm_mean_shift': float(np.nanmean(restored_bm) - np.mean(masked_truth)), 'bm_std_shift': float(np.nanstd(restored_bm) - np.std(masked_truth)), 'seconds': elapsed_seconds}  # 返回完整折级指标

imputation_results = [evaluate_imputation_fold(year, name, factory) for year in validation_years for name, factory in imputer_factories.items()]  # 完整运行三年乘三种方法
imputation_results_table = pd.DataFrame(imputation_results)  # 汇总逐折运行输出
baseline_rows = []  # 收集训练中位数和完整行 Ridge 基线
for validation_year in validation_years:  # 对同一时间外窗口计算简单基线
    training_rows = imputation_panel.query('year < @validation_year').dropna(subset=imputation_features + ['next_bm']).copy()  # 完整行训练样本
    validation_rows = imputation_panel.query('year == @validation_year').dropna(subset=imputation_features + ['next_bm']).copy()  # 完整行验证样本
    median_predictions = np.repeat(training_rows['next_bm'].median(), len(validation_rows))  # 训练期目标中位数基线
    complete_model = make_pipeline(StandardScaler(), Ridge(alpha=1.0)).fit(training_rows[imputation_features], training_rows['next_bm'])  # 只由完整行训练样本逐特征缩放并拟合 Ridge 基线
    complete_predictions = complete_model.predict(validation_rows[imputation_features])  # 生成完整行基线预测
    baseline_rows.append({'fold': validation_year, 'median_oos_mae': mean_absolute_error(validation_rows['next_bm'], median_predictions), 'complete_row_oos_mae': mean_absolute_error(validation_rows['next_bm'], complete_predictions)})  # 保存同窗基线指标
baseline_results = pd.DataFrame(baseline_rows)  # 形成逐折基线表

fold_results_display = imputation_results_table.rename(columns={'fold': '验证年', 'method': '方法', 'train_n': '训练', 'valid_n': '验证', 'masked_mae': '遮蔽MAE', 'oos_mae': '窗外MAE', 'bm_mean_shift': '均值移位', 'bm_std_shift': '标准差移位', 'seconds': '秒'})  # 缩短投影表头并保留全部九列原始含义
display(fold_results_display.round(5))  # 输出逐折三方法完整结果

imputation_summary = imputation_results_table.groupby('method').agg(
    mean_masked_mae=('masked_mae', 'mean'), mean_oos_mae=('oos_mae', 'mean'),
    worst_oos_mae=('oos_mae', 'max'), total_seconds=('seconds', 'sum')
).rename_axis('方法').rename(columns={
    'mean_masked_mae': '平均遮蔽误差', 'mean_oos_mae': '平均时间外误差',
    'worst_oos_mae': '最差时间外误差', 'total_seconds': '总耗时（秒）'
})
display(imputation_summary.round(5))

display(baseline_results.round(5))  # 输出训练中位数和完整行基线

mar_training = imputation_panel.query('year < 2023').copy()  # 固定2023年前记录为MAR训练样本
mar_features = mar_training[imputation_features].copy()  # 复制原始完整特征用于遮蔽
mar_probability = (mar_features['log_market_cap'].rank(pct=True) * 0.35).to_numpy()  # 令缺失概率仅依赖已观测市值排名
mar_random = np.random.default_rng(2026)  # 固定MAR遮蔽随机序列
mar_mask = mar_random.uniform(size=len(mar_features)) < mar_probability  # 按观测协变量生成缺失指示
mar_truth = mar_features.loc[mar_mask, 'book_to_market_ratio_lf'].copy()  # 保存人工遮蔽位置真值
mar_features.loc[mar_mask, 'book_to_market_ratio_lf'] = np.nan  # 只遮蔽B/M目标特征
mar_rows = []  # 收集三种方法的还原误差
for method_name, imputer_factory in imputer_factories.items():  # 在同一MAR遮蔽上比较候选方法
    mar_fitted_imputer = imputer_factory().fit(mar_features)  # 只用当前训练窗确定尺度、权重与插补参数
    mar_imputed = mar_fitted_imputer.transform(mar_features)  # 用同一确定变换还原遮蔽值并保持原量纲
    mar_rows.append({'method': method_name, 'masked_n': int(mar_mask.sum()), 'mar_recovery_mae': mean_absolute_error(mar_truth, mar_imputed[mar_mask, 0])})  # 计算遮蔽位置MAE
display(pd.DataFrame(mar_rows).round(5))  # 报告受控MAR恢复结果

reference_imputer = IterativeImputer(max_iter=10, random_state=42, sample_posterior=False).fit(mar_features)  # 以单次迭代插补建立MAR参考
reference_values = reference_imputer.transform(mar_features)[:, 0]  # 提取参考B/M插补值
delta_rows = []  # 收集不同不可识别偏移下的误差
for delta_bm in [-0.10, 0.00, 0.10]:  # 预设悲观、中性和乐观模式偏移
    adjusted_values = reference_values[mar_mask] + delta_bm  # 只调整缺失组的潜在B/M
    delta_rows.append({'delta_bm': delta_bm, 'recovery_mae': mean_absolute_error(mar_truth, adjusted_values), 'adjusted_mean': adjusted_values.mean()})  # 量化结论对MNAR偏移的敏感度
display(pd.DataFrame(delta_rows).round(5))  # 展示不可识别假设改变时的结果范围