/*
**                         FILE BPM14.ARC
**
** This file contains the bivariate probit model procedures.
** Part 14.
*/



/* Declaration of globals. */

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


/* Setting of libraries needed. */

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




/*
** BPROBSC
**
** FORMAT    {indv,info}=bprobsc(dataset,indvar,depvar,sig);
**
** INPUT      dataset  : String, the name of the dataset that
**                       contains the variables.
**            indvar   : An Lx1 or (L-1)x1 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 of chosen variables.
**            info     : rows(indv)x3 vector. Line j contains the log-
**                       likelihood with variables indv[1:j] included.
**                       info[j,2 3] contain the chisquared and P-value
**                       comparing the loglikelihood for variables 1 through
**                       j with the loglikelihood for variables 1 through
**                       j-1.
**
** 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     BPROBSC calls BPRS. 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 (2)= bprobsc(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,t,k,
      tijd1,tijd2,teta2,pending,stoppen,llt,step,chosen,fmt2,mask2,aantal,s,j,
      indexi,indexd,nvec,hulp1,i,hulp2,g1,g2,hulp,sortv,l,afg1,afg2,rval,sval,
      grad,hess,xx,yy,x,y,ll,var01,xh,yh,sh,minima,maxima,r,s,volgorde;


tijd1=hsec;

/* Creating/resetting the output file. */

_file=2;
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 :    Score test";
print "Significance level :" sig;
print "   Number of cases :" aantal;
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";




if dataset$=="";

   x=indvar;
   y=depvar;
   let depvar="y1" "y2";
   let indvar="constant" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;

else;

   {indvar,indexi}=indices(dataset,indvar);
   {depvar,indexd}=indices(dataset,depvar);
   nvec=setvars(dataset);
   open handle=^dataset for read;
   do until eof(handle);
      hulp1=readr(handle,100);
      i=1;
      do while i<=rows(indvar);
         hulp2=varget(indvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexi[i]];
         else;
            hulp2=hulp1[.,indexi[i]];
         endif;
         cnt=varput(hulp2,indvar[i]);
         i=i+1;
      endo;

      i=1;
      do while i<=rows(depvar);
         hulp2=varget(depvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexd[i]];
         else;
            hulp2=hulp1[.,indexd[i]];
         endif;
         cnt=varput(hulp2,depvar[i]);
         i=i+1;
      endo;

   endo;
   handle=close(handle);
   y=varget(depvar[1])~varget(depvar[2]);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;

   if not(x[.,1]==ones(rows(x),1));
      indvar="constant"|indvar;
   endif;

   nvec=setvars(dataset);

endif;


/* Checking for faulty data and deletion of rows with missing values. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=2);
   errorlog "y does not stand for a table of 0-1 variables.";
   end;
endif;

t=x~y;
_cases=rows(t);
g1=cols(x);
g2=cols(t);
t=packr(t);
_cases=_cases|rows(t);
y=t[.,g1+1:g2];
x=t[.,1:g1];


/* Determining the shape of the tables. */


minima=minc(y);
maxima=maxc(y);
i=minima[1];
rval=0;
do while i<=maxima[1];
   if not y[.,1]/=i;
      rval=rval|i;
   endif;
   i=i+1;
endo;
rval=rval[2:rows(rval)];
i=minima[2];
sval=0;
do while i<=maxima[2];
   if not y[.,2]/=i;
      sval=sval|i;
   endif;
   i=i+1;
endo;
sval=sval[2:rows(sval)];

r=rows(rval);
s=rows(sval);


/* Transforming the dependent variables into contingency tables. */

hulp=zeros(rows(y),r*s);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      hulp[.,(i-1)*s+j]=((y[.,1] .== rval[i]) .and (y[.,2] .== sval[j]));
      j=j+1;
   endo;
   i=i+1;
endo;

y=hulp;
clear hulp;


/* Multiple columns sort on the independent variable(s). */

g1=cols(x);
g2=g1+cols(y);
t=x~y;
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);


/* Summing over equal independent variable vectors. */

i=1;
k=1;
sh=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      sh[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      sh[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

sh=sh[1:k-1,.];


/*
** Separating dependent and independent variables and adding a column
** of ones to the matrix of independent variables if necessary.
*/


y=sh[.,g1+1:g2];
x=sh[.,1:g1];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;


/* Construction of auxiliary matrices containing only the (0,1)-covariates. */

var01=zeros(cols(x),1);
i=1;
do while i<=cols(x);
   var01[i]=sumc((x[.,i].==1)+(x[.,i].==0))==rows(x);
   i=i+1;
endo;

{xh,yh}=shrink(x.*var01',y);


/* Intercept. */

output off;
_restric={1 1,0 0};
if var01[2]==0;
   {xx,yy}=shrink(x[.,1 2],y);
else;
   {xx,yy}=shrink(xh[.,1 2],yh);
endif;
_restric=1;
{schat1,ll,schat0,llt,teta1,ll1,teta2,ll2}=
bprobrs(xx,yy,r,s,indvar[1 2],depvar);
_restric=1;
output on;


/* Scanning of individual independent variables. */

tel=2;
lvec=zeros(rows(indvar)-1,2);
schat0=schat0[.,1];
l=(rows(schat0)-r-s+3)/2;
schat1=schat0[1:r+l-2]|0|schat0[r+l-1:2*l+r+s-4]|0|schat0[2*l+r+s-3];
schat1=transoe(schat1',r,s)';
do while tel<=rows(indvar);
   if var01[tel]==0;
      {xx,yy}=shrink(x[.,1 tel],y);
   else;
      {xx,yy}=shrink(xh[.,1 tel],yh);
   endif;
   afg1=transeo(gradbprs(schat1,xx,yy,r,s)',r,s,1)';
   afg2=transeo(transeo(hessbprs(schat1,xx,yy,r,s),r,s,1)',r,s,1)';
   lvec[tel-1,1]=afg1'*inv(-afg2)*afg1;
   tel=tel+1;
endo;

tel=1;
do while tel<=rows(lvec);
   lvec[tel,2]=cdfchic(lvec[tel,1],2);
   tel=tel+1;
endo;


/* Ordering the independent variables. */

volgorde=seqa(2,1,rows(indvar)-1);
lvec=-lvec~volgorde;
lvec=sortc(lvec,1);
volgorde=lvec[.,3];
lvec=-lvec[.,1:2];


print "STEP  1";
print "=======";
print " ";
print " ";
mask=0~1~1;
let fmt[3,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4;
chosen="intercep"~llt~0~0;
print "Variables in model :";
print "--------------------";
format 10,4;
print $chosen[.,1]';
print " ";
print " ";
print "Candidates :";
print "------------";
print " ";
print "      variable      chisquare        p-value";
print "      --------      ---------        -------";
cnt=printfm(indvar[volgorde]~lvec[.,1]~(lvec[.,2]),mask,fmt);
print " ";
akkoord=(lvec[1,2]<=sig);
if akkoord;
   print $indvar[volgorde[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;
indvar=indvar[1]|indvar[volgorde];
x=x[.,1]~x[.,volgorde];
xh=xh[.,1]~xh[.,volgorde];
var01=var01[1]|var01[volgorde];
indv=indvar[1 2];
chosen=chosen|(indv[2]~0~lvec[1,.]);
if rows(indvar)>=3;
   pending=indvar[3:rows(indvar)];
   stoppen=not akkoord;
else;
   stoppen=1;
endif;
tel=1;
step=1;

do until stoppen;
   tel=1;
   step=step+1;
   output off;
   if var01[1:rows(chosen)]/=0;
      {xx,yy}=shrink(xh[.,1:rows(chosen)],yh);
   else;
      {xx,yy}=shrink(x[.,1:rows(chosen)],y);
   endif;
   {schat0,ll,tetra,llt,teta1,ll1,teta2,ll2}=
   bprobrs(xx,yy,r,s,chosen[.,1],depvar);
   output on;
   chosen[rows(chosen),2]=ll;
   w0=ll;
   lvec=zeros(rows(pending),2);
   format 3,0;
   print "STEP" step;
   print "=======";
   print " ";
   print "Variables in model :";
   print "--------------------";
   format 10,4;
   print $chosen[.,1]';
   print " ";
   print "Candidates :";
   print "------------";
   l=(rows(schat0)-r-s+3)/2;
   schat1=schat0[1:r+l-2,1]|0|schat0[r+l-1:2*l+r+s-4,1]|0|
   schat0[2*l+r+s-3,1];
   schat1=transoe(schat1',r,s)';
   do while tel<=rows(pending);
      if var01[1:rows(chosen) rows(chosen)+tel]/=0;
         {xx,yy}=shrink(xh[.,1:rows(chosen) rows(chosen)+tel],yh);
      else;
         {xx,yy}=shrink(x[.,1:rows(chosen) rows(chosen)+tel],y);
      endif;
      afg1=transeo(gradbprs(schat1,xx,yy,r,s)',r,s,1)';
      afg2=transeo(transeo(hessbprs(schat1,xx,yy,r,s),r,s,1)',r,s,1)';
      lvec[tel,1]=afg1'*inv(-afg2)*afg1;
      lvec[tel,2]=cdfchic(lvec[tel,1],2);
     tel=tel+1;
   endo;

   format 10,4;
   volgorde=seqa(rows(chosen)+1,1,rows(pending));
   lvec=-lvec~volgorde;
   lvec=sortc(lvec,1);
   volgorde=lvec[.,3];
   lvec=-lvec[.,1:2];
   x=x[.,1:rows(chosen)]~x[.,volgorde];
   xh=xh[.,1:rows(chosen)]~xh[.,volgorde];
   var01=var01[1:rows(chosen)]|var01[volgorde];
   volgorde=volgorde-rows(chosen);
   pending=pending[volgorde];

   print " ";
   print "      variable      chisquare        p-value";
   print "      --------      ---------        -------";
   cnt=printfm(pending~lvec[.,1]~lvec[.,2],mask,fmt);
   print " ";
   if rows(pending)==1;
      print " ";
   endif;
   akkoord=(lvec[1,2]<=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]~0~lvec[1,.]);
      indv=indv|pending[1];
      if rows(pending)==1;
         stoppen=1;
      else;
         pending=pending[2:rows(pending)];
      endif;
   else;
      stoppen=1;
   endif;

endo;

chosen[rows(chosen),2]=miss(chosen[rows(chosen),2],0);

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[1 2 2 3],fmt[1 2 2 3,.]);
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 " ";

output off;



if _final;
   print "We now present the analysis for selected independent variables.";
   print " ";
   if _sfile==1;
       lprint "\f";
       _file=1;
   else;
      _file=0;
   endif;
   {xx,yy}=shrink(x[.,1:rows(chosen)],y);
   {schat0,w0,tetra,llt,teta1,ll1,teta2,ll2}=
   bprobrs(xx,yy,r,s,chosen[.,1],depvar);
else;
   output off;
endif;

retp(chosen[.,1],chosen[.,2:4]);
endp;



/*
** SHRINK
**
** FORMAT  {xx,yy}=shrink(x,y);
**
** INPUT   x  : Matrix of covariates.
**         y  : Matrix of contingency tables.
**
** OUTPUT  xx : Shrunk matrix of covariates.
**         yy : Shrunk matrix of contingencies.
**
** REMARK  SHRINK is only meant for internal use as auxiliary procedure.
*/




proc (2)= shrink(x,y);

local t,s,g1,g2,sortv,i,k,j,cnt,test;


/* Multiple columns sort on the independent variable(s). */

t=x~y;
g1=cols(x);
g2=g1+cols(y);
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);


/* Summing over equal independent variable vectors. */

i=1;
k=1;
s=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      s[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      s[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

s=s[1:k-1,.];


/*
** Separating x and y.
*/

y=s[.,g1+1:g2];
x=s[.,1:g1];


retp(x,y);
endp;





/*
** BPROBSC2
**
** FORMAT    {indv,info}=bprobsc2(dataset,indvar,depvar,sig);
**
** INPUT      dataset  : String, the name of the dataset that
**                       contains the variables.
**            indvar   : An Lx1 or (L-1)x1 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 of chosen variables.
**            info     : rows(indvar)x3 vector. Line j contains the log-
**                       likelihood with variables indv[1:j] included.
**                       info[j,2 3] contain the chisquared and P-value
**                       comparing the loglikelihood for variables 1 through
**                       j with the loglikelihood for variables 1 through
**                       j-1.
**
** 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     BPROBSC2 is a version of BPROBSC only supporting the case of
**            2x2 tables, thus increasing the speed.
*/


proc (2)= bprobsc2(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,t,k,
      tijd1,tijd2,teta2,pending,stoppen,llt,step,chosen,fmt2,mask2,aantal,s,j,
      indexi,indexd,nvec,hulp1,i,hulp2,g1,g2,hulp,sortv,l,afg1,afg2,volgorde,
      grad,hess,xx,yy,x,y,ll,var01,xh,yh;


tijd1=hsec;

/* Creating/resetting the output file. */

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


/* 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 :    Score test";
print "Significance level :" sig;
print "   Number of cases :" aantal;
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";




if dataset$=="";

   x=indvar;
   y=depvar;
   let depvar="y1" "y2";
   let indvar="constant" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;

else;

   {indvar,indexi}=indices(dataset,indvar);
   {depvar,indexd}=indices(dataset,depvar);
   nvec=setvars(dataset);
   open handle=^dataset for read;
   do until eof(handle);
      hulp1=readr(handle,100);
      i=1;
      do while i<=rows(indvar);
         hulp2=varget(indvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexi[i]];
         else;
            hulp2=hulp1[.,indexi[i]];
         endif;
         cnt=varput(hulp2,indvar[i]);
         i=i+1;
      endo;

      i=1;
      do while i<=rows(depvar);
         hulp2=varget(depvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexd[i]];
         else;
            hulp2=hulp1[.,indexd[i]];
         endif;
         cnt=varput(hulp2,depvar[i]);
         i=i+1;
      endo;

   endo;
   handle=close(handle);
   y=varget(depvar[1])~varget(depvar[2]);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;

   if not(x[.,1]==ones(rows(x),1));
      indvar="constant"|indvar;
   endif;

   nvec=setvars(dataset);

endif;



/* Checking for faulty data and deletion of rows with missing values. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=2);
   errorlog "y does not stand for a table of 0-1 variables.";
   end;
endif;

t=x~y;
_cases=rows(t);
g1=cols(x);
g2=cols(t);
t=packr(t);
_cases=_cases|rows(t);
y=t[.,g1+1:g2];
x=t[.,1:g1];


if sumc(sumc((y[.,1 2] .ne 1) .and (y[.,1 2] .ne 0)))/=0;
   errorlog "The dependent variables are not of the desired form.";
   end;
endif;


/* Transforming the dependent variables into contingency tables. */

hulp=zeros(rows(y),4);
hulp[.,1]=((y[.,1] .== 0) .and (y[.,2] .== 0));
hulp[.,2]=((y[.,1] .== 1) .and (y[.,2] .== 0));
hulp[.,3]=((y[.,1] .== 0) .and (y[.,2] .== 1));
hulp[.,4]=((y[.,1] .== 1) .and (y[.,2] .== 1));

y=hulp;
clear hulp;
g2=g2+2;


/* Multiple columns sort on the independent variable(s). */

t=x~y;
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);


/* Summing over equal independent variable vectors. */

i=1;
k=1;
s=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      s[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      s[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

s=s[1:k-1,.];


/*
** Separating dependent and independent variables and adding a column
** of ones to the matrix of independent variables if necessary.
*/

y=s[.,cols(s)-3:cols(s)];
x=s[.,1:cols(s)-4];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;


/* Construction of auxiliary matrices containing only the (0,1)-covariates. */

var01=zeros(cols(x),1);
i=1;
do while i<=cols(x);
   var01[i]=sumc((x[.,i].==1)+(x[.,i].==0))==rows(x);
   i=i+1;
endo;

{xh,yh}=shrink(x.*var01',y);


/* Intercept. */

output off;
_restric={1 1,0 0};
if var01[2]==0;
   {xx,yy}=shrink(x[.,1 2],y);
else;
   {xx,yy}=shrink(xh[.,1 2],yh);
endif;
{schat0,llt,tetra,llt,teta1,ll1,teta2,ll2}=
bivprob(xx,yy,indvar[1 2],depvar);
_restric=1;
output on;




/* Scanning of individual independent variables. */

tel=2;
lvec=zeros(rows(indvar)-1,2);
schat0=schat0[.,1];
l=(rows(schat0)-1)/2;
do while tel<=rows(indvar);
   schat1=schat0[1:l]|0|schat0[l+1:2*l]|0|schat0[2*l+1];
   if var01[tel]==0;
      {xx,yy}=shrink(x[.,1 tel],y);
   else;
      {xx,yy}=shrink(xh[.,1 tel],yh);
   endif;
   afg1=gradbp(schat1,xx,yy);
   afg2=hessbp(schat1,xx,yy);
   lvec[tel-1,1]=afg1'*inv(-afg2)*afg1;
   tel=tel+1;
endo;


tel=1;
do while tel<=rows(lvec);
   lvec[tel,2]=cdfchic(lvec[tel,1],2);
   tel=tel+1;
endo;


/* Ordering the independent variables. */

volgorde=seqa(2,1,rows(indvar)-1);
lvec=-lvec~volgorde;
lvec=sortc(lvec,1);
volgorde=lvec[.,3];
lvec=-lvec[.,1:2];



print "STEP  1";
print "=======";
print " ";
print " ";
mask=0~1~1;
let fmt[3,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4;
chosen="intercep"~llt~0~0;
print "Variables in model :";
print "--------------------";
format 10,4;
print $chosen[.,1]';
print " ";
print " ";
print "Candidates :";
print "------------";
print " ";
print "      variable      chisquare        p-value";
print "      --------      ---------        -------";
cnt=printfm(indvar[volgorde]~lvec[.,1]~(lvec[.,2]),mask,fmt);
print " ";
akkoord=(lvec[1,2]<=sig);
if akkoord;
   print $indvar[volgorde[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;
indvar=indvar[1]|indvar[volgorde];
x=x[.,1]~x[.,volgorde];
xh=xh[.,1]~xh[.,volgorde];
var01=var01[1]|var01[volgorde];
indv=indvar[1 2];
chosen=chosen|(indv[2]~0~lvec[1,.]);
if rows(indvar)>=3;
   pending=indvar[3:rows(indvar)];
   stoppen=not akkoord;
else;
   stoppen=1;
endif;
tel=1;
step=1;

do until stoppen;
   tel=1;
   step=step+1;
   output off;
   if var01[1:rows(chosen)]/=0;
      {xx,yy}=shrink(xh[.,1:rows(chosen)],yh);
   else;
      {xx,yy}=shrink(x[.,1:rows(chosen)],y);
   endif;
   {schat0,ll,tetra,llt,teta1,ll1,teta2,ll2}=
   bivprob(xx,yy,chosen[.,1],depvar);
   output on;
   chosen[rows(chosen),2]=ll;
   w0=ll;
   lvec=zeros(rows(pending),2);
   format 3,0;
   print "STEP" step;
   print "=======";
   print " ";
   print "Variables in model :";
   print "--------------------";
   format 10,4;
   print $chosen[.,1]';
   print " ";
   print "Candidates :";
   print "------------";
   l=(rows(schat0)-1)/2;
   do while tel<=rows(pending);
      schat1=schat0[1:l,1]|0|schat0[l+1:2*l,1]|0|schat0[2*l+1,1];
      if var01[1:rows(chosen) rows(chosen)+tel]/=0;
         {xx,yy}=shrink(xh[.,1:rows(chosen) rows(chosen)+tel],yh);
      else;
         {xx,yy}=shrink(x[.,1:rows(chosen) rows(chosen)+tel],y);
      endif;
      afg1=gradbp(schat1,xx,yy);
      afg2=hessbp(schat1,xx,yy);
      lvec[tel,1]=afg1'*inv(-afg2)*afg1;
      lvec[tel,2]=cdfchic(lvec[tel,1],2);
     tel=tel+1;
   endo;

   format 10,4;
   volgorde=seqa(rows(chosen)+1,1,rows(pending));
   lvec=-lvec~volgorde;
   lvec=sortc(lvec,1);
   volgorde=lvec[.,3];
   lvec=-lvec[.,1:2];
   x=x[.,1:rows(chosen)]~x[.,volgorde];
   xh=xh[.,1:rows(chosen)]~xh[.,volgorde];
   var01=var01[1:rows(chosen)]|var01[volgorde];
   volgorde=volgorde-rows(chosen);
   pending=pending[volgorde];

   print " ";
   print "      variable      chisquare        p-value";
   print "      --------      ---------        -------";
   cnt=printfm(pending~lvec[.,1]~lvec[.,2],mask,fmt);
   print " ";
   if rows(pending)==1;
      print " ";
   endif;
   akkoord=(lvec[1,2]<=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]~0~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[1 2 2 3],fmt[1 2 2 3,.]);
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 " ";

output off;



if _final;
   print "We now present the analysis for selected independent variables.";
   print " ";
   if _sfile==1;
       lprint "\f";
       _file=1;
   else;
      _file=0;
   endif;
   {xx,yy}=shrink(x[.,1:rows(chosen)],y);
   {schat0,w0,tetra,llt,teta1,ll1,teta2,ll2}=
   bivprob(xx,yy,chosen[.,1],depvar);
else;
   output off;
endif;

retp(chosen[.,1],chosen[.,2:4]);
endp;





