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




/* Declaration of globals. */

declare _initval ?=1;
declare _itscrn ?=1;
declare _itfile ?=0;
declare _rhonorm ?=1;
declare _file ?=0;
declare _cases ?=0;
declare _teta ?=0;
declare _eps ?=0.001;
declare _maxit ?=30;
declare _transg ?=0;
declare _ordinar ?=1;
declare _rpar ?=0;
declare _restric ?=1;
declare _full ?=1;


/*
** BPROBRS
**
** FORMAT    {teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=
**                              bprobrs(x,y,r,s,indvar,depvar);
**
** 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        : Nx(r*s) matrix of contingencies.
**            r,s      : Dimensions of tables (e.g. 4x5 tables: r=4,s=5).
**            indvar   : Lx1 or (L-1)x1 character vector; names of independent
**                       variables.
**            depvar   : 2x1 character vector; names of dependent variables.
**
** OUTPUT     teta     : (L*(r+s-2)+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the bivariate parameters.
**            ll       : Bivariate log-likelihood.
**            tetra    : (r+s-1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the tetrachoric parameters.
**            llt      : Tetrachoric log-likelihood.
**            teta1    : (L*(r-1)+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the univariate parameters (first variable).
**            ll1      : Univariate log-likelihood (first variable).
**            teta2    : (L*(s-1)+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the univariate parameters (second variable).
**            ll2      : Univariate log-likelihood (second variable).
**
**
** OUTPUT FILE           bprobit.out
**
** GLOBALS    _initval : 0: The user is not asked for an opinion about the
**                          initial values.
**                       1: The user is asked for an opinion about the initial
**                          values (default).
**            _itscrn  : 0: A table with information about the current
**                          iteration does not appear on the screen.
**                       1: A table with information about the current
**                          iteration does appear on the screen (default).
**            _itfile  : 0: A table with information about the current
**                          iteration does not appear in the output file
**                          (default).
**                       1: A table with information about the current
**                          iteration does appear in the output file.
**            _rhonorm   0: Rho is re-parametrized into a variable stretching
**                          from -infty to +infty (inactive for rxs tables).
**                       1: Rho is not re-parametrized (default).
**            _file      0: Output device is file (default).
**                       1: Output device is printer.
**                       2: No output.
**            _eps        : Tolerance for the gradient (default 0.001).
**            _maxit      : Upper bound for the number of iterations
**                          (default 30).
**            _ordinar   1: Standard Probit Model (default).
**                       0: Extended model (generalized dependence on the
**                          bounds).
**              _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.
**
** REMARK     Newton-Raphson iteration method for the Bivariate Probit Model.
**            Initial values are calculated in INITRS. Gradients an Hessians
**            are computed using GRADBPRS and HESSBPRS respectively.
**            When r=2, s=2, you can also use BIVPROB and BPROBIT, as they
**            offer much more facilities that are not available in
**            the general case.
*/


proc (8)= bprobrs(x,y,r,s,indvar,depvar);


/* Declaration of variables. */

local hulp,kol1,teta0,behoud,modifica,modif,akkoord,ah,bh,rhot,sav,afg1,bad2,l,
      teta1,teta2,teta3,test,teller,poging,tetra0,afg1a,afg1b,q,afg2,infor,tal,
      tijd1,tijd2,tijd3,tijd4,bteta,j,i,langr,langs,stoppen,xtest,xtest1,ct,ll,
      a,b,rho,gem,afw,tetashow,mu,naam1,naam2,teta11,teta12,meanx,stdx,tval,p1,
      hoofd2,naamvec2,mask2,fmt2,pval,stde,cases,casesper,naam,good2,bad1,ssq2,
      yt1,yt2,y1,y2,teta01,ll1,eerste1,tweede1,teta02,ll2,eerste2,tweede2,ssq1,
      teta01s,teta02s,stde1,stde2,tval1,tval2,pval1,pval2,teta0a,eps,pvalt,cnt,
      aa,bb,good,bad,class,classper,tel,index,maxind,max,p,n,stdet,tvalt,trans,
      ytetr,tetra1,tetra2,eentetr,tweetetr,ssq,llt,meanx1,stdx1,transg,max1,p2,
      max2,maxind1,maxind2,class1,class2,classp1,classp2,good1,gcases,tetahulp,
      str,hulpa,hulpb,a1,b1,restro,restre,hulpc,rhoh,df;


format /rd 14,-log(_eps)+2;
disable;


/* Creating/resetting the output file. */

if _file==0;
   output file=bprobit.out reset;
elseif _file==1;
   output file="lpt1" reset;
endif;

if _file/=2;
   output on;
endif;


/* Initialising output values. */

teta1=0;
ll=0;
tetra1=0;
llt=0;
teta01s=0;
ll1=0;
teta02s=0;
ll2=0;


/* Check for faulty data. */

if (r<=1) or (s<=1);
   errorlog "r and s should at least be 2";
   end;
elseif (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=r*s);
   errorlog "y does not represent contingency tables.";
   end;
elseif not (y>=0);
   errorlog "y contains negative values.";
   end;
endif;

kol1=x[.,1];
if not(kol1==1);
   errorlog "The first colum of x does not consist of ones, we add one.";
   x=ones(rows(x),1)~x;
   indvar="constant"|indvar;
endif;
if cols(x)/=rows(indvar);
   errorlog "Vector of names of independent variables of the wrong size.";
   end;
elseif rows(depvar)/=2;
   errorlog "Vector of names of dependent variables of the wrong size.";
   end;
endif;


/* Reshaping the restriction matrix into an allowable form. */

l=cols(x);
_restric=vec(_restric);
if rows(_restric)/=1;
   if _rpar==0;
      hulp=2*l+r+s-3;
   else;
      hulp=3*l+r+s-4;
   endif;
   if rows(_restric)/=hulp;
      errorlog "Restriction vector of the wrong shape. We set it to default.";
      _restric=1;
   endif;
endif;

if rows(_restric)==1 and cols(_restric)==1;
   _restric=ones(2*l+r+s-3+_rpar*(l-1),1);
endif;

sav=_ordinar;
_ordinar=1;
restre=transoe(_restric',r,s)';
_ordinar=sav;
restro=transeo(restre',r,s,0)';


/* Dealing with missing values. */

if ismiss(x) or ismiss(y);
   hulp=y~x;
   _cases=rows(hulp);
   hulp=packr(hulp);
   _cases=_cases|rows(hulp);
   y=hulp[.,1:r*s];
   x=hulp[.,r*s+1:cols(hulp)];
   clear hulp;
elseif _cases==0;
   _cases=rows(x)|rows(x);
endif;

gcases=_cases;
_cases=0;


/* Creating marginal tables. */

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

y1=zeros(n,r);
y2=zeros(n,s);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      y1[.,i]=y1[.,i]+y[.,(i-1)*s+j];
      y2[.,j]=y2[.,j]+y[.,(i-1)*s+j];
      j=j+1;
   endo;
   i=i+1;
endo;


/* Means and standard errors of covariates. */

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

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

i=1;
do while i<=r;
   j=1;
   do while j<=s;
      if sumc(y[.,(i-1)*s+j])>0;
         meanx1[.,(i-1)*s+j]=
         sumc(x[.,2:cols(x)].*y[.,(i-1)*s+j])/sumc(y[.,(i-1)*s+j]);
      endif;
      if sumc(y[.,(i-1)*s+j])>1;
         stdx1[.,(i-1)*s+j]=
         sqrt(sumc(((x[.,2:cols(x)]-meanx1[.,(i-1)*s+j]')^2)
         .*y[.,(i-1)*s+j])/(sumc(y[.,(i-1)*s+j])-1));
      endif;
      j=j+1;
   endo;
   i=i+1;
endo;


/* Number of cases. */

cases=reshape(sumc(y),r,s);
casesper=100*(cases/sumc(sumc(cases)));


/*
** 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;
if _rpar==0;
   transg=eye(r+s-2).*.trans;
   transg=transg~zeros(rows(transg),1)|zeros(1,rows(transg))~1;
else;
   transg=eye(r+s-1).*.trans;
endif;
_transg=transg;

/*
transa=transls(transls(trans,_restric[.,1])',_restric[.,1])';
transb=transls(transls(trans,_restric[.,2])',_restric[.,2])';
*/



/* Printing the title and the names of te variables. */

if _full;
cls;
print " ";
print "            BIVARIATE  PROBIT  ESTIMATION  OF  PARAMETERS";
print "            =============================================";
print " ";
print " ";

print "Dependent variables :";
print "---------------------";
print $depvar';
print " ";
print " ";
print "Independent variables :";
print "-----------------------";
print $indvar';
print " ";
print " ";
if _rpar==1;
   print "The correlation depends on the independent variables.";
   print " ";
endif;
print "************************************************************"\
      "*****************";
print " ";
endif;


/* Calling INITRS for initial estimates. */

tijd1=hsec;

if _file/=2;
   output off;
endif;
sav=_rpar;
_rpar=0;
teta0=initrs(x,y,r,s);
_rpar=sav;
if _rpar;
   teta0=teta0|zeros(l-1,1);
endif;
if _file/=2;
   output on;
endif;

wijzigen:


a=teta0[1:l*(r-1)];
a=reshape(a',r-1,l)';
b=teta0[l*(r-1)+1:l*(r+s-2)];
b=reshape(b',s-1,l)';
if _rpar==0;
   rho=teta0[l*(r+s-2)+1];
else;
   rho=teta0[l*(r+s-2)+1:l*(r+s-1)];
endif;
tijd2=hsec;


/*
** If _initval is 1, the user will be asked if he agrees with the initial
** estimates. If not, he will be able to replace or to modify these values.
*/


if _initval;

   a=trans*a;
   b=trans*b;
   if _rpar;
      rho=trans*rho;
   endif;

   print "The initial estimates :";
   print a;
   print b;
   print " ";
   print rho;
   print " ";
   print "Do you want to keep them (1) or not (0) ?";
   behoud=con(1,1);

   if behoud;
      print "OK, we will go on with these values.";
   endif;
   if not behoud;
      print " ";
      print "Do you want to modify them (1) or to replace them (0) ?";
      modifica=con(1,1);
      if modifica;
         print " ";
         print "Please give in your modifications.";
         modif=con(l,r-1);
         a=a+modif;
         modif=con(l,s-1);
         b=b+modif;
         if _rpar==0;
            modif=con(1,1);
         else;
            modif=con(l,1);
         endif;
         rho=rho+modif;
         a=inv(trans)*a;
         b=inv(trans)*b;
         if _rpar;
            rho=inv(trans)*rho;
         endif;
         teta0=vec(a)|vec(b)|rho;
      else;
         print " ";
         print "Please give in your own vector of initial values.";
         a=con(l,r-1);
         b=con(l,s-1);
         if _rpar==0;
            rho=con(1,1);
         else;
            rho=con(l,1);
         endif;
         a=inv(trans)*a;
         b=inv(trans)*b;
         if _rpar;
            rho=inv(trans)*rho;
         endif;
         teta0=vec(a)|vec(b)|rho;
      endif;
   endif;

   if not behoud;
      print " ";
      print "The starting values now are :";
      print trans*a;
      print trans*b;
      print " ";
      if _rpar==0;
         print rho;
      else;
         print trans*rho;
      endif;
      print " ";
      print "Do you agree (1) or not (0) ?";
      akkoord=con(1,1);
      if akkoord;
         print "OK, we will go on with these values.";
      endif;
      if not akkoord;
         goto wijzigen;
      endif;
   endif;

print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";

endif;


if not _itscrn; screen off; endif;
if not _itfile and _file/=2; output off; endif;


/* Univariate Probit. */

tijd3=hsec;

cls;
format /rd 8,8;
print "Determination of Univariate Parameters for " $depvar[1];
print "===================================================";
print " ";

hulp=transls(teta0[1:l*(r-1)]',restre[1:l*(r-1)])';
{teta01,ll1,eerste1,tweede1}=
   uprobrs(hulp,transls(x,restre[1:l]),y1);

if ll1==0;
   errorlog "Unable to determine Univariate Parameters for first variable.";
endif;

format /rd 14,-log(_eps)+2;

cls;
format /rd 8,8;
print "Determination of Univariate Parameters for " $depvar[2];
print "===================================================";
print " ";

hulp=transls(teta0[l*(r-1)+1:l*(r+s-2)]',restre[l*(r-1)+1:l*(r+s-2)])';
{teta02,ll2,eerste2,tweede2}=
   uprobrs(hulp,transls(x,restre[l*(r-1)+1:l*r]),y2);

if ll2==0;
   errorlog "Unable to determine Univariate Parameters for second variable.";
endif;

format /rd 14,-log(_eps)+2;

if _ordinar==1;
   hulp=transsl(teta01',restro[1:l+r-2])';
else;
   hulp=transsl(teta01',restro[1:l*(r-1)])';
endif;
teta01s=trans*reshape(transoe1(hulp',r),r-1,l)';
teta01s=transls(teta01s',restre[1:l])';
if _ordinar==1;
   hulp=transsl(teta02',restro[l+r-1:2*l+r+s-4])';
else;
   hulp=transsl(teta02',restro[l*(r-1)+1:l*(r+s-2)])';
endif;
teta02s=trans*reshape(transoe1(hulp',s),s-1,l)';
teta02s=transls(teta02s',restre[l*(r-1)+1:l*r])';

if _ordinar;
   hulpa=eye(r-1)~(ones(r-1,l-1).*trans[1,2:l]);
   hulpa=hulpa|(zeros(l-1,r-1)~trans[2:l,2:l]);
   hulpa=transls(transls(hulpa,restro[1:l+r-2])',restro[1:l+r-2])';
   hulpb=eye(s-1)~(ones(s-1,l-1).*trans[1,2:l]);
   hulpb=hulpb|(zeros(l-1,s-1)~trans[2:l,2:l]);
   hulpb=transls(transls(hulpb,restro[l+r-1:2*l+r+s-4])',
   restro[l+r-1:2*l+r+s-4])';
endif;

if _ordinar;
   if ll1/=0;
      stde1=sqrt(diag(hulpa*inv(-tweede1)*hulpa'));
      tval1=transeo1(vec(teta01s)',r,0)'./stde1;
      pval1=2*(1-cdfn(abs(tval1)));
   endif;
   if ll2/=0;
      stde2=sqrt(diag(hulpb*inv(-tweede2)*hulpb'));
      tval2=transeo1(vec(teta02s)',s,0)'./stde2;
      pval2=2*(1-cdfn(abs(tval2)));
   endif;
else;
   if ll1/=0;
      hulp=eye(r-1).*.trans;
      hulp=transls(transls(hulp,restre[1:l*(r-1)])',restre[1:l*(r-1)])';
      stde1=sqrt(diag(hulp*inv(-tweede1)*hulp'));
      tval1=vec(teta01s)./stde1;
      pval1=2*(1-cdfn(abs(tval1)));
   endif;
   if ll2/=0;
      hulp=eye(s-1).*.trans;
      hulp=transls(transls(hulp,restre[l*(r-1)+1:l*(r+s-2)])',
      restre[l*(r-1)+1:l*(r+s-2)])';
      stde2=sqrt(diag(hulp*inv(-tweede2)*hulp'));
      tval2=vec(teta02s)./stde2;
      pval2=2*(1-cdfn(abs(tval2)));
   endif;
endif;

teta01s=vec(teta01s);
teta01s=transeo1(transsl(teta01s',restre[1:l*(r-1)]),r,0)';
if _ordinar==1;
   teta01s=transls(teta01s',restro[1:l+r-2])';
else;
   teta01s=transls(teta01s',restro[1:l*(r-1)])';
endif;
teta02s=vec(teta02s);
teta02s=transeo1(transsl(teta02s',restre[l*(r-1)+1:l*(r+s-2)]),s,0)';
if _ordinar==1;
   teta02s=transls(teta02s',restro[l+r-1:2*l+r+s-4])';
else;
   teta02s=transls(teta02s',restro[l*(r-1)+1:l*(r+s-2)])';
endif;


/* Tetrachoric Probit. */

if _full;
cls;
print "Determination of Tetrachoric Correlation.";
print "=========================================";
print " ";
ytetr=sumc(y)';

tetra1=zeros(r+s-1,1);
yt1=zeros(r,1);
yt2=zeros(s,1);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      yt1[i]=yt1[i]+ytetr[(i-1)*s+j];
      yt2[j]=yt2[j]+ytetr[(i-1)*s+j];
      j=j+1;
   endo;
   i=i+1;
endo;

yt1=yt1/sumc(yt1);
yt2=yt2/sumc(yt2);


i=1;
do while i<=r-1;
   tetra1[i]=cdfninv(sumc(yt1[1:i]));
   i=i+1;
endo;
j=1;
do while j<=s-1;
   tetra1[r-1+j]=cdfninv(sumc(yt2[1:j]));
   j=j+1;
endo;
tetra1[r+s-1]=0;

sav=_rpar;
_rpar=0;
poging=0;
herhaal:
test=1;
teller=0;
do until (test<=_eps or teller>=_maxit*2.5);
   teller=teller+1;
   eentetr=gradbprs(tetra1,1,ytetr,r,s);
   tweetetr=hessbprs(tetra1,1,ytetr,r,s);
   trap 1;
   infor=inv(tweetetr);
   if scalerr(infor);
      trap 0;
      poging=poging+1;
      if poging==1;
         /*tetra1=tetra2;*/
         errorlog "Tetrachoric Probit: Divergence. We start again.";
         errorlog "voorlopig is er nog geen tweede alternatief voorzien,";
         errorlog "hiertoe dient RT nog voor rxs tabellen aangepast.";
         goto herhaal;
      else;
         tetra1=0~0~0;
         llt=0;
         eentetr=0;
         tweetetr=0;
         errorlog "Tetrachoric Probit: Divergence again. We quit.";
         goto stoptetr;
      endif;
   endif;
   tetra1=tetra1-infor*eentetr;
   test=sumc(abs(eentetr));
   llt=loglikrs(tetra1,1,ytetr,r,s);
   if _rhonorm;
      if tetra1[rows(tetra1)]>=1;
         tetra1[rows(tetra1)]=0.9;
      elseif tetra1[rows(tetra1)]<=-1;
         tetra1[rows(tetra1)]=-0.9;
      endif;
   endif;
   cls;
   format /rd 5,0;
   print "TETRACHORIC PARAMETERS :  ITERATION NUMBER" teller;
   print "-----------------------------------------------";
   print " ";
   format /rd 14,-log(_eps)+2;
   print "Parameters/gradient at current iteration :";
   print " ";
   if _rhonorm;
      rhot=tetra1[r+s-1];
/*   else;
      mu=tetra1[3];
      rhot=transmr(mu);*/
   endif;
   print tetra1[1:r-1]'|eentetr[1:r-1]';
   print tetra1[r:r+s-2]'|eentetr[r:r+s-2]';
   print rhot|eentetr[r+s-1];
   print " ";
   print "Log-Likelihood at current iteration :" llt;
   print " ";
   print "*****************************************************************"\
   "***************";
   print " ";
endo;
_rpar=sav;


stdet=sqrt(diag(inv(-tweetetr)));
tvalt=(tetra1[1:r+s-2]|rhot)./stdet;
pvalt=2*(1-cdfn(abs(tvalt)));


endif; /*full*/

stoptetr:


/*
** Construction of a set of possible initial values for the bivariate probit
** estimation.
*/

if _rpar==0;
   bteta=teta0';
   if (ll1/=0) or (ll2/=0);
      if _ordinar==1;
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:2*l+r+s-4])~0,
         r,s);
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:2*l+r+s-4])~rho,
         r,s);
      else;
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:l*(r+s-2)])~0,
         r,s);
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:l*(r+s-2)])~rho,
         r,s);
      endif;
   endif;
   if llt/=0;
      bteta=bteta|(vec(tetra1[1:r+s-2]'|zeros(l-1,r+s-2))|tetra1[r+s-1])';
      bteta=bteta|(vec(tetra1[1:r+s-2]'|zeros(l-1,r+s-2))|0)';
   endif;
else;
   bteta=teta0';
   if (ll1/=0) or (ll2/=0);
      if _ordinar==1;
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:2*l+r+s-4])~
         zeros(1,l),r,s);
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:2*l+r+s-4])~rho'
         ,r,s);
      else;
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:l*(r+s-2)])~
         zeros(1,l),r,s);
         bteta=bteta|transoe(transsl(teta01'~teta02',restro[1:l*(r+s-2)])~rho'
         ,r,s);
      endif;
   endif;
   if llt/=0;
      bteta=bteta|(vec(tetra1[1:r+s-2]'|zeros(l-1,r+s-2))|tetra1[r+s-1])'~
      zeros(1,l-1);
      bteta=bteta|(vec(tetra1[1:r+s-2]'|zeros(l-1,r+s-2))|0)'~zeros(1,l-1);
   endif;
endif;

bteta=transsl(transls(bteta,restre),restre);


j=1;
ll=zeros(rows(bteta),1);
do while j<=rows(bteta);
   ll[j]=-loglikrs(bteta[j,.]',x,y,r,s);
   j=j+1;
endo;
bteta=ll~bteta;
bteta=sortc(bteta,1);
bteta=bteta[.,2:cols(bteta)]';


/* Bivariate Probit. */

cls;
print "Determination of Bivariate Parameters.";
print "======================================";
print " ";
eps=_eps;

poging=0;
teta1=bteta[.,1];
tetahulp=teta1;
teta1=transls(transeo(teta1',r,s,0),restro)';
herhaalb:
test=1;
teller=1;

do until ((test<eps) or (teller>_maxit));
   afg1=transls(transeo(gradbprs(tetahulp,x,y,r,s)',r,s,1),restro)';
   afg2=transls(transls(transeo(transeo(hessbprs(tetahulp,x,y,r,s),r,s,1)',
   r,s,1),restro)',restro)';
   trap 1;
   infor=inv(afg2);
   if scalerr(infor);
      trap 0;
      divbeh:
      poging=poging+1;
      if poging<=cols(bteta)-1;
         teta1=bteta[.,poging+1];
         tetahulp=teta1;
         teta1=transls(transeo(teta1',r,s,0),restro)';
         if eps<=0.1;
            eps=eps*10;
         endif;
         errorlog "Bivariate Probit: Divergence. We start again.";
         print "Trial number :" poging+1;
         goto herhaalb;
      else;
         errorlog "Bivariate Probit: All fail. We use a Monte Carlo method.";
         sav=_rhonorm;
         print "Monte Carlo is (nog?) niet aangepast voor rxs tabellen.";
         print "We zien ons dus verplicht van het programma af te breken.";
         goto einde;
         {teta1,ll,afg1,afg2}=bpmoncar(x,y);
         _rhonorm=sav;
         /*if _rhonorm==0;
            teta1[rows(teta1)]=transrm(teta1[rows(teta1)]);
         endif;*/
         goto stopbiv;
      endif;
   endif;

   teta2=teta1-infor*afg1;
   test=sumc(abs(afg1));
   teta1=teta2;
   if _rhonorm and _rpar==0;
      if teta1[rows(teta1)]>=1;
         teta1[rows(teta1)]=0.9999;
      elseif teta1[rows(teta1)]<=-1;
         teta1[rows(teta1)]=-0.9999;
      endif;
   endif;
   tetahulp=transoe(transsl(teta1',restro),r,s)';
   if _itfile or _itscrn;
      cls;
      format /rd 5,0;
      print "BIVARIATE PARAMETERS :  ITERATION NUMBER" teller;
      print "---------------------------------------------";
      format /rd 14,-log(_eps)+2;
      print "Parameters/gradient at current iteration :";
      a=reshape(tetahulp[1:l*(r-1)]',r-1,l)';
      b=reshape(tetahulp[l*(r-1)+1:l*(r+s-2)]',s-1,l)';
/*print "a" a;
print "b" b;*/
      if _rpar==0;
         rho=tetahulp[l*(r+s-2)+1];
      else;
         rho=tetahulp[l*(r+s-2)+1:l*(r+s-1)];
      endif;
      afg1a=transoe(transsl(afg1',restro),r,s)';
      afg1b=transoe(transsl(afg1',restro),r,s)';
      afg1a=reshape(afg1a[1:l*(r-1)],r-1,l)';
      afg1b=reshape(afg1b[l*(r-1)+1:l*(r+s-2)],s-1,l)';
      afg1a=transls(afg1a',restre[1:l])';
      afg1b=transls(afg1b',restre[l*(r-1)+1:l*r])';
      ah=a;
      bh=b;
      rhoh=rho;
      a=trans*a;
      b=trans*b;
      a=transls(a',restre[1:l])';
      b=transls(b',restre[l*(r-1)+1:l*r])';
      if _rpar;
         rho=trans*rho;
      endif;
      print a|afg1a;
      print b|afg1b;
      if _rpar==0;
         print rho|afg1[rows(afg1)];
      else;
         print transls(rho'|afg1[rows(afg1)-l+1:rows(afg1)]',
         restre[l*(r+s-2)+1:l*(r+s-1)]);
      endif;
      print " ";
      print "Log-likelihood for current parameters :";;
      ll=loglikrs(tetahulp,x,y,r,s);
      print ll;
      print "************************************************************"\
      "*****************";
   endif;
   teller=teller+1;
endo;

if teller>=_maxit or (llt/=0 and llt>ll) or (ll1/=0 and ll2/=0 and
ll1+ll2>ll);
   goto divbeh;
endif;

stopbiv:
tijd4=hsec;


if not _itscrn; screen on; endif;
if not _itfile and _file/=2; output on; endif;


/*
** Calculation of standard errors, t-values and p-values for the estimated
** parameters.
*/


if ll/=0;

if  _ordinar;
   trap 1;
   hulp=hulpa~zeros(rows(hulpa),cols(hulpb))|zeros(rows(hulpb),cols(hulpa))
   ~hulpb;
   if _rpar==0;
      hulp=hulp~zeros(rows(hulp),1)|zeros(1,cols(hulp))~1;
   else;
      hulpc=transls(transls(trans,restre[l*(r+s-2)+1:l*(r+s-1)])',
      restre[l*(r+s-2)+1:l*(r+s-1)])';
      hulp=hulp~zeros(rows(hulp),cols(hulpc))|zeros(rows(hulpc),cols(hulp))
      ~hulpc;
   endif;
   infor=hulp*inv(-afg2)*hulp';
   trap 0;
   if not scalerr(infor);
      stde=sqrt(diag(infor));
      if _rpar==1;
         rho=transls(rho',restre[l*(r+s-2)+1:l*(r+s-1)])';
      endif;
      tval=transls(transeo(transsl((vec(a)|vec(b)|rho)',restre),r,s,0),restro)'
      ./stde;
      pval=2*(1-cdfn(abs(tval)));
   else;
      hulp=ll;
      ll=0;
   endif;
else;
   trap 1;
   hulp=transls(transls(transg,restre)',restre)';
   infor=hulp*inv(-afg2)*hulp';
   trap 0;
   if not scalerr(infor);
      stde=sqrt(diag(infor));
      if _rpar==1;
         rho=transls(rho',restre[l*(r+s-2)+1:l*(r+s-1)])';
      endif;
      tval=(vec(a)|vec(b)|rho)./stde;
      pval=2*(1-cdfn(abs(tval)));
   else;
      hulp=ll;
      ll=0;
   endif;
endif;


/* Non-standardising the parameters. */

teta3=vec(ah)|vec(bh)|rhoh;
teta1=vec(a)|vec(b)|rho;
teta1=transls(transeo(transsl(teta1',restre),r,s,0),restro)';

if _full;

/* Classification table. */

aa=x*ah;
bb=x*bh;
if _rpar;
   rho=x*inv(trans)*transsl(rho',restre[l*(r+s-2)+1:l*(r+s-1)])';
endif;


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;

max=maxc(p');
maxind=zeros(n,1);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      maxind=maxind+((i-1)*s+j)*(max .==p[.,(i-1)*s+j]);
      j=j+1;
   endo;
   i=i+1;
endo;

class=zeros(r*s,r*s);

tel=1;
do while tel<=n;
   index=maxind[tel];
   class[index,.]=class[index,.]+y[tel,.];
   tel=tel+1;
endo;

classper=(class./maxc(sumc(class)'|ones(1,cols(class)))')*100;

good=sumc(diag((class/sumc(sumc(class)))*100));
bad=100-good;


/* Marginal classification tables. */


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

if ll1/=0;
max1=maxc(p1');
maxind1=zeros(n,1);
i=1;
do while i<=r;
   maxind1=maxind1+i*(max1 .==p1[.,i]);
   i=i+1;
endo;
class1=zeros(r,r);
tel=1;
do while tel<=n;
   index=maxind1[tel];
   class1[index,.]=class1[index,.]+y1[tel,.];
   tel=tel+1;
endo;
classp1=(class1./maxc(sumc(class1)'|ones(1,cols(class1)))')*100;
good1=sumc(diag((class1/sumc(sumc(class1)))*100));
bad1=100-good1;
endif;

if ll2/=0;
max2=maxc(p2');
maxind2=zeros(n,1);
i=1;
do while i<=s;
   maxind2=maxind2+i*(max2 .==p2[.,i]);
   i=i+1;
endo;
class2=zeros(s,s);
tel=1;
do while tel<=n;
   index=maxind2[tel];
   class2[index,.]=class2[index,.]+y2[tel,.];
   tel=tel+1;
endo;
classp2=(class2./maxc(sumc(class2)'|ones(1,cols(class2)))')*100;
good2=sumc(diag((class2/sumc(sumc(class2)))*100));
bad2=100-good2;
endif;


/* Comparison between predicted and observed frequencies. */

ssq=abs((y./sumc(y'))-p);
ssq=reshape(ssq,rows(ssq)*cols(ssq),1);
ssq=sumc(ssq)~meanc(ssq)~stdc(ssq);

ssq1=abs((y1./sumc(y1'))-p1);
ssq1=reshape(ssq1,rows(ssq1)*cols(ssq1),1);
ssq1=sumc(ssq1)~meanc(ssq1)~stdc(ssq1);

ssq2=abs((y2./sumc(y2'))-p2);
ssq2=reshape(ssq2,rows(ssq2)*cols(ssq2),1);
ssq2=sumc(ssq2)~meanc(ssq2)~stdc(ssq2);

endif; /*_full*/

endif;


/* Final output. */

if _teta/=1 and _full;

let hoofd2="name" "mean" "std.err";
mask2=0~1~1;
let fmt2[3,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4;
naamvec2=indvar[2:rows(indvar)];


print " ";
format /rd 5,0;
print "Computation time :" round((tijd2+tijd4-tijd1-tijd3)/100) "seconds";
format /rd 14,-log(_eps)+2;
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";

print "Means and standard errors of the independent variables :";
print "--------------------------------------------------------";
print $hoofd2';
print "      --------------------------------------";

cnt=printfm(naamvec2~meanx~stdx,mask2,fmt2);
print " ";
print " ";

print "Means and standard errors of the independent variables per class :";
print "------------------------------------------------------------------";
print " ";
print "class            |  ";;
format 8,3;
print $naamvec2';
print "--------------------";;
i=1;
do while i<=l-1;
   print "---------";;
   i=i+1;
endo;
print " ";
langr=(r<10)+2*(r>=10);
langs=(s<10)+2*(s>=10);
langr=maxc(langr|langs);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      format 0,0;
      naam="("$+ftocv(i,langr,0)$+","$+ftocv(j,langr,0)$+")";
      print $naam;;
      if langr==1; print "  ";; endif;
      print "    mean  |  ";;
      format /rd 8,3;
      print meanx1[.,(i-1)*s+j]';
      format 0,0;
      print "          st.e.  |  ";;
      format /rd 8,3;
      print stdx1[.,(i-1)*s+j]';
      j=j+1;
   endo;
   i=i+1;
endo;
print " ";
print " ";
print " ";


print "Number of cases :";
print "-----------------";
print " ";
print "    Total number of cases :" gcases[1];
print "            Cases deleted :" gcases[1]-gcases[2];;
print (100*((gcases[1]-gcases[2])/gcases[1])) " %";
print "Cases kept for processing :" gcases[2];;
print (100*(gcases[2]/gcases[1])) " %";
print " ";

stoppen=0;
tel=0;
do until stoppen;
format 8,0;
print "   Class |";;
   tal=minc((s-tel)|7);
   print seqa(tel+1,1,tal)';
   print "   ------";;
   j=tel+1;
   do while j<=tel+tal;
      print "---------";;
      j=j+1;
   endo;
   print " ";
   i=1;
   do while i<=r;
      format 8,0;
      print  i "|" cases[i,tel+1:tel+tal];
      format 8,2;
      print "       % |" casesper[i,tel+1:tel+tal];
      i=i+1;
   endo;
   print " ";
   tel=tel+7;
   if tel>=s;
      stoppen=1;
   endif;
endo;
print " ";

naam=zeros(r*s,1);
i=1;
do while i<=r;
   j=1;
   do while j<=s;
      naam[(i-1)*s+j]="("$+ftocv(i,langr,0)$+","$+ftocv(j,langr,0)$+")";
      j=j+1;
   endo;
   i=i+1;
endo;

if ll/=0;

print "Classification :";
print "----------------";
tel=0;
stoppen=0;
if langr==1;
   tal=10;
else;
   tal=8;
endif;
do until stoppen;
   if langr==1;
      format 5,0;
   else;
      format 7,0;
   endif;
   print " ";
   if langr==1;
      print "prd\\true |" $naam[tel+1:minc((tel+tal)|r*s)]';
   else;
      print "prd\\true   |" $naam[tel+1:minc((tel+tal)|r*s)]';
   endif;
   if langr==1;
      print "---------";;
   else;
      print "-----------";;
   endif;
   if langr==1;
      str="------";
   else;
      str="--------";
   endif;
   i=1;
   do while i<=minc(tal|(r*s)-tel);
      print str;;
      i=i+1;
   endo;
   print " ";
   i=1;
   do while i<=r;
      j=1;
      do while j<=s;
         if langr==1;
            format 5,0;
         else;
            format 7,0;
         endif;
         if not(class[(i-1)*s+j,tel+1:minc((tel+tal)|r*s)]==0);
            print $naam[(i-1)*s+j] "   |"
            class[(i-1)*s+j,tel+1:minc((tel+tal)|r*s)];
            if langr==1;
               format 5,1;
            else;
               format 7,1;
            endif;
            if langr==1;
               print "percent  |";;
            else;
               print "percent    |";;
            endif;
            classper[(i-1)*s+j,tel+1:minc((tel+tal)|r*s)];
         endif;
         j=j+1;
      endo;
      i=i+1;
   endo;
   if (r*s)<=tel+tal;
      stoppen=1;
   endif;
   tel=tel+tal;
endo;
print " ";
print " ";

format 6,2;
print "Correct classified :" good "%";
print "Misclassified      :" bad "%";
print " ";
print " ";

format 10,4;
print "Comparison between observed and predicted frequencies :";
print "-------------------------------------------------------";
print " ";
print "Total difference   :" ssq[1];
print "Average difference :" ssq[2];
print "Standard Error     :" ssq[3];
print " ";
print " ";

format /rd 8,8;
print "Marginal classification for :" $depvar[1];
print "-------------------------------------";
print " ";
stoppen=0;
tel=0;
do until stoppen;
   format 6,0;
   tal=minc(7|(r-tel));
   print "prd\\true |" seqa(tel+1,1,tal)';
   print "---------";;
   i=1;
   do while i<=tal;
      print "-------";;
      i=i+1;
   endo;
   print " ";
   i=1;
   do while i<=r;
      if not class1[i,.]==0;
         format 6,0;
         print i "  |" class1[i,tel+1:tel+tal];
         format 6,1;
         print "     %   |" classp1[i,tel+1:tel+tal];
      endif;
      i=i+1;
   endo;
   print " ";
   tel=tel+7;
   if tel>=r;
      stoppen=1;
   endif;
endo;
print " ";

format 6,2;
print "Correct classified :" good1 "%";
print "Misclassified      :" bad1 "%";
print " ";
print " ";


format /rd 8,8;
print "Comparison between obs. and pred. frequencies for :" $depvar[1];
print "-----------------------------------------------------------";
format 10,4;
print " ";
print "Total difference   :" ssq1[1];
print "Average difference :" ssq1[2];
print "Standard Error     :" ssq1[3];
print " ";
print " ";

format /rd 8,8;
print "Marginal classification for :" $depvar[2];
print "-------------------------------------";
print " ";
stoppen=0;
tel=0;
do until stoppen;
   format 6,0;
   tal=minc(7|(s-tel));
   print "prd\\true |" seqa(tel+1,1,tal)';
   print "---------";;
   i=1;
   do while i<=tal;
      print "-------";;
      i=i+1;
   endo;
   print " ";
   j=1;
   do while j<=s;
      if not class2[j,.]==0;
         format 6,0;
         print j "  |" class2[j,tel+1:tel+tal];
         format 6,1;
         print "     %   |" classp2[j,tel+1:tel+tal];
      endif;
      j=j+1;
   endo;
   print " ";
   tel=tel+7;
   if tel>=s;
      stoppen=1;
   endif;
endo;
print " ";

format 6,2;
print "Correct classified :" good2 "%";
print "Misclassified      :" bad2 "%";
print " ";
print " ";

format /rd 8,8;
print "Comparison between obs. and pred. frequencies for :" $depvar[2];
print "-----------------------------------------------------------";
format 10,4;
print " ";
print "Total difference   :" ssq2[1];
print "Average difference :" ssq2[2];
print "Standard Error     :" ssq2[3];
print " ";

print "************************************************************"\
      "*****************";
print " ";

endif;



print "T E T R A C H O R I C      P A R A M E T E R S ";
print "==============================================";
print " ";
print " ";

if llt/=0;
   format 10,4;
   print "    Name  Cutpoint |     Coeff   Std.err.    t-value    p-value";
   print " --------------------------------------------------------------";
   naam=depvar[1];
   i=1;
   do while i<=r-2;
      naam=naam|" ";
      i=i+1;
   endo;
   naam=naam|depvar[2];
   j=1;
   do while j<=s-2;
      naam=naam|" ";
      j=j+1;
   endo;
   naam=naam;
   tetra1=tetra1~stdet~tvalt~pvalt;
   i=1;
   do while i<=r-1;
      format 8,8; print  $naam[i];;
      format 6,0; print i "   |";;
      format 10,4; print tetra1[i,.];
      i=i+1;
   endo;
   j=1;
   do while j<=s-1;
      format 8,8; print  $naam[r-1+j];;
      format 6,0; print j "   |";;
      format 10,4; print tetra1[r-1+j,.];
      j=j+1;
   endo;
   print "    Correlation    |" tetra1[r+s-1,.];
   print " ";
   print " ";
   print "Log-likelihood for estimated parameters : ";;
   format 13,4;
   print llt;
   print "-------------------------------------------------------";
else;
   print "Unable to determine Tetrachoric Probit Parameters.";
   print " ";
   print " ";
endif;
print " ";
print " ";


print "************************************************************"\
      "*****************";
print " ";


print "U N I V A R I A T E      P R O B I T     P A R A M E T E R S";
print "============================================================";
print " ";
print " ";


if _ordinar;

   naam1=indvar[2:rows(indvar)];

   if ll1/=0;
      format /rd 8,8;
      print "Parameters for " $depvar[1] ":";
      print "-------------------------";
      print " ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print "---------------------------------------------------------------";
      teta01s=teta01s~stde1~tval1~pval1;

      i=1;
      tel=1;
      do while i<=r-1;
         if restro[i];
            format 7,0; print i "           |";;
            format 10,4; print teta01s[tel,.];
            tel=tel+1;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=l-1;
        if restro[r-1+i];
           format 8,0; print "          " $naam1[i] "|";;
           format 10,4; print teta01s[tel,.];
           tel=tel+1;
        endif;
        i=i+1;
      endo;

      print " ";
      print " ";
      print "Log-likelihood for estimated parameters : ";;
      format 13,4;
      print ll1;
      print "-------------------------------------------------------";
      print " ";
      print " ";
      print " ";
   endif;

   if ll2/=0;
      format /rd 8,8;
      print "Parameters for " $depvar[2] ":";
      print "-------------------------";
      print " ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print "---------------------------------------------------------------";
      teta02s=teta02s~stde2~tval2~pval2;

      i=1;
      tel=1;
      do while i<=s-1;
         if restro[r+l-2+i];
            format 7,0; print i "           |";;
            format 10,4; print teta02s[tel,.];
            tel=tel+1;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=l-1;
        if restro[r+s+l-3+i];
           format 8,0; print "          " $naam1[i] "|";;
           format 10,4; print teta02s[tel,.];
           tel=tel+1;
        endif;
        i=i+1;
      endo;

      print " ";
      print " ";
      print "Log-likelihood for estimated parameters : ";;
      format 13,4;
      print ll2;
      print "-------------------------------------------------------";
      print " ";
      print " ";
      print " ";
   endif;

else;

   if ll1/=0;
      format /rd 8,8;
      print "Parameters for " $depvar[1] ":";
      print "-------------------------";
      print "  ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";
      teta01s=teta01s~stde1~tval1~pval1;

      i=1;
      tel=1;
      do while i<=r-1;
         j=1;
         do while j<=l;
            if restre[(i-1)*l+j];
               format 7,0; print i;;
               format 10,8; print $indvar[j] "|";;
               format 10,4;print teta01s[tel,.];
               tel=tel+1;
            endif;
            j=j+1;
         endo;
         i=i+1;
      endo;

      print " ";
      print " ";
      print "Log-likelihood for estimated parameters : ";;
      format 13,4;
      print ll1;
      print "-------------------------------------------------------";
      print " ";
      print " ";
   endif;


   if ll2/=0;
      format /rd 8,8;
      print "Parameters for " $depvar[2] ":";
      print "-------------------------";
      print "  ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";
      teta02s=teta02s~stde2~tval2~pval2;

      i=1;
      tel=1;
      do while i<=s-1;
         j=1;
         do while j<=l;
            if restre[(r+i-2)*l+j];
               format 7,0; print i;;
               format 10,8; print $indvar[j] "|";;
               format 10,4;print teta02s[tel,.];
               tel=tel+1;
            endif;
            j=j+1;
         endo;
         i=i+1;
      endo;

      print " ";
      print " ";
      print "Log-likelihood for estimated parameters : ";;
      format 13,4;
      print ll2;
      print "-------------------------------------------------------";
      print " ";
      print " ";
   endif;

endif;


if ((ll1==0) and (ll2==0));
   print "Unable to determine Univariate Probit Parameters.";
   print " ";
   print " ";
endif;

print "************************************************************"\
      "*****************";
print " ";


print "B I V A R I A T E     P R O B I T     P A R A M E T E R S";
print "=========================================================";
print " ";
print " ";


if ll/=0;

   if _ordinar;

      teta1=teta1~(stde)~(tval)~(pval);
      naam1=indvar[2:rows(indvar)];
      format /rd 8,8;
      print "Parameters for " $depvar[1] ":";
      print "-------------------------";
      print " ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print "---------------------------------------------------------------";

      i=1;
      tel=1;
      do while i<=r-1;
         if restro[i];
            format 7,0; print i "           |";;
            format 10,4; print teta1[tel,.];
            tel=tel+1;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=l-1;
        if restro[r-1+i];
           format 8,0; print "          " $naam1[i] "|";;
           format 10,4; print teta1[tel,.];
           tel=tel+1;
        endif;
        i=i+1;
      endo;

      print " ";
      print " ";

      format /rd 8,8;
      print "Parameters for " $depvar[2] ":";
      print "-------------------------";
      print " ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print "---------------------------------------------------------------";

      i=1;
      do while i<=s-1;
         if restro[r+l-2+i];
            format 7,0; print i "           |";;
            format 10,4; print teta1[tel,.];
            tel=tel+1;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=l-1;
        if restro[r+s+l-3+i];
           format 8,0; print "          " $naam1[i] "|";;
           format 10,4; print teta1[tel,.];
           tel=tel+1;
        endif;
        i=i+1;
      endo;

      print " ";
      print " ";

   else;

      teta1=teta1~stde~tval~pval;
      format /rd 8,8;
      print "Parameters for " $depvar[1] ":";
      print "-------------------------";
      print "  ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";

      i=1;
      tel=1;
      do while i<=r-1;
         j=1;
         do while j<=l;
            if restre[(i-1)*l+j];
               format 7,0; print i;;
               format 10,8; print $indvar[j] "|";;
               format 10,4;print teta1[tel,.];
               tel=tel+1;
            endif;
            j=j+1;
         endo;
         i=i+1;
      endo;

      print " ";
      print " ";

      format /rd 8,8;
      print "Parameters for " $depvar[2] ":";
      print "-------------------------";
      print "  ";
      format 10,4;
      print "Cutpoint      Name |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";

      i=1;
      do while i<=s-1;
         j=1;
         do while j<=l;
            if restre[(r+i-2)*l+j];
               format 7,0; print i;;
               format 10,8; print $indvar[j] "|";;
               format 10,4;print teta1[tel,.];
               tel=tel+1;
            endif;
            j=j+1;
         endo;
         i=i+1;
      endo;
      print " ";
      print " ";

   endif;


   if _rpar==0;
      print "Correlation coefficient :";
      print "-------------------------";
      print " ";
      print "    Name           |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";
      print "    Correlation    |";;
      format 10,4;
      print teta1[rows(teta1),.];
      print "  ";
      print "  ";
   else;
      print "Parameters for the correlation coefficient :";
      print "--------------------------------------------";
      print " ";
      print "              Name |     Coeff   Std.err.    t-value    p-value";
      print " --------------------------------------------------------------";
      i=1;
      do while i<=l;
         if restre[l*(r+s-2)+i];
            format 7,0; print "        ";;
            format 10,8; print $indvar[i] "|";;
            format 10,4;print teta1[tel,.];
            tel=tel+1;
         endif;
         i=i+1;
      endo;
      print "  ";
      print "  ";
   endif;

   format 13,4;
   print "Log-likelihood for estimated parameters : ";;
   if _itfile or _itscrn;
      print ll;
   else;
      ll=loglik(teta3,x,y);
      print ll;
   endif;
   print "-------------------------------------------------------";
   print " ";
   print " ";

endif;


print "************************************************************"\
      "*****************";
print " ";


if (ll1/=0 and ll2/=0) or (llt/=0);
   print "S I G N I F I C A N C E     T E S T S";
   print "=====================================";
   print " ";
   if ll1/=0 and ll2/=0;
      print "Significance of correlation coefficient :";
      print "-----------------------------------------";
      format /rd 12,4;
      print "Chi-squared :";;
      print 2*(ll-ll1-ll2);
      if _rpar==0;
         print "    P-value :" cdfchic(abs(2*(ll-ll1-ll2)),1);
      else;
         df=sumc(restre[l*(r+s-2)+1:l*(r+s-1)]);
         print "    P-value :" cdfchic(abs(2*(ll-ll1-ll2)),df);
      endif;
      print " ";
   endif;
   if llt/=0;
      print "Significance of independent variables :";
      print "---------------------------------------";
      print "Chi-squared :";;
      print 2*(ll-llt);
      print "    P-value :";;
      print cdfchic(abs(2*(ll-llt)),rows(teta1)-rows(tetra1));
      print " ";
   endif;
   if (_rpar==0 and ll/=0 and _restric/=0);
      _rpar=1;
      tetahulp=tetahulp|zeros(l-1,1);
      afg1=transeo(gradbprs(tetahulp,x,y,r,s)',r,s,1)';
      afg2=transeo(transeo(hessbprs(tetahulp,x,y,r,s),r,s,1)',r,s,1)';
      trap 1;
      test=afg1'*invpd(-afg2)*afg1;
      trap 0;
      if scalerr(test);
         test=afg1'*inv(-afg2)*afg1;
      endif;
      if test>=0;
         print "Does the correlation depend on the covariates ?";
         print "-----------------------------------------------";
         print "Chi-squared :";;
         print test;
         print "    P-value :" cdfchic(abs(test),l-1);
         print " ";
      endif;
      _rpar=0;
   endif;
   print "************************************************************"\
         "*****************";
   print " ";
endif;

endif;

if _file/=2;
   output off;
endif;

format /rd 10,4;
if _file==1;
   lprint "\f";
endif;

_transg=0;
if ll==0;
   ll=hulp;
endif;

einde:
retp(teta1,ll,tetra1,llt,teta01s,ll1,teta02s,ll2);
endp;

