---
title: "Module 8 Discussion: ARIMA, Seasonal ARIMA, White Noise, Random Walk, and Drift"
author: "Adrian Aziza"
format:
html:
toc: true
code-fold: false
code-tools: true
execute:
echo: true
warning: false
message: false
---
```{r setup}
# Package setup ------------------------------------------------------------
pkgs <- c("forecast", "ggplot2", "dplyr", "tidyr", "tibble", "patchwork", "knitr")
to_install <- setdiff(pkgs, rownames(installed.packages()))
if (length(to_install) > 0) {
install.packages(to_install, repos = "https://cloud.r-project.org")
}
invisible(lapply(pkgs, library, character.only = TRUE))
theme_set(theme_minimal(base_size = 12))
fmt_arima <- function(fit) {
ord <- forecast::arimaorder(fit)
if (all(c("P", "D", "Q", "Frequency") %in% names(ord)) && ord[["Frequency"]] > 1) {
paste0("ARIMA(", ord[["p"]], ",", ord[["d"]], ",", ord[["q"]], ")(",
ord[["P"]], ",", ord[["D"]], ",", ord[["Q"]], ")[", ord[["Frequency"]], "]")
} else {
paste0("ARIMA(", ord[["p"]], ",", ord[["d"]], ",", ord[["q"]], ")")
}
}
```
# I. Manual ARIMA selection versus automatic ARIMA
For the non-seasonal example, I use the built-in annual `Nile` river-flow series. I hold out the last 10 years as the actual comparison period and fit models only on the historical training data.
```{r part1-data}
nile <- ts(as.numeric(Nile), start = 1871, frequency = 1)
nile_train <- window(nile, end = 1960)
nile_test <- window(nile, start = 1961)
nile_info <- tibble(
split = c("Training", "Holdout / actual"),
years = c("1871-1960", "1961-1970"),
n = c(length(nile_train), length(nile_test))
)
knitr::kable(nile_info)
```
The original series is not centered around a stable mean, so I difference it once. After first differencing, the ACF has its main meaningful spike at lag 1 and then weakens, while the PACF does not show a strong multi-lag cutoff. That pattern supports a low-order MA model on the differenced series, so my manual choice is **ARIMA(0,1,1)**: no AR term, one regular difference, and one MA term.
```{r part1-acf-pacf, fig.width=10, fig.height=4}
p_acf <- ggAcf(diff(nile_train), lag.max = 24) +
labs(title = "ACF of first-differenced Nile training series")
p_pacf <- ggPacf(diff(nile_train), lag.max = 24) +
labs(title = "PACF of first-differenced Nile training series")
p_acf + p_pacf
```
```{r part1-models}
manual_nile <- Arima(nile_train, order = c(0, 1, 1))
auto_nile <- auto.arima(
nile_train,
seasonal = FALSE,
stepwise = FALSE,
approximation = FALSE
)
model_compare <- tibble(
model = c("Manual", "Automatic"),
arima_form = c(fmt_arima(manual_nile), fmt_arima(auto_nile)),
AICc = c(manual_nile$aicc, auto_nile$aicc),
BIC = c(BIC(manual_nile), BIC(auto_nile))
)
knitr::kable(model_compare, digits = 2)
summary(manual_nile)
summary(auto_nile)
```
The automatic model is close to the manual model because both use first differencing and a low-order short-memory error structure. If `auto.arima()` adds one AR term, the difference is not conceptual: it means the algorithm found a small information-criterion improvement by allowing mild persistence in the differenced errors. I would report the automatic model for forecasting because it is chosen by AICc, but the ACF/PACF justification still correctly identified the needed differencing and the low-order ARMA structure.
```{r part1-forecast, fig.width=9, fig.height=5}
nile_fc <- forecast(auto_nile, h = length(nile_test), level = 95)
autoplot(nile_fc) +
autolayer(nile_test, series = "Actual holdout", linewidth = 0.8) +
labs(
title = "Nile annual flow: historical values, actual holdout, forecasts, and 95% prediction intervals",
x = "Year",
y = "Annual flow",
color = "Series"
)
```
```{r part1-forecast-table}
nile_forecast_table <- tibble(
year = 1961:1970,
actual = as.numeric(nile_test),
prediction = as.numeric(nile_fc$mean),
lo_95 = as.numeric(nile_fc$lower[, "95%"]),
hi_95 = as.numeric(nile_fc$upper[, "95%"])
)
knitr::kable(nile_forecast_table, digits = 1)
```
# II. Seasonal ARIMA
For the seasonal example, I use the built-in monthly `AirPassengers` series. It has a clear upward trend and multiplicative annual seasonality, so I use a log transformation through `lambda = 0`. The seasonal ARIMA model is:
\[
(1-B)(1-B^{12})\log(y_t)=(1+\theta_1B)(1+\Theta_1B^{12})\varepsilon_t.
\]
This is an **ARIMA(0,1,1)(0,1,1)[12]** model: one regular difference removes trend, one seasonal difference removes yearly seasonality, the non-seasonal MA(1) term captures short-run one-month error correction, and the seasonal MA(1) term captures same-month-last-year error correction.
```{r part2-seasonal-model}
ap <- AirPassengers
ap_sarima <- Arima(
ap,
order = c(0, 1, 1),
seasonal = c(0, 1, 1),
lambda = 0,
biasadj = TRUE
)
summary(ap_sarima)
coef_table <- tibble(
coefficient = names(coef(ap_sarima)),
estimate = as.numeric(coef(ap_sarima)),
interpretation = c(
"Non-seasonal MA(1): correction for the previous month's forecast error on the log-differenced scale.",
"Seasonal MA(1): correction for the same month in the previous year's seasonal forecast error."
)
)
knitr::kable(coef_table, digits = 3)
```
The negative MA coefficients mean that positive forecast errors tend to be partially corrected in the next relevant period. The regular MA(1) coefficient works month-to-month, while the seasonal MA(1) coefficient works across the 12-month seasonal lag. Because the model is estimated after one regular difference and one seasonal difference, the coefficients describe remaining dependence after trend and annual seasonality have been removed.
```{r part2-plot, fig.width=9, fig.height=5}
ap_fc <- forecast(ap_sarima, h = 24, level = 95)
autoplot(ap_fc) +
labs(
title = "Seasonal ARIMA forecast for AirPassengers",
x = "Year",
y = "Monthly passengers",
color = "Series"
)
```
```{r part2-residuals, fig.width=9, fig.height=5}
checkresiduals(ap_sarima)
```
# III. White noise, random walk, and random walk with drift
Let \(\varepsilon_t\) be independent random shocks with mean zero and constant variance \(\sigma^2\):
\[
\varepsilon_t \sim WN(0,\sigma^2).
\]
**White noise:**
\[
y_t=\varepsilon_t.
\]
Each value is just the current shock. There is no memory because the previous shock does not carry into the next value.
**Random walk:**
\[
y_t=y_{t-1}+\varepsilon_t.
\]
The current value equals the previous value plus a new shock. Differencing gives:
\[
\Delta y_t=y_t-y_{t-1}=\varepsilon_t,
\]
so the first difference of a random walk is white noise.
**Random walk with drift:**
\[
y_t=c+y_{t-1}+\varepsilon_t.
\]
Here \(c\) is a constant drift. Differencing gives:
\[
\Delta y_t=y_t-y_{t-1}=c+\varepsilon_t,
\]
so the first difference is white noise plus a constant.
```{r part3-simulate, fig.width=9, fig.height=7}
set.seed(7406)
n <- 160
t <- 1:n
sigma <- 1
drift <- 0.15
eps <- rnorm(n, mean = 0, sd = sigma)
sim_df <- tibble(
t = t,
`White noise` = eps,
`Random walk` = cumsum(eps),
`Random walk with drift` = cumsum(drift + eps)
)
sim_long <- sim_df |>
pivot_longer(-t, names_to = "process", values_to = "value")
ggplot(sim_long, aes(x = t, y = value)) +
geom_line(linewidth = 0.6) +
facet_wrap(~ process, scales = "free_y", ncol = 1) +
labs(
title = "Simulated white noise, random walk, and random walk with drift",
x = "Time",
y = "Value"
)
```
```{r part3-comparison-table}
comparison <- tibble(
process = c("White noise", "Random walk", "Random walk with drift"),
trend = c("None; fluctuates around a constant mean", "No deterministic trend, but it wanders", "Upward if c > 0; downward if c < 0"),
shock = c("Only affects the current period", "Changes the level from that point forward", "Changes the level permanently, plus drift moves the series"),
memory = c("None", "Permanent", "Permanent"),
stationarity = c("Yes, if shocks are iid with constant variance", "No", "No"),
forecast = c("Constant at the mean, usually 0", "Flat at the latest observed value, with widening intervals", "Trending line from the latest value, with widening intervals")
)
knitr::kable(comparison)
```
Visually, white noise looks like random jitter around a stable center. The random walk looks smoother but less stable because each shock accumulates into the level. The random walk with drift also accumulates shocks, but it has a systematic direction because the constant drift is added every period.
