Fit your own moment conditions by GMM

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

Usage

fit_gmm_moments(
  moment_conditions,
  initial,
  lower = NULL,
  upper = NULL,
  sample_size,
  jacobian = NULL,
  penalty = NULL,
  optimizer = "BFGS",
  strategy = "Iterative",
  max_gmm_iterations = 0L
)

Arguments

moment_conditions: the required function described above.
initial: numeric vector of starting values, one per parameter. Its length IS the parameter count.
lower, upper: numeric vectors of parameter bounds, the same length as initial, with every starting value inside them. Both are required: every optimizer this dispatches to takes a box, and the numerical Jacobian's step selection is bounds-aware.
sample_size: 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).
jacobian: an optional analytic Jacobian of the moment conditions, called as jacobian(parameters) and returning a q by p numeric matrix – one ROW per moment condition, one COLUMN per parameter. NULL, the default, uses the ported bounds-aware finite-difference Jacobian, which costs two extra moment_conditions calls per parameter per gradient.
penalty: an optional 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 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.
optimizer: one of "BFGS" (default, matching GeneralizedMethodOfMoments's own class default), "NelderMead", "Brent", "Powell", "DifferentialEvolution", "MultilevelSingleLinkage".
strategy: GMM estimation strategy: "Iterative" (default), "OneStep", or "TwoStep". "OneStep" is refused for an over-identified problem, matching the estimator.
max_gmm_iterations: maximum number of GMM iterations; 0 (default) keeps the estimator's own default cap.

Value

An object of class corehydro_fit with method == "GMM" – the same object fit_gmm() returns, so coef(), vcov(), print() 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 NA and confint() errors, as they do for fit_gmm(). $j_stat is Hansen's J and $j_stat_pval its p-value, which is NA whenever the fit is just-identified (as many moment conditions as parameters): zero degrees of freedom leaves no over-identifying restriction to test, and print() 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 NA 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. fit_diagnostics() and quantile_variance() are not available for this fit – both need the model a fit_gmm() fit carries.

The moment condition function

moment_conditions(parameters) is called with one numeric vector, as long as initial, and must return a list with two elements:

g: the sample mean of the moment conditions at those parameters: a numeric vector of length q, one entry per moment condition. GMM drives this towards zero.
s: their covariance: a q by q numeric matrix. 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:

avg <- function(v) {  # not mean(): see "Reproducing a run in Python" below
  total <- 0
  for (value in v) total <- total + value
  total / length(v)
}
moments <- function(p) {
  a <- x - p[1]
  b <- a * a - p[2]
  list(g = c(avg(a), avg(b)),
       s = matrix(c(avg(a * a), avg(a * b), avg(a * b), avg(b * b)), 2, 2))
}

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 NA – see the return value below.

Reproducing a run in Python

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

x <- c(4.1, 5.2, 4.8, 5.5, 4.9, 5.1, 5.3, 4.7)
# An explicit loop rather than mean(), so the same formula returns the same bits in Python:
# R's mean() and sum() accumulate in extended precision and Python's do not.
avg <- function(v) {
  total <- 0
  for (value in v) total <- total + value
  total / length(v)
}
moments <- function(p) {
  a <- x - p[1]
  b <- a * a - p[2]
  list(
    g = c(avg(a), avg(b)),
    s = matrix(c(avg(a * a), avg(a * b), avg(a * b), avg(b * b)), 2, 2)
  )
}
f <- fit_gmm_moments(moments,
  initial = c(5, 0.5), lower = c(0, 0.001), upper = c(10, 10),
  sample_size = length(x)
)
round(coef(f), 6) # the sample mean and the population variance
#>    p1    p2 
#> 4.950 0.165 
f$j_stat_pval # NA: just-identified, so there is nothing to test
#> [1] NA