fragility <- function(x, rollWindow = 504, nEigen = 2, method = c("sdchng", "raw"),
sc = 252, halfLife = TRUE)
{
# FRAGILITY calculates the rolling absorption ratio of a timeSeries object
#
# ARGUMENTS:
# x = timeSeries object of asset returns
# rollWindow = rolling window to calculate correlation matirx
# nEigen = number of critical eigenvalues for ratio
# method = 'sdchg' transforms 'raw' output into a standardized change or z-score
# sc = window for z-score calculation
# halfLife = TRUE - use exponential weighting for covariance calculation
#
# OUTPUT:
# timeSeries of ratio
#
# NOTES:
# Principal Components as a Measure of Systemic Risk (Kritzman ,et al)
if(nEigen >= NCOL(x)) stop("nEigen needs to be less than the number of asset columns")
if (is.unsorted(x) == TRUE) stop("asset returns need to be sorted oldest to newest")
n = NROW(x)
rawAR = rep(NA, n-rollWindow-1)
HL = function(x) exp(-(x-1:x)*log(2)/(x/2))
for(i in rollWindow:n){
rollReturns = x[(i-(rollWindow-1)):i, ]
if (halfLife == TRUE) rollReturns = rollReturns * HL(rollWindow)
p = princomp(rollReturns)
compVar = p$sdev^2
rawAR[i-rollWindow-1] = sum(compVar[1:nEigen]/sum(compVar))
}
rawAR = as.data.frame(rawAR)
rownames(rawAR) = rownames(x[(rollWindow+2):NROW(x), ])
rawAR = as.timeSeries(rawAR)
if (method[1] == "raw"){
rawAR
}
else{
ARsc <- array(dim = NROW(rawAR)-sc-1)
j <- 1
for(i in sc:(NROW(rawAR)-1)){
ARt = rawAR[(i-sc):i, ]
ARsc[j] = (rawAR[i+1, ] - mean(ARt)) / sd(ARt)
j = j + 1
}
ARsc = data.frame(ARsc)
rownames(ARsc) = rownames(rawAR)[(sc+1):NROW(rawAR)]
as.timeSeries(ARsc)
}
}
turbulence <- function(x, rollWindow = 504, method = c("sdchng", "raw"), sc = 252)
{
# TURBULENCE calculates rolling mahalanobis distance of a timeSeries object
#
# ARGUMENTS:
# x = timeSeries object
# rollWindow = rolling window for mahalanobis distance calculation
# method: "sdchng" - standadized change of metric (z-score), or
# "raw" - raw metric
# sc: window for z-score calculation
#
# OUTPUT:
# timeSeries object of rolling mahalobis distance
#
# NOTES:
# Skulls, Financial Turbulence, and Risk Management (Kritzman and Lee)
n = NROW(x)
rawTurb = rep(NA, n-rollWindow-1)
for(i in rollWindow:n){
rollReturns = x[(i-(rollWindow-1)):i, ]
rollTurb = mahalanobis(rollReturns, mean(rollReturns), cov(rollReturns))
rawTurb[(i-rollWindow-1)] = rollTurb[NROW(rollTurb)]
}
rawTurb = as.data.frame(rawTurb)
rownames(rawTurb) = rownames(x[(rollWindow+2):NROW(x), ])
rawTurb = as.timeSeries(rawTurb)
if (method[1] == "raw"){
rawTurb
}
else{
scTurb <- array(dim = NROW(rawTurb)-sc-1)
j <- 1
for(i in sc:(NROW(rawTurb)-1)){
Turbt = rawTurb[(i-sc):i, ]
scTurb[j] = (rawTurb[i+1, ] - mean(Turbt)) / sd(Turbt)
j = j + 1
}
scTurb = data.frame(scTurb)
rownames(scTurb) = rownames(rawTurb)[(sc+1):NROW(rawTurb)]
as.timeSeries(scTurb)
}
}
require(fImport)
require(xts)
require(ggplot2)
require(plotly)
require(reshape2)
require(fAssets)
require(PerformanceAnalytics)
# get equity ETF returns for 38 countries
tix = c('EWO','EWK','EWC','EDEN','EWQ','EFNL','EWG',
'EWH','EIRL','EIS','EWI','EWJ','EWN','ENZL',
'ENOR','EWS','EWP','EWD','EWL','EWU','EUSA',
'ECH','EWZ','MCHI','INDA','EIDO','EWY',
'EWM','EWW','EPU','EPHE','EPOL','ERUS','EZA',
'EWT','THD','TUR')
nms = c('Australia','Belgium','Canada','Denmark','France',
'Finland','Germany','Hong Kong','Ireland','Israel',
'Italy','Japan','Netherlands','New Zealand','Norway',
'Singapore','Spain','Sweden','Switzerland','UK','US',
'Chile','Brazi','China','India','Indonesia',
'South Korea','Malaysia','Mexico','Peru','Philippines',
'Poland','Russia','South Africa','Taiwan','Thailand','Turkey')
price = yahooSeries(tix[1], from = '1900-01-01')[,6]
ret = returns(price, method = 'discrete')
for (i in 2:length(tix))
{
price = yahooSeries(tix[i], from = '1900-01-01')[,6]
ret = cbind(ret, returns(price, method = 'discrete'))
}
# use common time period
ret = na.omit(ret)
# calculate absorption ratio and fragility
AR = fragility(ret)
AR = xts(AR, as.Date(rownames(AR)))
turb = turbulence(ret)
turb = xts(turb, as.Date(rownames(turb)))
d = data.frame(index(AR), AR, turb)
colnames(d) = c('date', 'AR', 'Turb')
d = melt(d, id = 'date')
g = ggplot(d, aes(x = date, y = value, col = variable)) +
geom_line() + labs(col = '') + scale_y_continuous(breaks = -100:100) +
geom_hline(yintercept = 2, col = 'red', lty = 2) +
geom_hline(yintercept = -2, col = 'green', lty = 2) +
ylab('z-score') +
theme_bw()
ggplotly(g)
colarrange = assetsArrange(ret)
ret_arrange = ret[,colarrange]
cormat = cor(ret_arrange)
rownames(cormat) = tix
colnames(cormat) = tix
d = melt(cormat)
g = ggplot(d, aes(x = Var1, y = Var2, fill = value)) +
geom_tile() + scale_fill_distiller(palette = 'RdYlBu') +
xlab('') + ylab('') + ggtitle('Country Equity Correlation') +
theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())
ggplotly(g)
ret_t63 = ret[(NROW(ret)-63):NROW(ret), ]
colarrange = assetsArrange(ret_t63)
ret_arrange = ret[,colarrange]
cormat = cor(ret_arrange)
rownames(cormat) = tix
colnames(cormat) = tix
d = melt(cormat)
g = ggplot(d, aes(x = Var1, y = Var2, fill = value)) +
geom_tile() + scale_fill_distiller(palette = 'RdYlBu') +
xlab('') + ylab('') + ggtitle('Country Equity Correlation Last 63 Trading Days') +
theme(axis.ticks.x = element_blank(), axis.text.x = element_blank())
ggplotly(g)
acwi = yahooSeries('acwi', from = '1900-01-01')[,6]
acwi = returns(acwi, method = 'discrete')
acwi = xts(acwi, as.Date(row.names(acwi)))
acwi_t = acwi['2015-02-10/']
dd = drawdowns(as.timeSeries(acwi_t))
dd = as.xts(dd)
rollVol = rollapply(acwi, width = 63, 'sd')
rollVol = rollVol * sqrt(252)
rollVol = rollVol['2015-02-10/']
d = data.frame(index(acwi_t), dd, rollVol)
colnames(d) = c('date', 'Loss', 'Vol')
d = melt(d, id = 'date')
g = ggplot(d, aes(x = date, y = value, col = variable)) +
geom_line() + labs(col = '') + scale_y_continuous(labels = scales::percent) +
ylab('') + xlab('') +
theme_bw()
ggplotly(g)