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

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.01                             # tolerance level
backN    <- 250                        # backtest sample
w_length <- 1000                       # rolling window length

###############################################################
# Traffic lights method (Bazel II); for p=0.01         #
###############################################################

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


varHS <- lag(varHS, -1)           # lags so that VaR(t+1) = quantile(t)
rr    <- r[index(varHS)]          # rr - realized returns
etaHS  <- tail(rr<=varHS, backN)  # VaR violations

# share of empirical violations (should be close to p)
nHS  <- sum(etaHS); nHS      # number of exceedances
piHS <- mean(etaHS); piHS    # share of exceedances

VaRplot(alpha=p, actual=tail(rr,backN), VaR=tail(varHS,backN))
title(main="Historical simulation: VaR violations")


## Parametric models: normal and t-Student

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

# VaR
varN <- MAmean + MAstd*qdist("norm", p);         varN <- lag(varN, -1)  
# VaR violations
etaN  <- tail(rr<=varN, backN)
# shares / number of violations
nN = sum(etaN); piN = mean(etaN)
# plot
VaRplot(alpha=p, actual=tail(rr,backN), VaR=tail(varN,backN))
title(main="Normal distribution: VaR violations")

# t-Student
varT <- MAmean + MAstd*qdist("std", p, shape=5); varT <- lag(varT, -1)
etaT  <- tail(rr<=varT, backN)
nT = sum(etaT); piT = mean(etaT)
VaRplot(alpha=p, actual=tail(rr,backN), VaR=tail(varT,backN))
title(main="t-Student distribution: VaR violations")

# Return to line 23 and change backtesting settings into  
backN    <- 1000                        # backtest sample
w_length <- 250                         # rolling window length
# Describe the changes?

## Sequential backtesting (number of violations for backN=250)
## Compare to "traffic ligts" methods
## This part needs settings from line 72-73


eta <- etaHS # etaN, etaT
var <- varHS # varN, varT

# rollling number of violations in window of 250 observations
empEx <- rollapply(eta, width=250, sum, align="right")

par(mfrow=c(2,1),oma=c(1,1,1,1),mar=c(1,1,1,1))
plot(c(zoo(NA, head(tail(index(eta), backN), 249)), empEx), xlab="",ylab="", xaxt="n", main="number of violations in window of 250 obs")
abline(h=4,col="green",lty=2)
abline(h=9,col="orange",lty=2)
# Rolling traffic light method
TRplot(actual=tail(rr,backN), VaR=tail(var,backN))

# Traffic lights method explanation
# Binomial distribution
require(binom)
require(knitr)
p       <- 0.01     # tolerance level
backN   <- 250      # backtest sample
lex     <- 0:12     # number of violations

Table = cbind(dbinom(lex,backN,p),pbinom(lex,backN,p))
colnames(Table) = c("pdf","cdf"); rownames(Table) = lex
kable(Table,digits=3)