/*
**                         FILE BPM5.ARC
**
** This file contains the bivariate probit model procedures.
** 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 _rpar ?=0;
declare _restric ?=1;
declare _tsiatis ?=0;
declare _signif ?=0;
declare _imat ?=0;


/*
** BIVPROB
**
** FORMAT    {teta,ll,tetra,llt,teta1,ll1,teta2,ll2}=
**                              bivprob(x,y,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        : Nx4 matrix of contingencies.
**            indvar   : Lx1 or (L-1)x1 character vector; names of independent
**                       variables.
**            depvar   : 2x1 character vector; names of dependent variables.
**
** OUTPUT     teta     : (2*L+1)x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the bivariate parameters.
**            ll       : Bivariate log-likelihood.
**            tetra    : 3x4 matrix, the columns containing
**                       estimates, standard errors, t-values and p-values
**                       for the tetrachoric parameters.
**            llt      : Tetrachoric log-likelihood.
**            teta1    : (L+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+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.
**                       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).
**            _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.
**            _restric   :  MxK matrix of booleans (M equals L or L-1; K equals
**                          2 if _rpar=0 and 3 if _rpar=1) to indicate which
**                          covariates should be taken into account for the
**                          first and second variable (and the correlation)
**                          respectively. By default _restric is scalar 1 (full
**                          model).
**            _tsiatis : 0: Tsiatis' Goodness-of-fit will not be performed.
**                          (default).
**                       1: Tsiatis' Goodness-of-fit will be peformed.
**            _signif  : 0: Significance tests will not be performed.
**                          (default).
**                       1: Significance tests will be peformed.
**            _imat    : 0: (default) Tsiatis' test uses program defined
**                          partition matrix.
**                       1: Tsiatis' test uses user defined partition matrix.
**
** REMARK     Newton-Raphson iteration method for the Bivariate Probit Model.
*/


proc (8)= bivprob(x,y,indvar,depvar);


/* Declaration of variables. */

local hulp,kol1,teta0,behoud,modifica,modif,akkoord,ah,bh,rhoh,rhot,sav,transa,
 l,teta1,teta2,teta3,test,teller,poging,tetra0,ch,tijd1,tijd2,tijd3,tijd4,mask,
 ll,afg1,afg2,infor,xtest,xtest1,bteta,j,a,b,rho,gem,afw,tetashow,mu,transb,y1,
 ct,meanx,stdx,meany,stdy,tval,pval,stde,cases,casesper,naam,scheid,fmt,eps,aa,
 hoofd1,hoofd2,naamvec1,naamvec2,naamvec3,mask1,mask2,fmt1,fmt2,cnt,pvalt,good,
 teta01,ll1,eerste1,tweede1,teta02,ll2,eerste2,tweede2,stdet,tvalt,transc,ssq1,
 teta01s,teta02s,stde1,stde2,tval1,tval2,pval1,pval2,teta0a,ssq,llt,tetra2,bad,
 bb,class,classper,tel,index,maxind,max,p,n,ytetr,tetra1,logl1,y2,naamvec4,alt,
 meanx1,meanx2,meanx3,meanx4,stdx1,stdx2,stdx3,stdx4,mask3,fmt3,eentetr,teta00,
 p1,p2,max1,max2,maxind1,maxind2,class1,class2,classp1,classp2,logl2,bad1,bad2,
 good1,good2,ssq2,trans,transg,gcases,tweetetr,teta00,cardit,tetahulp,la,lb,lc,
 imat;



format /rd 14,4;
disable;


/* Initialising the output variables. */

teta1=0;stde=0;tval=0;pval=0;ll=0;tetra1=0;stdet=0;tvalt=0;
pvalt=0;llt=0;teta01s=0;stde1=0;tval1=0;pval1=0;ll1=0;teta02s=0;
stde2=0;tval2=0;pval2=0;ll2=0;


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


/* Check for faulty data. */

if _rpar/=0 and _rpar/=1;
   errorlog "_rpar must be equal to either 0 or 1.  We change it to 0.";
endif;

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=4);
   errorlog "y does not represent 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;


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

if cols(_restric)/=1 or rows(_restric)/=1;
   if (_rpar==0 and cols(_restric)/=2) or (_rpar==1 and cols(_restric)/=3);
      errorlog "_restric has wrong number of columns. We set it to default.";
      _restric=ones(1,1);
   elseif rows(_restric)/=cols(x) and rows(_restric)/=cols(x)-1;
      errorlog "_restric has wrong number of rows.  We set it to default.";
      _restric=ones(1,1);
   elseif rows(_restric)==cols(x)-1;
      _restric=ones(1,cols(_restric))|_restric;
   elseif not(sumc(_restric)/=0);
      errorlog "Empty models are not allowed. We set _restric to default.";
      _restric=ones(1,1);
   endif;
endif;

if rows(_restric)==1 and cols(_restric)==1;
   _restric=ones(cols(x),2+_rpar);
endif;

la=sumc(_restric[.,1]);
lb=sumc(_restric[.,2]);
if _rpar==1;
   lc=sumc(_restric[.,3]);
endif;


/* 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:4];
   x=hulp[.,5:cols(hulp)];
   clear hulp;
elseif _cases==0;
   _cases=rows(x)|rows(x);
endif;

gcases=_cases;
_cases=0;
/*y=y[.,4]~y[.,3]~y[.,2]~y[.,1];*/


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

l=cols(x)-1;
ct=sumc(y);
meanx=sumc(x[.,2:cols(x)].*sumc(y'))/sumc(ct);
stdx=sqrt(sumc(((x[.,2:cols(x)]-meanx')^2).*sumc(y'))/(sumc(ct)-1));
meany=ct;
meany=((meany[2]+meany[4])|(meany[3]+meany[4]))/sumc(meany);
stdy=zeros(2,1);
stdy[1]=sqrt((ct[1]+ct[3])*(meany[1]^2)+(ct[2]+ct[4])*(1-meany[1])^2);
stdy[2]=sqrt((ct[1]+ct[2])*(meany[2]^2)+(ct[3]+ct[4])*(1-meany[2])^2);
stdy=stdy/sqrt(sumc(ct)-1);
if sumc(y[.,1])>0;
   meanx1=sumc(x[.,2:cols(x)].*y[.,1])/sumc(y[.,1]);
else;
   meanx1=zeros(l,1);
endif;
if sumc(y[.,1])>1;
   stdx1=sqrt(sumc(((x[.,2:cols(x)]-meanx1')^2).*y[.,1])/(sumc(y[.,1])-1));
else;
   stdx1=zeros(l,1);
endif;

if sumc(y[.,2])>0;
   meanx2=sumc(x[.,2:cols(x)].*y[.,2])/sumc(y[.,2]);
else;
   meanx2=zeros(l,1);
endif;
if sumc(y[.,2])>1;
   stdx2=sqrt(sumc(((x[.,2:cols(x)]-meanx2')^2).*y[.,2])/(sumc(y[.,2])-1));
else;
   stdx2=zeros(l,1);
endif;

if sumc(y[.,3])>0;
   meanx3=sumc(x[.,2:cols(x)].*y[.,3])/sumc(y[.,3]);
else;
   meanx3=zeros(l,1);
endif;
if sumc(y[.,3])>1;
   stdx3=sqrt(sumc(((x[.,2:cols(x)]-meanx3')^2).*y[.,3])/(sumc(y[.,3])-1));
else;
   stdx3=zeros(l,1);
endif;

if sumc(y[.,4])>0;
   meanx4=sumc(x[.,2:cols(x)].*y[.,4])/sumc(y[.,4]);
else;
   meanx4=zeros(l,1);
endif;
if sumc(y[.,4])>1;
   stdx4=sqrt(sumc(((x[.,2:cols(x)]-meanx4')^2).*y[.,4])/(sumc(y[.,4])-1));
else;
   stdx4=zeros(l,1);
endif;


/* Number of cases. */

cases=reshape(sumc(y),2,2);
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;
transa=transls(transls(trans,_restric[.,1])',_restric[.,1])';
transb=transls(transls(trans,_restric[.,2])',_restric[.,2])';
if _rpar==1;
   transc=transls(transls(trans,_restric[.,3])',_restric[.,3])';
endif;

transg=transa~zeros(rows(transa),cols(transb))|
zeros(rows(transb),cols(transa))~transb;
if _rpar==0;
   transg=transg~zeros(rows(transg),1)|zeros(1,rows(transg))~1;
elseif _rpar==1;
   transg=transg~zeros(rows(transg),cols(transc))|
   zeros(rows(transc),cols(transg))~transc;
endif;
_transg=transg;


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

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.";
   if not _rhonorm;
      print "The correlation coefficient is transformed.";
      print "You cannot compare the parameters for the correlation";
      print "coefficient to those obtained in case rho is not transformed,";
      print "as, in contrast with rho constant, you now have a different";
      print "model !";
   endif;
   print " ";
   print " ";
endif;
print "************************************************************"\
      "*****************";
print " ";


/* Calling INITVAL for initial estimates. */

tijd1=hsec;

if rows(_teta)==1;
   {teta00,logl1,alt,logl2}=initval(x,y);
   if logl2>logl1;
      print "\g\g";
      teta00=initrs(x,y,2,2);
      teta00=(-teta00[1:rows(teta00)-1])|teta00[rows(teta00)];
      l=l+1;
      teta00=teta00[l+1:2*l 1:l 2*l+1];
      l=l-1;
      logl1=loglik(teta00,x,y);
      alt=teta00;
      logl2=logl1;
   endif;
else;
   teta00=inv(transg)*_teta;
   alt=teta00;
   logl1=loglik(teta00,x,y);
   logl2=logl1;
endif;

wijzigen:
if logl1>=logl2;
   teta0=teta00;
else;
   teta0=alt;
endif;

l=(rows(teta0)-1)/2;

if _rpar==1;
   teta0=teta0|zeros(l-1,1);
   teta00=teta00|zeros(l-1,1);
   alt=alt|zeros(l-1,1);
endif;


if not _rhonorm;
   teta00[2*l+1]=transrm(teta00[2*l+1]);
   alt[2*l+1]=transrm(alt[2*l+1]);
endif;

if _rpar==0;
   if _rhonorm;
      rho=teta0[rows(teta0)];
   else;
      rho=teta0[rows(teta0)];
      mu=transrm(rho);
      teta0[rows(teta0)]=mu;
   endif;
elseif _rpar==1;
   if _rhonorm==0;
      teta0[2*l+1]=transrm(teta0[2*l+1]);
   endif;
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.
*/

ah=teta0[1];
bh=teta0[l+1];
if _rpar==1;
   rhoh=teta0[2*l+1];
endif;

if _initval;

   a=teta0[1:l];
   b=teta0[l+1:2*l];
   a=trans*a;
   b=trans*b;
   if _rpar==1;
      rho=teta0[2*l+1:3*l];
      rho=trans*rho;
   endif;
   tetashow=a|b|rho;


   print "The initial estimates :";
   print tetashow;
   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 a vector of modifications.";
         modif=con(rows(teta0),1);
         tetashow=tetashow+modif;
         a=tetashow[1:l];
         b=tetashow[l+1:2*l];
         a=inv(trans)*a;
         b=inv(trans)*b;
         if _rpar==0;
            rho=tetashow[2*l+1];
            if _rhonorm;
               teta0=a|b|rho;
            else;
               mu=transrm(rho);
               teta0=a|b|mu;
            endif;
         elseif _rpar==1;
            rho=tetashow[2*l+1:3*l];
            rho=inv(trans)*rho;
            teta0=a|b|rho;
         endif;
      else;
         print " ";
         print "Please give in your own vector of initial values.";
         tetashow=con(rows(tetashow),1);
         a=tetashow[1:l];
         b=tetashow[l+1:2*l];
         a=inv(trans)*a;
         b=inv(trans)*b;
         if _rpar==0;
            rho=tetashow[2*l+1];
            if _rhonorm;
               teta0=a|b|rho;
            else;
               mu=transrm(rho);
               teta0=a|b|mu;
            endif;
         elseif _rpar==1;
            rho=tetashow[2*l+1:3*l];
            rho=inv(trans)*rho;
            teta0=a|b|rho;
         endif;
      endif;
   endif;

   if not behoud;
      print " ";
      print "The starting vector now is :";
      print tetashow;
      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;
y1=(y[.,1]+y[.,3])~(y[.,2]+y[.,4]);
y2=(y[.,1]+y[.,2])~(y[.,3]+y[.,4]);


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

hulp=transls(teta0[1:l]',_restric[.,1])';
{teta01,ll1,eerste1,tweede1}=uniprob(hulp,transls(x,_restric[.,1]),y1);
if ll1==0;
   hulp=transls(teta00[1:l]',_restric[.,1])';
   {teta01,ll1,eerste1,tweede1}=uniprob(hulp,transls(x,_restric[.,1]),y1);
endif;
if ll1==0;
   hulp=transls(alt[1:l]',_restric[.,1])';
   {teta01,ll1,eerste1,tweede1}=uniprob(hulp,transls(x,_restric[.,1]),y1);
endif;

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

format /rd 14,4;

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

hulp=transls(teta0[l+1:2*l]',_restric[.,2])';
{teta02,ll2,eerste2,tweede2}=uniprob(hulp,transls(x,_restric[.,2]),y2);
if ll2==0;
   hulp=transls(teta00[l+1:2*l]',_restric[.,2])';
   {teta02,ll2,eerste2,tweede2}=uniprob(hulp,transls(x,_restric[.,2]),y2);
endif;
if ll2==0;
   hulp=transls(alt[l+1:2*l]',_restric[.,2])';
   {teta02,ll2,eerste2,tweede2}=uniprob(hulp,transls(x,_restric[.,2]),y2);
endif;

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

format /rd 14,4;

teta01s=transa*teta01;
teta02s=transb*teta02;

if ll1/=0;
   stde1=sqrt(diag(transa*inv(-tweede1)*transa'));
   tval1=teta01s./stde1;
   pval1=2*(1-cdfn(abs(tval1)));
endif;

if ll2/=0;
   stde2=sqrt(diag(transb*inv(-tweede2)*transb'));
   tval2=teta02s./stde2;
   pval2=2*(1-cdfn(abs(tval2)));
endif;


/* Tetrachoric Probit. */

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

tetra1=zeros(1,3);
tetra1[1]=cdfninv((ytetr[2]+ytetr[4])./sumc(ytetr'));
tetra1[2]=cdfninv((ytetr[3]+ytetr[4])./sumc(ytetr'));
tetra1[3]=teta00[2*l+1];

trap 1;
if not _rhonorm;
   tetra1[3]=transmr(tetra1[3]);
endif;
sav=_rhonorm;
_rhonorm=1;
rhot=rt(tetra1,ytetr);
_rhonorm=sav;
trap 0;
if not scalerr(rhot);
   if not _rhonorm;
      rhot=transrm(rhot);
   endif;
   tetra1=tetra1[1]|tetra1[2]|rhot;
else;
   tetra1=tetra1[1]|tetra1[2]|0;
endif;

tetra2=tetra1;
tetra1=-loglik(tetra1,1,ytetr)~tetra1';
tetra2[3]=0;
tetra2=-loglik(tetra2,1,ytetr)~tetra2';
tetra1=tetra1|tetra2;
tetra1=sortc(tetra1,1);
tetra2=tetra1[2,2:4]';
tetra1=tetra1[1,2:4]';


poging=0;
herhaal:
test=1;
teller=0;
do until (test<=_eps or teller>=_maxit*2.5);
   teller=teller+1;
   eentetr=gradbp(tetra1,1,ytetr);
   tweetetr=hessbp(tetra1,1,ytetr);
   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.";
         goto herhaal;
      else;
         tetra1=0~0~0;
         llt=0;
         eentetr=0;
         tweetetr=0;
         errorlog "Tetrachoric Probit: Divergence again. We quit.";
         goto stoptetr;
      endif;
   endif;
   tetra0=tetra1-infor*eentetr;
   test=sumc(abs(eentetr));
   tetra1=tetra0;
   llt=loglik(tetra1,1,ytetr);
   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,4;
   print "Parameters/gradient at current iteration :";
   if _rhonorm;
      rhot=tetra1[3];
   else;
      mu=tetra1[3];
      rhot=transmr(mu);
   endif;
   print tetra1[1:2]'~rhot|eentetr';
   print " ";
   print "Log-Likelihood at current iteration :" llt;
   print " ";
   print "*****************************************************************"\
   "***************";
   print " ";
endo;

if not _rhonorm;
   hulp=diagrv(eye(3),ones(2,1)|(((1+exp(-mu))^2)/(2*exp(-mu))));
   tweetetr=hulp*tweetetr*hulp;
endif;
stdet=sqrt(diag(inv(-tweetetr)));
tvalt=(tetra1[1 2]|rhot)./stdet;
pvalt=2*(1-cdfn(abs(tvalt)));

if not _rhonorm;
   tetra1[3]=transmr(tetra1[3]);
endif;


stoptetr:


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

if _rpar==0;

   bteta=teta00'|alt'|teta00[1:2*l]'~0|alt[1:2*l]'~0;
   if not(teta0==teta00) and not(teta0==alt);
      bteta=bteta|teta0';
      bteta=bteta|teta0[1:2*l]'~0;
   endif;
   if (ll1/=0) or (ll2/=0);
      bteta=bteta|transsl(teta01'~teta02',_restric)~teta00[2*l+1];
      bteta=bteta|transsl(teta01'~teta02',_restric)~0;
   endif;
   if llt/=0;
      bteta=bteta|tetra1[1]~zeros(1,l-1)~tetra1[2]~zeros(1,l-1)~tetra1[3];
      bteta=bteta|tetra1[1]~zeros(1,l-1)~tetra1[2]~zeros(1,l-1)~0;
   endif;

   bteta=transsl(transls(bteta,vec(_restric)|1),vec(_restric)|1);

   j=1;
   ll=zeros(rows(bteta),1);
   do while j<=rows(bteta);
      ll[j]=-loglik(bteta[j,.]',x,y);
      j=j+1;
   endo;

   bteta=ll~bteta;
   bteta=sortc(bteta,1);
   bteta=bteta[.,2:cols(bteta)]';

elseif _rpar==1;

   bteta=teta00'|alt'|teta00[1:2*l]'~zeros(1,l)|alt[1:2*l]'~zeros(1,l);
   if not(teta0==teta00) and not(teta0==alt);
      bteta=bteta|teta0';
      bteta=bteta|teta0[1:2*l]'~zeros(1,l);
   endif;
   if (ll1/=0) or (ll2/=0);
      bteta=bteta|transsl(teta01'~teta02',_restric[.,1 2])
      ~teta00[2*l+1]~zeros(1,l-1);
      bteta=bteta|transsl(teta01'~teta02',_restric[.,1 2])~zeros(1,l);
   endif;
   if llt/=0;
      bteta=bteta|tetra1[1]~zeros(1,l-1)~tetra1[2]~zeros(1,l-1)~tetra1[3]~
         zeros(1,l-1);
      bteta=bteta|tetra1[1]~zeros(1,l-1)~tetra1[2]~zeros(1,2*l-1);
   endif;

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

   j=1;
   ll=zeros(rows(bteta),1);
   do while j<=rows(bteta);
      ll[j]=-loglikr(bteta[j,.]',x,y);
      j=j+1;
   endo;

   bteta=ll~bteta;
   bteta=sortc(bteta,1);
   bteta=bteta[.,2:cols(bteta)]';

endif;


/* Bivariate Probit. */

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

poging=0;
teta1=bteta[.,1];
if _rpar==0;
   teta1=transls(teta1',vec(_restric)|1)';
elseif _rpar==1;
   teta1=transls(teta1',_restric)';
endif;
herhaalb:
test=1;
teller=1;


do until ((test<eps) or (teller>_maxit));
   trap 1;
   if _rpar==0;
      tetahulp=transsl(teta1',vec(_restric)|1)';
      afg1=gradbp(tetahulp,x,y);
      afg2=hessbp(tetahulp,x,y);
      afg1=transls(afg1',vec(_restric)|1)';
      afg2=transls(transls(afg2,vec(_restric)|1)',vec(_restric)|1)';
   elseif _rpar==1;
      tetahulp=transsl(teta1',_restric)';
      afg1=gradbpr(tetahulp,x,y);
      afg2=hessbpr(tetahulp,x,y);
      afg1=transls(afg1',_restric)';
      afg2=transls(transls(afg2,_restric)',_restric)';
   endif;
   infor=inv(afg2);
   trap 0;
   if scalerr(infor) or scalerr(afg1);
      divbeh:
      poging=poging+1;
      if poging<=cols(bteta)-1;
         teta1=bteta[.,poging+1];
         if _rpar==0;
            teta1=transls(teta1',vec(_restric)|1)';
         elseif _rpar==1;
            teta1=transls(teta1',_restric)';
         endif;
         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;
         {teta1,ll,afg1,afg2}=bpmoncar(x,y);
         _rhonorm=sav;
         if _rpar==0 and _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 _rpar==0;
      if _rhonorm;
         if teta1[rows(teta1)]>=1;
            teta1[rows(teta1)]=0.9999;
         elseif teta1[rows(teta1)]<=-1;
            teta1[rows(teta1)]=-0.9999;
         endif;
      endif;
   endif;
   a=teta1[1:la];
   b=teta1[la+1:la+lb];
   ah=a;
   bh=b;
   a=transa*a;
   b=transb*b;
   if _rpar==0;
       if _rhonorm;
         rho=teta1[la+lb+1];
      else;
         mu=teta1[la+lb+1];
         rho=transmr(mu);
      endif;
   elseif _rpar==1;
      rho=teta1[la+lb+1:la+lb+lc];
      ch=rho;
      rho=transc*rho;
   endif;
   if _itfile or _itscrn;
      cls;
      format /rd 5,0;
      print "BIVARIATE PARAMETERS :  ITERATION NUMBER" teller;
      print "---------------------------------------------";
      print " ";
      format /rd 14,4;
      print "Parameters/gradient at current iteration :";
      print (a|b|rho)'|afg1';
      print " ";
      print "Log-likelihood for current parameters :";
      print " ";
      if _rpar==0;
         ll=loglik(transsl(teta1',vec(_restric)|1)',x,y);
      elseif _rpar==1;
         ll=loglikr(transsl(teta1',_restric)',x,y);
      endif;
      print ll;
      print " ";
      print "************************************************************"\
      "*****************";
      print "  ";
   endif;
   teller=teller+1;
endo;

if teller>=_maxit or ll<llt;
   goto divbeh;
endif;

stopbiv:
cardit=teller-1;
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 _rpar==0 and not _rhonorm;
   hulp=diagrv(eye(la+lb+1),ones(la+lb,1)|(((1+exp(-mu))^2)/(2*exp(-mu))));
   afg2=hulp*afg2*hulp;
endif;
trap 1;
infor=transg*inv(-afg2)*transg';
trap 0;
if not scalerr(infor);
   stde=sqrt(diag(infor));
   tval=(a|b|rho)./stde;
   pval=2*(1-cdfn(abs(tval)));
else;
   hulp=ll;
   ll=0;
endif;


/* Non-standardising the parameters. */

teta3=teta1;
teta2=teta1;
a=teta1[1:la];
b=teta1[la+1:la+lb];
a=transa*a;
b=transb*b;
if _rpar==0;
   if _rhonorm;
      rho=teta1[la+lb+1];
   else;
      mu=teta1[la+lb+1];
      rho=transmr(mu);
   endif;
elseif _rpar==1;
   rho=teta1[la+lb+1:la+lb+lc];
   rho=transc*rho;
endif;
teta1=a|b|rho;

/* Classification table. */

n=rows(x);
aa=x*transsl(ah',_restric[.,1])';
bb=x*transsl(bh',_restric[.,2])';
if _rpar==1;
   rho=x*transsl(ch',_restric[.,3])';
   if not _rhonorm;
      rho=transmr(rho);
   endif;
endif;
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];
max=maxc(p');
maxind=(max .==p[.,1])+2*(max .==p[.,2])+3*(max .==p[.,3])+4*(max .==p[.,4]);
class=zeros(4,4);

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,2);
p1[.,1]=p[.,1]+p[.,3];
p1[.,2]=1-p1[.,1];
max1=maxc(p1');
maxind1=(max1 .==p1[.,1])+2*(max1 .==p1[.,2]);
class1=zeros(2,2);
tel=1;
do while tel<=n;
   index=maxind1[tel];
   class1[index,1]=class1[index,1]+y[tel,1]+y[tel,3];
   class1[index,2]=class1[index,2]+y[tel,2]+y[tel,4];
   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;

p2=zeros(n,2);
p2[.,1]=p[.,1]+p[.,2];
p2[.,2]=1-p2[.,1];
max2=maxc(p2');
maxind2=(max2 .==p2[.,1])+2*(max2 .==p2[.,2]);
class2=zeros(2,2);
tel=1;
do while tel<=n;
   index=maxind2[tel];
   class2[index,1]=class2[index,1]+y[tel,1]+y[tel,2];
   class2[index,2]=class2[index,2]+y[tel,3]+y[tel,4];
   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;


/* 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(((y[.,1]+y[.,3])~(y[.,2]+y[.,4]))./sumc(y')-p1);
ssq1=reshape(ssq1,rows(ssq1)*cols(ssq1),1);
ssq1=sumc(ssq1)~meanc(ssq1)~stdc(ssq1);

ssq2=abs(((y[.,1]+y[.,2])~(y[.,3]+y[.,4]))./sumc(y')-p2);
ssq2=reshape(ssq2,rows(ssq2)*cols(ssq2),1);
ssq2=sumc(ssq2)~meanc(ssq2)~stdc(ssq2);

endif;


/* Final output. */

if rows(_teta)/=1;
   goto geendruk;
endif;

mask3=0~0~1~1~1~1~1~1~1~1;
let fmt3[10,3]=
"*.*s "  7 7
"*.*s"   1 1
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2
"*.*lf " 7 2;

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

let hoofd1="name " "coeff" "std.err" "t-value" "p-value";
mask1=0~1~1~1~1;
let fmt1[5,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lE " 14  4
"*.*lf " 14  4
"*.*lf " 14  4;
naamvec1=indvar;
naamvec3="constant";


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

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

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

print "Means and standard errors of the independent variables per class :";
print "------------------------------------------------------------------";
print " ";
print "  name  |      (0,0)           (1,0)           (0,1)     "\
"      (1,1)";
print "        |  mean    st.e.   mean    st.e.   mean    st.e. "\
"  mean    st.e. ";
print "---------------------------------------------------------"\
"---------------";
naamvec4=naamvec1[2:rows(naamvec1)];
scheid=zeros(rows(naamvec4),1);
scheid=scheid $+ "|";
cnt=printfm(naamvec4~scheid~meanx1~stdx1~meanx2~stdx2~meanx3~stdx3~meanx4~stdx4,
mask3,fmt3);
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 " ";
naam="Class |      0     0 (%)  |      1     1 (%)";
print naam;
print "--------------------------------------------";
mask=0~0~1~1~0~1~1;
let fmt[7,3]=
"*.*s "   3 7
"*.*s "   3 2
"*.*lf "  8 0
"*.*lf "  7 2
"*.*s "   2 2
"*.*lf "  8 0
"*.*lf "  7 2;
let naam="0" "1";
let scheid= "|" "|";
cnt=printfm(naam~scheid~
cases[.,1]~casesper[.,1]~scheid~cases[.,2]~casesper[.,2],mask,fmt);
print " ";
print " ";

if ll/=0;

print "Classification :";
print "----------------";
print " ";
naam="prd\\true |       (0,0)           (1,0)           (0,1)"\
"           (1,1)";
print naam;
print "------------------------------------------------------"\
"--------------------";
scheid=scheid|scheid;
let naam="(0,0)" "(1,0)" "(0,1)" "(1,1)";
mask=0~0~1~1~1~1~1~1~1~1;
let fmt[10,3]=
"*.*s "   7 7
"*.*s "   2 2
"*.*lf "  7 0
"*.*lf "  7 2
"*.*lf "  7 0
"*.*lf "  7 2
"*.*lf "  7 0
"*.*lf "  7 2
"*.*lf "  7 0
"*.*lf "  7 2;
cnt=printfm(naam~scheid~class[.,1]~classper[.,1]~class[.,2]~classper[.,2]~
class[.,3]~classper[.,3]~class[.,4]~classper[.,4],mask,fmt);
print " ";
print "Correct classified :" good "%";
print "Misclassified      :" bad "%";
print " ";
print " ";

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 "-------------------------------------";
format /rd 14,4;
print " ";
naam="prd\\true |         0               1  ";
print naam;
print "------------------------------------------";
scheid=scheid|scheid;
let naam="0" "1";
mask=0~0~1~1~1~1;
let fmt[6,3]=
"*.*s "   7 7
"*.*s "   2 2
"*.*lf "  7 0
"*.*lf "  7 2
"*.*lf "  7 0
"*.*lf "  7 2;
cnt=printfm(naam~scheid[1:2]~class1[.,1]~classp1[.,1]~class1[.,2]~classp1[.,2],
mask,fmt);
print " ";
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 /rd 14,4;
print " ";
print "Total difference   :" ssq1[1];
print "Average difference :" ssq1[2];
print "Standard Error     :" ssq1[3];
print " ";


format /rd 8,8;
print "Marginal classification for :" $depvar[2];
print "-------------------------------------";
format /rd 14,4;
print " ";
naam="prd\\true |         0               1  ";
print naam;
print "------------------------------------------";
scheid=scheid|scheid;
let naam="0" "1";
cnt=printfm(naam~scheid[1:2]~class2[.,1]~classp2[.,1]~class2[.,2]~classp2[.,2],
mask,fmt);
print " ";
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 /rd 14,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;
   print $hoofd1';
   naam=depvar|"corr";
   print "      -------------------------------------------------------------"\
   "-------";
   cnt=printfm(naam~tetra1~stdet~tvalt~pvalt,mask1,fmt1);
   print " ";
   print " ";
   print "Log-likelihood for estimated parameters :";
   print "-----------------------------------------";
   print llt;
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 ll1/=0;
   format /rd 8,8;
   print "Parameters for " $depvar[1] ":";
   print "-------------------------";
   format /rd 14,4;
   print $hoofd1';
   print "      --------------------------------------------------------------"\
   "------";
   cnt=printfm(transls(naamvec1',_restric[.,1])'~teta01s~stde1~tval1~pval1,
   mask1,fmt1);
   print " ";
   print " ";
   print "Log-likelihood for estimated parameters :";
   print "-----------------------------------------";
   print ll1;
   print " ";
   print " ";
endif;

if ll2/=0;
   format /rd 8,8;
   print "Parameters for " $depvar[2] ":";
   print "-------------------------";
   format /rd 14,4;
   print $hoofd1';
   print "      --------------------------------------------------------------"\
   "------";
   cnt=printfm(transls(naamvec1',_restric[.,2])'~teta02s~stde2~tval2~pval2,
   mask1,fmt1);
   print " ";
   print " ";
   print "Log-likelihood for estimated parameters :";
   print "-----------------------------------------";
   print ll2;
   print " ";
   print " ";
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;

format 0,0;
print "Number of iterations: " cardit;
print " ";
print " ";
format /rd 8,8;
print "Parameters for " $depvar[1] ":";
print "-------------------------";
format /rd 14,4;
print $hoofd1';
print "      ----------------------------------------------------------------"\
"----";
cnt=printfm(transls(naamvec1',_restric[.,1])'~teta1[1:la]~stde[1:la]~
tval[1:la]~pval[1:la],mask1,fmt1);
print " ";
print " ";

format /rd 8,8;
print "Parameters for " $depvar[2] ":";
print "-------------------------";
format /rd 14,4;
print $hoofd1';
print "      ----------------------------------------------------------------"\
"----";
cnt=printfm(transls(naamvec1',_restric[.,2])'~teta1[la+1:la+lb]~
stde[la+1:la+lb]~tval[la+1:la+lb]~pval[la+1:la+lb],mask1,fmt1);
print " ";
print " ";

if _rpar==0;
   print "Correlation coefficient :";
   print "-------------------------";
   print $hoofd1';
   print "      --------------------------------------------------------------"\
   "------";
   cnt=printfm(naamvec3~teta1[la+lb+1]~stde[la+lb+1]~tval[la+lb+1]~
   pval[la+lb+1],mask1,fmt1);
   print "  ";
   print "  ";
elseif _rpar==1;
   print "Parameters for the correlation coefficient :";
   print "--------------------------------------------";
   print $hoofd1';
   print "      --------------------------------------------------------------"\
   "------";
   cnt=printfm(transls(naamvec1',_restric[.,3])'~teta1[la+lb+1:la+lb+lc]~
   stde[la+lb+1:la+lb+lc]~tval[la+lb+1:la+lb+lc]~pval[la+lb+1:la+lb+lc],
   mask1,fmt1);
   print "  ";
   print "  ";
endif;

print "Log-likelihood for estimated parameters :";
print "-----------------------------------------";
if _itfile or _itscrn;
   print ll;
else;
   if _rpar==0;
      ll=loglik(teta3,x,y);
   elseif _rpar==1;
      ll=loglikr(teta3,x,y);
   endif;
   print ll;
endif;
print " ";
print " ";

print "Gradient at estimated parameters :";
print "----------------------------------";
print afg1;
print " ";
print " ";

print "Covariance matrix at estimated parameters :";
print "-------------------------------------------";
format /re 14,4;
print infor;
format /rd 14,4;
print " ";
print "************************************************************"\
      "*****************";
print " ";

endif;


if ((ll1/=0 and ll2/=0) or (llt/=0) or (ll/=0 and _rpar==0)) and _signif;

   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 and ll/=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);
      elseif _rpar==1;
         print "    P-value :" cdfchic(abs(2*(ll-ll1-ll2)),lc);
      endif;
      print " ";
   endif;
   if _rpar==0;lc=0;endif;
   if (llt/=0) and ((_rpar==0 and la+lb-2>0) or (_rpar==1 and la+lb+lc-3>0))
   and (ll/=0);
      print "Significance of independent variables :";
      print "---------------------------------------";
      print "Chi-squared :";;
      print 2*(ll-llt);
      if _rpar==0;
         print "    P-value :" cdfchic(abs(2*(ll-llt)),la+lb-2);
      elseif _rpar==1;
         print "    P-value :" cdfchic(abs(2*(ll-llt)),la+lb+lc-3);
      endif;
      print " ";
   endif;
   if (_rpar==0 and ll/=0 and _teta==0);
      print "Does the correlation depend on the covariates ?";
      print "-----------------------------------------------";
      if _rhonorm==1;
         teta3[rows(teta3)]=transrm(teta3[rows(teta3)]);
      endif;
      sav=_rhonorm;
      _rhonorm=0;
      hulp=transsl(teta3',vec(_restric)|1|zeros(l-1,1))';
      afg1=transls(gradbpr(hulp,x,y)',vec(_restric)|ones(l,1))';
      afg2=transls(transls(hessbpr(hulp,x,y),vec(_restric)|ones(l,1))',
      vec(_restric)|ones(l,1))';
      test=afg1'*inv(-afg2)*afg1;
      print "Chi-squared :";;
      print test;
      print "    P-value :" cdfchic(abs(test),l-1);
      print " ";
      _rhonorm=sav;
   endif;
   print "************************************************************"\
         "*****************";
   print " ";

endif;

if (ll1/=0 or ll2/=0 or ll/=0) and rows(x)>2 and _tsiatis;

   print "G O O D N E S S - O F - F I T";
   print "=============================";
   print " ";
   if ll1/=0;
      format 8,8;
      print "Univariate : " $depvar[1] ":";
      print "-----------------------";
      hulp=transls(x,_restric[.,1]);
      if _imat;
         if _file /=2; output off;endif;
         print "Name of the imat-matrix (Number of rows should be ";;
         format 0,0;
         print rows(x);;
         print ") :";
         if _file /=2;output on;endif;
         naam=cons;
         imat=varget(naam);
      else;
         imat=imatmk(x,y,_restric[.,1]);
      endif;
      call tsiatis(teta01,hulp,y1,imat);
   endif;
   if ll2/=0;
      format 8,8;
      print "Univariate : " $depvar[2] ":";
      print "-----------------------";
      hulp=transls(x,_restric[.,2]);
      if _imat;
         if _file/=2;output off;endif;
         print "Name of the imat-matrix (Number of rows should be ";;
         format 0,0;
         print rows(x);;
         print ") :";
         if _file/=2;output on;endif;
         naam=cons;
         imat=varget(naam);
      else;
         imat=imatmk(x,y,_restric[.,2]);
      endif;
      call tsiatis(teta02,hulp,y2,imat);
   endif;
   if ll/=0;
      print "Bivariate :";
      print "-----------";
      if _imat;
         if _file/=2;output off;endif;
         print "Name of the imat-matrix (Number of rows should be ";;
         format 0,0;
         print rows(x);;
         print ") :";
         if _file/=2;output on;endif;
         naam=cons;
         imat=varget(naam);
      else;
         imat=imatmk(x,y,_restric[.,1].or _restric[.,2]);
      endif;
      call tsiatis2(teta2,x,y,imat);
   endif;
   print "************************************************************"\
         "*****************";
   print " ";

endif;

geendruk:

if _file/=2;
   output off;
endif;

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

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

retp(teta1~stde~tval~pval,ll,tetra1~stdet~tvalt~pvalt,llt,
     teta01s~stde1~tval1~pval1,ll1,teta02s~stde2~tval2~pval2,ll2);
endp;




