############################################################################
#
# GARCH(1,1) on STAN/RSTAN & bayesGARCH
#
############################################################################
#install.packages("fGarch")
#install.packages("rstan")
library("rstan")
options(mc.cores=4)
library(tidyquant)
library(ggplot2)
options(stringsAsFactors = FALSE)
library(fGarch)
library(bayesGARCH)

# Getting PBR data from August 10th 2000 up to last friday, February 26th 2021.
getSymbols("PBR",from='2000-08-10',to='2021-02-26',warnings=FALSE,auto.assign=TRUE)    

n        = nrow(PBR)
price    = as.numeric(PBR[,6])
ret      = log(price[2:n]/price[1:(n-1)])
ret      = ret-mean(ret)
ret.std0 = ret/sqrt(var(ret))

par(mfrow=c(1,2))
ts.plot(price,xlab="August 10th 2000 - February 26th 2021",ylab="Price")
ts.plot(ret,xlab="August 10th 2000 - February 26th 2021",ylab="Log returns")

par(mfrow=c(2,2))
acf(ret,lag=100,main="Log returns")
pacf(ret,lag=100,main="Log returns")
acf(ret^2,lag=100,main="Squared log returns")
pacf(ret^2,lag=100,main="Squared log returns")



#############################
# Running stan code
#############################
y      = ret
n      = length(y)
sigma1 = sqrt(var(y[1:100]))
model  = stan_model("GARCH11.stan")
fit    = sampling(model,list(n=n,y=y,sigma1=sigma1),iter=500,chains=4)
params = extract(fit)
qsig   = t(apply(params$sigma,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(2,4))
ts.plot(params$mu)
ts.plot(params$alpha0)
ts.plot(params$alpha1)
ts.plot(params$beta1)
hist(params$mu,prob=TRUE)
hist(params$alpha0,prob=TRUE)
hist(params$alpha1,prob=TRUE)
hist(params$beta1,prob=TRUE)

par(mfrow=c(1,1))
plot(abs(y),type="h",col=grey(0.8),xlab="August 10th 2000 - February 26th 2021",
     ylab="Log returns (absolute value)")
lines(qsig[,2])


#############################
# Runing bayesGARCH
#############################
M0    = 2000
M     = 2000
niter = M0+M
run   = bayesGARCH(ret,mu.alpha=c(0,0),Sigma.alpha=1000*diag(1,2),mu.beta=0,
                   Sigma.beta=1000,lambda=0.01,delta=2,
                   control=list(n.chain=1,l.chain=niter,refresh=100))
draws = run$chain1

range = (M0+1):niter 

par(mfrow=c(2,4)) 
ts.plot(draws[range,1],xlab="iterations",main=expression(alpha[0]),ylab="") 
ts.plot(draws[range,2],xlab="iterations",main=expression(alpha[1]),ylab="") 
ts.plot(draws[range,3],xlab="iterations",main=expression(beta),ylab="") 
ts.plot(draws[range,4],xlab="iterations",main=expression(nu),ylab="")
for (i in 1:4)
  hist(draws[range,i],prob=TRUE)

sig2s = matrix(0,M,n)
for (t in 2:n)
  sig2s[,t] = draws[range,1]+draws[range,2]*ret[t-1]^2+draws[range,3]*sig2s[,t-1]
q.sig = t(apply(sqrt(sig2s),2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(1,1))
ts.plot(q.sig[,2])
lines(qsig[,2],col=2)

####################################
# Classical GARCH(1,1) - garchFit
####################################
fit.garch  = garchFit(~garch(1,1),data=ret,trace=F,include.mean=FALSE)    

fit.garch

res.garch     = fit.garch@residuals
sigmat.garch  = fit.garch@sigma.t
ret.std.garch = ret/sigmat.garch

par(mfrow=c(1,1))
ts.plot(q.sig[,2])
lines(qsig[,2],col=2)
lines(sigmat.garch,col=3)

par(mfrow=c(1,1))
plot(abs(y[(n-500):n]),type="h",xlab="August 10th 2000 - February 26th 2021",
     ylab="Log returns (absolute value)")
lines(q.sig[(n-500):n,2],col=2,lwd=3)
lines(qsig[(n-500):n,2],col=4,lwd=3)
lines(sigmat.garch[(n-500):n],col=6,lwd=3)
legend("topleft",legend=c("bayesGARCH","stan","garchFit"),col=c(2,4,6),lty=1,bty="n",lwd=3)