Generate chart with 10y yield oscillator (10y yield minus 260d moving average), +/- 1 stdev dashed lines

10y-oscillator.R

recplot <- function(var, rec.ind, maintitle = "", ylab = "", ylim = NULL)
{
  # Give each recession a separate number, so we can split data by recession
  # and calculate the necessary values to pass to addPolygon() for shading
  r <- rle(as.integer(rec.ind))
  r$values[r$values > 0] <- seq_len(sum(r$values))
  rec <- xts(inverse.rle(r), index(rec.ind))

  # Merge variable and recession data, as a left-join so we do not have
  # observations that only exist in the recession indicator data
  a <- merge(var, rec, join = "left", fill = na.locf)

  # Calculate the y values needed for recession shading
  rng <- matrix(max(abs(var), na.rm = TRUE) * 1.1, 2, 2)
  rng[,2] <- -rng[,2]
  # Subset data to only include recession observations
  reconly <- a[,2] > 0
  # Split data by recession, and create an xts object with x/y pairs
  pl <- lapply(split(a[reconly,], a[reconly,2]),
               function(x) xts(rng, c(start(x), end(x))))

  # Plot data series. Store the plot in an object, so we can build it up
  # without plotting each intermediate plot.
  p <- plot(a[,1], main = maintitle, ylim = ylim, ylab = ylab,
            border = "grey", major.ticks = "years", minor.ticks = "months",
            grid.ticks.on = "years")
  # Add legend first (placement appears to be based on last series added)
  p <- addLegend("topright", c("10y minus its 1y mean", "+/- 1 stdev"),
                 col = "black", lty = c(1, 2))

  # Loop over each recession and add a shaded polygon behind the main series
  for (pr in pl) {
    p <- addPolygon(pr, col = "grey", on = -1)
  }

  # Add standard deviation bands behind main series
  dsd <- xts(rep(sd(var), nrow(a)), index(a))
  p <- addSeries( dsd, on = -1, lty = 2)
  p <- addSeries(-dsd, on = -1, lty = 2)

  # Return the plot object so you could continue to build the plot
  return(p)
}

library(quantmod)
getSymbols("DGS10;USRECD", src="FRED")

var <- na.omit(DGS10 - SMA(na.omit(DGS10), n=260))
recplot(var, USRECD)