proc cmlecx dylist rho srho vard svard logl conv ;
?
? ========================= CMLECX ===============================
?
?  Univariate Time Series Process 
?  Estimation by CMLE (Kruiniger 1999)
?  The model is
?
?       y(i,t) = a(i) + rho*y(i,t-1) + e(i,t)
?
?  Assume that overall year means d(t) have already been removed
?  from the data.
?
?  Hall-Mairesse paper March 2001.
?  This version uses the concentrated likelihood function
?  so the variance is derived, not estimated.
?
?  The user must have defined the following variables in the
?  main program:
?    t      Number of time periods before differencing.
?    nob    Number of individuals (units).
?    dylist A list with the difference series (T-1 of them,
?           each with N observations). 
?    
  set t1 = t-1 ; set nobt1 = nob*t1 ;
  set nobt12 = nobt1/2 ;
  set log2pi = log(2*3.14159) ;
  set const = -nobt12*(log2pi+1) ;  ? Log L constant plus trace (kernel).

?  The code below only needs to be executed once, so
?  we do not include it in the procedure passed to ML PROC.

  mform (nrow=t1,ncol=t) d ;
  mform (nrow=t,type=sym) vrho ;
  
  do i = 1 to t1 ;
    set vrho(i,i) = 1 ;
    set d(i,i) = 1 ;
    set i1 = i+1 ;
    set d(i,i1) = -1 ;
  enddo ;
  set vrho(t,t) = 1 ;
  
  smpl 1 nob ;
  mmake dy dylist ;
  mat zmat = dy'dy ;
  
? If srho is negative, then fix rho at the starting value.
  if srho<0 ; then ; do ;
    set vard = 1 ;
    ml (silent,maxit=100) cmlec vard ;
  enddo ;
  if srho>=0 ; then ; do ;
    ml (silent,maxit=100) cmlec rho vard ;
	unmake @ses srho svar ;
  enddo ;
  set svard = vard/sqrt(nobt12) ; ? Estimated std error of the variance.
  set logl = @logl ;			  ? Return the log likelihood.
  set conv = @ifconv ;			  ? Return the convergence status.
  
endproc cmlecx ;
?===================================================================
proc cmlec ;
?
?     This proc evaluates the likelihood function for the
?     AR1 panel data model with fixed effects (and constant variance).
?     The only parameter is rho, although the variance of the
?     disturbance is also evaluated.
?     The only difference between the evaluation of logl for
?     rho=(-1,1) and rho=1 is that in the latter case, the covariance
?     matrix has a different form (V(rho) is identity).
?
  if abs(rho)<1 ; then ; do ;    	? Branch for stationary rho.

?   Fill in the off-diagonals of the matrix V(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 ;
?   Evaluate the matrix phi and the variance of the disturbances.
    mat phi = (d*vrho*d')/(1-rho*rho) ;
    mat phiinv = yinv(phi) ;
    mat detphi = logdet(phi) ;

  enddo ;

  if abs(rho-1)<.000001 ; then ; do ; ? Branch for rho==1.

	mat phi = ident(t1) ;
?	mat detphi = logdet(phi) ;  ? Actually this is always just zero. 
	set detphi = 0 ;
?	mat phiinv = yinv(phi) ;	? And this is identity again.
	mat phiinv = phi ;
  enddo ;

  if rho>-1 & rho<=1 ; then ; do ;
    mat vard = tr(phiinv*zmat)/nobt1 ;
    set logl = const-nob*detphi/2-nobt12*log(vard) ;
	set @logl = logl ;
  enddo ;
   
  else ; do ;  ? Branch for rho out of range.
    set @logl = @miss ;
  enddo ;
endproc cmlec ;