Cálculo manual y backtesting de VaR y CVaR de DCC-GARCH R

Question

Estimé un ajuste GARCH para los rendimientos logarítmicos de tres series (CAC 40, un índice inmobiliario francés y una serie de rendimientos de bonos franceses T10) utilizando rugarch . A continuación, calculé manualmente y realicé pruebas retrospectivas de las medidas VaR y CVaR. También ajusté un DCC-GARCH(1,1) a los rendimientos logarítmicos de las 3 series utilizando rmgarch y ahora me gustaría hacer un backtest de las medidas VaR y CVaR de forma similar a como lo hice para los casos GARCH univariantes.

Antes de continuar, tendremos que especificar las siguientes funciones para el CVaR:

#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)

Los datos se pueden encontrar en Dropbox en los siguientes enlaces (uno para los datos del índice inmobiliario francés y el otro para los bonos franceses) los datos del CAC 40 se pueden descargar con 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

El código comentado es el siguiente:

#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))

Ahora quiero hacer los pasos de backtesting y trazado para las medidas de VaR y CVaR. Lo ideal sería realizar también la prueba de Kupiec y Christoffersen como en la función report del paquete rugarch . Estoy realmente perplejo, ya que he intentado encontrar una respuesta en Internet pero no he podido.

Answer 1

Aparentemente, los argumentos del VaR son simplemente ignorados por la función. Si intentas ejecutar algo absurdo, como especificar que calculate.VaR = 'banana', la función sigue ejecutándose. El objeto tiene el mismo tamaño sin importar los argumentos para calculate.VaR y VaR.alpha, y no he podido encontrar ningún rastro de los valores. El plot(DCCroll, que = 5) le dará un gráfico EWMA del 2,5% de VaR, pero no he podido encontrar ninguna forma de recuperar los valores reales para el backtesting. Actualmente me encuentro con el mismo problema y, por lo que veo, este documento ( https://repositorio.iscte-iul.pt/bitstream/10071/17461/1/Master_Vladimir_Krecmer.pdf ) recuperó la covarianza a través de DCC, la utilizó para pronosticar los rendimientos reales (no estoy seguro de cómo), y luego hizo manualmente el VaR. Si, por casualidad, alguna vez tienes una actualización, dilo. Estoy considerando trabajar una función para hacerlo, pero no creo que tenga las habilidades.