Stata: 动态面板门槛模型

发布时间：2020-03-13 阅读 19021

Stata 连享会主页 || 视频 || 推文

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

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

✨ 课程主页： https://gitee.com/arlionn/Course

严雅雪 (中南财经政法大学)
E-mail: 535743919@qq.com)

1、动态门槛效应回归模型背景简介

在线性回归中，所有观测样本回归函数相同，那如果所有观测样本是离散类，那回归函数的差异性将如何体现?

在回归分析中，我们常常关心系数估计值是否稳定，即如果将整个样本分成若干子样本分别进行回归，是否还能得到大致相同的估计系数?

如果样本的变量不是离散型变量，而是连续型变量，就需要给定一个标准，即 “门槛值” 。

传统的做法是，研究者主观的确定一个门槛值，然后将门槛值帮样本进行分割，即不对门槛值进行参数估计，也不对其显著性进行统计估计检验，所以很显然，这样得到的结果并不可靠。

Hansen 于 1996 年在《Econometrica》上发表文章《Inference when a nuisance parameter is not identified under the null hypothesis》，提出了时间序列门槛自回归模型 (TAR) 的估计和检验。之后，他在门槛模型上连续追踪，发表了几篇经典文章，尤其是 1999 年的《Threshold effects in non-dynamic panels: Estimation, testing and inference》，2000年的《Sample splitting and threshold estimation》和 2004 年与他人合作的《Instrumental Variable Estimation of a Threshold Model》。Hansen (1999) 提出的门槛效应模型是为存在个体固定效应非动态面板所设计的，门槛值的可以利用固定效应转换，通过最小二乘估计来得到回归斜率。利用一个非标准的渐近理论进行推断，并允许构建置信区间来进行假设测试。

Hansen (2000) 假设样本数据为

{y_{i t}, q_{i t}, x_{i t} : 1 \leq i \leq n, 1 \leq t \leq T}

表达式为

y_{i t} = μ_{i} + β_{1}^{'} x_{i t} I (q_{i t} \leq γ) + β_{2}^{'} x_{i t} I (q_{i t} > γ) + e_{i t}

代表划分样本的门槛变量，为待估计的门槛值，为外生解释变量，与扰动项不相关。I 为指示函数，即括号中的表达式为真，则取值为1；反之，取 0。写成分段函数为：

y_{i t} = {\begin{array}{l} μ_{i} + β_{1}^{'} x_{i t} + e_{i t}, & q_{i t} \leq γ \\ μ_{i} + β_{2}^{'} x_{i t} + e_{i t}, & q_{i t} > γ \end{array}

但在这些文章中，均基于的假设是协变量必须为强外生变量，否则会形成有偏估计。而在现实中较难实现这一假设，故无法扩展应用领域。

Caner and Hansen (2004) 年采用较为间接的方法解决了这个问题。他们研究了带有内生性的解释变量和一个外生门槛变量的面板门槛模型，只要利用简化型对内生变量进行一定的处理，然后用 2SLS (两阶段最小二乘法) 或者 GMM (广义矩估计) 对参数进行估计, 即放松了对内生解释变量的要求，但对于门槛变量的要求仍是必须强外生变量。

本文即采用这个方法来运用动态门槛效应模型进行表征。Hansen (1999) 中的具有阈值效应的面板模型已被广泛用于实证研究。

Wang (2015) 开发了 Stata 命令 xthreg，他的固定效应估计量要求协变量必须是强烈外生的才能使估计量保持一致。但是，强外生性可能在许多实际应用中受到限制。因此，Seo 和 Shin (2016) 将模型扩展为具有潜在内生阈值变量的动态面板模型。虽然先前的 xthreg 命令计算固定效果估计量，因此在此一般设置下不一致，但是的 xthenreg 命令产生一致且渐近的正态估计。Seo M H 和 Kim S (2019) 提出了一种计算上更有吸引力的自举算法来实现线性测试，而非 Seo 和 Shin (2016) 最初提出的非参数 i.i.d 自举算法。此外，我们还提出了约束 GMM 估计量，该估计量反映了近年来更为流行的扭力限制，如 Zhang 等人 (2017) 给出了其渐近方差公式和一致估计量。

2、用 STATA 运行动态门槛效应命令-基于 Wang (2015) 方法

2.1 命令的安装示例数据导入

在 Stata 命令窗口执行第一行代码即可完成对 xthreg 然后输入第二行命令下载中国工业经济面板数据 http://www.ciejournal.org/Magazine/Show?id=69845 数据.dta, 并执行第三行命令导入数据。

ssc install xthreg, replace //安装命令
net get xthreg  //下载命令附带的数据到当前工作路径下
use "数据.dta", replace

2.2 命令的语法

该命令的语法如下：

* xthreg -- Estimate fixed-effect panel threshold model
  xthreg depvar [indepvars] [if] [in], rx(varlist) qx(varname)    
             [thnum(#) grid(#) trim(numlist) bs(numlist) 
             thlevel(#) gen(newvarname) noreg nobslog
             thgiven options]

各个选项的含义说明如下：

depvar，必选项，为因变量。
indepvars，可选项，为自变量.
rx(varlist) 必选项，核心解释变量。允许使用时间序列运算符。
qx(varname) 必选项，是门槛变量。允许使用时间序列运算符。
thnum(#) 是阈值的数量。在当前版本 (Stata 13) 中，#必须等于或小于 3。默认值是 thnum(1)。
grid(#)是网格点的数量。使用 grid() 可以避免在计算大样本时花费太多时间。默认值是 grid(300)。
trim(numlist)是估计每个阈值的修剪比例。修剪比例的数量必须等于 thnum() 中指定的阈值的数量。所有阈值的缺省值都是 trim(0.01)。例如，为了适应三阈值模型，您可以设置 trim(0.01 0.01 0.05)。
bs(numlist) 是引导 Bootstrap 的数量。如果没有设置 bs()， xthreg 不会使用 bootstrap 进行阈值效应测试
thlevel(#) 指定阈值的置信区间的置信级别 (以百分比表示)。默认值是 thlevel(95)。

2.3 命令操作

2.3.1 拟合值估计

首先，本文将基于 xtdpdsys 命令和 nlswork.dta 数据对因变量 api 和核心解释变量 is 进行拟合值估计，代码如下：

gen api1=L. api
gen is1=L.is
xtdpdsys wage1 grade1 age race hours,twostep

运行结果如下：

System dynamic panel-data estimation            Number of obs     =      2,095
Group variable: city                            Number of groups  =        271
Time variable: year
                                                Obs per group:
                                                              min =          1
                                                              avg =   7.730627
                                                              max =          8

Number of instruments =     39                  Wald chi2(4)      =    1690.05
                                                Prob > chi2       =     0.0000
Two-step results
------------------------------------------------------------------------------
        api1 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        api1 |
         L1. |   .6248146   .0215764    28.96   0.000     .5825256    .6671036
             |
         is1 |   .0046254   .0028238     1.64   0.101    -.0009092      .01016
       lngrp |  -.1840629   .0725119    -2.54   0.011    -.3261836   -.0419422
       lnfdi |   .0001234   .0134864     0.01   0.993    -.0263095    .0265564
       _cons |   1.950786   .4309568     4.53   0.000     1.106127    2.795446
------------------------------------------------------------------------------
Warning: gmm two-step standard errors are biased; robust standard
         errors are recommended.
Instruments for differenced equation
        GMM-type: L(2/.).api1
        Standard: D.is1 D.lngrp D.lnfdi
Instruments for level equation
        GMM-type: LD.api1
        Standard: _cons```

然后我们进行 sargan 检验，和获取拟合值，代码如下：

estat sargan predict hapi

运行结果如下：

```stata

Sargan test of overidentifying restrictions
        H0: overidentifying restrictions are valid

        chi2(34)     =  83.05857
        Prob > chi2  =    0.0000

2.3.2代入拟合值估计门槛值和置信区间

首先，我们将科技强度 stei 作为门槛变量，考察在科技强度下，产业结构对空气质量的影响程度。代码如下：

xthreg api hapi lngrp pd, rx(is) qx(stei) thnum(3) grid(400) trim(0.01 0.01 0.01) bs(300 300 300)

运行结果如下：

Threshold estimator (level = 95):
-----------------------------------------------------
     model |    Threshold         Lower         Upper
-----------+-----------------------------------------
      Th-1 |       0.1315        0.1268        0.1353
     Th-21 |       0.1315        0.1268        0.1353
     Th-22 |       0.4464        0.4419        0.4491
      Th-3 |       1.0693        1.0650        1.0743
-----------------------------------------------------

Threshold effect test (bootstrap = 300 300 300):
-------------------------------------------------------------------------------
 Threshold |       RSS        MSE      Fstat    Prob   Crit10    Crit5    Crit1
-----------+-------------------------------------------------------------------
    Single |  335.0607     0.1655      17.45  0.1533  19.7967  24.1998  30.3118
    Double |  333.8384     0.1649       7.41  0.3567  11.1727  13.0692  17.5784
    Triple |  333.0669     0.1646       4.69  0.8467  14.7811  17.0482  23.5613
-------------------------------------------------------------------------------

Fixed-effects (within) regression               Number of obs      =      2032
Group variable: id                              Number of groups   =       254

R-sq:  within  = 0.1001                         Obs per group: min =         8
       between = 0.5008                                        avg =       8.0
       overall = 0.3400                                        max =         8

                                                F(7,1771)          =     28.14
corr(u_i, Xb)  = 0.4123                         Prob > F           =    0.0000

------------------------------------------------------------------------------
         api |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        hapi |   .3679484   .0368293     9.99   0.000     .2957149    .4401818
       lngrp |   .2459264   .0878522     2.80   0.005     .0736215    .4182313
          pd |  -.0000987   .0002223    -0.44   0.657    -.0005348    .0003373
             |
   _cat#c.is |
          0  |  -.0067473   .0037753    -1.79   0.074    -.0141519    .0006572
          1  |  -.0127042   .0033193    -3.83   0.000    -.0192143    -.006194
          2  |  -.0149214   .0032984    -4.52   0.000    -.0213905   -.0084522
          3  |  -.0135887    .003351    -4.06   0.000     -.020161   -.0070164
             |
       _cons |   .6807577   .6063483     1.12   0.262    -.5084758    1.869991
-------------+----------------------------------------------------------------
     sigma_u |  .53213522
     sigma_e |  .43366387
         rho |  .60090941   (fraction of variance due to u_i)
------------------------------------------------------------------------------
F test that all u_i=0: F(253, 1771) = 3.72                   Prob > F = 0.0000

3、用 STATA 运行动态门槛效应命令-基于 Seo M H 和 Kim S (2019)

3.1 命令的安装

ssc install xthenreg, replace //安装命令
net get xthenreg  //下载命令附带的数据到当前工作路径下

3.2 命令的语法

该命令的语法如下：

xthenreg depvar indepvars[if] [in][,endogenous(varlist) inst(varlist) kink staticgridnum (integer20) trimrate (real0.4) h0(real1.5) boost(real0)]

以下是注意事项：

xtset 应该进行在运行此命令之前，变量必须事先按照 (i) 面板变量和 (ii) 时间变量进行排序
需要高度平衡的面板数据
输入应为 y q x1 x2···，其中 q 为门槛变量, x1 x2 ··· 为其他自变量。
因为此命令使用 mm_quantilefunction，所以需要 moremata 外部命令包。
当存在内生自变量时应该设置内生选项。例如，如果 x1 是外生的，而 x2 是内生的，则输入必须为 y q x1，endo(x2).

各个主要选项的含义如下：

endogenous(varlist) 指定内生自变量。在逗号之前，必须将内生变量从自变量列表中排除。
inst(varlist) 指定其他工具变量的列表。
static 将模型设置为静态。与动态模型相反，静态模型不会自动包含 L.yas 独立变量。
Kink 设置模型 kink。
gridnum(integer) 确定网格点数以估计阈值 γ，默认值为 20。
trimrate(real) 在构造用于估计 r 的网格时确定修整率，默认值为 0.4。
h_0 (real) 用于内核估计的参数，默认值为 1.5。
boost(integer) 为线性测试的引导次数，默认值为 0。

3.3 命令的应用

首先，应用该方法用于评估肥胖对工人生产率的影响。肥胖是通过体重指数 (BMI) 来衡量的，体重指数 (公斤) 除以身高 (米) 平方。BMI 在 25 到 30 之间的个人被认为是超重的，且 BMI 为 30 或更高被视为肥胖。使用英国队列研究的数据和前面部分所述的方法，我们检查了 BMI 与工作时间的关系。有关更详细的讨论，请参见 Kim (2019)。在此示例中，我们使用以下模型结合 BMI 中的扭结来考虑男性工人的工作时间和 BMI。在这个例子中，我们使用下面的模型考虑男性工人的工作时间和身体质量指数。

y_{i t} = β_{0} + x_{i t} β_{1} + q_{i t} β_{2} + δ (q_{i t} - γ) I {q_{i t} \geq γ} + α_{i} + ε_{i t}

在这里，我们呈现出个体化的工作时间，例如：牙周病，家庭大小和体质指数的 BMI。我们有两个周期面板数据 (t = 1,2)，并采用如下的一阶差分法来去除与工作时间相关的各个时间不变特征 αi。

Δ y_{i 2} = Δ x_{i 2} β_{1} + Δ q_{i 2} β_{2} + δ (q_{i 2} - γ) 1 {q_{i 2} \geq γ} - (q_{i 1} - γ) δ 1 {q_{i 1} \geq γ} + Δ ε_{i 2}

为了实现 GMM 估计，我们在第一个差分模型中使用四个工具变量，出生体重 (bweight) 和工人自己的童年体重指数 (bmic) 和父母的体重指数 (bmim，bmid)，用于体重指数变量 qi2，qi2，qi1 在第一差分模型。xthenreg 的默认模型是一个动态模型。因为考虑静态模型，而不是动态模型，所以将使用静态选项，还通过使用扭结选项在模型中强加了一个扭结 (kink)。

表格前面的信息如下。N 是唯一受试者的总数，T 是时间段的数量 Number of moment conditions 取决于仪器的选择。在这个例子中，我们可以通过使用 exo 选项将所有外部变量收集到一个地方来获得相同的结果，如下所示：


xthenreg nour bmi, endo (bith)  exo (nsize) inst   ///
         (bwelght  bmic bmim bmid ) kink static

还可以使用 inst 选项来更改包含和排除工具的集合，根据力矩条件的数量相应地变化。

还可以在样本有限制的情况下估计模型，如下所示：

接下来，以下模型考察的是 BMI 效应的不连续性，而不考虑模型的扭曲。

y_{i t} = β_{0} + x_{i t} β_{1} + q_{i t} β_{2} + (δ_{0} + x_{i t} δ_{1} + q_{i t} δ_{2}) I {q_{i t} > γ} + α_{i} + ε_{i t}

通过取一阶差分，得到以下模型，并仅用统计方法进行估计。

Δ y_{i 2} = Δ x_{i 2} β_{1} + Δ q_{i 2} β_{2} + (δ_{0} + x_{i 2} δ_{1} + q_{i 2} δ_{2}) I {q_{i 2} > γ} - (δ_{0} + x_{i 1} δ_{1} + q_{i 1} δ_{2}) I {q_{i 1} > γ} + Δ ε_{i 2}

参考文献：

[1] Adam, C. S., and Bevan, D. L. (2005). “Fiscal deficits and growth in devel-oping countries,”Journal of Public Economics,89(4), 571-597.
[2] Giacomini, R., Politis, D. N., & White, H. (2013). “A warp-speed methodfor conducting Monte Carlo experiments involving bootstrap estimators,”Econometric theory, 29(3), 567-589.
[3] Hansen, B. E. (1999). “Threshold effects in non-dynamic panels: Estimation,testing, and inference,”Journal of econometrics, 93(2), 345-368.
[4] Hidalgo, J., Lee, J., and M.H. Seo (2019). “Robust Inference for ThresholdRegression Models,”Journal of Econometrics, to appear.
[5] Khan, M. S., and Ssnhadji, A. S. (2001). “Threshold effects in the relation-ship between inflation and growth,”IMF Staff papers, 48(1), 1-21.
[6] Kim, Y-J. (2019). “The effect of weight on work hours,” working paper.
[7] Seo, M. and Y. Shin (2016). “Dynamic panels with threshold effect andendogeneity,”Journal of Econometrics, 195: 169-186.
[8] Wang, Q. (2015). “Fixed-effect panel threshold model using Stata,”TheStata Journal, 15(1), 121-134.
[9] Zhang, Y., Zhou, Q., and Jiang, L. (2017). “Panel kink regression with anunknown threshold,”Economics Letters, 157, 116-121.

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata: 动态面板门槛模型

1、动态门槛效应回归模型背景简介

2、用 STATA 运行动态门槛效应命令-基于 Wang (2015) 方法

2.1 命令的安装示例数据导入

2.2 命令的语法

2.3 命令操作

2.3.1 拟合值估计

2.3.2代入拟合值估计门槛值和置信区间

3、用 STATA 运行动态门槛效应命令-基于 Seo M H 和 Kim S (2019)

3.1 命令的安装

3.2 命令的语法

3.2 命令的语法

3.3 命令的应用

参考文献：

相关课程

课程一览

关于我们

热门资讯