statspai.gmm

gmm

Dynamic panel GMM estimators for StatsPAI.

Provides: - Arellano-Bond (1991) first-differenced GMM - Blundell-Bond (1998) system GMM - Arellano-Bond test for serial correlation (AR(1)/AR(2)) - Hansen/Sargan test for overidentifying restrictions

xtabond

xtabond(data: DataFrame, y: str, x: Optional[Sequence[str]] = None, id: str = 'id', time: str = 'time', lags: int = 1, gmm_lags: Tuple[int, Optional[int]] = (2, None), method: str = 'difference', twostep: bool = False, robust: bool = True, alpha: float = 0.05, predetermined: Optional[Sequence[str]] = None, endogenous: Optional[Sequence[str]] = None, predetermined_lags: Optional[Tuple[Optional[int], Optional[int]]] = None, endogenous_lags: Optional[Tuple[Optional[int], Optional[int]]] = None, collapse: bool = False, time_dummies: bool = False, orthogonal: bool = False, cluster: Optional[str] = None, steps: Optional[object] = None, iter_tol: float = 1e-10, iter_maxiter: int = 100, ah_instrument: str = 'levels', constant: Optional[bool] = None) -> CausalResult

Arellano-Bond / Blundell-Bond dynamic panel GMM estimator.

Equivalent to Stata's xtabond / xtabond2.

Parameters:

Name	Type	Description	Default
data	DataFrame	Balanced or unbalanced panel in long format.	required
y	str	Dependent variable.	required
x	list of str	Strictly exogenous regressors. Entered in first differences both as regressors and as their own (standard) instruments. Accepts Stata lag-operator syntax — "w", "L.w", "L2.k", "L(0/2).k" — so the canonical Arellano-Bond (1991) employment equation is a single call: x=["l(0/1).w", "l(0/2).k"], lags=2.	None
id	str	Unit identifier.	'id'
time	str	Time period variable. Treated as an ordinal sequence: the sorted distinct values define consecutive periods, so a missing period (gap) is recognised, but non-integer / irregularly-spaced codes are collapsed to their rank order.	'time'
lags	int	Number of lags of Y to include (ρ₁ Y_{t-1} + ... + ρ_p Y_{t-p}).	1
gmm_lags	tuple(min, max)	Range of lags of Y (in levels) used as GMM instruments. min must be ≥ 2 (deeper lags are orthogonal to the differenced error). max=None uses all available deeper lags, matching Stata's xtabond default. Setting max caps the instrument count (Stata's maxldep()).	(2, None)
method	str	'difference' — Arellano-Bond first-differenced GMM; parity with Stata's xtabond ..., noconstant. 'system' — Blundell-Bond (1998) system GMM: the level equation is stacked alongside the transformed one and instrumented with lagged differences (E[Δy_{i,t-1}(α_i + ε_{it})] = 0), and an intercept becomes identified. Prefer it when the series is persistent: as ρ approaches 1 the lagged levels become weak instruments for the differences and difference GMM is badly biased (on abdata it returns ρ̂ = 1.02; system GMM returns 0.69). Parity with xtabond2 ... robust — one-step, two-step Windmeijer, and collapsed — to machine precision. The extra level moments are an additional assumption (the deviations of the initial conditions from the long-run mean must be uncorrelated with α_i). Report the difference-in-Hansen test for the level instruments before relying on them. 'ah' — Anderson-Hsiao (1981) simple IV: one pooled instrument for the differenced lagged dependent variable instead of the block-diagonal set, chosen by ah_instrument. Consistent but inefficient, and the natural robustness check on an Arellano-Bond fit: it uses so few moments that instrument proliferation cannot be driving the answer, so a large gap between the two is informative about the instrument set rather than about the data.	'difference'
twostep	bool	Use two-step GMM with the efficient weight matrix. When robust=True the Windmeijer (2005) finite-sample correction is applied to the two-step standard errors; with robust=False the conventional (downward-biased) two-step SEs are returned and a warning is issued.	False
robust	bool	Heteroskedasticity-robust standard errors (Windmeijer-corrected for two-step). When False, the classical one/two-step VCE.	True
alpha	float	Significance level.	0.05
predetermined	list of str	Regressors that are predetermined (weakly exogenous): E[x_{is} ε_{it}] = 0 for s ≤ t but not for s > t. Their own lags from predetermined_lags (default 1 and deeper) enter as block-diagonal GMM instruments rather than as a single Δx column. Accepts the same lag-operator syntax as x. Equivalent to Stata's xtabond ..., pre() / xtabond2 ... gmm(x, lag(1 .)).	None
endogenous	list of str	Regressors correlated with the contemporaneous error. Their lags from endogenous_lags (default 2 and deeper) are used as GMM instruments. Equivalent to xtabond ..., endogenous() / xtabond2 ... gmm(x, lag(2 .)).	None
predetermined_lags	tuple(min, max)	Instrument lag windows for the two classes above. None on either element means "class default" for the minimum and "all available deeper lags" for the maximum. Lags are absolute, i.e. counted from the equation period, which is xtabond2's gmm(x, lag(a b)) convention and the one the modern literature states. Stata's older xtabond counts further lags beyond the deepest lag of the variable that appears as a regressor, so xtabond ..., pre(w, lagstruct(p, .)) maps to predetermined=['l(0/p).w'], predetermined_lags=(p + 1, None) and endogenous(w, lagstruct(p, .)) to endogenous_lags=(p + 2, None). Both mappings are pinned by reference-parity tests.	None
endogenous_lags	tuple(min, max)	Instrument lag windows for the two classes above. None on either element means "class default" for the minimum and "all available deeper lags" for the maximum. Lags are absolute, i.e. counted from the equation period, which is xtabond2's gmm(x, lag(a b)) convention and the one the modern literature states. Stata's older xtabond counts further lags beyond the deepest lag of the variable that appears as a regressor, so xtabond ..., pre(w, lagstruct(p, .)) maps to predetermined=['l(0/p).w'], predetermined_lags=(p + 1, None) and endogenous(w, lagstruct(p, .)) to endogenous_lags=(p + 2, None). Both mappings are pinned by reference-parity tests.	None
collapse	bool	Collapse the block-diagonal instrument sets — one column per lag distance instead of one per (period, distance) pair. This is Roodman's (2009) remedy for instrument proliferation: the uncollapsed count grows as O(T²), which overfits the endogenous regressor, biases the estimate toward the within estimator, and pushes the Hansen p-value toward an uninformative 1.0. Matches xtabond2, collapse.	False
orthogonal	bool	Use Arellano-Bover (1995) forward orthogonal deviations instead of first differences for the transformed equation: subtract from each observation the mean of its available future ones, scaled so the transformed errors stay serially uncorrelated. Prefer it on gappy panels — first differencing destroys the equations on both sides of a hole, forward deviations only the one at the hole. Matches xtabond2, orthogonal (difference and system, one- and two-step).	False
time_dummies	bool	Add period dummies (first period dropped) as regressors and as their own standard instruments. Roodman (2009) recommends these as a default: they absorb common shocks, which is what makes the no-cross-sectional-dependence assumption behind the moment conditions plausible.	False
cluster	str	Column to cluster the standard errors on, instead of the panel unit. The moment conditions are summed within a unit by construction, so the cluster must be at least as coarse as the unit (industry, region, cohort) and constant within it; a finer variable raises. Only the meat of the sandwich re-groups — the one-step weight (Z'HZ)^{-1} stays a within-unit object because H encodes the serial structure the transform induces. Matches xtabond2, cluster().	None
steps	int or {iterated, cue}	Generalises twostep. 1 and 2 are the one- and two-step estimators (steps=2 is exactly twostep=True); an integer above 2 repeats the same recursion — re-estimate the moment covariance at the current residuals, re-solve — that many times. 'iterated' runs the recursion to a fixed point, at which the coefficient vector and the weight matrix it implies are mutually consistent, removing the arbitrariness of stopping at two. 'cue' is the continuously-updated estimator (Hansen, Heaton & Yaron 1996), which re-evaluates the weight inside the objective and so never depends on preliminary residuals at all — the dependence the Windmeijer correction exists to patch. CUE optimises a non-convex objective numerically; converged is reported in model_info. Passing both twostep and steps raises. On a heavily over-identified fit the iterated and CUE estimates can sit well away from the two-step one; that is information about the instrument set, not noise.	None
iter_tol	(float, int)	Convergence tolerance on the maximum coefficient change, and the iteration cap, for steps='iterated' / 'cue'. Failure to converge warns rather than silently returning the last iterate.	1e-10
iter_maxiter	(float, int)	Convergence tolerance on the maximum coefficient change, and the iteration cap, for steps='iterated' / 'cue'. Failure to converge warns rather than silently returning the last iterate.	1e-10
ah_instrument	(levels, differences)	Which Anderson-Hsiao instrument to use for method='ah': y_{t-2} in levels, or its first difference Δy_{t-2}. The differences variant costs one further period of data (it reaches back to y_{t-3}) and is usually the weaker instrument, but it is valid under slightly weaker assumptions about the initial conditions. Ignored for the other methods.	'levels'
constant	bool	Include an intercept. Defaults to True for method='system' (where the level equation identifies it, as in xtabond2) and False for method='difference' (where it differences away); requesting it for difference GMM raises NotImplementedError.	None

Returns:

Type

Description

CausalResult

estimate / se are the lagged-Y (ρ₁) coefficient. detail carries the per-coefficient table (lagged Y first, then exogenous regressors, then predetermined, endogenous, time dummies and — for system GMM — _cons). model_info holds the diagnostics: n_obs (transformed-equation rows for difference GMM, level-equation rows for system GMM — the counts xtabond and xtabond2 respectively print, with n_obs_diff / n_obs_level / n_obs_total always available), n_instruments, the AR(1)/AR(2) Arellano-Bond statistics, the Sargan test (valid under homoskedasticity) and the Hansen J (heteroskedasticity-robust, reported for one-step fits too). The Sargan statistic uses xtabond's scale σ̂² = ê*'ê* / (2 (N* − k)) over the transformed rows only (level residuals still contain α_i and carry no information about σ²). xtabond2 divides by 2 N* instead, so its Sargan sits a factor N*/(N* − k) higher; the Hansen J, which has no such free scale, matches xtabond2 exactly.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> rows = []
>>> for firm in range(40):  # build a dynamic panel
...     alpha = rng.normal(0, 1)
...     y_prev = rng.normal(0, 1)
...     for t in range(8):
...         capital = rng.normal(0, 1)
...         labor = rng.normal(0, 1)
...         y = (0.5 * y_prev + 0.3 * capital + 0.2 * labor
...              + alpha + rng.normal(0, 1))
...         rows.append({'firm': firm, 'year': 2000 + t, 'output': y,
...                      'capital': capital, 'labor': labor})
...         y_prev = y
>>> df = pd.DataFrame(rows)
>>> # Arellano-Bond (difference GMM)
>>> result = sp.xtabond(df, y='output', x=['capital', 'labor'],
...                     id='firm', time='year')
>>> print(result.summary())

>>> # Two-step with Windmeijer-corrected SEs
>>> result = sp.xtabond(df, y='output', x=['capital', 'labor'],
...                     id='firm', time='year', twostep=True)

>>> # Blundell-Bond system GMM (use when the series is persistent)
>>> result = sp.xtabond(df, y='output', x=['capital', 'labor'],
...                     id='firm', time='year', method='system',
...                     twostep=True, collapse=True)

>>> # Lag operators, a predetermined regressor, and collapsed instruments
>>> result = sp.xtabond(df, y='output', x=['l(0/1).capital'],
...                     predetermined=['labor'], id='firm', time='year',
...                     collapse=True)

Notes

Arellano-Bond (1991): First-differences the equation to remove fixed effects α_i, then uses lagged levels Y_{i,t-2}, Y_{i,t-3}, ... as a block-diagonal set of GMM instruments for ΔY_{i,t-1}.

No constant / time trend is included (unlike Stata's default xtabond, which adds a _cons via a level moment). This matches Stata's xtabond ..., noconstant; the reported ρ / β coefficients are identical to Stata's _cons run when the series has no drift.

Missing values are handled per variable, not per row. A covariate that is unobserved at period t costs only the equations that need it; y_{i,t} remains available as a GMM instrument and as a lag source. (Before v1.21 a listwise dropna deleted the whole row, which amputated the instrument set whenever a lagged regressor was used — on abdata that moved ρ̂ from 0.849 to 0.660. See MIGRATION.md.)

Balanced vs gapped panels. Coefficients, standard errors, the Sargan/Hansen tests and the AR(1)/AR(2) tests are validated to machine precision against Stata for balanced panels, ragged-but-gap-free panels, and panels with interior gaps alike. A warning is emitted on gapped panels, but it is an efficiency advisory only: first differencing loses two equations per hole where forward orthogonal deviations lose one, so orthogonal=True is usually preferable there.

When cross-checking a gapped panel against xtabond2, write the instrument set on the level — gmm(y, lag(a b)), matching gmm_lags=(a, b) — rather than on a lagged expression such as gmm(L.y, lag(a-1 b-1)). The two are the same moment set on a gap-free panel and different once the panel has holes, because Stata materialises L.y row by row before xtabond2 lags it again.

Key diagnostics: - AR(1) test: Should reject (expected in first differences). - AR(2) test: Should NOT reject (validates instrument exogeneity). - Sargan / Hansen test: Should NOT reject (overidentification). Sargan (one-step) is not robust to heteroskedasticity; prefer the two-step Hansen J when that is a concern. - Instrument count: reported in model_info['n_instruments']; a warning fires when it reaches the number of units.

See Roodman (2009, Stata Journal) for practical guidance.

References

Arellano, M. and Bond, S. (1991). Some tests of specification for panel data: Monte Carlo evidence and an application to employment equations. Review of Economic Studies. [@arellano1991some] Roodman, D. (2009). How to do xtabond2: An introduction to difference and system GMM in Stata. Stata Journal. [@roodman2009xtabond]

xtdpdsys

xtdpdsys(data: DataFrame, y: str, x: Optional[Sequence[str]] = None, id: str = 'id', time: str = 'time', lags: int = 1, gmm_lags: Tuple[int, Optional[int]] = (2, None), twostep: bool = False, robust: bool = True, alpha: float = 0.05, **kwargs) -> CausalResult

Blundell-Bond (1998) system GMM for dynamic panels.

Thin alias for sp.xtabond(..., method='system'), named after Stata's xtdpdsys so the estimator is discoverable under the name applied researchers already use. Every keyword of :func:xtabond is accepted.

Parameters:

Name	Type	Description	Default
data	DataFrame	As in :func:xtabond.	required
y	DataFrame	As in :func:xtabond.	required
x	DataFrame	As in :func:xtabond.	required
id	DataFrame	As in :func:xtabond.	required
time	DataFrame	As in :func:xtabond.	required
lags	DataFrame	As in :func:xtabond.	required
gmm_lags	DataFrame	As in :func:xtabond.	required
twostep	DataFrame	As in :func:xtabond.	required
robust	DataFrame	As in :func:xtabond.	required
alpha	DataFrame	As in :func:xtabond.	required
**kwargs		Forwarded to :func:xtabond — predetermined, endogenous, collapse, time_dummies, constant, and the instrument lag windows.	{}

Returns:

Type	Description
CausalResult	As :func:xtabond, including the _cons row the level equation identifies.

Notes

System GMM stacks the level equation alongside the first-differenced one and instruments it with lagged differences. It is the right default when the dependent variable is persistent — the case in which the lagged levels that difference GMM relies on are weak instruments and the estimate collapses toward (or past) unity.

The extra moments buy that power at the price of an extra assumption: the deviation of each unit's initial condition from its long-run mean must be uncorrelated with the fixed effect. Test it — the level instruments have their own difference-in-Hansen statistic.

Validated to machine precision against xtabond2 ..., robust (one-step, two-step Windmeijer, and collapsed) on the abdata panel; see tests/reference_parity/test_dynpanel_abdata_parity.py.

Examples:

>>> import statspai as sp
>>> res = sp.xtdpdsys(df, y='n', x=['w', 'k'], id='id',
...                   time='year', twostep=True, collapse=True)

References

Blundell, R. and Bond, S. (1998). Initial conditions and moment restrictions in dynamic panel data models. Journal of Econometrics. [@blundell1998initial] Roodman, D. (2009). How to do xtabond2. Stata Journal. [@roodman2009xtabond]

gmm

gmm(moment_fn: Callable[[ndarray, Optional[DataFrame]], Any], theta0: ndarray, data: Optional[DataFrame] = None, W: Optional[ndarray] = None, method: str = 'twostep', se: str = 'robust', maxiter: int = 200, tol: float = 1e-08, param_names: Optional[List[str]] = None, alpha: float = 0.05, jacobian: Optional[Callable[[ndarray, Optional[DataFrame]], Any]] = None, vcov: str = 'mds', cluster: Optional[Any] = None, hac_bandwidth: Optional[int] = None, center: bool = False) -> EconometricResults

General GMM estimator for arbitrary moment conditions.

Minimizes Q(theta) = gbar(theta)' W gbar(theta) where gbar(theta) = (1/n) sum_i g_i(theta).

Parameters:

Name	Type	Description	Default
moment_fn	callable	g(theta, data) -> ndarray of shape (n, q): the moment contribution of each observation.	required
theta0	ndarray	Starting values. Also the expansion point for the closed form when the moments are affine.	required
data	DataFrame	Passed through to moment_fn.	None
W	ndarray	Weighting matrix (q, q) for the first step. Defaults to the identity. Supplying one and asking for se='unadjusted' warns unless it happens to be efficient.	None
method	(onestep, twostep, iterative, cue)		'onestep'
se	(robust, unadjusted)	'robust' returns the sandwich, valid for any W. 'unadjusted' returns the efficient-GMM variance (D'WD)^{-1}/n, which describes the estimator only at the efficient weight; otherwise it warns.	'robust'
maxiter	int		200
tol	float		1e-8
param_names	list of str		None
alpha	float		0.05
jacobian	callable	D(theta, data) -> ndarray of shape (q, k), the derivative of the average moment. Supplying it removes finite-difference error from the standard errors and from the affine test.	None
vcov	(mds, iid, hac, cluster)	Estimator for the moment covariance S.	'mds'
cluster	array - like	Group labels, one per moment row. Required when vcov='cluster'.	None
hac_bandwidth	int	Bartlett bandwidth for vcov='hac'; the kernel vanishes at lag == hac_bandwidth.	None
center	bool	Centre the moments before forming S. False matches Stata, True matches R's gmm.	False

Returns:

Type	Description
EconometricResults	diagnostics carries J_stat / J_df / J_p, converged, and n_iter (0 when the closed form was used).

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> z1, z2, u = rng.normal(size=n), rng.normal(size=n), rng.normal(size=n)
>>> x1 = 0.7 * z1 + 0.5 * z2 + u + rng.normal(size=n)
>>> y = 1.0 + 2.0 * x1 + u + rng.normal(size=n)
>>> df = pd.DataFrame({'y': y, 'x1': x1, 'z1': z1, 'z2': z2})
>>>
>>> def moment_fn(theta, data):
...     y, X, Z = data['y'].values, data[['x1']].values, data[['z1', 'z2']].values
...     X_full = np.column_stack([np.ones(len(y)), X])
...     resid = y - X_full @ theta
...     Z_full = np.column_stack([np.ones(len(y)), Z])
...     return resid[:, np.newaxis] * Z_full
>>>
>>> result = sp.gmm(moment_fn, theta0=np.zeros(2), data=df,
...                 param_names=['_cons', 'x1'])
>>> bool(result is not None)
True
>>> result.diagnostics['n_iter']  # affine moments -> closed form
0

xtlsdvc

xtlsdvc(data: DataFrame, y: str, x: Optional[Sequence[str]] = None, id: str = 'id', time: str = 'time', initial: str = 'ab', bias_order: int = 2, alpha: float = 0.05, bootstrap: int = 0, seed: Optional[int] = None) -> CausalResult

Bias-corrected LSDV (LSDVC) estimator for dynamic panels.

Equivalent to Stata's xtlsdvc (Bruno 2005).

Parameters:

Name	Type	Description	Default
data	DataFrame	Panel in long format.	required
y	str	Dependent variable. Its first lag is added automatically — do not pass a hand-built lag in x.	required
x	list of str	Strictly exogenous regressors. Accepts the same lag-operator syntax as :func:~statspai.gmm.arellano_bond.xtabond ("l(0/1).w").	None
id	str	Unit and period identifiers.	'id'
time	str	Unit and period identifiers.	'id'
initial	('ab', 'ah', 'bb')	Consistent estimator used to evaluate the bias expression: Arellano-Bond, Anderson-Hsiao, or Blundell-Bond system GMM. The correction is only as good as this input, and the three can disagree materially on a persistent series — 'bb' is the safer choice there, for the same reason system GMM is.	'ab'
bias_order	(1, 2, 3)	How many terms of the Bun-Kiviet expansion to subtract: O(1/T), O(1/NT), or O(1/NT^2). Bun & Kiviet (2003) find diminishing returns beyond the second.	1
alpha	float	Significance level.	0.05
bootstrap	int	Number of parametric-bootstrap replications for the standard errors. The reported analytic standard errors are the LSDV ones and do not account for the bias correction, exactly as in Stata's xtlsdvc; set this to get honest ones. 0 skips the bootstrap and warns.	0
seed	int	Seed for the bootstrap.	None

Returns:

Type	Description
CausalResult	estimate / se are the corrected lagged-Y coefficient. detail holds the full corrected table, the uncorrected LSDV coefficients, and the subtracted bias.

Notes

The bias expression is the Bun-Kiviet/Bruno one, evaluated per unit on the balanced period grid with a selection matrix marking the observed rows (so unbalanced panels are handled as in Bruno 2005). Validated against Stata's xtlsdvc at all three bias orders and all three initial estimators.

When to prefer this over GMM. Small N with small T — the regime where the GMM instrument count rivals the number of units and the Hansen test stops being informative. LSDVC uses no instruments at all, so it cannot be undone by instrument proliferation; the price is that it assumes strict exogeneity and homoskedasticity, which GMM does not.

Examples:

>>> import statspai as sp
>>> res = sp.xtlsdvc(df, y='n', x=['w', 'k'], id='id',
...                  time='year', initial='bb', bootstrap=200)

References

Kiviet, J.F. (1995). Journal of Econometrics 68(1), 53-78. [@kiviet1995bias] Bruno, G.S.F. (2005). Economics Letters 87(3), 361-366. [@bruno2005approximating]