install.packages("tidyquant")
install.packages("fGarch")
install.packages("bayesGARCH")
library("tidyquant")
library("fGarch")
library("bayesGARCH")


getSymbols('PBR',src='yahoo',return.class='ts',from="2019-01-01")
n = nrow(PBR)

price = PBR[,6]

par(mfrow=c(1,1))
ts.plot(price,ylab="Price")

ret = 100*(price[2:n]/price[1:(n-1)]-1)

n = length(ret)

hist(ret,breaks=seq(min(ret),max(ret),length=50))


ts.plot(ret,ylab="Return")


fit.garch = garchFit(~garch(1,1),data=ret,trace=F,include.mean=FALSE,cond.dist=c("std"))    

fit.garch

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

par(mfrow=c(2,2))
range  = seq(min(ret.std.garch),max(ret.std.garch),length=30)
range1 = seq(min(ret.std.garch),max(ret.std.garch),length=100)
ts.plot(ret,ylab="Return")
ts.plot(ret.std.garch,xlab="Days",ylab="Std returns",main="")
acf(ret.std.garch,main="")
hist(ret.std.garch,breaks=range,prob=TRUE,main="",xlab="")
lines(range1,dnorm(range1,mean(ret.std.garch),sqrt(var(ret.std.garch))),col=2,lwd=2)

par(mfrow=c(1,1))
ts.plot(abs(ret),xlab="Days",ylab="Standard deviations",main="GARCH(1,1) model")
lines(sigmat.garch,col=2)


omega   = 0.15501
alpha1  = 0.04386
beta1   = 0.93748
shape   = 4.67726
par.mle = c(omega,alpha1,beta1,shape)

sigmas = rep(sigmat.garch[1]^2,n)
for (t in 2:n)
  sigmas[t] = omega+alpha1*ret[t-1]^2+beta1*sigmas[t-1]
sigmas = sqrt(sigmas)

ts.plot(sigmas)
lines(sigmat.garch,col=2)
cbind(sigmas,sigmat.garch)[1:10,]


M0 = 10000
M   = 1000
thin = 50
niter = M0+thin*M
time = system.time({
  fit.bayes = bayesGARCH(ret,control = list(n.chain = 2, l.chain = niter))
})
time
#    user  system elapsed 
# 126.658  24.870 151.496
# M0=10K, M=1K, thin=50

names = c("omega","alpha1","beta1","shape")
ind = seq(M0+1,niter,by=thin)
par(mfrow=c(2,2))
for (i in 1:4){
  ts.plot(fit.bayes$chain1[ind,i],ylab="",main=names[i])
  lines(fit.bayes$chain2[ind,i],col=2)
}

omega.b  = c(fit.bayes$chain1[ind,1],fit.bayes$chain2[ind,1])
alpha1.b = c(fit.bayes$chain1[ind,2],fit.bayes$chain2[ind,2])
beta1.b  = c(fit.bayes$chain1[ind,3],fit.bayes$chain2[ind,3])
shape.b  = c(fit.bayes$chain1[ind,4],fit.bayes$chain2[ind,4])

pars = cbind(omega.b,alpha1.b,beta1.b,shape.b)

par(mfrow=c(2,2))
for (i in 1:4)
  ts.plot(pars[,i],xlab="",main=names[i])

par(mfrow=c(2,2))
for (i in 1:4)
  acf(pars[,i],main="")

par(mfrow=c(2,2))
for (i in 1:4){
  hist(pars[,i],prob=TRUE,xlab="",main=names[i])
  abline(v=par.mle[i],col=2,lwd=2)
}

omega.m  = mean(c(fit.bayes$chain1[ind,1],fit.bayes$chain2[ind,1]))
alpha1.m = mean(c(fit.bayes$chain1[ind,2],fit.bayes$chain2[ind,2]))
beta1.m  = mean(c(fit.bayes$chain1[ind,3],fit.bayes$chain2[ind,3]))
shape.m  = mean(c(fit.bayes$chain1[ind,4],fit.bayes$chain2[ind,4]))

c(omega,alpha1,beta1,shape)
c(omega.m,alpha1.m,beta1.m,shape.m)

sigmas.m = rep(sigmat.garch[1]^2,n)
for (t in 2:n)
  sigmas.m[t] = omega.m+alpha1.m*ret[t-1]^2+beta1.m*sigmas.m[t-1]
sigmas.m = sqrt(sigmas.m)

par(mfrow=c(1,1))
ind1 = (n-500):n
ts.plot(sigmat.garch[ind1],ylim=range(sigmat.garch[ind1],sigmas.m[ind1]),ylab="Standard deviation")
lines(sigmas.m[ind1],col=2)
legend("topleft",legend=c("GARCH","Bayes GARCH"),col=1:2,bty="n",lty=1,lwd=2)


sigmas.b = matrix(sigmat.garch[1]^2,2*M,n)
for (t in 2:n)
  sigmas.b[,t] = omega.b+alpha1.b*ret[t-1]^2+beta1.b*sigmas.b[,t-1]
sigmas.b = sqrt(sigmas.b)

qsigmas.b = t(apply(sigmas.b,2,quantile,c(0.025,0.5,0.975)))

h = 400
ts.plot(qsigmas.b[(n-h):n,],col=c(3,2,3))
lines(sigmas.m[(n-h):n],col=4)