Monday, April 30, 2012

French Global Factors

I have said it already in multiple posts, but Kenneth French’s data library is one of the most generous and powerful contributions to the financial community.  To build on Systematic Investor’s series on factors, I thought I should run some basic analysis on the Global Factors maintained by Kenneth French.  I cannot imagine how long this would take without the data library and the incredible set of R packages available.

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

R code from GIST:

#get very helpful Ken French data
#for this project we will look at Global Factors
#http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/Global_Factors.zip
require(PerformanceAnalytics)
require(quantmod)
require(RColorBrewer)
#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/Global_Factors.zip"
#this will be the temp file set up for the zip file
my.tempfile<-paste(tempdir(),"\\frenchfactors.zip",sep="")
#my.usefile is the name of the txt file with the data
my.usefile<-paste(tempdir(),"\\Global_Factors.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_factor <- read.table(file=my.usefile,
header = TRUE, sep = "",
as.is = TRUE,
skip = 6, nrows=261)
#get dates ready for xts index
datestoformat <- rownames(french_factor)
datestoformat <- paste(substr(datestoformat,1,4),
substr(datestoformat,5,7),"01",sep="-")
#get xts for analysis
french_factor_xts <- as.xts(french_factor[,1:NCOL(french_factor)],
order.by=as.Date(datestoformat))
french_factor_xts <- french_factor_xts/100
#replace -0.9999 which means data does not exist
#know there is a better method to index, but can't find my easy approach
french_factor_xts[which(french_factor_xts < -0.99,arr.ind=TRUE)[,1],
unique(which(french_factor_xts < -0.99,arr.ind=TRUE)[,2])] <- 0
colnames(french_factor_xts) <- c("market","size","value","momentum","tbill")
chart.Correlation(french_factor_xts, main="Correlation of Ken French Global Market Factors")
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,line=4,cex=0.5, col="blue",adj=0)
chart.RollingCorrelation(french_factor_xts[,2:5],french_factor_xts[,1],
legend.loc="topleft",width=36,lwd=3,
main="Global Factor Rolling Correlation (3 Years)",
colorset=c("darkseagreen4","slateblue3","deepskyblue3","tan4"))
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,line=2.25,cex=0.5, col="blue",adj=0)
require(FactorAnalytics)
chart.RollingStyle(french_factor_xts[,1],french_factor_xts[,2:NCOL(french_factor_xts)],
width=12,
colorset=c("darkseagreen3","slateblue2","deepskyblue2","tan1"),
main="Global Market Rolling 1y French Factor Weights")
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,line=1,cex=0.5, col="blue",adj=0)
chart.Boxplot(french_factor_xts,main="Global Factors Return Distribution",
sort.by="",mean.symbol=19,symbol.color=c("gray60","darkseagreen4","slateblue3","deepskyblue3","tan3"),
colorset=c("gray60","darkseagreen4","slateblue3","deepskyblue3","tan3"))
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,cex=0.5, col="blue")
require(fPortfolio)
mycolors = c("gray60","darkseagreen4","slateblue3","deepskyblue3","tan4")
frontier <- portfolioFrontier(as.timeSeries(french_factor_xts))
pointsFrontier = frontierPoints(frontier, frontier = "both", auto=TRUE)
targetRisk = getTargetRisk(frontier@portfolio)[,1]
targetReturn = getTargetReturn(frontier@portfolio)[,1]
ans = cbind(Risk = targetRisk, Return = targetReturn)
colnames(ans) = c("targetRisk", "targetReturn")
rownames(ans) = as.character(1:NROW(ans))
#plot frontier points
#this method gives us much better control than frontierPlot
plot(ans,xlim=c(min(ans[,1]),max(ans[,1])+.01),ylim=c(0,0.0075),type="l",lwd=2, xlab=NA,ylab=NA)
minvariancePoints(frontier,pch=19,col="red")
tangencyPoints(frontier,pch=19,col="blue")
equalWeightsPoints(frontier,pch=15,col="grey")
singleAssetPoints(frontier,pch=19,cex=1.5,col=mycolors)
twoAssetsLines(frontier,lty=3,col="grey")
#label assets
stats <- getStatistics(frontier)
text(y=stats$mean,x=sqrt(diag(stats$Cov)),labels=names(stats$mean),pos=4,col=mycolors,cex=0.7)
#title(main="Efficient Frontier Small and Mid Since 1984")
#set up function from equalWeightsPoints to also label the point
equalLabel <- function (object, return = c("mean", "mu"), risk = c("Cov", "Sigma",
"CVaR", "VaR"), auto = TRUE, ...)
{
return = match.arg(return)
risk = match.arg(risk)
data = getSeries(object)
spec = getSpec(object)
constraints = getConstraints(object)
numberOfAssets = getNAssets(object)
setWeights(spec) = rep(1/numberOfAssets, times = numberOfAssets)
ewPortfolio = feasiblePortfolio(data, spec, constraints)
assets = frontierPoints(ewPortfolio, return = return, risk = risk,
auto = auto)
text(assets, labels = "Equal-Weight", pos=4,...)
invisible(assets)
}
equalLabel(frontier,cex=0.7,col="grey")
#title(main="Efficient Frontier 2000-October 2011",xlab="Risk(cov)",ylab="Monthly Return")
title(main="Efficient Frontier of French Global Factors",xlab="Risk(cov)",ylab="Monthly Return")
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,cex=0.5, col="blue",adj=0)
weightsPlot(frontier,col=mycolors)
mtext("Source: Kenneth French Data Library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/Data_Library",
side=1,cex=0.5, col="blue",adj=0)

Friday, April 27, 2012

Real Time Structural Break

Yesterday as I played with bfast I kept thinking “Yes, but this is all in hindsight.  How can I potentially use this in a system?”  Fortunately, one of the fine authors very generously commented on my post Structural Breaks (Bull or Bear?):

Jan Verbesselt Apr 27, 2012 02:01 AM

Nice application! you can also detect seasonal breaks. also check some new near real-time disturbance detection functionality using bfastmonitor() http://bfast.r-forge.r-project.org/Verbesselt+Zeileis+Herold-2012.pdf
cheers, Jan”

And away I went on an unexpected but very pleasant journey into bfastmonitor.  Please see the following paper for all the details.

Jan Verbesselt, Achim Zeileis, Martin Herold (2011). Near Real-Time Disturbance
  Detection in Terrestrial Ecosystems Using Satellite Image Time Series: Drought
  Detection in Somalia. Working Paper 2011-18. Working Papers in Economics and
  Statistics, Research Platform Empirical and Experimental Economics, Universitaet
  Innsbruck. URL http://EconPapers.RePEc.org/RePEc:inn:wpaper:2011-18

Doing a walk-forward test seemed like the best method of testing and illustration, so I chose the excruciating and incredibly volatile period from late 2008 to early 2009 as our example.  Amazingly, it picked with no optimizing or manual intervention March 2009 as the breakpoint. Of course, we would not know this until the end of March, but picking real-time with only a month lag is unbelievable to me.  Please try it out, and let me know your results.  Of course, I already have the 30 year bond bull in mind as a next trial.

Thanks to Yihui Xie who resurfaced again (see posts on knitr) with his animation package, which I used to create a good old-fashioned animated GIF.  I wish I had time to play more with the prettier and more robust options offered by the package.

animation

R code from GIST:

#analyze breakpoints in realtime with the R package bfast
#please read the paper
#Jan Verbesselt, Achim Zeileis, Martin Herold (2011). Near Real-Time Disturbance
#Detection in Terrestrial Ecosystems Using Satellite Image Time Series: Drought
#Detection in Somalia. Working Paper 2011-18. Working Papers in Economics and
#Statistics, Research Platform Empirical and Experimental Economics, Universitaet
#Innsbruck. URL http://EconPapers.RePEc.org/RePEc:inn:wpaper:2011-18
#http://scholar.google.com/scholar?cluster=9016488513865299942&hl=en&as_sdt=0,1
#install r-forge development version of bfast
#install.packages("bfast", repos="http://R-Forge.R-project.org")
require(bfast)
require(quantmod)
require(animation) #Yihui Xie surfaces again
getSymbols("^GSPC",from="1950-01-01")
#convert to log price
#for the sake of this example, let's assume we are in January 2009
#and we will see how this progresses with a for loop to go through
#the painful months of late 2008 and early 2009
#to see when our real time monitoring detects a break
#get a vector of the dates
evaldates <- c("2008-09","2008-10","2008-11","2008-12",
"2009-01","2009-02","2009-03","2009-04",
"2009-05","2009-06")
saveGIF(
for(i in 1:length(evaldates)) {
#notice this removes the foresight bias by only going to the current month
GSPC.monthly <- log(to.monthly(GSPC)[paste("1990::",evaldates[i],sep=""),4])
#need ts representation so do some brute force conversion
GSPC.ts <- ts(as.vector(GSPC.monthly),start=as.numeric(c(format(index(GSPC.monthly)[1],"%Y"),format(index(GSPC.monthly)[1],"%m"))),frequency=12)
#since we know the peak was October of 2010 let's use that
#as our start point for real time monitoring
mon5y <- bfastmonitor(GSPC.ts,
start=c(2007,10))
plot(mon5y)
mtext(evaldates[i],col="green",cex=1.5)
}
)
#really does not work, but always nice to have more examples
#let's get the minimum and maximum for our first experiment
GSPC.max <- GSPC.monthly[which(GSPC.monthly==max(last(GSPC.monthly,"10 years")))]
GSPC.min <- GSPC.monthly[which(GSPC.monthly==min(last(GSPC.monthly,"10 years")))]
#for evaluation let's pick whatever point is farthest from most recent close
GSPC.eval <- GSPC.monthly[GSPC.monthly == ifelse(as.numeric(last(GSPC.monthly))-as.numeric(GSPC.min) <
abs(as.numeric(last(GSPC.monthly)-as.numeric(GSPC.max))),GSPC.max,GSPC.min)]
#start monitoring from the max or min farthest away from current
mon5y <- bfastmonitor(GSPC.ts,
start=as.numeric(c(format(index(GSPC.eval),"%Y"),
format(index(GSPC.eval),"%m"))))

Thursday, April 26, 2012

Structural Breaks (Bull or Bear?)

When I spotted the bfast R package, I could not resist attempting to apply it to identify bull and bear markets.  For all the details that I do not understand, please see the references:

Jan Verbesselt, Rob Hyndman, Glenn Newnham, Darius Culvenor (2010). Detecting Trend and
  Seasonal Changes in Satellite Image Time Series. Remote Sensing of Environment, 114(1),
  106-115. doi:10.1016/j.rse.2009.08.014

Jan Verbesselt, Rob Hyndman, Achim Zeileis, Darius Culvenor (2010). Phenological Change
  Detection while Accounting for Abrupt and Gradual Trends in Satellite Image Time
  Series. Remote Sensing of Environment, 114(12), 2970 - 2980.
  doi:10.1016/j.rse.2010.08.003

I believe the result on the S&P 500 (even with a high h and only one iteration) is a fairly close marker for bull and bear markets, and I thoroughly enjoyed applying forestry techniques to financial time series.

From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio

R code from GIST:

#analyze breakpoints with the R package bfast
#please read the paper
#Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010)
#Detecting Trend and Seasonal Changes in Satellite Image Time Series.
#Remote Sensing of Environment, 114(1), 106–115.
#http://dx.doi.org/10.1016/j.rse.2009.08.014
require(bfast)
require(quantmod)
getSymbols("^GSPC",from="1950-01-01")
#convert to log price
GSPC.monthly <- log(to.monthly(GSPC)[,4])
#get monthly returns for the close price
#not necessary, leave in price form
#GSPC.return <- monthlyReturn(GSPC[,4])
#need ts representation so do some brute force conversion
GSPC.ts <- ts(as.vector(GSPC.monthly["1951-01::"]),start=c(1951,1),frequency=12)
#look at the stl Seasonal-Trend decomposition procedure already in R
GSPC.stl <- stl(GSPC.ts,s.window="periodic")
plot(GSPC.stl,main="STL Decomposition of S&P 500")
#get the results from bfast
#adjusting h lower will result in more breakpoints
GSPC.bfast <- bfast(GSPC.ts,h=0.2,max.iter=1,season="none")
plot(GSPC.bfast,type="components",ylim=c(3,max(GSPC.monthly)+1),main="S&P 500 with bfast Breakpoints and Components")
plot(GSPC.bfast,type="trend",ylim=c(3,max(GSPC.monthly)+1),main="S&P 500 with bfast Trend Breakpoints")
#see everything with type="all" but in bfast calculation set seasonal to "none"
#play away with this
#plot(GSPC.bfast,type="all")
#do some additional plotting
#[[1]] is used since for speed I only did one iteration
#could plot each iteration if I did multiple
plot(GSPC.bfast$Yt/GSPC.bfast$output[[1]]$Tt-1,
main="bfast Remainder as % of S&P 500 Price",
xlab=NA, ylab="remainder (% of price)",bty="l")
#add vertical line for the breakpoints
abline(v=breakdates(GSPC.bfast$output[[1]]$bp.Vt),col="gray70")
#add horizontal line at 0
abline(h=0,col="black",lty=2)
text(x=breakdates(GSPC.bfast$output[[1]]$bp.Vt),y=par("usr")[3]+.01,
labels=breakdates(GSPC.bfast$output[[1]]$bp.Vt,format.times=TRUE),
srt=90,pos=4,cex=0.75)

Monday, April 23, 2012

Drawdown Look at Frontier of Assets and Systems

In Efficient Frontier of Funds and Allocation Systems, I had hoped to start exploring how a frontier can potentially be created with only one asset, or how an even more efficient frontier could be created with assets and also systems on those assets.  I am obsessed with drawdown, so of course I need to extend that post with a look at drawdown to satisfy this obsession. If you have not read the original post, please read it before this, since I will only show the additional graph created.

From TimelyPortfolio

R code from GIST:

require(quantmod)
require(PerformanceAnalytics)
getSymbols("VFINX",from="1990-01-01",adjust=TRUE)
getSymbols("VBMFX",from="1990-01-01",adjust=TRUE)
perf <- na.omit(merge(monthlyReturn(VBMFX[,4]),monthlyReturn(VFINX[,4])))
colnames(perf) <- c("VBMFX","VFINX")
#get 8 month RSI; randomly picked 8; no optimization
rsi<- lag(merge(RSI(perf[,1],n=8),RSI(perf[,2],n=8)),k=1)
#allocate between vbmfx and vfinx based on highest RSI
rsi.perf <- ifelse(rsi[,1]>rsi[,2],perf[,1],perf[,2])
rsi.each <- as.xts(as.matrix(rsi>50) * as.matrix(perf),
order.by=index(perf))
#get cumulative returns for moving average
cumul <- as.xts(apply(perf+1,MARGIN=2,cumprod),order.by=index(perf))
#do 10 month Mebane Faber style system
ma <- lag(merge(runMean(cumul[,1],n=10),runMean(cumul[,2],n=10)),k=1)
#apply 50% allocation to each fund if they are > 10 month moving average
ma.perf <- as.xts(apply(as.matrix(cumul>ma) * as.matrix(perf)/2,
MARGIN=1,sum),
order.by=index(perf))
ma.each <- as.xts(as.matrix(cumul>ma) * as.matrix(perf),
order.by=index(perf))
#add omega as another allocation method
omega <- lag(merge(apply.rolling(perf[,1],width=6,by=1,FUN=Omega),
apply.rolling(perf[,2],width=6,by=1,FUN=Omega)),
k=1)
#if omega >= 1 then apply 50% allocation
omega.perf <- as.xts(apply(as.matrix(omega>=1) * as.matrix(perf)/2,
MARGIN=1,sum),
order.by=index(perf))
omega.each <- as.xts(as.matrix(omega>=1) * as.matrix(perf),
order.by=index(perf))
perf.all <- merge(perf,rsi.perf,rsi.each,ma.perf,ma.each,omega.perf,omega.each)
perf.all[is.na(perf.all)]<-0
colnames(perf.all) <- c(colnames(perf),paste(c(rep("rsi",3),rep("ma",3),rep("omega",3)),
c("",".VBMFX",".VFINX"),sep=""))
#now let's add two very basic systems
#and explore on Systematic Investor's efficient frontier
########################################################
#continue to highlight the very fine work of
#http://systematicinvestor.wordpress.com/
#adapted some of his code to provide
#a not-so-novel additional example for
#those that might be interested
#######################################################
# Load Systematic Investor Toolbox (SIT)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
#--------------------------------------------------------------------------
# Create Efficient Frontier
#--------------------------------------------------------------------------
ia = list()
#amend to use the funds and basic systems
ia$symbols = colnames(perf.all)
ia$n = len(ia$symbols)
#use PerformanceAnalytics tables to get return (geometric) and risk
#for the entire period
ia$expected.return = as.matrix(t(table.Stats(perf.all)[7,]))
ia$risk = as.matrix(t(table.Stats(perf.all)[14,]))
ia$correlation = cor(perf.all)
ia$cov = cov(perf.all)
n = ia$n
# 0 <= x.i <= 1
constraints = new.constraints(n, lb = 0, ub = 1)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
# create efficient frontier
ef.risk = portopt(ia, constraints, 50)
#I am getting an error here
#plot.ef(ia, ef.risk), transition.map=TRUE)
#know what is happening but not motivated to fix
#"Error in x$weight : $ operator is invalid for atomic vectors"
#will do manually plot
colors <- c("purple","indianred3","steelblue2","steelblue3","steelblue4",
"darkolivegreen2","darkolivegreen3","darkolivegreen4",
"chocolate2","chocolate3","chocolate4")
plot(ef.risk$return~ef.risk$risk,col="grey60",lwd=3,type="l",
xlim=c(min(ia$risk),max(ia$risk)+.01),
ylim=c(min(ia$expected.return),max(ia$expected.return)))
points(x=as.numeric(ia$risk),y=as.numeric(ia$expected.return),pch=19,
col=colors,cex=1.5)
text(x=as.numeric(ia$risk),y=as.numeric(ia$expected.return),
labels=ia$symbols,pos=4,col=colors)
title(main="Efficient Frontier of VBMFX and VFINX and Systematic Allocation",
adj=0,outer=TRUE,line=-1)
plot.transition.map(ef.risk,col=colors)
chart.CumReturns(perf.all,colorset=colors,
main="Growth of VBMFX and VFINX and Systematic Allocations",
legend.loc="topleft")
#I am sure there is a better way to do this
#this function will calculate the maxDrawdown for each of the point allocations
#from the systematic investor (SIT) frontier calculated above
#if there is not a better way then you're welcome
frontier.drawdown <- function(weight,perfxts) {
point.drawdown <- matrix(nrow=NROW(weight))
for (i in 1:NROW(weight)) {
weight.matrix <- matrix(rep(weight[i,],NROW(perfxts)),byrow=TRUE,nrow=NROW(perfxts),ncol=NCOL(weight))
point.drawdown[i] <- maxDrawdown(as.matrix(apply((weight.matrix * perfxts),MARGIN=1,sum)))
}
return(point.drawdown)
}
drawdown <- frontier.drawdown(weight=ef.risk$weight,perfxts=perf.all)
#set outer margin on top to get the title in a better spot
par(oma=c(0,1,2,0))
#do frontier style plot with drawdown as the risk measure
plot(-t(maxDrawdown(perf.all)),ia$expected.return,
col=colors,xlim=c(-0.6,0),pch=19,cex=1.25,
xlab="drawdown",ylab="monthly return",bty="u")
#add a right side axis since 0 is the origin and drawdown is negative
axis(side=4,lwd=0,lwd.ticks=1)
#add labels for each of the points
text(x=-t(maxDrawdown(perf.all)),y=ia$expected.return,labels=ia$symbols,pos=2,col=colors,cex=0.75)
#add the frontier line generated from normal optimization of mean and return
points(x=-drawdown,y=ef.risk$return,type="l",lwd=2,col="grey50")
title(main="Efficient Frontier with Drawdown as Risk Measure",adj=0,outer=TRUE)

Wednesday, April 18, 2012

Efficient Frontier of Funds and Allocation Systems

I did a very basic experiment in Efficient Frontier of Buy-Hold and Tactical System where I determined the efficient frontier of the S&P 500 with itself transformed by a Mebane Faber 10-month moving average tactical allocation.

The result was interesting, but I did not pursue further.  Now with some inspiration and tools by Systematic Investor, I thought I would extend the post. This time around we will use both the Vanguard U.S. Total Bond Market (VBMFX) and Vanguard U.S. S&P 500 (VFINX) combined with both portfolios determined by tactical methods (moving average, RSI, and omega) and those funds transformed individually by the same tactical methods.  I will as always warn you that this is not advice and large losses are almost guaranteed.  Also, I would like to note that I have checked the 10-month moving average every way possible (even manually in Excel), and it has been incredible on the VFINX since 1990.  Prior to 1990, results were good but nowhere near as amazing.  If I messed up, please let me know.

From TimelyPortfolio
From TimelyPortfolio
From TimelyPortfolio

R code from GIST:

require(quantmod)
require(PerformanceAnalytics)
getSymbols("VFINX",from="1990-01-01",adjust=TRUE)
getSymbols("VBMFX",from="1990-01-01",adjust=TRUE)
perf <- na.omit(merge(monthlyReturn(VBMFX[,4]),monthlyReturn(VFINX[,4])))
colnames(perf) <- c("VBMFX","VFINX")
#get 8 month RSI; randomly picked 8; no optimization
rsi<- lag(merge(RSI(perf[,1],n=8),RSI(perf[,2],n=8)),k=1)
#allocate between vbmfx and vfinx based on highest RSI
rsi.perf <- ifelse(rsi[,1]>rsi[,2],perf[,1],perf[,2])
rsi.each <- as.xts(as.matrix(rsi>50) * as.matrix(perf),
order.by=index(perf))
#get cumulative returns for moving average
cumul <- as.xts(apply(perf+1,MARGIN=2,cumprod),order.by=index(perf))
#do 10 month Mebane Faber style system
ma <- lag(merge(runMean(cumul[,1],n=10),runMean(cumul[,2],n=10)),k=1)
#apply 50% allocation to each fund if they are > 10 month moving average
ma.perf <- as.xts(apply(as.matrix(cumul>ma) * as.matrix(perf)/2,
MARGIN=1,sum),
order.by=index(perf))
ma.each <- as.xts(as.matrix(cumul>ma) * as.matrix(perf),
order.by=index(perf))
#add omega as another allocation method
omega <- lag(merge(apply.rolling(perf[,1],width=6,by=1,FUN=Omega),
apply.rolling(perf[,2],width=6,by=1,FUN=Omega)),
k=1)
#if omega >= 1 then apply 50% allocation
omega.perf <- as.xts(apply(as.matrix(omega>=1) * as.matrix(perf)/2,
MARGIN=1,sum),
order.by=index(perf))
omega.each <- as.xts(as.matrix(omega>=1) * as.matrix(perf),
order.by=index(perf))
perf.all <- merge(perf,rsi.perf,rsi.each,ma.perf,ma.each,omega.perf,omega.each)
perf.all[is.na(perf.all)]<-0
colnames(perf.all) <- c(colnames(perf),paste(c(rep("rsi",3),rep("ma",3),rep("omega",3)),
c("",".VBMFX",".VFINX"),sep=""))
#now let's add two very basic systems
#and explore on Systematic Investor's efficient frontier
########################################################
#continue to highlight the very fine work of
#http://systematicinvestor.wordpress.com/
#adapted some of his code to provide
#a not-so-novel additional example for
#those that might be interested
#######################################################
# Load Systematic Investor Toolbox (SIT)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
#--------------------------------------------------------------------------
# Create Efficient Frontier
#--------------------------------------------------------------------------
ia = list()
#amend to use the funds and basic systems
ia$symbols = colnames(perf.all)
ia$n = len(ia$symbols)
#use PerformanceAnalytics tables to get return (geometric) and risk
#for the entire period
ia$expected.return = as.matrix(t(table.Stats(perf.all)[7,]))
ia$risk = as.matrix(t(table.Stats(perf.all)[14,]))
ia$correlation = cor(perf.all)
ia$cov = cov(perf.all)
n = ia$n
# 0 <= x.i <= 1
constraints = new.constraints(n, lb = 0, ub = 1)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
# create efficient frontier
ef.risk = portopt(ia, constraints, 50)
#I am getting an error here
#plot.ef(ia, ef.risk), transition.map=TRUE)
#know what is happening but not motivated to fix
#"Error in x$weight : $ operator is invalid for atomic vectors"
#will do manually plot
colors <- c("purple","indianred3","steelblue2","steelblue3","steelblue4",
"darkolivegreen2","darkolivegreen3","darkolivegreen4",
"chocolate2","chocolate3","chocolate4")
plot(ef.risk$return~ef.risk$risk,col="grey60",lwd=3,type="l",
xlim=c(min(ia$risk),max(ia$risk)+.01),
ylim=c(min(ia$expected.return),max(ia$expected.return)))
points(x=as.numeric(ia$risk),y=as.numeric(ia$expected.return),pch=19,
col=colors,cex=1.5)
text(x=as.numeric(ia$risk),y=as.numeric(ia$expected.return),
labels=ia$symbols,pos=4,col=colors)
title(main="Efficient Frontier of VBMFX and VFINX and Systematic Allocation",
adj=0,outer=TRUE,line=-1)
plot.transition.map(ef.risk,col=colors)
chart.CumReturns(perf.all,colorset=colors,
main="Growth of VBMFX and VFINX and Systematic Allocations",
legend.loc="topleft")

knitr Performance Report-Attempt 2

Over the years I have changed my learning process from reading thoroughly first before proceeding to reading minimally and then applying immediately.  I very quickly see the gaps in my knowledge.  This method is far more painful but seems to quicken progress up the learning curve.  You will definitely see the process in its entirety with this series on latex, knitr, and performance reporting.  For those who are still hesitating, jump right in with me and let me know your result.

As I experimented with Attempt 2, I could see from others sample .Rnw files that LyX can serve as the IDE for Sweave/knitr documents, and I thought this would be a more pleasant route.  However it was not, and I eventually reverted back to old-school manual coding.  This post from R, Ruby, and Finance talks about TeXnic Center as an IDE.  Maybe I will try it for attempt 3. I am still a long way from a useable result (any bets on how many attempts it takes?), but the report is definitely an improvement

In Attempt 2, I used the echo=FALSE option to not output the R code on the final pdf report.  However, the code is still all there for public view within the .Rnw file.

I just remembered that ttrTests offers a Latex output option. Just what I need--another distraction. Look for that in future posts.

R code from GIST:
require(knitr)
knit2pdf("knitr performance 2.rnw")
%% LyX 2.0.2 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,nohyper,noae]{tufte-handout}
\usepackage{helvet}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{babel}
\usepackage[unicode=true,pdfusetitle,
bookmarks=true,bookmarksnumbered=true,bookmarksopen=true,bookmarksopenlevel=1,
breaklinks=true,pdfborder={0 0 0},backref=false,colorlinks=false]
{hyperref}
\usepackage{breakurl}
\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
\title{knitr Performance Summary Attempt 2}
\author{Timely Portfolio}
\providecommand{\LyX}{\texorpdfstring%
{L\kern-.1667em\lower.25em\hbox{Y}\kern-.125emX\@}
{LyX}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
<<echo=F>>=
if(exists(".orig.enc")) options(encoding = .orig.enc)
@
\makeatother
\begin{document}
\maketitle
\begin{abstract}
This time we will try to incorporate some \LyX{} brilliance to help
us format the end result. Unfortunately, this process (most likely
due to my ignorance) was not nearly as seamless as the first experiments
with knitr. I eventually went back to old-school manual coding in RStudio.
\end{abstract}
\section{Performance Overview}
We all know that performance reports generally blend both text tables
and graphical charts to communicate. In this example, we will use
two of the most popular Vanguard funds (vbmfx and vfinx) as the subjects
for our performance evaluation.
<<setup,eval=TRUE,echo=FALSE, warning=FALSE,message=FALSE, results='hide'>>=
require(quantmod)
require(PerformanceAnalytics)
getSymbols("VFINX",from="1990-01-01",adjust=TRUE)
getSymbols("VBMFX",from="1990-01-01",adjust=TRUE)
perf <- na.omit(merge(monthlyReturn(VBMFX[,4]),monthlyReturn(VFINX[,4])))
colnames(perf) <- c("VBMFX","VFINX")
@
<<lattice-chart,echo=FALSE,eval=FALSE,tidy=TRUE,warning=FALSE,message=FALSE>>=
require(lattice)
require(latticeExtra)
require(reshape2)
#note: table.CalendarReturns only handles monthly data
perf.annual<-as.data.frame(table.CalendarReturns(perf)[,13:14])
perf.annual<-cbind(rownames(perf.annual),perf.annual)
perf.annual.melt <- melt(perf.annual,id.vars=1)
colnames(perf.annual.melt)<-c("Year","Fund","Return")
p1 <- dotplot(Year~Return,group=Fund,data=perf.annual.melt,
pch=19,
lattice.opts=theEconomist.opts(),
par.settings = theEconomist.theme(box = "transparent"),
main="Annual Returns of VFINX and VBMFX",
auto.key=list(space="right"),
xlim=c(min(perf.annual.melt[,3]),max(perf.annual.melt[,3])))
p2 <- densityplot(~Return, group=Fund, data=perf.annual.melt,
lattice.opts=theEconomist.opts(),
par.settings = theEconomist.theme(box = "transparent"))
print(p1,position=c(0,0,0.6,1),more=TRUE)
print(p2+p1,position=c(0.6,0,1,1))
@
<<ref.label='lattice-chart',echo=FALSE,fig.height=4.75,fig.width=7,warning=FALSE,message=FALSE>>=
@
<<xtable-table,echo=FALSE,eval=TRUE,results='tex',tidy=TRUE,warning=FALSE,message=FALSE>>=
require(xtable)
print(xtable(t(last(table.CalendarReturns(perf)[,13:14],13))), floating=FALSE)
@
\newpage
\section{Risk and Return}
Although the summary and distribution of annual returns is a good first step, any real due diligence will require much more than just return. Let's do a very basic plot of risk and return. Of course, there are much more sophisticated methods, which we will explore in future versions.
\begin{figure}
<<chart-stats,echo=FALSE,eval=TRUE,tidy=TRUE,warning=FALSE,message=FALSE,fig.height=4.75,fig.width=7>>=
perf.stats <- table.Stats(perf)
#eliminate observations,na,skewness,and kurtosis
perf.stats <- perf.stats[3:(NROW(perf.stats)-2),]
perf.stats.melt <- melt(as.data.frame(cbind(rownames(perf.stats),perf.stats),stringsAsFactors=FALSE),id.vars=1)
colnames(perf.stats.melt)<-c("Statistic","Fund","Value")
barchart(Statistic~Value,group=Fund,data=perf.stats.melt,
origin=0,
lattice.opts=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
main="Risk and Return Statistics")
@
\end{figure}
\section{Diversification}
In today's markets with almost universally high positive correlations, diversification is much more difficult. However, most performance reports should include some analysis of both correlation and diversification.
\begin{figure}
<<chart-correlation,echo=FALSE,eval=TRUE,tidy=TRUE,warning=FALSE,message=FALSE,fig.height=4,fig.width=6>>=
chart.Correlation(perf, main="Correlation Analysis")
@
<<chart-rollingcorrelation,echo=FALSE,eval=TRUE,tidy=TRUE,warning=FALSE,message=FALSE,fig.height=4,fig.width=6>>=
chart.RollingCorrelation(perf[,1],perf[,2],main="VBMFX and VFINX Rolling 12 Month Correlation",xlab="Correlation")
@
\end{figure}
\begin{figure}
<<chart-diversification,echo=FALSE,eval=TRUE,tidy=TRUE,warning=FALSE,message=FALSE,fig.height=5,fig.width=6,fig.keep="last">>=
require(fPortfolio)
frontier <- portfolioFrontier(as.timeSeries(perf))
frontierPlot(frontier,pch=19,title=FALSE)
singleAssetPoints(frontier,col=c("steelblue2","steelblue3"),pch=19,cex=2)
title(main="Efficient Frontier with VBMFX and VFINX",adj=0)
@
\end{figure}
\end{document}

Friday, April 13, 2012

knitr Performance Report-Attempt 1

I get very excited about new R packages, but rarely is my excitement so fulfilled as with knitr.  Even with no skill, I have already been able to adapt the example Yihui Xie provides in his knitr Graphics Manual into a crude first version of a performance report that I could actually show clients and prospects.  Although this is far from production quality, two days of experimentation already has me to a level that assures me of the wonderful potential of combining the existing amazing R packages with knitr.  Before knitr, I do not believe I could have accomplished even this rough first draft.

If you have not played with MiKTeX before, you will need to use the Package Manager to install the xcolor and tufte-handout templates for this example to work properly.  MiKTeX is installed automatically with LyX as discussed in yesterday’s post Latex Allergy Cured by knitr.  If xcolor causes problems, then just change your repository to another repository.

imageimage

knitr weaved with a little PerformanceAnalytics, lattice, and latticeExtra provides this first draft.  As always, thanks so much to the brilliant contributors that have made this possible.

R code in GIST:

require(knitr)
knit2pdf("performance summary.rnw")
%% LyX 2.0.3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[nohyper,justified]{tufte-handout}
\usepackage[T1]{fontenc}
\usepackage{url}
\usepackage[unicode=true,pdfusetitle, bookmarks=true,bookmarksnumbered=true,bookmarksopen=true,bookmarksopenlevel=2, breaklinks=true,pdfborder={0 0 1},backref=false,colorlinks=false] {hyperref}
\hypersetup{pdfstartview=FitH}
\usepackage{breakurl}
\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
\title{Performance Reporting with knitr}
\author{Timely Portfolio}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\renewcommand{\textfraction}{0.05}
\renewcommand{\topfraction}{0.8}
\renewcommand{\bottomfraction}{0.8}
\renewcommand{\floatpagefraction}{0.75}
\makeatother
\begin{document}
% \SweaveOpts{fig.path='figure/graphics-', cache.path='cache/graphics-', fig.align='center', dev='pdf', fig.width=5, fig.height=5, fig.show='hold', cache=TRUE, par=TRUE}
<<setup, include=FALSE, cache=FALSE>>=
options(replace.assign=TRUE,width=60)
@
\maketitle
\begin{abstract}
This sample performance report will begin to highlight the ability of the \textbf{knitr} package to generate marketing materials and client-facing performance reports with a little help from the \textbf{PerformanceAnalytics} package.
\end{abstract}
Thanks again to Yihui Xie for not only his amazing \textbf{knitr} work but also his numerous examples. His {\url{http://yihui.name/knitr/demo/graphics/} (\textbf{knitr Graphics Manual}) will provide the initial template for this report. As I learn, hopefully I will not have to mimic his example so closely.
\section{Performance Summary}
For this first example we will use the prebuilt \textbf{charts.PerformanceSummary} function to visualize cumulative growth and drawdown of the EDHEC style indexes provided by data(edhec). Although \textbf{charts.PerformanceSummary} was primarily intended as an example or template, I hear that it has appeared unadulterated in live performance reports and marketing.
<<loadlibrary, eval=TRUE, results='hide', message=FALSE>>=
require(PerformanceAnalytics)
data(edhec)
@
\begin{figure}
<<perf, fig.width=8, fig.height=8, out.width='.9\\linewidth',dev='pdf',echo=TRUE,warning=FALSE>>=
charts.PerformanceSummary(edhec,main="Performance of EDHEC Style Indexes")
@
\caption{\textbf{charts.PerformanceSummary} provides a nice chart of my favorite measures: compounded return and drawdown.\label{fig:perf}}
\end{figure}
\newpage
\section{Improvement??}
With a little help from \textbf{lattice} and \textbf{latticeExtra}, maybe we can get something that might fit my style a little better.
<<perfbetter, eval=FALSE, warning=FALSE>>=
require(lattice)
require(latticeExtra)
require(reshape2)
#get cumulative growth of $1
edhec.cumul <- apply(edhec+1,MARGIN=2,cumprod)
#use melt so we can get in a format lattice or ggplot2 likes
edhec.cumul.melt <- melt(as.data.frame(cbind(index(edhec),edhec.cumul)),id.vars=1)
#name columns something more appropriate
colnames(edhec.cumul.melt) <- c("Date","Style","Growth")
#get dates in text form
edhec.cumul.melt[,1] <- as.Date(edhec.cumul.melt[,1])
colors <- c(brewer.pal(9,name="PuBuGn")[3:9],brewer.pal(9,"YlOrRd")[4:9])
#plot with lattice
xyplot(Growth~Date,groups=Style,data=edhec.cumul.melt,
type="l",lwd=2,col=colors,
par.settings=theEconomist.theme(box="transparent"),
axis = theEconomist.axis,
scales=list(x=list(alternating=1),
y=list(alternating=1)),
main="Cumulative Growth of EDHEC Style Indexes")+ layer(panel.text(x=as.Date(index(edhec)[NROW(edhec)]-1),y=round(edhec.cumul[NROW(edhec.cumul),],2),colnames(edhec),pos=0,cex=0.8,col=colors))
@
\begin{figure}
<<latticeperf,ref.label='perfbetter', fig.width=8, fig.height=8, out.width='.9\\linewidth',dev='pdf',echo=FALSE,warning=FALSE>>=
@
\caption{Still a mess but definitely closer to what I expect for production quality reporting. Keep following. I will get better.\label{fig:perf}}
\end{figure}
\end{document}

Wednesday, April 11, 2012

The Making of…R Fell on Alabama

Slightly off topic but if you are as easily distracted by amazing technology as I am then you might enjoy my short (all made on my Iphone 4s) “The Making of…R Fell on Alabama.”  Clearly I am not accustomed to doing this type of video.

Latex Allergy Cured by knitr

I have always known that at some point I would have to succumb to the power of Latex, but Latex has been uncharacteristically intimidating to me.  I finally found the remedy to my Latex allergy with the amazing and fantastic knitr package from Yihui Xie.  With very minimal effort, I ran my first experiment and now am extremely excited to incorporate it in production-quality performance reports (I plan to document the steps to get there in future posts).

For those starting from scratch on Windows, I think the easiest method to get up and running is to install LyX, which will also install MikTex.  If these are successfully installed, then you should be ready to experiment with knitr in R.

I will use knitr’s stich function, which is clearly not designed for the robust production use of knitr, but makes for a very easy first test.  stitch will open a very short script, apply a template, and generate a Sweave style .Rnw (can be changed).  knit2pdf converts the .Rnw file into a pdf, and with a couple lines of code you get a remarkable result.

 

 

R code from GIST (unbelievably only 7 lines):

require(knitr)
#use stitch to open and run the knitr performance.r script
# and then generate Rnw code
#knit2pdf converts the Rnw to a pdf
knit2pdf(stitch("knitr performance.r"))
view raw knitr stitch.r hosted with ❤ by GitHub
## title:First Attempt at a knitr Performance Report
## author:Timely Portfolio
require(PerformanceAnalytics)
require(xtable)
data(edhec)
table.CalendarReturns(edhec[,1:3])[13:15]
charts.PerformanceSummary(edhec[,1:3],
colorset=c("gray50","steelblue2","steelblue4"),
main="EDHEC Performance by Type")

Monday, April 9, 2012

Piggybacking and Hopefully Publicizing R Experts

I was inspired by the Revolution Analytics blog post http://blog.revolutionanalytics.com/2009/11/charting-time-series-as-calendar-heat-maps-in-r.html on the d3.js style calendar heat map that Paul Bleicher from Humedica developed in R.  In an effort to publicize such fine work, I wanted to put a slightly different spin on it to visualize a system’s time in the market.  The system is ridiculously basic (not investment advice and will lose money), but the visualization is very helpful.

From TimelyPortfolio

To continue with the theme, I would like to continue to highlight some fine R work from http://systematicinvestor.wordpress.com.  Although his work is far more sophisticated than this, I thought I would use his plota function to plot the German Dax (used in this example) with RSI below.

From TimelyPortfolio

Thanks so much to the amazing R coders who provided the code and functions for this post.

R code from GIST:

#quick example of d3.js style calendar heat map
#thanks so much Paul Bleicher, MD PhD who is Chief Medical Officer
# at Humedica, a next-generation clinical informatics company
# that provides novel business intelligence solutions
# to the health care and life science industries
#discussed in Revolution Analytics blog post
#http://blog.revolutionanalytics.com/2009/11/charting-time-series-as-calendar-heat-maps-in-r.html
con=url("http://blog.revolution-computing.com/downloads/calendarHeat.R")
source(con)
close(con)
require(quantmod)
#just as an extremely simple example
#enter when 20 day RSI > 50 and exit when 20 day RSI < 50
#weight will be 0 or 1
getSymbols("^GDAXI",from="1900-01-01",to=Sys.Date())
#get 20 day RSI and lag one day
signal <- lag(RSI(GDAXI[,4],n=20),k=1)
weights <- ifelse(signal>50,1,0)
calendarHeat(as.Date(index(weights["2009::",])),coredata(weights["2009::"]),varname="RSI System Weight",ncolors=2)
#######################################################
#continue to highlight the very fine work of
#http://systematicinvestor.wordpress.com/
#adapted some of his code to provide
#a not-so-novel additional example for
#those that might be interested
#######################################################
# Load Systematic Investor Toolbox (SIT)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
#set up two regions for graphs
#candlestick price data on top 2 of 3
#and rsi on the bottom 1 of 3
layout(c(1,1,2))
#plot candestick of the daily gdaxi data
plota(GDAXI, type = 'candle', plotX = F)
#add description of the data in top right and last value
plota.legend('German Dax', 'grey70', GDAXI)
#plot the rsi as line
plota(rsi<-RSI(GDAXI[,4],n=20), type = 'l')
#get some transparency for the highlighted regions
addalpha <- function(cols,alpha=180) {
rgbcomp <- col2rgb(cols)
rgbcomp[4] <- alpha
return(rgb(rgbcomp[1],rgbcomp[2],rgbcomp[3],rgbcomp[4],maxColorValue=255))
}
#since our entry is RSI>50 we will highlight the region > 50 in green
plota.y.highlight(col=addalpha("green",80),highlight=c(50,100))
#since our exit is RSI<50 we will highlight the region < 50 in red
plota.y.highlight(col=addalpha("red",80),highlight=c(0,50))
#plot a line our at special 50
abline(h = 50, col = 'gray20')
#give description of data in bottom (RSI) and color code to in (black) or out(red)
plota.legend('RSI(20)', text.col = ifelse(last(rsi)[]>50,'black','red'), fill = ifelse(last(rsi)[]>50,'black','red'), rsi)