#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)&&(i<N))
    mui=(mu*(1-rho)+rho*(hm[t,i-1]+hm[t,i+1]))/(1+rho^2)
  
  phi=.5+((1-2*exp(sigmai))/(1-exp(sigmai)))
  lambda=(phi-1)*exp(mui+.5*sigmai) + .5*y[i]^2

  #proposition
  z=lambda/rgamma(1,shape=phi)

  #selection
  prob=targ(z,y[i],mui,sigmai)*ppal(exp(hm[t-1,i]),phi,lambda)
  prob=prob/(targ(exp(hm[t-1,i]),y[i],mui,sigmai)*ppal(z,phi,lambda))

  #back to log sv
  if (runif(1)<prob)
    newhm=log(z)
  else
    newhm=hm[t-1,i]

  newhm
  }

#First estimates of the h_t's
hm[1,]=h #log(y^2)+1.26

#First estimates of the parameters
#under the prior mu,rho N(0,10 sigma)
#and 1/sigma Exp(1)

sighm[1]=var(hm[1,])
rhohm[1]=acf(hm[1,],lag=1,plot=F)$lag[2]

muhm[1]=mu  #simu(hm[1,],sighm[1],rhohm[1])
rhohm[1]=rho #siro(hm[1,],sighm[1],muhm[1])
sighm[1]=sigma #sisi(hm[1,],rhohm[1],muhm[1])

#Full round
hm[2,]=hm[1,]

for (t in 2:T){

  for (i in sample(1:N))
    hm[t,i]=simh(t,i,muhm[t-1],sighm[t-1],rhohm[t-1])

  muhm[t]=simu(hm[t,],sighm[t-1],rhohm[t-1])
  rhohm[t]=siro(hm[t,],sighm[t-1],muhm[t])
  sighm[t]=sisi(hm[t,],rhohm[t],muhm[t])

  if (t<T) hm[t+1,]=hm[t,]
  } 

#Graphs
indix=1:T
if (T>10^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)