Package {tvbounds}


Type: Package
Title: Sensitivity Analysis and Bounds under Total Variation Neighborhoods
Version: 0.1.1
Description: Implements the sensitivity analysis framework of Palomba (2026) "Sensitivity Analysis in Population Shares" https://filippopalomba.github.io/#jmp for randomized experiments with attrition, counterfactuals in structural models, and recentered instrumental variables. Computes and plots sensitivity bounds together with their confidence intervals and robustness summary measures.
License: MIT + file LICENSE
Encoding: UTF-8
Depends: R (≥ 4.1.0)
Imports: ggplot2 (≥ 3.4.0), stats, utils, withr
Suggests: JuliaCall, testthat (≥ 3.1.5), knitr, rmarkdown
SystemRequirements: For tvbounds_counterfactual(): Julia (>= 1.9) and the Artelys KNITRO solver with a valid license.
VignetteBuilder: knitr
RoxygenNote: 7.3.3
Config/testthat/edition: 3
NeedsCompilation: no
Packaged: 2026-08-30 17:10:06 UTC; fpalomba
Author: Filippo Palomba [aut, cre, cph]
Maintainer: Filippo Palomba <fpalomba@princeton.edu>
Repository: CRAN
Date/Publication: 2026-09-10 14:20:02 UTC

tvbounds: Sensitivity Analysis and Bounds under Total Variation Neighborhoods

Description

Implements the sensitivity analysis framework of Palomba (2026), "Sensitivity Analysis in Population Shares". The estimand is an expectation \mathbb{E}_{P}[g(Z;\theta)] of a known integrand g under a distribution P of the data Z. Rather than committing to a single baseline distribution P_{*}, the package lets P range over a robustness set

\mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta) = \{P : \mathbb{E}_{P}[m(Z;\theta)] \in \mathcal{M}(\rho),\ D_{\phi}(P \| P_{*}) \le \delta\},

collecting the distributions that remain compatible with the moment restrictions and lie within a divergence budget \delta of the baseline, and reports the resulting sensitivity bounds

\inf_{\theta \in \Theta}\ \inf_{P \in \mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta)} \mathbb{E}_{P}[g(Z;\theta)] \qquad \mathrm{and} \qquad \sup_{\theta \in \Theta}\ \sup_{P \in \mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta)} \mathbb{E}_{P}[g(Z;\theta)].

Details

The sensitivity parameter is the budget \delta, supplied through the delta argument. Under the total variation entropy \phi_{\mathsf{TV}}(s) = |s - 1| / 2 the divergence is the total variation distance \mathsf{TV}(P,P_{*}) and \delta \in [0,1] bounds the fraction of baseline probability mass that may be misspecified; under the contamination neighborhood \mathcal{C}_{\kappa}(P_{*}) = \{P = \kappa P_{*} + (1 - \kappa) R\} the perturbed distribution is a mixture of the baseline distribution and an arbitrary distribution R.

Applications

Three ready-made interfaces cover the paper's empirical applications:

Reporting

All estimators return a common tvbounds object carrying the bound paths \underline{\tau}(\delta) and \overline{\tau}(\delta) over a grid of budgets. tvbounds_plot() (or plot()) displays the bounds against the budget, and tvbounds_summary() (or summary()) computes the summary measures of the paper: the breakdown budget \delta_{b}(\tau_{\star}) and its certified counterpart, the shadow price of robustness \underline{\eta}(\delta), the robustness standard error \varsigma_{b} = \sigma(\delta_{b}) / \underline{\eta}(\delta_{b}), and the certification frontier n^{\star}(\delta;\alpha).

KNITRO requirement

tvbounds_counterfactual() relies on Julia and the commercial Artelys KNITRO solver, which requires a valid license; see that function's help page for details.

Author(s)

Maintainer: Filippo Palomba fpalomba@princeton.edu [copyright holder]

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Borusyak, K. and Hull, P. (2023). "Nonrandom Exposure to Exogenous Shocks." Econometrica, 91(6), 2155–2185.

Christensen, T. and Connault, B. (2023). "Counterfactual Sensitivity and Robustness." Econometrica, 91(1), 263–298.

Lee, D. S. (2009). "Training, Wages, and Sample Selection: Estimating Sharp Bounds on Treatment Effects." Review of Economic Studies, 76(3), 1071–1102.


Plot sensitivity bounds against the budget

Description

Displays the bounds stored in a tvbounds object as functions of the budget \delta of the total-variation (or contamination, or divergence) neighborhood, supplied through the delta column of x$bounds and drawn on the horizontal axis. The region between the lower bound path \underline{\tau}(\delta) and the upper bound path \overline{\tau}(\delta) is shaded, the outer confidence band (when the object carries one) is drawn as a lighter ribbon delimited by dashed lines, a dashed horizontal line marks the reference value \tau_{\star}, a point marks the baseline estimate at \delta = 0, and a dotted vertical line marks the plug-in breakdown budget \widehat{\delta}_{b} when it is interior to the budget grid.

Usage

## S3 method for class 'tvbounds'
plot(x, ...)

tvbounds_plot(
  x,
  bands = TRUE,
  baseline = TRUE,
  breakdown = TRUE,
  tau_star = 0,
  color = "#1F4E79",
  xlab = NULL,
  ylab = NULL,
  title = NULL,
  log_x = FALSE,
  ...
)

Arguments

x

A tvbounds object returned by tvbounds_attrition(), tvbounds_counterfactual(), or tvbounds_riv().

...

For tvbounds_plot(): currently unused, accepted for compatibility with the plot() generic. For the plot() method: further arguments forwarded to tvbounds_plot().

bands

Logical; draw the outer confidence band when the object carries one (columns ci_lower/ci_upper of x$bounds). Sides without a band are drawn without one. Default TRUE.

baseline

Logical; mark the baseline point estimate, that is, the value of the estimand under P_{*}, at \delta = 0. Ignored when log_x = TRUE, since \delta = 0 cannot be placed on a logarithmic axis. Default TRUE.

breakdown

Logical; draw a dotted vertical line, with a label, at the plug-in breakdown budget \widehat{\delta}_{b} when it is interior to the budget grid. Default TRUE.

tau_star

Reference value \tau_{\star} of the estimand against which robustness is judged; drawn as a dashed horizontal line (default 0).

color

Colour of the bounds, ribbons, and breakdown mark (default "#1F4E79", the paper's blue).

xlab, ylab, title

Axis labels and plot title. NULL (the default) uses expression(delta) for the horizontal axis, the object's estimand_label for the vertical axis, and no title.

log_x

Logical; use a logarithmic budget axis. Rows with delta <= 0 are dropped with a warning. Default FALSE.

Details

The breakdown budget drawn by breakdown = TRUE is the plug-in breakdown budget of tvbounds_summary(), the estimated counterpart of \delta_{b}(\tau_{\star}) = \inf\{\delta : \underline{\tau}(\delta) \le \tau_{\star} \le \overline{\tau}(\delta)\}: the first budget at which the bound path adjacent to \tau_{\star} reaches it — the lower path when the baseline point estimate exceeds tau_star, the upper path otherwise. The line is annotated with the value of \widehat{\delta}_{b}. No line is drawn when the breakdown budget is censored at the right endpoint of the grid, when it sits at the left endpoint, or when the baseline point estimate is missing.

Rows of x$bounds with missing bound values (e.g. censored or infeasible budgets) are omitted from the corresponding layer. For the counterfactual application, whose divergence budgets may exceed one, log_x = TRUE switches to a logarithmic budget axis.

Value

A ggplot object, which prints to the active graphics device and can be modified further with ggplot2 layers.

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

See Also

tvbounds_summary() for the numerical summary measures; plot() dispatches here for tvbounds objects.

Examples

set.seed(123)
n <- 400
d <- rbinom(n, 1, 0.5)
s <- rbinom(n, 1, plogis(0.5 + 0.5 * d))
y <- ifelse(s == 1, 0.3 * d + rnorm(n), NA)
dat <- data.frame(y = y, d = d, s = s)
fit <- tvbounds_attrition(dat, outcome = "y", treatment = "d",
                          response = "s", delta = seq(0, 1, by = 0.1),
                          B = 100, seed = 1)
tvbounds_plot(fit)

# Bounds only, no confidence band and no breakdown mark:
tvbounds_plot(fit, bands = FALSE, breakdown = FALSE)

Print a tvbounds object

Description

Compact display of a tvbounds object: the application, the neighborhood over which the sensitivity bounds were computed, the sample size n, the baseline point estimate — the value of the estimand under the baseline distribution P_{*}, reported at \delta = 0 — the grid of budgets \delta on which the bound paths \underline{\tau}(\delta) and \overline{\tau}(\delta) were evaluated, and whether inference is attached.

Usage

## S3 method for class 'tvbounds'
print(x, digits = 3, ...)

Arguments

x

A tvbounds object.

digits

Number of significant digits (default 3).

...

Further arguments; ignored.

Value

x, invisibly.

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Examples

set.seed(1)
n <- 200
d <- rbinom(n, 1, 0.5)
s <- rbinom(n, 1, plogis(0.4 + 0.4 * d))
y <- ifelse(s == 1, 0.3 * d + rnorm(n), NA)
fit <- tvbounds_attrition(data.frame(y = y, d = d, s = s),
                          outcome = "y", treatment = "d", response = "s",
                          delta = seq(0, 1, by = 0.1), bootstrap = FALSE)
print(fit)

Print a tvbounds summary

Description

Compact display of the summary measures computed by tvbounds_summary(): the plug-in and certified breakdown budgets \delta_{b}(\tau_{\star}) and \widehat{\delta}_{b}^{\,\mathsf{C}}(\alpha), the shadow price of robustness \underline{\eta}(\delta) = -\underline{\tau}'(\delta), the robustness standard error \varsigma_{b} = \sigma(\delta_{b})/\underline{\eta}(\delta_{b}), and the certification frontier n^{\star}(\delta;\alpha), followed by any notes on measures that could not be computed. The reference value of the estimand, \tau_{\star}, is shown as tau_star, matching the argument name.

Usage

## S3 method for class 'tvbounds_summary'
print(x, digits = 3, ...)

Arguments

x

A "tvbounds_summary" object.

digits

Number of significant digits (default 3).

...

Further arguments; ignored.

Value

x, invisibly.

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Examples

set.seed(1)
n <- 200
d <- rbinom(n, 1, 0.5)
s <- rbinom(n, 1, plogis(0.4 + 0.4 * d))
y <- ifelse(s == 1, 0.3 * d + rnorm(n), NA)
fit <- tvbounds_attrition(data.frame(y = y, d = d, s = s),
                          outcome = "y", treatment = "d", response = "s",
                          delta = seq(0, 1, by = 0.1), bootstrap = FALSE)
# Without bootstrap draws the certified measures degrade to NA, with a
# note explaining why:
summary(fit)

Summary measures for total-variation sensitivity bounds

Description

Computes the summary measures of Palomba (2026) for a tvbounds object: the plug-in and certified breakdown budgets, the shadow price of robustness, the robustness standard error, and the certification frontier. The budget delta is the sensitivity parameter \delta of the total-variation (or contamination, or divergence) neighborhood over which the bounds were computed.

Usage

## S3 method for class 'tvbounds'
summary(object, ...)

tvbounds_summary(
  object,
  delta = NULL,
  tau_star = 0,
  level = 0.95,
  cost_per_unit = 50,
  jump = 0.05,
  direction = c("auto", "lower", "upper")
)

Arguments

object

A tvbounds object returned by tvbounds_attrition(), tvbounds_counterfactual(), or tvbounds_riv().

...

For the summary() method: further arguments forwarded to tvbounds_summary().

delta

Optional evaluation budget \delta (a single nonnegative number). NULL (default) evaluates the measures at the plug-in breakdown budget \delta_{b}(\tau_{\star}), with the censoring fallback described in Details.

tau_star

Reference value \tau_{\star} of the estimand against which robustness is judged (default 0).

level

Confidence level 1 - \alpha for the two-sided critical value z_{1-\alpha/2} used by the normal-approximation breakdown budget and the certification frontier (default 0.95, that is, qnorm(0.975)). The certified breakdown itself is read off the band stored in the object.

cost_per_unit

Marginal cost of one additional sampled unit, used for the cost equivalent of the certification frontier (default 50, the paper's USD benchmark from the 3ie closed-grant portfolio).

jump

Increase in the certified budget that the frontier prices, so that n_star is evaluated at \delta + \mathrm{jump} (default 0.05).

direction

One of "auto", "lower", "upper"; see Details.

Details

Write \underline{\tau}(\delta) and \overline{\tau}(\delta) for the lower and upper sensitivity bounds on the estimand at budget \delta. The first is nonincreasing and the second nondecreasing in the budget, and both collapse at \delta = 0 to the value of the estimand under the baseline distribution P_{*}. Robustness is judged against a reference value \tau_{\star} of the estimand, supplied through tau_star, and typically zero when it is the sign of an effect rather than its magnitude that is of interest.

All measures are computed on the signed bound path adjacent to \tau_{\star}, namely \underline{\tau}(\delta) - \tau_{\star} when direction is "lower" and \tau_{\star} - \overline{\tau}(\delta) when it is "upper", so that in both cases the path starts positive when the conclusion holds at the baseline. The two directions are thus treated symmetrically, as the paper treats them by replacing the integrand g by -g; the displays below are written for the lower path. With direction = "auto" (the default) the lower path is used when the baseline point estimate exceeds tau_star and the upper path otherwise.

The breakdown budget is the smallest budget at which the bounds cease to exclude the reference value,

\delta_{b}(\tau_{\star}) = \inf\{\delta \in [0,1] : \underline{\tau}(\delta) \le \tau_{\star} \le \overline{\tau}(\delta)\},

with the convention that the infimum over the empty set equals one. That convention is exactly the package's censoring rule: when the estimated path never reaches \tau_{\star} on the supplied budget grid, the breakdown is reported at the right endpoint of the grid, which is one for the total-variation and contamination neighborhoods, and is flagged by censored = TRUE rather than recorded as an inequality.

Sampling uncertainty is accounted for by the certified breakdown budget,

\widehat{\delta}_{b}^{\,\mathsf{C}}(\alpha) = \inf\left\{\delta \in [0,1] : \widehat{\underline{\tau}}_{n}(\delta) - \frac{z_{1-\alpha/2}\,\widehat{\sigma}_{n}(\delta)}{\sqrt{n}} \le \tau_{\star}\right\},

the largest budget at which the conclusion survives sampling uncertainty. Here \widehat{\underline{\tau}}_{n}(\delta) estimates the lower bound path, \sigma(\delta) is the asymptotic standard deviation of that estimator in units of the estimand and \widehat{\sigma}_{n}(\delta) its estimator, and z_{1-\alpha/2} is the two-sided normal critical value at the confidence level level.

The shadow price of robustness is the marginal cost, in units of the estimand, of one further unit of budget, \underline{\eta}(\delta) = -\underline{\tau}'(\delta); it is nonnegative and nonincreasing, and the package estimates it by central differences on the budget grid. Dividing the sampling standard deviation of the bound by the shadow price converts it from units of the estimand into units of the budget, which yields the robustness standard error

\varsigma_{b} = \frac{\sigma(\delta_{b})}{\underline{\eta}(\delta_{b})}.

The certification frontier is the sample size at which the population counterpart of the certified-breakdown inequality just clears the reference value at budget \delta,

n^{\star}(\delta;\alpha) = \frac{z^{2}_{1-\alpha/2}\,\sigma^{2}(\delta)}{(\underline{\tau}(\delta) - \tau_{\star})^{2}}.

It is real-valued, so the smallest certifying integer sample size is \lfloor n^{\star}(\delta;\alpha)\rfloor + 1. The frontier diverges as the budget approaches the breakdown budget and is meaningful only below it: at and past the breakdown budget no sample size certifies the conclusion, and the frontier is reported as NA. Its semi-elasticity \mathrm{d}\log n^{\star}(\delta;\alpha)/\mathrm{d}\delta, the certification elasticity, gives the rate at which the required sample size grows with the budget; cost_per_pp prices a discrete version of it.

The one-row data frame measures reports the following, alongside the symbol each one corresponds to in the paper:

The evaluation budget for eta, se, and varsigma is the delta argument when supplied; when delta = NULL (the default) it is the plug-in breakdown budget, replaced by the certified breakdown budget when the plug-in breakdown is censored, as in the paper. The frontier is priced at the certified breakdown (default) or at the supplied delta.

Objects without inference (the recentered-IV and counterfactual applications carry none by design) degrade gracefully: the plug-in measures are reported, the certified and frontier measures are NA, and the notes field explains why.

Value

An object of class "tvbounds_summary": a list with a one-row data frame measures (columns label, direction, n, point, tau_star, delta_b, censored, delta_b_ci, censored_ci, delta_b_ci_norm, delta_eval, eta, se, varsigma, varsigma_sc, frontier_at, n_cur, n_star, delta_n, cost_per_pp) and metadata fields (application, estimand_label, neighborhood, divergence, direction, tau_star, delta_eval, level, band_level, zc, cost_per_unit, jump, has_se, has_ci, notes, call). Details gives the symbol of the paper each column of measures corresponds to. Printed compactly by print.tvbounds_summary().

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

See Also

tvbounds_plot() to display the bounds; summary() dispatches here for tvbounds objects.

Examples

set.seed(123)
n <- 400
d <- rbinom(n, 1, 0.5)
s <- rbinom(n, 1, plogis(0.5 + 0.5 * d))
y <- ifelse(s == 1, 0.3 * d + rnorm(n), NA)
dat <- data.frame(y = y, d = d, s = s)
fit <- tvbounds_attrition(dat, outcome = "y", treatment = "d",
                          response = "s", delta = seq(0, 1, by = 0.1),
                          B = 100, seed = 1)
tvbounds_summary(fit)

# Evaluate the measures at a chosen budget instead of the breakdown:
tvbounds_summary(fit, delta = 0.2)

Sensitivity bounds for randomized experiments with attrition

Description

Computes sensitivity bounds on the average treatment effect for the always-observed subpopulation of a randomized experiment with attrition, following Palomba (2026). Writing P_0 for the distribution of the data, the estimand is

\tau_0 := \mathbb{E}_{P_0}[Y(1) - Y(0) \mid S(0) = 1, S(1) = 1],

the average treatment effect on the units that respond under either arm. The bounds are indexed by a budget \delta \in [0, 1], supplied through delta, which caps the total variation distance between the outcome distribution of the compliers (units that respond only under treatment) and that of the always-observed units, \mathsf{TV}(P_{\mathsf{C}} \,\|\, P_{\mathsf{AO}}) \le \delta. At \delta = 0 the two distributions coincide and the bounds collapse to the baseline difference in means among respondents (the estimand under missingness completely at random); at \delta = 1 the restriction is vacuous and they equal the trimming bounds of Lee (2009). Optionally computes a nonparametric bootstrap (with clustering) for standard errors and a percentile confidence band, and covariate-pooled bounds that allocate a single budget optimally across covariate cells.

Usage

tvbounds_attrition(
  data,
  outcome,
  treatment,
  response,
  covariates = NULL,
  delta = seq(0, 1, by = 0.01),
  neighborhood = c("tv", "contamination"),
  bootstrap = TRUE,
  B = 1000,
  cluster = NULL,
  level = 0.95,
  min_obs = 5,
  seed = NULL,
  verbose = FALSE
)

Arguments

data

A data frame containing the columns named by outcome, treatment, response, and (optionally) covariates and cluster.

outcome

String; name of the numeric column holding the outcome Y. May be NA for non-respondents (and for the occasional respondent with item non-response, which is dropped from the outcome samples).

treatment

String; name of the binary 0/1 column holding the treatment D (1 = treated).

response

String; name of the binary 0/1 column holding the response indicator S (1 = outcome observed).

covariates

Optional character vector of column names to stratify on, forming the discrete covariate X; cells are formed by their interaction. The covariate columns must be free of missing values — recode missing values into an explicit category first. Default NULL (no stratification).

delta

Numeric vector of budget values \delta \in [0, 1] at which the bounds are evaluated; sorted and de-duplicated internally. The baseline (\delta = 0) and Lee (\delta = 1) endpoints are always computed internally even when absent from delta. Default seq(0, 1, by = 0.01).

neighborhood

Either "tv" (total variation, default) or "contamination"; see the Neighborhoods section.

bootstrap

Logical; compute bootstrap standard errors and the percentile confidence band. Default TRUE.

B

Number of bootstrap replications B. Default 1000.

cluster

Optional string; name of a cluster identifier column. When supplied, the bootstrap resamples whole clusters with replacement. Default NULL (units resampled independently).

level

Confidence level of the percentile band, 1 - \alpha in the notation of the Inference section. Default 0.95.

min_obs

Minimum number of observed outcomes per arm for a covariate cell to be retained. Default 5.

seed

Optional integer seed for the bootstrap. When non-NULL, the random-number-generator state is set locally and restored on exit, so the call has no side effect on the caller's RNG. Default NULL.

verbose

Logical; emit progress messages. Default FALSE.

Value

An object of class c("tvbounds_attrition", "tvbounds"): a list with components

Setup and estimand

Let D \in \{0, 1\} be the binary treatment and, for d \in \{0, 1\}, let Y(d) be the potential outcome and S(d) \in \{0, 1\} the potential response indicator. Their realized counterparts are

Y = D Y(1) + (1 - D) Y(0), \qquad S = D S(1) + (1 - D) S(0),

so that the observed data are (YS, S, D): the outcome is recorded only when S = 1. The response pattern (S(0), S(1)) partitions the population into always-observed units (S(0) = 1, S(1) = 1), never-observed units (S(0) = 0, S(1) = 0), compliers (S(0) = 0, S(1) = 1) and defiers (S(0) = 1, S(1) = 0); the estimand \tau_0 is the average treatment effect on the first of these groups.

Two assumptions are maintained. Random assignment enters as (S(0), S(1)) \perp\!\!\!\perp (Y(0), Y(1)), labelled (MCAR), and the monotonicity condition of Lee (2009), S(1) \ge S(0) almost surely, labelled (Mono) — treatment never causes a unit that would respond under control to attrit — rules out defiers. Under (Mono) the outcome distribution of the observed treated, P_{\mathsf{T}}, is a mixture of the complier and always-observed outcome distributions,

P_{\mathsf{T}} = \pi P_{\mathsf{C}} + (1 - \pi) P_{\mathsf{AO}}, \qquad \pi = 1 - \frac{r_0}{r_1},

where r_1 := P_0[S = 1 \mid D = 1] and r_0 := P_0[S = 1 \mid D = 0] are the arm-specific response rates and \pi is the complier share. The always-observed control mean is identified, \mu^{\mathsf{AO}}(0) = \mathbb{E}_{P_0}[Y \mid D = 0, S = 1], whereas the always-observed treated mean \mu^{\mathsf{AO}}(1) is only partially identified; the bounds on \tau_0 follow by subtracting the identified control mean.

The complier share is estimated by \widehat{\pi} = \max\{1 - \widehat{r}_0 / \widehat{r}_1, 0\} and returned as details$p_star, and the estimated response rates (\widehat{r}_0, \widehat{r}_1) as details$response_rate. A negative unconstrained estimate of \pi is sampling noise under (Mono) and is projected to zero, in which case the bounds collapse to the baseline difference in means at every budget.

Neighborhoods

Two robustness sets are available through neighborhood. Both are indexed by the budget \delta and both restrict the unobserved complier outcome distribution P_{\mathsf{C}} relative to the unobserved always-observed outcome distribution P_{\mathsf{AO}}:

Under total variation the rescaling identity \mathsf{TV}(P_{\mathsf{C}} \,\|\, P_{\mathsf{T}}) = (1 - \pi)\,\mathsf{TV}(P_{\mathsf{C}} \,\|\, P_{\mathsf{AO}}) recenters the restriction on the identified distribution P_{\mathsf{T}}, so that the candidate complier distributions Q, among which P_{\mathsf{C}} lies, range over the robustness set

\mathcal{Q}_{\mathsf{C}}(\delta) := \{Q \in \Delta(\mathcal{Y}) : \mathsf{TV}(Q \,\|\, P_{\mathsf{T}}) \le (1 - \pi)\delta, \ \pi Q \le P_{\mathsf{T}}\},

where \Delta(\mathcal{Y}) denotes the distributions on the outcome space and the second restriction is the mixture structure of the observed treated arm. The resulting sensitivity bounds \underline{\tau}(\delta) and \overline{\tau}(\delta) on \tau_0 are available in closed form as trimmed means of P_{\mathsf{T}}. Writing F^{-1}_{P_{\mathsf{T}}} for the quantile function of the observed treated outcomes and

s_{\mathsf{L}}(\delta) := F^{-1}_{P_{\mathsf{T}}}(\pi \delta), \qquad s_{\mathsf{U}}(\delta) := F^{-1}_{P_{\mathsf{T}}}(1 - (1 - \pi) \delta),

the upper bound is

\overline{\tau}(\delta) = \mathbb{E}_{P_{\mathsf{T}}}[Y \mathbf{1}\{s_{\mathsf{L}}(\delta) < Y < s_{\mathsf{U}}(\delta)\}] + \frac{1}{1 - \pi} \mathbb{E}_{P_{\mathsf{T}}}[Y \mathbf{1}\{Y \ge s_{\mathsf{U}}(\delta)\}] - \mu^{\mathsf{AO}}(0),

and the lower bound is obtained symmetrically, trimming at t_{\mathsf{L}}(\delta) := F^{-1}_{P_{\mathsf{T}}}((1 - \pi) \delta) and t_{\mathsf{U}}(\delta) := F^{-1}_{P_{\mathsf{T}}}(1 - \pi \delta).

Under contamination, combining the same mixture identity with P_{\mathsf{C}} \ge (1 - \delta) P_{\mathsf{AO}} pins the density w := \mathrm{d}P_{\mathsf{C}} / \mathrm{d}P_{\mathsf{T}} between (1 - \delta) / (1 - \delta \pi) and 1 / \pi, and the bounds are again trimmed means of P_{\mathsf{T}}, now with effective trimming mass \delta \pi. The contamination neighborhood is contained in the total variation one at every budget, so its bounds are weakly tighter, and the two families share the same endpoints: the baseline at \delta = 0 and, at \delta = 1, the Lee (2009) bounds \underline{\tau}_{\mathsf{Lee}} = \mathbb{E}_{P_{\mathsf{T}}}[Y \mid Y \le y_{1 - \pi}] - \mu^{\mathsf{AO}}(0) and \overline{\tau}_{\mathsf{Lee}} = \mathbb{E}_{P_{\mathsf{T}}}[Y \mid Y \ge y_{\pi}] - \mu^{\mathsf{AO}}(0), where y_u := F^{-1}_{P_{\mathsf{T}}}(u).

Both bound families are monotone in the budget by construction: the robustness sets are nested in \delta.

Covariates

When covariates is supplied, units are stratified on the interaction of the covariate columns, which plays the role of a discrete covariate X with support \mathcal{X}. Cells with fewer than min_obs observed outcomes in either arm are dropped with a warning, and the retained cells are weighted by their control-respondent shares, which under (Mono) are the covariate distribution of the always-observed population, P_{X \mid D = 0, S = 1} = P_{X \mid \mathsf{AO}}. Within a cell the complier share \pi(x) and the observed treated outcome distribution P_{\mathsf{T}}(x) are identified, and the cell-level construction is the one above.

Two ways of spending the budget across cells are distinguished. The within-stratum ("pointwise") restriction imposes \mathsf{TV}(P_{\mathsf{C}}(x) \,\|\, P_{\mathsf{AO}}(x)) \le \delta in every cell separately, giving one robustness set \mathcal{Q}^{\mathsf{pw}}_{\mathsf{C}}(\delta; x) per cell, whereas the pooled restriction caps only the average departure,

\int_{\mathcal{X}} \mathsf{TV}(P_{\mathsf{C}}(x) \,\|\, P_{\mathsf{AO}}(x)) \, \mathrm{d}P_{X \mid \mathsf{AO}}(x) \le \delta,

and so allows heterogeneity across cells inside the single robustness set \mathcal{Q}_{\mathsf{C},X}(\delta). The pooled restriction is the weaker of the two, so \underline{\tau}_{X}(\delta) \le \underline{\tau}^{\mathsf{pw}}_{X}(\delta) and \overline{\tau}^{\mathsf{pw}}_{X}(\delta) \le \overline{\tau}_{X}(\delta), with equality at \delta = 0 and at \delta = 1, where both collapse to the covariate Lee (2009) bounds \underline{\tau}_{\mathsf{Lee},X} and \overline{\tau}_{\mathsf{Lee},X}.

For neighborhood = "tv" the reported bounds are the pooled (joint) bounds \underline{\tau}_{X}(\delta) and \overline{\tau}_{X}(\delta): the budget allocation t : \mathcal{X} \to \mathbb{R}_+ subject to \mathbb{E}_{P_{\mathsf{T}}}[t(X)] \le (1 - \pi) \delta is solved exactly by a greedy fill over the stratum value functions V(t; x), which are piecewise linear in the cell budget and concave for the upper bound and convex for the lower one. The within-stratum bounds \underline{\tau}^{\mathsf{pw}}_{X}(\delta) and \overline{\tau}^{\mathsf{pw}}_{X}(\delta) are returned in details$pooled$pw for reference. For neighborhood = "contamination" the reported bounds impose the common budget delta within every retained cell and aggregate; they remain weakly inside the total variation bounds at every budget.

Inference

The bootstrap resamples the full observation W = (YS, S, D), together with the covariates, with replacement — whole clusters when cluster is supplied — and recomputes the entire bounds curve on each replicate, yielding \widehat{\overline{\tau}}^{(b)}(\delta), b = 1, \dots, B. Each replicate therefore redraws the arm-specific response rates and hence \widehat{\pi}, so the reported standard errors carry the estimation uncertainty in \pi, which the naive variance \widehat{\sigma}^2_{\mathsf{naive}}(\delta) omits by treating \widehat{\pi} as known; in the paper's influence function \psi_{\mathsf{full}}(W; \delta) this uncertainty is the term \varkappa(\delta) \psi_{\pi}(W). The reported upper_se is the standard deviation of the draws \widehat{\overline{\tau}}^{(b)}(\delta) across replicates, that is the paper's \widehat{\sigma}_{\mathsf{boot}}(\delta) divided by \sqrt{n} for a sample of size n, and lower_se is its counterpart for the lower bound. The reported confidence band is the percentile band: writing \alpha for the value of 1 - level, ci_lower is the \alpha / 2 quantile of the lower-bound draws and ci_upper the 1 - \alpha / 2 quantile of the upper-bound draws, the outer envelope of the identified set. Replicates on which the bounds cannot be computed are dropped and counted (a warning reports their number).

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Lee, D. S. (2009). "Training, Wages, and Sample Selection: Estimating Sharp Bounds on Treatment Effects." Review of Economic Studies, 76(3), 1071-1102.

Huber, P. J. (1964). "Robust Estimation of a Location Parameter." Annals of Mathematical Statistics, 35(1), 73-101.

See Also

tvbounds_plot() and tvbounds_summary() for reporting.

Examples

set.seed(123)
n <- 400
d <- rbinom(n, 1, 0.5)
s <- rbinom(n, 1, ifelse(d == 1, 0.9, 0.7))
y <- ifelse(s == 1, rnorm(n, mean = 0.3 * d), NA)
x <- rbinom(n, 1, 0.5)
dat <- data.frame(y = y, d = d, s = s, x = x)

## Total variation bounds with a small bootstrap
fit <- tvbounds_attrition(dat, outcome = "y", treatment = "d",
  response = "s", delta = seq(0, 1, by = 0.1), B = 50, seed = 1)
fit$bounds
fit$details$lee

## Contamination neighborhood, no inference: weakly tighter bounds
fit_c <- tvbounds_attrition(dat, outcome = "y", treatment = "d",
  response = "s", delta = seq(0, 1, by = 0.1),
  neighborhood = "contamination", bootstrap = FALSE)
all(fit_c$bounds$lower >= fit$bounds$lower - 1e-12)

## Covariate-pooled bounds
fit_x <- tvbounds_attrition(dat, outcome = "y", treatment = "d",
  response = "s", covariates = "x", delta = seq(0, 1, by = 0.1),
  bootstrap = FALSE)


Control options for the counterfactual solver

Description

Constructs the list of tuning options consumed by tvbounds_counterfactual(). The defaults reproduce the settings used for the counterfactual application in Palomba (2026).

Usage

tvbounds_control(
  maxsolves = 10,
  startptrange = 0.01,
  use_optim = FALSE,
  time_limit = 60,
  iterations = 100,
  outer_iterations = 3,
  inner_opt = NULL,
  outer_opt = NULL,
  knitro_options = list(),
  eta_min = 1e-120,
  lower_limit = -10,
  psi_tv_eps = 1e-04,
  tvac_tau = 0.001,
  tvmix_tau = 0.001,
  purekl_acap = 500,
  tvmix_kappa = NULL
)

Arguments

maxsolves

Integer, number of multi-start restarts of the outer optimization over the structural parameter \theta for each budget and side.

startptrange

Positive scalar; restarts after the first perturb the initial \theta uniformly on ⁠[-startptrange, startptrange]⁠ coordinate-wise (clipped to the parameter box).

use_optim

Logical; if TRUE the outer optimization uses Optim.jl (projected L-BFGS with the analytic envelope-theorem gradient) instead of a nested KNITRO solve. The Optim fallback avoids nested KNITRO contexts, which segfault with some KNITRO.jl versions.

time_limit

Positive scalar, wall-clock limit in seconds per outer Optim run (used only when use_optim = TRUE).

iterations

Positive integer, inner iteration limit per outer Optim run (used only when use_optim = TRUE).

outer_iterations

Positive integer, number of Fminbox outer iterations per Optim run (used only when use_optim = TRUE).

inner_opt

Path to a KNITRO option file for the inner (dual) problem, or NULL for the shipped default inner.opt.

outer_opt

Path to a KNITRO option file for the outer problem (over the structural parameter \theta), or NULL for the shipped default outer.opt. The package also ships outer_fast.opt (analytic envelope gradient, looser tolerances, suited to plotting grids) and outer_boot_tv.opt; point this argument at system.file("julia", "opt", "outer_fast.opt", package = "tvbounds") to use them.

knitro_options

Named list of individual KNITRO options (e.g. list(maxit = 500, outlev = 2)). These are merged into both the inner and the outer option files by writing merged copies to tempdir(); entries override options already present and are appended otherwise. Values must be length-one character, numeric, or logical (logicals are written as 0/1).

eta_min

Positive scalar, lower bound on the dual variable \eta, the multiplier pricing the divergence budget (kept strictly positive so the perspective function in the dual objective is well defined).

lower_limit

Scalar; inner (dual) objective values at or below this threshold are treated as unbounded below (the inner solver's infeasibility guard).

psi_tv_eps

Positive scalar, Huber smoothing scale for the kinked total-variation conjugate \phi^{*}_{\mathsf{TV}}. The smoothed conjugate lies above the exact one, so computed bounds remain outward-conservative.

tvac_tau

Positive scalar, soft-max temperature for the log-sum-exp term of the "TVac" divergence, that is, total variation restricted to distributions with P \ll P_{*}.

tvmix_tau

Positive scalar, soft-max temperature for the log-sum-exp term of the "TVmix" divergence, that is, total variation intersected with the mixture constraint P \geq \kappa P_{*}.

purekl_acap

Positive scalar, overflow clamp on the exponent of the conjugate of the pure Kullback-Leibler entropy \phi_{\mathsf{KL}}(s) = s \log s - s + 1, which increases exponentially in its argument.

tvmix_kappa

NULL or a scalar in [0,1]: the contamination weight \kappa of the mixture constraint P \geq \kappa P_{*} used by the "TVmix" and "TVmixC" divergences (the perturbed distribution must contain the baseline as a mixing component with weight \kappa). NULL ties \kappa to the budget as \kappa = 1 - \delta, matching Palomba (2026). An explicit value decouples the two parameters; the reduced "TVmix" program requires \kappa \geq 1 - \delta (so that the total-variation constraint is redundant), while "TVmixC" supports any \kappa in [0,1).

Details

The options are stated in the notation of the paper. A candidate distribution P is measured against the baseline P_{*} by the divergence D_{\phi}(P \| P_{*}) generated by an entropy function \phi, and the budget \delta caps it. At a fixed structural parameter \theta the inner problem is solved in its dual form, over the multipliers (\zeta,\eta,\lambda) attached respectively to the total-mass constraint, to the divergence budget – so that \eta is the shadow price of robustness – and to the moment restrictions. The dual objective integrates against the baseline the perspective of the convex conjugate \phi^{*},

(\phi^{*})^{\pi}\big(g(U;\theta) - \lambda^{\top} m(U;\theta) - \zeta, \; \eta\big),

with g the counterfactual criterion and m the moment function. Under total variation, \phi_{\mathsf{TV}}(s) = |s-1|/2, the conjugate \phi^{*}_{\mathsf{TV}} is piecewise linear and the perspective collapses to the kinked \max\{g(U;\theta) - \lambda^{\top} m(U;\theta) - \zeta, \; -\eta/2\}; the smoothing options below round those kinks off, always from above, so the computed bounds stay outward-conservative. The contamination weight \kappa of the mixture constraint P \geq \kappa P_{*} is set by tvmix_kappa.

Value

A list of class "tvbounds_control" with the (validated) options above.

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Christensen, T. and B. Connault (2023). "Counterfactual Sensitivity and Robustness." Econometrica, 91(1), 263-298.

See Also

tvbounds_counterfactual()

Examples

ctrl <- tvbounds_control()
ctrl$maxsolves

# a faster configuration for exploratory grids
tvbounds_control(maxsolves = 3, knitro_options = list(maxit = 200))


Sensitivity bounds for counterfactual predictions in structural models

Description

Computes the lower and upper sensitivity bounds \underline{\mathsf{k}}(\delta) and \overline{\mathsf{k}}(\delta) on a counterfactual \mathbb{E}_P[g(U;\theta)], when the distribution P of the latent variables U ranges over a divergence neighborhood of the simulated baseline P_{*} with budget delta, and the structural parameter \theta ranges over the values compatible with the moment conditions \mathbb{E}_P[m(U;\theta)] \in \mathcal{M}(\rho), following Palomba (2026) and Christensen and Connault (2023). This is the only function in the package that supports general \phi-divergences beyond total variation (spelled out at first use; "TV" below) and contamination.

Usage

tvbounds_counterfactual(
  moments,
  d,
  theta_lb,
  theta_ub,
  delta,
  divergence = c("KL_chi2", "KL", "chi2", "TV", "TVmix", "TVmixC", "TVac"),
  side = c("both", "lower", "upper"),
  U = NULL,
  M = 50000,
  u_dim = NULL,
  gamma = NULL,
  gradient = NULL,
  theta_init = NULL,
  control = tvbounds_control(),
  seed = NULL,
  verbose = TRUE
)

Arguments

moments

The moment/counterfactual function: a length-2 character vector c(file, fname), a single string naming a Julia function, or an R function. See Details.

d

Integer, the number of moment conditions \mathsf{d}_m (columns of G).

theta_lb, theta_ub

Numeric vectors of equal length: the box for the structural parameter \theta over which the outer problems optimize. Their common length fixes the dimension of \theta.

delta

Numeric vector of strictly positive budgets \delta (the sensitivity parameter of the neighborhood), without duplicates. For the total-variation family the budget must lie in (0,1]; for "KL_chi2", "KL", and "chi2" any positive value is allowed.

divergence

Divergence keyword; one of "KL_chi2" (default), "KL", "chi2", "TV", "TVmix", "TVmixC", "TVac". See Details.

side

"both" (default), "lower", or "upper": which bound problems to solve at each budget.

U

Optional ⁠M x u_dim⁠ numeric matrix of latent draws U^{(1)}, \dots, U^{(M)}, one per row. When NULL, scrambled-Halton uniforms are generated (see Details) and M, u_dim are required; when supplied, M and u_dim are taken from its dimensions.

M

Integer, the number M of simulated draws when U = NULL (default 50000, the setting of the paper).

u_dim

Integer, the dimension \mathsf{d}_z of the latent draw when U = NULL (at most 15).

gamma

Optional R list: an arbitrary payload forwarded to the moments function (as obj.gamma for Julia moments, as the third argument for R moments).

gradient

How to differentiate the moments with respect to \theta in the outer optimization: NULL (default; automatic differentiation for Julia moments, finite differences for R moments), the string "fd" (finite differences), the name of a Julia function jac(theta, U, obj), or an R function ⁠function(theta, U, gamma)⁠. A user-supplied gradient must return either the stacked ⁠(M*(d+1)) x l⁠ Jacobian of c(K, G) (K rows first, then G in column-major order) or a list with components K (⁠M x l⁠) and G (⁠M x d x l⁠), where l = length(theta_lb).

theta_init

Optional numeric vector, the initial structural parameter \theta for the outer optimization (defaults to the midpoint of the box). Set it to the baseline estimate of the model: the reported baseline point is the plug-in counterfactual \mathsf{k}(\theta;P_{*}) = \mathbb{E}_{P_{*}}[g(U;\theta)] at theta_init, and is only returned when theta_init is supplied.

control

A list created by tvbounds_control() with solver tuning options.

seed

Optional integer seed for the scrambled-Halton draws (and for R's RNG, which is saved and restored so the call has no side effect on the caller's random-number stream).

verbose

Logical; print solver progress (default TRUE).

Details

Setup. Let U \in \mathcal{Z} collect the latent variables of the structural model (taste shocks, unobserved heterogeneity, measurement errors) and let P_{*} be the baseline distribution the econometrician postulates for them. At a structural parameter \theta \in \Theta, a distribution P is compatible with the model when \mathbb{E}_P[m(U;\theta)] \in \mathcal{M}(\rho), with m the moment function and \mathcal{M}(\rho) the moment constraint set; the object of interest is the counterfactual \mathbb{E}_P[g(U;\theta)], the expectation of a known criterion g that is linear in P at fixed \theta. The robustness set collects the distributions that are compatible with the model and within budget of the baseline,

\mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta) := \{P : D_{\phi}(P \| P_{*}) \leq \delta, \; \mathbb{E}_P[m(U;\theta)] \in \mathcal{M}(\rho)\},

where D_{\phi}(P \| P_{*}) is the \phi-divergence selected by divergence and \delta the budget, and the reported bounds are the nested extrema

\underline{\mathsf{k}}(\delta) = \inf_{\theta \in \Theta} \inf_{P \in \mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta)} \mathbb{E}_P[g(U;\theta)], \qquad \overline{\mathsf{k}}(\delta) = \sup_{\theta \in \Theta} \sup_{P \in \mathcal{P}_{\phi}(\theta;\rho,P_{*},\delta)} \mathbb{E}_P[g(U;\theta)].

The optimization over P at a fixed \theta is the inner problem, solved in its dual form; the optimization over \theta is the outer problem.

Correspondence between the paper and the code. The paper writes the counterfactual integrand as g and the moment function as m. The interfaces below fill an array named K with the values of g and an array named G with the values of m: read K as g and G as m throughout. The solver takes the moment conditions in centered equality form, \mathbb{E}_P[m(U;\theta)] = 0, so a nonzero target \rho is absorbed by centering the moment function. The argument d is the number of moment conditions, \mathsf{d}_m; the common length of theta_lb and theta_ub is the dimension of \theta; and u_dim is the dimension \mathsf{d}_z of a single latent draw.

Moments specification. moments can be supplied in three forms:

  1. a length-2 character vector c(file, fname): file is the path of a Julia source file that is included into the session, and fname the name of a function defined there (at the top level, in Main) with the in-place signature ⁠moments!(K, G, theta, U, obj)⁠. The function must fill K (an M-vector holding the counterfactual values g(U^{(j)};\theta)) and G (an ⁠M x d⁠ matrix holding the moment functions m(U^{(j)};\theta), one row per draw) and may read the user payload as obj.gamma (an R list arrives in Julia as an ordered dictionary keyed by symbols, so ⁠obj.gamma[:name]⁠);

  2. a single string naming a Julia function with the same signature that is already defined in the session;

  3. an R function ⁠function(theta, U, gamma)⁠ returning a list with components K (numeric of length M) and G (numeric ⁠M x d⁠ matrix). This path is much slower (every objective evaluation crosses the R/Julia boundary), and because ForwardDiff cannot differentiate through R code the outer optimization needs either a user-supplied gradient or finite differences (the default for this path).

For Julia moments the outer envelope-theorem gradient differentiates the moments by automatic differentiation (ForwardDiff), so the Julia function should be written generically in the element type of theta; pass gradient = "fd" for a non-generic function.

Latent draws. U is the ⁠M x u_dim⁠ matrix of latent draws U^{(1)}, \dots, U^{(M)} that discretizes the baseline distribution P_{*}, row j holding the draw U^{(j)}. When U = NULL the package generates M scrambled-Halton points in the unit cube (0,1)^{\mathsf{d}_z} of dimension u_dim (Owen, 2017), seeded by seed; the moments function is then responsible for mapping the uniform coordinates into baseline draws (e.g. through quantile transforms). The Halton generator supports u_dim <= 15; supply U directly for higher-dimensional draws.

Divergences. divergence selects the entropy function \phi whose divergence D_{\phi}(P \| P_{*}) defines the neighborhood; the budget grid delta must be strictly positive, and must lie in (0,1] for the total-variation family:

The kinked total-variation conjugates \phi^{*}_{\mathsf{TV}} are Huber-smoothed and the per-draw maxima in the dual objective log-sum-exp-smoothed (scales in tvbounds_control()); both smoothings lie above the exact functions, so computed bounds are outward-conservative (wider, never narrower) at order 1e-3.

Optimization. For each budget (in increasing order) and each side, the solver runs control$maxsolves multi-start outer optimizations over \theta (KNITRO, or Optim.jl when control$use_optim = TRUE), warm-started at the previous budget's optimum; the reported bound is the inner (dual) value re-solved at the best candidate, which makes the bound curves monotone in the budget by construction. A degenerate box (theta_lb == theta_ub) skips the outer optimization and reports the bounds at the fixed \theta supplied through theta_init. Failed budgets are reported as NA (for "TVmix" an NA typically signals an infeasible moment condition at every \theta in the box, that is an empty robustness set at that budget).

Value

An object of class c("tvbounds_counterfactual", "tvbounds"): a list with the fields described in the package overview, in particular bounds (data frame with columns delta, lower, upper, holding \delta, \underline{\mathsf{k}}(\delta) and \overline{\mathsf{k}}(\delta); no standard-error or confidence-band columns, since this application currently carries no inference), point (the plug-in counterfactual \mathbb{E}_{P_{*}}[g(U;\theta)] at theta_init, or NA when theta_init was not supplied), divergence, and details, a list with:

solver

data frame of per-budget diagnostics: outer multi-start flags, inner KNITRO status codes, and timings for each side (-999 marks entries that do not apply, e.g. outer flags in fixed-theta mode).

theta_lower, theta_upper

⁠l x length(delta)⁠ matrices of the outer-optimal structural parameters \theta at each budget, for the lower and the upper bound respectively.

M, u_dim

the number M of simulated draws and their dimension \mathsf{d}_z.

theta_lb, theta_ub, theta_init, fixed_theta

the parameter box, the initial point, and whether the box was degenerate.

control

the resolved control list, including the option files actually used.

moments

a short description of the moments specification.

KNITRO requirement

This function relies on Julia (>= 1.9) and on the commercial Artelys KNITRO solver, accessed through the 'JuliaCall' package and KNITRO.jl. A valid KNITRO license is required (free academic trials are available from Artelys at https://www.artelys.com/solvers/knitro/). On the first call in each R session the package initializes Julia, instantiates its Julia environment, and checks that KNITRO.jl loads and that a KNITRO solver context can be created (which exercises the license); a one-time message reports the outcome, and the call stops with installation and license guidance when the check fails. The check runs once per R session.

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Christensen, T. and B. Connault (2023). "Counterfactual Sensitivity and Robustness." Econometrica, 91(1), 263-298.

Owen, A. B. (2017). "A randomized Halton algorithm in R." arXiv:1706.02808.

See Also

tvbounds_control(), tvbounds_plot(), tvbounds_summary()

Examples

# The full solver requires Julia and a licensed KNITRO installation,
# so a complete run cannot be executed on CRAN or in checks:
## Not run: 
# Toy model (shipped with the package): U ~ Uniform(0, 1), one moment
# condition m(U; theta) = U - theta and counterfactual g(U; theta) = U,
# so under the "TVmix" neighborhood the bounds equal the endpoints of
# the theta box. In the Julia file, m is written into G and g into K.
toy <- system.file("julia", "examples", "toy.jl", package = "tvbounds")
fit <- tvbounds_counterfactual(
  moments  = c(toy, "tvb_toy_moments!"),
  d        = 1,
  theta_lb = 0.4, theta_ub = 0.6,
  delta    = c(0.5, 1),
  divergence = "TVmix",
  M = 500, u_dim = 1,
  theta_init = 0.5,
  control = tvbounds_control(maxsolves = 2),
  seed = 1234)
fit$bounds

## End(Not run)

# The control constructor is pure R and always available:
tvbounds_control(maxsolves = 3)$maxsolves


Sensitivity bounds for recentered instrumental variables

Description

Computes sensitivity bounds on a recentered (formula) instrumental-variables estimate when the postulated distribution of the shocks is allowed to vary within a total variation or a contamination neighborhood of the baseline assignment distribution P_{*}, as in Palomba (2026). The leading use case is the recentered instruments of Borusyak and Hull (2023), whose validity rests on a researcher-postulated distribution for the shock process: the bounds quantify how far the estimate can move when up to a fraction delta of that postulated probability mass is misspecified.

Usage

tvbounds_riv(
  y,
  x,
  z,
  Fmat,
  p = NULL,
  controls = NULL,
  delta = seq(0, 1, by = 0.002),
  neighborhood = c("tv", "contamination"),
  tau_star = 0,
  verbose = FALSE
)

Arguments

y

Numeric n-vector, the outcome y_i.

x

Numeric n-vector, the endogenous regressor x_i.

z

Numeric n-vector, the realized (un-recentered) candidate instrument: z[i] is the formula of unit i evaluated at the realized shocks, z_i = f_i(v;w).

Fmat

Numeric n x S matrix of counterfactual instrument draws: Fmat[i, s] is the formula of unit i evaluated at the s-th counterfactual shock configuration, f_i(v^{(s)};w). These are the draws from the postulated assignment process used for recentering (in Borusyak and Hull (2023), permuted or re-simulated shock allocations).

p

Optional numeric S-vector holding the probabilities that the postulated assignment distribution P_{*} attaches to the columns of Fmat, that is to the configurations v^{(1)}, \dots, v^{(S)}. Defaults to NULL, meaning the uniform distribution 1/S on each draw, which is the postulated assignment distribution of Borusyak and Hull (2023). Must be nonnegative and sum to one.

controls

Optional numeric matrix or data frame of control variables (n rows) to be partialled out of y, x, z, and each column of Fmat. A constant is always included. Default NULL.

delta

Numeric vector of sensitivity budgets \delta \in [0,1] at which the bounds are traced. Default seq(0, 1, by = 0.002).

neighborhood

Either "tv" (the total variation ball \mathcal{P}^{\mathsf{FI}}_{\mathsf{TV}}(\delta), the default) or "contamination" (the contamination neighborhood \mathcal{P}^{\mathsf{FI}}_{\mathsf{cont}}(\delta)).

tau_star

Numeric scalar reference value \tau_{\star} for the breakdown budget \delta_b(\tau_{\star}): the smallest budget at which the bounds cover \tau_{\star}. The default 0 gives the breakdown budget for the sign of \beta.

verbose

Logical; if TRUE, print progress messages. Default FALSE.

Details

The exercise is conducted conditionally on the realized sample, so every bound is a deterministic function of the data and of the budget delta; accordingly, and by design (as in the paper), no standard errors or confidence bands are produced for this application.

The model and the recentered estimate. For units i \in [n] the structural equation is

y_i = \beta x_i + \varepsilon_i,

with \beta the parameter of interest, y_i the outcome, x_i the endogenous regressor and \varepsilon_i the unobserved residual. Let v denote the vector of exogenous shocks, taking values in a space \mathcal{V}, let w collect predetermined covariates, and let f_i(\cdot\,;w) : \mathcal{V} \to \mathbb{R} be the known formula that maps a shock configuration into the instrument of unit i, so that z_i := f_i(v;w) is the candidate instrument at the realized shocks. For a distribution P on \mathcal{V}, the expected instrument and the recentered instrument are

\mu_i(P) := \mathbb{E}_P[f_i(v;w) \mid w] = \int_{\mathcal{V}} f_i(v';w) \,\mathrm{d}P(v'), \qquad \widetilde{z}_i(P) := z_i - \mu_i(P).

Borusyak and Hull (2023) postulate an assignment distribution P_{*} for the shocks and recenter at it.

The two criterion functions. The formula enters only through the two sample aggregates

g_y(\cdot) := \sum_{i=1}^{n} y_i f_i(\cdot\,;w), \qquad g_x(\cdot) := \sum_{i=1}^{n} x_i f_i(\cdot\,;w),

whose recentered values are the reduced form G_y(P) := g_y(v) - \mathbb{E}_P[g_y] and the first stage G_x(P) := g_x(v) - \mathbb{E}_P[g_x]. The estimate at a candidate assignment distribution is the ratio

\widehat{\beta}(P) := \frac{G_y(P)}{G_x(P)} = \frac{\sum_{i=1}^{n} \widetilde{z}_i(P) \, y_i}{ \sum_{i=1}^{n} \widetilde{z}_i(P) \, x_i},

and the reported estimate is \widehat{\beta}_{*} = \widehat{\beta}(P_{*}).

How the arguments encode the shock space. The package represents P_{*} by S shock configurations v^{(1)}, \dots, v^{(S)}: entry Fmat[i, s] holds f_i(v^{(s)};w), the formula of unit i at the s-th configuration, p[s] holds the probability that P_{*} assigns to that configuration, and z[i] holds z_i = f_i(v;w). The shock space is therefore taken to be the finite set \mathcal{V} = \{v^{(1)}, \dots, v^{(S)}\}, over which the candidate distributions P range.

Neighborhoods. The bounds report the range of \widehat{\beta}(P) as P ranges over the chosen neighborhood of P_{*}, intersected with the set \mathcal{P}_{\neq 0} := \{P \in \Delta(\mathcal{V}) : G_x(P) \neq 0\} of distributions at which the ratio is defined:

Every member of \mathcal{P}^{\mathsf{FI}}_{\mathsf{cont}}(\delta) lies in \mathcal{P}^{\mathsf{FI}}_{\mathsf{TV}}(\delta), so the contamination bounds are weakly tighter at every budget.

First-stage breakdown. Once the budget is large enough that some distribution in the neighborhood makes the recentered first stage G_x(P) vanish, \widehat{\beta}(P) is no longer well defined over the whole neighborhood and the identified set is the entire real line. The smallest such budget is the first-stage breakdown budget, \delta^{\mathsf{TV}}_{\mathsf{FS}} for the total variation ball and \delta^{\mathsf{cont}}_{\mathsf{FS}} for the contamination neighborhood; since the contamination neighborhood is the smaller of the two, \delta^{\mathsf{TV}}_{\mathsf{FS}} \leq \delta^{\mathsf{cont}}_{\mathsf{FS}}. The one for the neighborhood in use is reported in details$delta_fs. Following the paper, the infimum over an empty set is set to 1, so a reported value of 1 carrying delta_fs_censored = TRUE means that the first stage never breaks down over the budget range, not that breakdown occurs at 1. Rows of bounds at budgets where the bounds are vacuous carry NA.

Controls. When controls is supplied, y, x, z, and every column of Fmat are residualized on the controls and a constant (when controls = NULL, on the constant alone) before the bounds are computed. By the Frisch–Waugh–Lovell theorem this leaves the just-identified two-stage least squares coefficient on x unchanged, because the projection matrix is idempotent; the baseline estimate then replicates the estimate from the full regression with controls.

Value

An object of class c("tvbounds_riv", "tvbounds"), a list with entries

References

Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working paper.

Borusyak, K. and Hull, P. (2023). "Nonrandom Exposure to Exogenous Shocks." Econometrica, 91(6), 2155–2185.

Examples

# A small simulated formula-instrument design: S counterfactual shock
# configurations, a realized instrument, a first stage, and an outcome.
set.seed(123)
n <- 80
S <- 50
Fmat <- matrix(rnorm(n * S), n, S)     # Fmat[i, s] = f_i(v^(s); w)
z <- rowMeans(Fmat) + rnorm(n)         # realized instrument z_i = f_i(v; w)
x <- z + 0.5 * rnorm(n)                # endogenous regressor
y <- 0.4 * x + rnorm(n)                # outcome

fit <- tvbounds_riv(y, x, z, Fmat, delta = seq(0, 1, by = 0.05))
fit$point                              # recentered IV estimate at P_*
head(fit$bounds)                       # bounds along the budget grid
fit$details$delta_fs                   # first-stage breakdown budget
fit$details$delta_breakdown            # breakdown budget for the sign

# Contamination neighborhood, with controls partialled out.
W <- data.frame(w1 = rnorm(n), w2 = rnorm(n))
fit_cont <- tvbounds_riv(y, x, z, Fmat, controls = W,
                         delta = seq(0, 1, by = 0.05),
                         neighborhood = "contamination")
fit_cont$point