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