Problem with moptimize() and gauss-newton alogrithm when implementing a GMM estimator: covariance matrix has missing values

[Christophe Kolodziejczyk]

I am trying to estimate a model by the GMM with moptimize() and the Gauss-Newton algorithm. Although I can get a solution the returned covariance matrix is filled with missing values.
I have tried to estimate a linear regression model by GMM to see if it was model specific, but I get a similar result. The coefficients are the same compared to the command regress but the covariance matrix has missing values. What am I doing wrong?
I've got this problem on Stata IC versions 12, 13 and 14 on Windows.
Best regards
Christophe

. sysuse auto, clear
(1978 Automobile Data)

.
. clear mata

. mata:
------------------------------------------------- mata (type end to exit) -----------------------------------------------------------------
:
: void Reg(M, todo, b, r, S)
>                 {
>                
>                                 real scalar phi1
>                                 real scalar phi2
>                                 real colvector w, xb, zg, e, iota, t, u
>                                 real scalar N
>                                 real matrix W, Z
>                                
>                                
>                                 xb = moptimize_util_xb(M,b,1)
>                                 N = rows(xb)
>
>                                
>                                 u = moptimize_util_depvar(M,1) :- xb
>                                 W = moptimize_util_userinfo(M, 1)
>                                
>                                 phi1 = quadcross(W, 1, u, 0)
>                                 r = (phi1)                      
>                
>                
>                 }
note: argument todo unused
note: argument S unused
note: variable phi2 unused
note: variable w unused
note: variable zg unused
note: variable e unused
note: variable iota unused
note: variable t unused
note: variable N set but not used
note: variable Z unused

:
: Y = st_data(.,"weight")

: X = st_data(.,"mpg")

: M = moptimize_init()

: moptimize_init_evaluatortype(M,"q0")

: moptimize_init_evaluator(M,&Reg())

: moptimize_init_technique(M,"gn")

: moptimize_init_depvar(M,1,Y)

: moptimize_init_eq_indepvars(M,1,X)

:
: moptimize_init_userinfo(M, 1, (X) )

: moptimize(M)
initial:       f(p) =  2.024e+13
alternative:   f(p) =  2.024e+13
rescale:       f(p) =  1.608e+12
Iteration 0:   f(p) =  1.608e+12  
Iteration 1:   f(p) =  5.449e-11  
Iteration 2:   f(p) =  2.921e-16  

: H = moptimize_result_Hessian(M)

: b = moptimize_result_coefs(M)

: N = rows(Y)

:
: sigma = quadcolsum( (Y-(X,J(N,1,1))*b'):^2)/(N-1)

: qrinv(quadcross(X, 1,X, 1))
                  1              2
    +-------------------------------+
  1 |   .0004092558   -.0087160428  |
  2 |  -.0087160428    .1991416689  |
    +-------------------------------+

: V = qrinv(quadcross(X, 1,X, 1))*sigma

: b',sqrt(diagonal(V))
                  1              2
    +-------------------------------+
  1 |  -108.4315547    9.281294329  |
  2 |   5328.758517    204.7350436  |
    +-------------------------------+

:
:
: qrinv(H)
                  1              2
    +-------------------------------+
  1 |   .0000760945   -.0017403548  |
  2 |  -.0017403548    .0398039838  |
    +-------------------------------+

: b
                  1              2
    +-------------------------------+
  1 |  -108.4315547    5328.758517  |
    +-------------------------------+

: H
[symmetric]
                 1             2
    +-----------------------------+
  1 |   1299565119                |
  2 |  56821055.74   2484419.577  |
    +-----------------------------+

:
: moptimize_result_V(M)
[symmetric]
       1   2
    +---------+
  1 |  .      |
  2 |  .   .  |
    +---------+

: moptimize_result_display(M)

                                                  Number of obs   =         74

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |  -108.4316          .        .       .            .           .
       _cons |   5328.759          .        .       .            .           .
------------------------------------------------------------------------------

:
:                
: end
-------------------------------------------------------------------------------------------------------------------------------------------

.
.
. reg weight mpg

      Source |       SS       df       MS              Number of obs =      74
-------------+------------------------------           F(  1,    72) =  134.62
       Model |  28728735.3     1  28728735.3           Prob > F      =  0.0000
    Residual |  15365443.1    72  213408.932           R-squared     =  0.6515
-------------+------------------------------           Adj R-squared =  0.6467
       Total |  44094178.4    73  604029.841           Root MSE      =  461.96

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -108.4316   9.345526   -11.60   0.000    -127.0615   -89.80159
       _cons |   5328.759   206.1519    25.85   0.000     4917.802    5739.715
------------------------------------------------------------------------------

sysuse auto, clear

clear mata
mata:

void Reg(M, todo, b, r, S)
        {
        
                real scalar phi1
                real scalar phi2
                real colvector w, xb, zg, e, iota, t, u
                real scalar N
                real matrix W, Z
                
                
                xb = moptimize_util_xb(M,b,1)
                N = rows(xb)

                
                u = moptimize_util_depvar(M,1) :- xb
                W = moptimize_util_userinfo(M, 1)
                
                phi1 = quadcross(W, 1, u, 0)
                r = (phi1)            
        
        
        }

Y = st_data(.,"weight")
X = st_data(.,"mpg")
M = moptimize_init()
moptimize_init_evaluatortype(M,"q0")
moptimize_init_evaluator(M,&Reg())
moptimize_init_technique(M,"gn")
moptimize_init_depvar(M,1,Y)
moptimize_init_eq_indepvars(M,1,X)

moptimize_init_userinfo(M, 1, (X) )
moptimize(M)
H = moptimize_result_Hessian(M)
b = moptimize_result_coefs(M)
N = rows(Y)

sigma = quadcolsum( (Y-(X,J(N,1,1))*b'):^2)/(N-1)
qrinv(quadcross(X, 1,X, 1))
V = qrinv(quadcross(X, 1,X, 1))*sigma
b',sqrt(diagonal(V))


qrinv(H)
b
H

moptimize_result_V(M)
moptimize_result_display(M)

        
end


reg weight mpg

[Jeff Pitblado (StataCorp)]

I think the problem is in the evaluator. In the documentation for moptimize() the type q evaluators are supposed to return an r vector and optionally an S matrix. r is supposed to be a vector of indepenent elements, think residuals that define the moment conditions. S is supposed to be matrix of scores.

For a linear regression, the elements are really just residuals. Consider the following q0 evaluator.

void myreg(M, todo, beta, r, S)
{
        r = moptimize_util_depvar(M) :- moptimize_util_xb(M, beta, 1)
}

Here is Christophe's example using this evaluator:

: O = moptimize_init()

: moptimize_init_evaluatortype(O,"q0")

: moptimize_init_evaluator(O,&myreg())

: moptimize_init_valueid(O,"Residual SS")

: moptimize_init_technique(O,"gn")

: moptimize_init_depvar(O,1,Y)

: moptimize_init_eq_indepvars(O,1,X)

: moptimize(O)
initial:       Residual SS =  7.188e+08
alternative:   Residual SS =  7.185e+08
rescale:       Residual SS =  1.139e+08
Iteration 0:   Residual SS =  1.139e+08  
Iteration 1:   Residual SS =   15365443  
Iteration 2:   Residual SS =   15365443  (backed up)
Iteration 3:   Residual SS =   15365443  (backed up)
Iteration 4:   Residual SS =   15365443  (backed up)

: moptimize_result_display(O)

                                                Number of obs     =         74

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |  -108.4316   9.345683   -11.60   0.000    -126.7488   -90.11435
       _cons |   5328.759   206.1657    25.85   0.000     4924.681    5732.836
------------------------------------------------------------------------------

The (backed up) messages are caused by a combination of the scale of weight and numerical derivatives. Using a currently undocument function of moptimize(), the q1 version of the above evaluator would be:

void myreg1(M, todo, beta, r, S)
{
        r = moptimize_util_depvar(M) :- moptimize_util_xb(M, beta, 1)
        if (todo) { 
                S = -moptimize_util_indepvars(M, 1)
        }
}

Here is the result using this evaluator:

: P = moptimize_init()

: moptimize_init_evaluatortype(P,"q1")

: moptimize_init_evaluator(P,&myreg1())

: moptimize_init_valueid(P,"Residual SS")

: moptimize_init_technique(P,"gn")

: moptimize_init_depvar(P,1,Y)

: moptimize_init_eq_indepvars(P,1,X)

: moptimize(P)
initial:       Residual SS =  7.188e+08
alternative:   Residual SS =  7.185e+08
rescale:       Residual SS =  1.139e+08
Iteration 0:   Residual SS =  1.139e+08  
Iteration 1:   Residual SS =   15365443  
Iteration 2:   Residual SS =   15365443  

: moptimize_result_display(P)

                                                Number of obs     =         74

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |  -108.4316   9.345526   -11.60   0.000    -126.7484   -90.11466
       _cons |   5328.759   206.1519    25.85   0.000     4924.708    5732.809
------------------------------------------------------------------------------

[Christophe Kolodziejczyk]

Dear Jeff
Thanks for your reply and comments. I wanted to write a longer post and discuss these issues in more detail, but I couldn't find the time.

I took the OLS as an example of a GMM. Minimizing the sum of squared residuals should give the same results compared to making the residuals orthogonals to the Xs, but in my opinion these are two different approaches. The way the Gauss-Newton algorithm is setup in moptimize() it will also be difficult to implement instrumental variables estimators for non-linear models or systems of non-linear models. But it is likely I am mistaken.
So far, my conclusion from your post is that maybe moptimize() with the gauss-netwon algortihm is not well suited for GMM estimation and optimize() would be a better solution.

More documentation and examples on this algorithm would also be welcome.

p.s. which undocumented function of moptimize() are you talking about?

[Jeff Pitblado (StataCorp)]

In myreg1(), I reference moptimize_util_indepvars() which is not a documented function in [M-5] moptimize().

[Christophe Kolodziejczyk]

This is an old post which I started years ago. I am using moptimize with a q evaluator and the Gauss-Newton algorihtm (gn). I wanted to compare my analytical gradient with the numerical one, but it seems that moptimize with q1debug and the option for tracing the gradient is not working properly.

sysuse auto, clear

clear mata
mata:
Y = st_data(.,"weight")
X = st_data(.,"mpg")

void myreg(M, todo, beta, r, S)
{
        r = moptimize_util_depvar(M) :- moptimize_util_xb(M, beta, 1)
}

void myreg1(M, todo, beta, r, S)
{
        r = moptimize_util_depvar(M) :- moptimize_util_xb(M, beta, 1)
        if (todo) { 
                S = -moptimize_util_indepvars(M, 1)
        }
}

P = moptimize_init()

moptimize_init_evaluatortype(P,"q1debug")

moptimize_init_evaluator(P,&myreg1())

moptimize_init_valueid(P,"Residual SS")

moptimize_init_technique(P,"gn")

moptimize_init_depvar(P,1,Y)

moptimize_init_eq_indepvars(P,1,X)
moptimize_init_trace_coefs(P, "on")
moptimize_init_trace_gradient(P, "on")

moptimize(P)
moptimize_result_display(P)

The previous code, where the trace for the gradient is on, will give the following log. The comparison of the two gradients is not printed by moptimize, where it will be the case with other techniques (nr).

: moptimize(P)
initial:       Residual SS =  7.188e+08
trying nonzero initial values ++
alternative:   Residual SS =  7.185e+08
rescaling entire vector .++++++++++++.
rescale:       Residual SS =  1.139e+08
------------------------------------------------------------------------------------------------------------------------------------
Iteration 0:
Parameter vector:
      eq1:   eq1:
       x1  _cons
r1      0   2048

q1debug:  Begin derivative-comparison report ---------------------------------------------------------------------------------------
q1debug:  mreldif(gradient vectors) =         0
q1debug:  End derivative-comparison report -----------------------------------------------------------------------------------------

                                                      Residual SS =  1.139e+08
------------------------------------------------------------------------------------------------------------------------------------
Iteration 1:
Parameter vector:
          eq1:       eq1:
           x1      _cons
r1  -108.4316   5328.759

q1debug:  Begin derivative-comparison report ---------------------------------------------------------------------------------------
q1debug:  mreldif(gradient vectors) =  .0000508
q1debug:  End derivative-comparison report -----------------------------------------------------------------------------------------

                                                      Residual SS =   15365443
------------------------------------------------------------------------------------------------------------------------------------