###########################################
# Topic 5                                 # 
# Univariate GARCH model                  #
###########################################

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

require(zoo)
require(xts)
require(moments)
require(forecast)
require(tseries)
require(rugarch)
require(knitr)
require(ggplot2)
source("Block2Functions.R")

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

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

# filename  <- "http://web.sgh.waw.pl/~mrubas/EFII/Dane/wig20.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)

# portfolio returns
r   <- zoo(dy%*%w,index(dy))
R   <- as.numeric(coredata(r))

# the value of investment in the portfolio
P   <- exp(cumsum(r/100))

par(mfrow=c(2,1), cex = 0.75, bty="l")
plot(P, main="price")
plot(r, main="log returns"); abline(h=0)

# 2. Some statistics and charts
###################################

# summary statistics
Nyear <- 365/as.numeric(mean(diff(dates)))
mu    <- mean(r)*Nyear
sig   <- sd(r)*sqrt(Nyear) 
mom <- as.data.frame(c(Nyear,mu,sig,min(r),max(r), skewness(r), kurtosis(r),jarque.bera.test(R)$stat))
rownames(mom) <- c("N","mu","sig","min","max","skew","kurt", "JB test"); colnames(mom)="value"
kable(mom, digits=3)

dev.off()
Acf(R, main="ACF of daily returns" )
Box.test(R, lag = 1, type = c("Ljung-Box"))

# Density plot
Rstar <- (R-mean(R))/sd(R)
bwdth     <- 0.05
dens_plot <- ggplot(data.frame(Rstar), aes(x = Rstar)) +
  geom_histogram(binwidth = bwdth, colour = "white", fill = "blue", size = 0.1) +
  theme_light() +
  labs(title="Empirical vs theoretical distribution", y="", x="Standarized returns", caption="") +
  stat_function(fun = function(x) ddist("norm",x)*length(Rstar)*bwdth, color = "red", size = 2) 
dens_plot

# adding t-Student
v = 4 + 6/(kurtosis(R)-3)
dens_plot <- dens_plot + stat_function(fun = function(x) ddist("std",x,shape=v)*length(Rstar)*bwdth, color = "black", size = 2)    # rozklad t-Studenta
dens_plot
dens_plot+ xlim(-7,-2)

# QQ plot
q            <- seq(1/length(Rstar),1-1/length(Rstar), 1/length(Rstar))
# Qteo         <- qnorm(q)        # rozklad normalny
# Qteo  <- qt(q,v)*sqrt((v-2)/v)             # for t-Student the variance is v/(v-2)
Qteo  <- qdist("std",p=q,shape=v)

Qemp         <- quantile(Rstar,q)         # our data
# Qemp         <- qnorm(q)                 
# Qemp         <- quantile(rnorm(10000),q)  # with some randomness


QQtab = data.frame(Qteo,Qemp)
ggplot(QQtab) + 
  geom_point(aes(x=Qemp, y=Qteo), size=2,shape=23, col="blue")+
  theme_light()+
  labs(title="QQ plot") +
  xlim(min(QQtab),max(QQtab))+ylim(min(QQtab),max(QQtab)) +
  geom_abline(intercept=0, slope=1, size=1)


# 3a. VaR/ES for 1-steap ahead horizon: 
# Normal, t-Student, Historical Simulation
######################################################

# normal distribution 
p   <- 0.05      # tolerance level
m   <- mean(R) 
s   <- sd(R)
VaRnorm     <- qnorm(p)*s + m             
ESnorm      <- m - s*dnorm(qnorm(p))/p    
qf          <- function(x) qdist("norm", p=x)
ESnorm1     <- m + s*(1/p * integrate(qf, 0, p)$value) 

# t-Student distribution 
# v = 4 + 6/(kurtosis(R)-3)               # GMM estimate
# v <- fitdist("std", R)$pars["shape"]    # ML estimate
v = 5
VaRt  <- m + s*qdist("std",shape=v,p=p)
qf    <- function(x) qdist("std", p=x, shape=v)
ESt   <- m + s*(1/p * integrate(qf, 0, p)$value) 

# Historical simulation
R0    <- sort(R)              
N0    <- floor(length(R0)*p)                            
VaRhs <- R0[N0]                      # compare to: quantile(R,p)
EShs  <- mean(R0[1:N0])              

# 3b. VaR/ES for H-steap ahead horizon: 
# Normal, t-Student, Historical Simulation
######################################################

H = 10

# normal distribution: square root of time
VaRHnorm <- sqrt(1:H)*qnorm(p)*s + (1:H)*m    
ESHnorm  <- (1:H)*m - sqrt(1:H)*s*dnorm(qnorm(p))/p    

# t-distribution: Monte Carlo simulations
Nsim    <- 10000
Rdraws  <- matrix(rdist(distribution="std", Nsim*H, mu = m, sigma = s, shape = v),H,Nsim) 
RdrawsC <- apply(Rdraws,2,cumsum)               # cumualted changes
VaRHt   <- ESHt <- rep(NaN,H)
M0  <- floor(Nsim*p)     # observation for p-th quantile 
for(h in 1:H){
  temp   = sort(RdrawsC[h,])
  VaRHt[h] = temp[M0]
  ESHt[h]  = mean(temp[1:M0])
}

# The same function in Block2Functions.R
# temp  <- RdrawsToVaRES(RDraws,p)
# VaRHt <- temp$VaR
# ESHt  <- temp$ES

# Historical Simulation: Bootstrap
Rdraws  <- matrix(sample(R, Nsim*H, replace=TRUE),H,Nsim) 
temp   <- RdrawsToVaRES(RDraws,p)
VaRHhs <- temp$VaR
ESHhs  <- temp$ES

kable(data.frame(h=1:H,VaRHhs,VaRHnorm, VaRHt),digits=2)

# 4. Volatility clustering
###########################################

Acf(R^2, main="ACF of squared daily returns" )
Box.test(R^2, lag = 20, type = c("Ljung-Box"))

# min-max spread
rollmin <- rollapply(r, width=Nyear, min, by=1)
rollmax <- rollapply(r, width=Nyear, max, by=1)
plot(merge(rollmin,rollmax),screens=c(1,1))

# rolling standard deviation 
sigROLL <- rollapply(r, width=Nyear, sd, by=1, align="right")
plot(sigROLL)

# exponentially weighted moving average - EWMA/riskMetrics
lambda      <- 0.94                          # smoothing parameter
temp        <- rep(0,length(R))
temp[1]     <- var(R)                        # starting point
for (t in 2:length(R)){             
  temp[t] = lambda * temp[t-1]  + (1-lambda) * R[t-1]^2
}
sigEWMA      <- zoo(temp^0.5,order.by=index(r))
plot(sigEWMA)
lines(sigROLL,col='red')

# EWMA in rugarch package
require(rugarch)
# help(ugarchspec)

EWMAspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), include.mean=FALSE),
                       variance.model=list(model="iGARCH", garchOrder=c(1,1)),
                       fixed.pars=list(alpha1=1-lambda, omega=0), distribution.model = "norm")
EWMAfit  <- ugarchfit(data = r, spec = EWMAspec)
# calibrated params.
round(coef(EWMAfit),3)

# comparison
tail(merge(as.zoo(sigma(EWMAfit)),sigEWMA))
plot(merge(as.zoo(sigma(EWMAfit)),sigEWMA), plot.type="single", col=c(1,2))

plot(EWMAfit, which =1)


# 5. GARCH(1,1) model
#######################

spec0 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(1,1)), 
                   mean.model=list(armaOrder=c(0,0), include.mean=TRUE),  
                   distribution.model="norm")

# help(ugarchfit) 
fit0 = ugarchfit(data=r, spec=spec0) # solver="nlminb"

fit0
slotNames(fit0)
sigGARCH <- as.zoo(sigma(fit0))

# equilibrium value for variance
sig2eq = fit0@fit$coef["omega"]/(1-fit0@fit$coef["alpha1"]-fit0@fit$coef["beta1"])

par(mfrow=c(1,1), cex = 0.7, bty="l")
plot(sigGARCH)
lines(sigEWMA,col='red')
lines(zoo(sqrt(sig2eq),index(sigGARCH)))

# standarized residuals
u   <- as.zoo(residuals(fit0)/sigma(fit0))

par(mfrow=c(1,2), cex=0.7, bty="l", lwd=1)
plot(u, main="standarized residuals", xlab="", ylab="")
Acf(coredata(u)^2, main="ACF for squarred standarized residuals", xlab="", ylab="")

Box.test(u^2, lag=12, type='Ljung')


# 6. Selecting the specificatio of GARCH(p,q)  
#############################################
# dist = c("norm", "snorm", "std", "sstd")

LagSel <- function(x, Pmax=4, Qmax=4, crit="SIC", dist="norm"){
  IC <- matrix(NA, Pmax, Qmax+1)
  for(p in 1:Pmax){
    for(q in 0:Qmax){
      
      spec = ugarchspec(variance.model=list(model="sGARCH", garchOrder=c(p,q)), 
                        mean.model=list(armaOrder=c(0,0), include.mean=TRUE),  
                        distribution.model=dist)
      fit  = ugarchfit(data=x, spec=spec)
      if(crit == "AIC"){IC[p,q+1] <- infocriteria(fit)[1] }
      if(crit == "SIC"){IC[p,q+1] <- infocriteria(fit)[2] }
      if(crit == "HQ"){	IC[p,q+1] <- infocriteria(fit)[4] }
    }
  }
  rownames(IC) <- paste('p=',1:Pmax, sep="")
  colnames(IC) <- paste('q=',0:Qmax, sep="")
  return(IC)
}

ICtab <- LagSel(r,4,4,crit="SIC", dist="std")
kable(ICtab,digits=3)

# optimal lags
pq   = c(1,1)
PQ   = c(0,0)
dist = "std"

spec1 = ugarchspec(variance.model=list(model="sGARCH", garchOrder=pq), 
                   mean.model=list(armaOrder=PQ, include.mean=TRUE),  
                   distribution.model=dist)
fit1 = ugarchfit(data=r, spec=spec1)
plot(fit1, which=12)


# 7. Extensions GARCH(1,1)  
######################################

# Leverage effect
Ccf(R^2,R, type="correlation", main="cross correlation of squarred returns and returns")

# assymetric models
spec.e   = ugarchspec(variance.model=list(model="eGARCH", garchOrder=pq), 
                      mean.model=list(armaOrder=PQ, include.mean=TRUE),  
                      distribution.model=dist)
spec.gjr = ugarchspec(variance.model=list(model="gjrGARCH", garchOrder=pq), 
                      mean.model=list(armaOrder=PQ, include.mean=TRUE),  
                      distribution.model=dist)

fit.e   = ugarchfit(data=r, spec=spec.e)    # solver="nlminb"
fit.gjr = ugarchfit(data=r, spec=spec.gjr)

# info criteria
IC <- cbind(infocriteria(fit1), infocriteria(fit.e), infocriteria(fit.gjr))
colnames(IC) <- c("GARCH", "eGARCH", "gjrGARCH")
kable(IC,digits=3)

mod = "gjrGARCH"

# GARCH-in-Mean #
spec.m = ugarchspec(variance.model=list(model=mod, garchOrder=pq), 
                    mean.model=list(armaOrder=PQ, include.mean=TRUE, archm = TRUE), distribution.model=dist)
fit.m  = ugarchfit(data=r, spec=spec.m)


IC <- cbind(infocriteria(fit1), infocriteria(fit.e), infocriteria(fit.gjr), infocriteria(fit.m))
colnames(IC) <- c("GARCH", "eGARCH", "gjrGARCH","GARCH-in-mean")
kable(IC,digits=3)

archM = TRUE

# best model
spec2 = ugarchspec(variance.model=list(model=mod, garchOrder=pq), 
                   mean.model=list(armaOrder=PQ, include.mean=TRUE, archm = archM), distribution.model=dist)
fit2  = ugarchfit(data=r, spec=spec2)


# 8. Analitical VaR/ES for h=1 from the best GARCH model
########################################################

p <- 0.05                        # tolerance level
v <- fit2@fit$coef["shape"]      # shape of t-distribution

# help(ugarchforecast)
fct2 <- ugarchforecast(fit2,data=r, n.ahead = 1)
sig <- sigma(fct2)
mu  <- fitted(fct2)

qf   <- function(x) qdist("std", p=x, shape=v)
#qf   <- function(x) qdist("norm", p=x)
VaRgarch <- mu + sig*qdist("std", p, shape=v)
ESgarch  <- mu + sig*(1/p * integrate(qf, 0, p)$value)  

# 9. Monte Carlo VAR/ES for longer horizons for the best GARCH model
####################################################################

H    = 10
Nsim = 10000
# help(ugarchsim)
sim2 <- ugarchsim(fit2, n.sim = H, n.start = 0, m.sim = Nsim, startMethod = "sample")
#names(sim1@simulation)
Rdraws  <- fitted(sim2)
temp   <- RdrawsToVaRES(RDraws,p)
VaRHgarch <- temp$VaR
ESHgarch <- temp$ES

# 10. Saving results for the next meetings
##############################################
VaRH <- cbind(VaRHnorm, VaRHt,VaRHhs,VaRHgarch)
ESH  <- cbind(ESHnorm, ESHt,ESHhs,ESHgarch)

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

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