R file for charting nominal GDP versus Treasury rates

GDP vs Treasury Rates

# R plot of nominal GDP vs 10-year rate using data from FRED
# Frozenquant.com
# Thanks to KLR for his posts on charting with xtsExtra

require(quantmod)
require(PerformanceAnalytics)
require(xtsExtra)
require(RColorBrewer)

# Get data for chart --------------------------------------------------------
symbols <- c("GDP","DGS10")
getSymbols(symbols, src='FRED')

# Went with an index class yearqtr from zoo
GDP2 <- xts(coredata(GDP), order.by=as.yearqtr(index(GDP)))
# yoy change
GDP.ch <-100*(GDP2/lag(GDP2, 4)-1)
# to.quarterly uses a yearqtr index class
TY.qtr <- Cl(to.quarterly(DGS10))
colnames(TY.qtr) <- "TY.qtr"
# Merge data series into one xts object
gdp.ty <- merge(GDP.ch, TY.qtr, join="left")["1983-12-31::"]
gdp.ty <- merge(gdp.ty, gdp.ty[,2]-gdp.ty[,1])
colnames(gdp.ty) <- c("GDP.chng", "TY", "Diff")
# Obviously this should be more elegant; get average diff for 
# non-recession periods.  Used for chart horizontal lines.
avg.real.rate <- rep(0,4)
avg.real.rate[1] <- mean((gdp.ty["1984-01/1990-07"])[,3])
avg.real.rate[2] <- mean((gdp.ty["1991-04/2001-03"])[,3])
avg.real.rate[3] <- mean((gdp.ty["2001-11/2007-12"])[,3])
avg.real.rate[4] <- mean((gdp.ty["2009-07/"])[,3])

# ----------------------------------------------------------------
# NBER recession dates
## 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"
)

# This is using KLR's idea of coding buy and sell dates
buydates = as.Date(c("1984-01-31",
             "1991-04-01",
             "2001-11-30",
             "2009-07-01"))
selldates = as.Date(c("1990-07-31",
              "2001-03-30",
              "2007-12-31",
              "2013-12-31"))  
# ----------------------------------------------------------------------------
# Custom panel with no black baseline, 200 day MA or skinny boxplot
custom.panel4 <- function(index,x,...) {
  default.panel(index,x,...)
  # horizontal baseline, lhs axis (side 2)
  abline(h=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=par("yaxp")[3]),col="gray60",lty=3)
  axis(side=2,col="gray60",col.axis="black",lwd=0,lwd.ticks=FALSE,las=1,
       at=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=abs(par("yaxp")[3])),
       labels=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=abs(par("yaxp")[3])))
}


# Custom panel with no 200 day MA with skinny boxplot and black baseline
custom.panel3 <- function(index,x,...) {
  default.panel(index,x,...)
  # horizontal baseline, lhs axis (side 2)
  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")
  axis(side=2,col="gray60",col.axis="black",lwd=0,lwd.ticks=FALSE,las=1,
       at=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=abs(par("yaxp")[3])),
       labels=pretty(c(par("yaxp")[1],par("yaxp")[2]),n=abs(par("yaxp")[3])))
  
  # Skinny box-whisker plot in RHS margin lines are 1.5 x IQR
  # Draw lines rather than segments to get outside of the plot area.
  # par("usr") = gives extremes of the user coordinates of the plotting region
  rt.base <- par("usr")[2] * 0.9975 + par("usr")[1] * 0.0025 # far dim - 1% x (far - near)
  # print(par("usr"))
  lines(rep(rt.base, 2),
        c(max(min(x), quantile(x, 0.25) - 1.5 * (quantile(x, 0.75)-quantile(x, 0.25))),
          min(max(x), quantile(x, 0.75) + 1.5 * (quantile(x, 0.75)-quantile(x, 0.25)))),
        lwd=1) 
  lines(rep(rt.base, 2),
        c(quantile(x, 0.25),
          quantile(x, 0.75)),
        lwd= 4, col="navajowhite4")  # gray25
  # Median point added - could add others
  points(x= rep(rt.base, 1),
         y=quantile(x, c(0.5)), 
         cex= c(1.2), pch=18,
         col= c("darkred")) 
 # Segment Averages - note changed to as.yearqtr from as.Date
 for(i in 1:4) {
   segments(x0=index[which(index(x) == as.yearqtr(buydates[i]))],
            x1=index[which(index(x) == as.yearqtr(selldates[i]))],
            y0=avg.real.rate[i],
            y1=avg.real.rate[i])
 } 
}

# xtsExtra plots --------------------------------------
col.bl <- brewer.pal(9,"Blues")[8]
col.gr <- brewer.pal(9,"Greens")[7]
col.red <- brewer.pal(9,"Reds")[8]
col.comb <- c(col.gr, col.bl, col.red)

plot.xts(gdp.ty,  #limit to Dec 2009 to Current so more easily visible
         col = col.comb,
         lwd = 2,                            #line width; will do 2
         las = 1,                            #do not rotate y axis labels
         screens = c(1,1,2),  
         layout.screens = c(1,1,1,1,2,2),
         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 = "gray94"),
         bty="n",
         auto.grid=FALSE,
         major.format="%Y",    # was "%b %Y"
         major.ticks="years",     
         minor.ticks=FALSE,
         #col.axis="transparent", # This is the thing that causes the invisible axis
         yax.loc="none",
         cex.axis=0.9,
         panel=c(custom.panel4,custom.panel3),
         main = NA)  #will do title later so we have more control

#define title separately so we have more control
title(main = "10-Year Treasury Rate - Nominal GDP",
      outer=TRUE,
      line=-2,
      adj=0.05,
      cex=0.8)
legend("topright",legend=c("Nominal GDP", "10-Yr Tsy"), lty=1, col=c(col.gr, col.bl), 
       bty="n", cex=0.6)