/*
** GMM FOR THTREE EQUATIONS
** Using Newey-West or Andrews Covariance Matrix
*/

NEW ;

@ LOADING DATA @

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

@ OPEN OUTPUT FILE @

    output file = gmm_3.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 in GMM  @

    n = rows(dat) ;

@  DEFINE EQUATIONS @

	proc eqq1(b) ;
	local equ1   ;

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

	retp(equ1) ;
	endp       ;

	proc eqq2(b) ;
	local equ2   ;

	   equ2 = lhwork-b[4]*lrate-b[5]*ones(n,1)-b[6]*kids5 ;

	retp(equ2) ;
	endp       ;

	proc eqq3(b) ;
	local equ3   ;

	   equ3 = lofinc-b[7]*lrate-b[8]*lhwork - b[9]*unempr  ;

	retp(equ3) ;
	endp       ;

@   DEFINE INSTRUMENTAL VARIABLES FOR EACH EQUATION  @

    iv1 = ones(n,1)~tenure~kids5~unempr ;
    iv2 = ones(n,1)~tenure~kids5~unempr ;
    iv3 = ones(n,1)~tenure~kids5~unempr ;

@   DEFINE INITIAL VALUE @

    bb = { 1, 1, 1, 1, 1, 1, 1, 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 = 0 ;


/*
**  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(omeo) ;
endp       ;


library optmum;
#include optmum.ext;
OPTSET ;

/* FIRST ROUND */

@   Creating weighting matrix @

    k1 = cols(iv1) ;
    k2 = cols(iv2) ;
    k3 = cols(iv3) ;

    w1 = (iv1'iv1)~zeros(k1,k2)~zeros(k1,k3) ;
    w2 = zeros(k2,k1)~(iv2'iv2)~zeros(k2,k3) ;
    w3 = zeros(k3,k1)~zeros(k3,k2)~(iv3'iv3) ;
    w  = invpd(w1|w2|w3) ;

proc  f(b) ;
local m    ;

    m = (iv1'eqq1(b))|(iv2'eqq2(b))|(iv3'eqq3(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 = (iv1'eqq1(b))|(iv2'eqq2(b))|(iv3'eqq3(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( (iv1'eqq1(b))|(iv2'eqq2(b))|(iv3'eqq3(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(iv1~iv2~iv3) - rows(b) ;
"" ;
"J TEST, DF, P-Val =" func  df  cdfchic(func,df) ;

end ; @ This means no longer using subroutines  @

output off ;

weight1:
    wee = (iv1.*eqq1(initb))~(iv2.*eqq2(initb)) ;
    wee = wee~(iv3.*eqq3(initb) ) ;
    w   = invpd( rows(iv1)*nw(wee,bandw) )      ;
    return ;

weight2:
    wee = (iv1.*eqq1(initb))~(iv2.*eqq2(initb))  ;
    wee = wee~(iv3.*eqq3(initb) ) ;
    w   = invpd( rows(iv1)*andrews(wee,bandw) )  ;
    return ;