Stata：合成控制法的预测区间-scpi

发布时间：2022-11-17 阅读 1156

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

作者：高净鹤 (东北财经大学)
邮箱：jh.gaook@gmail.com

编者按：本文主要摘译自下文，特此致谢！
Source：Cattaneo M D, Feng Y, Titiunik R. Prediction intervals for synthetic control methods[J]. Journal of the American Statistical Association, 2021, 116(536): 1865-1880. -PDF- -Supplemental Material- -R-

1. 简介

合成控制法 (Synthetic Control Methods，SCM) 是一种使用未受干预个体 (控制组) 的加权平均值，来近似在没有干预的情况下受干预个体 (处理组) 的反事实替身，以此来估计处理效应。合成控制法的实现通常涉及两个主要的估计步骤：

首先，通过回归方法 (通常是有约束的)，仅使用干预前数据将 treated unit(s) 与 control units 匹配；
其次，通过将干预前的匹配权重与干预后的 control units 相结合，获得 treated unit(s) 的反事实结果的预测。

因此，合成控制法提供了干预政策实施后，treated unit(s) (因果) 处理效应的预测或点估计。本文拓展了合成控制法的不确定性量化，提出了 SC 框架下的条件预测区间：将潜在的结果变量作为随机变量，并采用有限样本概率集中的方法，为 treated unit(s) 在干预后的反事实结果制定预测区间，从而提供了一种评估统计不确定性的 (有条件的) 替代推断方法。

2. 理论背景

本文提出的方法是通过 (有条件的) 预测区间来量化不确定性，因为在 SC 框架中，处理效应估计量基于干预前数据构建的 SC 估计权重，是从样本外预测问题中出现的随机变量。该方法推断出的不是通常意义上的置信区间 (即为感兴趣的非随机参数在参数空间中给出一个区域)，而是一个描述在随机变量的支持下可能观察到的新实现的区域的区间。

SC 预测的统计不确定性由两个不同的随机性来源控制：一个是样本内不确定性——由于干预前期 SC 权重的构造 (可能被错误指定)；另一个是样本外不确定性——由于干预后对处理效应分析时不可观测的随机误差。

因此，本文建议的预测区间是在考虑两个随机性来源的情况下构建的：

对于不确定性的第一个来源，提出了一种基于模拟的方法，该方法通过非渐近概率集中来证明其合理性，从而享有概率保证。该方法考虑了 SC 权重的具体构造。
对于第二个不确定性来源，即由于干预后的不可观测误差而导致的样本外预测，讨论了几种基于非参数概率近似和参数概率近似的方法作为原则灵敏度分析的框架。第二个不确定性源较难以非参数方式处理，因此应谨慎考虑其对总体预测区间的影响，文中的方法是采用不可知性敏感性分析。

2.1 设定

考虑具有单个 treated unit 和多个 control units 的标准合成控制框架，允许平稳和非平稳数据。数据可能只包括感兴趣的结果变量，或者包括感兴趣的结果变量加上其他变量。研究观测的是 $T_{0} + T_{1}$ 时期内 ( $t = 1, 2, \dots, T_{0}, T_{0} + 1, \dots, T_{0} + T_{1}$ ) $N + 1$ 个个体 ( $i = 1, 2, \dots, N, N + 1$ )。

在前 $T_{0}$ 时期内，所有的个体都没有接受干预。从 $T_{0} + 1$ 期开始，个体 1 ( $i = 1$ ) 接受干预而其他个体仍然保持不接受干预，此后个体是否接受干预状态保持不变，直到观测期结束。

每一个个体 $i$ 在时期 $t$ 都有两种潜在的结果：接受干预的结果 $Y_{i t} (1)$ 和未接受干预的结果 $Y_{i t} (0)$ 。即，设定包含两个隐形假设：无溢出 (个体 $i$ 的潜在结果只取决于其 $i$ 是否接受干预) 和无预期 ( $t$ 时期的潜在结果只取决于同时期是否接受干预)。我们关注的处理效应是接受干预的个体的结果与在不接受干预的情况下该个体应该出现的结果之间的差异，即：

τ_{t} = Y_{1 t} (1) - Y_{1 t} (0), t > T_{0} (1)

根据所考虑的框架， $τ_{t}$ 可以是随机或非随机的。在本文中，将 $τ_{t}$ 视为随机变量。对于每个个体，我们只能观测到与实际干预状态相对应的潜在结果，定义 $Y_{i t}$ 为实际观测到的结果：

Y_{i t} = {\begin{cases} Y_{i t} (0) & i f i = 2, \dots, N + 1 \\ Y_{i t} (0) & i f i = 1 a n d t \in {1, 2, \dots, T_{0}} \\ Y_{i t} (1) & i f i = 1 a n d t \in {T_{0} + 1, \dots, T_{0} + T_{1}} \end{cases} (2)

这意味着在 $τ_{t}$ 中，treated unit 的潜在结果 $Y_{1 t} (0)$ 对于所有的 $t > T_{0}$ 是无法观测到的。合成控制法的思想是将控制组干预后观察结果进行适当组合，以接近处理组的反事实干预后结果： $Y_{1 t} (0) f o r t > T_{0}$ 。

在 SC 框架中，选择一组权重 $w = (ω_{2}, ω_{3}, \dots, ω_{N + 1})^{'}$ 使给定的损失函数在约束条件下最小化。给定一组估计权重 $\hat{w}$ ，在 $t > T_{0}$ 时，处理组的反事实预测结果为 ${\hat{Y}}_{1 t} (0) = \sum_{i = 2}^{N + 1} {\hat{w}}_{i} Y_{i t} (0)$ 。加权平均 ${\hat{Y}}_{1 t} (0)$ 为被处理个体的合成控制，因为它代表了如何将未干预的个体组合起来，为干预后的个体提供最佳的反事实。

当数据只包含与感兴趣结果有关的信息时，选择 $w$ 使未干预个体结果的加权平均值与干预前一段时间内 treated unit 的结果轨迹很接近。也就是说，权重 $w$ 的选择满足：

\sum_{i = 2}^{N + 1} w_{i} Y_{i t} (0) \approx Y_{1 t} (0), for t = 1, 2, \dots, T_{0}

其中，符号 $\approx$ 的含义因所考虑的具体框架而异。一个主要例子的约束权重为非负且和为 1，并通过约束最小二乘估计 $w$ ：

{({\hat{w}}^{'}, \hat{r})}^{'} \in \underset{w \in W, r \in R}{\arg min} \sum_{t = 1}^{T_{0}} {(Y_{1 t} - Y_{2 t} w_{2} - \dots - Y_{(N + 1) t} w_{N + 1} - r)}^{2} (3)

其中 $r$ 表示截距， $W$ 和 $R$ 表示相应的约束集 (或可行性集)。

然而，SC 方法的这种只考虑结果的方式，不能保证所得到的合成控制个体将在除 (干预前) 结果之外的任何特征上与 treated unit 相似。在实证中往往还会获得其他特征，因此我们还希望确保合成控制在这些其他特征方面与 treated unit 接近。

SC 框架可以通过为这些额外的特性包含额外的方程来处理这种情况，并将综合损失最小化。在这种情况下，让 $l = 1, 2, \dots, M$ 定义为将被匹配产生权重的变量下标，上述最小化问题可推广为：

{({\hat{w}}^{'}, \hat{r})}^{'} \in \underset{w \in W, r \in R}{\arg min} \sum_{l = 1}^{M} \sum_{t = 1}^{T_{0}} v_{t, l} {(Y_{1 t, l} - Y_{2 t, l} w_{2} - \dots - Y_{(N + 1) t, l} w_{N + 1} - r_{l})}^{2} (4)

其中 $\hat{r} = {({\hat{r}}_{1}, \dots, {\hat{r}}_{M})}^{'}$ 和 ${v_{t, l}}_{1 \leq t \leq T_{0}, 1 \leq l \leq M}$ 是正常数，反映了不同方程和周期的相对重要性。

2.2 一般框架

现在讨论一个包含上述两个特殊例子的一般框架，并允许以统一的方式进行协变量调整和非平稳数据。

考虑为 treated unit 的 M 个特征同时构造的合成控制权重，用 $A_{l} = (a_{1, l}, \dots, a_{T_{0}, l})^{'} \in R^{T_{0}}$ 表示，其中 $l = 1, 2, \dots, M$ 。对于每一个特征 $l$ ，都存在 $J + K$ 个变量，可以用来预测或匹配 $T_{0}$ 维向量 $A_{l}$ 。这些 $J + K$ 个变量被分为两组，分别为 $B_{l} = (B_{1, l}, B_{2, l} \dots, B_{J, l}) \in R^{T_{0} \times J}$ 和 $C_{l} = (C_{1, l}, C_{2, l} \dots, C_{K, l}) \in R^{T_{0} \times K}$ 。

更准确的说，对于每一个 $j$ ， $B_{j, l} = (b_{j 1, l}, \dots, b_{j T_{0}, l})^{'}$ 对应于在干预前 $T_{0}$ 期间内观测到的第 $j$ 个个体的第 $l$ 个特征；对于每一个 $k$ ， $C_{k, l} = (c_{k 1, l}, \dots, c_{k T_{0}, l})^{'}$ 是另一个控制变量向量，也可能用于预测同一干预前时间跨度内的 $A_{l}$ 。为了方便，设 $d = J + K M$ 。

合成控制法的目标是在 M 个特征上搜寻一个公共权重向量 $w \in W \subseteq R^{J}$ 和一个系数向量 $r \in R \subseteq R^{K M}$ ，使得 $B_{l}$ 和 $C_{l}$ 的线性组合在所有 $1 \leq l \leq M$ 的情况下尽可能匹配 $A_{l}$ ，设目标通常通过以下优化问题来实现：

\hat{β} := {({\hat{w}}^{'}, {\hat{r}}^{'})}^{'} \in \underset{w \in W, r \in R}{\arg min} (A - B w - C r)^{'} (A - B w - C r) (5)

其中，

A = [\begin{matrix} A_{1} \\ A_{2} \\ ⋮ \\ A_{M} \end{matrix}], B = [\begin{matrix} B_{1} \\ B_{2} \\ ⋮ \\ B_{M} \end{matrix}], C = [\begin{array}{cccc} C_{1} & 0 & \dots & 0 \\ 0 & C_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & C_{M} \end{array}]

可行性集 $W$ 和 $R$ 捕捉了施加的限制。该框架包含了文献中多种先前的合成控制的形式，它们的不同之处在于是否包含额外的协变量，是否假设数据是平稳的，以及使用的约束集 $W$ 和 $R$ 的特定选择，以及其他可能性。

为了进一步理解本文的推断方法，定义相对于一个 $σ$ 域的伪真值 $w_{0}$ 和 $r_{0}$ ：

β_{0} := {(w_{0}^{'}, r_{0}^{'})}^{'} = \underset{w \in W, r \in R}{\arg min} E [(A - B w - C r)^{'} (A - B w - C r) ∣ H] (6)

因此，

A = B w_{0} + C r_{0} + U, w_{0} \in W, r_{0} \in R (7)

其中， $U = (u_{1, 1}, \dots, u_{T_{0}, 1}, \dots, u_{1, M}, \dots, u_{T_{0}, M})^{'} \in R^{T_{0} M}$ 是相对于 $σ$ 域的伪真值残差。也就是说， $w_{0}$ 和 $r_{0}$ 是与 $H$ 条件下 (可能受约束的) 最佳线性预测系数 $\hat{w}$ 和 $\hat{r}$ 相关的均方误差估计。

给定估计权重 $\hat{w}$ 和系数 $\hat{r}$ ， treated unit 干预后反事实结果为：

{\hat{Y}}_{1 T} (0) = x_{T}^{'} \hat{w} + g_{T}^{'} \hat{r} = p_{T}^{'} \hat{β}, p_{T} := {(x_{T}^{'}, g_{T}^{'})}^{'}, T > T_{0}

其中， $x_{T} \in R^{N}$ 是在时间 $T$ 上观察到的 control units 的预测因子向量， $g_{T} \in R^{K M}$ 是在时间 $T$ 上观测到的另一组用户指定的预测因子。包含在 $x_{T}$ 和 $g_{T}$ 中的变量不需要与 $B$ 和 $C$ 中的变量相同，但必须是 $σ$ 域 $H$ 上的一部分。分解 treated unit 的潜在结果：

Y_{1 T} (0) \equiv x_{T}^{'} w_{0} + g_{T}^{'} r_{0} + e_{T} = p_{T}^{'} β_{0} + e_{T}, T > T_{0} (8)

其中 $e_{T}$ 由结构定义。本文假设 $w_{0}$ 和 $r_{0}$ 可能是 $\hat{w}$ 和 $\hat{r}$ 分别集中在概率周围的随机元素，所以称它们为伪真值。处理组的估计处理效应与目标群体的距离为：

{\hat{τ}}_{T} - τ_{T} = (Y_{1 T} (1) - {\hat{Y}}_{1 T} (0)) - (Y_{1 T} (1) - Y_{1 T} (0)) = Y_{1 T} (0) - {\hat{Y}}_{1 T} (0) (9)

在合成控制框架内，将感兴趣的变量 $τ_{T}$ 视为一个随机变量，不称其为“参数”，称 ${\hat{τ}}_{T}$ 为 $τ_{T}$ 的预测，而不是其“估计量”，并着重建立预测区间而不是置信区间，以此表征 ${\hat{τ}}_{T}$ 的不确定性。

2.3 预测区间

{\hat{τ}}_{T} - τ_{T} = Y_{1 T} (0) - {\hat{Y}}_{1 T} (0) = e_{T} - p_{T}^{'} (\hat{β} - β_{0})

其中， $e_{T}$ 是由于错误设定以及干预后阶段 ( $T > T_{0}$ ) 中出现的任何额外噪声造成的样本外误差， $p_{T}^{'} (\hat{β} - β_{0})$ 是来自合成控制权重估计的样本内误差。本文的目标是分别找到这两项的概率边界，从而给出不确定性量化。

本文使用以下引理来构造有效的条件预测区间。假设存在 $M_{1, L}, M_{1, U}, M_{2, L}$ 和 $M_{2, U}$ ，可能取决于 $α_{1}, α_{2}, π_{1}, π_{2} \in (0, 1)$ 和条件 $σ$ 域 $H$ ，使得满足以下条件：

\begin{aligned} P {P [M_{1, L} \leq p_{T}^{'} (β_{0} - \hat{β}) \leq M_{1, U} ∣ H] \geq 1 - α_{1}} \geq 1 - π_{1}, and \\ P {P [M_{2, L} \leq e_{T} \leq M_{2, U} ∣ H] \geq 1 - α_{2}} \geq 1 - π_{2} \end{aligned}

那么，

P {P [{\hat{τ}}_{T} - M_{1, U} - M_{2, U} \leq τ_{T} \leq {\hat{τ}}_{T} - M_{1, L} - M_{2, L} ∣ H] \geq 1 - α_{1} - α_{2}} \geq 1 - π_{1} - π_{2}

该引理提供了一种简单的方法来构造具有 $(α, π)$ 有效的 $H$ 条件预测区间， $α = α_{1} + α_{2}$ 且 $π = π_{1} + π_{2}$ ：

I = [{\hat{τ}}_{T} - M_{1, U} - M_{2, U}, {\hat{τ}}_{T} - M_{1, L} - M_{2, L}]

一旦为 $τ_{T}$ 构建了 $H$ 条件下 $(α ， π)$ 有效的预测区间，treated unit 在干预后时期 $T$ 的反事实结果的类似预测区间 $Y_{1 T} (0)$ 也很容易得到，根据式 (1)，可见：

P {P [Y_{1 T} (1) - Y_{1 T} (0) \in I ∣ H] \geq 1 - α} \geq 1 - π

因此， $[M_{1, L} + M_{2, L} + {\hat{Y}}_{1 T} (0), M_{1, U} + M_{2, U} + {\hat{Y}}_{1 T} (0)]$ 是 $Y_{1 T} (0)$ 的条件有效预测区间。

3. Stata实操

3.1 准备工作

在 R、Python 和 Stata 中都有软件包可以实现本文提出的方法。读者可以自行在 https://nppackages.github.io/scpi/ 中获取关于 scpi 的各种资料。本文重点介绍 Stata 实现方法。Stata 的实现依赖于 Python ，因此需要提前安装 Python，并完成 Stata (16.0 以上) 与 Python (3.8 以上) 的连接。

安装 Python 有两种方式：

直接从 https://realpython.com/installing-python/ 下载安装 Python
下载安装 Anaconda Windows，macOS，Linux

Python 中安装 scpi 包：

pip install scpi_pkg

Stata 中安装 scpi 命令：

lxhget scpi.pkg, install replace // https://github.com/nppackages/scpi/tree/main/stata

3.2 命令介绍

3.2.1 数据准备：scdata

scpi 中的 scdata 为点估计和预测准备数据。该函数以 DataFrame 类对象为输入，输出 scpi 类数据对象，包含上述矩阵 $A$ 、 $B$ 、 $C$ 和干预后预测因子矩阵 $P = (p_{T_{0} + 1}, \dots, p_{T 0 + T})^{'}$ 。该命令允许使用者指定结果变量、待匹配 treated unit 的特征变量，以及逐个特征的协变量调整。

scdata features [if] [in] , id(idvar) time(timevar) outcome(outcomevar) 
  treatment(treatmentvar) dfname(string) 
  [anticipation(#) cointegrated constant pypinocheck]

其中，

id(idvar)：个体变量
time(timevar)：时间变量。
outcome(outcomevar)：结果变量。
treatment(treatmentvar)：处理变量，接受干预 (被处理) 取 1，否则取 0。
dfname(string)：指定保存并将传递给 scest 或 scpi 的 Python 对象的名称。
covadj(string)：用于调整每个特征的变量。如果想为所有特征指定相同的协变量集，则应采取的形式为 covadj("cov1,cov2")。反之，如果必须为每一个特征指定一组不同的协变量，应采取的形式为 covadj("cov1,cov2;cov1,cov3")。注意，在后种情况下，用 ; 分隔的子列表的数量必须等于特征的数量。字列表的顺序也很重要，第一个子列表是用于调整第一个特征的协变量集，以此类推。
anticipation(#)：潜在预期影响的时期数。默认值为 anticipation(0)。
cointegrated：如果该项被指定，意味着认为这些特征构成协整系统。
constant：该项被指定意味着包含跨特征的常数项。
pypinocheck：如果该项被指定，那么可以避免检查 Python 中 scpi_pkg的版本是否是 Stata 中 scdata 所需的版本，如果该项未指定，则执行检查并存储调用的宏，以避免多次检查。

3.2.2 点估计：scest

在创建了所有的设计矩阵 $A$ 、 $B$ 、 $C$ 和 $P$ 后，通过 scest 继续对感兴趣的反事实结果进行点估计/预测。scest 使用最小二乘回归、Lasso 回归、岭回归或单纯形约束来实现合成控制方法的估计程序。

scest, dfname(string) [p(#) direc(string) Q(#) lb(#) name(string) 
  V(string) opt(string) pypinocheck]

其中，

dfname(string)：使用 scdata 创建的已处理数据的 Python 对象的名称。
p(#)：设置要约束的范数类型。取 0 时不对权重的范数施加约束，取 1 时对权重的 L1 范数施加约束 (默认值)，取 2 时对权重的 L2 范数施加约束。
direc(string)：指定权重范数上约束的方向。选项包括：<= 权重范数上的约束是不等式约束，== 权重范数上的约束是等式约束 (默认值)。
Q(#)：指定权重范数上约束的大小。
lb(#)：指定权重的下限，默认是 lb(0)。
name(string)：指定使用的约束的名称：
- simplex：单纯形回归，经典合成控制估计，权重约束为非负且 L1 范数必须等于 1；
- lasso：使用 Lasso 型惩罚估计权重；
- ridge：使用 Ridge 型惩罚估计权重；
- L1-L2：使用单纯形约束和岭型惩罚估计权重；
- ols：使用最小二次惩罚无约束地估计权重。
V(string)：指定损失函数中使用的加权矩阵。默认值是单位矩阵 (V("separate"))，即所有观测值都具有相同的权重；另一种可能是 V("separate")，即集合拟合。
opt(string)：点估计的底层优化算法的终止条件，默认为非线性约束梯度优化 (SLSQP) 的序列二次规划 (SQP)，默认值为 opt("'maxeval' = 5000, 'xtol_rel' = 1e-8, 'xtol_abs' = 1e-8, 'ftol_rel' = 1e-12, 'ftol_abs' = 1e-12, 'tol_eq' = 1e-8, 'tol_ineq' = 1e-8"。当使用 Lasso 型约束，则使用不同的优化算法 (cvxpy)，并且不能更改终止条件。

3.2.3 预测区间：scpi

scpi 使用最小二乘、lasso、ridge 或单纯形类型约束来实现合成控制 (SC) 方法的预测区间估计。

scpi, dfname(string) [p(#) direc(string) Q(#) lb(#) V(string) name(string) u_missp) 
  u_sigma(string) u_order(#) u_lags(#) u_alpha(#) sims(#) e_method(string)
  e_order(#) e_lags(#) e_alpha(#) lgapp(string) rho(#) rho_max(#) opt_est(string) 
  opt_inf(string) pypinocheck]

其中，与上文中相同的选项都具有相同的含义，因此不再做过多的说明，下面介绍其他的选项。

样本内不确定性：

u_missp：如果被指定，则表示考虑了模型设置错误。
u_sigma(string)：估计伪残差的条件方差时要使用的方差-协方差估计量的类型。包括 HCO，HC1 (默认)，HC2，HC3。
u_order(#)：用于估计伪残差条件矩的预测中多项式的阶数。默认值为 u_order(1)。如果存在过度拟合的风险，该选项将自动设置为 0。
u_lags(#)：用于估计伪残差条件矩的预测的滞后时间。默认值为 u_lags(0)。如果存在过度拟合的风险，该选项将自动设置为 0。
u_alpha(#)：样本内不确定性的置信水平，即置信水平为 1-u_alpha。默认值为u_alpha(0.05)。
sims(#)：样本内不确定性的模拟次数。默认值为 sims(200)。

样本外不确定性：

e_method(string)：量化样本外不确定性的方式，选项包括 gaussian 条件亚高斯边界，ls location-scale 模型，qreg 分位数回归，all 上面三种都包括 (默认设置)。
e_order(#)：估计样本外误差的条件矩的预测中多项式的阶数。默认值为 e_order(1)。如果存在过度拟合的风险，该选项将自动设置为 0。
e_lags(#)：用于估计样本外误差条件矩的预测的滞后时间。默认值为 e_lags(0)。如果存在过度拟合的风险，该选项将自动设置为 0。
e_alpha(#)：样本外不确定性的置信水平，即置信水平为 1-e_alpha。默认值为 e_alpha(0.05)。

其他选项：

lgapp(string)：选择模拟中近似局部几何关系的方式，选项包括 generalized 和 linear，前一种方式适用于非线性约束，第二种方式适用于线性约束。
rho(#)：在估计的权重向量上施加稀疏性正则化参数。如果未指定，则根据优化不等式计算调整参数。
rho_max(#)：调优参数 rho 可以达到的最大值。
opt_est(string)：点估计的底层优化算法的终止条件，默认为非线性约束梯度优化 (SLSQP) 的序列二次规划 (SQP)，默认值为 opt("'maxeval' = 5000, 'xtol_rel' = 1e-8, 'xtol_abs' = 1e-8, 'ftol_rel' = 1e-12, 'ftol_abs' = 1e-12, 'tol_eq' = 1e-8, 'tol_ineq' = 1e-8"。
opt_inf(string)：点估计的底层优化算法的终止条件，默认为非线性约束梯度优化 (SLSQP) 的序列二次规划 (SQP)，默认值为 opt("'maxeval' = 5000, 'xtol_rel' = 1e-8, 'xtol_abs' = 1e-8, 'ftol_rel' = 1e-4, 'ftol_abs' = 1e-4, 'tol_eq' = 1e-8, 'tol_ineq' = 1e-8"。

3.2.4 绘图：scplot

最后一步，将估计的合成控制可视化，并利用 scplot 将其与观测到的 treated unit 的时间序列进行比较。

scplot, [scest uncertainty(string) joint dots_tr_col(colorstyle) 
  dots_tr_symb(symbolstyle) dots_tr_size(markersizestyle) dots_sc_col(colorstyle) 
  dots_sc_symb(symbolstyle) dots_sc_size(markersizestyle) line_tr_col(colorstyle) 
  line_tr_patt(linepatternstyle) line_tr_width(linewidthstyle) line_sc_col(colorstyle) 
  line_sc_patt(linepatternstyle) line_sc_width(linewidthstyle) spike_sc_col(colorstyle) 
  spike_sc_patt(linepatternstyle)spike_sc_width(linewidthstyle) gphoptions(string) 
  gphsave(string) savedata(dta_name) pypinocheck]

其中，

scest：如果指定该项，则必须在 scest 之后调用 scplot。否则，假定在 scpi 之后调用 scplot。
uncertainty(string)：绘制的预测区间的类型。选项包括：
- insample：预测区间只量化样本内不确定性。
- gaussian：预测区间使用条件亚高斯边界量化样本内和样本外不确定。
- ls：预测区间使用 location-scale 模型量化样本内和样本外不确定性。
- qreg：预测区间使用分位数回归量化样本内和样本外不确定性。
joint：如果该项被指定，则图中包含同步预测区间。
dots 相关选项：绘图时有关标记的选项 (颜色、大小、形式等)。
line 相关选项：绘图时有关线条的选项 (颜色、图案、宽度等)。
spike 相关选项：绘图时有关条形 (尖峰) 的选项。如果指定了 scest，则这些选项无效。
gphoptions(string)：修改绘图的其他选项。
gphsave(string)：指定命令保存的.gph 文件的路径和名称。
savedata(dta_name)：保存一个名为 dta_name.dta 的文件，包含用于生成绘图的已处理数据。

3.3 具体应用

使用德国的数据，关注 1990 年德国统一对西德人均 GDP 的因果影响。将西德统一后的结果与从 1960 年到 1990 年使用 16 个经合组织国家建立的合成控制个体的结果进行了比较。下面演示仅有一个 treated unit 的实操过程。

3.3.1 一个特征变量 (M=1：gdp)

**加载数据
. lxhuse "scpi_germany.dta", clear

**准备数据
. scdata gdp, dfname("python_scdata") id(country) outcome(gdp) time(year) ///
>     treatment(status) cointegrated constant

**点估计 (单纯形 simplex) 
. scest, dfname("python_scdata") name(simplex) 
. scest, dfname("python_scdata") p(1) q(1) direc("==") lb(0)

-------------------------------------------------------
Call: scest
Synthetic Control Estimation - Setup

Constraint Type:                          user provided
Constraint Size (Q):                                  1
Treated Unit:                              West Germany
Size of the donor pool:                              16
Features                                              1
Pre-treatment period                          1960-1990
Pre-treatment periods used in estimation:            31
Covariates used for adjustment:                       1

Synthetic Control Estimation - Results
Active donors: 6
Coefficients:
                          Weights
Treated Unit Donor               
West Germany Australia      0.000
             Austria        0.441
             Belgium        0.000
             Denmark        0.000
             France        -0.000
             Greece         0.000
             Italy          0.177
             Japan          0.014
             Netherlands    0.058
             New Zealand    0.000
             Norway         0.000
             Portugal       0.000
             Spain          0.000
             Switzerland    0.036
             UK             0.000
             USA            0.274

                       Covariates
                                
West Germany constant       0.158

**画图
. scplot, scest gphoptions("ytitle(GDP per capita) xtitle(Year)")

**预测区间 (单纯形 simplex，模型设定错误) 
. set seed 8894
. scpi, dfname("python_scdata") name(simplex) u_missp u_order(1) u_lags(0) ///
>   u_sigma("HC1") e_order(1) e_lags(0) e_method(qreg)

**画图
. scplot, gphoptions("ytitle(GDP per capita) xtitle(Year)")

3.3.2 多个特征变量 (M=2：gdp，trade)

**多个特征
. scdata gdp trade, dfname("python_scdata") id(country) outcome(gdp) ///
>   time(year) treatment(status) cointegrated
. scpi, dfname("python_scdata") name(simplex) u_missp u_order(1)    ///
>   u_lags(0) u_sigma("HC1") e_order(1) e_lags(0) e_method(qreg)
. scplot, gphoptions("ytitle(GDP per capita) xtitle(Year)") joint

* 多个特征和指定特征协变量调整 
scdata gdp infrate, dfname("python_scdata") id(country) outcome(gdp) ///
  time(year) treatment(status) cointegrated covadj("constant, trend; constant")

4. 总结

本文关注标准的合成控制框架下感兴趣的 $τ_{T} = Y_{1 t} (1) - Y_{1 t} (0)$ 的不确定性量化，即干预后时期内实际观测到的处理组的结果与在未接受干预情况下处理组的结果之间的差。本文将 $τ_{T}$ 视为随机变量，并且处理组只包含一个个体，因此提出条件预测区间，为反事实处理结果的实现提供了有限样本概率保证。

该方法以 SC 约束最小二乘优化方法为出发点，利用控制组在 T 处的特征的加权 (利用干预前数据估计的权重) 和误差项为处理组在时期 T 中的反事实结果建模。这种分解方法突出了两个不确定性来源，一个来自干预前对 SC 权重的样本内估计，另一个来自于 SC 方法中不可不免的样本外预测而产生的干预后误差。利用有限样本浓度边界，得到了包含两个不确定性来源的预测区间。

该过程可以通过 R、Python 和 Stata 的 scpi 软件包实现合成控制方法的点估计和预测区间、不确定性量化的过程，

5. 相关推文

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

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