Empirical Likelihood Covariate Adjustment for Regression Discontinuity Designs

Empirical Likelihood Covariate Adjustment for Regression Discontinuity Designs This version: April 22, 2024

Introduction

The RD design resembles a randomized experiment conducted near the cut-off of the score (forcing variable) and exploits the discontinuous variation in the probability of treatment to nonparametrically identify the local average treatment effect (LATE) at the cut-off under mild continuity assumptions on the latent variables.\footnote{In a recent study, hyytinen2018does confirmed that RD produces estimates that are in line with the results from a comparable experiment if inference is implemented with the method of calonico2014robust.} The transparent close-form identification (hahn2001identification) of the RD LATE calls for nonparametric estimation and inference methods as they avoid functional form assumptions. See Cattaneo2019 for a recent review of RD. In practical implementations, information from pre-treatment covariates (i.e., variables that have already been determined before the assignment of the treatment) is incorporated to enhance efficiency and compensate for the low accuracy of nonparametric methods. A widely used procedure is augmented local polynomial (LP) regression, where the covariates enter linearly. calonico2019regression formalize this augmented regression approach and derive its (first-order) asymptotic properties. CCFT shows that the augmented LP regression estimator consistently estimates the RD LATE under the covariate balance condition, i.e., the expectations of covariates coincide at both sides of the cut-off. Apart from CCFT, covariate adjustment for RD has received much attention in recent literature. See frolich2019including for an alternative approach that requires smoothing over covariates but allows for the potential failure of covariate balance. arai2021regression and Kreiss2022 extend CCFT's approach to control for a high-dimensional covariate vector by regularization. Noack2021 extend CCFT's linear regression adjustment to nonparametric adjustment with machine learning methods. See Cattaneo2021 for a recent review of covariate adjustment for RD.

This paper studies a novel and versatile approach based on (generalized) entropy balancing (EB) to incorporate covariates for RD. The recent literature on the estimation of the average treatment effect (ATE) under the unconfoundedness assumption and also broader causal inference literature (e.g., Doudchenko2016) flourishes with methods based on balancing. See BenMichael2021 for a review of this strand of literature. To the best of our knowledge, the balancing approach has not been investigated in the RD literature. In this paper, we follow CCFT to consider a potential outcome and covariate framework. Here, the covariate balance condition, which is a restriction on the population feature of the observed covariates, is directly implied by the predeterminedness (zero RD LATE on covariates) assumption and standard smoothness assumptions. Our balancing approach adjusts for covariates by using weights that achieve exact local covariate balance and have the least Kullback--Leibler (KL) divergence from the uniform weights. The EB estimator can be constructed in two intuitive steps: the first step computes the EB weights from a minimum relative entropy problem subject to the covariate balance constraints, and the second step replaces the uniform weights in the standard local polynomial RD estimator (without covariates) with the EB weights. The EB estimator can also be formulated as an empirical likelihood (EL) estimator, for which covariate balance translates to a set of over-identifying LP moment conditions and is used as “side information.” Therefore, our approach explicitly incorporates the covariate balance condition, which is treated as a maintained assumption in CCFT, into the estimation and inference procedure. We show in Theorem (ref) that the EB (EL) estimator is first-order equivalent to the regression adjustment estimator of CCFT. Although CCFT doubted whether covariate adjustment can always lead to asymptotic efficiency gain in RD estimation, it has been pointed out by Kreiss2022 that this is true. This paper provides an explanation of the asymptotic efficiency gain brought by covariate adjustment from the perspective of the generalized method of moments (GMM): the efficiency gain can be attributed to the efficient inclusion of covariate balance as side information (Remark (ref)). We also offer another explanation of the efficiency gain from the perspective of local randomization (Remark (ref)). Under CCFT's stronger version of covariate balance (see Page 446 of CCFT), incorporating functions of baseline covariates can further improve efficiency. Theorem (ref) shows that the asymptotic variance of the EB estimator incorporating basis functions of baseline covariates attains the lower bound derived in Noack2021, if the number of basis functions (i.e., the dimension of the corresponding linear sieve spaces) grows with the sample size.

Since the EB estimator can be formulated as an EL estimator, we expect that the favorable second-order properties (newey_smith_2004_higher) may also be shared by the EB estimator. Theorem (ref) shows that the EB estimator has a small (bounded) “nonlinearity bias”. Such a property is analogous to newey_smith_2004_higher. Then, we study covariate-adjusted EL inference for RD. A common advantage of EL inference is that it does not require calculating standard errors and explicit studentization. Theorem (ref) shows that the EL confidence set is a finite interval with probability approaching one. Theorem (ref) shows a new uniform-in-bandwidth extension of the standard Wilks theorem (i.e., the EL ratio is asymptotically $\chi^{2}$). Our uniform-in-bandwidth version adjusts for specification search over multiple bandwidths, known as bandwidth snooping (armstrong2017simple), and takes into account the effects from data-dependent bandwidths in a robust manner (Remarks (ref) and (ref)). It also provides a useful tool for sensitivity analysis in the sense of AK (Remark (ref)). By deriving distributional expansions, we investigate the second-order properties of our EL inference method and show that it enjoys a couple of nice properties in this setting. Theorem (ref) characterizes the leading coverage error term (i.e., the discrepancy between the nominal and finite-sample coverage probabilities; see, e.g., calonico2018optimal for Wald-type inference). We consider two choices of the LP order: one less than the assumed smoothness ($p$-th order) and exhausting the smoothness ($\left(p+1\right)$-th order). In the first case, the coverage optimal (CO) bandwidth, which is defined as the minimizer of this leading coverage error, has a simple closed form (Remark (ref)), which, to the best of our knowledge, cannot be obtained for Wald-type inference (calonico2018optimal). In both cases, the simple coverage expansion for the EL confidence sets makes analytical correction possible. The correction aims to remove the leading term in the coverage error and does not require resampling. The correction factor has a very simple form and thus can be estimated with good accuracy in finite samples. Remark (ref) proposes Analytically corrected likelihood ratio statistics and confidence sets for conducting covariate-adjusted RD inference. Remark (ref) combines the analytical correction and AK-type correction (Remark (ref)) and provides a more accurate uniform confidence band that is useful for sensitivity analysis and robust inference.

Theorem (ref) considers possible deviations from covariate balance and shows that the coverage accuracy of our proposed EL confidence set is highly insensitive to mild deviations (Remark (ref)), which we refer to as local imbalance in this paper. Failure of the covariate balance assumption may happen in a realistic situation when the balance condition holds for pre-treatment covariates in theory, but our sample observations on these covariates are contaminated (possibly due to measurement errors that occur after treatment) so that they are drawn from a perturbed population (kitamura2013robustness) that slightly violates the balance condition. When covariate balance does not hold exactly, the coverage accuracy of the EL confidence set stays relatively unaffected, while other inference methods may exhibit severe undercoverage (Remark (ref)). To the best of our knowledge, such a robustness property is novel in the literature.

Our balancing approach is versatile in dealing with covariate-adjustment estimation/inference for parameters and/or models beyond the standard RD, such as the treatment effect derivative (TED) of Dong2015 and nonlinear estimators for RD with limited outcome variables (e.g., xu2017regression,Xu2018). An algorithmic extension of CCFT's regression adjustment may not be straightforward in these scenarios. Indeed, applying our EB approach is about reweighting a sample-analogue-type estimator (without covariate) in the RD-related context using the EB weights that are fully determined by the covariate balance condition. It does not matter if the initial estimator (without covariate) involves derivative or nonlinear transformation. For this reason, our balancing approach serves as a useful complement to the regression adjustment. We consider the following example in this paper. In addition to the standard RD LATE parameter in the standard RD model, one may be interested in estimating other parameters that have important causal interpretations, such as the TED as a measurement of the external validity of RD. Theorem (ref) shows the efficiency gain of the simple TED estimator using our EB weights in place of uniform weights, for which the only assumption needed for consistency is covariate balance. Another class of problems that our approach can tackle is nonlinear estimators with limited outcome variables (e.g., xu2017regression,Xu2018). Estimators of xu2017regression,Xu2018 using the EB weights achieve desired properties (consistency and potential efficiency gain) under covariate balance. Lastly, various extensions to the standard RD model and estimation of the relevant causal parameters have been considered in the recent literature. Our approach has the potential to provide easy-to-implement covariate adjustment with clear causal interpretation. Further investigation is needed in a case-by-case manner.

Related literature. Our EB estimator resembles the method of Hainmueller2012,Chan2016 in the literature on balancing methods for estimating ATE under unconfoundedness. See Wong2017,Kallus2020,Hirshberg2021 for more recent development of this strand of literature. Graham2012 show that their balancing-type estimator enjoys a similar small nonlinearity bias property. EL and generalized EL (newey_smith_2004_higher) are popular alternatives to GMM, and they do not require first-step estimation of the efficient weighting matrix. See, e.g., kitamura_2006_el_review for a comprehensive review of EL and generalized EL. See, e.g., chen2000empirical,otsu2013estimation,otsu2015empirical,ma2019minimum for EL inference in the context of non-parametric curves. It was shown that EL has favorable properties relative to GMM. See, e.g., chen2007second,kitamura_2001_el_optimality,matsushita_otsu_2013_el,newey_smith_2004_higher,otsu_2010_bahadur,ma2017second among many others. In relation to the literature, otsu2015empirical proposed EL inference for RD without covariates. Their method was based on first-order conditions from standard local linear regression. This paper focuses on covariate adjustment and uses different moment conditions. In another related paper, ma2019minimum studied EL inference for the parameter of interest in the density discontinuity design (jales_yu_aie). Our paper uses a similar approach to covariate adjustment as Wu2011,Zhang2018 who formulated covariate balance in randomized experiments as moment conditions and proposed EL-type methods. We formulate local imbalance and study its impact on coverage accuracy by using standard local asymptotic analysis (e.g., the Pitman approach to local power analysis). Local imbalance can also be viewed as a special case of local misspecification in the GMM framework (see, e.g., armstrong2021sensitivity and references therein). However, the approach we take differs from those employed by papers in this strand of literature. Our approach follows Bravo:2003fj and is based on the second-order asymptotic expansion of the coverage probability under drifting data-generating processes (i.e., local imbalance).

Organization. Section (ref) quickly reviews the RD design. Section (ref) introduces our EB method for RD with covariates. Section (ref) provides results on the asymptotic properties of the EB estimator, including asymptotic normality with a discussion on the efficiency gain (Section (ref)), calculation of the nonlinearity bias (Section (ref)) and extension to balancing over sieve spaces (Section (ref)). In Section (ref), we consider inference using the likelihood ratio and show several properties, including a uniform-in-bandwidth Wilks theorem (Section (ref)), derivation of a simple analytical correction (Section (ref)), and sensitivity of the coverage probability to the covariate balance condition (Section (ref)). Section (ref) proposes a covariate-adjusted estimator of the TED and provides an asymptotic normality result that shows the efficiency gain. Sections (ref) and (ref) present results from simulation and empirical exercises. Section (ref) concludes. Proofs are collected in the online appendix (available at \url{ruc-econ.github.io/supplement_Rev_V12.pdf}).

Notation. $\sum_{i}$ is understood as $\sum_{i=1}^{n}$. “$a\coloneqq b$” means that $a$ is defined by $b$ and “$a\eqqcolon b$” means that $b$ is defined by $a$. For any $k$-times differentiable univariate function $f$, let $f^{\left(k\right)}$ denote the $k-$th order derivative. Let $\mathbbm{1}\left(\cdot\right)$ denote the indicator function. For a $d$-dimensional vector $x$, let $x^{\left(j\right)}$ denote its $j$-th coordinate, $x^{\top}$ denote its transpose, $x^{\otimes k}$ denote a vector of the distinct entries of $k$-th Kronecker power for $k=2,3,4$ ($x^{\otimes2}\coloneqq\mathrm{vech}\left(xx^{\top}\right)$, where $\mathrm{vech}\left(xx^{\top}\right)$ denotes the half vectorization of $xx^{\top}$, $x^{\otimes3}$ is the vector obtained by stacking $\left\{ x^{\left(j\right)}\mathrm{vech}\left(x_{j}x_{j}^{\top}\right):j=1,...,d\right\} $, where $x_{j}\coloneqq\left(x^{\left(j\right)},...,x^{\left(d\right)}\right)^{\top}$, and $x^{\otimes4}$ can be defined similarly) and $\left\Vert x\right\Vert $ denote its Euclidean norm. Let $\mathrm{I}_{K}$ denote the $K$-dimensional identity matrix. Let $0_{J\times K}$ denote the $J\times K$ matrix in which all elements are zeros. Let $0_{J}$ denote the $J$-dimensional vector in which all elements are zeros. $\mathrm{A}^{\left(jk\right)}$ denotes the $jk$-th element of a matrix $\mathrm{A}$. For a square matrix $\mathrm{A}$, let $\mathrm{tr}\left(\mathrm{A}\right)$ denote its trace and $\mathrm{mineig}\left(\mathrm{A}\right)$ and $\mathrm{maxeig}\left(\mathrm{A}\right)$ denote the smallest and the largest eigenvalues, respectively. For a real-valued function $f:\mathcal{X}\rightarrow\mathbb{R}$, let $\left\Vert f\right\Vert _{\infty}\coloneqq\mathrm{sup}_{x\in\mathcal{X}}\left|f\left(x\right)\right|$ denote the sup-norm. We write $a_{n}\asymp b_{n}$, if $a_{n}=O\left(b_{n}\right)$ and $b_{n}=O\left(a_{n}\right)$. Let $\mathrm{e}_{k,s}$ denote the $s$-th unit vector in $\mathbb{R}^{k}$.

Regression discontinuity designs

Let $X\in\mathbb{R}$ be a continuous score supported on $\left[\underline{x},\overline{x}\right]$. Let $f_{X}$ denote its density function. We normalize the cutoff point to zero (so that $0\in\left[\underline{x},\overline{x}\right]$ without loss of generality) for notational brevity. In this paper, we assume that $f_{X}$ is continuous at the cutoff. Denote $\varphi\coloneqq f_{X}\left(0\right)$ for simplicity. For a random vector (or matrix) $V$, denote $g_{V}\left(x\right)\coloneqq\mathrm{E}\left[V\mid X=x\right]$, $m_{V}\left(x\right)\coloneqq g_{V}\left(x\right)f_{X}\left(x\right)$ and $g_{V\mid ZX}\left(z,x\right)\coloneqq\mathrm{E}\left[V\mid Z=z,X=x\right]$. Denote $\mu_{V,-}^{\left(k\right)}\coloneqq\mathrm{lim}_{x\uparrow0}g_{V}^{\left(k\right)}\left(x\right)$ and $\psi_{V,-}^{\left(k\right)}\coloneqq\mathrm{lim}_{x\uparrow0}m_{V}^{\left(k\right)}\left(x\right)$. $\left(\mu_{V,+}^{\left(k\right)},\psi_{V,+}^{\left(k\right)}\right)$ are defined similarly with $\mathrm{lim}_{x\uparrow0}$ replaced by $\mathrm{lim}_{x\downarrow0}$. For simplicity, also denote $\mu_{V,\mathrm{s}}\coloneqq\mu_{V,\mathrm{s}}^{\left(0\right)}$, $\psi_{V,\mathrm{s}}\coloneqq\psi_{V,\mathrm{s}}^{\left(0\right)}$ ($\mathrm{s}\in\left\{ -,+\right\} $), $\mu_{V,\pm}\coloneqq\mu_{V,+}+\mu_{V,-}$, $\mu_{V,\dagger}\coloneqq\mu_{V,+}-\mu_{V,-}$, $\psi_{V,\pm}\coloneqq\psi_{V,+}+\psi_{V,-}$, $\psi_{V,\dagger}\coloneqq\psi_{V,+}-\psi_{V,-}$. Let $\mu_{V}$ ($\psi_{V}$) denote the common value if $\mu_{V,+}=\mu_{V,-}$ ($\psi_{V,+}=\psi_{V,-}$). For random vectors $V$ and $U$, $\mathrm{Var}_{\mid0^{+}}\left[U\right]$ is understood as $\mathrm{lim}_{x\downarrow0}\mathrm{Var}\left[U\mid X=x\right]=\mu_{UU^{\top},+}-\mu_{U,+}\mu_{U^{\top},+}$ and $\mathrm{Cov}_{\mid0^{+}}\left[V,U\right]$ is understood as $\mathrm{lim}_{x\downarrow0}\mathrm{Cov}\left[V,U\mid X=x\right]=\mu_{VU^{\top},+}-\mu_{V,+}\mu_{U^{\top},+}$. Similarly, $\mathrm{Var}_{\mid0^{-}}\left[U\right]\coloneqq\mathrm{lim}_{x\uparrow0}\mathrm{Var}\left[U\mid X=x\right]$ and $\mathrm{Cov}_{\mid0^{-}}\left[V,U\right]\coloneqq\mathrm{lim}_{x\uparrow0}\mathrm{Cov}\left[V,U\mid X=x\right]$. Also for notational simplicity, let $\mathrm{Var}_{\mid0^{\pm}}\left[U\right]\coloneqq\mathrm{Var}_{\mid0^{+}}\left[U\right]+\mathrm{Var}_{\mid0^{-}}\left[U\right]$ and $\mathrm{Cov}_{\mid0^{\pm}}\left[V,U\right]\coloneqq\mathrm{Cov}_{\mid0^{+}}\left[V,U\right]+\mathrm{Cov}_{\mid0^{-}}\left[V,U\right]$. $\mathrm{Var}_{\mid0}$ and $\mathrm{Cov}_{\mid0}$ are understood as $\mathrm{Var}\left[\cdot\mid X=0\right]$ and $\mathrm{Cov}\left[\cdot\mid X=0\right]$.

Let $Y\in\mathbb{R}$ denote the outcome variable, $D\in\left\{ 0,1\right\} $ be the binary treatment and $Z$ be pre-treatment covariates. Variables in $Z$ can be continuous, discrete or mixed. We observe $\left(Y,D,Z\right)$ and the score $X$. In an RD model, an incentive is assigned if $X\geq0$. In the sharp RD case $D=I\coloneqq\mathbbm{1}\left(X\geq0\right)$ (i.e., perfect compliance). The more general fuzzy RD model assumes $D\neq I$ but $g_{D}$ has a jump discontinuity at $x=0$ ($\mu_{D,+}\neq\mu_{D,-}$) due to the incentive. This is known as limited compliance in the literature. The RD model can be embedded in the potential outcome and treatment framework. Let $\left(Y\left(1\right),Y\left(0\right)\right)$ be the potential outcomes with or without treatment. Let $\left(D\left(1\right),D\left(0\right)\right)$ denote the potential treatments with or without incentives. The observed outcome $Y$ and treatment $D$ are determined by $Y=D\cdot Y\left(1\right)+\left(1-D\right)Y\left(0\right)$ and $D=I\cdot D\left(1\right)+\left(1-I\right)D\left(0\right)$ respectively. The complier group is defined to be individuals with $D\left(1\right)>D\left(0\right)$ (i.e., $\left(D\left(1\right),D\left(0\right)\right)=\left(1,0\right)$). We use “$\mathsf{co}$” to denote this event. Following CCFT, we let $\left(Z\left(1\right),Z\left(0\right)\right)$ denote potential covariates and then $Z=D\cdot Z\left(1\right)+\left(1-D\right)Z\left(0\right)$.\footnote{The RD design can be represented by a structural model. See dong2018alternative. $\left(Y,D,Z\right)$ are assumed to be generated by the structural model $Y=g\left(D,X,Z,\epsilon\right)$, $D=h\left(I,X,\eta\right)$ and $Z=m\left(D,X,\xi\right)$, where $\left(g,h,m\right)$ are unknown functions and $\left(\epsilon,\eta,\xi\right)$ are (potentially correlated) unobserved disturbances of unrestricted dimensionality. Then the potential outcomes, covariates and treatments are given by $Y\left(d\right)=g\left(d,X,Z,\epsilon\right)$, $D\left(d\right)=h\left(d,X,\eta\right)$ and $Z\left(d\right)=m\left(d,X,\xi\right)$.} Let $B\left(d\right)\coloneqq\left(Y\left(d\right),Z\left(d\right)\right)$, for $d\in\left\{ 0,1\right\} $. Denote $g_{dd'}\left(x\right)\coloneqq\mathrm{Pr}\left[D\left(1\right)=d,D\left(0\right)=d'\mid X=x\right]$ and $g_{B\left(j\right)\mid dd'}\left(x\right)\coloneqq\mathrm{E}\left[B\left(j\right)\mid D\left(1\right)=d,D\left(0\right)=d',X=x\right]$ for $\left(j,d,d'\right)\in\left\{ 0,1\right\} ^{3}$. The RD LATE (the average treatment effect for individuals with zero score in the complier group) is given by $\mathrm{E}\left[Y\left(1\right)-Y\left(0\right)\mid X=0,\mathsf{co}\right]$. The following assumption is sufficient for the identification in RD and is also imposed in CCFT.

It can be shown that under (a,b,c,d), the RD LATE is identified by the standard RD estimand $\vartheta\coloneqq\mu_{Y,\dagger}/\mu_{D,\dagger}$ (i.e., $\mathrm{E}\left[Y\left(1\right)-Y\left(0\right)\mid X=0,\mathsf{co}\right]=\vartheta$, see hahn2001identification,dong2018alternative and Arai2021 for testable implications of these assumptions), where $\vartheta$ is a population feature of the observed variables.\footnote{In the sharp RD model ($\mu_{D,+}=1$ and $\mu_{D,-}=0$ in this case) or under a stronger conditional independence assumption (hahn2001identification), a causal parameter that corresponds to a broader subpopulation (conditional average treatment effect) is identified by the same ratio: $\mathrm{E}\left[Y\left(1\right)-Y\left(0\right)\mid X=0\right]=\vartheta$ .} As in frolich2019including, the continuity assumption (a) can be viewed as an exclusion restriction. Intuitively, continuity of $g_{Y\left(j\right)\mid dd'}$ essentially requires that $Y\left(j\right)$ cannot depend on $I$ or (observed or unobserved) variables related to $I$ (so that their distributions change discontinuously at the cutoff). Since $Y\left(j\right)$ often depends on $Z$, continuity of $g_{Y\left(j\right)\mid dd'}$ also implicitly requires that the conditional distributions of $Z$ given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$ change smoothly at $x=0$. Since the distribution of $Z$ coincides with that of $Z\left(d\right)$ ($Z\left(d'\right)$), given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$ with $x\geq0$ ($x<0$), continuity of the conditional distribution of $Z$ given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$, for $\left(d,d'\right)\in\left\{ \left(1,1\right),\left(0,0\right),\left(1,0\right)\right\} $, holds if the conditional distributions of the potential covariates change smoothly at $x=0$ and the distribution of $Z\left(1\right)$ given $\mathsf{co}$ and $X=0$ is the same as that of $Z\left(0\right)$ given $\mathsf{co}$ and $X=0$. Following CCFT, we consider using weaker versions of these assumptions in (e,f). We consider using the strong versions in Section (ref). (e) essentially requires that the covariates satisfy the same exclusion restriction (not affected by $I$). It is clear from $g_{Z}\left(x\right)=\sum_{d,d'}g_{Z\mid dd'}\left(x\right)g_{dd'}\left(x\right)$, where “$\sum_{d,d'}$” is understood as “$\sum_{\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}}$”, that covariate balance $\mu_{Z,+}=\mu_{Z,-}$ holds as a testable implication for the population of the observed variables.\footnote{Indeed, $\mu_{Z,+}=\mu_{Z,-}$ is the null hypothesis of a popular falsification or placebo test for the RD model. See, e.g., lee2008randomized,Canay2017. Evidence against $\mu_{Z,+}=\mu_{Z,-}$ in the data (so that a hypothesis test of $\mu_{Z,+}=\mu_{Z,-}$ is rejected) casts doubts on the validity of the key identifying assumption of the RD design (i.e., Assumption (ref)(a)). While most empirical works conduct the balance test separately for each covariate, some researchers have noted that the problem of multiple testing may generate statistical imbalance of some covariates by chance. See, e.g., hyytinen2018does.}

Empirical entropy balancing

This section introduces the EB method. We quickly review the idea of entropy balancing and reweighting in the literature on ATE estimation under unconfoundedness (i.e., conditional independence of the potential outcomes and the treatment given the covariates). Then, we utilize the idea of EB to propose a new balancing-based method for covariate adjustment for RD.

In observational studies, because of the selection bias, the difference in the sample means corresponding to the treatment and control groups does not consistently estimate the ATE. The balancing weights satisfy the requirement that the weighted control (treatment) group sample moments of the covariates match the unweighted sample moments of the covariates of all units. Within all balancing weights, Hainmueller2012 defines the EB weights as those being as close as possible to the uniform weights in the sense of minimal relative entropy (KL divergence). Hainmueller2012 replaces the uniform weights used by the simple sample means with the EB weights. Chan2016 construct EB weights that equalize weighted and unweighted sample means of transformations of the covariates via basis functions. Chan2016 show that the estimator using these EB weights overcomes the selection bias under the unconfoundedness assumption if the number of basis functions of the covariates increases with the sample size.

From a GMM/EL perspective, in Hainmueller2012,Chan2016, the entropy balancing and reweighting approach uses weights under which some intentionally misspecified (biased) moment restrictions are satisfied to correct for the selection bias.\footnote{In observational data, the population moments of covariates in the control or treatment group may not be the same as the unconditional population moments, since the treatment status is not independent from the covariates.} In our RD case, the moment restrictions (balancing constraints) are correctly specified and entropy balancing and reweighting aim at enhancing efficiency (Section (ref)). In our case, the EB estimator can be formulated as a standard EL estimator (Section (ref)).

Entropy balancing for covariate adjustment in RD

Now we elaborate on the entropy balancing and reweighting approach to covariate adjustment in the RD context. Firstly, we introduce some notations. Let $K$ denote the kernel function and let $h$ denote the bandwidth. We assume that $h=h_{n}$ decreases with the sample size $n$. For notational simplicity, we suppress the dependence of $h$ on $n$. Let the data $\left\{ \left(Y_{i},D_{i},X_{i},Z_{i}\right)\right\} _{i=1}^{n}$ be i.i.d. copies of $\left(Y,D,X,Z\right)$. We drop the subscript $i$ when we refer to population-level estimands. Let $p\geq1$ be the integer-valued LP order. Let $r_{p}\left(t\right)\coloneqq\left(1,t,\ldots,t^{p}\right)^{\top}$ and let $\mathrm{H}$ be the $\left(p+1\right)\times\left(p+1\right)$ diagonal matrix with $\left(1,h,...,h^{p}\right)$ being on the diagonal. Denote

Let $\widehat{\Pi}_{p,+}$ be defined similarly by the right-hand side of ((ref)) with $\mathbbm{1}\left(X_{i}<0\right)$ replaced by $\mathbbm{1}\left(X_{i}>0\right)$. Let

Let $\widehat{W}_{p;+,i}$ be defined similarly by the right-hand side of ((ref)) with $\mathbbm{1}\left(X_{i}<0\right)$ and $\widehat{\Pi}_{p,-}$ replaced by $\mathbbm{1}\left(X_{i}>0\right)$ and $\widehat{\Pi}_{p,+}$.\footnote{We restrict the bandwidths on the left and the right of the cut-off to be the same. It is possible to extend all of the theorems in this paper to accommodate different bandwidths on different sides.}

Let $\widehat{W}_{p,i}\coloneqq\widehat{W}_{p;+,i}-\widehat{W}_{p;-,i}$. The standard LP regression estimator of $\vartheta$ is

where the numerator $\left(nh\right)^{-1}\text{\ensuremath{\sum_{i}\widehat{W}_{p,i}Y_{i}}}$ is the LP regression estimator of $\mu_{Y,\dagger}$ and the denominator is the LP regression estimator of $\mu_{D,\dagger}$.

Now we incorporate the covariate information to the standard LP estimator $\widehat{\vartheta}_{p}^{\mathsf{lp}}$ by reweighting its numerator and denominator using the EB weights computed from the covariate balance constraints . Denote $\bar{Z}_{i}\coloneqq\left(1,Z_{i}^{\top}\right)^{\top}$. We define EB weights $\left(w_{1}^{\mathsf{eb}},...,w_{n}^{\mathsf{eb}}\right)$ as the solution to the following minimum relative entropy problem:

where $\mathit{KL}\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)\coloneqq-\sum_{i}\mathrm{log}\left(n\cdot w_{i}\right)/n$ is the KL divergence from $\left(w_{1},...,w_{n}\right)$ to the uniform weights $\left(1/n,...,1/n\right)$. The construction of these balancing weights is similar to those in Hainmueller2012. By solving the minimization problem ((ref)), we find the set of weights with the least KL divergence from the uniform weights among these balancing weights. The uniform weights satisfy the important finite-sample property of $\sum_{i}\widehat{W}_{p,i}=0$. By requiring $\sum_{i}w_{i}

=0$ in the constraint of (\ref{eq:entropy balancing}), we require that the balancing weights satisfy the same property. The balancing weights should also satisfy $\sum_{i}w_{i}

=\sum_{i}w_{i}

$. This requires that the (kernel-weighted) local averages of the covariates on both sides of the thresholds coincide in finite samples under the new weights for the data points. To solve for the optimal weights, we use strong duality and concentrate out $\left(w_{1},...,w_{n}\right)$ to obtain the following dual characterization

where

is the Lagrangian multiplier. Computing the EB weights requires dealing with a well-understood convex optimization problem ((ref)) that can be solved by the Newton algorithm. The domain of its objective function is the convex set $\left\{ \lambda:1+\lambda^{\top}\left(\widehat{W}_{p,i}\bar{Z}_{i}\right)>0\textrm{ for all }i\right\} $. The algorithm should either take these constraints into account or use a modified objective function defined for all $\lambda\in\mathbb{R}^{d_{z}+1}$. ((ref)) has no solution if the origin is not an interior point of the convex hull of $\left\{ \widehat{W}_{p,1}\bar{Z}_{1},...,\widehat{W}_{p,n}\bar{Z}_{n}\right\} $.\footnote{By arguments similar to those in owen_2001_el_book, we can show that if covariate balance holds, the origin lies in the convex hull with probability approaching one.} See kitamura_2006_el_review and owen_2001_el_book for more algorithmic details.\footnote{In the “no solution” scenario, the Newton algorithm trying to solve ((ref)) returns a sequence of vectors with diverging lengths. In this scenario, after the algorithm terminates (either the gradient is sufficiently small or the maximal number of iterations is reached), we would get weights not summing up to one in the former case (owen_2001_el_book) or a large gradient in the latter case. In our simulation studies and computation for the empirical application, we use the Matlab code written by Kirill Evdokimov and Yuichi Kitamura (\url{https://kitamura.sites.yale.edu/matlabstata-codes-el}) and hardly see any “no solution” case.}

We propose the following empirical entropy balancing estimator by reweighting the numerator and denominator of $\widehat{\vartheta}_{p}^{\mathsf{lp}}$ in ((ref)) using the EB weights $w_{i}^{\mathsf{eb}}$ defined by ((ref)):

The reweighting form of our EB estimator $\widehat{\vartheta}_{p}^{\mathsf{eb}}$ has a clear causal interpretation. The continuity and predeterminedness assumptions in Assumption (ref) imply a restriction $\mu_{Z,+}=\mu_{Z,-}$ on the population distribution of the observed covariates. $\widehat{\vartheta}_{p}^{\mathsf{eb}}$ directly uses weights that explicitly exploit such information from the covariates.

Generalized entropy balancing

The entropy balancing approach looks for balancing weights closest to the uniform weights, where “closeness” is measured by the KL divergence. It is useful to consider the following extension. Let $\left(p_{1},...,p_{n}\right)$ and $\left(p_{1}',...,p_{n}'\right)$ be two sets of probability masses. For any $\varrho\in\mathbb{R}$, let

be the Cressie-Read divergence from $\left(p_{1},...,p_{n}\right)$ to $\left(p_{1}',...,p_{n}'\right)$. Taking $\varrho=0$ gives the KL divergence from $\left(p_{1},...,p_{n}\right)$ to $\left(p_{1}',...,p_{n}'\right)$. Taking $\varrho=-1$ gives the KL divergence from $\left(p_{1}',...,p_{n}'\right)$ to $\left(p_{1},...,p_{n}\right)$.\footnote{$D_{0}\left(p_{1},...,p_{n}\parallel p_{1}',...,p_{n}'\right)$ (or $D_{-1}\left(p_{1},...,p_{n}\parallel p_{1}',...,p_{n}'\right)$) is defined as the limit of the right hand side of ((ref)) as $\varrho\rightarrow0$ (or $\varrho\rightarrow-1$).} The generalized balancing estimator is based on the weights $\left(w_{\varrho,1}^{\mathsf{gb}},...,w_{\varrho,n}^{\mathsf{gb}}\right)$ that solve

We define the generalized balancing estimator $\widehat{\vartheta}_{\varrho,p}^{\mathsf{gb}}$ by the right hand side of ((ref)) with $w_{i}^{\mathsf{eb}}$ replaced by $w_{\varrho,i}^{\mathsf{gb}}$.

Interestingly, although CCFT takes an augmented regression approach to incorporate pre-treatment covariates, a slight modification of CCFT's estimator can be written as a generalized balancing estimator with $\varrho=-2$. It is easy to see that $D_{-2}\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)$ is proportional to the square of the Euclidean distance between $\left(w_{1},...,w_{n}\right)$ and the uniform weights. CCFT's covariate adjusted estimator $\widehat{\vartheta}_{Y,p}^{\mathsf{CCFT}}$ for $\mu_{Y,\dagger}$ is given by the regression coefficient of $I_{i}\coloneqq\mathbbm{1}\left(X_{i}\geq0\right)$ in

Similarly, CCFT's estimator $\widehat{\vartheta}_{D,p}^{\mathsf{CCFT}}$ for $\mu_{D,\dagger}$ is defined by the right hand side of the above equation with $Y_{i}$ replaced by $D_{i}$. Then CCFT's estimator of $\vartheta$ is $\widehat{\vartheta}_{p}^{\mathsf{CCFT}}\coloneqq\widehat{\vartheta}_{Y,p}^{\mathsf{CCFT}}/\widehat{\vartheta}_{D,p}^{\mathsf{CCFT}}$. CCFT shows that by the partitioned regression argument, the numerator $\widehat{\vartheta}_{Y,p}^{\mathsf{CCFT}}$ can be written as

where $\widehat{\gamma}_{Y}^{\mathsf{CCFT}}$ is a consistent estimator of $\gamma_{Y}\coloneqq\left(\mathrm{Var}_{\mid0^{\pm}}\left[Z\right]\right)^{-1}\mathrm{Cov}_{\mid0^{\pm}}\left[Z,Y\right]$. A similar result holds for $\widehat{\vartheta}_{D,p}^{\mathsf{CCFT}}$. To find the optimal weights that solve ((ref)) with $\varrho=-2$, we again apply the Lagrangian multiplier method. What differs from EB is that in this case the Lagrangian multiplier has an explicit form. Then we can see that $\sum_{i}w_{-2,i}^{\mathsf{gb}}\widehat{W}_{p,i}Y_{i}/h$ (or $\sum_{i}w_{-2,i}^{\mathsf{gb}}\widehat{W}_{p,i}D_{i}/h$) can be written in the form of ((ref)) with a slightly different estimator of $\gamma_{Y}$ (or $\gamma_{D}$). More details can be found in Section S9 of our online supplement. Therefore, our formulation provides information-theoretic and balancing interpretation of CCFT's estimator. The generalized balancing estimators are all first-order equivalent to the EB estimator in the sense that the conclusion of Theorem (ref) also holds for them.

Connection to empirical likelihood

We now show that the EB estimator can be formulated as an EL estimator which incorporates covariate balance as side information. The RD estimand $\vartheta$, which has causal interpretation under the identifying assumptions of the RD model, can be approximately identified by a moment condition. Note that \[ \underset{x\downarrow0}{\mathrm{lim}}\,\mathrm{E}\left[Y-\theta D\mid X=x\right]=\underset{x\uparrow0}{\mathrm{lim}}\,\mathrm{E}\left[Y-\theta D\mid X=x\right]\textrm{ if and only if }\theta=\vartheta. \] By the standard LP regression theory, we have \[ \frac{1}{nh}\sum_{i}\widehat{W}_{p,i}\left(Y_{i}-\theta D_{i}\right)\rightarrow_{p}\underset{x\downarrow0}{\mathrm{lim}}\,\mathrm{E}\left[Y-\theta D\mid X=x\right]-\underset{x\uparrow0}{\mathrm{lim}}\,\mathrm{E}\left[Y-\theta D\mid X=x\right], \] under standard assumptions. Under covariate balance, $\left(nh\right)^{-1}\sum_{i}\widehat{W}_{p,i}\bar{Z}_{i}\rightarrow_{p}0_{d_{z}+1}$. We consider the following EL-type criterion function:

Note that we have $2+d_{z}$ LP moment conditions that approximately identify one parameter of interest $\vartheta$. Note that the covariate balance condition provides a set of over-identifying moment restrictions. We can easily see that the EB estimator is also an EL estimator, defined as a minimizer of the EL criterion function. It is clear that $\ell_{p}^{\mathsf{\mathsf{el}}}\left(\theta\mid h\right)\geq-n^{-1}\sum_{i}\mathrm{log}\left(n\cdot w_{i}^{\mathsf{eb}}\right)$ for all $\theta$, since $-n^{-1}\sum_{i}\mathrm{log}\left(n\cdot w_{i}^{\mathsf{eb}}\right)$ is the minimum corresponding to a larger constraint set. Since the constraint set of ((ref)) with $\theta=\widehat{\vartheta}_{p}^{\mathsf{\mathsf{eb}}}$ contains the EB weights, we have $\ell_{p}^{\mathsf{\mathsf{el}}}\left(\widehat{\vartheta}_{p}^{\mathsf{\mathsf{eb}}}\mid h\right)\leq-n^{-1}\sum_{i}\mathrm{log}\left(n\cdot w_{i}^{\mathsf{eb}}\right)$. Therefore, $\widehat{\vartheta}_{p}^{\mathsf{\mathsf{eb}}}$ is a minimizer of $\ell_{p}^{\mathsf{\mathsf{el}}}\left(\cdot\mid h\right)$.

We consider replacing the kernel-dependent weight $\widehat{W}_{p,i}$ in ((ref)) by the weight from population-level LP fitting (or solving a minimum contrast (MC) problem) and define the MC-EL estimator as the minimizer. See bickel2015mathematical for more details about the construction of the population-level LP fitting. Denote $\mathrm{V}_{p;-}\coloneqq\int_{-1}^{0}r_{p}\left(t\right)r_{p}^{\top}\left(t\right)K\left(t\right)\mathrm{d}t$ and $\mathcal{K}_{p;-}\left(t\right)\coloneqq\mathrm{e}_{p+1,1}^{\top}\mathrm{V}_{p;-}^{-1}r_{p}\left(t\right)K\left(t\right)$. Let $\left(\mathrm{V}_{p;+},\mathcal{K}_{p;+}\right)$ be defined by the same equations with the integral range $\left[-1,0\right]$ replaced by $\left[0,1\right]$. $\left(\mathcal{K}_{p;+},\mathcal{K}_{p;-}\right)$ coincide with the “equivalent kernel” associated with the LP regression. See, e.g., Section S2.1 of AK. Let $W_{p;-,i}\coloneqq\mathbbm{1}\left(X_{i}<0\right)\mathcal{K}_{p;-}\left(X_{i}/h\right)$, $W_{p;+,i}\coloneqq\mathbbm{1}\left(X_{i}>0\right)\mathcal{K}_{p;+}\left(X_{i}/h\right)$ and $W_{p,i}\coloneqq W_{p;+,i}-W_{p;-,i}$. By Taylor expansion (see Jiang_2003), $\mathrm{E}\left[W_{p}\left(Y-\theta D\right)\right]=O\left(h^{p+2}\right)$ if and only if $\theta=\vartheta$ and $\mathrm{E}\left[W_{p}\cdot\bar{Z}\right]=O\left(h^{p+2}\right)$ under suitable smoothness assumptions. Let the MC-EL criterion function $\ell_{p}^{\mathsf{mc}}$ be defined by the right hand side of ((ref)) with $\widehat{W}_{p,i}$ replaced by $W_{p,i}$. The $p$-th order MC-EL estimator is given by $\widehat{\vartheta}_{p}^{\mathsf{mc}}\coloneqq\mathrm{argmin}_{\theta}\ell_{p}^{\mathsf{mc}}\left(\theta\mid h\right)$. Similar derivations show that the estimator can be written as a balancing-type estimator. $\widehat{\vartheta}_{p}^{\mathsf{mc}}$ is equal to the right hand side of ((ref)) with $\left(\widehat{W}_{p,i},w_{i}^{\mathsf{eb}}\right)$ replaced by $\left(W_{p,i},w_{i}^{\mathsf{mc}}\right)$, where $\left(w_{i}^{\mathsf{mc}},\lambda_{p}^{\mathsf{mc}}\right)$ are defined by the right hand sides of ((ref)) and ((ref)) with $\widehat{W}_{p,i}$ replaced by $W_{p,i}$. The MC-EL estimator has similar asymptotic properties as the EB estimator (see Theorems (ref) and (ref) ahead). The MC-EL criterion function is useful for constructing confidence sets for $\vartheta$ with favorable second-order properties (see Section (ref) ahead).

Properties of the empirical likelihood (balancing) estimator

In this section, we show several large-sample properties of the proposed empirical balancing (likelihood) estimator. Section (ref) gives the asymptotic normality result of the estimators proposed in the preceding section. We compare our result with that of CCFT's estimator and discuss the efficiency gain brought by the covariates. Section (ref) gives a result on the “nonlinearity bias” of the proposed estimators. We argue that our estimators have small nonlinearity biases, especially in the situation when a relatively large number of valid covariates satisfying the balance condition are available. In Section (ref), we consider the situation when strong continuity and predeterminedness assumptions hold. The main result in Section (ref) shows that an extension of our EB estimator whose weights balance functions of covariates in a sequence of linear sieve spaces achieves the variance lower bound derived in Noack2021.

Efficiency gain

This section shows asymptotic normality of the EB and MC-EL estimators, and gives the expression for the asymptotic mean square error (AMSE). We then compare our results with the asymptotic result from CCFT. Let $B\coloneqq\left(Y,Z^{\top}\right)^{\top}$, $\bar{B}\coloneqq\left(Y,D,Z^{\top}\right)^{\top}$ and $M\coloneqq Y-\vartheta D$. The following assumptions are imposed on the population distribution of the observed variables. Let $\mathbb{B}\subseteq\left[\underline{x},\overline{x}\right]$ denote a neighborhood around 0.

Assumption (ref) parallels Assumption SA-5 of CCFT. We will invoke it directly in the proofs. These assumptions are satisfied under suitable conditions imposed on the population distribution of the latent variables as in Assumption (ref). Since $B=D\left(1\right)B\left(1\right)+\left(1-D\left(1\right)\right)B\left(0\right)$ if $X\geq0$ and $B=D\left(0\right)B\left(1\right)+\left(1-D\left(0\right)\right)B\left(0\right)$ if $X<0$, by the law of iterated expectations (LIE), for any function $\varphi\left(\cdot,\cdot\right)$, we have

(a) is satisfied if $\left(g_{B\left(d\right)\mid dd'},g_{B\left(d'\right)\mid dd'},g_{dd'}\right)$ are $\left(p+1\right)$-times continuously differentiable on $\mathbb{B}$ with uniformly continuous derivatives for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$. The smoothness level in (a) is similar to that commonly assumed in the literature, i.e., the minimal smoothness level ($p+1$) such that that the leading smoothing bias term of the estimator (using $p$-th order LP) can be explicitly characterized. (b) is satisfied if we impose the additional condition that for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$, $\left(g_{B\left(d\right){}^{\otimes2}\mid dd'},g_{B\left(d'\right){}^{\otimes2}\mid dd'}\right)$ are uniformly continuous on $\mathbb{B}$. (a,c) also guarantee that $m_{B}=g_{B}f_{X}$ have uniformly continuous derivatives up to $\left(p+1\right)$-th order on the left and right neighborhoods of $0$. Existence of $\mathrm{Var}_{0^{+}}\left[\left(M,Z^{\top}\right)^{\top}\right]$ and $\mathrm{Var}_{\mid0^{-}}\left[\left(M,Z^{\top}\right)^{\top}\right]$ is guaranteed by (b). By the law of total variance (writing $\mathrm{Var}\left[\left(M,Z^{\top}\right)^{\top}\mid X\right]$ as the sum of $\mathrm{E}\left[\mathrm{Var}\left[\left(M,Z^{\top}\right)^{\top}\mid D\left(1\right),D\left(0\right),X\right]\mid X\right]$ and $\mathrm{Var}\left[\mathrm{E}\left[\left(M,Z^{\top}\right)^{\top}\mid D\left(1\right),D\left(0\right),X\right]\mid X\right]$) and (b), $\mathrm{Var}_{\mid0^{+}}\left[\left(M,Z^{\top}\right)^{\top}\right]$ (or $\mathrm{Var}_{\mid0^{-}}\left[\left(M,Z^{\top}\right)^{\top}\right]$) is guaranteed to be positive definite if $\mathrm{Var}\left[B\left(1\right)\mid X=0,\mathsf{co}\right]$ (or $\mathrm{Var}\left[B\left(0\right)\mid X=0,\mathsf{co}\right]$) is positive definite.

(a) is standard and also imposed in CCFT. (b) is also found in AK. (a) implies that $\mathcal{K}_{p;+}\left(t\right)=\mathcal{K}_{p;-}\left(-t\right)$ and therefore (b) also holds for $\mathcal{K}_{p;-}$. Denote $\omega_{p;+}^{j,k}\coloneqq\int_{0}^{1}t^{j}\mathcal{K}_{p;+}^{k}\left(t\right)\mathrm{d}t$ and $\omega_{p;-}^{j,k}\coloneqq\int_{-1}^{0}t^{j}\mathcal{K}_{p;-}^{k}\left(t\right)\mathrm{d}t$. It is easy to see that $\omega_{p;-}^{0,k}=\omega_{p;+}^{0,k}\eqqcolon\omega_{p}^{0,k}$. Let $\gamma_{M}\coloneqq\left(\mathrm{Var}_{\mid0^{\pm}}\left[Z\right]\right)^{-1}\mathrm{Cov}_{\mid0^{\pm}}\left[Z,M\right]$, $\epsilon\coloneqq M-Z^{\top}\gamma_{M}$ and $\sigma^{2}\coloneqq\mathrm{Var}_{\mid0^{\pm}}\left[M\right]-\mathrm{Cov}_{\mid0^{\pm}}\left[M,Z\right]\cdot\gamma_{M}=\mathrm{Var}_{\mid0^{\pm}}\left[\epsilon\right]$. Existence of these quantities is guaranteed by Assumption (ref)(b). Under Assumption (ref)(d), $\sigma^{2}$ is strictly positive. Under Assumption (ref), $\mu_{\epsilon,+}=\mu_{\epsilon,-}\eqqcolon\mu_{\epsilon}$. Assumption (ref)(a) guarantees that $g_{\epsilon}$ and $m_{\epsilon}$ admit continuous derivatives up to $\left(p+1\right)$-th order on the left and right neighborhoods of $0$ so that the leading bias terms can be characterized. For any $j\in\mathbb{N}$, $g_{\left\Vert B\right\Vert ^{j}}$ is bounded on $\mathbb{B}\setminus\left\{ 0\right\} $ if $\left(g_{\left\Vert B\left(d\right)\right\Vert ^{j}\mid dd'},g_{\left\Vert B\left(d'\right)\right\Vert ^{j}\mid dd'},g_{dd'}\right)$ are bounded on $\mathbb{B}$, for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$. The following result shows the asymptotic normality of the EB and MC-EL estimators.

Nonlinearity bias

This section carries out a higher-order analysis of the MC-EL and EB estimators $\widehat{\vartheta}_{p}^{\mathsf{mc}}$ and $\widehat{\vartheta}_{p}^{\mathsf{eb}}$. We apply the quadratic stochastic expansion (newey_smith_2004_higher) to the estimator and write it as the sum of a quadratic function of centered sample averages and a remainder term of a smaller order of magnitude. E.g., for CCFT's estimator, using the expression ((ref)), we simply write \[ \widehat{\vartheta}_{Y,p}^{\mathsf{CCFT}}=\frac{1}{nh}\sum_{i}\widehat{W}_{p,i}\left(Y_{i}-Z_{i}^{\top}\gamma_{Y}\right)-\left(\frac{1}{nh}\sum_{i}\widehat{W}_{p,i}Z_{i}\right)^{\top}\left(\widehat{\gamma}_{Y}^{\mathsf{CCFT}}-\gamma_{Y}\right). \] The first-order asymptotic analysis is based on the linear term $\left(nh\right)^{-1}\sum_{i}\widehat{W}_{p,i}\left(Y_{i}-Z_{i}^{\top}\gamma_{Y}\right)$. Using the second term on the right hand side of above equation and replacing $\widehat{\gamma}_{Y}^{\mathsf{CCFT}}-\gamma_{Y}$ by its linearization, we extract the quadratic terms. For $\widehat{\vartheta}_{p}^{\mathsf{mc}}$ and $\widehat{\vartheta}_{p}^{\mathsf{eb}}$, more complicated derivations are needed.

In our nonparametric context, we write the leading (linear and quadratic) terms as the sum of the first-order stochastic variability term, the first-order smoothing bias term, the second-order smoothing bias term, the second-order stochastic variability term, and a (smoothing) bias-variability interaction term. The first-order (second-order) stochastic variability term is a linear (quadratic) function of centralized sample averages. The first-order stochastic variability term is approximately distributed as $\mathrm{N}\left(0,\mathscr{V}_{p}/\left(nh\right)\right)$. The first-order smoothing bias has a leading term given by $\mathscr{B}_{p}^{\mathsf{eb}}h^{p+1}$ (or $\mathscr{B}_{p}^{\mathsf{mc}}h^{p+1}$). The expectation of the second-order stochastic variability term is referred to as the nonlinearity bias.\footnote{Such a bias is referred to as “higher-order bias” by newey_smith_2004_higher and Graham2012. We use terminology similar to Cattaneo2013 to distinguish such a bias incurred by (second-order) stochastic variability from smoothing bias in our nonparametric context.} The following theorem provides an asymptotic representation for the nonlinearity bias.

Balancing over functions in linear sieve spaces

The EB approach can incorporate information from not only the covariate balance conditions imposed on $Z$ but also on those imposed on functions of $Z.$ This improves efficiency relative to CCFT and can achieve the best attainable asymptotic variance derived in Noack2021. Let $\mathcal{Z}\subseteq\mathbb{R}^{d_{z}}$ denote the support of $Z$. We assume that $\mathcal{Z}$ is compact. Let $\left(b_{1},...,b_{k},...\right)$ be approximating basis functions defined on $\mathcal{Z}$ (typically, $b_{1}=1$). Denote $\rho\coloneqq\left(b_{1},...,b_{k}\right)^{\top}$. We assume that $k=k_{n}$ increases with the sample size $n$. For notational simplicity, we suppress the dependence of $k$ on $n$ and also the dependence of $\rho$ on $k$. Examples of such basis functions commonly used in econometrics include algebraic polynomials (and their transformations), trigonometric polynomials, and B-spline functions, among others. See, e.g., chen2007large and Belloni:2015fo for more details. Consider the following problem of balancing functions in the linear sieve space $\mathcal{M}_{k}\coloneqq\left\{ \rho^{\top}\gamma:\gamma\in\mathbb{R}^{k}\right\} $ and define the EB weights $\left(w_{1}^{\mathsf{sieve}},...,w_{n}^{\mathsf{sieve}}\right)$ as the solution to

The constraint in ((ref)) imposes the balancing constraint that $\sum_{i}w_{i}

=0$ for all $f\in\mathcal{M}_{k}$ for the balancing weights. Then, by the Lagrangian multiplier method, we get the optimal weights and the associated Lagrangian multiplier $\left(w_{i}^{\mathsf{sieve}},\lambda_{p}^{\mathsf{sieve}}\right)$ defined by the right-hand sides of (\ref{eq:entropy balancing weight definition}) and (\ref{eq:lambda_hat}) with $\bar{Z}_{i}$ replaced by $\rho\left(Z_{i}\right)$. Define the sieve EB estimator by

The weak continuity and predeterminedness assumptions (Assumption (ref)(e,f)) allow us to use only the margins in constructing the balancing weights. Now, in order to use functions in broader classes for further efficiency gain, we essentially need CCFT's strong predeterminedness assumption $F_{Z\left(1\right)\mid10}\left(\cdot\mid0\right)=F_{Z\left(0\right)\mid10}\left(\cdot\mid0\right)$ (see Section III of CCFT for discussion), where $F_{Z\left(j\right)\mid dd'}\left(\cdot\mid x\right)$ denotes the conditional cumulative distribution function (CDF) of $Z\left(j\right)$ given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$ and $F_{Z\mid X}$ denotes the conditional CDF of $Z$ given $X$. In addition, we need to replace Assumption (ref)(e) with the stronger assumption that the conditional distributions of $\left(Z\left(d\right),Z\left(d'\right)\right)$ given $\left(D\left(1\right)=d,D\left(0\right)=d',X=x\right)$ change smoothly around the threshold $x=0$. Under the strong continuity and predeterminedness assumptions, $F_{Z\mid X}$ changes smoothly ($\mathrm{lim}_{x\downarrow0}F_{Z\mid X}\left(z\mid x\right)=\mathrm{lim}_{x\uparrow0}F_{Z\mid X}\left(z\mid x\right)$ for all $z\in\mathcal{Z}$). Then $\left(g_{f\left(Z\left(d\right)\right)\mid dd'},g_{f\left(Z\left(d'\right)\right)\mid dd'}\right)$ are continuous at 0 for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$ if $f\in\mathcal{M}_{k}$ satisfies some mild conditions.\footnote{If for all $\left(j,d,d',x\right)\in\left\{ 0,1\right\} ^{3}\times\left[\underline{x},\overline{x}\right]$, the conditional distribution of $Z\left(j\right)$ given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$ admits a density $f_{Z\left(j\right)\mid dd'}\left(\cdot\mid x\right)$ with respect to the $\sigma$-finite dominating measure $\nu$ with $\nu\left(\mathcal{Z}\right)<\infty$ such that $f_{Z\left(j\right)\mid dd'}\left(z\mid\cdot\right)$ is continuous at 0 for all $z\in\mathcal{Z}$, we can write $g_{f\left(Z\left(j\right)\right)\mid dd'}\left(x\right)=\int f\left(z\right)f_{Z\left(j\right)\mid dd'}\left(z\mid x\right)\nu\left(\mathrm{d}z\right)$. If for all $x$ in an open neighborhood of 0, $f_{Z\left(j\right)\mid dd'}\left(\cdot\mid x\right)$ is uniformly bounded, this condition is satisfied if $\int_{\mathcal{Z}}\left|f\right|\mathrm{d}\nu<\infty$.} And the strong predeterminedness assumption implies that $g_{f\left(Z\left(1\right)\right)\mid10}\left(0\right)=g_{f\left(Z\left(0\right)\right)\mid10}\left(0\right)$. It is clear from $g_{f\left(Z\right)}\left(x\right)=\sum_{d,d'}g_{f\left(Z\right)\mid dd'}\left(x\right)g_{dd'}\left(x\right)$ that the covariate balance condition for $f\left(Z\right)$ (i.e., $\mu_{f\left(Z\right),+}=\mu_{f\left(Z\right),-}$) is fulfilled.

Theorem (ref) below shows that the sieve EB estimator given in ((ref)) achieves further efficiency gain relative to CCFT's estimator and our EB estimator in ((ref)). Interestingly, the asymptotic variance of the sieve EB estimator coincides with Noack2021's best attainable asymptotic variance of their LP estimator in which a flexible function is subtracted from the dependent variable. It is clear from ((ref)) that one can write CCFT's estimator for $\mu_{Y,\dagger}$ ($\mu_{D,\dagger}$) as a standard LP regression estimator using $Y_{i}-Z_{i}^{\top}\widehat{\gamma}_{Y}^{\mathsf{CCFT}}$ ($Y_{i}-Z_{i}^{\top}\widehat{\gamma}_{D}^{\mathsf{CCFT}}$) as the dependent variable. Noack2021 consider replacing the linear adjustment $Z_{i}^{\top}\widehat{\gamma}_{Y}^{\mathsf{CCFT}}$ in ((ref)) with a nonlinear transformation of the baseline covariates $Z_{i}$. Let $\left(\eta_{Y},\eta_{D}\right)$ be a real-valued adjustment functions defined on $\mathcal{Z}$ and we consider the standard LP regression estimator $\widehat{\vartheta}_{p}\left(\eta_{Y},\eta_{D}\right)$ using $Y_{i}-\eta_{Y}\left(Z_{i}\right)$ (or $D_{i}-\eta_{D}\left(Z_{i}\right)$) as the dependent variable in ((ref)). Such an estimator is consistent if $\mu_{\eta_{Y}\left(Z\right),+}=\mu_{\eta_{Y}\left(Z\right),-}$ and $\mu_{\eta_{D}\left(Z\right),+}=\mu_{\eta_{D}\left(Z\right),-}$. See Noack2021 for more discussion. Denote $\mu_{+}^{*}\left(z\right)\coloneqq\mathrm{lim}_{x\downarrow0}g_{M\mid ZX}\left(z,x\right)$, $\mu_{-}^{*}\left(z\right)\coloneqq\mathrm{lim}_{x\uparrow0}g_{M\mid ZX}\left(z,x\right)$, $\eta^{*}\left(z\right)\coloneqq\left(\mu_{+}^{*}\left(z\right)+\mu_{-}^{*}\left(z\right)\right)/2$, and $\epsilon^{*}\coloneqq M-\eta^{*}\left(Z\right)$. Noack2021 show that under some mild conditions on $\left(\eta_{Y},\eta_{D}\right)$, the asymptotic variance of $\widehat{\vartheta}_{p}\left(\eta_{Y},\eta_{D}\right)$ cannot be smaller than the best attainable asymptotic variance $\mathscr{V}_{p}^{\mathsf{opt}}\coloneqq\omega_{p}^{0,2}\sigma_{\mathsf{opt}}^{2}/\left(\varphi\mu_{D,\dagger}^{2}\right)$, where $\sigma_{\mathsf{opt}}^{2}\coloneqq\mathrm{Var}_{\mid0^{\pm}}\left[\epsilon^{*}\right]$.\footnote{ Let $\eta^{\dagger}\coloneqq\eta_{Y}-\vartheta\eta_{D}$. $\mathscr{V}_{p}^{\mathsf{opt}}$ is an asymptotic variance lower bound for all $\widehat{\vartheta}_{p}\left(\eta_{Y},\eta_{D}\right)$ with adjustment functions $\left(\eta_{Y},\eta_{D}\right)$ fulfilling the condition that $\mathrm{Cov}_{\mid0^{+}}\left[\eta^{\dagger}\left(Z\right),\mu_{\mathrm{s}}^{*}\left(Z\right)\right]=\mathrm{Cov}_{\mid0^{-}}\left[\eta^{\dagger}\left(Z\right),\mu_{\mathrm{s}}^{*}\left(Z\right)\right]$ for $\mathrm{s}\in\left\{ -,+\right\} $ and $\mathrm{Var}_{\mid0^{+}}\left[\eta^{\dagger}\left(Z\right)-\eta^{*}\left(Z\right)\right]$ and $\mathrm{Var}_{\mid0^{-}}\left[\eta^{\dagger}\left(Z\right)-\eta^{*}\left(Z\right)\right]$ exist. Under the conditions imposed on the densities in Footnote (ref) , this assumption is satisfied, if $\int_{\mathcal{Z}}\left(\eta^{\dagger}\right)^{2}\mathrm{d}\nu<\infty$ and $\int_{\mathcal{Z}}\left(\mu_{\mathrm{s}}^{*}\right)^{2}\mathrm{d}\nu<\infty$ for $\mathrm{s}\in\left\{ -,+\right\} $. Noack2021 show how to construct estimators that attain the optimal asymptotic variance.}

To show that ((ref)) is asymptotically normally distributed with the asymptotic variance $\mathscr{V}_{p}^{\mathsf{opt}}$, we impose the following assumption on the distribution of the observed variables, which we invoke directly in the proof of Theorem (ref). $\mathbb{B}$ is defined in Assumption (ref). Let $\bar{B}_{\eta}\coloneqq\left(Y,D,\eta^{*}\left(Z\right)\right)^{\top}$ and let $f_{X\mid Z}$ be the conditional PDF of $X$ given $Z$.

Sufficient and easy-to-interpret conditions can be imposed on the population distribution of the latent variables to guarantee that Assumption (ref) holds. Under the strong predeterminedness assumption $F_{Z\left(1\right)\mid10}\left(\cdot\mid0\right)=F_{Z\left(0\right)\mid10}\left(\cdot\mid0\right)$, (a) is satisfied if the assumptions on the densities in Footnote (ref) hold and $\eta^{*}$ satisfies the integrability condition in Footnote (ref). (b) is a mild regularity condition similar to Noack2021.\footnote{Under the existence of the densities of the latent variables defined in Footnote (ref), the conditional distribution of $Z$ given $X=x$ admits a density $f_{Z\mid X}\left(\cdot\mid x\right)$ with respect to $\nu$ as a mixture:

The assumption that $f_{Z\mid X}\left(z\mid\cdot\right)$ is $\left(p+1\right)$-times continuously differentiable on $\mathbb{B}\setminus\left\{ 0\right\} $ for all $z\in\mathcal{Z}$ and $\left(\partial/\partial x\right)^{j}f_{Z\mid X}\left(\cdot\mid x\right)$ is uniformly bounded for all $\left(x,j\right)\in\left(\mathbb{B}\setminus\left\{ 0\right\} \right)\times\left\{ 0,1,...,p+1\right\} $ is satisfied if for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$, (1) $\left(f_{Z\left(d\right)\mid dd'}\left(z\mid\cdot\right),f_{Z\left(d'\right)\mid dd'}\left(z\mid\cdot\right)\right)$ are $\left(p+1\right)$-times continuously differentiable on $\mathbb{B}$ with uniformly (in $z\in\mathcal{Z}$) bounded derivatives; (2) $g_{dd'}$ is $\left(p+1\right)$-times continuously differentiable on $\mathbb{B}$ with bounded derivatives. Then under these assumptions, we have $g_{\varepsilon_{n}\left(Z\right)}^{\left(p+1\right)}\left(x\right)=\int\varepsilon_{n}\left(z\right)\left(\left(\partial/\partial x\right)^{p+1}f_{Z\mid X}\left(z\mid x\right)\right)\nu\left(\mathrm{d}z\right)$ and Part (b) holds. } (c,d) are similar to Assumption (ref)(a,b). By ((ref)), these are satisfied under suitable smoothness assumptions on $\left(g_{B_{\eta}\left(d\right)\mid dd'},g_{B_{\eta}\left(d'\right)\mid dd'},g_{dd'}\right)$ and $\left(g_{B_{\eta}\left(d\right)^{\otimes2}\mid dd'},g_{B_{\eta}\left(d'\right)^{\otimes2}\mid dd'}\right)$ for $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$, where $B_{\eta}\left(d\right)\coloneqq\left(Y\left(d\right),\eta^{*}\left(Z\left(d\right)\right)\right)^{\top}$. The first part of Assumption (ref)(e) is satisfied if $f_{X\mid Z}\left(\cdot\mid z\right)$ are differentiable on $\mathbb{B}\setminus\left\{ 0\right\} $ with uniformly (in $z\in\mathcal{Z}$) bounded derivatives.\footnote{By the Bayes theorem, we can show that the first part is satisfied if the assumptions discussed in Footnote (ref) hold and $f_{X}$ is continuously differentiable on $\mathbb{B}$ with uniformly continuous derivatives.} The second part of Assumption (ref)(e) is satisfied if the conditional PDF $f_{M\mid ZX}\left(y\mid z,\cdot\right)$ of $M$ given $\left(Z,X\right)$ is differentiable on $\mathbb{B}\setminus\left\{ 0\right\} $ with derivatives that satisfy some dominance and integrability condition.\footnote{Assume for simplicity that the support $\mathcal{Y}$ of $Y$ is bounded. Let $f_{Y\left(j\right)Z\left(j\right)\mid dd'}\left(\cdot\mid x\right)$ denote the conditional joint density of $\left(Y\left(j\right),Z\left(j\right)\right)$ given $\left(D\left(1\right),D\left(0\right),X\right)=\left(d,d',x\right)$, for $\left(j,d,d',x\right)\in\left\{ 0,1\right\} ^{3}\times\left[\underline{x},\overline{x}\right]$. Then we can write the conditional joint density $f_{MZ\mid X}$ of $\left(M,Z\right)$ given $X$ as a mixture similar to ((ref)) and write $f_{M\mid ZX}=f_{MZ\mid X}/f_{Z\mid X}$. It is clear that the second part is satisfied, if for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$ (1) $\left(f_{Y\left(d\right)Z\left(d\right)\mid dd'}\left(y,z\mid\cdot\right),f_{Y\left(d\right)Z\left(d\right)\mid dd'}\left(y,z\mid\cdot\right)\right)$ are differentiable on $\mathbb{B}$ with uniformly (in $\left(y,z\right)\in\mathcal{Y}\times\mathcal{Z}$) bounded derivatives and similar assumptions hold for $\left(f_{Z\left(d\right)\mid dd'}\left(z\mid\cdot\right),f_{Z\left(d'\right)\mid dd'}\left(z\mid\cdot\right)\right)$ and $g_{dd'}$; (2) $\left(f_{Z\left(d\right)\mid dd'}\left(\cdot\mid\cdot\right),f_{Z\left(d'\right)\mid dd'}\left(\cdot\mid\cdot\right)\right)$ are bounded away from zero on $\mathcal{Z}\times\mathbb{B}$ and a similar assumption holds for $g_{dd'}$.} Assumption (ref)(f) is similar to Assumption (ref)(d). Under Assumption (ref)(c), it is satisfied as long as $\mathrm{Var}\left[Y\left(d\right)-\eta^{*}\left(Z\left(d\right)\right)\mid X=0,\mathsf{co}\right]>0$ for $d\in\left\{ 0,1\right\} $.

We also impose the following assumption on the basis functions. For notational simplicity, let $P\coloneqq\rho\left(Z\right)$.

(a) is satisfied by all commonly used basis functions, as long as the assumption in Footnote (ref) is fulfilled and $\left(\mu_{+}^{*},\mu_{-}^{*}\right)$ satisfy the integrability condition in Footnote (ref). (b) is a standard regularity condition imposing restrictions on the collinearity of the basis functions for which mild sufficient conditions are available (see Belloni:2015fo).\footnote{Part (b) is satisfied if (1) the conditional distribution of $Z$ given $X=x$ admits a Lebesgue density that is uniformly (for all $x\in\mathbb{B}\setminus\left\{ 0\right\} $) bounded above and away from zero; (2) the marginal distribution of $Z$ admits a Lebesgue density that is bounded above and away from zero; (3) the basis functions are orthonormal with respect to the Lebesgue measure. By ((ref)), Condition (1) is satisfied if for all $\left(d,d'\right)\in\left\{ 0,1\right\} ^{2}$, $\left(f_{Z\left(d\right)\mid dd'}\left(\cdot\mid\cdot\right),f_{Z\left(d'\right)\mid dd'}\left(\cdot\mid\cdot\right)\right)$ are bounded above and away from zero on $\mathcal{Z}\times\mathbb{B}$ and a similar assumption holds for $g_{dd'}$.} (c) is analogous to Assumption (ref)(c). (d) is satisfied under (b) and some other mild regularity conditions.\footnote{Under conditions in Footnote (ref), $g_{P}^{\left(p+1\right)}\left(x\right)=\int\rho\left(z\right)\left(\left(\partial/\partial x\right)^{p+1}f_{Z\mid X}\left(z\mid x\right)\right)\nu\left(\mathrm{d}z\right)$. Then, by Jensen's and Cauchy-Schwarz inequalities, \[ \left\Vert g_{P}^{\left(p+1\right)}\left(x\right)\right\Vert ^{2}\leq\mathrm{E}\left[\left(\frac{\left(\partial/\partial x\right)^{p+1}f_{Z\mid X}\left(Z\mid x\right)}{f_{Z\mid X}\left(Z\mid x\right)}\right)^{2}\right]\cdot\mathrm{E}\left[\left\Vert P\right\Vert ^{2}\mid X=x\right]. \] Under Part (b), $\mathrm{E}\left[\left\Vert P\right\Vert ^{2}\mid X=x\right]=\mathrm{tr}\left(\mathrm{E}\left[PP^{\top}\mid X=x\right]\right)\leq k\cdot\overline{\sigma}$. Part (d) holds if the first term is bounded. This holds if (1) $f_{Z\mid X}$ is bounded away from zero on $\mathcal{Z}\times\left(\mathbb{B}\setminus\left\{ 0\right\} \right)$; (2) $\left(\partial/\partial x\right)^{p+1}f_{Z\mid X}\left(\cdot\mid x\right)$ is uniformly bounded for all $x\in\mathbb{B}\setminus\left\{ 0\right\} $. Sufficient conditions for these assumptions are discussed in Footnotes (ref) and (ref).}

Let $\alpha_{k}\coloneqq\mathrm{inf}_{\gamma\in\mathbb{R}^{k}}\left\Vert \eta^{*}-\rho^{\top}\gamma\right\Vert _{\infty}$ and $\beta_{k}\coloneqq\mathrm{sup}_{x\in\mathcal{Z}}\left\Vert \rho\left(z\right)\right\Vert $. Bounds for $\alpha_{k}$ under commonly used basis functions are available from the approximation theory. E.g., if we take $\left(b_{1},...,b_{k},...\right)$ to be the algebraic polynomials and $\eta^{*}$ is $s$-smooth (see, e.g., chen2007large for its definition), then $\alpha_{k}$ is bounded by $k^{-s/d_{z}}$, up to a constant. Bounds for $\beta_{k}$ are also available in the literature for commonly used basis functions. For the algebraic polynomials, $\beta_{k}$ is bounded by $k$ up to a constant. If $\left(b_{1},...,b_{k},...\right)$ are B-splines, then an upper bound is $\sqrt{k}$. See chen2007large and Belloni:2015fo for results for other basis functions. In the statement of the following theorem, we impose Assumptions (ref) and (ref) in place of Assumption (ref).

Likelihood ratio based inference

In this section, we consider inference using the likelihood ratio statistics. Denote $M_{i}\left(\theta\right)\coloneqq Y_{i}-\theta D_{i}$, $U_{i}\left(\theta\right)\coloneqq\left(M_{i}\left(\theta\right),\bar{Z}_{i}^{\top}\right)^{\top}$ and $U_{i}\coloneqq U_{i}\left(\vartheta\right)$ for notational simplicity. Let $\tau\in\left(0,1\right)$ be the significance level. Let $F_{\chi_{1}^{2}}$ and $f_{\chi_{1}^{2}}$ denote the CDF and the PDF of a $\chi_{1}^{2}$ ($\chi^{2}$ with one degree of freedom) random variable respectively. Let $c_{\tau}\coloneqq F_{\chi_{1}^{2}}^{-1}\left(1-\tau\right)$ be the $\left(1-\tau\right)$ quantile of the $\chi_{1}^{2}$ distribution. The standard EL ratio statistic is given by $\mathit{LR}_{p}\left(\theta\mid h\right)\coloneqq2n\left(\ell_{p}^{\mathsf{mc}}\left(\theta\mid h\right)-\ell_{p}^{\mathsf{mc}}\left(\widehat{\vartheta}_{p}^{\mathsf{mc}}\mid h\right)\right)$, which is a function of $\theta$. An EL confidence set for $\vartheta$ with nominal coverage probability $1-\tau$ is $\mathit{CS}_{p,\tau}\left(h\right)\coloneqq\left\{ \theta:\mathit{LR}_{p}\left(\theta\mid h\right)\leq c_{\tau}\right\} $.\footnote{For fuzzy RD, as noack2019bias's method, the EL confidence set avoids a “delta method” argument used by the Wald-type inference of CCFT.} When $p=2$ is taken, our smoothness assumption and construction of $\mathit{CS}_{p,\tau}\left(h\right)$ parallel CCFT in that $\mathit{CS}_{p,\tau}\left(h\right)$ uses the same LP order as CCFT's inference method and Assumption (ref)(a) assumes the same (three-times differentiability) smoothness as CCFT's Assumption SRD.\footnote{CCFT proposes Wald-type inference using their local linear estimator with bias correction and standard errors that take into account estimation of the bias. calonico2014robust show that subtracting the $p$-th order LP estimator by the nonparametric estimator for the leading bias term with the same bandwidth is the same as a $\left(p+1\right)$-th order LP estimator. CCFT's bias-corrected local linear estimator (with common bandwidths) is the same as a local quadratic regression estimator.} By the Lagrangian multiplier method and strong duality, for fixed $\theta$,

By similar derivations as those in Section (ref),

Computation of $\mathit{LR}_{p}\left(\theta\mid h\right)$ only requires solving convex optimization problems. The right hand side of the second equality in ((ref)) can be $\infty$ in the “no solution” scenario discussed in Section (ref). If our algorithm finds a solution for the maximization problem in ((ref)), then we proceed to compute ((ref)) for fixed $\theta$ using a similar algorithm. The right hand side of ((ref)) is $\infty$ if the origin is not in the interior of the convex hull of $\left\{ W_{p,1}U_{1}\left(\theta\right),...,W_{p,n}U_{n}\left(\theta\right)\right\} $. In this case, the Newton algorithm would return a very large value for $\ell_{p}^{\mathsf{mc}}\left(\theta\mid h\right)$ and $\theta$ is excluded from the confidence set. We have the following result on the shape of $\mathit{CS}_{p,\tau}\left(h\right)$.

In the rest of this section, we give several large-sample properties of the EL inference method. Section (ref) establishes uniform-in-bandwidth (first-order) validity of the EL confidence set. Sections (ref) and (ref) are devoted to second-order properties. Section (ref) shows the distributional expansion for the likelihood ratio and proposes a simple analytical correction to improve coverage accuracy. Section (ref) considers a scenario in which covariate balance fails to hold and analyze the sensitivity of the coverage accuracy to this assumption. We derive the distributional expansion under local perturbation to the covariate balance condition.

Uniform-in-bandwidth Wilks theorem

The following theorem parallels the main result of AK and is a substantial extension of the standard Wilks theorem, which states that $\mathit{LR}_{p}\left(\vartheta\mid h\right)\rightarrow_{d}\chi_{1}^{2}$. Our result incorporates covariates and accommodates unbounded outcomes. The proof techniques we use differ from those employed by AK. Let $\ell^{\infty}\left(\left[1,\overline{h}/\underline{h}\right]\right)$ denote the space of all bounded real-valued functions defined on $\left[1,\overline{h}/\underline{h}\right]$ endowed with the sup-norm. Let $\mathbb{H}\coloneqq\left[\underline{h},\overline{h}\right]$ be a compact bandwidth set where $\underline{h}>0$ and $\overline{h}>0$ ($\underline{h}<\overline{h}$) are bandwidths that depend on the sample size.\footnote{As the main result of AK, Theorem (ref) assumes deterministic upper and lower bounds. Let $\left(\underline{h}^{*},\overline{h}^{*}\right)$ denote some deterministic bounds that some data-dependent bounds $\left(\underline{h},\overline{h}\right)$ capture. As argued by AK, the conclusion of Theorem (ref) still holds under data-dependent bounds if the orders of $\overline{h}/\overline{h}^{*}-1$ and $\underline{h}/\underline{h}^{*}-1$ are sufficiently small and $\left(\underline{h}^{*},\overline{h}^{*}\right)$ satisfy the assumptions of Theorem (ref).}

Analytical correction

This section provides coverage expansions of the EL confidence sets. Similar to calonico2018optimal, Theorem (ref) below considers two scenarios under the given smoothness assumption (Assumption (ref)(a)). The first scenario uses the LP order $p$ so that the leading bias term in the coverage error of the confidence set $\mathit{CS}_{p,\tau}\left(h\right)$ can be characterized. The second scenario exhausts the smoothness by setting LP order to $p+1$. The smoothing bias, in this case, is of a smaller order $O\left(h^{p+1+\mathfrak{h}}\right)$ but its leading term can not be explicitly characterized. The following mild assumption on the kernel function is used when establishing the validity of the Edgeworth expansions in the proofs of Theorems (ref) and (ref).

Since $K\left(\cdot\right)$ is assumed to be symmetric, an analogous property holds for $\left(1,\mathcal{K}_{p;-},\mathcal{K}_{p;-}^{2},\mathcal{K}_{p;-}^{3}\right)$ as functions on $\left(-1,0\right)$ under this assumption. It is clear that the assumption is satisfied if $\mathcal{K}_{p;+}$ is a non-constant polynomial on $\left[-1,1\right]$. If $p\geq1$, this condition is satisfied if $K\left(\cdot\right)$ is any of the commonly used kernel functions (triangular, biweight, triweight, etc.) including the uniform kernel.\footnote{Suppose that $p=1$ and $K$ is the uniform kernel, i.e., $K\left(t\right)=\mathbbm{1}\left(\left|t\right|\leq1\right)/2$. Then, by simple calculation, $\mathcal{K}_{p;+}\left(t\right)=\left(4-6t\right)\mathbbm{1}\left(\left|t\right|\leq1\right)$ and $\mathcal{K}_{p;-}\left(t\right)=\left(4+6t\right)\mathbbm{1}\left(\left|t\right|\leq1\right)$.} Denote $\Xi\coloneqq\mu_{UU^{\top},\pm}^{-1}$, $\Psi_{1}^{\mathsf{kl}}\coloneqq\mathrm{tr}\left(\Xi\cdot\mu_{U^{\left(\mathsf{k}\right)}U^{\left(\mathsf{l}\right)}UU^{\top},\pm}\right)$ and $\Psi_{2}^{\mathsf{kl}}\coloneqq\mathrm{tr}\left(\Xi\cdot\mu_{U^{\left(\mathsf{k}\right)}UU^{\top},\dagger}\Xi\cdot\mu_{U^{\left(\mathsf{l}\right)}UU^{\top},\dagger}\right)$. Let

Let $\mathscr{V}_{p}^{\ddagger}$ be defined by the same formula with $U$ replaced by $\bar{Z}$ and the range changed to $1,...,d_{z}+1$ accordingly. Let $\mathscr{V}_{p}^{\mathsf{LR}}\coloneqq\mathscr{V}_{p}^{\dagger}-\mathscr{V}_{p}^{\ddagger}$ and $\mathscr{B}_{p}^{\mathsf{LR}}\coloneqq\left(\mathscr{B}_{p}^{\mathsf{mc}}\right)^{2}/\mathscr{V}_{p}$.

Now we provide distributional expansions for both $\mathit{LR}_{p}\left(\vartheta\mid h\right)$ and $\mathit{LR}_{p+1}\left(\vartheta\mid h\right)$. Asymptotic expansions of the coverage probabilities follow from these results (e.g, $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p,\tau}\left(h\right)\right]=\mathrm{Pr}\left[\mathit{LR}_{p}\left(\vartheta\mid h\right)\leq c_{\tau}\right]$). The proof uses the method of calonico2018coverage and calculations in chen2007second.

Local imbalance

This section shows that the coverage performance of the EL confidence set is maintained even if the covariate balance assumption is slightly violated, a scenario we call “local imbalance”. Specially, we assume that the observed covariates $Z$ are subject to data contamination (measurement errors) that occurs after treatment. The contaminated covariates may not satisfy the predeterminedness assumption and can be drawn from some perturbed probability law that generates a slight imbalance (kitamura2013robustness). On the other hand, the genuine but unobserved predetermined covariates $Z^{\star}\in\mathbb{R}^{d_{z}}$, which typically affects $Y\left(d\right)$, still satisfy the balance condition. The continuity of $\left(g_{Y\left(d\right)\mid dd'},g_{Y\left(d'\right)\mid dd'}\right)$ remains to hold. In other words, the imbalance is caused by measurement errors that are known to be excluded from the data-generating process of $Y\left(d\right)$. In this case, the standard RD estimand $\vartheta$, which confidence sets try to cover, remains to identify a causal parameter of interest.\footnote{frolich2019including consider a different scenario where conditionally on $Z$, the continuity (exclusion) assumption is satisfied. This essentially assumes that there are no unobserved variables that both affect the potential outcomes and are affected by $I$. frolich2019including show that the RD LATE is still identifiable (under additional assumptions), but smoothing over $Z$ is required for estimation.}

Formally, let $\zeta\in\mathbb{R}^{d_{z}}$ denote the measurement errors realized after treatment. The measurement error $\zeta$ is nonclassical in the sense that it relates to $\left(D,X,Z^{\star}\right)$. Let $\left(Z^{\star}\left(1\right),Z^{\star}\left(0\right),\zeta\left(1\right),\zeta\left(0\right)\right)$ be potential covariates and measurement errors.\footnote{The structural model representation in Footnote (ref) can be extended to $Y=g\left(D,X,Z^{\star},\epsilon\right)$, $D=h\left(I,X,\eta\right)$, $Z^{\star}=m\left(D,X,\xi\right)$ and $\zeta=q\left(D,Z^{\star},X,\nu\right)$ for some unknown functions $\left(g,h,m,q\right)$ and unobserved disturbances $\left(\epsilon,\eta,\xi,\nu\right)$.} The contaminated potential covariates $\left(Z\left(1\right),Z\left(0\right)\right)$ are generated by $Z\left(d\right)=Z^{\star}\left(d\right)+\zeta\left(d\right)$ for $d=0,1$. And the observed contaminated covariates are $Z=D\cdot Z\left(1\right)+\left(1-D\right)Z\left(0\right)$. We assume that the true covariates satisfy the “predeterminedness” assumption $g_{Z^{\star}\left(1\right)\mid10}\left(0\right)=g_{Z^{\star}\left(0\right)\mid10}\left(0\right)$, but the measurement errors fail to satisfy it. As a result, local imbalance in essence assumes that $g_{Z\left(1\right)\mid10}\left(0\right)-g_{Z\left(0\right)\mid10}\left(0\right)$ approaches 0 at the rate of $\left(nh\right)^{-1/2}$. We are interested in the coverage probability $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)\right]$, which is expected to have a limit in $\left(0,1-\tau\right)$ and thus captures the phenomenon that covariate imbalance results in undercoverage.We set the bandwidth to $h\asymp n^{-1/\left(p+2\right)}$ as discussed in Remark (ref). Let $l_{n}\coloneqq n^{-\left(p+1\right)/\left(2p+4\right)}\asymp\left(nh\right)^{-1/2}$. The following assumption formalizes local imbalance.

Part (a) essentially assumes that no other variables depending on $I$ affect $\left(Z^{\star}\left(d\right),\zeta\left(d\right)\right)$. Under Assumption (ref)(a,b,c,d), the standard RD estimand still identifies the RD LATE. (a,b) imply that the true covariates that may affect $Y\left(d\right)$ still satisfy the balance condition $\mu_{Z^{\star},+}=\mu_{Z^{\star},-}$. (c) assumes that the RD LATE on $\zeta$ is $\delta\cdot l_{n}$, which generates the local imbalance in the observed covariate: $\mu_{Z,+}-\mu_{Z,-}\asymp\left(nh\right)^{-1/2}$. By using local asymptotic analysis, we analyze the performance of our EL confidence set under such a local imbalance condition, which is similar to using locally misspecified moment conditions in the sense of armstrong2021sensitivity. Our result differs from armstrong2021sensitivity and focuses on the coverage performance of the confidence set when $\delta$ is close to 0.\footnote{The approach of armstrong2021sensitivity specifies a set in which $\delta$ possibly lies and then adjusts the critical value to take into account the maximal misspecification bias. We take a very different approach in this paper.}

Let $N\coloneqq Z-\left(\delta\cdot l_{n}\right)D=D\cdot N\left(1\right)+\left(1-D\right)N\left(0\right)$, where $N\left(1\right)\coloneqq Z^{\star}\left(1\right)+\zeta\left(1\right)-\delta\cdot l_{n}$ and $N\left(0\right)\coloneqq Z^{\star}\left(0\right)+\zeta\left(0\right)$. It now follows that $g_{N\left(1\right)\mid10}\left(0\right)=g_{N\left(0\right)\mid10}\left(0\right)$ and $\mu_{N,+}=\mu_{N,-}$. Let $\gamma_{N}\coloneqq\left(\mathrm{Var}_{\mid0^{\pm}}\left[N\right]\right)^{-1}\mathrm{Cov}_{\mid0^{\pm}}\left[N,M\right]$ and $\mathscr{V}_{N}\coloneqq\left(\omega_{p+1}^{0,2}\mathrm{Var}_{\mid0^{\pm}}\left[M-N^{\top}\gamma_{N}\right]\right)/\left(\varphi\mu_{D,\dagger}^{2}\right)$. For simplicity, we assume that the distribution of $N$ does not vary with $n$.\footnote{E.g., this holds if the measurement errors are the following form: $\zeta\left(1\right)=\delta\cdot l_{n}+\zeta_{1}$ and $\zeta\left(0\right)=\zeta_{0}$, for some zero-mean $\left(\zeta_{0},\zeta_{1}\right)$ that are independent of other variables in the model. Relaxation of this assumption requires more complicated arguments and suitable modification of the assumptions.} CCFT shows that the covariate-adjusted estimator is inconsistent and the confidence interval fails to have asymptotically correct coverage probability under “global imbalance” $\mu_{Z,+}\neq\mu_{Z,-}$. Under local imbalance in Assumption (ref), CCFT's estimator and the generalized EB estimators are still consistent.\footnote{We can show that $\widehat{\gamma}_{Y}^{\mathsf{CCFT}}$ (see Section 6 of the online supplement of CCFT for its expression) in the representation ((ref)) converges in probability to $\left(\mathrm{Var}_{\mid0^{\pm}}\left[N\right]\right)^{-1}\mathrm{Cov}_{\mid0^{\pm}}\left[N,Y\right]$. Then, since $\mu_{N,+}=\mu_{N,-}$, we have \[ \widehat{\vartheta}_{Y,p}^{\mathsf{CCFT}}=\frac{1}{nh}\sum_{i}\widehat{W}_{p,i}\left(Y_{i}-N_{i}^{\top}\widehat{\gamma}_{Y}^{\mathsf{CCFT}}\right)-\left(\delta\cdot l_{n}\right)\left(\frac{1}{nh}\sum_{i}\widehat{W}_{p,i}D_{i}\right)=\mu_{Y,\dagger}+o_{p}\left(1\right). \] Under global imbalance, the EB estimator has a probabilistic limit different from that of CCFT's estimator (see Lemma 1 of CCFT). Neither of them is equal to $\vartheta$.} Inference suffers from the undercoverage problem, since the coverage probabilities of the confidence sets (CCFT's or the EL) converge to a limit in $\left(0,1-\tau\right)$.

We now consider $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)\right]$ as a function of $\delta$ under local imbalance. A measure of sensitivity of the coverage accuracy to local imbalance (i.e., how the coverage probability drops relative to that under $\delta=0$) is given by the slope of $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)\right]$ as a function of $\delta$ at $\delta=0$. We extend Theorem (ref) and derive a two-term asymptotic expansion $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)\right]=R\left(\delta\right)+o\left(l_{n}\right)$, where $R\left(\delta\right)$ is the sum of the leading terms as an approximation to $\mathrm{Pr}\left[\vartheta\in\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)\right]$ in finite samples. We show that $R\left(0\right)=1-\tau$ and the gradient $\nabla R\left(\delta\right)\coloneqq\left(\partial/\partial\delta\right)R\left(\delta\right)$ at $\delta=0$ is equal to 0, so that $R\left(\delta\right)$ is locally constant around $\delta=0$.

Let $F\left(\cdot\mid\iota\right)$ denote the CDF of a $\chi_{1}^{2}\left(\iota\right)$ (non-central $\chi^{2}$ with one degree of freedom and non-centrality parameter $\iota\geq0$) random variable. Let $F^{\left(k\right)}\left(x\mid\iota\right)\coloneqq\left(\partial/\partial\iota\right)^{k}F\left(x\mid\iota\right)$ be the $k$-times partial derivative of $F\left(x\mid\iota\right)$ with respect to $\iota$.

Covariate-adjusted estimation of the treatment effect derivative

The EB approach for covariate adjustment applies to parameters of interest other than the standard RD LATE parameter. This section applies EB to covariate-adjusted estimation of the treatment effect derivative (TED). To focus on the main ideas, we consider the sharp design first. Dong2015 propose using the TED defined as $\left.\left(\mathrm{d}/\mathrm{d}x\right)\mathrm{E}\left[Y\left(1\right)-Y\left(0\right)\mid X=x\right]\right|_{x=0}$ for evaluating the external validity of RD. A large TED suggests that the LATE would be quite different if the score changes slightly, raising more concern about external validity. The researcher can check whether the RD LATE is likely to have external validity by testing for zero TED. Under the assumption that $g_{Y\left(d\right)}$ is continuously differentiable on a neighborhood of 0 (Dong2015), the TED is identified: $\left.\left(\mathrm{d}/\mathrm{d}x\right)\mathrm{E}\left[Y\left(1\right)-Y\left(0\right)\mid X=x\right]\right|_{x=0}=\pi_{\mathsf{srd}}\coloneqq\mu_{Y,+}^{\left(1\right)}-\mu_{Y,-}^{\left(1\right)}$. This section proposes an EB estimator for the TED. An inferential procedure and standard errors can be found in Section S10 of the online supplement.

Let $\dot{W}_{p;-,i}$ be defined by the right-hand side of ((ref)) with $\mathrm{e}_{p+1,1}^{\top}$ replaced by $\mathrm{e}_{p+1,2}^{\top}$. Similarly, we define $\dot{W}_{p;+,i}$ and in addition, let $\dot{W}_{p,i}=\dot{W}_{p;+,i}-\dot{W}_{p;-,i}$. In our notation, the standard LP estimator proposed in Dong2015 for the TED is given by $\widehat{\pi}_{p}^{\mathsf{lp}}\coloneqq\left(nh^{2}\right)^{-1}\sum_{i}\dot{W}_{p,i}Y_{i}$ ($p\geq2$). As the EB estimator $\widehat{\vartheta}_{p}^{\mathsf{eb}}$ for the RD LATE proposed in Section (ref), the EB-based TED estimator $\widehat{\pi}_{p}^{\mathsf{eb}}$ with covariate adjustment can also be obtained by replacing the uniform weights with the EB weights $\left(w_{1}^{\mathsf{eb}},...,w_{n}^{\mathsf{eb}}\right)$ defined by ((ref)):

The above construction illustrates the convenience of EB-based covariate adjustment: One can start with the standard estimator (without covariates) for a parameter of interest in an RD-related context and then replace its standard uniform weights with the EB weights. The EB weights are computed using the covariates only, and are independent of the standard estimator. Such an adjustment strategy also works straightforwardly in other RD-related settings, such as the (nonlinear) estimators of xu2017regression,Xu2018 in the scenarios with limited outcome variables.

Let $\gamma_{\mathsf{ted}}\coloneqq\left(\mathrm{Var}_{\mid0^{\pm}}\left[Z\right]\right)^{-1}\left(\mathrm{Cov}_{\mid0^{+}}\left[Z,Y\right]-\mathrm{Cov}_{\mid0^{-}}\left[Z,Y\right]\right)$ and $\mathcal{\dot{K}}_{p;\mathrm{s}}\left(t\right)\coloneqq\mathrm{e}_{p+1,2}^{\top}\mathrm{V}_{p;\mathrm{s}}^{-1}r_{p}\left(t\right)K\left(t\right)$, for $\mathrm{s}\in\left\{ -,+\right\} $. One can easily verify that $\mathcal{\dot{K}}_{p;+}\left(t\right)=-\mathcal{\dot{K}}_{p;-}\left(-t\right)$ and $\int_{-1}^{0}\mathcal{K}_{p;-}\left(t\right)\dot{\mathcal{K}}_{p;-}\left(t\right)\mathrm{d}t=-\varpi_{p}$, where $\varpi_{p}\coloneqq\int_{0}^{1}\mathcal{K}_{p;+}\left(t\right)\dot{\mathcal{K}}_{p;+}\left(t\right)\mathrm{d}t$. Also denote $\dot{\omega}_{p;+}^{j,k}\coloneqq\int_{0}^{1}t^{j}\dot{\mathcal{K}}_{p;+}^{k}\left(t\right)\mathrm{d}t$ and $\dot{\omega}_{p;-}^{j,k}\coloneqq\int_{-1}^{0}t^{j}\dot{\mathcal{K}}_{p;-}^{k}\left(t\right)\mathrm{d}t$. It can be checked that $\dot{\omega}_{p;+}^{0,2}=\dot{\omega}_{p;-}^{0,2}\eqqcolon\dot{\omega}_{p}^{0,2}$. The following theorem shows the asymptotic normality of the EB estimator.

Monte Carlo simulations

We conduct simulations to evaluate the finite sample performance of the proposed EL-based inference methods for sharp RD designs with covariates. The data-generating process (DGP) of the outcome variable $Y_{i}$, the score $X_{i}$ and the first covariate $Z_{i}^{\left(1\right)}$ is based on the simulation design of CCFT. The incorporation of additional covariates $Z_{i}^{\left(2\right)},...,Z_{i}^{\left(l\right)}$ follows that of arai2021regression. Let

and

where with the coefficients $\gamma_{l}=0.22$, $\gamma_{r}=0.28$, $\theta_{l}=1.06$ and $\theta_{r}=0.61$, all following CCFT. Then, $Y_{i}=\mu_{y}\left(X_{i},Z_{i}^{\left(1\right)}\right)+\sum_{j=2}^{l}\pi^{j-1}Z_{i}^{\left(j\right)}+\varepsilon_{y,i}$ and $Z_{i}^{\left(1\right)}=\mu_{z}\left(X_{i}\right)+\varepsilon_{z,i}$. Error terms $\left(\varepsilon_{y,i},\varepsilon_{z,i}\right)$ are bivariate normal with mean $0$, standard deviation $1$ and correlation coefficient $\rho=0.269.$ Additional covariates $\left(Z_{i}^{\left(2\right)},...,Z_{i}^{\left(l\right)}\right)$ have a multivariate normal distribution with mean zero and covariance matrix given by $\mathrm{Cov}\left[Z_{i}^{\left(j\right)},Z_{i}^{\left(k\right)}\right]=0.5^{\left|j-k\right|},$ for all $j,k\geq2$. We take $\pi=0.2$. We consider three scenarios with $l=0,2,4$, corresponding to the total number of covariates $d_{z}=l+1$ being $1,3,5$. CCFT uses local linear regression with bias correction, equivalent to local quadratic regression. Our EL approach parallels CCFT in that the degree of the LP is set to be $p=2$. The sample sizes are $n=1000,2000$. The number of Monte Carlo replications is $5000$.

Table (ref) presents the bias, root mean square error (RMSE) of the MC-EL estimator $\widehat{\vartheta}_{p}^{\mathsf{mc}}$ defined in Section (ref), as well as the empirical coverage probability and the average length of the EL confidence sets $\mathit{CS}_{p,\tau}^{\mathsf{bc}}\left(h\right)$ and $\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h\right)$ defined in Remark (ref). Following Remark (ref), we select a bandwidth of the form $h=H\cdot n^{-1/\left(p+2\right)}$, replace $H$ with a consistent estimator $\widehat{H}$, and use the bandwidth $\widehat{h}\coloneqq\widehat{H}\cdot n^{-1/\left(p+2\right)}$. calonico2018optimal propose an approach that takes the estimated AMSE optimal bandwidth and rescales it to make it obey the coverage optimal rate (see Section IV(C) of CCFT). One choice of bandwidth $\widehat{h}$ is to follow this approach and use CCFT's bandwidth, denoted as CCFT in Table (ref). CCFT's bandwidth is computed from $\mathtt{R}$ function $\mathtt{rdrobust}$ with the options $\mathtt{p=1}$, $\mathtt{rho=1}$, and $\mathtt{bwselect}=$“$\mathtt{cerrd}$”.\footnote{The rate of CCFT's bandwidth is $n^{-1/4}$, which matches the rate of $\widehat{h}=\widehat{H}\cdot n^{-1/\left(p+2\right)}$ with $p=2$.} Another simpler choice is a rescaled rule of thumb (ROT) bandwidth that uses the constant part $\widehat{H}$ computed according to Hansen2021's ROT bandwidth. For comparison, Table (ref) also includes results from CCFT's method that uses the CCFT bandwidth and restricts $\rho=h/b=1$, where $b$ stands for the pilot bandwidth used for bias estimation. Table (ref) shows that both EL and CCFT approaches perform well for estimation and inference. A closer look reveals that EL with $p=2$ and using the rescaled ROT bandwidth yields similar bias and RMSE compared with CCFT, but slightly better coverage (especially for $d_{z}=3$ and $5$), and shorter confidence intervals. On the other hand, EL that uses CCFT's bandwidth, which amounts to half of the ROT bandwidth, yields smaller bias but larger RMSE and longer confidence intervals. In particular, the length of $\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(\widehat{h}\right)$ is longer than those of other confidence sets. In sum, all the methods we consider deliver satisfactory finite-sample performances.\footnote{The EL confidence intervals are also well-centered. E.g., the average center (across all $5,000$ simulation replications) of the EL intervals with ROT bandwidth ($\mathrm{EL}_{p}$ in Table 1), $d_{z}=5$ and $n=1,000$ is only $0.0015$ away from the true treatment parameter ($0.0494$).} Computing the EB weights (for the point estimator) and the EL likelihood ratio statistic (for the confidence set) only requires solving convex optimization problems (corresponding to the “inner loop” in the standard EL computation) and thus is very fast. E.g., computing the row $\mathrm{EL}_{p}$ in Table 1 with $d_{z}=5$, $n=1,000$ and CCFT bandwidth costs $0.06$ to $1.08$ seconds for one replication, with the average computation time per replication about $0.23$ second on an Intel Core i7 processor with 32 GB of RAM.

We also examine how the coverage performance of EL and CCFT confidence sets changes when the covariate balance condition is slightly violated. We consider the case with one covariate ($d_{z}=1$). The data-generating process for $\left(Y_{i},X_{i},Z_{i}^{\left(1\right)}\right)$ remains the same but the incorporated covariate is given by $\tilde{Z}_{i}^{\left(1\right)}\coloneqq Z_{i}^{\left(1\right)}+\mathbbm{1}\left(X_{i}<0\right)\delta$, so that the local covariate imbalance is measured by the perturbation $\delta$. Figure (ref) plots the simulated coverage probabilities of the EL and CCFT confidence sets as a function of $\delta\in\left[-0.3,0.3\right]$. We observe that the coverage probability of $\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(\widehat{h}\right)$ is less sensitive to the change of $\delta$, which parallels the discussion in Remark (ref).

We then investigate the performance of the EB approach in the covariate-adjusted estimation of TED. We consider the case with one covariate ($d_{z}=1$) and modify the coefficients $\gamma_{l}$, $\gamma_{r}$, $\theta_{l}$ and $\theta_{r}$ in the design in order to highlight two features of covariate-adjusted estimation of TED. First, the magnitude of efficiency gain from incorporating the single covariate is determined by $|\gamma_{l}-\gamma_{r}|$. We choose $\gamma_{l}=3$ and $\gamma_{r}=0$ to highlight the efficiency contribution of the covariate adjustment. Second, note that $\theta_{r}-\theta_{l}$ corresponds to $\mu_{Z,+}^{(1)}-\mu_{Z,-}^{(1)}$, which is required to be zero for CCFT's augmented regression estimator of TED (Remark (ref)). More specifically, the asymptotic bias of CCFT's estimator is proportional to $|\theta_{l}-\theta_{r}|$. We set $\theta_{l}=3$ and $\theta_{r}=0$ to highlight such a bias. Table (ref) reports the finite-sample performances of three TED estimators: the EB estimator with $p=2$ given by ((ref)) and its confidence interval constructed following the procedure in Section S10 in our online supplement, CCFT's TED estimator with $p=2$ in Remark (ref), and the standard local quadratic (LQ) TED estimator (Dong2015) without using covariate information, all three methods using calonico2014robust's bandwidth for the first derivative computed from the $\mathtt{R}$ function $\mathtt{rdrobust}$.\footnote{The rate of CCT's bandwidth here is still $n^{-1/4}$, which is also the CO rate for Wald-type inference on the TED. See calonico2018optimal.} As expected, Table (ref) shows that the CCFT's TED estimator leads to a substantial bias and undercoverage for TED inference, given that $\mu_{Z,+}^{(1)}-\mu_{Z,-}^{(1)}$ is away from zero. This problem can be solved by the EB estimator, which incorporates the correctly specified covariate balance condition $\mu_{Z,+}=\mu_{Z,-}$ rather than the misspecified condition $\mu_{Z,+}^{(1)}=\mu_{Z,-}^{(1)}$ . The standard LQ TED estimator without covariates remains valid but has a larger RMSE and yields a longer confidence interval than those from the EB estimator, which reflects the efficiency gain of EB from covariate adjustment.

Empirical illustration: Finnish municipal election data

We apply our estimation/inference method to analyze the individual incumbent advantage in Finnish municipal elections, which was first studied by hyytinen2018does. The outcome variable $Y$ indicates whether the candidate is elected in an election, and the score $X$ is the vote share margin in the previous election. Table (ref) presents the RD LATE point estimate $\widehat{\vartheta}$, the $p$-value for testing the null hypothesis $\vartheta=0$, the $95\%$ confidence intervals (CI), and the CI length. The first row of Table (ref) presents the standard LQ regression estimator that ignores the covariates. Then, we incorporate four covariates $Z$: candidates' age, gender, age squared, and $\mathrm{age}\times\mathrm{gender}$. EL estimation and inference use CCT's bandwidth ($h_{\mathsf{CCT}}=0.396$) and the rescaled ROT bandwidth (equal to $2.917$). The last row of Table (ref) reproduces the “experiment benchmark” reported originally by hyytinen2018does(see their Table 2, Column 4, the $p$-value is imputed by us).\footnote{The dataset includes 1351 candidates “for whom the (previous) electoral outcome was determined via random seat assignment due to ties in vote counts” (hyytinen2018does), which constitutes an experiment benchmark to evaluate the credibility of the RD treatment effect estimated from the non-experimental data (candidates with previous electoral ties are excluded from the RD sample).} Apparently, all RD estimates, with or without covariates, are small in magnitude and statistically insignificant, which agrees with the finding in the experiment benchmark. By comparing the covariate-adjusted estimates (EL and CCFT) with the standard LQ regression without covariates, we see that incorporating covariates helps to reduce the CI length for four out of five confidence intervals, except for $\mathit{CS}_{p+1,\tau}^{\mathsf{bc}}\left(h_{\mathsf{CCT}}\right)$. Among them, the EL confidence set $\mathit{CS}_{p,\tau}^{\mathsf{bc}}\left(h_{\mathsf{CCT}}\right)$ that uses the same bandwidth as the standard LQ regression is $7.2\%$ shorter than the standard method and is $5.5\%$ shorter than CCFT. Here, the efficiency improvement is moderate, probably because the election outcome is only weakly correlated with age and gender.

We then conduct a sensitivity analysis of the EL-based covariate-adjusted inference with respect to the bandwidth choice by plotting the confidence band (Remarks (ref) and (ref)). We consider the continuous range of bandwidths $h\in\left[\underline{h},\overline{h}\right]$ with the lower bound $\underline{h}=h_{\mathsf{CCT}}/3\approx0.13$ and the upper bound $\overline{h}=h_{\mathsf{CCT}}\times2\approx0.78$. The rate of $h_{\mathsf{CCT}}$ is $n^{-1/4}$, which satisfies the conditions for $\underline{h}$ and $\overline{h}$ in Theorem (ref). Using the $\mathtt{R}$ package $\mathtt{BWSnooping}$, we calculate the snooping corrected critical value $\mathrm{2.413}^{2}$ for the triangular kernel and bandwidth ratio $\overline{h}/\underline{h}=6$. In Figure (ref), the solid (or dotted) lines correspond to a $95\%$ uniform (or pointwise) confidence band. For small bandwidth (say, less than $0.2$), the uniform confidence band is wide. However, as long as the bandwidth is not so small, the confidence band appears stable. Moreover, the confidence band includes zero over the entire bandwidth range, demonstrating the robustness of the finding of no incumbency advantage with respect to the bandwidth choice.

Lastly, we evaluate the external validity by testing the null hypothesis that the TED is zero. It will tell us whether the RD estimate, which by design only applies to the “local” incumbents whose previous vote share margin resides at the $0$ cutoff, can be applied to incumbents whose previous vote share margins are slightly higher than $0$. When estimating the TED, we maintain the usual covariate balance condition $\mu_{Z,+}=\mu_{Z,-}$ but do not impose the balance condition $\mu_{Z,+}^{(1)}=\mu_{Z,-}^{(1)}$ for the derivatives, so the CCFT's augmented regression estimator for TED is not a proper choice, as discussed in Remark (ref). Our EB method gives a point estimate of TED equal to $-0.631$, and a $p$-value for testing a zero TED equal to $0.064$. In comparison, the standard estimate (without covariates) of TED is $-0.634$ with the $p$-value equal to $0.023$.\footnote{Both EB and the standard estimates use the CCT bandwidth for the first derivative, which equals $0.462$.} Therefore, both methods raise the concern of external validity of applying the RD estimate to incumbents with share margins above $0$, as the treatment effect is likely to significantly decrease in response to a marginal increase in the score.

Conclusion and further discussion

This paper proposes a balancing approach to covariate adjustment for RD. The covariate balance condition can be viewed as over-identifying restrictions, which the EB estimator incorporates when formulated as an EL estimator. By establishing the first-order equivalence between the EB estimator and CCFT's regression estimator, we show that the efficiency gain can be attributed to incorporating covariate balance as side information.

The EB problem ((ref)) can be cast in a more general framework under which several extensions can be considered. The construction follows BenMichael2021. Consider the following imbalance measure $\mathrm{imbalance}_{\mathcal{M}}\left(w_{1},...,w_{n}\right)\coloneqq\mathrm{sup}_{f\in\mathcal{M}}\left|\sum_{i}w_{i}\widehat{W}_{p,i}f\left(Z_{i}\right)\right|$ with respect to a function space $\mathcal{M}$. Let $m:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}\cup\left\{ \infty\right\} $ be an increasing and convex function. Let $\mathrm{complexity}\left(w_{1},...,w_{n}\right)$ denote some complexity (or dispersion) measure of the weights. Consider the following risk minimization problem similar to BenMichael2021:

for some tuning parameter $\varsigma>0$. Denote $\mathcal{M}_{0}\coloneqq\left\{ \mathbb{R}^{d_{z}}\ni z\mapsto a+z^{\top}b:\left|a\right|+\sum_{j=1}^{d_{z}}\left|b^{\left(j\right)}\right|\leq1\right\} $. We may also take the Cressie-Read divergence $D_{\varrho}\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)$ defined by ((ref)) as a complexity measure. It is clear that under $\mathrm{complexity}\left(w_{1},...,w_{n}\right)=D_{\varrho}\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)$ the generalized balancing problem ((ref)) can be written in the form ((ref)) with $\mathcal{M}=\mathcal{M}_{0}$ (see BenMichael2021) and $m\left(\cdot\right)$ taken to be $\mathbb{R}_{+}\ni x\mapsto\infty\cdot\mathbbm{1}\left(x>0\right)$, so that exact balance is required.\footnote{The risk minimization problem can now be written as

Clearly, the optimal weights do not depend on the choice of $\varsigma$ in this case. Relaxation of the “exact balance” constraint by using a strictly positive threshold (see BenMichael2021) is also straightforward.} The sieve balancing problem ((ref)) is also of the form ((ref)) with $\mathrm{complexity}\left(w_{1},...,w_{n}\right)=KL\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)$ and $\mathcal{M}$ taken to be the broader sieve space.

An alternative balancing scheme similar to Hirshberg2021 is based on solving ((ref)) with $\mathcal{M}$ taken to be the sieve space, $m\left(\cdot\right)$ taken to be $x\mapsto x^{2}$ and $\mathrm{complexity}\left(w_{1},...,w_{n}\right)$ taken to be the “square Euclidean” divergence given by $D_{-2}\left(w_{1},...,w_{n}\parallel1/n,...,1/n\right)$. Then, we expect to find a dual characterization of the optimal weights by using results from Hirshberg2021. An asymptotic normality result similar to Theorem (ref) is expected to hold under a suitable choice of tuning parameters $\left(h,k,\varsigma\right)$. With $\mathcal{M}$ taken to be a ball in a Reproducing Kernel Hilbert Space (RKHS), we get a balancing scheme similar to Kallus2020,Wong2017 (see Wong2017). We also expect an asymptotic normality result similar to Theorem ((ref)) holds under a suitable choice of the three tuning parameters ($\left(h,\varsigma\right)$ and the radius of the ball) and the assumption that the “optimal adjustment function” $\eta^{*}$ (see Section (ref)) lies in the RKHS.

Our EB approach avoids the selection of the additional tuning parameter. Another advantage is the favorable second-order properties developed in the EL literature carry over to our proposed method. These include a small nonlinearity bias of the point estimator and a simple analytical correction to improve coverage accuracy for the confidence set. We also show a uniform-in-bandwidth Wilks theorem, which can be used for sensitivity analysis and robust inference along the lines of AK. We also derive the distributional expansion for the EL ratio statistics under the local imbalance condition and analyze the sensitivity of the coverage performance to the balance assumption. Lastly, we demonstrate that our approach can address previously unsolved covariate adjustment problems in RD by deriving an EB-based covariate-adjusted estimator for the TED. We also expect the large-deviation optimality results for EL (e.g., otsu_2010_bahadur) to carry over. In the presence of high-dimensional covariates (arai2021regression,Kreiss2022), resorting to the dual characterization ((ref)), we apply appropriate penalization in ((ref)) (Chang2018) to reduce the effective number of covariates. Properties of the penalized EB are left for future investigation.