Instrumental variables and GMM

From restrictions to an estimator

Suppose a structural residual u_i(\theta) satisfies E[z_i u_i(\theta_0)]=0 even though a regressor is correlated with the residual. Define the observation-level vector g_i(\theta)=z_i u_i(\theta) and minimize

Q_n(\theta)=\bar g_n(\theta)'W\bar g_n(\theta),\qquad \bar g_n(\theta)=\frac1n\sum_i g_i(\theta).

The weighting matrix decides how competing sample restrictions are reconciled. Identity weighting is simple; efficient two-step GMM estimates a moment covariance at a preliminary fit and then uses its inverse. Neither weighting choice can repair invalid instruments or weak identification.

Trex takes a callback moment_cond(z, y, x, theta) returning one row per observation, shape (n, q). Do not return only the average: observation-level moments are needed for covariance estimation.

Endogenous regressor example

import numpy as np
import pandas as pd
import torch
from trex import LinearRegression
from trex.gmm import GMMEstimator

torch.set_num_threads(2)
rng = np.random.default_rng(807)
n = 1600
z = rng.normal(size=(n, 2))
v, e = rng.normal(size=(2, n))
x = .9*z[:, 0] + .5*z[:, 1] + v
y = .4 + 1.6*x + .8*v + e
X = np.column_stack([np.ones(n), x])
Z = np.column_stack([np.ones(n), z])

def iv_moments(z, y, x, theta):
    return z * (y - x @ theta)[:, None]

fits = {backend: GMMEstimator(iv_moments, backend=backend,
            weighting_matrix="optimal", device="cpu").fit(Z, y, X)
        for backend in ["scipy", "torch"]}
ols = LinearRegression(device="cpu").fit(X, y)
pd.DataFrame({"truth": [.4, 1.6], "OLS": ols.params["coef"].numpy(),
              "GMM / SciPy": fits["scipy"].theta_,
              "GMM / Torch": fits["torch"].theta_,
              "GMM SE": fits["torch"].std_errors_}, index=["intercept", "slope"]).round(4)

	truth	OLS	GMM / SciPy	GMM / Torch	GMM SE
intercept	0.4	0.3864	0.4133	0.4133	0.0325
slope	1.6	2.0045	1.6049	1.6049	0.0322

OLS absorbs the correlation induced by v; the instruments exclude that component. Both backends target the same sample criterion. SciPy uses numerical derivatives, while the Torch backend differentiates the actual callback. For Torch callbacks, avoid converting intermediate quantities to NumPy or detaching them from the parameter graph.

Weighting, dependence and inference

With G=\partial\bar g_n(\theta)/\partial\theta' and a consistent long-run covariance \widehat\Omega, the reported covariance is

\widehat V=\frac1n(G'WG)^{-1}G'W\widehat\Omega WG(G'WG)^{-1}.

Under independent sampling, Trex uses the uncentered second product n^{-1}\sum_i g_i g_i' (HC0 convention). With ordered, weakly dependent observations, iid=False uses Bartlett/Newey–West accumulation for both optimal weights and final inference:

\widehat\Omega=\widehat\Gamma_0+ \sum_{\ell=1}^{L}\left(1-\frac{\ell}{L+1}\right) (\widehat\Gamma_\ell+\widehat\Gamma_\ell'),\quad \widehat\Gamma_\ell=\frac1n\sum_{i=\ell+1}^n g_i g_{i-\ell}'.

Set max_lags explicitly when a particular bandwidth is part of the design. Reordering the rows changes this calculation. HAC is not a substitute for an arbitrary cluster structure.

m = fits["torch"]
print({"moment_mean": iv_moments(Z, y, X, m.theta_).mean(0).round(5).tolist(),
       "parameter_covariance_shape": list(m.covariance_.shape),
       "backend_coefficient_max_gap": float(np.max(np.abs(m.theta_-fits["scipy"].theta_)))})
assert np.allclose(m.covariance_, m.vtheta_/n)

{'moment_mean': [6e-05, -0.00674, 0.01349], 'parameter_covariance_shape': [2, 2], 'backend_coefficient_max_gap': 3.4866181053416767e-09}

Use covariance_ for estimated coefficient variance. vtheta_ retains the asymptotic, non-divided scale. A fixed positive-definite matrix is accepted as weighting_matrix; "identity" avoids the second weighting step. init_params determines parameter dimension when it differs from x.shape[1], including nonlinear moment models.

Rank-deficient moment covariances and locally unidentified parameters raise errors. Near-singularity can still produce unstable estimates before a rank check fails: examine instrument relevance, perturbations and uncertainty rather than taking optimizer success as identification.

API: GMMEstimator, Torch backend, SciPy backend. Background: Hansen (1982) and Newey–West (1987), linked in references.

Experimental scalar CUE: polynomial candidates

Moreira, Newey and Sharifvaghefi (2026) show that linear-IV CUE with one endogenous regressor and a fixed positive-definite joint reduced-form covariance is a polynomial ratio. Writing g(\beta)=a-\beta b and \Omega(\beta)=S_{11}-\beta(S_{12}+S_{21})+\beta^2 S_{22},

Q(\beta)=g(\beta)'\Omega(\beta)^{-1}g(\beta) =\frac{p(\beta)}{q(\beta)},\qquad h(\beta)=p'(\beta)q(\beta)-p(\beta)q'(\beta).

With k instruments, the degrees of p,q are at most 2k and the stationary polynomial has degree at most 4k-2. Evaluate the original matrix criterion at all real stationary roots and at the common infinite endpoint b'S_{22}^{-1}b. This targets CUE, not the fixed-weight two-step GMM estimator above.

from trex.gmm import ScalarIVCUE

# One endogenous column; demeaning handles the intercept.
cue_problem = ScalarIVCUE.from_data(z, y, x, weighting="reduced_form")
cue_result = cue_problem.minimize_polynomial()
print({"beta": cue_result.beta, "objective": cue_result.fun,
       "at_infinity": cue_result.at_infinity,
       "certified": cue_result.certified})

{'beta': 1.6047524911443598, 'objective': 0.1998541053005381, 'at_infinity': False, 'certified': False}

This opt-in CPU/float64 prototype uses Chebyshev interpolation and companion eigenvalues in two projective coordinate charts, then checks/polishes roots against derivatives of the matrix criterion. It supports HC0, Bartlett HAC and one-way CR0 cluster-score covariances. weighting="raw" instead gives structural-residual CUE; with HC0 this is the existing quadratic GEL profile up to its factor 1/(2n). Do not compare optimizers using different weight definitions and attribute the difference to optimization.

The polynomial solver currently permits at most 12 instruments and raises on failed interpolation or stationarity checks. Numerical checks are not an exact global certificate. An infinite result is reported as such, not turned into a large finite estimate. Numerically constant objectives are flagged and do not establish identification. Singular covariances, multiple endogenous regressors, CLR calculations and weak-IV standard errors are outside this implementation; existing GMM/GEL defaults are unchanged.

See the executed comparison notebook and implementation and reproduction notes.

API: ScalarIVCUE, PolynomialCUEResult.