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.
94,203 characters
An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls
\if11
{
\title{\bf An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls}
\author{Victor Chernozhukov\thanks{Massachusetts Institute of Technology; 50 Memorial Drive, E52-361B, Cambridge, MA 02142, USA; Email: [email removed]} \qquad Kaspar W\"uthrich\thanks{Department of Economics, University of California San Diego, 9500 Gilman Dr., La Jolla, CA 92093, USA; Email: [email removed]} \qquad Yinchu Zhu\thanks{Brandeis University; 415 South Street, Waltham, MA 02453, USA; Email: [email removed]}}
} \fi
\if01
{
\title{\bf An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls}
\date{ }
} \fi
\maketitle
\begin{abstract}
We introduce new inference procedures for counterfactual and synthetic control methods for policy evaluation. We recast the causal inference problem as a counterfactual prediction and a structural breaks testing problem. This allows us to exploit insights from conformal prediction and structural breaks testing to develop permutation inference procedures that accommodate modern high-dimensional estimators, are valid under weak and easy-to-verify conditions, and are provably robust against misspecification. Our methods work in conjunction with many different approaches for predicting counterfactual mean outcomes in the absence of the policy intervention. Examples include synthetic controls, difference-in-differences, factor and matrix completion models, and (fused) time series panel data models. Our approach demonstrates an excellent small-sample performance in simulations and is taken to a data application where we re-evaluate the consequences of decriminalizing indoor prostitution. Open-source software for implementing our conformal inference methods is available.
\bigskip
\noindent \textit{Keywords:} permutation inference, model-free validity, difference-in-differences, factor model, matrix completion, constrained Lasso
\vfill
\end{abstract}
\newpage
\spacingset{1.5}
\pagenumbering{arabic}
\section{Introduction}
We consider the problem of making inferences on the causal effect of a policy intervention in an aggregate time series setup with a single treated unit. The treated unit is observed for $T_0$ periods before and $T_\ast$ periods after the intervention occurs. Often, there is additional information in the form of possibly very many untreated units, which can serve as controls. Such settings are ubiquitous in applied research, and there are many different approaches for estimating the causal effect of the policy. Examples include difference-in-differences methods, synthetic control (SC) approaches, factor, matrix completion, and interactive fixed effects (FE) models, and times series models.\footnote{We refer to \citet{DI16}, \citet{gobillon2016regional}, and \citet{abadie2019jel} for excellent comparative overviews and reviews.}
We refer to these methods as counterfactual and synthetic control (CSC) methods.
This paper provides generic and robust procedures for making inferences on policy effects estimated by CSC methods. We propose a general counterfactual modeling framework that nests and generalizes many traditional and new methods for counterfactual analysis. We focus on methods that are able to generate mean-unbiased proxies, $P_t^N$, for the counterfactual outcomes of the treated unit in the absence of the policy intervention, $Y_{1t}^N$:
\[
Y_{1t}^N=P_t^N+u_t, \quad E\left(u_t \right)=0, \quad t=1,\dots,T_0+T_\ast.
\]
The policy effect in period $t$ is $\theta_t=Y_{1t}^I-Y_{1t}^N$, where $Y_{1t}^I$ is the counterfactual outcome of the treated unit with the policy intervention. We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, $\theta=\{\theta_t\}_{t=T_0+1}^{T_0+T_\ast}$. Specifically, we postulate a trajectory $\theta^0= \{ \theta^0_t\}_{t=T_0+1}^{T_0+T_\ast}$ and test the sharp null hypothesis that $\theta=\theta^0$. We also consider the problem of testing hypotheses about per-period effects $\theta_t$ and propose a simple algorithm for constructing pointwise confidence intervals via test inversion.
We recast the inference problem as a (counterfactual) prediction and a structural breaks testing problem. This allows us to build on the literature on conformal prediction \citep{vovk2005algorithmic} and end-of-sample stability testing \citep[][]{dufour1994generalized,andrews2003end} to construct inference procedures that are provably robust against misspecification and accommodate many classical and modern high-dimensional methods for estimating $P_t^N$.
The basic idea of our testing procedures is as follows. Suppose that there is only one post-treatment period and that $P_t^N$ is known. Under the sharp null that $\theta_{T_0+1}=\theta_{T_0+1}^0$, we can compute $Y_{1t}^N$ and $u_t=Y_{1t}^N-P_t^N$ for all time periods. If the stochastic shock process $\{u_t\}$ is stationary and weakly dependent, and its distribution is invariant under the intervention, the distribution of the error in the post-treatment period, $u_{T_0+1}$, should be the same as the distribution of the errors in the pre-treatment period, $\{u_t\}_{t=1}^{T_0}$. We operationalize this idea by proposing inference methods in which $p$-values are obtained by permuting blocks of estimated residuals across the time series dimension.
The proposed methods are valid under two different sets of conditions:
\begin{itemize} \setlength{\parskip}{0pt}
\item[\textbf{(i)}] \textbf{Estimator Consistency and Stationary Weakly Dependent Errors}
If the data exhibit dynamics, trends, and serial dependence but the stochastic shock sequence $\{u_t\}$ is stationary and weakly dependent, our inference procedures are approximately valid if the estimator of $P_t^N$ is consistent (pointwise and in prediction norm). Consistency can be verified for many different CSC methods. We provide concrete sufficient conditions for a representative selection of methods, including difference-in-differences, SC, factor models, matrix completion and interactive FE models, linear and nonlinear time series models, and fused time series panel data models.
\item[\textbf{(ii)}] \textbf{Estimator Stability and Stationary Weakly Dependent Data}
In practice, misspecification is an important concern. We show that even if the model for $P_t^N$ is misspecified and the estimator of $P_t^N$, $\hat{P}_t^N$, is inconsistent, our procedures are still valid, provided that the data are stationary and weakly dependent and $\hat{P}_t^N$ satisfies a \emph{stability} condition. This condition requires that $\hat{P}_t^N$ is stable under perturbations in a few observations. It is implied, for instance, if $\hat{P}_t^N$ is consistent for a pseudo-true parameter value but is shown to hold even in high-dimensional settings where consistency results under misspecification are often not available.
\end{itemize}
The main theoretical results in this paper are finite sample (non-asymptotic) bounds on the size accuracy of our methods; these bounds imply that our methods are exact as $T_0\rightarrow \infty$. Unlike traditional asymptotic results, which are only informative when the sample size is large enough, our non-asymptotic bounds show how different factors affect the finite sample performance. This feature is relevant in CSC applications where sample sizes are often small.
A key feature of our conformal inference methods is that $P_t^N$ is estimated under the null hypothesis based on data from all $T_0+T_\ast$ periods. Estimation under the null guarantees the exact finite sample validity of our procedures if the data are iid or exchangeable. Even when exchangeability fails, imposing the null for estimation is essential for a good performance in CSC applications where $T_0$ is often rather small. Figure \ref{fig:importance_H0} plots the empirical rejection probabilities for testing the null that $\theta_{T_0+1}=0$ when $T_0=19$, $J=50$ (as in our empirical application), $\{u_t\}$ is an AR(1) process, and $P_t^N$ is estimated using SC. The size properties of our method are excellent. By contrast, estimating $P_t^N$ based on the $T_0$ pre-treatment periods without imposing the null yields substantial size distortions. Figure \ref{fig:importance_H0} suggests that imposing the null continues to improve size accuracy even when exchangeability fails and that these improvements can be substantial in small samples.
\begin{center}
[Figure \ref{fig:importance_H0} around here.]
\end{center}
We make two additional contributions that may be of independent interest. First, we introduce the $\ell_1$-constrained least squares estimator or \emph{constrained Lasso} \citep[e.g.,][]{raskutti2011minimax} as an essentially tuning-free alternative to existing penalized regression estimators and study its theoretical properties. Constrained Lasso nests SC and difference-in-differences, providing a unifying approach for the regression-based estimation of the mean proxies $P_t^N$. Second, we obtain theoretical consistency results for SC estimators in settings with potentially very many control units.
We develop three extensions of our main results. First, we show that our method can be modified to test hypotheses about average effects over time. Second, we extend our method to settings with multiple treated units. Third, we propose easy-to-implement placebo tests for assessing the credibility of inferences based on our method.
Monte Carlo simulations suggest that our procedures exhibit excellent size properties and are robust to misspecification. We find that imposing additional constraints (e.g., using SC instead of the more general constrained Lasso) does not improve power when these additional restrictions are correct but can cause power losses when they are not.
Finally, we re-analyze the causal effect of decriminalizing indoor prostitution on sexually transmitted infections. Following \citet{cunningham2018decriminalizing}, we exploit the unanticipated decriminalization of indoor prostitution in Rhode Island in 2003. We find that decriminalizing indoor prostitution significantly decreased the incidence of female gonorrhea.
\paragraph{Related Literature.} We contribute to the literature on inference procedures for CSC methods with few treated units.
A popular method is the finite population permutation approach of \citet{abadie10sc}, see also \citet{firpo18synthetic} and \citet{abadie2019jel}. This approach permutes the policy assignment and relies on permutation distributions for inference. It corresponds to conventional randomization inference \citep{fisher1935design} under random assignment of the policy \citep[e.g.,][]{abadie10sc,abadie2019jel}. However, random assignment is not plausible in typical CSC applications, and assignment mechanisms are difficult to model and estimate when there are only few treated units \citep[e.g.,][]{abadie2019jel}. \citet{shaikh2019randomization} propose randomization tests for settings with staggered treatment adoption, which encompass the approach of \citet{abadie10sc}. The main assumption of \citet{shaikh2019randomization}'s approach is that policy adoption follows a Cox proportional hazards model. We do not model the assignment mechanism. Instead, we exploit stationarity and weak dependence of the errors across time in a repeated sampling framework. One advantage of exploiting the time series dimension is that we only require a suitable model for the potential outcome of the treated unit. By contrast, cross-sectional approaches often require estimating models for all units. That is, we only require a good ``local'' instead of a good ``global'' fit, which reduces the risk of model misspecification. On the other hand, our approach requires a large number of pre-treatment periods and relies on invariance of the error distribution under the intervention.
There is also an active literature on asymptotic inference methods for CSC models. Several papers focus on testing hypotheses about average or expected effects over time, requiring $T_0$ and $T_\ast$ to be large. \citet{li2017estimation}, \citet{carvalho2018arco}, \citet{chernozhukov2019ttest}, and \citet{li2020statistical} introduce inference methods based on penalized and constrained regression methods. \citet{arkhangelsky2018synthetic} propose inference methods for a version of SC with time and unit weights, which admits a weighted regression formulation. Asymptotic inference methods based on factor and interactive FE models are proposed by \citet{hsiao2012panel}, \citet{gobillon2016regional}, \citet{chan2016policy}, \citet{li2017estimation}, \citet{xu2017generalized}, and \citet{li2018inference}.
Here we focus on sharp null hypotheses and permutation distributions and provide non-asymptotic performance guarantees. Our approach is generic and valid with many different methods, including constrained regression, factor models, and interactive FE estimators. \citet{conley2011inference} propose inference methods for difference-in-difference settings with few treated units. They exploit the cross-sectional dimension, relying on weak dependence and stationarity of the error terms across units, which may be hard to justify in typical CSC settings. By contrast, our procedures rely on stationarity and weak dependence of the errors over time. On the other hand, exploiting the time series dimension, our approach requires $T_0$ to be large, whereas \citet{conley2011inference} allow $T_0$ to be fixed.
In related work, \citet{cattaneo2021prediction} provide prediction intervals for per-period effects estimated by SC methods. Their key observation is that there is both randomness from estimating the SC weights and from the prediction error. They propose a sampling-based inference method based on non-asymptotic probability bounds that accounts for both types of randomness and is valid with stationary and non-stationary data.
We recast the causal inference problem as a (counterfactual) prediction problem and build on the literature on conformal prediction \citep[e.g.,][]{vovk2005algorithmic,vovk2009online,lei2013distribution,lei2014distribution,lei2017distributionfree} and on the literature on permutation tests \citep[e.g.,][]{romano1990behavior,lehmann2005testing}, which was started by \citet{fisher1935design} in the context of randomization; see \citet{rubin1984bayesianly} for a Bayesian justification. On a more general level, our approach is also connected to transformation-based approaches to model-free prediction \citep[e.g.,][]{politis2015modelfree}. Let us discuss in more detail the relationship to \citet[][CWZ18 henceforth]{chernozhukov2018exact}, who extend classical conformal prediction to time series settings. Besides a different focus (prediction intervals for future outcome values vs.\ inference on policy effects), there are several important differences. First, we rely on permuting residuals, whereas CWZ18 permute blocks of data. Second, we theoretically analyze different types of permutations. In particular, we study the set of all permutations, which yields precise $p$-values in small samples. This set of permutations cannot be used in the framework of CWZ18 unless the data are iid. Third, we allow for non-stationary data, whereas the prediction methods in CWZ18 are strictly limited to stationary data. Forth, CWZ18 rely on abstract high-level conditions on the test statistics and do not provide any primitive conditions. By contrast, we develop transparent sufficient conditions that can be verified for many traditional and modern CSC methods, and we provide explicit primitive conditions for a large selection of popular approaches. Finally, we establish the validity of our methods with time series data under misspecification and stability, whereas the theoretical results for weakly dependent data in CWZ18 require correct specification and consistency.
Finally, we show that the problem of making inferences on policy effects can be recast as a structural breaks testing problem with a known break date. Therefore, we build on and contribute to the literature on testing for structural breaks and, in particular, to the literature on structural breaks testing using permutation approaches \citep[e.g.,][]{antoch2001permutation,zeileis2013toolbox}. Besides a different focus (inference on policy effects vs.\ testing for structural breaks), our paper differs from the existing literature in that we specifically focus on testing at the end of the sample, allow for a very general class of estimators, including modern high-dimensional methods, and provide non-asymptotic performance guarantees under correct specification and misspecification. Our paper is also related to tests for structural breaks at the end of the sample \citep[e.g.,][]{dufour1994generalized,andrews2003end}.\footnote{\citet{hahn2017synthetic} informally suggest applying a variant of \citet{andrews2003end}'s end-of-sample stability test in the context of SC, and \citet{ferman2019inference} use a version of this test in the context of difference-in-differences approaches with few treated groups.} Let us discuss the differences to \citet{andrews2003end}'s end-of-sample instability test based on subsampling in more detail. First, we focus on causal inference, whereas \citet{andrews2003end} is concerned with structural breaks testing. Second, our procedures are exactly valid under exchangeability, and we obtain finite sample bounds under weak conditions on the estimators, while the theoretical properties of \citet{andrews2003end}'s test rely on asymptotic analyses. Third, our methods are valid under misspecification, whereas \citet{andrews2003end} assumes correct specification. Forth, our results under correct specification only require stationarity and weak dependence of $\{u_t\}$, while \citet{andrews2003end}'s test assumes stationarity of the data.\footnote{ \citet{andrews2003end} briefly comments on page 1681 (comment 4) that his test can be shown to be asymptotically valid under stationary errors but does not provide a formal result.} Finally, our procedures work in conjunction with many modern high-dimensional estimators, whereas \citet{andrews2003end} focuses on low-dimensional GMM models.
\paragraph{Notation.} For $q\geq 1$, the $\ell_q$-norm of a vector is denoted by $\| \cdot\|_{q}$. We use $\|\cdot\|_0$ to denote the number of nonzero entries of a vector; $\|\cdot\|_{\infty}$ is used to denote the maximal absolute value of entries of a vector. We use the
notation $a\lesssim b$ to denote $a \leq cb$ for some constant $c > 0$ that does not depend on the sample size. We use the notation $a \asymp b $ to denote $a\lesssim b$ and $b\lesssim a$. For a set $A$, $|A|$ denotes the cardinality of $A$. For any $a\in\mathbb{R} $, we define $\lfloor a \rfloor =\max\{z\in\mathbb{Z}:z\leq a\} $ and $\lceil a \rceil =\min\{z\in\mathbb{Z}:z\geq a\} $, where $\mathbb{Z}$ is the set of integers. We use $\mathbb{N}$ to denote the set of natural numbers.
\section{A Conformal Inference Method}
\label{sec:conformal_inference}
\subsection{The Counterfactual Model}
We consider a time series of $T$ outcomes for a treated unit, labeled $j=1$. During the first
$T_0$ periods, the unit is not treated by a policy and, during the remaining $T-T_0=T_\ast$ periods, it is treated
by the policy. Extensions to more than one treated unit are discussed in the Appendix. Our typical setting is where $T_\ast$ is short compared to $T_0$. There may be other units that are not exposed to the policy, and they will be introduced below. We denote the observed outcome of the treated unit by $Y_{1t}$. We employ the potential (latent) outcomes framework \citep{neyman1923application,rubin1974estimating} and denote potential outcomes with and without the policy as $Y_{1t}^I$ and $Y_{1t}^N$. The effect of the policy intervention in period $t$ is $\theta_t=Y_{1t}^I-Y_{1t}^N$.
Our conformal inference method will rely on the following counterfactual modeling framework, which nests many traditional and new methods for counterfactual policy analysis; see Sections \ref{subsec:sc_panel}--\ref{subsec:ts_fused} for examples.
\begin{assumption}[Counterfactual Model]\label{ass:dgp}
Let $\{P_t^N\}$ be a given sequence of mean-unbiased predictors or proxies for the counterfactual outcomes $\{Y_{1t}^N\}$ in the absence of the policy intervention, that is $\{E\left( P_t^N\right)\} = \{E \left( Y_{1t}^N\right)\}$. Let $\{\theta_t\}$ be a fixed policy effect sequence with $\theta_t = 0 $ for $ t \leq T_0$, so that potential outcomes under the intervention are given by $\{Y_{1t}^I\} = \{Y_{1t}^N + \theta_t\}$.\footnote{In the Appendix, we consider an extension to random policy effects.}
In other words, potential outcomes can be written as
\begin{equation}\tag{CMF}
\begin{array}{l}
Y_{1t}^N = P_t^N + u_t\\
Y_{1t}^I = P_t^N + \theta_t + u_t \\
\end{array} \Bigg | \quad E( u_t) = 0, \quad t=1,\dots,T , \\
\end{equation}
where $\{u_t\}$ is a centered stationary stochastic process. Observed outcomes are related to potential outcomes as $Y_{1t}=Y_{1t}^{N}+D_{t}\left(Y_{1t}^{I}-Y_{1t}^{N}\right)$, where $D_{t}=1\left(t>T_0 \right)$.
\end{assumption}
Assumption \ref{ass:dgp} introduces the potential outcomes, but also postulates an identifying assumption in the form of the existence of mean-unbiased proxies $P^N_t$ such that $E \left( P_t^N\right) = E \left( Y_{1t}^N \right)$. Assumption \ref{ass:dgp} allows $\{P_t^N\}$ to be fixed or random and does not impose any restrictions on the dependence between $\{P_t^N\}$ and $\{u_t\}$. In Sections \ref{subsec:sc_panel}--\ref{subsec:ts_fused}, we will discuss specific panel data and time series models that postulate (and identify) what $P_t^N$ is under a variety of conditions. Additional assumptions on the stochastic shock process $\{u_t\}$ will be introduced later, in essence requiring $\{u_t\}$ to be either iid or, more generally, a stationary and weakly dependent process.
Assumption \ref{ass:dgp} also postulates that the stochastic shock sequence is invariant under the intervention. This is the fundamental identifying assumption. It requires that the timing of the policy intervention is independent of factors that change the distribution of $\{u_t\}$.\footnote{In principle, we can relax this assumption by specifying, for example, the scale and quantile shifts in the stochastic shocks that result from the policy, and then working with the resulting model; we leave this extension to future work.} If the policy changes the distribution of $\{u_t\}$, one can either interpret our method as a structural breaks test or view the policy effect as random in which case our method yields valid prediction sets; see the Appendix for details.
Often, there is additional information in the form of untreated units, which can serve as controls. Specifically, suppose that there are $J \geq 1$ control units, indexed by $j=2,\dots,J+1$. We assume that we observe all units for all $T$ periods, although this assumption can be relaxed. Let $Y_{jt}$ denote the observed outcome for these untreated units. This observed outcome is equal to the outcome in the absence of the policy intervention, i.e., $Y_{jt} = Y_{jt}^N$ for $2\le j\le J+1$ and $1\le t\le T$. For each unit, we may also observe a vector of covariates $X_{jt}$. This motivates a variety of strategies for modeling and identifying $P_t^N$ as discussed below.
\subsection{Hypotheses of Interest, Test Statistics, and $p$-Values}
\label{subsec:hypotheses}
We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, $\theta=\left(\theta_{T_0+1},\dots,\theta_{T}\right)'$. Our main hypothesis of interest is
\begin{eqnarray}\label{eq:H0}
H_0:~\theta=\theta^0,
\end{eqnarray}
where $\theta^0=\left(\theta^0_{T_0+1},\dots,\theta^0_{T}\right)'$ is a postulated policy effect trajectory. Hypothesis \eqref{eq:H0} is a sharp null hypothesis. It fully determines the value of the counterfactual outcome in the absence of the intervention in the post-treatment period since $Y_{1t}^N=Y^I_{1t}-\theta_t=Y_{1t}-\theta_t$. In the Appendix, we show that our method can also be used to test hypotheses about average effects.
To describe our procedure, we write the data under the null hypothesis as $\mathbf{Z}:=\mathbf{Z}(\theta^0)=(Z_1,\dots,Z_{T})'$, where
\[
Z_{t}=\begin{cases}
\left(Y^N_{1t},Y^N_{2t},\dots,Y^N_{J+1t},X'_{1t},\dots,X'_{J+1t}\right)', & t\le T_0\\
\left(Y_{1t}^I-\theta_t^{0},Y^N_{2t},\dots,Y^N_{J+1t},X'_{1t},\dots,X'_{J+1t}\right)', & t>T_0.
\end{cases}
\]
Using one of the methods described below, we will obtain a counterfactual proxy estimate, $\hat{P}^N_t$,
based on $\mathbf{Z}$, and compute the residuals $\hat{u}=\left(\hat{u}_1,\dots,\hat{u}_T \right)'$, where $\hat{u}_t=Y_{1t}^N-\hat{P}^N_t$ for $1\le t\le T$.
Since $P_t^N $ is computed using $\mathbf{Z}=\mathbf{Z}(\theta^0)$, $P_t^N$ is estimated under the null hypothesis, which is essential for a good small sample performance. In ``ideal'' settings where the data are iid, imposing the null guarantees the model-free exact finite sample validity of our method; see the Appendix for details. By contrast, when $P_t^N$ is estimated based on the pre-treatment data $\{Z_t\}_{t=1}^{T_0}$ without imposing the null, permuting blocks of residuals does not yield procedures with exact finite sample validity, not even with iid data.
\smallskip
\noindent \textbf{Definition of Test Statistic $\mathbf{S}$.} We consider the following test statistic:
\[
S (\hat u) = S_q(\hat u) = \left ( \frac{1}{\sqrt{T_{\ast}}} \sum_{t=T_0+1}^T|\hat{u}_t|^q \right)^{1/q}.
\]
\smallskip
Note that $S$ is constructed such that high values indicate rejection. Different choices of $q$ lead to power against different alternatives. For instance, if the intervention has a large but only temporary effect (i.e., if $|\theta_t|$ is large for few periods), choosing $q=\infty$ yields high power. On the other hand, if the intervention has a permanent effect (i.e., if $\theta_t$ is non-zero for many post-treatment periods), tests using $S_1$ or $S_2$ exhibit good power properties. In our application, we will be using $S_1$, which behaves well under heavy-tailed data. Throughout the paper, when the nature of the statistic is not essential, we write $S = S_q$.
\begin{remark}[Choice of Test Statistic] While we focus on $S_q$, other test statistics can be used as well. For example, when capturing deviations in the average effect $ T_{*}^{-1} \sum_{t=T_0+1}^T \theta_t$, it is useful to consider $S(\hat u) = T_\ast^{-1/2} \left |\sum_{t=T_0+1}^T \hat{u}_t \right |$. \qed
\end{remark}
We use (block) permutations to compute $p$-values. A permutation $\pi$ is a one-to-one mapping $\pi:\{1,\dots,T\}\mapsto\{1,\dots,T\}$. We denote the set of permutations under study as $\Pi$ and assume that $\Pi$ contains the identity map $\mathbb{I}$. We focus on two different sets of permutations: (i) the set of all permutations, which we call \emph{iid permutations}, $\Pi_{\text{all}}$, and (ii) the set of all (overlapping) \emph{moving block permutations}, $\Pi_{\to}$.\footnote{We can also consider other types of permutations; for example, the ``iid block'' permutations. Specifically, let $\{b_1, \dots, b_K\}$ be a partition of $\{1,\dots,T\}$, then we collect all the permutations $\pi$ of these blocks, forming the ``iid m-block'' permutations $\Pi_{mb}$. In our context, choosing $m=T_*$ is natural, though other choices should work as well, similarly to the choice of block size in the time series bootstrap. We refer to CWZ18 for more results on block permutations.} The elements of $\Pi_{\to}$ are indexed by $j \in \{0,1, \dots, T-1\}$, and the permutation $\pi_{j}$ is defined as
\[
\pi_{j}(i)=\begin{cases}
i+j & {\text{if}}\ i+j\leq T\\
i+j-T& {\text{otherwise}}.
\end{cases}
\]
Figure \ref{fig:illustration_permutation} presents a graphical illustration of $\Pi_{\text{all}}$ and $\Pi_{\to}$.
\begin{center}
[Figure \ref{fig:illustration_permutation} around here.]
\end{center}
The choice of $\Pi$ does not matter for the exact finite sample validity of our procedures if the residuals are exchangeable. However, $\Pi_{\text{all}}$ has more elements than $\Pi_{\to}$, allowing us to compute more precise $p$-values and to test at lower significance levels. For the asymptotic validity under estimator consistency, the choice of $\Pi$ depends on the assumptions that we are willing to impose on the stochastic shock sequence $\{u_t\}$ (cf. Section \ref{subsec:small_error}).
For each $\pi \in \Pi$, let $\hat{u}_\pi=(\hat{u}_{\pi(1)},\dots,\hat{u}_{\pi(T)})'$ denote the vector of permuted residuals.\footnote{If the estimator of $P_t^N$ is invariant under permutations of the data $\{Z_t\}$ across the time series dimension (which is the case for many estimators in Section \ref{subsec:sc_panel}), permuting the residuals $\{\hat{u}_t\}$ is equivalent to permuting the data $\{Z_t\}$.} The permutation $p$-value is defined as follows.
\smallskip
\noindent \textbf{Definition of $p$-Value.} The $p$-value is
\begin{equation}
\hat{p}=1- \hat{F}\left( S(\hat{u})\right), ~~\text{where}~ \hat{F}\left( x\right)=\frac{1}{|\Pi|}\sum_{\pi\in\Pi}\mathbf{1}\left\{ S\left(\hat{u}_\pi \right)< x\right\}.
\label{eq:p_value}
\end{equation}
We are often interested in testing pointwise hypotheses about $\theta_t$, $H_0:\theta_{t}=\theta_{t}^0$, and in constructing pointwise confidence intervals for $\theta_{t}$. Pointwise hypotheses can be tested by defining the data under the null as $\mathbf{Z}=\left(Z_1,\dots,Z_{T_0},Z_{t}\right)'$, provided that $P_t^N$ can be estimated based on $\mathbf{Z}$. Pointwise $(1-\alpha)$ confidence intervals for $\theta_{t}$ can be constructed via test inversion as described in Algorithm \ref{algo:pointwise_ci}.
\begin{algo}[Pointwise Confidence Intervals]\label{algo:pointwise_ci} (i) Choose a fine grid of $G$ candidate values $\tilde\Theta_t=\{\tilde\theta_{1t}^0,\dots,\tilde\theta_{Gt}^0\}$. (ii) For $\tilde\theta_{t}^0\in \tilde{\Theta}_t$, define $\mathbf{Z}$ for the null hypothesis $H_0:\theta_{t}=\tilde\theta_{t}^0$ and compute the corresponding $p$-value, $\hat p (\tilde\theta_{t}^0)$, using \eqref{eq:p_value}. (iii) Return the $(1-\alpha)$ confidence set $\mathcal{C}_{1-\alpha}(t)=\left\{\tilde\theta_{t}^0\in \tilde\Theta_t: ~\hat p (\tilde\theta_{t}^0)> \alpha \right\}$.
\end{algo}
\subsection{Models for Counterfactual Proxies via Synthetic Control and Panel Data}
\label{subsec:sc_panel}
The availability of control units motivates several strategies for modeling the counterfactual mean proxies $P_t^N$. We estimate $P_t^N$ based on the imputed data under the null hypothesis, $\mathbf{Z}(\theta^0)$, and write $Y_{1t}^N$ instead of $Y^I_{1t}-\theta_t^0$ to alleviate the exposition.
\subsubsection{Difference-in-Differences Methods}
\label{ex:did}
The difference-in-differences method postulates the following model for the counterfactual mean proxy \citep[e.g.,][Section 5.1]{DI16}:
$
P_t^N=\mu+\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}.
$
This model automatically embeds the identifying information. The counterfactual mean proxy can be estimated as
$
\hat{P}^N_t=\frac{1}{T}\sum_{s=1}^{T}\left(Y^N_{1s}-\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{js}\right)+\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}.
$
\subsubsection{Synthetic Control and Constrained Lasso}
\label{ex:sc}
The canonical SC method \citep[e.g.,][]{abadie2003economic,abadie10sc,abadie2015comparative} postulates the following model:
\begin{eqnarray}
P^N_t=\sum_{j=2}^{J+1}w_jY^N_{jt}, ~~\text{where}~~ w\ge 0 ~~\text{and}~~ \sum_{j=2}^{J+1} w_j=1. \label{eq:sc_model}
\end{eqnarray}
We need to impose an identification condition that allows us to identify the weights $w$, for example:\footnote{More generally, other exclusion restrictions and identifying assumptions could be used. See also \citet{abadie10sc}, \citet{ferman2019synthetic} and \citet{ferman2019properties}, who study the behavior of SC when the data are generated by a factor model.}
\begin{itemize}
\item[(SC)] Assume that the structural shocks $u_t$ for the treated unit are uncorrelated with
contemporaneous values of the outcomes, namely: $E \left( u_t Y^N_{jt} \right)= 0~\text{for}~ 2\leq j\leq J+1$.
\end{itemize}
The counterfactual is estimated as $\hat{P}^N_t=\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$.
We focus on the following canonical SC estimator for $w$:\footnote{This formulation of canonical SC
without covariates is due to \citet{DI16}, who refer to the estimator \eqref{eq:sc_problem} as ``constrained regression''. Note that unlike \citet{DI16}, we estimate $w$ under the null hypothesis based on all the data. We focus on the canonical problem \eqref{eq:sc_problem} for concreteness. \citet{abadie10sc,abadie2015comparative} consider a more general version that also includes covariates into the estimation of the weights. Our inference method also works in conjunction with more recently proposed modified versions of SC, such as the augmented SC estimator of \citet{benmichael2018augmented}.}
\begin{eqnarray}
\hat{w}=\arg\min_{w} \sum_{t=1}^{T}\left(Y^N_{1t}-\sum_{j=2}^{J+1}w_{j}Y^N_{jt}\right)^2~~ \text{s.t.}~~ w\ge 0 ~~\text{and}~~ \sum_{j=2}^{J+1}w_j=1\label{eq:sc_problem}.
\end{eqnarray}
As an alternative, we can consider the more flexible model\footnote{The idea to relax the non-negativity constraint on the weights is not new. It first appeared in \citet{hsiao2012panel}, who compared their factor model approach to SC, and also in \citet{valero2015synthetic}, who used the cross-validated Lasso to estimate the weights, and in \citet{DI16}, who used cross-validated Elastic Net for estimation of weights. They do not establish the formal properties of these estimators. Here we emphasize another version of relaxing SC, model \eqref{ex:cLasso}, which leads to constrained Lasso \eqref{eq:cLasso}. Constrained Lasso demonstrates an excellent theoretical and practical performance: it is tuning-free, performs very well empirically and in simulations, and we prove that it is consistent for dependent data without any sparsity conditions on the weights and that it satisfies the estimator stability condition required for validity under misspecification. We emphasize that this estimator generally differs from the cross-validated Lasso estimator.}
\begin{equation}\label{ex:cLasso}
P^N_t=\mu + \sum_{j=2}^{J+1}w_jY^N_{jt}, ~~ \text{where}~~ \| w\|_1 \leq 1,
\end{equation}
maintaining the same identifying assumption (SC). The counterfactual is estimated as $\hat{P}^N_t=\hat\mu+\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$
by the $\ell_1$-constrained least squares estimator, or constrained Lasso \citep[e.g.,][]{raskutti2011minimax}:
\begin{eqnarray}
(\hat\mu,\hat{w})=\arg\min_{(\mu,w)} \sum_{t=1}^{T}\left(Y^N_{1t}-\mu-\sum_{j=2}^{J+1}w_{j}Y^N_{jt}\right)^2 ~~\text{s.t.}~~ \| w\|_1\le 1.\label{eq:cLasso}
\end{eqnarray}
The advantage over other penalized regression methods discussed next is that constrained Lasso is essentially tuning-free, does not rely on any sparsity conditions, and is valid for dependent data under weak assumptions. Moreover, constrained Lasso encompasses both difference-in-differences and canonical SC as special cases (by setting $w=(1/J,\dots,1/J)'$ and $\mu=0,w\ge 0$, respectively) and, thus, provides a unifying approach for the regression-based estimation of $P^N_t$.
Section \ref{sub: low level SC} provides primitive conditions that guarantee that the SC and the constrained Lasso estimators are valid in our framework in settings with potentially many control units (large $J$). Finally, we note that it is straightforward to incorporate (transformations of) covariates $X_{jt}$ into the estimation problems \eqref{eq:sc_problem} and \eqref{eq:cLasso}.
\subsubsection{Penalized Regression Methods}
\label{ex:pen_reg}
Consider a linear model for $P_t^N$:
$
P^N_t=\mu+\sum_{j=2}^{J+1}w_jY^N_{jt}.
$
We maintain the identifying assumption (SC). The counterfactual is estimated by $\hat{P}^N_t=\hat\mu+\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$,
where
\begin{eqnarray}
(\hat\mu,\hat{w})=\arg\min_{(\mu,w)} \sum_{t=1}^{T}\left(Y^N_{1t}-\mu-\sum_{j=2}^{J+1}w_{j}Y^N_{jt}\right)^2 + \mathcal{P}(w)\label{eq:restr_reg},
\end{eqnarray}
and $\mathcal{P}(w)$ is a penalty function that penalizes deviations away from zero. If it is desired
to penalize deviations away from other focal points $w^0$, for example, $w^0 =(1/J, \dots , 1/J)$ used in the difference-in-differences approach, we may always use instead: $\mathcal{P}(w) \leftarrow \mathcal{P}(w- w^0)$. Note that it is straightforward to incorporate covariates $X_{jt}$ into the estimation problem \eqref{eq:restr_reg}.
Different variants of $\mathcal{P}(w)$ can be considered. Examples include: Lasso \citep{tibshirani96regression}, where $\mathcal{P}(w)=\lambda \|w\|_1$ and $\lambda$ is a tuning parameter; Elastic Net \citep{Zou2005}, where $\mathcal{P}(w)=\lambda \left( (1-\alpha)\|w\|_2^2+\alpha \|w\|_1 \right)$ and $\lambda$ and $\alpha$ are tuning parameters; Lava \citep{chernozhukov2017lava}, where $\mathcal{P}(w)=\inf_{a + b = w}\lambda \left((1-\alpha) \|a\|_2^2+\alpha \|b\|_1 \right)$ and $\lambda$ and $\alpha$ are tuning parameters.
In the context of CSC methods, Lasso was used by \citet{valero2015synthetic}, \citet{li2017estimation}, and \citet{carvalho2018arco}, while \citet{DI16} proposed to use Elastic Net. We will impose only weak requirements on the performance of the estimators (pointwise consistency and consistency in prediction norm), which implies that these estimators are valid in our framework under any set of sufficient conditions that exists in the literature.
\subsubsection{Interactive Fixed Effects, Factor, and Matrix Completion Models}
\label{ex:interactive_fe}
Consider the following interactive FE model for treated and untreated units:
\begin{eqnarray}
\begin{array}{l}
Y^N_{jt}=\lambda_{j}'F_{t}+ X_{jt}'\beta + u_{jt},\quad \text{for}\quad 1\leq j \leq J+1~\text{and}~ 1\leq t\leq T,
\end{array}\label{eq:interactive_feM}
\end{eqnarray}
where $F_t$ are unobserved factors, $\lambda_j$ are unit-specific factor loadings, and $\beta$ is a vector of common coefficients. Model \eqref{eq:interactive_feM} nests the classical factor model when $\beta=0$ and also covers the traditional linear FE model, in which $\lambda_j'F_t = \lambda_j+F_t$. Consider the following assumption.
\begin{itemize}
\item[(FE)] Assume that $u_{jt}$ is uncorrelated with $(X_{jt},F_{t},\lambda_j)$,
as well as other identification conditions in \citet{bai2009panel}.
\end{itemize}
The model leads to the following proxy:
\begin{eqnarray}
P_t^N=\lambda_1'F_t + X_{1t}'\beta.
\label{eq:interactive_fe}
\end{eqnarray}
Counterfactual proxies are estimated by $\hat{P}_t^N=\hat\lambda_1'\hat{F}_t + X_{1t}'\hat \beta$, where $\hat{\lambda}_1$ and $\hat{F}_t$, and $\hat \beta$ are obtained using the alternating least squares method applied to the model \eqref{eq:interactive_feM}; see, for example, \citet{bai2009panel} and \citet{hansen2019factor} for a version with high-dimensional covariates.
\citet{hsiao2012panel} appears to the be first work that proposed the use of factor models for predicting the (missing) counterfactual responses specifically in SC settings. \citet{gobillon2016regional} and \citet{xu2017generalized} employ \citet{bai2009panel}'s estimator in this setting, albeit provide no formally justified inference methods. Formal inference results for interactive FE and factor models in SC designs are developed in \citet{chan2016policy} and \citet{li2018inference} among others.\footnote{Factor models are widely used in macroeconomics for causal inference and prediction; see, for example, \citet{stock2016factor} and the references therein. In microeconometrics, factor models are used for estimation of treatment/structural effects; see, for example, \citet{hansen2019factor} who use interactive FE models to estimate the effect of gun prevalence on crime.}
Other recent applications to predicting counterfactual responses include \citet{amjad2018robust} and \citet{athey2018matrix} (using, respectively, singular value thresholding and the nuclear norm penalization).\footnote{Note that \citet{athey2018matrix}'s analysis applies to a broader collection of problems with general missing data patterns, nesting SC and difference-in-difference problems as special cases.} Our method delivers a way to perform valid inference for policy effects using any of the factor model estimators used in these proposals applied to the complete data under the null.\footnote{Note that in our case the sharp null allows us to impute the missing counterfactual response and apply any of the factor estimators to estimate the factor model
for the entire data, which is then used for conformal inference. Hence our inference approach does
not provide inference for the counterfactual prediction methods given in those papers. Indeed, there, the missing data entries are being predicted using factor models, whereas in our case the missing data entries are known under the null, and we use any form of low-rank approximation or interactive FE model to estimate the model for the entire data under the null hypothesis.} We shall be focusing on \citet{bai2009panel}'s alternating least squares estimator\footnote{We choose to focus on PCA/SVD and the alternating least squares estimator for the following reasons: (1) they are by far the most widely used in practice, (2) the alternating least squares estimator is computationally attractive and easily accommodates unbalanced data.} and on matrix completion via nuclear norm penalization when verifying our conditions.
\subsection{Models for Counterfactual Proxies via Time Series and Fused Models}
\label{subsec:ts_fused}
\subsubsection{Simple Time Series Models}
If no control units are available, one can use time series models for the single unit exposed to the intervention. For example, consider the following autoregressive model:\footnote{We can also add a moving average component for the errors, but we do not do so for simplicity.}
\begin{equation}
\begin{array}{l}
Y_{1t}^N - \mu = \rho (Y_{1{(t-1)}}^N -\mu ) + u_{t}\\
Y_{1t}^I - \mu = \rho (Y_{1{(t-1)}}^N - \mu ) + \theta_t + u_t
\end{array} \Bigg | \quad E( u_t) = 0, \quad \{u_t\} ~\text{iid}, \quad t=1,\dots,T. \label{eq:model_fused}
\end{equation}
In model \eqref{eq:model_fused}, the mean unbiased proxy is given by
$
P_t^N = \mu + \rho (Y_{1{(t-1)}}^N- \mu).
$
Note that the policy effect here is transitory, namely it does not feed-forward itself on the future values of $Y_{1t}^I$ beyond the current values.\footnote{We leave the model with persistent feed-forward effects,
$Y_{1t}^I = \rho (Y_{1(t-1)}^I) + \theta_t + u_t $, to future work. } Under the null hypothesis, we can impute the unobserved counterfactual as $Y_{1t}^N = Y_{1t} - \theta_t$ and estimate the model using traditional time series methods, and we can conduct inference by permuting the residuals.
The simplest form of the autoregressive model is the AR($K$) process, where the $\rho(\cdot)$ take the form: $\rho ( \cdot ) = \sum_{k=0}^K \rho_k \mathrm{L}^k (\cdot ),$ where $\mathrm{L}$ is the lag operator. There are many identifying conditions
for these models, see, for example, \citet{Hamilton1994} or \citet{brockwell2013time}. More generally, we can use a nonlinear function of lag operators, $ \rho (\cdot ) = m( \cdot, \mathrm{L}^1 (\cdot), \dots, \mathrm{L^k} (\cdot ))$, as, for example, when applying neural networks to time series data \citep[e.g.,][]{chen1999improved,chen2001semiparametric}, and we refer to the latter for identifying conditions.
\subsubsection{Fused Time-Series/Panel Models}
A simple and generic way to combine the insights from the panel data and time series models is as follows. Consider the system
of equations:
\begin{equation}\label{eq: fused model AR error}
\begin{array}{l}
Y_{1t}^N = C_t^N + \varepsilon_t \\
Y_{1t}^I = C_t^N + \theta_t + \varepsilon_t \\
\end{array} \Bigg |
\begin{array}{l}
\varepsilon_t = \rho (\varepsilon_{t-1}) + u_t, ~ \{u_t\} ~\text{iid},~ E( u_t) = 0, \\
\{u_t\} \textrm{ is independent of } \{C_t^N\},
\end{array}
\Bigg | \quad t=1,\dots,T,
\end{equation}
where $C_t^N$ is a panel model proxy for $Y_{1t}^N$, identified by one of the panel data methods. Note that the model has the autoregressive formulation:
$
Y_{1t}^N = C_t^N + \rho (Y_{1(t-1)}^N - C_{t-1}^N) + u_t,
$
thereby generalizing the previous model.
Here the mean unbiased proxy for $Y_{1t}^N$ is given by $P_t^N = C_t^N + \rho (\varepsilon_{t-1})$.
$P_t^N$ is a better proxy than $C_t^N$ because it provides an additional noise reduction through prediction of the stochastic shock by its lag. The model combines any favorite panel model $C_t^N$ for counterfactuals with a time series model for the stochastic shock model in a nice way: we can identify $C_t^N$ under the null by ignoring the time series structure, and then identify the time series structure of the residuals $Y_{1t}^N - C_t^N$. Estimation can proceed analogously. This approach will often improve the size accuracy of our inferential procedures.
\section{Theory}
\label{sec:theory}
When the data are iid (or exchangeable), our procedure is exactly valid in finite samples as shown in the Appendix. In this section, we establish the validity of our inference methods with time series data. Our results are non-asymptotic in nature and, hence, hold in {\it finite samples}. Finite sample bounds are provided for the size properties of our procedure; these bounds imply that our approach is exact as $T_0\rightarrow \infty$. In Section \ref{subsec:small_error}, we establish the validity of our procedure when the estimator of $P_t^N$ satisfies weak and easy-to-verify small error conditions (pointwise consistency and consistency in the prediction norm). This result accommodates non-stationary data and only requires stationarity and weak dependence of the stochastic shock process $\{u_t\}$. In Section \ref{subsec:stability}, we consider a setting that accommodates misspecification and inconsistent estimators. We show that if the data are stationary and weakly dependent, our procedure is valid, provided that the estimators are stable.
\subsection{Approximate Validity under Estimator Consistency}
\label{subsec:small_error}
The main condition underlying the results in this section is the following assumption on the stochastic shock process.
\begin{assumption}[Regularity of the Stochastic Shock Process]
\label{ass:u} Assume that the density function of $S(u)$ exists and is bounded, and that the stochastic process $\{u_t\}_{t=1}^T$ satisfies one of the following conditions.
\begin{enumerate} \setlength{\itemsep}{0pt} \setlength{\parskip}{0pt}
\item $\{u_t\}_{t=1}^T$ are iid, or \label{ass:u_iid}
\item $\{u_t\}_{t=1}^T$ are stationary, strongly mixing, with sum of mixing coefficient bounded by $M$. \label{ass:u_weak_dependent}
\end{enumerate}
\end{assumption}
Assumption \ref{ass:u} allows the data to be non-stationary and exhibit general dependence patterns. Assumption \ref{ass:u}.\ref{ass:u_iid} of iid shocks is our first sufficient condition. Under this condition, we will be able to use iid permutations, giving us a precise estimate of the $p$-value. The iid assumption can be replaced by Assumption \ref{ass:u}.\ref{ass:u_weak_dependent}, which holds for many commonly encountered stochastic processes such as ARMA and GARCH. It can be easily replaced by an even weaker ergodicity condition, as can be inspected in the proofs. Under this assumption, we will have to rely on the moving block permutations.
\begin{remark}[Heteroscedasticity] Assumption \ref{ass:u} does not rule out conditional heteroscedasticity in the stochastic shock process $\{u_t\}$. Unconditional heteroscedasticity is allowed in $\{Z_t\}$ but not in $\{u_t\} $. When we suspect unconditional heteroscedasticity in $\{u_t\}$, we can apply another filter or model to obtain ``standardized residuals'' from $ \{\hat{u}_t\}$. This will generally require another layer of modeling assumptions, leading to an overall procedure that reduces the data to ``fundamental'' shocks that are assumed to be stationary under the null. \qed
\end{remark}
We also impose the following condition on the estimation error under the null hypothesis. Let $P^N=(P^N_1,\dots,P^N_T)'$ and $\hat P^N=(\hat{P}^N_1,\dots,\hat{P}^N_T)'$.
\begin{assumption}[Consistency of the Counterfactual Estimators under the Null]\label{ass:high_level_est_error} Let there be sequences of constants $\delta_T$ and $\gamma_T$ converging to zero. Assume that with probability $1- \gamma_T$,
\begin{enumerate} \setlength{\itemsep}{0pt} \setlength{\parskip}{0pt}
\item the mean squared estimation error is small, $\| \hat P^N - P^N \|^2_{2}/T \leq\delta^2_{T}$;
\item for $T_{0}+1\leq t\leq T$, the pointwise errors are small, $|\hat P^N_t- P^N_t| \leq\delta_{T}$.
\end{enumerate}
\end{assumption}
Assumption \ref{ass:high_level_est_error} imposes weak and easy-to-verify conditions on the performance of the estimators $\hat P_t^N$ of the counterfactual
mean proxies $P_t^N$. These conditions are readily implied by the existing results for many estimators discussed in Section \ref{sec:conformal_inference}. In Section \ref{sec:primitive_conditions}, we provide explicit primitive conditions and references to primitive conditions implying Assumption \ref{ass:high_level_est_error}.\footnote{While our general results in this section are non-asymptotic, some of the analysis in Section \ref{sec:primitive_conditions} will not be non-asymptotic in nature.}
\begin{thm}[Approximate Validity under Consistent Estimation]
\label{thm:approximate_validity}
Assume that $T_*$ is fixed. Suppose that Assumptions \ref{ass:dgp} and \ref{ass:high_level_est_error} hold. Impose Assumption \ref{ass:u}.\ref{ass:u_iid} if $\Pi=\Pi_{\text{all}}$; impose Assumption \ref{ass:u}.\ref{ass:u_weak_dependent} if $\Pi=\Pi_{\to}$. Assume $S(u)$ has a density function bounded by $D$ under the null. Then, under the null hypothesis, the p-value is approximately unbiased in size:
\[
|P\left(\hat{p}\leq\alpha\right)- \alpha| \leq C ( \tilde \delta_T + \delta_T + \sqrt{\delta_T} + \gamma_T),
\]
where $\tilde \delta_T = (T_*/T_0)^{1/4}(\log T)$ and the constant $C$ depends on $T_*$, $M$ and $D$, but not on $T$.
\end{thm}
The above bound is non-asymptotic, allowing us to claim uniform validity with respect to a rich variety of data generating processes. Using simulations and empirical examples, we verify that our tests have good power and generate meaningful empirical results. There are other considerations that also affect power. For example, the better the model for $P_t^N$, the less variance the stochastic shocks will have, subject to assumed invariance to the policy. The smaller the variance of the shocks, the more powerful the testing procedure will be.
\subsection{Approximate Validity under Estimator Stability}
\label{subsec:stability}
Misspecification is an important practical concern, and consistency of the estimators of the counterfactual mean proxies $P_t^N$ may be questionable in certain settings. The classical analysis of misspecification focuses on convergence to pseudo-true values \citep[e.g.,][]{white1996estimation}. If it is possible to show that the estimator of the counterfactual mean proxy, $\hat{P}_t^N$, is consistent for some pseudo-true value $P_t^{N\ast}$ and that $\left\{Y_{1t}^N-P_t^{N\ast}\right\}_{t=1}^T$ is stationary and weakly dependent, the theoretical results in Section \ref{subsec:small_error} imply the validity of our procedure. Pseudo-true consistency can often be verified for low-dimensional models, but consistency results under misspecification remain elusive in high-dimensional settings. Therefore, we consider a notion of approximate exchangeability, which only requires the estimator to be \emph{stable} instead of consistent for a pseudo-true value. This stability condition does not require $\hat{P}_t^N$ to be consistent for anything, nor does it rely on correct specification of the counterfactual mean proxies. In the Appendix, we illustrate the difference between consistency and stability based on the analytically tractable example of Ridge regression.
The basic idea underlying the theoretical analysis here is as follows. If the estimators are non-random or independent of the data, then stationarity and weak dependence of the data would mean that $\hat{p}$ based on moving block permutations approximately has a uniform distribution under the null. This result follows from uniform laws of large numbers for dependent data. However, in practice, the estimators are computed using the data and are thus not independent of the data. Our key insight is that stable estimators are approximately independent of individual observations.
We now formalize the notion of stability of an estimator. To emphasize the dependence of $S(\hat{u}) $ on the estimator, with a slight abuse of notation, we write $S(\mathbf{Z},\beta)=\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\beta)$.
Let $\{\tilde{Z}_{t}\}_{t=1}^{T}$ be iid from the distribution
of $Z_{1}$ and independent of $\mathbf{Z}$. For any $H\subset\{1,\dots,T\}$,
let $Z_{t,H}=Z_{t}\mathbf{1}\{t\notin H\}+\tilde{Z}_{t}\mathbf{1}\{t\in H\}$,
and $\mathbf{Z}_{H}=\{Z_{t,H}\}_{t=1}^{T}$. Hence, $\mathbf{Z}_{H}$ is a perturbed
version of $\mathbf{Z}$ under $H$, i.e., $\mathbf{Z}$ with elements in $H$ replaced
by $\{\tilde{Z}_{t}\}_{t\in H}$.
By stability, we mean that the estimator computed using $\mathbf{Z}$ is similar to that computed using $\mathbf{Z}_H$ for $H\in \mathbb{H}$. Let $R\in\mathbb{N}$ and define $m=\left\lfloor T_{0}/R\right\rfloor $. The class $\mathbb{H}=\{\widetilde{H}_{1},\dots,\widetilde{H}_{R}\}$ contains $R$ members with $|\widetilde{H}_j|\leq 3m$ elements. The plan is to require stability under $R\asymp T_0/\log(T_0)$ (so $|\widetilde{H}_j|\asymp \log(T_0)$). Since $\log(T_0)\ll T_0$, swapping out $O(\log(T_0))$ out of $T_0+T_*$ data points should not cause a large change in the estimator for reasonable estimators.
We now give precise definitions of sets in $\mathbb{H}$. For $j\in\{1,\dots,R\}$, let $H_{j}=\{(j-1)m+1,\dots,jm\}$. Since the test statistic depends on $T_*$ data points after obtaining the estimator, defining $\mathbb{H}$ to be $\{H_1,\dots,H_R \}$ is not enough for technical arguments; we need a ``wedge'' to ensure that these $T_*$ data points do not cause a problem. To do so, we enlarge $H_j$ as follows. Let $k\in\mathbb{N}$
satisfy $T_{*}<k<m$. We let $\widetilde{H}_{j}$ denote the $k$-enlargement
of $H_{j}$, i.e., $\widetilde{H}_{j}=\{s: \min_{t\in H_{j}}|s-t|\leq k\}$.
Note that $\widetilde{H}_{j}=\{(j-1)m+1-k,\dots,jm+k\}$ for
$2\leq j\leq R-1$, $\widetilde{H}_{1}=\{1,\dots,m+k\}$ and $\widetilde{H}_{R}=\{(R-1)m+1-k,\min\{Rm+k,T\}\}$.
\begin{assumption}[Estimator Stability]
\label{assu: stability}Let $\Pi= \Pi_{\to}$. There exist non-decreasing functions $\varrho_{T}(\cdot)$
such that
$
P\left(\max_{\pi\in\Pi}\left|S\left(\mathbf{Z}^{\pi},\hat{\beta}(\mathbf{Z})\right)-S\left(\mathbf{Z}^{\pi},\hat{\beta}(\mathbf{Z}_{H})\right)\right|\leq\varrho_{T}(|H|)\right)\geq1-\gamma_{1,T}
$
and \\
$
P\left(\max_{\pi\in\Pi}\left|S\left((\dot{\mathbf{Z}})^{\pi},\hat{\beta}(\mathbf{Z})\right)-S\left((\dot{\mathbf{Z}})^{\pi},\hat{\beta}(\mathbf{Z}_{H})\right)\right|\leq\varrho_{T}(|H|)\right)\geq1-\gamma_{1,T}
$
for any $H\in\{\widetilde{H}_{1},\dots,\widetilde{H}_{R}\}$, where $\dot{\mathbf{Z}}\overset{d}{=}\mathbf{Z}$
and $\dot{\mathbf{Z}}$ is independent of $(\mathbf{Z},\{\tilde{Z}_{t}\}_{t=1}^{T})$.
\end{assumption}
Assumption \ref{assu: stability} specifies the estimator stability condition. It strengthens the perturb-one sensitivity of \citet[][Assumption A.3]{lei2017distributionfree}. When the model is misspecified, Assumption \ref{assu: stability} holds whenever the estimator $\hat{\beta}(\mathbf{Z})$ is consistent to a pseudo-true parameter value. However, it is more general in that the estimator $\hat{\beta}(\mathbf{Z})$ need not converge to any non-random quantity as long as it is stable under perturbations in a few observations. This feature is crucial in our setting as it allows us to accommodate high-dimensional CSC methods for many of which consistency results under misspecification are not available. Primitive sufficient conditions for Assumption \ref{assu: stability} are provided in the Appendix.
Let $\Psi(x;\beta)=P(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\beta)\leq x)$. Our strategy is to show that, under the null hypothesis,
$\hat{F}(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z})))$ is approximately uniform on $(0,1)$. We exploit the stability condition in Assumption \ref{assu: stability} and show that $\hat{F}(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z})))$ can be approximated by $\Psi\left(\phi(\bar{Z}_{T_{0}+1},\dots,\bar{Z}_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}}));\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}})\right) $, which has the uniform distribution on (0,1). Here $(\bar{Z}_{T_{0}+1},\dots,\bar{Z}_{T_{0}+T_{*}})$ has the same distribution as $(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}})$ and is independent of $\mathbf{Z}_{\widetilde{H}_{R}}$. This essentially confirms the above intuition that for stable estimators, $\hat{\beta}(\mathbf{Z})$ is almost independent of the last few observations $(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}})$.
We impose the following regularity conditions on the data.
\begin{assumption}[Regularity of the Data]
\label{assu: regularity resid}The data under the null, $\{Z_{t}\}_{t=1}^{T}$, are stationary
and $\beta$-mixing with coefficient $\beta_{{\rm mixing}}(\cdot)$ satisfying
$\beta_{{\rm mixing}}(i)\leq D_{1}\exp(-D_{2}i^{D_{3}})$ for some constants $D_{1},D_{2},D_{3}>0$.
For $1\leq j\leq R$, there exist sequences $\xi_{T}>0$ and $\gamma_{2,T}=o(1)$ such that
$P\left(\sup_{x\in\mathbb{R}}\left|\partial\Psi\left(x;\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}})\right)/\partial x\right|\leq\xi_{T}\right)\geq1-\gamma_{2,T}$.
\end{assumption}
Stationarity and $\beta$-mixing are commonly
imposed conditions on time series data. For a large class of Markov chains, GARCH and various
stochastic volatility models, $D_{3}=1$ \citep[cf.][]{carrasco2002mixing}.
Let $(\dot{Z}_{T_{0}+1},\ldots,\dot{Z}_{T_{0}+T_{*}})$ be an independent copy of $ (Z_{T_{0}+1},\ldots,Z_{T_{0}+T_{*}})$ and also independent of $(\mathbf{Z},\{\tilde{Z}\}_{t=1}^{T})$. The bounded derivative of $\Psi\left(x;\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}})\right)$ condition says that the density of $\phi(\dot{Z}_t,\ldots,\dot{Z}_{t+T_{*}-1};\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}}))$ conditional on $\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}})$ is bounded by $\xi_{T}$ with high probability.
The bounded density condition states that the distribution of the residual does not collapse into a degenerate one or one with point mass. In many cases, $\xi_{T}=O(1)$ for continuous distributions. For example, if $(Y_t,X_t)$ is jointly Gaussian and the variance of $Y_t$ given $X_t$ is bounded below by a constant, then for any $w$ satisfying the SC restrictions, the density of $Y_t-X_t'w$ is bounded by a constant that does not depend on $w$.
The following result states the approximate validity of our testing
procedure.
\begin{thm}[Approximate Validity under Estimator Stability]
\label{thm: approx exchange} Let $\Pi=\Pi_{\to}$. Suppose that Assumptions \ref{assu: stability}
and \ref{assu: regularity resid} hold. Then, under the null hypothesis, there exists a constant
$C_{1}>0$ depending only on $D_{1}$, $D_{2}$ and $D_{3}$ such
that for any $R$ with $k<\left\lfloor T_{0}/R\right\rfloor$ and $R<T_{0}/2$,
\begin{multline*}
\left|P\left(\hat{p}\leq\alpha\right)-\alpha\right|\leq C_{1}\sqrt{\xi_{T}\varrho_{T}(T_{0}/R+2k)}+C_{1}\left(T_{0}^{-1}R[\log(T_{0}/R)]^{1/D_{3}}\right)^{1/4}\\
+C_{1}\exp\left(-(k-T_{*}+1)^{1/D_{3}}\right)+C_{1}\sqrt{\gamma_{1,T}}+C_{1}\sqrt{\gamma_{2,T}}.
\end{multline*}
\end{thm}
In the theoretical arguments, we actually show a stronger result. The above bound holds for $E|P(\hat{p}\leq\alpha\mid \hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}}))-\alpha|$. Since the stability condition states that $\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}})\approx\hat{\beta}(\mathbf{Z})$, this means that $\hat{p} $ conditional on $\hat{\beta}(\mathbf{Z})$ almost has a uniform distribution on $(0,1)$; with iid or exchangeable data, $\hat{p} $ conditional on $\hat{\beta}(\mathbf{Z})$ has an exact uniform distribution. Therefore, we can view Theorem \ref{thm: approx exchange} as a result for approximate exchangeability.
Due to the exponential decay of $\beta_{{\rm mixing}}(\cdot)$, the bound in Theorem
\ref{thm: approx exchange} tends to zero if we choose
$k$ to be a slowly growing sequence and $T_{0}/R$ to be of the same order. For example, we can choose $k$ and $R$ such that $k\asymp T_{0}/R \asymp \log T_0$. Since $|\widetilde{H}_{j}|=\left\lfloor T_{0}/R\right\rfloor +2k $,
Assumption \ref{assu: regularity resid} only requires that the changes
to $S(\mathbf{Z}^{\pi},\hat{\beta}(\mathbf{Z}))$ are small if we replace only $\log T_{0}$
observations in computing $\hat{\beta}(\mathbf{Z})$. Under finite dependence, it suffices to choose $k$ and $T_0/R$ to be large enough constants. Note that $R$ is only needed in the theoretical arguments; we do not need to choose $R$ when implementing the proposed procedure.
The theoretical analysis in this section suggests that allowing for both unrestricted patterns of non-stationarity and misspecification is not possible in general. To obtain valid inferences with non-stationary data, one has to either rely on correct specification and consistency or impose assumptions on the particular structure of the non-stationarity, which allow for pre-processing the data to make them stationary.
\section{Sufficient Conditions for Consistent Estimation}
\label{sec:primitive_conditions}
In this section, we revisit the representative models of counterfactual proxies introduced in Section \ref{sec:conformal_inference}. Primitive conditions are provided to guarantee that the estimation of the counterfactual mean proxies is accurate enough for the asymptotic validity of the proposed procedure. In particular, these conditions can be used to verify Assumption \ref{ass:high_level_est_error}. The regularity conditions (e.g., bounded moments, weak serial dependence) for different models are stated in the Appendix and are commonly imposed in the literature for these models. The counterfactual mean proxies $P_t^N$ are estimated based on the imputed data under the null, $\mathbf{Z}(\theta^0)$, and we write $Y_{1t}^N$ instead of $Y^I_{1t}-\theta_t^0$ to alleviate the exposition. All the results in this section assume that $T_0\rightarrow \infty$ and $J\rightarrow \infty$ (if $J$ is present in the model).
\subsection{Difference-in-Differences} \label{sub: low level DiD}
In Section \ref{ex:did}, we have seen that the counterfactual mean proxies implied by the canonical difference-in-differences model are:
$
P_t^N=\mu+J^{-1}\sum_{j=2}^{J+1}Y^N_{jt}.
$
We consider the following estimator:
$
\hat{P}^N_t=\hat{\mu}+\frac{1}{J}\sum_{j=2}^{J+1}Y_{jt},$ where $ \hat{\mu}=\frac{1}{T}\sum_{t=1}^{T}\left(Y^N_{1t}-\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}\right)=\mu+ \frac{1}{T}\sum_{t=1}^{T} u_t.
$
Since $\hat{P}_t^N-P_t^N=\hat \mu -\mu$, Assumption \ref{ass:high_level_est_error} holds for the simple difference-in-differences model provided that $ T^{-1}\sum_{t=1}^{T} u_t = o_P(1)$, which is true under very weak conditions.
\subsection{Synthetic Control and Constrained Lasso} \label{sub: low level SC}
Several models in Section \ref{sec:conformal_inference} (including SC and constrained Lasso) imply a structure in which the counterfactual proxy is a linear function of observed outcomes of untreated units.
To provide a unified framework for these models, we use $Y$ to denote a generic vector of outcomes and $X$ to denote the design matrix throughout this section. For example, in Section \ref{sec:conformal_inference}, we set $Y=Y^N_1$ and $X=(Y^N_2,\ldots,Y^N_{J+1}) $, where $Y^N_j=(Y^N_{j1},\ldots,Y^N_{jT})' \in \mathbb{R}^T$ for $1\leq j\leq J+1 $. These models can be written as
\begin{equation}\label{eq: linear model sec 5}
Y=Xw+u,
\end{equation}
where $u=(u_1,\dots,u_{T})' \in \mathbb{R}^T$. Identification is achieved by requiring that $X$ and $u$ be uncorrelated (cf.\ Condition (SC)).
Under the framework in \eqref{eq: linear model sec 5}, different models correspond to different specifications for the weight vector $w$. For the SC model in Section \ref{ex:sc}, $w $ is an unknown vector whose elements are nonnegative and sum up to one. More generally, one can simply restrict $w$ to be any vector with bounded $\ell_1 $-norm. This is the constrained Lasso estimator.
Since $P_t^N $ is the $t$-th element of the vector $Xw$, the natural estimator is $\hat{P}_t^N $ being the $t$-th element of $X\hat{w} $, where $\hat{w}$ is an estimator for $w$. The estimation of $w$ depends on the specification. Let $\mathcal{W} $ be the parameter space for $w$. We consider the following version of the original SC estimator
\begin{equation}
\hat{w} = \arg\min_w \ \|Y-Xw\|_{2}\quad \text{ s.t. } w \in \mathcal{W} = \{ v \geq0, \| v\|_1 =1\} \label{eq:sc_estimator_compact}.
\end{equation}
The constrained Lasso estimator is
\begin{equation}
\hat{w}= \arg\min_w \|Y-Xw\|_{2}\quad \text{ s.t. } w \in \mathcal{W} = \{ v: \|v\|_{1}\leq K \}, \label{eq:classo_estimator_compact}
\end{equation}
where $K$ is bounded and $K>0$. In light of the estimator \eqref{eq:sc_estimator_compact}, a natural choice is $K=1$.
In general, we choose the parameter space $\mathcal{W} $ to be an arbitrary subset of an $\ell_1 $-ball with bounded radius. The following result gives very mild conditions under which the constrained least squares estimators are consistent and satisfy Assumption \ref{ass:high_level_est_error}.\footnote{To simplify the exposition, we do not include an intercept in Lemma \ref{lem: constrained LS}. Similar arguments could be used to prove an analogous result with an unconstrained intercept.}
\begin{lem}[Constrained Least Squares Estimators]
\label{lem: constrained LS} Consider
$\hat{w}=\arg\min_{v} \ \|Y-Xv\|_{2}$ s.t. $v\in\mathcal{W}$,
where $\mathcal{W}$ is a subset of $\{v:\|v\|_{1}\leq K\}$
and $K$ is bounded. Assume $w \in \mathcal{W}$, the data are $\beta$-mixing
with exponential speed, and other assumptions listed
at the beginning of the proof, including the identification condition (SC), then the estimator enjoys
the performance bounds stated in the proof, in particular:
$
\frac{1}{T}\sum_{t=1}^{T} (\hat{P}_{t}^{N}-P_{t}^{N} )^{2}=o_{P}(1)$ and $
\hat{P}_{t}^{N}-P_{t}^{N}=o_{P}(1)$, for any $ T_{0}+1\leq t\leq T.
$
\end{lem}
Lemma \ref{lem: constrained LS} provides several features that are important for counterfactual inference in our setup. First, we allow $J$ to be large relative to $T$. To be precise, we only require $\log J=o(T^c) $, where $c>0$ is a constant depending only on the $\beta$-mixing coefficients; see the Appendix for details. This is particularly relevant for settings in which the number of (potential) control units and the number of time periods have a similar order of magnitude as in our empirical application in Section \ref{sec:application}. Second, Lemma \ref{lem: constrained LS} does not rely on any sparsity assumptions on $w$, allowing for dense vectors. Third, compared to typical high-dimensional estimators (e.g., Lasso or Dantzig selector), our estimator does rely on tuning parameters that can be difficult to choose in times series settings. Finally, Lemma \ref{lem: constrained LS} provides new theoretical consistency results for the canonical SC estimator in settings with time series data and potentially very many control units.
\subsection{Models with Factor Structures} \label{sub: low level factor models}
The models for counterfactual proxies introduced in Section \ref{ex:interactive_fe} have factor structures. We provide estimation results for pure factor models (without regressors), factor models with regressors (interactive FE models), and matrix completion models. In this subsection, following standard notation, we let $N=J+1$.
\subsubsection{Pure Factor Models}
Recall from Section \ref{ex:interactive_fe} the standard factor model
$
Y_{jt}^{N}=\lambda_{j}'F_{t}+u_{jt},
$
where $F=(F_1,\ldots,F_{T})'\in \mathbb{R}^{T\times k} $ and $\Lambda=(\lambda_{1},\ldots,\lambda_{N})'\in \mathbb{R}^{N\times k}$ represent the $k$-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for $Y_{1t}^N $ is $P_t^N=\lambda_{1}'F_{t} $. We identify $P_t^N $ by imposing the condition that the idiosyncratic terms and the factor structure are uncorrelated (cf.\ Condition (FE)).
We use the standard principal component analysis (PCA) for estimating $P_t^N $.\footnote{Note that PCA amounts to singular value decomposition, which can be computed using polynomial time algorithms, \citep[e.g., ][Lecture 31]{trefethen1997numerical}. } Let $Y^{N}\in\mathbb{R}^{T\times N}$
be the matrix whose $(t,j)$ entry is $Y_{jt}^{N}$. We compute $\hat{F}=(\hat{F}_1,\ldots,\hat{F}_T)'\in\mathbb{R}^{T\times k}$
to be the matrix containing the eigenvectors corresponding to the
largest $k$ eigenvalues of $Y^{N}(Y^{N})'$ with $\hat{F}'\hat{F}/T=I_{k}$.
Let $\hat{\lambda}_{j}'$
denote the $j$-th row of $\hat{\Lambda}=(Y^{N})'\hat{F}/T$. Let $\hat{F}_{t}'$ denote the $t$-th row of $\hat{F}$. Our estimate for $P_{t}^N $ is $\hat{P}_t^N= \hat{\lambda}_{1}'\hat{F}_t $. The following lemma guarantees the validity of this estimator in our context under mild regularity conditions.
\begin{lem}[Pure Factor Model]\label{thm: low level pure factor}
Assume standard regularity conditions
given in \citet{bai2003inferential}, including the identification condition (FE). Consider the factor model and the principal component estimator. Then, for any $1\leq t\leq T$, as $N \to \infty$ and $T \to \infty$, we have $
\hat{P}^N_{t}-P^N_{t}=O_{P}(1/\sqrt{N} + 1/\sqrt{T} )
$
and $\frac{1}{T} \sum_{t=1}^T (\hat{P}^N_{t}-P^N_{t})^2=O_{P}(1/ N + 1/T)$.
\end{lem}
The only requirement on the sample size is that both $N$ and $T$ need to be large. Similar to Theorem 3 of \citet{bai2003inferential}, we do not restrict the relationship between $N$ and $T$. This is flexible enough for a wide range of applications in practice as the number of units is allowed to be much larger than, much smaller than, or similar to the number of time periods.
\subsubsection{Factor plus Regression Model: Interactive FE Model}
Now we study the general form of panel models with interactive FEs. Following Section \ref{ex:interactive_fe}, these models take the form
$
Y_{jt}^{N}=\lambda_{j}'F_{t}+X_{jt}'\beta+u_{jt},
$
where $X_{jt}\in\mathbb{R}^{k_{x}}$ are observed covariates and $F=(F_1,\ldots,F_{T})'\in \mathbb{R}^{T\times k} $ and $\Lambda=(\lambda_{1},\ldots,\lambda_{N})'\in \mathbb{R}^{N\times k}$ represent the $k$-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for $Y_{1t}^N $ is $P_t^N=\lambda_{1}'F_{t}+X_{1t}'\beta $. In this model, we identify the counterfactual proxy through the condition that the idiosyncratic terms are independent of the factor structure and the observed covariates (cf.\ Condition (FE)).
The two most popular estimators are the common correlated effects (CCE) estimator by \citet{pesaran2006estimation} and the iterative least squares estimator by \citet{bai2009panel}. We focus on the iterative least squares approach, but analogous results can be established for CCE estimators.
The notations for $F_t $, $\lambda_{j} $, $\hat{F}_t $ and $\hat{\lambda}_{j} $ are the same as before. We compute
$$
(\hat{F},\hat{\Lambda},\hat{\beta})=\underset{F,\Lambda,\beta}{\arg \min } \sum_{t=1}^T \sum_{j=1}^{N} (Y^N_{jt}-X_{jt}'\beta-F_t'\lambda_{j} )^2\quad \text{ s.t. }\quad F'F/T=I_k \quad \Lambda'\Lambda = \text{Diagonal}_k.
$$
The estimate for $P_t^N $ is $\hat{P}_t^N=\hat{\lambda}_{1}'\hat{F}_t+X_{1t}'\hat{\beta} $. The following result states the validity of applying this estimator in conjunction with our inference method.
\begin{lem}[Interactive FE Model]
\label{thm: low level interactive FE} Assume the standard conditions
in \citet{bai2009panel}, including the identification condition (FE).
Then, for any $1\leq t\leq T$,
$
\hat P^N_{t}-P^N_{t}=O_{P}(1/\sqrt{T} + 1/\sqrt{N})$ and $ \frac{1}{T}\sum_{t=1}^{T}(\hat{P}^N_{t}-P^N_{t})^{2}=O_{P}(1/T+ 1/N).
$
\end{lem}
Under the conditions in Theorem 3 of \citet{bai2009panel}, $N$ is of the same order as $T$ so that rate is really $T^{-1/2}$;
however, the stated bound should hold more generally.
\subsubsection{Matrix Completion via Nuclear Norm Regularization}
Suppose that
\begin{equation}
Y_{jt}^{N}=M_{jt}+u_{jt},\quad{\rm for}\ 1\leq j\leq J+1\ {\rm and}\ 1\leq t\leq T, \label{eq: matrix completion model}
\end{equation}
where $M_{jt}$ is the $(j,t)$-element of an unknown matrix $M\in\mathbb{R}^{(J+1)\times T}$ satisfying $\|M\|_{*}\leq K$, where $\|\cdot\|_{*}$ denotes the nuclear norm (the sum of singular values). We observe $Y_{jt}^N$ for $(j,t)\in \{1,\dots,T\}\times \{1,\dots,J+1\} \backslash \{(1,t):T_0+1\leq t\leq T\} $. The identifying condition is that $E(u\mid M)=0$ and that conditional on $M$, $\{u_{j}\}_{j=1}^{J+1}$ is independent across $j$, where $u_{j}=(u_{j1},\dots,u_{jT})'\in\mathbb{R}^{T}$. The counterfactual proxy is $P_{t}^N=M_{1t} $ for $1\leq t\leq T$.
The main challenge is to recover the entire matrix $M$ despite the missing entries $\{Y_{1t}^N:T_0+1\leq t\leq T\} $. The literature on matrix completion considers the model \eqref{eq: matrix completion model} under the assumption of missingness at random and exploits the assumption that the rank of $M$ is low.\footnote{See, for example, \citet{candes2009exact}, \citet{recht2010guaranteed}, \citet{candes2011tight}, \citet{koltchinskii2011nuclear}, \citet{negahban2011estimation}, \citet{rohde2011estimation}, and \citet{chatterjee2015matrix}.} Recently, \citet{athey2018matrix} introduce this method to study treatment effects in panel data models and point out the unobserved counterfactuals correspond to entries that are missing in a very special pattern, rather than at random. Assuming the usual low rank condition on $M$, they employ the nuclear norm penalized estimator and provide bounds on the estimation error in the typical setup of causal panel data models.
We take a different approach here since our main goal is hypothesis testing instead of estimation. The key observation is that under the null hypothesis, there are no missing entries in the data. By imposing the null hypothesis, we replace the missing entries with the hypothesized values and obtain a dataset that contains $\{Y_{jt}^N:1\leq j \leq J+1,\ 1\leq t \leq T \}$. The estimator for $M$ we examine here is closely related to existing nuclear norm regularized estimators and is defined as
\begin{align}
\hat{M}= & \underset{A\in\mathbb{R}^{N\times T}}{\arg\min}\sum_{t=1}^{T}\sum_{j=1}^{N}(Y_{jt}^{N}-A_{jt})^{2}\quad {\rm s.t. } \ \ \|A\|_{*}\leq K, \label{eq: constrained nuclear est}
\end{align}
where $K>0$ is the bound on the nuclear norm of the true matrix. In principle, it can be a sequence that tends to infinity. When $M$ represents a factor structure with strong factors, $K$ can be shown to grow at the rate $\sqrt{NT} $. Clear guidance on how to choose $K$ is still unavailable, but following \cite{athey2018matrix} one can use cross-validation.\footnote{The properties of cross-validation remain unknown in these settings.} Alternatively one can use a pilot thresholded SVD estimator to get a sense of what $K$ is, and use a somewhat larger value of $K$. The following result guarantees the validity of this estimator in our context under mild regularity conditions.
\begin{lem}
\label{lem: constrained nuclear}Consider the estimator $\hat{M}$
defined in (\ref{eq: constrained nuclear est}). Assume that $\|M\|_* \leq K $. Let the conditions listed at the beginning of the proof hold. Then, for any $T_{0}+1\leq t\leq T$,
$
\hat{P}^N_{t}-P^N_{t}=o_{P}(1)$ and $ \frac{1}{T}\sum_{t=K+1}^{T}\left(\hat{P}^N_{t}-P^N_{t}\right)^{2}=o_{P}(1).
$
\end{lem}
The result is notable because no sub-Gaussian assumptions are required. The estimator in (\ref{eq: constrained nuclear est}) does not explicitly require a low-rank condition on $M$. Instead, we impose a growth restriction on $K$. When $M$ is generated by a strong factor structure and the null hypothesis contains full information on the missing entries, we can choose $K\asymp \sqrt{NT}$ and our consistency result holds as long as $N,T\rightarrow \infty $ and $E(|u_{jt}|^{2+c}\mid M) $ is uniformly bounded for some $c>0$. In the case of weak factors, we can choose $K \ll \sqrt{NT} $ and obtain consistency.
\subsection{Time Series and Fused Models}
As pointed out in Section \ref{subsec:ts_fused}, time series models can be used to model counterfactual proxies with or without control units. We now discuss low-level conditions under which fitting these models yields estimates good enough for the purpose of our conformal inference approach.
\subsubsection{Autoregressive Models}
The linear autoregressive model with $K$ lags can be written as
$
Y_{1t}^{N}=\rho_{0}+\sum_{j=1}^{K}\rho_{j}Y_{1t-j}^{N}+u_{t},
$
where $\{u_{t}\}_{t=1}^{T}$ is an iid sequence with $E(u_{t})=0$.\footnote{Here
the model seems different, but Section \ref{subsec:ts_fused}'s model implies this one with $\rho_0 = \mu (1 - \sum_{j=1}^{K}\rho_{j})$.} The counterfactual proxy for $Y_{1t}^N $ is $P_t^N= \rho_{0}+\sum_{j=1}^{K}\rho_{j}Y_{1t-j}^{N}$. We write $P_t^N$ as $P_t^N = y_t'\rho$, where $y_{t}=(1,Y_{1t-1}^{N},Y_{1t-2}^{N},\dots,Y_{1t-K}^{N})'\in\mathbb{R}^{K+1}$
and $\rho=(\rho_0,\dots,\rho_K)'\in \mathbb{R}^{K+1}$. The coefficient vector $\rho$ can be estimated using least squares: $\hat{\rho}=\left(\sum_{t=K+1}^{T}y_{t}y_{t}'\right)^{-1}\left(\sum_{t=K+1}^{T}y_{t}Y_{1t}^{N}\right)$. The estimator for $P_t^N $ is $\hat P_t^N = y_t '\hat \rho$.
\begin{lem}[Linear AR Model]
\label{lem: low level AR}Suppose that $\{u_{t}\}_{t=1}^{T}$ is an iid sequence with $E(u_{1})=0$
and $E(u_{1}^{4})$ uniformly bounded and the roots of $1-\sum_{j=1}^{K}\rho_{j}L^{j}=0$
are uniformly bounded away from the unit circle. Then, for any $T_{0}+1\leq t\leq T$,
$
\hat{P}^N_{t}-P^N_{t}=o_{P}(1)$ and $\frac{1}{T}\sum_{t=K+1}^{T}(\hat{P}^N_{t}-P^N_{t})^{2}=o_{P}(1).
$
\end{lem}
As mentioned in Section \ref{subsec:ts_fused}, we can also apply nonlinear autoregressive models
$
Y_{1t}^N=\rho(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N)+u_{t},
$
where $\rho $ is a nonlinear function, in which case the counterfactual proxy is $P_t^N= \rho(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N)$.
Let $\hat{\rho} $ be an estimator for $\rho$ and $\hat{P}_{t}^N= \hat{\rho}(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N) $. This estimator can be parametric, semiparametric, or fully nonparametric and is only required to be consistent.
\begin{lem}[Nonlinear AR Model]
\label{lem: low level nonliear AR} Suppose that (1) $\|\hat{\rho}-\rho\|= O_P(r_T)$ with $r_T =o(1) $ for some appropriate norm $\|\cdot\| $ and $\max_{K+1 \leq t \leq T } |\hat{\rho}(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N) -\rho(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N)| \leq \ell_T \| \hat \rho - \rho \|$ for some $\ell_T r_T =o(1) $. Then, for any $T_{0}+1\leq t\leq T$,
$
\hat{P}^N_{t}-P^N_{t}=o_{P}(1)$ and $\frac{1}{T} \sum_{t=K+1}^{T}(\hat{P}^N_{t}-P^N_{t})^{2} =o_{P}(1).
$
\end{lem}
The primitive regularity conditions and the definitions of the neural network estimators possessing these properties can be found, for example, in \citet{chen1999improved} and \cite{chen2001semiparametric}.
\subsubsection{Fused Panel/Time Series Models with AR Errors}
Here we provide generic conditions for the fused panel/time series models described in Section \ref{subsec:ts_fused}.
In particular, AR models can be used to filter the estimated residuals and obtain
near iid errors. In Equation (\ref{eq: fused model AR error}) of Section \ref{subsec:ts_fused}, we introduce an autoregressive structure in the error terms:
$
Y_{1t}^N=C_t^N+\varepsilon_t$ and $\varepsilon_t=\rho(\varepsilon_{t-1})+u_t,
$
where $C_t^N $ can be specified as a panel data model discussed before. Due to the autoregressive structure in $\varepsilon_{t}$, the counterfactual proxy is $P_t^N=C_t^N+\rho(\varepsilon_{t-1}) $.
We estimate $P_t^N $ via a two-stage procedure. In the first stage, we estimate $C_t^N $ using the techniques we considered before and obtain say $\hat{C}_t^N $. In the second stage, we estimate $\rho(\varepsilon_{t-1}) $ by fitting an autoregressive model to the estimated residuals $\{\hat{\varepsilon}_t\}_{t=1}^T $, where $\hat{\varepsilon}_t=Y_{1t}^N-\hat{C}_t^N $. For simplicity, we consider a linear model in the second stage estimation. Analogous results can be obtained for more general models. To be specific, assume that $\varepsilon_{t}=x_{t}'\rho+u_{t}$, where $x_{t}=(\varepsilon_{t-1},\varepsilon_{t-2},\dots,\varepsilon_{t-K})'\in\mathbb{R}^{K}$
and $\rho=(\rho_{1},\rho_{2},\dots,\rho_{K})'\in\mathbb{R}^{K}$.
Given $\{\hat{\varepsilon}_t \}_{t=1}^T $ from the first-stage estimation, we define $\hat{x}_{t}=(\hat{\varepsilon}_{t-1},\hat{\varepsilon}_{t-2},\dots ,\hat{\varepsilon}_{t-K})'\in\mathbb{R}^{K}$
and $\hat{\rho}=\left(\sum_{t=K+1}^{T}\hat{x}_{t}\hat{x}_{t}'\right)^{-1}\left(\sum_{t=K+1}^{T}\hat{x}_{t}\hat{\varepsilon}_{t}\right)$. To compute the $p$-value, we use $\{\hat{u}_t\}_{t=K+1}^T $ with $\hat{u}_t=\hat{\varepsilon}_t-\hat{x}_t'\hat{\rho} $ in the permutation. By the following result, this procedure is valid under very mild conditions for the first-stage estimation.
\begin{lem}[AR Errors]
\label{lem: pre-whitening u hat}Suppose that $\{u_{t}\}_{t=1}^{T}$
is an iid sequence with $E(u_{t})=0$ and $E(u_{1}^{4})$ uniformly
bounded and the roots of $1-\sum_{j=1}^{K}\rho_{j}L^{j}=0$ are uniformly
bounded away from the unit circle. We assume that
(1) $\sum_{t=1}^{T}(\hat{C}^N_{t}-C^N_{t})^{2}=o_{P}(T)$,
and (2) $\hat{C}^N_{t}-C^N_{t}=o_{P}(1)$ for $T_{0}-K+1\leq t\leq T$. Then, for any $T_{0}+1\leq t\leq T$,
$
\hat{P}^N_{t}-P^N_{t}=o_{P}(1)$ and $\sum_{t=K+1}^{T}\left(\hat{P}^N_{t}-P^N_{t}\right)^{2}=o_{P}(T)
$
\end{lem}
Note that the conditions in Lemma \ref{lem: pre-whitening u hat} for the autoregressive part are the same as in Lemma \ref{lem: low level AR}. Consistency of $\hat{C}_t^N $ can be verified using existing results, for example, those in Sections \ref{sub: low level DiD}--\ref{sub: low level factor models}.
\section{Empirical Application}
\label{sec:application}
We revisit the analysis in \citet{cunningham2018decriminalizing} who study the impact of decriminalizing indoor prostitution.
They consider the case of Rhode Island, where a judge unanticipatedly decriminalized indoor sex work in July 2003 such that, until the recriminalization in November 2009, Rhode Island had decriminalized indoor and prohibited street prostitution.
We focus on the effect of legalizing indoor prostitution on female gonorrhea incidence. Our outcome of interest is log female gonorrhea incidence per 100,000. We use the data on gonorrhea cases from the Center for Disease Control (CDC)'s Gonorrhea Surveillance Program previously analyzed by \citet{cunningham2018decriminalizing}; see their Section 3 for a detailed description and descriptive statistics. The female gonorrhea series date back to 1985 such that $T_0=19$ and $T_\ast=6$. Figure \ref{fig:raw_data} displays the raw data for Rhode Island and the rest of the U.S. states.
\begin{center}
[Figure \ref{fig:raw_data} around here.]
\end{center}
We apply three different CSC methods: difference-in-differences, canonical SC, and constrained Lasso with $K=1$. Recall that constrained Lasso nests both difference-in-differences and SC. Following \citet{cunningham2018decriminalizing}, the set of potential control units includes all other U.S. states and the District of Columbia ($J=50$). We choose $S_1$ as our test statistic and report $p$-values computed based on moving block and iid permutations.\footnote{To keep computation tractable, we randomly sample 10,000 iid permutations with replacement.} All computations were performed in \texttt{R} \citep{R2020}.
Before turning to the main results, we use the placebo tests proposed in the Appendix to assess the plausibility of the underlying assumptions. Specifically, based on the pre-treatment data, we test $H_0: \theta_{2003-\tau+1}=\dots=\theta_{2003}=0$ for $\tau\in \{1,2,3\}$. Rejections of this null undermine the credibility of the assumptions underlying our procedure and the inferences on policy effects in the post-treatment period. Table \ref{tab:placebo_specification_tests} presents the results. Figure \ref{fig:placebo_graphical} complements the formal tests with plots of the residuals from fitting the three models to the pre-treatment data. The placebo tests and the residual plots provide evidence in favor of the credibility of our inference method in conjunction with SC and, especially, constrained Lasso, but suggest that the difference-in-differences results need to be interpreted with caution.
\begin{center}
[Table \ref{tab:placebo_specification_tests} around here.]
\end{center}
\begin{center}
[Figure \ref{fig:placebo_graphical} around here.]
\end{center}
Table \ref{tab:zero_effect} reports $p$-values from testing the null hypothesis of a zero effect:
\begin{equation}
H_0: \theta_{2004}=\theta_{2005}=\dots =\theta_{2009}=0 \label{eq:zero_effect}.
\end{equation}
The null hypothesis \eqref{eq:zero_effect} is rejected at the 10\% level based on both permutation schemes and all three methods.
\begin{center}
[Table \ref{tab:zero_effect} around here.]
\end{center}
Figure \ref{fig:ci_gonorrhoea} displays pointwise 90\% confidence intervals. The results are similar for all three methods. While the effect was not or only marginally significant during the first three years, legalizing indoor prostitution significantly decreased the incidence of female gonorrhea thereafter, corroborating the findings by \citet{cunningham2018decriminalizing}.
\begin{center}
[Figure \ref{fig:ci_gonorrhoea} around here.]
\end{center}
To investigate the robustness of our results, we perform a leave-one-out robustness check \citep[e.g.,][]{abadie2015comparative} to assess whether our findings are driven by a single control state. We iteratively exclude from the control group one of the states for which either the SC or constrained Lasso weights estimated based on the pre-treatment data are non-zero and compute the $p$-values for testing hypothesis \eqref{eq:zero_effect}. Figure \ref{fig:robustness} displays the distribution of the resulting $p$-values. Overall, our results are robust and not driven by a single control state: except for one specification, all results are significant at the 10\%-level.
\begin{center}
[Figure \ref{fig:robustness} around here.]
\end{center}
\if11
{
\section*{Acknowledgements}
We are grateful to Guido Imbens, Jacopo Diquigiovanni, Bruno Ferman, the Co-Editor (Matias Cattaneo), anonymous referees, and many seminar and conference participants for valuable comments. We would like to thank Scott Cunningham and Manisha Shah for sharing the data for the empirical application. W\"uthrich is also affiliated with CESifo and the Ifo Institute. Victor Chernozhukov gratefully acknowledges funding by the National Science Foundation. All errors are our own.
} \fi
\if01
{
} \fi
\spacingset{1.0}
\setlength{\bibsep}{2pt}
\bibliographystyle{apalike}
\bibliography{SC_biblio}
\section*{Figures Main Text}
\begin{figure}[H]
\caption{Small Sample Size Properties (Nominal Level: 10\%)}
\begin{center}
\includegraphics[width=0.6\textwidth,trim={0 1cm 0 2.5cm}]{imposing_null_intro}
\end{center}
{\footnotesize \textit{Notes:} Empirical rejection probability from testing $H_0:\theta_{T_0+1}=0$. The data are generated as $Y^N_{1t}=\sum_{j=2}^{J+1}w_jY_{jt}^N+u_t$, where $Y_{jt}^N\sim N(0,1)$ is iid across $(j,t)$, $\{u_t\}$ is a Gaussian AR(1) process, $(w_2,\dots,w_{J+1})'=(1/3,1/3,1/3,0,\dots,0)'$, $T_0=19$, and $J=50$. The weights are estimated using the canonical SC method (cf. Section \ref{ex:sc}).}
\label{fig:importance_H0}
\end{figure}
\begin{figure}[H]
\caption{ Graphical Illustration Permutations}
\label{fig:illustration_permutation}
\begin{center}
\begin{tikzpicture}[scale=0.65, transform shape]
\foreach \pt/\r/\ang in {1/3/0,2/3/45,3/3/90,4/3/135,5/3/180} {
\node[circle,draw,fill=white] (\pt) at (\ang:\r){\pt};
}
\foreach \pt/\r/\ang in {6/3/225,7/3/270, 8/3/315} {
\node[circle,draw,fill=gray] (\pt) at (\ang:\r){\pt};
}
\foreach \x/\y in {1/3, 2/6, 3/4, 4/7, 5/8, 6/5, 7/2, 8/1} {
\draw[->-] (\x) -- (\y);
}
\end{tikzpicture} \hspace{.3in} \begin{tikzpicture}[scale=0.65, transform shape]
\draw (0,0) circle [radius=3];
\foreach \pt/\r/\ang in {1/3/0,2/3/45,3/3/90,4/3/135,5/3/180} {
\node[circle,draw,fill=white] (\pt) at (\ang:\r){\pt};
}
\foreach \pt/\r/\ang in {6/3/225,7/3/270, 8/3/315} {
\node[circle,draw,fill=gray] (\pt) at (\ang:\r){\pt};
}
\foreach \ang in {25,70,115,160, 205, 250, 295, 340} {
\draw[->-] (\ang-1:3) -- (\ang+1:3);
}
\end{tikzpicture}
\hspace{.3in}
\end{center}
{\footnotesize \textit{Notes:} The left figure gives an example of an iid permutation of $\{1,2,3,4,5,6,7,8\}$. The right figure gives an example of a moving block permutation of $\{1,2,3,4,5,6,7,8\}$. $T_0=5$, $T_\ast=3$. Pre-treatment periods are white; post-treatment periods are gray.}
\end{figure}
\begin{figure}[H]
\caption{Raw Data}
\begin{center}
\includegraphics[width=0.65\textwidth,trim={0 1cm 0 1cm}]{gonorrhoea_data_raw.pdf}
\end{center}
\label{fig:raw_data}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. The figure shows the raw state-level data on log female gonorrhea cases per 100,000.}
\end{figure}
\begin{figure}[H]
\caption{Graphical Placebo Checks}
\begin{center}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{gonorrhoea_resid_pre_did.pdf}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{gonorrhoea_resid_pre_sc.pdf}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{gonorrhoea_resid_pre_classo.pdf}
\end{center}
\label{fig:placebo_graphical}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. The figure plots the pre-treatment residuals estimated using difference-in-differences, SC, and constrained Lasso.}
\end{figure}
\begin{figure}[H]
\caption{Pointwise Confidence Intervals}
\begin{center}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{ci_gonorrhoea_did.pdf}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{ci_gonorrhoea_sc.pdf}
\includegraphics[width=0.325\textwidth,trim={0 1cm 0 1cm}]{ci_gonorrhoea_classo.pdf}
\end{center}
\label{fig:ci_gonorrhoea}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. The figure plots pointwise 90\% confidence intervals computed using Algorithm \ref{algo:pointwise_ci}.}
\end{figure}
\begin{figure}[H]
\caption{Leave-one-out Robustness Checks}
\begin{center}
\includegraphics[width=0.45\textwidth,trim={0 1cm 0 1cm}]{robustness_mb.pdf}
\includegraphics[width=0.45\textwidth,trim={0 1cm 0 1cm}]{robustness_iid.pdf}
\end{center}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. This figure shows the distribution of $p$-values from testing null hypothesis \eqref{eq:zero_effect}, leaving-out one of the control states with non-zero weight at the time. The size of the circles is proportional to the number of $p$-values. DID: difference-in-differences; SC: synthetic control; CL: constrained Lasso.}
\label{fig:robustness}
\end{figure}
\section*{Tables Main Text}
\begin{table}[H]
\begin{center}
\caption{Placebo Specification Tests}
\label{tab:placebo_specification_tests}
\footnotesize
\begin{tabular}{lcccccc}
\\
\toprule
\midrule
&\multicolumn{3}{c}{Moving Block Permutations} & \multicolumn{3}{c}{iid Permutations} \\
\cmidrule(l{5pt}r{5pt}){2-4} \cmidrule(l{5pt}r{5pt}){5-7} \
$\tau$ &Diff-in-Diffs&Synth. Control&Constr. Lasso&Diff-in-Diffs&Synth. Control&Constr. Lasso \\
\midrule
1 & 0.11 & 0.32 & 1.00 & 0.11 & 0.31 & 1.00 \\
2 & 0.16 & 0.32 & 0.89 & 0.06 & 0.31 & 0.93 \\
3 & 0.11 & 0.26 & 1.00 & 0.03 & 0.25 & 0.95 \\
\midrule
\bottomrule
\end{tabular}
\end{center}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. Table shows $p$-values from testing $H_0: \theta_{2003-\tau+1}=\dots=\theta_{2003}=0$ for $\tau\in \{1,2,3\}$ based on the pre-treatment data.}
\normalsize
\end{table}
\begin{table}[H]
\begin{center}
\caption{Zero Effect Null Hypothesis}
\label{tab:zero_effect}
\footnotesize
\begin{tabular}{cccccc}
\\
\toprule
\midrule
\multicolumn{3}{c}{Moving Block Permutations} & \multicolumn{3}{c}{iid Permutations} \\
\cmidrule(l{5pt}r{5pt}){1-3} \cmidrule(l{5pt}r{5pt}){4-6} \
Diff-in-Diffs&Synth. Control&Constr. Lasso&Diff-in-Diffs&Synth. Control&Constr. Lasso \\
\midrule
0.08 & 0.04 & 0.08 & 0.01 & 0.03 & 0.01 \\
\midrule
\bottomrule
\end{tabular}
\end{center}
{\footnotesize \textit{Notes:} Data are from \citet{cunningham2018decriminalizing}. Table shows $p$-values from testing $H_0: \theta_{2004}=\theta_{2005}=\dots =\theta_{2009}=0$.}
\normalsize
\end{table}
\newpage