评注章节：Chapter 4. Least Squares Regression (HANSEN 2021 ,Chpt4, Sec4.21, pg128)。
评注要点：
群组分析中可能都会面临一个软件计算问题，也就是需要在群组之间进行循环（loop）计算。
理论过程
群组稳健协方差（Arellano法）
现在考虑 \(\widehat{\beta}\) 的协方差矩阵。
令 \(\Sigma_{g}=\mathbb{E}\left[\boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime} \mid \boldsymbol{X}_{g}\right]\) 表示 \(g^{t h}\) 群组内误差的 \(n_{g} \times n_{g}\) 条件协方差矩阵。
由于观测结果在群组之间是独立的，
\[
\begin{aligned}
\operatorname{var}\left[\left(\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \boldsymbol{e}_{g}\right) \mid \boldsymbol{X}\right]
&=\sum_{g=1}^{G} \operatorname{var}\left[\boldsymbol{X}_{g}^{\prime} \boldsymbol{e}_{g} \mid \boldsymbol{X}_{g}\right] \\
&=\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \mathbb{E}\left[\boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime} \mid \boldsymbol{X}_{g}\right] \boldsymbol{X}_{g} \\
&=\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \Sigma_{g} \boldsymbol{X}_{g} \\
& \stackrel{\text { def }}{=} \Omega_{n}
\end{aligned}
\tag{1}\]
进一步考察系数方差的计算公式
\[
\boldsymbol{V}_{\widehat{\beta}}=\operatorname{var}[\widehat{\beta} \mid \boldsymbol{X}]=\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \Omega_{n}\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} .
\tag{2}\]
由于群组内观测值之间的相关性，这与独立情况下的公式不同。差异的大小取决于群组内观测值之间的相关程度以及群组内观测值的数量。
例子 1 (简单情形计算) 要看到这一点，我们做一个简单情形的计算比较。
这一简单情形有如下假设：
所有群组对于 \(i \neq \ell\) 都有相同数量的观测值 \(n_{g}=N\)
\(\mathbb{E}\left[e_{i g}^{2} \mid \boldsymbol{X}_{g}\right]=\sigma^{2}\)，也即各组之间随机干扰项是同方差的
\(\mathbb{E}\left[e_{i g} e_{\ell g} \mid \boldsymbol{X}_{g}\right]=\sigma^{2} \rho\)，也即组内随机干扰项之间是相关的
并且回归元 \(X_{i g}\) 在群组内没有变化。
在这种情况下，OLS 估计量的精确方差等于 （经过一些计算）
\[
\boldsymbol{V}_{\widehat{\beta}}=\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \sigma^{2}(1+\rho(N-1))
\]
如果 \(\rho>0\) ，则正确方差大概是传统公式的 \(\rho N\) 倍。在肯尼亚学校的示例中，平均群组大小为 48 。如果 \(\rho=0.25\) 这意味着精确方差超出传统公式约 12 倍。在这种情况下，正确的标准误差（方差的平方根）大约是传统公式的三倍。这是一个巨大的差异，不应被忽视。
Arellano (1987) 提出了一种群组稳健协方差矩阵估计量，它是 White 估计量的扩展。回想一下，White 协方差估计量的见解是平方误差 \(e_{i}^{2}\) 对于 \(\mathbb{E}\left[e_{i}^{2} \mid X_{i}\right]=\sigma_{i}^{2}\) 是无偏的。类似地，由于群组依赖性，矩阵 \(\boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime}\) 对于 \(\mathbb{E}\left[\boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime} \mid \boldsymbol{X}_{g}\right]=\Sigma_{g}\) 是无偏的。这意味着 (4.51) 的无偏估计量是 \(\widetilde{\Omega}_{n}=\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime} \boldsymbol{X}_{g}\)。这是不可行的，但我们可以用 OLS 残差替换未知误差以获得 Arellano 估计量：
\[
\begin{aligned}
\widehat{\Omega}_{n} &=\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \widehat{\boldsymbol{e}}_{g} \widehat{\boldsymbol{e}}_{g}^{\prime} \boldsymbol{X}_{g}
&& \text{(a)}\\
&=\sum_{g=1}^{G} \sum_{i=1}^{n_{g}} \sum_{\ell=1}^{n_{g}} X_{i g} X_{\ell g}^{\prime} \widehat{e}_{i g} \widehat{e}_{\ell g}
&& \text{(b)} \\
&=\sum_{g=1}^{G}\left(\sum_{i=1}^{n_{g}} X_{i g} \widehat{e}_{i g}\right)\left(\sum_{\ell=1}^{n_{g}} X_{\ell g} \widehat{e}_{\ell g}\right)^{\prime}
&& \text{(c)}
\end{aligned}
\tag{3}\]
式 3 中的三个表达式给出了三个可用于计算 \(\widehat{\Omega}_{n}\) 的等效公式。最终表达式(c)根据群组和 \(\sum_{\ell=1}^{n_{g}} X_{\ell g} \widehat{e}_{\ell g}\) 写入 \(\widehat{\Omega}_{n}\)，这是我们的示例 \(\mathrm{R}\) 和下面所示的 MATLAB 代码的基础。
给定表达式 式 1 - 式 2，群组协方差矩阵的一个自然的估计量采用以下形式
\[
\widehat{\boldsymbol{V}}_{\widehat{\beta}}=a_{n}\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \widehat{\Omega}_{n}\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1}
\tag{4}\]
其中 \(a_{n}\) 是可能的有限样本调整。 Stata群组命令使用的调整因子为：
\[
a_{n}=\left(\frac{n-1}{n-k}\right)\left(\frac{G}{G-1}\right)
\]
因子 \(G /(G-1)\) 是由 Chris Hansen (2007) 在同等大小群组的背景下导出的，目的是在群组 \(G\) 数量较小时提高性能。调整因子 \((n-1) /(n-k)\) 是专门的一般化表达，它嵌套了HC1计算公式中使用的调整，因为 \(G=n\) 意味着简化 \(a_{n}=n /(n-k)\)。
群组稳健协方差（leave-one-out法）
可以使用群组级（cluster-level）预测误差构建替代的群组稳健协方差矩阵估计量，例如 \(\widetilde{\boldsymbol{e}}_{g}=\boldsymbol{Y}_{g}-\boldsymbol{X}_{g} \widehat{\beta}_{(-g)}\)，其中 \(\widehat{\beta}_{(-g)}\) 是省略群组 \(g\) 的最小二乘估计量。
如第 3.20 节所示，我们可以证明
\[
\widetilde{\boldsymbol{e}}_{g}=\left(\boldsymbol{I}_{n_{g}}-\boldsymbol{X}_{g}\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \boldsymbol{X}_{g}^{\prime}\right)^{-1} \widehat{\boldsymbol{e}}_{g}
\tag{5}\]
和
\[
\widehat{\beta}_{(-g)}=\widehat{\beta}-\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \boldsymbol{X}_{g}^{\prime} \widetilde{\boldsymbol{e}}_{g} .
\tag{6}\]
然后我们就有了稳健的协方差矩阵估计量
\[
\widehat{\boldsymbol{V}}_{\widehat{\beta}}^{\mathrm{CR} 3}=\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1}\left(\sum_{g=1}^{G} \boldsymbol{X}_{g}^{\prime} \widetilde{\boldsymbol{e}}_{g} \widetilde{\boldsymbol{e}}_{g}^{\prime} \boldsymbol{X}_{g}\right)\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1}
\tag{7}\]
标签“CR”是指“群组稳健性”，“CR3”是指 HC3 估计量的类似公式。
与异方差稳健情况类似，你可以证明 CR3 是 \(\boldsymbol{V}_{\widehat{\beta}}\) 的保守估计量，因为 \(\widehat{\boldsymbol{V}}_{\widehat{\beta}}^{\mathrm{CR} 3}\) 的条件期望超过 \(\boldsymbol{V}_{\widehat{\beta}}\)。然而，这种协方差矩阵估计量实现起来比较麻烦，因为群组级预测误差 式 5 无法通过简单的线性运算来计算，并且需要跨群组的循环来计算。
案例实证
案例背景
为了在肯尼亚学校教育示例的背景下进行说明，我们在学校级别的跟踪虚拟模型上展示了学生测试成绩的回归，并显示了两个标准误差。
\[
\begin{alignedat}{999}
TestScore_{i g}=&-
0.071 &&+ 0.138 \text { Tracking }_g +e_{i g}\\
&(0.019) &&(0.026) \\
& [0.054] &&[0.078]
\end{alignedat}
\tag{8}\]
其中：
我们可以看到，群组稳健标准误大约是传统稳健标准误的三倍。因此，系数的置信区间很大程度上受到选择的影响。
为了便于说明，我们在此列出了在 Stata 和R中生成具有群组标准误差的回归结果所需的命令。
Arellano计算法（R代码）
代码R文件路径为code-cluster-var-coef.R。
```{r}
#| label: code-ddk2011
#| code-line-numbers: true
# text source file: "chapt4.R"
# Load the data DDK2011 and create variables
data <- haven::read_dta(here("code/Econometrics-Programs/DDK2011.dta"))
y <- scale(as.matrix(data$totalscore)) # 标准化变换
n <- nrow(y)
x <- cbind(as.matrix(data$tracking),matrix(1,n,1))
schoolid <- as.matrix(data$schoolid)
k <- ncol(x)
xx <- t(x)%*%x
invx <- solve(xx)
beta <- solve(xx,t(x)%*%y)
xe <- x*rep(y-x%*%beta,times=k) # x矩阵的每列元素（n个）都乘以对应的残差元素（n个）
# Clustered robust standard error
xe_sum <- rowsum(xe, group = schoolid)
G <- nrow(xe_sum)
omega <- t(xe_sum)%*%xe_sum
scale <- G/(G-1)*(n-1)/(n-k)
V_clustered <- scale*invx%*%omega%*%invx
se_clustered <- sqrt(diag(V_clustered))
names(se_clustered) <- c("se_b2", "se_b1")
se_clustered
```
se_b2 se_b1
0.07723621 0.05439336
这里\(\boldsymbol{\widehat{\Omega}}_n\)计算实际上使用了 式 3 (c)的计算公式。
R代码解读
source(here::here("code/demo-five.R"))
kable(dt, align = "c", digits = c(0, 1,4))
表 1: 示例数据集
| 1 |
1.0 |
1.1669 |
| 1 |
1.5 |
-0.3617 |
| 1 |
2.0 |
1.2458 |
| 2 |
2.5 |
2.7125 |
| 2 |
3.0 |
2.2266 |
群组稳健性标准差计算代码如下：
# x矩阵的每列元素（n个）都乘以对应的残差元素（n个）
xe <- x*rep(e,times=k)
# same as above
#xe <- x*(e %*% matrix(1,1,k))
# Clustered robust standard error
id <- as.matrix(id)
xe_sum <- rowsum(xe, group = id)
G <- nrow(xe_sum)
omega <- t(xe_sum)%*%xe_sum
scale <- G/(G-1)*(n-1)/(n-k)
V_clustered <- scale*invx%*%omega%*%invx
se_clustered <- sqrt(diag(V_clustered))
names(se_clustered) <- c("se_b2", "se_b1")
se_clustered
se_b2 se_b1
0.17313569 0.07608076
xe <- x*rep(e,times=k)，等价于e %*% matrix(1,1,k)。代码中x为\(n \times 2\)的矩阵（matrix），而rep(e, times =k)则是一个长度为\(n \times 2\)的向量（vector），二者相乘的结果，实际上是X矩阵的每列元素（n个）都乘以对应的残差元素（n个）。
\[
\begin{aligned}
\begin{bmatrix}
X_{11} &1 \\
X_{12} & 1 \\
\vdots & \vdots\\
X_{1n} & 1
\end{bmatrix}
\odot
\begin{bmatrix}
e_{1} &e_1 \\
e_{2} & e_2 \\
\vdots & \vdots\\
e_{n} & e_n
\end{bmatrix}
\end{aligned}
\]
rep(e,times=k)
[1] 0.8075928 -1.2403359 -0.1522368 0.7951100 -0.2101301 0.8075928
[7] -1.2403359 -0.1522368 0.7951100 -0.2101301
等价于e %*% matrix(1,1,k)
[,1] [,2]
[1,] 0.8075928 0.8075928
[2,] -1.2403359 -1.2403359
[3,] -0.1522368 -0.1522368
[4,] 0.7951100 0.7951100
[5,] -0.2101301 -0.2101301
xe_sum <- rowsum(xe, group = schoolid)可以基于群组标准（这里是schoolid）进行群组内加总。
xe <- x*(e %*% matrix(1,1,k))
[,1] [,2]
[1,] 0.8075928 0.8075928
[2,] -1.8605038 -1.2403359
[3,] -0.3044736 -0.1522368
[4,] 1.9877749 0.7951100
[5,] -0.6303903 -0.2101301
xe_sum <- rowsum(xe, group = id)
[,1] [,2]
1 -1.357385 -0.5849799
2 1.357385 0.5849799
leave-one-out计算法
需要使用群组循环的方法。相关计算方法类似于 节删一法回归
参考
HANSEN, BRUCE. 2021. 《Econometrics》.