Discussion 4

library(tidyverse)
library(fpp3)
library(fredr)
library(ggthemes)

params <- list(
  series_id = c("UNRATENSA", "ICNSA"),
  frequency = c("m", "m")
)

unemp_df <- pmap_dfr(
  .l = params,
  .f = ~ fredr(series_id = .x, frequency = .y)
) |> 
  select(date:value) |> 
  filter(date >= "1970-01-01") |> 
  pivot_wider(names_from = "series_id", values_from = "value") |> 
  mutate(date = tsibble::yearmonth(date)) |> 
  as_tsibble()

add_rec_shade <- function(st_date, ed_date, shade_color = "darkgray") {
  library(fredr)
  library(ecm)
  library(ggplot2)
  # fredr_set_key("insert_your_api_key")

  recession <- fredr(series_id = "USRECD", observation_start = as.Date(st_date), observation_end = as.Date(ed_date))

  recession$diff <- recession$value - lagpad(recession$value, k = 1)
  recession <- recession[!is.na(recession$diff), ]
  recession.start <- recession[recession$diff == 1, ]$date
  recession.end <- recession[recession$diff == (-1), ]$date

  if (length(recession.start) > length(recession.end)) {
    recession.end <- c(recession.end, Sys.Date())
  }
  if (length(recession.end) > length(recession.start)) {
    recession.start <- c(min(recession$date), recession.start)
  }

  recs <- as.data.frame(cbind(recession.start, recession.end))
  recs$recession.start <- as.Date(as.numeric(recs$recession.start), origin = as.Date("1970-01-01"))
  recs$recession.end <- as.Date(recs$recession.end, origin = as.Date("1970-01-01"))
  if (nrow(recs) > 0) {
    rec_shade <- geom_rect(
      data = recs, inherit.aes = F,
      aes(xmin = recession.start, xmax = recession.end, ymin = -Inf, ymax = +Inf),
      fill = shade_color, alpha = 0.5
    )
    return(rec_shade)
  }
}

Unemployment chart

unemp_df |> 
  pivot_longer(where(is.numeric), names_to = "series_id") |> 
  filter(series_id == "UNRATENSA") |> 
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("5 years"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(2),
    labels = scales::label_percent(scale = 1)
  ) +
  add_rec_shade("1970-01-01", "2023-10-01") +
  labs(
    title = "Seasonally unadjusted unemployment rate",
    subtitle = "Since 1970, shading is recession indicated by NBER"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  theme_calc() +
    theme(axis.title = element_blank())

Initial claims chart

unemp_df |> 
  pivot_longer(where(is.numeric), names_to = "series_id") |> 
  filter(series_id == "ICNSA") |> 
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("5 years"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(500000),
    labels = scales::label_comma(scale = 1)
  ) +
  add_rec_shade("1970-01-01", "2023-10-01") +
  labs(
    title = "Seasonally unadjusted initial unemployment claims",
    subtitle = "Since 1970, shading is recession indicated by NBER"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  theme_calc() +
    theme(axis.title = element_blank())

Initialize model

train <- unemp_df |> 
  filter_index("1970 Jan" ~ "2013 Feb")
test <- unemp_df |> 
  filter_index("2013 Mar" ~ .)

Unemployment rate (NSA) as a function of initial claims

fit <- train |> 
  model(
    arimax = ARIMA(UNRATENSA ~ ICNSA, stepwise = FALSE, approximation = FALSE),
    seasonal_arimax = ARIMA(UNRATENSA ~ ICNSA + season(period = 12)),
  )

fit |> report()

# A tibble: 2 × 8
  .model          sigma2 log_lik   AIC  AICc   BIC ar_roots   ma_roots  
  <chr>            <dbl>   <dbl> <dbl> <dbl> <dbl> <list>     <list>    
1 arimax          0.0369    113. -212. -211. -182. <cpl [3]>  <cpl [13]>
2 seasonal_arimax 0.0385    116. -193. -192. -113. <cpl [12]> <cpl [16]>

fit |> 
  select(seasonal_arimax) |> 
  gg_tsresiduals()

Add new data

future_claims <- new_data(train, 12) |> 
  mutate(ICNSA = 340000)

forecast(fit, future_claims) |> 
  autoplot() +
  autolayer(unemp_df |> 
              filter(date >= yearmonth("2010 Jan"),
                     date <= yearmonth("2015 Jan")),
            show.legend = FALSE) +
  theme_calc() +
      theme(axis.title = element_blank())

model_fc <- forecast(fit, future_claims)

Accuracy metrics

accuracy(model_fc, unemp_df) |> 
  arrange(RMSE) |> 
  select(-.type) |> 
  gt::gt() |> 
  gt::cols_label(
    .model = "Model"
  )

Model	ME	RMSE	MAE	MPE	MAPE	MASE	RMSSE	ACF1
arimax	-0.3669498	0.4300769	0.3669498	-5.238117	5.238117	0.4499070	0.3799473	0.5772289
seasonal_arimax	-0.4851459	0.5584871	0.4851459	-6.924764	6.924764	0.5948239	0.4933900	0.6381919