Optimized Inference in Regression Kink Designs

Optimized Inference in Regression Kink Designs

\linespread{1.2}

\linespread{1.5}

\linespread{1.5}

Introduction

Research designs that warrant a causal interpretation of model parameter estimates based on observational data are becoming increasingly popular in the empirical social sciences. The quality of such observational studies is often assessed in terms of how well the assumptions that allow for the estimation of counterfactual quantities can be justified. In this context, methods that exploit discontinuities arising naturally from institutional rules are particularly attractive, as they allow the researcher to be precise about the source of variation utilized in the estimation of the parameter and discuss the assumptions that allow for a causal interpretation.

The Regression Kink Design (RKD) can be employed to estimate causal parameters when the variable of interest is a kinked function of an assignment variable. Analogous to the better-known Regression Discontinuity Design (RDD), which estimates the effect of a variable that changes its level discontinuously at a threshold, kink designs utilize discontinuous changes in the slope of the policy variable, effectively exposing units on each side of the threshold to different incentives. Such kinks arise naturally whenever marginal rates change discontinuously or benefit formulas involve maxima or minima and can be exploited to address important questions that are often difficult to be addressed experimentally, e.g. questions regarding the optimal design of unemployment insurance policies (see card2015effect; landais2015assessing; kolsrud2018optimal).

For example, a common feature of unemployment insurance systems is that unemployment benefits, typically a function of income in some base period, are capped at some maximum level. Consequently, the incentives for reemployment differ at each side of the threshold defined by the cap, and the (local) causal effect of the benefit level on unemployment duration can be estimated by comparing how the duration outcome changes with prior income at each side of the cutoff relative to the size of the kink. In practice, this amounts to estimating the jump in the first derivative of the conditional expectation function (CEF) of the duration outcome given the assignment variable at the kink point, and dividing the estimate by the kink size. If unobserved confounders vary smoothly at the discontinuity point, this corresponds up to scale to the (local) elasticity of unemployment duration with respect to the benefit level, a key parameter in dynamic labor supply models and corresponding welfare optimal benefit formulas.

Since the welfare effects of policy changes can often be expressed in terms of elasticities, kink designs can be used in a sufficient statistics approach to policy evaluation, avoiding the need for parametric assumptions and the estimation of structural primitives of the model, combining credible identification with the ability to make welfare predictions chetty2009sufficient.

A practical challenge for researchers utilizing discontinuity designs is to choose estimation and inference procedures that preserve the credibility of their findings stemming from the nonparametric identification of the causal parameter. The standard approach to inference in regression kink designs relies on local polynomial regression, that is fitting a polynomial model of order $p \geq 1$ using only observations within a prespecified window of length $2h$ around the threshold $c$ by weighted least squares. The theoretical properties of local polynomial estimators have been studied extensively (e.g. ruppert1994multivariate; fan1996local; fan1997local) and it is well understood that if the true CEF differs from a polynomial of order $p$ on $[c-h,c+h]$, the resulting estimator for the kink parameter is generally biased. The two standard approaches to inference in regression kink designs recognize this by combining an asymptotic normality result with an argument that adresses this smoothing bias.

One approach is to choose the bandwidth $h$ sufficiently small such that for a given sample size $n$ the resulting bias is (hopefully) negligible relative to the estimator's standard deviation. This strategy is referred to as undersmoothing in the nonparametric regression literature and uses plug-in estimates of bandwidth sequences that shrink faster than the asymptotic mean squared error (MSE) optimal sequence, eliminating the bias from the asymptotic approximation invoked for inference. In practice, the bandwidth is chosen by multiplying a regularized estimate of the pointwise, that is evaluated at the unknown true CEF, asymptotic MSE optimal bandwidth by $n^{- \delta}$ for some small $\delta > 0 $ (cf. imbens2008regression).

A second approach relies on estimating the leading bias term using a higher order polynomial and constructing the interval around a bias corrected point estimate, taking into account the additional variance introduced by bias estimation in the asymptotic approximation. This approach is referred to as robust bias correction (calonico2014robust) and is implemented using a pilot bandwidth tuned for bias estimation and a regularized estimate of the pointwise asymptotic MSE-optimal bandwidth for confidence interval (CI) construction.

In a recent paper, armstrong2020simple demonstrate that both default approaches to nonparametric confidence interval construction can lead to severe undercoverage in finite samples when implemented using bandwidth selectors justified by pointwise asymptotics, as is common. They attribute this finding to the pointwise statistical guarantees underlying these procedures in general, and in particular to the fact that, irrespective of the global curvature of the CEF, the pointwise asymptotic MSE-optimal bandwidth selector can be arbitrarily large if the $(p+1)$-th derivative of the CEF is close to zero at the threshold, resulting in large bandwidth choices for potentially highly nonlinear functions. While this problem was previously recognized, the regularization terms that are added to prevent the empirical bandwidth selectors from selecting too large bandwidths are sensitive to ad-hoc choices of tuning parameters that drive the finite sample coverage rates of the resulting intervals. As a solution, the authors propose to explicitly restrict the parameter space by placing a bound on the $(p+1)$-th derivative of the unknown CEF and select bandwidths according to a minimax criterion, avoiding the need for regularization and allowing for the construction of intervals that are honest, in the sense that they are valid uniformly over the considered parameter space. This approach is feasible, as the bound allows for the computation of the magnitude of the exact worst-case smoothing bias of the local polynomial estimator over the class of functions restricted by the bound, which is taken into account by adjusting critical values accordingly. It is shown that explicitly accounting for possible smoothing bias in this fashion can substantially sharpen inference relative to traditional methods.

This paper proposes confidence intervals for regression kink designs that leverage this potential. Our proposed method is optimal in the sense that it minimizes the interval length amongst all honest procedures that utilize linear estimators and thus improves upon the performance of intervals based on local polynomial regression in a minimax sense. The efficiency gains in terms of the length of two-sided 95% intervals relative to the uniform MSE-optimal local linear\footnote{We focus on local linear estimators as the relevant benchmark as they are a popular choice in applied work and have leading bias proportional to the second derivative of the unknown function.} intervals for continuous designs are approximately 6 percent, and increase when the assignment variable is discrete. It follows from the arguments in armstrong2020simple that the efficiency gains relative to undersmoothing and robust bias-correction under valid bandwidth choices are even larger, which we demonstrate in a simulation study and an empirical illustration. Our results suggest that the efficiency gains relative to traditional approaches to nonparametric inference in regression kink designs are substantial. These improvements are particularly valuable for applied work that utilizes RKDs, as power concerns are common and existing procedures that take into account the smoothing bias often yield uninformative confidence intervals despite graphical evidence suggesting the existence of kink effects card2017regression. While optimally tuned honest intervals based on local polynomial regression attain on par efficiency in continuous designs, the optimization approach to inference allows us to flexibly incorporate shape constraints that can further shrink the confidence set. We extend the optimization approach to fuzzy designs by employing the inversion strategy proposed in noack2019bias and state conditions that ensure that the optimized linear confidence intervals are honest in the sense of li1989honest. Finally, we illustrate the utility of the procedure in the context of inference on the elasticity of unemployment duration with respect to the benefit level and provide software that implements the procedure.

The key tuning parameter of our proposed method is the bound on the second derivative of the CEF, which governs how much curvature one plausibly allows for. In contrast to undersmoothing and robust bias correction, our method requires that this bound is specified explicitly. However, as pointed out in armstrong2020simple, these methods can not factually avoid this choice, as they must also implicitly restrict a derivative of order two or higher in order to maintain coverage over a given nonparametric class of functions. Once this bound is specified, the proposed method is fully data driven and avoids the choices associated with local polynomial regression regarding the polynomial order, the kernel or the bandwidth. These advantages come at the cost of a closed-form expression for the estimator that we center our interval around, which is defined as the solution to a linear minimax problem over the space of distributions defined by the bound and obtained by numerical convex optimization.

The study of linear minimax estimators for nonparametric regression problems (among other problems) goes back to donoho1994statistical, who showed that, under broad conditions, the ratio of the linear minimax risk to the general minimax risk is bounded by 1.25, and that the minimax linear estimator for linear functionals over convex parameter spaces can be obtained via convex optimization. In the context of inference in discontinuity designs, armstrong2018optimal first applied this result to RDDs over a class of functions proposed by sacks1978linear that restricts the approximation error of a Taylor expansion about the discontinuity point\footnote{In this class, the minimax linear estimator can be obtained in closed form. However, since discontinuity designs are predicated on the assumption of continuity of the CEF away from the threshold it is conceptually less appealing than the Hölder class considered in imbens2019optimized, armstrong2020simple and this paper. Moreover, permitting functions that are discontinous away from the threshold also leads to inference that is too conservative at smooth functions, as the worst-case bias is larger in this class.}. In the same context, imbens2019optimized impose a bound on the second derivative and propose a numerical optimization strategy to construct honest RDD intervals based on the linear minimax estimator under homoskedasticity, emphasizing the advantages of direct optimization in settings with discrete assignment variables and multivariate assignment rules.

The method proposed in this paper contributes to the methodological literature on inference in regression kink designs (calonico2014robust; card2015inference; ganong2018permutation; noack2019bias). In implementing our method, we rely on the discretization strategy of imbens2019optimized, which is employed to the dual of the linear minimax problem for kink estimation. Conceptually, this paper builds on the "bias-aware" approach to nonparametric inference in particular (armstrong2018optimal; kolesar2018inference; imbens2019optimized; ignatiadis2019bias; noack2019bias; rambachan2019honest; schennach2020bias; armstrong2020simple) and a large methodological literature on inference in discontinuity designs in general (cf. lee2010regression; cattaneo2019practical).

The remainder of this paper is structured as follows: Section 2 introduces our notation, sketches the identification result underlying regression kink designs and defines the parameter of interest. Section 3 formally states the objective of this paper, explains how the bound permits a bias characterization that facilitates the implementation of the procedure, states assumptions under which optimized linear confidence intervals attain honesty and discusses practical questions related to their implementation. In Section 4, we conduct a simulation study to investigate the performance of optimized linear intervals relative to a variety of methods based on local linear regression. In Section 5, we demonstrate the utility of our method in a sensitivity analysis to landais2015assessing. Section 6 concludes. The Appendix contains proofs and additional results.

Framework and Notation

\paragraph{2.1 Setup.} We observe a random sample of $n$ independent pairs $(X_i,Y_i)$, where $Y_i \in \mathbb{R}$ is the outcome of interest and $X_{i} \in \mathbb{R}$ is the assignment variable for unit $i$. We write the nonparametric regression model as

so that $\mu(X_i) = \operatorname{\mathbb{E}}[Y_i | X_i]$ denotes the conditional expectation function of the outcome given the assignment variable. The subscript $N$ indicates a column vector of length $n$ where the $i$-th element corresponds to unit $i$, such that e.g. $X_{N} = (X_{1},X_{2},\cdots, X_{n})^{T}$. We normalize the threshold $c$ to zero and define an indicator function $D(x) = \mathbbm{1}_{[x \geq 0]}$. For a general function $f(x)$ we write $f_{+}(x) = f(x) D(x)$ and $f_{-}(x) = f(x) (1-D(x))$ and denote the $j$-th derivative of $f(x)$ with respect to $x$ by $f^{j}(x)$, where $f^{0}(x)$ is understood to mean the function itself. Moreover, let $f^{j}_{\pm}(0) = \lim_{ \pm x \downarrow 0} f_{\pm}^{j}(x)$ so that in particular $\mu^{1}_{+}(0) = \lim_{x \downarrow 0} \mu^{1}_{+}(x)$ and $\mu^{1}_{-}(0) = \lim_{x \uparrow 0} \mu^{1}_{-}(x)$. Let $R_{X}$ denote the support of $X_{i}$, i.e. the smallest closed set in $\mathbb{R}$ such that $\Pr \lbrace X_{i} \in R_{X} \rbrace =1$. We assume that $R_{X}$ covers zero, is bounded and denote by $\underline{x}$ and $\bar{x}$ its minimal and maximal element. Let $\mathcal{X}_{-} = [\underline{x},0)$ and $\mathcal{X}_{+} = (0,\bar{x}]$. We assume that the CEF $\mu$ is a member of the class

which formalizes the notion that $\mu$ is two times differentiable on either side of the threshold, potentially discontinuous at the threshold, and has second derivative bounded by $L$ uniformly over $\mathcal{X}= \mathcal{X}_{+} \cup \mathcal{X}_{-}$. The bound $L$ effectively governs how much curvature one allows for, as values of $L$ close to zero imply that the members of $\mathcal{F}(L)$ are close to linear, while larger values of $L$ allow for increasing amounts of curvature. Given $L$ and $\sigma_{N}$, the method that is proposed in this paper is fully data-driven and we will assume throughout the derivation that both are known. The choice of $L$ and estimation of $\sigma_{N}$ are discussed separately in Section 3.7.

Since conditional expectation functions are unique only over the support of the conditioning variable, we follow kolesar2018inference in that, assuming $\mu \in \mathcal{F}(L)$ is understood to mean that there exist $\mu$ in $\mathcal{F}(L)$ such that $\Pr\lbrace\mu(X_i) = \operatorname{\mathbb{E}}[Y_i|X_i] \rbrace=1$. While the canonical regression kink design is predicated on continuous assignment variables, our setup therefore does not assume the assignment variable to be of any specific type. This is advantageous in settings in which we only have coarse measurements at our disposal, in particular if the support of the assignment variable does not contain an open neighborhood around the threshold. In such situations the bound $L$ allows for extrapolation that ensures meaningful partial identification of $\mu^{1}_{+}(0)$ and $\mu^{1}_{-}(0)$, as discussed in Section 3.6.

\paragraph{2.2 Parameter of Interest.} Regression Kink Designs consider structural models of the form

where the function $Y: \mathbb{R}^{3} \mapsto \mathbb{R}$ describes how the outcome is produced, the random variable $E \in \mathbb{R}$ aggregates unobserved influences potentially correlated with the assignment variable, and the policy function $T: \mathbb{R} \mapsto \mathbb{R}$ is differentiable away from a kink location normalized to zero. The parameter of interest is the average partial effect of the policy variable at the kink

which inherits its causal interpretation from the definition of $Y$. Let $f_{E|X}(e|x)$ denote the conditional probability density function of $E$ given $X=x$. Under regularity conditions, we can write the first derivative of the CEF in this framework as

Since $T$ is kinked at zero, this decomposition implies that, if the average partial effect of the assignment variable and the distribution of unobservables are continuous at zero, $\tau_{RKD}$ is identified by

card2015inference characterize models for which this is the case and discuss conditions under which $\tau_{RKD}$ is equivalent to the "treatment on the treated" parameter in florens2008identification or the "local average response" parameter in altonji2005cross, respectively. In the next section, we focus on the sharp RKD, that is we assume that the denominator of $\tau_{RKD}$ is known, such that inference on $\tau_{RKD}$ is solely concerned with the jump in the first derivative of $\mu(x)$, $\theta = \mu^{1}_{+}(0) - \mu^{1}_{-}(0)$, which we refer to as the kink parameter. In Appendix (ref), we extend our method to fuzzy designs by using a strategy recently proposed by noack2019bias.

Optimized Honest Confidence Intervals

\paragraph{3.1 Problem Statement.} We seek to construct efficient confidence intervals of the form $\mathcal{I}_{\alpha} = [ \underline{\theta}, \bar{\theta} ]$ which cover the kink parameter $\theta$ with at least probability $1-\alpha$ for some prespecified level $\alpha>0$ in large samples. Furthermore, we strengthen this requirement by demanding our confidence intervals to be honest in the sense of li1989honest with respect to the class $\mathcal{F}(L)$

The uniform requirement imposed by honesty disciplines our inference in the sense that it requires us to specify and take into account plausible adversarial distributions in our asymptotic approximation. In the present setting, this means to specify $L$ and guarantee coverage for the worst-case function in $\mathcal{F}(L)$. Honesty ensures that, for any tolerance level $\eta$, we can find a sample size $n_{\eta}$ such that for $n>n_{\eta}$ coverage of $\mathcal{I}_{\alpha}$ is above $1-\alpha -\eta$ for all $f \in \mathcal{F}(L)$. As discussed in armstrong2020simple, the requirement to explicitly specify $L$ is not a disadvantage of uniform procedures, as methods that rely on pointwise guarantees of the type

must implicitly restrict $L$ to justify that coverage is controlled over a given function class $\mathcal{F}$. This holds true, irrespective of the regularization problem that such procedures need to solve.

\paragraph{3.2 General Approach.}

We center the confidence interval around a linear estimator of the form $\hat{\theta} = \sum_{i=1}^{n} w_{i}Y_i$, with weights $w_{N}$ optimized for a uniform criterion. The estimator is linear in the sense that the weights depend only on $X_N$ and a non-random tuning parameter $\kappa>0$, which we keep implicit in our notation. This choice is motivated by the relative minimax-efficiency result in donoho1994statistical and the fact that the mean squared error of linear estimators, and related quantities governing the length of our confidence intervals, depend on the unknown function only through the bias. In order to ensure honesty of the enclosing interval, we compute the magnitude of the exact worst-case conditional bias $\bar{B}(w_{N}) = \sup_{ \mu \in \mathcal{F}(L)} \operatorname{\mathbb{E}} [ \hat{\theta} - \theta | X_{N}]$ of the estimator over $\mathcal{F}(L)$ during the optimization, and inflate the interval width accordingly. Following imbens2019optimized, we obtain the estimate for the kink parameter by numerical convex optimization, solving a version of the linear minimax problem under known variance $\sigma_{N}^{2}$

Note that, since the class $\mathcal{F}(L)$ is symmetric with respect to zero, we do not need an absolute value inside the supremum, as the worst-case negative and positive biases over $\mathcal{F}(L)$ have the same magnitude. For $\kappa=1$, the objective function thus corresponds to the uniform conditional MSE of $\hat{\theta}$ over $\mathcal{F}(L)$. In general, $\kappa$ governs the worst-case conditional bias-variance tradeoff and is either chosen to minimize the uniform MSE or the interval length. Since the solution of ((ref)) depends on the data only through $X_{N}$, the confidence interval obtained by optimizing the interval length over $\kappa$ provides the same statistical guarantee as those obtained for a fixed $\kappa$.

\paragraph{3.3 Bias Characterization.} In order to translate ((ref)) into a tractable optimization problem we rely on the restrictions defining $\mathcal{F}(L)$. For any $\mu\in \mathcal{F}(L)$ and $X_{i} \in R_{X}$ we can write

where $R_{+}$ and $R_{-}$ denote the remainders of the expansions. Note that by definition $R(0)= 0$, $R^{1}(0) = 0$, and $|R^{2}(x)| \leq L$ for all $\mu \in \mathcal{F}(L)$ and $x \in \mathcal{X}$. The worst-case conditional bias $\bar{B}(w_{N})$ of a general linear estimator of $\theta$ over $\mathcal{F}(L)$ is

Since the assumption $\mu \in \mathcal{F}(L)$ does not impose any constraints on $(\mu^{0},\mu^{1})$ at zero, this expression is infinite unless the following discrete moment conditions are satisfied

As a consequence, any solution to ((ref)) must satisfy these constraints, a fact that we utilize in the implementation of our estimator. We refer to a linear estimator and the corresponding set of weights with associated $X_{N}$ as "of the correct order" if these constraints are met, in which case the conditional bias of $\hat{\theta}$ is given by the weighted sum of approximation errors of the Taylor approximation to $\mu_{\pm}$ near zero. Since $\mu \in \mathcal{F}(L)$ implies that $\mu^{1}$ is absolutely continuous, an integration by parts argument allows further characterization of the conditional bias under the constraints, yielding

Applying an argument based on Fubini's Theorem then obtains

This representation of the conditional bias characterizes the choice of an an adversarial nature that needs to pick $\mu$ out of $\mathcal{F}(L)$ in response to $(w_{N},X_{N})$ for weights that satisfy the above constraints. Let $\bar{w}(t) = D(t) \sum_{i: X_{i} \geq t }w_{i,+} (X_i - t) - (1-D(t)) \sum_{i: X_{i} < t} w_{i,-} (X_i - t).$ In this notation, the conditional bias is $\int_{\mathbb{R}} \mu^{2}(t) \bar{w}(t) dt$, which an adversarial nature maximizes by setting $\mu^{2}(t) = \operatorname{sign} (\bar{w}(t)) L$, yielding

It follows that the worst-case conditional bias of a linear estimator, if it is finite, is proportional to $L$ and that the worst-case function in $\mathcal{F}(L)$ at which it is attained is a quadratic spline with piecewise constant second derivative of magnitude $L$. Another consequence of ((ref)) is that, for any given set of weights of the correct order and bound $L$, the computation of $\bar{B}(w_{N})$ amounts to finding the roots of $\bar{w}(t)$. This has practical value, as the sign of $\bar{w}(t)$ is known for local polynomial estimators under regularity conditions, a fact we utilize in our estimation strategy.

\paragraph{3.4 Estimation via Dual Optimization.} Our approach to implementing the optimized linear confidence intervals relies on a renormalization justified by the result in ((ref)) and the convexity of ((ref)). From equation ((ref)) it follows that for linear estimators with finite worst-case bias, it holds that $\sup_{ \mu \in \mathcal{F}(L)} \operatorname{\mathbb{E}} [ \hat{\theta} - \theta | X_{N}] = L \sup_{ R \in \bar{\mathcal{F}}(1)} \sum_{i=1}^{n} w_{i} R(X_i)$ with $\bar{\mathcal{F}}(1)$ defined as

The class $\bar{\mathcal{F}}(1)$ can be understood as the class of remainder functions corresponding to the conditional expectation functions in $\mathcal{F}(1)$, which is reflected by the additional constraints that correspond to the aforementioned properties of remainders. The renormalization allows us to equivalently state the optimization problem defining $\hat{\theta}$ as follows. For a fixed $\kappa$, we write the primal problem as

The weights solving problems ((ref)) and ((ref)) are equivalent since, at the optimal value of $r$, the objective functions are identical and any candidate solution to ((ref)) lies in the feasible set for $w_{N}$ of ((ref)). This reformulation is helpful, as we require the optimal weights as well as a sharp uniform upper bound on $\operatorname{\mathbb{E}} [ \hat{\theta} - \theta | X_{N}]$ to construct our confidence interval and, more importantly, rely on the additional constraints of $\mathcal{\bar{F}}(1)$ in our implementation. The Lagrangian of ((ref)) is given by

In order to solve ((ref)) we rely on a second equivalence result. In Appendix (ref), we show that the local polynomial weights with corresponding worst-case conditional bias lie in the feasible set of ((ref)), implying that a refined Slater's condition applies to ((ref)). As a consequence, strong duality holds and any primal optimal point is also a minimizer of $L(w_{N},r,\nu^{*},\lambda^{*})$ where $(\nu^{*},\lambda^{*})$ is the solution to the dual problem

In Appendix (ref), we show that we can interchange the order of the infimum and the supremum in the dual objective by applying a minimax theorem, yielding an inner convex quadratic minimization problem that is solved analytically. This results in closed-form expressions for the primal parameters as functions of the dual parameters and a remainder function $R \in \bar{\mathcal{F}}(1),$

as well as a simplified expression for the dual objective

Let $(w^{*},r^{*}) = (w(\nu^{*},\lambda^{*}),r(\nu^{*},\lambda^{*}))$ for the element of $\bar{\mathcal{F}}(1)$ that attains the supremum in $q(\nu,\lambda)$. It follows from strong duality and strict convexity of $L(w,r,\nu^{*},\lambda^{*})$ that $(w^{*},r^{*})$ is the solution of ((ref)). As a consequence, we can recover the weights solving ((ref)) as well as the associated worst-case conditional bias over $\mathcal{F}(L)$ by the solution of the simplified dual problem

in conjunction with the mapping ((ref)) between primal and dual parameters and the result ((ref)). This translates the primal problem of $n+1$ parameters into a problem over the space $\bar{\mathcal{F}}(1)$ and five dual parameters, which we solve numerically by discretization as described in Appendix (ref).

Remark 1. Abstracting from the constraint $R \in \bar{\mathcal{F}}(1)$, the dual problem ((ref)) is a standard quadratic program and the remaining challenge is to find a suitable approximation strategy to the the functional constraint. In our implementation, we approximate the function on an equidistant grid to permit approximation of the second order constraint via finite central differences, e.g. $R^{2}(x) = [R(x+h)-2R(x) + R(x+h)]/h^2 + O(h^2)$ (see Appendix (ref) for details). This approach is formally justified by Proposition 2 of imbens2019optimized, which states that, for assignment variables with compact and convex support, the optimal weights can be recovered with arbitrary small $L_{2}$- error under the proposed discretization strategy.

Remark 2. Under the proposed strategy, the optimization approach to bias-aware inference allows us to incorporate shape constraints in a simple fashion. As explained in more detail in Appendices (ref) and (ref), any additional constraint on the CEF that can be approximated in terms of finite differences of Taylor remainders can be utilized by modifying the feasible set of ((ref)).

\paragraph{3.5 Interval Construction.}

Given a solution $(\nu^{*},\lambda^{*},R^{*})$ to ((ref)) for a fixed value of $\kappa$, we recover the optimal weights $w_{N} \gets w^{*}$ and the corresponding worst-case bias magnitude $\bar{B}(w_{N}) \gets r^{*}L$ via the mapping ((ref)) and the result ((ref)) to construct the optimized linear interval. Intuitively, the construction relies on the following decomposition of $\hat{\theta} - \theta$

By definition, $\operatorname{\mathbb{E}}[\hat{\theta} - \theta|X_{N}]$ is bounded in absolute value by $\bar{B}(w_{N})= \sup_{ \mu \in \mathcal{F}(L)}E [ \hat{\theta} - \theta | X_{N}]$ uniformly over $\mathcal{F}(L)$. Let $s_{n}^{2} = \sum_{i=1}^{n} w_{i}^{2} \sigma_{i}^{2}$ denote the conditional variance of $\hat{\theta}$ given $X_{N}$ and define the conditional bias to standard deviation ratio $t_{n} = s_{n}^{-1} \operatorname{\mathbb{E}}[\hat{\theta} - \theta|X_{N}]$. The uniform bound implies that the t-statistic

is the sum of a term that is bounded in absolute value by $\bar{t}_{n} = \bar{B}(w_{N}) s_{n}^{-1}$ uniformly over $\mathcal{F}(L)$ and a term that, under conditions stated in the next section, converges to a standard normal distribution uniformly over $\mathcal{F}(L)$ by a suitable central limit theorem. Provided that this is the case, it follows that an honest $(1-\alpha)$ confidence interval for $\theta$ is given by

where $\text{cv}_{1-\alpha}$ denotes the $(1-\alpha)$ quantile of the folded normal distribution $|N(\bar{t}_{n},1)|$ with mean $\bar{t}_{n}$ and variance one, that is the distribution of the absolute value of a normal distribution with mean $\bar{t}_{n}$ and variance one. Intuitively, this construction works because the bias can not be negative and positive at the same time, which implies that a hypothetical interval that adds and substracts $\bar{B}(w_{N}) + z_{1-\alpha/2} s_{n}$ from $\hat{\theta}$ would be too conservative. For the sharp\footnote{In fuzzy discontinuity designs, the strategy underlying the construction of the interval ((ref)) can also be employed in principle. However, in this case the smoothing bias of the first stage estimator must be dealt with and additional problems arise. In Appendix (ref), we discuss how our implementation utilizes the strategy proposed in noack2019bias to extend the optimization approach to fuzzy discontinuity designs.} regression kink design, an honest $(1-\alpha)$ interval for $\tau_{RKD}$ is then immediately obtained by rescaling the upper and lower ends of $\mathcal{I}_{\alpha}$ by the inverse of the magnitude of the kink.

We construct two types of optimized linear confidence intervals according to ((ref)): Uniform MSE-optimal ($\kappa_{UMSE}=1$) and length-optimal ($\kappa_{LE}= \operatorname*{arg\,min}_{\kappa>0} s_{n} \text{cv}_{1-\alpha}(\bar{t}_{n})$) intervals. However, the same construction principle can be applied to any uniform performance criterion that specifies a worst-case bias-variance trade-off, as the guarantees of ((ref)) hold for any fixed $\kappa>0$.\\

Remark 3. In order to obtain the length-optimal interval, we search for the optimal value $\kappa_{LE}$ using a combination of golden section search and successive parabolic interpolation as implemented in standard derivative free optimization libraries. This is feasible at high accuracy in practice as the runtime of a single optimization iteration as implemented is low (approximately $0.129$ seconds for a sample size of $6000$), leading to an average total runtime of $3.59$ seconds for the same sample size on a standard desktop computer (see Appendix (ref) for more details). This could in principle be further improved by restricting attention to values of $\kappa$ smaller than the uniform MSE-optimal choice $\kappa=1$. This is because the length-optimal weights will "oversmooth" relative to the uniform MSE-optimal weights, which armstrong2020simple show for estimators in their regularity class (cf. Figure 1 therein).

\paragraph{3.6 Theoretical Properties.}

In order to discuss the statistical properties of confidence intervals constructed according to ((ref)) we impose the following assumptions.\\

Assumption 1 Let $(C,\delta,\sigma_{\min}, \sigma_{\max}) \in \mathbb{R}_{+}^{4}$ be some fixed vector.

Assumption 1 is sufficient for a central limit theorem to apply to $s_{n}^{-1} w_{N}^{T}u_{N} =s_{n}^{-1}\left(\hat{\theta} - \operatorname{\mathbb{E}}[\hat{\theta}|X_{N}]\right)$ uniformly over $\mathcal{F}(L)$ and ensures the consistency of $\hat{\theta}$ in the identified setting. Part (i) is the standard model for survey data and provides that the kink parameter is a well defined quantity. Part (ii) and (iii) ensure that the quadratic program defining $\hat{\theta}$ is strictly convex, which guarantees that the optimal weights are uniquely recovered by the dual optimal parameters, provided that the data contains atleast two distinct points on either side of the threshold.

Assumptions (iii) and (iv) guarantee the existence of and establish bounds on the second and $(2+\delta)$-th absolute conditional moment functions of the CEF error uniformly over the support of the assignment variable and the class of permitted CEFs. The two assumptions restrict the class of permitted distributions beyond the CEF constraint (ii) in that they require uniformly bounded and non-zero conditional variances as well as the existence of a strictly finite higher order moment function. Assumptions 1 (i)-(v) are sufficient to establish that Lyapunov's condition applies to each element of $\mathcal{F}(L)$, implying convergence of $s_{n}^{-1} w_{N}^{T}u_{N}$ to a standard normal variable uniformly over $\mathcal{F}(L)$.

Assumptions (v) restricts the limit behavior of the set of optimal weights and are difficult to derive from higher-level conditions. This is because, to the best of our knowledge, a closed form solution to (1) is not known and a general characterization of $w_{N}^{*}$ beyond the spline property derived in Section 3 is difficult. While this is unattractive from a theoretical point of view, the good news is that we can verify that the condition is approximately met in any given application, and our implementation reports the finite sample counterpart to part (v). Assumption (v) together with (iii) is sufficient for $s_{n}^{2} = o_{p}(1)$ uniformly over $\mathcal{F}(L)$ and implies consistency of $\hat{\theta}$ in the identified setting.

Roughly speaking, Assumption 1 rules out distributions such that, for some $X_{i} \in R_{X}$ and $\mu \in \mathcal{F}(L)$, in the limit $w_i u_i$ is "too large", in the sense that it dominates the behavior of the sequence $s_{n}^{-1} w_{N}^{T} u_{N}$. This rules out that only a "small" proportion of the data is driving the estimate under this function. Appendix (ref) contains a formal discussion of how the relevant components of Assumption 1 can be used to show that Lyapunov's condition holds conditionally on $X_{N}$ uniformly over $\mathcal{F}(L)$, which is the key ingredient in ensuring that the optimized interval attains honesty. Once uniform convergence of $s_{n}^{-1}[\hat{\theta} - \operatorname{\mathbb{E}}[\hat{\theta}|X_{N}] \overset{D}{\rightarrow} N(0,1)$ is established, it follows from the definitions of the worst-case bias to standard deviation ratio $\bar{t}_{n}$ and the critical value $\text{cv}_{1-\alpha}$ that, for a standard normal random variable $Z\sim N(0,1)$, it holds uniformly over $\mathcal{F}(L)$ that

Taken together, the two results imply that ((ref)) is honest, which we record in Proposition 1.\\

Proposition 1 Suppose that Assumption 1 holds. Then uniformly over $\mathcal{F}(L)$

and the interval $\mathcal{I}_{\alpha} = \left[ \hat{\theta} \pm s_{n} \text{cv}_{1-\alpha}(\bar{t}_{n}) \right]$ satisfies

for $t_{n}$, $\bar{B}(w_{N})$ and $\text{cv}_{1-\alpha}$ as defined in Section 3.5.\\

The relevant difference of the statistical guarantee given in Proposition 1 relative to those that pointwise approaches to nonparametric confidence interval construction rely on is as follows: It ensures that, for any tolerance level $\eta$, one can find a sample size $n_{\eta}$ such that for all $n>n_\eta$ coverage is at least $1-\alpha-\eta$ for all functions in $\mathcal{F}(L)$. In contrast, pointwise procedures can not generally ensure the existence of such a sample size for any given non-trivial tolerance level without restricting the curvature or a higher order derivative, since the true CEF is unknown. Thus, their coverage properties can theoretically be poor even in large samples. Consequently, the uniform guarantee provided by honesty is required for reliably good finite sample performance. Once such a restriction is imposed, it follows from the definition of optimized intervals that they are the minimax optimal choice (under known variances) amongst all linear intervals.

Another attractive property of ((ref)) and "bias-aware" intervals in general is that they remain valid, irrespective of whether the assignment variable has support arbitrarily close to the threshold or not, in the sense that the statistical guarantee of the interval remains the same. In settings in which this is not the case, the interval will have positive length in the limit, but not necessarily cover the whole identification interval, that is the interval of values for $\theta$ that are consistent with the distribution $(X,Y)$ and the restriction imposed by $\mu \in \mathcal{F}(L)$. Partial identification intervals of this type were proposed in imbens2004confidence and are also useful in other non-standard situations, e.g. if one wants to exclude data for reasons such as data entry errors or other institutional characteristics that could justify such a choice.\\

In order to discuss two potential threats to the quality of the approximation underlying the construction of ((ref)), we further introduce the following two assumptions.\\

Assumption 2 Let $C_{1} \in \mathbb{R_{+}}$ be fixed and $\hat{s}_{n}$ denote an estimator at our disposal.

Assumption 2 (i) allows us to clarify the role of the ratio $\bar{w}_{R} = \left[ \max_{i}|w_{i}|\right] \left[\sum_{i=1}^{n} |w_{i}|\right]^{-1}$ addressed by Assumption 1 (v) with respect to the quality of the normal approximation underlying ((ref)). It follows from the Berry-Essen Theorem (cf. Theorem 3 in deasymptotic) that under Assumption 2 (i)

where $\Phi$ denotes the standard normal CDF and the constant $D$ lies in $0.4097 < D \leq 0.56$. This illustrates the role of Assumption 1 (v) in ensuring the quality of the distributional approximation. In particular, under the maintained assumptions, one would expect that the finite sample coverage rate of ((ref)) at the worst-case function is close to nominal whenever the ratio $\bar{w}_{R}$ is small. We therefore report $\bar{w}_{R}^{2}$ as a diagnostic statistic in our implementation and recommend to verify that this is the case in practice. If the ratio is "large", in the sense that the weights concentrate on a small set of observations, it is recommended to modify $\kappa$. In doing so, one effectively trades the quality of the estimator resulting from the initial choice of $\kappa$ in terms of the respective performance criterion for an improvement in the quality of the distributional approximation.

Assumption 2 (ii) emphasizes that the honesty property of ((ref)) was derived under known variances and that, in principle, one needs an appropriate estimator for $\sigma_{N}$ in order to preserve honesty under estimated conditional variance. Assumption 2 (ii) provides that such a uniformly consistent estimator of $s_{n}$ is available. In this case, the feasible interval that replaces $s_{n}$ with $\hat{s}_{n}$ remains asymptotically uniformly valid. This is pointed out because commonly used estimators for the conditional variance have a leading bias that is proportional to $\mu^{1}$, which is unrestricted over $\mathcal{F}(L)$. noack2019bias propose an estimator for the conditional variance based on a regression adjusted version of the nearest-neighbor estimator of abadie2014inference that has leading bias proportional to $\mu^{2}$ and can thus preserve honesty under second order bounds. Our implementation contains their proposed estimator as well as standard estimators of $s_{n}$.

Note that the estimator solving ((ref)) is the finite-sample minimax linear estimator of $\theta$ over $\mathcal{F}(L)$ only under known conditional variances $\sigma_{N}^{2}$. If the conditional variances need to be estimated, the estimator is no longer guaranteed to achieve the minimax risk in finite samples. However, in the case of homoskedasticity, it suffices to estimate $\sigma^{2}_N$ by an efficient estimator for the conditional variance to obtain honest and asymptotically minimax optimal intervals. Under heteroskedasticity, a uniformly consistent estimator for $\sigma^{2}_{N}$ is required to maintain honesty as discussed above, and the minimax properties of the estimator solving ((ref)) depend on this choice.

\paragraph{3.7 Practical Implementation.} So far we have assumed that the curvature bound $L$ and the conditional variance $\sigma_{N}^{2}$ are known. In practice, it needs to be specified how $\sigma_{N}$ should be estimated and $L$ chosen. While the previous discussion gives some guidance on how to estimate $\sigma_{N}$, the most important choice in implementing our proposed method is the choice of the tuning parameter $L$, which, without additional assumption, can not be determined from the data without undermining the honesty of ((ref)) (cf. armstrong2018optimal and references therein). This is due to a result in low1997nonparametric, who shows that when $\mathcal{F}$ is a derivative smoothness class, it is, without further assumptions, not possible to adapt to $\mathcal{F}$ while maintaining uniform coverage at the same time.

\paragraph*{3.7.1. Choice of $L$.} As a consequence of Low's impossibility result, the curvature bound $L$ has to be chosen a priori and application-specific knowledge on what constitutes plausible amounts of curvature is required to obtain suitable values of $L$. We reiterate that this requirement is not unique to bias-aware approaches to inference, as confidence intervals based on pointwise procedures must implicitly restrict $L$ to be informative at any given tolerance level.

In the absence of reliable information on the magnitude of $L$, it is recommended to conduct a sensitivity analysis by considering a range of plausible bounds together with rule of thumb (ROT) estimates of $L$ based on modelling the CEF over the largest part of its domain that is plausibly informative. In our implementation, we consider three rules of thumb. The first was suggested by armstrong2020simple and is based on fitting a global quartic polynomial on each side of the threshold. The ROT estimate of $L$ is then obtained by computing the global maximum of the absolute value of the second derivative of the polynomial implied by the estimated coefficients. The second rule of thumb we consider was proposed in imbens2019optimized. It involves fitting a quadratic polynomial on each side of the cutoff and computing an estimate of $L$ by scaling the maximum second derivative magnitude by a factor of 2-4. Finally, we propose a third rule of thumb that is based on fitting a cubic smoothing spline with evenly spaced knots on each side of the threshold. The bound $\hat{L}_{ROT}$ is then estimated by the maximum magnitude of the implied second derivative. This approach is heuristically motivated by the fact that the worst-case function of our estimator is a quadratic spline. In our simulation study, we report results based on this approach.

While it is not possible to consistently recover $L$ from the data, we are aware of two methods that were proposed to guard against overly optimistic choices of $L$ and to gain intuition for what might constitute plausible degrees of curvature. The first method is due to kolesar2018inference, who propose a method to estimate a lower bound on $L$ based on the observations that any function in $\mathcal{F}(L)$ can, between any two points that are $\Delta$ units apart, not depart from a straight line by more than $\dfrac{L \Delta^{2}}{8}$. The second method is due to noack2019bias, who propose a graphical procedure based on the solution to a constrained least squares problem to visualize "extreme" elements of $\mathcal{F}(L)$. The idea is to plot this element while iteratively increasing the curvature bound until the resulting function become implausibly erratic. We generally recommend to combine subject knowledge, ROT estimates and such heuristic devices to gain intuition in any given application.

\paragraph*{3.7.2. Estimation of $\sigma_{N}$.} The discussion in the previous section implies that the estimator derived under known variances underlying the optimized linear confidence intervals can be understood as motivated by a homoskedastic model. In order to ensure that the inference based on ((ref)) is robust to heteroskedasticity, it is required to construct confidence intervals using an approriate estimator for the conditional variances $\sigma_{N}^2$, analogous to a regression analysis that uses ordinary least squares estimators but builds confidence intervals using Eicker–Huber–White standard errors. In our implementation, we initialize $\sigma_{N}$ by a naive homoskedastic estimate to obtain the weights, before building confidence intervals based on a function that implements different heteroskedasticity-robust estimators of the conditional variance, including estimators based on standard estimates of $\sigma_{i}$ relying on the residuals of linear regressions, nearest-neighbor estimates proposed and considered in abadie2006large and abadie2014inference, as well as the uniformly consistent modification proposed in noack2019bias. The results reported in the simulation study in the next section are obtained using the nearest-neighbor approach of abadie2014inference to estimate $\sigma_{N}$ based on 10 nearest-neighbor matches.

Comparison with other Methods

In this section, we compare the performance of optimized linear confidence intervals to a variety of procedures based on local linear regression in a simulation study. In order to be precise about the comparison we briefly introduce the considered methods.

\paragraph{4.1 Local Linear Methods.} The local linear estimate $\hat{\theta}_{LL}$ of $\theta$ with bandwidth $h$ is the coefficient on $X_i D_i$ in a weighted OLS regression of $Y_{i}$ on the vector $(1,X_i,D_i,X_iD_i)$, using only observations $i$ such that $|X_{i}| \leq h$, with weights determined by a kernel function. Under regularity conditions, it holds that if the density of the assignment variable $f_{X}(x)$ is continuous and bounded away from zero in an open neighborhood around the threshold, the MSE of the local linear estimator with bandwidth sequence $h_{n} \rightarrow 0$ evaluated at a function $\mu$ is

where $B \propto (\mu^{2}_{+}(0)-\mu^{2}_{-}(0))$ is the leading asymptotic bias and $V \propto (\sigma_{+}(0)^{2} + \sigma_{-}(0)^{2})f_{X}(0)^{-1}$ is the asymptotic variance. If $B \neq 0$, an asymptotic MSE-optimal bandwidth sequence is thus

In our simulation study, we compare the performance of uniform procedures to methods that rely on plug-in estimates of this quantity which add a regularization term to the denominator that shrinks with the sample size imbens2012optimal. We denote such estimates by $\hat{h}_{PMSE}$ and compute them using the plug-in estimators proposed by calonico2014robust. The undersmoothing bandwidths are computed relative to the obtained pointwise asymptotic MSE-optimal estimate $\hat{h}_{US} = n^{-1/20} \hat{h}_{PMSE}$. The two bandwidth choices required for the RBC intervals are either both set pointwise optimal $b= \hat{b}_{PMSE}$, $h=\hat{h}_{PMSE}$, where $\hat{b}_{PMSE}$ refers to the pointwise asymptotic MSE-optimal plug-in estimate for the local quadratic bias estimator, or both set to the pointwise asymptotic MSE-optimal estimate $\hat{h}_{PMSE}$ for the local linear estimator. In addition, we consider RBC bandwidth choices $h_{CE}$, $b_{CE}$ that optimize the pointwise asymptotic coverage error calonico2018coverage, which can be considered an intermediate form of undersmoothing and robust bias correction.

Under regularity conditions and restrictions on the rate of the bandwidth sequence $h_{n} \rightarrow 0$

where $B_{n} \overset{P}{\rightarrow} B$. The methods that we consider differ in whether and how they take into account the incorrect centering induced by the smoothing bias. In our simulation, we consider the following four types of two-sided local linear intervals for the above bandwidth choices:

where $z_{1-\alpha /2}$ denotes the $(1-\alpha /2)$ quantile of the standard normal distribution, $\hat{V}$ denotes an estimate of the respective asymptotic variance, and $\text{cv}_{1-\alpha}$ as well as $t_{n}$ are defined as in ((ref)) for the worst-case magnitude of the conditional bias of $\hat{\theta}_{LL}(\hat{h}_{FL})$.

The first two types, conventional and undersmoothed confidence intervals, essentialy assume the smoothing bias away. While the conventional method directly assumes that $h_{n} B_{n} s_{n}^{-1} \approx 0$, undersmoothed intervals rely on the asymptotic promise that $h_{US}/h_{PMSE} \underset{n \rightarrow \infty}{\rightarrow} 0$, implying that $\sqrt{nh_{n}^{3}} \left[ \hat{\theta}(h_{n}) - \theta - h_{n} B_{n} \right] = \sqrt{nh_{n}^{3}} \left[ \hat{\theta}(h_{n}) - \theta \right] + o_{p}(1) \overset{d}{\rightarrow} N(0,V)$, which is uninformative about the smoothing bias in a given sample.

The last two types explicitly address the smoothing bias. RBC intervals are centered around a bias-corrected point estimate, using a higher order local polynomial estimator to estimate the bias. They differ from traditional bias-corrected intervals, which are known to perform poorly in finite samples hall1992effect, in that they do not require $h_{RBC}/b_{RBC} \underset{n \rightarrow \infty}{\rightarrow} 0$. As a consequence, the standardized bias-correction term is not negligible asymptotically, leading to a different asymptotic variance $V_{RBC}$ that captures the additional uncertainty introduced by bias estimation. In our simulation, all RBC bias estimates are based on local quadratic estimators. The fixed length intervals are constructed in the same spirit as the optimized intervals discussed in Section 3 with weights defined by the local linear estimator. They take into account the exact magnitude of the worst-case conditional bias by inflating the critical value according to the ratio $\bar{t}_{n}$, and are therefore valid for any bandwidth choice. The bandwidths for the fixed length local linear intervals are obtained analogously to $\kappa_{UMSE}$ and $\kappa_{LE}$ by minimizing the finite sample uniform MSE $h_{UMSE}$ or the interval half-length $h_{HL}$. The variances for interval construction are estimated by the conditional variances of the estimators implied by the weights, using the nearest-neighbor approach of abadie2014inference based on 10 nearest neighbor matches.

\paragraph{4.2 Monte Carlo Setup.}

In order to investigate the performance of the optimized linear confidence intervals relative to the methods introduced above we conduct a simulation study. The setups differ in the conditional mean function, its degree of curvature and the distribution of the assignment variable, yielding a total of 8 different settings. The assignment variable is drawn from an equidistant uniform distribution with support $\lbrace -1, -1 + \frac{2}{K}, \cdots, 1- \frac{2}{K}, 1 \rbrace$, where the parameter $K$ controls the number of support points and $K_\infty$ means the continuous uniform distribution with support $[-1,1]$. The outcome data is generated according to

where the CEF error $\varepsilon_{i}$ is drawn from a mean zero normal distribution with $\sigma=0.1$ and

where $s_{+}^{2}(x) = D(x)x^{2}$ denotes square of the plus function. Note that both CEFs are second order splines with maximal second order derivative magnitude $L$ and thus elements of $\mathcal{F}(L)$. The function $\mu_{1}$ attains the second order bound only in $(-0.15,0.15)$ while $\mu_{2}$ attains the bound everywhere on $[-1,1]$ with alternating signs in each interval defined by the knots. In the simulation, we set $\theta=-0.5$ and consider bounds $L \in \lbrace{2,6 \rbrace}$. Figure (ref) displays the shape of the functions.\\

\paragraph{4.3 Monte Carlo Results.}

Tables (ref) and (ref) show the results of 5000 Monte Carlo runs\footnote{In Appendix (ref), we provide analogous results for 20.000 Monte Carlo runs that were conducted on a cluster due to the required computational resources.} with sample size $n=2000$ for $\mu_{1}$ and $\mu_{2}$ respectively. The top panel in each table displays the results for the case when the assignment variable is drawn from the continunous uniform distribution, while the bottom panel shows the results for 80 equidistant support points. We use a triangular kernel for all local linear methods. In the case that a bandwidth selector chooses a bandwidth for which the respective estimator is not defined, we manually adjust the bandwidth such that it covers three support points on either side of the cutoff.\footnote{This occured only in the discrete setting.} The left panel in each table reports the results for the low curvature version of the respective CEF, while the right panel reports the results for the high curvature variant. Columns 1 and 2 indicate the method and the tuning target. The curvature bound $L$ is either chosen by the rule of thumb $\hat{L}$ or fixed to 2 or 6, as indicated in the tuning subscript. We report the empirical coverage rate at nominal level 95% (Cov.), the average length relative to the optimal linear interval with correct curvature bound (RL), as well as the average tuning parameter choice ($h/\kappa$) of each method.

Unsurprisingly, conventional and undersmoothing confidence intervals show below nominal coverage in all designs, with undercoverage of undersmoothed intervals becoming more severe in the high curvature regime. The performance of robust bias-corrected intervals varies with the tuning target. While both, the default RBC method\footnote{This refers to the default in the authors' R implementation of RBC CIs.} that picks both bandwidths using the respective asymptotic MSE-optimal estimate and the coverage error optimized RBC interval undercover severly, the RBC interval obtained by setting both bandwidths to the local linear pointwise MSE bandwidth or the UMSE bandwidth under the ROT estimate of $L$ show close to nominal coverage and are insensitive to the true curvature. However, for the latter this comes at a cost in terms of their length relative to fixed-length and optimized intervals. In all designs, RBC intervals with tuning that attains close to nominal coverage are at least approximately twice as long as the infeasible length-optimal interval and 40% longer than feasible length optimized intervals under the ROT choice of $L$. Fixed-length and optimized intervals show above or close to nominal coverage under the correct curvature, with the lowest empirical coverage at 94.7% $(\mu_{1,2},K=80,L=6)$. As one would expect, they are conservative when the true curvature is lower than specified and undercover when the true curvature is higher than specified. The rule of thumb choice of $L$ tends to overestimate the curvature with the exception of $(\mu_{1},K_{\infty},L=6)$, leading to above nominal coverage of both uniform interval types in most

\setcounter{table}{0}

\newgeometry{top=0.75cm,bottom=1.25cm,footskip=0cm}

\restoregeometry

\newgeometry{top=0.75cm,bottom=1.25cm,footskip=0cm}

\restoregeometry

designs with lowest empirical coverage at 93.0% $(\mu_{1},K_{\infty},L=6)$. The resulting intervals thus tend to be conservative but are still substantially shorter than their pointwise counterparts under valid tuning. Interestingly, the length-optimal bias-aware intervals show higher coverage rates than their uniform MSE-optimal counterparts on average. Moreover, the length-optimal weights of both, fixed length and optimized estimators "oversmooth" relative to the uniform MSE-optimal choice of the respective tuning parameter. In all designs, the optimally tuned fixed-length intervals demonstrate performance on par with their optimized counterparts, indicating that the high minimax efficiency of local linear estimators under second order bounds demonstrated in armstrong2020simple for estimating the value of the CEF at a point also holds for the estimation of first derivatives.

\paragraph{4.4 Gains from Optimization.} The popularity of local linear estimators in empirical practice is motivated by their intuitive appeal and a range of attractive theoretical properties of local polynomials, in particular their asymptotic minimax efficiency over the Taylor class of functions (fan1993local, fan1997local, cheng1997automatic). armstrong2020simple show that local polynomial estimators can also attain high minimax efficiency in the Hölder class of functions defined by derivative bounds among a large class of estimators to which a central limit theorem applies and that have worst-case bias and standard deviation that scale as powers of a bandwidth parameter. This result relies on their observation that, for relevant performance criteria, the asymptotic minimax performance of two estimators in this class does not depend on the criterion but is solely governed by their worst-case biases, their standard deviations and their rate exponents $r = \gamma_{b}/(\gamma_{b}-\gamma_{s})$, where $\gamma_{b}$ and $\gamma_{s}$ denote the scaling exponents of the worst-case bias and the standard deviation respectively. Moreover, they show that the optimal worst-case bias to standard deviation ratio depends only on the criterion and $r$ (cf. Theorem 2.1 in armstrong2020simple).

Their analytic results provide us with guidance on what to expect with respect to the asymptotic efficiency gains of optimized linear confidence intervals and allow us to state a lower bound for the efficiency gain of our method relative to uniform MSE-optimal fixed length intervals. For the local linear estimator of the kink parameter, $r=0.4$ and their calculations imply an efficiency gain of approximately 6% at $\alpha= 0.05$ for moving from the uniform MSE-optimal to the length-optimal local linear interval (cf. Figure 3 in armstrong2020simple), which is consistent with our Monte-Carlo results. Thus, a lower bound for the efficiency gain of the optimized interval relative to the UMSE-optimal fixed length interval is 6%. This

implies a larger lower bound for the efficiency gain relative to undersmoothed and RBC intervals based on valid bandwidth choices. Figure (ref) shows the relative risk of fixed length length-optimal and uniform MSE-optimal estimators in terms of interval length (a) and uniform MSE (b) relative to their optimized linear counterparts for the respective criterion for the continuous and discrete ($K=40$) uniform design for different sample sizes. It illustrates that the efficiency gains from optimization increase as the discreteness of the assignment variable becomes more severe.\footnote{This point was previously made in imbens2019optimized. The convergence of the risk observed in Figure (ref) is due to manual bandwidth adjustments to ensure that the estimators are well defined.} This is because the shape of the optimized weighting function increasingly deviates from the local linear weights as the coarseness of the data becomes more severe, as shown in Figure (ref), which depicts uniform MSE-optimal local linear and optimized weights for different degrees of discreteness. In continuous designs fixed-length and the optimal linear kernels are nearly identical and the only advantage of the optimization based approach is that it can be easily modified to sharpen inference via shape constraints as discussed in Appendix (ref).

Empirical Illustration

We apply our method to the data of landais2015assessing, who estimates the effect of unemployment benefits on the duration of unemployment in a regression kink design. The paper exploits kinks in the schedule of unemployment benefits arising from a hard cap at a maximum benefit amount $b_{max}$. In the US, the weekly benefit amount $b$ received by an eligible unemployed is a fixed fraction $\gamma$ of a function of previous quarterly earnings $hqw$ in a base period up to the cap.

landais2015assessing reports estimates of $\tau_{RKD}$ for five US states: Louisiana, Idaho, Missouri, New Mexico and Washington. For the sake of exposition, we focus on the results for Louisiana, which serves as the leading example in the paper. Figure (ref) displays the benefit schedule for Louisiana for the time period covered by the data. Due to adjustments, the maximum benefit level changed over time, resulting in five distinct kinks. For the period under consideration the weekly benefit rate was fixed at $\gamma=0.04$, which corresponds to a constant replacement rate of \\

52% up to the respective kink, from where onwards the replacement rate decreases.

The paper utilizes data from the Continuous Wage and Benefit History (CWBH), a publicly available administrative UI data set for the US that contains the universe of unemployment spells and wage records for the five US states from the late 1970s to 1984, with different states starting the recording at different points in time. For Louisiana, the dataset contains $n=44702$ unemployment spells for the whole time period. See landais2015assessing Section II.A for a detailed description of the data. Figure (ref) plots the pooled data for all five time periods, where we have normalized the assignment variable (highest quarterly earnings) by the location of the respective kink and included only data within an 85% interval of the kink $[0.15,1.85]$. This corresponds to bandwidths in the range of 3000-4300 USD. In order to reduce noise, it shows the average unemployment durations, measured in weeks, in 30 equally wide bins. The paper reports RKD estimates of the effect of the benefit level on unemployment duration separately for each time period, with point estimates and standard errors rescaled to the 2010 USD price level. In the main specification, landais2015assessing estimates $\tau_{RKD}$ using local linear regression with a fixed bandwidth $h=2500$\footnote{In the paper's online Appendix, the author reports robustness checks for the pooled sample consisting of the last two periods using bandwidths $1500$ and $4500$.} and reports conventional 95% confidence intervals for $\tau_{RKD}$ based on Eicker-Huber-White standard errors. Table (ref) replicates\footnote{The point estimates in the top-left panel for period 3 and 4 deviate by $0.001$ from the results reported in the paper. We attribute this difference to rounding.} the results reported in Table 2 of the paper and additionally presents robust bias-corrected ($h=b=\hat{h}_{PMSE}$) and optimized linear $(L=\hat{L}_{ROT})$ point estimates and 95% confidence intervals computed on the same data.

The point estimates reported in Table (ref) correspond to the estimated effects of a 1 USD increase in weekly benefits on the duration of paid unemployment in weeks. For example, the point estimate for period 4 reported in landais2015assessing (top-left panel) suggest that a 1 dollar increase in weekly unemployment benefits leads to a $0.043$ weeks increase in the duration of unemployment at the kink. The estimates in the top-left panel correspond to (local) elasticities in the range of .2 and .7, suggesting that a 10% increase in the average weekly benefit amount\\

increases unemployment duration by 2 to 7% on average at the kink. As can be seen from the confidence intervals reported in the remaining panels, this finding is sensitive to potential smoothing biases for most of the time periods under consideration, with both RBC and optimized intervals covering zero at the 95% level for periods 1-4. However, while the RBC point estimates differ substantially from those reported in the top-left panel, with confidence intervals that are mostly uninformative for the question at hand, our procedure yields point estimates relatively close to those reported in landais2015assessing and lower bounds of 95% confidence intervals that are marginally below zero for most periods. This is because the estimated curvature is rather low relative to the estimators standard deviation in all data sets except for period 3, with $\hat{L}_{ROT} \times 10^5 = (0.0137, 0.058, 1.99, 0.477, 0.145)$. For the same reason, the difference between length-optimal and UMSE-optimal intervals is rather small in this data set as shown in the lower panel of Table (ref). Overall, these results indicate that the uncertainty associated with the estimated effect of unemployment benefits on unemployment duration is higher than suggested by the top-left panel unless one is certain about the linearity of the CEF in the domain specified by the bandwidth.

As mentioned earlier, the optimization approach to bias-aware RKD inference allows us to readily sharpen inference by imposing shape constraints on the CEF. This is often useful, as the plateaus in the schedules that are typically exploited in regression kink designs often give rise to empirically plausible concavity and convexity restrictions. In Appendices (ref) and (ref), we discuss this further and explain how shape constraints are implemented in practise. Table (ref) presents the counterparts of the optimized intervals reported in the lower panel of Table (ref) under the restriction that $\mu$ is a concave function, demonstrating that such constraints can substantially shrink our confidence set.\\

Conclusion

Motivated by the finite sample coverage problems of pointwise approaches to nonparametric inference, this paper proposes a robust and efficient alternative method for the construction of nonparametric confidence intervals in regression kink designs. Given a curvature bound, the method is fully data driven, easy to implement, and has excellent finite sample coverage and length properties due to its minimax construction that explicitly takes into account the worst-case smoothing bias in a given data set.\\