大標本理論に基づく仮説検定の有意水準を有限標本で評価

プログラム

#include<oxstd.h>
main(){
decl cK,cn,crepeat,i;
decl dc,dphi,dsigma2,ds2;
decl vb,vbeta,ve,vepsilon,veta,vreject,vtvalue,vx,vy;
decl mX,mSb;
// cn: sample size, crepeat:number of simulations
cn=30;crepeat=50000;
// true values
vbeta=<1;0.5>;dsigma2=1;dc=2;dphi=0.6;
cK=rows(vbeta);vreject=zeros(2,1);
// generate x
// x_0 => N(dc/(1-dphi), 1/(1-dphi^2)).
vx=zeros(cn,1);vx[0]=dc/(1-dphi)+rann(1,1)/sqrt(1-dphi^2);
// x_t=dc+dphi*x_{t-1}+eta_t, eta_t=> N(0,1)
veta=rann(cn-1,1);
for(i=0;i<cn-1;++i){vx[i+1]=dc+dphi*vx[i]+veta[i];}
//
vtvalue=zeros(crepeat,1);
for(i=0;i<crepeat-1;++i){
// Generate y_t=beta_1+beta_2*x_t+epsilon_t
//
// If epsilon_t => U(-0.5,-0.5), 
    vepsilon=ranu(cn,1)-0.5; 
// If epsilon_t => N(0,sigma2),
//  vepsilon=sqrt(dsigma2)*rann(cn,1);
 vy=vbeta[0]+vbeta[1]*vx+vepsilon;
 mX=ones(cn,1)~vx;          // X matrix
 vb=invert(mX'*mX)*mX'*vy;  // beta-hat
 ve=vy-mX*vb;               // residual
 ds2=ve'*ve/(cn-cK);        // estimate of sigma^2
 mSb=ds2*invert(mX'*mX);    // estimated covariance matrix of beta-hat 
// Calculate t ratio for H0: beta2=0
 vtvalue[i]=(vb[1]-vbeta[1])/sqrt(mSb[1][1]);
 if(tailt(fabs(vtvalue[i]),cn-cK) < 0.025){vreject[0]+=1;}
 if(tailn(fabs(vtvalue[i])) < 0.025){vreject[1]+=1;}
}
 println("Estimated Size of tests:", "    (Sample Size=",cn,")");
 println(
 "%r",{"Size"},
 "%c",{"True", "t", "Normal"},0.05~vreject'/crepeat);
}

計算結果

Estimated Size of tests:    (Sample Size=30)

               True            t       Normal
Size       0.050000     0.049580     0.058440