############################################################################################
#
# 
# Generalized autoregressive conditional heteroscedasticity (ARCH) model 
#
#  Model: y(t)    = sig(t)*e(t)                                
#         sig2(t) = a1+a2*y(t-1)^2+a3*sig2(t-1)
#
# where e(t) ~ N(0,1)    t=1,...,n.
#  
############################################################################################
#
#  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/
#
############################################################################################
q025 = function(x){quantile(x,0.025)}
q975 = function(x){quantile(x,0.975)}
like = function(sig2){
  llike = sum(dnorm(y,0,sqrt(sig2),log=TRUE))
}
sigma2s = function(a){
  sig2 = rep(0,n)
  sig2[1] = a[1]+a[2]*y20+a[3]*sig20
  for (t in 2:n)
    sig2[t] = a[1]+a[2]*y[t-1]^2+a[3]*sig2[t-1]
  return(sig2)	
}


# Simulating data
set.seed(1255)
n         = 300
a.true    = c(0.1,0.4,0.59)
sig20     = 0.1
y20       = 0.1
sig2      = rep(0,n)
y         = rep(0,n)
sig2[1]   = a.true[1]+a.true[2]*y20+a.true[3]*sig20
y[1]      = rnorm(1,0,sqrt(sig2[1]))
for (t in 2:n){
  sig2[t]   = a.true[1]+a.true[2]*y[t-1]^2+a.true[3]*sig2[t-1]
  y[t] = rnorm(1,0,sqrt(sig2[t]))
}

sig2.true = sigma2s(a.true)
  
pdf(file="garch-data.pdf",width=10,height=5)
par(mfrow=c(1,1))
plot(y,xlab="Time",ylab="",main=expression(y[t]),type="b")
dev.off()

# Prior hyperparameters
a0 = c(0,0,0)
C0 = rep(10,3)

# MCMC setup
M0    = 10000
M     = 1000
L     = 100
niter = M0+M*L
v     = 0.01
a     = a.true
sig2  = sigma2s(a)

pars  = a
set.seed(1235)
garchtime = system.time({
for (i in 1:niter){
  a.new = rnorm(3,a,v)
  sig2.new = sigma2s(a.new)
  if (sum(sig2.new>0)==n){
    alpha = min(0,like(sig2.new)-like(sig2))
    if (log(runif(1))<alpha){
      a = a.new
      sig2 = sig2.new
    }
  }
  pars = rbind(pars,a)
}
})
pars = pars[seq(M0+1,niter,by=L),]

true = a.true
names = c("a1","a2","a3")
pdf(file="garch-mcmc.pdf",width=10,height=7)
par(mfrow=c(3,3))
for (i in 1:3){
  ts.plot(pars[,i],xlab="iteration",ylab="",main=names[i])
  abline(h=true[i],col=2,lwd=3)
  acf(pars[,i],main="")
  hist(pars[,i],xlab="",ylab="",main="")
  abline(v=true[i],col=2,lwd=3)
}
dev.off()

pdf(file="garch-pars.pdf",width=10,height=7)
par(mfrow=c(2,2))
plot(pars[,1:2],xlab=expression(a[1]),ylab=expression(a[2]))
plot(pars[,c(1,3)],xlab=expression(a[1]),ylab=expression(a[3]))
plot(pars[,2:3],xlab=expression(a[2]),ylab=expression(a[3]))
hist(pars[,2]+pars[,3],xlab="",ylab="",main=expression(a[2]+a[3]))
abline(v=sum(a.true[2:3]),lwd=3,col=2)
dev.off()

sig2s = NULL
for (i in 1:M)
  sig2s = rbind(sig2s,sigma2s(pars[i,]))
  
  
msig2 = apply(sig2s,2,median)
lsig2 = apply(sig2s,2,q025)
usig2 = apply(sig2s,2,q975)

pdf(file="garch-vars.pdf",width=10,height=5)
par(mfrow=c(1,2))
ts.plot(cbind(lsig2,msig2,usig2),col=c(2,1,2),main="")
plot(sig2.true,msig2,xlab="TRUE",ylab="POSTERIOR MEAN",pch=16,cex=0.75)
abline(0,1,col=2,lwd=3)
dev.off()

# Forecasting
k     = 10
sig2f = cbind(sig2s,matrix(0,M,k))
y1    = cbind(y[n],matrix(0,M,k))
for (t in 1:k){
  sig2f[,n+t]= pars[,1]+pars[,2]*y1[,t]^2+pars[,3]*sig2f[,n+t-1]
  y1[,t+1] = rnorm(M,0,sqrt(sig2f[,n+t]))
}
sig2f  = sig2f[,(n+1):(n+k)]
msig2f = apply(sig2f,2,median)
lsig2f = apply(sig2f,2,q025)
usig2f = apply(sig2f,2,q975)

pdf(file="garch-vars-forecast.pdf",width=10,height=7)
par(mfrow=c(1,1))
ts.plot(cbind(lsig2f,msig2f,usig2f),col=c(2,1,2),main="")
dev.off()