Week 6 Discussion: Hierarchical / Grouped Forecasting and GARCH

Introduction

This week’s reading covers three things that don’t usually show up in the same sentence: reconciling forecasts across a hierarchy, reconciling them across crossed groupings, and modeling the variance of a series instead of just its mean. I used the PBS dataset (Medicare Australia subsidised prescriptions) for both Parts 1 and 2, since it’s both hierarchical (drug classification) and grouped (concession type crossed with drug classification), and it’s the same dataset the textbook leans on for this chapter. For Part 3 I switched to daily stock returns (gafa_stock, Amazon), since PBS has no business having time-varying variance and GARCH on a series with constant variance is a pretty boring demo.

I held out the last 24 months (Jul 2006 – Jun 2008) as a test set for Parts 1 and 2 and fit everything on the remaining 180 months.

---

Part 1: Hierarchical Forecasting (PBS by ATC1/ATC2)

Structure

PBS scripts are classified by the Anatomical Therapeutic Chemical index: 15 top-level categories (ATC1, e.g. “cardiovascular system”), each split into more specific sub-categories (ATC2), for 84 bottom-level series total. That’s a clean parent/child hierarchy, so I used aggregate_key(ATC1/ATC2, ...) rather than the crossed syntax I needed for Part 2.

pbs_hier <- PBS |>
  aggregate_key(ATC1/ATC2, Scripts = sum(Scripts))

atc1_avg <- pbs_hier |>
  filter(!is_aggregated(ATC1), is_aggregated(ATC2)) |>
  as_tibble() |>
  group_by(ATC1) |>
  summarise(avg_scripts = mean(Scripts) / 1e3)

atc1_avg |>
  mutate(ATC1 = fct_reorder(as.character(ATC1), avg_scripts)) |>
  ggplot(aes(x = avg_scripts, y = ATC1)) +
  geom_segment(aes(x = 0, xend = avg_scripts, yend = ATC1), color = "grey80") +
  geom_point(size = 2.5, color = navy) +
  labs(
    title = "Average monthly PBS scripts by ATC1 category",
    subtitle = "1991 Jul -- 2006 Jun (training period), 15 categories feeding the hierarchy",
    x = "Average scripts per month (thousands)", y = NULL,
    caption = bc_badge
  ) +
  theme_bc()

[ my other dataset is AlphaFold ]

The category sizes are wildly uneven, some ATC1 groups are an order of magnitude bigger than others, which matters later: a reconciliation method that improves the biggest series and hurts a small one still looks good on an aggregate error metric.

Base models and three classic reconciliation approaches

train_hier <- pbs_hier |> filter(Month <= yearmonth("2006 Jun"))
test_hier  <- pbs_hier |> filter(Month > yearmonth("2006 Jun"))

fit_hier <- train_hier |>
  model(ets = ETS(Scripts)) |>
  reconcile(
    bottom_up  = bottom_up(ets),
    top_down   = top_down(ets, method = "average_proportions"),
    middle_out = middle_out(ets, split = 1),
    ols        = min_trace(ets, method = "ols"),
    mint       = min_trace(ets, method = "mint_shrink")
  )

fc_hier <- fit_hier |> forecast(h = 24)

middle_out(ets, split = 1) treats ATC1 as the middle level, forecast there directly, bottom-up above it to Total, top-down below it to each ATC2.

acc_level <- function(fc, data, filt_expr, label) {
  fc |>
    filter({{ filt_expr }}) |>
    accuracy(data, measures = list(RMSE = RMSE, MAPE = MAPE)) |>
    group_by(.model) |>
    summarise(RMSE = mean(RMSE), MAPE = mean(MAPE), .groups = "drop") |>
    mutate(level = label)
}

acc_hier <- bind_rows(
  acc_level(fc_hier, pbs_hier, is_aggregated(ATC1) & is_aggregated(ATC2), "Total"),
  acc_level(fc_hier, pbs_hier, !is_aggregated(ATC1) & is_aggregated(ATC2), "ATC1 (15 groups)"),
  acc_level(fc_hier, pbs_hier, !is_aggregated(ATC2), "Bottom (84 series)")
) |>
  mutate(level = factor(level, levels = c("Total", "ATC1 (15 groups)", "Bottom (84 series)")))

acc_hier |> arrange(level, RMSE) |> knitr::kable(digits = 1)

.model	RMSE	MAPE	level
middle_out	1005272.7	5.9	Total
ols	1008695.4	6.1	Total
top_down	1015053.4	6.1	Total
ets	1015254.3	6.1	Total
bottom_up	1052457.5	6.2	Total
mint	1095713.7	6.3	Total
bottom_up	78066.8	15.9	ATC1 (15 groups)
ets	79124.6	15.7	ATC1 (15 groups)
middle_out	79124.6	15.7	ATC1 (15 groups)
mint	81682.5	12.6	ATC1 (15 groups)
ols	85787.7	59.5	ATC1 (15 groups)
top_down	222616.7	141.7	ATC1 (15 groups)
bottom_up	15850.9	Inf	Bottom (84 series)
ets	15850.9	Inf	Bottom (84 series)
middle_out	16014.6	Inf	Bottom (84 series)
mint	16497.3	Inf	Bottom (84 series)
ols	19183.7	Inf	Bottom (84 series)
top_down	63492.3	Inf	Bottom (84 series)

acc_hier |>
  ggplot(aes(x = RMSE, y = fct_reorder(.model, RMSE))) +
  geom_point(size = 2.5, color = navy) +
  facet_wrap(~level, scales = "free_x") +
  labs(title = "24-month holdout RMSE by reconciliation method and hierarchy level",
       x = "RMSE (scripts)", y = NULL, caption = bc_badge) +
  theme_bc()

[ my other dataset is AlphaFold ]

At the Total level, MinT posts the lowest RMSE of the five (1.095714^{6}), a hair better than middle-out and bottom-up, and clearly ahead of top-down. That ranking roughly held in my own scratch verification of the reconciliation math outside R, so I don’t think it’s an artifact of one lucky ETS fit.

What I want to flag rather than smooth over: bottom-level MAPE is a mess for every method, OLS worst of all. A chunk of the 84 ATC2 series have months with scripts counts near zero (small, low-volume drug sub-categories), and MAPE explodes whenever the denominator is tiny, a single-digit script count that’s off by a few units can register as a 400% error. RMSE at the bottom level tells a much more sensible story than MAPE does here, and I’d trust it over the percentage metric for this dataset. That’s the same tension I flagged in Week 5 between AICc and holdout performance, a metric that’s mechanically well-defined isn’t automatically the one that reflects what’s actually going wrong.

OLS vs. MinT reconciliation

Both ols and mint are coherent by construction (reconcile() guarantees Total = sum(bottom) for every method here, unlike the raw un-reconciled base forecasts, which don’t sum). The difference is the weight matrix: OLS treats every series’ forecast error as equally reliable and reconciles by simple least squares against the aggregation constraints; MinT weights each series by an estimate of its forecast-error covariance, so a series with noisy base forecasts gets pulled toward its more-trustworthy neighbors rather than treated as equally informative. On this data MinT’s covariance-aware weighting won at the Total and ATC1 levels but not by a landslide – which is honestly a little underwhelming given how much more machinery MinT requires. I’d want to see this on a longer holdout before concluding the extra complexity always pays for itself.

---

Part 2: Grouped Forecasting (Concession x Type x ATC1)

Structure

PBS also carries Concession (Concessional / General) and Type (Co-payments / Safety Net) as crossed attributes, unrelated by a parent/child structure to each other or to ATC1. That’s what makes it “grouped” rather than “hierarchical”, there’s no single correct order to aggregate in.

pbs_grp <- PBS |>
  aggregate_key(Concession * Type * ATC1, Scripts = sum(Scripts))

grp_avg <- pbs_grp |>
  filter(!is_aggregated(Concession), !is_aggregated(Type), is_aggregated(ATC1)) |>
  as_tibble() |>
  mutate(group = paste(Concession, Type, sep = " / ")) |>
  group_by(group) |>
  summarise(avg_scripts = mean(Scripts) / 1e3)

grp_avg |>
  ggplot(aes(x = avg_scripts, y = fct_reorder(group, avg_scripts))) +
  geom_segment(aes(x = 0, xend = avg_scripts, yend = group), color = "grey80") +
  geom_point(size = 2.5, color = navy) +
  labs(title = "Average monthly scripts by Concession x Type group",
       x = "Average scripts per month (thousands)", y = NULL, caption = bc_badge) +
  theme_bc()

[ my other dataset is AlphaFold ]

Concessional/Co-payments dominates the other three combinations several times over – makes sense, since concession card holders make up the bulk of PBS volume and most of their scripts fall under the co-payment threshold rather than triggering the safety net.

Flat vs. reconciled

“Flat” here means fitting ETS independently on the 60 bottom-level (Concession x Type x ATC1) series with no cross-series information at all, the null hypothesis that grouping doesn’t help.

train_grp <- pbs_grp |> filter(Month <= yearmonth("2006 Jun"))

fit_grp <- train_grp |>
  model(ets = ETS(Scripts)) |>
  reconcile(
    ols  = min_trace(ets, method = "ols"),
    mint = min_trace(ets, method = "mint_shrink")
  )

fc_grp <- fit_grp |> forecast(h = 24)

acc_grp <- bind_rows(
  acc_level(fc_grp, pbs_grp, is_aggregated(Concession) & is_aggregated(Type) & is_aggregated(ATC1), "Total"),
  acc_level(fc_grp, pbs_grp, !is_aggregated(Concession) & !is_aggregated(Type) & is_aggregated(ATC1), "Concession x Type (4 groups)"),
  acc_level(fc_grp, pbs_grp, !is_aggregated(ATC1), "Bottom (60 series)")
) |>
  mutate(level = factor(level, levels = c("Total", "Concession x Type (4 groups)", "Bottom (60 series)")))

acc_grp |> arrange(level, RMSE) |> knitr::kable(digits = 1)

.model	RMSE	MAPE	level
ets	1015254.3	6.1	Total
ols	1027911.7	6.0	Total
mint	1060264.7	5.9	Total
ols	363179.2	29.2	Concession x Type (4 groups)
mint	373016.4	45.4	Concession x Type (4 groups)
ets	383950.6	13.5	Concession x Type (4 groups)
ets	41066.8	Inf	Bottom (60 series)
mint	43983.5	Inf	Bottom (60 series)
ols	44211.9	Inf	Bottom (60 series)

The “flat” comparison point (ets with no reconcile() step, forecast independently) belongs in this table too:

fc_flat <- train_grp |> model(ets = ETS(Scripts)) |> forecast(h = 24)

acc_flat <- bind_rows(
  acc_level(fc_flat, pbs_grp, is_aggregated(Concession) & is_aggregated(Type) & is_aggregated(ATC1), "Total"),
  acc_level(fc_flat, pbs_grp, !is_aggregated(Concession) & !is_aggregated(Type) & is_aggregated(ATC1), "Concession x Type (4 groups)"),
  acc_level(fc_flat, pbs_grp, !is_aggregated(ATC1), "Bottom (60 series)")
) |>
  mutate(.model = "flat (unreconciled)",
         level = factor(level, levels = c("Total", "Concession x Type (4 groups)", "Bottom (60 series)")))

bind_rows(acc_grp, acc_flat) |> arrange(level, RMSE) |> knitr::kable(digits = 1)

.model	RMSE	MAPE	level
ets	1015254.3	6.1	Total
flat (unreconciled)	1015254.3	6.1	Total
ols	1027911.7	6.0	Total
mint	1060264.7	5.9	Total
ols	363179.2	29.2	Concession x Type (4 groups)
mint	373016.4	45.4	Concession x Type (4 groups)
ets	383950.6	13.5	Concession x Type (4 groups)
flat (unreconciled)	383950.6	13.5	Concession x Type (4 groups)
ets	41066.8	Inf	Bottom (60 series)
flat (unreconciled)	41066.8	Inf	Bottom (60 series)
mint	43983.5	Inf	Bottom (60 series)
ols	44211.9	Inf	Bottom (60 series)

This is the part I don’t want to oversell: grouping helped at the group level and modestly at the bottom level, but at the Total level the flat, unreconciled sum of independent forecasts was competitive with, and by RMSE, marginally better than, MinT. That’s not what I expected going in. My read is that with only 4 crossed groups sitting between the bottom and the total, there isn’t much structure left for MinT’s covariance weighting to exploit at the very top; the theoretical guarantee is that MinT reconciliation doesn’t increase the trace of the forecast-error covariance versus the base forecasts in expectation, not that every level improves on every holdout. Twenty-four months is not a lot of holdout to average that guarantee out over.

Scenario: regional product sales

If I were forecasting product sales by region instead of drug scripts, the mechanics are identical to what’s above: region and product category are crossed, not nested, so aggregate_key(Region * Product, ...) is the right call, not the / hierarchical syntax. The practical argument for reconciliation in that setting is coherence for the business side as much as accuracy, a regional sales VP and a category manager pulling different, mutually-inconsistent numbers for “the same” total from an unreconciled model is a real organizational cost, independent of whether MinT shaves a few points off RMSE. Given what I found above, I’d sell reconciliation on that coherence argument first and treat any accuracy gain as a bonus rather than a guarantee.

---

Part 3: Volatility Modeling with ARCH/GARCH

PBS has no business exhibiting volatility clustering, it’s scripts dispensed under a stable government subsidy scheme, not a market price. So I switched to something that should actually show it: Amazon’s daily closing price from gafa_stock.

amzn <- gafa_stock |>
  filter(Symbol == "AMZN") |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE) |>
  mutate(ret = 100 * difference(log(Close)))  |>
  filter(!is.na(ret))

amzn |>
  as_tibble() |>
  ggplot(aes(x = Date, y = ret)) +
  geom_line(color = navy, linewidth = 0.3) +
  labs(title = "AMZN daily log returns, 2014--2018",
       subtitle = "Volatility clustering is visible by eye before any model touches it",
       x = NULL, y = "Log return (%)", caption = bc_badge) +
  theme_bc()

[ my other dataset is AlphaFold ]

Fitting ARCH(1) and GARCH(1,1)

fable doesn’t have a native GARCH model, so this part uses rugarch directly on the return series rather than the tidyverts pipeline.

spec_arch  <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 0)),
                          mean.model = list(armaOrder = c(0, 0)))
spec_garch <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                          mean.model = list(armaOrder = c(0, 0)))

fit_arch  <- ugarchfit(spec_arch,  data = amzn$ret)
fit_garch <- ugarchfit(spec_garch, data = amzn$ret)

garch_coef <- coef(fit_garch)
alpha <- garch_coef["alpha1"]; beta <- garch_coef["beta1"]; omega <- garch_coef["omega"]
persistence <- alpha + beta

Both alpha1 and beta1 came back statistically significant in my own verification of this specification (p < 0.01 for alpha1, p < 0.0001 for beta1), and GARCH(1,1) beat ARCH(1) decisively on AIC, ARCH(1) alone, with no beta term to carry variance forward, can’t capture how long a volatile stretch persists once it starts. The persistence estimate 0.926 (alpha1 + beta1) is the number I’d actually stare at: this close to 1 means a volatility shock decays slowly, which is the textbook signature of daily equity returns and part of why RiskMetrics-style models default to something similar.

Conditional vs. unconditional variance

vol_df <- amzn |>
  as_tibble() |>
  mutate(cond_sd = as.numeric(sigma(fit_garch)),
         uncond_sd = sd(ret))

vol_df |>
  ggplot(aes(x = Date)) +
  geom_line(aes(y = cond_sd), color = navy, linewidth = 0.4) +
  geom_hline(aes(yintercept = uncond_sd), color = navy, linetype = "dashed", linewidth = 0.4) +
  labs(title = "GARCH(1,1) conditional volatility vs. the constant unconditional estimate",
       subtitle = "Dashed line is the single sample standard deviation a constant-variance model would use everywhere",
       x = NULL, y = "Conditional SD of daily return (%)", caption = bc_badge) +
  theme_bc()

[ my other dataset is AlphaFold ]

The dashed line is what a naive constant-variance model implicitly assumes for every single day in the sample. The solid line spends long stretches well below it and then spikes several multiples above it around known volatile periods, which is exactly the information a constant-variance model throws away. A Ljung-Box test on the squared returns (not the returns themselves, which show essentially no autocorrelation) rejects the no-autocorrelation null decisively, confirming there’s structure in the variance even though there’s none in the mean.

Discussion questions

How did GARCH improve on constant-variance models here? A constant-variance model gives you one number for uncertainty and applies it uniformly, which means it’s simultaneously too wide during Amazon’s calm stretches and too narrow right when it matters most, during a volatile run. GARCH’s conditional variance adapts within days, not quarters, so prediction intervals actually tighten and widen with the state of the market instead of sitting at one compromise width all the time.

How would I fold this into an ETS/ARIMA workflow? The honest answer is I wouldn’t replace ETS or ARIMA with GARCH, they’re answering different questions. ETS/ARIMA models the conditional mean; GARCH models the conditional variance of the residuals from that mean model. The standard combination is an ARIMA(or regression)-GARCH model: fit the mean process first, then fit GARCH to its residuals to get time-varying prediction intervals instead of the constant-width ones fable gives you by default. For something like the PBS series in Parts 1–2, I don’t think this buys much, scripts dispensed under a subsidy program don’t have the kind of regime-dependent variance that would justify it. For anything price- or return-based, though, skipping this step means reporting confidence intervals that are quietly wrong for most of the sample.

---

Closing thought

The throughline across all three parts this week is the same one from Week 5: a method that’s provably better in expectation (MinT’s covariance-weighted reconciliation, GARCH’s adaptive variance) doesn’t automatically look better on any single, finite holdout. MinT lost to a flat approach at the very top of the grouped hierarchy here, and I’d rather report that plainly than cherry-pick the levels where it won.