/*
**                         FILE BPMRS2.ARC
**
** This file contains the bivariate probit model procedures (r x s tables).
** Part 2.
*/



/* Declaration of globals. */

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




/*
** LOGLUPRS
**
** FORMAT      v=logluprs(teta,x,y);
**
** INPUT       teta   :  L*(r-1) parameter vector.
**             x      :  NxL matrix of independent variables. The first column
**                       has to consist completely of ones (to deal with
**                       the constant term).
**             y      :  Nxr matrix of contingencies.
**
** OUTPUT      v      :  Scalar, log-likelihood for specified parameters.
*/


proc logluprs(teta,x,y);

local  a,l,n,q,p,aa,i,r;


/* Definition of working variables. */

r=cols(y);
n=rows(x);
l=cols(x);
a=reshape(teta',r-1,l)';
aa=x*a;


/* Computation of cell probabilities. */

q=zeros(n,r+1);
p=zeros(n,r);

q[.,2:r]=cdfn(aa);
q[.,r+1]=ones(n,1);

i=1;
do while i<=r;
   p[.,i]=q[.,i+1]-q[.,i];
   i=i+1;
endo;

p=p.*(p.>0)+(1e-12)*(p.<=0);

retp(sumc(sumc(y.*ln(p))));
endp;






/*
** GRADUPRS
**
** FORMAT      grad=graduprs(teta,x,y);
**
** INPUT       teta   :  L*(r-1) parameter vector.
**             x      :  NxL matrix of independent variables. The first column
**                       has to consist completely of ones (to deal with
**                       the constant term).
**             y      :  Nxr matrix of contingencies.
**
** OUTPUT      grad   :  L*(r-1) gradient vector.
*/


proc graduprs(teta,x,y);

local r,n,a,l,q,aa,i,grad,p;

/* Definition of working variables. */

r=cols(y);
l=cols(x);
n=rows(x);
a=reshape(teta',r-1,l)';
aa=x*a;

/* Computation of cell probabilities. */

q=zeros(n,r+1);
p=zeros(n,r);

q[.,2:r]=cdfn(aa);
q[.,r+1]=ones(n,1);

i=1;
do while i<=r;
   p[.,i]=q[.,i+1]-q[.,i];
   i=i+1;
endo;

p=p.*(p .gt 1e-9)+(1e-9).*(p .le 1e-9);

/* Computation of gradient vector. */

q=(y[.,1:r-1]./p[.,1:r-1]-y[.,2:r]./p[.,2:r]).*phi1(aa);
grad=x'*q;

retp(vec(grad));
endp;



/*
** HESSUPRS
**
** FORMAT      hess=hessuprs(teta,x,y);
**
** INPUT       teta   :  L*(r-1) parameter vector.
**             x      :  NxL matrix of independent variables. The first column
**                       has to consist completely of ones (to deal with
**                       the constant term).
**             y      :  Nxr matrix of contingencies.
**
** OUTPUT      hess   :  (L*(r-1))x(L*(r+1)) Hessian matrix.
*/


proc hessuprs(teta,x,y);

local r,n,a,l,q,aa,i,grad,p,j,p2,hulp,hess,pp,qq;


/* Definition of working variables. */

r=cols(y);
l=cols(x);
n=rows(x);
a=reshape(teta',r-1,l)';
aa=x*a;


/* Computation of cell probabilities. */

q=zeros(n,r+1);
p=zeros(n,r);

q[.,2:r]=cdfn(aa);
q[.,r+1]=ones(n,1);

i=1;
do while i<=r;
   p[.,i]=q[.,i+1]-q[.,i];
   i=i+1;
endo;

p2=p.*p;

p=p.*(p .gt 1e-9)+(1e-9).*(p .le 1e-9);
p2=p2.*(p2 .gt 1e-9)+(1e-9).*(p2 .le 1e-9);


/* Computation of auxiliary matrix. */

hess=zeros(l*(r-1),l*(r-1));

i=1;
do while i<=r-1;
   j=1;
   do while j<=r-1;
      if i==j;
         hulp=-aa[.,i].*phi1(aa[.,i]).*(y[.,i]./p[.,i]-y[.,i+1]./p[.,i+1])
              -(y[.,i]./p2[.,i]+y[.,i+1]./p2[.,i+1])
              .*phi1(aa[.,i]).*phi1(aa[.,i]);
      elseif i+1==j;
         hulp=y[.,i+1]./p2[.,i+1].*phi1(aa[.,i]).*phi1(aa[.,i+1]);
      elseif i-1==j;
         hulp=y[.,i]./p2[.,i].*phi1(aa[.,i]).*phi1(aa[.,i-1]);
      endif;
      if abs(i-j)<=1;
         pp=1;
         do while pp<=l;
            qq=1;
            do while qq<=l;
               hess[(i-1)*l+pp,(j-1)*l+qq]=sumc(hulp.*x[.,pp].*x[.,qq]);
               qq=qq+1;
            endo;
            pp=pp+1;
         endo;
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;


retp(hess);
endp;






/*
** UPROBRS
**
** FORMAT     {teta1,logl,grad,hess}=
**                              uprobrs(teta0,x,y);
**
** INPUT      teta0      : L*(r-1)x1 vector of initial estimates.
**            x          : NxL vector of observations.
**            y          : Nxr matrix of contingencies.
**
** OUTPUT     teta1      : L*(r-1)x1 vector of estimated parameters.
**            logl       : Log-likelihood at estimated parameters.
**            grad       : L*(r-1)x1 gradient vector.
**            hess       : L*(r-1)x L*(r-1) Hessian matrix.
**
** GLOBALS    _eps       : Tolerance for the gradient (default 1e-3).
**            _maxit     : Upper bound for the number of iterations
**                         (default 30).
*/



proc (4)=uprobrs(teta0,x,y);

local teta1,teta2,teller,test,eerste,tweede,ll,infor,r,l,tetahulp;


/* Definition of working variables. */

r=cols(y);
l=cols(x);


/* Newton-Raphson iterations. */

teta1=teta0;
tetahulp=teta1;
if l>1;
   teta1=transeo1(teta1',r,0)';
endif;
teller=1;
test=1;

do until ((test<=_eps) or (teller>2.5*_maxit));
   if l>1;
      eerste=transeo1(graduprs(tetahulp,x,y)',r,1)';
      tweede=hessuprs(tetahulp',x,y);
      tweede=transeo1(transeo1(tweede,r,1)',r,1)';
   else;
      eerste=graduprs(tetahulp,x,y);
      tweede=hessuprs(tetahulp,x,y);
   endif;
   trap 1;
   infor=inv(tweede);
   if not(scalerr(infor));
      trap 0;
      teta2=teta1-infor*eerste;
   else;
      trap 0;
      locate 21,1;
      errorlog "Divergence of Newton-Raphson for Univariate Probit.";
      teta1=zeros(rows(teta1),1);
      eerste=-1;
      tweede=-1;
      ll=0;
      goto label;
   endif;
   test=sumc(abs(eerste));
   teta1=teta2;
   if l>1;
      tetahulp=transoe1(teta1',r)';
   else;
      tetahulp=teta1;
   endif;
   ll=logluprs(tetahulp,x,y);
   if _itscrn==1; cls; endif;
   format /rd 5,0;
   print "UNIVARIATE PARAMETERS :  ITERATION NUMBER" teller;
   print "----------------------------------------------";
   print " ";
   format /rd 14,-log(_eps)+2;
   print "Parameters/gradient at current iteration :";
   if l>1;
      print reshape(tetahulp',r-1,l)'|reshape(transoe1(eerste',r),r-1,l)';
   else;
      print reshape(tetahulp',r-1,l)'|reshape(eerste,r-1,l)';
   endif;
   print " ";
   print "Log-Likelihood at current iteration :" ll;
   print " ";
   print "*****************************************************************"\
   "***************";
   print " ";
   teller=teller+1;
endo;


label:
retp(teta1,ll,eerste,tweede);
endp;







