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