Confidence Intervals for Projections of Partially Identified Parameters

Confidence Intervals forProjections of Partially Identified Parameters

\thispagestyle{empty} \onehalfspacing

\pagenumbering{arabic}

\etocdepthtag.toc{mtchapter} \etocsettagdepth{mtchapter}{subsection} \etocsettagdepth{mtappendix}{none}

Introduction

This paper provides novel confidence intervals for projections and smooth functions of a parameter vector $\theta \in \Theta \subset \mathbb{R}^{d}$, $d<\infty$, that is partially or point identified through a finite number of moment (in)equalities. In addition, we develop a new algorithm for computing these confidence intervals and, more generally, for solving optimization problems with “black box” constraints, and obtain its rate of convergence.

Until recently, the rich literature on inference for moment (in)equalities focused on confidence sets for the entire vector $\theta$, usually obtained by test inversion as

where the test statistic $T_n(\theta)$ aggregates violations of the sample analog of the moment (in)equalities and the critical value $c_{1-\alpha}(\theta)$ controls asymptotic coverage, often uniformly over a large class of data generating processes (DGPs). However, applied researchers are frequently interested in a specific component (or function) of $\theta$, e.g., the returns to education. Even if not, they may simply want to report separate confidence intervals for components of a vector, as is standard practice in other contexts. Thus, consider inference on the projection $p^{\prime }\theta$, where $p$ is a known unit vector. To date, it is common to report as confidence set the corresponding projection of $\mathcal C_n(c_{1-\alpha})$ or the interval

which will miss any “gaps" in a disconnected projection but is much easier to compute. This approach yields asymptotically valid but typically conservative and therefore needlessly large confidence regions. The potential severity of this effect is easily appreciated in a point identified example. Given a $\sqrt{n}$-consistent estimator $\hat{\theta}_n \in \mathbb{R}^d$ with limiting covariance matrix equal to the identity matrix, the usual 95% confidence interval for $\theta_k$ equals $[\hat{\theta}_{n,k}-1.96,\hat{\theta}_{n,k}+1.96]$. Yet the analogy to $CI^{proj}_n$ would be projection of a 95% confidence ellipsoid, which with $d=10$ yields $[\hat{\theta}_{n,k}-4.28,\hat{\theta}_{n,k}+4.28]$ and a true coverage of essentially $1$.

Our first contribution is to provide a bootstrap-based {\em calibrated projection\/} method to largely anticipate and correct for the conservative effect of projection. The method uses an estimated critical level $\hat{c}_{n,1-\alpha}$ calibrated so that the projection of $\mathcal C_n(\hat{c}_{n,1-\alpha})$ covers $p'\theta$ (but not necessarily $\theta$) with probability at least $1-\alpha$. As a confidence region for the true $p'\theta$, one may report this projection, i.e.

or, for computational simplicity and presentational convenience, the interval

We prove uniform asymptotic validity of both over a large class of DGPs.

Computationally, calibration of $\hat{c}_{n,1-\alpha}$ is relatively attractive: We linearize all constraints around $\theta $, so that coverage of $p^\prime\theta$ can be calibrated by analyzing many linear programs. Nonetheless, computing the above objects is challenging in moderately high dimension. This brings us to our second contribution, namely a general method to accurately and rapidly compute confidence intervals whose construction resembles (ref). Additional applications within partial identification include projection of confidence regions defined in CHT, AS, or AShiECMA, as well as (with minor tweaking; see Appendix (ref)) the confidence interval proposed in BCS14_subv and further discussed later. In an application to a point identified setting, fre:rev17 use our method to construct uniform confidence bands for an unknown function of interest under (nonparametric) shape restrictions. They benchmark it against gridding and find it to be accurate at considerably improved speed. More generally, the method can be broadly used to compute confidence intervals for optimal values of optimization problems with estimated constraints.

Our algorithm (henceforth called E-A-M for Evaluation-Approximation-Maximization) is based on the response surface method, thus it belongs to the family of {\em expected improvement algorithms\/} jonesa2001,jonesefficient1998. Bull_Convergence_2011 established convergence of an expected improvement algorithm for unconstrained optimization problems where the objective is a “black box" function. The rate of convergence that he derives depends on the smoothness of the black box objective function. We substantially extend his results to show convergence, at a slightly slower rate, of our similar algorithm for constrained optimization problems in which the constraints are sufficiently smooth “black box" functions. Extensive Monte Carlo experiments (see Appendix (ref) and Section 5 of KMS_2017) confirm that the E-A-M algorithm is fast and accurate.

Relation to existing literature. The main alternative inference prodedure for projections -- introduced in Romano_Shaikh2008aJSPIWP and significantly advanced in BCS -- is based on profiling out a test statistic. The classes of DGPs for which calibrated projection and the profiling-based method of BCS (BCS-profiling henceforth) can be shown to be uniformly valid are non-nested.\footnote{See KMS_2017 for a comparison of the statistical properties of calibrated projection and BCS-profiling, summarized here at the end of Section (ref).}

Computationally, calibrated projection has the advantage that the bootstrap iterates over linear as opposed to nonlinear programming problems. While the “outer" optimization problems in (ref) are potentially intricate, our algorithm is geared toward them. Monte Carlo simulations suggest that these two factors give calibrated projection a considerable computational edge over profiling, though profiling can also benefit from the E-A-M algorithm. Indeed, in Appendix (ref) we replicate the Monte Carlo experiment of BCS and find that adapting E-A-M to their method improves computation time by a factor of about $4$, while switching to calibrated projection improves it by a further factor of about $17$.

In an influential paper, PakesPorterHo2011 also use linearization but, subject to this approximation, directly bootstrap the sample projection. This is valid only under stringent conditions.\footnote{The published version of PPHI, i.e. PPHI_ECMA, does not contain the inference part. KMS_2017 show that calibrated projection can be much simplified under the conditions imposed by PPHI.} Other related articles that explicitly consider inference on projections include ber:mol08, BontempsMagnacMaurin2012E, Kaido12, and KT15. None of these establish uniform validity of confidence sets. CCOT establish uniform validity of MCMC-based confidence intervals for projections, but aim at covering the projection of the entire identified region $\Theta_I(P)$ (defined later) and not just of the true $\theta$. GMO16 use our insight in the context of set identified spatial VARs.

Regarding computation, previous implementations of projection-based inference CilibertoTamer09,Grieco14,dic:mor16 reported the smallest and largest value of $p^\prime \theta$ among parameter values $\theta \in \mathcal C_n(c_{1-\alpha})$ that were discovered using, e.g., grid-search or simulated annealing with no cooling. This becomes computationally cumbersome as $d$ increases because it typically requires a number of evaluation points that grows exponentially with $d$. In contrast, using a probabilistic model, our method iteratively draws evaluation points from regions that are considered highly relevant for finding the confidence interval's end point. In applications, this tends to substantially reduce the number of evaluation points.

Structure of the paper. Section (ref) sets up notation and describes our approach in detail, including computational implementation of the method and choice of tuning parameters. Section (ref) establishes uniform asymptotic validity of $CI_n$, and Section (ref) shows that our algorithm converges at a specific rate which depends on the smoothness of the constraints. Section (ref) reports the results of an empirical application that revisits the analysis in KT15. Section (ref) draws conclusions. The proof of convergence of our algorithm is in Appendix (ref). Appendix (ref) shows that our algorithm can be used to compute BCS-profiling confidence intervals. Appendix (ref) reports the results of Monte Carlo simulations comparing our proposed method with that of BCS. All other proofs, background material for our algorithm, and additional results are in the Online Appendix.\footnote{Appendix (ref) provides convergence-related results and background material for our algorithm and describes how to compute $\hat{c}_{n,1-\alpha}(\theta)$. Appendix (ref) presents the assumptions under which we prove uniform asymptotic validity of $CI_n$. Appendix (ref) verifies, for a number of canonical partial identification problems, the assumptions that we invoke to show validity of our inference procedure and for our algorithm. Appendix (ref) contains the proof of Theorem (ref). Appendix (ref) collects Lemmas supporting this proof. }

Detailed Explanation of the Method

Setup and Definition of $CI_n$

Let $X_i\in\mathcal X\subseteq\mathbb R^{d_X}$ be a random vector with distribution $P$, let $\Theta\subseteq\mathbb R^{d}$ denote the parameter space, and let $m_j:\mathcal{X}\times\Theta\to\mathbb{R}$ for $j=1,\dots,J_1+J_2$ denote known measurable functions characterizing the model. The true parameter value $\theta$ is assumed to satisfy the moment inequality and equality restrictions

The identification region $\Theta_I(P)$ is the set of parameter values in $\Theta$ satisfying (ref)-(ref). For a random sample $\{X_i,i=1,..., n\}$ of observations drawn from $P$, we write

for the sample moments and the analog estimators of the population moment functions' standard deviations $\sigma_{P,j}$. The confidence interval in (ref) then is

with

and similarly for $(-p)$. Henceforth, to simplify notation, we write $\hat{c}_n$ for $\hat{c}_{n,1-\alpha}$. We also define $J\equiv J_1+2J_2$ moments, where $\bar m_{n,J_1+J_2+k}(\theta)=-\bar m_{J_1+k}(\theta)$ for $k=1,\dots,J_2$. That is, we treat moment equality constraints as two opposing inequality constraints.

For a class of DGPs $\mathcal P$ that we specify below, define the asymptotic size of $CI_n$ by\footnote{Here we focus on the confidence interval $CI_n$ defined in (ref). See Appendix (ref) for the analysis of the confidence region given by the mathematical projection in (ref).}

We next explain how to control this size and then how to compute $CI_n$.

Calibration of $\hat{c}_n(\theta)$

Calibration of $\hat{c}_n$ requires careful analysis of the moment restrictions' local behavior at each point in the identification region. This is because the extent of projection conservatism depends on (i) the asymptotic behavior of the sample moments entering the inequality restrictions, which can change discontinuously depending on whether they bind at $\theta$ or not, and (ii) the local geometry of the identification region at $\theta$, i.e. the shape of the constraint set formed by the moment restrictions. Features (i) and (ii) can be quite different at different points in $\Theta_I(P)$, making uniform inference challenging. In particular, (ii) does not arise if one only considers inference for the entire parameter vector, and hence is a new challenge requiring new methods.

To build an intuition, fix $P\in\mathcal P$ and $\theta\in\Theta_I(P)$. The projection of $\theta$ is covered when

Here, we first substituted $\vartheta=\theta+\lambda/\sqrt n$ and took $\lambda$ to be the choice parameter; intuitively, this localizes around $\theta$ at rate $1/\sqrt{n}$. We then make the event smaller by adding the constraint $\lambda \in \rho B^d$, with $B^d \equiv [-1,1]^d$ and $\rho \geq 0$ a tuning parameter. We motivate this step later.

Our goal is to set the probability of (ref) equal to $1-\alpha$. To ease computation, we approximate (ref) by linear expansion in $\lambda$ of the constraint set. For each $j$, add and subtract $\sqrt n E_P[m_j(X_i,\theta+\lambda/\sqrt n)]/\hat\sigma_{n,j}(\theta+\lambda/\sqrt n)$ and apply the mean value theorem to obtain

Here $\mathbb G_{n,j}(\cdot) \equiv \sqrt n (\bar m_{n,j}(\cdot)-E_P[m_j(X_i,\cdot)])/\sigma_{P,j}(\cdot)$ is a normalized empirical process indexed by $\theta\in \Theta$, $D_{P,j}(\cdot)\equiv\nabla_\theta \{E_P[m_j(X_i,\cdot)]/\sigma_{P,j}(\cdot)\}$ is the gradient of the normalized moment, $\gamma_{1,P,j}(\cdot)\equiv E_P(m_j(X_i,\cdot))/\sigma_{P,j}(\cdot)$ is the studentized population moment, and the mean value $\bar\theta$ lies componentwise between $\theta$ and $\theta+\lambda/\sqrt n$.\footnote{The mean value $\bar\theta$ changes with $j$ but we omit the dependence to ease notation.}

We formally establish that the probability of the last event in (ref) can be approximated by the probability that $0$ lies between the optimal values of two stochastic linear programs. The components that characterize these programs can be estimated. Specifically, we replace $D_{P,j}(\cdot)$ with a uniformly consistent (on compact sets) estimator, $\hat{D}_{n,j}(\cdot)$,\footnote{See Online Appendix (ref) for such estimators in some canonical moment (in)equality examples.} and the process $\mathbb G_{n,j}(\cdot)$ with its simple nonparametric bootstrap analog, $\mathbb G^b_{n,j}(\cdot)\equiv n^{-1/2}\sum_{i=1}^{n} (m_{j}(X_{i}^{b},\cdot)-\bar{m}_{n,j}(\cdot))/\hat{\sigma}_{n,j}(\cdot )$.\footnote{BCS approximate $\mathbb{G}_{n,j}(\cdot)$ by $n^{-1/2}\sum_{i=1}^{n} [(m_{j}(X_{i},\cdot)-\bar{m}_{n,j}(\cdot))/\hat{\sigma}_{n,j}(\cdot )]\chi_i$ with $\{\chi_i\sim N(0,1)\}_{i=1}^n$ i.i.d. This approximation is equally valid in our approach, and can be faster as it avoids repeated evaluation of $m_{j}(X^b_{i},\cdot)$.} Estimation of $\gamma_{1,P,j}(\theta)$ is more subtle because it enters (ref) scaled by $\sqrt{n}$, so that a sample analog estimator will not do. However, this specific issue is well understood in the moment inequalities literature. Following AS and others Bugni2009E,Canay2010JE,Stoye09, we shrink this sample analog toward zero, leading to conservative (if any) distortion in the limit. Formally, we estimate $\gamma_{1,P,j}(\theta)$ by $\varphi(\hat{\xi}_{n,j}(\theta))$, where $\varphi:\mathbb{R}^J_{[\pm\infty]} \mapsto \mathbb{R}^J_{[\pm\infty]}$ is one of the Generalized Moment Selection (GMS henceforth) functions proposed by AS,

and $\kappa_n\to \infty$ is a user-specified thresholding sequence.\footnote{A common choice of $\varphi$ is given component-wise by

Restrictions on $\varphi$ and the rate at which $\kappa_n$ diverges are imposed in Assumption (ref). While for concreteness here we write out the “hard thresholding" GMS function, Theorem (ref) below applies to all but one of the GMS functions in AS, namely to $\varphi^1-\varphi^4$, all of which depend on $\kappa_n^{-1}\sqrt n\bar m_{n,j}(\theta)/\hat\sigma_{n,j}(\theta)$. We do not consider GMS function $\varphi^5$, which depends also on the covariance matrix of the moment functions.} In sum, we replace the random constraint set in (ref) with the (bootstrap based) random polyhedral set\footnote{Here, we implicitly assume that $\Theta$ is a polyhedral set. If it is instead defined by smooth convex (in)equalities, these can be linearized too.}

The critical level $\hat{c}_n(\theta)$ to be used in (ref) then is

where $P^*$ denotes the law of the random set $\Lambda_n^b (\theta ,\rho ,c)$ induced by the bootstrap sampling process, i.e. by the distribution of $(X_1^b,\dots,X_n^b)$ conditional on the data. Expression (ref) uses convexity of $\Lambda_n^b (\theta, \rho ,c)$ and reveals that the probability inside curly brackets can be assessed by repeatedly checking feasibility of a linear program.\footnote{We implement a program in $\mathbb{R}^d$ for simplicity but, because $p'\lambda=0$, one could reduce this to $\mathbb{R}^{d-1}$.} We describe in detail in Online Appendix (ref) how we compute $\hat{c}_n(\theta)$ through a root-finding algorithm.

We conclude by motivating the “$\rho$-box constraint" in (ref), which is a major novel contribution of this paper. The constraint induces conservative bias but has two fundamental benefits: First, it ensures that the linear approximation of the feasible set in (ref) by (ref) is used only in a neighborhood of $\theta$, and therefore that it is uniformly accurate. More subtly, it ensures that coverage induced by a given $c$ depends continuously on estimated parameters even in certain intricate cases. This renders calibrated projection valid in cases that other methods must exclude by assumption.\footnote{In (ref), set $(\mathbb{G}_{n,1}^b(\cdot),\mathbb{G}_{n,2}^b(\cdot))\sim N(0,I_2)$, $p=\hat{D}_{n,1}=\hat{D}_{n,2}=(0,1)$, $\varphi_1(\cdot)=\varphi_2(\cdot)=0$, and $\alpha=.05$. Then simple algebra reveals that (with or without $\rho$-box) $\hat{c}_n(\cdot)=\Phi^{-1}(\sqrt{.95}) \approx 1.95$. If $\hat{D}_{n,1}=(0,1-\delta)$ and $\hat{D}_{n,2}=(0,1-\delta)$, then without $\rho$-box we have $\hat{c}_n(\cdot)=\Phi^{-1}(.95)/\sqrt{2}\approx 1.16$ for any small $\delta>0$, and we therefore cannot expect to get $\hat{c}_n(\cdot)$ right if gradients are estimated. With $\rho$-box, $\hat{c}_n(\cdot) \to 1.95$ as $\delta \to 0$, so the problem goes away. This stylized example is relevant because it resembles polyhedral identified sets where one face is near orthogonal to $p$. It violates assumptions in BCS and PPHI.}

Computation of $CI_n$ and of Similar Confidence Intervals

Projection based methods as in (ref) and (ref) have nonlinear constraints involving a critical value which in general is an unknown function, with unknown gradient, of $\theta$. Similar considerations often apply to critical values used to build confidence intervals for optimal values of optimization problems with estimated constraints. When the dimension of the parameter vector is large, directly solving optimization problems with such constraints can be expensive even if evaluating the critical value at each $\theta$ is cheap.

This concern motivates this paper's second main contribution, namely a novel algorithm for constrained optimization problems of the following form:

where $\theta^*$ is an optimal solution of the problem and $g_j(\cdot),j=1,...,J$ as well as $c(\cdot)$ are fixed functions of $\theta$. In our own application, $g_j(\theta)=\sqrt n\bar m_{n,j}(\theta)/\hat \sigma_{n,j}(\theta)$ and, for calibrated projection, $c(\theta)=\hat c_n(\theta)$.\footnote{We emphasize that, in analyzing the computational problem, we take the data, including bootstrap data, as given. Thus, while an econometrician would usually think of $\sqrt n\bar m_{n,j}(\theta)/\hat \sigma_{n,j}(\theta)$ and $\hat c_n(\theta)$ as random variables, for this section's purposes they are indeed just functions of $\theta$.}

The key issue is that evaluating $c(\cdot)$ is costly.\footnote{For simplicity and to mirror our motivating application, we suppose that $g_j(\cdot)$ is easy to compute. The algorithm is easily adapted to the case where it is not. Indeed, in Appendix (ref), we show how E-A-M can be employed to compute BCS-profiling confidence intervals, where the profiled test statistic itself is costly to compute and is approximated together with the critical value.} Our algorithm does so at relatively few values of $\theta$. Elsewhere, it approximates $c(\cdot)$ through a probabilistic model that gets updated as more values are computed. We use this model to determine the next evaluation point but report as tentative solution the best value of $\theta$ at which $c(\cdot)$ was computed, not a value at which it was merely approximated. Under reasonable conditions, the tentative optimal values converge to $p'\theta^*$ at a rate (relative to iterations of the algorithm) that is formally established in Section (ref).

After drawing an initial set of evaluation points that we set to grow linearly with $d$, the algorithm has three steps called E, A, and M below.

\noindentInitialization: Draw randomly (uniformly) over $\Theta$ a set $(\theta^{(1)},...,\theta^{(k)})$ of initial evaluation points. Evaluate $c(\theta^{(\ell)})$ for $\ell=1,...,k-1$. Initialize $L=k$.

\noindentE-Step: Evaluate $c(\theta^{(L)})$ and record the tentative optimal value

with $\bar g(\theta)=\max_{j=1,...,J}g_j(\theta)$.

\noindentA-step: Approximate $\theta\mapsto c(\theta)$ by a flexible auxiliary model. We use a Gaussian-process regression model (or kriging), which for a mean-zero Gaussian process $\zeta(\cdot)$ indexed by $\theta$ and with constant variance $\varsigma^2$ specifies

where $\Upsilon^{(\ell)}=c(\theta^{(\ell)})$ and $K_\beta$ is a kernel with parameter vector $\beta \in \bigtimes_{h=1}^d[\underline{\beta}_h,\overline{\beta}_h]\subset \mathbb{R}^d_{++}$; e.g., $K_\beta(\theta-\theta')=\exp(-\sum_{h=1}^d|\theta_h-\theta'_h|^{2}/\beta_h)$. The unknown parameters $(\mu,\varsigma^2)$ can be estimated by running a GLS regression of $\mathbf\Upsilon=(\Upsilon^{(1)},...,\Upsilon^{(L)})'$ on a constant with the given correlation matrix. The unknown parameters $\beta$ can be estimated by a (concentrated) MLE.

The (best linear) predictor of the critical value and its gradient at $\theta$ are then given by

where $\mathbf r_L(\theta)$ is a vector whose $\ell$-th component is $Corr(\zeta(\theta),\zeta(\theta^{(\ell)}))$ as given above with estimated parameters, $\mathbf Q_L(\theta)=\nabla_\theta \mathbf r_L(\theta)'$, and $\mathbf R_L$ is an $L$-by-$L$ matrix whose $(\ell,\ell')$ entry is $Corr(\zeta(\theta^{(\ell)}),\zeta(\theta^{(\ell')}))$ with estimated parameters. This surrogate model has the property that its predictor satisfies $c_L(\theta^{(\ell)})=c(\theta^{(\ell)}), \ell=1,..., L$. Hence, it provides an analytical interpolation, with analytical gradient, of evaluation points of $c(\cdot)$.\footnote{See details in jonesefficient1998. We use the DACE MATLAB kriging toolbox (\url{http://www2.imm.dtu.dk/projects/dace/}) for this step in our empirical application and Monte Carlo experiments.} The uncertainty left in $c(\cdot)$ is captured by the variance

\noindentM-step: With probability $1-\epsilon$, obtain the next evaluation point $\theta^{(L+1)}$ as

where $\mathbb{EI}_L(\theta)$ is the {\em expected improvement function\/}.\footnote{Heuristically, $\mathbb{EI}_L(\theta)$ is the expected improvement gained from analyzing parameter value $\theta$ for a Bayesian whose current beliefs about $c$ are described by the estimated model. Indeed, for each $\theta$, the maximand in (ref) multiplies improvement from learning that $\theta$ is feasible with this Bayesian's probability that it is.} This step can be implemented by standard nonlinear optimization solvers, e.g. MATLAB's \verb1fmincon1 or \verb1KNITRO1 (see Appendix (ref) for details). With probability $\epsilon$, draw $\theta^{(L+1)}$ randomly from a uniform distribution over $\Theta$. Set $L \leftarrow L+1$ and return to the E-step.

The algorithm yields an increasing sequence of tentative optimal values $p'\theta^{*,L},L=k+1,k+2,...$, with $\theta^{*,L}$ satisfying the {\em true\/} constraints in (ref) but the sequence of evaluation points leading to it obtained by maximization of expected improvement defined with respect to the {\em approximated\/} surface. Once a convergence criterion is met, $p'\theta^{*,L}$ is reported as the end point of $CI_n$. We discuss convergence criteria in Appendix (ref).

The advantages of E-A-M are as follows. First, we control the number of points at which we evaluate the critical value; recall that this evaluation is the expensive step. Also, the initial $k$ evaluations can easily be parallelized. For any additional E-step, one needs to evaluate $c(\cdot)$ only at a single point $\theta^{(L+1)}$. The M-step is crucial for reducing the number of additional evaluation points. To determine the next evaluation point, it trades off “exploitation” (i.e. the benefit of drawing a point at which the optimal value is high) against “exploration” (i.e. the benefit of drawing a point in a region in which the approximation error of $c$ is currently large) through maximizing expected improvement.\footnote{It is also possible to draw multiple points in each iteration Schonlau_Global_1998, as we do in our implementation of the method.} Finally, the algorithm simplifies the M-step by providing constraints and their gradients for program (ref) in closed form, thus greatly aiding fast and stable numerical optimization. The price is the additional approximation step. In the empirical application in Section (ref) and in the numerical exercises of Appendix (ref), this price turns out to be low.

Choice of Tuning Parameters

Practical implementation of calibrated projection and the E-A-M algorithm is detailed in KMST_code. It involves setting several tuning parameters, which we now discuss.

Calibration of $\hat{c}_n$ in (ref) must be tuned at two points, namely the use of GMS and the choice of $\rho$. The trade-offs in setting these tuning parameters are apparent from inspection of (ref). GMS is parameterized by a shrinkage function $\varphi$ and a sequence $\kappa_n$ that controls the rate of shrinkage. In practice, choice of $\kappa_n$ is more delicate. A smaller $\kappa_n$ will make $\Lambda_n^b$ larger, hence increase bootstrap coverage probability for any given $c$, hence reduce $\hat{c}_n$ and therefore make for shorter confidence intervals -- but the uniform asymptotics will be misleading, and finite sample coverage therefore potentially off target, if $\kappa_n$ is too small. We follow the industry standard set by AS and recommend $\kappa_n=\sqrt{\log n}$.

The trade-off in choosing $\rho$ is similar but reversed. A larger $\rho$ will expand $\Lambda_n^b$ and therefore make for shorter confidence intervals, but (our proof of) uniform validity of inference requires $\rho<\infty$. Indeed, calibrated projection with $\rho=0$ will disregard any projection conservatism and (as is easy to show) exactly recovers projection of the AS confidence set. Intuitively, we then want to choose $\rho$ large but not too large.

To this end, we heuristically calibrate $\rho$ based on how much conservative distortion one is willing to accept in well-behaved cases. This distortion -- denote it $\eta$, for which we suggest a numerical value of $0.01$ -- is compared against a bound on conservative distortion that is itself likely to be conservative but data free and trivial to compute. In particular, we set

The underlying heuristic is as follows: If all basic solutions (i.e., intersections of exactly $d$ constraints) that potentially define vertices of $\Lambda^b_n$ realize inside the $\rho$-box, then the $\rho$-box cannot affect the values in (ref) and hence not whether coverage obtains in a given bootstrap sample. Conversely, the probability that at least one basic solution realizes outside the $\rho$-box bounds from above the conservative distortion. This probability is, of course, dependent on unknown parameters. Our data free approximation imputes multivariate standard normal distributions for all basic solutions and Bonferroni adjustment to handle their covariation.\footnote{To reproduce the expression, recall that if $a\equiv\tbinom{J_1+J_2}{d}$ random variables in $\mathbb{R}^d$ are individually multivariate standard normal, then a Bonferroni upper bound on the probability that not all of them realize inside the $\rho$-box equals $a\bigl(1-\left(1-2\Phi(-\rho)\right)^d\bigr).$ Also, if Bonferroni is replaced with an independence assumption, the expression changes to $\rho=\Phi^{-1}\bigl(\tfrac{1}{2}+\tfrac{1}{2}(1-\eta)^{1/ad}\bigr)$. The numerical difference is negligible for moderate $J_1+J_2$.}

The E-A-M algorithm also has two tuning parameters. One is $k$, the initial number of evaluation points. The other is $\epsilon$, the probability of drawing $\theta^{(L+1)}$ randomly from a uniform distribution on $\Theta$ instead of by maximizing $\mathbb{EI}_L$. In calibrated projection use of the E-A-M algorithm there is a single “black box" function, $\hat{c}_n(\theta)$. We therefore suggest setting $k=10d+1$, similarly to the recommendation in jonesefficient1998. In our Monte Carlo exercises we experimented with larger values, e.g. $k=20d+1$, and found that the increased number had no noticeable effect on the computed $CI_n$. If a user applies our E-A-M algorithm to a constrained optimization problem with {\em many\/} “black box" functions to approximate, we suggest using a larger number of initial points.

The role of $\epsilon$ Bull_Convergence_2011 is to trade off the greediness of the $\mathbb{EI}_L$ maximization criterion with the overarching goal of global optimization. sut:bar98 explore the effect of setting $\epsilon=0.1$ and $0.01$ on different optimization problems, and find that for sufficiently large $L$, $\epsilon=0.01$ performs better. In our own simulations we have found that drawing {\em both\/} a uniform point and computing the value of $\theta$ for each $L$ (thereby sidestepping the choice of $\epsilon$) is fast and accurate, and that is what we recommend doing.

Theoretical Results

Asymptotic Validity of Inference

In this section we establish that $CI_n$ is uniformly asymptotically valid in the sense of ensuring that (ref) equals at least $1-\alpha$. The result applies to: (i) Confidence intervals for one projection; (ii) joint confidence regions for several projections, in particular confidence hyperrectangles for subvectors; (iii) confidence intervals for smooth nonlinear functions $f:\Theta \mapsto \mathbb{R}$. Examples of the latter extension include policy analysis and estimation of partially identified counterfactuals as well as demand extrapolation subject to rationality constraints.\footnote{In Appendix (ref), we show that the result actually applies to the mathematical projection in (ref).}

All assumptions can be found in Online Appendix (ref). Assumptions (ref) and (ref) are mild regularity conditions typical in the literature; see, e.g., Definition 4.2 and the corresponding discussion in BCS. Assumption (ref) is based on AS and constrains the GMS function $\varphi(\cdot)$ as well as the rate at which $\kappa_n$ diverges. Assumption (ref) requires normalized population moments to be sufficiently smooth and consistently estimable. Assumption (ref) is our key departure from the related literature. In essence, it requires that the correlation matrix of the moment functions corresponding to close-to-binding moment conditions has eigenvalues uniformly bounded from below.\footnote{Assumption (ref) allows for high correlation among moment inequalities that cannot cross. This covers equality constraints but also entry games as the ones studied in CilibertoTamer09.} Under this condition, we are able to show that in the limit problem corresponding to (ref) --where constraints are replaced with their local linearization using population gradients and Gaussian processes-- the probability of coverage increases continuously in $c$. If such continuity is directly assumed (Assumption (ref)), Theorem (ref) remains valid (Online Appendix (ref)). While the high level Assumption (ref) is similar in spirit to a key condition (Assumption A.2) in BCS, we propose Assumption (ref) due to its familiarity and ease of interpretation; a similar condition is required for uniform validity of standard point identified Generalized Method of Moments inference. In Online Appendix (ref) we verify that our assumptions hold in some of the canonical examples in the partial identification literature: mean with missing data, linear regression and best linear prediction with interval data (and discrete covariates), entry games with multiple equilibria (and discrete covariates), and semi-parametric binary regression models with discrete or interval valued covariates mag:mau08.

Assumptions (ref)-(ref) define the class of DGPs over which our proposed method yields uniformly asymptotically valid coverage. This class is non-nested with the class of DGPs over which the profiling-based methods of Romano_Shaikh2008aJSPIWP and BCS are uniformly asymptotically valid. KMS_2017 show that in well behaved cases, calibrated projection and BCS-profiling are asymptotically equivalent. They also provide conditions under which calibrated projection has lower probability of false coverage in finite sample, thereby establishing that the two methods' finite sample power properties are non-ranked.

Convergence of the E-A-M Algorithm

We next provide formal conditions under which the sequence $p'\theta^{*,L}$ generated by the E-A-M algorithm converges to the true end point of $CI_n$ as $L\to\infty$ at a rate that we obtain. Although $p'\theta^{*,L}=\max\{p'\theta^{(\ell)}:\ell \in\{1,..., L\}, \bar g(\theta)\le c(\theta^{(\ell)})\}$, so that $\theta^{*,L}$ satisfies the {\em true\/} constraints for each $L$, the sequence of evaluation points $\theta^{(\ell)}$ is mostly obtained through expected improvement maximization (M-Step) with respect to the {\em approximating\/} surface $c_L(\cdot)$. Because of this, a requirement for convergence is that the function $c(\cdot)$ is sufficiently smooth, so that the approximation error in $|c(\theta)-c_L(\theta)|$ vanishes uniformly in $\theta$ as $L\to\infty$.\footnote{As in Bull_Convergence_2011, our convergence result accounts for the fact that the parameters of the Gaussian process prior in (ref) are re-estimated for each iteration of the A-step using the “training data" $\{\theta^\ell,c(\theta^\ell)\}_{\ell=1}^L$.} We furthermore assume that the constraint set in (ref) satisfies a degeneracy condition introduced to the partial identification literature by CHT.\footnote{CHT impose the condition on the population identified set.} In our application, the condition requires that $\mathcal{C}_n(\hat{c}_n)$ has an interior and that the inequalities in (ref), when evaluated at points in a (small) $\tau$-contraction of $\mathcal{C}_n(\hat{c}_n)$, are satisfied with a slack that is proportional to $\tau$. Theorem (ref) below establishes that these conditions jointly ensure convergence of the E-A-M algorithm at a specific rate. This is a novel contribution to the literature on response surface methods for constrained optimization.

In the formal statement below, the expectation $E_{\mathbb Q}$ is taken with respect to the law of $(\theta^{(1)},...,\theta^{(L)})$ determined by the Initialization step and the M-step but conditioning on the sample. We refer to Appendix (ref) for a precise definition of $E_{\mathbb Q}$ and a proof of the theorem.

The requirement that $\Theta$ is a compact hyperrectangle with nonempty interior can be replaced by a requirement that $\Theta$ belongs to the interior of a closed hyperrectangle in $\mathbb{R}^d$. Assumption (ref) specifies the types of kernel to be used to define the correlation functional in (ref). Assumption (ref) collects requirements on differentiability of $g_j(\theta),j=1,\dots,J$, and smoothness of $c(\theta)$. Assumption (ref) is the degeneracy condition discussed above.

To apply Theorem (ref) to calibrated projection, we provide low level conditions (Assumption (ref) in Online Appendix (ref)) under which the map $\theta \mapsto \hat c_n(\theta)$ uniformly stochastically satisfies a Lipschitz-type condition. To get smoothness, we work with a mollified version of $\hat{c}_n$, denoted $\hat c_{n,\tau_n}$ in equation (ref), where $\tau_n=o(n^{-1/2})$.\footnote{For a discussion of mollification, see e.g. Rockafellar_Wets2005aBK.} Theorem (ref) in the Online Appendix shows that $\hat c_n$ and $\hat c_{n,\tau_n}$ can be made uniformly arbitrarily close, and that $\hat c_{n,\tau_n}$ yields valid inference as in (ref). In practice, we directly apply the E-A-M steps to $\hat c_n$.

The key condition imposed in Theorem (ref) is Assumption (ref). It requires that the GMS function used is Lipschitz in its argument,\footnote{This requirement rules out the GMS function in footnote (ref), but it is satisfied by other GMS functions proposed by AS.} and that the standardized moment functions are Lipschitz in $\theta$. In Online Appendix (ref) we establish that the latter condition is satisfied by some canonical examples in the moment (in)equality literature: mean with missing data, linear regression and best linear prediction with interval data (and discrete covariates), entry games with multiple equilibria (and discrete covariates), and semi-parametric binary regression models with discrete or interval valued covariates mag:mau08.\footnote{For these same examples we verify the differentiability requirement in Assumption (ref) on $g_j(\theta)$.}

The E-A-M algorithm is proposed as a method to implement our statistical procedure, not as part of the statistical procedure itself. As such, its approximation error is not taken into account in Theorem (ref). Our comparisons of the confidence intervals obtained through the use of E-A-M as opposed to directly solving problems (ref) through the use of MATLAB's fmincon in our empirical application in the next section suggest that such error is minimal.

Empirical Illustration: Estimating a Binary Game

We employ our method to revisit the study in KT15 of “what explains the decision of an airline to provide service between two airports." We use their data and model specification.\footnote{The data, which pertains to the second quarter of the year 2010, is downloaded from \url{http://qeconomics.org/ojs/index.php/qe/article/downloadSuppFile/371/1173}.} Here we briefly summarize the set-up and refer to KT15 for a richer discussion.

The study examines entry decisions of two types of firms, namely Low Cost Carriers ($LCC$) versus Other Airlines ($OA$). A market is defined as a trip between two airports, irrespective of intermediate stops. The entry decision $Y_{\ell,i}$ of player $\ell\in\{LCC,OA\}$ in market $i$ is recorded as a $1$ if a firm of type $\ell$ serves market $i$ and $0$ otherwise. Firm $\ell$'s payoff equals $Y_{\ell,i}(Z_{\ell,i}'\vartheta_\ell+\delta_i Y_{-\ell,i}+u_{\ell,i})$, where $Y_{-\ell,i}$ is the opponent's entry decision. Each firm enters if doing so generates non-negative payoffs. The observable covariates in the vector $Z_{\ell,i}$ include the constant and the variables $W_i^{size}$ and $W_{\ell,i}^{pres}$. The former is market size, a market-specific variable common to all airlines in that market and defined as the population at the endpoints of the trip. The latter is a firm-and-market-specific variable measuring the market presence of firms of type $\ell$ in market $i$ KT15. While $W_i^{size}$ enters the payoff function of both firms, $W_{LCC,i}^{pres}$ (respectively, $W_{OA,i}^{pres}$) is excluded from the payoff of firm $OA$ (respectively, $LCC$). Each of market size and of the two market presence variables are transformed into binary variables based on whether they realized above or below their respective median. This leads to a total of 8 market types, hence $J_1=16$ moment inequalities and $J_2=16$ moment equalities. The unobserved payoff shifters $u_{\ell,i}$ are assumed to be i.i.d. across $i$ and to have a bivariate normal distribution with $E(u_{\ell,i})=0$, $Var(u_{\ell,i})=1$, and $Corr(u_{LCC,i},u_{OA,i})=r$ for each $i$ and $\ell\in\{LCC,OA\}$, where the correlation $r$ is to be estimated. Following KT15, we assume that the strategic interaction parameters $\delta_{LCC}$ and $\delta_{OA}$ are negative, that $r\ge 0$, and that the researcher imposes these sign restrictions. To ensure that Assumption (ref) is satisfied,\footnote{This assumption, common in the literature on projection inference, requires that $D_{P,j}(\theta)$ are Lipschitz in $\theta$ and have bounded norm. But $\partial(\{E_P[m_j(X,\cdot)]/\sigma_{P,j}(\cdot)\})/\partial r$ includes a denominator equal to $(1-r^2)^2$. As $r\to 1$, this leads to a violation of the assumption and to numerical instability.} we furthermore assume that $r\le 0.85$ and use this value as its upper bound in the definition of the parameter space.

The results of the analysis are reported in Table (ref), which displays $95\%$ nominal confidence intervals (our $CI_n$ as defined in equations (ref)-(ref)) for each parameter. The output of the E-A-M algorithm is displayed in the accordingly labeled column. The next column shows a robustness check, namely the output of MATLAB's fmincon function, henceforth labelled “direct search," that was started at each of a widely spaced set of feasible points that were previously discovered by the E-A-M algorithm. We emphasize that this is a robustness or accuracy check, not a horse race: Direct search mechanically improves on E-A-M because it starts (among other points) at the point reported by E-A-M as optimal feasible. Using the standard MultiStart function in MATLAB instead of the points discovered by E-A-M produces unreliable and extremely slow results. In 10 out of 18 optimization problems that we solved, the E-A-M algorithm's solution came within its set tolerance ($0.005$) from the direct search solution. The other optimization problems were solved by E-A-M with a minimal error of less than $5\%$.

Table (ref) also reports computational time of the E-A-M algorithm, of the subsequent direct search, and the total time used to compute the confidence intervals. The direct search greatly increases computation time with small or negligible benefit. Also, computational time varied substantially across components. We suspect this might be due to the shape of the level sets of $\max_{j=1,\dots,J}\sqrt n \bar m_{n,j}(\theta)/\hat\sigma_{n,j}(\theta)$: By manually searching around the optimal values of the program, we verified that the level sets in specific directions can be extremely thin, rendering search more challenging.

Comparing our findings with those in KT15, we see that the results qualitatively agree. The confidence intervals for the interaction effects ($\delta_{LCC}$ and $\delta_{OA}$) and for the effect of market size on payoffs ($\vartheta^{size}_{LCC}$ and $\vartheta^{size}_{OA}$) are similar to each other across the two types of firms. The payoffs of $LCC$ firms seem to be impacted more than those of $OA$ firms by market presence. On the other hand, monopoly payoffs for $LCC$ firms seem to be smaller than for $OA$ firms.\footnote{Monopoly payoffs are those associated with a market with below-median size and below-median market presence (i.e., the constant terms).} The confidence interval on the correlation coefficient is quite large and includes our upper bound of 0.85.\footnote{Being on the boundary of the parameter space is not a problem for calibrated projection; indeed, it is accounted for in the calibration of $\hat{c}_n$ in equations (ref)-(ref).}

For most components, our confidence intervals are narrower than the corresponding 95% credible sets reported in KT15.\footnote{For the interaction parameters $\delta$, Kline and Tamer's upper confidence points are lower than ours; for the correlation coefficient $r$, their lower confidence point is higher than ours.} However, the intervals are not comparable for at least two reasons: We impose a stricter upper bound on $r$ and we aim to cover the projections of the true parameter value as opposed to the identified set.

Overall, our results suggest that in a reasonably sized, empirically interesting problem, calibrated projection yields informative confidence intervals. Furthermore, the E-A-M algorithm appears to accurately and quickly approximate solutions to complex smooth nonlinear optimization problems.

Conclusion

This paper proposes a confidence interval for linear functions of parameter vectors that are partially identified through finitely many moment (in)equalities. The extreme points of our {\em calibrated projection\/} confidence interval are obtained by minimizing and maximizing $p^\prime \theta$ subject to properly relaxed sample analogs of the moment conditions. The relaxation amount, or critical level, is computed to insure uniform asymptotic coverage of $p^\prime \theta$ rather than $\theta$ itself. Its calibration is computationally attractive because it is based on repeatedly checking feasibility of (bootstrap) linear programming problems. Computation of the extreme points of the confidence intervals is furthermore attractive thanks to an application of the response surface method for global optimization; this is a novel contribution of independent interest. Indeed, one key result is a convergence rate for this algorithm when applied to constrained optimization problems in which the objective function is easy to evaluate but the constraints are “black box" functions. The result is applicable to any instance when the researcher wants to compute confidence intervals for optimal values of constrained optimization problems. Our empirical application and Monte Carlo analysis show that, in the DGPs that we considered, calibrated projection is fast and accurate, and also that the E-A-M algorithm can greatly improve computation of other confidence intervals.

\ifx\undefined\leavevmode\rule[.5ex]{3em}{.5pt}\ \fi \ifx\undefined\textsc \let\tmpsmall\tmpsmall\sc \fi

\tmpsmall\sc

\numberwithin{equation}{section} \numberwithin{figure}{section} \numberwithin{table}{section}

\newgeometry{ letterpaper, total={210mm,297mm}, left=20mm, right=20mm, top=30mm, bottom=30mm, }

\onehalfspacing

\pagenumbering{arabic}

\setdisplayskipstretch{0.8}

\tmpsmall\sc