GMM with moment evaluator program reports no (too small) standard errors

[Adjmal Sirak]

Hello Statalist,
I am using Stata MP, version 13.1 on Windows 10. I have a problem with a GMM estimation using a moment evaluator program. Before I get into the details, let me provide you with a bit of backup information.

As part of my research project I have specified a model that could be called structural. It contains 137 parameters identified by 181 moment conditions. I am interested in quantities that are nonlinear functions of these parameters and I want to obtain bootstrapped confidence intervals. In order to speed up the bootstrap, I want to provide analytical derivatives. It is impossible to do so in the interactive mode of the gmm command, as the number of options that can be passed is limited to 70. This is why I use a GMM moment evaluator program.

I have set up the program, but when I pass it to the gmm command, I obtain point estimates of the parameters, but no standard errors. This suggests that something goes wrong in my moment evaluator program. I have looked into that and found out that the problem persists even in much more parsimonious settings, as I will demonstrate now.

A much cooked-down version of my problem could be the following. There are two random variables X1, X2 with a common mean mu. The task is to estimate mu. Hence there is one parameter and there are two moment equations.

E(X1-mu) = 0

E(X2-mu) = 0

The following Stata code implements the GMM estimation of mu, once in the interactive mode and once using a moment evaluator program. The user may set the true value of mu via the local ``mu'' and the number of random variables (two in the example given) may be set by the global ``num_vars''.

clear all
set seed 12345
set obs 1000

* set mean of random vars
global num_vars = 2
local mu = 1
foreach nn of numlist 1/$num_vars {
    quietly gen x`nn' = `mu' + rnormal()
}

* define moment evaluator program to estimate mu
capture program drop mymep
program define mymep
    version 13.1
    syntax varlist if, at(name) [derivatives(varlist)]
        
    foreach nn of numlist 1/$num_vars {
        * moment equations
        local mom_eq`nn' : word `nn' of `varlist'
        * calculate sample moments
        tempvar foo
        quietly egen double `foo' = mean(x`nn' - `at'[1,1]) `if'
        * replace in varlist
        quietly replace `mom_eq`=`nn''' = `foo'
    }
    
    * derivatives
    if "`derivatives'" == "" {
        exit
    }
    foreach nn of numlist 1/$num_vars {
        local dxdmu : word `nn' of `derivatives'
        quietly replace `dxdmu' = -1
    }
end

* GMM estimation of mu in the interactive mode
local eq_list = ""
local inst_list = ""
local deriv_list = ""
foreach nn of numlist 1/$num_vars {
    local eq_list = "`eq_list'" + "(eq`nn': (x`nn' - {mu})) "
    local inst_list = "`inst_list'" + "eq`nn' "
    local deriv_list = "`deriv_list'" + "deriv(eq`nn'/mu = -1) "
}

gmm `eq_list', inst(`inst_list': ) `deriv_list' winitial(unadjusted, independent)
est store gmm1
mat b_gmm1 = e(b)
mat V_gmm1 = e(V)

* GMM estimation of mu using the moment evaluator program
local eq_list = ""
foreach nn of numlist 1/$num_vars {
    local eq_list = "`eq_list'" + "eq`nn' "
}
gmm mymep, eq(`eq_list') inst(`eq_list': ) param(mu) ///
    winitial(unadjusted, independent) hasderivatives
est store gmm2
mat b_gmm2 = e(b)
mat V_gmm2 = e(V)

* compare point estimates
mat list b_gmm1
mat list b_gmm2
* compare VCEs
mat list V_gmm1
mat list V_gmm2

* compare ereturn
est restore gmm1
ereturn list
est restore gmm2
ereturn list

Here is the output from calling gmm in the interactive version:

. * GMM estimation of mu in the interactive mode
. local eq_list = ""

. local inst_list = ""

. local deriv_list = ""

. foreach nn of numlist 1/$num_vars {
  2.         local eq_list = "`eq_list'" + "(eq`nn': (x`nn' - {mu})) "
  3.         local inst_list = "`inst_list'" + "eq`nn' "
  4.         local deriv_list = "`deriv_list'" + "deriv(eq`nn'/mu = -1) "
  5. }
.
. gmm `eq_list', inst(`inst_list': ) `deriv_list' winitial(unadjusted, independent)

Step 1
Iteration 0:   GMM criterion Q(b) =  2.0713458  
Iteration 1:   GMM criterion Q(b) =  .00067259  
Iteration 2:   GMM criterion Q(b) =  .00067259  

Step 2
Iteration 0:   GMM criterion Q(b) =  .00062707  
Iteration 1:   GMM criterion Q(b) =  .00062701  

GMM estimation

Number of parameters =   1
Number of moments    =   2
Initial weight matrix: Unadjusted                     Number of obs  =    1000
GMM weight matrix:     Robust

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /mu |   1.017355   .0224233    45.37   0.000      .973406    1.061304
------------------------------------------------------------------------------
Instruments for equation 1: _cons
Instruments for equation 2: _cons

. est store gmm1

. mat b_gmm1 = e(b)

. mat V_gmm1 = e(V)

And here is the output from calling gmm using the moment evaluator program:

. * GMM estimation of mu using the moment evaluator program
. local eq_list = ""

. foreach nn of numlist 1/$num_vars {
  2.         local eq_list = "`eq_list'" + "eq`nn' "
  3. }

. gmm mymep, eq(`eq_list') inst(`eq_list': ) param(mu) ///
>         winitial(unadjusted, independent) hasderivatives

Step 1
Iteration 0:   GMM criterion Q(b) =  2.0713458  
Iteration 1:   GMM criterion Q(b) =  .00067259  
Iteration 2:   GMM criterion Q(b) =  .00067259  

Step 2
Iteration 0:   GMM criterion Q(b) =          1  
Iteration 1:   GMM criterion Q(b) =  9.163e-30  
Iteration 2:   GMM criterion Q(b) =  9.163e-30  

GMM estimation

Number of parameters =   1
Number of moments    =   2
Initial weight matrix: Unadjusted                     Number of obs  =    1000
GMM weight matrix:     Robust

------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /mu |   .9991766   1.76e-18  5.7e+17   0.000     .9991766    .9991766
------------------------------------------------------------------------------
Instruments for equation 1: _cons
Instruments for equation 2: _cons

. est store gmm2

. mat b_gmm2 = e(b)

. mat V_gmm2 = e(V)

Here are the comparisons of point estimates, VCEs, and the ereturn output:

. * compare point estimates
. mat list b_gmm1

symmetric b_gmm1[1,1]
           mu:
        _cons
y1  1.0173548

. mat list b_gmm2

symmetric b_gmm2[1,1]
           mu:
        _cons
y1  .99917656

. * compare VCEs
. mat list V_gmm1

symmetric V_gmm1[1,1]
                mu:
             _cons
mu:_cons  .0005028

. mat list V_gmm2

symmetric V_gmm2[1,1]
                 mu:
              _cons
mu:_cons  3.081e-36

.
. * compare ereturn
. est restore gmm1
(results gmm1 are active now)

. ereturn list

scalars:
               e(rank) =  1
                  e(N) =  1000
                  e(Q) =  .0006270145626651
                  e(J) =  .6270145626650803
               e(J_df) =  1
                e(k_1) =  1
                e(k_2) =  1
          e(converged) =  1
         e(has_xtinst) =  0
               e(type) =  1
               e(n_eq) =  2
                  e(k) =  1
          e(n_moments) =  2
              e(k_aux) =  1
         e(k_eq_model) =  0
               e(k_eq) =  1

macros:
         e(properties) : "b V"
             e(sexp_1) : "(x1 - {mu} )"
           e(params_1) : "mu"
             e(inst_1) : "_cons"
             e(sexp_2) : "(x2 - {mu} )"
           e(params_2) : "mu"
             e(inst_2) : "_cons"
             e(params) : "mu"
            e(vcetype) : "Robust"
                e(vce) : "robust"
            e(wmatrix) : "robust"
          e(estimator) : "twostep"
              e(winit) : "Unadjusted"
          e(technique) : "gn"
            e(eqnames) : "eq1 eq2"
       e(marginsnotok) : "_ALL"
            e(predict) : "gmm_p"
          e(estat_cmd) : "gmm_estat"
                e(cmd) : "gmm"
            e(cmdline) : "gmm (eq1: (x1 - {mu})) (eq2: (x2 - {mu})) , inst(eq1 eq2 : ) deriv(eq1/mu = -1) de.."
    e(_estimates_name) : "gmm1"

matrices:
                  e(b) :  1 x 1
                  e(V) :  1 x 1
                  e(W) :  2 x 2
                  e(S) :  2 x 2
       e(V_modelbased) :  1 x 1
               e(init) :  1 x 1

functions:
             e(sample)  

. est restore gmm2
(results gmm2 are active now)

. ereturn list

scalars:
               e(rank) =  0
                  e(N) =  1000
                  e(Q) =  9.16303687807e-30
                  e(J) =  9.16303687807e-27
               e(J_df) =  1
          e(converged) =  1
         e(has_xtinst) =  0
               e(type) =  2
               e(n_eq) =  2
                  e(k) =  1
          e(n_moments) =  2
              e(k_aux) =  1
         e(k_eq_model) =  0
               e(k_eq) =  1

macros:
         e(properties) : "b V"
             e(inst_1) : "_cons"
             e(inst_2) : "_cons"
             e(params) : "mu"
           e(evalprog) : "mymep"
            e(vcetype) : "Robust"
                e(vce) : "robust"
            e(wmatrix) : "robust"
          e(estimator) : "twostep"
              e(winit) : "Unadjusted"
          e(technique) : "gn"
            e(eqnames) : "eq1 eq2"
       e(marginsnotok) : "_ALL"
            e(predict) : "gmm_p"
          e(estat_cmd) : "gmm_estat"
                e(cmd) : "gmm"
            e(cmdline) : "gmm mymep, eq(eq1 eq2 ) inst(eq1 eq2 : ) param(mu) winitial(unadjusted, independen.."
    e(_estimates_name) : "gmm2"

matrices:
                  e(b) :  1 x 1
                  e(V) :  1 x 1
                  e(W) :  2 x 2
                  e(S) :  2 x 2
       e(V_modelbased) :  1 x 1
               e(init) :  1 x 1

functions:
             e(sample)  

.
end of do-file

Unlike my actual application the gmm reports standard errors when called with the moment evaluator program, but similar problems arise. As you can see, both calls of the gmm command yield point estimates that lie reasonably close to the true value of mu. However, when using the moment evaluator program there are some strange things going on:

1. Whereas in the interactive mode the standard errors are of the order of magnitude of 10^-2, they are of the order of magnitude of 10^-18 when using the moment evaluator program.
2. In the interactive mode the GMM criterion converges at an order of magnitude of 10^-4, whereas it is at 10^-30 when using the moment evaluator program.
3. When looking at the ereturn output, we find that the rank is reported as 1 in the interactive mode, but 0 when using the moment evaluator program.

In my actual application the diagonal elements of the VCE matrix are zero and no standard errors are reported at all.
Dear Statalist, can somebody please provide any insight into what might be the problem with my moment evaluator program or how to proceed with debugging? Your help will be gratefully appreciated.
Kind regards,
Adjmal

[Sebastian Kripfganz]

Hi Adjmal,

In your moment-evaluator program you need to replace the line

        quietly egen double `foo' = mean(x`nn' - `at'[1,1]) `if'

by

        quietly gen double `foo' = x`nn' - `at'[1,1] `if'

As an aside:
You might want to use version control before setting the seed to make sure that the same random numbers are drawn in different Stata versions (otherwise the results are differerent e.g. under Stata 14.2):

version 13.1
set seed 12345

Best wishes,
Sebastian

[Adjmal Sirak]

So the moment evaluator program is supposed to return observation-level output and not the sample moments, I see. Thank you very much for the quick reply.
Best regards,
Adjmal.