# Fourth implementation # with importance sampling for the h_t's # and fixed parameters (!!!) # Data producion # warning: h = log sv # ie y = exp(h/2) eps # Parameters N=100 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) #Importance simulation of h simh=function(t,i,beta,sigma,rho){ #transforms if (i>1) sigmai=1/(sigma^(-1)+.5*beta^(-2)*y[i]^2*exp(-rho*hm[t,i-1])) if (i==1) sigmai=1/(sigma^(-1)+.5*beta^(-2)*y[1]^2) if (i==1) mui=-.5*(1-beta^(-2)*y[1]^2)/(sigma^(-1)+.5*beta^(-2)*y[1]^2) if (i>1){ mui=rho*hm[t,i-1]-.5*(1-beta^(-2)*y[i]^2*exp(-rho*hm[t,i-1]))/(sigma^(-1) +.5*beta^(-2)*y[i]^2*exp(-rho*hm[t,i-1])) } #Generate h[t,i] as normal hti=sqrt(sigmai)*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(sigmai)) 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(sigmai)) 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)