proc cmlechsx dylist rho srho var svar logl conv ;
?
? ========================= CMLECHSX ===============================
?
?  AR(1) Model for a Panel Univariate Time Series Process
?  with individual-specific mean and variance.
?
?  Estimate by Concentrated CMLE (Kruiniger 1998, 1999 as enhanced
?  by Hall and Mairesse 2001). The model is
?
?       y(i,t) = a(i) + rho*y(i,t-1) + e(i,t)
?
?  (This assumes that year means d(t) have already been removed).
?  The errors are heteroskedastic with Var(e) = sigsq(i) sigsq(t).
?  This version of CMLECHSX has the option to estimate the variance(t)
?  only, or variance and rho. The likelihood is "concentrated" by
?  plugging in the equation for individual-specific variances. Note
?  that the score is likely to be biased in finite samples so
?  iteration behavior can probably be improved.
?
?  Hall-Mairesse paper revision September 2001.

  set t1 = t-1 ; set t1half = t1/2 ; set nobt1 = nob*t1 ;
  set log2pi = log(2*3.14159) ;
  set const = -(nobt1/2)*(log2pi+1) ;

  mform (nrow=t,type=sym) vrho ;
  do i = 1 to t ;
    set vrho(i,i) = 1 ;
  enddo ;
	
  mform (nrow=t1,ncol=t) d ;
  do i = 1 to t1 ;
    set d(i,i) = 1 ;
    set i1 = i+1 ;
    set d(i,i1) = -1 ;
  enddo ;

  dot 1-12 ;
   param sig. 1 ;
  enddot ;

  smpl 1 nob ;
  sigsq = 0 ;
  
? If srho is negative, then fix rho at the starting value.
  if srho<0 ; then ; do ;
    ml (silent,maxit=100) cmlechs sig3-sig12 ;
?	print sig2-sig12 ;
  enddo ;
  if srho>=0 ; then ; do ;
    ml (silent,maxit=100) cmlechs rho sig3-sig12 ;
?	print sig2-sig12 ;
	unmake @ses srho ssig3-ssig12 ;
  enddo ;
  set logl = @logl ;
  set conv = @ifconv ;
  smpl 1 nob ;
  logsig = log(sigsq) ;
  msd (silent) sigsq logsig ;		?  Get the sigsq (i) mean.
  unmake @mean var logvar ;
  set svar = var/sqrt(nobt1/2) ; ? Estimated std error of the variance.

?  hist sigsq ;
?  hist logsig ;
endproc cmlechsx ;
?===================================================================
proc cmlechs ;
?
?     Proc for log likelihood of AR1 panel data model
?     with fixed effects (and heteroskedastic errors).
?     Concentrated likelihood (s(i) removed).
?     Dy(i) ~ N[0,D S(i) V	S(i) D']
?           = N[0,sigsq(i) DPVPD']

?     where D is the differencing operator
?     sigsq(i) is the variance of the ith obs.
?     V(rho) is the AR1 covariance matrix
?     P is the matrix of the time varying sigmas.
?     The form of the variance of e(i,t) is sigsq(i)sigsq(t).
?     sigsq(t=1) is normalized at unity.
?
  local sigt dphid dsigt sigisq phi detphi phiinv ;

  set logl = const ;  ? This includes the trace term.
  
  if abs(rho)<1 ; then ; do ;	    ? Branch for stationary rho.
    do i = 2 to t ;
      set jend = t-i+1 ;
      do j = 1 to jend ;
	    set ii = j+i-1 ;
        set vrho(ii,j) = rho**(i-1) ;
	  enddo ;
    enddo ;
    mmake sigt sig1-sig12 ;
    mform (nrow=t,type=diag) sigt ;
	mat dsigt = d*sigt ;
	mat phi = (dsigt*vrho*dsigt')/(1-rho*rho) ;
  enddo ;

  if abs(rho-1)<.000001 ; then ; do ;		?  Branch for rho==1.
    mmake sigt sig2-sig12 ;
    mform (nrow=t1,type=diag) sigt ;
	mat phi = sigt*sigt ;
  enddo ;

  if rho>-1 & rho<=1 ; then ; do ;  
    mat phiinv = yinv(phi) ;
    mat detphi = logdet(phi) ;  
    set logl = logl-nob*detphi/2 ;
    do i = 1 to nob ;
      smpl i i ;
      mmake dy dylist ;
      mat sigisq = (dy*phiinv*dy')/t1 ;
      sigsq = sigisq ;		  ? Store the estimate of the ith var.
      set logl = logl - (t1half)*log(sigisq) ; 
    enddo ;
    set @logl = logl ;
  enddo ;
  
  else ; do ;	                          ? Branch for illegal rho.
    set @logl = @miss ;
  enddo ;

endproc cmlechs ;
?===================================================================