Friday, August 31, 2012

Russell 2000 Softail Fat Boy

If the Russell 2000 were a motorcycle, maybe it should be a Harley-Davidson Softail Fat Boy.  I have explored the exceptional case of the Russell 2000 in quite a few posts

but I still am not sure I have done a good job of clearly and simply explaining some of the unique characteristics of the Russell 2000 over the last  25 years.  The Russell 2000 has been extremely difficult to beat because the upside has been so good (fat boy) and the downside has been the same or better (softtail).

From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio

R code from GIST:

#explore exceptional case of the Russell 2000
require(quantmod)
require(PerformanceAnalytics)
require(xtsExtra)
getSymbols("^RUT", from = "1900-01-01")
getSymbols("^GSPC", from = "1900-01-01")
#do initial exploration of distribution
chart.QQPlot(na.omit(ROC(RUT[,4],type="discrete",n=250))["1987::"])
chart.QQPlot(na.omit(ROC(GSPC[,4],type="discrete",n=250))["1987::"])
#explore barplot.xts to do a chart of annual returns for both indexes
#merge prices
prices <- merge(GSPC[,4],RUT[,4])
#use endpoints to get annual returns
returns.annual <- as.xts(apply(
ROC(prices[endpoints(prices,"years")],type="discrete",n=1),
MARGIN = 2,
FUN = na.fill, fill = 0),
order.by = index(prices[endpoints(prices,"years")]))
#name columns something a little more clear
colnames(returns.annual) <- c("S&P 500","Russell 2000")
#using barplot.xts create the plot
#I made some subtle changes to barplot.xts to experiment so plot will be cosmetically different
barplot.xts(returns.annual,
stacked=FALSE,
box="transparent", #get rid of box surrounding the plot
ylim=c(-0.5,0.5),
ylab=NA,
border=c(brewer.pal(n=11,"BrBG")[c(4,9)]),
col=c(brewer.pal(n=11,"BrBG")[c(4,9)])) #deliberately trying some new colors
title(main="Annual Returns of S&P 500 and Russell 2000",
outer = TRUE,
adj=0.05, font.main = 1, cex.main = 1.25, line = -2)
require(latticeExtra)
require(reshape2)
roc <- na.omit(merge(ROC(GSPC[,4],type="discrete",n=250),ROC(RUT[,4],type="discrete",n=250)))
#name columns something a little more clear
colnames(roc) <- c("S&P 500","Russell 2000")
roc.melt <- melt(coredata(roc))
asTheEconomist(
densityplot(~value,data=roc.melt,groups=Var2,
#par.settings=theEconomist.theme(box="transparent"),
#lattice.options=theEconomist.opts(),
auto.key=list(space="right",col=(brewer.pal(n=11,"BrBG")[c(9,4)]),lines=FALSE),
col=(brewer.pal(n=11,"BrBG")[c(9,4)]),
ylab=NA,
main="Annual Returns (Rolling 250 Day) of S&P 500 and Russell 2000")
)

Thursday, August 30, 2012

Another Great Google Summer of Code 2012 R Project

Tradeblotter announced the very nice features that will be added to the PerformanceAnalytics package as a result of the Google Summer of Code (GSOC) 2012 project:

“…Matthieu commenced to produce dozens of new functions, extend several more existing ones, and add more than 40 pages of additional documentation (complete with formulae and examples) to PerformanceAnalytics. He’s included Bacon’s small data set and several new table.* functions for testing and demonstrating that the functions match the published results. All on plan, I would add.

He also wrote a very nice overview of the functions developed from Bacon (2008) that are included in PerformanceAnalytics, which should be helpful to readers or teachers of Bacon’s work. Matthieu’s summary document will also be distributed as a vignette in the package, accessible using vignette('PA-Bacon'). Additional detail is included in the documentation for each function, as well. I’ll highlight some of those functions in later posts…”

Although the Bacon vignette is very helpful, I highly recommend reading the book also.  For those of you on O’Reilly’s Safari, you can add the book to your bookshelf.

Preview

or buy it from Amazon.

Since I saw the announcement (been patiently waiting) when I was also reading this very good article “Tomatoes and the Low Vol Effect” by Research Affiliates, I was already working on some R to explore the Sharpe and Information Ratio between a buy/hold on the S&P 500 and a 10 month moving average strategy.  When you look at the code, you might enjoy some of the other charts and methods, but  I decided to change directions with my research and integrate a new PerformanceAnalytics function and continue to push the limits with another fine Google Summer of Code 2012 addition xtsExtra plot.xts.  I was not real selective and opted for ProspectRatio, since it is a formula I do not often hear discussed.  On page 20, the Bacon vignette describes ProspectRatio as

image

I basically tried to include absolutely everything possible (scattered, smothered, covered, chunked, and diced) from xtsExtra plot.xts in this chart.  The chart includes cumulative growth of buy/hold and 10 month moving average on the S&P 500 in the top panel, shows drawdown in the middle, and then adds a horizon plot of Prospect Ratio for the ma strategy – Prospect Ratio for the buy/hold.  Since plot.xts also allows for blocks, I thought it would be helpful to also shade those periods when the 36 month rolling Sharpe Ratio of the ma strategy outperforms buy/hold.  While this is probably a bit much for a single chart, I think it demonstrates clearly how powerful these tools are.

From TimelyPortfolio

R code from GIST (do raw for copy/paste):

#integrate 2 Google Summer of Code 2012 Projects
#plot.xts and PerformanceAnalytics both received very nice additions
#wish I knew the exact link but believe this section came from stackoverflow
#this allows you to work with the code straight from r-forge SVN
## If you want to source() a bunch of files, something like
## the following may be useful:
#path="C:\\Program Files\\R\\R-2.15.1\\sandbox\\svnsource\\returnanalytics\\pkg\\PerformanceAnalytics\\r"
#sourceDir <- function(path, trace = TRUE, ...) {
# for (nm in list.files(path, pattern = "\\.[RrSsQq]$")) {
# if(trace) cat(nm,":")
# source(file.path(path, nm), ...)
# if(trace) cat("\n")
# }
#}
#sourceDir(path)
require(RColorBrewer)
require(quantmod)
require(xtsExtra)
getSymbols("^GSPC",from="1890-01-01")
sp500.monthly <- GSPC[endpoints(GSPC, "months"),4]
roc <- ROC(sp500.monthly, type = "discrete", n = 1)
n = 10 #set n for number of periods; this is 10 for 10 months
roc.ma <- lag(ifelse(sp500.monthly > runMean(sp500.monthly, n = n), 1, 0), k = 1) * roc
returns <- merge(roc, roc.ma)
returns <- as.xts(apply(returns, MARGIN = 2, na.fill, fill = 0), order.by = index(returns))
colnames(returns) <- c("SP500.buyhold", "SP500.ma")
charts.PerformanceSummary(returns, ylog = TRUE)
sr <- SharpeRatio.annualized(returns)
ir <- InformationRatio(returns, returns[,1])
applyacross <- function(x,rollFun=Omega,width=12,by=1,Rb,...) {
if(missing(Rb)) { #add this so we can also use for InformationRatio
result <- apply.rolling(x,FUN=rollFun,width,by,...)
} else {
result <- apply.rolling(x,FUN=rollFun,width,by,Rb=Rb,...)
}
result <- merge(x,result)[,2] #merge to pad with NA and then get second column
colnames(result) <- colnames(x)
return(result)
}
stat = "SharpeRatio.annualized" #ES, KellyRatio, Omega, skewness, kurtosis, VaR, SharpeRatio, CalmarRatio
width = 36
sharpe.rolling <- as.xts(matrix(unlist(lapply(returns,FUN=applyacross,rollFun=get(stat),width=width)),byrow=FALSE,nrow=NROW(returns)), order.by = index(returns))
sharpe.rolling <- as.xts(apply(sharpe.rolling,MARGIN=2,FUN=na.fill,fill=0),order.by=index(sharpe.rolling))
colnames(sharpe.rolling) <- colnames(returns)
plot.xts(sharpe.rolling,screens=1)
plot.xts(sharpe.rolling[,2] - sharpe.rolling[,1],screens=1)
plot(x=coredata(sharpe.rolling[,1]),y=coredata(sharpe.rolling[,2]),pch=19,
col=ifelse(sharpe.rolling[,2]>0,"green","red"))
abline(lm(SP500.ma ~ SP500.buyhold, data=as.data.frame(sharpe.rolling)))
text(x=coredata(sharpe.rolling[,1]),y=coredata(sharpe.rolling[,2]),labels=format(as.Date(index(sharpe.rolling)),"%Y"),cex=0.5, pos = 2 )
plot(lm(SP500.ma ~ SP500.buyhold, data=as.data.frame(sharpe.rolling)),which=2)
stat = "InformationRatio"
#this allows multiple columns if we have more than one index or comparison
ir.rolling <- as.xts(matrix(unlist(lapply(returns,FUN=applyacross,rollFun=get(stat),width=width,Rb=returns[,1])),byrow=FALSE,nrow=NROW(returns)), order.by = index(returns))
colnames(ir.rolling) <- colnames(returns)
#could also do this if we only have one column
#ir.rolling <- apply.rolling(returns[,2],FUN=InformationRatio,width=width,by=1,Rb=returns[,1])
plot.xts(na.omit(merge(ir.rolling[,2],ROC(sp500.monthly,type="discrete",n=36))),screens=1)
plot(x=coredata(ROC(sp500.monthly,type="discrete",n=36)),y=coredata(ir.rolling[,2]),
pch=19,
col=ifelse(ir.rolling[,2]>0,"green","red"),
las=1)
abline(h=0)
abline(v=0)
plot(x = coredata(sharpe.rolling[,1]),y=coredata(ir.rolling[,2]),
pch = 19,
col = ifelse(ir.rolling[,2]>0,"green","red"),
las = 1)
abline(h=0)
abline(lm(coredata(ir.rolling[,2])~coredata(sharpe.rolling[,1])))
plot(lm(coredata(ir.rolling[,2])~coredata(sharpe.rolling[,1])))#,which=2)
#use the new ProspectRatio function in PerfomranceAnalytics
stat = "ProspectRatio" #ES, KellyRatio, Omega, skewness, kurtosis, VaR, SharpeRatio, CalmarRatio
width = 36
pr.rolling <- as.xts(matrix(unlist(lapply(returns,FUN=applyacross,rollFun=get(stat),width=width,MAR=0.025)),byrow=FALSE,nrow=NROW(returns)), order.by = index(returns))
pr.rolling <- as.xts(apply(pr.rolling,MARGIN=2,FUN=na.fill,fill=0),order.by=index(pr.rolling))
colnames(pr.rolling) <- colnames(returns)
#set up horizon plot functionality
horizon.panel <- function(index,x,...) {
#get some decent colors from RColorBrewer
#we will use colors on the edges so 2:4 for red and 7:9 for blue
require(RColorBrewer)
col.brew <- brewer.pal(name="RdBu",n=10)
#ease this reference later
n=NROW(x)
#clean up NA with either of the two methods below
#x[which(is.na(x),arr.ind=TRUE)[,1],
# unique(which(is.na(x),ar.ind=TRUE)[,2])] <- 0
x <- apply(x,MARGIN=2,FUN=na.fill,fill=0)
#get number of bands for the loop
#limit to 3
nbands = 3
#first tried this but will not work since each series needs to have same number of bands
#min(4,ceiling(max(abs(coredata(x)))/horizonscale))
par(usr=c(index[1],par("usr")[2],origin,horizonscale))
for (i in 1:nbands) {
#draw positive
polygon(
c(index[1], index, index[n]),
c(origin, coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[length(col.brew)-nbands+i-1],
border=NA
)
#draw negative
polygon(
c(index[1], index, index[n]),
c(origin, -coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[nbands-i+1],
border=NA
)
}
#delete trash drawn below origin that we keep so no overlap between positive and negative
polygon(
c(index[1], index, index[n]),
c(origin, -ifelse(coredata(x)==origin,horizonscale*5,abs(coredata(x))),origin),
col=par("bg"),
border=NA
)
#draw a line at the origin
abline(h=origin,col="black")
#draw line at top of plot or otherwise polygons will cover boxes
abline(h=par("usr")[4],col="black")
#mtext("ProspectRatio Difference", side = 3, adj = 0.02, line = -1.5, cex = 0.75)
}
horizonscale = 0.25
origin = 0
#trying this to color sections or color lines based on sharpe
rle.sharpe <- rle(as.vector(sharpe.rolling[,2]>sharpe.rolling[,1]))
dates <- index(returns)[cumsum(rle.sharpe$lengths)]
start.i=ifelse(na.omit(rle.sharpe$values)[1],2,1)
#png("plotxts with everything and ProspectRatio.png",height=600, width=640)
plot.xts(merge(log(cumprod(1+returns)),Drawdowns(returns),pr.rolling[,2] - pr.rolling[,1]),
screens = c(1,1,2,2,3), #since 2 columns for cumul and drawdown repeat screens
layout.screens = c(1,1,1,1,2,2,3), #make screen 1 4/7 of total 2 2/7 and 3 (horizon) 1/7
col = brewer.pal(9,"Blues")[c(5,8)], #get two blues that will look ok
lwd = c(1.5,2), #line width; will do smaller 1.5 for benchmark buy/hold
las = 1, #do not rotate y axis labels
ylim = matrix(c(0,5,-0.55,0,origin,horizonscale),byrow=TRUE,ncol=2), #plot.xts accepts ylim in matrix form; print matrix to see how it works
auto.legend = TRUE, #let plot.xts do the hard work on the legend
legend.loc = c("topleft",NA, NA), #just do legend on the first screen
legend.pars = list(bty = "n", horiz=TRUE), #make legend box transparent and legend horizontal
panel = c(default.panel,default.panel,horizon.panel), #specify panels for each screen
main = NA, #will do title later so we have more control
#log="y", #log scale does not work with blocks
blocks = list(start.time=dates[seq(start.i,NROW(dates),2)], #overlay blocks in which 36-mo sharpe ratio of ma exceeds buy/hold
end.time=dates[-1][seq(start.i,NROW(dates),2)],col="gray90")) #darkolivegreen2"))
title(main = "Strategy Comparison on S&P 500 - Buy Hold versus Moving Average", adj = 0.05, line = -1.5, outer = TRUE, cex.main = 1.1, font.main = 3)
#dev.off()

Monday, August 27, 2012

Horizon on ggplot2

SocialDataBlog’s kind reference in post Horizon plots with ggplot (not) motivated me to finish what the post started.  I knew that ggplot2 would be a little more difficult to use for the purpose of a horizon plot, but I felt compelled to provide at least one example of a horizon plot for each of the major R graphing packages.  I achieved a good result but the code is not as elegant or as flexible as I would like.  Readers more comfortable with ggplot2, please bash, fork, and improve.

From TimelyPortfolio

R code in GIST (do raw for copy/paste):

#first attempt at implementing horizon plots in ggplot2
#pleased with result but code sloppy and inflexible
#as always very open to improvements and forks
require(ggplot2)
require(reshape2)
require(quantmod)
require(PerformanceAnalytics)
require(xtsExtra)
data(edhec)
origin = 0
horizonscale = 0.1
#get 12 month rolling return of edhec indexes
roc <- as.xts(apply(cumprod(edhec+1),MARGIN=2,ROC,n=12,type="discrete"),order.by=index(edhec))
roc.df <- as.data.frame(cbind(index(roc),coredata(roc)))
roc.melt <- melt(roc.df,id.vars=1)
roc.melt[,1] <- as.Date(roc.melt[,1]) #convert back to a Date
horizon.panel.ggplot <- function(df, title) {
#df parameter should be in form of date (x), grouping, and a value (y)
colnames(df) <- c("date","grouping","y")
#get some decent colors from RColorBrewer
#we will use colors on the edges so 2:4 for red and 7:9 for blue
require(RColorBrewer)
col.brew <- brewer.pal(name="RdBu",n=10)
#get number of bands for the loop
#limit to 3 so it will be much more manageable
nbands = 3
#loop through nbands to add a column for each of the positive and negative bands
for (i in 1:nbands) {
#do positive
df[,paste("ypos",i,sep="")] <- ifelse(df$y > origin,
ifelse(abs(df$y) > horizonscale * i,
horizonscale,
ifelse(abs(df$y) - (horizonscale * (i - 1) - origin) > origin, abs(df$y) - (horizonscale * (i - 1) - origin), origin)),
origin)
#do negative
df[,paste("yneg",i,sep="")] <- ifelse(df$y < origin,
ifelse(abs(df$y) > horizonscale * i,
horizonscale,
ifelse(abs(df$y) - (horizonscale * (i - 1) - origin) > origin, abs(df$y) - (horizonscale * (i - 1) - origin), origin)),
origin)
}
#melt data frame now that we have added a column for each band
#this will fit ggplot2 expectations and make it much easier
df.melt <- melt(df[,c(1:2,4:9)],id.vars=1:2)
#name the columns for reference
#try to be generic
colnames(df.melt) <- c("date","grouping","band","value")
#use ggplot to produce an area plot
p <- ggplot(data=df.melt) +
geom_area(aes(x = date, y = value, fill=band),
#alpha=0.25,
position="identity") + #this means not stacked
scale_fill_manual(values=c("ypos1"=col.brew[7], #assign the colors to each of the bands; colors get darker as values increase
"ypos2"=col.brew[8],
"ypos3"=col.brew[9],
"yneg1"=col.brew[4],
"yneg2"=col.brew[3],
"yneg3"=col.brew[2])) +
ylim(origin,horizonscale) + #limit plot to origin and horizonscale
facet_grid(grouping ~ .) + #do new subplot for each group
theme_bw() + #this is optional, but I prefer to default
opts(legend.position = "none", #remove legend
strip.text.y = theme_text(),#rotate strip text to horizontal
axis.text.y = theme_blank(),#remove y axis labels
axis.ticks = theme_blank(), #remove tick marks
axis.title.y = theme_blank(),#remove title for the y axis
axis.title.x = theme_blank(),#remove title for the x axis
title = title, #add a title from function parameter
plot.title = theme_text(size=16, face="bold", hjust=0)) #format title
return(p)
}
horizon.panel.ggplot(roc.melt, "EDHEC Indexes Return (Rolling 1 Year)")

Thursday, August 23, 2012

Bonds Much Sharpe -r Than Buffett

Mebane Faber’s post Buffett’s Alpha points out Warren Buffett’s 0.76 Sharpe Ratio discussed in the similarly title paper Buffet’s Alpha.  I of course immediately think about the 8th Wonder of the World – the US Bond Market, whose Sharpe Ratio has trounced Buffett’s for the last 30 years.  What I like even better are all the tactical systems that employ US bonds in their backtests and make no adjustment for substantially different returns going forward.  For those of you who do not know, bonds with absolute certainty cannot achieve >8% annualized returns with max drawdown < 5% for the next 30 years with a starting yield to worst at 1.86% (Barclays Agg 8/23/2012). 

If anyone can show me where to get 8% annualized returns with max drawdown of 5% for the next 30 years, please let me know, and I will buy that with leverage and enjoy life. I’ll be happy to share my gains with whomever has the answer.

From TimelyPortfolio

In addition to an unbelievable Sharpe Ratio, bonds have exhibited a low/negative correlation with stocks during stocks’ bear market, which is also historically very anomalous.

From TimelyPortfolio

Just because I love horizon plots.

From TimelyPortfolio

For a longer perspective, here is rolling correlation since 1900 from the CSFB 2011 Yearbook.

image

So the experience can be dramatically different than what we have been fortunate enough to experience recently.

image

R code from GIST (install xtsExtra):

require(quantmod)
require(PerformanceAnalytics)
require(xtsExtra)
require(RColorBrewer)
#unfortunately don't feel like fighting IP lawyers so I cannot share this index data
portfolio <- read.csv("file.csv",stringsAsFactors=FALSE)
portfolio <- portfolio[2:NROW(portfolio),2:NCOL(portfolio)]
portfolio <- portfolio[,c(1,3,5)]
#since export has duplicate colnames we need to remove the .1 added
#colnames(portfolio) <- substr(colnames(portfolio),1,nchar(colnames(portfolio))-2)
len <- nchar(portfolio[,1])
xtsdate <- paste(substr(portfolio[,1],len-3,len),"-",
ifelse(len==9,"0",""),substr(portfolio[,1],1,len-8),"-01",sep="")
portfolio.xts <- xts(data.matrix(portfolio[,2:NCOL(portfolio)]),order.by=as.Date(xtsdate))
portfolio.xts <- portfolio.xts/100
portfolio.xts[1,]<-0
charts.RollingPerformance(portfolio.xts[,1], #Barclays Aggregate
width = 60, #rolling 5 year (60 months)
Rf = portfolio.xts[,2], #Barclays 3-month
main = "Barclays Aggregate Performance (Rolling 5 Year)")
getSymbols("SP500", src = "FRED")
getSymbols("GS10", src = "FRED")
require(RQuantLib)
GS10pricereturn<-GS10
GS10pricereturn[1,1]<-0
colnames(GS10pricereturn)<-"US10yPrice"
for (i in 1:(NROW(GS10)-1)) {
GS10pricereturn[i+1,1]<-FixedRateBondPriceByYield(yield=GS10[i+1,1]/100,issueDate=Sys.Date(),
maturityDate= advance("UnitedStates/GovernmentBond", Sys.Date(), 10, 3),
rates=GS10[i,1]/100,period=2)[1]/100-1
}
#total return will be the price return + yield/12 for one month
GS10totalreturn<-GS10pricereturn+lag(GS10,k=1)/12/100
colnames(GS10totalreturn)<-"US10yTotalReturn"
#thanks to PerformanceAnalytics package for this handy rolling correlation calculation
rollcorr <- rollapplyr(na.omit(merge(ROC(SP500,n=1,type="discrete"),
GS10totalreturn)),
width = 60, #5 year rolling
FUN = function(x) cor(x[, 1, drop = FALSE],
x[, 2, drop = FALSE]),
by = 1,
by.column = FALSE,
na.pad = TRUE)
rollcorr[is.na(rollcorr)] <- 0 #change leading NAs to 0
custom.panel <- function (index, x, col, lwd, ...) {
default.panel(index, ifelse(x > 0, x, 0), col = brewer.pal("Greens", n = 9)[7], lwd = 2, ...)
default.panel(index, ifelse(x < 0, x, 0), col = brewer.pal("Reds",n = 9)[7], lwd = 2, ...)
abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lwd=0.5)
abline (h = 0, col = "black", lwd = 2)
}
plot.xts(rollcorr,
auto.grid = FALSE, #turn off grid for more precise control
bty = "n", #turn off box because I like better
las = 1,
ylim = c(-1,1),
main = NA,
panel = custom.panel, #use our customized panel
type = "l",
major.format = "%Y",
minor.tick = FALSE,
yaxs = "i",
cex.axis = 0.8,
lwd = 2,
blocks = list(start.time = "2000-03-01", end.time = "2012-12-31"))
title(main = "Correlation of US 10y Bonds with S&P 500 (Rolling 5 Year)", outer = TRUE, line = -2, adj = 0.05)
mtext("Source: Federal Reserve Bank of St. Louis (FRED)", side = 1, cex = 0.7, font = 3, adj = 0.9, outer = TRUE, line = -2)
horizon.panel <- function(index,x,...) {
#get some decent colors from RColorBrewer
#we will use colors on the edges so 2:4 for red and 7:9 for blue
require(RColorBrewer)
col.brew <- brewer.pal(name="RdBu",n=10)
#ease this reference later
n=NROW(x)
#clean up NA with either of the two methods below
#x[which(is.na(x),arr.ind=TRUE)[,1],
# unique(which(is.na(x),ar.ind=TRUE)[,2])] <- 0
x <- apply(x,MARGIN=2,FUN=na.fill,fill=0)
#get number of bands for the loop
#limit to 3
nbands = 3
#first tried this but will not work since each series needs to have same number of bands
#min(4,ceiling(max(abs(coredata(x)))/horizonscale))
par(usr=c(index[1],index[n],origin,horizonscale))
for (i in 1:nbands) {
#draw positive
polygon(
c(index[1], index, index[n]),
c(origin, coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[length(col.brew)-nbands+i-1],
border=NA
)
#draw negative
polygon(
c(index[1], index, index[n]),
c(origin, -coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[nbands-i+1],
border=NA
)
}
#delete trash drawn below origin that we keep so no overlap between positive and negative
polygon(
c(index[1], index, index[n]),
c(origin, -ifelse(coredata(x)==origin,horizonscale*5,abs(coredata(x))),origin),
col=par("bg"),
border=NA
)
#draw a line at the origin
abline(h=origin,col="black")
}
horizonscale = 0.25
origin = 0
plot.xts(rollcorr,
auto.grid = FALSE, #turn off grid for more precise control
bty = "n", #turn off box because I like better
las = 1,
ylim = c(origin,horizonscale),
main = NA,
panel = horizon.panel, #use our customized panel
type = "l",
major.format = "%Y",
minor.tick = FALSE,
yax.loc = "none",
yaxt = "n", xaxt = "n",
cex.axis = 0.6,
lwd = 2)
title(main = "Correlation of US 10y Bonds with S&P 500 (Rolling 5 Year)", outer = TRUE, line = -2, adj = 0.05, cex = 0.7)
mtext("Source: Federal Reserve Bank of St. Louis (FRED)", side = 1, cex = 0.7, font = 3, adj = 0.9, outer = TRUE, line = -2)

Monday, August 20, 2012

plot.xts with Moving Average Panel

(for all plot.xts posts, see http://timelyportfolio.blogspot.com/search/label/plot.xts)

As another example of all that we can do with the new plot.xts, let’s try to do a price plot with a moving average overlays.  We will use the ETFs shown by Mebane Faber at http://www.mebanefaber.com/timing-model/.  With the panel functionality, it is very easy to specify a panel to draw the price line and then add the calculated moving average.  Notice how in all the examples, the recession block appears easily and very nicely.

From TimelyPortfolio

Also, if you wanted to specify some funky layouts, we have that option.  For this case, I do not think it makes much sense, but in the future I will demonstrate some more appropriate uses.

From TimelyPortfolio

Or as one more example, let’s change it up just slightly.

From TimelyPortfolio

R code in GIST (do raw for copy/paste and also install xtsExtra from r-forge):

#install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
require(quantmod)
require(RColorBrewer)
require(xtsExtra) #if you get an error, see first line and install from r-forge
#use Mebane Faber tickers VTI, VEU, IEF, VNQ, DBC
#as discussed in http://papers.ssrn.com/sol3/papers.cfm?abstract_id=962461
tckrs <- c("VTI", "VEU", "IEF", "VNQ", "DBC")
getSymbols(tckrs, from = "2000-01-01")
prices <- get(tckrs[1])[,6]
for (i in 2:length(tckrs)) {
prices <- na.omit(merge(prices, get(tckrs[i])[,6]))
}
ma.panel <- function (index, x, col, ...) {
#draw line for price
default.panel(index, x, col, ...)
#label each panel with first 3 characters of column name
mtext(substr(colnames(x), 1, 3), side = 3, cex = 0.8, line = -2.5, adj = 0.5, col = col)
#now get n=200 moving average
ma <- runMean(x, n = 200)
#add the moving average line
default.panel(index, ma, col="indianred1", ...)
abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lty=3,lwd=0.5)
#abline(h=par("yaxp")[1], col="black")
}
plot.xts(prices,
las = 1, #no rotate for y axis tickmark labels
panel = ma.panel, #do the moving average panel to overlay moving average lines
auto.grid = FALSE, #get some extra contol of grid
col = brewer.pal("Set1", n = 7)[c(2:5,7)],
main = NA,
blocks = list(start.time = "2007-12-01", end.time = "2009-06-01", col="lightblue"))
title(main = "ETFs (www.mebanefaber.com) and 200 day Moving Average", adj = 0.05, outer = TRUE, line = -2)
mtext("source: Yahoo! Finance", outer = TRUE, side = 1, adj = 0.94, font = 1, cex = 0.7, line = -1)
#not sure why you would want to do something like this
#but for example purposes
#we can change the layout radically with plot.xts
plot.xts(prices,
screens = c(1,2,4,3,5),
layout.screens = matrix(c(1,1,2,2,3,3,4,5), nrow = 2, byrow = FALSE), #to see how the matrix works print the matrix print(matrix(c(1,1,2,2,3,3,4,5), nrow = 2, byrow = FALSE))
ylim = matrix(c(range(prices[,c(1,2,4)]),range(prices[,c(1,2,4)]),range(prices[,c(1,2,4)]),range(prices[,3]),range(prices[,5])),
byrow=TRUE,ncol=2),
las = 1, #no rotate for y axis tickmark labels
panel = ma.panel, #do the moving average panel to overlay moving average lines
auto.grid = FALSE, #get some extra contol of grid
col = brewer.pal("Set1", n = 7)[c(2:5,7)],
main = NA,
minor.ticks = FALSE,
major.format = "%Y",
blocks = list(start.time = "2007-12-01", end.time = "2009-06-01", col="lightblue"))
title(main = "ETFs (www.mebanefaber.com) and 200 day Moving Average", adj = 0.05, outer = TRUE, line = -2)
mtext("source: Yahoo! Finance", outer = TRUE, side = 1, adj = 0.94, font = 1, cex = 0.7, line = -1)
plot.xts(prices,
screens = c(1,2,4,3,5),
layout.screens = matrix(c(1,1,2,2,3,3,4,5,4,5), nrow = 2, byrow = FALSE), #to see how the matrix works print the matrix print(matrix(c(1,1,2,2,3,3,4,5,4,5), nrow = 2, byrow = FALSE))
ylim = matrix(c(range(prices[,c(1,2,4)]),range(prices[,c(1,2,4)]),range(prices[,c(1,2,4)]),range(prices[,3]),range(prices[,5])),
byrow=TRUE,ncol=2),
las = 1, #no rotate for y axis tickmark labels
panel = ma.panel, #do the moving average panel to overlay moving average lines
auto.grid = FALSE, #get some extra contol of grid
col = brewer.pal("Set1", n = 7)[c(2:5,7)],
main = NA,
minor.ticks = FALSE,
major.format = "%Y",
blocks = list(start.time = "2007-12-01", end.time = "2009-06-01", col="lightblue"))
title(main = "ETFs (www.mebanefaber.com) and 200 day Moving Average", adj = 0.05, outer = TRUE, line = -2)
mtext("source: Yahoo! Finance", outer = TRUE, side = 1, adj = 0.94, font = 1, cex = 0.7, line = -1)

Friday, August 17, 2012

GARCH Panel in plot.xts

I’m clearly out of my realm of competence with most of the rugarch functions, but I thought it might be nice to provide an example combining plot.xts and uGARCHroll.

image

R code from GIST:

#install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
require(quantmod)
require(rugarch)
require(xtsExtra) #if you get an error, see first line and install from r-forge
getSymbols("DEXJPUS",src="FRED")
DEXJPUS<-1/to.weekly(DEXJPUS)
ugarch.panel <- function(index,x,type,cex,col,pch,...){
spec = ugarchspec(
variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(1,1), include.mean=T))
#get ugarchroll; I cannot say I completely understand what I'm doing here
ugr <- ugarchroll(spec,
data = na.omit(x),
forecast.length = floor(NROW(na.omit(x)) / 1000) * 1000,
n.ahead = 1,
refit.every = 50,
refit.window = "moving")
#get garch coefficients
ugr.var <- merge(x,as.data.frame(ugr,which="VaR"))
ugr.var <- as.xts(apply(ugr.var,MARGIN=2,na.fill,fill=c(0,"extend")),order.by=index(x))[,2:4]
print(tail(ugr.var))
default.panel(index,ugr.var[,3],type="h",cex,pch,col,...)
default.panel(index,ugr.var[,1],type="l",cex,pch,col="red",...)
default.panel(index,ugr.var[,2],type="l",cex,pch,col="gray50",...)
text(x=index[1],y=par("usr")[4],"VaR from rugarch ugarchroll",adj=c(0,1),cex=0.8)
}
plot.xts(na.omit(merge(DEXJPUS[,4],ROC(DEXJPUS[,4],n=1,type="discrete"))),
screens=c(1,2),
minor.ticks=FALSE, major.format="%Y",
panel=c(default.panel,ugarch.panel),
main="Japanese Yen (source: St. Louis Federal Reserve)")

Horizon Plots with plot.xts

Anyone who has read

  1. 48 Industries (Dendrogram Ordered) Over 50 Years
  2. 48 Industries Since 1963
  3. “Trend is Not Your Friend” Applied to 48 Industries
  4. Horizon Plots in Base Graphics
  5. More on Horizon Charts
  6. Application of Horizon Plots
  7. Horizon Plot Already Available
  8. Cubism Horizon Charts in R

should already know that I really like horizon charts.  Also, you would probably guess that I would be very excited if the new plot.xts could produce horizon charts.  So I am excited, since with its panel functionality, plot.xts very capably creates horizon charts.

From TimelyPortfolio

R code from GIST (also, please install xtsExtra from r-forge):

#plot.xts with horizons
#install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
require(PerformanceAnalytics)
require(quantmod)
require(xtsExtra) #if you get error, install xtsExtra from r-forge
horizon.panel <- function(index,x,...) {
#get some decent colors from RColorBrewer
#we will use colors on the edges so 2:4 for red and 7:9 for blue
require(RColorBrewer)
col.brew <- brewer.pal(name="RdBu",n=10)
#ease this reference later
n=NROW(x)
#clean up NA with either of the two methods below
#x[which(is.na(x),arr.ind=TRUE)[,1],
# unique(which(is.na(x),ar.ind=TRUE)[,2])] <- 0
x <- apply(x,MARGIN=2,FUN=na.fill,fill=0)
#get number of bands for the loop
#limit to 3
nbands = 3
#first tried this but will not work since each series needs to have same number of bands
#min(4,ceiling(max(abs(coredata(x)))/horizonscale))
for (i in 1:nbands) {
#draw positive
polygon(
c(index[1], index, index[n]),
c(origin, coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[length(col.brew)-nbands+i-1],
border=NA
)
#draw negative
polygon(
c(index[1], index, index[n]),
c(origin, -coredata(x) - (i-1) * horizonscale,origin),
col=col.brew[nbands-i+1],
border=NA
)
}
#delete trash drawn below origin that we keep so no overlap between positive and negative
polygon(
c(index[1], index, index[n]),
c(origin, -ifelse(coredata(x)==origin,horizonscale*5,abs(coredata(x))),origin),
col=par("bg"),
border=NA
)
#draw a line at the origin
abline(h=origin,col="black")
#labelling just for example purposes
#do label at top left and top right
text(x=index[1], y=horizonscale, colnames(x), adj=c(0,1))
text(x=index[n], y=horizonscale, colnames(x), adj=c(1,1))
#do label at center right for ease of reference
#text(x=index[n], y=horizonscale/2, colnames(x), pos=2)
#do label at bottom center
#text(x=index[n/2], y=origin, colnames(x), pos=3)
}
data(edhec)
#set these up for global access since I do not see how to pass into panel function
origin = 0
horizonscale = 0.1
#get 12 month rolling return of edhec indexes
roc <- as.xts(apply(cumprod(edhec+1),MARGIN=2,ROC,n=12,type="discrete"),order.by=index(edhec))
plot.xts(roc,layout.screens = 1:NCOL(roc),
ylim=c(origin,horizonscale),
panel=horizon.panel,bty="n",
auto.grid=FALSE,
yax.loc="none",
yaxt="n", #, xaxt="n")
main="Horizon Plot of EDHEC Indexes 12-month Rolling Return")

Thursday, August 16, 2012

plot.xts is wonderful

As mentioned in FOSS Trading post A New plot.xts yesterday

“The Google Summer of Code (2012) project to extend xts has produced a very promising new plot.xts function. Michael Weylandt, the project's student, wrote R-SIG-Finance to request impressions, feedback, and bug reports. The function is housed in the xtsExtra package of the xts project on R-Forge.

Please try xtsExtra::plot.xts and let us know what you think. A sample of the eye-candy produced by the code in Michael's email is below. Granted, this isn't a one-liner, but it's certainly impressive! Great work Michael!”

and as announced by author Michael Weylandt in https://stat.ethz.ch/pipermail/r-sig-finance/2012q3/010652.html

“As the community which makes the most heavy use of xts, I would like to draw your attention to a new set of plotting functions for xts objects available as part of Google Summer of Code 2012. This work represents a major overhaul of previously existing plot.xts and should provide you with the most comprehensive and flexible time series plotting available in R. Features include:

  • "automagic" layout construction and axis alignment
  • smart argument recycling
  • panel function abilities
  • more attractive candle and bar plots for OHLC objects
  • scatterplots to view the co-evolution of multiple series
  • event markers
  • regime highlighting
  • time-oriented barplots via barplot.xts [based on code by Peter Carl]
  • interoperability with all known R time series classes using the xts try/reclass paradigm

while retaining the same smart axis formatting and gridlines that plot.xts provided. We have made every effort to maintain complete compatibility with documented usages of the old plot.xts and to be 95% compatible with plot.zoo. My goal has been to craft a design which uses smart defaults to put attractive and informative graphics ever at your fingertips, while remaining flexible enough for "power-users" to craft every detail as they desire…”

Michael has done a fantastic job with this Google Summer of Code (GSOC) project, and I look forward to incorporating all the new features of plot.xts.  As a quick example very different from that shown in the post and mailing list announcement, I thought it would be fun to show how we can use plot.xts to replace the veteran chart.TimeSeries and charts.PerformanceSummary plots from the PerformanceAnalytics package.

For the chart.TimeSeries, I will almost exactly replicate the chart given in the documentation.

From TimelyPortfolio

Not really as a useful application but for the sake of demonstration, we can produce something like this to separate styles and also avoid event label collision (just implemented by Michael very late last night).

From TimelyPortfolio

plot.xts even allows plot.zoo and lattice panel type functionality which allows us to do very nice things like in this charts.PerformanceSummary style chart.

From TimelyPortfolio

I very much appreciate this very powerful contribution to R.  Thanks to everyone involved.

R code in GIST (note: install xtsExtra from R-forge) :

#install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
require(PerformanceAnalytics)
require(xtsExtra) #if you get error, please install xtsExtra from r-forge as shown in the top line
require(RColorBrewer)
#function add alpha or transparency to colors
addalpha <- function(cols,alpha=180) {
rgbcomp <- col2rgb(cols)
rgbcomp[4] <- alpha
return(rgb(rgbcomp[1],rgbcomp[2],rgbcomp[3],rgbcomp[4],maxColorValue=255))
}
#replicate example provided in chart.TimeSeries with plot.xts
# These are start and end dates, formatted as xts ranges.
## http://www.nber.org-cycles.html
cycles.begin.dates<-c("1857-06",
"1860-10",
"1865-04",
"1869-06",
"1873-10",
"1882-03",
"1887-03",
"1890-07",
"1893-01",
"1895-12",
"1899-06",
"1902-09",
"1907-05",
"1910-01",
"1913-01",
"1918-08",
"1920-01",
"1923-05",
"1926-10",
"1929-08",
"1937-05",
"1945-02",
"1948-11",
"1953-07",
"1957-08",
"1960-04",
"1969-12",
"1973-11",
"1980-01",
"1981-07",
"1990-07",
"2001-03",
"2007-12"
)
cycles.end.dates<-c("1858-12",
"1861-06",
"1867-12",
"1870-12",
"1879-03",
"1885-05",
"1888-04",
"1891-05",
"1894-06",
"1897-06",
"1900-12",
"1904-08",
"1908-06",
"1912-01",
"1914-12",
"1919-03",
"1921-07",
"1924-07",
"1927-11",
"1933-03",
"1938-06",
"1945-10",
"1949-10",
"1954-05",
"1958-04",
"1961-02",
"1970-11",
"1975-03",
"1980-07",
"1982-11",
"1991-03",
"2001-11",
"2009-06"
)
# Event lists - FOR BEST RESULTS, KEEP THESE DATES IN ORDER
risk.dates = c(
"1987-10-19",
"1994-02-01",
"1997-07-01",
"1998-08-17",
"1998-09-23",
"2000-07-01",
"2011-09-11")
risk.labels = c(
"Black Monday",
"Bond Crash",
"Asian Crisis",
"Russian Crisis",
"LTCM",
"Tech Bubble",
"Sept 11")
data(edhec)
R=edhec[,1:3,drop=FALSE]
Return.cumulative = cumprod(1+R) - 1
plot.xts(Return.cumulative, main="Default plot.xts Chart") #using all the defaults
#png("chartTimeSeries.png",width=640,height=567,units="px")
plot.xts(Return.cumulative, screens=1, #screens=1 probably most appropriate for this application
blocks = list( #get blocks to plot from above set dates in list form
start.time=paste(cycles.begin.dates,"-01",sep=""),
end.time=paste(cycles.end.dates,"-01",sep=""),
col = "lightblue"),
events = list( #get events to plot from above risk.date in list form
time = risk.dates,
label = risk.labels,
col = "red"),
lwd = 2,
legend.loc = "bottomright", auto.legend=TRUE,
main="EDHEC Style Indexes")
#dev.off()
#png("chartTimeSeries with extra.png",width=640,height=567,units="px")
plot.xts(Return.cumulative,
screens=c(1,2,2), #plot 1st series in 1st panel and 2nd and 3rd series in 2nd panel
layout.screens=c(1,2,2), #just as an example change the layout so 1st panel is 1/3 of area and 2nd is bottom 2/3
blocks = list( #set up blocks for recessions
start.time=paste(cycles.begin.dates,"-01",sep=""),
end.time=paste(cycles.end.dates,"-01",sep=""),
col = "lightblue"),
events = list( #add some event lines
time = risk.dates,
label = risk.labels,
offset = c(0.5,0.5,0.5,0.5,-0.5,0.5), #collision control can also use yadj as in next line
y.adj = c(0,0.5), #with recycling will alternate position for each label
col = "purple"), #don't know why you would use purple but in case you do
lwd = 2,
legend.loc = "bottomright", auto.legend=TRUE,
main="EDHEC Style Indexes")
#dev.off()
#for a little more advanced application we can start to do charts.PerformanceSummary style plot.xts
cumulreturn.panel <- function(...) {
mtext("Cumulative Return", side=1, adj=1, line=-3)
default.panel(...)
abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lty=3)
abline(h=0, col="black")
}
es.panel <- function(index,x,...) {
mtext("Expected Shortfall", side=1, adj=1, line=-3)
default.panel(index,x,...)
#silly to do this but if we wanted just certain points like every 4 months we could do something like this
#default.panel(index[seq(1,NROW(index),by=4)],coredata(x[seq(1,NROW(index),by=4)]),...)
#abline(h=0, col="black")
abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lty=3)
abline(h=par("yaxp")[1], col="black")
}
drawdown.panel <- function(index,x,...) {
mtext("Drawdown", side=1, adj=1, line=-2)
default.panel(index,x,...)
#silly to do this but if we wanted just certain points like every 4 months we could do something like this
#default.panel(index[seq(1,NROW(index),by=4)],coredata(x[seq(1,NROW(index),by=4)]),...)
#abline(h=0, col="black")
abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lty=3)
abline(h=par("usr")[3], col="black")
}
#get some risk measurements to add for Performance Summary style plot
Risk.drawdown <- Drawdowns(R)
Risk.es <- rollapplyr(R,FUN="ES",width=36,p=0.95,na.pad=TRUE)
#take care of NA with 0 at beginning and interpolation at end
Risk.es <- apply(Risk.es,MARGIN=2,FUN=na.fill,fill=c(0,"extend"))
#something wrong with returned value from apply.rolling so indexes don't work properly
data.to.plot <- as.xts(cbind(coredata(Return.cumulative),Risk.es,coredata(Risk.drawdown)),order.by=index(edhec))
#png("chartsPerformanceSummary.png",width=640,height=600,units="px")
plot.xts(data.to.plot,
lwd = c(2,1,1), #do this to show how arguments are recycled
col = brewer.pal(n=9,"PuBu")[c(8,4,6)],
auto.grid = FALSE, #usually auto.grid works just fine but turn off for example purposes
las = 1,yax.loc = "right", # yax.loc could also be flip or left in this case
screens = c(1,1,1,2,2,2,3,3,3), #3 series for each so first 3 in panel 1 second 3 in panel 2 and last 3 in panel 3
layout.screens = c(1,1,2,3), #panel 1 take up first 2 of 4 so 50% and panels 2 and 3 each 25%
bty = "n",
panel = c(cumulreturn.panel,es.panel,drawdown.panel), #c(first.panel,"auto"), #panel cycles through by panel rather than each series
ylab = NA, major.format = "%b %Y", minor.ticks = FALSE,
legend.loc = c("topleft",NA,NA), auto.legend = TRUE,
legend.pars = list(bty = "n", horiz=TRUE), #make legend box transparent
cex.axis = 0.9,
main = NA)
title(main = "Performance Summary of EDHEC Indexes", adj = 0, outer = TRUE, line = -1.5)
#dev.off()
#just playing around with type="h"
#think will only work if all colors very different
plot.xts(Risk.drawdown["2000::2001",],type="h",screens=1,lwd=c(10,25,15),
col=apply(as.matrix(brewer.pal(n=9,"PuBu")[c(8,4,6)]),MARGIN=1,FUN=addalpha,alpha=120))

Friday, August 10, 2012

Animated GIF Annual Correlation of 48 Industries for 50 Years

Inspired by http://blogs.reuters.com/felix-salmon/2012/08/06/chart-of-the-day-hft-edition/, I thought I would build on 48 Industries (Dendrogram Ordered) Over 50 Years, “Trend is Not Your Friend” Applied to 48 Industries, and 48 Industries Since 1963 with an animated GIF of annual correlation for the Kenneth French 48 Industry Data set.

From TimelyPortfolio

R code in GIST (do raw for copy/paste):

#get very helpful Ken French data
#for this project we will look at Industry Portfolios
#http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip
require(latticeExtra)
require(animation)
require(PerformanceAnalytics)
require(quantmod)
#my.url will be the location of the zip file with the data
my.url="http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip"
#this will be the temp file set up for the zip file
my.tempfile<-paste(tempdir(),"\\frenchindustry.zip",sep="")
#my.usefile is the name of the txt file with the data
my.usefile<-paste(tempdir(),"\\48_Industry_Portfolios_daily.txt",sep="")
download.file(my.url, my.tempfile, method="auto",
quiet = FALSE, mode = "wb",cacheOK = TRUE)
unzip(my.tempfile,exdir=tempdir(),junkpath=TRUE)
#read space delimited text file extracted from zip
french_industry <- read.table(file=my.usefile,
header = TRUE, sep = "",
as.is = TRUE,
skip = 9, nrows=12211)
#get dates ready for xts index
datestoformat <- rownames(french_industry)
datestoformat <- paste(substr(datestoformat,1,4),
substr(datestoformat,5,6),substr(datestoformat,7,8),sep="-")
#get xts for analysis
french_industry_xts <- as.xts(french_industry[,1:NCOL(french_industry)],
order.by=as.Date(datestoformat))
#divide by 100 to get percent
french_industry_xts <- french_industry_xts/100
#delete missing data which is denoted by -0.9999
french_industry_xts[which(french_industry_xts < -0.99,arr.ind=TRUE)[,1],
unique(which(french_industry_xts < -0.99,arr.ind=TRUE)[,2])] <- 0
#get a vector of the end of years
evaldates <- endpoints(french_industry_xts,"years")
saveGIF(
for(i in 2:length(evaldates)) {
#do correlation table
ca <- cor(french_industry_xts[evaldates[i-1]:evaldates[i],])
#replace na with 0
ca[which(is.na(ca),arr.ind=TRUE)[,]] <- 0
#get colors to use for heat map
brew <- brewer.pal(name="RdBu",n=5)
#get color ramp
cc.brew <- colorRampPalette(brew)
#apply color ramp
cc <- cc.brew(nrow(ca))
#do heatmap and sort by degree of correlation to VFINX (Vanguard S&P 500)
#heatmap(ca,symm=TRUE,Rowv=NA,Colv=NA,col=cc,RowSideColors=cc,main="")
#title(main=paste("Correlation Table\n",index(french_industry_xts)[evaldates[i]],sep=""),font.main=1,outer=TRUE,line=-2,cex.main=1.3)
#do with dendrogram ordering
heatmap(ca,symm=TRUE,col=cc,RowSideColors=cc,main="")
title(main=paste("Correlation Table (Dendrogram Ordered)\n",index(french_industry_xts)[evaldates[i]],sep=""),font.main=1,outer=TRUE,line=-3,cex.main=1.3,adj=0)
}
)

Thursday, August 9, 2012

48 Industries (Dendrogram Ordered) Over 50 Years

Thanks to reader AHWest for the comment on post 48 Industries Since 1963.

“I think it would be interesting to see the industries ordered by some sort of similarity of returns.”

I think this is a great suggestion, and I would like to see it also.  I tried the dendrogram plot technique from Inspirational Stack Overflow Dendrogram Applied to Currencies, but then I spotted the dendrogramGrob in the latticeExtra documentation.  This was much easier, and in a couple of lines, we are able to order and connect the 48 industries.

From TimelyPortfolio
R code from GIST (do raw for copy/paste):
require(fAssets)
require(latticeExtra)
require(quantmod)
require(PerformanceAnalytics)
#my.url will be the location of the zip file with the data
my.url="http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip"
#this will be the temp file set up for the zip file
my.tempfile<-paste(tempdir(),"\\frenchindustry.zip",sep="")
#my.usefile is the name of the txt file with the data
my.usefile<-paste(tempdir(),"\\48_Industry_Portfolios_daily.txt",sep="")
download.file(my.url, my.tempfile, method="auto",
quiet = FALSE, mode = "wb",cacheOK = TRUE)
unzip(my.tempfile,exdir=tempdir(),junkpath=TRUE)
#read space delimited text file extracted from zip
french_industry <- read.table(file=my.usefile,
header = TRUE, sep = "",
as.is = TRUE,
skip = 9, nrows=12211)
#get dates ready for xts index
datestoformat <- rownames(french_industry)
datestoformat <- paste(substr(datestoformat,1,4),
substr(datestoformat,5,6),substr(datestoformat,7,8),sep="-")
#get xts for analysis
french_industry_xts <- as.xts(french_industry[,1:NCOL(french_industry)],
order.by=as.Date(datestoformat))
#divide by 100 to get percent
french_industry_xts <- french_industry_xts/100
#delete missing data which is denoted by -0.9999
french_industry_xts[which(french_industry_xts < -0.99,arr.ind=TRUE)[,1],
unique(which(french_industry_xts < -0.99,arr.ind=TRUE)[,2])] <- 0
#get price series or cumulative growth of 1
french_industry_price <- cumprod(french_industry_xts+1)
#get 250 day rate of change or feel free to change to something other than 250
roc <- french_industry_price
#split into groups so do not run out of memory
for (i in seq(12,48,by=12)) {
roc[,((i-11):(i))] <- ROC(french_industry_price[,((i-11):(i))],n=250,type="discrete")
}
roc[1:250,] <- 0
# try to do http://stackoverflow.com/questions/9747426/how-can-i-produce-plots-like-this
# was much easier to use latticeExtra
# get dendrogram data from hclust
# backward and repetitive but it works
t <- assetsDendrogramPlot(as.timeSeries(french_industry_xts))
# thanks to the latticeExtra example
dd.row <- as.dendrogram(t$hclust)
row.ord <- order.dendrogram(dd.row)
xyplot(roc[,row.ord],
layout=c(1,48), ylim=c(0,0.25),
scales = list(tck = c(1,0), y = list(draw = FALSE,relation = "same")),
horizonscale=0.25,
origin = 0,
colorkey = TRUE,
#since so many industries, we will comment out grid
panel = function(x,y,...) {
panel.horizonplot(x,y,...) #feel free to change to whatever you would like)
# panel.grid(h=3, v=0,col = "white", lwd=1,lty = 3)
},
ylab = list(rev(colnames(roc[,row.ord])), rot = 0, cex = 0.7, pos = 3),
xlab = NULL,
par.settings=theEconomist.theme(box = "gray70"),
#use ylab above for labelling so we can specify FALSE for strip and strip.left
strip = FALSE,
strip.left = FALSE,
main = "French Daily 48 Industry (Dendrogram Ordered) 1963-2011\n source: http://mba.tuck.dartmouth.edu/pages/faculty/ken.french",
legend =
list(right =
list(fun = dendrogramGrob,
args =
list(x = dd.row, ord = row.ord,
side = "right",
size = 10))))

Wednesday, August 8, 2012

“Trend is Not Your Friend” Applied to 48 Industries

Please see previous post Crazy RUT in Academic Context Why Trend is Not Your Friend.

I’ll repeat the intro to the post mentioned above, so we can all get caught back up.

In response to Where are the Fat Tails?, reader vonjd very helpfully referred me to this paper The Trend is Not Your Friend! Why Empirical Timing Success is Determined by the Underlying’s Price Characteristics and Market Efficiency is Irrelevant by Peter Scholz and Ursula Walther. The authors conclude

“Our study on the basis of real data clearly confirms the hypothesis that the asset price characteristics of the underlying price process have a crucial impact on timing results. This allows us to forecast the timing success depending on the market's parameters. An OLS
regression analysis supports our predictions and verifies our assumption that the drift has the
strongest influence on timing success. By contrast, the higher moments (skewness, kurtosis)
seem not to have any significant impact on the timing result in the empirical sample. As we
presumed, the level of market development, and hence the degree of efficiency, does not play
any role. Trading worked coincidentally rather well in the developed world and quite poorly in
the emerging markets. The driving factor for the timing success is the parametric environment the trading system stumbles on…

Our study contributes to the discussion by providing a structured analysis of the relevance of the most important price process parameters. As a result, the traditional explanations for timing success can be abandoned: we find that it is very likely for the SMA trading rule to generate excess returns over its benchmark if the underlying price path exhibits negative drifts, high serial autocorrelation, low volatilities of returns, and highly clustered volatilities. Drift and autocorrelation of the underlying asset seem to have the largest impact, though.”

One of my initial ideas for extending the research was to incorporate a much larger set of indexes over a longer period of time.  As I was working on 48 Industries Since 1963, I decided 50 years of data on 48 different indexes would be a great dataset to apply the ideas and methods presented in the paper.

Using R and all its wonderful packages, it is surprisingly easy to accomplish.  Let’s see if we can test “drift and autocorrelation…have the largest impact” on excess returns with industries also.  I’ll try not to get too statistical.

In terms of drift or annualized return, we can see a linear inverse relationship between return and out(under)performance of the 200 day moving average system, so the better the performance of the industry, the less likely a moving average system is to outperform.

From TimelyPortfolio

In terms of the GARCH model effects on excess returns, a parallel coordinate chart will best start our exploration.  Lines are colored by the excess return of a moving average system on each industry.

From TimelyPortfolio

Mu and alpha1 seem to most heavily influence the ability of a moving average system to outperform.  Let’s isolate our chart to mu and alpha1 and add ar1 based on the authors’ findings.  I am very tempted to try to explain GARCH here, but for the sake of brevity, I’ll refrain.  This paper and this Portfolio Probe post offer a good introduction to GARCH.

From TimelyPortfolio

Now let’s look at a scatterplot of excess return versus each of the GARCH stats, starting with those with the most influence.

From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio

In conclusion, it seems the same effects observed by the authors also apply to US industry indexes.  In future posts, I’ll add a little more statistical rigor to the analysis and apply to other indexes.

Now, I just cannot resist using a horizon plot to evaluate the rolling 250 day excess returns of a moving average system over buy and hold.  As you can see, a bull market favors buy and hold.  The 70s and 2008-2009 were very kind to a moving average approach.

From TimelyPortfolio

R code in GIST (do raw for copy/paste):

#get very helpful Ken French data
#for this project we will look at Industry Portfolios
#http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip
require(latticeExtra)
require(PerformanceAnalytics)
require(quantmod)
#my.url will be the location of the zip file with the data
my.url="http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip"
#this will be the temp file set up for the zip file
my.tempfile<-paste(tempdir(),"\\frenchindustry.zip",sep="")
#my.usefile is the name of the txt file with the data
my.usefile<-paste(tempdir(),"\\48_Industry_Portfolios_daily.txt",sep="")
download.file(my.url, my.tempfile, method="auto",
quiet = FALSE, mode = "wb",cacheOK = TRUE)
unzip(my.tempfile,exdir=tempdir(),junkpath=TRUE)
#read space delimited text file extracted from zip
french_industry <- read.table(file=my.usefile,
header = TRUE, sep = "",
as.is = TRUE,
skip = 9, nrows=12211)
#get dates ready for xts index
datestoformat <- rownames(french_industry)
datestoformat <- paste(substr(datestoformat,1,4),
substr(datestoformat,5,6),substr(datestoformat,7,8),sep="-")
#get xts for analysis
french_industry_xts <- as.xts(french_industry[,1:NCOL(french_industry)],
order.by=as.Date(datestoformat))
#divide by 100 to get percent
french_industry_xts <- french_industry_xts/100
#delete missing data which is denoted by -0.9999
french_industry_xts[which(french_industry_xts < -0.99,arr.ind=TRUE)[,1],
unique(which(french_industry_xts < -0.99,arr.ind=TRUE)[,2])] <- 0
#get price series or cumulative growth of 1
french_industry_price <- cumprod(french_industry_xts+1)
#do 200 day moving average for initial testing
#just change n to test on other widths or rolling period
ma <- as.xts(apply(french_industry_price, MARGIN = 2, FUN = runMean, n=200), order.by = index(french_industry_price))
#do plot to test reasonableness
chart.TimeSeries(ma,ylog=TRUE)
#set up system to enter when price moves above moving average
#exit when below
ma.system <- lag(as.xts(
apply(french_industry_price > ma, MARGIN = 2, as.numeric),
order.by = index(french_industry_price)),
k=1) * french_industry_xts
#get returns cumulative and annualized for the entire period
ret.comp.cumul <- Return.cumulative(ma.system) - Return.cumulative(french_industry_xts)
ret.bh.ann <- Return.annualized(french_industry_xts)
ret.comp.ann <- Return.annualized(ma.system) - ret.bh.ann
rownames(ret.comp.ann) <- "Out(under)performance"
#get colors to use for heat map style coloring by out/under performance
brew <- brewer.pal(name="RdBu",n=5)
#get color ramp
cc.brew <- colorRampPalette(brew)
#apply color ramp
cc <- cc.brew(length(ret.comp.ann))
#do colors based on out/under performance but with gray so visible when labelling
cc.palette <- colorRampPalette(c(cc[1],"gray60",cc[length(cc)]))
cc.levpalette <- cc.palette(length(ret.comp.ann))
cc.levels <- level.colors(ret.comp.ann, at = do.breaks(c(-max(abs(ret.comp.ann)),max(abs(ret.comp.ann))),length(ret.comp.ann)),
col.regions = cc.levpalette)
#plot(ret.comp.cumul ~ Return.cumulative(french_industry_xts),pch=19)
#plot out/underperformance of ma system vs annualized compounded return of buyhold
plot(x=ret.bh.ann,y=ret.comp.ann,
pch=19, col=cc.levels,bty="l",
xlab = "Ann. Return of BuyHold", ylab = "Out(under)Performance of MA")
text(x=ret.bh.ann,y=ret.comp.ann, labels=colnames(ret.comp.ann), cex=0.7, pos=2, col=cc.levels)
#plot out/underperformance of ma system vs std. dev. of buyhold
plot(y=ret.comp.ann, x=apply(na.omit(french_industry_xts),MARGIN=2,sd),
pch=19,col=cc.levels,bty="l",
xlab = "StdDev of BuyHold", ylab = "Out(under)Performance of MA")
text(y=ret.comp.ann, x=apply(na.omit(french_industry_xts),MARGIN=2,sd), labels=colnames(ret.comp.ann), cex=0.7, pos=2, col=cc.levels)
#plot out/underperformance of ma system vs Omega of buyhold
plot(y=ret.comp.ann, x=Omega(french_industry_xts),
pch=19,col=cc.levels,bty="l",
xlab = "Omega of BuyHold", ylab = "Out(under)Performance of MA")
text(y=ret.comp.ann, x=Omega(french_industry_xts), labels=colnames(ret.comp.ann), cex=0.7, pos=2, col=cc.levels)
#using rugarch get garch stats similar to those explored in the research
require(rugarch)
spec = ugarchspec(
variance.model=list(garchOrder=c(1,1)),
mean.model=list(armaOrder=c(1,1), include.mean=T))
#set up function to get garch stats through apply function
gfNa <- function(data, spec) {
x <- na.omit(coredata(data))
gf <- suppressWarnings(ugarchfit(spec=spec, data=x))
stats <- coef(gf)
return(stats)
}
#do apply to get garch stats across all industries
gf.stats <- apply(french_industry_xts[,1:NCOL(french_industry_xts)],MARGIN=2,FUN=gfNa,spec=spec)
#plot return comparison by each garch stat
for (i in 1:NROW(gf.stats)) {
plot(y=ret.comp.ann, x=gf.stats[i,],
pch=19, col=cc.levels, bty="l",
xlab=rownames(gf.stats)[i], ylab="Out(under)Performance of MA")
text(y=ret.comp.ann, x=gf.stats[i,], labels=colnames(ret.comp.ann), cex=0.7, pos=2, col=cc.levels)
}
#thinking through linear model
#linmod = lm(as.vector(ret.comp.ann)~as.vector(gf.stats[1,]))
#plot(linmod)
#do parallel coordinates chart with color
require(MASS)
parcoord(cbind(t(ret.comp.ann),t(gf.stats)[,c(1,2,5)]),
col = cc.levels,lwd = 2,
main = "Out(under)Performance by GARCH Stat")
parcoord(cbind(t(ret.comp.ann),t(gf.stats)),
col = cc.levels,lwd = 2,
main = "Out(under)Performance by GARCH Stat")
#get rolling out(under) performance for horizon chart
#change na to 0 in ma.system returns
ma.system[which(is.na(ma.system),arr.ind=TRUE)[,1],
unique(which(is.na(ma.system),arr.ind=TRUE)[,2])] <- 0
ma_system_price <- cumprod(1+ma.system)
roc <- french_industry_price
#split into groups so do not run out of memory
for (i in seq(12,48,by=12)) {
#get difference in rolling performance
roc[,((i-11):(i))] <- ROC(ma_system_price[,((i-11):(i))],n=250,type="discrete") -
ROC(french_industry_price[,((i-11):(i))],n=250,type="discrete")
}
roc[1:250,] <- 0
#do a horizon plot of all 48 industries with horizonscale of 0.25
horizonplot(roc,
layout=c(1,48),
horizonscale=0.25, #feel free to change to whatever you would like
scales = list(tck = c(1,0), y = list(draw = FALSE,relation = "same")),
origin = 0,
colorkey = FALSE,
#since so many industries, we will comment out grid
# panel = function(x, ...) {
# panel.horizonplot(x, ...)
# panel.grid(h=3, v=0,col = "white", lwd=1,lty = 3)
# },
ylab = list(rev(colnames(roc)), rot = 0, cex = 0.7, pos = 3),
xlab = NULL,
par.settings=theEconomist.theme(box = "gray70"),
#use ylab above for labelling so we can specify FALSE for strip and strip.left
strip = FALSE,
strip.left = FALSE,
main = "Moving Averages System Performance on French Daily 48 Industry 1963-2011\n source: http://mba.tuck.dartmouth.edu/pages/faculty/ken.french")

Tuesday, August 7, 2012

48 Industries Since 1963

Please see http://timelyportfolio.blogspot.com/search/label/horizonplot for all horizon plot posts.

Once more thanks to Ken French for his data, we can accomplish something I think is fairly amazing.  In 640x800, we can see 250 day rollling returns for 48 U.S. industries since 1963.

From TimelyPortfolio

R code in GIST (do raw for copy/paste):

#get very helpful Ken French data
#for this project we will look at Industry Portfolios
#http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip
require(latticeExtra)
require(PerformanceAnalytics)
require(quantmod)
#my.url will be the location of the zip file with the data
my.url="http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/48_Industry_Portfolios_daily.zip"
#this will be the temp file set up for the zip file
my.tempfile<-paste(tempdir(),"\\frenchindustry.zip",sep="")
#my.usefile is the name of the txt file with the data
my.usefile<-paste(tempdir(),"\\48_Industry_Portfolios_daily.txt",sep="")
download.file(my.url, my.tempfile, method="auto",
quiet = FALSE, mode = "wb",cacheOK = TRUE)
unzip(my.tempfile,exdir=tempdir(),junkpath=TRUE)
#read space delimited text file extracted from zip
french_industry <- read.table(file=my.usefile,
header = TRUE, sep = "",
as.is = TRUE,
skip = 9, nrows=12211)
#get dates ready for xts index
datestoformat <- rownames(french_industry)
datestoformat <- paste(substr(datestoformat,1,4),
substr(datestoformat,5,6),substr(datestoformat,7,8),sep="-")
#get xts for analysis
french_industry_xts <- as.xts(french_industry[,1:NCOL(french_industry)],
order.by=as.Date(datestoformat))
#divide by 100 to get percent
french_industry_xts <- french_industry_xts/100
#delete missing data which is denoted by -0.9999
french_industry_xts[which(french_industry_xts < -0.99,arr.ind=TRUE)[,1],
unique(which(french_industry_xts < -0.99,arr.ind=TRUE)[,2])] <- 0
#get price series or cumulative growth of 1
french_industry_price <- cumprod(french_industry_xts+1)
#get 250 day rate of change or feel free to change to something other than 250
roc <- french_industry_price
#split into groups so do not run out of memory
for (i in seq(12,48,by=12)) {
roc[,((i-11):(i))] <- ROC(french_industry_price[,((i-11):(i))],n=250,type="discrete")
}
roc[1:250,] <- 0
#do a horizon plot of all 48 industries with horizonscale of 0.25
horizonplot(roc,
layout=c(1,48),
horizonscale=0.25, #feel free to change to whatever you would like
scales = list(tck = c(1,0), y = list(draw = FALSE,relation = "same")),
origin = 0,
colorkey = FALSE,
#since so many industries, we will comment out grid
# panel = function(x, ...) {
# panel.horizonplot(x, ...)
# panel.grid(h=3, v=0,col = "white", lwd=1,lty = 3)
# },
ylab = list(rev(colnames(roc)), rot = 0, cex = 0.7, pos = 3),
xlab = NULL,
par.settings=theEconomist.theme(box = "gray70"),
#use ylab above for labelling so we can specify FALSE for strip and strip.left
strip = FALSE,
strip.left = FALSE,
main = "French Daily 48 Industry 1963-2011\n source: http://mba.tuck.dartmouth.edu/pages/faculty/ken.french")

Friday, August 3, 2012

Horizon Plots in Base Graphics

for background please see prior posts More on Horizon Charts, Application of Horizon Plots, Horizon Plot Already Available, and Cubism Horizon Charts in R

There are three primary graphics routes in R (base graphics, lattice, and ggplot2), and each have their zealots.  Last time in More on Horizon Charts, I used lattice and latticeExtra. This time we will build horizon plots in base graphics, and I was pleased with the result. Unfortunately, there is one small issue in that the points of change from positive to negative overlap.  Please let me know if you have a solution. Thanks to helpful readers for their comments and code changes, this is now fixed.

From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio

Now, let’s implement a for loop and mirror the negative values.

From TimelyPortfolio

With horizon chart functionality in base graphics, hopefully we can add this type now to other packages.  Here is a potential example using quantmod.

From TimelyPortfolio

 

R code in GIST (do raw for copy/paste)

require(RColorBrewer)
require(quantmod)
require(PerformanceAnalytics)
data(managers)
#let's do managers from 2002 to 2004 to get positive and negative
x <- cumprod(1+managers["2002::2004"])[,1] - 1
#get some decent colors from RColorBrewer
#we will use colors on the edges so 2:4 for red and 7:9 for blue
col.brew <- brewer.pal(name="RdBu",n=10)
#get this to ease using it later
n<-nrow(x)
#set scale to be 10%
horizonscale=0.1
#remove space around chart
par(mar=c(2,1,1,1))
plot(index(x), coredata(x), type="n", bty="n", las=1, yaxt="n", xlab=NA, ylab=NA, ylim=c(-horizonscale,horizonscale))
#thanks http://stackoverflow.com/questions/9630014/polygon-for-xts-objects
#draw first positive band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) > 0,coredata(x), 0),0),
col=col.brew[7]
)
#draw first negative band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) < 0 ,coredata(x), 0),0),
col=col.brew[4]
)
#overlay second positive band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) > 0.1,coredata(x) - 0.1, 0),0),
col=col.brew[8]
)
#overlay second negative band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) < -0.1 ,coredata(x) + 0.1, 0),0),
col=col.brew[3]
)
#overlay third positive band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) > 0.2 ,coredata(x) - 0.2, 0),0),
col=col.brew[9]
)
#overlay third negative band
polygon(
index(x)[c(1,1:n,n)],
c(0,ifelse(coredata(x) < -0.2 ,coredata(x) + 0.2, 0),0),
col=col.brew[2]
)
#little touch up to get a line at extending left to right
abline(h=0,col="black")
#add a line at the bottom of the chart
abline(h=par("usr")[3],col="black")
#now let's do it with a loop and flip the negative up
nbands = ceiling(max(abs(coredata(x)))/horizonscale)
plot(index(x), abs(coredata(x)), type="n", bty="n", las=1, yaxt="n", xaxt="n", xlab=NA, ylab=NA, ylim=c(0,horizonscale))
#thanks to helpful reader A. Zolot
par(usr=c(index(x)[1],index(x)[n],0,horizonscale),mar=c(2,1,1,1)) # 0-margines
for (i in 1:nbands) {
#draw positive
polygon(
c(index(x)[1], index(x), index(x)[n]),
c(0, coredata(x) - (i-1) * horizonscale,0),
col=col.brew[length(col.brew)-nbands+i-1],
border=NA
)
#draw negative
polygon(
c(index(x)[1], index(x), index(x)[n]),
c(0, -coredata(x) - (i-1) * horizonscale, 0),
col=col.brew[nbands-i+1],
border=NA
)
}
abline(h=0,col="black")
axis.Date(side=1,x=index(x),pos=0)

Thursday, August 2, 2012

How do you say “We Will Do Whatever It Takes” in Thai?

As the market has already started to poke holes in Draghi’s promise, I thought it would be good to continue the series of posts that I began with the British version “We Will Do Whatever it Takes” with my favorite article written during the Asia Pacific crisis in 1997.

The Economist, “The baht spills over”: May 22, 1997

Too soon to crow

“Even if Thailand's troubles are in a regional class of their own, its neighbours still have reason to fear the knock-on effect of a successful onslaught on the baht. For now, Thailand has been crowing about a victory over the speculators. Mr Amnuay has even promised to bring interest rates, which have been kept high to protect the baht, down by two percentage points this year. It is as if a little guy, having survived one round with the big bullies, was dancing around with his fists up, saying, go on, hit me again.

They almost certainly will. The defence of the baht was achieved by a series of extraordinary measures that will be hard to sustain…”

A little graphical perspective might help us remember when this article was written.

image

For the entire set of Economist articles on Thailand during this period, go to http://www.economist.com/topics/thailand-1?page=33&%2524Version=0&%2524Path=%252F&%2524Domain=.economist.com.  Also, please read the incredibly thorough account of the 1997 crisis at The Nation blog.

If throughout financial market history, the central banks and politicians ever defeated the markets (greedy speculators), please let me know.  I am absolutely amazed at the confidence in central banks and politicians, especially in the United States.

 

R code:

#plot Thai baht
require(lattice)
require(latticeExtra)
require(quantmod)
getSymbols("DEXTHUS",src="FRED")
DEXTHUS <- 1/DEXTHUS["1995-01-01::1999-12-31",]
asTheEconomist(
  xyplot(DEXTHUS,
       panel = function (x,y,...) {
         panel.xyplot(x,y,...)
         panel.abline(v=index(DEXTHUS)[which(index(DEXTHUS)=="1997-05-22")],col="black", lty=3)
         panel.text(x=index(DEXTHUS)[which(index(DEXTHUS)=="1997-05-22")], y = 0.025, labels="Economist article The Baht Spills Over", col="indianred4", srt=90, pos=3)
       },
        main = "Thai Baht Before and After Crisis (1995-1999) \n source : Federal Reserve Bank of St. Louis")
)

Wednesday, August 1, 2012

More on Horizon Charts

for background please see prior posts Application of Horizon Plots, Horizon Plot Already Available, and Cubism Horizon Charts in R

Some feedback has led me to think that I might have been a little ambitious with my last post on horizon charts. I thought it might be helpful to quickly provide an overview on horizon charts. Since I am not a visualization expert, I will try to compile the best examples and tutorials that I have found. At the end, I will do a step by step walkthrough of the construction of a horizon chart in R.

Below is a great talk in its entirety by Mike Bostock.  For the discussion on horizon charts, skip to the 11:20 mark.

Mike Bostock @ Square talks about Time Series Visualization from Librato on Vimeo.

Fortunately, Mike also provides the interactive example from the video demonstrating the construction of a horizon chart from an area chart.  I have embedded it below (see https://bl.ocks.org/1483226 for the full example from the original source).

This quick YouTube clip from Panopticon all the way back in 2009 also does a very nice job explaining construction.  If you learn better from reading, a short paper from Hannes Reijner of Panopticon Software covers the same material.

Jeffrey Heer of Stanford University with Nicholas Kong and Maneesh Agrawala from University of California, Berkeley, investigate the effectiveness of horizon plots and make some recommendations for their use.  The chart below also offers another explanation of the steps in building a horizon plot.

image

They recommend and conclude

Layered Bands Are Beneficial As Chart Size Decreases
We found that dividing a chart into layered bands reliably increased estimation time and increased estimation error at constant chart heights. However, we also found that 2-band mirrored charts led to better estimation accuracies for chart heights less than 24 pixels (6.8 mm on our displays). For larger chart sizes, we advise scaling 1-band mirrored charts.  For smaller sizes, we advise adding layered bands.

Extending the horizon into healthcare, here is an interesting project using horizon graphs for visualization of diabetes care.

HorizonVis

Interactive Visual Exploration of Multivariate Medical Measurements in Diabetes Care


HorizonVis

Lead / Contact
Wolfgang Aigner

Team
Wolfgang Aigner, Vienna University of Technology
Michael Atanasov, HTBL Krems
Alexander Rind, Vienna University of Technology
Philipp Schindler, HTBL Krems
Reinhardt Wenzina, HTBL Krems

Partners
Vienna University of Technology, Institute of Software Technology & Interactive Systems
HTBL Krems, Department of Information Technology

 

Now For My Own Attempt at Explaining

If we use yesterday’s example of a 200 day moving average system where you enter when above the moving average and exit when below, we might like to see a standard time series plot like this one.

From TimelyPortfolio

In a simple world with only one asset or stock, this might be sufficient.  However, we probably will have multiple instruments that we would like to monitor, and dedicating this much height per instrument will require lots of space.  Ideally, we could condense each of these plots, so that we could see many at a time.

In the first step toward condensing, we could extract just the information that is most meaningful, which I consider to be the percent above or below the moving average.

From TimelyPortfolio

We might then try an area chart to better depict above or below 0.

From TimelyPortfolio

I’m sure you are wondering though when we will start reducing height and saving space.  We could start by changing all the negative values to positive values, and add color to represent positive or negative.  This means we cut our chart height by about 1/2.

From TimelyPortfolio

However, we need much more efficiency, so let’s separate the chart into bands.

From TimelyPortfolio

Is there any potential way of separating each of these bands and then recombining them to save space?  Let’s look at each band separately.

Band 1 (0 to 10%)

From TimelyPortfolio

Band 2 (10% to 20%)

From TimelyPortfolio

Band 3 (20% to 30%)

From TimelyPortfolio

Band 4 (30% to 40%) exceeding 3 bands is not recommended

From TimelyPortfolio

Notice that similar to a heat map, the band colors increase in intensity and opaqueness as their values increase. If we layer each band on top of each other, we get a horizon plot, and we can reduce the original chart by 1/6 or even more without significant loss of information.

From TimelyPortfolio

I hope this helps explain why we might use horizon plots and how to make them.  If nothing else, maybe you will have learned some lattice techniques in R.

 

R code in GIST (do raw for copy/paste):

#look at steps in constructing a horizon plot version
#of http://www.mebanefaber.com/timing-model/
#do horizon of percent above or below 10 month / 200 day moving average
require(lattice)
require(latticeExtra)
require(quantmod)
#since we are focused on the horizon plot, let's just look at one stock
tckrs <- "VTI"
getSymbols(tckrs, from = "2006-12-31")
#do horizon of percent above or below 10 month or 200 day moving average
prices <- get(tckrs[1])[,6]
#remove comments below if you would like to look at more than one symbol
#for (i in 2:length(tckrs)) {
# prices <- merge(prices,get(tckrs[i])[,4])
#}
colnames(prices) <- tckrs
#set n to desired moving average width; we'll do 200
n=200
ma <- runMean(prices, n = n)
colnames(ma) <- paste(tckrs, ".MovAvg", sep = "")
xyplot(merge(prices,ma),
col = c("black", "red"),
lty = c(1,3),
screens = 1,
scales = list(tck = c(1,0)),
xlab = NULL,
main = "VTI and 200-day Moving Average")
#but for timing system more interested in whether above or below
#get percent above or below
#we'll leave code to expand beyond one symbol
pctdiff <- (prices / apply(prices, MARGIN = 2, FUN = runMean, n = n) - 1)[n:NROW(prices),]
xyplot(pctdiff,
col.line = "steelblue4",
scales = list(tck = c(1,0)),
xlab = NULL)
xyplot(pctdiff,
border = NA,
col.line = "steelblue4",
scales = list(tck = c(1,0)),
xlab = NULL,
panel = function (...) {
panel.xyarea(origin=0, ...)
#draw horizontal lines at 10% to show where we will place bands
panel.abline(h = seq(-0.4, 0.4, 0.10), col = "white", lwd = 2)
#add black 0 axis back
panel.abline(h = 0, col = "black")
})
#takes a lot of height to represent so let's mirror the negative
#so we can cut height by 1/2
xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit to max of absolute value since we will mirror the negative
ylim = c(0,ceiling(max(abs(coredata(pctdiff)))*10)/10),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0, y ,0), col.line = "steelblue4", origin=0, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < 0, abs(y) ,0), col.line = "indianred3", origin=0, ...)
#draw horizontal lines at 10% to show where we will place bands
panel.abline(h = seq(-0.4, 0.4, 0.10), col = "white", lwd = 2)
#add black 0 axis back
panel.abline(h = 0, col = "black")
})
#do same chart as above but draw box around bands
xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit to max of absolute value since we will mirror the negative
ylim = c(0,ceiling(max(abs(coredata(pctdiff)))*10)/10),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0, y ,0), col.line = "steelblue4", origin=0, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < 0, abs(y) ,0), col.line = "indianred3", origin=0, ...)
#draw horizontal lines at 10% to show where we will place bands
panel.abline(h = seq(-0.4, 0.4, 0.10), col = "white", lwd = 2)
#add black 0 axis back
panel.abline(h = 0, col = "black")
panel.xblocks(x,abs(y)>0,height=0.128,col="white",border="black",alpha=0.3)
panel.text(x=x[1],y=0.11,labels="band1", pos=4)
panel.xblocks(x,abs(y)>0,height=0.228,col="white",border="black",alpha=0.25)
panel.text(x=x[1],y=0.21,labels="band2", pos=4)
panel.xblocks(x,abs(y)>0,height=0.328,col="white",border="black",alpha=0.2)
panel.text(x=x[1],y=0.31,labels="band3", pos=4)
panel.xblocks(x,abs(y)>0,height=0.428,col="white",border="black",alpha=0.15)
panel.text(x=x[1],y=0.41,labels="band4", pos=4)
})
#get 4 reds and 4 blues (one for each band)
reds <- brewer.pal("Reds", n=8)[4:8]
blues <- brewer.pal("Blues", n=8)[4:8]
#so let's start banding to use even less height
#for band1 so we will only show graph from 0 to 0.1
band1 <- xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
lattice.options = list(axis.padding = list(numeric = 0)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
x = list(col.line="black"),
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit y to band height; in this case 10%
ylim = c(0,0.1),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0, y ,0), col.line = blues[4], origin=0, alpha = 0.3, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < 0, abs(y) ,0), col.line = reds[4], origin=0, alpha = 0.3, ...)
#add black 0 axis back
panel.abline(h = 0, col = "black")
})
print(band1)
#we are missing all the values > 0.10
#we will draw band2 0.1 to 0.2 with a darker color
band2 <- xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
lattice.options = list(axis.padding = list(numeric = 0)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
x = list(draw = FALSE),
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit y to band height; in this case 10%
ylim = c(0, 0.1),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0.1, y - 0.1, 0), col.line = blues[4], origin=0, alpha = 0.5, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < -0.1, abs(y) - 0.1, 0), col.line = reds[4], origin=0, alpha = 0.5, ...)
})
print(band2)
#now we are missing all the values > 0.10 + 0.10
#we will draw band3 0.2 to 0.3 with a darker color
band3 <- xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
lattice.options = list(axis.padding = list(numeric = 0)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
x = list(draw = FALSE),
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit y to band height; in this case 10%
ylim = c(0, 0.1),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0.2, y - 0.2, 0), col.line = blues[4], origin=0, alpha = 0.7, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < -0.2, abs(y) - 0.2, 0), col.line = reds[4], origin=0, alpha = 0.7, ...)
})
print(band3)
#going to four bands is not recommended but for this example we will
band4 <- xyplot(pctdiff,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
lattice.options = list(axis.padding = list(numeric = 0)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
x = list(draw = FALSE),
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL,
#limit y to band height; in this case 10%
ylim = c(0, 0.1),
panel = function (x, y, ...) {
#do the positive values in blue
panel.xyarea(x, ifelse(y > 0.3, y - 0.3, 0), col.line = blues[4], origin=0, alpha = 1, ...)
#do the positive values in blue
panel.xyarea(x, ifelse(y < -0.3, abs(y) - 0.3, 0), col.line = reds[4], origin=0, alpha = 1, ...)
})
print(band4)
#combine all four bands/layers to one horizonplot
print(band1+band2+band3+band4)
#of course using the horizonplot function from latticeExtra
#makes this much easier
horizonplot(pctdiff,
strip.left=FALSE,
strip=FALSE,
#remove border around chart and axis lines
par.settings = list(axis.line = list(col = NA),
strip.border = list(col = NA),
strip.background = list(col = NA)),
lattice.options = list(axis.padding = list(numeric = 0)),
border = NA,
scales = list(tck = c(1,0), #remove tick lines on top
x = list(col.line="black"),
y = list(col.line="black", rot=0)), #make ticks black and not rotated
xlab = NULL
)