Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,058 characters · 8 sections · 103 citation commands
PLRD: Partially Linear Regression Discontinuity Inference
\allowdisplaybreaks \sloppy {5pt}
The fraction \begingroup \footnote{We are grateful to Evan Munro for sharing results on the behavior of existing methods for RDDs under generative adversarial simulation following athey2021using. We are also grateful to Tim Armstrong, P. Aronow, Max Farrell, Paul Goldsmith-Pinkham, Michal Koles\'ar, Christoph Rothe and Roc\'io Titiunik for insightful conversations and for comments on a previous version of this paper. This work was supported by the National Science Foundation under grant SES-2242876, the Office of Naval Research under grants N00014-17-1-2131 and N00014-19-1-2468, and a gift from Amazon. Software and replication files are available at \url{https://github.com/ghoshadi/plrd}.} \addtocounter{footnote}{-1} \endgroup
of empirical NBER working papers using regression discontinuity designs (RDDs) increased from nearly zero in 1995 to between 15% and 20% by 2015 currie2020technology.\footnote{goldsmith2024tracking maintains up-to-date numbers at \url{https://paulgp.com/econlit-pipeline}.} Despite the ubiquity of RDDs in applied work, however, there are still open questions as to econometric best practices for inference in RDDs. In particular, we find that widely used methods systematically fail to achieve nominal coverage.
We assess the performance of three commonly used methods by evaluating them on simulations calibrated to nine high-profile applications of RDDs. First, we consider conventional local linear regression \citep*{Hahn-et-al-2001} with bandwidth chosen as in imbens2012optimal. We also consider two more recent proposals by \citet*{CCT2014robustCI} and ArmstrongKolesar2018, as implemented in the R functions rdrobust and RDHonest respectively. These methods are representative of the current state of the art for inference in RDDs \citep*{stommes2023reliability}, and both offer automated end-to-end confidence interval constructions that can directly be applied to raw data. The method rdrobust by CCT2014robustCI is a dominant approach to regression-discontinuity inference in current econometric practice.\footnote{This is the only general purpose regression discontinuity package listed in the “Econometrics Task View” for highly used R packages on CRAN (\url{https://cran.r-project.org/web/views/Econometrics.html}); and both the R calonico2015rdrobust and Stata calonico2017rdrobust versions of the function are widely used and the original paper CCT2014robustCI cited is over 4500 times as of now.} The method starts by running a local linear regression as recommended by \citet*{Hahn-et-al-2001}, and then applies a 2nd-order bias correction to mitigate curvature bias.\footnote{The method \texttt{rdrobust} offers two options for building confidence intervals: A basic bias-corrected method which inherits the standard error estimator from the original local linear regression (and is justified by asymptotic arguments), and a robust method that incorporates a finite-sample variance inflation to account for the bias correction. The latter is the recommended one in the paper and this is the one we consider here.} Meanwhile, the method \texttt{RDHonest} by ArmstrongKolesar2018 explicitly inflates the width of local linear regression confidence intervals to accommodate worst-case plausible curvature bias (rather than trying to correct for this bias).\footnote{In principle, the method \texttt{RDHonest} asks the user to specify a worst-case bound on the curvature of the underlying conditional response functions up front. However, the method also has a default fall-back option where curvature bounds are derived using the rule-of-thumb of ArmstrongKolesar2020; in our experiments we use this default.}
To ensure that the evaluations of these methods are credible, we follow \citet*{athey2021using} and use generative adversarial neural networks goodfellow2014generative to generate simulation instances that are representative of a set of seven seminal regression discontinuity studies.\footnote{Two of the studies have two outcomes, so we have a total of nine applications.} The generative adversarial learning approach yields simulators whose data-generating distribution closely mimics the moments and idiosyncrasies of the original datasets; see (ref) for an illustration. The generated simulators do not use explicit parametric models and do not enforce any built-in smoothness assumptions, thus bringing us towards a neutral comparison ground for evaluating methods that were motivated under different assumptions.
This semi-synthetic simulation setup is an important part of our argument. It comes with two qualifications. First, it relies on the availability of datasets where regression discontinuity methods are applicable. We apply it to nine such data sets, but we encourage researchers to apply it to additional ones. Second, it relies on our generative adversarial nets yielding realistic data generating processes given the provided datasets. Well-known results in non-parametric inference imply that, fixing a sample size and given any dataset where standard regression-discontinuity estimators work well, one can construct another adversarial dataset that is statistically indistinguishable from the first but where the same estimators fail to get coverage low1997nonparametric. Our approach to evaluation implicitly relies on an assumption that the real datasets we are calibrated to aren't adversarial in this sense. Furthermore, our estimation process involves some choices of tuning parameters for which different researchers may make different choices; in our experiments we seek to mitigate concerns of researcher bias by using a single set of tuning parameters across all datasets.
We summarize results for 95% confidence intervals in (ref); see (ref) for details on the simulation design. A first observation is that the conventional local linear regression intervals systematically undercover the true parameter. The method rdrobust achieves better---although not perfect---coverage, with nominal rates in {4 out of the 9 settings}. However, it does so at the expense of wider confidence intervals than conventional local linear regression {(18% to 36% wider in our applications)}. RDHonest in turn achieves reliable coverage throughout all designs, but at the cost of still-wider intervals: The median ratio between the width of these intervals and {the rdrobust ones varies between 1.06 to 1.76 in our simulations, and between 1.25 to 2.11 relative to the conventional ones}. Furthermore, these confidence intervals are often excessively conservative, thus suggesting that some of this additional width is not strictly needed.\footnote{In several settings, RDHonest achieves 98% coverage or above---and this may come at a substantial cost in width. For example, in the classical setting with no bias adjustment, Gaussian 98% confidence intervals are 18.7% wider than the conventional 95% intervals.} These qualitative findings do not appear to be specific to our simulation designs; for example, in (ref) we observe similar coverage properties when applying the method to pure noise with no relationship between the outcome and the running variable.
These semi-synthetic simulation findings leave open a pragmatic question: Is it possible to design automated confidence intervals that reliably achieve nominal coverage in realistic settings, while remaining competitive with the conventional intervals in terms of interval width? We emphasize here the “realistic” qualifier. Even if it is not feasible to achieve uniformly better properties, it may be that most real-world applications are sufficiently smooth to allow for nominal coverage and narrower confidence intervals than the existing methods.
The main contribution of this paper is a method, partially linear regression discontinuity inference (PLRD), that provides an affirmative answer to this question---at least in the context of the numerical settings we consider. As described further below, algorithmically, the PLRD estimator adapts the method of ArmstrongKolesar2018 to accommodate a partially linear treatment effect model; conceptually, however, our full approach fits within the econometric framework of CCT2014robustCI. We implement our method in the package plrd for R r2023, which provides end-to-end functionality for regression discontinuity inference without relying on the user to specify hard-to-choose tuning parameters. As illustrated in (ref), our proposed plrd method achieves nominal and near-exact coverage throughout, and in nearly every instance yields the shortest intervals among those with valid coverage.
Consider the sharp RDD formalized using potential outcomes as in IL2008review. For each unit $i = 1, \, \ldots, \, n$, we observe IID-sampled pairs $(X_i,\, Y_i)$, where $X_i \in \mathbb{R}$ is the running variable and $Y_i \in \mathbb{R}$ is the outcome of interest. Following the sharp discontinuity design, we assume there exists a cutoff $c \in \mathbb{R}$ such that treatment is assigned as $W_i = \text{$\mathbf{1}$}{\left\{X_i\ge c\right\}}$. We posit potential outcomes $\{Y_i(0), \, Y_i(1)\}$ such that $Y_i = Y_i(W_i)$, and our goal is to estimate the conditional average treatment effect at the cutoff, $\tau_c:=\mu_{1}(c)-\mu_{0}(c)$ where $\mu_{w}(x) := \mathbb{E}[Y_i(w)|X_i = x]$.
It is well known that, provided the $\mu_{w}(x)$ are continuous in $x$ and that the running variable has continuous support around $c$, we can identify our target parameter as a discontinuity in the conditional response surface at $X_i = c$ Hahn-et-al-2001:
In order to obtain practical estimation or inference results for $\tau_c$, though, one needs to make further regularity assumptions on the $\mu_{w}(x)$; our PLRD method relies on the following two assumptions.
(ref) is a standard smoothness assumption; this is exactly the same as the main smoothness assumption used by CCT2014robustCI to justify rdrobust. (ref) implies that the treatment effect has a simpler dependence on the running variable than the baseline effect $\mu_{0}(x)$. This assumption is well in line with assumptions often made in the literature on heterogeneous treatment effect estimation wager2024causal; however, it has not typically been used in the RDD literature to date, which is why we emphasize this assumption in naming our method. Given these two assumptions, PLRD proceeds as follows:
Our full procedure---including some extra algorithmic details such as cross fitting---is detailed as (ref). As a pragmatic safeguard against failures of (ref), we start our procedure with a hypothesis test, and if we reject this assumption with high confidence then we proceed with a variant of our method that only relies on (ref); see (ref) for further discussion.
Our main formal result about the PLRD estimator is that it yields valid confidence intervals under our assumptions in the large-sample limit.\footnote{In our formal results, we assume a continuous running variable with support near the cutoff. In the case of discrete running variables, steps (ref) and (ref) of our procedure remain valid, i.e., given an appropriate $B$ they yield partial identification intervals for $\tau_c$ with a guarantee of the type (ref) optrdd2019,KolesarRothe2018. However, when $X_i$ has discrete support, step (ref) is no longer guaranteed to consistently recover $B$ under (ref) alone; and so further assumptions would be required to justify data-driven detection of $B$.}
We note that our algorithm does in principle require two tuning parameters, namely the pilot bandwidth $\ell$ used to run the initial cubic regression to choose $B$ and the lower bound $\varepsilon$ for $\wh B$. We note, however, that our algorithmic performance is largely insensitive to these choices (and simple defaults appear to provide good performance in practice); and similar objects are also required in CCT2014robustCI and imbens2012optimal. First, regarding $\ell$, this bandwidth can in practice be chosen simply to be “large”: We only need $\ell$ to decay faster than $n^{-1/12}$, whereas the mean-squared error optimal bandwidth for local linear regression in our setting (i.e., under (ref)) scales as $n^{-1/7}$ \citep*{Cheng-Fan-Marron}. In our implementation, we by default simply set $\ell$ to contain the full range of the data, and this choice is what we use in all our experiments.\footnote{This choice mirrors imbens2012optimal who seed their bandwidth-selection procedures by fitting global cubic polynomials to the full dataset; see also CCT2014robustCI.} As such, practitioners should be able to get reliable results out of our method unless the nature of the relationship between $X_i$ and $Y_i$ changes completely as we get far from the cutoff---and in settings like this it would be possible to pick $\ell$ by plotting the data. Meanwhile, $\varepsilon$ will have no effect in large samples as long as we choose it to be small enough that $\varepsilon < |\mu'''(c)|$.\footnote{One interesting question is what happens at the critical point where $\mu'''(c) = 0$; as, in this case, our constraint does necessarily bind and we fall back to using $\varepsilon$ as our curvature bound. The challenge is that, when $\mu'''(c) = 0$, our cubic polynomial estimator \smash{$\wh B$} (in step 3 of (ref)) converges to 0, but the exact way it does so is left underspecified under (ref). Our use of the lower bound $\varepsilon$ guarantees validity of our confidence intervals at the expense of some excess conservativeness when $\mu'''(c) = 0$. We note that existing methods, including CCT2014robustCI and imbens2012optimal, also do not offer guidance on how to potentially benefit from “superefficiency” in similar critical cases; i.e., our approach is again in line with current best practices.}
Regression discontinuity designs were originally introduced by thistlethwaite1960regression. The modern econometric literature goes back to Hahn-et-al-2001, who emphasized the identification result (ref) and recommended estimation $\tau_c$ via local linear regression,
where $h_n \rightarrow 0$ is a bandwidth that “localizes” the regression and $K(\cdot)$ is a kernel function with support on $[-1, \, 1]$. In the decade following Hahn-et-al-2001, local linear (or polynomial) regression remained the prevalent approach to inference in RDDs, with confidence intervals produced by applying off-the-shelf tools for heteroskedasticity-robust inference directly to the above regression IL2008review,lee2010regression. A limitation of this approach, however, is that it leads to contradictory guidance on the bandwidth choice under a non-parametric specification as given in, {\it e.g.}, (ref). The accuracy of \smash{$\wh\tau_\texttt{LLR}$} depends delicately on the bandwidth $h_n$, with a first systematic proposal given in imbens2012optimal. For choices of $h_n$ that optimize this accuracy the bias and standard error of \smash{$\wh\tau_\texttt{LLR}$} are of the same order---so off-the-shelf local linear regression inference will not have valid coverage for this choice of $h_n$. The practitioner is thus forced to invoke “undersmoothing”, whereby they voluntarily use a non-accuracy-optimizing choice of $h_n$ to justify validity of their inferential procedure.
In a major advance, CCT2014robustCI and ArmstrongKolesar2018 introduced paradigms for inference in RDDs that formally account for the bias of local linear regression. As discussed in stommes2023reliability, these methods offer material improvements in the credibility of regression-discontinuity analyses over the earlier practice. The main idea in CCT2014robustCI is to estimate and correct for the bias of local linear regression by leveraging higher-order smoothness; we refer to this approach as the “bias-corrected” approach.\footnote{When reporting results using the bias-correct approach, we use robust standard error estimates as recommended by CCT2014robustCI and as implemented in the package rdrobust calonico2015rdrobust,calonico2017rdrobust.} \citet*{CCFT2019} discuss the use of covariates under the bias-corrected approach to RDDs while cattaneo2022regression provide a recent review. In contrast, the “bias-aware” approach developed by ArmstrongKolesar2018 involves conservatively widening intervals to account for worst-case bias (without needing to actually estimate the bias); see KolesarRothe2018, optrdd2019 and noack2024bias for further applications of this paradigm.
Interestingly, not only do the above papers under the bias-corrected vs. bias-aware paradigms explore different methodological strategies for inference in RDDs; they also focus on proving formal results under different formal paradigms. CCT2014robustCI and follow-ups focus entirely on asymptotic validity results---as we do here. They fix a data-generating process, and guarantee that they eventually get coverage for this data-generating process as the sample size gets large (as is also done in (ref)). In contrast, most papers written under the bias-aware paradigm specify a quantitative smoothness class for the $\mu_w(x)$ ({\it e.g.}, all functions with a second derivative bounded by a pre-specified constant $B$), and then seek guarantees that hold in finite samples and uniformly over this whole class. Getting uniform guarantees in finite samples of course requires stronger assumptions than getting pointwise asymptotic guarantees; and for this reason papers associated with the bias-aware paradigm often start by making stronger assumptions than those associated with the bias-corrected paradigm. We emphasize, however, that this is due to the sought formal guarantees; and it's not that bias-aware methods inherently require stronger assumptions to work than bias-corrected methods.
Our proposed method, PLRD, straddles these two literatures. Our confidence intervals are rooted in the bias-aware paradigm of ArmstrongKolesar2018, and we use numerical convex optimization to construct our estimator as in optrdd2019. In contrast, our formal results are focused on pointwise asymptotics, and the assumptions and guarantee type given in our (ref) mirror those in Theorem 1 of CCT2014robustCI. One notable paper that also straddles these literatures in ArmstrongKolesar2020, who propose bias-aware confidence intervals for local linear regression with pointwise asymptotic guarantees. In particular, they show that if one runs local linear regression as in (ref) with a mean-square-error optimal bandwidth choice $h_n$, then under generic conditions we can build valid, bias-aware 95% confidence intervals for $\tau_c$ as $\wh\tau_\texttt{LLR} \pm 2.18 \times \text{standard errors}$ rather than the conventional $\wh\tau_\texttt{LLR} \pm 1.96 \times \text{standard errors}$. We provide a comparison of PLRD to this approach in (ref).
Wrapping up, PLRD can best be understood within the context of earlier developments in the bias-aware paradigm; and in fact {\it steps ((ref))} and {\it ((ref))} of our algorithmic outline fall squarely within the roadmap for minimax linear inference via numerical convex optimization as presented in optrdd2019. However, in a deviation from existing work:
Regarding this last point, the earlier focus on the $M$-bounded curvature class appears to have been largely guided by intuitive simplicity of this class rather than any conceptual or empirical arguments. Given the empirical performance of PLRD on calibrated simulations, we suggest that our $B$-Lipschitz curvature class with linear treatment effects may be a good default choice in many applications. We conduct an ablation analysis to examine the effect of different choices on the behavior of our estimator in (ref).
We now describe in more detail the experiment presented in the introduction ((ref)). As context for this experiment, (ref) shows a replication of a number of high-profile RDD applications using a variety of methods; citations for each application are given in the leftmost column and further description of each benchmark dataset is provided (ref).
A fundamental challenge in assessing methods for causal inference on real-world data is that we typically do not have access to a ground truth we can use for evaluation. In light of this challenge, athey2021using proposed using Wasserstein Generative Adversarial Networks (WGANs) to generate synthetic data that mimic benchmark datasets---and thus obtaining simulation instances that enable us to compare methods in settings that are as realistic as possible. Our WGAN evaluation pipeline involves the following three steps:
As illustrated in Figure (ref), a trained WGAN is able to simulate data whose distribution is very similar to that of the original data; thus suggesting that if a method gets reliable coverage on the simulator we can also trust it on the real data.
{Our WGAN training procedure consists of the following steps. First, apply marginal CDF transforms to the data to convert $(X_i, Y_i)$ to $(U_i, V_i)$, where $U_i=F_x(X_i)$ and $V_i = F_y(Y_i)$. For discrete data we break ties by randomization, and for additional stability we use a Gaussian transformation. We train two generators $\widehat{g}_w(u, z)$ ($w=0,1$), one for treated units and one for controls, that mimic the conditional distribution of the CDF-transformed outcome $V$ given the CDF-transformed running variable $U$. Finally, we invert the CDF-transforms to obtain $(\widetilde{X}_i,\widetilde{Y}_i)$ from the pair $(\widetilde{U}_i, \widetilde{V}_i)$ generated by the trained WGAN.
We calculate the ground truth $\widetilde{\tau}$ for the GAN estimated distribution as follows. We first calculate the cutoff $c'$ under the applied CDF transforms. We then sample random noise $Z_{0,q}, Z_{1,q}\in \mathbb{R}$ for $q=1,\dots,Q$, where $Q$ is a large number ($10^7$ in our simulations), generate CDF-transformed outcomes at the threshold $c'$, and invert them: $$\widetilde{\tau}:=Q^{-1}\sum_{q=1}^Q \left(F_y^{-1}(\widehat{g}_1(c', Z_{1,q}))-F_y^{-1}(\widehat{g}_0(c',Z_{0,q}))\right).$$ We generate $N=10^4$ simulated datasets indexed by $i=1,\dots,N$, and compute confidence intervals $[a_{i,k},b_{i,k}]$, for each method $k$. The empirical coverage for the $k$-th method is given by $\widehat{p}_k=N^{-1}\sum_{i=1}^n \mathbf{1}\{a_{i,k}\le \widetilde{\tau}\le b_{i,k}\}$. (ref) shows the obtained results. }
For completeness, we also benchmark the PLRD estimator against other existing methods on some standard simulation designs. First, we consider a pure noise setting where we generate the running variables as $X_i \sim \text{Unif}([-1, \, 1])$, the outcomes as $Y_i \sim \mathcal{N}(0,\, 1)$, and set treatment to $W_i = 1(\{X_i \geq 0\})$. Next, we consider simulation settings used in CCT2014robustCI. In each of these settings, we draw the running variable $X_i\sim 2 \, \text{Beta}(2, 4)-1$, and the outcome is given by $Y_i = m(X_i)+ \varepsilon_i$, where $\varepsilon_i$ are mean-zero Gaussian noise with variance $0.1295^2$. The form of $m(\cdot)$ varies for the different simulation settings; see CCT2014robustCI for details on Settings 1 through 3. Setting 4 uses the function $m=m_\text{quad}$ as in imbens2012optimal. We note that the linear treatment effect condition ((ref)) is not satisfied in several of these designs, so we expect our pragmatic safeguard (given in (ref)) to play an important role in delivering robust performance.
(ref) presents the empirical coverage and average width of 95% CIs produced by different methods across these simulation settings. As before, rdrobust falls short of nominal coverage; note that these numbers are in line with those originally reported by CCT2014robustCI. RDHonest delivers valid coverage at the expense of wider intervals, and plrd delivers valid coverage with intervals whose width is competitive with rdrobust.
Our point estimator, i.e., the $\widehat{\gamma}_i$-weighted average of outcomes $Y_i$ with weights as given in (ref) is, ignoring cross-fitting, equivalent to a minimax linear estimator for $\tau_c$ over the class $\mathcal{M}_{\widehat{B}}$ of all functions satisfying (ref) and for which $\mu''_0(x)$ is globally $\widehat{B}$-Lipschitz, i.e.,\footnote{The conditional variance of a linear estimator as in (ref) is \smash{$\sum_{i = 1}^n \widehat{\gamma}_i^2 \sigma_i^2$} with \smash{$\sigma_i^2 = \operatorname{Var}[Y_i \,|\, X_i]$}. If we had \smash{$\sigma_i^2 = \widehat{\sigma}^2$} for all $i = 1, \, \ldots, \, n$, then (ref) would exactly be the minimax linear estimation problem for $\tau_c$. In the case of heteroskedasticity, the term \smash{$\widehat{\sigma}^2 \lVert \gamma \rVert_2^2$} should be interpreted as a variance proxy that yields a practical estimator---and one that still maintains provable validity guarantees under heteroskedasticity; see optrdd2019 for further discussion.}
If $\widehat{B}$ were a constant chosen a-priori and $\mu''_0(x)$ were in fact is globally $\widehat{B}$-Lipschitz, we could use existing results for minimax linear regression-discontinuity inference ArmstrongKolesar2018,optrdd2019 to directly provide guarantees for $\widehat{\tau}_{\mathtt{PLRD}}$. But here, of course, this is not the case: Our choice of $\widehat{B}$ is learned in step 1 of our algorithm, and (ref) only implies local---not global---Lipschitz continuity in $\mu''_0(x)$ .
{The goal of this section is to provide technical tools to address these challenges, culminating in the proof of (ref).} We start by verifying validity of the cubic regression in the first step of our algorithm. The following result is a direct consequence of standard results for local polynomial regression fan1996framework, and implies that $\widehat{B}_{(1)}$ and $\widehat{B}_{(2)}$ in (ref) satisfy \smash{$\widehat{B}_{(k)}\rightarrow_p B$} as $n\to\infty$.
{
}ur core analytic result is the following central limit theorem for our PLRD estimator. Asymptotic validity of the PLRD confidence intervals as claimed in (ref) then follows immediately.
Define $\widehat{s}^2(\wh\gamma)=\sum_{i=1}^n \wh\gamma_i^2\wh\sigma_i^2$ as in (ref). It follows from white1980heteroskedasticity that $\widehat{s}(\wh\gamma)/s(\wh\gamma)\rightarrow_p 1$; see (ref) for details. Next, fix any $\delta\in(0,1)$. In view of (ref) and Slustky's lemma, $\mathbb{P}(A_n)\to 1$ where $A_n:=\{|b(\wh\gamma)|\le \widehat{b}+\delta\widehat{s}(\wh\gamma)\}$. We show in (ref) that for any $s$, $h$, $t$ and $\delta>0$, $$\inf_{|b|\,\le\, t+\delta s}\mathbb{P}(|b+s Z|\le h) \ge \inf_{|b|\,\le\, t}\mathbb{P}(|b+s Z|\le h) -\delta,$$ where $Z\sim \mathcal{N}(0,1)$. Using this with $s=\widehat{s}(\wh\gamma)$, $h=\widehat{h}_\alpha$ and $t=\widehat{b}$, we get, for every $n$,
where the last inequality follows from the definition of $\widehat{h}_\alpha$. On the other hand, (ref) and Slutsky's lemma imply that \smash{$Z_n:=(\wh\tau_\texttt{PLRD}-\tau_c-b(\wh\gamma))/\widehat{s}(\wh\gamma)\stackrel{d}{\longrightarrow} Z\sim \mathcal{N}(0,1)$}. Consequently,
Since $\delta\in(0,1)$ is arbitrary, we are through.
The optimization problem in (ref) is inherently related to its convex dual (see BoydVandenberghe for a definition). We show in (ref) that the dual problem is given by
where $\theta=({\rho},\lambda)$ varies in the set
and
In (ref), we use sample splitting to estimate $B$ (cf. (ref)) and the conditional variance $\sigma^2=\mathbb{E}[{\rm Var}(Y_i\mid X_i)]$. However, for ease of exposition, we assume throughout this section that we have $\widehat{B}\rightarrow_p B$ and $\wh\sigma^2\rightarrow_p\sigma^2$ independently of the data, and $$\wh\theta_n := \wh\theta_n(\wh\sigma^2,\widehat{B})=\mathop{\rm argmin}_{\theta\,\in\,\Theta} \widehat{L}_n(\theta;\wh\sigma^2,\widehat{B}).$$ We can recover the primal solution $(\wh\gamma_n,\widehat{t}_n\,)$ from the dual solution $\wh\theta_n$ as: $$\wh\gamma_{n,i} = -\frac{1}{2\wh\sigma_i^2} \, G(\wh\theta_n; \, X_i-c, W_i),\quad \widehat{t}_n = \frac{1}{2\widehat{B}^2}\widehat{\lambda}_{n,1}.$$ Next, define the population-level dual problem as:
Our proof strategy hinges on establishing fast-enough convergence of \smash{$\wh\theta_n$} to \smash{$\theta_n^*(\sigma^2,B)$} to enable us to read-off large-sample properties of \smash{$\widehat{\tau}_{\mathtt{PLRD}}$} from the solution to the population dual (ref).
Although the empirical minimization problem (ref) looks similar to a moment-estimation problem, it has some technical complications, e.g., (i) the dual parameter $\theta=(\rho,\lambda)$ involves a function $\rho$ and a vector $\lambda$, (ii) the set $\Theta$ in (ref) is not the Cartesian product of a function space and a vector space, and (iii) the population dual problem in (ref) involves the sample size $n$. To tackle these issues, we borrow tools from empirical process theory (see, {\it e.g.,} VW for a textbook treatment). The following result, proven in (ref), characterizes the rate of convergence of the empirical dual optimizer $\wh\theta_n$ to the population dual optimizer $\theta_n^*$; and implies ((ref)) that the weights underlying \smash{$\widehat{\tau}_{\mathtt{PLRD}}$} satisfy a Lindeberg-type condition.
(ref) immediately implies that the central limit theorem claimed in (ref) holds; see the proof of optrdd2019 for details. We delegate showing that $s^2(\wh\gamma)=\mathcal{O}_p(n^{-6/7})$ ((ref)) and $\widehat{b}/{s}(\wh\gamma)=\mathcal{O}_p(1)$ ((ref)) to (ref). Finally, to verify (ref), we note that (ref) yield $b(\widehat{\gamma}) = \sum_{i=1}^n \wh\gamma_i \mu_{W_i}(X_i) = B\sum_{i=1}^n\wh\gamma_i \rho(X_i-c)$, where $$\rho(x-c):=(\mu_0(x)-\mu_0(c)-\mu'_0(c)(x-c)-\mu''_0(c)(x-c)^2/2)/B.$$ It follows from (ref) (with cross-fitted $\widehat{B}$) that $|b(\wh\gamma)|\le B\widehat{t}$ where $\widehat{t}=\widehat{b}/\widehat{B}$. Consequently,
Furthermore, since $\widehat{B}\rightarrow_p B$ and $\widehat{b}/{s}(\wh\gamma)=\mathcal{O}_p(1)$, we have $$\lim_{n\to\infty} \frac{(B/\widehat{B}-1)\widehat{b}}{{s}(\wh\gamma)} =_p 0, $$ and so (ref) implies (ref).