iRBCrd: IMPROVED ROBUST BIAS CORRECTION FOR RD DESIGNS
Version 2.1 -- October 29, 2026
Website: https://irbc-packages.github.io/

AUTHORS AND CONTACTS
--------------------
Giuseppe Cavaliere
  University of Bologna and Exeter Business School
  giuseppe.cavaliere@unibo.it
Silvia Goncalves
  McGill University
  silvia.goncalves@mcgill.ca
Morten Orregaard Nielsen
  Aarhus University
  mon@econ.au.dk
Edoardo Zanelli
  Aarhus University
  ez@econ.au.dk

1. OVERVIEW
-----------
iRBCrd provides bias-corrected point estimates, analytic confidence intervals,
and bootstrap confidence intervals for:
  - Sharp regression discontinuity (RD).
  - Fuzzy RD.
  - Covariate-adjusted sharp RD.
  - Covariate-adjusted fuzzy RD, adjusting both outcome and treatment equations.
  - Sharp and fuzzy kink RD, with or without covariate adjustment.

Use iRBCrd() for every design. Supply fuzzy for a treatment variable and covs
for covariate adjustment. Add deriv = 1 for a kink design. The function selects
the appropriate implementation from these inputs; deriv = NULL or 0 selects level RD.

This is a distribution of R functions loaded with source(), not an installable
R package or a CRAN release. No license file is included. Bootstrap inference
uses a fixed-regressor wild bootstrap.

2. FILES AND LOADING
--------------------
iRBC_functions.R  Kernels, local-polynomial fits, variance calculations, and
                  bandwidth regularization helpers.
iRBCrd.R          Unified entry point, level RD, and kink RD implementations.
README.txt        This manual.
examples.R        Reproducible level and kink examples.
CHANGELOG.txt     Release overview.
smoke_test.R      Standalone interface and bootstrap checks.
kink_test.R       Independent kink formulas, bootstrap refits, and option checks.
output_test.R     Returned-field and bandwidth-reporting checks for all designs.
bootstrap_weights_test.R  Independent refits for bootstrap weighting choices.
VERIFICATION.txt  Summary of the development regression checks.
CONTACT_SOURCES.txt  Sources used to verify the institutional contacts.

Extract the ZIP into a folder. With R's working directory set to that folder:

  source("iRBC_functions.R")
  source("iRBCrd.R")

Both files are needed, in that order. No additional R packages are required
for estimation. rdrobust is optional for bandwidth selection or comparisons.
Loading these files does not attach or install any additional packages.

To run the examples or standalone checks after setting the working directory:

  source("examples.R")
  source("smoke_test.R")
  source("kink_test.R")
  source("output_test.R")
  source("bootstrap_weights_test.R")

The functions can be used in R, RStudio, or via Rscript. They do not require
RStudio or a compiler. This release was tested in R 4.5.2 on macOS.

3. QUICK START
--------------
  set.seed(17041998)
  n <- 400
  x <- runif(n, -1, 1)
  z <- 0.3 * x + rnorm(n)
  d <- rbinom(n, 1, 0.2 + 0.6 * (x >= 0))
  y <- 1 + x + x^2 + 0.8 * d + 0.5 * z + rnorm(n)
  Z <- cbind(z = z)

  sharp <- iRBCrd(y, x, h = 0.6)
  fuzzy_fit <- iRBCrd(y, x, h = 0.6, fuzzy = d)
  adjusted_sharp <- iRBCrd(y, x, h = 0.6, covs = Z)
  adjusted_fuzzy <- iRBCrd(y, x, h = 0.6, fuzzy = d, covs = Z)

The bandwidth above is illustrative. It is not an estimated optimal bandwidth.
The bootstrap and analytic reportable estimates are identical. Their intervals
can differ, and the bootstrap interval need not be symmetric around that
bias-corrected estimate.

4. FUNCTION INTERFACE
---------------------
  iRBCrd(y, x, h, c = NULL, p = NULL, kernel = NULL,
         fuzzy = NULL, deriv = NULL, covs = NULL,
         res = NULL, level = 95, bwcheck = .bwcheck_default,
         fast = NULL,
         boot.ci = NULL, boot.B = NULL, boot.seed = NULL,
         boot.wild = NULL, boot.keep = NULL,
         q.method = NULL, fuzzy.boot.weights = NULL,
         gamma.boot.weights = NULL)

Design and estimation:
  y, x       Numeric vectors in the same observation order. Prepare finite,
             complete-case data before calling the function.
  h          Required positive bandwidth, scalar or a pair ordered left, right.
  c          Cutoff; default 0. The right side includes x = c.
  p          Local-polynomial order; default 1 for levels and 2 for kinks.
             Use a nonnegative integer for levels, and p >= 1 for kinks.
  kernel     "tri" (default), "epa", or "uni". Long aliases "triangular",
             "epanechnikov", and "uniform" are accepted. A pair is allowed,
             ordered left, right, or named left/right or minus/plus.
  fuzzy      Treatment vector. NULL selects sharp RD. This is not a Boolean.
             Binary, discrete, and continuous treatments are accepted.
  deriv      Derivative order: 0 for discontinuities in levels, 1 for slope
             changes. NULL defaults to 0. Only numeric 0 and 1 are accepted.
  covs       Numeric covariate vector, matrix, or data frame; default NULL.
             One row per observation. Both equations are adjusted in fuzzy RD.
  res        "irbc-hc0" through "irbc-hc3", or "cct-hc0" through "cct-hc3".
             The default is pointwise iRBC-HC3 ("irbc-hc3").
  level      Confidence level in percent; default 95 gives a 95% interval.
             Supply a finite number strictly between 0 and 100.
  q.method   "sample" or "asymptotic". Default "sample" when fast = FALSE;
             "asymptotic" when fast = TRUE.
  fast       Standard-error shortcut; default FALSE. See Section 8.
  bwcheck    Minimum-neighbor regularization rule; default 21. NULL disables it.

Bootstrap:
  boot.ci    TRUE to add a bootstrap interval; default FALSE.
  boot.B     Number of draws; default 999 when boot.ci = TRUE. A supplied
             value must be a positive integer. At least two finite draws are
             needed for an interval; use substantially more for inference.
  boot.seed  Optional numeric seed. Default NULL uses the current R RNG state.
             A supplied seed is passed to set.seed() after integer conversion.
             The prior random-number state is not restored.
  boot.wild  "rademacher" (default), "mammen", or "normal".
  boot.keep  Retain bootstrap statistic draws and weight diagnostics; default FALSE.
  fuzzy.boot.weights
             "fixed" (default) or "bootstrap" in every fuzzy design,
             including kink and covariate-adjusted designs.
  gamma.boot.weights
             "fixed" (default) or "bootstrap" in every covariate-adjusted
             design, including sharp, fuzzy, and kink designs.

5. BANDWIDTHS AND DATA
----------------------
The function does not select an optimal bandwidth. Supply h, for example from
rdrobust::rdbwselect(), with the appropriate fuzzy and covariate inputs in the
selector. For a kink, also set deriv = 1 and the same p in the selector.
An unnamed bandwidth pair is left, right. Named pairs accept
left/right or minus/plus. A scalar applies to both sides. Supply bandwidths
only through h. Kernels follow the same convention through kernel.

  fit <- iRBCrd(y, x, h = c(left = 0.5, right = 0.7),
                kernel = c(left = "tri", right = "epa"))

When using rdrobust for a comparison with rho = 1, explicitly use the same
selected h and set b = h there. For an exact bandwidth comparison, consider
bwcheck = NULL here, provided all needed local fits remain nonsingular.

With bwcheck = k, a local bandwidth is enlarged if necessary to reach the kth
closest observation, or the farthest observation if fewer than k are present.
This counts observations, not distinct running-variable values. The initial
cutoff regularization uses pooled observations for scalar h,
and side-specific observations for paired bandwidth inputs.
Auxiliary fits can enlarge their own bandwidths as well.

For kink RD, cutoff regularization is side-specific even for scalar h. The
reported h.minus and h.plus are the actual slope-fitting bandwidths. These
effective bandwidths also start the side's auxiliary fits. pilot.h records
any further bandwidth enlargement in the level pilots used for the bias.
Kink bias, variance, and bootstrap calculations use the same coefficient maps,
including when a pilot neighborhood is enlarged.

Returned h is a named list of requested minus/plus bandwidths.
h.minus and h.plus are output fields, not inputs; they report the effective
cutoff bandwidths, not the bandwidth at every auxiliary fitting point. This
rule does not guarantee full rank with ties, collinearity, or zero endpoint
kernel weights. Do not trim the data to the cutoff h window before calling
iRBCrd(): auxiliary fits need observations outside that window. Bias and
variance calculations can involve observations up to 2h from the cutoff when
regularization is disabled; pointwise residual fits use their own neighbors.

6. EXTENSIONS
-------------
All extensions use iRBCrd() and support analytic and bootstrap inference.
The design inputs can be combined; bandwidth and kernel conventions remain
as in Section 5.

Kink RD:
  Set deriv = 1 to estimate the right slope minus the left slope at the
  cutoff. The default is local quadratic estimation (p = 2); p = 1 or a
  higher order can also be supplied. Sharp kink reports the slope difference
  in units of y/x, without dividing by a known policy-schedule slope change.
  For that target, rescale y by the known change before estimation. Bias
  correction uses slope-specific Q adjustments; residuals and bootstrap
  outcomes remain in levels.

    kink_fit <- iRBCrd(y, x, h, deriv = 1)

Fuzzy RD:
  Supply the treatment vector through fuzzy. The target is the outcome
  discontinuity divided by the treatment discontinuity, or the corresponding
  ratio of slope changes when deriv = 1. iRBC corrects the linearized ratio
  and retains the covariance between the outcome and treatment components.
  A near-zero first stage prevents ordinary ratio inference; a nonzero
  estimate alone does not establish a strong first stage. Bootstrap weighting
  choices are described in Section 7.

    fuzzy_fit <- iRBCrd(y, x, h, fuzzy = d)

Covariate adjustment:
  Supply covs to add common covariate slopes to the side-specific order-p
  weighted regressions. The outcome discontinuity is adjusted by the
  estimated covariate coefficients times the covariate discontinuities.
  Fuzzy designs adjust both outcome and treatment equations, using separate
  coefficient vectors and retaining their covariance. In kink designs, the
  adjustment uses covariate slope changes. Preserving the unadjusted target
  requires covariate balance at the cutoff: continuity for level RD and no
  slope change for kink RD, alongside the design's identification conditions.

    adjusted_fit <- iRBCrd(y, x, h, covs = Z)

7. BOOTSTRAP
-------------
The running variable, bandwidths, and Q adjustments are held fixed. For local
observations, bootstrap outcomes equal pointwise local-polynomial fitted
values plus HC-adjusted residuals multiplied by wild weights. The same draw
is used for all components of an observation in joint designs, preserving
their dependence. Observations outside the cutoff working window are left
unchanged; the bandwidth is not re-estimated in each bootstrap draw.

For each sharp component j, the bootstrap statistic is
  T_j* = Q_plus [ghat_j,plus*(c) - ghat_j,plus(c)]
       - Q_minus [ghat_j,minus*(c) - ghat_j,minus(c)].

For kink designs, replace the cutoff level fits in this expression with
cutoff slope fits. Bootstrap outcomes and residuals are still constructed
in levels, and Q uses the slope-specific adjustment.

The sharp, fuzzy, and covariate-adjusted statistics use the relevant component
linear combination. By default, fuzzy and covariate weights remain fixed.
The two choices can be changed independently in level and kink designs:
  fuzzy.boot.weights = "bootstrap" uses the bootstrap first-stage discontinuity
  (or slope change) in the denominator of the fuzzy weights. The original
  outcome-to-treatment ratio remains fixed, as in the linearized statistic.
  gamma.boot.weights = "bootstrap" refits the common covariate coefficients
  using bootstrap outcomes and covariates. Fuzzy designs also refit the
  treatment equation. The bootstrap first stage is adjusted using whichever
  covariate coefficients are selected by gamma.boot.weights.
Thus, covariate-adjusted fuzzy designs allow all four fixed/bootstrap
combinations. These options change the bootstrap statistic's combination
weights, not the original analytic estimate, variance, or Q adjustments.
Draws with singular covariate refits or undefined fuzzy weights are omitted
from the interval calculation.

Let se_us denote the standard error of the uncorrected Q-weighted component
combination, and mhat = se_irbc/se_us. The interval uses empirical quantiles
of the bootstrap statistic at
  pnorm(mhat * qnorm(c((100-level)/200, (100+level)/200))).
Quantiles use R's type = 7 interpolation and only finite draws. They are
inverted around the corresponding uncorrected discontinuity or ratio estimate.
The REPORTABLE bootstrap point estimate is nevertheless the same bias-corrected
estimate as for analytic iRBC. This distinction is intentional.

Rademacher weights take -1 or 1 with equal probability. Mammen weights use
the usual mean-zero, variance-one asymmetric two-point distribution. Normal
weights are also available. Bootstrap quantiles have simulation error,
especially at high confidence levels. B = 999 is a default, not a guarantee
of accurate tail quantiles. Check stability across seeds or larger B when
differences between intervals are small.

8. OTHER DESIGN CHOICES
-----------------------
Residuals:
  iRBC-HCk fits a p-th order local polynomial separately on each cutoff side,
  centered at each X_i. The raw residual is Y_i - ghat(X_i). There is no
  additional Q-weighted bias subtraction in this fitted value. The HC
  adjustment uses the pointwise fit's leverage and local sample size.

  CCT-HCk fits a (p+1)-th order local polynomial at the cutoff on each side.
  Evaluate it at X_i, subtract it from Y_i, and apply the HC adjustment using
  that cutoff fit's leverage and local sample size. Both residual choices
  use fits in levels, including for kink designs.

For raw residual e_i and leverage ell_i:
  HC0: e_i
  HC1: e_i / sqrt((n_local - df) / n_local)
  HC2: e_i / sqrt(1 - ell_i)
  HC3: e_i / (1 - ell_i)

Here df = p+1 for iRBC residuals and df = p+2 for CCT residuals; n_local counts
nonzero-kernel observations. The same residuals are used in the variance
estimator and bootstrap data-generating process. In joint designs, residuals
are evaluated componentwise and their cross-products enter the covariance.

"hc0" through "hc3" select iRBC-HCk; "rbc-hc0" through "rbc-hc3" select
CCT-HCk. The returned res field uses "irbc-hck" or "cct-hck", respectively,
with k replaced by 0, 1, 2, or 3.

res = "loo" uses unscaled residuals from the ordinary pointwise auxiliary fits,
without an HC adjustment.

Q and analytic inference:
  q.method = "sample" uses finite-sample Q_n values (the default).
  q.method = "asymptotic" uses kernel-specific population values; keep
  fast = FALSE to retain the full variance calculation. Asymptotic Q requires
  p = 1 for level RD and also supports higher orders for kink RD.
  fast = TRUE uses asymptotic Q and a standard-error shortcut for sharp designs,
  with slope-specific constants for kink RD; joint designs retain full
  covariance matrices. It cannot use sample Q.

The analytic interval is estimate_irbc +/- qnorm((100+level)/200) * se_irbc.

9. RESULTS
-----------
The primary iRBCrd() result fields are:
  estimate_irbc       Bias-corrected point estimate.
  se_irbc             Analytic standard error.
  ci_irbc             Analytic interval, lower then upper.
  level               Confidence level in percent.
  deriv               0 for level RD; 1 for kink RD.
  estimate_irbc_boot  Same point estimate, when boot.ci = TRUE.
  ci_irbc_boot        Bootstrap interval, lower then upper.
  h                  Requested bandwidths, as a named minus/plus list.
  h.minus, h.plus     Effective cutoff bandwidths.
  kernel, kernel.minus, kernel.plus
                     Kernel specifications used on the two sides.
  res                Residual specification.
  boot.B, boot.wild, boot.t
                     Bootstrap details; boot.t is omitted by default.

Set boot.keep = TRUE to retain bootstrap statistics and, where applicable,
component statistics, covariate coefficients, and combination weights for
individual draws. This does not change point estimates or confidence intervals.

Design-specific details include uncorrected discontinuities, bias components,
first-stage estimates, covariate slopes, covariance matrices, and sharp.*
component objects. Residual arrays are computational inputs over their local
support, not a full-sample diagnostic for fitted values.

Kink results also report deriv = 1, p, design, linear.weights, kink.weights,
and pilot.h. Component estimates are slope changes. Vectors in kink.weights
and the residuals follow the original observation order within each side.
Sharp kink returns tau.hat; fuzzy kink returns zeta.hat.

10. REFERENCES
---------------
Cavaliere, G., Goncalves, S., Nielsen, M. O., and Zanelli, E. (2026):
  Improved Inference for Nonparametric Regression and Regression-Discontinuity
  Designs. arXiv:2512.00566. https://arxiv.org/abs/2512.00566.
Calonico, S., Cattaneo, M. D., Farrell, M. H., and Titiunik, R. (2019):
  Regression Discontinuity Designs Using Covariates, Review of Economics
  and Statistics, 101(3), 442-451.

The code is supplied as R functions. No install.packages("iRBCrd") or
library(iRBCrd) step applies to this distribution. See CHANGELOG.txt for the
release overview and VERIFICATION.txt for tests performed before preparing
the ZIP.
