Article From: Department of Economics
Subject: GMMTRE.SET

Posted by: Information Provider
Phone Number: (202) 885-3770
E-Mail Address: econ@american.edu

Post Date: 28 Sep 1994
Expiration Date: 28 Oct 2004


               GMMTRE.SET (Hansen, Ogaki, & Heaton)

@ GMMTRE.SET @

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

 Last Revision:  12/15/87
------------------------------------------------------------------------------@
/*
   @@ This program sets procedure GMMTRE.  See GMMTRE.DOC for documentation.

      This program checks positive definiteness of W0.  If not, then this
      program uses the modified Durbin's method or Newey and West method.


    This program is for simultaneous estimation for time trends and nonlinear
    Euler equations.
*/
#include "partinv.prc";
#include "ywp.prc";
@ Globals @
clear bgm,const,const2,W0,gufn,w0flag,calwflag,nv,ordmav,ordard,seq11t,tend;
lend=15;
proc hu(b); @ This is a dummy definition, and not used by GMM programs.@
retp(b);
endp;

proc (1)=ufn(b);
local u,htp,vp;
u=hu(b);
htp=sumc(u);
vp=u[.,1:nv];
vp=seq11t.*vp;
retp(const*((sumc(vp)./tend)|htp));
endp;

proc startval;
retp(bgm);
endp;

proc (0)=gmmtre(&gradq,tend,nzv,rwv,nw,kgm,zero,maxitegm);
local gradq:proc;
local L,nz,maxrw,df,chi,differ,itergm,vbgm,sd,prob,wzp,omega22,tau,ometau,
      snzi,snzi0,igm,w,zero,maxitegm,yesno,ugm,omega11,omega12,vp,one3,ib,ie,
      i,nzpnv,rzw,rv,rvwz,ep,tma,xpx,ypx,bma,arcoef,omega2,differp,prob0;
? "lend=" lend;
? "rwv=" rwv';
one3=(1/3~0.5)|(0.5~1);
prob0=10;
seq11t=seqa(1,1,tend);
nz=sumc(nzv);
maxrw=maxc(rwv);
L=nz+2*nv;
nzpnv=nz+nv;
? "Initial values of the coefficients=" bgm';
format /m1 16,8;
df=L-kgm;
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(L);
 ? "Initial W0=I";
 ? "Scaling Constant used for first iteration=" const;
 itergm=1;
 bgm=minquad;
 save bgmi=bgm;
 goto lab100;
endif;
const=const2;
itergm=1;
if w0flag==3;
 vof=ofn(bgm);
 ? "W0 in the memory and initial bgm are used for the first GMM result";
 goto lab3200;
endif;
? "W0 in the memory is used as initial W0";
do until differ<zero;
     const=const2;
     bgm=minquad;
     save bgm;
lab3200:
output on;
 gufn=gradq(&ufn,bgm,ufnv,0);
 vbgm=tend*invpd(gufn'w0*gufn./(const^2));
 sd=sqrt(diag(vbgm));
 save  varb=vbgm;
    @ Save variance for bgm [not sqrt(T)*bgm] in VARB.FMT @
          @ calculate S.D. for (1/sqrt(T))*g(b) @
 vbgm=invpd(w0)-(gufn*vbgm*(gufn'))/((const^2)*tend);
 vbgm=sqrt(diag(vbgm));
? "
================= GMM ITERATION " itergm "====================";
 ? "   b=" bgm';
 ? "   s.d.=" sd';
   if df>0;
    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;
 elseif df==0;? "   Just Identified";
 else;? "   L<kgm; Not Identified";endif;
 ? "(1/sqrt(T))g(b)=" (ufnv')/(const*sqrt(tend));
 ? "           s.d.=" vbgm';
 ? "   scaling const used is" const;
 if itergm>maxitegm-0.5 or differp<zero;
  goto labe;
 endif;
 lab100:
 ugm=hu(bgm);
 vp=ugm[.,1:nv];
 wzp=ugm[.,(nv+1):nzpnv];
 clear ugm;
 omega11=vp'vp;
 clear omega12;
 rvwz=(omega11./tend)~zeros(nv,nz);  @ this is for Durbin's method @
 tau=1;
 do until tau>(lend-1);
  ometau=vp[(1+tau):tend,.]'vp[1:(tend-tau),.];
  rv=ometau./(tend-tau);
  omega11=omega11+(1.0-tau./lend).*(ometau+ometau');
  ometau=vp[(1+tau):tend,.]'wzp[1:(tend-tau),.];
  rvwz=rvwz~(rv~(ometau./(tend-tau)));
  omega12=omega12+(1.0-tau./lend).*ometau;
  tau=tau+1;
 endo;
 clear ometau,rv;
 omega11=omega11./tend;
 omega12=omega12./tend;
 omega22=wzp'wzp;
 omega2=omega22;  @ Used when omega is not p.d. @
 rzw=zeros(nz,nv)~(omega22./tend);@ for Durbin's method @
 tau=1;
 do until tau>=maxrw;
  ometau=wzp[(1+tau):tend,.]'wzp[1:(tend-tau),.];
  omega2=omega2+ometau+(ometau');
  clear snzi;
  igm=1;
  do until igm>nw;
   snzi0=snzi+1;
   snzi=snzi+nzv[igm,.];
   if tau>=rwv[igm,.];
    ometau[.,snzi0:snzi]=zeros(sumc(nzv),nzv[igm,.]);
   endif;
   igm=igm+1;
  endo;
  rzw=rzw~(zeros(nz,nv)~(ometau./(tend-tau)));
  omega22=omega22+ometau+ometau';
  tau=tau+1;
 endo;
 omega12=(0.5|1).*.omega12;
 omega22=omega22./tend;
 omega22=((one3.*.omega11)~omega12)|((omega12')~omega22);
 if calwflag==2;
  goto ldurb;
 elseif calwflag>2;
  goto lnw;
 endif;
 trap 1;
 w=invpd(omega22);
 trap 0;
 if not scalerr(w);
  ? "------- For next GMM iteration
  All the zero restrictions are successfully imposed on W0 -------";
  goto lnexite;
 endif;
 if minc(rwv) /= maxrw;
  ? "------- OMEGA is not p.d. Now some of the zero restirictions
     on OMEGA22 are neglected -------";
  omega2=((one3.*.omega11)~omega12)|((omega12')~(omega2./tend));
  trap 1;
  w=invpd(omega2);
  clear omega2,omega11,omega12;
  trap 0;
  if not scalerr(w);
   ? "-------- Zero restirictions for t>=maxrw successfully imposed ----";
   goto lnexite;
  endif;
  ? "------- OMEGA is still not p.d. ---------";
 endif;
 if calwflag==0;
 ? "OMEGA is not positive definite when OMEGA12 and OMEGA22 are
 calculated with zero restriction.  Now the modified Durbin's method
 is used to calculate W0.";
  goto ldurb;
 endif;
 goto lnw;

@ 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 for v(t)=" ordmav;
  tau=lend;
  do until tau>ordard;
   ometau=vp[(1+tau):tend,.]'vp[1:(tend-tau),.];
   rv=ometau./(tend-tau);
   ometau=vp[(1+tau):tend,.]'wzp[1:(tend-tau),.];
   rvwz=rvwz~(rv~(ometau./(tend-tau)));
   tau=tau+1;
  endo;
  arcoef=yw(rvwz[.,1:nzpnv*(ordard+1)]|
         (rzw~zeros(nz,nzpnv*(ordard+1-maxrw))));
  clear rvwz,rzw;
  @ 2nd Step: OLS to get MA coefficients @
  wzp=vp~wzp;
  clear vp;
    @ Get e(t)'s in ep@
  ep=wzp[(ordard+1),.]-reshape(wzp[1:ordard,.],1,nzpnv*ordard)*(arcoef');
  tau=2;
  do until tau>tend-ordard;
   ep=ep|(wzp[(ordard+tau),.]-
       reshape(wzp[tau:ordard+tau-1,.],1,nzpnv*ordard)*(arcoef'));
  tau=tau+1;
  endo;
  wzp=wzp[ordard+1:tend,.];
  tma=tend-ordard;
   @ A.Regress v(t) on lagged e(t)'s @
      @ Get X'X in xpx@
  igm=1;
  xpx=zeros(ordmav*nzpnv,ordmav*nzpnv);
  do until igm>ordmav; @row@
   tau=1;
   do until tau>ordmav; @column@
    xpx[(igm-1)*nzpnv+1:igm*nzpnv,(tau-1)*nzpnv+1:tau*nzpnv]=
       ep[ordmav+1-igm:tma-igm,.]'ep[ordmav+1-tau:tma-tau,.];
    tau=tau+1;
   endo;
   igm=igm+1;
  endo;
      @ Get Y'X in ypx @
  ypx=wzp[ordmav+1:tma,1:nv]'ep[ordmav:tma-1,.];
  tau=2;
  do until tau>ordmav;
   ypx=ypx~(wzp[ordmav+1:tma,1:nv]'ep[ordmav+1-tau:tma-tau,.]);
   tau=tau+1;
  endo;
      @ Get MA coefficients for v(t) @
  bma=ypx*invpd(xpx);
      @ Get residual for MA in vp @
  vp=wzp[ordmav+1:tma,1:nv]-ep[ordmav:tma-1,.]*(bma[.,1:nzpnv]');
  tau=2;
  do until tau>ordmav;
   vp=vp-ep[ordmav+1-tau:tma-tau,.]*(bma[.,(tau-1)*nzpnv+1:tau*nzpnv]');
   tau=tau+1;
  endo;
    @ Regress w(t)z(t) on e(t) @
  vp=vp~wzp[ordmav+1:tma,(nv+1):nzpnv];
        @ Regress wi(t)zi(t) onto e(t) for i=1,..,nw @
  i=1;
  ie=nv;
  do until i>nw;
  ib=ie+1;
  ie=ib+nzv[i,1]-1;
  if rwv[i,1]<1.5;
   bma=bma|zeros(nzv[i,1],ordmav*nzpnv);
   goto linex;
  endif;
      @ Get X'X in xpx@
  igm=1;
  xpx=zeros((rwv[i,1]-1)*nzpnv,(rwv[i,1]-1)*nzpnv);
  do until igm>(rwv[i,1]-1);  @row@
   tau=1;
   do until tau>(rwv[i,1]-1); @column@
    xpx[(igm-1)*nzpnv+1:igm*nzpnv,(tau-1)*nzpnv+1:tau*nzpnv]=
       ep[maxrw-igm:tma-igm,.]'ep[maxrw-tau:tma-tau,.];
    tau=tau+1;
   endo;
   igm=igm+1;
  endo;
      @ Get Y'X in ypx @
  ypx=wzp[rwv[i,1]:tma,ib:ie]'ep[(rwv[i,1]-1):tma-1,.];
  tau=2;
  do until tau>(rwv[i,1]-1);
   ypx=ypx~(wzp[rwv[i,1]:tma,ib:ie]'ep[rwv[i,1]-tau:tma-tau,.]);
   tau=tau+1;
  endo;
      @ Get MA coefficients for zw(t) @
  bma=bma|((ypx*invpd(xpx))~zeros(nzv[i,1],(ordmav-rwv[i,1]+1)*nzpnv));
  linex:
  i=i+1;
  endo;
       @ Get residual for MA in vp @
  vp[.,nv+1:nzpnv]=vp[.,nv+1:nzpnv]-ep[ordmav:tma-1,.]*
            (bma[nv+1:nzpnv,1:nzpnv]');
  tau=2;
  do until tau>(maxrw-1);
   vp[.,nv+1:nzpnv]=vp[.,nv+1:nzpnv]-
  ep[ordmav+1-tau:tma-tau,.]*(bma[nv+1:nzpnv,(tau-1)*nzpnv+1:tau*nzpnv]');
   tau=tau+1;
  endo;
  lomec:
  omega22=eye(nzpnv);
  tau=1;
  do until tau>ordmav;
   omega22=omega22+bma[.,(tau-1)*nzpnv+1:tau*nzpnv];
   tau=tau+1;
  endo;
  omega22=(omega22*vp'vp*(omega22'))/(tma-ordmav);
  omega11=one3.*.omega22[1:nv,1:nv];
  omega12=(0.5|1).*.omega22[1:nv,nv+1:nzpnv];
  omega22=(omega11~omega12)|((omega12')~omega22[nv+1:nzpnv,nv+1:nzpnv]);
  w=invpd(omega22);
 goto lnexite;

@ NEWEY and WEST Method @
 lnw:
 ? "-----For next GMM iteration Newey and West Method used for W0 ------";
  vp=vp~wzp;
  clear wzp;
  omega2=vp'vp;
  tau=1;
  do until tau>lend;
   ometau=vp[(tau+1):tend,.]'vp[1:(tend-tau),.];
   omega2=omega2+(1.0-tau./lend).*(ometau+(ometau'));
   tau=tau+1;
  endo;
  omega2=omega2./tend;
  omega11=one3.*.omega2[1:nv,1:nv];
  omega12=(0.5|1).*.omega2[1:nv,nv+1:nzpnv];
  omega2=(omega11~omega12)|((omega12')~omega2[nv+1:nzpnv,nv+1:nzpnv]);
  omega22[1:2*nv,.]=omega2[1:2*nv,.];
  omega22[.,1:2*nv]=omega2[.,1:2*nv];
  if calwflag==4;
   goto lnw4;
  endif;
  trap 1;
  w=invpd(omega22);
  trap 0;
  if scalerr(w);
   ? "OMEGA is not positive definite when OMEGA22 is calculated with zero
restrictions and OMEGA12 is calculated without zero restrictions.";
  lnw4:
  ? "No zero restrictions imposed for OMEGA12 and OMEGA22.";
   omega22=omega2;
   w=invpd(omega2);
  endif;

@ Next Iteration over W0 @
lnexite:
 differ=w-w0;differ=maxc(maxc(abs(differ)));
 ? "   max(|difference|)=" differ;
 w0=w;
 save w0;
 save omega=omega22;   @ weight before invert @
 itergm=itergm+1;
endo;
labe:
endp;
@ ------------------------ End of Program ------------------------------ @

************************************************************************
Department of Economics is entirely responsible for the information
provided above. Please direct any comments to Information Provider at:
          E-Mail Address: econ@american.edu
          Phone Number:   (202) 885-3770