Discussion 8
Load Data:
I. ACF/PACF, ARIMA:
This plot suggests the original series is nonstationary. The mean level changes a lot over time, especially with the strong run-up in the 2000s, the sharp collapse around 2008–2009, and then the rebound afterward. Because the series has clear trend/level shifts and changing volatility, it does not appear appropriate to model it directly with a stationary ARMA process, which supports taking a first difference.
The first-differenced series looks much closer to stationary in mean. After differencing, the long upward trend is largely removed, and the series fluctuates around zero rather than around a changing level, which supports setting \(d=1\). That said, there are still some large shocks and periods of higher volatility, especially in the late 2000s, so while differencing improves stationarity, the series is not perfectly stable in variance.
op |>
mutate(diff_oil = difference(value)) |>
gg_tsdisplay(diff_oil, plot_type = "partial")
In the differenced series, both the ACF and PACF show a strong spike at lag 1, while the remaining lags mostly taper off rather than cutting off sharply. This pattern suggests that both autoregressive and moving average behavior may be present at low order, so ARIMA(1,1,1) is a reasonable model choice.
fit_auto <- op |>
model(auto_arima = ARIMA(value))
report(fit_auto)Series: value
Model: ARIMA(1,1,0)
Coefficients:
ar1
0.3894
s.e. 0.0463
sigma^2 estimated as 12.21: log likelihood=-1072.87
AIC=2149.75 AICc=2149.78 BIC=2157.74The automatically selected ARIMA model was ARIMA(1,1,0), while my manual choice was ARIMA(1,1,1). These models are very similar, since both include one order of differencing and one autoregressive term. This suggests that my visual interpretation of the ACF and PACF was broadly correct, especially in identifying the need for \(d=1\) and a low-order AR structure. However, the automatic procedure preferred the simpler ARIMA(1,1,0) model because it produced the lowest AIC, AICc, and BIC. In other words, the MA(1) term in my manual model did not improve fit enough to justify the added complexity. This is consistent with how the ARIMA() function in fpp3 works: it uses repeated KPSS tests to select differencing and then minimizes AICc over candidate models using a stepwise search.
II. Seasonal ARIMA:
The retail turnover series has a clear upward trend and visible recurring seasonal peaks, so it does not appear stationary in its original form. The seasonal swings also seem to grow larger over time, suggesting increasing variance as the level rises. This indicates that both trend and seasonality should be accounted for when modeling the series.
fit_retail <- retail |>
model(auto_sarima = ARIMA(Turnover))
report(fit_retail)Series: Turnover
Model: ARIMA(0,1,2)(0,1,1)[12]
Coefficients:
ma1 ma2 sma1
-0.2759 -0.1023 -0.6761
s.e. 0.0491 0.0515 0.0350
sigma^2 estimated as 250: log likelihood=-1791.13
AIC=3590.27 AICc=3590.36 BIC=3606.5The fitted seasonal ARIMA model is ARIMA(0,1,2)(0,1,1)[12]. This means the series was differenced once to remove the overall trend and seasonally differenced once at lag 12 to remove the yearly seasonal pattern, which is appropriate for monthly retail data. The model includes two non-seasonal moving average terms and one seasonal moving average term, so the current month’s turnover is explained by recent forecast errors as well as forecast errors from the same season one year earlier. The estimated coefficients are all negative: ma1 = -0.2759, ma2 = -0.1023, and sma1 = -0.6761. The negative MA terms suggest that unexpected shocks tend to be partially corrected in the following month or two, rather than continuing in the same direction. The seasonal MA coefficient is the largest in magnitude, which indicates that seasonal shocks from 12 months earlier play an important role in explaining current turnover. Since each coefficient is several standard errors away from zero, they all appear to be statistically meaningful contributors to the model.
III. WN/RW/RW w/D:
White noise is the simplest process and represents pure randomness. It can be written as \(y_t = \varepsilon_t\), where \(\varepsilon_t\) is a random shock with mean \(E(\varepsilon_t)=0\), constant variance \(\mathrm{Var}(\varepsilon_t)=\sigma^2\), and \(\mathrm{Cov}(\varepsilon_t,\varepsilon_{t-k})=0\) for all \(k \neq 0\). In words, white noise has no trend, no persistence, and no predictable structure, so knowing the past tells us nothing about the future.
A random walk is formed by accumulating white noise over time. Its equation is \(y_t = y_{t-1} + \varepsilon_t\), where \(y_{t-1}\) is last period’s value and \(\varepsilon_t\) is a white noise shock. Each new value equals the previous value plus a random change, so shocks permanently affect the level of the series. Because the shocks build up over time, the variance of \(y_t\) grows with \(t\), which means a random walk is nonstationary. However, if we difference it, we get \(\Delta y_t = y_t - y_{t-1} = \varepsilon_t\), so the first difference of a random walk is just white noise.
A random walk with drift adds a constant average change each period. Its equation is \(y_t = \delta + y_{t-1} + \varepsilon_t\), where \(\delta\) is the drift term and represents a systematic upward or downward trend over time. If \(\delta > 0\), the series tends to rise on average; if \(\delta < 0\), it tends to fall. Differencing gives \(\Delta y_t = y_t - y_{t-1} = \delta + \varepsilon_t\). The differenced series is white noise plus a constant.
Overall, White noise is the stationary building block, the random walk is its cumulative sum, and the random walk with drift is the cumulative sum of white noise plus a fixed average change.
set.seed(111)
# Simulate 200 observations
n <- 200
wn <- tibble(
Time = 1:n,
y = rnorm(n, mean = 0, sd = 1)
) |>
as_tsibble(index = Time)
rw <- wn |>
mutate(y = cumsum(y))
rw_drift <- wn |>
mutate(y = cumsum(0.5 + y))
# Plot white noise
p1 <- wn |>
autoplot(y) +
labs(
title = "White Noise",
x = "Time",
y = "Value"
)
# Plot random walk
p2 <- rw |>
autoplot(y) +
labs(
title = "Random Walk",
x = "Time",
y = "Value"
)
# Plot random walk with drift
p3 <- rw_drift |>
autoplot(y) +
labs(
title = "Random Walk with Drift",
x = "Time",
y = "Value"
)
p1 / p2 / p3
The white noise series has no trend, since it fluctuates randomly around a constant level near zero rather than moving persistently upward or downward. Each shock is just a one-time random jump at time \(t\), and those shocks have no memory, meaning they do not affect future values at all. Because its mean, variance, and autocovariances stay constant over time, white noise is stationary. As a result, its forecast is flat at the mean, because there is no structure in the past to help predict anything else.
The random walk also has no deterministic trend, but visually it can drift up or down for long stretches because the shocks accumulate over time. A shock to a random walk does not disappear after one period. It rather becomes part of the level of the series, so the process has permanent memory. This is why the path appears to wander and why the variance grows over time, making the random walk nonstationary. Its forecast is flat in the sense that the best prediction is the current level carried forward, but that flat forecast is centered on the most recent value rather than on a fixed long-run mean.
The random walk with drift shows a clear upward trend in this simulation because each period includes both a random shock and a positive constant drift term. The shock still has a permanent effect, just as in the ordinary random walk, so the process again has permanent memory. It is also nonstationary, since both the level and the variance change over time as shocks accumulate. Its forecast is trending, because future values are projected by extending the current level forward while adding the average drift each period.