Australian Beer Production ARIMA Analysis
Meghma Banerjee
2026-07-30
options(repos = c(CRAN = "https://cloud.r-project.org"))
install.packages(
c("fpp3", "dplyr", "ggplot2", "tsibble", "tsibbledata",
"feasts", "fable", "fabletools", "tseries"),
dependencies = TRUE
)Step 1: Prepare the univariate time series
beer_data <- tsibbledata::aus_production |>
dplyr::select(Quarter, Beer) |>
dplyr::filter(!is.na(Beer))
print(beer_data)## # A tsibble: 218 x 2 [1Q]
## Quarter Beer
## <qtr> <dbl>
## 1 1956 Q1 284
## 2 1956 Q2 213
## 3 1956 Q3 227
## 4 1956 Q4 308
## 5 1957 Q1 262
## 6 1957 Q2 228
## 7 1957 Q3 236
## 8 1957 Q4 320
## 9 1958 Q1 272
## 10 1958 Q2 233
## # ℹ 208 more rowshead(beer_data)## # A tsibble: 6 x 2 [1Q]
## Quarter Beer
## <qtr> <dbl>
## 1 1956 Q1 284
## 2 1956 Q2 213
## 3 1956 Q3 227
## 4 1956 Q4 308
## 5 1957 Q1 262
## 6 1957 Q2 228summary(beer_data$Beer)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 213.0 379.2 422.0 415.4 465.8 599.0Plot the original series
# Check for trend, seasonality, and changes in variation
# ------------------------------------------------------------
beer_data |>
autoplot(Beer) +
labs(
title = "Quarterly Australian Beer Production",
subtitle = "Original Time Series",
x = "Quarter",
y = "Beer Production (Megalitres)"
) +
theme_minimal()
# ------------------------------------------------------------Check seasonality
feasts::gg_season(
beer_data,
Beer,
labels = "both"
) +
ggplot2::labs(
title = "Seasonal Plot of Australian Beer Production",
x = "Quarter",
y = "Beer Production"
) +
ggplot2::theme_minimal()
feasts::gg_subseries(
beer_data,
Beer
) +
ggplot2::labs(
title = "Seasonal Subseries Plot of Beer Production",
x = "Quarter",
y = "Beer Production"
) +
ggplot2::theme_minimal()
# ------------------------------------------------------------Plot the ACF of the original series
beer_data |>
ACF(Beer) |>
autoplot() +
labs(
title = "ACF of Original Beer Production Series"
) +
theme_minimal()
# ------------------------------------------------------------Determine recommended differencing
# Number of regular differences suggested
regular_differences <- beer_data |>
features(Beer, unitroot_ndiffs)
print(regular_differences)## # A tibble: 1 × 1
## ndiffs
## <int>
## 1 1# Number of seasonal differences suggested
seasonal_differences <- beer_data |>
features(Beer, unitroot_nsdiffs)
print(seasonal_differences)## # A tibble: 1 × 1
## nsdiffs
## <int>
## 1 1# ------------------------------------------------------------ADF test on the original series
# Null hypothesis:
# The series has a unit root and is nonstationary.
#
# If p-value < 0.05:
# Reject the null hypothesis; evidence supports stationarity.
#
# If p-value > 0.05:
# Fail to reject the null hypothesis; evidence of
# nonstationarity remains.
adf_original <- tseries::adf.test(
as.numeric(beer_data$Beer)
)
print(adf_original)##
## Augmented Dickey-Fuller Test
##
## data: as.numeric(beer_data$Beer)
## Dickey-Fuller = -1.2827, Lag order = 6, p-value = 0.877
## alternative hypothesis: stationaryApply seasonal differencing
# Quarterly data have a seasonal period of 4
# ============================================================
beer_diff <- beer_data
beer_diff$seasonal_diff <- c(
rep(NA_real_, 4),
diff(
as.numeric(beer_data$Beer),
lag = 4
)
)
# View the new column
head(beer_diff, 10)## # A tsibble: 10 x 3 [1Q]
## Quarter Beer seasonal_diff
## <qtr> <dbl> <dbl>
## 1 1956 Q1 284 NA
## 2 1956 Q2 213 NA
## 3 1956 Q3 227 NA
## 4 1956 Q4 308 NA
## 5 1957 Q1 262 -22
## 6 1957 Q2 228 15
## 7 1957 Q3 236 9
## 8 1957 Q4 320 12
## 9 1958 Q1 272 10
## 10 1958 Q2 233 5# Plot the seasonally differenced series
seasonal_difference_plot <- ggplot2::ggplot(
beer_diff,
ggplot2::aes(
x = Quarter,
y = seasonal_diff
)
) +
ggplot2::geom_line() +
ggplot2::labs(
title = "Seasonally Differenced Beer Production",
subtitle = "Beer(t) minus Beer(t − 4)",
x = "Quarter",
y = "Seasonal Difference"
) +
ggplot2::theme_minimal()
print(seasonal_difference_plot)
# ============================================================ADF test after seasonal differencing
seasonal_diff_values <- beer_diff$seasonal_diff
# Remove missing values
seasonal_diff_values <- seasonal_diff_values[
!is.na(seasonal_diff_values)
]
adf_seasonal <- tseries::adf.test(
seasonal_diff_values
)
print(adf_seasonal)##
## Augmented Dickey-Fuller Test
##
## data: seasonal_diff_values
## Dickey-Fuller = -4.6644, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary# ============================================================Apply regular and seasonal differencing
# First four values are already missing because of the
# seasonal difference. Applying diff() introduces one more
# missing value.
regular_seasonal_values <- diff(
seasonal_diff_values,
lag = 1
)
beer_diff$regular_seasonal_diff <- c(
rep(NA_real_, 5),
regular_seasonal_values
)
# Confirm the values were added
head(beer_diff, 10)## # A tsibble: 10 x 4 [1Q]
## Quarter Beer seasonal_diff regular_seasonal_diff
## <qtr> <dbl> <dbl> <dbl>
## 1 1956 Q1 284 NA NA
## 2 1956 Q2 213 NA NA
## 3 1956 Q3 227 NA NA
## 4 1956 Q4 308 NA NA
## 5 1957 Q1 262 -22 NA
## 6 1957 Q2 228 15 37
## 7 1957 Q3 236 9 -6
## 8 1957 Q4 320 12 3
## 9 1958 Q1 272 10 -2
## 10 1958 Q2 233 5 -5# Plot the regular and seasonally differenced series
final_difference_plot <- ggplot2::ggplot(
beer_diff,
ggplot2::aes(
x = Quarter,
y = regular_seasonal_diff
)
) +
ggplot2::geom_line() +
ggplot2::labs(
title = "Regular and Seasonally Differenced Beer Production",
subtitle = "One regular difference and one seasonal difference",
x = "Quarter",
y = "Differenced Beer Production"
) +
ggplot2::theme_minimal()
print(final_difference_plot)
# ============================================================ADF test after both differences
final_diff_values <- beer_diff$regular_seasonal_diff
# Remove missing values
final_diff_values <- final_diff_values[
!is.na(final_diff_values)
]
adf_final <- tseries::adf.test(
final_diff_values
)
print(adf_final)##
## Augmented Dickey-Fuller Test
##
## data: final_diff_values
## Dickey-Fuller = -10.649, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary# ============================================================Plot the ACF and PACF of the differenced series
# Create a separate tsibble containing the final
# differenced series
final_diff_data <- beer_diff[
!is.na(beer_diff$regular_seasonal_diff),
c("Quarter", "regular_seasonal_diff")
]
final_diff_data <- tsibble::as_tsibble(
final_diff_data,
index = Quarter
)
# ACF
final_acf <- feasts::ACF(
final_diff_data,
regular_seasonal_diff
)
final_acf_plot <- ggplot2::autoplot(
final_acf
) +
ggplot2::labs(
title = "ACF of Differenced Beer Production",
x = "Lag",
y = "Autocorrelation"
) +
ggplot2::theme_minimal()
print(final_acf_plot)
# PACF
final_pacf <- feasts::PACF(
final_diff_data,
regular_seasonal_diff
)
final_pacf_plot <- ggplot2::autoplot(
final_pacf
) +
ggplot2::labs(
title = "PACF of Differenced Beer Production",
x = "Lag",
y = "Partial Autocorrelation"
) +
ggplot2::theme_minimal()
print(final_pacf_plot)
# ============================================================Fit manual and automatic ARIMA models
beer_models <- fabletools::model(
beer_data,
# Manual nonseasonal ARIMA(1,1,1)
manual_arima = fable::ARIMA(
Beer ~
pdq(1, 1, 1) +
PDQ(0, 0, 0)
),
# Manual seasonal ARIMA(1,1,1)(0,1,1)[4]
manual_seasonal_arima = fable::ARIMA(
Beer ~
pdq(1, 1, 1) +
PDQ(0, 1, 1)
),
# Automatically selected ARIMA model
automatic_arima = fable::ARIMA(
Beer
)
)
print(beer_models)## # A mable: 1 x 3
## manual_arima manual_seasonal_arima automatic_arima
## <model> <model> <model>
## 1 <ARIMA(1,1,1) w/ drift> <ARIMA(1,1,1)(0,1,1)[4]> <ARIMA(1,1,2)(0,1,1)[4]># ============================================================View model estimates
# Report manual ARIMA
manual_report <- beer_models[
"manual_arima"
]
fabletools::report(manual_report)## Series: Beer
## Model: ARIMA(1,1,1) w/ drift
##
## Coefficients:
## ar1 ma1 constant
## -0.0395 -0.8538 0.8209
## s.e. 0.0721 0.0251 0.5150
##
## sigma^2 estimated as 2581: log likelihood=-1159.44
## AIC=2326.87 AICc=2327.06 BIC=2340.39# Report manual seasonal ARIMA
seasonal_report <- beer_models[
"manual_seasonal_arima"
]
fabletools::report(seasonal_report)## Series: Beer
## Model: ARIMA(1,1,1)(0,1,1)[4]
##
## Coefficients:
## ar1 ma1 sma1
## -0.3442 -0.5874 -0.7118
## s.e. 0.0801 0.0651 0.0511
##
## sigma^2 estimated as 244.6: log likelihood=-888.37
## AIC=1784.73 AICc=1784.93 BIC=1798.18# Report automatically selected ARIMA
automatic_report <- beer_models[
"automatic_arima"
]
fabletools::report(automatic_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.63# View coefficients in table format
model_coefficients <- fabletools::tidy(
beer_models
)
print(model_coefficients)## # A tibble: 10 × 6
## .model term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 manual_arima ar1 -0.0395 0.0721 -0.548 5.84e- 1
## 2 manual_arima ma1 -0.854 0.0251 -34.1 4.81e-89
## 3 manual_arima constant 0.821 0.515 1.59 1.12e- 1
## 4 manual_seasonal_arima ar1 -0.344 0.0801 -4.30 2.65e- 5
## 5 manual_seasonal_arima ma1 -0.587 0.0651 -9.02 1.12e-16
## 6 manual_seasonal_arima sma1 -0.712 0.0511 -13.9 9.08e-32
## 7 automatic_arima ar1 0.0495 0.196 0.253 8.01e- 1
## 8 automatic_arima ma1 -1.01 0.183 -5.53 9.49e- 8
## 9 automatic_arima ma2 0.375 0.153 2.45 1.52e- 2
## 10 automatic_arima sma1 -0.743 0.0502 -14.8 1.29e-34# ============================================================Compare information criteria
model_comparison <- fabletools::glance(
beer_models
)
# Keep the important comparison columns
model_comparison <- model_comparison[
,
c(
".model",
"sigma2",
"log_lik",
"AIC",
"AICc",
"BIC"
)
]
# Sort from lowest to highest AICc
model_comparison <- model_comparison[
order(model_comparison$AICc),
]
print(model_comparison)## # A tibble: 3 × 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 automatic_arima 241. -886. 1783. 1783. 1800.
## 2 manual_seasonal_arima 245. -888. 1785. 1785. 1798.
## 3 manual_arima 2581. -1159. 2327. 2327. 2340.# The model at the top has the lowest AICc.
# ============================================================Compare in-sample accuracy
accuracy_results <- fabletools::accuracy(
beer_models
)
print(accuracy_results)## # A tibble: 3 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 manual_arima Trai… -0.0950 50.3 42.0 -1.06 9.98 2.71 2.60 -0.0255
## 2 manual_seasonal… Trai… -0.0747 15.3 11.6 0.0123 2.76 0.750 0.794 -0.0199
## 3 automatic_arima Trai… -0.0853 15.2 11.7 -0.00516 2.77 0.752 0.787 0.00560# Compare measures such as:
# RMSE = Root Mean Squared Error
# MAE = Mean Absolute Error
# MAPE = Mean Absolute Percentage Error
#
# Lower values indicate better in-sample accuracy.
# ============================================================Check residual diagnostics
# Manual ARIMA residual diagnostics
manual_residual_plot <- feasts::gg_tsresiduals(
beer_models["manual_arima"]
) +
ggplot2::labs(
title = "Residual Diagnostics: Manual ARIMA(1,1,1)"
)
print(manual_residual_plot)
# Manual seasonal ARIMA residual diagnostics
seasonal_residual_plot <- feasts::gg_tsresiduals(
beer_models["manual_seasonal_arima"]
) +
ggplot2::labs(
title = paste(
"Residual Diagnostics:",
"Manual Seasonal ARIMA"
)
)
print(seasonal_residual_plot)
# Automatic ARIMA residual diagnostics
automatic_residual_plot <- feasts::gg_tsresiduals(
beer_models["automatic_arima"]
) +
ggplot2::labs(
title = "Residual Diagnostics: Automatic ARIMA"
)
print(automatic_residual_plot)
# Good residuals should:
# 1. Fluctuate around zero
# 2. Have approximately constant variance
# 3. Show no strong pattern
# 4. Have no major residual ACF spikes
# 5. Be approximately normally distributed
# ============================================================Run the Ljung–Box white-noise test
# Obtain model innovations/residuals
beer_residuals <- fabletools::augment(
beer_models
)
# Run the Ljung-Box test for each model
residual_tests <- fabletools::features(
beer_residuals,
.innov,
feasts::ljung_box,
lag = 8,
dof = 0
)
print(residual_tests)## # A tibble: 3 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 automatic_arima 4.06 0.851
## 2 manual_arima 543. 0
## 3 manual_seasonal_arima 8.07 0.427# Interpretation:
#
# H0: Residuals are independently distributed white noise.
#
# If p-value > 0.05:
# Fail to reject H0. Residuals are consistent with white noise.
#
# If p-value < 0.05:
# Reject H0. Significant autocorrelation remains, suggesting
# that the model may be inadequate.
# ============================================================Forecast using all models
beer_forecasts <- fabletools::forecast(
beer_models,
h = "3 years"
)
print(beer_forecasts)## # A fable: 36 x 4 [1Q]
## # Key: .model [3]
## .model Quarter
## <chr> <qtr>
## 1 manual_arima 2010 Q3
## 2 manual_arima 2010 Q4
## 3 manual_arima 2011 Q1
## 4 manual_arima 2011 Q2
## 5 manual_arima 2011 Q3
## 6 manual_arima 2011 Q4
## 7 manual_arima 2012 Q1
## 8 manual_arima 2012 Q2
## 9 manual_arima 2012 Q3
## 10 manual_arima 2012 Q4
## # ℹ 26 more rows
## # ℹ 2 more variables: Beer <dist>, .mean <dbl>all_forecasts_plot <- autoplot(
beer_forecasts,
beer_data
) +
ggplot2::labs(
title = "ARIMA Forecasts of Australian Beer Production",
subtitle = "Comparison of Manual and Automatic Models",
x = "Quarter",
y = "Beer Production (Megalitres)"
) +
ggplot2::theme_minimal()
print(all_forecasts_plot)
# ============================================================Forecast using the best model
# Start by looking at the AICc comparison again
print(model_comparison)## # A tibble: 3 × 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 automatic_arima 241. -886. 1783. 1783. 1800.
## 2 manual_seasonal_arima 245. -888. 1785. 1785. 1798.
## 3 manual_arima 2581. -1159. 2327. 2327. 2340.# This selects the model with the lowest AICc automatically
best_model_name <- as.character(
model_comparison$.model[1]
)
print(best_model_name)## [1] "automatic_arima"# Select the best model
best_model <- beer_models[
best_model_name
]
print(best_model)## # A mable: 1 x 1
## automatic_arima
## <model>
## 1 <ARIMA(1,1,2)(0,1,1)[4]># Forecast three years ahead
best_forecast <- fabletools::forecast(
best_model,
h = "3 years"
)
print(best_forecast)## # A fable: 12 x 4 [1Q]
## # Key: .model [1]
## .model Quarter
## <chr> <qtr>
## 1 automatic_arima 2010 Q3
## 2 automatic_arima 2010 Q4
## 3 automatic_arima 2011 Q1
## 4 automatic_arima 2011 Q2
## 5 automatic_arima 2011 Q3
## 6 automatic_arima 2011 Q4
## 7 automatic_arima 2012 Q1
## 8 automatic_arima 2012 Q2
## 9 automatic_arima 2012 Q3
## 10 automatic_arima 2012 Q4
## 11 automatic_arima 2013 Q1
## 12 automatic_arima 2013 Q2
## # ℹ 2 more variables: Beer <dist>, .mean <dbl># Plot the best forecast
best_forecast_plot <- autoplot(
best_forecast,
beer_data
) +
ggplot2::labs(
title = "Best ARIMA Forecast of Australian Beer Production",
subtitle = paste(
"Lowest-AICc model:",
best_model_name
),
x = "Quarter",
y = "Beer Production (Megalitres)"
) +
ggplot2::theme_minimal()
print(best_forecast_plot)
Part 2: Forecast Accuracy and Comparison with Exponential Smoothing
Create training and test sets
# Reserve the final 12 quarters (3 years) for forecast evaluation
n_test <- 12
# Identify the final quarter included in the training data
training_cutoff <- beer_data$Quarter[nrow(beer_data) - n_test]
beer_train <- beer_data |>
dplyr::filter(Quarter <= training_cutoff)
beer_test <- beer_data |>
dplyr::filter(Quarter > training_cutoff)
cat("Training observations:", nrow(beer_train), "\n")## Training observations: 206cat("Test observations:", nrow(beer_test), "\n")## Test observations: 12cat("Training ends at:", as.character(max(beer_train$Quarter)), "\n")## Training ends at: 2007 Q2cat("Test begins at:", as.character(min(beer_test$Quarter)), "\n")## Test begins at: 2007 Q3# Plot the training and test periods
beer_data |>
ggplot2::ggplot(ggplot2::aes(x = Quarter, y = Beer)) +
ggplot2::geom_line() +
ggplot2::geom_vline(
xintercept = as.numeric(training_cutoff),
linetype = "dashed"
) +
ggplot2::labs(
title = "Training and Test Periods",
subtitle = "The final 12 quarters are reserved for forecast evaluation",
x = "Quarter",
y = "Beer Production (Megalitres)"
) +
ggplot2::theme_minimal()
Fit ARIMA and exponential smoothing models
comparison_models <- beer_train |>
fabletools::model(
manual_arima = fable::ARIMA(
Beer ~ fable::pdq(1, 1, 1) + fable::PDQ(0, 1, 1)
),
automatic_arima = fable::ARIMA(Beer),
exponential_smoothing = fable::ETS(Beer)
)
print(comparison_models)## # A mable: 1 x 3
## manual_arima automatic_arima exponential_smoothing
## <model> <model> <model>
## 1 <NULL model> <ARIMA(1,1,2)(0,1,1)[4]> <ETS(M,A,M)># Display the selected model structures and coefficients
fabletools::report(comparison_models["manual_arima"])## Series: Beer
## Model: NULL model
## NULL modelfabletools::report(comparison_models["automatic_arima"])## Series: Beer
## Model: ARIMA(1,1,2)(0,1,1)[4]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.0394 -1.0187 0.3913 -0.7436
## s.e. 0.1954 0.1811 0.1486 0.0518
##
## sigma^2 estimated as 245.3: log likelihood=-838.13
## AIC=1686.26 AICc=1686.57 BIC=1702.78fabletools::report(comparison_models["exponential_smoothing"])## Series: Beer
## Model: ETS(M,A,M)
## Smoothing parameters:
## alpha = 0.2070047
## beta = 0.03209157
## gamma = 0.1994071
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 256.1253 0.6053282 1.183274 0.9110569 0.8615274 1.044142
##
## sigma^2: 0.0013
##
## AIC AICc BIC
## 2217.957 2218.875 2247.908Generate forecasts for the test period
comparison_forecasts <- comparison_models |>
fabletools::forecast(h = n_test)
print(comparison_forecasts)## # A fable: 36 x 4 [1Q]
## # Key: .model [3]
## .model Quarter Beer .mean
## <chr> <qtr> <dist> <dbl>
## 1 manual_arima 2007 Q3 NA NA
## 2 manual_arima 2007 Q4 NA NA
## 3 manual_arima 2008 Q1 NA NA
## 4 manual_arima 2008 Q2 NA NA
## 5 manual_arima 2008 Q3 NA NA
## 6 manual_arima 2008 Q4 NA NA
## 7 manual_arima 2009 Q1 NA NA
## 8 manual_arima 2009 Q2 NA NA
## 9 manual_arima 2009 Q3 NA NA
## 10 manual_arima 2009 Q4 NA NA
## # ℹ 26 more rows# Plot forecasts against the actual observations in the test set
comparison_forecasts |>
autoplot(beer_train) +
ggplot2::autolayer(
beer_test,
Beer,
series = "Actual test values"
) +
ggplot2::labs(
title = "ARIMA and Exponential Smoothing Forecasts",
subtitle = "Forecasts compared with the final 12 observed quarters",
x = "Quarter",
y = "Beer Production (Megalitres)"
) +
ggplot2::theme_minimal()
Evaluate forecast accuracy using RMSE and MAE
forecast_accuracy <- comparison_forecasts |>
fabletools::accuracy(beer_test) |>
dplyr::select(.model, RMSE, MAE, MAPE) |>
dplyr::arrange(RMSE)
print(forecast_accuracy)## # A tibble: 3 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 exponential_smoothing 9.62 8.92 2.13
## 2 automatic_arima 9.83 9.04 2.16
## 3 manual_arima NaN NaN NaN# The preferred forecasting model has the lowest test-set RMSE and MAE.Compare AIC and AICc
information_criteria <- comparison_models |>
generics::glance() |>
dplyr::select(.model, AIC, AICc, BIC) |>
dplyr::arrange(AICc)
print(information_criteria)## # A tibble: 2 × 4
## .model AIC AICc BIC
## <chr> <dbl> <dbl> <dbl>
## 1 automatic_arima 1686. 1687. 1703.
## 2 exponential_smoothing 2218. 2219. 2248.# AIC and AICc evaluate in-sample model fit while penalizing complexity.
# Lower values are preferred, but test-set RMSE and MAE provide the more
# direct measure of forecast performance on unseen observations.Combine the accuracy and information-criterion results
model_evaluation <- forecast_accuracy |>
dplyr::left_join(
information_criteria,
by = ".model"
) |>
dplyr::arrange(RMSE)
print(model_evaluation)## # A tibble: 3 × 7
## .model RMSE MAE MAPE AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 exponential_smoothing 9.62 8.92 2.13 2218. 2219. 2248.
## 2 automatic_arima 9.83 9.04 2.16 1686. 1687. 1703.
## 3 manual_arima NaN NaN NaN NA NA NAPart 2: Dynamic Regression Models (Chapter 10)
Prepare the Victorian electricity demand dataset
vic_data <- tsibbledata::vic_elec |>
dplyr::select(Time, Demand, Temperature, Holiday) |>
tidyr::drop_na()
print(vic_data)## # A tsibble: 52,608 x 4 [30m] <Australia/Melbourne>
## Time Demand Temperature Holiday
## <dttm> <dbl> <dbl> <lgl>
## 1 2012-01-01 00:00:00 4383. 21.4 TRUE
## 2 2012-01-01 00:30:00 4263. 21.0 TRUE
## 3 2012-01-01 01:00:00 4049. 20.7 TRUE
## 4 2012-01-01 01:30:00 3878. 20.6 TRUE
## 5 2012-01-01 02:00:00 4036. 20.4 TRUE
## 6 2012-01-01 02:30:00 3866. 20.2 TRUE
## 7 2012-01-01 03:00:00 3694. 20.1 TRUE
## 8 2012-01-01 03:30:00 3562. 19.6 TRUE
## 9 2012-01-01 04:00:00 3433. 19.1 TRUE
## 10 2012-01-01 04:30:00 3359. 19.0 TRUE
## # ℹ 52,598 more rowssummary(vic_data)## Time Demand Temperature Holiday
## Min. :2012-01-01 00:00:00 Min. :2858 Min. : 1.50 Mode :logical
## 1st Qu.:2012-09-30 22:52:30 1st Qu.:3969 1st Qu.:12.30 FALSE:51120
## Median :2013-07-01 22:45:00 Median :4635 Median :15.40 TRUE :1488
## Mean :2013-07-01 22:45:00 Mean :4665 Mean :16.27
## 3rd Qu.:2014-04-01 23:37:30 3rd Qu.:5244 3rd Qu.:19.40
## Max. :2014-12-31 23:30:00 Max. :9345 Max. :43.20Visualize electricity demand and external predictors
vic_data |>
feasts::autoplot(Demand) +
ggplot2::labs(
title = "Victorian Half-Hourly Electricity Demand",
x = "Time",
y = "Electricity Demand"
) +
ggplot2::theme_minimal()
vic_data |>
ggplot2::ggplot(
ggplot2::aes(x = Temperature, y = Demand)
) +
ggplot2::geom_point(alpha = 0.15) +
ggplot2::geom_smooth(method = "lm", se = FALSE) +
ggplot2::labs(
title = "Electricity Demand and Temperature",
x = "Temperature",
y = "Electricity Demand"
) +
ggplot2::theme_minimal()
vic_data |>
ggplot2::ggplot(
ggplot2::aes(x = Holiday, y = Demand)
) +
ggplot2::geom_boxplot() +
ggplot2::labs(
title = "Electricity Demand by Holiday Status",
x = "Holiday",
y = "Electricity Demand"
) +
ggplot2::theme_minimal()
Create training and test sets
# Use the final seven days as the test period.
# The data are measured every 30 minutes, so seven days contain 336 observations.
n_test_vic <- 7 * 48
n_total_vic <- nrow(vic_data)
vic_train <- vic_data |>
dplyr::slice_head(n = n_total_vic - n_test_vic)
vic_test <- vic_data |>
dplyr::slice_tail(n = n_test_vic)
print(vic_train)## # A tsibble: 52,272 x 4 [30m] <Australia/Melbourne>
## Time Demand Temperature Holiday
## <dttm> <dbl> <dbl> <lgl>
## 1 2012-01-01 00:00:00 4383. 21.4 TRUE
## 2 2012-01-01 00:30:00 4263. 21.0 TRUE
## 3 2012-01-01 01:00:00 4049. 20.7 TRUE
## 4 2012-01-01 01:30:00 3878. 20.6 TRUE
## 5 2012-01-01 02:00:00 4036. 20.4 TRUE
## 6 2012-01-01 02:30:00 3866. 20.2 TRUE
## 7 2012-01-01 03:00:00 3694. 20.1 TRUE
## 8 2012-01-01 03:30:00 3562. 19.6 TRUE
## 9 2012-01-01 04:00:00 3433. 19.1 TRUE
## 10 2012-01-01 04:30:00 3359. 19.0 TRUE
## # ℹ 52,262 more rowsprint(vic_test)## # A tsibble: 336 x 4 [30m] <Australia/Melbourne>
## Time Demand Temperature Holiday
## <dttm> <dbl> <dbl> <lgl>
## 1 2014-12-25 00:00:00 4042. 16.4 TRUE
## 2 2014-12-25 00:30:00 4053. 16.4 TRUE
## 3 2014-12-25 01:00:00 3821. 16.2 TRUE
## 4 2014-12-25 01:30:00 3624. 16.1 TRUE
## 5 2014-12-25 02:00:00 3470. 15.8 TRUE
## 6 2014-12-25 02:30:00 3318. 15.1 TRUE
## 7 2014-12-25 03:00:00 3202. 15 TRUE
## 8 2014-12-25 03:30:00 3113. 15.2 TRUE
## 9 2014-12-25 04:00:00 3056. 15.1 TRUE
## 10 2014-12-25 04:30:00 3032. 15.4 TRUE
## # ℹ 326 more rowsFit the dynamic regression model
dynamic_regression_model <- vic_train |>
fabletools::model(
dynamic_regression = fable::TSLM(
Demand ~ Temperature + Holiday + trend() + season()
)
)
print(dynamic_regression_model)## # A mable: 1 x 1
## dynamic_regression
## <model>
## 1 <TSLM>fabletools::report(dynamic_regression_model)## Series: Demand
## Model: TSLM
##
## Residuals:
## Min 1Q Median 3Q Max
## -1859.93 -667.19 -12.86 537.85 3772.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.144e+03 1.349e+01 307.257 <2e-16 ***
## Temperature 4.131e+01 6.451e-01 64.039 <2e-16 ***
## HolidayTRUE -7.069e+02 2.269e+01 -31.157 <2e-16 ***
## trend() -4.724e-03 2.417e-04 -19.544 <2e-16 ***
## season() -3.030e+00 7.282e+00 -0.416 0.677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 832.4 on 52267 degrees of freedom
## Multiple R-squared: 0.09162, Adjusted R-squared: 0.09155
## F-statistic: 1318 on 4 and 52267 DF, p-value: < 2.22e-16generics::tidy(dynamic_regression_model)## # A tibble: 5 × 6
## .model term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 dynamic_regression (Intercept) 4144. 13.5 307. 0
## 2 dynamic_regression Temperature 41.3 0.645 64.0 0
## 3 dynamic_regression HolidayTRUE -707. 22.7 -31.2 3.49e-211
## 4 dynamic_regression trend() -0.00472 0.000242 -19.5 9.39e- 85
## 5 dynamic_regression season() -3.03 7.28 -0.416 6.77e- 1Examine dynamic regression fitted values
dynamic_regression_model |>
fabletools::augment() |>
ggplot2::ggplot(
ggplot2::aes(x = Time)
) +
ggplot2::geom_line(
ggplot2::aes(y = Demand, linetype = "Actual")
) +
ggplot2::geom_line(
ggplot2::aes(y = .fitted, linetype = "Fitted")
) +
ggplot2::labs(
title = "Dynamic Regression: Actual and Fitted Demand",
x = "Time",
y = "Electricity Demand",
linetype = NULL
) +
ggplot2::theme_minimal()
Diagnose dynamic regression residuals
feasts::gg_tsresiduals(dynamic_regression_model)
dynamic_regression_residuals <- dynamic_regression_model |>
fabletools::augment()
dynamic_regression_ljung_box <- dynamic_regression_residuals |>
fabletools::features(
.innov,
feasts::ljung_box,
lag = 48,
dof = 4
)
print(dynamic_regression_ljung_box)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 dynamic_regression 529779. 0Fit the ARIMAX model
arimax_model <- vic_train |>
fabletools::model(
arimax = fable::ARIMA(
Demand ~ Temperature + Holiday
)
)
print(arimax_model)## # A mable: 1 x 1
## arimax
## <model>
## 1 <LM w/ ARIMA(2,1,1)(0,0,2)[2] errors>fabletools::report(arimax_model)## Series: Demand
## Model: LM w/ ARIMA(2,1,1)(0,0,2)[2] errors
##
## Coefficients:
## ar1 ar2 ma1 sma1 sma2 Temperature HolidayTRUE
## 1.8105 -0.8404 -0.9876 -0.2336 0.0358 2.3806 -17.1850
## s.e. 0.0028 0.0027 0.0018 0.0051 0.0049 0.6222 9.7583
##
## sigma^2 estimated as 8690: log likelihood=-311214.2
## AIC=622444.5 AICc=622444.5 BIC=622515.4generics::tidy(arimax_model)## # A tibble: 7 × 6
## .model term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 arimax ar1 1.81 0.00282 641. 0
## 2 arimax ar2 -0.840 0.00273 -308. 0
## 3 arimax ma1 -0.988 0.00178 -555. 0
## 4 arimax sma1 -0.234 0.00507 -46.1 0
## 5 arimax sma2 0.0358 0.00486 7.36 1.88e-13
## 6 arimax Temperature 2.38 0.622 3.83 1.30e- 4
## 7 arimax HolidayTRUE -17.2 9.76 -1.76 7.82e- 2Diagnose ARIMAX residuals
feasts::gg_tsresiduals(arimax_model)
arimax_residuals <- arimax_model |>
fabletools::augment()
arimax_ljung_box <- arimax_residuals |>
fabletools::features(
.innov,
feasts::ljung_box,
lag = 48,
dof = 4
)
print(arimax_ljung_box)## # A tibble: 1 × 3
## .model lb_stat lb_pvalue
## <chr> <dbl> <dbl>
## 1 arimax 49995. 0Fit both models together for comparison
dynamic_models <- vic_train |>
fabletools::model(
dynamic_regression = fable::TSLM(
Demand ~ Temperature + Holiday + trend() + season()
),
arimax = fable::ARIMA(
Demand ~ Temperature + Holiday
)
)
print(dynamic_models)## # A mable: 1 x 2
## dynamic_regression arimax
## <model> <model>
## 1 <TSLM> <LM w/ ARIMA(2,1,1)(0,0,2)[2] errors>Generate forecasts using known test-set predictors
dynamic_forecasts <- dynamic_models |>
fabletools::forecast(new_data = vic_test)
print(dynamic_forecasts)## # A fable: 672 x 6 [30m] <Australia/Melbourne>
## # Key: .model [2]
## .model Time
## <chr> <dttm>
## 1 dynamic_regression 2014-12-25 00:00:00
## 2 dynamic_regression 2014-12-25 00:30:00
## 3 dynamic_regression 2014-12-25 01:00:00
## 4 dynamic_regression 2014-12-25 01:30:00
## 5 dynamic_regression 2014-12-25 02:00:00
## 6 dynamic_regression 2014-12-25 02:30:00
## 7 dynamic_regression 2014-12-25 03:00:00
## 8 dynamic_regression 2014-12-25 03:30:00
## 9 dynamic_regression 2014-12-25 04:00:00
## 10 dynamic_regression 2014-12-25 04:30:00
## # ℹ 662 more rows
## # ℹ 4 more variables: Demand <dist>, .mean <dbl>, Temperature <dbl>,
## # Holiday <lgl>dynamic_forecasts |>
fabletools::autoplot(vic_train |>
dplyr::slice_tail(n = 7 * 48)) +
ggplot2::autolayer(
vic_test,
Demand,
series = "Actual test demand"
) +
ggplot2::labs(
title = "Dynamic Regression and ARIMAX Forecasts",
subtitle = "Forecasts for the final seven days",
x = "Time",
y = "Electricity Demand"
) +
ggplot2::theme_minimal()
Compare test-set forecast accuracy
dynamic_accuracy <- dynamic_forecasts |>
fabletools::accuracy(vic_test) |>
dplyr::select(.model, RMSE, MAE, MAPE) |>
dplyr::arrange(RMSE)
print(dynamic_accuracy)## # A tibble: 2 × 4
## .model RMSE MAE MAPE
## <chr> <dbl> <dbl> <dbl>
## 1 arimax 670. 575. 16.5
## 2 dynamic_regression 789. 716. 19.8Compare model information criteria
dynamic_information_criteria <- dynamic_models |>
generics::glance() |>
dplyr::select(.model, AIC, AICc, BIC) |>
dplyr::arrange(AICc)
print(dynamic_information_criteria)## # A tibble: 2 × 4
## .model AIC AICc BIC
## <chr> <dbl> <dbl> <dbl>
## 1 arimax 622444. 622444. 622515.
## 2 dynamic_regression 702995. 702995. 703049.Combine forecast accuracy and model-fit results
dynamic_model_comparison <- dynamic_accuracy |>
dplyr::left_join(
dynamic_information_criteria,
by = ".model"
) |>
dplyr::arrange(RMSE)
print(dynamic_model_comparison)## # A tibble: 2 × 7
## .model RMSE MAE MAPE AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 arimax 670. 575. 16.5 622444. 622444. 622515.
## 2 dynamic_regression 789. 716. 19.8 702995. 702995. 703049.#Discussions For the first part of this project, I analyzed quarterly Australian beer production using ARIMA and exponential smoothing models. The original time-series and seasonal plots showed a clear quarterly seasonal pattern, along with changes in the overall production level over time. Seasonal differencing was therefore important for making the series more stationary. After seasonal differencing, the Augmented Dickey-Fuller test produced a statistically significant result, supporting the conclusion that the differenced series was stationary.
The model comparisons showed that the seasonal ARIMA models performed much better than the nonseasonal manual ARIMA model. The automatic procedure selected an ARIMA(1,1,2)(0,1,1)[4] model, which had the lowest AICc of approximately 1,783. It also produced an in-sample RMSE of about 15.2, compared with 50.3 for the nonseasonal ARIMA model. This demonstrates that explicitly accounting for quarterly seasonality greatly improved the fit. On the test set, exponential smoothing performed slightly better than automatic ARIMA, with an RMSE of 9.62 and MAE of 8.92, compared with an RMSE of 9.83 and MAE of 9.04 for ARIMA. The difference was small, suggesting that both approaches captured the main trend and seasonal behavior effectively. ARIMA models seasonality through seasonal differencing and lagged relationships, while exponential smoothing updates the level, trend, and seasonal components directly.
For the second part, I modeled Victorian electricity demand using temperature and holiday indicators as external predictors. The dynamic regression model showed that temperature had a significant positive relationship with electricity demand, while holidays were associated with substantially lower demand. However, the residual plots and Ljung-Box test indicated strong remaining autocorrelation, meaning that the regression model did not fully capture the time-based structure in electricity usage.
The ARIMAX model improved the analysis by combining the predictors with ARIMA errors. It produced a lower test RMSE of approximately 670 and MAE of 575, compared with an RMSE of 789 and MAE of 716 for the dynamic regression model. Its AIC was also considerably lower. These results suggest that temperature and holidays are useful predictors, but electricity demand also depends heavily on its own recent values and recurring patterns. The ARIMAX model performed better because it accounted for both the external influences and the autocorrelation remaining in the regression residuals.