The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
103,048 characters
Synthetic Difference in Differences
\title{\textbf{Synthetic Difference in Differences}
\begingroup
\renewcommand{}\footnote{{\small We are grateful for
helpful comments and feedback from a co-editor and referees, as well as from Alberto Abadie, Avi Feller, Paul Goldsmith-Pinkham,
Liyang Sun, Yiqing Xu, Yinchu Zhu, and seminar participants at several venues. This research was generously supported
by ONR grant N00014-17-1-2131 and the Sloan Foundation.
The R package for implementing the methods developed here is available at https://github.com/synth-inference/synthdid.
The associated vignette is at https://synth-inference.github.io/synthdid/.}}
\endgroup
}
\author{Dmitry Arkhangelsky\thanks{{\small Associate Professor, CEMFI, Madrid, \texttt{[email removed]}.}}
\and
Susan Athey\thanks{{\small Professor of
Economics, Graduate School of Business, Stanford University, SIEPR, and NBER,
\texttt{[email removed]}. }}
\and David A.~Hirshberg\thanks{{\small Postdoctoral Fellow, Graduate School of Business and Department of Statistics, Stanford University, \texttt{[email removed]}.}}
\and Guido W.~Imbens\thanks{{\small Professor of Economics,
Graduate School of Business, and Department of Economics, Stanford University, SIEPR, and NBER,
\texttt{[email removed]}.}}
\and Stefan Wager\thanks{{\small Associate Professor of Operations, Information and Technology,
Graduate School of Business, and of Statistics (by courtesy), Stanford University, \texttt{[email removed]}.}} }
\date{Version \ifcase\month\or January\or February\or March\or April\or May\or June\or July\or August\or September\or October\or November\or December\fi \ \number \year\ \
}
\maketitle\thispagestyle{empty}
\linespread{1.5}
\begin{abstract}
\singlespacing
\noindent
We present a new estimator for causal effects with panel data that builds on insights behind the widely used difference in differences and synthetic control methods. Relative to these methods we find, both theoretically and empirically, that this ``synthetic difference in differences'' estimator has desirable robustness properties, and that it performs well in settings where the conventional estimators are commonly used in practice. We study the asymptotic behavior of the estimator when the systematic part of the outcome model includes latent unit factors interacted with latent time factors, and we present conditions for consistency and asymptotic normality.
\vspace{0.4\baselineskip}
\noindent \textbf{Keywords}: causal inference, difference in differences, synthetic controls, panel data
\end{abstract}
\baselineskip=20pt
\setcounter{page}{1}
\newpage
\section{Introduction}
Researchers are often interested in evaluating
the effects of policy changes using panel data, {\it i.e.}, using repeated observations of units across time, in a setting where
some units are exposed to the policy in some time periods but not others.
These policy changes are frequently not random---neither across units of analysis, nor across
time periods---and even unconfoundedness given observed covariates may not be credible ({\it e.g., } \citet{imbens2015causal}). In the absence of exogenous variation researchers have focused on statistical models that connect observed data to unobserved counterfactuals.
Many approaches have been developed for this setting but, in practice, a handful of methods
are dominant in empirical work. As documented by \citet*{currie2020technology}, Difference in Differences (DID)
methods have been widely used in applied economics over the last three decades; see also
\citet*{ashenfelter1985using}, \citet*{Bertrand2004did}, and \citet{angrist2008mostly}.
More recently, Synthetic Control (SC) methods, introduced in a
series of seminal papers by Abadie and coauthors \citep*{abadie2003,Abadie2010, abadie2014, abadie2016},
have emerged as an important alternative method for comparative case studies.
Currently these two strategies are often viewed as targeting different types of empirical applications.
In general, DID methods are applied in cases where we have a substantial number of units that are exposed
to the policy, and researchers are willing to make a ``parallel trends'' assumption which implies that we can
adequately control for selection effects by accounting for additive unit-specific and time-specific
fixed effects.
In contrast, SC methods, introduced in a setting with only a single (or small number)
of units exposed, seek to compensate for the lack of parallel trends
by re-weighting units to match their pre-exposure trends.
In this paper, we argue that although the empirical settings where DID and
SC methods
are typically used differ, the fundamental assumptions that justify both methods are closely related.
We then propose a new method, Synthetic Difference in Differences (SDID), that combines attractive features of both. Like
SC, our method re-weights and matches pre-exposure trends to weaken the reliance on parallel trend type
assumptions. Like DID, our method is invariant to additive unit-level shifts, and allows
for valid large-panel inference. Theoretically, we establish consistency and asymptotic normality of our
estimator. Empirically, we find that our method is competitive with (or dominates) DID
in applications where DID methods have been used in the past, and likewise is competitive with (or dominates)
SC in applications where SC methods have been used in the past.
To introduce the basic ideas,
consider a balanced panel with $N$ units and $T$ time periods, where the outcome for unit $i$ in period $t$ is denoted by $Y_{it}$,
and exposure to the binary treatment is denoted by $W_{it}\in\{0,1\}$. Suppose moreover that the first $N_{\rm co}$ (control) units are never
exposed to the treatment, while the last $N_{\rm tr} = N - N_{\rm co}$ (treated) units are exposed after time $T_{\rm pre}$.\footnote{Throughout the main part of our analysis, we focus on the block treatment assignment case where
$W_{it} = 1\left(\left\{i > N_{\rm co}, \,t > T_{\rm pre}\right\}\right)$. In the closely related staggered adoption case (\citet*{athey2021design}) where
units adopt the treatment at different times, but remain exposed after they first adopt the treatment, one can modify the methods
developed here.
See Section \ref{sec:staggered} in the Appendix for details.}
Like with
SC methods, we start by finding weights $\hat{\omega}^{\rm sdid}$ that align pre-exposure trends in the outcome of unexposed units with those
for the exposed units, {\it e.g.,}
$\sum_{i = 1}^{N_{\rm co}} \hat{\omega}^{\rm sdid}_i Y_{it} \approx N_{\rm tr}^{-1} \sum_{i = N_{\rm co} + 1}^{N} Y_{it}$ for all $t=1,\,\ldots,\,T_{\rm pre}$.
We also look for time weights $\hat{\lambda}^{\rm sdid}_t$ that balance pre-exposure time periods with
post-exposure ones (see Section \ref{sec:calif} for details).
Then we use these weights in a basic two-way fixed effects regression
to estimate the average causal effect of exposure (denoted by $\tau$):\footnote{This estimator
also has an interpretation as a difference-in-differences of weighted averages of observations.
See Equations~\ref{eq:tauhat_def}-~\ref{eq:delta_def} below.}
\begin{equation}
\label{main_sdid}
\left(\hat\tau^{\rm sdid}, \, \hat\mu, \, \hat\alpha, \, \hat\beta\right)=
\operatorname*{arg\,min}_{\tau,\mu,\alpha,\beta} \left\{ \sum_{i=1}^N\sum_{t=1}^T \Bigl( Y_{it}-\mu-\alpha_i-\beta_t-W_{it}\tau\Bigr)^2\hat{\omega}_i^{\rm sdid}\hat{\lambda}_t^{\rm sdid}\right\}.
\end{equation}
In comparison, DID estimates the effect of treatment exposure by solving the same two-way fixed effects regression problem without either time or unit weights:
\begin{equation}
\label{main_did}
\left(\hat\tau^{\rm did}, \, \hat\mu, \, \hat\alpha, \, \hat\beta\right)=
\operatorname*{arg\,min}_{\alpha,\beta,\mu,\tau} \left\{\sum_{i=1}^N\sum_{t=1}^T \Bigl( Y_{it}-\mu-\alpha_i-\beta_t-W_{it}\tau\Bigr)^2\right\}.
\end{equation}
The use of weights in the SDID estimator effectively makes the two-way fixed effect regression ``local,'' in that it emphasizes (puts more weight on) units that on average are similar
in terms of their past to the target (treated) units, and it emphasizes periods that are on average similar to the target (treated) periods.
This localization can bring two benefits relative to the standard DID estimator. Intuitively, using only similar units and similar periods makes the estimator more
robust. For example, if one is interested in estimating the effect of anti-smoking legislation on California (\citet*{Abadie2010}),
or the effect of German reunification on West Germany (\citet*{abadie2014}), or the effect of the Mariel boatlift on Miami (\citet*{cardmariel},
\citet*{peri2019labor}), it is natural to emphasize states, countries or cities that are similar to California, West Germany,
or Miami respectively relative to states, countries or cities that are not. Perhaps less intuitively, the use of the weights can also improve the estimator's precision by implicitly removing
systematic (predictable) parts of the outcome. However, the latter is not guaranteed: If there is little systematic heterogeneity in outcomes
by either units or time periods, the unequal weighting of units and time periods may worsen the precision of the estimators relative to the DID estimator.
Unit weights are designed so that the average outcome for the treated units is approximately parallel to the weighted average for control units. Time weights are designed so that
the average post-treatment outcome for each of the control units differs by a constant from the weighted average of the pre-treatment outcomes for the same control units.
Together, these weights make the DID strategy more plausible. This idea is not far from the current empirical practice. Raw data rarely exhibits parallel time trends for treated and control units, and researchers use different techniques, such as adjusting for covariates or selecting appropriate time periods to address this problem ({\it e.g.}, \citet{abadie2005semiparametric, callaway2020difference}). Graphical evidence that is used to support the parallel trends assumption is then based on the adjusted data. SDID makes this process automatic and applies a similar logic to weighting both units and time periods, all while retaining statistical guarantees. From this point of view, SDID addresses pretesting concerns recently expressed in \citet{roth2018pre}.
In comparison with the SDID estimator, the SC estimator omits the unit fixed effect and the time weights from the regression function:
\begin{equation}
\label{main_sc}
\left(\hat\tau^{\rm sc}, \, \hat\mu, \, \hat\beta\right)=
\operatorname*{arg\,min}_{\mu,\beta,\tau} \left\{\sum_{i=1}^N\sum_{t=1}^T \Bigl( Y_{it}-\mu-\beta_t-W_{it}\tau\Bigr)^2\hat{\omega}^{\rm sc}_i\right\}.
\end{equation}
The argument for including time weights in the SDID estimator is the same as the argument for including the unit weights presented earlier:
The time weight can both remove bias and
improve precision by eliminating the role of time periods that are very different from the post-treatment periods.
Similar to the argument for the use of weights, the argument for the inclusion of the unit fixed effects is twofold. First, by making
the model more flexible, we strengthen its robustness properties. Second, as demonstrated in the application
and simulations based on real data,
these unit fixed effects often explain much of the variation in
outcomes and can improve precision. Under some conditions, SC weighting can account for the
unit fixed effects on its own. In particular, this happens when the weighted average of the outcomes for the control
units in the pre-treatment periods is exactly equal to the average of outcomes for the treated
units during those pre-treatment periods. In practice, this equality holds only approximately,
in which case including the unit fixed effects in the weighted regression will remove some of the
remaining bias. The benefits of including unit fixed effects in the SC regression \eqref{main_sc} can also be obtained by applying the synthetic control method after centering the data by subtracting, from each unit's trajectory, its pre-treatment mean.
This estimator was previously suggested
in \citet{doudchenko2016balancing} and \citet{ferman2019synthetic}.
To separate out the benefits of allowing for fixed effects from those stemming from the use of time-weights, we include in our application and simulations this DIFP estimator.
\section{An Application}
\label{sec:calif}
To get a better understanding of how $\hat{\tau}^{\rm did}$, $\hat{\tau}^{\rm sc}$ and $\hat{\tau}^{\rm sdid}$ compare to each other,
we first revisit the California smoking cessation program example of \citet*{Abadie2010}.
The goal of their analysis was to estimate the effect of increased cigarette taxes on smoking in California.
We consider observations for 39 states (including California) from 1970 through 2000. California
passed Proposition 99 increasing cigarette taxes (i.e., is treated) from 1989 onwards.
Thus, we have $T_{\rm pre}=19$ pre-treatment periods, $T_{\rm post}=T-T_{\rm pre}=12$ post-treatment periods, $N_{\rm co}=38$ unexposed states,
and $N_{\rm tr}=1$ exposed state (California).
\subsection{Implementing SDID} \label{sec:implementing}
Before presenting results on the California smoking case, we discuss in detail how we choose the synthetic control
type weights $\hat{\omega}^{\rm sdid}$ and $\hat{\lambda}^{\rm sdid}$ used for our estimator as specified in \eqref{main_sdid}. Recall that,
at a high level, we want to choose the unit weights to roughly match pre-treatment trends of unexposed units with those for the exposed ones,
$\sum_{i = 1}^{N_{\rm co}} \hat{\omega}^{\rm sdid}_i Y_{it} \approx N_{\rm tr}^{-1} \sum_{i = N_{\rm co} + 1}^{N} Y_{it}$ for all $t=1,\,\ldots,\,T_{\rm pre}$,
and similarly we want to choose the time weights to balance pre- and post-exposure periods for unexposed units.
In the case of the unit weights $\hat{\omega}^{\rm sdid}$, we implement this by solving the optimization problem
\begin{equation}
\label{unit_weights}
\begin{aligned}
&\left(\hat{\omega}_0, \, \hat{\omega}^{\rm sdid}\right) = \operatorname*{arg\,min}_{ \omega_0 \in \mathbb{R}, \omega\in \Omega} \ell_{unit}(\omega_0, \omega) \quad \text{ where } \\
& \quad \quad \quad\ell_{unit}(\omega_0, \omega) = \sum_{t=1}^{T_{\rm pre}}
\left( \omega_0 + \sum_{i=1}^{N_{\rm co}} \omega_{i} Y_{it} -
\frac{1}{N_{\rm tr}} \sum_{i = N_{\rm co} + 1}^N Y_{it}\right)^2 + \zeta^2 T_{{\rm pre}} \left\|\omega\right\|_2^2, \\
& \quad \quad \quad \Omega = \left\{\omega \in \mathbb{R}_+^N : \sum_{i=1}^{N_{\rm co}} \omega_i=1, \, \omega_i=N_{\rm tr}^{-1} \text{ for all } i=N_{\rm co}+1, \, \ldots \, ,N\right\},
\end{aligned}
\end{equation}
where $\mathbb{R}_+$ denotes the positive real line. We set the regularization parameter $\zeta$ as
\begin{equation}
\label{eq:zeta_calif}
\begin{split}
&\zeta= (N_{\rm tr} T_{\rm post})^{1/4} \ \hat \sigma \ \text{ with }\ \hat \sigma^2 = \frac{1}{N_{\rm co} (T_{\rm pre}-1)}\sum_{i=1}^{N_{\rm co}}\sum_{t=1}^{T_{\rm pre}-1} \left(\Delta_{it}-\overline{\Delta}\right)^2, \\
&\mathrm{where}\ \ \Delta_{it}=Y_{i(t+1)}-Y_{it},\hskip0.5cm {\rm and}\ \ \overline{\Delta}=\frac{1}{N_{\rm co} (T_{\rm pre}-1)}\sum_{i=1}^{N_{\rm co}}\sum_{t=1}^{T_{\rm pre}-1} \Delta_{it}.
\end{split}
\end{equation}
That is, we choose the regularization parameter $\zeta$ to match the size of a typical one-period outcome change $\Delta_{it}$ for unexposed units in the pre-period, multiplied by a theoretically motivated scaling $(N_{\rm tr} T_{\rm post})^{1/4}$.
The SDID weights $\hat{\omega}^{\rm sdid}$ are closely related to the weights used in \citet*{Abadie2010}, with two minor differences.
First, we allow for an intercept term $\omega_0$, meaning that the weights $\hat{\omega}^{\rm sdid}$ no longer need to make
the unexposed pre-trends perfectly match the exposed ones; rather, it is sufficient that the weights make the trends parallel. The reason we
can allow for this extra flexibility in the choice of weights is that our use of fixed effects $\alpha_i$ will absorb any constant
differences between different units. Second, following \citet{doudchenko2016balancing}, we add a regularization penalty
to increase the dispersion, and ensure the uniqueness, of the weights. If we were to omit the intercept
$\omega_0$ and set $\zeta = 0$, then \eqref{unit_weights} would correspond exactly to a choice of weights discussed
in \citet{Abadie2010} in the case where $N_{\rm tr} = 1$.
We implement this for the time weights $\hat{\lambda}^{\rm sdid}$ by solving\footnote{The weights $\hat{\lambda}^{\rm sdid}$ may not be uniquely defined, as $\ell_{time}$ can have multiple minima. In principle our results hold for any argmin of $\ell_{time}$. These tend to be similar in the setting we consider, as they all converge to unique `oracle weights' $\tilde \lambda^{\rm sdid}$ that are discussed in Section~\ref{sec:oracle-and-adaptive}. In practice, to make the minimum defining our time weights unique, we add a very small regularization term $\zeta^2 N_{{\rm co}} \lVert\lambda\rVert^2$
to $\ell_{time}$, taking $\zeta = 10^{-6} \ \hat \sigma$ for $\hat \sigma$ as in \eqref{eq:zeta_calif}.
}
\begin{equation}
\label{time_weights}
\begin{aligned}
&\left(\hat{\lambda}_0, \, \hat{\lambda}^{\rm sdid}\right) = \operatorname*{arg\,min}_{ \lambda_0 \in \mathbb{R}, \lambda \in \Lambda} \ell_{time}(\lambda_0, \lambda) \quad \text{ where } \\
&\quad \quad \quad \ell_{time}(\lambda_0, \lambda) = \sum_{i=1}^{N_{\rm co}}
\left( \lambda_0 + \sum_{t=1}^{T_{\rm pre}} \lambda_{t} Y_{it} -
\frac{1}{T_{\rm post}} \sum_{t = T_{\rm pre} + 1}^T Y_{it}\right)^2 , \\
&\quad \quad \quad \Lambda=\left\{\lambda\in\mathbb{R}_+^T : \sum_{t=1}^{T_{\rm pre}} \lambda_t=1, \, \lambda_t=T_{\rm post}^{-1} \text{ for all } t=T_{\rm pre}+1,\, \ldots \, ,T\right\}.
\end{aligned}
\end{equation}
The main difference between \eqref{unit_weights} and \eqref{time_weights} is that we use regularization for the
former but not the latter. This choice is motivated by our formal
results, and reflects the fact we allow for correlated observations within time periods for the same unit, but not across units within
a time period, beyond what is captured by the systematic component of outcomes as represented by a latent factor model.
\RestyleAlgo{boxruled}
\LinesNumbered
\begin{algorithm}[t]
\KwData{$\boldsymbol{Y}, \boldsymbol{W}$}
\KwResult{Point estimate $\hat \tau^{\rm sdid}$}
Compute regularization parameter $\zeta$ using \eqref{eq:zeta_calif}\;
Compute unit weights $\hat{\omega}^{\rm sdid}$ via \eqref{unit_weights}\;
Compute time weights $\hat{\lambda}^{\rm sdid}$ via \eqref{time_weights}\;
Compute the SDID estimator via the weighted DID regression
\begin{equation*}
\label{main_sdid_cov}
\left(\hat\tau^{\rm sdid}, \, \hat\mu, \, \hat\alpha, \, \hat\beta\right)=
\operatorname*{arg\,min}_{\tau,\mu,\alpha,\beta} \left\{ \sum_{i=1}^N\sum_{t=1}^T \Bigl( Y_{it}-\mu-\alpha_i-\beta_t-W_{it}\tau\Bigr)^2\hat{\omega}_i^{\rm sdid}\hat{\lambda}_t^{\rm sdid}\right\};
\end{equation*}
\caption{Synthetic Difference in Differences (SDID)}
\label{alg:sdid}
\end{algorithm}
We summarize our procedure as Algorithm \ref{alg:sdid}.\footnote{Some applications feature time-varying
exogenous covariates $X_{it} \in \mathbb{R}^p$.
We can incorporate adjustment for these covariates by applying SDID to the residuals $Y_{it}^{\mathrm{res}} = Y_{it}-X_{it}\hat\beta$ of the regression of $Y_{it}$ on $X_{it}$.}
In our application and simulations we also report
the SC and DIFP estimators. Both of these use
weights solving \eqref{unit_weights} without regularization. The SC estimator also omits the intercept $\omega_0$.\footnote{Like the time weights $\hat\lambda^{{\rm sdid}}$, the unit weights
for the SC and DIFP estimators may not be uniquely defined. To ensure uniqueness in practice, we take $\zeta = 10^{-6} \ \hat \sigma$, not $\zeta=0$, in $\ell_{unit}$.
In our simulations, SC and DIFP with this minimal form of regularization outperform more strongly regularized variants with $\zeta$ as in \eqref{eq:zeta_calif}. We show this comparison in Table \ref{table_reg_vs_unreg}.}
Finally, we report results for the matrix completion (MC) estimator proposed by \citet{athey2017matrix}, which is based on imputing the missing $Y_{it}(0)$ using a low rank factor model with nuclear norm regularization.
\subsection{The California Smoking Cessation Program}
The results from running this analysis
are shown in Table \ref{estimates_smoking}. As argued in \citet{Abadie2010}, the assumptions underlying
the DID estimator are suspect here, and the -27.3 point estimate likely overstates the
effect of the policy change on smoking. SC provides a reduced (and generally considered more credible) estimate
of -19.6. The other methods, our proposed SDID, the DIFP and the MC estimator are all smaller than the DID estimator with the SDID and DIFP estimator substantially smaller than the SC estimator. At the very least, this difference in point estimates implies that the
use of time weights and unit fixed effects in \eqref{main_sdid} materially affects conclusions; and, throughout this
paper, we will argue that when $\hat{\tau}^{\rm sc}$ and $\hat{\tau}^{\rm sdid}$ differ, the latter is often more credible.
Next, and perhaps surprisingly, we see that the standard errors obtained for SDID
(and also for SCIFP, and MC) are smaller than those for DID, despite our method
being more flexible. This is a result of the local fit of SDID (and SC) being improved by the weighting.
\begin{table}[t]
\begin{center}
\begin{tabular}{|l|ccccc|}
\hline
& SDID & SC & DID & MC & DIFP \\
\hline
Estimate &
-15.6 & -19.6 &-27.3 & -20.2 & -11.1 \\
Standard error &
(8.4) & (9.9)& (17.7) & (11.5) & (9.5)
\\
\hline
\end{tabular}
\caption{Estimates for average effect of increased cigarette taxes on California per capita cigarette sales over twelve post-treatment years,
based on
synthetic difference in differences (SDID),
synthetic controls (SC),
difference in differences (DID),
matrix completion (MC), synthetic control with intercept (DIFP),
along with estimated standard errors. We use the `placebo method' standard error estimator discussed in Section~\ref{sec:inference}.}
\label{estimates_smoking}
\end{center}
\end{table}
To facilitate direct comparisons, we observe that each of the three estimators can be rewritten as a weighted average difference in adjusted outcomes $\hat{\delta}_i$ for appropriate sample weights $\hat{\omega}_i$:
\begin{equation}
\label{eq:tauhat_def}
\hat{\tau} = \hat{\delta}_{\rm tr} - \sum_{i = 1}^{N_{\rm co}} \hat{\omega}_i \hat{\delta}_i
\quad \text{ where }\quad \hat{\delta}_{\rm tr} = \frac{1}{N_{{\rm tr}}}\sum_{i=N_{\rm co}+1}^N \hat{\delta}_i.
\end{equation}
DID uses constant weights $\hat{\omega}_i^{\rm did} = N_{\rm co}^{-1}$, while the construction of SDID and SC weights is outlined in Section \ref{sec:implementing}.
For the adjusted outcomes $\hat{\delta}_i$, SC uses unweighted treatment period averages,
DID uses unweighted differences between average treatment period and pre-treatment outcomes,
and SDID uses weighted differences of the same.
\begin{equation}
\label{eq:delta_def}
\begin{split}
&\hat{\delta}_i^{\rm sc} = \frac{1}{T_{\rm post}} \sum_{t = T_{\rm pre} + 1}^T Y_{it}, \\
&\hat{\delta}_i^{\rm did} = \frac{1}{T_{\rm post}} \sum_{t = T_{\rm pre} + 1}^T Y_{it} - \frac{1}{T_{\rm pre}} \sum_{t = 1}^{T_{\rm pre}} Y_{it}, \\
&\hat{\delta}_i^{\rm sdid} = \frac{1}{T_{\rm post}} \sum_{t = T_{\rm pre} + 1}^T Y_{it} - \sum_{t = 1}^{T_{\rm pre}} \hat{\lambda}_t^{\rm sdid} Y_{it}.
\end{split}
\end{equation}
The top panel of Figure \ref{californiadiagrams} illustrates how each method operates.
As is well known \citep{ashenfelter1985using}, DID relies on the assumption that cigarette sales in different states would have evolved
in a parallel way absent the intervention. Here, pre-intervention trends are obviously not parallel, so the DID
estimate should be considered suspect. In contrast, SC re-weights
the unexposed states so that the weighted of outcomes for these states match California pre-intervention as close as possible, and then attributes
any post-intervention divergence of California from this weighted average to the intervention.
What SDID does here is to re-weight the unexposed control units to make their time trend
parallel (but not necessarily identical) to California pre-intervention, and then applies a DID
analysis to this re-weighted panel. Moreover, because of the time weights, we only
focus on a subset of the pre-intervention time periods when carrying out this last step. These time periods were
selected so that the weighted average of historical outcomes predict average treatment period outcomes for control units, up to a constant.
It is useful to contrast the data-driven SDID approach to selecting the time weights to both DID, where all pre-treatment periods are given equal weight, and to event studies where typically the last pre-treatment period is used as a comparison and so implicitly gets all the weight ({\it e.g.,} \citet{borusyak2016revisiting, freyaldenhoven2019pre}).
\begin{figure}[p]
\begin{center}
\begin{tabular}{cccc} & \
\makebox[0.3\textwidth]{Difference in Differences} &
\makebox[0.3\textwidth]{Synthetic Control} &
\makebox[0.3\textwidth]{Synthetic Diff.~in Differences}
\end{tabular}
\begin{minipage}{.02\textwidth}
\begin{tabular}{c}
\rotatebox{90}{\tiny \qquad\qquad cigarette consumption (packs/year)} \\
\rotatebox{90}{\tiny \quad difference in consumption (packs/year) \quad}
\end{tabular}
\end{minipage}
\begin{minipage}{.97\textwidth}
\includegraphics[width=\textwidth]{california-trends-plot-1.pdf}
\includegraphics[width=\textwidth]{california-units-plot-1.pdf}
\end{minipage}
\caption{A comparison between difference-in-differences, synthetic control, and synthetic differences-in-differences estimates
for the effect of California Proposition 99 on per-capita annual cigarette consumption (in packs/year).
In the first row, we show trends in consumption over time for California and the relevant weighted average of control states,
with the weights used to average pre-treatment time periods at the bottom of the graphs.
The estimated effect is indicated by an arrow.
In the second row, we show the state-by-state adjusted outcome difference $\hat{\delta}_{\rm tr} -\hat{\delta}_i$ as specified in \eqref{eq:tauhat_def}-\eqref{eq:delta_def}, with the weights $\hat{\omega}_i$ indicated by dot size and the weighted average of these differences --- the estimated effect --- indicated by a horizontal line. Observations with zero weight are denoted by an $\times$-symbol.
\label{californiadiagrams}}
\end{center}
\end{figure}
The lower panel of Figure \ref{californiadiagrams} plots $\hat{\delta}_{\rm tr} - \hat{\delta}_i$ for each method and for each unexposed state,
where the size of each point corresponds to its weight $\hat{\omega}_i$; observations with zero weight are denoted by
an $\times$-symbol. As discussed in \citet*{Abadie2010}, the
SC weights $\hat{\omega}^{\rm sc}$ are sparse. The SDID weights $\hat{\omega}^{\rm sdid}$ are also sparse---but
less so. This is due to regularization and the use of the intercept $\omega_0$, which allows greater flexibility in solving \eqref{unit_weights}, enabling more balanced weighting. Observe that both DID
and SC have some very high influence states, that is, states with large absolute values of $\hat{\omega}_i (\hat{\delta}_{\rm tr} - \hat{\delta}_i)$
({\it e.g.}, in both cases, New Hampshire). In contrast, SDID does not give any state particularly high influence,
suggesting that after weighting, we have achieved the desired ``parallel trends'' as illustrated in the top panel of Figure \ref{californiadiagrams} without inducing excessive variance in the estimator by using concentrated weights.
\section{Placebo Studies}
\label{sec:simu}
So far, we have relied on conceptual arguments to make the claim that SDID
inherits good robustness properties from both traditional DID and SC methods,
and shows promise as a method that can be used is settings where either DID and SC
would traditionally be used. The goal of this section is to see how these claims play out in realistic empirical settings. To this end, we consider
two carefully crafted simulation studies, calibrated to datasets representative of those typically used for panel data studies. The first simulation study mimics settings where DID would be used in practice (Section \ref{sec:cps}), while the second mimics settings suited to SC (Section \ref{sec:penn}). Not only do we base the outcome model of our simulation study on real datasets, we further ensure
that the treatment assignment process is realistic by seeking to emulate the distribution of real policy initiatives. To be specific,
in Section \ref{sec:cps}, we consider a panel of US states. We estimate several alternative treatment assignment models to create the hypothetical treatments, where the models are based on the state laws related to minimum wages,
abortion or gun rights.
In order to run such a simulation study, we first need to commit to an econometric specification
that can be used to assess the accuracy of each method. Here, we work with the following latent factor model
(also referred to as an ``interactive fixed-effects model'', \citet{xu2017generalized}, see also \citet{athey2017matrix}),
\begin{equation}
\label{basic_model_scalar}
Y_{it}=\boldsymbol{\gamma}_i \boldsymbol{\upsilon}_t^\top +\tau W_{it} +\varepsilon_{it},
\end{equation}
where $\boldsymbol{\gamma}_i$ is a vector of latent unit factors of dimension $R$, and $\boldsymbol{\upsilon}_t$ is a vector of latent time factors of dimension $R$. In matrix form, this can be written
\begin{equation}
\label{basic_model}
\boldsymbol{Y}=\boldsymbol{L}+\tau\boldsymbol{W}+\boldsymbol{E}
\quad \text{ where } \quad
\boldsymbol{L}=\boldsymbol{\Gamma \Upsilon}^\top.
\end{equation}
We refer to $\boldsymbol{E}$ as the idiosyncratic component or error matrix, and to $\boldsymbol{L}$ as the systematic\ component.
We assume that
the conditional expectation of the error matrix $\boldsymbol{E}$ given the assignment matrix $\boldsymbol{W}$ and the systematic\ component $\boldsymbol{L}$ is zero. That is,
the treatment assignment cannot depend on $\boldsymbol{E}$. However, the treatment assignment may in general depend on the
systematic component $\boldsymbol{L}$ ({\it i.e.}, we do not take $\boldsymbol{W}$ to be randomized). We assume that $\boldsymbol{E}_i$ is independent of $\boldsymbol{E}_{i'}$ for each pair of units $i,i'$, but we allow for correlation across time periods within a unit. Our goal is to estimate the treatment effect $\tau$.
The model \eqref{basic_model} captures several qualitative challenges that have received considerable attention
in the recent panel data literature. When the matrix $\boldsymbol{L}$ takes on an additive form, i.e., $L_{it} = \alpha_i + \beta_t$,
then the DID regression will consistently recover $\tau$. Allowing for interactions
in $\boldsymbol{L}$ is a natural way to generalize the fixed-effects specification and discuss inference in settings where
DID is misspecified \citep{bai2009panel,moon2015linear,moon2017dynamic}. In our formal results given in
Section \ref{sec:formal}, we show how, despite not explicitly fitting the model \eqref{basic_model}, SDID can consistently estimate $\tau$ in this design under reasonable conditions. Finally, accounting for correlation over time within observations of the same unit is widely considered to be an important ingredient to credible inference
using panel data \citep*{angrist2008mostly, Bertrand2004did}.
In our experiments, we compare DID, SC, SDID, and DIFP, all implemented exactly as in Section \ref{sec:calif}.
We also compare these four estimators to an alternative that
estimates $\tau$ by directly fitting both $\boldsymbol{L}$ and $\tau$ in \eqref{basic_model}; specifically, we consider
the matrix completion (MC) estimator recommended in \citet*{athey2017matrix} which uses nuclear norm
penalization to regularize its estimate of $\boldsymbol{L}$.
In the remainder of this section, we focus on comparing the bias and root-mean-squared error of the estimator.
We discuss questions around inference and coverage in Section \ref{sec:inference}.
\subsection{Current Population Survey Placebo Study}
\label{sec:cps}
Our first set of simulation experiments revisits the landmark placebo study of \citet*{Bertrand2004did} using
the Current Population Survey (CPS). The main goal of \citet{Bertrand2004did} was to study the behavior
of different standard error estimators for DID. To do so, they randomly assigned a subset of states in the CPS
dataset to a placebo treatment and the rest to the control group, and examined how well different approaches
to inference for DID estimators covered the true treatment effect of zero. Their main finding was that only methods that were robust to serial correlation of
repeated observations for a given unit ({\it e.g.,} methods that clustered observations by unit) attained valid coverage.
We modify the placebo analyses in \citet{Bertrand2004did} in two ways. First, we no longer assigned exposed states
completely at random, and instead use a non-uniform assignment mechanism that is inspired by different policy
choices actually made by different states. Using a non-uniformly random assignment is important because it allows us to differentiate between various estimators in ways that completely random assignment would not. Under completely random assignment, a number of methods, including DID, perform well because the presence of $\boldsymbol{L}$ in
\eqref{basic_model} introduces zero bias. In contrast, with a non-uniform random assignment ({\it i.e.},
treatment assignment is correlated with systematic effects), methods that
do not account for the presence of $\boldsymbol{L}$ will be biased.
Second, we simulate values for the outcomes based on a model estimated on the CPS data, in order to have more control over the data generating process.
\subsubsection{The Data Generating Process}
\label{sec:cps_dgp}
For the first set of simulations we use as the starting point data on wages for women with positive wages in the March outgoing rotation
groups in the Current Population Survey (CPS) for the years 1979 to 2019. We first transform these by taking logarithms and then average them by state/year cells.
Our simulation design has two components, an outcome model and an assignment model. We generate outcomes
via a simulation that seeks to capture the behavior of the average by state/year of the logarithm of wages for those
with positive hours worked in the CPS data as in \citet{Bertrand2004did}. Specifically, we simulate data using the
model \eqref{basic_model}, where the rows $\boldsymbol{E}_{i}$ of $\boldsymbol{E}$
have a multivariate Gaussian distribution $\boldsymbol{E}_{i} \sim \mathcal{N}(0, \Sigma)$, and
we choose both $\boldsymbol{L}$ and $\Sigma$ to fit the CPS data as follows. First, we fit a rank four factor model for $\boldsymbol{L}$:
\begin{equation}
\boldsymbol{L} :=\operatorname*{arg\,min}_{L: \text{rank}(L) = 4} \sum_{it}(Y^*_{it} - L_{it})^2,
\end{equation}
where $Y^*_{it}$ denotes the true state/year average of log-wage in the CPS data. We then estimate
$\Sigma$ by fitting an AR(2) model to the residuals of $Y^*_{it} - L_{it}$. For purpose of interpretation,
we further decompose the systematic component $\boldsymbol{L}$ into an additive (fixed effects) term $\boldsymbol{F}$ and
an interactive term $\boldsymbol{M}$, with
\begin{equation}
\begin{aligned}
&F_{it} = \alpha_i+\beta_t=\frac{1}{T}\sum_{l=1}^T L_{il} + \frac{1}{N}\sum_{j=1}^N L_{jt} - \frac{1}{NT} \sum_{it}L_{it},\\
& M_{it}=L_{it}-F_{it}.
\end{aligned}
\end{equation}
This decomposition of $\boldsymbol{L}$ into an additive two-way fixed effect component $\boldsymbol{F}$ and an interactive component $\boldsymbol{M}$ enables us to study the sensitivity
of different estimators to the presence of different types of systematic effects.
Next we discuss generation of the treatment assignment. Here, we are designing a ``null effect'' study, meaning that
treatment has no effect on the outcomes and all methods should estimate zero. However, to make this more challenging, we choose the treated units so that the assignment mechanism is correlated with the systematic component $\boldsymbol{L}$.
We set $W_{it} = D_i \boldsymbol{1}_{t > T_0}$, where $D_i$ is a binary exposure indicator generated as
\begin{equation}
\label{eq:assignment_lr}
D_i \,\big|\, \boldsymbol{E}_{i}, \alpha_i, \boldsymbol{M}_i \sim \text{Bernoulli}\left(\pi_i\right), \ \ \ \
\pi_i = \pi(\alpha_i, \boldsymbol{M}_i;\phi) = \frac{\exp(\phi_\alpha \alpha_i+\phi_M \boldsymbol{M}_i)}
{1+\exp(\phi_\alpha \alpha_i+\phi_M \boldsymbol{M}_i)}.
\end{equation}
In particular, the distribution of $D_i$ may depend on $\alpha_i$ and $\boldsymbol{M}_i$; however, $D_i$ is independent of $\boldsymbol{E}_i$, i.e., the assignment is strictly exogenous.\footnote{In the simulations below, we restrict the maximal number of treated units (either to $10$ or $1$). To achieve this, we first sample $D_i$ independently and accept the results if the number of treated units satisfies the constraint. If it does not, then we choose the maximal allowed number of treated units from those selected in the first step uniformly at random.}
To construct probabilities $\{\pi_i\}$ for this assignment model, we choose $\phi$ as the coefficient estimates from a logistic regression of an observed binary characteristic of the state $D_i$ on $\boldsymbol{M}_i$ and $\alpha_i$. We consider three different choices for $D_i$, relating to minimum wage laws, abortion rights, and gun control laws.\footnote{See Section~\ref{sec:placebo-study-details} in the appendix for details.} As a result, we get assignment probability models that reflect actual differences across states with respect to important economic variables. In practice the $\alpha_i$ and $\boldsymbol{M}_i$ that we construct predict a sizable part of variation in $D_i$, with $R^2$ varying from $15\%$ to $30\%$.
\begin{table}[t]
\begin{center}
\begin{adjustbox}{width=1\textwidth}
\begin{tabular}{|l|rrrr|rrrrr|rrrrr|}
\hline
& \multirow{2}{*}{$\frac{\|\boldsymbol{F}\|_F}{\sqrt{NT}}$} & \multirow{2}{*}{$\frac{\|\boldsymbol{M}\|_F}{\sqrt{NT}}$} & \multirow{2}{*}{$\sqrt{\frac{\operatorname{tr}(\Sigma)}{T}}$} & \multirow{2}{*}{AR(2)} &
\multicolumn{5}{c|}{RMSE}& \multicolumn{5}{c|}{Bias}\\
& && & &SDID & SC & DID &MC & DIFP & SDID & SC & DID &MC& DIFP \\ \hline
&&&& &&&& &&&& &&\\
Baseline & 0.992 & 0.100 & 0.098 & (.01,-.06) & \bf 0.028 & 0.037 & 0.049 & 0.035 & 0.032 & 0.010 & 0.020 & 0.021 & 0.015 & 0.007 \\ \hline
\multicolumn{15}{|l|}{\it Outcome Model}\\
No Corr& 0.992 & 0.100 & 0.098 & (.00, .00) & \bf 0.028 & 0.038 & 0.049 & 0.035 & 0.032 & 0.010 & 0.020 & 0.021 & 0.015 & 0.007 \\
No $\boldsymbol{M}$ & 0.992 & 0.000 & 0.098 & (.01, -.06) & 0.016 & 0.018 & \bf 0.014 & \bf 0.014 & 0.016 & 0.001 & 0.004 & 0.001 & 0.001 & 0.001 \\
No $\boldsymbol{F}$ & 0.000 & 0.100 & 0.098 & (.01, -.06) & 0.028 & \bf 0.023 & 0.049 & 0.035 & 0.032 & 0.010 & 0.004 & 0.021 & 0.015 & 0.007\\
Only Noise& 0.000 & 0.000 & 0.098 & (.01, -.06) & 0.016 & \bf 0.014 & \bf 0.014 & \bf 0.014 & 0.016 & 0.001 & 0.001 & 0.001 & 0.001 & 0.001 \\
No Noise & 0.992 & 0.100 & 0.000 & (.00, .00) & 0.006 & 0.017 & 0.047 & \bf 0.004 & 0.011 & 0.004 & 0.004 & 0.020 & 0.000 & 0.001 \\ \hline
\multicolumn{15}{|l|}{\it Assignment Process}\\
Gun Law &0.992 & 0.100 & 0.098 & (.01, -.06) &\bf 0.026 & 0.027 & 0.047 & 0.035 & 0.030 & 0.008 & -0.003 & 0.015 & 0.015 & 0.009 \\
Abortion & 0.992 & 0.100 & 0.098 & (.01, -.06) & \bf 0.023 & 0.031 & 0.045 & 0.031 & 0.027 & 0.004 & 0.016 & 0.003 & 0.003 & 0.001
\\ Random & 0.992 & 0.100 & 0.098 & (.01, -.06)& \bf 0.024 & 0.025 & 0.044 & 0.031 & 0.027 & 0.001 & -0.001 & 0.002 & 0.001 & -0.000 \\ \hline
\multicolumn{15}{|l|}{\it Outcome Variable}\\
Hours & 0.789 & 0.402 & 0.575 & (.06, .00) & 0.190 & 0.203 & 0.206 & \bf 0.185 & 0.197 & 0.111 & -0.049 & 0.085 & 0.100 & 0.099 \\
U-rate & 0.752 & 0.441 & 0.593 & (-.02, -.01) & 0.191 & \bf 0.184 & 0.353 & 0.247 & 0.187 & 0.100 & 0.080 & 0.304 & 0.187 & 0.078\\
\hline
\multicolumn{15}{|l|}{\it Assignment Block Size}\\
$T_{\rm post} =1$ & 0.992 & 0.100 & 0.098 & (.01, -.06) & \bf 0.050 & 0.059 & 0.070 & 0.051 & 0.054 & 0.019 & 0.017 & 0.038 & 0.021 & 0.012\\
$N_{\rm tr} = 1$ & 0.992 & 0.100 & 0.098 & (.01, -.06) & \bf 0.063 & 0.072 & 0.126 & 0.081 & 0.083 & 0.002 & 0.014 & 0.011 & 0.004 & -0.002 \\
$ T_{\rm post} = N_{\rm tr} = 1$ & 0.992 & 0.100 & 0.098 & (.01, -.06) & 0.112 & 0.124 & 0.153 & \bf 0.108 & 0.117 & 0.014 & 0.024 & 0.033 & 0.016 & 0.011 \\ \hline
\end{tabular}
\end{adjustbox}
\caption{Simulation Results for CPS Data.
The baseline case uses state minimum wage laws to simulate treatment assignment, and generates outcomes using the full
data-generating process described in Section \ref{sec:cps_dgp}, with $T_{\rm post}=10$ post-treatment periods and at most $N_{\rm tr}=10$
treatment states. In subsequent settings, we omit parts of the data-generating process (rows 2-6), consider different distributions
for the treatment exposure variable $D_i$ (rows 7-9), different distributions for the outcome variable (rows 10-11), and vary the number of treated cells (rows 12-14).
The full dataset has $N=50$, $T=40$, and outcomes are normalized to have mean zero and unit variance.
All results are based on 1000 simulation replications.
\label{table1}}
\end{center}
\end{table}
\subsubsection{Simulation Results}
Table \ref{table1} compares the performance of the four aforementioned estimators in the simulation design described above. We
consider various choices for the number of treated units and the treatment assignment distribution. Furthermore,
we also consider settings where we drop various components of the outcome-generating process, such as the fixed effects $\boldsymbol{F}$ or
the interactive component $\boldsymbol{M}$, or set the noise correlation matrix $\Sigma$ to be diagonal.
The magnitude of the $\boldsymbol{F}$, $\boldsymbol{M}$ and $\boldsymbol{E}$ components as well as the strength of the
autocorrelation effects in $\Sigma$ captured by the first two autoregressive coefficients are shown in the first four columns of Table \ref{table1}.
At a high level, we find that SDID has excellent performance relative to the benchmarks ---both in terms of bias
and root-mean squared error. This holds in the baseline simulation design and over a number of other designs where we vary the treatment assignment (from being based on minimum wage laws to gun laws, abortion laws, or completely random), the outcome (from average of log wages to average hours and unemployment rate), and the maximal number of treated units (from 10 to 1) and the number of exposed periods (from 10 to 1). We find that when the treatment assignment is uniformly random, all methods are essentially
unbiased, but SDID is more precise. Meanwhile, when the treatment assignment is not uniformly random, SDID
is particularly successful at mitigating bias while keeping variance in check.
In the second panel of
Table \ref{table1}
we provide some additional insights into the superior performance of the SDID estimator by
sequentially dropping some of the components of the model that generates the potential outcomes.
If we drop the interactive component $\boldsymbol{M}$ from the outcome model (``No $\boldsymbol{M}$''), so that the fixed effect specification is correct, the DID estimator performs best (alongside MC).
In contrast, if we drop the fixed effects component $\bf$ (``No $\boldsymbol{F}$'') but keep the interactive component, the SC estimator does best.
If we drop both parts of the systematic component, and there is only noise,
the superiority of the SDID estimator vanishes and all estimators are essentially equivalent.
On the other hand, if we remove the noise component so that there is only signal, the increased flexibility of the SDID estimator allows it (alongside MC) to outperform the SC and DID estimators dramatically.
\begin{figure}[t]
\begin{center}
\includegraphics[width=.7\textwidth]{error-density-plot-1.pdf}
\caption{Distribution of the errors of SDID, SC and DID in the setting of the ``baseline'' (i.e., with minimum wage)
and random assignment rows of Table \ref{table1}.}
\label{fig_mw_rand}
\end{center}
\end{figure}
Next, we focus on two designs of interest: One with the assignment probability model based on parameters estimated in the minimum
wage law model and one where the treatment exposure $D_i$ is assigned uniformly at random.
Figure \ref{fig_mw_rand} shows the errors of the DID, SC and SDID estimators in both settings, and reinforces our observations above.
When assignment is not uniformly random, the distribution of the DID errors is visibly off-center, showing the bias of the estimator. In contrast,
the errors from SDID are nearly centered. Meanwhile, when treatment assignment is uniformly random, both estimators are centered but the errors
of DID are more spread out. We note that the right panel of Figure \ref{fig_mw_rand} is closely related to the simulation
specification of \citet*{Bertrand2004did}. From this perspective, \citet{Bertrand2004did} correctly argue that the error
distribution of DID is centered, and that the error scale can accurately be recovered using appropriate robust estimators.
Here, however, we go further and show that this noise can be substantially reduced by using an estimator like SDID
that can exploit predictable variation by matching on pre-exposure trends.
Finally, we note that Figure \ref{fig_mw_rand} shows that the error distribution of SDID is nearly unbiased and Gaussian in both designs,
thus suggesting that it should be possible to use $\hat{\tau}^{\rm sdid}$ as the basis for valid inference. We postpone a discussion
of confidence intervals until Section \ref{sec:inference}, where we consider various strategies for inference based on SDID
and show that they attain good coverage here.
\subsection{Penn World Table Placebo Study}
\label{sec:penn}
The simulation based on the CPS is a natural benchmark for applications that traditionally rely on DID-type methods to estimate the policy effects. In contrast, SC methods are often used in applications where
units tend to be more heterogeneous and are observed over a longer timespan as in, e.g., \citet*{abadie2014}.
To investigate the behavior of SDID in this type of setting, we propose a second set of simulations based on the Penn World
Table. This dataset contains observations on annual real GDP for $N = 111$ countries for $T = 48$ consecutive years, starting from 1959;
we end the dataset in 2007 because we do not want the treatment period to coincide with the Great Recession.
We construct the outcome and the assignment model following the same
procedure outlined in the previous subsection. We select $\log(\mathrm{real\ GDP})$ as
the primary outcome. As with the CPS dataset, the two-way fixed effects explain most of the variation; however, the interactive
component plays a larger role in determining outcomes for this dataset than for the CPS data. We again derive treatment assignment via an exposure
variable $D_i$, and consider both a uniformly random distribution for $D_i$ as well as two non-uniform ones
based on predicting Penn World Table indicators of democracy and education respectively.
Results of the simulation study are presented in Table \ref{table2}. At a high level, these results mirror the ones above: SDID
again performs well in terms of both bias and root-mean squared error and across all simulation settings dominates the other estimators.
In particular, SDID is nearly unbiased, which is important for constructing confidence intervals with accurate coverage rates. The main difference
between Tables \ref{table1} and \ref{table2} is that DID does substantially worse here relative to SC than before. This appears to be due to
the presence of a stronger interactive component in the Penn World Table dataset, and is in line with the empirical practice of preferring
SC over DID in settings of this type. We again defer a discussion of inference to Section \ref{sec:inference}.
\begin{table}[t]
\begin{center}
\begin{adjustbox}{width=1\textwidth}
\begin{tabular}{|l|rrrr|rrrrr|rrrrr|}
\hline
& \multirow{2}{*}{$\frac{\|\boldsymbol{F}\|_F}{\sqrt{NT}}$} & \multirow{2}{*}{$\frac{\|\boldsymbol{M}\|_F}{\sqrt{NT}}$} & \multirow{2}{*}{$\sqrt{\frac{\operatorname{tr}(\Sigma)}{T}}$} & \multirow{2}{*}{AR(2)} &
\multicolumn{5}{c|}{RMSE}& \multicolumn{5}{c|}{Bias}\\
& && & &SDID & SC & DID &MC& DIFP & SDID & SC & DID &MC& DIFP\\ \hline
Democracy & 0.972 & 0.229 & 0.070 & (.91, -.22) & \bf 0.031 & 0.038 & 0.197 & 0.058 & 0.039 & -0.005 & -0.004 & 0.175 & 0.043 & -0.007\\
Education & 0.972 & 0.229 & 0.070 & (.91, -.22) & \bf 0.030 & 0.053 & 0.172 & 0.049 & 0.039 & -0.003 & 0.025 & 0.162 & 0.040 & -0.005 \\
Random & 0.972 & 0.229 & 0.070 & (.91, -.22) & \bf 0.037 & 0.046 & 0.129 & 0.063 & 0.045 & -0.002 & -0.011 & -0.006 & -0.004 & -0.004 \\
\hline
\end{tabular}
\end{adjustbox}
\caption{Simulation results based on the the Penn World Table dataset.
We use $\log(GDP)$ as the outcome, with $N_{\rm tr} = 10$ out of $N = 111$ treatment countries,
and $T_{\rm post} = 10$ out of $T = 48$ treatment periods.
In the first two rows we consider treatment assignment distributions based on democracy status and education metrics,
while in the last row the treatment is assigned completely at random.
All results are based on 1000 simulations.
\label{table2}}
\end{center}
\end{table}
\section{Formal Results}
\label{sec:formal}
In this section we discuss the formal results.
For the remainder of the paper, we assume that the data generating process follows a generalization of the latent factor model \eqref{basic_model},
\begin{equation}
\label{eq:model}
\boldsymbol{Y} = \boldsymbol{L} + \boldsymbol{W} \circ \boldsymbol{\tau} + \boldsymbol{E}, \qquad \mathrm{where} \qquad (\boldsymbol{W} \circ \boldsymbol{\tau})_{it} = \boldsymbol{W}_{it} \boldsymbol{\tau}_{it}.
\end{equation}
The model allows for heterogeneity in treatment effects $\tau_{it}$, as in
\citet{de2020two}. As above, we
assume block assignment $W_{it} = 1\left(\left\{i > N_{\rm co}, \, t > T_{\rm pre}\right\}\right)$, where
the subscript ``${\rm co}$'' stands for control group,
``${\rm tr}$'' stands for treatment group, ``${\rm pre}$'' stands for pre-treatment,
and ``${\rm post}$'' stands for post-treatment.
It is useful to characterize the systematic component $\boldsymbol{L}$ as a factor model $\boldsymbol{L}=\boldsymbol{\Gamma} \boldsymbol{\Upsilon}^{\top}$
as in \eqref{basic_model}, where we define factors $\boldsymbol{\Gamma}=\boldsymbol{UD}^{1/2}$ and $\boldsymbol{\Upsilon}^{\top}=\boldsymbol{D}^{1/2}\boldsymbol{V}^{\top}$ in terms of the singular value decomposition $\boldsymbol{L}=\boldsymbol{U D V}^{\top}$.
Our target estimand is the average treatment effect for the treated units during the periods they were treated,
which under block assignment is
\begin{equation}
\label{eq:estimand}
\tau = \frac{1}{N_{\rm tr} T_{\rm post}} \sum_{i=N_{\rm co}+1}^N\ \sum_{t=T_{\rm pre}+1}^T \boldsymbol{\tau}_{it}.
\end{equation}
For notational convenience, we partition the matrix $\boldsymbol{Y}$ as
\[ \boldsymbol{Y} =\left(
\begin{array}{cc}
\boldsymbol{Y}_{{\rm co},{\rm pre}} &
\boldsymbol{Y}_{{\rm co},{\rm post}}\\
\boldsymbol{Y}_{{\rm tr},{\rm pre}} &
\boldsymbol{Y}_{{\rm tr},{\rm post}}
\end{array}\right),\]
with $\boldsymbol{Y}_{{\rm co},{\rm pre}}$ a $N_{\rm co}\times T_{\rm pre}$ matrix,
$\boldsymbol{Y}_{{\rm co},{\rm post}}$ a $N_{\rm co}\times T_{\rm post}$ matrix,
$\boldsymbol{Y}_{{\rm tr},{\rm pre}}$ a $N_{\rm tr}\times T_{\rm pre}$ matrix, and
$\boldsymbol{Y}_{{\rm tr},{\rm post}}$ a $N_{\rm tr}\times T_{\rm post}$ matrix,
and similar for $\boldsymbol{L}$, $\boldsymbol{W}$, $\tau$, and $\boldsymbol{E}$.
Throughout our analysis, we will assume that the errors $\boldsymbol{E}_{i.}$ are homoskedastic across units (but not across time),
{\it i.e.}, that $\operatorname{Var}_{}\left[\boldsymbol{E}_{i.}\right] = \Sigma \in \mathbb{R}^{T \times T}$ for all units $i = 1, \, \ldots, \, n$. We partition $\Sigma$ as
\[ \Sigma=\left(
\begin{array}{cc} \Sigma_{{\rm pre},{\rm pre}} & \Sigma_{{\rm pre},{\rm post}}\\
\Sigma_{{\rm post},{\rm pre}} & \Sigma_{{\rm post},{\rm post}}\end{array}
\right).\]
Given this setting, we are interested in guarantees on how accurately SDID can recover $\tau$.
A simple, intuitively appealing approach to estimating $\tau$ in \eqref{eq:model} is to directly fit both $\boldsymbol{L}$ and $\tau$
via methods for low-rank matrix estimation, and several variants of this approach have been proposed in the literature
\citep*[e.g.,][]{athey2017matrix,bai2009panel,xu2017generalized, agarwal2019robustness}. However, our main interest is in $\tau$ and not in $\boldsymbol{L}$,
and so one might suspect that approaches that provide consistent estimation of $\boldsymbol{L}$ may rely on assumptions
that are stronger than what is necessary for consistent estimation of $\tau$.
Synthetic control methods address confounding bias without explicitly estimatin $\boldsymbol{L}$ in \eqref{eq:model}.
Instead, they take an indirect approach more akin to balancing as
in \citet{zubizarreta2015stable} and \citet*{athey2018approximate}. Recall that the SC weights
$\hat{\omega}^{\rm sc}$ seek to balance out the pre-intervention trends in $\boldsymbol{Y}$. Qualitatively, one might hope that doing so also leads us to balance out the unit-factors $\boldsymbol{\Gamma}$ from \eqref{basic_model}, rendering
$\sum_{i = N_{\rm co} + 1}^N \hat{\omega}_i^{\rm sc} \boldsymbol{\Gamma}_{i.} - \sum_{i = 1}^{N_{\rm co}} \hat{\omega}_i^{\rm sc} \boldsymbol{\Gamma}_{i.} \approx 0$.
\citet*{Abadie2010} provide some arguments for why this should be the case, and our formal analysis outlines
a further set of conditions under which this type of phenomenon holds. Then, if $\hat{\omega}^{\rm sc}$ in fact succeeds
in balancing out the factors in $\boldsymbol{\Gamma}$, the SC estimator can be approximated as $\hat{\tau}^{\rm sc} \approx \tau + \sum_{i = 1}^N (2W_i -1) \hat{\omega}_i^{\rm sc} \bar \varepsilon_{i}$ with
$\bar \varepsilon_i=
T_{\rm post}^{-1} \sum_{t=T_{pre}+1}^T \varepsilon_{it}$
; in words, SC weighting has succeeded in removing the bias associated with the systematic component $\boldsymbol{L}$
and in delivering a nearly unbiased estimate of $\tau$.
Much like the SC estimator, the SDID estimator seeks to recover $\tau$ in \eqref{eq:model} by reweighting to remove the bias associated with $\boldsymbol{L}$. However, the SDID estimator takes a
two--pronged approach. First, instead of only
making use of unit weights $\hat{\omega}$ that can be used to balance out $\boldsymbol{\Gamma}$, the estimator also incorporates time weights $\hat{\lambda}$
that seek to balance out $\boldsymbol{\Upsilon}$. This provides a type of double robustness property, whereby if
one of the balancing approaches is effective, the dependence on $\boldsymbol{L}$ is approximately removed.
Second, the use of two-way fixed effects in \eqref{main_sdid} and intercept terms in \eqref{unit_weights} and
\eqref{time_weights} makes the SDID estimator invariant to additive shocks to any row or column, i.e.,
if we modify $\boldsymbol{L}_{it} \leftarrow \boldsymbol{L}_{it} + \alpha_i + \beta_t$ for any choices $\alpha_i$ and $\beta_t$ the
estimator $\hat{\tau}^{\rm sdid}$ remains unchanged. The estimator shares this invariance property with DID
(but not SC).\footnote{More specifically, as suggested by \eqref{main_sc}, SC is invariant to shifts in $\beta_t$
but not $\alpha_i$. In this context, we also note that the DIFP estimator proposed by \citet{doudchenko2016balancing}
and \citet{ferman2019synthetic} that center each unit's trajectory before applying the synthetic control method is also invariant to shifts in $\alpha_i$.}
The goal of our formal analysis is to understand how and when the SDID weights
succeed in removing the bias due to $\boldsymbol{L}$. As discussed below, this requires assumptions
on the signal to noise ratio. The assumptions require that $\boldsymbol{E}$ does not incorporate too much serial correlation within units, so that we can
attribute persistent patterns in $\boldsymbol{Y}$ to patterns in $\boldsymbol{L}$; furthermore, $\boldsymbol{\Gamma}$ should be stable over time, particularly through the treatment periods. Of course, these are non-trivial assumptions.
However, as discussed further in Section \ref{sec:relworks}, they are considerably weaker than what is required in results of
\citet{bai2009panel} or \citet{moon2015linear,moon2017dynamic} for methods that require explicitly estimating $\boldsymbol{L}$ in \eqref{eq:model}.
Furthermore, these assumption are aligned with standard practice in the literature; for example, we can assess the claim that we balance
all components of $\boldsymbol{\Gamma}$ by examining the extent to which the method succeeds in balancing pre-intervention periods. Historical context may be needed to justify the assumption that that there were no
other shocks disproportionately affecting the treatment units at the time of the treatment.
\subsection{Weighted Double-Differencing Estimators}
\label{sec:double_diff}
We introduced the SDID estimator \eqref{main_sdid} as the
solution to a weighted two-way fixed effects regression. For the purpose of our formal results,
however, it is convenient to work with the alternative characterization described above
in Equation~\ref{eq:tau-weighted-differences}.
For any weights $\omega \in \Omega$ and $\lambda \in \Lambda$, we can define a weighted
double-differencing estimator\footnote{This weighted double-differencing structure plays a key
role in understanding the behavior of SDID. As discussed further in Section \ref{sec:relworks},
despite relying on a different motivation, certain specifications of the recently proposed ``augmented synthetic
control'' method of \citet*{ben2018augmented} also result in a weighted double-differencing estimator.}
\begin{equation}
\label{eq:tau-weighted-differences}
\hat{\tau}(\omega,\lambda) = \omega_{{\rm tr}}^\top \boldsymbol{Y}_{{\rm tr},{\rm post}} \lambda_{{\rm post}} - \omega_{{\rm co}}^\top \boldsymbol{Y}_{{\rm co},{\rm post}} \lambda_{{\rm post}} - \omega_{{\rm tr}}^\top \boldsymbol{Y}_{{\rm tr},{\rm pre}} \lambda_{{\rm pre}} +
\omega_{{\rm co}}^\top \boldsymbol{Y}_{{\rm co},{\rm pre}} \lambda_{{\rm pre}}.
\end{equation}
One can verify that the basic DID estimator is of the form \eqref{eq:tau-weighted-differences}, with
constant weights $\omega_{\rm tr} = 1/N_{\rm tr}$, etc. The proposed SDID estimator \eqref{main_sdid} can also be written as
\eqref{eq:tau-weighted-differences}, but now with weights $\hat{\omega}^{\rm sdid}$ and $\hat{\lambda}^{\rm sdid}$ solving \eqref{unit_weights}
and \eqref{time_weights} respectively. When there is no risk of ambiguity, we will omit the SDID-superscript from the weights
and simply write $\hat{\omega}$ and $\hat{\lambda}$.
Now, note that for any choice of weights $\omega \in \Omega$ and $\lambda \in \Lambda$, we have $\omega_{{\rm tr}} \in \mathbb{R}^{N_{\rm tr}}$
and $\lambda_{{\rm post}} \in \mathbb{R}^{T_{\rm post}}$ with all elements equal to $1/N_{\rm tr}$ and $1/T_{\rm post}$ respectively, and so
$\omega_{{\rm tr}}^\top \boldsymbol{\tau}_{{\rm tr},{\rm post}} \lambda_{{\rm post}} = \tau$. Thus, we can decompose the error of any weighted double-differencing
estimator with weights satisfying these conditions as the sum of a bias and a noise component:
\begin{equation}
\label{eq:sdid_error}
\begin{split}
\hat{\tau}(\omega,\lambda) - \tau
&= \underbrace{\omega_{{\rm tr}}^\top \boldsymbol{L}_{{\rm tr},{\rm post}} \lambda_{{\rm post}} - \omega_{\rm co}^\top \boldsymbol{L}_{{\rm co},{\rm post}} \lambda_{{\rm post}} - \omega_{{\rm tr}}^\top \boldsymbol{L}_{{\rm tr},{\rm pre}} \lambda_{\rm pre} + \omega_{\rm co}^\top \boldsymbol{L}_{{\rm co},{\rm pre}} \lambda_{\rm pre}}_{\text{bias}\ B(\omega,\lambda)} \\
&\quad\quad\quad\quad + \underbrace{\omega_{{\rm tr}}^\top \boldsymbol{E}_{{\rm tr},{\rm post}} \lambda_{{\rm post}} - \omega_{\rm co}^\top \boldsymbol{E}_{{\rm co},{\rm post}} \lambda_{{\rm post}} - \omega_{{\rm tr}}^\top\boldsymbol{E}_{{\rm tr},{\rm pre}} \lambda_{\rm pre} +
\omega_{\rm co}^\top \boldsymbol{E}_{{\rm co},{\rm pre}} \lambda_{\rm pre}.}_{\text{noise}\ \varepsilon(\omega,\lambda)}
\end{split}
\end{equation}
In order to characterize the distribution of $\hat{\tau}^{\rm sdid} - \tau$, it thus remains to carry out two tasks. First, we need
to understand the scale of the errors $ B(\omega,\lambda)$ and $\varepsilon(\omega,\lambda)$, and second, we need
to understand how data-adaptivity of the weights $\hat{\omega}$ and $\hat{\lambda}$ affects the situation.
\subsection{Oracle and Adaptive Synthetic Control Weights}
\label{sec:oracle-and-adaptive}
To address the adaptivity of the SDID weights $\hat{\omega}$ and $\hat{\lambda}$ chosen via \eqref{unit_weights}
and \eqref{time_weights}, we construct alternative ``oracle'' weights that have similar properties to $\hat{\omega}$ and $\hat{\lambda}$
in terms of eliminating bias due to $\boldsymbol{L}$, but are deterministic. We can then further decompose the error of $\hat{\tau}^{\rm sdid}$ into
the error of a weighted double-differencing estimator with the oracle weights and the difference between the oracle and feasible
estimators. Under appropriate conditions, we find the latter term negligible relative to the error of the
oracle estimator, opening the door to a simple asymptotic characterization of the error distribution of $\hat{\tau}^{\rm sdid}$.
We define such oracle weights $\tilde{\omega}$ and $\tilde{\lambda}$ by minimizing the expectation of the objective functions
$\ell_{unit}(\cdot)$ and $\ell_{time}(\cdot)$ used in \eqref{unit_weights} and \eqref{time_weights} respectively, and set
\begin{equation}
\label{eq:oracle_gen}
\begin{split}
&\left(\tilde{\omega}_0, \tilde{\omega}\right) = \operatorname*{arg\,min}_{\omega_0 \in \mathbb{R}, \omega \in \Omega} \mathbb{E}_{}\left[\ell_{unit}(\omega_0, \omega)\right], \ \ \ \
\left(\tilde{\lambda}_0, \, \tilde{\lambda}\right) = \operatorname*{arg\,min}_{\lambda_0 \in \mathbb{R}, \lambda \in \Lambda}\mathbb{E}_{}\left[\ell_{time}(\lambda_0, \lambda)\right].
\end{split}
\end{equation}
In the case of our model \eqref{eq:model} these weights admit a simplified characterization
\begin{align}
&\left(\tilde{\omega}_0, \tilde{\omega}\right)
= \operatorname*{arg\,min}_{\omega_0 \in \mathbb{R}, \omega \in \Omega} \left\lVert\omega_0 + \omega_{{\rm co}}^{\top} \boldsymbol{L}_{{\rm co},{\rm pre}} - \omega_{{\rm tr}}^{\top} \boldsymbol{L}_{{\rm tr}, {\rm pre}}\right\rVert^2_2 + \left(\operatorname{tr}(\Sigma_{{\rm pre},{\rm pre}}) + \zeta^2 T_{\rm pre}\right) \left\lVert\omega\right\rVert_2^2,
\label{eq:row_oracle} \\
&\left(\tilde{\lambda}_0, \, \tilde{\lambda}\right)
= \operatorname*{arg\,min}_{\lambda_0 \in \mathbb{R}, \lambda \in \Lambda} \left\lVert\lambda_0 + \boldsymbol{L}_{{\rm co},{\rm pre}} \lambda_{{\rm pre}} - \boldsymbol{L}_{{\rm co},{\rm post}}\lambda_{{\rm post}}\right\rVert_2^2 + \left\lVert\tilde\Sigma \lambda\right\rVert_2^2, \
\label{eq:col_oracle} \\
& \quad\quad\quad\quad\quad\quad \text{ where }\ \ \ \tilde\Sigma=\begin{pmatrix}\ \Sigma_{{\rm pre},{\rm pre}} & -\Sigma_{{\rm pre},{\rm post}} \\
-\Sigma_{{\rm post},{\rm pre}} & \ \ \ \Sigma_{{\rm post},{\rm post}}\end{pmatrix}. \nonumber
\end{align}
The error of the synthetic difference in differences estimator can now be decomposed as follows,
\begin{equation}
\label{eq:full_decomp}
\hat{\tau}^{\rm sdid} - \tau
= \underbrace{\varepsilon(\tilde{\omega},\tilde{\lambda})}_{\text{oracle noise}}
+ \underbrace{B(\tilde{\omega},\tilde{\lambda})}_{\text{oracle confounding bias}}
+\quad \underbrace{\hat{\tau}(\hat{\omega},\hat{\lambda}) - \hat{\tau}(\tilde{\omega},\tilde{\lambda}),}_{\text{deviation from oracle}}
\end{equation}
and our task is to characterize all three terms.
First, the oracle noise term tends to be small when the weights are not too concentrated,
{\it i.e.,} when $\lVert\tilde{\omega}\rVert_2$ and $\lVert\tilde{\lambda}\rVert_2$ are small, and we have a sufficient number of exposed units and time periods.
In the case with $\Sigma = \sigma^2 I_{T \times T}$, {\it i.e.,} without any cross-observation correlations, we note that
$\operatorname{Var}_{}\left[\varepsilon(\tilde{\omega},\tilde{\lambda})\right] = \sigma^2 \left(N_{\rm tr}^{-1} + \lVert\tilde{\omega}\rVert_2^2\right)\left(T_{\rm post}^{-1} + \lVert\tilde{\lambda}\rVert_2^2\right)$.
When we move to our asymptotic analysis below, we work under assumptions that make this oracle noise
term dominant relative to the other error terms in \eqref{eq:full_decomp}.
Second, the oracle confounding bias will be small either when the pre-exposure oracle row regression
fits well and generalizes to the exposed rows, i.e.,
$\tilde{\omega}_0 + \tilde{\omega}_{\rm co}^\top \boldsymbol{L}_{{\rm co},{\rm pre}} \approx \tilde{\omega}_{{\rm tr}}^\top \boldsymbol{L}_{{\rm tr},{\rm pre} }$ and
$\tilde{\omega}_0 + \tilde{\omega}_{\rm co}^\top \boldsymbol{L}_{{\rm co},{\rm post}} \approx \tilde{\omega}_{{\rm tr}}^{\top} \boldsymbol{L}_{{\rm tr},{\rm post}}$,
or when the unexposed oracle column regression fits well and generalizes
to the exposed columns,
$\tilde{\lambda}_0 + \boldsymbol{L}_{{\rm co},{\rm pre}}\tilde{\lambda}_{\rm pre} \approx \boldsymbol{L}_{{\rm co},{\rm post}}\tilde{\lambda}_{{\rm post}}$ and
$\tilde{\lambda}_0 + \boldsymbol{L}_{{\rm tr},{\rm pre}}\tilde{\lambda}_{\rm pre} \approx \boldsymbol{L}_{{\rm tr},{\rm post}} \tilde{\lambda}_{{\rm post}}$.
Moreover, even if neither model generalizes sufficiently well on its own, it suffices for one model to predict the generalization error of the other:
\begin{align*}
B(\omega,\lambda) &= (\omega_{{\rm tr}}^\top \boldsymbol{L}_{{\rm tr},{\rm post}} - \omega_{{\rm co}}^\top \boldsymbol{L}_{{\rm co},{\rm post}})\lambda_{{\rm post}} -
(\omega_{{\rm tr}}^\top \boldsymbol{L}_{{\rm tr},{\rm pre}} - \omega_{{\rm co}}^\top \boldsymbol{L}_{{\rm co},{\rm pre}})\lambda_{{\rm pre}} \\
&= \omega_{{\rm tr}}^\top (\boldsymbol{L}_{{\rm tr},{\rm post}}\lambda_{{\rm post}} - \boldsymbol{L}_{{\rm tr},{\rm pre}}\lambda_{{\rm pre}}) -
\omega_{{\rm co}}^\top(\boldsymbol{L}_{{\rm co},{\rm post}}\lambda_{{\rm post}} - \boldsymbol{L}_{{\rm co},{\rm pre}}\lambda_{\rm pre}).
\end{align*}
The upshot is even if one of the sets of weights fails to remove the bias from the presence of $\boldsymbol{L}$, the combination of weights $\tilde{\omega}$ and $\tilde{\lambda}$ can
compensate for such failures. This double robustness property is similar to that of the augmented inverse
probability weighting estimator, whereby one can trade off between accurate estimates of the outcome
and treatment assignment models \citep*{ben2018augmented, scharfstein1999adjusting}.
We note that although poor fit in the oracle regressions on the unexposed rows and columns of $\boldsymbol{L}$ will often be indicated by
a poor fit in the realized regressions on the unexposed rows and columns of $\boldsymbol{Y}$,
the assumption that one of these regressions generalizes to exposed rows or columns
is an identification assumption without clear testable implications.
It is essentially an assumption of no unexplained confounding:
any exceptional behavior of the exposed observations, whether due to exposure or not, can be ascribed to it.
Third, our core theoretical claim, formalized in our asymptotic analysis, is that the SDID estimator will be close to the oracle
when the oracle unit and time weights look promising on their respective training sets, i.e,
when $\tilde{\omega}_0 + \tilde{\omega}_{\rm co}^{\top} \boldsymbol{L}_{{\rm co},{\rm pre}} \approx \tilde{\omega}_{{\rm tr}}^{\top} \boldsymbol{L}_{{\rm tr},{\rm pre}}$ and $\lVert\tilde{\omega}\rVert_2$ is not too large and $\tilde{\lambda}_0 + \boldsymbol{L}_{{\rm co},{\rm pre}}\tilde{\lambda}_{\rm pre} \approx \boldsymbol{L}_{{\rm co},{\rm post}}\tilde{\lambda}_{{\rm post}}$ and $\lVert\tilde{\lambda}\rVert_2$ is not too large. Although the details differ, as described above
these qualitative properties are also criteria for accuracy of the oracle estimator itself.
Finally, we comment briefly on the behavior of the oracle time weights $\tilde{\lambda}$ in the presence of
autocorrelation over time. When $\Sigma$ is not diagonal, the effective regularization term in \eqref{eq:col_oracle}
does not shrink $\tilde{\lambda}_{{\rm pre}}$ towards zero, but rather toward an autoregression vector
\begin{equation}
\label{eq:autoregression-vector}
\psi = \operatorname*{arg\,min}_{v \in \mathbb{R}^{T_{\rm pre}}} \left\lVert\tilde\Sigma \begin{pmatrix}v \\ \lambda_{{\rm post}} \end{pmatrix}\right\rVert
= \Sigma_{{\rm pre},{\rm pre}}^{-1}\Sigma_{{\rm pre},{\rm post}}\lambda_{{\rm post}}.
\end{equation}
Here $\lambda_{{\rm post}}$ is the $T_{\rm post}$-component column vector with all elements equal to $1/T_{\rm post}$
and $\psi$ is the population regression coefficient in a regression of the average of the post-treatment errors on the pre-treatment errors.
In the absence of autocorrelation, $\psi$ is zero, but when autocorrelation is present, shinkage toward $\psi$
reduces the variance of the SDID estimator---and enables us to gain precision over the basic DID estimator \eqref{main_did} even when the two-way fixed effects model is correctly specified. This explains some of the behavior noted in the simulations.
\subsection{Asymptotic Properties}
\label{section:asym}
To carry out the analysis plan sketched above, we need to embed our problem into an asymptotic setting.
First, we require the error matrix $\boldsymbol{E}$ to satisfy some regularity properties.
\begin{assumption}\label{ass:noise}{\sc (Properties of Errors)}
The rows $\boldsymbol{E}_{i}$ of the noise matrix are independent and identically distributed Gaussian
vectors and the eigenvalues of its covariance matrix $\Sigma$ are bounded and bounded away from zero.
\end{assumption}
Next, we spell out assumptions about the sample size. At a high level, we want the panel to be large
({\it i.e.}, $N, \, T \rightarrow \infty$), and for the number of treated cells of the panel to grow to infinity but
slower than the total panel size. We note in particular that we can accommodate sequences where
one of $T_{\rm post}$ or $N_{\rm tr}$ is fixed, but not both.
\begin{assumption}\label{ass:sample_sizes} {\sc (Sample Sizes)}
We consider a sequence of populations where \\
$(i)$ the product $N_{\rm tr} \, T_{\rm post}$ goes to infinity, and both $N_{\rm co}$ and $T_{\rm pre}$ go to infinity,\\
$(ii)$ the ratio $T_{\rm pre}/N_{\rm co}$ is bounded and bounded away from zero,\\
$(iii)$ $N_{\rm co} / (N_{\rm tr} T_{\rm post} \max(N_{\rm tr}, T_{\rm post})\log^2(N_{\rm co})) \to \infty$.
\end{assumption}
We also need to make assumptions about the spectrum of $\boldsymbol{L}$; in particular, $\boldsymbol{L}$ cannot have too many large singular values,
although we allow for the possibility of many small singular values. A sufficient, but not necessary, condition for the assumption
below is that the rank of $\boldsymbol{L}$ is less than $\sqrt{\min(T_{\rm pre},N_{\rm co})}$. Notice that we do not assume any
lower bounds for non-zero singular values of $\boldsymbol{L}$; in fact can accommodate arbitrarily many non-zero but
very small singular values, much like, {\it e.g.}, \citet*{belloni2014inference} can accommodate arbitrarily many
non-zero but very small signal coefficients in a high-dimensional inference problem.
We need that the $\sqrt{\min(T_{\rm pre},N_{\rm co})}$th singular value of $\boldsymbol{L}_{{\rm co},{\rm pre}}$ is sufficiently small. Formally:
\begin{assumption}\label{rank}{\sc (Properties of $\boldsymbol{L}$)}
Letting $\sigma_1(\boldsymbol{\Gamma}), \sigma_2(\boldsymbol{\Gamma}), \ldots$ denote the singular values of the matrix $\boldsymbol{\Gamma}$ in decreasing order
and $R$ the largest integer less than $\sqrt{\min(T_{\rm pre},N_{\rm co})}$,
\begin{equation}
\sigma_{R}(\boldsymbol{L}_{{\rm co},{\rm pre}})/R = o\left(\min\Big\{N_{\rm tr}^{-1/2}\log^{-1/2}(N_{\rm co}), T_{\rm post}^{-1/2}\log^{-1/2}(T_{\rm pre})\Bigr\}\right)
\end{equation}
\end{assumption}
The last---and potentially most interesting---of our assumptions concerns the relation between the factor
structure $\boldsymbol{L}$ and the assignment mechanism $\boldsymbol{W}$. At a high level, it plays the role of an identifying assumption,
and guarantees that the oracle weights from \eqref{eq:row_oracle} and \eqref{eq:col_oracle} that are directly
defined in terms of $\boldsymbol{L}$ are able to adequately cancel out $\boldsymbol{L}$ via the weighted double-differencing strategy.
This requires that the optimization problems \eqref{eq:row_oracle} and \eqref{eq:col_oracle} accommodate
reasonably dispersed weights, and that the treated units and after periods not be too dissimilar from
the control units and the before periods respectively.
\begin{assumption}\label{weightss}{\sc (Properties of Weights and $\boldsymbol{L}$)}
The oracle unit weights
$\tilde\omega$ satisfy
\begin{equation}
\label{eq:weightss1}
\begin{aligned}
&\|\tilde\omega_{{\rm co}}\|_2 = o([(N_{\rm tr} T_{\rm post})\log(N_{\rm co})]^{-1/2}) \qquad \text{ and } \\
&\lVert\tilde\omega_0 + \tilde\omega_{{\rm co}}^{\top}\boldsymbol{L}_{{\rm co},{\rm pre}} - \tilde\omega_{{\rm tr}}^{\top}\boldsymbol{L}_{{\rm tr},{\rm pre}}\rVert_2 = o( N_{\rm co}^{1/4} (N_{\rm tr} T_{\rm post}\max(N_{\rm co},T_{\rm post}))^{-1/4} \log^{-1/2}(N_{\rm co})),
\end{aligned}
\end{equation}
the oracle time weights $\tilde\lambda$ satisfy
\begin{equation}
\label{eq:weightss2}
\begin{aligned}
&\|\tilde\lambda_{{\rm pre}}-\psi\|_2 = o([(N_{\rm tr} T_{\rm post})\log(N_{\rm co})]^{-1/2}) \qquad \text{ and } \\
&\lVert\tilde\lambda_0 + \boldsymbol{L}_{{\rm co},{\rm pre}}\tilde\lambda_{\rm pre} - \boldsymbol{L}_{{\rm co},{\rm post}}\tilde\lambda_{\rm post}\rVert_2 = o(N_{\rm co}^{1/4} (N_{\rm tr} T_{\rm post})^{-1/8}),
\end{aligned}
\end{equation}
and the oracle weights jointly satisfy
\begin{equation}
\label{eq:weightss3}
\tilde\omega_{\rm tr}^\top\boldsymbol{L}_{{\rm tr},{\rm post}}\tilde\lambda_{\rm post}
-\tilde\omega_{\rm co}^\top\boldsymbol{L}_{{\rm co},{\rm post}}\tilde\lambda_{\rm post}
-
\tilde\omega_{\rm tr}^\top\boldsymbol{L}_{{\rm tr},{\rm pre}}\tilde\lambda_{\rm pre}
+\tilde\omega_{\rm co}^\top\boldsymbol{L}_{{\rm co},{\rm pre}}\tilde\lambda_{\rm pre}
=o\left((N_{\rm tr} T_{\rm post})^{-1/2}\right).
\end{equation}
\end{assumption}
Assumptions \ref{ass:noise}-\ref{weightss} are substantially weaker than those used to establish asymptotic
normality of comparable methods.\footnote{In particular, note that our assumptions are
satisfied in the well-specified two-way fixed effect setting model. Suppose we have \smash{$L_{it}=\alpha_i+\beta_t$}
with uncorrelated and homoskedastic errors, and that the sample size restrictions in Assumption \ref{ass:sample_sizes} are satisfied.
Then Assumption \ref{ass:noise} is automatically satisfied, and
the rank condition on $\boldsymbol{L}$ from Assumption \ref{rank} is satisfied with $R=2$.
Next, we see that the oracle unit weights satisfy \smash{$\tilde\omega_{{\rm co},i}=1/N_{\rm co}$} so that
\smash{$\|\tilde\omega\|_2=1/\sqrt N_{\rm co}$}, and the oracle time weights satisfy \smash{$\tilde\lambda_{{\rm pre},i}=1/T_{\rm pre}$}
so that \smash{$\|\tilde\lambda-\psi\|_2=1/\sqrt N_{\rm co}$}. Thus if the restrictions on the rates at which the sample sizes increase
in Assumption \ref{ass:sample_sizes} are satisfied, then \eqref{eq:weightss1} and \eqref{eq:weightss2} are satisfied.
Finally, the additive structure of $\boldsymbol{L}$ implies that, as long as the weights for the controls sum to one,
\smash{$\tilde\omega_{\rm tr}^\top\boldsymbol{L}_{{\rm tr},{\rm post}}\tilde\lambda_{\rm post}
-\tilde\omega_{\rm co}^\top\boldsymbol{L}_{{\rm co},{\rm post}}\tilde\lambda_{\rm post}=0$},
and
\smash{$\tilde\omega_{\rm tr}^\top\boldsymbol{L}_{{\rm tr},{\rm pre}}\tilde\lambda_{\rm pre}
+\tilde\omega_{\rm co}^\top\boldsymbol{L}_{{\rm co},{\rm pre}}\tilde\lambda_{\rm pre}
=0$},
so that \eqref{eq:weightss3} is satisfied.}
We do not require that double differencing alone removes the individual and time effects
as the DID assumptions do. Furthermore, we do not require that unit comparisons alone are sufficient to remove the biases
in comparisons between treated and control units as the SC assumptions do. Finally, we do not require a low rank factor model
to be correctly specified, as is often assumed in the analysis of methods that estimate $\boldsymbol{L}$ explicitly
\citep[e.g.,][]{bai2009panel,moon2015linear,moon2017dynamic}. Rather, we only need the combination of the
three bias-reducing components in the SDID estimator, $(i)$ double differencing, $(ii)$ the unit weights, and $(iii)$ the time
weights, to reduce the bias to a sufficiently small level.
Our main formal result states that under these assumptions, our estimator is asymptotically normal.
Furthermore, its asymptotic variance is optimal, coinciding with the variance we would get
if we knew $\boldsymbol{L}$ and $\Sigma$ a-priori and could therefore estimate $\tau$ by a simple average
of $\tau_{it}$ plus unpredictable noise,
$N_{{\rm tr}}^{-1}\sum_{i=N_{{\rm co}}+1}^N [T_{{\rm post}}^{-1}\sum_{t=T_{\rm pre}+1}^T (\boldsymbol{\tau}_{it} + \varepsilon_{it}) - \boldsymbol{E}_{i,{\rm pre}}\psi]$.
\begin{theorem}
\label{theo:asymptotic-linearity}
Under the model \eqref{eq:model} with $\boldsymbol{L}$ and $\boldsymbol{W}$ taken as fixed, suppose that we run the
SDID estimator \eqref{main_sdid} with regularization parameter $\zeta$
satisfying
$(N_{{\rm tr}}T_{{\rm post}})^{1/2}
\log(N_{\rm co}) = o(\zeta^2)$.
Suppose moreover that Assumptions \ref{ass:noise}-\ref{weightss} hold. Then,
\begin{equation}
\label{eq:asymptotic-linearity}
\hat{\tau}^{\rm sdid}-\tau=\frac{1}{N_{\rm tr}}\sum_{i=N_{{\rm co}}+1}^N \left(\frac{1}{T_{\rm post}}\sum_{t=T_{\rm pre}+1}^T \varepsilon_{it} - \boldsymbol{E}_{i,{\rm pre}}\psi\right)
+ o_p\left((N_{\rm tr} T_{\rm post})^{-1/2}\right),
\end{equation}
and consequently
\begin{equation}
\label{eq:CLT1}
\left( \hat{\tau}^{\rm sdid} - \tau\right) \,\big/\, {V_{\tau}^{1/2}} \, \Rightarrow \, {\cal N}\left(0, \, 1\right), \ \ {\rm where}\ \
V_{\tau} = \frac{1}{N_{\rm tr}} \operatorname{Var}_{}\left[\frac{1}{T_{\rm post}}\sum_{t=T_{\rm pre}+1}^T \varepsilon_{it} - \boldsymbol{E}_{i,{\rm pre}} \psi\right].
\end{equation}
Here $V_\tau$ is on the order of $1/(N_{{\rm tr}}T_{{\rm post}})$, i.e., $N_{\rm tr} T_{\rm post} V_{\tau}$ is bounded and bounded away from zero.
\end{theorem}
\section{Large-Sample Inference}
\label{sec:inference}
The asymptotic result from the previous section can be used to motivate practical methods for large-sample inference using SDID. Under appropriate conditions, the estimator is asymptotically normal and zero-centered; thus, if these conditions hold and
we have a consistent estimator for its asymptotic variance $V_\tau$, we can use conventional confidence intervals
\begin{equation}
\label{eq:GCI}
\tau \in \hat{\tau}^{\rm sdid} \pm z_{\alpha/2}\sqrt{\widehat{V}_\tau}
\end{equation}
to conduct asymptotically valid inference. In this section, we discuss three approaches to variance estimation for use in confidence intervals
of this type.
\begin{algorithm}[t]
\KwData{$\boldsymbol{Y},\boldsymbol{W}, B$}
\KwResult{ Variance estimator $\widehat{V}^{cb}_{\tau}$ }
\For{$i \leftarrow 1$ \KwTo $B$}{
Construct a bootstrap dataset $(\boldsymbol{Y}^{(b)}, \boldsymbol{W}^{(b)})$ by sampling $N$ rows of \\ $(\boldsymbol{Y}, \boldsymbol{W})$ with replacement. \\
\If{the bootstrap sample has no treated units or no control units}{ Discard and resample (\textbf{go to 2})}
Compute the SDID estimator $\hat \tau^{(b)}$ based on $(\boldsymbol{Y}^{(b)}, \boldsymbol{W}^{(b)})$
}
Define $\widehat{V}^{b}_\tau=\frac{1}{B}\sum_{b=1}^{B} (\hat{\tau}^{(b)} - \frac{1}{B}\sum_{b=1}^B\hat{\tau}^{(b)})^2$\;
\caption{Bootstrap Variance Estimation}
\label{alg:boot}
\end{algorithm}
The first proposal we consider, described in detail in Algorithm \ref{alg:boot}, involves a clustered bootstrap
\citep{efron1979bootstrap} where we independently resample units. As argued in \citet*{Bertrand2004did},
unit-level bootstrapping presents a natural approach to inference with panel data when repeated observations
of the same unit may be correlated with each other. The bootstrap is simple to implement and, in our experiments,
appears to yield robust performance in large panels.
The main downside of the bootstrap is that it may be computationally costly as it involves running
the full SDID algorithm for each bootstrap replication, and for large datasets this can be prohibitively expensive.
\begin{algorithm}[t]
\KwData{$\hat\omega, \hat \lambda, \boldsymbol{Y}, \boldsymbol{W}, \hat \tau$}
\KwResult{ Variance estimator $\widehat{V}_{\tau}$ }
\For{$i \leftarrow 1$ \KwTo $N$}{
Compute
$\hat \tau^{(-i)}:\operatorname*{arg\,min}_{\tau, \{\alpha_j,\beta_t\}_{j \ne i, t}}\sum_{j\ne i,t}\left(\boldsymbol{Y}_{jt} - \alpha_j - \beta_t - \tau \boldsymbol{W}_{it}\right)^2\hat\omega_j\hat \lambda_t$
}
Compute $\widehat{V}^{\mathrm{jack}}_\tau=(N-1)N^{-1}\sum_{i=1}^{N} (\hat{\tau}^{(-i)} - \hat{\tau})^2$\;
\caption{Jackknife Variance Estimation}
\label{alg:jack}
\end{algorithm}
To address this issue we next consider an approach to inference that is more closely tailored to the SDID method and only involves
running the full SDID algorithm once, thus dramatically decreasing the computational burden.
Given weights $\hat{\omega}$ and $\hat{\lambda}$ used to get the SDID point estimate, Algorithm \ref{alg:jack} applies the jackknife
\citep{miller1974jackknife} to the weighted SDID regression \eqref{main_sdid}, with the weights treated as fixed.
The validity of this procedure is not implied directly by asymptotic linearity as in \eqref{eq:asymptotic-linearity}; however, as shown below,
we still recover conservative confidence intervals under considerable generality.
\begin{theorem}
\label{theo:jack}
Suppose that the elements of $\boldsymbol{L}$ are bounded.
Then, under the conditions of Theorem \ref{theo:asymptotic-linearity},
the jackknife variance estimator described in Algorithm \ref{alg:jack} yields conservative confidence intervals,
i.e., for any $0 < \alpha < 1$,
\begin{equation}
\label{eq:jack_cons}
\liminf \mathbb{P}_{}\left[\tau \in \hat{\tau}^{\rm sdid} \pm z_{\alpha/2}\sqrt{\widehat{V}^{\mathrm{jack}}_\tau}\right] \geq 1 - \alpha.
\end{equation}
Moreover, if the treatment effects $\boldsymbol{\tau}_{it} = \tau$ are constant\footnote{When
treatment effects are heterogeneous, the jackknife implicitly treats the estimand \eqref{eq:estimand} as
random whereas we treat it as fixed, thus resulting in excess estimated variance; see \citet{imbens2004}
for further discussion.}
and
\begin{equation}
\label{eq:regr_cons}
T_{{\rm post}} N_{\rm tr}^{-1} \left\lVert\hat{\lambda}_0 + \boldsymbol{L}_{{\rm tr},{\rm pre}} \hat{\lambda}_{{\rm pre}} - \boldsymbol{L}_{{\rm tr},{\rm post}}\hat{\lambda}_{{\rm post}}\right\rVert_2^2 \rightarrow_p 0,
\end{equation}
i.e., the time weights $\hat{\lambda}$ are predictive enough on the exposed units,
then the jackknife yields exact confidence intervals and \eqref{eq:jack_cons} holds with equality.
\end{theorem}
In other words, we find that the jackknife is in general conservative and is exact
when treated and control units are similar enough that time weights that fit the control units generalize
to the treated units.
This result depends on specific structure of the SDID estimator, and does not hold for
related methods such as the SC estimator. In particular, an analogue to Algorithm \ref{alg:jack} for SC would be
severely biased upwards, and would not be exact even in the well-specified fixed effects model. Thus, we do
not recommend (or report results for) this type of jackknifing with the SC estimator. We do report results for
jackknifing DID since, in this case, there are no random weights $\hat{\omega}$ or $\hat{\lambda}$ and so our jackknife
just amounts to the regular jackknife.
Now, both the bootstrap and jackknife-based methods discussed so far are designed with the setting of Theorem \ref{theo:asymptotic-linearity} in
mind, i.e., for large panels with many treated units. These methods may be less reliable when the number of treated units $N_{\rm tr}$
is small, and the jackknife is not even defined when $N_{\rm tr} = 1$. However, many applications of synthetic controls have $N_{\rm tr} = 1$, e.g.,
the California smoking application from Section \ref{sec:calif}. To this end, we consider a third variance estimator that is motivated by placebo evaluations as often considered in the literature
on synthetic controls \citep*{Abadie2010,abadie2014}, and that can be applied with $N_{\rm tr} =1$. The main idea of such placebo evaluations
is to consider the behavior of synthetic control estimation when we replace the unit that was exposed to the treatment with
different units that were not exposed.\footnote{Such a placebo test is closely connected to permutation tests in randomization inference;
however, in many synthetic controls applications, the exposed unit was not chosen at random, in which case placebo tests do not have
the formal properties of randomization tests \citep{firpo2018synthetic,hahn2016}, and so may need to be interpreted via a
more qualitative lens.}
Algorithm \ref{alg:placebo} builds on this idea, and uses placebo predictions using only the unexposed units to estimate
the noise level, and then uses it to get \smash{$\widehat{V}_\tau$} and build confidence intervals as in \eqref{eq:GCI}. See \citet{bottmer2021design} for a discussion of the properties of such placebo variance estimators in small samples.
\begin{algorithm}[t]
\KwData{$\boldsymbol{Y}_{{\rm co},\cdot}, N_{{\rm tr}}, B$}
\KwResult{ Variance estimator $\widehat{V}^{\mathrm{placebo}}_{\tau}$ }
\For{$b \leftarrow 1$ \KwTo $B$}{
Sample $N_{{\rm tr}}$ out of the $N_{{\rm co}}$ control units without replacement to `receive the placebo'\;
Construct a placebo treatment matrix $\boldsymbol{W}_{{\rm co},\cdot}^{(b)}$ for the controls\;
Compute the SDID estimator $\hat \tau^{(b)}$\
based on $(\boldsymbol{Y}_{{\rm co},\cdot}, \boldsymbol{W}_{{\rm co},\cdot}^{(b)})$
;
}
Define $\widehat{V}^{\mathrm{placebo}}_\tau=\frac{1}{B}\sum_{b=1}^{B} (\hat{\tau}^{(b)} - \frac{1}{B}\sum_{b=1}^B\hat{\tau}^{(b)})^2$\;
\caption{Placebo Variance Estimation}
\label{alg:placebo}
\end{algorithm}
Validity of the placebo approach relies fundamentally on homoskedasticity across units, because if the exposed and unexposed units
have different noise distributions then there is no way we can learn \smash{$V_\tau$} from unexposed units alone. We also
note that non-parametric variance estimation for treatment effect estimators is in general impossible if we only have one
treated unit, and so homoskedasticity across units is effectively a necessary assumption in order for inference to be possible
here.\footnote{In Theorem \ref{theo:asymptotic-linearity}, we also assumed homoskedasticity. In contrast to the case of placebo
inference, however, it's likely that a similar result would also hold without homoskedasticity; homoskedasticity
is used in the proof essentially only to simplify notation and allow the use of concentration inequalities which have been
proven in the homoskedastic case but can be generalized.}
Algorithm \ref{alg:placebo} can also be seen as an adaptation of the method of \citet{conley2011inference} for inference in DID models with
few treated units and assuming homoskedasticity, in that both rely on the empirical distribution of residuals for placebo-estimators
run on control units to conduct inference. We refer to \citet{conley2011inference} for a detailed analysis of this class of algorithms.
\begin{table}[t]
\begin{center}
\begin{tabular}{|l|rrr|rrr|rrr|}
\hline
&\multicolumn{3}{c|}{Bootstrap} & \multicolumn{3}{c|}{Jackknife}& \multicolumn{3}{c|}{Placebo}\\
& SDID & SC & DID & SDID & SC & DID& SDID & SC & DID \\ \hline
Baseline & 0.96 & 0.93 & 0.89 & 0.93 & --- & 0.92 & 0.95 & 0.89 & 0.96 \\
\hline
Gun Law & 0.97 & 0.96 & 0.93 & 0.94 & --- & 0.93 & 0.94 & 0.95 & 0.93 \\
Abortion & 0.96 & 0.94 & 0.93 & 0.93 & --- & 0.95 & 0.97 & 0.91 & 0.96 \\
Random & 0.96 & 0.96 & 0.92 & 0.93 & --- & 0.94 & 0.96 & 0.96 & 0.94 \\
Hours & 0.92 & 0.96 & 0.94 & 0.89 & --- & 0.95 & 0.91 & 0.89 & 0.96 \\
Urate & 0.91 & 0.90 & 0.57 & 0.86 & --- & 0.64 & 0.88 & 0.89 & 0.62 \\
\hline
$T_{post}= 1$ &0.93 & 0.94 & 0.84 & 0.92 & --- & 0.88 & 0.93 & 0.91 & 0.92\\
$N_{tr} = 1$ & --- & --- & --- & --- & --- & --- & 0.97 & 0.95 & 0.96 \\
$T_{post}= N_{tr} = 1 $ & --- & --- & --- & --- & --- & --- & 0.95 & 0.94 & 0.94\\
\hline
Resample, $N=200$ & 0.94 & 0.95 & 0.92 & 0.95 & --- & 0.93 & 0.96 & 0.94 & 0.94 \\
Resample, $N=400$ & 0.95 & 0.92 & 0.96 & 0.96 & --- & 0.95 & 0.96 & 0.91 & 0.96\\
\hline
Democracy & 0.93 & 0.96 & 0.55 & 0.94 & --- & 0.59 & 0.98 & 0.97 & 0.79 \\
Education &0.95 & 0.95 & 0.30 & 0.95 & --- & 0.34 & 0.99 & 0.90 & 0.94 \\
Random & 0.93 & 0.95 & 0.89 & 0.96 & --- & 0.91 & 0.95 & 0.94 & 0.91 \\
\hline
\end{tabular}
\caption{Coverage results for nominal 95\% confidence intervals in the CPS and Penn World Table simulation setting from Tables \ref{table1} and \ref{table2}.
The first three columns show coverage of confidence intervals obtained via the Placebo method.
The second set of columns show coverage from the jackknife method.
The last set of columns show coverage from the clustered bootstrap.
Unless otherwise specified, all settings have $N = 50$ and $T = 40$
cells, of which at most $N_{\rm tr} = 10$ units and $T_{\rm post} = 10$ periods are treated. In rows 7-9, we reduce the number
of treated cells. In rows 10 and 11, we artificially make the panel larger by adding rows, which makes the assumption that the number of treated units is small relative to the number of control units more accurate (we set $N_{\rm tr}$ to $10\%$ of the total number of units).
We do not report jackknife and bootstrap coverage rates for $N_{\rm tr} = 1$ because the estimators are not well-defined. We do not report jackknife coverage rates for SC because, as discussed in the text, the variance estimator is not well justified in this case.
All results are based on 400 simulation replications.}
\label{table_cps_cov}
\end{center}
\end{table}
Table \ref{table_cps_cov} shows the coverage rates for the experiments described in Section \ref{sec:cps} and \ref{sec:penn},
using Gaussian confidence intervals \eqref{eq:GCI} with variance estimates obtained as described above.
In the case of the SDID estimation, the bootstrap estimator performs particularly well, yielding nearly nominal $95\%$ coverage,
while both placebo and jackknife variance estimates also deliver results that are close to the nominal $95\%$ level.
This is encouraging, and aligned with our previous observation that the SDID estimator
appeared to have low bias. That being said, when assessing the performance of the placebo estimator, recall
that the data in Section \ref{sec:cps} was generated with noise that is both Gaussian and homoskedastic across
units---which were assumptions that are both heavily used by the placebo estimator.
In contrast, we see that coverage rates for DID and SC can be relatively low, especially in cases with significant bias
such as the setting with the state unemployment rate as the outcome. This is again in line with what one may have expected
based on the distribution of the errors of each estimator as discussed in Section \ref{sec:cps}, e.g., in
Figure \ref{fig_mw_rand}: If the point estimates $\hat{\tau}$ from DID and SC are dominated by bias, then we should not
expect confidence intervals that only focus on variance to achieve coverage.
\section{Related Work}
\label{sec:relworks}
Methodologically, our work draws most directly from the literature on SC methods, including
\citet{abadie2003}, \citet*{Abadie2010, abadie2014}, \citet{abadie2016}, \citet{doudchenko2016balancing}, and \citet*{ben2018augmented}.
Most methods in this line of work can be thought of as focusing on constructing unit weights that create comparable
(balanced) treated and control units, without relying on any modeling or weighting across time.
\citet*{ben2018augmented} is an interesting exception. Their augmented synthetic control estimator,
motivated by the augmented inverse-propensity weighted estimator of \citet*{robins1994estimation},
combines synthetic control weights with a regression adjustment for improved accuracy (see also \citet*{kellogg2020combining} which explicitly connects SC to matching).
They focus on the case of $N_{\rm tr} = 1$ exposed units and $T_{\rm post} = 1$ post-exposure periods, and their method involves fitting a
model for the conditional expectation $m(\cdot)$ for $Y_{iT}$ in terms of the lagged outcomes $\boldsymbol{Y}_{i,{\rm pre}}$, and
then using this fitted model to ``augment'' the basic synthetic control estimator as follows.
\begin{equation}
\label{eq:asc}
\begin{split}
\hat{\tau}_{\rm asc}
&= Y_{NT} - \left(\sum_{i=1}^{N-1} \hat{\omega}^{\rm sc}_i Y_{iT}+ \left(
\hat m(\boldsymbol{Y}_{N,{\rm pre}})-\sum_{i=1}^{N-1}\hat{\omega}^{\rm sc}_i \hat m(\boldsymbol{Y}_{i,{\rm pre}})\right)\right).
\end{split}
\end{equation}
Despite their different motivations, the augmented synthetic control and synthetic difference in differences methods share an
interesting connection: with a linear model $m(\cdot)$,
$\hat{\tau}_{\rm sdid}$ and $\hat{\tau}_{\rm asc}$ are very similar. In fact, had we fit $\hat{\omega}^{{\rm sdid}}$ without intercept,
they would be equivalent for $\hat m(\cdot)$ fit by least squares on the controls,
imposing the constraint that its coefficients are nonnegative and to sum to one, that is,
for $\hat m(\boldsymbol{Y}_{i,{\rm pre}}) = \hat{\lambda}_0^{{\rm sdid}} + \boldsymbol{Y}_{i,{\rm pre}}\hat{\lambda}_{{\rm pre}}^{{\rm sdid}}$.
This connection suggests that weighted two-way bias-removal methods are a natural way of working with
panels where we want to move beyond simple difference in difference approaches.
We also note recent work of \citet{roth2018pre} and \citet{rambachan2019honest},
who focus on valid inference in difference in differences settings when users look at past outcomes to
check for parallel trends. Our approach uses past data not only to check
whether the trends are parallel, but also to construct the weights to make them parallel. In this setting, we show
that one can still conduct valid inference, as long as $N$ and $T$ are large enough and the size of the treatment block is small.
In terms of our formal results, our paper fits broadly in the literature on panel models with interactive fixed effects and the matrix completion literature \citep{athey2017matrix, bai2009panel, moon2015linear,moon2017dynamic, robins1985comparison,xu2017generalized}.
Different types of problems of this form have a long tradition in the econometrics literature,
with early results going back to \citet*{ahn2001gmm}, \citet{chamberlain1992efficiency} and \citet*{holtz1988estimating} in the case
of finite-horizon panels (i.e., in our notation, under asymptotics where $T$ is fixed and only $N \rightarrow \infty$).
More recently, \citet{freyberger2018non} extended the work of \citet{chamberlain1992efficiency} to a setting that's
closely related to ours, and emphasized the role of the past outcomes for constructing moment restrictions in the fixed-$T$ setting.
\citet{freyberger2018non} attains identification by assuming that the errors $\boldsymbol{E}_{it}$ are uncorrelated,
and thus past outcomes act as valid instruments. In contrast, we allow for correlated errors within rows, and thus need to work in a large-$T$ setting.
Recently, there has considerable interest in models of type \eqref{basic_model} under asymptotics where both $N$ and $T$
get large. One popular approach, studied by \citet{bai2009panel} and \citet{moon2015linear,moon2017dynamic}, involves
fitting \eqref{basic_model} by ``least squares'', i.e., by minimizing squared-error loss while constraining $\widehat{\boldsymbol{L}}$ to
have bounded rank $R$. While these results do allow valid inference for $\tau$, they require strong assumptions. First,
they require the rank of $\boldsymbol{L}$ to be known a-priori (or, in the case of \citet{moon2015linear}, require a known upper bound
for its rank), and second, they require a $\beta_{\min}$-type condition whereby the normalized non-zero singular values of $\boldsymbol{L}$
are well separated from zero. In contrast, our results require no explicit limit on the rank of $\boldsymbol{L}$ and allow for $\boldsymbol{L}$ to have
to have positive singular values that are arbitrarily close to zero, thus suggesting that the SDID method may be more robust than the
least squares method in cases where the analyst wishes to be as agnostic as possible regarding properties of $\boldsymbol{L}$.\footnote{By
analogy, we also note that, in the literature on high-dimensional inference, methods that do no assume a uniform lower bound on
the strength of non-zero coefficients of the signal vector are generally considered more robust than ones that do
\citep*[e.g.,][]{belloni2014inference,zhang2014confidence}.}
\citet*{athey2017matrix}, \citet*{amjad2018robust}, \citet{moon2018nuclear} and \citet{xu2017generalized} build on this line of work, and replace the fixed-rank constraint with data-driven regularization on $\widehat{\boldsymbol{L}}$. This innovation is very helpful from a computational perspective; however, results for inference about
$\tau$ that go beyond what was available for least squares estimators are currently not available. We also note recent papers that
draw from these ideas in connection to synthetic control type analyses, including \citet{chan2020pcdid} and \citet{gobillon2016regional}.
Finally, in a paper contemporaneous to ours, \citet*{agarwal2019robustness} provide improved bounds from principal component
regression in an errors-in-variables model closely related to our setting, and discuss implications for estimation in synthetic
control type problems. Relative to our results, however, \citet{agarwal2019robustness} still require assumptions on the behavior
of the small singular values of $\boldsymbol{L}$, and do not provide methods for inference about $\tau$.
In another direction, several authors have recently proposed various methods that implicitly control for the systematic component
$\boldsymbol{L}$ in models of time \eqref{basic_model}. In one early example, \citet*{hsiao2012panel} start with a factor model similar to ours and
show that under certain assumptions it implies the moment condition
\begin{equation}
\label{eq:linear_model_stat}
Y_{Nt}= a+ \sum_{j=1}^{N-1}\beta_j Y_{jt} + \epsilon_{Nt}, \ \ \ \ \ \mathbb{E}_{}\left[\varepsilon_{Nt} \,\big|\, \{Y_{jt}\}_{j=1}^{N-1}\right]=0,
\end{equation}
for all $t = 1, \, \ldots, \, T$. The authors then estimate $\beta_j$ by (weighted) OLS. This approach is further refined by \cite{li2017estimation},
who additionally propose to penalizing the coefficients $\beta_j$ using the lasso \citep{tibshirani1996regression}.
In a recent paper, \citet*{chernozhukov2018inference} use the model \eqref{eq:linear_model_stat} as a starting point for inference.
While this line of work shares a conceptual connection with us, the formal setting is very different. In order to derive a representation of
the type \eqref{eq:linear_model_stat}, one essentially needs to assume a random specification for \eqref{basic_model} where both $\boldsymbol{L}$ and
$\boldsymbol{E}$ are stationary in time. \citet{li2017estimation} explicitly assumes that the outcomes $\boldsymbol{Y}$ themselves are weakly stationary, while
\citet*{chernozhukov2018inference} makes the same assumption to derive the results that are valid under general misspecification.
In our results, we do not assume stationarity anywhere: $\boldsymbol{L}$ is taken as deterministic and the errors $\boldsymbol{E}$ may be non-stationary.
Moreover, in the case of most synthetic control and difference in differences analyses, we believe stationarity to be a fairly restrictive
assumption. In particular, in our model, stationarity would imply that a simple pre-post comparison for exposed units would be an
unbiased estimator of $\tau$ and, as a result, the only purpose of the unexposed units would be to help improve efficiency.
In contrast, in our analysis, using unexposed units for double-differencing is crucial for identification.
\citet{ferman2019synthetic} analyze the performance of synthetic control estimator using essentially the same model as we do. They focus on the situations where $N$ is small, while $T_{\rm pre}$ (the number of control periods) is growing. They show that unless time factors have strong trends ({\it e.g.}, polynomial) the synthetic control estimator is asymptotically biased. Importantly \cite{ferman2019synthetic} focus on the standard synthetic control estimator, without time weights and regularization, but with an intercept in the construction of the weights.
Finally, from a statistical perspective, our approach bears some similarity to the work on ``balancing'' methods for
program evaluation under unconfoundedness, including \citet*{athey2018approximate}, \citet*{graham1},
\citet{hirshberg2017augmented}, \citet{imai2014covariate}, \citet{kallus2020generalized}, \citet{zhao2019covariate}
and \citet{zubizarreta2015stable}. One major result of this line of work is that, by algorithmically finding weights that balance
observed covariates across treated and control observations, we can derive robust estimators with good asymptotic properties
(such as efficiency). In contrast to this line of work, rather than balancing observed covariates, we here need to balance unobserved
factors $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$ in \eqref{basic_model} to achieve consistency; and accounting for this forces us to follow a different formal
approach than existing studies using balancing methods.
\bibliographystyle{plainnat}
\bibliography{references}
\newpage