Stata：基于Lasso的合成控制法-scul

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

作者：郭盼亭 (厦门大学)
邮箱：gpting2020@163.com

目录

1. 什么是合成控制法？

合成控制法 (Synthetic Control Method) 是一种用来估计因果处理效应的方法。在 Abadie、Diamond 和 Hainmueller (2010) (ADH，2010) 开发的典型应用中，研究人员观察了一个接受处理单元和许多未接受处理单元的时间序列结果。具体来看，

• 首先，通过尽量减少合成对照组和处理组在未接受处理时的差距，来选择合适的权重。未受处理序列的加权平均值被用作受处理序列的反事实估计，被称为合成对照组。
• 其次，通过计算处理发生之后，处理组和合成对照组之间的差异得到处理效应的估计值。
• 最后，通过安慰剂分析来进行统计推断。

1.1 识别假设

合成对照组可以进行外推的假设前提是：如果不接受处理，处理组和对照组之间的横截面相关性在处理发生后将保持相同，即合成控制法假设。具体分析如下图所示：

合成控制方法识别因果效应所需的两个关键假设：

• 条件独立性假设；
• 没有异质的时间趋势。

1.2 构建合成对照组

假设 $x_t=(y_{1t},\dots ,y_{St})$ 是 $t$ 时期的对照组的结果向量，$\omega ={({\pi }_1,\dots ,{\pi }_S)}^T$ 是一个 $S\times 1$ 的权重向量，有：

1.3 求解权重向量

其中，$V$ 是 $M$ X $M$ 的变量重要性权重矩阵。计算方法如下：

• 令 $X_1=(Z_1$, ${\overline{Y}}_i^{K_1}$,...,, ..., ${\overline{Y}}_i^{K_M}$ )' 为处理组干预前特征的一个 ($k\times 1$) 向量，$X_0$ 是控制组的 ($k\times J$) 阶矩阵；
• 向量 $W^{\ast }$ 的选择目标就是在权重约束下最小化 $||X_1—X_{0}W||$。

1.4 理论推导

或者

假设存在 R 个随时间变化的因子 $f_t$；每个单位 (如国家、参与者) 依据因子载荷 ${\lambda }_i$ 选取这些因子的某个固定线性组合；对于处理组和控制组，这些 “混杂变量” 在干预前的特征都得到呈现，因而我们可以使用这种信息来 “平衡” 混淆因素；令 $W=(w_2,...,w_{J+1}{)}^{\prime }$ 对于 $w_j>=0$，并且 $w_2+...+w_{J+1}=1$；令 ${\overline{Y}}_i^{K_1}$,...,, ..., ${\overline{Y}}_i^{K_M}$ 代表干预前结果变量的 M 个线性组合 $(M>R)$。

假设我们可以选择 $W^{\ast }$ 满足：

当 $T_0$ 足够大时，${\tau }_{1t}$ 的一个近似无偏估计量可以表示为：

其中，$Y_{1t}$ 代表控制组，$\sum _{j=2}^{J+1}{w}_j^{\ast }{Y}_{jt}$ 代表加权之后的控制组。

2. 什么是 SCUL？

SCUL (Synthetic Control Using Lasso) 顾名思义就是应用 Lasso 的合成控制法，是在 ADH (2010) 提出的经典合成控制法的基础上，创新性地使用 Lasso 回归，来求解构建合成对照组的权重的一种特殊的合成控制法。

可以看出，${\hat{\omega }}_{\text{lasso }}$ 和 ${\hat{\omega }}_{\text{Synth }}$ 的区别就在于加入了一个 $\lambda |\omega {|}_1$。其中，Lasso 回归是一种在机器学习中广泛使用的方法，应用了 Lasso 回归的 SCUL 可以更好地进行样本外预测。除此之外，SCUL 还在以下几个方面进行了创新性贡献：

• 阐明合成控制方法的识别假设；
• 为通常的临时选择提供建议；
• 使用机器学习方法，允许高维的对照组备选池和自动选择模型；
• 允许各种各样的对照组变量类型；
• 允许外推和负权重。

由上所述可知，SCUL 的估计流程和经典合成控制法的区别在于权重 $w$ 的估计，其余的步骤大致一样。

利用机器学习相关文献中经常提到的滚动原点交叉验证 (rolling-origin cross-validation) 方法，求出最优的 $\lambda$。这一方法可以更好地进行样本外预测，以及最小化自相关带来的麻烦。

3. 命令介绍

命令安装：

* 安装命令到 stata
net describe http://repec.org/bocode/s/scul
net install scul, replace
net get scul, replace

* 使用 scul 命令之前，还需要安装以下几个命令：
ssc inst tabstatmat, replace
ssc install distinct, replace
lxhget lassopack.pkg, install replace
cnssc install gtools, replace
ssc install coefplot, replace

*基本语法介绍
scul depvar, treated(variable)
   [ahead(number)] [placebos] [plat]
   [times(numlist)] [lambda(string)]
   [covs(varlist)] [cv(string)]
   [scheme(string)]  [rellab(numlist)]
   [obscol(string)] [cfcol(string)]
   [conf(string)] [legpos(integer)]
   [transform(string)] [q(real)]
   [donoradj(string)]
   [before(integer)] [after(integer)]

其中，depvar 是被解释变量，指定为结果变量，该变量不允许存在缺失值。treated 指定处理变量，该变量不允许有缺失值，必须为 0 或 1，并且在处理开始之后必须始终为 1。如果多于一个单元被处理，scul 将采用交错方法分析 (staggered adoption)。options 主要包含以下内容：

• ahead：指定用户希望在训练数据集中预测的提前周期数。有关此过程的详细信息。
• lambda：可选项。根据 Hastie 等 (2019) 第二章的讨论，用户可以指定 lopt(optimal lambda) 或 lse(standard error rule)。
• placebos：可选项。安慰剂检验。
• plat：可选项。需要设定时间，时间代表着用户希望用来做安慰剂的相对发生时间区间。例如，如果干预发生在 2000 年，并且用户指定 5 和 10，则干预将在 1995 年和 1990 年进行评估。
• covs：可选项。指定用来帮助估计反事实的协变量。
• scheme：可选项。指定想要的图形方案。
• legpos：可选项。指定在 Stata 的点钟表示法中，图例出现的位置。
• squerr：如果指定了安慰剂，则需要 Squerr。这一选项删除干预前 RMSPE 比同期处理组的 RMSPE 大一个数量级的安慰剂组的数据。
• obscol 和 cfcol：可选项。指定接受处理的和未接受处理单位线的颜色。
• q：可选项。当 q=0 时，采用 Ridge penalty；当 q=1 时，采用 LASSO。默认情况下，指定为 LASSO。
• conf：可选项。指定从 t 检验绘制置信区间图。
• transform：可选项。将结果变量归一化在干预发生前的时间。

staggered options：主要包含以下内容：

• donoradj：必选项。通过输入 et 或 nt 作为选项来指定如何调整对照组备选池。指定 et 意味着用户希望使用曾经接受过处理的单元，而 nt 意味着只使用从未接受过处理的单元。
• before and after：必选项。用户指定处理效应的平均相对事件时间周期。注意，这只保留了在周期前后有这么多周期的处理单元。假设一个用户有三个处理组的一年月度数据，一个在 1 月处理 (阶段 1)，一个在 3 月处理 (阶段 3)，另一个在 8 月处理 (阶段 8)。如果用户指定 before after(2)，那么第一个处理单元将被删除，因为它没有 2 个阶段的干预前数据。
• rellab：可选项。指定事件时间轴标记的周期。例如，如果指定 (-5(5)5)，则该图将显示干预前 5 个周期到干预后 5 个周期的平均处理效应。注意，该选项不会改变计算处理效应的周期，只会改变处理效应图的显示方式。

4. 具体示例

4.1 示例 1

. use "http://fmwww.bc.edu/repec/bocode/s/scul_basque.dta", clear
. qui xtset
. local lbl: value label `r(panelvar)'
. loc unit ="Basque Country (Pais Vasco)":`lbl'
. loc int_time = 1975
. qui xtset
. cls
. g treat = cond(`r(panelvar)'==`unit' & `r(timevar)' >= `int_time',1,0)
. scul gdpcap, ahead(3) treat(treat) obscol(black) cfcol("170 19 15") legpos(11)

-----------------------------------------------------------------------------------------
Algorithm: Synthetic LASSO, Single Unit Treated
-----------------------------------------------------------------------------------------
First Step: Data Setup
-----------------------------------------------------------------------------------------
Checking that setup variables make sense.
Setup successful!! All variables id (ID), year (Time) and gdpcap (Outcome) pass.
All are numeric, not missing and non-constant.
-----------------------------------------------------------------------------------------
Inspecting our treatment variable...
-----------------------------------------------------------------------------------------
Treatment is measured from 1975 to 1997 (20 pre-periods)
-----------------------------------------------------------------------------------------
               Treated Unit: Basque Country (Pais Vasco)
-----------------------------------------------------------------------------------------
               Control Units: 16 total donor pool units
-----------------------------------------------------------------------------------------
               Specifically: Andalucia, Aragon, Asturias, Baleares (Islas), Canarias, 
               Cantabria, Castilla Y Leon, Castilla-La Mancha, Cataluna, Comunidad Valenciana,
               Extremadura, Galicia, Madrid (Comunidad De), Murcia (Region de), Navarra
               (Comunidad Foral De), Rioja (La),
-----------------------------------------------------------------------------------------
Second Step: Data Reorganizing
-----------------------------------------------------------------------------------------
Reshaping...
Done!
-----------------------------------------------------------------------------------------
Third Step: Estimation
-----------------------------------------------------------------------------------------
Optimizing (optimal lambda)... This could take quite a while...
Warning: lopt is at the limit of the lambda range.
Optimization took 2.98 minutes
------------------------------------------------------
Relative  |
Time to   |
Event     |              gdpcap5        Counterfactual
----------+-------------------------------------------
      -20 |             3.853185             3.8409667
      -19 |             3.945658              3.943313
      -18 |             4.033562             4.0399026
      省略...
       20 |             9.440874             10.183964
       21 |             9.686518             10.439515
       22 |             10.17067             10.802309
------------------------------------------------------

Collection: default
      Rows: rowname
   Columns: colname
   Table 1: 9 x 1
--------------------------------------
                             | Weights
-----------------------------+--------
Asturias                     |   0.279
Canarias                     |  -0.188
Castilla Y Leon              |   1.127
Cataluna                     |   0.129
Comunidad Valenciana         |  -0.858
Extremadura                  |  -0.774
Madrid (Comunidad De)        |   0.302
Navarra (Comunidad Foral De) |   0.489
Rioja(La)                    |   0.272
--------------------------------------

4.2 示例 2

. loc int_time = 1989
. use "http://fmwww.bc.edu/repec/bocode/s/scul_p99_region", clear
. qui xtset
. local lbl: value label `r(panelvar)'
. loc unit ="California":`lbl'
. qui xtset
. g treat = cond(`r(panelvar)'==`unit' & `r(timevar)' >= `int_time',1,0)
. cls
. scul cigsale, ahead(1) treated(treat) obscol(black) cfcol(blue) legpos(7) cv(adaptive)

-----------------------------------------------------------------
Algorithm: Synthetic LASSO, Single Unit Treated
-----------------------------------------------------------------
First Step: Data Setup
-----------------------------------------------------------------
Checking that setup variables make sense.
Setup successful!! All variables region (ID), year (Time) and cigsale (Outcome) pass.
All are numeric, not missing and non-constant.
-----------------------------------------------------------------
Inspecting our treatment variable...
-----------------------------------------------------------------
Treatment is measured from 1989 to 2000 (19 pre-periods)
-----------------------------------------------------------------
               Treated Unit: California
-----------------------------------------------------------------
               Control Units: 8 total donor pool units
-----------------------------------------------------------------
               Specifically: East South Central, West South Central, 
               Mountain, New England, South Atlantic, East North
               Central, West North Central, Middle Atlantic,
-----------------------------------------------------------------
Second Step: Data Reorganizing
-----------------------------------------------------------------
Reshaping...
Done!
-----------------------------------------------------------------
Third Step: Estimation
-----------------------------------------------------------------
Optimizing (optimal lambda)... This could take quite a while...

Optimization took 0.20 minutes
------------------------------------------
Relative  |
Time to   |
Event     |     California  Counterfactual
----------+-------------------------------
      -19 |            123       119.44122
      -18 |            121       119.94731
      -17 |          123.5       126.71828
       省略...
        9 |           52.3       77.715147
       10 |           47.2       77.391106
       11 |           41.6       70.796067
------------------------------------------

Collection: default
      Rows: rowname
   Columns: colname
   Table 1: 2 x 1
----------------------------
                   | Weights
-------------------+--------
East South Central |  -0.242
Mountain           |   1.032
----------------------------

5. 相关推文

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

• 专题：Stata命令

  • Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？
• 专题：回归分析

  • Stata：拉索回归和岭回归-(Ridge,-Lasso)-简介
  • Stata Blogs - An introduction to the lasso in Stata (拉索回归简介)
• 专题：IV-GMM

  • Stata：IV-Lasso经典论文介绍
  • Lasso一下：再多的控制变量和工具变量我也不怕-T217
• 专题：合成控制法

  • Stata：合成控制法介绍-synth2
  • Stata：合成控制法的预测区间-scpi
  • 合成控制法简介
  • Stata：纠偏合成控制法介绍-allsynth
  • Synth_Runner命令：合成控制法高效实现
  • Stata：合成控制法程序分享
  • Stata：合成控制法-synth-命令无法加载-plugin-的解决办法
  • 合成控制法 (Synthetic Control Method) 及 Stata实现
• 专题：机器学习

  • 图解Lasso系列A：Lasso的变量筛选能力
• 专题：答疑-板书

  • FAQs答疑-2021寒假-Stata高级班-Day3-连玉君-RDD-合成控制法
• 专题：机器学习

  • Lasso：拉索中如何做统计推断
  • Stata：拉索开心读懂-Lasso入门

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版
• Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 连享会-主页 和 知乎专栏，700+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis

New！ lianxh 和 songbl 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh