Unemployment with "Doppler map"

### This script maps trend unemployment by state, with a "Doppler map" that
### shows whether unemployment is falling or rising.   In both instances, it
### uses a 6 month average of the data.
###
### 1.  Fetch unemployment data from each state from FRED
### 2.  Use mstl() to decompose data into trend/seasonal/remainder
### 3.  Measure the 6 month average of trend unemployment and d(trend)/dt
### 4.  Turn data into pretty map

### You need a FRED API key in order to pull data from FRED.
### You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
api_key_fred   <- "INSERT_YOUR_FRED_API_KEY_HERE"

### We need the following packages for this example.
packages <- c("tidyverse", "lubridate", "fredr", "ggplot2", "maps", "sf",
              "ggthemes", "tsibble", "dplyr", "broom", "ggfortify",
              "usmap", "forecast", "gridExtra", "cartogram", "RColorBrewer")

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

### Now set the FRED API key
fredr_set_key(api_key_fred)

### We need several years worth of data in order to filter out
### the seasonal component.  Ten years is a nice round number...
date_start <- as.Date("2010-01-01")
date_end   <- as.Date(now())

### Create an empty tibble to hold our data
trend_unemployment <- tibble(state        = character(),
                             unemployment = numeric(),
                             velocity     = numeric())

for (state in state.abb) {

  # Pull data from FRED
  data <- fredr(series_id = paste(state, "URN", sep = ""),
                observation_start = date_start,
                observation_end   = date_end,
                frequency = "m")

  date   <- data %>% as_tsibble(index = "date") %>% pull("date")
  values <- data %>% as_tsibble(index = "date") %>% pull("value")

  # Decompose the data and pluck out the trend component
  trend <- values %>%
           ts(frequency = 12, start = c(year(min(date)), month(min(date)))) %>%
           mstl() %>%
           trendcycle()

  # Compute the 6 month average of unemployment
  unemployment <- trend %>% tail(n = 6) %>% mean()

  # Compute the 6 month average of unemployment velocity
  velocity <- trend %>% diff() %>% tail(n = 6) %>% mean()

  # Append result to our tibble
  trend_unemployment <- trend_unemployment %>%
                        add_row(state, unemployment, velocity)

}

### use the "state.name" and "state.abb" databases to convert two-letter state
### codes to lowercase names (tennnessee, california, etc).   Add the resulting
### lowercase name to the tibble
match <- match(trend_unemployment$state, state.abb)
trend_unemployment$state_name <- state.name[match] %>% tolower()

### Load a map of states in the USA.   The maps have a list of lowercase
### state names, which will use to match against "pres" down below
us_map <- maps::map("state", plot = FALSE, fill = TRUE)

### Change the latitude/longitude data to a simple feature object
us_map <- sf::st_as_sf(us_map)

### Change the name of the "ID" column to "state_name"
names(us_map) <- c("geometry", "state_name")

### Remove the District of Colombia from our map
us_map <- us_map %>% filter(state_name != "district of columbia")

### Add our velocity data to the map data
us_map <- us_map %>% left_join(trend_unemployment, by = "state_name")

### Map the unemployment data
p1 <- ggplot(us_map, aes(fill = unemployment), col = "black") +

      theme_void() +
      theme(legend.position = "bottom") +

      # Title of our graph
      labs(title = "Trend Unemployment Rate by State, 6 month average") +

      # Draw our map
      geom_sf() +

      # Change the map projection to make it look nicer
      coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0) +

      # Use "heat" colors to color the map by unemployment
      scale_fill_gradientn(name = "Unemployment (%)",
                           colours = rev(brewer.pal(9, "Blues")),
                           trans = "log10")

p2 <- ggplot(us_map, col = "black") +

      theme_void() +
      theme(legend.position = "bottom") +

      # Title of our graph
      labs(title = "Trend Unemployment Rate Change by State, 6 month average") +

      # Color states using green for falling unemployment, red for rising
      # unemployment
      geom_sf(aes(fill = velocity)) +
      scale_fill_gradient2(name = "Unemployment (%)",
                           low  = "palegreen4",
                           mid  = "white",
                           high = "firebrick2") +

      # Change the map projection to make it look nicer
      coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0)

### Combine the two graphs
grid.arrange(p1, p2, nrow = 2)