/*
**                         FILE BPM11.ARC
**
** This file contains the bivariate probit model procedures.
** Part 11.
*/





declare _eps ?=0.001;
declare _maxit ?=30;



/*
** RT
**
** FORMAT       rho=rt(tetra0,y);
**
** INPUT        tetra0 :  3x1 vector of cutpoints and correlation for a 2x2
**                        table.  The cutpoints are considered to be exact.
**                        The correlation is an initial value to be optimised.
**                   y :  Contingency table, ordered as a 4x1 vector.
**
** OUTPUT       rho    :  Scalar, optimisation for the correlation coefficient.
**
** GLOBALS      _eps   :  Criterion to end the Newton-Raphson iterations.
**              _maxit :  Upper bound for the number of iterations.
**
** REMARK       Rt computes the tetrachoric correlation coefficient for a 2x2
**              contingency table, using Newton-Raphson iteration.
*/



proc rt(tetra0,y);

local p,teller,eenr,s12,pp,rho,vrho,afgr,vrr,tweer,aa,bb;


/* Initialising some values. */

aa=tetra0[1];
bb=tetra0[2];
rho=tetra0[3];
teller=1;
eenr=1;
s12=1~-1~-1~1;


/* Newton-Raphson iterations. */

do while (teller<=_maxit) and (sumc(abs(eenr)) > _eps);

   p=zeros(1,4);
   p[1]=cdfbvn(-aa,-bb,rho);
   p[2]=cdfn(-bb)-p[1];
   p[3]=cdfn(-aa)-p[1];
   p[4]=1-p[1]-p[2]-p[3];
   p=p.*(p .> 1e-12) + (1e-12).*(p .<=1e-12);

   pp=p^2;

   vrho=zeros(1,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(1,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(1,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;
   if rho>=1;
      rho=1-1e-4;
   elseif rho<=-1;
      rho=-1+1e-4;
   endif;
   teller=teller+1;

endo;

retp(rho);
endp;
