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 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.

4 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.

5 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.

6 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)

7 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

8 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	()