Inference with few treated units

Inference with few treated units

Keywords: Causal Inference, Small Sample Inference, Few Clusters, Synthetic Control, Difference-in-Differences, Randomization Inference

Introduction

In many causal inference applications, only one or a few units (or clusters of units) are treated. Examples include comparative case studies based on aggregate panel data, difference-in-differences (DiD) designs with few treated clusters, and randomized controlled trials (RCTs) with expensive treatments, among others. A key challenge in such settings is that standard inference methods based on asymptotic theory may be unreliable, even when the total number of units is large. For example, in simple treatment-control comparisons with few treated units (or clusters of units), the treated observations have a high leverage, which can lead to severe downward biases in conventional robust and cluster-robust standard errors Chesher,Carter,young2016improved,young_QJE,MacKinnonStata.

This survey aims to review and categorize the fast-growing literature proposing inference methods that are specifically designed to accommodate settings with few treated units. We consider both cross-sectional and panel data methods. We discuss the main assumptions that different inference methods rely on, explain the rationale behind them, and discuss their theoretical properties. We also emphasize connections and trade-offs between different approaches. As byproducts, we provide rigorous foundations for some heuristic approaches proposed in the literature, and we suggest slight modifications to improve the finite-sample performance of existing methods.

The existing inference methods can be first categorized into two main groups based on the source of uncertainty they consider: model-based and design-based inference methods. Model-based methods focus on the uncertainty coming from sampling the potential outcomes from an infinite super-population. Design-based methods focus on the uncertainty coming from the randomness in the treatment assignment.\footnote{ Our distinction between model-based and design-based approaches is consistent with usage of these terms in early statistical literature Sarndal1978, as well as more recent discussion in Econometrics baker2025. The term “sampling-based” is also used to cover our notion of “model-based” in some settings Abadie_finitepop. We can also consider methods that account for uncertainty in both potential outcomes and the treatment assignment Abadie_finitepop,Abadie2022. Since this approach also requires knowledge on the treatment assignment mechanism, we discuss these methods along with design-based approaches in Section (ref).}$^,$\footnote{We only cover frequentist inference, as this is more common in econometrics. However, we should note that some recent papers have suggested Bayesian approaches for inference in settings with few treated units Pang2022,benmichael2022,martinez2024. }

The main focus of this review is on model-based inference methods, as these methods are more prevalent in econometrics. We organize the presentation of the different model-based inference methods according to data availability. We first consider methods that are valid even in the extreme case in which there is only one treated unit and one treated period. The main challenge in this setting is that there is very limited information on the distribution of potential outcomes of the treated unit in the treated period. Therefore, the solutions proposed in the literature rely on extrapolating cross-sectional information from the untreated units or time-series information from untreated periods to assess uncertainty about the distribution of the treated potential outcomes. The choice between these alternatives should depend on the data availability and on the assumptions one is willing to make regarding the distribution of the errors in modeling the potential outcomes. Methods that rely on cross-sectional information generally allow for unrestricted time-series dependence, but rely on assumptions that restrict how the distribution of errors varies across treated and control units and often require a large number of control units. In contrast, methods that rely on time-series information allow for cross-sectional dependence and heteroskedasticity in the errors, but instead impose restrictions on the time-series properties of the errors (such as stationarity and weak dependence) and typically require a large number of untreated periods.

Given the limited information on the potential outcomes of the treated unit, methods that are valid even with one treated unit and period often restrict treatment effect heterogeneity either by assuming homogeneous treatment effects or by considering sensitivity analysis for the degree of treatment effects heterogeneity. We discuss how these assumptions relate to {a} sharp null hypothesis of no effect whatsoever and how we can reinterpret results that are valid under treatment effect homogeneity as inference on the realized treatment effect or in terms of prediction intervals (instead of confidence intervals).\footnote{In model-based settings, in which potential outcomes are treated as stochastic, we can think of a sharp null as $\Pr(Y(1) = Y(0))=1$, where $Y(1)$ and $Y(0)$ are potential outcomes. }

We then consider alternatives that are valid when there are few (but more than one) treated units, when there are many treated periods, or when there are many individual-level observations within each cluster. We discuss how the additional available information in each setting can be used to relax some of the strong assumptions required for valid inference in the extreme case with one treated unit and period. Nevertheless, these methods generally require alternative assumptions in other dimensions, and still typically require stronger assumptions than standard methods that are valid with many treated and many control units. For example, with multiple treated units, there are methods that can relax assumptions limiting heteroskedasticity and treatment effect heterogeneity while instead relying on symmetry assumptions on the errors and treatment effects. Moreover, while these methods control size with a fixed number of treated units, they may have low power when the number of treated units is very small. Therefore, methods designed specifically for settings with a single treated unit and period might be preferred due to power considerations when the number of treated units is very small (but greater than one).

Finally, we also discuss design-based approaches for inference. The theoretical justifications for design-based approaches are conceptually different from those for model-based approaches. Moreover, the focus is typically on different target parameters. This makes it difficult to directly compare them to model-based approaches. Nevertheless, we argue that design-based methods that are valid with few treated units have similar limitations as their model-based counterparts. We also argue that, even when the treatment is randomly assigned, settings with few treated units are prone to imbalances between treated and controls in the realized treatment assignment. Therefore, the longstanding debate on whether one should focus on unconditional inference or on inference conditional on imbalances between treated and controls is particularly relevant in these settings.

In categorizing different approaches for inference in settings with few treated, we also show asymptotic equivalence between some of the methods considered in the literature. In some cases, this provides a theoretical justification for methods that have been only heuristically justified. For example, we show that, in some settings, a wild bootstrap with the null imposed is asymptotically equivalent to an approximate randomization test based on sign-changes when the number of control units goes to infinity (keeping the number of treated units fixed).\footnote{Canay_wild_bootstrap derive conditions under which wild-cluster bootstrap is valid with a finite number of clusters. However, as we discuss in Section (ref), their results do not directly apply to treatment-control comparisons or DiD applications in which we use the wild-cluster bootstrap at the unit level.} Moreover, we propose slight modifications to improve the finite-sample performance some methods. The formal details are in the appendix.

This survey complements the excellent existing surveys and books on causal panel data methods where settings with few treated units are ubiquitous abadie2021using,dechaisemartin2022twoway,roth2023what,arkhangelsky2024causal,dechaisemartin2024book and the surveys and guides to practice regarding inference in regression models with clusters cameron2015practitioner,mackinnon2023cluster. Different from these references, we focus on the inferential challenges arising from the presence of few treated units (or clusters of units) across a wide variety of settings and methods. This allows us to highlight connections and common principles and, as a byproduct, yields some new formal justifications for existing heuristic procedures and variants of methods with improved finite sample performance. Many of the ideas and methods we discuss are generic and can be used in conjunction with a wide variety of panel and cross-sectional causal inference methods. We organize the literature based on conceptual aspects, such as the type of uncertainty, and practical aspects, such as the number of periods and treated units, to make the survey easy to navigate and useful for practitioners across a variety of fields and theoretical researchers alike.

In Section (ref) we start with a simple example that illustrates why standard methods that are asymptotically valid with many treated and many control units typically fail in settings with few treated units. This example also clarifies that the challenges for inference we discuss arise when the number of treated units is small and not necessarily when the share of treated units is small. Sections (ref) and (ref) present the model-based setting and inference methods. We then discuss in Section (ref) the use of design-based methods in settings with few treated units. Section (ref) concludes.

Why standard methods fail with few treated units: A simple example

We consider a simple example that illustrates why standard inference methods typically fail in settings with few treated units and the main challenges for inference in such settings. As we will see below, the insights from this simple example generalize to many other causal inference problems, including the estimation of treatment effects using regression and matching methods and panel data methods such as DiD, factor and interactive fixed effects models, and synthetic control.

We are interested in estimating the effect of a binary treatment $D_j$ on an outcome $Y_j$. We consider the regression model, $$ Y_j=\mu+\tau D_j+\eta_j,\quad \mathbb{E}[\eta_j | D_j]=0, $$ where $\tau$ is the treatment effect of interest.\footnote{To simplify the exposition, we interpret $\tau$ as the population OLS parameter of a regression of $Y_j$ on $D_j$ and a constant, and we do not define treatment effects in terms of potential outcomes, unlike in the following sections.} There are $N_1$ treated and $N_0$ control units. The outcomes of the treated (control) units are independently drawn from the distribution of $Y_j|D_j=1$ ($Y_j|D_j=0$).\footnote{See, for example, Assumption 3$'$ in AI_2006 for a similar setting.} To illustrate the issues with few treated units, we can consider an asymptotic regime where $N_1$ is fixed and $N_0\rightarrow \infty$. The OLS estimator of $\tau$ is equal to the difference between the treated and control mean,

where $\mathcal{J}_d\equiv \{j:D_j=d\}$ for $d=0,1$. Defining $\sigma^2_d = \mathbb{V}(\eta_j|D_j=d)$ for $d=0,1$, the variance of $\hat\tau$ is

We allow for heteroskedasticity, that is, $\sigma^2_0$ can be different from $ \sigma^2_1$.\footnote{As we discuss in detail in Section (ref), in a potential outcomes setting, heteroskedasticity can arise due to differences in the distribution of the untreated potential outcomes between the treated and control group and/or due to treatment effect heterogeneity.} Equation (ref) reflects the fact that we are estimating two means (for the treated and control group) to construct an estimator for $\tau$, so $\mathbb{V}(\hat \tau)$ is the sum of the variances of the treated and control averages.

The standard approach for making inferences on $\tau$ is to use a $t$-test based on the heteroskedasticity-robust variance estimator,

where $\hat \sigma^2_a = \frac{1}{N_a} \sum_{j \in \mathcal{J}_a} \hat \eta^2_j$, for $a \in \{0,1\}$, and $\hat \eta_j$ is the OLS residual. A key feature of $\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny heterosk}}$ is that it only uses information from the treated units to estimate $\sigma^2_1$ and information from the control units to estimate $\sigma^2_0$.

Under standard regularity conditions, as $\min\{N_1,N_0\} \rightarrow \infty$ (that is, we have many treated and many control observations), the OLS estimator $\hat \tau$ is asymptotically normal and

This result implies that we can conduct asymptotically valid inference for the null $H_0: \tau =0$ using the $t$-statistic $t_{\mbox{\tiny heterosk}}=\frac{\hat \tau}{\sqrt{\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny heterosk}}}}$ and standard normal critical values. We do not need to impose any assumptions on the relative magnitude of $\sigma_1^2$ and $\sigma_0^2$. Intuitively, this is because there is enough information to separately estimate $\sigma_1^2$ based on the treated units and $\sigma_0^2$ based on the control units when $N_1$ and $N_0$ are both large.

Consider now a setting in which $\min\{N_1,N_0\} \not \rightarrow \infty$. As an extreme case, suppose that $N_1 = 1$, while $N_0 \rightarrow \infty$. The first thing to notice in this case is that, while $\hat \tau$ remains unbiased, it will not be consistent and $\sqrt{N}$-asymptotically normal (even though $N=N_1+N_0 \rightarrow \infty$). In addition, the coefficient on the treatment dummy in this case would be such that the residual of the treated unit equals exactly zero. Therefore, $\hat \sigma^2_1=0$, implying $\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny heterosk}} = \frac{\hat \sigma^2_0}{N_0}$. Since the true variance of $\hat \tau$ is $\mathbb{V}(\hat \tau) =\sigma^2_1 + \frac{\sigma^2_0}{N_0}$, this means that $\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny heterosk}}$ severely underestimates $\mathbb{V}(\hat \tau)$. If $\sigma^2_0 = \sigma^2_1$, then the estimator for the variance of $\hat \tau$ is approximately $N$ times smaller than the true variance of $\hat \tau$. Intuitively, since there is only one treated observation, there is only enough information to estimate $\mathbb{E}[Y_j|D_j=1]$ (and, therefore, construct an estimator for $\tau$), but there is no information left to estimate $\sigma_1^2$ using only information from treated units.

In settings with more than one (but few) treated units, this problem is less severe, but the standard normal approximation for $t_{\mbox{\tiny heterosk}}$ can still be very inaccurate. Table (ref) provides an illustration based on a simple simulation study in which $Y_j \sim N(0,1)$ for all $j$. Note that this is a relatively favorable setting in that the estimator $\hat \tau$ is normally distributed even in finite samples. Moreover, we are in a setting in which errors are homoskedastic so that $\sigma_1^2=\sigma_0^2$. We set $N_0=100$ and show results for $N_1\in \{1,\dots,5\}$. The nominal level is 5%. While the problem is substantially more severe with $N_1 =1$, it is still relevant when $N_1>1$ but small. For example, with $N_1=5$, the true variance is still 32% larger than the expected value of the heteroskedasticity-robust variance estimator, which leads to a rejection rate of 15% instead of 5%. We stress that these simulations are only presented to illustrate the problem, and that these numbers could be different in other settings. In particular, there is no guarantee that the over-rejection would, in general, be limited to $15\%$ when $N_1 = 5$.

Interestingly, under this particular DGP in which errors are normal and homoskedastic, using homoskedastic variance estimators leads to tests with approximately correct size, even when $N_1=1$. The homoskedastic variance estimator is

where $\hat \sigma^2 = \frac{1}{N} \sum_{j=1}^N \hat \eta^2_j$. A notable difference relative to $\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny heterosk}}$ is that $\widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny homosk}}$ uses information from both the treated and the control units to estimate the variance of the error (which here is assumed to be the same for treated and control units). Under those simplifying assumptions, we have $\hat \tau \sim N\left( \tau, \sigma^2 + \frac{ \sigma^2}{N_0} \right)$, so $\frac{\hat \tau - \tau}{\sqrt{ \widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny homosk}}}} \overset{d}{\rightarrow} N(0,1)$ when $N_0 \rightarrow \infty$, even when $N_1=1$. Since $\sigma_0^2=\sigma_1^2$, we can use the control units to consistently estimate $\sigma_1^2$. In other words, we are extrapolating information from the control units to learn something about the distribution of the errors of the treated units. We emphasize that using $ \widehat{\mathbb{V}(\hat \tau)}_{\mbox{\tiny homosk}}$ does not yield valid inferences when the errors are heteroskedastic (even maintaining the normality assumption), because such extrapolation is not valid with heteroskedasticity. As we will see in Section (ref), while the assumption that $\eta_j$ has the same distribution for treated and control observations plays an important role, normality assumptions are typically not required for valid inference when $N_1=1$. That said, normality assumptions would be necessary for inference based on a $t$-test with homoskedastic standard errors.

Model-based setting

Notation and sources of uncertainty

We consider settings in which we observe $N$ units, indexed by $j=1,\dots,N$, over $T$ periods, indexed by $t=1,\dots,T$. We discuss cross-sectional settings where $T=1$ and panel data settings where $T>1$. Our goal is to identify and estimate the causal effect of a binary treatment $D_{j,t}$ on an outcome of interest $Y_{j,t}$. For ease of exposition, we focus on settings in which all treated units adopt the treatment in the same period and remain treated afterwards, so that we can define the treated and untreated potential outcomes as $Y_{j,t}(1)$ and $Y_{j,t}(0)$, respectively. However, the basic principles we discuss in the survey apply more generally. The observed outcome is related to the potential outcomes as $Y_{j,t}=D_{j,t}Y_{j,t}(1)+(1-D_{j,t})Y_{j,t}(0)$. We may also observe a vector of covariates $X_{j,t}$. The causal effect of $D_{j,t}$ for unit $j$ in period $t$ is $\tau_{j,t}=Y_{j,t}(1)-Y_{j,t}(0)$, where $\tau_{j,t}$ could be random or fixed. In the following, we suppress the index $t$ whenever we describe cross-sectional methods.

We first consider a model-based setting in which treatment assignment is fixed (or conditioned on), and the stochastic variation comes from the potential outcomes, $(Y_{j,t}(0),Y_{j,t}(1))$. Since treatment assignment is fixed and we are considering settings with no variation in treatment timing, we can define $\mathcal{J}_1$ ($\mathcal{J}_0$) as the set of treated (control) units, and $\mathcal{T}_1$ ($\mathcal{T}_0$) as the set of post- (pre-)treatment periods. We also define $N_1=|\mathcal{J}_1|$ ($N_0=|\mathcal{J}_0|$) as the number of treated (control) observations, and $T_1 = |\mathcal{T}_1|$ ($T_0=|\mathcal{T}_0|$) as the number of post- (pre-)treatment periods.

This model-based perspective is appropriate in applications where there is a well-defined large population from which the sample is drawn.\footnote{See, e.g., roth2023what for a discussion in the context of DiD methods.} The standard justification for this interpretation is that there is a large population of units, and we are sampling a negligible fraction of it. Alternatively, we can interpret this large population as a super-population representing different possible realizations of random variables that determine the outcomes of the units in the sample.

This alternative interpretation of sampling is not new: it dates back to at least since Haavelmo1944, who wrote:

and also appears in more recent discussions on the nature of uncertainty in structural analyses heckman2000causal,heckman2005scientific and in Biostatistics hernan2010causal.

This idea of considering a super-population representing different possible realization of potential outcomes (a “multiverse”) is actually prevalent in popular culture, and has been present in movies such as “Spider-Man: Into the Spider-Verse (2018),” “Doctor Strange in the Multiverse of Madness (2022),” “Spider-Man: Across the Spider-Verse (2023),” “The Flash (2023),” and “Deadpool & Wolverine (2024).”

This alternative interpretation of sampling is particularly relevant for justifying model-based approaches in applications with few treated units. For example, we often consider applications where a few states are treated, and we observe outcomes for all states. Moreover, we can even consider settings in which treatment would not be well-defined for the control units. For example, in the study of the economic impact of the German reunification by abadie2015comparative, it is difficult to imagine other countries reunifying abadie2021using. Still, in such setting, we can consider uncertainty coming from the fact that, conditional on the German reunification occurring in 1989, there were infinite possible realizations of the random variables that determine the potential outcomes of the treated and control countries, from which we observe a single realization.

Finally, model-based analyses can be interpreted as being conditional on the treatment assignment if the treatment assignment is stochastic. This interpretation allows for explicitly incorporating treatment assignment mechanisms into the analysis. More specifically, the characteristics of the treatment assignment mechanism are captured by the fact that we are considering the distribution of potential outcomes conditional on the treatment allocation. For example, if the units select into treatment based on some unobservable (to the econometrician) factors, then the distribution of these unobservables, conditional on the treatment allocation, depends on the treatment status Ferman_JASA.

Model, parameters of interest, and estimators

The literature typically focuses on treatment effects on the treated units in the treated period. A popular target parameter is the average treatment effect on the treated units (ATT),

where this expectation is taken over the distribution of the potential outcomes in the population (or super-population). Most of our discussions remain relevant if we consider alternative target parameters, such as the sequence $\left\{\mathbb{E}\left[\frac{1}{N_1}\sum_{j \in \mathcal{J}_1} \tau_{j,t} \right] \right\}_{t \in \mathcal{T}_1}$.\footnote{The methods discussed in Section (ref) are an exception in that they would not allow for inference on the sequence of treatment effects. } We discuss other alternatives, such as inference on the realized treatment effects, in Section (ref).

The key challenge when identifying and estimating parameters like the ATT is that the untreated potential outcomes are not observed for the treated units in the post-treatment period. By the definition of the potential outcomes, we can write the treatment effect $\tau_{j,t}$ for $j\in \mathcal{J}_1$ and $t\in \mathcal{T}_1$ as $$ \tau_{j,t}=Y_{j,t}(1)-Y_{j,t}(0)=Y_{j,t}-Y_{j,t}(0). $$ Thus, the key unknowns are the untreated potential outcomes for the treated units in the post-treatment period.

Most of the model-based methods we analyze in this review postulate models for the potential outcomes that can be written as

where $M_{j,t}$ is a mean-predictor of $Y_{j,t}(0)$ and $\epsilon_{j,t}$ is an error term, which is typically assumed to have mean zero. Given a model $M_{j,t}$, the estimator of the ATT will usually take the form

where $\hat M$ is an estimator for $\frac{1}{N_1} \frac{1}{T_1} \sum_{j \in \mathcal{J}_1} \sum_{t \in \mathcal{T}_1} M_{j,t}$. For methods in which we estimate treatment effects by imputing $\hat M_{j,t}$ as the unobserved counterfactual for unit $j$ at period $t$, the estimator takes the form $\hat M = \frac{1}{N_1} \frac{1}{T_1} \sum_{j \in \mathcal{J}_1} \sum_{t \in \mathcal{T}_1} \hat M_{j,t}$.

Given Equation (ref), the estimation error $\hat{\tau}-\tau^\ast$ can be decomposed into three terms,

The first term is the average of the errors in $Y_{j,t}(0)$ for the treated units in the post-treatment periods. The second term captures the heterogeneity in the treatment effects. Here we distinguish between two kinds of heterogeneity in the treatment effects, which are not mutually exclusive. We call stochastic heterogeneous treatment effects as the possibility that, for a given $j$ and $t$, $\tau_{j,t}$ is stochastic. This would be the case, for example, in settings in which the treatment has an effect on the variance of $Y_{j,t}(1)$. We call deterministic heterogeneous treatment effects when $\mathbb{E}[\tau_{j,t}]$ varies with $j$ and/or $t$. The third term captures the estimation error of the counterfactual.

Note that the expected value of the second term is zero by definition. We focus on cases in which the expected value of the first and third terms is also zero (or at least asymptotes to zero when $N_0$ and/or $T_0$ increases). Under these assumptions, the estimator is asymptotically unbiased. The main challenge for inference in settings with few treated units is that we generally cannot apply a Central Limit Theorem (CLT) to the first and second term.\footnote{Exceptions would be if $T_1 \rightarrow \infty$ or if we consider the case with many individual-level observations within cluster, and we impose weak dependence on the individual-level errors. We discuss those settings in Section (ref).} Moreover, it is not possible to consistently estimate the variances of $\epsilon_{j,t}$ and $\tau_{j,t}$ using only information from the treated units. As a result, standard inference methods that rely on asymptotic normality of the estimator and on normal approximations of studentized test statistics are not reliable with few treated units.

Examples

Our framework encompasses popular methods in cross-sectional analyses. For example:

Comparison of means/RCT. In this example, $M_{j}$ is equal to the mean of the untreated potential outcome, $M_{j}=\mathbb{E}[Y_{j}(0)]$. $M_{j}$ can be estimated as $\hat M_j = \frac{1}{N_0} \sum_{j \in \mathcal{J}_0} Y_j$.

Regression. Let $M_{j}=\mu_0+X_{j}'\beta_0$, for some observed variables $X_j$. In this case, $\beta_{0}$ is the population regression parameter from a regression of $Y_{j}$ on $X_j$ in the subpopulation with $D_j=0$. $M_{j}$ can be estimated using the sample analog of this regression as $\hat{M_{j}}=\hat\mu_0+X_{j}'\hat\beta_0$, where $(\hat{\mu}_0,\hat\beta_0)$ are obtained from a regression of $Y_{j}$ on $X_j$ in the subpopulation with $D_j=0$.\footnote{ Another way to implement an estimator for the ATT in this case would be to run a linear regression of $Y_j$ on $D_j$, $X_j$, and interactions of $D_j$ and $X_j - \bar X$ ImbensWoodlrdige2009. If we run a linear regression of $Y_j$ on $D_j$ and $X_j$, then the estimator would not necessarily have the format from Equation (ref). Still, the main takeaways regarding inference with few treated would remain valid in this case (see Section (ref) of Appendix (ref) for an example). }

Matching. Let $M_j= \mu_0(X_j)$, where $\mu_0(X_j)=\mathbb{E}[Y_{j}(0)|X_j]$. In this case, the $K-$nearest neighbor matching estimator would take the form $\hat M_j = \frac{1}{K} \sum_{q \in \mathcal{V}_j} Y_q$, where $\mathcal{V}_j$ is the set with the $K$ control units with value of $X_q$ closest to $X_j$. Other types of matching estimators could also be used Imbens_matching.

The framework also encompasses many popular panel data methods. For example:

Difference-in-differences. DiD methods are typically motivated by a two-way fixed effects model for the untreated potential outcome, $M_{j,t} = \lambda_t + \mu_j$, where $\lambda_t$ is a time fixed effect and $\mu_j$ is an unit fixed effect. If all treated units start treatment at the same period, a canonical DiD estimator is $\hat \tau = \frac{1}{N_1} \sum_{j \in \mathcal{J}_1} \left( \bar Y_{j,post} - \bar Y_{j,pre} \right) - \frac{1}{N_0} \sum_{j \in \mathcal{J}_0} \left( \bar Y_{j,post} - \bar Y_{j,pre} \right)$, where $\bar Y_{j,post} = \frac{1}{T_1} \sum_{t \in \mathcal{T}_1}Y_{j,t}$ and $\bar Y_{j,pre} = \frac{1}{T_0} \sum_{t \in \mathcal{T}_0}Y_{j,t}$. In this case, we have that $$\hat M_{j,t}=\bar Y_{j,pre} + \frac{1}{N_0} \sum_{j' \in \mathcal{J}_0} \left( Y_{j',t} - \bar Y_{j',pre} \right).$$ See dechaisemartin2022twoway and roth2023what for alternative DiD estimators and settings with more complicated adoption patterns.

Factor and interactive fixed effects approaches. The interactive fixed effects model is $M_{j,t} = X_{j,t}\beta+ \lambda_t '\mu_j$, where $\lambda_t$ is a vector of time-varying factors and $\mu_j$ is a vector of unit-specific loadings. It nests the standard factor model $M_{j,t} = \lambda_t '\mu_j$ as a special case. A natural estimator of $M_{j,t}$ is $\hat{M}_{j,t}=X_{j,t}\hat\beta+ \hat\lambda_t '\hat\mu_j$, where $(\hat\beta, \hat\lambda_t, \hat\mu_j)$ are estimated based on untreated periods and untreated groups gobillon2016regional,xu2017generalized. Alternatively, models with factor structures can be estimated using matrix completion techniques amjad2018robust,athey2021matrix.

Synthetic control. This method has been considered under different assumptions on $M_{j,t}$, such as linear factor models, low-rank matrices, or autoregresive models. Synthetic control estimators of $M_{j,t}$ typically take the form $\hat{M}_{j,t}=\sum_{j\in \mathcal{J}_0}\hat{w}_jY_{j,t}$, where the weights $\{\hat{w}_j\}_{j\in \mathcal{J}_0}$ are obtained based on the pre-treatment period abadie2021using.

Model-based inference methods

Methods that are valid even with $N_1=1$ and $T_1=1$

We consider first inference methods that are valid even for an extreme case in which $N_1=1$ and $T_1=1$. Therefore, in this section, we let $j=1$ be the treated unit, and $t=T$ be the post-treatment period. Importantly, though, these methods are usually also valid for settings with $N_1>1$ and/or $T_1 > 1$.

A distinctive feature of this extreme setting with $N_1=T_1=1$ is that there is only a single observation $Y_{j,t}$ that is treated. Therefore, we have one observation to estimate the treatment effect $\tau^\ast$, but not enough variation to quantify the uncertainty about the distribution of $Y_{1,T}(1)$ (which, following the decomposition in Equation (ref), encompasses the uncertainty on $\epsilon_{1,T}$ and $\tau_{1,T}$). The solutions that are valid even in such extreme settings then attempt to use information from the control units and/or the pre-treatment periods to learn about the distribution of $Y_{1,T}(1)$. Therefore, in all cases, we need to impose assumptions that are stronger than those required by standard inference methods for settings with many treated and many controls units. For example, since the untreated units do not provide any information on the distribution of treatment effects, restrictions on treatment effect heterogeneity are typically unavoidable.

In Sections (ref) and (ref), we consider settings with no stochastic treatment effect heterogeneity in which $\tau_{1,T}$ is treated as a fixed parameter. Therefore, in these cases we have that the ATT is given by $\tau^\ast = \tau_{1,T}$. In Section (ref), we focus on methods that exploit the cross-section variation for inference, while in Section (ref) we focus on methods that exploit time series variation. Then, in Section (ref), we discuss alternative interpretations and possible ways to relax the homogeneous treatment effects assumption.

Methods that exploit cross-sectional variation

conley20211inference provide a leading example of a model-based method that exploits variation in the cross-section. They propose an inference method for DiD settings that is valid when $N_1$ is fixed (including the extreme case where $N_1=1$) and $N_0 \rightarrow \infty$. Their main method is valid when (i) $\{\epsilon_{j,t}\}_{t \in \mathcal{T}}$ is iid across $j$ and (ii) the treatment effects are homogeneous.\footnote{conley20211inference also consider alternatives that relax both of these assumptions.} We discuss the interpretation of their method with stochastic heterogeneous treatment effects in Section (ref).

With $N_1=T_1=1$, the ATT from Equation (ref) is the (homogeneous) treatment effect for the treated unit in the post-treatment period, $\tau^\ast = \tau_{1,T}$, and the DiD estimator is given by $\hat \tau = \left[ Y_{1,T} - \bar Y_{1,pre}\right] - \frac{1}{N_0} \sum_{j \in \mathcal{J}_0}\left[Y_{j,T} - \bar Y_{j,pre} \right]$. Under standard regularity conditions, as $N_0 \rightarrow \infty$, we have that $\hat \tau = \tau^\ast + W_1 + \frac{1}{N_0} \sum_{j \in \mathcal{J}_0} W_j \overset{p}{\rightarrow} \tau^\ast + W_1$, where $W_j \equiv \epsilon_{j,T} - \frac{1}{T_0} \sum_{t \in \mathcal{T}_0} \epsilon_{j,t}$ is a linear combination of the errors. Therefore, under standard DiD identification assumptions ($\mathbb{E}[W_j] = 0$ for all $j$), the DiD estimator in this setting with one treated unit would be unbiased, but it would not be consistent. Moreover, the asymptotic distribution of $\hat \tau$ when $N_0 \rightarrow \infty$ will depend only on the errors of the treated unit, $W_1$.

The main intuition underlying conley20211inference's method is that the residuals of the controls $\{\widehat W_j \}_{j \in \mathcal{J}_0}$ asymptotically recover the distribution of $W_j$ in the control group. Then, under the assumption that $W_1$ has the same distribution as $W_j$ for $j \in \mathcal{J}_0$, this implies that we also recover the distribution of $W_1$. Therefore, we can construct confidence intervals for $\tau^\ast$ that are asymptotically valid when $N_0 \rightarrow \infty$ using the quantiles of the distribution of $\{\widehat W_j \}_{j \in \mathcal{J}_0}$, $I = [\hat \tau - \hat Q_W(1-\gamma/2), \hat \tau - \hat Q_W(\gamma/2)$], where $\hat Q_W(u)$ is the $u$ empirical quantile of $\{\widehat W_j \}_{j \in \mathcal{J}_0}$. Likewise, we can construct a p-value

which is asymptotically valid when $N_0 \rightarrow \infty$ for a two-sided test that $H_0: \tau^\ast = c$.

Crucially, since we rely on information on the control units to learn about the distribution of the errors of the treated unit, we need to rely on assumptions that link these two distributions. conley20211inference in their standard implementation assume that $W_1$ and $W_j$ for $j \in \mathcal{J}_0$ have the same distribution. Moreover, we need a large number of control units to consistently estimate the distribution of $W_j$ for the controls. This is achieved in their setting under the assumption that $W_j$ is iid across $j$, and considering an asymptotic approximation with $N_0 \rightarrow \infty$.

Interestingly, by exploiting cross-section variation for inference, conley20211inference do not have to restrict the time series properties of $\epsilon_{j,t}$. The main idea is that, to derive the distribution of $\hat \tau$ in this setting, we can essentially collapse the data into differences between post- and pre-treatment periods. Then this post- and pre-differences in the errors, $W_j$, already incorporates any serial correlation in $\epsilon_{j,t}$ that is relevant to derive the distribution of $\hat \tau$. Moreover, since this collapsed data is essentially a comparison of means, this approach can also be used as a model-based inference approach in RCTs with few treated units.

The method proposed by conley20211inference can also be used in settings with $N_1>1$ and/or $T_1>1$. In this case, we still cannot allow for stochastic treatment effects heterogeneity, but we can allow for deterministic treatment effects heterogeneity. Therefore, we continue to treat $\tau_{j,t}$ as deterministic parameters, but we allow them to vary with both $j$ and $t$, so that $\tau^\ast = \frac{1}{N_1} \frac{1}{T_1} \sum_{j \in \mathcal{J}_1}\sum_{t \in \mathcal{T}_1} \tau_{j, t}$ (see Section (ref) for more details on the distinction between these two types of treatment effects heterogeneity). Note that with $N_1>1$ and/or $T_1>1$ we have that $\hat \tau \overset{p}{\rightarrow} \tau^\ast + \frac{1}{N_1}\sum_{j \in \mathcal{J}_1} \dot W_j$, where, in this case, $\dot W_j \equiv \frac{1}{T_1} \sum_{t \in \mathcal{T}_1} \epsilon_{j,t} - \frac{1}{T_0} \sum_{t \in \mathcal{T}_0} \epsilon_{j,t}$. Therefore, the asymptotic distribution of $\hat \tau$ depends only on the errors $\dot W_j$ of the treated, regardless of whether $\tau_{j,t}$ varies with $j$ or $t$. Since, under the assumptions considered above, we can recover the distribution of $\dot W_j$ for the treated units by extrapolating from information from the control units, the inference method is valid in this case, even with a small number of treated units.\footnote{With $N_1>1$, it is also possible to consider an alternative using a studentized test statistic, as proposed by MACKINNON2020435. While with $N_1$ fixed it would not be possible to guarantee that the test is valid when there is treatment effect heterogeneity (even when $N_0 \rightarrow \infty$), this has important advantages in settings with $N_1 \rightarrow \infty$, as discussed in Appendix (ref) and Section (ref) . }

As conley20211inference noted, their inference procedure is asymptotically equivalent to a permutation test. We show in Appendix (ref) that a slightly modified version of their method is exactly equivalent to a permutation test, being valid with fixed $N_1$ and $N_0$ provided that $\{\epsilon_{j,t}\}_{t=1}^T$ is iid across $j$. The only additional assumption for validity with fixed $N_0$ is that we cannot have deterministic treatment effect heterogeneity in the cross-section when we have more than one treated unit (which is allowed asymptotically when $N_0 \rightarrow \infty$). We present details on that in Appendix (ref). The proposed alternative implementation of conley20211inference's method exhibits relevant advantages at essentially no cost. More specifically, it has the advantage of having finite-$N_0$ validity (under slightly stronger effect homogeneity assumptions), while being asymptotically equivalent to the original implementation when $N_0 \rightarrow \infty$. Therefore, we recommend the use of this alternative implementation. An important caveat of this approach when both $N_0$ and $N_1$ are small is that the reference distribution for the test will have a small number of support points. This implies that, in a setting with $N_1=1$, the p-value would always be, by construction, greater than $\frac{1}{N_0 + 1}$. Therefore, at least $N_0 \geq 10$ ($N_0 \geq 20$) would be required to reject the null at a 10% ($5\%$) significance level.

A number of papers propose similar approaches that exploit cross-sectional variation for inference with a single or few treated units under different sets of assumptions and for different settings. ferman2019inference consider a DiD setting in which $W_j$ may exhibit heteroskedasticity based on observed variables. As a leading example, when $Y_{j,t}$ represents state $\times$ time aggregates, we should expect $W_j$ for states with larger populations to have relatively lower variances. As a result, the standard method of conley20211inference over- (under-)rejects when the treated state is relatively small (large).\footnote{If we assume a setting in which treatment assignment is also stochastic, and all units have uniform probability of being the treated one, then the method proposed by conley20211inference would remain valid when we consider inference unconditional with respect to the treatment assignment. However, we would have size distortions for inference conditional on the populations of the treated states. See, e.g., ferman2019inference for some arguments in favor of considering conditional inference in these settings.} To address this issue, ferman2019inference propose to estimate this heteroskedasticity using control residuals. Then we can re-scale $\{\widehat W_j \}_{j \in \mathcal{J}_0}$ using the estimated heteroskedasticity and recover the distribution of $W_1$.\footnote{More specifically, ferman2019inference show that, under a wide range of structures for the within-state correlation, $\mathbb{V}(W_j)$ would have a parsimonious formula $A + B/X_j$, where $X_j$ is the population of unit $j$. Then one could run an OLS regression of $\widehat W_j^2$ on a constant and $1/X_j$ using the control states, and consider the re-scaled residuals $\widetilde W_j = \widehat W_j \sqrt{\frac{\widehat{\mathbb{V}(W_1)}}{\widehat{\mathbb{V}(W_j)}}}$ using the estimated parameters $\hat A$ and $\hat B$. By focusing on $W_j$, we bypass the problem of dealing with the serial correlation, and the fact that $\hat \epsilon_{j,t}$ is not a consistent estimator for $\epsilon_{j,t}$ when we have a panel with fixed $T$.} This approach is asymptotically valid in settings with $N_0 \rightarrow \infty$ and $N_1$ fixed (even when $N_1=1$) if we assume a scale-change model for the distribution of $W_j = \xi_j \times \sigma(X_j)$, where $\xi_j$ is iid across treated and control units, but $\mathbb{V}(W_j)$ depends on a set of observed covariates, $X_j$. In particular, in the example in which $Y_{j,t}$ represents state $\times$ time aggregates, this approach corrects for the fact that the errors of the treated state should have a higher variance than the errors of the control states, if the treated state is relatively smaller. Importantly, this approach still does not allow for {stochastic} heterogeneous treatment effects, and does not allow for heteroskedasticity in $W_j$ beyond the one that is estimable with the observed data. We also note that the method proposed by ferman2019inference relies crucially on $N_0 \rightarrow \infty$, which allows for consistent estimation of the heteroskedasticity function using the control units. As a result, the modifications in Appendix (ref) that ensure the finite $N_0$ validity of conley20211inference's method do not guarantee the finite $N_0$ validity of ferman2019inference's method.

Still considering DiD settings, alvarez2023extensions consider settings with variation in treatment timing, a topic that has been extensively studied in the recent DiD literature (see dechaisemartin2022twoway and roth2023what for surveys). They show that some recently proposed alternatives may also lead to substantial over-rejection when there are few treated units and extend the methods proposed by conley20211inference and ferman2019inference to accommodate variation in treatment timing. They also derive uniform confidence bands for dynamic DiD specifications that are valid with fixed $N_1$, including the case $N_1=1$). In another paper, alvarez2023inference relax the independence assumption across units. They show that the methods proposed by conley20211inference and ferman2019inference remain valid under weak dependence in the cross section when $N_1=1$, and propose alternatives when $N_1 > 1$.

The main idea underlying conley20211inference's method has also been used for inference in combination with other estimators than DiD. For example, synthetic_did (Algorithm 4) build on this idea to develop an approach for making inferences based on their Synthetic DiD estimator in applications with $N_1=1$, and ferman_matching develops a related approach for inference based on matching estimators with few treated and many control units.

There are also other alternative methods that are valid with $N_1=1$, even when $N_0$ is fixed. For example, Donald consider DiD settings with fixed $N_1$ and $N_0$, and derive the exact distribution of the $t$-statistic under normality and homoskedasticity assumptions on the errors (in contrast to conley20211inference, who do not require normality assumptions). Finally, hagemann2020inference proposes another alternative for settings with a finite number of heterogeneous clusters, where $N_1=1$. Each cluster is assumed to be large, so that the average error for each cluster is approximately Gaussian, and heteroskedasticity is allowed by imposing upper bounds on the ratio between the variance of the treated unit relative to the variance of the controls. In addition to dealing with heteroskedasticity, this solution can also deal with stochastic treatment effects heterogeneity (see Section (ref) for more details).

Methods that exploit time-series variation

Another alternative is to exploit the time-series variation. In such cases, we generally consider settings in which $T_0 \rightarrow \infty$, and we need to impose assumptions on the time-series of the errors (e.g., stationarity and weak dependence). Exploiting the time-series dimension allows for constructing inference methods that are typically valid with $N_0$ fixed. These methods typically only require models for the potential outcomes of the treated units, and allow for spatial correlation and richer heteroskedasticity in the cross-section. Therefore, these methods are complementary to the methods reviewed in Section (ref) in terms of the dimensions in which we need to impose strong assumptions and those in which we can allow for more flexibility.

As in Section (ref), we consider a setting with $N_1=T_1=1$ and no stochastic treatment effect heterogeneity, so that $\tau^\ast = \tau_{1,T}$ (we relax this assumption in Section (ref)). If the errors $\{\epsilon_{1,t}\}_{t=1}^T$ are stationary, then testing $H_0:\tau^\ast=0$ is akin to testing the null hypothesis of no structural break at the end of the sample. More specifically, under the null of no effect and stationarity of $\{\epsilon_{1,t}\}_{t=1}^T$, we have that $Y_{1,t} - M_{1,t} = \epsilon_{1,t}$ for $t=1,\dots,T$ is stationary. hahn2017synthetic suggest testing this implication using the end-of-sample instability test of andrews2003end in the context of the synthetic control method, while ferman2019inference propose to apply this test to make inferences in DiD applications with large $T_0$. The basic idea is to compare the post-treatment residual under the null $H_0:\tau^\ast=0$, $\hat\epsilon_{1,T}=Y_{1,T}-\hat{M}_{1,T}=Y_{1,T}(1)-\hat{M}_{1,T}=Y_{1,T}(0)-\hat{M}_{1,T}$, to the pre-treatment residuals $\{\hat\epsilon_{1,t}\}_{t\in \mathcal{T}_0}$ and reject if $\hat\epsilon_{1,T}$ is large relative to $\{\hat\epsilon_{1,t}\}_{t\in \mathcal{T}_0}$, where $\hat \epsilon_{1,t} = Y_{1,t} - \hat M_{1,t}$.\footnote{If we are interested in testing $H_0:\tau^\ast=c$, then the post-treatment residual (under the null) is $\hat\epsilon_{1,T}=Y_{1,T}-c-\hat{M}_{1,T}$.} For example, one can compute (two-sided) $p$-values for testing $H_0:\tau^\ast=0$ as

The $p$-value in (ref) is asymptotically valid under two conditions. First, since $\hat\epsilon_{1,t}-\epsilon_{1,t}=\hat{M}_{1,t}-M_{1,t}$, the estimation error $\hat{M}_{1,t}-M_{1,t}$ needs to be negligible so that the difference between the estimated errors $\hat\epsilon_{1,t}$ and the true errors $\epsilon_{1,t}$ is negligible. Second, the true errors need to be stationary and weakly dependent so that the (infeasible) $p$-value based on the true errors, $p=T^{-1}\sum_{t=1}^T\{|\epsilon_{1,t}|\ge |\epsilon_{1,T}|\}$, is valid. Importantly, however, because this approach only exploits the time-series dimension of the problem, it does not require assumptions on the cross-sectional heteroskedasticity.

Many modern approaches for estimating $M_{1,t}$ involve high-dimensional estimation problems. A leading example is the synthetic control method where researchers are estimating a weight for each control unit based on the pre-treatment data, so that there are $N_0$ parameters and $T_0$ data points. In many applications, $T_0$ is small or moderate, while $N_0$ is comparable to or even larger than $T_0$. In such applications, the estimation error $\hat{M}_{1,t}-M_{1,t}$ can be substantial and render inference methods based on asymptotic approximations inaccurate. To address this challenge, chernozhukov2021exact propose a conformal inference method that has a double-justification: it is valid in finite samples when the untreated potential outcomes are exchangeable in the time-series dimension (e.g., when the potential outcomes are iid across time), and it is asymptotically valid when the data exhibit dynamics and time-series dependence as $T_0\rightarrow \infty$.

The key idea of the conformal inference procedure is to estimate $\hat{M}_{1,t}$ using the data from all $T$ periods under the null hypothesis $H_0:\tau^\ast=0$. Estimation under the null hypothesis implies that $\{\hat\epsilon_{1,t}\}_{t=1}^T$ is exchangeable if the data are exchangeable, which implies that the method is exact in finite samples based on classical arguments for randomization tests hoeffding1952large,romano1990behavior, that is $P(\hat{p}\le \alpha)\le \alpha $ for any $(N,T)$. chernozhukov2021exact show that estimation under the null is also crucial for a good finite sample performance with dependent (non-exchangeable) data. Intuitively, this is because estimation under the null ensures that all residuals are affected equally by the estimation error. By contrast, if $M_{1,t}$ is estimated using only pre-treatment data, then $\{\hat\epsilon_{1,t}\}_{t=1}^{T_0}$ may be too small compared to $\hat\epsilon_{1,T}$ due to overfitting, which can lead to substantial size distortions. cattaneo2021prediction develop an alternative approach for dealing with the estimation error $\hat{M}_{1,T}-M_{1,T}$. Since their method is designed for constructing prediction intervals when the treatment effects are stochastic, we discuss it in Section (ref).

Finally, when exploiting the time-series dimension, non-stationary data are ubiquitous and lead to non-standard behavior of common estimators of $M_{j,t}$. To this end, masini2021counterfactual_jasa,masini2022counterfactual_jbes characterize the asymptotic properties of regression and Lasso estimators of $M_{1,t}$ when the data exhibit deterministic or stochastic trends. For settings with $T_1=1$ and $T_0\rightarrow \infty$, their theoretical results justify a residual-based inference method as discussed above when the errors are stationary.

Allowing for stochastic treatment effects heterogeneity

In Sections (ref) and (ref), we focused on methods that are valid under the assumption that there is no stochastic treatment effect heterogeneity. That is, the treatment effects $\tau_{j,t}$ are fixed parameters. This assumption is strong in many settings: it requires $Y_{j,t}(1)$ and $Y_{j,t}(0)$ only differ by a constant, so that the treatment only has a location-shift effect. In this section, we discuss alternative interpretations of the inference methods discussed in Sections (ref) and (ref) and approaches for assessing the sensitivity to the fixed treatment effects assumption.

\paragraph{Testing a sharp null.} An alternative interpretation of the methods discussion in Sections (ref) and (ref) is that they are testing the null hypothesis $H_0: \mathbb{P}[Y_{1,T}(1) = Y_{1,T}(0) + c]=1$. This is a sharp null that the treatment effect is homogeneous and equal to $c$ in the population. This has been considered by Chung2021 (see also Remark 4.11 of Bugni2018). The main idea is that, under this null, we have the homogeneity in treatment effects that is required by the methods revised in Sections (ref) and (ref). By setting $c=0$, note that testing this alternative null is similar in spirit to testing a sharp null of no effect whatsoever in design-based settings (see Section (ref)), though, in a model-based setting, “whatsoever” concerns the population (or super-population) from which the sample was drawn and not the sample at hand.

\paragraph{Inference on the realized treatment effects.} Another option is to allow for stochastic treatment effect heterogeneity and to make inferences conditional on the treatment effects. This corresponds to considering a repeated sampling framework in which $\tau_{j,t}$ is fixed, while uncertainty comes from different realizations of $\epsilon_{j,t}$ and other stochastic variables contained in $M_{j,t}$ (conditional on not only treatment assignment, but also on $\tau_{j,t}$ for $j \in \mathcal{J}_1$ and $t \in \mathcal{T}_1$). In this case, instead of considering a target parameter that is the expected treatment effect on the treated, we consider an alternative target parameter which is the realized treatment effect on the treated.

To make this distinction more concrete, consider an example from AF2025, in which the treatment effect may be $a>0$ if it rains in the region of the treated unit (which happens with probability $p$), and $0$ otherwise. In this case, the ATT is given by $\tau^\ast = p\times a$. With $N_1=T_1=1$, it would be essentially impossible to draw inference on $\tau^\ast$ without further assumptions, because we only observe $Y_{1,T}(1)$ either in a setting that rained, or in a setting that did not rain. However, we may still draw inference on the realized treatment effects (in this simple example, we would be conditioning on whether it rained or not).

Note that settings that allow deterministic treatment effect heterogeneity, but restrict stochastic treatment effect heterogeneity --- that is, $\tau_{j,t}$ is non-stochastic, but can vary across $j$ and $t$ --- can be seen as (implicitly) considering a model with stochastic treatment effect heterogeneity, but with the realized treatment effect as the target parameter.

For the methods in Sections (ref) and (ref) to be valid once we condition on the treatment effects, the underlying assumptions need to hold conditionally on the treatment effect. Thus, the key question then is whether those assumptions are reasonable once we condition on $\tau_{1,T}$. To illustrate, consider the method proposed by conley20211inference. Suppose that, unconditionally, $\{\epsilon_{j,t}\}_{t=1}^T$ is iid across $j$ (for both treated and controls). Then, a simple sufficient condition for the conditional validity of conley20211inference's method is independence of $\tau_{1,T}$ and $\{\epsilon_{j,t}\}_{t=1}^T$ for $j=1,\dots,N$. An alternative sufficient condition is that the distribution of $\epsilon_{1,t}$ changes conditional on $\tau_{1,T}$, but $\{\epsilon_{1,t}\}_{t=1}^T | \tau_{1,T} \overset{d}{=} \{\epsilon_{j,t}\}_{t=1}^T | \tau_{1,T}$ for all $j \in \mathcal{J}_0$. This alternative sufficient condition may be plausible if the randomness in the treatment effects is due to common shocks that affect all units. That is, if there are some unobserved shocks that affect not only $\tau_{1,T}$, but also $\{\epsilon_{j,t}\}_{t=1}^T$ for all units. By contrast, it would generally be implausible if the shocks that affect $\tau_{1,T}$ are specific to the treated unit. See AF2025 for further discussion on the validity of these methods for inference on the realized treatment effects.

\paragraph{Constructing prediction intervals.} An alternative to considering those relatively strong assumptions is to consider prediction intervals for treatment effects. A prediction interval for a scalar random variable $X$, defined in the same probability space as the sample, is a random set $\mathcal{I}_X$ such that $P[X \in \mathcal{I}_X] \geq 1-\alpha$ for $\alpha\in (0,1)$ Lei2013. When specialized to a treatment effect setting, e.g. when one considers $X = \tau_{1,T}$, the set $\mathcal{I}_X$ will be a function of the data that, over repeated sampling of the outcomes, contains the (stochastic) individual treatment effect $\tau_{1,T}$ with probability at least $1-\alpha$. Prediction intervals have a long tradition in the forecasting literature Brockwell1991, and have been recently considered in causal panel data settings with few treated units by cattaneo2021prediction,cattaneo2023uncertainty and chernozhukov2021exact.

Notice that, if both $M_{1,T}$ and the distribution of $\epsilon_{1,T}$ were known, one could construct a valid prediction interval for $Y_{1,T}(0)$ by taking $\mathcal{I}_{Y(0)} = [M_{1,T}+Q_{\epsilon_{1,T}}(\alpha/2), M_{1,T}+ Q_{\epsilon_{1,T}}(1-\alpha/2)]$, where $Q_{\epsilon_{1,T}}$ is the quantile function of $\epsilon_{1,T}$. Given this interval, and by noticing that $\tau_{1,T} = Y_{1,T}(1) - Y_{1,T}(0) = Y_{1,T} - Y_{1,T}(0)$, one could construct an interval for $\tau_{1,T}$ as $\mathcal{I}_{Y(1)} = Y_{1,T} - \mathcal{I}_{Y(0)}$. In practice, $M_{1,T}$ and $Q_{\epsilon_{1,T}}$ are unknown and have to be estimated. Moreover, methods that replace the unknown quantities in $\mathcal{I}_{Y(0)}$ by consistent estimators, whilst asymptotically justified, may exhibit relevant finite-sample distortions. cattaneo2021prediction propose a method for constructing prediction intervals that explicitly takes into account finite-sample estimation error. They consider a general class of synthetic control estimators $\hat{M}_{1,t}=X_{t}'\hat\beta$, where $X_t$ is a vector of features. They note the estimation error of the treatment effect estimator $\hat\tau_{1,T}=Y_{1,T}-X_{T}'\hat\beta$ can be decomposed into two parts, $ \hat\tau_{1,T}-\tau_{1,T}=X_{T}'(\beta-\hat\beta)+\epsilon_{1,T}, $ where $\beta$ is the population (constrained) linear prediction coefficient. Motivated by this decomposition, they propose procedures for constructing separate prediction intervals for each component. Denote these prediction intervals by $[L_{X_{T}'(\beta-\hat\beta)}, U_{X_{T}'(\beta-\hat\beta)}]$ and $[L_{\epsilon_{1,T}}, U_{\epsilon_{1,T}}]$. Based on these separate prediction intervals, the approximately valid prediction interval for the treatment effect $\tau_{1,T}$ is $\hat {\mathcal{I}}_{Y(1)}=[\hat\tau_{1,T}-U_{X_{T}'(\beta-\hat\beta)}-U_{\epsilon_{1,T}},\hat\tau_{1,T}-L_{X_{T}'(\beta-\hat\beta)}-L_{\epsilon_{1,T}}]$. By the union bound, the probability of miscoverage of $\hat {\mathcal{I}}_{Y(1)}$ is bounded above by the sum of the probability of miscoverage of each interval. cattaneo2023uncertainty extend these ideas to construct prediction intervals for a broader class of effects, including various types of average effects, in settings with staggered treatment adoption.

\paragraph{Relationships between inference on sharp nulls, inference on realized effects, and prediction intervals.} Methods for constructing tests in settings with nonstochastic treatment effects can sometimes be inverted to produce valid prediction intervals in a setting where treatment effects are random. For example, chernozhukov2021exact show that their conformal inference method yields valid prediction intervals when the treatment effects are random, while chernozhukov2024ttest, discussed in more detail in Section (ref) below, provide a similar result in the context of making inferences on average effects. More generally, AF2025 provide conditions under which any inference method that is valid for inference on the realized treatment effects assuming the distribution of treatment effects is independent of the distribution of errors, is also valid to construct prediction intervals even when we allow for arbitrary dependence between treatment effects and errors. This includes, for example, the methods proposed by conley20211inference and ferman2019inference. AF2025 also discuss how prediction intervals and tests of sharp nulls are intrinsically connected --- each being obtained from “inverting” the other under the same set of assumptions --- in a wide class of inference methods considered in settings with few treated units. Taken together, these results provide a unifying view on the interpretation of several approaches to inference available in the literature.

\paragraph{Sensitivity analysis.} Finally, if researchers are interested in $\tau^\ast$ and not willing to change the target parameter, they can perform an analysis to assess the sensitivity of the results with respect to the distribution of the heterogeneous treatment effects, or the distribution of $Y_{1,T}(1)$. hagemann2020inference introduces a valid inference method for a setting with a single treated cluster in an asymptotic framework where the number of observations within each cluster is large but the number of clusters is fixed. In hagemann2020inference's setting, if a normal approximation holds for the average of outcomes in each cluster, then a valid test of a null $\tau^\ast = c$ can be conducted, provided that the researcher provides an upper bound for $ \bar{\rho} \equiv \mathbb{V}[\bar{Y}_{1,T}(1)]/\min_{j \in \mathcal{J}_0}\mathbb{V}[\bar{Y}_{j,T}(0)]$, where $\bar{Y}_{j,T}(d)$ denotes the average potential outcome in cluster $j$, with the first cluster being the treated one.\footnote{hagemann2020inference's results also cover the more general case where one of the control clusters has zero variance. In this case, a bound on the ratio between the variance of the average outcome in the treated cluster and the second smallest variance among control clusters must be specified.} Note that hagemann2020inference's method can be alternatively recast as a sensitivity analysis: given a significance level $\alpha$, one can find the smallest-value of $\bar{\rho}$ compatible with not rejecting the null.

Methods exploiting $N_1>1$

The methods reviewed in Section (ref) typically also work in settings with $N_1$ fixed, but greater than one. However, settings with $N_1 > 1$ provide other alternatives for inference. Since we have more variation in the data on potential outcomes under treatment, we do not necessarily have to rely on information from the control units or from the pre-treatment periods to learn about the errors on the treated. This allows us to accommodate heteroskedasticity and stochastic treatment effect heterogeneity. Therefore, we consider the case in which $\tau_{j,t}$ may be stochastic, and we focus on the target parameter $\tau^\ast$ defined in Equation (ref).

While these methods can accommodate heteroskedasticity and stochastic treatment effect heterogeneity, they generally rely on assumptions that would not be required by the methods reviewed in Section (ref). Moreover, these methods usually still rely on stronger assumptions than standard inference methods that are asymptotically valid in settings with many treated and many control units. Finally, an important practical consideration is that these methods may have low or even trivial power when $N_1$ is very small. Thus, from a power perspective, the methods presented in Section (ref) might be preferable in some settings with $N_1>1$ when $N_1$ is very small. See Appendix (ref) for simulation evidence.

Sign-changes & Wild Bootstrap

\paragraph{Idea and implementation of sign-changes.} We start by considering randomization inference methods based on sign changes. To understand the main idea of these methods, consider a cross-sectional setting where $N_1$ is fixed and the estimator $\hat \tau$ can be written as $\hat \tau = \frac{1}{N_1} \sum_{j \in \mathcal{J}_1} \hat \tau_j$, with $\hat \tau_j \overset{p}{\rightarrow} \tau^\ast + (\tau_j - \tau^\ast) + \epsilon_j$ when $N_0 \rightarrow \infty$. This setting encompasses a comparison of means with $N_1$ fixed and $N_0 \rightarrow \infty$. In this case, if we define $\hat \tau_j = Y_j - \frac{1}{N_0} \sum_{j' \in \mathcal{J}_0} Y_{j'}$ for each $j \in \mathcal{J}_1$, then $\hat \tau_j \overset{p}{\rightarrow} \tau_j + \epsilon_j$ when $N_0 \rightarrow \infty$, if we can apply a law of large number for the average errors of the controls. Recall also that a DiD setting with no variation in treatment timing can be recast as a comparison of means, so this would be another example.

The main idea in this case is that, if we assume that $(\tau_j - \tau^\ast) + \epsilon_j$ are independent across $j$ and symmetric about zero, then the asymptotic distribution of $(\hat \tau - \tau^\ast)$ would be invariant to the group of transformations $\mathcal{G} \equiv \{-1,1\}^{N_1}$, meaning that, for any $g = (g_1,g_2,\ldots, g_{N_1})' \in \mathcal{G}$, the asymptotic distribution of $\hat \tau^g(\tau^*) =\frac{1}{N_1} \sum_{j \in \mathcal{J}_1} g_j (\hat \tau_j - \tau^\ast)$ would not depend on the choice of $g$. Therefore, we can apply the theory of randomization inference under approximate symmetry from Canay2017 (see also CaiCanay2023 for implementation details). To test the null $H_0: \tau^\ast = c$, we compute $\{\hat \tau^g(c)\}_{g \in \mathcal {G}}$, and compare $(\hat \tau - c)$ to this randomization distribution. The p-value would be given by $\hat p = \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}}{\mathbf{1}\{|\hat \tau -c | \leq |\hat \tau^g(c) | \}}$. Under the null, this p-value satisfies $\limsup_{N_0\to \infty} P[\hat p \leq \alpha] \leq \alpha$, for any $\alpha \in (0,1)$. Therefore, asymptotically, this test controls for size.

Intuitively, under the null hypothesis that $\tau^\ast = c$, the asymptotic distribution of $(\hat \tau_1 - c, \ldots, \hat \tau_{N_1} - c)$ is centered around zero. If we also assume symmetry, its asymptotic distribution is invariant to sign changes. That is, $(\hat \tau_1 - c, \ldots, \hat \tau_{N_1} - c)$ has the same asymptotic distribution as $(g_1(\hat \tau_1 - c), \ldots, g_{N_1}(\hat \tau_{N_1} - c))$ for any $g \in \mathcal{G}$. As a result, the probability that $\hat \tau$ is among the $k$ largest values of $\{\hat \tau^g(c)\}_{g \in \mathcal{G}}$, would (asymptotically) be $k/\operatorname{dim}(\mathcal{G})$, which guarantees that the test controls for size. Under the alternative hypothesis $\tau^\ast > c$, however, each $\hat \tau_j - c$ tends to be positive. In this case, flipping signs at random typically reduces the magnitude of the average, so most $\hat \tau^g(c)$ will be smaller than the observed $\hat \tau - c$. This shift in distribution leads to small p-values, giving the test power against alternatives where $\tau^\ast \neq c$.

Note that assuming that $(\tau_j - \tau^\ast) + \epsilon_j$ is symmetric allows for heteroskedasticity in the error term $\epsilon_j$, and for stochastic treatment effects heterogeneity $(\tau_j - \mathbb{E}[\tau_j])$, which was not generally allowed in the inference methods considered in Section (ref). However, these methods require assumptions that were not required by the inference methods considered in Section (ref), such as symmetry. This symmetry condition is often justified in settings in which outcomes $Y_j$ represent averages of individual level or time-series observations within group $j$, under assumptions that allow us to rely on a CLT within each unit Canay2017. Moreover, these methods do not allow for deterministic treatment effects heterogeneity, as in this case $\tau_j - \tau^\ast$ would not be symmetric about zero for all $j\in \mathcal{J}_1$. This is a more general feature of methods discussed in this section: they may allow for some stochastic heterogeneity, but, since there is only a finite number of treated units, there is not enough information to discern between nonstochastic location shifts and stochastic treatment effect heterogeneity. In contrast, methods in previous sections allowed for deterministic heterogeneity by extrapolating from controls or pre-treatment periods.

\paragraph{Power issues.} A practical limitation of inference methods based on sign-changes is that they never reject the null if $\alpha<1/\operatorname{dim}(\mathcal{G})$, and thus have trivial power in this case. In the extreme case with $N_1=1$, there are only two possible transformations with $|\hat \tau^g| = |\hat \tau|$ for $g \in \{-1,1\}$, so the p-value is equal to one and the test has trivial power for any $\alpha$. Likewise, with $N_1=2$ we would only have two distinct values for $|\hat \tau^g|$, so $\hat p$ could only be equal to 0.5 or 1. For a test at the 10% (5%) level to have non-trivial power, i.e. for the rejection probability to be greater than zero, there must be at least five (six) treated units CaiCanay2023. Therefore, the alternatives reviewed in Section (ref) might be preferred in settings with $N_1$ very small (even when $N_1>1$), due to power considerations. See Appendix (ref) for an illustration.

\paragraph{Applications of sign-changes.} Inference methods based on sign-changes method have been considered in DiD settings with few treated and many control units when we have uniform treatment timing Canay2017. In this case, we define $\epsilon_j$ as the post-pre average errors. These methods have also been considered for matching estimators with few treated and many control units ferman_matching. In this case, for $j \in \mathcal{J}_1$, we can define $\hat \tau_j$ as the difference between $Y_j$ and the average of its nearest neighbors. ferman_matching provides conditions under which $\hat \tau_j \overset{p}{\rightarrow} \tau^\ast + (\tau_j - \tau^\ast) + \epsilon_j - \xi_j$ as $N_0 \rightarrow \infty$, where $\xi_j$ is the average of the errors of the nearest neighbors of treated unit $j$. If units are independent and the probability that two treated units share the same nearest neighbor goes to zero, then dependence between $\hat \tau_j$'s would go to zero when $N_0 \rightarrow \infty$. Therefore, the sign-changes test would be asymptotically valid when $N_1$ is fixed and $N_0 \rightarrow \infty$ if $(\tau_j - \tau^\ast) + \epsilon_j - \xi_j$ is symmetric about zero. ferman_matching considers finite-sample corrections to take into account that, in finite samples, treated units may share the same nearest neighbor. ferman_matching also conjectures conditions under which the sign-changes test may be used in synthetic control applications with more than one treated unit. Randomization tests based on sign changes are also applicable, for example, in RCTs with few treated and many control units (so we consider an asymptotic approximation with $N_1$ fixed and $N_0 \rightarrow \infty$). Therefore, relative to the methods considered in Section (ref), sign-changes randomization tests provide an alternative (in a model-based framework) that allows for stochastic treatment effect heterogeneity at the expense of assuming symmetry.

s\paragraph{Equivalence to wild bootstrap with the null imposed.} In Appendix (ref), we show that the sign-changes test is asymptotically equivalent to a wild bootstrap with the null imposed when $N_1$ is fixed, $N_0 \rightarrow \infty$, $\hat \tau_j \overset{p}{\rightarrow} \tau^\ast + (\tau_j - \tau^\ast) + \epsilon_j$, and $\hat{\tau} = \frac{1}{N_1}\sum_{j\in \mathcal{J}_1} \hat{\tau}_j$ can be obtained from a (partially) linear regression. This result provides a theoretical justification for the use of the wild bootstrap with the null imposed in settings with few (but more than one) treated and many control clusters, which has been a common recommendation for settings with few treated clusters, at least since CGM.\footnote{Canay_wild_bootstrap provide a theoretical justification for the Wild Cluster Bootstrap with null imposed for inference in linear models in settings with a fixed number of clusters, when the number of observations in each cluster is large. While their main theory is not directly valid in settings in which we consider clusters at the treatment assignment level, they show validity of a wild cluster bootstrap when we consider clustering at a coarser level. In Appendix (ref), we show that this alternative is asymptotically equivalent to a wild bootstrap at the treatment assignment level when $N_0 \rightarrow \infty$ and $N_1$ is fixed.} We emphasize that imposing the null is crucial here: the wild bootstrap without imposing the null generally does not control size in settings with few treated units or clusters MacKinnon2018.

\paragraph{Sign-changes with finite $N_1$ and $N_0$ validity.} In some cases, it is also possible to consider a sign-changes method that is valid when both $N_1$ and $N_0$ are fixed. For example, Canay2017 discuss how their proposed sign-changes test may be adapted to settings with fixed $N_1$ and $N_0$, under additional assumptions. In this case, the idea is to cluster the units in $N_1$ disjoint groups, each one containing one treated and some of the control units. Suppose that the adopted estimator is “linear”, in the sense that $\hat \tau = \frac{1}{N_1} \sum_{j \in \mathcal{J}_1} \tilde \tau_j$, with $\tilde{\tau}_j$ being an estimator of the treatment effect of unit $j$, using data from cluster $j$. Suppose each $\tilde{\tau}_j$ may be decomposed as $\tilde \tau_j = \tau^\ast + (\tau_j - \tau^\ast) + \tilde \epsilon_j$, where $\tilde \epsilon_j$ is a function of the errors in the $j$-th cluster. If the $(\tau_j - \tau^\ast) + \tilde \epsilon_j$ are independent and symmetric about zero across $j$, the sign-changes test of the null $\tau^\ast = c$ based on $ \frac{1}{N_1}\sum_{j \in \mathcal{J}_1} g_j (\tilde \tau_j - c)$ and on the group of transformations $(g_1,\ldots,g_n) \in \mathcal{G}$, is valid with fixed $N_0$.

One detail in this case is that the p-value would depend on the division of the units into disjoint groups. As a solution, CaiCanay2023 recommend users “combine” the conclusions from different partitions by relying on the strategies proposed by DiCiccio2020 in the context of a general hypothesis testing problem. Alternatively, we show in Section (ref) of Appendix (ref) that another valid solution to this problem is to consider a number of different partitions, and then to construct an “aggregate” p-value across these partitions, following an idea similar to song2018orderingfree and Leung. By leveraging the special structure of the sign changes test and standard results on randomization tests, this alternative yields an inference procedure that is easier to implement and generallly more powerful than the solutions in DiCiccio2020. Moreover, as we remark in Appendix (ref), this aggregate p-value is asymptotically equivalent to the standard sign-changes test when $N_0 \rightarrow \infty$. Therefore, compared to the finite $N_0$ version of the sign-changes method proposed by Canay2017, this alternative has the advantage of not depending on the partition of the units, without requiring additional assumptions. At the same time, it is asymptotically equivalent to the standard implementation of the sign-changes method proposed by Canay2017 when $N_0$ is large.\footnote{ lau2025 provides an alternative approach that seeks to find a single partitioning scheme in order to maximize (local) power. This alternative is valid with $N_0$ fixed and (approximately) Gaussian $\tilde \epsilon_j$.}

Behrens-Fisher solutions

If we impose normality on the potential outcomes, then a treatment/control cross-sectional comparison or a DiD setting with uniform treatment timing collapses to a version of the classical Behrens-Fisher problem where the variance can be different between and within groups Bakirov1998. In theses cases, there are several alternatives available in the literature that could be used for settings with fixed $N_1$ and $N_0$. The assumption of normality may be justified as holding at least approximately in settings where each unit consists of an average of many observations, and it is plausible to consider a CLT for these within-unit averages. However, this justification for relying on normality approximations precludes the possibility of within-unit aggregate shocks.

Ibragimov2016 provide conditions on sample sizes and significance levels under which a degrees-of-freedom adjustment to an unequal variance $t$-test for the comparison of means of two populations is conservative when observations are independently normally distributed with possibly heterogeneous variances.\footnote{Bloom apply Ibragimov2016's (Ibragimov2016) approach to an RCT with $N_1 = 11$ and $N_0 = 7$. They estimate effects for unit outcomes that consist of a before-after comparisons with a large number of time periods, so the Gaussian approximation for each individual observation is appropriate.} {More recently, potscher2023 extended the results of Ibragimov2016 to a general Gaussian linear model. They provide sufficient conditions for the existence of modified critical values that ensure the heroskedasticity-robust t-test controls size in finite sample.\footnote{Their results also extend to some non-Gaussian settings, provided that the standardized errors of the model are spherically and symmetrically distributed with no point mass at zero.} Their critical values coincide with Ibragimov2016 in the comparison-of-means problem for the ranges of sample sizes and significance values in Ibragimov2016, but are also computable for other sample sizes and significance values. potscher2024 shows that the sufficient condition in potscher2023 is in fact necessary for finite-sample size control of tests based on a t-statistic with heteroskedasticity-consistent standard errors}; and preinerstorfer2021 provides computational code to find the modified critical values.

In a setting similar to the one from Ibragimov2016, Hagemann2023_permutation constructs a (generally conservative) permutation test for the null that the difference in means between the two populations is equal to some value $c$. His results cover different configurations of test statistics, sample sizes and significance levels than Ibragimov2016.

Still considering alternative inference procedures, Ibragimov2010 consider settings in which the parameter of interest can be estimated separately in a finite number of independent clusters. In DiD designs, this may require considering coarser clusters, containing both treated and control units, similarly to the discussion in Section (ref). In such settings, they provide conditions on sample size and significance levels for their procedure to be conservative. Under an additional assumption that estimators computed in each cluster have the same variance, we show in Appendix (ref) that their procedure in DiD designs collapses to a modified version of the method of Bester2011, which is exact regardless of sample size and significance level. Finally, hansen2024jackknife introduces a jackknife variance estimator for linear regressions that is never downward-biased, while Hansen2025 studies its application to clustered DiD designs. He shows through simulations that confidence intervals that rely on this variance estimator can have better coverage when compared to confidence intervals based on cluster-robust standard errors.

Overall, a key point is that a common feature of the normality-based methods discussed in this section is that they become invalid or have trivial power when $\min \{N_1,N_0\} = 1$, meaning they are only viable in settings with $N_1,N_0>1$. Moreover, even when $\min \{N_1,N_0\} > 1$, these methods can have lower power than those discussed in Section (ref) if the number of treated units is very small. However, they can exhibit non-trivial power with small $N_1$ in settings where the alternatives from Section (ref) would have trivial power. See Appendix (ref) for an illustration. We recall, though, that power comparisons for methods that rely on non-nested assumptions should be considered with caution.

Dual justification: stronger assumptions with $N_1$ fixed & weaker assumptions with $N_1,N_0 \rightarrow \infty$

A common theme for methods that are valid with fixed $N_1$ (whether with fixed $N_0$ or $N_0 \rightarrow \infty$) is that they rely on stronger assumptions relatively to standard methods that are valid with large $N_1$ and $N_0$. This motivates inference methods that are valid under stronger assumptions when $N_1$ (or $N_1$ and $N_0$) is fixed and under weaker assumptions when $N_1,N_0 \rightarrow \infty$ (in some cases, there are restrictions on the behavior of $N_1/N_0$, such as $N_1/N_0 \rightarrow 0$). One way to achieve that in methods that are based on randomization inference (such as those based on permutations or sign changes) is to consider studentized test statistics, as considered by, for example, Janssen1997, Chapter 15 of Lehmann2005, Chung2013, DiCiccio2017, MACKINNON2020435, ferman_matching, Canay_wild_bootstrap, Bertanha2023 and DHaultfuille2024. This is also a feature of the approaches considered by Ibragimov2016 and Bester2011, as it is well-known that $t$-tests with heteroskedasticity-robust (cluster-robust) standard errors are asymptotically valid under much weaker conditions when there are many treated and untreated units (clusters). Following a similar idea, Chaisemartin2022 propose a testing procedure in a DiD setting that is exact in finite samples under normality, homoskedasticity, and treatment effect homogeneity; and that remains asymptotically valid when the data are non-normal, heteroskedastic, and treatment effects are heterogeneous. Their construction relies on a $t$-test with a heteroskedasticity-robust standard error and critical values simulated under the normality, homoskedasticity, and treatment effect homogeneity assumptions. In large samples, these simulated critical values collapse to those obtained from the standard normal distribution, ensuring the asymptotic validity under weaker assumptions.

Higher-order improvements to standard asymptotics

Another related set of approaches consists in constructing inference procedures that are valid under standard asymptotics as $N_0,N_1\rightarrow \infty$, while achieving better performance in finite samples. This includes corrections to standard error formulae with an aim to remove higher-order bias terms and downweigh high-leverage data points Mckinnon1985,Cribari12000,Cribari2004,Cribari2007,mackinnon2023cluster, and approaches that, building on an idea originally due to Welch1951, seek to compute improved critical values for $t$-tests by relying on a t-distribution with degrees of freedom chosen in a way such that moments of the ratio between the variance estimator and its population counterpart mimic those of a chi-squared distribution bell2002bias,Imbens2016,young2016improved,Hansen2025. Studentized bootstrap methods that achieve higher-order improvements can also be placed in this category CGM,cameron2015practitioner,Djogbenou2019,mackinnon2023cluster. It is important to note that the theoretical justifications of these approaches still rely on asymptotic approximations where $N_1 \to \infty$, which may be inaccurate if $N_1$ is small.

Methods that exploit within-cluster Central Limit Theorems

In this section, we consider methods that exploit CLT approximations for averages of within-cluster observations. In Section (ref), we consider a setting where $N_1=1$ and there are enough post-treatment periods to justify asymptotic regimes where $T_1\rightarrow \infty$. In this case, it is possible to obtain asymptotically normal ATT estimators under stationarity and weak dependence assumptions on the errors, which motivate normality-based inference methods. In Section (ref), we consider a setting in which the number of observations within each treated cluster goes to infinity, but we do not have information on the dependence structure between observations (for example, when the number of individual-level observations within each treated cluster goes to infinity). In this case, under weak dependence assumptions on the errors, we may still have that the ATT estimator is asymptotically normal, and it is possible to conduct inference based on subsampling methods, while remaining agnostic to the dependence structure within clusters.

Central limit theorems exploiting time series ($T_1 \rightarrow \infty$)

Consider a setting where $N_1=1$, but $T_0$ and $T_1$ are large. The target parameter is the ATT for the treated unit over the entire post-treatment period, $\tau^\ast = \mathbb{E}\left[ T_1^{-1} \sum_{t \in \mathcal{T}_1}\tau_{1,t} \right]$, where $\tau^\ast$ may change with $T_1$, but we suppress such dependence for notational convenience. In this setting, the availability of many post-treatment periods makes it possible to develop inference methods based on CLTs.\footnote{For example, methods exploiting large $T_1$ in panel data settings have been proposed by, for example, li2017estimation,carvalho2018arco,li2020statistical,synthetic_did,masini2021counterfactual_jasa,masini2022counterfactual_jbes,li2023statistical,chernozhukov2024ttest. Some of those methods would nest classical DiD as a special case.} Ideally, one would like to establish asymptotic normality results, such as $\sqrt{T}_1(\hat\tau -\tau^\ast)\overset{d}\rightarrow N(0,\sigma_\tau^2)$ as $T_1\rightarrow \infty$, where the asymptotic variance $\sigma_\tau^2$ can be consistently estimated. However, there are at least four major challenges.

The first challenge is dealing with the error from estimating $M_{1,t}$. Since this term is rescaled by $\sqrt{T}_1$, it is not sufficient to just have a consistent estimator for $M_{1,t}$. Moreover, as discussed in Section (ref), some approaches for constructing $\hat{M}_{1,t}$ involve high-dimensional estimation problems (e.g., synthetic control with many control units). In such cases, the simple plug-in estimator $\hat \tau = \frac{1}{T_1} \sum_{t \in \mathcal{T}_1} (Y_{1,t} - \hat{M}_{1,t})$ will typically be biased, which motivates bias correction procedures chernozhukov2024ttest.

The second challenge is to estimate the asymptotic variance in the presence of serial correlation. The estimation error $\sqrt{T}_1(\hat\tau -\tau^\ast)$ typically contains terms like $\frac{1}{\sqrt{T_1}} \sum_{t \in \mathcal{J}_1} \epsilon_{1,t}$, which satisfy a CLT if $\{\epsilon_{1,t}\}$ is stationary and weakly dependent, in which case $\frac{1}{\sqrt{T_1}} \sum_{t \in \mathcal{J}_1} \epsilon_{1,t}\overset{d}\rightarrow N(0,\sigma_\epsilon^2)$. A standard approach to estimate $\sigma_\epsilon^2$ when $\{\epsilon_{1,t}\}$ exhibits serial correlation is to use newey1987simple variance estimators li2017estimation,carvalho2018arco. However, these estimators can perform poorly when the number of periods is small or moderate. To overcome this challenge, chernozhukov2024ttest propose an inference method that avoids the estimation of $\sigma_\epsilon^2$. The idea is to construct a test statistic that is “self-normalized”---a test statistic in which the denominator and the numerator are proportional to $\sigma_\epsilon$, which thus cancels out.

The third challenge is dealing with non-stationary data. As mentioned in Section (ref), non-stationarity is ubiquitous in applications where exploiting the time series dimension is useful for inference. One approach for dealing with non-stationarity is to impose explicit time series models for the non-stationarity and analyze the properties of specific estimators under these models li2020statistical,masini2021counterfactual_jasa,masini2022counterfactual_jbes. The drawback of this approach is that the resulting inference methods are typically not robust against violations of the underlying models of non-stationarity. Instead of imposing specific models for the non-stationarity, one can alternatively impose restrictions on the heterogeneity in the non-stationarity across units chernozhukov2024ttest. We view these two approaches as complementary.

The fourth challenge is that the estimator may not be asymptotically normal. For example, li2020statistical shows that synthetic control estimators may not be asymptotically normal in this case, given the restrictions on the synthetic control weights. In this case, she suggests a subsampling procedure to overcome the non-standard asymptotic distribution of the estimator.

In addition to these challenges, there is a trade-off between allowing for stochastic and deterministic treatment effect heterogeneity. If the treatment effect sequence $\{\tau_{1,t}\}_{t\in \mathcal{T}_1}$ is deterministic, we can typically allow for arbitrary effect heterogeneity over time carvalho2018arco,chernozhukov2024ttest.\footnote{As discussed in Section (ref), the confidence intervals obtained using the methods that treat $\{\tau_{1,t}\}_{t\in \mathcal{T}_1}$ as deterministic can often be reinterpreted as prediction intervals if $\{\tau_{1,t}\}$ is stochastic. See AF2025 for further discussion.} If the treatment effects are stochastic, stationarity and weak dependence assumptions on $\{\tau_{1,t}\}_{t\in \mathcal{T}_1}$ are typically required li2017estimation,li2020statistical,chernozhukov2024ttest, ruling out deterministic treatment effect heterogeneity.

On the one hand, settings with many post-treatment periods allow for developing inference methods based on asymptotic normality results. These methods typically are standard and easy to implement and communicate. Moreover, some of these methods can relax assumptions required by methods discussed in Sections (ref) and (ref), such as homoskedasticity and the absence of stochastic treatment effect heterogeneity. On the other hand, the validity of these methods relies on asymptotic approximations, which may be inaccurate in applications where $T_1$ is small or moderate.

Accounting for within-cluster dependence without a distance metric

The methods in the previous section relied on estimating the dependence structure between observations. This was possible because a metric along which the dependence between observations was assumed to decay -- in that case, time -- was naturally available. In other settings, however, treated observations may be arranged in clusters where no natural metric is available. For example, when we have many individual-level observations within each cluster. For these cases, Leung proposes a general inference method that remains agnostic about the dependence structure between observations. His approach remains valid even in settings with a single treated and a single control cluster, provided that the number of observations in each cluster ($N_j$) is large, and the conditions for the validity of a CLT on the within-cluster averages hold. Relatively to the methods from Section (ref) and (ref), this approach requires access to individual-level data, and it requires the use of a CLT within cluster (so it precludes, for example, cluster-level shocks, which are allowed by these other methods).\footnote{The discussion in Remark (ref) on inference conditional on aggregate shocks applies to this setting.}

In a comparison of means setting with one treated and one control cluster, Leung's approach to inference relies on a standard two-sample $t$-test statistic, with the average in each cluster being computed by randomly resampling $R_j$ units with replacement from each cluster $j$. Provided that $R_j \to \infty$ with $R_j/N_j \to 0$, the results in Leung ensure the test statistic converges in distribution to a standard normal, which enables researchers to construct tests with asymptotic validity. Importantly, the procedure remains agnostic about the dependence structure within each cluster (except for assuming that the dependence is such that a CLT within clusters is valid, which requires some sort of weak dependence within clusters). Intuitively, this is due to the restriction that $R_j/N_j \to 0$, which ensures that resampled draws are approximately independently distributed, thus ensuring convergence to a standard normal even in the presence of dependence. Indeed, as argued by Leung, the restriction that $R_j/N_j \to 0$ may be seen as a price to pay in order to be agnostic about the dependence structure of observations, since we effectively “lose” observations by working with $R_j < N_j$ draws. Presently, there is no general method to choose $R_j$, though Leung provides some guidance in specific settings.\footnote{The choice of $R_j$ is subject to a trade-off between statistical power and size control: a larger $R_j$ increases the power of the test, though possibly at the cost of a poorer quality of the normal approximation.}

Design-based inference

Notation and sources of uncertainty

An alternative approach is to consider design-based inference, in which (at least part of) the stochastic variation comes from the treatment assignment. This approach has a long tradition in the analysis of experiments (Neyman1990, Neyman1990 [1923]; Fisher1992, Fisher1992 [1926]). In design-based analyses, potential outcomes of the sample are often considered as fixed, $\{Y_j(1),Y_j(0) \}_{j=1}^N$, and uncertainty would come from the treatment assignment $\{D_j\}_{j=1}^N$.\footnote{To be consistent with the notation from Section (ref), we continue to consider $\{Y_j,D_j\}_{j=1}^N$ as the sample, which differs from Abadie_finitepop, who let $N$ be the number of observations in the finite population.} The target parameter can be the sample average treatment effect (SATE), $\tau_{\mbox{\tiny SATE}} \equiv \frac{1}{N} \sum_{j=1}^N (Y_j(1) - Y_j(0))$. If treatment is randomly assigned, then the standard difference in means estimator, $\hat \tau = \frac{1}{\sum_{j=1}^N D_j}\sum_{j=1}^N Y_j D_j - \frac{1}{N - \sum_{j=1}^N D_j}\sum_{j=1}^N Y_j (1-D_j)$ would be unbiased for the SATE, in this framework in which the only source of uncertainty comes from the assignment of treatment Imbens_Rubin_2015. Alternatively, we may consider the sample $\{Y_j,D_j\}_{j=1}^N$ as being drawn from a larger (finite or infinite) population, and define the target parameters as the Population Average Treatment Effect (PATE). In this case, uncertainty quantification should account for both randomness in the assignment mechanism, as well as sampling.

Note that the focus in such settings is generally on parameters related to average treatment effects given the realized potential outcomes of the sample (or of a finite population), while in Section (ref) the focus was on parameters related to the average treatment effects on the treated over different realizations of the potential outcomes.\footnote{An exception is Rambachan_designbased, who also define the expected average treatment effect on the treated, as an analog of the ATT in model-based settings. Still, this target parameter is also defined conditional on the realization of the potential outcomes, which differ from the target parameter in model-based approaches.} Therefore, from a conceptual perspective, a decision on whether to focus on model-based or design-based approaches for inference should reflect the target parameter of interest. From a pragmatic perspective, design-based approaches might be preferable when it is hard to conceive a super-population, while there is knowledge about the treatment assignment mechanism roth2023what. An extreme example is when the treatment is randomly assigned, in which case knowledge of the assignment mechanism can be exploited to construct exact inference procedures. Design-based approaches might also be helpful in settings in which it is possible to conceptualize a super-population, but modeling the outcome process is difficult. This is especially useful in settings where it is difficult to posit the dependence structure underlying sampling uncertainty, as design-based procedures remain agnostic about these Barrios2012,Adao2019. More recently, design-based approaches have also been considered in natural experiments, including in settings where treatment assignment probabilities are not equal across units Rambachan_designbased.

Inference methods in design-based settings

The use of heteroskedasticity-robust standard errors, or clustered standard errors, has been considered in this design-based framework for testing hypotheses regarding the SATE or PATE, with a justification that they conservatively quantify uncertainty. Indeed, Abadie_finitepop show that, in a cross-sectional treatment-control comparison, the two-sample unequal variance estimator of Neyman1990 (Neyman1990 [1923]) conservatively reflects estimation uncertainty of both the SATE and the PATE.\footnote{This variance estimator corresponds to relying on a heteroskedasticity-robust standard error with the degrees of freedom correction of Mckinnon1985 Imbens2016,Abadie_finitepop.} They then extend these results to a more general linear regression setting, where they show that heteroskedasticity-robust standard errors are generally asymptotically conservative in reflecting estimation uncertainty on weighted versions of the SATE and PATE. Abadie2022 provide similar results for cluster-robust standard errors, whereas Athey2022 and Roth2023 provide conservative variance estimators for estimators of weighted versions of the SATE in staggered DiD designs with randomized treatment adoption dates. See Adao2019 and Rambachan_designbased for related results in other quasi-experimental settings.

The previous results provide a theoretical justification for reporting conventional standard error formulae for estimators of the SATE or PATE. However, it is important to note that, for inference using $t$-statistics and normal critical values to be appropriate, these justifications rely on CLTs that presume that $N_1,N_0\rightarrow \infty$. In particular, if we are in a setting with few treated units, then methods based on asymptotic approximations can perform poorly.

Randomization tests

In settings where the treatment assignment mechanism is known (e.g., RCTs), randomization tests constitute a natural alternative for conducting design-based inference fisher1949design,Imbens_Rubin_2015,young_QJE. Randomization tests are exact in finite samples for testing sharp null hypotheses, such as $H_0: Y_j(1)=Y_j(0)$ for all $j$, and are therefore well-suited for applications with a small number of observations (or with few treated observations). They are valid even when there is only a single treated unit. To illustrate the main idea of randomization tests, consider the simple difference-in-means estimator.

In an RCT, the distribution of the treatment assignment vector $(D_1,\dots,D_N)$ is known, so that if we knew both potential outcomes for all $j=1,...,N$, we would know the distribution of $\hat \tau$. Sharp null hypotheses such as $H_0: Y_j(1)=Y_j(0)$ for all $j$ allow us to compute both potential outcomes for each unit as $Y_j=Y_j(1)=Y_j(0)$, so that the distribution of $\hat \tau$ is known under the null. For example, in a completely randomized experiment with $N_1$ treated and $N_0$ control units, this distribution can be computed by recalculating $\hat \tau$ under each of the ${N\choose N_1}$ possible treatment assignments with $N_1$ treated units. Hypotheses tests can then be conducted by comparing the actual estimate to the quantiles of this randomization distribution.

Randomization tests beyond experiments

Randomization test have also been considered in a variety of non-experimental settings. For example, abadie2010synthetic consider the idea of using a permutation test for the synthetic control estimator. They are very clear, though, that this should not be viewed as a formal test, since we usually cannot guarantee that treatment assignment probabilities are uniform across units. FirpoPossebom+2018 formalize this procedure for settings with uniform assignment, and consider sensitivity analysis for deviations from the uniform assignment probabilities benchmark. lei2024inference propose an alternative permutation test for synthetic controls under both uniform and nonuniform treatment assignment that improves power over existing alternatives when the pool of controls is small. Shaikh2021 consider randomization tests for the sharp null of no treatment effect in any unit or time period in a panel data setting with staggered treatment adoption, under the assumption that time-to-treatment follows a Cox proportional hazards model. In a setting with few treated units, the validity of their method requires knowledge of the parameters of the assignment mechanism.\footnote{If these parameters are unknown, a valid test can be conducted in a few treated setting by maximizing the p-value that treats the parameters as known over a valid confidence set for these parameters, and then adding a correction for one minus the confidence level of the set Berger1994. } Rosenbaum2002 discusses exact tests of sharp nulls in observational studies under a selection on observables assumption. His construction requires the propensity score to follow a logistic model, $\mathbb{P}[D_j=1|X_j] =\operatorname{Logit}(\psi'X_j)$, with $\psi$ possibly unknown.\footnote{Rosenbaum2002 also proposes a randomization test for observational studies with matched pairs, under the assumption that the probability of treatment assignment, within a pair, is the same.} Rosenbaum1996,Rosenbaum2002 and Imbens2005 propose randomization tests for sharp nulls in instrumental variable designs under knowledge of the instrument assignment mechanism (as in an experiment with imperfect compliance). CattaneoFrandsenTitiunik2015 and Cattaneo2017 introduce permutation tests for sharp nulls in regression discontinuity designs, under the assumption that the treatment is locally randomized in a neighborhood of the threshold. Borusyak2023 propose randomization tests for designs that leverage differential exposure to random shocks as a means of identifying causal effects. Their analysis nests shift-share designs as a particular case, and requires researchers to have at least partial knowledge of the shock assignment mechanism. Alvarez_shift_share propose randomization tests in shift-share designs that, by leveraging the special structure of the instrument in these settings, remain robust to misspecification of the assignment mechanism when the number of shocks is large.

Issues with randomization tests with few treated units

In the following, we highlight two important points regarding the use of randomization tests that are particularly relevant when we consider settings with few treated units.

\paragraph{Testing a sharp null.} An important point to notice is that this finite-sample justification of randomization tests does not consider testing null hypotheses regarding the SATE or PATE. Therefore, we may have that $\tau_{\mbox{\tiny SATE}}=0$, but a permutation test would reject at a rate greater than $\alpha$. This could happen, for example, if $\frac{1}{N} \sum_{j=1}^N Y_j(1) = \frac{1}{N} \sum_{j=1}^N Y_j(0)$, but $\frac{1}{N} \sum_{j=1}^N (Y_j(1) - \bar Y(1))^2 > \frac{1}{N} \sum_{j=1}^N (Y_j(0) - \bar Y(0))^2$, where $\bar Y(d) = \frac{1}{N}\sum_{j=1}^N Y_j(d)$ (that is, treatment does not affect the average of the potential outcomes, but affects their variances).

Another way to justify the use of permutation tests is to consider a test for the null $H_0: \tau_{\mbox{\tiny SATE}}=0$, but assume that $Y_j(1) = Y_j(0) + c$ for all $j$, for a constant $c$. Under this assumption, we have that the null $H_0: \tau_{\mbox{\tiny SATE}}=0$ implies $Y_j(1) = Y_j(0)$ for all $j$, which is the main building block to show that the permutation test is valid when we know the distribution of treatment assignment. Assuming that treatment effects are homogeneous in this finite-population setting is similar in spirit to assuming that there is no stochastic treatment effect heterogeneity in a model-based setting (e.g., conley20211inference). This should not be surprising, given the asymptotic equivalence between a permutation test and the method proposed by conley20211inference in model-based settings. This highlights that, whether we are in design- or model-based settings, permutation tests rely on the same kind of restrictions on the treatment effect heterogeneity for exact validity.

Similarly to the discussion in Section (ref), an interesting feature of randomization tests of sharp nulls in design-based settings is that, by properly studentizing the test statistic, it is possible to construct tests that are exact in finite samples for the sharp null, but that are at the same time asymptotically valid (albeit generally conservative) for a weaker null hypothesis. In completely randomized experiments, Wu2020 show that a permutation test of the sharp null $H_0: Y_j(0) = Y_j(1) + c$ that relies on an unequal variance $t$-statistic is exact in finite samples for this sharp null, and asymptotically conservative for the null $H_0: \tau_{\text{SATE}} = 0$ when both the number of treated and control units is large. They provide similar results for stratified experiments. Bugni2018 provide related results for the PATE in stratified experiments, whereas Young2024 considers the case of assessing treatment effect heterogeneity in the population. Alvarez_shift_share provide similar constructions for shift-share designs. Therefore, under a suitable choice of test statistic, the fact that randomization tests are generally designed for testing a sharp null is not a crucial issue when we are in a setting with many treated and many control observations.

\paragraph{Unconditional vs conditional inference.} Randomization tests are valid for unconditional inference. Now consider a setting in which we observe a characteristic $W_j \in \{0,1\}$ for the units in our sample (so that $\{Y_j(1),Y_j(0),W_j\}_{j=1}^N$ is treated as fixed). It might be that after running the experiment, the experimenter observes that there is an important imbalance in terms of the average of $W_j$ between the treated and the control group. Let $\bar W_1$ ($\bar W_0$) be the average of $W_j$ for the treated (control) group. In this case, while the randomization test would remain valid for unconditional tests of the sharp null, it may cease to be valid conditional on the fact that we had an imbalance in $(\bar W_1,\bar W_0)$. In this case, an alternative for valid conditional inference would be to consider only permutations with the same imbalance $(\bar W_1,\bar W_0)$ Hennessy2016. Whether considering unconditional or conditional inference in such settings has been subject to a longstanding debate Mutz2019,Johansson2022. While we do not intend to contribute to this debate, we highlight that this discussion is particularly relevant in settings with few treated units, which is the focus of our survey. In such settings, the probability of having relevant imbalances is higher than with many treated and many control units. In particular, if we consider a setting with $N_1=1$ and $W_j \in \{0,1\}$, then we would have that $\bar W_1 \in \{0,1\}$, so we would have imbalance in all realizations of the treatment assignment.

Conclusion

This survey considers the challenges with inference in applications with few treated units. We present a simple example that illustrates why standard inference methods typically fail in such applications and then discuss alternative inference methods that are valid with one or few treated units. We organize the survey according to the notion of uncertainty (model-based vs.\ design-based uncertainty) and the data availability. Two broad themes emerge. First, methods that are valid with few treated units typically rely on stronger assumptions than standard methods that are valid with many treated units. Second, the choice of method is highly context-specific. It depends, among other things, on the data availability, the target parameter of interest, the notion of uncertainty, and the assumptions one is willing to make.

For applied researchers, this survey first serves as a warning that standard inference methods may perform poorly in empirical applications with few treated units. It then provides guidance on how to choose a suitable alternative inference method, taking into account the specifics of the empirical application at hand.

For econometricians, this survey provides an overview of state-of-the-art inference methods for settings with few treated units, with an emphasis on the connections and differences between the different procedures. We hope that this survey will foster more research in this area. We see the development of novel inference methods relying on alternative sets of assumptions not previously considered in the literature as an important avenue for future research. A wider range of available methods will help empirical researchers choose inference methods that are more specifically tailored to their empirical applications.

{0pt}