CoVaR – BA

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/")
##=========================================================
## 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")
# attributes(out)
##columns are 95% CoVaR by EVT method, 99% CoVaR by EVT method and 95% CoVaR by empirical method
covar.KO <- out$CoVaR.KO
covar.BA <- out$CoVaR.BA
covar.DIS <- out$CoVaR.DIS
# covar.GE <- out$CoVaR.GE
## identify subsample where L[t+1]^j > VaR[t+1]^j
n=dim(covar.KO)[1]
end = nrow(r.KO) + 1
load("outKO.RDATA")
var.KO <- out.p$stEVT[,2:3]
loss.KO <- -r.KO$r[(end-n):(end-1)]
load("outBA.RDATA")
var.BA <- out.p$stEVT[,2:3]
loss.BA <- -r.BA$r[(end-n):(end-1)]
load("outDIS.RDATA")
var.DIS <- out.p$stEVT[,2:3]
loss.DIS <- -r.DIS$r[(end-n):(end-1)]
# covar.GE <- out$CoVaR.GE

Out-of-sample Test

We calcualte one-day ahead predicted VaR and CoVaR using past 500 days (rolloing window) and then match with loss series. Specifically, we use data from day 1 to 500 to fit Garch model and get residual series and predicted mu and sigmal. Then we calculate one-day ahead predicted VaR and CoVaR, that is day 501 VaR and CoVaR.

## identify subsample where L[t+1]^j > VaR[t+1]^j
n=dim(covar.KO)[1]
end = nrow(r.KO) + 1
setwd("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500/")
load("outKO.RDATA")
var.KO <- out.p$stEVT[,2:3]
loss.KO <- -r.KO$r[(end-n):(end-1)]
load("outBA.RDATA")
var.BA <- out.p$stEVT[,2:3]
loss.BA <- -r.BA$r[(end-n):(end-1)]
load("outDIS.RDATA")
var.DIS <- out.p$stEVT[,2:3]
loss.DIS <- -r.DIS$r[(end-n):(end-1)]

y <- loss.BA
rmat <- var.BA
covar <- covar.BA
yI <- as.vector(-rSP[(end-n):(end-1)])

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.147774 1.249532

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])
## 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]
  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]
  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.sub99>r.evt.sub99)/m99*100, sum(yI.sub99>r.emp.sub99)/m99*100, sum(yI.sub99>r.fp.sub99)/m99*100),2)

[1] 12.84 10.30  8.36  3.07 20.86  4.29

The EVT method does not perform well for 95% CoVaR estimation but relatively good for 99% CoVaR estimation. Simple CCT is consistent with this result, but comparative backtesting and TLM suggests that FP is slightly better than EVT. In summary, EVT and FP is significantly better than empirical method at 10% level for estimation 99% CoVaR, but EVT and FP are quite similar.

cct <- matrix(nrow=3, 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)
round(cct,3)

      [,1]  [,2]
[1,] 0.000 0.000
[2,] 0.000 0.000
[3,] 0.002 0.006

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)
round(cct,3)

      [,1]  [,2]
[1,] 0.129 0.298
[2,] 0.000 0.000
[3,] 0.041 0.120

## ONE-SIDED TEST of super-calibration as in Basel Accord
cct.1s <- matrix(nrow=3, 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)
round(cct.1s,3)

      [,1]  [,2]
[1,] 0.000 0.000
[2,] 0.000 0.000
[3,] 0.001 0.003

cct.1s <- matrix(nrow=3, 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)
round(cct.1s,3)

      [,1]  [,2]
[1,] 0.064 0.193
[2,] 0.000 0.000
[3,] 0.020 0.052

## Comparative backtesting
## q=.95
smat <- matrix(nrow=3, 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))
apply(smat,2,"rank")

     [,1]
[1,]    3
[2,]    1
[3,]    2

## q=.99
smat <- matrix(nrow=3, 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))
apply(smat,2,"rank")

     [,1]
[1,]    1
[2,]    3
[3,]    2

## TLM
rm=1
tlm95 <- TLMfn(r=rbind(r.evt.sub95,r.emp.sub95,r.fp.sub95), y=yI.sub95, rm=rm, lev=0.95,h=1)
plotTLM(tlm95$TLM10, rm=1, lev=0.95,method=c("evt","emp","st"))

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

In-Sample Test

However, we can show that movement in CoVaR is driven by the predict sigmal from Garch model, and the predicted sigmal shares the similar pattern of lagged loss series. The figure blow plots subsample with 100 days for SP500 loss series. As shown in the figure, the peaks of sigma_SP500 are in line with CoVaR from EVT and FP, but are one-period head of loss series (see x-axis 80-100). This result comes from the fact that predicted sigma is heavily affected by the most recent data. When there is a high reutrn (or loss) today, the predicted volatility will be very high tomorrow.

plot(covar[5800:5900,1], type = "l", col = "blue", ylim =  c(-0.1,0.2), lty = 3, ylab = "Return")
points(covar[5800:5900,5], type = "l", col = "black", lty = 1)

points(yI[5800:5900], type = "l", col = "red", lty = 2)
setwd("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500/")
load("outSP500.RDATA")
points(out.p$sigt.st[5800:5900], type = "l", col = 'green', lty = 4)

legend("topright", legend = c("EVT", "FP", "SP500","sigma_SP500"), col = c("blue", "black", "red","green"), lty = c(3,1,2,4))

Similarly, if we plot subsample in which BA loss is beyond 99% VaR, the peaks of SP00 and CoVaR from different methods are not aligned.

plot(r.fp.sub99, type = "l", col = "blue", ylim =  c(-0.03,0.3), 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))

If we shift Loss series one-day ahead (or shift CoVaR one day one-day backward), then the peaks are aligned. Essentially, this is in-sample analysis. We compare CoVaR estiamted using day 1 to day 500 data with day 500 loss.

## identify subsample where L[t+1]^j > VaR[t+1]^j
n=dim(covar.KO)[1]
end = nrow(r.KO)
loss.BA <- -r.BA$r[(end-n):(end-1)]
y <- loss.BA
yI <- as.vector(-rSP[(end-n):(end-1)])
plot(covar[5800:5900,1], type = "l", col = "blue", ylim =  c(-0.1,0.2), lty = 3, ylab = "Return")
points(covar[5800:5900,5], type = "l", col = "black", lty = 1)

points(yI[5800:5900], type = "l", col = "red", lty = 2)
setwd("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Natalia/out500/")
load("outSP500.RDATA")
points(out.p$sigt.st[5800:5900], type = "l", col = 'green', lty = 4)

legend("topright", legend = c("EVT", "FP", "SP500","sigma_SP500"), col = c("blue", "black", "red","green"), lty = c(3,1,2,4))

Now, we repeat all the analysis for in-sample CoVaR.

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] 3.8982417 0.5836139

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])
## 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]
  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]
  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.sub99>r.evt.sub99)/m99*100, sum(yI.sub99>r.emp.sub99)/m99*100, sum(yI.sub99>r.fp.sub99)/m99*100),2)

[1] 6.51 2.37 4.14 2.74 1.37 4.11

From above numbers, EVT looks work quite well for in-sample analysis. However, comparative backtesting and TLM seem to suggest that EVT is worse than empirical method. I do not understand why we get different resutls from different tests.

cct <- matrix(nrow=3, 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)
round(cct,3)

      [,1]  [,2]
[1,] 0.169 0.143
[2,] 0.000 0.000
[3,] 0.333 0.163

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)
round(cct,3)

      [,1]  [,2]
[1,] 0.365 0.239
[2,] 0.786 0.581
[3,] 0.186 0.090

## ONE-SIDED TEST of super-calibration as in Basel Accord
cct.1s <- matrix(nrow=3, 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)
round(cct.1s,3)

      [,1]  [,2]
[1,] 0.085 0.254
[2,] 1.000 1.000
[3,] 0.834 1.000

cct.1s <- matrix(nrow=3, 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)
round(cct.1s,3)

      [,1]  [,2]
[1,] 0.183 0.333
[2,] 0.393 1.000
[3,] 0.093 0.190

## Comparative backtesting
## q=.95
smat <- matrix(nrow=3, 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))
apply(smat,2,"rank")

     [,1]
[1,]    2
[2,]    1
[3,]    3

## q=.99
smat <- matrix(nrow=3, 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))
apply(smat,2,"rank")

     [,1]
[1,]    2
[2,]    1
[3,]    3

## TLM
rm=1
tlm95 <- TLMfn(r=rbind(r.evt.sub95,r.emp.sub95,r.fp.sub95), y=yI.sub95, rm=rm, lev=0.95,h=1)
plotTLM(tlm95$TLM10, rm=1, lev=0.95,method=c("evt","emp","st"))

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

plot(r.fp.sub99, type = "l", col = "blue", ylim =  c(-0.03,0.4), 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))

plot(yI, ylab="", ylim=c(-.12,0.4), col="gray",xlab="", main="BA", cex.main=0.8, lty = 1)
lines(time(yI),covar.BA[,1], col="red", lty = 2)

lines(time(yI),covar.BA[,5], col="black", lty = 3)
lines(time(yI),covar.BA[,3], col="green", lty = 4)

legend("topleft", legend = c("SP500", "EVT", "FP", "HS"), col = c("grey", "red", "black", "green"), lty = c(1,2,3,4))