Theoretical Background

This article gives the mathematical background for the bias-bound approach, based on Schennach (2020).

The bias-variance tradeoff

In nonparametric estimation we face a fundamental tradeoff:

Large bandwidth: low variance but high bias
Small bandwidth: low bias but high variance

Traditional approaches take one of two routes:

Undersmooth: use a smaller bandwidth to reduce bias, which inflates variance
Ignore bias: use the optimal MSE bandwidth but produce invalid confidence intervals

The bias-bound approach takes a third route: it bounds the bias. This permits inference at an MSE-oriented bandwidth, but the explicit bound can make the resulting interval conservative and wider in finite samples.

A Fourier view of the bias

For a sample \(X_1, \ldots, X_n\) from a density \(f\), the kernel density estimator is

\[\hat{f}_h(x) = \frac{1}{nh} \sum_{i=1}^{n} K\left(\frac{x - X_i}{h}\right)\]

where \(K\) is the kernel and \(h\) is the bandwidth. Its error splits into two parts:

\[\hat{f}_h(x) - f(x) = \underbrace{[\hat{f}_h(x) - E[\hat{f}_h(x)]]}_{\text{variance term}} + \underbrace{[E[\hat{f}_h(x)] - f(x)]}_{\text{bias term}}\]

The variance term is random with a known distribution. The bias term is deterministic but unknown, and the method must control it.

The bias-bound approach works with the Fourier representation of the bias. For a kernel estimator,

\[E[\hat{f}_h(x)] - f(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} [K^{FT}(h\xi) - 1] f^{FT}(\xi) e^{i\xi x} d\xi\]

where \(K^{FT}\) and \(f^{FT}\) are Fourier transforms. For the bounded-variation derivative class used by Schennach (2020), the Fourier transform satisfies the polynomial envelope

\[|f^{FT}(\xi)| \leq A |\xi|^{-r}\]

where \(A\) is an amplitude constant and \(r\) measures smoothness (a larger \(r\) means a smoother function). The package detects \((A, r)\) from the data by fitting the empirical Fourier transform:

Schennach’s formal class uses an integer \(r\in\{2,3,\ldots\}\), bounds the total variation of the \((r-1)\)st derivative by \(A\), and imposes integrability and boundary conditions. This is related in spirit to Sobolev smoothness because both constrain derivatives, but it is not an RKHS assumption and is not equivalent to a standard Sobolev ball. The bounded-variation condition is what gives the pointwise power-law Fourier envelope above.

# Generate two-fold-uniform data
set.seed(42)
x_data <- runif(500) + runif(500)

# Estimate density
fit <- bias_bound_density(x_data, bw = 0.08, kernel = "schennach")

# View detected smoothness parameters
coef(fit)
#>        A        r        h 
#> 4.262401 2.000000 0.080000

# Visualize Fourier transform fit
plot(fit, type = "ft")

The legend labels each line:

Empirical |phi|: the empirical Fourier transform magnitude
Fitted envelope: the fitted envelope \(A|\xi|^{-r}\), drawn over the selected window
Shaded band: the frequency range used for fitting

Constructing the bias bound

Given the smoothness envelope, the largest possible bias at \(x\) is

\[\bar{b}(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} |K^{FT}(h\xi) - 1| \min\{1,A|\xi|^{-r}\}\,d\xi,\]

which the package evaluates numerically. This \(\bar{b}\) is the worst-case bias consistent with the fitted smoothness envelope, and under the paper’s function class the true bias satisfies

\[|E[\hat{f}_h(x)] - f(x)| \leq \bar{b}(x).\]

A textbook confidence interval ignores the bias,

\[CI_{\text{naive}} = \hat{f}(x) \pm z_{1-\alpha/2} \hat{\sigma}(x),\]

so its coverage is wrong when the bias is non-negligible. The bias-bound interval adds the bound on both sides,

\[CI_{\text{bias-bound}} = [\hat{f}(x) - \bar{b}(x) - z_{1-\alpha/2}\hat{\sigma}(x), \quad \hat{f}(x) + \bar{b}(x) + z_{1-\alpha/2}\hat{\sigma}(x)],\]

accounting for the worst-case bias in both directions:

# The plot shows both bands
plot(fit)

The legend labels the two bands:

Bias bound: the bias range \([\hat{f} - \bar{b}, \hat{f} + \bar{b}]\)
95% CI: the full confidence interval including sampling uncertainty

Kernels

For the formal coverage result, an infinite-order kernel whose Fourier transform equals one on a neighborhood of the origin is required. Write the radius of that neighborhood as \(c_K>0\):

\[K^{FT}(\xi) = 1 \text{ for } |\xi| \leq c_K.\]

This removes the kernel contribution to the bias integrand for frequencies below \(c_K/h\). Four kernels are available:

Kernel	Order	Fourier Transform
Schennach2004	\(\infty\)	One through \(\|\xi\|=0.5\); smooth transition to zero at \(1.5\)
sinc	\(\infty\)	Sharp cutoff at \(\|\xi\|=\pi\) under the package convention
normal	2	Gaussian decay
epanechnikov	2	Finite support

library(gridExtra)

fit_sch <- bias_bound_density(x_data, kernel = "schennach")
fit_sinc <- bias_bound_density(x_data, kernel = "sinc")

grid.arrange(
  plot(fit_sch) + ggtitle("Schennach2004 (recommended)"),
  plot(fit_sinc) + ggtitle("Sinc kernel"),
  ncol = 1
)

Extension to regression

The same principles apply to regression. For \(E[Y \mid X = x]\) the Nadaraya-Watson estimator is

\[\hat{m}(x) = \frac{\sum_{i=1}^{n} K_h(x - X_i) Y_i}{\sum_{i=1}^{n} K_h(x - X_i)},\]

and its numerator and denominator biases are bounded separately. In particular, the numerator bound retains the same Fourier-normalization factor:

\[\bar b_{Y;X}(x) = \frac{1}{2\pi}\int_{-\infty}^{\infty} |K^{FT}(h\xi)-1|\min\{B,A|\xi|^{-r}\}\,d\xi.\]

The two bounds are then propagated through the ratio:

# Generate regression data
y_data <- sin(2 * pi * x_data) + rnorm(500, sd = 0.3)

# Estimate with bias bounds
fit_reg <- bias_bound_regression(x_data, y_data, bw = 0.1)

# View smoothness parameters
coef(fit_reg)
#>          A          r          B          h 
#> 22.9169202  2.0000000  0.6374611  0.1000000

Bandwidth selection

For density estimation, the bandwidth is chosen by ordinary least-squares (leave-one-out) cross-validation, the same MSE-oriented criterion used in standard kernel density estimation. This step is entirely separate from the Fourier analysis above: cross-validation selects the bandwidth, and the Fourier envelope is used only afterwards to bound the bias at that bandwidth.

\[h_{CV} = \arg\min_h \left[ \int \hat{f}_h(x)^2 \, dx \;-\; \frac{2}{n} \sum_{i=1}^{n} \hat{f}_{-i,h}(X_i) \right]\]

where \(\hat{f}_{-i,h}\) is the estimate computed without observation \(i\).

For regression, the selector instead minimizes leave-one-out squared Nadaraya–Watson prediction error. It therefore uses both the predictor and response; it does not reuse density cross-validation.

fit_cv <- bias_bound_density(x_data, bw = "cv", kernel = "schennach")
fit_silv <- bias_bound_density(x_data, bw = "silverman", kernel = "normal")
h_cv <- unname(coef(fit_cv)["h"])
h_silv <- unname(coef(fit_silv)["h"])

cat("CV bandwidth:", round(h_cv, 4), "\n")
#> CV bandwidth: 0.1687
cat("Silverman bandwidth:", round(h_silv, 4))
#> Silverman bandwidth: 0.1045

At an MSE-oriented (cross-validation) bandwidth the bias can be large enough to matter. The kernel rounds off peaks and other curvature, so a textbook confidence interval that uses only a standard error can under-cover. The bias bound addresses this directly: it widens the interval by the worst-case allowance and, under Schennach’s assumptions, supports coverage without changing the bandwidth.

The figure below makes this concrete at the optimal bandwidth. The narrow band is the naive interval (standard deviation only); the wider band is the bias-bound interval. Near the peak the true density (dashed line) leaves the naive band but stays inside the bias-bound band:

fit_opt <- bias_bound_density(x_data, bw = h_cv, kernel = "schennach")

# True density: the convolution of two U[0,1] is the triangular density on [0, 2]
tri_density <- function(x) ifelse(x >= 0 & x <= 2, 1 - abs(x - 1), 0)

z <- qnorm(0.975)
band_df <- data.frame(
  x        = fit_opt$x,
  estimate = fit_opt$density,
  truth    = tri_density(fit_opt$x),
  naive_lo = pmax(fit_opt$density - z * fit_opt$std_error, 0),
  naive_hi = fit_opt$density + z * fit_opt$std_error,
  bb_lo    = fit_opt$conf_int$lower,
  bb_hi    = fit_opt$conf_int$upper
)

ggplot(band_df, aes(x)) +
  geom_ribbon(aes(ymin = bb_lo, ymax = bb_hi, fill = "Bias-bound 95% CI")) +
  geom_ribbon(aes(ymin = naive_lo, ymax = naive_hi, fill = "Naive 95% CI (ignores bias)")) +
  geom_line(aes(y = estimate, color = "Estimate"), linewidth = 0.7) +
  geom_line(aes(y = truth, color = "True density"), linetype = "dashed", linewidth = 0.7) +
  scale_fill_manual(values = c("Bias-bound 95% CI" = "#9ECAE1",
                               "Naive 95% CI (ignores bias)" = "#FCAE91")) +
  scale_color_manual(values = c("Estimate" = "#08306B", "True density" = "black")) +
  labs(x = "x", y = "density", fill = NULL, color = NULL) +
  theme(legend.position = "top")

The classical route to valid inference is undersmoothing: shrink the bandwidth until the bias is negligible relative to the standard error. This increases variance. The bias-bound interval instead retains the selected bandwidth and pays an explicit price for bias. The finite-sample comparison is empirical, not guaranteed: across the benchmark cells their average pointwise coverage ranges from 0.991 to 1.000, while their intervals are generally wider, especially in regression, than conventional intervals at one-quarter of the selected bandwidth.

Frequency window and smoothness floor

The default frequency_window = "snr" setting uses a sustained signal-to-noise crossing of the empirical Fourier transform. It is a calibrated finite-sample heuristic; it replaces the original worst-case feasibility test, which often returns no usable frequency at realistic sample sizes. The default threshold is

\[\tau_n = \sqrt{2\log(e^2+n)}.\]

The default smoothness_floor = TRUE solves the theoretical mixed-integer problem with \(r\in\{2,3,\ldots\}\). Setting it to FALSE requests a continuous-slope diagnostic extension; Schennach’s 2020 coverage theorem does not cover that extension. The numerical search has the documented cap smoothness_max = 50. For regression, the package selects and fits the Fourier windows separately for the marginal object (\(Y=1\)) and the response-weighted numerator object, matching the two applications of the theory.

The exact frequency_window = "schennach" rule stops with guidance when no frequency passes its feasibility test. It does not silently substitute the theoretical upper cap. The default SNR heuristic is offered precisely because this empty-set outcome is common at practical sample sizes.

Summary

The bias-bound approach provides:

Bias-aware inference at MSE-oriented bandwidths
Automatic smoothness detection via Fourier analysis
Explicit bias accounting in confidence intervals
An explicit validity-versus-width tradeoff

References

Schennach, S. M. (2020). A Bias Bound Approach to Non-parametric Inference. The Review of Economic Studies, 87(5), 2439-2472. doi:10.1093/restud/rdz065