Midterm

Author

Sam Kuhn

Instructions

Find monthly time series data for at least five years.
Use 80% of the data (e.g., the first four years-the training set) to build three models + one ensemble model (average of the three forecasts) for the last year.
You will estimate appropriate model statistics on the with-held 20% test set.
You will compare these statistics across the models and suggest the best.
You will prepare a poster presentation with the attached template.

The poster will include an abstract summarizing the problem, the significance of the problem, the methods, and the results in about 200 words.
An introduction section with a discussion of three to five references associated with the topic (e.g., energy forecasting if that is what you are doing) and the selected methods (e.g., ETS, ARIMA, Regression, Ensemble).
A methods section discussing the chosen methods and formulation.
A results section that illustrates the model performance.
A discussion section that indicates the utility of the models.

You must include 3-5 graphs that provide descriptive and inferential information. The poster presentation format is attached.

Data set

Weekly initial unemployment claims since 1970, not seasonally adjusted

library(tidyverse)
library(fpp3)
library(fredr)
library(ggthemes)
library(ggiraph)
library(fable.prophet)
library(future)
options(scipen = 999)
initial_claims <- fredr(
  series_id = "ICNSA",
  observation_start = as.Date("1967-01-01"),
  observation_end = as.Date("2023-10-01"),
  frequency = "m",
  units = "lin"
) |>
  select(date, value) |>
  mutate(date = tsibble::yearmonth(date)) |>
  as_tsibble()

Plot

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

Major change points:

Height of Volcker-shock - 1981 through 1982
Height of GFC - 2008 through 2010
Peak of COVID pandemic - 2020

initial_claims |>
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  ggrepel::geom_label_repel(aes(date, value, label = scales::comma(value)),
                            show.legend = FALSE,
                            nudge_x = 100,
                            nudge_y = 100,
                            data = initial_claims |> 
                              filter(date == max(date))) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("5 years"),
    labels = scales::label_date("%Y")
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(1000000),
    labels = scales::label_comma()
  ) +
  labs(
    title = "Monthly unemployment insurance claims",
    subtitle = "Seasonally un-adjusted, since 1967",
    caption = "Source: Department of Labor and FRED \n Grey shading represents recessions according to the NBER"
  ) +
  ggthemes::theme_calc() +
  annotate(
    geom = "text",
    x = as.Date("2020-04-01"),
    y = 5000000,
    label = "Peak of COVID-19",
    fontface = "bold"
  ) +
    annotate(
    geom = "text",
    x = 14000,
    y = 1000000,
    label = "Peak of GFC",
    fontface = "bold"
  ) +
    annotate(
    geom = "text",
    x = as.Date("1982-01-01"),
    y = 1200000,
    label = "Volcker-shock",
    fontface = "bold"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  add_rec_shade("1967-01-01", "2023-10-01") +
  theme(axis.title = element_blank())

ggsave("ui_plot.png")

Remove COVID-19 distortion (2021 ~ )

initial_claims |>
  filter_index("2021 Jan" ~ .) |> 
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  ggrepel::geom_label_repel(aes(date, value, label = scales::comma(value)),
                            show.legend = FALSE,
                            nudge_x = 100,
                            nudge_y = 100,
                            data = initial_claims |> 
                              filter(date == max(date))) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("3 months"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(100000),
    labels = scales::label_comma()
  ) +
  labs(
    title = "Monthly unemployment insurance claims",
    subtitle = "Seasonally un-adjusted, since 2021",
    caption = "Source: Department of Labor and FRED \n Grey shading represents recessions according to the NBER"
  ) +
  ggthemes::theme_calc() +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  theme(axis.title = element_blank())

ggsave("covid_ui_plot.png")

Recovery from GFC (2010 ~ 2014)

initial_claims |>
  filter_index("2007 Dec" ~ "2014 Dec") |> 
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  geom_hline(yintercept = 321000, linetype = "dashed", color = "grey50") +
  ggrepel::geom_label_repel(aes(date, value, label = scales::comma(value)),
                            show.legend = FALSE,
                            nudge_x = 100,
                            nudge_y = 100,
                            data = initial_claims |> 
                              filter(date == yearmonth("2014 Dec"))) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("6 months"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(100000),
    labels = scales::label_comma()
  ) +
  labs(
    title = "Monthly unemployment insurance claims",
    subtitle = "Seasonally un-adjusted, recovery period from GFC, 2008-2014",
    caption = "Source: Department of Labor and FRED \n Grey shading represents recessions according to the NBER,
    dashed line is lowest month of UI claims during GFC"
  ) +
  ggthemes::theme_calc() +
  add_rec_shade("2008-01-01", "2015-01-01") +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  theme(axis.title = element_blank())

ggsave("gfc_ui_plot.png")

Volcker-shock 1979-1987

initial_claims |>
  filter_index("1978 Dec" ~ "1987 Dec") |> 
  ggplot(aes(date, value, color = "value")) +
  geom_line(linewidth = .75, show.legend = FALSE) +
  ggrepel::geom_label_repel(aes(date, value, label = scales::comma(value)),
                            show.legend = FALSE,
                            nudge_x = 100,
                            nudge_y = 100,
                            data = initial_claims |> 
                              filter(date == yearmonth("1987 Dec"))) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("1 year"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(100000),
    labels = scales::label_comma()
  ) +
  labs(
    title = "Monthly unemployment insurance claims",
    subtitle = "Seasonally un-adjusted, Volcker-shock, 1979-1987",
    caption = "Source: Department of Labor and FRED \n Grey shading represents recessions according to the NBER"
  ) +
  ggthemes::theme_calc() +
  add_rec_shade("1978-01-01", "1987-12-01") +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  theme(axis.title = element_blank())

ggsave("volcker_shock_ui.png")

STL Decomposition

decomp <- initial_claims |> 
  model(stl = STL(value))
decomp |> 
  components() |> 
  autoplot() +
  scale_x_yearmonth(
    breaks = scales::breaks_width("5 years"),
    labels = scales::label_date_short()
  ) +
  theme_calc()

ggsave("stl_decomp.png")

Build models

ARIMA
Prophet
Neural Net

Create training set

80% of time series ~ 546 observations

train <- initial_claims |> 
  filter_index("1967 Jan" ~ "2012 Jun")

test <- initial_claims |> 
  filter_index("2012 Jul" ~ .)

Warning: This model can take a long time to converge! Enable parallel processing if your compute engine works

model_train <- train |> 
  model(
    auto_arima = ARIMA(value, stepwise = FALSE, approximation = FALSE),
    lm = TSLM(value ~ season() + trend()),
    prophet = prophet(value ~ season("month", 4, type = "multiplicative"))
  ) |> 
  mutate(ensemble = (auto_arima + lm + prophet) / 3)

Report only works for single models

model_train |> 
  select(-ensemble) |> 
  report()

# A tibble: 3 × 19
  .model         sigma2 log_lik    AIC   AICc    BIC ar_roots ma_roots r_squared
  <chr>           <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <list>   <list>       <dbl>
1 auto_arima     6.10e8  -6161. 12335. 12336. 12365. <cpl>    <cpl>       NA    
2 lm             7.18e9  -6964. 12406. 12407. 12466. <NULL>   <NULL>       0.429
3 prophet       NA          NA     NA     NA     NA  <NULL>   <NULL>      NA    
# ℹ 10 more variables: adj_r_squared <dbl>, statistic <dbl>, p_value <dbl>,
#   df <int>, CV <dbl>, deviance <dbl>, df.residual <int>, rank <int>,
#   sigma <dbl>, changepoints <list>

Forecast against test set

model_fc <- model_train |> 
  forecast(test)

model_fc |> 
  autoplot(
    bind_rows(train, test),
    level = NULL
  ) +
  scale_y_continuous(
    labels = scales::label_comma()
  ) +
  labs(
    title = "Forecasted data againt test set",
    subtitle = "Monthly UI claims, not seasonally adjusted",
    caption = "Source: Department of Labor and FRED \n Grey shading represents recessions according to the NBER"
  ) +
  scale_color_manual(
    values = c("#8a100b", "#b29d6c", "#726158", "red")
  ) +
  add_rec_shade("1967-01-01", "2023-10-01") +
  ggthemes::theme_calc() +
  theme(axis.title = element_blank())

ggsave("training_models.png")

Performance metrics - ARIMA performing the best according to RMSE

training_tbl <- accuracy(model_fc, initial_claims) |> 
  arrange(RMSE) |> 
  select(-.type) |> 
  gt::gt() |> 
  gt::cols_label(
    .model = "Model"
  )
gt::gtsave(training_tbl, "training_data_results.png")

Forecast ahead on entire dataset

arima_fc <- initial_claims |> 
  model(auto_arima = ARIMA(value, stepwise = FALSE, approximation = FALSE))

Forecast 12 months ahead

arima_fc |> 
  forecast(h = 12) |> 
  autoplot(level = c(90, 95),
           linewidth = 1.25) +
  autolayer(initial_claims |> 
              filter_index("2021 Jan" ~ .),
            linewidth = 1.25,
            show_gap = FALSE) +
  scale_x_yearmonth(
    breaks = scales::breaks_width("3 months"),
    labels = scales::label_date_short()
  ) +
  scale_y_continuous(
    labels = scales::label_comma()
  ) +
  labs(
    title = "12 months ahead forecast of unemployment insurance claims",
    subtitle = "Monthly UI claims, not seasonally adjusted, ARIMA model",
    caption = "Source: Department of Labor and FRED\nShaded region is 90th and 95th % CI"
  ) +
  scale_color_manual(
    values = c("#8a100b")
  ) +
  ggthemes::theme_calc() +
    theme(axis.title = element_blank())

ggsave("arima_fc.png")