source("~/OneDrive - INSEAD/SkewEllipt/Code_Update/Final_function.R")
path = "~/OneDrive - INSEAD/SkewEllipt/Code_Update/CoVaR/"
k = 1
#### Read in Data
data = fread(paste(path,"ALL.csv",sep=""))
DJ = fread(paste(path,"Don.csv",sep=""))
DJ[,date := as.numeric(format(strptime(date,"%m/%d/%Y"),"%Y%m%d"))]
data = rbind(data, DJ)
# data <- data[ date>=20040101]
data <- data[!is.na(PRC)]
Ldata <- data.table(dcast(data,formula=date ~ PERMNO, value.var="PRC", drop = FALSE))
Ldata <- Ldata[!(rowSums(is.na(Ldata[,2:(ncol(Ldata)-1)])) == (ncol(Ldata)-2))]
Ldata[, date:=as.Date(as.character(date), "%Y%m%d")]
Dat = as.timeSeries(Ldata)
R_com=Return.calculate(Dat, method = "log")[-1,] * 100
num = ncol(R_com)
ndays = length(unique(Ldata$date)) - 1
nstock = length(unique(data$TICKER)) - 1
days = unique(Ldata$date)[-1]
days = as.numeric(days)
#### There are in total 99 financial institutions. The starting date of sample is 2002-01-02, and the end date is 2016-12-31. The longest sample has 3776 return data.
### ARMA-Garch Filter -- individual stock ###
INDEX_id = ncol(R_com)
loss_stock = -R_com[,k]
na.pos = which(is.na(loss_stock))
if(length(na.pos) == 0){
na.pos = 10^6
}
spec = ugarchspec(mean.model=list(armaOrder=c(1,0),include.mean=T), distribution.model="sstd")
garch11.fit = ugarchfit(spec = spec, data =loss_stock[-na.pos], solver.control=list(trace = 1),solver = "hybrid")
##
## Iter: 1 fn: 7473.8248 Pars: -0.04453 -0.07113 0.06168 0.10975 0.88111 1.02609 5.40679
## Iter: 2 fn: 7473.8248 Pars: -0.04453 -0.07113 0.06168 0.10976 0.88111 1.02610 5.40640
## solnp--> Completed in 2 iterations
### Calcualte VaR for individual stock
nu=coef(garch11.fit)["shape"]; ga=coef(garch11.fit)["skew"]
qmodel = qdist("sstd",p=0.99,mu=0,sigma=1,shape=nu,skew=ga)
# VaR95 = VaR.forecasts(fit=garch11.fit,qmodel=qmodel,alpha=0.95,ret_z = FALSE)
VaR99 = VaR.forecasts(fit=garch11.fit,qmodel=qmodel,alpha=0.99,ret_z = FALSE)
##### ARMA-Garch Filter -- Index ###
garch11.fit.Index = ugarchfit(spec = spec, data = -R_com[-na.pos,INDEX_id], solver.control=list(trace = 1),solver = "hybrid")
##
## Iter: 1 fn: 5976.3353 Pars: -0.05197 -0.08092 0.01545 0.10291 0.89270 1.08213 7.20901
## Iter: 2 fn: 5976.3353 Pars: -0.05197 -0.08092 0.01545 0.10290 0.89270 1.08212 7.20866
## solnp--> Completed in 2 iterations
#### DCC Filter ##########
library(bayesDccGarch)
### Get de-mean residuals from ARMA filter above
resid_s = residuals(garch11.fit); mut_s = fitted(garch11.fit)
resid_I = residuals(garch11.fit.Index); mut_I = fitted(garch11.fit.Index)
### Fit DCC model to residuals
sim_y = cbind(array(resid_s), array(resid_I)) #Loss
bayesDcc <- bayesDccGarch(sim_y, nSim = 1000)
## Maximizing the log-posterior density function.
## Done.
## One approximation for covariance matrix of parameters cannot be directly computed through the hessian matrix.
## Calibrating the standard deviations for simulation:
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.34 0.09 0.23 0.07 0.12 0.08 0.21 0.05 0.12 0.05
## phi_11
## 0.03
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.35 0.15 0.25 0.12 0.16 0.15 0.21 0.11 0.20 0.05
## phi_11
## 0.02
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.39 0.24 0.23 0.17 0.19 0.23 0.22 0.17 0.16 0.10
## phi_11
## 0.06
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.37 0.24 0.26 0.19 0.16 0.22 0.20 0.17 0.15 0.13
## phi_11
## 0.10
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.36 0.24 0.25 0.19 0.19 0.23 0.20 0.19 0.29 0.18
## phi_11
## 0.13
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate:
## phi_1 phi_2 phi_3 phi_4 phi_5 phi_6 phi_7 phi_8 phi_9 phi_10
## 0.35 0.22 0.22 0.21 0.18 0.24 0.21 0.22 0.28 0.23
## phi_11
## 0.28
## Computing the covariance matrix of pilot sample.
## Done.
## Calibrating the Lambda coefficient:
## lambda: 0.5
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.05
## lambda: 0.4
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.05
## lambda: 0.32
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.11
## lambda: 0.256
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.11
## lambda: 0.2048
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.11
## lambda: 0.16384
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Accept Rate: 0.19
## Done.
## Starting the simulation by one-block random walk Metropolis-Hasting algorithm.
## Simulation number = 100
## Simulation number = 200
## Simulation number = 300
## Simulation number = 400
## Simulation number = 500
## Simulation number = 600
## Simulation number = 700
## Simulation number = 800
## Simulation number = 900
## Simulation number = 1000
## Done.
H = bayesDcc$H ### Sigma function
Z = matrix(NA, nrow(H), 2)
Sigma.t = array(NA,dim=c(2,2,nrow(H)) )
### Calculated standarized residuals and store in Z
for (i in 1:nrow(H)){
Sigma.t[,,i] = matrix(H[i,],2,2)
tmp1 = matrix(sim_y[i,],1,2)
tmp2 = Sigma.t[,,i]%^%(-1/2)
Z[i,] = tmp1 %*% tmp2
}
# Since we use Baysian method, we drop first 50 observations to avoid non-converenge
library(skewt)
begin = 50
end = nrow(Z)
#### Fit EVT method ######
fit = New_cov_update(Z[begin:end,], prob = 0.85)
print(c("alpha1","alpha2","rho","nu","sd1","sd2"))
## [1] "alpha1" "alpha2" "rho" "nu" "sd1" "sd2"
print(fit$EVT)
## [1] -0.04345775 0.12461728 0.16131336 4.15373426 0.74021534 0.85045754
library(gld)
print(fit.gpd(apply(Z[begin:500,],1,function(x) sqrt(x[1]^2 + x[2]^2))))
## $estA
## alpha beta delta lambda
## 0.6657047 0.9903079 0.8168256 0.1489984
##
## $estB
## alpha beta delta lambda
## 2.1696639 13.0308369 0.3080824 4.2219394
##
## $param
## [1] "gpd"
##
## attr(,"class")
## [1] "GldFitMultiple"
#### Extract parameters for multivariate skew-t from DCC model
nu = summary(bayesDcc)$statistics[1,1]
gamma1 = summary(bayesDcc)$statistics[2,1]
gamma2 = summary(bayesDcc)$statistics[6,1]
gamma = c(gamma1, gamma2)
#### Calculate CoVaR ##########
CoVaR = function(x, q, fit, VaR, Sigma, mu_s, mu_I){ # x is CoVaR
z = (c(VaR, x) - c(mu_s, mu_I)) %*% (Sigma%^%(-1/2))
diag_inv_EVT = matrix(c(1/fit$EVT[5], 0,0,1/fit$EVT[6]),2,2)
tmp_value_EVT = diag_inv_EVT %*% t(z)
z_EVT = tmp_value_EVT[2] / tmp_value_EVT[1] - fit$EVT[3]
theorm_EVT = 1-Thrm4_1(z_EVT, fit$EVT[1], fit$EVT[2], fit$EVT[3], fit$EVT[4])
# return(theorm_EVT-q)
return( (theorm_EVT - q)^2*10^10 )
}
CoVaR_Skew <- function (x, VaR, Sigma, gamma, nu, q, marginal, mu_s, mu_I) { # x is CoVaR
z = (c(VaR, x)- c(mu_s, mu_I)) %*% (Sigma%^%(-1/2))
joint = hcubature(f=dsst,lowerLimit = z, upperLimit = c(10^6,10^6) , tol = 1e-8, gamma = gamma, nu = nu)$integral
return( (joint - q*marginal)^2*10^10 )
}
CoVaR_func <- function(fit, VaR95, VaR99, Sigma.t,pos, gamma, nu, mut_s, mut_I,Z){
##### Marginal VaR #########
# VaR_Z95 = as.numeric((VaR95$evt[pos] - VaR95$mut[pos])/VaR95$sigt[pos])
VaR_Z99 = as.numeric((VaR99$evt[pos] - VaR99$mut[pos])/VaR99$sigt[pos])
# marginal95 = sum(Z[,1]>VaR_Z95)/length(Z[,1]) #1 - pskt(VaR_Z95, nu, gamma = gamma[1])
marginal99 = 1 - pskt(VaR_Z99, nu, gamma = gamma[1]) #sum(Z[,1]>VaR_Z99)/length(Z[,1]) #
tmp95_Skew =NA; tmp99_Skew = NA; tmp95 = NA; tmp99 = NA;
##### Parametric method ###########
# try(tmp95_Skew <- as.numeric(uniroot(CoVaR_Skew, q = 0.05, VaR = as.numeric(VaR95$evt[pos]), Sigma = Sigma.t[,,pos], gamma = gamma, nu = nu, marginal = marginal95, mu_s = as.numeric(mut_s[pos]), mu_I = as.numeric(mut_I[pos]), lower = 0.1, upper = 10, tol = 8e-6)$root, print =1))
try(tmp99_Skew<-as.numeric(nlminb(start=5, objective=CoVaR_Skew,q = 0.01, VaR = as.numeric(VaR99$evt[pos]), Sigma = Sigma.t[,,pos], gamma = gamma, nu = nu, marginal = marginal99, mu_s = as.numeric(mut_s[pos]), mu_I = as.numeric(mut_I[pos]))$par))
##### EVT method ###############
# try(tmp99 <- as.numeric(uniroot(CoVaR, q = 0.01, fit = fit, VaR = as.numeric(VaR99$evt[pos]), Sigma = Sigma.t[,,pos], mu_s = as.numeric(mut_s[pos]), mu_I = as.numeric(mut_I[pos]), lower = max(0.5,(tmp99_Skew-2)), upper = (tmp99_Skew+2), tol = 8e-6)$root, print =1))
try(tmp99 <- as.numeric(nlminb(start=tmp99_Skew, objective=CoVaR,q = 0.01, fit = fit, VaR = as.numeric(VaR99$evt[pos]), Sigma = Sigma.t[,,pos], mu_s = as.numeric(mut_s[pos]), mu_I = as.numeric(mut_I[pos]))$par))
return(c(marginal99, tmp99, tmp99_Skew))
}
pos = 50
all_pos = matrix(1:ndays, ndays,1)
all_pos = all_pos[-na.pos]
tmp = CoVaR_func(fit, VaR95, VaR99, Sigma.t, pos, gamma, nu, mut_s, mut_I, Z)
df = data.table(date = days[all_pos[pos]], firm = colnames(R_com)[k], marginal99 = tmp[1], CoVaR_evt99 = tmp[2], CoVaR_fp99 = tmp[3], VaR99 = as.numeric(VaR99$evt[pos]), Loss_stock = as.numeric(-R_com[all_pos[pos],k]), Loss_Index = as.numeric(-R_com[all_pos[pos], INDEX_id]))
print(df)
## date firm marginal99 CoVaR_evt99 CoVaR_fp99 VaR99 Loss_stock
## 1: 11761 10138 0.02008714 5.025907 5.469257 4.162123 -0.5865119
## Loss_Index
## 1: -1.767393