-
Notifications
You must be signed in to change notification settings - Fork 11
Downside Risk Evaluation with the R Package GAS
Leopoldo Catania edited this page Sep 12, 2018
·
1 revision
# This code reproduces the results of:
# Ardia, Boudt, and Catania (2018) 'Downside Risk Evaluation with the R Package GAS'.
# The first part of the script reproduces the chuck codes reported in the paper.
# The second part of the code reproduces the results obtained in the empirical section.
# The code made use of parallel computation.
##############################################################
# INSTALLATION AND SESSION INFO #
##############################################################
# install.packages("GAS") # version 0.2.8
# R version 3.5.0 (2018-04-23)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 7 x64 (build 7601) Service Pack 1
#
# Matrix products: default
#
# locale:
# [1] LC_COLLATE=French_Switzerland.1252 LC_CTYPE=French_Switzerland.1252 LC_MONETARY=French_Switzerland.1252
# [4] LC_NUMERIC=C LC_TIME=French_Switzerland.1252
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] GAS_0.2.8
#
# loaded via a namespace (and not attached):
# [1] cubature_1.4 MASS_7.3-50 zoo_1.8-3 compiler_3.5.0 parallel_3.5.0 tools_3.5.0 xts_0.11-0
# [8] yaml_2.2.0 Rcpp_0.12.18 grid_3.5.0 truncnorm_1.0-8 numDeriv_2016.8-1 Rsolnp_1.16 lattice_0.20-35
##############################################################
# FIRST PART #
##############################################################
rm(list = ls())
options(digits = 7, max.print = 40)
library("GAS")
data("dji30ret", package = "GAS")
dji30ret <- tail(dji30ret, 2500)
GASSpec_N <- UniGASSpec(Dist = "norm", GASPar = list(scale = TRUE))
GASSpec_ST <- UniGASSpec(Dist = "std", GASPar = list(scale = TRUE))
GASSpec_SKST <- UniGASSpec(Dist = "sstd", GASPar = list(scale = TRUE))
library("parallel")
cluster <- makeCluster(2)
H <- 1000
Roll_N <- UniGASRoll(dji30ret[, "GE"], GASSpec_N, RefitEvery = 5,
cluster = cluster, ForecastLength = H)
alpha <- 0.01
VaR_N <- quantile(Roll_N, probs = alpha)
ES_N <- ES(Roll_N, probs = alpha)
VaRBacktest_N <- BacktestVaR(data = tail(dji30ret[, "GE"], H), VaR = VaR_N, alpha = alpha)
VaRBacktest_N$DQ
Roll_ST <- UniGASRoll(dji30ret[, "GE"], GASSpec_ST, RefitEvery = 5,
cluster = cluster, ForecastLength = H)
VaR_ST <- quantile(Roll_ST, probs = alpha)
VaRBacktest_ST <- BacktestVaR(data = tail(dji30ret[, "GE"], H), VaR = VaR_ST, alpha = alpha)
VaRBacktest_ST$DQ
round(VaRBacktest_ST$Loss$Loss / VaRBacktest_N$Loss$Loss, 2)
FZL <- FZLoss(data = tail(dji30ret[, "GE"], H), VaR = VaR_N, ES = ES_N, alpha = alpha)
stopCluster(cluster)
##############################################################
# SECOND PART #
##############################################################
rm(list = ls())
options(digits = 7, max.print = 40)
library("GAS")
library("parallel")
Ncluster <- 36
################## FOLDERS ##################
## Path where all the results are saved
sDir <- paste0(getwd(), "/")
sDir_Tables <- paste(sDir, "Tables/", sep = "")
dir.create(sDir_Tables, showWarnings = FALSE)
sDir_Empirical <- paste(sDir, "Empirical/", sep = "")
dir.create(sDir_Empirical, showWarnings = FALSE)
sDir_Figures <- paste(sDir, "Figures/", sep = "")
dir.create(sDir_Figures, showWarnings = FALSE)
################## DATA ##################
data("dji30ret", package = "GAS")
dji30ret <- tail(dji30ret, 2500)
vAsset <- colnames(dji30ret)
########## MODEL SPECIFICATION ###############
uGASSpec_sstd <- UniGASSpec(Dist = "sstd", ScalingType = "Identity",
GASPar = list(location = FALSE, scale = TRUE, shape = FALSE))
uGASSpec_std <- UniGASSpec(Dist = "std", ScalingType = "Identity",
GASPar = list(location = FALSE, scale = TRUE, shape = FALSE))
uGASSpec_norm <- UniGASSpec(Dist = "norm", ScalingType = "Identity",
GASPar = list(location = FALSE, scale = TRUE, shape = FALSE))
#### Forecast design
iS <- 1500 # length of the in sample period
iH <- nrow(dji30ret) - iS
RefitEvery <- 21 # length of refit window
RefitWindow <- "moving" # updating scheme
iN <- ncol(dji30ret)
#### clusters
if (!is.null(Ncluster)) {
cluster <- makeCluster(Ncluster)
} else {
cluster <- NULL
}
### Univariate Forecast # Note: This is computationally expensive
luGASRoll_sstd <- list()
luGASRoll_std <- list()
luGASRoll_norm <- list()
for (i in 1:iN) {
cat(i, "\n")
luGASRoll_sstd[[i]] <- UniGASRoll(data = dji30ret[, i], GASSpec = uGASSpec_sstd, ForecastLength = iH,
RefitEvery = RefitEvery, RefitWindow = RefitWindow, cluster = cluster)
luGASRoll_std[[i]] <- UniGASRoll(data = dji30ret[, i], GASSpec = uGASSpec_std, ForecastLength = iH,
RefitEvery = RefitEvery, RefitWindow = RefitWindow, cluster = cluster)
luGASRoll_norm[[i]] <- UniGASRoll(data = dji30ret[, i], GASSpec = uGASSpec_norm, ForecastLength = iH,
RefitEvery = RefitEvery, RefitWindow = RefitWindow, cluster = cluster)
}
names(luGASRoll_sstd) <- vAsset
names(luGASRoll_std) <- vAsset
names(luGASRoll_norm) <- vAsset
save(luGASRoll_sstd, file = paste(sDir_Empirical, "/luGASRoll_sstd.RData", sep = ""))
save(luGASRoll_std, file = paste(sDir_Empirical, "/luGASRoll_std.RData", sep = ""))
save(luGASRoll_norm, file = paste(sDir_Empirical, "/luGASRoll_norm.RData", sep = ""))
################## EMPIRICAL ##################
### VaR and ES Analysis
load(file = paste(sDir_Empirical, "/luGASRoll_sstd.RData", sep = ""))
load(file = paste(sDir_Empirical, "/luGASRoll_std.RData", sep = ""))
load(file = paste(sDir_Empirical, "/luGASRoll_norm.RData", sep = ""))
vAlpha <- c(0.05, 0.01)
lVaR <- list()
lES <- list()
vDist <- c("norm", "std", "sstd")
clusterEvalQ(cluster, library(GAS))
clusterExport(cluster, c("vDist", "vAlpha", "iH", paste("luGASRoll_", vDist, sep = "")))
lVaR <- parLapply(cluster, vAsset, function(sAsset) {
mOut <- array(NA, dim = c(iH, length(vDist), length(vAlpha)),
dimnames = list(1:iH, vDist, vAlpha))
for (sDist in vDist) {
obj_GAS <- get(paste("luGASRoll_", sDist, sep = ""))
for (dAlpha in vAlpha) {
mOut[, sDist, paste(dAlpha)] <- quantile(obj_GAS[[sAsset]], dAlpha)
}
}
return(mOut)
})
names(lVaR) <- vAsset
lES <- parLapply(cluster, vAsset, function(sAsset) {
mOut <- array(NA, dim = c(iH, length(vDist), length(vAlpha)),
dimnames = list(1:iH, vDist, vAlpha))
for (sDist in vDist) {
obj_GAS <- get(paste("luGASRoll_", sDist, sep = ""))
for (dAlpha in vAlpha) {
mOut[, sDist, paste(dAlpha)] <- ES(obj_GAS[[sAsset]], dAlpha)
}
}
return(mOut)
})
names(lES) <- vAsset
if (!is.null(Ncluster)) stopCluster(cluster)
save(lVaR, file = paste(sDir_Empirical, "/lVaR.RData", sep = ""))
save(lES, file = paste(sDir_Empirical, "/lES.RData", sep = ""))
load(file = paste(sDir_Empirical, "/lVaR.RData", sep = ""))
load(file = paste(sDir_Empirical, "/lES.RData", sep = ""))
## Loss VaR
lVaRBacktest <- list()
for (sAsset in vAsset) {
lVaRBacktest[[sAsset]] <- list()
for (sDist in vDist) {
lVaRBacktest[[sAsset]][[sDist]] = list()
obj_GAS <- get(paste("luGASRoll_", sDist, sep = ""))
vRealised <- tail(getObs(obj_GAS[[sAsset]]), iH)
for (dAlpha in vAlpha) {
lVaRBacktest[[sAsset]][[sDist]][[paste(dAlpha)]] <- BacktestVaR(data = vRealised,
VaR = lVaR[[sAsset]][, sDist, paste(dAlpha)], alpha = dAlpha)
}
}
}
save(lVaRBacktest, file = paste(sDir_Empirical, "/lVaRBacktest.RData", sep = ""))
## Loss VaR and ES
lVaRESBacktest <- list()
for (sAsset in vAsset) {
lVaRESBacktest[[sAsset]] <- list()
for (sDist in vDist) {
lVaRESBacktest[[sAsset]][[sDist]] = list()
obj_GAS <- get(paste("luGASRoll_", sDist, sep = ""))
vRealised <- tail(getObs(obj_GAS[[sAsset]]), iH)
for (dAlpha in vAlpha) {
lVaRESBacktest[[sAsset]][[sDist]][[paste(dAlpha)]] <- FZLoss(data = vRealised,
VaR = lVaR[[sAsset]][, sDist, paste(dAlpha)],
ES = lES[[sAsset]][, sDist, paste(dAlpha)],
alpha = dAlpha)
}
}
}
save(lVaRESBacktest, file = paste(sDir_Empirical, "/lVaRESBacktest.RData", sep = ""))
# Extract Losses
cLoss <- array(NA, dim = c(iH, length(vAsset), length(vDist), length(vAlpha), 2),
dimnames = list(1:iH, vAsset, vDist, vAlpha, c("QL", "FZL")))
for (sAsset in vAsset) {
for (sDist in vDist) {
for (dAlpha in vAlpha) {
cLoss[, sAsset, sDist, paste(dAlpha), "QL"] <- lVaRBacktest[[sAsset]][[sDist]][[paste(dAlpha)]]$Loss$LossSeries
cLoss[, sAsset, sDist, paste(dAlpha), "FZL"] <- lVaRESBacktest[[sAsset]][[sDist]][[paste(dAlpha)]]
}
}
}
# #
save(cLoss, file = paste(sDir_Empirical, "/cLoss.RData", sep = ""))
# #
################ TABLES ########################
# Testing difference in predictive ability between GAS specifications
### clear workspace
rm(list = ls()[-which(ls() %in% c("sDir", "sDir_Empirical", "sDir_Figures", "sDir_Tables"))])
####### VaR analysis
load(paste(sDir_Empirical, "/cLoss.RData", sep = ""))
vAlpha <- dimnames(cLoss)[[4]]
vAsset <- dimnames(cLoss)[[2]]
vDist <- dimnames(cLoss)[[3]]
## Losses VaR
vModels <- paste("GAS", vDist, sep = "-")
load(paste(sDir_Empirical, "/lVaRBacktest.RData", sep = ""))
load(paste(sDir_Empirical, "/lVaRESBacktest.RData", sep = ""))
vTestNames <- c("LRuc", "LRcc", "DQ", "Loss-QL", "Loss-FZL", "AE", "ADmean", "ADmax")
cBackTest <- array(NA, dim = c(length(vAsset), length(vTestNames),length(vModels), length(vAlpha)),
dimnames = list(vAsset, vTestNames, vModels, vAlpha))
for (sAsset in vAsset) {
for (sDist in vDist) {
for (dAlpha in vAlpha) {
cBackTest[sAsset, "LRuc", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$LRuc["Pvalue"]
cBackTest[sAsset, "LRcc", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$LRcc["Pvalue"]
cBackTest[sAsset, "DQ", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$DQ$pvalue
cBackTest[sAsset, "AE", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$AE
cBackTest[sAsset, "ADmean", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$AD["ADmean"]
cBackTest[sAsset, "ADmax", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$AD["ADmax"]
cBackTest[sAsset, "Loss-QL", paste("GAS", sDist, sep = "-"), dAlpha] <- lVaRBacktest[[sAsset]][[sDist]][[dAlpha]]$Loss$Loss
cBackTest[sAsset, "Loss-FZL", paste("GAS", sDist, sep = "-"), dAlpha] <- mean(lVaRESBacktest[[sAsset]][[sDist]][[dAlpha]])
}
}
}
mTab_DQ <- format(round(cbind(cBackTest[, "DQ", , "0.01"], cBackTest[, "DQ", , "0.05"]), 2), digits = 2, scientific = FALSE)
mTab_DQ_numeric <- cbind(cBackTest[, "DQ", , "0.01"], cBackTest[, "DQ", , "0.05"])
## Define shadow cells using the \usepackage[table]{xcolor} package
## In TeX define
# \newcommand{\grb}[1]{\cellcolor[gray]{0.6}{#1}}
mTab_DQ[mTab_DQ_numeric < 0.01] <- paste("\\grb{",mTab_DQ[mTab_DQ_numeric < 0.01], "}", sep = "")
mTab_DQ[, ncol(mTab_DQ)] <- paste(mTab_DQ[, ncol(mTab_DQ)], "\\\\")
mTab_DQ <- mTab_DQ[sort(rownames(mTab_DQ)), ]
write.table(mTab_DQ, file = paste(sDir_Tables, "/DQ.txt", sep = ""), sep = " & ", quote = FALSE,
row.names = TRUE, col.names = FALSE)
# QL Loss
mTab_Loss_numeric <- cbind(cBackTest[, "Loss-QL", , "0.01"]/cBackTest[, "Loss-QL", "GAS-norm", "0.01"],
cBackTest[, "Loss-QL", , "0.05"]/cBackTest[, "Loss-QL", "GAS-norm", "0.05"])
mTab_Loss <- format(round(cbind(cBackTest[, "Loss-QL", , "0.01"]/cBackTest[, "Loss-QL", "GAS-norm", "0.01"],
cBackTest[, "Loss-QL", , "0.05"]/cBackTest[, "Loss-QL", "GAS-norm", "0.05"]), 2),
digits = 2, scientific = FALSE)
mTab_Loss[mTab_Loss_numeric > 1.0] <- paste("\\grb{",mTab_Loss[mTab_Loss_numeric > 1.0], "}", sep = "")
mTab_Loss_QL = mTab_Loss[, -c(1, 4)]
## Loss VaR-ES
mTab_Loss_numeric <- cbind(cBackTest[, "Loss-FZL", , "0.01"]/cBackTest[, "Loss-FZL", "GAS-norm", "0.01"],
cBackTest[, "Loss-FZL", , "0.05"]/cBackTest[, "Loss-FZL", "GAS-norm", "0.05"])
mTab_Loss <- format(round(cbind(cBackTest[, "Loss-FZL", , "0.01"]/cBackTest[, "Loss-FZL", "GAS-norm", "0.01"],
cBackTest[, "Loss-FZL", , "0.05"]/cBackTest[, "Loss-FZL", "GAS-norm", "0.05"]), 2),
digits = 2, scientific = FALSE)
mTab_Loss[mTab_Loss_numeric > 1.0] <- paste("\\grb{",mTab_Loss[mTab_Loss_numeric > 1.0], "}", sep = "")
mTab_Loss_FZL = mTab_Loss[, -c(1, 4)]
mTab_Loss = cbind(mTab_Loss_QL, mTab_Loss_FZL)
mTab_Loss[, ncol(mTab_Loss)] <- paste(mTab_Loss[, ncol(mTab_Loss)], "\\\\")
mTab_Loss <- mTab_Loss[sort(rownames(mTab_Loss)), ]
write.table(mTab_Loss, file = paste(sDir_Tables, "/mTab_Loss.txt", sep = ""), sep = " & ", quote = FALSE,
row.names = TRUE, col.names = FALSE)
################ FIGURE ###########################
### clear workspace
rm(list = ls()[-which(ls() %in% c("sDir", "sDir_Empirical", "sDir_Figures", "sDir_Tables"))])
load(file = paste(sDir_Empirical, "/lVaR.RData", sep = ""))
vAsset <- sort(names(lVaR))
vAlpha <- dimnames(lVaR[[1]])[[3]]
vDist <- dimnames(lVaR[[1]])[[2]]
iH <- nrow(lVaR[[1]])
data("dji30ret", package = "GAS")
dji30ret_oos = tail(dji30ret, iH)
sAsset <- vAsset[11]
sAlpha <- "0.01"
vY <- dji30ret_oos[, sAsset]
cVaR <- lVaR[[sAsset]][, c("norm", "std"), sAlpha, drop = "FALSE"]
vDates <- as.Date(rownames(dji30ret_oos))
pdf(paste(sDir_Figures, "/VaR_", sAsset, ".pdf", sep = ""), width = 15, height = 10)
layout(matrix(1, 1, 1), heights = c(2, 2, 2, 2))
vLim <- c(min(cVaR, vY), max(cVaR, vY))
plot(vDates, rep(0, length(vDates)), type = "n", xlab = "", ylab = "",
xaxt = "n", cex.axis = 1.5, las = 1, ylim = vLim)
grid(nx = 10, ny = 10, col = "gray", lty = "dotted")
points(vDates, vY, col = "gray30", cex = 1.5)
lines(vDates, cVaR[, "norm", "0.01"], col = "red", lwd = 2)
lines(vDates, cVaR[, "std", "0.01"], col = "blue", lwd = 2, lty = 2)
legend("topleft", legend = c("GAS-N", "GAS-ST"), lwd = 2, lty = c(1, 2),
col = c("red", "blue"), cex = 1.5, bg = "white")
Label_Dates <- seq(min(vDates), max(vDates) + 10, "quarter")
Label_Dates_Q <- sapply(Label_Dates, function(x) substr(x, 1, 7))
axis.Date(1, at = Label_Dates, labels = Label_Dates_Q, cex.axis = 1.5)
axis.Date(1, at = seq(min(vDates), max(vDates) + 0, "month"), labels = FALSE, tcl = -0.2)
dev.off()