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.
102,098 characters
Continuously Updating GMM in Linear IV Models: A Polynomial Approach
\title{Continuously Updating GMM in Linear IV Models:\\
A Polynomial Approach}
\author{Marcelo J. Moreira\thanks{
FGV EPGE. Email: [email removed].} \and Whitney K. Newey\thanks{
MIT Department of Economics. Email: [email removed].} \and Mahrad
Sharifvaghefi\thanks{
University of Pittsburgh. Email: [email removed].}}
\date{\today }
\maketitle
\begin{abstract}
This paper characterizes the global minimum of the continuously updating
generalized method of moments (CU-GMM) objective in linear instrumental
variables models. We allow optimal weighting matrices under
heteroskedasticity, autocorrelation, or clustering. We show that the
objective is a ratio of polynomials. For one endogenous regressor,
stationary points of CU-GMM objective function are real eigenvalues of a
companion matrix. Comparing their objective values with the value at
infinity gives the global minimum, extending the classical eigenvector
approach to limited information maximum likelihood. With multiple endogenous
regressors, algebraic elimination and checks for real solutions identify the
minimum among finitely many candidate objective values, including boundary
values. Galois theory rules out general formulas by radicals even with one
endogenous regressor and two instruments, while numerical root finding
remains possible. Finding the global minimum allows us to compute
overidentification and likelihood ratio tests, including the conditional
likelihood ratio (CLR) test.
\end{abstract}
\noindent\textbf{Keywords:} continuously updating GMM, instrumental
variables, weak instruments, conditional likelihood ratio, polynomial roots,
companion matrix, Groebner basis, homotopy continuation.
\noindent\textbf{JEL codes:} C12, C26, C36.
\section{Introduction}
In the classical linear instrumental variables model, conditional
homoskedasticity means that the reduced-form errors have a common variance
matrix that does not depend on the instruments. When these errors are also
conditionally uncorrelated across observations, the variance matrix of the
reduced-form sample moments has a Kronecker product structure. Using the
usual reduced-form residual variance estimator, the continuously updating
GMM estimator (CUE) coincides with the limited information maximum
likelihood estimator (LIML) of \citet{AndersonRubin49}; see also
\citet{HansenHeatonYaron96}. The objective is a generalized Rayleigh
quotient. Its minimum is the smallest generalized eigenvalue, and a
corresponding eigenvector gives the coefficient estimates after
normalization; see \citet{MagnusNeudecker88}. Section S-2
in the online supplement presents this familiar characterization, including the
case in which the minimizing direction is at infinity.
The Kronecker product structure need not hold with heteroskedastic,
autocorrelated, or clustered errors. The GMM framework of \citet{Hansen1982}
allows more general variance matrices, and \citet{HansenHeatonYaron96}
introduce CUE by updating the weighting matrix at each candidate parameter
value. Without the Kronecker structure, the variance matrix of the moments
can change in both shape and scale as the parameter varies. The objective
need not be a generalized Rayleigh quotient, and local optimization may fail
to find its global minimum. This computational problem is particularly
relevant with weak or many instruments, when conventional asymptotic
approximations can be unreliable. \citet{StaigerStock97} study two-stage
least squares and LIML under weak instruments. Under homoskedasticity,
\citet{Bekker94} develops asymptotic approximations that allow the number of
instruments to grow with the sample size. \citet{NeweyWindmeijer09} analyze
continuously updating GMM with many weak moment conditions and
heteroskedastic errors. \citet{NeweySmith04} compare the higher-order
properties of GMM and generalized empirical likelihood estimators, including
CUE. These statistical analyses motivate reliable computation of the global
minimum that defines CUE.
This paper provides an algebraic characterization of that minimum. We first
consider a scalar structural coefficient and a fixed positive definite
estimate of the joint reduced-form variance matrix. The corresponding
variance matrix of the sample moments is quadratic in the coefficient, so
the continuously updating objective is a ratio of polynomials. We derive a
polynomial equation for its finite stationary points and bound its degree in
terms of the number of instruments. Evaluating the objective at the real
roots and comparing these values with its limit at infinity characterizes
the global minimum. This characterization does not require starting values
for local optimization or an arbitrarily bounded search region.
Computationally, finding these polynomial roots corresponds to finding the
eigenvalues of a companion matrix. This connects the general problem to the
classical LIML calculation, but the eigenvalues play different roles. In the
homoskedastic problem, the smallest generalized eigenvalue is the minimized
objective, and the estimator is recovered from its eigenvector. In the
companion problem, the real eigenvalues are candidate coefficient values. A
finite CUE is one of these eigenvalues, selected by evaluating the original
objective and comparing it with its limiting value at infinity.
The classical algebraic characterizations of constrained extrema in
\citet[Chapters~11 and~17]{MagnusNeudecker88} do not provide this global
minimization over the structural coefficient. Variational eigenvalue
formulas, Fischer's min--max theorem, and the trace and determinant extrema
obtained from Poincar\'{e}'s separation theorem concern fixed quadratic
forms under normalization or orthogonality restrictions. Their use requires
a reduction to that structure, which a general CU-GMM objective does not
provide. Matrix inequalities, including Cauchy--Schwarz, H\"{o}lder, and
Minkowski, and variational representations of trace powers and determinants
provide bounds or alternative representations. However, their equality
conditions do not by themselves select the structural coefficient that
globally minimizes CU-GMM. Our results supply an algebraic characterization
of this global minimum over the structural coefficient, together with
methods for computing it.
\citet{HansenLee21} study convergence of iterated GMM and inference under
misspecification. Each iteration minimizes a criterion with the weight held
fixed at the preceding estimate; the limit is a fixed point of this map.
Their population contraction result assumes a compact parameter space,
unique interior minima for the fixed-weight problems, smooth moments and
weights, uniformly positive definite weights, and misspecification
sufficiently small relative to curvature and weight sensitivity. With
independent observations and further moment and smoothness conditions, the
sample map is a contraction with probability tending to one. Their
asymptotic normality and inference results also require a nonsingular
identification matrix and do not cover weak identification.
This fixed-point characterization concerns a different estimator. CUE varies
the moments and the weight simultaneously, so its first-order conditions
generally contain terms absent from the iterated-GMM conditions.
Consequently, convergence of iterated GMM need not deliver CUE, even when
the contraction conditions hold. Our result characterizes the global CUE
minimum for the observed sample without requiring a contraction or strong
identification. The singular-variance extension discussed below also falls
outside their positive-definiteness conditions.
The polynomial representation also permits us to ask whether a general
formula using arithmetic operations and radicals can replace numerical root
finding. Polynomial equations of degree at most four admit such formulas,
whereas higher degrees require an analysis of the Galois group. Degree alone
does not establish impossibility for a particular family of objectives.
Proposition S.3 constructs a CU-GMM example with a scalar
coefficient and two instruments whose stationary polynomial has degree six,
is irreducible, and has the full symmetric Galois group. Every global
minimizer in this example is finite and cannot be expressed by radicals.
Thus no formula by radicals can produce CUE for every admissible variance
matrix, even in this small model. Exact descriptions by polynomial equations
and exact algorithms for isolating real roots remain available.
Numerical eigenvalue computation is already part of the classical approach.
In the homoskedastic model, the order of the generalized eigenvalue problem
is one greater than the number of endogenous regressors. In the full-rank
case, formulas by radicals therefore cover models with at most three
endogenous regressors. With four or more, the characteristic polynomial has
degree at least five, and no general formula by radicals is available for an
unrestricted eigenvalue problem of that order. The relevant dimension is
that of the structural coefficient, not the number of instruments. The
eigenvector characterization itself remains valid in every dimension.
Likewise, the value of the companion representation is that it reduces
computation to a well-studied problem with established numerical methods;
see \citet{AndersonEtAl1999}.
Finding the global minimum also allows inference. \citet{Moreira03} develops
the conditional likelihood ratio (CLR) test for linear IV models, and
\citet{AndrewsMoreiraStock06} study its power relative to the envelope for
invariant similar tests in the homoskedastic Gaussian model.
\citet{Kleibergen05} develops identification-robust GMM tests, while
\citet{AndrewsMikusheva16} give a conditional inference framework for
general moment condition models. For the likelihood criterion considered
here, the LR statistic subtracts the global minimum from the objective
evaluated at the null. The same minimization must be performed for the
simulated processes used to obtain conditional critical values.
\citet{MoreiraMoreira19} study conditional testing with heteroskedastic and
autocorrelated errors. In this setting, computing the full likelihood
criterion is relevant to the power comparisons in
\citet{MoreiraRidderSharifvaghefi26}. The consequences of replacing a
continuous search by a finite grid are studied by
\citet{MoreiraSharifvaghefi26}. Our results instead solves the minimization
step.
The same minimum is the overidentification statistic for testing whether the
moment restrictions hold for some structural coefficient. Under conventional
strong-identification conditions, its chi-square approximation has degrees
of freedom equal to the number of overidentifying restrictions; see
\citet{Hansen1982}. Weak identification can invalidate this calibration, but
we can use a conservative chi-square critical value equal to the number of
moments. A local minimum can overstate the overidentification statistic and
mislead the test to reject too often. Computing the global minimum therefore
matters for both parameter inference and specification testing.
The scalar characterization extends to singular variance matrices when the
objective is defined with the Moore--Penrose inverse. \citet{NeweySmith04}
already allow generalized inverses in their definition of CUE. We show how
to minimize the Moore--Penrose criterion globally when rank can vary with
the coefficient. On intervals of constant rank, the objective is rational.
Its values where rank falls must be evaluated separately, along with the
limit at infinity, since they need not agree with the limiting rational
expression; see Section \ref{sec:singular-univariate}.
We also establish an analogue of partialling out in homoskedastic LIML. For
a fixed positive definite full reduced-form variance estimate, minimizing
over the coefficients on included exogenous covariates gives exactly the
CU-GMM criterion based on projected moments, provided the variance estimate
is projected correspondingly. The two formulations yield the same finite
minimizing structural coefficients and the same infimum. Thus covariates do
not add dimensions to the remaining global search; see Section S-1.
The contribution is to connect an important, well-known econometric
optimization problem to theoretical and computational algebra. This
connection allows computational methods for CU-GMM to be assessed against
established methods for solving polynomial equations. The polynomial
representation identifies all candidates for the scalar minimum, companion
matrices provide a numerical implementation. Polynomial interpolation
recovers the coefficients without symbolic matrix inversion. With several
endogenous regressors, the first-order conditions become a system of
polynomial equations. We discuss methods for this extension, including the
treatment of stationary components and boundary directions. The multivariate
problem requires additional computation beyond the scalar root calculation.
Section \ref{sec:setup} introduces the model and its statistical
applications. Section \ref{sec:implementation} develops the polynomial
characterization and computation. Section \ref{sec:multivariate} treats
multiple endogenous regressors. Section \ref{sec:simulation} presents
numerical results. The online supplement contains the homoskedastic
benchmark, interpolation results, and proofs.
\section{Problem setup}
\label{sec:setup}
Let $y_1$ be the $n\times1$ outcome vector, $y_2$ be the $n\times1$ vector
of observations on an endogenous regressor, and $X$ be an $n\times p$ full-column-rank matrix
of observations on $p$ exogenous covariates. The
structural equation is
\begin{equation}
y_1=y_2\beta^*+X\delta^*+u. \label{eq:structural}
\end{equation}
The parameter of interest is $\beta^*$. Let $\widetilde{Z}$ be an $n\times k$
matrix of observations on $k$ instrumental variables (IVs). The reduced form
for $y_2$ is
\begin{equation}
y_2=\widetilde{Z}\pi+X\varphi+v_2. \label{eq:rf-y2}
\end{equation}
Substituting \eqref{eq:rf-y2} into \eqref{eq:structural} gives
\begin{equation}
y_1=\widetilde{Z} \pi \beta^{\ast} + X \gamma^{\ast}+v_1, \qquad
\gamma^{\ast}=\delta^{\ast}+\varphi \beta^{\ast}, \qquad
v_1=u+v_2\beta^{\ast}. \label{eq:rf-y1}
\end{equation}
Writing $Y=(y_1,y_2)$, $V=(v_1,v_2)$, $a^{\ast}=(\beta^{\ast},1)^{\prime}$,
and $\Gamma=(\gamma^{\ast}, \varphi)$, the reduced form is
\begin{equation}
Y=\widetilde{Z}\pi a^{*\prime}+X\Gamma+V. \label{eq:matrix-rf}
\end{equation}
Let
\begin{equation*}
Z=M_X\widetilde{Z}, \qquad
M_X=I_n-X\left(X^{\prime}X\right)^{-1}X^{\prime}.
\end{equation*}
Each column of $Z$ contains the variation in the
corresponding instrumental variable that is orthogonal
to the exogenous covariates in $X$. The instrumental-variable exclusion restrictions are $\mathbb{E}
[Z^{\prime}u]=0$. Because $Z^{\prime}X=0$, the sample moments for $\beta$
can be written as
\begin{align}
g_n(\beta) &=n^{-1/2}Z^{\prime}(y_1-\beta y_2) =n^{-1/2}Z^{\prime}Yb(\beta),
\label{eq:moment} \\
&=\left(b(\beta)^{\prime}\otimes I_k\right) \operatorname{vec}\left(n^{-1/2}Z^{
\prime}Y\right), \notag
\end{align}
where $b(\beta)=(1,-\beta)^{\prime}$.
Let $\Sigma_n$ denote the conditional variance matrix of
$\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$ given
$\widetilde{Z}$ and $X$. Then, the conditional variance matrix of
$g_n(\beta)$ is
\begin{equation}
\Omega_n(\beta)
=\left(b(\beta)^{\prime}\otimes I_k\right)\Sigma_n
\left(b(\beta)\otimes I_k\right).
\end{equation}
Let $\widehat{\Sigma}_{n}$ be a consistent estimator of $\Sigma_{n}$; see
\citet{MacKinnonWhite85} for variance estimation under heteroskedasticity,
\citet{NeweyWest1987} for HAC variance estimation, and
\citet{CameronGelbachMiller11} for multiway clustering.
Then, a consistent estimator of $\Omega_{n}(\beta)$ is
\begin{equation}
\widehat{\Omega}_n(\beta) =\left(b(\beta)^{\prime}\otimes I_k\right)
\widehat{\Sigma}_{n} \left(b(\beta)\otimes I_k\right). \label{eq:omega-hat}
\end{equation}
To simplify the presentation, we first assume
\begin{assumption}
\label{ass:sigma-pd} The matrix $\widehat{\Sigma}_n$ is symmetric positive
definite.
\end{assumption}
Section \ref{sec:singular-univariate} extends the analysis to symmetric
positive semidefinite variance matrices, using the Moore--Penrose inverse.
Assumption \ref{ass:sigma-pd} implies that $\widehat{\Omega}_n(\beta)$ is
positive definite for every $\beta\in\mathbb{R}$, as established by Lemma
S.2. Therefore, under this assumption, the CU-GMM
objective can be written as
\begin{equation} \label{eq:cue-objective}
Q_n(\beta) =g_n(\beta)^{\prime} \widehat{\Omega}_n(\beta)^{-1} g_n(\beta).
\end{equation}
Section S-1 in the online supplement shows
that minimizing this partialled-out CU-GMM objective gives the same finite
estimates of $\beta$ as jointly minimizing the corresponding full CU-GMM
objective over the endogenous and exogenous coefficients.
In the linear IV model, the CU-GMM objective function depends on $\beta$
only through the vector $b(\beta)$. So, we can write $Q_n(\beta)=Q_n\left(b(
\beta)\right)$. For $\beta\in\mathbb{R}$, let $c(\beta)=1/\left\|
b(\beta)\right\|>0$ and define
\begin{equation}
\overline{b}(\beta)=c(\beta)b(\beta)=\left(\frac{1}{\left\| b(\beta)\right\|}
,-\frac{\beta}{\left\| b(\beta)\right\|}\right)^{\prime}.
\end{equation}
Note that $Q_n\left(b(\beta)\right)$ is homogeneous of degree zero in $
b(\beta)$; that is, for any constant $c\neq0$, $Q_n\left(c\,b(\beta)
\right)=Q_n\left(b(\beta)\right)$. Therefore, $Q_n\left(b(\beta)\right)=Q_n
\left(\overline{b}(\beta)\right)$.
As $\beta\to\pm\infty$, $\overline{b}(\beta)\to\mp e_2$,
where $e_2=(0,1)^{\prime}$. Therefore, the normalized sample moments
converge to $\mp n^{-1/2}Z^{\prime}y_2$, and their
estimated variance matrix converges to
$\widehat{\Sigma}_{22,n}$, the lower-right $k\times k$ block
of $\widehat{\Sigma}_n$. Under Assumption \ref{ass:sigma-pd}, $\widehat{\Sigma}_{22,n}$ is positive definite.
By continuity of matrix inversion,
\[
\lim_{\beta\to-\infty}Q_n(\beta)
=\lim_{\beta\to+\infty}Q_n(\beta)
=\frac{y_2^{\prime}Z\widehat{\Sigma}_{22,n}^{-1}
Z^{\prime}y_2}{n}.
\]
We define
\begin{equation}\label{eq:cu-gmm-end_point-value}
Q_n(-\infty)=Q_n(+\infty)
=\frac{y_2^{\prime}Z\widehat{\Sigma}_{22,n}^{-1}
Z^{\prime}y_2}{n}.
\end{equation}
By Lemma S.2,
$\widehat{\Omega}_n(\beta)$ is positive definite for every
finite $\beta$. Therefore, continuity of the sample moments and matrix
inversion implies that $Q_n(\beta)$ is continuous
on $\mathbb{R}$. Together with the endpoint limits established
above, this gives a continuous extension to
$\overline{\mathbb{R}}=\mathbb{R}\cup\{-\infty,+\infty\}$.
Compactness of this space and the extreme value theorem
yield the following result.
\begin{lemma}
\label{lem:cue-existence} Under Assumption \ref{ass:sigma-pd}, the CU-GMM
objective $Q_n(\beta)$, with endpoint values defined by
\eqref{eq:cu-gmm-end_point-value}, is continuous on the extended real line $
\overline{\mathbb{R}}$ and attains its minimum there.
\end{lemma}
While the global minimum may occur at the infinite endpoints, this event has
probability zero if the conditional distribution of $\operatorname{vec}
\left(n^{-1/2}Z^{\prime}Y\right)$ given $\widehat{\Sigma}_n$ is absolutely
continuous with respect to Lebesgue measure almost surely; see Proposition
S.2 in the online supplement.
In the just-identified case, where $k=1$, the minimum value of the CU-GMM
objective $Q_n(\beta)$ over $\overline{\mathbb{R}}$ is zero. If $
Z^{\prime}y_2\neq0$, the continuously updating estimator (CUE) that
minimizes $Q_n(\beta)$ over $\overline{\mathbb{R}}$ is
\begin{equation*}
\widehat{\beta}_{\mathrm{CUE}}=\frac{Z^{\prime}y_1}{Z^{\prime}y_2},
\end{equation*}
so that $g_n\left(\widehat{\beta}_{\mathrm{CUE}}\right)=0$ and $Q_n\left(
\widehat{\beta}_{\mathrm{CUE}}\right)=0$. Since the objective is
nonnegative, this is its global minimum. If $Z^{\prime}y_2=0$, the minimum
is also zero because $Q_n(\pm\infty)=0$.
In the rest of the paper, we focus on the overidentified case, where $k>1$.
Finding the global minimum of the CU-GMM objective over $\overline{\mathbb{R}
}$ generally requires solving a non-convex optimization problem. The main
text considers a general setting that allows for heteroskedasticity,
clustering, or autocorrelation in the errors. Section S-2
in the online supplement discusses the simpler homoskedastic case.
The global minimum determines three statistical objects: the CUE,
overidentification statistic, and the LR statistic underlying the CLR test.
The CUE minimizes $Q_n(\beta)$ over $\overline{\mathbb{R}}$:
\begin{equation}
Q_n\left(\widehat{\beta}_{\mathrm{CUE}}\right) =\min_{\beta\in\overline{
\mathbb{R}}}Q_n(\beta). \label{eq:cue-estimator}
\end{equation}
When the CU-GMM objective uses the homoskedastic variance specification, CUE
coincides with LIML; see Section S-2 in the online
supplement. With a more general robust variance estimator, the
two estimators need not coincide, even when the underlying errors are
homoskedastic.
Under standard regularity conditions, with a fixed number of valid
instruments and strong identification, 2SLS, LIML, and CUE are consistent
and asymptotically normal. Their properties
differ when the number of instruments grows with the sample size. Under the
homoskedastic many-instrument asymptotics of \citet{Bekker94}, with $
k/n\to\alpha\in(0,1)$ and the maintained identification conditions, 2SLS can
be inconsistent, whereas LIML remains consistent and asymptotically normal.
Allowing for heteroskedasticity and many weak moment conditions,
\citet{NeweyWindmeijer09} establish consistency of CUE when $k^2/n\to0$, and
asymptotic normality when $k^3/n\to0$. These results also require their
regularity and identification conditions.
The $J$-statistic for overidentification tests can be written as
\begin{equation}
J_n=\min_{\beta\in\overline{\mathbb{R}}}Q_n(\beta). \label{eq:overid-stat}
\end{equation}
The minimum is attained by Lemma \ref{lem:cue-existence}. The null
hypothesis is that the moment restrictions hold for some finite coefficient $
\beta^*$. Under the usual strong-identification regularity conditions, with
a fixed number $k>1$ of instruments, $J_n\overset{d}{\longrightarrow}
\chi^2_{k-1}$. The usual overidentification test therefore compares $J_n$
with a $\chi^2_{k-1}$ critical value. The loss of one degree of freedom
reflects estimation of the scalar coefficient $\beta$; see \citet[Section
3]{NeweyWindmeijer09}.
Under weak identification, this $\chi^2_{k-1}$ approximation need not be
reliable. A conservative alternative compares $J_n$ with a $\chi^2_k$
critical value. To justify this comparison, suppose that the null sample
moments satisfy a central limit theorem with a positive definite limiting
variance matrix and that the variance estimator is consistent along the
sequence considered. For fixed $k$, these conditions give the
Anderson--Rubin null limit
\begin{equation*}
AR_n(\beta^*)=Q_n(\beta^*) \overset{d}{\longrightarrow}\chi^2_k,
\end{equation*}
without requiring strong identification. Since $J_n\leq Q_n(\beta^*)$,
letting $\chi^2_{k,1-\alpha}$ denote the $1-\alpha$ quantile of $\chi^2_k$,
we have
\begin{equation}
\Pr\left(J_n>\chi^2_{k,1-\alpha}\right)
\leq\Pr\left(Q_n(\beta^*)>\chi^2_{k,1-\alpha}\right) \longrightarrow\alpha.
\label{eq:overid-size-bound}
\end{equation}
Thus the test has asymptotic rejection probability at most $\alpha$ along
any null sequence for which this AR approximation holds. Uniform validity
over a class of null models requires a corresponding uniform AR
approximation. Under strong identification, the larger $\chi^2_k$ critical
value makes this test conservative relative to the usual $\chi^2_{k-1}$ test.
Suppose that, conditional on the instruments and exogenous
covariates, $\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$ is normally
distributed with known positive definite variance matrix $\Sigma_n$.
Using $\Sigma_n$ in the CU-GMM objective, the likelihood ratio
statistic underlying the CLR test of $H_0:\beta^*=\beta_0$ is
\begin{equation}\label{eq:lr-stat}
LR_n(\beta_0)=Q_n(\beta_0)-J_n.
\end{equation}
Under weak-instrument asymptotics, consider the sequence
$\pi_n=\mu/\sqrt{n}$, where $\mu\neq0$ is the local first-stage
parameter. The null distribution of $LR_n(\beta_0)$ generally
depends on $\mu$, which cannot be consistently estimated.
Let $a_0=(\beta_0,1)^{\prime}$ and define
\begin{equation}
T_n=
\left[
(a_0^{\prime}\otimes I_k)\Sigma_n^{-1}(a_0\otimes I_k)
\right]^{-1/2}
(a_0^{\prime}\otimes I_k)\Sigma_n^{-1}
\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right).
\end{equation}
Under $H_0$, $T_n$ is sufficient for $\mu$, so the conditional
null distribution of the LR statistic given $T_n$ does not
depend on $\mu$. Let $c_{\alpha}(t;\Sigma_n)$ denote the
$(1-\alpha)$ quantile of this distribution given $T_n=t$.
Following the conditioning argument of \citet{Moreira03},
this conditional critical value controls the size of the
CLR test at level $\alpha$.
The conditional critical value
$c_{\alpha}(T_{n};\Sigma_n)$
can be well approximated using the simulation method described
in Section \ref{sec: MC-CLR}. This requires repeated solution
of the global minimization problem for simulated data.
In general, $\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$ need not
be normally distributed in finite samples, and $\Sigma_n$
is unknown. Under standard regularity conditions,
$\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$ is asymptotically
normal, and replacing $\Sigma_n$ with a consistent estimator
$\widehat{\Sigma}_n$ gives a feasible CLR test with
asymptotic size control.
The global minimizers determine the possible CUE estimates, while the
minimum value determines the overidentification statistic and enters both
the observed LR statistic and the simulated statistics used to obtain its
conditional critical value. Evaluating the objective at a local minimum or
over a restricted grid can overstate $J_n$ and lead to additional rejections
of the overidentifying restrictions. The approximation need not be bounded
above by $Q_n(\beta^*)$, so the size guarantee in
\eqref{eq:overid-size-bound} does not automatically carry over.
In HAC IV models, \citet{MoreiraRidderSharifvaghefi26} show that the CLR
test attains the decision-theoretic power frontier, whereas LM and CQLR
procedures can have power arbitrarily close to size in designs where the
null and alternative distributions are well separated. At the same time,
\citet{MoreiraSharifvaghefi26} show that the grid search approximation of
the global minimum value of CU-GMM objective function does not preserve the
size and power properties of the CLR test uniformly over their class of
data-generating processes, even when the number of grid points grows
polynomially with the sample size.
Computing the CUE, the overidentification statistic, and the CLR test
requires a reliable algorithm for global minimization. With a general
variance matrix, $Q_n(\beta)$ need not be convex, so local optimization and
grid search do not in general guarantee the global minimum. The next section
shows that, with one endogenous regressor, $Q_n(\beta)$ is a ratio of
polynomials. Thus, its global minimum can be computed by finding all real
roots of the polynomial equation obtained from the first-order condition and
then evaluating the CU-GMM objective at those roots and at either infinite
endpoint. This procedure makes the computed CUE, overidentification
statistic, and CLR statistic coincide with the objects defined above, up to
numerical error in coefficient recovery and root finding.
\section{Global minimization of the CU-GMM objective}
\label{sec:implementation}
Section \ref{sec:poly-ratio-representation} expresses the CU-GMM objective
as a ratio of polynomials. Section \ref{sec:poly-approach} characterizes its
global minimum using the real roots of the first-order condition and the
objective value at infinity. Section \ref{sec:companion} explains why
finding these roots generally requires numerical computation and shows how
to obtain them through an eigenvalue problem. Section \ref
{sec:singular-univariate} extends these results to positive semidefinite
variance matrices using the Moore--Penrose inverse.
\subsection{Polynomial ratio representation}
\label{sec:poly-ratio-representation}
The vector of sample moments $g_n(\beta)=n^{-1/2}Z^{\prime}\left(y_1-\beta
y_2\right)$, defined in Section 2, is an affine function of $\beta$.
Moreover, by partitioning $\widehat{\Sigma}_n$ as
\begin{equation*}
\widehat{\Sigma}_n=
\begin{pmatrix}
\widehat{\Sigma}_{11,n} & \widehat{\Sigma}_{12,n} \\
\widehat{\Sigma}_{21,n} & \widehat{\Sigma}_{22,n}
\end{pmatrix}
,
\end{equation*}
where each block is $k\times k$, we can write the estimated variance matrix
of the sample moments as
\begin{equation} \label{eq:sample_moment_variance_expanded}
\widehat{\Omega}_n(\beta) =\widehat{\Sigma}_{11,n} -\beta\left(\widehat{
\Sigma}_{12,n}+\widehat{\Sigma}_{21,n}\right) +\beta^2\widehat{\Sigma}
_{22,n}.
\end{equation}
Each element of $\widehat{\Omega}_n(\beta)$ is a polynomial in $\beta$ of
degree at most two. Under Assumption \ref{ass:sigma-pd}, the adjugate
formula gives
\begin{equation}
\widehat{\Omega}_n(\beta)^{-1} =\frac{\operatorname{adj}\left\{\widehat{\Omega}
_n(\beta)\right\}} {\det\left\{\widehat{\Omega}_n(\beta)\right\}}.
\end{equation}
The CU-GMM objective can therefore be written as
\begin{equation} \label{eq:cue-objective-adjugate}
Q_n(\beta) =\frac{g_n(\beta)^{\prime}\operatorname{adj}\left\{\widehat{\Omega}
_n(\beta)\right\}g_n(\beta)} {\det\left\{\widehat{\Omega}_n(\beta)\right\}}.
\end{equation}
The determinant $\det\left\{\widehat{\Omega}_n(\beta)\right\}$ is a sum of
signed products of $k$ entries, so its degree is at most $2k$.
For $\beta\neq0$,
\[
\det\left\{\widehat{\Omega}_n(\beta)\right\}
=\beta^{2k}\det\left\{
\widehat{\Sigma}_{22,n}
-\frac{\widehat{\Sigma}_{12,n}+\widehat{\Sigma}_{21,n}}{\beta}
+\frac{\widehat{\Sigma}_{11,n}}{\beta^2}
\right\}.
\]
Thus,
\[
\lim_{\beta\to \pm \infty}
\frac{\det\left\{\widehat{\Omega}_n(\beta)\right\}}{\beta^{2k}}
=\det\left(\widehat{\Sigma}_{22,n}\right)>0.
\]
Since the determinant is a polynomial of degree at most $2k$,
this limit is its coefficient on $\beta^{2k}$. Therefore, its degree
is exactly $2k$.
Each entry of the adjugate is a signed determinant of a
$(k-1)\times(k-1)$ submatrix and has degree at most $2k-2$. Since $g_n(\beta)
$ is affine in $\beta$, the numerator in \eqref{eq:cue-objective-adjugate}
has degree at most $2k$.
By \eqref{eq:cu-gmm-end_point-value}, the endpoint value is strictly
positive when $Z^{\prime}y_2\neq0$, so the numerator and denominator must
both have degree $2k$. If $Z^{\prime}y_2=0$, the endpoint value is zero, and
the numerator has lower degree or is identically zero. Proposition \ref
{prop:ratio-formal} summarizes these findings.
\begin{proposition}
\label{prop:ratio-formal} Under Assumption \ref{ass:sigma-pd}, the CU-GMM
objective function in \eqref{eq:cue-objective-adjugate} can be written as
\begin{equation} \label{eq:ratio-formal}
Q_n(\beta)=\frac{p_n(\beta)}{q_n(\beta)},
\end{equation}
where $p_n(\beta)$ and $q_n(\beta)$ are the numerator and denominator in
\eqref{eq:cue-objective-adjugate}, respectively, and (i) $p_n(\beta)$ is a
polynomial of degree at most $2k$ and $q_n(\beta)$ is a polynomial of degree
exactly $2k$, (ii) $q_n(\beta)>0$ for every finite $\beta$, and (iii) both
polynomials have degree $2k$ if $Z^{\prime}y_2\neq0$.
\end{proposition}
The degrees in Proposition \ref{prop:ratio-formal} refer to the determinant
and adjugate representation. Common factors can reduce these degrees. Under
the homoskedastic variance structure, a factor of degree $2k-2$ cancels,
leaving a ratio of quadratic forms; see Section S-2 in
the online supplement.
\subsection{Global minimization through polynomial roots}
\label{sec:poly-approach} Under Assumption \ref{ass:sigma-pd}, the CU-GMM
objective $Q_n(\beta)$ is a continuously differentiable function of $\beta$
on $\mathbb{R}$. Using the representation $Q_n(\beta)=p_n(\beta)/q_n(\beta)$
from Section \ref{sec:poly-ratio-representation}, its derivative with
respect to $\beta$ is
\begin{equation}
\frac{\partial Q_n(\beta)}{\partial\beta} =\frac{\frac{\partial p_n(\beta)}{
\partial\beta}q_n(\beta) -p_n(\beta)\frac{\partial q_n(\beta)}{\partial\beta}
}{q_n(\beta)^2}.
\end{equation}
Since $q_n(\beta)>0$ for every finite $\beta$, the finite stationary points
of $Q_n(\beta)$ are exactly the solutions to
\begin{equation} \label{eq:foc}
h_n(\beta) =\frac{\partial p_n(\beta)}{\partial\beta}q_n(\beta) -p_n(\beta)
\frac{\partial q_n(\beta)}{\partial\beta} =0.
\end{equation}
Both $\frac{\partial p_n(\beta)}{\partial\beta}q_n(\beta)$ and $p_n(\beta)
\frac{\partial q_n(\beta)}{\partial\beta}$ are polynomials of degree at most
$4k-1$ in $\beta$. However, when $\deg\left(p_n\right)=\deg\left(q_n
\right)=2k$, the coefficients of the terms of degree $4k-1$ in their
difference cancel. If $p_n(.)$ has degree $d<2k$, both products have degree at
most $d+2k-1\leq4k-2$. If $p_n(.)$ is identically zero, then so is $h_n(\beta)$
. Therefore, $h_n(\beta)$ is either the zero polynomial or a polynomial of
degree at most $4k-2$.
If $h_n(\beta)$ is not identically zero, the stationary points of $Q_n(\beta)
$ are the real roots of a polynomial of degree at most $4k-2$. We know that
such a polynomial has at most $4k-2$ distinct real roots. Therefore, to find
the global minimum of $Q_n(\beta)$ over the extended real line $\overline{
\mathbb{R}}=\mathbb{R}\cup\{-\infty,+\infty\}$, we need to evaluate the
objective at no more than $4k-1$ points, that are at most $4k-2$ stationary
points and the endpoint $\beta=+\infty$. It is enough to evaluate one
endpoint because $Q_n(-\infty)=Q_n(+\infty)$.
If $h_n(\beta)$ is identically zero, the derivative of $Q_n(\beta)$ is zero
for every finite $\beta$. So, $Q_n(.)$ is constant on $\mathbb{R}$ and, by
continuity at the endpoints, on $\overline{\mathbb{R}}$. In this case,
evaluating the objective at any one point gives its global minimum. Theorem
\ref{thm:global-min-scalar} formalizes these findings.
\begin{theorem}
\label{thm:global-min-scalar} Suppose Assumption \ref{ass:sigma-pd} holds,
and let $h_n(\beta)$ be defined by \eqref{eq:foc}. If $h_n(\beta)$ is not
identically zero, let
\begin{equation*}
\mathcal{C}_n =\left\{\beta\in\mathbb{R}:h_n(\beta)=0\right\}\cup\left\{+
\infty\right\}.
\end{equation*}
Then $\mathcal{C}_n$ contains at most $4k-1$ points, and
\begin{equation} \label{eq:min-candidate-formal}
\min_{\beta\in\overline{\mathbb{R}}}Q_n(\beta) =\min_{\beta\in\mathcal{C}
_n}Q_n(\beta).
\end{equation}
If $h_n(\beta)$ is identically zero, then $Q_n(.)$ is constant on $\overline{
\mathbb{R}}$, and every point in $\overline{\mathbb{R}}$ is a global
minimizer.
\end{theorem}
\subsection{Computing polynomial roots using companion matrices}
\label{sec:companion}
Suppose $h_n(\beta)$ is not identically zero. Write
\begin{equation} \label{eq:foc-polynomial}
h_{n}(\beta)= \delta_{m,n} \beta^{m} + \delta_{m-1,n} \beta^{m-1}+\cdots+
\delta_{1,n} \beta + \delta_{0,n},
\end{equation}
where the coefficients $\delta_{i,n}$, $i = 0, 1, \ldots, m$, depend on the
observed data, $\delta_{m,n}\neq0$, and the degree $m$ is at most $4k-2$.
For polynomial equations of degree greater than four, Galois theory implies
that there is no general formula expressing their roots in terms of the
coefficients using arithmetic operations and radicals; see \citet[Section
14.7, Corollary 40]{DummitFoote2004}. Although particular polynomials may
admit such formulas, Proposition S.3 in the online
supplement gives a CU-GMM example with two instruments in which every global
minimizer is finite and cannot be expressed by radicals. This motivates a
numerical approach to computing the roots of $h_n(\beta)$. These roots can
be obtained as the eigenvalues of a companion matrix and computed using the
QR algorithm. To construct the companion matrix, consider the following
monic version of the polynomial $h_n(\beta)$:
\begin{equation*}
\overline{h}_{n}(\beta)= \beta^{m} + \overline{\delta}_{m-1,n}
\beta^{m-1}+\cdots+ \overline{\delta}_{1,n} \beta + \overline{\delta}_{0,n},
\end{equation*}
where, for $i = 0, 1, \ldots, m-1$,
\begin{equation}
\overline{\delta}_{i,n} = \frac{\delta_{i,n}}{\delta_{m,n}}.
\end{equation}
Define the companion matrix of $\overline{h}_{n}(\beta)$ as
\begin{equation} \label{eq:companion}
C_{\overline{h}}=
\begin{pmatrix}
0 & 1 & 0 & \cdots & 0 \\
0 & 0 & 1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1 \\
-\overline{\delta}_{0,n} & -\overline{\delta}_{1,n} & -\overline{\delta}
_{2,n} & \cdots & -\overline{\delta}_{m-1,n}
\end{pmatrix}
.
\end{equation}
Normalization does not change the roots, so the eigenvalues of $C_{\overline{
h}}$ are the roots of $h_n(\beta)$. Retaining the real eigenvalues gives the
following characterization of the finite stationary points.
\begin{proposition}
\label{prop:companion} Under Assumption \ref{ass:sigma-pd}, suppose $
h_n(\beta)$ is not identically zero. The finite stationary points of the
CU-GMM objective function $Q_{n}(\beta)$ given by \eqref{eq:cue-objective}
are exactly the real eigenvalues of the companion matrix $C_{\overline{h}}$
given by \eqref{eq:companion}.
\end{proposition}
Section S-2 in the online supplement gives the
generalized eigenvalue characterization under homoskedastic errors.
Proposition \ref{prop:companion} provides an eigenvalue approach for general
variance structures. Together with Theorem \ref{thm:global-min-scalar}, it
reduces global minimization to finding all real eigenvalues of the companion
matrix and comparing the objective values there and at infinity.
The roots of $h_n(\beta)$ can also be computed using other methods, such as
the Jenkins--Traub algorithm or Ehrlich--Aberth iterations. The important
point is not the specific root solver. It is that the computation requires
finding all real roots, rather than performing local optimization.
For implementation, the coefficients of $p_n(\beta)$ and $q_n(\beta)$ can be
computed symbolically from \eqref{eq:cue-objective-adjugate} or recovered by
polynomial interpolation as described in Section S-3 in the online supplement. Once the real roots
of $h_n(\beta)$ have been computed, evaluate the original CU-GMM objective
at these roots and at either infinite endpoint, as prescribed by Theorem \ref
{thm:global-min-scalar}.
\subsection{CU-GMM with singular variance matrices}
\label{sec:singular-univariate}
Throughout this subsection, we allow $\widehat{\Sigma}_n$ to be symmetric
positive semidefinite. As a result, $\widehat{\Omega}_n(\beta)$ is symmetric
positive semidefinite for every $\beta\in\mathbb{R}$, with rank between zero
and $k$. Its rank may vary with $\beta$. For every $\beta \in \mathbb{R} $,
define
\begin{equation} \label{eq:singular-objective}
Q_n(\beta)=g_n(\beta)^{\prime} \widehat{\Omega}_n(\beta)^{\dagger}g_n(\beta),
\end{equation}
where $\dagger$ denotes the Moore--Penrose inverse; see
\citet{Andrews87b} for conditions under which
test statistics using generalized inverse weighting matrices
keep their usual asymptotic chi-square distributions.\footnote{Alternatively, we can use the ridge-regularized variance estimator $\widehat{\Sigma}_n+\lambda_n I_{2k}$, where $\lambda_n>0$ does not depend on $\beta$. This matrix is positive definite, so the corresponding CU-GMM objective is a ratio of polynomials and the results of the preceding subsections apply. This regularization generally defines a different objective from \eqref{eq:singular-objective}.} For a symmetric positive semidefinite matrix, this inverse reciprocates its positive eigenvalues and leaves its zero eigenvalues equal to zero. It coincides with the ordinary inverse whenever the matrix is nonsingular.
Set $q_{0,n}(\beta)=1$. For $j=1,\ldots,k$, define
\begin{equation} \label{eq:singular-coefficients}
q_{j,n}(\beta)= \sum_{\substack{ I,J\subseteq\{1,\ldots,k\} \\ |I|=|J|=j}}
\det\left\{\left[\widehat{\Omega}_n(\beta)\right]_{I,J}\right\}^{2},
\end{equation}
where the subscripts $I,J$ select rows and columns of $\widehat{\Omega}
_n(\beta)$, respectively. Thus $q_{j,n}(\beta)$ is the sum of the squared $
j\times j$ minors. It is a polynomial of degree at most $4j$ and is positive
if and only if the rank is at least $j$. When the rank equals $r\geq1$,
Lemma \ref{lem:singular-inverse} gives
\begin{equation} \label{eq:singular-inverse}
\widehat{\Omega}_n(\beta)^{\dagger} =\frac{\displaystyle
\sum_{j=0}^{r-1}(-1)^j q_{r-1-j,n}(\beta) \widehat{\Omega}_n(\beta)^{2j+1}}{
q_{r,n}(\beta)}.
\end{equation}
\begin{lemma}
\label{lem:singular-inverse} For every $\beta\in\mathbb{R}$ at which $\operatorname{rank}
\left\{\widehat{\Omega}_n(\beta)\right\}=r\geq1$, $\widehat{\Omega}
_n(\beta)^{\dagger}$ is given by \eqref{eq:singular-inverse}, with $
q_{r,n}(\beta)>0$. If the rank is zero, $\widehat{\Omega}_n(\beta)^{
\dagger}=0$ and $Q_n(\beta)=0$.
\end{lemma}
Equation \eqref{eq:singular-inverse} supplies a rational formula at every
positive rank, including ranks below $k$. Let
\begin{equation*}
r_{\max}=\max_{\beta\in\mathbb{R}} \operatorname{rank}\left\{\widehat{\Omega}
_n(\beta)\right\}.
\end{equation*}
If $r_{\max}=0$, the objective is identically zero. Otherwise, define the
generic numerator and denominator by
\begin{align}
p_n(\beta)&=\sum_{j=0}^{r_{\max}-1}(-1)^j q_{r_{\max}-1-j,n}(\beta)
g_n(\beta)^{\prime} \widehat{\Omega}_n(\beta)^{2j+1}g_n(\beta),
\label{eq:singular-generic-p} \\
q_n(\beta)&=q_{r_{\max},n}(\beta) . \label{eq:singular-generic-q}
\end{align}
Both polynomials have degree at most $4r_{\max}$. These bounds need not be
sharp. When $r_{\max}=k$, the determinant and adjugate formula provides an
alternative representation with numerator and denominator degrees at most $2k
$.
The real zeros of $q_n(\beta)$ identify the rank-drop points. Away from
those points, the objective has the representation
\begin{equation} \label{eq:singular-generic-ratio}
Q_n(\beta)=\frac{p_n(\beta)}{q_n(\beta)} = \frac{\sum_{j=0}^{r_{
\max}-1}(-1)^j q_{r_{\max}-1-j,n}(\beta) g_n(\beta)^{\prime}\widehat{\Omega}
_n(\beta)^{2j+1}g_n(\beta) }{q_{r_{\max},n}(\beta)}
\end{equation}
The following proposition gives the corresponding partition.
\begin{proposition}
\label{prop:singular-intervals} Suppose $r_{\max}>0$. The polynomial $
q_n(\beta)$ given by \eqref{eq:singular-generic-q} is not identically zero.
Its distinct real roots $\beta_1<\cdots<\beta_m$, possibly none, are
exactly the points where the rank is below $r_{\max}$, with $
m\leq4r_{\max}\leq4k$. The open intervals
\begin{equation*}
\left(-\infty,\beta_1\right),\quad
\left(\beta_1,\beta_2\right),\quad\ldots,\quad
\left(\beta_m,+\infty\right),
\end{equation*}
together with the singleton sets $\{\beta_1\},\ldots,\{\beta_m\}$,
partition $\mathbb{R}$. On every interval, the rank is $r_{\max}$, and $
Q_n(\beta)$ has the representation \eqref{eq:singular-generic-ratio} and
is infinitely differentiable. If $m=0$, the rational representation holds
on all of $\mathbb{R}$.
\end{proposition}
At the rank-drop points, the objective values are computed directly from
\eqref{eq:singular-objective}. These points must be identified before
cancelling common factors in \eqref{eq:singular-generic-ratio}. Cancellation
can hide a rank-drop point whose actual objective value differs from the
simplified ratio. For example, if $\widehat{\Sigma}_n=\operatorname{diag}(0,1,1,0)$
and $g_n(\beta)=\left(\beta,0\right)^{\prime}$, then $\widehat{\Omega}
_n(\beta)=\operatorname{diag}(\beta^2,1)$. The generic objective is $
\beta^2/\beta^2=1$ for $\beta\neq0$, whereas $Q_n(0)=0$.
On all open intervals in Proposition \ref{prop:singular-intervals},
\begin{equation} \label{eq:singular-foc}
\frac{\partial}{\partial \beta}Q_n(\beta)=\frac{h_n(\beta)}{q_n(\beta)^2}
,\qquad h_n(\beta)= \frac{\partial p_n(\beta)}{\partial \beta} q_n(\beta) -
p_n(\beta) \frac{\partial q_n(\beta)}{\partial \beta} .
\end{equation}
Rank-drop points and infinite endpoints are evaluated separately.
We extend the objective to the infinite endpoints by defining
\begin{equation} \label{eq:singular-endpoints}
Q_n(-\infty)=\lim_{\beta\to-\infty}Q_n(\beta),\qquad
Q_n(+\infty)=\lim_{\beta\to+\infty}Q_n(\beta).
\end{equation}
Outside finitely many points where the rank of $\widehat{\Omega}_n(\beta)$
falls below its maximum, the objective equals the same nonnegative rational
function. Comparing its leading terms shows that the limits at positive and
negative infinity coincide, either at a finite nonnegative value or at $
+\infty$. When $\widehat{\Sigma}_{22,n}$ is positive definite, the common
limit is
\begin{equation} \label{eq:singular-endpoints-pd}
Q_n(-\infty)=Q_n(+\infty) =\frac{y_2^{\prime}Z\widehat{\Sigma}
_{22,n}^{-1}Z^{\prime}y_2}{n}.
\end{equation}
When $\widehat{\Sigma}_{22,n}$ is singular, replacing its inverse in this
expression by its Moore--Penrose inverse need not give the limit. If $
r_{\max}>0$, the endpoint values are obtained from the limits of the
rational expression in \eqref{eq:singular-generic-ratio}. If $r_{\max}=0$,
both endpoint values are zero.
The following theorem shows that although $Q_n(\beta)$ may be discontinuous
at rank-drop points, it remains lower semicontinuous. Its value at each
finite rank-drop point is therefore no larger than either one-sided limit.
Lower semicontinuity and compactness of the extended real line ensure that
the extended objective attains its minimum. Defining the endpoint values by
\eqref{eq:singular-endpoints} makes this minimum equal to $\inf_{\beta\in
\mathbb{R}}Q_n(\beta)$.
\begin{theorem}
\label{thm:singular-global} Suppose that $\widehat{\Sigma}_n$ is symmetric
positive semidefinite, and define $Q_n(\beta)$ on $\overline{\mathbb{R}}$
by \eqref{eq:singular-objective} and \eqref{eq:singular-endpoints}. Then
\begin{itemize}
\item[(i)] The extended objective is lower semicontinuous and attains a
finite minimum equal to $\inf_{\beta\in\mathbb{R}}Q_n(\beta)$.
\item[(ii)] If $r_{\max}>0$ and $h_n(\beta)$ is not identically zero, the
set
\begin{equation} \label{eq:singular-candidates}
\mathcal{C}_n= \left\{\beta\in\mathbb{R}:h_n(\beta)=0,\ q_n(\beta)>0\right\}
\cup\left\{\beta_1,\ldots,\beta_m\right\} \cup\left\{-\infty,+\infty\right\}
\end{equation}
is finite and contains every global minimizer. In particular,
\begin{equation} \label{eq:singular-global-min}
\inf_{\beta\in\mathbb{R}}Q_n(\beta) =\min_{\beta\in\overline{\mathbb{R}}
}Q_n(\beta) =\min_{\beta\in\mathcal{C}_n}Q_n(\beta).
\end{equation}
\item[(iii)] If $r_{\max}>0$ and $h_n(\beta)$ is identically zero, the
generic objective has the same constant value on all the open intervals.
Replacing the first set in \eqref{eq:singular-candidates} by any one finite
point with $q_n(\beta)>0$ still gives \eqref{eq:singular-global-min} and at
least one global minimizer.
\item[(iv)] If $r_{\max}=0$, every point of $\overline{\mathbb{R}}$ is a
minimizer and the minimum is zero.
\end{itemize}
\end{theorem}
\section{Linear IV models with multiple endogenous variables}
\label{sec:multivariate}
In the linear IV model with $q$ endogenous regressors, let $
Y_2=\left(y_2,\ldots,y_{q+1}\right)$. The exclusion restrictions associated
with $k>q$ instruments are
\begin{equation}
\mathbb{E}\left[Z^{\prime}\left(y_{1} - Y_{2} \beta^{\ast}\right)\right]=0,
\label{eq:general-population-moments}
\end{equation}
and the corresponding sample moments are
\begin{align}
g_n(\beta) &=n^{-1/2}Z^{\prime}\left(y_{1} -Y_{2} \beta \right)
=n^{-1/2}Z^{\prime} Y b(\beta), \label{eq:general-sample-moments} \\
&=\left(b(\beta)^{\prime}\otimes I_k\right) \operatorname{vec}\left(n^{-1/2}Z^{
\prime}Y\right), \notag
\end{align}
where $\beta\in\mathbb{R}^q$, $Y=\left(y_1,Y_{2}\right)$ and $
b(\beta)=\left(1,-\beta^{\prime}\right)^{\prime}$. As in Section \ref
{sec:setup}, let $\widehat{\Sigma}_n$ be a consistent estimator of the
asymptotic variance matrix of $\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$.
The corresponding estimator of the variance matrix of $g_n(\beta)$ is
\begin{equation}
\widehat{\Omega}_n(\beta) =\left(b(\beta)^{\prime}\otimes I_k\right)
\widehat{\Sigma}_n \left(b(\beta)\otimes I_k\right).
\label{eq:general-induced-variance}
\end{equation}
We impose Assumption \ref{ass:sigma-pd} on the enlarged $k(q+1)\times k(q+1)$
matrix $\widehat{\Sigma}_n$. This assumption implies that $\widehat{\Omega}
_n(\beta)$ is positive definite for every $\beta \in \mathbb{R}^{q}$.
Therefore, the CU-GMM objective is well defined by
\begin{equation}
Q_n(\beta)=g_n(\beta)^{\prime} \widehat{\Omega}_n(\beta)^{-1}g_n(\beta).
\label{eq:general-cue-objective}
\end{equation}
When $\widehat{\Sigma}_n$ is symmetric positive semidefinite, the CU-GMM
objective can also be defined using the Moore--Penrose inverse, as in
Section \ref{sec:singular-univariate}.
As in Section \ref{sec:setup}, the objective depends on $\beta$ only through
$b(\beta)$ and is homogeneous of degree zero in $b(\beta)$, $
Q_n(cb(\beta))=Q_n(b(\beta))$ for every constant $c\neq0$. Thus, defining
\begin{equation}
\overline{b}(\beta) =\frac{b(\beta)}{\left\| b(\beta)\right\|} =\frac{
\left(1,-\beta^{\prime}\right)^{\prime}} {\sqrt{1+\left\|\beta\right\|^2}},
\label{eq:general-normalized-direction}
\end{equation}
we have $Q_n(\beta)=Q_n\left(\overline{b}(\beta)\right)$.
Let $\widehat{\Sigma}_{Y_{2}Y_{2},n}$ denote the positive definite
lower-right $kq\times kq$ block of $\widehat{\Sigma}_n$, corresponding to
the estimated variance matrix of $\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y_2
\right)$.
For $q>1$, the limit at infinity can depend on the direction in which $\beta$
diverges. Let $d\in\mathbb{R}^q$ satisfy $\left\| d\right\|=1$. Along any
sequence with $\left\|\beta\right\|\to\infty$ and $\beta/\left\|\beta\right
\|\to d$, we have $\overline{b}(\beta)\to\left(0,-d^{\prime}\right)^{\prime}$
. The normalized moments and variance matrix consequently converge to $
-n^{-1/2}Z^{\prime}Y_2d$ and $\left(d^{\prime}\otimes I_k\right)\widehat{
\Sigma}_{Y_{2}Y_{2},n} \left(d\otimes I_k\right)$, respectively. The latter
matrix is positive definite because $\widehat{\Sigma}_{Y_{2}Y_{2},n}$ is
positive definite and $d\neq0$. Therefore, we define the objective at
infinity in direction $d$ by
\begin{equation}
\begin{split}
Q_{n}(\infty_{d}) &:=\lim_{\substack{ \left\|\beta\right\|\to\infty \\
\beta/\left\|\beta\right\|\to d}}Q_n(\beta) \\
&=\frac{1}{n}d^{\prime}Y_2^{\prime}Z \left[\left(d^{\prime}\otimes I_k\right)
\widehat{\Sigma}_{Y_{2}Y_{2},n} \left(d\otimes I_k\right)\right]^{-1}
Z^{\prime}Y_2d.
\end{split}
\label{eq:general-boundary-objective}
\end{equation}
Here, $\infty_{d}$ denotes a boundary point representing the direction $d$.
Opposite directions have the same objective value, that is $
Q_{n}(\infty_{d})=Q_{n}\left(\infty_{-d}\right)$.
To keep these directions, define the extended parameter space as
\begin{equation}
\overline{\mathbb{R}}^{q} =\mathbb{R}^q\cup \left\{\infty_{d}:d\in\mathbb{R}
^q,\ \left\| d\right\|=1\right\}. \label{eq:general-extended-space}
\end{equation}
To define the topology on this extended space, map each finite $\beta$ to $
\beta/\sqrt{1+\left\|\beta\right\|^2}$ in the open unit ball and each
direction at infinity $\infty_d$ to the boundary point $d$. With the
topology associated to this mapping, the extended space is homeomorphic to
the closed unit ball in $\mathbb{R}^q$ and therefore it is compact.
The normalized vector $\overline{b}(\beta)$ ranges over the part of the unit
sphere with a strictly positive first coordinate. Its limits at infinity
fill the boundary where that coordinate is zero. On the part of the unit
sphere with a nonnegative first coordinate, the associated variance matrix
is positive definite. The moments and variance matrix are continuous
functions of the normalized vector. Since matrix inversion is continuous on
the set of positive definite matrices, the objective extends continuously to
$\overline{\mathbb{R}}^{q}$. This extended space is compact, so the
objective attains its minimum. The following lemma summarizes this result.
\begin{lemma}
\label{lem:general-cue-existence} Under Assumption \ref{ass:sigma-pd} for
the enlarged matrix $\widehat{\Sigma}_n$, the objective in
\eqref{eq:general-cue-objective}, with the boundary values in
\eqref{eq:general-boundary-objective}, is continuous on the compact space $
\overline{\mathbb{R}}^{q}$ defined in \eqref{eq:general-extended-space}. It
attains its minimum, and
\begin{equation*}
\min_{\beta\in\overline{\mathbb{R}}^{q}}Q_n(\beta) =\inf_{\beta\in\mathbb{R}
^q}Q_n(\beta).
\end{equation*}
\end{lemma}
We define the CUE over this extended space by
\begin{equation}
Q_n\left(\widehat{\beta}_{\mathrm{CUE}}\right) =\min_{\beta\in\overline{
\mathbb{R}}^{q}}Q_n(\beta). \label{eq:general-cue-estimator}
\end{equation}
The minimum may occur at a finite parameter vector or at a direction at
infinity. Opposite boundary directions are distinct points of the extended
space but have the same objective value.
As in the case with one endogenous regressor, the probability that a global
minimizer occurs at a direction at infinity is zero if the conditional
distribution of $\operatorname{vec}\left(n^{-1/2}Z^{\prime}Y\right)$ given $\widehat{
\Sigma}_n$ is absolutely continuous with respect to Lebesgue measure on $
\mathbb{R}^{k(q+1)}$ almost surely; see Proposition S.4 in the online supplement.
\subsection{Polynomial ratio representation}
\label{subsec:general-polynomial-ratio}
We now follow Section \ref{sec:poly-ratio-representation} to write the
CU-GMM objective as a ratio of polynomials in $\beta_1,\ldots,\beta_q$.
Throughout this subsection, degree means total degree, that is the degree of
a monomial is the sum of its exponents, and the degree of a nonzero
polynomial is the largest degree among its monomials with nonzero
coefficients.
Partition $\widehat{\Sigma}_n$ into $k\times k$ blocks $\widehat{\Sigma}
_{ij,n}$, for $i,j=1,\ldots,q+1$, in the order of the columns of $
Y=\left(y_1,Y_2\right)$. Then
\begin{equation}
\begin{split}
\widehat{\Omega}_n(\beta) ={}&\widehat{\Sigma}_{11,n} -\sum_{j=1}^{q}\beta_j
\left(\widehat{\Sigma}_{1,j+1,n}+\widehat{\Sigma}_{j+1,1,n}\right) \\
&+\sum_{i=1}^{q}\sum_{j=1}^{q}\beta_i\beta_j \widehat{\Sigma}_{i+1,j+1,n}.
\end{split}
\label{eq:general-variance-expanded}
\end{equation}
So, each entry of $\widehat{\Omega}_n(\beta)$ is a polynomial of degree at
most two. Under Assumption \ref{ass:sigma-pd}, this matrix is positive
definite for every finite $\beta$. Therefore, the adjugate formula gives
\begin{equation}
Q_n(\beta) =\frac{g_n(\beta)^{\prime} \operatorname{adj}\left\{\widehat{\Omega}
_n(\beta)\right\}g_n(\beta)} {\det\left\{\widehat{\Omega}_n(\beta)\right\}}
\label{eq:general-polynomial-ratio}
\end{equation}
For every nonzero direction $d\in\mathbb{R}^q$, the coefficient on $\tau^2$
in $\widehat{\Omega}_n( \tau d)$ is the positive definite matrix $
\left(d^{\prime}\otimes I_k\right)\widehat{\Sigma}_{Y_{2}Y_{2},n}
\left(d\otimes I_k\right)$. Therefore, its determinant is the strictly
positive coefficient on $\tau^{2k}$ in $\det\left\{\widehat{\Omega}_n(\tau
d)\right\}$. Thus, $\det\left\{\widehat{\Omega}_n(\beta)\right\}$ has degree
exactly $2k$.
Each entry of the adjugate has total degree of at most $2(k-1)$. Since the
sample moments are linear in $\beta$, the numerator of $Q_n(\beta)$ has
total degree at most $2k$.
As in the scalar case, the values at infinity determine whether the
numerator reaches this degree. If $Z^{\prime}Y_2\neq0$, there is a unit
vector $d$ such that $Z^{\prime}Y_2d\neq0$, and
\eqref{eq:general-boundary-objective} gives $Q_n(\infty_{d})>0$. Along the
ray $\beta= \tau d$, the denominator has degree $2k$ in $\tau$, so the
numerator must also have degree $2k$ to give this positive limit. This
establishes the following extension of Proposition \ref{prop:ratio-formal}.
\begin{proposition}
\label{prop:general-ratio} Under Assumption \ref{ass:sigma-pd} for the
enlarged matrix $\widehat{\Sigma}_n$, the CU-GMM objective can be written as
\begin{equation}
Q_n(\beta)=\frac{p_n(\beta)}{q_n(\beta)},
\end{equation}
where $p_n(\beta)$ and $q_n(\beta)$ are the numerator and denominator in
\eqref{eq:general-polynomial-ratio}, respectively, and (i) $p_n(\beta)$ is a
polynomial of degree at most $2k$, including the possibility that it is
identically zero, and $q_n(\beta)$ is a polynomial of degree exactly $2k$,
(ii) $q_n(\beta)>0$ for every $\beta\in\mathbb{R}^q$, and (iii) both
polynomials have degree $2k$ if $Z^{\prime}Y_2\neq0$.
\end{proposition}
The degree bounds depend on the number of instruments $k$, but not on the
number of endogenous regressors $q$. Nevertheless, the number of polynomial
coefficients and the difficulty of solving the resulting systems can
increase with the number of endogenous regressors.
As in the scalar case, common factors can reduce the degrees in Proposition
\ref{prop:general-ratio}. The cancellation under homoskedastic errors
described in Section S-2 in the online supplement also
applies here, with $Y=(y_1,Y_2)$ and $b(\beta)=(1,-\beta^{\prime})^{\prime}$
, and reduces the objective to a ratio of quadratic forms.
\subsection{Global minimization through polynomial systems}
\label{subsec:general-polynomial-systems}
Under Assumption \ref{ass:sigma-pd} for the enlarged matrix $\widehat{\Sigma}
_n$, the CU-GMM objective is continuously differentiable on $\mathbb{R}^q$.
Using the representation in \eqref{eq:general-polynomial-ratio}, we have
\begin{equation}
\nabla Q_n(\beta) =\frac{h_n(\beta)}{q_n(\beta)^2}, \qquad
h_n(\beta)=q_n(\beta)\nabla p_n(\beta) -p_n(\beta)\nabla q_n(\beta).
\label{eq:general-foc}
\end{equation}
Since $q_n(\beta)>0$ for every finite $\beta$, the finite stationary points
are exactly the real solutions of the polynomial system
\begin{equation}
h_n(\beta)=0. \label{eq:general-stationary-system}
\end{equation}
Each element of the $q \times 1$ vector $h_{n}(\beta)$ has total degree of
at most $4k-1$. For $q>1$, the cancellation that gives the bound $4k-2$ in
Section \ref{sec:poly-approach} need not occur. When the scalar numerator
and denominator both have degree $2k$, their highest-degree terms are
multiples of the same monomial. In the multivariate case, their
highest-degree homogeneous parts need not be proportional, so their partial
derivatives need not give the same cancellation.
Unlike a nonzero polynomial in one variable, a multivariate polynomial
system can have infinitely many real solutions. Thus, even when the
objective is not constant, \eqref{eq:general-stationary-system} need not
give a finite set of stationary points. If all $h_{j,n}(\beta)$ are zero for
every $\beta \in \mathbb{R}^{q}$, the gradient is zero everywhere on $
\mathbb{R}^q$. The objective is then constant on $\mathbb{R}^q$ and, by
continuity, on the extended parameter space.
We must also consider minima at infinity. For $q>1$, the boundary contains a
sphere of directions, and its objective value can vary with the direction.
Let $p_{n,\infty}(d)$ and $q_{n, \infty}(d)$ be the homogeneous parts of
degree $2k$ of $p_n(\beta)$ and $q_n(\beta)$, evaluated at $d$.
Equivalently,
\begin{equation}
p_{n,\infty}(d)=\lim_{t\to\infty}\frac{p_n(td)}{t^{2k}}, \qquad
q_{n,\infty}(d)=\lim_{t\to\infty}\frac{q_n(td)}{t^{2k}}.
\label{eq:general-leading-polynomials}
\end{equation}
Here, $p_{n,\infty}(d)$ is identically zero if $p_n(\beta)$ has degree less
than $2k$. By Proposition \ref{prop:general-ratio} and its proof,
\begin{equation*}
q_{n,\infty}(d) =\det\left\{\left(d^{\prime}\otimes I_k\right) \widehat{
\Sigma}_{Y_{2}Y_{2},n}\left(d\otimes I_k\right)\right\}>0 \quad\text{for }
d\neq0.
\end{equation*}
Therefore, the boundary objective in \eqref{eq:general-boundary-objective}
can be written as
\begin{equation}
Q_{n}(\infty_{d})=\frac{p_{n,\infty}(d)}{q_{n,\infty}(d)}, \qquad \left\|
d\right\|=1. \label{eq:general-boundary-ratio}
\end{equation}
At every stationary point of the boundary objective on the unit sphere,
\begin{equation*}
\nabla\left\{\frac{p_{n,\infty}(d)}{q_{n,\infty}(d)}\right\} =\lambda d,
\end{equation*}
where $\lambda$ is a Lagrange multiplier. Since $d^{\prime}d=1$,
\begin{equation*}
d^{\prime}\nabla \left\{\frac{p_{n,\infty}(d)}{q_{n,\infty}(d)}\right\}
=\lambda.
\end{equation*}
For $d\neq0$, the ratio is homogeneous of degree zero, so for every $\tau>0$
,
\begin{equation*}
\frac{p_{n,\infty}(\tau d)}{q_{n,\infty}(\tau d)} =\frac{p_{n,\infty}(d)}{
q_{n,\infty}(d)}.
\end{equation*}
Therefore,
\begin{equation*}
\left. \frac{\partial}{\partial\tau} \left\{\frac{p_{n,\infty}(\tau d)}{
q_{n,\infty}(\tau d)}\right\} \right|_{\tau=1} = d^{\prime}\nabla \left\{
\frac{p_{n,\infty}(d)}{q_{n,\infty}(d)}\right\}=0.
\end{equation*}
Consequently, $\lambda=0$, and the gradient of the ratio with respect to $d$
is zero at every stationary point on the unit sphere. Thus, the stationary
boundary directions are exactly the real solutions of
\begin{equation}
h_{n,\infty}(d)=q_{n,\infty}(d)\nabla p_{n,\infty}(d) -p_{n,\infty}(d)\nabla
q_{n,\infty}(d) =0, \qquad \text{ and } \qquad d^{\prime }d =1.
\label{eq:general-boundary-system}
\end{equation}
Each element of the $q\times 1$ vector $h_{n,\infty}(d)$ is either
identically zero or a homogeneous polynomial of degree $4k-1$. The following
theorem extends Theorem \ref{thm:global-min-scalar}.
\begin{theorem}
\label{thm:general-global-min} Suppose Assumption \ref{ass:sigma-pd} holds
for the enlarged matrix $\widehat{\Sigma}_n$. Define
\begin{equation}
\mathcal{C}_n =\left\{\beta\in\mathbb{R}^q:h_n(\beta)=0\right\} \cup
\left\{\infty_{d}: \left\| d\right\|=1, \ h_{n,\infty}(d)=0\right\}.
\label{eq:general-candidate-set}
\end{equation}
Every global minimizer belongs to $\mathcal{C}_n$, and
\begin{equation}
\min_{\beta\in\overline{\mathbb{R}}^{q}}Q_n(\beta) =\min_{\beta\in\mathcal{C}
_n} Q_n(\beta) \label{eq:general-min-candidates}
\end{equation}
If $h_{n}(\beta)$ is zero for all $\beta$, then $Q_n(\beta)$ is constant on
the extended parameter space and every point is a global minimizer.
\end{theorem}
Thus, global minimization can be characterized by two polynomial systems:
one for finite stationary points and one for stationary directions at
infinity. Opposite boundary directions give the same objective value, so
only one representative of each pair is needed. For $q=1$, the boundary
ratio is constant on the two directions $d=\pm1$, and the finite system
reduces to the scalar polynomial equation in Section \ref{sec:poly-approach}.
Even when there are infinitely many stationary points, there are only
finitely many distinct stationary objective values. The real solution sets
of \eqref{eq:general-stationary-system} and
\eqref{eq:general-boundary-system} are defined by polynomial equations and
hence are semialgebraic. Therefore, they have finitely many connected
components, each semialgebraic; see \citet[Theorem
2.4.5]{BochnakCosteRoy1998}. Any two points in a component can be joined by
a continuous semialgebraic path in that component; see \citet[Proposition
2.5.13]{BochnakCosteRoy1998}. Moreover, the path is smooth on each open
subinterval of a finite partition of its parameter interval; see
\citet[Proposition 2.9.10 and its proof]{BochnakCosteRoy1998}. Along each
smooth piece of a path $\beta(\tau)$ in the finite stationary set,
\begin{equation}
\frac{\partial}{\partial \tau}Q_n\left(\beta(\tau)\right) =\nabla
Q_n\left(\beta(\tau)\right)^{\prime} \frac{\partial}{\partial \tau}
\beta(\tau)=0. \label{eq:general-constant-along-path}
\end{equation}
Thus, the objective is constant on each smooth piece. Continuity implies
that these constants agree across the partition points, so the objective is
constant on each connected component of the finite stationary set. The same
argument applies to the boundary ratio in \eqref{eq:general-boundary-ratio},
whose tangent derivatives are zero at every stationary direction. Thus, a
connected component can contain infinitely many stationary points but
contributes only one objective value.
\begin{proposition}
\label{prop:general-finite-critical-values} Under the assumptions of Theorem
\ref{thm:general-global-min}, the finite stationary set in
\eqref{eq:general-stationary-system} and the stationary direction set in
\eqref{eq:general-boundary-system} each have finitely many connected
components. The corresponding objective is constant on each component.
Consequently,
\begin{equation*}
\left\{Q_n(\beta):\beta\in\mathcal{C}_n\right\}
\end{equation*}
is a finite, nonempty set, and its smallest element is the global minimum of
the extended CU-GMM objective.
\end{proposition}
The proposition reduces the global minimization problem to a finite
comparison of objective values, although finding a representative from each
real connected component can still be challenging.
\subsection{Computing stationary objective values}
\label{subsec:general-solving-systems}
Proposition \ref{prop:general-finite-critical-values} shows that global
minimization requires comparing finitely many stationary objective values,
even when there are infinitely many stationary points. Therefore, we can
target these values directly rather than recover all stationary points. Both
the finite stationary system \eqref{eq:general-stationary-system} and the
boundary system \eqref{eq:general-boundary-system} must be considered.
Introduce a scalar $t$ representing the objective value and add
\begin{equation}
p_n(\beta)-t q_n(\beta)=0 \label{eq:general-critical-value-equation}
\end{equation}
to \eqref{eq:general-stationary-system}. These are $q+1$ polynomial
equations in $\beta$ and $t$. While $q_n(\beta)>0$ for every real $\beta$,
it can be equal to zero at complex points. Algebraic elimination considers
complex solutions, so we also introduce an auxiliary variable $s$ and impose
\begin{equation}
s q_n(\beta)-1=0. \label{eq:general-denominator-exclusion}
\end{equation}
This equation excludes denominator zeros and preserves every real stationary
point, with $s=1/q_n(\beta)$.
Substituting $p_n(\beta)=tq_n(\beta)$ into \eqref{eq:general-foc} gives, on
this level set,
\begin{equation*}
h_{n}(\beta) =q_n(\beta)\left( \nabla p_n(\beta) -t \nabla q_n(\beta)
\right).
\end{equation*}
Since \eqref{eq:general-denominator-exclusion} ensures $q_n(\beta)\neq0$, we
can replace the stationary equations by
\begin{equation}
\nabla p_n(\beta) -t \nabla q_n(\beta) = 0.
\label{eq:general-value-derivatives}
\end{equation}
We use \eqref{eq:general-critical-value-equation},
\eqref{eq:general-denominator-exclusion}, and
\eqref{eq:general-value-derivatives} for elimination. The derivative
equations have degree at most $2k-1$ in $\beta$ and are at most linear in $t$
, rather than degree at most $4k-1$ in $\beta$ for the original stationary
equations. Their degree in $(\beta,t)$ is at most $2k$.
A Gr\"{o}bner basis is a finite generating set for the ideal of polynomial
combinations of the equations that permits systematic polynomial division.
With a lexicographic order placing $\beta$ and $s$ before $t$, it permits
elimination of $\beta$ and $s$; see \citet[Theorem 1.4.1 and Corollary
1.4.2]{Hibi2013}. For this augmented stationary system, the resulting
elimination ideal contains a nonzero univariate polynomial
\begin{equation}
\psi_n(t)=0. \label{eq:general-value-eliminant}
\end{equation}
Every real stationary objective value is a root of $\psi_n(t)$. The
justification, given in Section S-5 in the online supplement, uses the fact that there
are finitely many complex stationary values as well, that is the rational
objective is constant on each connected component of the complex stationary
set where its denominator is nonzero. This additional argument is needed
because Proposition \ref{prop:general-finite-critical-values} concerns real
solutions. If the augmented system has no complex solution, its elimination
ideal contains $1$, indicating that there are no candidates.
For boundary values, apply the same construction to
\eqref{eq:general-boundary-system}, including $d^{\prime}d-1=0$, and add
\begin{equation}
p_{n,\infty}(d)-t q_{n,\infty}(d)=0, \qquad s q_{n,\infty}(d)-1=0.
\label{eq:general-boundary-value-equations}
\end{equation}
As in the finite case, replace the boundary stationary equations by
\begin{equation}
\nabla p_{n,\infty}(d) -t\nabla q_{n,\infty}(d)=0.
\label{eq:general-boundary-value-derivatives}
\end{equation}
These equations have degree at most $2k-1$ in $d$ and are at most linear in $
t$. Together with $d^{\prime}d=1$ and
\eqref{eq:general-boundary-value-equations}, they define the same augmented
solution set as the original boundary stationary equations. Eliminating $d$
and $s$ similarly gives a nonzero univariate polynomial whose roots include
every real stationary boundary value.
Collect the distinct real roots of the two elimination polynomials and order
them from smallest to largest. A real candidate value can arise from complex
stationary points without being attained at any real point. For each
candidate $t_i$, starting with the smallest, check for real solutions
according to its source. For a candidate from the finite stationary system,
check whether $p_n(\beta)-t_i q_n(\beta)=0$ admits a solution $\beta\in
\mathbb{R}^q$. For a candidate from the boundary stationary system, check
whether $p_{n,\infty}(d)-t_i q_{n,\infty}(d)=0$ admits a real solution
satisfying $d^{\prime}d=1$. No attained objective value can be below the
global minimum, and the global minimum is among the candidate values.
Therefore, the first candidate admitting a real solution is the global
minimum; the stationary equations need not be included in these checks.
Cylindrical algebraic decomposition provides an exact method for checking
whether these equations have real solutions; see \citet{Collins1975} or
\citet{BasuPollackRoy2006}. For a candidate value $t_i$ from the finite
stationary system, it partitions $\mathbb{R}^q$ into finitely many cells on
which $p_n(\beta)-t_i q_n(\beta)$ has constant sign and provides a sample
point in each cell. The equation has a real solution if and only if this
polynomial is equal to zero at one of these sample points. For a candidate
from the boundary stationary system, apply the same procedure to $
p_{n,\infty}(d)-t_i q_{n,\infty}(d)$ and $d^{\prime}d-1$, requiring both
polynomials to be equal to zero at the same sample point. At the smallest
candidate value for which either check yields a real solution, the
corresponding sample point provides a finite CUE or a minimizing boundary
direction. The computational cost of this procedure depends on the number of
variables and the structure of the polynomial systems.
\section{Numerical studies}\label{sec:simulation}
In this section, we first examine the finite-sample behavior of CUE
using the simulation designs in \citet{ChaoHasmanNeweysSansonWoutersen14}
and their extension to two endogenous regressors. We compare the
proposed polynomial method with quasi-Newton methods using various
starting points. The results show that failure by quasi-Newton
methods to find the global minimum can distort the reported
finite-sample behavior of CUE. We also compare CUE computed using
the polynomial method with 2SLS, LIML, two-step GMM, and iterated
GMM. In the reported designs, CUE generally has smaller absolute
median bias than these estimators and less dispersion than LIML.
Motivated by the theoretical results of \citet{NeweyWindmeijer09}
under many weak moment conditions and heteroskedasticity, we next
compare these estimators across five applications selected for
their prominence in the IV literature. Three are discussed in
the review article by \citet{AndrewsStockSun19}: the public
spending multiplier in \citet{AcconciaCorsettiSimonelli14},
returns to schooling in \citet{StephensYang14}, and the elasticity
of worker efficacy with respect to employment shares in
\citet{Young14}. We also consider two commonly referenced
applications: returns to schooling in \citet{AngristKrueger91}
and the intertemporal elasticity of substitution in \citet{Yogo04}.
Finally, we examine the CLR test using the simulation designs of
\citet{MoreiraRidderSharifvaghefi26}, which allow for heteroskedasticity
and autocorrelation. We compare implementations of the CLR test that
use the polynomial method or quasi-Newton methods to compute the
minimum of the CU-GMM objective.
The results show that the numerical method can change rejection
decisions under both the null and alternative hypotheses, with
potential consequences for the size and power of the CLR test.
\subsection{Simulation Studies: Estimation}
We closely follow the simulation designs of
\citet{ChaoHasmanNeweysSansonWoutersen14}, and consider
\begin{equation*}
y_{1i}=y_{2i}\beta^{\ast}+u_i
\qquad\text{and}\qquad
y_{2i}=z_{1i}\pi+v_{2i},
\end{equation*}
where $z_{1i}$ and $v_{2i}$ are independent standard normal
variables. The structural parameter is set to $\beta^{\ast} = 5$. Furthermore,
\begin{equation*}
u_i=\rho v_{2i}
+\sqrt{1-\rho^2}
\left(w\xi_{1i}+\sqrt{1-w^2}\xi_{2i}\right),
\qquad
\xi_{1i}\mid z_{1i}\sim\mathcal{N}(0,z_{1i}^2),
\qquad
\xi_{2i}\sim\mathcal{N}(0,1).
\end{equation*}
Conditional on $z_{1i}$, the variables $v_{2i}$, $\xi_{1i}$,
and $\xi_{2i}$ are independent, and observations are independent
across $i$.
We set the sample size to $n=800$. The parameter $\rho$ controls
endogeneity: when $\rho=0$, the regressor is exogenous, whereas
larger values of $|\rho|$ increase the absolute correlation between
$v_{2i}$ and $u_i$. We consider $\rho=0.3,0.6,0.9$.
The parameter $w$ controls heteroskedasticity. When $w=0$,
the structural error is homoskedastic, whereas larger values
of $w$ increase the variation in its conditional variance.
We consider $w=0,1/2,\sqrt{2}/2,\sqrt{3}/2,1$.
We control the first-stage signal through the concentration
parameter $\mu^2=n\pi^2$, taking values $1$, $4$, and $16$.
\footnote{The normalization uses $\operatorname{Var}(z_{1i})=\operatorname{Var}(v_{2i})=1$.
Under heteroskedasticity, $\mu^2$ alone does not fully characterize
the information available for estimating $\beta^{\ast}$.}
We consider $k=1,4,10,30,60$ instruments. When $k=1$,
$z_{1i}$ is the only instrument. As discussed in
Section \ref{sec:setup}, in this just-identified case, CUE,
LIML, 2SLS, two-step GMM, and iterated GMM coincide whenever
the sample first-stage coefficient is nonzero. We include
this case as a diagnostic benchmark for the numerical procedures.
For $k=4$, the instruments are the first four powers of $z_{1i}$.
For $k=10,30,60$, we also include interactions with independent
Bernoulli variables:
\[
Z_i=\left(
z_{1i},z_{1i}^2,z_{1i}^3,z_{1i}^4,
D_{1i}z_{1i},\ldots,D_{(k-4)i}z_{1i}
\right)^{\prime},
\qquad
D_{ji}\sim\operatorname{Bernoulli}\left(\frac{1}{2}\right).
\]
The dummy variables are independent of the other underlying
variables.
The reported results are based on 10,000 Monte Carlo replications.
In every replication, we treat the intercept as a deterministic
exogenous variable and partial it out by demeaning $y_{1i}$,
$y_{2i}$, and each component of $Z_i$. The tables report median
bias and, in parentheses, the 90\% interquantile range, defined
as the difference between the 95th and 5th percentiles of the
estimator's empirical distribution. To conserve space, the main
text reports results only for $k=1,10,60$, $\rho=0.6$, and
$w=\sqrt{2}/2$. The full set of simulation results is reported
in the online supplement.
\begin{table}[htbp]
\centering
\caption{Numerical computation of CUE:
median bias and 90\% interquantile range}
\label{tablelocalopt}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{ccccccc}
\toprule\toprule
& & Polynomial & \multicolumn{4}{c}{Quasi-Newton method with starting value:} \\
\cmidrule(lr){4-7}
$\mu^2$ & $k$ & method
& $\widehat{\beta}_{\mathrm{2SLS}}$ & $5$ & $4.5$ & $0$ \\
\midrule
\multirow{6}{*}{$1$}
& \multirow{2}{*}{$1$} & 0.2450 & 0.2450 & 0.0026 & -0.1525 & -0.8488 \\
& & (9.7000) & (9.7000) & (640.6557)
& (1064.2449) & (1993.0712) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & 0.3960 & 0.4015 & -0.2092 & -0.6690 & -3.8133 \\
& & (17.8233) & (18.9915) & (1334.1182)
& (1895.1619) & (4134.2534) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & 0.5016 & 0.5265 & -0.9370 & -2.0738 & -11.0427 \\
& & (25.2496) & (125.4843) & (1604.0900)
& (2311.2732) & (5924.6761) \\
\midrule
\multirow{6}{*}{$4$}
& \multirow{2}{*}{$1$} & 0.0286 & 0.0286 & -0.0004 & -0.0136 & -0.1201 \\
& & (3.3584) & (3.3584) & (4.0734)
& (4.8540) & (1556.6287) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & 0.1081 & 0.1180 & -0.1091 & -0.2737 & -1.0649 \\
& & (12.2227) & (10.7686) & (818.5090)
& (1359.1899) & (3945.2074) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & 0.1847 & 0.3174 & -0.7154 & -1.4356 & -5.1030 \\
& & (22.4477) & (48.8783) & (1340.0942)
& (2063.8933) & (5884.8722) \\
\midrule
\multirow{6}{*}{$16$}
& \multirow{2}{*}{$1$} & 0.0025 & 0.0025 & 0.0025 & 0.0025 & -0.0003 \\
& & (1.1681) & (1.1681) & (1.1681)
& (1.1681) & (1.1821) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & -0.0257 & -0.0211 & -0.0340 & -0.0510 & -0.1219 \\
& & (2.4633) & (2.2577) & (2.3959)
& (2.7748) & (2438.0664) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & -0.1694 & -0.0716 & -0.3475 & -0.5579 & -1.1984 \\
& & (12.0223) & (9.5672) & (375.6969)
& (1076.6775) & (5481.0188) \\
\bottomrule
\end{tabular}
\medskip
\begin{minipage}{\linewidth}
\footnotesize
\textit{Notes:}
Entries report median bias, with the difference between the
95th and 5th percentiles in parentheses.
Results are based on 10,000 replications with $n=800$,
$\beta^{\ast}=5$, $\rho=0.6$, and $w=\sqrt{2}/2$.
\end{minipage}
\end{table}
Table \ref{tablelocalopt} compares CUE computed using the polynomial
method with estimates obtained by the quasi-Newton method,
\footnote{We implement the quasi-Newton method in MATLAB using
\texttt{fminunc}.}
starting from the 2SLS estimator and the fixed values $5$, $4.5$,
and $0$. Starting at the true value $\beta^{\ast}=5$ provides an
infeasible benchmark for assessing whether knowledge of the true
parameter resolves the numerical difficulty. The other fixed
values assess sensitivity to starting points near and farther
from the true parameter.
Even in the just-identified case $k=1$, where an explicit solution
is available as discussed in Section \ref{sec:setup}, the
quasi-Newton procedure can fail to recover CUE. Starting at
the true parameter does not eliminate this problem. For example,
when $\mu^2=1$ and $k=1$, the 90\% interquantile range is
$640.66$ when starting at $5$, compared with $9.70$ for CUE
computed using the polynomial method. Starting at 2SLS recovers
CUE when $k=1$, but discrepancies arise when $k=10$ or $60$.
For $\mu^2=1$ and $k=60$, the corresponding ranges are $125.48$
and $25.25$.
These results show that numerical implementation can substantially
alter the reported finite-sample behavior of CUE.
\footnote{In the simulations reported in the online supplement,
the objective value obtained by the polynomial method is never
larger than that obtained by the quasi-Newton methods and can
be substantially smaller.}
The theoretical properties established for CUE, including those
under the many weak moment conditions of
\citet{NeweyWindmeijer09}, concern the estimator defined by
global minimization. They do not automatically extend to
numerical approximations that fail to recover that estimator.
\begin{table}[htbp]
\centering
\caption{Estimator comparison:
median bias and 90\% interquantile range}
\label{tableestimators}
\small
\setlength{\tabcolsep}{7pt}
\begin{tabular}{@{}ccccccc@{}}
\toprule\toprule
$\mu^2$ & $k$ & 2SLS & LIML & Two-step & Iterated & CUE \\
\midrule
\multirow{6}{*}{$1$}
& \multirow{2}{*}{$1$} & 0.2450 & 0.2450 & 0.2450 & 0.2450 & 0.2450 \\
& & (9.7000) & (9.7000) & (9.7000) & (9.7000) & (9.7000) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & 0.5451 & -0.0872 & 0.5373 & 0.5392 & 0.3960 \\
& & (1.3657) & (33.9866) & (1.4725) & (1.4799) & (17.8233) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & 0.5912 & -0.9389 & 0.5909 & 0.5907 & 0.5016 \\
& & (0.4693) & (63.1004) & (0.5636) & (0.5645) & (25.2496) \\
\midrule
\multirow{6}{*}{$4$}
& \multirow{2}{*}{$1$} & 0.0286 & 0.0286 & 0.0286 & 0.0286 & 0.0286 \\
& & (3.3584) & (3.3584) & (3.3584) & (3.3584) & (3.3584) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & 0.4211 & -0.5080 & 0.4077 & 0.4084 & 0.1081 \\
& & (1.1999) & (24.1870) & (1.2788) & (1.2792) & (12.2227) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & 0.5636 & -1.6394 & 0.5629 & 0.5624 & 0.1847 \\
& & (0.4605) & (57.7575) & (0.5564) & (0.5569) & (22.4477) \\
\midrule
\multirow{6}{*}{$16$}
& \multirow{2}{*}{$1$} & 0.0025 & 0.0025 & 0.0025 & 0.0025 & 0.0025 \\
& & (1.1681) & (1.1681) & (1.1681) & (1.1681) & (1.1681) \\
\addlinespace[3pt]
& \multirow{2}{*}{$10$} & 0.2163 & -0.3357 & 0.2043 & 0.2065 & -0.0257 \\
& & (0.8410) & (3.6623) & (0.8617) & (0.8602) & (2.4633) \\
\addlinespace[3pt]
& \multirow{2}{*}{$60$} & 0.4748 & -1.7711 & 0.4755 & 0.4759 & -0.1694 \\
& & (0.4329) & (30.9108) & (0.5170) & (0.5167) & (12.0223) \\
\bottomrule
\end{tabular}
\medskip
\begin{minipage}{\linewidth}
\footnotesize
\textit{Notes:}
Entries report median bias, with the difference between the
95th and 5th percentiles in parentheses.
Results are based on 10,000 replications with $n=800$,
$\beta^{\ast}=5$, $\rho=0.6$, and $w=\sqrt{2}/2$.
CUE is computed using the polynomial method.
\end{minipage}
\end{table}
Table \ref{tableestimators} compares CUE with 2SLS, LIML,
two-step GMM, and iterated GMM. For two-step GMM, we first
minimize the GMM objective using the identity weight matrix.
We then evaluate the estimated variance matrix at this preliminary
estimate and minimize the GMM objective using its inverse as
the fixed weight matrix. For iterated GMM, we repeat the
weight-update and minimization procedure 100 times.
In the reported overidentified designs, CUE has a larger
90\% interquantile range than 2SLS, two-step GMM, and iterated
GMM, but a smaller range than LIML. Except for one comparison
with LIML, CUE has smaller absolute median bias than all these
estimators. The combination of larger absolute median bias and
smaller dispersion shows that 2SLS, two-step GMM, and iterated
GMM can be tightly concentrated away from the true parameter.
For 2SLS, this pattern is consistent with the inconsistency
established under the many-instrument asymptotics of
\citet{Bekker94}. Two-step and iterated GMM display a similar
finite-sample pattern and produce results close to 2SLS;
iterating the weight updates does not generally recover CUE.
Table \ref{tab:multivariate_case} considers an extension of the
designs in \citet{ChaoHasmanNeweysSansonWoutersen14} to two
endogenous regressors:
\begin{align*}
y_{1i}
&=y_{2i}\beta_1^{\ast}+y_{3i}\beta_2^{\ast}+u_i,\\
y_{2i}&=z_{1i}\pi_1+v_{2i},
\qquad
y_{3i}=z_{2i}\pi_2+v_{3i},
\end{align*}
where $z_{1i}$, $z_{2i}$, $v_{2i}$, and $v_{3i}$ are independent
standard normal variables. The structural error is
\[
u_i=\rho_1v_{2i}+\rho_2v_{3i}
+\sqrt{1-\rho_1^2-\rho_2^2}
\left(w\xi_{1i}+\sqrt{1-w^2}\xi_{2i}\right),
\]
with
\[
\xi_{1i}\mid z_{1i},z_{2i}
\sim\mathcal{N}\left(0,\frac{z_{1i}^2+z_{2i}^2}{2}\right),
\qquad
\xi_{2i}\sim\mathcal{N}(0,1).
\]
Conditional on the instruments, $v_{2i}$, $v_{3i}$,
$\xi_{1i}$, and $\xi_{2i}$ are independent, and observations
are independent across $i$.
We set $n=800$, $\pi_1=\pi_2=\pi$ with $n\pi^2=4$,
$\beta_1^{\ast}=\beta_2^{\ast}=2$,
$\rho_1=\rho_2=0.6\sqrt{2}/2$, and $w=\sqrt{2}/2$.
The five instruments are $Z_i=\left(z_{1i},z_{1i}^2,z_{2i},z_{2i}^2,
z_{1i}z_{2i}\right)^{\prime}.$
The multivariate polynomial method requires solving systems
of polynomial equations, as discussed in
Section \ref{subsec:general-solving-systems}.
The results broadly agree with the univariate findings.
In the reported simulations, the objective value obtained
using the polynomial method is no larger than that obtained
using the quasi-Newton methods. Panel A of
Table \ref{tab:multivariate_case} shows that the numerical
method and starting point can substantially change the reported
median bias and dispersion of CUE. Starting from 2SLS produces
results close to those of the polynomial method in this design.
However, as the univariate results in Table \ref{tablelocalopt} show,
starting from 2SLS does not always recover CUE,
Panel B compares CUE with the other estimators. For both
coefficients, CUE has the smallest absolute median bias.
Its 90\% interquantile range is smaller than that of LIML,
but larger than those of 2SLS, two-step GMM, and iterated GMM.
\begin{table}[htbp]
\centering
\caption{Multivariate estimation:
median bias and 90\% interquantile range}
\label{tab:multivariate_case}
\small
\setlength{\tabcolsep}{6pt}
\begin{tabular}{@{}cccccc@{}}
\toprule \toprule
\multicolumn{6}{@{}l}{\textit{Panel A: Numerical methods for computing CUE}} \\
\addlinespace[4pt]
& \multicolumn{4}{c}{Quasi-Newton method with starting values:} & Polynomial \\
\cmidrule(lr){2-5}
Parameter & $\widehat{\beta}_{\mathrm{2SLS}}$
& $(2,2)^{\prime}$ & $(1.5,1.5)^{\prime}$ & $(0,0)^{\prime}$ & method \\
\midrule
$\beta_1$ & 0.0533 & -0.0351 & -0.1270 & -0.3624 & 0.0608 \\
& (12.4382) & (165.0510) & (74777.2079)
& (246259.2421) & (13.3977) \\
\addlinespace[4pt]
$\beta_2$ & 0.0552 & -0.0326 & -0.1134 & -0.3965 & 0.0527 \\
& (11.8000) & (61.8557) & (67203.5305)
& (269375.2046) & (12.5874) \\
\midrule
\multicolumn{6}{@{}l}{\textit{Panel B: Comparison across estimators}} \\
\addlinespace[4pt]
Parameter & 2SLS & Two-step & Iterated & LIML & CUE \\
\midrule
$\beta_1$ & 0.1951 & 0.1679 & 0.1737 & -0.1387 & 0.0608 \\
& (2.7174) & (2.6391) & (2.5934) & (23.0126) & (13.3977) \\
\addlinespace[4pt]
$\beta_2$ & 0.1930 & 0.1698 & 0.1694 & -0.1564 & 0.0527 \\
& (2.6814) & (2.6573) & (2.6254) & (22.1940) & (12.5874) \\
\bottomrule
\end{tabular}
\medskip
\begin{minipage}{\linewidth}
\footnotesize
\textit{Notes:}
Entries report median bias, with the difference between the
95th and 5th percentiles in parentheses.
The true coefficients are $\beta_1^{\ast}=\beta_2^{\ast}=2$.
In Panel A, $\widehat{\beta}_{\mathrm{2SLS}}$ denotes the
vector of 2SLS estimates used to initialize the quasi-Newton
method. In both panels, the final column reports CUE computed
using the polynomial method.
\end{minipage}
\end{table}
\subsection{Empirical applications}
\label{subsec:empirical-applications}
We compare CUE with 2SLS, LIML, and iterated GMM across the five
applications introduced above, retaining the specifications considered
in each application. Our comparison includes 20 specifications from
\citet{AcconciaCorsettiSimonelli14}, each with two instruments;
30 from \citet{StephensYang14}, each with three instruments; and
44 from \citet{Yogo04}, each with four instruments. We also consider
40 specifications from \citet{Young14}, with 21--59 instruments,
and 24 from \citet{AngristKrueger91}, with 3--180 instruments across
different census samples. The purpose is to assess how estimator
choice affects the reported coefficients, rather than to reassess
the validity of the instruments or the specifications.
Figure \ref{pctdiff} compares the percentage differences between CUE
and the other estimators across specifications. Each difference is
scaled by the average of the two estimates and expressed as a
percentage. We apply the real cube-root transformation to compress
large differences while preserving their signs and ordering.
The left panel plots the transformed differences for 2SLS against
those for LIML; the right panel compares 2SLS with iterated GMM.
Points near the origin indicate agreement with CUE, whereas points
near the 45-degree line indicate similar transformed percentage
differences for the two estimators. Because the measure divides by
the average of the estimates, large percentage differences can also
arise when that average is close to zero.
Blue diamonds denote \citet{AcconciaCorsettiSimonelli14}, red
upward triangles denote \citet{StephensYang14}, green squares denote
\citet{Young14}, orange circles denote \citet{Yogo04}, and magenta
downward triangles denote \citet{AngristKrueger91}.
\begin{figure}[tbp]
\centering
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[trim = 0mm 0mm 0mm 0mm, clip, width=\linewidth]{pctdiff_2.pdf}
\end{subfigure}
\begin{subfigure}[t]{0.49\textwidth}
\centering
\includegraphics[trim = 0mm 0mm 0mm 0mm, clip, width=\linewidth]{pctdiff_1.pdf}
\end{subfigure}
\caption{Cube-root-transformed percentage differences between CUE and alternative estimators}
\label{pctdiff}
\end{figure}
In the left panel, LIML estimates are close to CUE in many
specifications from \citet{Young14} and \citet{AngristKrueger91},
while 2SLS estimates show larger percentage differences. Other
specifications lie near the diagonal but away from the origin,
indicating similar percentage differences for 2SLS and LIML relative
to CUE. The right panel shows that iterated GMM often produces
estimates close to 2SLS. As in the simulations, repeated weight
updates need not recover CUE.
\begin{table}
\centering
\caption{Quartiles of absolute percentage differences between CUE and alternative estimators} \label{pdifftable}
\begin{tabular}{ccccc}
\hline\hline
Application & Estimator & 25th percentile & Median & 75th percentile \\
\hline
& 2SLS & 0.0800 & 0.3107 & 2.3135 \\
\cite{AcconciaCorsettiSimonelli14} & LIML & 0.0808 & 0.3482 & 2.1771 \\
& Iterated GMM & 0.0068 & 0.0282 & 0.1479 \\ \hline
& 2SLS & 5.1138 & 16.6242 & 41.5516 \\
\cite{StephensYang14} & LIML & 4.7907 & 13.5528 & 23.1492 \\
& Iterated GMM & 1.8198 & 9.8726 & 34.0177 \\ \hline
& 2SLS & 38.5661 & 84.4015 & 178.7024 \\
\cite{Young14} & LIML & 0.0003 & 0.0017 & 0.0163 \\
& Iterated GMM & 38.5604 & 84.4022 & 178.6905 \\ \hline
& 2SLS & 12.3428 & 36.3257 & 136.0933 \\
\cite{Yogo04} & LIML & 30.3433 & 98.6778 & 160.5497 \\
& Iterated GMM & 8.9277 & 84.7788 & 183.7710 \\ \hline
& 2SLS & 6.2673 & 18.6199 & 35.1441 \\
\cite{AngristKrueger91} & LIML & 0.5423 & 2.0802 & 3.1129 \\
& Iterated GMM & 4.4153 & 17.5725 & 34.6890 \\ \hline
\end{tabular}
\end{table}
Table \ref{pdifftable} reports the quartiles of the absolute
percentage differences within each application, before the cube-root
transformation. The median differences are below $0.4\%$ for all
three estimators in \citet{AcconciaCorsettiSimonelli14}. For
\citet{StephensYang14}, they range from $9.87\%$ to $16.62\%$,
and for \citet{Yogo04}, from $36.33\%$ to $98.68\%$.
Thus, substantial differences occur even in applications with few
instruments. In \citet{Young14}, the median difference is approximately
$84.4\%$ for both 2SLS and iterated GMM, compared with $0.0017\%$
for LIML. In \citet{AngristKrueger91}, the corresponding medians are
$18.62\%$, $17.57\%$, and $2.08\%$. These comparisons show that
estimator choice can materially affect the reported coefficients;
they do not by themselves identify which estimator is closer to the
true parameter or establish the source of the differences.
\subsection{Simulation studies: testing}\label{sec: MC-CLR}
We examine how the computation of the global minimum of the CU-GMM
objective affects the implementation of the CLR test. Following
\citet{MoreiraRidderSharifvaghefi23}, we simulate the standardized
statistic $R=\left(Z^{\prime}Z\right)^{-1/2}Z^{\prime}Y$ directly
from the Gaussian experiment
\begin{equation}\label{eq: dgp_r}
\operatorname{vec}(R)\sim\mathcal{N}\left(
\operatorname{vec}\left(\mu a^{\ast\prime}\right),\widetilde{\Sigma}\right),
\end{equation}
where $a^{\ast}=\left(\beta^{\ast},1\right)^{\prime}$,
$\mu\in\mathbb{R}^{k}$, and the variance matrix
\[
\widetilde{\Sigma}
=\left(I_2\otimes\left(Z^{\prime}Z/n\right)^{-1/2}\right)
\Sigma
\left(I_2\otimes\left(Z^{\prime}Z/n\right)^{-1/2}\right)
\]
is treated as known. Partition $\widetilde{\Sigma}$ into $k\times k$
blocks. We use the HAC design in \citet{MoreiraRidderSharifvaghefi23}, with
\[
\widetilde{\Sigma}_{11}=c_{11}I_k,\qquad
\widetilde{\Sigma}_{12}=\widetilde{\Sigma}_{21}=c_{12}J_k,\qquad
\widetilde{\Sigma}_{22}=c_{22}I_k,
\]
where $J_k$ has ones on its anti-diagonal and zeros elsewhere.
We set $c_{11}=1$, $c_{12}=100$, and
$c_{22}=c_{12}^{2}+c_{12}^{-3}$.
With $b(\beta)=\left(1,-\beta\right)^{\prime}$, the CU-GMM objective is
\begin{equation}\label{eq:CU-GMM_obj_r}
\begin{split}
Q_n(\beta)
&=\frac{1}{n}b(\beta)^{\prime}Y^{\prime}Z
\left[\left(b(\beta)^{\prime}\otimes I_k\right)\Sigma
\left(b(\beta)\otimes I_k\right)\right]^{-1}
Z^{\prime}Yb(\beta)\\
&=b(\beta)^{\prime}R^{\prime}
\left[\left(b(\beta)^{\prime}\otimes I_k\right)\widetilde{\Sigma}
\left(b(\beta)\otimes I_k\right)\right]^{-1}Rb(\beta).
\end{split}
\end{equation}
For testing $H_0:\beta^{\ast}=\beta_0$ against
$H_1:\beta^{\ast}\neq\beta_0$, the LR statistic is
\[
LR_n(\beta_0)=Q_n(\beta_0)
-\min_{\beta\in\overline{\mathbb{R}}}Q_n(\beta).
\]
Let $a_0=\left(\beta_0,1\right)^{\prime}$ and $b_0=b(\beta_0)$, and define
\begin{align*}
T&=\left[\left(a_0^{\prime}\otimes I_k\right)
\widetilde{\Sigma}^{-1}\left(a_0\otimes I_k\right)\right]^{-1/2}
\left(a_0^{\prime}\otimes I_k\right)\widetilde{\Sigma}^{-1}\operatorname{vec}(R),\\
S&=\left[\left(b_0^{\prime}\otimes I_k\right)
\widetilde{\Sigma}\left(b_0\otimes I_k\right)\right]^{-1/2}
\left(b_0^{\prime}\otimes I_k\right)\operatorname{vec}(R).
\end{align*}
The mapping between $R$ and $\left(S,T\right)$ is one-to-one.
Under $H_0$, $T$ is sufficient for $\mu$, while
$S\sim\mathcal{N}\left(0,I_k\right)$ is independent of $T$.
To approximate the conditional null distribution described in
Section \ref{sec:setup}, we hold $T=t$ fixed and generate 1,000 draws
of $S$. For each draw, we reconstruct $R$ from $\left(S,t\right)$
and compute $LR_n(\beta_0)$. The empirical 95th percentile gives
the conditional critical value for a test at the 5\% level.
We compare the LR statistics, conditional critical values, and rejection
decisions obtained using the polynomial method with those obtained using
quasi-Newton methods. We use $\beta_0$ as a starting point and also
consider additional starting points. As discussed by
\citet{MoreiraRidderSharifvaghefi20}, multiple starting points and a compact
parameterization can improve numerical performance, although they do not
guarantee global minimization. Specifically, setting $\beta=\tan\theta$
and $\overline{b}(\theta)=\left(\cos\theta,-\sin\theta\right)^{\prime}$ gives
\[
Q_n\left(\tan\theta\right)
=\overline{b}(\theta)^{\prime}R^{\prime}
\left[\left(\overline{b}(\theta)^{\prime}\otimes I_k\right)
\widetilde{\Sigma}
\left(\overline{b}(\theta)\otimes I_k\right)\right]^{-1}
R\overline{b}(\theta),
\qquad -\pi/2 \leq \theta \leq \pi/2.
\]
We denote the quasi-Newton implementations by
$\mathrm{QN}(n_d,n_r)$, where $n_d$ and $n_r$ are the numbers of
deterministic and random starting points, respectively.
Table \ref{tab:clr_test} reports results for $\beta_0=0$ and
$\mu=\sqrt{\lambda/k}\left(1,\ldots,1\right)^{\prime}$, so that
$\mu^{\prime}\mu=\lambda$. We vary $\beta^{\ast}$ through the
noncentrality parameter of the AR statistic,
\[
\Delta_{\mathrm{AR}}^{2}
=\left(\beta^{\ast}-\beta_0\right)^2\mu^{\prime}
\left[\left(b_0^{\prime}\otimes I_k\right)\widetilde{\Sigma}
\left(b_0\otimes I_k\right)\right]^{-1}\mu.
\]
We consider $\Delta_{\mathrm{AR}}=0$ under the null and
$\Delta_{\mathrm{AR}}=1.5$ under the alternative.
Failure to find the global minimum overstates the minimized objective
and therefore understates the LR statistic. That is, if a numerical
procedure returns $\widetilde{m}\geq
\min_{\beta\in\overline{\mathbb{R}}}Q_n(\beta)$, then
\[
Q_n(\beta_0)-\widetilde{m}
\leq Q_n(\beta_0)
-\min_{\beta\in\overline{\mathbb{R}}}Q_n(\beta).
\]
When the same conditional simulation draws are used, this ordering also
implies that the empirical critical value cannot exceed the one based
on global minimization. Since both the statistic and its critical value
can be understated, the ordering alone does not determine the direction
of changes in rejection decisions.
Table \ref{tab:clr_test} shows differences in the LR statistics in
39.0\% to 91.2\% of replications and differences in the conditional
critical values in every replication. Rejection decisions differ in
3.6\% to 8.5\% of replications under the null and 10.2\% to 20.1\%
under the alternative. These are frequencies of disagreement between
implementations, rather than rejection probabilities. Under the
alternative, the disagreement increases with the number of instruments
for each quasi-Newton implementation. Additional starting points reduce
the frequency of differences in the LR statistics, but do not eliminate
differences in rejection decisions.
\begin{table}[htbp]
\centering
\small
\setlength{\tabcolsep}{5pt}
\caption{Differences from the polynomial implementation of the CLR test}
\label{tab:clr_test}
\begin{threeparttable}
\begin{tabular}{@{}llrrrrrr@{}}
\toprule\toprule
& & \multicolumn{3}{c}{Null: $\Delta_{\mathrm{AR}}=0$}
& \multicolumn{3}{c}{Alternative: $\Delta_{\mathrm{AR}}=1.5$}\\
\cmidrule(lr){3-5}\cmidrule(l){6-8}
Method & Quantity & $k=2$ & $k=5$ & $k=10$ & $k=2$ & $k=5$ & $k=10$\\
\midrule
\multirow{3}{*}{$\mathrm{QN}(1,0)$}
& LR statistic & 75.6 & 78.9 & 85.3 & 82.0 & 88.8 & 91.2\\
& Critical value & 100.0 & 100.0 & 100.0 & 100.0 & 100.0 & 100.0\\
& Rejection decision & 6.6 & 8.5 & 7.8 & 10.4 & 14.2 & 20.1\\
\addlinespace[0.6em]
\multirow{3}{*}{$\mathrm{QN}(1,10)$}
& LR statistic & 49.7 & 57.7 & 67.9 & 42.5 & 53.1 & 61.4\\
& Critical value & 100.0 & 100.0 & 100.0 & 100.0 & 100.0 & 100.0\\
& Rejection decision & 4.5 & 6.5 & 8.4 & 11.1 & 14.9 & 19.7\\
\addlinespace[0.6em]
\multirow{3}{*}{$\mathrm{QN}(1,50)$}
& LR statistic & 45.7 & 52.5 & 61.9 & 39.0 & 50.6 & 59.5\\
& Critical value & 100.0 & 100.0 & 100.0 & 100.0 & 100.0 & 100.0\\
& Rejection decision & 3.6 & 4.7 & 7.2 & 10.2 & 13.0 & 16.7\\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item \textit{Notes:} Entries are percentages of replications in which
an implementation differs from the polynomial method. Critical values
are conditional 95th percentiles. $\mathrm{QN}(n_d,n_r)$ uses $n_d$
deterministic and $n_r$ random starting points.
\end{tablenotes}
\end{threeparttable}
\end{table}
\section{Conclusion}
\label{sec:conclusion}
This paper gives an algebraic characterization of global CU-GMM minimization
in linear IV models. For a scalar structural coefficient and a positive
definite reduced-form variance estimate, the minimum is obtained by
comparing the objective at real stationary roots with its limit at infinity.
Companion matrices express the root calculation as an eigenvalue problem,
providing a counterpart to the classical eigenvector characterization of
LIML. We extend the scalar characterization to the Moore--Penrose criterion,
allowing the variance matrix to be singular and its rank to change. With
several endogenous regressors, polynomial systems and directional limits at
infinity replace the scalar root calculation. Galois theory rules out a
general formula by radicals even with a scalar coefficient and two
instruments. Exact algebraic algorithms remain available.
These results suggest extensions to nonlinear GMM and minimum distance (MD)
estimation. Examples include GMM with polynomial moment restrictions and MD
that matches polynomial structural moments to sample moments. More general
models, including the Euler-equation models of \citet{HansenHeatonYaron96},
could be approximated by polynomial, rational, or spline functions. When the
resulting criterion is a ratio of polynomials on finitely many intervals,
the scalar approach suggests finding stationary points within each interval
and comparing their objective values with admissible boundary values and
one-sided limits, including limits at infinity. Changes of formula and
denominator zeros would be treated explicitly. Such approximations must
preserve the variance restrictions required by the criterion. Uniform
approximation of criterion values controls the error in the infimum;
approximating its minimizers also requires identification.
Under the variance specification considered here, these calculations provide
a way to compute the estimator studied by \citet{NeweyWindmeijer09} and the
minima needed for overidentification and likelihood ratio statistics.
Conditional likelihood ratio procedures, following \citet{Moreira03},
\citet{AndrewsMikusheva16} and \citet{MoreiraMoreira19}, also require
minimization of simulated criteria. Extending these procedures to nonlinear
models requires the appropriate distributional and variance estimation
assumptions. For approximated criteria, one must also control the effect of
approximation on critical values, as the analysis of discrete approximation
by \citet{MoreiraSharifvaghefi26} illustrates.
\bibliographystyle{ecta}
\bibliography{References_combined}