Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,542 characters · 15 sections · 38 citation commands
Post-Matching Two-Way Fixed Effects Estimation
Keywords: two-way fixed effects, difference-in-differences, matching \\ \noindentJEL Classification: C14, C21, C23
When estimating treatment effects with two-ways fixed effects (2WFE) regressions, researchers often use matching as a pre-processing step. In such settings, each treated unit is matched to an untreated unit using pre-treatment observed covariates. In a second step, a 2WFE specification is estimated in the matched sample, reweighting the untreated units by the number of times they are used as a match for the treated units. This approach is typically justified under a conditional parallel trends assumption, by which units that are similar in their observed covariates are more likely to have followed the same outcome trends in the absence of treatment. While several alternative methods are available for estimating treatment effects under a conditional parallel trends assumption, matching-based methods have the advantage of being intuitively appealing, easy to implement, computationally stable and fully nonparametric, not requiring the functional form assumptions often invoked in inverse-probability weighting or regression adjustment methods. This combination of matching and 2WFE is common in empirical work: Table (ref) lists 22 papers published between 2013 and 2025 in top general interest and field journals in Economics and Political Science that use some version of this approach.
In this paper, we formally analyze this commonly employed strategy and highlight two problems that generally invalidate its findings. First, when different treatment cohorts enter treatment in different time periods, the post-matching 2WFE estimator that pools all treated cohorts converges to a linear combination of average treatment effects on the treated (ATTs) plus an asymptotic bias, which we refer to as “pooling bias”, even when treatment effects are constant across units and over time. To explain this result, we show that the post-matching 2WFE estimator can be decomposed into a linear combination of 2x2 (that is, two-group-two-period) reweighted difference-in-differences (DiD) estimators. These 2x2 estimators compare each treated cohort to the never-treated cohort, and each treatment cohort to other treatment cohorts, with a specific set of weights. Because each treated unit is matched to a never-treated unit in the first step, comparisons between treated and reweighted never-treated units consistently estimate ATTs for different cohorts at different time periods. When comparing units across different treatment cohorts, however, units are not matched to each other, and thus the group used as comparison is not reweighted. As a result, the two groups will generally exhibit different covariate distributions. When the parallel trends assumption holds conditionally, these unweighted comparisons between different treated cohorts are invalid, even when the comparison group consists of not-yet-treated units.
To address this issue, we propose post-matching difference-in-differences (DiD) estimators that compare each treated cohort to the never-treated units separately, instead of pooling all treated cohorts. We show that under the conditional parallel trends assumption and mild regularity conditions, these pairwise estimators are consistent for well-defined averages of ATTs over time.
The second problem is that matching introduces variability that needs to be accounted for when conducting inference, and failing to account for this additional variability results in invalid standard errors. We characterize the probability limit of the naive cluster-robust variance estimator for the post-matching DiD estimator and show that failing to account for the matching step creates an asymptotic bias in the variance estimator, which can be positive or negative depending on the data generating process. We rely on recent advancements on measures of high-order Voronoi cells by Chen-Han_2024_wp to provide a closed-form representation of the correct limiting variance of our proposed estimators, a result previously unavailable in the literature, and propose consistent variance estimators. Based on these results, we also compare the true limiting variance to the semiparametric efficiency bound for a general dimension of the covariate vector. We show that the post-matching DiD estimators do not attain the bound when the number of matches is fixed, but approach it as the number of neighbors diverges under some conditions.
We illustrate our results with simulations and with a reanalysis of the National Supported Work (NSW) Demonstration. Using the experimental benchmark from the randomized evaluation following Lalonde_1986_AER, we find that the post-matching 2WFE estimator yields estimates that are much closer to the experimental estimates than the unmatched 2WFE, and accounting for the matching step increases the standard errors.
Finally, while our main results focus on nearest-neighbor matching with replacement, we discuss how they extend to matching on discrete covariates, nearest-neighbor matching without replacement and propensity score matching.
Our paper contributes to a large and growing literature on DiD and 2WFE methods in staggered designs, that is, designs where treated units exhibit variation in their treatment timing Chaise-Dhaut_2022_EJ,Roth-etal_2023_JoE. In particular, studies have found that in staggered designs, 2WFE estimators generally recover linear combinations of ATTs with weights that are hard to interpret and possibly negative when the treatment effects are heterogeneous Athey-Imbens_2021_JoE,Goodman-Bacon_2021_JoE,Chaise-Dhaut_2020_AER,Imai-Kim_2021_PA,Sun-Abraham_2021_JoE,Borusyak-Jaravel-Spiess_2024_ReStud. These studies focus on unmatched 2WFE estimators, mainly under the unconditional version of the parallel trends assumption. We complement this literature by showing that these issues remain for post-matching 2WFE estimators when the parallel trends assumption holds conditionally. Furthermore, we show that pooling different treatment cohorts in the matching step creates an additional bias. Specifically, our Theorem (ref) decomposes the probability limit of the post-matching 2WFE estimator into a linear combination of cohort-specific ATTs and three types of bias: one due to the heterogeneity of treatment effects over time, one due to the heterogeneity of treatment effects across cohorts and a bias due to differences in covariate distributions across different treated cohorts. Because the third bias, which is specific to our setting, is caused by differences in covariate distributions and not by treatment effect heterogeneity, this bias remains even when the treatment effect is constant.
Another closely related strand of the literature is the one analyzing DiD and 2WFE under a conditional parallel trends assumption Abadie_2005_ReStud,SantAnna-Zhao_2020_JoE,Callaway-SantAnna_2021_JoE. These studies focus on inverse-probability weighting and doubly-robust estimators. We build on their framework to analyze alternative estimators that control for covariates employing fully non-parametric nearest-neighbor matching methods.
Our paper also contributes to the literature on cross-sectional matching methods Heckman-Ichimura-Todd_1998_ReStud,Abadie-Imbens_2006_ECMA,Abadie-Imbens_2008_ECMA,Abadie-Imbens_2011_JBES,Abadie-Imbens_2012_JASA,Abadie-Imbens_2016_ECMA,Abadie-Spiess_2022_JASA by generalizing the results to 2WFE estimators with panel data. Specifically, our Propositions (ref) and (ref) show that while the post-matching 2WFE that pools all treated cohorts is a complicated linear combination of 2x2 DiD estimators, the estimator that compares one treated cohort to the never-treated can be recast as a cross-sectional matching estimator using a simple transformation of the outcome of interest, and thus existing matching methods can be used to characterize the probability limit and asymptotic distribution of this estimator. In addition, our closed-form characterization of the asymptotic variance based on the results by Chen-Han_2024_wp allows us to compare it to the semiparametric efficiency bound for a general dimension of the covariate vector, a result previously available only for the scalar case. In line with existing results for the average treatment effect (ATE) in cross-sectional settings Abadie-Imbens_2006_ECMA,Lin-Ding-Han_2021_ECMA, we find that the asymptotic variance is larger than the semiparametric efficiency bound, but approaches it under some conditions as the number of neighbors diverges.
Finally, among the papers at the intersection of matching and panel data methods, Heckman-Ichimura-Todd_1997_ReStud propose a DiD-matching estimator in a two-period setting and derive its limiting distribution when matching is conducted using kernel and local linear regression methods. We complement their findings by considering nearest-neighbor matching with multiple periods and a staggered design, a setting that is very common in empirical work. Imai-Kim-Wang_2023_AJPS propose matching methods for panel data, matching units based on their treatment history. They also show that their estimators can be written as weighted 2WFE regressions, but their inference is conducted conditionally on the matching weights and does not account for the variability introduced by the matching step.
The remainder of the paper is organized as follows. Section (ref) describes the setup, notation and provides identification conditions for the parameters of interest. Section (ref) characterizes the estimators of interest and their probability limits. Section (ref) studies the asymptotic distribution of the estimators, provides a closed-form formula for their limiting variance, discusses efficiency and the inconsistency of the naive cluster-robust variance estimator that ignores the matching step. Section (ref) contains simulations studies and Section (ref) illustrates our results using data from the NSW program. Section (ref) discusses extensions of our results to other types of matching, and Section (ref) provides concluding remarks.
Consider panel data following $n$ units $i=1,\ldots,n$ over $T$ time periods $t=1,2,\ldots,T$. Let $t_i^*$ be a random variable with support $\mathcal{S}=\{s\in\{2,\ldots,T,\infty\}:\mathbb{P}[t_i^*=s]>0\}$ that indicates the period in which unit $i$ receives the treatment for the first time. The values of $t_i^*$ identify different treatment cohorts. We assume that there is at least one pure pre-treatment period (and possibly more) in which no unit is treated, $\mathbb{P}[t_i^*=1]=0$, and denote the never-treated units by $t_i^*=\infty$. We denote the population proportion of each treatment cohort by $p_s=\mathbb{P}[t_i^*=s]$. We assume that the treatment is an absorbing state, so that once unit $i$ is treated, it remains treated through period $T$. Let $D_{it}=\mathbbm{1}(t\ge t_i^*)$ be the treatment indicator in each period.
For $t'>t$, let $\mathbf{d}_{t:t'}=(d_t,d_{t+1},\ldots,d_{t'})$ be a vector of treatments between periods $t$ and $t'$, where $d_{t}\in\{0,1\}$ and $d_{t+1}\ge d_t$ for all $t$, and let $\mathbf{0}_{t:t'}$ and $\mathbf{1}_{t:t'}$ be vectors with all their entries being equal to zero or one, respectively. We use the shorthand notation $\mathbf{d}_{1:t}=\mathbf{d}_t$. Without further assumptions, the potential outcome for each unit in each time period may be a function of the whole treatment path, $Y_{it}(\mathbf{d}_T)$. We introduce the following assumption, standard in the DiD literature, which ensures that potential outcomes do not depend on future treatments.
We use the shorthand notation $Y_{it}(\mathbf{0}_t)=Y_{it}(0)$ to denote the potential outcome of a unit not receiving treatment up to period $t$ and, for any $t\ge t_i^*$, we write $Y_{it}(t_i^*)=Y_{it}(\mathbf{0}_{t_i^*-1},\mathbf{1}_{t^*_i:t})$, which is the potential outcome in period $t$ for a unit that enters treatment in period $t_i^*$. The unit-level treatment effect is defined as $\tau_{it}=Y_{it}(t_i^*)-Y_{it}(0)$. The main parameters of interest are the cohort-specific ATTs, $\mathbb{E}[\tau_{it}|t_i^*=s]$.
Let $X_i$ be a time-invariant vector of covariates of dimension $q\ge 1$, measured at a pure pre-treatment period $t<\min\{t_i^*\}$. To justify the matching procedure, we assume that the parallel trends assumption holds conditional on covariates in the following way.
This assumption states that for each value of $X_i$, all treatment cohorts experience the same outcome trend in the absence of treatment.
Finally, let $e_s(X_i)=\mathbb{P}[t_i^*=s|X_i]$ be the cohort-specific propensity score, that is, the probability of entering treatment in period $s$ given covariates $X_i$. We impose the following conditions which guarantees a non-zero mass of never treated and common support.
The following identification result is standard in the literature, and we include it here for reference.
Consider a sample $(Y_{i1},\ldots,Y_{iT},t_i^*,X_i')_{i=1}^n$. Let $N_s=\sum_i \mathbbm{1}(t_i^*=s)$ be the number of units in treatment cohort $t_i^*=s<\infty$ and $N_0=\sum_i\mathbbm{1}(t_i^*=\infty)$ be the number of never-treated units in the sample, with corresponding sample proportions $\hat{p}_s=N_s/n$ and $\hat{p}_\infty=N_0/n$. We assume observations are drawn independently from the same distribution.
We study the following two-step estimation strategy commonly used in empirical practice. In the first step, each treated unit is matched with replacement to $M\ge1$ units from the pool of never-treated units based on their covariates $X_i$. In the second step, treatment effects are estimated using a 2WFE regression with weights defined by the matching procedure, as described in detail below.
To analyze the matching step, we rely on the setup from Abadie-Imbens_2006_ECMA. As in their paper, we consider the case in which all the covariates are continuous, and discuss matching on discrete covariates in Section (ref). Let $\left\lVertx\right\rVert=\sqrt{x'x}$. For each unit $i$, we define $j_m(i,s')$ as the index $j\in\{1,\ldots,n\}$ of units in comparison cohort $s'\ne t_i^*$ such that: \[\sum_{l:t_l^*=s'}\mathbbm{1}\left\{\left\lVertX_l-X_i\right\rVert\le \left\lVertX_j-X_i\right\rVert\right\}=m,\] that is, the index of the $m$-th closest unit to $i$ in covariate values among the comparison cohort $s'$. Let $\mathcal{J}_M(i,s')=\{j_1(i,s'),j_2(i,s'),\ldots,j_M(i,s')\}$ be the set of indices for the $M$ nearest neighbors in cohort $s'$ for unit $i$. Finally, let: \[K_M(i,s)=\sum_{l:t_l^*=s}\mathbbm{1}\{i\in\mathcal{J}_M(l,t_i^*)\}\] be the number of times unit $i$ is used as a match for units in cohort $s$. We also let: \[K_M(i)=\sum_{s\ne t_i^*}K_M(i,s)\] the total number of times unit $i$ is used as a match for units in any other cohort.
Once the matching is conducted, treatment effects are estimated using a 2WFE regression on the matched sample that includes the treated units and their nearest neighbors. We consider the following general weighted 2WFE estimator. Let $w_i$ be (possibly random) unit-level weights, collected in a vector $w=(w_1,\ldots,w_n)$. For each time period $t=1,\ldots,T$, let $\lambda_t\in\{0,1\}$ be a nonrandom indicator of whether the period is included in the analysis, and collect these indicators in a vector $\lambda=(\lambda_1,\ldots,\lambda_T)$ where $\sum_t\lambda_t\ge 2$ so there are at least two periods. The weighted 2WFE estimator is defined as the estimator of $\tau$ from the 2WFE regression:\footnote{We assume that the weights $(w,\lambda)$ are chosen in such a way that this estimator is well defined, which requires avoiding perfect collinearity between the regressor $D_{it}$ and the unit and time effects. This rules out obvious cases where, for example, no included cohort experiences a treatment change between the included periods, or all included cohorts experience the same treatment change.}
using weights $w_i\lambda_t$, and can be written as:
where \[\bar{D}_i=\frac{\sum_t\lambda_tD_{it}}{\sum_t\lambda_t},\quad \tilde{D}_t=\frac{\sum_iw_iD_{it}}{\sum_iw_i},\quad \bar{D}=\frac{\sum_i\sum_tw_i\lambda_tD_{it}}{\sum_iw_i\sum_t\lambda_t}.\] Lemma (ref) in the appendix provides a general representation of this estimator for an arbitrary set of weights.
In what follows, we focus on the empirical approach described above, which assigns a weight equal to one to all treated units, pooling all treatment cohorts, and reweight the never-treated units by the number of times they are used as a match, scaled by the number of neighbors $M$. We refer to this estimator as the full-sample matched 2WFE estimator. The following result characterizes this estimator.
Proposition (ref) shows that the full-sample matched 2WFE estimator can be decomposed as a linear combination of all possible 2x2 estimators, comparing each cohort $s$ to the never-treated cohort in different post- and pre-treatment periods, $\hat\tau^s_{s'}(t,t')$, and comparing each treated cohort $s$ to other eventually treated cohorts, $\hat\tau^s_{s'}(t,t')$ for $s'\ne\infty$. As pointed out in the literature on unmatched 2WFE estimators Chaise-Dhaut_2020_AER,Athey-Imbens_2021_JoE,Goodman-Bacon_2021_JoE,Imai-Kim_2021_PA,Sun-Abraham_2021_JoE,Borusyak-Jaravel-Spiess_2024_ReStud, the second set of comparisons is problematic because it involves some comparison groups that are already treated, often known as “forbidden comparisons”. Importantly, while the 2x2 estimators using the never-treated, $\hat\tau^s_\infty(t,t')$, use the matching weights, the 2x2 estimators comparing different eventually treated cohorts, $\hat\tau^s_{s'}(t,t')$ for $s'\ne\infty$, are unweighted DiD comparisons. These unweighted comparisons generally bias the estimator, as we show next. To characterize the probability limit of the estimator, we introduce the following regularity conditions.
Theorem (ref) shows that in a staggered design, the post-matching 2WFE estimator does not generally recover a causally interpretable parameter. The probability limit of this estimator can be decomposed into four components. The first component, in the first line of the expression for $\hat\tau(w^\mathsf{pool},\lambda)$ above, is a weighted average of ATTs across cohorts and over the included time periods, where the weights are non-negative and sum up to one, $\sum_{s<\infty}\phi_1(s,\lambda)+\sum_{s<\infty}\sum_{s<s'<\infty}\phi_2(s,s',\lambda)=1$. Although this is a proper weighted average, the weights may be highly unequal across ATTs, and this term is generally hard to interpret.
In addition to this weighted average of ATTs, the full-sample matched 2WFE estimator contains three bias terms. The first bias, $B^\mathsf{het}_\mathsf{cohort}$, involves differences in ATTs in a given time period between different cohorts, $\mathbb{E}[\tau_{it}|t_i^*=s]-\mathbb{E}[\tau_{it}|t_i^*=s']$. These terms are non-zero whenever treatment effects are heterogeneous across cohorts. Notice that because this term consists of differences in ATTs, some ATTs enter negatively. This is an illustration of the negative weights issue pointed out in the literature on unmatched DiD with staggered adoption.
The second bias, $B^\mathsf{het}_\mathsf{time}$, involves the change in the ATTs for a given cohort over time, $\mathbb{E}[\tau_{it}-\tau_{it'}|t_i^*=s']$. This term appears because the 2WFE estimator uses already-treated units as comparison units, resulting in some ATTs entering with a negative sign in the probability limit of the estimator. Notice that the bias terms $B^\mathsf{het}_\mathsf{cohort}$ and $B^\mathsf{het}_\mathsf{time}$ have been discussed in existing literature analyzing unmatched DiD estimators. Theorem (ref) shows that these issues persist when using matched 2WFE estimators.
The third bias, $B_\mathsf{pool}$, on the other hand, is specific to the post-matching setting when the matching is conducted as explained above. This “pooling bias” involves differences in untreated outcome trends across eventually-treated cohorts, $B_{s'}^s(t,t')$. To understand the intuition behind this bias, notice that these terms can be rewritten as: \[B_{s'}^s(t,t')=\mathbb{E}\left\{\mathbb{E}\left[Y_{it}(0)-Y_{it'}(0)|t_i^*=s,X_i\right]\left(\frac{f_s(X_i)-f_{s'}(X_i)}{f(X_i)}\right)\right\}\] where $f(x)$ is the unconditional density of $X_i$ and $f_s(x)$ and $f_{s'}(x)$ are the conditional densities of $X_i$ for cohorts $s$ and $s'$, respectively. Thus, this bias is driven by the difference in covariate distributions across eventually-treated cohorts. As shown in Proposition (ref), the full-sample matched 2WFE estimator involves comparisons between different treated cohorts. Although each treated unit is matched in the first step to a never-treated unit, ensuring that the distribution of covariates for the treated is identical to the reweighted distribution of the never-treated, the different treated cohorts are not matched to each other, and thus their covariate distributions will differ in general. Because the parallel trends assumption only holds conditionally on $X_i$, different treated cohorts may have different unconditional trends in their untreated outcomes, which creates this bias. In other words, when the parallel trends assumption holds conditionally, any unmatched comparison between eventually treated cohorts is invalid, even those involving not-yet-treated units as comparisons.
When the treatment effect is constant both across units and over time $\tau_{it}=\tau$, $\mathbb{E}[\tau_{it}|t_i^*=s]-\mathbb{E}[\tau_{it}|t_i^*=s']=\mathbb{E}[\tau_{it}-\tau_{it'}|t_i^*=s']=0$, so the biases due to treatment effect heterogeneity equal zero, $B^\mathsf{het}_\mathsf{time}=B^\mathsf{het}_\mathsf{cohort}=0$. However, the bias due to the difference in covariate distributions, $B_\mathsf{pool}$, remains regardless of whether treatment effects are homogeneous or heterogeneous. This bias equals zero when either the parallel trends assumption holds unconditionally, $\mathbb{E}\left[Y_{it}(0)-Y_{it'}(0)|t_i^*=s,X_i\right]=\mathbb{E}\left[Y_{it}(0)-Y_{it'}(0)|t_i^*=s\right]$ for all $s$, so matching is unnecessary, or when the distribution of covariates across cohorts are the same, $f_s(x)=f_{s'}(x)$, which is generally not true even when the treatment effect is constant.
Thus, the commonly used full-sample matched 2WFE estimator requires stringent homogeneity conditions on treatment effects and covariate distributions to recover causally interpretable parameters when there is variation in treatment timing. To avoid this problem, a staggered design can be split into multiple pairwise designs, as suggested by Callaway-SantAnna_2021_JoE,Sun-Abraham_2021_JoE,Steigerwald-VazquezBare-Meier_2021_JAERE, among others, in unmatched panel designs. We now consider an estimator that assigns weight equal to one to units in a specific cohort $t_i^*=s$, weights equal to $K_M(i,s)/M$ to never-treated units, and weight equal to zero to all other cohorts. We refer to these estimators as pairwise matched DiD estimators. The following proposition characterizes this estimator.
Proposition (ref) shows that the pairwise matched DiD estimator can be recast as cross-sectional matching estimator like the ones studied by Abadie-Imbens_2006_ECMA, where the outcome is the difference in average post- and pre-treatment periods. As an illustration, for a given $s$, when $\lambda_t^*=\mathbbm{1}(t\in\{\ell',\ell\})$ for a pair of periods such that $\ell'<s\le \ell$,
which is a 2x2 (i.e. two-group-two-period) matched DiD estimator.
The pairwise matched DiD estimator is consistent for an average of the ATTs for cohort $s$ over the included post-treatment periods $t\ge s$, as we show below.
As an illustration, if $\lambda_t=\mathbbm{1}(t\in\{\ell',\ell\})$ with $\ell'<s\le \ell$, we have that $\hat\tau(w^s,\lambda)\to_\mathbb{P} \mathbb{E}[\tau_{i\ell}|t_i^*=s]$ which is the ATT for cohort $s$ in period $\ell$. On the other hand, when all periods are included, so that $\lambda_t=1$ for all $t$, $\hat\tau(w^s,\lambda)\to_\mathbb{P}\sum_{t\ge s}\mathbb{E}[\tau_{it}|t_i^*=s]/(T+1-s)$ which is a simple average of ATTs over periods $t\ge s$.
These results demonstrate that the inconsistency of the full-sample matched 2WFE estimator can be avoided by comparing each treated cohort to the never treated, which consistently estimates an average of ATTs over time. When conducting inference based on the pairwise matched DiD estimators, the variance needs to account for the variability introduced by the matching step, and failing to account for this variability generally results in invalid standard error estimators. The next section studies the asymptotic distribution of the pairwise matched DiD estimators, proposes valid standard error estimators and characterizes the bias of the naive variance estimators that ignore the matching step.
In this section, we study the asymptotic distribution and valid inference procedures for the pairwise estimators (ref). In what follows, for a given $s\in\mathcal{S}$, let $\Delta \bar{Y}_i^s(\lambda)=\bar{Y}^{\mathsf{post},s}_i(\lambda)-\bar{Y}^{\mathsf{pre},s}_i(\lambda)$, $\mu_\lambda(s,x)=\mathbb{E}[\Delta \bar{Y}_i^s(\lambda)|t_i^*=s,X_i=x]$ and $\sigma^2_\lambda(s,x)=\mathbb{V}[\Delta \bar{Y}_i^s(\lambda)|t_i^*=s,X_i=x]$ where we use the notation $\mu_\lambda(0,x)$ and $\sigma^2_\lambda(0,x)$ when $s=\infty$. Also let: \[\bar{\tau}_i^{\mathsf{post},s}(\lambda)=\frac{\sum_{t\ge s}\lambda_t\tau_{it}}{\sum_{t\ge s}\lambda_t}\] be the average unit-level treatment effect over the included post-treatment periods $t\ge s$ with $s<\infty$. The following result characterizes the asymptotic distribution of the pairwise matched DiD estimator.
Theorem (ref) shows that the pairwise matched DiD estimator is asymptotically normal after appropriate centering and scaling, which means that valid inference can be conducted using standard procedures. In general, the estimator needs to be bias-corrected, as initially pointed out by Abadie-Imbens_2006_ECMA. We refer to this bias $B_\lambda(s,M)$ as the “matching bias”. The matching bias $B_\lambda(s,M)$ is a finite-sample bias due to the fact that, when covariates are continuous, it is not possible to find perfect matches for treated units in any particular sample. This type of finite sample bias is typical in nonparametric estimators. It can be shown that under the assumptions in the theorem, $B_\lambda(s,M)=O_\mathbb{P}(n^{-1/q})$ where $q$ is the dimension of $X_i$. This implies that $B_\lambda(s,M)=o_\mathbb{P}(1)$ so the matching bias does not appear in the probability limit of the estimator, but it needs to be accounted for in the asymptotic distribution because it is multiplied by the convergence rate.
In some cases, the bias correction is not needed. For example, when matching on a single covariate, $\sqrt{n}B_\lambda(s,M)=o_\mathbb{P}(1)$. Abadie-Imbens_2006_ECMA also provide an alternative sampling scheme in which the bias is asymptotically negligible when the never-treated group is “sufficiently larger” than the treated group; see Section (ref) for further discussion. Finally, when matching on discrete covariates, matches are exact, so the bias is zero, as discussed in Section (ref).
For these reasons, the nature of the matching bias is very different from that of the pooling bias of the full-sample matched 2WFE estimator characterized in Theorem (ref), $B_\mathsf{pool}$. Specifically, the pooling bias is a large-sample bias that appears in the probability limit of the estimator, and does not decrease with the sample size.
Because the pairwise estimator is asymptotically normal, valid inference can be conducted using standard methods after estimating the bias and the asymptotic variance that accounts for the matching step. Proposition (ref) shows that the pairwise matched DiD estimator can be written as a cross-sectional estimator, and therefore the methods available in the literature for bias and variance estimation like the ones in Abadie-Imbens_2006_ECMA,Abadie-Imbens_2011_JBES can be applied in our setting after transforming the data appropriately.
The analytic formula for the asymptotic variance we derive in Theorem (ref) allows us to study the efficiency of the pairwise matched DiD estimator for a general dimension of the covariate vector. Under Assumptions (ref), (ref) and (ref), SantAnna-Zhao_2020_JoE derive the semiparametric efficiency bound for the ATT with panel data under the conditional parallel trends assumption. When comparing cohort $s$ against the never-treated, the bound is given by:
Then we have the following result.
Proposition (ref) shows that the matched DiD estimator does not attain the semiparametric efficiency bound, and that the difference between the asymptotic variance and the bound depends on the number of nearest neighbors $M$ and on the dimension of the covariates $q$ through the function $\alpha(M,q)$. As an illustration, with a single covariate $q=1$, $\alpha(M,1)=M(2M+1)/2$ Chen-Han_2024_wp and thus:
In this case it is clear that the difference between the variance and the bound vanishes as $M\to\infty$, in line with the results in Abadie-Imbens_2006_ECMA and Lin-Ding-Han_2021_ECMA for the ATE with cross-sectional data under unconfoundedness. Because the closed form of $\alpha(M,q)$ has not been characterized for a general $M$ and $q$, it is hard to extend this statement for the general case. However, the tabulations in Chen-Han_2024_wp show that for $M\in\{2,3,\ldots10\}$ and $d\in\{2,3,\ldots,10\}$, $\alpha(M,q)/M^2-1<0.15$, and that $\alpha(M,q)/M^2-1$ is positive and decreasing in $M$ for each of these values of $M$ and $q$.
Because the pairwise matched DiD estimator is constructed in two steps, inference needs to account for the first-step estimation in the second-step asymptotic variance, as done in Theorem (ref). In practice, however, researchers often rely on the cluster-robust variance estimator from the second-stage 2WFE regression, failing to account for the variability introduced by the matching step, which results in inconsistent variance estimators. With cross-sectional data under unconfoundedness, Abadie-Imbens_2016_ECMA compare the asymptotic variance of the oracle propensity-score matching ATT estimator with the variance of the ATT estimator using the estimated propensity score, and show that the latter can be smaller or larger depending on the data generating process. On the other hand, Abadie-Spiess_2022_JASA consider post-matching regression estimators using cross-sectional data and matching without replacement and show that standard errors that fail to account for the matching step are inconsistent when the regression of interest is misspecified. In this section, we derive the probability limit of the cluster-robust 2WFE variance estimator that ignores the matching step, which we refer to as the naive cluster-robust variance estimator.
By Proposition (ref), the cluster-robust 2WFE estimator of the variance of $\hat\tau(w^s,\lambda)$ can be written as:
where:
is the regression residual and the dependence on $\lambda$ is left implicit to reduce notation. The following result characterizes the probability limit of this estimator.
Proposition (ref) shows that if $\mathbb{E}[\Delta\bar{Y}_i^s(0,\lambda)|t_i^*=s,X_i]=\mathbb{E}[\Delta\bar{Y}_i^s(0,\lambda)|t_i^*=s]$, so that covariates do not affect the average counterfactual trend, the naive cluster-robust variance estimator that ignores the matching step is consistent for the true variance. This case, however, only occurs when the matching step is unnecessary, because counterfactual trends are parallel unconditionally. More generally, the naive cluster-robust variance estimator is inconsistent for the true asymptotic variance, with an asymptotic bias given by:
The direction and magnitude of this bias depend on the data generating process. The first and third and terms are always positive, whereas the covariance can be positive or negative. For example, if the conditional ATT, $\mathbb{E}[\tau_{it}|t_i^*=s,X_i]$, is constant over $X_i$, the covariance is zero and thus the naive variance estimator is upward biased. In such cases, ignoring the matching step results in overly conservative standard errors. On the other hand, when the covariance between the conditional ATT and the conditional trend of the untreated outcome is negative and large enough, the bias may be negative, resulting in standard errors that underestimate the true variability of the estimator. We further illustrate these cases with simulations in the next section.
We now present two simulation studies to illustrate the problems with post-matching 2WFE estimators analyzed in previous sections. In the first study, we consider a staggered treatment adoption setting to illustrate the asymptotic bias. In the second study, we show how failing to account for the matching step when conducting inference results in invalid standard errors, considering cases where the unadjusted standard errors are too large and too small.
Our first setup consists of panel data with a total of $n=1,000$ units and four time periods, $t=1,2,3,4$. Units belong to one of three cohorts, $t_i^*\in\{2,3,\infty\}$, so the “early adopters” enter treatment in period $t=2$ and the “late adopters” enter treatment in period $t=3$. Units are assigned to cohorts based on a multinomial rule: \[\mathbb{P}[t_i^*=\infty|X_i]=\frac{1}{1+\exp(X_i)+\exp(2X_i)},\quad \mathbb{P}[t_i^*=s|X_i]=\frac{\exp((s-1)X_i)}{1+\exp(X_i)+\exp(2X_i)},\quad s=2,3\] where $ X_i$ is a scalar covariate with $X_i\sim \mathrm{Uniform}(-1,2)$. The untreated potential outcome is defined as: \[Y_{it}(0)=\alpha_i+(t-1)+5X_i(t-1)+u_{it},\quad \alpha_i\sim\mathcal{N}(0,1),\quad u_{it}\sim \mathcal{N}(0,1)\] with $(\alpha_i,u_{it},X_i)$ mutually independent. This implies that $\mathbb{E}[Y_{it}(0)-Y_{it-1}(0)|t_i^*=s,X_i]=1+5X_i$ for all $s$, so the conditional parallel trends assumption holds (the unconditional parallel trends assumption does not because $\mathbb{E}[X_i|t_i^*=s]$ differs across $s$). To focus on the bias due to the matching procedure and abstract from the other biases characterized in Theorem (ref), we assume that the treatment effect is constant, $\tau=5$. The model for the treated potential outcome is: \[Y_{it}(1)=5 + \alpha_i+(t-1)+5X_i(t-1)+v_{it},\quad v_{it}\sim \mathcal{N}(0,1).\] We estimate $\tau$ by matching each treated unit to a never-treated unit in the first step and running Equation (ref) in the reweighted data. The results are displayed in the first column of Table (ref). The results show that the bias of the full sample matched 2WFE is about 0.7 or 15% of the true ATT. In addition, the naive cluster-robust standard error severely overestimates the true variability.
Our second simulation exercise illustrates how failing to account for the matching step leads to invalid inference. We consider panel data with a total of $n=1,000$ units and four time periods, $t=1,2,3,4$. Units belong to one of two cohorts, $t_i^*\in\{3,\infty\}$, so there is a treated cohort that enters treatment in period 3, and a never-treated cohort. Units are assigned to treatment based on a logistic model: \[\mathbb{P}[t_i^*=3|X_i]=\frac{\exp(X_i)}{1+\exp(X_i)}\] where $X_i$ is a scalar covariate with $X_i\sim\mathrm{Uniform}[-1/2,1/2]$. We consider two designs. The first design has a constant treatment effect equal to 5 and the potential outcomes are given by:
where as before $u_{it}\sim \mathcal{N}(0,1)$, $v_{it}\sim \mathcal{N}(0,1)$ and $(\alpha_i,u_{it},X_i)$ mutually independent. Because in this setup the treatment effect is constant, Proposition (ref) shows that the naive cluster-robust standard error is biased upwards.
In the second design, the potential outcomes are given by:
which implies that the conditional ATT is $\mathbb{E}[Y_{it}(1)-Y_{it}(0)|t_i^*=3,X_i]=7X_i(t-1)$ which changes across $X_i$ and over time. On the other hand, $\mathbb{E}[\Delta Y_{it}(0)|X_i,t_i^*=s]=1-2X_i$ and $\mathbb{C}\mathrm{ov}(\mathbb{E}[Y_{it}(1)-Y_{it}(0)|t_i^*=3,X_i],\mathbb{E}[\Delta Y_{it}(0)|X_i,t_i^*=s]|t_i^*=s)=-14(t-1)\mathbb{V}[X_i]$ which in this case results in downward-biased standard errors.
The results from these simulations are shown Table (ref). Column (2) reports the result for the constant treatment effect scenario. As shown in Proposition (ref), when the treatment effect is constant, the naive variance estimator overestimates the true variability. In this case, the Monte Carlo (MC) standard error is 0.086, while the average of the naive standard error across simulations is 0.257, about three times larger. This upward bias results in overcoverage of confidence intervals and loss of power. On the other hand, the standard error estimator that accounts for the variability of the matching procedure is 0.085, almost identical to the MC standard error, and reaches correct coverage.
Column (3) reports the results for the case with a heterogeneous treatment effect and a negative covariance between the conditional ATT and the conditional outcome trend under no treatment, the third term in Proposition (ref). In this case, the negative covariance term yields a downward-biased standard error of 0.209, about 13% smaller than the MC standard error (0.239) and the naive CI undercovers the true parameter. As before, the standard error that accounts for the matching procedure gives is very close to the MC standard error and reaches correct coverage.
We now illustrate our results using data from the National Supported Work (NSW) Demonstration. Using the NSW experimental data as a benchmark, Lalonde_1986_AER compared different non-experimental methods to determine whether they can replicate the experimental findings. This study sparked an extensive debate about the credibility of non-experimental methods Imbens-Xu_2025_JEL.
Following Lalonde_1986_AER, we use data from the Current Population Survey (CPS) to construct a comparison group for the experimental treatment group using nearest-neighbor matching. The outcome of interest is earnings, which is measured in 1974, 1975 (before treatment) and 1978 (after treatment). Our list of covariates includes indicators of age at baseline, indicators of years of education and indicators of being married, Black and Hispanic.
We compare the following specifications using earnings data from 1975 and 1978, which constitutes a two-period panel: (1) a difference in mean earnings between treated and controls in the experimental sample (“Experimental”), (2) an unmatched 2WFE regression (“2WFE”), (3) a naive matched 2WFE regression where standard errors do not account for the matching step (“Naive matched 2WFE”), (4) a matched 2WFE regression calculating valid standard errors (“Matched 2WFE”), (5) a matched 2WFE regression with valid standard errors and bias correction using a linear specification for the covariates (“Matched 2WFE - BC”) and (6) a 2WFE regression on the matched sample, that is, on the sample excluding any untreated units that are not matched to any treated unit, but without accounting for the number of times each untreated is used as a match (“2WFE in matched sample”). All standard errors are clustered at the individual level.
The results are shown in Column (1) of Table (ref). The experimental estimate of the ATT is \$886.3 with a p-value of 0.069. We use this value as a benchmark. The unmatched 2WFE regression yields an estimate of about \$1,714, almost twice as big as the experimental benchmark, with a p-value close to zero. The weighted 2WFE that uses matching in the first step gives an estimate of around \$930 (p-value 0.117), much closer to the experimental benchmark. As shown in Proposition (ref), the standard error estimator in this specification is inconsistent because it does not account for the matching step. When adjusting for the matching first step, the standard error increases from \$593 to \$635. The bias-corrected matched 2WFE specification reports a slightly lower estimate of \$809.21 with a very similar standard error. Finally, the 2WFE regression dropping unmatched units reveals an estimate of about \$418, approximately half the experimental benchmark.
Column (2) of Table (ref) repeats the same analysis for a placebo outcome. In the experimental sample, this corresponds to a difference in mean earnings in 1975 (before the treatment) between treated and controls. For all other specifications, the placebo test uses data on earnings in 1974 and 1975, and can be interpreted as a pre-trends test. Because all the data in these specifications is pre-treatment, all estimated effects should be close to zero and insignificant. While all the estimates are insignificant at the usual levels, the one from the unweighted 2WFE and the one from the regression that ignores repeated matches show the largest magnitudes and lowest p-values.
When the vector $X_i$ is a low-dimensional vector of discrete covariates, the matching can be done exactly, which eliminates the matching bias. Without loss of generality, suppose that $X_i$ is a scalar taking values in a finite set $\mathcal{X}$. For such cases, a consistent and asymptotically normal estimator is given by:
where \[\hat{e}_s(X_i)=\hat{\mathbb{P}}[t_i^*=s|X_i=x]=\frac{\sum_i \mathbbm{1}(t_i^*=s)\mathbbm{1}(X_i=x)}{\sum_i\mathbbm{1}(X_i=x)}.\] This estimator that reweights never-treated units by the ratio of the estimated propensity scores is equivalent to the matching DiD estimator (ref) that matches each treated unit to all never-treated units in the same covariate cell $X_i=x$. Because in this setting the number of units in each covariate cell increases as $n\to\infty$, this estimator can be seen as a matching estimator with a diverging number of matches, similar to the ATE estimator considered in a cross-sectional setting under unconfoundedness by Lin-Ding-Han_2021_ECMA. This estimator attains the semiparametric efficiency bound (ref), as shown in the following result.
Matching is sometimes carried out without replacement. In such cases, whenever an untreated unit is matched to a treated unit, that unit is removed from the pool of untreated. This implies that each treated unit is used at most once, $K_M(i)\in\{0,1\}$, and the size of the matched sample is $(M+1)N_s$. We now discuss how to adapt our results to the case of matching without replacement, building on the results by Abadie-Imbens_2012_JASA and Abadie-Spiess_2022_JASA. For simplicity we consider the case where the population consists of only two cohorts, $t_i^*\in\{s,\infty\}$. In this case, the post-matching 2WFE estimator using matching without replacement is:
The following proposition summarizes the behavior of the ATT estimator under matching without replacement.
The conditional sampling scheme in Proposition (ref) and the requirement that $N_s=O(N_0^{1/r})$, which implies that the untreated group is “sufficiently larger” than the treated group, ensures that the bias of the matching estimator is $o_\mathbb{P}(n^{-1/2})$ Abadie-Imbens_2006_ECMA and thus bias correction is not required for asymptotic normality. This is not specific to matching without replacement, and can also be applied to the results in the previous sections of the paper.
When the number of (continuous) covariates is large, matching directly on the vector of covariates can be computationally cumbersome. In such cases, researchers often resort to matching on the propensity score, which reduces the problem to matching on a single covariate. Rosenbaum-Rubin_1983_BIO showed that the propensity score is a balancing score, and thus matching on the propensity score instead of the full vector of covariates is sufficient for eliminating confounding bias under a selection on observables assumption. In practice, the propensity score is unknown and needs to be estimated, which introduces additional variability that needs to be accounted for when conducting inference. Abadie-Imbens_2016_ECMA derive the asymptotic distribution of the cross-sectional propensity score matching estimator. For the sake of brevity, we do not elaborate on the technical aspects of this case, but we note that the results in Abadie-Imbens_2016_ECMA can be adapted to our setting to derive the asymptotic distribution of our pairwse matched DiD estimator when matching on the estimated propensity score.
We analyze the commonly employed empirical strategy of estimating 2WFE regressions in a reweighted sample where treated units are matched to never-treated units in a first step. This approach is often used when the parallel trends assumption is believed to hold conditionally on covariates. By formally characterizing the post-matching 2WFE estimators and their asymptotic distribution, we highlight two problems with this strategy that generally invalidates its findings. First, when treated units exhibit variation in treatment timing, the post-matching 2WFE estimator that pools all treated units contains an asymptotic bias that remains even when the treatment effect is constant. The reason for this bias is that the 2WFE estimator implicitly involves comparisons between different treatment cohorts, using one of them as a comparison group. Because different treated cohorts are not matched to each other, they generally have different covariate distributions, and thus comparing their outcome evolutions over time is invalid when the parallel trends assumption holds conditionally, even when the comparison group is not yet treated.
Second, ignoring the matching step when conducting inference based on the post-matching 2WFE estimators results in inconsistent variance estimators, which can be biased upwards or downwards depending on the data generating process. Specifically, the naive variance estimator is upward-biased (conservative) when the treatment effect is constant, but can be downward-biased when there is sufficient treatment effect heterogeneity.
Then, we provide conditions under which the post-matching 2WFE estimator that compares each treated cohort to the never treated separately is consistent for an average of cohort-specific ATTs over time, and derive a closed-form representation of its asymptotic variance which can be used to conduct valid inference post-matching. We also show that this estimator does not generally attain the semiparametric efficiency bound, but it approaches it as the number of neighbors grows. We also illustrate our results using simulations and a reanalysis of the data from the NSW job training program.
While different methods are available to incorporate covariates when the parallel trends assumption holds conditionally, matching-based methods have the advantage of being intuitively appealing, easy to implement, computationally stable and fully nonparametric. Our findings enable empirical researchers to exploit these methods and conduct valid inference on the parameters of interest. Moreover, because we show that post-matching 2WFE estimators can be recast as cross-sectional matching estimators, existing statistical software for matching estimators such as the teffects nnmatch package in Stata can be used after appropriately transforming the data.