/* This program is for the generalized method of moments (GMM) procedure.
See "Generalized Instrumental Variable Estimation of Nonlinear Rational
Expectations Models," by L.P. Hansen and K.J. Singleton (Econometrica
1982, Vol.50, 1269-1286).  This example file is for Hansen and Singleton
model, but uses different data than theirs.

Model:   E[h(t,b0)]=0:  nmr moments restrictions, h(t,b0) is nmr by 1

Algorithm:
   This program calculates and iterates the following:
     ititial W0 given.

  1. g(b)=C*[h(1,b)+,...,+h(T,b)]
     b=argmin{g(t,b)'W0*g(t,b)}

  2. Rzw(j)=(1/T)[h(1,b)h(1,b)'+,..,+h(T,b)h(T-j,b)']
                       j=0,1,..,mas


*/
new ;                               @clear workspace and program space@
Library optmum;
#include optmum.ext ;
optset ; @reset OPTMUM globals @
__design = 0;

          @---choose weighting matrix estimation method------@
calwflag = 1;     @  1: Newey-West  2: quadr spctrl  3: c:\gauss\w0  @

pw=2 ;        @-----choice of prewhitening of weightmatrix----@
if pw==1 ; "VAR(1) prewhitening " ;
elseif pw==2; "individual prewhitening " ;
else ; "no prewhitening" ;
endif ;

outfile= "hansin.out";
output file = ^outfile reset;
? "outfile: ";; outfile;
output on;
datestr(date) ;;" " ;;timestr(time) ;

 @ initialization of constants and matrices @
clear alpq1,alpq2,S_Tnw,S_Tqs,w0,tsm,n,nmom,dsm,vk,prob;
clear noconv,test,pd_err,nwtest ;

@ ----------- READING IN DATA --------------------------- @
tt=334 ;   @no of obs@
load x[335,3]=d:\teach\hansin.dat;
c=x[1:tt+1,1];     @ c(t)/c(t-1) @
re=x[1:tt+1,2];    @ vwr @
rf=x[1:tt+1,3];    @ T-Bill rate @
clear x,nwtest;

  @------------------- Constructing instrumental vars --------@

  zp=ones(tt,1)~c[1:tt,.]~re[1:tt,.];


S_Tnw = .5*TT^(1/3) ;    @number of lags in Newey-West weight scheme@
maxnw = 10*tt^(1/3);

@-------------- initial param values --------------------@
bgm=1|1 ;

nmom=3 ;              @------no moment restrictions-------@

w0= eye(nmom) ;       @initial estimation uses id weight  @

"" ;
"** GMM Results * first round ";
     _opmiter = 50; @max 50 iter for first round @
call gmmq;

"" ;
"** initiating second round ** ";
"";
w0 = wocalc(calwflag);
 if nwtest ==1 ;
   pd_err = 1 ;  ndpclex ;
   goto pdcount ;
 endif ;

   if isinf(w0) ==1 ;
   pd_err = 1 ; ndpclex ;
   goto pdcount ;
endif ;


     _opmiter = 200; @max 200 iter for 2nd round @
 call gmmq;

 if noconv>0 ;
    "non convergence - return code = " ;; noconv ;
    "this estimation discarded" ;
    n_noconv = n_noconv+1 ;
    repeat = 1 ;
    goto endloop ; @mostly so that failure is only counted once@
  endif ;


pdcount:
   if pd_err == 1;
        ndpclex ;
        "non-pd error estimation discarded" ;
        n_pderr = n_pderr + 1 ;
        repeat = 1 ;
        goto endloop ;
   endif ;

endloop:

"" ;

@ ------ END OF PROGRAM ----------------- @

datestr(date) ;;" " ;;timestr(time) ;
output off;

proc (0)=gmmq;
    local tm,gort,df,chi,vbgm,g,h,vof,qd,
          oldval,__btol,__gtol,np,sd;

? "Initial values =" bgm';
format /m1 12,7;
df=nmom-rows(bgm);

    vof = ofn(bgm);
    ? "initial vof=";; vof;
/*    _opditer = 100 ;                @maximum number of iterations@ */
     __btol = .0000001;               @convergence tolerance @
     __gtol = .0000001;
     _opstep = 4 ;                   @4 is brent - the manual is wrong here@
     _opmbkst = 30 ;                 @max number of backsteps @
/*  _opgdprc=&derivat ;            @analytic deriv if you have them@    */
    output off;              @or output file will fill w garbage@
    {bgm,vof,g,h} = optmum(&ofn,bgm) ;
    output on;
    noconv= h ;
    np = rows(bgm);
     ? " parameter estimate";; bgm';

     gort=gradcd(&ort,bgm);

    if isindef(gort) == 1 ;
      pd_err = 1 ; ndpclex ;
      goto pdout ;
    endif ;
    if isinf(gort) == 1;
      "infty reached" ; pause(1) ;
      pd_err = 1 ; ndpclex ;
      goto pdout ;
    endif ;

/* the following method to trap the non-pd error
   has been copies from the GAUSS command reference manual
   see un TRAP */
oldval=trapchk(1) ;
trap 1,1 ;
    vbgm=invpd(gort'w0*gort);
trap oldval,1 ;
if scalerr(vbgm) ;
   pd_err = 1 ; ndpclex ;
   goto pdout ;
else ;
   pd_err = 0 ;
endif ;

    sd=sqrt(diag(vbgm)/tt);
     ? "   s.d.=";; sd';
   @  ? "obs:"; qd; @
   chi=vof*tt ;
   prob=cdfchic(chi,df);
     ? "vof:" ;; vof ;;  ? " chi square=" chi "(" prob ")";
pdout:
endp;            @ this is the end of the GMMQ routine @

                      @returns the time average of the ortho cond@
proc ort(bgm) ;
local zw1 ;
zw1=ortho(bgm) ;
retp(sumc(zw1)/tt) ;
endp ;

proc ortho(b);        @*** orthogonality condition **@
 local ps,w1,w2,zw1,zw2;

   ps=b[1].*(c[2:tt+1,.]^(-b[2]));
   w1=ps.*re[2:tt+1,.]-1;
   zw1=zp.*w1;

 retp(zw1);
endp;


proc ofn(b);              @**** quadr object fuct *@
  local np,fval,ortg;
  np = rows(b);
  ortg=ort(b);
  fval = ortg'w0*ortg;
retp(fval);
endp;



proc(1) = wocalc(calwflag);
    local w,nm;
if calwflag == 1 ;
  w=newest ;
elseif calwflag == 2 ;
    w=qspec;
elseif calwflag == 3;
    loadm w = c:\gauss\w06;
elseif calwflag == 4;
    nm = rows(nmom);
    w = eye(nm);
endif;
retp(w);
endp;


proc(1)= newest;                      @ NEWEY-WEST W0 **@
@"calculating NW weight matrix" ;@
local tau,omega22,ometau,w,gikgodt,oldval,wzp1,wzp,AA;

nwtest=0 ;
wzp1 = ortho(bgm) ; ;
if pw == 1 ;
   {wzp,AA} = prewhit1(wzp1) ;
 @  "VAR(1) pre-whitening" ;@
elseif pw == 2 ;
   {wzp,AA} = prewhit2(wzp1) ;
   "individual AR(1) pre-whitening" ;
else ;
  wzp = wzp1 ;
endif ;
omega22=wzp'wzp;

tau=1;
call alcalcpw(wzp) ;
do until tau > (S_Tnw-1);
 /*  "doing lag no :"; tau; */
   ometau = wzp[tau+1:tt-1,.]'wzp[1:tt-1-tau,.];
   omega22 = omega22 + (1 - tau/S_Tnw) * (ometau+ometau');
   tau=tau+1;
endo;

omega22=omega22/(tt-1);

if (pw ==1) .or (pw==2) ;
@  "re-color" ;  pause(1) ;@
  omega22= recolor(omega22,AA) ;
endif ;

/* the following method to trap the non-pd error
   has been copies from the GAUSS command reference manual
   see un TRAP */
oldval=trapchk(1) ;
trap 1,1 ;
    w=invpd(omega22);
trap oldval,1 ;
if scalerr(w) ;
   "Problem inverting W" ; ndpclex ;   nwtest=1 ;
else ;
   "Inversion of W OK" ;
endif ;

@"NW done " ;@
retp(w) ;
endp ;

proc(0) = alcalcpw(wzp) ;
local roa,sigma2a,resa,alpdenom;
@"calculating Andrews alpha " ;@
roa = diag(wzp[2:tt-1,.]'wzp[1:tt-2,.]) ;
roa = roa./diag(wzp'wzp) ;
resa = wzp[2:tt-1,.] - wzp[1:tt-2,.]*roa ;
sigma2a = diag(resa'resa)/(tt-2) ;
alpdenom = sumc(sigma2a^2./(1-roa)^4) ;
alpq1 = sumc( 4*(roa^2).*(sigma2a^2)./(((1-roa)^6).*((1+roa)^2)) )/alpdenom ;
alpq2 = sumc( 4*(roa^2).*(sigma2a^2)./((1-roa)^8) )/alpdenom ;
S_Tnw = round( 1.1447*(alpq1*(tt-1))^(1/3) ) ; "S_Tnw" ;; S_Tnw ;
S_Tqs = round( 1.3221*(alpq2*(tt-1))^(1/5) ) ; "S_Tqs" ; S_Tqs ;
if ((S_Tnw > maxnw) OR (S_Tnw==0vfff8000000000000 )) ;
    S_Tnw = maxnw;
    "S_Tnw forced to maxnw" ;
endif;
if S_Tqs > 25;
    S_Tqs = 25;
endif;
endp;


proc kx(x) ;
local kx ;
kx = 25 * (sin(6*pi*x/5)/(6*pi*x/5)-cos(6*pi*x/5)) / (12*pi^2*x^2) ;
"kx " ; kx ;
retp(kx) ;
endp ;

proc qspec ;                     @**** QS W0 **@
"calculating QS weight matrix" ;
 local omega22,ometau,tau,w,tsm,wzp,aa;
local oldval;
tsm=ortho(bgm) ;
tau=1;
if pw == 1 ;
   {wzp,AA} = prewhit1(tsm) ;
 @  "VAR(1) pre-whitening" ;@
elseif pw == 2 ;
   {wzp,AA} = prewhit2(tsm) ;
   "individual AR(1) pre-whitening" ;
else ;
  wzp = tsm ;
endif ;
omega22=wzp'wzp;
call alcalcpw(wzp) ;
do until tau > 3 *  S_Tqs ;
  ometau = wzp[tau+1:tt,.]'wzp[1:tt-tau,.];
  omega22=omega22 + kx(tau/S_Tqs)*(ometau+ometau');
   tau=tau+1;
endo;
omega22=omega22/tt;
w = invpd(omega22);


if (pw ==1) .or (pw==2) ;
@  "re-color" ;  pause(1) ;@
  omega22= recolor(omega22,AA) ;
endif ;

/* the following method to trap the non-pd error
   has been copies from the GAUSS command reference manual
   see un TRAP */
oldval=trapchk(1) ;
trap 1,1 ;
    w=invpd(omega22);
trap oldval,1 ;
if scalerr(w) ;
   "Problem inverting W" ; ndpclex ;   nwtest=1 ;
else ;
   "Inversion of W OK" ;
endif ;

"QS done " ;

retp(w) ;
endp ;


proc(2) = prewhit1(wzp) ;
local aa ;
aa = wzp[2:tt,.]'wzp[1:tt-1,.]/(wzp[1:tt-1,.]'wzp[1:tt-1,.]) ;
"singular values for aa:  " ;; svd(aa) ;
pause(2) ;
@Note: parametrization of aa corresponds to y_t = aa*y_{t-1}, where
y_t is k times 1  as in Andrews and Monahan '92                      @

retp(wzp[2:tt,.] - wzp[1:tt-1,.]*aa',aa) ;
endp ;

proc(2) = prewhit2(wzp) ;
local aa,aaa,ii ;
aaa=ones(nmom,1) ;
ii = 1 ;
do until ii> nmom ;
 aaa[ii] = wzp[2:tt,ii]'wzp[1:tt-1,ii]/(wzp[1:tt-1,ii]'wzp[1:tt-1,ii]) ;
 ii = ii+1 ;
endo ;
   if (maxc(abs(aaa)) .> .99)
       OR (maxc( abs(aaa)).==0vfff8000000000000 )==1 ;
      aaa=zeros(nmom,1) ;
      "wild values for prewhitening - prewhitening abandoned " ;
   endif ;
@"values for aaa:  " ;; aaa' ;@
aa = diagrv(zeros(nmom,nmom),aaa) ;
pause(1) ;
@Note: parametrization of aa corresponds to y_t = aa*y_{t-1}, where
y_t is k times 1  as in Andrews and Monahan '92                      @

retp(wzp[2:tt,.] - wzp[1:tt-1,.]*aa',aa) ;
endp ;

proc recolor(w,aa) ;
local don_d ;
don_d = inv(eye(nmom) - aa) ;
retp(don_d*w*don_d') ;
endp ;

   @   -------------   PROC  ISINDEF   ---------------   @

proc isindef(x);
    local indef;
    indef = 0vfff8000000000000 ;
    retp(not x $/= indef);
endp;


   @   -------------   PROC  ISINF   ---------------   @

proc isinf(x);
    local plus,minus;
    plus = 0v7ff0000000000000 ;
    minus = 0vfff0000000000000 ;
    retp(not x $/= plus or not x $/= minus) ;
endp;