SpS5EBV9

1### Create a map showing trend unemployment velocity - Mathew Binkley
2
3### You need a FRED API key in order to pull data from FRED.
4### You may request an API key at:
5### https://research.stlouisfed.org/useraccount/apikeys
6api_key_fred   <- "INSERT_YOUR_FRED_API_KEY_HERE"
7
8### We need the following packages for this example.
9packages <- c("tidyverse", "lubridate", "fredr", "ggplot2", "maps", "sf",
10              "ggthemes", "tsibble", "dplyr", "broom", "ggfortify", "forecast")
11
12### Install packages if needed, then load them quietly
13new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
14if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
15invisible(lapply(packages, "library",
16                 quietly = TRUE,
17                 character.only = TRUE,
18                 warn.conflicts = FALSE))
19
20### Now set the FRED API key
21fredr_set_key(api_key_fred)
22
23### We need several years worth of data in order to filter out
24### the seasonal component.  Ten years is a nice round number...
25date_start <- as.Date("2010-01-01")
26date_end   <- as.Date(now())
27
28trend_velocity <- tibble(state = "", velocity = "")
29
30for (state in state.abb) {
31
32  data <- fredr(series_id = paste(state, "URN", sep = ""),
33                observation_start = date_start,
34                observation_end   = date_end,
35                frequency = "m")
36
37  # Extract the trend
38  date   <- data %>% as_tsibble(index = "date") %>% pull("date")
39  values <- data %>% as_tsibble(index = "date") %>% pull("value")
40
41  # Decompose the data and pluck out the trend component
42  trend <- values %>%
43           ts(frequency = 12, start = c(year(min(date)), month(min(date)))) %>%
44           mstl() %>%
45           trendcycle()
46
47  # Take a 6-month average of unemployment velocity.
48  velocity <- trend %>% diff() %>% tail(n = 6) %>% mean()
49
50  trend_velocity <- trend_velocity %>%
51                    add_row(state = state, velocity = velocity)
52
53}
54
55### use the "state.name" and "state.abb" databases to convert two-letter state
56### codes to lowercase names (tennessee, california, etc).   Add the resulting
57### lowercase name to the tibble
58trend_velocity$state_name <- tolower(state.name[match(trend_velocity$state, state.abb)])
59
60### Filter out an empty value
61trend_velocity <- trend_velocity %>% filter(!(state == ""))
62
63### Load a map of states in the USA.   The maps have a list of lowercase
64### state names, which will use to match against "pres" down below
65us_map <- maps::map("state", plot = FALSE, fill = TRUE)
66
67### Change the latitude/longitude data to a simple feature object
68us_map <- sf::st_as_sf(us_map)
69
70### Change the name of the "ID" column to "state_name"
71names(us_map) <- c("geometry", "state_name")
72
73### Remove the District of Colombia from our map
74us_map <- us_map %>% filter(state_name != "district of columbia")
75
76### Add our velocity data to the map data
77us_map <- us_map %>% left_join(trend_velocity, by = "state_name")
78
79ggplot(us_map, aes(fill = as.numeric(velocity) <= 0), col = "black") +
80   geom_sf(aes(alpha = abs(as.numeric(velocity)))) +
81   coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0) +
82   scale_fill_manual(values = c("TRUE" = "darkgreen", "FALSE" = "red")) +
83   scale_alpha(range = c(0.1, 1)) +
84   theme_void() +
85   theme(legend.position = "none")