/*
**                         FILE BPM4.ARC
**
** This file contains the bivariate probit model procedures.
** Part 4.
*/



/* Declaration of globals. */

declare _rhonorm ?=1;



/*
** LOGLIK
**
** FORMAT    v=loglik(teta,x,y);
**
** INPUT     teta   :  (2*L+1)x1 vector of parameters.
**           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    v      :  Scalar, computed log-likelihood.
*/


proc loglik(teta,x,y);

local a,b,rho,l,aa,bb,p,n,result,mu;



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

l=rows(teta);
n=rows(x);
if _rhonorm;
   rho=teta[l];
else;
   mu=teta[l];
   rho=transmr(mu);
endif;
l=(l-1)/2;
a=teta[1:l];
b=teta[l+1:2*l];
aa=x*a;
bb=x*b;

/* Computation of log-likelihood. */

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];

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

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



/*
** GRADBP
**
** FORMAT    v=gradbp(teta,x,y);
**
** INPUT     teta   :  (2*L+1)x1 vector of parameters.
**           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    v      :  (2*L+1)x1 gradient vector of the log-likelihood function
**                     at the parameter vector.
**
** REMARK    Explicit formulas are used to calculate the derivatives of the
**           cumulative bivariate normal distribution function.
*/


proc gradbp(teta,x,y);

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


/* Initialising some values. */

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


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

l=rows(teta);
if _rhonorm;
   rho=teta[l];
else;
   mu=teta[l];
   rho=transmr(mu);
endif;
l=(l-1)/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);
vrho=zeros(n,4);

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

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

vrho=eerster(aa,bb,rho);
vrho=vrho~(-vrho)~(-vrho)~vrho;
vrho=y.*vrho;

if not _rhonorm;
   vrho=vrho*2*exp(-mu);
   vrho=vrho/((1+exp(-mu))^2);
endif;

zrho=vrho./p;
clear vrho;


/* 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;

drho=sumc(sumc(zrho));


eersteaf=da|db|drho;

retp(eersteaf);
endp;




/*
** HESSBP
**
** FORMAT    v=hessbp(teta,x,y);
**
** INPUT     teta   :  (2*L+1)x1 vector of parameters.
**           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    v      :  (2*L+1)x(2*L+1) hessian matrix of the log-likelihood
**                     function at the parameter vector.
**
** REMARK    Explicit formulas are used to calculate the first and second
**           order derivatives of the
**           cumulative bivariate normal distribution function.
*/




proc hessbp(teta,x,y);

local s1,s2;
local a,b,l,n,aa,bb,p,pp,rho,mu;
local afgx,afgy,afgr,vaa,vbb,vrr,vab,var,vbr,vrr2;
local ddaa,ddbb,ddrr,ddab,ddar,ddbr,tweedeaf;
local i,j,xhulp,hulp;


/* Initialising some values. */

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


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

l=rows(teta);
if _rhonorm;
   rho=teta[l];
else;
   mu=teta[l];
   rho=transmr(mu);
endif;
l=(l-1)/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]=afgx[.,1];
afgx[.,3]=eerstex(aa,bb,rho);
afgx[.,4]=afgx[.,3];

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

afgr=zeros(n,4);
afgr=eerster(aa,bb,rho);
afgr=afgr~afgr~afgr~afgr;

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

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

vrr=zeros(n,4);
hulp=tweederr(aa,bb,rho);
vrr[.,1]=y[.,1].*(p[.,1].*hulp-(afgr[.,1])^2);
vrr[.,2]=y[.,2].*(p[.,2].*(-hulp)-(afgr[.,2])^2);
vrr[.,3]=y[.,3].*(p[.,3].*(-hulp)-(afgr[.,3])^2);
vrr[.,4]=y[.,4].*(p[.,4].*hulp-(afgr[.,4])^2);
vrr=vrr./pp;

if not _rhonorm;
   vrr=vrr*2*exp(-mu);
   vrr2=zeros(n,4);
   hulp=eerster(aa,bb,rho);
   vrr2[.,1]=-y[.,1].*hulp*(1-exp(-2*mu));
   vrr2[.,2]=y[.,2].*hulp*(1-exp(-2*mu));
   vrr2[.,3]=y[.,3].*hulp*(1-exp(-2*mu));
   vrr2[.,4]=-y[.,4].*hulp*(1-exp(-2*mu));
   vrr2=vrr2./p;
   vrr=vrr+vrr2;
   clear vrr2;
   vrr=vrr*2*exp(-mu);
   vrr=vrr/((1+exp(-mu))^4);
endif;


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

var=zeros(n,4);
hulp=tweedexr(aa,bb,rho);
var[.,1]=y[.,1].*(p[.,1].*(-hulp)-afgx[.,1].*afgr[.,1]);
var[.,2]=y[.,2].*(p[.,2].*hulp-afgx[.,2].*afgr[.,2]);
var[.,3]=y[.,3].*(p[.,3].*(-hulp)-afgx[.,3].*afgr[.,3]);
var[.,4]=y[.,4].*(p[.,4].*hulp-afgx[.,4].*afgr[.,4]);
if not _rhonorm;
   var=2*var*exp(-mu)/((1+exp(-mu))^2);
endif;
var=var./pp;

vbr=zeros(n,4);
hulp=tweedeyr(aa,bb,rho);
vbr[.,1]=y[.,1].*(p[.,1].*(-hulp)-afgr[.,1].*afgy[.,1]);
vbr[.,2]=y[.,2].*(p[.,2].*(-hulp)-afgr[.,2].*afgy[.,2]);
vbr[.,3]=y[.,3].*(p[.,3].*hulp-afgr[.,3].*afgy[.,3]);
vbr[.,4]=y[.,4].*(p[.,4].*hulp-afgr[.,4].*afgy[.,4]);
if not _rhonorm;
   vbr=2*vbr*exp(-mu)/((1+exp(-mu))^2);
endif;
vbr=vbr./pp;


clear afgx,afgy,afgr;


/* 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;

ddrr=sumc(sumc(vrr));

clear vrr;

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;

ddar=zeros(l,1);

i=1;
do while i<=l;
   xhulp=x[.,i];
   xhulp=(s2*xhulp')';
   ddar[i]=sumc(sumc(var.*xhulp));
   i=i+1;
endo;

clear var;

ddbr=zeros(l,1);

i=1;
do while i<=l;
   xhulp=x[.,i];
   xhulp=(s1*xhulp')';
   ddbr[i]=sumc(sumc(vbr.*xhulp));
   i=i+1;
endo;

clear vbr;


tweedeaf=(ddaa~ddab~ddar)|(ddab'~ ddbb~ddbr)|(ddar'~ddbr'~ddrr);

retp(tweedeaf);
endp;





/*
** LOGLIKR
**
** FORMAT    v=loglikr(teta,x,y);
**
** INPUT     teta   :  (3*L)x1 vector of parameters.
**           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    v      :  Scalar, calculated log-likelihood.
*/



proc loglikr(teta,x,y);

local a,b,c,l,aa,bb,cc,p,n,result;


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

l=rows(teta);
n=rows(x);

l=l/3;

a=teta[1:l];
b=teta[l+1:2*l];
c=teta[2*l+1:3*l];

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

if not _rhonorm;
   cc=transmr(cc);
endif;

/* Computing the log-likelihood. */

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

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

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



/*
** GRADBPR
**
** FORMAT    v=gradbpr(teta,x,y);
**
** INPUT     teta   :  (3*L)x1 vector of parameters.
**           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    v      :  (3*L)x1 gradient vector of the log-likelihood function at
**                     the parameter vector.
**
** REMARK    Explicit formulas are used to calculate the derivatives of the
**           cumulative bivariate normal distribution function.
*/



proc gradbpr(teta,x,y);

local s1,s2,l,n,nr,k,eersteaf,mu;
local a,b,c,aa,bb,cc,p,za,zb,zc,va,vb,vc;
local xak,xbk,xck,da,db,dc;


/* 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/3;
a=teta[1:l];
b=teta[l+1:2*l];
c=teta[2*l+1:3*l];

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

if not _rhonorm;
   mu=cc;
   cc=transmr(cc);
endif;


cc=cc.*(cc.>=-0.9999 .and cc.<=0.9999)
    -0.9999.*(cc.<-0.9999)+0.9999.*(cc.>0.9999);


/* Matrices to calculate first order derivatives. */

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

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

va[.,1]=y[.,1].*eerstex(-aa,-bb,cc);
va[.,2]=y[.,2].*eerstex(aa,-bb,-cc);
va[.,3]=y[.,3].*eerstex(-aa,bb,-cc);
va[.,4]=y[.,4].*eerstex(aa,bb,cc);

vb[.,1]=y[.,1].*eerstey(-aa,-bb,cc);
vb[.,2]=y[.,2].*eerstey(aa,-bb,-cc);
vb[.,3]=y[.,3].*eerstey(-aa,bb,-cc);
vb[.,4]=y[.,4].*eerstey(aa,bb,cc);

vc[.,1]=y[.,1].*eerster(-aa,-bb,cc);
vc[.,2]=y[.,2].*eerster(aa,-bb,-cc);
vc[.,3]=y[.,3].*eerster(-aa,bb,-cc);
vc[.,4]=y[.,4].*eerster(aa,bb,cc);
if not _rhonorm;
   vc=vc.*(2.*(exp(-mu)./((1+exp(-mu))^2)));
endif;


za=va./p;
zb=vb./p;
zc=vc./p;
clear va,vb,vc,p;


/* Calculating components of gradient vector. */

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

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

eersteaf=da|db|dc;

retp(eersteaf);
endp;





/*
** HESSBPR
**
** FORMAT    v=hessbpr(teta,x,y);
**
** INPUT     teta   :  (3*L)x1 vector of parameters.
**           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    v      :  (3*L)x(3*L) hessian matrix of the log-likelihood function
**                     at the parameter vector.
**
** REMARK    Explicit formulas are used to calculate the first and second
**           order derivatives of the
**           cumulative bivariate normal distribution function.
*/


proc hessbpr(teta,x,y);

local l,n,i,j,p,pp,xhulp;
local s1,s2;
local a,b,c,aa,bb,cc,mu;
local afgx,afgy,afgr,vaa,vbb,vrr,vab,var,vbr,vrr2;
local ddaa,ddbb,ddrr,ddab,ddar,ddbr,tweedeaf;


/* 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/3;
a=teta[1:l];
b=teta[l+1:2*l];
c=teta[2*l+1:3*l];

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

if not _rhonorm;
   mu=cc;
   cc=transmr(cc);
endif;



cc=cc.*(cc.>=-0.9999 .and cc.<=0.9999)
    -0.9999.*(cc.<-0.9999)+0.9999.*(cc.>0.9999);


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

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

pp=p^2;

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

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

afgr=zeros(n,4);
afgr[.,1]=eerster(-aa,-bb,cc);
afgr[.,2]=eerster(aa,-bb,-cc);
afgr[.,3]=eerster(-aa,bb,-cc);
afgr[.,4]=eerster(aa,bb,cc);

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

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

vrr=zeros(n,4);
vrr[.,1]=y[.,1].*(p[.,1].*tweederr(-aa,-bb,cc)-(afgr[.,1])^2);
vrr[.,2]=y[.,2].*(p[.,2].*tweederr(aa,-bb,-cc)-(afgr[.,2])^2);
vrr[.,3]=y[.,3].*(p[.,3].*tweederr(-aa,bb,-cc)-(afgr[.,3])^2);
vrr[.,4]=y[.,4].*(p[.,4].*tweederr(aa,bb,cc) - (afgr[.,4])^2);
vrr=vrr./pp;

if not _rhonorm;
   vrr=vrr.*(2.*exp(-mu));
   vrr2=zeros(n,4);
   vrr2[.,1]=-y[.,1].*eerster(-aa,-bb,cc).*(1-exp(-2*mu));
   vrr2[.,2]=y[.,2].*eerster(aa,-bb,-cc).*(1-exp(-2*mu));
   vrr2[.,3]=y[.,3].*eerster(-aa,bb,-cc).*(1-exp(-2*mu));
   vrr2[.,4]=-y[.,4].*eerster(aa,bb,cc).*(1-exp(-2*mu));
   vrr2=vrr2./p;
   vrr=vrr+vrr2;
   clear vrr2;
   vrr=vrr.*(2.*exp(-mu));
   vrr=vrr./((1+exp(-mu))^4);
endif;

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

var=zeros(n,4);
var[.,1]=y[.,1].*(p[.,1].*tweedexr(-aa,-bb,cc)-afgx[.,1].*afgr[.,1]);
var[.,2]=y[.,2].*(p[.,2].*tweedexr(aa,-bb,-cc)-afgx[.,2].*afgr[.,2]);
var[.,3]=y[.,3].*(p[.,3].*tweedexr(-aa,bb,-cc)-afgx[.,3].*afgr[.,3]);
var[.,4]=y[.,4].*(p[.,4].*tweedexr(aa,bb,cc) - afgx[.,4].*afgr[.,4]);
var=var./pp;
if not _rhonorm;
   var=2*var.*exp(-mu)./((1+exp(-mu))^2);
endif;

vbr=zeros(n,4);
vbr[.,1]=y[.,1].*(p[.,1].*tweedeyr(-aa,-bb,cc)-afgr[.,1].*afgy[.,1]);
vbr[.,2]=y[.,2].*(p[.,2].*tweedeyr(aa,-bb,-cc)-afgr[.,2].*afgy[.,2]);
vbr[.,3]=y[.,3].*(p[.,3].*tweedeyr(-aa,bb,-cc)-afgr[.,3].*afgy[.,3]);
vbr[.,4]=y[.,4].*(p[.,4].*tweedeyr(aa,bb,cc) - afgr[.,4].*afgy[.,4]);
vbr=vbr./pp;
if not _rhonorm;
   vbr=2*vbr.*exp(-mu)./((1+exp(-mu))^2);
endif;

clear afgx,afgy,afgr,pp;


/* 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;


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;


ddrr=zeros(l,l);

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



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;

ddar=zeros(l,l);

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

ddbr=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*xhulp')';
         ddbr[i,j]=sumc(sumc(vbr.*xhulp));
      else;
         ddbr[i,j]=ddbr[j,i];
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;


tweedeaf=(ddaa~ddab~ddar)|(ddab'~ ddbb~ddbr)|(ddar'~ddbr'~ddrr);


retp(tweedeaf);
endp;




