/*
**                         FILE BPM2.ARC
**
** This file contains the bivariate probit model procedures.
** Part 2.
*/



/* Declaration of globals. */

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


/*
** LOGLIKUP
**
** FORMAT      waarde=loglikup(teta,x,y);
**
** INPUT       teta   :  Vector of parameters (Lx1).
**             x      :  NxL matrix of independant variables.
**             y      :  Nx2 matrix of contingencies.
**
** OUTPUT      waarde :  Scalar, log-likelihood for specified parameters.
**
*/


proc loglikup(teta,x,y);

local p,aa,n;

n=rows(x);
aa=x*teta;
p=zeros(n,2);
p[.,1]=cdfn(-aa);
p[.,2]=1-p[.,1];
p=p.*(p.>0)+(1e-12)*(p.<=0);
retp(sumc(sumc(y.*ln(p))));
endp;




/*
** GRADUP
**
** FORMAT      eerste=gradup(teta,x,y);
**
** INPUT       teta   :  Lx1 vector of parameters.
**             x      :  NxL vector of independant variables.
**             y      :  Nx2 vector of contingencies.
**
** OUTPUT      eerste :  Lx1 vector: gradient of log-likelihood function.
*/



proc gradup(teta,x,y);

local xak,n,l,s,p,aa,d1,q,v1,k;

n=rows(x);
s=-1~1;
l=rows(teta);
p=zeros(n,2);
aa=x*teta;
d1=zeros(l,1);

p[.,1]=cdfn(-aa);
p[.,2]=1-p[.,1];

q=phi1(aa);

v1=(y.*q)./p;

k=1;
do while (k<=l);
   xak=x[.,k]*s;
   d1[k]=sumc(sumc(v1.*xak));
   k=k+1;
endo;

retp(d1);
endp;




/*
** HESSUP
**
** FORMAT      tweede=hessup(teta,x,y);
**
** INPUT       teta   :  Lx1 vector of parameters.
**             x      :  NxL vector of observations.
**             y      :  Nx2 vector of contingencies.
**
** OUTPUT      tweede :  LxL matrix: Hessian of log-likelihood function.
*/



proc hessup(teta,x,y);

local n,s,l,p,aa,d2,q,v2,i,j,xhulp;

n=rows(x);
s=-1~1;
l=rows(teta);
p=zeros(n,2);
aa=x*teta;
d2=zeros(l,l);
p[.,1]=cdfn(-aa);
p[.,2]=1-p[.,1];
q=phi1(aa);
v2=-y.*q.*(s.*p.*aa+q)./(p^2);

i=1;
do while (i<=l);
   j=1;
   do while (j<=l);
      if (i<=j);
         xhulp=x[.,i].*x[.,j];
         d2[i,j]=sumc(sumc(v2.*xhulp));
      else;
         d2[i,j]=d2[j,i];
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;


retp(d2);
endp;





/*
** UNIPROB
**
** FORMAT     {teta1,loglik1,grad,hess}=
**                              uniprob(teta0,x,y);
**
** INPUT      teta0      : Lx1 vector of initial estimates.
**            x          : NxL vector of observations.
**            y          : Nx2 vector of contingencies.
**
** OUTPUT     teta1      : Lx1 vector of estimates for te parameters.
**            loglik1    : Log-likelihood at parameters.
**            grad       : Lx1 vector: gradient at parameter vector.
**            hess       : LxL matrix: hessian at parameter vector.
**
** GLOBALS    _eps       : Tolerance for the gradient (default 1e-3).
**            _maxit     : Upper bound for the number of iterations
**                         (default 30).
*/



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

local teta1,teta2,teller,test,eerste,tweede,ll,infor;

teta1=teta0;
teller=1;
test=1;

do until ((test<=_eps) or (teller>2.5*_maxit));
   eerste=gradup(teta1,x,y);
   tweede=hessup(teta1,x,y);
   trap 1;
   infor=inv(tweede);
   if not(scalerr(infor));
      trap 0;
      teta2=teta1-infor*eerste;
   else;
      trap 0;
      errorlog "Divergence of Newton-Raphson for Univariate Probit.";
      teta1=zeros(rows(teta0),1);
      eerste=-1;
      tweede=-1;
      ll=0;
      goto label;
   endif;
   test=sumc(abs(eerste));
   teta1=teta2;
   ll=loglikup(teta1,x,y);
   if _itscrn; 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 :";
   print teta1'|eerste';
   print " ";
   print "Log-Likelihood at current iteration :" ll;
   print " ";
   print "*****************************************************************"\
   "***************";
   print " ";
   teller=teller+1;
endo;

if (test>_eps); errorlog "The uniprob-test is not satisfied.";endif;

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








