/*
** GMM FOR SINGLE EQUATION
** Using Newey-West and Andrews Covariance Matrix
*/

new ;

@ LOADING DATA @

    load dat[923,17] = mwemp.db ; @ Data on employed married women @

@ OPEN OUTPUT FILE @

    output file = gmm_1.out reset ;

@ DEFINE VARIABLES @

    lrate   = dat[.,1]  ;
    edu     = dat[.,2]  ;
    urb     = dat[.,3]  ;
    minor   = dat[.,4]  ;
    age     = dat[.,5]  ;
    tenure  = dat[.,6]  ;
    expp    = dat[.,7]  ;
    regs    = dat[.,8]  ;
    occw    = dat[.,9]  ;
    occb    = dat[.,10] ;
    indumg  = dat[.,11] ;
    indumn  = dat[.,12] ;
    unionn  = dat[.,13] ;
    unempr  = dat[.,14] ;
    lofinc  = dat[.,15] ;
    hwork   = dat[.,16] ;
    kids5   = dat[.,17] ;
    lhwork  = ln(hwork+1) ;

@ Define # of observations @

   n = rows(dat) ;

@  DEFINE EQUATION @

proc eqq(b)  ;
local equ    ;

   equ = lhwork - b[1]*ones(n,1) - b[2]*lrate ;

retp(equ) ;
endp      ;

@   DEFINE INSTRUMENTAL VARIABLES  @

    iv = ones(n,1)~age~expp~edu;

@   DEFINE INITIAL VALUE @

    bb = {1,1} ;

@   Choose np = 1 for Newey-West;
           np = 0 for Andrews        @

    np = 1 ;

@   Specify Bandwidth @

    bandw = 0 ; @ if bandw = 0, White-Correction @

@   iter = 0 for two-step GMM ; iter = 1 for iterative GMM @

    iter = 1 ;

/*
**  DO NOT CHANGE FROM HERE
*/

@   Procedure for Newey-West   @

proc nw(uv,l) ;
local j,omeo,n,w,ome  ;

    n    = rows(uv)  ;
    omeo = uv'uv/n   ;

    j = 1; do while j <= l  ;

        w    = 1 - j/(l+1)       ;
        ome  = uv[j+1:n,.]'uv[1:n-j,.]/n ;
        omeo = omeo + w*(ome+ome')       ;

    j = j + 1 ; endo ;

retp(omeo) ;
endp       ;

@   Procedure for Andrews (1991) @

proc quadk(z) ;
local rr ;

  rr = ( 3/(6*pi*z/5)^2 )
        *( 5*sin(6*pi*z/5)/(6*pi*z) - cos(6*pi*z/5) ) ;

retp(rr) ;
endp     ;

proc andrews(uv,l) ;
local j,omeo,n,w,ome  ;

    n    = rows(uv)  ;
    omeo = uv'uv/n   ;

    j = 1; do while j <= n-1 ;

        w    = j/(l+1)                    ;
        ome  = uv[j+1:n,.]'uv[1:n-j,.]/n  ;
        omeo = omeo + quadk(w)*(ome+ome') ;

    j = j + 1 ; endo ;

retp((n/(n-rows(bb)))*omeo) ;
endp       ;


library optmum;
#include optmum.ext;
optset ;

/* FIRST ROUND */

@   Creating weighting matrix @

    w = invpd(iv'iv) ;

proc  f(b) ;
local m    ;

    m = iv'eqq(b) ;

RETP(  m'w*m  ) ;
ENDP  ;

@ starting values @

b0 = bb  ;

__title = "FIRST-STAGE GMM ";

_opgtol = 1e-4;
_opstmth = "BFGS HALF";
__output = 0 ;

{b,func,grad,retcode} = optprt(optmum(&f,b0)) ;

initb = b ;

/* SECOND ROUND */

cre = 10 ; do while cre >= 0.1e-4 ;
initb = initb ;

@ CREATING WEIGHTING MATRICES @

    if np == 1 ; gosub weight1 ;
    else       ; gosub weight2 ;
    endif      ;

library optmum;
#include optmum.ext;
optset ;

/* Second Round GMM */

proc  f1(b) ;
local m    ;

    m = iv'eqq(b) ;

RETP(  m'w*m ) ;
ENDP  ;

@ starting values @

B0 = initb ;

__title = "SECOND-STAGE GMM ";
_opgtol = 1e-4;
_opstmth = "BFGS HALF";
__output = 0 ;

{b,func,grad,retcode} = optprt(optmum(&f1,b0)) ;

cre = iter*(initb-b)'(initb-b) ;
initb = b ;
endo      ;


/* CONPUTING STANDARD ERRORS */

format /rd 12,4 ;

proc mom(b)     ;
retp(iv'eqq(b)) ;
endp            ;

gra = gradp(&mom,b) ;
cov = invpd(gra'w*gra) ;
se  = sqrt(diag(cov))  ;
tst = b./se  ;

"" ;
"GMM RESULTS" ;
"" ;
"      coeff.      std. err.      t-st " ;
b~se~tst     ;

/*  HANSEN TEST */

df = cols(iv) - rows(b) ;
"" ;
"J TEST, DF, P-Val =" func  df  cdfchic(func,df) ;

end ; @ This means no longer using subroutines  @

output off ;

weight1:
    w  = invpd( n*nw(iv.*eqq(initb),bandw) ) ;
    return ;

weight2:
    w  = invpd( n*andrews(iv.*eqq(initb),bandw) ) ;
    return ;