Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,849 characters · 16 sections · 30 citation commands
Semiparametric Triple Difference Estimators
The triple difference framework for causal inference extends the well-known difference-in-differences (DiD) framework ashenfelter1984using,card1990impact, card1994minimum,heckman1997matching, card2000minimum,abadie2005semiparametric by incorporating data from an auxiliary domain. Specifically, it considers a setting with two domains: a target domain, in which the causal parameter of interest is defined, and a reference domain. The triple difference framework is designed for settings where the DiD assumptions do not necessarily hold in the target domain. It characterizes a set of requirements on the relationship between the two domains under which identifiability is attainable by fusing data from the two domains olden2022triple,zhuang2024way,berck2016note.
The triple difference framework relaxes the parallel trends assumption, which is fundamental to the canonical DiD, and thus allows for a more flexible identification scheme. Given its flexibility and practical relevance, the triple difference framework has been widely adopted in economics and policy evaluation---see olden2022triple for a survey. However, formal identification and estimation theory for triple difference has received little attention. Although an identification formula for triple difference appears in wooldridge2020introductory and frohlich2019impact, others have settled on incorporating the high-level ideas without technical details lechner2011estimation, angrist2009mostly. Recently, olden2022triple conducted the first formal study of identification in the triple difference framework. They outlined the necessary assumptions for the identification of the average treatment effect on the treated (ATT) in the target domain based on outcome regression. Beyond that, however, key aspects such as weighting-based identification, and the design of robust and efficient estimators for ATT or other causal quantities have remained largely unexplored. This work aims to fill these gaps. We consider both identification and estimation aspects of the problem, as described below.
From the identification point of view, we contribute to the triple difference literature by (i) formalizing the outcome regression-based identification for the repeated cross-sections setting (described in (ref)); and (ii) establishing weighting methods to identify the ATT, which were missing in the literature, for both panel data (also known as repeated outcomes) and repeated cross-sections settings. In the repeated cross-sections setting, we avoid the widely adopted assumption of no compositional changes, i.e., the assumption that the composition of the units is time-invariant. Such an assumption requires us to sample observations from the same population, which is unrealistic in most scenarios hong2013measuring, sant2025difference. We show that this assumption, although commonly made in the literature, is not required for identification.
From the estimation point of view, we develop semiparametric estimators for the triple difference framework in both panel data and repeated cross-sections settings. We derive influence function-based estimators that attain the semiparametric efficiency bound. We demonstrate that they are doubly robust, in the sense that they remain consistent as long as either the estimator for outcome regression or the one for treatment assignment is consistent, but not necessarily both. Furthermore, we characterize the conditions under which our proposed estimators are $\sqrt{n}$-consistent and asymptotically normal even if the nuisance function estimators do not converge at $\sqrt{n}$ rate. In the repeated cross-sections setting, we propose estimators for when compositional changes assumption holds as well as when it is violated. Additionally, we discuss the implications of this assumption in terms of the robustness properties that our estimators can guarantee.
\paragraph{Related literature.} The triple difference framework is an extension of the canonical DiD framework. A large body of work has contributed to refining the latter. abadie2005semiparametric developed weighting estimators and introduced semiparametric approaches to estimation of ATT in the DiD framework. More recently, sant2020doubly proposed a set of doubly robust estimators, which are robust for inference under parametric assumptions. Others have proposed methods to combine DiD with synthetic controls arkhangelsky2021synthetic, explored heterogeneous treatment effects nie2019nonparametric, provided efficient estimation methods for discrete-valued outcomes li2019double, analyzed DiD with continuous treatments callaway2024difference, and extended DiD to nonlinear models athey2006identification, torous2024optimal. Yet, similar developments for triple difference remain limited zhuang2024way, akbari2024triple. Our work extends identification and estimation techniques for the triple difference framework building on the aforementioned literature.
The triple difference framework can be naturally interpreted as a data fusion framework bareinboim2016causal, degtiar2023review,colnet2024causal,yang2024causal, where the goal is to identify causal effects in a target domain using additional information from an auxiliary or reference domain. Despite the conceptual similarity, the data fusion literature has largely focused on conditional independence-type assumptions on potential outcome variables and has paid little attention to parallel trends-type assumptions. As such, existing data fusion frameworks do not directly address the type of identification problems encountered in the triple difference setting.
As mentioned earlier, in the literature of canonical DiD, when repeated cross-sectional data are available, it is commonly assumed that the covariates and the treatment assignment are time-invariant, known as the {no compositional changes} assumption heckman1997matching, abadie2005semiparametric, sant2020doubly, callaway2021difference. To the best of our knowledge, hong2013measuring was the first to analyze DiD without imposing this assumption, at the expense of replacing the parallel trends requirement with a selection on observables condition. zimmert2018efficient showed that ATT can be identified under a set of weaker independence assumptions and, only recently, sant2025difference proved that identification is possible without making such an assumption altogether. We generalize the viewpoint of sant2025difference to the triple difference setting.
\paragraph{Organization.} The rest of this paper is organized as follows: we formally introduce the triple difference setup and our parameters of interest in (ref). We provide the identification results in (ref), and in (ref), we describe our proposed estimators. In all sections, we present the results for the panel data setting first, followed by those for the repeated cross-sections setting; proofs of all results are provided in the Appendix. In (ref), we evaluate our proposed estimators on synthetic data. In (ref), we apply our proposed methodology to study the effect of mandated maternity benefits on the hourly wages of women of childbearing age.
We assume that data are collected from two domains indicated by $D\in\{0,1\}$, where $D=0$ and $D=1$ represent the reference and target domains, respectively. Units in each domain are divided into two groups, indicated by $G\in\{0,1\}$. Members of the groups $G=0$ and $G=1$ are ineligible and eligible to receive a specific treatment, respectively. For instance, if only certain people, such as those in lower quantiles of income, are eligible for a government financial assistance, the lower quantiles of income correspond to $G=1$. Let $A\in\{0,1\}$ denote a binary treatment (or policy) which is assigned only to the eligible units in the target domain. We denote by $X$ a set of observed pre-treatment covariates of the units. In what follows, we briefly describe the panel data and repeated cross-sections settings.
In this setting, data are collected from the same population at two time points indicated by the index $t\in\{0,1\}$. Let $Y_0, Y_1\in\mathbb{R}$ represent the observed outcomes at time $t=0$ and $t=1$, respectively. Moreover, for $t\in\{0,1\}$, let $Y_t^{(0)}$, $Y_t^{(1)}\in\mathbb{R}$ denote the potential outcomes of $Y_t$, if (potentially contrary to the fact) the treatment was set to $A=0$ (control) and $A=1$ (treatment), respectively. The treatment is administered after the outcomes are measured at time $t=0$, only in the target domain, and to all eligible units, represented by $(D=1,G=1)$. We observe i.i.d. samples $\{O_i=(G,D,X,Y_0,Y_1)_i\}_{i=1}^n$. Our parameter of interest is \[\tau_\mathrm{pd}\coloneqq \ex{}\left[Y_1^{(1)}-\Y{1}\:\big|\: G=1, D=1\right].\] Note that the parameter $\tau_\mathrm{pd}$ represents the average treatment effect at time $t=1$ in the population corresponding to $(D=1,G=1)$, which is the treated population. Therefore, the parameter $\tau_\mathrm{pd}$ is the ATT in the target domain.
In this setting, cross-sectional data are available at two time points. We denote by $T\in\{0,1\}$ the (random) collection time of samples, where $T=0$ and $T=1$ correspond to pre- and post-treatment periods, respectively. Let $Y$ denote the observed outcome variable, and $Y^{(0)}, Y^{(1)}\in\mathbb{R}$ represent the potential outcomes under control and treatment, respectively. The treatment is administered to all units indicated by $(G=1,D=1,T=1)$, that is, the eligible group in the target domain at time $T=1$. We observe i.i.d. samples $\{O_i=(G,D,T,X,Y)_i\}_{i=1}^n$. Note that unlike most work in the literature, we do not require the no compositional changes assumption, which requires $p(X,G,D\mid T=0)=p(X,G,D\mid T=1)$. The parameter of interest is defined as \[ \tau_\mathrm{rc}\coloneqq \ex{}\left[Y^{(1)}-\Y{}\:\big|\: G=1,D=1,T=1\right], \] which is the ATT in the target domain.
In the panel data setting, we track the same individuals (or units) across the two time periods. This allows researchers to observe changes within individuals over time. In contrast, in the repeated cross-sections setting, we are given different samples at the two time points. Each cross-section represents a snapshot of one of the populations at a given moment, but the individuals in each sample are not necessarily the same. Therefore, the repeated cross-sections setting can introduce additional complications when composition of the units changes across the two populations. Specifically, the researcher should disentangle whether observed differences between time periods are due to actual changes within individuals or shifts in the composition of the population. We will address the complications due to compositional changes in the task of identification and estimation in Sections (ref) and (ref), respectively.
In this section, we present our identification results, first for the panel data setting, followed by the repeated cross-sections setting.
We require consistency and no anticipation effects assumptions in our setting. The former links the potential outcomes to the observed outcomes, while the latter requires that receiving treatment has no effect before its actual implementation. Both of these assumptions are commonly made in the literature of the DiD framework.
We also require positivity, also known as overlap, in our setting, which is another commonly made assumption in the literature of causal inference. This assumption ensures that (i) at least some units are treated in the target domain, to make the parameter of interest well-defined, and (ii) in each stratum of $X$, there is a non-zero probability of being assigned to each group.
We say strict positivity holds if $\epsilon>0$. As we shall see, positivity is sufficient for identification, whereas inference requires strict positivity. Finally, we make the following assumption.
(ref) is a relaxation of the conditional parallel trends assumption in the canonical DiD framework, in the sense that when data from one domain are available, say $D=1$, conditional parallel trends asserts that the left hand side of Equation (ref) is equal to $0$. Rather than assuming that the left hand side term is zero, here, we assume that it is equal to its counterpart in the reference domain $D=0$.
We next present our nonparametric identification result for the ATT based on outcome regression (OR).
The first statement in (ref) provides the identification result for the conditional ATT only for $x$ where $p(G,D\mid X=x)$ is positive. (ref) ensures that this condition always holds for all $x$ observed in the subpopulation of treated units, $(G=1,D=1)$, resulting in the identification of ATT as given by the second statement in (ref). Note that $\psi_\mathrm{pd}$ is a functional of the observational law. In Section (ref), we will consider $\psi_\mathrm{pd}$ as our parameter of interest.
We next present an identification result for ATT in the triple difference framework based on inverse propensity weighting. As in the case of OR-based identification, we first identify the conditional ATT in the target domain, and subsequently use it to identify the ATT.
As in the case of panel data setting, we begin by presenting the necessary assumptions for identification. Specifically, (ref) establishes a link between the potential outcomes and the observed outcomes, (ref) rules out anticipation effects, (ref) ensures that the parameter of interest is well-defined, and (ref) ensures the identifiability of the ATT. These assumptions are the counterparts of Assumptions (ref) through (ref) in the repeated cross-sections setting.
We will say strict positivity holds when $\epsilon>0$.
Similar to the previous subsection, we first present the OR-based identification result.
As in the previous case, $\psi_\mathrm{rc}$ is a functional of the observational law, and will be considered as the parameter of interest in (ref). We next turn our focus to weighting-based identification for the repeated cross-sections setting.
Importantly, (ref) allows for compositional changes. That is, the composition of the units is allowed to change arbitrarily across time points. However, the commonly made no compositional changes assumption simplifies the identification. Below, we formally present this assumption and discuss its implications for identification of the ATT.
Under Assumption (ref), $\phi_0(\cdot)$ in (ref) can be further simplified to: \[
\] where $\rho_0(\cdot)$ is defined in Equation (ref). That is, the identification will now require the less complex nuisance function $\pi_{g,d}(\cdot)$, as opposed to $\pi_{g,d,t}(\cdot)$. Similarly, the identification functional of the ATT estimand then simplifies to \[ \tau_\mathrm{rc}=\ex{}\left[\frac{p(G=1,D=1\mid X)}{\ex{}[G\cdot D]}\cdot\rho_0(X,G,D)\cdot \frac{T-\ex{}[T]}{\ex{}[T](1-\ex{}[T])}\cdot Y\right]. \] Note how the latter resembles the identification functional for the panel data setting ((ref)) after replacing $(Y_1-Y_0)$ by $\frac{T-\ex{}[T]}{\ex{}[T](1-\ex{}[T])}\cdot Y$. This aligns well with analogous results in the canonical DiD literature where no compositional changes assumption is made---see, for instance, abadie2005semiparametric.
In this section, we first present two estimation strategies for each setting based on the identification results of (ref). We discuss the potential issues that these estimation strategies can face. Then we propose efficient and robust estimators based on the influence functions of our parameters of interest to address these issues. Similar to the previous sections, we present our results first for the panel data setting, followed by those for the repeated cross-sections setting.
Based on the definition of the parameter $\psi_{pd}$, one can estimate it using the following plug-in estimator: \[ \hat{\psi}^\mathrm{or}_\mathrm{pd} = \frac{\ex{n}\left[\big(Y_1-Y_0-\hat{\mu}_{0,1,\Delta}(X)-\hat{\mu}_{1,0,\Delta}(X)+\hat{\mu}_{0,0,\Delta}(X)\big)\cdot \ind{G=1,D=1}\right]}{\ex{n}[\ind{G=1,D=1}]}, \] where $\hat{\mu}_{g,d,\Delta}(\cdot)$ is an estimator of the outcome regression function $\mu_{g,d,\Delta}(\cdot)$, and $\ex{n}$ represents empirical mean. Alternatively, one can estimate the parameter of interest based on the identification functional given by (ref). Specifically, letting $\hat{\pi}_{g,d}(\cdot)$ be an estimator of $\pi_{g,d}(\cdot)$, $\psi_\mathrm{pd}$ can be estimated using the following weighting estimator:
where $\hat{e}_\mathrm{pd}$ is an estimate of $1/\ex{}[G\cdot D]$, and
Both $\hat{\psi}_\mathrm{pd}^\mathrm{or}$ and $\hat{\psi}_\mathrm{pd}^\mathrm{ipw}$ provide consistent estimates of the ATT under correct model specification. The OR estimator ($\hat{\psi}_\mathrm{pd}^\mathrm{or}$) relies on correctly modeling the outcome regressions, while the IPW estimator ($\hat{\psi}_\mathrm{pd}^\mathrm{ipw}$) depends on the correct specification of the propensity score models. If either model is misspecified, the respective estimator may be biased. To avoid model misspecifications, one can use nonparametric nuisance estimators. However, the convergence rate of such estimators are often slower than $\sqrt{n}$.
To address the aforementioned issues, we next propose an estimation strategy based on the influence function of $\psi_\mathrm{pd}$. Our estimator is (i) robust to model misspecifications (as described in (ref)), (ii) achieves the semiparametric efficiency bound, and (iii) can achieve parametric convergence rates even if the nuisance function estimators do not converge at $\sqrt{n}$ rate. We start by deriving the efficient influence function of $\psi_\mathrm{pd}$ in the following result.
Based on the influence function in (ref), we propose the following procedure to estimate the parameter $\psi_\mathrm{pd}$. We use the cross-fitting approach chernozhukov2018double to separate the estimation of nuisance parameters from that of the parameter of interest. In particular, we partition the samples into $L$ equally-sized folds of size $m$ indexed by $\{1,\dots,L\}$. For each $\ell\in\{1,\dots,L\}$, let $\hat{\mu}^\ell_{g,d,\Delta}(X)$ and $\hat{\pi}^\ell_{r,g,d}(X)$ be the estimators of $\ex{}[Y_1-Y_0\mid X, G=g, D=d]$ and $\pi_{r,g,d}(X)$, respectively, using the data in all but the $\ell$-th fold. Let $\ee^\ell_m[\cdot]$ represent the empirical mean in the $\ell$-th fold of data. Furthermore, let $\hat{e}_\mathrm{pd}^\ell$ be the estimator of $\frac{1}{\ex{}[G\cdot D]}$ using the data in the $\ell$-th fold, defined as
where we take the maximum in the denominator so that $\hat{e}_\mathrm{pd}^\ell$ is well-defined even if $\ee^{\ell}_m[G\cdot D]=0$. Our estimator for the parameter of interest is:
where
and
Below, we show that our proposed estimator is robust against misspecifications.
Based on (ref), using the estimator $\hat{\psi}^\mathrm{dr}_\mathrm{pd}$ provides the researcher with two opportunities for consistent estimation of the ATT. Specifically, consistency holds as long as for all $\ell,g,d$, either the outcome estimator $\hat{\mu}_{g,d,\Delta}(X)$, or the propensity score ratio estimator $\hat{\pi}_{r,g,d}(X)$ is $L^2(p)$-consistent for the true nuisance function, but not necessarily both.
Next, we show that under certain regularity and convergence rate conditions, our estimator will be $\sqrt{n}$-consistent and asymptotically normal (CAN). Importantly, the requirement is imposed on the product of pairs of convergence rates, rather than on the individual rates themselves. Consequently, none of the nuisance functions are required to converge at the $\sqrt{n}$ rate.
As a corollary of (ref), we can use the influence function $\psi^{1}_\mathrm{pd}(O)$ to obtain confidence intervals for the parameter of interest, $\psi_\mathrm{pd}$. Specifically, for every $\ell\in\{1,\dots,L\}$, we estimate the variance of $\psi^{1}_\mathrm{pd}(\cdot)$ in the $\ell$-th fold as \[ \hat{\sigma}^2_\ell = \ee_m^\ell\left[\Big(\hat{e}^{-\ell}_\mathrm{pd}\cdot\eta_\mathrm{pd}(O; \{\hat{\mu}^\ell_{g,d}\}_{g,d}, \{\hat{\pi}^\ell_{r,g,d}\}_{g,d}) - \hat{e}^{-\ell}_\mathrm{pd}\cdot G\cdot D\cdot \hat{\psi}_\mathrm{pd}^\mathrm{dr}\Big)^2\right], \] where $\hat{e}^{-\ell}_\mathrm{pd}$ is the estimator of $1/\ex{}[G\cdot D]$ using data from all but the $\ell$-th fold of data, defined similarly to (ref). Then we define $\hat{\sigma}^2=\frac{1}{L}\sum_{\ell=1}^L\hat{\sigma}^2_\ell$. Using this estimated variance, the $100(1-\alpha)\%$ confidence interval of $\psi_\mathrm{pd}$ can be obtained as \[ \hat{\psi}_\mathrm{pd}^\mathrm{dr} \pm z_{1-\alpha/2}\frac{\hat{\sigma}}{\sqrt{n}}, \] where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the standard normal distribution.
Similar to the panel data setting, based on the definition of $\psi_\mathrm{rc}$, it can be estimated using the following plug-in estimator:
where $\hat{\mu}_{1,1,0}(\cdot)$ and $\hat{\mu}_{g,d,\Delta}(\cdot)$ are estimators of $\mu_{1,1,0}(\cdot)$, and $(\mu_{g,d,1}-\mu_{g,d,0})(\cdot)$, respectively. Alternatively, (ref) suggests the following weighting-based estimator: \[ \hat{\psi}^\mathrm{ipw}_\mathrm{rc}=\ex{}\big[\hat{e}_\mathrm{rc}\cdot\hat{\pi}_{1,1,1}(X)\cdot\hat{\phi}(X,G,D,T)\cdot Y\big], \] where $\hat{e}_\mathrm{rc}$ is an estimate of $1/\ex{}[G\cdot D\cdot T]$, $\hat{\pi}_{g,d,t}(\cdot)$ is an estimator of the propensity score function $\pi_{g,d,t}(\cdot)$, and \[ \hat{\phi} (X,G,D,T) =-\sum_{g,d,t\in\{0,1\}}\frac{(1-g-G)(1-d-D)(1-t-T)}{\hat{\pi}_{g,d,t}(X)}. \]
As in the panel data setting, these two estimators ($\hat{\psi}^\mathrm{or}_\mathrm{rc}$ and $\hat{\psi}^\mathrm{ipw}_\mathrm{rc}$) can be biased if the corresponding models are misspecified. One can use nonparametric nuisance function estimators to avoid misspecifications, but such estimators often have slow convergence rates. We therefore propose an influence function-based estimator which is (i) robust to misspecifications (see (ref)), (ii) achieves the semiparametric efficiency bound, and (iii) can attain parametric convergence rates even if the nuisance functions converge slower than the $\sqrt{n}$ rate. We begin by deriving the efficient influence function of $\psi_\mathrm{rc}$ in the following result.
The influence function in (ref) suggests the following estimation strategy based on cross-fitting chernozhukov2018double. We partition the data into $L$ folds of size $m$, indexed by $\{1,\dots,L\}$. For each $\ell\in\{1,\dots,L\}$, we let $\hat{\mu}^\ell_{g,d,t}(X)$ and $\hat{\pi}^\ell_{r,g,d,t}(X)$ be estimators of $\ex{}[Y\mid X, G=g, D=d, T=t]$ and $\pi_{r,g,d,t}(X)$, respectively, using the data in all but $\ell$-th fold. Additionally, let $\hat{e}^\ell_\mathrm{rc}$ be the estimator of $(1/\ex{}[G\cdot D\cdot T])$ using the data in the $\ell$-th fold, defined as
Our proposed estimator for $\psi_\mathrm{rc}$ is: \[ \hat{\psi}_\mathrm{rc,1}^\mathrm{dr}= \frac{1}{L}\sum_{\ell=1}^L\hat{e}^\ell_\mathrm{rc}\cdot\ee_m^\ell \Big[ \eta_\mathrm{rc,1}(O; \{\hat{\mu}^\ell_{g,d,t}\}_{g,d,t}, \{\hat{\pi}^\ell_{r,g,d,t}\}_{g,d,t}) \Big], \] where \[ \eta_\mathrm{rc,1}(O;\{\hat{\mu}^\ell_{g,d,t}\}_{g,d,t}, \{\hat{\pi}^\ell_{r,g,d,t}\}_{g,d,t}) \coloneqq \sum_{g,d,t\in\{0,1\}}(-1)^{(g+d+t)}\hat{\omega}^\ell_{g,d,t}(X,G,D,T)\cdot\big(Y-\hat{\mu}^\ell_{g,d,t}(X)\big), \] and, \[ \hat{\omega}_{g,d,t}(X,G,D,T) = G\cdot D\cdot T - \hat{\pi}_{r,g,d,t}(X)\cdot\mathbbm{1}\{G=g,D=d,T=t\}. \] Below, we show that $\hat{\psi}_\mathrm{rc,1}^\mathrm{dr}$ is robust against misspecifications.
(ref) demonstrates that $\hat{\psi}^\mathrm{dr}_\mathrm{rc,1}$ is consistent as long as for all $\ell,g,d,t$, either the outcome regression estimator $\hat{\mu}_{g,d,t}$ or the propensity score ratio estimator $\hat{\pi}_{r,g,d,t}$ (but not necessarily both) is $L^2(p)$-consistent for the true nuisance function. Next, we show that under certain regularity and convergence rate conditions, our estimator will be $\sqrt{n}$-consistent and asymptotically normal (CAN). Similarly to the panel data setting, the requirement is imposed on the product of pairs of convergence rates, rather than on the individual rates themselves. Consequently, none of the nuisance functions are required to converge at the $\sqrt{n}$ rate.
As a corollary of (ref), one can use the influence function $\psi^1_\mathrm{rc,1}(\cdot)$ to obtain confidence intervals for $\psi_\mathrm{rc}$. Specifically, for every $\ell\in\{1,\dots,L\}$, we estimate the variance of $\psi^1_\mathrm{rc,1}(O)$ in the $\ell$-th fold as
where $\hat{e}^{-\ell}_\mathrm{rc}$ is the estimator of $1/\ex{}[G\cdot D\cdot T]$ using data in all but $\ell$-th fold, defined similarly to (ref). Then we define the estimated variance $\hat{\sigma}^2=\frac{1}{L}\sum_{\ell=1}^L\hat{\sigma}^2_\ell$, and, the $100(1-\alpha)\%$ confidence interval of $\psi_\mathrm{rc}$ can be obtained as \[ \hat{\psi}_\mathrm{rc,1}^\mathrm{dr} \pm z_{1-\alpha/2}\frac{\hat{\sigma}}{\sqrt{n}}, \] where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of the normal distribution.
Both (ref) and (ref) allow for compositional changes. In (ref), we showed how (ref) simplifies weighting-based identification. Here, we show that this assumption also has significant implications for estimation. To this end, we first provide the efficient influence function of $\psi_\mathrm{rc}$ under (ref).
Let $\hat{\mu}^\ell_{g,d,t}(X)$, and $\hat{\pi}^\ell_{r,g,d}(X)$ be estimators of $\ex{}[Y\mid X, G=g, D=d, T=t]$, and $\pi_{r,g,d}(X)$ respectively, using the data in all but $\ell$-th fold. Also, define \[ \hat{e}^\ell_\mathrm{rc,2}=\frac{1}{\max\{\frac{1}{m},\ee_m^{\ell}[G\cdot D]\}}. \] Here, based on the influence function in (ref), we propose the following estimator for $\psi_\mathrm{rc}$ under (ref): \[
\] where \[
\] and
The following result demonstrates the robustness of $\hat{\psi}_\mathrm{rc,2}^\mathrm{dr}$ under (ref) against misspecifications.
As evident from (ref), under (ref), a stronger robustness can be achieved. In particular, $\hat{\psi}_{\mathrm{rc},2}^\mathrm{dr}$ is consistent under similar assumptions to the panel data setting; it suffices to have access to consistent estimators of either $\mu_{g,d,1}-\mu_{g,d,0}$, or $\pi_{r,g,d}$. However, without this assumption, one needs consistent estimators of both $\mu_{g,d,1}$ and $\mu_{g,d,0}$ (as opposed to only their contrast), or consistent estimators of both $\pi_{r,g,d,1}$ and $\pi_{r,g,d,0}$ (as opposed to only $\pi_{r,g,d}$). Under (ref), the required assumptions for $\sqrt{n}$-consistency and asymptotic normality of $\hat{\psi}^\mathrm{dr}_\mathrm{rc,2}$ are weaker than those for $\hat{\psi}^\mathrm{dr}_\mathrm{rc,1}$. In particular, convergence rate conditions are imposed only on $(\hat{\mu}_{g,d,1}-\hat{\mu}_{g,d,0})$ and $\hat{\pi}_{r,g,d}$, rather than on $\hat{\mu}_{g,d,1},\hat{\mu}_{g,d,0},\hat{\pi}_{r,g,d,1}$, and $\hat{\pi}_{r,g,d,0}$ (which was required in (ref)).
Once again, a corollary of (ref) is that confidence intervals for the parameter of interest can be obtained by estimating the variance of the influence function $\psi^1_\mathrm{rc,2}(O)$.
In this section, we present simulation studies to assess the performance of our proposed methodology. We adapted the data-generating mechanism of kang2007demystifying to the triple-difference setting, with the key modification of introducing an unobserved confounder $U$, drawn from a standard normal distribution. Four observed covariates $(X_1, X_2, X_3, X_4)$ were constructed from nonlinear transformations of latent normal variables, while treatment eligibility and domain assignments $(G,D)$ were jointly generated from a multinomial logit model that depends on both observed covariates and $U$. Outcomes were generated from nonlinear functions of $(X_1,\dots,X_4)$ and $U$, with potential outcomes defined to satisfy the conditional parallel difference in trends assumption ((ref) and (ref) in the panel data and repeated cross-sections settings, respectively). In the repeated cross-sections setting, the time indicator $T$ was sampled from an expit model with a parameter that depends on both the covariates and the treatment and domain indicators, thereby violating no-compositional changes ((ref)). A detailed description of our data generating process is included in (ref).
In our evaluations, we considered four different variants of estimators $\hat{\psi}^\mathrm{dr}_\mathrm{pd}$ and $\hat{\psi}^\mathrm{dr}_\mathrm{rc,1}$: (i) with both outcome regression (OR) functions and propensity scores (PS) correctly specified, (ii) with correctly specified OR functions but misspecified PS, (iii) with misspecified OR functions but correctly specified PS, and (iv) with all nuisance functions misspecified. To help ensure that the nuisance functions were adequately captured and the estimator can be regarded as correctly specified, we used fully connected neural networks, designed with sufficient depth and width to flexibly model complex, non-linear nuisance functions. For OR functions, we used networks with three hidden layers, while for PS we used networks with two hidden layers. For the case of misspecified OR function, we used a ridge regression estimator that takes only $X_1$, the first observed covariate, into account. For the case of misspecified PS, we used a logistic regression estimator that works only with $X_1$, the first observed covariate. In all four cases, cross-fitting with three folds was employed. The complete implementation details and reproducing python code are available at \href{https://github.com/SinaAkbarii/triplediff}{https://github.com/SinaAkbarii/triplediff}.
(ref) illustrates the relative bias of each of the variants of $\hat{\psi}_\mathrm{pd}^\mathrm{dr}$ in the panel data setting with sample sizes ranging from $1000$ to $10000$. Relative bias is defined as $(\hat{\psi}/\psi)-1$, where $\hat{\psi}$ is the estimate and $\psi$ is the true value of the parameter of interest, which was $10$ in our data generating mechanism. Each boxplot summarizes results from 1,000 Monte Carlo replications. The results show that when at least one set of the OR functions or PS is correctly specified, the relative bias converges to zero, confirming the double-robustness property. In contrast, when all nuisance functions are misspecified, the estimates converge to a value that differs from the true parameter, causing the relative bias to converge to a negative value rather than vanish. (ref) presents the numerical values of average bias and root mean squared error (RMSE).
(ref) and (ref) show the relative bias, average bias, and RMSE of variants of $\hat{\psi}_\mathrm{rc,1}^\mathrm{dr}$ estimators in the repeated cross-sections setting. Overall, the patterns closely mirror those observed in the panel data setting. However, for a fixed sample size, these estimators exhibit greater variance than their panel data counterparts. This is expected, since each panel data sample contains both $Y_0$ and $Y_1$ and therefore provides more information, whereas each repeated cross-sections sample includes only one outcome. To make the comparison clearer, in (ref), we report the bias, RMSE, and empirical coverage of 95% confidence intervals of our estimators $\hat{\psi}_\mathrm{pd}^\mathrm{dr}$ and $\hat{\psi}_\mathrm{rc,1}^\mathrm{dr}$, when the nuisance functions are correctly specified. Moreover, the results of this table support the asymptotic normality of these estimators, as coverage exceeds $95\%$ with increasing sample size.
We applied our method to examine the impact of mandated maternity benefits on wages of women of childbearing age in the United States. We used the May Current Population Survey (CPS) of census_cps_may_2023 for the years 1974--1975, representing the pre-treatment period, and 1977--1978, representing the post-treatment period. The treatment of interest ($A$) is the adoption of state-level laws in the late 1970s requiring employers to offer health insurance policies covering childbirth costs and maternity benefits. The outcome of interest ($Y$) is the logarithm of the hourly wages. All wage measures were adjusted using the Consumer Price Index (CPI) for the corresponding year: 49.3 (1974), 53.8 (1975), 60.6 (1977), and 65.2 (1978). The design of our study follows that of gruber1994maternity. The analysis was restricted to eight states: Connecticut, Illinois, Indiana, Massachusetts, New Jersey, New York, North Carolina, and Ohio, where Illinois, New Jersey and New York adopted the policy and were considered as the target domain ($D=1$), whereas the other states did not do so and served as the reference domain ($D=0$). Eligibility for the maternity mandate is defined by an indicator variable $G=1$ for married women aged 20--40 (childbearing age) and zero otherwise. We used the following demographic and occupational covariates ($X$): education, age, sex, race (white/non-white), marital status, union status, and white-collar occupation status. We further restricted the sample to individuals aged 20--65 and followed Gruber's exclusion criteria, removing single or divorced women aged 20--40 as well as married men in the same age range. Individuals with hourly wages below $1$ or above $100$ US dollars were also excluded.
The primary distinction between our analysis and that of gruber1994maternity lies in the estimation strategy. While Gruber considered a simple parametric model with fixed effects for the outcome such that one of the parameters will become the causal parameter of interest, we used our semiparametric influence function-based estimator $\hat{\psi}_\mathrm{rc,1}^\mathrm{dr}$, where nuisance functions were estimated using fully connected neural networks as described in (ref). The results of our analysis is presented in (ref). For the standard error and p-value (for testing $H_0$: effect $=0$), we report results based on both the influence function and bootstrap through 1000 replications.\footnote{The python implementations to reproduce these results can be found at \url{https://github.com/SinaAkbarii/triplediff}.} We note that mandated maternity benefits increase labor costs for employers, which are likely passed on to employees in the form of lower wages. In particular, it is implausible that mandated maternity benefits would raise wages for married women of childbearing age. For this reason, we report a one-sided p-value. Our analysis concluded a point estimate of $\hat{\psi}_\mathrm{rc,1}^\mathrm{dr} = -0.02633$ with standard error $0.02223$ for change in the log of hourly wages. This implies that the mandates resulted in a drop of $2.6\%$ in the hourly wages. This is slightly lower than Gruber's conclusion, which was a $4.2\%$ drop in hourly wages. For evaluating the significance, the one-sided p-values are $0.118$ and $0.109$, obtained from the influence function and bootstrap approaches, respectively. Taken together, these results provide evidence that mandated maternity benefits reduced the wages of married women of childbearing age in adopting states.
We studied the identification and estimation of the average treatment effect on the treated within the triple difference framework, focusing on both panel data and repeated cross-sections settings. From the identification standpoint, we presented the first weighting-based identification results, notably, while allowing for compositional changes in the repeated cross-sections setting. From an estimation standpoint, we proposed semiparametric estimators for the triple difference framework in both panel data and repeated cross-section settings. These estimators employ a cross-fitting approach, allowing flexible machine learning methods to estimate the nuisance functions. We identified the conditions under which our estimators are efficient, doubly robust, $\sqrt{n}$-consistent, and asymptotically normal. Our work lays a foundation for further analysis of the triple difference framework. Given the practical relevance and growing use of this framework in empirical research, there is ample scope for future studies to extend such analyses, especially by exploring settings with multiple time periods, dynamic treatment effects, and continuous treatments. As an application of our proposed methodology, we assessed the effect of mandated maternity benefits on the hourly wages of women of childbearing age and found that these mandates result in a $2.6\%$ drop in hourly wages. Our estimated 2.6% reduction in hourly wages is somewhat smaller than the 4.2% reported by gruber1994maternity; however, while Gruber relied on a simple parametric model, our approach employs nonparametric machine learning methods, which are more flexible and reliable. Therefore, the effect of mandated maternity benefit policies might be slightly smaller than previously estimated.