Prepare Data:

rm(list = ls())
library("rugarch")
library("fGarch")
library("MASS")
library("ismev")
# library("lmom")
# library("QRM")
library("skewt")
library("zoo")
library("plotrix")

source("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Rfns.R")
setwd("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500_new/")
##=========================================================
## Backtesting CoVaR
##=========================================================

dat = fread("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/Ret_data.csv") # raw data file (1962/07/02 - 2017/07/31)
dat = dat[,`:=`(RET=as.numeric(as.character(RET)), SP500 = as.numeric(as.character(SP500)))]
dat = dat[!is.na(RET)]
dat = unique(dat, by = c("TICKER", "date") )
setnames(dat, "RET", "r")
r.KO = dat[TICKER == 'KO', list(date, TICKER, r, SP500)]
r.GE = dat[TICKER == 'GE', list(date, TICKER, r, SP500)]
r.BA = dat[TICKER == 'BA', list(date, TICKER, r, SP500)]
r.DIS = dat[TICKER == 'DIS', list(date, TICKER, r, SP500)]

dat.indices <- unique(data.frame(dat$date,dat$Ret_M,dat$SP500))
colnames(dat.indices)<-c("date","retM","sp500")
rM <- zoo(x=dat.indices$retM,order.by = strptime(dat.indices$date,"%Y%m%d"))
rSP<- zoo(x=dat.indices$sp500,order.by = strptime(dat.indices$date,"%Y%m%d"))

load("CoVaR.RDATA")
covar.BA <- out$CoVaR.BA
covar.DIS <- out$CoVaR.DIS
covar.KO <- out$CoVaR.KO
# covar.GE <- out$CoVaR.GE

Out-of-sample Test

## identify subsample where L[t+1]^j > VaR[t+1]^j
n=dim(covar.BA)[1] 
end = n+1

var.BA <- covar.BA[, 25:26]# out$stEVT[,2:3]#
loss.BA <- -r.BA$r[(end-n):(end-1)] 

var.DIS <- covar.DIS[, 25:26]#out$stEVT[,2:3]#
loss.DIS <- -r.DIS$r[(end-n):(end-1)] 

var.KO <- covar.KO[, 25:26]#out$stEVT[,2:3]
loss.KO <- -r.KO$r[(end-n):(end-1)]

# var.GE <- out.p$stEVT[,2:3]
# loss.GE <- -r.GE$r[(end-n):(end-1)] 
# market index data
yI <- as.vector(-rSP[(end-n):(end-1)]) * 100
####
y <- loss.DIS * 100
rmat <- var.DIS
covar <- covar.DIS

n_end = dim(var.DIS)[1]
n_start = 1
y = y[n_start:n_end]
rmat = rmat[n_start:n_end,]
covar = covar[n_start:n_end,]
yI = yI[n_start:n_end]
n = length(y)

We check the performance of predicted VaR:

k=dim(rmat)[2] # number of methods/cases
ymat = matrix(rep(y, k), ncol = k, byrow=FALSE)
hmat = (ymat > rmat) 

apply(hmat,2,"sum",na.rm=T)/n*100
## [1] 5.693977 1.301908

We check the performance of predicted CoVaR:

r.evt.sub95 <- subset(covar[,1],hmat[,1]) 
r.evt.sub99 <- subset(covar[,2],hmat[,2])
r.emp.sub95 <- subset(covar[,3],hmat[,1]) 
r.emp.sub99 <- subset(covar[,4],hmat[,2])
r.fp.sub95 <-  subset(covar[,5],hmat[,1])
r.fp.sub99 <-  subset(covar[,6],hmat[,2])
r.AFG.sub95 <-  subset(covar[,7],hmat[,1])
r.AFG.sub99 <-  subset(covar[,8],hmat[,2])

## Checking % violations of CoVaR by the index losses
yI.sub95 <- as.numeric(subset(yI, hmat[,1]))
yI.sub99 <- as.numeric(subset(yI, hmat[,2]))

# Remove NA values
r.sub95 = rbind(r.evt.sub95,r.emp.sub95, r.fp.sub95)
ind.na95 = which( (is.na(apply(r.sub95,2,sum))) )

r.sub99 = rbind(r.evt.sub99,r.emp.sub99, r.fp.sub99)
ind.na99 = which((is.na(apply(r.sub99,2,sum))) )

if (length(ind.na95)>0){
  r.evt.sub95 <- r.evt.sub95[-ind.na95]
  r.emp.sub95 <- r.emp.sub95[-ind.na95]
  r.fp.sub95 <- r.fp.sub95[-ind.na95]
  r.AFG.sub95 <- r.AFG.sub95[-ind.na95]
  yI.sub95 <- yI.sub95[-ind.na95]
}
if (length(ind.na99)>0){
  r.evt.sub99 <- r.evt.sub99[-ind.na99]
  r.emp.sub99 <- r.emp.sub99[-ind.na99]
  r.fp.sub99 <- r.fp.sub99[-ind.na99]
  r.AFG.sub99 <- r.AFG.sub99[-ind.na99]
  yI.sub99 <- yI.sub99[-ind.na99]
}



m95 <- length(yI.sub95) #sub-sample size for q=0.95
m99 <- length(yI.sub99) #sub-sample size for q=0.99

round(c(sum(yI.sub95>r.evt.sub95)/m95*100,sum(yI.sub95>r.emp.sub95)/m95*100, sum(yI.sub95>r.fp.sub95)/m95*100, sum(yI.sub95>r.AFG.sub95)/m95*100, sum(yI.sub99>r.evt.sub99)/m99*100, sum(yI.sub99>r.emp.sub99)/m99*100, sum(yI.sub99>r.fp.sub99)/m99*100, sum(yI.sub99>r.AFG.sub99)/m99*100 ),2)
## [1] 14.73 10.17  5.76 15.26  0.60 20.24  1.19  2.38
cct <- matrix(nrow=4, ncol=2) # simple & general CCT

tmp <- cct.2s.VaR(x=yI.sub95, r=r.evt.sub95,lev=0.95); cct[1,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub95, r=r.emp.sub95,lev=0.95); cct[2,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub95, r=r.fp.sub95,lev=0.95); cct[3,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub95, r=r.AFG.sub95,lev=0.95); cct[4,]<-c(tmp$pv.avg,tmp$pv.cond)
round(cct,3)
##       [,1] [,2]
## [1,] 0.000 0.00
## [2,] 0.000 0.00
## [3,] 0.375 0.37
## [4,] 0.000 0.00
tmp <- cct.2s.VaR(x=yI.sub99, r=r.evt.sub99,lev=0.99); cct[1,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub99, r=r.emp.sub99,lev=0.99); cct[2,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub99, r=r.fp.sub99,lev=0.99); cct[3,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.2s.VaR(x=yI.sub99, r=r.AFG.sub99,lev=0.99); cct[4,]<-c(tmp$pv.avg,tmp$pv.cond)
round(cct,3)
##       [,1]  [,2]
## [1,] 0.496 0.487
## [2,] 0.000 0.000
## [3,] 0.820 0.357
## [4,] 0.242 0.066
## ONE-SIDED TEST of super-calibration as in Basel Accord
cct.1s <- matrix(nrow=4, ncol=2) # simple & general one-sided CCT of super-calibration
tmp <- cct.1s.VaR(x=yI.sub95, r=r.evt.sub95,lev=0.95); cct.1s[1,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub95, r=r.emp.sub95,lev=0.95); cct.1s[2,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub95, r=r.fp.sub95,lev=0.95); cct.1s[3,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub95, r=r.AFG.sub95,lev=0.95); cct.1s[4,]<-c(tmp$pv.avg,tmp$pv.cond)
round(cct.1s,3)
##       [,1] [,2]
## [1,] 0.000 0.00
## [2,] 0.000 0.00
## [3,] 0.188 0.47
## [4,] 0.000 0.00
cct.1s <- matrix(nrow=4, ncol=2) # simple & general one-sided CCT of super-calibration
tmp <- cct.1s.VaR(x=yI.sub99, r=r.evt.sub99,lev=0.99); cct.1s[1,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub99, r=r.emp.sub99,lev=0.99); cct.1s[2,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub99, r=r.fp.sub99,lev=0.99); cct.1s[3,]<-c(tmp$pv.avg,tmp$pv.cond)
tmp <- cct.1s.VaR(x=yI.sub99, r=r.AFG.sub99,lev=0.99); cct.1s[4,]<-c(tmp$pv.avg,tmp$pv.cond)
round(cct.1s,3)
##       [,1]  [,2]
## [1,] 0.752 1.000
## [2,] 0.000 0.000
## [3,] 0.410 0.615
## [4,] 0.121 0.234
## Comparative backtesting
## q=.95
smat <- matrix(nrow=4, ncol=1)
smat[1,1] <- mean(sfVaR(r=r.evt.sub95,x=yI.sub95,a=0.95,h=1))
smat[2,1] <- mean(sfVaR(r=r.emp.sub95,x=yI.sub95,a=0.95,h=1))
smat[3,1] <- mean(sfVaR(r=r.fp.sub95,x=yI.sub95,a=0.95,h=1))
smat[4,1] <- mean(sfVaR(r=r.AFG.sub95,x=yI.sub95,a=0.95,h=1))
apply(smat,2,"rank")
##      [,1]
## [1,]    3
## [2,]    1
## [3,]    2
## [4,]    4
## q=.99
smat <- matrix(nrow=4, ncol=1)
smat[1,1] <- mean(sfVaR(r=r.evt.sub99,x=yI.sub99,a=0.99,h=1))
smat[2,1] <- mean(sfVaR(r=r.emp.sub99,x=yI.sub99,a=0.99,h=1))
smat[3,1] <- mean(sfVaR(r=r.fp.sub99,x=yI.sub99,a=0.99,h=1))
smat[4,1] <- mean(sfVaR(r=r.AFG.sub99,x=yI.sub99,a=0.99,h=1))
apply(smat,2,"rank")
##      [,1]
## [1,]    2
## [2,]    1
## [3,]    4
## [4,]    3
## TLM
rm=1
tlm95 <- TLMfn(r=rbind(r.evt.sub95,r.emp.sub95,r.fp.sub95,r.AFG.sub95), y=yI.sub95, rm=rm, lev=0.95,h=1)
plotTLM(tlm95$TLM10, rm=1, lev=0.90,method=c("evt","emp","st","AFG"))

tlm99 <- TLMfn(r=rbind(r.evt.sub99,r.emp.sub99,r.fp.sub99,r.AFG.sub99), y=yI.sub99, rm=rm, lev=0.99,h=1)
plotTLM(tlm99$TLM10, rm=1, lev=0.90,method=c("evt","emp","st","AFG"))

Plots

Time series plot of SP500 losses and its CoVaRs from EVT and FP methods

plot(covar[,1], type = "l", col = "blue",  lty = 3, ylab = "Return")
points(covar[,5], type = "l", col = "black", lty = 1)
points(yI, type = "l", col = "red", lty = 2)
legend("topright", legend = c("EVT", "FP", "Loss"), col = c("blue", "black", "red"), lty = c(3,2,1))

Subsample of SP500 losses and CoVaR when individual stock is above 99% VaR

plot(r.fp.sub99, type = "l", col = "blue", ylim =  c(-3,25), lty = 3, ylab = "Return")
points(yI.sub99, type = "l", col = "black", lty = 1)
points(r.evt.sub99, type = "l", col = "red", lty = 2)
legend("topright", legend = c("FP", "EVT", "Loss"), col = c("blue", "red", "black"), lty = c(3,2,1))