rm(list=ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
# packages used
listofpackages <- c("tidyverse","dygraphs", "rugarch", "forecast","dplyr","ellipse","reshape2","ggplot2","xts","xlsx","readxl")

for (j in listofpackages){
  if(sum(installed.packages()[, 1] == j) == 0) {
    install.packages(j)
  }
  library(j, character.only = T)
}

# loading the  databases
load("VaRdata.Rdata")
ggplot(data_simplot, aes(x = Date)) +
  geom_line(aes(y = BA_ret, color = "BA"), size = 1, linetype = "solid") +
  geom_line(aes(y = var_BA_capm, color = "Var CAPM"), size = 1, linetype = "solid") +
  geom_line(aes(y = var_BA_garch, color = "VaR GARCH"), size = 1, linetype = "solid") +
  labs(x = "Date", y = "Returns and VaR") +
  ylim(-0.30, 0.20) +
  theme_minimal() +
  theme(
    legend.position = c(0.15, 0.95),  # Set the legend position (top-left)
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 10),
    plot.title = element_text(size = 12, hjust = 0.5)
  ) +
  scale_color_manual(
    values = c("blue", "green", "red"),
    labels = c("BA", "VaR CAPM", "Var GARCH")
  )
## --------------------------------------------------------------------------------------------------
# VaR tail
alpha <- 0.1 


## --------------------------------------------------------------------------------------------------
violations <- (data_simplot$BA_ret - data_simplot$var_BA_capm) < 0
table(violations)
plot(y = violations*runif(length(violations), min = 0.99, max = 1.01), x = data_simplot$Date, main = "VaR violations", 
     ylab = "violations", xlab = "time") # adding jitter to make sure that adjacent observations don't overlap

## testing
### Unconditional coverage
p <- alpha
T1 <- sum(violations)
T0 <- length(violations) - sum(violations)

x <- T1/(T1+T0) # violations as fraction of sample length, which is also the test statistic estimate

L_p <- (1-p)^T0 * p^T1
L_x <- (1-x)^T0 * x^T1

LR_uc <- -2 * log(L_p/L_x) # test statistic 
critical <- qchisq(p = 0.95, df = 1) # the 5% critical value

LR_uc; critical
LR_uc > critical 

### Independence testing
temp1 <- abs(diff(violations)) # to identify moments of change
temp2 <- violations[2:length(violations)] # to identify the ending points 
T01 <- sum(temp2 * temp1) # those that finish with 1 and had a change
T11 <- sum(temp2 * (1-temp1)) # those that finish with 1 and had no change
T10 <- sum((1 - temp2) * temp1) # those finishing with 0 and having a change
T00 <- sum((1-temp2) * (1-temp1)) # finishing with 0 and no change

xhat <- x # from before
x00 <- T00/(T00 + T01)
x11 <- T11/(T11 + T10)

L_x_different <- x00^T00 * (1-x00)^T01 * (1-x11)^T10 * x11^T11
L_x_equal <- (1-xhat)^T00 * xhat^T01 * (1-xhat)^T10 * xhat^T11

LR_ind <- -2 * log(L_x_equal/L_x_different)
critical <- qchisq(p = 0.95, df = 1) # the 5% critical value
LR_ind; critical
LR_ind > critical