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



/* Declaration of globals. */

declare _rhonorm ?=1;
declare _rpar ?=0;


/*
** LOGLIKRS
**
** FORMAT    v=loglikrs(teta,x,y,r,s);
**
** INPUT     teta   :  L*(r+s-2)+1 or L*(r+s-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      :  Nx(r*s) matrix of contingencies.
**           r,s    :  Dimensions of tables (e.g. 4x5 tables: r=4,s=5).
**
** OUTPUT    v      :  Scalar, calculated log-likelihood.
**
** GLOBALS   _rpar  : 0: The correlation coefficient is a constant to be
**                       estimated (default).
**                    1: The correlation depends on the independent variables
**                       in the same way as the variables do.
*/



proc loglikrs(teta,x,y,r,s);

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


/* Definition of working variables. */

n=rows(x);
l=cols(x);
a=teta[1:l*(r-1)];
a=reshape(a',r-1,l)';
b=teta[l*(r-1)+1:l*(r+s-2)];
b=reshape(b',s-1,l)';
rho=teta[l*(r+s-2)+1:rows(teta)];

aa=x*a;
bb=x*b;
if _rpar;
   rho=x*rho;
endif;


/* Computation of log-likelihood. */

q=zeros(n,(r+1)*(s+1));

i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      q[.,i*(s+1)+j+1]=cdfbvn(aa[.,i],bb[.,j],rho);
      j=j+1;
   endo;
   i=i+1;
endo;

i=1;
do while i<=r-1;
   q[.,(i+1)*(s+1)]=cdfn(aa[.,i]);
   i=i+1;
endo;

j=1;
do while j<=s-1;
   q[.,r*(s+1)+j+1]=cdfn(bb[.,j]);
   j=j+1;
endo;

q[.,(r+1)*(s+1)]=ones(n,1);


p=zeros(n,r*s);

i=1;
do while i<=r;
   j=1;
   do while j<=s;
      p[.,(i-1)*s+j]=q[.,i*(s+1)+j+1]-q[.,i*(s+1)+j]
         -q[.,(i-1)*(s+1)+j+1]+q[.,(i-1)*(s+1)+j];
      j=j+1;
   endo;
   i=i+1;
endo;


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


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





/*
** GRADBPRS
**
** FORMAT    grad=gradbprs(teta,x,y,r,s);
**
** INPUT     teta   :  L*(r+s-2)+1 or L*(r+s-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      :  Nx(r*s) matrix of contingencies.
**           r,s    :  Dimensions of tables (e.g. 4x5 tables: r=4,s=5).
**
** OUTPUT    grad   :  L*(s+r-2)+1 or L*(r+s-1) gradient vector.
**
** GLOBALS   _rpar  : 0: The correlation coefficient is a constant to be
**                       estimated (default).
**                    1: The correlation depends on the independent variables
**                       in the same way as the variables do.
*/



proc  gradbprs(teta,x,y,r,s);


local n,l,q,p,i,j,aa,bb,va,vb,grada,gradb,hulp,pp,k,za,zb,vaa,vbb,vr,qr,gradr,
      a,b,rho;


/* Initialising some values. */

n=rows(x);
l=cols(x);
a=teta[1:l*(r-1)];
a=reshape(a',r-1,l)';
b=teta[l*(r-1)+1:l*(r+s-2)];
b=reshape(b',s-1,l)';
rho=teta[l*(r+s-2)+1:rows(teta)];

aa=x*a;
bb=x*b;
if _rpar;
   rho=x*rho;
endif;


/* Matrices to compute first order derivatives. */

q=zeros(n,(r+1)*(s+1));
qr=zeros(n,(r+1)*(s+1));

i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      q[.,i*(s+1)+j+1]=cdfbvn(aa[.,i],bb[.,j],rho);
      qr[.,i*(s+1)+j+1]=eerster(aa[.,i],bb[.,j],rho);
      j=j+1;
   endo;
   i=i+1;
endo;

i=1;
do while i<=r-1;
   q[.,(i+1)*(s+1)]=cdfn(aa[.,i]);
   i=i+1;
endo;

j=1;
do while j<=s-1;
   q[.,r*(s+1)+j+1]=cdfn(bb[.,j]);
   j=j+1;
endo;

q[.,(r+1)*(s+1)]=ones(n,1);


p=zeros(n,r*s);
va=zeros(n,r*s);
vb=zeros(n,r*s);
vr=zeros(n,r*s);


i=1;
do while i<=r;
   j=1;
   do while j<=s;
      p[.,(i-1)*s+j]=q[.,i*(s+1)+j+1]-q[.,i*(s+1)+j]
         -q[.,(i-1)*(s+1)+j+1]+q[.,(i-1)*(s+1)+j];
      vr[.,(i-1)*s+j]=qr[.,i*(s+1)+j+1]-qr[.,i*(s+1)+j]
         -qr[.,(i-1)*(s+1)+j+1]+qr[.,(i-1)*(s+1)+j];
      j=j+1;
   endo;
   i=i+1;
endo;

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

i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      va[.,(i-1)*s+j]=eerstex(aa[.,i],bb[.,j],rho);
      vb[.,(i-1)*s+j]=eerstey(aa[.,i],bb[.,j],rho);
      j=j+1;
   endo;
   i=i+1;
endo;

i=1;
do while i<=r-1;
   va[.,i*s]=phi1(aa[.,i]);
   i=i+1;
endo;

j=1;
do while j<=s-1;
   vb[.,(r-1)*s+j]=phi1(bb[.,j]);
   j=j+1;
endo;

vaa=va;
i=1;
do while i<=r;
   j=2;
   do while j<=s;
      vaa[.,(i-1)*s+j]=va[.,(i-1)*s+j]-va[.,(i-1)*s+j-1];
      j=j+1;
   endo;
   i=i+1;
endo;
va=vaa;
clear vaa;


vbb=vb;
j=1;
do while j<=s;
   i=2;
   do while i<=r;
      vbb[.,(i-1)*s+j]=vb[.,(i-1)*s+j]-vb[.,(i-2)*s+j];
      i=i+1;
   endo;
   j=j+1;
endo;
vb=vbb;
clear vbb;


za=zeros(n,r*s);
zb=zeros(n,r*s);

i=1;
do while i<=r-1;
   j=1;
   do while j<=s;
      za[.,(i-1)*s+j]=(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
         y[.,i*s+j]./p[.,i*s+j]).*va[.,(i-1)*s+j];
      j=j+1;
   endo;
   i=i+1;
endo;

j=1;
do while j<=s-1;
   i=1;
   do while i<=r;
      zb[.,(i-1)*s+j]=(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
         y[.,(i-1)*s+j+1]./p[.,(i-1)*s+j+1]).*vb[.,(i-1)*s+j];
      i=i+1;
   endo;
   j=j+1;
endo;

vr=vr./p;


/* Computation of gradients. */

grada=zeros(l,r-1);
gradb=zeros(l,s-1);

pp=1;
do while pp<=l;
   i=1;
   do while i<=r-1;
      j=1;
      hulp=0;
      do while j<=s;
         hulp=hulp+za[.,(i-1)*s+j];
         j=j+1;
      endo;
      hulp=hulp.*x[.,pp];
      grada[pp,i]=sumc(hulp);
      i=i+1;
   endo;
   pp=pp+1;
endo;

pp=1;
do while pp<=l;
   j=1;
   do while j<=s-1;
      i=1;
      hulp=0;
      do while i<=r;
         hulp=hulp+zb[.,(i-1)*s+j];
         i=i+1;
      endo;
      hulp=hulp.*x[.,pp];
      gradb[pp,j]=sumc(hulp);
      j=j+1;
   endo;
   pp=pp+1;
endo;

if _rpar==0;
   gradr=sumc(sumc(y.*vr));
else;
   hulp=y.*vr;
   gradr=zeros(l,1);
   i=1;
   do while i<=l;
      gradr[i]=sumc(sumc(x[.,i].*hulp));
      i=i+1;
   endo;
endif;


retp(vec(grada)|vec(gradb)|gradr);
endp;







/*
** HESSBPRS
**
** FORMAT    hess=hessbprs(teta,x,y,r,s);
**
** INPUT     teta   :  L*(r+s-2)+1 or L*(r+s-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      :  Nx(r*s) matrix of contingencies.
**           r,s    :  Dimensions of tables (e.g. 4x5 tables: r=4,s=5).
**
** OUTPUT    hess   :  (L*(s+r-2)+1) x (L*(s+r-2)+1)) or
**                     (L*(s+r-1) x L*(s+r-1)) Hessian matrix.
**
** GLOBALS   _rpar  : 0: The correlation coefficient is a constant to be
**                       estimated (default).
**                    1: The correlation depends on the independent variables
**                       in the same way as the variables do.
*/



proc hessbprs(teta,x,y,r,s);


local n,l,q,p,i,j,aa,bb,va,vb,hulp,hulp2,pp,k,za,zb,vaa,vbb,vr,qr,hessar,
      a,b,rho,vrr,qrr,hessrr,p2,var,vbr,zar,zbr,hessbr,hess,vab,m,
      hulp3,hessaa,hessbb,hessab,i1,qq,j1;


/* Initialising some values. */

n=rows(x);
l=cols(x);
a=teta[1:l*(r-1)];
a=reshape(a',r-1,l)';
b=teta[l*(r-1)+1:l*(r+s-2)];
b=reshape(b',s-1,l)';
rho=teta[l*(r+s-2)+1:rows(teta)];

aa=x*a;
bb=x*b;
if _rpar;
   rho=x*rho;
endif;

if _rpar==0;
   hess=zeros(l*(r+s-2)+1,l*(r+s-2)+1);
else;
   hess=zeros(l*(r+s-1),l*(r+s-1));
endif;


/* Matrices to compute first order derivatives. */

q=zeros(n,(r+1)*(s+1));
qr=zeros(n,(r+1)*(s+1));
qrr=zeros(n,(r+1)*(s+1));

i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      q[.,i*(s+1)+j+1]=cdfbvn(aa[.,i],bb[.,j],rho);
      qr[.,i*(s+1)+j+1]=eerster(aa[.,i],bb[.,j],rho);
      qrr[.,i*(s+1)+j+1]=tweederr(aa[.,i],bb[.,j],rho);
      j=j+1;
   endo;
   i=i+1;
endo;

i=1;
do while i<=r-1;
   q[.,(i+1)*(s+1)]=cdfn(aa[.,i]);
   i=i+1;
endo;

j=1;
do while j<=s-1;
   q[.,r*(s+1)+j+1]=cdfn(bb[.,j]);
   j=j+1;
endo;

q[.,(r+1)*(s+1)]=ones(n,1);


p=zeros(n,r*s);
va=zeros(n,r*s);
vb=zeros(n,r*s);
vr=zeros(n,r*s);
vaa=zeros(n,r*s);
vbb=zeros(n,r*s);
vab=zeros(n,r*s);
var=zeros(n,r*s);
vbr=zeros(n,r*s);
vrr=zeros(n,r*s);


i=1;
do while i<=r;
   j=1;
   do while j<=s;
      p[.,(i-1)*s+j]=q[.,i*(s+1)+j+1]-q[.,i*(s+1)+j]
         -q[.,(i-1)*(s+1)+j+1]+q[.,(i-1)*(s+1)+j];
      vr[.,(i-1)*s+j]=qr[.,i*(s+1)+j+1]-qr[.,i*(s+1)+j]
         -qr[.,(i-1)*(s+1)+j+1]+qr[.,(i-1)*(s+1)+j];
      vrr[.,(i-1)*s+j]=qrr[.,i*(s+1)+j+1]-qrr[.,i*(s+1)+j]
         -qrr[.,(i-1)*(s+1)+j+1]+qrr[.,(i-1)*(s+1)+j];
      j=j+1;
   endo;
   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);


i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      va[.,(i-1)*s+j]=eerstex(aa[.,i],bb[.,j],rho);
      vb[.,(i-1)*s+j]=eerstey(aa[.,i],bb[.,j],rho);
      vaa[.,(i-1)*s+j]=tweedexx(aa[.,i],bb[.,j],rho);
      vbb[.,(i-1)*s+j]=tweedeyy(aa[.,i],bb[.,j],rho);
      vab[.,(i-1)*s+j]=tweedexy(aa[.,i],bb[.,j],rho);
      var[.,(i-1)*s+j]=tweedexr(aa[.,i],bb[.,j],rho);
      vbr[.,(i-1)*s+j]=tweedeyr(aa[.,i],bb[.,j],rho);
      j=j+1;
   endo;
   i=i+1;
endo;


i=1;
do while i<=r-1;
   va[.,i*s]=phi1(aa[.,i]);
   vaa[.,i*s]=-aa[.,i].*va[.,i*s];
   i=i+1;
endo;

j=1;
do while j<=s-1;
   vb[.,(r-1)*s+j]=phi1(bb[.,j]);
   vbb[.,(r-1)*s+j]=-bb[.,j].*vb[.,(r-1)*s+j];
   j=j+1;
endo;

hulp=va;
hulp2=var;
hulp3=vaa;
i=1;
do while i<=r;
   j=2;
   do while j<=s;
      hulp[.,(i-1)*s+j]=va[.,(i-1)*s+j]-va[.,(i-1)*s+j-1];
      hulp2[.,(i-1)*s+j]=var[.,(i-1)*s+j]-var[.,(i-1)*s+j-1];
      hulp3[.,(i-1)*s+j]=vaa[.,(i-1)*s+j]-vaa[.,(i-1)*s+j-1];
      j=j+1;
   endo;
   i=i+1;
endo;
va=hulp;
var=hulp2;
vaa=hulp3;


hulp=vb;
hulp2=vbr;
hulp3=vbb;
j=1;
do while j<=s;
   i=2;
   do while i<=r;
      hulp[.,(i-1)*s+j]=vb[.,(i-1)*s+j]-vb[.,(i-2)*s+j];
      hulp2[.,(i-1)*s+j]=vbr[.,(i-1)*s+j]-vbr[.,(i-2)*s+j];
      hulp3[.,(i-1)*s+j]=vbb[.,(i-1)*s+j]-vbb[.,(i-2)*s+j];
      i=i+1;
   endo;
   j=j+1;
endo;
vb=hulp;
vbr=hulp2;
vbb=hulp3;
clear hulp hulp2 hulp3;


za=zeros(n,r*s);
zb=zeros(n,r*s);
zar=zeros(n,r*s);
zbr=zeros(n,r*s);

i=1;
do while i<=r-1;
   j=1;
   do while j<=s;
      za[.,(i-1)*s+j]=(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
         y[.,i*s+j]./p[.,i*s+j]).*va[.,(i-1)*s+j];
      zar[.,(i-1)*s+j]=-(y[.,(i-1)*s+j].*vr[.,(i-1)*s+j]./p2[.,(i-1)*s+j]-
         y[.,i*s+j].*vr[.,i*s+j]./p2[.,i*s+j]).*va[.,(i-1)*s+j]
         +(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-y[.,i*s+j]./p[.,i*s+j])
         .*var[.,(i-1)*s+j];
      j=j+1;
   endo;
   i=i+1;
endo;


j=1;
do while j<=s-1;
   i=1;
   do while i<=r;
      zb[.,(i-1)*s+j]=(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
       y[.,(i-1)*s+j+1]./p[.,(i-1)*s+j+1]).*vb[.,(i-1)*s+j];
      zbr[.,(i-1)*s+j]=-(y[.,(i-1)*s+j].*vr[.,(i-1)*s+j]./p2[.,(i-1)*s+j]-
       y[.,(i-1)*s+j+1].*vr[.,(i-1)*s+j+1]./p2[.,(i-1)*s+j+1]).*vb[.,(i-1)*s+j]
       +(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
       y[.,(i-1)*s+j+1]./p[.,(i-1)*s+j+1]).*vbr[.,(i-1)*s+j];
      i=i+1;
   endo;
   j=j+1;
endo;


/* Computing parts of Hessian. */

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

i=1;
do while i<=r-1;
   i1=1;
   do while i1<=r-1;
      hulp=0;
      if i==i1;
         j=1;
         do while j<=s;
            hulp=hulp+(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-y[.,i*s+j]./p[.,i*s+j]).*
            vaa[.,(i-1)*s+j]
            -(y[.,(i-1)*s+j]./p2[.,(i-1)*s+j]+y[.,i*s+j]./p2[.,i*s+j]).*
            va[.,(i-1)*s+j].*va[.,(i-1)*s+j];
            j=j+1;
         endo;
      elseif i-1==i1;
         j=1;
         do while j<=s;
            hulp=hulp+y[.,(i-1)*s+j]./p2[.,(i-1)*s+j].*va[.,(i-1)*s+j].*
            va[.,(i-2)*s+j];
            j=j+1;
         endo;
      elseif i+1==i1;
         j=1;
         do while j<=s;
            hulp=hulp+y[.,i*s+j]./p2[.,i*s+j].*va[.,(i-1)*s+j].*
            va[.,i*s+j];
            j=j+1;
         endo;
      endif;
      if abs(i-i1)<=1;
         pp=1;
         do while pp<=l;
            qq=1;
            do while qq<=l;
               hessaa[(i-1)*l+pp,(i1-1)*l+qq]=
               sumc(hulp.*x[.,pp].*x[.,qq]);
               qq=qq+1;
            endo;
            pp=pp+1;
         endo;
      endif;
      i1=i1+1;
   endo;
   i=i+1;
endo;


hessbb=zeros(l*(s-1),l*(s-1));

j=1;
do while j<=s-1;
   j1=1;
   do while j1<=s-1;
      hulp=0;
      if j==j1;
         i=1;
         do while i<=r;
            hulp=hulp+(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-
            y[.,(i-1)*s+j+1]./p[.,(i-1)*s+j+1]).*vbb[.,(i-1)*s+j]
            -(y[.,(i-1)*s+j]./p2[.,(i-1)*s+j]+
            y[.,(i-1)*s+j+1]./p2[.,(i-1)*s+j+1]).*
            vb[.,(i-1)*s+j].*vb[.,(i-1)*s+j];
            i=i+1;
         endo;
      elseif j-1==j1;
         i=1;
         do while i<=r;
            hulp=hulp+y[.,(i-1)*s+j]./p2[.,(i-1)*s+j].*vb[.,(i-1)*s+j].*
            vb[.,(i-1)*s+j-1];
            i=i+1;
         endo;
      elseif j+1==j1;
         i=1;
         do while i<=r;
            hulp=hulp+y[.,(i-1)*s+j+1]./p2[.,(i-1)*s+j+1].*vb[.,(i-1)*s+j].*
            vb[.,(i-1)*s+j+1];
            i=i+1;
         endo;
      endif;
      if abs(j-j1)<=1;
         pp=1;
         do while pp<=l;
            qq=1;
            do while qq<=l;
               hessbb[(j-1)*l+pp,(j1-1)*l+qq]=
               sumc(hulp.*x[.,pp].*x[.,qq]);
               qq=qq+1;
            endo;
            pp=pp+1;
         endo;
      endif;
      j1=j1+1;
   endo;
   j=j+1;
endo;


hessab=zeros(l*(r-1),l*(s-1));

i=1;
do while i<=r-1;
   j=1;
   do while j<=s-1;
      hulp=(y[.,(i-1)*s+j]./p[.,(i-1)*s+j]-y[.,i*s+j]./p[.,i*s+j]
      -y[.,(i-1)*s+j+1]./p[.,(i-1)*s+j+1]+y[.,i*s+j+1]./p[.,i*s+j+1])
      .*vab[.,(i-1)*s+j]
      -y[.,(i-1)*s+j]./p2[.,(i-1)*s+j].*va[.,(i-1)*s+j].*vb[.,(i-1)*s+j]
      +y[.,(i-1)*s+j+1]./p2[.,(i-1)*s+j+1].*va[.,(i-1)*s+j+1].*vb[.,(i-1)*s+j]
      +y[.,i*s+j]./p2[.,i*s+j].*va[.,(i-1)*s+j].*vb[.,i*s+j]
      -y[.,i*s+j+1]./p2[.,i*s+j+1].*va[.,(i-1)*s+j+1].*vb[.,i*s+j];
      pp=1;
      do while pp<=l;
         qq=1;
         do while qq<=l;
            hessab[(i-1)*l+pp,(j-1)*l+qq]=
            sumc(hulp.*x[.,pp].*x[.,qq]);
            qq=qq+1;
         endo;
         pp=pp+1;
      endo;
      j=j+1;
   endo;
   i=i+1;
endo;


if _rpar==0;
   hessar=zeros(l,r-1);
   i=1;
   do while i<=r-1;
      j=1;
      hulp=0;
      do while j<=s;
         hulp=hulp+zar[.,(i-1)*s+j];
         j=j+1;
      endo;
      pp=1;
      do while pp<=l;
         hessar[pp,i]=sumc(hulp.*x[.,pp]);
         pp=pp+1;
      endo;
      i=i+1;
   endo;
else;
   hessar=zeros(l*(r-1),l);
   i=1;
   do while i<=r-1;
      j=1;
      hulp=0;
      do while j<=s;
         hulp=hulp+zar[.,(i-1)*s+j];
         j=j+1;
      endo;
      pp=1;
      do while pp<=l;
         m=1;
         do while m<=l;
            hessar[(i-1)*l+pp,m]=sumc(hulp.*x[.,pp].*x[.,m]);
            m=m+1;
         endo;
         pp=pp+1;
      endo;
      i=i+1;
   endo;
endif;


if _rpar==0;
   hessbr=zeros(l,s-1);
   j=1;
   do while j<=s-1;
      i=1;
      hulp=0;
      do while i<=r;
         hulp=hulp+zbr[.,(i-1)*s+j];
         i=i+1;
      endo;
      pp=1;
      do while pp<=l;
         hessbr[pp,j]=sumc(hulp.*x[.,pp]);
         pp=pp+1;
      endo;
      j=j+1;
   endo;
else;
   hessbr=zeros(l*(s-1),l);
   j=1;
   do while j<=s-1;
      i=1;
      hulp=0;
      do while i<=r;
         hulp=hulp+zbr[.,(i-1)*s+j];
         i=i+1;
      endo;
      pp=1;
      do while pp<=l;
         m=1;
         do while m<=l;
            hessbr[(j-1)*l+pp,m]=sumc(hulp.*x[.,pp].*x[.,m]);
            m=m+1;
         endo;
         pp=pp+1;
      endo;
      j=j+1;
   endo;
endif;


if _rpar==0;
   hessrr=(p.*vrr-vr.*vr)./p2;
   hessrr=sumc(sumc(y.*hessrr));
else;
   hessrr=zeros(l,l);
   hulp=y.*(p.*vrr-vr.*vr)./p2;
   i=1;
   do while i<=l;
      j=1;
      do while j<=l;
         if i<=j;
            hessrr[i,j]=sumc(sumc(x[.,i].*x[.,j].*hulp));
         else;
            hessrr[i,j]=hessrr[j,i];
         endif;
         j=j+1;
      endo;
      i=i+1;
   endo;
endif;


/* Patching parts of Hessian together. */

hess[1:l*(r-1),1:l*(r-1)]=hessaa;
hess[l*(r-1)+1:l*(r+s-2),l*(r-1)+1:l*(r+s-2)]=hessbb;
hess[1:l*(r-1),l*(r-1)+1:l*(r+s-2)]=hessab;
hess[l*(r-1)+1:l*(r+s-2),1:l*(r-1)]=hessab';
if _rpar==0;
   hess[1:l*(r+s-2),l*(r+s-2)+1]=vec(hessar)|vec(hessbr);
   hess[l*(r+s-2)+1,1:l*(r+s-2)]=(vec(hessar)|vec(hessbr))';
   hess[rows(hess),cols(hess)]=hessrr;
else;
   hess[1:l*(r+s-2),l*(r+s-2)+1:l*(r+s-1)]=hessar|hessbr;
   hess[l*(r+s-2)+1:l*(r+s-1),1:l*(r+s-2)]=hessar'~hessbr';
   hess[l*(r+s-2)+1:l*(r+s-1),l*(r+s-2)+1:l*(r+s-1)]=hessrr;
endif;


retp(hess);
endp;

