/*
**                         FILE BPM8.ARC
**
** This file contains the bivariate probit model procedures.
** Part 8.
*/



/* Declaration of globals. */

declare _rpar ?=0;
declare _sfile ?=0;
declare _final ?=1;
declare __miss !=1;
declare _eps ?=0.001;
declare _maxit ?=30;
declare _file ?=0;
declare __miss !=1;
declare __output !=0;
declare _full ?=1;


/* Setting of libraries needed. */

library a:\bpm,stat,pgraph,dutil;




/*
** BPROBSEL
**
** FORMAT    {indv,teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=
**                          bprobsel(dataset,indvar,depvar,sig);
**
** INPUT      dataset  : String, the name of the dataset that
**                       contains the variables.
**            indvar   : An Lx1 character vector containing
**                       the names of the independent variables.
**                       case a unit column is appended and L-1 is incremented).
**            depvar   : A 2x1 character vector containing the names
**                       of the dependent variables.
**            sig      : Bound for chi-square test.
**
** OUTPUT     indv     : Character vector, names of the independent variables
**                       selected.
**            teta     : (2*L+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the bivariate parameters.
**            ll       : Bivariate log-likelihood.
**            tetra    : 3x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the tetrachoric parameters.
**            llt      : Tetrachoric log-likelihood.
**            teta1    : (L+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the univariate parameters (first variable).
**            ll1      : Univariate log-likelihood (first variable).
**            teta2    : (L+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the univariate parameters (second variable).
**            ll2      : Univariate log-likelihood (second variable).
**
** GLOBALS    _sfile  0: Output device is file (default).
**                    1: Output device is printer.
**            _final  0: No final analysis done.
**                    1: Final analysis done (default).
**            _eps     : Tolerance for the gradient (default 0.001).
**            _maxit   : Upper bound for the number of iterations
**                       (default 30).
**
** REMARK     BPROBSEL calls BPROBIT. It starts with one independent variable,
**            according to a criterion, other independent variables are added.
**            The global _rpar, controlling wich algorithm is called by BPROBIT
**            (BIVPROB, BIVPROBR, PROFLIK) should be 0 or 1. PROFLIK should
**            not be called by BPROBSEL, as in that case no information is
**            returned to feed the test subprocedure of BPROBSEL.
*/


proc (9)= bprobsel(dataset,indvar,depvar,sig);

local nin,inlbl,nout,outlbl,outfile,indv,tel,schat0,w0,een0,twee0,test,handle,
      akkoord, schat1,w1,een1,twee1,lvec,mask,fmt,cnt,ll1,ll2,tetra,teta1,
      tijd1,tijd2,teta2,pending,stoppen,llt,step,chosen,fmt2,mask2,aantal;


tijd1=hsec;

/* Creating/resetting the output file. */

if _sfile==0;
   output file=bprobsel.out reset;
else;
   output file="lpt1" reset;
endif;


/* Avoiding that PROFLIK is called from within BPROBSEL. */

label:
if _rpar==2;
   errorlog "You shouldn't try to call PROFLIK from BPROBSEL. Do you want";
   errorlog "rho parametrized (1) or not (0) ?";
   _rpar=con(1,1);
   goto label;
endif;


/*
** Defining the working data set, and creating a dataset with the variables
** passed in indvar and depvar.
*/

outfile="bprobsel";
{nin,inlbl,nout,outlbl}=dselvar(dataset,outfile,indvar|depvar);
open handle=^outfile for read;
aantal=rowsf(handle);
handle=close(handle);


/* Printing the title and the names of te variables. */

cls;
print " ";
print "            BIVARIATE  PROBIT  ESTIMATION :  SELECTION";
print "            ==========================================";
print " ";
print " ";

format 14,4;
print "Dependent variables :";
print "---------------------";
print $depvar';
print " ";
print " ";
print "Candidate predictors :";
print "----------------------";
print $indvar';
print " ";
print " ";
format 10,4;
print "         Criterion :    Likelihood ratio";
print "Significance level :" sig;
print "   Number of cases :" aantal;
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";


/* Scanning of individual independent variables. */

_file=2;
tel=1;
lvec=zeros(rows(indvar),3);
do while tel<=rows(indvar);
   output off;
   _full=0;
   {schat0,lvec[tel,1],tetra,llt,teta1,ll1,teta2,ll2}=
   bprs(outfile,indvar[tel],depvar);
   _full=1;
   output on;
   tel=tel+1;
endo;

tel=1;
do while tel<=rows(indvar);
   lvec[tel,2]=lvec[tel,1]-llt;
   if _rpar==0;
      lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),2);
   elseif _rpar==1;
      lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),3);
   endif;
   tel=tel+1;
endo;



/* Ordering the independent variables. */

lvec=-lvec~indvar;
lvec=sortc(lvec,1);
indvar=lvec[.,4];
lvec=-lvec[.,1:3];


print "STEP  1";
print "=======";
print " ";
mask=0~1~1~1;
let fmt[4,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4
"*.*lf " 14  4;
fmt2=fmt[1 2,.];
mask2=0~1;
chosen="intercept"~llt~0~0;
print "Variables in model :";
print "--------------------";
print $chosen[.,1]';
print " ";
print " ";
print "Candidates :";
print "------------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";
cnt=printfm(indvar~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
print " ";
if lvec[1,3]<sig;
   print $indvar[1] "will be kept in the model.";
else;
   print "No more variables will be taken in the model.";
endif;
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";



/* Selection procedure. */

llt=0;
akkoord=(lvec[1,3]<=sig);
indv=indvar[1];
chosen=chosen|(indv~lvec[1,.]);
pending=indvar[2:rows(indvar)];
tel=1;
step=1;

stoppen=not akkoord;
do until stoppen;
   tel=1;
   step=step+1;
   w0=lvec[1,1];
   lvec=zeros(rows(pending),3);
   format 3,0;
   print "STEP" step;
   print "=======";
   print " ";
   print "Variables in model :";
   print "--------------------";
   format 8,0;
   print $chosen[.,1]';
   print " ";
   print "Candidates :";
   print "------------";
   do while tel<=rows(pending);
      output off;
      _full=0;
      {schat0,lvec[tel,1],tetra,llt,teta1,ll1,teta2,ll2}=
      bprs(outfile,indv|pending[tel],depvar);
      _full=1;
      output on;
      lvec[tel,2]=lvec[tel,1]-w0;
      if _rpar==0;
         lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),2);
      elseif _rpar==1;
         lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),3);
      endif;
      tel=tel+1;
   endo;

   lvec=-lvec~pending;
   lvec=sortc(lvec,1);
   pending=lvec[.,4];
   lvec=-lvec[.,1:3];

   print " ";
   print "      variable        log-lik      chisquare        p-value";
   print "      --------        -------      ---------        -------";
   cnt=printfm(pending~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
   print " ";
   if rows(pending)==1;
      print " ";
   endif;
   akkoord=(lvec[1,3]<=sig);
   if akkoord;
      print $pending[1] "will be kept in the model.";
   else;
      print "No more variables will be taken in the model.";
   endif;
   print " ";
   print "************************************************************"\
         "*****************";
   print " ";

   if akkoord;
      chosen=chosen|(pending[1]~lvec[1,.]);
      indv=indv|pending[1];
      if rows(pending)==1;
         stoppen=1;
      else;
         pending=pending[2:rows(pending)];
      endif;
   else;
      stoppen=1;
   endif;

endo;


output on;

print "Summary :";
print "---------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";

if rows(chosen)==2 and chosen[2,4]>sig;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
else;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
   print " ";
   cnt=printfm(chosen[2:rows(chosen),.],mask,fmt);
endif;
print " ";
print " ";
print "************************************************************"\
      "*****************";


print " ";
tijd2=hsec;
tijd2=round((tijd2-tijd1)/100);
format /rd 5,0;
print "Computation time :" round(tijd2/60) "minutes" tijd2%60 "seconds.";
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";

if _final;
   print "We now present the analysis for selected independent variables.";
   print " ";
   if _sfile==1;
       lprint "\f";
       _file=1;
   else;
      _file=0;
   endif;
   {schat0,w0,tetra,llt,teta1,ll1,teta2,ll2}=bprs(outfile,indv,depvar);
else;
   output off;
endif;

retp(indv,schat0,w0,tetra,llt,teta1,ll1,teta2,ll2);
endp;


