Discussion 2: Trends, Seasonality, and Simple Benchmarks

Victorian Liquor Retailing

Data and why I picked it

I’m working with the aus_retail dataset from the fpp3 package: monthly retail turnover for liquor retailing in Victoria, April 1982 through December 2018 (441 observations, $AUD millions, originally from the ABS Retail Trade Survey).

Last week my MAUISA series turned out to be pre-seasonally-adjusted by the BEA, so its seasonal component was basically an artifact (seasonal strength around 0.01). This week I deliberately went the other direction and picked a raw, unadjusted retail series where I expected the seasonality to be real and loud. Liquor sales around the holidays seemed like a safe bet.

library(fpp3)
library(forecast)
library(ggplot2)
library(gridExtra)
library(grid)

vic <- aus_retail |>
  filter(State == "Victoria", Industry == "Liquor retailing") |>
  arrange(Month)
y <- ts(vic$Turnover, start = c(1982, 4), frequency = 12)

cat("Series length:", length(y), "obs | Range:", min(y), "-", max(y), "$M AUD")

## Series length: 441 obs | Range: 17.3 - 342.9 $M AUD

What the plots show

df <- data.frame(t = as.numeric(time(y)), v = as.numeric(y))
p1 <- ggplot(df, aes(t, v)) +
  geom_line(color = INK, linewidth = 0.45) +
  annotate("text", x = 2016.2, y = 336, hjust = 1, size = 3, color = "grey30",
           label = "Dec 2018: $336.8M") +
  annotate("point", x = 2018 + 11/12, y = 336.8, size = 1.4, color = INK) +
  labs(title = "Victorian Liquor Retailing Turnover, Apr 1982 - Dec 2018",
       subtitle = "aus_retail (fpp3) | ABS Retail Trade Survey, $AUD millions",
       x = NULL, y = "Turnover ($M)") +
  theme_bc()
show_badged(p1)

Figure 1. Monthly turnover, Victorian liquor retailing. Point marks the series maximum.

The time plot shows turnover climbing from $17.3M to $336.8M over 37 years, roughly a 20x increase, with a sawtooth pattern layered on top. The seasonal plot and monthly distributions make the source of the sawtooth obvious: December turnover averages 1.62 times the overall mean, November runs elevated, and January dips as households apparently recover from the holidays.

df$yr <- floor(df$t); df$mo <- cycle(y)
pa <- ggplot(df, aes(mo, v, group = yr, color = yr)) +
  geom_line(linewidth = 0.35) +
  scale_color_gradient(low = "grey85", high = INK, name = "Year",
                       breaks = c(1985, 2000, 2015)) +
  scale_x_continuous(breaks = 1:12, labels = month.abb) +
  labs(title = "Seasonal plot", subtitle = "one line per year; darker = more recent",
       x = NULL, y = "Turnover ($M)") +
  theme_bc() + theme(legend.position = "right", legend.key.height = unit(0.5, "cm"))
mo_means <- aggregate(v ~ mo, df, mean)
pb <- ggplot(df, aes(factor(mo, labels = month.abb), v)) +
  geom_point(color = INK, alpha = 0.25, size = 0.7,
             position = position_jitter(width = 0.12, height = 0)) +
  geom_point(data = data.frame(mo = factor(mo_means$mo, labels = month.abb), v = mo_means$v),
             aes(x = mo, y = v), color = INK, size = 2.6, shape = 18) +
  labs(title = "Monthly distribution", subtitle = "dots = individual months; diamond = mean",
       x = NULL, y = "Turnover ($M)") +
  theme_bc()
grid.arrange(pa, pb, ncol = 2, bottom = badge_grob())

Figure 2. Left: seasonal plot, darker ink = more recent year. Right: turnover distribution by month; diamonds mark monthly means.

cat("Mean December turnover / overall mean:",
    round(tapply(y, cycle(y), mean)[12] / mean(y), 2), "\n")

## Mean December turnover / overall mean: 1.62

early <- window(y, end = c(1989, 12)); late <- window(y, start = c(2012, 1))
cat("Seasonal range 1982-89:", round(max(early) - min(early), 1), "$M |",
    "2012-18:", round(max(late) - min(late), 1), "$M\n")

## Seasonal range 1982-89: 34 $M | 2012-18: 203.7 $M

What I found more interesting is that the seasonal swings grow with the level of the series. The within-year range was about $34M in the 1980s and about $204M in 2012–2018. That is the signature of multiplicative seasonality (Hyndman & Athanasopoulos, 2021, sec. 3.2), and it drove two choices downstream: a multiplicative classical decomposition, and a log transformation before STL.

Classical vs. STL decomposition

dec_classical <- decompose(y, type = "multiplicative")
show_badged(autoplot(dec_classical) +
  labs(title = "Classical Multiplicative Decomposition", x = NULL) + theme_bc())

Figure 3. Classical multiplicative decomposition.

dec_stl <- stl(log(y), s.window = 13, robust = TRUE)
show_badged(autoplot(dec_stl) +
  labs(title = "STL Decomposition of log(Turnover), s.window = 13, robust",
       x = NULL) + theme_bc())

Figure 4. STL decomposition of log(turnover), s.window = 13, robust.

Tt <- dec_stl$time.series[, "trend"]
St <- dec_stl$time.series[, "seasonal"]
Rt <- dec_stl$time.series[, "remainder"]
F_trend    <- max(0, 1 - var(Rt) / var(Tt + Rt))
F_seasonal <- max(0, 1 - var(Rt) / var(St + Rt))
cat("STL trend strength Ft =", round(F_trend, 3),
    "| seasonal strength Fs =", round(F_seasonal, 3), "\n")

## STL trend strength Ft = 0.996 | seasonal strength Fs = 0.893

St_mat <- matrix(St[1:432], ncol = 12, byrow = TRUE)
cat("STL December seasonal effect (log scale), first 5 yrs vs last 5 yrs:",
    round(mean(St_mat[1:5, 9]), 3), "vs", round(mean(St_mat[32:36, 9]), 3), "\n")

## STL December seasonal effect (log scale), first 5 yrs vs last 5 yrs: 0.575 vs 0.479

cat("Classical Dec seasonal index (fixed):", round(dec_classical$figure[9], 3), "\n")

## Classical Dec seasonal index (fixed): 1.581

On the surface the two methods tell a similar story: dominant trend, regular seasonality, modest remainder. The STL-based strength measures came out at Ft = 0.996 and Fs = 0.893, almost exactly the mirror image of last week’s MAUISA (0.89 and 0.01).

The real difference shows up in how each method treats the seasonal component. Classical decomposition assumes one fixed seasonal index for the whole sample: December is locked at 1.581 whether it is 1985 or 2018. STL lets the pattern evolve, and here that matters. The December effect on the log scale shrinks from about 0.575 in the first five years to 0.479 in the last five. In plain terms, the holiday spike is still huge but proportionally fading; Victorians appear to be spreading their liquor spending more evenly across the year. Classical decomposition is structurally incapable of detecting that, which to me is the strongest argument for STL, beyond the usual textbook points about robustness to outliers and not losing observations at the ends (the classical method drops six months at each end to its 2x12 moving average).

rem_classical <- na.omit(dec_classical$random)
cat("Ljung-Box on classical remainder: p =",
    format.pval(Box.test(rem_classical, lag = 24, type = "Ljung-Box")$p.value, digits = 3), "\n")

## Ljung-Box on classical remainder: p = 2.31e-11

cat("Ljung-Box on STL remainder:       p =",
    format.pval(Box.test(Rt, lag = 24, type = "Ljung-Box")$p.value, digits = 3), "\n")

## Ljung-Box on STL remainder:       p = 3.45e-10

One humbling similarity: neither method produces a clean remainder. A Ljung-Box test on the classical remainder gives p ≈ 2e-11, and the STL remainder is barely better at p ≈ 3e-10. Same lesson as last week: decomposition describes the series, it does not fully explain it.

The forecaster’s toolbox: 70/30 split

I trained on the first 308 observations (through November 2007) and held out the remaining 133 (December 2007 – December 2018), then ran the four benchmarks from section 5.2: mean, naïve, seasonal naïve, and drift.

n     <- length(y)
n_trn <- floor(0.70 * n)
train <- window(y, end   = time(y)[n_trn])
test  <- window(y, start = time(y)[n_trn + 1])
h     <- length(test)

fc_mean   <- meanf(train, h = h)
fc_naive  <- naive(train, h = h)
fc_snaive <- snaive(train, h = h)
fc_drift  <- rwf(train, h = h, drift = TRUE)

acc_row <- function(fc, label) {
  e  <- as.numeric(test) - as.numeric(fc$mean)
  pe <- 100 * e / as.numeric(test)
  data.frame(Method = label, ME = round(mean(e), 2), MPE = round(mean(pe), 2),
             MAPE = round(mean(abs(pe)), 2), MSE = round(mean(e^2), 1),
             MAE = round(mean(abs(e)), 2))
}
acc <- rbind(acc_row(fc_mean, "Mean"), acc_row(fc_naive, "Naive"),
             acc_row(fc_snaive, "Seasonal naive"), acc_row(fc_drift, "Drift"))
knitr::kable(acc, caption = "Table 1. Hold-out accuracy, h = 133 months.")

Table 1. Hold-out accuracy, h = 133 months.

Method	ME	MPE	MAPE	MSE	MAE
Mean	122.60	70.11	70.11	17229.0	122.60
Naive	60.92	31.24	33.21	5910.2	62.75
Seasonal naive	69.08	38.04	38.24	6086.2	69.30
Drift	40.88	20.09	22.31	3263.0	42.96

p5 <- autoplot(window(y, start = c(2000, 1)), color = "grey55", linewidth = 0.4) +
  autolayer(fc_mean$mean, series = "Mean") +
  autolayer(fc_naive$mean, series = "Naive") +
  autolayer(fc_snaive$mean, series = "Seasonal naive") +
  autolayer(fc_drift$mean, series = "Drift") +
  scale_color_manual(values = c("Mean" = "#B23A48", "Naive" = "#E0A100",
                                "Seasonal naive" = "#5B8C5A", "Drift" = "#1F3864")) +
  labs(title = "Benchmark Forecasts vs. Actuals (test: Dec 2007 - Dec 2018)",
       subtitle = "grey = observed turnover; trained on Apr 1982 - Nov 2007 (70%)",
       x = NULL, y = "Turnover ($M)", colour = "Method") +
  theme_bc()
show_badged(p5)

Figure 5. The four benchmark forecasts against the observed test data (grey).

Two things surprised me. First, drift wins comfortably (MAPE 22.3%), which makes sense in hindsight: with trend strength at 0.996, the only method that extrapolates any growth at all is going to dominate over an 11-year horizon. Second, the seasonal naïve lost to the plain naïve (38.2% vs. 33.2% MAPE) despite seasonal strength of 0.89. I stared at that for a while before it clicked. Every method except drift forecasts a flat line or a flat repeating pattern, so far out into the test set the error is almost entirely trend bias. Snaïve repeats 2007’s seasonal shape at 2007’s level, and by 2018 actual turnover has nearly doubled; seasonal fidelity does not matter when the level is wrong by 40%. All four ME values are strongly positive, confirming systematic underprediction across the board.

r <- residuals(fc_snaive)
ra <- autoplot(r, color = INK, linewidth = 0.35) +
  labs(title = "Seasonal naive training residuals", x = NULL, y = "Residual ($M)") + theme_bc()
rb <- ggAcf(r, lag.max = 36) + labs(title = "ACF of residuals", x = "Lag", y = "ACF") +
  theme_bc() + theme(plot.title = element_text(size = 11))
rc <- ggplot(data.frame(r = as.numeric(na.omit(r))), aes(r)) +
  geom_histogram(bins = 30, fill = INK, color = "white", linewidth = 0.2) +
  labs(title = "Histogram of residuals", x = "Residual ($M)", y = "Count") + theme_bc()
grid.arrange(ra, arrangeGrob(rb, rc, ncol = 2), nrow = 2, bottom = badge_grob())

Figure 6. Residual diagnostics for the seasonal naive method (training residuals).

Box.test(residuals(fc_snaive), lag = 24, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals(fc_snaive)
## X-squared = 601.3, df = 24, p-value < 2.2e-16

Honestly, I am a little skeptical of my results/measures. A 133-month holdout is not a forecasting scenario anyone would face; MAPE at that horizon is mostly a referendum on whether a method has a trend term, and percentage errors shrink mechanically as the denominator grows in a trending series. It also does not help that my training window ends one month before the global financial crisis. Section 5.10’s time series cross-validation seems like a more honest evaluation design, and I would expect the seasonal naïve to look far better at h = 12 than at h = 133.

The obvious next step is a method that carries both components, Holt-Winters or a seasonal ARIMA, since the decomposition says this series is essentially trend plus multiplicative seasonality with a little noise left over.

Question for the class: Did anyone else’s benchmark ranking flip depending on where they cut the split? Has anyone tried tsCV() instead? I would like to hear how it compared.

References

Australian Bureau of Statistics. (2019). Retail trade, Australia (Cat. No. 8501.0) [Data set]. Accessed via the fpp3 R package (aus_retail). https://www.abs.gov.au/

Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal-trend decomposition procedure based on loess. Journal of Official Statistics, 6(1), 3–73.

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