################################################
# 0. Simulate the bt input                        #
################################################
x<-rnorm(250) #realised portfolio P&L
var1<- rep(-qnorm(0.01),250) # projected risk (estimated or true) - true VAR 1% is close to true ES 2.5%
y<-x+var1 # secured position [Realised portfolio P&L + projected portfolio risk]

###IMPORTANT NOTE: If risk profile is changing, or the portfolio size is changing one should normalize it:
## instead of y=P&L+RISK one can set y=(P&L+RISK)/RISK, i.e. so called 'overshooting ratio' as set in TRIM 
#[we resize position so that the 'estimated' risk for it would be equal to 1]

################################################
# 1. VaR & ES backtest statistics definition   #
################################################
VAR.BT<-function(y){sum(       sort(y) < 0)}
ES.BT<- function(y){sum(cumsum(sort(y))< 0)}

sort(y)[1:10]
VAR.BT(y)
ES.BT(y)
################################################
# 2. Exemplary backtests on theoretical data   #
################################################
x<-  rnorm(250)  #realised P&L sample
var1<- -qnorm(0.01) #true VAR 1% risk
es1 <-  dnorm(qnorm(0.025))/(0.025) #true ES 2.5% risk
var2<-0.9*var1; es2<-0.9*es1 #minor underestimated VAR/ESrisks
var3<-0.8*var1; es3<-0.8*es1 #major underestimated VAR/ESrisks

VAR.BT(x+var1); VAR.BT(x+var2); VAR.BT(x+var3) #VAR standard breach bt
ES.BT(x+es1); ES.BT(x+es2); ES.BT(x+es3) #ES cumulative breach bt

# Traffic-zone specification:
# VAR (250obs): 0-4 Green zone; 5-9 Yellow zone; 10+ Red Zone
# ES  (250obs): 0-11 Green zone; 12-24 Yellow zone; 25+ Red zone

#############################################################
# 3. VAR BT/ES BT Consistency check using previous example  #
#############################################################
z1<-cbind(rep(0,10000),rep(0,10000)); z2<-cbind(rep(0,10000),rep(0,10000)); z3<-cbind(rep(0,10000),rep(0,10000))
#z1,z2,z3 is 10000 x 2 matrix storing T_n and G_n values for each run

for(i in 1:10000){
  x<-rnorm(250)
  #var1,var2,var3 and es1,es2,es3 as before
  z1[i,]<-c(min(VAR.BT(x+var1),14),min(ES.BT(x+es1),37)) #upper threshold if more than 14 or 37 let set to 14 or 37 - it's in red zone anyway
  z2[i,]<-c(min(VAR.BT(x+var2),14),min(ES.BT(x+es2),37))
  z3[i,]<-c(min(VAR.BT(x+var3),14),min(ES.BT(x+es3),37))
  }
library(MASS)
par(mfrow=c(1,3))
#True risk consistency
  image(kde2d(z1[,1],z1[,2]),col=gray((32:0)/32),xlim=c(-1,15),ylim=c(-1,38),main="Proper risk (10k strong MC)",xlab="VAR BT",ylab="ES BT")
  abline(v=4.5,lty=2);abline(v=9.5,lty=2); abline(h=11.5,lty=2);abline(h=24.5,lty=2)
#Minor underestimation consistency
  image(kde2d(z2[,1],z2[,2]),col=gray((32:0)/32),xlim=c(-1,15),ylim=c(-1,38),main="10% risk underestimation (10k strong MC)",xlab="VAR BT",ylab="ES BT")
  abline(v=4.5,lty=2);abline(v=9.5,lty=2); abline(h=11.5,lty=2);abline(h=24.5,lty=2)
#Major underestimation consistency 
  image(kde2d(z3[,1],z3[,2]),col=gray((32:0)/32),xlim=c(-1,15),ylim=c(-1,38),main="20% risk underestimation (10k strong MC)",xlab="VAR BT",ylab="ES BT")
  abline(v=4.5,lty=2);abline(v=9.5,lty=2)
  abline(h=11.5,lty=2);abline(h=24.5,lty=2)
  

########################################################
# 4. Test statistic distribution check (for normal)    #
########################################################
library(dplyr)
par(mfrow=c(1,1))  
  
## VAR
  #VAR true risk
  plot(y=table(z1[,1]),x=as.numeric(names(table(z1[,1]))),type="h",xlim=c(0,14),lwd=4,main="VAR BT",col="green",xlab="breaches",ylab="nr.runs")
  abline(v=4.5,lty=2);abline(v=9.5,lty=2)
  #adding z2 (mior underestimation)
  lines(y=table(z2[,1]),x=as.numeric(names(table(z2[,1])))+0.2,type="h",lwd=4,col="orange")
  #adding z3 (major underestimation)
  lines(y=table(z3[,1]),x=as.numeric(names(table(z3[,1])))+0.4,type="h",lwd=4,col="red")
  
## ES
  ##ES true risk
  plot(y=table(z1[,2]),x=as.numeric(names(table(z1[,2]))),type="h",xlim=c(0,25),lwd=4,main="ES BT",col="green")
  abline(v=11.5,lty=2);abline(v=24.5,lty=2)
  #adding z2 (mior underestimation)
  lines(y=table(z2[,2]),x=as.numeric(names(table(z2[,2])))+0.2,type="h",lwd=4,col="orange")
  #adding z3 (major underestimation)
  lines(y=table(z3[,2]),x=as.numeric(names(table(z3[,2])))+0.4,type="h",lwd=4,col="red")
  
#############################################################
# 5. Checking the same for t-student (or any other distr)   #
#############################################################
  #It's very easy - you need to slightly modify 3.
  
#define risk specifications for t-student  
df<-5
var1<- -qt(0.01,df=df) #true VAR risk
es1<- dt(qt(0.025,df=df),df=df) / (0.025) * (df+qt(0.025,df=df)^2)/(df-1) #true ES risk
var2<-0.85*var1; es2<-0.85*es1 #minor underestimated VAR/ESrisks
var3<-0.7*var1; es3<-0.7*es1 #major underestimated VAR/ESrisks
  
#Modify input in 3 to sample from t-student, i.e.
#x<-rt(250,df=df)

#You can check other null specification. We did it for multiple other choices including GARCH specifications (e.g. with skew t-student innovations)