############################################################################################
# Stochastic volatility
# Individual move algorithm for log-volatilities
# Random walk Metropolis
############################################################################################
#
# HEDIBERT FREITAS LOPES
# Associate Professor of Econometrics and Statistics
# The University of Chicago Booth School of Business
# 5807 South Woodlawn Avenue
# Chicago, Illinois, 60637
# Email : hlopes@ChicagoGSB.edu
# URL: http://faculty.chicagobooth.edu/hedibert.lopes/research/
#
############################################################################################
svol.rw = function(y,mu,phi,tau2,h0,h,vh){
  eta2  = c(rep(tau2/(1+phi^2),n-1),tau2)
  eta   = sqrt(eta2)
  coef1 = (1-phi)/(1+phi^2)*mu
  coef2 = phi/(1+phi^2)
  for (t in 1:n){
    if (t==n){
      mut  = mu + phi*h[t-1]
    }else{
      if (t==1){
        mut = coef1 + coef2*(h[2]+h0)
      	}else{
        mut = coef1 + coef2*(h[t+1]+h[t-1])
      }
    }
    ht   = rnorm(1,h[t],vh)
    num  = dnorm(ht,mut,eta[t],log=TRUE)+dnorm(y[t],0,exp(ht/2),log=TRUE)
    den  = dnorm(h[t],mut,eta[t],log=TRUE)+dnorm(y[t],0,exp(h[t]/2),log=TRUE)
    lacc = min(0,num-den)
    if (log(runif(1))<lacc){
      h[t] = ht
    }
  }
  return(h)
}

#-------------------------------------------------------
# Stochastic volatility
# Individual move algorithm for log-volatilities
#-------------------------------------------------------
svolsingle = function(y,mu,phi,tau2,h0,h){
  eta2  = c(rep(tau2/(1+phi^2),n-1),tau2)
  eta   = sqrt(eta2)
  coef1 = (1-phi)/(1+phi^2)*mu
  coef2 = phi/(1+phi^2)
  for (t in 1:n){
    if (t==n){
      mut  = mu + phi*h[t-1]
    }else{
      if (t==1){
        mut = coef1 + coef2*(h[2]+h0)
      	}else{
        mut = coef1 + coef2*(h[t+1]+h[t-1])
      }
    }
    mut1 = mut + 0.5*eta2[t]*(y[t]^2*exp(-mut)-1)     
    ht   = rnorm(1,mut1,eta[t])
    num  = dnorm(ht,mut,eta[t],log=TRUE)+
           dnorm(y[t],0,exp(ht/2),log=TRUE)-
           dnorm(ht,mut1,eta[t],log=TRUE)
    den  = dnorm(h[t],mut,eta[t],log=TRUE)+
           dnorm(y[t],0,exp(h[t]/2),log=TRUE)-
           dnorm(h[t],mut1,eta[t],log=TRUE)
    lacc = min(0,num-den)
    if (log(runif(1))<lacc){
      h[t] = ht
    }
  }
  return(h)
}

#-------------------------------------------------------
# Univariate FFBS: 
# y(t)     ~ N(alpha(t)+F(t)*theta(t);V(t))        
# theta(t) ~ N(gamma+G*theta(t-1);W)               
#-------------------------------------------------------
ffbsu = function(y,F,alpha,V,gamma,G,W,m0,C0){
  n = length(y)
  if (length(F)==1){F = rep(F,n)}
  if (length(alpha)==1){alpha = rep(alpha,n)}
  if (length(V)==1){V = rep(V,n)}
  a = rep(0,n)
  R = rep(0,n)
  m = rep(0,n)
  C = rep(0,n)
  B = rep(0,n-1)
  H = rep(0,n-1)
  a1 = gamma + G*m0
  R1 = G^2*C0 + W
  # time t=1
  a[1] = a1
  R[1] = R1
  f    = alpha[1]+F[1]*a[1]
  Q    = R[1]*F[1]**2+V[1]
  A    = R[1]*F[1]/Q
  m[1] = a[1]+A*(y[1]-f)
  C[1] = R[1]-Q*A**2
  # forward filtering
  for (t in 2:n){
    a[t] = gamma + G*m[t-1]
    R[t] = C[t-1]*G**2 + W
    f    = alpha[t]+F[t]*a[t]
    Q    = R[t]*F[t]**2+V[t]
    A    = R[t]*F[t]/Q
    m[t] = a[t]+A*(y[t]-f)
    C[t] = R[t]-Q*A**2
    B[t-1] = C[t-1]*G/R[t]
    H[t-1] = sqrt(C[t-1]-R[t]*B[t-1]**2)
  }
  theta = rep(0,n)
  # backward sampling
  theta[n] = rnorm(1,m[n],sqrt(C[n]))
  for (t in (n-1):1)
    theta[t] = rnorm(1,m[t]+B[t]*(theta[t+1]-a[t+1]),H[t])
  return(theta)
}

#-------------------------------------------------------
# y    = X*beta + u    u ~ N(0,sig2*I_n)
# beta|sig2 ~ N(b,sig2*inv(A))
#      sig2 ~ IG(v/2,v*lam/2)
#-------------------------------------------------------
fixedpar = function(y,X,b,A,v,lam){
  n     = length(y)
  k     = ncol(X)
  par1  = (v+n)/2
  var   = solve(crossprod(X,X)+A)
  mean  = matrix(var%*%(crossprod(X,y)+crossprod(A,b)),k,1)
  par2  = v*lam + sum((y-crossprod(t(X),mean))^2)
  par2  = (par2 + crossprod(t(crossprod(mean-b,A)),mean-b))/2
  sig2  = 1/rgamma(1,par1,par2)
  var   = var*sig2
  mean  = mean + crossprod(t(chol(var)),rnorm(2))
  return(c(mean,sig2))
}

#----------------------------------------------------------------------------
# Sample Z from 1,2,...,k, with P(Z=i) proportional to q[i]N(mu[i],sig2[i])
#----------------------------------------------------------------------------
ncind = function(y,mu,sig,q){
  w = dnorm(y,mu,sig)*q
  return(sample(1:length(q),size=1,prob=w/sum(w)))
}

#------------------------------------------------
# Quantile statistics
#------------------------------------------------
quant05 = function(x){quantile(x,0.05)}
quant95 = function(x){quantile(x,0.95)}

#------------------------------------------------
# Sampling the log-volatilities in a 
# standard univariate stochastic volatility model
#------------------------------------------------
svu = function(logy2,lambda,gamma,G,W,m0,C0){
  mu   = c(-11.40039,-5.24321,-9.83726,1.50746,-0.65098,0.52478,-2.35859)
  sig2 = c(5.795960,2.613690,5.179500,0.167350,0.640090,0.340230,1.262610)
  q    = c(0.007300,0.105560,0.000020,0.043950,0.340010,0.245660,0.257500)
  sig  = sqrt(sig2)
  z    = sapply(logy2-lambda,ncind,mu,sig,q)
  ffbsu(logy2,1.0,mu[z],sig2[z],gamma,G,W,m0,C0)
}