MBAS 921: Week 6 Tutorial
This tutorial covers ARIMA forecasting.
i. First, we load the required packages and prepare the dataset. The data used in
this tutorial comes from the us_employment dataset. We filter the data so that it
includes observations from the year 2000 onwards and only focuses on the
'Retail Trade' sector.
library(fpp3)
library(readxl)
emp_ts <- us_employment %>%
filter(year(Month) >= 2000, Title == "Retail Trade") |>
select(-Series_ID)
emp_ts
ii. We begin by visualising the time series to observe general patterns such as trend
and seasonality in retail employment.
emp_ts %>%
autoplot(Employed) +
labs(title = "US retail employment time series")
iii. Seasonality can be explored using seasonal plots and subseries plots. These
visualisations help identify recurring patterns within each year.
These plots are very good tools to observe seasonality and trend in the data.
They indeed show both on the same plot, one emphasising on trend and the
other on seasonality.
emp_ts %>%
gg_season(Employed) +
labs(title = "Seasonal plot")
emp_ts %>%
gg_subseries(Employed) +
labs(title = "Subseries plot")
iv. We next examine the autocorrelation and partial autocorrelation functions of the
original series. These plots help identify the presence of correlation between
observations at different lags .
Autocorrelation is a fundamental concept in time series. Understanding it is
therefore important. Plotting ACF for the original series shows whether there are
seasonality or trend in the data.
emp_ts %>%
ACF(Employed) %>%
autoplot() +
labs(title = "ACF of original series")
emp_ts %>%
PACF(Employed) %>%
autoplot() +
labs(title = "PACF of original series")
v. To make the time series stationary, differencing can be applied. Below we
compute the first difference, the seasonal difference with lag 12, and a
combination of both.
emp_ts %>%
mutate(diff_emp = difference(Employed)) %>%
autoplot(diff_emp) +
labs(title = "First differenced series",
y = "First difference")
emp_ts %>%
mutate(season_diff = difference(Employed, lag = 12)) %>%
autoplot(season_diff) +
labs(title = "Seasonally differenced series",
y = "Seasonal difference")
emp_ts %>%
mutate(double_diff = difference(difference(Employed), lag = 12)) %>%
autoplot(double_diff) +
labs(title = "First + seasonal differenced series",
y = "Double difference")
vi. After differencing, we re-examine the ACF and PACF plots. These help guide the
selection of ARIMA model parameters (ACF helps identify MA terms, and PACF
help identify AR terms), although we can go with automatic parameter selection.
We need to do the differencing to the point that differenced data show no
autocorrelation, i.e., all values are within the dashed blue lines.
emp_ts %>%
mutate(double_diff = difference(difference(Employed), lag = 12)) %>%
ACF(double_diff) %>%
autoplot() +
labs(title = "ACF of differenced series")
emp_ts %>%
mutate(double_diff = difference(difference(Employed), lag = 12)) %>%
PACF(double_diff) %>%
autoplot() +
labs(title = "PACF of differenced series")
vii. We divide the dataset into training and test sets. The training set contains
observations up to 2016, while the remaining data is used to evaluate forecasting
performance.
emp_train <- emp_ts %>%
filter(year(Month) <= 2016)
emp_test <- emp_ts %>%
filter(year(Month) >= 2017)
viii. We estimate an automatic ARIMA model using the training dataset.
Although we did the differencing manually above, we did not use it here! We are
using the automatic selection of parameters of ARIMA. Next week we try to do it
manually.
auto_arima <- emp_train %>%
model(ARIMA(Employed))
report(auto_arima)
ix. Forecasts are then generated for the test dataset.
fc <- auto_arima %>%
forecast(new_data = emp_test)
fc %>%
autoplot(emp_ts) +
labs(title = "ARIMA forecasts",
y = "Retail employment")
x. The forecasting accuracy can be measured by comparing the forecasts with the
actual values in the test dataset.
accuracy(fc, emp_test)
xi. Residual diagnostics are used to assess whether the ARIMA model has
adequately captured the structure of the time series.
This is usually done to verify that assumptions are met. Assumptions are that
residuals are independent, so, there should be no autocorrelation, and they are
normally distributed.
aug <- auto_arima %>%
augment()
autoplot(aug, .innov) +
labs(title = "Residuals from ARIMA model",
y = "Residuals")
aug %>%
ACF(.innov) %>%
autoplot() +
labs(title = "ACF of residuals")
aug %>%
ggplot(aes(x = .innov)) +
geom_histogram(bins = 30) +
labs(title = "Histogram of residuals",
x = "Residuals",
y = "Frequency")
aug %>%
features(.innov, ljung_box, lag = 24, dof = 0)
xii. Finally, we compare the ARIMA model with some exponential smoothing
methods.
This step is good as we understand that ARIMA, with all complexities that it has
compared to smoothing methods, is inferior for this data. So, we don’t always get
better forecast if we use more complex methods.
fit <- emp_train %>%
model(
Naive = NAIVE(Employed),
SNaive = SNAIVE(Employed),
ETS_ANN = ETS(Employed),
)
fc <- fit %>%
forecast(h = "3 years")
accuracy(fc, emp_ts)