/*
**                         FILE BPM7.ARC
**
** This file contains the bivariate probit model procedures.
** Part 7.
*/



/* Declaration of globals. */

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



/* Setting of libraries needed. */

library a:\bpm,pgraph,dutil,stat;



/*
** GRADPL
**
** FORMAT      grad=gradpl(teta,rho,x,y);
**
** INPUT       teta  :  (2*L)x1 vector of parameters.
**             rho   :  Scalar, correlation coefficient.
**             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      grad  :  (2*L)x1 gradient vector.
**
** REMARK      The maximum of the likelihood function is calculated with
**             given rho as a boundary function. Thus, the gradient is only
**             taken with respect to the parameters for the first and second
**             variable.
*/


proc gradpl(teta,rho,x,y);

local a,b,l,n,aa,bb,p,za,zb,nr,k;
local xak,xbk,da,db,pp;
local eersteaf,rho,s1,s2;
local va,vb;


/* Initialising some values. */

n=rows(x);
let s1=-1,1,-1,1;
let s2=-1,-1,1,1;


/* Splitting  the parameter vector into its components. */

l=rows(teta);
l=l/2;
a=teta[1:l];
b=teta[l+1:2*l];

aa=x*a;
bb=x*b;


/* Matrices to calculate first order derivatives. */

p=zeros(n,4);
p[.,1]=cdfbvn(-aa,-bb,rho);
p[.,2]=cdfn(-bb)-p[.,1];
p[.,3]=cdfn(-aa)-p[.,1];
p[.,4]=ones(n,1)-p[.,1]-p[.,2]-p[.,3];

va=zeros(n,4);
vb=zeros(n,4);

va[.,1]=y[.,1].*eerstex(-aa,-bb,rho);
va[.,2]=y[.,2].*eerstex(aa,-bb,-rho);
va[.,3]=y[.,3].*eerstex(-aa,bb,-rho);
va[.,4]=y[.,4].*eerstex(aa,bb,rho);
za=va./p;
clear va;

vb[.,1]=y[.,1].*eerstey(-aa,-bb,rho);
vb[.,2]=y[.,2].*eerstey(aa,-bb,-rho);
vb[.,3]=y[.,3].*eerstey(-aa,bb,-rho);
vb[.,4]=y[.,4].*eerstey(aa,bb,rho);
zb=vb./p;
clear vb;


/* Calculating components of gradient vector. */

xak=zeros(l,1);
xbk=zeros(l,1);
da=zeros(l,1);
db=zeros(l,1);

k=1;
do while (k<=l);
   xak=s1*(x[.,k]');
   xbk=s2*(x[.,k]');
   da[k]=sumc(sumc(za.*xak'));
   db[k]=sumc(sumc(zb.*xbk'));
   k=k+1;
endo;

clear za,zb;

eersteaf=da|db;

retp(eersteaf);
endp;




/*
** HESSPL
**
** FORMAT      hess=hesspl(teta,rho,x,y);
**
** INPUT       teta  :  (2*L)x1 vector of parameters.
**             rho   :  Scalar, correlation coefficient.
**             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      hess  :  (2*L)x(2*L) gradient vector.
**
** REMARK      The maximum of the likelihood function is calculated with
**             given rho as a boundary function. Thus, the hessian is only
**             taken with respect to the parameters for the first and second
**             variable.
*/



proc hesspl(teta,rho,x,y);

local s1,s2;
local a,b,l,n,aa,bb,p,pp,rho;
local afgx,afgy,vaa,vbb,vab;
local ddaa,ddbb,ddab,tweedeaf;
local i,j,xhulp;


/* Initialising some values. */

n=rows(x);
let s1=-1,1,-1,1;
let s2=-1,-1,1,1;


/* Splitting the parameter vector into its components. */

l=rows(teta);
l=l/2;
a=teta[1:l];
b=teta[l+1:2*l];

aa=x*a;
bb=x*b;


/* Matrices to calculate first and second order derivatives. */

p=zeros(n,4);
p[.,1]=cdfbvn(-aa,-bb,rho);
p[.,2]=cdfn(-bb)-p[.,1];
p[.,3]=cdfn(-aa)-p[.,1];
p[.,4]=ones(n,1)-p[.,1]-p[.,2]-p[.,3];

pp=p^2;

afgx=zeros(n,4);
afgx[.,1]=eerstex(-aa,-bb,rho);
afgx[.,2]=eerstex(aa,-bb,-rho);
afgx[.,3]=eerstex(-aa,bb,-rho);
afgx[.,4]=eerstex(aa,bb,rho);

afgy=zeros(n,4);
afgy[.,1]=eerstey(-aa,-bb,rho);
afgy[.,2]=eerstey(aa,-bb,-rho);
afgy[.,3]=eerstey(-aa,bb,-rho);
afgy[.,4]=eerstey(aa,bb,rho);

vaa=zeros(n,4);
vaa[.,1]=y[.,1].*(p[.,1].*tweedexx(-aa,-bb,rho) - (afgx[.,1])^2);
vaa[.,2]=y[.,2].*(p[.,2].*tweedexx(aa,-bb,-rho)-(afgx[.,2])^2);
vaa[.,3]=y[.,3].*(p[.,3].*tweedexx(-aa,bb,-rho)-(afgx[.,3])^2);
vaa[.,4]=y[.,4].*(p[.,4].*tweedexx(aa,bb,rho)-(afgx[.,4])^2);
vaa=vaa./pp;

vbb=zeros(n,4);
vbb[.,1]=y[.,1].*(p[.,1].*tweedeyy(-aa,-bb,rho) - (afgy[.,1])^2);
vbb[.,2]=y[.,2].*(p[.,2].*tweedeyy(aa,-bb,-rho)-(afgy[.,2])^2);
vbb[.,3]=y[.,3].*(p[.,3].*tweedeyy(-aa,bb,-rho)-(afgy[.,3])^2);
vbb[.,4]=y[.,4].*(p[.,4].*tweedeyy(aa,bb,rho)-(afgy[.,4])^2);
vbb=vbb./pp;

vab=zeros(n,4);
vab[.,1]=y[.,1].*(p[.,1].*tweedexy(-aa,-bb,rho) - afgx[.,1].*afgy[.,1]);
vab[.,2]=y[.,2].*(p[.,2].*tweedexy(aa,-bb,-rho)-afgx[.,2].*afgy[.,2]);
vab[.,3]=y[.,3].*(p[.,3].*tweedexy(-aa,bb,-rho)-afgx[.,3].*afgy[.,3]);
vab[.,4]=y[.,4].*(p[.,4].*tweedexy(aa,bb,rho)-afgx[.,4].*afgy[.,4]);
vab=vab./pp;

clear afgx,afgy;


/* Calculation of ddaa,ddbb,ddab,ddrr,ddar,ddbr. */

ddaa=zeros(l,l);

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

clear vaa;

ddbb=zeros(l,l);

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

clear vbb;

ddab=zeros(l,l);

i=1;
do while i<=l;
   j=1;
   do while (j<=l);
      if (i<=j);
         xhulp=x[.,i].*x[.,j];
         xhulp=((s1.*s2)*xhulp')';
         ddab[i,j]=sumc(sumc(vab.*xhulp));
      else;
         ddab[i,j]=ddab[j,i];
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;

clear vab;

tweedeaf=(ddaa~ddab)|(ddab'~ ddbb);

retp(tweedeaf);
endp;




/*
** PROFLIK
**
** FORMAT   {rho,par,pl,tvec,pvec}=
**            proflik(x,y,indvar,depvar,stap,min,max);
**
** INPUT    x        : NxL or Nx(L-1) matrix of independent variables.
**                     (L-1 if the first column isn't a unit column; in that
**                     case a unit column is appended and L-1 is incremented).
**          y        : Nx4 matrix of contingencies.
**          indvar   : Lx1 or (L-1)x1 character vector; names of independent
**                     variables.
**          depvar   : 2x1 character vector; names of dependent variables.
**          stap     : Scalar,stepsize for rho.
**          min      : Minimum value for rho.
**          max      : Maximum value for rho.
**
** OUTPUT   rho      : Vector of rho values used (Kx1; K=2/stap-1).
**          par      : Kx4*L matrix of parameters and standard errors.
**          pl       : Kx1 vector of log-likelihoods.
**          tvec     : Kx2*L matrix of t-values.
**          pvec     : Kx2*L matrix of p-values.
**
** GLOBALS  _draw : 0: Pictures are made by the program.
**                  1: Pictures are not made by the program (default).
**          _eps     : Tolerance for the gradient (default 0.001).
**          _maxit   : Upper bound for the number of iterations
**                     (default 30).
**
** REMARK   A vector of rho values is created between -1 and +1 with steps
**          given by the variable stap.  For each rho the log-likelihood
**          function is maximalized. The log-likelihood and the parameters
**          are plotted against rho.
*/



proc (5)= proflik(x,y,indvar,depvar,stap,min,max);

local kol1,hulp,ct,meanx,stdx,meany,stdy,gem,afw,xtest,xtest1,trans,teta0,pl,
 teta1,tvec,pvec,par,rho,teller,test,afg1,afg2,teta2,stde,tval,pval,rhonum,a,
 b,ll,l,ind1,ind2,titel,l1,l2,teta0a,str,titel,transg,infor;


format /rd 10,4;
str="   * ";
titel="    Number       Rho    Log-Lik.    *     Number       Rho    Log-Lik.";

/* Checking for faulty data. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=4);
   errorlog "y does not represent a table of contingencies.";
   end;
elseif not (y>=0);
   errorlog "y contains negative values.";
   end;
endif;

kol1=x[.,1];
if kol1/=ones(rows(x),1);
   errorlog "The first colum of x does noet consist of ones, we add one.";
   x=ones(rows(x),1)~x;
endif;

if stap<=0 or stap>=2;
   errorlog "Stepsize must be in the interval ]0,2[";
   end;
endif;


/* Dealing with missing values. */

if ismiss(x) or ismiss(y);
   hulp=y~x;
   hulp=packr(hulp);
   y=hulp[.,1:4];
   x=hulp[.,5:cols(hulp)];
   clear hulp;
endif;


/* Means and standard errors of the variables. */

ct=sumc(y);
meanx=sumc(x[.,2:cols(x)].*sumc(y'))/sumc(ct);
stdx=sqrt(sumc(((x[.,2:cols(x)]-meanx')^2).*sumc(y'))/(sumc(ct)-1));
meany=ct;
meany=((meany[2]+meany[4])|(meany[3]+meany[4]))/sumc(meany);
stdy=zeros(2,1);
stdy[1]=sqrt((ct[1]+ct[3])*(meany[1]^2)+(ct[2]+ct[4])*(1-meany[1])^2);
stdy[2]=sqrt((ct[1]+ct[2])*(meany[2]^2)+(ct[3]+ct[4])*(1-meany[2])^2);
stdy=stdy/sqrt(sumc(ct)-1);


/*
** Replacing x by a standardised copy; creating values to transform
** the standardised parameters thus found into non-standardised ones.
*/

gem=meanc(x[.,2:cols(x)]);
afw=stdc(x[.,2:cols(x)]);
xtest=(x .==1) + (x .==0);
xtest1=rows(xtest);
xtest=sumc(xtest);
xtest=(xtest .== xtest1);
gem=gem.*(xtest[2:rows(xtest)] .==0);
afw=afw.*(xtest[2:rows(xtest)] .==0) + (xtest[2:rows(xtest)]  .==1);
x[.,2:cols(x)]=(x[.,2:cols(x)]-gem')./afw';
trans=zeros(rows(gem),rows(gem));
trans=diagrv(trans,1./afw);
trans=zeros(rows(trans),1)~trans;
trans=(1~(-gem./afw)')|trans;
transg=trans~zeros(rows(trans),cols(trans))|zeros(rows(trans),cols(trans))~
trans;


/* Calling INITVAL for initial estimates. */

screen off;
{teta0,l1,teta0a,l2}=initval(x,y);
screen on;
teta0=teta0[1:rows(teta0)-1];
l=rows(teta0)/2;


/* Newton-Raphson iterations. */

print titel;
if min==-1; min=-0.95; endif;
if max==1; max=0.95;endif;
teta1=teta0;
pl=0;
tvec=0;
pvec=0;
rhonum=1;
par=zeros(1,4*l);
rho=min;

do while (rhonum<=(round((max-min)/stap)+1));
   teller=1;
   test=1;

print rhonum rho;;

   do until ((test<_eps) or (teller>_maxit));
      afg1=gradpl(teta1,rho,x,y);
      afg2=hesspl(teta1,rho,x,y);
      teta2=teta1-inv(afg2)*afg1;
      test=sumc(abs(afg1));
      teta1=teta2;
      teller=teller+1;
   endo;


   a=teta1[1:l];
   b=teta1[l+1:2*l];
   a=trans*a;
   b=trans*b;
   infor=transg*inv(-afg2)*transg';
   stde=sqrt(diag(infor));
   tval=(a|b)./stde;
   pval=2*(1-cdfn(abs(tval)));
   ll=loglik(teta1|rho,x,y);
   if (rhonum%2==0);
      print ll;
   else;
      print ll str;;
   endif;
   pl=pl|ll;
   par=par|(a'~b'~stde');
   tvec=tvec|tval;
   pvec=pvec|pval;
   rho=rho+stap;
   rhonum=rhonum+1;

endo;

pl=pl[2:rows(pl)];
tvec=tvec[2:rows(tvec)];
pvec=pvec[2:rows(pvec)];
par=par[2:rows(par),.];
rho=seqa(min,stap,round((max-min)/stap)+1);


/* Graphics section. */


if _draw;

   clear x,y;


   title("Profile Log-Likelihood");
   xlabel("Correlation Coefficient");
   ylabel("Log-Likelihood at alpha(rho)");
   print "we gaan pl tegen rho uitzetten";
   xy(rho,pl);

   title("Evolution for a Parameter");
   xlabel("Rho");
   ylabel("Parameter");


   ind1=1;
   do while (ind1<=2);
      ind2=1;
      do while (ind2<=l);
         titel="Coeff of ";
         titel=titel$+indvar[ind2]$+" for "$+depvar[ind1];
         title(titel);
         teller=(ind1-1)*l+ind2;
         xy(rho,par[.,teller]);
         ind2=ind2+1;
      endo;
      ind1=ind1+1;
   endo;

endif;


retp(rho,par,pl,tvec,pvec);
endp;



