\section{Introduction}

Difference-in-Differences (DID) models present a series of challenges for inference. There is a large number of inference methods for DID, but the effectiveness of different solutions depend crucially on the set of assumptions one is willing to make on the errors, and on many features of the empirical design, such as the number of treated and control units.\footnote{A non-exhaustive list of papers that proposed and/or analyzed different inference methods in the DID settings include  \cite{Arellano}, \cite{Bertrand04howmuch}, \cite{cameron2008bootstrap},  \cite{IZA},  \cite{CT}, \cite{Donald}  \cite{FP},   \cite{Canay},  \cite{MW}, and  \cite{ferman2019simple}. } We consider a common setting in which  a satisfactory solution is not yet available, when (i) there is a small number of treated units, (ii) the number of periods is fixed, and (iii) errors are possibly spatially correlated, but there is no available distance metric across units. Throughout, we refer to ``spatial correlation'' as any correlation in the cross section, not necessarily related to a geographical distance metric.
Importantly, as we explain in detail below, our findings are relevant even for settings in which the number of treated units has been considered as large enough to rely on inference methods based on asymptotic approximations with a large number of both treated and control units.

We first derive conditions in which the inference methods proposed by \cite{CT} (henceforth, CT) and \cite{FP} (henceforth, FP) remain asymptotically valid in the presence of spatial correlation when there is a single treated and many control units. The main assumptions are that (i)  the post-pre difference in average errors for each unit has the same marginal distribution for all units --- we can relax this assumption by allowing for heteroskedasticity depending on observed covariates, as considered by FP ---, and (ii) the cross-section distribution of this post-pre difference in average errors is strongly mixing for the control units. Under these conditions, the asymptotic distribution of the DID estimator depends only on the post-pre difference in average errors of the treated unit (as originally shown by CT), and the residuals of the control units asymptotically recover the distribution of  the errors of the treated unit even when there is spatial correlation (in contrast, both CT and FP assume independence across units for this result).\footnote{{CT} have in their appendix an example model that allows spatial dependence. However, their model requires an observed distance metric across units, and correct specification of the spatial correlation structure. We consider instead a setting in which such distance metric is not available, as is common in DID applications.}  While a standard permutation test, as proposed by \cite{fisher}, would be valid under any structure on the potential outcomes \emph{if treatment is randomly assigned}, our result does not require random assignment. In particular, we find conditions in which the inference method proposed by FP is valid even when a standard permutation test would not be valid.

However, when there is more than one  treated unit, the inference methods proposed by CT and FP will not be asymptotically valid if there is spatial correlation.  The intuition is the following: when we (mistakenly) assume errors are independent across clusters, we underestimate the volatility of the average of the errors for the treated units.  We propose an alternative inference method that is asymptotically valid, though generally conservative in this setting. We also show that in common settings in which the outcome variable is the aggregate of unit $\times$ time individual-level observations, it is possible to leverage on such structure to derive an asymptotically valid test that is more powerful.

While the methods proposed by CT and FP have been usually considered as good alternatives when the number of treated units is very small, our conservative tests may provide an interesting alternative even for settings in which applied researchers have been confidently using  alternatives such as the cluster-robust variance estimator (CRVE). The rationale is the following: if we assume errors are independent across groups, then there is a trade-off between relying on CRVE relative to  CT and FP. While CRVE relies on weaker assumptions on the errors, it relies on an asymptotic theory that depends on a large number of both treated and control units. With, say, 20 or more treated units, such asymptotic approximation required by CRVE is generally considered good enough so that the benefits of relying on weaker assumptions dominate.\footnote{\cite{ferman2019simple} proposes an assessment that can be used to  evaluate such trade offs.} However, when  we allow for spatial correlation (and considering the case with more than one treated unit), all of those methods will generally be invalid. In contrast, the inference methods we propose are robust to spatial correlation, even in settings in which a distance metric is not available.

We analyze a Monte Carlo simulation based on the spatial correlation structure of the American Community Survey. Finally, we revisit the work by \cite{Sommers} and \cite{Kaestner}, who analyzed the effects of the Massachusetts 2006 health reform, in light of our results. This empirical application highlights some common misunderstandings in the use of randomization inference in DID applications, and illustrates how our results  can be used to provide proper inference in this and other more complex settings.

\section{Setting}

Let $Y_{it}(0)$ ($Y_{it}(1)$) be the potential outcome of unit $i$ at time $t$ when this unit is untreated (treated) at this period. We consider that potential outcomes are given by
\begin{eqnarray} \label{simple_did_equation}
\begin{cases} Y_{it}(0) =   \theta_i + \gamma_t + \eta_{it}  \\  Y_{it}(1) = \alpha_{it} +Y_{it}(0),  \end{cases}
\end{eqnarray}
where $\alpha_{it}$ is the (possibly heterogeneous) treatment effects for unit $i$ at time $t$, $\theta_i$ are time-invariant unobserved effects, and $\gamma_t$ are group-invariant unobserved effects. The error term $\eta_{it}$ represents unobserved determinants of $Y_{it}(0)$ that are not captured by the fixed effects. This structure for the potential outcomes, together with Assumption \ref{As_W}(i) below, is similar to assuming  parallel trends for all periods. We can extend our results to consider alternative parallel trends assumptions  \citep{Marcus} and alternative estimators.

This implies  a DID  model
\begin{eqnarray} \label{baseline}
Y_{it} = \alpha d_{it} + \theta_i + \gamma_t  + (\alpha_{it} - \alpha)d_{it}+ \eta_{it},
\end{eqnarray}
where $Y_{it}$ is the observed outcome variable for unit $i$ at time $t$, $d_{it}$ is an indicator variable equal to one if unit $i$ is treated at time $t$, and zero otherwise.
We focus on the case in which $d_{it}$ changes to 1 for all treated units starting after date $t^\ast$ (see Remark \ref{staggering} for the case  with variation in adoption time). Therefore, we can let $D_i$ be an indicator variable equal to one if unit $i$ is treated.   There are $N_1$ treated units, $N_0$ control units, and $T$ time periods. Let $\mathcal{I}_1$ ($\mathcal{I}_0$) be the set of indices for treated (control) units, while  $\mathcal{T}_1$ ($\mathcal{T}_0$) be the set of indices for post- (pre-) treatment periods.  In this case, we define $\alpha = \frac{1}{N_1} \frac{1}{T_1} \sum_{i \in \mathcal{I}_1}\sum_{t \in \mathcal{T}_1} \alpha_{it}$, where   $\{ \alpha_{i1},...,\alpha_{iT} \}_{i \in \mathcal{I}_1}$ are fixed parameters.\footnote{CT and FP consider a setting with homogeneous treatment effects, so not only they consider $\{ \alpha_{i1},...,\alpha_{iT} \}_{i \in \mathcal{I}_1}$ as fixed, but also that $\alpha_{it} =\alpha$ for all $i$ and $t$ in their settings. Therefore, we consider a more general setting than the settings considered by CT and FP.  Since we consider a setting with $N_1$ fixed, it would not be possible to assess the distribution of  $\alpha_{it}$, if we considered this term as stochastic. Therefore, it would not be possible to provide valid inference if we consider $\alpha_{it}$ as stochastic and define the estimand as its expected value. }

 Following FP, for a generic variable $A_t$, define $\nabla A = \frac{1}{T-t^\ast} \sum_{t \in \mathcal{T}_1} A_{t} -  \frac{1}{t^\ast} \sum_{t \in \mathcal{T}_0} A_{t}$. In particular,  we consider $W_i = \nabla \eta_{it}$, which is the post-pre difference in average errors for each unit $i$. In this  case in which treatment starts at the same period for all treated units, the DID estimator is numerically equivalent to the two-way fixed effects (TWFE) estimator of $\alpha$, which is given by
\begin{eqnarray}  \label{alpha}
\widehat \alpha  &=& \frac{1}{N_1} \sum_{i \in \mathcal{I}_1} \nabla Y_{i}  - \frac{1}{N_0} \sum_{t \in \mathcal{I}_0} \nabla Y_{i}= \alpha +  \frac{1}{N_1} \sum_{i \in \mathcal{I}_1} W_i - \frac{1}{N_0} \sum_{i \in \mathcal{I}_0} W_i.
\end{eqnarray}

Our results are also relevant and can be extended for alternative estimators.\footnote{For example, we can consider the first differences estimator, or the estimators proposed by \cite{chaisemartin2018twoway} and \cite{Pedro}.}

We consider  a repeated sampling framework over the distribution of $\{ W_i\}_{i \in \mathcal{I}_0 \cup \mathcal{I}_1 }$, in which treatment assignment is pre-determined (alternatively, we can consider that the analysis is conditional on the vector of treatment assignment).
For simplicity, we consider first the homoskedastic case, as analyzed by CT, and we do not include covariates in model (\ref{baseline}). In Appendix \ref{Appendix_extensions}  we allow for heteroskedasticity based on a set of covariates $Z_i$, as considered by FP, and for  a set of covariates $X_{it}$ in model (\ref{baseline}). In this simpler case, we impose the following assumptions on the distributions of $W_i$.
\begin{ass}{(distribution of $W_i$)} \label{As_W}
\normalfont
(i) $ \mathbb{E} \left[ W_i  \right] =0$ for $i \in \mathcal{I}_1$ and for $i \in \mathcal{I}_0$, (ii) $W_i$ has the same marginal distribution for all $i \in \mathcal{I}_1 \cup \mathcal{I}_0$, with finite second moment and continuous distribution, and   (iii)  $\{W_i\}_{i \in \mathcal{I}_0}$ is such that $\frac{1}{N_0} \sum_{i \in \mathcal{I}_0} W_i \buildrel p \over \rightarrow 0 $  and $\frac{1}{N_0} \sum_{i \in \mathcal{I}_0} \mathbbm{1} \left\{ W_i \leq c \right\} \buildrel p \over \rightarrow F(c) $ when $N_0 \rightarrow \infty$, where  $\mathbbm{1} \left\{ . \right\}$ is an indicator function and $F(c)$ is the cdf of $W_i$.

\end{ass}

Assumption \ref{As_W}(i) is a standard DID identification assumption. Assumption \ref{As_W}(ii)  restricts the marginal distribution of the treated and control units to be the same. This is essentially an homoskedasticity assumption, as considered by CT. Importantly, we can relax this assumption to allow for heteroskedasticity based on a set of covariates $Z_i$ (see  Appendix \ref{Appendix_extensions}).  Finally, Assumption \ref{As_W}(iii) allows for spatially correlated shocks, but restricts such dependence so that we can apply a law of large numbers when we consider the control units. This will be satisfied, for example, if we assume strong mixing conditions in the cross section (see Theorem 3 from \cite{JENISH200986}). For simplicity, we refer in the text  to Assumption \ref{As_W}(iii) as a strong mixing condition, even though strong mixing is not a necessary condition for this assumption. Importantly, Assumption \ref{As_W} allows for  unrestricted spatial correlation among the treated units.

Since we focus on settings in which the researcher does not have information on what generates the spatial correlation or he/she is not willing to assume such structure, we do not need to model in detail the sources of spatial correlation. As a concrete  example, we can considered a setting in which the $N_1$ treated units are closely located geographically, but we have a larger number of control units in different locations. In this case, if spatial correlation goes to zero when geographical distance increases, then we would have  Assumption \ref{As_W}(iii) satisfied, even though we may have arbitrarily strong spatial correlation among the $N_1$ treated units. The spatial correlation may also arise from other mechanisms not related to geographical distance (for example, units with similar industry shares may have more correlated errors). See also Appendix A2 from   \cite{ferman2019inference} for an example of a spatial correlation model in a DID setting that would satisfy this strong mixing condition. Importantly, we consider a setting in which the applied researcher may be unaware or may not have information on the relevant distance metrics in the cross section, which is common in DID applications.

Moreover, while we consider a setting in which the treatment allocation is fixed, we can think about an underlying model in which treatment assignment and errors are stochastic, and units that experience correlated shocks are more likely to receive the same treatment status. In this case, once we condition on the treatment allocation, we would have a setting in which the errors of the treated units are correlated.

From Proposition 1 from CT,  it follows from equation (\ref{alpha}) that the limiting distribution of $\widehat \alpha$ when $N_1$ is fixed and $N_0 \rightarrow \infty$ is given by $\widehat \alpha \buildrel p \over \rightarrow \alpha + \widetilde{W}$, where $\widetilde{W} \equiv {N_1}^{-1} \sum_{i \in \mathcal{I}_1} W_i$. The intuition is that, given the strong mixing condition (Assumption \ref{As_W}(iii)), the average of the control errors  converges in probability to zero as $N_0 \rightarrow \infty$. However, since the number of treated units remain fixed, the distribution of $\widehat \alpha$ will still depend on the errors of the treated units even asymptotically. In this case, the estimator is unbiased, but not consistent.

\section{Inference with spatially correlated shocks} \label{inference}

CT propose an interesting inference method in this setting by noting that the residuals $\widehat W_i$ of the control units may be informative about the distribution of $W_i$ and, therefore, about the distribution of $\widehat \alpha$. The main idea is to estimate the distribution of $W_i$ using the empirical distribution of  $\{\widehat W_i \}_{i \in \mathcal{I}_0}$. Intuitively, this works because, as $N_0 \rightarrow \infty$, $\widehat W_i \buildrel p \over \rightarrow W_i$ for $i \in \mathcal{I}_0$. Therefore, hypothesis tests and confidence intervals can be constructed based on the empirical distribution of $\{\widehat W_i \}_{i \in \mathcal{I}_0}$. They show validity of this inference method  when $W_i$ is iid (Assumption 2 in CT).\footnote{CT have in their appendix an example model that allows spatial dependence. However, their model requires an observed distance metric across units, and correct specification of the spatial correlation structure. We consider instead a setting in which such distance metric is not available, as is common in DID applications. }

 It is easy to understand the importance of this independence assumption when $N_1>1$.  Consider a simple case in which $N_1=2$, and $\{W_1,W_2\}$ is multivariate normally distributed with  correlation $\rho$.  In this case, $\widehat \alpha \buildrel p \over \rightarrow \alpha + \widetilde{W}$, where the variance of $\widetilde{W}$ is given by $2^{-1} (1 + \rho) var(W_i)$. However, given Assumption \ref{As_W}(iii), if we consider two random draws from  $\{\widehat W_i \}_{i \in \mathcal{I}_0}$ to recover the distribution of $\widetilde{W}$, the correlation between these draws would converge to zero when $N_0 \rightarrow \infty$. As a consequence, the approach proposed by CT would recover a distribution for $\widetilde{W}$ that has a variance given by $2^{-1} var(W_i)$. Therefore, the variance of $\widehat \alpha$ would be underestimated if $\rho>0$, leading to over-rejection. The same problem applies for the inference method proposed by FP.

When $N_1=1$, however, we show that the inference methods proposed by CT and FP can remain valid even if we allow for spatial correlation. The main intuition is that, under Assumption \ref{As_W}, the asymptotic distribution of $\widehat \alpha$ depends only on $W_1$, and the distribution of $W_1$ is still asymptotically approximated by the empirical distribution of $\{\widehat W_i \}_{i \in \mathcal{I}_0}$, even when errors are spatially correlated. Of course, the  strong mixing condition is crucial for this argument. To formalize that,  let $\widehat F (c) = {N_0}^{-1} \sum_{i \in \mathcal{I}_0} \mathbbm{1} \left\{ \widehat W_i \leq c \right\}$.

\begin{prop} \label{Prop_N1}
Under Assumption \ref{As_W}, as $N_0 \rightarrow \infty$, $\widehat F (c)$ converges in probability to $F(c)$ uniformly on any compact subset of the support of $W$.
\end{prop}

The proof is similar to the  proof of Proposition 2 from CT, but allowing for spatial correlation. We present details in Appendix \ref{Proof_N1}. It immediately follows that the inference method proposed by CT remains valid for the case with $N_1=1$, even when we may have spatial correlation.

\begin{cor} \label{Cor_N1}
Consider the problem of testing the null hypothesis $H_0: \alpha =\alpha_0$, for some $\alpha_0 \in \mathbb{R}$. Under Assumption \ref{As_W}, if $N_1=1$, then the inference method proposed by CT is asymptotically valid  when $N_0 \rightarrow \infty$.
\end{cor}

We show in Appendix \ref{Appendix_extensions} that these results remain valid when we include covariates in model (\ref{baseline}). We also show in Appendix \ref{Appendix_extensions} that the inference method proposed by FP remains asymptotically valid, in this setting with $N_1=1$ and spatial correlation,  when  we relax the assumption that $W_i$ has the same marginal distribution for all $i$. Therefore, we also have results that are valid when errors are heteroskedastic.

{While a standard permutation test, as proposed by \cite{fisher}, would be valid under any structure on the potential outcomes \emph{if treatment is randomly assigned}, note that Assumptions \ref{As_W} and \ref{As_W_appendix} do not require random assignment. For example, Assumption \ref{As_W_appendix} allows for a setting in which the error of the treated unit has larger variance than the error of the control units.   In this case, we show that FP would be valid, while a standard permutation test would lead to over-rejection.

As explained above, however, the inference methods proposed by CT and FP would not be valid when $N_1>1$ and errors are spatially correlated. Given this potential problem, we consider the case in which the errors of all treated units are perfectly correlated as a worst-case scenario to provide an asymptotically valid, though conservative, inference method in this setting. More specifically, from Proposition \ref{Prop_N1}, $\widehat F(c)$ asymptotically recovers the distribution of $W_i$ even when we allow for spatial correlation. Given that, instead of considering independent draws from this distribution to put in place of the treated units, we sample only one $\widehat W_i$ and place it for all  treated units, which would asymptotically recover the distribution of $\widehat \alpha$ if $\{W_i\}_{i \in \mathcal{I}_1}$ were perfectly correlated.  When treated units are not perfectly correlated, though, we would recover a distribution for $\widehat \alpha$ that has a \emph{higher} variance relative to the true distribution of $\widehat \alpha$. To formalize this idea, we consider the following high-level assumption.

\begin{ass}{(regularity condition)} \label{As_reg}
\normalfont
For any $\tau \in (0,1)$, let $c_{\tau}$ be the $\tau-$quantile of $W_i$, and define $\widetilde{W} = N_1^{-1} \sum_{i \in \mathcal{I}_1}W_i$. Then   $\tilde p \equiv Pr\left( \{ \widetilde W  \leq c_{\tau/2} \} \cup \{ \widetilde W  > c_{1-\tau/2} \} \right) \leq \tau $.

\end{ass}

What Assumption \ref{As_reg} states is that, regardless of the spatial correlation among the treated units, the probability of having extreme values for the average of the treated units, $\widetilde W$, is weakly smaller than the probability of having extreme values for a single draw of $W_i$.  While we can guarantee that $var(\widetilde W) \leq var(W_i)$ irrespectively of the spatial correlation among the treated units, this is not sufficient to guarantee Assumption \ref{As_reg} for a general distribution of $W_i$. Still, Assumption \ref{As_reg} is intuitively plausible, and it is satisfied when $\{W_i\}_{i \in \mathcal{I}_1}$ is multivariate normal. In this case, $\widetilde W$ will also be normal, and the fact that $var(\widetilde W) \leq var(W_i)$ guarantees that Assumption \ref{As_reg} holds. This is valid regardless of how strong the spatial correlation is. This assumption can also be justified in settings in which $Y_{it}$ represents the aggregate of individual-level observations in unit $i$ and time $t$, and we can apply a central limit theorem to these  averages. In this case, $\{W_i\}_{i \in \mathcal{I}_1}$ would be approximately multivariate normal. To provide more intuition on this assumption,  we present in Appendix \ref{Appendix_reg} an example in which this assumption may not hold. To engineer such example,  we had to consider a distribution of $W_i$ that is bimodal, and such that the two peaks are very far apart, which is not something we should expect in common empirical applications. In this example, when the two peaks are not very far apart, then Assumption \ref{As_reg} would still hold.

\begin{prop} \label{Prop_conservative}
Consider testing the null hypothesis $H_0: \alpha = \alpha_0$ with a significance level $\tau$ by contrasting $\widehat \alpha - \alpha_0$ with the $\tau/2$ and $1-\tau/2$ quantiles of $\widehat F(c)$.  Under Assumptions \ref{As_W} and \ref{As_reg}, as $N_0 \rightarrow \infty$, this test is asymptotically level $\tau$ for any $\tau \in (0,1)$.
\end{prop}

This result follows directly from Proposition \ref{Prop_N1} and Assumption \ref{As_reg}. Therefore, this conservative test provides a viable alternative in settings where $N_1 << N_0$ that is robust to strongly mixing spatial correlation. While we guarantee that this modified test does not over-reject (asymptotically), it will generally be conservative when $N_1 > 1$, unless the errors of treated units are perfectly correlated.
 In Appendix \ref{Appendix_extensions}, we extend the result from  Proposition \ref{Prop_conservative} for settings in which we have heteroskedasticity, as considered by FP.  Therefore, a valid test is possible even when we relax the assumption that $W_i$ has the same marginal distribution for all $i$. We also allow for covariates in model (\ref{baseline}).

This conservative test may be an interesting alternative even when $N_1$ is such that applied researchers would usually rely on large $N_1$/large $N_0$ methods, such as CRVE. Once we relax the assumption that clusters are independent, CRVE would generally be invalid even when $N_1$ is reasonably large \citep{ferman2019inference}.  In contrast, the method we propose remains valid when we have a strongly mixing spatial correlation. Of course, the cost of being robust to spatial correlation is that the conservative test will generally have lower power.

If a distance metric is available and the researcher is willing to assume such distance metric is the relevant one for the spatial correlation, then other  available alternatives might present better power (for example, the inference method proposed in the appendix of CT). In this case, there would be a trade-off between a test that is  more powerful but requires correct specification of the spatial correlation, versus a test that is less powerful, but does not require correct specification of the spatial correlation. We also recall that the test we propose does not require observation of a distance metric in the cross section.  In Section \ref{Conservative2} we consider an alternative test that will generally be more powerful in common DID applications relative to the test presented in this section, without relying on knowledge of the spatial correlation structure.

\begin{rem}
\normalfont

The assumption  that $W_i$ is strongly mixing would not be satisfied if there are unobserved shocks that affect all treated units differently from the control units. In this case, the DID residuals $\widehat W_i$ would not capture these group-specific shocks, and $\widehat F(c)$ would underestimate the dispersion of the marginal distribution of $W_i$.  As a consequence, even the conservative test may over-reject. Importantly, however, the over-rejection of the conservative test in this case would be no larger than the over-rejection of CT (or FP). In contrast,  this  mixing assumption may be reasonable  when units closer in some distance metric have more correlated errors, but such correlation goes to zero when this distance increases. Such distance metric does not need to coincide with geographical distance (for example, we may have that units with similar industry shares have more correlated errors), and we do not even need to have information on the relevant distance metrics.

\end{rem}

\begin{rem}
\normalfont \label{staggering}
Constructing a conservative test  becomes substantially more complicated if treated units start treatment at different periods, as we discuss in  Appendix \ref{Appendix_staggered}.

\end{rem}

\begin{rem}
\normalfont
Since CT and FP are valid when $N_1=1$ (Corollaries \ref{Cor_N1} and \ref{Cor_Appendix}), an alternative for the case with $N_1>1$ is to consider $N_1$ separate DID regressions --- one for each treated unit ---, and adjust the inference for multiple  testing \citep{Romano2005}. In this case, however, we would be considering $N_1$ different hypotheses related to $\alpha_i = \frac{1}{T_1}\sum_{t \in \mathcal{T}_1} \alpha_{it}$ for $i \in \mathcal{I}_1$. In contrast, the inference methods we propose in Sections \ref{inference} and \ref{Conservative2} are valid for a single hypothesis related to the aggregate parameter $\alpha = \frac{1}{N_1} \frac{1}{T_1} \sum_{i \in \mathcal{I}_1}\sum_{t \in \mathcal{T}_1} \alpha_{it}$.

\end{rem}

\section{A less conservative test}
 \label{Conservative2}

While the conservative test presented in Section \ref{inference} is valid even when we have strongly mixing spatial correlation, such test may be too conservative, implying  low power. We now show that it is possible to construct a test that  is asymptotically valid and less conservative in common settings in which  $Y_{it}$ represents averages of individual-level observations in unit $i$ and time $t$. We can consider either the case in which we observe only the aggregates $Y_{it}$ or the individual-level data. This is the leading example considered by FP, where they show that variation in unit sizes generally generate heteroskedastiticy in $W_i$. More specifically, let $M_{it}$ be the number of observations used to construct the aggregate observations for unit $i$ at time $t$. For simplicity, assume that $M_{it} = M_i$ for all $t$. FP show that, under a wide range of structures for the within-unit correlation of the individual-level errors, the variance of $W_i$ conditional on $M_i$ is given by  $var(W_i | M_i) = A + B{M_i}^{-1}$ for constants $A,B \geq 0$.

The main idea of FP is that we can use the residuals from the control units to first estimate the parameters $A$ and $B$. Then, we can use $\hat A$ and $\hat B$ to re-scale the residuals from the control units, so that they become informative about the distribution of the errors of the treated units. More specifically, instead of considering the empirical distribution of $\{\widehat W_j\}_{j \in \mathcal{I}_0}$ to approximate the distribution of $W_i$ for a treated unit  $i \in \mathcal{I}_1$, FP consider the empirical distribution of $\{\widehat W_j \times {(\hat A + \hat B M_i^{-1})^{1/2}/(\hat A + \hat B M_j^{-1})^{1/2}} \}_{j \in \mathcal{I}_0}$. FP provide conditions in which this re-scaled empirical distribution asymptotically recovers the distribution of a treated unit $W_i$ with $M_i$ observations. The problem, however, is that FP consider the case in which $W_i$ is independent across units, so their procedure would not be valid if we consider a setting with spatial correlation, following the same rationale as the one presented in Section \ref{inference}.

Consider a simple setting with $N_1=2$, where the aggregate data of the treated units  are constructed from  $M_1$ and $M_2$ individual-level observations per period, and consider a TWFE estimator using $M_i$ as sampling weights. This would lead to the same estimator as if we used the individual-level data. Note that assuming that both treated units start treatment after $t^\ast$, and given the assumption that $M_{it} = M_i$, the estimand of the TWFE will be a well-defined weighted average of the (average across time) treatment effects of units 1 and 2.  In this case, a naive version of the inference method proposed by FP ignoring spatial correlation would assume that $W_1$ and $W_2$ are independent, implying that it would asymptotically recover a distribution for $\widehat \alpha$ with variance equal to
\begin{eqnarray}
\sigma^2_{\mbox{\tiny naive}} = \left( \frac{M_1}{M_1+M_2} \right)^2 \left[ A + B{M_1}^{-1} \right] +  \left( \frac{M_2}{M_1+M_2} \right)^2 \left[ A + B{M_2}^{-1} \right].
\end{eqnarray}

As we discussed in Section \ref{inference}, this would generally lead to over-rejection if the treated units are spatially correlated. We can consider instead the  conservative method proposed in Appendix \ref{Appendix_extensions}, which is an extension of the method proposed in Section \ref{inference}. We refer to this test as ``conservative test 1.'' Asymptotically, this test recovers a distribution with variance
\begin{eqnarray}
\sigma^2_{\mbox{\tiny cons1}} = \left[ \left( \frac{M_1}{M_1+M_2} \right)\left( A + B{M_1}^{-1}  \right)^{\frac{1}{2}} +  \left( \frac{M_2}{M_1+M_2} \right) \left( A + B{M_2}^{-1}  \right)^{\frac{1}{2}}    \right]^2.
\end{eqnarray}

It is clear that $\sigma^2_{\mbox{\tiny cons1}} \geq \sigma^2_{\mbox{\tiny naive}}$. Under the assumptions considered in Appendix \ref{Appendix_extensions}, this would lead to a valid, but generally conservative test.

An alternative approach in this particular setting would be to consider that we have only one treated unit with $M_1 + M_2$ observations. In this case, we would asymptotically recover a distribution with variance equal to
\begin{eqnarray}
\sigma^2_{\mbox{\tiny cons2}} =  A + B\left({M_1+M_2}\right)^{-1} .
\end{eqnarray}

We have that  $\sigma^2_{\mbox{\tiny cons1}} \geq \sigma^2_{\mbox{\tiny cons2}} \geq \sigma^2_{\mbox{\tiny naive}}$. The main question is whether $\sigma^2_{\mbox{\tiny cons2}}$ is  conservative for the true variance of $\widetilde W = \frac{M_1}{M_1+M_2} W_1 +\frac{M_2}{M_1+M_2} W_2$. This will be the case if two individuals in the same unit are more correlated than two individuals in different treated units, which is a plausible assumption.

Intuitively, this approach would consider a scenario for the spatial correlation between individual-level observations of two different units by assuming they have the same spatial correlation as if we considered two  individual-level observations in the same unit. Assuming individuals in the same unit are weakly more spatially correlated than individuals in different units, this would be a worst-case scenario for the spatial correlation, which would be much less stringent than the one considered in the conservative test 1. Importantly, we do not need individual-level data for this approach. We only need information on the number of observations used to construct the unit $\times$ time aggregates. We refer to this test as ``conservative test 2'', and we formalize this idea in Appendix \ref{Appendix_cons2}.

Note that, if $A=0$, then $\sigma^2_{\mbox{\tiny cons2}} = \sigma^2_{\mbox{\tiny naive}}$. Therefore, in a setting in which there is no spatial correlation (whether we consider two individual-level observations in the same unit or in different units), then the conservative test 2 would be asymptotically equivalent to FP. There would not be a loss in power in this case. In contrast, the conservative test 1 would have lower power relative to FP even in this setting.

When $M_{it}$ is constant across $t$, we can implement the inference method above by generating a single treated unit that is a weighted average of the $N_1$ treated units,  $\widetilde{y}_{t} = \sum_{i \in \mathcal{I}_1}\frac{M_i}{\sum_{i' \in \mathcal{I}_1} M_{i'}}y_{it}$, with a number of observations $\widetilde{M} = \sum_{i \in \mathcal{I}_1}M_i$. In this case, the DID estimator of $\widehat{\alpha}$ using this collapsed treated unit with $\widetilde M$ and $\{M_i\}_{i \in \mathcal{I}_0}$ as weights is numerically the same as the one using the original data with $\{M_i\}_{i \in \mathcal{I}_0 \cup  \mathcal{I}_1}$ as weights. Then we can apply FP for this DID model with a single treated unit. We also consider the case in which $M_{it}$ varies across $t$ in Appendix \ref{Appendix_cons2}

\section{Simulations with Real Datasets} \label{simulations}

We analyze the spatial correlation problem, and the proposed  conservative tests, in simulations based on the  American Community Survey (ACS) \citep{ipums}. We restrict the sample to women between the ages 25 and 50, and consider log wages as outcome variable. We construct a DGP with variance/covariance matrix estimated from  this data. More specifically, with data from 2015 to 2018, we estimate the residuals $\widehat W_i$ at the PUMA level, considering a setting in which that treatment starts after 2016.\footnote{We compute the residuals from a PUMA $\times$ year regression using PUMA and time fixed effects, and then calculate  $\widehat W_i = (\hat e_{i,2018} + \hat e_{i,2017})/2 - (\hat e_{i,2016} + \hat e_{i,2015})/2$.} Given $\widehat W_i$ for all PUMAs, we estimate a model for the spatial correlation in which the covariance between two PUMAs may depend on their distance, on whether they belong to the same state, and on their population sizes (details in Appendix \ref{Appendix_MC}). If we consider all pairs of PUMAs, the correlation between their errors in this estimated model is greater than $0.05$ in only $0.42\%$ of the cases. Therefore, the strong mixing condition we consider seems reasonable in this setting. Still, we have PUMAs with spatial correlation as strong as $0.30$. Therefore, we may have a setting in which the $N_1$ treated PUMAs exhibit strong spatial correlation.

In line with our theoretical model, we consider first a setting in which we fix one PUMA as the treated one, and generate multivariate normal draws of $W_i$ with this estimated spatial dependence. For the results presented in Table \ref{Table_MC}, we selected a treated PUMA that has a strong spatial correlation with some other PUMAs. In Appendix \ref{Appendix_MC} we discuss how our results vary when we consider alternative treated PUMAs.  When we consider the test proposed by FP,  we find a rejection rate of approximately $5\%$ when we consider a $5\%-$level test when the null is true, despite the fact that PUMAs are spatially correlated (Panel A of Table \ref{Table_MC}). This is consistent with Corollary \ref{Cor_N1}. Note also that the treated PUMA is fixed and has a heterogeneous variance, so our setting is not consistent with random assignment of the treatment.

\begin{table}[h!]

 \begin{center}
\caption{{\bf MC simulations}} \label{Table_MC}

\begin{tabular}{lccc}
\hline
\hline

 & Original FP   & Conservative 1 & Conservative 2   \\
 & (1) & (2) & (3) \\ \hline

\multicolumn{4}{c}{\textit{Panel A: $N_1=1$} } \\

$\alpha=0$ & 0.054 & - & - \\
 \\
$\alpha=0.1 \sigma$ & 0.434 & - & - \\
 \\
$\alpha=0.15 \sigma$ & 0.759 & - & - \\

\\

\multicolumn{4}{c}{\textit{Panel B: $N_1=5$} } \\

$\alpha=0$ & 0.170 & 0.003 & 0.027 \\
 \\
$\alpha=0.1 \sigma$ & 0.933 & 0.440 & 0.700 \\
 \\
$\alpha=0.15 \sigma$ & 0.999 & 0.898 & 0.970 \\

\hline

\end{tabular}

 \end{center}
{\small{ Notes: This table presents rejection rates for 5\% tests using the FP inference method (column 1), the conservative test presented in Section \ref{inference} (column 2), and the conservative test presented in Section \ref{Conservative2} (column 3). Details of the simulations are presented in Appendix \ref{Appendix_MC}. In Panel A, we present results when there is a single treated unit, while in Panel B we consider the case in which there are five treated units. The true treatment effect $\alpha$ is based on the standard deviation of log wages in the individual-level data from 2018.    }}

\end{table}

We then consider a setting in which four other PUMAs that are spatially correlated with the originally treated PUMA are also treated (so $N_1=5$). In this case, the test proposed by FP leads to relevant over rejection when the null is true ($\approx 17\%$), which is consistent with the conclusions from Section \ref{inference} (Panel B of Table \ref{Table_MC}).  The two alternative tests, in contrast, control for size, although they are conservative (with rejection rates of $0.5\%$ and $2.7\%$). As expected from the discussion in Section \ref{Conservative2}, the conservative test 2 is less conservative than the conservative test 1. Moreover, the conservative test 2 is more powerful than the conservative test 1, although both tests present large power when the treatment effect is large.

Overall, these simulations highlight the main messages of the paper. First, that the inference methods proposed by CT and FP remain valid when we have a single treated unit, even when errors are spatially correlated. Second, that the proposed conservative tests control well for size. Finally, these simulations also  highlight that it is possible to explore the structure of common empirical applications in which units are aggregates of individual-level observations to construct a more powerful test (even if we do not have information on the individual-level observations). This is important, because we inherently have to accept a loss in power so that we can be robust to spatial correlation when $N_1>1$.

Of course, if we have information on the relevant distance metrics in the cross section, then it would be possible to correct for spatial correlation with a more powerful test. Still, these simulations illustrate that it is possible to construct valid tests even when such information is unavailable or the researcher is not willing to impose a structure on the spatial correlation.

\section{Empirical Illustration} \label{illustration}

We illustrate our findings analyzing the effects of the Massachusetts 2006 health care reform. This reform was analyzed by \cite{Sommers} using a DID design comparing 14 Massachusetts counties with  513 control counties that were selected based on a propensity score to be more similar with the treated counties.\footnote{The propensity score used age distribution, sex, race/ethnicity, poverty rate, median income, unemployment, uninsured rate, and baseline annual mortality as predictors. We take this first selection step as given in our analysis. We find similar results if we consider a DID regression using all counties, so that there is no pre-selection of control counties. } \cite{Sommers} find a reduction of 2.9\%-4.2\% in mortality in Massachusetts relative to the controls after the reform (depending on whether covariates are included).

\cite{Sommers} reported standard errors clustered at the state level,  and also considered standard errors clustered at the county level in their online appendix. Given the recent findings that clustered standard errors may perform badly when there are few treated clusters, \cite{Kaestner}  re-analyzed the findings from \cite{Sommers} using  randomization inference techniques, which is essentially the same as the inference method proposed by CT. He found substantially larger p-values, ranging from 0.22 to 0.78, concluding that there is no evidence that the reform caused significant reductions in mortality. \cite{Kaestner} considered both randomization inference at the county and at the state levels.

Overall, standard errors clustered at the state level can indeed lead to large distortions when there is only one treated state, while standard errors clustered at the county level could lead to distortions due to spatial correlation (see \cite{ferman2019inference} for a thorough discussion on that). In contrast, a permutation test in the DID setting is essentially the same as the method proposed by CT, which is based on homoskedastic and independent assumptions on the errors. As we show, however,  spatial correlation can lead to substantial over-rejection if we consider CT at the county level. If we consider randomization inference at the state level, since we would have only one treated unit, spatial correlation would not generate problems if errors are strongly mixing (Corollary \ref{Cor_N1}). However, in both cases we may have distortions given that counties/states vary by population sizes, which would naturally lead to heteroskedasticity.

We revisit the use of the method proposed by CT in this empirical application in light of our findings. We consider a DID estimator based on OLS TWFE regression with no covariates. Point estimates are slightly different than reported in the original paper because we use the publicly available data set (so we do not have death counts for cells with fewer than 10 deaths), and because we weight observations by population mean across years to avoid problems with TWFE estimand highlighted in recent papers.\footnote{See, for example, \cite{Bacon}, \cite{chaisemartin2018twoway}, and \cite{Pedro}. We also restrict to counties with information non-missing information for all years. We end up with 485 control counties, as compared to 513 from the original study. } We find a p-value of $0.437$ when we use CT method at the county level, which is similar to the findings from \cite{Kaestner} based on permutations, who found a p-value of $0.384$.\footnote{\cite{Sommers} report a p-value of 0.04 for the OLS unadjusted regression using CRVE at the state level (see their Appendix Table 4). } Note, however, that the average population of the treated counties is around 280,000, while the average population of the control counties is around 91,000. Therefore, following FP, we should expect these p-values to be very conservative, as the method proposed by CT would recover a distribution with larger variance relative to the distribution of the DID estimator.  Correcting for that using the method proposed by FP, we find a p-value of $0.148$, which is much lower than the one using CT method.

However, the FP p-value above does not take spatial correlation into account. In this setting, if we assume errors are independent across states, then considering an aggregate  model at the state level would solve the problem. In this case, we find a p-value of 0.49 from CT and 0.42 from FP.\footnote{\cite{Kaestner} does not report p-value for the state-level OLS unadjusted regression. However, he finds a p-value of 0.404 for a state-level generalized linear model using a negative binomial distribution and log link, which is close to the p-value we found.  } Comparing the p-values from FP at the county and state levels suggests that spatial correlation is a significant problem when we consider the county-level model. In Appendix  \ref{test} we propose a test for spatial correlation, which reinforces  that  spatial correlation is important in this setting.

 When we consider the effects on amenable mortality rates, the differences in p-values become even more striking. We find a p-value of $0.314$ for CT at the county level, which goes down to $0.032$ when we correct for heteroskedasticity using FP. Considering FP at the state level, however, the p-value is up again to $0.238$. Overall, these results suggest that spatial correlation within states is important and, after taking that into account, it is not possible to reject the null that the reform had no effect on mortality.

Finally, we also consider the conservative test 2, proposed in Section \ref{Conservative2}. While in this setting it would be possible to solve the problem by aggregating the data at the state level, our method would remain feasible if we had treated counties from different states, or even in applications in which there is no distance metric across units. In this case, we find a p-value of $0.381$ ($0.200$ considering amenable mortality rates), which are remarkably close to the FP p-values at the state level. This reinforces that the conservative test 2 would have worked well even if we did not have information on a distance metric, and despite the fact that we have spatial correlation across counties. Therefore, this test provides an interesting alternative, particularly for  more complex applications in which it would not be possible to solve the inference problem by aggregating the data at the state level.

\section{Conclusion}

We consider the problem of inference in DID models when there are few treated units and errors are spatially correlated. We first show that, when there is a single  treated unit, the inference methods proposed by CT and FP, which were  designed for settings with few treated and many control units, remain asymptotically valid when errors are strongly mixing. This extends the set of possible applications in which the tests proposed by CT and FP can be reliably used when there is only a single treated unit. However, these methods can lead to over-rejection with more than one treated unit. We propose two  alternative inference methods that are asymptotically valid, though generally conservative, in the presence of spatial correlation. These tests provide  interesting alternatives when spatial correlation is likely relevant, but the researcher does not have information on a distance metric or is not willing to impose a structure on the spatial correlation.

\pagebreak