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

The program is the actual program written by Hansen's RA's (Heaton
and Ogaki), I believe. New ways of estimating variance matrices
have been added, as well as ways of "capturing" crashes */

/* The main structure of the program is as follows:

a) you set variables that guides choices of weighting matrix and
how it is calculated

b) you read in the data

c) you choose intial values and initial weighting matrix (typically
the identity matrix)

d) you call gmmq procedure, which estimates the parameters

e) you calculate the weighting matrix using the `residuals' from
the first calculation

f) you estimate the model again

The gmmq procedure calls the optimization procedure "feeding" it
the objective function and after optimization calculates std errors
using the numerical gradient

The objective function is the procedure ofn

The ofn procedure calls the procedure ort which is the time average
of the period t "residuals"

the procedure ort calls the procedure ortho, which contains the
economic residuals multiplied by period t instruments---the reason
for having these routines is that the period-by-period output is
needed for calculating the variance matrix, even if only time
averaged output from ort goes into the objective function

the procedure newest calculates the variance matrix --- in the
process it calls alcalcpw which calculates the optimal "bandwidth"
S_Tnw (an alternative procedure for this is qspec, but understand
the logic of newest first

***************
 Note: This program is set up so that it can be used
in a Monte Carlo study, this obscure subroutine isinf check for
number that has gone to infinity and would make the program crash
if I don't trap it. I also check if the weighting matrix is
positive definite, but you can safely skip all that if you
translate the program to Matlab.

*/


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=0 ;

@-----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]=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;