# Second implementation # with slice sampling for the h_t's # Data producion # warning: h = log sv # ie y = exp(h/2) eps # Parameters N=2500 rho=0.5 #Natural mu mu=2 #True mu mu=mu*(1-rho) #True sigma^2 sigma=0.1 # Data y=rnorm(N) h=rnorm(N)*sqrt(sigma) for (t in 2:N) h[t]=mu+h[t-1]*rho+h[t] y=exp(h/2)*y #Gibbs sampler T=10^4 hm=matrix(0,ncol=N,nrow=T) sighm=rep(1,T) muhm=0*sighm rhohm=0*sighm #Simulation of mu simu=function(h,sig,rho){ rnorm(1)*sqrt(sig/(N-.9))+(sum(h[-1])-rho*sum(h[-N]))/(N-.9) } #Simulation of rho siro=function(h,sig,mu){ rnorm(1)*sqrt(sig/(sum(h[-N]^2)+.1))+sum((h[-1]-mu)*h[-N])/(sum(h[-N]^2)+.1) } #Simulation of sigma sisi=function(h,rho,mu){ .5*(2+sum((h[-1]-mu-rho*h[-N])^2)+.1*mu^2+.1*rho^2)/rgamma(1,.5*N+2) } #Slice simulation of h[t,i] simh=function(t,i,mu,sigma,rho){ #transforms sigmai=sigma/(1+rho^2) if (i==1) mui=(mu*(1-rho)+rho*hm[t,i+1]-sigma)/(1+rho^2) if (i==N) mui=(mu*(1-rho)+rho*hm[t,i-1]-sigma)/(1+rho^2) if ((i>1)&&(i10^3) #subsampling indix=seq(1,T,le=10^3) par(mfrow=c(2,2),mar=c(2,1,1,1)) hist(muhm,ncl=65,col="midnightblue",main=expression(mu),xlab="") par(new=T) plot(muhm[indix],type="l",axes=F,xlab="",ylab="") hist(rhohm,ncl=65,col="midnightblue",main=expression(rho),xlab="") par(new=T) plot(rhohm[indix],type="l",axes=F,xlab="",ylab="") hist(sighm,ncl=65,col="midnightblue",main=expression(sigma^2),xlab="") par(new=T) plot(sighm[indix],type="l",axes=F,xlab="",ylab="") plot(apply(hm,2,mean),type="l",ylim=range(hm),col="sienna3",main=expression(h[t]),xlab="") for (t in indix) lines(hm[t,],col="grey50") lines(apply(hm,2,mean),col="sienna3",lwd=2) lines(h,col="forestgreen",lwd=2)