

/********************************************/
@  quick OLS on first differenced stuff @
@         {b,varhat}=modelD2a(f1) ;@

  proc (2)=modeld2a(f1);

    local nextobs,foo,kd;
    local ynew,xnew,oldsize,yrvals,tranmtx,nobs0 ;

    xx=0;xy=0;yy=0;nobs0=0;ysum=0 ;      @ computed by makemom @
    totobs=rowsf(f1);

    nextobs=1; tranmtx=0;oldsize=0;

    foo=seekr(f1,1);              @ rewind @

    do until scalmiss(nextobs) ;

@  Grab a panel. @
        {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
             grab_d2(oldsize,nextobs,tranmtx,xid);

        if scalmiss(y);                  @ if empty, skip to next panel @
           continue ;
        endif;

@ Add to moment,  etc. matrices.  @

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

    endo;

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

@ compute b, other stuff (returned to global vars ) @
     if compols(xx,xy,yy,ysum,nobs,df)==0;
        retp(0,0);
     else ;
        retp(b,varhat);
     endif;

 endp;


/***************************/
@ do MLE of serial 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.

 NOte: should never get here if grid search is requested.

@

  proc modeld2c(f1,parname1);

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

    l0=miss(1,1); u0=l0; f0=u0;

    header=" 1st diff, Rho = P_e + P_y^Dt, model, concentrated loglikehood ";

    dispname="RHO_E"|"RHO_Y"|"LAMBDA" ;
    if not(scalmiss(yearid));
       dispnam1=dispname ;
       btmp=rho_e|rho_y|lambda ;
    else ;
       dispnam1="RHO_E"|"LAMBDA" ;
       btmp=rho_e|lambda ;
    endif;

    smalldat=0;                       @force read of file @



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

        if usestart==1 ;
ask1: ;
             bbs=userset(dispnam1,btmp);?;?;?;
             if scalmiss(yearid);
                bbs=bbs[1]|0|bbs[2] ;
             endif;
             tmp=bbs[1]+bbs[2] ;
             if tmp ge 1.0;
                " Inadmissible values ";
                goto ask1 ;
             endif;
             if infinitp==1 ;
              if (bbs[3]> ( 1/(1-tmp^2))) or (bbs[3]< ( 1/(1-(bbs[1]^2))) ) ;
                  " Internally inconsistent value for Lambda ";
                  " Min value = " (1/(1-bbs[1]^2)) ;
                  " Max value = " (1/(1-tmp^2)) ;
                  " Enter ESC to use'em anyways .... ";;
                 if keypause(2,27) ne 1; goto ask1; endif ;
              endif;
             endif ;
        else ;
            let bbs=0 0 1  ;  @ use arbitrary values @
        endif;

        " Loading data ..(concentrated 1st diff)... ";
        {smalldat,numpanls,nobs,xid,zid}=setiod2(f1,parname1);
        if numpanls==0 ; retp(0)    ; endif ;
        if scalmiss(xid); retp(0)     ; endif ;

@ Iterated (1-by-1) concentrated likelihood uses two tricks:
  1) Compute beta's given 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 figd2c automatically (using compbd2).
   The second item is handled below.
       The basic procedure is :
       1   Set rho_y and lambda to prior values
              Solve for rho_e (0<=rho_e<1)
       2   Set rho_e and lambda to prior values (with rho_e prior value being
           the number computed just above), and solve for rho_y
               (0<=rho_y<(1-rho_e)
       3   Set rho_e and rho_y, solve for lambda
               (1/(rho_e+rho_Y)^2)<lambda<infinity.
       4   Solve all parameters simultantously.
       4 Check change in coeffs and loglik, and if convergence met exit;
           else goto step 1
@

oldepsb=varget("_mlepsb") ;
oldepsg=varget("_mlepsg") ;
oldmaxi=varget("_maxiter");


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

   if mleconc==1 ;                     @ 1-by-1 concentrated @

       call varput(_max1itr,"_maxiter");
       call varput(oldepsb*100,"_mlepsb");
       call varput(oldepsg*100,"_mlepsg");

@ STEP 1 -- solving for rho_e @
        if prntit==1 ;
           ?; output on ;  " Solving for rho_e " ; output off;
        endif;

         {cnstlist,bbs}=setulfd2(bbs,yearid,1,infinitp,negok,dispnam1);
         foo=varput(1,"_mlatmax");    @ just 1 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", rho_e step ",
                               &figd2i,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 ;
        if prntit==1 ;
           ?; output on ;  " Solving for rho_y " ; output off;
        endif;
         {cnstlist,bbs}=setulfd2(bbs,yearid,2,infinitp,negok,dispnam1);
         foo=varput(1,"_mlatmax");    @ just 1 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", rho_y step ",
                               &figd2i,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 3 -- solving for lambda @
step3: ;
        if prntit==1 ;
          ?;  output on ;  " Solving for lambda " ; output off ;
        endif;
         {cnstlist,bbs}=setulfd2(bbs,yearid,3,infinitp,negok,dispnam1);
         foo=varput(1,"_mlatmax");    @ just 1 criteria @
         {loglik,bbs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", lambda step ",
                               &figd2i,dispnam1);
         b=bbs ;

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

        iter=varget("_mliter")+iter ;

@ Step 4 -- all 3 parameters  @
       if prntit==1 ;
          ?;    output on ;  " Solving for all_3 parameters  " ; output off ;
        endif;
        {cnstlist,bbs}=setulfd2(bbs,yearid,0,infinitp,negok,dispnam1);
        foo=varput(2,"_mlatmax");    @ 2 criteria @
        {loglik,rs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", all_3 parameters ",
                               &figd2i,dispnam1);
         b=bbs ;

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

        iter=varget("_mliter")+iter ;


@ 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-Iterative concentrated (simultaneous) @
    else ;
           {cnstlist,bbs}=setulfd2(bbs,yearid,0,infinitp,negok,dispnam1);
           foo=varput(2,"_mlatmax");    @ 2 criteria @
           {loglik,rs,freeindx}=maxlikc(bbs,cnstlist,ditype,steptype,
                                        prntit,
                               header $+ ", all 3 parameters ",
                               &figd2i,dispnam1);

          iter=varget("_mliter")+iter ;
     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 @

      kb=rows(bconc);
      var=bconc[kb] ; beta=bconc[1:kb-1] ;
      b=beta|rs|var ;       @ starting values @

@ Use non-concentrated Hessian for VC matrix@
     " Computing covariance matrix ... ";
      {vc,hmtx}=compcvc(b,vctype,miss(1,1),&figd2);
      stderr=sqrt(diag(vc));
      tstat=b./stderr ;
      nobs2=nobs-numpanls-kb ;
      parnames=parname1|dispname|"åý";
      output on ;
      nlstprnt(_titles,dataset,depvar,nobs,nobs2,loglik,
              " Concentrated Loglik 1d_Serial ",
               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(oldmaxi,"_maxiter");   @ reset @
     call varput(oldepsb,"_mlepsb");
     call varput(oldepsg,"_mlepsg");

     retp(1);                                   @ got an answer @


notgot: ;
     call varput(oldmaxi,"_maxiter");   @ reset @
     call varput(oldepsb,"_mlepsb");
     call varput(oldepsg,"_mlepsg");

     retp(0);                                   @ no answer @

    endp;

/***************************/
@ do MLE of serial correlated irregular timed panel  -- all parameters free @
@ GRID SEARCH METHOD ONLY !!! @

  proc modeld2(f1,parname1);

  local cv,logl,status,tmpe,tmpl,bbs,g1s,li;
  local l0,u0,f0,prnt,header,kb,beta,nobs2 ;
  local hmtx,iter;


    l0=miss(1,1); u0=l0; f0=u0;

    header=" 1st diff, Rho = P_e + P_y^Dt, model";

    dispname=parname1|"RHO_E"|"RHO_Y"|"LAMBDA"|"åý" ;
    smalldat=0;                       @force read of file @


@ Use this if you want grid search @
    if dogrid==1 ;

         ?; ?;" Loading data .(grid search method).... ";?;
         {smalldat,numpanls,nobs,xid,zid}=setiod2(f1,parname1);
         if numpanls==0; retp(0); endif ;
         if scalmiss(xid); retp(0); endif ;
         {loglik,beta,rho_e,rho_y,lambda,var}=
               doiter(b,varhat,header,_useeps,_prntall);
         if scalmiss(loglik);
            " Regression aborted ";
            retp(1);
         endif ;
         b=beta|rho_e|rho_y|lambda|var ;
        " Computing covariance matrix ... ";
         {vc,hmtx}=compcvc(b,vctype,miss(1,1),&figd2);
         stderr=sqrt(diag(vc));
         tstat=b./stderr ;
         nobs2=nobs-numpanls-rows(b);
         iter=varget("_mliter") ;
         output on ;
         parnames=dispname;
         nlstprnt(_titles,dataset,depvar,nobs,nobs2,loglik,
               " 1d_SERIAL Grid search Loglik",
               dispname,b,stderr,tstat,cdftc(abs(tstat),nobs2),hmtx,iter);
          output off ;
          retp(1);

     endif;


/**************** NON-concentrated methods no longer supported  ****/

 errorlog "% Non-concentrated method no longer supported: MODELD2 " ;
   retp(0);       @ should NEVER get here @
  endp;


/*****************************/
@ concentrated log likelihood, 1st diff serial correlated model @
@ 1 params at at time @
@ Note use of global parameter ditype @

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

    local lambda,ry,re,b,ls,var,kb;
    local lisum,g1sum,h,tmp,dx,dl ;
 @ bs= b re ry lambda var @

    if ditype ne 3 ;
       {lisum,g1sum,h}=figd2c(bs,narg);
       retp(lisum,g1sum,h);
    endif;

@ else, use simple gradient (ditype==3)@

    kb=rows(bs);
    lambda=bs[3] ;
    ry=bs[2] ;
    re=bs[1];

   h=eye(kb);  lisum=0; g1sum=0;

    if narg ne 2 ;
      {lisum,var,b}=compbd2(re,ry,lambda,1);  @ 1==compute loglik @

    endif ;

    if narg==0 or narg==1 or narg==4 or narg==5 ;

       if thisfree==1 ;
          dx=maxc(abs(re*10e-5)|0.00001) ;
          re=re+dx;
       elseif thisfree==2 ;
          dx=maxc(abs(ry*10e-5)|0.00001) ;
          ry=ry+dx ;
       else ;
          dx=maxc(abs(lambda*10e-5)|0.00001) ;
          lambda=lambda+dx ;
       endif;
       {tmp,var,b}=compbd2(re,ry,lambda,1);  @ 1==compute loglik @
       dl=(tmp-lisum) ;
       if thisfree==1 ;
          g1sum=(dl/dx)|0|0 ;
       elseif thisfree==2 ;
          g1sum=0|(dl/dx)|0 ;
       else ;
          g1sum=0|0|(dl/dx) ;
       endif;

    endif;

    retp(lisum,g1sum,h);
  endp;



/*****************************/
@ concentrated log likelihood, 1st diff serial correlated model @

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

    local lambda,ry,re,b,ls,var,kb;
    local lisum,g1sum,h,tmp ;
 @ bs= b re ry lambda var @


    kb=rows(bs);
    lambda=bs[3] ;
    ry=bs[2] ;
    re=bs[1];

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

    if narg==3  ;
      {lisum,var,b}=compbd2(re,ry,lambda,1);  @ 1==compute loglik @
    elseif (narg==1) or (narg==2) ;
      " Hessian unavailable ." ;
      stop ;

    else ;                     @ 0 ,4 or 5 @

       {lisum,g1sum,h}=figd2cb(re,ry,lambda,narg);

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



/************************************/
@ compute log like, gradient, and g'g  for concentrated first diff model,
  given parameters rho_e,rho_y,lambda.
  ALSO -- set global parameters Bconc = beta|var
  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 var and beta,
  and the second pass to compute gradients wrt rho_e,rho_y and lambda.

@

    proc (3)=figd2cb(rho_e,rho_y,lambda,narg);

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

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

    dldlam=0;dldre=0;dldry=0;

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


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

    {foo,var,beta}=compbd2(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. @

      {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
           grab_d2(oldsize,nextobs,tranmtx,xid);

       if scalmiss(ynew); continue ; endif ;

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

@ Note: we don't need deriv wrt var @
      {foo,dcvdlam,dcvdre,dcvdry,cv}=
         gradcvd2(rho_e,rho_y,yrvals,lambda,var);   @ deriv of cv @

      cv=cv*var ;
      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) )   ;

 @ bs= b re ry lambda var @
       g1=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-numpanls)*ln(1/sqrt(2*pi)) );
    endif;

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

    retp(lik,gtot',ggtot);

  endp;



/************************************/
@ Non-concentrated . Used ONLY for hessian for VC matrix
 (uses figlikd2 and figd2b )
@

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


/*** NOTE: used only for VC matrix ***/


    local lambda,ry,re,b,ls,var,kb;
    local lisum,g1sum,h ;
    kb=rows(bs);                  @ bs= b re ry lambda var @
    var=bs[kb];
    lambda=bs[kb-1] ;
    ry=bs[kb-2] ;
    re=bs[kb-3];
    b=bs[1:kb-4];

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

    if narg==3  ;           @ return -loglik @
       lisum=figlikd2(f1,b,var,re,ry,lambda);

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

    else ;                     @ 0 ,4 or 5 @
       {lisum,g1sum,h}=figd2b(b,var,re,ry,lambda,narg);
    endif;
     retp(lisum,g1sum,h);

  endp;




/************************************/
@ compute log like, gradient, and g'g  for first diff model,
  given parameters b,rho_e,rho_y,lambda
  Returns li,g,g'g -- dependent on narg:
    narg=0  compute all 3
         4  G'G only
         5  Li and G
NOTE: USED BY FIGD2 -- ONLY FOR HESSIAN FOR VC MATRIX
@
    proc (3)=figd2b(b,var,rho_e,rho_y,lambda,narg);

    local inpanel,cv,lik;
    local foo,atline,atcol,nextobs,dta,oldsize,bbs;
    local tranmtx,xnew,ynew,z,cvdet,cvinv,yrvals;
    local dldvar,dldre,dldry,dldlam,dldb,dcvdvar,dcvdlam,dcvdre,dcvdry;
    local g1,gtot,ggtot ;

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

    dldvar=0;dldb=0;dldlam=0;dldre=0;dldry=0;

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

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

@  Grab a panel. @

      {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
             grab_d2(oldsize,nextobs,tranmtx,xid);

       if scalmiss(ynew); continue ; endif ;

@ else, okay panel @

      {dcvdvar,dcvdlam,dcvdre,dcvdry,cv}=
         gradcvd2(rho_e,rho_y,yrvals,lambda,var);   @ deriv of cv @

      cv=cv*var ;
      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;

   dldvar=dldvar+ (sumc(diag(cvinv*dcvdvar)) - (z'(cvinv*dcvdvar*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 re ry lambda var @
       g1=dldb~dldre~dldry~dldlam~dldvar ; 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-numpanls)*ln(1/sqrt(2*pi)) );
    endif;
    retp(lik,gtot',ggtot);

  endp;


/************************************/

@ compute  true log like for first diff model,
  given parameters b,rho_e,rho_y,lambda@

    proc figlikd2(f1,b,var,rho_e,rho_y,lambda);

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

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


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


@  Grab a panel. @

      {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
                  grab_d2(oldsize,nextobs,tranmtx,xid);
       if scalmiss(ynew);
          continue ;
       endif ;
@ else, got a valid  panel @

         cv=mkcvd2(rho_e,rho_y,yrvals,lambda)*var;

       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-numpanls)*ln(1/sqrt(2*pi));
    retp(lik+aconst);

  endp;


/***********************************/
@ there are 4 elements to take derive of the CV matrix with respect to
   sigma, rho_e, rho_y, and lambda.
retp {dcvdvar,dcvdlam,dcvdre,dcvdry,cv};
@

 proc (5)=gradcvd2(rho_e,rho_y,years,lambda,var);

local sth,tth,aprod,aprod2,tmp,tmp0,tmp1,tmp2;
local cv,dcvdvar,dcvdre,dcvdry,dcvdlam ;
local tr,dt,p,psq,pm1,pm1sq,tmp4;
local dpdre,dpdry,dpm1sqre,dpm1sqry,daprodre,daprodry ;
local prodt,dprodret,dprodryt,prods,dprodres,dprodrys ;
local j,k,m,n;

    cv=mkcvd2(rho_e,rho_y,years,lambda);
    Tr=rows(years);

    dcvdre=zeros(tr-1,tr-1); dcvdry=dcvdre; dcvdlam=dcvdry ;  @initialiaze@

    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 @
    psq=p.*p ;
    pm1=p-1;
    pm1sq=pm1.*pm1 ;

    dpdry=dt.*(rho_y.^(dt-1));
    dpdry=1|dpdry;   @ needed below @
    dpdre=ones(rows(dpdry),1);

    dpm1sqre=2*pm1 ;
    dpm1sqry=2*(pm1.*dpdry) ;


@ First, set dcvdvar .... @
  dcvdvar=cv ;                   @ dCV/dVAR is simple ! @

@ constuct dcvdre,dcvdry, and dcvdlam element by element . @

@ Create diagonal elements.  Note that element 1,1 is composed from
  t=2 data @

@ tth=2 @
  tth=2 ;
  dcvdlam[tth-1,tth-1]=pm1sq[tth] ;
  dcvdre[tth-1,tth-1]=dpm1sqre[tth]*lambda ;
  dcvdry[tth-1,tth-1]=dpm1sqry[tth]*lambda ;


@ tth=3..T @
  if tr<3; goto done ; endif ;
  tth=3 ;
  do until tth>tr ;
     tmp1=0;tmp0=0;
  {aprod,daprodre,daprodry}=dprodm(psq,2*(p.*dpdre),2*(p.*dpdry),0,tth-3,tth-1);

     dcvdlam[tth-1,tth-1]=pm1sq[tth]*aprod ;

     tmp0=tmp0+(dpm1sqre[tth]*aprod*lambda) ;
     tmp0=tmp0+(pm1sq[tth]*daprodre*lambda) ;
     tmp1=tmp1+(dpm1sqry[tth]*aprod*lambda) ;
     tmp1=tmp1+(pm1sq[tth]*daprodry*lambda) ;

     n=0;               @ when tth=3, this is skipped @
     do until n>tth-4 ;
      {aprod,daprodre,daprodry}=dprodm(psq,2*(p.*dpdre),2*(p.*dpdry),0,n,tth-1);
       tmp0=tmp0+(dpm1sqre[tth]*aprod) ;
       tmp0=tmp0+(pm1sq[tth]*daprodre);
       tmp1=tmp1+(dpm1sqry[tth]*aprod) ;
       tmp1=tmp1+(pm1sq[tth]*daprodry);
       n=n+1 ;
    endo;

    dcvdre[tth-1,tth-1]=tmp0+dpm1sqre[tth] ;
    dcvdry[tth-1,tth-1]=tmp1+dpm1sqry[tth] ;

    tth=tth+1;
  endo;

s1: ;
@ sth=2 and tth=3 @
   dcvdlam[1,2]=pm1[3]*pm1[2]*p[2] ;
   dcvdlam[2,1]=dcvdlam[1,2] ;

   dcvdre[1,2]=pm1[2]*p[2]*lambda ;
   dcvdry[1,2]=pm1[2]*p[2]*lambda*dpdry[3] ;

   dcvdre[1,2]=dcvdre[1,2]+ (pm1[3]*p[2]*lambda);
   dcvdry[1,2]=dcvdry[1,2]+ (pm1[3]*p[2]*lambda*dpdry[2]) ;

   dcvdre[1,2]=dcvdre[1,2]+ (pm1[3]*pm1[2]*lambda);
   dcvdry[1,2]=dcvdry[1,2]+ (pm1[3]*pm1[2]*lambda*dpdry[2]);

   dcvdre[1,2]=dcvdre[1,2]+1;
   dcvdry[1,2]=dcvdry[1,2]+dpdry[3];

   dcvdre[2,1]=dcvdre[1,2] ;
   dcvdry[2,1]=dcvdry[1,2] ;


@ sth=2, t=4...T @
   if tr==3 ; goto done ; endif ;
   tth=4 ;

   do until tth>tr ;

      {aprod,daprodre,daprodry}=dprod(p,dpdre,dpdry,2,tth-1);
      dcvdlam[1,tth-1]=pm1[tth]*pm1[2]*aprod  ;
      dcvdlam[tth-1,1]=dcvdlam[1,tth-1]  ;

      tmp=pm1[2]*aprod ; ;
      tmp1=tmp ;
      tmp2=tmp*dpdry[tth] ;

      tmp=pm1[tth]*aprod ;
      tmp1=tmp1+ tmp ;
      tmp2=tmp2+ (tmp*dpdry[2]) ;

      tmp1=tmp1+(pm1[tth]*pm1[2]*daprodre) ;
      tmp2=tmp2+(pm1[tth]*pm1[2]*daprodry) ;


      tmp1=tmp1*lambda ;
      tmp2=tmp2*lambda ;

      {aprod,daprodre,daprodry}=dprod(p,dpdre,dpdry,3,tth-1);
      tmp1=tmp1+(pm1[tth]*daprodre)+(aprod) ;
      tmp2=tmp2+(pm1[tth]*daprodry)+(aprod*dpdry[tth]) ;

      dcvdre[1,tth-1]=tmp1;
      dcvdre[tth-1,1]=tmp1;
      dcvdry[1,tth-1]=tmp2;
      dcvdry[tth-1,1]=tmp2;

      tth=tth+1;
   endo ;



@ s=3, t=4 @
    if tr<4 ; goto s2; endif ;

    dcvdlam[2,3]=pm1[4]*pm1[3]*p[3]*psq[2] ;
    dcvdlam[3,2]=dcvdlam[2,3];

    tth=4;
    aprod=(psq[2]*lambda)+1;
    tmp1=pm1[3]*p[3]*aprod ;
    tmp2=tmp1*dpdry[4] ;
    tmp=pm1[4]*p[3]*aprod ;
    tmp1=tmp1+tmp ;
    tmp2=tmp2+(tmp*dpdry[3]);
    tmp=pm1[4]*pm1[3]*aprod;
    tmp1=tmp1+tmp ;
    tmp2=tmp2+(tmp*dpdry[3]);
    tmp=pm1[4]*pm1[3]*p[3]*(p[2]*2*lambda);
    tmp1=tmp1+tmp ;
    tmp2=tmp2+(tmp*dpdry[2]);

    tmp1=tmp1+1;
    tmp2=tmp2+dpdry[4];

   dcvdre[2,tth-1]=tmp1;
   dcvdre[tth-1,2]=tmp1;
   dcvdry[2,tth-1]=tmp2;
   dcvdry[tth-1,2]=tmp2;

s2: ;
@ s=t-1, t=5..T @

     tth=5 ;
     do until tth>tr ;
         sth=tth-1;
         {aprod,daprodre,daprodry}=
              dprodm(psq,2*(p.*dpdre),2*(p.*dpdry),0,tth-4,tth-2);
         dcvdlam[tth-1,tth-2]=pm1[tth]*pm1[sth]*p[sth]*aprod;
         dcvdlam[tth-2,tth-1]=dcvdlam[tth-1,tth-2] ;

        tmp=pm1[sth]*p[sth]*aprod ;
        tmp1=tmp ;
        tmp2=tmp*dpdry[tth];
        tmp=pm1[tth]*p[sth]*aprod ;
        tmp1=tmp1+tmp;
        tmp2=tmp2+(tmp*dpdry[sth]);
        tmp=pm1[tth]*pm1[sth]*aprod ;
        tmp1=tmp1+tmp;
        tmp2=tmp2+(tmp*dpdry[sth]);

        tmp1=tmp1+(pm1[tth]*pm1[sth]*p[sth]*daprodre);
        tmp2=tmp2+(pm1[tth]*pm1[sth]*p[sth]*daprodry);

        tmp1=tmp1*lambda ;
        tmp2=tmp2*lambda ;

        n=0;
        do until n>tth-5 ;
            {aprod,daprodre,daprodry}=
                dprodm(psq,2*(p.*dpdre),2*(p.*dpdry),0,n,tth-2);

            tmp=pm1[sth]*p[sth]*aprod ;
            tmp1=tmp1+tmp ;
            tmp2=tmp2+(tmp*dpdry[tth]);
            tmp=pm1[tth]*pm1[sth]*aprod ;
            tmp1=tmp1+tmp;
            tmp2=tmp2+(tmp*dpdry[sth]);
            tmp=pm1[tth]*p[sth]*aprod ;
            tmp1=tmp1+tmp;
            tmp2=tmp2+(tmp*dpdry[sth]);

            tmp1=tmp1+(pm1[tth]*pm1[sth]*p[sth]*daprodre);
            tmp2=tmp2+(pm1[tth]*pm1[sth]*p[sth]*daprodry);

            n=n+1;
        endo ;

        tmp=pm1[sth]*p[sth] ;
        tmp1=tmp1+tmp ;
        tmp2=tmp2+(tmp*dpdry[tth]);
        tmp=pm1[tth]*p[sth] ;
        tmp1=tmp1+tmp;
        tmp2=tmp2+(tmp*dpdry[sth]);
        tmp=pm1[tth]*pm1[sth] ;
        tmp1=tmp1+tmp;
        tmp2=tmp2+(tmp*dpdry[sth]);

        tmp1=tmp1+1 ;
        tmp2=tmp2+dpdry[tth] ;

       dcvdre[tth-2,tth-1]=tmp1;
       dcvdre[tth-1,tth-2]=tmp1;
       dcvdry[tth-2,tth-1]=tmp2;
       dcvdry[tth-1,tth-2]=tmp2;

       tth=tth+1;
     endo ;                             @ tth @

@ s=3..t-2, t=5..T @
     tth=5 ;
     do until tth>tr ;
       sth=3;
       do until sth>tth-2 ;

@ first term @
          {prodt,dprodret,dprodryt}=
               dprodm(p,dpdre,dpdry,0,tth-3,tth-1);
          {prods,dprodres,dprodrys}=
               dprodm(p,dpdre,dpdry,0,sth-3,sth-1);

           dcvdlam[sth-1,tth-1]=pm1[tth]*pm1[sth]*prodt*prods ;
           dcvdlam[tth-1,sth-1]=dcvdlam[sth-1,tth-1] ;

           tmp1=pm1[tth]*prodt*prods*lambda ;
           tmp2=tmp1*dpdry[sth] ;

           tmp=pm1[sth]*prodt*prods*lambda ;
           tmp1=tmp1+tmp ;
           tmp2=tmp2+(tmp*dpdry[tth]) ;

           tmp0=pm1[tth]*pm1[sth]*lambda;
           tmp1=tmp1+(tmp0*dprodret*prods);
           tmp1=tmp1+(tmp0*prodt*dprodres);
           tmp2=tmp2+(tmp0*dprodryt*prods);
           tmp2=tmp2+(tmp0*prodt*dprodrys);

@ second term @
          {prodt,dprodret,dprodryt}=
               dprodm(p,dpdre,dpdry,0,tth-sth-1,tth-1);

           tmp=pm1[sth]*prodt ;
           tmp1=tmp1+tmp ;
           tmp2=tmp2+(tmp*dpdry[tth]);
           tmp=pm1[tth]*prodt ;
           tmp1=tmp1+tmp ;
           tmp2=tmp2+(tmp*dpdry[sth]);

           tmp0=pm1[tth]*pm1[sth];
           tmp1=tmp1+(tmp0*dprodret);
           tmp2=tmp2+(tmp0*dprodryt);

@ fourth term @

         {aprod,daprodre,daprodry}=
              dprodm(p,dpdre,dpdry,0,tth-sth-2,tth-1);

           tmp1=tmp1+aprod ;
           tmp2=tmp2+(aprod*dpdry[tth]);
           tmp1=tmp1+(pm1[tth]*daprodre);
           tmp2=tmp2+(pm1[tth]*daprodry);


@ Third term. if s=3, skip this term @
         if sth==3; goto nos3 ; endif ;

         {aprod,daprodre,daprodry}=
              dprodm(p,dpdre,dpdry,0,tth-sth-1,tth-1);
          n=tth-sth ;
          do until n>tth-4 ;
            {prodt,dprodret,dprodryt}=
                dprodm(psq,2*(p.*dpdre),2*(p.*dpdry),tth-sth,n,tth-1);

             tmp=pm1[sth]*aprod*prodt ;
             tmp1=tmp1+tmp;
             tmp2=tmp2+(tmp*dpdry[tth]);
             tmp=pm1[tth]*aprod*prodt ;
             tmp1=tmp1+tmp;
             tmp2=tmp2+(tmp*dpdry[sth]);

             tmp0=pm1[sth]*pm1[tth] ;

             tmp1=tmp1+(tmp0*daprodre*prodt);
             tmp2=tmp2+(tmp0*daprodry*prodt);
             tmp1=tmp1+(tmp0*aprod*dprodret);
             tmp2=tmp2+(tmp0*aprod*dprodryt);


             n=n+1;
         endo;


nos3: ;         @ skip here if sth=3 @

         dcvdre[sth-1,tth-1]=tmp1 ;
         dcvdre[tth-1,sth-1]=tmp1 ;
         dcvdry[sth-1,tth-1]=tmp2 ;
         dcvdry[tth-1,sth-1]=tmp2 ;

         sth=sth+1;
     endo;
     tth=tth+1;
   endo;

done:;
    dcvdre=dcvdre*var ;
    dcvdry=dcvdry*var ;
    dcvdlam=dcvdlam*var ;

    retp(dcvdvar,dcvdlam,dcvdre,dcvdry,cv);
endp;

/******************************************/
@ create the CV matrix for the first-differenced serial-correlated
  irregular-timing model.  User supplies:
  rho_e = event correlation coefficient
  rho_y = year correlation coefficient)
  years     = Tx1 vector of 'years' of event occurences
  lambda = the past history correction factor
@

 proc mkcvd2(re,ry,years,lambda);


@  there are a number of special cases.  NOte that we derive the
   lower diagonal, and then copy this to the upper diagonal

Note: dimension of CV matrix will be (Tr-1)x(Tr-1), with Tr=rows(Years)
@
    local
      tr, dt, p, psq, pm1, pm1sq, cvmtx, tmp, it, tmp1, tmp2, n, is, tmp4,
      tmp3, tmp3a ;

     Tr=rows(years);

@ vector of P correlation coefficients, and related values @

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

@ useful stuff used below @
     psq=p.*p ;
     pm1=p-1;
     pm1sq=pm1.*pm1 ;

     cvmtx=zeros(tr-1,tr-1);   @ vc matrix is T-1xT-1 @

    tmp=(lambda*pm1sq[2]) + 1 ;                         @ 1,1 element @

@ 2 years, CV is 1x1 @
     if tr==2 ;
        retp(tmp);
     else ;
         cvmtx[1,1]=tmp ;
     endif ;

@ s=t=3  @
    tmp=(psq[2]*lambda)+1 ; tmp=tmp*pm1sq[3];
    cvmtx[2,2]=tmp+1;


@ s=2,t=3 @
    tmp=(pm1[2]*p[2]*lambda)+1 ;
    cvmtx[1,2]=tmp*pm1[3]; cvmtx[2,1]=cvmtx[1,2];

    if tr==3 ; retp(cvmtx); endif ;


@ > 3 years @

@ s=t=4..T @
    it=4 ;
    do until it gt tr ;
       tmp1=pm1sq[it] * prodc(psq[2:(it-1)]) * lambda  ;
       tmp2=0;
       n=0;
       do until n gt (it-4) ;
         tmp2=tmp2+prodc(psq[(it-1-n):(it-1)]);
         n=n+1;
       endo;
       tmp2=tmp2*pm1sq[it];
       cvmtx[it-1,it-1]=tmp1+tmp2+pm1sq[it]+1 ;
       it=it+1;
     endo ;

@ s=2, t=4..T @
    it=4 ;
    do until it gt tr ;
      tmp=pm1[2]*prodc(p[2:(it-1)])*lambda ;
      tmp=tmp+  prodc(p[3:(it-1)]);
      cvmtx[it-1,1]=pm1[it]*tmp ; cvmtx[1,it-1]=cvmtx[it-1,1];
      it=it+1 ;
    endo ;

@ s=3,t=4 @
    tmp=pm1[3]*pm1[4]*p[3] * ( (psq[2]*lambda)  + 1 ) ;
    cvmtx[2,3]=tmp+pm1[4] ; cvmtx[3,2]=cvmtx[2,3] ;


    if tr==4  ; retp(cvmtx); endif ;

@ > 4 years @

@ do tridiagonal elements @

   it=5;
   do until it gt tr ;
       is=it-1;
       tmp1=prodc(psq[2:(it-2)]) * lambda;
       n=0;  tmp2=0;
       do until n gt (it-5) ;
          tmp2= tmp2 + prodc(psq[(it-n-2):(it-2)]);
          n=n+1;
       endo;
       tmp=(tmp1+tmp2+1)*pm1[is]*p[is] ;
       cvmtx[it-1,is-1]=(tmp+1)*pm1[it] ;
       cvmtx[is-1,it-1]=cvmtx[it-1,is-1];
       it=it+1;


   endo ;


@ do all other elements @
   it=5 ;
   do until it gt tr ;
     is=3 ;
     do until is gt it-2 ;
        tmp1=pm1[is]* prodc(p[2:(it-1)]) * prodc(p[2:(is-1)]) * lambda ;
        tmp2=pm1[is]*prodc(p[is:(it-1)]) ;
        tmp4=prodc(p[(is+1):(it-1)]);

        tmp3=0;
        if is>3 ;                       @ eps_3* eps_t is special @
          tmp3=pm1[is]*prodc(p[is:(it-1)]);
          n=it-is ; tmp3a=0;
          do until n gt it-4 ;
             tmp3a=tmp3a+prodc(psq[(it-n-1):(it-(it-is)-1)]) ;
             n=n+1;
          endo;
          tmp3=tmp3a*tmp3 ;
        endif;
        cvmtx[is-1,it-1]=(tmp1+tmp2+tmp3+tmp4)*pm1[it] ;
       cvmtx[it-1,is-1]=cvmtx[is-1,it-1];
        is=is+1;
     endo;
     it=it+1;
   endo;


  retp(cvmtx);
endp;



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

   proc mkcvd2a(rho,years);

   local inpanel,a,b,c,pmtx,kol,ro,nd,nd1;
   local rhond,rhond1;

   inpanel=rows(years);

   PMTX=eye(inpanel-1);

   A=(rho^2-rho^3 ) / (1+rho) ;
   b= (rho*rho)-rho;
   c= rho-1 ;

   pmtx=pmtx*(a+(c*c)+1) ;

   ro =0;
   do until ro == inpanel-1;
     ro =ro +1;
     KOL=0;
     do until Kol== ro;
        Kol=Kol+1;
        nd=ABS(ro -Kol) ;
        if nd ne 0;
          ND1=ND-1; RHOND=RHO^ND; RHOND1=RHO^ND1;
          pmtx[ro,Kol]=(rhond*A) + (rhoNd1*B*C) + (rhond1*C) ;
          pmtx[kol,ro]=pmtx[ro,kol] ;
        endif;          @ else, leave as is @
     ENDO;
    ENDO;

    retp(pmtx);

   endp;




/*************************/
@set bounds for first diff model @
@   vars  = rho_e|rho_y|lambda  --- concentrated @

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

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

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


  if usec==0 ;          @ concentrated @

      let cnstlist= 0 0 -1 1 -1e-9 ;   @ lambda constraint @
      if infinitp==1 ; cnstlist[5]=-1 ; endif ;  @ lambda gt 1 @
      if scalmiss(yearid) ;             @ no rhoy @
          let tmp1= 0 1 0 2  0 ;   @ rhoy=0 @
          let tmp2= 1 0 0 1  0.9999 ;   @ rhoe<1 @
          let tmp3= -1 0 0 1 0 ; @ rhoe>0@
          cnstlist=(cnstlist'|tmp1') ;
      else ;                          @ rhoy and rhoe @
          let tmp1= 1 1 0 1 0.999 ; @ re+ry<1 @
          let tmp2= -1 -1 0 1 0 ;    @ re+ry>0 @
          let tmp3= 1 0 0 1 0.999 ; @ rhoe < 1 @
          let tmp4= 0 1 0 1 0.999 ;   @ rhoy < 1 @
          cnstlist=(cnstlist'|tmp1'|tmp2'|tmp3'|tmp4'|tmp5'|tmp6') ;
      endif;

  else ;
       cnstlist=eye(3);
       cnstlist=cnstlist~(2*ones(3,1));
       cnstlist=cnstlist~beta;

      if usec==1 ;                @ only rho_e param free @
             let tmp3=1 0 0 1 0 ; tmp3[5]=0.9999-beta[2]; @ rhoe+rhoy<1@
             cnstlist[1,.]=tmp3' ;
             tmp=0.9999-beta[2] ;
             beta[1]=minc(beta[1]|tmp);
             cnstlist=cnstlist|tmp5' ;

      elseif usec==2 ;           @ only rho_y free @
             let tmp3= 0 1 0 1 0 ; tmp3[5]=0.9999-beta[1]; @ rhoe+rhoy<1@
             cnstlist[2,.]=tmp3' ;
             tmp=0.99999-beta[1]  ;
             beta[2]=minc(beta[2]|tmp);
             cnstlist=cnstlist|tmp6' ;

       elseif usec==3 ;           @ only lambda free @
          let tmp3=0 0 1 1 1e26 ;   @ lambda lt infinity @
          if infinitp==1 ;        @ process begins at time -infinity @
             tmp=beta[1]+beta[2] ;       @which implies bounds on lambda@
             tmp= ( 1/(1-tmp^2)) ;
             tmp3[5]=tmp ;                 @ upper bound not infinity @
             cnstlist[3,.]=tmp3' ;

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

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

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

          else ;
             cnstlist[3,.]=tmp3' ;
             tmp3[3]=-1 ; tmp3[5]=-1e-9 ;
             cnstlist=cnstlist|tmp3' ;
          endif;

      endif;

    endif ;

    thisfree=usec ;      @ GLOBAL used by figd2i @

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



      retp(cnstlist,beta);



  endp;



/************************************/
@ user grid search  to find best likelihood.. d2 model

  There are K+4 parameters to estimate.  K+1 of them (betas and var)
   are strictly
   dependent on the other 3, these 3 being rho_e, rho_y, and lambda
   -- which are used in a GLS estimator for beta, var.

   A grid search is performed over these 4 elements. First, rho_e is set
   (allowed to vary between 0 and 1.0 ).  Then, rho_y is set (allowed to
   vary between 0 and 1-rho_e).  Then lambda is set (allowed to vary
   between 0 (or 1-(1-(rho_E+rho_y)^2) and infinity.   Lastly, var is set.
   Then, the log-likelihood is found.

        L = SUM,i { ln(det(CV[i]) + z[i]'inv(CV[i])*z[i]  }

   Noting that CV[i] = CV1[i]*var, we can simplify L to

        L = SUM,i {ln(det(CV1[i]))}
              +  N*ln(var)
              +  Sum,i {Z[i]'inv(CV1[i])*Z[i]} / var

    so    ==> var=J/N
    where J equals Sum,i {Z[i]'inv(CV1[i])*Z[i]}

   A grid search is  over rho_e, rho_y.  For lambda, a form of
   golden search is used.
-----------------------------------------------------@

    proc (6)=doiter(b,var,header,eps,doall);

    local ry,re,aninc,tmp,lm0,li,varuse,beta,arow,acol;
    local buse,liuse,likmax,ryuse,reuse,lmuse,lmgot,vargot ;
    local re0,re1,ry0,ry1,oldbest,ryloop,iter,lm0a ;
    local _mliter,liworse,foo ;

    beta=miss(1,1);
    liuse=miss(1,1); buse=0; reuse=0;ryuse=0;lmuse=0;varuse=0;

@ Use _useinc as initial step size @
      aninc=_useinc ;

      ?;?; ?; arow=csrlin-1 ; acol=csrcol ;

      re=_usere0 ;     likmax=-1e90 ; oldbest=0; iter=0;

      do until fge(re,minc(1.0|_usere2)) ;
        ry=_usery0 ; ryloop=0;
        liworse=miss(1,1);      @ skip ONLY in rhoy direction @

        do until ryloop==1;
           lm0a=1/(1-(re*re)) ;          @ minimum allowable lambda @
           if infinitp==1 ; lm0=lm0a ; else ; lm0=_uselmin ; endif ;

           if (abortit2(arow,acol,re,ry))==1; goto abort2 ; endif ;

           if ry eq 0 ;
               lmgot=lm0a;
               {li,vargot,beta}=compbd2(re,ry,lmgot,1);  @returns true loglik@
           else;
               if infinitp==0 ;
                  tmp=_uselm2 ;
                else ;
                  tmp=1-( (re+ry)^2); tmp=1/tmp ;
                  tmp=minc(tmp|_uselm2);           @ maximum allowable lambda@
                endif ;
                {lmgot,vargot,beta,li}=goldd2(re,ry,lm0,0.001,tmp);
           endif;
           if li>likmax ;               @ maximize likelihood @
               buse=beta; liuse=li ; reuse=re; ryuse=ry;
               varuse=vargot ; lmuse=lmgot ; likmax=li ;
               liworse=miss(1,1);
           else ;                   @ not an improvement @
               if _usequik==1 ;         @stop if worsening ? @
                  if scalmiss(liworse) ;      @ first "worse" @
                     liworse=li ;
                  else ;                @ we have a prior "worse" @
                     if li<liworse ;   @ continued decrease? @
                          ryloop=1 ;     @ so skip out @
                     endif ;
                  endif;
               endif;          @ usequik @
            endif;

            if (li>likmax) or (doall==1) ;
               if li<likmax ; foo=arow+4 ; else ; foo=arow+10 ; endif ;
               call doprint(prntit,re,ry,lmgot,vargot,li,beta,foo);
            endif;

            if scalmiss(yearid);
              ryloop=1;
            else;
               ry=ry+aninc ;
               if ry > _usery2; ryloop=1 ; endif;
               if fge(ry+re,1.0) ; ryloop=1; endif;
             endif;

        endo;                   @ ryloop @

        re=re+aninc ;
      endo;

loop1: ;
  iter=iter+1;
  if (aninc < eps )   ;                 @ close enuf ?@
       output on ;
        _mliter=iter;
        retp(liuse,beta,reuse,ryuse,lmuse,varuse);
    endif;

@ else, focus in.  Note if no year var, then ry0=0 @
   re0=maxc(_usere0|(reuse-aninc));
   re1=minc(_usere2|(reuse+aninc));
   ry0=maxc(_usery0|(ryuse-aninc));
   ry1=minc(_usery2|(ryuse+aninc));
   aninc=aninc/5 ;
   re=re0;
   liworse=miss(1,1);

   do until fge(re,re1) ;          @ re loop @
        ry=ry0 ; ryloop=0;

        do until ryloop==1;
           lm0a=1/(1-(re*re)) ;          @ minimum allowable lambda @
           if infinitp==1 ; lm0=lm0a ; else ; lm0=_uselmin ; endif ;

           if (abortit2(arow,acol,re,ry))==1; goto abort2 ; endif ;

@ get loglik, etc. @
           if ry eq 0 ;
               lmgot=lm0a;
               {li,vargot,beta}=compbd2(re,ry,lmgot,1);
           else;
               if infinitp==0 ;
                  tmp=_uselm2 ;
                else ;
                  tmp=1-( (re+ry)^2); tmp=1/tmp ;
                  tmp=minc(tmp|_uselm2);           @ maximum allowable lambda@
                endif ;
               {lmgot,vargot,beta,li}=goldd2(re,ry,lm0,0.001,tmp);
           endif;
           if li>likmax ;                       @ maximize true loglik @
               buse=beta; liuse=li ; reuse=re; ryuse=ry;
               varuse=vargot ; lmuse=lmgot ; oldbest=likmax;likmax=li ;
               liworse=miss(1,1);
           else;
               if _usequik==1 ;         @stop if worsening ? @
                  if scalmiss(liworse) ;      @ first decrease @
                     liworse=li ;
                  else ;                @ we have a prior "worse" @
                     if li<liworse ;   @ continued decrease! @
                          ryloop=1 ;     @ so skip out @
                     endif ;
                  endif;
               endif;          @ usequik @
            endif;

            if (li>likmax) or (doall==1) ;
               if li<likmax ; foo=arow+4 ; else ; foo=arow+10; endif ;
               call doprint(prntit,re,ry,lmgot,vargot,li,beta,foo);
            endif;

            if scalmiss(yearid);
               ryloop=1 ;
            else;
                ry=ry+aninc ;
                if fge(ry+re,1) or fge(ry,ry1); ryloop=1; endif ;
            endif;
        endo;

        re=re+aninc ;
   endo;

   goto loop1 ;      @ keep on trucking @


abort2: ;?;?;?;?;?;?;?;?;?; " Regression aborted.";
    _mliter=iter;
    retp(liuse,buse,reuse,ryuse,lmuse,varuse);

  endp;


/*********************************/
@ golden algorithim to minimize -loglik of serial correlated first diff,
  over lambda.

        {lambda,var,beta,val}=goldd2(re,ry,lm0,eps,lmupper);

         INPUTS
           re,ry -- set parameters
           lm0   -- min allowable lambda
           lmupper -- max allowable lambda
           eps   -- search tolerance (min step size tolerated)
         OUTPUTS
            lambda,var,beta = maximizing parameters for re,ry
           val  = log likelihood at stepsize*di

        IF val=missing, problem or abort.

@

        proc (4)=goldd2(re,ry,lm0,eps,lmupper);

local model:proc ;
local maxtries,  v1, v2, cd, cs, vr0, vr1, vr2, a0, b00, g0, g1, g2,
l, i, rd, rs, del, sig, s,t;
local ktht,diff0,diff1,foo,maxst,di;
local var,beta,amiss,var0,b0 ;
local lm1,var1,var2,b1,b2,lm2,TMP ;


  if lm0>lmupper ;
          retp(miss(1,1),0,0,0);
  endif;

/* --------------- initializations ------------------- */
   maxtries=100;    @ Will stop trying to find min after this many tries @
   v1=(3-sqrt(5))/2; v2=(sqrt(5)-1)/2;
   cd=1; cs=1;    @ if 1, compute rd, rs, respectively, in seven below @
   t=1;
   di=lm0;
   amiss=miss(1,1);
   maxst=1000000;
/* --------------- test to see if [0,t] brackets minimum ---------------- */
step1: ;
    {vr0,var0,b0}=compbd2(re,ry,lm0,1); vr0=-vr0;
    if badmtx(Vr0);        @ shouldn't happen, but .. @
        retp(0,0,0,amiss);
    endif;
step2: ;
    lm1=minc((lm0+(t*di))|lmupper);
    {vr1,var1,b1}=compbd2(re,ry,lm1,1);vr1=-vr1 ;
    if badmtx(vr1) ;           @try smaller initial step @
       t=t/20;
       if t< eps ;
          retp(lm0,var0,b0,vr0);
       else ;
          goto step2;
       endif;
    endif;


g0=0; g1=t; ktht=1; lm2=lm1; var2=var1;b2=b1;vr2=vr1 ;
if vr0 <= vr1; a0=0; b00=t; diff0=abs(vr1-vr0); goto seven; endif;


/* -- if not, expand interval by increments of t until min is bracketed -- */
four:  g2=g1+(t*(1.5^ktht));     @bigger expansions at each iteration @
    ktht=ktht+1;

    if g2>maxst  ;                      @ maxed out at maxstep @
         retp(lm2,var2,b2,vr2);
    endif;

    IF (LM0+(G2*DI)) GE LMUPPER ;
        G2=(LMUPPER-LM0)/DI ;
    ENDIF;
    lm2=lm0+(g2*di);
    {vr2,var2,b2}=compbd2(re,ry,lm2,1); vr2=-vr2 ;

    if badmtx(vr2);                   @ cant bracket minimum, so return vr1@
         retp(lm1,var1,b1,vr1);
    endif ;

    if vr2 >= vr1;
        a0=g0; b00=g2; diff0=abs(vr2-vr0);
        goto seven;
    else;
         if FGE(lm2,lmupper) ;             @ maxed out at upper bound @
              retp(lm2,var2,b2,vr2);
         endif;
        vr0=vr1 ; vr1=vr2; lm1=lm2 ; b1=b2; var1=var2 ;
        g0=g1; g1=g2; goto four;
    endif;

/* ---------------- search to isolate minimum ---------------------------- */
seven: i=1;
if eps<1e-05; eps=1e-05; endif;         @ just in case of wierdness@

do until i > maxtries;
    l=b00-a0;
    if abs(l) < eps; goto fin; endif;      @step size small? enuf @
    if cd; del=a0+v1*l;
       TMP=LM0+(DI*DEL);
       {rd,var2,b2}=compbd2(re,ry,TMP,1); rd=-rd ;
    endif;
    if cs; sig=a0+v2*l;
       TMP=LM0+(DI*SIG);
       {rs,var2,b2}=compbd2(re,ry,lm0+(di*sig),1); rs=-rs ;
    endif;
    if badmtx(rd|rs);            @if bad values, assume aborted @
        b00=a0;
        goto aborted ;
    endif;

    diff1=abs(rd-rs) ;                          @dmh addition @
    if abs(diff1) < eps ; goto fin ; endif ;
    if diff1/abs(rd) < eps ;goto fin ; endif ;
    if rd < rs; b00=sig; rs=rd; sig=del; cs=0; cd=1;
          else; a0=del; rd=rs; del=sig; cs=1; cd=0;
    endif;
i=i+1;
endo;

fin:
s=(a0+b00)/2;
{rs,var2,b2}=compbd2(re,ry,lm0+(di*s),1);
retp(lm0+(di*s),var2,b2,rs);

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

endp;

/****************************/
@ grab needed info  (note: only extract xid vars )
   {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
         grab_d2(oldsize,nextobs,tranmtx,xid)

@

 proc (6)=grab_d2(oldsize,nxt,tranmtx,xid);

   local dta,xnew,ynew,inpanel,nxtpnl,a1,xt,yt,tmp,yrs,pat ;

    nxt=missrv(nxt,0);
    if nxt==0;
           retp(miss(1,1),0,0,0,miss(1,1),0) ;
    endif;

    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) ;             @ panelens is global @
             nxtpnl=miss(1,1);
        else;
             nxtpnl=nxtpnl+1 ;
        endif;
        inpanel=rows(yt);
        xt=xt[1:inpanel-1,xid]; yt=yt[1:inpanel-1];
        retp(yt,xt,yrs,0,nxtpnl,0);

     else ;                             @ else, read from datafile @
         {yt,xt,yrs,nxt}=grabpanl(f1,nxt);

         if rows(y)==1;                  @ if empty, skip to next panel @
           retp(miss(1,1),0,0,0,nxt,0) ;
         endif;

         inpanel=rows(y);
         if inpanel/=oldsize ;                     @dont use prior matrix @
            tranmtx=chosemtx(inpanel);
            oldsize=inpanel;
         endif;                          @else, use prior value of tranmtx@
         xnew=tranmtx*xt; ynew=tranmtx*yt;
         retp(ynew,xnew[.,xid],yrs,oldsize,nxt,tranmtx);

      endif;

endp;


/*********************************/
@ see if we can fit it all into memory.  Sets a bunch of globals
    {smalldat,numpanls,nobs,xid,zid}=setiod2(f1,parname1)
@
   proc (5)=setiod2(f1,b);

    local nk,vcmem,memavail,nr,vn,x,y;
    local ynew,xnew,yrvals,oldsize,nextobs,tranmtx,amiss,tmp,iat,ry ;
    local novary,smalldt0,xid,zidnobs,numpanls,aa ;

    novary=1;                           @ for pointing to panel invariant vars@
    amiss=miss(1,1);
    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 @

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

       if (numpanls % 25)==1 ; format /ld 5,0; locate aa,1 ;
          numpanls " read ... ";
       endif;
       {y,x,yrvals,nextobs}=grabpanl(f1,nextobs);

       if scalmiss(y);
            continue ;          @ not a valid panel @
       else;                     @ process it @
           numpanls=numpanls+1;
           nobs=nobs+rows(yrvals);
           novary=novary.* dotfeq(minc(x),maxc(x)) ;

           if smalldt0==1 ;               @ store it @
              if rows(y) ne oldsize ;
                 oldsize=rows(y);
                 tranmtx=chosemtx(oldsize);
               endif;
               xnew=tranmtx*x ; ynew=tranmtx*y ;
@ add a row to keep consistent with yrvals @
               xnew=xnew|ones(1,cols(xnew));
               tmp=xnew~(Ynew|1);
               ry=rows(yrvals)-1;
               if scalmiss(bigdata);
                   bigdata=tmp~yrvals; panelens=iat~(iat+ry);
               else;
                   bigdata=bigdata|(tmp~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. ";
    endif;


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

endp;


/**********************************/
@ GIven re,ry,lambda -- compute beta and var that minimizes log like.
 For serial correlated first diff model
 Use global var f1
    { lik,var,b}=compbd2(re,ry,lambda)
@

   proc (3)=compbd2(re,ry,lambda,dolik);


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

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

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

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

        {ynew,xnew,yrvals,oldsize,nextobs,tranmtx}=
             grab_d2(oldsize,nextobs,tranmtx,xid);
         if scalmiss(ynew);
            continue ;
         endif;

@ else, got a panel @

        if ry==0;               @ use simple form @
            cv=mkcvd2a(re,yrvals);
        else;
            cv=mkcvd2(re,ry,yrvals,lambda);
        endif;
        trap 1 ;
        cvinv=invpd(cv);
        if scalerr(cvinv);
            output on ; format 12,6 ;?;
           " Can't take inverse of :" cv ;
           " Re ry lambda " re ry lambda ;
           " Yrvals " yrvals'; ?;
           output off ;
           goto aborted ;
        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 @
    bgls=invpd(xcvx)*xcvy ;

@ and var @
    var=(ycvy- (xcvy)'bgls)/(nobs-numpanls) ;

@ and lik @
    if dolik==1 ;
       lik=figlikd2(f1,bgls,var,re,ry,lambda);
    else;
      lik=miss(1,1);
    endif;

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

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

  endp;



