####################################################################################################
#
# Multivariate time-varying covariance modeling:
# Dynamic conditional correlation (Frequentist and Bayesian) and 
# Factor stochastic volatility (Bayesian)
#
# References:
#
# 1.Robert Engle (2002) Dynamic Conditional Correlation, JBES, 20(3), 339-350.
#
# 2.Fioruci, Ehlers and Andrade-Filho (2014) Bayesian multivariate GARCH models with DCC 
#   and asymmetric error distributions, Journal of Applied Statistics, 41(2), 320-331.
#   https://arxiv.org/pdf/1412.2967.pdf
#
# 3.Kastner, Fruehwirth-Schnatter and Lopes (2017) Efficient Bayesian inference for multivariate 
#   SV models, JCGS, 26, 905-917.
#   https://arxiv.org/pdf/1906.12123.pdf
#
# 4. I've also used this slides from Eric Zivot to understand the R package rmgarch
#    https://faculty.washington.edu/ezivot/econ589/DCCgarchPowerpoint.pdf
#
#
#  By Hedibert F. Lopes
#  March 31st 2021
#  www.hedibert.org
#
####################################################################################################


#install.packages("factorstochvol")
#install.packages("rmgarch")
#install.packages("~/Downloads/bayesDccGarch_2.1.tar", repos = NULL, type="source")
library(tidyquant)
library(rmgarch)
library("bayesDccGarch")
library(factorstochvol)


# Downloading the financial time series
# -------------------------------------
getSymbols("AAPL",from='2006-01-01',to='2021-03-25',warnings=FALSE,auto.assign=TRUE)# Apple Inc.
getSymbols("MSFT",from='2006-01-01',to='2021-03-25',warnings=FALSE,auto.assign=TRUE)# Microsoft Corp.
getSymbols("PBR",from='2006-01-01',to='2021-03-25',warnings=FALSE,auto.assign=TRUE)# Petrobras
getSymbols("VALE",from='2006-01-01',to='2021-03-25',warnings=FALSE,auto.assign=TRUE)# Vale

names = c("apple","microsoft","petrobras","vale")
data  = as.matrix(cbind(AAPL$AAPL.Adjusted,MSFT$MSFT.Adjusted,PBR$PBR.Adjusted,VALE$VALE.Adjusted))
n     = nrow(data)             
p     = ncol(data)
ret   = matrix(0,n-1,p)
for (i in 1:p)
 ret[,i] = log(data[2:n,i]/data[1:(n-1),i])
n = nrow(ret)

pdf(file="fsv-dcc-bayesdcc-US-BR.pdf",width=15,height=8)
             
par(mfrow=c(2,4))
for (i in 1:p)             
  ts.plot(data[,i],main=names[i],ylab="Price")
for (i in 1:p)             
  ts.plot(ret[,i],main=names[i],ylab="Log returns")


# Multivariate exploratory analysis
# ---------------------------------
L    = 25
ts   = (L+1):(n-L)
nt   = length(ts)
corr = array(0,c(nt,p,p))
cov  = array(0,c(nt,p,p))
for (i in 1:nt){
  ind = (ts[i]-L):(ts[i]+L)
  corr[i,,] = cor(ret[ind,])
  cov[i,,] = cov(ret[ind,])
}

par(mfrow=c(2,2))
for (i in 1:p)
  plot(ts,sqrt(cov[,i,i]),ylab="Standard deviation (rolling window)",main=names[i],type="l",xlab="Time")

rcorr = range(corr)
par(mfrow=c(2,3))
for (i in 1:(p-1))
for (j in (i+1):p){
  plot(ts,corr[,i,j],type="l",xlab="Time",ylab="Correlation (rolling window)",ylim=rcorr)
   title(paste(names[i]," x ",names[j],sep=""),cex=0.75)
}


# Multivariate GARCH modeling
# ---------------------------
garch11.spec = ugarchspec(mean.model=list(armaOrder=c(0,0)),variance.model=list(garchOrder=c(1,1),model = "sGARCH"),distribution.model="norm")

dcc.garch11.spec = dccspec(uspec=multispec(replicate(p,garch11.spec)),dccOrder=c(1,1),distribution="mvt")

dcc.fit = dccfit(dcc.garch11.spec,data=ret)


class(dcc.fit)

slotNames(dcc.fit)

names(dcc.fit@mfit)

names(dcc.fit@model)

dcc.cov = rcov(dcc.fit)
dcc.cor = rcor(dcc.fit)

par(mfrow=c(2,2))
for (i in 1:p){
  range = range(sqrt(dcc.cov[i,i,]),sqrt(cov[,i,i]))
  ts.plot(sqrt(dcc.cov[i,i,]),main=names[i],ylab="Standard deviation",ylim=range)
  lines(ts,sqrt(cov[,i,i]),col=2)
}
legend("topright",legend=c("DCC","Rolling window"),col=2:1,bty="n",lty=1)

  
par(mfrow=c(2,3))
for (i in 1:(p-1))
for (j in (i+1):p){
  plot(ts,corr[,i,j],ylim=range(dcc.cor,corr),ylab="Correlation",type="l",xlab="Time")
  lines(dcc.cor[i,j,],col=2)
  title(paste(names[i]," x ",names[j],sep=""),cex=0.75)
}
legend("bottomright",legend=c("DCC","Rolling window"),col=2:1,bty="n",lty=1)



# Multivariate SV modeling
# ------------------------
fsv.fit   = fsvsample(ret,factors=2,draws=5000,burnin=5000,runningstore=6)
fsv.cov   = array(0,c(n,p,p))
fsv.cor   = array(0,c(n,p,p))
l   = 0
for (j in 1:p){
  for (i in j:p){
    l = l + 1 
    fsv.cov[,i,j] = fsv.fit$runningstore$cov[,l,1]
    fsv.cov[,j,i] = fsv.cov[,i,j]
  }
}
for (i in 1:p)
  for (j in 1:p)
    fsv.cor[,i,j] = fsv.cov[,i,j]/sqrt(fsv.cov[,i,i]*fsv.cov[,j,j])


par(mfrow=c(2,2))
for (i in 1:p){
  range = range(sqrt(dcc.cov[i,i,]),sqrt(fsv.cov[,i,i]),sqrt(cov[,i,i]))
  plot(ts,sqrt(cov[,i,i]),main=names[i],ylab="Standard deviation",ylim=range,xlab="Time",type="l")
  lines(sqrt(dcc.cov[i,i,]),col=2)
  lines(sqrt(fsv.cov[,i,i]),col=3)  
}
legend("topright",legend=c("Rolling window","DCC","FSV"),col=1:3,bty="n",lty=1)
  

par(mfrow=c(2,3))
for (i in 1:(p-1))
for (j in (i+1):p){
    plot(ts,corr[,i,j],ylim=c(-0.5,1),ylab="",xlab="Time",type="l")
    lines(dcc.cor[i,j,],col=2)
    lines(fsv.cor[,i,j],col=3)
    title(paste(names[i]," x ",names[j],sep=""),cex=0.75)
  }
legend("bottomright",legend=c("Rolling window","DCC","FSV"),col=1:3,bty="n",lty=1)




# Bayesian DCC-GARCH modeling
# ---------------------------
bayesdcc = bayesDccGarch(ret,nSim=10000,errorDist=1)
bayesdcc = window(bayesdcc,start=5000)
summary(bayesdcc)

bdcc.cov = array(0,c(n,p,p))
bdcc.cor = array(0,c(n,p,p))
l = 0
for (i in 1:p){
  for (j in 1:p){
    l = l + 1
    bdcc.cov[,i,j] = bayesdcc$H[,l]
    bdcc.cor[,i,j] = bayesdcc$R[,l]
  }
}

par(mfrow=c(2,2))
for (i in 1:p){
  range = range(sqrt(dcc.cov[i,i,]),sqrt(cov[,i,i]),sqrt(bdcc.cov[,i,i]))
  ts.plot(sqrt(dcc.cov[i,i,]),main=names[i],ylab="Standard deviation",ylim=range)
  lines(ts,sqrt(cov[,i,i]),col=2)
  lines(sqrt(bdcc.cov[,i,i]),col=3)
}
legend("topright",legend=c("DCC","Bayes DCC","Rolling window"),col=c(2,3,1),bty="n",lty=1)

  
par(mfrow=c(2,3))
for (i in 1:(p-1))
for (j in (i+1):p){
  range = range(dcc.cor[i,j,],corr[,i,j],bdcc.cor[10:n,i,j])
  plot(ts,corr[,i,j],ylim=range,ylab="Correlation",type="l",xlab="Time")
  lines(dcc.cor[i,j,],col=2)
  lines(10:n,bdcc.cor[10:n,i,j],col=3)
  title(paste(names[i]," x ",names[j],sep=""),cex=0.75)
}
legend("bottomleft",legend=c("DCC","Bayes DCC","Rolling window"),col=c(2,3,1),bty="n",lty=1)



par(mfrow=c(2,1),mai=c(0.4,0.8,0.05,0.05))
ind = (n-1000):n
i=3
range = range(sqrt(dcc.cov[i,i,ind]),sqrt(bdcc.cov[ind,i,i]),sqrt(fsv.cov[ind,i,i]))
ts.plot(sqrt(dcc.cov[i,i,ind]),xlab="",main="",ylab="",ylim=range,lwd=2)
mtext(side=2,cex=1,line=2,"Standard deviation")
lines(sqrt(bdcc.cov[ind,i,i]),col=2,lwd=2)
lines(sqrt(fsv.cov[ind,i,i]),col=3,lwd=2)
legend("topleft",legend=c("DCC","Bayes DCC","FSV"),col=1:3,bty="n",lty=1,lwd=2)
legend("topright",legend=c(names[i],"last 1000 obs."))


i=3
j=4
ind = (n-1000):n
range = range(dcc.cor[i,j,ind],bdcc.cor[ind,i,j],fsv.cor[ind,i,j])
ts.plot(dcc.cor[i,j,ind],ylim=range,ylab="",xlab="")
mtext(side=2,cex=1,line=2,"Correlation")
lines(bdcc.cor[ind,i,j],col=2)
lines(fsv.cor[ind,i,j],col=3)
legend("topright",c(paste(names[i]," x ",names[j],sep=""),"last 1000 obs."),cex=0.75)
legend("topleft",legend=c("DCC","Bayes DCC","FSV"),col=1:3,bty="n",lty=1,lwd=2)

dev.off()