Stata：工具变量的秩检验-bootrantest

发布时间：2022-02-28 阅读 1888

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下载 - 推文合集

作者：卫明昊 (中山大学)
邮箱：weimh7@mail2.sysu.edu.cn

编者按：本文主要摘译自下文，特此致谢！
Source：Chen Q, Fang Z, Huang X. Implementing an Improved Test of Matrix Rank in Stata[J]. arXiv preprint arXiv:2108.00511, 2021. -PDF-

1. 简介

1.1 工具变量

OLS 能成立的最重要条件是解释变量与扰动项不相关，否则 OLS 估计量将是不一致的，即无论样本容量多大，OLS 估计量都不会收敛至真实的总体参数。在实证研究中我们往往通过引入工具变量来解决这一问题。

对于一个内生变量，我们可以找一个变量 $z$ ，只要能满足以下两个条件，就能将其当做 $x$ 的一个工具变量：

相关性： $z$ 与 $x$ 相关，即 $Cov (z, x) \neq 0$ ；
外生性： $z$ 与 $u$ 不相关，即 $Cov (z, u) = 0$ 。

接下来，我们推导使用工具变量法时回归模型中各参数的估计值。假设回归模型为：

y_{i} = β_{1} x_{i 1} + β_{2} x_{i 2} + . . . + β_{k} x_{i k} + u_{i}

其中 $x_{k}$ 是内生变量，因此 OLS 估计量是不一致的。假设 $w$ 是 $x_{k}$ 的有效工具变量，同时由于其他解释变量 $x_{1}, x_{2}, . . ., x_{k - 1}$ 是外生的，故可以把自己作为自身的工具变量。

记解释向量 $x_{i} x_{i} = (x_{i 1}, x_{i 2}, . . ., x_{i k})^{T}$ ，则原模型为 $y_{i} = x_{i}^{T} x_{i}^{T} β β + u_{i}$ 。记工具向量 $z_{i} z_{i} = (z_{i 1}, z_{i 2}, . . ., z_{k - 1}, w_{i})^{T}$ 。由于工具变量与扰动项不相关，即 $Cov (z_{i} z_{i}, u_{i}) = E z_{i} z_{i} u_{i} - E z_{i} z_{i} E u_{i} = 0$ ，又有 $E u_{i} = 0$ ，则 $E z_{i} z_{i} u_{i} = 0$ ，这被称为 “总体矩条件” 或 “正交条件”。由此可得：

\begin{aligned} E z_{i} u_{i} = 0 & \Rightarrow E [z_{i} (y_{i} - x_{i}^{T} β)] = 0 \\ \Rightarrow E z_{i} y_{i} = (E z_{i} x_{i}^{T}) β \\ \Rightarrow β = {(E z_{i} x_{i}^{T})}^{- 1} E z_{i} y_{i} \end{aligned}

以样本矩代替总体矩，即可得到工具变量估计量：

\hat{β β_{I V}} = (\frac{1}{n} \sum_{i = 1}^{n} z_{i} x_{i}^{T} z_{i} x_{i}^{T})^{- 1} (\frac{1}{n} \sum_{i = 1}^{n} z_{i} z_{i} y_{i})

1.2 秩条件与阶条件

细心的读者可以发现，在上一部分推导中，我们没有特别说明矩阵 $E z_{i} x_{i}^{T} z_{i} x_{i}^{T}$ 是否可逆。事实上如果这个矩阵不可逆，则无法得到工具变量的估计量 $\hat{β β_{I V}}$ 。

矩阵 $E z_{i} x_{i}^{T} z_{i} x_{i}^{T}$ 可逆等价于该矩阵满秩，即 $r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k$ ，称此条件为秩条件。在满足秩条件的情况下，可以推导出在一定的正则条件下， $\hat{β β_{I V}}$ 是 $β β$ 的一致估计量，且 $\hat{β β_{I V}}$ 服从渐近正态分布。

从直观上理解，秩条件 $r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k$ 成立意味着，工具变量 $w$ 与内生变量 $x_{k}$ 相关。以一元回归为例，此时， $k = 2$ ， $x_{i} x_{i} = (1, x_{i 2})^{T}$ ， $z_{i} z_{i} = (1, w_{i})^{T}$ ，则：

\begin{aligned} E z_{i} x_{i}^{T} z_{i} x_{i}^{T} & ＝ E [\begin{array}{c} (\begin{matrix} 1 \\ w_{i} \end{matrix}) (\begin{matrix} 1 & x_{i 2} \end{matrix}) \end{array}] \\ = E [\begin{array}{c} 1 & x_{i 2} \\ w_{i} & w_{i} x_{i 2} \end{array}] \\ = [\begin{array}{c} 1 & E x_{i 2} \\ E w_{i} & E w_{i} x_{i 2} \end{array}] \end{aligned}

因此，

\begin{aligned} r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k = 2 & \Leftrightarrow | \begin{array}{c} 1 & E x_{i 2} \\ E w_{i} & E w_{i} x_{i 2} \end{array} | \neq 0 \\ \Leftrightarrow E w_{i} x_{i 2} - E w_{i} E x_{i 2} \neq 0 \\ \Leftrightarrow C o v (w_{i}, x_{i 2}) \neq 0 \end{aligned}

即 $w$ 与 $x_{k}$ 相关。

显然，满足秩条件的必要条件是在 $z_{i} z_{i}$ 中至少包含 $k$ 个变量，即不在方程中出现的工具变量个数不能少于方程中内生解释变量的个数，称此条件为阶条件。根据是否满足阶条件可分为三种情况：

不可识别：工具变量个数小于内生变量个数；
恰好识别：工具变量个数等于内生变量个数；
过度识别：工具变量个数大于内生变量个数。

1.3 工具变量的秩检验

使用工具变量法的前提之一是秩条件成立，即 $r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k$ (满列秩)。其中 $z_{i} z_{i} = (z_{i 1}, z_{i 2}, . . .,, z_{l})^{T}$ ( $l$ 个工具变量)， $x_{i} x_{i} = (x_{i 1}, x_{i 2}, . . ., x_{i k})^{T}$ ( $k$ 个解释变量)， $z_{i} z_{i}$ 与 $x_{i} x_{i}$ 可有重叠元素，且 $l ≧ k$ (满足阶条件)。

对于秩条件是否成立，可进行 “不可识别检验”。其原假设及备择假设为：

\begin{aligned} H_{0} & : r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k - 1 \\ H_{1} & : r a n k (E z_{i} x_{i}^{T} z_{i} x_{i}^{T}) = k \end{aligned}

在同方差假设下，可以使用 Anderson LM 统计量 (Anderson，1951)，其渐进分布为 $χ^{2} (L - K + 1)$ 。如果允许存在异方差，则应使用 Kleibergen-Paap rk LM 统计量 (Kleibergen-Paap，2006)，其渐进分布同样为 $χ^{2} (L - K + 1)$ 。在 Stata 中使用 ranktest 命令即可执行工具变量秩检验。

1.4 传统秩检验方法的缺陷

可以看出，上文提及的秩检验方法忽略了 $rank(E\pmb{z_ix_i^T})

2. 改良秩检验方法

考虑到现有秩检验方法的缺陷，Chen 和 Fang (2019) 提出了一种改良秩检验方法。由于这部分涉及较多数学推导，读者可能感觉晦涩难懂。因此，我们将先用较为通俗的语言简单介绍这一方法。具体如下：

首先，对于一个未知矩阵 $Π_{0} \in M M^{m \times k}$ ( $m ≧ k$ ， $M M^{m \times k}$ 代表所有 $m \times k$ 矩阵的集合)，通过修改检验的原假设，将 $rank(\Pi_0)
然后，对 $Π_{0}$ 进行奇异值分解。由于 $Π_{0}$ 有 $r a n k (Π_{0})$ 个非零奇异值，这一步将秩检验简化为检验矩阵的非零奇异值个数 (同理，也可以检验等于零的奇异值个数)；
接下来，对于 $Π_{0}$ 的估计量 ${\hat{Π}}_{n}$ ，可以构造一个统计量 $n ϕ ({\hat{Π}}_{n})$ 。这个统计量的直观意义是 ${\hat{Π}}_{n}$ 最小的 $t$ 个 ( $t$ 的取值依赖于原假设) 奇异值是否足够接近零。如果 $n ϕ ({\hat{Π}}_{n})$ 很小，就不能拒绝原假设；
最后，通过 bootstrap 方法得到统计量的分布，并根据给定的置信水平确定临界值。这样才能判断 $n ϕ ({\hat{Π}}_{n})$ 是 “大” 还是 “小”。

下面是具体的数学推导。

2.1 修改原假设

Chen 和 Fang 将原假设和备择假设修改为：

\begin{aligned} H_{0} & : r a n k (Π_{0}) ≦ r \\ H_{1} & : r a n k (Π_{0}) > r \end{aligned}

其中， $Π_{0} \in M M^{m \times k}$ ( $m ≧ k$ ) 是一个未知矩阵，在秩检验中，它代表 2SLS 方法一阶段回归中内生变量对外生变量回归的系数矩阵。 $r$ 一般取 $k - 1$ 。

2.2 奇异值分解

对 $Π_{0}$ 进行奇异值分解 (SVD)，得到：

\begin{aligned} Π_{0} = P_{0} Σ_{0} Q_{0}^{T} (1) \end{aligned}

其中 $P_{0} \in M M^{m \times m}$ 和 $Q_{0} \in M M^{k \times k}$ 都是对称矩阵， $Σ_{0} \in M M^{m \times k}$ 是一个对角矩阵，对角线上从大到小排列着 $Π_{0}$ 的奇异值，也就是 $Π_{0}^{T} Π_{0}$ 的特征值的平方根。

由于 $P_{0}$ 和 $Q_{0}$ 都是可逆矩阵， $Π_{0}$ 的秩就等于 $Σ_{0}$ 的秩，因此只需基于 $Σ_{0}$ 进行秩检验。最后，记 $r_{0} = r a n k (Π_{0})$ ，也就是说 $r_{0}$ 是 $Π_{0}$ 真实的秩。对 (1) 式等号右边三个矩阵展开得到：

\begin{aligned} [\begin{array}{c} P_{0, 1}, P_{0, 2} \end{array}] [\begin{array}{c} σ_{0, 1} & \dots & 0 & 0 & \dots & 0 \\ ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & σ_{0, r_{0}} & 0 & \dots & 0 \\ 0 & \dots & 0 & σ_{0, r_{0} + 1} & \dots & 0 \\ ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 & 0 & \dots & σ_{0, k} \\ 0 & \dots & 0 & 0 & \dots & 0 \\ ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 & 0 & \dots & 0 \end{array}] {[\begin{array}{c} Q_{0, 1}, Q_{0, 2} \end{array}]}^{T} (2) \end{aligned}

其中 $σ_{0, j} = σ_{j} (Π_{0})$ 代表 $Π_{0}$ 第 $j$ 大的奇异值， $P_{0, 1}$ 包括 $P_{0}$ 的前 $r_{0}$ 列向量，分别对应着 $Π_{0}$ 的非零奇异值， $P_{0, 2}$ 则是 $P_{0}$ 的后 $k - r_{0}$ 列向量。同理， $Q_{0, 1}$ 和 $Q_{0, 2}$ 代表着对应的向量。

此时，原假设成立，当且仅当最小的 $k - r$ 个 $Π_{0}$ 的奇异值为 0。也就是说，原假设和备择假设等价于：

\begin{aligned} H_{0} & : ϕ (Π_{0}) = 0 \\ H_{1} & : ϕ (Π_{0}) > 0 \end{aligned}

其中， $ϕ (Π_{0}) = \sum_{j = r + 1}^{k} σ_{j}^{2} (Π_{0})$ ，表示最小的 $k - r$ 个 $Π_{0}$ 的奇异值平方之和。

2.3 构造统计量

对于 $Π_{0}$ 的估计量 ${\hat{Π}}_{n}$ ，可以构造一个统计量 $n ϕ ({\hat{Π}}_{n})$ 。如果该统计量大于临界值，就可以拒绝原假设。在一阶段回归中， ${\hat{Π}}_{n}$ 就是 OLS 估计量。为了得到临界值，Chen 和 Fang 证明了：

\begin{aligned} n ϕ ({\hat{Π}}_{n}) \overset{L}{\to} \sum_{j = r - r_{0} + 1}^{k - r_{0}} σ_{j}^{2} (P_{0, 2}^{T} M Q_{0, 2}) (3) \end{aligned}

其中 $\overset{L}{\to}$ 表示依分布收敛， $M \in M^{m \times k}$ 是 ${\hat{Π}}_{n}$ 的渐进分布，即 $\sqrt{n} ({\hat{Π}}_{n} - Π_{0}) \overset{L}{\to} M$ 。因此，就像我们用标准正态分布的分位数作为 $t$ 检验的临界值，可以把 (3) 式中分布的分位数当做 $n ϕ ({\hat{Π}}_{n})$ 的临界值。问题是 $M, P_{0, 2}, Q_{0, 2}$ 和 $r_{0}$ 都是未知的。但 (3) 式同时说明可以用它们的估计量代替这些未知量。

2.4 得到统计量的分布并确定临界值

$M$ 的分布的估计可以通过 bootstrap 来完成。我们直接用 ${\hat{M}}_{n}^{*} = \sqrt{n} ({\hat{Π}}_{n}^{*} - Π_{0})$ 的分布来代替即可。其中， ${\hat{Π}}_{n}^{*}$ 是基于 bootstrap 样本的一阶段回归的 OLS 估计量。

Chen 和 Fang 提出了两种方法来得到 $r_{0}$ 的估计量。其中一种方法是统计 ${\hat{Π}}_{n}$ 中 “不等于” 零的奇异值个数。为此，我们要设置一个参数 $κ_{n} = n^{- 1 / 2}$ (或者 $n^{- 1 / 4}$ )，则 $r_{0}$ 的估计量为：

\begin{aligned} {\hat{r}}_{n} = max {j = 1, \dots, r : σ_{j} ({\hat{Π}}_{n}) ≧ κ_{n}} (4) \end{aligned}

在得到 $r_{0}$ 的估计量 ${\hat{r}}_{n}$ 后，对 ${\hat{Π}}_{n}$ 进行奇异值分解就能得到 $P_{0, 2}$ 和 $Q_{0, 2}$ 的估计量 ${\hat{P}}_{n, 2}$ 和 ${\hat{Q}}_{n, 2}$ 。

最后，对于给定的显著性水平 $α \in (0, 1)$ ，我们将下面这个分布的 $1 - α$ 分位数作为临界值，记为 ${\hat{c}}_{n, 1 - α}$ 。

\begin{aligned} \sum_{j = r - {\hat{r}}_{n} + 1}^{k - {\hat{r}}_{n}} σ_{j}^{2} ({\hat{P}}_{n, 2}^{T} {\hat{M}}_{n}^{*} {\hat{Q}}_{n, 2}) (5) \end{aligned}

如果统计量 $n ϕ ({\hat{Π}}_{n}) > {\hat{c}}_{n, 1 - α}$ ，则拒绝原假设。我们把这种方法称为分析法 (analytic approach) 。

另一种得到 ${\hat{r}}_{n}$ 的方法是基于 Kleibergen-Paap rk LM 统计量指定一个置信水平 $β \in (0, α)$ 。此时的估计值有 $1 - β$ 的概率与 $Π_{0}$ 真实的秩 $r_{0}$ 一致。

如果 ${\hat{r}}_{n} > r$ ，就拒绝原假设，否则就将 ${\hat{c}}_{n, 1 - α + β}$ 作为临界值。如果 $n ϕ ({\hat{Π}}_{n}) > {\hat{c}}_{n, 1 - α + β}$ ，则拒绝原假设。我们把这种方法称为两步法 (two-step approach)。

3. 命令介绍

基于上部分介绍的改良秩检验方法，Chen 等 (2021) 开发了 bootranktest 命令，下面对这一命令进行介绍。

bootranktest 命令可以直接通过 ssc insatll bootranktest, replace 进行安装。(注：由于该论文目前还在 The Stata Journal 的 RR 中，命令暂时还不能通过 SSC 进行安装，请读者耐心等待。)

bootranktest 命令语法如下：

bootranktest (varlist1) (varlist2) [weight] [if exp] [in range]
             [, rank(#) allrank numboot(#) beta(#) kappan(#)
             blocksize(#) partial(varlist3) cluster(varname) noconstant cfa]

varlist1：工具变量；
varlist2：内生变量；
varlist3：回归方程中的其他外生变量；
rank：假设矩阵的秩为 $r$ ，默认 $r = k - 1$ ，注意 $r$ 必须小于 $k$ ；
allrank：对于 $r = 0, \dots, k - 1$ ，执行命令并报告结果；
numboot：bootstrap 次数，默认为 1000；
beta：设置两步法参数 $β$ ，默认 $β = 0.05 / 10$ ；
kappan：设置分析法参数 $κ_{n}$ ，默认 $κ_{n} = n^{- 1 / 4}$ ；
blocksize：对时间序列数据进行 bootstrap 时，每次抽样数据时间段长度；
partial(varlist3)：指定回归方程中其他非常数项外生变量；
cluster：指定聚类依据的变量；
nonconstant：一阶段回归方程中没有常数项；
cfa：报告分析法和两步法的具体结果。

4. 案例演示

下面使用 Stata 自带的 klein 数据集对 bootranktest 命令进行演示。klein 数据集包含 22 个时间序列观测值 (1920-1941 年)。主要变量为：

消费 consumption，私人收益 profits，美国工资总额 wagetot，政府支出 govt，间接巴士税加净出口 taxnetx，年份减去 1931 year，政府工资总额 wagegovt，股本的滞后项 capital1 和总需求 totinc。

计量模型为：

c o n s u m p t i o n_{t} = β_{0} + β_{1} p r o f i t s_{t - 1} + β_{2} p r o f i t s_{t} + β_{3} w a g e t o t_{t} + e_{t}

我们假设 profits 的滞后项为外生变量，profits 和 wagetot 为内生变量。工具变量是 govt，taxnetx，year，wagegovt，capital1 和 totinc 的滞后项。一阶段回归方程为：

[\begin{matrix} {profits}_{t} \\ {wagetot}_{t} \end{matrix}] = Π_{0}^{⊤} [\begin{matrix} {govt}_{t} \\ {taxnetx}_{t} \\ {year}_{t} \\ {wagegovt}_{t} \\ {capital1}_{t} \\ {totinc}_{t - 1} \end{matrix}] + Γ_{0}^{⊤} [\begin{matrix} 1 \\ {profits}_{t - 1} \end{matrix}] + u_{t},

我们的目标是检验工具变量的秩条件是否成立，即 $Π_{0}$ 的秩是否为 2 (满列秩)。原假设为：

H_{0} : r a n k (Π_{0}) ≦ 1

接下来使用 bootranktest 命令进行秩检验，得到如下结果：

在 5% 置信水平下，使用两步法得到的 $p$ 值为 0.03，拒绝原假设；使用分析法得到的