因果面板数据模型的矩阵补全方法

发布时间：2022-04-24 阅读 1598

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：王文茂 (中山大学)
邮箱：wangwm8@mail2.sysu.edu.cn

编者按：本文主要摘译自下文，特此致谢！
Source: Athey S, Bayati M, Doudchenko N, et al. Matrix completion methods for causal panel data models[J]. Journal of the American Statistical Association, 2021, 116(536): 1716-1730. -PDF-

1. 引言

相关性并不意味着因果关系。在进行因果推断时，我们总是希望能够找到一个合理的反事实对照组。已有文献在估计面板数据的平均处理效应 (ATE) 时，为了得到该对照组，采取了非混淆性和合成控制两种方法。其中，

非混淆性可理解为，使用与前一时期观察结果值相似的对照单元观察结果，来推算处理单元丢失的潜在对照结果；
合成控制可理解为，使用对照单元的加权平均结果，来计算处理单元的潜在对照结果，并选择权重，以便对照单元的加权滞后结果与处理单元的滞后结果相匹配。

本文在前人的基础上提出一种新的方法来解决这一问题。借鉴关于因子模型和交互固定效应的计量经济学文献，以及关于矩阵补全的计算机科学和统计学文献，作者提出核范数矩阵补全估计法。在基于真实数据的模拟中，证明该方法的性能优于现有的方法。

本文的贡献主要有以下三方面：

一是推广了矩阵补全文献中的结果。修改了矩阵补全和因子模型文献中的估计量 (提出核范数矩阵补全估计法)，以考虑非正则化的单位和时间固定效应；
二是证明了非混淆性和合成控制，以及核范数最小化方法都可以看作是基于矩阵分解的矩阵补全方法。它们都具有基于 Fröbenius 范数的目标函数。根据这个共同的目标函数，非混淆性和合成控制方法对矩阵分解中的因素施加了不同的限制。核范数最小化方法不施加任何限制，而是使用正则化来表示估计量；
三是将该方法应用于两个真实数据集。人为地指定某些单元和时间段缺失的结果，然后比较不同补全估计量的表现。核范数矩阵补全估计量与另外的方法相比，在一系列情况下表现更为良好。

2. 基本模型设置

在 $T$ 个时期内观察到 $N$ 个单位的面板数据中，每个时期每个单位都有两个潜在的结果： $Y_{i t}$ (0) 和 $Y_{i t} (1)$ 。 $W_{i t} = 1$ 表示该单元被处理，否则 $W_{i t} = 0$ 。对于每个单元和时期，我们观察结果为：

Y_{i t} = Y_{i t} (W_{i t})

文章集中于估计被处理组的平均效应 (ATE)：

τ = \sum_{(i, t) : W_{i t} = 1} [Y_{i t} (1) - Y_{i t} (0)] / \sum_{i, t} W_{i t}

要估计这样的平均处理效应，一种方法是计算缺失的潜在结果。由于已经观察到了 $Y_{i t} (1)$ 的所有相关值，我们只需要为 $W_{i t} = 1$ 的受试者计算其相对应 $Y (0)$ 矩阵的缺失值。这样就可以使 $Y_{i t} (1)$ 拥有与其相对应的对照项 $Y_{i t} (0)$ ，从而能够计算其效应。

如下所示， $Y$ 、 $W$ 均是 $N \times T$ 的矩阵。同时，将 $M$ 定义为对应 $Y$ 中 $W_{i t} = 1$ 对应的 $(i, t)$ 中的项集合，将 $O$ 定义为对应 $Y$ 中 $W_{i t} = 0$ 对应的 $(i, t)$ 中的项集合。文章主要就是找到关于在 $Y$ 中缺失值的统计问题 (即公式中的 $?$ )，一旦这些值被输入，我们就可以估计平均因果效应 $τ$ 。

Y = (\begin{array}{ccccc} Y_{11} & Y_{12} & ? & \dots & Y_{1 T} \\ ? & ? & Y_{23} & \dots & ? \\ Y_{31} & ? & Y_{33} & \dots & ? \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ Y_{N 1} & ? & Y_{N 3} & \dots & ? \end{array}),

W = (\begin{array}{ccccc} 0 & 0 & 1 & \dots & 0 \\ 1 & 1 & 0 & \dots & 1 \\ 0 & 1 & 0 & \dots & 1 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 1 & 0 & \dots & 1 \end{array}),

其中，

W_{i t} = {\begin{cases} 1 & if (i, t) \in M \\ 0 & if (i, t) \in O \end{cases}

3. 矩阵形式和分类

本节中将讨论矩阵 $Y$ 和 $W$ 的一些特殊形式并分类。矩阵补全文献主要关注 $W$ 是完全随机的情况，如上节所示的情形，并且 $Y$ 和 $W$ 的维度都很大。首先，考虑缺失数据的模式，即 $W$ 的分布不同于完全随机的分布。其次，考虑矩阵 $Y$ 的不同形状，其中维度 $N$ 和 $T$ 的相对大小可能非常不同。第三，考虑计量经济学文献中的一些具体分析，这些分析侧重于缺失数据模式和矩阵形状的特定组合。

3.1 数据缺失模式

数据缺失模式主要分为两种，块状结构和交错处理结构。其中，块状结构如下所示：

Y_{N \times T} = (\begin{array}{cccccc} ✓ & ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ✓ & ? & \dots & ? \\ ✓ & ✓ & ✓ & ? & \dots & ? \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ ✓ & ✓ & ✓ & ? & \dots & ? \end{array})

块状结构有两种特殊情况。在非混淆性 (Imbens 和 Rubin，2015) 下，估计平均处理效应的文献大多集中在 $T_{0} = T - 1$ 的情况下，因此仅有的处理单元在最后一个周期，称为单处理周期块结构 (Single-TreatedPeriod-Block Structure)。相反，合成控制文献 (Abadie 等，2010；Abadie，2019) 主要关注具有从时段 $T_{0} + 1$ 开始的多个时段处理的单个处理单元的情况，称为单处理单元块结构，如下所示：

Y = (\begin{array}{cccccc} ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ & ? \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ ✓ & ✓ & ✓ & \dots & ✓ & ? \\ ↑ \\ treated period \end{array})

Y = (\begin{array}{llllll} ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ ✓ & ✓ & ✓ & \dots & ✓ \\ ✓ & ✓ & ? & \dots & ? & \leftarrow treated unit \end{array})

对于仅仅缺少单个 $Y_{i t}$ 的情况，上面提到的两种方法均适用：

Y = (\begin{array}{cccccc} ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ ✓ & ✓ & ✓ & \dots & ✓ & ✓ \\ ✓ & ✓ & ✓ & \dots & ✓ & ? \end{array})

交错处理结构是指不同的单位第一次接触处理的时间可能不同 (Athey 和 Imbens，2018；Shaikh 和 Toulis，2019)，具体如下所示：

Y_{N \times T} = (\begin{array}{ccccccc} ✓ & ✓ & ✓ & ✓ & \dots & ✓ & (never adopter) \\ ✓ & ✓ & ✓ & ✓ & \dots & ? & (late adopter) \\ ✓ & ✓ & ? & ? & \dots & ? \\ ✓ & ✓ & ? & ? & \dots & ? & (medium adopter) \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ ✓ & ? & ? & ? & \dots & ? & (early adopter) \end{array})

3.2 瘦矩阵和胖矩阵

第二种分类涉及矩阵 $Y$ 的形状。由于 $N 、 T$ 的相对大小变化，矩阵的形状就会产生不同。例如，

当 $N ≫ T$ 时， $Y$ 是瘦矩阵：

\begin{aligned} Y = (\begin{array}{ccc} ? & ✓ & ? \\ ✓ & ? & ✓ \\ ? & ? & ✓ \\ ✓ & ? & ✓ \\ ? & ? & ? \\ ⋮ & ⋮ & ⋮ \\ ? & ? & ✓ \end{array}) (瘦) \end{aligned}

当 $N ≪ T$ 时， $Y$ 是胖矩阵：

Y = (\begin{array}{ccccccc} ? & ? & ✓ & ✓ & ✓ & \dots & ? \\ ✓ & ✓ & ✓ & ✓ & ? & \dots & ✓ \\ ? & ✓ & ? & ✓ & ? & \dots & ✓ \end{array}) (胖)

当 $N \approx T$ 时， $Y$ 是近似正方形矩阵：

Y = (\begin{array}{cccccc} ? & ? & ✓ & ✓ & \dots & ? \\ ✓ & ✓ & ✓ & ✓ & \dots & ✓ \\ ? & ✓ & ? & ✓ & \dots & ✓ \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ ? & ? & ✓ & ✓ & \dots & ✓ \end{array}) (近似正方形)

3.3 水平回归和垂直回归

依据数据缺失模式和矩阵形状分类的两种特殊组合值得特别关注，因为它们是很多文献的关注点。一是水平回归和非混淆性文献。非混淆性文献 (Rosenbaum 和 Rubin，1983；Rubin，2006；Imbens 和 Wooldridge，2009；Abadie 和 Cattneo，2018) 主要关注具有瘦矩阵 ( $N ≫ T$ )、且具有大量处理和控制单元的单处理周期块结构的矩阵。

\begin{aligned} Y = (\begin{array}{ccc} ? & ✓ & ? \\ ✓ & ? & ✓ \\ ? & ? & ✓ \\ ✓ & ? & ✓ \\ ? & ? & ? \\ ⋮ & ⋮ & ⋮ \\ ? & ? & ✓ \end{array}) \end{aligned}

非混淆性方法的一个简化处理是将上一阶段的结果与滞后的结果进行回归，并使用估计回归来预测缺失的潜在结果，将其称为水平回归。预测结果如下式：

{\hat{Y}}_{i T} = {\hat{β}}_{0} + \sum_{s = 1}^{T - 1} {\hat{β}}_{s} Y_{i s}

其中，

\hat{β} = \underset{β}{\arg min} \sum_{i : (i, T) \in O} {(Y_{i T} - β_{0} - \sum_{s = 1}^{T - 1} β_{s} Y_{i s})}^{2}

二是垂直回归和合成控制文献。合成控制文献 (Abadie 等，2010) 主要关注具有相对较胖 ( $N ≪ T$ ) 或近似正方形矩阵 ( $T \approx N$ )、且具有大量预处理周期的单处理单元块结构的矩阵。

\begin{aligned} Y = (\begin{array}{ccccccc} ? & ? & ✓ & ✓ & ✓ & \dots & ? \\ ✓ & ✓ & ✓ & ✓ & ? & \dots & ✓ \\ ? & ✓ & ? & ✓ & ? & \dots & ✓ \end{array}) \end{aligned}

Doudchenko 和 Imbens (2016)、Ferman 和 Pinto (2019) 展示了运用合成控制方法进行回归分析，将其称为垂直回归。预测结果如下式：

{\hat{Y}}_{N t} = {\hat{γ}}_{0} + \sum_{i = 1}^{N - 1} {\hat{γ}}_{i} Y_{i t},

其中，

\hat{γ} = \underset{γ}{\arg min} \sum_{t : (N, t) \in O} {(Y_{N t} - γ_{0} - \sum_{i = 1}^{N - 1} γ_{i} Y_{i t})}^{2}

3.4 固定效应和因子模型

面板数据中通常会产生个体固定效应和时间固定效应。水平回归中存在时间固定效应，垂直回归中存在个体固定效应。二者结合会产生双向固定效应。常见的双向固定效应模型如下所示：

Y_{i t} = γ_{i} + δ_{t} + ϵ_{i t}

更一般的因子模型可以写成：

Y_{i t} = \sum_{r = 1}^{R} u_{i r} v_{t r} + ε_{i t}, or Y = {U V}^{⊤} + ε

其中， $U$ 是 $N \times R$ 矩阵， $V$ 是 $R \times T$ 矩阵。

早期文献 (Anderson，1958；Goldberger，1972) 主要研究瘦矩阵，渐进估计是基于时期数目固定而单元数目增加。在现在文献中 (Bai，2003；Pesaran，2006；Moon 和 Weidner，2015；Bai 和 Ng，2017)，研究人员考虑了当 $N$ 和 $T$ 都增大时更复杂的渐近性，在归一化处理之后，允许一致地估计因子 $V$ 和载荷 $U$ 。

在这类文献中，通常假设因子 $R$ 大小是固定的，Bai 和 Ng (2002)、Moon 和 Weidner (2015) 讨论了估计秩 $R$ 的方法。Xu (2017) 将这种交互式固定效应方法应用于具有块状分配特点的矩阵补全问题。Gobillon 和 Magnac (2016)、Kim 和 Oka (2014)、以及 Hsiao 等 (2012) 均对此进行了应用。这种分块结构极大地简化了固定秩估计量的计算。但是，在更复杂的数据缺失模式中，这种方法效率不高，在计算上也不简单。

在机器学习和统计领域关于矩阵补全的文献中 (Srebro 等，2005；Candès 和 Recht，2009；Candès 和 Tao，2010；Gross，2011；Rohde 和 Tsybakov，2011)，研究人员从一个不能完全观察到的矩阵 $Y$ 出发，提出了低秩矩阵模型作为矩阵补全的基础。重点不是估计 $U$ 和 $V$ ，而是 $Y$ 的缺失元素。这些估计量依赖于正则化，特别是核范数正则化。

4. 核范数最小化估计下的矩阵补全

在不考虑协变量的情况下，将数据矩阵 $Y (N \times T)$ 模型化为：

Y = L^{*} + ε, where E [ε ∣ L^{*}] = 0

其中 $ε_{i t}$ 为测量误差。

假设 1： $ε$ 是独立于 $L^{*}$ ，且 $ε$ 相互独立且服从亚高斯分布。

此时目标就是估计矩阵 $L^{*}$ ，注意此时固定效应都包含在 $L^{*}$ 里。为了便于描述估计量，定义两个矩阵：

\begin{aligned} P_{O} (A)_{i t} = {\begin{cases} A_{i t} & if (i, t) \in O, \\ 0 & if (i, t) \notin O, \end{cases} \\ P_{O}^{⊥} (A)_{i t} = {\begin{cases} 0 & if (i, t) \in O, \\ A_{i t} & if (i, t) \notin O . \end{cases} \end{aligned}

为了估计 $L^{*}$ ，文中进一步引入矩阵范数和机器学习中的惩罚回归模型，并引入 $Γ \in R^{N \times 1} 和 Δ \in R^{T \times 1}$ 来避免对个体和时间的固定效应正则化，以此来构造估计量。

4.1 估计量

估计量的一般形式为：

L^{*} = \hat{L} + \hat{Γ} 1_{T}^{⊤} + l_{N} {\hat{Δ}}^{⊤}

其中，

(\hat{L}, \hat{Γ}, \hat{Δ}) = \arg min_{L, Γ, Δ} {\frac{1}{| O |} {∥ P_{O} (Y - L - Γ 1_{T}^{⊤} - 1_{N} Δ^{⊤}) ∥}_{F}^{2} + λ ∥ L ∥_{*}}

上式被称为核范数最小化的矩阵补全估计量 (MC-NNM)。式中 $λ$ 为规则参数 (惩罚因子)。

第一项为损失函数，用 Fröbenius 范数形式表示；
第二项为规则项，用矩阵 $L$ 的核范数来表示。文中指出了其他形式的范数由于不同的原因在这里并不适用，故选用核范数。

核范数 (Nuclear Norm) 是指矩阵奇异值的和。使用核范数的主要优点是可以使用快速凸优化程序来计算所得到的估计量，它可以用来约束低秩的。

秩可以度量相关性，而矩阵的相关性实际上有带有了矩阵的结构信息。如果矩阵之间各行的相关性很强，那么就表示这个矩阵实际可以投影到更低维的线性子空间，也就是用几个向量就可以完全表达了，它就是低秩的。

低秩矩阵每行或每列都可以用其他的行或列线性表出，可见它包含大量的冗余信息。利用这种冗余信息，可以对缺失数据进行恢复，也可以对数据进行特征提取。因为低秩矩阵是非凸的函数，在优化问题里面很难求解，那么就需要寻找它的凸近似来近似它了，而核范数就是它的凸近似。

4.2 计算估计量

首先不考虑固定效应，定义矩阵收缩算子:

{shrink}_{λ} (A) = S \tilde{Σ} R^{⊤}

以及矩阵 $L$ ：

L_{k + 1} (λ, O) = {shrink}_{\frac{λ ∣ O}{2}} {P_{O} (Y) + P_{O}^{⊥} (L_{k} (λ, O))}

通过交叉验证选择 $λ$ 的最优值。然后可以计算 $\hat{L}$ ，通过平方项中的一阶条件可以得到 $\hat{Γ}$ 和 $\hat{Δ}$ 。

5. 水平和垂直回归的关系

本节讨论了矩阵补全回归与水平回归 (无混淆性)、垂直回归 (合成控制) 和双重差分方法之间的关系，这是本文的第二个贡献。

将矩阵分块后发现，矩阵补全回归、水平回归、垂直回归、合成控制回归、弹性网估计和双重差分回归估计量之间的关系非常密切，它们都可以被看做是基于完全相同的目标函数，只是在正则化和对目标函数的参数的附加限制方面有所不同。定义目标函数：

Q (Y; R, A, B, γ, δ) = \frac{1}{| O |} {∥ P_{O} (Y - {A B}^{⊤} - γ 1_{T}^{⊤} - 1_{N} δ^{⊤}) ∥}_{F}^{2}

矩阵补全回归 (MC-NNM)：

\begin{aligned} (R^{m c - n n m}, A_{λ}^{m c - n n m}, B_{λ}^{m c - n n m}, γ_{λ}^{m c - n n m}, δ_{λ}^{m c - n n m}) = \\ \underset{R, A, B, γ, δ}{argmin} {Q (Y; R, A, B, γ, δ) + \frac{λ}{2} ∥ A ∥_{F}^{2} + \frac{λ}{2} ∥ B ∥_{F}^{2}} \end{aligned}

水平回归：

(R^{h r}, A^{h r}, B^{h r}, γ^{h r}, δ^{h r}) = \underset{R, A, γ, δ}{argmin} Q (Y; R, A, B, γ, δ),

subject to \begin{aligned} R = T - 1, A = (\begin{array}{c} Y_{0} \\ y_{2}^{⊤} \end{array}), γ = 0 \\ δ_{1} = δ_{2} = \dots = δ_{T - 1} = 0 \end{aligned}

垂直回归：

(R^{v t}, A^{v t}, B^{v t}, γ^{v t}, δ^{v t}) = \underset{R, A, B, γ, δ}{argmin} Q (Y; R, A, B, γ, δ),

subject to \begin{aligned} R = N - 1, B = (\begin{array}{c} Y_{0}^{⊤} \\ y_{1}^{⊤} \end{array}), \\ γ_{1} = γ_{2} = \dots = γ_{N - 1} = 0, δ = 0, \end{aligned}

合成控制：

\begin{matrix} (R^{s c - a d h}, A^{s c - a d h}, B^{s c - a d h}, γ^{s c - a d h}, δ^{s c - a d h}) \\ = \underset{R, A, B, γ, δ}{argmin} Q (Y; R, A, B, γ, δ) \end{matrix}

subject to \begin{aligned} R = N - 1, B = (\begin{array}{c} Y_{0}^{⊤} \\ y_{1}^{⊤} \end{array}), δ = 0, γ = 0, \\ \forall i, A_{i T} \geq 0, \sum_{i = 1}^{N - 1} A_{i T} = 1, \end{aligned}

弹性网络垂直回归：