/*
**                         FILE BPM3.ARC
**
** This file contains the bivariate probit model procedures.
** Part 3.
*/





/*
** INITVAL
**
** FORMAT     {teta,loglik1,alt,loglik2}=initval(x,y);
**
** INPUT      x         : NxL matrix of independent variables. The first column
**                        has to consist completely out of ones (to deal with
**                        the constant term).
**            y         : Nx4 matrix of contingencies.
**
** OUTPUT     teta      : Estimated initial value to start the Newton-Raphson
**                        iteration.
**            loglik1   : Log-likelihood of initial values.
**            alt       : Alternative initial estimates.
**            loglik2   : Log-likelihood of alternative initial estimates.
**
** REMARK     The initial estimate for the coefficients of the variables
**            are determined by OLS through the probabilities of the
**            observed frequencies. The initial estimate for the correlation
**            is taken to be the maximum of the log-likelihood function with
**            the initial estimates of the coefficients already determined.
**            If this last step fails, initial estimates are
**            determined alternatively by fitting correlations to the
**            observed frequencies.
*/


proc (4)= initval(x,y);

local prop1,prop2,i,n,aa,bb,prop,test,teta;
local vnam,m,a,stb,rc,stderr,sigma,cx,rsq,resid,dwstat,b;
local rho,rhoa,rhob,rhoc,rhod,rhogem,rho1,rho2,rhohulp;
local ll,logl1,logl,logla,loglb,loglc,logld;
local rhovec,loglvec,propvec,naamvec,fmt,mask,selectie,cnt;
local teller,rho3,s12,p,vrho,afgr,eenr,pp,vrr,tweer;
local loglik1,loglik2,teta1,teta2,alt,teta3,spring;
clearg __output;


output off;


/* Initial value for a and b by means of OLS */

print "Determining starting point for initial values of a and b.";
n=rows(y);
prop1=(y[.,2]+y[.,4])./sumc(y');
prop2=(y[.,3]+y[.,4])./sumc(y');

prop1=cdfninv(prop1);
prop2=cdfninv(prop2);

if rows(x)>cols(x);
__output=0;
{vnam,m,a,stb,rc,stderr,sigma,cx,rsq,resid,dwstat}=ols("",prop1,x[.,2:cols(x)]);
{vnam,m,b,stb,rc,stderr,sigma,cx,rsq,resid,dwstat}=ols("",prop2,x[.,2:cols(x)]);
else;
   a=inv(x)*prop1;
   b=inv(x)*prop2;
endif;


/* Optimization of initial values. */

print "Optimisation of initial values of a and b.";
teta1=a|b|0;
loglik1=loglik(teta1,x,y);
teta2=teta1/2;
loglik2=loglik(teta2,x,y);
/*cls;
print "Current initial estimates :";
print "---------------------------";
print "teta1";
print teta1';
print "loglik1 :" loglik1;
print " ";
print "teta2";
print teta2';
print "loglik2 :" loglik2;*/
do while loglik1<loglik2;
   teta1=teta2;
   loglik1=loglik2;
   teta2=teta1/2;
   loglik2=loglik(teta2,x,y);
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "teta1";
   print teta1';
   print "loglik1 :" loglik1;
   print " ";
   print "teta2";
   print teta2';
   print "loglik2 :" loglik2;*/
endo;

print "Further optimisation of initial values for a and b.";
loglik1=loglik(teta1,x,y);
loglik2=loglik(teta2,x,y);
/*cls;
print "Current initial estimates :";
print "---------------------------";
print "teta1";
print teta1';
print "loglik1 :" loglik1;
print " ";
print "teta2";
print teta2';
print "loglik2 :" loglik2;*/

do while abs(loglik1-loglik2)>5;
   teta3=(teta1+teta2)/2;
   if loglik1>loglik2;
      teta2=teta3;
      loglik2=loglik(teta2,x,y);
   else;
      teta1=teta3;
      loglik1=loglik(teta1,x,y);
   endif;
   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "teta1";
   print teta1';
   print "loglik1 :" loglik1;
   print " ";
   print "teta2";
   print teta2';
   print "loglik2 :" loglik2;*/
endo;
if loglik2>loglik1;
   teta1=teta2;
endif;

a=teta1[1:rows(a)];
b=teta1[rows(b)+1:2*rows(b)];




/*
** Initial value for the correlation coefficient rho :
** Maximalisation of the log-likelihood function with a and b fixed and
** rho=0 as starting value.
*/

print "Determining an initial value for the correlation coefficient.";

aa=x*a;
bb=x*b;

rho=0;
eenr=1;
teller=1;
spring=0;
do while (teller<=20) and (sumc(abs(eenr)) > 1e-6);

   s12=ones(n,4);
   s12[.,2:3]=-s12[.,2:3];

   p=zeros(n,4);
   p[.,1]=cdfbvn(-aa,-bb,rho);
   p[.,2]=cdfn(-bb)-p[.,1];
   p[.,3]=cdfn(-aa)-p[.,1];
   p[.,4]=ones(n,1)-p[.,1]-p[.,2]-p[.,3];

   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;

   /*cls;
   print "Current initial estimates :";
   print "---------------------------";
   print "teta";
   print (a|b|rho)';
   print "loglik :" loglik(a|b|rho,x,y);*/

   teller=teller+1;

   if abs(rho)<1;
      spring=0;
   elseif ((rho>=1) and (spring==0));
      rho=0.9;
      spring=1;
      print "Upper truncation needed.";
   elseif ((rho<=-1) and (spring==0));
      rho=-0.9;
      spring=1;
      print "Lower truncation needed.";
   elseif spring==1;
      print "We will proceed in an alternative way.";
      goto label;
   endif;

endo;


/*
** If the estimate for rho as just determined is not satisfying, we choose
** an initial estimate for rho in another way, i.e. by choosing, for every
** entry in y, the rho that best fits to the observed frequency, and then
** taking the weighted average over every group.
*/


if (eenr>1e-6);

label:

print "Alternative determination of an initial estimate for the correlation"\
"coefficient";

   /* contingency 1 */

   prop=y[.,1]./sumc(y');
   rhoa=zeros(n,1);
   test=cdfbvn(-aa,-bb,rhoa);
   rho2=(test .lt prop);
   rho1= -(test .ge prop);

   i=1;
   do while (i<=10);
      rhohulp=(rho1+rho2)/2;
      test=cdfbvn(-aa,-bb,rhohulp);
      rho2=rho2.*(test .le prop) + rhohulp.*(test .gt prop);
      rho1=rho1.*(test .gt prop) + rhohulp.*(test .le prop);
      i=i+1;
   endo;

   rhoa=(meanc(rho1)+meanc(rho2))/2;


   /* contingency 2 */

   prop=y[.,2]./sumc(y');
   rhob=zeros(n,1);
   test=cdfbvn(aa,-bb,-rhob);
   rho2=(test .ge prop);
   rho1= -(test .lt prop);

   i=1;
   do while (i<=10);
      rhohulp=(rho1+rho2)/2;
      test=cdfbvn(aa,-bb,-rhohulp);
      rho2=rho2.*(test .gt prop) + rhohulp.*(test .le prop);
      rho1=rho1.*(test .le prop) + rhohulp.*(test .gt prop);
      i=i+1;
   endo;

   rhob=(meanc(rho1)+meanc(rho2))/2;


   /* contingency 3 */

   prop=y[.,3]./sumc(y');
   rhoc=zeros(n,1);
   test=cdfbvn(-aa,bb,-rhoc);
   rho2=(test .ge prop);
   rho1= -(test .lt prop);

   i=1;
   do while (i<=10);
      rhohulp=(rho1+rho2)/2;
      test=cdfbvn(-aa,bb,-rhohulp);
      rho2=rho2.*(test .gt prop) + rhohulp.*(test .le prop);
      rho1=rho1.*(test .le prop) + rhohulp.*(test .gt prop);
      i=i+1;
   endo;

   rhoc=(meanc(rho1)+meanc(rho2))/2;


   /* contingency 4 */

   prop=y[.,4]./sumc(y');
   rhod=zeros(n,1);
   test=cdfbvn(aa,bb,rhod);
   rho2=(test .lt prop);
   rho1= -(test .ge prop);

   i=1;
   do while (i<=10);
      rhohulp=(rho1+rho2)/2;
      test=cdfbvn(aa,bb,rhohulp);
      rho2=rho2.*(test .le prop) + rhohulp.*(test .gt prop);
      rho1=rho1.*(test .gt prop) + rhohulp.*(test .le prop);
      i=i+1;
   endo;

   rhod=(meanc(rho1)+meanc(rho2))/2;


   rhogem=((rhoa~rhob~rhoc~rhod)*sumc(y))/sumc(sumc(y));


   /* Log-likelihood for the different rho-values. */

   logl=loglik(a|b|rhogem,x,y);
   logla=loglik(a|b|rhoa,x,y);
   loglb=loglik(a|b|rhob,x,y);
   loglc=loglik(a|b|rhoc,x,y);
   logld=loglik(a|b|rhod,x,y);


   /* Values to place in a table. */

   rhovec=rhoa|rhob|rhoc|rhod|rhogem;
   loglvec=logla|loglb|loglc|logld|logl;
   propvec=(sumc(y)/sumc(sumc(y)))|1;
   let naamvec= "rhoa  " "rhob  " "rhoc  " "rhod  " "rhogem";


   /* Selection of rho with best likelihood. */

   prop=1;
   rho=rhogem;
   logl1=logl;
   if logl1<logla;
     rho=rhoa;
     logl1=logla;
     prop=propvec[1];
   endif;
   if logl1<loglb;
     rho=rhob;
     logl1=loglb;
     prop=propvec[2];
   endif;
   if logl1<loglc;
     rho=rhoc;
     logl1=loglc;
     prop=propvec[3];
   endif;
   if logl1<logld;
     rho=rhod;
     logl1=logld;
     prop=propvec[4];
   endif;


   /* Table. */

   mask=0~1~1~1;
   let fmt[4,3]=
   "*.*s "  7 7
   "*.*lf " 10 4
   "*.*lf " 17 4
   "*.*lf " 8 4;
   print "Information about calculated rho-values :";
   print "-----------------------------------------";
   cnt=printfm(naamvec~rhovec~loglvec~propvec,mask,fmt);
   print " ";

   /*
   ** In case the rho with the best likelihood is based on to little data,
   ** one chooses the average rho as initial value.
   */

   if (prop<0.05); rho=rhogem;logl1=logl; endif;

   clear rhoa,rhob,rhoc,rhod,rhohulp,rhogem;

endif;


/* Final choice of returns. */

if rows(y)>1;
   alt=(a[1]/abs(a[1]))|zeros(rows(a)-1,1)|
      (b[1]/abs(b[1]))|zeros(rows(a)-1,1)|rho;
   loglik2=loglik(alt,x,y);
   loglik1=loglik(a|b|rho,x,y);
else;
   alt=a|b|rho;
   loglik1=loglik(a|b|rho,x,y);
   loglik2=loglik1;
endif;

/*print "a|b|rho" (a|b|rho)';
print "alt" alt';
print "loglik1" loglik1;
print "loglik2" loglik2;*/

if _file/=2; output on; endif;

retp(a|b|rho,loglik1,alt,loglik2);

endp;




