Local Polynomial Estimation at a Regression Discontinuity Cutoff

Why Estimate Two Limits?

In a sharp regression discontinuity (RD) design, treatment changes deterministically when a running variable crosses a known cutoff. The causal comparison is not between all treated and untreated units. It is the difference between the outcome limits for units just above and just below the threshold. Since the observations do not generally lie exactly at the cutoff, estimating those limits requires fitting the outcome relationship locally on each side and evaluating each fit at the boundary.

A global straight-line regression can obscure the feature of interest: outcomes may vary nonlinearly with the running variable, and the relationship need not have the same slope on both sides. Local polynomial estimation instead uses observations near the cutoff, weights them by proximity, and allows separate shapes on the two sides. The resulting intercepts estimate the two limits; their difference estimates the sharp RD treatment effect under the continuity-based identification conditions.

Definition: Let $X_i$ be unit $i$'s running variable, $c$ the cutoff, and $D_i=1\{X_i\ge c\}$ the sharp treatment assignment. Write $Y_i(1)$ and $Y_i(0)$ for the potential outcomes under treatment and control. Define $m_d(x)=E[Y_i(d)\mid X_i=x]$, for $d\in\{0,1\}$. If these conditional mean functions are continuous at $c$, the treatment effect at the cutoff is $$ \tau(c)=m_1(c)-m_0(c) =\lim_{x\downarrow c}E[Y_i\mid X_i=x]-\lim_{x\uparrow c}E[Y_i\mid X_i=x]. $$

Here the right-hand limit is from above and the left-hand limit is from below. The equality uses the sharp assignment rule: just above the cutoff, observed outcomes are treated potential outcomes; just below, they are control potential outcomes. Continuity of the potential-outcome conditional means lets the corresponding limits identify their values at $c$. The target is a local causal effect at the threshold, not an effect for every value of $X_i$.

Separate Kernel-Weighted Polynomial Fits

Choose a bandwidth $h>0$, a polynomial order $p$, and a nonnegative kernel function $K$, typically zero outside a bounded interval. The bandwidth restricts attention to observations near $c$; the kernel gives relatively more weight to some observations within that neighborhood. For example, the triangular kernel is $K(u)=1-|u|$ when $|u|\le 1$, and zero otherwise. The scaled distance $u_i=(X_i-c)/h$ expresses proximity in bandwidth units.

For each side $s\in\{-,+\}$, retain observations with $X_i<c$ for $s=-$ and $X_i\ge c$ for $s=+$. Let $r_p(u)=(1,u,\ldots,u^p)^\mathsf{T}$. Estimate the polynomial coefficient vector by minimizing the weighted sum of squared residuals on that side:

\widehat{\beta}_{s,p} =\arg\min_{\beta\in\mathbb{R}^{p+1}} \sum_{i:\,X_i\text{ on side }s} K\!\left(\frac{X_i-c}{h}\right) \left[Y_i-r_p\!\left(\frac{X_i-c}{h}\right)^\mathsf{T}\beta\right]^2.

The fitted value at the cutoff is the intercept, because $r_p(0)=(1,0,\ldots,0)^\mathsf{T}$. Thus the estimated jump is

\widehat{\tau}_p(c)=\widehat{\beta}_{+,p,0}-\widehat{\beta}_{-,p,0}.

The minimization has a unique solution when the weighted design matrix on that side has full column rank. In practical terms, there must be enough positively weighted observations at distinct running-variable values to identify the chosen polynomial. A higher order requires more local information; it does not create information that the data lack. The procedure also requires a neighborhood with observations on both sides and a kernel that assigns positive weight to usable observations there.

Using the scaled distance $u_i$ makes the intercept directly interpretable as the cutoff fit and helps express the local neighborhood consistently. One could use powers of the unscaled distance $X_i-c$ instead; rescaling changes slope and higher-order coefficients, but not the fitted value at $c$.

Local Linear Estimation and Its Derivation

For $p=1$, each side gets its own local line. On one side, write $w_i=K((X_i-c)/h)$ and $x_i=X_i-c$. Minimizing $\sum_i w_i[Y_i-a-bx_i]^2$ gives the weighted normal equations. Define $S_k=\sum_i w_i x_i^k$, $T_0=\sum_i w_iY_i$, and $T_1=\sum_i w_ix_iY_i$, with each sum restricted to that side. The equations are

\begin{pmatrix}S_0&S_1\\S_1&S_2\end{pmatrix} \begin{pmatrix}\widehat a\\\widehat b\end{pmatrix} = \begin{pmatrix}T_0\\T_1\end{pmatrix}.

Solving the two equations gives the intercept

\widehat a =\frac{S_2T_0-S_1T_1}{S_0S_2-S_1^2}, \qquad \widehat b =\frac{S_0T_1-S_1T_0}{S_0S_2-S_1^2}.

The denominator is positive when the positively weighted running-variable values vary on that side. The intercept $\widehat a$ is the local linear estimate of that side's outcome limit. Repeating the calculation separately on the other side and subtracting the intercepts gives the estimated jump.

This derivation also shows what the local linear intercept does: it adjusts the weighted outcomes for a fitted local slope rather than treating their weighted average as the boundary value. Because the observations are one-sided, a line fitted near the cutoff is extrapolated to $x_i=0$. That extrapolation is deliberate: the estimand is the limit at the boundary, not an average outcome across the entire neighborhood.

Worked Example: A Triangular-Kernel Local Linear Jump

Suppose hypothetical sharp RD data use cutoff $c=0$, bandwidth $h=2$, and the triangular kernel. The running-variable distances and outcomes are $(-1.5,6),(-1,8),(-0.5,9)$ below the cutoff and $(0.5,12),(1,14),(1.5,17)$ at or above it. Outcomes are measured in points. All six distances are within the kernel's support, and each side has three distinct running-variable values, enough to fit a line.

The kernel weights at distances $0.5,1,1.5$ are $0.75,0.5,0.25$, respectively. On the right, the weighted sums are $S_0=1.5$, $S_1=1.25$, $S_2=1.25$, $T_0=20.25$, and $T_1=17.875$. Therefore,

\widehat a_+ =\frac{(1.25)(20.25)-(1.25)(17.875)} {(1.5)(1.25)-(1.25)^2} =\frac{2.96875}{0.3125} =9.5\text{ points}.

On the left, using the signed distances $-1.5,-1,-0.5$, the corresponding sums are $S_0=1.5$, $S_1=-1.25$, $S_2=1.25$, $T_0=12.25$, and $T_1=-9.625$. Thus,

\widehat a_- =\frac{(1.25)(12.25)-(-1.25)(-9.625)} {(1.5)(1.25)-(-1.25)^2} =\frac{3.28125}{0.3125} =10.5\text{ points}.

The estimated jump is $9.5-10.5=-1.0$ point. Under the sharp RD continuity conditions, this is an estimate of the treatment effect at the cutoff. It is not evidence that treatment lowers outcomes for units at all running-variable values; the conclusion concerns the local effect at $c=0$.

Why the Boundary Creates Bias

At an interior point, a symmetric neighborhood can balance observations on both sides of that point. At the RD cutoff, each fit uses observations from only one side. Even a symmetric kernel does not restore two-sided balance once the sample is restricted to $X_i<c$ or $X_i\ge c$. A local constant fit, $p=0$, is just a kernel-weighted average. If the conditional mean has a slope, that average generally differs from the value at the boundary: all the right-side observations are to the right, or all the left-side observations are to the left.

A local linear fit estimates both a level and a slope, allowing it to extrapolate back to the cutoff. Under sufficient smoothness, its leading error at the boundary is typically of order $h^2$, rather than the order-$h$ error of a local constant fit. More generally, a local polynomial of degree $p$ reproduces polynomials through degree $p$; with sufficient smoothness, its leading boundary bias is typically of order $h^{p+1}$. These are bias orders for the fitted limit under regularity conditions, not guarantees that every finite-sample estimate has small error. Particular designs can also make a leading bias term vanish.

The logic can be seen directly by considering a quadratic conditional mean, $m(t)=a+bt+qt^2$, where $t=X_i-c$. A local constant fit averages the slope and curvature contributions instead of removing them. A local linear fit absorbs the constant and linear components, but curvature can still affect its intercept. The one-sided weights determine the size and direction of that remaining curvature bias.

Worked Example: Comparing Boundary Bias on One Side

Use the hypothetical right-side distances and triangular weights from the previous example: $t=0.5,1,1.5$, with weights $0.75,0.5,0.25$. Suppose the conditional mean is $m(t)=10+2t+0.2t^2$, measured in outcome points. The true limit is $m(0)=10$ points. The total weight is $1.5$, the weighted mean distance is $1.25/1.5=0.8333$, and the weighted mean squared distance is also $1.25/1.5=0.8333$. The local constant fitted value is therefore

\widehat m_{\mathrm{constant}} =10+2(0.8333)+0.2(0.8333) =11.8333\text{ points}.

Its difference from the true limit is $11.8333-10=1.8333$ points for this weighted design. For the local linear fit, the linear term contributes exactly its own slope to the fitted line. To find the intercept contribution from $t^2$, the weighted sums are $S_0=1.5$, $S_1=1.25$, $S_2=1.25$, and $S_3=1.4375$. The intercept from regressing $t^2$ on a constant and $t$ is

\frac{S_2S_2-S_1S_3}{S_0S_2-S_1^2} =\frac{1.5625-1.796875}{0.3125} =-0.75.

Consequently the local linear intercept is $10+0.2(-0.75)=9.85$ points. It is closer to the true limit in this example, though it still has curvature bias. This calculation illustrates why fitting a slope addresses first-order boundary bias without implying that a local linear fit removes all bias.

Higher-Order Fits and Practical Interpretation

A quadratic fit ($p=2$) adds a local curvature term. On each side, the fitted polynomial is evaluated at $u=0$, so its intercept remains the estimated outcome limit. If the conditional mean is truly quadratic over the fitted neighborhood and the weighted design has full rank, the quadratic fit reproduces that function exactly. In applications, however, the true function is not known, and a higher-order polynomial can be sensitive to sparse or unevenly distributed observations.

Worked Example: A Quadratic Fit Recovers the Cutoff Jump

Consider hypothetical observations with cutoff $c=0$, bandwidth $h=2$, and the same triangular kernel. Suppose the right-side outcomes at distances $t=0.5,1,1.5$ are $11.75,15,19.75$. These values lie exactly on $m_+(t)=10+2t+3t^2$. Suppose the left-side outcomes at $t=-0.5,-1,-1.5$ are $8,9,11$, lying exactly on $m_-(t)=8+t+2t^2$. All weights are positive, and each side has three distinct observations, so its quadratic design matrix is full rank.

A quadratic fit on the right recovers intercept $10$, since $m_+(0)=10+2(0)+3(0)^2=10$. The separate quadratic fit on the left recovers intercept $8$, since $m_-(0)=8+0+2(0)^2=8$. Subtracting gives $\widehat{\tau}_2(0)=10-8=2$ outcome points. The exact recovery follows here because each side's data lie precisely on its fitted quadratic; it is not a general claim that a quadratic model removes all finite-sample error.

A polynomial order is a modeling choice about local shape. Local linear estimation is often useful because it allows a slope and has favorable boundary-bias behavior without adding curvature parameters. A higher order can represent more curvature, but requires a better-supported local design. Kernel choice changes the relative influence of observations; it does not change the estimand, which remains the difference between the two cutoff limits. Bandwidth choice also matters because it determines the neighborhood used to approximate those limits. The appropriate choice and its consequences require explicit consideration of the bias-variance tradeoff.

Key takeaway: Fit the outcome relationship separately on each side of the cutoff, evaluate each fitted polynomial at the boundary, and subtract the left intercept from the right. The causal interpretation belongs to this cutoff jump and depends on the sharp assignment rule and continuity of the potential-outcome conditional means.

Common Mistakes

Using one global line without separate side behavior. A single slope can force the two sides to share a shape and may fail to represent the jump. Fit each side separately so the two limits are allowed to differ.
Calling a one-sided weighted average the cutoff limit. The local constant average generally includes slope-related boundary bias. A local linear fit estimates a slope and evaluates its fitted line at the cutoff.
Interpreting the jump as an effect everywhere. The RD estimand is local to the threshold. Do not extend it to units far from the cutoff without additional assumptions.
Assuming higher polynomial order automatically improves estimation. A higher-order fit needs enough positively weighted, distinct running-variable values. It may be unstable when the local design is sparse.
Ignoring the kernel's support or the side restriction. Only observations with nonzero kernel weight contribute, and each fit must use the correct side of the threshold.

Check Your Understanding

Answer each question using the local polynomial framework.

Why is the sharp RD effect at the cutoff expressed as the right-hand conditional outcome limit minus the left-hand limit?
For a local linear fit, what does the intercept represent, and why is the slope included?
What condition on the weighted running-variable values is needed for a unique local linear fit on each side?
Why can a one-sided local constant fit be biased for the outcome limit even when the kernel is symmetric?
If a quadratic fit has an intercept of 14 points on the right and 11 points on the left, what cutoff jump does it estimate, and what is the scope of its causal interpretation?