### Create a map showing trend unemployment velocity - Mathew Binkley

### 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", "forecast")

### 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_velocity <- tibble(state    = character(),
                         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()

  # Take a 6-month average of unemployment velocity.
  velocity <- trend %>% diff() %>% tail(n = 6) %>% mean()

  # Append the result to our tibble
  trend_velocity <- trend_velocity %>%
                    add_row(state = state, velocity = velocity)

}

### use the "state.name" and "state.abb" databases to convert two-letter state
### codes to lowercase names (tennessee, california, etc).   Add the resulting
### lowercase name to the tibble
match <- match(trend_velocity$state, state.abb)
trend_velocity$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_velocity, by = "state_name")

ggplot(us_map, aes(fill = velocity <= 0), col = "black") +
   geom_sf(aes(alpha = abs(velocity))) +
   coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0) +
   scale_fill_manual(values = c("TRUE" = "darkgreen", "FALSE" = "red")) +
   scale_alpha(range = c(0.1, 1)) +
   theme_void() +
   theme(legend.position = "none")