Ljung-BoxのQ統計量, Breusch-Godfrey　Test

#include<oxstd.h>
main(){
decl mData;
decl vyear,vmonth,vpai1,vpai3,vtb1,vtb3,vcpi;
decl vpai,vr,vdr;
decl cnlag,cnobs,cK;
decl file;
decl mX,vy,vyhat,vb,vpvalue,vse,vtvalue,ve;
decl dr2,ds2,drobustt;
decl mSb,mSxxinv,mS;
decl vrho,vg,vq,vqpv,i;
// Read Data
// Sample period 1950/2-1990/12
// year month pai1 pai3 tb1 tb3 cpi
file=fopen("mishkin.asc");fscan(file,"%#m",491,7,&mData);fclose(file);
vyear=mData[][0];vmonth=mData[][1];vpai1=mData[][2];
vpai3=mData[][3];vtb1=mData[][4];vtb3=mData[][5];vcpi=mData[][6];
//
format("%15.7g");
//format(200);println("%#13.7g",vyear~vmonth~vpai1~vpai3~vtb1~vtb3~vcpi);
println("Sample period: 1953/1-1971/7 (36th-258th. 233 obs.)");
vpai=((vcpi[35:257]./vcpi[34:256]).^(12)-1)*100;
cnobs=rows(vpai);vr=vtb1[35:257]-vpai;
// Ljung-Box Statistic for R
vdr=vr-meanc(vr);
cnlag=12;
vg=vrho=zeros(cnlag+1,1);
for(i=0;i<cnlag+1;++i){
vg[i]=(1/cnobs)*lag0(vdr,i)'*vdr;
}
vrho=(1/vg[0])*vg;
println("Autocorrelations for R and Std. Err.):",
        vrho[1:12]~(1/sqrt(cnobs))*ones(cnlag,1));
vq=zeros(cnlag,1);vq[0]=cnobs*(cnobs+2)*vrho[1]^2/(cnobs-1);
for(i=1;i<cnlag;++i){
vq[i]=vq[i-1]+cnobs*(cnobs+2)*vrho[i+1]^2/(cnobs-(i+1));
}
vqpv=tailchi(vq,cumsum(ones(cnlag,1),1));
println("Ljung-Box Statistics and p-values:",vq~vqpv);

// Regression: PAI = beta1+beta2*R
vy=vpai; 
// 
println("Mean of dependent variable(PAI)=",meanc(vy));
mX=ones(cnobs,1)~vtb1[35:257];
vb=invert(mX'*mX)*mX'*vy;
vyhat=mX*vb;
ve=vy-vyhat;cK=columns(mX);
ds2=ve'*ve/(cnobs-cK);
dr2=sumsqrc(vyhat-meanc(vyhat))/sumsqrc(vy-meanc(vy));
println("R-squared: ",dr2);
println("SER (Standard error of regression)=",sqrt(ds2));
mSxxinv=cnobs*invert(mX'*mX);
mS=zeros(cK,cK);
for(i=0;i<cnobs;++i){
mS+=ve[i]^2*mX[i][]'*mX[i][];
}
mS=(1/cnobs)*mS;
mSb=(1/cnobs)*mSxxinv*mS*mSxxinv;
vse=sqrt(diagonal(mSb))';
vtvalue=vb./vse;
vpvalue=2*tailt(fabs(vtvalue),cnobs-cK);
println(
 "%r",{"B1","B2"},
 "%c",{"Est. Coef.","Std. Err.","t-value", "p-value"},
       vb~vse~vtvalue~vpvalue);
// Robut t-ratio 
drobustt=(vb[1]-1.0)/vse[1];
println("Robut t-ratio for H0:beta2=1 vs H1:beta != 1: ",drobustt);
println("p-value: ",2*tailn(fabs(drobustt)), " (using normal) ",
2*tailt(fabs(drobustt),cnobs-cK),  " (using t)");
//
// Breusch-Godfrey Test
for(i=1;i<=cnlag;++i){mX=mX~lag0(ve,i);}
vy=ve;vb=invert(mX'*mX)*mX'*vy;vyhat=mX*vb;
dr2=sumsqrc(vyhat-meanc(vyhat))/sumsqrc(vy-meanc(vy));
println("Breusch-Godfrey Test of serial correaltion using p=12:");
println("nR^2: ",cnobs*dr2,"p-value: ",tailchi(cnobs*dr2,cnlag));
}