/*************size bcf paper this calculates size, power sbcfrh, heteroskedastic case ***********/ /***for power change line 34,61 only****/ /****The current version is Table 12, Row3, This corresponds to correlation of 0.5 between the instrument and the st, error, setup2, power***/ /****we have 1 instrument, 1 endog regressor, no controls****/ /****H0: theta0=0, so u=y*****/ cls; n=200; om={1 .3 0, .3 1 0.5, 0 .5 1}; /****cov. matrix is ready 1,2, 2,1 elements are cov zu******/ ci=1000;/*******iterations ********/ nbb=8;/****no. of power alternatives=8 when power, in size =3****/ rej=zeros(ci,nbb); rej1=zeros(ci,nbb); jkk=1; do while jkk<=nbb; /****for setups 1-3, for size, batch program****/ czuw=.5; czu=czuw; /***************/ /* for size only czu=czuw[jkk]; */ /****if setup 1, use czu, if setup 3 use czu2, setup 2: close both setups1,3, and then run see om matrix above*******/ om[1,2]=czu; om[2,1]=czu; /********************************/ /****cov. matrix is ready 1,2, 2,1 elements are cov zu******/ k=1; do while k<=ci; l=1;/********no. of instruments************/ dat=rndn(n,l+2); dat=dat*chol(om); u=dat[.,l+1];/********since bo=0, u=y, st.eqn*******/ tt={-2, -1.5, -1, -0.5, 0.5, 1, 1.5,2}; /********true val. for theta, 0 for siZe , and others for power sucgh as -2, -1.5*********/ v=dat[.,l+2]; z=dat[.,1:l];/********instruments**********/ p=ones(l,1)*2;/*******controls weakness of instruments******/ x=z*p+v;/**********red. form eqn*************/ u=u.*abs(z); y=x*tt[jkk]+u; /********st.eqn, put tt[jkk] for power*************/ /**********dgp is provided , now test stat for emp. likel.*******/ om1=z.*y; om11=(om1'om1/n); ar=(y'z)*inv(om11)*(z'y)/n; f=1/2-1.5/sqrt(n);/****automatic fraction choice for finite samples*****/ b=ceil(n*f); /***********resampling*********/ nb=1000;/***iterations to get resample distribution****/ arb=zeros(nb,1);/***********arb is the test stat at subsample**********/ i=1; do while i<=nb; kk4=zeros(b,1); ss=ceil(seqa(n,-1,b).*rndu(b,1))+seqa(0,1,b); xx=seqa(1,1,n); iii=0; do while iiicvalb); k=k+1; endo; "Rejection percentage of null, when theta =0, size, when theta not equal to 0, power:" meanc(rej[.,jkk]); jkk=jkk+1; endo;