An Identification and Dimensionality Robust Test for Instrumental Variables Models

An Identification-and Dimensionality-Robust Test for Instrumental Variables Models

\thispagestyle{empty}

Introduction

When instruments are suspected to be weak, researchers may want to test hypotheses about structural parameters using testing procedures that are robust to identification strength. These procedures all rely on some conditions on the rate of growth of the number of instruments, \(d_z\), in relation to the sample size, \(n\). The initial identification robust tests developed in StaigerStock-1997, Moreira_2003, and Kleibergen-2005 are shown by Andrews-Stock-2007 to control size under heteroskedasticity when the number of instruments cubes is small relative to the sample size, \(d_z^3/n\to 0\). Meanwhile, recent and interesting “many-instrument” tests (crudu_mellace_sandor_2021, Mikusheva_Sun_2022, Matsushita_Otsu_2022, lim2022conditional) allow the number of instruments to be proportional to the sample size, \(d_z/n\to \varrho \in [0,1)\), but require that the number of instruments itself be large, \(d_z \to\infty\).

In practice, these conditions can be difficult to interpret and the variety of tests available under alternate regimes may make it difficult for the researcher to know which test, if any, should be applied in her exact setting. As examples, consider settings such as that of derenoncourt-2022, where \(d_z = 9\) and \(n = 130\), Paravisini-2014, where \(d_z = 10\) and \(n = 5{,}995\), and gilchrist-glassberg-2016 where \(d_z = 52\) and \(n = 1{,}671\). In all three cases, the number of instruments cubed, \(d_z^3 = 729\), \(1{,}000\) and \(140{,}608\), respectively, cannot be treated as negligible relative to the sample size. Indeed, all of these papers use post-LASSO estimates of the first-stage, suggesting concern about the large number of instruments by the authors. However as the number of instruments is only moderately large in all of these settings, asymptotic approximations that rely on \(d_z \to \infty\) may not accurately resemble the finite sample distribution of the test statistic. This can make size control of many-instrument tests questionable.\footnote{Mikusheva_Sun_2022 point out that, when errors are homoskedastic, under fixed \(d_z\) their proposed test statistic weakly converges to a distribution whose 95th quantile is close to that of the standard normal. However, it is unclear what the behavior of the test would be when errors are heteroskedastic.}

By comparison, the test proposed in this paper can be applied in any of the settings listed above as well as when the instruments are high-dimensional, \(d_z \gg n\), a setting for which there has been little progress on identification robust testing to this point. The main problem in these settings has been that the exact limiting behaviors of the regularized first-stage estimators used when \(d_z \gg n\) are difficult to analyze and typically unknown. Existing analyses sidestep this issue by assuming strong identification and exploiting a certain orthogonality, termed “Neyman Orthogonality” by CCDDHNR-2018, of the structural parameter estimate to the first-stage estimation error BCCH-2012. This approach is explicitly not applicable under weak identification where the first-stage estimation error is on the same or lesser order as the signal from the instruments and thus relevant to asymptotic analysis. As such, the few existing proposals for identification robust testing that allow \(d_z \gg n\) (BCCH-2012, gautier2021highdimensional, Mikusheva-2023-ManyTimeSeries) either fail to incorporate first-stage information or rely on sample splitting, both of which may reduce power.

To construct the test statistic I borrow an idea from Kleibergen_2002,Kleibergen-2005 and leverage a conditional slope parameter, which can be simply estimated using regularized methods even when \(d_z \gg n\), to partial out the structural error from the endogenous variable. I then use this partialled out version of the endogenous variable to construct first-stage estimates. The key idea is that, if variables in the model followed a jointly Gaussian distribution, we could exploit the fact that uncorrelated jointly Gaussian random variables are independent to show that these first-stage estimates are independent of the structural error under the null. Thus, even under weak identification where their behavior becomes relevant to the distribution of the proposed test statistic, one could easily derive the limiting \(\chi^2\) distribution of the proposed test statistic by conditioning on these estimates. Deriving the limiting distribution of the test statistic in the general model then reduces to showing that one can treat all observations as if they followed a jointly Gaussian distribution.

When the number of instruments is small, this can be justified in large samples through central limit and continuous mapping theorems. However, when the number of instruments is large these standard tools cannot be applied. Instead, my asymptotic analysis proceeds using modifications of Lindeberg's interpolation argument Lindeberg1922, which roughly proceeds by showing to be negligible the total change in distribution from a one-by-one replacement of each observation in the expression of the test statistic with a Gaussian version. The modifications of Lindeberg's original interpolation method are non-trivial and required to deal with the fact that derivatives of the test statistic with respect to individual observations may be unbounded generally, an issue arising from the “divide-by-zero” problem of weak identification ASS-2019.\footnote{Existing interpolation results require these derivatives to be bounded while in this setting they may not even have finite moments. Given this, these modified approaches may be of independent interest to a growing literature on direct Gaussian approximation techniques (Chatterjee-2006, cck2013, Pouzo-Quadratic, CMY-2020).} This interpolation argument requires some conditions on the first-stage estimates, in particular I require that the first-stage estimates take on a “jackknife-linear” form and suggest using a jackknife ridge regression in practice to allow \(d_z \gg n\). Interestingly, though, these first-stage estimates are not required to converge to limiting values which allows the researcher some flexibility if she wanted to use a different approach.

Through the Gaussian approximation result, I examine the power properties of my proposed testing procedure in local neighborhoods of the null. These local neighborhoods are characterized by a bounded local power index. In the case of a single endogenous variable I show that, under an additional regularity condition, the local power index diverging implies that the test is consistent. Unfortunately, the process of partialling out the structural error can introduce bias into the first-stage estimate under the alternative hypothesis. Against certain alternatives, this bias can be particularly pronounced and erase the first-stage signal from the instruments, a problem pointed out by Andrews-2016 in the context of the Kleibergen-2005 K-statistic. To address this, I propose a simple combination with the sup-score test of BCCH-2012, which can also allow \(d_z \gg n\). As with the Anderson-Rubin test, while the sup-score statistic does not incorporate first-stage information it does not face a power decline against particular alternatives.

Identification-robust testing procedures may be of particular interest in high-dimensional settings due to a lack of clarity on how to pretest for weak identification. Current empirical practice when using post-LASSO estimates of the first-stage appears to be to use standard t-test inference if the first stage F-statistic on the LASSO selected variables is larger than 10 (Paravisini-2014, gilchrist-glassberg-2016, derenoncourt-2022). Using a simple numerical demonstration, I argue that first stage F-statistics on LASSO selected variables may not be reliable indicators of identification strength. Given uncertainty about the strength of identification I apply the newly proposed testing procedures to the data of gilchrist-glassberg-2016 and generate weak instrument-robust confidence intervals for the effect of social spillovers on movie consumption. The newly proposed confidence intervals are smaller than those obtained by inverting the sup-score and many instrument tests across all specifications. To verify these results, I also generate confidence intervals in the data of AK-1991, again finding that the newly proposed methods generate tighter confidence bands than the many instruments methods. In (ref), I argue these improvements in power may be explained by my proposed test statistic's use of higher quality first-stage estimates and individual scores that are uncorrelated across observations.

Finally, I examine the applicability of the theoretical results in this paper through a simulation study. While existing tests seem to face size distortions in alternate regimes, the test based on my proposed test statistic has nearly exact size in a variety of settings. While this new test may have diminished power against certain alternatives, this deficiency is ameliorated through combination with the sup-score test. After combining the new test with the sup-score statistic, I find the newly proposed testing is also have favorable power properties compared to the many instruments and sup-score test in my simulation setup.

The outline of this paper is as follows. (ref) formally defines the model considered and introduces the jackknife K-statistic and it's combination with the sup-score test. (ref) demonstrates the usefuleness of the proposed testing procedures in two empirical applications while (ref) provides evidence from simulation study. (ref) provides an overview of the Gaussian approximation approach and characterizes the limiting behavior of the test statistic. (ref) uses this characterization to examine the power properties of the test and establishes the validity of the combination test. Proofs of the main results as well as the presentation of some auxiliary results are deferred to the Appendix.

\paragraph{Notation.} For any \(n \in \SN\) let \([n]\) denote the set \(\{1,\dots,n\}\). I work with a sequence of probability measures \(P_n\) on the data \(\{(y_i,x_i,z_i): i \in [n]\}\). To accommodate independent but not identically distributed observations, let \(\E_n[f_i] = n^{-1}\sum_{i=1}^n f_i\) denote the empirical expectation and \(\bar\E[f] = \E_n[\E[f_i]]\) denote the average expectation operator.

Model and Setup

Though the analysis below allows for exogenous regressors, to simplify the exposition I follow Mikusheva_Sun_2022 and assume that they have already been partialed out of both the outcome, \(y_i\), and the endogenous regressors, \(x_i\). As the controls are assumed to be of fixed dimension, this is without loss of generality.\footnote{For discussion refer to (ref).} Along with the first stage, the IV model can then be written as a system of simultaneous equations:

The researcher observes the outcome \(y_i \in \SR\), the endogenous variable \(x_i \in \SR^{d_x}\), and the instruments \(z_i \in \SR^{d_z}\) but neither the structural error \(\varepsilon_i \in \SR\) nor the first-stage errors \(v_i \in \SR^{d_x}\). The structural error is assumed to be conditional-mean independent of the instruments, \(\E[\varepsilon_i|z_i] = 0\). I denote \(\E[x_i|z_i]\) as \(\Pi_i \coloneqq \E[x_i|z_i]\) and make no assumptions about the functional form of the conditional expectation so the instruments are allowed to affect the endogenous variable in a nonlinear fashion.

The random variables \(\{(z_i,\varepsilon_i, v_i)\}_{i=1}^n\) are assumed to be independent across observations. Observations need not be identically distributed but the errors are assumed to have a common covariance structure conditional on the instruments \(z_i\): \[ \Var((\varepsilon_i, v_i)'|z_i) \coloneqq \Omega(z_i) =

\in \SR^{(1 + d_x) \times (1 + d_x)} \] As \(\Omega(z_i)\) is otherwise left unrestricted, the errors are allowed to be heteroskedastic. All results in this paper hold conditionally on a realization of the instruments \(\vz := (z_1',\dots, z_n') \in \SR^{n\,\times\,d_z}\) so from this point forth they are treated as fixed and all expectations can be understood as conditional on the instruments.

Under this setup, the researcher wishes to test a two-sided restriction on the structural parameter: \[ H_0: \beta = \beta_0\vsbox H_1: \beta \neq \beta_0 \] I am interested in constructing powerful tests for this null-alternate pair that are asymptotically valid under arbitrarily weak identification and with minimal restrictions on the number of instruments \(d_z\). To this end, define the null errors \(\varepsilon_i(\beta_0) \coloneqq y_i - x_i'\beta_0\). Using these, I construct a variable, \(r_i\), that is a “partialed-out” version of the endogenous variable satisfying \(\Cov(r_i,\eps_i(\beta_0))= 0\):

Each element of the nuisance parameter \(\rho(z_i)\), \(\rho_\ell(z_i)\) for \(\ell = 1,\dots,d_x\), can be interpreted as the (conditional) slope coefficient from a simple linear regression of \(x_{\ell i}\) on \(\eps_i(\beta_0)\). Thus, if \(\rho_\ell(\cdot)\) falls in some function class \(\Phi\) it can be estimated directly under \(H_0\) by solving empirical analogs of:\footnote{Under \(H_1\), \(\rho_\ell(z_i)\) can be estimated directly by solving empirical analogs of \(\rho_\ell(z_i) = \arg\min_{\phi \in \Phi} \E[(x_{\ell i} - \eta_i(\beta_0)\phi(z_i))^2]\) where \(\eta_i(\beta_0) = \eps_i(\beta_0) - \E[\eps_i(\beta_0)|z_i]\). This requires an initial estimate of \(\E[\eps_i(\beta_0)|z_i]\), however.} \[ \rho_\ell(z_i) = \arg\min_{\varphi\in\Phi} \bar\E[(x_{\ell i} - \eps_i(\beta_0)\varphi(z_i))^2]. \] While other estimators of \(\rho(\cdot)\) are, in principle, possible (see (ref)), I will focus on \(\ell_1\)-penalized/LASSO estimators. These estimators are consistent under the assumption that \(\rho(z_i)\) has an approximately sparse representation in some basis \(b(z_i) \coloneqq (b_1(z_i),\dots,b_{d_b}(z_i))' \in \SR^{d_b}\), that is \(\rho_\ell(z_i) = b(z_i)'\phi_\ell + \xi_{\ell i}\) where \(\xi_{\ell i}\) represents an approximation error that tends to zero with the sample size and \(\phi_\ell\) is sparse in the sense that many of its coefficients are zero. This allows for nesting of the low-dimensional case, where the number of instruments is fixed, and the high dimensional case, where the number of instruments is potentially much larger than the sample size, under a unified estimation procedure. Under homoskedasticity, \(\rho_\ell(z_i)\) is a constant function and thus has a spare representation in any basis that contains a constant term. In general, the aproximate sparsity assumption can either be interpreted as an assumption that there are only a few instruments that are important for explaining variation in the covariance matrix \(\Omega(z_i)\) or as an assumption that the function \(\rho(z_i)\) can be accurately approximated using only a smaller set of basis terms in \(b(z_i)\).

The parameter \(\phi_\ell\) can be estimated via LASSO:

or via post-LASSO, refitting an unpenalized version of (ref) using only the basis terms associated with nonzero coeffecients in the inital LASSO regression. The estimating procedure in (ref) is a simple \(\ell_1\)-penalized regression of \(x_{\ell i}\) against \(\eps_i(\beta_0)b(z_i)\). It can be easily implemented using out-of-the-box software available on most platforms. Under standard conditions, this leads to a consistent estimate of \(\rho_\ell(z_i)\) as long as the sparsity condition \(s^2\log^M(d_bn)/n \to 0\) where \(s\) is the number of nonzero elements of \(\phi_\ell\) and \(M\) is a positive constant that depends on the moment bounds imposed. The estimation procedure is discussed in more detail in (ref). With \(\hat\rho(z_i) \coloneqq b(z_i)'\hat\phi_\ell\), I construct the estimated version of \(r_{\ell i}\), \(\hat r_{\ell i} \coloneqq x_i - \hat\rho(z_i)\eps_i(\beta_0)\) for each \(\ell \in [d_x]\).

Test Statistic

The test statistic is based on an arbitrary jackknife-linear estimate of the first stage, \[ \widehat\Pi_{\ell i} = \sum_{j\neq i} h_{ij}\hat r_{\ell j},\;\; \ell \in [d_x] \] for some “hat” matrix \(H = [h_{ij}] \in \SR^{n\times n}\). The phrase “hat matrix” is borrowed from ordinary least squares (OLS) where the projection matrix, \(\vz(\vz'\vz)^{-1}\vz'\), is sometimes referred to as the hat matrix in the sense that \(\hat x = \vz(\vz'\vz)^{-1}\vz' x\). In practice, the hat matrix, \(H\), can be any matrix that depends only on \(\vz\). It is important to note that while \(\widehat\Pi_{\ell i}\) does not depend on \(\hat r_{\ell i}\), it may depend on \(z_i\) through the hat matrix \(H\). This gives the test power against alternatives where \(\E[\eps_i(\beta_0)z_i] \neq 0\). For technical reasons, I will assume that \(h_{ii} = 0\) for each \(i \in [n]\) so that \(\widehat\Pi_{\ell i}\) can be written as \(\widehat\Pi_{\ell i} = \sum_{j=1}^n h_{ij} r_{\ell j}\).

Formally, the only structure I require on the hat matrix \(H\) is a balanced-design condition described in (ref). However, for reasons explained in (ref) it may be optimal to introduce some regularization in estimating the first-stage models \(\widehat\Pi_{\ell i}\) so I suggest using a jackknife ridge regression procedure setting:

where \(\hat\pi_{(-i)}(\lambda)\) is the coefficient from a ridge regression of \(\hat r\) on \(\vz\), leaving out observation \(i\) and with penalty parameter set equal to \(\lambda\). \footnote{A ridge regression coefficient estimate from a regression of \(\tilde\vy \in \SR^n\) on \(\tilde\vx\in \SR^{k\times n}\) with penalty parameter \(\lambda \in \SR_+\) solves \(\hat\pi \in \arg\min_{\pi\in\SR^K} \|\tilde\vy - \pi'\tilde\vx\|_2^2 + \lambda\|\pi\|_2^2\). nyquist1988applications shows how the fitted values from a jackknife ridge regression can be calculated without having to recompute \(\hat\pi_{(-i)}\) for each observation. Angrist_Krueger_JIVE provides a similar analysis for jackknife OLS.} Following recommendations in harrell_2015 and lecture_notes_ridge, the penalty parameter \(\lambda^\star\) is set so that the effective degrees of freedom is no more than a fraction of the sample size: \[ \lambda^\star = \inf \{\lambda \geq 0: \trace(\vz(\vz'\vz + \lambda I_{d_z})^{-1}\vz') \leq n/5\} \] The jackknife ridge first-stage estimate has the benefit of being well defined even when the number of instruments is larger than the sample size. I stress, though, that the \(\widehat\Pi_{\ell i}\) estimators are not required to be consistent and the researcher may use any other hat matrices that she believes will lead to plausible first-stage estimates. Other possible choices of hat matrix include the jackknife OLS hat matrix of Angrist_Krueger_JIVE, the deleted diagonal projection matrix introduced in Chao-2012-JIVE and successfully used in KSS-2018,crudu_mellace_sandor_2021,Mikusheva_Sun_2022, and Matsushita_Otsu_2022, or hat matrices based on selecting instruments via some preliminary unsupervised technique such as principal component analysis (PCA). (ref) below discusses how the balanced-design condition may be verified for arbitrary choices of hat matrices.

For each \(i = 1,\dots,n\), define \(\widehat\Pi_i = (\widehat\Pi_{1i},\dots,\widehat\Pi_{d_xi})\in \SR^{d_x}\) and \(\widehat\Pi_{\eps i} = \eps_i(\beta_0)\widehat\Pi_i\). Collect these in the matrices

The jackknife K-statistic can then be defined

I will show that, under appropriate moment bounds and conditions on the hat matrix, \(H\), the limiting distribution of \(\JK(\beta_0)\) under \(H_0\) is \(\chi^2_{d_x}\). For exposition, I will largely focus on the case where \(d_x = 1\), in which case the form of the test statistic simplifies to \(\JK(\beta_0) = \big(\sum_{i=1}^n \eps_i(\beta_0)\widehat\Pi_i\big)^2/\sum_{i=1}^n \eps_i^2(\beta_0)\widehat\Pi_i^2\). The extension to \(d_x > 1\) is not immediate but is possible under strengthened moment conditions.

Combination with Sup-Score Test

As will be discussed further in (ref), the test based on the \(\JK(\beta_0)\) statistic can have deficient power against certain alternatives. This loss of power is similar to that faced by tests based on the K-statistics of Kleibergen_2002,Kleibergen-2005 and derives from the fact that the process of partialling out the null errors, \(\eps(\beta_0)\), from the endgenous variables introduces bias into the first stage estimates, \(\widehat\Pi\), under the alternative hypothesis. Against certain alternatives, this bias can “erase” the first stage signal from the instruments.

The power deficiency in tests based on the K-statistic is typically adressed by combining the K-statistic with the Anderson-Rubin statistic based on a constructed conditioning variable. Prominent examples of such combinations include the celebrated conditional likelihood test of Moreira_2003, GMM-M test of Kleibergen-2005, and minimax regret tests of Andrews-2016. I take a similar approach here in combining the newly proposed tests with tests based on the sup-score statistic of BCCH-2012,

which have correct asymptotic size even when the instruments is much larger than the sample size, \(d_z \gg n\). A level \(\alpha \in (0,1)\) test based on the sup-score statistic rejects whenever \(S(\beta_0) > c^S_{1-\alpha}\) where, for \(e_1,\dots,e_n\) iid standard normal and generated independently of the data, \(c^S_{1-\alpha}\) is the simulated multiplier bootstrap critical value:\footnote{This conditional quantile can also be approximated using an empirical bootstrap procedure as demonstrated by deng2020beyond.} \[ c^S_{1-\alpha} \coloneqq (1-\alpha)\text{ quantile of } \sup_{1 \leq \ell \leq d_z} \bigg|\frac{\sum_{i=1}^n e_i \eps_i(\beta_0)z_{\ell i}}{(\sum_{i=1}^n z_{\ell i}^2)^{1/2}}\bigg|\text{ conditional on }\{(y_i,x_i,z_i)\}_{i=1}^n. \] As with the Anderson-Rubin test, tests based on the sup-score statistic may have suboptimal power properties in overidentified models as it does not incorporate first-stage information. However, the sup-score statistic does retain the benefit of directing power evenly in all directions, avoiding pitfalls of tests based on \(\JK(\beta_0)\) against certain alternatives.

The combination test decides which test to run based on an attempt to detect whether the alternative \(\beta\) is such that \(\E[\widehat\Pi_{\ell,i}^I] = 0\) for all \(i \in [n]\) and some \(\ell \in [d_x]\). When this happens, tests based on the \(\JK(\beta_0)\) statistic have trivial power against deviations in the \(\ell\textsuperscript{th}\) coordinate of \(\beta\), so in local neighborhoods of these values of \(\beta\), tests based on the sup-score statistic may be preferable. Detection of whether \(\E[\widehat\Pi_{\ell,i}^I] = 0\) is based on the conditioning statistic:

Under the assumption that \(\E[\widehat\Pi_i^I] = 0\) for all \(i \in [n]\), quantiles of the conditioning statistic can be simulated analogously to the sup-score critical value. For a new set of \(\{(e_{\ell1},\dots,e_{\ell n}): \ell \in [d_x]\}\) iid standard normal and generated independently of the data, and for any \(\theta\in (0,1)\), define the conditional quantile

The thresholding test decides which test to run by comparing the conditioning statistic \(C\) to a treshhold value \(\tau\),

In principle, the thresholding statistic has correct size for any (preset) choice of parameter \(\tau\). In practice, however, I find that setting \(\tau = c^C_{0.75}\) leads to a reasonable balance of power between local and distant alternatives.

Empirical Application

I apply the testing procedures proposed in this paper to the data of gilchrist-glassberg-2016, who examine the effect of social spillovers in movie ticket sales, and to the data of AK-1991, who examine the returns to education. In both studies the number of instruments, \(d_z = 52\) and \(d_z = 180\), respectively, cannot be treated as negligible relative to the sample size (\(n = 1{,}671\) and \(n = 329{,}509\), respectively). To deal with the large number of instruments, gilchrist-glassberg-2016 employ a post-LASSO estimate of the first stage. This strategy is also shown to work well in the data of AK-1991 by angrist2022machine. Using a simple simulation study, I demonstrate that the first-stage F-statistics on LASSO selected variables typically reported by authors can be misleading indicators of identification strength. When revisting the initial analyses using identification robust testing procedures, the confidence intervals constructed by inverting the tests proposed in (ref) are consistently narrower than those constructed from inverting the sup-score and many-instrument testing procedures.

Application to Social Spillovers in Movie Consumption

The gilchrist-glassberg-2016 sample consists of 1,671 opening weekend days between January 1, 2002 and January 1, 2012. For each opening weekend, the authors observe gross ticket sales for movies wide released in theaters in the United States with a run in theaters of at least six weeks.\footnote{An opening weekend day is a Friday, Saturday, or Sunday of opening weekend and a wide released movie is any movie that ever shows on 600 or more screens.} The data are obtained through Box Office Mojo, a subsidiary of the Internet Movie Database (IMDb).

The outcome variables of interest are gross ticket sales of movies that opened in a given weekend in the second through sixth weeks of their run, while the endogenous variable is the gross ticket sales of a movie in its opening weekend. To control for seasonal periodicity in both the supply of and demand for movies, a vector of date controls are included. Formally, the authors are interested in the parameters \(\beta_w\), \(w = 2,\dots,7\) from the linear IV model(s):

where, for \(w = 1,\dots,6\), \(\text{Sales}_{wi}^\perp\) represents gross national ticket sales, after the partialing out of date controls and a constant, \(7(w-1)\) days after day \(i\), of movies that opened on the opening weekend of \(i\). The variable \(\text{Sales}_{7i}^\perp = \sum_{w = 1}^6 \text{Sales}_{wi}^\perp\) denotes the cumulative national ticket sales in the second through sixth running weekends of movies who opened in weekend \(i\), after the partialing out of date controls and a constant. The parameter \(\beta_w\) represents the social spillover effect of strong opening weekend sales on sales in later weeks.

To instrument for sales on opening weekend the authors employ a vector of nationally aggregated weather measures. These weather measures include the proportion of movie theaters experiencing maximum temperatures in \(5^\circ\) Fahrenheit bins on the interval \([10^\circ, 100^\circ]\), the proportion of movie theaters experiencing precipitation levels in 0.25 inch per hour increments on the interval \([0,1.5]\), and the proportions of theaters experiencing any type of snow and of theaters experiencing any type of rain. Since unusually poor weather may cause people to substitute away from outdoor activities and into watching a movie, these measures provide a source of exogenous variation in opening weekend sales that can be used to identify the effect of social spillovers.

Putting together the nationally aggregated weather measures leaves gilchrist-glassberg-2016 with 48 linearly independent instrumental variables. \footnote{There are 52 instruments in total, but four linearly dependent ones are ignored in the following.} To handle the large number of instruments, the authors employ a post-LASSO estimate of the first stage BCCH-2012; they set the first-stage penalty parameter so that the number of instrument selected is one, two, or three. The resulting first-stage F-statistics using the selected instrument(s), 38.80, 25.86, and 20.95, respectively, seem to indicate strong identification. However, the first-stage F-statistic on the full set of instrumental variables is only 3.80. Moreover, since the LASSO objective is an \(\ell_1\) penalized version of the OLS loss, using the variables selected by LASSO may mechanically lead to higher F-statistics even if the underlying relationship between the instruments and the endogenous variables is weak.

(ref) provides evidence from a simple simulation experiment to demonstrate this. For the simulation experiment I generate an iid sample of size \(n = 1000\). For each \(i \in [n]\), I generate 10 mutually independent instruments \(Z_{ki} \sim N(0,1)\) for \(1 \leq k \leq 10\). The endogenous variable is generated to only have a weak relationship with the instruments, \(X_i = \frac{2}{\sqrt{n}}\sum_{\ell = 1}^{10} Z_{\ell i}+ v_i\), and the outcome is generated \(Y_i = X_i + \eps_i\) where \((\eps_i,v_i)\) are independent standard normals. From the initial set of 10 instrumental variables I generate an additional 55 technical instruments by squaring and taking all interactions between variables in the initial set. These generated instruments are correlated with the initial instruments but do not directly enter the first stage.

I then set the LASSO penalty so that only a certain number of instruments are chosen and report the resulting average first stage F-statistics and 95% confidence interval coverage over one thousand simulations. As a comparasion I also report the average first-stage F-statistics and 95% confidence interval coverage from the oracle estimator, which only uses the relevant 10 initial instruments. Despite the fact that the first-stage F-statstic on selected instruments is more than double the first-stage F-statistic using the oracle first stage estimator, the coverage rate of 95% confidence intervals based on LASSO selected instruments is significantly degraded compared to both the nominal coverage probability and the coverage probability using the oracle first-stage estimator.

Given a lack of clarity on the strength of identification, I seek to validate the results of gilchrist-glassberg-2016 using the weak identification testing procedures proposed in this paper. The setting is particularly suitable for weak IV testing using the jackknife K-statistic. With 48 instruments and a sample size of 1671, \(d_z^3 = 110{,}592 \gg n \), making the tests of Moreira_2003,Moreira_2009,Kleibergen-2005, and Andrews-2016 inapplicable. On the other hand, it is unclear whether asymptotic approximations based on \(d_z \to \infty\) will accurately describe the finite-sample distribution of test statistics with 48 instruments. Moreover, since fluctuations in movie theater attendance seem to be largely driven by either particularly cold or particularly hot weather (see Figure 4 in gilchrist-glassberg-2016), the nuisance parameter \(\rho(z_i)\) is plausibly approximately sparse.

I compare the 95% confidence intervals based on the jackknife K-statistic to those based on the jackknife Lagrange-Multiplier (JLM) statistic of Matsushita_Otsu_2022 and the sup-score statistic of BCCH-2012. Confidence intervals based on the jackknife AR statistic of Mikusheva_Sun_2022 are not reported as they were empty for all specifications, a result that could indicate misspecification of the linear model. Similarly, confidence intervals based on the thresholding statistic, implemented as recommended in (ref), are also not reported as they always align with the those of jackknife K-statistic. The narrower confidence bands of the jackknife K-statistic in specifications where the sup-score confidence interval is non-empty indicates higher power from the jackknife K-test in this setting, so the combination test suggesting its use may be expected. The jackknife K-statistic is implemented using the jackknife OLS hat matrix of Angrist_Krueger_JIVE, that is by setting \(\lambda = 0\) in (ref).

Tables (ref) reports the \(95\%\) confidence intervals for \(\beta_1,\dots,\beta_7\) generated by weak-instrument robust confidence intervals for three sets of instruments: the first is the initial set of 48 instruments in gilchrist-glassberg-2016, the second set includes only the temperature instruments for \(d_z = 36\), and the final includes the initial instruments as well as all interactions between the temperature instruments and the remaining instruments for \(d_z = 524\). For reference, I also provide point estimates and standard errors for \(\beta_1,\dots,\beta_7\) from gilchrist-glassberg-2016, Table 2. To facilitate comparison, these point estimates and standard errors come from a specification that uses all the instruments in the first stage of a 2SLS procedure.

Qualitatively, the results from the weak-instrument robust confidence intervals are similar to that of the author's original analysis; indeed the gilchrist-glassberg-2016 point estimates are always in the identification robust confidence intervals when using either the initial instrument set (\(d_z = 48\)) or the reduced instrument set (\(d_z = 36\)). However, when using the larger instrument set of \(d_z = 524\) we obtain confidence bands using the jackknife K-test that rule out the author's initial point estimates for the parameters \(\beta_5, \beta_6\) and \(\beta_7\) and suggest somewhat smaller social spillover effects in movie consumption. Across all specifications and instrument sets, confidence intervals obtained from inverting the jackknife K-test are consistently smaller than those obtained from inverting both the sup-score and jackknife Lagrange Multiplier test. These reductions in confidence interval length are most noticeable when using the complete set of interactions, across all parameters the jackknife K-confidence intervals are nearly half the length of their jackknife Lagrange Multiplier counterparts. As with the jackknife Anderson-Rubin test, the sup-score confidence intervals are also often empty which again could suggest misspecification of the linear IV model.

Application to Returns to Education

For additional comparison, I revisit the setting of AK-1991, who study the effect of education on log-wages, instrumenting for education using various combinations of quarter-of-birth (QOB), year-of-birth (YOB), and place-of-birth indicators (POB). This dataset has the benefit of being used in the empirical studies of both Mikusheva_Sun_2022 and Matsushita_Otsu_2022, facilitating an easy comparison of the performance of the newly proposed testing procedure to these existing “many-instrument” methods. Moreover, angrist2022machine report improved power in this dataset when using post-LASSO estimates in the first stage, suggesting that high-dimensional or machine-learning techniques can be useful in this setting.

(ref) displays the resulting identification robust confidence intervals using two set of instruments. The first set, which contains all QOB \(\times\) YOB and all QOB \(\times\) POB interactions for \(d_z = 180\), corresponds to the specification in Table VII of AK-1991. The second set of instruments, \(d_z = 1{,}530\), contains all QOB \(\times\) YOB \(\times\) POB interactions. The confidence intervals for the jackknife Anderson-Rubin (JAR) and jackknife Lagrange-Multiplier (JLM) tests are taken directly from Mikusheva_Sun_2022 and Matsushita_Otsu_2022, respectively. To implement the jackknife K-test I follow a similar procedure to that of the prior empirical exercise, setting \(\lambda = 0\) when constructing the \(\widehat\Pi_i\) values. However, for computational reasons I opt to split the data into 11 pieces when estimating \(\widehat\Pi_i\) rather than using a true “leave-one-out” jackknife approach. Confidence intervals from the combination test are not reported as the sup-score confidence interval is empty in both specifications.

The results in (ref) are similar to those in (ref). In both specifications the confidence intervals obtained from inverting the jackknife K-test are substantially smaller than those obtained by inverting the “many-instrument” JAR and JLM tests. Indeed, as before, the confidence intervals obtained from inverting the jackknife K-statistic are nearly half the length of those obtained from inverting the JLM test and less than a quarter the length of those obtained from inverting the JAR test. It is notable that the range of plausible returns to education values implied by the jackknife K-confidence interval with \(d_z = 1{,}530\) are strictly below that implied by the jackknife-K confidence interval with \(d_z = 180\). This may be related to misspecification of the linear model as indicated by the empty sup-score confidence interval. “Downward bias” of the confidence intervals when using the larger instrument set is also seen for the JAR and JLM confidence intervals, though to a lesser extent. (ref), below, provides reasoning for why confidence intervals based on the \(\JK(\beta_0)\) statistic may be tighter than those based on the JAR and JLM statistics.

As always, these are just particular empirical examples, and my results should not be interpreted as a critique of BCCH-2012, Mikusheva_Sun_2022, or Matsushita_Otsu_2022, whose prior work I relied upon and was inspired by.

Simulation Study

In a simple simulation study, I examine the performance of tests based on the \(\JK(\beta_0)\) statistic and compare it with that of other tests that may be used in settings where the number of instruments is nonnegligible as a fraction of sample size. I consider a reduced-form data-generating process (DGP) similar to that of Matsushita_Otsu_2022. The outcome variable, \(y_i\), and endogenous variable, \(x_i\), are generated according to

where \(\Pi_i = \frac{1}{r_n}\sum_{k=1}^5 \frac{3}{4}\bar z_{ki} + \frac{1}{4}\bar z_{ki}^2 + \frac{1}{4}\bar z_{ki}^3\) is a (dense) transformation of an initial set of instruments \(\bar z_i \in \SR^{15}\) generated as described below. The value of \(r_n\) varies depending on the strength of identification considered; under weak identification, \(r_n = n^{-1/2}\) while for power curves I consider an intermediate identification strength, \(r_n = n^{-1/3}\).\footnote{Intermediate identification strength is considered to let the power curves come up to one at the boundaries of the considered range of \(\beta\) values. Power curves with \(r_n = n^{-1/2}\) look similar, but shrunk towards zero.} To model heteroskedasticity, the errors \((e_i,v_i)\) are generated \(\eps_i = (1 + \varrho_1(\bar z_{1i}^2 + \bar z_{2i}^2 + \bar z_{2i}\bar z_{3i}))e_{1i}\), and \(v_i = \varrho_2(1 + \bar z_{1i})\eps_i + (1 - \varrho_2)^2e_{2i}\) where \(e_{1i}\) and \(e_{2i}\) are generated independently of each other and other variables in the model according to a Laplace distribution with location parameter \(\mu = 0\) and scale parameter \(b = 1\). Since jackknife K-statistic has an exact \(\chi^2\) distribution when the errors are jointly Gaussian and \(\rho(z_i)\) is known, I purposefully avoid normally distributed errors to investigate the quality of asymptotic approximations. The parameters \(\varrho_1\) and \(\varrho_2\) control the degree of heteroskedasticity and endogeneity, respectively.

I examine the size of the test under three different instrument regimes. In all three regimes, I begin with an initial set of instruments \(\bar z_i = (\bar z_{1i},\dots,\bar z_{15i})'\) generated independently across indices according to a multivariate Gaussian distribution with Toeplitz covariance structure, \(\Cov(\bar z_{\ell i},\bar z_{ki}) = 2^{-|\ell-k|}\). In the first regime, the instruments only include the initial set, \(z_i = \bar z_i\) so that \(d_z = 15\). In the second regime, the full set of instruments \(z_i\) additionally includes all quadratic and cubic terms, \((z_{\ell i}^2,z_{\ell i}^3)\), \(\ell = 1,\dots,15\) so that in total \(d_z = 45\). In the third regime, the full set of instrument includes the initial set of instruments, \(\bar z_i\), all quadratic and cubic terms (30 additional terms) and interactions of the initial set of instruments (\(\binom{15}{2} = 105\) additional terms), so that in total \(d_z = 150\). Under each regime, the full set of instruments is passed to the test statistics with no indication about which instruments correspond to the initial set, and thus no indication about which instruments are relevant to the DGP.

In constructing the jackknife K-statistic, I opt to use the deleted diagonal ridge matrix, \(H = [h_{ij}]\) where \(h_{ij} = [\vz(\vz'\vz + \lambda I_{d_z})\vz']_{ij}\bm{1}\{i \neq j\}\), instead of a true jackknife ridge procedure in order to make the simulations compuationally tractable. Following reccomendations in lecture_notes_ridge, the penalty paramter \(\lambda\) is set so that effective degrees of freedom of the resulting hat matrix is no more than \(n/5\).\footnote{Precisely, the penalty parameter is set \(\lambda = \max{0, (n/5)^{-1}d_1^2(\vz)}(d_z - n/5)\), where \(d_1^2(\vz)\) is the square of the maximum singular value of the design matrix \(\vz\).} To estimate the parameter \(\rho(z_i)\), I implement the default cross-validated \(\ell_1\)-penalized procedure of (ref) via the glmnet package in R glmnet. I use the full vector of instruments as the basis to approximate \(\rho(z_i)\).

I compare the simulated size of the jackknife K test and to the performance of the sup-score test, \(S(\beta_0)\), of BCCH-2012, the thresholding test introduced in (ref), the standard Anderson-Rubin (A.Rbn.) test of Anderson-Rubin-1949 and StaigerStock-1997, the jackknife AR test (JAR) of crudu_mellace_sandor_2021 and Mikusheva_Sun_2022, and the jackknife LM test (JLM) of Matsushita_Otsu_2022. Critical values of the sup-score and conditioning statistic are simulated with the procedures described in (ref) with 1000 bootstrap replications. For the combination test cutoff, I report results using two different quantiles of the conditioning statistic under the assumption that \(\E[\widehat\Pi_i^I] = 0\) for all \(i \in [n]\); \(\tau_{0.3}\) corresponding to the 30\textsuperscript{th} quantile and \(\tau_{0.75}\) corresponding to the 75\textsuperscript{th} quantile.

(ref) reports the simulated size for all tests under weak and strong identification, respectively. One can see that the \(\JK(\beta_0)\) statistic has nearly exact size in almost all the setups considered. In contrast, both the jackknife AR and jackknife LM test all exhibit moderate size distortions in various regimes. The jackknife AR test in particular appears to overreject in nearly all setups considered. This property also observed in the simulation study of Matsushita_Otsu_2022 and so may be driven by the similarity of this simulation setup to theirs. The jackknife LM statistic appears to have good size properties when \(d_z = 10\) and \(d_z = 150\) but is consistently (though only moderately) undersized in the intermediate regime where \(d_z = 45\). This size distortion does not improve when moving from \(n = 200\) to \(n = 500\), suggesting that the requirement that \(d_z \to \infty\) is important for the quality of finite-sample approximation by its limiting distribution. Though the good performance of the jackknife LM statistic when \(d_z = 10\) is notable, it should also be remarked that this is the setup with the least amount of correlation between the instruments.

The sup-score and Anderson Rubin test seem to be undersized in all regimes, with the Anderson-Rubin test nearly never rejecting when \(d_z = 150\). However, and in line with the theory, both of their size properties appear to improve when increasing the sample size from \(n = 200\) to \(n = 500\). It is possible that the size properties of the sup-score test could be improved by using an empirical bootstrap based approach to simulate the critical value, as proposed by deng2020beyond, however I do not consider that approach here. The threshholding test appears to inherit the conservative nature of the sup-score test, though to a lesser degree due to the combination with the \(\JK(\beta_0)\) test.

In addition to examining the size properties of the given tests I also investigate the power properties of the tests in this setting. (ref) plot calibrated local power curves under an intermediate identification strength where the first stage is in a \(n^{-1/3}\) neighborhood of zero for \(n \in \{200, 500\}\), the number of instrument is plausible large, \(d_z = 150\), \(\varrho_1 \in \{0.2,0.5\}\) and \(\varrho_2 \in \{0.3,0.6\}\). The critical value of each test is set to simulated 95\textsuperscript{th} quantile of the distribution of the corresponding test-statistic under \(H_0\). I compare the calibrated local power curves of the \(\JK(\beta_0)\) test, the combination test with cutoff \(\tau_{0.75}\), the jackknife AR test, the Jackknife LM test, and the sup-score test.\footnote{For both the jackknife AR and jackknife LM tests, I use cross-fit estimates of test statistic variances proposed and shown to improve power by Mikusheva_Sun_2022.}

In all setups considered, the jackknife K and combination tests have substantially stronger power than the jackknife AR, jackknife LM, and sup-score tests in local neighborhoods of the null as well as for negative values of \((\beta_0 - \beta)\). For values of \((\beta_0 - \beta)\) larger than 1.5, tests based on the jackknife K-statistic suffer from a loss of power as described in (ref). This power decline is largerly ameliorated by combining the jackknife K-statistic with the sup-score statistic and the thresholding test has good power properties over all alternatives considered. However, tests based on the jackknife AR or jackknife LM statistic can still provide better power than the threshholding test for very positive values of \((\beta_0 - \beta)\) in some setups.

Limiting Behavior of the Test Statistic

The limiting behavior of the test statistic is analyzed via a direct Gaussian approximation technique. In this section, I describe the approach and characterize the limiting behavior of the test statistic under local alternatives to \(H_0\). This direct approach has the advantage of not relying on any particular central limit theorem, which allows a great deal of flexibility in the choice of hat matrix \(H\) used to construct the first stage estimates. When there is only a single endogenous variable, \(d_x = 1\), the approach can be considerably simplified. The extension to \(d_x > 1\) requires a more involved argument which relies on stronger moment conditions.

Formally, I show that quantiles of the jackknife K-statistic can be approximated by analogous quantiles of the Gaussian statistic:

where, for each \(i \in [n]\), \((\tilde\eps_i(\beta_0),\tilde r_i')'\) are generated independently across indices following a Gaussian distribution with the same mean and covariance matrix as \((\eps_i(\beta_0),r_i')'\). Further, for \(\tilde\Pi_{\ell i} = \sum_{j\neq i} h_{ij} \tilde r_{\ell j}\) define \(\tilde\Pi_i \coloneqq (\tilde\Pi_{1i},\dots,\tilde\Pi_{d_xi})'\in\SR^{d_x}\), \(\tilde\Pi_{\eps i} \coloneqq (\E[\eps_i^2(\beta_0)])^{1/2}\tilde\Pi_i\),

As uncorrelated jointly Gaussian random variables are independent, under \(H_0\) the vector \(\tilde\eps(\beta_0)\) is mean zero and independent of \((\tilde\Pi, \tilde\Pi_\eps)\). Conditional on any realization of \((\tilde\Pi,\tilde\Pi_\eps)\) the \(\JK_G(\beta_0)\) statistic then follows a \(\chi^2_{d_x}\) distribution, and thus, its unconditional distribution is also \(\chi^2_{d_x}\).

To outline the approach, consider functions \(\varphi_\gamma(\cdot) \in C_b^3(\SR)\) that approximate the indicators \(\bm{1}\{\cdot \leq a\}\), where \(a \in \SR\) is arbitrary and as \(\gamma \to 0\) the quality of the approximation improves but the derivatives of \(\varphi_\gamma\) become larger in magnitutde. A primary goal is to show, for a sequence \(\gamma_n\) tending to zero, that

for a version of the test statistic, \(\JK_I(\beta_0)\), that could be constructed if \(\rho(\cdot)\) was known to the researcher. In order to establish (ref), existing interpolation methods cannot be applied as they require the derivative of the test statistic with respect to individual obserations are bounded. In this setting, the derivative of the test statistics with respect to terms in the denominator matrix, \(\widehat\Pi_\eps'\widehat\Pi_\eps\), may be as large as the inverse of the minimum eigenvalue of the denominator matrix. When identification is sufficiently weak, the eigenvalues of the denominator matrix can be arbitrarily close to zero and thus inverse of its minimum eigenvalue may not even have finite moments.

To get around this, I modify the argument by considering a “data-dependent” choice of approximation parameter \(\gamma_n\). This choice of approximation parameter inversely scales with the determinant of the denominator matrix and thus, since the determinant is the product of the eigenvalues, inversely scales with the minimum eigenvalue.\footnote{The determinant has the benefit of being a smooth function of elements of the matrix. This makes it nicer to work with than the minimum eigenvalue itself, which loses differentiability when the dimension of its eigenspace is larger than one.} Geometrically, this approach can be thought of as “stretching out” the function \(\varphi_{\gamma_n}(\cdot)\) in directions where the minimum eigenvalue of the denominator matrix is close to zero. Through the chain rule, this allows for control of the overall derivative of \(\varphi_{\gamma_n}(\JK_I(\beta_0))\) with respective to an individual observation.

High Level Assumptions

I now detail the assumptions needed for the argument, starting with the assumptions that are commmon to the cases with \(d_x = 1\) and \(d_x > 1\). In what comes below \(c > 1\) can be considered an arbitrary constant that may be updated upon each use but that does not depend on sample size \(n\).

(ref)(i) requires that the average second moment of the infeasible first-stage estimators be on the same order as the maximum first-stage estimator second moment. This is imposed mainly to rule out hat matrices that are all zeroes or nearly all zeros so that the effective number of observations used to test the null is growing with the sample size. (ref) provides further intuition for this assumption and below discusses how this assumption and (ref)(ii) may be verified in practice. (ref) compares this balanced design assumption to that in the many-instruments literature crudu_mellace_sandor_2021,Mikusheva_Sun_2022,Matsushita_Otsu_2022,lim2022conditional, noting that their balanced design neither implies nor is implied by the one in this paper.

(ref)(ii) requires that the maximum leverage of any observation be bounded. When \(H\) is symmetric, it is automatically satisfied.\footnote{To see this for \(d_x = 1\) notice that \(s_n^{-2} = \max_i\E[(\widehat\Pi_i^I)^2] \geq \max_i \Var(\widehat\Pi_i^I) = \max_i \sum_{j\neq i} h_{ij}^2 \Var(r_j)\), while \(\Var(r_j)\) will be assumed bounded from below by \(c^{-1}\). Inverting this chain of inequalities yields that \(s_n^2 \sum_{j\neq i} h_{ij}^2\) is bounded from above uniformly over all \(i \in [n]\).} (ref)(iii) can be viewed as a mild technical requirement that there be more than one “effective” instrument in the hat matrix.\footnote{In the case of a standard projection matrix (no deleted diagonal), (ref)(iii) would be satisfied whenever \(\rank(z(z'z)^{-1}z) > 1\), which occurs whenever there are at least two linearly independent instruments.} This condition can be easily verified in practice by examining the eigenvalues of \(HH'\).

Next, I make a high level assumption that the estimation error in \(\hat \rho(\cdot)\) can be treated as negligible in both the numerator and denominator. Later, I will verify this assumption for the particular choice of \(\hat \rho(\cdot)\) described in (ref). For each \(\ell \in [d_x]\) define \(\widehat\Pi^I_{\ell,i} \coloneqq \sum_{j\neq i} h_{ij}r_{\ell j}\), the version of the first stage estimates that could be constructed if \(\rho(\cdot)\) was known to the researcher. Using these, define the magnitude of estimation error in the numerator and denominator as

Showing that \((\Delta_N, \Delta_D)\to_p 0\) implies that estimation error can be treated as negligible for the test statistic \(\JK(\beta_0)\) requires some care. In a standard approach, this would straightforwardly follow from application of the continous mapping theorem. However, this approach requires that the scaled numerator and denominator each have well defined distributional limits, something that is not required by the direct Gaussian approximation. Instead, I establish and make use of anticoncentration bounds to show that the basic results of the continuous mapping theorem still apply even when the numerator and denominator do not have weak limits.

Finally, in addition to characterizing the limiting distribution of \(\JK(\beta_0)\) under \(H_0\), I also examine the behavior of \(\JK(\beta_0)\) in local neighborhoods of the null. These local neighborhoods are characterized by the local power index \(P\), defined below, as well as an additional regularity condition that restricts the size of \(\E[\eps_i(\beta_0)]\) relative to \(\E[r_{\ell i}]\).

Under \(H_0\), (ref) is trivially satisfied since \((\beta - \beta_0) = 0\) and \(\sum_{j\neq i} s_{\ell,n}^2h_{ji}^2 \leq c\) for each \(\ell \in [d_x]\). The local power index is the second moment of the scaled numerator of the test statistic is a measure of the association between the true first stage \(\Pi_i\) and the first-stage estimates \(\widehat\Pi_i\). In (ref), I discuss how the strength of this association is related to the power of the test under local alternatives. (ref)(ii) can be roughly interpreted as requiring the local neighborhoods of \(H_0\) considered to be those in which the means of \((\eps_1(\beta_0),\dots,\eps_n(\beta_0))\) are of the same or lesser order than the means of \((r_1,\dots,r_n)\).

Limiting Behavior of the Test Statistic

Having detailed the assumptions needed for the argument, I now present results showing that, in local neighborhoods of \(H_0\), the distribution of the test statistic can be uniformly approximated by the distribution of the Gaussian statistic described in (ref). As mentioned, the argument can be simplified to require lighter moment bounds when \(d_x = 1\). These moment bounds will be made on the model primitives, \(\eta_i \coloneqq (\beta - \beta_0)v_i + \eps_i = \eps_i(\beta_0) - \E[\eps_i(\beta_0)]\) and \(\zeta_{\ell i} \coloneqq v_{\ell i} - \rho_\ell(z_i)\eta_i = r_{\ell i} - \E[r_{\ell i}]\) for each \(\ell \in [d_x]\).

In the case of \(d_x = 1\), I additionally show that the test based on the \(\JK_I(\beta_0)\) statistic is consistent whenever the power index diverges, \(P\to\infty\), and (ref)(ii) holds.

The dependence of the consistency result on (ref)(ii) is a nontrivial restriction because of the bias taken on in constructing \(r_i\). In particular, against certain alternatives it is possible that \(\E[\widehat\Pi_i^I] = 0\) for all \(i \in [n]\) even under strong identification. This is an extreme case, however. In general, bias in \(\E[r_i]\) does not imply a violation of (ref)(ii), which requires only that the size of \(\E[r_i]\) be of a weakly greater order than that of \(\E[\eps_i(\beta_0)]\). Moreover, as discussed in (ref), (ref) does not necessarily rule out consistency when \(P \to \infty\) but (ref)(ii) fails.

Regardless, bias taken on in constructing \(r_i\) has consequences for the power of the test in finite samples. The combination with the sup-score test described in (ref) is an attempt to rectify this. While this attempt is not perfect, it appears to work well both in the empirical application to the data of gilchrist-glassberg-2016 and in the simulation study of (ref).

The argument when \(d_x > 1\) is considerably more involved than the case where \(d_x = 1\) and requires strengthened moment condition on the variables \(\eta_i\) and \(\zeta_i\). Given a random variable \(X\) and \(\upsilon > 0\) the Orlicz (quasi-)norm is defined \[ \|X\|_{\psi_\upsilon} \coloneqq \inf\{t > 0: \E\exp(|X|^\upsilon/t^\upsilon) \leq 2\} \] Random variables with a finite Orlicz norm for some \(\upsilon \in (0,1] \cup\{2\}\) are termed \(\alpha\)-sub-exponential random variables Gotze-subexponential-polynomial,sambale2022notes. This class of encompasses a wide range of potential distributions including all bounded and sub-Gaussian random random variables (with \(\upsilon = 2\)), all sub-exponential random variables such as Poisson or noncentral \(\chi^2\) random variables (with \(\upsilon = 1\)), as well as random variables with “fatter” tails such as Weibull distributed random variables with shape parameter \(\upsilon \in (0,1]\). Thus, while imposing that the variables \(\eta_i\) and \(\zeta_i\) are \(\alpha\)-sub-exponential is notably stronger than the finite sixth moments required by (ref), it may still be plausible in a wide range of empirical settings.

While \(\JK_G(\beta_0)\) does not have a fixed distribution, examining its behavior is still tractable and allows for insight into the power properties of the jackknife K-test. In (ref), I use this result to analyze the local power of the proposed test. To improve power against certain alternatives, I suggest a combination with the sup-score statistic of BCCH-2012.

Controlling Estimation Error

The final step is to show that estimation error in \(\hat \rho\) can be treated as negligible, that is verify (ref). Establishing that this holds even when identification is weak makes use of the fact that estimation error in \(\hat \rho(\cdot)\) enters the test statistic only through its interaction with the implied error \(\eps_i(\beta_0)\). When identification is weak, the implied error is nearly conditional mean independent of the instruments and thus the product, \(\eps_i(\beta_0)\hat\rho(\cdot)\), is insensitive to small pertubations in \(\hat\rho(\cdot)\).\footnote{In the langauge of CCDDHNR-2018, this is termed “Neyman Orthogonality.” Under strong identification, Neyman orthogonality allows for the use of many machine learning techniques to estimate the first stage. This (approximate) orthogonality also holds under structural parameter sequences that are local to the null.}

In theory, this orthogonality could be combined with cross-fitting to allow for the use of other machine learning methods to estimate \(\hat\rho(\cdot)\), as in CCDDHNR-2018. This possibility is explored briefly in (ref). The use of other machine learning methods may be useful if sparsity of \(\rho(\cdot)\) is not a plausible assumption in the researcher's empirical setting as estimators such as random forests and neural networks can be consistent in high dimensional settings under alternate assumptions.\footnote{chi2022asymptotic provide results for random forests under a “sufficient impurity decrease” assumption while farrellNeuralNetwork2021 and Schmidt-Hieber-2020-neural-networks provide results for neural networks under a generalized heirarchical model assumption.} However, I focus on the \(\ell_1\)-penalized procedure proposed in (ref) both for expositional simplicity and because the sparsity assumption required for the consistency of this procedure mirrors that needed for the popular post-Lasso first-stage estimator.

(ref)(i) ensures that \(\eta_i\) has finite exponential moments (i.e has a well defined moment generating function), which is required to allow the number of basis terms used to approximate \(\rho(\cdot)\) to grow at a near exponential rate compared to the sample size (\(d_b \gg n\)). When using fewer basis terms this assumption may be relaxed. (ref)(ii) is a standard condition in \(\ell_1\)-penalized estimation. At the cost of extra notation, it can be relaxed and the sup-norm of the basis terms can be allowed to grow slowly with the sample size to accommodate bases such as normalized b-splines or wavelets. (ref)(iii) is a bound on the rate of decay of the approximation error, similar to the approximate sparsity condition of BCCH-2012.

(ref)(iv) is a high-level condition on the rate of consistency of the parameter estimate \(\hat\phi\) in the \(\ell_1\) norm. This can be verified under approximate sparsity for both the LASSO estimator in (ref) or post-LASSO procedures based on refitting an unpenalized version of (ref) only using the basis terms selected in a LASSO first stage. See BCCH-2012,vanDerGreer2016,Tan-2017, and CS-2021 for references under various choices of penalty parameter. This condition allows for the dimensionality of the basis terms, \(d_b\), to grow near exponentially as a function of the sample size. Following the analysis of Tan-2017 this may be satisfied as long as \(s^2\log^{2(\upsilon + 1)/\upsilon}(d_b n)/n\to 0\), where the sparsity index \(s\) denotes the number of nonzero elements of \(\phi\).

(ref)(v) is a strengthening of the definition of local neighborhoods and can be interpreted similarly to (ref)(ii). Since the moment conditions in (ref)(va,vb) hold with \(b_\ell(z_j)\eps_j(\beta_0)\) replaced with \(r_j\), (ref)(v) can be interpreted as requiring that \(|\E[\sum_{j\neq i} h_{ij}b_\ell(z_j)\eps_j(\beta_0)]|\) is on the same order as \(|\E[\sum_{j\neq i} h_{ij}r_j]|\) for all \(i =1,\dots,n\) and \(\ell = 1,\dots,d_b\). As with (ref)(ii), it is trivially satisfied under \(H_0\) or, using the fact that \(\max_i \sum_{j\neq i} s_n^2 h_{ij}^2 \leq c\), whenever \(\E[\eps_i(\beta_0)] = \Pi_i(\beta - \beta_0)\) is in a \(\sqrt{n}\)-neighborhood of zero.

Local Power and Improvements

Using the characterization of the limiting behavior of the test statistic derived in (ref), I analyze the local power properties of the test. Unfortunately, against certain alternatives the test statistic may have trivial power, a deficiency shared with the K-statistics of Kleibergen_2002,Kleibergen-2005. Finally, I describe how the the combination with the sup-score statistic, described in (ref), attempts to remedy this and formally establish its validity.

Local Power Properties

For exposition, I focus on the case where \(d_x = 1\). In local neighborhoods of \(H_0\), as defined in (ref), (ref) implies that the limiting behavior of \(\JK(\beta_0)\) can be analyzed by examining the behavior of the Gaussian analog statistic, \(\JK_G(\beta_0)\). Conditional on the vector \(\tilde r = (\tilde r_1,\dots,\tilde r_n)\), the distribution of \(\JK_G(\beta_0)\) is nearly non-central \(\chi^2_1\) with noncentrality parameter \(\mu(\tilde r)\), \(\JK_G(\beta_0) |\tilde r \sim A^2(\tilde r) \cdot \chi^2_1(\mu(\tilde r))\):

Under local alternatives, the terms \(\Pi_i^2(\beta - \beta_0)^2 \to 0\) so that \(A(\tilde r) \to 1\) and \(|\mu^2(\tilde r) - \mu_\infty^2(\tilde r)|\to 0\), where

The numerator of \(\mu^2_\infty(\tilde r)\) suggests that power is maximized when the first-stage estimate \(\tilde \Pi_i\) is close to the true first stage value \(\Pi_i\). Indeed, when errors are homoskedastic \(\mu_\infty^2(\tilde r)\) is maximized by setting \(\tilde \Pi_i = \Pi_i\) reflecting the classical result of Chamberlain-1987. The denominator of \(\mu^2_\infty(\tilde r)\) suggests that having first-stage estimates \(\tilde \Pi_i\) with low second moments may increase power. This guides the recommendation for the use of \(\ell_2\)-regularization in constructing the hat matrix, \(H\).

Unfortunately, estimators of \(\Pi_i\) based on \(r_i = x_i - \rho(z_i)\eps_i(\beta_0)\) may not be close to \(\Pi_i\) under \(H_1\). This is because the mean of \(r_i\) will in general differ from \(\Pi_i\) \[ \E[r_i] = \Pi_i - \rho(z_i)\Pi_i(\beta - \beta_0) \] This deficiency is inherited from the similarity of the \(\JK(\beta_0)\) statistic to the K-statistic. As pointed out by Moreira-2001, this need not be an issue as long as there is a fixed constant \(C\neq 0\) such that \(\E[r_i] = C\Pi_i\) for all \(i \in [n]\). However, in general, this will introduce bias into the first-stage estimates \(\widehat\Pi_i\) under \(H_1\). The power implications of this bias are particularly pronounced when \(\rho(z_i)\) is a constant \((\beta - \beta_0) = 1/\rho(z_i)\). In this case, \(\E[r_i]\), and thus \(\E[\tilde\Pi_i]\), will equal zero for each \(i\in[n]\), and the \(\JK(\beta_0)\) statistic will select a direction completely at random to direct power into.\footnote{AMS-2006 and Andrews-2016 point out this deficiency in the context of the K-statistics of Kleibergen_2002,Kleibergen-2005.}

A Simple Combination Test

To combat this loss of power for tests based on the K-statistic, a common strategy is to combine the K-statistic with the Anderson-Rubin statistic based on a conditioning statistic. While the Anderson-Rubin statistic does not have optimal power on its own, it has the benefit of directing power equally in all directions avoiding the pitfalls of the K-statistic which lacks power in certain directions. Prominent examples of such tests are the conditional likelihood ratio test of Moreira_2003, the GMM-M test of Kleibergen-2005, and the minimax regret tests of Andrews-2016. These combinations make use of the fact that the Anderson-Rubin statistic is asymptotically independent of both the K-statistic and the conditioning statistic.

Unfortunately, the asymptotic validity of these tests under heteroskedasticity is based on the assumption that \(d_z^3/n\to 0\), which may not reasonably describe many settings discussed above. Instead, to improve the power of tests based on the jackknife K-statistic, I consider a simple combination with the sup-score statistic of BCCH-2012. The test based on the sup-score statistic (ref) is similar in spirit to the Anderson-Rubin test but controls size even when \(d_z\) grows near exponentially as a function of the sample size.

A size \(\theta \in (0,1)\) test based on the sup-score statistic rejects whenever \(S(\beta_0) > c^S_{1-\theta}\) where, for \(e_1,\dots,e_n\) iid standard normal and generated independently of the data, \(c^S_{1-\theta}\) is the simulated multiplier bootstrap critical value: \[ c^S_{1-\theta} \coloneqq (1-\theta)\text{ quantile of } \sup_{1 \leq \ell \leq d_z} \bigg|\frac{\sum_{i=1}^n e_i \eps_i(\beta_0)z_{\ell i}}{(\sum_{i=1}^n z_{\ell i}^2)^{1/2}}\bigg|\text{ conditional on }\{(y_i,x_i,z_i)\}_{i=1}^n. \] As with the Anderson-Rubin test, tests based on the sup-score statistic may have suboptimal power properties in overidentified models as it does not incorporate first-stage information. However, the sup-score statistic does retain the benefit of directing power evenly in all directions, avoiding pitfalls of tests based on \(\JK(\beta_0)\) against certain alternatives.

The combination test will be based on an attempt to detect whether the alternative \(\beta\) is such that \(\E[\widehat\Pi_{\ell,i}^I] = 0\) for all \(i = 1,\dots, n\) and some \(\ell \in [d_x]\). When this is the case, the researcher would prefer to test the null hypothesis using the sup-score statistic. As mentioned in (ref), detection of whether \(\E[\widehat\Pi_{\ell,i}^I] = 0\) for some \(i \in [n]\) is based on the conditioning statistic:

Under the assumption that \(\E[\widehat\Pi_i^I] = 0\) for all \(i \in [n]\), quantiles of the conditioning statistic can be simulated analogously to the sup-score critical value. For a new set of \(e_1,\dots,e_n\) iid standard normal and generated independently of the data, and for any \(\theta\in (0,1)\), define the conditional quantile

Depending on the value of the conditioning statistic, the thresholding test decides whether the test based on \(\JK(\beta_0)\) or one based on \(S(\beta_0)\) should be run.

for some cutoff \(\tau\), which I take in the simulation study and empirical exercise to be the 75\textsuperscript{th} quantile of the distribution of \(C\) under the assumption that \(\E[\widehat\Pi_i^I] = 0, \forall i \in [n]\).

To show that the thresholding test controls size, I compare the rejection probability to that of a Gaussian analog. In addition to \(\JK_G(\beta_0)\), defined in (ref), define the Gaussian analogs of \(\S(\beta_0)\) and the conditioning statistic \(C\):

where, as in (ref), \((\tilde\eps_i(\beta_0),\tilde r_i)'\) are generated independently of each other and the data following a Gaussian distribution with the same mean and covariance matrix as \((\eps_i(\beta_0),r_i)\). Since \(\Cov(\tilde\eps_i(\beta_0),\tilde r_i) = 0\) under \(H_0\), the statistics \(C_G\) and \(S_G(\beta_0)\) are independent under the null. Similarly, the null distribution of \(\JK_G(\beta_0)\) is the same conditional on any realization of \((\tilde r_1,\dots,r_n)\); it is also independent of \(C_G\) under the null. The Gaussian analog thresholding test decides whether the researcher should run a test based on \(S_G(\beta_0)\) or \(\JK_G(\beta_0)\) depending on the value of \(C_G\) as in (ref).

The test statistics \(\JK_G(\beta_0)\) and \(S_G(\beta_0)\) are only marginally independent of the conditioning statistic \(C_G\) under the null. This limits the ways in which the test statistics can be combined using the conditioning statistic while still controlling size. This marginal independence in the Gaussian limit is enough, however, for the asymptotic validity of the thresholding test, \(T(\beta_0;\tau)\). To establish that the behavior of the pairs \((C,\JK(\beta_0))\) and \((C,S(\beta_0))\) can be approximated by the behavior of \((C_G,\JK_G(\beta_0))\) and \((C_G,S_G(\beta_0))\), respectively, I rely on the following assumption:

(ref)(i) is a strengthening of the moment bound on \(r_i\) similar to that of (ref)(i). As discussed, while more restrictive than the condition in (ref), this still allows for a wide range of potential distributions for \(r_i\). (ref)(ii) requires that the number of observations used to test \(\E[\widehat\Pi_i] = 0\) via the conditioning statistic and the number of observations used to test the null hypothesis via the sup-score test are both growing with the sample size. It can be verified by looking at the hat matrix \(H\) and the instruments. Finally, (ref)(iii) is a light requirement on the number of instruments \(d_z\) needed for the validity of the sup-score test. It allows the number of instruments to grow near exponentially as a function of sample size.

The first part of (ref) establishes the asymptotic validity of the thresholding test \(T(\beta_0;\tau)\) for any choice of cutoff \(\tau\). While not explicitly stated in the statment of the theorem, this result is uniform in the choice of \(\tau\); for any sequence \(\{\tau_n\}\subset \SR_+\) the sequnce of testing procedures \(T(\beta_0;\tau_n)\) will also have asymptotic size \(\alpha\). The proof of this statement follows the logic outlined above. The second part of (ref) establishes the validity of the multiplier bootstrap procedure to approximate quantiles of the conditioning statistic. It follows directly from results in belloni2018highdimensional after verifying that the conditions needed for error taken on from estimation of \(\rho(z_i)\) can treated as negligible under (ref).

In the case of a single endogenous variable, \(d_x = 1\), (ref) could be established under the lighter conditions of (ref) along with (ref). However, for brevity, I do not seperate the two cases here.

Conclusion

I propose a new test for the structural parameter in a linear instrumental variables model. This test is based on a jackknife version of the K-statistic and the limiting behavior of the test is analyzed via a novel direct Gaussian approximation argument. I show that, as long as an auxiliary parameter can be consistently estimated, the test is robust to both the strength of identification and the number of instruments; the limiting distribution of the test statistic does not depend on either of these factors. Consistency of the auxiliary parameter can be achieved under approximate sparsity using simple-to-implement \(\ell_1\)-penalized methods.

I characterize the behavior of the jackknife K-statistic in local neighborhoods of the null. To address a power deficiency that tests based on jackknife K-statistic inherit from their non-jackknife namesakes, I propose a testing procedure that decides whether the researcher should run a test via the jackknife K-statistic or one via the sup-score statistic based on the value of a conditioning statistic. While this combination may not fully address the power decline, I show that it works well in a simulation study and leave further refinements to future work.

\startappendix

Proof of (ref)

(ref) follows from the following two main technical lemmas, the proofs of which will comprise the majority of this appendix section. Let \(\JK_I(\beta_0)\) be the version of the test statistic that could be constructed if \(\rho(\cdot)\) was known to the researcher, defined in more detail shortly.

Proof of (ref)

Before proceeding, we will introduce some notation. Let \(\tilde H = s_n H\) and \(\tilde h_{ij} = s_n h_{ij}\), where \(s_n\) is as in (ref). Recall that \(\tilde h_{ii} = 0\) and define

where \((\tilde\eps_i(\beta_0),\tilde r_i)\) are jointly Gaussian with the same mean and covariance matrix as \((\eps_i(\beta_0),r_i)\) and \(\kappa_i^2(\beta_0) = \E[\eps_i^2(\beta_0)]\). Under this notation we can write \(\JK_I(\beta_0) = \frac{N^2}{D}\bm{1}_{\{D > 0\}}\) and \(\JK_G(\beta_0) = \frac{\tilde N^2}{\tilde D}\). Dealing with these forms of the statistics is difficult for the interpolation argument, since the denominator is random. Instead, we will notice that since \(D = 0 \implies N = 0\) and \(\Pr(\tilde D > 0) = 1\), for any \(a \geq 0\) we can rewrite the events

With this in mind define

Showing (ref) is then equivalent to showing that \(\sup_a|\Pr(\JK^a \leq 0) - \Pr(\tilde \JK^a \leq 0)|\to 0\). The statement \(\sup_{a < 0}\big|\Pr(\JK_I(\beta_0) \leq a) - \Pr(\JK_G(\beta_0) \leq a)\big| = 0\) is immediate since both \(\JK_I(\beta_0)\) and \(\JK_G(\beta_0)\) are always weakly positive. It thus suffices to show \[ \sup_{a \geq 0} \big|\Pr(\JK_I(\beta_0) \leq a) - \Pr(\JK_G(\beta_0) \leq a)\big| \to 0 \] We do so in a few lemmas, the final result being shown in (ref) at the bottom of this subsection.

Proof of (ref)

Proof of (ref)

Proof of (ref)

I provide the proof for the case where \(d_x = 1\) with the general case following from symmetric logic. For any \(j = 1,\dots, d_b\) define the matrix \(B_j = \diag(b_j(z_1),\dots,b_j(z_n))\) and collect observations \(\eps(\beta_0) = (\eps_1(\beta_0),\dots,\eps_n(\beta_0))'\in \SR^n\), \(r = (r_1,\dots,r_n)'\in\SR^n\), \(\hat r = (\hat r_1,\dots, \hat r_n)'\in\SR^n\), and \(\xi = (\xi_1,\dots,\xi_n)' \in \SR^n\). In addition, collect \(b_\eps = (b_{\eps1},\dots,b_{\eps n})\in \SR^{d_b\times n}\) where \(b_{\eps i} = \eps_i(\beta_0)b(z_i) \in \SR^{d_b}\). Finally, let \(\mH = \frac{s_n}{\sqrt{n}}H\), \(\tilde H = s_n H\) and \(\tilde h_{ij} = s_n h_{ij}\).

Step 1: \(\Delta_N\to_p 0\). To show that \(\Delta_N \to_p 0\) write

To bound \(\vA\) we move to apply (ref) to the quadratic form \(\eps(\beta_0)'(\mH B_j)\eps(\beta_0)\). First notice that, under (ref)(v), we have \[ \|\E[\mH b_j\eps(\beta_0)]\|_2 = \frac{1}{n}\sum_{i=1}^n (\E[s_n\sum_{j\neq i} h_{ij} b(z_j)\eps_j(\beta_0)])^2 \leq c^2 \] In the notation of (ref) this give us an upper bound on \(\|\E f^{(1)}(X)\|_{\text{HS}}\). Next, (ref) gives us that the Frobenius norm of \(\mH = \frac{s_n}{\sqrt{n}}\mH\) is bounded, since the rows of \(s_n H\) are square summable, \(\sum_{j\neq i} (s_n h_{ij})^2 \leq c\) for all \(i = 1,\dots, n\). In the notation of (ref) this gives us an upper bound on \(\|\E f^{(2)}(X)\|_{\text{HS}}\). Applying (ref) and a union bound then gives us that

Since \(\max_{1 \leq j \leq d_b}|\E[\eps(\beta_0)'\mH B_j\eps(\beta_0)]| \leq c\) under (ref)(v), (ref) gives that \[ \max_{1 \leq j \leq d_b} |\eps(\beta_0)'\mH B_j \eps(\beta_0)| = O_p(\log^{2/a}(d_b)) \] Since \(\log^{2/a}(d_b)\|\widehat\gamma - \gamma\|_1 \to_p 0\) by assumption, this yields that \(\vA \to_p 0\).

To bound \(\vB\) see that \(\|\eps(\beta_0)'\mH\|_2 = \frac{s_n^2}{n}\sum_{i=1}^n (\sum_{j\neq i} h_{ij}\eps_i(\beta_0))^2 = O_p(1)\) under (ref)(ii) while under (ref) \(\|\xi\|_2 = o(1)\).

Step 2: \(\Delta_D\to_p 0\). Notice that \(a^2 - b^2 = 2b(a - b) + (a - b)^2\) and bound:

Since both \(\vE = O_p(1)\) and \(\vF = O_p(1)\) under the moment bounds of (ref) and (ref), it suffices to show that \[ \max_i |\sum_{j\neq i} \tilde h_{ij} (\widehat r_j - r_j)| \to_p 0 \] To do so write

To bound \(\vA\), note that by (ref)(v) \(\max_{i,j}|\E[\sum_{j\neq i} \tilde h_{ij} b(z_j) \eps_j(\beta_0)| \leq c\). Under (ref)(ii), \(\max_{i,j} \sum_{j\neq i} \tilde h_{ij}^2 b^2(z_j) \leq c^2\) so we can apply (ref) and a union bound to obtain that \[ \max_{\substack{1 \leq i \leq n \\ 1 \leq j \leq d_b}} \big|\sum_{j\neq i} \tilde h_{ij} b(z_j)\eps_j(\beta_0)\big| = O_p(\log^{1/a}(d_bn)) \] Along with the implied rate on \(\|\hat\gamma - \gamma\|_1\) from (ref)(iv) this shows that \(\vA \to_p 0\).

To show that \(\vB \to 0\), use Cauchy-Schwarz, \(\sum_{j\neq i} \tilde h_{ij}^2b^2(z_j) \leq c\) for any \(i,j\) by (ref)(ii), and \(\sum_{i=1}^n \xi_i^2 = o(1)\) by (ref)(iii).

Proof of (ref)

Throughout this section, define the scaled elements of the infeasible and gaussian numerators and denominators

Collect these in \(N = (N_1,\dots N_{d_x})' \in \SR^{d_x}\), \(\tilde N = (\tilde N_1,\dots, \tilde N_{d_x})' \in \SR^{d_x}\), \(D = [D_{\ell k}]_{\ell, k \in [d_x]} \in \SR^{d_x \times d_x}\), and \(\tilde D = [\tilde D_{\ell k}]_{\ell, k \in [d_x]} \in \SR^{d_x \times d_x}\). After multiplying by scaling matrix \(\diag(s_{1,n},\dots,s_{d_x,n})\) and the inverse of the scaling matrix we rewrite the infeasible and gaussian test statistics

As with (ref), the result (ref) follows directly from combining the following lemmas. The first is the main technical lemmma, and shows that the distribution of the infeasible statistic \(\JK_I(\beta_0)\) can be uniformly approximated by that of \(\JK_G(\beta_0)\). The proof of this technical lemma is involved and deferred to (ref). The second lemma establishes that estimation error can be treated as negligible. As with (ref), the main difficulty hear is in dealing with the fact that neither the numerator vector nor denominator matrix of the \(\JK(\beta_0)\) statistic may have stable limiting distributions.

Proofs of Results in (ref)

The statement of (ref) relies on showing

In particular, since \((\JK_G(\beta_0) \perp C_G)\) and \((\S_G(\beta_0)\perp C_G)\) under \(H_0\), showing the above will imply the test based on \(T(\beta_0;\tau)\) has asymptotic size \(\alpha\) for any choice of cutoff \(\tau\). The second line in the above display follows imediately from (ref) after verifying (ref), below.

The first line in the top display relies on a joint interpolation of the infeasible \(\JK_I(\beta_0)\) test statistic and the infeasible conditioning statistic \(C_I\), which could be constructed if \(\rho(z_i)\) was known to the researcher.

This joint interpolation argument is rather involved however, and deferred to (ref). The interpolation argument for the conditioning statistic very closely follows the results in cck2013. The results of (ref) rely on showing that the difference between \(C\) and \(C_I\) can be treated as negligible. This in turn reduces to verifying (ref), which is done in (ref), below.

Proof of (ref)

The first result in (ref) with \(\JK(\beta_0)\) and \(C\) replaced with their infeasible analogs \(\JK_I(\beta_0)\) and \(C_I\) follows from the argument in (ref). After verifying that \(|\JK(\beta_0) - \JK(\beta_0)| \to_p 0\) via (ref) and that (ref) is satisfied via (ref) follow the same steps as in the proof of belloni2018highdimensional, {Theorem 2.1} to see that approximation result holds for the feasible \(\JK(\beta_0)\) and \(C\).

For the second statement, I show that the conditions of (ref) are satisfied. To see that (ref)(i,ii) is satisfied under the moment assumptions of (ref) use (i) the definition of \(\dot h_{ij} = h_{ij}\big/(n^{-1}\sum_{j=1}^n h_{ij}^2)^{1/2}\); (ii) that the variance of each \(r_{j}\) is bounded away from zero and (iii) that the fourth moments of \(r_j\) are bounded from above. (ref)(iii) is satisfied with \(B_n = \log^{1/\upsilon}(n)\) by (ref)(i,iii) and (ref). Finally (ref) is satisfied by applying (ref). Apply (ref) to conclude.

Joint Gaussian Approximation of \(\JK(\beta_0)\) and \(C\)

(ref) rely on a joint interpolation of the conditioning and testing statistics as well as a joint interpolation of the conditioning and testing statistics. The joint interpolation of \(\JK(\beta_0)\) and the conditioning statistic \(C\) is given in (ref) after introducing some notation in (ref). The joint gaussian approximation of \(S(\beta_0)\) and \(C\) follows immediately from results in \cites{belloni2018highdimensional,CCK-2018-HDCLT}. The result is presented below for the general form of the \(\JK(\beta_0)\) statistic under \(H_0\) however the proof strategy is very similar when using the decomposed form of \(\JK(\beta_0)\) when \(d_x = 1\). This proof is available on request.

Notation

\paragraph{Jackknife Statistic Definitions.} Define \(\tilde h_{\ell, ij} = s_{n,\ell} h_{ij}\) for each \(\ell = 1,\dots,d_x\) and the scaled leave-one-out quasi-numerator and denominators

where \(\dot\eps_j(\beta_0)\) is equal to \(\tilde\eps_j(\beta_0)\) if \(j < i\) and equal to \(\eps_j(\beta_0)\) if \(j > i\), \(\dot r_{\ell j}\) is equal to \(\tilde r_{\ell j}\) if \(j < i\) and equal to \(r_j\) if \(j > i\), and \(\ddot \eps_j(\beta_0)\) is equal to \(\E[\eps_j^2(\beta_0)]\) if \(j < i\) and equal to \(\eps_j(\beta_0)\) if \(j > i\). As in the proof of (ref) while the definitions of \(\dot\eps_j(\beta_0), \dot r_{\ell j}\), and \(\ddot\eps_j(\beta_0)\) depend on \(i\) this dependence is suppressed to conslidate notation and since we only consider one step deviations at a time.

Also define the one step deviations

where

Notice that in this notation we can write the test statistic and gaussian test statistics, after scaling by \(\diag(s_{n,1},\dots,s_{n,d_x})\), as

In this proof we will use these representations for the test statistics. Finally define

\paragraph{Conditioning Statistic Definitions.} Let \(h_{\ell, ii} = 0\) for any \(\ell = 1,\dots,d_x\) and \(i= 1,\dots,n\). Define \(\tilde h_{\ell,ij} = h_{\ell,ij}/\omega_{\ell i}\) for \(\omega_{\ell i} = n^{-1}\sum_{j=1}^n |h_{\ell,ij}|\). Also define the one-step deviations:

Notice that \(C = \max_{1 \leq \iota \leq 2nd_x} (C_{-1} + \frac{1}{\sqrt{n}}\Delta_{C1})_\iota\) while \(\tilde C = \max_{1 \leq \iota \leq 2nd_x}(C_{-n} + \Delta_{Cn})_\iota \).

\paragraph{Function Definitions.} As in cck2013 consider the “smooth max” function, \(F_\beta:\SR^p \to \SR\) defined \[ F_\beta(z) = \beta^{-1}\log\bigg(\sum_{i=1}^n \exp(\beta z_i)\bigg) \] which satisfies \[ 0 \leq F_\beta(z) - \max_{1 \leq i \leq n} z_i \leq \beta^{-1}\log p. \] (ref) notes some useful properties of the smooth max function which we will use in the joint interpolation argument. In addition let \(\varphi(\cdot) \in C_b^3(\SR)\) be such that \(\varphi(x) = 1\) if \(x \leq 0\), \(\varphi'(x) < 0\) for \(x \in (0,1)\), and \(\varphi(x) = 0\) for \(x \geq 1\). For any \(\gamma > 0\) and \(a = (a_1,a_2)'\in \SR^2\) define the function \(\tilde \varphi(\cdot,\cdot,\cdot): \SR^{d_x} \times \vec(\SR^{d_x \times d_x}) \times \SR^{2nd_x} \to \SR\) via

where

The function \(\tilde\varphi_{\gamma,a}(\cdot,\cdot,\cdot)\) is meant to approximate the indicator function \(\bm{1}\{K(\beta_0) \leq a_1\}\bm{1}\{C \leq a_2\}\) with \(\gamma\) governing the quality of approximation. Where it is obvious, we will supress the subscripts \(\gamma,a\) from our notation.

Main Argument

\newblankpage

Incorporating Exogenous Controls

In this section, I analyze the model with exogeneous controls. To this end, define the vector \(z_2 = (z_{21}',\dots,z_{2n}')' \in \SR^{n \times d_c}\). Let \(P_2 = z_2(z_2'z_2)^{-1}z_2' \in \SR^{n\times n}\) denote the projection onto the column space of \(z_2\) and \(M_2 = I_n - P_2\) denote the projection onto to orthocomplement of the column space. Focus will be on the case where \(d_x = 1\) to simplify notation, but the basic concepts apply generally to \(d_x > 1\).

For \(y \coloneqq (y_1,\dots,y_n)' \in \SR^n\) and \(x \coloneqq (x_1',\dots,x_n')'\in \SR^{n\times}\) define \(y^\perp \coloneqq M_2y\) and \(x^\perp \coloneqq M_2 x\) as the “partialled out” versions of \(y\) and \(x\), respectively. Let \(y^\perp_i\) be the \(i\)\textsuperscript{th} element of \(y^\perp\) and \(x^\perp_i\) be the \(i\)\textsuperscript{th} element of \(x^\perp\). From here we can define \(\eps(\beta_0) \coloneqq y - x\beta_0 \), \(\eps^\perp(\beta_0) = M_2\eps(\beta_0)\) and \(r^\perp \coloneqq M_2r\) where as in the main text \(r = (r_1,\dots,r_n)'\) is constructed \(r_i = x_i - \rho(z_i)\eps_i(\beta_0)\). The definition of \(\rho(z_i)\) does not change after partialling out \(z_2\) since all expectations are understood to be conditional on the instruments \(z\). Notice that \(\eps^\perp(\beta_0)\) is mean zero. Finally I assume that the controls have been partialled out of hat matrix so that the effective hat matrix is \(M_2H\) and the vector \(\widehat\Pi \in \SR^n\) is defined \(\widehat\Pi = (M_2 H)(M_2 r)\). This does not make a difference for the numerator of the \(\JK(\beta_0)\) statistic but does affect the denominator slightly. When this is not done, inference may be conservative.

Using matrix notation in the numerator to make things clear, we can write the version of the \(\JK(\beta_0)\) statistic with the partialled out vectors, \(\eps^\perp(\beta_0)\) and \(r^\perp\), in terms of the original vectors, \(\eps(\beta_0)\) and \(r\),

where \(\vh_{ij} = [M_2 \tilde HM_2]_{ij}\), \(\tilde H = s_n H\), and \(m_{ij} = [M_2]_{ij}\). I seek to characterize the limiting distribution of \(\JK(\beta_0)\) under \(H_0\). To do so, we show that quantiles \(\JK(\beta_0)\) can be approximated by quantiles of the gaussian analog statistic

where \((\tilde\eps_i, \tilde \eps_i(\beta_0), \tilde r_i)\) are generated gaussian independent of the data and with the same mean and covariance as \((\eps_i, \eps_i(\beta_0),r_i)\). Since \(\Var(\tilde\eps(\beta_0)) = \Var(\eps_i)\) under \(H_0\), \(\E[\tilde\eps(\beta_0)'M_2] = 0\), and \(\tilde r \perp \tilde\eps(\beta_0)\), this gaussian analog statistic has a \(\chi^2_1\) distribution conditional on any realization of \(\tilde r\) and thus its unconditional distribution is also \(\chi^2_1\).

Showing that quantiles of \(\JK(\beta_0)\) can be approximated by quantiles of \(\tilde\JK(\beta_0)\) proceeds in two steps. In the first step, we show that \(\JK(\beta_0)\) converges in probability to an intermediate statistic.

We will then show that quantiles of this intermediate statistic can be approximated by quantiles of \(\tilde\JK(\beta_0)\). In view of (ref), it suffices to show for the first step that \(\Delta_D \to_p 0\), where

To do this, notice that under \(H_0\) we can write \(\eps_i^\perp(\beta_0) = \eps_i + z_{2i}'(\widehat\Gamma - \Gamma)\) where \(\widehat\Gamma = (z_2'z_2)^{-1}z_2\eps(\beta_0)\) is a \(\sqrt{n}\)-consistent estimate of \(\Gamma\). Exploiting this fact we get

Both of these terms will tend to zero by the consistency \(\widehat\Gamma\) to \(\Gamma\), giving that \(\Delta_D \to_p 0\).

In our second step, we argue that quantiles of \(\JK^{\text{int}}(\beta_0)\) can be approximated by quantiles of \(\JK_G(\beta_0)\). To make this comparasion, we can follow almost exactly the same steps as in (ref). The only difference between analysis in this case and analysis in the original case is that the partialling out of controls leads the test statistic to not strictly have a jackknife form; the effective hat matrix \(M_2HM_2\) no longer has a deleted diagonal. However, as I will argue below, this will not make a difference in the interpolation argument since the diagonal terms of \([P_2]_{ii}\) are small in the sense that they sum to \(d_c\).

The (ref) analog one step deviations for the numerator are given

where as (ref), a dotted variable is equal to the gaussian analog if \(j > i\) but equal to the standard version otherwise. The first and second moments of the first two terms in \(\Delta_{1i}\) can be matched with their gaussian analog terms as in the proof of (ref). While we cannot match seconds moments of the third term in the one step deviation, this sum of all these third terms can be treated as negligible after scaling by \(1/\sqrt{n}\) as \(\sum_{i=1}^n |\vh_{ii}| \lesssim d_c\). This is because \(M_2\tilde H M_2 = \tilde H - P_2 \tilde H - \tilde H P_2 - P_2 \tilde H P_2\). The matrix \(\tilde H\) has zeros on it's diagonal. Meanwhile

where the final inequality comes because the matrix \(P_2\) is symmetric and idempotent and since \(\Big(\sum_{j\neq i} H_{ji}^2\Big) \lesssim 1\) by (ref)(ii). A similar argument can be used to show that \([P_2 \tilde H P_2]_{ii}^2 \lesssim [P_2]_{ii} \). Since \(P_2\) is a projection matrix we must have that \(\|P_2He_j\| \leq \|He_j\|\) for any basis vector \(e_j \in \SR^n\). Thus \(\sum_{j=1}^n [P_2H]_{ji}^2 \leq \sum_{j=1}^n [H]_{ji}^2\). Finally, we can use the fact that the trace of \(P_2\) is equal to its rank to show that \(\sum_{i=1}^n|\vh_{ii}| \lesssim d_c\)

The one step deviations in the denominator can be bounded using the same logic. These one step deviations are given

where \(\ddot\eps_j\) is equal to \(\Var(\eps_j)\) if \(j < i\) and equal to \(\eps_j\) if \(j > i\). The first three terms in this expansion are can be dealt with exactly as in the proof of (ref). The fourth term is new, however summing over the fourth terms and scaling by \(1/n\) will be negligible as \(\sum_{i=1}^n |\vh_{ii}| \lesssim d_c\). After showing the lindeberg interpolation step, the rest of the proof follows exactly as in (ref).

Alternative Construction of Test Statistic via Cross Fitting

To accomodate a general class of estimators for \(\hat\rho(z_i)\) I propose a cross-fit form of the jackknife K-statistic. In this section I detail the cross-fitting procedure and present high level conditions needed for estimation error to be treated as negligible. These conditions can be satisfied by a large class of machine learning estimators under alternate conditions.

Cross Fit Test Statistic

To construct the cross fit test statistic, evenly (and randomly) split the sample into two subsets, \(\calI_1\) and \(\calI_2\) such that \(\calI_1 \cap \calI_2 = \emptyset\) and \(\calI_1 \cup \calI_2 = [n]\). For each \(k = 1,2\) construct an estimator \(\hat\rho^{(k)}(\cdot)\) using only the observations in \(\calI_k\). For observations in \(i \in \calI_1\) form the first-stage estimates \[ \widehat\Pi_i = \sum_{j \in \calI_1\setminus\{i\}} h_{ij}^{(1)}(x_j - \eps_j(\beta_0)\hat\rho^{(2)}(z_j)) \] where the weights \(h_{ij}^{(1)}\) come from a hat matrix \(H^{(1)}\) that only depends on the observations in \(\calI_1\), for example the ridge regression hat matrix of (ref) using only the instruments \((z_i)_{i\in\calI_1}\). First stage estimates for observations in \(\calI_2\) symmetrically, using the estimator \(\hat\rho^{(1)}(\cdot)\). The test statistic is then constructed as before \[ \JK(\beta_0) = \frac{\big(\sum_{i=1}^n \eps_i(\beta_0)\widehat\Pi_i\big)^2}{\sum_{i=1}^n \eps_i^2(\beta_0)(\widehat\Pi_i)^2} \] Notice that while only half the sample is being used to construct each first stage estimate, the full sample is still being used to test the null hypothesis.

Controlling Estimation Error

After writing the overall hat matrix as \[ H =

, \] analysis of the infeasible statistic proceeds as before. Thus, to show that the crossfit test statistic has a limiting \(\chi^2_1\) distribution by (ref) it suffices to state high-level conditions under which \((\Delta_N, \Delta_D)' \to_p 0\).

Relevant Moment Bounds

Moment Bounds for (ref)

Here I provide some lemmas that are useful in the proof of \Crefrange{lemma:lindeberg-interpolation}{lemma:approximate-distribution}

Useful Properties of Smooth Max

Moment Bounds for (ref)

Matrix Derivative Lemmas

The purpose of this section is largely to establish some matrix derivative expressions that will be useful for the Lindeberg interpolation in

Technical Lemmas

Probability Lemmas

Matrix Lemmas

Miscellaneous Lemmas

The following is a standard bound, but it is used a lot so it is restated here.

Assorted Results from Literature

Concentration Inequalities and Tail Bounds

Anticoncentration Bounds

Let \(\xi \in \SR^n\) follow a normal distribution on \(\SR^n\) with mean zero and covariance matrix \(\Sigma_\xi\). Order the eigenvalues of \(\Sigma_\xi\) in non-increasing order \(\lambda_{1\xi} \geq \lambda_{2\xi} \geq ... \geq \lambda_{n\xi}\). Define the quantities \[ \Lambda_{k\xi}^2 = \sum_{j=k}^\infty \lambda_{j\xi}^2,\;\; k =1,2 \]

We use the following anticoncentration lemma from Nazarov2003 noted in CCK-2018-HDCLT.

Gaussian Comparasions and Approximations

We also use the following gaussian approximation results from \cites{belloni2018highdimensional,CCK-2018-HDCLT}. Let \(X_1,\dots,X_n \in \SR^p\) be independent, mean zero, random vectors and let \(Y_1,\dots,Y_n \in \SR^p\) be independent random vectors such that \(Y_i \sim N(0, \E[X_iX_i'])\). Suppose that the researcher does not directly observe \(X_1,\dots,X_n\) but instead observes noisy estimates \(\widehat X_1,\dots,\widehat X_n \in \SR^p\).

Define the sums

Let \(\calA^{\text{re}}\) be the class of all hyperrectangles in \(\SR^p\); that is, \(\calA^{\text{re}}\) consists of all sets \(A\) of the form

for some \(-\infty \leq a_j \leq b_j \leq \infty\), \(j = 1,\dots, p\). Define \[ \rho_n(\calA^{\text{re}}) \coloneqq \sup_{A \in \calA^{\text{re}}}\big|\Pr(S_n^X \in A) - \Pr(S_n^Y \in A)\big| \] Bounding \(\rho_n(\calA^{\text{re}})\) relies on the following moment conditions:

as well as the following bounds on the estimation error

Let \(e_1,\dots,e_n \overset{\text{iid}}{\sim}N(0,1)\) be generated independently of the data. A gaussian bootstrap draw is defined \[ S_n^{X,\star} \coloneqq \frac{1}{\sqrt{n}}\sum_{i=1}^n e_i \widehat X_i \]