Stata 倍分法: 不满足平行趋势假设咋办？

发布时间：2020-04-17 阅读 14400

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

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

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

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

作者：黄欣怡 (中山大学) 邮箱：eugenehuangcheeks@163.com

「Source：Weights to address non-parallel trends in panel difference-in-differences models 」

1. 引言

双重差分法 ( Difference in differences，简称为 DID 或 DD ) 是用于估计处理效应的最为常见的计量方法，广泛应用于评估政策实施、项目执行及其他外生事件产生的效应。DID 法通过衡量处理组与控制组在被处理前后发生的变化，将处理组的变动减去控制组的变动，得到处理效应。由于在估计过程中采用了两次差分，故名为双重差分法。

采用 DID 法的重要前提假设为平行趋势假设，即假定处理组在未受到处理情况下的变化趋势与控制组相同。该假设在具体实证研究中转化为检验处理组与控制组在接受处理前 ( pre-treatment period ) 是否具有相同的变化趋势。关于平行趋势假设的检验在连享会之前的推文「多期 DID：平行趋势检验图示」和「Stata: 多期倍分法 (DID) 详解及其图示」中均有涉猎。

实证研究中，当个体接受处理的可能性或接受处理产生的效应与该个体随时间变动的因素相关时，可能不再满足平行趋势假设。在违背平行趋势假设的情况下，不能直接采用 DID 估计处理效应。部分研究采用了增加控制变量的方法，在回归中进一步控制了交互固定效应 ( interactive fixed effects ) 或随时间变化的混合因素 ( time-varing confounding factors ) 。

与此同时，部分学者通过改变样本的权重，使处理组与新构造的控制组在接受处理前具有相似的变化趋势。这类方法包括合成控制法 ( synthetic control method ) 、倾向得分匹配法 ( propensity score matching ) 、逆概率加权法 ( inverse probability weighting ) 等。

Ahlfeldt (2018) 指出，这类计量方法只适用于估计单一处理效应的情形；在需要区分多类相互关联的处理效应的研究情境中，或面对连续型 DID 的情形时，这些方法不再适用。Ahlfeldt (2018) 由此提出 Weighted Parallel Trends DD ( 以下简称 WPT DD ) 的估计方法。下文主要对 Ahlfeldt (2018) 提出的网格搜索算法 ( Grid Search ) 和迭代法 ( Iterative Approach ) 进行介绍。

2. 模拟数据的生成

2.1 双重差分模型

我们参照如下模型生成模拟数据：

y_{i t} = \sum_{n} α_{n} P_{t} D_{i}^{n} + μ_{i} + φ_{t} + ω_{i} f (t) + ϵ_{i t} (1)

其中， $y_{i t}$ 是个体 $i = 1, . . ., I$ 在时期 $t = 1, . . ., T$ 可观测到的结果变量， $D_{i}^{n} (n = 1, . . ., N)$ 代表所考察的 N 个处理组虚拟变量， $P_{t}$ 为时间虚拟变量 ( 当 $t \geq z$ 时， $P_{t} = 1$ ，否则取值为 0 ) 。 $μ_{i}$ 和 $φ_{t}$ 分别为个体固定效应与时间固定效应。 $ω_{i}$ 代表研究者难以观测的个体变化趋势， $f (t) = Σ_{0} Γ_{0} t^{0}$ 表示时间趋势，两者的交乘项 $ω_{i} f (t)$ 则表示因个体而异的时间趋势。

在模型 (1) 中， $α_{n} (n = 1, . . ., N)$ 是我们关注的 N 类处理效应。当 $ω_{i}$ 与 $D_{i}^{n}$ 不具有相关性时，遗漏 $ω_{i}$ 不会导致估计结果有偏；然而当 $ω_{i}$ 与 $D_{i}^{n}$ 相关时，遗漏变量将导致估计结果出现偏差。我们假设 $ω_{i}$ 与 $D_{i}^{n}$ 具有如下的相关性：

ω_{i} = ϑ_{i} + \frac{1}{N} H_{i} \sum_{n} D_{i}^{n} (2)

其中， $H_{i}$ 是我们难以观测的变量， $ϑ_{i}$ 是随机变量。此时，当我们在 DID 估计中遗漏 $ω_{i} f (t)$ 时，我们估计所得的处理效应为：

$E (\frac{\partial y_{i, t > z}}{\partial D_{i}^{n}} - \frac{\partial y_{i, t < z}}{\partial D_{i}^{n}}) = α_{n} + \frac{1}{N} E (H_{i}) ({\bar{f}}^{t \geq z} - {\bar{f}}^{t < z}) (3)$

其中， ${\bar{f}}^{t \geq z}$ 和 ${\bar{f}}^{t < z}$ 分别指样本在接受处理前期间 ( pre-treatment ) 与接受处理后期间 ( post-treatment ) 的平均时间趋势。可以看出，若 $E (H_{i}) \neq 0$ ，则处理效应的估计结果有偏。

遵循这一思路，Ahlfeldt (2018) 首先生成总体 $j = 1, . . ., J$ 。该总体满足 $E (H_{j}) = 0$ 。随后，依照抽样比率 $F_{j} = 1 / \sum_{m} r_{m} b_{j}^{m}$ 从总体中非随机地抽取样本 $i = 1, . . ., I$ 。我们假定 $F_{j}$ 与 $H_{j}$ 具有相关性，且 $H_{j} = - M / 2 + \sum_{m} b_{j}^{m}$ 。其中， $b_{j}^{m}$ 是不随时间变化的个体特征变量。令抽样权重 $S_{i} = F_{i}^{- 1}$ 。则在后续研究中，我们试图根据可观测样本 $i = 1, . . ., I$ 寻找 $S_{i}$ ，以纠正非随机抽样产生的偏差，从而使 $E (H_{j}) = 0$ 。

下文将利用样本接受处理前 ( pre-treatment period ) 的数据对权重 $S_{i}$ 进行搜索。若能够正确识别 $S_{i}$ ，则可消除 $Δ y_{j, p r e}$ 与 $D_{i}^{n}$ 的相关关系。由于 $b_{j}^{m}$ 是不随时间变化的特征变量，则 $H_{j}$ 亦不随时间变化，因此识别得到的 $S_{i}$ 也能够消除 $Δ y_{j, p o s t}$ 与 $D_{i}^{n}$ 的相关关系。

2.2 模拟数据生成过程

我们采用 Ahlfeldt (2018) 提供的 [数据和程序] 进行模拟。模拟数据由 1000 个体 ×10 年构成 ( $I = 1000, T = 10$ )，是含有 10000 个观测值的平行面板数据。

在模拟数据生成过程中，我们设定 $α_{n} = 1$ ，引入连续型处理变量 T1、T2，令样本接受处理的时点为 t = 5 。我们考察 2 类处理效应 $D^{n = 1}$ 和 $D^{n = 2}$ ( N=2 ) ；设定 M = 3，即 $H_{j}$ 包含 3 个个体特征变量： $b^{m = 1}, b^{m = 2}, b^{m = 3}$ 。我们仅就时间趋势服从一阶函数分布 ( $f (t) = t$ ) 的情形进行讨论。

* 导入模拟数据
cd "E:\Data\GA__PTW_REP"
u "DATA/INPUT/MCDATA.dta", replace

* 生成不随时间变化的个体特征变量，即 b1, b2, b3
qui gen h1 = runiform(0,1)
qui gen h2 = runiform(0,1)
qui gen h3 = runiform(0,1)

* 生成模型(2)中对应的变量 Hi
qui gen H =  (-1.5+h1+h2+h3)

* 生成两组 treatment variable (连续型 DID)
qui gen T1 = runiform(0,1)
qui gen T2 = runiform(0,1)

* 将 TREATMENT x H 加入至 TREND 变量中，从而生成模型 (2) 对应的个体变化趋势 OMEGA
qui replace TREND = TREND + (0.5 * T1 + 0.5 * T2) * H

* 生成抽样比率 Fi 与抽样权重 Si
* 生成 r1, r2, r3
gen RW = runiform(0,1)
* 生成抽样权重 Si
qui gen MODELW =  (RW[1]*h1+RW[2]*h2+(RW[3])*h3)
* 对抽样权重作标准化处理
qui sum MODELW
qui replace MODELW = MODELW /r(mean)
* 生成抽样权重的倒数：抽样比率 Fi
qui gen IVW = 1/MODELW

* 依据抽样比率 Fi 从总体 j = 1,...,J 中抽取样本 i = 1,...,I
qui gsample `obs' [w=IVW]

* 将数据由宽型数据转换为长型数据
        qui foreach num of numlist 1/10 {
        preserve
        keep ID LE YE Y`num' h*  TREND T1 T2 MODELW IVW
        gen PERIOD = `num'
        ren Y`num' Y
        gen panelID = _n
        save "TEMP/re_temp`num'.dta", replace
        restore
        }
        qui u "TEMP/re_temp1.dta", clear
        qui erase "TEMP/re_temp1.dta"
        qui foreach num of numlist 2/10 {
        append using "TEMP/re_temp`num'.dta"
        erase "TEMP/re_temp`num'.dta"
        }


* 定义面板数据
qui xtset panelID PERIOD

* 将f(T) x OMEGA 加入至结果变量中
qui replace Y = Y + TREND*PERIOD

* 将处理效应 TREATMENT 加入至结果变量中 ，以t=5为接受处理的时点
qui replace  Y = Y + T1 if PERIOD >=6
qui replace  Y = Y + T2 if PERIOD >=6

* 生成 TREATMENT x POST DID 变量，其系数即我们考察的处理效应
qui gen T1_POST = T1 * (PERIOD >=6)
qui gen T2_POST = T2 * (PERIOD >=6)

* 在结果变量中加入随机扰动项
qui replace Y = Y + rnormal(0,0.1)
* 存储数据以进行下一步运算
save "TEMP/temp", replace

3. DID 估计与 Stata 的实现

我们首先采用 OLS DD 估计模拟数据的处理效应。

** 估计 OLS DD
reg Y T1_POST T2_POST i.PERIOD, abs(panelID)
------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     T1_POST |   .6451472   .0786551     8.20   0.000     .4909652    .7993292
     T2_POST |   .3899474   .0813273     4.79   0.000     .2305274    .5493675
             |
      PERIOD |
          2  |   .2779982   .0528507     5.26   0.000     .1743988    .3815976
          3  |   .8783501   .0528507    16.62   0.000     .7747507    .9819495
          4  |   1.648978   .0528507    31.20   0.000     1.545378    1.752577
          5  |   1.991826   .0528507    37.69   0.000     1.888226    2.095425
          6  |   2.832213   .0769469    36.81   0.000      2.68138    2.983046
          7  |   2.538575   .0769469    32.99   0.000     2.387742    2.689409
          8  |   3.607691   .0769469    46.89   0.000     3.456857    3.758524
          9  |   3.472193   .0769469    45.12   0.000      3.32136    3.623027
         10  |   4.313866   .0769469    56.06   0.000     4.163033      4.4647
             |
       _cons |   1.315558   .0373711    35.20   0.000     1.242303    1.388814
------------------------------------------------------------------------------

其次，我们已知生成样本数据所采用的抽样比率 $F_{i}$ ，由此可通过 WPT DD 估计处理效应：

** 采用抽样权重估计 WPT DD
reg Y T1_POST T2_POST i.PERIOD [w=MODELW], abs(panelID)

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     T1_POST |   .8356594   .0766697    10.90   0.000     .6853693    .9859495
     T2_POST |     1.0463   .0769499    13.60   0.000     .8954609    1.197139
             |
      PERIOD |
          2  |   .3825006   .0507956     7.53   0.000     .2829295    .4820716
          3  |   1.083468   .0507956    21.33   0.000     .9838966    1.183039
          4  |   1.959798   .0507956    38.58   0.000     1.860227    2.059369
          5  |   2.404149   .0507956    47.33   0.000     2.304578     2.50372
          6  |   2.931054   .0742346    39.48   0.000     2.785538    3.076571
          7  |   2.739205   .0742346    36.90   0.000     2.593688    2.884722
          8  |   3.910568   .0742346    52.68   0.000     3.765051    4.056084
          9  |   3.879056   .0742346    52.25   0.000     3.733539    4.024573
         10  |   4.824459   .0742346    64.99   0.000     4.678942    4.969976
             |
       _cons |   1.420531   .0359179    39.55   0.000     1.350124    1.490938
------------------------------------------------------------------------------

下文我们将分别通过网格搜索法 ( Grid Search ) 与迭代法 ( Iterative Approach ) 搜索 $S_{i}$ ，进而对处理效应作 WPT DD 估计。

3.1 网格搜索 ( Grid Search )

我们假定抽样权重 $\hat{S_{i}} = \sum_{m} q_{m} b_{m}^{i}$ ，并且研究者已知 $b_{i}^{m}$ ，因此我们尝试通过网格搜索的方法确定 $Q = (q_{m = 1}, q_{m = 2}, q_{m = 3})$ 。若将步长设定为 0.1，则搜索范围为 $q_{m} = 0, 0.1, 0.2, \dots, 1 (m = 1, 2, 3)$ 。在迭代搜索的过程中将采用如下模型估计：

Δ {\dot{y}}_{i 1} = c_{Q}^{0} + \sum_{n} c_{Q}^{n} {\dot{D}}_{i}^{n} + η_{Q i} (4)

$Δ y_{i 1}$ 是结果变量自第 1 期至第 2 期的变动幅度，变量上加点表明该变量经标准化处理； $c_{Q}^{n}$ 衡量了样本在该期间的标准化处理效应。我们试图寻找参数 $Q$ 使 $Δ y_{i, p r e}$ 与 $D_{i}^{n}$ 的相关性最小化，即使 $c_{Q}^{n}$ 最小化。

可采用的目标函数共分为 3 类，分别为：

叠加法 (additive approach)：最小化目标方程 $B^{A} = \sum_{n} ({\hat{c}}_{Q}^{n})^{2}$ 。
连乘法 (multiplicative approach)：最小化目标方程 $B^{m} = \prod_{n} ({\hat{c}}_{Q}^{n})^{2}$ 。
极值法 (min-max approach)：最小化目标方程 $B^{M} = m a x (| {\hat{c}}_{Q}^{n} |)$ 。

在网格搜索中，我们选择叠加法进行演算，即最小化目标方程 $B^{A} = \sum_{n} ({\hat{c}}_{Q}^{n})^{2}$ 。

* 导入数据并调整数据格式
qui u "TEMP/temp", clear
    keep if PERIOD ==2
    keep ID Y
    ren Y Y2
    save "TEMP/ptwtemp", replace
qui u "TEMP/temp", clear
    keep if PERIOD ==1
    merge m:m ID using "TEMP/ptwtemp.dta"
    drop _m
    ren Y Y1
    erase "TEMP/ptwtemp.dta"

* 生成结果变量自第1期至第2期的变化值
    gen DELTA = Y2- Y1

* 将个体特征变量标准化
    foreach num of numlist 1/3 {
    gen VAR`num'= h`num'
    sum VAR`num',d
    gen wa_VAR`num' = VAR`num'- r(min)
    replace wa_VAR`num' = wa_VAR`num'/r(sd)
    }

* 将TREATMENT VARIABLE标准化
    sum T1
    replace T1 = T1/ r(sd)
    sum T2
    replace T2 = T2 / r(sd)

* 生成网格搜索法中需要用到的暂时性变量
    gen wa = .
    gen OBJ = .
    gen OBJmin = .
    gen VAR1min = .
    gen VAR2min = .
    gen VAR3min = .
    gen wa_min = .

* 确定搜索范围与搜索步长
    local min = 0
    local max = 1
    local step = 0.1

* 运行网格搜索

* 开启对q1的循环
local om_VAR1 = `min'
while `om_VAR1' <= `max' {
    * 开启对q2的循环
    local om_VAR2 = `min'
    while `om_VAR2' <= `max' {
        * 开启对q3的循环
        local om_VAR3 = `min'
        while `om_VAR3' <= `max' {

            * 生成权重Si
            qui replace wa = `om_VAR1'*wa_VAR1+`om_VAR2'*wa_VAR2+`om_VAR3'*wa_VAR3
            * 确保抽样权重Si>0
            qui sum wa
            if r(mean) > 0 {
            * 对Si作标准化处理
            qui replace wa = wa / r(mean)
            * 依照模型(4)进行加权回归估计处理效应
            qui reg DELTA T1 T2 [iw=wa]
            * 计算目标方程
            qui replace OBJ = _b[T1]^2+_b[T2]^2
            * 若 qm 使目标方程较上一轮循环更小，则存储参数
            qui replace VAR1min = `om_VAR1' if OBJ < OBJmin
            qui replace VAR2min = `om_VAR2' if OBJ < OBJmin
            qui replace VAR3min = `om_VAR3' if OBJ < OBJmin
            * 若目标方程较上一轮循环更小，存储计算所得的 WPT 权重 Si
            qui replace wa_min = wa if OBJ < OBJmin
            qui replace OBJmin  = OBJ if OBJ < OBJmin
            display "VAR1 "`om_VAR1' " VAR2 " `om_VAR2' " 'VAR3 " `om_VAR3'
                    }
                    else {
                    }
            * 结束对q3的循环
            local om_VAR3 = `om_VAR3'+`step'
            }
        * 结束对q2的循环
        local om_VAR2 = `om_VAR2'+`step'
        }
    * 结束对q1的循环
    local om_VAR1 = `om_VAR1'+`step'
    }

* 结束循环并存储样本权重
    ren wa_min WA_MC
    keep ID WA_MC
    save "TEMP/WA_MC.dta", replace

* 将计算所得的样本权重与原模拟数据进行合并
    qui u "TEMP/temp", clear
    qui merge m:m ID using "TEMP/WA_MC.dta"
    erase "TEMP/WA_MC.dta"
    qui tab _m
    qui drop _m

* 估计 WPT DD
    reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)

* 去除权重变量
    qui drop WA_MC

我们得到 WPT DD 的估计结果如下：

reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     T1_POST |   .8385112    .077283    10.85   0.000     .6870188    .9900035
     T2_POST |   1.090782   .0778193    14.02   0.000     .9382385    1.243326
             |
      PERIOD |
          2  |   .3796894   .0510905     7.43   0.000     .2795404    .4798385
          3  |   1.079326   .0510905    21.13   0.000     .9791772    1.179475
          4  |   1.952109   .0510905    38.21   0.000      1.85196    2.052258
          5  |   2.392124   .0510905    46.82   0.000     2.291975    2.492273
          6  |   2.893577   .0748875    38.64   0.000      2.74678    3.040373
          7  |    2.69863   .0748875    36.04   0.000     2.551833    2.845426
          8  |   3.865206   .0748875    51.61   0.000     3.718409    4.012002
          9  |   3.832881   .0748875    51.18   0.000     3.686084    3.979677
         10  |    4.77509   .0748875    63.76   0.000     4.628293    4.921886
             |
       _cons |   1.417661   .0361264    39.24   0.000     1.346845    1.488477
------------------------------------------------------------------------------

3.2 迭代法 ( Iterative Approach )

迭代法中，我们以 ${\hat{S_{i}}}^{s = 1} = 1$ 作为迭代循环的起点 ( s = 1 表示第一轮循环)。在每一轮循环中首先估计如下模型：

Δ {\dot{y}}_{i 1} = c_{s}^{0} + \sum_{n} c_{s}^{n} {\dot{D}}_{i}^{n} + Ψ_{s i} (5)

其中 $c_{s}^{n}$ 是我们关注的处理效应。随后在模型 (5) 中加入处理变量与个体特征变量的交乘项 ${\dot{D}}_{i}^{n} \times b_{i}^{m}$ 进行估计：

Δ {\dot{y}}_{i 1} = b_{s}^{0} + \sum_{n} b_{s}^{n} {\dot{D}}_{i}^{n} + \sum_{m} \sum_{n} b_{s}^{n, m} {\dot{D}}_{i}^{n} \times b_{i}^{m} + ξ_{s i} (6)

在模型 (6) 中，估计得到的边际处理效应为

{\hat{π}}_{i}^{n} = \partial {\dot{y}}_{i 1} / \partial {\dot{D}}_{i}^{n} = {\hat{b}}_{s}^{n} + \sum_{m} \sum_{n} {\hat{b}}_{s}^{n, m} \times b_{i}^{m}

迭代的具体思路为：每一轮循环对权重 ${\hat{S}}_{i}^{s}$ 的计算均基于上一轮循环所得权重 ${\hat{S}}_{i}^{s - 1}$ 以及本轮循环估计模型 (6) 得到的边际处理效应 ${\hat{π}}_{i}^{n}$ 进行调整，即 ${\hat{S}}_{i}^{s, n} = g ({\hat{S}}_{i}^{s - 1, n}, {\hat{π}}_{i}^{n})$ 。

当模型 (5) 估计所得 $c_{s}^{n} < 0$ 时，我们在 ${\hat{S}}_{i}^{s - 1}$ 的基础上赋予 ${\hat{π}}_{i}^{n} > 0$ 的个体更高权重，此时 $g_{S} > 0, g_{π} > 0$ ；当模型 (5) 估计所得 $c_{s}^{n} > 0$ 时，我们在 ${\hat{S}}_{i}^{s - 1}$ 的基础上赋予 ${\hat{π}}_{i}^{n} < 0$ 的个体更高权重，此时 $g_{S} > 0, g_{π} < 0$ 。如此一来，经过多次迭代后可使得 $∣ {\hat{c}}_{s}^{n} ∣$ 最小化。

迭代法中我们选取目标方程是 $B^{M} = m a x (| {\hat{c}}^{1} |, | {\hat{c}}_{s}^{2} |)$ ，当目标方程函数值小于我们设定的门槛值 0.005 时，即 $B^{M} < 0.005$ 时，循环结束。

* 存储数据
save "TEMP/temp", replace

* 清除内存
set more off
estimates drop _all
program drop _all

* 设定迭代次数为100
local MAX = 100
* 设定目标方程的目标值为0.005
local TARGET = 0.005

* 导入数据并整理数据格式
u "TEMP/temp.dta", clear
keep if PERIOD ==2
keep ID Y
ren Y Y2
save "TEMP/ptwtemp.dta", replace
u "TEMP/temp.dta", clear
keep if PERIOD ==1
merge m:m ID using "TEMP/ptwtemp.dta"
drop _m
ren Y Y1
erase "TEMP/ptwtemp.dta"

* 生成结果变量自第1期至第2期的变化值
    gen DELTA = Y2- Y1

* 对结果变量和处理变量作标准化处理
    foreach var of varlist DELTA T1 T2 {
        qui sum `var'
        qui gen S_`var' = `var'/r(sd)
        }
* 生成处理变量T1 T2与个体特征变量bm的交乘项
    qui foreach var1 of varlist S_T1 S_T2 {
        foreach var2 of varlist h1 h2 h3 {
        gen `var1'_`var2' = `var1'*`var2'
        }
        }
* 生成迭代所需要的暂时性变量
    gen W = 1
    gen W1 = 1
    gen W2 = 1
    gen WOLD = 1
    qui gen ME_T1 = .
    qui gen ME_T2 = .
    gen WBEST = 1
    gen bBEST = 100
    scalar TH = 0.00
    gen R = .

* 定义迭代法所使用的程序 ************************************************
    program ALGO
        * 估计模型(5)
            qui replace WOLD = W
            qui reg S_DELTA S_T1 S_T2 [w=W]
        * 存储估计所得的处理效应 (模型(5)中的c_s^n)
            qui scalar bT1 = _b[S_T1]
            qui scalar bT2 = _b[S_T2]
        * 存储目标方程的函数值
            qui scalar bsum = abs(bT1)+abs(bT2)

        * 估计模型(6)
            qui reg S_DELTA S_T1 S_T2 S_T1_h* S_T2_h* [w=W]
        * 存储估计所得的处理效应 (模型(6)中的b_s^n,m)
                qui replace ME_T1 = _b[S_T1]
                qui replace ME_T2 = _b[S_T2]
                qui foreach num of numlist 1/3 {
                qui     replace ME_T1 = ME_T1 + _b[S_T1_h`num']*h`num'
                qui     replace ME_T2 = ME_T2 + _b[S_T2_h`num']*h`num'
                }

        * 生成调整后的权重Si

            * TREATMENT 1 (第一类处理效应)
                * 若模型(5)中估计所得的处理效应为负
                qui if bT1 < -1*TH {
                display "neg"
                * 增加处理效应为正的个体在样本中的权重
                replace W1 = W + W*ME_T1 if ME_T1 > 0
                }
                else {
                * 若模型(5)中估计所得的处理效应为正
                if bT1 > TH {
                        display "pos"
                 * 增加处理效应为负的个体在样本中的权重
                         replace W1 = W - W*ME_T1 if ME_T1 < 0
                        }
                        else {
                        }
                    }
                        else {
                        }

            * TREATMENT 2 (第二类处理效应)
                * 若模型(5)中估计所得的处理效应为负
                qui if bT2 < -1*TH {
                display "neg"
                * 增加处理效应为正的个体在样本中的权重
                replace W2 = W + W*ME_T2 if ME_T2 > 0
                }
                else {
                * 若模型(5)中估计所得的处理效应为正
                if bT2 > TH {
                        display "pos"
                 * 增加处理效应为负的个体在样本中的权重
                        replace W2 = W - W*ME_T2 if ME_T2 < 0
                        }
                        else {
                        }
                    }
                        else {
                        }

        * 根据两类处理效应的估计结果对权重进行调整
        * 对于绝对值较大的处理效应给予更大幅度的调整
            qui replace W =[exp(abs(bT1))* W1+exp(abs(bT2))*W2] / [exp(abs(bT1))+exp(abs(bT2))]
        * 基于新生成的权重估计模型(5)
            qui reg S_DELTA S_T1 S_T2 [w=W]
        * 将处理效应存储为单值
            qui scalar bT1new = _b[S_T1]
            qui scalar bT2new = _b[S_T2]
        * 计算新的目标方程函数值
            qui scalar bsumnew = abs(_b[S_T1])+abs(_b[S_T2])
        * 基于目标方程的函数值判定是否对权重进行更新
            qui replace bBEST = bsumnew if  bsumnew<bBEST
            qui replace WBEST = W if  bsumnew<bBEST

        * 若上一轮循环得到的权重更优，则结合本次循环和上一轮循环得到的权重进行调整
            if bsumnew > bBEST + 0.025 {
                replace W = 0.8*WBEST+0.2*W
                }
                else {
                }

    end
* 程序ALGO定义结束    ****************************************************


* 迭代循环过程 *************************************************************
    local it = 1

    while `it' < `MAX'  {
        * 运行以上定义的程序 ALGO
        qui ALGO
        * 基于本轮循环得到的权重估计处理效应
        qui reg S_DELTA S_T1 S_T2 [w=W]
        * 存储估计所得处理效应
        qui local T1CORR = _b[S_T1]
        qui local T2CORR = _b[S_T2]

        * 展示迭代结果
        display "Iteration "`it' " T1 CORR " `T1CORR' " T2 CORR " `T2CORR'

        * 检验是否达到目标方程函数值小于0.005的目标
        if abs(`T1CORR') < `TARGET' & abs(`T2CORR') < `TARGET' {
            local it = 1000
            }
            else {
            local it = `it'+1
            }
        }
* 存储最终得到的pre-treatment period的处理效应
    reg  S_DELTA S_T1 S_T2 [w=W]
    scalar T1CORR = _b[S_T1]
    scalar T2CORR = _b[S_T2]

* 存储最终得到的最优权重
    qui sum W
    gen WA_MC = W / r(mean)
    keep ID WA_MC
    save "TEMP/WA_MC.dta", replace

* 将原模拟数据与权重数据合并
    qui u "TEMP/temp", clear
    qui merge m:m ID using "TEMP/WA_MC.dta"
    erase  "TEMP/WA_MC.dta"
    qui tab _m
    qui drop _m
* 估计 WPT DD
    reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)

* 去除权重变量
    qui drop WA_MC
    erase "TEMP/temp.dta"

我们得到 WPT DD 的估计结果如下：

reg Y T1_POST T2_POST i.PERIOD [w=WA_MC], abs(panelID)

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     T1_POST |   .8385112    .077283    10.85   0.000     .6870188    .9900035
     T2_POST |   1.090782   .0778193    14.02   0.000     .9382385    1.243326
             |
      PERIOD |
          2  |   .3796894   .0510905     7.43   0.000     .2795404    .4798385
          3  |   1.079326   .0510905    21.13   0.000     .9791772    1.179475
          4  |   1.952109   .0510905    38.21   0.000      1.85196    2.052258
          5  |   2.392124   .0510905    46.82   0.000     2.291975    2.492273
          6  |   2.893577   .0748875    38.64   0.000      2.74678    3.040373
          7  |    2.69863   .0748875    36.04   0.000     2.551833    2.845426
          8  |   3.865206   .0748875    51.61   0.000     3.718409    4.012002
          9  |   3.832881   .0748875    51.18   0.000     3.686084    3.979677
         10  |    4.77509   .0748875    63.76   0.000     4.628293    4.921886
             |
       _cons |   1.417661   .0361264    39.24   0.000     1.346845    1.488477
------------------------------------------------------------------------------

4. WPT DD 的估计结果解读

Ahlfeldt (2018) 进行了总计 1000 次蒙特卡洛模拟实验，图 1 呈现了蒙特卡洛模拟计算所得的处理效应 T1 和 T2 的分布，表 1 展示了处理效应分布的均值、中位数和标准差。

由图表可以看出，OLS DD 估计所得处理效应明显下偏，低于我们实际设定的处理效应 $α_{n} = 1$ 。通过网格搜索 ( Grid Search ) 与迭代法 ( Iterative Approach ) 调整所得的处理效应则集中于我们的设定值 1 附近，经过多次模拟计算所得的处理效应均值 ${\hat{α}}_{n}$ 亦更接近于 1。与网格搜索 ( Grid Search ) 相比，迭代法 ( Iterative Approach ) 对应的估计结果更加离散，但估计准确度仍然显著优于 OLS DD 。

由此可以看出，WPT DD 有利于改善因违背平行趋势假设导致 DID 估计结果有偏的问题。

参考文献

Ahlfeldt, G. M., 2018, Weights to address non-parallel trends in panel difference-in-differences models, CESifo Economic Studies, 64 (2): 216-240. [PDF]
多期 DID：平行趋势检验图示
Stata: 多期倍分法 (DID) 详解及其图示

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