aH0nsnMg

1### Render a map of US unemployment by state, and a "Doppler map" showing
2### how quickly unemployment is falling/growing
3### by <Mathew.Binkley@Vanderbilt.edu>
4
5### You need a FRED API key in order to pull data from FRED.
6### You may request an API key at:
7### https://research.stlouisfed.org/useraccount/apikeys
8api_key_fred   <- "INSERT_YOUR_FRED_API_KEY_HERE"
9
10### We need the following packages for this example.
11packages <- c("tidyverse", "lubridate", "fredr", "ggplot2", "maps", "sf",
12              "ggthemes", "tsibble", "dplyr", "broom", "ggfortify",
13              "usmap", "forecast", "gridExtra", "cartogram")
14
15### Install packages if needed, then load them quietly
16new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
17if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
18invisible(lapply(packages, "library",
19                 quietly = TRUE,
20                 character.only = TRUE,
21                 warn.conflicts = FALSE))
22
23### Now set the FRED API key
24fredr_set_key(api_key_fred)
25
26### We need several years worth of data in order to filter out
27### the seasonal component.  Ten years is a nice round number...
28date_start <- as.Date("2010-01-01")
29date_end   <- as.Date(now())
30
31### Create an empty tibble to hold our data
32trend_unemployment <- tibble(state        = character(),
33                             unemployment = numeric(),
34                             velocity     = numeric())
35
36for (state in state.abb) {
37
38  # Pull data from FRED
39  data <- fredr(series_id = paste(state, "URN", sep = ""),
40                observation_start = date_start,
41                observation_end   = date_end,
42                frequency = "m")
43
44  date   <- data %>% as_tsibble(index = "date") %>% pull("date")
45  values <- data %>% as_tsibble(index = "date") %>% pull("value")
46
47  # Decompose the data and pluck out the trend component
48  trend <- values %>%
49           ts(frequency = 12, start = c(year(min(date)), month(min(date)))) %>%
50           mstl() %>%
51           trendcycle()
52
53  # Compute the 6 month average of unemployment
54  unemployment <- trend %>% tail(n = 6) %>% mean()
55
56  # Compute the 6 month average of unemployment velocity
57  velocity <- trend %>% diff() %>% tail(n = 6) %>% mean()
58
59  # Append result to our tibble
60  trend_unemployment <- trend_unemployment %>%
61                        add_row(state, unemployment, velocity)
62
63}
64
65### use the "state.name" and "state.abb" databases to convert two-letter state
66### codes to lowercase names (tennnessee, california, etc).   Add the resulting
67### lowercase name to the tibble
68match <- match(trend_unemployment$state, state.abb)
69trend_unemployment$state_name <- state.name[match] %>% tolower()
70
71### Load a map of states in the USA.   The maps have a list of lowercase
72### state names, which will use to match against "pres" down below
73us_map <- maps::map("state", plot = FALSE, fill = TRUE)
74
75### Change the latitude/longitude data to a simple feature object
76us_map <- sf::st_as_sf(us_map)
77
78### Change the name of the "ID" column to "state_name"
79names(us_map) <- c("geometry", "state_name")
80
81### Remove the District of Colombia from our map
82us_map <- us_map %>% filter(state_name != "district of columbia")
83
84### Add our velocity data to the map data
85us_map <- us_map %>% left_join(trend_unemployment, by = "state_name")
86
87### Map the unemployment data
88p1 <- ggplot(us_map, aes(fill = unemployment), col = "black") +
89
90      # Title of our graph
91      labs(title = "Trend Unemployment by State, 6 month average") +
92
93      # Draw our map
94      geom_sf() +
95
96      # Change the map projection to make it look nicer
97      coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0) +
98  
99      # Use "heat" colors to color the map by unemployment
100      scale_fill_gradientn(name = "Unemployment",
101                           colours = heat.colors(7),
102                           trans = "log10") +
103      theme_void()
104
105p2 <- ggplot(us_map, col = "black") +
106
107      # Title of our graph
108      labs(title = "Trend Unemployment Change by State, 6 month average") +
109
110      # Color states using green for falling unemployment, red for rising
111      # unemployment
112      geom_sf(aes(fill = velocity)) +
113      scale_fill_gradient2(name = "Unemployment",
114                           low  = "palegreen4",
115                           mid  = "white",
116                           high = "firebrick2") +
117
118      # Change the map projection to make it look nicer
119      coord_sf(crs = "+proj=aea +lat_1=25 +lat_2=50 +lon_0=-100", ndiscr = 0) +
120
121      # Disable theming
122      theme_void()
123
124### Combine the two graphs
125grid.arrange(p1, p2, nrow = 2)