/*
**                         FILE BPM1.ARC
**
** This file contains the bivariate probit model procedures.
** Part 1.
*/


/* Declaration of globals. */

declare _ordinar ?=1;
declare _rpar ?=0;



/*
** EERSTEX   EERSTEY   EERSTER
**
** FORMAT   v=eerste.(x,y,rho);
**
** INPUT    x,y : two variables (Nx1).
**          rho : correlation coefficient (1x1 or Nx1).
**
** OUTPUT   Nx1 vector of first derivatives of the Bivariate Normal Integral.
*/


proc eerstex(x,y,rho);


local result;
result=phi1(x).*cdfn((y-rho.*x)./sqrt(1-rho.*rho));
retp(result);
endp;




proc eerstey(x,y,rho);

local result;
result=phi1(y).*cdfn((x-rho.*y)./sqrt(1-rho.*rho));
retp(result);
endp;




proc eerster(x,y,rho);

local result;
retp(phi2(x,y,rho));
endp;




/*
** TWEEDEXX   TWEEDEXY   TWEEDEYY   TWEEDEXR   TWEEDEYR   TWEEDERR
**
** FORMAT   v=tweede..(x,y,rho);
**
** INPUT    x,y : two variables (Nx1).
**          rho : correlation coefficient (1x1 or Nx1).
**
** OUTPUT   Nx1 vector of second order derivatives of the Bivariate Normal
**          Integral.
*/



proc tweedexx(x,y,rho);

local result;
result=rho./(sqrt(1-rho.*rho));
result=result.*phi1((-rho.*x+y)./sqrt(1-rho.*rho));
result=result+x.*cdfn((-rho.*x+y)./sqrt(1-rho.*rho));
result=-result.*phi1(x);
retp(result);
endp;




proc tweedeyy(x,y,rho);

local result;
result=rho./(sqrt(1-rho.*rho));
result=result.*phi1((-rho.*y+x)./sqrt(1-rho.*rho));
result=result+y.*cdfn((-rho.*y+x)./sqrt(1-rho.*rho));
result=-result.*phi1(y);
retp(result);
endp;




proc tweedexy(x,y,rho);

retp(phi2(x,y,rho));
endp;




proc tweedexr(x,y,rho);

local result;
result=(rho.*y-x)./(1-rho.*rho);
result=result.*phi2(x,y,rho);
retp(result);
endp;




proc tweedeyr(x,y,rho);

local result;
result=(rho.*x-y)./(1-rho.*rho);
result=result.*phi2(x,y,rho);
retp(result);
endp;




proc tweederr(x,y,rho);

local result;
result=rho.*(1-rho.*rho)+x.*y.*(1+rho.*rho)-(x.*x+y.*y).*rho;
result=result./((1-rho.*rho)^2);
result=result.*phi2(x,y,rho);
retp(result);
endp;








/*
** PHI1
**
** FORMAT    t=phi1(x);
**
** INPUT     x: Nx1 vector of variables.
**
** OUTPUT    t: Nx1 vector of densities.
**
** REMARK    Calculates the univariate normal density function.
*/

proc phi1(x);

local result;
result=1/sqrt(2*pi);
result=result*exp(-.5*x.*x);
retp(result);
endp;




/*
** PHI2
**
** FORMAT     t=phi2(x,y,rho);
**
** INPUT      x   : Nx1 vector of first variables.
**            y   : Nx1 vector of second variables.
**            rho : Scalar of Nx1 vector of correlation coefficients.
**
** OUTPUT     t   : Nx1 vector of densities.
**
** REMARK     Calculates the bivariate normal density function.
*/


proc phi2(x,y,rho);

local result,exponent;
result=1/(2*pi);
result=result/(sqrt(1-rho.*rho));
exponent=((-x.*x+2*rho.*x.*y-y.*y)*0.5)./(1-rho.*rho);
result=result.*exp(exponent);
retp(result);
endp;





/*
** TRANSRM
**
** FORMAT  mu=transrm(rho);
**
** INPUT   rho : Scalar or column vector of correlation coefficient.
**
** OUTPUT  mu  : Scalar or column vector of transformed correlations.
**
** REMARK  Rho is transformed into a real variable, stretching from minus
**         infinity to plus infinity, to avoid the correlation surpassing
**         the interval from -1 to +1.
*/


proc transrm(rho);

retp(-ln((1-rho)./(1+rho)));
endp;





/*
** TRANSMR
**
** FORMAT  rho=transmr(mu);
**
** INPUT   mu  : Scalar or column vector of transformed correlations.
**
** OUTPUT  rho : Scalar or column vector of correlation coefficients.
**
** REMARK  Rho is transformed into a real variable, stretching from minus
**         infinity to plus infinity, to avoid the correlation surpassing
**         the interval from -1 to +1. This function is the inverse of transrm.
*/



proc transmr(mu);

retp((1-exp(-mu))./(1+exp(-mu)));
endp;





/*
** TRANSLS
**
** FORMAT  teta0=transls(teta1,restric);
**
** INPUT   teta1   :   Nx(M*K) matrix.
**         restric :   MxK matrix of booleans indicating which columns
**                     of teta1 are to be kept.
**
** OUTPUT  teta0   :   Nx(sumc(sumc(restric))), selection of teta1, according
**                     to restric.
**
*/



proc transls(teta1,restric);

local teta0,i;

restric=vec(restric);
teta0=zeros(rows(teta1),1);
i=1;
do while i<=cols(teta1);
   if restric[i];
      teta0=teta0~teta1[.,i];
   endif;
   i=i+1;
endo;
teta0=teta0[.,2:cols(teta0)];

retp(teta0);
endp;






/*
** TRANSSL
**
** FORMAT  teta1=transsl(teta0,restric);
**
** INPUT   teta0   :   Nx(sumc(sumc(restric))) matrix.
**         restric :   MxK matrix of booleans, indicating in which columns
**                     of teta1 the columns of teta0 are to be placed.
**
** OUTPUT  teta1   :   Nx(M*K) matrix, enlarged version of teta0, following
**                     restric.
*/



proc transsl(teta0,restric);

local teta1,i;

restric=vec(restric);
teta1=zeros(rows(teta0),rows(restric));

i=1;
do while i<=cols(teta1);
   if restric[i];
      teta1[.,i]=teta0[.,sumc(restric[1:i])];
   endif;
   i=i+1;
endo;

retp(teta1);
endp;







/*
** TRANSEO
**
** FORMAT  teta0=transeo(teta1,r,s,trans);
**
** INPUT   teta1   :   Mx(L*(r+s-2)+1) matrix.
**         r,s     :   Format of tables.
**         trans   :   Indicates wether the shortening of teta1 into teta0
**                     has to be done either by replacing the various
**                     coefficients of the covariates by a single copy
**                     (trans==0) or by summing corresponding components
**                     (trans==1). The first alternative typically is meant
**                     for the parameters, the second one is meant for
**                     gradients en Hessian matrices.
**
** OUTPUT  teta0   :   Mx(r+s+2*l-3), shortening of teta1.
**
** REMARK  The name of the procedure indicates a TRANSformation of the
**         Extended model to the Ordinary model.
*/



proc transeo(teta1,r,s,trans);

local teta0,i,l,selec,tel;


if _ordinar;

   if _rpar<=0;
      l=(cols(teta1)-1)/(r+s-2);
   elseif _rpar>=1;
      l=cols(teta1)/(r+s-1);
   endif;

   if l>1;

   if trans==0;
      selec=seqa(1,l,r-1);
      teta0=teta1[.,selec];
      teta0=teta0~teta1[.,2:l];
      selec=seqa(1,l,s-1);
      selec=selec+l*(r-1);
      teta0=teta0~teta1[.,selec];
      teta0=teta0~teta1[.,(r-1)*l+2:r*l];
      if _rpar<=0;
         teta0=teta0~teta1[.,(r+s-2)*l+1];
      elseif _rpar>=1;
         teta0=teta0~teta1[.,(r+s-2)*l+1:(r+s-1)*l];
      endif;
   elseif trans==1;
      selec=seqa(1,l,r-1);
      teta0=teta1[.,selec];
      teta0=teta0~teta1[.,2:l];
      tel=2;
      do while tel<=r-1;
         teta0[.,cols(teta0)-l+2:cols(teta0)]=
         teta0[.,cols(teta0)-l+2:cols(teta0)]+teta1[.,(tel-1)*l+2:tel*l];
         tel=tel+1;
      endo;
      selec=seqa(1,l,s-1);
      selec=selec+l*(r-1);
      teta0=teta0~teta1[.,selec];
      teta0=teta0~teta1[.,(r-1)*l+2:r*l];
      tel=2;
      do while tel<=s-1;
         teta0[.,cols(teta0)-l+2:cols(teta0)]=
       teta0[.,cols(teta0)-l+2:cols(teta0)]+teta1[.,(r+tel-2)*l+2:(r+tel-1)*l];
         tel=tel+1;
      endo;
      if _rpar<=0;
         teta0=teta0~teta1[.,(r+s-2)*l+1];
      elseif _rpar>=1;
         teta0=teta0~teta1[.,(r+s-2)*l+1:(r+s-1)*l];
      endif;
   endif;

   elseif l==1;

      teta0=teta1;

   endif;

   retp(teta0);

else;

   retp(teta1);

endif;

endp;






/*
** TRANSOE
**
** FORMAT  teta1=transoe(teta0,r,s);
**
** INPUT   teta0   :   Mx(r+s+2*l-3), shortening of teta1.
**         r,s     :   Format of tables.
**
** OUTPUT  teta1   :   Mx(L*(r+s-2)+1) matrix, teta0 prolonged.
**
** REMARK  The name of the procedure indicates a TRANSformation of the
**         Ordinary model to the Extended model.
*/



proc transoe(teta0,r,s);

local teta1,l,tel;

if _ordinar;

   if _rpar<=0;
      l=(cols(teta0)+3-r-s)/2;
   elseif _rpar>=1;
      l=(cols(teta0)+4-r-s)/3;
   endif;

   if l>1;

   if _rpar<=0;
      teta1=zeros(rows(teta0),l*(r+s-2)+1);
   elseif _rpar>=1;
      teta1=zeros(rows(teta0),l*(r+s-1));
   endif;
   tel=1;
   do while tel<=r-1;
      teta1[.,(tel-1)*l+1]=teta0[.,tel];
      teta1[.,(tel-1)*l+2:tel*l]=teta0[.,r:r+l-2];
      tel=tel+1;
   endo;
   tel=1;
   do while tel<=s-1;
      teta1[.,(r+tel-2)*l+1]=teta0[.,l+r+tel-2];
      teta1[.,(r+tel-2)*l+2:(r+tel-1)*l]=teta0[.,r+l+s-2:r+s+2*l-4];
      tel=tel+1;
   endo;
   if _rpar<=0;
      teta1[.,l*(r+s-2)+1]=teta0[.,2*l+r+s-3];
   elseif _rpar>=1;
      teta1[.,l*(r+s-2)+1:l*(r+s-1)]=teta0[.,2*l+r+s-3:3*l+r+s-4];
   endif;

   elseif l==1;

      teta1=teta0;

   endif;

   retp(teta1);

else;

   retp(teta0);

endif;

endp;






/*
** TRANSEO1
**
** FORMAT  teta0=transeo1(teta1,r,trans);
**
** INPUT   teta1   :   Mx(L*(r-1)) matrix.
**         r       :   Format of tables.
**         trans   :   Indicates wether the shortening of teta1 into teta0
**                     has to be done either by replacing the various
**                     coefficients of the covariates by a single copy
**                     (trans==0) or by summing corresponding components
**                     (trans==1). The first alternative typically is meant
**                     for the parameters, the second one is meant for
**                     gradients en Hessian matrices.
**
** OUTPUT  teta0   :   Mx(r+l-2), shortening of teta1.
**
** REMARK  The name of the procedure indicates a TRANSformation of the
**         Extended model to the Ordinary model for the univariate case(1).
*/



proc transeo1(teta1,r,trans);

local teta0,i,l,selec,tel;

if _ordinar;

   l=cols(teta1)/(r-1);

   if trans==0;
      selec=seqa(1,l,r-1);
      teta0=teta1[.,selec];
      teta0=teta0~teta1[.,2:l];
   elseif trans==1;
      selec=seqa(1,l,r-1);
      teta0=teta1[.,selec];
      teta0=teta0~teta1[.,2:l];
      tel=2;
      do while tel<=r-1;
         teta0[.,cols(teta0)-l+2:cols(teta0)]=
         teta0[.,cols(teta0)-l+2:cols(teta0)]+teta1[.,(tel-1)*l+2:tel*l];
         tel=tel+1;
      endo;
   endif;

   retp(teta0);

else;

   retp(teta1);

endif;

endp;






/*
** TRANSOE1
**
** FORMAT  teta1=transoe1(teta0,r);
**
** INPUT   teta0   :   Mx(r+2*l-2), shortening of teta1.
**         r       :   Format of tables.
**
** OUTPUT  teta1   :   Mx(L*(r-1)) matrix, teta0 prolonged.
**
** REMARK  The name of the procedure indicates a TRANSformation of the
**         Ordinary model to the Extended model, univariate case (1).
*/



proc transoe1(teta0,r);

local teta1,l,tel;

if _ordinar;

   l=cols(teta0)+2-r;

   teta1=zeros(rows(teta0),l*(r-1));
   tel=1;
   do while tel<=r-1;
      teta1[.,(tel-1)*l+1]=teta0[.,tel];
      teta1[.,(tel-1)*l+2:tel*l]=teta0[.,r:r+l-2];
      tel=tel+1;
   endo;
   teta1[.,cols(teta1)]=teta0[.,cols(teta0)];

   retp(teta1);

else;

   retp(teta0);

endif;

endp;




/*
** CDFNINV
**
** FORMAT   y=cdfninv(x);
**
** INPUT    x  : Vector of values between 0 and 1.
**
** OUTPUT   y  : Vector of corresponding inverse cumulative distribution
**               function values.
**
*/



proc cdfninv(x);


local w0,w1,i;


/* Testing for faulty input. */

x=x.*(x.>0.000000001 .and x.<0.99999999)+0.000000001.*(x.<=0.000000001)
  +0.99999999.*(x.>=0.99999999);


/* Newton-iteration. */

i=1;
w0=-5;
w1=0;
do while i<=10 and maxc((abs(w0-w1)))>=1e-10;
   w0=w1;
   w1=w0+(x-cdfn(w0))./phi1(w0);
   i=i+1;
endo;

retp(w1);
endp;


