1
$\begingroup$

I estimated a GARCH fit to the log returns of three series (CAC 40, a french real estate index and french T10 bond yield series) using rugarch. I then manually calculated and backtested the VaR and CVaR measures. I also fitted a DCC-GARCH(1,1) to the log returns of the 3 series using rmgarch and now I would like to backtest the VaR and CVaR measures in a similar way as I did for the univariate GARCH cases.

We'll need to specify the following functions for the CVaR before we proceed:

#This function calculates the CVaR at a certain position gdist list
cvar <- function(p=0.05, s = "CAC", dist_params = gdist_var, pos = l, v = df, dist = "jsu"){

  ES <- abs((integrate(qdist, lower = 0, upper = p, distribution = dist, mu = gdist_var[[s]][, 'Mu'][pos], sigma = gdist_var[[s]][, 'Sigma'][pos],
                       shape = gdist_var[[s]][, 'Shape'][pos], skew = gdist_var[[s]][, 'Skew'][pos])$value)/p * v[nrow(v),s])

  return(ES)
}

#This function calculates the CVaR given the arguments
cvar_df <- function(p=0.01, dist = "jsu", mu = Mu, sigma = Sigma, shape = Shape, skew = Skew){

  ES <- (integrate(qdist, lower = 0, upper = p, distribution = dist, mu = mu, sigma = sigma, shape = shape, skew = skew)$value)/p 

  return(ES)
}

#This function is a vectorized form of the above
vcvar_df <- Vectorize(cvar_df)

The data can be found on dropbox under the following links (one for french real estate index data and the other for french bonds) the CAC 40 data is downloadable with quantmod:

https://www.dropbox.com/s/vy8sl88fs5opmi3/IEIF%20SIIC%20FRANCE_quote_chart.csv?dl=0

https://www.dropbox.com/s/xljxk5izy6pt1ds/entre_obligations.csv?dl=0

The commented code is the following:

#loading libraries and data:

require(tidyquant)
require(reshape2)
require(astsa)
require(GGally)
require(forecast)
source("functions.R", local = T)

# #https://www.banque-france.fr/statistiques/taux-et-cours/les-indices-obligataires

obli_10 <- read.csv("entre_obligations.csv", sep = ";", na.strings = "-", stringsAsFactors = F) %>%
  rename(Date = 1) %>%
  mutate(Date = dmy(Date)) %>%
  mutate_at(vars(-Date), funs(gsub("\\,", ".", .))) %>%
  mutate_at(vars(-Date), funs(as.numeric)) %>%
  dplyr::select(c(1,2))

# #https://live.euronext.com/en/product/indices/QS0010980447-XPAR/quotes indices nu de:
# #https://www.ieif.fr/wp-content/plugins/aa-indices/datas/histo/index.php?IndiceNu=SIICNu&IndiceNet=SIICNet&IndiceBrut=SIICBrut&Indice=Euronext%20IEIF%20SIIC%20France
reit <- read.csv("IEIF SIIC FRANCE_quote_chart.csv", stringsAsFactors = F) %>%
  dplyr::select(1,2) %>%
  rename(Date = 1) %>%
  mutate(Date = substr(Date, 1, 10)) %>%
  mutate(Date = ymd(Date))

cac <- as.data.frame(Ad(getSymbols("^FCHI", src = "yahoo", adjust = T, auto.assign = FALSE)))

cac <- cac %>%
  mutate(Date = rownames(.)) %>%
  mutate(Date = ymd(Date)) %>%
  dplyr::select(Date, everything())

#Calculate the log returns

lr_df <- as.data.frame(sapply(df[2:ncol(df)], function(x) diff(log(x))))

lr_df <-cbind(df$Date[2:nrow(df)], lr_df) %>%
  dplyr::rename(Date = !!names(.)[1])

#Specification of GARCH models

cac_egarch_spec <- ugarchspec(mean.model = list(armaOrder = c(3, 3), include.mean = T, archm = F, archpow = 1),
                              variance.model = list(model = "eGARCH", garchOrder = c(2, 1)),
                              distribution.model="jsu")

reit_egarch_spec <- ugarchspec(mean.model = list(armaOrder = c(3, 1), include.mean = T, archm = F, archpow = 1),
                               variance.model = list(model = "eGARCH", garchOrder = c(2, 1)),
                               distribution.model="nig")

obli_apgarch_spec <- ugarchspec(mean.model = list(armaOrder = c(2, 1), include.mean = T, archm = F, archpow = 1),
                                 variance.model = list(model = "apARCH", garchOrder = c(1, 1)),
                                 distribution.model="jsu")

#Get VaR and CVaR

cac_roll <- ugarchroll(cac_egarch_spec, lr_df[,2],n.start = 750, refit.every = 50, refit.window = "moving",
                                   solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = c(0.01, 0.025, 0.05), keep.coef = T,
                                   fit.control = list(scale = 1))

reit_roll <- ugarchroll(reit_egarch_spec, lr_df[,3],n.start = 750, refit.every = 50, refit.window = "moving",
                       solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = c(0.01, 0.025, 0.05), keep.coef = T,
                       fit.control = list(scale = 1))

obli_roll <- ugarchroll(obli_apgarch_spec, lr_df[,4],n.start = 750, refit.every = 50, refit.window = "moving",
                        solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = c(0.01, 0.025, 0.05), keep.coef = T,
                        fit.control = list(scale = 1))

gdist_var <- list()

gdist_var[["CAC"]] <- as.data.frame(cac_roll, which = 'density') 
gdist_var[["REIT"]] <- as.data.frame(reit_roll, which = 'density') 
gdist_var[["OBLI_10"]] <- as.data.frame(obli_roll, which = 'density')

#VaR and CVaR calculations
p <- c(0.05, 0.025, 0.01)
l <- nrow(gdist_var[["CAC"]])

for(j in p){
  for(i in 1:3){
    print(paste("VaR", names(gdist_var)[i], 1-j))
    print(abs(qdist(dg[[i]], p=j, mu=gdist_var[[i]]$Mu[l], sigma=gdist_var[[i]]$Sigma[l], skew=gdist_var[[i]]$Skew[l], shape=gdist_var[[i]]$Shape[l]))*df[nrow(df),i+1])
  }
}

for(j in p){
  for(i in 1:3){
    print(paste("CVaR", names(gdist_var)[i], 1-j))
    print(cvar(p = j, s = names(gdist_var)[i], dist_params = gdist_var, pos = l, v = df, dist = dg[[i]]))
  }
}

#VaR plots for cac only but will be done for others the same way
var_cac <- gdist_var$CAC
var_cac <- cbind.data.frame(tail(lr_df[,c("Date","CAC")],2438), var_cac) %>%
  dplyr::select(-`Shape(GIG)`, -Realized) %>%
  dplyr::mutate(VaR_99 = qdist("jsu", p = 0.01, mu = Mu, sigma = Sigma, skew = Skew, shape = Shape)) %>%
  dplyr::select(-Mu, -Sigma, -Skew, -Shape)
var_cac <- melt(var_cac, id.vars = "Date")
ggplot(data = var_cac, aes(x = Date, value)) + geom_line(aes(colour = variable)) +
  ggtitle("Series with 1% 1D VaR Limit") +
  theme(plot.title = element_text(hjust = 0.5))

#VaR backtesting reports using report function
report(cac_roll, type = "VaR", VaR.alpha = 0.05, conf.level = 0.95)
report(reit_roll, type = "VaR", VaR.alpha = 0.05, conf.level = 0.95)
report(obli_roll, type = "VaR", VaR.alpha = 0.05, conf.level = 0.95)

#CVaR plots for CAC only but will be done for others
cvar_cac <- gdist_var$CAC
cvar_cac <- cbind.data.frame(tail(lr_df[,c("Date","CAC")],2438), cvar_cac) %>%
  dplyr::select(-`Shape(GIG)`, -Realized) %>%
  dplyr::mutate(CVaR_99 = vcvar_df(p = 0.01, dist = "jsu", mu = Mu, sigma = Sigma, shape = Shape, skew = Skew)) %>%
  dplyr::select(-Mu, -Sigma, -Skew, -Shape)

mcvar_cac <- melt(cvar_cac, id.vars = "Date")
ggplot(data = mcvar_cac, aes(x = Date, value)) + geom_line(aes(colour = variable)) +
  ggtitle("Series with 1% 1D CVaR Limit") +
  theme(plot.title = element_text(hjust = 0.5))

#Bactesting CVaRby calculating nuber of times CVaR crossed
cvar_cac <- gdist_var$CAC
cvar_cac <- cbind.data.frame(tail(lr_df[,c("Date","CAC")],2438), cvar_cac) %>%
  dplyr::select(-`Shape(GIG)`, -Realized) %>%
  dplyr::mutate(CVaR_99 = vcvar_df(p = 0.01, dist = "jsu", mu = Mu, sigma = Sigma, shape = Shape, skew = Skew)) %>%
  dplyr::mutate(CVaR_975 = vcvar_df(p = 0.025, dist = "jsu", mu = Mu, sigma = Sigma, shape = Shape, skew = Skew)) %>%
  dplyr::mutate(CVaR_95 = vcvar_df(p = 0.05, dist = "jsu", mu = Mu, sigma = Sigma, shape = Shape, skew = Skew)) %>%
  mutate(depasse_99 = case_when(CVaR_99 >= .[[2]] ~ 1, TRUE ~ 0)) %>%
  mutate(depasse_975 = case_when(CVaR_975 >= .[[2]] ~ 1, TRUE ~ 0)) %>%
  mutate(depasse_95 = case_when(CVaR_95 >= .[[2]] ~ 1, TRUE ~ 0)) %>%
  mutate(sum_99 = sum(depasse_99)) %>%
  mutate(sum_975 = sum(depasse_975)) %>%
  mutate(sum_95 = sum(depasse_95))

#DCC GARCH of GARCH models above:
require(rmgarch)

dcc_garch <- multispec(c(cac_egarch_spec, reit_egarch_spec, obli_apgarch_spec))
dcc_multfit <- multifit(dcc_garch, lr_df[,2:ncol(lr_df)]) #fitting many univariate models
dcc_spec <- dccspec(uspec = dcc_garch, dccOrder = c(1,1), distribution = "mvnorm")
dcc_fit <- dccfit(dcc_spec, lr_df[,2:ncol(lr_df)], fit.control = list(eval.se = TRUE), fit = dcc_multfit) #fit = dcc_multfit not really necessary but more robust
dcc_roll <- dccroll(dcc_spec, lr_df[,2:4],n.start = 750, refit.every = 50, refit.window = "moving",
                       solver = "solnp", calculate.VaR = TRUE, VaR.alpha = c(0.01, 0.025, 0.05), keep.coef = T,
                       fit.control = list(scale = 1))

Now I want to do the backtesting and plotting steps for both VaR and CVaR measures. Ideally I would also conduct the Kupiec and Christoffersen test just like in the function report of the package rugarch. I am realy stumped as I tried to find an answer online but couldn't.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
0
$\begingroup$

Apparently, the VaR arguments are simply ignored by the function. If you try running something absurd, such as specifying that calculate.VaR = 'banana', the function still runs. The object has the same size no matter the argumentos for calculate.VaR and VaR.alpha, and I was unable to find any trace of the values. The plot(DCCroll, which = 5) will give you a EWMA plot of 2.5% VaR, but I've been unable to find any way of recovering the actual values for backtesting. I'm currently facing the same problem, and, as far as I can tell, this paper (https://repositorio.iscte-iul.pt/bitstream/10071/17461/1/Master_Vladimir_Krecmer.pdf) recovered the covariance through DCC, used it to forecast the actual returns (I'm not sure how), then manually did the VaR. If, by any chance, you ever get an update, do tell. I'm considering working a function for doing it, but I don't think I have the skills.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.