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.
54,491 characters · 14 sections · 71 citation commands
Fast, Robust Inference for Linear Instrumental Variables Models using Self-Normalized Moments
Instrumental variables are widely used in applied econometrics, yet computationally tractable robust inference remains challenging. It is well known that inference robust to weak instruments can be conducted by inverting robust tests. However, to our knowledge there do not exist inference methods which are simultaneously robust to weak instruments, many instruments/regressors (e.g., larger than the sample size), potentially invalid instruments (i.e., an exclusion restriction may be violated) and potentially endogenous instruments, nor which can offer coverage guarantees in small samples. Moreover, inverting robust tests can be computationally challenging andrews19 because their non-rejection regions are not convex mikusheva10 and a grid search is infeasible with just a handful of regressors (see andrews16con, supplementary material).
We address both the statistical and computational challenges. First, we provide an approach to inference based on self-normalized sample moment conditions, which we refer to as Self-Normalized Instrumental Variables (SNIV). Due to the minimal assumptions it requires, SNIV simultaneously allows for weak instruments and conditional heteroskedasticity, and can be applied equally to the standard low-dimensional setting (large sample size, few regressors, few instruments) and to the high-dimensional setting in which the number of instruments and/or regressors can be large, possibly much larger than the sample size. For example, in a model with a single endogenous regressor, SNIV could be used to construct a confidence interval which is simultaneously robust to many and to weak instruments.\footnote{To further motivate our framework, there are several reasons to expect models with multiple endogenous regressors to become increasingly popular. Potential applications include demand systems with many goods and endogenous expenditure gautier21 or prices belloni22, as well as models of peer effects with unknown peer relationships rose18. More broadly, due to increasing availability of rich datasets and the potential to allow for heterogeneous treatment effects by using interactions of the treatment with individual characteristics, applied research has recently considered models with multiple exogenous variables of interest (e.g., farrell20). A natural extensionis to allow for endogenous treatment belloni22.} We extend SNIV to settings in which one or more instrument may be `invalid' (i.e., the exclusion restriction fails, see kolesar15,kang16) or endogenous without requiring a pre-test, and propose the use of an a-priori upper bound on the number of invalid/endogenous instruments for settings in which the set of identifiable parameters would otherwise be unbounded.
Second, we provide a computational implementation of SNIV which we show can also be applied to rapidly invert other robust tests. This is because test inversion can typically be cast as a semi-algebraic optimization problem, hence we can apply methods from the literature on semi-algebraic optimization. The basic idea to deal with computational intractability is to attempt to solve a hiearchy of semidefinite optimization problems. Solving each optimization problem delivers an outer bound on the confidence region. As we proceed up the hierarchy, the bounds become sharper at the expense of greater computational burden. This allows the researcher to effectively trade off sharpness with available computational resources. In practice, the bounds obtained towards the beginning of the hierarchy are often exact. A simple diagnostic informs the researcher if exact bounds have been attained. Similar computational approaches have been applied by gautier18, lee20 and auerbach22.
In contrast to a grid search, our approach can be applied to settings with multiple regressors. In contrast to heuristic/local optimization methods, our approach guarantees an outer bound on the confidence region, hence does not risk compromising the coverage guarantee. We illustrate how our computational approach can be used to rapidly invert existing robust tests by combining it with the results of guggenberger12, guggenberger19 and guggenberger21 to obtain weak instrument robust Anderson-Rubin (AR) confidence intervals which can be computed near instantaneously, even when a grid search is infeasible. We also show that our approach can be applied to rapidly invert other robust tests such as the Lagrange-multiplier (LM) test kleibergen02,moreira02 and the Conditional Likelihood Ratio (CLR) test moreira03.
We conduct a Monte-Carlo experiment in which we demonstrate SNIV and AR confidence regions are both easily implemented in a setting in which a grid search is computationally intractable. SNIV has similar coverage to the AR test in designs with either strong or weak instruments. With many instruments, SNIV maintains coverage close to the nominal level but AR does not. We also show that SNIV can be applied to conduct informative inference with invalid instruments and endogenous instruments in challenging designs, to which existing approaches cannot be applied.
Our work is related to the literature on many instruments (e.g., bekker94,angrist99,donald01,anderson05,chao05,stock05,hansen08,ackerberg09,van10,chao12,hausman12,anatolyev13,hansen14,kolesar18; see anatolyev19 for a recent review) and weak instruments (e.g., anderson49,kleibergen02,moreira02,moreira03,mikusheva10,guggenberger12,andrews16con,guggenberger19,guggenberger21; see andrews19 for a recent review), but allows simultaneously for weak instruments and for the number of regressors and/or instruments to be large, possibly much larger than the sample size. This is because we conduct inference based on moderate deviations of self-normalized sample moments (e.g., pinelis94,bertail08,jing03) instead of using a Central Limit Theorem.
The most closely related papers are gautier11, belloni12, gold20, gautier21 and belloni22, all of which also consider the linear instrumental variables model in a potentially high-dimensional setting. gautier11 and gautier21 suggest to combine a point estimator with lower bounds on its sensitivity characteristics to perform robust inference, but their confidence region is larger than ours and the authors do not provide a disciplined way to trade off computational complexity and sharpness when implementing their approach. belloni12 discuss inference based on self-normalization, but do not propose a practical computational solution, nor allow for invalid/endogenous instruments. gold20, gautier21 and belloni22 propose confidence regions for a subset of parameters of interest (e.g., confidence intervals) but these rely on stronger assumptions than we use below. For example, these papers propose methods which are not robust to weak instruments and do not allow for potentially invalid nor endogenous instruments. We view our work as complementary to gold20, gautier21 and belloni22, providing the applied researcher with a more robust alternative, but one which may sometimes yield wider bounds in practice.
Finally, our extensions of SNIV are related to the literature on invalid and endogenous instruments. Regarding invalid instruments (e.g., kolesar15,kang16), we allow for a setting with multiple endogenous regressors, potentially weak instruments, and for the number of instruments and/or regressors to be larger than the sample size. Regarding endogenous instruments (e.g., sargan58,hansen82,anatolyev11,lee12,chao14), we do not perform specification tests, but instead perform inference directly, accounting for potential endogeneity of an unknown subset of instruments.
We proceed as follows. Section (ref) sets out our model. Section (ref) defines SNIV, establishes its coverage guarantee, and provides extensions to potentially invalid and endogenous instruments. Section (ref) presents our computational method, which we apply to invert existing robust tests in Section (ref). Section (ref) presents a Monte-Carlo experiment and Section (ref) concludes. All proofs are gathered in the appendix.
To simplify the exposition we consider an i.i.d. sample of size $n$. The i.i.d. setting is not critical for our results and can be relaxed by using an alternative choice of $r_n$ below, for which we provide appropriate references. The population model comprises an outcome $Y$, regressors $X\in{\mathbb R}^{d_X}$, and instrumental variables $Z\in{\mathbb R}^{d_Z}$ of joint distribution $\mathbb{P}$. ${\mathbb E}$ is the expectation under $\mathbb{P}$ and ${\mathbb E}_n$ is its sample counterpart. Our results apply to a sequence of models indexed by $n$. For simplicity of exposition we do not make this explicit, but we occasionally note that certain objects can depend on $n$. To allow for high-dimensional data, the relative magnitudes of $n,d_X$ and $d_Z$ are unrestricted, and both $d_X$ and $d_Z$ can grow with $n$. For $b\in{\mathbb R}^{d_X}$, $U(b)\triangleq Y-X^{\top}b$, $\mathbb{P}(b)$ is the distribution of $\left(X,Z,U(b)\right)$ implied by $\mathbb{P}$. For $S\subseteq [d]\triangleq \{1,2,...,d\}$, $|S|$ is its cardinality and $S^c$ its complement. For $\Delta\in\mathbb{R}^d$, $S(\Delta)\triangleq\{k\in[d]:\ \Delta_{k}\ne0\}$ and $|\Delta|_p$ is the $\ell^p$-norm of $\Delta$. For a polynomial $p$, deg$(p)$ is its degree. We use $\mathbf{M}\succcurlyeq \mathbf{0}$ to say that the matrix $\mathbf{M}$ is positive semidefinite.
The linear instrumental variables model is
where $\mathcal{B}\subseteq \mathbb{R}^{d_X}$ is the parameter space and $\mathcal{P}$ is a nonparametric class. The set $\mathcal{I}$ collects the vectors which satisfy (ref)-(ref). As will be made clear below, our results are for all $\beta\in\mathcal{I}$, hence for the true value $\beta^*$. We use the class $\mathcal{P}$ to permit the use of results on moderate deviations of self-normalized sums for inference. We consider four classes, including\\
{\bf Class 1.} There exists $\delta$ in $(0,1]$ and $\mu_{2+\delta}>0$ such that $$\displaystyle\left|\left( \left({\mathbb E}\left[|Z_lU(\beta)|^{2+\delta}\right]\right) \left({\mathbb E}\left[Z_l^2U(\beta)^2\right]\right)^{-(2+\delta)/2}\right)_{l\in[{d_Z}]}\right|_{\infty} \le\mu_{2+\delta},$$ and ${d_Z}\le \alpha/(2\Phi(-n^{1/2-1/(2+\delta)}\mu_{2+\delta}^{-1/(2+\delta)}))$,
where $\alpha\in(0,1)$ is a confidence level and $\Phi$ the normal CDF;\\
{\bf Class 2.} $\exists \mu_4>0:$ $\max_{l\in[{d_Z}]}\mathbb{E}[Z_l^4U(\beta)^4](\mathbb{E}[Z_l^2U(\beta)^2])^{-2}\le \mu_{4}$, ${d_Z}<\alpha\exp\left(n/\mu_4\right)/(2e+1)$ and $n-\mu_4\log({d_Z}(2e+1)/\alpha)\geq n/2$;\\
{\bf Class 3.} $Z_lU(\beta)$ is symmetric for all $l\in[d_Z]$ and ${d_Z}<9\alpha/\left(4e^3\Phi\left(-\sqrt{n}\right)\right)$.\\
Classes 1-2 require mild bounds on ratios of moments, whereas Class 3 requires no bounds but uses symmetry. Further classes allowing for dependence and non i.d. data can be found in chen16 and references therein. In Section (ref) we consider a fourth class based on Gaussian approximation rather than self-normalization.
The $1-\alpha$ SNIV confidence set is
where $\mathbf{D}(\beta)$ is the $d_Z\times d_Z$ positive, diagonal matrix with $l^{\text{th}}$ diagonal element $\mathbb{E}_n[Z_l^2U(\beta)^2]^{-1/2}$ used to self-normalize the $d_Z$ moments and $r_n$ depends on the class. Under Class 1 we set $r_n=-\Phi^{-1}\left(\alpha/(2{d_Z})\right)/\sqrt{n}$. Under Class 2 we set $r_n=2\sqrt{\log({d_Z}(2e+1)/\alpha)/n}$. Under Class 3 we set $r_n=-\Phi^{-1}(9\alpha/(4d_Ze^3))/\sqrt{n}.$
Proposition (ref) shows that the coverage of SNIV is at least the nominal level uniformly over the identifiable parameters and the distributions of the data they imply. Beyond the class used, no further assumptions are needed. Classes 1-3 allow for conditional heteroscedasticity, do not restrict the joint distribution of $X$ and $Z$ (hence are robust to weak instruments) and have very mild requirements on the relative magnitudes of $n$ and $d_Z$ (hence are robust to many regressors and/or instruments). For Class 1 the coverage guarantee is asymptotic in $n$ in such a way that $d_X$ and $d_Z$ can grow with (and be much larger than) $n$. For Classes 2-3 the coverage guarantee is for any $n$.
The SNIV confidence set collects vectors for which the $\ell_\infty$ deviation from zero of the self-normalized sample moment is at most $r_n$. The core components which deliver uniformity, finite sample validity and robustness to identification are the $\ell_\infty$-norm and self-normalization of the moments. The $\ell_\infty$-norm is crucial so as to allow for $d_Z$ larger than $n$ because it permits $r_n$ to be of the order $\log(d_Z)/\sqrt{n}$. This means that the SNIV confidence set can be small even when the number of instruments is much larger than the sample size.\footnote{For Class 3, $r_n\le2\log\left(4{d_Z}e^3/(9\alpha)\right)/\sqrt{n},$ $\forall\alpha\in[0,1],d_Z\ge1$(because $\Phi^{-1}(a)\ge2\log(a)$ if $0<a\le\exp(-1/(4\pi))$).} Note also that $d_X$ can be arbitrarily large with respect to $n$.
The SNIV confidence set may be conservative when the instruments are strongly correlated with one another because $r_n$ is based on a union bound over $d_Z$ self-normalized sample moments. gautier21 propose an alternative to self-normalization based around the multiplier bootstrap of chernozhukov13, which we implement here. We modify the SNIV confidence set by replacing $\mathbf{D}(\beta)$ by ${\mathbb E}_n[U(\beta)]^{-1/2}\mathbf{D_Z}$ and $r_n$ by the $1-\alpha$ quantile of $|\mathbf{D_Z}{\mathbb E}_n[ZW]|_\infty$ (computed by simulation) in its definition, where $\mathbf{D_Z}$ is a $d_Z\times d_Z$ diagonal matrix with $l^{th}$ diagonal element ${\mathbb E}_n[Z_l^2]^{-1/2}$ and $W$ is a standard normal random variable which is independent of $Z$. The corresponding class is\\
{\bf Class 4.} There exist constants $C$ and $c$, and $B_n$ such that, for all $(\beta,\mathbb{P})$: $\beta\in\mathcal{I}$, $U(\beta)\perp Z$; $|Z|_\infty\leq B_n$ (a.s.); ${\mathbb E}[U(\beta)^4]\leq C$; and $B_n^4\log(d_Zn)^7/n\leq Cn^{-c}$,\\
which delivers the following coverage guarantee.
Further classes based on Gaussian approximation but allowing for non i.d. and dependent data can be found in zhang17 and references therein.
In the high-dimensional setting with $d_X$ larger than $n$, a natural and commonly used restriction is that $\beta$ is sparse, meaning that it has many elements exactly equal to zero but the researcher does not know which ones. Sparsity implies that there exists an underlying parsimonious model which is unknown to the researcher. It can be used to motivate $\ell_1$ penalized estimators such as the LASSO of tibshirani96 for regression or the STIV of gautier21 for instrumental variables. In the instrumental variables context, sparsity can be interpreted as imposing exclusion restrictions of unknown locations. As explained below, this is particularly useful in the underidentified case with $d_Z<d_X$, which can arise, for example, when there is uncertainty as to which candidate instruments can be excluded, implying that some instruments may be invalid kolesar15,kang16.
The SNIV confidence set can easily accommodate sparsity. We define $S_Q\subseteq [d_X]$ as the indices of the regressors of questionable relevance (i.e., whose entry of $\beta^*$ may be zero). We denote by $d_Q\triangleq|S_Q|$ and modify $\mathcal{I}$ and $\widehat{\mathcal{C}}$ to include the restriction
where $s\in[d_Q]$ is an upper bound chosen by the researcher, and recalling that $S(\beta)\subseteq[d_X]$ is the support of $\beta$. Though we do not make it explicit, both $s$ and $S_Q$ can depend on $n$. Since the choice of $s$ provides a guarantee on the sparsity, we refer to it as a sparsity certificate. When the sparsity certificate $s$ is used, we use the notation $\mathcal{I}_s$ and $\widehat{\mathcal{C}}_s$ in place of $\mathcal{I}$ and $\widehat{\mathcal{C}}$. Clearly, $\mathcal{I}_{d_Q}=\mathcal{I}$ and $\widehat{\mathcal{C}}_{d_Q}=\widehat{\mathcal{C}}$.
Interestingly, $\mathcal{I}_s$ can be a singleton even when $\mathcal{I}$ is not. This means that sparsity can lead to point identification even in `underidentified' models (i.e., when $d_Z<d_X$). In general, $\mathcal{I}_s$ is a singleton if there is a solution for only one of the ${d_Q \choose s}$ overdetermined systems based on (ref)-(ref) and it is unique. For example, $\mathcal{I}_s$ can be a singleton when $s+d_X-d_Q<d_Z<d_X$ and sparsity implies that some exogenous regressors have a zero coefficient (i.e., they are excluded, see kang16). The basic idea is that excluded exogenous regressors can serve as instruments for included endogenous regressors, but we need not necessarily know which regressors are excluded. Finally, if $S_Q=[d_X]$, $\mathcal{I}_s$ is a singleton if all matrices formed from $2s$ columns of ${\mathbb E}[ZX^\top]$ have rank $2s$ candes07. The corresponding order condition is $s\leq d_Z/2$, which does not depend on $d_X$.
However, SNIV does not require that $\mathcal{I}_s$ be a singleton. For example, $\mathcal{I}_s$ can comprise a finite union of singletons. Figure (ref) depicts such an example with $d_Z=1$, $d_X=2$, $\mathcal{B}={\mathbb R}^{d_X}$, $d_Q=[d_X]$ and $s=1$, in which case $\mathcal{I}_s$ is the intersection of the line ${\mathbb E}[Zy]={\mathbb E}[ZX^\top]\beta$ with the set $\{\beta\in{\mathbb R}^2:\beta_1=0\text{ or } \beta_2=0\}$. The SNIV confidence set allows for such partially identified cases due to the uniformity over $s$ and $\mathcal{I}_s$ in the coverage guarantee, which is obtained by replacing $\inf_{(\beta,\mathbb{P}):\beta\in\mathcal{I}}\mathbb{P}(\beta\in\widehat{\mathcal{C}})$ by $\min_{s\in[d_Q]}\inf_{(\beta,\mathbb{P}):\beta\in\mathcal{I}_s}\mathbb{P}(\beta\in\widehat{\mathcal{C}}_s)$ in Propositions (ref) and (ref).
Testing instrument exogeneity is a classical problem to which our framework can be applied. Introducing $\theta\in{\mathbb R}^{d_Z}$ to account for the possible failure of exogeneity, we replace (ref)-(ref) by
where $\theta_l\ne0$ means that $Z_l$ is endogenous, $\mathbb{P}\left(b,t\right)$ is the distribution of $\left(X,Z,ZU(b)-t\right)$ implied by $\mathbb{P}$ and $\Theta\subseteq{\mathbb R}^{d_Z}$ encodes restrictions on $\theta$. For example, $\Theta$ may be such that the sign of the correlation of a regressor and the structural error is known. An important restriction encoded by $\Theta$ is $\theta_{ S_{\perp}}=0$ for $S_{\perp}\subseteq[d_Z]$, which indexes the instruments known to be exogenous. The remaining instruments are potentially endogenous. We can use a sparsity certificate to place an upper bound on the number of endogenous instruments, given by
for a given $\widetilde{s}\in[\widetilde{d}_Q]$, where $\widetilde{d}_Q\triangleq d_Z-|S_\perp|$. Thus, though the identities of the endogenous instruments may not be known, their number can be restricted. The counterpart of $\mathcal{I}_s$, denoted by $\mathcal{I}_{s,\widetilde{s}}$, collects the vectors which satisfy (ref)-(ref) and the sparsity restrictions in (ref) and (ref). Under Classes 1-3,\footnote{Class 4 is not applicable with possibly endogenous instruments.} SNIV is
where $\mathbf{D}(\beta,\theta)$ is the $d_Z\times d_Z$ positive, diagonal matrix with $l^{\text{th}}$ diagonal element $\mathbb{E}_n[(Z_lU(\beta)-\theta_l)^2]^{-1/2}$. This set allows one to simutaneously perform inference on $\beta^*$ and $\theta^*$, without requiring, for example, a pilot estimator and subsequent test of instrument exogeneity. The coverage guarantee is obtained by replacing $\inf_{(\beta,\mathbb{P}):\beta\in\mathcal{I}}\mathbb{P}(\beta\in\widehat{\mathcal{C}})$ by $\min_{\widetilde{s}\in [\widetilde{d}_Q]}\min_{s\in[d_Q]}\inf_{(\beta,\theta,\mathbb{P}):(\beta,\theta)\in\mathcal{I}_{s,\widetilde{s}}}\mathbb{P}((\beta,\theta)\in\widehat{\mathcal{C}}_{s,\widetilde{s}})$ in Proposition (ref). Both $\widetilde{S}_Q$ and $\widetilde{s}$ can depend on $n$.
To implement SNIV the researcher needs some way to summarize the vectors which lie in the confidence set. belloni12 propose to use a grid for a confidence set with no sparsity constraints nor potentially endogenous instruments. This involves checking whether the inequalities in the definition of the SNIV confidence set are verified for every $\beta$ on a grid over $\mathcal{B}$, and is a practical solution when $d_X$ is small. However, a grid search quickly becomes infeasible for moderate $d_X$. In a low-dimensional setting, one can first partial-out a small number of exogenous regressors (see Remark (ref)) so that $d_X$ is the number of endogenous regressors, which may be sufficiently small so as to use a grid. Otherwise we require an alternative.
We propose a method based on solving convex optimization problems. For a given direction $\mathbf{u}\in{\mathbb R}^{d_X+d_Z}$ normalized to satisfy $|\mathbf{u}|_2=1$ and the function $f_{\mathbf{u}}(\beta,\theta)\triangleq\mathbf{u}^\top(\beta^\top,\theta^\top)^\top$ we seek to compute
which is the support function of $\widehat{\mathcal{C}}_{s,\widetilde{s}}$. By solving (ref) for all directions $\mathbf{u}\in\{\mathbf{u}\in{\mathbb R}^{d_X+d_Z}:|\mathbf{u}|_2=1\}$, we obtain the convex envelope of $\widehat{\mathcal{C}}_{s,\widetilde{s}}$ defined by the inequalities $f_\mathbf{u}(\beta,\theta)\geq f^*(\mathbf{u})$ for all $\mathbf{u}\in\{\mathbf{u}\in{\mathbb R}^{d_X+d_Z}:|\mathbf{u}|_2=1\}$.\footnote{In practice, we consider only a finite number of directions.}
If $\widehat{\mathcal{C}}_{s,\widetilde{s}}$ is convex, solving (ref) is straightforward. In general $\widehat{\mathcal{C}}_{s,\widetilde{s}}$ is not convex because none of the inequalities in its definition define a convex set. This is unavoidable because $\mathcal{I}_{s,\widetilde{s}}$ need not be convex. We now show that $\widehat{\mathcal{C}}_{s,\widetilde{s}}$ is a semi-algebraic set (i.e., a set defined by polynomial inequalities) and apply methods in semi-algebraic optimization to problem (ref).
The requirement that $\mathcal{B}$ and $\Theta$ are semi-algebraic is mild. For example, $\mathcal{B}=\mathbb{R}^{d_x}$ is semi-algebraic. The additional parameter $\gamma$ is required to model the sparsity constraints in (ref) and (ref). Proposition (ref) implies that
which can be computed by solving a polynomial optimization problem. Due to non-convexity, exact computation of $f^*(\mathbf{u})$ is NP-hard. Instead, we focus on solving convex relaxations of (ref). Convex relaxation is routinely used to construct computationally tractable estimators. For example, LASSO uses an $\ell_1$ penalty as a convex relaxation of a sparsity constraint such as (ref). We solve a sequence of convex relaxations, delivering a hierarchy of convex optimization problems. Following the seminal paper of lasserre01, such hierarchies have attracted much attention in the optimization literature in recent years. We first provide a general summary of the approach, then explain the specific hierarchy we propose.
The most important feature of a hierarchy is that it is disciplined, meaning that it delivers a monotone sequence of lower bounds converging to $f^*(\mathbf{u})$. If $f^*_h(\mathbf{u})$ is the optimal value obtained by solving the $h^{\text{th}}$ convex optimization problem in the hierarchy, we have $f^*_{h}(\mathbf{u})\leq f_{h+1}^*(\mathbf{u})\leq f^*(\mathbf{u})$ for all $h\in\mathbb{N}$ and $f^*_h(\mathbf{u})\to f^*(\mathbf{u})$ as $h\to\infty$. As $h$ increases, though convex, the optimization problems become more computationally intensive. It is also generically the case that there exists finite $h^*$ such that $f^*_{h^*}(\mathbf{u})=f^*(\mathbf{u})$, and that the researcher can identify when such $h^*$ has been encountered.
Monotonicity of the sequence of lower bounds on $f^*(\mathbf{u})$ is crucial. This is because it allows us to construct bounds on the convex envelope of the SNIV confidence set defined by the linear inequalities $f_\mathbf{u}(\beta,\theta)\geq f^*_{h}(\mathbf{u})$ for all $\mathbf{u}\in\mathcal{U}$ and some $h\in\mathbb{N}$, where $\mathcal{U}$ is a finite collection of directions. Since we construct a superset, the coverage guarantee cannot fall below $1-\alpha$. The larger is $h$, the closer the superset becomes to the SNIV confidence set. Hence, by varying $h$ and $\mathcal{U}$, we can trade off the computational burden with the quality of the approximation without compromising the coverage guarantee. Such a trade-off cannot be achieved by local or heuristic optimization methods nor by adjusting the spacing of a grid, both of which may compromise the coverage guarantee. Figure (ref) illustrates the SNIV confidence set and its outer approximations for a partially identified model with $d_Z=1$, $d_X=2$, $\mathcal{B}={\mathbb R}^{d_X}$, $d_Q=[d_X]$, $s=1$ and all instruments known to be exogenous.
The method we propose can also be used to compute bounds on a polynomial function of interest $p(\beta^*,\theta^*)$. For example, to obtain a confidence interval for $\beta_1^*$, we can solve (ref) for $f_{\mathbf{u}_1}$ and $f_{\mathbf{u}_2}^\top$, where $\mathbf{u}_1=(1,0,...,0)^\top$ and $\mathbf{u}_2=-\mathbf{u}_1$ (i.e., use the projection method). More generally one can consider a vector $\mathbf{p}$ of functions of interest. If the dimension of $\mathbf{p}$ is small relative to $d_X$, it is well known that the projection method can be conservative. A leading case with $d_X={\rm dim}(\mathbf{p})=1$ is a confidence interval in a model with one endogenous regressor. In this case, the projection method is not conservative and SNIV is simultaneously robust to weak instruments and to $d_Z$ much larger than $n$.
Since it is straightforward to implement, we present our application of the seminal hierarchy first proposed by lasserre01. This hierarchy is sufficiently computationally tractable to deal with problems of size likely to be encountered in empirical work. Recent advances allowing for even larger problems are provided by lasserre17 and weisser18.
To simplify the exposition, we denote the decision variable in problem (ref) by $\delta\triangleq(\beta^\top,\theta^\top,\gamma^\top)^\top$ of size $d_\delta\triangleq d_X+d_Z+d_Q+\widetilde{d}_Q$. The hierarchy uses the decision variable $\mu$, each entry of which represents a monomial of $\delta$. For example, if $d_\delta=2$ then $\mu=(1,\delta_1,\delta_2,\delta_1^2,\delta_1\delta_2,\delta_2^2,...)^\top$, so the polynomial $p(\delta)=\delta_2+2\delta_1^2$ is equivalently expressed as $\mu_3+2\mu_4$. This allows us to define the Riesz linear functional of $p$ as $L_\mu(p)=\mu_3+2\mu_4$.\footnote{We present the case in which he order of the entries of $\mu$ is a graded lexicographic order. Other orderings are possible. None of our results depend on the ordering used.} Now let the vector $\mathbf{m}_e(\delta)$ comprise all monomials of $\delta$ of degree no larger than $e$. For example, $\mathbf{m}_1(\delta)=(1,\delta_1,\delta_2)^\top$. Then we can define the moment matrix $\mathbf{M}_e(\mu)\triangleq L_\mu(\mathbf{m}_e(\delta)\mathbf{m}_e(\delta)^\top)$. For example, if $d_\delta=2$ and $e=1$, we have
Given another polynomial $q(\delta)$, we can similarly define the localizing matrix $\mathbf{M}_e(q \mu)\triangleq L_\mu(q(\delta)\mathbf{m}_e(\delta)\mathbf{m}_e(\delta)^\top)$.
At level $h$ of the hierachy we solve the semidefinite program
where $e_j$ is the smallest integer which is at least as large as deg$(\widehat{\mathbf{g}}_j)/2$ for $j\in[d_g]$. This program has a linear objective function and $d_g+1$ semidefinite constraints. The semidefinite constraint on $\mathbf{M}_h(\mu)$ arises because $\mathbf{m}_h(\delta)\mathbf{m}_h(\delta)^\top$ has rank 1. In principle we would like to impose that $\mathbf{M}_h(\mu)$ has rank 1. However, the set of rank 1 matrices is not convex. To obtain a convex problem, we use instead the set of positive semidefinite matrices. The intuition is the same for the other $d_g$ semidefinite constraints because the polynomials $\widehat{\mathbf{g}}$ are restricted to be non-negative.
Corollary (ref) follows from Proposition (ref) due to Theorem 4.2 of lasserre01. The only assumption beyond the class $\mathcal{P}$ is a technical assumption requiring that the parameter space be compact. Compactness is useful because it allows us to find $B$ sufficiently large such that the redundant polynomial constraint
holds. In practice, we augment the constraints $\widehat{\mathbf{g}}(\beta,\theta,\gamma)\geq 0$ to include (ref) prior to applying the semidefinite hierarchy.
Though compactness is a common technical assumption, in practice we may often not have compact $\mathcal{B}$ and $\Theta$. For example, we may have $\mathcal{B}={\mathbb R}^{d_X}$. In this case we suggest increasing $B$ until (ref) ceases to bind at the solution. If the SNIV confidence set is unbounded in direction $\mathbf{u}$, (ref) will always bind. In practice this is of little consequence since there is little distinction between $f^*(\mathbf{u})$ being $-\infty$ or an arbitrarily small finite constant. Thus, when the parameter space is not compact, our approach characterizes the intersection of the SNIV confidence set with an arbitrarily large $\ell_2$ ball.
The intuition for the result that $f^*_{h}(\mathbf{u}) \leq f^*(\mathbf{u})$ for all $h\in\mathbb{N}$ comes from convex relaxation. By replacing rank 1 constraints for the moment and localizing matrices by positive semidefinite constraints, we minimize over a larger set, hence it must be that we obtain a lower bound on the optimal value. The intuition for $f^*_{h}(\mathbf{u}) \leq f_{h+1}(\mathbf{u})$ for all $h\in\mathbb{N}$ is that increasing $h$ reduces the size of the set over which we minimize, hence must always deliver a larger optimal value. The computational trade-off is also clear from the form of problem (ref) because the dimension of the moment matrix $\mathbf{M}_h(\mu)$ is ${d_\delta+h \choose h}$, which is increasing in $h$. Similarly, the dimensions of the localizing matrices are combinatorically increasing in $h$. Thus, increasing $h$ delivers a tigher bound but at increased computational cost.
We implement the hierarchy using the following algorithm proposed by lasserre15.
The basic idea is to begin with the most computationally tractable semidefinite program with $h=1$ and continue to increase $h$ until either we know that $f^*_h(\mathbf{u})=f^*(\mathbf{u})$ (step 3) or we hit the largest computationally feasible level of the hierarchy ($\overline{h}$). In practice, $\overline{h}$ is determined by the size of the problem and the available computational resources. In our Monte-Carlo experiment we use $\overline{h}=2$ on a standard desktop machine. Step 3 provides a stopping criterion which can be used to establish finite convergence of the hierarchy. For brevity, we do not provide technical conditions under which finite convergence is possible, which can be found in lasserre15 (see Theorem 6.5) and involve standard Karusch-Kuhn-Tucker conditions for an optimal solution to be a local minimizer of a nonlinear program. In fact, these conditions imply that finite convergence is achieved generically lasserre15 (see Theorem 7.6), though there is no guarantee that it is achieved for small values of $h$. In our Monte-Carlo experiment we achieve finite convergence with high frequency in some designs but with low frequency in others.
In the standard linear instrumental variables setting with large $n$ and fixed $d_Z\geq d_X$, robust inference can be conducted by inverting robust tests andrews19. In this section we show that the computational approach of Section (ref) can be applied to do so. There are myriad such tests, including but not limited to, the Anderson Rubin (AR) test anderson49, Lagrange-multiplier (LM) test kleibergen02,moreira02 and the Conditional Likelihood Ratio (CLR) test moreira03. All of these tests have a non-rejection region (i.e., a confidence set) of the form
where $\widehat{p}$ and $\widehat q_\alpha$ are polynomials and the coefficients of $\widehat q_\alpha$ depend on the confidence level $\alpha$. As with the SNIV confidence set, $\widetilde{\mathcal{C}}$ is semi-algebraic whenever $\mathcal{B}$ is. The polynomial inequality in the definition of $\widetilde{\mathcal{C}}$ can be degree 2 (AR test, CLR test with $d_X=1$) or larger (LM test, CLR test with $d_X>1$). For example, the AR test under homoskedasticity uses $\widehat{p}(\beta)=\widehat{p}_{AR}(\beta)$ and $\widehat q_\alpha(\beta)= \widehat q_{AR,\alpha}(\beta)$, where
$\mathbf{y}$ is the $n\times 1$ vector of outcomes, $\mathbf{X}$ is the $n\times d_X$ matrix of regressors, $\mathbf{U}(\beta)\triangleq\mathbf{y}-\mathbf{X}\beta$, $\mathbf{Z}$ is the $n\times d_Z$ matrix of instruments and $C_{\alpha}(d)$ is the $1-\alpha$ quantile of the $\chi^2_{d}$ distribution. None of the above tests yield a convex confidence set mikusheva10, making a grid search computationally demanding (see andrews16con, supplementary material). Alternatives to a grid search (e.g., mikusheva10 for the CLR test) can also be computationally intensive for moderate $d_X$. Hierarchies of semidefinite optimization problems provide a practical alternative.
Sometimes the object of interest may be a function of $\beta^*$ of dimension smaller than $d_X$ (see Section (ref)). For example, one may be interested in a sub-vector of $\beta^*$. To obtain coverage probability $1-\alpha$ for a confidence set for a sub-vector (e.g., a confidence interval), we need to adjust $\widehat{q}_\alpha$. guggenberger12, guggenberger19 and guggenberger21 provide appropriate adjustments for the AR test. Decomposing $\beta=(\beta_1^\top,\beta_2^\top)^\top$, the results of guggenberger12 imply that an asymptotic $1-\alpha$ confidence set for $\beta_1^*$ under homoskedasticity is
where $d_{X_1}$ is the dimension of $\beta_1$.\footnote{This decomposition is without loss of generality because the regressors can be reordered.} Thus, for a given direction $\mathbf{u}_1\in{\mathbb R}^{d_{X_1}}$ normalized to have $|\mathbf{u}_1|_2=1$ we seek to compute $\inf_{\beta_1\in\widetilde{\mathcal{C}}_{AR}}\mathbf{u}_1^\top\beta_1$, equivalently expressed as
where $\mathbf{u}=(\mathbf{u}^\top_1,\mathbf{0}^\top)^\top$ is $d_X\times 1$. This is a polynomial optimization problem, hence we can apply convex hierarchies to find a monotonic sequence of lower bounds. For the special case of a confidence interval we have $d_{X_1}=1$ hence only need consider $\mathbf{u}_1=\pm 1$. guggenberger19 replace $ C_{\alpha}(d_Z-d_X+d_{X_1})$ with an alternative which delivers a less conservative confidence set in a finite sample and guggenberger21 extend the approach to allow for conditional heterokskedasticity. Thus, we can apply hierarchies of semidefinite optimization problems to (ref) in order to obtain AR confidence intervals which can be rapidly computed.
To illustrate our approach, we consider a setting with $d_X=10$ endogenous regressors. We choose this design because $d_X$ is large enough to render a grid search infeasible yet small enough to permit many replications of our experiment on a standard desktop machine within a reasonable timeframe.
We consider an i.i.d. sample of size $n=2000$ satisfying (ref)-(ref). The instruments are related to the regressors according to $\mathbb{E}[ZV(\Pi)]=0$ where $V(\Pi)\triangleq X-\Pi Z$ and $\Pi$ is $d_X\times d_Z$. We set $\beta^*=(1,-1,0,...,0)^\top$ and vary $\Pi^*$ by design, as explained below. The instruments follow $\mathcal{N}(0,I_{d_Z})$ and the error terms verify $(U(\beta^*),V(\Pi^*)^\top)^\top\sim\mathcal{N}(0,\Omega)$, where $\Omega_{11}=1$ (homoskedasticity) or $\Omega_{11}=Z_1^2$ (conditional heteroskedasticity), $\Omega_{1j}=(-1)^j(1-\pi^{*})/5$, $\Omega_{jj}=1-\pi^{*}$ for $j>1$ and all other entries are equal to zero. The parameter $\pi^*\in[0,1]$ determines the fraction of the variance of each regressor which is due to the instruments. In all designs, the variances of each regressor and the structural error are equal to 1.
To compute the SNIV confidence set we choose $r_n$ using Class 1 and Class 3 with $\alpha=0.05$. To implement the hierarchy, we use $\overline{h}=2$ and $B=1000$. We compute the coverage probability for the SNIV confidence set, and, for designs in which it is feasible, the AR confidence set and confidence intervals. For the SNIV and AR confidence sets, we also report the coverage probability for their outer approximations obtained by solving hierarchies of semidefinite optimization problems, defined by $f_\mathbf{u}(\beta)\geq f^*_h(\mathbf{u})$ for all $\mathbf{u}\in\mathcal{U}$, where, for $\mathcal{U}$ we use a grid of 1600 points over the surface of an $\ell_2$ ball of radius 1.\footnote{In practice we parallelize over $\mathbf{u}\in\mathcal{U}$.}
Our results are collected in Table (ref). We focus the discussion of SNIV on Class 1, which, identically to AR, provides an asymptotic coverage guarantee. Class 3 provides a finite guarantee, hence a larger confidence set in all designs. Nevertheless, we find that whenever Class 1 provides an informative confidence set, so does Class 3.\% Class 3 is also slower to compute, though still sufficiently fast to conduct informative inference using a standard machine. \\
{Classical design.} We set $d_Z=d_X=10$, $\pi^*=0.3$ and $\Pi^*=\sqrt{\pi^*}I_{d_Z}$ and do not impose any sparsity constraint. The SNIV and AR confidence sets have similar coverage, both of which are marginally below the nominal level. The SNIV confidence set is marginally narrower than the AR confidence set. The coverage of the AR confidence intervals are almost exactly equal to the nominal level, and their width is narrower than either of the confidence sets, as expected. Almost all optimization problems solved yielded an exact global optimum (i.e., $f^*_h(\mathbf{u})=f^*(\mathbf{u})$). The time taken to solve an optimization problem is a little over one second for SNIV and the AR confidence set/interval. In the design with conditional heteroskedasticity the SNIV confidence set is marginally wider than under homoskedasticity and attains the nominal coverage.\\
{Many instruments.} We take the classical design and add 1989 redundant instruments, all drawn from the standard normal distribution. This yields $d_Z=1999$ instruments and $n=2000$ observations. The SNIV confidence sets have coverage almost identical to the nominal level but are wider than the classical design. However, they remain sufficiently narrow as to be informative on the sign of the nonzero entries of $\beta^*$. In contrast, the AR confidence sets and intervals do not have the correct coverage and are too wide so as to be informative. We also consider an identical design but with $d_Z=2100$. The SNIV confidence sets are similar to the case of $d_Z=1999$, whereas AR confidence sets and intervals are not defined. Almost all optimization problems solved yielded a global optimum. The SNIV optimization problems are solved more slowly than the classical design, taking around 4 seconds on average. In the designs with conditional heteroskedasticity the SNIV confidence set has nominal coverage but is marginally narrower than under homoskedasticity, likely because a greater fraction of the optimization problems yielded an exact global optimum. \\
{Weak instruments.} We take the classical design and set $\pi^*=0.03$. The AR confidence sets are narrower than SNIV but have coverage further from the nominal level. The AR confidence intervals have coverage slightly larger than the nominal level. Almost all optimization problems solved yielded a global optimum. In the design with conditional heteroskedasticity the SNIV confidence set is marginally wider than under homoskedasticity with coverage close to the nominal level. Computation timings are similar to the classical design.\\
{Invalid instruments.} We take the classical design with one endogenous regressor and $d_Z=9$ instruments, all of which are included as regressors (i.e., $X_k$=$Z_{k-1}$ for $k=2,3,...,d_X$). Hence there are $d_X=10$ regressors, but only the first is endogenous. We set $\Pi^*=(0,0,...,0,\sqrt{\pi^*}/2,-\sqrt{\pi^*}/2)^\top$ so that only the final two instruments are correlated with the endogenous regressor. We suppose that $S_Q=[d_X]$ (i.e., the relevance of all regressors is questionable) and consider the sparsity certificates $s\in\{2,3\}$ (recalling that $\beta^*$ has two nonzero entries). This design is such that $\mathcal{I}_2$ is a singleton, $\mathcal{I}_3$ is not a singleton but is bounded, and $\mathcal{I}_s$ is unbounded for $s>3$. When $s=3$, though $\beta^*_1=1$ we have $\min_{\beta\in\mathcal{I}_3} \beta_1=0$ and $\max_{\beta\in\mathcal{I}_3} \beta_1=1$. The AR confidence sets and intervals cannot be computed.
The SNIV confidence set has coverage slightly larger than the nominal level. For $s=2$, SNIV is sufficiently narrow so as to be informative on the sign of the nonzero entries of $\beta^*$. For $s=3$, the width of SNIV for $\beta_1$ is around 1.25 on average, which is not sufficiently narrow so as to be informative on the sign of $\beta^*_1$. This is expected because $\beta^*_1$ is not point identified (the identified set has width 1), as explained in the previous paragraph. In the design with conditional heteroskedasticity the SNIV confidence set is marginally wider than under homoskedasticity. Each optimization problem is solved in around 17-27 seconds depending on the sparsity certificate used. This is likely due to the non-convex nature of the sparsity constraint. Nevertheless, the problems are sufficiently tractable so as to allow informative inference on a standard machine.\\
{Endogenous instruments.} We take the classical design and but add the instruments $Z_{11}=X_1$ and $Z_{12}=X_2$ (i.e., include the first two regressors as additional instruments). This results in $d_Z=12$ instruments, two of which are endogenous, with $\theta_{11}=\Omega_{1,2}$ and $\theta_{12}=\Omega_{1,3}$. We suppose that the researcher questions the exogeneity of the two endogenous instruments and the final three exogenous instruments, hence $S_\perp=[7]$. This implies that are seven instruments known to be exogenous, whereas $d_X=10$, hence the model using only the instruments known to be exogenous is underidentified. Thus, classical tests of overidentifying restrictions are infeasible.
Since the identified set is otherwise unbounded, we restrict the number of endogenous instruments using the sparsity certificate $\widetilde{s}$ on $\theta^*$. We do not make any sparsity restriction on $\beta^*$. Using $\widetilde{s}=2$ corresponds to the case where we assume that there are ten exogenous instruments, but we do not know all of their identities. This design is such that $\mathcal{I}_{d_X,2}$ is a singleton. In contrast, $\mathcal{I}_{d_X,\widetilde{s}}$ is unbounded for $\widetilde{s}>2$. We compute the SNIV confidence set for $\widetilde{s}=2$.
The SNIV confidence set has coverage almost exactly equal to the nominal level and is sufficiently narrow so as to be informative on the sign of the nonzero entries of $\beta^*$. Moreover, the SNIV confidence set allows the null hypotheses of $\theta_{S_\perp^c}=0$ to be (correctly) rejected with probability 0.84. In the design with conditional heteroskedasticity the SNIV confidence set is marginally wider than under homoskedasticity. Each optimization problem is solved more slowly than under the classical design, taking around 36 seconds on average. Nevertheless, the problems are sufficiently tractable so as to allow informative inference on a standard machine.
We use self-normalization of sample moments to conduct robust, computationally tractable inference in linear instrumental variables models. We also show that our computational approach is not unique to self-normalzation, and can be applied to perform fast inversion of other tests. In our view there are two avenues for future work.
First, though SNIV requires minimal assumptions and has desirable statistical and computational properties, when $d_X$ is large it can be conservative when the object of interest is low dimensional (e.g., a single treatment effect). Though this is not an issue in the leading case in which $d_X$ is small (e.g., one endogenous regressor, possibly after partialing our a small number of exogenous regressors) and $d_Z$ may be large (with possibly weak instruments), future work may seek to adapt our approach to perform robust inference directly on a sub-vector of parameters of interest.
Second, we believe that our computational approach is applicable beyond the instrumental variables context. An obvious setting to which our results may be applied is that of inference in partially identified models, which is often based on solving programming problems such as (ref). A simple example in which the optimization problem is semi-algebraic (hence our approach is applicable) is the $2\times2$ entry game considered by kaido19.