/***************************/
@ Do MLE of serial correlated panel , ignoring timing @

   proc modelre0(f1,parname1);

  local cv,logl,status,tmpe,tmpl,bbs,g1s,li,tmp,oldepsb,oldepsg,oldmaxi ;
  local l0,u0,f0,prnt,header,kb,beta,nobs2,dispnam1 ;
  local rs,tmp1,tmp2,btmp,foo ;
  local rhotmp ;
  local hmtx,iter;
   local bbs0,loglik0,freeindx,cnstlist ;

  header=" RE serial model, concentrated loglikehood";

    if usestart== 0;
       " Loading data (RE serial model).. ";
       {smalldat,numpanls,nobs,xid,zid}=setiore(f1,parname1);
       if numpanls==0 ; goto notgot ; endif ;
    endif ;

     dispnam1="åîý"|"åàý"|"RHO_E"|"RY"|"LAM";
     btmp=varfe|vargrp|rho_e ;


@ Use this if you want maximum likelihood @
@ bbs= b varfe vargrp rho_e rho_y lambda @

     if usestart==1 ;
ask1: ;
             bbs=userset(dispnam1[1:3],b);?;?;?;
             if bbs[3]<0 or bbs[3]>1 ; goto ask1 ; endif ;
             " Loading data (RE serial model).. ";
             {smalldat,numpanls,nobs,xid,zid}=setiore(f1,parname1);
             if numpanls==0 ; goto notgot ; endif ;

    else ;                          @ use a RE regression @

         bbs=modelre2(f1,parname1,xid); format /lz 10,5 ;
         if scalmiss(bbs);
            goto notgot ;
         endif;
         bbs=bbs[rows(bbs)-1:rows(bbs)];
         output on;
         " Starting values (VARFE, VARGRP): " bbs';
         output off ;
         bbs=bbs|0;

    endif;


   let cnstlist[6,7]=-1 0  0     0 0     1 -1e-9
                      0 -1 0     0 0     1 -1e-9
                      0 0  -1    0 0     1 0
                      0 0  1     0 0     1  0.9999
                      0 0  0     1 0     2  0
                      0 0  0     0 1     2  1 ;
    if prntit==1 ;
         output on ;
         call shwconst(cnstlist,dispnam1);
         output off ;
         call keypause(10,-1) ;
    endif ;

    infinitp=1;

    {loglik,rs,freeindx}=maxlikc(bbs|0|0,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", all params",
                               &figrec,dispnam1);
    iter=varget("_mliter") ;

    if scalmiss(loglik) ;goto notgot ;  endif;


@ Compute final results, and output it @
@ Values of Beta and Var, given rho_e,rho_y and  lambda, were computed
  by figd2c -- pick em up @

      beta=bconc;
      b=beta|rs ;              @ the final solution @
      kb=rows(b);

      rhotmp=rs[3];
      b[kb-1]=0 ;            @ make sure rhoY is zero @
      b[kb]=1/(1-(rhotmp*rhotmp));    @ set correct lambda @

      FREEINDX=FREEINDX+ROWS(BETA);         @ MODIFY TO INCLUDE BETAS @
      FREEINDX=SEQA(1,1,ROWS(BETA))|FREEINDX ;
     " Computing covariance matrix ... ";
      {vc,hmtx}=compcvc(b,vctype,FREEINDX,&figre);

      kb=kb-2;       @ drop ry,lambda @
      b=b[1:kb] ;

      vc=vc[1:kb,1:kb] ;
      stderr=sqrt(diag(vc));
      tstat=b./stderr ;
      nobs2=nobs-numpanls-kb ;

      output on ;

      parnames=parname1|dispnam1[1:3] ;            @parnames is global @
      nlstprnt(_titles,dataset,depvar,nobs,nobs2,loglik,
              " Concentrated Loglik RE_Serial ",
               parnames,
               b,stderr,tstat,cdftc(abs(tstat),nobs2),hmtx,iter);
      ?;

      "  0 ó RHO_E < 1.0 " ;
      output off ;

     retp(1);

notgot: ;

     retp(0);

    endp;


/***************************/
@ do MLE of sporadic correlated irregular timed panel  -- concentrated
   likelihood.  Two variants: solve for all CV parameters simultaneously, of
solve singly and iterate.  For both variants, the betas are solved via GLS
using current values of rho_e,rho_y and lambda.
@

  proc modelrec(f1,parname1);

  local cv,logl,status,tmpe,tmpl,bbs,g1s,li,tmp,oldepsb,oldepsg,oldmaxi ;
  local l0,u0,f0,prnt,header,kb,nobs2,dispnam1 ;
  local rs,tmp1,tmp2,btmp,foo ;
  local hmtx,iter;
   local bbs0,loglik0,freeindx,cnstlist ;

    l0=miss(1,1); u0=l0; f0=u0;
    oldepsb=varget("_mlepsb") ;
    oldepsg=varget("_mlepsg") ;
    oldmaxi=varget("_maxiter");

    call varput(1.0,"_useallb");    @loosened for 1-by-1 @

    if usestart== 0;
       " Loading data (RE sporadic model).. ";
       {smalldat,numpanls,nobs,xid,zid}=setiore(f1,parname1);
       if numpanls==0 ; goto notgot ; endif ;
    endif ;

  header=" RE sporadic, Rho=P_e + P_y^Dt,  concentrated likelihood";


     dispnam1="åîý"|"åàý"|"RHO_E"|"RHO_Y"|"LAMBDA";
     btmp=varfe|vargrp|rho_e|rho_y|lambda ;


@ Use this if you want maximum likelihood @
@ bbs= b varfe vargrp rho_e rho_y lambda @

     if usestart==1 ;
ask1: ;
             bbs=userset(dispnam1,b);?;?;?;
             if scalmiss(yearid);
                bbs=bbs[1:3]|0|bbs[4] ;
             endif;
             tmp=bbs[3]+bbs[4] ;
             if tmp ge 1.0 ;
                " Inadmissible values ";
                goto ask1 ;
             endif;
             if infinitp==1 ;
               if (bbs[5]> ( 1/(1-tmp^2))) or (bbs[5]< ( 1/(1-(bbs[3]^2))) ) ;
                   " Internally inconsistent value for Lambda ";
                   " Min value = " (1/(1-bbs[3]^2)) ;
                   " Max value = " (1/(1-tmp^2)) ;
                   " Enter ESC to use'em anyways .... ";;
                  if keypause(4,27) ne 1; goto ask1; endif ;
               endif;
             endif ;
             " Loading data (RE sporadic model).. ";
             {smalldat,numpanls,nobs,xid,zid}=setiore(f1,parname1);
             if numpanls==0 ; goto notgot ; endif ;

    else ;                          @ use a RE regression @

         bbs=modelre2(f1,parname1,xid); format /lz 10,5 ;
         if scalmiss(bbs);
            goto notgot ;
         endif;
         bbs=bbs[rows(bbs)-1:rows(bbs)];
         output on;
         " Starting values (VARFE, VARGRP): " bbs';
         output off ;
         bbs=bbs|0|0|1;

    endif;


@ Either 1-by-1 or simultaneous ...... @

  if mleconc==1 ;

@ Iterated (1-by-1) concentrated likelihood uses two tricks:
  1) Compute beta's given varfe,vargrp,rho_E,rho_y, and lambda (using GLS)
  2) Estimate with MLE only one parameter at a time, and iterate until
    convergence criteria for all parameters is met.

   The first item is handled by figrec automatically (using compre2).
   The second item is handled below.
       The basic procedure is :
       1   Solve for rho_e (0<=rho_e<1)
       2   Solve for rho_y (0<=rho_y<(1-rho_e)
       3   Solve for lambda    (1/(rho_e+rho_Y)^2)<lambda<infinity.
       4   Solve for varfe
       5   Solve for vargrp
       6   Solve all 5_parameters
       6 Check full "gradient" at current values, and if convergence met exit.
           else goto step 1
@


      call varput(oldepsb*100,"_mlepsb");    @loosened for 1-by-1 @
      call varput(oldepsg*100,"_mlepsg");
      call varput(_max1itr,"_maxiter");   @ set in main program @


@  create bounds and do mle for each step  @
       bbs0=bbs ;         @ value before starting this step of iter@
       loglik0=-1e20 ; iter=0;

iter1: ;


@ STEP 1 -- solving for rho_e @
        if prntit==1 ;
           ?; output on ;  " Solving for rho_e " ; output off;
         endif;
         {cnstlist,bbs}=setulfre(bbs,3,infinitp,yearid,negok,dispnam1);
         foo=varput(2,"_mlatmax");    @ just 2 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", rho_e step ",
                               &figrec,dispnam1);
        b=bbs ;

        if scalmiss(loglik) ;
           " Hit SPACE to abort completely, ESC to try rho_y step ";;
           foo=1-keypause(12,-27);
           if foo==1 ;
              goto notgot ;
            endif;
        endif ;

        iter=varget("_mliter")+iter ;


@ STEP 2 -- solving for rho_y  -- assuming year is selected @
        if scalmiss(yearid); goto step3 ; endif ;  @ year not selected @
        if prntit==1 ;
           ?; output on ;  " Solving for rho_y " ; output off;
         endif;
         {cnstlist,bbs}=setulfre(bbs,4,infinitp,yearid,negok,dispnam1);
         foo=varput(2,"_mlatmax");    @ just 2 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", rho_y step ",
                               &figrec,dispnam1);
         b=bbs ;
        iter=varget("_mliter")+iter ;

        if scalmiss(loglik) ;
           " Hit SPACE to abort completely, ESC to try rho_y step ";;
           foo=1-keypause(12,-27);
           if foo==1 ;
              goto notgot ;
            endif;
        endif ;




@ STEP 3 -- solving for lambda @
step3: ;
       if prntit==1 ;
             ?;  output on ;  " Solving for lambda " ; output off ;
        endif;
         {cnstlist,bbs}=setulfre(bbs,5,infinitp,yearid,negok,dispnam1);
         foo=varput(2,"_mlatmax");    @ just 2 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", lambda step ",
                               &figrec,dispnam1);
        b=bbs ;
        iter=varget("_mliter")+iter ;
        if scalmiss(loglik) ;
           " Hit SPACE to abort completely, ESC to try rho_y step ";;
           foo=1-keypause(12,-27);
           if foo==1 ;
              goto notgot ;
            endif;
        endif ;


@ STEP 4 -- solving for varfe @
        if prntit==1 ;
           ?; output on ;  " Solving for varfe " ; output off;
         endif;
         {cnstlist,bbs}=setulfre(bbs,1,infinitp,yearid,negok,dispnam1);
         foo=varput(1,"_mlatmax");    @ just 1 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", varfe  step ",
                               &figrec,dispnam1);
        b=bbs ;
        iter=varget("_mliter")+iter ;

        if scalmiss(loglik) ;
           " Hit SPACE to abort completely, ESC to try rho_y step ";;
           foo=1-keypause(12,-27);
           if foo==1 ;
              goto notgot ;
            endif;
        endif ;

@ STEP 5-- solving for vargrp @
        if prntit==1 ;
           ?; output on ;  " Solving for vargrp " ; output off;
         endif;
         {cnstlist,bbs}=setulfre(bbs,2,infinitp,yearid,negok,dispnam1);
         foo=varput(1,"_mlatmax");    @ just 1 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", vargrp  step ",
                               &figrec,dispnam1);
        b=bbs ;
        iter=varget("_mliter")+iter ;
        if scalmiss(loglik) ;
           " Hit SPACE to abort completely, ESC to try all_5 step ";;
           foo=1-keypause(12,-27);
           if foo==1 ;
              goto notgot ;
            endif;
        endif ;

@ STEP 6 -- solve for all_5 @
          if prntit==1 ;
               ?; output on ;  " Solving for all_5 parameters  " ; output off ;
           endif;
           {cnstlist,bbs}=setulfre(bbs,0,infinitp,yearid,negok,dispnam1);
            foo=varput(2,"_mlatmax");    @ 2 criteria @
           {loglik,rs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", all_5 ",
                               &figrec,dispnam1);
           b=bbs ;
          iter=varget("_mliter")+iter ;

          if scalmiss(loglik) ;
              " Hit SPACE to abort completely, ESC to try rho_e step ";;
               foo=1-keypause(12,-27);
              if foo==1 ;
                 goto notgot ;
              endif;
          endif ;

@ Now, see if a signficant change has occured.  This will be the case if
  changes in b and loglik are sufficiently small @

        if iter>oldmaxi ;                       @ too many iterations @
           output on ; ? ;
           " MAXIMUM # of TOTAL (all steps) iterations exceeded = " oldmaxi ;
            " Current coefficient values will be used. ";
           output off ;
        else ;
           tmp1=abs(bbs-bbs0) ; tmp2=abs(loglik-loglik0);
           tmp1=tmp1./(abs(bbs0)+0.00001) ; tmp2=tmp2/(abs(loglik0)+0.00001);
           if (maxc(tmp1|tmp2)>oldepsb) ;            @ a significant change @
               bbs0=bbs ;
               loglik0=loglik ;
              goto iter1 ;
           endif ;                         @ Else, success @
         endif ;

@ if here, done estimating one at a time @


@ ___ NON ITERATED >>>> @
 else ;                 @ NON Iterated @
           {cnstlist,bbs}=setulfre(bbs,0,infinitp,yearid,negok,dispnam1);
           foo=varput(2,"_mlatmax");    @ 2 criteria @
           {loglik,rs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", all params",
                               &figrec,dispnam1);
          iter=varget("_mliter") ;
 endif;



 if scalmiss(loglik) ;goto notgot ;  endif;



@ Compute final results, and output it @
@ Values of Beta and Var, given rho_e,rho_y and  lambda, were computed
  by figd2c -- pick em up @

      beta=bconc;
      b=beta|rs ;              @ the final solution @
      kb=rows(b);

     " Computing covariance matrix ... ";
      {vc,hmtx}=compcvc(b,vctype,miss(1,1),&figre);
      stderr=sqrt(diag(vc));
      tstat=b./stderr ;
      nobs2=nobs-numpanls-kb ;

      output on ;

      parnames=parname1|dispnam1 ;            @parnames is global @
      nlstprnt(_titles,dataset,depvar,nobs,nobs2,loglik,
              " Concentrated Loglik RE_Sporadic ",
               parnames,
               b,stderr,tstat,cdftc(abs(tstat),nobs2),hmtx,iter);
      ?;

      if negok==1 ; " 0 ó RHO_E + RHO_Y < 1.0 " ;;
          else ; " 0 ó RHO_E, RHO_Y < 1.0 " ;;
      endif ;
      if infinitp==1 ; ", 1/(1-RHO_Eý) < Lambda < 1/(1-((RHO_E+RHO_Y)ý) ";
          else ; ", Lambda > 0 ";
      endif ;
      output off ;

     call varput(oldepsb,"_mlepsb");
     call varput(oldepsg,"_mlepsg");
     call varput(oldmaxi,"_maxiter");

beta=b ;  @ global var @
     retp(1);

notgot: ;

     call varput(oldepsb,"_mlepsb");
     call varput(oldepsg,"_mlepsg");
     call varput(oldmaxi,"_maxiter");
     retp(0);

    endp;

/***************************/
@ RE model. retp b|varfe|vargrp @

 proc modelre2(f1,parname1,xid);
    local xx,xy,yy,ysum,foo,x,y,nextobs,yrvals,i1,tranmtx,oldsize;
    local atline,plen,ith,xnew,ynew,bre;
    local kd,ssr,suminvt,bfe,varfe,varbet,vargrp,df,bbet,psi2,psi,t1 ;

    foo=seekr(f1,1);
    nextobs=1 ; oldsize=0;
    xx=0;yy=0;xy=0;ysum=0;foo=0;ith=0; suminvt=0;
     " Computing within estimates .. ";?;
    atline=csrlin-1 ;
    do until scalmiss(nextobs);

       if abortreg(ith,atline,1)==1;
           retp(miss(1,1));
       endif;

       {x,y,yrvals,nextobs}=grab_re(nextobs,smalldat);
       if scalmiss(y); continue; endif;
       ith=ith+1;
       plen=rows(y);
       if plen ne oldsize;
         i1=ones(plen,1);
         tranmtx=eye(plen)-(ones(plen,plen)/plen) ;
         oldsize=plen ;
        endif;
        xnew=tranmtx*x[.,xid] ; ynew=tranmtx*y ;
        {xx,xy,yy,ysum,foo}=makemom(xnew,ynew,xx,xy,yy,ysum,foo);
        suminvt=suminvt+ (1/plen);       @ used by RE @
    endo;

@ Compute Df @
    kd=rows(xid);
    df=nobs-kd-numpanls ;
    if df<1;
         errorlog " Not enough degrees of freedom .. ";
        retp(miss(1,1));
    endif;

    trap 1 ;
    bfe=solpd(xy,xx);
    trap 0 ;
    if scalmiss(bfe);
     errorlog
"ERROR: There is perfect collinearity among the independent variables,
and the moment matrix is singular.";
        call keypause(1.2,0);
        retp(miss(1,1));
     endif;

    output on ; " FE estimates: ";
    format 8,8; $parname1[xid]';
    format /lz 10,5 ; bfe' ;?; output off ;

     ssr=yy-bfe'xy;
     varfe=ssr/df;
     format /lo 15,9;
     output on ;
     "     Variance of iid portion of error, îit  =" varfe ; output off ;

    ? ; " Computing between estimates .. ";?;
    atline=csrlin-1 ;
@ Now do between estimator @
    foo=seekr(f1,1);
    nextobs=1 ; oldsize=0;
    xx=0;yy=0;xy=0;ysum=0;foo=0;ith=0;
    do until scalmiss(nextobs);

       if abortreg(ith,atline,1)==1;
           retp(miss(1,1));
       endif;

       {x,y,yrvals,nextobs}=grab_re(nextobs,smalldat);
       if scalmiss(y); continue; endif;
       ith=ith+1;
       plen=rows(y);
       if plen ne oldsize;
           tranmtx=ones(1,plen)./plen;     @reduces to 1 obs per panel@
           oldsize=plen ;
        endif;
        xnew=tranmtx*x ; ynew=tranmtx*y ;
        {xx,xy,yy,ysum,foo}=makemom(xnew,ynew,xx,xy,yy,ysum,foo);
    endo;

    bbet=solpd(xy,xx);
    trap 0 ;
    if scalmiss(bbet);
       errorlog
"ERROR: There is perfect collinearity among the between variables,
and the moment matrix is singular.";
        call keypause(1.2,0);
        retp(miss(1,1));
     endif;
     df=numpanls-rows(xx);
     ssr=yy-bbet'xy;
     varbet=ssr/df;
     t1=nobs/numpanls;            @ average panel size @
     format /lz 15,9 ;

    output on ;
     "      Variance of between estimator (VARBET) = "  varbet ;
     "      ä(1/Ti) = " suminvt ;;       @suminvt computed in FE @
             "    Average panel size= " t1 ;
    vargrp=varbet- (suminvt*varfe/numpanls) ;
    "      Variance of  random effect, ài, (VARGRP) = " vargrp ;
    if vargrp<=0;
           "               ::    VARGRP set to 0 ";
        vargrp=0.01;
    endif;

    psi=varfe/(varfe+(t1*vargrp));
  "         Psi==  Varfe/(Varfe+(T*Vargrp)) ; at average T=" psi;
  "                NB: Psi = 0 Í> FE    :: PSI= 1 Í> OLS   " ;
    output off ;

@ Note: for each panel, a term of the form:
         psi=varfe/(varfe + (plen*vargrp))
  Will be used to quasi demean the observation @

@ Do the Quasi demean (RE) regression, using results from above@

     ?;?;" Computing RE estimates ... ";
     atline=csrlin ;
     xx=0;xy=0;yy=0;nobs=0;ysum=0 ;      @ computed by makemom @
     nextobs=1; oldsize=0;
     foo=seekr(f1,1);    ith=0;
     do until scalmiss(nextobs);

          if abortreg(ith,atline,1)==1;
              retp(miss(1,1));
          endif;

@  Grab a panel. @
       {x,y,yrvals,nextobs}=grab_re(nextobs,smalldat);
       if scalmiss(y); continue; endif;
       ith=ith+1;
       plen=rows(y);

@ get transformation matrix, if new one needed@
       if plen/=oldsize ;                           @dont use prior matrix @
           psi=varfe/(varfe + (plen*vargrp));
           psi2=1-sqrt(psi) ;
           tranmtx=(eye(plen)-(psi2.*ones(plen,plen)./plen) );
           oldsize=plen;
        endif;                             @else, use prior value of tranmtx@
        xnew=tranmtx*x; ynew=tranmtx*y;

       {xx,xy,yy,ysum,nobs}=makemom(xnew,ynew,xx,xy,yy,ysum,nobs);
    endo;


@ Compute Df @
    kd=rows(parname1);
    df=nobs-kd ;
    if df<1;
         errorlog " Not enough degrees of freedom .. ";
        retp(miss(1,1));
    endif;

    trap 1 ;
    bre=solpd(xy,xx);
    trap 0 ;

    if scalmiss(bre);
        errorlog "ERROR: Quasi-demean step: moment matrix is singular";
        call keypause(1.2,0);
        retp(miss(1,1));
     endif;

     retp(bre|varfe|vargrp);

 endp;



/************************************/
@ compute  true log like for re model,
  given parameters b,rho_e,rho_y,lambda varfe vargrp @

    proc figlikre(f1,b,varfe,vargrp,rho_e,rho_y,lambda);

    local inpanel,cv,lik;
    local foo,atline,atcol,nextobs,dta ;
    local xnew,ynew,znew,cvdet,cvinv,yrvals,aconst,x,y;

    atline=csrlin;atcol=csrcol;  @for writing current status @


    nextobs=1;   lik=0;
    foo=seekr(f1,1);              @ rewind @
    do until scalmiss(nextobs);


@  Grab a panel. @

       {xnew,ynew,yrvals,nextobs}=grab_re(nextobs,smalldat);
       if scalmiss(ynew);
          continue ;
       endif ;

@ else, got a valid  panel @
       cv=mkcvre(varfe,vargrp,rho_e,rho_y,yrvals,lambda,miss(1,1));

       znew=ynew-(xnew*b);
       trap 1 ;
       cvinv=invpd(cv);
         if scalerr(cv);
         " Non psd CV "; cv ; stop ;
         endif;
       trap 0;
       cvdet=detl ;

       lik=lik+ (ln(cvdet) + (znew'cvinv*znew)) ;

    endo;
    lik=-lik/2 ;
    aconst=nobs*ln(1/sqrt(2*pi));
    retp(lik+aconst);

  endp;

@  Non-concentrated RE sporadic corr model.  Used only for hessian for VC
@

  proc (3)=figre(bs,narg);

    local lambda,ry,re,b,ls,varfe,vargrp,kb;
    local lisum,g1sum,h ;
 @ bs= b varfe vargrp re ry lambda  @
    kb=rows(bs);

    lambda=bs[kb] ;
    ry=bs[kb-1] ;
    re=bs[kb-2];
    vargrp=bs[kb-3];
    varfe=bs[kb-4];
    b=bs[1:kb-5];


     lisum=0; h=zeros(kb,kb); g1sum=0;

    if narg==3  ;           @ return -loglik @
       lisum=figlikre(f1,b,varfe,vargrp,re,ry,lambda);

    elseif narg==1 or narg==2 ;
      " Hessian unavailable ." ; stop ;

    else ;                     @ 0 ,4 or 5 @
       {lisum,g1sum,h}=figreb(b,varfe,vargrp,re,ry,lambda,narg);
    endif;
     retp(lisum,g1sum,h);
  endp;



/************************************/
@ compute log like, gradient, and g'g  for re sporadic corr model,
  given parameters b,rho_e,rho_y,lambda,varfe,vargrp
  Returns li,g,g'g -- dependent on narg:
    narg=0  compute all 3
         4  G'G only
         5  Li and G

@

    proc (3)=figreb(b,varfe,vargrp,rho_e,rho_y,lambda,narg);

    local inpanel,cv,lik;
    local foo,atline,atcol,nextobs,dta,bbs;
    local xnew,ynew,z,cvdet,cvinv,yrvals,var;
    local dldvarfe,dldvargr,dldre,dldry,dcvdvarf,dcvdvarg;
    local dldlam,dldb,dcvdvar,dcvdlam,dcvdre,dcvdry;
    local g1,gtot,ggtot ;

    atline=csrlin;atcol=csrcol;  @for writing current status @

    dldvarfe=0;dldvargr=0;dldb=0;dldlam=0;dldre=0;dldry=0;

    nextobs=1;
    gtot=0;ggtot=0;lik=0;

    foo=seekr(f1,1);              @ rewind @
    do until eof(f1) or scalmiss(nextobs) ;

@  Grab a panel. @

       {xnew,ynew,yrvals,nextobs}=grab_re(nextobs,smalldat);

       if scalmiss(ynew); continue ; endif ;

@ else, okay panel @

      {dcvdvarf,dcvdvarg,dcvdre,dcvdry,dcvdlam,cv}=
         gradcvre(rho_e,rho_y,yrvals,lambda,varfe,vargrp);   @ deriv of cv @

       z=ynew-(xnew*b);

      trap 1;
      cvinv=invpd(cv);
      if scalerr(cvinv);
          " NON-psd CV ";
          stop ;
      endif;
      trap 0 ;
      if narg ne 4 ;            @  compute likelihood @
           cvdet=detl ;
           lik=lik+(ln(cvdet) + (z'cvinv*z)) ;    @ - loglik @
      endif;

   dldvarfe=dldvarfe+ (sumc(diag(cvinv*dcvdvarf)) -
                       (z'(cvinv*dcvdvarf*cvinv)*z) );
   dldvargr=dldvargr + (sumc(diag(cvinv*dcvdvarg)) -
                    (z'(cvinv*dcvdvarg*cvinv)*z) );
   dldlam=dldlam+ (sumc(diag(cvinv*dcvdlam)) - (z'(cvinv*dcvdlam*cvinv)*z) ) ;
   dldre=dldre +  (sumc(diag(cvinv*dcvdre)) - (z'(cvinv*dcvdre*cvinv)*z) )  ;
   dldry=dldry +  (sumc(diag(cvinv*dcvdry)) - (z'(cvinv*dcvdry*cvinv)*z) )   ;

       dldb=dldb-( (2*cvinv*z)'xnew) ;

 @ bs= b varfe vargrp re ry lambda  @
       g1=dldb~dldvarfe~dldvargr~dldre~dldry~dldlam ; g1=-g1 ;

       if narg ne 4 ;
          gtot=gtot+g1;
       endif ;

       if narg ne 5 ;
          ggtot=ggtot+(g1'g1) ;
       endif;
    endo ;

    if narg ne 4 ;
        lik=(-lik/2)+ ( (nobs)*ln(1/sqrt(2*pi)) );
    endif;
    retp(lik,gtot',ggtot);

  endp;

/*****************************/
@ concentrated log likelihood, re  sporadic correlated model
 all 5 params at once @

  proc (3)=figrec(bs,narg);

    local lambda,ry,re,b,ls,kb,varfe,vargrp;
    local lisum,g1sum,h,tmp ;

 @ bs= b varfe vargrp re ry lambda @
    kb=rows(bs);
    lambda=bs[5] ;
    ry=bs[4] ;
    re=bs[3];
    varfe=bs[1];
    vargrp=bs[2];
     lisum=0; h=zeros(kb,kb); g1sum=0;

    if narg==3  ;
      {lisum,b}=compre2(varfe,vargrp,re,ry,lambda,1);
      if scalmiss(lisum) ; retp(lisum,0,0); endif;

     elseif (narg==1) or (narg==2) ;
        " Hessian unavailable ." ;
        stop ;

    else ;                     @ 0 ,4 or 5 @

       {lisum,g1sum,h}=figrecb(varfe,vargrp,re,ry,lambda,narg);

    endif;
    retp(lisum,g1sum,h);
  endp;



/************************************/
@ compute log like, gradient, and g'g  for concentrated RE
  sporadic and serial model,
  given parameters varfe,vargrp, rho_e,rho_y,lambda.

  ALSO -- set global parameters Bconc = beta
  Returns li,g,g'g -- dependent on narg:
    narg=0  compute all 3
         4  G'G only
         5  Li and G

  Must make 2 passes through data. First pass to compute beta,
  and the second pass to compute gradients wrt
  varfe vargrp rho_e,rho_y and lambda.

@

    proc (3)=figrecb(varfe,vargrp,rho_e,rho_y,lambda,narg);

    local cv,lik,foo;
    local atline,atcol,nextobs,dta,bbs;
    local nmtx,xnew,ynew,z,cvdet,cvinv,yrvals,var;
local dldre,dldry,dldlam,dldvarf,dldvarg ;
local dcvdvarf,dcvdvarg,dcvdlam,dcvdre,dcvdry;
    local g1,gtot,ggtot ;
    local beta ;

    atline=csrlin;atcol=csrcol;  @for writing current status @

    dldlam=0;dldre=0;dldry=0;dldvarf=0;dldvarg=0;

    nextobs=1;
    gtot=0;ggtot=0;lik=0;


@ first pass -------- @

    {foo,beta}=compre2(varfe,vargrp,rho_e,rho_y,lambda,0);  @ 0==ignore loglik
@


@ second pass -- compute likelihood and gradient. Note we don't compute
  likelihood above simply to save time @

    foo=seekr(f1,1);                  @ rewind @
    do until eof(f1) or scalmiss(nextobs) ;

@  Grab a panel. @

       {xnew,ynew,yrvals,nextobs}=grab_re(nextobs,smalldat);

       if scalmiss(ynew); continue ; endif ;

@ else, okay panel.  Build grad, etc. a panel at a time @

      {dcvdvarf,dcvdvarg,dcvdre,dcvdry,dcvdlam,cv}=
         gradcvre(rho_e,rho_y,yrvals,lambda,varfe,vargrp);   @ deriv of cv @

      z=ynew-(xnew*beta);
      trap 1;
      cvinv=invpd(cv);
      if scalerr(cvinv);
          " NON-psd CV ";
          stop ;
      endif;
      trap 0 ;
      if narg ne 4 ;            @  compute likelihood @
           cvdet=detl ;
           lik=lik+(ln(cvdet) + (z'cvinv*z)) ;    @ - loglik @
      endif;

   dldlam=dldlam+ (sumc(diag(cvinv*dcvdlam)) - (z'(cvinv*dcvdlam*cvinv)*z) ) ;
   dldre=dldre +  (sumc(diag(cvinv*dcvdre)) - (z'(cvinv*dcvdre*cvinv)*z) )  ;
   dldry=dldry +  (sumc(diag(cvinv*dcvdry)) - (z'(cvinv*dcvdry*cvinv)*z) )   ;
   dldvarf=dldvarf +  (sumc(diag(cvinv*dcvdvarf))
                - (z'(cvinv*dcvdvarf*cvinv)*z) )   ;
   dldvarg=dldvarg +  (sumc(diag(cvinv*dcvdvarg))
                - (z'(cvinv*dcvdvarg*cvinv)*z) )   ;

 @ concentrated bs= varfe vargrp re ry lambda  @
       g1=dldvarf~dldvarg~dldre~dldry~dldlam ;
       g1=-g1 ;    @ true grad @

       if narg ne 4 ;
          gtot=gtot+g1;
       endif ;

       if narg ne 5 ;
          ggtot=ggtot+(g1'g1) ;
       endif;
    endo ;

    if narg ne 4 ;
      lik= -lik/2 ;            @  variable part of true likelihood @
      lik=lik + ( (nobs)*ln(1/sqrt(2*pi)) );
    endif;

@ set global pecuiliar to this concentrated likelihood function @
    bconc=beta ;

    retp(lik,gtot',ggtot);

  endp;

/******************************/
@ set bounds for RE model
 Note: always do as a concentrated likelihood (perhaps as iterated conc)

            {cnstlist,bbs}=setulfre(bbs,1,infinitp,yearid,negok,pnames);
@

 proc (2)=setulfre(beta,usec,infinitp,yearid,negok,pnames);

      local ff,tmp,cnstlist,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp9 ;

@ Vars are: variid,vargrp,rhoe,rhoy,lambda @

@ note: the idea is that 1,1 element of CV matrix in sporadic correlation
        equals 1/(1-p**2).  Given 0<p<1, this fraction is at least 1.
        Furthermore, if p=rhoe+rhoy**dt, then dt=1 yields max value of p=
        rhoe+rhoy, so fraction is maximized at 1/(1-(rhoe+rhoy)**2).
        However, if infinitp==0, then we ARE NOT assuming that the error
        process was generate in this fashion, so we allow lambda to be
        anything (> 0 and < infinity.

        Note: in unconstrained model, we are unable to set max value of
        lambda using linear constraints.  We hope it works.
@

@ used below @
      if negok==1 ;        @ constrain 0<rhoe+rhoy<1 @
            let tmp5=0 0 -1 0 0 1 -0.999 ; @ rhoe > -1 @
            let tmp6=0 0 0 -1 0 1 -0.999 ;   @ rhoy > -1 @
       else;               @ negative rhoe,y not allowed @
            let tmp5=0 0 -1 0 0 1 0 ; @ rhoe > 0 @
            let tmp6=0 0 0 -1 0 1 0 ;   @ rhoy >0 @
      endif ;

@ Not 1-by-1 ???? @
   if usec==0 ;             @ fully free @
@ First, variid and vargrp are always > 0, lambda > 0 @
        let cnstlist[3,7]=-1 0 0 0 0 1 -1e-9
                           0 -1 0 0 0 1 -1e-9
                           0 0 0 0 -1 1 -1e-9 ;
        if infinitp==1 ; cnstlist[3,7]=-1 ; endif ;  @ lambda gt 1 @

         let tmp1=0 0 1 1 0 1 0.999 ; @ re+ry<1 @
         let tmp2=0 0 -1 -1 0 1 0 ;    @ re+ry>0 @
         let tmp3=0 0 1 0 0 1 0.999 ; @ rhoe < 1 @
         let tmp4=0 0 0 1 0 1 0.999 ;   @ rhoy < 1 @
         cnstlist=cnstlist|tmp1'|tmp2'|tmp3'|tmp4'|tmp5'|tmp6' ;

@ else, 1-by-1 concentrated logl @
   else ;
@ Constrain all variables to be their current values @
       cnstlist=eye(5);
       cnstlist=cnstlist~(2*ones(5,1));
       cnstlist=cnstlist~beta;
       ff=usec;
       if ff==1 ;                    @ varfe free @
           cnstlist[1,1]=-1 ; cnstlist[1,6]=1 ; cnstlist[1,7]=-1e-9 ;
       elseif ff==2 ;
           cnstlist[2,2]=-1; cnstlist[2,6]=1 ; cnstlist[2,7]=-1e-9 ;
       elseif ff==3 ;                @ only rho_e param free @
             let tmp3=0 0 1 0 0 1 0 ; tmp3[7]=0.9999-beta[4]; @ rhoe+rhoy<1@
             cnstlist[3,.]=tmp3' ;
             tmp=0.9999-beta[4] ;
             beta[3]=minc(beta[3]|tmp);
             cnstlist=cnstlist|tmp5' ;
       elseif ff==4 ;           @ only rho_y free @
             let tmp3=0 0 0 1 0 1 0 ; tmp3[7]=0.9999-beta[3]; @ rhoe+rhoy<1@
             cnstlist[4,.]=tmp3' ;
             tmp=0.99999-beta[3]  ;
             beta[4]=minc(beta[4]|tmp);
             cnstlist=cnstlist|tmp6' ;
       elseif ff==5 ;           @ only lambda free @
          let tmp3=0 0 0 0 1 1 1e26 ;   @ lambda lt infinity @
          let tmp9=0 0 0 0 -1 1 -1e-9 ;    @ lambda > 0 @
          if infinitp==1 ;        @ process begins at time -infinity @
             tmp=beta[3]+beta[4] ;       @which implies bounds on lambda@
             tmp= ( 1/(1-tmp^2)) ;
             tmp3[7]=tmp ;                 @ upper bound not infinity @
             cnstlist[5,.]=tmp3' ;

             beta[5]=minc(beta[5]|(tmp-1e6));

             tmp= (1/(1-(beta[3]^2)));     @ lower bound as if rhoy=0 @
             tmp3[7]=-tmp ;
             tmp3[5]=-1 ;
             cnstlist=cnstlist|(tmp3')|(tmp9') ;

             beta[5]=maxc(beta[5]|tmp);

          else ;
             cnstlist[5,.]=tmp3' ;
             cnstlist=cnstlist|(tmp9') ;
          endif;
       endif;


   endif;         @ 1 param at a time @

   if prntit==1 ;
         output on ;
           " beta " beta' ;
           call shwconst(cnstlist,pnames);
         output off ;
         call keypause(10,-1) ;
   endif ;

   retp(cnstlist,beta);

 endp;

/***************************/
 proc (4)=grab_re(nxt,smalldat);
   local dta,xt,yt,yrs,nxtpnl,inpanel,pat,a1,tmp ;

    if smalldat==1 ;                            @smalldat is global @
        nxtpnl=nxt;
        pat=panelens[nxtpnl,.];
        a1=bigdata[pat[1]:pat[2],.] ;
        tmp=cols(a1);
        xt=a1[.,1:tmp-2] ;
        yt=a1[.,tmp-1];
        yrs=a1[.,tmp];
        if nxtpnl==rows(panelens) ;
             nxtpnl=miss(1,1);
        else;
             nxtpnl=nxtpnl+1 ;
        endif;
        retp(xt,yt,yrs,nxtpnl);

     else ;                             @ else, read from datafile @
         {y,x,yearvals,nxt}=grabpanl(f1,nxt);   @f1 is global @
         retp(x,y,yearvals,nxt);             @ x,y,yearvals are global @
     endif;

endp;




/*********************************/
@ see if we can fit it all into memory, for RE.  @
@ retp smalldat,numpanls,nobs,xid,zid @
   proc (5)=setiore(f1,b);

    local nk,vcmem,memavail,nr,vn,aa;
    local ynew,xnew,yrvals,nextobs,amiss,tmp,iat,ry ;
    local novary,smalldt0 ;

    amiss=miss(1,1);
    novary=1;                           @ for pointing to panel invariant vars@

    nk=rows(b)+4;
    vcmem=nk*nk*5*8;
    memavail=coreleft-vcmem;

    if memavail>500000 ;
       nr = ceil(minc(
               coreleft/((nk+1)*8*5)|100000/((nk+1)+1)));
@nb: use 100000 instead of maxvec ! @

    else ;                  @ probably 2.0 @

       vn=minc(64*5|(memavail/1024));
       nr=floor(((vn/5)*128-2)/(nk+1));

    endif;
    if nr < 1;
 " ERROR: Insufficient memory  -- reduce # of variables in dataset. " ;
      retp(0,0,0,0,0) ;
   elseif nr < rowsf(f1);
       smalldt0=0;
   else ;
       smalldt0=1 ;
   endif;

@ Get some bookkeeping info, and load into memory if possiblwe @

    nextobs=1;bigdata=amiss ;iat=1;
    vn=seekr(f1,1); numpanls=0; nobs=0; ?; aa=csrlin-1;
    do until scalmiss(nextobs);

       {ynew,xnew,yrvals,nextobs}=grabpanl(f1,nextobs);
       if scalmiss(ynew);
            continue ;          @ not a valid panel @
       else;                  @ store in bigdata @
           numpanls=numpanls+1;         @ GLOBAL @
           if (Numpanls % 25)==1 ;
              locate aa,1; format 4,0; numpanls " read .. ";
           endif;
           nobs=nobs+rows(yrvals);      @ GLOBAL @
           novary=novary.* dotfeq(minc(xnew),maxc(xnew)) ;

            if smalldt0==1 ;
               ry=rows(yrvals)-1;
               if scalmiss(bigdata);
                   bigdata=xnew~ynew~yrvals; panelens=iat~(iat+ry);
               else;
                   bigdata=bigdata|(xnew~ynew~yrvals) ;
                   panelens=panelens|(iat~(iat+ry));
               endif;
               iat=iat+ry+1;
            endif;                    @ scalmiss(bigdata) @
       endif;                          @ smalldat0=1 @
    endo;                       @ through dataset @

   {xid,zid}=getxzid(novary);       @ GLOBALS @

   if scalmiss(xid);
        " No panel variant variables. "; call keypause(1.2,0);
    endif;
    if not(scalmiss(zid));
      " Note: the following are panel invariant:";; format 1,8;
        $parname1[zid] ; call keypause(1.4,0);
    endif;

  retp(smalldt0,numpanls,nobs,xid,zid);

endp;


/******************************************/
@ create the CV matrix for the RE sporadic model
User supplies:
  varfe  = IID component of variancec
  vargrp = RE component of variance
  rho_e = event correlation coefficient
  rho_y = year correlation coefficient)
  years     = Tx1 vector of 'years' of event occurences
  lambda = the past history correction factor

Basically two components:

C = i*i' * vargrp + varfe*E(U*U')  -- with U the (normalized) serial
correleated component

@

  proc mkcvre(varfe,vargrp,re,ry,years,lambda,p);
   local
      cv, tr, dt, cvmtx, psq, t, pr1, tmp1, n, s, pr2, k ;


@ is it serial (not sporadic) model, or a special case of sporadic ? @
   if (ry==0 and infinitp==1) ;         @ use simple form @
           cv=mkcvrea(varfe,vargrp,re,years);
           retp(cv);
   endif;

   Tr=rows(years);

@ special case @
     if tr==1 ;
       retp(vargrp+(lambda*varfe));
     endif;

@ vector of P correlation coefficients, and related values @
@ note: may already have been computed @
     if scalmiss(p);
       dt=years[2:tr]-years[1:tr-1] ;
       dt=dt+(dt.==0);                  @ set minimum dt to 1 @
       P=re + (ry.^dt) ;
       p=1|p ;   @ use any arbitrary value, since year[0] is unknown @
    endif ;

@ another special case @
    if tr==2 ;
       cvmtx=eye(2);
       cvmtx[1,1]=lambda ;
       cvmtx[2,2]=(p[2]*p[2]*lambda) + 1 ;
       cvmtx[1,2]=p[2]*lambda;
       cvmtx[2,1]=cvmtx[1,2];
       retp((varfe*cvmtx)+vargrp);
    endif;

   psq=p.*p ;           @ used below @

@ not special case @
   cv=zeros(tr,tr);

   cv[1,1]=lambda ;                     @ by definition @

@ compute diagonals @
   t=1;
   do until t==tr ;

     t=t+1;

     pr1=prodc(psq[2:t])*lambda;

     tmp1=1+pr1 ;

     if t>2 ;
        n=0;
        do until n gt t-3 ;
          tmp1=tmp1+prodc(psq[(t-n):t]);
          n=n+1;
        endo ;
     endif ;

     cv[t,t]=tmp1 ;

   endo;

@ now do off diagonals @
  t=2 ;
  do until t gt tr ;                  @ t is rows, > s @
    s=1 ;
    do until s eq t ;                 @ s is columns @

@ compute outer term @
        pr2=prodc(p[(s+1):t]);

@ compute inner term @
        if s==1 ;
            pr1=1 ;
        else ;
           pr1=prodc(psq[2:s]);
        endif ;
        tmp1=(pr1*lambda) ;

        if s>1 ;
          tmp1=tmp1+1 ;
          k=3 ;
          do until k>s ;
            tmp1=tmp1+prodc(psq[k:s]);
            k=k+1 ;
          endo;
        endif;
        cv[t,s]=pr2*tmp1 ;
        cv[s,t]=cv[t,s];
        s=s+1;
    endo;                      @ s @
    t=t+1;
  endo;                        @ t @

@ incorporate variance components @
  cv=(cv*varfe)+vargrp ;

  retp(cv);

endp ;


/****************************/
@ make cv for RE serial corr model when ry=0 @

   proc mkcvrea(varfe,vargrp,rho,years);

   local tmp,cv,cv1,i,tmp1,ny;

   ny=rows(years);              @ years used to determine panel size @
   tmp1=seqa(0,1,ny);
   cv1=(rho.^tmp1)' ; cv=cv1 ;
   i=1 ;
   do until i> (ny-1);
     tmp=shiftr(cv1,i,0);
     cv=cv|tmp;
     i=i+1;
   endo;
   tmp=zeros(ny,1);
   tmp=diagrv(cv,tmp);
   cv=cv+tmp';
   tmp=varfe/(1-rho^2) ;
   retp((cv*tmp)+vargrp);

endp;

/*******************/
   proc mkdcva(rho,years);

   local tmp,cv,cv1,i,tmp1,ny;

   ny=rows(years);
   tmp1=seqa(0,1,ny);
   if rho==0 ;
      cv1=zeros(1,ny);cv1[2]=1;
   else ;
      cv1=tmp1'.*(rho.^(tmp1-1))' ;
   endif;
   cv=cv1 ;
   i=1 ;
   do until i> (ny-1);
     tmp=shiftr(cv1,i,0);
     cv=cv|tmp;
     i=i+1;
   endo;
   tmp=zeros(ny,1);
   tmp=diagrv(cv,tmp);
   cv=cv+tmp' ;
   retp(cv);

endp;


/***************************/
@{dcvdvarf,dcvdvarg,dcvdre,dcvdry,dcvdlam,cv}=@
proc (6)= gradcvre(rho_e,rho_y,years,lambda,varfe,vargrp);
  local
      tr, dcvdre, dcvdry, dcvdlam, cv, dt, p, dcvdvarg, dcvdvarf, amtx,
      dlam, dcva, psq, dpdry, dpdre, dpsqdre, dpsqdry, t, trm1, dtrm1,
      dtrm1y, s, term2, dterm2, dterm2y, term1, dterm1, dterm1y, k, cvtmp,
      n ;

@ Initialize some vars @

    Tr=rows(years);
    dcvdre=zeros(tr,tr); dcvdry=dcvdre ; dcvdlam=dcvdre ; cv=dcvdre ;


    dt=years[2:tr]-years[1:tr-1] ;
    dt=dt+(dt .== 0);           @ minimum of 1 year @
    p = rho_e + ( rho_y.^dt) ;
    p=1|p ;   @ use any arbitrary value, since year[0] is unknown @

@ Check for special case of ry=0 or serial (not sporadic) correlation @
    if (rho_y==0 and infinitp==1) ; @ use simple form @

          dcvdvarg=ones(tr,tr);

@ COMPUTE deriv wrt to varfe @
          dcvdvarf=mkcvrea(1,0,rho_e,years);

@ COMPUTE CV @
          cv=(dcvdvarf*varfe)+vargrp ;

@ dcvdry and dcvdlam are non-existent (so we might as well return zeros) @
@ Compute dcvdre here @
            lambda=1/(1-(rho_e*rho_E)) ;

            amtx=dcvdvarf.*(1/lambda) ;
            dlam=(2*rho_e)/((1-(rho_e*rho_e))^2) ;
            dcva=mkdcva(rho_e,years);
            dcvdre=((lambda.*dcva)+(amtx.*dlam))*varfe;

            retp(dcvdvarf,dcvdvarg,dcvdre,dcvdry,dcvdlam,cv);
    endif ;



@ if here, not simple case @
    cv[1,1]=lambda ;
    dcvdlam[1,1]=1;

    psq=p.*p ;
    dpdry=dt.*(rho_y.^(dt-1)) ; dpdry=1|dpdry ;

    dpdre=ones(tr,1);

    dpsqdre=2*p ;
    dpsqdry=2*dpdry.*p;

    t=2 ;
    do until t > tr ;

@ s=1 @
       {trm1,dtrm1,dtrm1y}=dprodx(p,dpdre,dpdry,2,t);
       dcvdre[t,1]=lambda*dtrm1;
       dcvdre[1,t]=dcvdre[t,1];

       dcvdry[t,1]=lambda*dtrm1y;
       dcvdry[1,t]=dcvdry[t,1];

       dcvdlam[t,1]=trm1 ; dcvdlam[1,t]=dcvdlam[t,1] ;

       cv[t,1]=lambda*trm1 ; cv[1,t]=cv[t,1] ;

       s=2 ;

       do until s==t ;

@  d term 2 @
            {term2,dterm2,dterm2y}=dprodx(p,dpdre,,dpdry,s+1,t);

 @ d term1 @
           {term1,dterm1,dterm1y}=dprodx(psq,dpsqdre,dpsqdry,2,s);

            dcvdlam[t,s]=term1*term2 ;
            dcvdlam[s,t]=dcvdlam[t,s] ;

             term1=(lambda*term1)+1;

             dterm1=lambda*dterm1;
             dterm1y=lambda*dterm1y;

             k=3 ;
             do until k> s ;
                {trm1,dtrm1,dtrm1y}=dprodx(psq,dpsqdre,dpsqdry,k,s);
                term1=term1+trm1; dterm1=dterm1+dtrm1 ;
                dterm1y=dterm1y+dtrm1y ;
                k=k+1;
             endo;

@ d cv = term1 * dterm2 + dterm1 * term2 @

             dcvdry[s,t]=(term1*dterm2y)+(dterm1y*term2);
             dcvdry[t,s]=dcvdry[s,t];

             dcvdre[s,t]=(term1*dterm2)+(dterm1*term2);
             dcvdre[t,s]=dcvdre[s,t];

             cv[s,t]=term1*term2 ; cv[t,s]=cv[s,t] ;

             s=s+1;
          endo ;                @ IS @

@ diagonal term @
          {term1,dterm1,dterm1y}=dprodx(psq,dpsqdre,dpsqdry,2,t);
          dcvdlam[t,t]=term1 ;
          dterm1=dterm1*lambda ;
          dterm1y=dterm1y*lambda ;
          cvtmp=(term1*lambda)+1 ;
          n=0;
          do until n gt t-3 ;
            {trm1,dtrm1,dtrm1y}=dprodx(psq,dpsqdre,dpsqdry,t-n,t);
            cvtmp=cvtmp+trm1 ;
            dterm1=dterm1+dtrm1 ;
            dterm1y=dterm1y+dtrm1y ;
            n=n+1;
          endo;

          dcvdre[t,t]=dterm1 ;
          dcvdry[t,t]=dterm1y ;
          cv[t,t]=cvtmp ;

          t=t+1;
      endo ;

@ DcvDvarf and CV matrix @
    dcvdvarf=cv ;
    cv=(cv*varfe)+vargrp ;

@ COMPUTE deriv wrt vargrp  (dcvdvarg is ones matrix @
   dcvdvarg=ones(tr,tr);

@{dcvdvarf,dcvdvarg,dcvdre,dcvdry,dcvdlam,cv}=@
retp(dcvdvarf,dcvdvarg,dcvdre*varfe,dcvdry*varfe,dcvdlam*varfe,cv);

endp;



/************************************/
@ GIven varfe,vargrp,re,ry,lambda -- compute beta that minimizes log
like.  For sporadic and serial correlated RE model
 Use global var f1
    { lik,b}=compre2(varfe,vargrp,re,ry,lambda,dolik)
@

   proc (2)=compre2(varfe,vargrp,re,ry,lambda,dolik);


    local cv,lik;
    local foo,npanels,atline,atcol,nextobs,dta ;
    local xnew,ynew,znew,cvdet,cvinv,tmp,yrvals;
    local xcvx,xcvy,ycvy ;
    local sumcvdet,sumzcvz,bgls ;

    atline=csrlin;atcol=csrcol;  @for writing current status @

    nextobs=1;
    sumcvdet=0;xcvx=0;xcvy=0; ycvy=0;

    foo=seekr(f1,1);              @ rewind @
    do until eof(f1) or scalmiss(nextobs);

       {xnew,ynew,yrvals,nextobs}=grab_re(nextobs,smalldat);
         if scalmiss(ynew);
            continue ;
         endif;

@ else, got a panel @

@ get VC matrix @
        cv=mkcvre(varfe,vargrp,re,ry,yrvals,lambda,miss(1,1));

        trap 1 ;
        cvinv=invpd(cv);
        if scalerr(cvinv);
output on ; ? ;?;
           " Problem with invpd(cv) ,varfe,vargrp,re,ry,lambda" ;
           varfe vargrp re ry lambda ; yrvals';
?; output off ;
            goto aborted ;
            cvinv=inv(cv);
        endif;
        trap 0;
        xcvx=xcvx+ (xnew'cvinv*xnew);
        xcvy=xcvy+(xnew'cvinv*ynew);
        ycvy=ycvy+(ynew'cvinv*ynew);

        cvdet=detl ;
        sumcvdet=sumcvdet+ln(cvdet);
    endo;

@ comute b gls @
clearg fooxcvx ; fooxcvx=xcvx ;
    bgls=invpd(xcvx)*xcvy ;

@ and lik @
    if dolik==1 ;
       lik=figlikre(f1,bgls,varfe,vargrp,re,ry,lambda);
    else;
      lik=miss(1,1);
    endif;

    retp(lik,bgls);                @ return loglik, etc. @

aborted: ;
   retp(miss(1,1),0);

  endp;






