@ Include file for SERPANEL.PRG @

  clear yearname,yearid ; yearid=miss(1,1);

@ year variable @
  clear yearvals;
  clear sortopt;                  @flag: sort data set on yearvals @
  sortopt="OFF";

@ value of  panelid variable @
  clear panlvalu ;

@ model selected @
  clear model;

@ DO not ask for weight var @
  popname=miss(1,1);

@ for batch processing @
  CLEAR COMMANDS ;
  outname="SERPANEL.OUT";

@ regression vars @
  clear x,y,qconst;
  clear dta;                  @ make global for debugging purposes @
  clear parname1,parnames ;   @ parnames is what was last printed by PRNTOLS @
                              @ it can be used by SAVE @

@ Z variables are panel invariant @
  clear zid ;
@ X vars are panel variand @
  clear xid ;

@ for mle gradient, and grid search @
clear xuse,yuse ;
@ for io @
clear bigdata,panelens,smalldat ;

@ Coefficients, etc computed. Make global so that SAVE will work @
   clear xx,xy,yy,ysum;                         @moment matrices @
   clear ssr,varhat,xxinv,vc,b,ybar,sst,loglk;    @ computed in  compols @
   clear numpanls,nobs,df,totobs ;                  @ tabulated in modelXX @
   clear  sigma,stderr,t,pvt;                   @ computed in prntols1@
   clear li,b,varfe,vargrp,rho_e,rho_y,lambda,var ;
   clear var,loglik,tstat ;           @  results @

  clear crlf ;                  @ for building multiline string @
  crlf=chrs(10~13);

  clear diffmtx;                        @ global 1st diff matrix @
@ Used to first difference in small panels @
let  diffmtx[6,7]=
                   -1  1  0  0  0  0  0
                    0 -1  1  0  0  0  0
                    0  0 -1  1  0  0  0
                    0  0  0 -1  1  0  0
                    0  0  0  0 -1  1  0
                    0  0  0  0  0 -1  1 ;


  clear f1;             @ other stuff @

clear freeg,th;
clear prntit,b ;

clear dogrid ;          @ use grid search @
clear infinitp ;       @ assume NOT infinite process (lambda is free ) @
clear negok ;
clear mleconc,bconc;   @ use concentrated likelihood @
clear _max1itr ; _max1itr=3 ;     @ # of iterations in 1 step algorithim @

clear dispname,parnames ;
clear thisfree ;                        @ used by figd2i   @

output file=^outname ; output off ;

clear _useinc,_useeps,_usere0,_usere2,_usery0,_usery2,_uselm2,_prntall ;
clear _uselmin,_usequik ;
clear gridparm ;
let gridparm=
 _useinc _useeps _usere0 _usere2 _usery0 _usery2 _uselmin _uselm2 _prntall
 _usequik ;

clear _usems ;
_usems=make_mtx("_use",11);  @ grid search paramter changer @
           _use01=" Initial step size " ;
           _use02=" Minimum step size  "  ;
           _use03=" Minimum allowable rho_e  " ;
           _use04=" Maximum allowable rho_e  " ;
           _use05=" Minimum allowable rho_y  " ;
           _use06=" Maximum allowable rho_y  " ;
           _use07=" Minimum allowable lambda " ;
           _use08=" Maximum allowable lambda " ;
           _use09=" 1=Print ALL results "  ;
           _use10=" 1=If likelihood increases, skip to next rho_e ";
           _use11=" Select a grid search parameter to change:";

clear steptype,ditype,vctype;
ditype=2; steptype=2 ;  vctype=2;

 let spopts=RUN VARIABLE MODEL ALG_SET STARTá YEAR LOGLIK SORT CHOSEOBS
        STATS TITLE DATAFILE PRNT OUTPUT SAVE LISTVARS DOS EXIT QUIT ;
 spopts2=figvlite(spopts)|(0~"?"~0)|("?"~0~0);

  clear prompts ;
  prompts=make_mtx("SRHELP",rows(spopts)+2) ;
  srhelp01=" Estimate coefficients. ";
  srhelp02=" Select independent, dependent, etc. variables ";
  srhelp03=" Select model (RE, FE, etc) to estimate ";
  srhelp04=" Change MLE technique and related control parameters ";
  srhelp05=" Toggle user-selection of starting values ";
  srhelp06=" Select 'YEAR' variable ";
  srhelp07=" Compute log-likelihood ";
  srhelp08=" Toggle sort of panels on event year ";
  srhelp09=" Enter observation selection criteria ";
  srhelp10=" Compute basic statistics ";
  srhelp11=" Change title written to output ";
  srhelp12=" Select GAUSS datafile ";
  srhelp13=" Toggle printing of intermediate output ";
  srhelp14=" Change output file ";
  srhelp15=" Save/change global parameters/coefficient vectors ";
  srhelp16=" List currently chosen variables ";
  srhelp17=" Temporary exit to DOS ";
  srhelp18=" Exit to GAUSS ";
  srhelp19=" Exit to DOS ";
  srhelp20=" More help ... ";
  srhelp21=" Enter SERPANEL Option: ";
  prompts=prompts~(3*ones(rows(prompts),1));
  prompts[rows(prompts),2]=2 ;

clear mnames,mname0;
    let mnames=OLS  "1ST_DIFF" RE_SER RE_SPOR   KIEFER  "1D_SPOR" ;
    let mname0=O D RE0 RE K D2 ;




/********************************************/
@ OLS @

  proc modelo(f1,parname1,doout);

     local atline,atcol,nextobs,status,foo,kd;

    xx=0;xy=0;yy=0;nobs=0;ysum=0 ;      @ computed by makemom @
    totobs=rowsf(f1);
    xid=miss(1,1); zid=miss(1,1);

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

    numpanls=0; nextobs=1; smalldat=0;
    foo=seekr(f1,1);              @ rewind @
    do until scalmiss(nextobs);    @ ----- LOOP THRU DATA -------- @

       if abortreg(numpanls,atline,atcol)==1;
           retp(0);
       endif;

@  Grab a panel. @

       {y,x,yearvals,nextobs}=grabpanl(f1,nextobs);

        if scalmiss(y);                  @ if empty, skip to next panel @
           goto endloop;
        else;
           numpanls=numpanls+1;
        endif;

@ For OLS, no transformations required. So just  as is, to add to moment,
  etc. matrices @

         {xx,xy,yy,ysum,nobs}=makemom(x,y,xx,xy,yy,ysum,nobs);

endloop: ;                      @ jump here if no valid obs in this panel @
    endo;

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

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

@ display results (the global vars ) @
     if doout==1;
       status=PRNTOLS1(parname1," OLS Regression, no transformations.");
     endif;
     retp(1);

 endp;

/*************************************************/
@ First difference regression @

   proc modelD(f1,parname1,prntit);

     local atline,atcol,nextobs,status,foo;
     local inpanel,oldsize,xnew,ynew,tranmtx;
     local novary,kd;


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

    xx=0;xy=0;yy=0;nobs=0;ysum=0 ;      @ computed by makemom @
    novary=1;                           @ for pointing to panel invariant vars@
    totobs=rowsf(f1);

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

       if abortreg(numpanls,atline,atcol)==1;
           retp(0);
       endif;

@  Grab a panel. @

       {y,x,yearvals,nextobs}=grabpanl(f1,nextobs);
        if rows(y)==1 ;         @ if empty, skip to next panel @
           continue;
        endif;

        inpanel=rows(y);
        if  inpanel==1 ;         @ need > 1 obs for 1st diff@
           continue     ;
        else;
           numpanls=numpanls+1;
        endif;

    novary=novary.* dotfeq(minc(x),maxc(x)) ; @ used to create xid,zid @

@ get transformation matrix, if new one needed@
   if inpanel/=oldsize ;                           @dont use prior matrix @
      tranmtx=chosemtx(inpanel);
      oldsize=inpanel;
   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;

@ create pointers to panel variant  & panel invariant vars @
   {xid,zid}=getxzid(novary);

   if scalmiss(xid);
        " No panel variant variables. ";
         retp(0);
    endif;

@  extract elements of xx,xy that are panel variant @
    xx=xx[xid,xid];
    xy=xy[xid];

@ Compute Df @
    kd=rows(xid);
    df=nobs-kd ;            @ note that nobs is reduced by first differencing@
    if df<1;
         " Not enough degrees of freedom .. ";
        retp(0);
    endif;

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

@ display results (the global vars ) @
    if prntit==1 ;
       status=PRNTOLS1(parname1[xid]," First Differencing Regression.");
    endif;

    retp(1);

 endp;



/***************************/
@ do Kiefer model @

 proc modelk(f1,parname1);

   local nextobs,tranmtx,oldsize,foo,novary,maxsize,xx,xy,ysum,inpanel;
   local xnew,ynew,df,kd,b,x,y ;
   local yytot,qvyar,vuse,vuseinv,qyvar ;
   local d1,d2,omat,fmt,mask,atline,atcol ;
   local pvt ;

    atline=csrlin ; atcol=csrcol ;
    nextobs=1; tranmtx=0;oldsize=0;
    foo=seekr(f1,1);                  @ rewind @
    novary=1;                         @ for pointing to panel invariant vars@

    numpanls=0;  maxsize=0;
    xx=0;xy=0;yy=0;ysum=0; nobs=0;     @ used by makemom @
    totobs=rowsf(f1); foo=0;

    do until scalmiss(nextobs) ;

       if abortreg(foo,atline,atcol)==1;
           retp(0);
       endif;

       {y,x,foo,nextobs}=grabpanl(f1,nextobs);
        if ismiss(y)==1 ;                  @ if empty, skip to next panel @
           continue;
        endif;
        inpanel=rows(y);
        numpanls=numpanls+1;
        novary=novary.* dotfeq(minc(x),maxc(x)) ; @ used to create xid,zid @
        if inpanel/=oldsize ;                       @dont use prior matrix @
           tranmtx=eye(inpanel)- (ones(inpanel,inpanel)/inpanel);
           oldsize=inpanel;
           maxsize=maxc(inpanel|maxsize);
        endif;                             @else, use prior value of tranmtx@
        xnew=tranmtx*x; ynew=tranmtx*y; foo=foo+1;
        {xx,xy,yy,ysum,nobs}=makemom(xnew,ynew,xx,xy,yy,ysum,nobs);
    endo;

   {xid,zid}=getxzid(novary);           @ use panel variant variables @
   if scalmiss(xid);
        " No panel variant variables. ";
         retp(0);
    endif;

    xx=xx[xid,xid];
    xy=xy[xid];

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

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

     output on ;
     ? ;
    " FE estimates: ";
    format /lz 10,5 ; b' ;?;
    output off ;

@ We have initial, consistent beta.  Now compute SIGMA, the estimate of the
VC of Qu (Q being the demeaner matrix).  Note that we assume that for panels of
differing sizes, all common elements are the same.  In other words, the 3x3
upper left sub matrix of a 5x5 matrix (from a 5 observation panel) will have
the same elements as the 3x3 matrix from a 3 observation panel. @

   nextobs=1; tranmtx=0;oldsize=0;
   foo=seekr(f1,1);              @ rewind @
   yytot=zeros(maxsize,maxsize) ;
   foo=0;

   do until eof(f1) or scalmiss(nextobs) ;

       if abortreg(foo,atline,atcol)==1;
           retp(0);
       endif;

        {y,x,foo,nextobs}=grabpanl(f1,nextobs);
        if ismiss(y)==1 ;                  @ if empty, skip to next panel @
           continue;
        endif;
        inpanel=rows(y); foo=foo+1;
        if inpanel/=oldsize ;                       @dont use prior matrix @
           tranmtx=eye(inpanel)- (ones(inpanel,inpanel)/inpanel);
           oldsize=inpanel;
        endif;                             @else, use prior value of tranmtx@
        yy=moment((y-x[.,xid]*b)',0) ;
        yytot[1:inpanel,1:inpanel]=yytot[1:inpanel,1:inpanel]+yy ;
     endo;
     qyvar=yytot/numpanls ;

@ Now do GLS using QYVAR.  Again, we use submatrices of QYVAR when panel
  size is less than largest panel @

   nextobs=1; tranmtx=0;oldsize=0;
   foo=seekr(f1,1);              @ rewind @
   xx=0;xy=0;
   foo=0;
   do until eof(f1) or scalmiss(nextobs) ;

       if abortreg(foo,atline,atcol)==1;
           retp(0);
       endif;
        {y,x,foo,nextobs}=grabpanl(f1,nextobs);
        if rows(y)==1 ;                  @ if empty, skip to next panel @
           continue;
        endif;

        inpanel=rows(y);
        foo=foo+1;
        if inpanel/=oldsize ;                       @dont use prior matrix @
           tranmtx=eye(inpanel)- (ones(inpanel,inpanel)/inpanel);
           oldsize=inpanel;
           vuse=eye(inpanel);
           vuse[1:inpanel-1,1:inpanel-1]=qyvar[1:inpanel-1,1:inpanel-1] ;
           vuseinv=invpd(vuse);
           vuseinv[inpanel,inpanel]=0;
        endif;                             @else, use prior value of tranmtx@
        xx=((tranmtx*x)'*vuseinv*(tranmtx*x))+xx ;
        xy=((tranmtx*x)'*vuseinv*(tranmtx*y))+xy ;

     endo;
     vc=invpd(xx[xid,xid]);
     b=vc*xy[xid,1] ;
     stderr=sqrt(diag(vc));
     tstat=b./stderr ;
     pvt=cdftc(abs(tstat),df);
output on ;
?;
"------- Kiefer  Results---- ";
call wrtetitl(_titles);

"                 ____________________________________";
"                    Data set          :  " $dataset;
"                    Dependent variable:  " $depvar;
"                 ------------------------------------";?;
     d1=datestr(0); d2=timestr(0);
     "Date        :" $d1 " , " $d2 ;
     FORMAT /ld 8,0;
     "Observations used:" nobs   "  (in dataset:" totobs ")";
     " # of panels:   " numpanls "   Degrees of freedom:" df;

/*
     format /lo 11,7;
     "R-squared   :" rsq;;  "     Rbar-squared       :" rbsq;
     format /lo 11,7;
     "Residual SS :" ssr;;  "   Std error of est (å) :" sigma ;
     "Total SS    :" sst;;
     format 1,0; "     F(" k "," df ")=";; format 6,4; fstat;;
     format 4,2; /rd "  P-value=" pvf;
     format /lo 14,9;  "Log-Likelihood:  " loglk;
*/

     print;
     "Var             Coef           Std. Error     t-Stat         P-Value";
     "--------------------------------------------------------------------";
omat=parname1[xid]~b~stderr~tstat~pvt;
mask=0~1~1~1~1;
let fmt[5,3]="-*.*s" 10 8
             "*.*lg " 14 6
             "*.*lg " 14 6
             "*.*lg " 14 6
             "*.*lg " 11 3;
call printfm(omat,mask,fmt);
print;
output off ;

clearg beta;
parnames=parname1;  beta=b ;
    retp(1);
 endp;



/**********/
@ proc to create pointers to elements that are panel variant (Xid) and to
  elements that are panel invariant (Zid) @

   proc (2)=getxzid(novary);
     local zid,xid;

@ create pointer to vars that are panel invariant @
      zid=seqa(1,1,rows(novary)).*novary;
      zid=packr(miss(zid,0));

@ create pointer to vars that are panel variant @
      xid=seqa(1,1,rows(novary)).* .not(novary) ;
      xid=packr(miss(xid,0));

      retp(xid,zid);
  endp;

/*****************/
 proc chosemtx(inpanel);
 local tmp,ii,psi2;


 if inpanel<=1;      @only 1 obs in this panel: no transformation possible@
    retp(1);            @ so leave panel as is @
 endif;

    if inpanel < 8;
        retp(diffmtx[1:inpanel-1,1:inpanel]);
    else;                               @ big panel, make a diffmtx@
        tmp=-eye(inpanel);
        ii=0; do until ii==inpanel-1;
          ii=ii+1; tmp[ii,ii+1]=1;
        endo;
        retp(tmp[1:inpanel-1,1:inpanel]);
    endif;


endp;
/********************************/
/********************************
NLSTPRNT -- Proc to print the final results for non-linear statistical
            estimation. From MAXPROC3.ARC

            Call as:
          call nlstprnt(atitle,dataset,depvar,nobs,df,vof,voftitle,parnames,
              b,stderr,t,pvt,cvid,iter);

            INPUTS
              atitle -- title, used by wrtetitl
            dataset -- the name of the data set used; if null, will not
                       print anything for this
             depvar -- the name of the dependent variable; if null will not
                       print anything for this
               nobs -- number of observations
                 df -- degress of freedom
                vof -- value of the objective function at the optimum
           voftitle -- a string to be printed that describes the objective
                       function (e.g. VOFTITLE="LOG LIKELIHOOD FUNCTION";)
           parnames -- a character vector containing the names to be
                       associated with the parameters.
                  b -- the kx1 vector of estimated parameters
             stderr -- the kx1 vector of standard errors of the estimated
                       parameters.
                  t -- the kx1 vector of t-ratios of the parameters.
                pvt -- the kx1 vector of p-values of the t-ratios.
                cvid  -- Type of cv matrix
                iter  -- # of iterations
............................................................................*/

proc nlstprnt(atitle,dataset,depvar,nobs,df,vof,voftitle,parnames,
              b,stderr,t,pvt,cvid,iter);
local depflag, tme, hr, min, dt, yr, mo, day, dte, tmv, omat, dfmt, tfmt,
mask, fmt;
clear depflag;

/* convert to string */
depvar=""$+depvar;

/* convert to upper case */
dataset=upper(dataset);
depvar=upper(depvar);
parnames=upper(parnames);

print; format /ro 8,8;
call wrtetitl(atitle);
?;
"                 ____________________________________";
"                    Data set          :  " $dataset;
"                    Dependent Variable:  " $depvar;
"                 ------------------------------------";?;

     call datetime;
  format /rd 6,0 ;
 " # iterations:" iter ;;

 "  # Observations:" nobs;; "     D.F. : " df;
 format 10,6;  print;

  "         ***   " $voftitle vof "  *** ";  ?;


"Var          Coef       Std. Error    t-Stat  P-Value ";
@    Mean     Elasticity@
"---------------------------------------------------------------------";
omat=parnames~b~stderr~t~pvt;        @ form output matrix to be printed    @
mask=0~1~1~1~1 ;    @  mask to print first column as strings @

/**
omat=omat~bars~elasts ;
mask=0~1~1~1~1~1~1   @mask to print first column as strings @
*/

let fmt[5,3]="-*.*s" 9 8    @ format info for variable names column @
             "*.*lg " 12 7   @ format info for coef column           @
             "*.*lg " 11 6   @ format info for std. error column     @
             "*.*lg " 9 4   @ format info for T-stat column         @
             "*.*lg " 9 4 ; @ format info for P-value column        @

/*
             "*.*lg " 10 5  @ format info for means column        @
             "*.*lg " 9 4;  @ format info for elasticies column        @
*/

call printfm(omat,mask,fmt); @ print omat using above format         @

if cvid==1;
      " Hessian used for CV  ";
elseif cvid==2;
      " BHHH matrix used for CV ";
elseif cvid==3;
      " 'White heteroscedastic consistent' matrix used for CV ";
else;
     " CV matrix unobtainable: identity matrix used ";
endif;


retp("");
endp;







/******************* PRoc PRNTOLS1 *****************/
@ PRNTOLS1 -- Print results of OLS regression.
  Uses global variables, and computes other variables.
  Status=1 if okay, 0 otherwise

  Global vars used are:     ssr,varhat,xxinv,vc,b,ybar,sst,df,nobs
  Globals set are:      sigma,stderr,t,pvt;
  @

proc PRNTOLS1(parname1,olstitle);

    local rsq,fstat,rm,pvf,k,rbsq;
    local d1,d2;
    local fmt,omat,mask;

      parnames=parname1;                @ save global version @
      k=rows(parname1);
      sigma=sqrt(varhat);
      stderr=sigma.*sqrt(diag(xxinv));
      t=b./stderr;
      pvt=2*cdftc(abs(t),df);

      rsq=1-(ssr/sst);           @ R square stat @
      rbsq=1-(ssr/df)/(sst/(nobs-1));


      if qconst == 1 and k>1;             @ compute f stat @
         fstat=((sst-ssr)/(k-1))/(ssr/df);
      else;
        rm=eye(k);
        fstat=b'inv(xxinv)*b*df/(ssr*k);
      endif;

     if qconst==1 and k==1;     @ prob value of f stat @
        pvf=miss(1,1);
     else;
         if fstat> 0 ;
         pvf=cdffc(fstat,k,df);
         else ; pvf=miss(1,1);
         endif;
     endif;


@ Display title, ancillary info @

"------- OLS  Results---- ";
call wrtetitl(_titles);
olstitle ;                        @ model title @


"                 ____________________________________";
"                    Data set          :  " $dataset;
"                    Dependent variable:  " $depvar;
"                 ------------------------------------";?;

     d1=datestr(0); d2=timestr(0);
     "Date        :" $d1 " , " $d2 ;
     FORMAT /ld 8,0;
     "Observations used:" nobs   "  (in dataset:" totobs ")";
     " # of panels:   " numpanls "   Degrees of freedom:" df;
     format /lo 11,7;
     "R-squared   :" rsq;;  "     Rbar-squared       :" rbsq;
     format /lo 11,7;
     "Residual SS :" ssr;;  "   Std error of est (å) :" sigma ;
     "Total SS    :" sst;;
     format 1,0; "     F(" k "," df ")=";; format 6,4; fstat;;
     format 4,2; /rd "  P-value=" pvf;
     format /lo 14,9;  "Log-Likelihood:  " loglk;
     print;
     "Var             Coef           Std. Error     t-Stat         P-Value";
     "--------------------------------------------------------------------";
omat=parname1~b~stderr~t~pvt;
mask=0~1~1~1~1;
let fmt[5,3]="-*.*s" 10 8
             "*.*lg " 14 6
             "*.*lg " 14 6
             "*.*lg " 14 6
             "*.*lg " 11 3;
call printfm(omat,mask,fmt);
print;


retp(0);
endp;

/**********/
@ Compute ols coefficients, etc.
 Variables returned to globals:  ssr,varhat,xxinv,vc,b,ybar,sst,loglk
@

 proc compols(xx,xy,yy,ysum,nobs,df);

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

@ compute other stuff @
     ssr=yy-b'xy;
     varhat=ssr/df;
     xxinv=invpd(xx);
     vc=varhat*xxinv;
     ybar=ysum/nobs;
     sst=yy-nobs*ybar*ybar;

     loglk=((-nobs/2)*ln(2*pi))-((nobs/2)*ln(varhat))-(ssr/(2*varhat));

     retp(1);

endp;

/*******      LISTVARS
@ List variables to output (all of them are globals @

proc listvars(dataset);

        local mask,fmt;
        local allnames,foo,i,goo ;
        format 1,8;?;

         if scalmiss(yindex);
             " No variables selected !! ";

         else;
            output on ;
             " ------------" ;
             " Displaying all selected variables: ";
             "   DATASET: " $dataset ;
             "   Dependent variable: " $depvar ;;
             if not(scalmiss(yearid));
                "   Year variable: " $yearname ;
             endif;
             IF NOT(SCALMISS(PANELID));
                IF ISNUMBER(PANLNAME);
                   " Balanced Panels of size: " $panlname ;
                else ;
                   " Panel Variable: " $panlname ;
                endif;
             endif;
             If qconst ; " Constant included " ;
                 else ; " No constant " ; endif;
             " Number of independent variables: ";; format /ld 5,0; nindep ;
             if nindep>0; format 1,8;
                 " Independent vars:" $xname';
                 ALLNAMES=xname;
             endif;

        retp(0);
   endp;

*****/


/***********/
@ Add to moment, etc matrices, given x,y and old values @

     proc (5)=makemom(x,y,xx,xy,yy,ysum,nobs);
@ retp new values of xx,xy,xy,ysum,nobs @
        retp(xx+(x'x),xy+(x'y),yy+(y'y),ysum+sumc(y),nobs+rows(y));

     endp;


/****************************/
@ Little proc to drive generic MKPARNMS @
  proc mkparnm2;
    local parname1;
    parname1=mkparnms(nindep,qconst,xname,dumnames,polyname,miss(1,1));
 retp(parname1);
endp;


/*****************************/
@ compute product from j=0 to t-m of P[j-lag].
  Vector returned will correspond to values for t=m..rows(p) @

  proc mkppm(p,m,lag);
  local tr,tth,ppm;
  tr=rows(p);
  ppm=ones(rows(p),1);

  tth=m;
  do until tth>tr ;
     ppm[tth]=prodc(p[m-lag:tth-lag]);
     tth=tth+1;
  endo;
  retp(ppm);
  endp;


/**********************************/
@ Read next panel. If balanced, just get next -panelid obs.
 If unbalanced, continue until id changes.
       {ynew,xnew,yrvals,nextobs}=grabpanl(f1,nextobs);
  This is used only by setiore and setiod2
@
 proc (4)=grabpanl(f1,nextobs);

local id,a1,tmp,plen,pat,dta,thisp,m,rws,thisread,iread;

     if nextobs== 0;
       retp(miss(1,1),0,0,nextobs);          @ miss = not valid panel @
     endif;

    if panelid<0 ;                @ panelid is global. balanced panel? @
         nextobs=nextobs-panelid ;
         dta=readr(f1,-panelid);
         plen=-panelid ;
         if eof(f1); nextobs=miss(1,1); endif;
   else ;                           @  unbalanced @
        dta=miss(1,1);
        iread=5;        @ could be increased if bigger panels on average @
morepanl:;                          @ jump here on long panels @
        nextobs=seekr(f1,nextobs);
        a1=readr(f1,iread);            @ probably enuf @
        thisread=rows(a1);
        if scalmiss(dta);
          id=a1[1,panelid];           @ panel defined by first obs @
        endif;
        thisp=seqa(1,1,rows(a1)).*(a1[.,panelid].==id) ;
        thisp=packr(miss(thisp,0));

        if scalmiss(thisp);       @ this can happen if 2nd read @
                 ;       @e.g. -- if panel of iread length read on 1st read @

        else ;
           rws=rows(thisp);
           if rws==1 ;             @ 1 element panel @
               m=1; thisp=1;
           else;                   @ more than 1, just take consecutive ones@
              tmp=thisp[2:rws]-thisp[1:rws-1]; m=rows(thisp);
              if maxc(tmp)>1 ;    @ uh oh, not all members of panel together?@
                 m=1 ;         @ then assume more than one panel with same id @
                 do until 1==2 ;  @ break should always work @
                    if tmp[m]>1 ; break  ; endif;
                    m=m+1;               @nb: first obs always in panel @
                 endo;
                 thisp=thisp[1:m];
              endif;                @ maxc(tmp) @
           endif;                  @ rws==1 @
           if scalmiss(dta);
               dta=a1[thisp,.];
           else ;
              dta=dta|a1[thisp,.];
           endif;
           nextobs=nextobs+m;
           plen=rows(dta);
           if m eq thisread ;     @ uh oh, panel may be longer then 25 elements@
              if not(eof(f1));
                 goto morepanl ;
              else;
                 nextobs=miss(1,1);
              endif;
           endif;
       endif;               @ scalmiss(thisp) @
   endif;                     @ balanced @


   call defvar2(dta,sortopt);             @ sortopt is global @

   if rows(y) eq plen ;              @ y,x,yearvals are global @
      retp(y,x,yearvals,nextobs);
   else ;
      retp(miss(1,1),0,0,nextobs);          @ miss = not valid panel @
   endif;

   endp;


/****************************************************************/
@ Pull out wanted variables. Assume dta contains 1 panel @
@ Note: selection on DEFVAR to assure that years, etc. match up @
@ Note: sort, if desired @

    proc defvar2(dta,sortopt);

       local m1,tmpc,tmp,nn,ifoo;

       x=miss(1,1); y=x;

       dta=packr(doselct(dta,critlist));          @select obs @
       if scalmiss(dta);
          retp(0);
       endif;

       if sortopt$=="ON" and not(scalmiss(yearid));            @ sort panel @
              dta=sortc(dta,yearid);
       endif;

       nn=rows(dta);

       tmpc=critlist;                    @critlist is global @
       critlist=miss(1,1);                @a hack to speed up defvar@
       call defvar(dta);                 @ pull out generic vars @
       critlist=tmpc;                    @ restore @


       if panelid>0;                     @ create panelid variable @
          panlvalu=dta[1,panelid];
       else;
          panlvalu=miss(1,1);           @ create later @
       endif;

       if not(scalmiss(yearid));         @ create year var @
           yearvals=dta[1:nn,yearid];
       else;
          yearvals=seqa(1,1,rows(dta));
       endif;

       yearvals=substute(yearvals,yearvals.==0,1);

      retp(1);

    endp;


/*********/
@ display info, possibly abort regression @
   proc abortreg(ithpan,atline,atcol);
       local akey;

       if ithpan % 25 == 0 ;
          locate atline,atcol ;
           format /ld 6,0;
           output off;
           " # of Panels read=" ithpan  " ( ESC to abort regression ) ";
            akey=key;
            if (akey==27) ;
                   " Aborting regression ..... ";
                   output on ;
                   retp(1);
             endif;
             output on ;
             retp(0);
        endif;
        retp(0);
   endp;



@ mle procs for serpanel program @


/*************/
@ 2 procs used by gradcvd2 @

proc (3)=dprodm(ps,dps1,dps2,a,b,sub);

  local aprod,l1,l2,tmp1,tmp2,pz ;
  local tmpp,j1;

  l1=sub-b ;
  l2=sub-a ;

  pz=sumc(ps[l1:l2].==0);

  if pz gt 1  ;
      retp(0,0,0);
  endif;

  j1=l1 ; tmp1=0;tmp2=0;
  do until j1>l2 ;
      tmpp=ps ; tmpp[j1]=1 ;
      aprod=prodc(tmpp[l1:l2]);
      tmp1=tmp1+(aprod*dps1[j1]);
      tmp2=tmp2+(aprod*dps2[j1]);
      j1=j1+1;
  endo;
  aprod=prodc(ps[l1:l2]) ;

  retp(aprod,tmp1,tmp2);

  endp;

proc (3)=dprod(ps,dps1,dps2,a,b);

  local aprod,tmp,l1,l2,tmp1,tmp2,j1,pz,tmpp ;

  pz=sumc(ps[a:b].==0);

  if pz gt 1 ;
      retp(0,0,0);
  endif;

  j1=a ; tmp1=0;tmp2=0;
  do until j1>b ;
      tmpp=ps ; tmpp[j1]=1 ;
      aprod=prodc(tmpp[a:b]);
      tmp1=tmp1+(aprod*dps1[j1]);
      tmp2=tmp2+(aprod*dps2[j1]);
      j1=j1+1;
  endo;
  aprod=prodc(ps[a:b]) ;


  retp(aprod,tmp1,tmp2);

  endp;

/*********************************/
@ compute 2 derivs of a PROD x, given x and dx and dx2, and start/end in x,etc.
  Note that we assume each element of x is a function of some parameter
  z, so that we are compute d (prod(x)) / d z --
  dx and dx2 contains a T vector, with element t = d X[t]/ d z.

Also, return prod(x) (since its calcluated anyways @

   proc (3)=dprodX(x,dx,dx2,i1,i2);

      local j2,amtx,prodxn,tmp,amtx2,amtx1 ;

      j2=rows(x);

      if ( rows(x) ne rows(dx)) or (i1<0) or (i2>j2) or (i1>i2) ;
          " % Bad Input to DprodX " ; stop ;
      endif ;

@ extract subvector @
      x=x[i1:i2];
      dx=dx[i1:i2] ;
      dx2=dx2[i1:i2] ;
      j2=rows(x);

      if j2==1 ;        @ simple case @
        retp(x,dx,dx2);
      endif;

@ The form of deriv is
      d1 * prod,not_1  + d2* prod,not_2 + .. + dj2 * prod,not_j2
  where dn is the derivative of the nth term of p = dp[n]
  and prod,not_2 is the product of p, not including n @

       amtx=ones(1,j2).*x ;
       amtx1=diagrv(amtx,dx);
       amtx2=diagrv(amtx,dx2);

       prodxn=prodc(amtx1~amtx2~x);

       retp(prodxn[(2*j2)+1],sumc(prodxn[1:j2]),sumc(prodxn[(j2+1):(2*j2)]) );
  endp;

