$\texttt{rdid}$ and $\texttt{rdidstag}$: Stata commands for robust difference-in-differences

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.

54,699 characters · 31 sections · 38 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

rdid and rdidstag: Stata commands for robust difference-in-differences

\inserttype[st0001]{article} { Kyunghoon Ban\\Rochester Institution of Technology\\Rochester, NY\\ [email removed] \and D\'esir\'e K\'edagni\\UNC-Chapel Hill\\Chapel Hill, NC\\[email removed] }

abstractThis article provides a Stata package for the implementation of the robust difference-in-differences (RDID) method developed in bk2023rdid. It contains three main commands: rdid, rdid_dy, rdidstag, which we describe in the introduction and the main text. We illustrate these commands through simulations and empirical examples.

Introduction

In this article, we present the rdid, rdid_dy, rdidstag commands for estimation and inference on robust difference-in-differences (DID) bounds as developed by bk2023rdid. These commands, summarized in Table (ref), allows one to analyze the average treatment effects on the treated (ATT) under the canonical $2 \times 2$ DID setting where observational data consist of two different groups (a treatment group and a control group) and two time periods (pre-treatment and post-treatment) as well as the staggered adoption design where there are multiple cohorts with different timings in the treatment adoption.

table[table omitted — 884 chars of source]

Our article contributes to the expanding suite of publicly available software within the DID design community, including villa2016diff, houngbedji2016abadie, de2019fuzzy, rios2022drdid, rios2023csdid, among others. Specifically, our software is more closely aligned with research exploring alternative identifying assumptions within the DID framework; for instance, mora2015didq examines DID estimations under alternative assumptions and assesses their robustness, while bravo2022honestdid discusses methods to bound post-treatment parallel trends (PT) violations using pre-treatment PT violations.

The rest of this article is structured as follows. In section 2, we recap the underlying frameworks of the robust DID bounds by bk2023rdid. In sections 3, 4, and 5, we describe the rdid, rdid_dy, rdidstag commands, respectively. In section 6, we illustrate the performance of each command through a series of Monte Carlo simulations.

Framework and methodology

Consider the following potential outcome model:

eqnarray[eqnarray omitted — 53 chars of source]

where $(Y_t,D)$, $t\in \mathcal T_0 \cup \{1\}$ represents the observed data, while the vector $(Y_1(0), Y_1(1))$ is latent. The variables $Y_{t_0}, Y_1 \in \mathcal Y$ are respectively the observed outcomes in the baseline period $t_0 \in \mathcal T_0$ and the follow-up period 1, while $D\in \left\{0,1\right\}$ is the observed treatment that occurred between periods 0 and 1, $Y_1(0)$ and $Y_1(1)$ are the potential outcomes that would have been observed in period 1 had the treatment $D$ been externally set to 0 and 1, respectively. $\mathcal T_0$ denotes the set of all pre-treatment periods. Model ((ref)) assumes that there is no anticipatory effect of the treatment, so that $Y_{t_0}(1)=Y_{t_0}(0)$ for all baseline periods $t_0 \in \mathcal T_0$.

The standard DID estimand is defined as the difference between the OLS estimand in period 1 and the selection bias in period 0: $$\theta_{DID} \equiv \theta_{OLS}-SB_0$$ where $\theta_{OLS}$ is the coefficient of a regression of $Y_1$ on the binary $D$, $$ \theta_{OLS} = \mathbb{E}[Y_1 \vert D=1 ] - \mathbb{E}[Y_1 \vert D=0 ], $$ and $SB_t$ is the selection bias in period $t$,

align*[align* omitted — 82 chars of source]

for $t\in \mathcal T_0 \cup \{1\}$. Recall that $SB_0$ is identified with the no anticipatory effect but $SB_1$ is an unobserved counterfactual.

Note that the average treatment effect on the treated (ATT), $$ATT\equiv \mathbb E[Y_1(1)-Y_1(0)\vert D=1],$$ is identified by $\theta_{DID}$ under the parallel trend (PT) assumption, $$ \mathbb{E}[Y_1(0) - Y_0(0) \vert D=1] = \mathbb{E}[Y_1(0) - Y_0(0) \vert D=0], $$ or $SB_1 = SB_0$ (bias stability).

Hence, the identification of ATT can be re-framed as imputing the counterfactual selection bias $SB_1$. This observation inspires us to introduce a robust version of the DID estimand.

Define the baseline information set $\mathcal I_0\equiv \left\{(t_0,x_{t_0}): x_{t_0} \in \mathcal X_{t_0}, t_0 \in \mathcal T_0\right\}$, where $\mathcal X_{t_0}$ is the support of a baseline covariate $X_{t_0}$ in period $t_0$, and $\mathcal T_0$ is the set of pre-treatment periods. For simplicity, suppose there are no baseline covariates so that $\mathcal I_0=\mathcal T_0$.

definitionGiven the baseline information set $\mathcal I_0=\mathcal T_0$, we define the robust difference-in-differences (RDID) estimand as \begin{eqnarray} \theta_{RDID}\equiv \theta_{OLS}-\operatorname{Conv}\big(SB(\mathcal T_0)\big), \end{eqnarray} where $\operatorname{Conv}(A)$ is the convex hull of a generic set $A$, $SB$ is a function defined as \begin{eqnarray*} SB\colon \quad \mathcal I_0& \to & \mathbbm R,\\ \iota_0& \mapsto & SB(t_0)=SB_{t_0}, \end{eqnarray*} and $SB(\mathcal T_0)$ is the image of the information set $\mathcal I_0=\mathcal T_0$ through the function $SB$.

When baseline covariates exists, we define the function $SB$ as $$

array[array omitted — 221 chars of source]

In the above definition, if there is only one element in the information, i.e., $\mathcal I_{0}=\mathcal T_0=\{0\}$, then $\operatorname{Conv}\big(SB(\mathcal T_0)\big)=SB_0$, and the robust DID estimand is the same as the standard DID estimand.

Bias set stability

Let us first consider the simple case with no covariates in the model.

assumption[Bias set stability] \begin{eqnarray*} SB_1 \in \left[ \inf_{\iota_0\in \mathcal I_0} SB(\iota_0), \sup_{\iota_0 \in \mathcal I_0} SB(\iota_0)\right] \equiv \Delta_{SB(\mathcal I_0)}, \end{eqnarray*} where $SB(\iota_0)$ is defined above.

Assumption (ref) is weaker than the standard “parallel/common trends” assumption. Indeed, if $\mathcal{I}_0$ is the singleton of a single baseline information $I_0=\{0\}$, then Assumption (ref) is equivalent to $SB_1= SB_0,$ which is equivalent to the parallel trends assumption.

Under Assumption (ref), robust DID estimand yields the following bounds on the ATT.\footnote{Note that the RDID can be generalized by using other correspondences $f(SB(\mathcal I_0))$ than $\operatorname{Conv}(SB(\mathcal I_0))$ associated with other assumptions bk2023rdid.}

propositionSuppose that model ((ref)) along with Assumption (ref) holds. Then, the following bounds hold for the ATT: \begin{eqnarray*} ATT \in \left[ \theta_{OLS}- \sup_{\iota_0 \in \mathcal I_0} SB(\iota_0), \theta_{OLS}- \inf_{\iota_0 \in \mathcal I_0} SB(\iota_0) \right] \equiv \Theta_I. \end{eqnarray*} These bounds are sharp, and $\Theta_I$ is the identified set for the ATT.

The bounds in Proposition (ref) are never empty, as they always contain the standard DID estimand under the parallel trends assumption. However, they may not contain the OLS estimand in period 1, $\theta_{OLS}$, as 0 may not lie within the set $\Delta_{SB(\mathcal I_0)}$. If all pre-treatment periods selection biases are equal, i.e., $SB(\iota_0)=SB_0$ for all $\iota_0$, then our bounds collapse to a point, the standard DID estimand. In case the information set $\mathcal I_0$ is the set of pre-treatment periods, the above bounds are robust to violations of PT that can be captured in the pre-treatment periods.

PO-RDID: Best predictor of $SB_1$ based on a loss function

Consider a random variable $I_0$ representing the baseline information with a distribution $\mu_{I_0}$. Provided $\mu_{I_0}$ is known, we are going to assume that the decision maker will choose the selection bias $SB_1$ in such a way that a loss function is minimized. By plugging such an optimal selection bias $SB_1$ into the definition of the robust DID, we obtain what we call a policy-oriented robust difference-in-differences (PO-RDID) estimand. This estimand may not have a causal interpretation, but it may help the policy-maker in her decision making process.

assumptionLet $\mathcal L(b; SB(\mathcal I_0), \mu_{I_0})$ be the decision maker's loss function when she observes a (random) baseline information $I_0$ and chooses a selection bias $b$. The decision maker chooses $b$ (or $SB_1$) to minimize the loss $\mathcal L(b)$: $SB_1^{opt}=\arg\min_b \mathcal L(b; SB(\mathcal I_0), \mu_{I_0})$.

In this paper, we consider the class of $p$-norm losses defined as: $$\mathcal L_p(b; SB(\mathcal I_0), \mu_{I_0})=\left(\mathbb E_{\mu_{I_0}}\left[\vert b - SB(I_0)\vert ^p \right]\right)^{1/p},$$ where $1 \leq p \leq \infty$. We are going to derive the optimal selection bias $SB_1^{opt}$ for $p\in \left\{1,2,\infty\right\}$. We consider those special loss functions because the solutions to the optimization problem have closed-form expressions. Other loss functions can also be considered.

The command rdid with an option rdidtype(1) can be used, and the distribution of $SB_0(I_0)$ is estimated by the proportion of observations at each level of $I_0$.

$L1$ loss: Mean absolute error (MAE)

Note that $\mathcal L_1(b; SB(\mathcal I_0), \mu_{I_0})=\mathbb E_{\mu_{I_0}}\left[\vert b - SB(I_0)\vert \right]$. Given this $L1$ loss function, under Assumption (ref), the decision maker solves the following optimization problem:

eqnarray*[eqnarray* omitted — 90 chars of source]

The optimal decision is to set the selection $SB_1$ to be equal to the median selection bias in the baseline period, i.e., $SB_1^{opt}=Med_{\mu_{I_0}}\big(SB(I_0)\big)$. In such a case, the policy-oriented robust DID estimand is given by $$\theta_{PO-RDID}^{\mathcal{L}_{1}}=\theta_{OLS}-Med_{\mu_{I_0}}\big(SB(I_0)\big).$$

$L2$ loss: Root mean square error (RMSE)

We have $\mathcal L_2(b; SB(\mathcal I_0), \mu_{I_0})=\left(\mathbb E_{\mu_{I_0}}\left[\vert b - SB_0(I_0)\vert ^2 \right]\right)^{1/2}$, and minimizing the RMSE is equivalent to minimizing the mean square error (MSE). Therefore, under Assumption (ref), the decision maker solves the following optimization problem:

eqnarray*[eqnarray* omitted — 94 chars of source]

This yields an optimal decision for the selection $SB_1$ to be set equal to the average selection bias in the baseline period, i.e., $SB_1^{opt}=\mathbb E_{\mu_{I_0}}[SB(I_0)]$. Hence, we have $$\theta_{PO-RDID}^{\mathcal{L}_{2}}=\theta_{OLS}-\mathbb E_{\mu_{I_0}}[SB(I_0)].$$

$L\infty$ loss: Maximal regret

Note that

align*[align* omitted — 213 chars of source]

This optimization problem with the $L\infty$ loss is equivalent to a minimax criterion, and yields the mid-point of the bounds on $SB_1$ stated in Assumption (ref). Hence, the PO-RDID estimand is given by $$\theta_{PO-RDID}^{\mathcal{L}_{\infty}}=\theta_{OLS}-\frac{1}{2}\bigg(\inf_{\iota_0 \in \mathcal I_0}SB(\iota_0)+\sup_{\iota_0 \in \mathcal I_0}SB(\iota_0)\bigg).$$

Note that in all cases, if the information in the baseline period is a singleton, then the optimal $SB_1^{opt}$ is the selection bias in the baseline period $SB_0$, which is equivalent to the parallel trends assumption. Unlike the PO-RDID estimand obtained from $L1$ and $L2$ loss functions, that obtained from the $L\infty$ loss function does not require the knowledge of the distribution $\mu_{I_0}$ of the information $I_0$ but only its support and is easy to compute. However, when the distribution of $SB(I_0)$ is uniform over $\left[ \inf_{\iota_0\in \mathcal I_0} SB(\iota_0), \sup_{\iota_0 \in \mathcal I_0} SB(\iota_0)\right]$, then the optimal selection bias $SB_1^{opt}$ is the same in all three cases.

Forecasting $SB_1$ when the baseline information is ordered

Suppose that the baseline information set $\mathcal I_0$ is ordered (e.g., a set of multiple pre-treatment periods $\mathcal I_0=\{-T_0, -T_0+1, \ldots, -1, 0\}$ or a continuous baseline covariate $X_0$ like age). We can regress $SB(I_0)$ on $\{I_0, I_0^2, \ldots\}$ and use this regression to predict $SB_1$.

For instance, if $\mathcal I_0=\{-T_0, -T_0+1, \ldots, -1, 0\}$ and the selection bias $SB(I_0)$ is increasing over time, Assumption (ref) may not hold. The researcher could instead use this increasing trend information about the selection bias to forcast the next period selection bias $\widehat{SB}_1$.

The command rdid with an option rdidtype(2) can be used for a linear predicition model, and an option peval(\#) can be used to specify the evaluation point of the prediction.

Identification with multiple periods

We generalize our analysis to a setting where the treatment receipt occurs at multiple periods. We consider the following multiple treatment periods potential outcome model:

eqnarray[eqnarray omitted — 155 chars of source]

for $t\in \mathcal T_0 \cup \{1,\ldots,T\}$, where $Y_t$ denotes the observed outcome in period $t$, $D_t$ is the observed treatment status in period $t$ with $D_0=0$ by definition, while $Y_t(0,d_1, \ldots, d_T)$ is the potential outcome when the treatment path $(D_0,D_1, \ldots, D_T)$ is externally set to $(0,d_1,\ldots, d_T)$.\footnote{See Robins1986, Robins1987, and Han2021 for a similar definition of the potential outcome model.} Under the no-anticipation assumption, we have $Y_{t_0}(0,d_1,\ldots, d_T) =Y_{t_0}(0)$ for all $(d_1,\ldots,d_T) \in \{0,1\}^T$ and $t_0 \in \mathcal T_0$. That is, we assume that individuals do not anticipate any effects of the treatment before it occurs for the first time in $t=1$. However, we allow the individuals to anticipate the effects of the treatment for the rest of the period. This assumption is less restrictive than the commonly used no-anticipatory effects assumption.

In the above framework, the parameter of interest is the average treatment effect in period $t$ on the treated group following the path $\bm d^0 \equiv (0,d_1^0,\ldots, d_T^0)$ to $\bm d^1 \equiv (0,d_1^1,\ldots, d_T^1)$, which is defined as:

eqnarray*[eqnarray* omitted — 166 chars of source]

Similarly to what we have in the one post-treatment setting, we can write the difference-in-means estimand $(\theta_{DIM}^t)$ between the two groups $\bm d^0$ and $\bm d^1$ in period $t$ as:

align*[align* omitted — 263 chars of source]

where $SB_t[\bm d^0\rightarrow \bm d^1]\equiv \mathbb E\left[Y_t(\bm d^0) \vert (D_0,D_1, \ldots, D_T)=\bm d^1\right]-\mathbb E\left[Y_t(\bm d^0) \vert (D_0,D_1, \ldots, D_T)=\bm d^0\right]$. We extend Assumption (ref) to the current setting.

assumption[Extended bias set stability] For each $t$, \begin{align*} SB_t[\bm d^0\rightarrow \bm d^1] &\in \left[ \inf_{t_0\in \mathcal T_0} SB_{t_0}[\bm d^0\rightarrow \bm d^1], \sup_{t_0 \in \mathcal T_0} SB_{t_0}[\bm d^0\rightarrow \bm d^1]\right] \\ &\equiv \Delta_{SB(\mathcal T_0)}[\bm d^0\rightarrow \bm d^1], \end{align*} where $SB_{t_0}[\bm d^0\rightarrow \bm d^1]\equiv \mathbb E[Y_{t_0}(0)\vert (D_0,D_1, \ldots, D_T)=\bm d^1] - \mathbb E[Y_{t_0}(0)\vert (D_0,D_1, \ldots, D_T)=\bm d^0]$ is the baseline selection bias with respect to the treatment status in period $t_0 \in \mathcal T_0$ when the information set $\mathcal I_0$ is equal to $\mathcal T_0$.

Assumption (ref) is a generalization of Assumption (ref), and we have the following identification results for $ATT_t$ with different treatment paths.

propositionSuppose that model ((ref)) along with Assumption (ref) holds. Then, the following bounds are valid for $ATT_t$: \begin{align*} &ATT_t[\bm d^0 \rightarrow \bm d^1] \\ &\qquad \in \left[ \theta_{DIM}^t[\bm d^0\rightarrow \bm d^1] - \sup_{t_0 \in \mathcal T_0} SB_{t_0}[\bm d^0\rightarrow \bm d^1], \theta_{DIM}^t[\bm d^0\rightarrow \bm d^1]- \inf_{t_0 \in \mathcal T_0} SB_{t_0}[\bm d^0\rightarrow \bm d^1] \right], \\ & \qquad \equiv \Theta_I^t[\bm d^0 \rightarrow \bm d^1]. \end{align*} These bounds are sharp, and $\Theta_I^t[\bm d^0 \rightarrow \bm d^1]$ is the identified set for $ATT_t[\bm d^0 \rightarrow \bm d^1]$.

In a staggered design setting, the average treatment effect in period $t$ on units who are treated for the first time in period $g$ could be an interesting parameter, as considered in CallawaySantAnna2021: $$ATT(g, t) \equiv ATT_t[(0,\ldots,0,d_g=0,0,\ldots, 0)\rightarrow (0,\ldots,0,d_g=1,1,\ldots, 1)].$$

By defining $G$ as $$ G = \left\{

array[array omitted — 178 chars of source]

\right. , $$ where $G = \infty$ denotes the never-treated cohort as our control group, we can write

align*[align* omitted — 398 chars of source]

for $t \in \mathcal T_0 \cup \{1, \ldots, T \}$, by abusing the notations of $g$ and $\infty$ for paths.

Hence, Assumption (ref) can be re-written as

align*[align* omitted — 237 chars of source]

for each $t \in \{1, \ldots, T \}$, and $SB_{t_0}[\infty \rightarrow g] \equiv \mathbb E[Y_{t_0}(0)\vert G = g] - \mathbb E[Y_{t_0}(0)\vert G = \infty ]$ for $t \in \mathcal T_0$. Thus, we have

align*[align* omitted — 263 chars of source]

The rdid command

Syntax

The rdid command estimates bounds on the ATT under the assumptions discussed above. The syntax for the rdid command is

stsyntaxrdid \depvar\ \optindepvars\ \optif\ \optin , \underbar{treat}name(\varname) \underbar{post}name(\varname) \underbar{info}name(\varname) \optional{ rdidtype(\num) peval(\num) \underbar{lev}el(\num) \underbar{fig}ure \underbar{b}rep(\num) \underbar{cl}ustername(\varname)}

Options

\hangpara {\tt treatname(\varname)} specifies the variable name of the treatment indicator. treatname() is required.

\hangpara {\tt postname(\varname)} specifies the variable name of the post-treatment period indicator. postname() is required.

\hangpara {\tt infoname(\varname)} specifies the variable name of the information index. infoname() is required.

\hangpara {\tt rdidtype(\num)} specifies the type of RDID estimator; 0 for the simple RDID, 1 for the policy-oriented (PO) RDID, and 2 for the RDID with linear predictions. The default is {\tt rdidtype(0)}.

\hangpara {\tt level(\num)} specifies the confidence level, as a percentage, for confidence intervals (CIs). The default is {\tt level(95)}. Three different CIs are provided with {\tt rdidtype(0)}; 1) CI for the bounds ye2023negative, 2) CI for the ATT ye2023negative, and 3) CI for the bounds (union bounds). Bootstrapped CIs are provided with {\tt rdidtype(1)} and {\tt rdidtype(2)}.

\hangpara When {\tt figure(\ststring)} is specified, rdid saves a scatter plot of estimated selection biases and the information index in the pre-treatment periods as {\tt figure(\ststring).png}.

\hangpara {\tt brep(\num)} specifies the number of bootstrap replicates. The default is {\tt brep(500)}.

\hangpara {\tt clustername(\varname)} specifies the variable name that identifies resampling clusters in bootstrapping.

Stored results

rdid stores the following in e(). When indepvars is specified, $\tau^{DR}$ is estimated using post-treatment periods (see bk2023rdid, bk2023rdid) and stored in {\tt e(DR)}. Otherwise, $\theta_{OLS}$ is estimated and stored in {\tt e(OLS)}. Also, stored scalars depend on {\tt rdidtype(\num)} specification (referred as `type' below). In the following, LB (UB) stands for lower bound (upper bound).

stresults\stresultsgroup{Scalars (common)} \\ \stcmd{e(DR)} & $\widehat{\tau}^{DR}$ estimate & \stcmd{e(OLS)} & $\widehat{\theta}_{OLS}$ estimate \\ \stcmd{e(N)} & number of observations & & \\ \stresultsgroup{Scalars (type=0)} \\ \stcmd{e(SB_LB)} & LB of $\widehat{\Delta}_{SB(\mathcal I_0)}$ & \stcmd{e(SB_UB)} & UB of $\widehat{\Delta}_{SB(\mathcal I_0)}$ \\ \stcmd{e(RDID_LB)} & LB of $\widehat{\Theta}_I$ & \stcmd{e(RDID_UB)} & UB of $\widehat{\Theta}_I$ \\ \stcmd{e(CI1_LB)} & LB of CI-1 & \stcmd{e(CI1_UB)} & UB of CI-1 \\ \stcmd{e(CI2_LB)} & LB of CI-2 & \stcmd{e(CI2_UB)} & UB of CI-2 \\ \stcmd{e(CI3_LB)} & LB of CI-3 & \stcmd{e(CI3_UB)} & UB of CI-3 \\ \stresultsgroup{Scalars (type=1)} \\ \stcmd{e(L1_CI_LB)} & LB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_1}$ & \stcmd{e(L1_CI_LB)} & UB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_1}$ \\ \stcmd{e(L2_CI_LB)} & LB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_2}$ & \stcmd{e(L2_CI_LB)} & UB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_2}$ \\ \stcmd{e(Linf_CI_LB)} & LB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_{\infty}}$ & \stcmd{e(Linf_CI_LB)} & UB of CI for $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_{\infty}}$ \\ \stcmd{e(L1_PE)} & $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_1}$ estimate & \stcmd{e(L2_PE)} & $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_2}$ estimate \\ \stcmd{e(Linf_PE)} & $\widehat{\theta}_{PO-RDID}^{\mathcal{L}_{\infty}}$ estimate & & \\ \stresultsgroup{Scalars (type=2)} \\ \stcmd{e(SB_hat)} & $\widehat{SB}_1$ estimate & \stcmd{e(proj_PE)} & $\widehat{\theta}_{RDID}$ estimate \\ \stcmd{e(CI_LB)} & LB of CI for $\widehat{\theta}_{RDID}$ & \stcmd{e(CI_UB)} & UB of CI for $\widehat{\theta}_{RDID}$ \\ \stresultsgroup{Macros} \\ \stcmd{e(cmd)} & \stcmd{rdid} & \stcmd{e(indep)} & name of independent variable \\ \stcmd{e(depvar)} & name of dependent variable(s) & \stcmd{e(clustvar)} & name of cluster variable \\ \stcmd{e(level)} & confidence level & & \\ \stresultsgroup{Matrices} \\ \stcmd{e(results)} & displayed summary matrix & & \\

The rdid_dy command

Syntax

The rdid_dy command is a wrapper command that implements RDID estimation separately for each levels of a time variable in the post-treatment period and collects the results. The syntax for the rdid_dy command is

stsyntaxrdid_dy \depvar\ \optindepvars\ \optif\ \optin , \underbar{treat}name(\varname) \underbar{post}name(\varname) \underbar{info}name(\varname) \underbar{t}name(\varname) \optional{ rdidtype(\num) peval(\num) \underbar{lev}el(\num) \underbar{fig}ure(\ststring) \underbar{b}rep(\num) \underbar{cl}ustername(\varname) \underbar{ci}type(\num) \underbar{l}osstype(\num)}

Options

\hangpara {\tt tname(\varname)} specifies the variable name for the time. tname() is required.

\hangpara {\tt peval(\num)} specifies the evaluation point for the RDID with linear predictions. The default is the mean value of {\tt infoname(\varname)} in the post-treatment period. In particular, if {\tt infoname(\varname)} is the same as {\tt tname(\varname)}, each RDID estimate is obtained using each level of {\tt tname(\varname)} as the evaluation point. This option is ignored when {\tt rdidtype(\num)} is not equal to 2.

\hangpara When {\tt figure(\ststring)} is specified, rdid_dy saves line plots of RDID estimates and confidence intervals over the post-treatment periods as {\tt figure(\ststring).png}.

\hangpara {\tt citype(\varname)} specifies the type of confidence intervals to collect for {\tt rdidtype(0)}; 1 yields CI for the RDID bounds ye2023negative, 2 yields CI for the ATT ye2023negative, and 3 yields CI for the RDID bounds (Union Bounds). The default is {\tt citype(1)}. This option is ignored when {\tt rdidtype(\num)} is not equal to 0.

\hangpara {\tt losstype(\num)} specifies the type of loss function for {\tt rdidtype(1)} (PO-RDID); 1 for $\mathcal{L}_1$, 2 for $\mathcal{L}_2$, and 0 for $\mathcal{L}_{\infty}$. The default is {\tt losstype(1)}. This option is ignored when {\tt rdidtype(\num)} is not equal to 1.

Stored results

rdid_dy stores the following in e(). Stored scalars depend on {\tt rdidtype(\num)} specification; In the following, LB (UB) stands for lower bound (upper bound) and (t) denotes each level of {\tt tname(\varname)}.

stresults\stresultsgroup{Scalars} \\ \stcmd{e(N)} & \qquad number of observations & \stcmd{e(RDID_PE_(t))} & \qquad $\widehat{\theta}_{RDID}$ estimate for (t) \\ \stcmd{e(RDID_LB_(t))} & \qquad LB of $\widehat{\Theta}_{I}$ for (t) & \stcmd{e(RDID_UB_(t))} & \qquad UB of $\widehat{\Theta}_{I}$ for (t) \\ \stcmd{e(RDID_LB_(t))} & \qquad LB of CI for (t) & \stcmd{e(RDID_UB_(t))} & \qquad LB of CI for (t) \\ \stresultsgroup{Macros} \\ \stcmd{e(cmd)} & \qquad \stcmd{rdid_dy} & \stcmd{e(indep)} & \qquad name of independent \\ &&& \qquad variable \\ \stcmd{e(depvar)} & \qquad name of dependent & \stcmd{e(clustvar)} & \qquad name of cluster variable \\ & \qquad variable(s)&& \\ \stcmd{e(level)} & \qquad confidence level & & \\ \stresultsgroup{Matrices} \\ \stcmd{e(results)} & \qquad displayed summary matrix & & \\

The rdidstag command

Syntax

The rdidstag command estimates bounds on the ATT over time for different cohorts in the staggered adoption design under the assumptions discussed above. The syntax for the rdidstag command is

stsyntaxrdidstag \depvar\ \optindepvars\ \optif\ \optin , \underbar{g}name(\varname) \underbar{t}name(\varname) \optional{ \underbar{post}name(\varname) \underbar{info}name(\varname) \underbar{lev}el(\num) \underbar{fig}ure \underbar{cl}ustername(\varname)}

Options

\hangpara {\tt gname(\varname)} specifies the variable name of the cohort index, $g \in \mathcal{G}$, where $g=0$ is assigned for the never-treated group. gname() is required.

\hangpara {\tt tname(\varname)} specifies the variable name of time index, $t$. tname() is required.

\hangpara {\tt postname(\varname)} specifies the variable name of the post-treatment period indicator. The default is using units with $t \geq \min(\mathcal{G} \setminus \{0\})$ as post-period units.

\hangpara {\tt infoname(\varname)} specifies the variable name of the information index. The default is using pre-treatment periods; i.e., $\mathcal{I}_0 = \{t \mid t < \min(\mathcal{G} \setminus \{0\}) \}$.

\hangpara {\tt level(\num)} specifies the confidence level, as a percentage, for confidence intervals. The default is {\tt level(95)}.

\hangpara When {\tt figure} is specified, rdidstag saves line plots of RDID bounds and confidence intervals over the post-treatment periods for each $g \in \mathcal{G}$ as {\tt figure(\ststring).png}.

\hangpara {\tt clustername(\varname)} specifies the variable name that identifies resampling clusters in bootstrapping.

Stored results

rdidstag stores the following in e(). In the following, LB (UB) stands for lower bound (upper bound). Also, (g) and (t) denote each level of {\tt gname(\varname)} and {\tt tname(\varname)}, respectively.

stresults\stresultsgroup{Scalars} \\ \stcmd{e(RDID_LB_(g)_(t))} & \qquad LB of $\widehat{\Theta}(g, t)$ & \stcmd{e(RDID_UB_(g)_(t))} & \qquad UB of $\widehat{\Theta}(g, t)$ \\ \stcmd{e(CI_LB_(g)_(t))} & \qquad LB of CI for $\widehat{\Theta}(g, t)$ & \stcmd{e(CI_UB_(g)_(t))} & \qquad UB of CI for $\widehat{\Theta}(g, t)$ \\ \stcmd{e(N)} & \qquad number of observations & & \\ \stresultsgroup{Macros} \\ \stcmd{e(cmd)} & \qquad \stcmd{rdidstag} & \stcmd{e(indep)} & \qquad name of independent \\ &&& \qquad variable \\ \stcmd{e(depvar)} & \qquad name of dependent & \stcmd{e(clustvar)} & \qquad name of cluster variable \\ & \qquad variable(s)&& \\ \stcmd{e(level)} & \qquad confidence level & & \\ \stresultsgroup{Matrices} \\ \stcmd{e(results)} & \qquad displayed summary matrix & & \\

Applications

We illustrate usages of rdid, rdid_dy, and rdidstag commands using two applications. We first demonstrate rdid and rdid_dy commands using an application of cai2016 in Section (ref) and then rdidstag command using data from dube2016 in Section (ref).

Effects of insurance provision on tobacco production

cai2016 investigates the impact of insurance provision on tobacco production using a household-level panel dataset provided by the Rural Credit Cooperative (RCC), the main rural bank in China.

Using rdid command

In the following runs, the dependent variable is area_tob and the covariates are hhsize, educ_scale, and age, yielding the doubly-robust estimate bk2023rdid. Also, treatment variable and the post-treatment period indicator is is treatment and policy2, respectively, and the information set is specified as pre-treatment periods by using \texttt{year}. We use 1,000 bootstrap replicates and cluster the standard errors using \texttt{hhno}. First, the following is yielding the RDID bounds with the three different confidence intervals.

stlog. use "analysis.dta" . drop if sector != 1 (20,519 observations deleted) . rdid area_tob hhsize educ_scale age, treat(treatment) post(policy2) /// > info(year) b(1000) cl(hhno) **** RDID Estimation version 1.8 **** Y name: area_tob D name: treatment X name(s): hhsize educ_scale age \HLI{16}{\TOPT}\HLI{19} Y: area_tob {\VBAR} LB UB \HLI{16}{\PLUS}\HLI{19} RDID {\VBAR} 0.8086 0.9479 CI_1 {\VBAR} 0.6203 1.1257 CI_2 {\VBAR} 0.6235 1.1222 CI_3 {\VBAR} 0.6187 1.1257 \HLI{16}{\BOTT}\HLI{19} * RDID: Point estimates for RDID bounds * CI_1: Confidence interval for the bounds (Ye et al.) * CI_2: Confidence interval for the ATT (Ye et al.) * CI_3: Confidence interval for the bounds (Union Bounds) \nullskip

Second, by specifying the option {\tt rdidtype(1)}, we can have PO-RDID estimates and their confidence intervals under the three loss functions.

stlog. rdid area_tob hhsize educ_scale age, treat(treatment) post(policy2) /// > info(year) rdidtype(1) b(1000) cl(hhno) **** RDID Estimation version 1.8 **** Y name: area_tob D name: treatment X name(s): hhsize educ_scale age \HLI{16}{\TOPT}\HLI{23} Y: area_tob {\VBAR} PE CI_LB CI_UB \HLI{16}{\PLUS}\HLI{23} L1 {\VBAR} 0.8770 0.6732 1.0806 L2 {\VBAR} 0.8778 0.6911 1.0753 Linf {\VBAR} 0.8783 0.6978 1.0710 \HLI{16}{\BOTT}\HLI{23} \nullskip

Lastly, by specifying the option {\tt rdidtype(2)}, we can have an RDID estimate based on a linear prediction with an ordered information set, pre-treatment periods. The predicted value is evaluated at the mean of the post-treatment periods, 2005.5, with the option {\tt peval(2005.5)}.

stlog. rdid area_tob hhsize educ_scale age, treat(treatment) post(policy2) /// > info(year) rdidtype(2) peval(2005.5) b(1000) cl(hhno) **** RDID Estimation version 1.8 **** Y name: area_tob D name: treatment X name(s): hhsize educ_scale age \HLI{16}{\TOPT}\HLI{39} {\VBAR} PE CI_LB CI_UB TAU_DR SB_hat \HLI{16}{\PLUS}\HLI{39} Y {\VBAR} area_tob {\VBAR} 0.7239 0.4760 0.9725 1.9154 1.1914 \HLI{16}{\BOTT}\HLI{39} \nullskip

Using rdid_dy command

Recall that rdid_dy command is a wrapper command that implements RDID estimation separately for each levels of a time variable in the post-treatment period and collects the results. In the following runs, the same options are specified as before, but the option {\tt tname(year)} is used to denote the variable for the post-treatment periods. Also, the following run produces Figure (ref) that summarizes the output.

stlog. use "analysis.dta" . drop if sector != 1 (20,519 observations deleted) . rdid_dy area_tob hhsize educ_scale age, treat(treatment) post(policy2) /// > info(year) tname(year) fig b(1000) cl(hhno) file{\bftt{ rdid_dy.png}} saved as PNG format \HLI{16}{\TOPT}\HLI{39} T: year {\VBAR} RDID_LB RDID_UB CI_LB CI_UB \HLI{16}{\PLUS}\HLI{39} 2003 {\VBAR} 0.1829 0.3221 0.0239 0.4855 2004 {\VBAR} 0.6181 0.7573 0.3988 0.9545 2005 {\VBAR} 0.8452 0.9844 0.6401 1.1692 2006 {\VBAR} 0.7897 0.9289 0.5738 1.1377 2007 {\VBAR} 1.1269 1.2661 0.8829 1.4996 2008 {\VBAR} 1.2372 1.3764 0.9328 1.6571 \HLI{16}{\BOTT}\HLI{39} * Confidence intervals are obtained for the bounds (Ye et al.) \nullskip

Note that rdid_dy command can be used for the option {\tt rdidtype(1)} or {\tt rdidtype(2)}, and the following code runs for collecting PO-RDID with $\mathcal{L}_1$ (default), producing Figure (ref) that summarizes the output.

stlog. rdid_dy area_tob hhsize educ_scale age, rdidtype(1) treat(treatment) /// > post(policy2) info(year) tname(year) fig b(1000) cl(hhno) file{\bftt{ rdid_dy.png}} saved as PNG format \HLI{16}{\TOPT}\HLI{29} T: year {\VBAR} RDID_PE CI_LB CI_UB \HLI{16}{\PLUS}\HLI{29} 2003 {\VBAR} 0.2513 0.0605 0.4273 2004 {\VBAR} 0.6865 0.4497 0.8893 2005 {\VBAR} 0.9136 0.6951 1.1206 2006 {\VBAR} 0.8581 0.6253 1.0898 2007 {\VBAR} 1.1953 0.9517 1.4561 2008 {\VBAR} 1.3056 0.9999 1.5860 \HLI{16}{\BOTT}\HLI{29} * RDID estimates are obtained as PO-RDID estimates (L1) \nullskip

figure[figure omitted — 177 chars of source]

Effects of minimum wage increases on teen employment

Following CallawaySantAnna2021 and dube2016, we used the Quarterly Workforce Indicators (QWI) data to collect the first quarter teen employment as our outcome variable. CallawaySantAnna2021 considered 7 years of periods between 2001 and 2007 where the federal minimum wage did not change over time, and 3 different control groups of $g=2004, 2006,$ and $2007$ with states that raised their minimum wage in or right before the beginning of years 2004, 2006, and 2007, respectively. The specific timing of the raise can be found in CallawaySantAnna2021, and it should be noted that there is some heterogeneity in the size of the minimum wage increase within each group. The control group consists of states that did not raise their minimum wage during this period.

The following code specifies the option {\tt gname(G)} for the cohort index variable, which has four levels: 0, 2004, 2006, and 2007. Note that because the options {\tt postname(\varname)} and {\tt infoname(\varname)} are not specified, all periods from 2004 onward are treated as post-treatment periods, and all periods before 2004 are treated as informational elements.

stlog. use "qwi_minwageclean2.dta" . destring years, replace years: all characters numeric; replaced as int . rdidstag lnemp_0A01_BS lnpop lnteenpop lnemp_0A00_BS, gname(G) tname(years) cl(countyfips) fig **** RDID Estimation version 1.7 **** Y name: lnemp_0A01_BS G name: G X name(s): lnpop lnteenpop lnemp_0A00_BS postname is not specified: using units with T >= min(G) as post-period units. infoname is not specified: using pre-period T < min(G) as the information set. \HLI{16}{\TOPT}\HLI{39} ATT(G/T) {\VBAR} RDID_LB RDID_UB 95CI_LB 95CI_UB \HLI{16}{\PLUS}\HLI{39} ATT(2004/2004) {\VBAR} -0.0250 -0.0183 -0.0646 0.0164 ATT(2004/2005) {\VBAR} -0.0551 -0.0484 -0.0973 -0.0111 ATT(2004/2006) {\VBAR} -0.1086 -0.1018 -0.1504 -0.0639 ATT(2004/2007) {\VBAR} -0.1153 -0.1086 -0.1617 -0.0653 ATT(2006/2004) {\VBAR} 0.0186 0.0578 -0.0117 0.0858 ATT(2006/2005) {\VBAR} 0.0408 0.0800 0.0077 0.1121 ATT(2006/2006) {\VBAR} 0.0230 0.0622 -0.0130 0.0974 ATT(2006/2007) {\VBAR} -0.0227 0.0165 -0.0585 0.0520 ATT(2007/2004) {\VBAR} 0.0065 0.0422 -0.0301 0.0738 ATT(2007/2005) {\VBAR} 0.0013 0.0370 -0.0329 0.0676 ATT(2007/2006) {\VBAR} -0.0346 0.0012 -0.0692 0.0383 ATT(2007/2007) {\VBAR} -0.0512 -0.0154 -0.0846 0.0207 \HLI{16}{\BOTT}\HLI{39} - Information elements: 2001 2002 2003 - Post-periods: 2004 2005 2006 2007 - Groups: 2004 2006 2007 file{\bftt{ rdidstag_g2004.png}} saved as PNG format file{\bftt{ rdidstag_g2006.png}} saved as PNG format file{\bftt{ rdidstag_g2007.png}} saved as PNG format \nullskip

figure[figure omitted — 153 chars of source]

Monte-Carlo Simulations

Modified Example 2: Ashenfelter1978's (Ashenfelter1978) dip

Consider

eqnarray*[eqnarray* omitted — 272 chars of source]

for $t \in \mathcal T_0 \cup \{1\}$ where $\mathcal I_0=\mathcal T_0=\{-2,-1,0\}$.

In this model, $SB_t=(1+\vert t \vert + t^2)(\alpha_1-\alpha_0)$ where $\alpha_1=\frac{\phi(1)}{1-\Phi(1)}\approx 1.53$ and $\alpha_0=-\frac{\phi(1)}{\Phi(1)} \approx -0.29$. We have $SB_0= \alpha_1-\alpha_0 \neq 3(\alpha_1-\alpha_0) =SB_{-1} \neq SB_{-2}=7(\alpha_1-\alpha_0)$ and $SB_0=\alpha_1-\alpha_0 \neq 3(\alpha_1-\alpha_0)=SB_1$. So, the standard parallel trends assumption does not hold. However, the selection bias $SB_1$ in period 1 belongs to the convex hull of all selection biases in period 0, i.e., $SB_1 \in [\min\{SB_0,SB_{-1},SB_{-2}\}, \max\{SB_0,SB_{-1},SB_{-2}\}]=[\alpha_1-\alpha_0,7 (\alpha_1-\alpha_0)]$. Hence, our identifying assumption holds.

We have $\theta_{OLS}=3(\alpha_1-\alpha_0)+\theta$, and $\Theta_I=[\theta-4(\alpha_1-\alpha_0),\theta+2(\alpha_1-\alpha_0)]$. The true $ATT=\theta$, and the DID estimand is $\theta_{DID}=\theta_{OLS}-SB_0=\theta+2(\alpha_1-\alpha_0)$. Thus, the DID estimand is upward biased and the bias is equal to $2(\alpha_1-\alpha_0)$. The bounds are generally informative about the magnitude of the average treatment effect on the treated (ATT) and can identify the sign of the ATT in some circumstances. For example, when the true ATT is equal to $-5$ or $9$, our bounds as well as the standard DID estimand identify the correct sign. On the other hand, when the true ATT is equal to $-1$, our bounds do not identify any sign, as they contain zero. But, the standard DID estimand identifies a wrong sign, it shows that the ATT is positive while it is actually negative.

The following table shows a series of Monte-Carlo simulations for this DGP for three different confidence intervals (CIs): (1) ye2023negative's (ye2023negative) bootstrap CI for $\Theta_I$, (2) ye2023negative's (ye2023negative) bootstrap CI for ATT, (3) union bounds for $\Theta_I$.. The number of bootstrap replicates is 300, and the number of simulations is 1,000 for each case of $N \in \{200, 500, 1000, 5000, 10000, 50000\}$.

table[table omitted — 882 chars of source]

Example 1 from the article: with a covariate

Consider

eqnarray*[eqnarray* omitted — 189 chars of source]

where $X \rotatebox[origin=c]{90}{$\models$} U$ and $\mathcal I_0=\mathcal X=\{0,1 \}$.

We have $SB_0(x)=(1+x)(\alpha_1-\alpha_0)$, and $SB_1(x)=(1+0.5x)(\alpha_1-\alpha_0)$ where $\alpha_1=\frac{\phi(1)}{1-\Phi(1)}\approx 1.53$ and $\alpha_0=-\frac{\phi(1)}{\Phi(1)} \approx -0.29$. We have $SB_0(x) \in [\alpha_1-\alpha_0, 2(\alpha_1-\alpha_0)]$ and $SB_1(x) \in \left[\alpha_1-\alpha_0, 1.5(\alpha_1-\alpha_0)\right] \subseteq [\alpha_1-\alpha_0, 2(\alpha_1-\alpha_0)]\equiv \Delta_{SB(\mathcal X)}$. So, the standard parallel trends assumption does not hold as $SB_0(x)\neq SB_1(x)$. However, the selection bias $SB_1(x) $ in period 1 belongs to the convex hull of all selection biases in period 0, i.e., $SB_1(x) \in \Delta_{SB(\mathcal X)}$. Hence, our identifying assumption holds. We have $\theta_{OLS}(x)=(1+ 0.5 x)(\alpha_1-\alpha_0)+\theta x$, and our new bounds $ \Theta_I $ are obtained as $ATT(x) \in [\theta x - (1-0.5 x)(\alpha_1-\alpha_0),\theta x + 0.5 x(\alpha_1-\alpha_0)]$. The actual conditional ATT function is $ATT(x)=\theta x$, but the standard conditional DID estimand is $\theta_{DID}(x)=\theta x - 0.5 x (\alpha_1 - \alpha_0)$.

Moreover, using the distribution of $X$, we have an identified set for ATT $$ \Theta_I = \left[ \left(\frac{1}{2}p - 1 \right)(\alpha_1-\alpha_0) + p\theta, \frac{1}{2}(\alpha_1-\alpha_0) + p\theta \right] $$ where $ATT = p\theta$. The simulation is implemented using $p=0.5$, and the results are presented in the following table.

table[table omitted — 864 chars of source]

The modified staggered adoption design

Consider

eqnarray*[eqnarray* omitted — 221 chars of source]

where $U_i \rotatebox[origin=c]{90}{$\models$} (\{\varepsilon_{it}, \theta_{it} \}_{t=1}^T )$, $\theta_{it} \sim \mathcal{N}(\frac{1+t^2}{2}, 1)$, $\mathcal{I}_0 = \mathcal T_0 = \{-T, \ldots, 0 \}$, $\varepsilon_{it} \sim \mathcal{N}(t^2, 1)$, $U_i \sim \mathcal{U}_{[0, 2]}$, and $T=4$.

Note that the PT is violated as

align*[align* omitted — 215 chars of source]

for $t > 0$ when $g \neq g'$, especially with $g' = \infty$ for the never-treated cohorts being the control units.

However, note that

align*[align* omitted — 207 chars of source]

and

align*[align* omitted — 250 chars of source]

for all $t=1, \ldots, T$.

In the case of $g' = \infty$, we have

align*[align* omitted — 156 chars of source]

and

align*[align* omitted — 173 chars of source]

Hence, we have

align*[align* omitted — 98 chars of source]

where $\mu(g) \equiv \mathbbm{E}[U_{i} \mid G_i = g] - \mathbbm{E}[U_{i} \mid G_i = \infty]$.

Under this setup, we implemented a series of Monte-Carlo simulations with various $N$ values, and the following tables summarize the coverage probability ($CP_{inf}$) and the average length of the 95% confidence interval for each case.

table[table omitted — 977 chars of source]

table[table omitted — 958 chars of source]

table[table omitted — 953 chars of source]

Conclusion

In this article, we introduce the rdid, rdid_dy, and rdidstag commands to implement the RDID framework by bk2023rdid for both the canonical $2 \times 2$ DID settings and staggered adoption designs. By allowing users to specify their information set as pre-treatment periods, the commands enhance robustness to potential violations of the parallel trends assumption in both settings. We anticipate that continued enhancements and community feedback will further increase its effectiveness and adoption in empirical research.