4つのHeteroschedasticity-consistentな推定方法の比較 (HCTYPE=0,1,2,3)

#include<oxstd.h>
main(){
decl mData;
decl vyear,vmonth,vpai1,vpai3,vtb1,vtb3,vcpi;
decl vpai,vr,vdr;
decl cnobs,cK;
decl file,i;
decl mX,vy,vyhat,vb,vpvalue,vse,vtvalue,ve;
decl mSb0,mSb1,mSb2,mSb3,mSxxinv,mS0,mS1,mS2,mS3,mH;
// 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];
//
// Use 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;
// Regression: PAI = beta1+beta2*R
vy=vpai;mX=ones(cnobs,1)~vtb1[35:257];
vb=invert(mX'*mX)*mX'*vy;
vyhat=mX*vb;
ve=vy-vyhat;cK=columns(mX);
mSxxinv=cnobs*invert(mX'*mX);
mS0=mS1=mS2=mS3=zeros(cK,cK);
mH=mX*invert(mX'*mX)*mX';
for(i=0;i<cnobs;++i){
mS0+=(1/cnobs)*ve[i]^2*mX[i][]'*mX[i][];
mS1+=(1/(cnobs-cK))*ve[i]^2*mX[i][]'*mX[i][];
mS2+=(1/cnobs)*(ve[i]^2/(1-mH[i][i]))*mX[i][]'*mX[i][];
mS3+=(1/cnobs)*(ve[i]^2/(1-mH[i][i])^2)*mX[i][]'*mX[i][];
}
mSb0=(1/cnobs)*mSxxinv*mS0*mSxxinv;
mSb1=(1/cnobs)*mSxxinv*mS1*mSxxinv;
mSb2=(1/cnobs)*mSxxinv*mS2*mSxxinv;
mSb3=(1/cnobs)*mSxxinv*mS3*mSxxinv;
vse=sqrt(diagonal(mSb0))|sqrt(diagonal(mSb1))
    |sqrt(diagonal(mSb2))|sqrt(diagonal(mSb3));
vtvalue=vb'./vse;
vpvalue=2*tailt(fabs(vtvalue),cnobs-cK);
format("%10.6f");
println("Estimated Beta1 is: ",vb[0]);
println(
 "%r",{"HCTYPE=0","HCTYPE=1","HCTYPE=2","HCTYPE=3"},
 "%c",{"Std. Err.","t-value", "p-value"},
       vse[][0]~vtvalue[][0]~vpvalue[][0]);
println("Estimated Beta2 is: ",vb[1]);
println(
 "%r",{"HCTYPE=0","HCTYPE=1","HCTYPE=2","HCTYPE=3"},
 "%c",{"Std. Err.","t-value", "p-value"},
       vse[][1]~vtvalue[][1]~vpvalue[][1]);
}