第 09 章 主成分分析与因子模型

本章对应 Tsay 原书第 9 章。第 08a 章的 VAR 在维数稍高时参数就会爆炸：\(k=30\) 只股票的协方差矩阵已有 465 个参数，而月度数据往往只有几十到几百期。本章的思路是降维：假设众多资产的收益由少数共同因子驱动，用因子结构压缩协方差矩阵。这正是商业多因子风险模型（Barra、Axioma）、Fama–French 因子研究和 PCA 统计套利的共同理论基础。第 05 册第 17d 章从预测角度讲动态因子模型与主成分回归，第 05 册第 19c 章给出主成分的矩阵推导，可对照阅读。

学习目标

1. 掌握一般因子模型 \(\boldsymbol r_t=\boldsymbol\alpha+\boldsymbol\beta\boldsymbol f_t+\boldsymbol\epsilon_t\) 及其协方差结构 \(\boldsymbol\Sigma_r=\boldsymbol\beta\boldsymbol\Sigma_f\boldsymbol\beta'+\boldsymbol D\)，能说出它把参数从 \(k(k+1)/2\) 压缩到多少。
2. 分清三类因子模型的"已知 / 未知"：宏观因子模型（\(\boldsymbol f_t\) 已知，时间序列回归求 \(\boldsymbol\beta\)）、基本面因子模型（BARRA：\(\boldsymbol\beta\) 已知，截面 WLS 求 \(\boldsymbol f_t\)；Fama–French：排序构造组合得 \(\boldsymbol f_t\)，再求 \(\boldsymbol\beta\)）、统计因子模型（两者都未知）。
3. 会实现 BARRA 两步 GLS 估计，理解 GLS 权重就是"因子模拟组合"。
4. 掌握 PCA 的理论（特征分解、方差分解、碎石图）和经验用法，知道股票第一主成分≈市场、债券前两个主成分≈水平与斜率。
5. 掌握正交因子模型、主成分法与极大似然估计、LR 检验定因子数、varimax 旋转；了解 \(k>T\) 时的渐近主成分分析和 Bai–Ng 准则。
6. 会用全局最小方差组合（GMVP）比较不同协方差估计，理解因子模型在组合优化中的价值。

读前导读

这一章在解决什么问题。 做组合优化需要协方差矩阵，但资产一多，样本协方差就既估不准、又可能不可逆：30 只股票有 465 个协方差参数，而你手里可能只有 60 个月的数据。本章的办法是假设"所有股票的共同运动只来自少数几个因子"。你在 CFA 里学过的单指数模型就是最简单的例子：\(r_i=\alpha_i+\beta_ir_m+\epsilon_i\)，于是任意两只股票的协方差只能通过市场传递，\(\mathrm{Cov}(r_i,r_j)=\beta_i\beta_j\sigma_m^2\)。本章把它推广到多个因子，并回答三种情形下怎么估计：因子已知（宏观因子、市场收益）、暴露已知（Barra 的行业与风格）、两者都不知道（统计因子：PCA 和因子分析）。

你熟悉的 APT、Fama–French 三因子、Barra 风险报告里的"因子风险 + 特质风险"分解，在本章都会落到同一个公式 \(\boldsymbol\Sigma_r=\boldsymbol\beta\boldsymbol\Sigma_f\boldsymbol\beta'+\boldsymbol D\) 上。PCA 部分则回答另一个问题："如果不预设任何因子，数据自己会告诉我们哪几个方向最重要？"你在固定收益里听过的"水平、斜率、曲率"三因子，正是国债收益率 PCA 的前三个主成分。

需要先想起来的数学。

• 组合方差的矩阵写法。 组合 \(\boldsymbol w\) 的方差是 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w\)，两个组合的协方差是 \(\boldsymbol w_1'\boldsymbol\Sigma\boldsymbol w_2\)；更一般地，\(\mathrm{Cov}(\boldsymbol A\boldsymbol x)=\boldsymbol A\,\mathrm{Cov}(\boldsymbol x)\boldsymbol A'\)。例：两资产各 50%，方差 0.04、0.09、协方差 0.03，组合方差 \(=0.25\times0.04+0.25\times0.09+2\times0.25\times0.03=0.0475\)。见 第 00 册第 06 章 线性代数速成。
• 对称矩阵的特征分解。 协方差矩阵对称，可写成 \(\boldsymbol\Sigma=\sum_i\lambda_i\boldsymbol e_i\boldsymbol e_i'\)，特征向量两两正交、长度为 1，特征值非负。直观上：把风险拆成 \(k\) 个互不相关的"方向"，第 \(i\) 个方向上的方差是 \(\lambda_i\)。例：\(\begin{bmatrix}2&1\\1&2\end{bmatrix}\) 的特征值 3、1，对应方向 \((1,1)/\sqrt2\)（同涨同跌）和 \((1,-1)/\sqrt2\)（一涨一跌）。见第 00 册第 06 章。
• 迹与行列式。 迹 \(\mathrm{tr}(\boldsymbol\Sigma)\) 是对角元之和，也等于全部特征值之和，即"总方差"；行列式等于特征值之积。见第 00 册第 06 章。
• 拉格朗日乘子。 求 \(f(\boldsymbol w)\) 在约束 \(g(\boldsymbol w)=c\) 下的极值：构造 \(f-\lambda(g-c)\)，对 \(\boldsymbol w\) 求偏导令其为零。你在 CFA 里见过的有效前沿、最小方差组合就是这样解出来的。本章 GMVP、因子模拟组合、PCA 的推导都用它。见 第 00 册第 05 章 多元微积分与优化。
• 加权最小二乘（WLS/GLS）。 误差方差不相等时，给噪声大的观测较小的权重：\(\hat{\boldsymbol b}=(\boldsymbol X'\boldsymbol W\boldsymbol X)^{-1}\boldsymbol X'\boldsymbol W\boldsymbol y\)，\(\boldsymbol W\) 取误差方差的倒数。这是 CFA 二级"异方差修正"的标准做法。

怎么读这一章。 核心必读：9.1（因子模型与协方差结构）、9.2.2（单因子市场模型与 GMVP）、9.3.1（Barra 截面回归与因子模拟组合）、9.4（PCA）。9.3.2 Fama–French 你大概已经熟悉，快速浏览即可。9.5 统计因子分析中，正交因子模型的假设和主成分法要读；ML 估计的 LR 检验、varimax 旋转第一次可以只看结论和例子。9.6 渐近主成分适合需要处理"股票数多于期数"的读者，第一次可跳过。建议顺序：9.1 → 9.2 → 9.3 → 9.4 → 示例二、三 → 9.5、9.6。

9.1 一般因子模型

9.1.1 模型与假设

设有 \(k\) 个资产、\(T\) 期，

\(f_{jt}\) 是 \(m\) 个共同因子（common factors），\(\beta_{ij}\) 是资产 \(i\) 对因子 \(j\) 的因子载荷（factor loading），\(\epsilon_{it}\) 是特质因子（specific factor）。假设：

• \(\boldsymbol f_t\) 是 \(m\) 维平稳过程，\(E(\boldsymbol f_t)=\boldsymbol\mu_f\)，\(\mathrm{Cov}(\boldsymbol f_t)=\boldsymbol\Sigma_f\)；
• \(E(\epsilon_{it})=0\)，\(\epsilon\) 与所有因子不相关；
• 特质因子之间、不同时点之间都不相关：\(\mathrm{Cov}(\epsilon_{it},\epsilon_{js})=\sigma_i^2\) 当 \(i=j,t=s\)，否则为 0。

注意共同因子之间不必不相关。因子分析通常还假设因子无序列相关；若收益有序列相关，先用第 08a 章的模型去除。

写成截面形式：

于是

这一行就是所有因子风险模型的核心。参数个数从 \(k(k+1)/2\) 降为 \(km+m(m+1)/2+k\)：\(k=500\)、\(m=40\) 时，从约 12.5 万降为约 2.1 万，而且估计更稳定。

推导拆解：\(\boldsymbol r_t\) 由两部分相加：\(\boldsymbol\beta\boldsymbol f_t\) 和 \(\boldsymbol\epsilon_t\)（\(\boldsymbol\alpha\) 是常数，不影响协方差）。 第一步，两部分不相关（假设 2），所以总协方差 = 各自协方差之和，没有交叉项。 第二步，\(\mathrm{Cov}(\boldsymbol\beta\boldsymbol f_t)=\boldsymbol\beta\,\mathrm{Cov}(\boldsymbol f_t)\boldsymbol\beta'=\boldsymbol\beta\boldsymbol\Sigma_f\boldsymbol\beta'\)，用的是"线性变换的协方差 = 左乘 + 右乘转置"。 第三步，\(\mathrm{Cov}(\boldsymbol\epsilon_t)=\boldsymbol D\) 是对角阵（假设 3）。 逐元素看：\(\mathrm{Cov}(r_i,r_j)=\boldsymbol\beta_i'\boldsymbol\Sigma_f\boldsymbol\beta_j\)（\(i\ne j\)），即两只股票的协方差只通过共同因子产生。单因子时就是 CFA 里的 \(\beta_i\beta_j\sigma_m^2\)。

金融直觉：组合 \(\boldsymbol w\) 的风险 \(\boldsymbol w'\boldsymbol\Sigma_r\boldsymbol w=(\boldsymbol\beta'\boldsymbol w)'\boldsymbol\Sigma_f(\boldsymbol\beta'\boldsymbol w)+\sum_iw_i^2\sigma_i^2\)。第一项只依赖组合的因子暴露 \(\boldsymbol\beta'\boldsymbol w\)，第二项是特质风险，随持仓分散按 \(1/k\) 的速度变小。这正是风险报告里"因子风险 + 特质风险"的来源，也解释了为什么充分分散的组合只剩因子风险。

9.1.2 时间序列形式与多元线性回归

对第 \(i\) 个资产，把 \(T\) 期堆叠起来：\(\boldsymbol R_i=\alpha_i\boldsymbol 1_T+\boldsymbol F\boldsymbol\beta_i'+\boldsymbol E_i\)，\(\boldsymbol F\) 为 \(T\times m\) 因子矩阵。令 \(\boldsymbol g_t=(1,\boldsymbol f_t')'\)，\(\boldsymbol\xi=[\boldsymbol\alpha,\boldsymbol\beta]\)，所有资产合并为

\(\boldsymbol R\) 为 \(T\times k\)，\(\boldsymbol G\) 为 \(T\times(m+1)\)。因子可观测时这就是多元线性回归的特例（特殊之处在于误差协方差是对角的）。

9.2 宏观经济因子模型

9.2.1 估计

因子可观测（宏观变量或市场收益），对 (9.4) 做最小二乘：

\(\hat\sigma_i^2\) 取 \(\hat{\boldsymbol E}'\hat{\boldsymbol E}/(T-m-1)\) 的第 \(i\) 个对角元，\(R_i^2=1-[\hat{\boldsymbol E}'\hat{\boldsymbol E}]_{ii}/[\boldsymbol R'\boldsymbol R]_{ii}\)（\(\boldsymbol R\) 已去均值）。LS 没有施加"特质因子互不相关"的约束，一般不是有效估计，但施加约束很繁琐，通常忽略。模型检验：看残差相关矩阵 \(\hat{\boldsymbol E}'\hat{\boldsymbol E}/(T-m-1)\) 的非对角元是否接近 0。

9.2.2 单因子市场模型

市场模型（Sharpe 1970）：

\(r_{it}\)、\(r_{mt}\) 为超额收益，\(\beta_i\) 就是通常说的股票 β。

原书例子。 13 只美股（AA、AGE、CAT、F、FDX、GM、HPQ、KMB、MEL、NYT、PG、TRB、TXN）1990–2003 年月度超额收益，\(k=13\)，\(T=168\)，S&P 500 为市场。部分结果（\(\hat\beta\) / \(\hat\sigma_i\) / \(R^2\)）：

金融股和高科技股 β 与 \(R^2\) 较高，KMB、PG 较低。\(R^2\) 在 0.09–0.41 之间：市场只能解释不到一半的个股波动。模型隐含相关矩阵 \(\hat\sigma_m^2\hat{\boldsymbol\beta}\hat{\boldsymbol\beta}'+\hat{\boldsymbol D}\) 的元素多在 0.1–0.4，而样本相关矩阵中 Cor(AA,CAT)、Cor(F,GM)、Cor(HPQ,TXN) 都约 0.6——同行业的高相关单因子模型捕捉不到。残差相关也印证这一点：Cor(CAT,AA)=0.45、Cor(GM,F)=0.48，遗漏了行业因子。

用 GMVP 比较协方差估计。 全局最小方差组合

只依赖协方差矩阵，是比较协方差估计的实用工具。原书中模型 GMVP 与样本 GMVP 都重仓 KMB、NYT、PG（低 β 低波动股），但对 TRB 的权重差异明显（0.14 对 0.017）。

推导拆解：GMVP 公式的来历。 第一步，拉格朗日函数 \(\mathcal L=\boldsymbol\omega'\boldsymbol\Sigma\boldsymbol\omega-2\lambda(\boldsymbol\omega'\boldsymbol 1-1)\)（乘 2 只是为了让后面的式子干净）。 第二步，对 \(\boldsymbol\omega\) 求偏导：二次型 \(\boldsymbol\omega'\boldsymbol\Sigma\boldsymbol\omega\) 的梯度是 \(2\boldsymbol\Sigma\boldsymbol\omega\)（\(\boldsymbol\Sigma\) 对称），线性项 \(\boldsymbol\omega'\boldsymbol 1\) 的梯度是 \(\boldsymbol 1\)。令梯度为零：\(2\boldsymbol\Sigma\boldsymbol\omega-2\lambda\boldsymbol 1=\boldsymbol 0\)，得 \(\boldsymbol\omega=\lambda\boldsymbol\Sigma^{-1}\boldsymbol 1\)。 第三步，代入约束 \(\boldsymbol 1'\boldsymbol\omega=1\) 求 \(\lambda\)：\(\lambda\boldsymbol 1'\boldsymbol\Sigma^{-1}\boldsymbol 1=1\)，于是 \(\lambda=1/(\boldsymbol 1'\boldsymbol\Sigma^{-1}\boldsymbol 1)\)，代回即得。 为什么用它比较协方差估计？GMVP 不需要预期收益，避开了收益预测的巨大误差；同时它要用到 \(\boldsymbol\Sigma^{-1}\)，对协方差估计误差非常敏感——估计越差，样本外实现的波动越大。所以"样本外 GMVP 波动"是协方差估计质量的直接考题。

9.2.3 多因子宏观模型

Chen, Roll & Ross (1986) 用宏观变量的意外变化（surprises）作因子。为什么要用意外？因为宏观变量本身高度可预测（有强自相关），而市场只对未被预期的部分做出反应。意外序列可以取 VAR 模型的残差。

原书例子。 CPI 增长率与就业增长率（季调，1975–2003），BIC 选 VAR(3)，取 1990–2003 年的残差作为两个因子，对 13 只股票做回归。所有股票对 CPI 意外的 β 都为负（约 −14 到 −2），经济上合理；但 \(R^2\) 全都低于 0.06，模型隐含相关矩阵接近单位阵，残差相关矩阵与原始收益相关矩阵几乎相同。教训：宏观因子经济含义清楚，但对个股月收益的解释力通常很弱。

9.3 基本面因子模型

用可观测的公司属性——行业、市值、账面市值比、动量、波动率等——构造共同因子。两种方法的"已知 / 未知"恰好相反。

9.3.1 BARRA 方法：暴露已知，截面回归求因子收益

BARRA 方法（Bar Rosenberg 创立，Grinold & Kahn 2000）把观测到的基本面直接当作载荷 \(\boldsymbol\beta\)，在每个时点用截面回归估计因子实现值 \(\boldsymbol f_t\)。设超额收益已去均值，

这是 \(k\) 个观测、\(m\) 个未知数的回归。由于 \(\mathrm{Cov}(\boldsymbol\epsilon_t)=\boldsymbol D\) 异方差，有效估计是加权最小二乘：

白话解释：这里的回归方向与 CFA 里习惯的相反。平常做 beta 回归，是对一只股票、沿时间做回归，估出 beta；Barra 是在某一天、沿股票截面做回归——"因变量"是这一天 \(k\) 只股票的收益，"自变量"是它们已知的暴露（行业哑变量、市值得分等），估出来的"系数"\(\boldsymbol f_t\) 就是这一天各因子的收益。例如只有一个"科技行业"哑变量时，回归系数就是"科技股比其他股票多赚了多少"。 为什么要加权：特质波动大的股票（如小盘科技股）收益里噪声多，对"这一天科技因子赚了多少"提供的信息少，所以权重取 \(1/\sigma_i^2\)，即 \(\boldsymbol D^{-1}\)。这就是 (9.7) 中 \(\boldsymbol D^{-1}\) 的作用。

\(\boldsymbol D\) 未知，于是两步法：

1. 每个 \(t\) 做 OLS：\(\hat{\boldsymbol f}_{t,o}=(\boldsymbol\beta'\boldsymbol\beta)^{-1}\boldsymbol\beta'\tilde{\boldsymbol r}_t\)（一致但非有效），残差 \(\hat{\boldsymbol\epsilon}_{t,o}\)；合并所有时点估计 \(\hat{\boldsymbol D}_o=\mathrm{diag}\{\frac1{T-1}\sum_t\hat{\boldsymbol\epsilon}_{t,o}\hat{\boldsymbol\epsilon}_{t,o}'\}\)。
2. 代入得可行 GLS 估计

再由 GLS 残差得 \(\hat{\boldsymbol D}_g\)，由 \(\{\hat{\boldsymbol f}_{t,g}\}\) 的样本协方差得 \(\hat{\boldsymbol\Sigma}_f\)，最终 \(\widehat{\mathrm{Cov}}(\boldsymbol r_t)=\boldsymbol\beta\hat{\boldsymbol\Sigma}_f\boldsymbol\beta'+\hat{\boldsymbol D}_g\)。

行业因子模型示例。 10 只股票 1990–2003 年月度超额收益：金融（AGE、C、MWD、MER）、科技（DELL、HPQ、IBM）、其他（AA、CAT、PG）。三个因子代表三个行业，\(\beta_{ij}=1\) 当资产 \(i\) 属于行业 \(j\)，否则为 0（如 IBM 的 \(\boldsymbol\beta_i=(0,1,0)'\)）。由于 β 是哑变量，OLS 因子估计就是各时点行业内的等权平均超额收益，特质收益就是个股减去行业平均。GLS 则不再等权：金融因子对 AGE、C、MWD、MER 的权重为 0.187、0.255、0.259、0.300；科技因子对 DELL、HPQ、IBM 为 0.227、0.402、0.371——特质波动大的 DELL 被降权。

模型隐含相关矩阵中，行业内相关普遍高于样本值（金融股 0.7–0.8；CAT–PG 样本仅 0.1，模型为 0.6）。原因是"其他"这个大杂烩行业把 CAT 和 PG 强行视为同一因子驱动——行业划分的质量直接决定模型质量。

因子模拟组合（factor mimicking portfolio）。单因子情形下，WLS 估计有一个漂亮的组合解释。求解

拉格朗日条件给出 \(\boldsymbol\omega'=(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)^{-1}\boldsymbol\beta'\boldsymbol D^{-1}\)，所以 \(\hat f_t=\boldsymbol\omega'\boldsymbol r_t\) 恰好是这个组合的收益：在所有对该因子暴露为 1 的组合中，特质风险最小的那一个。多因子时，GLS 权重矩阵 \(\boldsymbol H=(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)^{-1}\boldsymbol\beta'\boldsymbol D^{-1}\) 满足 \(\boldsymbol H\boldsymbol\beta=\boldsymbol I\)：第 \(j\) 行是对因子 \(j\) 暴露为 1、对其他因子暴露为 0 的组合，也就是实务中说的纯因子组合。实践中超额收益样本均值常与 0 无显著差异，拟合前可不去均值。

推导拆解：拉格朗日条件的细节（单因子情形）。 第一步，\(\mathcal L=\tfrac12\boldsymbol\omega'\boldsymbol D\boldsymbol\omega-\lambda(\boldsymbol\omega'\boldsymbol\beta-1)\)，对 \(\boldsymbol\omega\) 求偏导令其为零：\(\boldsymbol D\boldsymbol\omega=\lambda\boldsymbol\beta\)，所以 \(\boldsymbol\omega=\lambda\boldsymbol D^{-1}\boldsymbol\beta\)——每只股票的权重正比于 \(\beta_i/\sigma_i^2\)。 第二步，代入约束 \(\boldsymbol\beta'\boldsymbol\omega=1\)：\(\lambda=1/(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)\)。 第三步，比较 (9.7)：单因子时 \(\hat f_t=\boldsymbol\beta'\boldsymbol D^{-1}\tilde{\boldsymbol r}_t/(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)\)，恰好等于 \(\boldsymbol\omega'\tilde{\boldsymbol r}_t\)。 多因子时 \(\boldsymbol H\boldsymbol\beta=(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)^{-1}(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)=\boldsymbol I\)，直接相乘就能验证。

金融直觉：于是"截面回归估出的因子收益"和"一个可以真实持有的组合的收益"是同一个数。比如价值因子收益，就是一个"价值暴露为 1、行业和其他风格暴露都为 0、特质风险尽量小"的多空组合当天的收益。这让因子收益可以被交易、被归因，也让你能检查它是否需要过高的杠杆或换手。

9.3.2 Fama–French 方法：构造组合得因子收益，时间序列回归求暴露

Fama & French (1992, 1993) 反过来做：

1. 按某个基本面（如账面市值比 B/M）对股票排序，做多排名最高的一组、做空最低的一组，构成对冲组合（hedge portfolio），其收益就是该因子在 \(t\) 时的实现值 \(f_{jt}\)；
2. 对每个资产做时间序列回归，估计其对各因子的 β。

三因子为：市场超额收益；SMB（small minus big，小盘减大盘，市值按中位数分两组）；HML（high minus low，高 B/M 减低 B/M，按 30%/40%/30% 分三组）。（原书称"top quintile"，与 FF 原文的分组不符，以 FF 原文为准。）

"因子个数"要谨慎。 FF 的三个因子是三个基本面；若把它们线性组合成一个新属性，形式上也可以称为单因子模型。统计因子模型中的因子个数则相对明确。

9.3.3 三类模型对照

9.4 主成分分析

9.4.1 理论

PCA 用少数线性组合解释 \(k\) 维收益 \(\boldsymbol r\) 的协方差结构。设 \(y_i=\boldsymbol w_i'\boldsymbol r\)，\(\boldsymbol w_i'\boldsymbol w_i=1\)（若 \(\boldsymbol r\) 是简单收益，\(y_i\) 就是一个组合的收益）。则

定义：第 1 主成分在 \(\boldsymbol w_1'\boldsymbol w_1=1\) 下最大化 \(\mathrm{Var}(y_1)\)；第 \(i\) 主成分在单位长度且与前 \(i-1\) 个主成分不相关的约束下最大化方差。

结论 9.1：设 \(\boldsymbol\Sigma_r\) 的特征对为 \((\lambda_i,\boldsymbol e_i)\)，\(\lambda_1\ge\cdots\ge\lambda_k\ge0\)，则第 \(i\) 主成分 \(y_i=\boldsymbol e_i'\boldsymbol r\)，\(\mathrm{Var}(y_i)=\lambda_i\)，不同主成分互不相关，并且

证明思路：拉格朗日函数 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w-\lambda(\boldsymbol w'\boldsymbol w-1)\) 对 \(\boldsymbol w\) 求导得 \(\boldsymbol\Sigma\boldsymbol w=\lambda\boldsymbol w\)，此时方差 \(=\boldsymbol w'\boldsymbol\Sigma\boldsymbol w=\lambda\)，所以取最大特征值对应的特征向量；加入与 \(\boldsymbol e_1\) 正交的约束后重复，依次得到其余特征向量（对称矩阵的特征向量天然正交，见第 01 册）。

于是第 \(i\) 主成分解释总方差的比例为 \(\lambda_i/\sum_j\lambda_j\)。用相关矩阵时 \(\mathrm{tr}(\boldsymbol\rho_r)=k\)，比例为 \(\lambda_i/k\)。副产品：若 \(\lambda_k=0\)，则 \(\boldsymbol e_k'\boldsymbol r\) 是常数，分量间存在精确线性关系，可以降维。

推导拆解：把证明思路逐步展开。 第一步，目标：在 \(\boldsymbol w'\boldsymbol w=1\) 下最大化 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w\)。约束是必要的，否则把 \(\boldsymbol w\) 放大一倍方差就放大四倍，问题无界。 第二步，对 \(\mathcal L=\boldsymbol w'\boldsymbol\Sigma\boldsymbol w-\lambda(\boldsymbol w'\boldsymbol w-1)\) 求梯度：\(2\boldsymbol\Sigma\boldsymbol w-2\lambda\boldsymbol w=\boldsymbol 0\)，即 \(\boldsymbol\Sigma\boldsymbol w=\lambda\boldsymbol w\)。所以极值点只能是特征向量，乘子 \(\lambda\) 就是特征值。 第三步，在特征向量处，方差 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w=\boldsymbol w'(\lambda\boldsymbol w)=\lambda\boldsymbol w'\boldsymbol w=\lambda\)。要方差最大，就选最大的特征值。 第四步，(9.13)：\(\mathrm{tr}(\boldsymbol\Sigma)\) 按定义是对角元（各资产方差）之和；线性代数的一个基本事实是迹等于特征值之和。所以"各资产方差之和"="各主成分方差之和"，PCA 只是把总方差重新分配到互不相关的方向上，没有增减。 数值例子（练习 4）：\(\boldsymbol\Sigma=\begin{bmatrix}2&1\\1&2\end{bmatrix}\)，迹 4 = 3 + 1。第一主成分 \((r_1+r_2)/\sqrt2\) 方差 3，占 75%；第二主成分 \((r_1-r_2)/\sqrt2\) 方差 1，占 25%。

金融直觉：如果 \(\boldsymbol r\) 是收益，每个主成分就是一个"组合"（权重向量长度为 1，不是权重和为 1）。第一主成分是"波动最大的那个方向"，在股票里几乎总是所有股票同向持有的"市场组合"；后面的主成分则是多空组合，例如"做多科技、做空金融"。

9.4.2 经验 PCA

用样本协方差 \(\hat{\boldsymbol\Sigma}_r=\frac1{T-1}\sum(\boldsymbol r_t-\bar{\boldsymbol r})(\boldsymbol r_t-\bar{\boldsymbol r})'\) 或样本相关矩阵 \(\hat{\boldsymbol\rho}_r=\hat{\boldsymbol S}^{-1}\hat{\boldsymbol\Sigma}_r\hat{\boldsymbol S}^{-1}\) 做特征分解。

协方差还是相关？ 协方差矩阵 PCA 会被高波动资产主导；各资产波动差异大或量纲不同时应使用相关矩阵。主成分的符号是任意的（特征向量乘 −1 仍是特征向量），解释时要留意。

例 9.1（五只股票）。 IBM、HPQ、INTC、JPM、BAC 月度对数收益，1990–2008，228 个观测。样本相关矩阵（下三角）为

相关阵 PCA 中，\(\boldsymbol e_1=(0.428,0.460,0.451,0.479,0.416)'\) 近似等权——市场成分；\(\boldsymbol e_2=(0.341,0.356,0.385,-0.469,-0.623)'\) 是科技股减金融股——行业成分。两者共解释约 74%。

碎石图（scree plot）：画 \(\hat\lambda_i\) 对 \(i\)，找"肘部"——之后的特征值较小且相近。例 9.1 取 2 个主成分。除非被舍弃的特征值恰好为 0，取前几个主成分只是近似。

9.5 统计因子分析

9.5.1 正交因子模型

当因子和载荷都不可观测时，

这已不是回归模型。正交因子模型假设：\(E(\boldsymbol f_t)=\boldsymbol 0\)，\(\mathrm{Cov}(\boldsymbol f_t)=\boldsymbol I_m\)；\(\mathrm{Cov}(\boldsymbol\epsilon_t)=\boldsymbol D\) 对角；\(\boldsymbol f_t\) 与 \(\boldsymbol\epsilon_t\) 独立。于是

逐元素看，\(\mathrm{Var}(r_{it})=c_i^2+\sigma_i^2\)，其中 \(c_i^2=\sum_j\beta_{ij}^2\) 称为共同度（communality），\(\sigma_i^2\) 称为特殊方差（specific variance / uniqueness）。

不唯一性。 对任意 \(m\times m\) 正交阵 \(\boldsymbol P\)，\(\boldsymbol\beta^*=\boldsymbol\beta\boldsymbol P\)、\(\boldsymbol f_t^*=\boldsymbol P'\boldsymbol f_t\) 满足全部假设且给出同一个 \(\boldsymbol\Sigma_r\)。这既是缺点（载荷含义任意），也是优点（可以旋转出更好解释的载荷）。

推导拆解：两处一行就带过的结论。 (a) \(\mathrm{Cov}(\boldsymbol r_t,\boldsymbol f_t)=E[(\boldsymbol\beta\boldsymbol f_t+\boldsymbol\epsilon_t)\boldsymbol f_t']=\boldsymbol\beta E[\boldsymbol f_t\boldsymbol f_t']+E[\boldsymbol\epsilon_t\boldsymbol f_t']=\boldsymbol\beta\boldsymbol I_m+\boldsymbol 0=\boldsymbol\beta\)。所以在正交因子模型里，载荷就是"资产与因子的协方差"；若资产也标准化为单位方差，载荷就是相关系数。 (b) 正交阵 \(\boldsymbol P\)（满足 \(\boldsymbol P\boldsymbol P'=\boldsymbol I\)，几何上是旋转或翻转）：\(\boldsymbol\beta^*\boldsymbol f_t^*=\boldsymbol\beta\boldsymbol P\boldsymbol P'\boldsymbol f_t=\boldsymbol\beta\boldsymbol f_t\)，收益完全不变；\(\mathrm{Cov}(\boldsymbol f_t^*)=\boldsymbol P'\boldsymbol I\boldsymbol P=\boldsymbol I\)，假设仍成立。

金融直觉：两个因子"市场"和"科技减金融"，与旋转后的"科技板块"和"金融板块"两个因子，能解释完全相同的协方差矩阵。数据本身无法区分这两种说法，选哪种取决于哪种更好解释、更好用。这就是 9.5.3 旋转的意义。

9.5.2 估计

主成分法。 不需要正态假设，也不需要预设因子数。取样本协方差（或相关）矩阵前 \(m\) 个特征对，

\(\hat\sigma_i^2=\hat\sigma_{ii}-\sum_j\hat\beta_{ij}^2\)。近似误差 \(\hat{\boldsymbol\Sigma}_r-(\hat{\boldsymbol\beta}\hat{\boldsymbol\beta}'+\hat{\boldsymbol D})\) 的元素平方和不超过 \(\hat\lambda_{m+1}^2+\cdots+\hat\lambda_k^2\)。增加 \(m\) 时已有载荷不变。

白话解释：为什么载荷是 \(\sqrt{\hat\lambda_j}\hat{\boldsymbol e}_j\)？由特征分解 \(\hat{\boldsymbol\Sigma}=\sum_{j=1}^k\hat\lambda_j\hat{\boldsymbol e}_j\hat{\boldsymbol e}_j'\)，只保留前 \(m\) 项，并把每一项写成 \((\sqrt{\hat\lambda_j}\hat{\boldsymbol e}_j)(\sqrt{\hat\lambda_j}\hat{\boldsymbol e}_j)'\)，前 \(m\) 项之和正好是 \(\hat{\boldsymbol\beta}\hat{\boldsymbol\beta}'\)。开平方是因为正交因子模型要求因子方差为 1：主成分 \(y_j\) 的方差是 \(\lambda_j\)，把它除以 \(\sqrt{\lambda_j}\) 变成单位方差，这个尺度就转移到了载荷上。 剩下被丢掉的 \(\sum_{j>m}\hat\lambda_j\hat{\boldsymbol e}_j\hat{\boldsymbol e}_j'\) 中，对角线部分由 \(\hat{\boldsymbol D}\) 补上，非对角部分就是近似误差。被丢掉的特征值越小，误差越小。

极大似然法。 若 \(\boldsymbol f_t,\boldsymbol\epsilon_t\) 联合正态，则 \(\boldsymbol r_t\sim N(\boldsymbol\mu,\boldsymbol\beta\boldsymbol\beta'+\boldsymbol D)\)，在唯一性约束 \(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta=\) 对角阵下极大化似然。需预设 \(m\)，并可用修正 LR 统计量检验 \(m\) 个因子是否足够：

自由度是"协方差矩阵的自由参数"减去"因子模型的自由参数"。

推导拆解：自由度的点数。无约束的 \(k\times k\) 协方差矩阵有 \(k(k+1)/2\) 个自由参数。因子模型有 \(\boldsymbol\beta\) 的 \(km\) 个和 \(\boldsymbol D\) 的 \(k\) 个，但旋转不唯一性意味着其中有 \(m(m-1)/2\) 个是"多余"的（唯一性约束 \(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta\) 为对角阵正好施加这么多个条件）。相减： \(\frac{k(k+1)}2-\left[km+k-\frac{m(m-1)}2\right]=\frac{(k-m)^2-k-m}2\)。 验算例 9.4：\(k=10\)，\(m=2\) 得 \((64-12)/2=26\)；\(m=3\) 得 \((49-13)/2=18\)，与正文一致。 读法：LR 检验的原假设是"\(m\) 个因子够了"，p 值大（如例 9.4 的 0.089）表示不能拒绝，即 \(m\) 个因子可以接受。从小到大逐个试，第一个不被拒绝的 \(m\) 就是选择。

9.5.3 因子旋转

Varimax 准则（Kaiser 1958）：令 \(\tilde\beta^*_{ij}=\beta^*_{ij}/c_i\)，选正交阵 \(\boldsymbol P\) 最大化

即让每一列载荷平方的"方差"尽量大，使每个因子只在少数资产上有大载荷、其余接近 0。旋转不改变共同度与特殊方差，只是辅助解释的工具，有时有用有时无信息；还有 quartimax 等其他准则。

9.5.4 应用

先检验无序列相关假设；若有，可先拟合 VARMA 再对残差做因子分析。很多时候残差相关矩阵与原始数据很接近，可直接分析。

例 9.2（同例 9.1 的五只股票）。 \(Q_5(1)=39.99\)、\(Q_5(5)=160.60\)、\(Q_5(10)=293.04\)，p 值 0.029、0.017、0.032，1% 水平下忽略。相关阵 ML 估计、2 因子：

两因子共解释约 60%。旋转后科技股重载 \(f_1^*\)，金融股重载 \(f_2^*\)，行业区分一目了然。IBM 共同度仅 0.387，特质成分大。

例 9.3（美国国债指数，期限 30、20、10、5、1 年）。 主成分法 2 因子：\(f_1\) 载荷 0.952、0.954、0.956、0.955、0.800，近似相等——水平因子；\(f_2\) 载荷 0.253、0.240、0.140、−0.142、−0.585，随期限单调且和约为 0——斜率因子（长债减短债）。两因子解释 95.7% 的变异。旋转后变成"长端因子"和"短端因子"。ML 法结果类似（合计 90.5%）。

例 9.4（10 只股票，同 9.3.1 节）。 2 因子被 LR 检验拒绝（\(LR(2)=72.96\)，\(\chi^2_{26}\)，p≈0）；3 因子 \(\chi^2=26.48\)、自由度 18、p=0.089，5% 水平可接受。载荷：因子 1 集中在金融股（0.68–0.82），因子 2 在科技股和 AA，因子 3 主要是 CAT（0.970）与 AA。拟合相关矩阵比 BARRA 行业模型更接近样本值（如 CAT–PG 0.1）。注意 CAT 的特殊方差估为 0，这是 Heywood 情形（ML 估计落在参数边界），提示该因子几乎由 CAT 一只股票定义，解释时要谨慎。

9.6 渐近主成分分析（\(k>T\)）

9.6.1 方法

前面的 PCA 需要 \(k<T\)，否则 \(k\times k\) 样本协方差矩阵奇异。Connor & Korajczyk (1986, 1988) 的渐近主成分分析（APCA）利用 \(k\to\infty\) 的渐近理论，转而分解 \(T\times T\) 矩阵

前 \(m\) 个特征向量就是因子 \(\boldsymbol f_t\) 的估计（第 \(t\) 个分量对应第 \(t\) 期）。直觉：\(T\times T\) 矩阵 \(\boldsymbol X\boldsymbol X'\) 与 \(k\times k\) 矩阵 \(\boldsymbol X'\boldsymbol X\) 有相同的非零特征值（SVD 的对偶性，见第 01 册），\(k\gg T\) 时分解小矩阵快得多。仿照 BARRA 可以精化：用初始因子对每个资产做时间序列回归得 \(\hat\sigma_i^2\)，把收益按 \(\hat{\boldsymbol D}^{-1/2}\) 缩放后重新分解。

推导拆解：为什么两个矩阵特征值相同、特征向量能互相换算？设 \(\boldsymbol X'\boldsymbol X\boldsymbol v=\lambda\boldsymbol v\)（\(\boldsymbol v\) 是 \(k\) 维，"资产权重方向"）。两边左乘 \(\boldsymbol X\)：\(\boldsymbol X\boldsymbol X'(\boldsymbol X\boldsymbol v)=\lambda(\boldsymbol X\boldsymbol v)\)。所以 \(\boldsymbol u=\boldsymbol X\boldsymbol v\) 是 \(\boldsymbol X\boldsymbol X'\) 的特征向量，特征值同为 \(\lambda\)。而 \(\boldsymbol X\boldsymbol v\) 的第 \(t\) 个分量正是"按权重 \(\boldsymbol v\) 持有的组合在第 \(t\) 期的收益"，也就是主成分的时间序列。 白话解释：传统 PCA 先找"权重"再算"因子时间序列"；APCA 反过来，直接在 \(T\times T\) 的小矩阵上找出因子的时间序列，再对每只股票回归得到载荷。\(k=3000\) 只股票、\(T=60\) 个月时，前者要分解 3000×3000 的矩阵（而且它的秩最多只有 60，大部分特征值是 0），后者只要分解 60×60 的矩阵。

9.6.2 因子个数

• Connor & Korajczyk (1993)：若 \(m\) 正确，从 \(m\) 增加到 \(m+1\) 时特质误差的截面方差不应显著下降。
• Bai & Ng (2002)：令 \(\hat\sigma^2(m)\) 为用 \(m\) 个 APCA 因子回归后残差方差的截面平均，\(M\) 为预设上限，

\(P_{kT}=\min(\sqrt k,\sqrt T)\)，取使准则最小的 \(m\)。

原书实例。 40 只 NASDAQ/NYSE 最活跃股票，2001–2003 月度收益，\(k=40>T=36\)。Connor–Korajczyk 法选 1 个因子（回归 \(R^2\) 中位数 0.487），Bai–Ng 法选 6 个（\(R^2\) 中位数 0.695，6 因子累计解释约 89.4%）。两种方法结论差异很大，\(T=36\) 时尤其不稳定，实务中要结合经济解释和样本外表现。

量化实战

应用场景

多因子风险模型。 BARRA 方法就是 Barra USE/CNE 等商业风险模型的原型：行业哑变量 + 风格暴露（市值、动量、波动率、估值、流动性等，截面标准化后）作为 \(\boldsymbol\beta\)，每期截面 WLS 得到因子收益，再估 \(\boldsymbol\Sigma_f\)（通常用 EWMA，见第 10 章）和特质方差 \(\boldsymbol D\)。组合风险 \(\boldsymbol w'(\boldsymbol\beta\boldsymbol\Sigma_f\boldsymbol\beta'+\boldsymbol D)\boldsymbol w\) 可以分解为各因子贡献和特质贡献，用于风险预算与归因。实务中截面回归的权重常取市值平方根而不是 \(1/\hat\sigma_i^2\)，以兼顾稳健性和对大市值股票的关注；同时加入国家（市场）因子时需对行业因子收益施加"市值加权和为零"的约束以消除共线性。

因子研究。 Fama–French 排序法是构造 SMB、HML 以及各类异象因子收益的标准做法；截面回归得到的"纯因子收益"（因子模拟组合）剔除了其他因子的干扰，常用于因子绩效归因，与 IC 分析互补。

组合优化。 \(k\) 接近或大于 \(T\) 时样本协方差矩阵病态甚至奇异，直接用于均值–方差优化会放大估计误差。因子模型协方差是最常用的"结构化收缩"，GMVP 的样本外波动是检验协方差估计质量的简单基准。

统计套利。 PCA / APCA 提取的统计因子可用来构建市场中性组合：把个股收益对前 \(m\) 个主成分回归，残差（特质收益）累积成"特质价格"，再像第 08b 章的价差一样做均值回复交易。Bai–Ng 准则可用于确定剔除多少个因子。

固定收益。 债券收益的前两个主成分即利率曲线的水平与斜率（例 9.3），用于把久期/DV01 风险按主成分分解并对冲。

Python 示例

用一个统一的模拟器生成月度超额收益（单位 %）：市场因子 + 3 个行业因子 + 规模风格因子 + 特质噪声，每个行业 10 只股票。先把下面这段代码保存为工作目录下的 c09_sim.py：后面的示例一、二、三都以 from c09_sim import simulate 调用它，不保存会报 ModuleNotFoundError。

# 文件名：c09_sim.py（示例一至示例三共用）
import numpy as np

def simulate(T=240, n_ind=3, per_ind=10, seed=0):
    """月度超额收益（%）：市场 + 行业 + 规模风格 + 特质"""
    rng = np.random.default_rng(seed)
    k = n_ind * per_ind
    ind = np.repeat(np.arange(n_ind), per_ind)            # 行业归属
    beta_m = rng.uniform(0.6, 1.5, k)                      # 市场 beta
    size = rng.standard_normal(k); size = (size - size.mean()) / size.std()  # 规模暴露（标准化）
    sig_e = rng.uniform(4, 9, k)                           # 特质波动
    f_mkt = 0.6 + 4.5 * rng.standard_normal(T)
    f_ind = 3.0 * rng.standard_normal((T, n_ind))
    f_size = -0.2 + 2.0 * rng.standard_normal(T)           # 规模因子：负暴露（小盘）占优
    eps = rng.standard_normal((T, k)) * sig_e
    R = np.outer(f_mkt, beta_m) + f_ind[:, ind] + np.outer(f_size, size) + eps
    return dict(R=R, ind=ind, beta_m=beta_m, size=size, sig_e=sig_e,
                f_mkt=f_mkt, f_ind=f_ind, f_size=f_size)

示例一：市场模型、残差相关、GMVP 与宏观意外因子

import numpy as np
from statsmodels.tsa.api import VAR
from c09_sim import simulate

d = simulate(T=240, seed=1)
R, mkt, ind = d["R"], d["f_mkt"], d["ind"]
T, k = R.shape

# ---------- 1. 单因子市场模型：时间序列 OLS，式 (9.5) ----------
G = np.column_stack([np.ones(T), mkt])                    # G = [1, f_t]
xi = np.linalg.lstsq(G, R, rcond=None)[0]                 # (G'G)^{-1} G'R，2 x k
E = R - G @ xi
D_hat = np.diag(E.T @ E / (T - 2))
Rc = R - R.mean(0)
r2 = 1 - np.diag(E.T @ E) / np.diag(Rc.T @ Rc)
print("beta 估计 vs 真值（前 5 只）:", xi[1, :5].round(2), d["beta_m"][:5].round(2))
print(f"R^2 范围: {r2.min():.2f} ~ {r2.max():.2f}, 中位数 {np.median(r2):.2f}")

# 残差相关：同行业 vs 跨行业（检验"特质因子互不相关"假设）
C = np.corrcoef(E.T)
same = (ind[:, None] == ind[None, :]) & ~np.eye(k, dtype=bool)
diff = ind[:, None] != ind[None, :]
print(f"残差平均相关：同行业 {C[same].mean():.3f}，跨行业 {C[diff].mean():.3f}")

# ---------- 2. 全局最小方差组合（GMVP）：模型协方差 vs 样本协方差 ----------
def gmvp(S):
    w = np.linalg.solve(S, np.ones(len(S)))
    return w / w.sum()
S_model = np.var(mkt, ddof=1) * np.outer(xi[1], xi[1]) + np.diag(D_hat)
S_samp = np.cov(R.T)
w_m, w_s = gmvp(S_model), gmvp(S_samp)
print("GMVP 权重相关性（模型 vs 样本）:", np.corrcoef(w_m, w_s)[0, 1].round(3),
      " 卖空权重个数:", (w_m < 0).sum(), (w_s < 0).sum())

# ---------- 3. 宏观因子：用 VAR 残差作为"意外"（Chen-Roll-Ross 思路） ----------
rng = np.random.default_rng(5)
macro = np.zeros((T + 100, 2))                            # (通胀增速, 就业增速)
A = np.array([[0.6, 0.1], [0.0, 0.5]])
for t in range(1, T + 100):
    macro[t] = 0.2 + A @ macro[t-1] + rng.standard_normal(2) * [0.3, 0.2]
macro = macro[100:]
surprise = VAR(macro).fit(maxlags=4, ic="bic").resid      # 去掉可预测部分
Tm = len(surprise)
# 让股票对"通胀意外"有负暴露（真值 -4），对就业意外无暴露
Rm = R[-Tm:] + np.outer(surprise[:, 0], -4.0 * np.ones(k))
Gm = np.column_stack([np.ones(Tm), surprise])
bm = np.linalg.lstsq(Gm, Rm, rcond=None)[0]
Em = Rm - Gm @ bm
r2m = 1 - np.diag(Em.T @ Em) / np.diag((Rm - Rm.mean(0)).T @ (Rm - Rm.mean(0)))
se = np.sqrt(np.diag(np.linalg.inv(Gm.T @ Gm)))[1:, None] * np.sqrt((Em**2).sum(0) / (Tm - 3))
print(f"通胀意外 beta 均值 {bm[1].mean():.2f}（真值 -4），平均 t 值 {np.mean(bm[1]/se[0]):.2f}；"
      f"就业意外 beta 均值 {bm[2].mean():.2f}，平均 t 值 {np.mean(bm[2]/se[1]):.2f}")
print(f"宏观因子模型 R^2 中位数 {np.median(r2m):.3f}")

输出：

beta 估计 vs 真值（前 5 只）: [0.94 1.51 0.69 1.56 0.76] [1.06 1.46 0.73 1.45 0.88]
R^2 范围: 0.08 ~ 0.54, 中位数 0.27
残差平均相关：同行业 0.174，跨行业 -0.001
GMVP 权重相关性（模型 vs 样本）: 0.933  卖空权重个数: 12 11
通胀意外 beta 均值 -2.11（真值 -4），平均 t 值 -1.06；就业意外 beta 均值 -0.11，平均 t 值 -0.05
宏观因子模型 R^2 中位数 0.007

单因子模型的 \(R^2\) 中位数 0.27，与原书 0.09–0.41 的范围相当；同行业残差平均相关 0.17、跨行业约 0——和原书 Cor(GM,F)=0.48 一样，说明遗漏了行业因子。宏观因子部分即使设定了真实暴露 −4，估计出的 β 也只有 −2.1、平均 t 值约 −1，\(R^2\) 不到 1%：通胀意外的方差相对个股特质波动太小，这就是原书中宏观因子模型"方向合理、解释力极弱"的原因。

示例二：BARRA 两步 GLS、因子模拟组合与 Fama–French 排序

import numpy as np
from c09_sim import simulate

d = simulate(T=240, seed=1)
R, ind, size = d["R"], d["ind"], d["size"]
T, k = R.shape

# ---------- 1. BARRA：暴露已知（3 个行业哑变量 + 规模），逐期截面回归求因子收益 ----------
B = np.column_stack([(ind == j).astype(float) for j in range(3)] + [size])   # k x m
Rt = R - R.mean(0)                                         # 去均值超额收益（式 9.6）

# 第一步：逐期 OLS
F_ols = np.linalg.solve(B.T @ B, B.T @ Rt.T).T             # T x m
E_ols = Rt - F_ols @ B.T
D_o = E_ols.var(0, ddof=1)                                 # 合并全部时点估计特质方差
# 第二步：GLS / WLS，式 (9.8)
H = np.linalg.solve(B.T @ (B / D_o[:, None]), (B / D_o[:, None]).T)   # (B'D^-1B)^-1 B'D^-1，m x k
F_gls = Rt @ H.T
E_gls = Rt - F_gls @ B.T
D_g = E_gls.var(0, ddof=1)

true_ind = d["f_ind"] + np.outer(d["f_mkt"], [d["beta_m"][ind == j].mean() for j in range(3)])
for j, name in enumerate(["行业1", "行业2", "行业3"]):
    print(f"{name}因子：GLS 估计与(行业+平均beta×市场)真值相关 "
          f"{np.corrcoef(F_gls[:, j], true_ind[:, j])[0, 1]:.3f}")
print(f"规模因子：OLS 相关 {np.corrcoef(F_ols[:, 3], d['f_size'])[0, 1]:.3f}，"
      f"GLS 相关 {np.corrcoef(F_gls[:, 3], d['f_size'])[0, 1]:.3f}")

# GLS 权重 = 因子模拟组合：对自身因子暴露 1，对其他因子暴露 0
print("H @ B（应为单位阵）:\n", (H @ B).round(6) + 0.0)
w_ind1 = H[0, ind == 0]
print("行业1因子组合在本行业 10 只股票上的权重:", w_ind1.round(3),
      "\n  对应特质波动:", np.sqrt(D_o[ind == 0]).round(1))

# 模型协方差与样本协方差：同行业平均相关
S_f = np.cov(F_gls.T)
S_model = B @ S_f @ B.T + np.diag(D_g)
def avg_corr(S, mask):
    c = S / np.sqrt(np.outer(np.diag(S), np.diag(S)))
    return c[mask].mean()
same = (ind[:, None] == ind[None, :]) & ~np.eye(k, dtype=bool)
print(f"同行业平均相关：样本 {avg_corr(np.cov(R.T), same):.3f}，BARRA 模型 {avg_corr(S_model, same):.3f}")

# ---------- 2. Fama-French 方法：按规模排序构造对冲组合，再做时间序列回归 ----------
order = np.argsort(size)
small, big = order[:k // 3], order[-(k // 3):]
smb = R[:, small].mean(1) - R[:, big].mean(1)              # small minus big
print(f"SMB 与真实规模因子的相关 {np.corrcoef(smb, d['f_size'])[0, 1]:.3f}"
      f"（SMB 做多小盘，故为负相关）")
G = np.column_stack([np.ones(T), d["f_mkt"], smb])
coef = np.linalg.lstsq(G, R, rcond=None)[0]
print("SMB 载荷与规模暴露的截面相关:", np.corrcoef(coef[2], size)[0, 1].round(3))

输出：

行业1因子：GLS 估计与(行业+平均beta×市场)真值相关 0.935
行业2因子：GLS 估计与(行业+平均beta×市场)真值相关 0.941
行业3因子：GLS 估计与(行业+平均beta×市场)真值相关 0.947
规模因子：OLS 相关 0.851，GLS 相关 0.865
H @ B（应为单位阵）:
 [[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
行业1因子组合在本行业 10 只股票上的权重: [0.139 0.108 0.103 0.086 0.055 0.123 0.136 0.094 0.103 0.053] 
  对应特质波动: [5.2 4.8 6.5 6.8 6.7 5.2 4.9 6.5 6.8 8.3]
同行业平均相关：样本 0.384，BARRA 模型 0.439
SMB 与真实规模因子的相关 -0.806（SMB 做多小盘，故为负相关）
SMB 载荷与规模暴露的截面相关: -0.957

• 只用行业哑变量（没有单独的市场因子）时，行业因子吸收了"行业因子 + 行业平均 β × 市场"，估计与之相关 0.94 左右。
• \(\boldsymbol H\boldsymbol\beta=\boldsymbol I\) 验证了 GLS 权重就是纯因子组合：行业 1 因子组合对本行业暴露 1、对其他行业和规模暴露 0。权重随特质波动递减（特质波动 8.3 的股票只得 0.053），同时被规模暴露微调以保持规模中性——与原书 DELL 被降权的现象一致。
• BARRA 模型的同行业平均相关（0.44）高于样本值（0.38），重现了原书"行业模型高估行业内相关"的现象：真实数据中个股对市场的 β 不同，而行业哑变量假设同行业股票对行业因子的暴露都是 1。
• Fama–French 排序构造的 SMB 与真实规模因子相关 −0.81（SMB 做多小盘，而模拟中规模暴露为正表示大盘），时间序列回归得到的 SMB 载荷与真实暴露的截面相关 −0.96。

示例三：PCA、因子分析（LR 定阶 + varimax）、Bai–Ng 与样本外 GMVP

import numpy as np
from sklearn.decomposition import FactorAnalysis
from scipy.stats import chi2
from c09_sim import simulate

# ---------- 1. PCA（相关阵）：特征值、解释比例、第一主成分 ≈ 市场 ----------
d = simulate(T=240, seed=1)
R, ind = d["R"], d["ind"]
T, k = R.shape
corr = np.corrcoef(R.T)
lam, vec = np.linalg.eigh(corr)
lam, vec = lam[::-1], vec[:, ::-1]
print("前 6 个特征值:", lam[:6].round(2))
print("累计解释比例:", (np.cumsum(lam) / k)[:6].round(3))
e1 = vec[:, 0] * np.sign(vec[:, 0].sum())
print("第一特征向量（全部同号≈市场）: 最小 %.3f 最大 %.3f" % (e1.min(), e1.max()))
for j in (1, 2):
    print(f"第{j+1}特征向量按行业的平均载荷:",
          [round(float(vec[ind == g, j].mean()), 3) for g in range(3)])

# ---------- 2. 统计因子分析（ML）：LR 检验定因子数（式 9.20），再做 varimax 旋转 ----------
Z = (R - R.mean(0)) / R.std(0)
logdet_S = np.linalg.slogdet(np.corrcoef(Z.T))[1]
for m in range(2, 7):
    fa = FactorAnalysis(n_components=m, random_state=0).fit(Z)
    Sig = fa.get_covariance()
    LR = -(T - 1 - (2 * k + 5) / 6 - 2 * m / 3) * (logdet_S - np.linalg.slogdet(Sig)[1])
    dfree = ((k - m) ** 2 - k - m) / 2
    print(f"m={m}: LR={LR:7.1f}, df={dfree:.0f}, p={chi2.sf(LR, dfree):.3f}")
fa = FactorAnalysis(n_components=4, rotation="varimax", random_state=0).fit(Z)
L = fa.components_.T                                       # k x m 载荷
print("varimax 旋转后：各因子按行业的平均载荷 | 与规模暴露、市场 beta 的截面相关")
for j in range(4):
    print(f"  因子{j+1}:", [round(float(L[ind == g, j].mean()), 2) for g in range(3)],
          "|", round(float(np.corrcoef(L[:, j], d["size"])[0, 1]), 2),
          round(float(np.corrcoef(L[:, j], d["beta_m"])[0, 1]), 2))
comm = (L**2).sum(1)
print(f"共同度范围 {comm.min():.2f} ~ {comm.max():.2f}")

# ---------- 3. k > T：渐近主成分（APCA）+ Bai-Ng 准则选因子数 ----------
d2 = simulate(T=60, n_ind=4, per_ind=25, seed=2)           # k=100 只股票，T=60 个月
R2 = d2["R"]; T2, k2 = R2.shape
X = R2 - R2.mean(0)
Omega = X @ X.T / k2                                       # T x T 矩阵
ev, F = np.linalg.eigh(Omega); F = F[:, ::-1]
M = 10
def sigma2(m):
    if m == 0:
        return (X**2).mean()
    Fm = F[:, :m]
    resid = X - Fm @ np.linalg.lstsq(Fm, X, rcond=None)[0]
    return (resid**2).mean()
s2M = sigma2(M)
pen = (k2 + T2) / (k2 * T2)
Cp1 = [sigma2(m) + m * s2M * pen * np.log(k2 * T2 / (k2 + T2)) for m in range(M + 1)]
Cp2 = [sigma2(m) + m * s2M * pen * np.log(min(k2, T2)) for m in range(M + 1)]
print(f"k={k2}, T={T2}: Bai-Ng Cp1 选 m={int(np.argmin(Cp1))}，Cp2 选 m={int(np.argmin(Cp2))}"
      f"（真实共同因子：市场+4行业+规模，秩为 5）")

# ---------- 4. 样本外检验：不同协方差估计下 GMVP 的实现波动 ----------
d3 = simulate(T=360, seed=3)
R3 = d3["R"]; k3 = R3.shape[1]
def gmvp(S):
    w = np.linalg.solve(S, np.ones(len(S))); return w / w.sum()
def pca_cov(X, m):
    S = np.cov(X.T); l, v = np.linalg.eigh(S); l, v = l[::-1], v[:, ::-1]
    low = (v[:, :m] * l[:m]) @ v[:, :m].T
    return low + np.diag(np.maximum(np.diag(S - low), 1e-6))
def single_index(X):
    mkt = X.mean(1)
    b = np.cov(np.column_stack([mkt, X]).T)[0, 1:] / mkt.var(ddof=1)
    e = X - np.outer(mkt, b)
    return mkt.var(ddof=1) * np.outer(b, b) + np.diag(e.var(0, ddof=1))
res = {"样本协方差": [], "单指数模型": [], "PCA 4 因子": []}
win = 60                                                   # 只用 60 个月估计 30x30 协方差
for t0 in range(win, 360, 12):
    Xw, Xo = R3[t0 - win:t0], R3[t0:t0 + 12]
    for name, S in (("样本协方差", np.cov(Xw.T)), ("单指数模型", single_index(Xw)),
                    ("PCA 4 因子", pca_cov(Xw, 4))):
        res[name].append(Xo @ gmvp(S))
for name, r in res.items():
    r = np.concatenate(r)
    print(f"{name}: 样本外 GMVP 月波动 {r.std():.3f}%")

输出：

前 6 个特征值: [10.14  2.31  1.93  1.58  1.03  0.98]
累计解释比例: [0.338 0.415 0.479 0.532 0.567 0.599]
第一特征向量（全部同号≈市场）: 最小 0.115 最大 0.243
第2特征向量按行业的平均载荷: [-0.039, 0.187, -0.122]
第3特征向量按行业的平均载荷: [0.22, -0.015, -0.189]
m=2: LR=  794.1, df=376, p=0.000
m=3: LR=  507.0, df=348, p=0.000
m=4: LR=  347.2, df=321, p=0.150
m=5: LR=  305.8, df=295, p=0.321
m=6: LR=  265.5, df=270, p=0.566
varimax 旋转后：各因子按行业的平均载荷 | 与规模暴露、市场 beta 的截面相关
  因子1: [0.47, 0.17, 0.2] | 0.78 0.28
  因子2: [-0.19, -0.54, -0.2] | 0.4 -0.07
  因子3: [-0.16, -0.18, -0.57] | -0.23 -0.43
  因子4: [0.3, 0.12, 0.17] | -0.72 0.18
共同度范围 0.20 ~ 0.80
k=100, T=60: Bai-Ng Cp1 选 m=5，Cp2 选 m=5（真实共同因子：市场+4行业+规模，秩为 5）
样本协方差: 样本外 GMVP 月波动 5.064%
单指数模型: 样本外 GMVP 月波动 4.227%
PCA 4 因子: 样本外 GMVP 月波动 4.167%

• PCA：第一特征值 10.14（解释 34%），对应特征向量全部同号——市场；第二、三特征向量在不同行业上符号相反——行业对比，与例 9.1 的结构相同。真实共同因子有 5 个，但后两个特征值（1.03、0.98）已经混入噪声，碎石图的肘部在 4 附近。
• LR 检验：\(m=2,3\) 被拒绝，\(m=4\) 在 5% 水平被接受（p=0.15）。较弱的规模因子（波动 2%，特质波动 4–9%）在 \(T=240\) 时难以与噪声区分，这与原书例 9.4 中"统计上只能识别 3 个因子"类似。
• Varimax：因子 2、3 分别清晰对应行业 2、3（符号任意）；因子 1、4 是行业 1 与规模的混合。旋转能改善解释，但不保证得到"真实"因子——这正是 9.5.1 节不唯一性的体现。
• Bai–Ng：\(k=100\)、\(T=60\) 时，两个准则都选出 5 个因子，与真实秩一致。
• 样本外 GMVP：用 60 个月估计 30×30 协方差，样本协方差的 GMVP 样本外月波动 5.06%，单指数模型 4.23%，PCA 4 因子 4.17%。结构化的因子协方差比样本协方差降低了约 17% 的实现风险，这就是因子模型在组合优化中的价值。

本章小结

因子模型假设收益由少数共同因子和互不相关的特质成分构成，协方差简化为 \(\boldsymbol\beta\boldsymbol\Sigma_f\boldsymbol\beta'+\boldsymbol D\)。三类模型的区别在于什么是已知的：宏观因子模型已知因子、时间序列回归求载荷，解释力通常很弱；BARRA 已知载荷（行业、风格暴露）、逐期截面 WLS 求因子收益，GLS 权重就是纯因子组合；Fama–French 用排序对冲组合得到因子收益再求载荷；统计因子模型两者都未知，用 PCA 或极大似然估计，可用 LR 检验或碎石图定因子数，用 varimax 旋转辅助解释，\(k>T\) 时用 APCA 与 Bai–Ng 准则。实证上，股票的第一主成分是市场、第二主成分是行业对比，债券的前两个主成分是水平与斜率；单因子模型遗漏行业相关，行业哑变量模型可能高估行业内相关。在量化实务中，因子模型是风险模型、因子研究、组合优化和统计套利的共同基础。

练习

基础

1. \(k=100\)、\(m=10\) 的因子模型有多少个参数？与样本协方差矩阵相比少了多少？ 提示：\(100\times10+55+100=1155\)，样本协方差 5050 个。
2. 证明单因子 BARRA 模型中，\(\min\frac12\boldsymbol\omega'\boldsymbol D\boldsymbol\omega\) s.t. \(\boldsymbol\omega'\boldsymbol\beta=1\) 的解为 \(\boldsymbol\omega=\boldsymbol D^{-1}\boldsymbol\beta/(\boldsymbol\beta'\boldsymbol D^{-1}\boldsymbol\beta)\)。
3. 证明行业哑变量模型中，OLS 因子估计等于行业内等权平均收益。 提示：\(\boldsymbol\beta'\boldsymbol\beta=\mathrm{diag}(n_1,\dots,n_m)\)。
4. 对 \(\boldsymbol\Sigma=\begin{bmatrix}2&1\\1&2\end{bmatrix}\) 做 PCA，求两个主成分及其解释比例。 提示：特征值 3、1，特征向量 \((1,1)'/\sqrt2\)、\((1,-1)'/\sqrt2\)，比例 75%、25%。
5. 为什么宏观因子模型要用 VAR 残差（意外）而不是宏观变量本身？
6. 解释"Heywood 情形"，以及它在例 9.4 中提示了什么。

进阶

7. 证明主成分法中近似误差矩阵的元素平方和不超过 \(\sum_{i>m}\hat\lambda_i^2\)。 提示：\(\hat{\boldsymbol\Sigma}-\hat{\boldsymbol\beta}\hat{\boldsymbol\beta}'=\sum_{i>m}\hat\lambda_i\hat{\boldsymbol e}_i\hat{\boldsymbol e}_i'\)，其 Frobenius 范数平方为 \(\sum_{i>m}\hat\lambda_i^2\)；再减去对角阵 \(\hat{\boldsymbol D}\) 只会把对角元置零。
8. 修改示例二：加入"市场"因子（所有股票暴露为 1）与 3 个行业哑变量，观察 \(\boldsymbol\beta'\boldsymbol\beta\) 奇异；再施加"行业因子收益按股票数加权和为 0"的约束求解，比较结果。
9. 修改示例三：把估计窗口从 60 个月改为 36、120 个月，比较三种协方差估计的样本外 GMVP 波动，并加入 Ledoit–Wolf 收缩（sklearn.covariance.LedoitWolf）作对照。
10. 实现一个简化的 PCA 统计套利：对 \(k=100\) 只股票的日收益滚动提取 5 个主成分，计算每只股票的残差累积序列及其 z 分数，按 \(|z|>1.5\) 开仓、\(|z|<0.5\) 平仓，评估模拟数据上的表现。

原书推荐习题：9.4（BARRA 行业因子模型完整实现，风险模型的核心练习）、9.3（市场模型 + GMVP 比较协方差估计）、9.2(c) 与 9.6（统计因子分析、LR 检验、主成分法与 ML 法）、9.5（PCA 与碎石图）、9.7（用 VAR 残差构造宏观意外因子）。

原书对照

（原书印刷页码约等于 PDF 页码减 20。）