結果の診断プログラム

#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>
#include <oxdraw.h>
//Calculate variance of time series 
fTsvar(const x, const Bm){
decl n,sp;
// Calculate Spectral density function using
// Bandwidth=Bm at 2 points (0, pi) and Parzen  Window
n=rows(x);sp=periodogram(x,Bm,2,1)/n;
return M_2PI*sp[0];
}
// Calculate p-value for Convergence (Geweke's Method)
// H0: Convergence, H1: No Convergence
fGeweke(const x, const Bm){
decl n,n1,n2,sp1,sp2,var1,var2,x1,x2,x1bar,x2bar,z;
 n=rows(x);n1=floor(0.1*n);n2=floor(0.5*n);
 x1=x[:n1-1];x2=x[n-n2:];
 x1bar=meanc(x1);x2bar=meanc(x2);
 var1=fTsvar(x1,Bm);var2=fTsvar(x2,Bm);
 z=(x1bar-x2bar)/sqrt(var1/n1+var2/n2);
return 2*tailn(fabs(z));
}
// *************************************************
// Main Part 
// *************************************************
main(){
decl file,cm,crepeat,dbm,mx;
decl x,x1,x2,z,pvalue,i,j,y1,y2;
// rx: Number of Samples, bm:Bandowdth for Parzen Window
// ci: Quantiles for Credible Intervals
// dim_b: dimension of coeff. b
// n_acf: # of points for ACF plot

crepeat=10000;dbm=0.1*crepeat;cm=5;
// read beta 
 file = fopen("beta.txt");fscan(file,"%#m",crepeat,cm,&mx);fclose(file);
// 
// Summary Statistics
println(
"%r",{"Const"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mx[][0])~(varc(mx[][0])^0.5)~quantilec(mx[][0],<0.025,0.975>)'~fGeweke(mx[][0],dbm),"\n");
//
println(
"%r",{"PRIV"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mx[][1])~(varc(mx[][1])^0.5)~quantilec(mx[][1],<0.025,0.975>)'~fGeweke(mx[][1],dbm),"\n");

println(
"%r",{"SCHOOL"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mx[][2])~(varc(mx[][2])^0.5)~quantilec(mx[][2],<0.025,0.975>)'~fGeweke(mx[][2],dbm),"\n");

println(
"%r",{"LOGINC"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mx[][3])~(varc(mx[][3])^0.5)~quantilec(mx[][3],<0.025,0.975>)'~fGeweke(mx[][3],dbm),"\n");

println(
"%r",{"PTCON"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke"},
meanc(mx[][4])~(varc(mx[][4])^0.5)~quantilec(mx[][4],<0.025,0.975>)'~fGeweke(mx[][4],dbm),"\n");

decl n_acf=50;
//ACF
 DrawCorrelogram(0,mx[][0]',{"Const"},n_acf);
 DrawCorrelogram(1,mx[][1]',{"PRIV"},n_acf);
 DrawCorrelogram(2,mx[][2]',{"SCHOOL"},n_acf);
 DrawCorrelogram(3,mx[][3]',{"LOGINC"},n_acf);
 DrawCorrelogram(4,mx[][4]',{"PTCON"},n_acf);
 SaveDrawWindow("acf.ps");
 CloseDrawWindow();
//Path
 DrawTMatrix(0,mx[][0]',{"Const"},1,1,1);
 DrawTMatrix(1,mx[][1]',{"PRIV"},1,1,1);
 DrawTMatrix(2,mx[][2]',{"SCHOOL"},1,1,1);
 DrawTMatrix(3,mx[][3]',{"LOGINC"},1,1,1);
 DrawTMatrix(4,mx[][4]',{"PTCON"},1,1,1);
 SaveDrawWindow("path.ps");
 CloseDrawWindow();
//Posterior density
 DrawDensity(0,mx[][0]',{"Const"},1,0,0);
 DrawDensity(1,mx[][1]',{"PRIV"},1,0,0);
 DrawDensity(2,mx[][2]',{"SCHOOL"},1,0,0);
 DrawDensity(3,mx[][3]',{"LOGINC"},1,0,0);
 DrawDensity(4,mx[][4]',{"PTCON"},1,0,0);
 SaveDrawWindow("density.ps");
 CloseDrawWindow();

}