fit_gmm_moments

fit_gmm_moments(
    moment_conditions,
    initial,
    lower=None,
    upper=None,
    sample_size=None,
    jacobian=None,
    penalty=None,
    optimizer='BFGS',
    strategy='Iterative',
    max_gmm_iterations=0,
)

Fit your own moment conditions by the generalized method of moments.

Estimates parameters by GMM against moment conditions you write. This is the constructor the upstream C# library itself exposes, GeneralizedMethodOfMoments(momentConditionFunction, ...): :func:~corehydropy.fit_gmm can only fit a :func:~corehydropy.model_bulletin17c model (the single implementation of the IGMMModel interface it takes), while this function takes any moment conditions you can write down.

The moment condition function

moment_conditions(parameters) is called with one list of numbers, as long as initial, and must return two things: the tuple (g, s), or a dict with keys "g" and "s".

g is the sample mean of the moment conditions at those parameters: a sequence of q numbers, one per moment condition. GMM drives this towards zero. With one moment condition it may be written as the bare number, which is R’s spelling for it – except when s is also written bare, where (number, number) cannot be told apart from a flat [g0, g1] and is refused by name.
s is their covariance: a q by q matrix, written as a sequence of ROWS (a list of lists, or a 2-D numpy array). In two-step and iterative GMM the optimal weighting matrix is its inverse.

Both are required and both are checked by name, because returning the wrong thing here is the likeliest mistake on this surface. A two-parameter method-of-moments fit of a Normal, whose answer is the sample mean and the population variance::

def moments(p):
    a = [xi - p[0] for xi in x]
    b = [ai * ai - p[1] for ai in a]
    n = len(x)
    mean = lambda v: sum(v) / n
    return ([mean(a), mean(b)],
            [[mean([ai * ai for ai in a]), mean([ai * bi for ai, bi in zip(a, b)])],
             [mean([ai * bi for ai, bi in zip(a, b)]), mean([bi * bi for bi in b])]])

q (the number of moment conditions) is measured by calling your function once at initial, so there is no argument to get wrong. When q equals the number of parameters the fit is just-identified and .j_stat_pval comes back None – see Returns.

Parameters

Name	Type	Description	Default
moment_conditions	callable	The required function described above.	required
initial	sequence of float	Starting values, one per parameter. Its length IS the parameter count.	required
lower	sequence of float	Parameter bounds, the same length as `initial`, with every starting value inside them. Both are required despite the `None` default: every optimizer this dispatches to takes a box, and the numerical Jacobian’s step selection is bounds-aware.	`None`
upper	sequence of float	Parameter bounds, the same length as `initial`, with every starting value inside them. Both are required despite the `None` default: every optimizer this dispatches to takes a box, and the numerical Jacobian’s step selection is bounds-aware.	`None`
sample_size	int	The number of observations behind the moment conditions. Required, and only you know it: your function hands over averages, not data. The sandwich covariance divides by it, so the standard errors scale as `1 / sqrt(sample_size)`.	`None`
jacobian	callable	An analytic Jacobian of the moment conditions, called as `jacobian(parameters)` and returning a q by p matrix – one ROW per moment condition, one COLUMN per parameter. `None`, the default, uses the ported bounds-aware finite-difference Jacobian, which costs two extra `moment_conditions` calls per parameter per gradient.	`None`
penalty	callable	A penalty added to the GMM objective, called as `penalty(parameters)` and returning one number. Ridge-type regularization, and the only way to fit a model with more parameters than moment conditions (which is otherwise refused as under-identified). Return `0.0` for no penalty. Note the ported half-quadratic convention: with a penalty the objective is `0.5 * g'Wg + penalty`, so a penalty should carry its own `1/2`.	`None`
optimizer	str	One of `"BFGS"` (matching `GeneralizedMethodOfMoments`’s own class default), `"NelderMead"`, `"Brent"`, `"Powell"`, `"DifferentialEvolution"`, `"MultilevelSingleLinkage"`.	`"BFGS"`
strategy	('Iterative', 'OneStep', 'TwoStep')	GMM estimation strategy. `"OneStep"` is refused for an over-identified problem, matching the estimator.	`"Iterative"`
max_gmm_iterations	int	Maximum number of GMM iterations; `0` keeps the estimator’s own default cap.	`0`

Returns

Name	Type	Description
	Fit	The same object :func:`~corehydropy.fit_gmm` returns, with `.method == "GMM"`, so `.parameters`, `.covariance`, `repr()` and `.summary()` behave identically. Parameters are named `p1`, `p2`, … since your moment conditions name nothing. Method of moments computes no likelihood surface, so `.log_likelihood`, `.aic` and `.bic` are `None` and `.confint()` raises, as they do for :func:`~corehydropy.fit_gmm`. `.j_stat` is Hansen’s J and `.j_stat_pval` its p-value, which is `None` whenever the fit is just-identified (as many moment conditions as parameters): zero degrees of freedom leaves no over-identifying restriction to test, and `.summary()` says so rather than showing a figure. `.j_stat` itself is not a goodness-of-fit number you can read there either. The residual covariance it is scaled by is singular – it has rank q - p, so it is exactly zero when the fit is just-identified – and inverting it amplifies the optimizer’s convergence tolerance rather than any property of your data. The result varies by many orders of magnitude and in sign between optimizers, and between this package and the C# library it ports, on fits whose parameters agree to ten significant figures. Sometimes it cannot be computed at all, and then it comes back `nan` rather than failing the fit. Over-identifying the model restores the p-value but not `.j_stat`, since the rank deficiency only shrinks from q to q - p. `.degree_of_freedom` and `.number_of_moment_conditions` carry q - p and q. :func:`~corehydropy.fit_diagnostics` and :func:`~corehydropy.quantile_variance` are not available for this fit – both need the model a `fit_gmm` fit carries.

Notes

There is no random number generator anywhere in this fit, so a repeated call returns the identical numbers. Across languages the guarantee is the usual one for this surface: the optimizer arithmetic all happens in the shared C++ core, but g and s are computed by your own Python code, and R and Python do not guarantee identical rounding for the same formula. Arithmetic (+ - * /) is IEEE-deterministic and does reproduce; log, exp and friends come from each platform’s math library. R’s sum() and mean() accumulate in extended precision where Python’s do not, so an explicit loop is the portable spelling.

Examples

>>> import corehydropy as ch
>>> x = [4.1, 5.2, 4.8, 5.5, 4.9, 5.1, 5.3, 4.7]
>>> def moments(p):
...     n = len(x)
...     a = [xi - p[0] for xi in x]
...     b = [ai * ai - p[1] for ai in a]
...     saa = sbb = sab = 0.0
...     for ai, bi in zip(a, b):
...         saa += ai * ai
...         sab += ai * bi
...         sbb += bi * bi
...     return ([sum(a) / n, sum(b) / n],
...             [[saa / n, sab / n], [sab / n, sbb / n]])
>>> f = ch.fit_gmm_moments(moments, initial=[5.0, 0.5], lower=[0.0, 0.001],
...                        upper=[10.0, 10.0], sample_size=len(x))
>>> round(f.parameters["p1"], 6), round(f.parameters["p2"], 6)
(4.95, 0.165)
>>> f.j_stat_pval is None  # just-identified, so there is nothing to test
True