The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
105,827 characters
Matched Triple-Differences: A Framework for Covariate Adjustment
\maketitle
\vspace{-1em}
\begin{abstract}
A common empirical strategy in triple-differences (DDD) is to include either covariate trends or covariate levels to a three-way fixed effects (3WFE) regression. This strategy is typically motivated by the conditional parallel gaps assumption which assumes that deviations from parallel trends are similar among units with comparable observed covariates. We formally study both 3WFE specifications and show that, in general, neither consistently estimates the average treatment effect on the treated (ATT). Our diagnosis has two parts. First, we show that the OLS estimands of both specifications fail to satisfy the covariate balancing condition that is sufficient for them to equal the ATT. Second, we characterize the additional assumptions for each OLS estimand to equal the ATT. To address these limitations, we propose a matched triple-differences framework that accommodates a general class of matching procedures, including nearest-neighbor matching and kernel matching. Within this framework, we construct a class of consistent matching estimators, establish their asymptotic normality, and provide a consistent variance estimator that accounts for the variability introduced by the matching step. Our empirical application shows that the proposed matched triple-differences estimator can yield estimates and qualitative conclusions that differ from those obtained using the 3WFE regression.
\end{abstract} \vspace{1em}
\noindent \textbf{Keywords:} Triple-Differences; Difference-in-Differences-in-Differences; Parallel Gaps Assumption; Matching \\\noindent\textbf{JEL Classification:} C14, C21, C23
\newpage
\bigskip
\noindent
\section{Introduction}
Difference-in-differences (DiD) is one of the most widely used causal inference methods in empirical economics, and its identification relies on the parallel trends assumption. In practice, researchers may worry that this assumption is too restrictive to provide a reasonable approximation to reality. This concern is particularly relevant when treatment requires a unit to satisfy two criteria simultaneously: being geographically exposed to a policy and belonging to the subpopulation eligible for it. In a DiD comparing the treated subgroup (units that are both exposed and eligible) with the unexposed eligible subgroup, geography-specific shocks may violate parallel trends (top two boxes in Figure \ref{figure:foursubgroupintro}). Similarly, subgroup-specific shocks may invalidate a DiD comparing the treated subgroup with the exposed ineligible subgroup (left two boxes in Figure \ref{figure:foursubgroupintro}). Moreover, both types of shocks may affect units differently depending on their characteristics. These concerns motivate triple-differences (DDD) designs with covariates. Their identifying assumption, the conditional parallel gaps assumption, accommodates both geography- and subgroup-specific shocks and allows the effects of these shocks to vary with covariates.
To implement these designs, researchers typically incorporate covariates through one of two three-way fixed effects (3WFE) regression specifications: one includes covariate-specific trends, and the other includes covariate levels. We refer to these as covariate-trend and covariate-level 3WFE regressions, respectively. These specifications are widely used in empirical work: Table \ref{table:ddd_literature} lists $18$ papers published between $2010$ and $2026$ in top economics journals (e.g., AER and QJE) that use at least one of the two specifications. Yet these specifications have received limited attention in the methodological and theoretical literature. In this paper, we formally analyze both specifications, characterize their limitations, and propose valid alternative procedures for estimation and inference based on matching.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[scale=0.9, transform shape,
box/.style={draw, rounded corners, minimum width=4cm, minimum height=1.5cm,
align=center, font=\small},
treated/.style={box, very thick, fill=gray!15, font=\small\bfseries},
lab/.style={font=\small\bfseries},
did/.style={<->, >=stealth, thick, shorten <=2pt, shorten >=2pt},
shock/.style={font=\footnotesize\itshape, text=red!70!black, align=center}
]
\node[lab] at (-3.8, 3.05) {Exposed};
\node[lab] at ( 3.8, 3.05) {Unexposed};
\node[lab, rotate=90] at (-6.5, 1.85) {Eligible};
\node[lab, rotate=90] at (-6.5, -1.85) {Ineligible};
\node[treated] (A) at (-3.8, 1.85) {Exposed Eligible\\(Treated)};
\node[box] (B) at ( 3.8, 1.85) {Unexposed Eligible};
\node[box] (C) at (-3.8, -1.85) {Exposed Ineligible};
\node[box] (D) at ( 3.8, -1.85) {Unexposed Ineligible};
\draw[did] (A) -- node[above, lab] {DiD 1}
node[below, shock] {geography shock} (B);
\draw[did] (A) -- node[left, lab] {DiD 2}
node[right, shock] {subgroup shock} (C);
\end{tikzpicture}
\caption{Limitations of the two DiD comparisons}
\label{figure:foursubgroupintro}
\end{figure}
Throughout, we consider a triple-differences design with multiple periods and simultaneous adoption of treatment, and the parameters of interest are the average treatment effects on the treated subgroup (ATT) in each post-treatment period. Our starting point is a weighting representation of the identification result: the ATT can be expressed as a weighting estimand in which the treated subgroup is left unweighted and each of the three untreated subgroups is reweighted by a density ratio. These density ratios align the covariate distribution of each untreated subgroup with that of the treated subgroup. This representation implies a sufficient covariate balancing condition: a weighting estimand equals the ATT if its weighting profile (i)~balances the covariate distribution of \emph{each} untreated subgroup with that of the treated subgroup, and (ii)~attains this balance while leaving the treated subgroup unweighted.
We then characterize the limitations of the commonly used covariate-trend and covariate-level 3WFE regressions. For the covariate-trend specification, the corresponding OLS estimand can be expressed as a weighting estimand that reweights all four subgroups with weights proportional to linear projection residuals, known in the literature as ``implied regression weights" \citep{Chatt_Zubi_2023_biometrika}. This weighting profile fails to satisfy the covariate balancing condition in three repects. First, the pooling issue, it targets covariate balance between the treated subgroup and a pooled untreated group merging the three untreated subgroups, rather than balancing each untreated subgroup separately with the treated subgroup. Second, the mean-only balance issue, the implied regression weights attain only mean balance rather than balance of the full covariate distribution: two subgroups can still have different covariate distributions even when their covariate means are balanced. Third, the treated-reweighting issue, these weights also reweight the treated subgroup. The covariate-trend OLS estimand equals the ATT if treatment effects are homogeneous in the covariates and, in addition, either the conditional mean functions of the three untreated subgroups are linear with a common slope or the generalized propensity scores take a specific inverse-linear form.
For the covariate-level specification, the corresponding OLS estimand equals the unconditional triple-differences estimand, which assigns a weight of one to all four subgroups. This weighting profile also fails to satisfy the covariate balancing condition because it leaves all three untreated subgroups unweighted. The covariate-level OLS estimand equals the ATT if and only if the parallel gaps assumption holds unconditionally, in which case there is no need to include covariates in the triple-differences design.
To address these limitations, we propose a matched triple-differences framework for valid estimation and inference. The estimation part of our framework is a two-step procedure. The first step conducts three pairwise matchings, each between the treated subgroup and one of the three untreated subgroups. The second step estimates a weighted 3WFE regression using the weights generated in the matching step. The resulting matched triple-differences estimator is designed to satisfy the covariate balancing condition: matching reweights each untreated subgroup toward the covariate distribution of the treated subgroup while leaving the treated subgroup unweighted. The inference part of our framework is based on a set of novel regularity conditions on the matching procedures. These conditions allow us to abstract from any specific matching method and establish asymptotic theory for a general class of matching estimators. Notably, this class includes nearest-neighbor matching, with either a fixed or a diverging number of neighbors, as well as kernel matching. To account for the uncertainty introduced by the matching step, we propose a unified variance estimator that is consistent for every matching procedure covered by our framework.
We illustrate our results with simulations and an empirical application. In our empirical application, we revisit \cite{Cai_2016_AEJ}, who studies the effect of a tobacco insurance program on the borrowing and saving behavior of the treated households. \cite{Cai_2016_AEJ} uses triple-differences as the main empirical design, and her main regression specification is similar to the covariate-level 3WFE regression. We find that our matched triple-differences estimates can differ in magnitude from those produced by her regression specification and lead to qualitatively different results for some outcomes of interest. Specifically, the estimates from her regression specification imply that the program’s effect on net saving switches sign over time (positive in years $0$ and $3$ and negative in year $5$), whereas our matched triple-differences estimates show no
significant effect in any post-treatment year.
\subsection{Related Literature and Contribution}
This paper contributes to four strands of literature. First, it contributes to the literature on triple-differences designs and the 3WFE method. Triple-differences was first proposed by \citet{Gruber_1994_AER}. \citet{Olden_Moen_2022_EJ} formalize the triple-differences estimator in a two-period setting without covariates. Recent work has begun to examine the validity of 3WFE methods and to study triple-differences designs with covariates. Three closely related papers are \citet{Strezhnev_2023_wp}, \citet{Leventer_2025_wp}, and \citet{Ortiz-Villavicencio_SantAnna_2025_working}. \citet{Strezhnev_2023_wp} studies the limitations of 3WFE regressions \emph{without} covariates under staggered adoption, where identification relies on an unconditional parallel gaps assumption. Our paper instead focuses on 3WFE regressions with covariates, where identification relies on a conditional parallel gaps assumption. We show that these regressions can fail even under simultaneous adoption, a setting in which 3WFE without covariates is valid under the unconditional parallel gaps assumption.
Both \citet{Leventer_2025_wp} and \citet{Ortiz-Villavicencio_SantAnna_2025_working} study triple-differences designs with covariates. \citet{Leventer_2025_wp} examines an empirical practice distinct from the 3WFE approach studied in our paper. Specifically, the paper shows that treating triple differences with covariates as the difference between two conditional DiD comparisons is invalid. \citet{Ortiz-Villavicencio_SantAnna_2025_working} also examine the practice of differencing two conditional DiD comparisons. In addition, they provide Monte Carlo evidence that the two 3WFE specifications considered in our paper can be biased in a two-period setting. Our paper differs from and complements these studies in two respects. First, we provide an analytical diagnosis that characterizes the distinct limitations of the two commonly used 3WFE specifications with covariates. This diagnosis explains the specific reasons why each specification can fail and characterizes the restrictions under which it is valid, thereby complementing the simulation evidence in \citet{Ortiz-Villavicencio_SantAnna_2025_working}. Second, the estimators proposed by \citet{Leventer_2025_wp} and \citet{Ortiz-Villavicencio_SantAnna_2025_working} are based on inverse-probability weighting and doubly robust estimation, whereas our proposed estimator is based on fully nonparametric matching methods.
Second, this paper contributes to the literature on cross-sectional matching methods. Earlier work focuses on asymptotic theory and inference for specific matching methods, such as nearest-neighbor matching (\citeinb{Abadie-Imbens_2006_ECMA}; \citeinb{Lin-Ding-Han_2023_ECMA}; \citeinb{Abadie-Spiess_2022_JASA}) and kernel matching (\citeinb{Heckman-etal_1998_ECMA}). More recently, \cite{Lin_Han_2025_JOE} and \cite{Meng_et_al_2026_wp} have developed general asymptotic theory and inference procedures that cover a class of matching methods. As in these two papers, our asymptotic framework also abstracts from a specific matching method but differs from them in two aspects. First, our framework is based on regularity conditions on the matching procedure that are non-nested with those in \cite{Lin_Han_2025_JOE} and \cite{Meng_et_al_2026_wp}. Second, we develop a different variance estimator. Our variance estimator builds on that of \cite{Abadie-Imbens_2006_ECMA} but extends it to the triple-differences setting and to a class of matching procedures, which requires a different proof strategy.
Third, this paper adds to the line of work that combines matching with panel data, which has primarily focused on the difference-in-differences setting. \citet{Heckman-Ichimura-Todd_1997_ReStud} propose a difference-in-differences matching estimator for a two-period setting based on kernel matching and local linear regression; \citet{Liu_Vazquez_Bare_2026_wp} study nearest-neighbor matching combined with difference-in-differences designs in settings with multiple periods and staggered adoption. We extend this literature to triple-differences designs: we lay out how to combine matching with triple-differences designs and develop asymptotic theory and inference for a class of matching methods.
Finally, our work is broadly related to the literature on interpreting OLS estimands in regression specifications with covariates \citep{Sloczynski_2022_Restat,Sloczynski_2026_Restud,Blandol_et_al_2026_Restud}. Specifically, the implied regression weights of the covariate-trend specification are closely related to those studied by \citet{Chatt_Zubi_2023_biometrika} in the cross-sectional setting and \citet{Caetano-Callaway_2026_wp} in the difference-in-differences setting. However, we show that the implied regression weights involve a pooling issue that has no counterpart in either the cross-sectional or the simultaneous-adoption DiD setting.
The remainder of the paper is organized as follows. Section \ref{section:econometricssetup} presents the motivating
example, introduces the setup and notation, and establishes the identification result and the covariate balancing condition. Section \ref{section:empiricaldiagnose} diagnoses the two 3WFE regression specifications commonly used in empirical practice. Section \ref{section:matchedtriplediffereces} introduces the matched triple-differences framework, establishes the consistency and asymptotic normality of the resulting estimators, and develops the variance estimator for inference. Section \ref{section:simulation} presents simulation studies, Section \ref{section:empiricalapplication} illustrates our methods by revisiting \citet{Cai_2016_AEJ}, and Section \ref{section:conclusion} concludes.
\section{Setup, Estimand, and Identification}\label{section:econometricssetup}
\subsection{Motivating Example}
As a motivating example, we use the study by \citet{Cai_2016_AEJ} to illustrate a triple-differences design and explain how conditioning on covariates can make its identifying assumption more credible. \citet{Cai_2016_AEJ} uses household-level panel data from $12$ tobacco-producing counties in Jiangxi Province, China, to study how the provision of a tobacco insurance program affects rural households' borrowing and saving behavior. In these counties, some households (hereafter, tobacco households) rely on tobacco as a main source of income, but tobacco production is vulnerable to adverse weather events. In 2003, one county (hereafter, the exposed county) collaborated with the People's Insurance Company of China (PICC) to design and launch a tobacco production insurance program that protected tobacco households against losses from weather-related disasters.
The key feature of this program is that a household is treated if and only if it satisfies two criteria simultaneously: exposure and eligibility. A household is exposed if it resides in the county offering the insurance program, and eligible if it grows tobacco. Households outside the exposed county and non-tobacco households are untreated.
To estimate the effect of the insurance program on tobacco households in the exposed county, a triple-differences design is more appealing than a difference-in-differences design because it accounts for both county- and industry-specific shocks, such as changes in local economic conditions and in the tobacco market, respectively. To see this, consider a difference-in-differences comparison of tobacco households inside and outside the exposed county. This comparison accounts for industry-specific shocks common to tobacco households across counties. However, it may be confounded by county-specific shocks that affect households' financial decisions differently in the exposed and unexposed counties, in which case the parallel trends assumption fails. A triple-differences design adds a second difference-in-differences comparison, between non-tobacco households inside and outside the exposed county, to capture these county-specific shocks. Differencing the two comparisons removes this bias caused by the violation of the parallel trends assumption\footnote{Alternatively, one could start from a difference-in-differences comparison of tobacco and non-tobacco households within the exposed county, which accounts for county-specific shocks. This comparison, however, may be confounded by industry-specific shocks. The triple-differences design removes this bias symmetrically, by differencing out the analogous comparison in the unexposed counties.}. Figure \ref{figure:fourhouseholds} illustrates the two difference-in-differences comparisons underlying this triple-differences design.
Moreover, county- and industry-specific shocks are likely to affect households differently depending on their characteristics. For example, households of different sizes may face different budget constraints and therefore adjust their borrowing and saving differently in response to the same shock. Consequently, deviations from parallel trends may depend on household covariates. In this setting, conditioning on covariates makes the triple-differences identification strategy more credible, because households with similar covariates tend to experience similar deviations from parallel trends.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[
box/.style={
draw,
rounded corners,
minimum width=4.2cm,
minimum height=1.8cm,
align=center,
font=\small
},
treated/.style={
box,
very thick,
fill=gray!15
}
]
\node[font=\small\bfseries] at (-3.5,3.25) {Exposed};
\node[font=\small\bfseries] at ( 3.5,3.25) {Unexposed};
\node[font=\small\bfseries, rotate=90] at (-7.0, 1.5) {Tobacco};
\node[font=\small\bfseries, rotate=90] at (-7.0,-1.5) {Non-tobacco};
\node[treated] (A) at (-3.5,1.5)
{\textcolor{red}{\textbf{Exposed Tobacco}}\\
\textcolor{red}{\textbf{(Treated)}}};
\node[box] (B) at (3.5,1.5)
{Unexposed Tobacco\\ (Untreated)};
\node[box] (C) at (-3.5,-1.5)
{Exposed Non-tobacco\\ (Untreated)};
\node[box] (D) at (3.5,-1.5)
{Unexposed Non-tobacco\\ (Untreated)};
\draw[<->, thick]
(A.east) -- node[above, font=\small]
{$\mathbf{DiD_1}$} (B.west);
\draw[<->, thick]
(C.east) -- node[below, font=\small]
{$\mathbf{DiD_2}$} (D.west);
\node[align=center, font=\small] at (0,0)
{$\mathbf{DDD}=\mathbf{DiD_1}-\mathbf{DiD_2}$};
\end{tikzpicture}
\caption{Four types of households and the two DiD comparisons}
\label{figure:fourhouseholds}
\end{figure}
\newpage
\noindent
\subsection{Econometric Setup and Identification}
Consider a balanced panel that follows $n$ units over $T$ periods, $t=1,\dots,T$. The treatment is implemented simultaneously in a single period $t^*$, with $2\leq t^*\leq T$, so that $t<t^*$ are the pre-treatment periods and $t\geq t^*$ are the post-treatment periods. Each unit is characterized by two time-invariant binary indicators: an exposure indicator \(G\) and an eligibility indicator \(S\). Specifically,
\begin{itemize}
\item the \emph{exposure} indicator \(G \in \{0,1\}\) indicates whether the unit belongs to a group (e.g., a county) that is exposed to the treatment in the post-treatment periods;
\item the \emph{eligibility} indicator \(S \in \{0,1\}\) indicates whether the unit belongs to a subpopulation (e.g., households growing tobacco) that is eligible for the treatment.
\end{itemize}
These two indicators partition the population into four cross-sectional subgroups displayed in Figure \ref{figure:foursubgroup}. The key feature of the triple-differences design is that only units satisfying both criteria are treated: a unit must belong to the exposed group and to the eligible subpopulation. The treatment is an absorbing state: once a unit gets treated, it remains treated through period $T$. Let $D\in\{0,1\}$ indicate whether the unit is treated in the post-treatment periods, and let \(D_{t} \in \{0,1\}\) indicate whether the unit is treated in period $t$, so that
\[D=GS, \quad D_{t}=GS\mathbbm{1}(t\geq t^*) \]
In Figure~\ref{figure:foursubgroup}, the treated subgroup is the top-left cell, \((G=1,S=1)\); we refer to the remaining three cells as the \emph{untreated subgroups}.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[
box/.style={
draw,
rounded corners,
minimum width=4.2cm,
minimum height=1.8cm,
align=center,
font=\small
},
treated/.style={
box,
very thick,
fill=gray!15
}
]
\node[font=\small\bfseries] at (-3.5,3.25) {Exposed};
\node[font=\small] at (-3.5,2.75) {\((G=1)\)};
\node[font=\small\bfseries] at (3.5,3.25) {Unexposed};
\node[font=\small] at (3.5,2.75) {\((G=0)\)};
\node[font=\small\bfseries, rotate=90] at (-7.0,1.5) {Eligible};
\node[font=\small, rotate=90] at (-6.5,1.5) {\((S=1)\)};
\node[font=\small\bfseries, rotate=90] at (-7.0,-1.5) {Ineligible};
\node[font=\small, rotate=90] at (-6.5,-1.5) {\((S=0)\)};
\node[treated] (A) at (-3.5,1.5) {\textbf{Treated subgroup}};
\node[box] (B) at ( 3.5,1.5) {Unexposed Eligible};
\node[box] (C) at (-3.5,-1.5) {Exposed Ineligible};
\node[box] (D) at ( 3.5,-1.5) {Unexposed Ineligible};
\end{tikzpicture}
\caption{Four cross-sectional subgroups}
\label{figure:foursubgroup}
\end{figure}
We adopt the potential-outcome framework in which the potential outcome depends on the entire treatment sequence. Let $Y_{t}(d_{1},\dots,d_{T})$ denote the potential outcome in period $t$ under treatment sequence $(d_{1},\dots,d_{T})$, and let $\mathbf{0}_s$ and $\mathbf{1}_s$ denote $s$-dimensional vectors of zeros and ones, respectively. Because adoption is simultaneous and the treatment is absorbing, only two treatment paths can occur in our setting: treated units follow $(\mathbf{0}_{t^*-1},\mathbf{1}_{T-t^*+1})$ and untreated units follow $\mathbf{0}_{T}$. We therefore abbreviate
\[ Y_{t}(1)\equiv Y_{t}(\mathbf{0}_{t^*-1},\mathbf{1}_{T-t^*+1}), \quad Y_{t}(0)\equiv Y_{t}(\mathbf{0}_{T}) \]
As is common in the causal inference literature on panel data (\citeinb{Chaise-Dhaut_2022_EJ}; \citeinb{Roth-etal_2023_JoE}), we assume SUTVA and consistency to link the observed and potential outcomes. As a result,
\begin{align*}
Y_{t}= Y_{t}(1)D_t+Y_{t}(0)(1-D_t)
\end{align*}
Our parameters of interest are the average treatment effects on the treated subgroup (ATT) $l$ periods after the initial treatment period,
\[
\tau_l=\mathbbm{E}\!\left[Y_{t^*+l}(1)-Y_{t^*+l}(0)\mid D=1\right] = \mathbbm{E}\!\left[Y_{t^*+l}(1)-Y_{t^*+l}(0)\mid G=1,\ S=1\right], \quad l=0,1,\dots,T-t^*
\]
To identify $\tau_l$, we make the following three assumptions:
\begin{assumption}\label{assumption:noanticipation}
(No Anticipation) For all $t=1,\dots,T,$ $Y_{t}(d_1,\dots,d_T)=Y_{t}(d_1,\dots,d_t)$
\end{assumption}
The no-anticipation assumption requires that future treatments do not affect current outcomes: the potential outcome in period $t$ depends only on the treatment path up to period $t$.
Let $X \in \mathbbm{R}^q$ denote a vector of time-invariant, pre-treatment covariates of dimension $q$.
\begin{assumption} \label{assumption:overlap}
(Strong Overlap) There exists \(\epsilon>0\) such that
\[
\Pr(G=g,S=s\mid X=x)>\epsilon
\]
for every \((g,s)\in\{0,1\}^2\) and almost every \(x\in \mathcal{X}\), where \(\mathcal X\) denotes the support of \(X\).
\end{assumption}
The strong overlap assumption requires that, at almost every covariate value $x$, each of the four subgroups occurs with probability bounded away from zero: within each covariate stratum, the population contains units from all four subgroups of Figure~\ref{figure:foursubgroup}.\footnote{Assumption \ref{assumption:overlap} is stronger than what identification of $\tau_l$ requires: it would suffice to impose overlap almost everywhere on the support of $X$ conditional on the treated subgroup, rather than on the whole support $\mathcal{X}$. We maintain the stronger version for simplicity. }
The key identifying assumption of the triple-differences design replaces the familiar parallel trends assumption of the difference-in-differences designs with a parallel gaps assumption, following the recent triple-differences literature (\citeinb{Olden_Moen_2022_EJ}; \citeinb{Ortiz-Villavicencio_SantAnna_2025_working}).
\begin{assumption}\label{assumption:parallelgap}
(Conditional Parallel Gaps) For all $t=2,\dots,T$,
\begin{align*}
&\mathbbm{E}[Y_{t}(0)-Y_{t-1}(0)\mid G=1,\ S=1,\ X] - \mathbbm{E}[Y_{t}(0)-Y_{t-1}(0)\mid G=0,\ S=1,\ X]\\
&= \mathbbm{E}[Y_{t}(0)-Y_{t-1}(0)\mid G=1,\ S=0,\ X] - \mathbbm{E}[Y_{t}(0)-Y_{t-1}(0)\mid G=0,\ S=0,\ X],\ a.s.
\end{align*}
\end{assumption}
Figure~\ref{figure:ddd_parallel_gaps} illustrates the intuition behind the conditional parallel gaps assumption. The assumption allows each of the two difference-in-differences comparisons used in the triple-differences design to violate the conditional parallel trends assumption, and allows the size of this violation to vary with the covariates $X$. It requires only that, at each covariate value $X=x$, the violation be the same in both comparisons. In our motivating example, the conditional parallel gaps assumption allows for both county- and industry-specific shocks, and allows these shocks to affect households' borrowing and saving behavior differently depending on their characteristics. It requires, however, that county-specific shocks generate the same deviation from parallel trends for tobacco and non-tobacco households with similar covariates, or, equivalently, that industry-specific shocks generate the same deviation from parallel trends for households inside and outside the exposed county with similar covariates.
\begin{figure}[htbp]
\centering
\begin{subfigure}[t]{0.48\textwidth}
\centering
\caption*{\large\bfseries Eligible $(S=1)$}
\begin{tikzpicture}[x=0.6cm,y=0.7cm]
\draw[axis] (0,0) -- (8.6,0) node[right] {\scriptsize Time};
\draw[axis] (0,0) -- (0,6.2) node[above] {\scriptsize Outcome};
\draw (1,-0.08) -- (1,0.08) node[below=3pt, font=\scriptsize] {Pre};
\draw (5,-0.08) -- (5,0.08) node[below=3pt, font=\scriptsize] {Post};
\coordinate (C0) at (1,1.2);
\coordinate (C1) at (5,2.4);
\coordinate (T0) at (1,3.0);
\coordinate (TB1) at (5,4.2);
\coordinate (TA1) at (5,5.7);
\draw[control] (C0) -- (C1);
\draw[treated] (T0) -- (TA1);
\draw[benchmark] (T0) -- (TB1);
\fill[blue] (C0) circle (1.5pt);
\fill[blue] (C1) circle (1.5pt);
\fill[red] (T0) circle (1.5pt);
\fill[red] (TA1) circle (1.5pt);
\fill[red!75!black] (TB1) circle (1.3pt);
\draw[bias] ($(TB1)+(0.35,0.05)$) -- ($(TA1)+(0.35,-0.05)$);
\node[font=\scriptsize, inner sep=1pt, anchor=west]
at ($(TB1)!0.5!(TA1)+(0.55,0)$) {$\delta_1(x)$};
\node[anchor=west, font=\scriptsize] at ($(TA1)+(0.55,0)$) {Exposed $(G{=}1,S{=}1)$};
\node[anchor=west, font=\scriptsize] at ($(TB1)+(0.55,0)$) {Trend under CPT};
\node[anchor=west, font=\scriptsize] at ($(C1) +(0.55,0)$) {Unexposed $(G{=}0,S{=}1)$};
\end{tikzpicture}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.48\textwidth}
\centering
\caption*{\large\bfseries Ineligible $(S=0)$}
\begin{tikzpicture}[x=0.6cm,y=0.7cm]
\draw[axis] (0,0) -- (8.6,0) node[right] {\scriptsize Time};
\draw[axis] (0,0) -- (0,6.2) node[above] {\scriptsize Outcome};
\draw (1,-0.08) -- (1,0.08) node[below=3pt, font=\scriptsize] {Pre};
\draw (5,-0.08) -- (5,0.08) node[below=3pt, font=\scriptsize] {Post};
\coordinate (C0) at (1,1.0);
\coordinate (C1) at (5,1.5);
\coordinate (T0) at (1,2.5);
\coordinate (TB1) at (5,3.0);
\coordinate (TA1) at (5,4.5);
\draw[control] (C0) -- (C1);
\draw[treated] (T0) -- (TA1);
\draw[benchmark] (T0) -- (TB1);
\fill[blue] (C0) circle (1.5pt);
\fill[blue] (C1) circle (1.5pt);
\fill[red] (T0) circle (1.5pt);
\fill[red] (TA1) circle (1.5pt);
\fill[red!75!black] (TB1) circle (1.3pt);
\draw[bias] ($(TB1)+(0.35,0.05)$) -- ($(TA1)+(0.35,-0.05)$);
\node[font=\scriptsize, inner sep=1pt, anchor=west]
at ($(TB1)!0.5!(TA1)+(0.55,0)$) {$\delta_0(x)$};
\node[anchor=west, font=\scriptsize] at ($(TA1)+(0.55,0)$) {Exposed $(G{=}1,S{=}0)$};
\node[anchor=west, font=\scriptsize] at ($(TB1)+(0.55,0)$) {Trend under CPT};
\node[anchor=west, font=\scriptsize] at ($(C1) +(0.55,0)$) {Unexposed $(G{=}0,S{=}0)$};
\end{tikzpicture}
\end{subfigure}
\caption{Conditional parallel gaps}
\label{figure:ddd_parallel_gaps}
\vspace{0.3em}
\begin{minipage}{0.95\textwidth}
\footnotesize
\textit{Note:} This figure illustrates the conditional parallel gaps assumption at a fixed covariate value \(X=x\). In each panel, the dotted line is the trend the exposed units would follow if conditional parallel trends (CPT) held. In the left panel, \(\delta_1(x)\) is the deviation from conditional parallel trends in untreated outcomes between eligible exposed units, \((G=1,S=1)\), and eligible unexposed units, \((G=0,S=1)\). In the right panel, \(\delta_0(x)\) is the corresponding deviation between ineligible exposed units, \((G=1,S=0)\), and ineligible unexposed units, \((G=0,S=0)\). The assumption allows conditional parallel trends to be violated in the two DiD comparisons used in the triple-differences design, but requires the violation to be the same across the two DiD comparisons: \(\delta_1(x)=\delta_0(x)\).
\end{minipage}
\end{figure}
Although the conditional parallel gaps assumption allows violations of conditional parallel trends, we emphasize that the former is not a weaker version of the latter. Mathematically, neither condition implies the other. A triple-differences design should therefore not be interpreted as a mechanical robustness check on a DiD design: the two designs rest on distinct identifying assumptions, and their credibility depends on the empirical context. Introducing the second difference-in-differences comparison improves credibility when the second comparison is subject to the same parallel-trends violation as the first, but it can introduce bias conditional parallel trends hold but conditional parallel gaps fail.
\begin{theorem}
\label{theorem:identification}
Under Assumptions \ref{assumption:noanticipation},
\ref{assumption:overlap}, and \ref{assumption:parallelgap}, $\tau_l$ is identified for every $l=0,1,\dots,T-t^*$,
\begin{equation*}
\begin{aligned}
\tau_l
&=
\mathbbm{E}[\Delta Y_l\mid G=1,S=1]
-\mathbbm{E}\!\left[
\frac{f_{11}(X)}{f_{01}(X)}
\Delta Y_l
\mid G=0,S=1
\right] \\
&\quad
-\mathbbm{E}\!\left[
\frac{f_{11}(X)}{f_{10}(X)}
\Delta Y_l
\mid G=1,S=0
\right]
+\mathbbm{E}\!\left[
\frac{f_{11}(X)}{f_{00}(X)}
\Delta Y_l
\mid G=0,S=0
\right] ,
\end{aligned}
\end{equation*}
where $\Delta Y_l\equiv Y_{t^*+l}-Y_{t^*-1}$ and \(f_{gs}(x)\) denotes the conditional density of \(X\) given
$(G,S)=(g,s)$.
\end{theorem}
Theorem \ref{theorem:identification} gives a weighting representation of $\tau_l$. The treated subgroup is left unweighted, and each untreated subgroup $(g,s)$ is reweighted by the density ratio $f_{11}(X)/f_{gs}(X)$. The density ratio is a balancing weight: it aligns the covariate distribution of the untreated subgroup $(g,s)$ with that of the treated subgroup. The identification result thus pins down a balancing weight profile:
\begin{equation}
\label{eq:balancing-weights}
w_{11}(X)=1,
\qquad
w_{gs}(X)=\frac{f_{11}(X)}{f_{gs}(X)}
\quad\text{for }(g,s)\in\{(0,1),(1,0),(0,0)\}.
\end{equation}
This profile yields a sufficient condition for a weighting estimand to equal the ATT, which we call the
\emph{covariate balancing condition}: the estimand’s weighting profile (i)balances the covariate distribution of each untreated subgroup with that ofthe treated subgroup, and (ii) attains this balance while leaving the treated subgroup unweighted. Consequently, an estimator that converges to such a weighting estimand is consistent for the ATT. This condition plays a key role in both our diagnosis of empirical practice in Section \ref{section:empiricaldiagnose} and our construction of the matching estimator in Section \ref{section:matchedtriplediffereces}.
\section{Diagnosis of Empirical Practice}\label{section:empiricaldiagnose}
\noindent
In this section, we diagnose the two leading empirical practices in the applied literature. Consider a sample $\{Y_{i1},Y_{i2},\dots, Y_{iT}, G_i,S_i,X_i\}_{i=1}^n$, drawn independently from the same distribution.
\begin{assumption}\label{assumption:iid}
(Random Sampling) $\{Y_{i1},Y_{i2},\dots, Y_{iT}, G_i,S_i,X_i\}_{i=1}^n$ is a sample of independent and identically distributed observations from a population $(Y_{1},Y_{2},\dots,Y_{T},G,S,X)$.
\end{assumption}
As documented in our survey of $18$ empirical papers (Table \ref{table:ddd_literature}), empirical practice in the literature largely follows two types of 3WFE regression specifications. In the first specification, which we call the covariate-trend 3WFE regression, the covariates $X_i$ enter with covariate specific trends, allowing for time varying coefficients:
\begin{equation}
\label{eq:spec1}
Y_{it}
= \lambda_{G_iS_i}
+ \lambda_{G_i t}
+ \lambda_{S_i t}
+ \sum_{l\neq -1}
\tau_l^{CT}\, G_iS_i\mathbbm{1}(t-t^*=l)
+ \sum_{r=1}^T X_i'\theta_r\mathbbm{1}(t=r)
+ u_{it}.
\end{equation}
In the second specification, which we call the covariate-level 3WFE regression, the covariates $X_i$ enter only in levels, with time-invariant coefficients:
\begin{equation}
\label{eq:spec2}
Y_{it}
= \lambda_{G_iS_i}
+ \lambda_{G_i t}
+ \lambda_{S_i t}
+ \sum_{l\neq -1}
\tau_l^{CL}\, G_iS_i\mathbbm{1}(t-t^*=l)
+ X_i'\theta
+ u_{it}.
\end{equation}
In both specifications, $\lambda_{G_iS_i}$, $\lambda_{G_i t}$, and $\lambda_{S_i t}$ are the subgroup, exposure-by-time, and eligibility-by-time fixed effects, respectively. The OLS coefficients $\widehat \tau_l^{CT}$ and $\widehat \tau_l^{CL}$ on the triple interaction term $ G_iS_i\mathbbm{1}(t-t^*=l)$ are often interpreted as estimates of $\tau_l$, the ATT $l$ periods after the treatment. In a balanced panel, replacing the subgroup fixed effect $\lambda_{G_iS_i}$ with an individual fixed effect $\lambda_i$ yields identical OLS estimates in both the covariate-trend and covariate-level specifications. Because more than two-thirds the papers in our survey use subgroup fixed effects, we focus our diagnosis on the specifications with subgroup fixed effects; since the estimates are identical, the diagnosis below also applies to the specifications with individual fixed effects.
For each specification, our diagnosis consists of two parts. First, we rewrite the corresponding OLS estimand as a weighting estimand and determine whether its weighting profile satisfies the covariate balancing condition. Second, if it does not, we characterize the additional assumptions under which the OLS estimand equals the ATT. Our results show that common empirical practice generally fails to consistently estimate the ATT: the induced weighting estimands do not correctly align the covariate distributions of the untreated subgroups with that of the treated subgroup, and they require additional restrictive assumptions to equal the ATT.
\subsection{Covariate-trend 3WFE Regression}
Recall that $D=GS$ is the treatment indicator. Let $L(D\mid G,S,X)$ denote the population linear projection of $D$ onto $(1,G,S,X)$, let $\kappa\equiv\mathbbm{E}[(D-L(D\mid G,S,X))^2]$ denote the mean squared projection residual, and let $p_{gs}=P(G=g,S=s)$ denote the subgroup proportions. Let $\tau_l(X)=\mathbbm{E}\!\left[Y_{t^*+l}(1)-Y_{t^*+l}(0)\mid G=1,\ S=1,X\right] $ denote the conditional ATT, and let $ L_0(\Delta Y_l|G,S,X)$ denote the conditional linear projection of $\Delta Y_l$ on $(1,G,S,X)$ in the untreated population $D=0$.
\begin{lemma}\label{lemma:covariatetrend}
\textbf{(I)} Under Assumptions~\ref{assumption:overlap} and \ref{assumption:iid},
\begin{align*}
\widehat \tau_l^{CT} &\xrightarrow{p}\tau_{l}^{CT}\\
&=\mathbbm{E}\!\left[\frac{p_{11}\left(1-L(D\mid G,S,X)\right)}{\kappa}\, \Delta Y_l \;\middle|\; G=1,S=1\right]
-\mathbbm{E}\!\left[\frac{p_{01}\,L(D\mid G,S,X)}{\kappa}\, \Delta Y_l \;\middle|\; G=0,S=1\right] \\
&\quad-\mathbbm{E}\!\left[\frac{p_{10}\,L(D\mid G,S,X)}{\kappa}\, \Delta Y_l \;\middle|\; G=1,S=0\right]
+\mathbbm{E}\!\left[\frac{-p_{00}\,L(D\mid G,S,X)}{\kappa}\, \Delta Y_l \;\middle|\; G=0,S=0\right].
\end{align*}
\textbf{(II)} $\tau_l^{CT}=\tau_l$ if the treatment effect is \emph{homogeneous} in $X$, i.e., $\tau_l(X)=\tau_l$ and
\emph{either} of the following holds:
\begin{enumerate}
\item Linear conditional means with a common slope.
\[
\mathbbm{E}[\Delta Y_l \mid G=g,S=s,X]
= L_0(\Delta Y_l|g,s,X) = \gamma_{gs} + X'\gamma,
\qquad (g,s)\in\{(0,0),(0,1),(1,0)\}.
\]
\item An inverse linear functional form for generalized propensity scores. Define $e_{gs}(X)=P(G=g,S=s|X)$.
\begin{align*}
e_{gs}(X)&=\frac{[(-1)^{g+s}L(D|g,s,X)]^{-1}}{\sum_{g=0}^1\sum_{s=0}^1[gs-(-1)^{g+s}L(D|g,s,X)]^{-1}}\\
&=
\frac{
[ \lambda_{gs} +(-1)^{g+s}X'\lambda ]^{-1}
}{
\displaystyle
\sum_{g=0}^{1}\sum_{s=0}^{1}
[ \lambda_{gs} +(-1)^{g+s}X'\lambda ]^{-1}
}.
\end{align*}
\end{enumerate}
\end{lemma}
Part (I) of Lemma \ref{lemma:covariatetrend} shows that the OLS estimator from the covariate-trend 3WFE regression converges to a weighting estimand, in which each of the four subgroups is reweighted by weights proportional to the linear projection residual $D-L(D\mid G,S,X)$. This weighting profile is the triple-differences analogue of the implied regression weights studied by \citet{Chatt_Zubi_2023_biometrika} in cross-sectional settings and by \citet{Caetano-Callaway_2026_wp} in DiD settings. In those two settings, the implied regression weights have an appealing property: they balance the covariate means of the treated and untreated groups. In the triple-differences setting, however, the implied regression weights induced by specification \eqref{eq:spec1} fail the covariate balancing condition in three respects.
First, they target the wrong balance structure (pooling issue). Lemma~\ref{lemma:regressionweight} in the appendix shows that the implied regression weights balance covariate means between the treated subgroup ($D=1$) and the pooled untreated group ($D=0$). In a triple-differences design this is the wrong balance structure because the untreated group ($D=0$) now pools three distinct untreated subgroups, whereas the covariate balancing condition requires aligning the covariate distribution of \emph{each} of them with the treated subgroup separately. Second, the implied regression weights attain only mean balance rather than the balance of the full covariate distribution required by the covariate balancing conditions (mean only balance issue): Two subgroups can still have different covariate distributions even when their covariate means are balanced. Third, the implied regression weights also reweight the treated subgroup (treated-reweighting issue). The covariate balancing condition leaves the treated subgroup unweighted while the implied regression weights profile reweights the treated subgroup by the weight $p_{11}(1-L(D\mid G,S,X))/\kappa$.
Part (II) of lemma \ref{lemma:covariatetrend} provides the additional sufficient assumptions for $\tau_l^{CT}$ to equal the ATT. The first is that the treatment effect is homogeneous with respect to the covariates. The second is that either the conditional mean functions of the three untreated subgroups are linear in $X$ with a common slope, differing only in their intercepts, or the generalized propensity scores admit the specific inverse linear functional form. Both requirements are restrictive: the first rules out treatment effect heterogeneity in $X$, and the second imposes a parametric form on either the conditional mean functions or the generalized propensity scores.
\begin{remark}
Specification \eqref{eq:spec1} can be enriched in two ways: (i) by further interacting the covariates with the subgroup indicators $G_i$ and $S_i$, and (ii) by including higher-order terms of the covariates. The first enrichment can attain mean balance between each untreated subgroup and the treated subgroup, addressing the pooling issue. The second enrichment can attain balance of higher-order moments of the covariates, partly addressing the mean-only balance issue. However, the enriched specifications still do not satisfy the covariate balancing condition: first, balance of higher-order moments among the four subgroups does not imply balance of the full covariate distribution, and second, the implied regression weights continue to reweight the treated subgroup. Hence, the enriched specifications do not consistently estimate $\tau_l$ in general.
\end{remark}
\subsection{Covariate-level 3WFE Regression}
\begin{lemma}\label{lemma:covariatelevel}
\textbf{(I)} Under Assumptions~\ref{assumption:overlap} and \ref{assumption:iid},
\begin{align*}
\widehat\tau_l^{CL}&\xrightarrow{p}\tau_l^{CL}\\
&=\mathbbm{E}[\Delta Y_l \mid G=1,S=1] -\mathbbm{E}[\Delta Y_l \mid G=0,S=1]- \mathbbm{E}[\Delta Y_l \mid G=1,S=0] + \mathbbm{E}[\Delta Y_l \mid G=0,S=0]
\end{align*}
\textbf{(II)} $\tau_l^{CL}=\tau_l$ \emph{if and only if} the parallel gaps assumption holds \textbf{unconditionally}:
\begin{align*}
&\mathbbm{E}[\Delta Y_l(0)\mid G=1,\ S=1] - \mathbbm{E}[\Delta Y_l(0)\mid G=0,\ S=1]\\
&= \mathbbm{E}[\Delta Y_l(0)\mid G=1,\ S=0] - \mathbbm{E}[\Delta Y_l(0)\mid G=0,\ S=0]
\end{align*}
\end{lemma}
Lemma \ref{lemma:covariatelevel} shows that the OLS estimator from the covariate-level 3WFE regression specification converges to the unconditional triple-differences estimand, $\tau_l^{CL}$. This estimand uses the constant weight profile $\widetilde w_{gs}(X)=1$ for all four subgroups and fails the sufficient covariate balancing condition because none of the three untreated subgroups is reweighted. The necessary and sufficient condition for $\tau_l^{CL} $ to equal the ATT is that the parallel gaps assumption holds unconditionally, in which case there is no need to include covariates in the triple-differences design.
In addition, although this 3WFE regression still produces OLS coefficient on $X$, in a balanced panel $\widehat{\tau}_l^{CL}$ is numerically identical to the estimator from the 3WFE regression without covariates.\footnote{In the alternative specification with individual fixed effects, discussed at the beginning of this section, the time-invariant covariates are absorbed by the fixed effects.} Therefore, including covariates in this way neither helps the identification nor improves the efficiency of estimation.
\section{Matched Triple-Differences: Estimation and Inference}\label{section:matchedtriplediffereces}
The diagnosis in Section \ref{section:empiricaldiagnose} shows that the OLS estimators from the two commonly used 3WFE specifications are, in general, inconsistent for the ATT. In this section, we propose a matched triple-differences framework for valid estimation and inference. The core of our framework consists of two parts: (i) an estimation procedure for constructing matching estimators in the triple-differences design; and (ii) a general asymptotic theory and inference that cover a class of matching procedures, such as the commonly used nearest-neighbor matching with a fixed or diverging number of neighbors and kernel matching.
\subsection{Estimation}
Estimation in our framework follows a two-step procedure. Let $\mathcal I_{gs}\equiv\left\{i: G_i=g,\, S_i=s\right\}$, and $n_{gs}=\left\lvert\mathcal I_{gs}\right\rvert$, $(g,s)\in\{0,1\}^2$, denote the four subgroups and their sizes. In the first step, we conduct three pairwise covariate matchings; in each, the treated subgroup $\mathcal{I}_{11}$ is matched with one of the three untreated subgroups $\mathcal{I}_{gs}$, $(g,s)\neq(1,1)$. The goal of this step is to obtain a matching weight for every unit in each of the three untreated subgroups. In Figure \ref{figure:matching_workflow_proposed}, we visualize the three pairwise matchings and the resulting matching weights.
To be specific, consider the pairwise covariate matching between the treated subgroup $\mathcal{I}_{11}$ and the untreated subgroup $\mathcal{I}_{gs}$. Each treated unit $i$ is matched to its \emph{matched set} $\mathcal{M}(i)$ which consists of units in the untreated subgroup $\mathcal{I}_{gs}$ with covariates similar to $X_i$. Each untreated unit $j\in \mathcal{M}(i)$ receives a non-negative \emph{pairwise weight} $w_{ij}(\mathcal I_{11},\mathcal I_{gs})$, which measures the contribution of untreated unit $j$ to the matched set of treated unit $i$. The pairwise weights are computed by a prespecified matching algorithm (e.g., nearest-neighbor or kernel matching). An untreated unit $j$ may appear in the matched sets of several treated units and can therefore receive several pairwise weights. Summing all the pairwise weights that unit $j$ receives yields the \emph{matching weight} $w_{j}(\mathcal I_{11},\mathcal I_{gs})$, which measures how intensely untreated unit $j$ is used in the pairwise matching:
\[
w_{j}(\mathcal I_{11},\mathcal I_{gs}) \equiv \sum_{i\in\mathcal I_{11}} w_{ij}(\mathcal I_{11},\mathcal I_{gs}),
\qquad j\in\mathcal I_{gs},
\]
In Figure \ref{figure:pairwisematching}, we illustrate, through a simplified example, how the matching weights are obtained in the pairwise matching. Carrying out all three pairwise matchings gives each untreated unit $j\in \mathcal{I}_{gs}$ a matching weight $w_{j}(\mathcal I_{11},\mathcal I_{gs})$ and leads to the following weighting profile:
\begin{align*}
w_i=
\begin{cases}
1 & i \in \mathcal I_{11}, \\[2pt]
w_i\!\left(\mathcal I_{11},\,\mathcal I_{gs}\right) & i \in \mathcal I_{gs},\ (g,s)\neq(1,1) \\
\end{cases}
\end{align*}
In the second step, we run a weighted 3WFE regression using the weighting profile we obtained in the first step:
\begin{equation}
\label{eq:ddd-regression}
\sqrt{w_i}Y_{it}
=
\sqrt{w_i}\left(\lambda_{G_iS_i}
+\lambda_{G_i t}
+\lambda_{S_i t}
+\sum_{l\neq-1}
\tau_l^{match} G_iS_i\mathbbm{1}(t-t^*=l)
+u_{it}\right).
\end{equation}
The coefficient on the triple interaction term, $\widehat{\tau}_l^{match}$, is our matched triple-differences estimator. It can then be equivalently written as
\begin{align*}
\widehat \tau_l^{match}&= \frac{1}{n_{11}}
\sum_{i\in\mathcal{I}_{11}}
\Delta Y_{il}
-\frac{1}{n_{11}}
\sum_{i\in\mathcal{I}_{01}}w_i\!\left(\mathcal I_{11},\,\mathcal I_{01}\right)\Delta Y_{il}\\
&-\left(\frac{1}{n_{11}}
\sum_{i\in\mathcal{I}_{10}}w_i\!\left(\mathcal I_{11},\,\mathcal I_{10}\right)\Delta Y_{il}
-\frac{1}{n_{11}}
\sum_{i\in\mathcal{I}_{00}}w_i\!\left(\mathcal I_{11},\,\mathcal I_{00}\right)\Delta Y_{il}\right).
\end{align*}
Our matched triple-differences estimator is closely connected to the covariate balancing condition in Section \ref{section:econometricssetup}. To see this connection, note that in $\widehat{\tau}_l^{match}$, the treated subgroup is unweighted, while the matching weights $w_i\!\left(\mathcal I_{11},\,\mathcal I_{gs}\right)$ approximate the density ratio, so that each untreated subgroup is reweighted to have the approximately the same covariate distribution as the treated subgroup.
\begin{figure}[htbp]
\centering
\savebox{\figbox}{
\begin{tikzpicture}[
>={Latex[length=2.4mm,width=2mm]},
node distance=1.6cm and 3.4cm,
box/.style={
draw, line width=0.6pt, rounded corners=3pt,
minimum width=4.8cm, minimum height=1.7cm,
align=center, font=\small
},
treated/.style={box, line width=1.2pt, fill=gray!20},
header/.style={font=\small\bfseries, align=center},
match/.style={->, line width=0.7pt},
lab/.style={font=\footnotesize, inner sep=2.5pt}
]
\node[treated] (A) {Treated subgroup $(\mathcal{I}_{11})$\\[2pt] $w_i = 1$};
\node[box, right=of A] (B) {Unexposed Eligible $(\mathcal{I}_{01})$\\[2pt] $w_i = w_i(\mathcal{I}_{11},\mathcal{I}_{01})$};
\node[box, below=of A] (C) {Exposed Ineligible $(\mathcal{I}_{10})$\\[2pt] $w_i = w_i(\mathcal{I}_{11},\mathcal{I}_{10})$};
\node[box, below=of B] (D) {Unexposed Ineligible $(\mathcal{I}_{00})$\\[2pt] $w_i = w_i(\mathcal{I}_{11},\mathcal{I}_{00})$};
\node[header, above=0.45cm of A] {Exposed\\ $(G_i = 1)$};
\node[header, above=0.45cm of B] {Unexposed\\ $(G_i = 0)$};
\node[header, rotate=90, anchor=south] at ([xshift=-0.45cm]A.west) {Eligible\\ $(S_i = 1)$};
\node[header, rotate=90, anchor=south] at ([xshift=-0.45cm]C.west) {Ineligible\\ $(S_i = 0)$};
\draw[match] (A) -- node[lab, above] {pairwise matching 1} (B);
\draw[match] (A) -- node[lab, left] {pairwise matching 2} (C);
\draw[match] (A) -- node[lab, above, sloped] {pairwise matching 3} (D);
\end{tikzpicture}}
\begin{adjustbox}{max width=\linewidth}\usebox{\figbox}\end{adjustbox}
\caption{The matched triple-differences procedure}
\label{figure:matching_workflow_proposed}
\vspace{0.35em}
\setlength{\notewidth}{\wd\figbox}
\ifdim\notewidth>\linewidth \setlength{\notewidth}{\linewidth}\fi
\begin{minipage}{\notewidth}
\footnotesize
\raggedright
\emph{Note.} The matched triple-differences procedure performs three separate pairwise matchings. In each pairwise matching, every unit in the treated subgroup \((\mathcal{I}_{11})\) is matched to similar units in one untreated subgroup, based on observed covariates $X_i$.
\end{minipage}
\end{figure}
\begin{figure}[htbp]
\centering
\sbox{\figbox}{
\resizebox{!}{0.25\textheight}{
\begin{tikzpicture}[
x=1cm,y=1cm,
every node/.style={font=\footnotesize},
unit/.style={
circle,draw=black!60,fill=white,
minimum size=6mm,inner sep=0pt,
font=\normalsize
},
pair/.style={
-{Latex[length=1.8mm,width=1.3mm]},thick
},
wt/.style={
fill=white,inner sep=1.5pt,
font=\footnotesize
},
eq/.style={
anchor=west,inner sep=0pt,
font=\footnotesize
}
]
\fill[myblue!7]
(5,.625) ellipse (.85cm and 1.05cm);
\fill[red!7,opacity=.60]
(5,-.625) ellipse (.85cm and 1.05cm);
\draw[myblue,thick]
(5,.625) ellipse (.85cm and 1.05cm);
\draw[red!80!black,thick,dashed]
(5,-.625) ellipse (.85cm and 1.05cm);
\node[font=\footnotesize\bfseries,align=center,anchor=north]
at (0,2.95) {Treated\\$\mathcal I_{11}$};
\node[font=\footnotesize\bfseries,align=center,anchor=north]
at (5,2.95) {Untreated\\$\mathcal I_{gs}$};
\node[unit,draw=myblue,thick]
(i1) at (0,1.25) {$i_1$};
\node[unit,draw=red!80!black,thick]
(i2) at (0,-1.25) {$i_2$};
\node[unit] (j1) at (5,1.25) {$j_1$};
\node[unit,draw=black!80,thick]
(j2) at (5,0) {$j_2$};
\node[unit] (j3) at (5,-1.25) {$j_3$};
\draw[pair,myblue] (i1) -- (j1);
\draw[pair,red!80!black,dashed] (i2) -- (j3);
\draw[pair,myblue] (i1) --
node[wt,text=myblue,pos=.5]
{$w_{i_1j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})$} (j2);
\draw[pair,red!80!black,dashed] (i2) --
node[wt,text=red!80!black,pos=.5]
{$w_{i_2j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})$} (j2);
\node[eq] at (6.3,0) {
$\begin{aligned}
w_{j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})
&=\textcolor{myblue}{w_{i_1j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})}\\[-1pt]
&\quad+\textcolor{red!80!black}{w_{i_2j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})}
\end{aligned}$
};
\coordinate (leg) at ([yshift=-.5cm]current bounding box.south);
\node[anchor=east,inner sep=0pt] (leg1)
at ([xshift=-.5cm]leg) {Matched set for $i_1$};
\draw[myblue,thick]
([xshift=-.65cm]leg1.west) -- ([xshift=-.15cm]leg1.west);
\draw[red!80!black,thick,dashed]
([xshift=.5cm]leg) -- ([xshift=1cm]leg);
\node[anchor=west,inner sep=0pt]
at ([xshift=1.15cm]leg) {Matched set for $i_2$};
\end{tikzpicture}
}
}
\begin{adjustbox}{max width=\linewidth}\usebox{\figbox}\end{adjustbox}
\caption{Construction of the matching weights}
\label{figure:pairwisematching}
\vspace{0.35em}
\setlength{\notewidth}{\wd\figbox}
\ifdim\notewidth>\linewidth \setlength{\notewidth}{\linewidth}\fi
\begin{minipage}{\notewidth}
\footnotesize
\raggedright
\emph{Note.} This figure illustrates how to obtain the matching weight $w_{j_2}(\mathcal{I}_{11},\mathcal{I}_{gs})$ for untreated unit $j_2\in\mathcal{I}_{gs}$ in the pairwise matching between the treated subgroup $\mathcal{I}_{11}$ and untreated subgroup $\mathcal{I}_{gs}$.
\end{minipage}
\end{figure}
\subsection{Asymptotic Theory}
The key novelty of our asymptotic theory is a set of regularity conditions on the pairwise covariate matching used in estimation. These conditions allow our framework to abstract from any specific matching method and to provide a general asymptotic theory for a class of estimators constructed from different matching methods. To avoid clutter, we suppress the dependence of the weights on the pair of subgroups, writing $w_{ij}$ and $w_j$ and we write $\widehat{\tau}_l$ for $\widehat{\tau}_l^{match}$ when no confusion arises.
\begin{assumption}[Regularity of the pairwise covariate matching]
\label{assumption:matching}
For each $(g,s)\in\{(0,1),(1,0),(0,0)\}$, the pairwise covariate matching between the treated subgroup $\mathcal I_{11}$ and untreated subgroup $\mathcal I_{gs}$ satisfies the following conditions.
\begin{enumerate}[label=(\roman*)]
\item The pairwise weights $\left\{w_{ij}:i\in\mathcal I_{11},\,j\in\mathcal I_{gs} \right\}$ are constructed using only the subgroup indicators and observed covariates of the treated and untreated subgroups, $(G_i,S_i,X_i)_{i\in \mathcal{I}_{11}\cup \mathcal{I}_{gs}}$. They do not depend on the order in which the observations are supplied: permuting the observations within $\mathcal I_{11}$ and $\mathcal I_{gs}$ permutes the corresponding pairwise weights accordingly. Moreover, for each $i\in\mathcal I_{11}$,
\[
w_{ij}\geq 0 \quad\text{for every }j\in\mathcal I_{gs},
\qquad
\sum_{j\in\mathcal I_{gs}} w_{ij} = 1.
\]
\item \textbf{(Locality)}
\[
\mathbbm{E}\!\left[
\sum_{j\in\mathcal I_{gs}} w_{ij}\left\lVert X_j-X_i\right\rVert
\,\middle|\, i\in\mathcal I_{11}
\right] = o(1).
\]
\item \textbf{(Bounded Moments)} For some $\delta>0$,
\[
\mathbbm{E}\!\left[ w_{j}^{4+\delta} \,\middle|\, j\in\mathcal I_{gs} \right] =O(1).
\]
\item \textbf{(Stability)} Let $w_{j}^{(1)}$ denote the matching weight of unit $j$ computed after replacing the first observation $(Y_{11},\dots, Y_{1T}, G_1,S_1,X_1 ) $ with an independent copy $(Y_{1'1},\dots, Y_{1'T}, G_{1'},S_{1'},X_{1'} )$.
\[
\mathbbm{E}\!\left[\Bigg(\sum_{j\in\mathcal I_{gs}\setminus\{1\}}
\bigl|\,w_{j}-w_{j}^{(1)}\bigr|\Bigg)^{\!4}\right]
= O(1).
\]
\end{enumerate}
\end{assumption}
Condition (i) restricts the information available to the matching procedure: the pairwise weights depend only on the subgroup indicators and observed covariates, and they are invariant to the order in which observations are supplied. The non-negativity and adding-up requirements imply that each treated unit is matched to a convex combination of untreated units.
Condition (ii) imposes locality: the covariate discrepancy between a treated unit and its matched untreated units, averaged using the pairwise weights, vanishes as the sample grows. In other words, each treated unit is asymptotically matched to untreated units with nearly identical covariate values.
Condition (iii) restricts how intensively an untreated unit may be reused across matched sets. The $4+\delta$ moment bound rules out matching procedures that heavily reuse only a handful of untreated units.
Condition (iv) is a stability requirement on the matching procedure: it rules out procedures in which a single observation can exert a large influence on the overall allocation of the matching weights. This condition cannot simply be removed if we want to establish unconditional asymptotic normality, because under conditions (i)--(iii) alone, asymptotic normality of the matching estimators is guaranteed only conditionally on the subgroup indicators and the covariates. In Appendix \ref{section:example}, we provide an example of a matching procedure that satisfies conditions (i)--(iii) but not (iv), for which the resulting matched triple-differences estimator fails to be asymptotically normal unconditionally.
Assumption~\ref{assumption:matching} is satisfied by several commonly used matching procedures: Lemma~\ref{lemma:matchingprocedure} in the Appendix \ref{appendix:matchingprocedure} shows that it holds for nearest-neighbor matching with replacement, with either a fixed or a diverging number of neighbors and for kernel matching. The assumption does, however, exclude some procedures used in the empirical literature. Nearest-neighbor matching \emph{without} replacement violates condition (i), because its pairwise weights depend on the order in which the observations are supplied; local linear matching violates condition (i) as well, because its pairwise weights can be negative.
We also impose the following regularity conditions, which are standard in the matching literature (e.g., \citeinb{Abadie-Imbens_2006_ECMA}).
\begin{assumption}[Regularity conditions]\label{assumption:regularity}
For every $(g,s)\in\{0,1\}^2$, $d\in\{0,1\}$, $t=1,\dots,T$, and event time $l=0,1,\dots,T-t^*$:
\begin{enumerate}[label=(\roman*)]
\item the conditional density $f_{gs}$ of $X$ given $(G,S)=(g,s)$ is continuous, has convex and compact support, and is bounded and bounded away from zero on its support;
\item $\mathbbm{E}[Y_{t}(d)\mid G=g,S=s,X=x]$ is Lipschitz continuous in $x$;
\item $\mathrm{Var}[\Delta Y_{l}(d)\mid G=g,S=s,X=x]$ is Lipschitz continuous in $x$ and uniformly bounded above and away from zero, that is, $0<\underline{\sigma}^2\leq \mathrm{Var}[\Delta Y_{l}(d)\mid G=g,S=s,X=x]\leq \overline{\sigma}^2$ ;
\item $\mathbbm{E}[Y_{t}(d)^4\mid G=g,S=s,X=x]$ is uniformly bounded.
\end{enumerate}
\end{assumption}
For notational simplicity, define $\mu_l^{gs}(x)=\mathbbm{E}[\Delta Y_l|G=g,S=s,X=x]$ and $\sigma^2_{gs,l}(x)=Var(\Delta Y_l|G=g,S=s,X=x)$.
\begin{theorem}\label{theorem:asymptotic}
Under Assumptions \ref{assumption:noanticipation},
\ref{assumption:overlap}, \ref{assumption:parallelgap}, \ref{assumption:iid}, \ref{assumption:matching}, and \ref{assumption:regularity},
\begin{align*}
\widehat\tau_l-\tau_l\xrightarrow{p} 0, \quad
(V_{\tau(X)}+V_E)^{-1/2} \sqrt{n_{11}}\left(\widehat{\tau}_l-\tau_l-B \right)\xrightarrow{d}\mathcal{N}(0,1)
\end{align*}
where
\begin{align*}
V_{\tau(X)}& =\mathbbm{E}\left[\left. \left(\mu_l^{11}(X_i)-\mu_l^{01}(X_i)-\mu_l^{10}(X_i)+\mu_l^{00}(X_i)-\tau_l \right)^2\right\vert G_i=1,S_i=1\right] \\
V_E&= \frac{1}{n_{11}}\left(\sum_{i\in\mathcal{I}_{11}}\sigma^2_{11,l}(X_{i})+\sum_{i\in\mathcal{I}_{01}} w_{i}^2\sigma^2_{01,l}(X_{i})+\sum_{i\in\mathcal{I}_{10}} w_{i}^2\sigma^2_{10,l}(X_{i}) +\sum_{i\in\mathcal{I}_{00}} w_{i}^2\sigma^2_{00,l}(X_{i})\right)\\
B&=\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{01})w_i\right]\mu_l^{01}(X_i)+\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{10})w_i\right]\mu_l^{10}(X_i)\\
&-\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{00})w_i\right]\mu_l^{00}(X_i)
\end{align*}
\end{theorem}
Theorem \ref{theorem:asymptotic} establishes consistency and asymptotic normality for the class of matched triple-differences estimators that fall within our framework. It extends classical asymptotic results for matching estimators, such as those of \cite{Abadie-Imbens_2006_ECMA}, both to the triple-differences setting and to a class of matching procedures.
Theorem \ref{theorem:asymptotic} also shows that our matched triple-differences estimator carries a finite-sample bias term $B$. The bias term does not affect consistency, but it can affect inference: $\sqrt{n_{11}}B$ need not be asymptotically negligible, so the matched triple-differences estimator may require bias correction for valid inference. This issue parallels the one first pointed out by \cite{Abadie-Imbens_2006_ECMA} for cross-sectional matching estimators. In some cases, bias correction is not required. For example, if covariates are discrete, we can match on them exactly, and matching bias is zero. In addition, the results of \cite{Viel_2025_wp} imply that, for nearest-neighbor matching and kernel matching with properly chosen tuning parameters, the matching bias is root-$n$ negligible when the dimension of the continuous covariate is at most $3$, so bias correction is not required\footnote{Let $r$ denote the number of continuous covariates and $(g,s)\in\{(0,1),(1,0),(0,0)\}$. The matching bias is root-$n$ negligible if (i) for nearest-neighbor matching, $M=o(n_{gs}^{2/3})$ when $r=1$, $M=o(n_{gs}^{1/2})$ when $r=2$, and $M=o(n_{gs}^{1/4})$ when $r=3$; and (ii) for kernel matching, $h=n_{gs}^{-a}$ with $1/4<a<1$ when $r=1$, $1/4<a<1/2$ when $r=2$, and $1/4<a<1/3$ when $r=3$.}.
When bias correction is required, we can correct our matched triple-differences estimator by adapting the bias-correction procedure in the cross-sectional setting \citep{Abadie-Imbens_2011_JBES,Lin_Han_2025_JOE} to each of the three terms of the bias $B$ in Theorem \ref{theorem:asymptotic}. To conduct the bias correction, we estimate the three unknown conditional mean functions $\mu_l^{gs}$, $(g,s)\in\{(0,1),(1,0),(0,0)\}$, nonparametrically by $\widehat{\mu}_l^{gs}$, using off-the-shelf methods such as series or spline estimators. The bias-corrected matched triple-differences estimator then takes the form:
\begin{align*}
\widehat{\tau}_l^{bc}&= \widehat{\tau}_l -\widehat{B} \\
&=\widehat{\tau}_l-\left(\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{01})w_i\right]\widehat\mu_l^{01}(X_i)+\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{10})w_i\right]\widehat\mu_l^{10}(X_i)\right.\\
&\left.-\frac{1}{n_{11}} \sum_{i=1}^{n} \left[\mathbbm{1}(i\in \mathcal{I}_{11})-\mathbbm{1}(i\in \mathcal{I}_{00})w_i\right]\widehat\mu_l^{00}(X_i)\right)
\end{align*}
\begin{theorem}\label{theorem:biascorrection}
Under Assumptions \ref{assumption:noanticipation},
\ref{assumption:overlap}, \ref{assumption:parallelgap}, \ref{assumption:iid}, \ref{assumption:matching}, \ref{assumption:regularity}, and \ref{assumption:biascorrection},
\begin{align*}
\sqrt{n_{11}}(\widehat{B}-B)\xrightarrow{p} 0,\quad
(V_{\tau(X)}+V_E)^{-1/2} \sqrt{n_{11}}\left(\widehat{\tau}_l-\tau_l-\widehat B \right)\xrightarrow{d}\mathcal{N}(0,1).
\end{align*}
\end{theorem}
Theorem \ref{theorem:biascorrection} shows that the bias can be estimated with an error of order smaller than $n_{11}^{-1/2}$. Consequently, the bias-corrected matched triple-differences estimator is $\sqrt{n_{11}}$-consistent and asymptotically normal, and the correction leaves the asymptotic variance unchanged, so the variance estimator in Section \ref{subsection:varianceestimation} still applies to our bias-corrected estimator.
\begin{remark}
Although our bias-corrected estimator is motivated by removing the matching bias, equation (\ref{eq:AIPW}) shows that it has the structure of an augmented inverse probability weighting (AIPW) estimator. The first term is the outcome regression part and the remaining three terms form the weighting (IPW) part, with the inverse propensity score weights replaced by the matching weights. This parallels the novel insight of \cite{Lin-Ding-Han_2023_ECMA} and \cite{Lin_Han_2025_JOE}, who point out that the bias-corrected matching estimator in the cross-sectional setting has the structure of an AIPW estimator and provide conditions under which it enjoys the desirable properties of AIPW estimators such as double robustness and attainment of the semiparametric efficiency bound. We conjecture that, under suitable assumptions, our bias-corrected estimator $\widehat{\tau}_{l}^{bc}$ is also doubly robust and semiparametrically efficient.
\begin{equation}
\label{eq:AIPW}
\begin{split}
\widehat{\tau}_l^{bc}&=
\frac{1}{n_{11}}\sum_{i\in\mathcal{I}_{11}} \left[
\Delta Y_{il}
-\widehat{\mu}_l^{01}(X_i)
-\widehat{\mu}_l^{10}(X_i)
+\widehat{\mu}_l^{00}(X_i)
\right]\\
&-\frac{1}{n_{11}}\sum_{i\in\mathcal{I}_{01}}
w_i\left[\Delta Y_{il}-\widehat{\mu}_l^{01}(X_i)\right]-\frac{1}{n_{11}}\sum_{i\in\mathcal{I}_{10}}
w_i\left[\Delta Y_{il}-\widehat{\mu}_l^{10}(X_i)\right]+\frac{1}{n_{11}}\sum_{i\in\mathcal{I}_{00}}
w_i\left[\Delta Y_{il}-\widehat{\mu}_l^{00}(X_i)\right].
\end{split}
\end{equation}
\end{remark}
\subsection{Variance Estimation}\label{subsection:varianceestimation}
Because our estimator is obtained from a weighted 3WFE regression, it is tempting to use the cluster-robust standard error from that regression. Previous studies (\citeinb{Abadie-Spiess_2022_JASA}; \citeinb{Liu_Vazquez_Bare_2026_wp}) show, however, that this standard error is inconsistent because it treats the matching weights as fixed. In this subsection, we propose a unified variance estimator that accounts for the matching step and is valid for all matching procedures within our framework.
For each untreated subgroup $(g,s)\in \{(0,1),(1,0),(0,0)\}$, let $\widehat{\Delta Y_{il}^{gs}}\equiv\sum_{j\in \mathcal{I}_{gs}}w_{ij}\Delta Y_{jl}$ denote the matched outcome change from subgroup $(g,s)$ for treated unit $i$, and let $\widetilde w_j^2\equiv\sum_{i\in \mathcal{I}_{11}}w_{ij}^2$ denote the sum of the squared pairwise weights received by untreated unit $j$. We construct our variance estimator by estimating each of the two components of the asymptotic variance in Theorem~\ref{theorem:asymptotic}, the effect-heterogeneity component, $V_{\tau(X)}$, and the conditional-variance component, $V_E$.
\begin{align*}
\widehat{V}_{\tau(X)} & = \frac{1}{n_{11}}\sum_{i\in \mathcal{I}_{11}}\left(\Delta Y_{il}- \widehat {\Delta Y_{il}^{01}}-\widehat {\Delta Y_{il}^{10}} + \widehat {\Delta Y_{il}^{00}}-\widehat{\tau}_l \right)^2, \\
\widehat{V}_E & =\frac{1}{n_{11}}\sum_{j\in \mathcal{I}_{01}} \left( w_j^2-\widetilde w_j^2\right)\widehat{\sigma}_{01,l}^2(X_j) + \frac{1}{n_{11}}\sum_{j\in \mathcal{I}_{10}} \left( w_j^2-\widetilde w_j^2 \right)\widehat{\sigma}_{10,l}^2(X_j) + \frac{1}{n_{11}}\sum_{j\in \mathcal{I}_{00}} \left( w_j^2-\widetilde w_j^2 \right)\widehat{\sigma}_{00,l}^2(X_j).
\end{align*}
Here, $\widehat{\sigma}_{gs,l}^2(X_j)$ is an estimator of the unknown conditional variance function, constructed following the strategy proposed by \cite{Abadie-Imbens_2006_ECMA}. For each untreated subgroup, we conduct an additional $K$-nearest-neighbor matching within the same subgroup and denote by $k_m(j)$ the index of the $m$-th nearest neighbor of unit $j$. The estimator $\widehat{\sigma}_{gs,l}^2(X_j)$ is then defined as:
\[
\widehat{\sigma}_{gs,l}^2(X_j)
=
\frac{K}{K+1}
\left(
\Delta Y_{jl}-\frac{1}{K}\sum_{m=1}^K \Delta Y_{k_m(j)l}
\right)^2.
\]
Theorem~\ref{theorem:varianceestimation} shows that our variance estimator is consistent under the same assumptions as Theorem~\ref{theorem:asymptotic}.
\begin{theorem}\label{theorem:varianceestimation}
Under Assumptions \ref{assumption:noanticipation},
\ref{assumption:overlap}, \ref{assumption:parallelgap}, \ref{assumption:iid}, \ref{assumption:matching}, and \ref{assumption:regularity},
\begin{align*}
\frac{\widehat V_{\tau(X)}+\widehat V_E}{V_{\tau(X)}+V_E } \xrightarrow{p} 1.
\end{align*}
\end{theorem}
Our variance estimation approach builds on that of \citet{Abadie-Imbens_2006_ECMA} for the cross-sectional setting but extends it in two directions. First, our variance estimator is tailored to the data structure and estimation procedure of the triple-differences design. This structure generates additional interaction terms between the residuals of different untreated subgroups. Our proof must handle these terms, which have no counterpart in the two-group cross-sectional setting. Second, the consistency proof for the variance estimator in \citet{Abadie-Imbens_2006_ECMA} relies on properties specific to nearest-neighbor matching with a fixed number of neighbors, so their proof cannot be applied directly here. In stead, we develop a refined proof strategy that is valid for the class of matching procedures characterized by Assumption~\ref{assumption:matching}, which includes nearest-neighbor matching with a fixed number of neighbors as a special case.
\section{Simulation}\label{section:simulation}
In this section, we provide the Monte Carlo evidence on (i) the inconsistency of the OLS estimators from the two 3WFE regression specifications studied in Section \ref{section:empiricaldiagnose} and (ii) the finite-sample performance of our proposed matched triple-differences estimators, their bias-corrected versions, and the proposed variance estimator.
We simulate a panel of $n=5,000$ units observed over four periods, $t=1,2,3,4$. Each unit belongs to one of four subgroups indexed by $(G,S)\in\{0,1\}^2$. The treatment is implemented simultaneously in period $t^*=3$ and only the exposed eligible subgroup, $(G,S)=(1,1)$, is treated.
Following the data generating process (DGP) of \cite{Otsu_Rai_2017_JASA}, we generate each unit's $q$-dimensional covariate vector $X_i=(X_{i1},\dots, X_{iq})'$:
\begin{align*}
X_{ij}&=\xi_i\frac{\left\lvert\zeta_{ij}\right\rvert}{\left\lVert\zeta_i\right\rVert},\quad \text{if } j=1,\dots,q \\
\xi_i&\sim \text{Uniform}(0,1), \quad \zeta_i\sim \mathcal{N}(0,I_q)
\end{align*}
Let $\left\lVertX\right\rVert$ denote the Euclidean norm of $X$. Under our DGP, $\left\lVertX\right\rVert$ follows a uniform distribution over $(0,1)$.
Each unit is assigned to one of the four subgroups according to a multinomial logit model:
\begin{align*}
P(G_i=g,S_i=s|X_i)=\frac{\exp(\alpha_{gs}+\gamma_{gs}(\left\lVertX_i\right\rVert-0.5))}{\sum_{(g,s)}\exp(\alpha_{gs}+\gamma_{gs}(\left\lVertX_i\right\rVert-0.5)) }
\end{align*}
The intercepts $\alpha_{gs}$ control the relative sizes of the four subgroups. We consider two sets of intercepts. The first gives the four subgroups approximately equal sizes: $\alpha_{11}=\alpha_{10}=\alpha_{01}=\alpha_{00}=0$. The second makes the exposed ineligible subgroup $(1,0)$ smaller than the other three subgroups: $\alpha_{11}=0, \alpha_{10}=-0.776, \alpha_{01}=0.228,\alpha_{00}=0.270 $. The targeted shares are approximately $25\%$, $11\%$, $30\%$, and $34\%$ for subgroups $(1,1)$, $(1,0)$, $(0,1)$, and $(0,0)$, respectively. This second set of intercepts is designed to mimic the subgroup sizes in our empirical application in Section \ref{section:empiricalapplication}.
The slopes $\gamma_{gs}$ determine how selection into the four subgroups depends on the covariates. We fix the slopes at $\gamma_{11}=0,\gamma_{10}=-2 ,\gamma_{01}=-0.5, \gamma_{00}=-3$. We choose these values so that selection on the covariate is stronger along the eligibility dimension than along the exposure dimension. This mirrors many triple-differences designs, in which the the policy is typically implemented along the exposure dimension (e.g., at the county or state level). Because the policy is usually implemented quasi-randomly, selection on individual covariates is weaker long the exposure dimension; along the eligibility dimension (e.g., growing tobacco, or being pregnant or not), selection on individual covariates is usually stronger.
The potential outcomes are defined as
\begin{align*}
Y_{it}(0)&=\alpha_i + 2 \beta_0 G_i(\left\lVertX_i\right\rVert-0.5)t +\beta_0 S_i (\left\lVertX_i\right\rVert-0.5)t +v_{it}(0)\\
Y_{it}(1)&=\tau + \alpha_i + 2 \beta_1 G_i(\left\lVertX_i\right\rVert-0.5)t +\beta_1 S_i (\left\lVertX_i\right\rVert-0.5)t +v_{it}(1)\\
\alpha_i&\sim \mathcal{N}(0,1), \quad v_{it}(0)\sim\mathcal{N}(0,1), \quad v_{it}(1)\sim \mathcal{N}(0,1)
\end{align*}
with $\alpha_i$, $v_{it}(0)$, $v_{it}(1)$, $X_i$ mutually independent. This specification satisfies the conditional parallel gaps assumption but violates the conditional parallel trends assumption. We set $\tau=1$ and consider two choices of $\beta_0$ and $\beta_1$. The first sets $\beta_0=\beta_1=2$, so the treatment effect is homogeneous with respect to $X$. The second sets $\beta_0=2,\beta_1=6$, so the treatment effect is heterogeneous with respect to $X$.
Based on the discussion above, we consider four main simulation designs. The first and second designs have equal subgroup sizes with a homogeneous/heterogeneous treatment effect. The third and fourth designs have a small $(1,0)$ subgroup size with a homogeneous/heterogeneous treatment effect.
We report results for $\tau_0$, the treatment effect in the implementation period $t^*=3$, based on $1,000$ repetitions with covariate dimension $q=2$. For each design, we report the OLS estimators from the two 3WFE regression specifications studied in Section \ref{section:empiricaldiagnose}, the matched triple-differences estimators under three matching methods, and their corresponding bias-corrected estimators. The three matching methods and their tuning parameters are: nearest-neighbor matching with one neighbor; nearest neighbor matching with a diverging number of neighbors, set to $n_{gs}^{1/3}$; and kernel matching with the Epanechnikov kernel and a bandwidth rule based on Stata's kmatch implementation, which sets the bandwidth to $1.5$ times the $90$th percentile of the non-zero pairwise distances from one-to-one nearest-neighbor matching with replacement. Bias correction uses a second-order series estimator. For each of the eight estimators, we report the Monte Carlo bias, the root mean squared error (RMSE), the Monte Carlo standard deviation, the average standard error computed by our proposed variance estimator, and the empirical coverage rate and length of the corresponding $95\%$ confidence interval.
The results reported in Tables \ref{table:design1} to \ref{table:design4} confirm our theoretical findings in Sections \ref{section:empiricaldiagnose} and \ref{section:matchedtriplediffereces}. In all four designs, the OLS estimators from the two 3WFE regression specifications exhibit substantial bias, ranging from $14\%$ to $81\%$ of the true treatment effect. By contrast, our proposed estimators have negligible bias in all four designs. For nearest-neighbor matching with a diverging number of neighbors, the bias-corrected estimator has smaller bias than the uncorrected one in all four designs. Our proposed variance estimator performs well: in almost all combinations of the four designs and six matching estimators, the average standard error is close to the Monte Carlo standard deviation, and the coverage rate of the $95\%$ confidence interval is close to the nominal rate. Nearest-neighbor matching with a diverging number of neighbors and the kernel matching are more efficient than nearest-neighbor matching with one neighbor, with smaller standard errors and shorter confidence intervals.
\begin{table}[htbp]
\centering
\caption{Balanced Subgroup Size, Homogeneous Treatment Effect}
\label{table:design1}
\small
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{RMSE} & \multicolumn{1}{c}{MCSD} & \multicolumn{1}{c}{Avg.~SE} & \multicolumn{1}{c}{Coverage} & \multicolumn{1}{c}{CI length} \\
\midrule
Cov. Trend ($X'\theta_t$) & $0.8084$ & $0.8131$ & $0.0879$ & $0.0874$ & $0.000$ & $0.3426$ \\
Cov. Level ($X'\theta$) & $0.7008$ & $0.7079$ & $0.1001$ & $0.0988$ & $0.000$ & $0.3875$ \\
\midrule
NN, $M = 1$ & $0.0060$ & $0.1203$ & $0.1202$ & $0.1225$ & $0.954$ & $0.4801$ \\
\quad bias-corrected & $0.0025$ & $0.1203$ & $0.1203$ & $0.1225$ & $0.953$ & $0.4801$ \\
\addlinespace
NN, $M = n_{gs}^{1/3}$ & $0.0238$ & $0.0964$ & $0.0935$ & $0.0946$ & $0.944$ & $0.3709$ \\
\quad bias-corrected & $0.0052$ & $0.0944$ & $0.0943$ & $0.0946$ & $0.950$ & $0.3709$ \\
\addlinespace
Kernel & $0.0086$ & $0.1029$ & $0.1026$ & $0.1040$ & $0.954$ & $0.4077$ \\
\quad bias-corrected & $0.0037$ & $0.1027$ & $0.1027$ & $0.1040$ & $0.954$ & $0.4077$ \\
\bottomrule
\end{tabular}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize
\emph{Notes:} This table reports the simulation results for design 1 with balanced subgroup size and homogeneous treatment effect. We run $1,000$ Monte Carlo replications, each with $n = 5,000$ units, $q = 2$ covariates, and $4$ periods. The treatment effect is $\tau_0 = 1$. The two OLS estimators are estimated using the two regression specifications studied in section \ref{section:empiricaldiagnose}. The three matching estimators use (i) nearest-neighbor matching with one neighbor, (ii) nearest-neighbor matching with a diverging number of neighbors, set to $n_{gs}^{1/3}$, and (iii) the kernel matching estimator with the Epanechnikov kernel with a bandwidth rule based on Stata's kmatch implementation, which is set to be $1.5$ times the $90$th percentile of non-zero pairwise distances from matching with replacement. The three bias-corrected matching estimators use the corresponding matching methods, and the bias correction is performed using a second-order series estimator.
\end{minipage}
\end{table}
\begin{table}[htbp]
\centering
\caption{Balanced Subgroup Size, Heterogeneous Treatment Effect}
\label{table:design2}
\small
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{RMSE} & \multicolumn{1}{c}{MCSD} & \multicolumn{1}{c}{Avg.~SE} & \multicolumn{1}{c}{Coverage} & \multicolumn{1}{c}{CI length} \\
\midrule
Cov. Trend ($X'\theta_t$) & $1.1133$ & $1.1472$ & $0.2768$ & $0.2768$ & $0.016$ & $1.0851$ \\
Cov. Level ($X'\theta$) & $0.6864$ & $0.7601$ & $0.3268$ & $0.3302$ & $0.454$ & $1.2943$ \\
\midrule
NN, $M = 1$ & $-0.0085$ & $0.2900$ & $0.2901$ & $0.2993$ & $0.957$ & $1.1731$ \\
\quad bias-corrected & $-0.0120$ & $0.2901$ & $0.2900$ & $0.2992$ & $0.958$ & $1.1727$ \\
\addlinespace
NN, $M = n_{gs}^{1/3}$ & $0.0093$ & $0.2839$ & $0.2839$ & $0.2894$ & $0.949$ & $1.1343$ \\
\quad bias-corrected & $-0.0092$ & $0.2834$ & $0.2834$ & $0.2888$ & $0.951$ & $1.1320$ \\
\addlinespace
Kernel & $-0.0058$ & $0.2856$ & $0.2857$ & $0.2921$ & $0.950$ & $1.1451$ \\
\quad bias-corrected & $-0.0107$ & $0.2856$ & $0.2856$ & $0.2920$ & $0.948$ & $1.1446$ \\
\bottomrule
\end{tabular}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize
\emph{Notes:} This table reports the simulation results for design 2 with balanced subgroup size and heterogeneous treatment effect. We run $1,000$ Monte Carlo replications, each with $N = 5,000$ units, $q = 2$ covariates, and $4$ periods. The treatment effect is $\tau_0 = 4.87$. The two OLS estimators are estimated using the two regression specifications studied in section \ref{section:empiricaldiagnose}. The three matching estimators use (i) nearest-neighbor matching with one neighbor, (ii) nearest-neighbor matching with a diverging number of neighbors, set to $n_{gs}^{1/3}$, and (iii) the kernel matching estimator with the Epanechnikov kernel with a bandwidth rule based on Stata's kmatch implementation, which is set to be $1.5$ times the $90$th percentile of non-zero pairwise distances from matching with replacement. The three bias-corrected matching estimators use the corresponding matching methods, and the bias correction is performed using a second-order series estimator.
\end{minipage}
\end{table}
\begin{table}[htbp]
\centering
\caption{Small $(1,0)$ Subgroup, Homogeneous Treatment Effect}
\label{table:design3}
\small
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{RMSE} & \multicolumn{1}{c}{MCSD} & \multicolumn{1}{c}{Avg.~SE} & \multicolumn{1}{c}{Coverage} & \multicolumn{1}{c}{CI length} \\
\midrule
Cov. Trend ($X'\theta_t$) & $0.8118$ & $0.8192$ & $0.1102$ & $0.1102$ & $0.000$ & $0.4318$ \\
Cov. Level ($X'\theta$) & $0.7087$ & $0.7196$ & $0.1253$ & $0.1242$ & $0.001$ & $0.4869$ \\
\midrule
NN, $M = 1$ & $0.0134$ & $0.1537$ & $0.1532$ & $0.1496$ & $0.951$ & $0.5866$ \\
\quad bias-corrected & $0.0011$ & $0.1541$ & $0.1542$ & $0.1496$ & $0.947$ & $0.5865$ \\
\addlinespace
NN, $M = n_{gs}^{1/3}$ & $0.0461$ & $0.1386$ & $0.1307$ & $0.1257$ & $0.922$ & $0.4927$ \\
\quad bias-corrected & $-0.0003$ & $0.1345$ & $0.1346$ & $0.1257$ & $0.942$ & $0.4927$ \\
\addlinespace
Kernel & $0.0206$ & $0.1409$ & $0.1394$ & $0.1353$ & $0.942$ & $0.5303$ \\
\quad bias-corrected & $0.0020$ & $0.1407$ & $0.1408$ & $0.1353$ & $0.941$ & $0.5303$ \\
\bottomrule
\end{tabular}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize
\emph{Notes:} This table reports the simulation results for design 3 with a small $(1,0)$ subgroup and homogeneous treatment effect. We run $1,000$ Monte Carlo replications, each with $n = 5,000$ units, $q = 2$ covariates, and $4$ periods. The treatment effect is $\tau_0 = 1$. The two OLS estimators are estimated using the two regression specifications studied in section \ref{section:empiricaldiagnose}. The three matching estimators use (i) nearest-neighbor matching with one neighbor, (ii) nearest-neighbor matching with a diverging number of neighbors, set to $n_{gs}^{1/3}$, and (iii) the kernel matching estimator with the Epanechnikov kernel with a bandwidth rule based on Stata's kmatch implementation, which is set to be $1.5$ times the $90$th percentile of non-zero pairwise distances from matching with replacement. The three bias-corrected matching estimators use the corresponding matching methods, and the bias correction is performed using a second-order series estimator.
\end{minipage}
\end{table}
\begin{table}[htbp]
\centering
\caption{Small $(1,0)$ Subgroup, Heterogeneous Treatment Effect}
\label{table:design4}
\small
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{1}{c}{Bias} & \multicolumn{1}{c}{RMSE} & \multicolumn{1}{c}{MCSD} & \multicolumn{1}{c}{Avg.~SE} & \multicolumn{1}{c}{Coverage} & \multicolumn{1}{c}{CI length} \\
\midrule
Cov. Trend ($X'\theta_t$) & $1.2256$ & $1.2709$ & $0.3361$ & $0.3394$ & $0.055$ & $1.3304$ \\
Cov. Level ($X'\theta$) & $0.6947$ & $0.7501$ & $0.2830$ & $0.2883$ & $0.336$ & $1.1300$ \\
\midrule
NN, $M = 1$ & $-0.0005$ & $0.2714$ & $0.2715$ & $0.2708$ & $0.949$ & $1.0615$ \\
\quad bias-corrected & $-0.0128$ & $0.2722$ & $0.2720$ & $0.2705$ & $0.946$ & $1.0604$ \\
\addlinespace
NN, $M = n_{gs}^{1/3}$ & $0.0322$ & $0.2613$ & $0.2595$ & $0.2590$ & $0.954$ & $1.0152$ \\
\quad bias-corrected & $-0.0142$ & $0.2611$ & $0.2608$ & $0.2579$ & $0.954$ & $1.0111$ \\
\addlinespace
Kernel, & $0.0067$ & $0.2632$ & $0.2633$ & $0.2632$ & $0.951$ & $1.0316$ \\
\quad bias-corrected & $-0.0119$ & $0.2640$ & $0.2638$ & $0.2629$ & $0.953$ & $1.0304$ \\
\bottomrule
\end{tabular}
\par\smallskip
\begin{minipage}{\linewidth}\footnotesize
\emph{Notes:} This table reports the simulation results for design 4 with a small $(1,0)$ subgroup and heterogeneous treatment effect. We run $1,000$ Monte Carlo replications, each with $n = 5,000$ units, $q = 2$ covariates, and $4$ periods. The treatment effect is $\tau_0 = 3.13$. The two OLS estimators are estimated using the two regression specifications studied in section \ref{section:empiricaldiagnose}. The three matching estimators use (i) nearest-neighbor matching with one neighbor, (ii) nearest-neighbor matching with a diverging number of neighbors, set to $n_{gs}^{1/3}$, and (iii) the kernel matching estimator with the Epanechnikov kernel with a bandwidth rule based on Stata's kmatch implementation, which is set to be $1.5$ times the $90$th percentile of non-zero pairwise distances from matching with replacement. The three bias-corrected matching estimators use the corresponding matching methods, and the bias correction is performed using a second-order series estimator.
\end{minipage}
\end{table}
\newpage
\section{Empirical Application}\label{section:empiricalapplication}
In this section, we revisit our motivating example from \cite{Cai_2016_AEJ}, who studies how the introduction of an agricultural insurance program affects rural households' borrowing and saving behavior in $12$ tobacco-producing counties in Jiangxi Province, China. In these counties, some rural households rely on tobacco cultivation as their main source of income (hereafter, tobacco households), while others rely on other crops such as rice (hereafter, non-tobacco households). In 2003, one of the 12 counties, Guangchang County (hereafter, the exposed county), collaborated with the People's Insurance Company of China (PICC) to design and launch a tobacco production insurance program. Under the program, a household was eligible for an insurance payout if a covered weather disaster reduced its tobacco yield by at least $30\%$. Only tobacco households in the exposed county were offered the insurance; non-tobacco households in the exposed county and households in the other counties were not eligible to purchase it. Because purchase was compulsory, all tobacco households in the exposed county were covered.
\cite{Cai_2016_AEJ} uses a household-level panel of 5,746 households from the 12 counties, including more than 3,000 tobacco households (approximately $1,200$ of them in the exposed county) and about $2,200$ non-tobacco households. The dataset has two parts. The first part is an administrative dataset on household saving and borrowing information from the Rural Credit Cooperative (RCC), the main rural bank in China. The second part is RCC annual survey data on the background information of those households such as demographic information about household heads, agricultural land allocation information, and revenue from the crop production. The panel covers the years 2001 to 2008; because the insurance program was implemented in 2003, it contains two pre-treatment years (2001--2002), and six post-treatment years (2003--2008).
\cite{Cai_2016_AEJ} uses a triple-differences design to estimate the effect of the program on household borrowing and saving. As discussed in our motivating example in Section \ref{section:econometricssetup}, a triple-differences design accounts for both the county-specific and industry-specific shocks, either of which could confound a standard difference-in-differences design. In our empirical application, we focus on the dynamic effects of insurance provision on the saving behavior of tobacco households in the exposed county. We use the same three saving outcomes as \cite{Cai_2016_AEJ}: net saving, defined as the annual increase in total savings; the saving rate, defined as the ratio of net saving to current household income; and the flexible-term saving ratio, defined as the ratio of net saving in checking accounts to total net saving.
We first analyze the three outcomes using an event-study regression similar to that in \cite{Cai_2016_AEJ}, which corresponds to the covariate-level 3WFE specification in Section \ref{section:empiricaldiagnose}:
\begin{align*}
Y_{it}
= \lambda_{G_iS_i}
+ \lambda_{G_i t}
+ \lambda_{S_i t}
+ \sum_{l\neq -1}
\tau_l^{CL}\, G_iS_i\mathbbm{1}(t-t^*=l)
+ X_i'\theta
+ u_{it},
\end{align*}
where $\lambda_{G_iS_i}$, $\lambda_{G_i t}$, and $\lambda_{S_i t}$ are subgroup, exposed-county-by-time, and tobacco-by-time fixed effects, respectively, and $X_i$ is a vector of time-invariant household characteristics, including the age and education level of the household head and household size.
We then implement our matched triple-differences estimator and construct confidence intervals using the variance estimator proposed in Section \ref{section:matchedtriplediffereces}. We match exactly on the education level of the household head and treat the head's age and household size as continuous variables because they take many values. Our matching method is bias-corrected kernel matching with an Epanechnikov kernel. The bandwidth follows the rule implemented in Stata's kmatch: $1.5$ times the $90$th percentile of the non-zero pairwise distances from one-to-one nearest-neighbor matching with replacement. The bias correction uses a second-order series estimator. For both the 3WFE regression and the matching estimator, the baseline year is $2002$, the year before the program was implemented.
Figure \ref{figure:empiricalcaikernelbc} plots the event-study estimators from the 3WFE regression and kernel matching\footnote{In appendix \ref{appendix:additionalempirical}, we also report the estimates from the kernel matching estimator without bias correction and bias-corrected nearest-neighbor matching estimator with a diverging number of neighbors. The main results are not sensitive to the choice of matching methods and whether the bias-correction is applied.}. For net saving (Panel A), the two estimators lead to qualitatively different conclusions. Specifically, the 3WFE regression finds an sign-reversing pattern that the program increases net saving at years $0$ and $3$ and then decreases it at year $5$. In contrast, the kernel matching estimator finds no statistically significant effect on net saving in any post-treatment years, which is more plausible. Recall from Section \ref{section:empiricaldiagnose} that this 3WFE specification is consistent if and only if the parallel gaps assumption holds unconditionally. This suggests that, for the net saving, the parallel gaps assumption may be more plausible after conditioning on household characteristics. For the other two outcomes (Panels B and C), the kernel matching estimator and the 3WFE regression estimator yield similar results, although the kernel matching estimator finds larger effects on the flexible-term saving ratio in years $2$ and $3$. Overall, the matched triple-differences estimates can differ from the 3WFE estimates both in magnitude and in their qualitative conclusions.
\begin{figure}[p]
\centering
\setstretch{1}
\setlength{\parindent}{0pt}
\caption{Event-study estimates: 3WFE regression and matched triple-differences estimator}
\label{figure:empiricalcaikernelbc}
\vspace{0.04in}
\begin{tikzpicture}[every node/.style={font=\footnotesize}]
\begin{axis}[
width=0.78\linewidth, height=0.24\textheight,
xmin=-2.35, xmax=5.35, ymin=-0.40027801463779, ymax=0.34504899909701,
xtick={-2,-1,0,1,2,3,4,5},
xlabel={Event time (year $-$ 2003)}, ylabel={Treatment effect},
axis lines=left, axis line style={black,thin},
tick align=outside, tick style={black,thin},
tick label style={font=\footnotesize},
label style={font=\footnotesize},
xlabel style={at={(axis description cs:0.5,-0.17)},anchor=north},
ylabel style={at={(axis description cs:-0.11,0.5)},anchor=south},
scaled y ticks=false,
yticklabel style={/pgf/number format/fixed,/pgf/number format/precision=3},
legend style={at={(1.02,0.5)},anchor=west,draw=none,
font=\footnotesize,inner sep=1pt},
legend cell align=left, legend columns=1,
clip=true,
]
\addplot[black,thin,no marks,forget plot] coordinates {(-2.35,0) (5.35,0)};
\addplot[black,dashed,thin,no marks,forget plot] coordinates {(-1,-0.40027801463779) (-1,0.34504899909701)};
\addplot[
only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-2.065,0.03278396501074) += (0,0.033038353351182) -= (0,0.033038353351182)
(-0.065,0.07723661597553) += (0,0.051545595431263) -= (0,0.051545595431263)
(0.935,-0.043751965344329) += (0,0.040673623325812) -= (0,0.040673623325812)
(1.935,0.032237012146934) += (0,0.053347809186021) -= (0,0.053347809186021)
(2.935,0.17860639235518) += (0,0.099240990749347) -= (0,0.099240990749347)
(3.935,0.010624681265614) += (0,0.084083444079536) -= (0,0.084083444079536)
(4.935,-0.19722942837161) += (0,0.12107520493823) -= (0,0.12107520493823)
};
\addlegendentry{Regression}
\addplot[forget plot,only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153}] coordinates {(-1.065,0)};
\addplot[
only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-1.935,-0.0072926863016265) += (0,0.038654032950374) -= (0,0.038654032950374)
(0.065,0.0033408336817603) += (0,0.052404924595985) -= (0,0.052404924595985)
(1.065,-0.0019253958329706) += (0,0.047644232644645) -= (0,0.047644232644645)
(2.065,0.015527154522128) += (0,0.063152640502042) -= (0,0.063152640502042)
(3.065,-0.012531710609) += (0,0.090321288477436) -= (0,0.090321288477436)
(4.065,-0.056798746382658) += (0,0.078135846351844) -= (0,0.078135846351844)
(5.065,-0.075814057075813) += (0,0.11710117087035) -= (0,0.11710117087035)
};
\addlegendentry{Matched DDD}
\addplot[forget plot,only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44}] coordinates {(-0.935,0)};
\end{axis}
\end{tikzpicture}\par
\vspace{-0.02in}
{\footnotesize Panel A: Net saving}\par
\vspace{0.05in}
\begin{tikzpicture}[every node/.style={font=\footnotesize}]
\begin{axis}[
width=0.78\linewidth, height=0.24\textheight,
xmin=-2.35, xmax=5.35, ymin=-0.13464693760165, ymax=0.066716038222654,
xtick={-2,-1,0,1,2,3,4,5},
xlabel={Event time (year $-$ 2003)}, ylabel={Treatment effect},
axis lines=left, axis line style={black,thin},
tick align=outside, tick style={black,thin},
tick label style={font=\footnotesize},
label style={font=\footnotesize},
xlabel style={at={(axis description cs:0.5,-0.17)},anchor=north},
ylabel style={at={(axis description cs:-0.11,0.5)},anchor=south},
scaled y ticks=false,
yticklabel style={/pgf/number format/fixed,/pgf/number format/precision=3},
legend style={at={(1.02,0.5)},anchor=west,draw=none,
font=\footnotesize,inner sep=1pt},
legend cell align=left, legend columns=1,
clip=true,
]
\addplot[black,thin,no marks,forget plot] coordinates {(-2.35,0) (5.35,0)};
\addplot[black,dashed,thin,no marks,forget plot] coordinates {(-1,-0.13464693760165) (-1,0.066716038222654)};
\addplot[
only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-2.065,-0.0083879405501005) += (0,0.023449453918368) -= (0,0.023449453918368)
(-0.065,0.023063620532447) += (0,0.022853104228507) -= (0,0.022853104228507)
(0.935,0.014396292948598) += (0,0.025536221754226) -= (0,0.025536221754226)
(1.935,-0.032938609492962) += (0,0.030063560432264) -= (0,0.030063560432264)
(2.935,-0.024354776486304) += (0,0.028987026928152) -= (0,0.028987026928152)
(3.935,-0.04001878515996) += (0,0.03117962236484) -= (0,0.03117962236484)
(4.935,-0.065892770571912) += (0,0.030702135926055) -= (0,0.030702135926055)
};
\addlegendentry{Regression}
\addplot[forget plot,only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153}] coordinates {(-1.065,0)};
\addplot[
only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-1.935,-0.0059915426586062) += (0,0.023468253643622) -= (0,0.023468253643622)
(0.065,0.0148217072547) += (0,0.02394563832467) -= (0,0.02394563832467)
(1.065,0.01115741914584) += (0,0.026431346625161) -= (0,0.026431346625161)
(2.065,-0.039257319286807) += (0,0.03176410827417) -= (0,0.03176410827417)
(3.065,-0.02801425477034) += (0,0.031174009477146) -= (0,0.031174009477146)
(4.065,-0.04102811553772) += (0,0.031813313816065) -= (0,0.031813313816065)
(5.065,-0.066929290403756) += (0,0.031769118306759) -= (0,0.031769118306759)
};
\addlegendentry{Matched DDD}
\addplot[forget plot,only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44}] coordinates {(-0.935,0)};
\end{axis}
\end{tikzpicture}\par
\vspace{-0.02in}
{\footnotesize Panel B: Saving rate}\par
\vspace{0.05in}
\begin{tikzpicture}[every node/.style={font=\footnotesize}]
\begin{axis}[
width=0.78\linewidth, height=0.24\textheight,
xmin=-2.35, xmax=5.35, ymin=-0.12680145612051, ymax=0.24111493392685,
xtick={-2,-1,0,1,2,3,4,5},
xlabel={Event time (year $-$ 2003)}, ylabel={Treatment effect},
axis lines=left, axis line style={black,thin},
tick align=outside, tick style={black,thin},
tick label style={font=\footnotesize},
label style={font=\footnotesize},
xlabel style={at={(axis description cs:0.5,-0.17)},anchor=north},
ylabel style={at={(axis description cs:-0.11,0.5)},anchor=south},
scaled y ticks=false,
yticklabel style={/pgf/number format/fixed,/pgf/number format/precision=3},
legend style={at={(1.02,0.5)},anchor=west,draw=none,
font=\footnotesize,inner sep=1pt},
legend cell align=left, legend columns=1,
clip=true,
]
\addplot[black,thin,no marks,forget plot] coordinates {(-2.35,0) (5.35,0)};
\addplot[black,dashed,thin,no marks,forget plot] coordinates {(-1,-0.12680145612051) (-1,0.24111493392685)};
\addplot[
only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-2.065,-0.031466247576672) += (0,0.040632916047508) -= (0,0.040632916047508)
(-0.065,0.006780677822692) += (0,0.04200750436697) -= (0,0.04200750436697)
(0.935,0.030510667282213) += (0,0.04027750711002) -= (0,0.04027750711002)
(1.935,0.049415648042741) += (0,0.048055316651295) -= (0,0.048055316651295)
(2.935,0.045013756686105) += (0,0.045749248829915) -= (0,0.045749248829915)
(3.935,0.070427588164144) += (0,0.052713895580232) -= (0,0.052713895580232)
(4.935,0.14103572272288) += (0,0.053576951374135) -= (0,0.053576951374135)
};
\addlegendentry{Regression}
\addplot[forget plot,only marks,mark=*,mark size=2.1pt,color={rgb,255:red,0;green,76;blue,153}] coordinates {(-1.065,0)};
\addplot[
only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44},
error bars/.cd,y dir=both,y explicit,
error bar style={line width=0.65pt},error mark options={rotate=90,mark size=2.3pt,line width=0.65pt}
] coordinates {
(-1.935,-0.026180834651252) += (0,0.047448326311303) -= (0,0.047448326311303)
(0.065,0.0055140263956653) += (0,0.05200625198073) -= (0,0.05200625198073)
(1.065,0.034105963281522) += (0,0.043579512138134) -= (0,0.043579512138134)
(2.065,0.082572994502872) += (0,0.050482823855256) -= (0,0.050482823855256)
(3.065,0.073966511152735) += (0,0.052421800101064) -= (0,0.052421800101064)
(4.065,0.081444440576871) += (0,0.055482952567311) -= (0,0.055482952567311)
(5.065,0.13609657334295) += (0,0.056392364553265) -= (0,0.056392364553265)
};
\addlegendentry{Matched DDD}
\addplot[forget plot,only marks,mark=triangle*,mark size=2.1pt,color={rgb,255:red,190;green,52;blue,44}] coordinates {(-0.935,0)};
\end{axis}
\end{tikzpicture}\par
\vspace{-0.02in}
{\footnotesize Panel C: Flexible-term saving ratio}\par
\vspace{0.05in}
\vspace{0.02in}
\begin{minipage}{\linewidth}
\scriptsize
\textit{Notes:} This figure plots event-study estimates from the covariate-level 3WFE regression and from the matched triple-differences estimator with bias-corrected kernel matching. Event time $-1$ (2002) is the reference period. Bars show 95\% confidence intervals. The covariates included in the analysis are: the age and education level of the household head and household size. Cross-sectional sample sizes for panels A, B, and C are 3,736 (29,888 observations), 3,332 households (26,656 observations), and 3,545 households (28,360 observations), respectively. Matching uses an Epanechnikov kernel and bias correction uses the second-order series estimator.
\par
\vspace{0.25\baselineskip}
\end{minipage}
\end{figure}
\newpage
\section{Conclusion}\label{section:conclusion}
We formally analyze two commonly used empirical strategies in triple-differences designs: augmenting the 3WFE regression specification with either covariate trends or covariate levels. These strategies are often used when the parallel gaps assumption is believed to hold conditionally on covariates. Our analysis shows that both generally lead to inconsistent estimators of the ATT. First, we show that the weights of the regression estimands for both specifications fail to align the covariate distribution of each untreated subgroup with that of the treated subgroup so neither satisfies the covariate balancing condition under which the estimands would equal the ATT. Second, we characterize the additional assumptions under which each estimand equals the ATT. For the covariate-trend specification, the sufficient conditions require the treatment effect to be homogeneous in the covariates, together with a strong parametric assumption on either the conditional mean function or the propensity score. For the covariate-level specification, the necessary and sufficient condition is that the parallel gaps assumption holds unconditionally, in which case there is no need to add covariates.
To address these limitations, we propose a matched triple-differences framework for valid estimation and inference. Our two-step estimation procedure first conducts three separate pairwise matchings to satisfy the covariate balancing condition and then runs a weighted 3WFE regression, which yields the matched triple-differences estimator. We then develop a novel asymptotic framework to establish the consistency and asymptotic normality of our estimator for a class of matching procedures, including kernel matching and nearest neighbor matching with a fixed or diverging number of neighbors. In addition, we propose a bias-corrected matched triple-differences estimator and a valid variance estimator that accounts for the matching step. Monte Carlo simulations confirm our theoretical results and show that the proposed estimators and variance estimator perform well in finite samples. Revisiting \cite{Cai_2016_AEJ}, we find that, compared with the regression specification used in that paper, our matched triple-differences estimators can yield treatment effect estimates of a different magnitude and lead to different qualitative conclusions.
For future work, we envision extending our proposed matching estimators to settings with staggered adoption and to more flexible sampling schemes, such as repeated cross sections. Another interesting direction is to adapt recent developments in the matching literature, such as \cite{Lin-Ding-Han_2023_ECMA} and \cite{Lin_Han_2025_JOE}, to establish the double robustness of our proposed bias-corrected matching estimator and to characterize when it attains the semiparametric efficiency bound.
\bibliography{matching_DDD_references}
\newpage
\section*{Appendix}