############################################################################
#
# ARCH(1) and GARCH(1,1) model
#
# Hedibert Lopes
# March 1st 2021
#
# Econometrics III - Time Series
# PhD in Business Economics
# Insper Institute of Education and Research
#
#
############################################################################
#install.packages("fGarch")
library(tidyquant)
library(ggplot2)
options(stringsAsFactors = FALSE)
library(fGarch)

# 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")

# Fitting ARCH(1,1) & GARCH(1,1)
fit.arch   = garchFit(~garch(1,0),data=ret,trace=F,include.mean=FALSE)    
fit.garch  = garchFit(~garch(1,1),data=ret,trace=F,include.mean=FALSE)    

fit.arch
fit.garch

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

par(mfrow=c(2,3))
ts.plot(ret,xlab="08/10/2020 - 02/26/2021",ylab="Log returns",main="Fitting a ARCH(1) model")
lines(2*sigmat.arch,col=2)
lines(-2*sigmat.arch,col=2)
legend("bottomleft",legend=c("Log returns","+/- 2 standard deviations"),col=1:2,lwd=2,lty=1,bty="n")
ts.plot(ret.std.arch,xlab="08/10/2020 - 02/26/2021",ylab="Standardized log-returns",main="")
breaks = seq(min(ret.std.arch),max(ret.std.arch),length=50)
hist(ret.std.arch,xlab="",main="Standardized log returns\n ARCH(1) model",prob=TRUE,breaks=breaks)

ts.plot(ret,xlab="08/10/2020 - 02/26/2021",ylab="Log returns",main="Fitting a GARCH(1,1) model")
lines(2*sigmat.garch,col=2)
lines(-2*sigmat.garch,col=2)
legend("bottomleft",legend=c("Log returns","+/- 2 standard deviations"),col=1:2,lwd=2,lty=1,bty="n")
ts.plot(ret.std.garch,xlab="08/10/2020 - 02/26/2021",ylab="Standardized log-returns",main="")
breaks = seq(min(ret.std.garch),max(ret.std.garch),length=50)
hist(ret.std.garch,xlab="",main="Standardized log returns\n GARCH(1,1) model",prob=TRUE,breaks=breaks)

par(mfrow=c(2,3))
breaks = seq(min(ret.std.arch,ret.std.garch,ret.std0),max(ret.std.arch,ret.std.garch,ret.std0),length=50)
xxx = seq(breaks[1],breaks[50],length=100)
hist(ret.std0,xlab="",main="Standardized log returns\n Constant variance",prob=TRUE,breaks=breaks)
lines(xxx,dnorm(xxx),col=2)
hist(ret.std.arch,xlab="",main="Standardized log returns\n ARCH(1) model",prob=TRUE,breaks=breaks)
lines(xxx,dnorm(xxx),col=2)
hist(ret.std.garch,xlab="",main="Standardized log returns\n GARCH(1,1) model",prob=TRUE,breaks=breaks)
lines(xxx,dnorm(xxx),col=2)
qqnorm(ret.std0,ylim=range(breaks))
qqline(ret.std0, col = 2)
qqnorm(ret.std.arch,ylim=range(breaks))
qqline(ret.std.arch, col = 2)
qqnorm(ret.std.garch,ylim=range(breaks))
qqline(ret.std.garch, col = 2)