/*
**                         FILE BPM9.ARC
**
** This file contains the bivariate probit model procedures.
** Part 9.
*/



/* Declaration of globals. */

declare _draw !=1;
declare _initval !=0;


/*
** PRBLIK
**
** FORMAT   {t,rho,par,pl,tvec,pvec}=prblik(dataset,indvar,depvar,sig);
**
** INPUT    indvar   : Lx1 or (L-1)x1 character vector; names of independent
**                     variables.
**          depvar   : 2x1 character vector; names of dependent variables.
**          sig      : Significance level for the chi-square test (e.g.
**                     0.95 or 0.98).
**
** OUTPUT   t        : 3x3 vector of rho values at level sig and the maximum
**                     together with likelihood and bound
**          rho      : Vector of rho values between bounds.
**          par      : Kx4*L matrix of parameters and standard errors.
**          pl       : Vector of log-likelihoods.
**          tvec     : Matrix of t-values.
**          pvec     : Matrix of p-values.
**
*/


proc (6)= prblik(dataset,indvar,depvar,sig);

local x,y,teta1,waarde1,eerste1,tweede1,test,test1,test2,tetra,llt;
local i,j,k,g1,g2,cnt,s,test,t,hulp,sortv,rhomax,teta1,ll1,teta2,ll2;
local indexi,indexd,nvec,handle,hulp1,i,hulp2,cnt,choice;
local rho,par,pl,tvec,pvec,step;
local aantal,onder,boven,min,max,g1,g2,g3,g4,min1,min2,max1,max2,hoog;
local ll1,ll2,ll3,ll4,llmax;

disable;

step=0.1;

/*
** 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;
g1=cols(x);
g2=cols(t);
t=packr(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;


/* Calling  PROFLIK  */

print "DETERMINING THE MAXIMUM";
screen off;
{rhomax,llmax,tetra,llt,teta1,ll1,teta2,ll2}=bivprob(x,y,indvar,depvar);
screen on;
clear eerste;
clear tweede;
rhomax=rhomax[rows(rhomax),1];

{rho,par,pl,tvec,pvec}=proflik(x,y,indvar,depvar,step,-1,1);
format /rd 10,4;

test=2*(llmax-pl);
test=1-cdfchic(test,1);
g2=1;
do while test[g2]>sig;
   g2=g2+1;
endo;
g3=rows(pl);
do while test[g3]>sig;
   g3=g3-1;
endo;
g1=g2-1;
g4=g3+1;

min1=rho[g1];
max1=rho[g2];
min2=rho[g3];
max2=rho[g4];


/* Lower bound. */

cls;
print "Bounds/Tests";
print (min1~max1~min2~max2)|(test[g1]~test[g2]~test[g3]~test[g4]);
print "**********************************************************************";
aantal=1;
do while aantal<=3;
   step=(max1-min1)/20;
   {rho,par,pl,tvec,pvec}=proflik(x,y,indvar,depvar,step,min1,max1);
   format /rd 10,4;
   test1=2*(llmax-pl);
   test1=1-cdfchic(test1,1);
   onder=1;
   do while test1[onder]>sig;
      onder=onder+1;
   endo;
   boven=onder;
   onder=onder-1;
   if onder<1; onder=1; endif;
   if boven>rows(rho); boven=rows(rho);endif;
   min1=rho[onder];
   max1=rho[boven];
   locate 1,1;
   print "Bounds/Tests";
   print (min1~max1~min2~max2)|(test1[onder]~test1[boven]~test[g3]~test[g4]);
   print "*******************************************************************"\
"***";
   aantal=aantal+1;
endo;

g1=onder;
g2=boven;
ll1=pl[g1];
ll2=pl[g2];

/* Upper bound. */

aantal=1;
do while aantal<=3;
   step=(max2-min2)/20;
   {rho,par,pl,tvec,pvec}=proflik(x,y,indvar,depvar,step,min2,max2);
   format /rd 10,4;
   test2=2*(llmax-pl);
   test2=1-cdfchic(test2,1);
   onder=1;
   do while test2[onder]<sig;
      onder=onder+1;
   endo;
   boven=onder;
   onder=onder-1;
   if onder<1; onder=1; endif;
   if boven>rows(rho); boven=rows(rho);endif;
   min2=rho[onder];
   max2=rho[boven];
   locate 1,1;
   print "Bounds/Tests";
   print (min1~max1~min2~max2)|(test1[g1]~test1[g2]~test2[onder]~test2[boven]);
   print "*******************************************************************"\
"***";
   aantal=aantal+1;
endo;

g3=onder;
g4=boven;
ll3=pl[g3];
ll4=pl[g4];


/* Choice of returns for bounds. */

t=min1~ll1~test1[g1]|max1~ll2~test1[g2];
t=t|min2~ll3~test2[g3]|max2~ll4~test2[g4];

choice=abs(sig-t[.,3]);
if choice[1]<choice[2];
   choice[1]=1;
else;
   choice[1]=2;
endif;
if choice[3]<choice[4];
   choice[2]=3;
else;
   choice[2]=4;
endif;
choice=choice[1 2];

t[3,.]=t[choice[2],.];
t[1,.]=t[choice[1],.];
t[2,.]=rhomax~llmax~0;

t=t[1:3,.];


/* Profile likelihood between bounds. */

print "Final Profile Likelihood Computation.                                  ";
print "-------------------------------------                                  ";
step=(t[3,1]-t[1,1])/20;
{rho,par,pl,tvec,pvec}=proflik(x,y,indvar,depvar,step,t[1,1],t[3,1]);


retp(t,rho,par,pl,tvec,pvec);

endp;


