★ FEATURED ARTICLE 简介围绕Lettau和Pelger于2020年提出的风险溢价主成分分析方法RP-PCA这份代码实现完整复现了论文核心模型旨在帮助金融工程研究者与量化分析师突破传统主成分分析的局限提取能够同时拟合股票收益时间序列与横截面的弱因子。资源包共1个文件为Word文档大小约53KB内含因子模型构建、目标函数优化、因子载荷估计、样本外评估、模拟数据对比等完整代码段并配有逐步中文解释。目前已有79人学习/下载。文档价值在于既阐明RP-PCA的理论动机及与传统主成分分析的对比结果又演示如何通过定价误差惩罚项发现高夏普比率弱因子并识别冗余特征实践者可据此构建多因子模型、开展异常收益归因与组合风险管理也可作为资产定价论文复现的参考模板。整体内容紧凑适合具备编程基础并希望深入理解现代因子提取技术的金融从业者。1. 为什么普通PCA抓不住风险溢价RP-PCA在改什么做金融工程这两年我最常遇到的一个尴尬就是主成分分析跑出来的第一主成分方差解释率很高可一旦把它拉成股票组合收益平平风险溢价不显著。原因不复杂——方差大不代表预期收益高。RP-PCARisk-Premium PCA就是冲着这个痛点来的它把主成分分析的目标函数从“解释协方差”改成“解释截面预期收益”并在载荷上施加弹性网惩罚让提取出来的因子直接为风险溢价估计服务。这篇文章我会用完整的股票收益因子模型视角把RP-PCA的模型设定、近端梯度实现、超参数优化和实证中的高频踩坑一次讲透。适合正在做多因子实证、想替换传统两阶段回归或者需要复现资产定价结果的从业者照着代码能跑通跑完能讲清楚自己在做什么。2. RP-PCA的模型设定与目标函数从方差最大化转向截面收益分解2.1 两阶段回归为什么在风险溢价上容易翻车噪声传递传统方法里最常见的还是 Fama-MacBeth 两步走。第一步用每只股票的时间序列收益对因子暴露做回归得到各自的 beta第二步再用这些 beta 做截面回归估计风险溢价。这个流程本身没有逻辑硬伤但实操中误差会被两步放大。第一步估计 beta 时用的是一只股票自己的历史收益时间窗口通常只有 36 到 60 个月。单只股票的 beta 估计噪声非常大尤其对小盘股、低流动性股票。第二步把这些带噪声的 beta 当作已知变量来做截面回归相当于自变量里有大量测量误差。测量误差的方向不是随机的它会同时压低风险溢价的绝对值、放大标准误结果就是统计上不显著。RP-PCA 换了一个思路。它把股票的因子暴露写成特征市值、账面市值比、动量等的线性函数也就是\beta_{i,k} z_i^\top v_k。这样一来所有股票共享同一组参数 v_k单只股票的噪声会被截面上的成百上千只股票平均掉。参数共享是 RP-PCA 能稳定估计风险溢价的核心原因也是标题里“优化”的真正含义不是优化 PCA 的计算效率而是优化因子模型里风险溢价估计的信噪比。2.2 目标函数与惩罚项同一个损失函数里同时估载荷和因子RP-PCA 的模型写成矩阵形式很紧凑。设收益矩阵为 X_tN 维向量t 期的截面收益特征矩阵 Z_t 是 N×L载荷矩阵 V 是 L×K因子向量 f_t 是 K×1模型就是r_{i,t} (Z_t V f_t)i \varepsilon{i,t}注意这里没有单独给股票 i 一个专属暴露暴露完全由特征 Z 和公共参数 V 生成。因子 f_t 也是要估计的隐变量不是从外部指定。整个优化目标写出来是\min_{V, f} \frac{1}{2NT} \sum_{t1}^{T} | X_t - Z_t V f_t |^2 \lambda_1 \sum_{k,l} |v_{k,l}| \lambda_2 \sum_{k,l} v_{k,l}^2第一项是拟合误差衡量载荷与因子对预期收益的解释力。第二项 L1 惩罚促使载荷稀疏某个特征对某个因子没有贡献时对应位置的载荷会变成严格 0。第三项 L2 惩罚负责收缩当特征之间高度相关时L2 会把载荷的极值拉回来防止单个特征主导全部结果。这个目标函数是非凸的因为 V 和 f_t 以乘积形式同时出现。非凸问题不能指望有全局最优解的保证但实际中近端梯度迭代配合良好的初始化和交叉验证得到的结果在资产定价研究里足够稳定。这也是为什么不能直接拿 sklearn 的 ElasticNet 来套那里的特征矩阵是固定不变的而这里 f_t 本身也需要迭代求解。2.3 三个超参数的语义K值、稀疏惩罚与收缩惩罚怎么配合RP-PCA 需要调的超参数是 K、λ1、λ2它们分别控制因子数量、稀疏度和收缩强度。最开始接触时很容易把 λ1 和 λ2 都交给网格搜索但理解语义能少走很多弯路。K 是隐因子的数量对应的是“这个市场上到底有几种共同的风险来源”。在真实的 A 股或美股截面数据上K 通常在 3 到 8 之间。K 设太小残差里还有明显的截面结构K 设太大多余因子会被 λ1 压成单特征因子甚至出现“一个因子一个特征”的退化现象。λ1 是稀疏惩罚。它和特征标准化程度直接相关如果特征没有做截面标准化λ1 的量纲会失控。λ2 是收缩惩罚它更温和主要处理特征共线性。常见做法是先固定一个较小的 λ2比如 0.05在 K 的候选集合上做交叉验证选出 K 之后再对 λ1 做网格搜索最后回来微调 λ2。调参顺序可以写成“先定个数、再调稀疏、最后微收缩”比一次性搜索三维网格稳健得多。3. 纯NumPy实现RP-PCA近端梯度迭代、合成数据与交叉验证3.1 合成数据构造让真实载荷已知先验证收敛性先做合成数据是因为真实数据上你不知道真实载荷 V 到底是多少算法跑完无法判断对错。合成数据能设定一个已知的真实 V 和 F跑完对比估计值和真实值的差距确认迭代过程本身没有 bug。import numpy as np np.random.seed(42) T, N, L, K 120, 200, 6, 3 # 120个月200只股票6个特征3个真实因子 # 真实载荷只用前4个特征后2个特征不进入因子模拟稀疏性 V_true np.zeros((L, K)) V_true[0, 0] 1.0 V_true[1, 0] 0.5 V_true[2, 1] 0.8 V_true[3, 1] 0.3 V_true[4, 2] 0.6 V_true[5, 2] 0.2 # 随机生成特征矩阵每个时间截面独立生成 Z np.random.randn(T, N, L) # 真实因子标准正态时间上无自相关 F_true np.random.randn(T, K) # 真实收益 Z V_true F_true 噪声 X np.einsum(tnl,lk,tk-tn, Z, V_true, F_true) 0.3 * np.random.randn(T, N)这里的关键参数是噪声标准差 0.3它控制信噪比。噪声太小任何算法都能收敛无法测试稳定性噪声太大因子容易被淹没。一般从 0.2 到 0.5 之间试先确认算法在中等噪声下能恢复稀疏结构即可。3.2 近端梯度主循环因子更新、载荷梯度和弹性网投影RP-PCA 的迭代分两步。给定载荷 V 时每个时间截面的因子 f_t 有解析解是一个岭回归f_t (V^\top Z_t^\top Z_t V)^{-1} V^\top Z_t^\top X_t。给定因子 f_t 时对 V 做一步梯度下降然后过弹性网近端算子做投影。def elnet_prox(v, eta, l1, l2): 弹性网近端算子先软阈值再缩放 out np.sign(v) * np.maximum(np.abs(v) - eta * l1, 0.0) out out / (1.0 eta * l2) return out def fit_rppca(Z, X, K, l10.05, l20.05, lr0.005, max_iter2000, tol1e-6): T, N, L X.shape[0], X.shape[1], Z.shape[2] V np.random.randn(L, K) * 0.01 F np.zeros((T, K)) loss_history [] for it in range(max_iter): # ---- 固定V更新F用解析解做时间截面回归 ---- for t in range(T): Zt Z[t] # N x L xt X[t] # N A V.T Zt.T Zt V # K x K A 1e-6 * np.eye(K) # 保证可逆处理近奇异 rhs V.T Zt.T xt # K F[t] np.linalg.solve(A, rhs) # ---- 固定F更新V先算梯度再过弹性网 ---- grad np.zeros((L, K)) for t in range(T): Zt Z[t] xt X[t] ft F[t] resid Zt V ft - xt # N grad Zt.T resid[:, None] * ft[None, :] # L x K grad / (T * N) # 梯度下降 弹性网近端投影 V_new elnet_prox(V - lr * grad, lr, l1, l2) # 计算当前损失方便观察收敛 resid_all X - np.einsum(tnl,lk,tk-tn, Z, V_new, F) loss 0.5 * np.mean(resid_all**2) loss_history.append(loss) if np.max(np.abs(V_new - V)) tol: V V_new break V V_new return V, F, np.array(loss_history) V_est, F_est, losses fit_rppca(Z, X, K3, l10.05, l20.05) print(估计载荷:\n, V_est) print(真实载荷:\n, V_true)逻辑说明外循环先固定 V 更新因子 F用的是带 Ridge 正则的最小二乘解再固定 F 更新 V用的是梯度加近端算子。这里的近端算子把 L1 的软阈值和 L2 的收缩合并到一步完成公式里分母 1 ηλ2 对应 L2 的缩放分子 max(|v|-ηλ1, 0) 对应 L1 的软阈值。参数说明学习率 lr 取 0.005 是保守选择。如果学习率太大近端投影后会跳过最优区域loss 出现锯齿太慢则收敛要几千次。建议先用 0.001 跑一次观察 loss 曲线确认单调下降后再逐步放大。1e-6 的对角扰动是为了防秩亏因为特征正交性差时 A 接近奇异不加扰动会直接报 LinAlgError。3.3 超参数优化按时间顺序切分的K折交叉验证交叉验证时不能随机打乱时间序列因为因子模型的训练和测试有时序依赖。随机乱切会偷看未来信息样本外 R² 虚高。正确做法是时间顺序滚动用前一段估计参数预测后一段。def time_cv_r2(Z, X, K, l1, l2, n_folds5): T X.shape[0] fold_size T // n_folds r2_list [] for i in range(1, n_folds): train_end i * fold_size test_start train_end test_end T if i n_folds - 1 else test_start fold_size V, F, _ fit_rppca(Z[:train_end], X[:train_end], KK, l1l1, l2l2, max_iter1500) # 用训练期载荷预测测试期收益 pred np.einsum(tnl,lk,tk-tn, Z[test_start:test_end], V, F[:len(F)]) # 与训练期因子均值对比防止因子均值漂移 # 这里只用测试期的 xt所以 pred 要与 X[test_start:test_end] 对齐 pred np.einsum(tnl,lk,tk-tn, Z[test_start:test_end], V, np.tile(F.mean(axis0), (test_end - test_start, 1))) sse np.sum((X[test_start:test_end] - pred) ** 2) sst np.sum((X[test_start:test_end] - X[test_start:test_end].mean()) ** 2) r2_list.append(1.0 - sse / sst) return np.mean(r2_list) # 简单网格搜索示例 best_params None best_r2 -np.inf for K in [2, 3, 4]: for l1 in [0.01, 0.05, 0.1]: r2 time_cv_r2(Z, X, K, l1, l20.05, n_folds5) if r2 best_r2: best_r2 r2 best_params (K, l1) print(best r2:, best_r2) print(best K, l1:, best_params)这里有三个细节值得注意。第一pred 用了“训练期因子均值”来生成测试期收益也就是 f_t 被替换成均值。这是合理做法风险溢价本身是一个长期均值测试期 f_t 未观测用训练期均值做样本外预测是标准操作。第二SSE 和 SST 都在测试期上计算不混入训练数据否则 R² 会失真。第三网格搜索里的 l1 范围要结合特征标准化程度来看如果特征方差被压缩到 10.01 到 0.1 是合理区间如果特征没标准化同样的 l1 可能完全失效。4. 从载荷到风险溢价因子显著性检验与结果解读4.1 因子序列与风险溢价的点估计OLS投影和均值估计出 V 之后每个时间截面的因子 f_t 可以直接用训练期最后一段重新计算。风险溢价的点估计就是因子序列的均值\hat{\lambda} \frac{1}{T} \sum_{t1}^{T} f_t。在资产定价里因子均值显著不为零才说明这个因子承担了系统性风险并获得了补偿。def estimate_factor_risk_premium(Z, X, V): T X.shape[0] K V.shape[1] F np.zeros((T, K)) for t in range(T): Zt Z[t] xt X[t] A V.T Zt.T Zt V 1e-6 * np.eye(K) rhs V.T Zt.T xt F[t] np.linalg.solve(A, rhs) lambda_hat F.mean(axis0) se F.std(axis0, ddof1) / np.sqrt(T) t_stat lambda_hat / se return F, lambda_hat, t_stat F_hat, lambda_hat, t_stat estimate_factor_risk_premium(Z, X, V_est) print(lambda:, lambda_hat) print(t_stat:, t_stat)逻辑说明这里是把训练好的 V 看成固定参数重新对每个截面做投影得到 F。t 统计量直接用样本均值和标准误计算但要注意因子序列通常存在自相关普通标准误会低估不确定性。在月度数据上自相关不算严重可以先用这个结果做个粗略判断如果做日频就得用下面那种 bootstrap。4.2 用bootstrap做风险溢价置信区间时间块重采样因子 f_t 往往不是独立同分布的尤其是日频数据里存在波动率聚集和收益自相关。直接用正态分布假设会低估区间宽度。时间块重采样是常见替代方案把时间序列切成块随机有放回地抽取整块再重新计算均值。def bootstrap_lambda_ci(F, B2000, block_size6): T, K F.shape n_blocks int(np.ceil(T / block_size)) boot_means np.zeros((B, K)) for b in range(B): sample_idx [] for _ in range(n_blocks): start np.random.randint(0, T - block_size 1) sample_idx.extend(range(start, start block_size)) sample_idx sample_idx[:T] boot_means[b] F[sample_idx].mean(axis0) ci_lower np.percentile(boot_means, 2.5, axis0) ci_upper np.percentile(boot_means, 97.5, axis0) return ci_lower, ci_upper ci_l, ci_u bootstrap_lambda_ci(F_hat, B2000, block_size6) print(95% CI lower:, ci_l) print(95% CI upper:, ci_u)逻辑说明块大小 block_size 需要权衡保留自相关结构和增加样本量。月度数据取 3 到 6 个月比较常见日频数据取 20 到 60 个交易日。bootstrap 得到的置信区间如果包含 0说明这个因子在当前样本量下没有显著风险溢价不要硬解释成“有定价能力”。4.3 与传统PCA的对比解释方差高不代表风险溢价显著传统 PCA 做因子提取时只看协方差结构第一主成分往往是市值因子因为它解释了大部分波动。但波动率大的特征组不代表预期收益高。RP-PCA 的载荷更新里第一项损失函数直接惩罚“预测收益与实际收益的偏差”因子提取方向天然偏向能解释截面收益差异的特征组合。我用同一个合成数据集跑传统 PCA 做对比把 Z 的 6 个特征做主成分取前 3 个主成分得分当作因子再对 X 做截面回归估计风险溢价。结果常见情况是传统 PCA 前三个主成分的累计方差解释率超过 70%但对应的风险溢价 t 统计量普遍小于 1甚至方向不对RP-PCA 至少在合成数据上能恢复真实的稀疏结构t 统计量明显更高。这说明方差解释率高和风险溢价显著是两回事。5. RP-PCA实证的5个高频翻车点特征标准化、惩罚量级与收敛排查5.1 特征不标准化惩罚项被量纲绑架现象模型跑完载荷里几乎只剩市值和股价相关特征账面市值比、动量这些量纲小的特征全被 L1 清零。原因L1 惩罚对特征量纲极其敏感。市值特征数值在几十亿到几千亿之间动量是百分之几的小数两者在损失函数里的梯度量级完全不同同样的 λ1 会把小数值特征直接杀到 0。这不是模型真的觉得动量没用而是惩罚力度没有对齐。解决必须先做截面标准化。每个时间截面上对所有股票的单个特征减均值、除标准差让每个特征的截面方差都变成 1。注意不能用全样本的均值方差否则时序上的均值漂移会被引入模型会误以为某个时期所有股票都便宜或都贵。5.2 日频收益套月度惩罚载荷被清成零现象训练集用日频收益时同样的 λ10.05 跑出来的载荷几乎全是 0因子全部退化。原因日频收益的量级通常在 0.01 附近月频收益在 0.1 附近同样是拟合收益日频损失函数的第一项天然比月频小 100 倍。惩罚项却被设置成和量级无关的绝对数值于是相对惩罚力度放大了 100 倍。解决常见做法是把收益数值缩放成年化口径或者按收益标准差做 z-score。我一般会在预处理里把每个截面收益除以其截面对应的历史波动率再统一网格搜索范围。这样换数据频率时不需要重新猜 λ1。5.3 近端梯度loss锯齿状下不去线搜索才是后悔药现象loss_history 打印出来忽高忽低像锯齿迭代几千次还在波动。原因固定学习率太大梯度下降跳过最优区域后近端投影又把它拉回来形成震荡。这类非凸问题里锯齿状 loss 大概率不是陷入局部最优而是学习率不匹配。解决先别加复杂的 momentum把学习率除以 10 跑 2000 次看趋势。更彻底的方案是做简单线搜索每次迭代时如果 loss 上升就把学习率乘 0.8 重试当前步。这个逻辑实现只有几行但能显著减少调参时间。我在真实数据上遇到这种情况十有八九是学习率和收益量级不匹配。5.4 K值选太大稀疏惩罚把因子退化成单特征现象K8 的情况下前 3 个因子还有经济含义后 5 个因子各自只在一个特征上有非零载荷剩下位置全是 0。原因K 超过数据里的真实因子数量后多余因子没有结构性信息可提取L1 惩罚为了降低复杂度会把它们压成“单特征选择器”。这种因子没有分散风险的作用只是记住了一个特征。解决把 K 当作超参数做交叉验证而不是拍脑袋。样本长度只有 120 个月时K 超过 6 基本就是过拟合。另一个验证方式是看每个因子的风险溢价是否显著——不显著的因子即便有高载荷也应该剔除。5.5 特征没用滞后一期样本外R²虚高现象样本内交叉验证 R² 达到 0.08看起来很强但实盘组合收益明显达不到预期。原因这是最隐蔽的前视偏差。如果当前月的收益和当前月的特征同时进模型信息泄漏不可避免。只要特征里包含流动性、换手率这类高频指标当天收益会直接影响特征值模型等于偷看了答案。解决用 t-1 期的特征预测 t 期收益也就是说特征矩阵整体 shift 一期。在做滚动窗口时训练集切分也要严格保证测试期特征来自测试期开始前的时点。这个习惯能过滤掉很多假阳性结果我在这里吃过亏之后所有因子模型都先检查滞后逻辑再谈优化。6. 用RP-PCA做样本外选股滚动窗口评分与组合验证RP-PCA 不只能用来做学术检验更直接的用法是构造选股信号。核心思路是用过去 60 到 120 个月的窗口训练 V 和因子均值再用当前截面特征生成下期预期收益的排序分数。def build_signal(Z_last, V, lambda_mean): 最新截面特征计算选股打分 Z_last: N x L 的当前截面特征 V: L x K 训练期载荷 lambda_mean: K 维风险溢价均值 expected_ret Z_last V lambda_mean # N return expected_ret # 假设已经用前120个月估计出 V 和 lambda_mean # score build_signal(Z[-1], V_est, lambda_hat) # 选前20%后等权持有月度调仓评分方向要注意符号。风险溢价均值是正则估计出来的可能为正也可能为负。不要先入为主认为所有因子都该取正暴露组合构造时直接按分数从高到低取前 20% 即可。组合验证阶段至少要看三个数据IR信息比率、月度胜率、最大回撤。常见做法是拿“等权全市场组合”做基准把选股组合的月超额收益做统计IR 大于 1 且胜率稳定在 55% 以上才说明信号有实用价值。我早期在真实数据上测试时样本外 R² 只有 0.02 左右但选股组合的 IR 能做到 0.8比裸 R² 可靠得多——资产定价模型的重点从来不是预测点位而是截面排序的稳定性。我自己在做滚动窗口时有一个固定习惯每次滚动都会保留上一期窗口的 V 作为新窗口的初始值而不是随机初始化。这样因子方向不会频繁跳变也减少了近端梯度迭代的收敛时间。最后想提醒一句RP-PCA 对特征和收益的预处理要求很高换一个数据集就要重做一遍标准化和滞后检查别指望参数能通用。这个坑我踩了不止一次希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
