data step1;

** This is the simulation program for comparing various interval estimators;

** of the gamma correlation in Table 2 under stratified sampling;

array obp1(i) ob1-ob12;

array obp2(i) ob13-ob24;

array obp3(i) ob25-ob36;

array obp4(i) ob37-ob48;

array obp5(i) ob49-ob60;

array obp6(i) ob61-ob72;

array obp7(i) ob73-ob84;

array obp8(i) ob85-ob96;

array obp9(i) ob97-ob108;

array obp10(i) ob109-ob120;

array orow(j) obp1-obp10;

array tforow(i) tfob1-tfob12;

array torow(i) tob1-tob12;

array n(j) size1-size10;

nsimul=10000;

za=1.96;

do gamma=-0.50,0.0,0.50;

do numstra=2,5,10;

do strasize=50,100,200,500;

sumcov1=0; sumcov2=0; sumcov3=0; sumcov4=0; sumcov5=0; sumcov6=0; sumcov7=0;

sumlen1=0; sumlen2=0; sumlen3=0; sumlen4=0; sumlen5=0; sumlen6=0; sumlen7=0;

sflag1=0; sflag2=0; sflag3=0; sflag4=0; sflag5=0; sflag6=0; sflag7=0;

do simul=1 to nsimul;

flag1=0; flag2=0; flag3=0; flag4=0; flag5=0; flag6=0; flag7=0;

do i = 1 to 12;

tforow=0;

end;

do j= 1 to numstra;

do i=1 to 12;

orow=0;

end; ** do i = 1;

end; ** do j=1;

sumwtwls=0; ** sum of weights for WLS;

summhnum=0; ** sum of MH estimator in the numberator;

sumwlspt=0; ** sum og WLS point estimator in the numerator;

summhden=0; ** sum of MH estimator in the denominator;

sumwtr=0; ** sum of weight for wls of the log point estimator of the ratio;

sumwtrpt=0; ** sum of the wls of the log point estimator of the ratio;

a2=0; b2=0; c2=0; tn=0;

do s=1 to numstra;

which=rantbl(1432171,0.25,0.25,0.25,0.25);

if gamma eq -0.5 and which eq 1 then do;

p11=0.01; p12=0.02; p13=0.03; p14=0.14;

p21=0.05; p22=0.02; p23=0.05; p24=0.10;

p31=0.14; p32=0.01; p33=0.383; p34=0.047; end;

if gamma eq -0.5 and which eq 2 then do;

p11=0.02; p12=0.05; p13=0.05; p14=0.20;

p21=0.02; p22=0.01; p23=0.02; p24=0.05;

p31=0.20; p32=0.02; p33=0.248; p34=0.112; end;

if gamma eq -0.5 and which eq 3 then do;

p11=0.03; p12=0.03; p13=0.03; p14=0.11;

p21=0.03; p22=0.03; p23=0.05; p24=0.15;

p31=0.25; p32=0.02; p33=0.192; p34=0.078; end;

if gamma eq -0.5 and which eq 4 then do;

p11=0.01; p12=0.05; p13=0.07; p14=0.32;

p21=0.02; p22=0.05; p23=0.02; p24=0.05;

p31=0.10; p32=0.02; p33=0.184; p34=0.106; end;

if gamma eq 0.0 and which eq 1 then do;

p11=0.05; p12=0.05; p13=0.05; p14=0.05;

p21=0.07; p22=0.07; p23=0.07; p24=0.07;

p31=0.13; p32=0.13; p33=0.13; p34=0.13; end;

if gamma eq 0.0 and which eq 2 then do;

p11=0.10; p12=0.05; p13=0.05; p14=0.10;

p21=0.15; p22=0.05; p23=0.05; p24=0.15;

p31=0.05; p32=0.05; p33=0.179; p34=0.021; end;

if gamma eq 0.0 and which eq 3 then do;

p11=0.10; p12=0.05; p13=0.15; p14=0.03;

p21=0.10; p22=0.05; p23=0.15; p24=0.03;

p31=0.10; p32=0.05; p33=0.167; p34=0.023; end;

if gamma eq 0.0 and which eq 4 then do;

p11=0.10; p12=0.05; p13=0.07; p14=0.32;

p21=0.02; p22=0.05; p23=0.02; p24=0.05;

p31=0.01; p32=0.05; p33=0.092; p34=0.168; end;

if gamma eq 0.5 and which eq 1 then do;

p11=0.10; p12=0.05; p13=0.05; p14=0.05;

p21=0.05; p22=0.20; p23=0.03; p24=0.02;

p31=0.05; p32=0.02; p33=0.20; p34=0.18; end;

if gamma eq 0.5 and which eq 2 then do;

p11=0.15; p12=0.03; p13=0.02; p14=0.01;

p21=0.05; p22=0.05; p23=0.10; p24=0.20;

p31=0.01; p32=0.02; p33=0.22; p34=0.14; end;

if gamma eq 0.5 and which eq 3 then do;

p11=0.20; p12=0.020; p13=0.02; p14=0.02;

p21=0.03; p22=0.10; p23=0.10; p24=0.10;

p31=0.02; p32=0.05; p33=0.326; p34=0.014; end;

if gamma eq 0.5 and which eq 4 then do;

p11=0.10; p12=0.05; p13=0.03; p14=0.02;

p21=0.03; p22=0.05; p23=0.10; p24=0.15;

p31=0.02; p32=0.05; p33=0.166; p34=0.234; end;

nx=ranpoi(1542151,strasize);

do while (nx le 5);

nx=ranpoi(1542151,strasize);

end; ** end of do while;

tn+nx;

j=s;

n=nx;

do numpat=1 to n;

x=rantbl(1432151,p11,p12,p13,p14,p21,p22,p23,p24,p31,p32,p33,p34);

i=x;

orow+1/n;

tforow+1;

end; ** end of do numpat=1;

p11c=0; p12c=0; p13c=0; p14c=0; p21c=0; p22c=0; p23c=0; p24c=0; p31c=0; p32c=0; p33c=0; p34c=0;

do i = 6,7,8,10,11,12;

p11c+orow; end;

do i = 7,8,11,12;

p12c+orow; end;

do i=8,12;

p13c+orow; end;

do i=10,11,12;

p21c+orow; end;

do i =1,11,12;

p22c+orow; end;

do i=1,2,12;

p23c+orow; end;

do i=1,2,3;

p24c+orow; end;

do i = 1, 5;

p32c+orow; end;

do i=1, 2, 5,6;

p33c+orow; end;

do i=1,2,3,5,6,7;

p34c+orow; end;

i=1; ep11=orow;

i=2; ep12=orow;

i=3; ep13=orow;

i=4; ep14=orow;

i=5; ep21=orow;

i=6; ep22=orow;

i=7; ep23=orow;

i=8; ep24=orow;

i=9; ep31=orow;

i=10; ep32=orow;

i=11; ep33=orow;

i=12; ep34=orow;

pic=ep11*p11c+ep12*p12c+ep13*p13c+ep21*p21c+ep22*p22c+ep23*p23c+ep24*p24c+ep32*p32c+ep33*p33c+ep34*p34c;

p11d=0; p12d=0; p13d=0; p14d=0; p21d=0; p22d=0; p23d=0; p24d=0; p31d=0; p32d=0; p33d=0; p34d=0;

do i = 5, 9;

p12d+orow; end;

do i = 5, 6, 9, 10;

p13d+orow; end;

do i = 5,6,7,9,10,11;

p14d+orow; end;

do i = 2, 3, 4;

p21d+orow; end;

do i = 3, 4, 9;

p22d+orow; end;

do i = 4, 9, 10;

p23d+orow; end;

do i = 9, 10, 11;

p24d+orow; end;

do i = 2,3,4,6,7,8;

p31d+orow; end;

do i = 3,4,7,8;

p32d+orow; end;

do i=4, 8;

p33d+orow; end;

pid=ep12*p12d+ep13*p13d+ep14*p14d+ep21*p21d+ep22*p22d+ep23*p23d+ep24*p24d+ep31*p31d+ep32*p32d+ep33*p33d;

if round(pic,0.00000000001) gt 0 and round(pid,0.00000000001) gt 0 then do;

egamma=(pic-pid)/(pic+pid);

vgamma=16*((p11c*pid-p11d*pic)**2*ep11+(p12c*pid-p12d*pic)**2*ep12+

(p13c*pid-p13d*pic)**2*ep13+(p14c*pid-p14d*pic)**2*ep14+

(p21c*pid-p21d*pic)**2*ep21+

(p22c*pid-p22d*pic)**2*ep22+(p23c*pid-p23d*pic)**2*ep23+

(p24c*pid-p24d*pic)**2*ep24+

(p31c*pid-p31d*pic)**2*ep31+(p32c*pid-p32d*pic)**2*ep32+

(p33c*pid-p33d*pic)**2*ep33+(p34c*pid-p34d*pic)**2*ep34)/(n*(pic+pid)**4);

ratio=pic/pid;

vlratio=4*((p11c-ratio*p11d)**2*ep11+(p12c-ratio*p12d)**2*ep12+

(p13c-ratio*p13d)**2*ep13+(p14c-ratio*p14d)**2*ep14+

(p21c-ratio*p21d)**2*ep21+

(p22c-ratio*p22d)**2*ep22+(p23c-ratio*p23d)**2*ep23+

(p24c-ratio*p24d)**2*ep24+

(p31c-ratio*p31d)**2*ep31+(p32c-ratio*p32d)**2*ep32+

(p33c-ratio*p33d)**2*ep33+(p34c-ratio*p34d)**2*ep34)/(n*pic**2);

if vgamma gt 0 then do;

wt=1/vgamma;

sumwtwls+wt;

sumwlspt+wt*egamma;

end; ** end of if vgamma gt 0;

if vlratio gt 0 then do;

sumwtr+1/vlratio;

sumwtrpt+log(ratio)/vlratio;

end; ** end of if vration gt 0;

summhnum+n*(pic-pid);

summhden+n*(pic+pid);

a2+za**2*4*n**2*((p11c+p11d)**2*ep11+(p12c+p12d)**2*ep12+

(p13c+p13d)**2*ep13+(p14c+p14d)**2*ep14+

(p21c+p21d)**2*ep21+

(p22c+p22d)**2*ep22+(p23c+p23d)**2*ep23+

(p24c+p24d)**2*ep24+

(p31c+p31d)**2*ep31+(p32c+p32d)**2*ep32+

(p33c+p33d)**2*ep33+(p34c+p34d)**2*ep34)/n;

b2+za**2*4*n**2*((p11c-p11d)*(p11c+p11d)*ep11+

(p12c-p12d)*(p12c+p12d)*ep12+

(p13c-p13d)*(p13c+p13d)*ep13+

(p14c-p14d)*(p14c+p14d)*ep14+

(p21c-p21d)*(p21c+p21d)*ep21+

(p22c-p22d)*(p22c+p22d)*ep22+

(p23c-p23d)*(p23c+p23d)*ep23+

(p24c-p24d)*(p24c+p24d)*ep24+

(p31c-p31d)*(p31c+p31d)*ep31+

(p32c-p32d)*(p32c+p32d)*ep32+

(p33c-p33d)*(p33c+p33d)*ep33+

(p34c-p34d)*(p34c+p34d)*ep34)/n;

c2+za**2*4*n**2*((p11c-p11d)**2*ep11+(p12c-p12d)**2*ep12+

(p13c-p13d)**2*ep13+(p14c-p14d)**2*ep14+

(p21c-p21d)**2*ep21+

(p22c-p22d)**2*ep22+(p23c-p23d)**2*ep23+

(p24c-p24d)**2*ep24+

(p31c-p31d)**2*ep31+(p32c-p32d)**2*ep32+

(p33c-p33d)**2*ep33+(p34c-p34d)**2*ep34)/n;

end; ** end of if pic gt 0 and pid gt 0;

end; ** end of do numstrat;

do i= 1 to 12;

torow=tforow/tn;

end;

if round(sumwtwls,0.000000000001) eq 0 then do; flag1=1; end;

if round(summhnum,0.000000000001) eq 0 or round(summhden,0.000000000001) eq 0 then flag3=1;

if round(summhnum,0.000000000001) ne 0 and round(summhden,0.000000000001) ne 0 then mhpt=summhnum/summhden;

if sumwtwls gt 0 then do;

wlspt=sumwlspt/sumwtwls;

low1=wlspt-za/sqrt(sumwtwls);

if low1 le -1 then low1=-1;

up1=wlspt+za/sqrt(sumwtwls);

if up1 ge 1 then up1=1;

if low1 le gamma le up1 then sumcov1+1;

sumlen1+(up1-low1);

end; ** end of if sumwtwls gt 0;

agamma=summhden**2-a2;

bgamma=summhnum*summhden-b2;

cgamma=summhnum**2-c2;

if agamma le 0 or (bgamma**2-agamma*cgamma) le 0 then flag4=1;

if agamma gt 0 and (bgamma**2-agamma*cgamma) gt 0 then do;

low4=(bgamma-sqrt(bgamma**2-agamma*cgamma))/agamma;

if low4 le -1 then low4=-1;

up4=(bgamma+sqrt(bgamma**2-agamma*cgamma))/agamma;

if up4 ge 1 then up4=1;

if low4 le gamma le up4 then sumcov4+1;

sumlen4+(up4-low4);

end; ** end of if agamma gt 0 and (bgamma**2-agamma*cgamma) gt 0;

if round(sumwtr,0.000000000001) le 0 then flag5=1;

if round(sumwtr,0.000000000001) gt 0 then do;

rwlsl=exp(sumwtrpt/sumwtr-za/sqrt(sumwtr));

rwlsu=exp(sumwtrpt/sumwtr+za/sqrt(sumwtr));

low5=(rwlsl-1)/(rwlsl+1);

up5=(rwlsu-1)/(rwlsu+1);

if low5 le gamma le up5 then sumcov5+1;

sumlen5+(up5-low5);

end; ** end of if sumwtr gt 0;

tp11c=0; tp12c=0; tp13c=0; tp14c=0; tp21c=0; tp22c=0; tp23c=0; tp24c=0; tp31c=0;

tp32c=0; tp33c=0; tp34c=0;

do i = 6,7,8,10,11,12;

tp11c+torow; end;

do i = 7,8,11,12;

tp12c+torow; end;

do i=8,12;

tp13c+torow; end;

do i=10,11,12;

tp21c+torow; end;

do i =1,11,12;

tp22c+torow; end;

do i=1,2,12;

tp23c+torow; end;

do i=1,2,3;

tp24c+torow; end;

do i = 1, 5;

tp32c+torow; end;

do i=1, 2, 5,6;

tp33c+torow; end;

do i=1,2,3,5,6,7;

tp34c+torow; end;

i=1; tep11=torow;

i=2; tep12=torow;

i=3; tep13=torow;

i=4; tep14=torow;

i=5; tep21=torow;

i=6; tep22=torow;

i=7; tep23=torow;

i=8; tep24=torow;

i=9; tep31=torow;

i=10; tep32=torow;

i=11; tep33=torow;

i=12; tep34=torow;

tpic=tep11*tp11c+tep12*tp12c+tep13*tp13c+tep21*tp21c+tep22*tp22c+tep23*tp23c+tep24*tp24c+tep32*tp32c+tep33*tp33c+tep34*tp34c;

tp11d=0; tp12d=0; tp13d=0; tp14d=0; tp21d=0; tp22d=0; tp23d=0; tp24d=0; tp31d=0; tp32d=0; tp33d=0; tp34d=0;

do i = 5, 9;

tp12d+torow; end;

do i = 5, 6, 9, 10;

tp13d+torow; end;

do i = 5,6,7,9,10,11;

tp14d+torow; end;

do i = 2, 3, 4;

tp21d+torow; end;

do i = 3, 4, 9;

tp22d+torow; end;

do i = 4, 9, 10;

tp23d+torow; end;

do i = 9, 10, 11;

tp24d+torow; end;

do i = 2,3,4,6,7,8;

tp31d+torow; end;

do i = 3,4,7,8;

tp32d+torow; end;

do i=4, 8;

tp33d+torow; end;

tpid=tep12*tp12d+tep13*tp13d+tep14*tp14d+tep21*tp21d+tep22*tp22d+tep23*tp23d+tep24*tp24d+tep31*tp31d+tep32*tp32d+tep33*tp33d;

if round(tpic,0.00000000001) le 0 or round(tpid,0.000000000001) le 0 then do; flag6=1; flag7=1; end;

if round(tpic,0.00000000001) gt 0 and round(tpid,0.00000000001) gt 0 then do;

tegamma=(tpic-tpid)/(tpic+tpid);

tvgamma=16*((tp11c*tpid-tp11d*tpic)**2*tep11+(tp12c*tpid-tp12d*tpic)**2*tep12+

(tp13c*tpid-tp13d*tpic)**2*tep13+(tp14c*tpid-tp14d*tpic)**2*tep14+

(tp21c*tpid-tp21d*tpic)**2*tep21+

(tp22c*tpid-tp22d*tpic)**2*tep22+(tp23c*tpid-tp23d*tpic)**2*tep23+

(tp24c*tpid-tp24d*tpic)**2*tep24+

(tp31c*tpid-tp31d*tpic)**2*tep31+(tp32c*tpid-tp32d*tpic)**2*tep32+

(tp33c*tpid-tp33d*tpic)**2*tep33+(tp34c*tpid-tp34d*tpic)**2*tep34)/(tn*(tpic+tpid)**4);

low6=tegamma-za*sqrt(tvgamma);

if low6 le -1 then low6=-1;

up6=tegamma+za*sqrt(tvgamma);

if up6 ge 1 then up6=1;

if low6 le gamma le up6 then sumcov6+1;

sumlen6+(up6-low6);

tratio=tpic/tpid;

tvlratio=4*((tp11c-tratio*tp11d)**2*tep11+(tp12c-tratio*tp12d)**2*tep12+

(tp13c-tratio*tp13d)**2*tep13+(tp14c-tratio*tp14d)**2*tep14+

(tp21c-tratio*tp21d)**2*tep21+

(tp22c-tratio*tp22d)**2*tep22+(tp23c-tratio*tp23d)**2*tep23+

(tp24c-tratio*tp24d)**2*tep24+

(tp31c-tratio*tp31d)**2*tep31+(tp32c-tratio*tp32d)**2*tep32+

(tp33c-tratio*tp33d)**2*tep33+(tp34c-tratio*tp34d)**2*tep34)/(tn*tpic**2);

rlow7=tratio*exp(-1*za*sqrt(tvlratio));

rup7=tratio*exp(za*sqrt(tvlratio));

low7=(rlow7-1)/(rlow7+1);

up7=(rup7-1)/(rup7+1);

if low7 le gamma le up7 then sumcov7+1;

sumlen7+(up7-low7);

end; ** if tpic gt 0 and tpid gt 0;

vmh=0;

vlmhr=0;

if flag3 eq 0 then do;

do j = 1 to numstra;

p11c=0; p12c=0; p13c=0; p14c=0; p21c=0; p22c=0; p23c=0; p24c=0; p31c=0; p32c=0; p33c=0; p34c=0;

do i = 6,7,8,10,11,12;

p11c+orow; end;

do i = 7,8,11,12;

p12c+orow; end;

do i=8,12;

p13c+orow; end;

do i=10,11,12;

p21c+orow; end;

do i =1,11,12;

p22c+orow; end;

do i=1,2,12;

p23c+orow; end;

do i=1,2,3;

p24c+orow; end;

do i = 1, 5;

p32c+orow; end;

do i=1, 2, 5,6;

p33c+orow; end;

do i=1,2,3,5,6,7;

p34c+orow; end;

i=1; ep11=orow;

i=2; ep12=orow;

i=3; ep13=orow;

i=4; ep14=orow;

i=5; ep21=orow;

i=6; ep22=orow;

i=7; ep23=orow;

i=8; ep24=orow;

i=9; ep31=orow;

i=10; ep32=orow;

i=11; ep33=orow;

i=12; ep34=orow;

pic=ep11*p11c+ep12*p12c+ep13*p13c+ep21*p21c+ep22*p22c+ep23*p23c+ep24*p24c+ep32*p32c+ep33*p33c+ep34*p34c;

p11d=0; p12d=0; p13d=0; p14d=0; p21d=0; p22d=0; p23d=0; p24d=0; p31d=0; p32d=0; p33d=0; p34d=0;

do i = 5, 9;

p12d+orow; end;

do i = 5, 6, 9, 10;

p13d+orow; end;

do i = 5,6,7,9,10,11;

p14d+orow; end;

do i = 2, 3, 4;

p21d+orow; end;

do i = 3, 4, 9;

p22d+orow; end;

do i = 4, 9, 10;

p23d+orow; end;

do i = 9, 10, 11;

p24d+orow; end;

do i = 2,3,4,6,7,8;

p31d+orow; end;

do i = 3,4,7,8;

p32d+orow; end;

do i=4, 8;

p33d+orow; end;

pid=ep12*p12d+ep13*p13d+ep14*p14d+ep21*p21d+ep22*p22d+ep23*p23d+ep24*p24d+ep31*p31d+ep32*p32d+ep33*p33d;

vmh+4*n**2*((p11c-p11d-mhpt*(p11c+p11d))**2*ep11+(p12c-p12d-mhpt*(p12c+p12d))**2*ep12+

(p13c-p13d-mhpt*(p13c+p13d))**2*ep13+(p14c-p14d-mhpt*(p14c+p14d))**2*ep14+

(p21c-p21d-mhpt*(p21c+p21d))**2*ep21+

(p22c-p22d-mhpt*(p22c+p22d))**2*ep22+(p23c-p23d-mhpt*(p23c+p23d))**2*ep23+

(p24c-p24d-mhpt*(p24c+p24d))**2*ep24+

(p31c-p31d-mhpt*(p31c+p31d))**2*ep31+(p32c-p32d-mhpt*(p32c+p32d))**2*ep32+

(p33c-p33d-mhpt*(p33c+p33d))**2*ep33+(p34c-p34d-mhpt*(p34c+p34d))**2*ep34)/(n*summhden**2);

end; ** end of do numstrat;

tmhpt=0.5*log((1+mhpt)/(1-mhpt));

low3=tanh(tmhpt-za*sqrt(vmh/(1-mhpt**2)**2));

up3=tanh(tmhpt+za*sqrt(vmh/(1-mhpt**2)**2));

if low3 le gamma le up3 then sumcov3+1;

sumlen3+(up3-low3);

end; ** end of if flag3 eq 0;

end; ** end of do nsimul;

sflag1+flag1;

sflag2+flag2;

sflag3+flag3;

sflag4+flag4;

sflag5+flag5;

sflag6+flag6;

sflag7+flag7;

cov1=sumcov1/(nsimul-sflag1);

cov3=sumcov3/(nsimul-sflag3);

cov4=sumcov4/(nsimul-sflag4);

cov5=sumcov5/(nsimul-sflag5);

cov6=sumcov6/(nsimul-sflag6);

cov7=sumcov7/(nsimul-sflag7);

len1=sumlen1/(nsimul-sflag1);

len3=sumlen3/(nsimul-sflag3);

len4=sumlen4/(nsimul-sflag4);

len5=sumlen5/(nsimul-sflag5);

len6=sumlen6/(nsimul-sflag6);

len7=sumlen7/(nsimul-sflag7);

fail1=sflag1/nsimul;

fail3=sflag3/nsimul;

fail4=sflag4/nsimul;

fail5=sflag5/nsimul;

fail6=sflag6/nsimul;

fail7=sflag7/nsimul;

/* put gamma 1-5 2 numstra 6-10 strasize 11-15 cov1 16-23 3 cov3 24-31 3 cov4 32-39 3 cov5 40-47 3 cov6 48-55 3

cov7 56-63 3/

@18 "(" len1 19-23 3 @24 ")" @26 "(" len3 27-31 3 @32 ")"

@34 "(" len4 35-39 3 @40 ")" @42 "(" len5 43-47 3 @48 ")"

@50 "(" len6 51-55 3 @56 ")" @58 "(" len7 59-63 3 @64 ")"

/

@18 "[" fail1 19-23 3 @24 "]" @26 "[" fail3 27-31 3 @32 "]"

@34 "[" fail4 35-39 3 @40 "]" @42 "[" fail5 43-47 3 @48 "]"

@50 "[" fail6 51-55 3 @56 "]" @58 "[" fail7 59-63 3 @64 "]"

; */

put gamma 1-5 2 numstra 6-10 strasize 11-15 cov1 16-23 3 cov3 24-31 3 cov4 32-39 3 cov5 40-47 3 cov6 48-55 3

cov7 56-63 3/

@18 "(" len1 19-23 3 @24 ")" @26 "(" len3 27-31 3 @32 ")"

@34 "(" len4 35-39 3 @40 ")" @42 "(" len5 43-47 3 @48 ")"

@50 "(" len6 51-55 3 @56 ")" @58 "(" len7 59-63 3 @64 ")"

;

end; ** do strasize;

end; ** do numstrat;

end; ** do gamma;