Week 6 Discussion
Matthew Gonzalez
2026-07-13
library(fpp3)## Warning: package 'fpp3' was built under R version 4.5.3## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.3 ──## ✔ tibble 3.3.1 ✔ tsibble 1.2.0
## ✔ dplyr 1.1.4 ✔ tsibbledata 0.4.1
## ✔ tidyr 1.3.2 ✔ ggtime 0.2.0
## ✔ lubridate 1.9.4 ✔ feasts 0.5.0
## ✔ ggplot2 4.0.2 ✔ fable 0.5.0## Warning: package 'tibble' was built under R version 4.5.2## Warning: package 'tidyr' was built under R version 4.5.2## Warning: package 'ggplot2' was built under R version 4.5.2## Warning: package 'tsibble' was built under R version 4.5.3## Warning: package 'tsibbledata' was built under R version 4.5.3## Warning: package 'ggtime' was built under R version 4.5.3## Warning: package 'feasts' was built under R version 4.5.3## Warning: package 'fabletools' was built under R version 4.5.3## Warning: package 'fable' was built under R version 4.5.3## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date() masks base::date()
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tsibble::setdiff() masks base::setdiff()
## ✖ tsibble::union() masks base::union()library(rugarch)## Warning: package 'rugarch' was built under R version 4.5.3## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:fabletools':
##
## reportlibrary(tidyverse)## Warning: package 'tidyverse' was built under R version 4.5.2## Warning: package 'readr' was built under R version 4.5.2## Warning: package 'purrr' was built under R version 4.5.2## Warning: package 'stringr' was built under R version 4.5.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.1 ✔ readr 2.1.6
## ✔ purrr 1.2.1 ✔ stringr 1.6.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ tsibble::interval() masks lubridate::interval()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::reduce() masks rugarch::reduce()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsPart 1
Create the Hierarchiary
tourism_h <- tourism %>%
aggregate_key(
State / Region,
Trips = sum(Trips)
)Visualize
tourism_h %>%
filter(
!is_aggregated(State),
is_aggregated(Region)
) %>%
autoplot(Trips) +
facet_wrap(vars(State), scales = "free_y") +
labs(
title = "Tourism Trips by State",
x = "Quarter",
y = "Trips"
)
Create training/testing data
tourism_train <- tourism_h %>%
filter(year(Quarter) <= 2015)
tourism_test <- tourism_h %>%
filter(year(Quarter) > 2015)fit models
tourism_fit <- tourism_train %>%
model(base = ETS(Trips)) %>%
reconcile(
bottom_up = bottom_up(base),
top_down = top_down(
base,
method = "average_proportions"
),
middle_out = middle_out(
base,
split = 1
),
ols = min_trace(
base,
method = "ols"
),
mint = min_trace(
base,
method = "mint_shrink"
)
)plot forecasts
tourism_fc <- tourism_fit %>%
forecast(h = "2 years")
tourism_fc %>%
filter(
is_aggregated(State),
is_aggregated(Region)
) %>%
autoplot(tourism_h, level = NULL) +
labs(
title = "Hierarchical Forecast Comparison",
x = "Quarter",
y = "Trips"
)
compare forecast accuracy
tourism_fc %>%
accuracy(tourism_test) %>%
group_by(.model) %>%
summarise(
RMSE = mean(RMSE, na.rm = TRUE),
MAPE = mean(MAPE, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(RMSE)## # A tibble: 6 × 3
## .model RMSE MAPE
## <chr> <dbl> <dbl>
## 1 ols 90.1 16.9
## 2 middle_out 94.7 16.1
## 3 base 96.2 16.7
## 4 mint 101. 15.9
## 5 top_down 108. 21.6
## 6 bottom_up 113. 16.8comparisons
tourism_fc %>%
accuracy(tourism_test) %>%
mutate(
Level = case_when(
is_aggregated(State) ~ "Australia",
is_aggregated(Region) ~ "State",
TRUE ~ "Region"
)
) %>%
group_by(Level, .model) %>%
summarise(
RMSE = mean(RMSE, na.rm = TRUE),
MAPE = mean(MAPE, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(Level, RMSE)## # A tibble: 18 × 4
## Level .model RMSE MAPE
## <chr> <chr> <dbl> <dbl>
## 1 Australia base 1721. 5.22
## 2 Australia top_down 1721. 5.22
## 3 Australia ols 1761. 5.37
## 4 Australia middle_out 2027. 6.48
## 5 Australia mint 2143. 7.07
## 6 Australia bottom_up 2514. 8.73
## 7 Region middle_out 47.0 16.9
## 8 Region ols 47.0 17.9
## 9 Region mint 48.7 16.7
## 10 Region base 52.6 17.5
## 11 Region bottom_up 52.6 17.5
## 12 Region top_down 61.4 22.8
## 13 State ols 291. 8.70
## 14 State middle_out 307. 9.83
## 15 State base 307. 9.83
## 16 State mint 339. 10.0
## 17 State top_down 347. 12.6
## 18 State bottom_up 389. 11.3The tourism data shows differences in both scale/movement for the Australian states. With that, most states generally increased near the end of the series. The reconciliation methods also performed differently depending on the level of hierarchy. Overall, OLS had the lowest RMSE, but MinT had the lowest MAPE. At the national level, the base/top down forecasts performed the best while middle out had the lowest regional RMSE and OLS had the lowest state-level RMSE. This shows that no single method was the unanimous best, but the reconciliation ensured that everything stayed consistent.
Part 2
Create dataset
industries <- aus_retail %>%
distinct(Industry) %>%
slice(1:3) %>%
pull(Industry)
retail_grouped <- aus_retail %>%
filter(Industry %in% industries) %>%
aggregate_key(
State * Industry,
RetailSales = sum(Turnover)
)visualize
retail_grouped %>%
filter(
!is_aggregated(State),
!is_aggregated(Industry)
) %>%
autoplot(RetailSales) +
facet_wrap(
vars(State, Industry),
scales = "free_y"
) +
labs(
title = "Retail Sales by State and Industry",
x = "Month",
y = "Retail Sales"
) +
theme(legend.position = "none")
retail_grouped %>%
filter(
!is_aggregated(State),
is_aggregated(Industry)
) %>%
autoplot(RetailSales) +
facet_wrap(vars(State), scales = "free_y") +
labs(
title = "Total Retail Sales by State",
x = "Month",
y = "Retail Sales"
) +
theme(legend.position = "none")
training sets
last_month <- max(retail_grouped$Month)
retail_train <- retail_grouped %>%
filter(Month <= last_month - 24)
retail_test <- retail_grouped %>%
filter(Month > last_month - 24)fit models
retail_fit <- retail_train %>%
model(
ETS = ETS(RetailSales)
)Reconcile
retail_reconciled <- retail_fit %>%
reconcile(
BottomUp = bottom_up(ETS),
MinTrace = min_trace(ETS, method = "mint_shrink")
)forecast
retail_forecasts <- retail_reconciled %>%
forecast(h = 24)plot
retail_forecasts %>%
filter(
!is_aggregated(State),
!is_aggregated(Industry)
) %>%
autoplot(
retail_grouped %>%
filter(Month > last_month - 60)
) +
facet_wrap(
vars(State, Industry),
scales = "free_y"
) +
labs(
title = "Reconciled Retail Sales Forecasts",
x = "Month",
y = "Retail Sales"
)
accuracy - grouped
grouped_accuracy <- retail_forecasts %>%
accuracy(retail_test) %>%
filter(
!is_aggregated(State),
!is_aggregated(Industry)
) %>%
select(.model, RMSE, MAE, MAPE)
grouped_accuracy## # A tibble: 72 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 BottomUp 2.72 2.18 5.55
## 2 BottomUp 3.33 2.74 4.10
## 3 BottomUp 2.47 2.05 9.72
## 4 BottomUp 61.7 56.2 7.47
## 5 BottomUp 29.9 21.7 1.66
## 6 BottomUp 17.0 13.0 2.48
## 7 BottomUp 2.22 1.58 8.43
## 8 BottomUp 4.30 3.76 9.71
## 9 BottomUp 0.805 0.737 8.61
## 10 BottomUp 41.1 36.2 9.34
## # ℹ 62 more rowsfit flat models
retail_flat <- aus_retail %>%
filter(Industry %in% industries) %>%
select(Month, State, Industry, Turnover)
flat_train <- retail_flat %>%
filter(Month <= last_month - 24)
flat_test <- retail_flat %>%
filter(Month > last_month - 24)
flat_fit <- flat_train %>%
model(
FlatETS = ETS(Turnover)
)
flat_forecasts <- flat_fit %>%
forecast(h = 24)evaluate flat forecast
flat_accuracy <- flat_forecasts %>%
accuracy(flat_test) %>%
select(.model, RMSE, MAE, MAPE)
flat_accuracy## # A tibble: 24 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 FlatETS 2.72 2.18 5.55
## 2 FlatETS 3.33 2.74 4.10
## 3 FlatETS 2.47 2.05 9.72
## 4 FlatETS 61.7 56.2 7.47
## 5 FlatETS 29.9 21.7 1.66
## 6 FlatETS 17.0 13.0 2.48
## 7 FlatETS 2.22 1.58 8.43
## 8 FlatETS 4.30 3.76 9.71
## 9 FlatETS 0.805 0.737 8.61
## 10 FlatETS 41.1 36.2 9.34
## # ℹ 14 more rowscompare
bind_rows(
grouped_accuracy %>%
mutate(Type = "Grouped and Reconciled"),
flat_accuracy %>%
mutate(Type = "Flat")
) %>%
group_by(Type, .model) %>%
summarise(
RMSE = mean(RMSE, na.rm = TRUE),
MAE = mean(MAE, na.rm = TRUE),
MAPE = mean(MAPE, na.rm = TRUE),
.groups = "drop"
)## # A tibble: 4 × 5
## Type .model RMSE MAE MAPE
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Flat FlatETS 16.5 13.9 7.18
## 2 Grouped and Reconciled BottomUp 16.5 13.9 7.18
## 3 Grouped and Reconciled ETS 16.5 13.9 7.18
## 4 Grouped and Reconciled MinTrace 16.6 14.2 6.97The retail data showed a strong upward trend across most states, but the size and growth of sales differed by both state and industry. The flat ETS, base ETS, and bottom-up forecasts produced the same average accuracy. MinTrace had a slightly higher RMSE and MAE but the lowest MAPE, so the grouped forecasts did not clearly outperform the flat models. In a regional sales setting, one region’s forecast can affect the reconciled totals that connect all regions and industries. Reconciliation improves reliability by adjusting the forecasts so that individual state and industry values add correctly to the overall totals, preventing conflicting forecasts across categories.
Part 3
create dax returns
data("EuStockMarkets")
dax_prices <- EuStockMarkets[, "DAX"]
dax_returns <- diff(log(dax_prices)) * 100
dax_returns <- as.numeric(na.omit(dax_returns))plot
returns_data <- tibble(
Time = seq_along(dax_returns),
Return = dax_returns
)
ggplot(returns_data, aes(x = Time, y = Return)) +
geom_line() +
labs(
title = "Daily DAX Returns",
x = "Observation",
y = "Percent Return"
)
plot squared returns
ggplot(returns_data, aes(x = Time, y = Return^2)) +
geom_line() +
labs(
title = "Squared DAX Returns",
x = "Observation",
y = "Squared Return"
)
specify model
garch_spec <- ugarchspec(
variance.model = list(
model = "sGARCH",
garchOrder = c(1, 1)
),
mean.model = list(
armaOrder = c(0, 0),
include.mean = TRUE
),
distribution.model = "norm"
)fit model
garch_fit <- ugarchfit(
spec = garch_spec,
data = dax_returns
)
garch_fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(0,0,0)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.065353 0.021576 3.0290 0.002454
## omega 0.047563 0.012813 3.7121 0.000206
## alpha1 0.068454 0.014975 4.5713 0.000005
## beta1 0.887569 0.023897 37.1417 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.065353 0.022151 2.9503 0.003175
## omega 0.047563 0.034132 1.3935 0.163465
## alpha1 0.068454 0.025102 2.7270 0.006391
## beta1 0.887569 0.045481 19.5153 0.000000
##
## LogLikelihood : -2594.796
##
## Information Criteria
## ------------------------------------
##
## Akaike 2.7959
## Bayes 2.8078
## Shibata 2.7959
## Hannan-Quinn 2.8003
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1997 0.6550
## Lag[2*(p+q)+(p+q)-1][2] 0.3475 0.7699
## Lag[4*(p+q)+(p+q)-1][5] 0.7977 0.9035
## d.o.f=0
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.1255 0.7231
## Lag[2*(p+q)+(p+q)-1][5] 0.3366 0.9797
## Lag[4*(p+q)+(p+q)-1][9] 0.4989 0.9985
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.01578 0.500 2.000 0.9000
## ARCH Lag[5] 0.32054 1.440 1.667 0.9348
## ARCH Lag[7] 0.37485 2.315 1.543 0.9883
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 2.5476
## Individual Statistics:
## mu 0.7035
## omega 0.1271
## alpha1 0.2864
## beta1 0.1054
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.07 1.24 1.6
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.4183 0.1563
## Negative Sign Bias 0.8066 0.4200
## Positive Sign Bias 0.4291 0.6679
## Joint Effect 4.2421 0.2365
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 97.12 1.776e-12
## 2 30 132.02 4.270e-15
## 3 40 130.52 8.456e-12
## 4 50 195.57 2.093e-19
##
##
## Elapsed time : 0.127902view params
coef(garch_fit)## mu omega alpha1 beta1
## 0.06535253 0.04756287 0.06845367 0.88756875garch_fit@fit$matcoef## Estimate Std. Error t value Pr(>|t|)
## mu 0.06535253 0.02157585 3.028967 2.453911e-03
## omega 0.04756287 0.01281283 3.712129 2.055231e-04
## alpha1 0.06845367 0.01497478 4.571263 4.847931e-06
## beta1 0.88756875 0.02389685 37.141671 0.000000e+00calc volatility
parameters <- coef(garch_fit)
persistence <- parameters["alpha1"] + parameters["beta1"]
persistence## alpha1
## 0.9560224compare the variance
conditional_variance <- sigma(garch_fit)^2
unconditional_variance <- var(dax_returns)
variance_data <- tibble(
Time = seq_along(conditional_variance),
ConditionalVariance = conditional_variance,
UnconditionalVariance = unconditional_variance
)
ggplot(variance_data, aes(x = Time)) +
geom_line(aes(y = ConditionalVariance)) +
geom_line(
aes(y = UnconditionalVariance),
linetype = "dashed"
) +
labs(
title = "Conditional and Unconditional Variance",
x = "Observation",
y = "Variance"
)## Don't know how to automatically pick scale for object of type <xts/zoo>.
## Defaulting to continuous.
forecast
garch_forecast <- ugarchforecast(
garch_fit,
n.ahead = 20
)
garch_forecast##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: sGARCH
## Horizon: 20
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=1975-02-03]:
## Series Sigma
## T+1 0.06535 1.527
## T+2 0.06535 1.509
## T+3 0.06535 1.491
## T+4 0.06535 1.475
## T+5 0.06535 1.458
## T+6 0.06535 1.442
## T+7 0.06535 1.427
## T+8 0.06535 1.412
## T+9 0.06535 1.398
## T+10 0.06535 1.384
## T+11 0.06535 1.371
## T+12 0.06535 1.358
## T+13 0.06535 1.346
## T+14 0.06535 1.334
## T+15 0.06535 1.322
## T+16 0.06535 1.311
## T+17 0.06535 1.300
## T+18 0.06535 1.290
## T+19 0.06535 1.280
## T+20 0.06535 1.270plot volatility
forecast_variance <- as.numeric(
sigma(garch_forecast)^2
)
forecast_data <- tibble(
Day = 1:20,
ForecastVariance = forecast_variance
)
ggplot(
forecast_data,
aes(x = Day, y = ForecastVariance)
) +
geom_line() +
labs(
title = "GARCH Variance Forecast",
x = "Forecast Day",
y = "Forecast Conditional Variance"
)
The DAX returns clearly show volatility clustering because periods of relatively small movements are followed by periods with much larger positive and negative returns. The GARCH model captured this changing variance better than a constant-variance model. The alpha estimate shows that recent shocks affected current volatility, while the beta estimate shows that volatility was highly persistent and faded slowly over time. The Ljung-Box and ARCH tests had large p-values, indicating that little serial correlation or remaining ARCH behavior was left after fitting the model. The forecast also shows conditional volatility gradually declining toward its longer-run level. In a broader forecasting workflow, ETS or ARIMA could model the expected value of the series, while GARCH could model changing uncertainty and produce more realistic prediction intervals.