Plot & Transform

Source Data

We use federal unemployment rate data located at FRED - https://fred.stlouisfed.org/series/UNRATE.

This is US Monthly Unemployment Data.

library(ggplot2)
library(tsibble)
library(tidyverse)
library(lubridate)
library(forecast)
library(tseries)
library(fabletools)
library(fable)
library(feasts)
library(kableExtra)
library(dplyr)


Import data

raw_data_path<-("https://raw.githubusercontent.com/amedina613/Data624-ARIMA-group/refs/heads/main/UNRATE.csv")

unrate_raw <- read.csv(raw_data_path)
head(unrate_raw)
##         DATE UNRATE
## 1 1948-01-01    3.4
## 2 1948-02-01    3.8
## 3 1948-03-01    4.0
## 4 1948-04-01    3.9
## 5 1948-05-01    3.5
## 6 1948-06-01    3.6
# Convert the date column to date format
#unrate_raw$DATE <- as.Date(unrate_raw$DATE, format="%Y-%m-%d")
# Convert data to tsibble
#unrate_tsibble <- unrate_raw |>
#  as_tsibble(index = DATE)

#head(unrate_tsibble)

unrate_raw$monthFormat <- yearmonth(unrate_raw$DATE)
unrate_tsibble <- unrate_raw |> select(-DATE) |> as_tsibble(index = monthFormat)
autoplot(unrate_tsibble, UNRATE) + 
  labs(title = "Unemployment Rate Time Series",
       x = "Month", y = "Unemployment Rate")+
  theme_minimal()

This graph clearly shows noticeable fluctuations and peaks.


Stabilize variance with Box-Cox

#find the optimal lambda value using the guerrero feature
lambda <- unrate_tsibble %>%
  features(UNRATE, features = guerrero) %>%
  pull(lambda_guerrero)

print(lambda)
## [1] -0.8999268
# apply box-cox transformation
unrate_transformed <- unrate_tsibble %>%
  mutate(UNRATE_trans = box_cox(UNRATE, lambda))

# plot the transformed series using autoplot
autoplot(unrate_transformed, UNRATE_trans) + 
  labs(title = "Box-Cox Transformed Unemployment Rate",
       x = "Month", y = "Transformed Unemployment Rate")+
  theme_minimal()


Test for Stationarity

Augmented Dickey-Fuller (ADF) test
# perform ADF test
adf_test <- adf.test(unrate_transformed$UNRATE_trans)
print(adf_test)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  unrate_transformed$UNRATE_trans
## Dickey-Fuller = -4.4331, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

The results of the ADF test show a p-value of 0.01, which is less than 0.05. This could indicate that it’s stationary. For further inspection, we can inspect the autocorrelation function


Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test

Sometimes, ADF and KPSS tests can show differing results.

We conduct a KPSS test as well.

unrate_tsibble |> features(UNRATE, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1      1.02        0.01

We get a p value of 0.01, indicating that we reject the null hypothesis that the data is stationary.

We conduct a KPSS test on differenced data:

unrate_tsibble |> mutate(diff_unrate =difference(UNRATE)) |>
  features(diff_unrate, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0405         0.1

The p-value of 0.1 indicates that the data appears stationary after one order of differencing.



Fit & Check

For the ARIMA model below, we are using the original unemployment rate values (unrate_tsibble) to demonstrate the non seasonal ARIMA model.


Plot the ACF & PACF

We plot the ACF and PACF to determine appropriate values for the p (auto-regressive) and q (moving average) parameters of an ARIMA model after differencing the time series to achieve stationarity. The differenced residuals are used to evaluate whether the data still exhibits patterns that need modeling or if it resembles white noise.

Both the ACF and PACF show minimal significant spikes, indicating a lack of strong correlations at different lags. This suggests that the appropriate model could be ARIMA(0,1,0), as there are no significant auto regressive or moving average components required after first-order differencing (d = 1). The residuals indicate that the time series is adequately modeled without additional AR or MA terms. The ARIMA(0,1,0) is equivalent to a random walk.

unrate_tsibble |> select(UNRATE) |>
gg_tsdisplay(difference(UNRATE), plot_type = "partial")


Identify Model Candidates

We fit multiple ARIMA models to the UNRATE time series to determine the best model for forecasting using fable.

  • The arima010 model fits an ARIMA(0, 1, 0) model (based on a manual check).
  • The step-wise model automatically selects the best ARIMA parameters using a step-wise algorithm.
  • The search model performs an search to find the optimal ARIMA configuration.


ARIMA Candidates

The model function checks the manually selected ARIMA(0, 1, 0) model for fit as well surfacing two other well fitting models: ARIMA (1, 1, 1) and seasonal ARIMA (1,1,1)(2,0,0)[12]

unrate_fit <- unrate_tsibble |> 
  model(
    arima010 = ARIMA(UNRATE ~ pdq(0, 1,0)),
    stepwise = ARIMA(UNRATE),
    search = ARIMA(UNRATE, stepwise = FALSE))

# Create table
unrate_fit %>%
  kable("html", caption = "ARIMA Model Candidates") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "responsive"))
ARIMA Model Candidates
arima010 stepwise search
<ARIMA(0,1,0)> <ARIMA(1,1,1)> <ARIMA(1,1,1)(2,0,0)[12]>


Check Model Fit & Select Model

To determine the best model for forecasting, we focus on the AICc metric to select the model. Both Sigma2, which measures variance across the model, and log-likelihood metric, which measures how well model fits the data, are very similar, the AICc points to the best candidate.

In this case, the ARIMA (1, 1, 1) model identified by the stepsister algorithm appears to be the best fit.

# Convert the model summary to a table, dropping the ar_roots and ma_roots columns
unrate_summary <- glance(unrate_fit) %>%
  select(-ar_roots, -ma_roots) %>%
  arrange(AICc)

# Create table
unrate_summary %>%
  kable("html", caption = "ARIMA Model Fits for Unemployment Rate") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "responsive"))
ARIMA Model Fits for Unemployment Rate
.model sigma2 log_lik AIC AICc BIC
stepwise 0.1736121 -498.4630 1002.926 1002.952 1017.396
search 0.1738165 -498.0089 1006.018 1006.084 1030.134
arima010 0.1746028 -502.0661 1006.132 1006.137 1010.955


Check Fit Residuals

We check the residuals of an ARIMA model to ensure that they resemble white noise, indicating that the model has captured all the underlying patterns in the time series and that there is no remaining structure.

The residuals appear mostly random with no clear pattern, and the ACF shows no significant spikes beyond the confidence bounds, suggesting no significant autocorrelation left in the residuals. The histogram also indicates that the residuals are centered around zero. These observations suggest that the selected ARIMA model provides an adequate fit for the data, as the residuals do not exhibit systematic patterns.

The significant spike / outlier in the residual occurs during the 2020 COVID pandemic.

unrate_fit |>
  select(stepwise)|>
  gg_tsresiduals()


Pormanteau Test

The Portmanteau test is used to determine whether the residuals of an ARIMA model are independently distributed, testing for remaining autocorrelation.

In the result shown, the p-value is 0.9849195, which is high. We fail to reject the null hypothesis, indicating that there is no significant autocorrelation present in the residuals. This implies that the ARIMA model has captured all the key patterns in the time series, and the residuals resemble white noise.

portm <- augment(unrate_fit) |>
  filter(.model == "stepwise") |>
  features (.innov, ljung_box, lag = 10, dof = 2)


portm %>%
  kable("html", caption = "ARIMA Model Candidates") %>%
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "responsive"))
ARIMA Model Candidates
.model lb_stat lb_pvalue
stepwise 1.863305 0.9849195


Forecast & Measure

Train & Test the ARIMA(1,1,1) model

We perform a train-test forecast on the ARIMA(1,1,1) model allows to evaluate its effectiveness at capturing the underlying patterns in the time series and providing reliable forecasts.

In case, the model’s predictive power forecasts well of a roughly 2 year period (24 monthly lags) before losing it predictive power and converging at a single value. The ARIMA(1,1,1) works well over a period of about 2 years – but the unemployment rate data is likely better suited for a rolling forecast using incremental intervals of 24 monthly lags.

# Split the data 80/20
n <- nrow(unrate_tsibble)
train_size <- floor(0.8 * n)
unrate_train <- unrate_tsibble[1:train_size, ]
unrate_test <- unrate_tsibble[(train_size + 1):n, ]

# Fit the ARIMA model to the training set
unrate_fit_train <- unrate_train |>
  model(stepwise = ARIMA(UNRATE ~ pdq(1,1,1)))

# Forecast using the stepwise fitted model
forecast_h <- n - train_size
unrate_forecast <- unrate_fit_train |>
  select(stepwise) |>
  forecast(h = forecast_h)

# Convert tsibble for plotting
unrate_train_tbl <- unrate_train |>
  as_tibble() |>
  mutate(type = "Training")

unrate_test_tbl <- unrate_test |>
  as_tibble() |>
  mutate(type = "Actual")

unrate_forecast_tbl <- unrate_forecast |>
  as_tibble() |>
  mutate(type = "Forecast")

# Combine training, actual, and forecast data
combined_data <- bind_rows(
  unrate_train_tbl %>%
    select(monthFormat, UNRATE, type),
  unrate_test_tbl %>%
    select(monthFormat, UNRATE, type),
  unrate_forecast_tbl %>%
    select(monthFormat, UNRATE = .mean, type))

# Plot
ggplot(combined_data, aes(x = monthFormat, y = UNRATE, color = type)) +
  geom_line() +
  labs(title = "ARIMA(1,1,1) Model Forecast on Actual Data",
       x = "Month",
       y = "Unemployment Rate",
       color = "Legend") +
  theme_minimal()


ARIMA Forecast

unrate_fit_train |> select(stepwise) |> forecast(unrate_test) |> autoplot(unrate_test)

While our model did a decent job for the first couple of years of the test set, the model resulted in an unemployment rate that stabilized around 7% while the actual unemployment rate dropped below 5% before seeing an unprecedented spike due to the COVID pandemic.

unrate_fit_train |>  accuracy()
## # A tibble: 1 × 10
##   .model   .type         ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <chr>    <chr>      <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 stepwise Training 0.00527 0.197 0.142 0.0612  2.66 0.164 0.165 -0.146

Wide confidence intervals make the point forecast somewhat less reliable. The real unemployment rate values mostly stayed inside the 80% confidence level except for the COVID spike.