@This example program is a GAUSS program to calculate
the empirical size and power of the t-test for H0: b=0,
where e follows t-distribution with 3 degrees of freedom.
The power is calculate for the case when b=0.5 etc. @

@ Modified from Ogaki, Jang and Lim (2004)@

RNDSEED 382974;

tend=25; @the sample size@
nor=1000; @the number of replications@
df=3; @ d.f. for the t-distribution of X@
i=1;
tn=zeros(nor,1); @used to store t-values under H0@
ta=zeros(nor,1); @used to store t-values under H1@
Z = RNDN(tend,1);
a = 2.5;
b = 0.1;
do until i>nor;
nrv=RNDN(tend,df+1); @normal r.v.'s@
crv=nrv[.,2:df+1]^2; @chi square r.v.'s@
e=nrv[.,1]./sqrt(sumc(crv')/df); @t distribution: used under H0@
@e=rndn(tend,1);@ 
X0 = a + e;
X1 = a + b*Z + e;
ZZ = ones(tend,1)~Z; 

b0 = invpd(ZZ'ZZ)*(ZZ'X0);
e0 = X0 - ZZ*b0;
sighat0 = (e0'e0)/(tend-2); @simgahat under H0@
b1 = invpd(ZZ'ZZ)*(ZZ'X1); 
e1 = X1 - ZZ*b1;
sighat1 = (e1'e1)/(tend-2); @sigmahat under H1@

vb0 = sighat0*invpd(ZZ'ZZ);
vb0 = vb0[2:2,2:2];

vb1 = sighat1*invpd(ZZ'ZZ);
vb1 = vb1[2:2,2:2];

tn[i]=b0[2]/sqrt(vb0); @t-value under H0@
ta[i]=b1[2]/sqrt(vb1); @t-value under H1@
i=i+1;
endo;
 "***** When H0 is true *****";
 "The estimated size with the nominal critical value";
 meanc(abs(tn).>1.96);
 "The estimated true 5-percent critical value";
sorttn=sortc(abs(tn),1);
etcv=sorttn[int(nor*0.95)];
 etcv;
 "***** When H1 is true *****";
 "The estimated power with the nominal critical value";
 meanc(abs(ta).>1.96);
 "The estimated size corrected power";
 meanc(abs(ta).>etcv);