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

new ;

@ LOADING DATA @

    load dat[201,12] = econmac.txt ; 
    @ Data on US macro data: 1953:2 - 2003:2 @

@ DEFINE VARIABLES @

    gdp    = dat[.,1]  ;
    rgdp   = dat[.,2]  ;
    pce    = dat[.,3]  ;
    rpce   = dat[.,4]  ;
    wage   = dat[.,5]  ;
    cpi    = dat[.,6]  ;
    profit = dat[.,7]  ;
    rexpert= dat[.,8]  ;
    rimport= dat[.,9]  ;
    fgdep  = dat[.,10] ;
    tb1    = dat[.,11] ;
    tb2    = dat[.,12] ;


@ OPEN OUTPUT FILE @

    output file = gmm_nw_t.out reset ;

@ Define lags to use @

    lag = 9 ;
    
@ Define # of observations @

   nn = rows(dat) ;
   n  = nn - lag  ;

@ Define variables to use @

    cpi_0 = cpi[10:nn,.];
    cpi_1 = cpi[9:nn-1,.];
    cpi_2 = cpi[8:nn-2,.];
    cpi_3 = cpi[7:nn-3,.];
    cpi_4 = cpi[6:nn-4,.];
    cpi_5 = cpi[5:nn-5,.];
    cpi_6 = cpi[4:nn-6,.];
    cpi_7 = cpi[3:nn-7,.];
    cpi_8 = cpi[2:nn-8,.];
    cpi_9 = cpi[1:nn-9,.];
    
    inf4_0 = ln(cpi_0./cpi_4);
    inf4_4 = ln(cpi_4./cpi_8);
    inf4_5 = ln(cpi_5./cpi_9);

    gret_0  = 1 + tb1[10:nn,.]/100 ;
    gret_1  = 1 + tb1[9:nn-1,.]/100 ;
    gret_2  = 1 + tb1[8:nn-2,.]/100 ;
    gret_3  = 1 + tb1[7:nn-3,.]/100 ;
    gret_4  = 1 + tb1[6:nn-4,.]/100 ;
    gret_5  = 1 + tb1[5:nn-5,.]/100 ;
    gret_6  = 1 + tb1[4:nn-6,.]/100 ;
    gret_7  = 1 + tb1[3:nn-7,.]/100 ;
    gret_8  = 1 + tb1[2:nn-8,.]/100 ;
    gret_9  = 1 + tb1[1:nn-9,.]/100 ;
 
    rgret_0 = gret_0 - inf4_0 ;
    rgret_4 = gret_4 - inf4_4 ;
    rgret_5 = gret_5 - inf4_5 ;
        

    rpce_0  = rpce[10:nn,.] ;
    rpce_1  = rpce[9:nn-1,.] ;
    rpce_2  = rpce[8:nn-2,.] ;
    rpce_3  = rpce[7:nn-3,.] ;
    rpce_4  = rpce[6:nn-4,.] ;
    rpce_5  = rpce[5:nn-5,.] ;
    rpce_6  = rpce[4:nn-6,.] ;
    rpce_7  = rpce[3:nn-7,.] ;
    rpce_8  = rpce[2:nn-8,.] ;
    rpce_9  = rpce[1:nn-9,.] ;
    
    rgc_0   = rpce_0./rpce_4;
    rgc_4   = rpce_4./rpce_8;
    rgc_5   = rpce_5./rpce_9;

@  DEFINE EQUATION @

proc eqq(b)  ;
local equ    ;

   equ = b[1]*(rgret_0).*rgc_0^(-b[2]) - ones(n,1)  ;

retp(equ) ;
endp      ;

  vnp = {"discp", "riskp" } ;


@   DEFINE INSTRUMENTAL VARIABLES  @

    iv = ones(n,1)~rgret_4~rgret_5~rgc_4~rgc_5 ;

@   DEFINE INITIAL VALUE @

    bb = {1,7} ;

@   Specify Bandwidth @

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


/*
**  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  ;

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

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

initb = b ;

/* SECOND ROUND */

gg = iv.*eqq(initb);
w  = invpd(n*nw(gg,bandw));

proc  f1(b) ;
local m    ;

    m = iv'eqq(b) ;

RETP(  m'w*m ) ;
ENDP  ;

@ starting values @

B0 = initb ;

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

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


/* 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  ;

econ = vnp~b~se~(b./se);

let mask[1,4] = 0 1 1 1;
let fmt[4,3] =
    "-*.*s" 8 8
    "*.*lf" 10 4
    "*.*lf" 10 4
    "*.*lf" 10 4;

format /rd 10,4 ;

"" ;
"GMM Estimation Results" ;
"------------------------" ;

"para         coeff.  std. err.   t-st " ;
yyprin = printfm(econ,mask,fmt); 
"" ;

/*  HANSEN TEST */

df = cols(iv) - rows(b) ;

"J TEST, DF, P-Val =" func  df  cdfchic(func,df) ;
"# of observations =" n ;

output off ;