The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
82,296 characters
Uniformly Valid Inference Under Interactive and High-Dimensional Constraints
\doublespacing
\title{Uniformly Valid Inference Under Interactive and High-Dimensional Constraints }
\author[$\dagger$]{Joseph Fry}
\affil[$\dagger$]{\footnotesize
Rutgers University, 75 Hamilton St, 08901, New Brunswick, New Jersey, United States\par\nopagebreak
\textit{E-mail address:} \texttt{[email removed]}}
\date{\today}
\maketitle
\begin{abstract}
Asymptotic normality approximations often fail to hold for extremum estimators when the true value of the parameter is at or close to the boundary of a parameter space. I analyze and develop tests using a quasi-unconstrained estimator, which is asymptotically normal even when the true parameter vector is near or at the boundary. These results generalize previous work with this estimator by allowing for more types of constraints and showing how the method can naturally be modified when a nuisance parameter is also high-dimensional. I show that variations of Wald, Likelihood Ratio, and Lagrange Multiplier tests can control size in a uniform sense, provided the initial constrained estimator is sufficiently accurate. Lastly, I apply the method to an application involving network estimation with panel data.
\end{abstract}
\par \textbf{Keywords: }\textit{Asymptotic Normality, Size Control, Boundary, Spatial Models, High Dimensional}
\section{Introduction} \label{s: Intro}
General asymptotic normality results for unconstrained extremum estimators are well-established in the literature (see, e.g., \cite{Newey-McFadden}). Similarly, when the true value of the parameter is fixed at some interior point of a constrained parameter space, then the constraints are non-binding with probability approaching one and the constrained estimator is asymptotically equivalent to the unconstrained estimator. However, \cite{Andrews1999} and \cite{ConstrainedMEstimation} show that when the parameter is fixed at a boundary point of the parameter space, the asymptotic distribution of the constrained estimator is equal to the distribution of a normal random vector that has been projected onto the tangent cone of the feasible set at that boundary point. The fact that the asymptotic distribution is discontinuous in the parameter results in it providing poor approximations in finite samples when the parameter is close to the boundary of the parameter space relative to the sample size, as can be illustrated by examining the asymptotics with sequences of the parameter drifting towards the boundary. Here, I generalize the results of \cite{KETZ2018}, which builds on earlier work using quadratic approximations by \cite{Andrews1999}, to show that conducting inference using a quasi-unconstrained estimator can allow us to control size in a uniform sense, meaning that the procedure can control size in large samples regardless of how close to the boundary the parameter is.
\par
\cite{KETZ2018} uses what they refer to as a quasi-unconstrained estimator to calculate test statistics. This quasi-unconstrained estimator is calculated by taking some $\sqrt{n}$-consistent constrained estimator and performing a single Newton–Raphson-like iteration. If the true value of the parameter is asymptotically a local minimum, then this quasi-unconstrained estimator will also be $\sqrt{n}$-consistent. Furthermore, since it satisfies a set of first-order conditions asymptotically, it is asymptotically normal by similar reasoning as for an unconstrained extremum estimator. \cite{KETZ2018} then uses this quasi-unconstrained estimator to conduct inference on subvectors using a Conditional Likelihood Ratio test based on the conditioning principle of \cite{Moreira2003}. They show that this test will control asymptotic size in a uniform sense, but only when the parameter space meets specific conditions. They assume the parameter space is a product space, so each element of the parameter vector is constrained to lie in some closed interval (e.g., $\theta_j \in [0,c]$ for some constant $c$ or $\theta_j \in [0, \infty)$). While this covers important applications, such as estimating parameters subject to sign restrictions, it excludes cases where constraints interact. For example, autoregressive models that are restricted to be stationary, the control weights for a synthetic control unit that are constrained to be convex combinations, a random coefficient model with correlated coefficients, and the constraints placed on the network parameters in the spillover model in section \ref{s: Network application}. Using the results of \cite{ANDREWS2020}, I generalize the work of \cite{KETZ2018} to show that the versions of Wald, Likelihood Ratio, and Lagrange Multiplier test statistics that use this quasi-unconstrained estimator can be used to conduct inference that will asymptotically control size in a uniform sense, even when much weaker assumptions are placed on the parameter space.
\par
A substantial literature has developed on inference when the parameter lies on or near the boundary of the parameter space. \cite{TestingOnTheBoundary} provides methods for when the parameters fixed by the null hypothesis may lie on the boundary of the parameter space and there may be a nuisance parameter that appears under the alternative hypothesis, but not under the null. \cite{LRTestOnTheBoundary} give the asymptotic distribution of the Likelihood Ratio test when the parameter is fixed on the boundary of the parameter space. However, in practice, it is generally unknown whether the parameter is on the boundary or close to it, and they do not address the case where the parameter is close to the boundary relative to the sample size. Standard resampling methods for inference also suffer from limitations when the parameter is close to or at the boundary. \cite{BootstrapInconsistency} shows that the standard bootstrap can be inconsistent when the parameter is on the boundary of the parameter space, and \cite{ProblemWithSubampling} show that subsampling and $m$ out of $n$ bootstrap methods do not control asymptotic size in a uniform sense when the parameter may be close to the boundary. Intuitively, this is because the shape of the finite-sample distribution depends on how close the parameter is to the boundary relative to the sample size, which differs when using the full sample versus subsamples. It is possible to modify existing tests to ensure asymptotic size control and/or to increase local asymptotic power by using adaptive critical values, but these modifications typically are computationally expensive (e.g., \cite{MCCLOSKEY2017}) or specialized to specific types of constraints (e.g., \cite{KetzMcCloskey2023} assume only the sign of parameters is constrained and the critical values for the Likelihood Ratio test of \cite{cavaliere2025uniformcriticalvalueslikelihood} and the bootstrap method of \cite{CAVALIERE2022} assume parameters that may be near the boundary lie in a product space). \cite{FAN2023} provide a method that weakens the product space assumption by allowing the parameter space to be defined by a set of linear inequality constraints. The proximal bootstrap method of \cite{LI2025} allows the parameter space not to be a product space. Also, the proximal bootstrap method does not have the limitation of \cite{KETZ2018}'s method that the initial constrained estimator must be $\sqrt{n}$-consistent. On the other hand, they do not apply their method in cases where the nuisance parameter is high-dimensional, and when the orthogonalized version of the method is used, the $\sqrt{n}$-consistent conditions can be weakened for the high-dimensional nuisance parameter.
Neyman orthogonality (introduced by \cite{Neyman1959}) has been used in many recent papers to achieve asymptotic normality for a low-dimensional parameter of interest after estimation of a high-dimensional nuisance parameter, including \cite{ChernozhukovDoubleDebiased},
\cite{ChernozhukovHansenSpindler2015-General}, \cite{ChernozhukovHansenSpindler2015}, \cite{BelloniChernozhukovHansen2012}, \cite{BelloniChernozhukovHansen2014}, and \cite{belloni2018highdimensional}. Once initial estimates of the parameters have been obtained and the moment conditions have been orthogonalized with respect to the high-dimensional nuisance parameter, these papers obtain an estimate of the parameter of interest either by minimizing the orthogonalized moment conditions directly or by applying a one-step correction to an initial estimate using those conditions. This second approach can be viewed as applying a Newton–Raphson step to a General Method of Moments (GMM) criterion function which uses the orthogonalized moments. As a result, it can be viewed as an orthogonalized version of the quasi-unconstrained estimator in \cite{KETZ2018}. In this sense, the orthogonalized quasi-unconstrained estimator is not itself a new approach to estimation. What has gone unnoticed, to the best of my knowledge, is its connection to \cite{KETZ2018} and, more broadly, the benefit this construction provides when the parameter of interest is itself constrained. While some of these applications allow the high-dimensional nuisance parameter to be constrained, they assume the parameter of interest is interior, and so the debiasing step is understood only as a technique for being insensitive to nuisance parameters rather than for handling the boundary. Yet the same reasoning that makes \cite{KETZ2018}'s estimator asymptotically normal at or near the boundary applies to the debiased estimate. As long as the first-stage constrained estimate satisfies certain rate-of-convergence conditions, the constraints on the parameter of interest can be dropped from the second-stage estimate, and the resulting orthogonalized quasi-unconstrained estimator of the parameter of interest will be asymptotically normal even when it is on or close to the boundary of its parameter space.
\par
In section \ref{s: MandE}, I first introduce the quasi-unconstrained estimator and discuss when \cite{KETZ2018}'s conditions for its asymptotic normality will hold. In section \ref{s: Inference}, I show that this estimator can be used to conduct uniformly valid full-vector inference. In section \ref{s: InferenceForSubvectors}, I extend this to methods of conducting uniformly valid subvector inference, including subvector inference with a high-dimensional nuisance parameter. While controlling size is the central focus of the paper, I also discuss challenges in establishing the methods' power properties in section \ref{s: Power}. Lastly, in section \ref{s: Network application}, I show how the method can be applied to the network model of \cite{dePaulaRasulSouza2024}, where panel data is used to simultaneously estimate the spillover effect and the latent network over which spillover occurs. In this setting, constraints on the parameter space can arise both for the identification of the latent network and for the interpretability of the estimated spillover effects. I replicate the tax evasion application of \cite{dePaulaRasulSouza2024} and \cite{BesleyCase1995}, and then run simulations using a network model fitted to this data.\footnote{All data used is from \cite{dePaulaRasulSouza2024}.}
\section{Model and Estimator} \label{s: MandE}
Generally, the parameter to be estimated $\theta$ does not fully specify the distribution of the data, so I use $\omega$ to denote a, usually infinite-dimensional, parameter so that $\gamma := (\theta, \omega)$ fully specifies the distribution of the data. The corresponding parameter space satisfies the conditions given in Assumption \hyperref[A1]{1}.
\textbf{Assumption 1} \phantomsection\label{A1} Let $\Gamma = \{ \gamma = (\theta, \omega) \;|\; \theta \in \Theta,\; \omega \in \Omega(\theta)\}$ be compact, where
\begin{enumerate}
\item $\Theta \subseteq \mathbb{R}^J$ is compact and non-empty.
\item $\Omega(\theta) \subseteq \Omega$ for all $\theta \in \Theta$ for some compact metric space $\Omega$.
\end{enumerate}
Because $\gamma$ specifies the distribution of the data, where I use $n$ to denote the number of observations, it also fully specifies the distribution of the stochastic objective function $\hat{Q}$ being minimized to estimate the value of $\theta$. The true value of $\theta$ is defined to be the unique minimum over $\Theta$ of the function $Q_n$ to which $\hat{Q}$ converges:
\begin{equation}\label{eq: thetaDefined}
\theta_n = \operatorname*{arg\,min}_{\theta \in \Theta} Q_n(\theta)
\end{equation}
I use the subscript $n$ to indicate that I consider drifting sequences of the parameter, so $\theta_n$ and $Q_n$ may change with the sample size. In order to use the method of \cite{KETZ2018}, we first need some initial estimate $\hat{\theta}$, which is obtained by minimizing $\hat{Q}$ over $\Theta$. However, it is not necessary for $\hat{\theta}$ to exactly minimize $\hat{Q}$ over $\Theta$, but rather it has to at least approximately achieve the minimum value. More specifically, I assume that
\begin{equation}\label{eq: ConstrainedDefined}
\hat{Q} (\hat{\theta}) = \inf_{\theta \in \Theta} \hat{Q} (\theta) + o_p(1/n).
\end{equation}
\par
\cite{cox2022} characterizes the asymptotic distribution of a constrained extremum estimator under drifting sequences of the parameter $\{ \theta_n \}_{n \in \mathbb{N}}$ with $\theta_n \in \Theta \subseteq \mathbb{R}^J$ in the case where the sequences of sets $\sqrt{n}(\Theta - \theta_n)$ converge to some set $S \subseteq \mathbb{R}^J$. In this case, they show that the results of \cite{Andrews1999} can be generalized so that the asymptotic distribution of $\sqrt{n}(\hat{\theta} - \theta_n)$ is equal to that of a normal random vector projected onto the set $S$. However, the sequence $\sqrt{n}(\Theta - \theta_n)$ will not converge under all drifting sequences $\{ \theta_n \}_{n \in \mathbb{N}}$ so $\sqrt{n}(\hat{\theta} - \theta_n)$ will not always converge in distribution and what it does converge to is significantly different depending on $S$. Because of this, instead of using $\hat{\theta}$ directly for inference, I focus on using the quasi-unconstrained estimator since it is asymptotically normal under quite general conditions on the drifting sequence $\{ \theta_n \}_{n \in \mathbb{N}}$. As a result, I do not rely on $\sqrt{n}(\Theta - \theta_n)$ converging to a particular set and therefore allow for more general conditions of $\Theta$.\footnote{For example, if $\theta_n$ is fixed at some boundary point $\theta^*$, we do not need to impose Chernoff Regularity on $\Theta$ to ensure that $\sqrt{n}(\Theta - \theta_n)$ is converging to the tangent cone of $\Theta$ at $\theta^*$.}
\par
The quasi-unconstrained estimator is obtained by taking the constrained estimator and performing one Newton-Raphson step using the quadratic approximation of $\hat{Q}$ at $\hat{\theta}$. I use $D Q_n$ and $D^2 Q_n$ to denote the generalized first-order and second-order partial derivatives of $Q_n$, respectively. The partial derivatives of $Q_n$ may not exist for all $\theta \in \Theta$. For example, $Q_n$ and $\hat{Q}$ may not be defined outside $\Theta$, such as in a random coefficient model. In this case, $D Q_n$ and $D^2 Q_n$ would denote the first-order and second-order left or right partial derivatives. Also, $Q_n$ may not be smooth but only “stochastically differentiable" (see \cite{Pollard1985}), which can be the case in quantile estimation (e.g., \cite{PakesPollard1989}). Also following \cite{KETZ2018}, I allow for the case that these derivatives are not known exactly and must be approximated numerically using $\hat{D \hat{Q}}$ and $\hat{D^2 \hat{Q}}$. In subsection \ref{s: HighDimensional}, I show that the idea of using approximations of the derivatives can be extended for the orthogonalized version of the method, where part of the "approximation" is using the estimated rather than true values of additional high-dimensional parameters beyond $\theta$. The quasi-unconstrained estimator can then be defined as:
\begin{equation} \label{eq: Quasi-unconstrainedDefined}
\Tilde{\theta} = \hat{\theta} - (\hat{D^2 \hat{Q}}(\hat{\theta}))^{-1} \hat{D \hat{Q}} (\hat{\theta}).
\end{equation}
The compactness conditions in Assumption \hyperref[A1]{1} ensures that for any sequences $\{\gamma_n \}_{n \in \mathbb{N}}$ with $\gamma_n \in \Gamma$, there is a convergent subsequence. I denote these convergent sequences as $\Gamma(\gamma^*) = \{ \{\gamma_n\}_{n \in \mathbb{N}} \;|\; \gamma_n \rightarrow \gamma^* \}$ for some $\gamma^* := (\theta^*, \omega^*) \in \Gamma$. Using the results of \cite{ANDREWS2020}, the asymptotic distribution of the test statistics under these convergent subsequences is used to show that the tests control size uniformly. In order to show the asymptotic behavior of the test statistics under these convergent drifting sequences, I impose conditions in Assumption \hyperref[A2]{2} that guarantee that $\Tilde{\theta}$ is asymptotically normal under these convergent drifting sequences.
\textbf{Assumption 2} \phantomsection\label{A2} Under $\{ (\theta_n, \omega_n) \}_{n \in \mathbb{N}} = \{ \gamma_n \}_{n \in \mathbb{N}} \in \Gamma(\gamma^*)$ for any $\theta^* \in \Theta$,
\begin{enumerate}
\item $\sqrt{n}(\hat{\theta} - \theta_n) = O_p(1)$.
\item $\sqrt{n} D \hat{Q}(\theta_n) \overset{d}{\rightarrow} N(0, V(\gamma^*))$, where $V(\gamma^*)$ is symmetric and positive definite.
\item $D^2 \hat{Q} (\theta_n) \overset{p}{\rightarrow} \mathcal{J}(\gamma^*)$, where $\mathcal{J}(\gamma^*)$ is symmetric and non-singular.
\item $|| \hat{D \hat{Q}}(\hat{\theta}) - D \hat{Q}(\theta_n) - D^2 \hat{Q} (\theta_n)(\hat{\theta} - \theta_n)||_2 = o_p(1/\sqrt{n})$ and $|| \hat{D^2 \hat{Q}}(\hat{\theta}) - D^2 \hat{Q} (\theta_n)||_2 = o_p(1)$.
\item $\hat{V}(\hat{\theta}) \overset{p}{\rightarrow} V(\gamma^*)$.
\end{enumerate}
Under Assumption \hyperref[A2]{2}, by Theorem 1 of \cite{KETZ2018}, under $\{ \gamma_n \}_{n \in \mathbb{N}} \in \Gamma(\gamma^*)$, $$\sqrt{n}(\Tilde{\theta} - \theta_n) \overset{d}{\rightarrow} \mathcal{Z} \sim N(0_J, \Sigma(\gamma^*))$$ where $\Sigma(\gamma^*) = \mathcal{J}(\gamma^*)^{-1}V(\gamma^*)\mathcal{J}(\gamma^*)^{-1}$ and $\hat{\Sigma} = (\hat{D^2 \hat{Q}}(\hat{\theta)})^{-1} \hat{V}(\hat{\theta})(\hat{D^2 \hat{Q}}(\hat{\theta}))^{-1} \overset{p}{\rightarrow} \Sigma(\gamma^*)$. In other words, $\Tilde{\theta}$ is asymptotically normal and its asymptotic variance can be consistently estimated.\footnote{For notational simplicity, I suppress the dependence of $\Sigma$ on $\gamma^*$ below.}
\par
I show in Appendix B in the Supplemental Materials how Assumption \hyperref[A2]{2.2} is satisfied when $\sqrt{n}D\hat{Q}$ converges weakly to a Gaussian process on $\Theta$. Since $D \hat{Q} (\theta)$ often involves a sample average, in many cases this can be shown by applying the proper Uniform Central Limit Theorem. The asymptotic normality condition in Assumption \hyperref[A2]{2.2} provides the basis for the asymptotic normality of $\Tilde{\theta}$, because the quasi-unconstrained estimator will asymptotically satisfy a first-order condition. In order for Assumption \hyperref[A2]{2.3} to hold, we want that $\sup_{\theta \in \Theta} |D^2 \hat{Q} (\theta) - D^2 Q_n(\theta)| \overset{p}{\rightarrow} 0$. When $D \hat{Q}$ and $D^2 \hat{Q}$ can be found analytically, Assumption \hyperref[A2]{2.4} is implied by Assumption \hyperref[A2]{2.1} and $\hat{Q}$ being twice continuously fully differentiable. If $D \hat{Q}$ and $D^2 \hat{Q}$ must be approximated numerically, Assumption \hyperref[A2]{2.4} requires that these approximations converge to the actual first and second generalized derivatives at a sufficiently fast rate as the sample size grows. Assumption \hyperref[A2]{2.5} simply requires that the variance of the derivative at $\theta_n$ can be consistently estimated using $\hat{\theta}$. Assumption \hyperref[A2]{2.1} imposes that $\hat{\theta}$ is $\sqrt{n}$-consistent without requiring that it is asymptotically normal. Given that the other conditions in Assumption \hyperref[A2]{2} hold, only one weaker additional condition is needed to show that any consistent estimator is also $\sqrt{n}$-consistent.
\textbf{Lemma 1 ($\sqrt{n}$-Consistency)}\phantomsection\label{L1} Suppose that Assumptions \hyperref[A1]{1} and \hyperref[A2]{2.2}-\hyperref[A2]{2.4} hold. Furthermore, suppose that $\hat{\theta} - \theta_n = o_p(1)$ and $\hat{\theta}$ satisfies equation \eqref{eq: ConstrainedDefined}, $\mathcal{J}(\gamma^*)$ is positive-definite, and for any $\epsilon > 0$, there exists $\kappa_\epsilon$, such that $$P(\sup_{\theta : ||\theta - \theta_n||_2 \leq \kappa_\epsilon} ||D^2 \hat{Q}(\theta) - D^2 \hat{Q}(\theta_n)||_2 > \epsilon) \rightarrow 0.$$
Then $\sqrt{n}(\hat{\theta} - \theta_n) = O_p(1)$.
\par
This additional convergence condition will hold when $D^2 \hat{Q}$ converges uniformly over some sequence of neighborhoods of $\theta_n$, and it does not impose a rate of convergence for $D^2 \hat{Q}$. Therefore, under fairly weak conditions, even when $\hat{\theta}$ is not asymptotically normal due to $\theta_n$ being on or close to the boundary, it is still $\sqrt{n}$-consistent.
\section{Full Vector Inference} \label{s: Inference}
I now show that tests using this quasi-unconstrained estimator can control the size uniformly over $\Theta$ asymptotically. Let $CP_n(\gamma)$ denote the coverage probability for $\theta$ of a given test when the sample size is $n$ (i.e., the probability of failing to reject the null hypothesis that $\theta_n = \theta$). Following \cite{ANDREWS2020}, let $$AsySz = \liminf_{n \rightarrow \infty} \inf_{\gamma \in \Gamma} CP_n(\gamma)$$ be the asymptotic size and let $$AsyMaxCP = \limsup_{n \rightarrow \infty} \sup_{\gamma \in \Gamma} CP_n(\gamma)$$ be the asymptotic maximum coverage probability. When the nominal size of a test is equal to $\alpha$, a test provides asymptotic size control in a uniform sense when $AsySz \geq 1 - \alpha$. In this case, we can be confident that the test does not over-reject when the sample size is large, even when the parameter is on or near the boundary. Uniform size control is shown in the proofs using the general technique of \cite{ANDREWS2020} of showing that for any sequence $\{\gamma_n\}_{n \in \mathbb{N}}$, there is a subsequence under which size is being controlled. \cite{KETZ2018} also uses the technique of \cite{ANDREWS2020}, but focuses specifically on when $\Theta$ is a product space. When $\Theta$ is a product space, there are no interactions between when one element is close to the boundary and another element is. As a result, the problem of considering all sequences where $\theta_n$ converges to a point on the boundary can be reduced to considering only a finite number of cases based on whether each element of $\theta$ is converging to one of its bounds. The key additional insight used in the proof is that because $\Theta \subset \mathbb{R}^J$, the sequence of sets $\sqrt{n}(\Theta - \theta_n)$ must have a convergent subsequence, and I am able to focus on analyzing the asymptotic distribution of the test statistic under these subsequences. I first consider the Wald test statistic, which is simply equal to
\begin{equation*}
(\Tilde{\theta} - \theta)'(\hat{\Sigma}/n)^{-1}(\Tilde{\theta} - \theta).
\end{equation*}
I show that it not only controls size, but also satisfies a stronger asymptotic similarity condition. Namely, that $AsySz = AsyMaxCP = 1 - \alpha$.
\textbf{Proposition 1}\phantomsection\label{P1} (\textbf{Wald Statistic}) Let
\begin{equation*}
CP_n(\bar\gamma) = P_{\bar\gamma}((\Tilde{\theta} - \bar\theta)'(\hat{\Sigma}/n)^{-1}(\Tilde{\theta} - \bar\theta) \leq cv_{1-\alpha}(\chi^2_J)),
\end{equation*}
where $cv_{1-\alpha}(\chi^2_J)$ denotes the $1-\alpha$ critical value of a Chi-Squared distribution with $J$ degrees of freedom. If Assumptions \hyperref[A1]{1} and \hyperref[A2]{2} hold, then $AsySz = AsyMaxCP = 1 - \alpha$.
\par
The asymptotic similarity condition in Proposition \hyperref[P1]{1} shows that the Wald test is also not unnecessarily conservative asymptotically. However, I also consider a Likelihood Ratio (LR) test because it may have superior power properties. The LR test statistic is given by,
\begin{equation*}
LR(\theta,\Tilde{\theta}) = (\Tilde{\theta} - \theta)'(\hat{\Sigma}/n)^{-1}(\Tilde{\theta} - \theta) - \inf_{b \in \Theta} (\Tilde{\theta} - b)'(\hat{\Sigma}/n)^{-1}(\Tilde{\theta} - b),
\end{equation*}
\par
where here the likelihood is based on $\Tilde{\theta}$ being asymptotically normal. Note that if the estimator was unconstrained, so $\Theta = \mathbb{R}^J$, then this would be equivalent to the Wald statistic. Relatedly, we could also consider the Lagrange Multiplier test statistic, which is typically defined as $\nabla_{\theta,n}'\hat{\Sigma}\nabla_{\theta,n}$ where $\nabla_{\theta,n}$ is the gradient of the log-likelihood function at the maximum value under the null hypothesis. However, if we again base the likelihood on the asymptotic normality of $\Tilde{\theta}$, then this will again give a test statistic proportional to the Wald test statistic. This is unsurprising given that these three test statistics are known to be identical when the log-likelihood is simply a quadratic function of $\theta$ and $\theta$ is unconstrained (see, e.g., Lemma 1 of \cite{ENGLE1984}). I therefore don't provide separate formal results for the LM test statistic in this section, but do separately consider a Conditional Lagrange Multiplier (CLM) statistic for subvector inference in section \ref{s: InferenceForSubvectors}. One difference between my LR test statistic and the Conditional Likelihood Ratio Test (CLR) used by \cite{KETZ2018}, is the feasible set. Since they assume $\Theta$ is the Cartesian product of intervals of the form $\theta_j \in [0,c]$ or $\theta_j \in [-c,c]$ for some large constant $c$, they set the feasible set in the CLR statistic to be equal to the Cartesian product of the extensions of these intervals (i.e., $\theta_j \in [0,\infty)$ or $\theta_j \in (-\infty, \infty)$). They do this because it makes sure that only boundaries at zero are taken into account, and they focus on drifting sequences where elements of $\theta_n$ may be drifting towards zero. In my approach, I use $\Theta$ as the feasible set to be agnostic about which part of the boundary the parameter is close to. $\Theta$ is used not only in calculating the test statistic, but also in calculating the critical values. For the Wald test, calculating p-values and confidence intervals is straightforward because the critical values are always calculated using a chi-squared distribution. For the LR test, the critical values come from the distribution of $LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta)$ where $\mathcal{Z} \sim N(0_J, \hat{\Sigma})$. Specifically, critical values are given by $$cv_{1-\alpha}(LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta)) = \inf \{q \in \mathbb{R} \;|\: P( LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta) \leq q ) \geq 1- \alpha).$$ Therefore, the critical values depend on $\theta$ and $\Theta$ and in general do need to be calculated numerically.
\textbf{Proposition 2 (Likelihood Ratio)} \phantomsection\label{P2} Let
\begin{equation*}
CP_n(\bar\gamma) = P_{\bar\gamma}(LR(\bar\theta, \Tilde{\theta}) \leq cv_{1-\alpha}(LR(\bar\theta, n^{-1/2}\hat{\mathcal{Z}} + \bar\theta))).
\end{equation*}
If Assumptions \hyperref[A1]{1} and \hyperref[A2]{2} hold, then $AsySz \geq 1 - \alpha$.
\par
While both the Wald and LR statistics provide asymptotic size control, the stronger similarity condition of $AsySz = AsyMaxCP = 1- \alpha$ does not hold for the LR test in general. This is related to the fact that the Wald test is pivotal but the LR test is not. For the Wald test, its critical values always come from a chi-squared distribution so they always come from a continuous distribution, but the distribution of $LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta)$ does have discontinuities in some cases. For example, in some cases there is a positive probability of the two terms that compose the LR statistic being equal to each other, so the distribution function of $LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta)$ does have a discontinuity at zero.\footnote{This is true in a case discussed by \cite{KETZ2018} where $\Theta = [0, \infty)^J$ and $\theta$ is equal to the zero vector.} Generally, the distribution of $LR(\theta, n^{-1/2}\hat{\mathcal{Z}} + \theta)$ must be calculated numerically by repeatedly sampling from $\hat{\mathcal{Z}}$. When $\theta$ is low-dimensional, this is generally not computationally intensive, although the critical values depend on the null value, so this must be done repeatedly to invert the test for confidence sets. I show in Appendix C of the Supplemental Materials that, depending on the form of the parameter space $\Theta$, we can increase the computational speed of projecting $n^{-1/2}\hat{\mathcal{Z}} + \theta$ using a change of variable and approximations for the new feasible set after that change of variable. When the full vector is high-dimensional, but the parameter of interest is low-dimensional, one extra benefit of using the orthogonalized version of this method will be helping to retain computational tractability for conducting inference as the resampling and projecting will only need to be done in a low-dimensional space.
\section{Subvector Inference} \label{s: InferenceForSubvectors}
I next consider the case in which the null hypothesis specifies only the value of a subvector of $\theta$. Following \cite{KETZ2018}, let $\theta = (\beta , \delta)$ where $\beta \in \mathbb{R}^K$ is the parameter of interest and $\delta \in \mathbb{R}^L$ is a nuisance parameter. Let $\Bar{\theta} = (\Bar{\beta}, \Bar{\delta})$ and let
\begin{equation*}
\begin{bmatrix}
\Sigma_{\beta \beta} & \Sigma_{\beta \delta} \\
\Sigma_{\delta \beta} & \Sigma_{\delta \delta}
\end{bmatrix}
\end{equation*}
denote the conformable submatrices of $\Sigma$. We are interested in testing
\begin{equation*}
H_0 \;:\: \Bar{\beta} = \beta_0, \; \Bar{\delta} \in D_{\beta_0} \quad vs.\quad H_1 \;:\: \Bar{\beta} \ne \beta_0,\; \Bar{\beta} \in B = \{\beta \in \mathbb{R}^K | (\beta,\delta) \in \Theta \text{ for some } \delta\}
\end{equation*}
where $D_{\beta_0} := \{\delta \in \mathbb{R}^L \;|\; (\beta_0, \delta) \in \Theta \}$. For the special case considered by \cite{KETZ2018} where $\Theta$ is equal to the Cartesian Product of closed, possibly infinite, intervals, there is no interaction between the true value of the subvector of interest and the feasible set for the nuisance parameter. However, in general $D_{\beta_0}$ depends on $\beta_0$.
\par
The results from the Wald statistics can be directly extended to this case, where the test statistic is now equal to
\begin{equation*}
(\Tilde{\beta} - \beta_0)'(\hat{\Sigma}_{\beta \beta}/n)^{-1}(\Tilde{\beta} - \beta_0).
\end{equation*}
This is because the point-wise asymptotic distribution of the Wald statistic for $\Bar{\beta}$ does not depend on $\Bar{\delta}$ and is therefore still fully specified under the null hypothesis. However, for the LR and Lagrange Multiplier test, the true value $\Bar{\delta}$ does affect the asymptotic behavior of the test statistic. Luckily, this can be addressed using the conditioning principle of \cite{Moreira2003} by performing CLR and CLM tests. Under the null hypothesis, one can use $X_n := \Tilde{\delta} - \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1} \Tilde{\beta}$ as a sufficient statistic for the unknown nuisance parameter $\Bar{\delta}$ that is asymptotically independent of $\Tilde{\beta}$, so the distribution of the CLR statistic conditional on $X_n$ is asymptotically nuisance-parameter-free. Using this sufficient statistic, the CLR statistic is given by
\begin{equation*}
CLR(\beta_0, \Tilde{\beta},\hat{\Sigma}/n, X_n) =
\inf_{d \in D_{\beta_0}}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
\Tilde{\delta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
\Tilde{\delta} - d
\end{pmatrix} -
\inf_{(b, d) \in \Theta}
\begin{pmatrix}
\Tilde{\beta} - b \\
\Tilde{\delta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - b \\
\Tilde{\delta} - d
\end{pmatrix}
\end{equation*}
\begin{equation*}
= \inf_{d \in D_{\beta_0}}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1} \Tilde{\beta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1} \Tilde{\beta} - d
\end{pmatrix}
\end{equation*}
\begin{equation*}
-
\inf_{(b, d) \in \Theta}
\begin{pmatrix}
\Tilde{\beta} - b \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1} \Tilde{\beta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - b \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1} \Tilde{\beta} - d
\end{pmatrix}.
\end{equation*}
The critical values for the CLR test are based on a normal distribution projected onto a set which now depends on $X_n$. More precisely, the critical values come from
\begin{equation*}
CLR(\beta_0, n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0, \hat{\Sigma}/n, x) =
\end{equation*}
\begin{equation*}
\inf_{d \in D_{\beta_0}}
\begin{pmatrix}
n^{-1/2}\hat{\mathcal{Z}}_b \\
x + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}(n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0) - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
n^{-1/2}\hat{\mathcal{Z}}_b \\
x + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}(n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0) - d
\end{pmatrix}
\end{equation*}
\begin{equation*}
-
\inf_{(b,d) \in \Theta}
\begin{pmatrix}
n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0 - b \\
x + (\hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}(n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0)) - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0 - b \\
x +(\hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}(n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0)) - d
\end{pmatrix},
\end{equation*}
where $x$ is the realized value of $X_n$ and $\hat{\mathcal{Z}}_b$ with $\hat{\mathcal{Z}}_b \sim N (0_K, \hat{\Sigma}_{\beta \beta})| X_n = x$.
So the critical value at the $1- \alpha$ level is given by
\begin{equation*}
cv_{1-\alpha}(CLR(\beta_0, n^{-1/2}\hat{\mathcal{Z}}_b + \beta_0,\Sigma, x)|x = X_n,\Sigma = \hat{\Sigma}/n)
\end{equation*}
\begin{equation*}
= \inf \{q \in \mathbb{R}\;|\: P(CLR(\beta, n^{-1/2}\hat{\mathcal{Z}}+\beta_0,\Sigma, x) \leq q \;|\;x = X_n, \Sigma = \hat{\Sigma}/n) \geq 1 - \alpha\}.
\end{equation*}
\textbf{Proposition 3 (Conditional Likelihood Ratio)} \phantomsection\label{P3} Let
\begin{equation*}
CP_n(\bar\gamma) = P_{\bar\gamma}(CLR(\bar\beta, \Tilde{\beta},\hat{\Sigma}/n,X_n) \leq cv_{1-\alpha}(CLR(\bar\beta, n^{-1/2}\hat{\mathcal{Z}}_b +\beta,\hat{\Sigma}/n,X_n))).
\end{equation*}
If Assumptions \hyperref[A1]{1} and \hyperref[A2]{2} hold, then $AsySz \geq 1 - \alpha$.
As mentioned in section \ref{s: Inference}, the Lagrange Multiplier test statistic is typically equal to $\nabla_{\theta, n}' \hat{\Sigma} \nabla_{\theta, n}$. However, because $\Tilde{\theta}$ is asymptotically normal, this is proportional to
\begin{equation*}
LM(\beta_0, \Tilde{\beta}) = \inf_{d \in D_{\beta_0}} \begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
\Tilde{\delta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
\Tilde{\delta} - d
\end{pmatrix}.
\end{equation*}
Note that this is equal to the first term of the CLR statistic, so the same approach can be used for the CLM statistic. Hence, we can express the test statistic as
\begin{equation*}
CLM(\beta_0, \Tilde{\beta}, \hat{\Sigma}/n, X_n) = \inf_{d \in D_{\beta_0}}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}\Tilde{\beta} - d
\end{pmatrix}'
(\hat{\Sigma}/n)^{-1}
\begin{pmatrix}
\Tilde{\beta} - \beta_0 \\
X_n + \hat{\Sigma}_{\delta \beta} \hat{\Sigma}_{\beta \beta}^{-1}\Tilde{\beta} - d
\end{pmatrix}.
\end{equation*}
Then, by the same method as for the CLR, the CLM can be shown to control size uniformly over $\Theta$.
\textbf{Proposition 4 (Conditional Lagrange Multiplier)} \phantomsection\label{P4} Let
\begin{equation*}
CP_n(\bar\gamma) = P_{\bar\gamma}(CLM(\bar\beta, \Tilde{\beta},\hat{\Sigma}/n,X_n) \leq cv_{1-\alpha}(CLM(\bar\beta, n^{-1/2}\hat{\mathcal{Z}}_b +\bar\beta,\hat{\Sigma}/n,X_n))).
\end{equation*}
If Assumptions \hyperref[A1]{1} and \hyperref[A2]{2} hold, then $AsySz \geq 1 - \alpha$.
Both the CLR and CLM fail to achieve the uniform asymptotic similarity condition for the same reason as the unconditional LR test. Because the distribution used to calculate critical values involves projecting a normal distribution onto a set, there may be discontinuities in that distribution. As a result, it is possible that the threshold cannot be chosen so that the probability of exceeding it is exactly equal to $\alpha$ under the null.
\subsection{Low Dimensional Subvector of a High Dimensional Vector}\label{s: HighDimensional}
All results presented so far, as well as those in \cite{KETZ2018}, assume that the dimension of the estimated parameters is fixed, so they may only provide a good approximation in finite samples where the number of estimated parameters is small relative to $n$. However, it is not uncommon to have applications in which the parameter vector is both constrained and high-dimensional. For example, the control weights for a synthetic control unit are usually constrained and often have a dimension comparable to the number of time periods. Also, some applications of the random coefficient model of \cite{BLP} have many markets or many products so $\theta$ is both high-dimensional and has some of its elements constrained (see, for example, \cite{Armstrong2016}). I explore another example of where this arises in the empirical application in section \ref{s: Network application}.
\par
The strategy in this section is to first estimate the high-dimensional nuisance parameter, possibly using regularization, and then estimate the remaining
parameters using a Neyman orthogonalized score. The reason this is useful is that the orthogonalized score, with the estimated nuisance parameters plugged in, can be treated as the criterion function to which the results of Sections \ref{s: MandE} through \ref{s: InferenceForSubvectors} apply. To make that correspondence transparent, it is convenient to slightly change the notation. From this point forward, $\theta \in \Theta \subseteq \mathbb{R}^{J}$ denotes only the parameters we do not orthogonalize with respect to, and $\psi \in \Psi_n \subseteq \mathbb{R}^{P}$ denotes the high-dimensional nuisance parameter we do, with $J$ fixed and $P$ permitted to grow with $n$. The criterion function before orthogonalization, therefore, depends on the pair $(\theta,\psi)$, and it is this pair that corresponds to what was called $\theta$ previously.\footnote{Note that $\psi$ and the nuisance subvector $\delta$ of Section \ref{s: InferenceForSubvectors} are distinct objects handled by different devices, orthogonalization in the former case and conditioning in the latter.} After orthogonalizing, the resulting moment conditions are a function of $\theta$ alone.
\par
The population orthogonalized score can be defined as
\begin{equation}\label{Neyman Orthogonal Score}
M(\theta, \psi, \eta) = D_\theta Q_n(\theta,\psi) - \eta D_\psi Q_n(\theta,\psi),
\end{equation}
where $\eta \in \mathbb{R}^{J \times L}$ is an additional nuisance parameter and the true value of this matrix $\eta_n$ satisfies $$D_{\theta \psi} Q_n(\theta_n) - \eta_n D_{\psi \psi}(\theta_n) = 0_{J \times L}.$$
Using some estimate of this matrix $\hat{\eta}$, we can define the sample orthogonalized score as
\begin{equation}\label{Sample Neyman Orthogonal Score}
\hat{M}(\theta,\hat{\psi},\hat{\eta}) = \hat{D_\theta \hat{Q}}(\theta,\hat{\psi}) - \hat{\eta} \hat{D_\psi \hat{Q}}(\theta,\hat{\psi}).
\end{equation}
Under suitable rates of convergence conditions on $\hat{\eta}$ and $\hat{\psi}$ and smoothness conditions on $Q_n$,
\begin{equation}\label{eq: adaptivity}
\sqrt{n}(\hat{M}(\theta_n,\hat{\psi},\hat{\eta}) - \hat{M}(\theta_n,\psi_n,\eta_n)) \overset{p}{\rightarrow} 0,
\end{equation}
and as a result, an estimator using the moment conditions $\hat{M}(\theta,\hat{\psi},\hat{\eta})$ can be asymptotically equivalent to one using $\hat{M}(\theta,\psi_n,\eta_n)$. There is a direct parallel between the adaptivity condition and the conditions imposed on the numerically approximated derivatives in Assumption \hyperref[A2]{2.4}. There, $\hat{D\hat{Q}}$ and $\hat{D^{2}\hat{Q}}$ were approximations to $D\hat{Q}$ and $D^{2}\hat{Q}$, and what was required was that the approximation error vanish quickly enough that $\sqrt{n}\hat{D\hat{Q}}$ and $\hat{D^{2}\hat{Q}}$ behave asymptotically as $\sqrt{n}D\hat{Q}$ and $D^{2}\hat{Q}$ do. In the orthogonalized version of the method, plugging the estimates $\hat{\psi}$ and $\hat{\eta}$ into the orthogonalized moments plays exactly this role. Because $\hat{\theta}$ is $\sqrt{n}$-consistent, the asymptotic distribution of $\tilde{\theta}$ will depend on the limiting behavior of $\sqrt{n}\hat{M}(\theta_n,\hat{\psi},\hat{\eta})$ and
$\partial_{\theta}\hat{M}(\theta_n,\hat{\psi},\hat{\eta})$, and what we require is that these be asymptotically equivalent to $\sqrt{n}\hat{M}(\theta_n,\psi_n,\eta_n)$ and
$\partial_{\theta}\hat{M}(\theta_n,\psi_n,\eta_n)$. The rates required are the same as in Assumption \hyperref[A2]{2.4} since the level must agree to order $o_p(1/\sqrt{n})$, which is exactly the adaptivity condition in equation \eqref{eq: adaptivity}, while the derivative needs only to agree to order $o_p(1)$, which Assumptions \hyperref[A3]{3.2} and \hyperref[A3]{3.5} impose. In this sense the estimation error in $\hat{\psi}$ and $\hat{\eta}$ can be treated as one more source of approximation error in the derivatives of the criterion function, and the reasoning of Section \ref{s: MandE} applies with $\hat{M}(\theta,\hat{\psi},\hat{\eta})$ in place of $\hat{D\hat{Q}}(\theta)$. This also means the two sources of approximation can be combined. If the derivatives of the orthogonalized criterion must themselves be computed numerically, it is enough that the total error from both numerical approximation and nuisance parameter estimation satisfies these rates.
\par
Rather than taking one Newton-Raphson step using the actual criterion function, we can instead use a GMM objective function:
\begin{equation}\phantomsection\label{eq: GMM Estimator}
\hat{M}(\theta,\hat{\psi},\hat{\eta})' \hat{M}(\theta,\hat{\psi},\hat{\eta}).
\end{equation}
Therefore, the quasi-unconstrained estimator is now equal to
\begin{equation}
\Tilde{\theta} = \hat{\theta} - (\partial_\theta \hat{M}(\hat{\theta},\hat{\psi},\hat{\eta})' \partial_\theta \hat{M}(\hat{\theta},\hat{\psi},\hat{\eta}))^{-1}\partial_\theta \hat{M}(\hat{\theta},\hat{\psi},\hat{\eta})' \hat{M}(\hat{\theta},\hat{\psi},\hat{\eta}).
\end{equation}
\par
This estimator can be considered as minimizing the quadratic approximation of the GMM objective function at the initial estimate $\hat{\theta}$. Note that in the special case where $\hat{M}(\theta,\hat{\psi},\hat{\eta})$ is linear in $\theta$, this is equivalent to setting the orthogonalized moment conditions equal to zero, so $\Tilde{\theta}$ be defined by
$$ \hat{M}(\Tilde{\theta},\hat{\psi},\hat{\eta}) = 0.$$
\textbf{Assumption 3}\phantomsection\label{A3} Under $\{ (\theta_n, \omega_n) \}_{n \in \mathbb{N}} = \{ \gamma_n \}_{n \in \mathbb{N}} \in \Gamma(\gamma^*)$ for any $\theta^* \in \Theta$,
\begin{enumerate}
\item $\sqrt{n}(\hat{\theta} - \theta_n) = O_p(1)$.
\item $\hat{Q}$ is twice continuously differentiable in $\theta$ and for or all $\epsilon >0$, there exists $\kappa_\epsilon$ such that $$P(\sup_{\theta \in \Theta_n, \psi \in \Psi_n : ||\theta - \theta_{n}||_1,||\psi - \psi_n||_1 < \kappa_\epsilon}||\hat{D_{\theta \theta} \hat{Q}}(\theta,\psi) - \hat{D_{ \theta \theta} \hat{Q}}(\theta_{n},\psi_n)||_2 > \epsilon ) \rightarrow 0$$
$$\text{ and } P(\sup_{\theta \in \Theta_n,\psi \in \Psi_n : ||\theta - \theta_{n}||_1,||\psi - \psi_n||_1 < \kappa_\epsilon}||\hat{D_{\theta \psi} \hat{Q}}(\theta,\psi) - \hat{D_{\theta \psi} \hat{Q}}(\theta_{n},\psi_{n})||_2 > \epsilon ) \rightarrow 0.$$
\item $\sqrt{n}\hat{M}(\theta_{n},\psi_{n},\eta_{n}) \overset{d}{\rightarrow} N(0,V_M)$ for some positive definite matrix $V_M$.
\item $\partial_\theta \hat{M}(\theta_{n},\psi_{n},\eta_{n}) \overset{p}{\rightarrow} M_\theta$ for some matrix $M_\theta$ with $rank(M_\theta) = J$.\footnote{Here, the dependence of $M_\theta$ and $V_M$ on $\gamma^*$ have been suppressed for notational convenience.}
\item $||\hat{\psi} - \psi_n||_1 = o_p(1)$ and $||vec(\hat{\eta}) - vec(\eta_n)||_1 = o_p(1)$.
\end{enumerate}
In Appendix C, I give examples of sufficient conditions for the adaptivity condition in equation \eqref{eq: adaptivity}.\footnote{A number of existing papers provide conditions under which this will hold, including \cite{fry2024}, \cite{belloni2018highdimensional}, and others.} In order to achieve this condition as well as the consistency condition in Assumption \hyperref[A3]{3.5}, it is common for $\psi_{n}$ and $\eta_{n}$ to be sparse and for them to be estimated by either penalizing or constraining their L1 norm. The empirical application in section \ref{s: Network application}, where $\psi$ corresponds to the weights in a sparse unobserved network, provides an example where this is the case.
\textbf{Proposition 5}\phantomsection\label{P5} Suppose under $\{\gamma_n\}_{n \in \mathbb{N}} \in \Gamma(\theta^*)$ for any $\theta^* \in \Theta$, $\theta_{n}$, $\psi_n$, and $\eta_{n}$ satisfies orthogonality equation for each $n$, the adaptivity condition in equation \eqref{eq: adaptivity} holds, and Assumptions \hyperref[A1]{1} and \hyperref[A3]{3} hold. Then as $n \rightarrow \infty$ with $J$ fixed and $L \rightarrow \infty$,
$$\sqrt{n}(\Tilde{\theta} - \theta_{n}) \overset{d}{\rightarrow} N(0,\Sigma(\gamma^*)),$$
where $\Sigma(\gamma^*) = (M_\theta' M_\theta)^{-1} M_\theta' V_M M_\theta (M_\theta' M_\theta)^{-1}$.
\par
The earlier results then apply with $\theta$ read as the full parameter vector, with $\psi_n$ and $\eta_n$ absorbed into the infinite-dimensional component $\omega_n$ of $\gamma_n$, so that Assumption \hyperref[A1]{1} and the convergent subsequences $\Gamma(\gamma^{*})$ retain their meaning without modification. Which results apply depends on what $\theta$ contains. If every element of $\theta$ is of interest, the full-vector methods of Section \ref{s: Inference} can be used directly. If $\theta$ itself splits into a parameter of interest and a remaining low-dimensional nuisance parameter, the subvector methods of Section \ref{s: InferenceForSubvectors} apply to that split, with the conditioning argument now operating within $\theta$ rather than on the original full vector.
The proof follows similar reasoning as Theorem 1 of \cite{KETZ2018}. As a result, size can be controlled in a uniform sense using the methods in section \ref{s: Inference} and \ref{s: InferenceForSubvectors} provided we also have a consistent estimator of the asymptotic variance. However, I have implicitly assumed in this section that the parameter space for $\theta$ does not depend on $\psi$, which is an important caveat for using the LR and LM tests in this setting.
\textbf{Remark 1 (Interactive Constraints Between High and Low Dimension Parameters)} One complication of the orthogonalization approach that the conditional approach of Section \ref{s: InferenceForSubvectors} does not share is accommodating constraints involving both $\psi$ and $\theta$. In the CLR and CLM statistics, the first term is minimized over all $\delta \in D_{\beta_0}$ and the second over all $\theta \in \Theta$, so the statistic is computed without imposing any restriction on how the constraints on the parameter of interest and those on the nuisance parameter interact. Here, by not having $\Theta$ depend on $\psi$, we are not allowing interactive constraints as we were when allowing $D$ to depend on $\beta$. Because we do not want to rely on $\sqrt{n}$-consistency and asymptotic normality of the estimate of the $P$-dimensional vector $\psi$, and because minimizing over a $P$-dimensional set is computationally burdensome, we now form the likelihood using only $\tilde{\theta}$ and $\hat{\Sigma}$. The feasible set entering the LR and LM statistics is then a set for $\theta$ alone. Writing $\mathcal{S}_n \subseteq \mathbb{R}^{J} \times \mathbb{R}^{P}$ for the joint constraint set on $(\theta,\psi)$, the relevant set is
\[
\Theta_{\psi_n} := \{\theta \in \mathbb{R}^{J} : (\theta,\psi_n) \in
\mathcal{S}_n\},
\]
which depends on the unknown $\psi_n$. For the statistic and its critical values to be computable, the feasible set for $\theta$ must therefore not vary with the value of the nuisance parameter, which is to say we require a product structure between the two blocks:
\[
\mathcal{S}_n = \Theta \times \Psi_n, \qquad \Theta \subseteq \mathbb{R}^{J},
\quad \Psi_n \subseteq \mathbb{R}^{P}.
\]
This is weaker than the product-space assumption in \cite{KETZ2018}. Interactions among the elements of $\theta$ are still permitted, as are interactions among the elements of $\psi$, including the simplex
constraints on the network weights in Section \ref{s: Network application}. What is ruled out is a constraint that links the two blocks. Lastly, this requirement applies only to the LR- and LM-based tests. The Wald statistic uses the parameter space neither in the statistic nor in its critical values, so the subvector Wald test remains available with no condition of this kind.\footnote{The proof of Proposition \hyperref[P5]{5} does not rely on this product space assumption, so the only change that would be required in this case is to modify Assumption \hyperref[A3]{3.2}'s convergence conditions to hold with the supremum taken over $(\theta,\psi) \in \mathcal{S}_n: ||\theta -\theta_n||,||\psi-\psi_n||_1 < \kappa_\epsilon$.}
\section{Asymptotic Power}\label{s: Power}
\cite{KETZ2018} analyzes the point-wise-in-$\beta$ power properties of their version of the CLR test when $\Theta$ is a product space and the subvector being tested $\beta$ is a scalar. They show that when the CLR test is point-wise similar, it is also admissible and extended Weighted Average Power (WAP) maximizing in the class of all tests with correct asymptotic size.\footnote{Extended-WAP-similar is defined \cite{montiel_olea_2020} to be a similar test that, for any tolerance $\epsilon > 0$, there exist some weighting of alternatives such that the test is within $\epsilon$ of being WAP-maximizing in the set of similar tests.} The core of their argument relies on showing that the test is similar with convex acceptance regions and using an argument based on the complete class result of \cite{MatthesTruax1967} to show that this means the CLR is admissible. This can then be connected to being WAP-maximizing using the results of \cite{montiel_olea_2020}, who shows that a test is similar and admissible if and only if it is similar and essentially WAP-maximizing (i.e., the limit of a sequence of WAP-maximizing tests). Extending \cite{KETZ2018}'s power results beyond a product parameter space is more delicate than extending the size results, largely because they are able to treat testing for a boundary point as a single special case (because the only boundary case is $H_0: \beta = 0$).\footnote{More specifically, for the boundary case, they must add a condition that $\alpha \leq 1/2$ and there is positive covariance of $\beta$ with the elements of $\delta$ which are also constrained to be non-negative.}
\par
To see why it is necessary to handle boundary points separately, first note that when testing points on the boundary, admissibility is not guaranteed by convex acceptance sections and similarity alone. This can be straightforwardly shown using the Wald test. Consider the class of tests $\mathcal{C}$ that asymptotically have size less than or equal to $\alpha$ under a data generating processes that induce $\sqrt{n}\hat{\Sigma}^{-1/2}(\Tilde{\theta} - \theta_0) \overset{d}{\rightarrow} N(0_J,I_J)$. Based on \cite{Muller2011}, for any test that is a function of this quasi-unconstrained estimator $\phi(\sqrt{n}\hat{\Sigma}^{-1/2}(\Tilde{\theta} - \theta_0))$, its WAP converges to the WAP of $\phi$ in the limiting problem. That is, its asymptotic WAP can be found by considering $\phi(\mathcal{Z})$ where $\mathcal{Z} \sim N(0_J,I_J)$. As a result of this, it is possible to show whether a test of this form is asymptotically WAP-maximizing in the class $\mathcal{C}$ by showing whether it is WAP-maximizing in this class for the limiting problem. For the Wald test, it is easy to find counterexamples where it is not WAP-maximizing in the limiting problem when the null value is on the boundary. Consider the following example: Let $\Theta = [0,\infty)$, $H_0: \; \theta = 0$, and the true value of $\theta$ is some $\theta_0 \geq 0$. The power of the Wald test in the limiting problem is $1- P(|\Tilde{\theta}| \leq cv_{1-\alpha/2}(N(0,\Sigma)))$ where $\Tilde{\theta} \sim N(\theta_0,\Sigma)$. If instead a right-sided test was used, then the power would be $1- P(\Tilde{\theta} \leq cv_{1-\alpha}(N(0,\Sigma)))$. Then the difference in the right-sided and the two-sided power is given by
\begin{equation*}
P( cv_{1-\alpha}(N(0,\Sigma)) \leq \Tilde{\theta} \leq cv_{1-\alpha/2}(N(0,\Sigma))) - P(\Tilde{\theta} \leq cv_{\alpha/2}(N(0,\Sigma))),
\end{equation*}
which is positive for all $\theta_0 > 0$ and equal to zero when $\theta_0 = 0$. Therefore, the Wald test is inadmissible among the class $\mathcal{C}$. As the Wald test illustrates, a similar test with convex acceptance sections can still be dominated when the tested value lies on the boundary of the parameter of interest because a one-sided competitor exploits information about $\Theta$ when determining the rejection region that a two-sided statistic ignores. The CLM test can suffer from a similar limitation when $\theta$ is constrained. While the statistic does take into account $D_{\beta_0}$, it does not take into account how surprising $\Tilde{\theta}$ is under alternatives $\beta \in B$. Therefore, a analogous counterexample could be constructed where $\Theta = [0,\infty) \times \mathbb{R}$ and $H_0 : \beta = 0$. Similarly to the Wald statistic, the CLM test will reject this null when the magnitude of $\Tilde{\beta}$ is too large, but this is dominated by a right-sided test that only rejects when $\Tilde{\beta}$ takes on a large positive value. On the other hand, as \cite{KETZ2018} highlights, in this setting the CLR test for $H_0 : \beta = 0$ is equivalent to the right-sided test mentioned above as long as $\alpha \leq 1/2$, which is why it is admissible in this case.
\par
A second complication in boundary cases is that, due to discontinuities in the asymptotic distribution of the test statistic, the test may not be similar. In the case of \cite{KETZ2018}, they are able to handle the special case of testing $\beta = 0$ because the only discontinuity occurs at probability one-half, which does not present an issue as long as $\alpha \leq 1/2$. For a general parameter space, no such simple characterization of the discontinuities is available, so ensuring the test is similar becomes correspondingly more difficult, although it can still be handled in other special cases. For example, \cite{andrews1996} shows that when $\Theta$ is a positively homogeneous set (i.e., $\theta \in \Theta$ implies $\tau \theta \in \Theta$ for all $\tau > 0$), then the LR test is similar and admissible.\footnote{Relatedly, when $\Theta$ is a closed convex cone, the LR tests asymptotic distribution will be a mixture of $\chi^2$ distributions by Theorem 3.4.2 of \cite{ConstrainedInferenceTextbook}. This distribution is continuous on $(0,\infty)$ and has an atom at the origin of mass $\omega_0$, corresponding to realizations in which the two terms composing the statistic coincide. The $1-\alpha$ quantile is therefore a continuity point whenever $\omega_0 < 1- \alpha$, in which case the LR test is point-wise similar.} They show this for the finite-sample Gaussian model, but as noted above, \cite{Muller2011} results imply there is a direct connection between being similar and admissible in the limiting experiment and being asymptotically similar and admissible.\footnote{One route to recovering \cite{KETZ2018}'s result when $\Theta$ is not a product space is working instead in the limit experiment, where the localized parameter ranges over the tangent cone and considering the case where that tangent cone has a product structure, with its nuisance section equal to a linear subspace or a product of half-spaces and lines. This holds when a product space is used and when the constraints binding at the tested point restrict the parameter of interest while leaving the nuisance parameter unconstrained. However, I do not focus on this extension, as it fails to cover the main constraint geometries of interest, such as simplex constraints and positive semi-definite constraints.} I leave a general characterization of when the CLR and LR tests are similar to future work.
\par
Lastly, it is worth considering the power properties of tests using the orthogonalized QUE. Subsection \ref{s: HighDimensional} focuses on settings with high-dimensional nuisance parameters, where admissibility results are going to be harder to establish. While the technique could also be used to handle low-dimensional nuisance parameters, it does not offer the same benefit when the full vector is $\sqrt{n}$-consistent and asymptotically normal, since we can already construct asymptotically pivotal test statistics in this case. When using the orthogonalized estimator in the low-dimensional setting, the Wald test will be inadmissible in some cases by the same reasoning discussed above (for example, just consider a setting with $\beta$ non-negativity constrained and the criterion function being defined using a set of moment conditions that are already orthogonal with respect to $\delta$, so the orthogonalized and non-orthogonalized QUE are equivalent). This means that we are testing for a scalar subvector of a low-dimensional vector in a product space, we should likely prefer the CLR approach to handling nuisance parameters. While the CLR's power properties are unclear when nuisance parameters interact with the constraints on the parameter of interest, as mentioned in subsection \ref{s: HighDimensional}, it is also harder to implement the orthogonalized version of the LR test in this case. Therefore, I would generally recommend only using the orthogonalization technique on a high-dimensional vector of nuisance parameters.
\section{Network Estimation Application}\label{s: Network application}
A common problem in the estimation of treatment effects is the presence of spillover effects, and this problem is especially challenging when the network over which spillovers occur is unobserved. There is a growing literature on using panel data to estimate these latent networks, including \cite{dePaulaRasulSouza2024}, \cite{Manresa2013}, \cite{ChernozhukovKarl2021}, \cite{MIAO2023}, \cite{BarigozziBrownless2019}, and \cite{Maung2022}, when the number of time periods is large. These methods generally use the covariance of the outcome $y$ across entities over time and the covariance of the covariates $x$ with $y$ across units to identify the network. Recently, \cite{ChernozhukovHuangWang2026} proposed a method for network estimation with panel data that leverages Neyman orthogonality and uses a multiplier bootstrap method for inference. However, their work assumes a linear model and focuses on inference for the latent network, rather than on the model's low-dimensional parameters as I will do here.
\par
\cite{dePaulaRasulSouza2024} show that when researchers have panel data with a large number of time periods, it is possible to point-identify the latent network, and they provide sufficient conditions under which both endogenous and exogenous social effect parameters are globally identified if the network is constant over time. This uses the model:
\[
y_{it} = \rho_0 W_{i,} y_t + \sum_{p=1}^P \alpha_{0,p} x_{it,p} + \sum_{p=1}^P \gamma_{0,p} W_{i,} x_{t,p} + u_i + \lambda_t + e_{it},
\]
and let $\theta = (\rho, \alpha, \gamma)$ and $\psi = vec(W)$. Here, the term involving $\rho$ captures the endogenous spillovers, $\alpha$ captures the direct effect of the $P$ covariates, and the term involving $\gamma$ captures the exogenous spillovers for each of these covariates. For simplicity, I focus on the case where there is a single exogenous covariate, so we have $ N^2$ moment conditions based on each entity's covariate $x_{it}$ being correlated with each entity's errors $\epsilon_{jt}$. Expressing them in a way that removes the fixed effects, the $Ni + j$ sample moment condition is equal to:
$$g_{Ni + j}(\theta,\psi) = \sum_{t=1}^T \ddot x_{it}( \ddot y_{jt} - \rho W_{j,} \ddot y_t - \alpha \ddot x_{jt} - \gamma W_{j,} \ddot x_{t})/T$$
where $\ddot x_{it,p}$ and $\ddot y_{it}$ denotes the double demean values of the data.\footnote{For the setting with an unconstrained $(\theta,\psi)$, \cite{basu2023} extends the Yule-Frisch-Waugh-Lovell Theorem to show that this GMM estimator will be identical to one that uses moment conditions including the fixed effects. Here, they may not be numerically identical due to the constraints, but the moment conditions will still identify $\theta$ under the same conditions.}
\par
\cite{dePaulaRasulSouza2024} provide identification results where they show that under several different restrictions on the parameter space, all parameters are point identified. In particular, in Corollary 3 they show that in a version of the model without the fixed effect the parameters are identified if $0 < \rho_0 < 1$ and $W_{j,i} \geq 0$ with at least one row of $W$ normalized to be one, but they also show that when time fixed effects are included, it is necessary to normalized all rows of $W$ to sum to one. These constraints on $W$ also have the nice feature that they allow $\rho$ and $\gamma$ to be interpreted as the total endogenous and exogenous spillover effect entities experience. The constraint that the endogenous spillovers are positive is also plausible in many contexts.\footnote{Note, however, that the constraint on $\rho$ was not originally imposed by \cite{dePaulaRasulSouza2024} during estimation. While they note that their estimated values satisfy their identification constraints, imposing those constraints during estimation is useful for achieving consistency.} I therefore use a GMM estimator with the constraints the $0 \leq \rho \leq 1$ and $W_{j,} \geq 0$ and $\sum_{i=1}^N W_{j,i} = 1$ for each $j$.
\par
While the focus of \cite{dePaulaRasulSouza2024} is primarily on identification, they also discuss an estimation method based on the adaptive elastic net GMM estimator of \cite{CanerZhang2014}. However, \cite{CanerZhang2014}'s asymptotic normality result for their elastic net GMM estimator requires the number of parameters and the number of moment conditions to grow more slowly than the sample size.\footnote{Also, if we impose the constraints on $W$ throughout the estimation, then its L1 norm is already fixed, so the elastic net estimator's L1 penalty does not affect the estimate. In \cite{dePaulaRasulSouza2024}, they implement the elastic net estimator as a multi-step estimator in which the L1 penalty does affect the solution, with the L1 norm of $W$ not constrained in the initial steps.} Here, I focus on using a constrained one-step GMM estimator as the initial constrained estimator. The simplex constraints tend to induce sparsity in the weights, similar to other contexts in which they are used such as synthetic controls (see e.g., \cite{Abadie2021}). Relatedly, \cite{dePaulaRasulSouza2024} note that small links in the network can be estimated as zeros due to the non-negativity constraint, and this, combined with the row-sum normalization, can lead to other links being overestimated. These errors affect the estimation of $(\rho,\alpha,\gamma)$ and can lead the GMM estimates for them to fail to be asymptotically normal as well. To correct this issue, I apply the method from \ref{s: HighDimensional} by orthogonalizing the moment with respect to the network's estimated weights and using these to obtain the quasi-unconstrained estimates $(\Tilde{\rho},\Tilde{\beta},\Tilde{\gamma})$. In Appendix D of the Supplementary Materials, I discuss the implementation of the constrained and quasi-unconstrained estimators in more detail.
\par
I apply the method by replicating the tax competition between U.S. states, as in \cite{BesleyCase1995} and \cite{dePaulaRasulSouza2024}, using \cite{dePaulaRasulSouza2024}'s extended dataset covering the 48 mainland US states ($N=48$) over the years 1962 to 2015 ($T=53$). The dependent variable, $\Delta\tau_{it}$, is the change in state tax liabilities between year $t$ and $t-2$. These tax liabilities include state income, sales, and corporate taxes. The covariate variable I include is income per capita, denoted $x_{it}$. \cite{dePaulaRasulSouza2024} additionally include variables for the unemployment rate and the proportions of young and elderly in the state's population, but I find similar estimates for $\rho$ when excluding these. There are two main theoretical reasons to expect endogenous spillovers in this context: factor mobility, in which labor and capital move in response to differential tax rates, and yardstick competition, in which voters compare their state's taxes to others to evaluate their own politicians' quality.
\par
\begin{table}[htpb]
\centering
\caption{US State Tax Spillovers}
\label{tab:empirical_results}\scalebox{.75}{
\begin{tabular}{lccc}
\hline\hline
& Endogenous Spillover ($\rho$) & Income Direct Effect ($\alpha_1$) & Income Exogenous Spillover ($\gamma_1$) \\
Constrained Estimate (CE) & 0.4860 & 0.0590 & 0.0240 \\
Orthogonalized Estimate (QUE) & 0.4780 & 0.0850 & 0.0070 \\
QUE CLR 95\% CI & [0.460, 0.498] & [0.007, 0.163] & [-0.032, 0.046] \\
QUE CLM 95\% CI & [0.459, 0.496] & [0.001, 0.163] & [-0.035, 0.046] \\
QUE-Wald 95\% CI & [0.459, 0.498] & [0.001, 0.168] & [-0.035, 0.049] \\
CE Wald 95\% CI & [0.460, 0.511] & [0.034, 0.083] & [0.000, 0.048] \\
M-Bootstrap 95\% CI & [0.459, 0.497] & [0.010, 0.167] & [-0.034, 0.048] \\
\hline
\end{tabular}}
\vspace{1ex}
\raggedright
\end{table}
Table \ref{tab:empirical_results} presents the empirical estimates, comparing confidence intervals constructed with a Wald test using the constrained estimate with those constructed using the proposed CLR, CLM, and Wald statistic. I also include a confidence interval constructed with a Multiplier Bootstrap using the orthogonalized moment conditions. We can see that the constrained and quasi-unconstrained estimates of the endogenous spillover are quite similar, and the inference methods agree that this effect is statistically significant.\footnote{Estimating the full specification using the Elastic Net estimator \cite{dePaulaRasulSouza2024} finds a similar estimate for the preferred endogenous spillover effect with $\hat{\rho} = 0.402 $ when treating the covariates as exogenous. When using instruments based on lagged values of geographic neighbors' covariates, they also obtain a similar estimate with $\hat{\rho} = 0.452$.} For the direct effect and exogenous spillover effect of income, whether the unorthogonalized, constrained estimate has more relevance for the statistical significance of the effects. With the direct effect, for most of the CI the lower bound is right by zero, indicating the estimates are only marginally significant, while this is not the case for the traditional Wald test, as the constrained estimate is about 4.6 times the size of the standard error. On the other hand, the traditional Wald test yields a marginally significant result for the exogenous spillover effect, whereas zero is well within the CI from the quasi-unconstrained inference methods.
\subsection{Simulations}
To evaluate the finite-sample size control of the proposed tests in a high-dimensional, constrained setting, I conduct Monte Carlo simulations calibrated to the empirical tax competition application. Since the fixed effects are removed by double demeaning the data, for the simulations, I simply remove the fixed effects and generate the data using the double-demeaned values of changes in tax rates and income per capita. I let each state's covariate data $x_{it}$ and errors $\epsilon_{it}$ be sampled independently and identically within a state. The changes in tax rates $\Delta \tau_{it}$ are then generated by solving the system of $N$ linear equations:
\begin{equation*}
\begin{pmatrix}
\Delta\tau_{1t} \\
... \\
\Delta\tau_{Nt}
\end{pmatrix}
=
\begin{pmatrix}
\rho_0 \sum_{j=1}^N W_{1j} \Delta\tau_{jt} + \alpha_{0} x_{1t} + \gamma_{0} \sum_{j=1}^N W_{1j} x_{jt} + \epsilon_{1t} \\
... \\
\rho_0 \sum_{j=1}^N W_{Nj} \Delta\tau_{jt} + \alpha_{0} x_{Nt} + \gamma_{0} \sum_{j=1}^N W_{Nj} x_{jt} + \epsilon_{Nt}
\end{pmatrix}
,
\end{equation*}
for each time period.
I use the values of $(\rho_0,\alpha_0,\gamma_0)$ estimated by the quasi-unconstrained estimator in the empirical application. For $W_0$, I also use the network estimated in the empirical application, but because I'm interested in how the performance changes as we change the size of $N$ and I want the results to not be dependent on a particular choice of network, I conduct the simulations with multiple different values of $N$ and in each draw of the simulations I choose $W_0$ by selecting a $N$ by $N$ subgraph of the estimated network uniformly as random.\footnote{The involves selecting $N$ states uniformly at random, selecting the subgraph for those states, and normalizing the edges so that each state's incoming edges still sum to one. If none of a selected state's incoming neighbors are included, I assign all of its weight to a randomly selected included state.}
\par
I focus on comparing the sizes of the standard Wald and Multiplier Bootstrap tests against the proposed quasi-unconstrained Wald, CLR, and CLM tests. I also examine the power of the tests in the same setting, with in both cases the nominal sizes of the tests being 0.05.\footnote{The null value is the true value for the sizes results, and the null value is 0 for the power results.} I use $T = 53$, like in the sample of \cite{dePaulaRasulSouza2024}. For the number of states, I have them equal to 5, 10, and 20, so smaller than in the empirical application, but still large enough to have the number of parameters and moments be significantly larger than the number of time periods.
\begin{table}[htb]\centering\caption{Main Size Results \label{TableSize}}\scalebox{.75}{
\begin{threeparttable}
\begin{tabular}{l *{5}{c}} \hline
& \multicolumn{1}{c}{QUE CLR} & \multicolumn{1}{c}{QUE CLM} & \multicolumn{1}{c}{QUE Wald} & \multicolumn{1}{c}{CE Wald} & \multicolumn{1}{c}{Multiplier Bootstrap} \\
\toprule
\multicolumn{6}{l}{\textbf{Testing $H_0: \rho = \rho_0$}} \\
\midrule
$N = 5$ & 0.147 & 0.146 & 0.148 & 0.215 & 0.160 \\
$N = 10$ & 0.072 & 0.072 & 0.072 & 0.116 & 0.076 \\
$N = 20$ & 0.063 & 0.062 & 0.061 & 0.256 & 0.061 \\
\hline
\multicolumn{6}{l}{\textbf{Testing $H_0: \alpha = \alpha_0$}} \\
\midrule
$N = 5$ & 0.100 & 0.101 & 0.098 & 0.137 & 0.106 \\
$N = 10$ & 0.083 & 0.081 & 0.083 & 0.098 & 0.089 \\
$N = 20$ & 0.087 & 0.088 & 0.087 & 0.088 & 0.088 \\
\hline
\multicolumn{6}{l}{\textbf{Testing $H_0: \gamma = \gamma_0$}} \\
\midrule
$N = 5$ & 0.045 & 0.044 & 0.044 & 0.081 & 0.050 \\
$N = 10$ & 0.052 & 0.049 & 0.051 & 0.136 & 0.058 \\
$N = 20$ & 0.079 & 0.079 & 0.078 & 0.230 & 0.075 \\
\hline
\end{tabular}
\begin{tablenotes}
\small
\item Notes: All simulations are conducted with a thousand replications.
\end{tablenotes}
\end{threeparttable}}
\end{table}
Looking at the size results in Table \ref{TableSize} for $\rho$, we can see that all three of the QUE inference methods, as well as the Multiplier Bootstrap, get close to the nominal size when $N = 10$ and $N = 20$, while the traditional Wald test has fairly substantial over-rejection. In general, the presence of constraints may lead tests that use the constrained estimate to either over- or under-reject, depending on the context. In the simplest case of a one-dimensional parameter constrained to lie on an interval like with the constraints on $\rho$, the asymptotic distribution of the constrained estimator will generally be a censored normal, leading to the traditional Wald test being overly conservative as discussed in section \ref{s: Power}. However, because of constraints on the network weights, the distribution of the constrained estimator of $\rho$ is not well approximated by a truncated normal distribution.
\begin{table}[htb]\centering\caption{Main Power Results \label{TablePower}}\scalebox{.75}{
\begin{threeparttable}
\begin{tabular}{l *{5}{c}} \hline
& \multicolumn{1}{c}{QUE CLR} & \multicolumn{1}{c}{QUE CLM} & \multicolumn{1}{c}{QUE Wald} & \multicolumn{1}{c}{CE Wald} & \multicolumn{1}{c}{Multiplier Bootstrap} \\
\toprule
\multicolumn{6}{l}{\textbf{Testing $H_0: \rho = 0$}} \\
\midrule
$N = 5$ & 0.782 & 0.707 & 0.705 & 0.868 & 0.716 \\
$N = 10$ & 0.905 & 0.867 & 0.868 & 0.962 & 0.871 \\
$N = 20$ & 0.985 & 0.966 & 0.966 & 0.997 & 0.967 \\
\hline
\multicolumn{6}{l}{\textbf{Testing $H_0: \alpha = 0$}} \\
\midrule
$N = 5$ & 0.931 & 0.928 & 0.931 & 1.000 & 0.930 \\
$N = 10$ & 0.991 & 0.991 & 0.991 & 1.000 & 0.991 \\
$N = 20$ & 0.998 & 0.997 & 0.997 & 1.000 & 0.998 \\
\hline
\multicolumn{6}{l}{\textbf{Testing $H_0: \gamma = 0$}} \\
\midrule
$N = 5$ & 0.038 & 0.038 & 0.037 & 0.066 & 0.041 \\
$N = 10$ & 0.034 & 0.036 & 0.037 & 0.112 & 0.044 \\
$N = 20$ & 0.055 & 0.057 & 0.054 & 0.151 & 0.057 \\
\hline
\end{tabular}
\begin{tablenotes}
\small
\item Notes: All simulations are conducted with a thousand replications.
\end{tablenotes}
\end{threeparttable}}
\end{table}
For the direct effect, all the methods are consistent is experiencing a modest amount of over-rejection. The similar performance of the constrained and quasi-unconstrained methods for $\alpha$ may be due to the direct effect being least sensitive to the estimated spillover network, so it benefits the least from orthogonalizing with respect to the network. For $\gamma$, the QUE methods and the Multiplier Bootstrap generally control size well with only some modest over-rejection when $N = 20$, while the over-rejection of the traditional Wald test gradually gets worse as $N$ grows. For the three parameters, the similarity in performance of the CLR, CLM, and Wald test with the QUE is unsurprising, but the fact that the Multiplier Bootstrap test has very similar performance as well suggests the other methods that make use of orthogonalization and don't impose constraints when using the orthogonalized moment conditions can also be well-suited to contexts with both high-dimensional nuisance parameters and constraints.
The power results for testing the parameters are equal to zero are in Table \ref{TablePower}. The different QUE tests and the multiplier bootstrap generally have very similar power for testing $H_0: \alpha = 0$ and $H_0: \gamma = 0$, where for $H_0: \rho = 0$ they also have similar power, but with the CLR test having a slight edge. This pattern is unsurprising given that $\alpha$ and $\gamma$ are unconstrained and $\rho$ is constrained but not especially close to the boundary. The power of the tests generally increases with $N$, although this pattern is hardly detectable for $\gamma$, likely because of how close to zero $\gamma_0$ is.
\subsection{Limitations and Open Questions}\label{s: Discussion}
While the simulations contain parameters close to or on the boundary due to the network weights, the non-negativity constraint for $\rho$ is only binding a small percentage of the simulations. It is worth highlighting this because while the method is relatively well-suited to handling boundary cases of the network constraints, the same cannot be said for the constraint on $\rho$. Table \ref{TableWeakSize} contains results for testing $\rho = 0$ when this hypothesis is true under otherwise the same conditions as the previous simulations. We can see that all inference methods now perform poorly, with over-rejection still the worst for the traditional Wald test.
\begin{table}[htb]\centering\caption{Weakly-Identified Size Results \label{TableWeakSize}}\scalebox{.9}{
\begin{threeparttable}
\begin{tabular}{l *{5}{c}} \hline
& \multicolumn{1}{c}{QUE CLR} & \multicolumn{1}{c}{QUE CLM} & \multicolumn{1}{c}{QUE Wald} & \multicolumn{1}{c}{CE Wald} & \multicolumn{1}{c}{Multiplier Bootstrap} \\
\toprule
\multicolumn{6}{l}{\textbf{Testing $H_0: \rho = 0$}} \\
\midrule
$N = 5$ & 0.207 & 0.134 & 0.135 & 0.233 & 0.149 \\
$N = 10$ & 0.456 & 0.328 & 0.327 & 0.604 & 0.338 \\
$N = 20$ & 0.815 & 0.708 & 0.709 & 0.916 & 0.714 \\
\hline
\end{tabular}
\begin{tablenotes}
\small
\item Notes: All simulations are conducted with a thousand replications.
\end{tablenotes}
\end{threeparttable}}
\end{table}
This is due to the limitation of the method that the initial estimator $\hat{\theta}$ must be $\sqrt{n}$-consistent. When Neyman orthogonality is used, $\sqrt{n}$-consistency is not required for the entire vector, but the nuisance-parameter estimator must satisfy rate conditions tied to the sensitivity of the derivatives of the criterion function in order for the adaptivity condition in equation \eqref{eq: adaptivity} to hold. This typically requires imposing structure on the nuisance parameter, for example, sparsity of the spillover network. These convergence rate conditions can also fail when proximity to the boundary coincides with weak identification. If $\rho$ and $\gamma$ are both close to zero, the latent network $W$ becomes weakly identified and therefore cannot be consistently estimated.
\par
Importantly, these size distortions are not simply due to a parameter being close to the boundary. To see this, Table \ref{TableNonNegativeSize} contains results under the same conditions as the previous subsection, except with all parameters constrained to be non-negative and the true value of the direct effect being zero. In this case, the QUE tests for $H_0 : \alpha = 0$ actually control size better than in the previous results, and the traditional Wald test also controls size, although it is slightly conservative due to the left tail being truncated. Together, these results help illustrate that the method is well-suited to handling boundary cases for both low-dimensional and high-dimensional parameters, except when the boundary indicates a point of identification failure. In such cases, we would ideally like a method that accommodates nuisance parameters that are weakly identified, constrained, and high-dimensional, and this is an interesting direction for future work.
\begin{table}[htb]\centering\caption{Non-Negative Direct Effect Size Results \label{TableNonNegativeSize}}\scalebox{.9}{
\begin{threeparttable}
\begin{tabular}{l *{5}{c}} \hline
& \multicolumn{1}{c}{QUE CLR} & \multicolumn{1}{c}{QUE CLM} & \multicolumn{1}{c}{QUE Wald} & \multicolumn{1}{c}{CE Wald} & \multicolumn{1}{c}{Multiplier Bootstrap} \\
\toprule
\multicolumn{6}{l}{\textbf{Testing $H_0: \alpha = 0$}} \\
\midrule
$N = 5$ & 0.036 & 0.038 & 0.038 & 0.022 & 0.041 \\
$N = 10$ & 0.047 & 0.047 & 0.049 & 0.025 & 0.048 \\
$N = 20$ & 0.065 & 0.066 & 0.067 & 0.034 & 0.063 \\
\hline
\end{tabular}
\begin{tablenotes}
\small
\item Notes: All simulations are conducted with a thousand replications.
\end{tablenotes}
\end{threeparttable}}
\end{table}
\par
Several other theoretical questions remain open, including characterizing when the CLR test is similar and admissible beyond the cases discussed in Section \ref{s: Power}. Also, when using the orthogonalized QUE, one could either orthogonalize with respect to all nuisance parameters and conduct inference using one of the methods from section \ref{s: Inference} or orthogonalize with respect to only the high-dimensional nuisance parameter vectors and one of the methods from section \ref{s: InferenceForSubvectors} could be employed. Here, I focused on the second approach, where you simply have a single orthogonalized estimate you can use to conduct inference for multiple parameters (e.g., orthogonalizing with respect to the elements of $W$ to conduct inference for $\rho$, $\alpha$, and $\gamma$). In addition to being simpler, one might expect the sufficient-statistic method to deliver superior power to the orthogonalization approach when handling a low-dimensional nuisance parameter, since it includes more information in the likelihood. Showing this formally would be an interesting area for future investigation.
\newpage
\noindent
\textbf{Acknowledgments}: I am grateful for helpful comments and feedback from Adam McCloskey, Carlos Martins-Filho, Philipp Ketz, and Xiaodong Liu.
\footnotesize
\bibliographystyle{econ-aea}
\bibliography{references}
\newpage