Sahm Rule (R Impementation)

### Sahm's rule in R - by <Mathew.Binkley@Vanderbilt.edu>
### based on Claudia Sahm's paper at:
### https://www.hamiltonproject.org/assets/files/Sahm_web_20190506.pdf
###
### Sahm's law is a lagging indicator of a recession.  It is based on the
### difference between the three month moving average of unemployment and
### the lowest unemployment over the last 12 months.  If the difference
### is > 0.5 (at the Federal level), then the US is or is about to be in
### recession.

### Set your FRED API key here.  You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
api_key_fred <- "WHATEVER_YOUR_FRED_API_KEY_IS"

### Specifiy the start/end date of the graph.  1948-01-01 is the
### earliest possible time for UNRATE
date_start <- "1948-01-01" %>% as.Date()
date_end   <- Sys.Date()   %>% as.Date()

####################################################################
### Load necessary R packages and set the FRED API key
####################################################################
packages <- c("fredr", "lubridate", "tidyverse", "tsibble", "TTR")

### Install packages if needed, then load them quietly
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
invisible(lapply(packages, "library", quietly = TRUE,
                 character.only = TRUE, warn.conflicts = FALSE))

fredr_set_key(api_key_fred)

####################################################################
### Add recession bars to ggplot graphs
####################################################################
geom_recession_bars <- function (date_start, date_end) {

  date_start <- date_start %>% as.Date(origin = "1970-01-01")
  date_end   <- date_end   %>% as.Date(origin = "1970-01-01")

  recessions_tibble <- tibble(

    peak   = c("1857-06-01", "1860-10-01", "1865-04-01", "1869-06-01",
               "1873-10-01", "1882-03-01", "1887-03-01", "1890-07-01",
               "1893-01-01", "1895-12-01", "1899-06-01", "1902-09-01",
               "1907-05-01", "1910-01-01", "1913-01-01", "1918-08-01",
               "1920-01-01", "1923-05-01", "1926-10-01", "1929-08-01",
               "1937-05-01", "1945-02-01", "1948-11-01", "1953-07-01",
               "1957-08-01", "1960-04-01", "1969-12-01", "1973-11-01",
               "1980-01-01", "1981-07-01", "1990-07-01", "2001-03-01",
               "2007-12-01") %>% as.Date(),

    trough = c("1858-12-01", "1861-06-01", "1867-12-01", "1870-12-01",
               "1879-03-01", "1885-05-01", "1888-04-01", "1891-05-01",
               "1894-06-01", "1897-06-01", "1900-12-01", "1904-08-01",
               "1908-06-01", "1912-01-01", "1914-12-01", "1919-03-01",
               "1921-07-01", "1924-07-01", "1927-11-01", "1933-03-01",
               "1938-06-01", "1945-10-01", "1949-10-01", "1954-05-01",
               "1958-04-01", "1961-02-01", "1970-11-01", "1975-03-01",
               "1980-07-01", "1982-11-01", "1991-03-01", "2001-11-01",
               "2009-06-01") %>% as.Date()
  )

  recessions_trim <- recessions_tibble %>%
    filter(peak >= min(date_start) & trough <= max(date_end))

  if (nrow(recessions_trim) > 0) {
      recession_bars = geom_rect(data        = recessions_trim,
                                 inherit.aes = F,
                                 fill        = "darkgray",
                                 alpha       = 0.25,
                                 aes(xmin = as.Date(peak,   origin="1970-01-01"),
                                     xmax = as.Date(trough, origin="1970-01-01"),
                                     ymin = -Inf, ymax = +Inf))
  } else {
    recession_bars = geom_blank()
  }
}

##########################################################################
### Fetch unemployment data from FRED and calculate Sahm Rule
##########################################################################
data <- fredr(series_id = "UNRATE", frequency = "m") %>%
        as_tsibble(index = "date")

date         <- data %>% pull("date") %>% as.Date()
unemployment <- data %>% pull("value")

### Compute the moving average of unemployment over past 3 months
unemployment_ma <- unemployment %>% SMA(n = 3)

### Compute the minimum unemployment over the past 12 months
unemployment_min <- vector(length = length(unemployment))
for (i in (length(unemployment):(1 + 12))) {
  unemployment_min[i] <- min(unemployment[(i - 13):(i)])
}

### Compute the difference between the 3 month average
### and the 12 month minimum unemployment
sahm_rule <- vector(length = (length(unemployment)))
for (i in seq_len(length(unemployment))) {
  sahm_rule[i] <- unemployment_ma[i] - unemployment_min[i]
}

### Strip off the first 14 months of data since it's junk
### 2 months come from the 3-month moving average
### 12 months come from the 12-month minimum window
sahm_rule <- sahm_rule[(14):(length(sahm_rule))]
date_diff <- date[(14):(length(date))]

sahm_rule_df <- data.frame(date = date_diff, values = sahm_rule) %>%
                filter(date >= date_start & date <= date_end)

###################################################################
### Graph 1:  Search results with trend and seasonality
###################################################################
c1 <- "U.S. Bureau of Labor Statistics, Unemployment Rate [UNRATE]\n"
c2 <- "retrieved from FRED, Federal Reserve Bank of St. Louis\n"
c3 <- "https://fred.stlouisfed.org/series/UNRATE\n"
c4 <- paste("Data retrieved on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))

title    <- "Sahm Rule for the United States"
subtitle <- "Recessions marked with vertical bars"
xlab     <- "Year"
ylab     <- "Percent"
caption  <- paste(c1, c2, c3, c4)

p <- ggplot(sahm_rule_df, aes(x = date, y = values)) +
     theme_bw() +
     theme(legend.position = "none") +
     geom_line(data = sahm_rule_df, size = 1.3, color = "darkblue",
               aes(y = sahm_rule, color = "Sahm Rule")) +
     geom_recession_bars(min(sahm_rule_df$date), max(sahm_rule_df$date)) +
     geom_hline(yintercept = 0.5, size = 1, linetype = "dotted",
                color = "darkred", alpha = 0.5) +
     labs(title = title, subtitle = subtitle, caption = caption,
          x = xlab, y = ylab)
print(p)