US Unemployment

####################################################################
### Graph US Unemployment data from FRED
###    - <Mathew.Binkley@Vanderbilt.edu>
####################################################################

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

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

### Set my FRED API key to access the FRED database.
### You may request an API key at:
### https://research.stlouisfed.org/useraccount/apikeys
api_key_fred <- "GO_GET_YOUR_OWN"
fredr_set_key(api_key_fred)

### We only want to plot the last ~3 years of data
date_start <- as.Date("2017-01-01")
date_end   <- as.Date(now())


####################################################################
### Fetch unemployment data from FRED
####################################################################
data   <- fredr(series_id = "UNRATENSA", frequency = "m") %>%
  as_tsibble(index = "date")

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


####################################################################
### Deseasonalize using X13-ARIMA-SEATS
####################################################################
values_ts <- values %>% ts(frequency = 12,
                           start = c(year(min(date)), month(min(date))))

x13arima <- seas(values_ts, estimate.maxiter = 10000, x11 = "")
x13arima_fit <- x13arima$data %>% as_tibble()
x13arima_trend <- x13arima_fit$final
x13arima_seasonal <- x13arima_fit$seasonal

data_df <- data.frame(date     = date,
                      values   = values,
                      trend    = x13arima_trend,
                      seasonal = x13arima_seasonal)

data_subset <- data_df %>% filter(date >= date_start & date <= date_end)

####################################################################
### Graph:  Deseasonalized unemployment compared to raw data
####################################################################
c1 <- "U.S. Bureau of Labor Statistics, Unemployment Rate in US [UNRATENSA]\n"
c2 <- "Retrieved from FRED, Federal Reserve Bank of St. Louis;\n"
c3 <- "https://fred.stlouisfed.org/series/UNRATENSA\n"
c4 <- paste("Data retrieved on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))
caption  <- paste(c1, c2, c3, c4)

title    <- paste("US Unemployment Rate\nbased on unemployment data from",
                  month.name[month(last(date))], year(last(date)))

s1 <- "Raw unemployment rate in grey, "
s2 <- "deseasonalized unemployment rate in gold\n"
s3 <- "Seasonal signal in dotted blue\n"

subtitle <- paste(s1, s2, s3, sep = "")
xlab     <- "Year"
ylab     <- "Unemployment Rate"
caption  <- paste(c1, c2, c3, c4)

p <- ggplot(data = data_subset, aes(x = date, y = values)) +

  theme_economist() +

  geom_point(color = "grey60", shape = 1) +

  geom_line(data = data_subset, aes(x = date, y = values),
            size = 1.3, color = "black", alpha = 0.1) +

  geom_line(data = data_subset, aes(x = date, y = trend),
            size = 1.3, color = "goldenrod2", alpha = 0.8) +

  geom_line(data = data_subset,
            aes(x = date,
                y = (max(trend) + min(trend)) / 1.5 + seasonal),
            size = 1.3, color = "steelblue2",
            alpha = 0.5, linetype = "dotted") +

  labs(title = title, subtitle = subtitle, caption = caption,
       x = xlab, y = ylab)
print(p)