data step1;
  array typeb(jx) error1-error3;
  array fab(jx) fx1-fx3;
  array types(jx) error4-error6;
  array fas(jx) fx4-fx6;
  za=1.645;
  zam=probit(1-0.05/4);
  zm=sqrt(cinv(0.9,3));
  nsimul=10000;
  do kmult=0.1, 1.0, 5.0;
  do margin=log(0.80),log(0.50);
   marg=exp(margin);
   alpha1=-1; beta1=margin;
   alpha2=0;  beta2=margin;
   alpha3=1;  beta3=margin;
  do n=50,100,200;
   if n = 50 then jx=1;
   if n =100 then jx=2;
   if n= 200 then jx=3;
   sumrej1=0; sumflag1=0;
   sumrej2=0; sumflag2=0;
   sumrej3=0; sumflag3=0;
   sumrej4=0; sumflag4=0;
    do i = 1 to nsimul;
   p1r10=0;
   p2r10=0;
   p3r10=0;
   p1r01=0;
   p2r01=0;
   p3r01=0;
   p1r102s10=0;
   p1r102s01=0;
   p1r103s10=0;
   p1r103s01=0;
   p2r103s10=0;
   p2r103s01=0;
   p1r012s10=0;
   p1r012s01=0;
   p1r013s10=0;
   p1r013s01=0;
   p2r013s10=0;
   p2r013s01=0;
   flag1=0;
   flag2=0;
   flag3=0;
   flag4=0;
    do size=1 to n;
   gamma=kmult*rannor(1432151);
   rp1re=exp(alpha1+beta1+gamma)/(1+exp(alpha1+beta1+gamma));
   rp2re=exp(alpha2+beta2+gamma)/(1+exp(alpha2+beta2+gamma));
   rp3re=exp(alpha3+beta3+gamma)/(1+exp(alpha3+beta3+gamma));
   rp1rs=exp(alpha1+gamma)/(1+exp(alpha1+gamma));
   rp2rs=exp(alpha2+gamma)/(1+exp(alpha2+gamma));
   rp3rs=exp(alpha3+gamma)/(1+exp(alpha3+gamma));

   y1=ranbin(1353151,1,rp1re);
   y2=ranbin(1353151,1,rp2re);
   y3=ranbin(1353151,1,rp3re);
   y4=ranbin(1353151,1,rp1rs);
   y5=ranbin(1353151,1,rp2rs);
   y6=ranbin(1353151,1,rp3rs);
 
   if y1 eq 1 and y4 eq 0 then p1r10+1/n;
   if y1 eq 0 and y4 eq 1 then p1r01+1/n;
   if y2 eq 1 and y5 eq 0 then p2r10+1/n;
   if y2 eq 0 and y5 eq 1 then p2r01+1/n;
   if y3 eq 1 and y6 eq 0 then p3r10+1/n;
   if y3 eq 0 and y6 eq 1 then p3r01+1/n;
 
   if y1 eq 1 and y2 eq 1 and y4 eq 0 and y5 eq 0 then p1r102s10+1/n;
   if y1 eq 1 and y2 eq 0 and y4 eq 0 and y5 eq 1 then p1r102s01+1/n;
   if y1 eq 1 and y3 eq 1 and y4 eq 0 and y6 eq 0 then p1r103s10+1/n;
   if y1 eq 1 and y3 eq 0 and y4 eq 0 and y6 eq 1 then p1r103s01+1/n;
   if y2 eq 1 and y3 eq 1 and y5 eq 0 and y6 eq 0 then p2r103s10+1/n;
   if y2 eq 1 and y3 eq 0 and y5 eq 0 and y6 eq 1 then p2r103s01+1/n;
  
   if y1 eq 0 and y2 eq 1 and y4 eq 1 and y5 eq 0 then p1r012s10+1/n;
   if y1 eq 0 and y2 eq 0 and y4 eq 1 and y5 eq 1 then p1r012s01+1/n;
   if y1 eq 0 and y3 eq 1 and y4 eq 1 and y6 eq 0 then p1r013s10+1/n;
   if y1 eq 0 and y3 eq 0 and y4 eq 1 and y6 eq 1 then p1r013s01+1/n;
   if y2 eq 0 and y3 eq 1 and y5 eq 1 and y6 eq 0 then p2r013s10+1/n;
   if y2 eq 0 and y3 eq 0 and y5 eq 1 and y6 eq 1 then p2r013s01+1/n;

end;  ** end of do size;
    if p1r10 le 0 or p1r01 le 0 then do;
     p1r10=(n*p1r10+0.50)/(n+1);
     p1r01=(n*p1r01+0.50)/(n+1);
     flag1=1;
    end;
    if p2r10 le 0 or p2r01 le 0 then do;
     p2r10=(n*p2r10+0.50)/(n+1);
     p2r01=(n*p2r01+0.50)/(n+1);
     flag2=1;
    end;
    if p3r10 le 0 or p3r01 le 0 then do;
     p3r10=(n*p3r10+0.50)/(n+1);
     p3r01=(n*p3r01+0.50)/(n+1);
     flag3=1;
    end;
   if p1r10 gt 0 and p1r01 gt 0 and p2r10 gt 0 and p2r01 gt 0 and 
      p3r10 gt 0 and p3r01 gt 0 then do; 
   if flag1 eq 0 then v1r=1/(n*p1r10)+1/(n*p1r01);
   if flag1 eq 1 then v1r=1/((n+1)*p1r10)+1/((n+1)*p1r01);
   z1=(log(p1r10/p1r01)-(margin))/sqrt(v1r);
   if flag2 eq 0 then v2r=1/(n*p2r10)+1/(n*p2r01);
   if flag2 eq 1 then v2r=1/((n+1)*p2r10)+1/((n+1)*p2r01);
   z2=(log(p2r10/p2r01)-(margin))/sqrt(v2r);
   if flag3 eq 0 then v3r=1/(n*p3r10)+1/(n*p3r01);
   if flag3 eq 1 then v3r=1/((n+1)*p3r10)+1/((n+1)*p3r01);
   z3=(log(p3r10/p3r01)-(margin))/sqrt(v3r);
   c121010=(p1r102s10-p1r10*p2r10)/(n*p1r10*p2r10);
   c121001=(p1r102s01-p1r10*p2r01)/(n*p1r10*p2r01);
   c120110=(p1r012s10-p1r01*p2r10)/(n*p1r01*p2r10);
   c120101=(p1r012s01-p1r01*p2r01)/(n*p1r01*p2r01);
   c12=c121010-c121001-c120110+c120101;
   
   c131010=(p1r103s10-p1r10*p3r10)/(n*p1r10*p3r10);
   c131001=(p1r103s01-p1r10*p3r01)/(n*p1r10*p3r01);
   c130110=(p1r013s10-p1r01*p3r10)/(n*p1r01*p3r10);
   c130101=(p1r013s01-p1r01*p3r01)/(n*p1r01*p3r01);
   c13=c131010-c131001-c130110+c130101;
  
   c231010=(p2r103s10-p2r10*p3r10)/(n*p2r10*p3r10);
   c231001=(p2r103s01-p2r10*p3r01)/(n*p2r10*p3r01);
   c230110=(p2r013s10-p2r01*p3r10)/(n*p2r01*p3r10);
   c230101=(p2r013s01-p2r01*p3r01)/(n*p2r01*p3r01);
   c23=c231010-c231001-c230110+c230101;

  w1=1/3; w2=1/3; w3=1/3;
  lc1=w1*log(p1r10/p1r01)+w2*log(p2r10/p2r01)+w3*log(p3r10/p3r01);
  vlc1=w1**2*v1r+w2**2*v2r+w3**2*v3r+2*w1*w2*c12+2*w1*w3*c13+2*w2*w3*c23;
  if vlc1 gt 0 then z4=(lc1-(margin))/sqrt(vlc1);
   if z1 gt zam or z2 gt zam or z3 gt zam or z4 gt zam then sumrej1+1;
  if z1 gt zm or z2 gt zm or z3 gt zm or z4 gt zm 
  then sumrej2+1; 
   end; ** end of if p1r10 gt 0 and p2r01 gt 0 and;
end; ** end of do nsimul;
   err1=sumrej1/nsimul;
   err2=sumrej2/nsimul;
   typeb=err1;
   types=err2;
 if jx eq 3 then 
put kmult 1-5 1 marg 6-10 2
    error1 15-21 3 error4 24-30 3 error2 33-39 3 error5 42-48 3 error3 51-57 3 error6 60-66 3;
end;  ** end of do n;
end;  ** end of do lmargin;
end; ** end of kmult;