/*
**                         FILE BPM16.ARC
**
** This file contains the bivariate probit model procedures.
** Part 16.
*/



/* Declaration of globals. */

declare _rpar ?=0;
declare _cases ?=0;
declare _file ?=0;
declare _eps ?=0.001;
declare _maxit ?=30;
declare xg ?=0;
declare yg ?=0;




/*
** MCORR
**
** FORMAT       {rho,eenr,tweer}=mcorr(teta0,x,y);
**
** INPUT        teta0 :   (2*L)x1 vector of cutpoints.
**                        The cutpoints are considered to be exact.
**                  x :   Matrix of covariates.
**                  y :   Contingency table, ordered as a Nx4 matrix.
**
** OUTPUT       rho   :   Scalar, optimisation for the correlation coefficient.
**              eenr  :   First order derivative.
**              tweer :   Second order derivative.
**
** GLOBALS      _eps  :   Criterion to end the Newton-Raphson iterations.
**              _maxit:   Upper bound for the number of iterations.
**
*/



proc (3)= mcorr(teta0,x,y);

local p,teller,eenr,s12,pp,rho,vrho,afgr,vrr,tweer,aa,bb,a,b,l,n;


/* Initialising some values. */

l=rows(teta0)/2;
a=teta0[1:l];
b=teta0[l+1:2*l];
aa=x*a;
bb=x*b;
rho=0;
teller=1;
eenr=1;
s12=1~-1~-1~1;
n=rows(x);


/* Newton-Raphson iterations. */

do while (teller<=_maxit) and (sumc(abs(eenr)) > _eps);

   p=zeros(n,4);
   p[.,1]=cdfbvn(-aa,-bb,rho);
   p[.,2]=cdfn(-bb)-p[.,1];
   p[.,3]=cdfn(-aa)-p[.,1];
   p[.,4]=1-p[.,1]-p[.,2]-p[.,3];
   p=p.*(p .> 1e-12) + (1e-12).*(p .<=1e-12);

   pp=p^2;

   vrho=zeros(n,4);
   vrho[.,1]=y[.,1].*eerster(-aa,-bb,rho);
   vrho[.,2]=y[.,2].*eerster(aa,-bb,-rho);
   vrho[.,3]=y[.,3].*eerster(-aa,bb,-rho);
   vrho[.,4]=y[.,4].*eerster(aa,bb,rho);

   vrho=vrho./p;
   eenr=sumc(sumc(vrho.*s12));

   afgr=zeros(n,4);
   afgr[.,1]=eerster(-aa,-bb,rho);
   afgr[.,2]=eerster(aa,-bb,-rho);
   afgr[.,3]=eerster(-aa,bb,-rho);
   afgr[.,4]=eerster(aa,bb,rho);

   vrr=zeros(n,4);
   vrr[.,1]=y[.,1].*(p[.,1].*tweederr(-aa,-bb,rho)-(afgr[.,1])^2);
   vrr[.,2]=y[.,2].*(p[.,2].*tweederr(aa,-bb,-rho)-(afgr[.,2])^2);
   vrr[.,3]=y[.,3].*(p[.,3].*tweederr(-aa,bb,-rho)-(afgr[.,3])^2);
   vrr[.,4]=y[.,4].*(p[.,4].*tweederr(aa,bb,rho) - (afgr[.,4])^2);
   vrr=vrr./pp;
   tweer=sumc(sumc(vrr));
   rho=rho-inv(tweer)*eenr;
   if rho>=1;
      rho=1-1e-4;
   elseif rho<=-1;
      rho=-1+1e-4;
   endif;
   cls;
   format /rd 5,0;
   print "CORRELATION  COEFFICIENT :  ITERATION NUMBER" teller;
   print "----------------------------------------------";
   print " ";
   format /rd 14,-log(_eps)+2;
   print "Parameter/gradient at current iteration :";
   print rho|eenr;
   print " ";
   print "*****************************************************************"\
   "***************";
   print " ";
   teller=teller+1;

endo;

retp(rho,eenr,tweer);
endp;







/*
** MPROBIT
**
** FORMAT     {teta,corr}=mprobit(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 Mx1 character vector containing the names
**                       of the dependent variables or, in case dataset is the
**                       null string, an NxM matrix of dependent 0-1 variables.
**
** OUTPUT     teta     : Matrix of estimated parameters for covariates, the
**                       meanings of the columns: estimated parameter,
**                       standard error, t-value and p-value.
**            corr     : Matrix of estimated correlations for the pairs
**                       (i,j) of dependent variables.
**
** 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.
**
*/



proc (2)= mprobit(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,i1,i2,
      yy,xx,teta0,corr,tel,grad,hess,pval,tval,stde,rho;



disable;

if _file==0;
   output file=mprobit.out reset;
elseif _file==1;
   output file="lpt1";
endif;
output off;


/*
** 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" "y3" "y4" "y5" "y6" "y7" "y8" "y9" "y10";
   let indvar="constant" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9" "x10";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;
   depvar=depvar[1:cols(y)];

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);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;
   y=varget(depvar[1]);
   i=2;
   do while i<=rows(depvar);
      y=y~varget(depvar[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;
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 .ne 1) .and (y .ne 0)))/=0;
   errorlog "The dependent variables are not of the desired form.";
   end;
endif;


/* To compute loglikelihood in 3-dimensional case. */


if cols(y)==3;
   {xg,yg}=conting3(x,y);
endif;

{xg,yg}=conting4(x,y);stop;


/* Initializing some values. */

corr=zeros(1,4);
if not(x[.,1]==1);
   teta0=zeros(1,4);
else;
   teta0=zeros(1,4);
endif;


/* CALLING THE UNIVARIATE PROCEDURES. */

i1=1;
do while i1<=cols(y);
   xx=x;
   yy=y[.,i1];


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

   hulp=zeros(rows(yy),2);
   hulp[.,1]=(yy .== 0);
   hulp[.,2]=(yy .== 1);

   yy=hulp;
   clear hulp;


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

   t=xx~yy;
   sortv=seqa(1,1,g1);
   t=sortmc(t,sortv);


   /* Summing over equal independent variable vectors. */

   i=1;
   k=1;
   s=zeros(rows(t),cols(t));
   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:cols(s)]=sumc(t[i:j,g1+1:cols(t)])';
         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.
   */

   yy=s[.,cols(s)-1:cols(s)];
   xx=s[.,1:cols(s)-2];
   if not(xx[.,1]==1);
      xx=ones(rows(xx),1)~xx;
   endif;



   {teta1,ll1,grad,hess}=uniprob(zeros(cols(xx),1),xx,yy);
   stde=zeros(cols(xx),1);
   tval=zeros(cols(xx),1);
   pval=zeros(cols(xx),1);
   if ll1/=0;
      stde=sqrt(diag(inv(-hess)));
      tval=teta1./stde;
      pval=2*(1-cdfn(abs(tval)));
   endif;

   teta0=teta0|(teta1~stde~tval~pval);


   i1=i1+1;
endo;

teta0=teta0[2:rows(teta0),.];
cls;


/*
** CALLING THE BIVARIATE PROCEDURES FOR ALL COMBINATIONS OF DEPENDENT
** VARIABLES.
*/

g2=g1+4;
tel=1;
i1=1;
do while i1<=cols(y)-1;
   i2=i1+1;
   do while i2<=cols(y);
      format 0,0;
      print "Determining correlation ( " i1 ", " i2 ")";
      format 10,4;
      xx=x;
      yy=y[.,i1 i2];

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

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

      yy=hulp;
      clear hulp;


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

      t=xx~yy;
      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.
      */

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




      /* Calling MCORR */

      teta1=teta0[cols(xx)*(i1-1)+1:cols(xx)*i1,1];
      teta2=teta0[cols(xx)*(i2-1)+1:cols(xx)*i2,1];
      {rho,grad,hess}=mcorr(teta1|teta2,xx,yy);
      stde=sqrt(diag(inv(-hess)));
      tval=rho./stde;
      pval=2*(1-cdfn(abs(tval)));
      corr=corr|(rho~stde~tval~pval);


      i2=i2+1;
   endo;
   i1=i1+1;
endo;

corr=corr[2:rows(corr),.];


/* Output routines. */

if _file/=2;
   output on;
endif;

cls;
print " ";
print "            MULTIVARIATE  PROBIT  ESTIMATION  OF  PARAMETERS";
print "            ================================================";
print " ";
print " ";

print "Dependent variables :";
print "---------------------";
print $depvar';
print " ";
print " ";
print "Independent variables :";
print "-----------------------";
print $indvar';
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";

print "U N I V A R I A T E      P R O B I T     P A R A M E T E R S";
print "============================================================";
print " ";
print " ";


print "  Dep.Var.   Ind.Var.  |    Coeff    Std.err.    t-value    p-value";
print "  -----------------------------------------------------------------";
format 10,4;
tel=1;
i=1;
do while i<=rows(depvar);
   j=1;
   do while j<=rows(indvar);
      if j==1;
         print $depvar[i];;
      else;
         print "           ";;
      endif;
      print $indvar[j];;
      print " |";;
      print teta0[tel,.];
      tel=tel+1;
      j=j+1;
   endo;
   print " ";
   i=i+1;
endo;


print " ";
print " ";
print "B I V A R I A T E      P R O B I T     C O R R E L A T I O N S";
print "==============================================================";
print " ";
print " ";


print "  Dep.Var1   Dep.Var2  |    Coeff    Std.err.    t-value    p-value";
print "  -----------------------------------------------------------------";
format 10,4;
tel=1;
i=1;
do while i<=rows(depvar)-1;
   j=i+1;
   do while j<=rows(depvar);
      print $depvar[i j]';;
      print " |";;
      print corr[tel,.];
      tel=tel+1;
      j=j+1;
   endo;
   i=i+1;
endo;

if _file /=2;
   output off;
endif;

if _file==1;
   lprint "\f";
endif;

retp(teta0,corr);
endp;



/*
** CONTING3
**
** FORMAT   {xc,yc}=conting3(x,y);
**
** INPUT    x :  NxL matrix, representing L covariates.
**          y :  Nx3 matrix, representing 3 dichotomous variables.
**
** OUTPUT   xc:  MxL matrix, representing L covariates.
**          yc:  Mx8 matrix, representing the 3 dichotomous variables,
**               in the form of contingency tables.
**
*/



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

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


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

y=hulp;


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

t=x~y;
g1=cols(x);
g2=cols(x)+8;
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)-7:cols(s)];
x=s[.,1:cols(s)-8];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;



retp(x,y);
endp;




/*
** CONTING4
**
** FORMAT   {xc,yc}=conting3(x,y);
**
** INPUT    x :  NxL matrix, representing L covariates.
**          y :  Nx3 matrix, representing 4 dichotomous variables.
**
** OUTPUT   xc:  MxL matrix, representing L covariates.
**          yc:  Mx16 matrix, representing the 4 dichotomous variables,
**               in the form of contingency tables.
**
*/



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

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


hulp=zeros(rows(y),16);
hulp[.,1]=((y[.,1] .== 1) .and (y[.,2] .== 1) .and (y[.,3] .== 1) .and
           (y[.,4] .== 1));
hulp[.,2]=((y[.,1] .== 0) .and (y[.,2] .== 1) .and (y[.,3] .== 1) .and
           (y[.,4] .== 1));
hulp[.,3]=((y[.,1] .== 1) .and (y[.,2] .== 0) .and (y[.,3] .== 1) .and
           (y[.,4] .== 1));
hulp[.,4]=((y[.,1] .== 0) .and (y[.,2] .== 0) .and (y[.,3] .== 1) .and
           (y[.,4] .== 1));
hulp[.,5]=((y[.,1] .== 1) .and (y[.,2] .== 1) .and (y[.,3] .== 0) .and
           (y[.,4] .== 1));
hulp[.,6]=((y[.,1] .== 0) .and (y[.,2] .== 1) .and (y[.,3] .== 0) .and
           (y[.,4] .== 1));
hulp[.,7]=((y[.,1] .== 1) .and (y[.,2] .== 0) .and (y[.,3] .== 0) .and
           (y[.,4] .== 1));
hulp[.,8]=((y[.,1] .== 0) .and (y[.,2] .== 0) .and (y[.,3] .== 0) .and
           (y[.,4] .== 1));

hulp[.,9]=((y[.,1] .== 1) .and (y[.,2] .== 1) .and (y[.,3] .== 1) .and
           (y[.,4] .== 0));
hulp[.,10]=((y[.,1] .== 0) .and (y[.,2] .== 1) .and (y[.,3] .== 1) .and
           (y[.,4] .== 0));
hulp[.,11]=((y[.,1] .== 1) .and (y[.,2] .== 0) .and (y[.,3] .== 1) .and
           (y[.,4] .== 0));
hulp[.,12]=((y[.,1] .== 0) .and (y[.,2] .== 0) .and (y[.,3] .== 1) .and
           (y[.,4] .== 0));
hulp[.,13]=((y[.,1] .== 1) .and (y[.,2] .== 1) .and (y[.,3] .== 0) .and
           (y[.,4] .== 0));
hulp[.,14]=((y[.,1] .== 0) .and (y[.,2] .== 1) .and (y[.,3] .== 0) .and
           (y[.,4] .== 0));
hulp[.,15]=((y[.,1] .== 1) .and (y[.,2] .== 0) .and (y[.,3] .== 0) .and
           (y[.,4] .== 0));
hulp[.,16]=((y[.,1] .== 0) .and (y[.,2] .== 0) .and (y[.,3] .== 0) .and
           (y[.,4] .== 0));

y=hulp;


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

t=x~y;
g1=cols(x);
g2=cols(x)+16;
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)-15:cols(s)];
x=s[.,1:cols(s)-16];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;



retp(x,y);
endp;



