@ GMM.SET @
@------------------------------------------------------------------------------
Prepared by Lars P. Hansen, Masao Ogaki, and John C. Heaton
    Financial support from the National Science Foundation
                 Grant numbers: SES-8512371 and SES-9213930

 Last Revision:  5/18/93
------------------------------------------------------------------------------@
/*
This program has been used and seem to be free of errors. However,
we (Lars P. Hansen, John C. Heaton, and Masao Ogaki) do not assume
responsibility for any remaining errors.  In no event shall we be liable to
for any damages whatsoever arising out ot the use of or inability to use
this program.

   @@ This program sets procedure GMM.
    GMM.EXP is an example file.
*/
#include "partinv.prc";
#include "ywp.prc";
@ Globals @
clear bgm,const,const2,W0,gufn,w0flag,calwflag,nv,ordard;
lend=15;
proc hu(b); @ This is a dummy definition, and not used by GMM programs.@
retp(b);
endp;

proc (1)=ufn(b);
retp(const*sumc(hu(b)));
endp;

proc startval;
retp(bgm);
endp;

proc gmm(&gradq,tend,nmr,mas,kgm,zero,maxitegm,st,wav,maxd,bst);
local gradq:proc;
local nz,df,chi,differ,itergm,vbgm,sd,prob,wzp,omega22,tau,ometau,
      snzi,snzi0,igm,w,vp,prob0,differp,
      bmai,i,ib,ie,rzw,rv,rvwz,ep,tma,xpx,ypx,bma,arcoef,omega2,
      tendm, wvstr,omegas,tenda,als,bls,clsc,va,dells,indm,iwsws,rhoa,
      ea,siga4,alpha2,alpha2d,v2512,v65,xtau,dinv,lamdas,omegat;
clear omega22;
? "Initial values of the coefficients=" bgm';
? "mas=" mas;
format /m1 16,8;
df=nmr-rows(bgm);
tendm=tend-kgm;   @a degree of freedom adjustment@
prob0=100;
clear chi;
differ=1e+12;
differp=differ;
if w0flag==1;
 clear w0;
 ? "initial values of parameters are used to caluculate initial W0";
 itergm=0;
 goto lab100;
elseif w0flag==0;
 w0=eye(nmr);
 ? "Initial W0=I";
 ? "Scaling Constant used for first iteration=" const;
 itergm=1;
 bgm=minquad;
 save bgmi=bgm;
 goto lab100;
endif;
? "W0 in the memory is used as initial W0";
itergm=1;
const=const2;
if w0flag==3;
 vof=ofn(bgm);
? "W0 and bgm in the memory is used for the first GMM result";
 goto lab3200;
endif;
? "W0 in the memory is used for the first GMM result";
do until differ<zero;
     save omega=omega22;
     const=const2;
     bgm=minquad;
     save bgm;
     save w0;
lab3200:
output on;
 gufn=gradq(&ufn,bgm,ufnv,0);
 vbgm=gufn/tend;
 save gammagm=vbgm;   @gamma@
 vbgm=tend*invpd(gufn'w0*gufn./(const^2));  @varb@
 sd=sqrt(diag(vbgm));
 save varb=vbgm;
   @ save variance for bgm [not sqrt(T)*bgm] in VARB.FMT @
? "
================= GMM ITERATION " itergm "====================";
 ? "   b=" bgm';
 ? "   s.e.=" sd';
   if df>0;
     @ S.D. for (1/sqrt(T))*g(b) @
     vbgm=invpd(w0)-(gufn*vbgm*(gufn'))/((const^2)*tend);
     vbgm=sqrt(diag(vbgm));

    chi=(vof)./((const^2).*tend);
    prob=cdfchic(chi,df);
    differp=abs(prob-prob0);
    prob0=prob;
 ? "   chi square=" chi "(" prob ")";
 ? "   d.f.=" df;
 ? "   difference in prob. value from last GMM iteration=" differp;
 ? "(1/sqrt(T))g(b)=" (ufnv')/(const*sqrt(tend));
 ? "           s.e.=" vbgm';
 elseif df==0;? "   Just Identified";
 else;? "   nmr<rows(bgm); Not Identified";endif;
 ? "   scaling const used is" const;
 if itergm>maxitegm-0.5;  @ or differp<zero;@
  goto labe;
 endif;
 lab100:
 wzp=hu(bgm);
 omega22=wzp'wzp;
 rzw=(omega22./tend);@ for Durbin's method @
 tau=1;
 do until tau>mas+0.5;
  ometau=wzp[(1+tau):tend,.]'wzp[1:(tend-tau),.];
  rzw=rzw~(ometau./(tend-tau));
  omega22=omega22+ometau+(ometau');
  tau=tau+1;
 endo;
 omega22=omega22./tendm;
 if calwflag==3;
  goto ldurb;
 elseif calwflag==4;
  goto lnq;
 endif;
 trap 1;
 w=invpd(omega22);
 trap 0;
 if not scalerr(w);
  ? "------- For next GMM iteration
  All the zero restirictions are successfully imposed on W0 -------";
  goto lnexite;
 endif;
 if mas==0;
  ? "!!!!! ERROR !!!!! Since mas=0, OMEGA should be p.s.d.
But OMEGA is not p.d.";
  output off;
  stop;
 endif;
 ? "OMEGA is not positive definite when
 calculated with zero restrictions.";
 if calwflag==0;
  goto ldurb;
 elseif calwflag==1;
  goto lnq;
 endif;

@ Modified Durbin's Method @
   @ 1st Step: Solving Yule-Walker equations to get AR coefficients @
  ldurb:
 ? "------ For next GMM iteration Durbin's Method used for W0 -------";
 ? "Order of AR used=" ordard;
 ? "Order of MA used=" mas;
  arcoef=yw(rzw~zeros(nmr,nmr*(ordard-mas)));
  clear rzw;
  @ 2nd Step: OLS to get MA coefficients @
    @ Get e(t)'s in ep@
  ep=wzp[ordard+1:tend,.];
  tau=1;
  do until tau>ordard;
   ep=ep-wzp[ordard+1-tau:tend-tau,.]*(arcoef[.,(tau-1)*nmr+1:tau*nmr]');
  tau=tau+1;
  endo;
  wzp=wzp[ordard+1:tend,.];
  tma=tend-ordard;

    @ Regress wz(t) on e(t)  @
      @ Get X'X in xpx@
  xpx=ep[mas:tma-1,.];
  tau=2;
  do until tau>mas;
   xpx=xpx~ep[mas+1-tau:tma-tau,.];
   tau=tau+1;
  endo;
  xpx=xpx'xpx;
      @ Get Y'X in ypx @
   ypx=wzp[mas+1:tma,.]'ep[mas:tma-1,.];
   tau=2;
   do until tau>mas;
    ypx=ypx~(wzp[mas+1:tma,.]'ep[mas+1-tau:tma-tau,.]);
    tau=tau+1;
   endo;
      @ Get MA coefficients @
   bma=ypx*invpd(xpx);    @nmr by mas*nmr matrix@
      @ Get residual for MA in vp @
  vp=wzp[mas+1:tma,.]-ep[mas:tma-1,.]*(bma[.,1:nmr]');
  tau=2;
  do until tau>mas;
    vp=vp-ep[mas+1-tau:tma-tau,.]*(bma[.,(tau-1)*nmr+1:tau*nmr]');
   tau=tau+1;
  endo;
  omega22=eye(nmr);
  tau=1;
  do until tau>mas;
   omega22=omega22+bma[.,(tau-1)*nmr+1:tau*nmr];
   tau=tau+1;
  endo;
  omega22=(omega22*(vp'vp)*(omega22'))/(tma-mas);
  w=invpd(omega22);
 goto lnexite;

@ The QS kernel @
 lnq:
if maxd==0;
 ? "-----For the next GMM iteration, the QS kernel (nonprewhitened) is
used for estimation of omega-----";
 wvstr=wzp;
 omegas=(wvstr'wvstr)/tendm;
 tenda=tend;
else;
 ? "-----For the next GMM iteration, the QS kernel (prewhitened) is
used for estimation of omega-----";
 tenda=tend-1;
 als=wzp[2:tend,.]'wzp[1:tend-1,.]*invpd(wzp[1:tend-1,.]'wzp[1:tend-1,.]);
 {va,bls}=eigrs2(als*(als'));
 {va,clsc}=eigrs2(als'als);
 dells=bls'als*clsc;
"ALS and DELTALS for prewhitening=" als dells;
"maxd=" maxd;
 indm=0;
 i=1;
 do until i>rows(als);
  if dells[i,i]>maxd;
   ? i "-th element of DeltaLS=" dells[i,i] "is replaced by" maxd;
   indm=1;
   dells[i,i]=maxd;
  elseif dells[i,i]<-maxd;
   ? i "-th element of DeltaLS=" dells[i,i] "is replaced by" -maxd;
   indm=1;
   dells[i,i]=-maxd;
  endif;
  i=i+1;
 endo;
 if indm==1;
  als=bls*dells*(clsc');
 endif;
 wvstr=wzp[2:tend,.]-wzp[1:tend-1,.]*(als');
 omegas=(wvstr'wvstr)/tendm;
endif;
if st==0;
? "***** Automatic Bandwidth Estimator is used";
 clear alpha2,alpha2d;
 i=1;
 do until i>cols(wzp);
  if wav[i,1]>0;
   iwsws=invpd(wvstr[1:tenda-1,i]'wvstr[1:tenda-1,i]);
   rhoa=iwsws*(wvstr[1:tenda-1,i]'wvstr[2:tenda,i]);
   ea=wvstr[2:tenda,i]-wvstr[1:tenda-1,i]*rhoa;
   siga4=((ea'ea)/(tenda-1))^2;
   alpha2=alpha2+4*wav[i,1]*(rhoa^2)*(siga4)/((1-rhoa)^8);
   alpha2d=alpha2d+wav[i,1]*siga4/((1-rhoa)^4);
  endif;
  i=i+1;
 endo;
 alpha2=alpha2/alpha2d;
 st=1.3221*((alpha2*(tend))^0.2);
 if st>bst;
  "(estimated st)=" st "is > bst, bst is used";
   st=bst;
 endif;
endif;
 ? "The bandwidth parameter used" st;
 v2512=25/(12*(pi^2));
 v65=6*pi/5;
 clear lamdas;
 tau=1;
 do until tau>tenda-1;           @change from >= 5/14/91@
  omegat=(wvstr[1:tenda-tau,.]'wvstr[1+tau:tenda,.])/tendm;
  xtau=tau/st;
  omegas=omegas+v2512*(sin(v65*xtau)/(v65*xtau)-cos(v65*xtau))
                     *(omegat+omegat')/(xtau^2);
  tau=tau+1;
 endo;
if maxd==0;
 omega22=omegas;
else;
 dinv=inv(eye(rows(als))-als);
 omega22=dinv*omegas*(dinv');
endif;

  w=invpd(omega22);
goto lnexite;

@ Next Iteration over W0 @
lnexite:
 differ=w-w0;differ=maxc(maxc(abs(differ)));
 ? "   max(|difference|)=" differ;
 w0=w;
 itergm=itergm+1;
endo;
labe:
retp(chi);
endp;
@ ------------------------ End of Program ------------------------------ @