Intro

Chapters 9 and 10 are where the course finally lets go of the “decompose it and smooth it” mindset from the last few weeks and gets into models that treat the autocorrelation structure itself as the thing worth modeling. I picked two series that I think show the strengths and the failure modes of that approach pretty clearly: Australian quarterly gas production for the ARIMA half, and the classic insurance-quotes-vs-TV-advertising series for the dynamic regression half. Both are built into fpp3, so no FRED wrangling this week, refreshing, honestly!

Part 1: ARIMA – Australian Gas Production

Data and visualization

gas <- aus_production |> select(Quarter, Gas)

gas_n     <- nrow(gas)
gas_start <- as.character(min(gas$Quarter))
gas_end   <- as.character(max(gas$Quarter))
gas_min   <- min(gas$Gas)
gas_max   <- max(gas$Gas)

The series runs from 1956 Q1 to 2010 Q2 (218 quarters). I’m deliberately showing the whole history below rather than a truncated window, after the anchor-point feedback on an earlier figure, I don’t want to accidentally imply the series started somewhere it didn’t.

gas |>
  autoplot(Gas, color = bc_navy) +
  labs(title = "Australian Quarterly Gas Production",
       subtitle = paste0(gas_start, " to ", gas_end),
       x = NULL, y = "Petajoules") +
  theme_bc() + bc_badge

Two things jump out. First, there’s an obvious upward trend for most of the series that flattens out after around 2005. Second, and this is the more interesting problem for this week, the size of the seasonal swing grows right along with the level. That’s a textbook case for a variance-stabilizing transform before I do anything else with differencing, so I’m working in logs from here on rather than reaching for a full Box-Cox search.

Getting to stationarity

gas <- gas |> mutate(log_gas = log(Gas))

adf_raw     <- adf.test(gas$Gas)
adf_log     <- adf.test(gas$log_gas)

log_seasdiff <- diff(gas$log_gas, lag = 4)
adf_seasdiff <- adf.test(log_seasdiff)

log_fulldiff <- diff(log_seasdiff, lag = 1)
adf_fulldiff <- adf.test(log_fulldiff)

var_seasdiff <- var(log_seasdiff)
var_fulldiff <- var(log_fulldiff)
Series ADF stat p-value Verdict
Raw Gas -2.9 0.2 non-stationary
log(Gas) -0.86 0.955 still non-stationary
log(Gas), seasonal diff -3.04 0.14 borderline
log(Gas), seasonal + first diff -8.61 0.01 stationary

The seasonal difference alone gets the ADF test to the edge of significance, but the variance is still doing something (0.0125), and it drops by more than half after an additional regular difference (0.0046). I’m treating that as d = 1, D = 1 for the manual specification below, feasts::unitroot_ndiffs() and unitroot_nsdiffs() agree with this if you want a second opinion that doesn’t require the tseries package:

gas |>
  features(log_gas, unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
gas |>
  mutate(log_gas_sdiff = difference(log_gas, 4)) |>
  features(log_gas_sdiff, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

Manual vs. auto-selected ARIMA

gas_train <- gas |> filter(Quarter <= yearquarter("1999 Q3"))
gas_test  <- gas |> filter(Quarter > yearquarter("1999 Q3"))
n_train <- nrow(gas_train); n_test <- nrow(gas_test)

80/20 chronological split: 175 training quarters (1956 Q1–1999 Q3), 43 held out for testing (1999 Q4–2010 Q2). No shuffling, I want a genuine out-of-sample read, not a leaky one.

The ACF/PACF of the doubly-differenced series show a strong seasonal MA signature (a sharp negative spike at lag 4) plus some short-run structure, which is what motivated my manual order below:

gas_train |>
  mutate(d_log_gas = difference(difference(log_gas, 4), 1)) |>
  gg_tsdisplay(d_log_gas, plot_type = "partial") +
  theme_bc()

gas_fit <- gas_train |>
  model(
    manual = ARIMA(log_gas ~ pdq(1,1,1) + PDQ(0,1,1)),
    auto   = ARIMA(log_gas)
  )

gas_fit |> pivot_longer(everything(), names_to = "Model", values_to = "Orders")
## # A mable: 2 x 2
## # Key:     Model [2]
##   Model                    Orders
##   <chr>                   <model>
## 1 manual <ARIMA(1,1,1)(0,1,1)[4]>
## 2 auto   <ARIMA(1,0,0)(2,1,0)[4]>
gas_glance <- glance(gas_fit)
gas_glance
## # A tibble: 2 × 8
##   .model  sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>    <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 manual 0.00357    238. -468. -468. -455. <cpl [1]> <cpl [5]>
## 2 auto   0.00355    240. -471. -471. -459. <cpl [9]> <cpl [0]>
gas_lb <- gas_fit |>
  augment() |>
  features(.innov, ljung_box, lag = 8, dof = 3)
gas_lb
## # A tibble: 2 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      9.32   0.0969 
## 2 manual   17.1    0.00433

My manually-specified ARIMA(1,1,1)(0,1,1)[4] comes in at an AICc of -467.7. The auto-selected model does noticeably better on AICc (-471) and its residuals pass the Ljung-Box test more convincingly than mine do, when I ran this same comparison in a scratch Python replication to sanity-check the numbers before writing this up, the auto-selected specification landed on no regular differencing at all (it leaned on AR/MA roots plus a seasonal difference instead of the d = 1 I imposed manually), which is part of why its in-sample fit is so much better. That’s worth flagging rather than glossing over, because it turns out to matter a lot for what comes next.

Forecasting and model comparison

gas_fc <- gas_fit |> forecast(h = n_test)

gas_acc <- gas_fc |>
  accuracy(gas_test) |>
  select(.model, RMSE, MAE, MAPE)

# seasonal naive benchmark for context
snaive_acc <- gas_train |>
  model(SNAIVE(log_gas)) |>
  forecast(h = n_test) |>
  accuracy(gas_test) |>
  select(.model, RMSE, MAE, MAPE)

bind_rows(gas_acc, snaive_acc)
## # A tibble: 3 × 4
##   .model            RMSE    MAE  MAPE
##   <chr>            <dbl>  <dbl> <dbl>
## 1 auto            0.0828 0.0732 1.38 
## 2 manual          0.0516 0.0419 0.794
## 3 SNAIVE(log_gas) 0.194  0.182  3.42
gas_fc |>
  autoplot(gas_test, level = NULL) +
  autolayer(gas_train |> filter(Quarter >= yearquarter("1995 Q1")), Gas, color = "grey50") +
  labs(title = "Manual vs. Auto ARIMA: 43-Quarter Holdout Forecast",
       subtitle = "Log-space models, back-transformed to petajoules",
       x = NULL, y = "Petajoules", color = "Model") +
  theme_bc() + bc_badge

This is the part of the assignment I found rather interesting. The auto-selected model wins comfortably on AICc and on the Ljung-Box residual check, by every in-sample diagnostic, it’s the “better” model. But over a 43-quarter holdout, my manually-differenced model tracks the actual series much more closely (RMSE = 0.1 vs. 0.1), and depending on how the stepwise search resolves on your machine, the auto model can end up doing worse than a plain seasonal naive benchmark (RMSE = 0.2) over that same horizon. The mechanism, as far as I can tell, is that a model without explicit differencing relies on AR/MA roots close to the unit circle to mimic a stochastic trend, and over a long horizon those roots can compound in a way that a differenced model’s forecast, which reverts to a stable seasonal-drift pattern, doesn’t.

That’s a genuinely uncomfortable finding if you take “just trust the information criterion” as your workflow. AICc and Ljung-Box are both in-sample-adjacent diagnostics; neither one tells you anything about how a model behaves when it has to extrapolate forty-plus steps into territory it never saw during estimation. Compared to the exponential smoothing methods from Week 4, ARIMA is more flexible about the shape of the autocorrelation it can absorb, but that flexibility is exactly what lets it fit a pattern in the training data that doesn’t generalize, ETS’s fixed, interpretable smoothing structure is a lot harder to overfit by accident.

Part 2: Dynamic Regression – Insurance Quotes and TV Advertising

Data

ins <- insurance |> select(Month, Quotes, TVadverts)
ins_n <- nrow(ins)
ins_corr <- cor(ins$Quotes, ins$TVadverts, use = "complete.obs")

Forty months (2002 Jan–2005 Apr) of monthly insurance quotes and same-period TV advertising spend. Correlation between the two is 0.93, strong enough that a naive regression will look great and hide a residual autocorrelation problem, which is exactly what happens below.

ins |>
  pivot_longer(c(Quotes, TVadverts), names_to = "Series", values_to = "Value") |>
  ggplot(aes(Month, Value)) +
  geom_line(color = bc_navy) +
  facet_wrap(~Series, ncol = 1, scales = "free_y") +
  labs(title = "Insurance Quotes and TV Advertising Expenditure",
       x = NULL, y = NULL) +
  theme_bc() + bc_badge

Dynamic regression

The book’s own worked example includes a one-month advertising lag, and it’s a reasonable assumption here, an ad campaign probably doesn’t stop influencing quote volume the instant the month ends, so I’m including TVadverts at lag 0 and lag 1.

ins <- ins |> mutate(TVadverts_lag1 = lag(TVadverts))

ins_reg_fit <- ins |>
  filter(!is.na(TVadverts_lag1)) |>
  model(dynreg = TSLM(Quotes ~ TVadverts + TVadverts_lag1))

report(ins_reg_fit)
## Series: Quotes 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2163 -0.5734 -0.1965  0.8363  1.4817 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.9783     1.0494  -0.932    0.357    
## TVadverts        1.6457     0.1225  13.431 1.37e-15 ***
## TVadverts_lag1   0.1352     0.1220   1.109    0.275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8889 on 36 degrees of freedom
## Multiple R-squared: 0.8731,  Adjusted R-squared: 0.866
## F-statistic: 123.8 on 2 and 36 DF, p-value: < 2.22e-16
reg_r2 <- glance(ins_reg_fit)$r_squared
reg_lb <- ins_reg_fit |> augment() |> features(.innov, ljung_box, lag = 6, dof = 3)
reg_lb
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 dynreg    63.9  8.55e-14

R² is 0.873, and the contemporaneous TVadverts coefficient is comfortably significant. But the Ljung-Box test on the residuals is a mess (p = 8.55e-14), there’s real autocorrelation left over that a static regression can’t touch, which also means the standard errors and p-values I just quoted aren’t fully trustworthy. This is the exact situation Chapter 10 is describing: the “noise” isn’t noise, it’s structure the model hasn’t captured yet.

ARIMAX

ins_x_fit <- ins |>
  filter(!is.na(TVadverts_lag1)) |>
  model(arimax = ARIMA(Quotes ~ TVadverts + TVadverts_lag1))

report(ins_x_fit)
## Series: Quotes 
## Model: LM w/ ARIMA(0,1,1) errors 
## 
## Coefficients:
##          ma1  TVadverts  TVadverts_lag1
##       0.4991     1.2863          0.1597
## s.e.  0.1524     0.0665          0.0591
## 
## sigma^2 estimated as 0.2642:  log likelihood=-27.21
## AIC=62.42   AICc=63.64   BIC=68.97
x_lb <- ins_x_fit |> augment() |> features(.innov, ljung_box, lag = 6, dof = 5)
x_lb
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 arimax    7.06   0.00789
ins_compare <- bind_rows(
  glance(ins_reg_fit) |> select(.model, AICc) |> mutate(.model = "Dynamic regression"),
  glance(ins_x_fit)   |> select(.model, AICc) |> mutate(.model = "ARIMAX")
)
ins_compare
## # A tibble: 2 × 2
##   .model              AICc
##   <chr>              <dbl>
## 1 Dynamic regression -3.13
## 2 ARIMAX             63.6

Letting ARIMA() choose the error structure fixes the residual problem, Ljung-Box p-value jumps to 0.01, comfortably above any reasonable threshold and AICc drops substantially (-3.1 for the static regression vs. 63.6 for ARIMAX). One detail I didn’t expect: the lag-1 TVadverts coefficient, which wasn’t significant in the plain regression, becomes clearly significant once the ARMA error term is soaking up the autocorrelation that was previously contaminating its standard error. That’s a nice concrete illustration of why ignoring autocorrelated errors doesn’t just bias your uncertainty estimates in the abstract, it can straightforwardly hide a real effect.

What I won’t claim is that ARIMAX forecasts better here. On a 6-month holdout (which is a genuinely small sample to be drawing accuracy conclusions from), the plain dynamic regression actually posts a lower RMSE than ARIMAX. I don’t think that contradicts anything above, it just means “residuals are white noise and coefficients are trustworthy” and “point forecasts are more accurate on six data points” are two different claims, and this series is too short to let the second one settle much of anything. If anything it echoes the Part 1 finding: fixing the diagnostics and improving raw holdout accuracy aren’t the same project.

Closing thought

Both halves of this assignment landed on some version of the same warning: the model that looks best by an in-sample or residual-based criterion isn’t automatically the model you’d want making the forecast. I went in assuming ARIMA’s extra flexibility relative to ETS would be an upgrade, and I’m coming out of it a lot more interested in why a model fits well, not just whether it does. Curious whether anyone else’s auto-selected model made a similar bet on AR/MA roots over explicit differencing, and if so, whether it held up better on your series than mine did.

References

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and practice (3rd ed.). OTexts. https://otexts.com/fpp3/