############################################################################
#
# Modeling time-varying variances via GARCH(1,1) vs SV-AR(1) models with
# Students' t errors.
#
############################################################################
#install.packages("fGarch")
library(tidyquant)
library(ggplot2)
options(stringsAsFactors = FALSE)
library(fGarch)
library(stochvol)
library(bayesGARCH)

pdf(file="GARCH11-SV1.pdf",width=12,height=8)

# Getting PBR data from August 10th 2000 up to last friday, March 5th 2021.
# -------------------------------------------------------------------------
getSymbols("PBR",from='2000-08-10',to='2021-03-05',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   = ret/sqrt(var(ret))
n     = length(ret)

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

# Fitting standard GARCH(1,1) with Gaussian errors
# ------------------------------------------------
fit.garch     = garchFit(~garch(1,1),data=ret,cond.dist="std",trace=F,include.mean=FALSE)
res.garch     = fit.garch@residuals
sigmat.garch  = fit.garch@sigma.t
ret.std.garch = ret/sigmat.garch

        # Estimate  Std. Error  t value Pr(>|t|)    
# omega   0.014364    0.003164    4.539 5.65e-06 ***
# alpha1  0.073478    0.009216    7.973 1.55e-15 ***
# beta1   0.910218    0.011001   82.739  < 2e-16 ***
# shape   7.025885    0.621874   11.298  < 2e-16 ***

# Fitting Bayesian GARCH(1,1) with Student's t errors
# ---------------------------------------------------
M0    = 2000
M     = 2000
niter = M0+M
range = (M0+1):niter 
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,]

m   = apply(draws,2,mean)
se  = sqrt(apply(draws,2,var))
tab = cbind(m,se,m/se,2*(1-pnorm(m/se)))
colnames(tab) = c("Estimate","Std. Error","t value","Pr(>|t|)")
tab

         # Estimate Std. Error   t value     Pr(>|t|)
# alpha0 0.01951128 0.00389684  5.006949 5.529955e-07
# alpha1 0.08608605 0.01035219  8.315736 0.000000e+00
# beta   0.89217674 0.01264973 70.529322 0.000000e+00
# nu     7.19102807 0.62112666 11.577394 0.000000e+00

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

ret.std.bgarch = ret/q.sig[,2]


par(mfrow=c(1,2))
ts.plot(q.sig[100:n,2],xlab="August 10th 2000 - March 5th 2021",ylab="Standard deviation")
lines(sigmat.garch[100:n],col=2)
title("GARCH(1,1) with Student's t errors\nClassical and Bayesian")
ts.plot(q.sig[(n-500):n,2],xlab="Last 500 days",ylab="Standard deviation")
lines(sigmat.garch[(n-500):n],col=2)
legend("topleft",legend=c("GARCH(1,1)","Bayesian GARCH(1,1)"),col=2:1,lty=1,bty="n")


# Fitting SV-AR(1)
fit.svar1 = svtsample(ret,draws=M,burnin=M0,thin=1)
params    = fit.svar1$para[[1]][,1:4]

m   = apply(params,2,mean)
se  = sqrt(apply(params,2,var))
tab = cbind(m,se,m/se,2*(1-pnorm(m/se)))
colnames(tab) = c("Estimate","Std. Error","t value","Pr(>|t|)")
tab

        # Estimate  Std. Error    t value     Pr(>|t|)
# mu    -0.4241572 0.134251199  -3.159429 1.998419e+00
# phi    0.9843182 0.003957667 248.711728 0.000000e+00
# sigma  0.1384504 0.015680593   8.829409 0.000000e+00
# nu    12.8571814 2.789885690   4.608498 4.055891e-06

q.vol = t(apply(exp(fit.svar1$latent[[1]]/2),2,quantile,c(0.025,0.5,0.975)))
ret.std.bsvar1 = ret/q.vol[,2]

par(mfrow=c(1,2))
ts.plot(q.sig[10:n,2],xlab="August 10th 2000 - March 5th 2021",ylab="Standard deviation")
lines(sigmat.garch[10:n],col=2)
lines(q.vol[10:n,2],col=4)
title("GARCH(1,1) with Student's t errors\nClassical and Bayesian")
ts.plot(q.sig[(n-1000):n,2],xlab="Last 1000 days",ylab="Standard deviation")
lines(sigmat.garch[(n-1000):n],col=2)
lines(q.vol[(n-1000):n,2],col=4)
legend("topleft",legend=c("GARCH(1,1)","Bayesian GARCH(1,1)","Bayesian SV-AR(1)"),col=c(2,1,4),lty=1,bty="n")


par(mfrow=c(3,3),mai=c(0.6,0.6,0.2,0.2))
ts.plot(ret.std.garch,ylab="Standardized returns",main="Student's t GARCH(1,1)")
ts.plot(ret.std.bgarch,ylab="Standardized returns",main="Bayesian Student's t GARCH(1,1)")
ts.plot(ret.std.bsvar1,ylab="Standardized returns",main="Bayesian Student's t SV-AR(1)")
breaks = seq(-10,10,length=30)
hist(ret.std.garch,prob=TRUE,breaks=breaks,xlab="",main="")
curve(dnorm,-10,10,add=TRUE,col=2,n=1000)
hist(ret.std.bgarch,prob=TRUE,breaks=breaks,xlab="",main="")
curve(dnorm,-10,10,add=TRUE,col=2,n=1000)
hist(ret.std.bsvar1,prob=TRUE,breaks=breaks,xlab="",main="")
curve(dnorm,-10,10,add=TRUE,col=2,n=1000)
qqnorm(ret.std.garch,xlim=c(-6,6))
qqline(ret.std.garch,col=2)
qqnorm(ret.std.bgarch)
qqline(ret.std.bgarch,col=2)
qqnorm(ret.std.bsvar1)
qqline(ret.std.bsvar1,col=2)

dev.off()