Stata：控制函数法及其Stata应用

发布时间：2022-12-31 阅读 1890

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

作者：张弛 (中山大学)
邮箱：zhangch369@mail2.sysu.edu.cn

1. 何时选择控制函数法？

在对政策效果进行评估时，选择参与到一个项目的主体可能在某些方面与选择不参与到这个项目的主体存在系统性差异，并且这样的差异很难被观察到，这会导致内生性问题。

例如，Sapelli and Vial (2002) 在研究智利上公立学校和私立学校对二年级学生成绩影响时，不同的父母教育程度、不同的家庭经济状况等因素，一方面会直接影响到学生的成绩，另一方面还会影响到学生是选择读公立学校还是私立学校。因此，“选择学校类型”这一变量本身就存在内生性。

控制函数法的核心思想就在于，通过对内生解释变量 (Endogenous Explanary Variables，下文简称 EEV) 的内生性来源进行提取，然后加入原始模型进行估计，来获得对政策效果的无偏估计。

2. 常系数线性模型

2.1 EEV 为连续型变量

2.1.1 原理简介

假设 $y_{1}$ 为被解释变量， $y_{2}$ 为 EEV， $z$ 是 $1 \times L$ 的外生变量 (其中， $z_{1} = 1 ， z_{1}$ 是 $z$ 的子向量)，我们估计的原始方程为：

y_{1} = z_{1} δ_{1} + α_{1} y_{2} + u_{1} (1)

外生性假设为：

E (z^{'} u_{1}) = 0 (2)

对 $y_{2}$ 的内生性假设：

y_{2} = z π_{2} + v_{2}, E (z^{'} v_{2}) = 0 (3)

其中， $π_{2}$ 是 $L \times 1$ 的向量。此时，定义残差项之间的关系：

u_{1} = ρ_{1} v_{2} + e_{1} (4)

ρ_{1} = E (v_{2} u_{1}) / E (v_{2}^{2}) (5)

其中， $ρ_{1}$ 为总体回归系数。将 (4) 代入 (1)，得：

y_{1} = z_{1} δ_{1} + α_{1} y_{2} + ρ_{1} v_{2} + e_{1} (6)

此时，该方程所估计的 $α_{1}$ 即与残差项无关。因为我们将 $y_{2}$ 中与残差项有关的，存在内生性的部分 $v_{2}$ 已单独提取出，作为一个新变量加入我们的估计模型中。

2.1.2 检验变量是否为 EEV

y_{i 1} = z_{i 1} δ_{1} + α_{1} y_{i 2} + ρ_{1} {\hat{v}}_{i 2} + {error}_{i} (7)

原假设 $H_{0} : ρ_{1} = 0$ 相当于 Hausman 检验中

H_{0} : Cov (y_{2}, u_{1}) = 0

在这样的线性方法中，对 $δ_{1}$ 和 $α_{1}$ 的估计和 2SLS 方法是一样的。注意 IV 方法和 CF 方法不一定等价。例如，原方程加入 $y_{2}$ 的平方项

y_{1} = z_{1} δ_{1} + α_{1} y_{2} + γ_{1} y_{2}^{2} + u_{1} (8)

若我们使用IV方法，我们需要加入 $z_{2}^{2}$ 作为 $y_{2}^{2}$ 的工具变量，此时工具变量向量为 $(z_{1}, z_{2}, z_{2}^{2})$ 。但是若使用 CF 方法，则需要施加更强的假设：

E (u_{1} ∣ z, y_{2}) = E (u_{1} ∣ v_{2}) = ρ_{1} v_{2} (9)

在 CF 方法中，第一阶段回归为：

y_{2} = z π_{2} + v_{2} (10)

在前述条件下，我们有新的估计方程 (11)，此时，与 2SLS 的估计方法不一致。

E (y_{1} ∣ z, y_{2}) = z_{1} δ_{1} + α_{1} y_{2} + γ_{1} y_{2}^{2} + ρ_{1} v_{2} (11)

用第一阶段回归所估计得到的 ${\hat{v}}_{2}$ 代入上式，用 OLS 方法估计。

这样的假设会使得 CF 方法更有效，但是稳健性不如 IV 方法，因为 IV 方法只需要满足 $E (u_{1} ∣ z) = 0$ ，这个假设更弱。

2.1.3 估计方法和 Stata 实现

第一步：将 $y_{2}$ 对 $z$ 做回归，获得简化方程的残差项 ${\hat{v}}_{2}$ 。第二步：将 $y_{1}$ 对 $z_{1}$ ， $y_{2}$ ， ${\hat{v}}_{2}$ 做回归，对 (6) 式的 OLS 估计就是 Control Function 的估计方法。在第二步的估计中， $δ_{1}$ ， $α_{1}$ ， $ρ_{1}$ 都是一致的。

. bcuse wage2 //一个关于教育回报的数据集，作为演示的范例
. gen exper2=exper^2 
. sum educ
. gen educ_average=r(mean)

. // 控制函数法stage1 获取存在内生性的残差项
. bootstrap, reps(1000): reg educ meduc feduc black
. predict v2, residual  //残差项吸收了内生性

. // 教育回报率方程其他变量的设置
. gen black_educ=black*(educ-educ_average)

. // stage2 将内生的残差项加入原始模型回归
. bootstrap, reps(1000): reg lwage educ exper exper2 black tenure black_educ v2

Linear regression                                       Number of obs =    722
                                                        Replications  =  1,000
                                                        Wald chi2(7)  = 165.62
                                                        Prob > chi2   = 0.0000
                                                        R-squared     = 0.1824
                                                        Adj R-squared = 0.1744
                                                        Root MSE      = 0.3811
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
       lwage | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        educ |      0.120      0.015     7.96   0.000        0.090       0.149
       exper |      0.008      0.016     0.50   0.620       -0.024       0.040
      exper2 |      0.000      0.001     0.66   0.508       -0.001       0.002
       black |     -0.149      0.057    -2.61   0.009       -0.261      -0.037
      tenure |      0.009      0.003     3.18   0.001        0.004       0.015
  black_educ |     -0.030      0.026    -1.19   0.234       -0.080       0.020
          v2 |     -0.056      0.016    -3.59   0.000       -0.087      -0.026
       _cons |      4.946      0.230    21.48   0.000        4.495       5.398
------------------------------------------------------------------------------

2.2 EEV 为离散型变量

2.2.1 原理方法

CF 方法相比于 IV 方法的优势在于，它可以针对非线性的 EEV 进行一致性估计，例如 EEV 可以是一个二值变量。当 $y_{2}$ 为非线性形式时，我们无法使用 IV 估计，CF 方法是更好的选择。例如，待估模型为：

y_{1} = z_{1} δ_{1} + γ_{1} y_{2} + u_{1} (12)

在该模型中，我们假设 $y_{2}$ 为二值变量，满足以下方程：

y_{2} = 1 [z δ_{2} + e_{2} \geq 0] (13)

运用 CF 方法，需施加以下假设：

e_{2} \sim Normal (0, 1) (14)

此时， $E (u_{1} ∣ z, y_{2})$ 取决于 $(u_{1}, y_{2})$ 的联合分布，且 $(u_{1}, y_{2})$ 独立于 $z$ ， $E (u_{1} ∣ e_{2}) = ρ_{1} e_{2}$ ，此时我们估计的方程为：

E (y_{1} ∣ z, y_{2}) = z_{1} δ_{1} + α_{1} y_{2} + ρ_{1} [y_{2} λ (z δ_{2}) - (1 - y_{2}) λ (- z δ_{2})] (15)

并且

E (u_{1} ∣ z, y_{2}) = ρ_{1} [y_{2} λ (z δ_{2}) - (1 - y_{2}) λ (- z δ_{2})] (16)

其中， $λ (\cdot)$ 是逆米尔斯比率 (inverse Mills ratio，IMR)。

2.2.2 估计步骤和 Stata 实现

第一步：估计 $P (y_{2} = 1 ∣ z) = Φ (z δ_{2})$ ，获得广义残差：

{\hat{r}}_{i 2} \equiv y_{i 2} λ (z_{i} {\hat{δ}}_{2}) - (1 - y_{i 2}) λ (- z_{i} {\hat{δ}}_{2}), i = 1, \dots, N

第二步：将 $y_{i 1}$ 对 $z_{i 1}$ ， $y_{i 2}$ ， ${\hat{r}}_{i 2}$ 进行 OLS 回归，获得对 $λ_{1}$ ， $α_{1}$ ， $ρ_{1}$ 的一致估计。

. webuse drugexp, clear
. etregress lndrug chron age lninc, treat(ins=age married lninc work) poutcomes cfunction

Linear regression with endogenous treatment              Number of obs = 6,000
Estimator: Control function
------------------------------------------------------------------------------
             |               Robust
             | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
lndrug       |
       chron |      0.467      0.032    14.61   0.000        0.405       0.530
         age |      0.102      0.003    34.98   0.000        0.096       0.108
       lninc |      0.055      0.023     2.45   0.014        0.011       0.099
       1.ins |     -0.860      0.348    -2.47   0.014       -1.543      -0.177
       _cons |      1.666      0.253     6.59   0.000        1.170       2.161
-------------+----------------------------------------------------------------
ins          |
         age |      0.021      0.002     9.21   0.000        0.017       0.026
     married |      0.085      0.036     2.35   0.019        0.014       0.155
       lninc |      0.102      0.023     4.55   0.000        0.058       0.146
        work |      0.288      0.037     7.75   0.000        0.215       0.361
       _cons |     -0.623      0.109    -5.73   0.000       -0.836      -0.410
-------------+----------------------------------------------------------------
    /athrho0 |      0.404      0.172     2.34   0.019        0.066       0.742
   /lnsigma0 |      0.316      0.050     6.31   0.000        0.218       0.414
    /athrho1 |      0.793      0.299     2.66   0.008        0.208       1.378
   /lnsigma1 |      0.187      0.061     3.04   0.002        0.066       0.307
-------------+----------------------------------------------------------------
        rho0 |      0.383      0.147                         0.065       0.630
      sigma0 |      1.372      0.069                         1.243       1.513
     lambda0 |      0.525      0.226                         0.082       0.969
        rho1 |      0.660      0.169                         0.205       0.881
      sigma1 |      1.205      0.074                         1.069       1.359
     lambda1 |      0.795      0.251                         0.303       1.288
------------------------------------------------------------------------------
Wald test of indep. (rho0 = rho1 = 0): chi2(2) =    8.88  Prob > chi2 = 0.0118

. // 在该估计中，内生解释变量为ins，它受到二值变量age，married，lninc，work的影响
. // poutcomes：表示用潜在结果模型估计第一阶段回归，如果不加该选项，即表示第一阶段回归为线性回归
. // cfunction：表示用控制函数法估计处理效应

2.2.3 EEV 是非线性的情况

假设 $y_{2}$ 和其他外生变量是交互的：

y_{1} = z_{1} δ_{1} + y_{2} z_{1} γ_{1} + u_{1} (17)

如果 $y_{2} z_{1}$ 中， $y_{2}$ 是线性的，可以按照最原始的方法，如果是二值变量，可以使用离散型随机变量的方法，均为以上两种情况的简单变形。

3. 相关随机系数模型

3.1 模型简介

有时候， $y_{2}$ 的效果也会受到不可观测的因素影响。例如，就读于教会学校的教育回报率，可能会因为入学学生的个体差异而不同，且这些不同是无法观测的。

相关随机系数 (correlated random coefficient，CRC)，指允许一个随机的系数影响解释变量。这个概念最早由 Heckman and Vytlac 在1998年提出。在研究处理效应时， CRC 模型允许异质性的处理效应和自选择结合。

3.2 模型设定

假设估计模型：

y_{i 1} = z_{i 1} δ_{1} + g_{i 1} y_{i 2} + u_{i 1} (18)

其中， $y_{i}$ 与 $u_{i 1}$ ， $g_{i 1}$ 相关，并且 EEV 系数满足以下两个条件：

g_{i 1} = γ_{1} + v_{i 1} (19)

γ_{1} = E (g_{i 1}) (20)

代入 (18) 得：

y_{i} = z_{i 1} δ_{1} + γ_{1} y_{i 2} + v_{i 1} y_{i 2} + u_{i 1} (21)

与此同时，针对 (21) 式，我们还需要施加如下的外生性假设：

E (u_{i 1} ∣ z_{i}) = 0, E (v_{i 1} ∣ z_{i}) = 0 (22)

注意，此处不能运用 2SLS 估计，因为此时估计出来的 $γ_{1}$ 是不一致的，残差为 $e_{i 1} = v_{i 1} y_{i 2} + u_{i 1}$ ，仍然内生。

假设 $y_{i 2}$ 的内生性设定如下：

y_{i 2} = z_{i} π_{2} + v_{i 2}, E (v_{i 2} ∣ z_{i}) = 0 (23)

存在不可观测因素 $u_{i 1}$ , $v_{i 1}$ ，他们和 $v_{i 2}$ 满足以下线性关系：

E (u_{i 1} ∣ v_{i 2}) = η_{1} v_{i 2}, E (v_{i 1} ∣ v_{i 2}) = ψ_{1} v_{i 2} (24)

并且，我们认为，所有不可观测的变量都与 $z_{i}$ 独立，所以，我们实际上需要估计的式子是：

E (y_{i 1} ∣ z_{i}, y_{i 2}) = E (y_{i 1} ∣ z_{i 1}, y_{i 2}, v_{i 2}) = z_{i 1} δ_{1} + γ_{1} y_{i 2} + η_{1} v_{i 2} + ψ_{1} v_{i 2} y_{i 2} (25)

接着，我们将 $y_{i 1}$ 对 $z_{i 1}, y_{i 2}, {\hat{v}}_{i 2}, {\hat{v}}_{i 2} y_{i 2}$ 做 OLS 回归，得到CF估计。此处 ${\hat{v}}_{i 2} y_{i 2}$ 的系数，衡量的就是随机系数 $v_{i 1}$ 对 $y_{i 2}$ 的影响。但是这里必须注意，如果 $v_{i 2}$ 的系数不为 0，对交互项 ${\hat{v}}_{i 2} y_{i 2}$ 的 $t$ 检验在第一阶段的回归是不可靠的，需要在第二步使用 bootstrap 去调整方差。

在这里，为了方便，我们是对交互项和残差项 ( ${\hat{v}}_{i 2}, {\hat{v}}_{i 2} y_{i 2}$ ) 做联合检验，检验是否存在外生性。

3.3 操作步骤和 Stata 实现

操作步骤：与线性常系数模型中EEV连续的情况基本类似，只是在第二阶段回归的时候加入 EEV 和残差的交互项。在这里，我们依然考虑一个讨论教育回报率的模型。

. bcuse wage2, clear
. gen exper2=exper^2 
. sum educ
. gen educ_average=r(mean)

. // 控制函数法 stage1 获取存在内生性的残差项
. bootstrap, reps(1000): reg educ meduc feduc black
. predict v2, residual  // 残差项吸收了内生性

. // 教育回报率方程其他变量的设置
. gen v2_educ=v2*educ  // 构造残差和 eev 的交互项
. gen black_educ=black*(educ-educ_average)

. // stage2 将内生的残差项加入原始模型回归
. bootstrap,reps(1000): reg lwage educ exper exper2 black tenure black_educ v2 v2_educ

Linear regression                                       Number of obs =    722
                                                        Replications  =  1,000
                                                        Wald chi2(8)  = 187.97
                                                        Prob > chi2   = 0.0000
                                                        R-squared     = 0.1866
                                                        Adj R-squared = 0.1775
                                                        Root MSE      = 0.3803
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
       lwage | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        educ |      0.124      0.014     8.56   0.000        0.095       0.152
       exper |      0.007      0.016     0.42   0.672       -0.024       0.037
      exper2 |      0.001      0.001     0.83   0.407       -0.001       0.002
       black |     -0.160      0.057    -2.82   0.005       -0.272      -0.049
      tenure |      0.010      0.003     3.20   0.001        0.004       0.016
  black_educ |     -0.036      0.026    -1.37   0.170       -0.087       0.015
          v2 |      0.040      0.050     0.80   0.423       -0.058       0.137
     v2_educ |     -0.007      0.003    -2.02   0.043       -0.013      -0.000
       _cons |      4.922      0.215    22.86   0.000        4.500       5.344
------------------------------------------------------------------------------

3.4 模型拓展

(28) 式的模型也可以拓展到非线性的 ( $z_{i 1}$ , $y_{i 2}$ )，包括多次方项和交互项。更一般化的模型如下：假设存在一个 $1 \times K_{1}$ 的 $x_{i 1}$ ，存在这样一个方程 $g_{1} (z_{i 1}, y_{i 2})$ ， $x_{i 1}$ 的随机系数为 $b_{i 1} = β_{1} + v_{i 1}$ 。此时，我们估计的方程为：

y_{i 1} = α_{1} + x_{i 1} β_{1} + u_{i 1} + x_{i 1} v_{i i}

在第一阶段回归获得残差 ${\hat{v}}_{i 2}$ 后，我们将 $y_{i 1}$ 对 $1$ ， $x_{i 1}$ ， ${\hat{v}}_{i 2}$ ， $x_{i 1}$ ， ${\hat{v}}_{i 2}$ 做 OLS 回归，并使用非参数的 bootstrap 方法调整标准差，从而得到有效的估计。如果 ${\hat{Ψ}}_{1}$ 是 $x_{i 1} {\hat{v}}_{i 2}$ 的 $K_{1}$ 维的向量，我们可以估计 $E (b_{i 1} ∣ v_{i 2}) = \hat{β_{1}} + \hat{v_{i 2}} \hat{Ψ_{1}}$ 。

类似的，对于 $y_{2}$ 是二值变量的情况，该方法也一样成立。由前所示， $y_{2}$ 需满足前面的条件：

y_{2} = 1 [z δ_{2} + e_{2} > 0]

e_{2} \sim Normal (0, 1)

然后，获取标准化的残差：

{\hat{r}}_{i 2} \equiv y_{i 2} λ (z_{i} {\hat{δ}}_{2}) - (1 - y_{i 2}) λ (- z_{i} {\hat{δ}}_{2}), i = 1, \dots, N

最后，将 $y_{i 1}$ 对 $z_{i 1}, y_{i 2}, y_{i 2} \cdot (z_{i 1} - {\bar{z}}_{1}), {\hat{r}}_{i 2}, y_{i 2} {\hat{r}}_{i 2}$ 做回归，中心化的 $z_{i 1}$ 帮助我们求出政策的平均效应。注意，针对 (26) 式的估计，我们不直接讨论 $y_{2} = 1$ 和 $y_{2} = 0$ 两种情况下分别的效果，ATE (Average Treatment Effect) 在这个模型下就是 $γ_{1}$ 。

E (y_{i 1} ∣