options double crt nwidth=17,signif=9;
name grsur 'Grunfeld data - SUR (ML) by firms';
?  Iterated SUR (ML) on all 10 firms
?  same coefficients in all equations
?  I.e. unrestricted firm vcov, but assuming
?  no autocorrelation across time.
?  This allows a non-diagonal heteroskedasticity.

?  Using the SBIC to compare different models for the Grunfeld
?  data/model, SUR without AR(1) fits better than diagonal het
?  without AR(1).  However, diagonal het with AR(1) (SBIC=874)
?  fits better than SUR with AR(1) and better than either of the above.
?
?  Before any conclusions are drawn, these het models should
?  also be fit with individual fixed effects to see if SBIC can
?  be lowered even further.  Models with an explicit lagged
?  dependent variable should also be considered, although this
?  makes the SBIC comparison more difficult.  That is, a common
?  factor test should be run to see if the AR(1) model is just
?  proxying for missing lags of Y and X.  (AR(1) could also be
?  proxying for missing fixed effects).

?  Usually (when N > T) the SUR/FIML model is applied with one
?  equation per time period, instead of one equation per firm.
?  This allows for unrestricted residual autocorrelation.
?  However, when N < T as in the Grunfeld data, this is infeasible.
?  This model normally assumes no individual heteroskedasticity.
?  However, a SIGI equation could be used with FIML like in
?  grhetar.tsp  to allow for this.

?  by Clint Cummins 11/1999

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 the data from the text data file, with variables numbered
? by firm, so that one equation can be written for each firm.
freq a;
smpl 35,54;
dot 1-10;
  read(file='grunfeld.dat') firm year i. f. k.;
  frml eq. i. = b0 + b1*f. + b2*k.;
enddot;
param b0 b1 b2;
supres regout covt w;

title 'Iterated SUR (ML) estimates';
? Proposed Benchmark:
? 1. SUR/FIML - by firm, no autocorrelation
?                 TSP
? F           .0387618741
? (SEHess)   (.001802790101)
? (SEBHHH)   (.001828408693)
? K           .0735533461
?            (.002102205745)
?            (.002393898222)
? Constant  -1.566557384
?            (.1808473092)
?            (.1876072270)
? LogL    -798.8881471
? SBIC     949.890192094
? dfk       57
?  Note:  see  grunsur2.tsp  for a published benchmark of SUR on
?  just 2 of the Grunfeld firms.

lsq(tol=1e-11,step=bardb,maxit=30) eq1-eq10;  ? same as FIML estimates
set dfk = 3 + (10*11/2) - 1; ? 57
csbic @logl dfk nt;
Proc csbic logl dfk nobneq;
  set sbic = -logl + dfk*log(nobneq)/2;
  print sbic;
endproc;

? SUR + AR(1) model via exact ML / FIML
? alternative ways of doing exact ML (which can be extended)
d1 = (year=ystart);
dot 1-10;
  frml u.  i. - (b0 + b1*f. + b2*k.);
  frml e.  d1*u.*sqrt(1-rho**2) + (1-d1)*(u. - rho*u.(-1));
  eqsub e. u.;
enddot;
title 'FIML reproduces iterated SUR, although with BHHH SEs';
const rho,0;
fiml(endog=(i1-i10),nodropmiss,tol=1e-7) e1-e10;

title 'FIML/SUR with AR(1) - exact ML';
? Proposed Benchmark:
? 2. SUR/FIML - by firm, with AR(1) exact ML
?                 TSP
? F           .0362672946
? (SEBHHH)   (.004083578417)
? K           .0609158857
?            (.007680004072)
? Constant   -.6899939688
?            (.5229315423)
? RHO         .6912210325
?            (.0677486515)
? LogL    -778.0180211
? SBIC     931.669224736
? dfk       58

param rho;
fiml(endog=(i1-i10),nodropmiss,tol=1e-7) e1-e10;
set dfk = 4 + (10*11/2) - 1; ? 58
csbic @logl dfk nt;