###########################################
# Topic 8                                 # 
# Backtesting                             # 
###########################################

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

require(zoo)
require(xts)
require(MASS)
require(rugarch)
require(forecast)
require(car)
require(knitr)
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)

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


# 2. Rolling VaR/ES 
#######################################

# Normal distribution / vs t-Distribution / vs historical simulations / vs EWMA

p     <- 0.05         # VaR/ES tolerance level
M     <- 250          # number of observations for backtesting
R     <- 500          # rolling window length (here we apply rolling forecasting scheme)


temp     <- NORMfct(r, M, R, p)
ESrnorm  <- temp$ES
VaRrnorm <- temp$VaR

temp   <- HSfct(r, M, R, p)
ESrhs  <- temp$ES
VaRrhs <- temp$VaR

temp   <- Tfct(r, M, R, p, v=5)
ESrt   <- temp$ES
VaRrt  <- temp$VaR

temp      <- EWMAfct(r, M, R, p, lambda=0.94)
ESrewma   <- temp$ES
VaRrewma  <- temp$VaR


#temp      <- GARCHfct(r, M, R, p, v=5)
#ESgarch  <- temp$ES
#VaRgarch  <- temp$VaR

#temp     <- DCCGARCHfct(dy,w, M, R, p, v=5)
#ESrdcc   <- temp$ES
#VaRrdcc  <- temp$VaR


realized <- tail(r,M)
VaR      <- VaRrewma
ES       <- ESrewma

plot(realized,main="Portfolio returns")
abline(h=0)
lines(VaR, col="red",lwd=3)
lines(ES,  col="blue",lwd=3)

# graph in rugarch
VaRplot(p, actual=realized, VaR=VaR)
abline(h=0)
lines(ES,  col="blue",lwd=3)

# 3. Backtesting
#################

# expected number of exceedances
require(knitr)
K = 25
x = data.frame(n=1:K,pdf=dbinom(1:K,size=M,prob=p),cdf=pbinom(1:K,size=M,prob=p))
kable(x,digits=3)

# actual number of exceedances
HR     <- as.numeric(realized < VaR) # hit ratio: VaR exceedances
a      <- mean(HR)

# a. Unconditional coverage: Kupiec test
# H0: alpha = p
# alpha - expected value of VaR exceedance share
# see: https://www.value-at-risk.net/backtesting-coverage-tests/
####################################################################

# regression approach, HR[t] = alpha + epsilon[t]
require(car)
HR_reg = lm(HR~1)
summary(HR_reg)
linearHypothesis(HR_reg, hypothesis.matrix=1, rhs=p, test="Chisq")

# Kupiec test
n  = length(HR)
n1 = sum(HR)
n0 = n - n1

LR_uc   = (p/a)^n1 *((1-p)/(1-a))^n0
stat_uc = -2*log(LR_uc)
prob_uc = 1 - pchisq(stat_uc,1, lower.tail=TRUE)
prob_uc

# The same in rugarch package
help(VaRTest)
temp <- VaRTest(alpha=p, coredata(realized), coredata(VaR))
temp$expected.exceed # p*M
temp$actual.exceed   # a*M = n1 
temp$uc.H0           # alpha = p
temp$uc.LRstat       # stat_k
temp$uc.LRp          # prob_k
temp$uc.Decision     

# b. Christofersen independence test) 
# H0: p(HR[t]=1|HR[t-1]=1) = p(HR[t]=1|HR[t-1]=0)      
###########################################################

# Autocorrelation function
require(forecast)
Acf(HR)
Box.test(HR,1)

# regression approach
# HR[t] = alpha + rho HR[t-1] + epsilon[t]
HR0 <- HR[2:M]
HR1 <- HR[1:(M-1)]

HR_reg = lm(HR0~HR1)
summary(HR_reg)
# H0: rho = 0
linearHypothesis(HR_reg, hypothesis.matrix=c(0,1), rhs=0, test="Chisq")

# Christofersen independence test
n00 = sum(!HR1 & !HR0) # no exceedance after no exceedance
n01 = sum(!HR1 & HR0)  # exceedance after no exceedance
n10 = sum(HR1 & !HR0)  # no exceedance after exceedance
n11 = sum(HR1 & HR0)   # exceedance after exceedance

n0  = n00 + n10
n1  = n01 + n11

pi01 = n01 / (n00+n01) # prob. of exceedance after no exceedance 
pi11 = n11 / (n10+n11) # prob. of exceedance after exceedance 
pi   = (n01+n11) / (n00+n01+n10+n11) # podobna wartosc do alpha

LR_ind   = pi^n1*(1-pi)^n0 / (pi01^n01 * pi11^n11 * (1-pi01)^n00  * (1-pi11)^n10)
stat_ind = -2*log(LR_ind)
prob_ind = 1 - pchisq(stat_ind,1)
prob_ind


# c. Christofersen conditional coverage test                       
# H0: p(HR[t]=1|HR[t-1]=1) = p(HR[t]=1|HR[t-1]=0)=p      
#################################################################

# regression approach
HR_reg = lm(HR0~HR1)
summary(HR_reg)
# H0: alpha = ap, rho = 0
linearHypothesis(HR_reg, diag(1,2), c(p,0), "Chisq")

# Christofersen conditional coverage test
pi01 = n01 / (n00+n01) # prob. of exceedance after no exceedance 
pi11 = n11 / (n10+n11) # prob. of exceedance after exceedance 

# LR_cc   = (p/pi01)^n01 * (p/pi11)^n11 * ((1-p)/(1-pi01))^n00 * ((1-p)/(1-pi11))^n10
LR_cc=LR_ind*LR_uc
stat_cc = -2*log(LR_cc)
prob_cc = 1 - pchisq(stat_cc,2)
prob_cc

# The same in rugarch package
temp <- VaRTest(alpha=p, coredata(realized), coredata(VaR))
temp$cc.H0            
temp$cc.LRstat       # stat_cc
temp$cc.LRp          # prob_cc
temp$cc.Decision      

# d. McNeil and Frey test for well calibrated Expected Shartfall  
################################################

# only when VaR was exceeded
ind <- which(HR==1)
mean(VaR[ind])
mean(ES[ind])
mean(realized[ind])

temp      <- realized[ind]-ES[ind]
m         <- mean(temp)
s         <- sd(temp)
N         <- length(temp)
stat_mf   <- mean(temp)/ (sd(temp)/sqrt(N))
prob_mf   <- 1- pnorm(abs(stat_mf))
prob_mf 

# the same in rugarch
temp <- ESTest(alpha = p, realized, ES, VaR)
temp$p.value