Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment