/*
**                         FILE BPM10.ARC
**
** This file contains the bivariate probit model procedures.
** Part 10.
*/


/* Declaration of globals. */

declare _rhonorm ?=0;
declare _transg ?=1;
declare _maxit ?=30;
declare _rpar ?=0;



/*
** BPMONCAR
**
** FORMAT     {teta,loglik,grad,hess}=bpmoncar(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      : (2*L+1)x1 vector; estimates for the parameters.
**            loglik    : Scalar; log-likelihood at estimated parameters.
**            grad      : Gradient at parameter vector.
**            hess      : Hessian at parameter vector.
**
** GLOBALS    _maxit    : Upper bound for the number of iterations.
**            _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.
**
** REMARK     BPMONCAR provides an alternative to determine bivariate probit
**            parameters. It is useful in case Newton-Raphson fails.
**            BPMONCAR uses Monte Carlo methods. It is considerably slower
**            than BIVPROB.
*/



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


local tetamc,factor,max,llmc,t,max,tel,rhomc,ll1,ll2,max2,l,ll,max0,test,
 afg1,rho,mu,tetashow,i,gedaan,sav,uitgev;


if _rpar/=0 and _rpar/=1;
   errorlog "_rpar should be either 0 or 1. We change it to 0.";
   _rpar=0;
endif;

tel=0;
if _rpar==0;
   sav=_rhonorm;
   _rhonorm=0;
endif;
gedaan=0;


/* Determining initial values. */

{max,ll1,max2,ll2}=initval(x,y);
l=rows(max);
l=(l-1)/2;

/* poging om startwaarden van dit proces te verbeteren */

local y1,y2,teta01,teta02,eerste,tweede;

y1=(y[.,1]+y[.,3])~(y[.,2]+y[.,4]);
y2=(y[.,1]+y[.,2])~(y[.,3]+y[.,4]);
screen off;
output off;
{teta01,ll1,eerste,tweede}=uniprob(max[1:l],x,y1);
{teta02,ll2,eerste,tweede}=uniprob(max[l+1:2*l],x,y2);
screen on;
if _file/=2;output on;endif;

max[1:2*l]=teta01|teta02;


/* Preparing values to start iteration process. */

if rows(_transg)/=rows(max);
   _transg=1;
endif;
if rows(_transg)==1 and cols(_transg)==1;
   if _rpar==0;
      _transg=eye(2*l+1);
   elseif _rpar==1;
      _transg=eye(3*l);
   endif;
endif;

factor=2;

max=max';
if _rhonorm==0;
   max[2*l+1]=transrm(max[2*l+1]);
   mu=max[2*l+1];
endif;

if _rpar==0;
   afg1=gradbp(max',x,y);
   ll=loglik(max',x,y);
elseif _rpar==1;
   max=max~zeros(1,l-1);
   afg1=gradbpr(max',x,y);
   ll=loglikr(max',x,y);
endif;

test=sqrt(sumc((abs(afg1))^2));


/* Iterations. */

if _itfile==0;
   output off;
endif;

do until test<0.25*l*rows(x)/20 or tel>=_maxit;
   slot:
   if _rpar==0;
      tetamc=factor*rndn(50,2*l+1);
      if abs(mu)>1.5;
         tetamc[.,2*l+1]=rndn(50,1);
      endif;
   elseif _rpar==1;
      tetamc=factor*rndn(50,3*l);
      if _rhonorm==0;
         if not abs(max[2*l+1:3*l])<=1.5;
print "aanpassingsvector";
print (abs(max[2*l+1:3*l]) .>1.5);
           tetamc[.,2*l+1:3*l]=tetamc[.,2*l+1:3*l].*(abs(max[2*l+1:3*l]) .<=1.5)
               +rndn(50,l).*(abs(max[2*l+1:3*l]) .>1.5);
         endif;
      endif;
   endif;
   tetamc=tetamc+max;
   uitgev=0;
   if tel>2;
      i=1;
      do while i<=2*l and not(gedaan);
         if abs(afg1[i])<=(test/(2.5*cols(max)));
            uitgev=uitgev~i;
            tetamc[.,i]=max[i]*ones(50,1);
         endif;
         i=i+1;
      endo;
      if not(gedaan) and abs(afg1[2*l+1])<=(test/(20*(2*l+1)));
         uitgev=uitgev~(2*l+1);
         tetamc[.,2*l+1]=max[2*l+1]*ones(50,1);
      endif;
   else;
      uitgev=seqa(0,1,2*l+1)';
      tetamc[.,1:2*l]=max[1:2*l].*ones(50,2*l);
   endif;
   if cols(uitgev)>1;
      print " ";
      format /rd 4,0;
      print "Will be kept fixed :";; uitgev[2:cols(uitgev)];
   endif;
   llmc=zeros(50,1);
   t=1;
   if _rpar==0;
      do while t<=50;
            llmc[t]=loglik(tetamc[t,.]',x,y);
         t=t+1;
      endo;
   elseif _rpar==1;
      do while t<=50;
            llmc[t]=loglikr(tetamc[t,.]',x,y);
         t=t+1;
      endo;
   endif;
   tetamc=sortc(-ll~max|-llmc~tetamc,1);
   max0=max;
   max=tetamc[1,2:cols(tetamc)];
   ll=-tetamc[1,1];
   tel=tel+1;
   cls;
   format /rd 5,0;
   print "BIVARIATE PARAMETERS (MONTE  CARLO) :  ITERATION NUMBER" tel;
   print "------------------------------------------------------------";
   format /rd 14,4;
   print " ";
   if max==max0;
      print "Not changed.";
      factor=factor*(0.5);
   else;
      print "Optimized.";
/*      factor=factor*0.65;*/
      if _rpar==0;
         afg1=gradbp(max',x,y);
      elseif _rpar==1;
         afg1=gradbpr(max',x,y);
      endif;
      test=sqrt(sumc((abs(afg1))^2));
   endif;
   print " ";
   print "Parameters/gradient at current iteration :";
   if _rpar==0;
      mu=max[2*l+1];
      rho=transmr(mu);
      tetashow=max[1:2*l]'|rho;
      tetashow=_transg*tetashow;
   elseif _rpar==1;
      tetashow=_transg*max';
   endif;
   print (tetashow)'|afg1';
   print " ";
   if _rpar==0;
      print "                     Corresponding mu :" mu;
   endif;
   print "Log-likelihood for current parameters :" ll;
   print "                        Testing value :" test;
   print "                               Factor :" factor;
   print " ";
   print "************************************************************"\
   "*****************";
   print "  ";
   if test<1 and not(gedaan);
      gedaan=1;
      goto slot;
   endif;
endo;

if _file/=2;
   output on;
endif;

_transg=1;

if _rpar==0;
   _rhonorm=sav;
endif;
if _rpar==0;
   retp(max[1:2*l]'|transmr(max[2*l+1]),ll,afg1,hessbp(max',x,y));
elseif _rpar==1;
   retp(max',ll,afg1,hessbpr(max',x,y));
endif;

endp;

