evt.risk.07feb2017.R

# R code re: CapitalSpecator.com post for analyzing tail risk with extreme value theory in R
# "Tail-Risk Analysis In R: Part II - Extreme Value Theory"
# http://www.capitalspectator.com/tail-risk-analysis-in-r-part-ii-extreme-value-theory/
# 7 Feb 2017
# By James Picerno
# http://www.capitalspectator.com/
# (c) 2017 by Beta Publishing LLC


# load packages
library(xts)
library(timeSeries)
library(tseries)
library(TTR)
library(quantmod)
library(PerformanceAnalytics)
library(fExtremes)

# download data, generate returns
symbols <-c("VFINX", "VFITX")
getSymbols(symbols, src='yahoo', from = '1991-12-31', to = "2017-02-06")

for(symbol in symbols) {
  x <- get(symbol)
  indexFormat(x) <- '%Y-%m-%d'
  colnames(x) <- gsub("x", symbol, colnames(x))
  x <- x[,6]
  assign(symbol, x)
}

prices <- do.call(merge, lapply(symbols, get))
colnames(prices) <-c(symbols)

ret <-ROC(prices, 1, "discrete", na.pad = F)

# generate portfolio data
w = c(0.6, 0.4) # asset weights

port <-Return.portfolio(ret,
                        rebalance_on = "years",
                        weights = w,
                        wealth.index = TRUE,
                        verbose = TRUE)

port.ret <-ROC(port$wealthindex, 1, "discrete", na.pad = F)

# Density histogram: full sample 
hist(port.ret,
     freq = F,
     breaks = 100,
     col = 'gray99',
     xlab = "return",
     main = "Distribution of 1-Day Portfolio Returns",
     panel.first = grid() )
curve(dnorm(x,
            mean(port.ret),
            sd(port.ret)),
      add = TRUE,
      col = "red",
      lwd = 2)
legend("topleft",
       c("empirical", "normal"), 
       fill = c("gray99", "red"))
mtext(side = 3,
      "density histogram based on daily returns for 60/40 strategy",
      line = 0.2)
box()

# Density histogram: left tail
hist(port.ret,
     xlim = c(-0.05, -0.01),
     ylim = c(0, 2),
     xlab = "return",
     main = "Left Tail: Distribution of 1-year Portfolio Returns",
     freq = F, breaks = 100, col = 'gray90',
     panel.first = grid() )
curve(dnorm(x, mean(port.ret), 
            sd(port.ret) ),
      from = -0.05,
      to = -0.01,
      add = TRUE, 
      col = "red", 
      lwd = 3)
legend("topleft", c("empirical", "normal"), 
       fill = c("gray99", "red") )
mtext(side = 3,
      "density histogram based on daily returns for 60/40 strategy",
      line = 0.2)
box()

# Empirical cdf plot: full sample
plot.ecdf(as.numeric(port.ret),
          main = "Empirical Cumulative Distribution Function",
          xlab = "return",
          lwd = 2,
          panel.first = grid() )
curve(pnorm(x, mean(port.ret), sd(port.ret) ), 
      add = TRUE,
      lwd = 2,
      col='red')
mtext(side = 3,
      "based on daily returns for 60/40 strategy",
      line = 0.2)
legend("topleft",c("empirical","normal"), 
       fill = c("black","red")) 
box()

# Empirical cdf plot: left tail
plot.ecdf(as.numeric(port.ret),
          main = "Left Tail: Empirical Cumulative Distribution Function",
          xlab = "return",
          xlim = c(-0.05, -0.01), 
          ylim = c(0, 0.04),
          lwd = 2,
          panel.first = grid() )
curve(pnorm(x, mean(port.ret), sd(port.ret) ), 
      add = TRUE, 
      lwd = 2,
      col = 'red')
mtext(side = 3,
      "based on daily returns for 60/40 strategy",
      line = 0.2)
legend("topleft",c("empirical","normal"), 
       fill = c("black", "red") )
box()

# Estimate threshold

# Q-Q plot 
stats::qqnorm(port.ret, panel.first = grid() )
stats::qqline(port.ret, col = 'red')
abline(a = -0.01, b = 0, col = "blue")
mtext(side = 3,
        "based on daily returns for 60/40 strategy",
        line = 0.2)

# Compute threshold at 1% quantile
threshold.1 <- quantile(port.ret, 0.01)

# Fit evt model
gpd.fit.1 <-gpdFit(as.numeric(port.ret) * -1, 
                   u = threshold.1 * -1, 
                   type = c("mle") )

# Calculate risk metrics
var.es.mod <-round(tailRisk(gpd.fit.1, prob=c(0.95, 0.99, 0.999) ), 3)

var.his <-round(c(PerformanceAnalytics::VaR(port.ret, p = 0.95, method = c("historical") ),
                  PerformanceAnalytics::VaR(port.ret, p = 0.99, method = c("historical")),
                  PerformanceAnalytics::VaR(port.ret, p = 0.999, method = c("historical") )), 3)

es.his <-round(c(PerformanceAnalytics::ETL(port.ret, p = 0.95, 
                                           method=c("historical")),
                 PerformanceAnalytics::ETL(port.ret, p = 0.99, method = c("historical") ),
                 PerformanceAnalytics::ETL(port.ret, p = 0.999, method = c("historical") )), 3)

var.es.all <-cbind(var.es.mod[,1], var.es.mod[,-1] * -1, var.his, es.his) * 100

colnames(var.es.all) <-c("Confidence", "VaR.mod", "ES.mod", "VaR.his", "ES.his")

# Empirical cdf plot with model data
plot.ecdf(as.numeric(port.ret),
          main = "Left Tail: Empirical Cumulative Distribution Function",
          xlim = c(-0.05, -0.01), 
          ylim = c(0, 0.04),
          xlab = "return",
          lwd = 2,
          panel.first = grid() )
plot.ecdf(gpd.fit.1@data$x, add=TRUE,col='green3')
curve(pnorm(x, mean(port.ret), sd(port.ret)), add = TRUE, col = 'red')
legend("topleft", c("empirical", "normal", "model"),
       fill = c("black", "red", "green3") )
mtext(side = 3,
        "based on daily returns for 60/40 strategy",
        line = 0.2)
box()

# Similuate random data from GPD model
set.seed(215)
gpd.sim.1 <-fExtremes::rgpd(1000000, 
                            xi = gpd.fit.1@fit$par.ests[1], 
                            mu=0, 
                            beta = gpd.fit.1@fit$par.ests[2] )

# Histogram of simulated data
hist(gpd.sim.1,
     freq = F,
     main = paste("Simulated Left-Tail Density Histogram"),
     ylim = c(0, 0.01),
     xlim = c(0, 0.35),
     xlab = c("estimated daily loss"),
     breaks = 100,
     col = "cornflowerblue",
     panel.first = grid() )
mtext(side = 3,"based on 1 million simulations", line = 0.2)
box() 

# END