目录

中文EN

理论篇：因子模型与风格

1、单因子模型代码2、多因子模型理论基础3、经典风格因子4、因子版图扩展

操作篇：因子生产与管理

5、因子构建与中性化6、因子测试框架代码7、因子组合优化与归因8、因子择时与拥挤

拓展篇：宏观量化与资产配置

9、宏观数据与因子10、宏观状态识别11、大类资产配置代码12、宏观量化十年笔记代码

6、因子测试框架

因子测试是因子研究中的重要环节，其目的在于系统评估因子的预测能力、稳定性以及实际可交易性。

在因子正式用于组合构建或实盘应用之前，有必要对其进行充分检验，判断其是否具有持续有效的截面解释力，以及在交易成本约束下是否仍具备可实施性。若测试框架设计不完整，回测结果即使表面上表现良好，也可能难以在样本外或实盘环境中保持一致。

部分研究者在回归中见到t值大于2，即认为找到了有效因子。然而，回归表现良好并不等同于因子具有稳定的选股能力——曾出现某因子回归结果理想，但分组回测时第三组与第四组收益乱序的情况，说明其可能只是统计上的巧合。因此，分层回测与IC分析是检验因子的关键方法。

本章将围绕四个核心模块展开：分层回测、IC/IR分析、因子收益序列分析与换手率分析。上述四个维度相互补充，构成因子有效性评估的基本框架。推荐工作顺序为：先做IC分析快速筛除无效因子，再做分层回测与收益序列验证，最后用换手率评估可交易性。

6.1 分层回测

分层回测是因子初筛中最常用的方法之一。其基本思路是：按照因子值对股票进行排序，并将其划分为若干组，进而考察各组未来收益的差异。

若因子具有有效性，分组收益通常应呈现一定的单调性，例如高因子组收益较高、低因子组收益较低，且中间各组大致依次排列。反之，若分组收益缺乏规律性，则该因子的有效性通常较弱。

6.1.1 分组构建方法

分组逻辑如下：每月末将全市场股票按因子值排序，再等分为N组。

实践中常用10组（十分位法），既能观察单调性，又可避免每组样本过少。若研究高频因子，5组亦够用。核心原则是，分组必须严格基于因子值，不能引入未来信息。

具体步骤如下：

1. 排序：每月末，对所有股票按因子值从大到小排序
2. 分组：等分成N组，第1组是因子值最大的，第N组是最小的
3. 计算收益：每组内股票等权或市值加权，计算下个月的收益
4. 滚动：每月重复，得到每组的时间序列收益

需注意：若因子值分布不均匀（如极端值较多），等分法可能导致各组股票数量差异较大。建议先做MAD去极值，再进行分组。

分层回测结果也容易受到市值、行业等共同风险暴露的影响。若某一组样本主要集中于小市值股票，而另一组主要由大市值股票构成，则分组收益差异可能反映的是市值效应，而非因子本身的预测能力。因此，在分层回测中应尽可能控制相关风险暴露。

可进一步对市值中性化后的因子进行分层回测，即先将因子对市值进行回归，取残差作为新的因子，再进行分组测试。这样有助于剔除市值干扰，更准确地识别因子的独立解释力。

若分组后某一组样本数量过少，例如少于10只股票，则该组统计结果的可靠性将显著下降。实践中通常要求每组至少包含20只股票；若难以满足该条件，则可考虑合并相邻组，以提高分组检验的稳定性。

分组标签可用如下方式生成：

import pandas as pd
import numpy as np

# Python伪代码：十分位分组
def factor_quantile(factor_df, n_groups=10):
 """
 factor_df: 日期为index, 股票代码为columns的因子值矩阵
 返回分组标签矩阵
 """
 # 每月排序分组
 groups = factor_df.apply(
 lambda x: pd.qcut(x.rank, q=n_groups, labels=False) + 1,
 axis=1
 )
 return groups

完整的分层回测框架示例如下：

def分层回测(factor_df, return_df, n_groups=10):
 """
 factor_df: 因子值，index为日期，columns为股票代码
 return_df: 下期收益，index为日期，columns为股票代码
 """
 results = {}
 for date in factor_df.index:
 # 获取当期因子值
 factors = factor_df.loc[date].dropna
 # 按因子值排序并分组
 ranked = factors.rank
 group_labels = pd.qcut(ranked, n_groups, labels=False)
 # 计算每组下期收益
 returns = return_df.loc[date]
 for g in range(n_groups):
 mask = group_labels == g
 group_returns = returns[mask.index[mask]].mean
 results.setdefault(g, []).append(group_returns)
 return pd.DataFrame(results)

  File "<ipython-input-2-ff26dcdeb80a>", line 1
    def分层回测(factor_df, return_df, n_groups=10):
                                              ^
SyntaxError: invalid syntax

6.1.2 多空组合收益

多空组合，就是做多第1组（因子值最大），做空第N组（因子值最小）。

该收益反映：若纯粹依据该因子选股，可获得多少超额收益。

考察多空组合的原因在于剥离市场整体走势的影响。例如，若市场上涨20%、第1组上涨25%、第5组上涨23%，因子创造的增量有限；而多空收益为2%，才更能反映因子的真实贡献。

从经验来看，多空收益年化低于5%的因子，通常不宜纳入多因子模型，扣除交易成本后剩余空间有限。

计算多空收益时，要注意两个细节：

• 换仓频率：月度换仓是最常见的，但高频因子可能需要周度
• 交易成本：多空组合换手率通常很高，一定要扣除千分之一到千分之二的单边成本

# 计算多空组合收益
long_short_ret = group_returns.iloc[:, 0] - group_returns.iloc[:, -1]
# 扣除交易成本（假设单边0.15%）
turnover = estimate_turnover(group_assignments)
long_short_ret_net = long_short_ret - turnover * 0.0015 * 2

6.1.3 单调性检验

一个有效的因子，分组收益应该是单调的。也就是说，从第1组到第N组，收益应该逐渐下降（或上升）。

如果第3组收益比第2组还高，那就说明因子在某些区间失效了。

常用的检验方法有三种：

方法	说明	适用场景
收益单调性	每组平均收益是否严格单调	初步筛选，快速排除明显无效因子
单调性统计量	计算组间收益差的正负比例	正式报告，量化单调程度
Patton-Timmermann检验	统计检验分组收益是否单调	学术论文，严谨性要求较高时使用

单调性不等于线性。因子收益可以是凸的（两端收益高，中间低），也可以是凹的。只要方向一致，就算单调。

曾出现某因子前5组收益单调、第6组突然跳升的情况，后查明系第6组中含个别异常个股所致。因此，单调性检验需结合异常值分析。

6.1.4 分组净值曲线

数值再精确，也不及图表直观。分组净值曲线将各组的累计收益可视化呈现。

好的因子，其净值曲线通常具备以下特征：

• 第1组在最上面，第N组在最下面
• 各组之间没有交叉
• 多空组合的净值曲线稳步向上

如果曲线交叉了，说明因子在某些时间段失效了。如果多空净值曲线大起大落，说明因子风险很高。

绘制净值曲线时，可采用对数坐标，以更清晰地观察不同时间段的收益差异，尤其是早期与后期。

以下总结实战中需特别注意的问题：

• 幸存者偏差：应使用全样本，包括退市股票。若仅包含存续股票，因子表现可能被高估——差质股票已退市。
• 前视偏差：分组所用因子值必须是调仓日收盘后可获得的。若使用未来数据，则结果无效。
• 行业集中度：若第1组以某一行业为主，因子反映的可能是行业效应而非选股能力。建议做行业中性化处理。
• 小市值效应：许多因子在小市值股票上表现更好。可按市值分层检验，确认因子在各市值区间均有效。
• 市值与行业控制：可对市值中性化后的因子再做分层回测，或在分组时控制行业暴露，以识别因子的独立解释力。

分层回测是因子投资的基础环节。分组构建需严谨，多空收益需扣除成本，单调性需统计检验，净值曲线需直观呈现。完成上述步骤后，因子价值基本可判定。

6.2 IC/IR分析

分层回测关注的是分组收益表现，而IC/IR分析则从相关系数的角度评估因子的预测能力。IC分析可快速筛除大量无效因子，通常作为因子投资的第一道防线。

6.2.1 IC的定义与计算

IC（Information Coefficient，信息系数），衡量因子值对下期收益的预测能力。换言之：因子值较高的股票，下期是否表现更优？

IC有两种主流计算方式：

• Spearman秩相关系数：实践中通常优先采用。该方法仅关注排序、不依赖具体数值，对极端值不敏感。在A股等易出现极端行情的市场中，秩相关较普通Pearson相关更为稳健。
• Pearson相关系数：要求线性关系，对异常值敏感。除非已做严格正态化处理，否则不建议直接使用。

计算IC的步骤如下：

1. 在每个调仓日，计算所有股票的因子值
2. 计算这些股票在下一期的收益（比如下一天、下一周）
3. 计算因子值和下期收益的相关系数

代码实现如下：

from scipy.stats import spearmanr

def calc_ic(factor_series, forward_return_series):
 """
 计算截面IC
 factor_series: 某一天的因子值，index为股票代码
 forward_return_series: 对应的下期收益
 """
 # 去掉缺失值
 valid = factor_series.notna & forward_return_series.notna
 f = factor_series[valid]
 r = forward_return_series[valid]

 # Spearman秩相关
 ic, p_value = spearmanr(f, r)
 return ic, p_value

# 示例：某一交易日的因子值与下期收益
# ic_value, p_val = calc_ic(factor_data['2024-01-05'],
# forward_ret['2024-01-05'])
# print(f"IC = {ic_value:.4f}, p-value = {p_val:.4f}")

亦可写成更简洁的单期接口：

def calculate_ic(factor_series, return_series):
 """计算单期IC（Spearman秩相关）"""
 combined = pd.concat([factor_series, return_series], axis=1).dropna
 if len(combined) < 30:
 return np.nan
 ic, p_value = spearmanr(combined.iloc[:, 0], combined.iloc[:, 1])
 return ic

  File "<ipython-input-4-dd7623c68873>", line 5
    return np.nan
    ^
IndentationError: expected an indented block after 'if' statement on line 4

关键点：IC值范围在[-1, 1]之间。正IC表示因子与未来收益正相关，负IC表示负相关。绝对值越大，预测能力越强。

6.2.2 IC序列的统计检验

单日IC缺乏说服力：若今日IC为0.05、次日变为-0.03，并不能据此判定因子有效。需要考察IC序列的统计显著性。

通常可进行以下分析：

• 计算IC的均值：看长期平均IC是否显著不为零
• 计算IC的t统计量：检验均值是否统计显著
• 看IC的正向比例：有多少比例的IC是正的？如果只有50%，那基本等于抛硬币

需注意：IC序列存在自相关。若直接使用普通t检验，会高估显著性。正确做法是采用Newey-West调整的标准误。

from statsmodels.stats.stattools import durbin_watson
from statsmodels.regression.linear_model import OLS

def ic_statistical_test(ic_series):
 """
 对IC序列做统计检验
 """
 # 基本统计量
 mean_ic = ic_series.mean
 std_ic = ic_series.std
 t_stat = mean_ic / (std_ic / np.sqrt(len(ic_series)))

 # Newey-West调整（处理自相关）
 X = np.ones((len(ic_series), 1))
 model = OLS(ic_series.values, X).fit(cov_type='HAC', cov_kwds={'maxlags': 5})
 nw_t_stat = model.tvalues[0]

 # 正向比例
 positive_ratio = (ic_series > 0).mean

 return {
 'mean_ic': mean_ic,
 'std_ic': std_ic,
 't_stat': t_stat,
 'nw_t_stat': nw_t_stat,
 'positive_ratio': positive_ratio
 }

通常要求IC均值至少大于0.03，且t统计量大于2（对应95%置信水平）。若IC均值仅为0.01，即使统计显著，实际交易中也难以覆盖交易成本。

关键指标参考：

指标	优秀	良好	一般	较差
IC均值	>0.05	0.03~0.05	0.01~0.03	<0.01
IR	>0.5	0.3~0.5	0.1~0.3	<0.1
IC正比例	>60%	55%~60%	50%~55%	<50%

一种误区是，仅依据IC均值判断因子的优劣。实际上，若IC均值为0.04，但标准差高达0.15，则对应的IR仅为0.27，说明该因子的预测能力并不稳定。因此，在评估因子有效性时，应更加重视IR（或ICIR）而非单纯依赖IC均值。

6.2.3 IC衰减分析

IC衰减指因子预测能力随时间推移而下降的速度。产生衰减的原因包括：市场结构变化、套利行为增加、因子被广泛使用导致超额收益摊薄等。

通常可考察不同持有期的IC：

• 持有1天的IC（日频IC）
• 持有5天的IC（周频IC）
• 持有20天的IC（月频IC）
• 持有60天的IC（季频IC）

如果因子只在第1天有效，第5天就衰减到零，说明这是个短期反转因子。如果因子在20天内都保持稳定IC，那就是趋势因子。

例如，某动量因子的月度IC均值为0.05，t值为3.2，表面上看具有一定有效性；但进一步的IC衰减分析表明，该因子仅对未来1个月收益具有解释力，至第2个月时IC已转为负值。这说明该因子更适用于短期持有，而不适合长期配置。

def decay_analysis(factor_df, forward_returns_dict):
 """
 分析IC衰减
 factor_df: 每日因子值，index为日期，columns为股票
 forward_returns_dict: 不同持有期的收益，如{'1d': ..., '5d': ...}
 """
 decay_results = {}
 for horizon, ret_df in forward_returns_dict.items:
 ics = []
 for date in factor_df.index:
 if date in ret_df.index:
 ic, _ = calc_ic(factor_df.loc[date], ret_df.loc[date])
 ics.append(ic)
 decay_results[horizon] = np.mean(ics)

 return decay_results

# 输出示例
# {'1d': 0.042, '5d': 0.038, '20d': 0.025, '60d': 0.008}

  File "<ipython-input-6-08a0b285071e>", line 9
    ics = []
    ^
IndentationError: expected an indented block after 'for' statement on line 8

曾出现某因子日频IC高达0.08、5日IC降至0.01的情况。此类快速衰减的因子往往仅适用于高频策略，难以在日频及以上频率覆盖交易成本。

6.2.4 ICIR：信息系数与信息比率

ICIR（Information Coefficient Information Ratio），是IC的均值除以IC的标准差。它衡量的是因子预测能力的稳定性。文中亦常简称IR。

公式如下：

$$ ICIR = mean(IC) / std(IC) $$

ICIR的重要性可通过下例说明：

• 因子A：IC均值0.05，标准差0.02，ICIR = 2.5
• 因子B：IC均值0.08，标准差0.10，ICIR = 0.8

因子A的IC虽较低，但稳定性高；因子B的IC较高，但波动较大。实际交易中，因子A往往优于因子B——因其可分配更高权重，且不易出现某月突然失效的情况。

ICIR筛选的常见标准如下：

ICIR范围	评价	建议
> 2.0	优秀	可直接使用
1.0 - 2.0	良好	需结合其他因子
0.5 - 1.0	一般	谨慎使用，注意风险
< 0.5	差	建议放弃

ICIR比IC本身更为重要。ICIR为2.0、IC均值为0.03的因子，通常优于ICIR为0.5、IC均值为0.08的因子。稳定性是关键指标。

说明：上表为基于IC序列稳定性的ICIR标准；前文关键指标参考中的IR阈值（如优秀>0.5）多对应月度频率下的经验区间，二者量纲与样本频率不同，实践中宜按同一频率统一比较。

以下列举实战中常见的问题：

1. IC为正但策略亏损：可能是因子在极端行情下失效。需检查IC时间序列，观察牛熊转换期IC是否转负。
2. ICIR很高但实盘表现差：可能存在过拟合。若样本内ICIR高达3.0、样本外仅0.5，则可能是数据挖掘的结果。
3. IC衰减过快：可能是因子已被市场广泛认知，或交易成本过高。建议检查换手率；若换手率超过500%，衰减较快属正常现象。

开发新因子时，通常先做IC分析，再做分层回测。IC分析可快速筛除大量无效因子。IC分析是因子投资的第一道防线。

6.3 因子收益序列分析

分层回测和IC分析主要从截面维度评估因子表现，而因子收益序列分析则进一步从时间序列维度考察因子的有效性。

具体而言，可以构建因子模拟组合（Factor Mimicking Portfolio），每月做多因子值最高的组合，做空因子值最低的组合，并据此观察多空组合的收益表现（其构建逻辑与6.1.2节多空组合一致，此处侧重时间序列统计与风险特征）。

该收益序列可以用于判断以下几个方面：

• 因子的累计收益曲线是否具有持续上升特征；
• 最大回撤的幅度及其发生时点；
• 收益波动是否处于可接受范围；
• 是否存在明显的风格暴露，尤其是在特定市场环境下的失效风险。

例如，曾有一个估值因子的多空组合年化收益率达到8%，但最大回撤亦达到25%。进一步分析发现，回撤主要集中于2015年牛市后期和2017年白马股行情阶段。这表明，该因子在特定风格环境下可能失效。若用于全市场选股，其适用性尚可；但若策略仅覆盖某一类股票，则需谨慎评估其稳定性。

因子收益序列的关键统计指标包括：年化收益率、年化波动率、夏普比率、最大回撤、Calmar比率（年化收益/最大回撤）以及月度胜率。

6.4 换手率分析

换手率分析是因子测试中不可忽视的重要环节。即使因子在收益表现上具有优势，若其换手率过高，交易成本仍可能显著侵蚀收益，甚至使策略的实际可行性大幅下降。

换手率的计算方式如下：

$$ 换手率 = （买入金额 + 卖出金额）/ 组合总市值 $$

对于多空组合，需分别统计多头端和空头端的换手率。通常可重点关注以下指标：

• 单边换手率：多头或空头每期的换手率；
• 双边换手率：多头与空头合计的总换手率；
• 年化换手率：单边换手率×调仓频率；
• 收益换手比：年化收益/年化换手率，用于衡量单位换手所对应的收益水平。

换手率分析应结合实际交易成本进行评估。A股市场的双边交易成本（包括佣金、印花税及冲击成本）通常约为0.1%~0.3%。若因子换手率过高，即便IC较高，实际收益也可能被交易成本显著削弱。

在实践中，某高频反转因子的原始换手率曾高达800%，年化收益率为12%。扣除交易成本后，净收益降至3%。随后通过加入筛选条件将换手率降至200%，虽然年化收益率下降至9%，但扣除成本后净收益提升至6%。这一结果表明，在部分情形下，降低换手率反而有助于提升策略的实际收益。

实操中，一般会设定一个换手率上限。对于月度调仓的因子，年化换手率超过300%就要警惕；对于周度调仓的因子，年化换手率超过1000%基本不可行。具体阈值取决于交易成本和策略容量。

换手率分析的代码示例：

def calculate_turnover(weights_prev, weights_curr):
 """
 计算单期换手率
 weights_prev: 上期持仓权重
 weights_curr: 本期持仓权重
 """
 # 买入量：本期权重增加的部分
 buy = (weights_curr - weights_prev).clip(lower=0).sum
 # 卖出量：本期权重减少的部分
 sell = (weights_prev - weights_curr).clip(lower=0).sum
 # 换手率取买入和卖出的平均值
 turnover = (buy + sell) / 2
 return turnover

最后，对因子测试框架作一简要总结。该框架主要包括以下四个模块，各自侧重点不同：

• 分层回测：用于检验因子分组收益的单调性，并识别其是否具有稳定的截面解释力（含多空收益、单调性检验、净值曲线与偏差控制）；
• IC/IR分析：用于评估因子的预测能力及其稳定性（含统计检验、衰减分析与ICIR）；
• 因子收益序列：用于考察因子的实际收益表现及风险特征；
• 换手率分析：用于衡量因子的交易成本与实际可行性。

推荐工作顺序：先做IC分析快速筛除无效因子，再做分层回测与收益序列验证，最后用换手率评估可交易性。上述四个维度均通过检验后，因子方可视为初步具备有效性，并进一步进入多因子组合构建阶段。

需要说明的是，因子测试并非一次性完成的工作。随着市场环境持续变化，因子的有效性也可能发生调整。因此，宜定期对测试框架进行复核，以判断因子是否仍具备稳定表现。大量因子在样本外阶段失效属于常见现象，应将其视为因子研究中的正常风险。

6.5 规范组合排序：单维与双维

分层回测与 IC 分析直观，但在学术复现与严格检验中，通常还需要明确：断点如何设定、是否控制另一维度（如市值）、多空收益如何定义。规范的组合排序（portfolio sorts）是实证资产定价中的标准流程。

6.5.1 单维组合排序（Univariate Portfolio Sorts）

基本步骤

1. 在每个再平衡时点（如每月末），按因子值对股票排序；
2. 依据分位数断点划分为 Q 组（常见为5组或10组）；
3. 组内按等权或市值加权计算持有期收益；
4. 构造高减低（或低减高）多空组合，形成因子收益时间序列；
5. 对多空收益做均值、t 检验（常用 Newey-West 调整）。

与前文分层回测一致的是“分组—持有—比较”；学术写法更强调：断点样本（是否仅用特定交易所样本计算断点）、滞后对齐（避免前视偏差）以及多空定义的符号方向。

实现要点

一套完整的排序流程通常包括：断点计算、组合标签分配、组内收益聚合，以及高减低多空收益构造。可在自有代码中模块化实现上述步骤。

在 A 股应用时，面板字段一般包括股票代码、交易日、收益、市值与因子值，并在样本过滤中加入 ST、停牌、涨跌停与上市天数等规则。

# 示意：单维排序的最小可运行骨架（演示数据）
import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
dates = pd.date_range("2020-01-31", periods=24, freq="ME")
rows = []
for dt in dates:
 n = 300
 factor = rng.normal(size=n)
 ret = 0.02 * factor + rng.normal(0, 0.05, size=n)
 mktcap = np.exp(rng.normal(22, 1, size=n))
 rows.append(pd.DataFrame({
 "date": dt, "id": np.arange(n), "factor": factor, "ret": ret, "mktcap": mktcap,
 }))
panel = pd.concat(rows, ignore_index=True)

def univariate_sort(df, q=5, weight="equal"):
 out = []
 for dt, g in df.groupby("date"):
 g = g.copy
 g["port"] = pd.qcut(g["factor"], q=q, labels=False, duplicates="drop") + 1
 if weight == "equal":
 port_ret = g.groupby("port")["ret"].mean
 else:
 def wavg(x):
 w = x["mktcap"] / x["mktcap"].sum
 return np.sum(w * x["ret"])
 port_ret = g.groupby("port").apply(wavg, include_groups=False)
 long_short = port_ret.loc[q] - port_ret.loc[1]
 out.append({"date": dt, "long_short": long_short})
 return pd.DataFrame(out).set_index("date")

res = univariate_sort(panel, q=5)
print(res[["long_short"]].mean)
print(
 "t-stat (iid approx):",
 res["long_short"].mean / (res["long_short"].std(ddof=1) / np.sqrt(len(res))),
)

  File "<ipython-input-8-a86b468f3358>", line 21
    g = g.copy()
    ^
IndentationError: expected an indented block after 'for' statement on line 20

6.5.2 双维组合排序（Bivariate Sorts）

单维排序无法回答：因子收益是否只是市值、行业或其他风格的代理。双维排序的常见做法是：

1. 先按控制变量（如市值）分组；
2. 再在每组内按目标因子分组；或采用独立排序形成二维网格；
3. 在控制维度上平均，得到“控制后”的因子组合收益。

Fama-French 的 Size-BM 2×3 或 5×5 组合即属此类。

实务含义：若目标因子在市值中性双维排序后多空收益仍显著，其增量信息才更可能独立于规模效应。

6.6 Fama-MacBeth 回归

Fama-MacBeth（1973）回归用于估计特征/暴露的风险溢价，与组合排序互补：排序看经济意义与单调性，FM 看统计定价与多变量稳健性。

第一步（截面）：对每个时期 t，估计

$$ R_{i,t}= \gamma_{0,t}+\gamma_{1,t}X_{i,t-1}+\varepsilon_{i,t} $$

其中 (X_{i,t-1}) 为滞后一期的因子暴露或公司特征。

第二步（时间序列）：对 ({\gamma_{1,t}}) 取均值得到风险溢价估计，并用时间序列标准误（常用 Newey-West）判断是否显著异于 0。

方法	主要问题	优点	局限
组合排序	高暴露组合是否赚得更多	直观、可呈现非线性	难同时控制多个变量
Fama-MacBeth	控制其他变量后暴露是否被定价	可多变量控制	依赖线性设定，对异常值敏感

# 示意：手写两步 Fama-MacBeth（沿用上一单元格的 panel）
import statsmodels.api as sm

def fama_macbeth(panel, y="ret", x_cols=("factor",), date_col="date"):
 gammas = []
 for dt, g in panel.groupby(date_col):
 if len(g) < 30:
 continue
 X = sm.add_constant(g.loc[:, list(x_cols)])
 model = sm.OLS(g[y], X).fit
 gammas.append({"date": dt, **model.params.to_dict})
 gdf = pd.DataFrame(gammas).set_index("date").sort_index
 summary = pd.DataFrame({
 "risk_premium": gdf.mean,
 "std_error_iid": gdf.std(ddof=1) / np.sqrt(len(gdf)),
 })
 summary["t_iid"] = summary["risk_premium"] / summary["std_error_iid"]
 return gdf, summary

gdf, summary = fama_macbeth(panel, y="ret", x_cols=("factor",))
print(summary)

  File "<ipython-input-9-9068e354d5e3>", line 7
    if len(g) < 30:
    ^
IndentationError: expected an indented block after 'for' statement on line 6

多变量控制与标准误

实务中常在同一截面回归中同时放入多个特征，例如市场 beta、账面市值比与对数市值，以识别目标因子的增量定价信息。对第二步得到的系数序列，优先使用 Newey-West 等 HAC 标准误。

多变量示例（示意）：

# 截面期 t：ret_excess ~ factor + beta + bm + log_mktcap
# 再对各期系数做时间序列均值与 Newey-West t 检验

实施要点

1. 滞后：解释变量必须对收益严格先验可知，避免前视偏差；
2. 标准化/缩尾：截面回归对极端值敏感，宜先做去极值与标准化（见第5章）；
3. 行业或市值控制：可在回归中加入控制变量，或先中性化再回归（见第5章）；
4. 标准误：系数序列常存在自相关，优先报告 Newey-West t 值；
5. 与 IC 的关系：Rank IC 关注次序相关性；FM 关注条件风险溢价，二者不可互相替代。

6.7 面板回归中的聚类标准误（要点）

在将多期截面合并为面板，或估计个股时序回归时，残差常在个股维度或时间维度上相关。若仍使用普通 OLS 标准误，显著性容易被高估。

面板或混合回归中，残差相关结构会影响标准误。常见处理包括：

• 个体聚类：允许同一股票残差跨期相关；
• 时间聚类：允许同一时期残差在股票间相关；
• 双向聚类：同时考虑两种相关（实证资产定价中较常见）。

与 FM 的关系：经典 FM 通过对各期截面系数再做时间序列推断，本身已部分处理截面相关；若改用面板一次性回归，则更需要显式聚类。Python 中可使用 linearmodels 等库实现固定效应与聚类标准误。

课程建议：因子定价检验优先 FM + Newey-West；面板特征回归至少报告一维聚类标准误，重要结论尽量做双向聚类稳健性检查。多重检验与排序设定中的 p-hacking 见第4章（因子版图扩展）。

6.8 A股复现：单维组合排序（Univariate Portfolio Sorts）

本节完整复现 tidyfinance 章节 Univariate Portfolio Sorts 的核心流程，并把数据替换为 A 股月频样本。组合排序是实证资产定价中最常用的方法之一：按某个可观测特征把股票分成若干组合，比较各组未来收益，从而检验该特征是否具有横截面预测能力。

单维排序只使用一个排序变量 $x_{i,t-1}$（下标 $t-1$ 表示投资者在 $t$ 期初已可观测）。本节以第 2.6 节估计的市场 beta 作为排序变量，评估其与下期超额收益 $r_{i,t}$ 的关系。

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

# 中文字体（与第2章绘图设置一致）
_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("DATA_DIR =", DATA_DIR)
print("beta exists =", (DATA_DIR / "beta_ashare.parquet").exists)
print("monthly cache exists =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)

DATA_DIR = D:\A_Topics\202607_03_TidyFinanceAShare\202608_02_传统量化\Data
beta exists = True
monthly cache exists = True

6.8.1 数据准备

对应原文 Data Preparation。资产宇宙使用 A 股月频个股超额收益与上期流通市值；市场超额收益用于计算 CAPM alpha。beta 取自第 2.6 节月频滚动估计结果。

crsp_monthly = (
    pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet")
    [["permno", "date", "ret_excess", "mkt_excess", "mktcap_lag"]]
    .copy
)
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)

# 市场因子：每月一条 mkt_excess（与原文 factors_ff3_monthly 角色相同）
factors_mkt_monthly = (
    crsp_monthly[["date", "mkt_excess"]]
    .drop_duplicates("date")
    .sort_values("date")
    .reset_index(drop=True)
)

beta = (
    pd.read_parquet(DATA_DIR / "beta_ashare.parquet")
    .query("return_type == 'monthly'")
    [["permno", "date", "beta"]]
    .copy
)
beta["date"] = pd.to_datetime(beta["date"])
beta["permno"] = beta["permno"].astype(str).str.zfill(6)

print("crsp_monthly:", crsp_monthly.shape)
print("beta:", beta.shape)
print("months:", factors_mkt_monthly["date"].nunique)
display(crsp_monthly.head(3))
display(beta.head(3))

crsp_monthly: (890537, 5)
beta: (623256, 3)
months: 319

	permno	date	ret_excess	mkt_excess	mktcap_lag
0	000001	2000-01-01	0.060035	0.158982	NaN
1	000001	2000-02-01	-0.013189	0.120168	19843822.88
2	000001	2000-03-01	0.000873	0.054070	19618933.36

	permno	date	beta
0	000001	2003-12-01	0.945177
1	000001	2004-01-01	0.956954
2	000001	2004-02-01	0.978056

6.8.2 按市场 Beta 排序

对应原文 Sorting by Market Beta。用滞后一期的 beta 作为排序变量，保证建仓时信息已可知。做法是把 beta 的日期整体平移一个月再与收益表按 (permno, date) 合并，而不是简单 groupby.shift(1)——后者在存在隐式缺失月份时容易错位。

beta_lag = beta.copy
beta_lag["date"] = beta_lag["date"] + pd.DateOffset(months=1)
beta_lag = beta_lag.rename(columns={"beta": "beta_lag"}).dropna(subset=["beta_lag"])

data_for_sorts = crsp_monthly.merge(
    beta_lag, on=["permno", "date"], how="inner"
).dropna(subset=["ret_excess", "mktcap_lag", "beta_lag"])

print("data_for_sorts:", data_for_sorts.shape)
print(
    "date range:",
    data_for_sorts["date"].min.date,
    "->",
    data_for_sorts["date"].max.date,
)
display(data_for_sorts.head(3))

data_for_sorts: (614523, 6)
date range: 2004-01-01 -> 2026-07-01

	permno	date	ret_excess	mkt_excess	mktcap_lag	beta_lag
0	000001	2004-01-01	0.088847	0.073093	11993670.32	0.945177
1	000001	2004-02-01	0.114744	0.069838	13078879.04	0.956954
2	000001	2004-03-01	0.027323	0.030699	14600989.96	0.978056

组合排序的第一步是计算断点。先用滞后 beta 的中位数把股票分成 low / high 两组，再计算各组市值加权超额收益（权重为 mktcap_lag）。

def value_weighted_ret(g: pd.DataFrame) -> float:
    w = g["mktcap_lag"]
    return float((g["ret_excess"] * w).sum / w.sum)

rows = []
for dt, g in data_for_sorts.groupby("date", sort=True):
    g = g.copy
    # 中位数断点：两组
    g["portfolio"] = pd.qcut(
        g["beta_lag"].rank(method="first"),
        q=2,
        labels=["low", "high"],
    )
    for port, pg in g.groupby("portfolio", observed=True):
        if pg["mktcap_lag"].sum <= 0:
            continue
        rows.append(
            {
                "date": dt,
                "portfolio": str(port),
                "ret": value_weighted_ret(pg),
            }
        )

beta_portfolios_2 = pd.DataFrame(rows).sort_values(["date", "portfolio"])
print(beta_portfolios_2.head)
print("n months:", beta_portfolios_2["date"].nunique)

        date portfolio       ret
1 2004-01-01      high  0.099191
0 2004-01-01       low  0.056309
3 2004-02-01      high  0.087797
2 2004-02-01       low  0.068211
5 2004-03-01      high  0.035140
n months: 271

6.8.3 绩效评价（两组多空）

对应原文 Performance Evaluation。构造多空组合：做多高 beta、做空低 beta。忽略摩擦时，多空头寸对市场净暴露可接近于零（此处先看原始多空收益）。

beta_longshort_2 = (
    beta_portfolios_2.pivot(index="date", columns="portfolio", values="ret")
    .sort_index
    .dropna(subset=["low", "high"])
)
beta_longshort_2["long_short"] = beta_longshort_2["high"] - beta_longshort_2["low"]
beta_longshort_2 = beta_longshort_2.reset_index

display(beta_longshort_2.head)
print(
    "mean long_short =",
    round(beta_longshort_2["long_short"].mean, 6),
    "std =",
    round(beta_longshort_2["long_short"].std, 6),
)

portfolio	date	high	low	long_short
0	2004-01-01	0.099191	0.056309	0.042882
1	2004-02-01	0.087797	0.068211	0.019586
2	2004-03-01	0.035140	0.025616	0.009524
3	2004-04-01	-0.103090	-0.114164	0.011074
4	2004-05-01	-0.014376	-0.028354	0.013978

mean long_short = -0.001431 std = 0.041392

检验多空组合平均超额收益是否显著异于零。资产定价文献常用 Newey–West（HAC）$t$ 统计量处理自相关；原文默认滞后 6 个月。这里用 statsmodels 的 cov_type="HAC" 实现。

y = beta_longshort_2["long_short"]
X = np.ones((len(y), 1))
model_fit = sm.OLS(y, X, missing="drop").fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("Long-short ~ 1  (Newey-West, lags=6)")
print(model_fit.summary)

Long-short ~ 1  (Newey-West, lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                       nan
Date:                Sat, 08 Aug 2026   Prob (F-statistic):                nan
Time:                        21:15:11   Log-Likelihood:                 479.01
No. Observations:                 271   AIC:                            -956.0
Df Residuals:                     270   BIC:                            -952.4
Df Model:                           0                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0014      0.002     -0.722      0.470      -0.005       0.002
==============================================================================
Omnibus:                       37.393   Durbin-Watson:                   1.942
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              187.790
Skew:                          -0.363   Prob(JB):                     1.67e-41
Kurtosis:                       7.013   Cond. No.                         1.00
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

若 CAPM 成立，高 beta 组合应有更高期望收益，因而“多高空低”的平均超额收益应显著为正。若结果不显著甚至为负，则与 CAPM 的直觉相冲突——这正是后文要展开的 low-beta / betting-against-beta 线索。

6.8.4 函数化：任意分组数的排序

对应原文 Functional Programming for Portfolio Sorts。把排序逻辑封装成函数，便于把股票分到任意 $N$ 组。对聚类严重、断点重复或早期样本成分股过少的情形，需要额外谨慎（空组合、极端权重等）。

def assign_portfolio(series: pd.Series, n_portfolios: int) -> pd.Series:
    # 按分位数把排序变量映射到 1..n_portfolios
    return pd.qcut(
        series.rank(method="first"),
        q=n_portfolios,
        labels=[str(i) for i in range(1, n_portfolios + 1)],
        duplicates="drop",
    )

def form_value_weighted_portfolios(
    data: pd.DataFrame,
    sorting_variable: str = "beta_lag",
    n_portfolios: int = 10,
) -> pd.DataFrame:
    # 每月按 sorting_variable 分成 n 组，计算市值加权超额收益
    rows = []
    for dt, g in data.groupby("date", sort=True):
        g = g.copy
        g["portfolio"] = assign_portfolio(g[sorting_variable], n_portfolios)
        for port, pg in g.groupby("portfolio", observed=True):
            if pg["mktcap_lag"].sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio": str(port),
                    "ret": value_weighted_ret(pg),
                }
            )
    out = pd.DataFrame(rows)
    # 有序类别，方便后续按 1..10 排序作图
    cats = [str(i) for i in range(1, n_portfolios + 1)]
    out["portfolio"] = pd.Categorical(out["portfolio"], categories=cats, ordered=True)
    out = out.merge(factors_mkt_monthly, on="date", how="left")
    return out.sort_values(["date", "portfolio"]).reset_index(drop=True)

beta_portfolios = form_value_weighted_portfolios(
    data_for_sorts, sorting_variable="beta_lag", n_portfolios=10
)
print(beta_portfolios.head)
print("portfolios:", sorted(beta_portfolios["portfolio"].dropna.unique.tolist))

        date portfolio       ret  mkt_excess
0 2004-01-01         1  0.024797    0.073093
1 2004-01-01         2  0.066696    0.073093
2 2004-01-01         3  0.054822    0.073093
3 2004-01-01         4  0.070948    0.073093
4 2004-01-01         5  0.082065    0.073093
portfolios: ['1', '10', '2', '3', '4', '5', '6', '7', '8', '9']

6.8.5 十分组绩效：Alpha、组合 Beta 与平均收益

对应原文 More Performance Evaluation。对每个 beta 组合估计 CAPM：

$$ r_{p,t}=\alpha_p+\beta_p\,r_{m,t}+\varepsilon_{p,t}, $$

并汇总 $\hat\alpha_p$、$\hat\beta_p$ 与平均超额收益。

summary_rows = []
for port, g in beta_portfolios.groupby("portfolio", observed=True):
    g = g.dropna(subset=["ret", "mkt_excess"])
    if len(g) < 24:
        continue
    X = sm.add_constant(g["mkt_excess"])
    m = sm.OLS(g["ret"], X).fit
    summary_rows.append(
        {
            "portfolio": str(port),
            "alpha": float(m.params["const"]),
            "beta": float(m.params["mkt_excess"]),
            "ret": float(g["ret"].mean),
            "n": int(len(g)),
        }
    )

beta_portfolios_summary = pd.DataFrame(summary_rows)
cats = [str(i) for i in range(1, 11)]
beta_portfolios_summary["portfolio"] = pd.Categorical(
    beta_portfolios_summary["portfolio"], categories=cats, ordered=True
)
beta_portfolios_summary = beta_portfolios_summary.sort_values("portfolio").reset_index(drop=True)
display(beta_portfolios_summary)

	portfolio	alpha	beta	ret	n
0	1	0.003172	0.715856	0.009123	271
1	2	0.001743	0.922710	0.009414	271
2	3	0.000964	1.012170	0.009378	271
3	4	0.000680	1.041719	0.009340	271
4	5	0.000177	1.070959	0.009080	271
5	6	0.000239	1.095577	0.009346	271
6	7	-0.000868	1.146563	0.008663	271
7	8	-0.000416	1.182178	0.009411	271
8	9	-0.001407	1.215726	0.008699	271
9	10	-0.006544	1.278707	0.004086	271

下图展示十分组的 CAPM alpha。若出现“低 beta 组 alpha 偏高、高 beta 组 alpha 偏低”，则与 CAPM（风险调整后 alpha 应接近 0，且收益随 beta 上升）相矛盾。

fig, ax = plt.subplots(figsize=(9, 4.5))
ax.bar(
    beta_portfolios_summary["portfolio"].astype(str),
    beta_portfolios_summary["alpha"],
    color="#4C78A8",
    edgecolor="white",
)
ax.axhline(0, color="black", linewidth=0.8)
ax.set_xlabel("Portfolio (1=low beta … 10=high beta)")
ax.set_ylabel("CAPM alpha (monthly)")
ax.set_title("CAPM alphas of beta-sorted portfolios (A-share)")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
plt.tight_layout
plt.show

6.8.6 证券市场线（SML）与 Beta 组合

对应原文 The Security Market Line and Beta Portfolios。CAPM 预测各组合应落在 SML 上，斜率等于市场风险溢价。下图同时画出：

1. 以平均市场超额收益为斜率、过原点的理论 SML（实线）；
2. 用十分组的（组合 beta，平均收益）拟合的经验直线（虚线）。

# 理论 SML：过原点，斜率为样本期平均市场超额收益
mkt_premium = float(factors_mkt_monthly["mkt_excess"].mean)

# 经验直线：ret ~ a + b * beta
X = sm.add_constant(beta_portfolios_summary["beta"])
sml_fit = sm.OLS(beta_portfolios_summary["ret"], X).fit
intercept_hat = float(sml_fit.params["const"])
slope_hat = float(sml_fit.params["beta"])

fig, ax = plt.subplots(figsize=(7.5, 5.5))
ax.scatter(
    beta_portfolios_summary["beta"],
    beta_portfolios_summary["ret"],
    c=np.arange(len(beta_portfolios_summary)),
    cmap="viridis",
    s=60,
    zorder=3,
)
for _, r in beta_portfolios_summary.iterrows:
    ax.annotate(
        str(r["portfolio"]),
        (r["beta"], r["ret"]),
        textcoords="offset points",
        xytext=(5, 4),
        fontsize=9,
    )

x_line = np.linspace(0, max(2.0, beta_portfolios_summary["beta"].max * 1.05), 100)
ax.plot(x_line, mkt_premium * x_line, color="black", linewidth=1.5, label="Theoretical SML")
ax.plot(
    x_line,
    intercept_hat + slope_hat * x_line,
    color="crimson",
    linestyle="--",
    linewidth=1.5,
    label="Fitted line on portfolios",
)
ax.set_xlim(0, max(2.0, beta_portfolios_summary["beta"].max * 1.05))
ymax = max(mkt_premium * 2, beta_portfolios_summary["ret"].max * 1.2, 0.02)
ax.set_ylim(min(0, beta_portfolios_summary["ret"].min * 1.2), ymax)
ax.set_xlabel("Beta")
ax.set_ylabel("Average excess return")
ax.set_title("Average portfolio excess returns and beta estimates (A-share)")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
ax.legend(frameon=False)
ax.grid(True, alpha=0.3)
plt.tight_layout
plt.show

print(f"market premium (mean mkt_excess) = {mkt_premium:.4%}")
print(f"fitted intercept = {intercept_hat:.4%}, slope = {slope_hat:.4f}")

market premium (mean mkt_excess) = 0.6856%
fitted intercept = 1.4026%, slope = -0.0050

再构造十分组下的极端多空：做多最高 beta 组、做空最低 beta 组，并做 Newey–West 均值检验与 CAPM alpha 检验。

wide = (
    beta_portfolios.pivot(index="date", columns="portfolio", values="ret")
    .sort_index
)
# 组合标签为 '1'..'10'
low_col, high_col = "1", "10"
beta_longshort = wide[[low_col, high_col]].dropna.copy
beta_longshort = beta_longshort.rename(columns={low_col: "low", high_col: "high"})
beta_longshort["long_short"] = beta_longshort["high"] - beta_longshort["low"]
beta_longshort = beta_longshort.reset_index.merge(
    factors_mkt_monthly, on="date", how="left"
)

print("long-short months:", len(beta_longshort))
display(beta_longshort.head)

# 平均收益是否为 0
y = beta_longshort["long_short"]
model_mu = sm.OLS(y, np.ones((len(y), 1))).fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("\nLong-short ~ 1  (NW lags=6)")
print(model_mu.summary)

# CAPM alpha
X = sm.add_constant(beta_longshort["mkt_excess"])
model_capm = sm.OLS(beta_longshort["long_short"], X, missing="drop").fit(
    cov_type="HAC", cov_kwds={"maxlags": 6}
)
print("\nLong-short ~ 1 + mkt_excess  (NW lags=6)")
print(model_capm.summary)

long-short months: 271

	date	low	high	long_short	mkt_excess
0	2004-01-01	0.024797	0.138344	0.113547	0.073093
1	2004-02-01	0.041176	0.122426	0.081249	0.069838
2	2004-03-01	0.008092	0.012647	0.004555	0.030699
3	2004-04-01	-0.144966	-0.124091	0.020875	-0.100944
4	2004-05-01	-0.036202	-0.007386	0.028816	-0.024572

Long-short ~ 1  (NW lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                       nan
Date:                Sat, 08 Aug 2026   Prob (F-statistic):                nan
Time:                        21:15:13   Log-Likelihood:                 324.30
No. Observations:                 271   AIC:                            -646.6
Df Residuals:                     270   BIC:                            -643.0
Df Model:                           0                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0050      0.004     -1.419      0.156      -0.012       0.002
==============================================================================
Omnibus:                       29.855   Durbin-Watson:                   1.980
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              152.214
Skew:                          -0.128   Prob(JB):                     8.85e-34
Kurtosis:                       6.663   Cond. No.                         1.00
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

Long-short ~ 1 + mkt_excess  (NW lags=6)
                            OLS Regression Results                            
==============================================================================
Dep. Variable:             long_short   R-squared:                       0.346
Model:                            OLS   Adj. R-squared:                  0.343
Method:                 Least Squares   F-statistic:                     63.87
Date:                Sat, 08 Aug 2026   Prob (F-statistic):           3.93e-14
Time:                        21:15:13   Log-Likelihood:                 381.75
No. Observations:                 271   AIC:                            -759.5
Df Residuals:                     269   BIC:                            -752.3
Df Model:                           1                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0097      0.004     -2.618      0.009      -0.017      -0.002
mkt_excess     0.5629      0.070      7.992      0.000       0.425       0.701
==============================================================================
Omnibus:                       60.286   Durbin-Watson:                   1.841
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              262.919
Skew:                          -0.826   Prob(JB):                     8.09e-58
Kurtosis:                       7.534   Cond. No.                         13.1
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and without small sample correction

若多空组合的平均收益不显著，但在控制市场后出现显著为负的 CAPM alpha，则与 CAPM 预测不一致，并与文献中的 betting-against-beta（做多低 beta、做空高 beta）方向一致 [@Frazzini2014]。需注意：Frazzini–Pedersen 的因子构造并非标准十分组市值加权，其异常收益是否稳健仍有争议 [@NovyMarx2022]。

6.8.7 年度收益分解

对应原文年度柱状图：观察低 beta、高 beta 与多空组合在各年的累计收益，检查是否存在少数年份主导全样本结论。

tmp = beta_longshort.copy
tmp["year"] = tmp["date"].dt.year

def ann_comp(s: pd.Series) -> float:
    return float(np.prod(1.0 + s.dropna) - 1.0)

annual = (
    tmp.groupby("year")
    .agg(
        low=("low", ann_comp),
        high=("high", ann_comp),
        long_short=("long_short", ann_comp),
    )
    .reset_index
)

plot_df = annual.melt(
    id_vars="year",
    value_vars=["low", "high", "long_short"],
    var_name="name",
    value_name="value",
)

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
for ax, name in zip(axes, ["low", "high", "long_short"]):
    sub = plot_df[plot_df["name"] == name]
    ax.bar(sub["year"], sub["value"], color="#4C78A8", width=0.8)
    ax.axhline(0, color="black", linewidth=0.8)
    ax.set_ylabel(name)
    ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
    ax.grid(True, axis="y", alpha=0.3)
axes[0].set_title("Annual returns of beta portfolios (A-share)")
axes[-1].set_xlabel("")
plt.tight_layout
plt.show

display(annual.tail(10))

	year	low	high	long_short
13	2017	0.287904	-0.157718	-0.353551
14	2018	-0.169761	-0.321624	-0.194517
15	2019	0.242932	0.247311	0.012042
16	2020	0.097370	0.157694	0.056586
17	2021	0.018111	0.079606	0.055313
18	2022	-0.081621	-0.204390	-0.125512
19	2023	0.121660	-0.132788	-0.236273
20	2024	0.281498	0.026509	-0.206064
21	2025	0.093440	0.269349	0.149724
22	2026	0.010071	-0.081584	-0.173771

6.8.8 要点小结

• 单维组合排序用于检验单一特征（如滞后市场 beta）能否预测未来超额收益。
• 每月按分位数断点分组，并用市值加权计算组合收益，更接近可投资策略。
• 以 beta 排序的多空组合往往难以稳定获得显著为正的平均超额收益，这与“高 beta 应有高期望收益”的 CAPM 直觉冲突。
• 低 beta 组合常呈现相对更高的 CAPM alpha，高 beta 组合相反，构成 betting-against-beta / 低 beta 异象的证据线索。
• 把 assign_portfolio 与市值加权收益写成函数后，可把同一套流程迁移到任意排序变量（价值、动量、质量等）。

练习（对应原文 Exercises，改为 A 股）

1. 比较中位数两组多空与十分组高低多空的收益差异：平均收益、波动、Sharpe，以及 Newey–West $t$ 值是否明显不同？
2. 对十分组的 CAPM alpha 逐个做显著性检验（可用 HAC）。哪些组合的 alpha 显著？
3. 把排序变量换成第 2.6 节的日频 beta（beta_ashare.parquet 中 return_type=="daily"），重复本节分析，指出与月频 beta 结果的差异。
4. 基于本节结果，设计一个有望获得正 alpha 的多空方向（提示：与“多高空低”相反），绘制其累计超额收益，并与市场超额收益对比。

6.9 A股复现：规模排序与 p-hacking（Size Sorts and p-Hacking）

本节复现 tidyfinance 章节 Size Sorts and p-Hacking：在单维排序框架下，把排序变量换成公司规模（流通市值），得到经典的规模溢价（做多小市值、做空大市值）；并展示断点样本、加权方式、分组数、样本期等研究设计选择如何改变溢价估计——这正是排序情境下的 non-standard errors / p-hacking 风险。

import warnings
from itertools import product
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("monthly cache exists =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)

monthly cache exists = True

6.9.1 数据准备

对应原文 Data Preparation。规模用上期流通市值 mktcap_lag（建仓时已知）；收益用 ret_excess。板块由 permno 规则映射，作为美股 exchange 的 A 股替代。

def map_board(permno: str) -> str:
    # 由证券代码推断上市板块（简化规则）
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "科创板"
    if code.startswith(("300", "301")):
        return "创业板"
    if code.startswith("6"):
        return "上交所主板"
    if code.startswith(("000", "001", "002", "003")):
        return "深交所主板"
    return "其他"

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)

# 排序与加权需要上期市值；缺失则当月无法纳入
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

print(crsp_monthly.shape)
print(crsp_monthly["board"].value_counts)
display(crsp_monthly.head(3))

(884560, 11)
board
上交所主板    344321
深交所主板    318516
创业板      143454
其他        44021
科创板       34248
Name: count, dtype: int64

	permno	month	ret	mktcap	ret_m	rf	date	ret_excess	mkt_excess	mktcap_lag	board
1	000001	2000-02	-0.011333	19618933.36	0.122024	0.001856	2000-02-01	-0.013189	0.120168	19843822.88	深交所主板
2	000001	2000-03	0.002729	19672478.48	0.055926	0.001856	2000-03-01	0.000873	0.054070	19618933.36	深交所主板
3	000001	2000-04	0.037017	20400692.17	0.013014	0.001856	2000-04-01	0.035161	0.011158	19672478.48	深交所主板

6.9.2 规模分布

对应原文 Size Distribution。先看市值有多集中：最大 1%/5%/10%/25% 公司占全市场流通市值的比例。若头部集中度很高，则市值加权组合会强烈被大票主导。

def top_cap_share(g: pd.DataFrame, q: float) -> float:
    thr = g["mktcap"].quantile(q)
    total = g["mktcap"].sum
    if total <= 0:
        return np.nan
    return float(g.loc[g["mktcap"] >= thr, "mktcap"].sum / total)

rows = []
for dt, g in crsp_monthly.groupby("date", sort=True):
    rows.append(
        {
            "date": dt,
            "Largest 1%": top_cap_share(g, 0.99),
            "Largest 5%": top_cap_share(g, 0.95),
            "Largest 10%": top_cap_share(g, 0.90),
            "Largest 25%": top_cap_share(g, 0.75),
        }
    )
conc = pd.DataFrame(rows).sort_values("date")
plot_df = conc.melt(id_vars="date", var_name="name", value_name="value")

fig, ax = plt.subplots(figsize=(10, 4.5))
for name, sub in plot_df.groupby("name"):
    ax.plot(sub["date"], sub["value"], label=name, linewidth=1.6)
ax.set_title("最大市值公司占总流通市值比例（A股）")
ax.set_ylabel("份额")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
ax.legend(frameon=False, ncol=2)
ax.grid(True, alpha=0.3)
plt.tight_layout
plt.show

display(conc.tail(3))

	date	Largest 1%	Largest 5%	Largest 10%	Largest 25%
315	2026-05-01	0.300931	0.529262	0.648327	0.809577
316	2026-06-01	0.300993	0.539631	0.662786	0.823773
317	2026-07-01	0.320657	0.552181	0.668810	0.822761

再看各上市板块的市值份额随时间变化（类比原文按 NYSE/NASDAQ/AMEX 堆叠面积图）。

share = (
    crsp_monthly.groupby(["date", "board"], as_index=False)["mktcap"]
    .sum
    .rename(columns={"mktcap": "board_cap"})
)
share["total"] = share.groupby("date")["board_cap"].transform("sum")
share["share"] = share["board_cap"] / share["total"]

pivot = share.pivot(index="date", columns="board", values="share").fillna(0).sort_index
# 固定堆叠顺序
cols = [c for c in ["上交所主板", "深交所主板", "创业板", "科创板", "其他"] if c in pivot.columns]
pivot = pivot[cols]

fig, ax = plt.subplots(figsize=(10, 4.5))
ax.stackplot(pivot.index, [pivot[c].values for c in cols], labels=cols, alpha=0.85)
ax.set_title("各上市板块占总流通市值份额（A股）")
ax.set_ylabel("份额")
ax.set_ylim(0, 1)
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.0f}%"))
ax.legend(loc="upper left", ncol=3, frameon=False, fontsize=9)
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout
plt.show

最新月份各板块市值的描述统计：主板公司平均规模通常更大；创业板/科创板公司数量多但体量更分散。这也是“只用主板样本算断点”会影响小市值组厚度的原因。

latest = crsp_monthly["date"].max
snap = crsp_monthly[crsp_monthly["date"] == latest].copy

def board_summary(df: pd.DataFrame) -> pd.DataFrame:
    def one(g):
        s = g["mktcap"]
        return pd.Series(
            {
                "count": len(s),
                "mean": s.mean,
                "std": s.std,
                "min": s.min,
                "p5": s.quantile(0.05),
                "p50": s.quantile(0.50),
                "p95": s.quantile(0.95),
                "max": s.max,
            }
        )

    by = df.groupby("board").apply(one, include_groups=False)
    overall = one(df)
    overall.name = "Overall"
    out = pd.concat([by, overall.to_frame.T])
    return out

print("snapshot month:", latest.date)
summary = board_summary(snap)
display(summary.round(0))

snapshot month: 2026-07-01

	count	mean	std	min	p5	p50	p95	max
上交所主板	1699.0	29093131.0	119091033.0	722864.0	1825229.0	6953583.0	104098370.0	2.202785e+09
其他	396.0	1315870.0	2552979.0	16781.0	78636.0	716231.0	4035889.0	3.429679e+07
创业板	1397.0	10177802.0	56703424.0	24043.0	1193492.0	3487125.0	26520174.0	1.684010e+09
深交所主板	1494.0	14654251.0	35785984.0	37467.0	1543550.0	5414599.0	56215101.0	6.028877e+08
科创板	608.0	16589284.0	46424082.0	209392.0	1581280.0	5953805.0	53572655.0	6.948920e+08
Overall	5594.0	17187780.0	75890924.0	16781.0	861935.0	4741815.0	56425655.0	2.202785e+09

6.9.3 灵活断点的单维规模排序

对应原文 Univariate Size Portfolios with Flexible Breakpoints。关键扩展：断点可以只在子集上计算（如仅上交所主板），再应用到当月全部股票——类似美股“NYSE breakpoints”。

def assign_portfolio(
    data: pd.DataFrame,
    boards: list[str],
    sorting_variable: str,
    n_portfolios: int,
) -> pd.Series:
    # 用 boards 子集计算分位断点，再映射到全样本
    subset = data[data["board"].isin(boards)]
    if len(subset) < n_portfolios:
        # 退化：全体同一组
        return pd.Series(1, index=data.index, dtype=int)

    quantiles = np.linspace(0, 1, n_portfolios + 1)
    breakpoints = subset[sorting_variable].quantile(quantiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    # cut 需要内部断点；两端扩展为 -inf/inf
    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

# 示例：某月用上交所主板断点 vs 全市场断点，比较落入“小市值组”的股票数
demo_date = crsp_monthly["date"].max
g = crsp_monthly[crsp_monthly["date"] == demo_date].copy
g["port_all"] = assign_portfolio(g, g["board"].unique.tolist, "mktcap_lag", 2)
g["port_sh"] = assign_portfolio(g, ["上交所主板"], "mktcap_lag", 2)
print("month:", demo_date.date)
print("small-group count | all-board breakpoints:", int((g["port_all"] == 1).sum))
print("small-group count | SH-main breakpoints:", int((g["port_sh"] == 1).sum))

month: 2026-07-01
small-group count | all-board breakpoints: 2797
small-group count | SH-main breakpoints: 3408

6.9.4 组合加权：市值加权 vs 等权

对应原文 Weighting Schemes。市值加权更接近被动可投资组合；等权需每月再平衡，实务成本更高，但常放大小市值溢价。规模溢价定义为：最小市值组收益 − 最大市值组收益（small minus big）。

def compute_portfolio_returns(
    n_portfolios: int = 10,
    boards: list[str] | None = None,
    value_weighted: bool = True,
    data: pd.DataFrame | None = None,
) -> float:
    # 返回全样本平均规模溢价（每月 small-big 的时间序列均值）
    if data is None:
        data = crsp_monthly
    if boards is None:
        boards = sorted(data["board"].dropna.unique.tolist)

    premia = []
    for dt, g in data.groupby("date", sort=True):
        g = g.copy
        g["portfolio"] = assign_portfolio(g, boards, "mktcap_lag", n_portfolios)
        # 各组收益
        rets = {}
        for port, pg in g.groupby("portfolio"):
            if value_weighted:
                w = pg["mktcap_lag"]
                if w.sum <= 0:
                    continue
                rets[int(port)] = float((pg["ret_excess"] * w).sum / w.sum)
            else:
                rets[int(port)] = float(pg["ret_excess"].mean)
        if len(rets) < 2:
            continue
        pmin, pmax = min(rets), max(rets)
        premia.append(rets[pmin] - rets[pmax])

    if not premia:
        return np.nan
    return float(np.mean(premia))

boards_all = ["上交所主板", "深交所主板", "创业板", "科创板"]
ret_all = compute_portfolio_returns(
    n_portfolios=2, boards=boards_all, value_weighted=True, data=crsp_monthly
)
ret_sh = compute_portfolio_returns(
    n_portfolios=2, boards=["上交所主板"], value_weighted=True, data=crsp_monthly
)

cmp = pd.DataFrame(
    {
        "断点样本": ["全部主要板块", "仅上交所主板"],
        "平均规模溢价(月)": [ret_all, ret_sh],
        "折年约%": [ret_all * 12 * 100, ret_sh * 12 * 100],
    }
)
display(cmp.round(4))
print("说明：仅用上交所主板断点时，全市场小市值组通常更‘厚’，溢价估计往往变化明显。")

	断点样本	平均规模溢价(月)	折年约%
0	全部主要板块	0.0081	9.6901
1	仅上交所主板	0.0068	8.1625

说明：仅用上交所主板断点时，全市场小市值组通常更‘厚’，溢价估计往往变化明显。

6.9.5 p-hacking 与 non-standard errors

对应原文 P-Hacking and Non-Standard Errors。排序至少要选：分组数、断点样本、等权/市值加权，以及是否截断样本期或剔除某些板块。这些选择都可以在顶级期刊文献中找到先例；它们不一定“错”，但会带来 non-standard errors（研究者选择带来的变异）。若只报告最显著的一种设定，就滑向 p-hacking。

下面用设定网格扫描多种组合，观察规模溢价分布。为控制运行时间，网格比原文略收敛，但仍覆盖关键决策节点。

n_portfolios_grid = [2, 5, 10]
boards_grid = [
    ["上交所主板"],
    ["上交所主板", "深交所主板", "创业板", "科创板"],
]
value_weighted_grid = [True, False]

# 数据切片：全样本 / 剔除双创（只留主板） / 2010前 / 2010及以后
data_main = crsp_monthly[crsp_monthly["board"].isin(["上交所主板", "深交所主板"])].copy
data_pre = crsp_monthly[crsp_monthly["date"] < "2010-01-01"].copy
data_post = crsp_monthly[crsp_monthly["date"] >= "2010-01-01"].copy
data_grid = [
    ("全样本", crsp_monthly),
    ("仅主板", data_main),
    ("2010前", data_pre),
    ("2010及以后", data_post),
]

setup = list(product(n_portfolios_grid, boards_grid, value_weighted_grid, data_grid))
print("specifications:", len(setup))

records = []
for n_port, boards, vw, (tag, df) in setup:
    if len(df) < 1000 or df["date"].nunique < 24:
        prem = np.nan
    else:
        prem = compute_portfolio_returns(
            n_portfolios=n_port,
            boards=boards,
            value_weighted=vw,
            data=df,
        )
    records.append(
        {
            "n_portfolios": n_port,
            "breakpoint_boards": "+".join(boards) if len(boards) == 1 else "多板块",
            "value_weighted": vw,
            "sample": tag,
            "size_premium": prem,
        }
    )

p_hacking_results = pd.DataFrame(records).dropna(subset=["size_premium"])
p_hacking_results = p_hacking_results.sort_values("size_premium", ascending=False)
print(p_hacking_results.head(10))
print("...")
print(p_hacking_results.tail(5))
print(
    "mean=",
    round(p_hacking_results["size_premium"].mean, 5),
    "std=",
    round(p_hacking_results["size_premium"].std, 5),
    "min=",
    round(p_hacking_results["size_premium"].min, 5),
    "max=",
    round(p_hacking_results["size_premium"].max, 5),
)

specifications: 48
    n_portfolios breakpoint_boards  value_weighted sample  size_premium
38            10             上交所主板           False  2010前      0.021159
46            10               多板块           False  2010前      0.020778
34            10             上交所主板            True  2010前      0.017018
45            10               多板块           False    仅主板      0.016913
42            10               多板块            True  2010前      0.016659
44            10               多板块           False    全样本      0.016365
36            10             上交所主板           False    全样本      0.015902
30             5               多板块           False  2010前      0.015687
37            10             上交所主板           False    仅主板      0.015563
22             5             上交所主板           False  2010前      0.014876
...
    n_portfolios breakpoint_boards  value_weighted   sample  size_premium
10             2               多板块            True    2010前      0.007132
3              2             上交所主板            True  2010及以后      0.006864
0              2             上交所主板            True      全样本      0.006802
2              2             上交所主板            True    2010前      0.006699
1              2             上交所主板            True      仅主板      0.006532
mean= 0.01224 std= 0.00363 min= 0.00653 max= 0.02116

6.9.6 规模溢价的设定变异

对应原文 Size-Premium Variation。直方图展示不同设定下的月均规模溢价；竖虚线为一种“基准设定”的溢价（十分组、上交所主板断点、市值加权、全样本），扮演原文中 FF-SMB 均值参照线的角色。

# 基准设定（类比文献常用口径，并非唯一正确）
benchmark_premium = compute_portfolio_returns(
    n_portfolios=10,
    boards=["上交所主板"],
    value_weighted=True,
    data=crsp_monthly,
)

fig, ax = plt.subplots(figsize=(9, 4.5))
vals = p_hacking_results["size_premium"].values
ax.hist(vals, bins=min(40, max(10, len(vals) // 2)), color="#4C78A8", edgecolor="white")
ax.axvline(benchmark_premium, color="crimson", linestyle="--", linewidth=1.8, label="基准设定")
ax.set_title("不同排序设定下的规模溢价分布（A股）")
ax.set_xlabel("平均月度规模溢价（small - big）")
ax.set_ylabel("设定个数")
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{100*x:.2f}%"))
ax.legend(frameon=False)
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout
plt.show

print(f"benchmark monthly size premium = {benchmark_premium:.4%}  (ann. ~ {benchmark_premium*12:.2%})")

# 哪个决策节点影响最大？看分组后的溢价离散度
impact = []
for col in ["n_portfolios", "breakpoint_boards", "value_weighted", "sample"]:
    grp = p_hacking_results.groupby(col)["size_premium"].mean
    impact.append({"choice": col, "range_across_levels": float(grp.max - grp.min)})
impact_df = pd.DataFrame(impact).sort_values("range_across_levels", ascending=False)
print("\n各决策节点：不同水平下平均溢价的极差（越大说明该选择越能左右结果）")
display(impact_df)

benchmark monthly size premium = 1.4274%  (ann. ~ 17.13%)

各决策节点：不同水平下平均溢价的极差（越大说明该选择越能左右结果）

	choice	range_across_levels
0	n_portfolios	0.007701
3	sample	0.002164
2	value_weighted	0.001620
1	breakpoint_boards	0.000731

6.9.7 要点小结

• 规模是重要排序变量；小市值相对大市值的多空收益即规模溢价，也是 FF 三因子中 SMB 的经验来源之一。
• 组合收益对断点样本、加权方式、分组数、样本期/板块过滤等设计选择敏感。
• 这种因合理备选流程而产生的变异，属于 non-standard errors；只报告最显著设定则接近 p-hacking。
• 稳健做法：预先登记主设定，并报告多种合理设定下的结果区间（见第4章 4.1.4）。

练习（对应原文 Exercises，改为 A 股）

1. 画出上交所主板、深交所主板、创业板的平均市值随时间变化，解读结构变迁。
2. 画出十分组规模排序的断点（或极端组断点差）时序，检查是否存在异常跳变。
3. 把 compute_portfolio_returns 改成输出多空组合的 CAPM alpha（相对 mkt_excess），比较与平均超额收益结论是否一致。
4. 在设定网格中，用回归或分组统计找出对溢价影响最大的单个决策变量。
5. 哪一种设定最接近“主板断点 + 市值加权 + 多空 small-big”的常见 SMB 口径？它在分布中处于什么位置？

6.10 A股复现：价值与双维排序（Value and Bivariate Sorts）

本节复现 tidyfinance 章节 Value and Bivariate Sorts：在单维排序之上，引入账面市值比（BM），并完成 Size × BM 的独立双维排序与依赖（条件）双维排序；最后比较两种协议下 5×5 组合的股票数与市值份额分布。

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"
BE_CACHE = CACHE_DIR / "book_equity_ashare.parquet"

print("crsp cache =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE cache   =", BE_CACHE.exists)

crsp cache = True
BE cache   = True

6.10.1 数据准备

对应原文 Data Preparation。月度收益/市值来自既有面板；账面权益若本地无缓存，则按年报截止日（YYYY1231）调用 akshare.stock_zcfz_em 拉取并写入 parquet（首次较慢，之后直接读缓存）。

def map_board(permno: str) -> str:
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "科创板"
    if code.startswith(("300", "301")):
        return "创业板"
    if code.startswith("6"):
        return "上交所主板"
    if code.startswith(("000", "001", "002", "003")):
        return "深交所主板"
    return "其他"

def load_or_download_book_equity(
    cache_path: Path,
    years: range | list[int] = range(1999, 2026),
) -> pd.DataFrame:
    """Load cached annual book equity; download via akshare if missing."""
    if cache_path.exists:
        be = pd.read_parquet(cache_path)
        be["permno"] = be["permno"].astype(str).str.zfill(6)
        be["datadate"] = pd.to_datetime(be["datadate"])
        return be

    import time
    import akshare as ak

    print("缓存不存在，开始下载年报股东权益（首次约十余分钟）…")
    frames: list[pd.DataFrame] = []
    partial = cache_path.with_suffix(".partial.parquet")
    done: set[int] = set
    if partial.exists:
        old = pd.read_parquet(partial)
        frames.append(old)
        done = set(pd.to_datetime(old["datadate"]).dt.year.unique.tolist)
        print("resume years:", sorted(done))

    for y in years:
        if y in done:
            continue
        date = f"{y}1231"
        print("fetch", date)
        try:
            raw = ak.stock_zcfz_em(date=date)
        except Exception as exc:  # noqa: BLE001
            print("  FAIL", date, exc)
            time.sleep(2)
            continue
        cols = list(raw.columns)
        code_col = [c for c in cols if "代码" in str(c)][0]
        be_col = [c for c in cols if "股东权益" in str(c)][0]
        out = pd.DataFrame(
            {
                "permno": raw[code_col].astype(str).str.zfill(6),
                "datadate": pd.to_datetime(date),
                "be": pd.to_numeric(raw[be_col], errors="coerce"),
            }
        )
        out = out.dropna(subset=["be"])
        out = out[out["be"] > 0].copy
        # 东财单位：元；CSMAR Msmvosd：千元
        out["be"] = out["be"] / 1000.0
        frames.append(out)
        done.add(y)
        pd.concat(frames, ignore_index=True).to_parquet(partial, index=False)
        time.sleep(0.5)

    be = pd.concat(frames, ignore_index=True)
    be = be.drop_duplicates(["permno", "datadate"], keep="last")
    cache_path.parent.mkdir(parents=True, exist_ok=True)
    be.to_parquet(cache_path, index=False)
    if partial.exists:
        partial.unlink
    print("saved", cache_path, be.shape)
    return be

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

book_equity = load_or_download_book_equity(BE_CACHE)
book_equity = book_equity[book_equity["be"] > 0].dropna(subset=["be"]).copy
# 与原文一致：会计日对齐到月（年报用 12-01 表示该月）
book_equity["date"] = book_equity["datadate"].dt.to_period("M").dt.to_timestamp

print("crsp:", crsp_monthly.shape, "BE rows:", len(book_equity),
      "years:", book_equity["datadate"].dt.year.min, "-", book_equity["datadate"].dt.year.max)
display(book_equity.head(3))

crsp: (884560, 11) BE rows: 88343 years: 1999 - 2025

	permno	datadate	be	date
0	601963	1999-12-31	4.083935e+05	1999-12-01
1	601229	1999-12-31	4.182000e+06	1999-12-01
2	601187	1999-12-31	3.140239e+05	1999-12-01

6.10.2 账面市值比（Book-to-Market）

对应原文 Book-to-Market Ratio。核心是避免前视偏差：

1. 规模滞后 1 个月（sorting_date = date + 1M）；
2. 在财报月把 be 与当月 mktcap 匹配得到 BM，再整体滞后 6 个月；
3. 将 BM 按股票向前填充，并丢弃会计信息超过 12 个月的观测。

高 BM 为价值股，低 BM 为成长股。

size = (
    crsp_monthly.assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=1))
    .rename(columns={"mktcap": "size"})
    [["permno", "sorting_date", "size"]]
)

bm = (
    book_equity.merge(
        crsp_monthly[["permno", "date", "mktcap"]],
        on=["permno", "date"],
        how="inner",
    )
    .assign(
        bm=lambda d: d["be"] / d["mktcap"],
        sorting_date=lambda d: d["date"] + pd.DateOffset(months=6),
        accounting_date=lambda d: d["date"] + pd.DateOffset(months=6),
    )
    [["permno", "sorting_date", "accounting_date", "bm"]]
)
bm = bm.replace([np.inf, -np.inf], np.nan).dropna(subset=["bm"])
bm = bm[bm["bm"] > 0].copy

data_for_sorts = (
    crsp_monthly.merge(
        bm,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
        suffixes=("", "_bm"),
    )
    .drop(columns=["sorting_date"], errors="ignore")
    .merge(
        size,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
    )
    .drop(columns=["sorting_date"], errors="ignore")
)

data_for_sorts = data_for_sorts.sort_values(["permno", "date"]).copy
data_for_sorts["bm"] = data_for_sorts.groupby("permno")["bm"].ffill
data_for_sorts["accounting_date"] = data_for_sorts.groupby("permno")["accounting_date"].ffill
data_for_sorts["threshold_date"] = data_for_sorts["date"] - pd.DateOffset(months=12)
data_for_sorts = data_for_sorts[
    data_for_sorts["accounting_date"] > data_for_sorts["threshold_date"]
].copy
data_for_sorts = data_for_sorts.dropna(
    subset=["ret_excess", "mktcap_lag", "size", "bm", "board"]
).copy
data_for_sorts = data_for_sorts.drop(columns=["threshold_date", "accounting_date"])

print(data_for_sorts.shape)
print(
    "date range:",
    data_for_sorts["date"].min.date,
    "→",
    data_for_sorts["date"].max.date,
)
display(data_for_sorts[["permno", "date", "size", "bm", "ret_excess", "board"]].head(3))

(703738, 13)
date range: 2001-06-01 → 2026-07-01

	permno	date	size	bm	ret_excess	board
16	000001	2001-06-01	22568621.18	0.173894	-0.056794	深交所主板
17	000001	2001-07-01	21328740.14	0.173894	-0.098525	深交所主板
18	000001	2001-08-01	19266915.49	0.173894	-0.082839	深交所主板

断点函数与 6.9 相同：在指定板块子集上算分位，再映射到当月全样本。

def assign_portfolio(
    data: pd.DataFrame,
    boards: list[str],
    sorting_variable: str,
    n_portfolios: int,
) -> pd.Series:
    subset = data[data["board"].isin(boards)]
    if len(subset) < n_portfolios:
        return pd.Series(1, index=data.index, dtype=int)

    quantiles = np.linspace(0, 1, n_portfolios + 1)
    breakpoints = subset[sorting_variable].quantile(quantiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

BREAKPOINT_BOARDS = ["上交所主板"]
print("default breakpoint boards:", BREAKPOINT_BOARDS)

default breakpoint boards: ['上交所主板']

6.10.3 独立双维排序（Independent Sorts）

对应原文 Independent Sorts。每月分别按 Size、BM 用主板断点各分 5 组，交叉得到 25 个组合；组内按 mktcap_lag 市值加权。价值溢价 = 高 BM 五组等权平均 − 低 BM 五组等权平均（在规模维上先等权合成）。

def monthly_independent_returns(df: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_bm"] = assign_portfolio(g, BREAKPOINT_BOARDS, "bm", 5)
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        for (p_bm, p_sz), pg in g.groupby(["portfolio_bm", "portfolio_size"]):
            w = pg["mktcap_lag"]
            if w.sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio_bm": int(p_bm),
                    "portfolio_size": int(p_sz),
                    "ret": float((pg["ret_excess"] * w).sum / w.sum),
                }
            )
    return pd.DataFrame(rows)

value_portfolios_ind = monthly_independent_returns(data_for_sorts)

# 每月：先在各 BM 组内对 5 个规模组合等权，再 high-low
vp = (
    value_portfolios_ind.groupby(["date", "portfolio_bm"], as_index=False)["ret"]
    .mean
)
prem_ind = []
for dt, g in vp.groupby("date"):
    hi = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    lo = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    prem_ind.append(hi - lo)
value_premium_ind = float(np.mean(prem_ind))

print(
    f"independent value premium (monthly) = {value_premium_ind:.4%}  "
    f"(ann. ~ {value_premium_ind * 12:.2%})"
)
display(value_portfolios_ind.head(3))

independent value premium (monthly) = 0.1419%  (ann. ~ 1.70%)

	date	portfolio_bm	portfolio_size	ret
0	2001-06-01	1	1	0.022497
1	2001-06-01	1	2	-0.014334
2	2001-06-01	1	3	-0.009245

6.10.4 依赖双维排序（Dependent Sorts）

对应原文 Dependent Sorts。先按 Size 分 5 组，再在每个规模组内部按 BM 分 5 组——BM 断点随规模组变化。依赖排序在“断点取自全市场”时往往使各组股票数更均匀；若断点仍只用主板，组内数量仍可能不均，但 BM 维的相对均匀性通常优于独立排序。

def monthly_dependent_returns(df: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        parts = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, BREAKPOINT_BOARDS, "bm", 5)
            parts.append(sg)
        g2 = pd.concat(parts, ignore_index=False)
        for (p_bm, p_sz), pg in g2.groupby(["portfolio_bm", "portfolio_size"]):
            w = pg["mktcap_lag"]
            if w.sum <= 0:
                continue
            rows.append(
                {
                    "date": dt,
                    "portfolio_bm": int(p_bm),
                    "portfolio_size": int(p_sz),
                    "ret": float((pg["ret_excess"] * w).sum / w.sum),
                }
            )
    return pd.DataFrame(rows)

value_portfolios_dep = monthly_dependent_returns(data_for_sorts)
vp_d = (
    value_portfolios_dep.groupby(["date", "portfolio_bm"], as_index=False)["ret"]
    .mean
)
prem_dep = []
for dt, g in vp_d.groupby("date"):
    hi = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    lo = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    prem_dep.append(hi - lo)
value_premium_dep = float(np.mean(prem_dep))

cmp = pd.DataFrame(
    {
        "排序方式": ["独立排序", "依赖排序"],
        "月均价值溢价": [value_premium_ind, value_premium_dep],
        "折年约%": [value_premium_ind * 12 * 100, value_premium_dep * 12 * 100],
    }
)
display(cmp.round(4))

	排序方式	月均价值溢价	折年约%
0	独立排序	0.0014	1.7027
1	依赖排序	0.0013	1.5334

6.10.5 组合构成：股票数与市值份额

对应原文 Portfolio Composition。收益之外，看 5×5 网格里股票如何分布、市值如何集中——独立 vs 依赖排序的差异主要体现在 BM 维是否“挤在角落”。

def assign_independent(df: pd.DataFrame) -> pd.DataFrame:
    parts = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        g["portfolio_bm"] = assign_portfolio(g, BREAKPOINT_BOARDS, "bm", 5)
        g["sorting_method"] = "Independent"
        parts.append(g)
    return pd.concat(parts, ignore_index=True)

def assign_dependent(df: pd.DataFrame) -> pd.DataFrame:
    parts = []
    for dt, g in df.groupby("date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, BREAKPOINT_BOARDS, "size", 5)
        sub = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, BREAKPOINT_BOARDS, "bm", 5)
            sub.append(sg)
        g2 = pd.concat(sub, ignore_index=False)
        g2["sorting_method"] = "Dependent"
        parts.append(g2)
    return pd.concat(parts, ignore_index=True)

assignments = pd.concat(
    [assign_independent(data_for_sorts), assign_dependent(data_for_sorts)],
    ignore_index=True,
)

monthly_char = (
    assignments.groupby(
        ["sorting_method", "date", "portfolio_size", "portfolio_bm"], as_index=False
    )
    .agg(n_stocks=("permno", "size"), mktcap=("mktcap_lag", "sum"))
)

portfolio_characteristics = (
    monthly_char.groupby(
        ["sorting_method", "portfolio_size", "portfolio_bm"], as_index=False
    )
    .agg(n_stocks=("n_stocks", "mean"), mktcap=("mktcap", "mean"))
)
portfolio_characteristics["mktcap_share"] = portfolio_characteristics[
    "mktcap"
] / portfolio_characteristics.groupby("sorting_method")["mktcap"].transform("sum")

display(portfolio_characteristics.head(6))

	sorting_method	portfolio_size	portfolio_bm	n_stocks	mktcap	mktcap_share
0	Dependent	1	1	106.258278	1.883282e+08	0.006357
1	Dependent	1	2	129.678808	2.398458e+08	0.008097
2	Dependent	1	3	132.099338	2.389889e+08	0.008068
3	Dependent	1	4	154.298013	2.474387e+08	0.008353
4	Dependent	1	5	144.781457	2.081802e+08	0.007028
5	Dependent	2	1	108.443709	3.549448e+08	0.011982

def plot_heatmap(df: pd.DataFrame, value_col: str, title: str, fmt) -> None:
    methods = ["Independent", "Dependent"]
    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)
    for ax, method in zip(axes, methods):
        sub = df[df["sorting_method"] == method]
        mat = sub.pivot(
            index="portfolio_bm", columns="portfolio_size", values=value_col
        ).sort_index(ascending=True)
        # 显示时 size 1→5 从左到右，bm 1→5 从下到上更直观：这里保持矩阵行=bm
        im = ax.imshow(mat.values, aspect="auto", origin="lower", cmap="Blues")
        ax.set_xticks(range(mat.shape[1]))
        ax.set_yticks(range(mat.shape[0]))
        ax.set_xticklabels(mat.columns.tolist)
        ax.set_yticklabels(mat.index.tolist)
        ax.set_xlabel("Size portfolio")
        ax.set_ylabel("Book-to-market portfolio")
        ax.set_title(method)
        for i in range(mat.shape[0]):
            for j in range(mat.shape[1]):
                val = mat.values[i, j]
                if np.isnan(val):
                    continue
                ax.text(j, i, fmt(val), ha="center", va="center", color="black", fontsize=8)
        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    fig.suptitle(title, y=1.02)
    plt.tight_layout
    plt.show

plot_heatmap(
    portfolio_characteristics,
    "n_stocks",
    "Average number of stocks per portfolio (A-share)",
    lambda v: f"{v:.0f}",
)
plot_heatmap(
    portfolio_characteristics,
    "mktcap_share",
    "Market-cap share per portfolio (A-share)",
    lambda v: f"{100 * v:.1f}%",
)

解读要点（与原文同构，换成 A 股语境）：

• 断点取自上交所主板时，小市值组往往吸收大量深市/双创股票，股票数向小市值列倾斜。
• 独立排序对所有规模组使用同一套 BM 断点，某一规模列内 BM 组股票数可不均；依赖排序在组内重算 BM 断点，BM 维更均匀。
• 市值份额通常集中在大市值 × 低 BM（成长）一角：少数大票主导市值加权收益——这提醒我们：看收益时也要看组合构成。

6.10.6 要点小结

• 双维排序用两个特征（此处 Size 与 BM）构造投资组合，是 Fama–French 类价值/规模因子的经验基础。
• 独立排序分别赋组再交叉；依赖排序在第一维组内再分第二维。
• 会计数据必须滞后并对齐；BM 的 6 个月滞后 + 12 个月有效期是控制前视偏差的常用保守做法。
• 价值溢价对独立/依赖、断点样本、加权方式等选择敏感——与 6.9 的 non-standard errors 叙事一脉相承。
• 组合构成热力图说明：股票数≠市值份额；市值加权会被大盘成长股主导。

练习（对应原文 Exercises，改为 A 股）

1. 分别用“全部主要板块”与“仅上交所主板”算断点，比较 Size×BM 各格的平均股票数差异。
2. 仿照 6.9，画出账面权益与 BM 的截面分布（或十分位断点时序）。
3. 写函数输出 25 个独立/依赖组合上 Size、BM 的组内均值（特征画像，而不只是收益）。
4. 对照文献常见 HML 口径（如 2×3 排序），做一个小设定网格，观察价值溢价如何随设计选择变化（呼应 6.9 的 p-hacking）。

6.11 A股复现：Fama–French 因子构造（Replicating Fama-French Factors）

本节复现 tidyfinance 章节 Replicating Fama-French Factors：按 FF 官方协议（每年 6 月定组、7 月起持有一年；Size×BM 的 2×3 排序）构造 SMB / HML，并与 CSMAR 公布的月度三因子对照评价复现质量。若本地已缓存盈利与投资变量，再扩展到五因子中的 RMW / CMA。

import warnings
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display
from matplotlib import font_manager

warnings.filterwarnings("ignore", category=UserWarning)

_cn_candidates = ["Microsoft YaHei", "SimHei", "PingFang SC", "Noto Sans CJK SC", "Arial Unicode MS"]
_available = {f.name for f in font_manager.fontManager.ttflist}
for _name in _cn_candidates:
    if _name in _available:
        plt.rcParams["font.sans-serif"] = [_name]
        break
plt.rcParams["axes.unicode_minus"] = False

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"
BE_CACHE = CACHE_DIR / "book_equity_ashare.parquet"
ANN_CACHE = CACHE_DIR / "compustat_annual_ashare.parquet"  # be/at/op/inv（五因子用）
FF3_XLSX = DATA_DIR / "Fama_French_Factor_Monthly.xlsx"

# CSMAR：沪深A股+创业板+科创板；流通市值加权
CSMAR_MARKET = "P9714"

print("crsp =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE   =", BE_CACHE.exists)
print("FF3  =", FF3_XLSX.exists)
print("ANN  =", ANN_CACHE.exists, "(optional for FF5)")

crsp = True
BE   = True
FF3  = True
ANN  = True (optional for FF5)

6.11.1 数据准备

对应原文 Data Preparation。读取月度面板、年报账面权益，以及 CSMAR 三因子；板块映射用于主板断点。

def map_board(permno: str) -> str:
    code = str(permno).zfill(6)
    if code.startswith("688"):
        return "科创板"
    if code.startswith(("300", "301")):
        return "创业板"
    if code.startswith("6"):
        return "上交所主板"
    if code.startswith(("000", "001", "002", "003")):
        return "深交所主板"
    return "其他"

def load_csmar_ff3(path: Path, market: str = CSMAR_MARKET) -> pd.DataFrame:
    """Load CSMAR monthly FF3; keep float-weighted factors for one market."""
    raw = pd.read_excel(path)
    # 前两行：中文名 / 单位
    df = raw.iloc[2:].copy
    df.columns = ["MarkettypeID", "TradingMonth", "mkt_excess", "mkt_excess_total",
                  "smb", "smb_total", "hml", "hml_total"]
    df = df[df["MarkettypeID"].astype(str) == market].copy
    df["date"] = pd.to_datetime(df["TradingMonth"].astype(str) + "-01")
    for c in ["mkt_excess", "smb", "hml"]:
        df[c] = pd.to_numeric(df[c], errors="coerce")
    out = df[["date", "mkt_excess", "smb", "hml"]].dropna(subset=["smb", "hml"]).sort_values("date")
    return out.reset_index(drop=True)

crsp_monthly = pd.read_parquet(CACHE_DIR / "crsp_monthly_ashare.parquet").copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly["board"] = crsp_monthly["permno"].map(map_board)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap", "mktcap_lag"]).copy
crsp_monthly = crsp_monthly[crsp_monthly["mktcap_lag"] > 0].copy

book_equity = pd.read_parquet(BE_CACHE).copy
book_equity["permno"] = book_equity["permno"].astype(str).str.zfill(6)
book_equity["datadate"] = pd.to_datetime(book_equity["datadate"])
book_equity = book_equity[book_equity["be"] > 0].copy

factors_ff3_monthly = load_csmar_ff3(FF3_XLSX, CSMAR_MARKET)

print("crsp:", crsp_monthly.shape)
print("BE years:", book_equity["datadate"].dt.year.min, "-", book_equity["datadate"].dt.year.max)
print("CSMAR FF3:", factors_ff3_monthly["date"].min.date, "→", factors_ff3_monthly["date"].max.date,
      "| market=", CSMAR_MARKET)
display(factors_ff3_monthly.tail(3))

crsp: (884560, 11)
BE years: 1999 - 2025
CSMAR FF3: 1991-07-01 → 2026-07-01 | market= P9714

	date	mkt_excess	smb	hml
416	2026-05-01	-0.003817	-0.009226	-0.022075
417	2026-06-01	0.002287	-0.008796	-0.035981
418	2026-07-01	-0.082799	0.004447	0.204287

与 6.10 不同，FF 官方时间轴是：

1. 规模：用 (t) 年 6 月市值，从 7 月起持有到次年 6 月；
2. BM：用 (t-1) 年报账面权益 ÷ (t-1) 年 12 月市值，同样在 (t) 年 7 月投入使用（会计信息最多可滞后约 18 个月）。

下面用 sorting_date（每年 7 月）把规模与 BM 对齐到同一套年度分组键。

# 6 月市值 → 7 月 sorting_date
size = (
    crsp_monthly.loc[crsp_monthly["date"].dt.month == 6, ["permno", "board", "date", "mktcap"]]
    .assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=1))
    .rename(columns={"mktcap": "size"})
    [["permno", "board", "sorting_date", "size"]]
)

# 12 月市值 → +7 个月 = 次年 7 月（与年报 BE 对齐）
market_equity = (
    crsp_monthly.loc[crsp_monthly["date"].dt.month == 12, ["permno", "date", "mktcap"]]
    .assign(sorting_date=lambda d: d["date"] + pd.DateOffset(months=7))
    .rename(columns={"mktcap": "me"})
    [["permno", "sorting_date", "me"]]
)

book_to_market = book_equity.copy
book_to_market["sorting_date"] = pd.to_datetime(
    dict(
        year=book_to_market["datadate"].dt.year + 1,
        month=7,
        day=1,
    )
)
book_to_market = book_to_market.merge(market_equity, on=["permno", "sorting_date"], how="inner")
book_to_market["bm"] = book_to_market["be"] / book_to_market["me"]
book_to_market = book_to_market.replace([np.inf, -np.inf], np.nan)
book_to_market = book_to_market.dropna(subset=["bm"])
book_to_market = book_to_market[book_to_market["bm"] > 0]

sorting_variables = (
    size.merge(
        book_to_market[["permno", "sorting_date", "me", "bm"]],
        on=["permno", "sorting_date"],
        how="inner",
    )
    .dropna
    .drop_duplicates(["permno", "sorting_date"], keep="first")
)

print(sorting_variables.shape)
print(
    "sorting years:",
    sorting_variables["sorting_date"].dt.year.min,
    "-",
    sorting_variables["sorting_date"].dt.year.max,
)
display(sorting_variables.head(3))

(63650, 6)
sorting years: 2001 - 2026

	permno	board	sorting_date	size	me	bm
0	000001	深交所主板	2001-07-01	21328740.14	20228171.57	0.173894
1	000001	深交所主板	2002-07-01	21140429.48	17264684.07	0.210121
2	000001	深交所主板	2003-07-01	15629824.19	14784207.01	0.266529

6.11.2 组合排序（2×3）

对应原文 Portfolio Sorts。规模按主板中位数分 Small / Big；BM 按主板 30% / 70% 分位分 Low / Mid / High。交叉得 6 个组合。

def assign_portfolio(
    data: pd.DataFrame,
    sorting_variable: str,
    percentiles: list[float],
    boards: list[str] | None = None,
) -> pd.Series:
    if boards is None:
        boards = ["上交所主板"]
    subset = data[data["board"].isin(boards)]
    if len(subset) < max(2, len(percentiles) - 1):
        return pd.Series(1, index=data.index, dtype=int)

    breakpoints = subset[sorting_variable].quantile(percentiles).to_numpy(dtype=float)
    breakpoints = np.unique(breakpoints)
    if len(breakpoints) < 2:
        return pd.Series(1, index=data.index, dtype=int)

    inner = breakpoints[1:-1]
    labels = list(range(1, len(inner) + 2))
    try:
        port = pd.cut(
            data[sorting_variable],
            bins=[-np.inf, *inner, np.inf],
            labels=labels,
            include_lowest=True,
        )
    except ValueError:
        return pd.Series(1, index=data.index, dtype=int)
    return port.astype(int)

parts = []
for dt, g in sorting_variables.groupby("sorting_date", sort=True):
    g = g.copy
    g["portfolio_size"] = assign_portfolio(g, "size", [0, 0.5, 1])
    g["portfolio_bm"] = assign_portfolio(g, "bm", [0, 0.3, 0.7, 1])
    parts.append(g[["permno", "sorting_date", "portfolio_size", "portfolio_bm"]])
portfolios_annual = pd.concat(parts, ignore_index=True)

# 月度收益：1–6 月用上一年 7 月的分组；7–12 月用本年 7 月分组
crsp = crsp_monthly.copy
y = crsp["date"].dt.year
m = crsp["date"].dt.month
crsp["sorting_date"] = pd.to_datetime(
    dict(year=np.where(m <= 6, y - 1, y), month=7, day=1)
)
portfolios = crsp.merge(portfolios_annual, on=["permno", "sorting_date"], how="inner")
print("monthly portfolio-months:", portfolios.shape)
print(portfolios.groupby(["portfolio_size", "portfolio_bm"]).size.unstack(fill_value=0))

monthly portfolio-months: (689085, 14)
portfolio_bm         1       2       3
portfolio_size                        
1               122372  179872  109645
2               106402  101274   69520

6.11.3 三因子：SMB 与 HML

对应原文 Fama-French Three-Factor Model。各 2×3 格子内按 mktcap_lag 市值加权；
SMB = 三个小市值组合收益的等权平均 − 三个大市值组合等权平均；
HML = 两个高 BM 组合等权平均 − 两个低 BM 组合等权平均。

def vw_ret(g: pd.DataFrame) -> float:
    w = g["mktcap_lag"]
    return float((g["ret_excess"] * w).sum / w.sum) if w.sum > 0 else np.nan

rows = []
for keys, g in portfolios.groupby(["date", "portfolio_size", "portfolio_bm"], sort=True):
    rows.append(
        {
            "date": keys[0],
            "portfolio_size": int(keys[1]),
            "portfolio_bm": int(keys[2]),
            "ret": vw_ret(g),
        }
    )
port_rets = pd.DataFrame(rows).dropna(subset=["ret"])

factors_replicated = []
for dt, g in port_rets.groupby("date"):
    small = g.loc[g["portfolio_size"] == 1, "ret"].mean
    big = g.loc[g["portfolio_size"] == 2, "ret"].mean
    high = g.loc[g["portfolio_bm"] == g["portfolio_bm"].max, "ret"].mean
    low = g.loc[g["portfolio_bm"] == g["portfolio_bm"].min, "ret"].mean
    factors_replicated.append(
        {"date": dt, "smb_replicated": small - big, "hml_replicated": high - low}
    )
factors_replicated = pd.DataFrame(factors_replicated).sort_values("date")

print(
    factors_replicated[["smb_replicated", "hml_replicated"]]
    .agg(["mean", "std", "count"])
    .round(4)
)
display(factors_replicated.tail(3))

       smb_replicated  hml_replicated
mean          -0.0027         -0.0018
std            0.0434          0.0308
count        318.0000        318.0000

	date	smb_replicated	hml_replicated
315	2026-05-01	0.005685	0.000339
316	2026-06-01	0.023147	-0.039969
317	2026-07-01	-0.140000	0.111628

6.11.4 复现评价（对照 CSMAR）

对应原文 Replication Evaluation：把自建因子对 CSMAR 官方因子做时序回归
replicated ~ α + β · official。理想情形下 α≈0、β≈1、调整 (R^2) 接近 1。

def replication_regression(y: pd.Series, x: pd.Series, name: str) -> dict:
    df = pd.concat([y, x], axis=1).dropna
    df.columns = ["y", "x"]
    if len(df) < 24:
        return {"factor": name, "n": len(df), "alpha": np.nan, "beta": np.nan, "adj_r2": np.nan}
    model = sm.OLS(df["y"], sm.add_constant(df["x"])).fit
    return {
        "factor": name,
        "n": int(model.nobs),
        "alpha": float(model.params["const"]),
        "alpha_t": float(model.tvalues["const"]),
        "beta": float(model.params["x"]),
        "beta_t": float(model.tvalues["x"]),
        "adj_r2": float(model.rsquared_adj),
        "corr": float(df["y"].corr(df["x"])),
    }

test = factors_replicated.merge(factors_ff3_monthly, on="date", how="inner")
# 对齐到 4 位小数，减轻浮点噪音（与原文 round 一致）
test["smb_replicated"] = test["smb_replicated"].round(4)
test["hml_replicated"] = test["hml_replicated"].round(4)

rows = [
    replication_regression(test["smb_replicated"], test["smb"], "SMB"),
    replication_regression(test["hml_replicated"], test["hml"], "HML"),
]
eval_df = pd.DataFrame(rows)
display(eval_df.round(4))

fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
for ax, lab, ycol, xcol in [
    (axes[0], "SMB", "smb_replicated", "smb"),
    (axes[1], "HML", "hml_replicated", "hml"),
]:
    ax.scatter(test[xcol], test[ycol], s=12, alpha=0.55, color="#4C78A8")
    lims = [
        np.nanmin([test[xcol].min, test[ycol].min]),
        np.nanmax([test[xcol].max, test[ycol].max]),
    ]
    ax.plot(lims, lims, "r--", linewidth=1, label="45°")
    ax.set_xlabel(f"CSMAR {lab}")
    ax.set_ylabel(f"Replicated {lab}")
    ax.set_title(lab)
    ax.grid(True, alpha=0.3)
    ax.legend(frameon=False)
plt.suptitle(f"A-share FF3 replication vs CSMAR ({CSMAR_MARKET})", y=1.02)
plt.tight_layout
plt.show

print(
    f"overlap months = {len(test)}; "
    f"SMB corr = {test['smb_replicated'].corr(test['smb']):.3f}; "
    f"HML corr = {test['hml_replicated'].corr(test['hml']):.3f}"
)

	factor	n	alpha	alpha_t	beta	beta_t	adj_r2	corr
0	SMB	318	-0.0080	-9.7005	0.8133	50.0896	0.8878	0.9424
1	HML	318	-0.0036	-3.6164	0.7315	25.6705	0.6749	0.8221

overlap months = 318; SMB corr = 0.942; HML corr = 0.822

若 β 明显偏离 1 或 (R^2) 不高，常见原因包括：样本股票池与 CSMAR MarkettypeID 不完全一致、ST/金融业过滤差异、市值口径（流通 vs 总市值）、BE 定义与财报时点差异。可改用 P9706（仅沪深主板）等编码做稳健性对照。

6.11.5 五因子扩展：RMW 与 CMA

对应原文 Fama-French Five-Factor Model。在 Size 组内再对 盈利（OP）、投资（INV） 做 30/70 分位，构造：

• RMW = 高盈利 − 低盈利；
• CMA = 低投资 − 高投资；
• 五因子口径下的 SMB 常取 Size×BM / Size×OP / Size×INV 三套小减大的平均。

需要年报变量 op、inv（缓存 compustat_annual_ashare.parquet）。若缓存不存在，下面会提示跳过；你也可以用 akshare 利润表/资产负债表自行补齐。CSMAR 本文件仅含三因子，故五因子不做官方对照。

def load_annual_fundamentals(path: Path) -> pd.DataFrame | None:
    if not path.exists:
        print("未找到", path, "→ 跳过五因子构造。")
        print("可先运行 Data 下载脚本生成 be/at/op/inv 后再执行本小节。")
        return None
    ann = pd.read_parquet(path).copy
    ann["permno"] = ann["permno"].astype(str).str.zfill(6)
    ann["datadate"] = pd.to_datetime(ann["datadate"])
    need = {"be", "op", "inv"}
    if not need.issubset(ann.columns):
        print("缓存缺少列", need - set(ann.columns), "→ 跳过五因子。")
        return None
    return ann

ann = load_annual_fundamentals(ANN_CACHE)
if ann is not None:
    other = ann.copy
    other["sorting_date"] = pd.to_datetime(
        dict(year=other["datadate"].dt.year + 1, month=7, day=1)
    )
    other = other.merge(market_equity, on=["permno", "sorting_date"], how="inner")
    other["bm"] = other["be"] / other["me"]
    other = other.replace([np.inf, -np.inf], np.nan)

    sorting5 = (
        size.merge(
            other[["permno", "sorting_date", "me", "be", "bm", "op", "inv"]],
            on=["permno", "sorting_date"],
            how="inner",
        )
        .dropna(subset=["size", "bm", "op", "inv"])
        .drop_duplicates(["permno", "sorting_date"], keep="first")
    )
    # 依赖排序：先 Size，再组内 BM/OP/INV
    parts = []
    for dt, g in sorting5.groupby("sorting_date", sort=True):
        g = g.copy
        g["portfolio_size"] = assign_portfolio(g, "size", [0, 0.5, 1])
        sub = []
        for _, sg in g.groupby("portfolio_size"):
            sg = sg.copy
            sg["portfolio_bm"] = assign_portfolio(sg, "bm", [0, 0.3, 0.7, 1])
            sg["portfolio_op"] = assign_portfolio(sg, "op", [0, 0.3, 0.7, 1])
            sg["portfolio_inv"] = assign_portfolio(sg, "inv", [0, 0.3, 0.7, 1])
            sub.append(sg)
        parts.append(pd.concat(sub, ignore_index=False))
    port_ann5 = pd.concat(parts, ignore_index=True)[
        ["permno", "sorting_date", "portfolio_size", "portfolio_bm", "portfolio_op", "portfolio_inv"]
    ]
    port5 = crsp.merge(port_ann5, on=["permno", "sorting_date"], how="inner")

    def factor_from_grid(df, char_col, high_minus_low=True):
        rows = []
        for keys, g in df.groupby(["date", "portfolio_size", char_col], sort=True):
            rows.append(
                {
                    "date": keys[0],
                    "portfolio_size": int(keys[1]),
                    char_col: int(keys[2]),
                    "ret": vw_ret(g),
                }
            )
        rets = pd.DataFrame(rows).dropna(subset=["ret"])
        out = []
        for dt, g in rets.groupby("date"):
            hi = g.loc[g[char_col] == g[char_col].max, "ret"].mean
            lo = g.loc[g[char_col] == g[char_col].min, "ret"].mean
            out.append({"date": dt, "ret": (hi - lo) if high_minus_low else (lo - hi)})
        return pd.DataFrame(out)

    hml5 = factor_from_grid(port5, "portfolio_bm", True).rename(columns={"ret": "hml_replicated"})
    rmw = factor_from_grid(port5, "portfolio_op", True).rename(columns={"ret": "rmw_replicated"})
    cma = factor_from_grid(port5, "portfolio_inv", False).rename(columns={"ret": "cma_replicated"})

    # SMB：三套 size 网格的小减大平均
    smb_parts = []
    for col in ["portfolio_bm", "portfolio_op", "portfolio_inv"]:
        rows = []
        for keys, g in port5.groupby(["date", "portfolio_size", col], sort=True):
            rows.append(
                {
                    "date": keys[0],
                    "portfolio_size": int(keys[1]),
                    "ret": vw_ret(g),
                }
            )
        rets = pd.DataFrame(rows).dropna(subset=["ret"])
        for dt, g in rets.groupby("date"):
            s = g.loc[g["portfolio_size"] == 1, "ret"].mean
            b = g.loc[g["portfolio_size"] == 2, "ret"].mean
            smb_parts.append({"date": dt, "ret": s - b})
    smb5 = (
        pd.DataFrame(smb_parts).groupby("date", as_index=False)["ret"].mean
        .rename(columns={"ret": "smb_replicated"})
    )

    factors5 = (
        smb5.merge(hml5, on="date", how="outer")
        .merge(rmw, on="date", how="outer")
        .merge(cma, on="date", how="outer")
        .sort_values("date")
    )
    print("FF5 replicated — monthly means:")
    display(factors5.filter(like="_replicated").agg(["mean", "std", "count"]).T.round(4))
    display(factors5.tail(3))

FF5 replicated — monthly means:

	mean	std	count
smb_replicated	-0.0021	0.0441	318.0
hml_replicated	-0.0019	0.0311	318.0
rmw_replicated	0.0022	0.0241	318.0
cma_replicated	-0.0029	0.0163	318.0

	date	smb_replicated	hml_replicated	rmw_replicated	cma_replicated
315	2026-05-01	0.002664	-0.003790	0.032230	-0.026152
316	2026-06-01	0.015789	-0.045325	-0.012644	-0.055543
317	2026-07-01	-0.119145	0.112453	0.058547	0.047182

6.11.6 要点小结

• FF 三因子在 CAPM 市场因子之外加入 SMB（规模） 与 HML（价值）；五因子再加入 RMW（盈利） 与 CMA（投资）。
• 官方构造强调：6/7 月年度再平衡、特定会计滞后、以及（美股）NYSE /（A 股）主板断点与 2×3 排序。
• 复现质量用对官方序列的时序回归评价：关注 α、β 与调整 (R^2)。
• 6.9/6.10 展示了排序自由度；本节展示“锁定一套 FF 协议”后如何得到可对照的标准因子。

练习（对应原文 Exercises，改为 A 股）

1. 用自建市值加权市场超额收益复现 RiskPremium1，并与 CSMAR mkt_excess 做同样的回归评价。
2. 分别用 P9706（不含双创）与 P9714（含双创+科创）作对照，比较 SMB/HML 的 β 与 (R^2)。
3. 把断点从“仅上交所主板”改为全市场分位，观察复现相关系数如何变化。
4. 在缓存齐备后，补全五因子，并查找/下载 CSMAR 五因子做对照（若有权限）。

6.12 A股复现：Fama–MacBeth 回归（Fama-MacBeth Regressions）

本节复现 tidyfinance 章节 Fama-MacBeth Regressions：用个股作为检验资产，估计与 Fama–French 三因子相关特征（市场 beta、对数市值、账面市值比）的风险溢价。对应章内概念节 6.6；此处给出与原文同构的完整 A 股流水线。

Fama–MacBeth 是两步法：

1. 截面回归：每个月 (t)，用当期（已知）特征解释下月超额收益 (r_{i,t+1})；
2. 时序聚合：对每月得到的 (\hat\lambda^f_t) 取时间平均，并做 (t) 检验（常配合 Newey–West）。

线性设定示意（用特征代理风险暴露）：

$$ r_{i,t+1}=\alpha_t+\lambda^{M}_t\beta^{M}_{i,t}+\lambda^{\mathrm{Size}}_t\log(\mathrm{ME})_{i,t}+\lambda^{\mathrm{BM}}_t\mathrm{BM}_{i,t}+\epsilon_{i,t+1}. $$

import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import statsmodels.api as sm
from IPython.display import display

warnings.filterwarnings("ignore", category=UserWarning)

DATA_DIR = Path(r"D:/A_Topics/202607_03_TidyFinanceAShare/202608_02_传统量化/Data")
CACHE_DIR = DATA_DIR / "_cache_beta"

print("crsp =", (CACHE_DIR / "crsp_monthly_ashare.parquet").exists)
print("BE   =", (CACHE_DIR / "book_equity_ashare.parquet").exists)
print("beta =", (DATA_DIR / "beta_ashare.parquet").exists)

crsp = True
BE   = True
beta = True

6.12.1 数据准备

对应原文 Data Preparation。读取月度收益与市值、正的账面权益，以及月频 CAPM beta。

crsp_monthly = pd.read_parquet(
    CACHE_DIR / "crsp_monthly_ashare.parquet",
    columns=["permno", "date", "ret_excess", "mktcap"],
).copy
crsp_monthly["date"] = pd.to_datetime(crsp_monthly["date"])
crsp_monthly["permno"] = crsp_monthly["permno"].astype(str).str.zfill(6)
crsp_monthly = crsp_monthly.dropna(subset=["ret_excess", "mktcap"])
crsp_monthly = crsp_monthly[crsp_monthly["mktcap"] > 0].copy

compustat_annual = pd.read_parquet(CACHE_DIR / "book_equity_ashare.parquet").copy
compustat_annual["permno"] = compustat_annual["permno"].astype(str).str.zfill(6)
compustat_annual["datadate"] = pd.to_datetime(compustat_annual["datadate"])
compustat_annual = compustat_annual[compustat_annual["be"] > 0].copy

beta = pd.read_parquet(DATA_DIR / "beta_ashare.parquet").copy
beta["permno"] = beta["permno"].astype(str).str.zfill(6)
beta["date"] = pd.to_datetime(beta["date"])
beta = beta.loc[beta["return_type"] == "monthly", ["permno", "date", "beta"]].copy

print(crsp_monthly.shape, compustat_annual.shape, beta.shape)
display(crsp_monthly.head(2))

(890537, 4) (88343, 3) (623256, 3)

	permno	date	ret_excess	mktcap
0	000001	2000-01-01	0.060035	19843822.88
1	000001	2000-02-01	-0.013189	19618933.36

用年报 BE 与同月市值构造 BM 与 (\log(\mathrm{ME}))，并并上当月 beta；特征整体滞后 6 个月后并入收益面板，再按股票向前填充。最后把收益前移一期，得到 ret_excess_lead（用 (t) 期特征预测 (t+1) 期收益）。

# 年报日对齐到月（12 月）
compustat_annual = compustat_annual.assign(
    date=lambda d: d["datadate"].dt.to_period("M").dt.to_timestamp
)

characteristics = (
    compustat_annual.merge(
        crsp_monthly[["permno", "date", "mktcap"]],
        on=["permno", "date"],
        how="left",
    )
    .merge(beta, on=["permno", "date"], how="left")
    .assign(
        bm=lambda d: d["be"] / d["mktcap"],
        log_mktcap=lambda d: np.log(d["mktcap"]),
        sorting_date=lambda d: d["date"] + pd.DateOffset(months=6),
    )
    [["permno", "bm", "log_mktcap", "beta", "sorting_date"]]
)
characteristics = characteristics.replace([np.inf, -np.inf], np.nan)

data_fama_macbeth = (
    crsp_monthly.merge(
        characteristics,
        left_on=["permno", "date"],
        right_on=["permno", "sorting_date"],
        how="left",
    )
    .drop(columns=["sorting_date"], errors="ignore")
    .sort_values(["permno", "date"])
)
for col in ["beta", "bm", "log_mktcap"]:
    data_fama_macbeth[col] = data_fama_macbeth.groupby("permno")[col].ffill

# 构造 lead 收益：把 date 减 1 个月后与当期对齐 => 当期特征配下月收益
lead = data_fama_macbeth[["permno", "date", "ret_excess"]].copy
lead["date"] = lead["date"] - pd.DateOffset(months=1)
lead = lead.rename(columns={"ret_excess": "ret_excess_lead"})

data_fama_macbeth = (
    data_fama_macbeth.merge(lead, on=["permno", "date"], how="left")
    [["permno", "date", "ret_excess_lead", "beta", "log_mktcap", "bm"]]
    .dropna
    .copy
)

print(data_fama_macbeth.shape)
print(
    "date range:",
    data_fama_macbeth["date"].min.date,
    "→",
    data_fama_macbeth["date"].max.date,
)
display(data_fama_macbeth.head(3))

(502694, 6)
date range: 2004-06-01 → 2026-06-01

	permno	date	ret_excess_lead	beta	log_mktcap	bm
53	000001	2004-06-01	-0.058612	0.945177	16.29989	0.366434
54	000001	2004-07-01	0.004530	0.945177	16.29989	0.366434
55	000001	2004-08-01	-0.001635	0.945177	16.29989	0.366434

6.12.2 截面回归

对应原文 Cross-Sectional Regression。每个月对全市场个股估计

ret_excess_lead ~ beta + log_mktcap + bm

得到该月的风险溢价估计 (\hat\lambda_t)。

def estimate_cross_section(group: pd.DataFrame) -> pd.Series | None:
    # 样本过少或共线时跳过
    if len(group) < 50:
        return None
    y = group["ret_excess_lead"]
    X = sm.add_constant(group[["beta", "log_mktcap", "bm"]], has_constant="add")
    try:
        model = sm.OLS(y, X).fit
    except Exception:  # noqa: BLE001
        return None
    out = model.params.rename({"const": "Intercept"})
    out["date"] = group["date"].iloc[0]
    return out

rows = []
for _, g in data_fama_macbeth.groupby("date", sort=True):
    res = estimate_cross_section(g)
    if res is not None:
        rows.append(res)

risk_premiums = pd.DataFrame(rows).sort_values("date").reset_index(drop=True)
risk_premiums = risk_premiums[["date", "Intercept", "beta", "log_mktcap", "bm"]]
print(risk_premiums.shape)
display(risk_premiums.tail(3))

(265, 5)

	date	Intercept	beta	log_mktcap	bm
262	2026-04-01	-0.173627	-0.010647	0.009697	-0.006927
263	2026-05-01	-0.325958	0.029290	0.015854	-0.015826
264	2026-06-01	-0.066273	-0.111225	0.007892	0.031225

6.12.3 时序聚合

对应原文 Time-Series Aggregation。对 (\hat\lambda_t) 取时间平均作为风险溢价，并计算普通 (t) 统计量：

[ t=\frac{\overline{\hat\lambda}}{\mathrm{sd}(\hat\lambda)/\sqrt{T}}. ]

long = risk_premiums.melt(
    id_vars="date", var_name="factor", value_name="estimate"
)

price_of_risk = (
    long.groupby("factor", as_index=False)
    .agg(
        risk_premium=("estimate", "mean"),
        t_statistic=(
            "estimate",
            lambda s: s.mean / s.std(ddof=1) * np.sqrt(len(s)),
        ),
    )
    .sort_values("factor")
)
display(price_of_risk.round(4))

	factor	risk_premium	t_statistic
0	Intercept	0.0565	2.9072
1	beta	-0.0018	-1.1494
2	bm	-0.0003	-0.5250
3	log_mktcap	-0.0027	-2.4319

实务中报告风险溢价时，常对 (\hat\lambda_t) 序列做 Newey–West（HAC） 调整以处理自相关。下面用“对常数项回归 + HAC 标准误”（滞后 6 阶）计算 NW (t) 统计量，并与普通 (t) 并表——对应原文用 pyfixest / tidyfinance 的做法。

def estimate_newey_west(group: pd.DataFrame, lags: int = 6) -> pd.Series:
    g = group.sort_values("date")
    y = g["estimate"].to_numpy(dtype=float)
    X = np.ones((len(y), 1))
    fit = sm.OLS(y, X).fit(cov_type="HAC", cov_kwds={"maxlags": lags})
    return pd.Series(
        {
            "factor": g["factor"].iloc[0],
            "t_statistic_newey_west": float(y.mean / fit.bse[0]),
        }
    )

nw_rows = [estimate_newey_west(g) for _, g in long.groupby("factor")]
price_of_risk_newey_west = pd.DataFrame(nw_rows)

fm_table = (
    price_of_risk.merge(price_of_risk_newey_west, on="factor")
    .assign(
        risk_premium=lambda d: d["risk_premium"].round(3),
        t_statistic=lambda d: d["t_statistic"].round(3),
        t_statistic_newey_west=lambda d: d["t_statistic_newey_west"].round(3),
    )
)
display(fm_table)

	factor	risk_premium	t_statistic	t_statistic_newey_west
0	Intercept	0.057	2.907	3.300
1	beta	-0.002	-1.149	-1.146
2	bm	-0.000	-0.525	-0.672
3	log_mktcap	-0.003	-2.432	-2.824

6.12.4 结果解读与一键封装

解读口径与原文一致（符号方向取决于样本）：

• BM 系数为正 → 价值溢价（高 BM 预期收益更高）；
• log_mktcap 为负 → 规模溢价（小市值预期收益更高）；
• beta 往往不显著或很弱 → 与前文单维 beta 排序、SML 平坦等结果相呼应。

下面给出与 tidyfinance.estimate_fama_macbeth 同角色的封装函数，便于复用（默认报告 Newey–West (t)）。

def estimate_fama_macbeth(
    data: pd.DataFrame,
    model: str = "ret_excess_lead ~ beta + bm + log_mktcap",
    vcov_lags: int = 6,
) -> pd.DataFrame:
    """Two-step Fama-MacBeth with HAC t-stats (A-share helper)."""
    # 解析极简公式：y ~ x1 + x2 + ...
    left, right = [s.strip for s in model.split("~")]
    x_cols = [c.strip for c in right.split("+")]
    y_col = left

    prem_rows = []
    for dt, g in data.groupby("date", sort=True):
        if len(g) < 50:
            continue
        y = g[y_col]
        X = sm.add_constant(g[x_cols], has_constant="add")
        try:
            fit = sm.OLS(y, X).fit
        except Exception:  # noqa: BLE001
            continue
        params = fit.params.rename({"const": "Intercept"})
        params["date"] = dt
        prem_rows.append(params)
    prem = pd.DataFrame(prem_rows)
    if prem.empty:
        return prem

    long_ = prem.melt(id_vars="date", var_name="factor", value_name="estimate")
    out_rows = []
    for fac, g in long_.groupby("factor"):
        y = g.sort_values("date")["estimate"].to_numpy(dtype=float)
        fit = sm.OLS(y, np.ones((len(y), 1))).fit(
            cov_type="HAC", cov_kwds={"maxlags": vcov_lags}
        )
        se = float(fit.bse[0])
        mu = float(y.mean)
        out_rows.append(
            {
                "factor": fac if fac != "Intercept" else "intercept",
                "risk_premium": mu,
                "n": len(y),
                "standard_error": se,
                "t_statistic": mu / se if se > 0 else np.nan,
            }
        )
    return pd.DataFrame(out_rows).sort_values("factor").reset_index(drop=True)

display(
    estimate_fama_macbeth(
        data_fama_macbeth,
        model="ret_excess_lead ~ beta + bm + log_mktcap",
        vcov_lags=6,
    ).round(4)
)

	factor	risk_premium	n	standard_error	t_statistic
0	beta	-0.0018	265	0.0016	-1.1460
1	bm	-0.0003	265	0.0005	-0.6723
2	intercept	0.0565	265	0.0171	3.3003
3	log_mktcap	-0.0027	265	0.0010	-2.8240

6.12.5 要点小结

• Fama–MacBeth 通过对截面回归系数的时序平均估计特征风险溢价，与组合排序互补。
• 常用于检验规模、价值、市场 beta 等是否在收益截面中被定价。
• 预先估计的 beta 存在测量误差；大样本下偏误趋于消失，实务上至少应对 (\lambda_t) 使用 Newey–West 等 HAC 标准误；个股回归也可考虑按市值 WLS（Hou et al., 2020）。
• 本节用 statsmodels 手写两步法，并用 estimate_fama_macbeth 封装，对应原文的手动流程与 tidyfinance 一键接口。

练习（对应原文 Exercises，改为 A 股）

1. 用滚动窗口、相对自建或 CSMAR 的 SMB/HML 估计个股的规模/价值因子载荷，替代特征 log_mktcap/bm 再跑 FM，比较系数差异。
2. 若有行业组合收益，用行业组合代替个股估计滚动 beta/规模/价值暴露，再做 FM，看结论是否稳健。
3. 按 Hou et al. (2020) 用市值加权 WLS 做截面回归；另做“每月剔除最小 20% 市值公司后的 OLS”。比较三种设定下的风险溢价。