rm(list = ls())
require(zoo)
require(rugarch)
require(car)
source("MRFzR_FunkcjeBlok2.R")

source("LoadFundData.R")
# source("LoadWigData.R")

par(mfrow=c(2,1), cex = 0.7, bty="l")
plot(P, main="Level", xlab="", ylab="")
plot(r, main="returns", xlab="", ylab="")

#####################################
# VaR model settings                #
#####################################

T  <- length(r)                        # full sample, nobs
N  <- 1250                             # we shorten the sample to the last 5 years
r  <- tail(r,N)                        # log-returns
R  <- coredata(r)
p  <- 0.05                             # tolerance level
backN    <- 250                        # backtest sample
w_length <- 1000                       # rolling window length

################################################
### Backtesting for ES - McNeil and Frey test  #
################################################

## Historical simulation
##########################################
varHS <- rollapply(r, width=w_length,
                   function(w){quantile(w,p)},
                   by=1, align="right")

esHS  <- rollapply(r, width=w_length,
                  function(w){mean(w[w < quantile(w,p)])},
                  by=1, align="right")

varHS <- lag(varHS, -1)           # lags so that VaR(t+1) = quantile(t)
esHS  <- lag(esHS, -1)

rr     <- r[index(varHS)]          # rr - realized returns
etaHS  <- tail(rr<=varHS, backN)   # VaR violations

# Kupiec and Cristoffersen test for VaR
VaRTest(alpha=p, actual=coredata(rr), VaR=coredata(varHS))

# McNeil and Frey tests
mean(varHS[etaHS])
mean(esHS[etaHS])
mean(rr[etaHS])
ESVaRplot(alpha=p, actual=rr, VaR=varHS, ES=esHS); title(main="Historical simulation")

temp      <- rr[etaHS]-esHS[etaHS]
stat_MF   <- mean(temp)/ (sd(temp)/sqrt(length(temp)))
prob_MF   <- 1- pnorm(abs(stat_MF))
prob_MF 

# The same in rugarch
ESTest(alpha=p, actual=coredata(rr), ES=coredata(esHS), VaR=coredata(varHS))

##########################
## Parametric models     #
##########################

MAmean  <- rollapply(r, width=w_length, mean, by=1, align="right")
MAstd   <- rollapply(r, width=w_length, sd,   by=1, align="right")

# Normal distribution
# VaR
varN <- MAmean + MAstd*qdist("norm", p);          
# VaR violations
etaN <- tail(rr<=varN, backN)
# ES
qf   <- function(x) qdist("norm", p=x)
esN  <- MAmean + MAstd*(1/p * integrate(qf, 0, p)$value)

varN <- lag(varN, -1)           
esN  <- lag(esN, -1)
etaN <- tail(rr<=varN, backN)  

ESVaRplot(alpha=p, actual=rr, VaR=varN, ES=esN); title(main="Normal distribution")
ESTest(alpha=p, actual=coredata(rr), ES=coredata(esN), VaR=coredata(varN))

# t-Student 
# VaR
varT <- MAmean + MAstd*qdist("std", p, shape=5)         
# VaR violations
etaT <- tail(rr<=varT, backN)
# ES
qf       <- function(x) qdist("std", p=x, shape=5)
esT  <- MAmean + MAstd*(1/p * integrate(qf, 0, p)$value)

varT <- lag(varT, -1)           
esT  <- lag(esT, -1)
etaT <- tail(rr<=varT, backN)  

ESVaRplot(alpha=p, actual=rr, VaR=varT, ES=esT); title(main="t-Student distribution")
ESTest(alpha=p, actual=coredata(rr), ES=coredata(esT), VaR=coredata(varT))

#################
## EWMA model   #
#################

lambda <- 0.94

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")
                      #fixed.pars=list(alpha1=1-lambda, omega=0,shape=5), distribution.model = "std")

## quantile function
qf     <- function(x) qdist("norm", p=x)
#qf     <- function(x) qdist("std", p=x, shape=5)
eqt <-  integrate(qf, 0, p)$value

varesEWMA <- rollapply(r, width=w_length,
                        function(w) {
                            frc <- ugarchforecast(EWMAspec,data=w, n.ahead=1)
                            var <- quantile(frc, p)
                            sigma <- sigma(frc)
                            es <- sigma*eqt/p
                            return(c(var, es))
                        },
                        by=1, align="right")
varesEWMA <- lag(varesEWMA, -1)
varEWMA <- varesEWMA[,1]
esEWMA  <- varesEWMA[,2]
etaEWMA <- (rr<varEWMA)

ESVaRplot(alpha=p, actual=rr, VaR=varEWMA, ES=esEWMA); title(main="Model EWMA")
ESTest(alpha=p, actual=coredata(rr), ES=coredata(esEWMA), VaR=coredata(varEWMA))

#######################
## GARCH  model #
#######################

GARCHspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
                       variance.model=list(model="sGARCH", garchOrder=c(1,1)),
                       distribution.model = "norm")
                       #fixed.pars=list(shape=5), distribution.model = "std")

qf     <- function(x) qdist("norm", p=x)
#qf     <- function(x) qdist("std", p=x, shape=5)
eqt <-  integrate(qf, 0, p)$value

varesGARCH <- rollapply(r, width=w_length,
                           function(w) {
                             fit <- ugarchfit(data=w, spec=GARCHspec, solver="hybrid")
                             frc <- ugarchforecast(fit, n.ahead=1)
                             var <- quantile(frc, p)
                             sigma <- sigma(frc)
                             mu <- fitted(frc)
                             es <- mu + sigma*eqt/p
                             return(c(var, es))
                           },
                           by=1, align="right")
varesGARCH <- lag(varesGARCH, -1)
varGARCH <- varesGARCH[,1]
esGARCH  <- varesGARCH[,2]
etaGARCH <- (rr<varGARCH)

ESVaRplot(alpha=p, actual=rr, VaR=varGARCH, ES=esGARCH); title(main="GARCH model")
ESTest(alpha=p, actual=coredata(rr), ES=coredata(esGARCH), VaR=coredata(varGARCH))


#############################################################
#        Backtesting entire distribution                    #
#                   Berkowitz test                          #
#############################################################

###########################
# Parametric model       #
###########################

MAmean  <- rollapply(r, width=w_length, mean, by=1, align="right"); MAmean=lag(MAmean,-1)
MAstd   <- rollapply(r, width=w_length, sd,   by=1, align="right"); MAstd=lag(MAstd,-1)

## rstar must be N(0,1))
rstar <- (rr-MAmean)/ MAstd
# PIT~U(0,1)
PIT   <- pdist("norm",rstar)
# PIT   <- pdist("std",rstar,shape=5) # t-Student
rstar <- qnorm(PIT)

regB <- lm(rstar[-1]~lag(rstar,-1))
summary(regB)
linearHypothesis(regB, hypothesis.matrix=diag(1,2), rhs=c(0,0), test="Chisq")

# the same in rugarch (test LR instead of LM)
testDist <- BerkowitzTest(data = rstar, lags = 1, significance = 0.05)
testDist
testDist$LRp 
testDist$H0  
testDist$Decision  

##########
## GARCH #
##########

GARCHspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
                         variance.model=list(model="sGARCH", garchOrder=c(1,1)),
                         fixed.pars=list(shape=5), distribution.model = "norm")
                         # fixed.pars=list(shape=5), distribution.model = "std")

temp <- rollapply(r, width=w_length,
                        function(w) {
                            fit <- ugarchfit(data=w, spec=GARCHspec, solver="hybrid")
                            frc <- ugarchforecast(fit, n.ahead=1)
                            sigma <- sigma(frc)
                            mu <- fitted(frc)
                            return(c(mu, sigma))
                        },
                        by=1, align="right")
temp <- lag(temp, -1)
MAmean = temp[,1]; MAstd = temp[,2];
rstar <- (rr-MAmean)/ MAstd
# PIT~U(0,1)
PIT   <- pdist("norm",rstar)
# PIT   <- pdist("std",rstar,shape=5) # t-Student
rstar <- qnorm(PIT)

BerkowitzTest(data = rstar, lags = 1, significance = 0.05)