Empirical

Read in Data

rm(list = ls())
source("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Rfns.R")
library("data.table")
library("ggplot2")
library("psych")
library("lintools")
source("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Rfns.R")

path = "/Users/jinyuanzhang/OneDrive - INSEAD/SkewEllipt/Code_Update/CoVaR/"
folders= 2:99
folders = folders[-71]
load(paste(path,"s",as.character(1),"/COVAR.RDATA",sep=""))  
DF = out$df
fit = c(out$df$firm[1],out$fit$EVT)
for(folder in folders){
  load(paste(path,"s",as.character(folder),"/COVAR.RDATA",sep=""))
  DF = rbind(DF, out$df)
  fit = rbind(fit, c(out$df$firm[1],out$fit$EVT))
}
DF = data.table(DF)
fit = data.table(fit)
fit[,  2:7 := lapply(.SD, as.numeric), .SDcols = 2:7]
colnames(fit) = c("firm", "alpha1", "alpha2", "rho", "nu", "sd1", "sd2")
#### Trim data at 0.5% tails
DF[, CoVaR_evt99 := winsor(CoVaR_evt99, 0.005)]
DF[, CoVaR_fp99 := winsor(CoVaR_fp99, 0.005)]
DF[, VaR99 := winsor(VaR99, 0.005)]
DF = DF[!(CoVaR_evt99==max(CoVaR_evt99))]
DF = DF[!(CoVaR_fp99==max(CoVaR_fp99))]
DF = DF[!(VaR99==max(VaR99))]

DF[, N := .N, by = .(firm)]

key <- fread("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Risk Contagion/key.csv")
colnames(key)[1] <- "firm"
key$firm = as.character(key$firm)
DF = merge(DF, key, by = "firm")
fit = merge(fit, key, by = "firm")

report.1 <- function(dat){
  c(length(dat), mean(as.matrix(dat), na.rm=TRUE),
    sqrt(mean( (as.matrix(dat) - mean(as.matrix(dat), na.rm=T))^2, na.rm =T)),
    quantile(as.matrix(dat), na.rm=T, c(0.05,0.25,0.50,0.75,0.95)))
}

Summary statistics for estiamted parameters

report = t(apply(fit[,2:7], 2, report.1))
colnames(report) = c("Count","Mean", "Std", "5%", "25%", "50%", "75%", "95%")
report

##        Count        Mean        Std          5%         25%          50%
## alpha1    98 -0.03314823 0.25036200 -0.08715549 -0.03819582 -0.005212298
## alpha2    98  0.09275409 0.04074092  0.02173147  0.06488612  0.097669435
## rho       98  0.10576065 0.07335300  0.03082969  0.06014987  0.094194268
## nu        98  4.17695280 0.50476527  3.29764704  3.92147810  4.175215573
## sd1       98  0.72740855 0.09604552  0.63929725  0.67956889  0.718663307
## sd2       98  0.82950944 0.10866660  0.69571860  0.81442924  0.835132940
##               75%        95%
## alpha1 0.02852873 0.06840548
## alpha2 0.12330919 0.14861774
## rho    0.13239995 0.23476957
## nu     4.45621180 5.03889515
## sd1    0.74827677 0.82289698
## sd2    0.87643267 0.94451352

Comparison between EVT and FP mthods

##### Empirical estimator for violation: Loss_stock > VaR_stock
hist(DF[, list(sum( (Loss_stock > VaR99)/N )), by = .(firm)]$V1)

##### Empirical estimator for violatoin: Loss_Index > CoVaR_Index given Loss_stock > VaR_stock
DF_CoVaR = DF[Loss_stock > VaR99]
DF_CoVaR[, CoVaR_N := .N, by = .(firm)]
# DF_CoVaR = DF_CoVaR[!(firm %in% c(10562,66157,86868))] #10562 66157 86868
DF_CoVaR = DF_CoVaR[CoVaR_N >= 5]
print(length(unique(DF_CoVaR$firm)))

## [1] 96

### EVT method 
mean(DF_CoVaR[, list(sum((Loss_Index > CoVaR_evt99)/CoVaR_N)), by = .(firm)]$V1)

## [1] 0.06310318

### FP method 
mean(DF_CoVaR[, list(sum((Loss_Index > CoVaR_fp99)/CoVaR_N)), by = .(firm)]$V1)

## [1] 0.03654653

#### LRuc #########
LRuc_pvalue = DF_CoVaR[, list(fp_LRuc = LRuc(0.01, Loss_Index > CoVaR_fp99)$pvalue, evt_LRuc = LRuc(0.01, Loss_Index > CoVaR_evt99)$pvalue), by = .(firm, Group)]
LRuc_pvalue[, `:=`(p10_evt = evt_LRuc > 0.1, p05_evt = evt_LRuc > 0.05, p01_evt = evt_LRuc > 0.01, p10_fp = fp_LRuc > 0.1, p05_fp = fp_LRuc > 0.05, p01_fp = fp_LRuc > 0.01)]

out = apply(LRuc_pvalue[, 5:10],2, sum) #cases with pvalue > 0.1, 0.05, 0.01
out = rbind(out, apply(LRuc_pvalue[Group == "depository", 5:10],2, sum)) 
out = rbind(out, apply(LRuc_pvalue[Group == "other", 5:10],2, sum))
out = rbind(out, apply(LRuc_pvalue[Group == "insurance", 5:10],2, sum))
out = rbind(out, apply(LRuc_pvalue[Group == "broker", 5:10],2, sum))
rownames(out) = c("All", "Depositories", "Others", "Insurance", "Broker-dealers")
count = c(nrow(LRuc_pvalue), nrow(LRuc_pvalue[Group == "depository"]),  nrow(LRuc_pvalue[Group == "other"]),  nrow(LRuc_pvalue[Group == "insurance"]),  nrow(LRuc_pvalue[Group == "broker"]) )
ave_p_fp = c(mean(LRuc_pvalue$fp_LRuc), mean(LRuc_pvalue[Group == "depository"]$fp_LRuc),  mean(LRuc_pvalue[Group == "other"]$fp_LRuc),  mean(LRuc_pvalue[Group == "insurance"]$fp_LRuc),  mean(LRuc_pvalue[Group == "broker"]$fp_LRuc) )
ave_p_evt = c(mean(LRuc_pvalue$evt_LRuc), mean(LRuc_pvalue[Group == "depository"]$evt_LRuc),  mean(LRuc_pvalue[Group == "other"]$evt_LRuc),  mean(LRuc_pvalue[Group == "insurance"]$evt_LRuc),  mean(LRuc_pvalue[Group == "broker"]$evt_LRuc) )
cbind(count, out[,1:3], ave_p_evt, out[,4:6], ave_p_fp)

##                count p10_evt p05_evt p01_evt ave_p_evt p10_fp p05_fp
## All               96      51      56      71 0.2433203     64     70
## Depositories      29      10      12      13 0.1642987     17     19
## Others            26      16      17      24 0.2825238     18     19
## Insurance         32      20      21      28 0.2717759     24     27
## Broker-dealers     9       5       6       6 0.2835152      5      5
##                p01_fp  ave_p_fp
## All                79 0.3213886
## Depositories       21 0.2891741
## Others             23 0.3274383
## Insurance          30 0.3576532
## Broker-dealers      5 0.2787727

#### cct.twosided #####
#### pv.cond
fp_cct =  DF_CoVaR[, (cct.2s.VaR(r = CoVaR_fp99, x = Loss_Index , lev=0.99, hac = FALSE)$pv.cond), by = .(firm, Group)]
colnames(fp_cct)[3] = "fp_cct"
evt_cct =  DF_CoVaR[, (cct.2s.VaR(r = CoVaR_evt99, x = Loss_Index , lev=0.99, hac = FALSE)$pv.cond), by = .(firm, Group)]
colnames(evt_cct)[3] = "evt_cct"
cct_pvalue = merge(fp_cct, evt_cct, by = c("firm","Group"))
cct_pvalue[, `:=`(p10_evt = evt_cct > 0.1, p05_evt = evt_cct > 0.05, p01_evt = evt_cct > 0.01, p10_fp = fp_cct > 0.1, p05_fp = fp_cct > 0.05, p01_fp = fp_cct > 0.01)]

out = apply(cct_pvalue[, 5:10],2, sum) #cases with pvalue > 0.1, 0.05, 0.01
out = rbind(out, apply(cct_pvalue[Group == "depository", 5:10],2, sum)) 
out = rbind(out, apply(cct_pvalue[Group == "other", 5:10],2, sum))
out = rbind(out, apply(cct_pvalue[Group == "insurance", 5:10],2, sum))
out = rbind(out, apply(cct_pvalue[Group == "broker", 5:10],2, sum))
rownames(out) = c("All", "Depositories", "Others", "Insurance", "Broker-dealers")
count = c(nrow(cct_pvalue), nrow(cct_pvalue[Group == "depository"]),  nrow(cct_pvalue[Group == "other"]),  nrow(cct_pvalue[Group == "insurance"]),  nrow(cct_pvalue[Group == "broker"]) )
ave_p_fp = c(mean(cct_pvalue$fp_cct), mean(cct_pvalue[Group == "depository"]$fp_cct),  mean(cct_pvalue[Group == "other"]$fp_cct),  mean(cct_pvalue[Group == "insurance"]$fp_cct),  mean(cct_pvalue[Group == "broker"]$fp_cct) )
ave_p_evt = c(mean(cct_pvalue$evt_cct), mean(cct_pvalue[Group == "depository"]$evt_cct),  mean(cct_pvalue[Group == "other"]$evt_cct),  mean(cct_pvalue[Group == "insurance"]$evt_cct),  mean(cct_pvalue[Group == "broker"]$evt_cct) )
cbind(count, out[,1:3], ave_p_evt, out[,4:6], ave_p_fp)

##                count p10_evt p05_evt p01_evt  ave_p_evt p10_fp p05_fp
## All               96      44      55      67 0.15019290     36     43
## Depositories      29      10      17      23 0.12875793     13     16
## Others            26      12      14      15 0.13040299     11     12
## Insurance         32      19      20      23 0.20008314     11     11
## Broker-dealers     9       3       4       6 0.09904439      1      4
##                p01_fp   ave_p_fp
## All                51 0.12441521
## Depositories       19 0.15319362
## Others             13 0.13049516
## Insurance          13 0.10648460
## Broker-dealers      6 0.07787378

#### pv.avg
fp_cct =  DF_CoVaR[, (cct.2s.VaR(r = CoVaR_fp99, x = Loss_Index , lev=0.99, hac = FALSE)$pv.avg), by = .(firm, Group)]
colnames(fp_cct)[3] = "fp_cct"
evt_cct =  DF_CoVaR[, (cct.2s.VaR(r = CoVaR_evt99, x = Loss_Index , lev=0.99, hac = FALSE)$pv.avg), by = .(firm, Group)]
colnames(evt_cct)[3] = "evt_cct"
cct_pvalue = merge(fp_cct, evt_cct, by = c("firm","Group"))
cct_pvalue[, `:=`(p10_evt = evt_cct > 0.1, p05_evt = evt_cct > 0.05, p01_evt = evt_cct > 0.01, p10_fp = fp_cct > 0.1, p05_fp = fp_cct > 0.05, p01_fp = fp_cct > 0.01)]

out = apply(cct_pvalue[, 5:10],2, sum) #cases with pvalue > 0.1, 0.05, 0.01
out = rbind(out, apply(cct_pvalue[Group == "depository", 5:10],2, sum)) 
out = rbind(out, apply(cct_pvalue[Group == "other", 5:10],2, sum))
out = rbind(out, apply(cct_pvalue[Group == "insurance", 5:10],2, sum))
out = rbind(out, apply(cct_pvalue[Group == "broker", 5:10],2, sum))
rownames(out) = c("All", "Depositories", "Others", "Insurance", "Broker-dealers")
count = c(nrow(cct_pvalue), nrow(cct_pvalue[Group == "depository"]),  nrow(cct_pvalue[Group == "other"]),  nrow(cct_pvalue[Group == "insurance"]),  nrow(cct_pvalue[Group == "broker"]) )
cbind(count, out)

##                count p10_evt p05_evt p01_evt p10_fp p05_fp p01_fp
## All               96      49      56      60     37     43     45
## Depositories      29      14      17      20     12     16     17
## Others            26      14      14      14     12     12     12
## Insurance         32      19      21      21     10     11     11
## Broker-dealers     9       2       4       5      3      4      5

ave_p_fp = c(mean(cct_pvalue$fp_cct), mean(cct_pvalue[Group == "depository"]$fp_cct),  mean(cct_pvalue[Group == "other"]$fp_cct),  mean(cct_pvalue[Group == "insurance"]$fp_cct),  mean(cct_pvalue[Group == "broker"]$fp_cct) )
ave_p_evt = c(mean(cct_pvalue$evt_cct), mean(cct_pvalue[Group == "depository"]$evt_cct),  mean(cct_pvalue[Group == "other"]$evt_cct),  mean(cct_pvalue[Group == "insurance"]$evt_cct),  mean(cct_pvalue[Group == "broker"]$evt_cct) )
cbind(count, out[,1:3], ave_p_evt, out[,4:6], ave_p_fp)

##                count p10_evt p05_evt p01_evt  ave_p_evt p10_fp p05_fp
## All               96      49      56      60 0.14973607     37     43
## Depositories      29      14      17      20 0.11927437     12     16
## Others            26      14      14      14 0.14304709     12     12
## Insurance         32      19      21      21 0.20062208     10     11
## Broker-dealers     9       2       4       5 0.08628608      3      4
##                p01_fp   ave_p_fp
## All                45 0.11431501
## Depositories       17 0.13690477
## Others             12 0.12583795
## Insurance          11 0.09379921
## Broker-dealers      5 0.08118234

## Comparative backtesting
com1 =  DF_CoVaR[, as.list(TLMfn.biv(r = CoVaR_evt99, rS = CoVaR_fp99,  y= Loss_Index , lev=0.99, h = 1)), by = .(firm, Group)]
com0 =  DF_CoVaR[, as.list(TLMfn.biv(r = CoVaR_evt99, rS = CoVaR_fp99,  y= Loss_Index , lev=0.99, h = 0)), by = .(firm, Group)]
colnames(com0)[3:6] = c("TLM01_0", "TLM05_0", "TLM10_0", "Tn_0")
colnames(com1)[3:6] = c("TLM01_1", "TLM05_1", "TLM10_1", "Tn_1")
com = merge(com0, com1, by = c("firm", "Group"))

#### All  ###
red.count <- apply((com[,c(3:5, 7:9)]==0),2,"sum") #number of 0's in each column
yellow.count <- apply((com[,c(3:5, 7:9)]==1),2,"sum")  #number of 1's in each column
green.count <-  apply((com[,c(3:5, 7:9)]==2),2,"sum")  #number of 2's in each column
round(rbind(red.count,yellow.count,green.count)/nrow(com)*100) ## percentage in each category

##              TLM01_0 TLM05_0 TLM10_0 TLM01_1 TLM05_1 TLM10_1
## red.count         20      30      40      12      24      31
## yellow.count      75      62      48      80      68      57
## green.count        5       7      12       7       8      11

#### depository ### 
red.count <- apply((com[Group == "depository",c(3:5, 7:9)]==0),2,"sum") #number of 0's in each column
yellow.count <- apply((com[Group == "depository",c(3:5, 7:9)]==1),2,"sum")  #number of 1's in each column
green.count <-  apply((com[Group == "depository",c(3:5, 7:9)]==2),2,"sum")  #number of 2's in each column
round(rbind(red.count,yellow.count,green.count)/nrow(com[Group == "depository"])*100) ## percentage in each category

##              TLM01_0 TLM05_0 TLM10_0 TLM01_1 TLM05_1 TLM10_1
## red.count         10      17      21       0      10      14
## yellow.count      83      72      66      90      79      76
## green.count        7      10      14      10      10      10

#### other ### 
red.count <- apply((com[Group == "other",c(3:5, 7:9)]==0),2,"sum") #number of 0's in each column
yellow.count <- apply((com[Group == "other",c(3:5, 7:9)]==1),2,"sum")  #number of 1's in each column
green.count <-  apply((com[Group == "other",c(3:5, 7:9)]==2),2,"sum")  #number of 2's in each column
round(rbind(red.count,yellow.count,green.count)/nrow(com[Group == "other"])*100) ## percentage in each category

##              TLM01_0 TLM05_0 TLM10_0 TLM01_1 TLM05_1 TLM10_1
## red.count         15      38      50       8      31      46
## yellow.count      77      50      35      81      54      38
## green.count        8      12      15      12      15      15

#### insurance ### 
red.count <- apply((com[Group == "insurance",c(3:5, 7:9)]==0),2,"sum") #number of 0's in each column
yellow.count <- apply((com[Group == "insurance",c(3:5, 7:9)]==1),2,"sum")  #number of 1's in each column
green.count <-  apply((com[Group == "insurance",c(3:5, 7:9)]==2),2,"sum")  #number of 2's in each column
round(rbind(red.count,yellow.count,green.count)/nrow(com[Group == "insurance"])*100) ## percentage in each category

##              TLM01_0 TLM05_0 TLM10_0 TLM01_1 TLM05_1 TLM10_1
## red.count         34      41      56      31      38      44
## yellow.count      66      59      41      69      62      53
## green.count        0       0       3       0       0       3

#### broker ### 
red.count <- apply((com[Group == "broker",c(3:5, 7:9)]==0),2,"sum") #number of 0's in each column
yellow.count <- apply((com[Group == "broker",c(3:5, 7:9)]==1),2,"sum")  #number of 1's in each column
green.count <-  apply((com[Group == "broker",c(3:5, 7:9)]==2),2,"sum")  #number of 2's in each column
round(rbind(red.count,yellow.count,green.count)/nrow(com[Group == "broker"])*100) ## percentage in each category

##              TLM01_0 TLM05_0 TLM10_0 TLM01_1 TLM05_1 TLM10_1
## red.count         11      11      11       0       0       0
## yellow.count      78      78      56      89      89      67
## green.count       11      11      33      11      11      33

#### Compare scores #####
c(sum((com$Tn_0<0)),sum((com$Tn_1<0)))

## [1] 32 36

Some plots

############## Correlation of CoVaR and VaR ################
corraltion = DF[, .(corr = cor(CoVaR_evt99, VaR99)), by = .(firm)]
hist(corraltion$corr)

#### Plot ########################
DF_ts = DF[, list(CoVaR_EVT=mean(CoVaR_evt99, na.rm=T), CoVaR_FP=mean(CoVaR_fp99, na.rm=T), VaR99 = mean(VaR99, na.rm=T), Loss_stock = mean(Loss_stock, na.rm=T), Loss_DJUSFN = mean(Loss_Index, na.rm=T)), by = "date"]
list=colnames(DF_ts)[2:6]
DF_ts_long <- data.table(reshape(DF_ts, varying = list, times = list,  timevar = "var", v.names = "value",  direction = "long" ))
ggplot(DF_ts_long[var%in%c("CoVaR_EVT", "CoVaR_FP","Loss_DJUSFN")], aes(as.Date.numeric(date,origin = "1970-01-01"), value, color = var)) + geom_line() +  theme(legend.title = element_blank(), legend.position ="bottom") + xlab('Year') + ylab("Losses(%)")