options double crt nwidth=17,signif=9;
name ghac;
?  HAC SEs for Panel Data in many patterns
?  same coefficients in all equations
?  Benchmarks with Grunfeld data/model
?  written by Clint Cummins  7/7/00

const n,10 T,20 ystart,1935;
freq(panel,n=n,t=T,id=firm,time=year,start=ystart) a;
set nT = n*T;
smpl 1,nT;
read(file='grunfeld.dat') firm year i f k;

olsq(terse) i c f k;
?select .not.miss(@res);  ? in case there are missing values
?mmake x @rnms;
?mat xpxi = (x'x)";
mat xpxi = @vcov/@s2;

title 'diagonal het-robust SEs - types 0-3';
olsq(hctype=0,terse) i c f k;
olsq(hctype=1,terse) i c f k;
olsq(hctype=2,terse) i c f k;
olsq(hctype=3,terse) i c f k;

title 'block diagonal het-robust SEs with df adj like hctype=1';
set dfadj = @nob/(@nob-@ncoef);

title 'het-auto robust SEs';
?  Omega = diag( Sigma_t )    -- each T x T    (or T_i x T_i if unbalanced)
?  Reference:  the middle term xpOmegax is the same as @VCOV for
?    GMM(HET) with one equation per time period (all parameters equal
?    across equations), plus it works for unbalanced.
mat xpOmegax = 0*xpxi;   ?  x'Omega*x = x1'Omega1*x1 + x2'Omega2*x2 + ...
supres smpl;
do ii=1,n;
  ?select (firm=ii) .and. .not.miss(@res);  ? in case there are missing values
  select firm=ii;
  mmake xi @rnms;
  mat xipui = xi'@res;
  mat xpOmegax = xpOmegax + xipui*xipui';   ? xipui*xipui' = xi'(ui*ui')*xi
enddo;
select 1;
mat vcovha = (xpxi'xpOmegax*xpxi)*dfadj;
tstats(names=@rnms) @coef vcovha;

title 'het-contemporaneous correlation robust SEs';
?  Omega = diag( Sigma_i )    -- each n x n    (or n_t x n_t if unbalanced)
?  (here we implicitly reorder the data to have all obs for year1,
?   then all obs for year2, etc.)
?  Reference:  the middle term xpOmegax is the same as @VCOV for
?    GMM(HET) with one equation per individual, plus it works for unbalanced.
mat xpOmegax = 0*xpxi;
supres smpl;
set yend = ystart+T-1;
do jj=ystart,yend;
  ?select (year=jj) .and. .not.miss(@res);  ? in case there are missing values
  select year=jj;
  mmake xt @rnms;
  mat xtput = xt'@res;
  mat xpOmegax = xpOmegax + xtput*xtput';   ? xtput*xtput' = xt'(ut*ut')*xt
enddo;
select 1;
mat vcovhcc = (xpxi'xpOmegax*xpxi)*dfadj;
tstats(names=@rnms) @coef vcovhcc;

? Note:  it might seem possible to handle a completely general
? Omega with this method, such as:
?   select 1;
?   mmake x @rnms;
?   mat xpu = x'@res;
?   mat xpOmegax = xpu*xpu';
? However, xpu and xpOmegax will be entirely zeros, since xpu are
? the orthogonality conditions for OLS.
?   Similarly, if the model includes coefficients which vary by individual
? or by time period (such as time dummies or individual fixed effects),
? parts of xpu and xpOmegax will be zeroes due to orthogonality
? conditions.  This sort of limitation with het-robust SEs has been
? observed before with variables that are single-observation dummies.

title 'Beck-Katz (non-het) contemporaneous correlation robust SEs';
?  Omega = Sigma_i # Ident(t)
?  Sigma_i formed as average over all individuals, so assumes no
?   conditional heteroskedasticity is present.
?  assumes balanced data  (T obs. for each individual)
? (vs. the above methods, which work for unbalanced data)
?  Reference:  Beck, Nathaniel L., and Katz, Jonathan N.,
?    "What to Do (and Not to Do) with Time-Series-Cross-Section Data",
?    American Political Science Review 89, 1995, p.634-647.
?  Note:  the middle term xpOmegax is the same as @VCOV for
?    GMM(NOHET) with one equation per individual.
?  Unbalanced:  Could be extended to unbalanced by using pairwise
?    version of (E_i'E_j)/T_ij, although this might result in a
?    non-positive-definite Sigma.  Another option would be to
?    form E_s'E_s/T_s from a smaller set of time periods T_s < T,
?    where all individuals are non-missing.
?    The term xt'sigma*xt would need to use the appropriate submatrix
?    of sigma when the number of individuals non-missing for a given
?    year is less than n.  The above methods do not involve any
?    of these complications.
mform(nrow=t,ncol=n) E = @res;
mat sigma = (E'E)/T;  ? n x n
mat xpOmegax = 0*xpxi;
supres smpl;
set yend = ystart+T-1;
do jj=ystart,yend;
  ?select (year=jj) .and. .not.miss(@res);  ? in case there are missing values
  select year=jj;
  mmake xt @rnms;
  mat xpOmegax = xpOmegax + xt'sigma*xt;  ? If unbalanced, xt would not be n x k
enddo;
select 1;
mat vcovbk = (xpxi'xpOmegax*xpxi)*dfadj;
tstats(names=@rnms) @coef vcovbk;