? ============================ cov2step =========================
proc cov2step neq ninst mask ivs an2 ;

?  Procedure to compute the estimated two-step covariance matrix for 
?  GMM estimation (after consistent estimation of the residuals); 
?  see Arellano and Bond (REStud 91 for details). 
?  Benchmarked with Steve Bond using DPD and these give estimates that
?  match his.  TSP's default GMM estimator does not match the cov2step
?  results - the default GMM is more like a 3-step estimator since it
?  does 2SLS, 3SLS, and then GMM. 

?  This version is based on a program by Benoit Mulkay. 
?  by Bronwyn H. Hall, Jan. 97;  minor revisions by Clint Cummins 3/97
?  Revised by C.C. 8/97 to handle discontinuous SMPL.
?  requires TSP 4.3 or later   (nooptcov option in GMM, when an2 is used)

?  neq   - number of equations (time periods T-2).
?  ninst - number of instruments per equation.
?  mask  - a ninst by neq matrix of zeros and ones that specifies 
?          which instruments are used for which equation. 
?  ivs   - a list of instrumental variables (ninst long). 
?  an1   - the covariance matrix estimate (returned). 

   local noc i1 mz ineq idt i z zi zm zj reslist lsmp ifcont smpful obs ;
   copy @smpl savsmp ;              ? Save the current sample.
   mat lsmp = nrow(savsmp) ;
   set ifcont = lsmp .eq. 2;
   copy @nob nob ;

?   Get residuals. 

   list(last=neq,prefix=r) reslist ;
   local reslist ;
   unmake @res reslist ;
     
?  Now compute the estimate of An by looping over the panel and 
?  summing over the individual estimate (this version is 
?  heteroskedastic-consistent). 
  
  mat mz = ident(neq) # vech(ident(ninst)) ;
  mat idt = ident(neq) ;
  mat ineq = vech(idt) ;  ? neq x 1  unit vector
    
?  Make sure mask is a row vector.
  set noc = neq*ninst ; 
  mform (nrow=noc,ncol=1) mask ;

  mform(nrow=noc,type=sym) an2 = 0 ;
  if .not.ifcont; then; do;
    set i = savsmp(1);
    set i1 = savsmp(lsmp);
    mmake smpful i i1;
    smpl smpful;
    obs = 0;
    smpl savsmp;
    trend obs;
    smpl smpful;
  enddo;
  do i = 1 nob ;
    if ifcont; then;
      smpl i i ;
    else;
      select obs = i ;
    mmake z ivs ;
    mmake nu reslist ;
    mat zi = ineq # z' ;
    mat zm = zi % mask ;
    mform (nrow=noc,ncol=1) zm ;
    mat zj = mz % (ineq'#zm) ;
    mat an2 = an2+(zj*nu'nu*zj') ;
?  Second method of computation - more space but simpler. 
?  May be easier to follow. It gives identical results.
?  mat zs = (idt # z)*diag(mask);
?  mat an2 = an2 + (zs'nu'nu*zs) ;
  enddo ;
  mat an2 = an2/nob ;
  smpl savsmp ;
endproc cov2step ;
? ==================================================================