options double crt nwidth=17,signif=9;
name gar1rei;
?  AR1, Random Effects Individual (exact ML) on all 10 firms.
?  Follows Baltagi and Li, Journal of Econometrics 48, 1991, p.385-393.
?  Parameterized as rho_ic, sigma2.
?  This allows for a negative rho_ic (intra-class correlation coefficient)
?  if necessary.

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;
ar1(tscs) i c f k;

d1 = (year=ystart);
d2T = 1-d1;
dT = (year=(ystart+T-1));
frml equ u = i - (b_0 + b_f*f + b_k*k);
param b_0 b_f b_k rho;
frml eqe d1*u*sqrt(1-rho**2) + d2T*(u - rho*u(-1));
frml ar1 logl = -log(sigma2)/2 + lnorm(eqe/sqrt(sigma2)) + d1*log(1-rho**2)/2;
param sigma2,1;
eqsub ar1 eqe equ;
ml(hiter=n,hcov=nb,nodropmiss,tol=1e-7,maxit=30) ar1;

set pi = 4*atan(1);
set lr2pi = log(sqrt(2*pi));
frml ar1rei logl = -log(sig_e) + (d1/2)*(log(1-rho**2) + log(theta)) - lr2pi
                    - upSiu/2 ;
frml etheta theta = sig2_e/sig2_a;
frml upSiu ystars**2/sig2_e;
frml ystars eqe - thetaA*b*(d1*alpha + d2T);
frml ethetaA thetaA = 1 - sig_e/sig_a;
frml esig_e sig_e = sqrt(sig2_e);
frml esig_a sig_a = sqrt(sig2_a);
frml esig2_a sig2_a = d2*sig2_mu*(1-rho)**2 + sig2_e;
frml ed2 d2 = alpha**2 + T-1;
frml ealpha alpha = sqrt((1+rho)/(1-rho));
frml esig2_e sig2_e = sigma2*(1-rho_ic);
frml esig2_mu sig2_mu = sigma2*rho_ic;
frml eb b = bn/d2;
frml ebn bn = (1+rho)*u_1 + us2 - rho*us1;
dot _1 s1 s2;
  frml u. i. - (b_0*c. + b_f*f. + b_k*k.);
enddot;
eqsub ar1rei
   etheta upSiu ystars ethetaA esig_a esig2_a eb ebn ed2 ealpha
   esig_e esig2_e esig2_mu
   eqe equ u_1 us1 us2;
? Create i_1 is1 is2 f_1 fs1 fs2 k_1 ks1 ks2 (sums of data)
supres smpl;
dot i c f k;
    ._1 = .;  ? first obs.
    .s1 = .;  ? Ultimately, sum from 1 to T-1
    .s2 = 0;  ? Ultimately, sum from 2 to T
    select d2t;  ? from 2 to T
      genr(silent) .s1 = .s1(-1) + . ;  ? sum from 1 to T
    select dT;  ? obs T
      mmake mtmpsum .s1;  ? extract obs T - the sum from 1 to T
      mmake mtmp_T  .;    ? extract obs T
    select d1;  ? obs 1 - put sums into first obs for replication
      unmake mtmpsum tmpsum;
      unmake mtmp_T tmp_T;
      .s1 = tmpsum - tmp_T;  ? sum from 1 to T-1 = sum(1,T) - y_T
      .s2 = tmpsum - ._1;    ? sum from 2 to T   = sum(1,T) - y_1
    select d2t;  ? from 2 to T - replicate first obs to all obs
      genr(silent) ._1 = ._1(-1);
      genr(silent) .s1 = .s1(-1);
      genr(silent) .s2 = .s2(-1);
    select 1;
enddot;
const rho_ic,0;
ml(hiter=n,hcov=nb,nodropmiss,tol=1e-7) ar1rei;

title 'with RHO=0 fixed, rho_ic estimated';
const rho,0;
param rho_ic,.5;
ml(hiter=n,hcov=nb,nodropmiss,tol=1e-7) ar1rei;
eqsub esig_e esig2_e;
frml esig_mu sig_mu = sqrt(sig2_mu);
eqsub esig_mu esig2_mu;
analyz esig_e esig_mu;

title 'full model';
param rho;
ml(hiter=n,hcov=nb,nodropmiss,tol=1e-7,maxit=50) ar1rei;

title 'hardcoded version';
ar1(rei,tol=1e-7) i c f k;