data step1;
  za=1.645;
  nsimul=10000;
  do pattern = 1,2,3;
  put pattern 1-30;
  do kmult=0.1, 1.0, 5.0;
  do margin=log(0.50);
   marg=exp(margin);
   if pattern eq 1 then do; beta1=0.60; beta2=0.65; beta3=0.70; end;
   if pattern eq 2 then do; beta1=0.70; beta2=0.65; beta3=0.60; end;
   if pattern eq 3 then do; beta1=0.70; beta2=0.70; beta3=0.70; end;
   alpha1=-1;
   alpha2=0; 
   alpha3=1; 
  do n=50,100,200;
   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 n*p1r10 le 0 or n*p1r01 le 0 then do;
   v1r=1/(n*p1r10+0.50)+1/(n*p1r01+0.50);
   z1=(log((p1r10+0.50)/(n*p1r01+0.50))-margin)/sqrt(v1r);
   if z1 gt za then sumrej1+1;
end; 
if n*p1r10 gt 0 and n*p1r01 gt 0 then do;
   v1r=1/(n*p1r10)+1/(n*p1r01);
   z1=(log(p1r10/p1r01)-(margin))/sqrt(v1r);
   if z1 gt za then sumrej1+1;
end;
if n*p2r10 le 0 or n*p2r01 le 0 then do;
   v2r=1/(n*p2r10+0.50)+1/(n*p2r01+0.50);
   z2=(log((p2r10+0.50)/(n*p2r01+0.50))-margin)/sqrt(v2r);
   if z2 gt za then sumrej1+1;
end; 
if n*p2r10 gt 0 and n*p2r01 gt 0 then do;
   v2r=1/(n*p2r10)+1/(n*p2r01);
   z2=(log(p2r10/p2r01)-(margin))/sqrt(v2r);
   if z2 gt za then sumrej2+1;
end;
if n*p3r10 le 0 or n*p3r01 le 0 then do;
   v3r=1/(n*p3r10+0.50)+1/(n*p3r01+0.50);
   z3=(log((p3r10+0.50)/(n*p3r01+0.50))-margin)/sqrt(v3r);
   if z3 gt za then sumrej1+1;
end;
if n*p3r10 gt 0 and n*p3r01 gt 0 then do;
   v3r=1/(n*p3r10)+1/(n*p3r01);
   z3=(log(p3r10/p3r01)-(margin))/sqrt(v3r);
   if z3 gt za then sumrej3+1;
end;
   if p1r10 le 0 or p1r01 le 0 then do;
    p1r10=(n*p1r10+0.50)/(n+1);
    p1r01=(n*p1r01+0.50)/(n+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);
   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);
   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; 
   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 z4 gt za then sumrej4+1;
   end; ** end of if p1r10 gt 0 and p2r01 gt 0 and;
end; ** end of do nsimul;
   err1=sumrej1/nsimul;
   err2=sumrej2/nsimul;
   err3=sumrej3/nsimul;
   err4=sumrej4/nsimul;
put kmult 1-5 1 marg 6-10 2 n 11-20 err1 26-30 3 err2 36-40 3 err3 46-50 3 err4 56-60 3;
end;  ** end of do n;
end;  ** end of do lmargin;
end; ** end of kmult;
end;  ** end of pattern;