GMM

Callback-driven generalized method of moments

from _api_doc_utils import *

1 Where it fits

Group: Estimation interfaces

GMM solves moment restrictions of the form

\[ \mathbb E[g_i(\theta)] = 0. \]

The user supplies a Python callback returning the per-observation moment matrix. In exactly identified cases the class can solve by Gauss-Newton; in overidentified cases it can use identity or two-step weighting and report sandwich covariance estimates.

2 Criterion and weighting

The moment callback returns rows \(g_i(\theta)'\in\mathbb R^m\), and

\[ \bar g(\theta)=\frac1n\sum_i g_i(\theta). \]

For a fixed positive-semidefinite weight matrix \(W\), the implemented criterion is

\[ Q(\theta;W) = \frac12\bar g(\theta)'W\bar g(\theta). \]

With mean-moment Jacobian \(D=\partial\bar g(\theta)/\partial\theta'\), one Gauss-Newton direction is

\[ \Delta = (D'WD+\rho I)^{-1}D'W\bar g, \]

where \(\rho\) is the configured numerical ridge. A backtracking line search updates \(\theta\leftarrow\theta-a\Delta\) only when the criterion decreases. Convergence is based on the step norm or criterion change; exhausting the iteration budget raises.

Identity weighting uses \(W=I_m\). Automatic weighting uses identity when \(m=p\) and otherwise performs two steps: it first obtains \(\tilde\theta\) under identity weighting, computes

\[ \hat\Omega_1=\frac1n\sum_i g_i(\tilde\theta)g_i(\tilde\theta)', \qquad W=(\hat\Omega_1+\rho I)^{-1}, \]

and re-optimizes from \(\tilde\theta\). This is iid two-step GMM; the fit-stage weight cannot be HAC or clustered.

If no Jacobian callback is supplied, \(D\) is computed by central differences with step \(h_j=\text{fd eps}\times\max\{|\theta_j|,1\}\). A supplied callback must return the Jacobian of the mean moments, with shape \(m\times p\), not per-observation derivatives.

3 Implementation walkthrough

The callback bridge, numerical differentiation, Gauss-Newton loop, two-step weighting, sketch, and covariance assembly are package-owned.

The first callback evaluation at \(\theta_0\) establishes \(n\) and \(m\) and must return a nonempty two-dimensional float array. The fit rejects \(m<p\). Callback outputs are copied from NumPy into owned Rust arrays on every evaluation.
Without a Jacobian callback, the solver first evaluates the base moment shape, then calls moment_fn at \(\theta\pm h_je_j\) for every parameter and central-differences the column means. With a callback, the returned matrix is used directly and checked against \((m,p)\) after any sketch projection.
One Gauss-Newton iteration evaluates moments and Jacobian, forms \(D'WD\), adds ridge to its diagonal, explicitly inverts it, and computes \(\Delta=(D'WD+\rho I)^{-1}D'W\bar g\). A step norm below tolerance returns before updating.
Otherwise the solver tries \(\theta-\alpha\Delta\) for \(\alpha=1,1/2,1/4,\ldots\) down to \(10^{-8}\). Every trial calls the Python moment function again. Only a strict criterion decrease is accepted; failure to find one raises rather than returning the current iterate.
An accepted criterion improvement below tolerance counts as convergence. The iteration counter increments only after an accepted step. Reaching max_iterations raises and no partial fit is installed.
Identity weighting runs this loop once. Two-step weighting computes the uncentered iid outer product of first-step moment rows, adds the numerical ridge before inversion, and starts a second complete Gauss-Newton solve at the first-step estimate. Reported iterations add both stages and first_step_theta is retained internally.
Summary reevaluates moments and the Jacobian at the stored estimate. It reconstructs \(A=D'WD\), applies the same diagonal ridge to its inverse, builds the requested moment covariance, and then forms either \(A^{-1}/n\) or the full sandwich. Wald tests duplicate this covariance path rather than caching it.

The line search makes the local quadratic method robust to a poor full step, but it also magnifies Python callback traffic. The diagonal ridge is a numerical device in fitting, weight inversion, and reported covariance; a successful inverse is not evidence of strong identification.

4 Moment sketch

The sketch does not reduce the number of observations. It draws a dense Rademacher projection \(R\in\mathbb R^{m\times s}\) with entries \(\pm1/\sqrt{s}\) and replaces

\[ g_i(\theta)\ \text{by}\ R'g_i(\theta), \qquad D\ \text{by}\ R'D. \]

The sketch size must satisfy \(p\leq s\leq m\). The Python callback still constructs the full \(n\times m\) moment matrix on every evaluation before projection, so sketching reduces weight-matrix and Jacobian linear algebra but not callback cost or full-moment allocation. Inference and the overidentification statistic apply to the projected moments.

5 Inference

Let \(A=D'WD\) and let \(\hat\Omega\) be the selected covariance of \(g_i(\hat\theta)\). The vanilla option returns

\[ \widehat V_{\mathrm{vanilla}}=\frac1nA^{-1}, \]

which assumes the fitted \(W\) is the inverse moment covariance. The sandwich option returns

\[ \widehat V_{\mathrm{sandwich}} = \frac1n A^{-1}D'W\hat\Omega WD A^{-1}. \]

Available \(\hat\Omega\) estimators are uncentered iid moment outer products, Bartlett Newey-West autocovariances, and cluster sums. The cluster version has no finite-cluster correction. Default Newey-West lags are \(\max\{1,\lfloor4(n/100)^{2/9}\rfloor\}\). The numerical ridge is also used in covariance inversions.

For \(m>p\), the summary reports

\[ J=n\,\bar g(\hat\theta)'W\bar g(\hat\theta) \]

with \(m-p\) degrees of freedom, but does not compute a \(p\)-value. Wald tests use the requested covariance.

6 Performance and numerical behavior

Each numerical-Jacobian evaluation calls the Python moment function at least \(2p+1\) times, and each line-search trial calls it again. With dense moments, weight construction and fitting use \(m\times m\) matrices and parameter systems of size \(p\). Two-step fitting roughly doubles optimization work. A callback Jacobian is therefore important when moments are expensive. Weak identification appears as an ill-conditioned \(D'WD\); the ridge may permit a numerical answer without establishing statistical identification.

7 Python API

Constructor: cm.GMM

Construct with GMM(moment_fn, jacobian_fn=None, max_iterations=100, tolerance=1e-6, ridge=1e-8, fd_eps=1e-6). fit(data, theta0, weighting='auto') stores the fitted parameters and raises when the iteration budget is exhausted without meeting the convergence tolerance. fit_sketch(...) projects the moment columns to a smaller dimension. summary(vcov='sandwich', omega='iid', lags=None, clusters=None) controls inference.

print(inspect.signature(cm.GMM))

(moment_fn, jacobian_fn=None, max_iterations=100, tolerance=1e-06, ridge=1e-08, fd_eps=1e-06)

cls = cm.GMM
display(HTML(html_table(["Public method"], public_methods(cls))))

Public method
`fit(self, /, data, theta0, weighting='auto')`
`fit_sketch(self, /, data, theta0, sketch_size, weighting='auto', seed=None)`
`summary(self, /, vcov='sandwich', omega='iid', lags=None, clusters=None)`
`wald_test(self, /, r, q=None, vcov='sandwich', omega='iid', lags=None, clusters=None)`

8 Minimal example

def moments(theta, data):
    resid = data['y'] - data['x'] * theta[0]
    return data['z'] * resid[:, None]
def jac(theta, data):
    return -(data['z'].T @ data['x'][:, None]) / data['x'].shape[0]
rng = np.random.default_rng(20)
n = 300
z = rng.normal(size=(n, 3))
v = rng.normal(size=n)
x = z @ np.array([0.9, 0.4, -0.3]) + v
y = 1.2 * x + 0.5 * v + rng.normal(size=n) * 0.3
model = cm.GMM(moments, jacobian_fn=jac, max_iterations=200)
model.fit({'x': x, 'y': y, 'z': z}, np.array([0.0]), weighting='identity')
print(model.summary()['coef'])
print(model.summary()['j_stat'])

[1.16176096]
0.7350528644433015

9 `summary()` contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

def moments(theta, data):
    resid = data['y'] - data['x'] * theta[0]
    return data['z'] * resid[:, None]
def jac(theta, data):
    return -(data['z'].T @ data['x'][:, None]) / data['x'].shape[0]
rng = np.random.default_rng(120)
n = 120
z = rng.normal(size=(n, 3))
v = rng.normal(size=n)
x = z @ np.array([0.9, 0.4, -0.3]) + v
y = 1.2 * x + 0.5 * v + rng.normal(size=n) * 0.3
model = cm.GMM(moments, jacobian_fn=jac, max_iterations=200)
model.fit({'x': x, 'y': y, 'z': z}, np.array([0.0]), weighting='identity')
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))

summary() key	shape
`coef`	`(1,)`
`se`	`(1,)`
`vcov`	`(1, 1)`
`criterion`	`()`
`nit`	`()`
`converged`	`()`
`weighting`	`()`
`vcov_type`	`()`
`omega_type`	`()`
`weight_matrix`	`(3, 3)`
`nobs`	`()`
`n_moments`	`()`
`original_n_moments`	`()`
`sketch_size`	`()`
`j_stat`	`()`
`j_df`	`()`