/*                         FILE BPM17.ARC
**
** This file contains the bivariate probit model procedures.
** Part 17.
*/





/* Declaration of globals. */

declare _rhonorm ?=1;
declare _maxit ?=30;
declare _eps ?=0.001;
declare _restric ?=1;
declare _svdtol !=1e-5;
declare _itscrn ?=1;


/*
** IMATMK
**
** FORMAT    imat=imatmk(x,y,restric);
**
** INPUT     x       : Matrix of independent variables.
**           y       : Matrix of dependent variables.
**           restric : Vector of restrictions on the independent variables.
**
** OUTPUT    imat    : Matrix indicating in which part of the partitioned
**                     covariate space a particular covariate vector lies.
**
** REMARK    The purpose of IMAT is to provide with a partition on the
**           covariate space to perfom Tsiatis' Test.
*/


proc imatmk(x,y,restric);

local imat,xhulp,nul,i,cut,klasse;

xhulp=restric.*seqa(1,1,rows(restric));
xhulp=sortc(xhulp,1);
nul=sumc(xhulp.==0);
xhulp=xhulp[nul+1:rows(xhulp)];
xhulp=x[.,xhulp];
cut=(sumc(xhulp[.,2:cols(xhulp)].*sumc(y')))/sumc(sumc(y));
klasse=xhulp[.,2:cols(xhulp)].<=cut';
klasse=(klasse*seqm(1,2,rows(cut)))+1;
imat=zeros(rows(xhulp),2^rows(cut));
i=1;
do while i<=rows(xhulp);
   imat[i,klasse[i]]=1;
   i=i+1;
endo;

retp(imat);
endp;





/*
** TSIATIS
**
** FORMAT    test=tsiatis(teta,x,y,imat);
**
** INPUT     teta    : Parameter vector.
**           x       : NxL matrix of independent variables.
**           y       : Nx2 matrix of dependent variable.
**           imat    : NxK matrix, indicating the partitioning of the
**                     covariate space.
**
** OUTPUT    test    : 1x3 vector of chi-squared value, degrees of freedom
**                     en p-value.
**
** REMARK    TSIATIS performs Tsiatis' Test in the Univariate case.
*/




proc tsiatis(teta,x,y,imat);

local xak,n,l,s,p,aa,d1,q,v1,k,d2,v2,i,j,xhulp,xmat,vmat,amat,bmat,cmat,m,
      df,tmat,test;



n=rows(x);
s=-1~1;
l=rows(teta);
p=zeros(n,2);
aa=x*teta;
m=cols(imat);
xmat=zeros(m,1);
amat=zeros(m,m);
bmat=zeros(m,l);
cmat=zeros(l,l);


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


v1=(y.*q)./p;
v2=-y.*q.*(s.*p.*aa+q)./(p^2);


k=1;
do while (k<=m);
   xak=s.*imat[.,k];
   xmat[k]=sumc(sumc(v1.*xak));
   k=k+1;
endo;


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



i=1;
do while (i<=m);
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,j];
         amat[i,j]=sumc(sumc(v2.*xhulp));
      else;
         amat[i,j]=amat[j,i];
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;


i=1;
do while (i<=m);
   j=1;
   do while (j<=l);
      xhulp=imat[.,i].*x[.,j];
      bmat[i,j]=sumc(sumc(v2.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


amat=-amat;
bmat=-bmat;
cmat=-cmat;


vmat=amat-bmat*inv(cmat)*bmat';

tmat=xmat'*pinv(vmat)*xmat;

df=rank(vmat);
test=cdfchic(tmat,df);

print " ";
format 10,3;
print "       Chi-squared : " tmat;
format 6,0;
print "Degrees of freedom : " df;
format 10,3;
print "    Prob > Chi-sq. : " test;
print " ";

retp(tmat~df~test);
endp;





/*
** TSIATIS2
**
** FORMAT    test=tsiatis2(teta,x,y,imat);
**
** INPUT     x       : NxL matrix of independent variables.
**           y       : Nx4 matrix of dependent variables.
**           imat    : NxK matrix, indicating the partitioning of the
**                     covariate space.
**
** OUTPUT    test    : 1x3 vector of chi-squared value, degrees of freedom
**                     en p-value.
**
** REMARK    TSIATIS2 performs Tsiatis' Test in the Bivariate case, where
**           also a change in the correlation coefficient is considered.
*/




proc tsiatis2(teta,x,y,imat);

local a,b,l,n,aa,bb,p,za,zb,zrho,nr,k,mu,m,xmat,amat,bmat,cmat,tmat,vmat,dfa,
  xbk,da,db,pp,rho,dr,xrk,drho,rho,s1,s2,va,vb,vrho,afgx,afgy,afgr,vaa,vbb,i,
  xak,vrr,vab,var,vbr,vrr2,ddaa,ddbb,ddrr,ddab,ddar,ddbr,xhulp,hulp,tweedeaf,
  j,dfb,dfr,dga,dgb,dgr,dha,dhb,dhr,df,test,sav,mu,vrr2,nul,cut,klasse;



/* Taking _restric into account. */

teta=transsl(teta',vec(_restric)|1)';


/* Initialising some values. */

n=rows(x);
let s1=-1,1,-1,1;
let s2=-1,-1,1,1;
l=rows(teta);
if _rhonorm;
   rho=teta[l];
   mu=transrm(rho);
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;
m=cols(imat);
sav=_rhonorm;
_rhonorm=0;


/* Matrix of probabilities. */

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;




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=y.*vrho;
vrho=vrho*2*exp(-mu);
vrho=vrho/((1+exp(-mu))^2);

zrho=vrho./p;
clear vrho;



/* Calculating components of xmat. */

xak=zeros(m,1);
xbk=zeros(m,1);
xrk=zeros(m,1);
da=zeros(m,1);
db=zeros(m,1);
dr=zeros(m,1);

k=1;
do while (k<=m);
   xak=s1*(imat[.,k]');
   xbk=s2*(imat[.,k]');
   xrk=(s1.*s2)*(imat[.,k]');
   da[k]=sumc(sumc(za.*xak'));
   db[k]=sumc(sumc(zb.*xbk'));
   dr[k]=sumc(sumc(zrho.*xrk'));
   k=k+1;
endo;

clear za,zb;

xmat=da|db|dr;


/* Second order matrices. */

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;
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;
vrr=vrr*2*exp(-mu);
vrr=vrr/((1+exp(-mu))^4);


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]);
var=2*var*exp(-mu)/((1+exp(-mu))^2);
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]);
vbr=2*vbr*exp(-mu)/((1+exp(-mu))^2);
vbr=vbr./pp;


clear afgx,afgy,afgr;


/* Calculation of cmat. */

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=sumc(sumc(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;


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;


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;

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




/* Calculation of amat. */

ddaa=zeros(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if i<=j;
         xhulp=imat[.,i].*imat[.,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(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if i<=j;
         xhulp=imat[.,i].*imat[.,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(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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;



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





/* Calculation of bmat. */


dfa=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=imat[.,i].*x[.,j];
      dfa[i,j]=sumc(sumc(vaa.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dfb=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=(s1.*s2)'.*(imat[.,i].*x[.,j]);
      dfb[i,j]=sumc(sumc(vab.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dga=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=(s1.*s2)'.*(imat[.,i].*x[.,j]);
      dga[i,j]=sumc(sumc(vab.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dgb=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=imat[.,i].*x[.,j];
      dgb[i,j]=sumc(sumc(vbb.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dha=zeros(m,l);

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



dhb=zeros(m,l);

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


dfr=zeros(m,1);

i=1;
do while i<=m;
   xhulp=imat[.,i];
   xhulp=(s2*xhulp')';
   dfr[i]=sumc(sumc(var.*xhulp));
   i=i+1;
endo;


dgr=zeros(m,1);

i=1;
do while i<=m;
   xhulp=imat[.,i];
   xhulp=(s1*xhulp')';
   dgr[i]=sumc(sumc(vbr.*xhulp));
   i=i+1;
endo;



dhr=zeros(m,1);

i=1;
do while i<=m;
   xhulp=imat[.,i];
   dhr[i]=sumc(sumc(vrr.*xhulp));
   i=i+1;
endo;

bmat=dfa~dfb~dfr|dga~dgb~dgr|dha~dhb~dhr;




/* Final calculations. */

amat=-amat;
bmat=-bmat;
cmat=-cmat;

bmat=transls(bmat,vec(_restric)|1);
cmat=transls(transls(cmat,vec(_restric)|1)',vec(_restric)|1)';


vmat=amat-bmat*inv(cmat)*bmat';
tmat=xmat'*pinv(vmat)*xmat;


df=rank(vmat);
test=cdfchic(abs(tmat),df);

print " ";
format 10,3;
print "       Chi-squared : " tmat;
format 6,0;
print "Degrees of freedom : " df;
format 10,3;
print "    Prob > Chi-sq. : " test;
print " ";

_rhonorm=sav;

retp(tmat~df~test);
endp;




/*
** TSIATIS3
**
** FORMAT    test=tsiatis3(teta,x,y,imat);
**
** INPUT     x       : NxL matrix of independent variables.
**           y       : Nx4 matrix of dependent variables.
**           imat    : NxK matrix, indicating the partitioning of the
**                     covariate space.
**
** OUTPUT    test    : 1x3 vector of chi-squared value, degrees of freedom
**                     en p-value.
**
** REMARK    TSIATIS3 performs Tsiatis' Test in the Bivariate case, where
**           the correlation is not considered.
*/




proc tsiatis3(teta,x,y,imat);

local a,b,l,n,aa,bb,p,za,zb,zrho,nr,k,mu,m,xmat,amat,bmat,cmat,tmat,vmat,dfa,
  xbk,da,db,pp,rho,dr,xrk,drho,rho,s1,s2,va,vb,vrho,afgx,afgy,afgr,vaa,vbb,i,
  xak,vrr,vab,var,vbr,vrr2,ddaa,ddbb,ddrr,ddab,ddar,ddbr,xhulp,hulp,tweedeaf,
  j,dfb,dfr,dga,dgb,dgr,dha,dhb,dhr,df,test,sav,mu,vrr2,nul,cut,klasse;



/* Taking _restric into account. */

teta=transsl(teta',vec(_restric)|1)';


/* Initialising some values. */

n=rows(x);
let s1=-1,1,-1,1;
let s2=-1,-1,1,1;
l=rows(teta);
if _rhonorm;
   rho=teta[l];
   mu=transrm(rho);
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;
m=cols(imat);
sav=_rhonorm;
_rhonorm=0;


/* Matrix of probabilities. */

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;

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=y.*vrho;
vrho=vrho*2*exp(-mu);
vrho=vrho/((1+exp(-mu))^2);

zrho=vrho./p;
clear vrho;


/* Calculating components of xmat. */

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

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

clear za,zb;

xmat=da|db;


/* Second order matrices. */

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;
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;
vrr=vrr*2*exp(-mu);
vrr=vrr/((1+exp(-mu))^4);


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]);
var=2*var*exp(-mu)/((1+exp(-mu))^2);
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]);
vbr=2*vbr*exp(-mu)/((1+exp(-mu))^2);
vbr=vbr./pp;


clear afgx,afgy,afgr;


/* Calculation of cmat. */

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=sumc(sumc(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;


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;


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;

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




/* Calculation of amat. */

ddaa=zeros(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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(m,m);

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


ddab=zeros(m,m);

i=1;
do while i<=m;
   j=1;
   do while (j<=m);
      if (i<=j);
         xhulp=imat[.,i].*imat[.,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;

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



/* Calculation of bmat. */


dfa=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=imat[.,i].*x[.,j];
      dfa[i,j]=sumc(sumc(vaa.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dfb=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=(s1.*s2)'.*(imat[.,i].*x[.,j]);
      dfb[i,j]=sumc(sumc(vab.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dga=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=(s1.*s2)'.*(imat[.,i].*x[.,j]);
      dga[i,j]=sumc(sumc(vab.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;


dgb=zeros(m,l);

i=1;
do while i<=m;
   j=1;
   do while (j<=l);
      xhulp=imat[.,i].*x[.,j];
      dgb[i,j]=sumc(sumc(vbb.*xhulp));
      j=j+1;
   endo;
   i=i+1;
endo;

dfr=zeros(m,1);

i=1;
do while i<=m;
   xhulp=imat[.,i];
   xhulp=(s2*xhulp')';
   dfr[i]=sumc(sumc(var.*xhulp));
   i=i+1;
endo;

dgr=zeros(m,1);

i=1;
do while i<=m;
   xhulp=imat[.,i];
   xhulp=(s1*xhulp')';
   dgr[i]=sumc(sumc(vbr.*xhulp));
   i=i+1;
endo;


bmat=dfa~dfb~dfr|dga~dgb~dgr;




/* Final calculations. */

amat=-amat;
bmat=-bmat;
cmat=-cmat;


bmat=transls(bmat,vec(_restric)|1);
cmat=transls(transls(cmat,vec(_restric)|1)',vec(_restric)|1)';


vmat=amat-bmat*inv(cmat)*bmat';
tmat=xmat'*pinv(vmat)*xmat;

df=rank(vmat);
test=cdfchic(abs(tmat),df);

if _itscrn;
   print " ";
   format 10,3;
   print "       Chi-squared : " tmat;
   format 6,0;
   print "Degrees of freedom : " df;
   format 10,3;
   print "    Prob > Chi-sq. : " test;
   print " ";
endif;

_rhonorm=sav;

retp(tmat~df~test);
endp;

