# Fourth implementation # with importance sampling for the h_t's # and fixed parameters (!!!) # based on another approximation # Data producion # warning: h = log sv # ie y = exp(h/2) eps # Parameters N=200 rho=0.99 #Natural mu mu=-2 #True mu mu=mu*(1-rho) #True sigma^2 sigma=0.01 #True beta beta = exp(mu*(1-rho)/2) # Data y=rnorm(N) h=rnorm(N)*sqrt(sigma) for (t in 2:N) h[t]=h[t-1]*rho+h[t] y=beta*exp(h/2)*y #Importance sampler T=10^4 hm=matrix(0,ncol=N,nrow=T) # Values derived from the moments of log(N(0,1)) xhat=log(y^2)-log(beta^2)+1.27 tau=4.92 siginv=1/(sigma^(-1)+tau^(-1)) #Importance simulation of h simh=function(t,i,beta,sigma,rho){ #transforms if (i==1) mui=tau^(-1)*xhat[1]*siginv if (i>1) mui=(sigma^(-1)*rho*hm[t,i-1]+tau^(-1)*xhat[i])*siginv #Generate h[t,i] as normal hti=sqrt(siginv)*rnorm(1)+mui if (i>1) wei=exp(-.5*(hti-rho*hm[t,i-1])^2/sigma)*exp(-y[i]^2*beta^(-2)*exp(-hti)-.5*hti)/ dnorm(hti,mean=mui,sd=sqrt(siginv)) if (i==1) wei=exp(-.5*hti^2/sigma)*exp(-y[1]^2*beta^(-2)*exp(-hti)-.5*hti)/ dnorm(hti,mean=mui,sd=sqrt(siginv)) list(h=hti,w=wei) } iswe=rep(1,T) #Full round for (t in 1:T){ for (i in (1:N)){ resus=simh(t,i,beta,sigma,rho) hm[t,i]=resus\$h iswe[t]=iswe[t]*resus\$w } } #Graphs indix=sample(1:T,10^3,iswe,rep=T) plot(apply(hm,2,weighted.mean,w=iswe/sum(iswe)), 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)