Summary
GeneralizedMethodOfMoments.PostProcess computes Hansen's J statistic as g' V^-1 g by calling V.Inverse() on the moment residual covariance built by GetMomentResidualCovariance. That matrix is singular by construction, so the inverse amplifies the optimizer's convergence tolerance and the resulting J describes which optimizer ran rather than how well the model fits.
V = S - D(D'S^-1 D)^-1 D'. Multiplying both sides by S^(-1/2) gives S^(-1/2) V S^(-1/2) = I - P where P is a projection of rank p, so rank(V) = q - p exactly. V is singular for every q and p. Over-identifying moves the rank from 0 to 1; it does not make V invertible.
Because V.Inverse() is an LU inverse of a rank-deficient matrix, what it amplifies is the first-order-condition residual left by the optimizer, which is a convergence tolerance around 1e-8, not floating-point rounding.
Steps to reproduce
Fit any over-identified model (q > p) through the delegate-based GeneralizedMethodOfMoments constructor and read JStat / JStatPval off the result. A three-moment fit of a two-parameter Normal is enough. Note that this is not reachable through Bulletin17CDistribution, which is always just-identified, so it only appears with user-supplied moment conditions.
Measured at the fitted parameters of one such fit:
- V has singular values
[1.9e-02, 9.8e-19, 2.3e-19], numerical rank 1, condition number 8.1e16.
- Perturbing S or D by 1e-11 swings J across
+892, -1061, +668, -820.
- Across four optimizers on parameters agreeing to 1e-5, J spans
-129.46, 1268.6, 1.2e+08, 3.8e+06.
- Replacing the inverse with a Moore-Penrose pseudo-inverse gives
g' V^+ g = 2.3466, which matches the textbook n * g' S^-1 g = 2.3466 on the same fit. S, D and g are therefore all correct; only the inversion is at fault.
I found this while validating a C++ port of Numerics and RMC.BestFit against the real libraries. The port is line-identical to GetMomentResidualCovariance and also uses an LU inverse, and on the same fit the C# library returns J = 214.59 with a p-value of 0 where the port returns J = -129.46 with a p-value of 1, from parameters that agree to 2e-11. Two implementations of the same formula disagreeing that far on inputs that close is the symptom of the rank deficiency rather than of either implementation.
Suggested fix
Compute J through a pseudo-inverse of V, or equivalently as n * g' S^-1 g, in place of the V.Inverse() call in PostProcess. Either form is stable against the rank deficiency V carries by construction, and both are reproducible across compilers and optimizers.
The just-identified case (q == p, rank(V) = 0) is the degenerate end of the same fact. There JStatPval is already NaN because the degrees of freedom are zero, so nothing downstream changes; it may still be worth declining to report JStat itself at zero degrees of freedom, since no value there is meaningful.
Environment
RMC-BestFit c2e6192 (v2.0.0) and Numerics 2a0357a (v2.1.4), library rather than desktop app. Source references: RMC.BestFit/Estimation/GeneralizedMethodOfMoments.cs, GetMomentResidualCovariance around lines 945-975 and the var Vinv = V.Inverse(); in PostProcess around line 2511.
Summary
GeneralizedMethodOfMoments.PostProcesscomputes Hansen's J statistic asg' V^-1 gby callingV.Inverse()on the moment residual covariance built byGetMomentResidualCovariance. That matrix is singular by construction, so the inverse amplifies the optimizer's convergence tolerance and the resulting J describes which optimizer ran rather than how well the model fits.V = S - D(D'S^-1 D)^-1 D'. Multiplying both sides byS^(-1/2)givesS^(-1/2) V S^(-1/2) = I - PwherePis a projection of rankp, sorank(V) = q - pexactly. V is singular for everyqandp. Over-identifying moves the rank from 0 to 1; it does not make V invertible.Because
V.Inverse()is an LU inverse of a rank-deficient matrix, what it amplifies is the first-order-condition residual left by the optimizer, which is a convergence tolerance around 1e-8, not floating-point rounding.Steps to reproduce
Fit any over-identified model (
q > p) through the delegate-basedGeneralizedMethodOfMomentsconstructor and readJStat/JStatPvaloff the result. A three-moment fit of a two-parameter Normal is enough. Note that this is not reachable throughBulletin17CDistribution, which is always just-identified, so it only appears with user-supplied moment conditions.Measured at the fitted parameters of one such fit:
[1.9e-02, 9.8e-19, 2.3e-19], numerical rank 1, condition number 8.1e16.+892,-1061,+668,-820.-129.46,1268.6,1.2e+08,3.8e+06.g' V^+ g = 2.3466, which matches the textbookn * g' S^-1 g = 2.3466on the same fit. S, D and g are therefore all correct; only the inversion is at fault.I found this while validating a C++ port of Numerics and RMC.BestFit against the real libraries. The port is line-identical to
GetMomentResidualCovarianceand also uses an LU inverse, and on the same fit the C# library returnsJ = 214.59with a p-value of 0 where the port returnsJ = -129.46with a p-value of 1, from parameters that agree to 2e-11. Two implementations of the same formula disagreeing that far on inputs that close is the symptom of the rank deficiency rather than of either implementation.Suggested fix
Compute J through a pseudo-inverse of V, or equivalently as
n * g' S^-1 g, in place of theV.Inverse()call inPostProcess. Either form is stable against the rank deficiency V carries by construction, and both are reproducible across compilers and optimizers.The just-identified case (
q == p,rank(V) = 0) is the degenerate end of the same fact. ThereJStatPvalis alreadyNaNbecause the degrees of freedom are zero, so nothing downstream changes; it may still be worth declining to reportJStatitself at zero degrees of freedom, since no value there is meaningful.Environment
RMC-BestFit
c2e6192(v2.0.0) and Numerics2a0357a(v2.1.4), library rather than desktop app. Source references:RMC.BestFit/Estimation/GeneralizedMethodOfMoments.cs,GetMomentResidualCovariancearound lines 945-975 and thevar Vinv = V.Inverse();inPostProcessaround line 2511.