因子投资-实验04-实现FamaMacBeth回归

QUANT RESEARCH NOTE / LAB 04

实现一个最小 Fama–MacBeth 特征回归：第一阶段每月对股票横截面做 OLS，第二阶段对月度系数求时间均值，并用 Newey–West 标准误做推断。

类型动手实验

模块研究实验

阅读约 12 分钟

材料9 节 · 3 段代码

问题

实现一个最小 Fama–MacBeth 特征回归：第一阶段每月对股票横截面做 OLS，第二阶段对月度系数求时间均值，并用 Newey–West 标准误做推断。

最终接口：

fama_macbeth(frame, y, xs) → pd.DataFrame

返回每个系数的 mean、hac_se、t_stat 和月份数；月度系数保存在返回表的 attrs 中，便于诊断。

输入数据

使用 8 个月、每月 6 只股票。两个当期特征是 size 和 value，未来一期收益 future_return 按以下过程生成：

$$R_{i,t+1}=0.002+\lambda_{size,t}size_i+ \lambda_{value,t}value_i+\varepsilon_{i,t+1}$$

size 风险价格的 8 个月真值均值为 1%，value 为 0.5%。噪声被构造成与当月常数、size、value 三列正交，因此每月 OLS 能精确恢复生成系数；这只用于测试代数实现。

真实输入长表至少包含 return_date、asset、future_return 和所有 x；信号形成日期应另列，并断言早于 return_date。本实验省略交易实现，只检验回归。

逐步实现

1. 每月独立横截面回归

PYTHON 可复现实验

# kb-skip
for date, block in frame.groupby("return_date"):
    X = sm.add_constant(block[xs], has_constant="add")
    fit = sm.OLS(block[y], X).fit()
    monthly_rows.append({"return_date": date, **fit.params.to_dict()})

一月有 N 个股票观测，得到一组当月 λ。每个月重新估计允许风险价格随时间变化，也避免把股票数量误当作时间证据。

2. 对月度系数做均值回归

PYTHON 可复现实验

# kb-skip
fit_mean = sm.OLS(values, np.ones((len(values), 1))).fit(
    cov_type="HAC",
    cov_kwds={"maxlags": lags},
)

只对常数回归的系数就是时间均值；HAC bse 是该均值在允许有限自相关下的标准误。lags 必须由持有期和误差结构决定。

3. 保存中间结果

最终表只含摘要会隐藏不稳定月份。将 monthly 放在 attrs，实际研究还应落盘保存每月股票数、R²、条件数、缺失比例和分位数。

完整代码

PYTHON 可复现实验

import numpy as np
import pandas as pd
import statsmodels.api as sm


def fama_macbeth(frame, y, xs, lags=1):
    required = {"return_date", "asset", y, *xs}
    missing = required.difference(frame.columns)
    if missing:
        raise ValueError(f"missing columns: {sorted(missing)}")
    if frame.duplicated(["return_date", "asset"]).any():
        raise ValueError("return_date and asset must be unique")
    if frame[list(required)].isna().any().any():
        raise ValueError("input contains missing values")

    names = ["const", *xs]
    monthly_rows = []
    for date, block in frame.groupby("return_date", sort=True):
        if len(block) <= len(names):
            raise ValueError("each month needs more assets than coefficients")
        X = sm.add_constant(block[xs], has_constant="add")
        fit = sm.OLS(block[y], X).fit()
        monthly_rows.append({
            "return_date": date,
            **{name: float(fit.params[name]) for name in names},
        })

    monthly = pd.DataFrame(monthly_rows).sort_values("return_date")
    summary_rows = []
    for name in names:
        values = monthly[name].to_numpy()
        mean_fit = sm.OLS(values, np.ones((len(values), 1))).fit(
            cov_type="HAC",
            cov_kwds={"maxlags": lags},
        )
        mean = float(mean_fit.params[0])
        se = float(mean_fit.bse[0])
        summary_rows.append({
            "variable": name,
            "mean": mean,
            "hac_se": se,
            "t_stat": mean / se if se > 0 else np.nan,
            "n_months": len(values),
        })

    summary = pd.DataFrame(summary_rows).set_index("variable")
    summary.attrs["monthly"] = monthly
    return summary


assets = list("ABCDEF")
dates = pd.date_range("2026-01-31", periods=8, freq="ME")
size = np.array([-1.0, -0.5, 0.0, 0.2, 0.7, 1.0])
value = np.array([0.3, -0.8, 1.0, -0.4, 0.5, -0.2])
lambda_size = np.array([0.010, 0.020, -0.005, 0.015, 0.000, 0.025, 0.010, 0.005])
lambda_value = np.array([0.005, -0.010, 0.015, 0.000, 0.010, 0.005, -0.005, 0.020])

design = np.column_stack([np.ones(6), size, value])
raw_noise = np.array([0.0010, -0.0010, 0.0005, -0.0005, 0.0008, -0.0008])
projection = design @ np.linalg.lstsq(design, raw_noise, rcond=None)[0]
orthogonal_noise = raw_noise - projection

rows = []
for month, date in enumerate(dates):
    future = (
        0.002
        + lambda_size[month] * size
        + lambda_value[month] * value
        + (1 + 0.1 * month) * orthogonal_noise
    )
    for asset, s, v, ret in zip(assets, size, value, future):
        rows.append({
            "return_date": date,
            "asset": asset,
            "size": s,
            "value": v,
            "future_return": ret,
        })

sample = pd.DataFrame(rows)
summary = fama_macbeth(sample, "future_return", ["size", "value"])
monthly = summary.attrs["monthly"]

first = sample.loc[sample["return_date"] == dates[0]]
direct = sm.OLS(
    first["future_return"],
    sm.add_constant(first[["size", "value"]], has_constant="add"),
).fit()

assert np.allclose(
    monthly.loc[0, ["const", "size", "value"]].to_numpy(dtype=float),
    direct.params.to_numpy(),
)
assert np.isclose(summary.loc["size", "mean"], lambda_size.mean())
assert np.isclose(summary.loc["value", "mean"], lambda_value.mean())
assert (summary["n_months"] == 8).all()

print(summary.round(6).to_string())
print(monthly.round(6).to_string(index=False))

预期输出

核心均值：

预期输出 3 行 × 3 列

variable	mean	n_months
const	0.002	8
size	0.010	8
value	0.005	8

每月 size、value 系数应分别与代码中的 lambda_size、lambda_value 相同。hac_se 和 t_stat 由这 8 个时间值及 maxlags=1 决定。8 个月样本过短，数值只验证算法，不应解释为经济证据。

断言

1. return_date、asset 联合唯一。
2. 必需字段无缺失。
3. 每月股票数大于回归系数数，避免零残差自由度。
4. 函数第一月系数与直接 statsmodels OLS 完全一致。
5. 第二阶段均值等于生成时 λ 的算术平均。
6. 每个摘要系数使用恰好 8 个月。

真实点时研究应再加入 signal_date<return_date、每月秩、条件数、极端杠杆点、行业覆盖和股票池变化断言。

结果解释

第一阶段回答：“在这个月的股票横截面中，size 或 value 增加一个单位，对应多少未来收益差？”第二阶段回答：“这种月度关系的时间平均是否稳定偏离零？”

Fama–MacBeth 的有效时间证据主要来自 8 个月，而不是 48 个股票月份。股票数增加能提高每月截面系数精度，却不能把共同宏观冲击复制成独立月份。

本例直接使用可观察特征。如果 x 是第一步时序回归估出的 Beta，第二步还会继承 Beta 的测量误差，需要 Shanken 型修正、联合估计或与设计匹配的自助法。

HAC 只改变第二阶段均值的标准误，不改变月度系数或均值。若通过调整 maxlags 来最大化 t 值，推断已经被研究者选择污染。

常见故障

• 把 return_date 当信号日，导致 future_return 与特征同期或反向错位。
• 全部股票月份池化 OLS 后用普通标准误。
• 某月股票数过少、特征常数或高度共线，仍保留无意义系数。
• 每月标准化使用全样本均值和标准差。
• 缺失值按未来是否有收益决定删除。
• Newey–West 滞后未经说明，或根据结果反复改变。
• 只保存摘要，无法检查系数是否由单个极端月份驱动。
• 把特征斜率直接叫风险溢价，没有证明特征代表可定价暴露。

练习

1. 把 lambda_size 后四个月全部设为负，观察全样本均值怎样掩盖结构变化。
2. 将持有期改为 3 个月并构造重叠收益，比较 maxlags=0、2、6 的标准误。
3. 某个月把 value 设为 size 的两倍，添加条件数检查让函数失败。
4. 在每月先对 size、value 做当期 z-score，解释系数单位变化。
5. 加入行业哑变量，说明 const、特征系数和自由度如何变化。
6. 用两维聚类面板回归做对照，解释它和 Fama–MacBeth 各允许什么相关结构。