第9章 合成控制法

在上一章中，你了解了面板数据对于因果识别的优势。也就是说，你不仅可以将个体相互之间进行比较，还可以将它们与过去的自身进行比较，这使得你能够在更合理的假设下估计反事实结果 \(Y_0\)。你还了解了双重差分法\(（DID）\)—— 以及它的许多变体 —— 这是众多利用面板数据的因果推断工具之一。通过依赖处理组和对照组之间相似的平行的增长轨迹，即使处理组和对照组之间 \(Y_0\) 的水平存在差异，双重差分法也能够识别出处理效应。在本章中，你将学习面板数据集的另一种常用技术：合成控制法\(（SC）\)。

虽然当你拥有的个体数量 \(N\) 相对于时间段 \(T\) 较多时，双重差分法\(（DID）\)效果很好，但当情况相反时，它就不够用了。相比之下，合成控制法旨在用于处理个体非常少，甚至只有一个处理个体的情况。其背后的思路相当简单：将控制个体组合起来，构建一个合成控制组，这个合成控制组在没有处理的情况下，能够近似处理个体的行为。通过这样做，它避免了做出平行趋势假设，因为精心构建的合成控制组，不会仅仅是平行的，而是会与反事实 \(E[Y_0|D = 1]\) 完全重叠。

在本章末尾，你还将学习如何将 \(DID\) 和 \(SC\) 结合起来。这种组合估计量不仅非常强大，而且最重要的是，它将让你对双重差分法、特别是合成控制法，以及总体上的面板数据方法有一个全新的认识。

在线营销数据集

作为合成控制法的一个用例，你将使用一个在线营销数据集。与线下营销相比，在线营销能实现更好的追踪，但这并不意味着它在因果推断方面就没有挑战。例如，在线营销确实能实现更优的归因：你能知道客户是否通过某个付费营销链接接触到了你的产品。但这并不意味着你知道，如果客户没有看到你的在线广告，他们会有怎样的行为。或许客户只是因为看到了广告才来的，在这种情况下，广告带来了额外的客户。但也有可能不管怎样客户都会来，而他们通过付费链接进入只是因为那个链接在页面顶部。

由于归因和增量性并非同一回事，而且你无法对谁能看到你的广告进行随机化处理，所以像上一章那样，对整个地理区域进行处理并开展某种面板数据分析，对于在线营销来说也是个不错的想法。因此，你这里的数据与上一章看到的数据并没有太大差异。同样，你以城市为单位，以日期为时间维度，有一个 \(treat\) 列，用于标记该城市最终是否被处理，还有一个 \(post-treatment\) 列，用于标记干预后的时期。你也有一些辅助列，比如该城市的人口（记录于 \(2013\) 年，所以是固定不变的）以及州（信息）。

from toolz.curried import *


import pandas as pd
import numpy as np

import statsmodels.formula.api as smf

import seaborn as sns
from matplotlib import pyplot as plt
import matplotlib

from cycler import cycler

color=['0.0', '0.4', '0.8']
default_cycler = (cycler(color=color))
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']

plt.rc('axes', prop_cycle=default_cycler)

import pandas as pd
import numpy as np

df = (pd.read_csv("./data/online_mkt.csv").astype({"date":"datetime64[ns]"}))

df.head()

##    app_download  population       city      state       date  post  treated
## 0        3066.0    12396372  sao_paulo  sao_paulo 2022-03-01     0        1
## 1        2701.0    12396372  sao_paulo  sao_paulo 2022-03-02     0        1
## 2        1927.0    12396372  sao_paulo  sao_paulo 2022-03-03     0        1
## 3        1451.0    12396372  sao_paulo  sao_paulo 2022-03-04     0        1
## 4        1248.0    12396372  sao_paulo  sao_paulo 2022-03-05     0        1

在这里，结果变量是每日应用下载量 \(app\_download\)，处理变量是为该城市开启营销活动。干预在被处理的城市同时实施，这意味着你采用了一种简单的区组设计。这里的问题在于，被处理的单位数量要少得多 —— 只有 \(3\) 个城市：

treated = list(df.query("treated==1")["city"].unique())
treated

## ['sao_paulo', 'porto_alegre', 'joao_pessoa']

如果你留意数据框中的人口列，可能会注意到其中一个城市 —— 圣保罗\(（sao \space paulo）\)，人口多达 \(1200\) 多万。事实上，它是世界上最大的城市之一！这也意味着圣保罗的应用下载量会比其他城市大得多，这带来了一些挑战。很难把其他城市组合起来，构建一个能与圣保罗下载量相匹配的合成控制组。这个问题在这里更为突出，但一般来说，整个市场规模各不相同，这使得在它们之间进行比较很困难。因此，一种常见的方法是按市场规模对结果进行标准化。这意味着将应用下载量除以城市人口，以创建结果的标准化版本。这个新的结果：\(app\_download\_pct\)，表示每日下载量占市场规模的百分比：

df_norm = df.assign(app_download_pct = 100*df["app_download"]/df["population"])

df_norm.head()

##    app_download  population       city  ... post treated  app_download_pct
## 0        3066.0    12396372  sao_paulo  ...    0       1          0.024733
## 1        2701.0    12396372  sao_paulo  ...    0       1          0.021789
## 2        1927.0    12396372  sao_paulo  ...    0       1          0.015545
## 3        1451.0    12396372  sao_paulo  ...    0       1          0.011705
## 4        1248.0    12396372  sao_paulo  ...    0       1          0.010067
## 
## [5 rows x 8 columns]

继续进行探索性分析，你会发现一场在线营销活动于 \(2022年5月1日\) 在那些城市启动。该活动在分析时间窗口的剩余时间内也持续进行着：

tr_period = df_norm.query("post==1")["date"].min()
tr_period

## Timestamp('2022-05-01 00:00:00')

现在是回顾你将要使用的一些面板数据符号的好时机。回想一下，为了避免混淆，我会用 \(D\) 表示处理变量，用 \(t\) 表示时间。\(T\) 表示时期数，其中 \(T_{pre}\) 是干预前的时期数，\(T_{post}\) 是干预后的时期数。因此，当 \(D = 1\) 且 \(t > T_{pre}\) 时，处理发生。为了简化，我有时会用一个 \(post\) 虚拟变量来表示 \(t > T_{pre}\)。处理组和处理后时期的组合将用 \(W_{it} = D_i * Post_t\) 来表示。

为了让你了解这些数据的样子，下面的图表展示了 \(3\) 个处理城市的平均结果的变化情况，背景中还有一组对照城市的样本（以浅灰色显示）。处理后时期的开始处用一条水平虚线标记：

看这个图表，你大致能看到处理组个体在干预后的结果有所上升，但并非 100% 清晰。要更精确的话，你必须估计反事实情况，并将其与观察到的结果进行比较，以得到处理组的平均处理效应（ATT）：

\(ATT = E[Y|D = 1, Post = 1] - E[Y_0|D = 1, Post = 1]\)

这就是合成控制法发挥作用的地方。它是一种极其巧妙的方法，利用（但不依赖于）过去的结果来估计 \(E[Y_0|D = 1, Post = 1]\)。

np.random.seed(123)

df_sc = df_norm.pivot(index="date", columns="city", values="app_download_pct")
ax = df_sc[treated].mean(axis=1).plot(figsize=(15,6), label="$E[Y|D=1]$")
ax.vlines(tr_period, 0, 0.07, ls="dashed", label="$t=T_{pre}+1$")
ax.legend()
# leave the control group
df_sc.drop(columns=treated).sample(frac=0.2, axis=1).plot(color="0.5", alpha=0.2, legend=False, ax=ax)
plt.ylabel("app_download_pct")

看这个图表，你大致能看到处理组个体在干预后的结果有所上升，但并非 \(100\%\) 清晰。要更精确的话，你必须估计反事实情况，并将其与观察到的结果进行比较，以得到处理组的平均处理效应：

\(ATT = E[Y|D = 1, Post = 1] - E[Y_0|D = 1, Post = 1]\)

这就是合成控制法发挥作用的地方。它是一种极其巧妙的方法，利用但不依赖于过去的结果来估计 \(E[Y_0|D = 1, Post = 1]\)。

矩阵表示

在上一章中，我向你展示过一幅将面板数据表示为矩阵的图像，其中一个维度是时间段，另一个维度表示个体单位。合成控制法会明确使用这个矩阵，所以有必要回顾一下。假设矩阵的行是时间段，矩阵的列是城市（个体）。你可以用四个块来表示处理分配：

\(W = \begin{bmatrix} 0_{pre,co} & 0_{pre,tr} \\ 0_{post,co} & 1_{post,tr} \end{bmatrix}\)

矩阵中的第一个块（左上角）对应处理期之前的控制组个体；第二个块（右上角）对应处理期之前的处理组个体；第三个块（左下角）包含处理期之后的控制组个体；第四个块（右下角）是处理期之后的处理组个体。处理指示变量 \(w_{ti}\) 除了处理期之后的处理组个体所在的块（右下角）外，其他地方均为零。

这个分配矩阵会得到以下观测到的潜在结果矩阵：

\(Y = \begin{bmatrix} Y(0)_{pre,co} & Y(0)_{pre,tr} \\ Y(0)_{post,co} & Y(1)_{post,tr} \end{bmatrix}\)

再次注意，处理后时期在下方，处理组个体在右侧。你的目标是估计 \(ATT = Y(1)_{post,tr} - Y(0)_{post,tr}\)。为此，你需要以某种方式估计缺失的潜在结果 \(Y(0)_{post,tr}\)，它是未被观测到的。换句话说，你需要知道，如果处理组个体在处理后时期没有被处理，会发生什么情况。为了实现这一点，理想情况下你会利用可供你使用的其他三个块：\(Y(0)_{pre,co}\)、\(Y(0)_{pre,tr}\) 和 \(Y(0)_{post,co}\)。在我向你展示合成控制法是如何做到这一点之前，让我们先创建一个函数，以这种矩阵格式来表示数据。

以下代码使用 \(.pivot()\) 方法来重塑数据框，以便最终每个时间段日期对应一行，每个城市对应一列，同时结果变量成为矩阵的值。然后，它将矩阵划分为处理组个体和控制组个体。进一步将它们划分为干预前时期和干预后时期：

def reshape_sc_data(df: pd.DataFrame,geo_col: str, time_col: str,y_col: str,tr_geos: str,tr_start: str):
    
    df_pivot = df.pivot(index=time_col, columns=geo_col, values=y_col)
    y_co = df_pivot.drop(columns=tr_geos)
    y_tr = df_pivot[tr_geos]
    
    y_pre_co = y_co[df_pivot.index < tr_start]
    y_pre_tr = y_tr[df_pivot.index < tr_start]
    
    y_post_co = y_co[df_pivot.index >= tr_start]
    y_post_tr = y_tr[df_pivot.index >= tr_start]
    
    return y_pre_co, y_pre_tr, y_post_co, y_post_tr

在本章中，你将一直使用这种四块矩阵表示法。如果你忘记了自己正在处理的内容，只需回到这个函数。为了了解它的工作原理，将 \(df\_norm\) 传递给 \(rshape\_sc\_date\) 会以矩阵格式返回这些 \(Y\) 值。以下是 \(y\_pre\_tr\) 的前五行：

y_pre_co,y_pre_tr,y_post_co,y_post_tr = reshape_sc_data(
  df_norm,
  geo_col="city",
  time_col="date",
  y_col="app_download_pct",
  tr_geos=treated,
  tr_start=str(tr_period)
)
y_pre_tr.head()

## city        sao_paulo  porto_alegre  joao_pessoa
## date                                            
## 2022-03-01   0.024733      0.004288     0.022039
## 2022-03-02   0.021789      0.008107     0.020344
## 2022-03-03   0.015545      0.004891     0.012352
## 2022-03-04   0.011705      0.002948     0.018285
## 2022-03-05   0.010067      0.006767     0.000000

作为水平回归的合成控制

合成控制法背后的核心思想相当简单。利用处理前时期，你要找到一种方法，将控制组个体组合起来，以近似处理组个体的平均结果。从数学角度而言，这可以构建为一个优化问题：你要寻找个体权重 \(\omega_i\)（不要与 \(w_{it} = Post_t * D_i\) 混淆），使得当你将每个权重与其对应个体的结果相乘（即 \(\omega_i y_i\)）时，得到的结果与处理组个体的结果相似：

\(\widehat{\omega}^{sc} = \underset{\omega}{\text{argmin}} \ \Vert \overline{y}_{pre,tr} - Y_{pre,co} \omega_{co} \Vert^2\)

这段内容主要是关于用合成控制法估计反事实结果以得到平均处理效应，并将合成控制与线性回归类比来解释。

首先，为了估计 \(E[Y(0)|D = 1, Post = 1]\) 并得到 \(ATT\) 估计值，你可以使用合成控制 \(Y_{post,co}\omega_{co}\)。

如果这看起来有点晦涩，或许一个好的替代解释是将合成控制与更熟悉的工具 —— 线性回归进行比较。回想一下，回归也可以表示为一个优化问题，其目标是最小化结果与协变量 X 的线性组合之间的（平方）差异：

\(\beta^* = \underset{\beta}{argmin} \ \Vert Y_i - X_i'\beta \Vert^2\)

---

结果建模

在这里，你可以将合成控制法与潜在结果建模进行类比，潜在结果建模是你在第\(5\)章阅读关于双重稳健估计的内容时看到的。在那里，你也必须建立一个在对照组中估计的回归模型。然后，你使用该模型来填补那些被处理个体的缺失潜在结果 \(Y_0\)。这里的思路几乎是一样的。

---

如你所见，两个目标是完全相同的！这意味着合成控制法只不过是一种回归，它使用对照组的结果作为特征，试图预测处理组个体的平均结果。巧妙之处在于，它仅使用干预前时期的数据来进行回归，从而估计 \(E[Y_0|D = 1]\)。

事实上，为了证明我的观点，现在让我们用 \(OLS\) 来构建一个合成控制组。你所要做的就是把 \(y\_pre\_co\) 当作协变量矩阵 \(X\)，把 \(y\_pre\_tr\) 的列均值当作结果 \(y\)。一旦你拟合了这个模型，就可以用 \(.coef\_\) 提取权重：

from sklearn.linear_model import LinearRegression

model = LinearRegression(fit_intercept=False)
model.fit(y_pre_co, y_pre_tr.mean(axis=1))

LinearRegression(fit_intercept=False)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# extract the weights
weights_lr = model.coef_
print(weights_lr.round(3))

## [-0.65  -0.058 -0.239  0.971  0.03  -0.204  0.007  0.095  0.102  0.106
##   0.074  0.079  0.032 -0.5   -0.041 -0.154 -0.014  0.132  0.115  0.094
##   0.151 -0.058 -0.353  0.049 -0.476 -0.11   0.158 -0.002  0.036 -0.129
##  -0.066  0.024 -0.047  0.089 -0.057  0.429  0.23  -0.086  0.098  0.351
##  -0.128  0.128 -0.205  0.088  0.147  0.555  0.229]

---

如你所见，每个对照城市都有一个权重。通常，回归用于当你有大量个体（大 \(N\)）的情况，这使得你可以将个体作为行，协变量作为列。但合成控制法旨在用于个体相对较少，但干预前时间跨度 \(T_{pre}\) 较大的情况。为了实现这一点，合成控制法实际上是把数据倒过来，将个体当作协变量来使用。这就是为什么合成控制法也被称为水平回归（见\(图9-1\)）。

一旦你估计出了回归参数或权重，你就可以用它们来预测 \(E[Y_0|D = 1]\) 的情况，不仅是在干预前时期，而且是在整个时间范围内：

# same as y0_tr_hat = model.predict(y_post_co)
y0_tr_hat = y_post_co.dot(weights_lr)
print(y0_tr_hat)

## date
## 2022-05-01    0.013356
## 2022-05-02    0.018682
## 2022-05-03    0.016760
## 2022-05-04    0.014669
## 2022-05-05    0.012445
##                 ...   
## 2022-06-26    0.019992
## 2022-06-27    0.029849
## 2022-06-28    0.024227
## 2022-06-29    0.028336
## 2022-06-30    0.029129
## Length: 61, dtype: float64

在这里，\(y0\_tr\_hat\) 可以被视为一个合成控制组：它是对照组个体的一种组合，这些个体结合在一起，近似模拟了如果处理组个体未被处理时，其平均表现的情况。

---

合成控制组的平均值

或者，不是寻找一个合成控制来复制处理组单元的平均结果，你也可以为每个处理组单元单独拟合一个合成控制，然后对这些合成控制取平均值：

# model = LinearRegression(fit_intercept=False)
# model.fit(y_pre_co, y_pre_tr)
# y0_tr_hat = model.predict(y_co).mean(axis=1)

---

如果你把这个合成控制组和观测到的结果一起绘制出来，你会得到这样的图：

plt.figure(figsize=(10,4))

y_co = pd.concat([y_pre_co, y_post_co])
y_tr = pd.concat([y_pre_tr, y_post_tr])

plt.plot(y_co.index, model.predict(y_co), label="Synthetic Control", ls="-.")
plt.plot(y_tr.mean(axis=1), label="Treated Avg.")
plt.vlines(pd.to_datetime("2022-05-01"), 0, 0.04, ls="dashed", label="$t=T_{pre}+1$")
plt.xticks(rotation=45)

## (array([19052., 19066., 19083., 19097., 19113., 19127., 19144., 19158.,
##        19174.]), [Text(19052.0, 0, '2022-03-01'), Text(19066.0, 0, '2022-03-15'), Text(19083.0, 0, '2022-04-01'), Text(19097.0, 0, '2022-04-15'), Text(19113.0, 0, '2022-05-01'), Text(19127.0, 0, '2022-05-15'), Text(19144.0, 0, '2022-06-01'), Text(19158.0, 0, '2022-06-15'), Text(19174.0, 0, '2022-07-01')])

plt.ylabel("app_download_pct")
plt.legend(fontsize=14)

注意预测值（合成控制）是如何低于处理组单元的实际结果的。这意味着，倘若没有进行处理（干预），观察到的结果要高于你所估计的结果。这表明在线营销活动产生了积极的营销效果。你可以通过将观察到的结果与合成控制进行对比来计算平均处理效应的估计值：

att = y_post_tr.mean(axis=1) - y0_tr_hat

该图表呈现出几个有趣的方面。首先，它表明这种影响需要一些时间才能达到峰值，之后才会逐渐减弱。在营销领域，这种逐渐上升的情况很常见，因为人们通常在看到广告后需要时间来采取行动。此外，影响逐渐消失往往可以归因于新奇效应，这种效应会随着时间的推移逐渐消退。

第二件有趣的事是干预前时期平均处理效应的大小。在那个时间段内，\(ATT\) 可以简单地解释为 \(OLS\) 模型的残差（样本内误差）。你可能会认为它接近 \(0\) 是件好事；毕竟，你不希望在处理之前就看到效果（预期效应）。但实际上还有更多情况。干预前误差极低这一事实也可能意味着 \(OLS\) 模型很可能存在过拟合。因此，本应估计 \(E[Y_0|D = 1, Post = 1]\) 的样本外预测可能会不准确。

这就是为什么简单回归不常被用作构建合成控制组的方法。由于列数（对照城市的数量）相对较多，它往往会过拟合，无法推广到干预后时期。正因如此，最初的合成控制方法不是简单的回归，而是一种施加了一些合理且直观约束的方法。

标准合成控制法

标准的合成控制法公式对回归模型施加了两个约束：

1. 所有权重均为正。

2. 权重之和为 1。

从数学角度看，优化目标变为：

\(\widehat{\omega}^{sc} = \underset{\omega}{\text{argmin}} \ \Vert \overline{\mathbf{y}}_{pre,tr} - \mathbf{Y}_{pre,co}\omega_{co} \Vert^2\)

约束条件：\(\sum \omega_i = 1\) 且 \(\omega_i > 0 \ \forall \ i\)

这些约束背后的理念是迫使合成控制组成为处理组个体的凸组合，避免外推。这意味着，如果处理组个体的结果高于（或低于）所有控制组个体的结果，这种标准公式将无法构建出一个合成控制组来还原 \(E[Y_0|D = 1]\)。你可以将这视为一种局限性，但实际上它是一种防护措施。这相当于在表明，你试图重建的处理组个体与对照组中的个体差异非常大，因此你甚至不应该尝试用这种方法。

要编写标准版本的合成控制代码，你可以使用凸优化软件，比如 \(cvxpy\) 。\(cvxpy\) 允许你使用 \(cp.Minmize\) 来定义优化目标。对于合成控制，你想要最小化平方误差，所以用 \(cp.Minimize(cp.sum_squares(y\_co\_pre@\omega - y\_tr\_pre))\) 。它还允许你传递优化约束。在这里，你希望所有的 \(\omega\) 都是非负的，并且 \(np.sum(\omega)==1\) 。

在下面的代码中，我按照 \(scikit-learn\) 的样板来构建合成控制模型。要做到这一点，你可以扩展 \(BaseEstimator\) 和 \(RegressorMixin\)，并定义一个 \(.fit\) 和一个 \(.predict\) 方法。其余的代码，比如 \(check\_X\_y\)、\(check\_array\) 和 \(check\_is\_fitted\)，只是一些标准检查，你无需担心：

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.utils.validation import (check_X_y, check_array,check_is_fitted)
import cvxpy as cp

class SyntheticControl(BaseEstimator, RegressorMixin):

    def __init__(self,):
        pass

    def fit(self, y_pre_co, y_pre_tr):

        y_pre_co, y_pre_tr = check_X_y(y_pre_co, y_pre_tr)
        # 待求解的权重变量
        w = cp.Variable(y_pre_co.shape[1])
        
        objective = cp.Minimize(cp.sum_squares(y_pre_co@w - y_pre_tr))
        constraints = [cp.sum(w) == 1, w >= 0]
        
        problem = cp.Problem(objective, constraints)
        
        self.loss_ = problem.solve(verbose=False)
        self.w_ = w.value
        
        self.is_fitted_ = True
        return self
        
        
    def predict(self, y_co):

        check_is_fitted(self)
        y_co = check_array(y_co)
        
        return y_co @ self.w_

在定义了 \(SyntheticControl\) 类之后，你可以像之前使用 \(LinearRegression\) 那样来使用它。注意，我存储了估计模型的最终损失。如果你想要在模型中纳入协变量，这会派上用场，正如你很快会看到的那样。此外，模型拟合后，你可以通过 \(.w\_\)来获取权重：

model = SyntheticControl()
model.fit(y_pre_co, y_pre_tr.mean(axis=1))

SyntheticControl()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

注意另一个有趣的点：你施加的凸性约束为优化问题提供了一个稀疏解。只有少数几个城市被用来构建最终的合成控制组。从这一点往后，就和之前完全一样了。你可以对整个数据集进行预测，以获得 \(E[Y_0|D = 1]\) 的合成控制估计值，并利用该估计值得到 \(ATT\) 的估计值：

将这个平均处理效应图与你之前使用无约束回归得到的图进行比较。现在，训练（处理前）误差略大，但 \(ATT\) 没那么嘈杂了。这就是正则化在起作用。

---

正则化回归

一旦你意识到合成控制法其实就是水平回归，你就可以找到其他对其进行正则化的方法。例如，你可以使用 \(Lasso\) 回归或 \(Ridge\) 回归。不过，不允许出现负权重仍然是很有意义的，尤其是当个体是不同地域时，因为这些地域的结果往往呈正相关。

---

---

合成控制法的假设

就像双重差分法一样，在使用合成控制法时，你也必须假设不存在对处理的预期，且不存在溢出效应。这两种方法之间的主要区别在于对潜在结果模型的参数假设。在双重差分法中，你必须假设处理组的 \(Y_0\) 趋势与对照组的 \(Y_0\) 趋势是平行的。另一方面，合成控制法对该潜在结果允许一个更灵活（但也更复杂）的模型：向量自回归模型或线性因子模型。

对于因子模型：

\(Y_{0it} = \lambda_{it}\hat{\epsilon}_{mt}f_t + e_{it},\)

如果你令 \(\lambda_1 = 1\)、\(f_1 = \beta_t\)，且 \(\lambda_1 = \alpha_i\)、\(f_1 = 1\)，你会发现它变成了双向固定效应模型 \(Y_{0it} = \alpha_i + \beta_t + e_{it}\) 的一种推广形式。

在《使用合成控制法：可行性、数据要求与方法学方面》这篇论文中，\(Alberto \space Abadie\) 指出，如果潜在结果遵循向量自回归模型或线性因子模型，并且合成控制组与处理组个体相匹配，那么合成控制法会对平均处理效应产生无偏估计。在实际情况中，合成控制组只能近似模拟处理组个体，因此存在一定偏差是可以预料的。

---

含协变量的合成控制法

通常，合成控制法仅使用控制组个体的处理前结果变量作为特征来预测 \(\overline{Y}_{tr}\)。这是因为这些特征往往是你可用的最具预测性的特征。不过，如果你认为某些额外的协变量具有良好的预测能力，你可能希望将它们纳入模型。不过，这种情况相当少见，所以如果时间紧张，你可以跳过这一部分。

假设你以某种方式获取了主要竞争对手的每日下载量数据，并且你也已经按市场规模对其进行了标准化，即 \(comp\_down\_pct\) 。你认为这个协变量是 \(\overline{Y}_{tr}\) 的良好预测因子，所以你希望将其纳入合成控制模型：

df_norm_cov = pd.read_csv("./data/online_mkt_cov.csv").astype({"date":"datetime64[ns]"})

df_norm_cov.head()

##    app_download  population  ... app_download_pct comp_download_pct
## 0        3066.0    12396372  ...         0.024733          0.026280
## 1        2701.0    12396372  ...         0.021789          0.023925
## 2        1927.0    12396372  ...         0.015545          0.018930
## 3        1451.0    12396372  ...         0.011705          0.015858
## 4        1248.0    12396372  ...         0.010067          0.014548
## 
## [5 rows x 9 columns]

用数学符号表示的话，你想要构建一个合成控制组，使得权重 \(w_i\) 不仅与 \(y_{co}\) 相乘，还与这个额外的协变量 \(x_{co}\) 相乘，以此来近似 \(\overline{Y}_{tr}\)。这里的问题在于，\(x_{co}\) 和 \(y_{co}\) 可能处于完全不同的尺度，或者其中一个可能比另一个更具预测性，这就是为什么在求解 \(SC\) 优化问题之前，你需要将每个协变量（包括 \(y_{co}\)）乘以一个缩放因子 \(v\)。为了考虑到这一点，你可以用协变量 \(X\) 来重写目标函数，将 \(y_{co}\) 仅仅视为另一个协变量：

\(\widehat{\omega}^{sc} = \underset{\omega}{\text{argmin}} \ \left\Vert \overline{\mathbf{y}}_{pre,tr} - \sum v_k^*\mathbf{X}_{k,pre,co}\omega_{co} \right\Vert^2\)

约束条件：\(\sum \omega_i = 1\) 且 \(\omega_i > 0 \ \forall \ i\)

然而，这个目标函数并没有告诉你如何找到最优的 \(v\)。为了做到这一点，你必须将整个合成控制过程纳入到另一个优化目标中。如果这听起来很复杂，别担心。在代码中理解它要容易得多。首先，让我们为协变量 \(comp\_download\_pct\)、\(app\_download\_pct\) 以及 \(Y_{pre,co}\) 创建 \(X\) 矩阵：

from toolz import partial

reshaper = partial(reshape_sc_data,df=df_norm_cov,geo_col="city",time_col="date",tr_geos=treated,tr_start=str(tr_period))

y_pre_co, y_pre_tr, y_post_co, y_post_tr = reshaper(y_col="app_download_pct")
x_pre_co, _, x_post_co, _ = reshaper(y_col="comp_download_pct")

接下来，让我们编写一个函数，当给定一个 \(v\) 的列表（每个协变量对应一个 \(v\)）时，返回合成控制权重和优化损失。记住，你可以通过拟合好的 \(SyntheticControl\) 模型的 \(.loss\_\) 来获取目标损失。

为了检查这是否有效，你可以传入 \([1,0]\) 作为 \(v\) 的值，并传入 \([y\_pre\_co, x\_pre\_co]\) 作为协变量列表。在这种情况下，由于额外的协变量权重为 \(0\)，你应该得到原来的合成控制结果：

def find_w_given_vs(vs, x_co_list, y_tr_pre):
    X_times_v = sum([x*v for x, v in zip(x_co_list, vs)])
    
    model = SyntheticControl()
    model.fit(X_times_v, y_tr_pre)
    
    return {"loss": model.loss_, "w": model.w_} 

find_w_given_vs([1, 0],[y_pre_co, x_pre_co],y_pre_tr.mean(axis=1)).get("w").round(3)

## array([-0.   , -0.   , -0.   , -0.   , -0.   , -0.   ,  0.076,  0.037,
##         0.083,  0.01 , -0.   , -0.   , -0.   , -0.   , -0.   , -0.   ,
##         0.061,  0.123,  0.008,  0.074, -0.   ,  0.   , -0.   , -0.   ,
##        -0.   , -0.   , -0.   , -0.   , -0.   ,  0.   , -0.   ,  0.092,
##        -0.   , -0.   ,  0.   ,  0.046,  0.089,  0.   ,  0.067,  0.061,
##         0.   , -0.   , -0.   ,  0.088,  0.   ,  0.086, -0.   ])

最后，你可以将 \(find\_w\_given\_vs\) 封装到一个只接收 \(v\) 的数组并返回优化损失的函数中。然后，你可以将这个函数传递给 \(scipy\) 的 \(minimize\) 函数，该函数会迭代地寻找最佳的 \(v\) 并返回给你：

from scipy.optimize import minimize

def v_loss(vs):
    return find_w_given_vs(vs, [y_pre_co, x_pre_co],y_pre_tr.mean(axis=1)).get("loss")

v_solution = minimize(v_loss, [0, 0], method='L-BFGS-B')
v_solution.x

## array([1.80162695, 0.00311651])

有了最优的 \(v\) 值，你可以回到 \(find\_w\_given\_vs\) 来获取考虑了协变量的合成控制权重。不过，有一点需要注意，最终的解决方案与不包含协变量的方案并没有太大差异。这并不奇怪，因为针对 \(comp\_download\_pct\) 协变量的最优 \(v\) 是一个非常小的数，而且它的尺度与 \(app\_download\_pct\) 相比并没有大很多：

w_cov = find_w_given_vs(v_solution.x,[y_pre_co, x_pre_co],y_pre_tr.mean(axis=1)).get("w").round(3)

w_cov

## array([-0.   , -0.   ,  0.   , -0.   , -0.   , -0.   ,  0.086,  0.   ,
##         0.034, -0.   , -0.   ,  0.035, -0.   ,  0.   , -0.   ,  0.   ,
##         0.016,  0.051,  0.031,  0.01 , -0.   ,  0.   ,  0.   ,  0.053,
##        -0.   ,  0.   , -0.   , -0.   ,  0.   ,  0.   , -0.   ,  0.048,
##         0.084,  0.007,  0.   ,  0.   ,  0.147,  0.   ,  0.023,  0.009,
##        -0.   ,  0.16 ,  0.   ,  0.036,  0.   ,  0.169, -0.   ])

利用这些权重，你可以对 \(Y(0)\) 进行预测，并得到考虑了协变量的最终 \(ATT\) 估计值：

y0_hat = sum([x*v for x, v in zip([y_post_co, x_post_co], v_solution.x)]).dot(w_cov)
att = y_post_tr.mean(axis=1) - y0_hat

下面的图表展示了所得的 \(ATT\)，以及不包含协变量的 \(SC\) 所得到的 \(ATT\) 估计值。如你所见，两者非常相似：

尽管纳入协变量并不困难，但需要相当多的额外复杂性。出于这个原因，再加上 \(Y_{pre,co}\) 往往足以预测 \(Y_{tr}\) 这一事实，我通常不会费心去添加协变量。不过，也许你能找到非常有预测性的特征，从而证明添加协变量是合理的。

---

通用水平回归

添加协变量的一种更简单的方法是，只需将任何你认为有价值的额外时间序列作为一列拼接到 \(\mathbf{Y}_{pre,co}\) 中。这相当于在水平回归中添加额外的协变量：

\(\left[ \mathbf{Y}_{pre,co} \vert \mathbf{X}_{pre,co} \right] \omega\)

从严格意义上来说，这并不是合成控制，因为你现在是用控制组个体和那些额外的时间序列来估计 \(E[Y(0)_{tr}]\)。因此，最终你得到的权重不仅是针对个体的，也是针对你所拼接的每一个额外列的。

---

合成控制法的去偏

与强大的机器学习模型非常相似，这些预测技术容易出现过拟合，尤其是在处理前时期的数量 \(T_{pre}\) 较小时。即使对标准合成控制法施加的约束也无法完全解决这个问题。因此，合成控制法是有偏的。为了理解这一点，让我们将 \(ATT\) 重新定义为干预后时期内的时间平均值：

\(ATT = \frac{1}{T1} \sum_{t=T0+1}^{T} \left( \overline{Y}_{1t} - \overline{Y}_{0t} \right)\)

其中，\(T0\) 和 \(T1\) 分别是干预前和干预后时期的时长，\(\overline{Y}_{dt}\) 是处理组个体的平均潜在结果。这只是将每个单独干预后时期的 \(ATT\) 平均为一个单一数值，使其更易于处理。现在，为了检查 \(SC\) 方法中的偏差，你可以将这个单一数值与其估计值进行比较。\(图9-2\) 展示了该偏差。

我正在按照合成控制法的设定模拟大量数据 —— 其中处理组个体是一些控制组个体的加权组合。这里，\(N = 16\)，\(T_{pre} = 15\)，所以在水平回归中列数多于行数。而且，真实的 \(ATT\) 为 \(0\)。尽管如此，用合成控制法得到的 \(ATT\) 估计值的分布并非以 \(0\) 为中心，这表明它确实存在偏差。

from statsmodels.tsa.arima_process import ArmaProcess

def sim_sc_data(N=16, T0=15, T1=28):
    Y_co = np.random.normal(
        np.array([np.concatenate([np.ones(3), np.zeros(N-3)])*2 for _ in range(T0+T1)])
    )
    
    w = np.concatenate([np.ones(3), np.zeros(N-3)])*1/3
    Y_tr = Y_co.dot(w) + ArmaProcess(0.8).generate_sample(T0+T1)
    
    model = SyntheticControl().fit(
        y_pre_co=Y_co[:T0],
        y_pre_tr=Y_tr[:T0].flatten()
    )
        
    return (Y_tr[T0:].flatten() - model.predict(Y_co[T0:])).mean()


np.random.seed(123)

atts = [sim_sc_data() for _ in range(500)]

plt.figure(figsize=(10,4))

## <Figure size 1000x400 with 0 Axes>

plt.hist(atts, alpha=0.2, density=True, bins=20)

## (array([0.05774333, 0.11548667, 0.25022111, 0.53893778, 0.73141555,
##        0.86615   , 1.00088444, 1.38584   , 1.00088444, 0.71216778,
##        0.94314111, 0.59668111, 0.42345111, 0.42345111, 0.23097333,
##        0.11548667, 0.11548667, 0.03849556, 0.01924778, 0.05774333]), array([-0.61660423, -0.51269613, -0.40878804, -0.30487994, -0.20097184,
##        -0.09706374,  0.00684436,  0.11075246,  0.21466056,  0.31856866,
##         0.42247676,  0.52638486,  0.63029296,  0.73420106,  0.83810916,
##         0.94201726,  1.04592536,  1.14983346,  1.25374156,  1.35764965,
##         1.46155775]), <BarContainer object of 20 artists>)

plt.vlines(np.mean(atts), 0, 1.5, label="Simulation Avg.")

## <matplotlib.collections.LineCollection object at 0x000001CFABB838E0>

plt.vlines(0, 0, 1.5, label="Correct ATT", ls="dashed")

## <matplotlib.collections.LineCollection object at 0x000001CFABBAD030>

plt.title("500 SC Simulations")

## Text(0.5, 1.0, '500 SC Simulations')

plt.legend()

## <matplotlib.legend.Legend object at 0x000001CFABB82FB0>

幸运的是，在第 \(7\) 章中，当你学习双重 / 去偏机器学习时，已经了解了如何解决过拟合偏差。答案在于交叉拟合。思路是将干预前时期划分为 \(K\) 个块，每个块的大小为 \(min\{T_{pre}/K, T_{post}\}\)（这个 \(min\) 函数的原因很快就会清楚。现在，先相信我）。然后，将每个块视为留出集，并在 \(Y_{pre,co}^{-k}\) 和 \(Y_{pre,tr}^{-k}\) 上拟合合成控制模型，其中 \(-k\) 表示从训练中排除该块。这一步会得到权重 \(\widehat{\omega}^{-k}\)。接下来，你将使用这些权重，利用留出数据 \(Y_{pre,co}^{k}\) 进行样本外预测。预测值与留出数据中观测值的平均差异就是对偏差的估计：

\(\widehat{\text{Bias}}^{k} = avg(\mathbf{Y}_{pre,tr}^{k} - \mathbf{Y}_{pre,co}^{k}\widehat{\omega}^{-k})\)

这意味着你可以用它来调整 \(ATT\) 的估计值：

\(\widehat{ATT}^{k} = \mathbf{Y}_{post,tr} - \mathbf{Y}_{post,co}\widehat{\omega}^{-k} - \widehat{\text{Bias}}^{k}\)

注意，这会给你 \(K\) 个不同的 \(ATT\)。你可以对它们取平均，得到最终的 \(ATT\) 估计值。

现在，让我们把这个用 \(Python\) 代码实现。这里最棘手的部分是定义这些块。所以，为了更清楚地说明，让我们看一个简单的例子，如 \(图9-3\) 所示。

假设你有 \(5\) 个干预前时期和 \(3\) 个干预后时期，并且你想要构建 \(K = 2\) 个块。块的大小是 \(2.5\)，这不是一个整数，所以你必须向下取整为 \(2\)。这意味着你要从干预前时期中取出 \(2\) 个大小为 \(2\) 的块。\(2×2\) 会给你 \(4\) 个时期，但你有 \(5\) 个时期。所以，我选择从干预前时期的末尾取出这些块，这会导致第一个时期永远不会被移除。这相当随意，但对整个过程不会有太大影响。你也可以选择修剪干预前时期，使其能被 \(K\) 整除。

最后，你将使用权重和偏差估计值来对干预后时期的 \(ATT\) 进行估计。

尽管描述起来有点复杂，但用 \(Numpy\) 获取这些块是相当容易的。首先，你要对干预前时期的末尾进行索引，通过 \(y\_pre\_tr.index[-K*block\_size:]\) 来获取恰好包含 \(K\) 个块的索引。然后，你可以使用 \(np.split\) 将这些索引拆分成 \(K\) 个块。这会返回一个数组，其中有 \(K\) 行，每一行包含在每次迭代中你想要移除的索引。一旦你有了这些块，就可以遍历它们，拟合 \(SC\) 模型，估计偏差以及干预后时期的 \(ATT\)。结果可以存储在一个数据框中，以便于展示：

def debiased_sc_atts(y_pre_co, y_pre_tr, y_post_co, y_post_tr, K=3):
    # 计算每个验证块的大小：取干预前数据可分块数量与干预后数据长度的最小值
    block_size = int(min(np.floor(len(y_pre_tr)/K), len(y_post_tr)))
    # 从干预前数据的末尾截取 K*block_size 个时间点，并分成 K 个块
    blocks = np.split(y_pre_tr.index[-K*block_size:], K)
    
    def fold_effect(hold_out):
        # 1. 用排除当前验证块的干预前数据拟合合成控制模型
        model = SyntheticControl()
        model.fit(
            y_pre_co.drop(hold_out),  # 控制组数据排除验证块时间
            y_pre_tr.drop(hold_out)   # 处理组数据排除验证块时间
        )
        
        # 2. 估计偏差：验证块时间内，处理组实际值与合成控制预测值的平均差异
        bias_hat = np.mean(y_pre_tr.loc[hold_out]  # 处理组验证块实际值
                           - model.predict(y_pre_co.loc[hold_out]))  # 合成控制对验证块的预测值
        
        # 3. 用拟合好的模型预测干预后处理组的反事实结果（无干预时的预测值）
        y0_hat = model.predict(y_post_co)
        
        # 4. 计算去偏后的处理效应：(实际结果 - 反事实预测) - 估计的偏差
        return (y_post_tr - y0_hat) - bias_hat
    
    
    # 对每个验证块计算去偏效应，结果转置后返回DataFrame
    return pd.DataFrame([fold_effect(block) for block in blocks]).T

要将这个函数应用到（已经透视过的）营销数据上，你只需要记住对处理组个体进行平均。

结果是一个包含所有平均处理效应（ATT）估计值的数据框。它有 K 列，且每个干预后时期对应一行：

deb_atts = debiased_sc_atts(y_pre_co,y_pre_tr.mean(axis=1),y_post_co,y_post_tr.mean(axis=1),K=3)
deb_atts.head()

##                    0         1         2
## date                                    
## 2022-05-01  0.003314  0.002475  0.003162
## 2022-05-02  0.003544  0.002844  0.003368
## 2022-05-03  0.004644  0.003698  0.004759
## 2022-05-04  0.004706  0.002866  0.003603
## 2022-05-05  0.000134 -0.000541  0.000173

要得到每个干预后时期的最终 \(ATT\) 估计值，你可以对列进行平均，即 \(deb\_atts.mean(axis=1)\) ；或者，如果你想要整个时期的单一 \(ATT\) 值，只需对所有值取平均：\(deb\_atts.mean(axis=1).mean()\) 。将去偏后的 \(ATT\) 与你之前得到的 \(SC\) 的 \(ATT\) 绘在一起，也能显示出，在大部分情况下，去偏提高了 \(ATT\) 估计值，尽管幅度不大：

# 实例化合成控制模型
model_biased = SyntheticControl()
# 用控制组和处理组干预前的数据拟合模型（处理组取均值，可能是多处理组的平均）
model_biased.fit(y_pre_co, y_pre_tr.mean(axis=1))

## SyntheticControl()

# 计算有偏的处理效应：处理组干预后实际均值 - 合成控制组的预测值
att_biased = y_post_tr.mean(axis=1) - model_biased.predict(y_post_co)

# 创建画布
plt.figure(figsize=(10,4))

## <Figure size 1000x400 with 0 Axes>

# 绘制去偏后的ATT（虚线）
plt.plot(att, label="ATT Debiased", ls="-.")

## [<matplotlib.lines.Line2D object at 0x000001CFABC06980>]

# 绘制有偏的ATT（实线，默认）
plt.plot(att_biased, label="ATT Biased")

## [<matplotlib.lines.Line2D object at 0x000001CFABC2BE20>]

# 添加图例、参考线、标签和旋转x轴刻度
plt.legend()  # 显示图例区分两条线

## <matplotlib.legend.Legend object at 0x000001CFABC49A50>

plt.hlines(0, att.index.min(), att.index.max())  # 添加y=0的水平线（基准线，效应为0的参考）

## <matplotlib.collections.LineCollection object at 0x000001CFABC48FA0>

plt.ylabel("ATT")  # y轴标签：处理效应值

## Text(0, 0.5, 'ATT')

plt.xticks(rotation=45);  # x轴刻度旋转45度，避免重叠（通常x轴是时间）

## (array([19113., 19120., 19127., 19134., 19144., 19151., 19158., 19165.,
##        19174.]), [Text(19113.0, 0, '2022-05-01'), Text(19120.0, 0, '2022-05-08'), Text(19127.0, 0, '2022-05-15'), Text(19134.0, 0, '2022-05-22'), Text(19144.0, 0, '2022-06-01'), Text(19151.0, 0, '2022-06-08'), Text(19158.0, 0, '2022-06-15'), Text(19165.0, 0, '2022-06-22'), Text(19174.0, 0, '2022-07-01')])

在你的营销数据中很难看出差异，但为了说明去偏的重要性，我可以重新进行之前的模拟，不过现在使用去偏过程。现在，模拟得到的 \(ATT\) 分布的均值为 \(0\)，正如它本应的那样：

def sim_deb_sc_data(N=16, T0=15, T1=28):
    Y_co = np.random.normal(np.array([np.concatenate([np.ones(3), np.zeros(N-3)])*2 for _ in range(T0+T1)]))
    w = np.concatenate([np.ones(3), np.zeros(N-3)])*1/3
    Y_tr = Y_co.dot(w) + ArmaProcess(0.8).generate_sample(T0+T1)
    df_co = pd.DataFrame(Y_co)
    df_tr = pd.DataFrame(Y_tr)
    atts = debiased_sc_atts(df_co.iloc[:T0],df_tr.iloc[:T0].mean(axis=1),df_co.iloc[T0:],df_tr.iloc[T0:].mean(axis=1), K=2)
    return atts.mean(axis=0).mean()


np.random.seed(123)
atts = [sim_deb_sc_data() for _ in range(500)]

plt.figure(figsize=(10,4))

## <Figure size 1000x400 with 0 Axes>

plt.hist(atts, alpha=0.2, density=True, bins=20)

## (array([0.01193041, 0.        , 0.        , 0.02386081, 0.11930407,
##        0.16702569, 0.25053854, 0.5368683 , 0.73968522, 0.78740684,
##        0.94250213, 0.68003318, 0.65617237, 0.51300749, 0.22667773,
##        0.10737366, 0.13123447, 0.01193041, 0.03579122, 0.02386081]), array([-1.75077237, -1.58313349, -1.41549461, -1.24785574, -1.08021686,
##        -0.91257798, -0.7449391 , -0.57730022, -0.40966134, -0.24202247,
##        -0.07438359,  0.09325529,  0.26089417,  0.42853305,  0.59617193,
##         0.7638108 ,  0.93144968,  1.09908856,  1.26672744,  1.43436632,
##         1.6020052 ]), <BarContainer object of 20 artists>)

plt.vlines(np.mean(atts), 0, 1.0, label="Sample Avg.")

## <matplotlib.collections.LineCollection object at 0x000001CFABCFF850>

plt.vlines(0, 0, 1, label="Correct ATT", ls="dashed")

## <matplotlib.collections.LineCollection object at 0x000001CFABCFEE30>

plt.title("500 SC (Debiased) Simulations")

## Text(0.5, 1.0, '500 SC (Debiased) Simulations')

plt.legend()

## <matplotlib.legend.Legend object at 0x000001CFABD2AA40>

推断

去偏过程本身就很有用，但它还有第二个有趣的原因，那就是可以为合成控制的 \(ATT\) 估计值设定置信区间。用合成控制进行推断被证明是一项艰巨的任务，主要是因为通常控制组个体非常少，甚至只有一个。你在第 \(8\) 章中学到的 块自助抽样 在这里行不通，因为很多自助抽样样本会排除所有处理组个体，导致 \(ATT\) 无定义。

合成控制法的推断是一个活跃的研究领域，并且发展迅速。在过去的几年里，已经提出了许多方法。它们中的大多数依赖于对时间维度进行置换，因为对个体进行 自助法似乎存在问题。在这里，我选择了其中一种我认为实现起来相当简单且计算效率很高的方法。尤其是如果你已经处理过去偏部分的话，因为它将去偏部分作为起点。只是作为回顾，要记得去偏为 \(K\) 个折中的每一个以及每个干预后时间周期都提供了一个 \(ATT\) 估计值，这些估计值在下面的数据框中以列的形式呈现：

deb_atts.head()

##                    0         1         2
## date                                    
## 2022-05-01  0.003314  0.002475  0.003162
## 2022-05-02  0.003544  0.002844  0.003368
## 2022-05-03  0.004644  0.003698  0.004759
## 2022-05-04  0.004706  0.002866  0.003603
## 2022-05-05  0.000134 -0.000541  0.000173

现在，假设你想要为干预后时期的整体 \(ATT\) 估计值设定一个置信区间。要做到这一点，首先你需要的是 \(\widehat{ATT}\) 本身。你可以对这个数据框的行进行平均，这会为 \(K\) 个折中中的每一个得到一个单一的 \(ATT\)。然后，你可以对这些值取平均：

atts_k = deb_atts.mean(axis=0).values
att = np.mean(atts_k)

print("atts_k:", atts_k)

## atts_k: [0.00414872 0.00260513 0.00317687]

print("ATT:", att)

## ATT: 0.003310235951728967

现在，来看推断部分。这里的思路是基于每个 \(ATT^{k}\) 构建一个标准误估计值：

\(\widehat{\sigma} = \sqrt{1 + \frac{BlockSize * K}{T_{post}}} * \sqrt{\frac{1}{K - 1} \sum_{k=1}^{K} (ATT^{k} - ATT)}\)

\(\widehat{SE} = \widehat{\sigma} / \sqrt{K}\)

在编写代码时，你只需要注意使用样本标准差，这意味着要向 \(np.std\) 传递 \(ddof=1\)：

K = len(atts_k)  # 交叉验证的折数（等于去偏效应的数量）
T0 = len(y_pre_co)  # 干预前的时间点总数
T1 = len(y_post_co)  # 干预后的时间点总数
# 每个验证块的大小：取"干预前可分块大小"和"干预后时间长度"的最小值（与去偏函数中的逻辑一致）
block_size = min(np.floor(T0/K), T1)
# 标准误计算公式
se_hat = np.sqrt(1 + ((K * block_size) / T1)) * np.std(atts_k, ddof=1) / np.sqrt(K)
print("SE:", se_hat)

## SE: 0.0006345693279981852

有了这个标准误，你可以构建检验统计量 \(\widehat{ATT} / \widehat{SE}\)，在原假设 \(H0: ATT = 0\) 下，该统计量渐近服从自由度为 \(K - 1\) 的 t 分布。这意味着你可以利用它，通过 \(t\) 分布来构建置信区间。例如，下面是构建 \(90\%\) 置信区间（\(\alpha = 0.1\)）的方法：

from scipy.stats import t

alpha = 0.1
[att - t.ppf(1-alpha/2, K-1)*se_hat,att + t.ppf(1-alpha/2, K-1)*se_hat]

## [0.001457302664238376, 0.005163169239219558]

你可能会看到标准误公式分母中的 \(K\)，并想把它设为一个很大的数。然而，天下没有免费的午餐。\(K\) 值越大，置信区间就越小，但代价是这些区间的覆盖率会降低。对于大的 \(K\) 值，\(1 - \alpha\) 的置信区间包含真实 \(ATT\) 的次数会少于 \(1 - \alpha\)，尤其是在处理前时期的数量较少时。在这种情况下，\(K\) 合理的选择是 \(3\)。当 \(T_0\) 与 \(N\) 相比非常大时，你可以尝试更大的 \(K\) 值来缩短置信区间的长度。

另一个重要的点是，这种方法不适用于处理效应的轨迹。也就是说，它不能用于逐期推断，因为其理论要求 \(T0\) 和 \(T1\) 都相对较大。

---

另见

这种推断方法是由 \(Victor \space Chernozhukov\) 等人在论文《A T-Test for Synthetic Controls》中提出的。如果你想要进行逐期推断，同一作者还有一篇补充论文，提出了用于合成控制的共形推断方法：《An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls》。

---

合成双重差分法

本章结束之际，我想为你提供关于合成控制法的另一种视角，即它与双重差分法的关联。通过这样做，你也将学习如何把这两种方法整合为一个单一的合成双重差分\(（SDID）\)估计量。这里的思路十分简单。首先，构建一个合成控制组。然后，在双重差分法的设定中，将其用作控制单元。最终的结果比两种方法简单相加要有趣得多。首先，双重差分法所需的平行趋势假设会变得更具合理性，因为你正在为 \(E[Y(0)_{ti}|D = 1]\) 构建一个合成控制组。其次，由于你使用了双重差分法，合成控制组可以只专注于捕捉处理组的趋势，因为它的 \(Y(0)\) 水平可以不同。不过，首先，在深入探讨合成双重差分法之前，让我们先回顾一些双重差分法的理论。

双重差分法回顾

在其标准形式中，假设有一个控制组（从未接受处理）和一个处理组（在同一时期全部接受处理），你可以用双向固定效应来表示双重差分法：

\(\widehat{\tau}^{did} = \underset{\mu, \alpha_i, \beta_t, \tau}{\operatorname{argmin}} \left\{ \sum_{n=1}^{N} \sum_{t=1}^{T} \left( Y_{it} - \mu - \alpha_i -\beta_t - \tau W_{it} \right)^2 \right\},\)

其中，\(\tau\) 是你所关注的 \(ATT\)，\(\alpha_i\) 是个体固定效应，\(\beta_t\) 是时间固定效应。在这种表述中，个体效应捕捉了每个个体截距的差异，而时间效应捕捉了处理组和控制组共同的总体趋势。双重差分法背后的主要假设是，处理组和未处理组具有相同的 \(Y_0\) 趋势：

\(\Delta Y(d)_i \perp D\)

再论合成控制法

接下来，让我们看看如何将合成控制估计量重新表述为类似于前面 \(DID\) 的形式。非常有趣的是，你可以将 \(SC\) 估计量写成求解以下优化问题的形式：

\(\widehat{\tau}^{sc} = \underset{\beta, \tau}{\operatorname{argmin}} \left\{ \sum_{n=1}^{N} \sum_{t=1}^{T} \left( Y_{it} - \beta_t - \tau W_{it} \right)^2 \widehat{\omega}_i^{sc} \right\},\)

其中，控制单元的权重 \(\widehat{\omega}_i^{sc}\) 是通过优化本章开头你看到的合成控制目标函数得到的。由于前面 \(SC\) 目标函数的表述是针对所有单元而不仅仅是控制单元定义的，你还需要考虑处理单元的权重。在这里，因为你关注的是 \(ATT\)，所以它们只是\(N_{tr}/N\)（均匀加权）。

为了验证这个新表述确实与你之前所学的表述等价，让我们对两者进行比较。首先，按照你到目前为止的做法估计合成控制，并计算 \(ATT\) ：

sc_model = SyntheticControl()
sc_model.fit(y_pre_co, y_pre_tr.mean(axis=1))

## SyntheticControl()

(y_post_tr.mean(axis=1) - sc_model.predict(y_post_co)).mean()

## 0.003327040979396121

接下来，在进行矩阵重塑之前，让我们把这些合成控制权重添加到原始的营销数据框中。要做到这一点，你可以创建一个单元权重数据框，将每个控制城市映射到其对应的权重：

unit_w = pd.DataFrame(zip(y_pre_co.columns, sc_model.w_),columns=["city", "unit_weight"])
unit_w.head()

##                    city  unit_weight
## 0            ananindeua    -0.000008
## 1  aparecida_de_goiania    -0.000001
## 2               aracaju    -0.000008
## 3                 belem    -0.000012
## 4          belford_roxo    -0.000006

然后，你可以使用城市作为键，将这个单元权重数据框合并到原始的营销数据框中。这会使处理组单元的权重为 \(NaN\)（缺失值）。你可以用处理组虚拟变量的平均值（即\(N_{tr}/N\)）来填充这些缺失值。

我也会借此机会，通过将\(D_i\)乘以\(Post_t\)来创建\(W_{it}\)变量：

df_with_w = (df_norm
             .assign(tr_post = lambda d: d["post"]*d["treated"])
             .merge(unit_w, on=["city"], how="left")
             .fillna({"unit_weight": df_norm["treated"].mean()}))
             
df_with_w.head()

##    app_download  population       city  ... app_download_pct tr_post  unit_weight
## 0        3066.0    12396372  sao_paulo  ...         0.024733       0         0.06
## 1        2701.0    12396372  sao_paulo  ...         0.021789       0         0.06
## 2        1927.0    12396372  sao_paulo  ...         0.015545       0         0.06
## 3        1451.0    12396372  sao_paulo  ...         0.011705       0         0.06
## 4        1248.0    12396372  sao_paulo  ...         0.010067       0         0.06
## 
## [5 rows x 10 columns]

最后，你可以按照之前提到的另一种合成控制公式，进行带有时间固定效应的加权普通最小二乘法\(（WLS）\)回归。只是要注意删除权重非常小的行，否则在尝试运行这个回归时，你可能会遇到一些错误：

mod = smf.wls(
    "app_download_pct ~ tr_post + C(date)",
    # 指定权重列 unit_weight，即每个观测值的权重。WLS 会根据这些权重调整回归，权重越高的样本对回归结果影响越大
    data=df_with_w.query("unit_weight>=1e-10"),
    weights=df_with_w.query("unit_weight>=1e-10")["unit_weight"]
)

mod.fit().params["tr_post"]

## 0.0033293800074678482

事实上，这里得到的 \(ATT\) 与你之前得到的完全相同，这表明两种合成控制公式是等价的。但更重要的是，新的 \(SC\) 公式更容易与 \(DID\) 的双向固定效应公式进行比较。首先，看起来合成控制法有时间固定效应，但没有个体固定效应。与此同时，双重差分法同时具有时间和个体固定效应，但没有个体权重。这表明可以将这两个模型合并，纳入来自合成控制法和双重差分法的元素：

\(\widehat{\tau}^{sdid} = \underset{\mu, \alpha, \beta, \tau}{\operatorname{argmin}} \left\{ \sum_{n=1}^{N} \sum_{t=1}^{T} \left( Y_{it} - \left( \mu + \alpha_i + \beta_t + \tau D_{it} \right) \right)^2 \widehat{\omega}_i \right\}\)

而且，既然都做到这一步了，为什么只对单元个体进行加权呢？时刻记住这里的最终目标：估计 \(E[Y_0 | Post = 1, Tr = 1]\)。单元权重的目的是利用控制单元来近似处理单元。但这里还有一个时间维度，这意味着你也可以在干预前时期使用权重，以更好地近似干预后时期。这会得到以下合成双重差分\(（SDID）\)公式：

\(\widehat{\tau}^{sdid} = \underset{\mu, \alpha, \beta, \tau}{\operatorname{argmin}} \left\{ \sum_{n=1}^{N} \sum_{t=1}^{T} \left( Y_{it} - \left( \mu + \alpha_i + \beta_t + \tau D_{it} \right) \right)^2 \widehat{\omega}_i \widehat{\lambda}_t \right\},\)

其中 \(\widehat{\lambda}_t\) 是时间权重。

估计时间权重

回想一下，为了得到单元权重，你是如何在干预前时期，将处理单元的平均结果对控制单元的结果进行回归的？

\(\begin{aligned} \widehat{\omega}_i^{sc} &= \underset{\omega}{\operatorname{argmin}} \left\| \overline{\mathbf{y}}_{pre,tr} - \mathbf{Y}_{pre,co} \omega_{co} \right\|^2 \\ \text{s.t } &\sum \omega_i = 1 \text{ and } \omega_i > 0 \ \forall \ i \end{aligned}\)

那么，要得到时间权重，你只需要对\(\mathbf{Y}_{pre,co}\)进行转置，然后将其对控制单元在干预后时期的平均结果进行回归：

\(\begin{aligned} \widehat{\lambda}_t^{sc} &= \underset{w}{\operatorname{argmin}} \left\| \overline{\mathbf{y}}_{pre,co}' - \mathbf{Y}_{pre,co}' \lambda_{pre} \right\|^2 \\ \text{s.t } &\sum \lambda_i = 1 \text{ and } \lambda_i > 0 \ \forall \ i \end{aligned}\)

但这里还有一个额外的问题。还记得 \(SC\) 不允许外推吗？如果结果存在某种趋势，这就会成为一个问题。要是那样的话，干预后时期的平均结果会比所有干预前时期的结果更高或更低，而你需要进行外推才能得到良好的拟合。因此，\(SDID\) 公式在寻找时间权重时，允许存在一个截距项偏移 \(\lambda_0\)：

\(\begin{aligned} \widehat{\lambda}_t^{sc} &= \underset{w}{\operatorname{argmin}} \left\| \overline{\mathbf{y}}_{pre,co}' - \left( \mathbf{Y}_{pre,co}' \lambda_{pre} + \lambda_0 \right) \right\|^2 \\ \text{s.t } &\sum \lambda_i = 1 \text{ and } \lambda_i > 0 \ \forall \ i \end{aligned}\)

幸运的是，修改 \(SyntheticControl\) 代码以选择性地拟合截距相当容易，可使用一个 \(fit\_itercept\) 参数。首先，你要创建一个截距列，如果 \(fit\_intercept=True\)，该列始终为 \(1\)，否则为 \(0\)。你可以利用 \(Python\) 中 \(True*1=1\) 这一特性。然后，你会将这一列前置到 \(y\_pre\_co\) 中，并在目标函数中使用它。另外，在构建约束条件时，你不希望包含截距。最后，你会去掉截距参数，只返回单元的权重。

from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.utils.validation import check_X_y, check_array, check_is_fitted
import cvxpy as cp

class SyntheticControl(BaseEstimator, RegressorMixin):

    def __init__(self, fit_intercept=False):
        # 控制是否允许合成控制组与处理组之间存在一个固定的偏移量
        self.fit_intercept = fit_intercept

    def fit(self, y_pre_co, y_pre_tr):
        y_pre_co, y_pre_tr = check_X_y(y_pre_co, y_pre_tr)
        
        intercept = np.ones((y_pre_co.shape[0], 1))*self.fit_intercept
        # 约束条件：
        # 1. 控制组权重之和为1（截距项权重不参与求和）
        # 2. 控制组权重非负（确保是合理的加权组合）
        X = np.concatenate([intercept, y_pre_co], axis=1)
        w = cp.Variable(X.shape[1])
        
        # 目标函数：最小化合成组与处理组的平方误差和
        objective = cp.Minimize(cp.sum_squares(X@w - y_pre_tr))
        constraints = [cp.sum(w[1:]) == 1, w[1:] >= 0]
        # 求解优化问题
        problem = cp.Problem(objective, constraints)
        
        self.loss_ = problem.solve(verbose=False)
        # 只保留控制组的权重（去掉截距项权重）
        self.w_ = w.value[1:]
        
        self.is_fitted_ = True
        return self
        
        
    def predict(self, y_co):

        check_is_fitted(self)
        y_co = check_array(y_co)
        
        return y_co @ self.w_

一旦你处理完那个（问题 / 部分），就可以继续估计时间权重了：

time_sc = SyntheticControl(fit_intercept=True)
time_sc.fit(y_pre_co.T,y_post_co.mean(axis=0))

## SyntheticControl(fit_intercept=True)

time_w = pd.DataFrame(zip(y_pre_co.index, time_sc.w_),columns=["date", "time_weight"])
time_w.tail()

##          date  time_weight
## 56 2022-04-26    -0.000011
## 57 2022-04-27     0.071965
## 58 2022-04-28    -0.000002
## 59 2022-04-29     0.078350
## 60 2022-04-30     0.000002

我也已经把权重存储在一个数据框中了，你之后会用到它。你也可以将这些权重绘制成图，看看干预前的时间段是如何被用来重构干预后时期控制组的平均结果的：

plt.figure(figsize=(10,4))

## <Figure size 1000x400 with 0 Axes>

plt.bar(time_w["date"], time_w["time_weight"])

## <BarContainer object of 61 artists>

plt.ylabel("Time Weights")

## Text(0, 0.5, 'Time Weights')

plt.xticks(rotation=45) 

## (array([19052., 19059., 19066., 19073., 19083., 19090., 19097., 19104.,
##        19113.]), [Text(19052.0, 0, '2022-03-01'), Text(19059.0, 0, '2022-03-08'), Text(19066.0, 0, '2022-03-15'), Text(19073.0, 0, '2022-03-22'), Text(19083.0, 0, '2022-04-01'), Text(19090.0, 0, '2022-04-08'), Text(19097.0, 0, '2022-04-15'), Text(19104.0, 0, '2022-04-22'), Text(19113.0, 0, '2022-05-01')])

合成控制法与双重差分法

好的。现在你已经有了干预前时期的权重以及所有单元的权重。剩下要做的就是把这些部分整合到最终的估计量中。你可以从之前定义的名为 \(df\_with\_w\) 的数据框开始，然后用 \(date\) 作为键，合并进时间权重数据框。由于 \(time\_w\) 只有干预前时期的权重，你需要用 \(T_{post}/T\)（同样对它们进行均匀加权）来填充干预后时期的时间权重。最后，将\(\lambda_t\)和\(\omega_i\)相乘，这样就可以了：

scdid_df = (
    df_with_w
    .merge(time_w, on=["date"], how="left")
    .fillna({"time_weight": df_norm["post"].mean()})
    .assign(weight = lambda d: (d["time_weight"]*d["unit_weight"]))
)

现在你可以使用 \(scdid\_df\) 数据和加权回归来运行 \(DID\)。与 \(W_{it} = D_i * Post_t\) 相关联的参数估计值就是你所关注的 \(ATT\) 估计值：

did_model = smf.wls(
    "app_download_pct ~ treated*post",
    data=scdid_df.query("weight>1e-10"),
    weights=scdid_df.query("weight>1e-10")["weight"]).fit()

did_model.params["treated:post"]

## 0.004086196404101933

为了理解 \(SDID\) 的作用，你可以绘制处理单元的双重差分线，以及通过使用合成控制组的趋势并将其投射到处理单元的基线上而得到的反事实趋势（虚线）。两者之间的差异就是你刚刚估计的 \(ATT\)。我还在第二个图上绘制了时间权重。你可以看到合成双重差分法主要使用更接近干预后时期的时间段：

这个 \(SDID\) 估计值略高，但与 \(SC\) 的 \(ATT\) 没有太大差异。那么，为什么合成双重差分法很有趣呢？估计量中的 \(SC\) 部分使得 \(DID\) 的平行趋势假设更具合理性。如果你首先构建一个合成控制组来模拟处理单元，要得到平行趋势就容易得多。因此，与双重差分法和合成控制法相比，合成双重差分法往往具有更低的偏差。其次，合成双重差分法的方差也往往比这两种方法都低。如果感兴趣的话，原始论文中有模拟结果可以证明这一点。

---

原始的合成双重差分法

这个合成双重差分（SDID）估计量是对原始 SDID 估计量的简化，原始估计量由德米特里・阿尔汉格尔斯基等人在论文《Synthetic Difference in Differences》中提出。该论文针对单元权重提出了一个略有不同的优化目标：

\(\begin{aligned} \widehat{\omega}_i^{sdid} &= \underset{\omega}{\operatorname{argmin}} \left\| \overline{\mathbf{y}}_{pre,tr} - \left( \mathbf{Y}_{pre,co} \omega_{co} + \omega_0 \right) \right\|_2^2 + \zeta^2 T_{pre} \left\| \omega_{co} \right\|_2^2 \\ \text{s.t } &\sum \omega_i = 1 \text{ and } \omega_i > 0 \ \forall \ i \end{aligned}\)

首先，这个目标函数也允许存在截距项偏移。原因在于，合成控制组不需要完全匹配处理单元，只需要匹配其趋势即可，因为之后会将其代入 \(DID\) 模型中。其次，他们在单元权重上添加了一个 \(L2\) 惩罚项，其中包含了这个新的\(\zeta\)项：

\(\zeta = \left( N_{tr} * T_{post} \right)^{1/4} \sigma \left( Y_{it} - Y_{i(t-1)} \right)\)

关于这个\(\zeta\)，存在一个复杂的理论原因，我就不深入探讨了，但额外的 \(L2\) 惩罚背后的主要思想是确保没有单个单元获得不成比例的大权重。因此，与标准的合成控制方法相比，这些权重往往更分散地分布在更多的单元上。

这篇论文还提出了一种专门为 \(SDID\) 设计的新的推断程序。如果这还不足以成为你去研究它的理由，那我真不知道还有什么理由了。

---

当然，天下没有免费的午餐。通过允许截距项偏移，\(SDID\) 去掉了 \(SC\) 的凸性保障。根据具体情况，你可以认为这是好的，因为 \(SC\) 允许更多的灵活性；也可以认为这是不好的，因为它也允许危险的外推。

关键思想

如果说我希望你们从本章中掌握的一个要点，那就是可以采用基于模型的方法来估计 \(E[Y(0)_t|D = 1, Post = 1]\)：只需拟合一个模型，用一组完全相同的处理前时间序列来预测受处理对象的处理前结果，然后使用该模型对处理后时期进行预测。通常，这些时间序列是对照组单元的结果，而这种方法相当于在干预前时期对受处理结果与对照结果进行横向回归：

\(\widehat{\omega}^{sc} = \underset{\omega}{\text{argmin}} \ \Vert \overline{\mathbf{y}}_{pre,tr} - \mathbf{Y}_{pre,co}\omega_{co} \Vert_2^2\)

结果是一组权重，当这些权重与对照组单元相乘时，会产生一个合成控制组：这是对照组单元的一种组合，至少在干预前时期能近似模拟受处理单元的行为。如果这种近似效果良好，并且能推广到干预后时期，你就可以用它来估计 \(ATT\)：

\(\widehat{ATT} = \mathbf{Y}_{post,tr} - \mathbf{Y}_{post,co}\widehat{\omega}_{co}\)

这就是基本思想。你可以在此基础上进行拓展。例如，在典型的合成控制情境中，你会添加约束条件 \(\sum \omega_i = 1\) 且 \(\omega_i > 0\) 对所有 \(i\) 成立，或者使用类似套索回归的方法。但基本思想仍然是：利用干预前时期，将干预前的结果变量对预测性很强的时间序列进行回归，然后将这些预测延伸到干预后时期。

---

另见

\(徐轶青\)有很多关于推广合成控制法的论文，也有实现这些方法的软件。举几个例子，论文《广义合成控制方法：具有交互固定效应模型的因果推断》还融合了双重差分法和合成控制法的特点，将后者推广到了可变的处理时期（ staggered adoption design，交错采用设计）。在《比较案例研究中合成控制的贝叶斯替代方法》中，作者们也提出了一个贝叶斯模型来估计反事实情况。

---

---

实用示例

因果影响

谷歌的研究团队利用合成控制法背后的核心思想构建了 \(causalimpact\) 库。他们使用贝叶斯结构时间序列模型，通过其他自身未受处理影响的时间序列来估计 \(E[Y(0)|D = 1]\) 的反事实时间序列。该方法是贝叶斯方法这一事实，也让他们能够很自然地给出不确定性指标。

---