Week 7 discussion
Part 1: Hierarchical Time Series Forecasting
Exploring Hierarchical Structures
The tourism dataset records quarterly Australian domestic trips from 1998 Q1 to 2017 Q4, with a clean nested hierarchy: Total -> State -> Region. Every region belongs to exactly one state, so summing regions gives states, and summing states gives the national total.
tdf <- as.data.frame(tourism) %>% select(Quarter, Region, State, Purpose, Trips)
tdf$Quarter <- as.numeric(unclass(tdf$Quarter))
agg <- tdf %>% group_by(Quarter, State, Region) %>% summarise(Trips = sum(Trips), .groups = "drop")
quarters <- sort(unique(agg$Quarter))
n_q <- length(quarters)
qtr_dates <- seq(as.Date("1998-01-01"), by = "quarter", length.out = n_q)
region_wide <- agg %>% select(Quarter, Region, Trips) %>%
pivot_wider(names_from = Region, values_from = Trips) %>% arrange(Quarter)
region_names <- setdiff(names(region_wide), "Quarter")
n_bottom <- length(region_names)
region_state_map <- agg %>% distinct(Region, State)
states <- sort(unique(region_state_map$State))
n_state <- length(states)The hierarchy has 8 states and 76 regions across 80 quarters of data.
region_state_map %>%
count(State, name = "n_regions") %>%
ggplot(aes(x = reorder(State, n_regions), y = n_regions)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = n_regions), hjust = -0.3, size = 3.2) +
coord_flip() +
labs(title = "Tourism Hierarchy: Number of Regions per State", x = NULL, y = "Number of Regions")
Number of regions per state
tourism %>% summarise(Trips = sum(Trips)) %>%
autoplot(Trips) +
labs(title = "Top of the Hierarchy: Total Australian Domestic Trips", y = "Trips ('000s)", x = "Quarter")
Total trips
tourism %>% group_by(State) %>% summarise(Trips = sum(Trips)) %>%
autoplot(Trips) +
labs(title = "Middle of the Hierarchy: Trips by State", y = "Trips ('000s)", x = "Quarter", colour = "State")
Trips by state
The hierarchy is unbalanced. Victoria alone contains 21 of the 76 regions, while the ACT is a single region. This matters for reconciliation: methods that pool information across many regions (Bottom-Up, MinT) have more to work with in states like Victoria or NSW than in the ACT. The total and state series share the same seasonal pattern and a trend that dips after 2008 and accelerates from 2013 onward, this shared structure is exactly what reconciliation exploits, since independently-forecast series will not, in general, sum back up to a sensible total.
Modeling and Forecasting Hierarchical Data
An independent ETS model is fit to every series in the hierarchy,
holding out the last 2 years (8 quarters) for testing. From these base
forecasts, forecasts are built four ways Bottom-Up, Top-Down
(historical-proportion split), Middle-Out (forecast States, aggregate
up, split down) and two reconciliation methods, OLS and MinT (diagonal),
computed from the summing matrix S.
h <- 8
train_idx <- 1:(n_q - h); test_idx <- (n_q - h + 1):n_q
region_mat <- as.matrix(region_wide[, region_names])
train_region <- region_mat[train_idx, , drop = FALSE]
test_region <- region_mat[test_idx, , drop = FALSE]
state_mat <- sapply(states, function(s) {
rs <- region_state_map$Region[region_state_map$State == s]
rowSums(region_mat[, rs, drop = FALSE])
})
total_vec <- rowSums(state_mat)
train_state <- state_mat[train_idx, , drop = FALSE]; test_state <- state_mat[test_idx, , drop = FALSE]
train_total <- total_vec[train_idx]; test_total <- total_vec[test_idx]
make_ts <- function(x) ts(x, start = c(1998, 1), frequency = 4)
fit_fc <- function(x_train) {
fit <- ets(make_ts(x_train))
list(fc = as.numeric(forecast(fit, h = h)$mean), resid = as.numeric(residuals(fit)))
}
fit_total <- fit_fc(train_total)
fit_states <- setNames(lapply(states, function(s) fit_fc(train_state[, s])), states)
fit_regions <- setNames(lapply(region_names, function(r) fit_fc(train_region[, r])), region_names)
base_total <- fit_total$fc
base_state <- sapply(fit_states, function(z) z$fc)
base_region <- sapply(fit_regions, function(z) z$fc)
# Bottom-Up
bu_state <- sapply(states, function(s) {
rs <- region_state_map$Region[region_state_map$State == s]
rowSums(base_region[, rs, drop = FALSE])
})
bu_total <- rowSums(bu_state)
# Top-Down (historical proportions)
prop_state <- colSums(train_state) / sum(train_total)
prop_region <- colSums(train_region) / sum(train_total)
td_state <- outer(base_total, prop_state); colnames(td_state) <- states
td_region <- outer(base_total, prop_region); colnames(td_region) <- region_names
# Middle-Out
mo_total <- rowSums(base_state)
prop_region_in_state <- sapply(region_names, function(r) {
s <- region_state_map$State[region_state_map$Region == r]
sum(train_region[, r]) / sum(train_state[, s])
})
mo_region <- sapply(region_names, function(r) {
s <- region_state_map$State[region_state_map$Region == r]
base_state[, s] * prop_region_in_state[r]
})
colnames(mo_region) <- region_names
# OLS / MinT reconciliation
S_state_rows <- sapply(states, function(s) {
rs <- region_state_map$Region[region_state_map$State == s]
as.numeric(region_names %in% rs)
})
S <- rbind(rep(1, n_bottom), t(S_state_rows), diag(n_bottom))
rownames(S) <- c("Total", states, region_names)
base_all <- cbind(Total = base_total, base_state, base_region)
resid_all <- c(list(Total = fit_total$resid),
lapply(fit_states, function(z) z$resid),
lapply(fit_regions, function(z) z$resid))
min_len <- min(sapply(resid_all, length))
resid_mat <- sapply(resid_all, function(r) tail(r, min_len))
reconcile <- function(base_h, W_diag = NULL) {
if (is.null(W_diag)) { P <- solve(t(S) %*% S) %*% t(S) }
else { Winv <- diag(1 / W_diag); P <- solve(t(S) %*% Winv %*% S) %*% t(S) %*% Winv }
t(S %*% P %*% t(base_h))
}
ols_all <- reconcile(base_all)
mint_all <- reconcile(base_all, W_diag = apply(resid_mat, 2, var))
extract <- function(M) list(total = M[, "Total"], state = M[, states], region = M[, region_names])
ols <- extract(ols_all); mint <- extract(mint_all)acc_row <- function(name, total_f, state_f, region_f) {
data.frame(Method = name,
RMSE_Total = rmse(test_total, total_f), MAPE_Total = mape(test_total, total_f),
RMSE_State = mean(sapply(states, function(s) rmse(test_state[, s], state_f[, s]))),
MAPE_State = mean(sapply(states, function(s) mape(test_state[, s], state_f[, s]))),
RMSE_Region = mean(sapply(region_names, function(r) rmse(test_region[, r], region_f[, r]))),
MAPE_Region = mean(sapply(region_names, function(r) mape(test_region[, r], region_f[, r])))
)
}
results1 <- rbind(
acc_row("Base ETS (incoherent)", base_total, base_state, base_region),
acc_row("Bottom-Up", bu_total, bu_state, base_region),
acc_row("Top-Down", base_total, td_state, td_region),
acc_row("Middle-Out", mo_total, base_state, mo_region),
acc_row("OLS Reconciliation", ols$total, ols$state, ols$region),
acc_row("MinT Reconciliation", mint$total, mint$state, mint$region)
)
results1[,-1] <- round(results1[,-1], 1)
results1| Method | RMSE_Total | MAPE_Total | RMSE_State | MAPE_State | RMSE_Region | MAPE_Region |
|---|---|---|---|---|---|---|
| Base ETS (incoherent) | 1720.7 | 5.2 | 306.8 | 9.8 | 52.6 | 17.5 |
| Bottom-Up | 2513.5 | 8.7 | 388.6 | 11.3 | 52.6 | 17.5 |
| Top-Down | 1720.7 | 5.2 | 344.9 | 12.5 | 61.3 | 22.7 |
| Middle-Out | 2027.0 | 6.5 | 306.8 | 9.8 | 60.6 | 19.3 |
| OLS Reconciliation | 1760.5 | 5.4 | 290.9 | 8.7 | 47.0 | 17.9 |
| MinT Reconciliation | 1720.7 | 5.2 | 274.1 | 8.9 | 49.0 | 19.4 |
fc_dates <- qtr_dates[test_idx]
data.frame(
Date = rep(fc_dates, 6),
Trips = c(test_total, base_total, bu_total, base_total, mo_total, mint$total),
Series = rep(c("Actual","Base ETS","Bottom-Up","Top-Down","Middle-Out","MinT"), each = h)
) %>%
ggplot(aes(Date, Trips, colour = Series, linetype = Series)) +
geom_line(linewidth = 0.9) +
labs(title = "Part 1: Total-Level Forecasts vs Actual (last 8 quarters)", y = "Trips ('000s)", x = NULL)
Total-level forecasts vs actual
method_abbrev <- c("Base ETS (incoherent)"="Base ETS","Bottom-Up"="Bottom-Up","Top-Down"="Top-Down",
"Middle-Out"="Middle-Out","OLS Reconciliation"="OLS","MinT Reconciliation"="MinT")
results1 %>%
mutate(Method = recode(Method, !!!method_abbrev),
Method = factor(Method, levels = rev(unname(method_abbrev)))) %>%
select(Method, RMSE_Total, RMSE_State, RMSE_Region) %>%
pivot_longer(-Method, names_to = "Level", values_to = "RMSE") %>%
mutate(Level = factor(sub("RMSE_", "", Level), levels = c("Total","State","Region"))) %>%
ggplot(aes(Method, RMSE, fill = Method)) +
geom_col(show.legend = FALSE) +
geom_text(aes(label = round(RMSE, 0)), hjust = -0.15, size = 3) +
facet_wrap(~ Level, scales = "free_x") +
coord_flip() + expand_limits(y = 0) +
labs(title = "Part 1: RMSE by Method and Hierarchy Level", x = NULL) +
theme(strip.text = element_text(face = "bold"))
RMSE by method and level
MinT and Top-Down tie the Base ETS model at the Total level (1720.7 RMSE, 5.22% MAPE), which makes sense because both essentially preserve the top-level ETS fit rather than degrading it with noisier bottom-level information. Reconciliation is further down: MinT wins at the State level (274.1 vs 306.8 for the unreconciled base, roughly an 11% RMSE reduction), and OLS wins at the Region level (47.0 vs 52.6, also roughly an 11% reduction). Bottom-Up is worst almost everywhere (2513.5 Total RMSE, 8.72% MAPE). Summing 76 independently-noisy region forecasts loses badly to modeling the total directly, confirming that bottom-level series are usually too noisy to trust in isolation.
Discussing Reconciliation Techniques
Three traditional strategies were tested before formal reconciliation: Bottom-Up forecasts every region independently and sums upward, Top-Down forecasts only the Total and splits it down by historical share, Middle-Out forecasts States directly and both aggregates up and disaggregates down (a compromise). OLS reconciliation finds the coherent forecast set closest, in a least-squares sense, to the independent base forecasts, treating every series equally. MinT (diagonal) does the same but weights each series by the inverse of its own forecast-error variance, so noisier series are trusted less and pulled more toward the better-behaved aggregate. Both guarantee forecasts sum correctly across the hierarchy which is something none of Bottom-Up, Top-Down, or Middle-Out achieve at every level simultaneously which is why reconciliation balances accuracy by borrowing strength from correlated series against coherence (numbers that always add up), rather than treating them as separate goals.
Part 2: Grouped Time Series Forecasting
Exploring Grouped Data Structures
aus_retail is restricted to 4 states x 4 industries (16 bottom-level groups), monthly turnover. Unlike Part 1, State and Industry are non-nested grouping dimensions meaning no single tree describes the data. The summing matrix instead encodes two independent sets of aggregation constraints.
rdf <- as.data.frame(aus_retail) %>% select(Month, State, Industry, Turnover)
rdf$Month <- as.numeric(unclass(rdf$Month))
sel_states <- c("New South Wales", "Victoria", "Queensland", "Western Australia")
sel_ind <- c("Cafes, restaurants and catering services", "Clothing retailing",
"Supermarket and grocery stores", "Newspaper and book retailing")
sub <- rdf %>% filter(State %in% sel_states, Industry %in% sel_ind)
months <- sort(unique(sub$Month)); n_m <- length(months)
month_dates <- seq(as.Date("1982-04-01"), by = "month", length.out = n_m)sub %>% group_by(State, Industry) %>% summarise(AvgTurnover = mean(Turnover), .groups = "drop") %>%
ggplot(aes(Industry, State, fill = AvgTurnover)) +
geom_tile() + scale_fill_viridis_c() +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(title = "Grouped Structure: Average Turnover by State x Industry", fill = "Avg Turnover ($m)", x = NULL, y = NULL)
Average turnover heatmap
sub %>% ggplot(aes(Month, Turnover, colour = Industry)) +
geom_line() + facet_wrap(~ State, scales = "free_y") +
labs(title = "Grouped Series: Turnover over Time by State x Industry", x = "Month index (1982 M4 = 1)", y = "Turnover ($m)")
Grouped series over time
Turnover levels differ a lot by industry (supermarkets dominate everywhere) and by state (NSW and Victoria are the largest markets), with neither dimension nested inside the other. Supermarket turnover shows a much stronger upward trend in every state than the other three industries.
Modeling and Forecasting Grouped Data
sub$Group <- paste(sub$State, sub$Industry, sep = " || ")
wide <- sub %>% select(Month, Group, Turnover) %>% pivot_wider(names_from = Group, values_from = Turnover) %>% arrange(Month)
group_names <- setdiff(names(wide), "Month"); n_bottom2 <- length(group_names)
bottom_mat <- as.matrix(wide[, group_names])
group_key <- as.data.frame(do.call(rbind, strsplit(group_names, " \\|\\| ")))
colnames(group_key) <- c("State", "Industry")
state_mat2 <- sapply(sel_states, function(s) rowSums(bottom_mat[, group_key$State == s, drop = FALSE]))
ind_mat2 <- sapply(sel_ind, function(i) rowSums(bottom_mat[, group_key$Industry == i, drop = FALSE]))
total_vec2 <- rowSums(bottom_mat)
h2 <- 24
train_idx2 <- 1:(n_m - h2); test_idx2 <- (n_m - h2 + 1):n_m
make_ts2 <- function(x) ts(x, start = c(1982, 4), frequency = 12)
fit_series <- function(x_train, method = "ets") {
ts_x <- window(make_ts2(x_train), start = c(2000, 1))
fit <- if (method == "ets") ets(ts_x) else auto.arima(ts_x)
list(fc = as.numeric(forecast(fit, h = h2)$mean), resid = as.numeric(residuals(fit)))
}
train_bottom <- bottom_mat[train_idx2, , drop = FALSE]; test_bottom <- bottom_mat[test_idx2, , drop = FALSE]
train_state2 <- state_mat2[train_idx2, , drop = FALSE]; test_state2 <- state_mat2[test_idx2, , drop = FALSE]
train_ind2 <- ind_mat2[train_idx2, , drop = FALSE]; test_ind2 <- ind_mat2[test_idx2, , drop = FALSE]
train_total2 <- total_vec2[train_idx2]; test_total2 <- total_vec2[test_idx2]
fit_bottom_ets <- setNames(lapply(group_names, function(g) fit_series(train_bottom[, g], "ets")), group_names)
fit_bottom_arima <- setNames(lapply(group_names, function(g) fit_series(train_bottom[, g], "arima")), group_names)
base_bottom_ets <- sapply(fit_bottom_ets, function(z) z$fc)
base_bottom_arima <- sapply(fit_bottom_arima, function(z) z$fc)
gr_state <- sapply(sel_states, function(s) rowSums(base_bottom_ets[, group_key$State == s, drop = FALSE]))
gr_ind <- sapply(sel_ind, function(i) rowSums(base_bottom_ets[, group_key$Industry == i, drop = FALSE]))
gr_total <- rowSums(base_bottom_ets)
gr_state_ar <- sapply(sel_states, function(s) rowSums(base_bottom_arima[, group_key$State == s, drop = FALSE]))
gr_ind_ar <- sapply(sel_ind, function(i) rowSums(base_bottom_arima[, group_key$Industry == i, drop = FALSE]))
gr_total_ar <- rowSums(base_bottom_arima)
flat_state <- sapply(sel_states, function(s) fit_series(train_state2[, s], "ets")$fc)
flat_ind <- sapply(sel_ind, function(i) fit_series(train_ind2[, i], "ets")$fc)
flat_total <- fit_series(train_total2, "ets")$fc
S2 <- rbind(rep(1, n_bottom2),
t(sapply(sel_states, function(s) as.numeric(group_key$State == s))),
t(sapply(sel_ind, function(i) as.numeric(group_key$Industry == i))),
diag(n_bottom2))
rownames(S2) <- c("Total", sel_states, sel_ind, group_names)
base_all2 <- cbind(Total = flat_total, flat_state, flat_ind, base_bottom_ets)
resid_all2 <- c(list(Total = fit_series(train_total2, "ets")$resid),
lapply(sel_states, function(s) fit_series(train_state2[, s], "ets")$resid),
lapply(sel_ind, function(i) fit_series(train_ind2[, i], "ets")$resid),
lapply(fit_bottom_ets, function(z) z$resid))
min_len2 <- min(sapply(resid_all2, length))
resid_mat2 <- sapply(resid_all2, function(r) tail(r, min_len2))
reconcile2 <- function(base_h, W_diag = NULL) {
if (is.null(W_diag)) { P <- solve(t(S2) %*% S2) %*% t(S2) }
else { Winv <- diag(1 / W_diag); P <- solve(t(S2) %*% Winv %*% S2) %*% t(S2) %*% Winv }
t(S2 %*% P %*% t(base_h))
}
ols_all2 <- reconcile2(base_all2)
mint_all2 <- reconcile2(base_all2, W_diag = apply(resid_mat2, 2, var))
get_cols <- function(M, nms) M[, nms, drop = FALSE]
ols_state2 <- get_cols(ols_all2, sel_states); ols_ind2 <- get_cols(ols_all2, sel_ind); ols_total2 <- ols_all2[, "Total"]
mint_state2 <- get_cols(mint_all2, sel_states); mint_ind2 <- get_cols(mint_all2, sel_ind); mint_total2 <- mint_all2[, "Total"]acc2 <- function(name, total_f, state_f, ind_f, bottom_f) {
data.frame(Method = name,
RMSE_Total = rmse(test_total2, total_f), MAPE_Total = mape(test_total2, total_f),
RMSE_State = mean(sapply(sel_states, function(s) rmse(test_state2[, s], state_f[, s]))),
MAPE_State = mean(sapply(sel_states, function(s) mape(test_state2[, s], state_f[, s]))),
RMSE_Industry = mean(sapply(sel_ind, function(i) rmse(test_ind2[, i], ind_f[, i]))),
MAPE_Industry = mean(sapply(sel_ind, function(i) mape(test_ind2[, i], ind_f[, i]))),
RMSE_Bottom = mean(sapply(group_names, function(g) rmse(test_bottom[, g], bottom_f[, g]))),
MAPE_Bottom = mean(sapply(group_names, function(g) mape(test_bottom[, g], bottom_f[, g])))
)
}
results2 <- rbind(
acc2("Flat ETS (non-grouped)", flat_total, flat_state, flat_ind, base_bottom_ets),
acc2("Grouped Bottom-Up (ETS)", gr_total, gr_state, gr_ind, base_bottom_ets),
acc2("Grouped Bottom-Up (ARIMA)", gr_total_ar, gr_state_ar, gr_ind_ar, base_bottom_arima),
acc2("OLS Reconciliation", ols_total2, ols_state2, ols_ind2, base_bottom_ets),
acc2("MinT Reconciliation", mint_total2, mint_state2, mint_ind2, base_bottom_ets)
)
results2[,-1] <- round(results2[,-1], 2)
results2| Method | RMSE_Total | MAPE_Total | RMSE_State | MAPE_State | RMSE_Industry | MAPE_Industry | RMSE_Bottom | MAPE_Bottom |
|---|---|---|---|---|---|---|---|---|
| Flat ETS (non-grouped) | 148.39 | 1.16 | 61.44 | 1.81 | 55.82 | 4.55 | 25.07 | 6.44 |
| Grouped Bottom-Up (ETS) | 140.24 | 1.06 | 66.55 | 2.02 | 64.55 | 4.94 | 25.07 | 6.44 |
| Grouped Bottom-Up (ARIMA) | 253.97 | 2.00 | 88.56 | 3.37 | 74.60 | 4.60 | 31.33 | 7.57 |
| OLS Reconciliation | 126.95 | 0.98 | 60.78 | 1.96 | 60.05 | 7.09 | 25.07 | 6.44 |
| MinT Reconciliation | 138.31 | 1.08 | 63.72 | 2.13 | 69.12 | 10.13 | 25.07 | 6.44 |
fc_dates2 <- month_dates[test_idx2]
data.frame(
Date = rep(fc_dates2, 4),
Turnover = c(test_total2, flat_total, gr_total, mint_total2),
Series = rep(c("Actual","Flat ETS","Grouped Bottom-Up","MinT"), each = h2)
) %>%
ggplot(aes(Date, Turnover, colour = Series, linetype = Series)) +
geom_line(linewidth = 0.9) +
labs(title = "Part 2: Total-Level Forecasts vs Actual (last 24 months)", y = "Turnover ($m)", x = NULL)
Total-level forecasts vs actual
results2 %>%
select(Method, RMSE_Total, RMSE_State, RMSE_Industry, RMSE_Bottom) %>%
pivot_longer(-Method, names_to = "Level", values_to = "RMSE") %>%
mutate(Level = factor(sub("RMSE_", "", Level), levels = c("Total","State","Industry","Bottom"))) %>%
ggplot(aes(Method, RMSE, fill = Method)) +
geom_col(show.legend = FALSE) +
facet_wrap(~ Level, scales = "free_x") +
coord_flip() +
labs(title = "Part 2: RMSE by Method and Group Level", x = NULL)
RMSE by method and group level
OLS reconciliation gives the best Total-level accuracy (126.9 RMSE, 0.98% MAPE, versus 148.4 RMSE / 1.17% MAPE for the Flat model which is roughly a 14% RMSE improvement) and ties for best at the State level (60.8 vs 61.4). The Flat model still wins at the Industry level (55.8 RMSE, 4.55% MAPE), and its MAPE advantage there is larger still, since OLS and especially MinT’s industry-level MAPE balloons (7.09% and 10.13%, respectively) even though RMSE is comparable these methods appear to get pulled off on the smaller-turnover industries. ETS clearly beats ARIMA at the bottom level (25.1 vs 31.3 RMSE). The ARIMA-bottom variant is worst across nearly every level (254.0 Total RMSE, 2.00% MAPE). Because OLS, MinT, and Grouped Bottom-Up (ETS) all reconcile from the same ETS bottom-level forecasts, they share an identical Bottom-level RMSE (25.1) by construction reconciliation only reallocates forecasts at the aggregate levels, it doesn’t change the bottom-level inputs. Grouped/reconciled forecasting outperforms the flat model where it matters most (the Total), but doesn’t uniformly dominate at every level. Reconciliation trades some accuracy at certain levels for coherence and improved accuracy elsewhere.
Scenario Analysis
In a scenario where you are forecasting product sales across regions and product categories instead of retail turnover. A demand shock in one region-category cell (a local promotion, a supply disruption) shows up as a forecast error at the bottom level. Left unreconciled, that error propagates inconsistently. The region’s total tells one story while the product-category total tells another. MinT-style reconciliation pools information across all related series, so a noisy or biased bottom-level forecast gets partially corrected using the more stable aggregate signal, and the corrected forecasts are guaranteed to add up consistently across both the regional and product hierarchies which is exactly what’s needed for downstream decisions like regional inventory allocation that depend on multiple categories agreeing with each other.
Part 3: Volatility Modeling with ARCH and GARCH
Exploring Volatility Clustering
library(rugarch)
library(quantmod)
getSymbols("AAPL", src = "yahoo", from = "2014-01-01", to = "2018-12-31")## [1] "AAPL"px <- Ad(AAPL)
ret <- na.omit(diff(log(px)) * 100)
uncond_var <- var(as.numeric(ret))Number of return observations: 1256. Mean return: 0.0615%, unconditional SD: 1.5081%.
data.frame(Date = index(px), Price = as.numeric(px)) %>%
ggplot(aes(Date, Price)) + geom_line(colour = "steelblue") +
labs(title = "AAPL Adjusted Close, 2014-2018", y = "Price ($)", x = NULL)
AAPL adjusted close
data.frame(Date = index(ret), Return = as.numeric(ret)) %>%
ggplot(aes(Date, Return)) + geom_line(colour = "darkred", linewidth = 0.3) +
labs(title = "AAPL Daily Log Returns (%) - Volatility Clustering",
subtitle = "Note the alternating calm and turbulent periods", y = "Return (%)", x = NULL)
AAPL daily returns
forecast::ggAcf(as.numeric(ret)^2, lag.max = 30) +
labs(title = "ACF of Squared Returns", subtitle = "Significant autocorrelation = evidence of ARCH effects")
ACF of squared returns
Returns show classic volatility clustering: large moves cluster together for example 2014, mid-2015, early 2016 while calmer stretches persist elsewhere, even though the returns themselves show little autocorrelation. The ACF of squared returns confirms this, significant autocorrelation at many lags is direct evidence of ARCH effects and time-varying variance, motivating a GARCH-family model rather than a constant-variance assumption.
Modeling Conditional Variance
spec_arch1 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 0)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE))
fit_arch1 <- ugarchfit(spec = spec_arch1, data = ret)
spec_garch11 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(0, 0), include.mean = TRUE))
fit_garch11 <- ugarchfit(spec = spec_garch11, data = ret)
coefs <- coef(fit_garch11)
persistence <- as.numeric(coefs["alpha1"] + coefs["beta1"])
cond_sd <- sigma(fit_garch11)
aic_arch1 <- infocriteria(fit_arch1)["Akaike",]
aic_garch11 <- infocriteria(fit_garch11)["Akaike",]| Model | mu | omega | alpha1 | beta1 | AIC |
|---|---|---|---|---|---|
| ARCH(1) | 0.300 | 4.472 | 0.0000 (not sig., p≈1.0) | – | 3.862 |
| GARCH(1,1) | 0.121 | 0.256 | 0.125 | 0.769 | 3.583 |
All GARCH(1,1) coefficients are statistically significant (p < 0.001 for alpha1 and beta1). Persistence (alpha1 + beta1) = 0.894. Conditional SD range: 1.078% - 3.284%, versus unconditional SD 1.508%.
mu <- coefs["mu"]
data.frame(Date = index(ret), Dev = as.numeric(ret) - mu, SD = as.numeric(cond_sd)) %>%
ggplot(aes(Date)) +
geom_line(aes(y = Dev), colour = "steelblue", linewidth = 0.3) +
geom_line(aes(y = SD), colour = "darkred", linewidth = 0.6) +
geom_line(aes(y = -SD), colour = "darkred", linewidth = 0.6) +
geom_hline(yintercept = 0) +
labs(title = "AAPL Daily Returns with +/- 1 Conditional SD Bands (GARCH(1,1))", y = "Percent", x = NULL)
Returns with conditional SD bands
data.frame(Date = index(ret), CondSD = as.numeric(cond_sd)) %>%
ggplot(aes(Date, CondSD)) +
geom_line(colour = "darkred") +
geom_hline(yintercept = sqrt(uncond_var), linetype = "dashed") +
labs(title = "GARCH(1,1) Conditional SD vs Unconditional SD",
subtitle = "Dashed line = constant unconditional SD", y = "Conditional SD (%)", x = NULL)
Conditional vs unconditional SD
ARCH(1) essentially fails – its alpha1 coefficient collapses to 0.0000 and is statistically insignificant (p ≈ 1.0), meaning a single lagged squared-shock term cannot explain the variance dynamics on its own. GARCH(1,1) fixes this decisively: both alpha1 (0.125) and beta1 (0.769) are highly significant, and AIC drops from 3.86 to 3.58. Persistence of 0.894 (close to but below 1) means volatility shocks decay slowly, a turbulent day raises expected volatility for weeks, not just the next day. The conditional-SD plot makes the payoff obvious: it swings well above the flat 1.51% unconditional line during turbulent stretches (peaking near 3.28%) and drops as low as 1.08% during calm periods – a single constant-variance number represents neither regime well.
Discussion Questions
1. How did ARCH/GARCH improve understanding of volatility compared to constant-variance models? A constant-variance model gives one number for typical risk over the whole sample. GARCH(1,1) instead produces a full time series of risk estimates that rises ahead of and during noisy periods and falls during calm ones, showing exactly when and how much risk increased rather than averaging it away.
2. How would volatility modeling fit into a forecasting workflow alongside ETS/ARIMA? ETS/ARIMA model the conditional mean. GARCH models the conditional variance of the residuals from that mean model. In practice, fit ARIMA/ETS for the point forecast, then fit GARCH on its residuals to get time-varying prediction intervals, point forecasts stay the same, but the uncertainty bands widen and narrow with market conditions, which is far more realistic than the fixed-width intervals ETS/ARIMA produce alone.
Reflection
Of the three approaches, hierarchical reconciliation produced the clearest, most consistent win: MinT and OLS improved on the incoherent base forecasts at both the State and Region levels while guaranteeing coherence, and the nested tree structure of tourism is exactly the setting reconciliation was designed for. Grouped reconciliation was more of a mixed bag, it won at the Total level but lost to the simpler flat model at the Industry level, reflecting that non-nested grouped structures give reconciliation less to work with than a clean hierarchy. Volatility modeling answered a fundamentally different question, not how do we get more accurate point forecast but how much can we trust today’s forecast and GARCH(1,1)’s clear AIC and diagnostic improvement over both ARCH(1) and a constant-variance assumption made it the most decisive result of the three in its own domain. In a real forecasting workflow these aren’t competitors: reconciliation improves point-forecast coherence and accuracy across a hierarchy, while GARCH adds risk-aware uncertainty bands around whatever point-forecast method is used.