The Financial Crisis on Tape Part II

FinCrisisPartII.R

rm(list=ls())
#install the superb quantmod library
#we will use it to download the data and compute returns
library(quantmod)
#install various other libraries for the data manipulation and graphing.
library(ggplot2)
library(scales)
library(RColorBrewer)
library(reshape2)
library(grid)
library(gridExtra)

# load historical prices from Yahoo Finance
# I use a set that I saw used by Systematic Investor
symbols = c('SPY','QQQ','EEM','IWM','EFA','TLT','IYR','GLD')
symbols.names = c('S&P 500,Nasdaq 100,Emerging Markets,Russell 2000,EAFE,20 Year
                  Treasury,U.S. Real Estate,Gold')

#Download the data from yahoo (default choice of getSymbols)
#It loads directly to the environment which depending on your R background may seem surprising
#One can use the 'get' function below to obtain the results.
getSymbols(symbols, src = 'yahoo', from = '2005-01-01')

#obtain the daily closing price for each and form into a data frame
#and compute the returns of the adjusted prices
hist.returns = 
  do.call(cbind,
          lapply(symbols, function(symbol){
            symbol.data = get(symbol) #get yahoo data from the environment
            symbol.data.adj = Ad(symbol.data) #extract the adjusted prices
            symbols.data.adjret = ROC(symbol.data.adj, n=1, type="continuous", na.pad=TRUE) #compute the simple returns
            symbols.data.adjret
          })
  )

#we will also need to get the plain adjusted prices which we will use later
hist.prices = 
  do.call(cbind,
          lapply(symbols, function(symbol){
            symbol.data = get(symbol) #get yahoo data from the environment
            symbol.data.adj = Ad(symbol.data) #extract the adjusted price
            symbol.data.adj
          })
  )

#a little prep for plotting with ggplot later
hist.prices = data.frame(this.date = index(hist.prices),hist.prices)
colnames(hist.prices) = gsub(".Adjusted","",colnames(hist.prices))
hist.prices.melt = melt(hist.prices,id.vars="this.date")

#define preferred order of the assets in the graphs.
pref.order = c("SPY","QQQ","EEM","IWM","EFA","IYR","GLD","TLT")
hist.prices.melt$variable  = factor(hist.prices.melt$variable,pref.order)

###########################################################
##############Compute Rolling Correlations#################
###########################################################
#there are many ways to apply a function to a rolling historical window
#I'm going to roll a very simple one of my own

#Function to compute the correlation matrix from an dataframe of returns:
DataFrameCorOutput<-function(hist.returns){
  require(reshape2)
  correls = melt(cor(as.matrix(na.omit(hist.returns)))) #r's built in correlation function returns the correlation matrix
  colnames(correls) = c("Var1","Var2","Correl") #label the melted correlation matrix

  #rename some of the series to make obtaining a pretty plot a little easier!
  correls$Var1 = factor(gsub(".Adjusted","",correls$Var1),rev(pref.order))
  correls$Var2 = factor(gsub(".Adjusted","",correls$Var2),rev(pref.order))
  
  list(correls = correls,min.date = min(index(hist.returns)),max.date = max(index(hist.returns))) #return all three objects
}

#choose the window length over which we will compute the correlation
window.length = 120 #~6M of business days

#calculate the dates on which we have sufficient data to compute the correlation (i.e. a window.length of history)
dates.to.do = tail(index(hist.returns),-(window.length - 1)) #these are the dates to right of the first window

#calculate the correlation for an N day window before each of these dates
rolling.correls = lapply(dates.to.do,FUN<-function(this.date){
                       return(DataFrameCorOutput(tail(hist.returns[index(hist.returns)<=this.date,],window.length)))
})
names(rolling.correls) = dates.to.do

#function to extract the correlation matrix from rolling.correls 
#and plot it using the ggplot package.
PlotCorrelsForThisDate<-function(this.date){
  #extract the relevant correlmatrix from the list of correl matrices
  correlmatrix = rolling.correls[[as.character(this.date)]][["correls"]]
  
  #I want to use a bespoke colour pallete see Part I of "The Financial Crisis on Tape" for an explanation.
  my.col.vec = c(
    rev(brewer.pal(8,"Blues")),
    colorRampPalette(
      c(
      colorRampPalette(c("white",brewer.pal(9,"Set3")[2],
                         brewer.pal(9,"Set1")[c(6,5,1)]))(10),
      rev(brewer.pal(9,"PRGn")[1])
      )
    )(9)
  )
  
  #plot a simple heat-map using ggplot
  correl.plot<-ggplot(correlmatrix, aes(x = Var1, y = Var2, fill = Correl)) + 
    geom_tile() + 
    scale_fill_gradientn(colours = my.col.vec
                         ,limits = c(-1,1)) +
    theme_bw()+
    theme(panel.grid.major.x = element_blank(),
          panel.grid.major.y = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x=element_text(angle=90))  
  
  return(correl.plot)
}

#The below is a funtion that given an index, j, produces the graphic for the date dates.to.do[j]
ProduceFilmStill<-function(j){
  #the main reason I use j rather than lapplying to a list of dates
  #is I wish to call the output Plot_j.png, so I need to record j.
  this.date = dates.to.do[j]
  
  #use the previous function to produce the correlation heatmap
  correlplot = PlotCorrelsForThisDate(this.date)

  #the strip.data highlights the period of time used to compute the correlation
  strip.data = data.frame(xmin=rolling.correls[[as.character(this.date)]][["min.date"]], xmax=this.date,ymin=-Inf,ymax=Inf)
  
  #Plot the timeseries plots:
  #the idea of these is that the dot shows the current point in time whilst the period used is shaded grey.
  #this gives context to the the correlation calculation
  timeseries.plot = ggplot(hist.prices.melt,aes(x=this.date,y=value,colour=variable))+
    geom_line()+
    geom_rect(data=strip.data,aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),fill="grey", alpha=0.2, inherit.aes = FALSE)+
    geom_point(data = hist.prices.melt[hist.prices.melt$this.date == this.date,],size=3.5)+
    theme_bw() +
    scale_y_continuous("Price")+
    scale_colour_brewer("",palette="Set1")+
    facet_grid(variable~.,scales="free_y")+
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          plot.margin = unit(c(0.53,0.2,0.75,0), "cm"))+
    guides(colour=FALSE)
  
  library(gtable)
  #Combine the two ggplots using gtable into a single graphic.
  #thanks to numerous stack-overflow posters
  #for your help with learning about gtable!
  gA <- ggplot_gtable(ggplot_build(timeseries.plot))
  gB <- ggplot_gtable(ggplot_build(correlplot))
  maxWidth = grid::unit.pmax(gA$widths[2:3], gB$widths[2:3])
  gA$widths[2:3] <- as.list(maxWidth)
  gB$widths[2:3] <- as.list(maxWidth)
  
  gt <- gtable(widths = unit(c(0.8, 1), "null"), height = unit(1, "null"))
  gt <- gtable_add_grob(gt, gA, 1, 1)
  gt <- gtable_add_grob(gt, gB, 1, 2)
  
  #save the graphic as png
  #this will form a single frame in the video
  base.graphpath = "~\\FilmStills_"
  graphpath = paste(base.graphpath,j,".png",sep="")
  
  png(graphpath,
      width = 1500,
      height = 750,
      res = 150
      )
    grid.newpage()
    grid.draw(gt)
  dev.off()

  return (graphpath)
}

#produce a graphic for each date in the date range
FilmStillPaths = lapply(1:length(dates.to.do),ProduceFilmStill)

#use the ffmpeg program to combine the pngs to make an mpeg.
#ffmpeg -y -r 50  -s 1500X750 -f image2 -i "FilmStills_%d.png" FinancialCrisisOnTape.mpg