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.
97,870 characters · 21 sections · 97 citation commands
Semiparametric Bayesian Difference-in-Differences
{ Keywords: Difference-in-differences, conditional parallel trends, semiparametric Bayesian inference, Bernstein–von Mises theorem, double robustness. \ }
The Difference-in-Differences (DiD) method is widely used in causal inference. It is particularly effective for evaluating policy interventions while accounting for unobserved time-invariant heterogeneity. The primary parameter of interest in this context is the average treatment effect on the treated (ATT). One of its key identifying conditions is the (conditional) parallel trends assumption, i.e., treated and control groups would exhibit similar trends absent treatment after adjusting for covariates abadie2005,sant2020doubly. While the related literature is largely frequentist, this paper introduces a Bayesian framework under conditional parallel trends, avoiding parametric assumptions on model primitives. Our approach yields point estimates and credible sets in a unified manner.
We propose two novel Bayesian methods for inference on the ATT in the DiD framework. First, we propose the Bayesian procedure using Gaussian process priors. This method places the Gaussian process prior on the conditional mean function for the control group and a Dirichlet process prior on the remaining part of the distribution in the likelihood. This avoids the need to impose Bayesian modeling on either the conditional mean of the treated group or the propensity score. Our method can be viewed as the Bayesian counterpart of heckman1997matching, with the added advantage of enabling automatic uncertainty quantification through the posterior distribution. We show that this Bayesian method satisfies the Bernstein-von Mises (BvM) theorem under regularity conditions and is therefore asymptotically equivalent to semiparametric efficient frequentist estimators. While its asymptotic BvM property does not hold under double robust smoothness conditions, the Bayesian method performs well empirically when the number of continuous covariates is moderate and in scenarios where the overlap assumption is nearly violated. This robustness stems, in part, from the Gaussian process prior being specified solely on the conditional mean function for the control group.
We also provide an extension of our Bayesian procedure, which incorporates robustification via estimated propensity scores and is particularly suited for more complex models, either due to a larger number of continuous covariates or when the underlying conditional mean functions are not smooth. Our Double Robust Bayesian procedure adjusts the prior and posterior distributions by incorporating an efficient influence function. By doing so, we leverage the rich frequentist literature on double-robust estimation, specifically sant2020doubly in the DiD framework, without sacrificing many of the desirable properties of the Bayesian approach. Under double-robust smoothness conditions, our robust Bayesian procedure satisfies the semiparametric Bernstein–von Mises (BvM) theorem, albeit with a “bias term" in the posterior. Specifically, the resulting posterior distribution depends on the unknown true conditional mean and propensity score functions. Our double-robust Bayesian approach addresses this “bias term" by incorporating an explicit posterior correction. Both the prior adjustment and the posterior correction are derived from functional forms closely associated with the efficient influence function of estimating the ATT.
In our Monte Carlo simulations, we find that our methods result in improved empirical coverage probabilities while maintaining competitive confidence interval lengths compared to existing frequentist methods. This finite sample advantage is also observed in low dimensional cases for our Bayesian method that does not involve prior or posterior corrections.\footnote{In contrast, a Bayesian method without prior correction performs poorly for the average treatment effect (ATE) as shown by BLY2022.} This can be explained by the construction of the Bayesian procedure for the ATT, which involves only a prior specification for the conditional mean function in the control arm. In particular, we note that our approach leads to more accurate uncertainty quantification and is less sensitive to estimated propensity scores that are close to boundary values. Our Bayesian methodology requires a prior specification through a likelihood function for the control arm, for which we impose an exponential family structure. In the Gaussian case, for instance, this leads to a computationally efficient procedure with the posterior being multivariate Gaussian, which avoids computational demanding methods like MCMC. We stress that the misspecification of this structure does not have serious consequences for estimation of the ATT. First, as shown in kleijn2006misInf, nonparametric Bayesian methods possess the same robustness to misspecification as the frequentist M-estimation using least squares. Second, we provide finite sample evidence through simulations, where we find that our Bayesian procedures are not sensitive to misspecifications of the likelihood functions.\\
In the related literature, nonparametric Bayesian causal inference has recently received considerable interest; see, for example, the numerous applications in daniels2024BNP. ray2020causal develop the comprehensive theory for establishing the BvM theorem in the missing data framework, employing Gaussian process priors for the conditional mean function. Extending their methodology to the ATT would require nonparametric Bayesian modeling of both the propensity score and the conditional mean function for the treated group, as discussed in Remark (ref). A key innovation of our proposal is to circumvent this route by building our Bayesian procedure on a reparametrization that is particularly convenient for the analysis of the ATT.
ray2020causal also propose the novel prior adjustment to the conditional mean, which makes use of the estimated propensity score. Building on this prior adjustment, BLY2022 introduce a debiasing step to further correct the posterior and establish the BvM theorem for the average treatment effect (ATE) under double robustness. Although they outline the extension to general semiparametric models where the parameter of interest can be written as the linear functional of conditional means, this approach does not cover the case of the ATT, because of its ratio form. To address the random denominator that estimates the proportion of treated individuals, we apply a new conditional Slutsky lemma introduced by yiu2023corrected in the Bayesian context. Also, the correction steps in our Bayesian method share the same motivation as in BLY2022, but differ in their functional form. This is in line with the well known subtle differences between the cases for ATE and ATT; see hahn1998role. Additionally, we find that the bias term in the BvM for ATT is substantially simpler than that for ATE, which also explains the favorable finite-sample behavior of Bayesian ATT estimators even without posterior corrections. In contrast, yiu2023corrected suggest a different type of posterior correction that assumes stronger regularity conditions in the context of the ATT. To the best of our knowledge, our proposed double robust BvM theorem for the ATT is the first to relax the Donsker property assumption for the conditional mean (see also Remark (ref) for a detailed comparison).
Our paper is also connected to the broader literature on robustifying standard Bayesian procedures in econometrics. Regarding Bayesian inference methods for partially or weakly identified models, we refer readers to chen2018MC,giacomini2020robust,andrews2022gmm. Under local misspecification of parametric models, muller2024 establish a novel BvM result utilizing the efficient influence function. There are also scattered results exploring Bayesian methodology to the study of ATT or DiD. In an earlier paper, chib2002treat developed a semiparametric Bayesian model for the ATT in both cross-sectional and panel data settings. Their semiparametric model differs from our setup in that covariates enter the outcome equation linearly, while the error term is modeled using flexible Dirichlet process mixtures. Recently, in the context of assessing sensitivity to the parallel trends assumption, kwon2024empirical proposed a Bayesian approach.\\
The remainder of this paper is organized as follows. Section (ref) presents the setup and introduces the Bayesian framework in the DiD setup. Section (ref) outlines our Bayesian methods. In Section (ref), we establish inference via semiparametric BvM theorems for our first method. In Section (ref), we derive a doubly robust, semiparametric BvM theorem for our second method. Section (ref) provides BvM results under primitive conditions when using squared expontential process priors. Section (ref) presents finite sample results via simulations and an empirical illustration. Towards the end, we outline two extensions of our methodology: Section (ref) provides an extension to the staggered intervention with multiple time periods and repeated cross-sectional data. Proofs of main theoretical results are collected in Appendix (ref). Supplementary Appendices (ref)--(ref) provide additional technical results and further simulation evidence.
This section provides the main setup of the average treatment effect on the treated (ATT) in the difference-in-differences (DiD) design. We first provide standard conditions for the identification of the ATT and introduce additional notations under the Bayesian formulation of the problem.
We focus on the canonical DiD design case, where there are two treatment periods and two treatment groups. Let $Y_{i t}$ be the outcome of interest for unit $i$ at time $t$. We assume that researchers have access to outcome data in a pre-treatment period $t=1$ and in a post-treatment period $t=2$. Let $D_{i t}=1$ if unit $i$ is treated before time $t$ and $D_{i t}=0$ otherwise. Note that $D_{i 1}=0$ for every $i$ and thus we may write $D_i=D_{i 2}$. Using the potential outcome notation, $Y_{it}(0)$ or $Y_{it}(1)$ denotes the outcome of unit $i$ at time $t$ if it does not receive or receives treatment by time $t$, respectively. Thus, the realized outcome for unit $i$ at time $t=1$ is $Y_{i 1}= Y_{i 1}(0)$, and at time $t=2$ it is $Y_{i 2}=D_i Y_{i 2}(1)+\left(1-D_i\right) Y_{i 2}(0)$. Below, $P_0$ denotes the frequentist distribution generating the observed data.
A vector of $p$-dimensional pre-treatment covariates $X_i$ is also available, with cumulative distribution function denoted by $F_X$.\footnote{If $X_i$ does not have a density we can simply consider the conditional density of $(\Delta Y_i, D_i)$ given $X_i=x$ instead of the joint density of $(\Delta Y_i, D_i, X_i)$.} Let $\pi_0(x)=P_0(D_i=1\mid X_i=x)$ denote the propensity score, $\pi_0=P_0(D_i=1)$ the proportion, and $m_0(x)= \mathbb E_0[\Delta Y_i\mid D_i=0, X_i=x]$ the conditional mean of the differenced outcome across two periods, where $\Delta Y_{i}:=Y_{i2}-Y_{i1}$ and where $\mathbb{E}_0[\cdot]$ denotes the expectation under $P_0$. The researcher observes an independent and identically distributed (i.i.d.) observations of $(Y_{i1},Y_{i2},D_i,X_i^\top)^\top$, $i=1,\dots,n$. In addition to this canonical panel data setup, we discuss how our results translate to repeated cross-sections data in Section (ref), while an extension to staggered intervention is provided in Section (ref). In addition to this canonical panel data setup, we provide an extension to staggered interventions in Section (ref) and discuss how our results translate to repeated cross-sectional data in Section (ref). For notational simplicity, we will henceforth suppress the unit index $i$.\\
Regarding the causal effect in the canonical DiD setup, the related literature primarily focuses on the average treatment effect on the treated (ATT) given by
For its identification, we impose the no anticipation assumption, conditional parallel trends (PTA) given covariates $X$, and the weak overlap conditions as follows.
Under Assumption (ref) the ATT is identified by
One can construct an estimator that replaces the conditional mean function $m_0$ with an estimator, known as the outcome regression approach, as described in heckman1997matching. As noted by abadie2005, plug-in estimators based on standard nonparametric estimators of the conditional mean function $m_0$ can face significant challenges due to the curse of dimensionality. This is where we can capitalize the strength of Bayesian estimation, which allows us to incorporate rich covariate information in the prior distribution. It has also been noted in the recent literature that, in the presence of heterogeneous treatment effects in $X$, i.e., when $\mathbb E_0[Y_{2}(1)-Y_{2}(0)\mid X=x, D=1]$ varies with $x$, the two-way fixed effect estimator (TWFE) is in general not consistent for the ATT. See also Remark 1 of sant2020doubly for an explicit discussion.
We now provide the formal Bayesian setup to the ATT in the DiD context. We consider a family of probability distributions $\{P_\eta:\eta\in\mathcal H\}$ for some parameter space $\mathcal H$. The (possibly infinite dimensional) parameter $\eta$ characterizes the probability model. Let $\eta_0$ be the true value of the parameter and denote $P_0=P_{\eta_0}$, which corresponds to the frequentist distribution generating the observed data. Under $P_\eta$ where $\eta=(\pi,f_{X}, f_{\Delta Y|D,X})$, the joint density function of $Z=(\Delta Y,D, X^\top)^\top$ can thus be written as
where $f_{\Delta Y \mid D,X}(y \mid d,x)$ for $d = 1$ is the unrestricted conditional density of $\Delta Y$ given $(D, X)$, while for $d = 0$, we impose the exponential family condition as in (ref). Here, $f$ denotes the joint density of $(D\Delta Y, D, X^\top)^\top$ under $P_\eta$, and the corresponding cumulative distribution function is denoted by $F$. Importantly, specifying only the prior distribution on density function $f$ and the conditional mean function $m$ is sufficient for identifying the ATT parameter of interest. Specifically, there is no need to additionally parameterize the propensity score $\pi$, the marginal density of $X$, or the conditional density for the control group, $f_{\Delta Y| D, X}(y\mid 0,x)$, which is a deterministic function of $m(x)$ due to the exponential family assumption imposed in (ref) below.
We consider the following reparametrization of $(m, f)$ given by $\eta=(\eta^m,\eta^{f})$. A central insight of this paper is to show that a nonparametric process prior specification on the conditional mean function $m$ and the density $f$ is sufficient for the Bayesian inference on the ATT parameter. We index the probability model by $P_{\eta}$, where
for some known, invertible function $q(\cdot)$, which we specify below. We can write the ATT depending on a hyperparameter $\eta$ as
where $\mathbb{E}_\eta$ denotes the expectation under $P_\eta$ and in this case, is the integral with respect to the density $f_\eta$.
As we saw above, we only need to impute the conditional mean of the outcome in the control group, making it unnecessary to impose a model for the treated group. We assume that the distribution of $\Delta Y$, conditional on $D=0$ and $X$, belongs to the “single-parameter" exponential family, where the unknown parameter is the nonparametric conditional mean function $m(x)=\mathbb{E}[\Delta Y\mid D=0,X=x]$. Specifically, we assume that the conditional density function is given by
where $A(m)= \log\int c(y)\exp\left[q(m) ay\right]dy$, some constant $a>0$, and the function $q(\cdot)$ links the conditional mean to the “natural parameter” of the exponential family. We also restrict the sufficient statistic to be linear in $y$. The exponential family assumption implies the conditional mean equation $\mathbb{E}[ \Delta Y\mid D=0, X=x]=A'(m(x))/(a\, q'(m(x)))$, which corresponds to generalized regression models. Interestingly, in the exponential family examples provided below, the previous equation implies $\mathbb{E}[\Delta Y\mid D=0, X=x]=m(x)$ and hence the exponential family assumption (ref) does not impose functional form assumptions on the conditional mean function $m$ in these cases and, in particular, does not restrict the ATT parameter.
The family ((ref)) allows for counting and continuous outcomes. For instance, when $a=1$, the Poisson distribution corresponds to the choices $c(y)= 1/(y!)$, $q(m)=\log m$, and $A(m)=m$, while the exponential distribution is represented by $c(y)=1$, $q(m)=-1/m$, and $A(m)=\log m$. Furthermore, the normal distribution with $\text{Var}(\Delta Y|D=0,X)=\sigma^2$ for some $\sigma>0$, is captured by $c(y)=\exp(-y^2/(2\sigma^2))/\sqrt{2\pi\sigma^2}$, $q(m)=m/\sigma$, $A(m)=m^2 / (2\sigma^2)$, and $a=1/\sigma$. For the normal case, we treat $\sigma$ as a hyperparameter and estimate it together with hyperparameters in Gaussian process prior by maximizing marginal likelihood. We note that while a generalization to multinomial outcomes, as in BLY2022, is possible, we do not consider this case explicitly in this paper.
A high-level assumption in our BvM theorem requires the posterior contraction of $m_{\eta}$ to the true conditional mean $m_0$. There are cases that this holds even if the exponential family is misspecified. Generally, the posterior of $m_{\eta}$ will contract on near the point (pseudo-true value) in the support of the prior that minimize the Kullback–Leibler (KL) divergence with respect to the true data generating probability. Related posterior contraction results and cases where the pseudo-truth concides with $m_0$ can be found in kleijn2006misInf. This aligns with our finite sample results, which are not sensitive to deviations from exponential family distributions, as shown in Appendix (ref). Beyond the exponential family, one can also consider the flexiable nonparametric Bayesian approach in norets2022 to model the conditional distribution of the outcome for the control group.
We now present two Bayesian procedures that build on flexible prior processes, enabling semiparametric inference on the parameter of interest. The first corresponds to a nonparametric Bayesian approach based on standard Gaussian process priors. The second involves Bayesian methods with frequentist modifications, incorporating an adjustment to the prior along with a posterior correction.
We first consider nonparametric Bayesian inference, which builds on a standard Gaussian process prior for the conditional mean function combined with an independent Dirichlet process prior for the conditional expectation in (ref). The proposed method does not include a propensity score adjustment, which prevents it from achieving double robustness, as we see in the next section. In our simulation results, however, we find that the proposed method is robust even in cases of near overlap failure.
The use of Gaussian process priors for the conditional mean has the following motivation. The mode of a posterior stemming from Gaussian process priors can be derived by a minimization problem involving the corresponding norm of a so-called reproducing kernel Hilbert space (RKHS). Gaussian process (GP) priors share close ties with spline estimation wahba1990spline, a connection that—along with their strong finite-sample performance—has fueled their popularity in machine learning rassmusen2006gaussian,murphy2023pml. For other notable applications in econometrics, see kasy2018tax, chib2018moment and florens2019gaussian.
The Dirichlet process is default prior on spaces of probability measures. By the definition of the ATT $\tau_{\eta}$, we assign a Dirichlet process prior to model the distribution $F_{\eta}$, which induces the so-called Bayesian bootstrap when the base measure of the Dirichlet process is taken to be zero; see rubin1981bayesian and also chamberlain2003bayesian.
The Bayesian Algorithm (ref) allows for simultaneous point estimation and uncertainty quantification. Our $100\cdot(1-\alpha)\%$ credible set $\mathcal{C}_n(\alpha)$ for the ATT parameter $\tau_0$ is computed by
where $q(a)$ denotes the $a$ quantile of $\{\tau_\eta^s:s=1,\ldots,B\}$. We also obtain the Bayesian point estimator (the posterior mean) by averaging the simulation draws: $\overline{\tau}_{\eta}=B^{-1}\sum_{s=1}^B \tau_\eta^s$.
For the choice of the prior process $W^m$, we use a Gaussian process with mean $\mu$ and the squared exponential (SE) covariance function $K\left(\cdot,\cdot\right)$ rassmusen2006gaussian given by
where the hyperparameter $\nu^2$ is the kernel variance and $a_{1n},\ldots,a_{pn}$ are rescaling parameters that reflect the relevance of each covariate in predicting $\eta^m$. In practice, the hyperparameters $\mu$, $\nu$, and $a_{1n},\ldots,a_{pn}$ can be chosen by maximizing the marginal likelihood. When the exponential family specification in ((ref)) takes the Gaussian form, Step (a) of posterior computation in Algorithm (ref) is analytically tractable and computationally very efficient, see Supplementary Appendix (ref) for details. For non-Gaussian cases, one can use Laplace approximation or Monte Carlo sampling for Step (a).
Our Bayesian approach relies on prior correction via inverse propensity score weighting (IPW) in the least favorable direction, as specified by the efficient influence function. In contrast to abadie2005, we do not incorporate IPW directly. Importantly, we make use of IPW for the prior and posterior correction of our Bayesian procedure. This resembles sant2020doubly who combine OR and IPW to achieve doubly robust estimation in the frequentist setting. Following hahn1998role,hirano2003efficient or, in the DiD setup sant2020doubly, the efficient influence function for the ATT is given by
with its Riesz representer $\gamma_\eta$ given by
We show in the Supplemental Appendix (ref) that the Riesz representer $\gamma_\eta$ determines the least favorable direction associated with the Bayesian submodel with the largest variance. Our prior adjustment using this Riesz representer provides exact invariance under shifts in nonparametric components along this direction. This extends the work of ray2020causal on unconditional average treatment effects to the ATT case, where the least favorable function, given by the efficient influence function, takes a different functional form. In a similar vein to BLY2022, we use the Riesz representer to correct for posterior bias under double robust smoothness conditions.
Our prior and posterior adjustments depend on a preliminary estimator of $\gamma_0$. A pilot estimator for the propensity score $\pi_0(\cdot)$ is denoted by $\widehat{\pi}(\cdot)$, based on an auxiliary sample, as is the estimator of the treated proportion $\pi_0$, which is taken to be the sample mean of the treatment indicators. We also make use of a pilot estimator $\widehat{m}$ for the conditional mean function $m_0$. We consider a plug-in estimator for the Riesz representer $\gamma_0$ given by
One could also consider an estimation proportion based on the sample information of $D_i$ alone, yet this would complicate the theoretical analysis without improving the finite sample performance of our procedure. The use of an auxiliary data for the estimation of unknown functional parameters simplifies the technical analysis and is common in the related Bayesian literature; see ray2020causal for propensity score adjusted priors in the case of missing data. In practice, we use the full data twice and do not split the sample, as we have not observed any over-fitting or loss of coverage thereby.
Algorithm (ref) describes our double robust Bayesian procedure that approximates the posterior distribution of $\tau_{\eta}$ given in equation (ref). Let $n_c$ denote the number of observations in the control arm. Based on simulations, we recommend the following choices of the pilot estimators. The initial estimator $\widehat \gamma$, given in ((ref)), is implemented based on logistic regression for the propensity scores $\pi(x)$ and the sample average of the treated proportion $\pi$. The pilot estimator for $m(x)$ is implemented by $\widehat m(x)=\sum_{s=1}^B m_\eta^{s}(x)/B$, where $m_\eta^{s}$ is obtained in Step (a) of the posterior computation in Algorithm (ref), and $B$ denotes the number of posterior draws.
Algorithm (ref) also leads to simultaneous point estimation and uncertainty quantification. The $100\cdot(1-\alpha)\%$ credible set $\mathcal{C}_n(\alpha)^{DR}$ for the ATT parameter $\tau_0$ is as in (ref), but here $q(a)$ denotes the $a$ quantile of $\{ \check{\tau}_\eta^s:s=1,\ldots,B\}$. The Bayesian point estimator by $\overline{\tau}_{\eta}^{DR}=B^{-1}\sum_{s=1}^B \check{\tau}_\eta^{s}$.
In this section, we establish a Bernstein-von Mises Theorem using standard Gaussian process priors as considered in our Bayesian procedure in Algorithm (ref).
We now provide additional notations used for the derivation of our semiparametric Bernstein-von Mises Theorem. Recall that we restrict the joint density for the control arm only, imposing the exponential family restriction as in (ref). We denote the observed data corresponding to the treated part as $Z_{\text{Treat}}^{(n)}:=(X_i,D_i,D_i\Delta Y_i)$. We express the posterior as follows:
where the conditional density $f_{\Delta Y|D,X}$ is a function of the conditional mean $m$ by the exponential family restriction given in (ref). Here, we used the fact that independent priors are placed on the conditional mean $m$ and the distribution function $F$.
We first introduce assumptions, which are high-level, and discuss primitive conditions for those in the next section. Below, we consider some measurable sets $\mathcal H^m_n$ of functions $\eta^m$ is understood only for the control arm such that $\Pi(\eta^m\in\mathcal{H}^m_n\mid Z^{(n)})\to_{P_0} 1$. To abuse the notation for convenience, we also denote $\mathcal{H}_n=\{\eta:\eta^m\in\mathcal{H}_n^m\}$ when we index the conditional mean function $m_{\eta}$ by its subscript $\eta$. We write the expression $\|\phi\|_{2,F_0}:= \sqrt{\int \phi^2(z)\mathrm{d}F_0(z)}$ for all $\phi\in L^2(F_0):=\{\phi:\|\phi\|_{2,F_0}<\infty\}$. When we consider the conditional moment function $m$ below, the integral simplifies to one that depends only on the marginal distribution of $X$ under $P_0$.
The posterior contraction rate for the conditional mean can be derived by modifying the classical results of ghosal2000rates. In the related literature, the requirement $\varepsilon_n=o(n^{-1/4})$ is stated explicitly in order to eliminate second-order remainder terms; see Condition (C) in castillo2012gaussian. This also aligns with the usual cut-off rate of the nonparametric components in frequentist semiparametric models newey1994var. Note that the ATT $\tau_{\eta}$ is linear in $m_{\eta}$, so that we do not need to deal with these second-order terms. Nevertheless, the posterior contraction rate also plays a crucial role in the next two assumptions related to the stochastic equicontinuity and prior stability. For the concrete example involving the H\"older class for the conditional mean function, we need to impose sufficient smoothness so that this contraction rate indeed satisfies $\sqrt n\varepsilon^2_n=o(1)$. If the exponential family structure (ref) is misspecified, the posterior contracts to the point in the support of the prior that is closest to the true distribution (as measured by the Kullback-Leibler divergence). Specifically, if one starts with the Gaussian model, the contraction result required by Assumption (ref) can be established utilizing Theorem 4.1 of kleijn2006misInf, provided that the true conditional mean lies within the prior's support.
We adopt the standard empirical process notation as follows. For a function $h$ of a random vector $Z=(Y,D, X^\top)^\top$ that follows distribution $P$, we let $P[h]=\int h(z)\mathrm{d}P(z)$, $\mathbb{P}_n[h]=n^{-1}\sum_{i=1}^{n}h(Z_i)$, and $\mathbb{G}_n[h]=\sqrt n\left(\mathbb{P}_n-P\right)[h]$. The next set of assumptions restrict the complexity of the conditional mean functions. The first part requires the class $\{m_\eta:\eta\in\mathcal{H}_n\}$ to be Glivenko-Cantelli, plus some mild moment conditions on its envelope function. The second part imposes the stochastic equicontinuity, which is holds when the conditional mean function belong to a Donsker class. For the H\"older class considered in Section (ref), this enforces the sufficient smoothness of those functions relative to the dimensionality of covariates.
The next assumption concerns the prior stability condition, which is common to semiparametric Bayesian inference ghosal2017fundamentals. This facilitates the technical proof for which we need to consider the perturbation along the least favorable direction. For standard parametric models, the absolute continuity of the prior density suffices. However, for nonparametric priors, the very notion of a Radon-Nikodym density is non-trivial, and one needs to apply the Cameron-Martin theorem; see Proposition I.20 in ghosal2017fundamentals. For that purpose, we introduce some necessary terminologies related to the general Gaussian process. Such a process determines a so-called reproducing kernel Hilbert space (RKHS) $\left(\mathbb{H}^m,\|\cdot\|_{\mathbb{H}^m}\right)$.
Our Bayesian method based on standard Gaussian process priors in Algorithm (ref) does not include a correction involving the Riesz representer $\gamma_0$ as defined in (ref). Yet to establish prior stability, an approximation condition for $\gamma_0$ is imposed, requiring sufficient regularity of the propensity score $\pi_0(\cdot)$. We introduce the ball in $\mathbb{H}^m$ centered at the true Riesz representer $\gamma_0$ given by
for some rate $r_n$, where $\Vert \cdot\Vert_{\infty}$ denotes the supremum norm.
Assumption (ref) imposes an approximation condition to the Riesz representer $\gamma_0$ via the restriction $\overline\gamma_n\in\mathbb{H}^m(\zeta_n)$. Based on this assumption, we provide the proof of this prior stability in Supplementary Appendix (ref). In comparison, the prior correction weakens the requirement with the help of a pilot estimator of the propensity score, pioneered by ray2020causal. Under propensity score adjusted priors analyzed in the next section, BLY2022 the approximation condition even holds under double robustness.
We now establish a Bernstein-von Mises Theorem for our nonparametric Bayesian method based on standard Gaussian process priors. When it comes to the centering point of the posterior, we consider an asymptotically efficient estimator $\widehat{\tau}$ with the following linear representation:
where $\widetilde{\tau}_0=\widetilde{\tau}_{\eta_0}$ is the efficient influence function given in (ref). Below, we write $\mathcal{L}_{\Pi}(\sqrt{n}(\tau_\eta-\widehat{\tau})|Z^{(n)})$ for the marginal posterior law of $\sqrt{n}(\tau_\eta-\widehat{\tau})$.
This asymptotic equivalence result is established using the so called bounded Lipschitz distance. For two probability measures $P,Q$ defined on a metric space $\mathcal{Z}$, we define the bounded Lipschitz distance as
where $BL(1)=\left\{f:\mathcal{Z}\mapsto\mathbb{R}, \sup_{z\in\mathcal{Z}}|f(z)|+\sup_{z\neq z'}\frac{|f(z)-f(z')|}{\|z-z'\|_{\ell_2}}\leq 1 \right\}$. Here, $\|\cdot\|_{\ell_2}$ denotes the vector $\ell_2$ norm. Below is our main statement about the asymptotic behavior of the posterior distribution of $\tau_{\eta}$, that is derived from the Bayes rule given the prior specification and the observed data $Z^{(n)}$. As in the modern Bayesian paradigm, the exact posterior is rarely of closed-form, and one needs to rely on certain Monte Carlo simulations, such as the implementation procedure in Section (ref), to approximate this posterior distribution, as well as the resulting point estimator and credible set.
Theorem (ref) establishes the BvM result for our Bayesian procedure using standard Gaussian process priors. The entropy condition uniformly over $\eta\in\mathcal{H}_n$ is satisfied if $m_\eta$ is sufficiently smooth, that is, if $m_\eta$ belong to a fixed $F_0$-Donsker class and, in particular, rules out double robustness. On the other hand, note that the asymptotic equivalence is obtained without any adjustment of prior or correction to posterior distributions, so the full Bayesian flavor is preserved.
In this section, we show that the Bayesian procedure in Algorithm (ref), which employs prior and posterior adjustments, satisfies the Bernstein-von Mises Theorem under double robust smoothness conditions. Herein, we clarify the notion of double robustness. In the earlier development, the focus is typically on developing working parametric models for either the propensity score or the conditional mean function, and the double robust estimation hedges against the risk of model misspecification. However, implausible parametric assumptions on the data generating process are of limited applicability to complex phenomena in economics. Recent advances in the double machine learning literature have led to a number of important developments in causal inference, utilizing flexible nonparametric or machine learning algorithms. In this context, double robustness means the possibility to trade off the estimation accuracy between nuisance functions.
Below, we present the assumptions that enable double-robust inference through our propensity score adjustments to the prior and posterior distributions.
Assumption (ref) imposes sufficiently fast convergence rates for the estimators for the conditional mean function $m_0$ and the propensity score $\pi_0$. The posterior convergence rate for the conditional mean can be derived by modifying the classical results of ghosal2000rates by accommodating the propensity score-adjusted prior, in the same spirit of ray2020causal. We refer to BLY2022 who showed that this assumption allows for double robustness under H\"older type smoothness assumptions.
Assumption (ref) restricts the functional class $\mathcal{H}_n$ to form a $P_0$-Glivenko-Cantelli class; see Section 2.4 of van1996empirical and imposes a stochastic equicontinuity condition on a product structure involving $\widehat\gamma$ and $m_\eta$. Hence, the complexity of the functional class $(m_{\eta}-m_0)$ can be compensated by certain high regularity of the corresponding Riesz representer and vice versa. This condition adapts the complexity requirement of BLY2022 by only restricting the control arm.
Recall the propensity score-dependent prior on $m$ given in (ref), i.e., $m(\cdot) = q^{-1}\left(W^m(\cdot) + \lambda\widehat \gamma(\cdot)\right)$. Below, we restrict the behavior for $\lambda$ through its hyperparameter $\varsigma_n>0$. For two sequences $\{a_n\}$ and $\{b_n\}$ of positive numbers, we write $a_n \lesssim b_n$ if $\limsup_{n\to\infty} (a_n / b_n)<\infty$, and $a_n \sim b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n$.
Assumption (ref) incorporates Conditions (3.9) and (3.10) from Theorem 2 in ray2020causal, and it is imposed to establish the stability property of the adjusted prior distribution. We will provide sufficient conditions for Assumption (ref) in Section (ref).
We now establish a semiparametric Bernstein–von Mises theorem for our double robust Bayesian procedure given in Algorithm (ref).
Theorem (ref) shows that, under double-robust smoothness conditions, the BvM theorem holds only up to a “bias term" $b_{0,\eta}$, which depends on the unknown conditional mean $m_0$. This biased posterior makes the BvM not feasible in practice. We also emphasize that the derivation of this result is different to the BvM results in BLY2022, as we need to control the denominator in the asymptotic expansions.
The next result is an immediate implication of Theorem (ref). Specifically, it provides a Bernstein-von Mises Theorem for Bayesian procedures that do not rely on posterior correction. This can be achieved if the bias term is asymptotically negligible uniformly over the class of hyperparameters, which requires more restrictive smoothness conditions on the conditional mean function $m_0$.
While Corollary (ref) allows for arbitrarily low regularity of propensity scores, it requires the conditional mean function to be sufficiently smooth; specifically, the smoothness of $m$ must be greater than or equal to $\dim(X_i)/2$ (also referred to as the Donsker property). This condition is also called single robustness by ray2020causal, and indeed, this corollary extends their findings to the inference on the ATT. Also, as they point out, propensity score adjusted priors (ref) relax the uniformity condition $\sup_{\eta\in\mathcal H_n}|\mathbb{G}_n\left[m_\eta-m_0\right]|=o_{P_0}(1)$ used in Theorem (ref) under standard Gaussian process priors.
Under double robust assumptions, however, the Bayesian procedure that achieves the BvM equivalence in Theorem (ref) is not feasible, because it depends on the term $b_{0,\eta}$, which is a function of the unknown conditional mean $m_0$. Our objective is to maintain double robust conditions, while considering pilot estimators for the unknown functional parameters in $b_{0,\eta}$. The correction term $\widehat{b}_{\eta}$, as introduced in (ref), results in a feasible Bayesian procedure that satisfies the BvM theorem, as demonstrated below.
Theorem (ref) shows that the Bayesian method proposed in Algorithm (ref), $\check\tau_\eta = \tau_{\eta} - \widehat{b}_\eta$, achieves the BvM result under double robust smoothness conditions. The following remark clarifies the relationship when considering posterior correction alone, in which case BvM results are available only under more restrictive smoothness assumptions on the propensity score and the conditional mean function.
In this section, we provide primitive conditions for the assumptions used to derive the BvM Theorems. To do so, we focus on squared exponential process priors as an example of Gaussian process priors. Moreover, we consider specific smoothness classes to derive the explicit regularity conditions implied by our high-level assumptions.
A Gaussian process (GP) is completely characterized by its mean and covariance functions rassmusen2006gaussian. Below we consider a GP prior, which has mean zero and the covariance function specified by $\mathbb{E}[W(s)W(t)]=\exp(-\Vert s-t\Vert_{\ell_2}^2)$. This so-called squared exponential process prior, which is one of the most commonly used priors in applications; see rassmusen2006gaussian and murphy2023pml. Following BLY2022, we consider a rescaled Gaussian process $\big(W(a_nt):\,t\in [0,1]^p\big)$. Intuitively, $a_n^{-1}$ can be thought as a bandwidth parameter. For a large $a_n$, the prior sample path $t\mapsto W(a_nt)$ is obtained by shrinking the long sample path $t\mapsto W(t)$. Thus, it incorporates more randomness and becomes suitable as a prior model for less regular functions, see van2008gaussian,van2009adaptive.
Below, $\mathcal{C}^{s_m}([0,1]^p)$ denotes a H\"older space with the smoothness index $s_m>0$. Specifically, we illustrate our theory with the case where $m_0\in \mathcal{C}^{s_m}([0,1]^{p})$. Given such a H\"older-type smoothness condition, we choose
The particular choice of $a_n$ mimics the corresponding kernel bandwidth based on any kernel smoothing method. Note that the minimax posterior contraction rate for the conditional mean function $m_\eta$ given by $\varepsilon_n=n^{-s_m/(2s_m+p)}(\log n)^{s_m(1+p)/(2s_m+p)}$; see Section 11.5 of ghosal2017fundamentals.
Proposition (ref) makes explicit the smoothness requirements for the BvM Theorem to hold when standard Gaussian process priors are placed on the conditional mean function $m$. This result shows that the smoothness of both the conditional mean function and the propensity score function must exceed $\dim(X)/2$. Conversely, in situations where one is confident that these regularity conditions are met, no additional modifications to the Bayesian procedures are necessary to achieve the BvM result.
Proposition (ref) requires $\sqrt{s_\pi \, s_m}>p/2$, which represents a trade-off between the smoothness requirement for $m_0$ and $\pi_0$. This corresponds to the double robustness property; i.e., a lack of smoothness of the conditional mean function $m_0$ can be mitigated by exploiting the regularity of the propensity score$\pi_0$, and vice versa.
This section investigates the finite sample performance of the proposed Bayesian estimation/inference approaches and then apply them to the well-known DiD study of card1994minimum.
We now present Monte Carlo simulation results to compare our proposed semiparametric Bayesian methods with existing frequentist approaches. Consider the following data generating process (DGP) for observed variables $(Y_{i1}, Y_{i2}, D_i, X_{i}^\top)^\top$ given by
where the covariance matrix $\Sigma=(\Sigma_{jk})_{1\leq j,k,\leq p}$ is determined by $\Sigma_{jk}=0.5^{|k-j|}$. We generate outcomes in two periods:
We consider the following four different designs based on different specifications of the functions $g$ and $h$:
The fixed effect $\alpha_i$ and the error terms are standard normal, with $(\alpha_i, \epsilon_{i1}, \epsilon_{i2}(0), \epsilon_{i2}(1))^\top \sim \mathcal{N}(0, I_4)$, where $I_4$ denotes the four-dimensional identity matrix.\footnote{Tables (ref) and (ref) in the Supplementary Appendix (ref) present additional simulation results for cases where the error terms follow chi-squared distributions or include heteroskedasticity. The finite-sample performance in these cases is similar to that observed in Table (ref) for standard normal errors.} The true ATT is zero for all cases. Our DGPs follow the structure of DGP1 in sant2020doubly but allow for more covariates and a possibly nonlinear function $h(\cdot)$. In our simulations, we analyze the effect of varying the dimension of covariates $p\in\{5,10, 20\}$ and varying sample sizes $n\in\{500, 1000\}$. Throughout our simulations, the number of Monte Carlo replication is set to $1000$.
Our nonparametric Bayesian (hereinafter Bayes) and the double robust Bayesian (DR Bayes) methods are implemented following Algorithms (ref) and (ref) in Section (ref), using the MATLAB package $\mathtt{GPML}$ to draw posteriors. Both Bayesian methods are implemented based on $B=5000$ posterior draws. Here, we apply the exponential family specification in ((ref)) to the Gaussian case. The resulting posterior distribution of the conditional mean function is available in closed form (see Supplementary Appendix (ref) for details), eliminating the need for computationally costly Monte Carlo samplers like MCMC.\footnote{If the conditional density function in ((ref)) belongs to another distribution in the exponential family, the posterior of the conditional mean can also be approximated using an analytical approximation, such as the Laplace method, see riihimaki2014laplace.} The tuning parameter $\varsigma_n$ for DR Bayes, which corresponds to the standard deviation of the adjusted prior, is set according to the prior specification step in Algorithm (ref). In Supplementary Appendix (ref), Table (ref) demonstrates that the performance of DR Bayes is stable with respect to the value of $\varsigma_n$. DR Bayes in Table (ref) uses the full sample twice in computing the prior/posterior adjustments and the posteriors of the conditional mean function. As shown in Table (ref) in Supplementary Appendix (ref), results from sample splitting are comparable to those in Table (ref).
We also compare the Bayesian methods to several frequentist DiD estimators. DR corresponds to the improved doubly robust DiD estimator proposed by sant2020doubly. OR, the outcome regression approach, refers to the sample analog of ((ref)) where the conditional mean $m_0$ is estimated by a linear regression of $\Delta Y_i$ on $X_i$ using the sample of the control arm. Two types of inverse propensity score weighted (IPW) estimators are considered: IPW$^{\text{HT}}$ refers to the IPW estimator in abadie2005, which is of the horvitz1952generalization type. IPW$^{\text{H\'{a}jek}}$ refers to the hajek1971discussion type IPW estimator that normalizes the weights to sum up to one.\footnote{The expression for the H\'{a}jek--type IPW DiD estimator is given in equation (4.1) of sant2020doubly. } TWFE corresponds to the standard two-way-fixed effect model that regresses $Y_{it}$ on $D_i$, $t$, the interaction $D_i\times t$ and $X_i$. DML corresponds to the double/debiased machine learning ATT estimator of chernozhukov2017double or chang2020double, where the nuisance function $\pi_0$ is estimated by logistic LASSO and $m_0$ estimated by random forests.\footnote{We apply random forest to estimate $m_0$ to cope with the nonlinear function forms in Designs II and IV, and to match DML with our Bayesian procedures that estimate $m_0$ nonparametrically. Frequentist DiD estimators, except for DML, are implemented using the R package $\mathtt{DRDID}$, while DML is implemented using the R package $\mathtt{DoubleML}$.} Table (ref) presents the finite sample (mean) bias of the point estimator, coverage probability (CP) and the average length (CIL) of the $95\%$ credible/confidence interval for the Bayesian and frequentist methods mentioned above.
Concerning the Bayesian DiD for estimating the ATT, Table (ref) shows that the nonparameteric Bayes performs well in Design I, but undercovers by $9\%$ to $14\%$ in Design II when $p=10$ and $20$. DR Bayes improves the coverage probability of nonparameteric Bayesian inference in these cases and performs well across both designs, different dimensions $p$, and sample sizes $n$. The point estimator produced by DR Bayes also leads to a smaller bias than the nonparametric Bayes. On the other hand, the nonparametric Bayes yields shorter confidence intervals. We also see that for large values of $p$, nonparametric Bayes tends to undercover.
In Table (ref), the frequentist DiD estimators DR, two types of IPW, and DML -- each of which is double robust or at least robust to misspecification in the conditional mean function -- exhibit good coverage performance in both designs. Among them, DR produces slightly longer CIs than our DR Bayes in most cases; IPW estimators yield longer CIs than most of other methods, including both Bayesian procedures; and DML yields slightly shorter CIs than DR Bayes in Design I but noticeably longer CIs in Design II. Unsurprisingly, OR suffers from severe undercoverage in Design II. TWFE performs poorly in both designs, where the time trend is linearly correlated with the covariates $X$, as also documented in sant2020doubly.
Table (ref) presents the finite sample performance of aforementioned DiD estimators when the function $g(\cdot)$, used to specification of the propensity score, is nonlinear. DR Bayes, DR, and IPW estimate the propensity score using logistic regression, while DML uses logistic LASSO. Table (ref) illustrates the impact of misspecifying the propensity score on these methods. DR Bayes maintains reasonably good performance in Design III. Although its performance deteriorates in Design IV, particularly as dimensionality of $X$ increases, it still outperforms frequentist estimators, including double-robust methods like DR and DML. Nonparametric Bayes, using standard Gaussian process priors, avoids estimating the propensity score under misspecification and hence performs well in Design III while outperforming all other methods in Design IV.
We apply Bayesian DiD methods to the well-known minimum wage study of card1994minimum.\footnote{The data is available on \href{https://davidcard.berkeley.edu/data_sets}{https://davidcard.berkeley.edu/data_sets}. } The outcome variables $Y_{1i}$ and $Y_{2i}$ are full time equivalent (FTE) employment of fast-food stores in New Jersey and Pennsylvania before and after New Jersey's raise of minimum wage. The treatment variable takes one for fast-food stores in New Jersey and zero otherwise. The set of covariates $X$ includes the twelve store characteristics surveyed before the minimum wage change: indicator for company ownership, three chain type dummies, numbers of managers, cash registers and hours open per weekday, time to the first wage raise, indicator for offering recruitment bonus, item prices of medium soda, small french fries and a main course. We would like to see whether the findings of card1994minimum which considers as store characteristics the company ownership indicator and chain type dummies in their regression-adjusted model would change if more covariates are included and a flexible functional form of $m_0(x)$ is allowed.
The sample size is $307$. Since the data contains a non-negligible proportion of units with propensity score estimates very close to $1$, we follow crump2009dealing and discard observations with the estimated propensity score outside the range $(0, 1-t]$, with the trimming threshold $t\in\{0.05, 0.01\}$.
Table (ref) presents the ATT estiamtes for Bayesian and frequentist methods. As we see, all methods produce positive but insignificant ATT estimates for the impact of minimum wage, which is in line with the findings of card1994minimum. For example, nonparametric Bayes and DR Bayes yield ATT point estimates ranging from $1.907$ to $2.024$ and confidence intervals covering $0$ with the length from $5.514$ to $6.683$. Bayesian methods also provide shorter confidence intervals than most of the frequentist methods including the widely-used TWFE estimator, except that the credible interval produced by semiparametric Bayes is slightly longer than the confidence interval of H\'{a}jek--type IPW when $t=0.05$.
If we do not trim the propensity score, the failure of the overlap condition prevents us from using estimators that involve the inverse propensity score. Among estimators that do not use the propensity score, nonparametric Bayes gives an ATT estimate of $1.935$, with a 95% confidence interval of $[-0.460, 4.341]$, for the full sample without any trimming ($t=0,\bar n_t=249, \bar n_c=58$). OR yields an estimated ATT of $3.351$, with a 95% CI of $[-1.233, 7.936]$. TWFE provides an estimated ATT of $2.635$, with a 95% CI of $[-0.622, 5.891]$. It turns out that our semiparametric Bayesian method continue to yield stable results when the overlap condition is nearly violated. In sum, our Bayesian methods, which allows a flexible form of the conditional mean function $m_0(x)$ as well as a rich set of covariate, generate comparable ATT estimate with the original findings in card1994minimum.\footnote{When covariates are not included in the model, card1994minimum report the difference-in-difference estimate of $2.76$ (standard error $1.36$), and the regression adjusted model with controls for chain and ownership dummies yields an estimate of $2.30$ (standard error $1.20$).} Therefore, our Bayesian DiD methods confirms the robustness of findings in the classic literature against model specifications.
We now provide extensions to the canonical DiD panel data setup and show that our Bayesian DiD methods, described in Section (ref), can be conveniently extended to cases such as multiple periods with staggered entry and repeated cross sections.
The Bayesian DiD methods described in Section (ref) can be conveniently extended to the cases with multiple periods and staggered intervention de2020two,callaway2021difference,sun2021estimating,borusyak2024revisiting. The related literature focuses on the identification of disaggregated causal parameters and some proper aggregation of these parameters. This section extends our Bayesian method for inference on the disaggregated ATT, specifically the group-time ATT proposed by callaway2021difference.
Suppose the available panel data consists of $T$ periods indexed by $t=1,\dots,T$ and the earliest treatment intervention occurs at period $S$. We assume that the treatment intervention remains once a unit get treated. As a result, the entire path of treatment assignment for each unit can be summarized by his/her first treated period (cohort), denoted by the cohort variable $G_i\in\{S,\dots,T,\infty\}$, where $G_i=\infty$ means the unit $i$ never gets treated. Let the cohort indicators $D_{ig}$ denote whether unit $i$ first receives treatment in period $g \in \{S, \ldots, T, \infty\}$, where $D_{i\infty}=1$ indicates that unit $i$ never receives treatment.
We assume that never-treated units exist. The potential outcomes depend on cohorts and thus are denoted as $Y_{it}(g)$ for $g\in\{S,\dots,T\}$ and $Y_{it}(0)$ for $G_i=\infty$. Obviously, $\sum_{g=S}^T D_{ig} + D_{i\infty}=1$. The realized outcome for unit $i$ at time $t$ is $ Y_{it}=Y_{it}(0)+\sum_{g=S}^{T}D_{ig}\left(Y_{it}(g)-Y_{it}(0)\right)$.
We focus on the analysis of treatment effect heterogeneity by allowing the ATT to vary with the cohort $g$ ($g \neq \infty$) and the time period $t \geq g$:
Suppose a vector of pre-treatment covariates $X_i$ is also available, a vector of dimension $p$, with the distribution $F_0$ and the density $f_0$. The researcher observes independent and identically distributed observations of $(Y_{i1},\dots,Y_{iT}, D_{iS},\dots,D_{iT},D_{i\infty}, X_i)$, $i=1,\ldots, n$.
Applying the identification strategy in callaway2021difference, the ATT parameters $\tau_0^{g,t}$ for $g\in\mathcal G:=\{S,\dots ,T\}$ and $t=g,\dots,T$ can be identified under Assumption (ref) below.\footnote{callaway2021difference propose two identification strategies, depending on the whether the parallel trend assumption is imposed on the never-treated cohort or “Not-Yet-Treated" cohorts. Here we consider the former version.} For the identification of the ATT, we follow the setup by callaway2021difference and impose the following conditions, which correspond to their Assumptions 3, 4, and 6 (in their $\delta = 0$ case).
Assumption (ref)(i) is a “no anticipation" assumption, Assumption (ref)(ii) is a conditional parallel trend assumption based on the never-treated cohort, and Assumption (ref)(iii) is an overlap restriction. Under Assumption (ref), callaway2021difference show that the ATT in the staggered entry case is identified by
where the difference operator $\Delta_g$ is defined by $\Delta_g Y_t := Y_t - Y_{g-1}$ and the conditional mean function $m_0^{g,t}(x):= \mathbb{E}_0\left[\Delta_g Y_t\mid D_{\infty}=1, X=x\right]$.
The identification result in ((ref)) uses the cohort $g$ (i.e., $D_{g}=1$) as the treated group and the “never treated" cohort ($D_{\infty}=1$) as the control group. Using the transformed cross-sectional data $\left(\Delta_g Y_{it}, D_{ig}, X_{i} \right)$ for $i=1,\dots, n$ and following the notation in Section (ref), we can write ATT for a given pair $(g,t)$ under a family of probability distributions $\{P_\eta:\eta\in\mathcal H\}$ as
where $\mathbb{E}_{\eta}$ denotes the expectation with respect to the distribution of $\left(\Delta_g Y_t , D_g, X\right)$. The Bayesian DiD procedures in Section (ref) can be applied in the staggered DiD case to obtain the posterior draws $\left\{(\tau_{\eta}^{g,t})^s: s=1,\ldots,B \right\}$. Specifically, this can be achieved by replacing $\Delta Y_i$, $D_i$, $m_{\eta}(\cdot)$, $\pi_{\eta}$, and $\pi_{\eta}(\cdot)$ in Algorithm (ref) or (ref) by $\Delta Y_{it}$, $D_{ig}$, $ m_{\eta}^{g,t}(\cdot)$, $\pi_{\eta}^g:=\mathbb{E}_{\eta}[D_g]$ and $\pi_{\eta}^g(\cdot):=P_\eta\left(D_g=1 \mid D_g + D_{\infty}=1, X=\cdot \right)$, respectively, as defined in this section.
The first resulting Bayesian estimator is denoted by $\tau_\eta^{g,t}$, while the second, double-robust Bayesian method is denoted by $\check{\tau}_\eta^{g,t}$ for a cohort $g \in \mathcal{G}$. The next result is an immediate implication of Theorem (ref) and Theorem (ref), and its proof is thus omitted.
Corollary (ref) pertains to inference on cohort-specific ATTs and establishes BvM results for our two Bayesian methods, employing either standard Gaussian process priors or prior/posterior adjustments via the cohort-specific propensity score. We note that extending this framework to aggregate ATTs is highly non-trivial, as it requires a joint modeling of outcome variables across different cohorts and time periods. Hence, distinct prior and likelihood specifications in the Bayesian methodology, as well as prior/posterior adjustments of the double robust version, are needed. A thorough investigation is therefore left for future research.
Our method also allows for repeated cross-sections following abadie2005, as also considered by sant2020doubly. In this case, we consider a dummy variable $T_i$ that takes the value two if observation $i$ is only observed in the post-treatment period, and one if observation $i$ is only observed in the pre-treatment period. Define $Y_i = (T_i-1) Y_{i2} + \left(2 - T_i\right) Y_{i1}$. The available data is $\left\{Y_i, D_i, T_i, X_i\right\}_{i=1}^n$. Let $n_2$ and $n_1$ be the sample sizes for the post-treatment and pre-treatment periods, respectively, such that $n = n_2 + n_1$; let $\mathbb{P}(T = 2) \in (0,1)$. The following assumption restates Assumption 3.3 of abadie2005.
Under Assumptions (ref), we can write
Then using Assumption (ref), we can identify ATT as
where $m_0 (x,t)\equiv \mathbb{E}_0[Y\mid D=0, X=x,T=t]$ for $t=1,2$.
With an analogous reparametrization as in the panel data case, we obtain
Interestingly, the analysis of the last conditional expectation involves a difference of conditional moment function as for the average treatment effect and can be analyzed similarly to BLY2022 in absence of prior corrections. For our double robust method in repeated cross-sections, we emphasize that the efficient influence function takes a different functional form (see sant2020doubly). This translates to a modified prior and posterior adjustments of our double robust Bayesian procedure in Algorithm (ref). While this procedure would be analogous to our double robust method, a full derivation of its asymptotic properties lies beyond the scope of this paper.
This paper introduces new semiparametric Bayesian procedures that satisfy the Bernstein-von Mises results in the DiD setup. Our first proposal, based on standard Gaussian process priors, provides a Bayesian analog to the outcome regression in heckman1997matching. Through simulations, we show that it performs well in models that are not overly complex and, since no propensity score specification is required, it is not sensitive to the overlap issues. Our second, double robust proposal incorporates prior/posterior corrections based on estimated propensity scores. In simulations it works well for complex models, i.e., when the number of covariates is large. Overall, our Bayesian methods exhibit remarkable finite sample performance, while adapting to the functional form of the conditional mean function. Although our focus is primarily on the DiD panel data case, we also discuss extensions to the repeated cross-section case and staggered interventions.