Load package

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()

Get/prep Data

monthly_elec <- vic_elec %>%
  index_by(Month = yearmonth(Time)) %>%
  summarise(
    Demand = sum(Demand),
    Temperature = mean(Temperature)
  )

monthly_elec
## # A tsibble: 36 x 3 [1M]
##       Month   Demand Temperature
##       <mth>    <dbl>       <dbl>
##  1 2012 Jan 7241049.        21.8
##  2 2012 Feb 6874543.        21.4
##  3 2012 Mar 6746560.        18.6
##  4 2012 Apr 6401078.        16.7
##  5 2012 May 7375177.        13.2
##  6 2012 Jun 7388456.        11.0
##  7 2012 Jul 7568114.        11.0
##  8 2012 Aug 7492261.        11.3
##  9 2012 Sep 6567196.        13.7
## 10 2012 Oct 6680705.        15.2
## # ℹ 26 more rows

Visualize Time Series

monthly_elec %>%
  autoplot(Demand) +
  labs(
    title = "Monthly Electricity Demand in Victoria",
    x = "Month",
    y = "Electricity Demand"
  )

check for seasonality

monthly_elec %>%
  gg_season(Demand) +
  labs(
    title = "Monthly Seasonal Pattern in Electricity Demand",
    x = "Month",
    y = "Electricity Demand"
  )

sub series plot

monthly_elec %>%
  gg_subseries(Demand) +
  labs(
    title = "Monthly Subseries Plot of Electricity Demand",
    x = "Month",
    y = "Electricity Demand"
  )

seasonality w/ regular differencing

monthly_elec %>%
  mutate(Diff_Demand = difference(Demand)) %>%
  autoplot(Diff_Demand) +
  labs(
    title = "Differenced Monthly Electricity Demand",
    x = "Month",
    y = "Differenced Demand"
  )
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

seasonal diff

monthly_elec %>%
  mutate(Seasonal_Diff_Demand = difference(Demand, lag = 12)) %>%
  autoplot(Seasonal_Diff_Demand) +
  labs(
    title = "Seasonally Differenced Monthly Electricity Demand",
    x = "Month",
    y = "Demand Difference from Same Month Last Year"
  )
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

KPSS tests

monthly_elec %>%
  features(Demand, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.224         0.1
monthly_elec %>%
  mutate(Diff_Demand = difference(Demand)) %>%
  features(Diff_Demand, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0416         0.1
monthly_elec %>%
  mutate(Seasonal_Diff_Demand = difference(Demand, lag = 12)) %>%
  features(Seasonal_Diff_Demand, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1     0.173         0.1

train/test splits

train <- monthly_elec %>%
  filter(Month <= max(Month) - 12)

test <- monthly_elec %>%
  filter(Month > max(Month) - 12)

Fit the manual ARIMA and auto ARIMA models

fit_arima <- train %>%
  model(
    Manual_ARIMA = ARIMA(Demand ~ pdq(1,1,1)),
    Auto_ARIMA = ARIMA(Demand)
  )

report(fit_arima)
## Warning in report.mdl_df(fit_arima): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
## # A tibble: 2 × 8
##   .model              sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>                <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 Manual_ARIMA 164688625764.   -330.  665.  666.  668. <cpl [1]> <cpl [1]>
## 2 Auto_ARIMA   154319973701.   -342.  690.  692.  694. <cpl [1]> <cpl [0]>

Check Auto ARIMA residuals

fit_arima %>%
  select(Auto_ARIMA) %>%
  gg_tsresiduals()

Check Manual ARIMA residuals

fit_arima %>%
  select(Manual_ARIMA) %>%
  gg_tsresiduals()

Ljung-Box test for Auto ARIMA

fit_arima %>%
  select(Auto_ARIMA) %>%
  augment() %>%
  features(.innov, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model     lb_stat lb_pvalue
##   <chr>        <dbl>     <dbl>
## 1 Auto_ARIMA    22.0    0.0381

Ljung-Box test for Manual ARIMA

fit_arima %>%
  select(Manual_ARIMA) %>%
  augment() %>%
  features(.innov, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 Manual_ARIMA    25.7    0.0119

Forecast the ARIMA models

fc_arima <- fit_arima %>%
  forecast(h = 12)

fc_arima %>%
  autoplot(train, level = NULL) +
  autolayer(test, Demand, color = "black") +
  labs(
    title = "Monthly ARIMA Forecasts for Electricity Demand",
    x = "Month",
    y = "Electricity Demand"
  )

Check ARIMA accuracy

accuracy(fc_arima, test)
## # A tibble: 2 × 10
##   .model       .type       ME    RMSE     MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>        <chr>    <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Auto_ARIMA   Test  -104194. 432882. 386436. -1.92  5.75   NaN   NaN 0.349
## 2 Manual_ARIMA Test   -91303. 431130. 382556. -1.72  5.69   NaN   NaN 0.346

Check ARIMA AIC values

glance(fit_arima)
## # A tibble: 2 × 8
##   .model              sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>                <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 Manual_ARIMA 164688625764.   -330.  665.  666.  668. <cpl [1]> <cpl [1]>
## 2 Auto_ARIMA   154319973701.   -342.  690.  692.  694. <cpl [1]> <cpl [0]>

Fit the dynamic regression model

fit_dynamic <- train %>%
  model(
    Dynamic_Regression = TSLM(Demand ~ Temperature)
  )

report(fit_dynamic)
## Series: Demand 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -566777 -348389   71960  279699  656064 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7650759     341379  22.411   <2e-16 ***
## Temperature   -48814      20553  -2.375   0.0267 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 381800 on 22 degrees of freedom
## Multiple R-squared: 0.2041,  Adjusted R-squared: 0.1679
## F-statistic: 5.641 on 1 and 22 DF, p-value: 0.026686

Check dynamic regression residuals

fit_dynamic %>%
  gg_tsresiduals()

Ljung-Box test for dynamic regression

fit_dynamic %>%
  augment() %>%
  features(.innov, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model             lb_stat lb_pvalue
##   <chr>                <dbl>     <dbl>
## 1 Dynamic_Regression    18.9    0.0917

Forecast the dynamic regression model

future_xreg <- test %>%
  select(Month, Temperature)

fc_dynamic <- fit_dynamic %>%
  forecast(new_data = future_xreg)

fc_dynamic %>%
  autoplot(train, level = NULL) +
  autolayer(test, Demand, color = "black") +
  labs(
    title = "Monthly Dynamic Regression Forecast",
    x = "Month",
    y = "Electricity Demand"
  )

Check dynamic regression accuracy

accuracy(fc_dynamic, test)
## # A tibble: 1 × 10
##   .model             .type       ME    RMSE    MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>              <chr>    <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Dynamic_Regression Test  -113164. 372349. 3.25e5 -1.97  4.86   NaN   NaN 0.190

Fit the ARIMAX model

fit_arimax <- train %>%
  model(
    ARIMAX = ARIMA(Demand ~ Temperature)
  )

report(fit_arimax)
## Series: Demand 
## Model: LM w/ ARIMA(0,0,0) errors 
## 
## Coefficients:
##       Temperature  intercept
##         -48814.17  7650759.0
## s.e.     19681.16   326907.4
## 
## sigma^2 estimated as 1.458e+11:  log likelihood=-341.48
## AIC=688.95   AICc=690.15   BIC=692.49

Check ARIMAX residuals

fit_arimax %>%
  gg_tsresiduals()

Ljung-Box test for ARIMAX

fit_arimax %>%
  augment() %>%
  features(.innov, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMAX    18.9    0.0917

Forecast the ARIMAX model

fc_arimax <- fit_arimax %>%
  forecast(new_data = future_xreg)

fc_arimax %>%
  autoplot(train, level = NULL) +
  autolayer(test, Demand, color = "black") +
  labs(
    title = "Monthly ARIMAX Forecast",
    x = "Month",
    y = "Electricity Demand"
  )

Check ARIMAX accuracy

accuracy(fc_arimax, test)
## # A tibble: 1 × 10
##   .model .type       ME    RMSE     MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>    <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMAX Test  -113164. 372349. 325272. -1.97  4.86   NaN   NaN 0.190

Compare all of the forecasts

fc_compare <- bind_rows(
  fc_arima,
  fc_dynamic,
  fc_arimax
)

fc_compare %>%
  autoplot(train, level = NULL) +
  autolayer(test, Demand, color = "black") +
  labs(
    title = "Monthly Forecast Comparison: ARIMA, Dynamic Regression, and ARIMAX",
    x = "Month",
    y = "Electricity Demand"
  )

Compare all model accuracy

accuracy(fc_compare, test)
## # A tibble: 4 × 10
##   .model             .type       ME    RMSE    MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>              <chr>    <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMAX             Test  -113164. 372349. 3.25e5 -1.97  4.86   NaN   NaN 0.190
## 2 Auto_ARIMA         Test  -104194. 432882. 3.86e5 -1.92  5.75   NaN   NaN 0.349
## 3 Dynamic_Regression Test  -113164. 372349. 3.25e5 -1.97  4.86   NaN   NaN 0.190
## 4 Manual_ARIMA       Test   -91303. 431130. 3.83e5 -1.72  5.69   NaN   NaN 0.346

Compare model information values

glance(fit_arima)
## # A tibble: 2 × 8
##   .model              sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>                <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 Manual_ARIMA 164688625764.   -330.  665.  666.  668. <cpl [1]> <cpl [1]>
## 2 Auto_ARIMA   154319973701.   -342.  690.  692.  694. <cpl [1]> <cpl [0]>
glance(fit_dynamic)
## # A tibble: 1 × 15
##   .model   r_squared adj_r_squared  sigma2 statistic p_value    df log_lik   AIC
##   <chr>        <dbl>         <dbl>   <dbl>     <dbl>   <dbl> <int>   <dbl> <dbl>
## 1 Dynamic…     0.204         0.168 1.46e11      5.64  0.0267     2   -341.  621.
## # ℹ 6 more variables: AICc <dbl>, BIC <dbl>, CV <dbl>, deviance <dbl>,
## #   df.residual <int>, rank <int>
glance(fit_arimax)
## # A tibble: 1 × 8
##   .model        sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 ARIMAX 145794570454.   -341.  689.  690.  692. <cpl [0]> <cpl [0]>

ets model

fit_ets <- train %>%
  model(
    ETS = ETS(Demand)
  )

fc_ets <- fit_ets %>%
  forecast(h = 12)

fc_ets %>%
  autoplot(train, level = NULL) +
  autolayer(test, Demand, color = "black") +
  labs(
    title = "ETS Forecast for Monthly Electricity Demand",
    x = "Month",
    y = "Electricity Demand"
  )

accuracy(fc_ets, test)
## # A tibble: 1 × 10
##   .model .type      ME    RMSE     MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>   <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ETS    Test  317583. 515704. 370265.  4.38  5.23   NaN   NaN 0.389

total forecasts

model_accuracy <- bind_rows(
  accuracy(fc_arima, test),
  accuracy(fc_dynamic, test),
  accuracy(fc_arimax, test),
  accuracy(fc_ets, test)
)

model_accuracy %>%
  select(.model, RMSE, MAE, MAPE)
## # A tibble: 5 × 4
##   .model                RMSE     MAE  MAPE
##   <chr>                <dbl>   <dbl> <dbl>
## 1 Auto_ARIMA         432882. 386436.  5.75
## 2 Manual_ARIMA       431130. 382556.  5.69
## 3 Dynamic_Regression 372349. 325272.  4.86
## 4 ARIMAX             372349. 325272.  4.86
## 5 ETS                515704. 370265.  5.23

For this discussion, I used the vic_elec dataset. This dataset includes electricity data from Victoria, Australia, along with temperature data. I chose this dataset because demand is a good example of something affected by both seasonal patterns and external influences. Since the original dataset is half-hourly, I converted it to monthly data. This will make it easier to see if there is a trend by month. After plotting, the monthly seasonal plot showed that electricity demand varies throughout the year; it appears to rise during warmer months and fall during cooler ones. This supports treating temperature as an external predictor, since electricity use is likely affected by seasonal weather. Before fitting ARIMA models, I checked for stationarity using regular and seasonal differencing and the KPSS test. Regular differencing was performed to examine changes between months. Seasonal differencing only compared each month to the same month in previous years. This helped determine whether the series needed stabilization. I then fit both a manual ARIMA model and an automatic ARIMA model. The manual model was kept simple as ARIMA(1,1,1) because the monthly dataset has a limited number of observations. The automatic ARIMA model was useful because it selected a model structure based on the data. I compared the models using residual diagnostics, Ljung-Box tests, AIC values, and forecast accuracy. For dynamic regression, I used temperature to predict electricity demand. This helped explain demand using an external variable, but it did not fully account for autocorrelation. Because of that, I also fit an ARIMAX model with ARIMA(Demand ~ Temperature). This model combines temperature as a predictor with ARIMA errors. Overall, I found that ARIMAX is the strongest approach because electricity demand depends on both temperature and past time-series patterns. The accuracy result showed that dynamic regressions and ARIMAX models performed the best overall by having lower RMSE, MAE and MAPE values. This suggests adding temperature to the forecast improved it than only using demand patterns. ETS performed better than the ARIMA models, but only by MAPE. The main limitation is that the monthly dataset is small, so the results may not be entirely accurate.