# 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)