Stata：纠偏合成控制法介绍-allsynth

发布时间：2022-01-19 阅读 3232

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

作者：伍站锋 (中山大学)
邮箱：wuzhf25@mail2.sysu.edu.cn

编者按：本文主要参考自下文，特此致谢！
Source：Abadie A, L’Hour J. A penalized synthetic control estimator for disaggregated data[J]. Journal of the American Statistical Association, 2021: 1-18. -PDF-

与 PSM-DID 相似，合成控制法也是用于反事实因果推断。它的应用场景是某项政策只针对部分地区 (称之为处理组) ，而其他地区并未受到政策影响 (称之为控制组)。文献中常见的做法是使用倾向性匹配得分 (PSM) 的方法，选取协变量消除处理组和控制组在其他诸多方面的系统性差异。

但是，在很多情况下，因为各个地区之间的经济、政治、文化等方面的差异本身就很大，即使在使用 PSM 方法之后，配对的地区在政策出台前依然有较大差异。因此，在这种场景下，使用合成控制法更为合适。

1. 加州控烟的政策效应

Abadie 和 Gardeazabal (2003) 最早在研究西班牙巴斯克地区恐怖袭击的经济后果时，使用了单变量的合成控制法。随后，Abadie 等 (2010) 在研究 1989 年开始实施的加州控烟法案对香烟销量的影响时，进一步提出了有多个协变量的合成控制法：即对控制组赋予一个权重矩阵，根据权重矩阵形成的线性组合形成一个 “合成加州”，而实际的加州称为 “真实加州”。

如果在政策实施前，“合成加州” 和 “真实加州” 的香烟销量能在统计意义上相等，而在政策实施之后，“真实加州” 的香烟销量在统计意义上比 “合成加州” 显著下降，那么可以认为加州控烟法案取得了效果。

接下来，我们将简单介绍合成控制法的估计原理以及它存在的问题，并且给出解决这一问题的方法，即纠偏合成控制法。

2. 传统合成控制法介绍

2.1 传统合成控制法的估计方法

在加州控烟问题中，我们可以用下面的公式表示 “合成加州” 的香烟销量：

{\hat{Y}}_{i, t}^{N} = \sum_{j = n_{1} + 1}^{n} {\hat{W}}_{i, j} \times Y_{j, t}

一般来说在 $n$ 个个体中，我们认为有 $n_{1}$ 个处理组个体， $(n - n_{1})$ 个控制组个体。具体到加州控烟案例中，等式左边的是反事实估计下的加州香烟销量，下标 $i$ 表示的是处理组中的个体 (在上述式子中是合成的)，等式右边的下标 $j$ 表示的是控制组中的个体，在加州控烟案例中是美国其他 38 个州。基于此，我们可以估计出某个时点 $t$ 上的平均处理效应：

{\hat{γ}}_{t} = \frac{1}{n_{1}} \sum_{i = 1}^{n_{1}} (Y_{i, t} - \sum_{j = n_{1} + 1}^{n} {\hat{W}}_{i, j} \times Y_{j, t})

如果我们令 $T$ 代表样本观测区间， $T_{0}$ 则指的是政策开始发生的时间，时间点 $t$ 满足 $T_{0} < t < T$ 。

整个区间上的 $\hat{γ}$ (包括政策前和政策后，我们称之为 "gap")，描述了在整个研究的时间段上处理组个体和合成个体的差异。换而言之，"gap" 在政策实施后时间段的数据，表示政策处理效应。"gap" 在整个研究区间上的时间趋势图正是合成控制法最关心的数据变化。

合成控制法的核心在于计算出权重矩阵 $W$ 的估计值 $\hat{W}$ 。合成控制法采用的是最小均方误差的估计办法，也就是我们希望在政策发生前，处理组和控制组的差异 "gap" 尽可能小。具体的数学表达式如下：

m i n ∥ X_{i, t} - \sum_{j = n_{1} + 1}^{n} W_{i, j} \times X_{j, t} ∥^{2}

s . t . W_{i, j} \geq 0,

\sum_{j = n_{1} + 1}^{n} W_{i, j} = 1

其中， $t < T_{0}$ 。矩阵 $X$ 包括了回归方程中的协变量 (控制变量) 和被解释变量，最优化条件可以表达为：寻找一个权重矩阵 $\hat{W}$ ，使得合成加州和真实加州在政策发生前的各项指标最为接近。

以上过程可以用 Stata 命令 synth 来实现。synth 命令的基本语法格式如下：

synth depvar predictorvars(x1 x2 x3),
      trunit(#) trperiod(#)
      [counit(numlist) xperiod(numlist) mspeperiod()
      resultsperiod() nested allopt unitnames(varname)
      figure keep(file) customV(numlist) optsettings ]

具体到这一案例中，加州控烟案例研究的代码如下：

lxhuse smoking.dta, clear
xtset state year
synth cigsale retprice lnincome age15to24 beer      ///
    cigsale(1975) cigsale(1980) cigsale(1988),      ///
    trunit(3) trperiod(1989) xperiod(1980(1)1988)   ///
    figure nested keep(smoking_synth)

其中，

必选项 trunit(#) 用于指定处理地区。在原有数据集中加州的对应数字是 3，因此代码为 trunit(3)；
必选项 trperiod(#) 用于指定政策干预开始的时期。加州控烟案例的政策开始时间是 1989 年；
选择项 counit(numlist) 用于指定潜在的控制地区，此处默认为数据中的所有控制组；
选择项 xperiod(numlist) 用于指定将预测变量进行平均的期间，默认为政策干预开始之前的所有时期；
选择项 mspeperiod() 用于指定最小化均方预测误差 (MSPE) 的时期，默认为政策干预开始之前的所有时期；
选择项 figure 绘制 “真实加州” 与 “合成加州” 的结果变量时间趋势图；
选择项 nested 表示使用嵌套的数值方法寻找最优的合成控制 (推荐使用此选项)，这比默认方法更费时间，但可能更精确。如果使用选择项 nested allopt，则比单独使用 nested 耗时更长，但精确度可能更高。
选择项 keep(filename) 将估计结果 (比如，合成控制的权重和结果变量) 存为 Stata 数据集 (smoking_synth.dta)，以便进行后续计算。

上述代码的运行结果如下：

Unit Weights:
----------------------------
         Co_No | Unit_Weight
---------------+------------
       Alabama |           0
      Arkansas |           0
      Colorado |        .146
   Connecticut |        .062
      Delaware |           0
       (output omitted)
     Wisconsin |           0
       Wyoming |           0
----------------------------

Predictor Balance:
------------------------------------------------------
                               |   Treated  Synthetic 
-------------------------------+----------------------
                      retprice |  89.42222   89.42084 
                      lnincome |  10.07656   9.853375 
                     age15to24 |  .1735324   .1735971 
                          beer |     24.28   24.22762 
                 cigsale(1975) |     127.1   127.0591 
                 cigsale(1980) |     120.2   120.2444 
                 cigsale(1988) |      90.1    91.4341 
------------------------------------------------------

事实上，合成控制法下大部分控制组的权重为 0，这一特征我们称之为：稀疏性。在论文写作中则往往需要呈现一张 "gap" 变化趋势图，用来直观地反映政策效应。在 synth 命令下，我们需要手动生成这张图片，具体命令如下：

use smoking.dta,clear
xtset state year

// 假设每个州都是处理组, 逐个遍历进行合成控制
forval i=1/39{
    qui synth cigsale beer(1984(1)1988) lnincome retprice age15to24 cigsale(1975) ///
    cigsale(1980) cigsale(1988), trunit (`i') trperiod(1989) bcorrect(replace) keep(synth_`i') rep
}
// 逐一生成每个州作为处理组时的 gap 文件, 存储为 dta 数据集
forval i=1/39{
    use synth_`i', clear
    rename _time years
    gen gap_`i' = _Y_treated - _Y_synthetic
    keep years gap_`i'
    drop if missing(years)
    save synth_`i', replace
}
// 合并上述的 dta 文件
use synth_1, clear
	forval i=2/39{
	qui merge 1:1 years using synth_`i', nogenerate
}
save synth.dta,replace
// 生成一个暂元, 存储图像绘制
use synth.dta,replace
local lp1
forval i=1/39 {
   local lp1 `lp1' line gap_`i' years, lpattern(dash) lcolor(gs8) ||
}
// 最终展示
twoway `lp1' || line gap_3 years, scale(1.3) lcolor(black) legend(pos(6) ///
    order(1 "其他控制组" 40 "加州")) xline(1989, lpattern(dash))

从图 1 中，我们可以观察到在政策发生之前，真实加州与合成加州的香烟销量差距很小，即在 gap=0 附近上下波动，这证明根据以上合成控制法合成的 “加州” 效果较好；而在政策发生之后，gap 逐渐下降，表明真实加州的香烟销量比合成加州的预测香烟销量更少，加州控烟法案或许有一定成效。

2.2 传统合成控制法存在的问题

然而，上述过程看似严谨，但是从代数的角度来剖析，待估的权重矩阵 $\hat{W}$ 存在解不唯一的问题。我们将加州控烟案例暂时简化为只有香烟销量数据的单变量问题。具体解释如下 ( $\dots$ 表示省略矩阵中具体的元素)：

设控制组矩阵 X_{0, 19 \times 38} = [\begin{array}{lll} s a l e s_{2, 1970} & s a l e s_{3, 1970} & s a l e s_{4, 1970} & \dots \\ s a l e s_{2, 1971} & s a l e s_{3, 1971} & s a l e s_{4, 1971} & \dots \\ s a l e s_{2, 1972} & s a l e s_{3, 1972} & s a l e s_{4, 1972} & \dots \\ \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots \end{array}]

设待估的 “ 合成加州 ” 向量 {\hat{X}}_{1, 19 \times 1} = [\begin{array}{lll} {\hat{s a l e s}}_{1, 1970} \\ {\hat{s a l e s}}_{1, 1971} \\ {\hat{s a l e s}}_{1, 1972} \\ \dots \\ \dots \end{array}]

当然，待估的 ${\hat{X}}_{1, 19 \times 1}$ 还需要满足以下条件：

m i n ∥ X_{1, 19 \times 1} - {\hat{X}}_{1, 19 \times 1} ∥^{2}, X_{1, 19 \times 1} 是真实加州向量

其中矩阵 $X_{0}$ 描述了加州以外 38 个州 (用下标 2、3、4 等来表示) 从 1970-1988 年这 19 年间 (即控烟法案出台之前) 历年的香烟销量，矩阵 $X_{1}$ 描述了法案出台前加州历年的香烟销量。在实际研究中，矩阵 $X_{0}$ 和矩阵 $X_{1}$ 还包括其他协变量。所以，最关键的在于推导出权重矩阵的估计值 $\hat{W}$ ：

X_{0, 19 \times 38} \times {\hat{W}}_{38 \times 1} = {\hat{X}}_{1, 19 \times 1}

由于每个州的时间序列之间很难存在线性相关关系，根据线性代数的知识，这就会导致系数矩阵 $X_{0}$ 不是满秩矩阵，造成线性方程的解不唯一。尤其是当控制组内的个体数量非常多时，解不唯一的问题会更加严重。如果权重矩阵解不唯一，那么基于合成控制法反事实因果推断在实操上会不严谨，甚至无法得到一个唯一的 “合成加州”，最后得到的结论就没有说服力。

纠偏合成控制法的重要目的就是解决这一的问题。接下来，将主要以图像展示的方式来解释 Abadie 和 L'Hour (2021) 提出的方法，以帮助读者更直观地理解纠偏合成控制法。

3. 纠偏合成控制法介绍

3.1 纠偏合成控制法的估计方法

Abadie 和 L'Hour (2021) 提出的方法就是在最小均方误差估计中添加一个 “惩罚项” (Penalized)，构造一个 “带惩罚项的合成控制” (penalized version of the synthetic control)。具体方法如下：

m i n {∥ X_{i, t} - \sum_{j = n_{1} + 1}^{n} W_{i, j} X_{j, t} ∥}^{2} + λ \sum_{j = n_{1} + 1}^{n} W_{i, j} {∥ X_{i, t} - X_{j, t} ∥}^{2}

s . t . W_{i, j} \geq 0,

\sum_{j = n_{1} + 1}^{n} W_{i, j} = 1

其中， $t < T_{0}$ ， $λ > 0$ 。最小化条件的右侧多项式就是惩罚项，我们暂时称之为 “惩罚项多项式”。

那么这里的 “惩罚项” 又有何含义呢？我们首先来理解矩阵 (或者说矩阵中的向量) 的含义。以加州控烟案例为例子，矩阵中的每一列向量包含了每个州的不同方面的信息 (香烟销量、啤酒销量、州内青少年人口等等)， $∥ X_{i, t} - X_{j, t} ∥^{2}$ 其实衡量了处理组 (加州) 和控制组的信息差异。

合成控制法的目的是合成一个在政策实施前和 “真实加州” 各个方面非常接近的 “合成加州”。所以，我们希望用于合成的控制组和真实加州 (处理组) 在各个维度上的信息含量比较接近。因此，对于 $∥ X_{i, t} - X_{j, t} ∥^{2}$ 的数值较大的个体，我们应当赋予它们较低的权重；反之，对于 $∥ X_{i, t} - X_{j, t} ∥^{2}$ 的数值较小的个体，我们应当赋予它们较高的权重。

所以，我们不妨先把 “惩罚项多项式” 中 $∥ X_{i, t} - X_{j, t} ∥^{2}$ 数值较大项的权重赋为 0。在加州控烟案例中，也就是把和加州在各个维度上的信息差异最大的州先赋予 0 的权重，在剩下的 38 个州里面重复上述步骤，直到系数矩阵 $X_{0}$ 成为为满秩矩阵为止。

细心的读者会提出疑问：那这样一来，很多控制组个体的权重不就是 0 了吗？这不就损失了数据的信息含量吗？其实我们不妨回忆在本文第一部分提到的合成控制法的重要特征 “稀疏性”。事实上，合成控制法的确会出现很多权重为 0 的情况，这是合成控制法得到的结果的基本特征。因此，“惩罚项” 方法没有违背合成控制法的基本思想，在操作中也符合合成控制法的特征。

另一方面，“合成加州” 不仅仅需要在政策发生前和 “真实加州” 的差异很小，还需要对政策发生后的的数据变化有较好地预测。如果 “合成加州” 将所有控制组个体的信息都包含在内，那么预测中就包含了许多不重要的信息，会产生过拟合问题。惩罚项方法在机器学习领域经常使用，目的也是为了解决过拟合问题。纠偏合成控制法的 “惩罚项” 也有这一作用。

本文特地使用图像来更简洁地表达这一过程：

图 2 是以 4 个控制组个体、1 个处理组个体的情况为例子。4 个绿色的向量 $A_{1}^{'}$ $A_{2}^{'}$ $A_{3}^{'}$ $A_{4}^{'}$ 代表了四个控制组个体所包含的信息，这四个向量构成了一个平面 $l$ ；蓝色向量 $B^{'}$ 代表处理组个体所包含的信息。 $∥ X_{i, t} - \sum_{j = n_{1} + 1}^{n} W_{i, j} \times X_{j, t} ∥$ 的几何含义就是向量 $B^{'}$ 到向量 $A^{'}$ (包含了 $A_{1}^{'}$ $A_{2}^{'}$ $A_{3}^{'}$ $A_{4}^{'}$ ) 合成的向量 $C^{'}$ 的距离 $d$ ，在图中表示为紫色虚线 $d$ 。

而最优化条件 $m i n ∥ X_{i, t} - \sum_{j = n_{1} + 1}^{n} W_{i, j} \times X_{j, t} ∥^{2}$ 的几何含义就是使得距离 $d$ 最小化。绿色向量上的红色部分是每个向量合成向量 $C^{'}$ 时提供的权重，他们可以构成权重矩阵(向量)的估计值 $\hat{W}$ 。

通过图 3 与图 2 对比，我们可以观察到，红色部分的权重向量在数值上进行调整后，依然可以合成和原来相同的向量 $C^{'}$ 。比如，在图 2 中，我们赋予了向量 $A_{4}^{'}$ 较高的权重；但是在调整过后，我们赋予向量 $A_{4}^{'}$ 更低的权重，同时赋予向量 $A_{3}^{'}$ 和向量 $A_{1}^{'}$ 更高的权重，但最终都可以合成同一个向量 $C^{'}$ 。

具体表现在图 3 中，就是向量 $A_{4}^{'}$ 的红色部分变短，而向量 $A_{3}^{'}$ 和向量 $A_{1}^{'}$ 的红色部分变长，但是他们合成的向量 $C^{'}$ 可以保持长度和位置不变。所以，我们可以认为在系数矩阵 $X_{j, t}$ 不满秩的情况下，权重矩阵不唯一。

解决方法如图 4，在采用 “惩罚项” 方法过后，我们将与向量 $B^{'}$ 差异最大 (距离最远) 的向量 $A_{4}^{'}$ 的权重给定为 0 (在图中表现为 $A_{4}$ 的红色部分长度为 0)，那么实际上权重向量 $\hat{W}$ 的未知元素只有 3 个。这样一来，控制组矩阵 $X_{j, t}$ 的秩不会超过 3，能够唯一确定一个权重矩阵 $\hat{W}$ 。

在加州控烟案例中，控烟法案出台前有 19 年的观测区间，但是却有 39 个州作为个体。按照纠偏合成控制法的思想，我们从 38 个控制组个体中寻找与加州各方面差异最大的州，对应的权重赋为 0，再在剩下 37 个控制组个体中寻找差异最大的州，赋权重为 0。这样一直进行下去，直到只有 19 个控制组个体没有被我们主动赋权重为 0 为止，再用传统的合成控制法估计权重矩阵。

需要读者们注意的是，虽然事实上大部分控制组个体的权重都是 0，但是差异较大的控制组个体是我们在估计之前主动赋值的，而其他为 0 的个体是通过最小均方误差估计出来的。

除了图示法，我们还可以从 “惩罚项” 的定义出发去理解权重矩阵被唯一确定的过程。

惩罚实际上是一种 “权衡” (trade off)。在合成控制法中，如果某个控制组个体的的信息与处理组个体的信息差异过大，那么为了 “惩罚” 这一差距，我们需要赋予这些差异过大的个体一个较低的权重。而根据定义，权重矩阵的范围是 $[0, 1)$ ，那么差异最大的个体应该赋予最低的权重，也就是 0，这也就照应了前文的图解法。

3.2 纠偏合成控制法的惩罚项估计

在了解了 “惩罚项” 的原理之后，细心的读者会发现，还有一个至关重要的问题没有解决：那就是最优化条件中，我们给惩罚项添加了一个正的系数 $λ$ ，所以此时求出来的权重矩阵 $\hat{W}$ 还是一个 $λ$ 相关的矩阵，即：

\hat{W} = \hat{W} (λ)

那么问题就是：如何通过实际数据判定惩罚项系数 $λ$ 的大小呢？

联系 Bootstrap 和 Jackknife 的思想，Abadie (2021) 使用了 "Pre-Intervention Holdout Validation" (政策前时间点有放回验证) 和 "Leave-One-Out Cross-Validation" (去一法交叉验证)。具体方法如下：

3.2.1 政策前时间点有放回验证

这一思想和 Bootstrap 有异曲同工之妙。Bootstrap 的基本思想是 “有放回地抽样”，最后根据抽样得到的数据估计总体的特征。"Pre-Intervention Holdout Validation" 则是针对政策实施前的一段时间 (这段时间开始的时间点记为 $k$ )，有放回地抽取不同时间点上处理组和控制组 (比如，在 1980-1988 这 9 年间有放回地随机取 9 个年份作为一组样本，重复多次)，计算出每次抽取得到的 “合成加州” 和 “真实加州” 的差异。

具体的数学表达式如下，首先需要计算抽出的样本中真实加州和合成加州的差异 gap：

\hat{τ} (λ) = Y_{i, t} - \sum_{j = n 1 + 1}^{n} {\hat{W}}_{i, j} (λ) Y_{j, t} (k < t < T_{0})

其中， $k$ 是政策出台时间点 $T_{0}$ 之前的任意一个时间点。注意，时间段 $(T_{0} - k)$ 必须是连续的，时间段 $(T_{0} - k)$ 必须紧跟着政策处理时间点 $T_{0}$ 。而为了让政策实施前加州和合成加州的差异尽可能小， $\hat{τ}$ 则需要满足以下条件：

m i n \sum_{i = 1}^{n_{1}} \sum_{t = T_{0} - k + 1}^{T_{0}} (\hat{τ} (λ))^{2}

我们随机抽得几个政策发生前时间点上的数据，最小化抽取样本中处理组和控制组的差异，既在最大程度上保证了 Bootstrap 抽样数据的完整性，又能够得到更广泛的样本数据。经过上述优化问题我们可以求出 $λ$ 的值，代入权重矩阵 $\hat{W}$ 的计算过程中，我们可以最终计算出 $\hat{W}$ ，完成合成控制法。

3.2.2 去一法交叉验证

这一思想和 Jackknife 方法很相似。Jackknife 的基本思想是 “轮流去除一个地抽样”，最后根据抽样得到的数据估计总体的特征。"Leave-One-Out Cross Validatio" 则是针对政策实施后一直到研究结束的时间段内，轮流挑选一个控制组个体，记该个体为地区 $m$ ，该地区与处理组的一个个体 $i$ 互换，于是地区 $i$ 将被视为控制组个体，而地区 $m$ 将被视为处理组个体，这样维持处理组个体和控制组个体的数量不变。

之后计算出每次抽取得到的合成地区 $m$ 的各项数据，再计算真实地区 $m$ 和合成地区 $m$ 在政策实施之后的差异，因为地区 $m$ 本身没有受到政策影响，所以我们期望它们没有差异。具体的数学表达式如下，首先计算的是政策发生后真实地区 $m$ 和合成地区 $m$ 的差异：

\hat{τ} (λ) = Y_{m, t} - \sum_{j = n_{1} + 1, j \neq m}^{n} {\hat{W}}_{i, j} (λ) Y_{j, t}, （ T_{0} < t < T ）

因为控制组地区实际上没有受到政策冲击，因此，我们需要让政策实施后被挑选出的地区 $m$ 和合成地区 $m$ 的差异尽可能小，基于此 $\hat{τ}$ 则需要满足以下条件：

m i n \sum_{m = n_{1} + 1}^{n} \sum_{t = T_{0} + 1}^{T} (\hat{τ} (λ))^{2}

其中 $T_{0}$ 指的是政策发生的当期。我们通过随机抽得政策发生后的一个控制组地区，并使得控制组地区 $m$ 和合成地区 $m$ 的差异最小化。最终通过这种方法能够计算出 $λ$ ，代入权重矩阵 $\hat{W}$ 中，我们可以计算出 $\hat{W}$ ，完成对权重矩阵的估计。

3.3 纠偏合成控制法下的 gap 估计

值得一提的是，在 allsynth 命令下，全时间段的 gap 估计采用的是 Abadie 和 Imbens (2011) 的估计方法。估计方法如下：

{\hat{γ}}_{b c, t} = \frac{1}{n_{1}} \sum_{i = 1}^{n_{1}} [(Y_{i, t} - {\hat{μ}}_{i, t}) - \sum_{j = n_{1} + 1}^{n} {\hat{W}}_{i, j} (Y_{j, t} - {\hat{μ}}_{j, t})]

这种估计方法被称为 bias-correct (简写为 bc) 估计，这也正是 allsynth 命令被称之为 “纠偏合成控制法” 的原因。其中,

{\hat{μ}}_{i, t} = E [Y | X = x_{i, t}]

{\hat{μ}}_{j, t} = E [Y | X = x_{j, t}]

其中， $i$ 代表处理组个体， $j$ 代表控制组个体。因此， ${\hat{μ}}_{i, t}$ 表示处理组个体的被解释变量的估计值， ${\hat{μ}}_{j, t}$ 表示控制组个体被解释变量的估计值，估计方法一般采用 OLS 回归或者岭回归。又由于权重矩阵是同时用 $Y$ 和 $X$ 的信息估计出来的，因此在估计 gap 时既需要考虑直接用控制组个体的被解释变量 $Y$ 加权，又需要考虑用解释变量 $X$ 估计出的被解释变量加权，二者之间还存在一个 bias，在估计最终的 gap 时需要考虑这一点。这便是用 bias-correct 方法计算 gap 的基本原理。

4. 纠偏合成控制法的 Stata 实操

4.1 代码解析

纠偏合成控制法会用到外部命令 allsynth，安装方式如下。建议读者安装时将下述两条命令写到 do 文档同时执行。

net from https://justinwiltshire.com/s
net install allsynth, replace

另外，在使用 allsynth 命令时，还需要安装额外的外部包。

ssc install elasticregress, replace

allsynth 的命令语法如下：

allsynth depvar,trunit(#) trperiod(#)
      [counit(numlist) xperiod(numlist) mspeperiod()
      resultsperiod() nested allopt
      unitnames(varname)
      figure keep(file) bcorrect(string)
      gapfigure(string)
      pvalues placeboskeep customV(numlist) optsettings]

allsynth 命令大部分语法与 synth 是相同的，如果我们仿造前文提到的 synth 命令，我们仅仅把 synth 命令修改为 allsynth 命令，看看会发生什么结果？

lxhuse smoking.dta, clear
xtset state year
allsynth cigsale beer(1984(1)1988) lnincome retprice      ///
    age15to24 cigsale(1988) cigsale(1980) cigsale(1975),  ///
    trunit(3) trperiod(1989)

这时，输出窗口会和 synth 命令出一致，但是在采用纠偏合成控制法之后，会有更多的州获得一定的权重，并且此时结果窗口会有一行提示：

If bcorrect() was specified but no option was input inside
   the parentheses, bias-correction was not applied.
   bcorrect() must contain one of -nosave-, -merge-,
   or -replace- for bias-correction to be applied

这一提示告诉我们，在使用 allsynth 的时候需要注意使用 borrect()，其参数为 nosave、merge 和 replace 三种选项中的其中一个。具体说明如下：

如果填入了 merge 和 replace，则 Stata 会保存纠偏合成控制法计算出的 gap。因此，需要在随后的代码中添加选项 keep()，并且指定保存的文件名称填入其中；
如果 borrect() 中指定的选项是 merge，那么 Stata 会同时将传统的 gap 和纠偏合成控制法计算出的 gap 保存在文件中，并且存为变量 gap；如果指定的选项是 replace，那么 Stata 仅仅将纠偏合成控制法控制法计算出的 gap 保存在指定文件中，存储为变量 gap；如果指定的是 nosave，那么 Stata 将不会保存相关文件。

相对于 synth，allsynth 命令其他的特有选项含义如下：

gapfigure()：绘图命令，可以绘制 classic (传统合成控制法) 计算出的 gap、bcorrect (纠偏合成控制法) 计算出的 gap 以及 placebos (安慰剂检验) 计算出的 gap，不过最多绘制两个方法计算出的 gap。另外，需要指定绘制线条的形状，常用为 lineback (实心曲线)。
pvalues()：计算每次安慰剂检验时，政策发生后的预测均方误差对应的 P 值，这一值越小，表明政策效应能够较好地通过安慰剂检验；
placeboskeep：保存安慰剂检验的结果，与 keep() 同时使用；
keep()：将 bcorrect() 和placeboskeep 需要保存的结果存入当前路径下的一个 dta 数据集中，注意数据集名称不能合同一路径下的其他文件重名；
关于是否有唯一解的问题：值得一提的是，使用新开发的 allsynth 命令时，如果权重矩阵的取值不唯一，那么 Stata 会出现 warning 提示。

4.2 效果展示

首先，比较重要的结果是每个控制组个体分别赋予了多少权重，在实际的论文写作中可以以此来验证现实直觉。其次，我们需要观察的数据是，在观察政策处理前，合成加州和真实加州的数据差异是否明显，我们期望是不明显。由于这里 allsynth 命令结果与 synth 类似，这里不再展示。最后，是研究者们关心的 gap 时间变化趋势图。具体代码如下：

lxhuse smoking.dta, clear
xtset state year

// 遍历每个州进行纠偏合成控制, 每个州的数据都会存为一个 dta 数据集
forval i=1/39{
    qui allsynth cigsale beer(1984(1)1988) lnincome retprice age15to24 cigsale(1975) ///
    cigsale(1980) cigsale(1988), trunit(`i') trperiod(1989) bcorrect(replace) keep(allsynth_`i') rep
}
// 生成 gap
forval i=1/39{
   use allsynth_`i', clear
   rename _time years
   gen gap_`i' = _Y_treated - _Y_synthetic
   label variable gap_`i' gap
   keep years gap_`i'
   drop if missing(years)
   save allsynth_`i', replace
}
// 合并所有数据集
use allsynth_1, clear
   forval i=2/39{
   qui merge 1:1 years using allsynth_`i', nogenerate
}
save allsynth.dta,replace

use allsynth.dta, replace  
// 因为发现第 22 个州和第 29 个州的数据有异常, 因此我们选择剔除
local lp1
forval i=1/21 {
   local lp1 `lp1' line gap_`i' years, lpattern(dash) lcolor(gs8) ||
}
local lp2
forval i=23/28 {
   local lp2 `lp2' line gap_`i' years, lpattern(dash) lcolor(gs8) ||
}
local lp3
forval i=30/39 {
   local lp3 `lp3' line gap_`i' years, lpattern(dash) lcolor(gs8) ||
}

// 绘图, 设定颜色、图例位置、显示的图例
twoway `lp1' `lp2' `lp3'|| line gap_3 years, scale(1.3) lcolor(black) ///
    legend(pos(6) order(1 "其他控制组" 38 "加州")) xline(1989, lpattern(dash))
// 删去了两个州, 因此黑色线的加州是第 38 条

从图 5 中我们观察到，位于黑色实线加州以下的虚线在政策发生后只有 2-3 条。我们可以计算出政策效应是由于偶然原因的犯错概率约为 $\frac{2}{39} \approx 0.05$ 以及 $\frac{3}{39} \approx 0.08$ 。在 10% 的显著水平下，我们可以认为加州控烟法案的政策效应通过了安慰剂检验。

最后，Stata 会自动生成一个数据集，数据集的名称由 keep() 中填入的名称决定。数据集中包含的重要的变量有：gap、gap_bc、rmpse (预测均方根误差，政策开始后的年份才有计算结果)、P 和 P_bc (预测均方根误差 RMPSE 的 P 值)，读者在论文写作时可以根据自身需求展示结果和数据。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 合成控制 bootstrap jackknife, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]