############################################################################
#
# Joint modelling rvol and sp500
#
# S&P 500 1-Month Realized Volatility Index
#
############################################################################
#install.packages("fGarch")
#install.packages("rstan")
library("rstan")
options(mc.cores=4)
library(tidyquant)
library(ggplot2)
options(stringsAsFactors = FALSE)
library(fGarch)

data = read.csv("rvol-sp500.csv",header=TRUE)

sp500  = data[,3]
n      = nrow(data)
rvol   = data[2:n,2]

ret    = log(sp500[2:n]/sp500[1:(n-1)])
ret    = (ret - mean(ret))/sqrt(var(ret))
n = n-1


pdf(file="rvol-sp500.pdf",width=12,height=9)
par(mfrow=c(1,2))
ts.plot(ret,ylab="",main="S&P 500 log returns",xlab="Feb 1 2011 - Feb 26 2021")
ts.plot(rvol,ylab="",main="S&P 500 1-Month Realized Volatility Index",xlab="Feb 1 2011 - Feb 26 2021")

par(mfrow=c(1,1))
ts.plot(abs(ret),ylab="",main="S&P 500 log returns vs realized volatility")
lines(sqrt(rvol),col=2)


#########################################
# GARCH(1,1) + RVOL
#########################################
sigma1 = sqrt(var(ret[1:100]))

model  = stan_model("rvol-sp500.stan")

fit    = sampling(model,list(n=n,ret=ret,rvol=rvol,sigma1=sigma1),iter=1000,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,xlab="Iterations",ylab="",main=expression(mu))
ts.plot(params$alpha0,xlab="Iterations",ylab="",main=expression(alpha[0]))
ts.plot(params$alpha1,xlab="Iterations",ylab="",main=expression(alpha[1]))
ts.plot(params$beta1,xlab="Iterations",ylab="",main=expression(beta[1]))
ts.plot(params$gamma0,xlab="Iterations",ylab="",main=expression(gamma[0]))
ts.plot(params$gamma1,xlab="Iterations",ylab="",main=expression(gamma[1]))
ts.plot(params$tau,xlab="Iterations",ylab="",main=expression(tau))

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

#########################################
# GARCH(1,1)
#########################################
model1  = stan_model("GARCH11.stan")
fit1    = sampling(model1,list(n=n,y=ret,sigma1=sigma1),iter=1000,chains=4)
params1 = extract(fit1)
qsig1   = t(apply(params1$sigma,2,quantile,c(0.025,0.5,0.975)))


par(mfrow=c(2,4))
ts.plot(params$mu,xlab="Iterations",ylab="",main=expression(mu))
ts.plot(params$alpha0,xlab="Iterations",ylab="",main=expression(alpha[0]))
ts.plot(params$alpha1,xlab="Iterations",ylab="",main=expression(alpha[1]))
ts.plot(params$beta1,xlab="Iterations",ylab="",main=expression(beta[1]))
hist(params1$mu,prob=TRUE,xlab="",main="")
hist(params1$alpha0,prob=TRUE,xlab="",main="")
hist(params1$alpha1,prob=TRUE,xlab="",main="")
hist(params1$beta1,prob=TRUE,xlab="",main="")

par(mfrow=c(2,2))
plot(density(params$mu),xlim=range(params$mu,params1$mu),main=expression(mu),xlab="")
lines(density(params1$mu),col=2)
plot(density(params$alpha0),xlim=range(params$alpha0,params1$alpha0),main=expression(alpha[0]),xlab="")
lines(density(params1$alpha0),col=2)
plot(density(params$alpha1),xlim=range(params$alpha1,params1$alpha1),main=expression(alpha[1]),xlab="")
lines(density(params1$alpha1),col=2)
plot(density(params$beta1),xlim=range(params$beta1,params1$beta1),main=expression(beta[1]),xlab="")
lines(density(params1$beta1),col=2)
legend("topleft",legend=c("GARC(1,1)","GARCH(1,1)+Rvol"),col=2:1,lty=1,bty="n")

par(mfrow=c(1,2))
hist(params1$alpha1+params1$beta1,xlim=c(0.9,1.0),
     xlab="GARCH(1,1)",main=expression(alpha[1]+beta[1]))
hist(params$alpha1+params$beta1,xlim=c(0.9,1.0),
     xlab="GARCH(1,1)+Rvol",main=expression(alpha[1]+beta[1]))


par(mfrow=c(1,1))
plot(abs(ret),type="h",col=grey(0.8),xlab="August 10th 2000 - February 26th 2021",
     ylab="Log returns (absolute value)",ylim=c(0,5))
lines(qsig[,2],col=1)
lines(qsig1[,2],col=2)
legend("topleft",legend=c("GARC(1,1)","GARCH(1,1)+Rvol"),col=1:2,lty=1,bty="n",lwd=2)


ind = (n-499):n
par(mfrow=c(1,1))
plot(abs(ret[ind]),type="h",col=grey(0.8),xlab="August 10th 2000 - February 26th 2021",
     ylab="Log returns (absolute value)",lwd=3,ylim=c(0,8))
lines(qsig[ind,2],col=1,lwd=3)
lines(qsig1[ind,2],col=2,lwd=3)
legend("topleft",legend=c("GARC(1,1)","GARCH(1,1)+Rvol"),col=2:1,lty=1,bty="n",lwd=3)


dev.off()