@ GMMP2S.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:  05/29/91
------------------------------------------------------------------------------@
/*
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 GMMP.  This program can be used in a
   similar way as GMMQ.SET See GMMQ.EXP@
INPUT:
nendv    a global variable a column vector of nend(j) for j=1,...,q
nz       the number of orthogonal conditions for each j
kgmv     the mubmer of the parameters estimated in each village.
         this is used for small sample correction of st. error only.
             (this feature is added at 21:45 on 11/13/91)

Define a procedure

hu(b)   sumc(nendv) by nz vector of disturbance with first nendv[1] for j=1,
       the next nendv[2] for j=2, etc.

*/
@ Globals @
clear bgm,const,const2,W0,gufn,w0flag,calwflag,nv,nendv;

proc hu(b); @ This is a dummy definition, and not used by GMM programs.@
retp(b);
endp;

proc sohu2(b12); @ This is a dummy definition, and not used by GMM programs.@
retp(b12);
endp;

proc (1)=ufn(b);
local u,suv,sib,j;
u=const*hu(b);
suv=sumc(u[1:nendv[1],.]);
sib=nendv[1]+1;
j=2;
do until j>rows(nendv);
 suv=suv|sumc(u[sib:sib+nendv[j]-1,.]);
 sib=sib+nendv[j];
 j=j+1;
endo;
retp(suv);
endp;

proc startval;
retp(bgm);
endp;

proc gmmp2s(&gradg,bgm1,del1,hu1,nendn,nz,kgmv,zero,maxitegm);
local gradg:proc;  @ used for gamma21 and gamma22, the global, gradname
     controls gradient calculation in minquad.set@
local df,chi,differ,itergm,vbgm,sd,prob,wzp,omep,tau,ometau,
      snzi,snzi0,igm,w,vp,prob0,differp,o11,o11j,o12,o12j,o22,o22j,
      bmai,i,ib,ie,rzw,rv,rvwz,ep,tma,xpx,ypx,bma,arcoef,omega2,
      nod,sub,sib,j,gamstr,wj,sie,kgm,kgm1,gam2,gam21,gam22,gdi,
      omestr,nzp,nz1,nendnk;
clear omep;
nod=rows(nendv);
kgm=rows(bgm);
kgm1=rows(bgm1);
nzp=round(nz*nod);
nz1=cols(hu1);
nendnk=nendn-kgmv;
? "1st step bgm1=" bgm1';
? "Initial values of the coefficients=" bgm';
? "kgmv=" kgmv;
format /m1 16,8;
df=nzp-kgm;
prob0=100;
clear chi;
differ=1e+12;
differp=differ;
if w0flag==1;
 clear w0;
 ? "initial values of parameters are used to caluculate initial W0";
 gam2=gradg(&sohu2,bgm1|bgm,sohu2(bgm1|bgm),0);
 gam21=gam2[.,1:kgm1]/nendn;
 itergm=0;
 goto lab100;
elseif w0flag==0;
 w0=eye(nzp);
 ? "Initial W0=I";
 ? "Scaling Constant used for first iteration=" const;
 itergm=1;
 bgm=minquad;
 gam2=gradg(&sohu2,bgm1|bgm,sohu2(bgm1|bgm),0);
 gam21=gam2[.,1:kgm1]/nendn;
 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;
do until differ<zero;
     const=const2;
     bgm=minquad;
lab3200:
output on;
     save bgm;
     save w0;
 gam2=gradg(&sohu2,bgm1|bgm,sohu2(bgm1|bgm),0);
 gam21=gam2[.,1:kgm1]/nendn;
 gam22=gam2[.,kgm1+1:kgm1+kgm];
 vbgm=nendn*invpd(gam22'w0*gam22);
 sd=sqrt(diag(vbgm));
 save varb=vbgm;

   @ save variance for bgm [not sqrt(T)*bgm] in VARB.FMT @
? "
================= GMM ITERATION " itergm "====================";
 ? "   b=" bgm';
 ? "   s.d.=" sd';
   if df>0;
    chi=(vof)./((const^2).*nendn);
    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;? "   nzp<kgm; Not Identified";endif;
 ? "   scaling const used is" const;
 if itergm>maxitegm-0.5;  @ or differp<zero;@
  goto labe;
 endif;
 lab100:
 wzp=hu(bgm);
 o11j=(hu1[1:nendv[1],.]'hu1[1:nendv[1],.])/nendnk;
 o12j=(hu1[1:nendv[1],.]'wzp[1:nendv[1],.])/nendnk;
 o22j=(wzp[1:nendv[1],.]'wzp[1:nendv[1],.])/nendnk;
 sie=nendv[1];
 o11=o11j;
 o12=o12j;
 o22=o22j;
 j=2;
 do until j>rows(nendv);
   o11j=((hu1[sie+1:sie+nendv[j],.]'hu1[sie+1:sie+nendv[j],.])/nendnk);
  o11=(o11~zeros(nz1*(j-1),nz1))|(zeros(nz1,nz1*(j-1))~o11j);
   o12j=((hu1[sie+1:sie+nendv[j],.]'wzp[sie+1:sie+nendv[j],.])/nendnk);
  o12=(o12~zeros(nz1*(j-1),nz))|(zeros(nz1,nz*(j-1))~o12j);
   o22j=((wzp[sie+1:sie+nendv[j],.]'wzp[sie+1:sie+nendv[j],.])/nendnk);
  o22=(o22~zeros(nz*(j-1),nz))|(zeros(nz,nz*(j-1))~o22j);
  sie=sie+nendv[j];
  j=j+1;
 endo;
 omep=(o11~o12)|((o12')~o22);
 gdi=(-gam21*del1)~eye(nzp);
 omestr=gdi*omep*(gdi');
 trap 1;
 w=invpd(omestr);
 trap 0;
 if not scalerr(w);
  goto lnexite;
 endif;
  ? "!!!!! ERROR !!!!! OMEGA should be p.s.d. but is not p.d.";
  output off;
  stop;

@ 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 ------------------------------ @