# MC for stochastic volatility model
# 03/06/2005

creastovol=function(n,beta2,sigma2,phi)
{

# Simulation of the data
# t=1,....,n,  y_t=beta*exp(alpha_t/2)*u_t
# alpha_1 ~ N(0,sigma2/(1-phi^2))
# t=2....,n, alpha_t=phi*alphz_{t-1}+sigma*v_t
# u_t ~ N(0,1) and v_t ~ N(0,1)

alpha=rep(0,n)
alpha[1]=rnorm(1,0,sqrt(sigma2/(1-phi^2)))
for (i in 2:n) alpha[i]=rnorm(1,phi*alpha[i-1],sqrt(sigma2))

y=rnorm(n,0,sqrt(beta2*exp(alpha)))
list(y=y,alpha=alpha)
}

#------------------------------------------------------------------------------#

lposshep=function(y,alpha,phi,sigma2,beta2)
{

# Log-posterior distribution up to a multiplicative constant
# f(y_1,....,y_n,alpha_1,....,alpha_n,beta2,phi,sigma2)

n=length(y)
0.5*log(1-phi^2)-0.5*(n+1)*log(beta2)-0.5*(n+1)*log(sigma2)-
0.5/sigma2*((1-phi^2)*alpha[1]^2+sum((alpha[2:n]-phi*alpha[1:(n-1)])^2))-0.5*sum(alpha)-
0.5/beta2*sum(y^2*exp(-alpha))
}

#------------------------------------------------------------------------------#

updalphashep=function(y,alphap,phi,sigma2,beta2)
{

# Proposal distribution for updating the latent vector alpha

n=length(y)
alpha=alphap
lopro=0
for (i in 1:n)
{
if (i==1)
{
mu=phi*alphap[2]
sigp2=1/sigma2+0.5*y[1]^2/exp(-mu)/beta2
mean=mu/sigma2+0.5*y[1]^2*exp(-mu)*(1+mu)/beta2-0.5
mean=mean/sigp2
}
if (i!=1 & i!=n)
{
mu=phi*(alpha[i-1]+alphap[i+1])/(1+phi^2)
sigp2=(1+phi^2)/sigma2+0.5*y[i]^2*exp(-mu)/beta2
mean=mu*(1+phi^2)/sigma2+0.5*y[i]^2*exp(-mu)*(1+mu)/beta2-0.5
mean=mean/sigp2
}
if (i==n)
{
mu=phi*alpha[n-1]
sigp2=1/sigma2+0.5*y[n]^2/exp(-mu)/beta2
mean=mu/sigma2+0.5*y[n]^2*exp(-mu)*(1+mu)/beta2-0.5
mean=mean/sigp2
}
alpha[i]=rnorm(1,mean,sqrt(1/sigp2))
lopro=lopro+dnorm(alpha[i],mean,sqrt(1/sigp2),log=T)
}
list(alpha=alpha,lopro=lopro)
}

#------------------------------------------------------------------------------#

loproalphashep=function(alpha,y,alphap,phi,sigma2,beta2)
{

# Weight associated to the proposal distribution of the latent vector alpha

n=length(y)
sigp2=rep(0,n)
mu=sigp2
mean=sigp2
mu[1]=phi*alphap[2]
mu[2:(n-1)]=phi*(alpha[1:(n-2)]+alphap[3:n])/(1+phi^2)
mu[n]=phi*alpha[n-1]
sigp2[1]=1/sigma2+0.5*y[1]^2/exp(-mu[1])/beta2
sigp2[2:(n-1)]=(1+phi^2)/sigma2+0.5*y[2:(n-1)]^2*exp(-mu[2:(n-1)])/beta2
sigp2[n]=1/sigma2+0.5*y[n]^2/exp(-mu[n])/beta2
mean[1]=mu[1]/sigma2+0.5*y[1]^2*exp(-mu[1])*(1+mu[1])/beta2-0.5
mean[2:(n-1)]=mu[2:(n-1)]*(1+phi^2)/sigma2+0.5*y[2:(n-1)]^2*exp(-mu[2:(n-1)])*(1+mu[2:(n-1)])/beta2-0.5
mean[n]=mu[n]/sigma2+0.5*y[n]^2*exp(-mu[n])*(1+mu[n])/beta2-0.5
mean=mean/sigp2
sum(dnorm(alpha,mean,sqrt(1/sigp2),log=T))
}

#------------------------------------------------------------------------------#

updparshep=function(y,alpha,phip)
{

# Proposal (Gibbs type) distribution for updating the parameters: sigma2, phi and beta2

n=length(y)
suffi1=0.5*sum(y^2*exp(-alpha))
beta2=1/rgamma(1,0.5*(n-1),rate=suffi1)
lopro=-suffi1/beta2-(0.5*(n+1))*log(beta2)+0.5*(n-1)*log(suffi1)
suffi2=0.5*((1-phip^2)*alpha[1]^2+sum((alpha[2:n]-phip*alpha[1:(n-1)])^2))
sigma2=1/rgamma(1,0.5*(n-1),rate=suffi2)
lopro=lopro-suffi2/sigma2-(0.5*(n+1))*log(sigma2)+0.5*(n-1)*log(suffi2)
suffi3=sum(alpha[2:(n-1)]^2)
mean=sum(alpha[2:n]*alpha[1:(n-1)])/suffi3
var=sigma2/suffi3
u=runif(1)
phi=mean+sqrt(var)*qnorm((u*(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var)))+pnorm((-1-mean)/sqrt(var))))
lopro=lopro+dnorm(phi,mean,sqrt(var),log=T)-log(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var)))
list(beta2=beta2,sigma2=sigma2,phi=phi,lopro=lopro)
}

#------------------------------------------------------------------------------#

loproparshep=function(beta2,sigma2,phi,y,alpha,phip)
{

# Weight associated to the proposal distribution of the parameters: sigma2, phi and beta2

n=length(y)
suffi1=0.5*sum(y^2*exp(-alpha))
suffi2=0.5*((1-phip^2)*alpha[1]^2+sum((alpha[2:n]-phip*alpha[1:(n-1)])^2))
lopro=-suffi1/beta2-(0.5*(n+1))*log(beta2)+0.5*(n-1)*log(suffi1)
lopro=lopro-suffi2/sigma2-(0.5*(n+1))*log(sigma2)+0.5*(n-1)*log(suffi2)
suffi3=sum(alpha[2:(n-1)]^2)
mean=sum(alpha[2:n]*alpha[1:(n-1)])/suffi3
var=sigma2/suffi3
lopro=lopro+dnorm(phi,mean,sqrt(var),log=T)-log(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var)))
lopro
}

#------------------------------------------------------------------------------#

hmstovol=function(niter,y)
{

# Hasting-Metropolis algorithm
# niter : number of iterations

n=length(y)
alpha=rep(0,n)
phi=rep(0,niter)
sigma2=rep(0,niter)
beta2=rep(0,niter)

phi[1]=2*runif(1)-1
sigma2[1]=abs(rnorm(1,0.5,1))
beta2[1]=abs(rnorm(1,0.5,1))

for (t in 1:(niter-1))
{
for (i in 1:n)
{

#Gibbs parameters for alpha[i]
if (i==1)
{
mu=phi[t]*alpha[2]
sigp2=1/sigma2[t]+0.5*y[1]^2*exp(-mu)/beta2[t]
mean=mu/sigma2[t]+0.5*y[1]^2*exp(-mu)*(1+mu)/beta2[t]-0.5
mean=mean/sigp2
}

if (i!=1 & i!=n)
{
mu=phi[t]*(alpha[i-1]+alpha[i+1])/(1+phi[t]^2)
sigp2=(1+phi[t]^2)/sigma2[t]+0.5*y[i]^2*exp(-mu)/beta2[t]
mean=mu*(1+phi[t]^2)/sigma2[t]+0.5*y[i]^2*exp(-mu)*(1+mu)/beta2[t]-0.5
mean=mean/sigp2
}

if (i==n)
{
mu=phi[t]*alpha[n-1]
sigp2=1/sigma2[t]+0.5*y[n]^2*exp(-mu)/beta2[t]
mean=mu/sigma2[t]+0.5*y[n]^2*exp(-mu)*(1+mu)/beta2[t]-0.5
mean=mean/sigp2
}

prop=rnorm(1,mean,sqrt(1/sigp2))

if (i==1)
{
d1=exp(-0.5*(prop^2-2*phi[t]*prop*alpha[2])/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[1]^2/beta2[t])
d2=exp(-0.5*(alpha[1]^2-2*phi[t]*alpha[2]*alpha[1])/sigma2[t]-0.5*alpha[1]-0.5*exp(-alpha[1])*y[1]^2/beta2[t])
}

if (i!=1 & i!=n)
{
d1=exp(-0.5*(prop^2*(1+phi[t]^2)-2*phi[t]*prop*(alpha[i-1]+alpha[i+1]))/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[i]^2/beta2[t])
d2=exp(-0.5*(alpha[i]^2*(1+phi[t]^2)-2*phi[t]*alpha[i]*(alpha[i-1]+alpha[i+1]))/sigma2[t]-0.5*alpha[i]-0.5*exp(-alpha[i])*y[i]^2/beta2[t])
}

if (i==n)
{
d1=exp(-0.5*(prop^2-2*phi[t]*prop*alpha[n-1])/sigma2[t]-0.5*prop-0.5*exp(-prop)*y[n]^2/beta2[t])
d2=exp(-0.5*(alpha[n]^2-2*phi[t]*alpha[n-1]*alpha[n])/sigma2[t]-0.5*alpha[n]-0.5*exp(-alpha[n])*y[n]^2/beta2[t])
}

prob=(d1*dnorm(alpha[i],mean,sqrt(1/sigp2)))/(d2*dnorm(prop,mean,sqrt(1/sigp2)))
if (runif(1)<prob) alpha[i]=prop
}

#Gibbs moves on sigma & beta
suffi1=0.5*sum(y^2*exp(-alpha))
beta2[t+1]=1/rgamma(1,0.5*(n-1),rate=suffi1)
suffi2=0.5*((1-phi[t]^2)*alpha[1]^2+sum((alpha[2:n]-phi[t]*alpha[1:(n-1)])^2))
sigma2[t+1]=1/rgamma(1,0.5*n,rate=suffi2)
suffi3=sum(alpha[2:(n-1)]^2)
mean=sum(alpha[2:n]*alpha[1:(n-1)])/suffi3
var=sigma2[t+1]/suffi3

#Proposal for phi
u=runif(1)
prop=mean+sqrt(var)*qnorm((u*(pnorm((1-mean)/sqrt(var))-pnorm((-1-mean)/sqrt(var)))+pnorm((-1-mean)/sqrt(var))))
prob=sqrt(1-prop^2)/sqrt(1-phi[t]^2)
if (runif(1)<prob) phi[t+1]=prop
else phi[t+1]=phi[t]

print(c(t,beta2[t],sigma2[t],phi[t]))
}

list(sigma2=sigma2,phi=phi,beta2=beta2)
}