Recession Indicator using d(Trend Unemployment)/dt

### Leading recession indicator using unemployment.   The algorithm is a close
### cousin of the Sahm Rule, which I derived after experimenting with her method.
###
### The Sahm Rule is a lagging indicator, while this offers a *leading*
### indicator of recessions.   It works perfectly across the entire 70 year
### UNRATENSA dataset, with no false positives or negatives.   So given its
### accuracy and its similar level of simplicity to the Sahm Rule, I believe
### that the rule is doing something very useful.   And I have not found anyone
### using this particular technique, so I believe the technique may be novel.
###
### What my rule does (and how it differs from the Sahm rule) are:
###
### * I use UNRATENSA unemployment data from FRED
###   (Sahm uses seasonally-adjusted UNRATE instead)
###
### * I use mstl() to extract the trend component of unemployment
###   (Sahm uses a moving average which should work somewhat similarly
###    on deseasonalized data)
###
### * I graph d(trend unemployment)/dt
###   (Sahm graphs the difference between the 3-month average and
###    the 12 month minimum)
###
### * If d(trend)/dt is above a small threshold, a recession is
###   extremely likely.  (Sahm requires her rule to exceed a
###   similar threshold of 0.5 before calling a recession).
###
### For this indicator, the important feature is *not* the peak, but rather the
### point where unemployment velocity turns positive.  When the velocity rises
### above the threshold, it is a reliable sign that a recession is coming a few
### months down the road.
###
### written by /u/MetricT (PM me for my actual email)

### 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"

####################################################################
### Load necessary R packages, set the FRED API key,  and set
### start/end dates for graph
####################################################################

### We need the following packages for this example.
packages <- c("fredr", "lubridate", "fredr", "ggplot2", "forecast",
              "ggthemes", "tsibble", "dplyr", "magrittr", "broom", "scales")

### 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)

### 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   <- now()        %>% as.Date()

####################################################################
### 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 data from FRED and compute unemployment rate
####################################################################
unemploy_data  <- fredr(series_id = "UNEMPLOY",  frequency = "m") %>% as_tsibble(index = "date")
clf160v_data   <- fredr(series_id = "CLF16OV",   frequency = "m") %>% as_tsibble(index = "date")

date     <- unemploy_data %>% pull("date")

unemploy_trend <- unemploy_data %>% as_tsibble(index = "date") %>% pull("value") %>%
                  ts(frequency = 12, start = c(year(min(date)), month(min(date)))) %>%
                  mstl() %>% trendcycle()

clf160v_trend  <- clf160v_data %>% as_tsibble(index = "date") %>% pull("value") %>%
                  ts(frequency = 12, start = c(year(min(date)), month(min(date)))) %>%
                  mstl() %>% trendcycle()

unemployment   <- 100 * unemploy_trend / clf160v_trend

####################################################################
### Calculations to derive unemployment "velocity", d(unemployment)/dt
####################################################################
diff_trend <- unemployment %>% diff()
diff_date  <- date %>% tail(n = length(diff_trend))

### Use Friedman's SuperSmoother to eliminate "stairstepping" in the
### derivative.   Just comment out the two lines below if you want to
### use raw data instead of the smoothed data.   It works the same and
### arrives at the same conclusions, but the unsmoothed data yields
### uglier graphs.
tmp <- supsmu(decimal_date(diff_date), diff_trend, span = 3/length(diff_trend))
diff_trend <- tmp$y

velocity_df <- data.frame(diff_date, diff_trend) %>%
               filter(diff_date >= date_start &
                      diff_date <= date_end)

####################################################################
### Find locations where the unemployment velocity line grows
### through x = 0
####################################################################
axis_crossing_locations <- function(date, value) {

  # The threshold necessary to eliminate the two false positives from the 1960's
  threshold <- 0.0105665

  # Shift the value by the threshold
  adj_value <- value - threshold

  # See if the sign changes before/after, which indicates it has
  # passed through x = 0
  sign     <- value     %>% sign()
  adj_sign <- adj_value %>% sign()

  # Computer time derivative of sign changes
  diff_sign     <- sign     %>% diff()
  diff_adj_sign <- adj_sign %>% diff()

  # Find dates when rising unemployment velocity exceeds the
  # threshold.   We want to add one month to the date this
  # gives, thus the "%m+% months(1)" bit.
  locs <- date[diff_adj_sign == 2] %m+% months(1)
  locs
}

axiscrossings <- axis_crossing_locations(diff_date, diff_trend) %>% as.Date()

####################################################################
### Graph:  Trend Unemployment Velocity
####################################################################
c1 <- "U.S. Bureau of Labor Statistics, Unemployment Level [UNEMPLOY]\n"
c2 <- "U.S. Bureau of Labor Statistics, Civilian Labor Force Level [CLF16OV]\n"
c3 <- "Retrieved from FRED, Federal Reserve Bank of St. Louis;\n"
c4 <- paste("Data retrieved on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))

s1 <- "Indicator (vertical dotted line) shows imminent recession when velocity rises above 0\n"
s2 <- "The farther above/below 0, the faster unemployment is growing/falling\n"
s3 <- "Recessions marked with vertical bars\n"

title    <- "Trend Unemployment Rate Velocity"
subtitle <- paste(s1, s2, s3, sep = "")
xlab     <- "Year"
ylab     <- "Percent/Month"
caption  <- paste(c1, c2, c3, c4, sep = "")

p <- ggplot(velocity_df, aes(x = diff_date, y = diff_trend)) +

  ### Plot unemployment velocity
  geom_line(size = 1.3, color = "red3") +

  ### Add a line at velocity = 0.  Above this line, unemployment is rising, and
  ### the farther above the line, the faster it rises.   Similarly, below the
  ### line unemployment is falling, and the further below the line, the faster
  ### it falls
  geom_hline(yintercept = 0, size = 1) +

  ### Draw the recession indicators when the unemployment velocity has
  ### has grown beyond x = 0
  geom_vline(xintercept = axiscrossings, linetype = "dashed") +

  # Add recession bars
  geom_recession_bars(min(velocity_df$diff_date), max(velocity_df$diff_date)) +

  ### Misc graph stuff
  theme_economist() +
  scale_x_date(breaks = pretty_breaks(10), limits = c(date_start, date_end)) +
  labs(title = title, subtitle = subtitle, caption = caption,
       x = xlab, y = ylab)

print(p)