#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter5: Section5.9.3
#
# Sensitivity analysis of Value at Risk
#
# Author : Tomohiro Ando
#
################################################################################

library(bayesGARCH)


# Nikkei 225 index data

x <- scan("Nikkei225-index.txt")

n <- length(x)


#MCMC estimation of GARCH with student-t innovations

MCMC <- bayesGARCH(x,control=list(n.chain=1,l.chain=10000,thin=10,start=10001))


#Summarize results

A <- as.matrix(MCMC$chain1)

Var <- rep(0,len=10000)

for(i in 1:10000){

a0 <- as.numeric(A[i,1])
a1 <- as.numeric(A[i,2])
b  <- as.numeric(A[i,3])
nu <- as.numeric(A[i,4])

h <- a0+a1*x[1]^2
for(j in 2:(n+1)){h <- a0+a1*x[j-1]^2+b*h}

Var[i] <- sqrt(h*(nu-2)/nu)*qt(0.01,df=nu)

}

mean(Var)

plot(hist(Var),main="",xlab=expression(Var[99]),xlim=c(-3,-1.5))
lines(c(mean(Var),mean(Var)),c(-10,10000),lwd=3)






#MCMC estimation of GARCH with (approximate) normal innovations

MCMC <- bayesGARCH(x,lambda = 0.01, delta = 50, control=list(n.chain=1,l.chain=10000,thin=10,start=10001))


#Summarize results

A <- as.matrix(MCMC$chain1)

Var <- rep(0,len=10000)

for(i in 1:10000){

a0 <- as.numeric(A[i,1])
a1 <- as.numeric(A[i,2])
b  <- as.numeric(A[i,3])
nu <- as.numeric(A[i,4])

h <- a0+a1*x[1]^2
for(j in 2:(n+1)){h <- a0+a1*x[j-1]^2+b*h}

Var[i] <- sqrt(h*(nu-2)/nu)*qnorm(0.01,0,1)

}

mean(Var)

plot(hist(Var),main="",xlab=expression(Var[99]),xlim=c(-3,-1.5))
lines(c(mean(Var),mean(Var)),c(-10,10000),lwd=3)