Week 6 Discussion
Juan Zabal Alonso
library(fpp3)
library(tseries) 1 Part 1: ARIMA Models
For this week’s discussion I decided to use aus_production which shows the quarterly Australian beer production. The series is from 1956 Q1 to 2010 Q2, n = 218.
beer <- aus_production %>%
select(Quarter, Beer) %>%
filter(!is.na(Beer))1.1 Data Visualization
beer %>%
autoplot(Beer) +
labs(title = "Australian Quarterly Beer Production",
y = "Megalitres")
By looking at the graph you can observe a rising trend through the 1970s-80s, a peak around 1990, then a gradual decline plus a strong, stable quarterly seasonal pattern. You can clearly see trend and seasonality are present here. This tells us we will need regular and seasonal differencing.
beer %>% gg_season(Beer) +
labs(title = "Seasonal Plot: Beer Production")
beer %>% gg_subseries(Beer) +
labs(title = "Subseries Plot: Beer Production by Quarter")
1.2 Differencing to make it stationarity
beer %>%
features(Beer, unitroot_nsdiffs) # seasonal differences ## # A tibble: 1 × 1
## nsdiffs
## <int>
## 1 1beer %>%
features(Beer, unitroot_ndiffs) # regular differences ## # A tibble: 1 × 1
## ndiffs
## <int>
## 1 1beer <- beer %>%
mutate(
diff_seasonal = difference(Beer, lag = 4),
diff_both = difference(diff_seasonal, lag = 1)
)
beer %>% autoplot(diff_seasonal) +
labs(title = "Seasonally Differenced Beer Production (lag = 4)")
beer %>% autoplot(diff_both) +
labs(title = "Seasonally + First Differenced Beer Production")
# ADF test on raw series
adf.test(na.omit(beer$Beer))##
## Augmented Dickey-Fuller Test
##
## data: na.omit(beer$Beer)
## Dickey-Fuller = -1.2827, Lag order = 6, p-value = 0.877
## alternative hypothesis: stationary# ADF test on differenced series
adf.test(na.omit(beer$diff_both))##
## Augmented Dickey-Fuller Test
##
## data: na.omit(beer$diff_both)
## Dickey-Fuller = -10.649, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryunitroot_nsdiffs recommends 1 seasonal difference and unitroot_ndiffs recommends another regular difference after that. The ADF test on the raw series gives DF = -1.28, p = 0.877 (non-stationary), while the ADF test on the seasonally + regularly differenced series gives DF = -10.65, p < 0.01 (stationary). This matches the ndiffs/nsdiffs recommendation exactly: one seasonal + one regular difference is enough to reach stationarity.
1.3 Fitting the ARIMA models
fit_arima <- beer %>%
model(
manual = ARIMA(Beer ~ pdq(2,1,1) + PDQ(0,1,1)), # plausible manual guess from ACF/PACF
auto = ARIMA(Beer) # let auto_arima search
)
fit_arima %>%
pivot_longer(everything(), names_to = "Model", values_to = "Orders") %>%
print()## # A mable: 2 x 2
## # Key: Model [2]
## Model Orders
## <chr> <model>
## 1 manual <ARIMA(2,1,1)(0,1,1)[4]>
## 2 auto <ARIMA(1,1,2)(0,1,1)[4]>glance(fit_arima) %>%
select(.model, AIC, AICc, BIC, log_lik)## # A tibble: 2 × 5
## .model AIC AICc BIC log_lik
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 manual 1784. 1784. 1800. -887.
## 2 auto 1783. 1783. 1800. -886.fit_arima %>% select(auto) %>% report()## Series: Beer
## Model: ARIMA(1,1,2)(0,1,1)[4]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.0495 -1.0091 0.3746 -0.7434
## s.e. 0.1959 0.1826 0.1530 0.0502
##
## sigma^2 estimated as 241.3: log likelihood=-886.41
## AIC=1782.82 AICc=1783.11 BIC=1799.63The auto-selected model is ARIMA(1,1,2)(0,1,1)[4] (AIC = 1782.8, AICc = 1783.1, BIC = 1799.6), which beats my manually-specified ARIMA(2,1,1)(0,1,1)[4] (AIC = 1784, AICc = 1784, BIC = 1800) on all three criteria, though the margin is fairly small. Interestingly the auto model’s AR term is tiny and not significant (ar1 = 0.05, se = 0.20), so most of the structure is really being captured by the two MA terms plus the seasonal MA term. Whereas my manual guess of AR(2) turned out to be more complex than needed. This is a good show of why auto_arima/ARIMA’s stepwise search is worth running even when you have a plausible manual guess from the ACF/PACF.
1.4 Residual diagnostics
fit_arima %>%
select(auto) %>%
gg_tsresiduals() +
labs(title = "Residual Diagnostics: Auto ARIMA")
#test for autocorrelation in residuals
fit_arima %>%
select(auto) %>%
augment() %>%
features(.innov, ljung_box, lag = 8, dof = 3)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 auto 4.06 0.540The Ljung-Box on the auto ARIMA residuals gives Q* = 4.06, df = 5, p = 0.540, well above 0.05, so we fail to reject “no autocorrelation.” The residuals look like white noise, meaning the model has captured the available structure in the series (the manual model also passes, with p = 0.315).
1.5 Forecasting and model comparison
# Last 8 quarters as test set
test_start <- max(beer$Quarter) - 7
train <- beer %>% filter(Quarter < test_start)
test <- beer %>% filter(Quarter >= test_start)
fit_compare <- train %>%
model(
auto_arima = ARIMA(Beer),
manual_arima = ARIMA(Beer ~ pdq(2,1,1) + PDQ(0,1,1)),
ets = ETS(Beer)
)
fc_compare <- fit_compare %>% forecast(new_data = test)
fc_compare %>%
autoplot(beer %>% filter(Quarter >= test_start - 8), level = NULL) +
labs(title = "Forecast Comparison: ARIMA (manual/auto) vs ETS")
accuracy(fc_compare, beer) %>%
select(.model, RMSE, MAE, MAPE)## # A tibble: 3 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 auto_arima 14.5 12.2 2.86
## 2 ets 12.8 11.3 2.65
## 3 manual_arima 14.1 12.3 2.88Despite the ARIMA models having the better in-sample AIC/BIC, ETS actually produced the most accurate out-of-sample forecasts on this holdout. That’s a useful finding because it shows AIC/BIC measure in-sample fit relative to model complexity, not genuine forecast accuracy, and the two don’t always agree. For a series like beer production which is dominated by a clean, stable trend plus seasonal pattern with modest residual autocorrelation, ETS’s smoothing-based structure can match or beat ARIMA’s explicit AR/MA terms. This doesn’t mean ARIMA “looses” in general; on series with strong short-term autocorrelation, ARIMA typically wins. The main takeaway is to always validate with an out-of-sample test rather than picking a model on AIC alone.
2 Part 2: Dynamic Regression Models
For this I decided to use vic_elec which is the half-hourly Victorian electricity demand, which includes temperature as a natural external predictor, aggregated to daily to keep the model manageable. The daily data spans from 2012-01-01 to 2014-12-31 with a total n of 1096.
2.1 Data prep and dynamic regression model
vic_daily <- vic_elec %>%
index_by(Date = as_date(Time)) %>%
summarise(
Demand = sum(Demand),
Temperature = mean(Temperature),
Holiday = any(Holiday)
) %>%
mutate(Day_Type = case_when(
Holiday ~ "Holiday",
wday(Date) %in% c(1,7) ~ "Weekend",
TRUE ~ "Weekday"
))
vic_daily %>%
autoplot(Demand) +
labs(title = "Daily Electricity Demand, Victoria")
vic_daily %>%
ggplot(aes(x = Temperature, y = Demand)) +
geom_point(alpha = 0.4) +
labs(title = "Demand vs Temperature (U-shape: heating + cooling load)")
The scatterplot shows a clear U-shape which means that demand is lowest around 15-20C and rises at both colder and hotter extremes, consistent with heating and cooling load.
fit_dynreg <- vic_daily %>%
model(
dynreg = ARIMA(Demand ~ Temperature + I(Temperature^2) + Day_Type)
)
report(fit_dynreg)## Series: Demand
## Model: LM w/ ARIMA(2,1,1)(2,0,0)[7] errors
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 Temperature I(Temperature^2)
## 0.6208 0.1056 -0.9673 0.2658 0.2069 -11726.2619 362.9539
## s.e. 0.0364 0.0334 0.0120 0.0308 0.0310 341.9217 8.7085
## Day_TypeWeekday Day_TypeWeekend
## 31174.180 -1611.768
## s.e. 1002.141 1188.094
##
## sigma^2 estimated as 40636913: log likelihood=-11142.38
## AIC=22304.76 AICc=22304.97 BIC=22354.75The positive Temperature^2 coefficient confirms the U-shape: demand falls as temperature rises from cold toward mild (heating demand drops), then rises again as it gets hot (cooling/air-conditioning demand kicks in). It is exactly what the scatter plot of Demand vs Temperature shows. The Day_Type coefficients are relative to the baseline (Holiday): weekdays show substantially and significantly higher demand than holidays (about 31,174 units higher, se = 1,002), consistent with more commercial/industrial activity on work days, while weekends are only modestly and not clearly significantly lower than holidays (about -1,612 units, se = 1,188, well under 2 standard errors from zero).
2.2 Diagnostics for the dynamic regression model
fit_dynreg %>% gg_tsresiduals() +
labs(title = "Residual Diagnostics: Dynamic Regression")
fit_dynreg %>%
augment() %>%
features(.innov, ljung_box, lag = 14, dof = 6)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 dynreg 34.9 0.0000284Even with ARIMA(2,1,1)(2,0,0)[7] errors already fitted, Ljung-Box (Q* = 34.9, df = 8, p < 0.001) still shows significant leftover autocorrelation which is limiting.Daily-aggregated demand has a strong weekly cycle for example: Monday consistently differs from Friday, that a 3-level Weekday/Weekend/Holiday factor doesn’t fully capture. A next step worth trying and worth would be adding day-of-week dummies or Fourier terms for a period-7 seasonal cycle, or working directly with the half-hourly data and multiple seasonal periods (daily +weekly) via fable’s multiple-seasonality support.
2.3 ARIMAX: comparing to a pure regression (no ARIMA errors)
train_vic <- vic_daily %>% filter(Date < max(Date) - 29)
test_vic <- vic_daily %>% filter(Date >= max(Date) - 29)
fit_vic_compare <- train_vic %>%
model(
plain_reg = TSLM(Demand ~ Temperature + I(Temperature^2) + Day_Type),
arimax = ARIMA(Demand ~ Temperature + I(Temperature^2) + Day_Type)
)
glance(fit_vic_compare) %>% select(.model, AIC, AICc, BIC)## # A tibble: 2 × 4
## .model AIC AICc BIC
## <chr> <dbl> <dbl> <dbl>
## 1 plain_reg 19740. 19740. 19770.
## 2 arimax 21684. 21685. 21734.# Need future predictor values to forecast forward
fc_vic <- fit_vic_compare %>%
forecast(new_data = test_vic)
# Zoomed to 30 days of history + the 30-day test window, same reasoning
# as the beer plot above.
fc_vic %>%
autoplot(vic_daily %>% filter(Date >= max(Date) - 29 - 30), level = NULL) +
labs(title = "ARIMAX vs Plain Regression Forecasts")
accuracy(fc_vic, vic_daily) %>%
select(.model, RMSE, MAE, MAPE)## # A tibble: 2 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 arimax 13560. 9325. 4.83
## 2 plain_reg 14508. 10086. 5.20The plain regression’s residuals show heavy leftover autocorrelation which is completely expected, since forcing ARIMA(0,0,0) errors is equivalent to assuming i.i.d. residuals, which clearly isn’t true for a daily demand series. ARIMAX, which lets auto.arima search for the right ARIMA error structure, actually has a much higher in-sample AIC than the plain regression (21684 vs. 19740) — though that’s not really an apples-to-apples comparison, since the two models differ in their order of differencing and so aren’t on the same likelihood scale. Where ARIMAX does clearly win is out-of-sample: it forecasts more accurately, with RMSE improving by about 6.5% and MAE by about 7.5%. This is a clean illustration of the core idea behind dynamic regression / ARIMAX: once you’ve included your external predictors, don’t just assume the leftover errors are white noise. Model them explicitly if they’re not, and you’ll typically get better forecasts, even if in-sample AIC comparisons between differently-specified models aren’t reliable.
3 Reflection
Based on my results, the two datasets tell somewhat different stories, which is itself the interesting finding.
For beer production: a series with no obvious external driver, plain ETS beat both ARIMA specifications out-of-sample (RMSE 12.8 vs. 14.1-14.5), even though the auto-selected ARIMA(1,1,2)(0,1,1)[4] had the better in-sample AIC. This is a reminder that AIC/BIC reward in-sample fit, not necessarily forecast accuracy, and that a simpler smoothing-based model can be perfectly competitive when a series is dominated by clean trend/seasonal structure.
For electricity demand, where a real causal driver (temperature, via heating/cooling load) and a real calendar effect (weekday vs. weekend vs. holiday) are available, ARIMAX still won out-of-sample, even though its in-sample AIC was actually higher than the plain regression’s (21684 vs. 19740, not a fair comparison given the different differencing orders): it cut test RMSE by about 6.5% and MAE by about 7.5% relative to a plain regression, by explicitly modeling the autocorrelation the plain regression’s residuals still contained. Even so, ARIMAX didn’t fully eliminate residual autocorrelation a leftover weekly cycle remained showing that external predictors help most when they capture a real, causal driver of the series, but they don’t excuse you from also getting the time-series (ARIMA) error structure right.
My overall takeaway: reach for dynamic regression/ARIMAX when you have predictors with real explanatory power and known/forecastable future values; when you don’t, a well-chosen univariate model (ARIMA or ETS,chosen by out-of-sample validation rather than AIC alone) may already be hard to beat.