Code
print(inspect.signature(cm.GMM))(moment_fn, jacobian_fn=None, max_iterations=100, tolerance=1e-06, ridge=1e-08, fd_eps=1e-06)GMM
Callback-driven generalized method of moments
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 = (A+\rho\operatorname{diag}(A))^{-1}D'W\bar g,\qquad A=D'WD, \]
where \(\rho\) is dimensionless solver damping. The implementation normalizes by the square roots of \(A\)’s diagonal before solving and rejects numerical rank deficiency. A backtracking line search accepts only criterion decreases. Convergence requires both the undamped relative Newton step and a moment-scale-normalized first-order residual to meet tolerance; a small damped step or absolute criterion change is insufficient. 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^{-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_fnat \(\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.Each iteration evaluates finite moments and Jacobian with stable shapes. It computes the diagonally normalized damped direction and a separate undamped Newton direction. Convergence checks the latter relative to \(1+|\theta_j|\), together with the normalized first-order residual
\[ \frac{\{\sum_j (D'W\bar g)_j^2/A_{jj}\}^{1/2}} {\{n^{-1}\sum_i g_i'Wg_i\}^{1/2}}. \]
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.
The iteration counter increments after an accepted step. Both convergence checks are evaluated at the next iterate before checking the budget. Failure clears the old fit and does not install a partial solution.
Two-step weighting uses the undamped inverse of the first-step iid moment outer product. Damping does not change the statistical weighting matrix or covariance bread. Reported iterations add both stages.
Fitted moments and their Jacobian are owned snapshots. Summary and Wald inference use these arrays without calling Python again, reconstruct the undamped covariance bread, and form either \(A^{-1}/n\) or the sandwich. No arbitrary Python data deep copy is attempted.
The line search magnifies Python callback traffic. Numerical rank checks reject singular systems but do not establish strong statistical identification. Callback-based optimization is not detached from the GIL.
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. This option requires two-step iid weighting unless assume_optimal_weighting=True is explicitly supplied to summary() or wald_test(). The assertion does not repair a nonoptimal weight. The default 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\}\). Covariance inversion is undamped.
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. j_test_valid flags the intended two-step iid weighting regime, while j_test_assumptions lists additional unverified regularity assumptions. A raw j_stat under identity or dependent-data weighting is not automatically chi-square calibrated. Wald tests use the requested covariance.
6 Performance and numerical behavior
Each numerical-Jacobian evaluation calls Python at least \(2p+1\) times, and each line-search trial calls it again. Dense fitting uses \(m\times m\) weights and \(p\times p\) parameter systems; snapshots retain \(n\times m\) fitted moments. A callback Jacobian matters when moments are expensive. Damping is not a remedy for weak identification, and rank checks are not a statistical identification test.
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.
Code
cls = cm.GMM
display(HTML(html_table(["Public method"], public_methods(cls))))| Public method |
|---|
GMM(moment_fn, jacobian_fn=None, max_iterations=100, tolerance=1e-06, ridge=1e-08, fd_eps=1e-06) |
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, *, assume_optimal_weighting=False) |
wald_test(self, /, r, q=None, vcov='sandwich', omega='iid', lags=None, clusters=None, *, assume_optimal_weighting=False) |
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.73505286444329139 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.
Code
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 |
() |
termination_reason |
() |
inference_data |
() |
assume_optimal_weighting |
() |
j_test_valid |
() |
j_test_assumptions |
() |
weighting |
() |
vcov_type |
() |
omega_type |
() |
weight_matrix |
(3, 3) |
nobs |
() |
n_moments |
() |
original_n_moments |
() |
sketch_size |
() |
j_stat |
() |
j_df |
() |