scripts/sizeLoss.R

In 3schwartz/SpecialeScrAndFun: Script and functions used for master thesis.

library(sfsmisc)
library(qrmtools)
library(Data)
library(scrAndFun)

outputFig = "C:/Users/sch/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/"
#### VaR and ES not always bigger for t-distribution than normal ####
V <- 1000
sig <- 0.2 / sqrt(365.25)
nu <- 4

alpha <- 1 - 2^seq(log(1-0.001, base = 2), -10, length.out = 256)
scale = sig * V

VaR.norm <- scale * qnorm(alpha)
VaR.t <- scale * sqrt((nu - 2) / nu) * qt(alpha, df = nu)
ES.norm <- (scale/(1 - alpha)) * dnorm(qnorm(alpha))
ES.t <- (scale*sqrt((nu-2)/nu)/(1 - alpha)) * dt(qt(alpha, df = nu), df = nu) * (nu + qt(alpha, df = nu)^2)/(nu - 1)

## Plot
pdf(file = paste0(outputFig, "VaRESexample.pdf"),
    width = 8, height = 6)
plot(1-alpha, ES.t, type = "l", ylim = range(VaR.norm, VaR.t, ES.norm, ES.t), log = "x",
     col = "maroon3", xlab = expression(1-alpha), ylab = "")
lines(1-alpha, VaR.t, col = "darkorange2")
lines(1-alpha, ES.norm,  col = "royalblue3")
lines(1-alpha, VaR.norm, col = "black")
legend("topright", bty = "n", lty = rep(1,4), col = c("maroon3", "darkorange2",
                                                      "royalblue3", "black"),
       legend = c(substitute(ES[alpha]~~"for "*italic(t[nu.])*" model", list(nu. = nu)),
                  substitute(VaR[alpha]~~"for "*italic(t[nu.])*" model", list(nu. = nu)),
                  expression(ES[alpha]~~"for normal model"),
                  expression(VaR[alpha]~~"for normal model")))
dev.off()

coinList.xts <- list(
  BitCoinYahoo.xts = xts::xts(BitCoinYahoo[,2:7], order.by = BitCoinYahoo[,1]),
  BitCashYahoo.xts = xts::xts(BitCashYahoo[,2:7], order.by = BitCashYahoo[,1]),
  EthereumYahoo.xts = xts::xts(EthereumYahoo[,2:7], order.by = EthereumYahoo[,1]),
  RippleYahoo.xts = xts::xts(RippleYahoo[,2:7], order.by = RippleYahoo[,1])
)

MSCI.xts <- list(
  days = xts::xts(MSCI_DP[,2], order.by = MSCI_DP[,1]),
  weeks = xts::xts(MSCI_WP[,2], order.by = MSCI_WP[,1]),
  months = xts::xts(MSCI_MP[,2], order.by = MSCI_MP[,1])
)

dailyExcess <- Log_Period_Vol(listCoins = coinList.xts, listIndex = MSCI.xts,
                              period = "days",
                              delta = 60/61, annu = 365.25,
                              dateRange = list(
                                BitCoin = "2013/",
                                BitCash = "2017-08-01/",
                                Ethereum = "2016-02-01/",
                                Ripple = "2014-01-01/",
                                Cum = "2013/"
                              ))

BTC_RFC <-  dailyExcess$BitCoin$excessReturn
BCH_RFC <-  dailyExcess$BitCash$excessReturn
XRP_RFC <-  dailyExcess$Ripple$excessReturn
ETH_RFC <-  dailyExcess$Ethereum$excessReturn

loss_Es_VaR(BTC_RFC, name = "Bitcoin", outputFig = outputFig)
loss_Es_VaR(BCH_RFC, name = "Bitcash", outputFig = outputFig)
loss_Es_VaR(XRP_RFC, name = "Ripple", outputFig = outputFig)
loss_Es_VaR(ETH_RFC, name = "Ethereum", outputFig = outputFig)

n <- 2500
B <- 1000

th <- 2
set.seed(271)
L <- rPar(n, shape = th)
name = "Pareto"
alpha <- 0.99
(VaR. <- VaR_np(L, level = alpha))
(ES.  <-  ES_np(L, level = alpha, verbose = TRUE))

pdf(file = paste0(outputFig, "VaRESalpha_",name,".pdf"),
    width = 8, height = 6)
plot(L, ylab = c("Loss"),
     main = bquote(atop(~widehat(VaR)[alpha]~"and"~widehat(ES)[alpha]~ " for "~ alpha~"="~.(alpha),
                        .(name))))
abline(h = c(VaR., ES.), lty = 2, col = c("royalblue3", "maroon3"))
legend("topleft", bty = "n", lwd = c(2,2),
       lty = c("dashed", "dashed"),
       col = c("royalblue3", "maroon3"),
       legend = c(
         expression(widehat(VaR)[alpha]),
         expression(widehat(ES)[alpha])),
       cex = 1)
dev.off()

alpha <- 1-10^-seq(0.5, 5, by = 0.05)
stopifnot(0 < alpha, alpha < 1)
VaR. <- VaR_np(L, level = alpha)
ES.  <-  ES_np(L, level = alpha, verbose = TRUE)
if(FALSE)
  warnings()

VaR.Par. <- VaR_Par(alpha, shape = th)
ES.Par.  <-  ES_Par(alpha, shape = th)

pdf(file = paste0(outputFig, "VaRESofAlpha_",name,".pdf"),
    width = 8, height = 6)
plot(1-alpha, ES.Par., type = "l", ylim = ran, log = "xy",
     xlab = expression(1-alpha),
     ylab = expression("True"~VaR[alpha]*","~ES[alpha]~"and estimates"),
     main = "Pareto simulated losses")
lines(1-alpha, ES., type = "l", col = "maroon3")
lines(1-alpha, VaR.Par., type = "l")
lines(1-alpha, VaR., type = "l", col = "royalblue3")
legend("topright", bty = "n", y.intersp = 1.2, lty = rep(1, 3),
       col = c("black", "maroon3", "royalblue3"),
       legend = c(expression(ES[alpha]~"and"~VaR[alpha]),
                  expression(widehat(ES)[alpha]),
                  expression(widehat(VaR)[alpha])))
dev.off()

set.seed(271)
VaR.boot <- bootstrap(L, B = B, level = alpha)
stopifnot(all(!is.na(VaR.boot)))

VaR.boot.    <- rowMeans(VaR.boot)
VaR.boot.var <- apply(VaR.boot, 1, var)
VaR.boot.CI  <- apply(VaR.boot, 1, CI)

set.seed(271)
system.time(ES.boot <- bootstrap(L, B = B, level = alpha, method = "ES"))
if(FALSE)
  warnings()

isNaN <- is.nan(ES.boot)
percNaN <- 100 * apply(isNaN, 1, mean)

pdf(file = paste0(outputFig, "numberOfNaN_ES_",name,".pdf"),
    width = 8, height = 6)
plot(1-alpha, percNaN, type = "l", log = "x",
     main = "Pareto",
     xlab = expression(1-alpha), ylab = expression("% of NaN when estimating "~ES[alpha]))
dev.off()

na.rm <- TRUE
ES.boot.    <- rowMeans(ES.boot, na.rm = na.rm)
ES.boot.var <- apply(ES.boot, 1, var, na.rm = na.rm)
ES.boot.CI  <- apply(ES.boot, 1, CI, na.rm = na.rm)
ran <- range(VaR.Par.,
             VaR.,
             VaR.boot.,
             VaR.boot.var,
             VaR.boot.CI,
             ES.Par., ES., ES.boot., ES.boot.var, ES.boot.CI, na.rm = TRUE)
pdf(file = paste0(outputFig, "VaRESBootAlpha_",name,".pdf"),
    width = 8, height = 6)
plot(1-alpha, VaR.Par., type = "l", log = "xy", xaxt = "n", yaxt = "n",
     ylim = ran, xlab = expression(1-alpha), ylab = "",
     main = "Pareto")
lines(1-alpha, VaR., lty = "dashed", lwd = 2, col = "royalblue3")
lines(1-alpha, VaR.boot., lty = "solid", col = "royalblue3")
lines(1-alpha, VaR.boot.var, lty = "dotdash", lwd = 1.4, col = "royalblue3")
lines(1-alpha, VaR.boot.CI[1,], lty = "dotted", col = "royalblue3")
lines(1-alpha, VaR.boot.CI[2,], lty = "dotted", col = "royalblue3")
lines(1-alpha, ES.Par.)
lines(1-alpha, ES., lty = "dashed", lwd = 2, col = "maroon3")
lines(1-alpha, ES.boot., lty = "solid", col = "maroon3")
lines(1-alpha, ES.boot.var, lty = "dotdash", lwd = 1.4, col = "maroon3")
lines(1-alpha, ES.boot.CI[1,], lty = "dotted", col = "maroon3")
lines(1-alpha, ES.boot.CI[2,], lty = "dotted", col = "maroon3")
eaxis(1)
eaxis(2)
mtext(substitute(B == B.~~"replications of size"~~n == n.~~"from Par("*th.*")",
                 list(B. = B, n. = n, th. = th)), side = 4, line = 1, adj = 0)
legend("bottomleft", bty = "n", lwd = c(1, rep(c(2, 1, 1.4, 1), 2)),
       lty = c("solid", rep(c("dashed", "solid", "dotdash", "dotted"), times = 2)),
       col = c("black", rep(c("royalblue3", "maroon3"), each = 4)),
       legend = c(
         expression("True"~VaR[alpha]~"and"~ES[alpha]),
         expression(widehat(VaR)[alpha]),
         expression("Bootstrapped"~~widehat(VaR)[alpha]),
         expression("Bootstrapped"~~Var(widehat(VaR)[alpha])),
         "Bootstrapped 95% CIs",
         expression(widehat(ES)[alpha]),
         expression("Bootstrapped"~~widehat(ES)[alpha]),
         expression("Bootstrapped"~~Var(widehat(ES)[alpha])),
         "Bootstrapped 95% CIs"),
       cex = 0.7)
dev.off()