/*
**                         FILE BPMRS3.ARC
**
** This file contains the bivariate probit model procedures (rxs tables).
** Part 3.
*/


declare __output ?=0;




/*
** INITRS
**
** FORMAT     teta=initrs(x,y,r,s);
**
** INPUT      x         : NxL matrix of independent variables. The first column
**                        has to consist completely out of ones (to deal with
**                        the constant term).
**            y         : Nx(r*s) matrix of contingencies.
**
** OUTPUT     teta      : Estimated  L*(r+s-2) vector of estimates using the
**                        marginal distribution.
*/



proc initrs(x,y,r,s);

local n,a,b,i,j,teta1,teta2,loglik1,loglik2,teta3,a1,a2,b1,b2,loglik3,rho1,
   rho2,rho3,step,maxstep,y1,y2;


/* Working variables. */

n=rows(x);
y1=zeros(n,r);
y2=zeros(n,s);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      y1[.,i]=y1[.,i]+y[.,(i-1)*s+j];
      y2[.,j]=y2[.,j]+y[.,(i-1)*s+j];
      j=j+1;
   endo;
   i=i+1;
endo;


/* Initial coefficients. */

maxstep=3+(cols(x)-1)*2;
step=1;
a=zeros(cols(x),r-1);
b=zeros(cols(x),s-1);

if r>=3;
   a[1,.]=seqa(-10,20/(r-2),r-1)';
else;
   a[1,1]=0;
endif;

if s>=3;
   b[1,.]=seqa(-10,20/(s-2),s-1)';
else;
   b[1,1]=0;
endif;


/* Optimisation of cutpoints for first variable. */

format 0,0;
print "Step " step "of " maxstep;
format 10,4;
step=step+1;
loglik1=logluprs(vec(a),x,y1);
loglik2=logluprs(vec(a)/2,x,y1);
do while loglik1<loglik2;
   a=a/2;
   loglik1=loglik2;
   loglik2=logluprs(vec(a)/2,x,y1);
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "a";
   print a;
   print "loglik1 :" loglik1;
   print "loglik2 :" loglik2;*/
endo;


a1=a;
a2=a;
a1[1,.]=a1[1,.]-3;
a2[1,.]=a2[1,.]+3;
loglik1=logluprs(vec(a1),x,y1);
loglik2=logluprs(vec(a2),x,y1);
i=1;
do while i<=8;
   if loglik1<loglik2;
      a1=(a1+a2)/2;
      loglik1=logluprs(vec(a1),x,y1);
   else;
      a2=(a1+a2)/2;
      loglik2=logluprs(vec(a2),x,y1);
   endif;
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "a1";
   print a1;
   print "loglik1 :" loglik1;
   print " ";
   print "a2";
   print a2;
   print "loglik2 :" loglik2;*/
i=i+1;
endo;
if loglik1>loglik2;
   a=a1;
else;
   a=a2;
endif;



/* Optimisation of cutpoints for second variable. */

format 0,0;
print "Step " step "of " maxstep;
format 10,4;
step=step+1;
loglik1=logluprs(vec(b),x,y2);
loglik2=logluprs(vec(b)/2,x,y2);
do while loglik1<loglik2;
   loglik1=loglik2;
   b=b/2;
   loglik2=logluprs(vec(b)/2,x,y2);
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "b";
   print b;
   print "loglik1 :" loglik1;
   print "loglik2 :" loglik2;*/
endo;

b1=b;
b2=b;
b1[1,.]=b1[1,.]-3;
b2[1,.]=b2[1,.]+3;
loglik1=logluprs(vec(b1),x,y2);
loglik2=logluprs(vec(b2),x,y2);
i=1;
do while i<=8;
   if loglik1<loglik2;
      b1=(b1+b2)/2;
      loglik1=logluprs(vec(b1),x,y2);
   else;
      b2=(b1+b2)/2;
      loglik2=logluprs(vec(b2),x,y2);
   endif;
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "b1";
   print b1;
   print "loglik1 :" loglik1;
   print " ";
   print "b2";
   print b2;
   print "loglik2 :" loglik2;*/
i=i+1;
endo;
if loglik1>loglik2;
   b=b1;
else;
   b=b2;
endif;


/* Optimisation of covariate-parameters for first variable. */

i=2;
do while i<=cols(x);
   format 0,0;
   print "Step " step "of " maxstep;
   format 10,4;
   step=step+1;
   a1=a;
   a2=a;
   a1[i,.]=-5*ones(1,cols(a1));
   a2[i,.]=5*ones(1,cols(a2));
   j=1;
   do while j<=8;
      loglik1=logluprs(vec(a1),x,y1);
      loglik2=logluprs(vec(a2),x,y1);
      if loglik1>loglik2;
         a2=(a1+a2)/2;
         loglik2=logluprs(vec(a2),x,y1);
      else;
         a1=(a1+a2)/2;
         loglik1=logluprs(vec(a1),x,y1);
      endif;
      j=j+1;
      /*cls;
      print "Current initial estimates :";
      print "---------------------------";
      print "a1";
      print a1;
      print "loglik1 :" loglik1;
      print " ";
      print "a2";
      print a2;
      print "loglik2 :" loglik2;*/
   endo;
i=i+1;
if loglik1<loglik2;
   a=a2;
else;
   a=a1;
endif;
endo;


/* Optimisation of covariate-parameters for second variable. */

i=2;
do while i<=cols(x);
   format 0,0;
   print "Step " step "of " maxstep;
   format 10,4;
   step=step+1;
   b1=b;
   b2=b;
   b1[i,.]=-5*ones(1,cols(b1));
   b2[i,.]=5*ones(1,cols(b2));
   j=1;
   do while j<=8;
      loglik1=logluprs(vec(b1),x,y2);
      loglik2=logluprs(vec(b2),x,y2);
      if loglik1>loglik2;
          b2=(b1+b2)/2;
         loglik2=logluprs(vec(b2),x,y2);
      else;
           b1=(b1+b2)/2;
         loglik1=logluprs(vec(b1),x,y2);
      endif;
      j=j+1;
      /*cls;
      print "Current initial estimates :";
      print "---------------------------";
      print "b1";
      print b1;
      print "loglik1 :" loglik1;
      print " ";
      print "b2";
      print b2;
      print "loglik2 :" loglik2;*/
   endo;
i=i+1;
if loglik1<loglik2;
   b=b2;
else;
   b=b1;
endif;
endo;



/* Optimisation of correlation coefficient. */

format 0,0;
print "Step " step "of " maxstep;
format 10,4;
step=step+1;
teta1=vec(a)|vec(b);
rho1=-0.9;
rho2=0.9;
loglik1=loglikrs(teta1|rho1,x,y,r,s);
loglik2=loglikrs(teta1|rho2,x,y,r,s);
i=1;
do while i<=7;
   rho3=(rho1+rho2)/2;
   loglik3=loglikrs(teta1|rho3,x,y,r,s);
   if loglik1>loglik2;
      rho2=rho3;
      loglik2=loglik3;
   else;
      rho1=rho3;
      loglik1=loglik3;
   endif;
   i=i+1;
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "teta1";
   print teta1'~rho1;
   print "loglik1 :" loglik1;
   print " ";
   print "teta2";
   print teta1'~rho2;
   print "loglik2 :" loglik2;*/
endo;

if loglik1>loglik2;
   teta1=teta1|rho1;
else;
   teta1=teta1|rho2;
endif;


retp(teta1);
endp;

