[R-SIG-Finance] Manually calculating and backtesting VaR and CVaR from DCC-GARCH

[Eliot Tabet]

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 1Day VaR and
1Day 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.

-- 
Eliot TABET | Structurer/Quantitative Analyst



Les informations contenues dans ce message sont confidentielles.
Si vous receviez ce message par erreur, merci de prévenir l'expéditeur
et de le supprimer ainsi que ses pièces jointes.

N'imprimez ce mail que si c'est nécessaire

	[[alternative HTML version deleted]]