Doubly Robust Semiparametric Difference-in-Differences Estimators with High-Dimensional Data

Doubly Robust Semiparametric Difference-in-Differences Estimators with High-Dimensional Data

Introduction

This paper proposes a doubly robust two-stage semiparametric difference-in-difference estimator for estimating heterogeneous treatment effects conditional on high-dimensional covariates. The difference-in-difference (DiD) design has been widely adopted in policy evaluation from academia to industry when a real experiment is expensive or infeasible. When a policy/feature only affects a fraction of the population, DiD design can identify the average treatment effect on the treated (ATT) based on observational data. It is based on the simple idea of comparing the difference in pre and post-treatment outcome of those individuals who are affected and those who are not affected by the policy/feature of interest.

A key identification assumption for the classical DiD design is the parallel trend assumption. It requires that the outcome variables for treated and non-treated individuals would have followed parallel paths over time in the absence of treatment. However, this assumption ignores the potential selection problem due to individual heterogeneity. For example, a company might want to evaluate the effect of an email marketing campaign (advertisement through email with an embedded promo link). A researcher can compare the customers' conversion rate (whether a purchase was made) before and after the campaign for a group of treated customers (click into the link) and a group of non-treated customers (did not click into the link). If existing customers are more likely to click into the link and also more likely to purchase again even without the campaign intervention, the classical DiD estimator will lead to a positive bias and exaggerate the effect of the campaign. To account for such case, Abadie2005 proposed a two-stage semiparametric estimator with the so-called conditional parallel trend assumption. In this framework, a propensity score is estimated in the first stage to explicitly account for any observed confounders that may affect both the treatment take-up as well as the outcome growth trend.

While the semiparametric DiD (semi-DiD) estimator is comprehensively used by researchers in academia and industry, three major challenges arise in practice. First, the semi-DiD estimator becomes difficult to implement when there exist too many covariates. Following the previous example, researchers may also observe customers' browsing history and may suspect customers who visited certain (unknown) websites are more likely to click into the link while also more likely to make the purchase. However, the semi-DiD estimator cannot be implemented if the number of attributes exceeds the number of observations. Therefore, researchers may be forced to select covariates based on their intuition or insights and may lead to further biases (Belloni2017). Even when the number of observations is larger than the number of covariates, the semi-DiD estimator may still contain a large bias when too many covariates are included (Matias2019). Second, the semi-DiD estimator is sensitive to the choice of specification for propensity score estimation. This is a similar problem for the inverse propensity score weighted (IPW) estimator and it becomes a more subtle problem if machine learning methods (e.g. random forrest, neural network, etc.) are used to predict propensity score. Third, conditional or heterogeneous treatment effects on treated (ATT) is often of interest to practitioners. While semi-DiD framework provides a way to estimate conditional or heterogeneous ATT under a vector of low-dimensional covariates, it is not clear how this framework can be extended to the high-dimensional case, (e.g. how to develop estimation and inference methods and theory).

In this paper, we propose a new estimator to solve the above three problems. Our doubly robust DiD (Dr-DiD) estimator is robust to model miss-specifications under high-dimensional covariates. We show the desired rate of convergence of our estimator can be achieved as long as either the propensity score function or the outcome equation can be approximated asymptotically at a moderate rate. Thus, a general set of machine learning methods can be used in our framework. Although diff-in-diff design is an ATT estimator, we show that the semi-DiD estimator can be extended to an augmented inverse propensity score weighted (AIPW) estimator. We show that the extended AIPW form still preserves the doubly robustness property under the parallel trend assumption.

To further incorporate high-dimensional covariates and heterogeneous treatment effects, we consider a partially linear specification in the potential outcome estimation. The partially linear form is composed of a nonparametric specification from a set of low-dimensional covariates as well as a linear parametric specification from a set of high-dimensional covariates, which provides a flexible functional form to model the potential outcome. This is a very useful specification in real world application. For example, researchers may be interested in the nonlinear relationship between the outcome variable and a set of covariates while also facing a large number of indicator variables such as age, gender and region.

We derive the rate of convergence for our estimator as well as a de-bias procedure for inference. We show that the high-dimensional linear part of the estimator can achieve the oracle rate of convergence, while the nonparametric part maintains the nonparametric rate of convergence. With bias correction, the high-dimensional linear part can achieve normality at $\sqrt{n}$-rate, while the nonparametric part can achieve the normality at the nonparametric rate. Finally, we demonstrate the finite sample performance of our estimator in a simulation study and apply our estimator to study the effect of the Fair Minimum Wage Act on the unemployment rate using the data collect by CallawayLi2020. We show that the heterogeneity in the effect of this policy can be explained by variations in demographics. More specifically, counties with larger population and higher median income are more likely to suffer from an increase in the unemployment rate. These findings coincide with the canonical economic theory on unemployment rate. For example, a higher median income level implies a higher substitution cost for workers currently at minimum wage and thus leads to an increase in the unemployment rate when minimum wage rises. On the other hand, regions with larger population sizes have more labor supply and thus a minimum wage raise can also lead to a surplus.

In summary, the main contributions of this work are as follows: first, we propose a doubly robust approach to estimate heterogeneous ATT conditional on covariates for DiD models that allows either the propensity scores or the model for ATT to be misspecified. Second, we propose a regularized two-stage estimation procedure for DiD models that allows (i) suitable machine learing tools to estimate the first-stage propensity socres and (ii) high-dimensional covariates and nonparametric specification for the heterogeneous ATT in the second-stage. Third, we provide a novel approach to simultaneously correct the biases due to both stages and provide a novel statistical inference procedure based on the de-biased estimator. Finally, as a useful byproduct, we derive novel estimation and inference methods for a partially linear model for both the high-dimensional parametric parameter and the nonparametric function.

Related Literature

This paper is related to the vast literature on robust estimation and inference for treatment effects models; see for example, robins1994estimation, tan2006distributional, chen2008semiparametric, graham2012inverse, OkuiTan2012doubly, farrell2015robust, vermeulen2015bias, ogburn2015doubly, Belloni2017, lee2017doubly, Chernozhukov2016, sloczynskiWooldridge2018, kennedy2019robust and Tan2020 among many others. Our work is particularly closely related to a recent work independently developed by SantAnnaZhao2019. Both are based on the seminal framework proposed in Abadie2005. Our estimator complements theirs as we focus on estimation and inference for heterogeneous ATT conditional on covariates in a high-dimensional setting while SantAnnaZhao2019 focus on efficient estimation of ATT when the dimension of covariates is fixed and is much smaller than the sample size.

This paper also contributes to the literature by connecting the widely used DiD estimator with the machine learning/ high-dimensional statistic literature. The DiD estimator has been an active research field in the economic literature, e.g. Card1994, Abadie2005, Athey2006, Imai2019, Athey2019, Callaway2019 among others. Our paper proposes a specific DiD estimator so that high-dimensional/machine learning tools can be applied. This paper also contributes to a set of works that apply machine learning tools to casual inference. This includes chernozhukov2016locally, Belloni2017, semenova2017estimation, Chernozhukov2016, Syrgkanis2019, fan2020estimation, Tan2020, etc. Our work distinguish this literature in the following two ways. First, we propose a doubly robust diff-in-diff estimator in the high-dimensional/machine learning setting that has not been studied. Second, to our best knowledge, the doubly-robust estimators in these papers use various high-dimensional set of covariates and machine learning methods to deal with the selection into treatment. However, the ultimate parameter of interest in the second-stage is a low-dimensional subset of the covariates so traditional nonparametric estimator can apply. By contrast, our parameter of interest in the second-stage contains both high-dimensional covariates in the parametric part and an unknown function, which brings substantial challenges for estimation and inference. We construct a new Neyman orthogonal moment condition (chernozhukov2016locally) and propose de-biased estimators for both the parametric parameters and the nonparametric function in the second-stage to construct valid confidence intervals.

Moreover, as useful by-products, we provide an inference method for a partially linear model for both the parametric parameter and the nonparametric function when the linear part contains high-dimensional covariates. Thus, this work is related to recent discussion in muller2015partial, MaHuang2016, YLC2019, ZhuYuGuang2019high, among others. Our paper departs from the existing papers in the following three aspects. First, our partially linear form is in the second-stage outcome equation so estimation and inference results have to take the first-stage estimators into consideration, while the existing papers focus on a one-stage regression problem. Second, the above papers propose estimators with penalized estimation in functional space. As is pointed out in shen1997methods, this approach often leads to undesirable properties of the estimates, such as inconsistency and roughness. Moreover, such an optimization procedure is difficult to implement in practice. Therefore, we consider the extension to sieve estimation in our estimator by approximating the nonparametric function with sieves so that we carry out optimization within a dense subset of the infinite- dimensional space, which is finite-dimensional and therefore easy to work with. Finally, the existing literature is concerned with asymptotic theories and inference procedures for the parametric parameters only. The nonparametric function is profiled out as an infinite-dimensional nuisance parameter. This paper considers the joint asymptotic theory and inference methods when the parameters of interest are not only the parametric parameter, but also the nonparametric function. To the best of our knowledge, these results are new to the literature. We show that the parametric parameter converges to a normal distribution with a $\sqrt{n}$-rate. The parametric estimator achieves the semiparametric efficiency bound when the error term is homoskedastic, while the functional of a nonparametric function converges to a normal distribution with a nonparametric rate. We observe that the marginal asymptotic variance for the nonparametric component is, in general, different from those derived without the high-dimensional parametric parameter, i.e., Newey1997, BCCK2015 and ChenChristensen2015. This result may be of independent interest to the readers. \\

Organization of the Paper

The paper is organized as follows. The estimator is proposed in Section 2. Rate of convergence of the estimator and inference theory are developed in Sections 3 and 4, respectively. Section 5 presents extensive simulation results to evaluate the finite sample performance. An empirical study on the effect of the Fair minimum Wage Act on the unemployment rate is presented in Section 6. Section 7 concludes. We defer the proofs to the Appendices.

Notation

For a vector ${\bf x}=(x_{1},\dots,x_{d})^{\top}\in\mathbb{R}^{d}$ and $1\leq q\leq\infty$, let $\|{\bf x}\|_{q}=\left(\sum_{i=1}^{d}|x_{i}|^{q}\right)^{1/q}$, $\|{\bf x}\|_{\infty}=\max_{1\leq i\leq d}|x_{i}|$, $\|{\bf x}\|_{0}=\vert{\rm supp}({\bf x})\vert,$ where ${\rm supp}({\bf x})=\{j:x_{j}\neq0\}$ and $|a|$ is the cardinality of a set $a$. For a symmetric matrix $A$, let $\Lambda_{\max}(A)$ and $\Lambda_{\min}(A)$ be the maximum and minimum eigenvalues of $A$. For a matrix $B=[B_{jk}]$, let $\Vert B\|_{\max}=\max_{jk}\vert B_{jk}\vert,$ $\|B\|_{1}=\sum_{jk}\vert B_{jk}\vert$, $\|B\|_{2}=\sqrt{\Lambda_{\max}(B^{\top}B)}$ and $\|B\|_{\ell_{\infty}}=\max_{j}\sum_{k}|B_{jk}|$. For any function $f:\mathcal{Z}\rightarrow\mathbb{R},$ let $\|f\|_{\infty}=\sup_{z\in\mathcal{Z}}\vert f(z)\vert$, $\|f\|_{P,2}=\sqrt{\mathbb{E}f^{2}(z)}$ and $\|f\|_{n}=\sqrt{n^{-1}\sum_{i=1}^{n}f^{2}(Z_{i})}$. We denote $I_{d}$ as the $d\times d$ identity matrix. For a set $S\subseteq\{1,\dots,d\}$, let ${\bf x}_{S}=\{x_{j}:j\in S\}$ and $S^{c}$ be the complement of $S$. Let $S_0$ be the set of all non-zero components of $\beta_0$ and $s_0=|S_0|$. We use $\bigtriangledown_{S}f({\bf x})$ to denote the gradient of $f({\bf x})$ with respect to ${\bf x}_{S}$. Given $a,b\in\mathbb{R}$, let $a\vee b$ and $a\wedge b$ denote the maximum and minimum of $a$ and $b$. For two positive sequences $a_{n}$ and $b_{n}$, let $a_{n}\asymp b_{n}$ denote $C\leq a_{n}/b_{n}\leq C'$ for some $C,C'>0$; let $a_{n}\lesssim b_{n}$ denote $a_{n}\leq Cb_{n}$ for some constant $C>0$. Also, we write $a_{n}=O(b_{n})$ if $|a_{n}|\leq C|b_{n}|$. We use $X_{n}\rightarrow_{p}a$ for some constant $a$ if a sequence of random variables $X_{n}$ converges in probability to $a$. Similarly, if $X_{n}$ converges weakly to $X$ we write $X_{n}\rightsquigarrow X$ for some random variable $X$. For notational simplicity, we use $C$, $C'$ and $C''$ to denote generic constants, whose values can change from line to line. Let $\mathbb{E}_{n}f=\frac{1}{n}\sum_{i=1}^{n}f(X_i)$ and $\mathbb{G}_{n}f=\mathbb{E}_{n}f-\mathbb{E}f$.

A random variable $X$ is called sub-exponential if there exists some positive constant $K_{1}$ such that $\mathbb{P}(|X|>t)\leq\exp(1-t/K_{1})$ for all $t\geq0$. The sub-exponential norm of $X$ is defined as $\|X\|_{\psi_{1}}=\sup_{q\geq1}q^{-1}(\mathbb{E}|X|^{q})^{1/q}$. Similarly, a random variable $X$ is called sub-Gaussian if there exists some positive constant $K_{2}$ such that $\mathbb{P}(|X|>t)\leq\exp(1-t^{2}/K_{2}^{2})$ for all $t\geq0$. And the sub-Gaussian norm of $X$ is defined as $\|X\|_{\psi_{2}}=\sup_{q\geq1}q^{-1/2}\left(\mathbb{E}|X|^{q}\right)^{1/q}.$

Doubly Robust DiD Estimator

Denote $Y^0(i, t)$ as the potential outcome of individual $i$ at time $t$ being not treated and $Y^1(i, t)$ as the potential outcome of individual $i$ at time $t$ being treated. We cannot observe both $Y^0(i,t)$ and $Y^1(i,t)$ for the same individual, but we observe the realized outcome for individual $i$ at time $t$ as \[Y(i, t) = D_iY^1(i, t) + (1-D_i)Y^0(i, t)\] where $D_i$ is the treatment status at time $t=1$. For some observed covariates $W_i = (X_i, Z_i)$, we want to learn the heterogeneous treatment effect on the treated conditional on the covariates $W_i$ such that

Our parameter of interest is different from ATT (e.g. SantAnnaZhao2019), which is defined as

While it is useful to know ((ref)), a doubly robust estimator for ((ref)) conditional on the covariates could also be relevant and important in empirical applications when the parameter of interest is the heterogeneous treatment effects conditional on the covariates.

As pointed out in Abadie2005, the conventional DiD estimator is based on the strong assumption that outcomes for treated and non-treated groups or individuals would have followed parallel paths over time in the absence of treatment. That assumption can be easily violated when differences in observed characteristics create non-parallel outcome dynamics between treated and non-treated populations. Abadie2005 generalizes this assumption by allowing the parallel trend assumption to hold after conditioning on the covariates as follows:

In addition, a full support assumption will guarantee the existence of the propensity score function.

Together with Assumptions (ref) and (ref), the Abadie2005 estimand can be defined as

Defining $\Delta Y_{i} := Y(i, 1) - Y(i, 0)$, we then have

It is easy to see that Equation ((ref)) is in the form of Horvitz-Thompson estimator (Horvitz1952). As a natural extension to the IPW form estimator, we study whether a doubly robust form exists under the DiD setting and this leads to our parameters of interest as follows. Define \[ \Delta Y_{1i} :=Y^1(i, 1) - Y^0(i, 0), \qquad \Phi_1(W_i) := \mathbb{E}[\Delta Y_{1i}|W_i, D_i=1], \] \[\Delta Y_{0 i} : = Y^0(i, 1) - Y^0(i, 0), \qquad \Phi_0(W_i):= \mathbb{E}[\Delta Y_{0i}|W_i, D_i=0].\] Let \[\rho_0 = \frac{D_i - \mathbb{ P}(D_i=1|W_i)}{\mathbb{ P}(D_i=1|W_i)\mathbb{E}(D_i=0|W_i)},\] our doubly robust estimand is defined as

where $\mathbb{E}[D_i=0|W_i]$, $\mathbb{E}[D_i=1|W_i]$, $\Phi_0(W_i)$ and $\Phi_1(W_i)$ are nuisance functions to be estimated from the first-stage.

Lemma (ref) shows that, with Assumptions (ref) and (ref), we can have a doubly robust estimator for $\tau_0(W_i)$ when either the regression models $\Phi_0(\cdot)$ and $\Phi_1(\cdot)$ are misspecified or the propensity score $\mathbb{P}(D_i=1|X_i)$ is misspecified.

To model $\tau_0(\cdot)$, we consider a class of flexible high-dimensional partially linear model such that

where the linear part contains the parametric Euclidean vector $\beta_0\in\mathcal{B}\subseteq \mathbb{R}^{p}$ with $p>n$, and the nonparametric part contains an unknown function $f(\cdot): \mathcal{Z}\rightarrow \mathbb{R}$, where $\mathcal{Z}$ is a compact subset of $\mathbb{R}^{d_{z}}$. We will assume that the unknown function belongs to a smoothed function class defined in Section (ref).

Compare with the definition of equation (11) in Abadie2005, we define our estimand in equation (ref) in a partial linear form rather than approximate it with a best linear predictor. The semi-parametric structure is slightly stronger as equation (11) in Abadie2005 is satisfied if we plug in equation (ref) and allow $g(X_k, \theta)$ to admit a partial linear specification.

On the other hand, the partially linear specification in ((ref)) provides a flexible functional form while still allowing us to maintain the Neyman orthogonality condition when designing the estimator under the high dimensionally covariates. Theoretical properties of the semiparametric partially linear model when the dimension of $X$ is fixed and smaller than $n$ have been thoroughly discussed in the econometrics literature (engle1986semiparametric, Robinson1988root, ahn1993semiparametric, donald1994series, linton1995second, fan1999root, to mention only a few; see LiRacine for a review). We complement the literature by providing new estimation and inference methods and theory when $X$ is high-dimensional.

As a result of Lemma (ref) and equation (ref), if Assumptions (ref) and (ref) hold, we have

We are going to construct a two-step estimator of $(\beta_0,f_0)$ based on the sample analogue of ((ref)), where the first-step estimator estimates $\rho_0$ and the second-step estimator estimates $(\beta_0,f_0)$. We allow the propensity score, hence $\rho_0$, to be estimated by any suitable machine learning methods as long as certain conditions in Section (ref) are satisfied.

Estimation

Let $\hat{\pi}(\cdot)$, $\hat{\Phi}_1(\cdot)$ and $\hat{\Phi}_0(\cdot)$ be nonparametric or machine learning estimators of $\pi(W_i)$, $\Phi_1(W_i)$ and $\Phi_0(W_i)$, respectively. We propose the following two-stage estimator such that

where $f_n(\cdot)=\psi^{k_n}(\cdot)^\top\gamma_{n}$ is a sieve approximation of the unknown function $f(\cdot)\in\mathcal{F}$ with \[f_{0}(Z_i) = \sum_{j = 1}^{k_n} \psi_j(Z_i)\gamma_{j, n0} + r_{n}(Z_i):=f_{n0}(Z_i) +r_{ni}, \] where $r_{ni}:=r_{n}(Z_i),i=1,\dots,n$ is a sieve approximation error that depends on the smoothness of $f_{0}$ and the sample size $n$.

For $\alpha>0$ and any vector $\tau=(\tau_1,\dots,\tau_{d_z})$ of $d_z$ integers, define the differential operator $D^\tau=\partial^{\tau.}/\partial z^{\tau_1}_1...z^{\tau_{d_z}}_{d_z}$, where $\tau.=\sum_{l=1}^{d_{z}}\tau_l$. For a function $g:\mathcal{Z}\rightarrow \mathbb{R}$, let

Let $C^{\alpha}_M(\mathcal{Z})$ be the set of all continuous functions $g:\mathcal{Z}\rightarrow \mathbb{R}$ with $\|g\|_{\infty,\alpha}\leq M$. We assume that $\mathcal{F}\subseteq C^{\alpha}_M(\mathcal{Z})$. Let $\psi^{k_{n}}(Z_{i})=\left(\psi_{1}(Z_{i}),\dots,\psi_{k_{n}}(Z_{i})\right)^\top$ be a $k_n\times1$ vector of basis functions, and we use the notation $\mathcal{F}_{n}$ to represent the space of sieve functions. Define the projection of $X_{ij}, i=1,\dots,n$ onto $\mathcal{F}_{n}$ as \[\Pi_n(X_{ij}|Z)= \arg\min_{h^* \in \mathcal{F}_{n}}\|X_ij-h^*\|^2_{n} =\psi^{k_n}(Z)(\psi^{k_n}(Z)^\top\psi^{k_n}(Z))^{-1}\psi^{k_n}(Z)^\top X_ij, j=1,\dots,p\] and $\Pi_{n,X_i|Z_i}:=(\Pi_n(X_{i1}|Z_i),\dots, \Pi_n(X_{ip}|Z_i))^{\top}, i=1,\dots,n.$

Next, define \[\hat S_i := \hat{\rho}(W_i) \Big(\Delta Y_i -(1-\hat{\pi}(W_i))\hat{\Phi}_1(W_i) - \hat{\pi}(W_i)\hat{\Phi}_0(W_i) \Big),\] \[S_i := \rho_0(W_i) \Big(\Delta Y_i -(1-\pi(W_i))\Phi_1(W_i) - \pi(W_i)\Phi_0(W_i) \Big),\] $\Pi_{n, X|Z} :=P_{Z}\bm{X}$, $P_{Z}:=\Psi_{n}\left(\Psi_{n}^{\top}\Psi_{n}\right)^{-1}\Psi_{n}^{\top}$ and $\tilde{\bm{X}} := \bm{X} - \Pi_{n, X|Z} $, where $\Psi_{n}:=\Psi_{n}(\bm{Z})=\left(\psi^{k_{n}}(Z_{1})^{\top},\dots,\psi^{k_{n}}(Z_{n})^{\top}\right)^{\top}$ is a $n\times k_n$ matrix. Let $\bar{\bm{X}}=\bm{X}-\Pi_{X|Z}=(X_1^{\top}-\Pi^{\top}_{X_1|Z_1},...,X_n^{\top}-\Pi^{\top}_{X_n|Z_n})^{\top},$ where $\Pi_{X|Z}=\mathbb{E}[X_i|Z_i]$.

Define $\bm{\eta}_{n}=\bm{S}-\bm{X}\beta_{0}-\Psi_{n}\gamma_{n0}=\bm{\epsilon}+\bm{r_{n}}$, where $\bm{\epsilon}=(\epsilon_{1},\dots,\epsilon_{i},\dots,\epsilon_{n})^{\top}$, $\epsilon_{i}=\frac{D_{i}}{\pi_{i}}\epsilon_{1i}+\frac{1-D_{i}}{1-\pi_{i}}\epsilon_{0i}$, $\epsilon_{0i}=\Delta Y_{0i}-\Phi_{0}(W_i)$, $\epsilon_{1i}=\Delta Y_{1i}-\Phi_{1}(W_i)$, $\bm{r_{n}} = (r_{n1},\dots, r_{ni}, \dots, r_{nn})^{\top}$. We have the following decomposition: \[\|\bm{X}\beta_0+f_{n0}\|_{n}^2 = \|\tilde{\bm{X}}\beta_0\|_{n}^2 + \|\Pi_{n,X|Z}\beta_0 + f_{n0}\|_{n}^2.\]

In Assumption (ref), (i) is Assumption 3.3 in Abadie2005. We follow the same sampling scheme to consider repeated cross sections; (ii) can be relaxed if we add a continuous nonnegative weight function in the definition of $\|\cdot\|_{\infty,\alpha}$ in ((ref)) freyberger2019practical; and (iii) and (iv) are standard identification conditions for a partially linear model. Assumption (ref) allows the error terms to be non-identically distributed and be conditionally heteroskedastic. One can replace it by the stronger sub-Gaussian assumption often used in the literature.

As is standard in the literature for high-dimensional data, we introduce the restricted eigenvalue condition for as \[ \Lambda^{2}_{\bar{X}}(s_0):=\min_{\delta\in\mathbb{R}^{p}\backslash\{0\},\|\delta_{S_{0}^{C}}\|_{1}\leq3\sqrt{s_{0}}\|\delta_{S_{0}}\|_{2}}\frac{\delta^{\top}\mathbb{E}\left[\bar{X}_{i}\bar{X}_{i}^{\top}\right]\delta}{\|\delta_{S_{0}}\|_{2}^{2}}>0 \] Let $\Sigma_{\bar{X}}=\mathbb{E}[\bar{X}_{i}\bar{X}_{i}^{\top}]$, $\Sigma_{\tilde{X}}=\mathbb{E}[\tilde{X}_{i}\tilde{X}_{i}^{\top}]$ and $\Sigma_{\Pi}=\mathbb{E}[\Pi_{X_{i}|Z_{i}}\Pi_{X_{i}|Z_{i}}^{\top}]$.

Assumption (ref) is standard in the literature: (i) can be relaxed if we replace it by some uniform moment conditions discussed in caner2018asymptotically; (ii) and (iii) restrict the eigenvalues. In particular, (ii) is a restricted eigenvalue condition.

Assumption (ref) is Assumption A.3 in BCCK2015. Note that $\mathcal{F}$ is a set of functions $f$ in $C^{\alpha}_M(\mathcal{Z})$, thus, $\|f\|_{\infty,\alpha}$ is bounded from above uniformly over all $f\in \mathcal{F}$. Then for instance, $c_{k_n}=O(k^{-\alpha/d_{z}})$ for the polynomial series and $c_{k_n}=O(k^{-(\alpha\wedge \alpha_0)/d_{z}})$ for splines with order $\alpha_0$. Assumption (ref) (i) and the first two conditions in (ii) are also standard in the literature. Given this assumption, it is without loss of generality to normalize $\mathbb{E}[\psi^{k_n}(Z_{i})\psi^{k_n}(Z_{i})^{\top}]=I_{k_n}$. The condition $k_{n}\xi_{0}(k_{n})^{2}\log p/n=O_{p}(1)$ in Assumption (ref)(ii) is new. It is a mild condition on the relationship between $p$ and $k_n$. Assumption (ref)(iii) is a smoothness condition on the approximation error of the projection $\Pi_{n,X|Z}$ to $\Pi_{X|Z}$.

Assumption (ref) imposes moderate conditions on the first stage approximations of the nuisance functions ${\pi}(W_i)$, ${\Phi}_1(W_i)$ and ${\Phi}_0(W_i)$. Only the interaction terms between ${\pi}(W_i)$ and ${\Phi}_1(W_i)$ or ${\pi}(W_i)$ and ${\Phi}_0(W_i)$ are required to converge at a mild rate $O_p(\log p/n \vee k_n\log k_n/n)$. This demonstrates the double robustness properties of our estimator such that when either ${\pi}(W_i)$ or $({\Phi}_1(W_i), {\pi}(W_i))$ are correctly specified, the desired rate of convergence in Theorem (ref) can be achieved. As pointed out in Chernozhukov2016, the benefit of using sample-splitting is that it makes the entropy condition become very weak, allowing machine learning methods (e.g. random forest, boosted trees, deep neural nets, and their aggregated and hybrid versions) to be applied to estimate the functions $\hat{\pi}(W_i)$, $\hat{\Phi}_1(W_i)$ and $\hat{\Phi}_0(W_i)$. One can provide more primitive conditions to verify these rates for each given machine learning method of chosen.

Theorem (ref) establishes the rate of convergence for our estimator. We show that for the parametric estimator, similar to the one for the high-dimensional linear regressors (e.g., Theorem 6.1 in buhlmann2011statistics), its convergence rate depends on the rate of the tuning parameter $\lambda$ and the level of sparsity $s_0$. When $\lambda=O(\sqrt{\log p/n})$, we have $\|\hat{\beta}-\beta_0\|_{1}=O_p(s_0\sqrt{\log p/n})$, which is the same rate in lasso regression for high-dimensional linear models without unknown functions. For the nonparametric estimator, the convergence rate maintains the same rate as the one obtained in nonparametric regressor models (e.g., Theorem 4.1 in BCCK2015), which depends on the order of basis function $k_{n}$ and the approximation error. Unlike the results in the literature of semiparametric partially linear model when the dimension of $X$ is much smaller than sample size $n$, the convergence rate of the parametric estimator is slower than $\sqrt{n}$ due to high dimensionality of the model. It makes the inference problem challenging. As we will show in Section (ref), the asymptotic variance of the nonparametric estimator will contain a projection term that reflects the effect of the high-dimensional parametric estimation.

Asymptotic Inference

In many applications, practitioners are not only interested in the estimation of the treatment effect but also the uncertainty quantification of the estimated treatment effect. The latter provides the confidence of the treatment effect estimation and is a routine procedure in most causal inference problems. While the inferential properties under high-dimensional linear/generalized linear models have been extensively investigated in the recent literature zhang2011confidence,javanmard2013confidence,van2013asymptotically,belloni2013honest,ning2014general,ning2017likelihood,cai2015confidence,neykov2018unified,gold2020inference, the asymptotic inference under the DiD design has not been studied, especially in the partially linear model specification. In this section, we consider how to construct confidence intervals for the parametric component $\beta_0$ and the nonparametric component $f_0(z)$ for given $z\in \mathcal{Z}$.

Consider the inference problem for a linear combination of $\beta_0$, say $\xi^{\top}\beta_0$, for a known vector $\xi\in \mathbb{R}^p$. For instance, if we take $\xi$ as the unit basis vector $e_j=(0,...,0,1,0,...0)$ with the $j$th position being 1 and 0 otherwise, then the linear functional reduces to $\xi^{\top}\beta_0=(\beta_0)_j$, which is the $j$th component of the regression coefficient. Similarly, if we are interested in the prediction for a given test sample $X=x_0$ and $Z=z_0$, then the parameter of interest becomes $x_0^{\top}\beta+f(z_0)$. Thus, the inference problems can be decomposed into two problems: the inference on $x_0^{\top}\beta$ and the inference on $f(z_0)$. The former is again a linear combination of $\beta_0$ with $\xi=x_0$. The inference on $f(z_0)$ will be studied later in this section. To construct the confidence intervals for $\xi^{\top}\beta_0$, we extend the de-biasing approach to the DiD design under the partially linear model specification. Given the Lasso estimator $\hat\beta$, we propose the following de-biased Lasso estimator:

where

with $\hat{\Sigma}_{\tilde{X}}=\frac{1}{n}\sum_{i=1}^n \tilde X_i \tilde X_i^{\top}$ and $\lambda'$ as a tuning parameter. We will show that $\hat w$ is a consistent estimator of $w_{0}=\Sigma^{-1}_{\bar{X}}\xi$. Let $s_{w} = |\{k: w_{0k} \neq 0\}|$ as the size of non-zero elements in $w_{0}$.

Assumption (ref) is a stronger version of Assumption (ref), which is required for constructing the asymptotic normality.

Let $\sigma^{2}_i=\mathbb{E}[\epsilon^{2}_{i}|X_i]$, $V_\beta=\Sigma^{-1}_{\bar{X}}\Omega_\beta\Sigma^{-1}_{\bar{X}}$ with $\Sigma_{\bar{X}}=\mathbb{E}\left[\bar{X}_{i}\bar{X}_{i}^{\top}\right]$ and $\Omega_{\beta}:=\mathbb{E}\left[\sigma_{i}^{2}\bar{X}_{i}\bar{X}_{i}^{\top}\right]$. Let $\hat{V}_\beta=\hat{w}^{\top}\hat{\Omega}_{\beta}\hat{w}$ with $\hat{\Omega}_{\beta}:=\mathbb{E}_{n}\left[\hat{\sigma}_{i}^{2}\tilde{X}_{i}\tilde{X}_{i}^{\top}\right]$.

Theorem (ref) implies that we can construct an asymptotic $(1-\alpha)$ confidence interval for $\xi^{\top}\beta_0$ as $(\hat T-z_{1-\alpha/2}(\hat w^{\top} \hat{V}_\beta\hat w)^{1/2}, \hat T+z_{1-\alpha/2}(\hat w^{\top} \hat{V}_\beta\hat w)^{1/2}$, where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard normal random variable. Note that by constructing the Neyman orthogonality condition and by using an de-biased estimator, the asymptotic variance of the parametric parameter coincides the one in low-dimensional partially linear models Robinson1988root. In particular, with homoskedasticity $\mathbb{E}[\epsilon_{i}^{2}|X_{i}]=\sigma^{2}$, the asymptotic variance achieves the semiparametry efficiency bound $V_{\beta}=\sigma^{2}\mathbb{E}[\tilde{X}_{i}\tilde{X}_{i}^{\top}]^{-1}$. \\

In the following, we extend the de-biasing approach to construct the confidence intervals for $f(z)$ for any given $z\in \mathcal{Z}$, where we assume $d_z$ is much smaller than $n$ to avoid the curse of dimentionality problem for nonparametric estimation. Recall that $f(z)$ can be approximated in the sieve space by $\psi^{k_n}(z)^{\top}\gamma_{n0}$. To construct the confidence interval for $f(z)$, it suffices to apply the debias approach to the parameter $\gamma_{n}$, in which the parameter $\beta$ is treated as a high-dimensional nuisance parameter. To this end, we first derive the score function for $\gamma_n$ as $$ U_{\gamma_{n}}(\beta,M)=\mathbb{E}_n\left\{\left(\hat{\rho}_i \Big(\Delta Y_i -(1-\hat{\pi}_i)\hat{\Phi}_1(W_i) - \hat{\pi}_i\hat{\Phi}_0(W_i) \Big)-X_i^\top\beta - \psi^{k_n}(Z_i)^{\top}\gamma_n \right)(\psi^{k_n}(Z_i)- MX_i)\right\}, $$ where $M=\mathbb{E}(\psi^{k_n}(Z_i)X_i^{\top})\{\mathbb{E}X_i^{\otimes 2}\}^{-1}$ is a $k_n\times p$ matrix. One key property of the score function is that $U_{\gamma_{n}}(\beta,M)$ is insensitive to the unknown high-dimensional nuisance parameters $\beta$ and $M$. In fact, we will show below that $U_{\gamma_{n}}(\hat\beta,\hat M)=U_{\gamma_{n}}(\beta,M)+o_p(n^{-1/2})$ for some suitable estimators $\hat\beta$ and $\hat M$ to be defined later.

Given this score function $U_{\gamma_{n}}(\beta,M)$, we can define the one-step updated de-biased estimator as $\bar{f}(z):=\psi^{k_{n}}(z)^{\top}\bar{\gamma}_{n}$, where $$ \bar{\gamma}_{n}:=\hat{\gamma}_{n}-\hat{\Sigma}^{-1}_{f}\mathbb{E}_n\left\{\left(\hat{\rho}_i \Big(\Delta Y_i -(1-\hat{\pi}_i)\hat{\Phi}_1(W_i) - \hat{\pi}_i\hat{\Phi}_0(W_i) \Big)-X_i^\top\hat\beta - \hat f(Z_i) \right)(\psi^{k_n}(Z_i)-\hat MX_i)\right\} $$ where $\hat{\Sigma}_f = \mathbb{E}_n\left\{(\psi^{k_n}(Z_i)-\hat MX_i)\psi^{k_n}(Z_i)^{\top}\right\}$ and $\hat M=[\hat M_1,...,\hat M_j,...,\hat M_{k_n}]^{\top}$ with

Let $s_{m}=\left|k:\left\{ \left(\mathbb{E}[\psi^{k_{n}}(Z_{i})X_{i}^{\top}]\right)\left(\mathbb{E}\left[X_{i}X_{i}^{\top}\right]^{-1}\right)\right\} _{k}\neq0\right| $ as the size of non-zeros elements in $\left(\mathbb{E}[\psi^{k_{n}}(Z_{i})X_{i}^{\top}]\right)\left(\mathbb{E}\left[X_{i}X_{i}^{\top}\right]^{-1}\right)$.

Let $\sigma_{z}^{2}=\psi^{k_{n}}(z)^{\top}V_{f}\psi^{k_{n}}(z)$ with $V_{f}=\Sigma_{f}^{-1}\Omega_{f}\Sigma_{f}^{-1}$, $\Sigma_{f}=\mathbb{E}\left[\left(\psi^{k_{n}}(Z_{i})-MX_{i}\right)\psi^{k_{n}}(Z_{i})^{\top}\right]$ and $\Omega_{f}=\mathbb{E}\left[\sigma_{i}^{2}\psi^{k_{n}}(Z_{i})\psi^{k_{n}}(Z_{i})^{\top}\right]-M\mathbb{E}\left[\sigma_{i}^{2}X_{i}X_{i}^{\top}\right]M^{\top}$.

We define the sample analogs similarly. Let $\hat{\sigma}_{z}^{2}=\psi^{k_{n}}(z)^{\top}\hat{V}_{f}\psi^{k_{n}}(z)$ with $\hat{V}_{f}=\hat{\Sigma}_{f}^{-1}\hat{\Omega}_{f}\hat{\Sigma}_{f}^{-1}$, $\hat{\Sigma}_{f}=\mathbb{E}\left[\left(\psi^{k_{n}}(Z_{i})-\hat{M}X_{i}\right)\psi^{k_{n}}(Z_{i})^{\top}\right]$ and $\hat{\Omega}_{f}= \mathbb{E}_n\left[\hat{\sigma}^{2}_i\psi^{k_n}(Z_i)\psi^{k_n}(Z_i)^{\top}\right] - \hat M\mathbb{E}_n\left[\hat{\sigma}^{2}_iX_iX_i^{\top}\right]\hat M^{\top}$.

Theorem (ref) provides asymptotic theory that can be used to construct the confidence intervals for $f(z)$ for any $z$. Unlike the standard results in the nonparametric literature, we need to construct a de-biased $\bar{f}(z)$ estimator that corrects bias caused by estimating the high-dimensional parametric component of the partially linear model. Since the Lasso estimator of the parametric linear part has a convergence rate that is slower than $\sqrt{n}$, the asymptotic variance of $\bar{f}(z)$ contains a projection term that reflects the effect of the parametric component on the nonparametric component. We can construct an asymptotic $(1-\alpha)$ confidence interval for $f_{0}(z)$ as $\left[ \bar{f}(z)-z_{1-\alpha/2}\hat{\sigma}_{z}, \bar{f}(z)+z_{1-\alpha/2}\hat{\sigma}_{z}\right]$, where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard normal random variable. Combining the results in Theorems (ref) and (ref), one can easily construct the confidence interval for the heterogeneous ATT denoted as $\tau_0(w)$.

Simulation

We compare the finite sample performance of the doubly robust estimator proposed in this paper with the semiparametric DiD estimator in Abadie2005 when the latter is applicable. We consider two data generating processes. In the first setting (DGP1), we allow $Y^0(i, 0)$ and $Y^1(i,0)$ to follow standard normal distribution, where $Y^1(i, 1)$ and $Y^0(i, 1)$ are defined as follows:

where $X_i$ and $Z_i$ are generated from independently standard normal distributions. The errors $\epsilon_{i1}$ and $\epsilon_{i0}$ are independently generated from standard normal distributions. In the second setting (DGP2), we define \[Y(i, 0) = \tilde\epsilon_{i}\cdot (1/\sqrt{2} \cdot Z_i + 1/\sqrt{2} \cdot X_{i1})\] where $\tilde\epsilon_{i}$ is generated from a standard normal distribution and $X_i \sim N(0, \Sigma)$ where $\Sigma_{jk} = \rho^{|j-k|}$. This allows both heteroskedasticity in the error term as well as correlation among regressors. In both DGP1 and DGP2, we set $\beta^1_i = 2/i$ and $\beta^0_i = 1/i$ for $i \leq 15$ and $f(Z_i) = \exp(Z_i)$. The treatment assignment probability is based on a logistic distribution with \[\mathbb{P}(T_i=1) = 1-(1+\exp(X_i^\top\theta_0))^{-1},\] where $\theta_{0i} = 1/i$ for $i \leq 10$. In both setting, we use 8th degree trigonometric polynomial basis for the non-parametric estimation.

Table (ref) and (ref) summarize the results for the two settings. We report the average bias, average standard errors, average mean squared errors, average coverages for a $90\%$ confidence intervals as well as the average lengths for this confidence intervals separately for both the linear coefficients and the nonparametric coefficients. To compare with the parametric part, we report the coverages for the linear combination of the nonparametric coefficients. Divided by the standard error, it also converges to standard normal with the same condition in (ref). The “Dr-DiD" columns represent the results for the doubly robust diff-in-diff estimator and “semi-DiD" columns represent the results for the Abadie2005 estimator. We present results with $n$ varying from 200, 500 and 1000 and the dimension for linear specification $p$ varying from 10, 50 , 500 and 1000. Notice that the Abadie2005 estimator is infeasible when $n\leq p$ so we omit to report the “semi-DiD" results when $p =500, 1000$ and denote them as “-" in the tables. Furthermore, the variance for the Semi-Did estimator becomes large when $n$ is relatively large compared to $p$ (e.g. $n=200$ and $p=50$). Although Semi-DiD estimator can still be computed when $n=1000$ and $p=500$, we choose not to report this result because of its large variance .

As shown in both tables, the Dr-DiD estimator has a smaller standard error, RMSE and confidence interval length in both linear and nonparametric specifications. When $p=50$, the Semi-DiD estimator becomes too conservative and produces larger standard errors. On the other hand, the Dr-DiD estimator is also more robust comparing with the Semi-DiD estimator when switching from homoskedastic errors to heteroskedastic errors. More importantly, our experiments show that in finite sample, the Dr-DiD estimator can deliver reasonable estimates under high-dimensional settings.

Empirical Application

We use the proposed method to study the effect of increasing the minimum wage on unemployment rates at the county level. We use the same dataset collected by CallawayLi2020, which contains the county level unemployment rates from 2005 to 2007 before the Fair Minimum Wage Act was enacted in all states on May 25, 2007. Eleven states increased their minimum wage by the first quarter of 2007, while the other states did not increase their minimum wage until the federal minimum wage increased in July of 2007.\footnote{New Hampshire and Pennsylvania are dropped for the same reason as in CallawayLi2020 }

We explore the variation in adopting the minimum wage policy among different states to evaluate its impact on county level unemployment rates. CallawayLi2020 consider identification and estimation of the quantile treatment effect on the treated under a distributional extension of the mean difference in differences assumption with fixed-dimensional covariates. Differing from the work in CallawayLi2020, this work focuses on studying the impact of covariates on heterogeneous ATT in this DiD design. Our proposed method allows us to weaken the parallel trend assumption to the conditional parallel trend assumption by conditioning on a large amount of potential confounders. For example, states with smaller populations may have higher variation in the unemployment rates, thus moving at different trends as compared to states with larger populations. Our method also allows us to derive marginal effect given a specific covariate of interest. As a result, customized policy recommendations can be designed based on those results.

Figure (ref) plots the simple difference for the 2005 to 2007 difference in unemployment rate by median income (panel a) and population (panel b). We separate the counties in the control and treated states by red and blue color. The solid lines on the graphs are the local means for the control and treated groups. There is a general decrease in the unemployment rate from 2005 to 2007 across all counties as the change in unemployment rate is centered below 0. The difference between the red and blue lines is the standard DiD estimator under the unconditional parallel trend assumption. The decrease in the unemployment rate for the treated counties is lower than the control counties at the low income region, while not much difference between treated and control is observed at high income region. On the other hand, the decrease in the unemployment rate for the treated counties is consistently lower than the control counties regardless of population size.\\ Figure (ref) compares the semi-parametric diff-in-diff estimator (Semi-DiD: blue) with the doubly robust estimator (Dr-DiD: red) on the nonparametric component $f$. Both methods use 4th degree trigonometric polynomial basis to approximate $f$ and have partially linear forms. The Semi-DiD estimator controls only the underlying covariate of interest ($Z_i$), while the Dr-DiD estimator controls not only the underlying covariate of interest ($Z_i$) but also 703 covariates ($X_i$) in a linear additive form. These covariates include 38 county level characteristics, as well as all interactions between them. It is also worth pointing out that when computing the marginal effect for median income, population is used as a confounder in the linear part and vice versa when computing the marginal effect for population. The dashed lines are 95% confidence intervals for the estimator.

We find that both estimators show that regions with high median income levels and larger population sizes may suffer from an increase in unemployment rate due to the minimum wage policy, while no significant effect of the policy is detected for regions with lower median income levels and smaller population sizes. These findings coincide with the canonical economic theory on unemployment rate. For example, a higher median income level implies a higher substitution cost for a worker currently at minimum wage and thus leads to an increase in unemployment rate when the minimum wage rises. On the other hand, a region with a larger population size means more labor supply and thus a raise in minimum wage can lead to a surplus. The difference in the general direction of the results predicted by Figure (ref) and Figure (ref) indicates the potential severeness of confounding problems in this design.

Next, while both the Semi-DiD estimator and Dr-DiD estimator show no significant impact of the policy at low median income and thin population regions, the Dr-DiD shows a larger impact in regions with a higher median income level and a denser population than the Semi-DiD estimator. The Dr-DiD estimated effects are also more significant in those regions. This is due to the controlling of additional covariates that further alleviates the concern of confoundedness as well as reducing the uncertainly in the model to yield more accurate estimation.

Discussion

In this paper, we propose a new doubly robust two-stage difference-in-differences estimator that allows for, but does not require, the number of potential confounding covariates to be greater than the number of observations. Our estimator is robust to model miss-specification and a general set of machine learning tools can be used in our estimation procedure to estimate the propensity score. The outcome equation is modeled as a flexible partially linear form. The rate of convergence is derived for the new estimator and a novel de-bias procedure is proposed for inference. This allows the user to construct confidence intervals for the heterogeneous treatment effects. A simulation study shows promising finite sample performance of our estimator under different data generation processes. Our method is applied to study the effect of the Fair Minimum Wage Act on local unemployment rates and show heterogenous effects could rise due to the differences in demographics. Moreover, an R package for implementing the proposed method is available on Github. More work remains to be done. For example, it will be interesting to consider a similar estimation and inference strategy for panel data, or develop estimators for quantile treatment effect on the treated with high-dimensional covariates. We leave these topics for future studies.