data step;
array sapow(k) sfinp1-sfinp3;
array an(k) finn1-finn3;
za=1.96;
zb=0.842;
desired=0.80;
nsimul=100000;
do mean=0.50,1.0,3;
do rho=0.10,0.50,0.80;
beta=rho/(1-rho);
alpha=mean/beta;
do eta=log(0.50),log(1.2),log(1.5);
rmt=exp(eta);
do gamma=log(0.90),log(1.0),log(1.1);
rmp=exp(gamma);
a=1+exp(eta+gamma);
b=exp(eta)+exp(gamma);
p1=exp(eta+gamma)/(1+exp(eta+gamma));
p2=exp(gamma)/(exp(gamma)+exp(eta));
pp=(exp(eta+gamma)+exp(gamma))/(1+exp(eta+gamma)+exp(eta)+exp(gamma));
v0=(1/(mean*a)+1/(mean*b))/(4*pp*(1-pp));
va=(1/(mean*a*p1*(1-p1))+1/(mean*b*p2*(1-p2)))/4;
n=ceil((za*sqrt(v0)+zb*sqrt(va))**2/eta**2);
sumrej=0;
do sim=1 to nsimul;
sx11=0;
sx21=0;
sx12=0;
sx22=0;
stot1=0;
stot2=0;
do num=1 to n;
mu=beta*rangam(1433151,alpha);
x11=ranpoi(1452151,mu);
sx11+x11;
x21=ranpoi(1452141,mu*exp(eta+gamma));
sx21+x21;
stot1+x11+x21;
mu=beta*rangam(1433151,alpha);
x12=ranpoi(1452151,mu*exp(eta));
sx12+x12;
x22=ranpoi(1452151,mu*exp(gamma));
sx22+x22;
stot2+x12+x22;
end; ** end of do num=1 to n;
if sx11 eq 0 or sx21 eq 0 or sx12 eq 0 or sx22 eq 0 then do;
sx11=sx11+0.50; sx21=sx21+0.50; sx12=sx12+0.50; sx22=sx22+0.50;
stot1=stot1+1; stot2=stot2+1;
end;
if sx11 gt 0 and sx21 gt 0 and sx12 gt 0 and sx22 gt 0 then do;
vartm=(1/stot1+1/stot2)/(4*(sx21+sx22)/(stot1+stot2)*(stot1+stot2-(sx21+sx22))/(stot1+stot2));
vart=(1/sx11+1/sx21+1/sx12+1/sx22)/4;
ztest=0.5*log((sx21*sx12)/(sx22*sx11));
z=ztest/sqrt(vartm);
if z gt 1.96 or z lt -1.96 then sumrej+1;
end; ** end of xdot1 gt 0 and xdot2 gt 0;
end; ** end of do sim=1 to nsimul;
spower=sumrej/nsimul;
if gamma=log(0.90) then do; k = 1; an=n; sapow=spower; end;
if gamma=log(1) then do; k = 2; an=n; sapow=spower; end;
if gamma=log(1.1) then do; k = 3; an=n; sapow=spower; end;
if k eq 3 then do;
put mean 1-5 1 rho 6-10 2 rmt 11-15 1 finn1 16-20 @21 "(" sfinp1 22-26 3 @27 ")"
finn2 36-40 @41 "(" sfinp2 42-46 3 @47 ")"
finn3 56-60 @61 "(" sfinp3 62-66 3 @67 ")";
end; ** end of if k eq 3;
end; ** end of do gamma;
end; ** end of do eta;
end; ** end of do rho;
end; ** end of do mean;