Posterior Simulation Comapriosn (flat prior)

プログラム

#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];
}
fTsvar_batch(const x, const n){
// N=n*m
decl i,m,N,y,xbar;
N=rows(x);m=N/n;xbar=meanc(x);y=zeros(m,1);
for(i=0;i<m;++i){y[i]=meanc(x[i*n:(i+1)*n-1]);}
return n/(m-1)*sumsqrc(y-xbar);
}
// 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(){
decl i,j,nobs,phi,phi_new,phi_s,mu_phi,var_phi,var_e,y;
decl burn,Bm,file,n,n_acf;
nobs=1001;y=zeros(nobs,1);phi=0.9;var_e=1.0;
// Generate Dataset to analyze
y[0]=rann(1,1)*sqrt(var_e/(1-phi^2));
for(i=0;i<nobs-1;++i){
  y[i+1]=phi*y[i]+rann(1,1)*sqrt(var_e);
  }
// Dataset
// Prior of phi ~ Uniform (-1,1), f(phi)=0.5 on (-1,1)
n=100000;burn=0.1*n;phi_s=zeros(n,1);
for(i=-burn;i<n;++i){
        mu_phi=(y[0:nobs-2]'*y[1:nobs-1])/sumsqrc(y[1:nobs-2]);
        var_phi=var_e/sumsqrc(y[1:nobs-2]);
        // Generate phi ~ TN(mu_phi,var_phi) on (-1,1)
        do{
        phi_new=mu_phi+rann(1,1)*sqrt(var_phi);
        }while(fabs(phi_new)>=1);
        // MH
        if(ranu(1,1)<= (0.5*sqrt(1-phi_new^2))/(0.5*sqrt(1-phi^2))){phi=phi_new;}
   // Add following Step to check the simulation procedure
     y[0]=rann(1,1)*sqrt(var_e/(1-phi^2));
     for(j=0;j<nobs-1;++j){
     y[j+1]=phi*y[j]+rann(1,1)*sqrt(var_e);
     }
   //
        if(i>=0){phi_s[i]=phi;}
 }

Bm=0.1*n;
//Bm=0;
println(
"%r",{"PHI"},
"%c",{"Mean","Stdev","95%L","95%U","Geweke","MC SE"},
meanc(phi_s)~varc(phi_s)^0.5~quantilec(phi_s,<0.025,0.975>)'~fGeweke(phi_s,Bm)~sqrt(fTsvar(phi_s,Bm)/n),"\n");
//
 n_acf=10000;
 //ACF
 DrawCorrelogram(0,phi_s',{"PHI"},n_acf);
 DrawTitle(0,"Sample ACF");
 //Path
 DrawTMatrix(1,phi_s',{"PHI"},1,1,1);
 DrawTitle(1,"Sample Path");
 //Posterior density
 DrawDensity(2,phi_s',{"Phi"},0,1,0);
 DrawTitle(2,"Posterior Density");
 SaveDrawWindow("psc-flat-prior.ps");
 CloseDrawWindow();
//
 file = fopen("phi.txt","w");fprint(file,"phi","%15.10f",phi_s);fclose(file);
 }

結果

               Mean        Stdev         95%L         95%U       Geweke
PHI       0.0046834      0.58563     -0.94853      0.95812      0.27859
              MC SE
PHI        0.049105