/*
**                         FILE BPM6.ARC
**
** This file contains the bivariate probit model procedures.
** Part 6.
*/



/* Declaration of globals. */

declare _rpar ?=0;
declare _cases ?=0;



/*
** BPROBIT
**
** FORMAT     {teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=
**                           bprobit(dataset,indvar,depvar);
**
**            or
**
**            {rho,par,pl,tvec,pvec,stepsize,min,max}=
**                           bprobit(dataset,indvar,depvar);
**
** INPUT      dataset  : Either a string, the name of the dataset that
**                       contains the variables or the null string.
**            indvar   : Either an Lx1 or (L-1)x1 character vector containing
**                       the names of the independent variables or,
**                       in case dataset is the null string, an
**                       NxL or Nx(L-1) matrix of independent variables.
**                       (L-1 if the first column isn't a unit column; in that
**                       case a unit column is appended and L-1 is incremented).
**            depvar   : Either a 2x1 character vector containing the names
**                       of the dependent variables or, in case dataset is the
**                       null string, an Nx2 matrix of dependent 0-1 variables.
**
** OUTPUT     If _rpar=0 or 1 :
**
**            teta     : (2*L+1)x4 or (3*L)*4 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).
**
**            If _rpar=2 :
**
**            rho      : Kx1 vector of rho-values used.
**            par      : (4*L)xK vector of parameters and standard errors.
**            pl       : Kx1 vector of log-likelihoods at K rho values.
**            tvec     : Kx1 vector of t-values at K rho values.
**            pvec     : Kx1 vector of p-values at K rho values.
**            stepsize : The interval between rho-values.
**            min      : Minimum value for rho.
**            max      : Maximum value for rho.
**
** GLOBALS    _rpar:  0: The correlation coefficient is a constant to be
**                       estimated (default).
**                    1: The correlation depends on the independent variables
**                       in the same way as the variables do.
**                    2: A profile likelihood picture is designed. In this case
**                       the inputs are the same but the outputs are different.
**                       The user will be requested upon stepsize for rho.
**
** REMARK     BPROBIT transforms input to a form, tractible by BIVPROB
**            and PROFLIK.  Depending on _rpar, one of these
**            procedures is called. Each of them can be called directly
**            (see specifications for these procedures), but using BPROBIT
**            first is more flexible and will save time in general.
*/



proc (8)= bprobit(dataset,indvar,depvar);

local x,y,teta,ll,teta1,ll1,teta2,ll2,tetra,llt,max,min,i,j,k,g1,g2,cnt,s,
      indexi,indexd,nvec,handle,hulp1,i,hulp2,cnt,test,t,hulp,sortv;


disable;

/*
** If the input is on a dataset, it is processed in the subsequent, if not
** it is read from the matrices passed in indvar and depvar.
*/

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;

/*xg=x;yg=y;*/


/* Calling BIVPROB (_rpar=0 or 1) or PROFLIK (_rpar=2). */

if _rpar==0 or _rpar==1;
   {teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=bivprob(x,y,indvar,depvar);
   retp(teta,ll,tetra,llt,teta1,ll1,teta2,ll2);
elseif _rpar==2;
   local stepsize,rho,par,pl,tvec,pvec;
   print "What step size do you want for rho ?";
   stepsize=con(1,1);
   min=-1+stepsize/2;
   max=-min;
   {rho,par,pl,tvec,pvec}=proflik(x,y,indvar,depvar,stepsize,min,max);
   retp(rho,par,pl,tvec,pvec,stepsize,min,max);
   end;
endif;


endp;


