3. 最小二乘法:线性代数观点#

《九章算术》第八章：方程

今有卖牛二、羊五，以买十三豕，有馀钱一千。卖牛三、豕三，以买九羊，钱适足。卖羊六、豕八，以买五牛，钱不足六百。问牛、羊、豕价各几何?

\[\begin{split} \left[\begin{array}{r} 2 & 5 & -13\\ 3 & -9 & 3\\ -5 & 6 & 8 \end{array} \right] \left[\begin{array}{c} 牛\\ 羊\\ 猪 \end{array} \right] = \left[\begin{array}{r} 1000\\ 0\\ -600 \end{array} \right] \end{split}\]
A = matrix(c(2, 3, -5, 5, -9, 6, -13, 3, 8), nrow = 3)
b = c(1000, 0, -600)
solve(A,b)

最小二乘法 (OLS) 是计量经济学中最基本的估计方法，它简单透明，易于理解。了解最小二乘法，有助于我们研究更复杂的线性估计方法。此外，许多非线性估计量在真实值附近与线性估计量的行为是类似的。在本讲义中，我们将从线性代数的运算讲起，学习一系列的知识。

套用Leopold Kronecker的名言“上帝创造了整数，其他都是人的作品”，我想说“高斯创造了最小二乘法，其他都是应用研究者的作品”。在科学界，最小二乘法的受欢迎程度远远超出了人们的想象。然而，最小二乘法是一种纯统计学的方法，或者说是一种监督式机器学习的方法，因此它揭示的是相关关系，而非因果关系。相反，经济理论假设因果关系的存在，然后我们收集数据来检验理论或量化效果。

数学标记: \(y_{i}\) 是标量。\(x_{i}=\left(x_{i1},\ldots,x_{iK}\right)'\) 是 \(K\times1\) 的向量。\(Y=\left(y_{1},\ldots,y_{n}\right)'\) 是 \( n\times1\) 的向量。

\[\begin{split} X=\left[\begin{array}{c} x_{1}'\\ x_{2}'\\ \vdots\\ x_{n}' \end{array}\right]=\left[\begin{array}{cccc} x_{11} & x_{12} & \cdots & x_{1K}\\ x_{21} & x_{22} & \cdots & x_{2K}\\ \vdots & \vdots & \ddots & \vdots\\ x_{n1} & x_{22} & \cdots & x_{nK} \end{array}\right] \end{split}\]

是 \(n\times K\) 的矩阵。\(I_{n}\) 是 \(n\times n\) 的单位矩阵。

3.1. 最小二乘法#

给定 \(n\) 个观测值 \(\left(y_{i},x_{i}\right)_{i=1}^{n}\) ， 我们想利用 \(X\) 的某个线性组合 \(Xb\) 尽量逼近 \(Y\)。说到逼近，我们需要定义两个向量 \(Y\) 和 \(Xb\) 之间的距离。欧几里得范数(Euclidean norm) \(\Vert Y - Xb \Vert = \sqrt{\sum_{i=1}^{n}\left(y_{i}-x_{i}'b\right)^{2}}\) 是常用的测量向量长度的指标。 定义一个目标函数

\[ Q\left(b\right)=\frac{1}{2n}\sum_{i=1}^{n}\left(y_{i}-x_{i}'b\right)^{2}=\frac{1}{2n}\left(Y-Xb\right)'\left(Y-Xb\right)=\frac{1}{2n}\left\Vert Y-Xb\right\Vert ^{2}. \]

残差平方和 (sum of squared residuals) \(\sum_{i=1}^{n}\left(y_{i}-x_{i}'b\right)^{2}\) 是 \(\Vert Y - Xb \Vert\) 的平方。 因为平方 \((\cdot)^2\) 在 \(\mathbb{R}^+\) 上严格单调递增，所以用 \(b\) 来最小化 \(Q(b)\) 等价于最小化欧氏距离。系数 \(\frac{1}{2n}\) 与 \(b\) 无关，不影响该最优化问题的解。

最小化 \(Q\left(b\right)\) 的一阶必要条件是

\[\begin{split} \frac{\partial}{\partial b}Q\left(b\right)=\left[\begin{array}{c} \partial Q\left(b\right)/\partial b_{1}\\ \partial Q\left(b\right)/\partial b_{2}\\ \vdots\\ \partial Q\left(b\right)/\partial b_{K} \end{array}\right]=-\frac{1}{n}X'\left(Y-Xb\right)=0. \end{split}\]

对上式移项可得

\[ (X'X)b = X'Y. \]

在最小二乘法的推导中， 我们假定 \(X\) 中的 \(K\) 列是线性独立的， 也就是说不存在 \(K\times1\) 的向量 \(b\ (b\neq0_{K})\) 使得 \(Xb=0_{n}\) 。同时，这也意味着 \(n\geq K\) ，并且 \(X'X/n\) 可逆。 我们记该最优解

\[ \widehat{\beta}=\left(X'X\right)^{-1}X'y. \]

另外，二阶条件

\[\begin{split} \frac{\partial^{2}}{\partial b\partial b'}Q\left(b\right)=\left[\begin{array}{cccc} \frac{\partial^{2}}{\partial b_{1}^{2}}Q\left(b\right) & \frac{\partial^{2}}{\partial b_{2}\partial b_{2}}Q\left(b\right) & \cdots & \frac{\partial^{2}}{\partial b_{K}\partial b_{1}}Q\left(b\right)\\ \frac{\partial^{2}}{\partial b_{1}\partial b_{2}}Q\left(b\right) & \frac{\partial^{2}}{\partial b_{2}^{2}}Q\left(b\right) & \cdots & \frac{\partial^{2}}{\partial b_{K}\partial b_{2}}Q\left(b\right)\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial^{2}}{\partial b_{1}\partial b_{K}}Q\left(b\right) & \frac{\partial^{2}}{\partial b_{2}\partial b_{K}}Q\left(b\right) & \cdots & \frac{\partial^{2}}{\partial b_{K}^{2}}Q\left(b\right) \end{array}\right]=\frac{1}{n}X'X \end{split}\]

表明 \(Q\left(b\right)\) 是关于 \(b\) 的凸函数，因为 \(X'X/n\) 是半正定矩阵。(如果 \(X'X/n\) 是正定矩阵，那么 \(Q\left(b\right)\) 是关于 \(b\) 的严格凸函数。)

评注 3.1

如果某些回归因子满足完全共线性 (perfectly collinear)，就违反了 \(X'X/n\) 可逆的线性独立的假设。 例如，使用虚拟变量来表示分类变量并将所有这些类别放入回归模型中时，通常计量经济学软件会自动检测并指出完全共线性问题。然而，软件难以侦测不完全共线性 (nearly collinear)，即 \(X'X/n\) 的最小特征值接近于0，而不等于0。我们将在渐进理论一章中讨论不完全共线性的后果。

3.1.1. 线性代数性质#

下面列举一些与最小二乘法估计量有关的定义与性质。

  • 拟合值 (Fitted value): \(\widehat{Y}=X\widehat{\beta}\).

  • 投影矩阵 (Projection matrix): \(P_{X}=X\left(X'X\right)^{-1}X\) ：残差生成矩阵 (Residual maker matrix)： \(M_{X}=I_{n}-P_{X}\).

  • \(P_{X}X=X\); \(X'P_{X}=X'\).

  • \(M_{X}X=0_{n\times K}\); \(X'M_{X}=0_{K\times n}\).

  • \(P_{X}M_{X}=M_{X}P_{X}=0_{n\times n}\).

  • \(P_{X}\) 与 \(M_{X}\) 都是 幂等矩阵 。

  • 如果 \(AA=A\) ， 那么 \(A\) 被称作 幂等矩阵 (idempotent matrix)。幂等矩阵的特征值只能是1或者0。

  • \(\mathrm{rank}\left(P_{X}\right)=K\) ， \(\mathrm{rank}\left(M_{X}\right)=n-K\).

  • 残差: \(\widehat{e}=Y-\widehat{Y}=Y-X\widehat{\beta}=Y-X(X'X)^{-1}X'Y=(I_{n}-P_{X})Y=M_{X}Y.\)

  • \(X'\widehat{e}=0_{K}\).

  • 若 \(x_{i}\) 中有一个为常数，则 \(\sum_{i=1}^{n}\widehat{e}_{i}=0\) 。

(因为 \(X'\widehat{e}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1\\ \heartsuit & \heartsuit & \cdots & \heartsuit\\ \cdots & \cdots & \ddots & \vdots\\ \heartsuit & \heartsuit & \cdots & \heartsuit \end{array}\right]\left[\begin{array}{c} \widehat{e}_{1}\\ \widehat{e}_{2}\\ \vdots\\ \widehat{e}_{n} \end{array}\right]=\left[\begin{array}{c} 0\\ 0\\ \vdots\\ 0 \end{array}\right]\) ，第一行运算说明了 \(\sum_{i=1}^{n}\widehat{e}_{i}=0\)。”\(\heartsuit\)” 代表与我们计算无关的值。)

We will produce a graph here

consider replace this simulation example with a real data example

例子 3.1

下面是一个简单的数值模拟案例，我们以此来说明最小二乘法估计量的特性。给定 \(\left(x_{1i},x_{2i},x_{3i},e_{i}\right)^{\prime}\sim N\left(0_{4},I_{4}\right)\) ， 因变量 \(y_{i}\) 的生成式为

\[ y_{i}=0.5+2\cdot x_{1i}-1\cdot x_{2i}+e_{i} \]

在不知道 \(x_{3i}\) 是多余的情况下，我们将 \(y_{i}\) 对 \(\left(1,x_{1i},x_{2i},x_{3i}\right)\) 进行回归。

最后得到的参数估计值与真实值接近。当然，由于样本容量较小，答案不是十分准确。

library(magrittr); set.seed(2022-6-15)
n = 20 # sample size 
K = 4 # number of paramters
b0 = as.matrix( c(0.5, 2, -1, 0) ) # the true coefficient
X = cbind(1, matrix( rnorm(n * (K-1)), nrow = n ) ) # the regressor matrix 
e = rnorm(n) # the error term
Y = X %*% b0 + e # generate the dependent variable
bhat = solve(t(X) %*% X, t(X) %*% Y ) %>% as.vector() %>% print()

我们继续用上面的数值例子来验证一些代数性质。

ehat = Y - X %*% bhat 
as.vector( t(X) %*% ehat ) %>% print()
MX = diag(n) - X %*% solve( crossprod(X) ) %*% t(X)

data.frame(e = e, ehat = ehat, MXY = MX%*%Y ) %>% head()
cat("The mean of the residual is ", mean(ehat), ".\n")
cat("The mean of the true error term is", mean(e), ".")

注意到，上例中 \(e_i\) 是用来生成 \(y_i\) 的扰动项，它和 \(\widehat{e}_i\) 不是一回事。\(e_i\) 只是出现在我们的数值模拟当中，在实际回归中我们只观测到 \((y_i,x_i)\) ，而 \(e_i\) 无从观测。

3.1.2. 几何意义#

\(\mathcal{X}=\left\{ Xb:b\in\mathbb{R}^{K}\right\}\) 是一个由 \(X=\left[X_{\cdot1},\ldots,X_{\cdot K}\right]\) 中的 \(K\) 列生成的线性空间; 如果这些列是线性独立的，那么 \(X\) 就是 \(K\) 维的。最小二乘法估计量 \(\widehat \beta\) 使得 \(\left\Vert Y-Xb\right\Vert\) 最小化。换言之， \(X\widehat{\beta}\) 是 \(\mathcal{X}\) 中与 \(Y\) 最接近的那个点。

等式 \(Y=X\widehat{\beta}+\widehat{e}\) 将 \(Y\) 分解为两个垂直的分量， \(X\widehat{\beta}\) 与 \(\widehat{e}\) ， 因为\(\left\langle X\widehat{\beta},\widehat{e}\right\rangle =\widehat{\beta}'X'\widehat{e}=0_{K}^{\prime}\) ， 其中\(\left\langle \cdot,\cdot\right\rangle\) 是向量的内积运算。那么， \(X\widehat{\beta}\) 就是 \(Y\) 在 \(\mathcal{X}\) 上的 投影 (projection)， \(\widehat{e}\) 则是对应的 投影残差 (projection residuals)，根据勾股定理，自然有

\[ \left\Vert Y\right\Vert ^{2}=\Vert X\widehat{\beta}\Vert^{2}+ \Vert \widehat{e}\Vert^{2}. \]

3.2. 子向量 (Subvector)#

有时候，我们并不关心整个向量 \(\widehat{\beta}\) 的数值，而是关心某个变量的系数，比如一个子向量 \(\widehat{\beta}_1\) 的数值。Frish-Waugh-Lovell (FWL) 定理 揭示最小二乘法子向量的代数性质。要得到FWL定理，我们需要分块矩阵求逆的知识。对称实数正定矩阵 \(A=\begin{pmatrix}A_{11} & A_{12}\\ A_{12}' & A_{22} \end{pmatrix}\) 逆矩阵可以写成

\[\begin{split} A^{-1}=\begin{pmatrix}\left(A_{11}-A_{12}A_{22}^{-1}A_{12}'\right)^{-1} & -\left(A_{11}-A_{12}A_{22}^{-1}A_{12}'\right)^{-1}A_{12}A_{22}^{-1}\\ -A_{22}^{-1}A_{12}'\left(A_{11}-A_{12}A_{22}^{-1}A_{12}'\right)^{-1} & \left(A_{22}-A_{12}'A_{11}^{-1}A_{12}\right)^{-1} \end{pmatrix}. \end{split}\]

将此性质运用至最小二乘法估计。记 \(X=(X_1, X_2)\) ，则

(3.1)#\[\begin{split} \begin{aligned} \begin{pmatrix}\widehat{\beta}_{1}\\ \widehat{\beta}_{2} \end{pmatrix} & =\widehat{\beta}=(X'X)^{-1}X'Y\\ & =\left(\begin{pmatrix}X_{1}'\\ X_{2}' \end{pmatrix}\begin{pmatrix}X_{1} & X_{2}\end{pmatrix}\right)^{-1}\begin{pmatrix}X_{1}'Y\\ X_{2}'Y \end{pmatrix}\\ & =\begin{pmatrix}X_{1}'X_{1} & X_{1}'X_{2}\\ X_{2}'X_{1} & X_{2}'X_{2} \end{pmatrix}^{-1}\begin{pmatrix}X_{1}'Y\\ X_{2}'Y \end{pmatrix}\\ & =\begin{pmatrix}\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1} & -\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'X_{2}\left(X_{2}'X_{2}\right)^{-1}\\ \heartsuit & \heartsuit \end{pmatrix}\begin{pmatrix}X_{1}'Y\\ X_{2}'Y \end{pmatrix}. \end{aligned} \end{split}\]

\(\widehat{\beta}\) 的子向量 \(\widehat{\beta}_{1}\) 可写成

\[\begin{split} \begin{aligned} \widehat{\beta}_{1} & =\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'Y-\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'X_{2}\left(X_{2}'X_{2}\right)^{-1}X_{2}'Y\\ & =\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'Y-\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'P_{X_{2}}Y\\ & =\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}\left(X_{1}'Y-X_{1}'P_{X_{2}}Y\right)\\ & =\left(X_{1}'M_{X_{2}}'X_{1}\right)^{-1}X_{1}'M_{X_{2}}Y.\end{aligned} \end{split}\]

以上 \(\widehat{\beta}_{1}\) 的表达式就是 FWL 定理。该表达式意味着我们可以用三步方法得到 \(\widehat{\beta}_{1}\) :

1.将 \(Y\) 对 \(X_{2}\) 做回归，得到残差 \(\tilde{Y}\) ；

2.将 \(X_{1}\) 对 \(X_{2}\) 做回归，得到残差 \(\tilde{X}_{1}\) ；

3.将 \(\tilde{Y}\) 对 \(\tilde{X}_{1}\) 做回归，得到最小二乘法估计值 \(\widehat{\beta}_{1}\) 。

同样的方法也适用于总体线性投影 (population linear projection)，参考 Hansen (2020) [E] Chapter 2.22-23。

# verify FWL theorem

X1 = X[,1:2];X2 = X[,3:4]
PX2 = X2 %*% solve( t(X2) %*% X2) %*% t(X2) 
MX2 = diag(rep(1,n)) - PX2

bhat1 <- (solve(t(X1)%*% MX2 %*% X1, t(X1) %*% MX2 %*% Y )) %>%
 as.vector() %>% print()

ehat1 = MX2 %*% Y - MX2 %*% X1 %*% bhat1 
data.frame(ehat = ehat, ehat1 = ehat1) %>% head() %>% print()

3.3. 适配度 (Goodness of Fit)#

考虑一个带截距项的最小二乘法估计

(3.2)#\[ Y=\widehat{Y}+\widehat{e}=(X_{1}\widehat{\beta}_{1}+\widehat{\beta}_{2})+\widehat{e} \]

应用FWL定理，令 \(X_{2}=\iota\) ，其中希腊字母 \(\iota\) (iota)是一个所有位置都为1的 \(n\times1\) 向量。那么\(M_{X_{2}}=M_{\iota}=I_{n}-\frac{1}{n}\iota\iota'\) 。其中， \(M_{\iota}z=z-\bar{z}\) 。也就是说，我们在原本的向量 \(z\) 中减去其均值 \(\bar{z}=\frac{1}{n}\sum_{i=1}^{n}z_{i}\) 。

那么，上述的三步方法变为

1.将 \(Y\) 对 \(\iota\) 做回归，得到残差 \(M_{\iota}Y\) ；

2.将 \(X_{1}\) 对 \(\iota\) 做回归，得到残差 \(M_{\iota}X_{1}\) ；

3.将 \(M_{\iota}Y\) 对 \(M_{\iota}X_{1}\) 做回归，根据 (3.1) 得到最小二乘法估计值 \(\widehat{\beta}_{1}\) 。

我们将最后一步进行分解

(3.3)#\[ M_{\iota}Y=M_{\iota}X_{1}\widehat{\beta}_{1}+\tilde{e}, \]

应用勾股定理得到

\[ \left\Vert M_{\iota}Y\right\Vert ^{2}=\Vert M_{\iota}X_{1}\widehat{\beta}_{1}\Vert^{2}+ \Vert \widehat{e} \Vert^{2}. \]

练习 3.1

证明：(3.2) 式中的 \(\widehat{e}\) 与 (3.3) 式中的 \(\tilde{e}\) 相等。

在线性回归中， 决定系数 ( \(R^2\) ) 是一个衡量适配度的指标。样本内决定系数(in-sample \(R^2\))可写作

\[ R^{2}=\frac{\Vert M_{\iota}X_{1}\widehat{\beta}_{1}\Vert^{2}}{\left\Vert M_{\iota}Y\right\Vert ^{2}}=1-\frac{\left\Vert \tilde{e}\right\Vert ^{2}}{\left\Vert M_{\iota}Y\right\Vert ^{2}}. \]

当归回包含截距项时， \(R^2\) 才有定义。

练习 3.2

在 (3.2) 式的分解下，证明 \(R^{2}\) 即为 \(\widehat{Y}\) 的样本方差与 \(Y\) 的样本方差之比：

\[ R^{2}=\frac{\widehat{Y}'M_{\iota}\widehat{Y}}{Y'M_{\iota}Y}=\frac{\sum_{i=1}^{n}\left(\widehat{y_{i}}-\overline{y}\right)^{2}}{\sum_{i=1}^{n}\left(y_{i}-\overline{y}\right)^{2}}. \]

决定系数 \(R^2\) 的大小在不同的实际问题当中差别很大。\(R^2\) 大于90%的情况在带有滞后效应的宏观模型中并不罕见。然而，在截面数据回归 (cross sectional regressions) 中， \(R^2\) 的值经常是低于20%。

练习 3.3

考虑一个较“短”的回归 “将 \(y_{i}\) 对 \(x_{1i}\) 做回归” 与一个较“长”的回归 “将 \(y_{i}\) 对 \(\left(x_{1i},x_{2i}\right)\) 做回归”: 给定相同的数据集 \(\left(Y,X_{1},X_{2}\right)\) ， 试说明“短”回归的 \(R^2\) 不比“长”回归的 \(R^2\) 大。也就是说，通过增加更多的回归变我们总能使得 \(R^2\) 不会变小。)

在传统意义上，回归因子的数量 \(K\) 总是远小于样本量 \(n\) 。然而，在大数据时代，回归变量的数量 \(K\) 有可能大于样本量 \(n\) 。

练习 3.4

证明: 当 \(K\geq n\) 时， \(R^{2}=1\)。注意，若 \(K>n\) ， 则矩阵 \(X'X\) 是缺秩 (rank deficient) 的。在这种情况下，我们可以将最小二乘法的定义扩展， \(\hat \beta\) 仍然是使得 \(\left\Vert Y-Xb\right\Vert ^{2}\) 最小化的参数值，但这个值不是唯一的。

以下是本练习的数值实例。

n = 5; K = 6; 
Y = rnorm(n)
X = matrix( rnorm(n*K), n)
summary( lm(Y~X) )

对于一个新的数据集 \(\left(Y^{\mathrm{new}},X^{\mathrm{new}}\right)\) ， 样本外决定系数(OOS \(R^2\)) 可写作

\[ OOS\ R^{2}=\frac{\widehat{\beta}^{\prime}X^{\mathrm{new}\prime}M_{\iota}X^{\mathrm{new}}\widehat{\beta}}{Y^{\mathrm{new}\prime}M_{\iota}Y^{\mathrm{new}}}. \]

从原始数据集中得到参数估计值以后， \(OOS\ R^{2}\) 衡量这一参数估计值在新数据集上的适配度。在金融市场的短期预测模型中，如果某策略能够系统性地获得2%的 \(OOS\ R^{2}\) ，那么这就是一个很好的赚钱策略。

3.4. 总结#

本章完全是线性代数操作，无关数据背后的概率模型。只要 \(X'X/n\) 可逆，以上线性代数性质对任何有限样本 \(n\) 都成立。

历史趣闻

Carl Friedrich Gauss (1777–1855) 宣称他在1795年就发明了最小二乘法，那时他用三个数据点预测了1801年矮行星谷神星的位置。虽然 Gauss 在1809年才发表了相关文章，而 Adrien-Marie Legendre (1752–1833) 在1805年就公开提出了此方法，但今天人们仍倾向于将最小二乘法的发明归功于 Gauss。因为像 Gauss 这样的数学巨人没有必要通过撒谎来窃取 Legendre 的成果。