#Data producion #warning: h = log sv #ie y = exp(h/2) eps # Parameters N=1000 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=5*10^3 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) } #target function for simh #warning: h is now regular sv #ie exp(previous h) targ=function(h,y,mu,sigma){ h^(-1.5)*exp(-(y^2/(2*h))-(log(h)-mu)^2/(2*sigma) ) } #proposal for simh: inverted gamma #warning: h is now regular sv ppal=function(x,phi,lambda){ dgamma(lambda/x,shape=phi)/x^2 } #Simulation of h[t] #warning: h is now log sv 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])/(1+rho^2) if (i==N) mui=(mu*(1-rho)+rho*hm[t,i-1])/(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)