###########################################
# Topic 6                                 # 
# Multivariate GARCH model                #
###########################################

rm(list = ls())
# setwd(.)

require(zoo)
require(xts)
require(knitr)
require(ggplot2)
require(rugarch)
require(rmgarch)
Sys.setlocale("LC_ALL","English")
source("Block2Functions.R")

# 1. Loading data for WIG20 stocks
####################################

filename <- "wig20.Rdata"
load(filename)

# filename  <- "http://web.sgh.waw.pl/~mrubas/EFII/Dane/wig.csv"
# data      <- read.csv(filename, row.names = "X")
# data      <- zoo(data,as.Date(rownames(data)))

dates     <- index(data) 
names(data)

# selecting the sample
startDate <- as.Date("2005-01-01")
endDate   <- as.Date("2050-01-01")
y         <- window(data, start=startDate, end=endDate)

# selecting a pair of stocks
y          <- na.omit(y[,c("pko","kgh")])

# log returns of stocks
dy  <- 100*diff(log(y))

# porfolio weights
w   <- c(0.5,0.5)

# 2. Some charts
#############################
# rolling moments

plot(dy)

wd <- 52
Rsd1  <- rollapply(dy[,1], width=wd, sd, align="right")
Rsd2  <- rollapply(dy[,2], width=wd, sd, align="right")
Rcor  <- rollapply(dy, width=wd, function(x) cor(x[,1],x[,2]), by.column=FALSE, align="right")
Rcov  <- rollapply(dy, width=wd, function(x) cov(x[,1],x[,2]), by.column=FALSE, align="right")

#windows()
par(mfrow=c(2,2), cex=0.7, bty="l")
plot(Rsd1, main="Rolling SD for asset 1", xlab="", ylab="")
plot(Rsd2, main="Rolling SD for asset 2", xlab="", ylab="")
plot(Rcor, main="Rolling correlation", xlab="", ylab=""); abline(h=0)
plot(Rcov, main="Rolling covariance", xlab="", ylab=""); abline(h=0)

############################################################
# 3. Model GO-GARCH                                        #
# Generalized Orthogonal GARCH model                       #
# x_t  = Z y_t, y_t ~ N(0,H_t)                             #
# H_t  = diag(h_1t,...,h_nt)                               #
# h_it = w_i + a_i y_{i,t-1}^2 + b_i h_{i,t-1}             #
# Var(x_t) = ZH_tZ'                                 	     #
# E(H_t) = I, so that unconditional variance of y_t is ZZ' #
# Spectral decomposition of uncond. variance:              #
# VL^(0.5)U U'L^(0.5)V' = ZZ', U - ortonormal matrix       #
############################################################

require(rmgarch)
# help(gogarchspec)
mod  = "gjrGARCH" # "sGARCH", "iGARCH"
pq   = c(1,1)

specGO <- gogarchspec(mean.model        = list(model = "constant", external.regressors = NULL),
                      variance.model     = list(model = mod, garchOrder = pq), 
                      distribution.model = c("mvnorm"))

# help(gogarchfit)
fitGO <- gogarchfit(specGO, dy)  # solver = "solnp"

# coefficients
fitGO@mfit$arcoef
fitGO@mfit$garchcoef

# some matrices and time series
resX <- fitGO@mfit$residuals      # orignal series (residuals)
resY <- fitGO@mfit$Y              # unobserved factors
U    <- fitGO@mfit$U              # orthonormal matrix U%*%t(U)=I
L    <- fitGO@mfit$D              # eigenvalues of ZZ': eigen(var(resX))
V    <- fitGO@mfit$E              # eigenvectors
Z    <- V%*%L^0.5%*%U             # check resX - resY%*%t(Z) = 0

# conditional sd and cor
# plot of conditional sd, corr, cov
Ht     <- rcov(fitGO)
dates  <- as.Date(names(Ht[1,1,]))
GOsd1  <- zoo(Ht[1,1,]^0.5,order.by = dates)
GOsd2  <- zoo(Ht[2,2,]^0.5,order.by = dates)
GOcov  <- zoo(Ht[2,1,],order.by = dates)
Rt     <- rcor(fitGO)
GOcor  <- zoo(Rt[2,1,],order.by = dates)

#windows()
par(mfrow=c(2,2), cex=0.75, bty="l", mar=c(3,3,3,3))
plot(GOsd1, main="conditional SD, asset 1", xlab="", ylab=""); 
plot(GOsd2, main="conditional SD, asset 2", xlab="", ylab=""); 
plot(GOcor, main="conditional correlation", xlab="", ylab=""); abline(h=0)
plot(GOcov, main="conditional covariance", xlab="", ylab=""); 

# GO-GARCH forecast
#####################

# help(gogarchforecast)
H     <- 10
fctGO <- gogarchforecast(fitGO, n.ahead = H) 

# the values of our forecasts
fitted(fctGO)  
rcov(fctGO)    
rcor(fctGO)    

HFt      <- rcov(fctGO)
datesF   <- as.Date(names(HFt)) + seq(1:H)         # for simplicity I don't adjust for working days calendar (it is doable in bizdays)
GOsd1F   <- zoo(HFt[[1]][1,1,]^0.5,order.by = datesF)
GOsd2F   <- zoo(HFt[[1]][2,2,]^0.5,order.by = datesF)
GOcovF   <- zoo(HFt[[1]][2,1,],order.by = datesF)
RFt      <- rcor(fctGO)
GOcorF   <- zoo(RFt[[1]][2,1,],order.by = datesF)

# plot
B      <- 100
x_lim  <- c(datesF[1]-B,datesF[H])

#windows()
par(mfrow=c(2,2), cex=0.75, bty="l", mar=c(3,3,3,3))
plot(GOsd1[dates>datesF[1]-B], main="conditional SD, asset 1", xlab="", ylab="", xlim=x_lim); lines(GOsd1F, col="red")
plot(GOsd2[dates>datesF[1]-B], main="conditional SD, asset 2", xlab="", ylab="", xlim=x_lim); lines(GOsd2F, col="red")
plot(GOcor[dates>datesF[1]-B], main="conditional correlation", xlab="", ylab="", xlim=x_lim); lines(GOcorF, col="red"); abline(h=0)
plot(GOcov[dates>datesF[1]-B], main="conditional covariance", xlab="", ylab="", xlim=x_lim); lines(GOcovF, col="red"); abline(h=0)

# Simulations and VaR/ES for the portolio
###########################################

# help(gogarchsim)
H      <- 10      # horizon
Nsim   <- 10000   # number of draws
p      <- 0.05    # tolerance level
simGO  <- gogarchsim(fitGO, n.sim = H, m.sim = Nsim, startMethod = c("sample"))

draws  <- simGO@msim$seriesSim # draws for returns

Rdraws <- matrix(NA,H,Nsim) # return from investing in our portfolio of 2 assets
for (m in 1:Nsim){
  Rdraws[,m]  = draws[[m]]%*%w    # w should be a vector of weights
}

temp     <- RdrawsToVaRES(RDraws,p)
VaRHgo   <- temp$VaR
ESHgo    <- temp$ES



###################################################################
# 4. DCC-GARCH                                                    #
#                                                                 #
# y_t = mu + e_t, e_t\sim N(0,\Sigma_t)                           #
# \Sigma_t = D_t*R_t*D_t,                                         # 
# e_t - D_t*z_t                                                   #
# R_t - conditional correlation matrix                            #
# D_t - diag(h_1t^0.5;...;h_Nt^0.5)                               #
# D_t - conditional variance                                      #
# h_t = a_0 + Ae_{t-1}.*e_{t-1} + Bh_{t-1}                        #
# R_t = (Q_t.* I)^{-0.5}Q_t(Q_t.* I)^{-0.5} scaling               #
# Q_t = (1-\alpha-\beta)Q + \alpha(z_{t-1}z_{t-1}' + \betaQ_{t-1} # 
###################################################################
require(rmgarch)

# help(dccspec)
# a. specification of univariate GARCH, separately for each individual asset

mod  = "gjrGARCH" # "sGARCH", "iGARCH"
pq   = c(1,1)
PQ   = c(0,0)
dist = "std"      # "norm"

specU = ugarchspec(variance.model=list(model=mod, garchOrder=pq), 
                   mean.model=list(armaOrder=PQ, include.mean=TRUE), 
                   distribution.model=dist)
mspec = multispec(c(specU, specU))


specDCC <- dccspec(mspec, dccOrder = c(1,1), model = "DCC")
fitDCC  <- dccfit(specDCC,dy)  # start.pars = list())
fitDCC

# some matrices and time series
Qbar  <- fitDCC@mfit$Qbar   
Qt    <- fitDCC@mfit$Q      # correlation matrix before scalling
Rt    <- fitDCC@mfit$R      # correlation matrix after scaling
Ht    <- fitDCC@mfit$H      # conditional covariance matrix

# compare average values of conditional covariance with unconditional covariance
apply(Ht,c(1,2),mean)
cov(dy)

# plot of conditional sd, corr, cov
Ht     <- rcov(fitDCC)
dates  <- as.Date(names(Ht[1,1,]))
DCCsd1  <- zoo(Ht[1,1,]^0.5,order.by = dates)
DCCsd2  <- zoo(Ht[2,2,]^0.5,order.by = dates)
DCCcov  <- zoo(Ht[2,1,],order.by = dates)
Rt     <- rcor(fitDCC)
DCCcor  <- zoo(Rt[2,1,],order.by = dates)


ylimC = c(min(GOcor,DCCcor)-0.01,max(GOcor,DCCcor)+0.01)
#windows()
par(mfrow=c(2,2), cex=0.75, bty="l",mar=c(3,3,3,3))
plot(DCCsd1, main="conditional SD, asset 1", xlab="", ylab=""); lines(GOsd1,col=2)
plot(DCCsd2, main="conditional SD, asset 2", xlab="", ylab=""); lines(GOsd2,col=2)
plot(DCCcor, main="conditional correlation", xlab="", ylab="",ylim=ylimC); abline(h=0); lines(GOcor,col=2)
plot(DCCcov, main="conditional covariance", xlab="", ylab=""); abline(h=0); lines(GOcov,col=2)

# Forecast
# help(dccforecast)
############################################

H      <- 10
fctDCC <- dccforecast(fitDCC, n.ahead = H) 

# the values of our forecasts
fitted(fctDCC)  # or fctDCC@mforecast$mu
rcov(fctDCC)    # or fctDCC@mforecast$H
rcor(fctDCC)    # or fctDCC@mforecast$R

HFt      <- rcov(fctDCC)
datesF   <- as.Date(names(HFt)) + seq(1:H)         # for simplicity I don't adjust for working days calendar
DCCsd1F  <- zoo(HFt[[1]][1,1,]^0.5,order.by = datesF)
DCCsd2F  <- zoo(HFt[[1]][2,2,]^0.5,order.by = datesF)
DCCcovF  <- zoo(HFt[[1]][2,1,],order.by = datesF)
RFt      <- rcor(fctDCC)
DCCcorF  <- zoo(RFt[[1]][2,1,],order.by = datesF)

# plot
B      <- 100
x_lim  <- c(datesF[1]-B,datesF[H])

#windows()
par(mfrow=c(2,2), cex=0.75, bty="l",mar=c(3,3,3,3))
plot(DCCsd1[dates>datesF[1]-B], main="conditional SD, asset 1", xlab="", ylab="", xlim=x_lim); lines(DCCsd1F, col="red")
plot(DCCsd2[dates>datesF[1]-B], main="conditional SD, asset 2", xlab="", ylab="", xlim=x_lim); lines(DCCsd2F, col="red")
plot(DCCcor[dates>datesF[1]-B], main="conditional correlation", xlab="", ylab="", xlim=x_lim); lines(DCCcorF, col="red"); abline(h=0)
plot(DCCcov[dates>datesF[1]-B], main="conditional covariance", xlab="", ylab="", xlim=x_lim); lines(DCCcovF, col="red"); abline(h=0)


# Simulations and VaR/ES for the portolio
###########################################

# help(dccsim)
H      <- 10      # horizon
p <- 0.05         # tolerance level for VaR/ES
Nsim   <- 10000   # number of draws
simDCC <- dccsim(fitDCC, n.sim = H, m.sim = Nsim, startMethod = c("sample"), rseed = 7)

draws  <- simDCC@msim$simX # draws for returns

Rdraws <- matrix(NA,H,Nsim) # return from investing in our portfolio of 2 assets
for (m in 1:Nsim){
  Rdraws[,m]  = draws[[m]]%*%w    # w should be a vector of weights
}

temp      <- RdrawsToVaRES(RDraws,p)
VaRHdcc   <- temp$VaR
ESHdcc    <- temp$ES

# 5. Saving the results
#########################

load("VaRESH.RData")

VaRH <- cbind(VaRH, VaRHdcc,VaRHgo)
ESH  <- cbind(ESH, ESHdcc,ESHgo)

save("VaRH","ESH", file = "VaRESH.RData")

require(knitr)
kable(cbind(1:H,-VaRH), digits=2)
kable(cbind(1:H,-ESH), digits=2)