Uniform inference for value functions

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

97,110 characters · 14 sections · 100 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\thispagestyle{empty}

abstract\begin{spacing}{1} We propose a method to conduct uniform inference for the (optimal) value function, that is, the function that results from optimizing an objective function marginally over one of its arguments. Marginal optimization is not Hadamard differentiable (that is, compactly differentiable) as a map between the spaces of objective and value functions, which is problematic because standard inference methods for nonlinear maps usually rely on Hadamard differentiability. However, we show that the map from objective function to an $L_p$ functional of a value function, for $1 \leq p \leq \infty$, are Hadamard directionally differentiable. As a result, we establish consistency and weak convergence of nonparametric plug-in estimates of Cram\'er-von Mises and Kolmogorov-Smirnov test statistics applied to value functions. For practical inference, we develop detailed resampling techniques that combine a bootstrap procedure with estimates of the directional derivatives. In addition, we establish local size control of tests which use the resampling procedure. Monte Carlo simulations assess the finite-sample properties of the proposed methods and show accurate empirical size and nontrivial power of the procedures. Finally, we apply our methods to the evaluation of a job training program using bounds for the distribution function of treatment effects. \end{spacing}

Keywords: uniform inference, value function, bootstrap, delta method, directional differentiability

JEL codes: C12, C14, C21

Introduction

\setcounter{page}{1}

The (optimal) value function, that is, a function defined by optimizing an objective function over one of multiple arguments, is used widely in economics, finance, statistics and other fields. The statistical behavior of estimated value functions appears, for instance, in the study of distribution and quantile functions of treatment effects under partial identification (see, e.g., FirpoRidder08, FirpoRidder19 and FanPark10,FanPark12). More recently, investigating counterfactual sensitivity, ChristensenConnault19 construct interval-identified counterfactual outcomes as the distribution of unobservables varies over a nonparametric neighborhood of an assumed model specification. Statistical inference for value functions can be analyzed at a point as in Roemisch06 or CarcamoCuevasRodriguez19 and inference can be carried out using the techniques of HongLi18 or FangSantos19.

This paper extends these pointwise inference methods to uniform inference for value functions that are estimated with nonparametric plug-in estimators. Uniform inference has an important role for the analysis of statistical models.\footnote{For example, AngristChernozhukovVal06 discuss the importance of uniform inference vis-\`a-vis pointwise inference in a quantile regression framework. Uniform testing methods allow for multiple comparisons across the entire distribution without compromising confidence levels, and uniform confidence bands differ from corresponding pointwise confidence intervals because they guarantee coverage over the family of confidence intervals simultaneously.} However, uniform inference for the value function can be quite challenging. In order to use standard statistical inference procedures, one must show that the map from objective function to value function is Hadamard (i.e., compactly) differentiable as a map between two functional spaces, which may require restrictive regularity conditions. As a result, the delta method, as the conventional tool used to analyze the distribution of nonlinear transformations of sample data, may not apply to marginal optimization as a map between function spaces due to a lack of differentiability. We propose a solution to this problem that allows for the use of the delta method to conduct uniform inference for value functions.

Because in general the map of objective function to value function is not differentiable, our solution is to directly analyze statistics that are applied to the value function. This is in place of a more conventional analysis that would first establish that the value function is well-behaved and as a second step use the continuous mapping theorem to determine the distribution of test statistics. In settings involving a chain of maps, this can be seen as choosing how to define the “links” in the chain to which the chain rule applies. The results are presented for $L_p$ functionals for $1 \leq p \leq \infty$, that are applied to a value function, which should cover many cases of interest. In particular, this family of functionals includes Kolmogorov-Smirnov and Cram\'er-von Mises tests.

By considering $L_p$ functionals of the value function, we may bypass the most serious impediment to uniform inference. However, these $L_p$ functionals are only Hadamard directionally differentiable. A directional generalization of Hadamard (or compact) differentiability was developed for statistical applications in Shapiro90 and Duembgen93. There is a more recent growing literature of applications using Hadamard directional differentiability. See, for example, SommerfeldMunk17, ChetverikovSantosShaikh18, ChoWhite18, HongLi18, FangSantos19, ChristensenConnault19, MastenPoirier20, and CattaneoJanssonNagasawa20. We separately analyze $L_p$ functionals for $1 \leq p < \infty$ and the supremum norm, as well as one-sided variants used for tests that have a sense of direction. As an intermediate step towards showing the differentiability of supremum-norm statistics, we show directional differentiability of a minimax map applied to functions that are not necessarily convex-concave.

Directional differentiability of $L_p$ functionals provides the minimal amount of regularity needed for feasible statistical inference. We use directional differentiability to find the asymptotic distributions of test statistics estimated from sample data. The distributions generally depend on features of the objective function, and we show how to combine estimation with resampling to estimate the limiting distributions, as was proposed by FangSantos19. This technique can be tailored to impose a null hypothesis and is simple to implement in practice. We provide three examples throughout the paper to illustrate the proposed methods --- dependency bounds, stochastic dominance, and cost and profit functions. Moreover, we establish local size control of the resampling procedure in general and uniform size control for some functionals.

As a practical illustration of the developed framework, we consider bounds for the distribution function of a treatment effect. We use this model in both the Monte Carlo simulations and empirical application. Treatment effects models have provided a valuable analytic framework for the evaluation of causal effects in program evaluation studies. When the treatment effect is heterogeneous, its distribution function and related features, for example its quantile function, are often of interest. Nevertheless, it is well known (e.g. HeckmanSmithClements97, AbbringHeckman07) that features of the distribution of treatment effects beyond the average are not point identified unless one imposes significant theoretical restrictions on the joint distribution of potential outcomes. Therefore it has became common to use partial identification of the distribution of such features, and to derive bounds on their distribution without making any unwarranted assumptions about the copula that connects the marginal outcome distributions.\footnote{Bounds for the CDF and quantile functions at one point in the support of the treatment effects distribution were developed in Makarov82, Rueschendorf82, FrankNelsenSchweizer87, WilliamsonDowns90. These bounds are minimal in the sense that they rely only on the marginal distribution functions of the observed outcomes, not the joint distribution of (potential) outcomes.}

Monte Carlo simulations (reported in a supplementary appendix) assess the finite-sample properties of the proposed methods. The simulations suggest that Kolmogorov-Smirnov statistics used to construct uniform confidence bands have accurate empirical size and power against local alternatives. A Cram\'er-von Mises type statistic is used to test a stochastic dominance hypothesis using bound functions, similar to the breakdown frontier tests of MastenPoirier20. This test has accurate size at the breakdown frontier and power against alternatives that violate the dominance relationship. All results are improved when the sample size increases and are powerful using a modest number of bootstrap repetitions.

We illustrate the methods empirically using job training program data set first analyzed by LaLonde86 and subsequently by many others, including HeckmanHotz89, DehejiaWahba99, SmithTodd01, SmithTodd05, Imbens03, and Firpo07, without making any assumptions on dependence between potential outcomes in the experiment. This experimental data set has information from the National Supported Work Program used by DehejiaWahba99. We document strong heterogeneity of the job training treatment across the distribution of earnings. The uniform confidence bands for the treatment effect distribution function are imprecisely estimated at some parts of the earnings distribution. These large uniform confidence bands may in part be attributed to the large number of zeros contained in the data set, but are also inherent to the fact that the distribution function is everywhere only partially identified.

The remainder of the paper is organized as follows. Section (ref) defines the statistical model of interest for the value function. Section (ref) establishes directional differentiability for the functionals of interest. Inference procedures are established in Section (ref). An empirical application to job training is discussed in Section (ref). Finally, Section (ref) concludes. Monte Carlo simulation descriptions and results and all proofs are relegated to the supplementary appendix.

Notation

For any set $T \subseteq \mathbb{R}^d$ let $\ell^\infty(T)$ denote the space of bounded functions $f: T \rightarrow \mathbb{R}$ and let $\mathcal{C}(T)$ denote the space of continuous functions $f: T \rightarrow \mathbb{R}$, both equipped with the uniform norm $\| f \|_\infty = \sup_{t \in T} |f(t)|$. Given a sequence $\{f_n\}_n \subset \ell^\infty(T)$ and limiting random element $f$ we write $f_n \leadsto f$ to denote weak convergence in $(\ell^\infty(T), \| \cdot \|_\infty)$ in the sense of Hoffmann-J\o rgensen vanderVaartWellner96. Let $[x]_+ = \max\{x, 0\}$ and $[x]_- = \min\{x, 0\}$. A set-valued map (or correspondence) $S$ that maps elements of $X$ to the collection of subsets of $Y$ is denoted $S: X \rightrightarrows Y$. For set-valued map $S: X \rightrightarrows Y$, let $\textnormal{gr} S$ denote the graph of $S$ as a set in $X \times Y$.

The model

Consider two sets $U \subseteq \mathbb{R}^{d_U}$ and $X \subseteq \mathbb{R}^{d_X}$, a set-valued map $A: X \rightrightarrows U$ that serves as a choice set, and an objective function $f \in \ell^\infty(\textnormal{gr} A)$. The objective function and the value function are linked by a marginal optimization step that maps one functional space to another. Let $\psi: \ell^\infty(\textnormal{gr} A) \rightarrow \ell^\infty(X)$ map the objective function to the value function obtained by marginally optimizing the objective $f$ with respect to $u \in A(x)$ for each $x \in X$. Without loss of generality, we consider only marginal maximization with respect to $u$ for each $x$:

equation[equation omitted — 71 chars of source]

The value function $\psi(f)$ is the object of interest and we would like to conduct uniform inference using a plug-in estimator $\psi(f_n)$. First we introduce examples of the value function model represented in equation (ref).

exm[Dependency bounds] Dependency bounds describe bounds on the distribution function of the sum of two random variables $X$ and $Y$, $F_{X+Y}$, using only the marginal distribution functions $F_X$ and $F_Y$ and no information from the joint distribution $F_{XY}$. FrankNelsenSchweizer87 showed that bounds derived for $F_{X+Y}$ by Makarov82 are related to copula bounds and can be extended to other binary operations. Upper and lower dependency bounds for the CDFs $F_{X + Y}$, $F_{X - Y}$, $F_{X \cdot Y}$ and $F_{X / Y}$ and their associated quantile functions are provided in WilliamsonDowns90. WilliamsonDowns90 show that for any $x$, $L_{X+Y}(x) \leq F_{X+Y}(x) \leq U_{X+Y}(x)$, where \begin{align} L_{X+Y}(x) &= \left[ \sup_{u \in \mathbb{R}} (F_X(u) + F_Y(x - u) - 1) \right]_+ \\ U_{X+Y}(x) &= 1 - \left[ \inf_{u \in \mathbb{R}} (F_X(u) + F_Y(x - u)) \right]_-. \end{align} Similarly, letting $f^{-1}$ denote the generalized inverse of an increasing function, bounds for the quantile function $F_{X+Y}^{-1}(\tau)$ may be considered for $0 < \tau < 1$. The functions above may be inverted to find that for each quantile level $\tau$, $U_{X+Y}^{-1}(\tau) \leq F_{X+Y}^{-1}(\tau) \leq L_{X+Y}^{-1}(\tau)$, and the quantile function bounds are \begin{align} U_{X+Y}^{-1}(\tau) &= \sup_{u \in (0, \tau)} \left\{ F_X^{-1}(u) + F_Y^{-1}(\tau - u) \right\} \\ L_{X+Y}^{-1}(\tau) &= \inf_{u \in (\tau, 1)} \left\{ F_X^{-1}(u) + F_Y^{-1}(\tau - u + 1) \right\}. \end{align} The functions $U_{X+Y}^{-1}$ and $L_{X+Y}^{-1}$ are value functions where marginal optimization takes place over a set-valued map that varies with $\tau$. All these related functions may be written as maps that depend on a value function in some part of a chain of maps. For example, take $L_{X+Y}$. For a pair of functions $f = (f_x, f_y) \in (\ell^\infty(\mathbb{R}))^2$, let $\Pi^+: (\ell^\infty(\mathbb{R}))^2 \rightarrow \ell^\infty(\mathbb{R}^2)$ \begin{equation} \Pi^+(f)(s, t) = f_x(s) + f_y(t - s) - 1 \end{equation} and let $\psi$ be as in (ref) with $A(x) = \mathbb{R}$ for all $x$. Then $L_{X+Y} = \psi(\Pi^+(F) \vee 0)$. Similarly, consider $U_{X+Y}^{-1}$. For a pair of CDFs $f = (f_x, f_y)$ where $f_x, f_y: \mathbb{R} \rightarrow [0, 1]$, define $\tilde{\Pi}^+: (\ell^\infty(\mathbb{R}))^2 \rightarrow \ell^\infty([0, 1]^2)$ by\footnote{We have simplified the problem by excluding $\{0, 1\}$ from the definition of $\tilde{\psi}$ for ease of exposition --- see p.115-117 of WilliamsonDowns90 for more complex cases which include these points.} \begin{equation} \tilde{\Pi}^+(f)(s, t) = f_x^{-1}(s) + f_y^{-1}(t - s) \end{equation} and let $\psi$ use $A(\tau) = (0, \tau)$ for each $\tau \in (0, 1)$. Then $U_{X+Y}^{-1} = \psi(\tilde{\Pi}^+(F))$. We will focus on constructing uniform confidence bands for $L_{X+Y}$ and $U_{X+Y}^{-1}$, since the calculations for the other functions are analogous.

exm[Stochastic dominance] Inference for features of $F_{X-Y}$ at a given $x$, when $X - Y$ may be considered a treatment effect distribution, has been considered in FanGuerreZhu17 and the references cited therein. Consider testing for first order stochastic dominance without assuming point identification. Suppose that $(X_A, X_B, X_0)$ have marginal CDFs $(G_A, G_B, G_0) \in (\ell^\infty(\mathbb{R}))^3$, and let $F_{\Delta_A}$ and $F_{\Delta_B}$ be the distribution functions of the treatment effects $\Delta_A = X_A - X_0$ and $\Delta_B = X_B - X_0$. Without point identification we cannot test the hypothesis $H_0: F_{\Delta_A} \succeq_{FOSD} F_{\Delta_B}$, which is equivalent to the condition that $F_{\Delta_A}(x) - F_{\Delta_B}(x) \leq 0$ for all $x$. However, the distributions of $\Delta_A$ and $\Delta_B$ may be bounded. For example, using references from the previous example it can be shown that for each $x$, \begin{align*} F_{\Delta_A}(x) \geq L_A(x) &= \left[ \sup_{u \in \mathbb{R}} (G_A(u) - G_0(u - x)) \right]_+, \\ F_{\Delta_B}(x) \leq U_B(x) &= 1 + \left[ \inf_{u \in \mathbb{R}} (G_B(u) - G_0(u - x)) \right]_-. \end{align*} Under the null hypothesis, $L_A \leq F_{\Delta_A} \leq F_{\Delta_B} \leq U_B$, so that $L_A(x) - U_B(x) \leq F_{\Delta_A} - F_{\Delta_B} \leq 0$ for all $x$, and the null hypothesis may be rejected whenever there exists an $x'$ such that $L_A(x') > U_B(x')$. This holds regardless of the correlations between treatment outcomes. Let $\Pi^-: (\ell^\infty(\mathbb{R}))^2 \rightarrow \ell^\infty(\mathbb{R}^2)$ be, for $f = (f_0, f_1)$ \begin{equation*} \Pi^-(f)(u, x) = f_1(u) - f_0(u - x) \end{equation*} and let $\psi$ be defined by (ref) with $A(x) = \mathbb{R}$ for all $x$. The function used to indicate a violation of dominance maps three distribution functions into a function in $\ell^\infty(\mathbb{R})$: for each $x \in \mathbb{R}$, \begin{align} L_A(x) - U_B(x) &= \sup_{u \in \mathbb{R}} (G_A(u) - G_0(u - x)) - 1 - \inf_{u \in \mathbb{R}} (G_B(u) - G_0(u - x)) \\ &= \psi(\Pi^-(G_0, G_A))(x) - 1 + \psi(-\Pi^-(G_0, G_B))(x). \end{align} We would like to test the hypothesis that $F_A$ dominates $F_B$, without knowledge of the dependence between $X_0$, $X_A$ and $X_B$, by looking for $x$ where $L_A(x) - U_B(x) > 0$, a case in which a one-sided functional $f \mapsto \|[f]_+\|_p$ is preferred for testing.

exm[Cost and profit functions] Suppose a firm produces outputs $x \in \mathbb{R}^n_+$ for $n \geq 1$. Given cost function $c: \mathbb{R}^n \rightarrow \mathbb{R}$, the firm's profit at given product prices $p \in \mathbb{R}^n$ is $\pi(p) = \sup_{x \in \mathbb{R}^n_+} (\langle p, x \rangle - c(x))$. Suppose that we have an accurate model of costs and would like to conduct inference for the profit curve, for example by putting a confidence band around an estimated profit curve. The Legendre-Fenchel transform can be used to find the dual conjugate function of a convex function $f$ or the convex hull of a non-convex function $f$ RockafellarWets98. For convex $f: \mathbb{R}^n \rightarrow \mathbb{R}$, the conjugate $f^\ast: \mathbb{R}^n \rightarrow \mathbb{R}$ is $f^\ast(x) = \sup_{u \in \mathbb{R}^n} \left( \langle u, x \rangle - f(u) \right)$, that is, the smallest function such that $f^\ast(x) + f(u) \geq \langle u, x \rangle$. We can recast this as an optimization map by defining the Legendre transform as a map $\mathcal{L}: \ell^\infty(\mathbb{R}^n) \rightarrow \ell^\infty(\mathbb{R}^n)$ defined by \begin{equation*} \mathcal{L}(f)(x) = \sup_{u \in \mathbb{R}^n} \left( \langle u, x \rangle - f(u) \right). \end{equation*} Then the profit function can be seen as the conjugate dual function of the cost function: $\pi(p) = \mathcal{L}(c)(p) = \sup_{x \in \mathbb{R}^n} \left( \langle p, x \rangle - f(x) \right)$. As discussed in the text, it is possible (i.e., when it is not possible to directly calculate $\pi$ from the cost function) that $\pi$ is not Hadamard differentiable. Instead, consider the supremum of the absolute value of transformed functions $\check{\lambda}(f) = \|\mathcal{L}(f)\|_\infty$. Since $\pi = \mathcal{L}(c)$, we can find a confidence band for $\pi$ using $\check{\lambda}(c)$ using a sample analog $\check{\lambda}(c_n)$.

In Example (ref), the bounds on the quantile function are an example where $A(x)$ varies with $x$. It would be notationally simpler to consider optimization problems over the rectangle $\mathcal{X} \times \mathcal{U}$. However, given an $x$-varying set of constraints we would need to use “indicator functions” as they are used in the optimization literature, e.g. we would change the optimization problem to $\tilde{f}(u, x) = f(u, x) + \delta_A(u, x)$, where $\delta_A(u, x) = 0$ for $u \in A(x)$ and (for maximization) $\delta_A(u, x) = -\infty$ for $u \notin A(x)$. This takes us out of the realm of bounded functions and towards the consideration of epi- or hypographs, which may be more complicated. For example, see BuecherSegersVolgushev14. The literature on weak convergence in the space of bounded functions is voluminous and more familiar to researchers in economics.

To analyze the asymptotic properties of the above examples with plug-in estimates, one would typically rely on the delta method because $\psi$ is a nonlinear map from objective function to value function. In the next section we discuss the difficulties encountered with the delta method for value functions and a solution for the purposes of uniform inference.

Differentiability properties

We use the delta method to establish uniform inference methods for value functions. This depends on the notion of Hadamard (or compact) differentiability to linearize the map near the population distribution (see, for example, vanderVaartWellner96), and has the advantage of dividing the analysis into a deterministic part and a statistical part. We first consider Hadamard derivatives without regard to sample data, before explicitly considering the behavior of sample statistics in the next section.

Directional differentiability

The delta method between two metric spaces usually relies on (full) Hadamard differentiability vanderVaartWellner96. Shapiro90, Duembgen93 and FangSantos19 discuss Hadamard directional differentiability and show that this weaker notion also allows for the application of the delta method.\footnote{The derivative $\phi_f'$ in the definition is also called a semiderivative RockafellarWets98.}

defn[Hadamard directional differentiability] Let $\mathbb{D}$ and $\mathbb{E}$ be Banach spaces and consider a map $\phi: \mathbb{D}_\phi \subseteq \mathbb{D} \rightarrow \mathbb{E}$. $\phi$ is Hadamard directionally differentiable at $f \in \mathbb{D}_\phi$ tangentially to a set $\mathbb{D}_0 \subseteq \mathbb{D}$ if there is a continuous map $\phi'_f: \mathbb{D}_0 \rightarrow \mathbb{E}$ such that \begin{equation*} \lim_{n \rightarrow \infty} \left\| \frac{\phi(f + t_n h_n) - \phi(f)}{t_n} - \phi_f'(h) \right\|_\mathbb{E} = 0 \end{equation*} for all sequences $\{h_n\} \subset \mathbb{D}$ and $\{t_n\} \subset \mathbb{R}_+$ such that $h_n \rightarrow h \in \mathbb{D}_0$ and $t_n \rightarrow 0^+$ as $n \rightarrow \infty$ and $f + t_n h_n \in \mathbb{D}_\phi$ for all $n$.

In case the derivative map $\phi'_f$ is linear and $t_n$ may approach $0$ from both negative and positive sides, the map is fully Hadamard differentiable. The delta method and chain rule can be applied to maps that are either Hadamard differentiable or Hadamard directionally differentiable.

To discuss the differentiability of $\psi$ in equation (ref), for any $\epsilon \geq 0$ let $U_f: X \rightrightarrows U$ define the set-valued map of $\epsilon$-maximizers of $f(\cdot, x)$ in $u$ for each $x \in X$:

equation[equation omitted — 129 chars of source]

Results on the directional derivatives of the optimization map $f \mapsto \sup_X f$ date to Danskin67. For example, CarcamoCuevasRodriguez19 show that, in our notation, for fixed $x \in X$, $\psi$ is directionally differentiable in $\ell^\infty(A(x) \times \{x\})$, and for directions $h(\cdot, x)$

equation[equation omitted — 125 chars of source]

Assuming the stronger conditions that $A$ is continuous and compact-valued and that $f$ is continuous on the graph of $A$ ($\textnormal{gr} A$), the maximum theorem implies that $U_f$ is non-empty, compact-valued and upper hemicontinuous AliprantisBorder06, and CarcamoCuevasRodriguez19 show that tangentially to $\mathcal{C}(U \times \{x\})$, the derivative of $\psi(f)(x)$ simplifies to $\psi'_f(h)(x) = \sup_{u \in U_f(x, 0)} h(u, x)$. Further, when (for fixed $x$) $U_f(x, 0)$ is a singleton set, $\psi(f)(x)$ is fully differentiable, not just directionally so.

The pointwise directional differentiability of $\psi(f)(x)$ at each $x$ might lead one to suspect that $\psi$ is differentiable more generally as a map from $\ell^\infty(\textnormal{gr} A)$ to $\ell^\infty(X)$. However, this is not true. The following example illustrates a case where $\psi$ is not Hadamard directionally differentiable as a map from $\ell^\infty(\textnormal{gr} A)$ to $\ell^\infty(X)$, although $f$ and $h$ are both continuous functions.

exmLet $U \times X = [0, 1] \times [-1, 1]$ and define $f: U \times X \rightarrow \mathbb{R}$ by \begin{equation*} f(u, x) = \begin{cases} 0 & x \in [-1, 0] \\ x(u^2 - u) & x \in (0, 1] \end{cases}. \end{equation*} Let $h(u, x) = -u^2 + u$. Then $\psi(f)(x) = 0$ for all $x$ and for any $t > 0$, \begin{equation*} \psi(f + th)(x) = \begin{cases} t / 4 & x \in [-1, 0] \\ (t - x) / 4 & x \in (0, t] \\ 0 & x \in (t, 1] \end{cases}. \end{equation*} Therefore for each $x$, \begin{equation*} \lim_{t \rightarrow 0^+} \frac{\psi(f + th)(x) - \psi(f)(x)}{t} = \frac{1}{4} I(x \leq 0). \end{equation*} However, for any $t > 0$, \begin{equation*} \sup_{x \in [-1, 1]} \left| (\psi(f + th)(x) - \psi(f)(x)) / t - (1/4) I(x \leq 0) \right| = 1 / 4. \end{equation*} This implies that no derivative exists as an element of $\ell^\infty([-1, 1])$. The functions $f$ and $h$ are well behaved, but because the candidate derivative $I(x \leq 0) / 4$ is not continuous in $x$, uniform convergence of the quotients to the candidate fails.

Directional differentiability of norms of value functions

The lack of uniformity of the convergence to $\psi_f'(h)(\cdot)$ appears to preclude the delta method for uniform inference, because the typical path of analysis would use a well-behaved limit of $r_n(\psi(f_n) - \psi(f))$ and the continuous mapping theorem to find the distribution of statistics applied to the limit. However, real-valued statistics are used to find uniform inference results, so we examine maps that include not only the marginal optimization step but also a functional applied to the resulting value function.

Consider $L_p$ norms (for $1 \leq p \leq \infty$) applied to $\psi(f)$ or $[\psi(f)]_+$. Letting $m$ denote the dominating measure, define $\lambda_j: \ell^\infty(\textnormal{gr} A) \rightarrow \mathbb{R}$ for $j = 1, \ldots, 4$ by

align[align omitted — 450 chars of source]

The maps $\lambda_2$ and $\lambda_4$ are included because one-sided comparisons may also be of interest, and the map $f \mapsto [f]_+$ is pointwise directionally differentiable but not differentiable as a map from $\ell^\infty(X)$ to $\ell^\infty(X)$. We illustrate its use in a stochastic dominance example below.

We make the following definitions in order to construct well-defined derivatives. First, let $\mu(f)$ be the supremum of $f$ over $\textnormal{gr} A$ and let $\sigma(f)$ be the maximin value of $f$ over $\textnormal{gr} A$:

align[align omitted — 161 chars of source]

Define $A_f^\mu \subseteq \textnormal{gr} A$ for any $\epsilon > 0$ to be the set of $\epsilon$-maximizers of $f$ in $(u, x)$ over $\textnormal{gr} A$. Next, let $X_f^\sigma \subseteq X$ collect the $\epsilon$-maximinimizers of $f$ in $x$ over $X$. That is, define

align[align omitted — 260 chars of source]

Also, note that $U_{(-f)}(x, \epsilon)$ is the set of $\epsilon$-minimizers of $f(u, x)$ in $u$ for given $x$.

CarcamoCuevasRodriguez19 show that the map (ref) has directional derivative

equation[equation omitted — 120 chars of source]

and Lemma (ref) in the \hyperref[appn]{Appendix} shows that (ref) has directional derivative

equation[equation omitted — 193 chars of source]

To analyze derivatives of the one-sided $\lambda_2(f)$ and $\lambda_4(f)$, define the contact set

equation[equation omitted — 86 chars of source]

The following two theorems describe the form of the four directional derivatives generally and under the condition that the value function is zero everywhere or on some subset of $X$.

thmSuppose $f, h \in \ell^\infty(\textnormal{gr} A)$, where $A: X \rightrightarrows U$ is non-empty-valued for all $x \in X$. Define $\mu$ by (ref) and $\sigma$ by (ref), and their derivatives $\mu_f'$ and $\sigma_f'$ as in (ref) and (ref). Then: \begin{enumerate} • $\lambda_1(f)$ is Hadamard directionally differentiable and \begin{equation*} \lambda'_{1f}(h) = \begin{cases} \mu_f'(h), & \mu(f) > \sigma(-f) \\ \max \left\{ \mu_f'(h), \sigma_{(-f)}'(-h) \right\}, & \mu(f) = \sigma(-f) \\ \sigma_{(-f)}'(-h), & \mu(f) < \sigma(-f) \end{cases}. \end{equation*} • If $\|\psi(f)\|_\infty = 0$, $\lambda_1(f)$ is Hadamard directionally differentiable and \begin{align*} \lambda'_{1f}(h) &= \max \left\{ \lim_{\epsilon \rightarrow 0^+} \sup_{(u, x) \in U_f(X, \epsilon)} h(u, x), \lim_{\epsilon \rightarrow 0^+} \sup_{x \in X} \inf_{u \in U_f(x, \epsilon)} (-h(u, x)) \right\} \\ &= \lim_{\epsilon \rightarrow 0^+} \sup_{x \in X} \left| \sup_{u \in U_f(x, \epsilon)} h(u, x) \right|. \end{align*} • $\lambda_2(f)$ is Hadamard directionally differentiable and \begin{equation*} \lambda'_{2f}(h) = \begin{cases} \mu_f'(h), & \mu(f) > 0 \\ \left[ \mu_f'(h) \right]_+, & \mu(f) = 0 \\ 0 & \mu(f) < 0 \end{cases}. \end{equation*} • If $\|[\psi(f)]_+\|_\infty = 0$, $\lambda_2(f)$ is Hadamard directionally differentiable and \begin{equation*} \lambda'_{2f}(h) = \begin{cases} \lim_{\epsilon \rightarrow 0^+} \sup_{(u, x) \in U_f(X_0, \epsilon)} \left[ h(u, x) \right]_+, & \mu(f) = 0 \\ 0, & \mu(f) < 0 \end{cases}. \end{equation*} \end{enumerate}

The first part of Theorem (ref) essentially takes advantage of the identity $|x| = \max\{x, -x\}$, and breaks the evaluation of the supremum of the absolute value of the value function into two parts that are well-defined, a maximization and a minimization problem. The second part of the theorem is simpler because the $x \mapsto [x]_+$ map simplifies the evaluation of the supremum.

The next result shows that the $L_p$ functionals for $1 \leq p < \infty$ defined in (ref), are Hadamard directionally differentiable. The form of the derivatives for these maps also change depending on whether $\lambda_j(f) = 0$ or $\lambda_j(f) \neq 0$.

thmSuppose that $f, h \in \ell^\infty(\textnormal{gr} A)$ where $A: X \rightrightarrows U$ is non-empty-valued for all $x \in X$ and $m(X) < \infty$. Recall $\psi'_f$ defined in (ref). Then: \begin{enumerate} • If $\|\psi(f)\|_p \neq 0$, then $\lambda_3$ is Hadamard directionally differentiable and \begin{equation*} \lambda_{3f}'(h) = \| \psi(f) \|_p^{1-p} \int_X \operatorname*{sgn}(\psi(f)(x)) |\psi(f)(x)|^{p-1} \psi_f'(h)(x) d m(x). \end{equation*} • If $\|\psi(f)\|_p = 0$, then $\lambda_3$ is Hadamard directionally differentiable and \begin{equation*} \lambda_{3f}'(h) = \left( \int_X \left| \psi_f'(h)(x) \right|^p d m(x) \right)^{1/p}. \end{equation*} • If $\|[\psi(f)]_+\|_p \neq 0$, then $\lambda_4$ is Hadamard directionally differentiable and \begin{equation*} \lambda_{4f}'(h) = \| [\psi(f)]_+ \|_p^{1-p} \int_X [\psi(f)(x)]_+^{p-1} \psi_f'(h)(x) d m(x). \end{equation*} • If $\|[\psi(f)]_+\|_p = 0$, then $\lambda_4$ is Hadamard directionally differentiable and \begin{equation*} \lambda_{4f}'(h) = \left( \int_{X_0} \left[ \psi_f'(h)(x) \right]_+^p d m(x) \right)^{1/p}. \end{equation*} \end{enumerate}

Theorem (ref) is related to the results of ChenFang19, who dealt with a squared $L_2$ statistic. It is interesting to note here that using an $L_2$ statistic instead of its square results in first-order (directional) differentiability of the map, unlike the squared $L_2$ norm, which has a degenerate first-order derivative but nondegenerate second-order derivative.

We conclude this section by revisiting Examples (ref), (ref), and (ref) and describing the computation of the directional derivatives in the corresponding cases.

excont{ex:bounds}{Dependency bounds} To construct a uniform confidence band for $L_{X+Y}$ we can invert a test for the hypothesis that $L_{X+Y} = L_0$. Recalling the map $\Pi^+$ defined in (ref), consider the statistic $\lambda: (\ell^\infty(\mathbb{R}))^2 \rightarrow \mathbb{R}$ defined by \begin{equation} \lambda(F) = \sup_{x \in \mathbb{R}} |\max\{ \sup_{u \in \mathbb{R}} \Pi^+(F)(u, x), 0\} - L_0(x)|. \end{equation} This is challenging because of the pointwise maximum inside the absolute value. Lemma (ref) in the appendix, which uses calculations similar to those in the proof of Theorem (ref), shows that, letting \begin{align*} X_0 &= \{x \in \mathbb{R} : \sup_{u \in \mathbb{R}} (F_X(u) + F_Y(x - u)) = 1\} \\ X_+ &= \{x \in \mathbb{R} : \sup_{u \in \mathbb{R}} (F_X(u) + F_Y(x - u)) > 1\} \\ U_F(x, \epsilon) &= \{u \in \mathbb{R} : F_X(u) + F_Y(x - u) \geq \sup_{u \in \mathbb{R}} (F_X(u) + F_Y(x - u)) - \epsilon\}, \end{align*} $\lambda$ is Hadamard directionally differentiable and its derivative is, for direction $h = (h_X, h_Y)$, \begin{multline} \lambda_F'(h) = \max \Bigg\{ \lim_{\epsilon \rightarrow 0^+} \sup_{x \in X_0} \sup_{u \in U_F(x, \epsilon)} [h_X(u) + h_Y(x - u)]_+, \\ \lim_{\epsilon \rightarrow 0^+} \sup_{x \in X_+} \left| \sup_{u \in U_F(x, \epsilon)} (h_X(u) - h_Y(x - u)) \right| \Bigg\}. \end{multline} Similarly, we may work with the bound $U_{X+Y}^{-1}$ on the quantile function $F_{X+Y}^{-1}$. When a CDF $F$ is continuously differentiable on a set $\mathcal{T} = [F^{-1}(p) - \epsilon, F^{-1}(q) + \epsilon]$ for $0<p<q<1$ and $\epsilon > 0$, it is fully Hadamard differentiable in $\ell^\infty(\mathcal{T})$ with derivative in direction $h$ at any $\tau \in \mathcal{T}$ given by $-h(F^{-1}(\tau)) / f(F^{-1}(\tau)$, where $f$ is the density associated with $F$ vanderVaartWellner96. Recall the definition of $\tilde{\Pi}^+$ in (ref) and consider, for $F = (F_X, F_Y)$, \begin{equation} \tilde{\lambda}(F) = \sup_{\tau \in \mathcal{T}} |\sup_{u \in (0, \tau)} \tilde{\Pi}^+(F)(u, x) - U_0^{-1}(x)|. \end{equation} Under differentiability of the CDFs, we have \begin{equation*} (\tilde{\Pi}^+)'_F(h)(s, t) = - \frac{h_X(F_X^{-1}(s))}{f_X(F_X^{-1}(s))} - \frac{h_Y(F_Y^{-1}(t-s))}{f_Y(F_Y^{-1}(t-s))}. \end{equation*} Then this fact, Theorem (ref) and the chain rule imply that, letting $\tilde{U}(\tau, \epsilon) = \{u \in (0, \tau) : F_X^{-1}(u) + F_Y^{-1}(\tau - u) \geq \sup_{u \in (0, \tau)} (F_X^{-1}(u) - F_Y^{-1}(\tau - u)) - \epsilon \}$, \begin{equation*} \tilde{\lambda}'_F(h) = \lim_{\epsilon \rightarrow 0} \sup_{\tau \in \mathcal{T}} \left| \sup_{u \in \tilde{U}(\tau, \epsilon)} \left( - \frac{h_X(F_X^{-1}(u))}{f_X(F_X^{-1}(u))} - \frac{h_Y(F_Y^{-1}(\tau-u))}{f_Y(F_Y^{-1}(\tau-u))} \right) \right|. \end{equation*}

excont{ex:stoc_dominance}{Stochastic dominance} Similarly to the previous example, we may test the necessary condition for stochastic dominance using the test statistic \begin{equation} \Lambda = \left( \int [ (L_{A}(x) - U_{B}(x))]_+^2 d m(x) \right)^{1/2}. \end{equation} This test statistic is zero when $L_A \leq U_B$ and greater than zero otherwise. The map $\Pi^-$ is fully differentiable with derivative $(\Pi)_f'(h) = h_1(u) - h_0(u - x)$. Then $\Lambda$ is directionally differentiable and, extending Theorem (ref) to two functions inside the integral, \begin{multline} \Lambda'_G(h) = \bigg( \int_{L_A = U_B} \Big[ \lim_{\epsilon \rightarrow 0^+} \sup_{u \in U_{\Pi^-(G_0, G_A)}(\epsilon, x)} \Pi^-(h_0, h_A) (u, x) \\ - \lim_{\epsilon \rightarrow 0^+} \inf_{u \in U_{-\Pi^-(G_0, G_B)}(\epsilon, x)} \Pi^-(h_0, h_B)(u, x) \Big]_+^2 d m(x) \bigg)^{1/2}. \end{multline}

excont{ex:legendre}{Cost and profit functions} This example shows how some of the technical difficulties encountered in the previous two examples may be alleviated with assumptions on the economic model. Recall that a uniform confidence band for $\pi = \mathcal{L}(c)$ depends on the derivative of $\check{\lambda}(f) = \sup_x |\mathcal{L}(f)(x)|$. Suppose that $(u, x) \in \mathcal{U} \times \mathcal{X}$, a compact set. More generally, the assumptions that firms will not operate with negative profits and that the cost function is 1-coercive (that is, that $\lim_{\|x\| \rightarrow \infty} c(x) / \|x\| = +\infty$) would ensure that $\pi$ is continuous. Given the existence of optimizers and a maximum, write \begin{align*} \check{\lambda}(c) &= \max_p \left| \max_x (p \cdot x - c(x)) \right| \\ &= \max \left\{ \max_{p,x} (p \cdot x - c(x)), \max_p \min_x (c(x) - p \cdot x) \right\}. \end{align*} Since $c$ is estimated and $c$ only depends on $x$, it is sufficient to perturb the above statistic in a direction that does not depend on $p$, i.e., $h(x)$, to calculate derivatives. Define the derivative of $\check{\mu}(c) = \max_{p,x} (p \cdot x - c(x))$ by \begin{equation*} \check{\mu}'_c(h) = \max_{\operatorname*{argmax} (p \cdot x - c(x))} (-h(x)), \end{equation*} and the derivative of $\check{\sigma}(c) = \max_p \min_x (c(x) - p \cdot x)$ by \begin{equation*} \check{\sigma}'_c(h) = \max_{\operatorname*{argmax}_p (\operatorname*{argmin}_x (c(x) - p \cdot x))} \min_{\operatorname*{argmin}_x (c(x) - p \cdot x)} h(x). \end{equation*} With strict convexity of $c$ the argmin expressions in $\check{\mu}'$ and $\check{\sigma}'$ would be single-valued, and these would be simpler. The directional derivative of $\lambda(c)$ is \begin{equation} \lambda_c'(h) = \begin{cases} \check{\mu}'_c(h), & \check{\mu}(c) > \check{\sigma}(c), \\ \check{\mu}'_c(h) \vee \check{\sigma}'_c(h), & \check{\mu}(c) = \check{\sigma}(c), \\ \check{\sigma}'_c(h) & \check{\mu}(c) < \check{\sigma}(c) \end{cases}. \end{equation}

Inference

In this section we derive asymptotic distributions for uniform statistics applied to value functions, and propose bootstrap estimators of their distributions that can be used for practical inference.

Asymptotic distributions

The derivatives developed in the previous subsections can be used in a functional delta method. We make a high-level assumption on the way that observations are related to estimated functions, but conditions under which such convergence holds are well-known, and examples will be shown below.

enumerate[label=A\arabic*., ref=A\arabic*] • Assume that for each $n$ there is a random sample $\{Z_i\}_{i=1}^n$ with $P$ denoting their joint probability distribution. Assume there is a map $\{Z_i\}_{i=1}^n \mapsto f_n$ where $f_n \in \ell^\infty(U \times X)$. In case $\lambda_3$ or $\lambda_4$ is used, also assume $m(X) < \infty$. Furthermore, for some sequence $r_n \rightarrow \infty$, $r_n (f_n - f) \leadsto \mathcal{G}_P$, where $\mathcal{G}_P$ is a tight random element of $\ell^\infty(\textnormal{gr} A)$.

In case an $L_p$ statistic is used (with $p < \infty$), we restrict the space of allowable functions. In order to ensure that $L_p$ statistics are well-defined we make the assumption that the measure of $X$ is finite. While stronger than the assumption used for the supremum statistic, it is sufficient to ensure that the $L_p$ statistics of the value function are finite. This assumption also has the advantage of providing an easy-to-verify condition based on the objective function $f$, rather than the more direct, but less obvious condition $\psi(f) \in L_p(m)$. Lifting this restriction would require some other restriction on the objective function $f$ that is sufficient to ensure the value function is $p$-integrable.\footnote{We attempted to show $p$-integrability by assuming that $f$ is bounded and integrable (note that this would imply $f$ is $p$-integrable in $\textnormal{gr} A$), but were unable to show that this implies that the value function is integrable (in $X$). Note that if one were able to show that $f$ is bounded and integrable, then the elegant results in Kaji19 would apply for the purposes of verifying Assumption (ref).}

Given assumption (ref) and the derivatives of the last section, the asymptotic distribution of test statistics applied to the value function is straightforward.

thmUnder Assumption (ref), \begin{equation*} r_n \left(\lambda_j(f_n) - \lambda_j(f) \right) \leadsto \lambda_{jf}'(\mathcal{G}_P), \end{equation*} where for $j \in \{1, \ldots 4\}$, $\lambda_j$ are defined in (ref) and $\lambda_{jf}'$ are defined in Theorems (ref) and (ref).

This theorem is abstract and hides the fact that the limiting distributions may depend on features of $\lambda_j$ and $f$. Therefore it is only indirectly useful for inference. A resampling scheme is the subject of the next section.

The delta method as used in Theorem (ref) is tailored to arguments that have weak limits in the space $\ell^\infty(\textnormal{gr} A)$, which is an intrinsic feature of the delta method. That implies some limitations on the type of functions that may be considered. For example, functions estimated using kernel methods do not converge to a weak limit as assumed in Assumption (ref). One can in general only show the stochastic order of the sequence $r_n(f_n - f)$, and methods for uniform inference must shown using something other than the delta method. We leave uniform inference for value functions based on less regular estimators to future research.

Our bootstrap technique ahead is designed for resampling under the null hypothesis that $\lambda_j(f) = 0$. We assume that when $j = 1$ or $j = 3$ (that is, conventional two-sided statistics are used), the statistics are used to test the null and alternative hypotheses

align[align omitted — 142 chars of source]

Meanwhile, when $j = 2$ or $j = 4$ (the one-sided cases) we assume the hypotheses are

align[align omitted — 144 chars of source]

The probability measure $P$ is assumed to belong to a collection $\mathcal{P}$ which describes the set of probability measures allowed by the model. A few subcollections of $\mathcal{P}$ serve to organize the asymptotic results below. When $P$ is such that $\psi(f) \equiv 0$, we label $P \in \mathcal{P}_{00}^E$. For functional inequalities, the behavior of test statistics under the null is more complicated. If $P$ is such that $\psi(f) \leq 0$ everywhere, then label $P \in \mathcal{P}_0^I$. When $P \in \mathcal{P}_0^I$ makes $\psi(f)(x) = 0$ for at least one $x \in X$, then we label $P \in \mathcal{P}_{00}^I$.

The following corollary combines the derivatives from Theorem (ref) and Theorem (ref) with the result of Theorem (ref) for these distribution classes. Recall that the set $X_0$, defined in (ref), represents the subset of $X$ where $\psi(f)$ is zero.

corUnder Assumption (ref), for $j \in \{1, \ldots 4\}$, \begin{enumerate} • For each $P \in \mathcal{P}_{00}^E$, $r_n (\lambda_1(f_n) - \lambda_1(f)) \leadsto \lim_{\epsilon \rightarrow 0^+} \sup_{x \in X} \left| \sup_{u \in U_f(x, \epsilon)} \mathcal{G}_P(u, x) \right|$. • \begin{enumerate} • For each $P \in \mathcal{P}_{00}^I$, $r_n(\lambda_2(f_n) - \lambda_2(f)) \leadsto \lim_{\epsilon \rightarrow 0^+} \sup_{(u, x) \in U_f(X_0, \epsilon)} \left[ \mathcal{G}_P(u, x) \right]_+.$ • For each $P \in \mathcal{P}_0^I \backslash \mathcal{P}_{00}^I$, $r_n(\lambda_2(f_n) - \lambda_2(f)) \stackrel{p}{\rightarrow} 0$. \end{enumerate} • For each $P \in \mathcal{P}_{00}^E$, $r_n(\lambda_3(f_n) - \lambda_3(f)) \leadsto \left( \int_X \left| \psi_f'(\mathcal{G}_P)(x) \right|^p \textnormal{d} m(x) \right)^{1/p}$. • \begin{enumerate} • For each $P \in \mathcal{P}_{00}^I$, $r_n(\lambda_4(f_n) - \lambda_4(f)) \leadsto \left( \int_{X_0} \left[ \psi_f'(\mathcal{G}_P)(x) \right]_+^p \textnormal{d} m(x) \right)^{1/p}$. • For each $P \in \mathcal{P}_0^I \backslash \mathcal{P}_{00}^I$, $r_n(\lambda_4(f_n) - \lambda_4(f)) \stackrel{p}{\rightarrow} 0$. \end{enumerate} \end{enumerate}

The distributions described in the above corollary under the assumption that $P \in \mathcal{P}_{00}^E$ or $\mathcal{P}_{00}^I$ are those that we emulate using resampling methods in the next section.

Resampling

In the previous section we established asymptotic distributions for the test statistics of interest. However, for practical inference we turn to resampling to estimate their distributions. This section suggests the use of a resampling strategy that was proposed by FangSantos19, and it combines a standard bootstrap procedure with estimates of the directional derivatives $\lambda_{jf}'$. Hadamard directional differentiability also implies that the numerical approximation methods of HongLi18 may be used to estimate these directional derivatives, an approach that we address in simulations (in the supplementary appendix).

The resampling scheme described below is designed to reproduce the null distribution under the assumption that the statistic is equal to zero, or in other terms, that $P \in \mathcal{P}_{00}^E$ or $P \in \mathcal{P}_{00}^I$. This is achieved by restricting the form of the estimates $\hat{\lambda}_{jn}'$. The bootstrap routine is consistent under more general conditions, as in Theorem 3.1 of FangSantos19. However, we discuss behavior of bootstrap-based tests under the assumption that one of the null conditions described in the previous section holds.

All the derivative formulas in Corollary (ref) require some form of an estimate of the near-maximizers of $f$ in $u$, that is, of the set $U_f(x, \epsilon)$ defined in (ref) for various $x$. The set estimators we use are similar to those used in LintonSongWhang10, ChernozhukovLeeRosen13 and LeeSongWhang18 and depend on slowly decreasing sequences of constants to estimate the relevant sets in the derivative formulas. For a sequence $a_n \rightarrow 0^+$ we estimate $U_f(x, \epsilon)$ with the plug-in estimator $U_{f_n}(x, a_n)$. The sequence $a_n$ should decrease more slowly than the rate at which $f_n$ converges uniformly to $f$, an assumption which will be formalized below.

The one-sided estimates $\hat{\lambda}_{2n}'$ and $\hat{\lambda}_{4n}'$ use a second estimate. For another sequence $b_n \rightarrow 0^+$, estimate the contact set $X_0$ defined in (ref) with

equation*[equation* omitted — 233 chars of source]

This is a method used by LintonSongWhang10 of achieving uniformity in the convergence of estimators that depend on contact set estimation under the null hypothesis $P \in \mathcal{P}_{00}^I$.

Define estimators of $\lambda_{1f}'$ and $\lambda_{2f}'$ by

align[align omitted — 262 chars of source]

and

equation[equation omitted — 110 chars of source]

These formulas impose the condition that $P \in \mathcal{P}_{00}^E$ or $\mathcal{P}_{00}^I$ on the derivative estimates, following the forms shown in Corollary (ref). For more general estimation we would need to devise more complex estimators to deal with the variety of limits discussed in Theorems (ref) and (ref), and also the technique of HongLi18 could be used in that case. However, inference is usually of primary interest, and generating a bootstrap distribution that respects the null hypothesis can improve the performance of bootstrap inference procedures. When $P \in \mathcal{P}_{00}^E$ the imposition amounts to the assumption that $\psi(f)(x) = 0$ for each $x$, so that the set of population $\epsilon$-maximizers in $\textnormal{gr} A$ is $U_f(X, \epsilon)$, and the estimator $\hat{\lambda}_{1n}'$ emulates that. The condition $P \in \mathcal{P}_{00}^I$ implies that $X_0$ is not empty, and $\hat{X}_0$ is meant to estimate this contact set, while using all of $X$ when the distribution is not as is hypothesized.

The estimators $\hat{\lambda}_{3n}'$ and $\hat{\lambda}_{4n}'$ are defined similarly to the supremum-norm estimators: let

equation[equation omitted — 144 chars of source]

and

equation[equation omitted — 174 chars of source]

We use an exchangeable bootstrap, which depends on a set of weights $\{W_i\}_{i=1}^n$ that put probability mass $W_i$ at each observation $Z_i$. This type of bootstrap describes many well-known bootstrap techniques including sampling with replacement from the observations vanderVaartWellner96.

We make the following assumptions to ensure the bootstrap is well-behaved.

enumerate[label=A\arabic*., ref=A\arabic*] \setcounter{enumi}{1} • Assume $a_n, b_n \rightarrow 0^+$, $a_n r_n \rightarrow \infty$ and $b_n r_n \rightarrow \infty$. • Suppose that for each $n$, $W$ is independent of the data $Z$ and there is a map $\{Z_i, W_i\}_{i=1}^n \mapsto f_n^*$ where $f_n^* \in \ell^\infty(\textnormal{gr} A)$. Suppose $r_n(f_n^* - f_n)$ is asymptotically measurable; for all continuous and bounded $g$, $g(r_n(f_n^* - f_n))$ is a measurable function of $\{W_i\}$ outer almost surely in $\{Z_i\}$; and $r_n(f_n^* - f_n) \leadsto \mathcal{G}_P$, where $\mathcal{G}_P$ was defined in Assumption (ref).

We remark that it is difficult to propose a general method for choosing how to set $a_n$ and $b_n$ in practice. We numerically investigated their effect in simulations. We apply similar rules to those used in LintonSongWhang10 and LeeSongWhang18. However, in general situations, where the variance of the objective function changes significantly over its domain, one would need to estimate the pointwise standard deviation of the function. For example, FanPark10 provide pointwise asymptotic distributions for estimated dependency bound functions for $F_{X-Y}$. Nevertheless, we leave this extension for future research.

Assumption (ref) is used to ensure consistency of the $\epsilon$-maximizer set estimators used in the bootstrap algorithm. Assumption (ref) is a condensed version of Assumption 3 of FangSantos19 to ensure that $r_n(f_n^* - f_n)$ behaves asymptotically like $r_n(f_n - f)$.

Resampling routine to estimate the distribution of $r_n(\lambda_j(f_n) - \lambda_j(f))$

enumerate• Estimate $\hat{\lambda}_{jn}'$ using sample data and formulas (ref)-(ref) above.

Then repeat steps (ref)-(ref) for $r = 1, \ldots R$:

enumerate\setcounter{enumi}{1} • Use an exchangeable bootstrap to construct $r_n(f_n^* - f_n)$. • Calculate the resampled test statistic $\lambda^*_r = \hat{\lambda}_{jn}'(r_n(f_n^* - f_n))$ using the estimate from step 1.

Finally,

enumerate\setcounter{enumi}{3} • Let $\hat{q}_{\lambda^{*}}(1-\alpha)$ be the $(1-\alpha)$-th sample quantile from the bootstrap distribution of $\{\lambda_r^*\}_{r=1}^R$, where $\alpha \in (0, 1)$ is the nominal size of the test. Reject the null hypothesis if $r_n \lambda_j(f_n)$ is larger than $\hat{q}_{\lambda^{*}}(1-\alpha)$.

The consistency of this resampling procedure under the null hypothesis is summarized in the following theorem. To discuss weak convergence it is easiest to use the space of bounded Lipschitz functions, which indicates weak convergence. That is, defining $BL_1(\mathbb{R}) = \{g \in C_b(\mathbb{R}) : \sup_x |g(x)| \leq 1, \sup_{x \neq y} |g(x) - g(y)| \leq |x - y|\}$, $X_n$ converges weakly to $X$ if and only if $\sup_{g \in BL_1(\mathbb{R})} |\textnormal{E} \left[ g(X_n) \right] - \textnormal{E} \left[ g(x) \right]| \rightarrow 0$ vanderVaartWellner96.

thmUnder Assumptions (ref)-(ref), if either $j \in \{1, 3\}$ and $P \in \mathcal{P}_{00}^E$, or $j \in \{2, 4\}$ and $P \in \mathcal{P}_{00}^I$, then for any $\epsilon > 0$, \begin{equation*} \limsup_{n \rightarrow \infty} P \left\{ \sup_{g \in BL_1(\mathbb{R})} \left| E \left[ g \left( \hat{\lambda}_{jn}' \left( r_n (f_n^* - f_n) \right) \right) \big| \{Z_i\}_{i=1}^n \right] - E \left[ g \left( \lambda_{jf}'(\mathcal{G}_P) \right) \right] \right| > \epsilon \right\} = 0. \end{equation*}

We conclude this section by revisiting Examples (ref), (ref), and (ref) and describing the resampling procedures in the corresponding cases.

excont{ex:bounds}{Dependency bounds} Recall that tests using the statistic $\lambda(F)$ defined in (ref) may be inverted to find a uniform confidence band for $L_{X+Y}$ (here we focus on the lower dependency bound for the CDF $F_{X+Y}$ only). Suppose that we observe two independent random samples, $\{X_i\}_{i=1}^{n_X}$ and $\{Y_i\}_{i=1}^{n_Y}$ and $\mathbb{F}_n = (\mathbb{F}_{Xn}, \mathbb{F}_{Yn})$ be their empirical CDFs. Then \begin{equation} \lambda(\mathbb{F}_n) = \sup_{x \in \mathbb{R}} | \mathbb{L}_n(x) - L_0(x) |, \end{equation} where $\mathbb{L}_n(x) = \mathbb{L}_n(x, \mathbb{F}_n) = \max \left\{ \sup_{u \in \mathbb{R}} (\mathbb{F}_{Xn}(u) + \mathbb{F}_{Yn}(x - u)) - 1, 0 \right\}$. Standard conditions ensure that, with $n = n_X + n_Y$, $\sqrt{n}(\mathbb{F}_n - F) \leadsto \mathcal{G}_F$, where $\mathcal{G}_F$ is a Gaussian process in $(\ell^\infty(\mathbb{R}))^2$. Under the hypothesis that $L_0$ is the true lower bound function, $\sqrt{n} \lambda(\mathbb{F}_n) \leadsto \lambda'_F(\mathcal{G}_F)$, where $\lambda_F'$ is defined in (ref). This derivative needs to be estimated as in (ref), estimating the set-valued map $U_F(x, \epsilon)$ with the plug-in estimate \begin{equation} U_{\mathbb{F}_n}(x, a_n) = \left\{ u \in \mathbb{R} : \max\{ \mathbb{F}_{Xn}(u) + \mathbb{F}_{Yn}(x - u), 0 \} \geq \mathbb{L}_n(x) - a_n \right\}, \end{equation} where $a_n$ is a sequence that converges slowly to zero, and with an additional slowly-decreasing sequence $b_n$ that is used to determine the regions where $\Pi^+(\mathbb{F}_n)$ falls in the interval $(-b_n, b_n)$ or is above $b_n$ to estimate $\hat{X}_0$ and $\hat{X}_+$. Labeling the estimated derivative $\hat{\lambda}_F'$, slight modifications of Corollary (ref) and Theorem (ref) imply that for all $q \geq 0$, \begin{equation*} P \left\{ \hat{\lambda}_F' (\sqrt{n} (\mathbb{F}_n^* - \mathbb{F}_n)) \leq q \; \Big| \{X_i\}_{i=1}^{n_X}, \{Y_i\}_{i=1}^{n_Y} \right\} \stackrel{p}{\rightarrow} P \left\{ \lambda_F' (\mathcal{G}_F) \leq q \right\}. \end{equation*} We can find a critical value of the asymptotic distribution using the bootstrap by estimating \begin{equation*} \hat{q}_{\lambda^*}(1-\alpha) = \min \left\{ q : P \left\{ \hat{\lambda}_F' (\sqrt{n} (\mathbb{F}_n^* - \mathbb{F}_n)) \leq c \; \Big| \{X_i\}_{i=1}^{n_X}, \{Y_i\}_{i=1}^{n_Y} \right\} \geq 1 - \alpha \right\} \end{equation*} by simulation. For $r = 1, \ldots R$, let $\lambda_r^* = \hat{\lambda}_F' (\sqrt{n} (\mathbb{F}_n^* - \mathbb{F}_n))$, and let $\hat{q}_{\lambda^*}(1-\alpha)$ be the $(1-\alpha)$-th quantile of the bootstrap sample $\{\lambda^*_r\}_{r=1}^R$. Simulation evidence presented in the supplementary appendix verifies that the coverage probability of these intervals is accurate in moderate sample sizes by examining the empirical rejection probabilities of the associated tests.

excont{ex:stoc_dominance}{Stochastic dominance} Suppose that three independent random samples $\{X_{0i}\}_{i=1}^n$, $\{X_{Ai}\}_{i=1}^n$ and $\{X_{Bi}\}_{i=1}^n$ are observed. Let $\mathbb{G}_n = (\mathbb{G}_{0n}, \mathbb{G}_{An}, \mathbb{G}_{Bn})$ be their empirical distribution functions, and define the sample statistic \begin{equation} \hat{\Lambda}(\mathbb{G}_n) = \left( \int_\mathbb{R} \left[ \mathbb{L}_{An}(x) - \mathbb{U}_{Bn}(x) \right]_+^2 d m(x) \right)^{1/2}, \end{equation} where $\mathbb{L}_{An}$ is the plug-in estimate of $L_A$ using $\mathbb{G}_{0n}$ and $\mathbb{G}_{An}$ and $\mathbb{U}_{Bn}$ is the plug-in estimate of $U_B$ using $\mathbb{G}_{0n}$ and $\mathbb{G}_{Bn}$. Under standard conditions, with $n = n_0 + n_A + n_B$, $\sqrt{n}(\mathbb{G}_n - G) \leadsto \mathcal{G}_G$, where $\mathcal{G}_G$ is a Gaussian process in $(\ell^\infty(\mathbb{R}))^3$. The theory above can be extended in a straightforward way to show that $\hat{\Lambda}(\sqrt{n} \mathbb{G}_n) \leadsto \Lambda_G'(\mathcal{G}_G)$, where $\Lambda_G'(h)$ was defined in (ref). To estimate the distribution of $\Lambda_G'(\mathcal{G}_G)$, some estimates of the derivative are required. Given sequence $\{b_n\}$, let \begin{equation*} \hat{X}_0 = \{x \in \mathbb{R}: |\psi(\Pi^-(\mathbb{G}_0, \mathbb{G}_A))(x) - \psi(-\Pi^-(\mathbb{G}_0, \mathbb{G}_B))(x)| \leq b_n\}, \end{equation*} and estimate the near-maximizers in $u$ for each dependency bound and all $x$ as in the previous example. Estimate the distribution by calculating resampled statistics, for $r = 1, \ldots, R$, \begin{multline*} \Lambda^*_r = \sqrt{n} \Bigg( \int_{\hat{X}_0} \bigg[ \sup_{u \in U_{\Pi^-(\mathbb{G}_{0n}, \mathbb{G}_{An})}(x, a_n)} \Pi(\mathbb{G}^*_{0n} - \mathbb{G}_{0n}, \mathbb{G}^*_{An} - \mathbb{G}_{An})(u, x) \\ + \sup_{u \in U_{-\Pi^-(\mathbb{G}_{0n}, \mathbb{G}_{Bn})}(x, a_n)} (-\Pi)(\mathbb{G}_{0n}^* - \mathbb{G}_{0n}, \mathbb{G}_{Bn}^* - \mathbb{G}_{Bn}) (u, x) \bigg]_+^2 d m(x) \Bigg)^{1/2}. \end{multline*} A test can be conducted by comparing $\sqrt{n}\hat{\Lambda}(\mathbb{G}_n)$ to the $(1-\alpha)$-th quantile of the bootstrap distribution $\{\Lambda^*_r\}_{r=1}^R$. A simulation in the supplementary appendix illustrates the accurate size and power of this testing strategy.

excont{ex:legendre}{Cost and profit functions} Now consider estimating a confidence band for a conjugate function. Suppose that we observe $\{X_i\}_{i=1}^n$ and have a parametric model $c(x, \theta)$, estimated using $c(x, \hat{\theta})$. Supposing that $\theta \in \mathbb{R}^p$, this is a VC-subgraph class of functions and $\sqrt{n}(c(\cdot, \hat{\theta}) - c(\cdot, \theta)) \leadsto \mathcal{G}_c$ where $\mathcal{G}_c$ is a Gaussian process vanderVaart98. To estimate the directional derivative $\lambda'(\mathcal{G}_c)$ where $\lambda'$ is defined in (ref), we need to estimate several quantities. Let $\widehat{\check{\mu}} = \max_{p, x} (px - c(x, \hat{\theta}))$ and $\widehat{\check{\sigma}} = \max_p \min_x (c(x) - px)$. For a slowly decreasing sequence $a_n$, decide whether $\widehat{\check{\mu}} + a_n < \widehat{\check{\sigma}}$, $|\widehat{\check{\mu}} - \widehat{\check{\sigma}}| \leq a_n$ or $\widehat{\check{\mu}} > \widehat{\check{\sigma}} + a_n$. Given this determination, estimate the derivative $\check{\lambda}'(\mathcal{G}_c)$ by $\widehat{\check{\sigma}}'(\mathcal{G}_c^*)$, $\max\{ \widehat{\check{\mu}}'(\mathcal{G}_c^*), \widehat{\check{\sigma}}'(\mathcal{G}_c^*)\}$ or $\widehat{\check{\mu}}'(\mathcal{G}_c^*)$ respectively, where $\mathcal{G}_c^* = \sqrt{n}(c(\cdot, \hat{\theta}^*) - c(\cdot, \hat{\theta}))$. Given $\check{q}_{\lambda^*}(1-\alpha/2)$, a quantile from the bootstrap sample $\{\check{\lambda}^*_r\}_{r=1}^R$ calculated as just described, we could estimate a symmetric level $(1-\alpha)$ uniform confidence band for $\mathcal{L}(c)$ using $\mathcal{L}(c(\cdot, \hat{\theta})) \pm \check{q}_{\lambda^*}(1-\alpha/2)$.

Local size control

It is also of interest to examine how these tests behave under sequences of distributions local to distributions that satisfy the null hypothesis. We consider sequences of local alternative distributions $\{P_n\}$ such that for each $n$, $\{Z_i\}_{i=1}^n$ are distributed according to $P_n$, and $P_n$ converges towards a limit $P$ that satisfies the null hypothesis. To describe this process, for $t \geq 0$ define a path $t \mapsto P_t$, where $P_t$ is an element of the space of distribution functions $\mathcal{P}$, such that

equation[equation omitted — 191 chars of source]

where the score function $h \in L_2(F)$ satisfies $\textnormal{E} \left[ h \right] = 0$ and $P_0$ is a distribution that satisfies the null hypothesis. The direction that the sequence approaches the null is described asymptotically by the score $h$. We assume that by letting $t = c / r_n$ for $c \in \mathbb{R}$, we can parameterize distributions that are local to $P_0$ and for $t \geq 0$ denote $f(P_t)$ as the function $f$ under distribution $f(P_t)$ so that the unmarked $f$ described above can be rewritten $f = f(P_0)$. See, e.g., vanderVaartWellner96 for more details.

The following assumption ensures that $f$ remains suitably regular under such local perturbations to the null distribution.

enumerate[label=A\arabic*., ref=A\arabic*] \setcounter{enumi}{3} • For all $c \in \mathbb{R}$, \begin{enumerate} • For given $c$, there exists some $f'(c) \in \ell^\infty(\textnormal{gr} A)$ such that $\|r_n (f(P_{c/r_n}) - f(P_0)) - f'(c)\|_\infty \rightarrow 0$, where $P_{c/r_n}$ satisfy (ref). • $r_n(f_n - f(P_{c/r_n})) \leadsto \mathcal{G}_{P_0}$ in $\ell^\infty(\textnormal{gr} A)$, where for each $n$, $\{Z_i\}_{i=1}^n \sim P_{c/r_n}$. \end{enumerate}

Both parts of Assumption (ref) ensure that $f_n$ behave regularly as distributions drift towards $P_0 \in \mathcal{P}_{00}^E$ or $P_0 \in \mathcal{P}_{00}^I$. These additional regularity conditions allow us to describe the size and local power properties of the test statistics.

thmUnder Assumptions (ref)-(ref): \begin{enumerate} • Suppose $j = 1$ or $3$, $P_0 \in \mathcal{P}_{00}^E$ and for each $n$, the observations $\{Z_i\}_{i=1}^n$ are distributed according to $P_n = P_{c / r_n} \in \mathcal{P}_{00}^E$. Suppose that the CDF of $\lambda_{jf}'(\mathcal{G}_{P_0})$ is continuous and increasing at $q_{1-\alpha}$, its $(1-\alpha)$-th quantile. Then \begin{equation*} \limsup_{n \rightarrow \infty} P_n \left\{ r_n \lambda_j(f_n) > \hat{q}_{\lambda^*}(1-\alpha) \right\} \geq P_0 \left\{ \lambda_{jf}'( \mathcal{G}_{P_0} + f'(c) ) > q_{1-\alpha} \right\}. \end{equation*} This holds with equality if $q_{1-\alpha}$ is a continuity point of the CDF of $\lambda_{jf}'( \mathcal{G}_{P_0} + f'(c) )$. If $c = 0$ then the limiting rejection probability is equal to $\alpha$. • Suppose $j = 2$ or $4$, $P_0 \in \mathcal{P}_{00}^I$ and for each $n$, the observations $\{Z_i\}_{i=1}^n$ are distributed according to $P_n = P_{c/r_n} \in \mathcal{P}_0^I$. Suppose that the CDF of $\lambda_{jf}'(\mathcal{G}_{P_0})$ is continuous and increasing at its $(1-\alpha)$-th quantile. Then \begin{equation*} \limsup_{n \rightarrow \infty} P_n \left\{ r_n \lambda_j(f_n) > \hat{q}_{\lambda^*}(1-\alpha) \right\} \leq \alpha. \end{equation*} \end{enumerate}

Theorem (ref) shows that the size of tests can be controlled locally to the null region (and the nominal rejection probability matches the intended probability) in only some cases. In particular, the tests for the null that $P_0 \in \mathcal{P}_{00}^E$ cannot be shown to control local size without more information about the direction of the local alternatives. This is because we must make assumptions about the underlying objective functions without being able to make similar assumptions about the corresponding value functions --- note that Assumption (ref) is made with regard to $f(P)$, not $\psi(f(P))$. To see where the problem lies we may take $\lambda_1$ as an example. The result of Theorem (ref) shows that when $P_n \in \mathcal{P}_{00}^E$ for each $n$,

equation*[equation* omitted — 95 chars of source]

Although $P \in \mathcal{P}_{00}^E$ implies $\lambda_f'(f'(c)) = 0$, as shown in the proof of Theorem (ref), it does not imply that $f'(c) \equiv 0$. For example, if $f$ are CDFs corresponding to a location-shift family of distributions, i.e., for all $c \in \mathbb{R}$, $G_{c}(x) = G_0(x - c)$ with differentiable densities that have square-integrable derivatives. Then since the $h$ in (ref) is $cg'_0(x) / g_0(x)$, the $f'(c)$ of assumption (ref) would be $-cg_0(x)$. This non-zero derivative can cause problems when optimizing it interacts with the absolute value function. Generally $|\sup_u (f(u) + g(u))| \not\leq |\sup_u f(u)| + |\sup_u g(u)|$ (take $f$ and $g$ strictly less than zero), and we do not achieve a simplification that would imply the probability on the right-hand side is less than or equal to $\alpha$. We cannot assume regularity of $\psi(f)$ because this would require the existence of a derivative $\psi_f'(\cdot) \in \ell^\infty(X)$, but this derivative cannot generally exist. It may be that size control could be shown for special cases (for example when theory implies a shape for the value function), but not in general. On the other hand, for the same reason, the size of one-sided test statistics is as intended. That is, when the value function is found through maximization, the positive-part map in the definitions of $\lambda_2$ and $\lambda_4$ interacts with it in a way that maintains size control. In contrast, the negative-part map would not guarantee size control. A stronger result in this vein will be shown in the next section.

The inequality in the second part of this theorem results from the one-sidedness of the test statistics $\lambda_2$ and $\lambda_4$. This is related to a literature in econometrics on moment inequality testing. Tests may exhibit size that is lower than nominal for local alternatives that are from the interior of the null region. A few possible solutions to this problem have been proposed. For example, one might evaluate the region where the moments appear to hold with equality, which leads to contact set estimates like in the bootstrap routine described above LintonSongWhang10. Alternatively, we may alter the reference distribution by shifting it in the regions where equality does not seem to hold AndrewsShi17.

Uniform size control

Uniformity of inference procedures in the data distribution was introduced in GineZinn91 and SheehyWellner92, and uniformity has become a topic of great interest in the econometrics literature (see, for example, AndrewsGuggenberger10, LintonSongWhang10, RomanoShaikh12 or HongLi18 for an application like ours). Under stronger assumptions that ensure the underlying $f_n$ converge weakly to a limiting process uniformly over the set of possible data distributions, some of the above results can be extended from size control under local alternative distributions to size control that holds uniformly over the class of distributions that satisfy the null. To ensure uniformity we assume the following regularity conditions are satisfied.

enumerate[label=A\arabic*., ref=A\arabic*] \setcounter{enumi}{4} • Suppose that $r_n (f_n - f) \leadsto \mathcal{G}_P$ uniformly over $P \in \mathcal{P}$: \begin{equation*} \limsup_{n \rightarrow \infty} \sup_{P \in \mathcal{P}} \sup_{g \in BL_1(\ell^\infty(gr A))} \left| E \left[ g(r_n(f_n - f(P))) \right] - E \left[ g(\mathcal{G}_P) \right] \right| = 0. \end{equation*} • Suppose that the limiting process $\mathcal{G}_P$ is tight uniformly over $P \in \mathcal{P}$: \begin{equation*} \lim_{M \rightarrow \infty} \sup_{P \in \mathcal{P}} P \left\{ \|\mathcal{G}_P\|_\infty \geq M \right\} = 0. \end{equation*} • Suppose that on $\mathcal{P}_{00}^I$, the processes $\mathcal{G}_P$ are regular for $\lambda_{2f}'$ or $\lambda_{4f}'$ in the following sense: for any $0 < \alpha < 1/2$, there exists a $\underline{q} > 0$ depending only on $\alpha$ such that \begin{equation*} \sup_{P \in \mathcal{P}_{00}^I} P \left\{ \lambda_{jf}'(\mathcal{G}_P) < q \right\} < 1 - \alpha \end{equation*} and for any $s > 0$, \begin{equation*} \limsup_{\eta \rightarrow 0^+} \sup_{P \in \mathcal{P}_{00}^I} P \left\{ \left| \lambda_{jf}'(\mathcal{G}_P) - s \right| \leq \eta \right\} = 0. \end{equation*} In addition, suppose that over $\mathcal{P}_0^I \backslash \mathcal{P}_{00}^I$, the processes $\mathcal{G}_P$ are regular for statistics related to the entire domain $X$ instead of $X_0$ --- letting $\tilde{\lambda}_{2f}'(h) = \lim_{\epsilon \rightarrow 0^+} \sup_{(u, x) \in U_f(X, \epsilon)} [h(u, x)]_+$ and $\tilde{\lambda}_{4f}'(h) = (\int_X \psi_f'(h)(x) \textnormal{d} m(x))^{1/p}$, assume that the above two displays are satisfied with $\tilde{\lambda}_{jf}'(\mathcal{G}_P)$ in the place of $\lambda_{jf}'(\mathcal{G}_P)$. • Suppose that conditional on the observations $\{Z_i\}_{i=1}^n$, $r_n(f_n^* - f_n) \leadsto \mathcal{G}_P$ uniformly over $P \in \mathcal{P}$: for all $\epsilon > 0$, \begin{equation*} \limsup_{n \rightarrow \infty} \sup_{P \in \mathcal{P}} P \left\{ \sup_{g \in BL_1(\ell^\infty(\textnormal{gr} A))} \left| \textnormal{E} \left[ g(r_n(f_n^* - f_n)) | \{Z_i\}_{i=1}^n \right] - \textnormal{E} \left[ g(\mathcal{G}_{P}) \right] \right| > \epsilon \right\} = 0. \end{equation*}

The above assumptions help define the collection $\mathcal{P}$ on which we may assert that inference is valid uniformly for $P \in \mathcal{P}$. Assumptions (ref) and (ref) require that the sample objective functions converge to well-behaved limits uniformly on $\mathcal{P}$. Assumption (ref) requires that also the limiting distributions of the sample statistics $r_n(\lambda_j(f_n) - \lambda_j(f))$ are regular enough over $\mathcal{P}_0^I$ that their CDFs do not have extremely large atoms (i.e., extending above the $(1-\alpha)$-th quantile of the distribution) and are strictly increasing at the relevant critical value. This assumption may be simplified when the $\mathcal{G}_P$ are Gaussian since the functionals $\lambda_2$ and $\lambda_4$ are convex DavydovLifshitsSmorodina98. Assumption (ref) requires that the bootstrap objective functions also converge uniformly, as the original sample functions do. This last assumption may be implied by Assumption (ref), see Lemma A.2 of LintonSongWhang10, for example.

thmUnder Assumptions (ref)-(ref) and (ref)-(ref), for $j = 2$ or $4$, \begin{equation*} \limsup_{n \rightarrow \infty} \sup_{P \in \mathcal{P}_0^I} P \left\{ r_n \lambda_j(f_n) > \hat{q}_{\lambda^*}(1-\alpha) \right\} \leq \alpha. \end{equation*}

Uniformity for $P \in \mathcal{P}$ is maintained by the one-sided statistics when optimization matches the “direction” of the test. That is, when the value function matches the positive-part map in $\lambda_2$ and $\lambda_4$, we have uniformity over $P \in \mathcal{P}$, and analogously uniform inference would be possible for minimization and the negative-part map. Uniformity can be lost when optimization and the direction of the test do not match, which may be the case with two-sided tests.

In the next section we illustrate the usefulness of our results by providing details on the construction of uniform confidence bands around bound functions for the CDF of a treatment effect distribution. See the first section of the supplemental appendix for numerical simulations evaluating the finite-sample performance of the tests developed in Examples (ref) and (ref).

The treatment effect distribution of job training on wages

This section illustrates the inference methods with an evaluation of a job training program. We construct upper and lower bounds for both the distribution and quantile function of the treatment effects, confidence bands for these bound function estimates, and describe a few inference results. This application uses an experimental job training program data set from the National Supported Work (NSW) Program, which was first analyzed by LaLonde86 and later by many others, including HeckmanHotz89, DehejiaWahba99, SmithTodd01, SmithTodd05, Imbens03, and Firpo07.

Recent studies in statistical inference for features of the treatment effects distribution in the presence of partial identification include, among others, FirpoRidder08, FirpoRidder19, FanPark09, FanPark10, FanWu10, FanPark12, FanShermanShum14, FanGuerreZhu17. These studies have concentrated on distributions of finite-dimensional functionals of the distribution and quantile functions, including these functions themselves evaluated at a point. Additional work includes, among others, GautierHoderlein11, ChernozhukovLeeRosen13, Kim14 and ChesherRosen15. Each of these papers provides pointwise inference methods for bounds on the distribution or quantile functions, and often for more complex objects. We contribute to this literature by applying the general results in this paper to provide uniform inference methods for the bounds developed by Makarov and others, and hope that it may indicate the direction that pointwise inference for bounds in more involved models may be extended to be uniformly valid.

The data set we use is described in detail in LaLonde86. We use the publicly available subset of the NSW study used by DehejiaWahba99. The program was designed as an experiment where applicants were randomly assigned into treatment. The treatment was work experience in a wide range of possible activities, such as learning to operating a restaurant or a child care center, for a period not exceeding 12 months. Eligible participants were targeted from recipients of Aid to Families With Dependent Children, former addicts, former offenders, and young school dropouts. The NSW data set consists of information on earnings and employment (outcome variables), whether treated or not.\footnote{The data set also contains background characteristics, such as education, ethnicity, age, and employment variables before treatment. Nevertheless, since we only use the experimental part of the data we refrain from using this portion of the data.} We consider male workers only and focus on earnings in 1978 as the outcome variable of interest. There are a total of 445 observations, where 260 are control observations and 185 are treatment observations. Summary statistics for the two parts of the data are presented in Table (ref).

table[table omitted — 898 chars of source]

To provide a more complete overview of the data, we also compute the empirical CDFs of the treatment and control groups in Figure (ref). From this figure we note that the empirical treatment CDF stochastically dominates the empirical control CDF, and that there are a large number of zeros in each sample. In particular, $\mathbb{F}_{tr,n}(0) \approx 0.24$ and $\mathbb{F}_{co,n}(0) \approx 0.35$.

figure[figure omitted — 429 chars of source]

Suppose a binary treatment is independent of two potential outcomes $(X_{co}, X_{tr})$, where $X_{co}$ denotes outcomes under a control regime and $X_{tr}$ denotes outcomes under a treatment, and $X_{co}$ and $X_{tr}$ have marginal distribution functions $F_{co}$ and $F_{tr}$ respectively. Suppose that interest is in the distribution of the treatment effect $\Delta = X_{tr} - X_{co}$ but we are unwilling to make any assumptions regarding the dependence between $X_{co}$ and $X_{tr}$. In this section we study the relationship between the identifiable functions $F_{co}$ and $F_{tr}$ and functions that bound the distribution function of the unobservable random variable $\Delta$.

$F_\Delta(\cdot)$ is not point-identified because the full bivariate distribution of $(X_{co}, X_{tr})$ is unidentified and the analyst has no knowledge of the dependence between potential outcomes. However, $F_\Delta$ can be bounded. Suppose we observe samples $\{X_{ki}\}_{i=1}^{n_k}$ for $k \in \{co, tr\}$. Define the empirical lower and upper bound functions

align[align omitted — 268 chars of source]

These are plug-in estimates of analogous population bounds $L_\Delta$ and $U_\Delta$, and are similar to Example (ref) because $F_\Delta = F_{(X_{tr}-X_{co})}$. We make the following assumptions on the observed samples of treatment and control observations.

enumerate[label=B\arabic*., ref=B\arabic*] • The observations $\{X_{co,i}\}_{i=1}^{n_{co}}$ and $\{X_{tr,i}\}_{i=1}^{n_{tr}}$ are iid samples and independent of each other and are distributed with marginal distribution functions $F_{co}$ and $F_{tr}$ respectively. Refer to the pair of distribution functions and their empirical distribution function estimates as $F = (F_{co}, F_{tr})$ and $\mathbb{F}_n = (\mathbb{F}_{co,n}, \mathbb{F}_{tr,n})$. • The sample sizes $n_{co}$ and $n_{tr}$ increase in such a way that $n_k / (n_{co} + n_{tr}) \rightarrow \nu_k$ as $n_{co}, n_{tr} \rightarrow \infty$, where $0 < \nu_k < 1$ for $k \in \{co, tr\}$. Define $n = n_{co} + n_{tr}$.

Under Assumptions (ref) and (ref), it is a standard result vanderVaart98 that for $k \in \{co, tr\}$, $\sqrt{n_k} (\mathbb{F}_{kn} - F_k) \leadsto \mathcal{G}_k$, where $\mathcal{G}_{co}$ and $\mathcal{G}_{tr}$ are independent $F_{co}$- and $F_{tr}$-Brownian bridges, that is, mean-zero Gaussian processes with covariance functions $\rho_k(x, y) = F_k(x \wedge y) - F_k(x) F_k(y)$. This implies in turn that $ \sqrt{n} (\mathbb{F}_n - F) \leadsto \mathcal{G}_F = ( \mathcal{G}_{co} / \sqrt{\nu_{co}}, \mathcal{G}_{tr} / \sqrt{\nu_{tr}} )$, where $\mathcal{G}_F$ is a mean-zero Gaussian process with covariance process $\rho_F(x, y) = \text{diag}\{\rho_k(x, y) / \nu_k\}$.

Now we focus on the calculation of a uniform confidence band for $L$ only. Because the bootstrap algorithm was described previously we only verify that the regularity conditions for this plan hold. Assumptions (ref) and (ref), along with the above discussion of the weak convergence of $\sqrt{n}(\mathbb{F}_n - F)$ imply that Assumption (ref) is satisfied. The class of functions $\{ (I(X \leq x), I(Y \leq y), x,y \in \mathbb{R} \}$ is uniform Donsker and the sequences $a_n$ and $b_n$ can be chosen to satisfy Assumption (ref) by making them converge to zero more slowly than $n^{-1/2}$. Assuming the weights are independent of the observations, assumption (ref) is satisfied by Lemma A.2 of LintonSongWhang10, which implies that the bootstrap algorithm described above is consistent, as described in Theorem (ref). It is also straightforward to verify that, under the high-level assumption (ref), both parts of Assumption (ref) are satisfied --- in the language of the assumption, $f(P_{c/\sqrt{n}})$ are the pair $(F_{co}^n, F_{tr}^n)$ under the local probability distribution $P_{c/\sqrt{n}}$, and $f' = (\int_{-\infty}^\cdot h_{co}^c, \int_{-\infty}^\cdot h_{tr}^c)$) for direction $h^c$ indexed by $c$ and $\sqrt{n}(\mathbb{F}_n - F^n) \leadsto \mathcal{G}_F$ vanderVaartWellner96.

The main objective is to provide uniform confidence bands for the CDF of the treatment effects distribution. We calculate the lower and upper bounds for the distribution function using the control and treatment samples as in equations (ref)--(ref). The bounds are computed on a grid with increments of \$100 dollars along the range of the common support of the bounds, which is roughly from -\$40,000 to \$60,000. We focus on the region between -\$40,000 and \$40,000, which contains almost all the observations (there is a single \$60K observation in the treated sample). The results are presented in Figure (ref) and given by the black solid lines in the picture. The most prominent feature is that, as expected, the upper bound for the CDF of treatment effects stochastically dominates the corresponding lower bound.

Next, we compute the uniform confidence bands as described in the text. They are shown in Figure (ref) as the dashed lines around the corresponding solid lines.

figure[figure omitted — 326 chars of source]

Due to the large number of zero outcomes in both samples, these bounds have some interesting features that we discuss further. First, we note that for any $\epsilon$ greater than zero but smaller than the next smallest outcome (about \$45), $\mathbb{F}_{tr,n}(0) - \mathbb{F}_{co,n}(-\epsilon) = 0.24$, which explains the jump in the lower bound estimate near zero (it is really for a point in the grid just above zero). Likewise, for the same $\epsilon$, $\mathbb{F}_{tr,n}(-\epsilon) - \mathbb{F}_{co,n}(0) = -0.35$, which explains the jump in the upper bound just below zero. Without these point masses at zero, both bounds would more smoothly tend towards 0 or 1. Second, the point masses at zero imply another feature of the bounds that can be discerned in the picture. The upper bound to the left of 0 is the same as $1 - \mathbb{F}_{co,n}(-x)$ and the lower bound to the right of zero is the same as $\mathbb{F}_{tr,n}(x)$. Taking the lower bound as an example, for each $x > 0$ find the closest observation from the control sample $y_{co,i^*}$, and set $x^*(x) = X_{co,i^*} + \epsilon$, leading to the supremum $\mathbb{F}_{tr,n}(X_{co,i} + \epsilon)$ at every point where $X_{co,i} + \epsilon < X_{tr,j}$ for all $j$ in the treated sample. It is identical to the empirical treatment CDF for the entire positive part of the support because the treatment first-order stochastically dominates the control. The situation would be different if there were a jump in the empirical control CDF at least as large as the jump in the empirical treatment CDF at zero. Because the opposite is the case for the upper bound, it does change slightly above the zero mark, tending from $1 - \mathbb{F}_{tr,n}(0) + \mathbb{F}_{co,n}(0)$ to 1 as $x$ goes from 0 to the right.

We can also use the confidence bands for the bound functions to construct a confidence band for the true distribution function of treatment effects. This is shown in Figure (ref). A $1 - \alpha$ confidence band can be constructed by using the upper $1 - \alpha / 2$ limit of the upper bound confidence band and the lower $\alpha / 2$ limit of the lower confidence band. This band is a uniform asymptotic confidence band for the true CDF, and uniform over correlation between the potential outcomes between samples. In other words, if $\mathcal{P}$ is the collection of bivariate distributions that have marginal distributions $F_{tr}$ and $F_{co}$, then

equation*[equation* omitted — 164 chars of source]

This confidence band is likely conservative, since

align*[align* omitted — 502 chars of source]

See Kim14 and FirpoRidder19 for a more thorough discussion of the sense in which these bounds are not uniformly sharp for the treatment effect distribution function. We leave more sophisticated, potentially tighter confidence bands for future research. Note that the technique of ImbensManski04 cannot be used to tighten these bounds, because the parameter, a function, could violate the null hypothesis at both sides of the confidence band.

figure[figure omitted — 1,033 chars of source]

We plot some other features in this figure for context. First, the dotted vertical line is positioned at $y = 0$, and it can be seen that we (just) reject the null hypothesis $H_0: P \left\{ \Delta = 0 \right\} = 1$. This hypothesis is closest to non-rejection, and it is clearer that one should reject the null that the treatment effect distribution is degenerate at any other point besides zero. This supports the notion that treatment effect heterogeneity is an important feature of these observations, especially because this band is completely agnostic about the form of the joint distribution. On the other hand, by examining the bands at horizontal levels, it can be seen that for the median effect and a wide interval in the center of the distribution, the hypothesis of zero treatment effect cannot be rejected (although these are uniform bounds and not tests of individual quantile levels).

The final feature in the figure is the dashed curve that represents the estimate that one would make under the assumption of comonotonicity (or rank invariance) --- the assumption that, had an individual been moved from the treatment to the control group, their rank in the control would be the same as their observed rank. Under this strong assumption the quantile treatment effects are the quantiles of the treatment distribution and they can be inverted to make an estimate. Clearly, the estimate under this assumption is just one point-identified treatment effect distribution function of many.

figure[figure omitted — 593 chars of source]

Finally, to provide context for the uniform confidence band results within the literature on inference for bounds, we compare the proposed uniform bands to the pointwise confidence intervals suggested by FanPark10. We used BickelSakov08's automatic procedure to choose subsample size and constructed confidence intervals for each individual point in the grid of the treatment effect support. This collection of pointwise confidence intervals are plotted along with uniform confidence bands in Figure (ref). The results show that the uniform bands are farther from the bound estimates than the set of pointwise confidence intervals.

Conclusion

This paper develops uniform statistical inference methods for optimal value functions, that is, functions constructed by optimizing an objective function in one argument over all values of the remaining arguments. Value functions can be seen as a nonlinear transformation of the objective function. The map from objective function to value function is not Hadamard differentiable, but statistics used to conduct uniform inference are Hadamard directionally differentiable. We establish the asymptotic properties of nonparametric plug-in estimators of these uniform test statistics and develop a resampling technique to conduct practical inference. Examples involving dependency bounds for treatment effects distributions are used for illustration. Finally, we provide an application to the evaluation of a job training program, estimating a conservative uniform confidence band for the distribution function of the program's treatment effect without making any assumptions about the dependence between potential outcomes.

Acknowledgements

We would like to thank Timothy Christensen, Zheng Fang, Hidehiko Ichimura, Lynda Khalaf, Alexandre Poirier, Pedro H.C. Sant'Anna, Andres Santos, Rami Tabri, Brennan Thompson, Tiemen Woutersen, and seminar participants at the University of Kentucky and New York Camp Econometrics XIV for helpful comments and discussions. We would also like to thank the referees who helped improve the paper, especially the comments of one referee that improved subsection 4.4 considerably. All the remaining errors are ours. This research was enabled in part by support from the Shared Hierarchical Academic Research Computing Network (SHARCNET: \url{https://www.sharcnet.ca}) and Compute Canada (\url{http://www.computecanada.ca}). Code used to implement the methods described in this paper and to replicate the simulations can be found at \url{https://github.com/tmparker/uniform_inference_for_value_functions}.