事後シミュレーション比較の結果の診断プログラム

#include <oxstd.h>
#include <oxprob.h>
#include <oxfloat.h>
#include <oxdraw.h>
//Calculate variance of time series 
fTsvar(const mX, const dBm){
// mX: T x m matrix (T: number of obs, m: number of series) 
decl cT,msp;
// Calculate Spectral density function using
// Bandwidth=Bm at 2 points (0, pi) and Parzen  Window
cT=rows(mX);msp=periodogram(mX,dBm,2,1)/cT;
// 1st row of msp is spectral density evaluated at 0
return M_2PI*(msp[0][])';
}
// Calculate p-value for Convergence (Geweke's Method)
// H0: Convergence, H1: No Convergence
fGeweke(const mX, const dBm){
// mX: T x m matrix (T: number of obs, m: number of series) 
decl cT,n1,n2,sp1,sp2,var1,var2,mX1,mX2,x1bar,x2bar,dz;
 cT=rows(mX);n1=floor(0.1*cT);n2=floor(0.5*cT);
 mX1=mX[:n1-1][];mX2=mX[cT-n2:][];
 x1bar=meanc(mX1)';x2bar=meanc(mX2)';
 var1=fTsvar(mX1,dBm);var2=fTsvar(mX2,dBm);
 dz=(x1bar-x2bar)./(var1/n1+var2/n2).^(0.5);
return 2*tailn(fabs(dz));
}
fMomentTest1(const mX, const vExpected, const dBm){
decl cT,vxbar,vvar,dz;
 cT=rows(mX);
 vxbar=meanc(mX)';vvar=fTsvar(mX,dBm);
 dz=(vxbar-vExpected)./(vvar/cT).^(0.5);
return 2*tailn(fabs(dz));

}
// *************************************************
// Main Part 
// *************************************************
main(){
decl file,cm,crepeat,dbm,mx,mz;
decl x,x1,x2,z,pvalue,i,j,y1,y2;
decl mstat;
// 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

// dbm=500 since autocorrelation vanishes after the lag 500.
crepeat=100000;dbm=500;cm=5;
// read beta 
 file = fopen("beta_psc.txt");fscan(file,"%#m",crepeat,cm,&mx);fclose(file);
// 
// Summary Statistics
mstat=meanc(mx)'~(varc(mx).^0.5)'~quantilec(mx,<0.025,0.975>)'
      ~fGeweke(mx,dbm);
println("%r",{"Const","PRIV","SCHOOL","LOGINC","PTCON"},
        "%c",{"Mean","Stdev","95%L","95%U","Geweke"},
        mstat);

// Posterior simulation comparison
//  Since b~N(0,I), we know that E(b_i)=0, E(b_i^2)=1, E(b_ib_j)=0
//  (so we don't have to sample from N(0,I).)
//  We will test them below. 

println("p-values of Posterior Simulation Comparison:");
decl dpvalue;
for(i=0;i<cm;i++){
 // Check the first moments
 // E(b_i)=0
 dpvalue=fMomentTest1(mx[][i],0,dbm);
 println("H_0:E(B",i,")=0",dpvalue);
 if(dpvalue <=0.01){println("  *Check errors in the MCMC program.*");} 
 // Check the second moments
  // E(b_ib_j)=0 if (i != 0)
  for(j=0;j<i;j++){
  dpvalue=fMomentTest1(mx[][i].*mx[][j],0,dbm);
  println("H_0:E(B",i,"*B",j,")=0",dpvalue);
  if(dpvalue <=0.01){println("  *Check errors in the MCMC program.*");} 
 }
  // E(b_i^2)=1
  dpvalue=fMomentTest1(mx[][i].^2,1,dbm);
  println("H_0:E(B",i,"^2)=1",dpvalue);
  if(dpvalue <=0.01){println("  *Check errors in the MCMC program.*");} 
 }
                
decl n_acf=1000;
//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_psc.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_psc.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_psc.ps");
 CloseDrawWindow();

}