Executive summary

The goal of this project was to forecast cash withdrawals from four ATMs (part A) and residential power usage (part B). For each time series, I first visualized the data and checked for missing observations, outliers, and non-zero values, which is a requirement for Box-Cox transformation and multiplicative ETS models. Missing data and outliers were interpolated using an ARIMA model. Then, I characterized the time series (eg, trend, seasonality) using STL decomposition. I also assessed whether variance was non-constant to determine the need for Box-Cox transformation and performed KPSS unit root tests to determine whether differencing was needed to achieve stationarity.

For all time series except ATM3, which only had 3 observations, I split the data into training and test sets and compared the performance of ETS and ARIMA models with a baseline naive model using four out-of-sample accuracy metrics (RMSE, MAE, MPE, and MAPE). I used residual plots and Ljung-Box tests to determine whether the residuals were white noise, which indicates that the model fit the data well. For part B, finding an appropriate model required two iterations of this process.

Key operations and the best model for each dataset are summarized below:

  • ATM1: Box-Cox transformation, first-order seasonal differencing, \(ETS(A,A,A)\) model

  • ATM2: Box-Cox transformation, first-order seasonal differencing, \(ARIMA(0,0,0)(1,1,1)_7 + drift\) model

  • ATM3: \(Mean\) model

  • ATM4: Box-Cox transformation, \(ARIMA(2,0,0)(2,0,0)_7 + mean\) model

  • Residential power: First-order seasonal differencing, \(ARIMA(1,0,1)(2,1,0)_{12} + drift\) model

Finally, I refit each model to the entire time series (ie, training + test) to forecast future values and saved the forecasts as Excel files.

A. ATM forecasts

Data input and wrangling

filename <- 'ATM624Data.xlsx'
url <- paste0(domain, path, filename)

atm_df <- read.xlsx(url, sheet = 1)

The raw dataset is a tibble (dataframe) with 3 variables and 1474 observations, which is composed of 365 observations for each of the 4 ATMs.

class(atm_df)
## [1] "data.frame"
glimpse(atm_df)
## Rows: 1,474
## Columns: 3
## $ DATE <dbl> 39934, 39934, 39935, 39935, 39936, 39936, 39937, 39937, 39938, 39…
## $ ATM  <chr> "ATM1", "ATM2", "ATM1", "ATM2", "ATM1", "ATM2", "ATM1", "ATM2", "…
## $ Cash <dbl> 96, 107, 82, 89, 85, 90, 90, 55, 99, 79, 88, 19, 8, 2, 104, 103, …
table(atm_df$ATM)
## 
## ATM1 ATM2 ATM3 ATM4 
##  365  365  365  365

The date variable is numeric and requires conversion to a date-time data type.

# Origin date for Excel
# https://www.r-bloggers.com/2024/02/taming-excel-dates-in-r-from-numbers-to-meaningful-dates/
atm_df <- atm_df %>%
  mutate(
    Date = as.Date(DATE, origin = '1899-12-30'),
    .after = DATE
  ) %>%
  select(-c(DATE))

Finally, I subset the main dataframe into 4 dataframes by ATM number and converted each to a tsibble.

# Split main dataframe into list of dataframes by ATM number
atm_df_list <- split(atm_df, atm_df$ATM)

# Convert each dataframe to tsibble
atm_ts_list <- lapply(atm_df_list, function(df) as_tsibble(df, index = Date))

# Remove ATM variable (no longer needed) from each tsibble
atm_ts_list <- lapply(atm_ts_list, function(ts) select(ts, -ATM))

# Create individual objects for each element in list
# This creates 4 tsibbles (ATM1, ATM2, ATM3, ATM4)
invisible(list2env(atm_ts_list, envir = .GlobalEnv))

Glimpsing each ATM dataframe revealed that the Cash values for ATM4 are floats, whereas the values for the other ATMs are whole numbers. For consistency, and because ATMs do not dispense coins, I rounded down all Cash values for ATM4 to the nearest whole number. However, in the real world, it would be best to contact the data provider to understand why the format of the Cash values are different.

glimpse(ATM1)
## Rows: 365
## Columns: 2
## $ Date <date> 2009-05-01, 2009-05-02, 2009-05-03, 2009-05-04, 2009-05-05, 2009…
## $ Cash <dbl> 96, 82, 85, 90, 99, 88, 8, 104, 87, 93, 86, 111, 75, 6, 102, 73, …
glimpse(ATM2)
## Rows: 365
## Columns: 2
## $ Date <date> 2009-05-01, 2009-05-02, 2009-05-03, 2009-05-04, 2009-05-05, 2009…
## $ Cash <dbl> 107, 89, 90, 55, 79, 19, 2, 103, 107, 118, 75, 111, 25, 16, 137, …
glimpse(ATM3)
## Rows: 365
## Columns: 2
## $ Date <date> 2009-05-01, 2009-05-02, 2009-05-03, 2009-05-04, 2009-05-05, 2009…
## $ Cash <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
glimpse(ATM4)
## Rows: 365
## Columns: 2
## $ Date <date> 2009-05-01, 2009-05-02, 2009-05-03, 2009-05-04, 2009-05-05, 2009…
## $ Cash <dbl> 776.99342, 524.41796, 792.81136, 908.23846, 52.83210, 52.20845, 5…
# Change ATM4 cash values to whole numbers
ATM4 <- ATM4 %>%
  mutate(
    Cash = floor(Cash)
  )

Initial observations

The facet plot below shows that ATM3 and ATM4 have different data characteristics than ATM1 and ATM2 that affect overall data analysis strategies.

Both ATM1 and ATM2 show time-dependent changes in cash withdrawal activity, which can be analyzed using either simple (eg, naive) or advanced forecasting methods (eg, ETS, ARIMA).

ATM3 only had activity during the last three days of the observation period (perhaps it is a recently installed ATM). Because of this, there is insufficient data to use exponential smoothing or ARIMA, as well as to split the data into training and test sets to evaluate model performance. So, only simple methods such as the mean or naive method can be used for forecasting.

filter(ATM3, Cash > 0)
## # A tsibble: 3 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2010-04-28    96
## 2 2010-04-29    82
## 3 2010-04-30    85

ATM4 has at least one extreme outlier. In practice, it would be a good idea to double check the data source (ie, perhaps there was a typo). However, for this exercise, the outlier(s) can be removed and replaced with interpolated values using an ARIMA model before model development.

# Combined ATM data and coerce ATM as factor
ATM_combined <- bind_rows(atm_df_list)
ATM_combined$ATM <- factor(ATM_combined$ATM, levels = c('ATM1', 'ATM2', 'ATM3', 'ATM4'))

# Facet plot
ggplot(ATM_combined, aes(x = Date, y = Cash)) +
  geom_line() +
  facet_wrap(~ ATM, scales = 'free_y', ncol = 2, nrow = 2) + 
  labs(
    y = 'Amount (USD, x100)',
    title = 'ATM cash withdrawals, May 2009 to Apr 2010'
  )

ATM1

Data checks

Missing data

There were three missing data points, all from June 2009.

(ATM1_missing <- filter(ATM1, is.na(Cash)))
## # A tsibble: 3 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2009-06-13    NA
## 2 2009-06-16    NA
## 3 2009-06-22    NA

Outliers

There were no outliers that exceeded 3 interquartile ranges.

# Boundaries for outlier detection
lower_fence <- quantile(ATM1$Cash, 0.25, na.rm = TRUE) - 3 * IQR(ATM1$Cash, na.rm = TRUE)
upper_fence <- quantile(ATM1$Cash, 0.75, na.rm = TRUE) + 3 * IQR(ATM1$Cash, na.rm = TRUE)

ATM1_outliers <- ATM1 %>%
  filter(Cash < lower_fence | Cash > upper_fence)

ATM1_outliers
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Negative or zero values

All Cash values were positive, which satisfies the requirement for Box-Cox transformation if it is needed.

filter(ATM1, Cash <= 0)
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Visualize time series

The plot below shows the ATM1 time series with missing observations and boundaries for outlier detection.

ATM1_plot <- ggtime::autoplot(ATM1, Cash) + 
  # Highlight outlier boundaries
  geom_hline(yintercept = lower_fence, linetype = 'dashed', color = 'goldenrod') +
  geom_hline(yintercept = upper_fence, linetype = 'dashed', color = 'goldenrod') +  
  labs(
    y = 'Amount (USD, x100)',
    title = 'ATM1 cash withdrawals, May 2009 to Apr 2010',
    subtitle = 'Raw data has 3 missing observations and no outliers'    
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(face = 'bold'),
    plot.title = element_text(face = 'bold')
  )

# Add vertical lines to highlight missing data points
for (i in 1:nrow(ATM1_missing)) {
  ATM1_plot <- ATM1_plot + 
    geom_vline(xintercept = ATM1_missing[[i, 'Date']], 
               linetype = 'dashed', color = 'steelblue')
}

ATM1_plot

Data interpolation

I replaced the missing data with estimates interpolated with an ARIMA model.

ATM1_interpolated <- ATM1 %>%
  # Fit ARIMA model to data containing missing values
  model(ARIMA(Cash)) %>%
  # Estimate Cash for all time points
  interpolate(ATM1)

STL decomposition

STL decomposition of the interpolated time series showed that cash withdrawals from ATM1 have a downward trend that appears to increase with the level of the series as well as seasonality with a weekly period. Together, this indicates that the time series is not stationary, and that an additive or multiplicative ETS or ARIMA model could be suitable.

Of note, the remainder component shows an increase in noise after mid-February 2010. This suggests that the trend and seasonal components do not fully account for changes in the data after this point. This may also indicate a need for a multiplicative error component in an ETS model.

ATM1_decomp <- ATM1_interpolated %>%
  model(STL(Cash ~ season(window = 'periodic'), robust = TRUE)) %>%
  components()

ATM1_decomp %>%
  ggtime::autoplot() +
    # Unomment next line to view a more limited time period  
    # coord_cartesian(xlim = as.Date(c('2010-01-01', '2010-04-30'))) +
    ggtitle('STL decomposition of interpolated ATM1 withdrawals')

Variance and autocorrelation

  • The time series plot (top) shows non-constant variance and the cash withdrawal amounts are all positive, so Box-Cox transformation could be helpful.

  • The ACF plot (lower left) appears to be somewhat sinusoidal and has three significant spikes that decrease in magnitude as lag increases, which is consistent with seasonality and non-stationarity.

  • The partial ACF plot (lower right) has several significant spikes up to lag=14, which suggests that a higher-order seasonal ARIMA model may be needed.

ATM1_interpolated %>%
  ggtime::gg_tsdisplay(Cash, plot_type = 'partial') +
    # I plotted a limited time period to make it easier to see changes in variance
    # To view the entire time series, comment out the next line
    coord_cartesian(xlim = as.Date(c('2009-10-01', '2010-01-31'))) +
    ggtitle('ATM1 withdrawals (interpolated)')

Box-Cox transformation

The Guerrero method indicated that a Box-Cox transformation with \(\lambda = 0.26\) would help stabilize the variance.

# Determine optimal value of lambda for Box-Cox transformation
ATM1_lambda <- ATM1_interpolated %>%
  features(Cash, features = guerrero) %>%
  pull(lambda_guerrero)

ATM1_lambda
## [1] 0.2622969

Differencing

The KPSS unit root test indicated that first-order seasonal differencing of the transformed data is sufficient to achieve stationarity.

ATM1_interpolated %>%
  features(box_cox(Cash, ATM1_lambda), unitroot_nsdiffs) %>%
  pull(nsdiffs)
## [1] 1

Models

Development

The data characteristics described above suggest that additive or multiplicative ETS and seasonal ARIMA models can be fit to the transformed data. I compared the performance of these models with the performance of a seasonal naive model as baseline.

# Create training and test datasets to develop models and assess performance
# I split by year for convenience. By observations, the split ratio is ~65:35.
ATM1_train <- ATM1_interpolated %>%
  filter(year(Date) == 2009)

ATM1_test <- ATM1_interpolated %>%
  filter(year(Date) == 2010)
# ETS
ATM1_train_ets_fit <- ATM1_train %>%
  # Automatically determine the best method from the specified options
  # The transformation is specified inside the model so fable will automatically
  # reverse it during forecasting
  model(ETS = ETS(box_cox(Cash, ATM1_lambda) ~ 
                            error(method = c('A', 'M')) +
                            trend(method = c('A', 'Ad', 'M')) +
                            season(method = c('A', 'M'))))

# ARIMA
ATM1_train_arima_fit <- ATM1_train %>%
  # specify first-order difference (D), automatic determination of P and Q
  model(ARIMA = ARIMA(box_cox(Cash, ATM1_lambda) ~ PDQ(D = 1)))   

# Baseline seasonal naive model
ATM1_train_snaive_fit = ATM1_train %>%
  model(SNAIVE = SNAIVE(box_cox(Cash, ATM1_lambda)))
model_list <- list(ATM1_train_ets_fit, 
                   ATM1_train_arima_fit,
                   ATM1_train_snaive_fit)

compare_forecast_accuracy(model_list, ATM1_test) %>%
  arrange(RMSE)
## # A tibble: 3 × 5
##   model   RMSE   MAE   MPE  MAPE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 ETS     44.0  32.9 -211.  261.
## 2 ARIMA   44.1  30.1 -283.  314.
## 3 SNAIVE  86.3  73.1 -331.  368.

Comparison

The \(ETS(A,A,A)\) model had the lowest RMSE, MPE, and MAPE, while the \(ARIMA(0,0,1)(0,1,1)_7\) model had the lowest MAE. Due to the different model classes, comparisons using AIC or BIC are not valid.1

Details of each model are shown below.

report(ATM1_train_arima_fit)
## Series: Cash 
## Model: ARIMA(0,0,1)(0,1,1)[7] 
## Transformation: box_cox(Cash, ATM1_lambda) 
## 
## Coefficients:
##          ma1     sma1
##       0.1331  -0.9243
## s.e.  0.0652   0.0453
## 
## sigma^2 estimated as 1.486:  log likelihood=-390.57
## AIC=787.13   AICc=787.23   BIC=797.55
report(ATM1_train_ets_fit)
## Series: Cash 
## Model: ETS(A,A,A) 
## Transformation: box_cox(Cash, ATM1_lambda) 
##   Smoothing parameters:
##     alpha = 0.05922859 
##     beta  = 0.0001000784 
##     gamma = 0.0001362675 
## 
##   Initial states:
##      l[0]         b[0]      s[0]      s[-1]    s[-2]     s[-3]     s[-4]
##  7.933366 -0.002664309 -4.160148 0.00941567 1.206887 0.5304338 0.9595673
##      s[-5]    s[-6]
##  0.8724516 0.581393
## 
##   sigma^2:  1.4368
## 
##      AIC     AICc      BIC 
## 1449.346 1450.691 1491.361
report(ATM1_train_snaive_fit)
## Series: Cash 
## Model: SNAIVE 
## Transformation: box_cox(Cash, ATM1_lambda) 
## 
## sigma^2: 2.6367

Residual diagnostics

The residual plots of the \(ETS(A,A,A)\) model shows several negative spikes, but they do not appear to have any pattern. All of the ACF spikes are non-significant. The distribution of residuals is close to normal (although there are a few outliers), and the mean is near zero. In addition, the Ljung-Box test P-value is not significant, so the null hypothesis that the residuals are white noise is not rejected. Collectively, these results suggest that this model is a good fit to the data.

ggtime::gg_tsresiduals(ATM1_train_ets_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a weekly period, so l = 2m = 2*7 = 14.
# Degrees of freedom = number of parameters in model, so dof = 3
ATM1_train_ets_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 14, dof = 3)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ETS       4.26     0.962

The residual plots of the \(ARIMA(0,0,1)(0,1,1)_7\) model shows several negative spikes, but they do not appear to have any pattern. All of the ACF spikes are non-significant. The distribution of residuals is close to normal (although there are a few outliers), and the mean is approximately zero. In addition, the Ljung-Box test P-value is not significant, so the null hypothesis that the residuals are white noise is not rejected. Collectively, these results suggest that this model is a good fit to the data.

ggtime::gg_tsresiduals(ATM1_train_arima_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a weekly period, so l = 2m = 2*7 = 14.
# Degrees of freedom = number of parameters in model, so dof = 2
ATM1_train_arima_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 14, dof = 2)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     6.54     0.887

Forecast future values

Although both the ETS and ARIMA models were deemed appropriate, I selected the seasonal \(ETS(A,A,A)\) model as the final model because it had the lowest of three accuracy metrics (including RMSE, which is a key metric) whereas the ARIMA model only had one. In addition, the ETS model is simpler.

To forecast future values, I refit this model to the entire dataset (ie, train + test), and then used it to forecast cash withdrawals 1 month into the future.

# Refit the training model to the full dataset
ATM1_final_ets_fit <- ATM1_train_ets_fit %>%
  fabletools::refit(ATM1_interpolated)
# Forecast future values
ATM1_forecast <- ATM1_final_ets_fit %>%
  forecast(h = '1 month')

The forecast looks reasonable, although the upper range of the 95% prediction interval is a little high relative to the historical data. This may reflect the increased noise in the data near the end of the observation period, which was noted in the STL decomposition. More advanced techniques, such as machine learning, may be needed to improve the forecast.

# Plot forecast
ATM1_forecast %>%
  ggtime::autoplot(ATM1_interpolated) +
    labs(
      y = 'Amount (USD, x100)',
      title = 'ATM1 cash withdrawals, May 2009 to Apr 2010 + Forecast',
      subtitle = 'ETS(A,A,A) model'    
    ) +
    theme_classic() +
    theme(
      axis.title = element_text(face = 'bold'),
      plot.title = element_text(face = 'bold')
    )

Save forecast

# Note this block only creates a dataframe with forecasts;
# these will be consolidated and exported at the end of the project
ATM1_forecast <- ATM1_forecast %>%
  as.data.frame() %>%
  select(Date, Cash_forecast = `.mean`) %>%
  # Round forecasts since ATMs don't dispense coins
  mutate(
    Cash_forecast = floor(Cash_forecast)
  )

# Workaround because Excel doesn't read date field correctly
ATM1_forecast$Date <- as.character(ATM1_forecast$Date)

ATM2

Data checks

Missing data

There were two missing data points, both from June 2009.

(ATM2_missing <- filter(ATM2, is.na(Cash)))
## # A tsibble: 2 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2009-06-18    NA
## 2 2009-06-24    NA

Outliers

There were no outliers that exceeded 3 interquartile ranges.

# Boundaries for outlier detection
lower_fence <- quantile(ATM2$Cash, 0.25, na.rm = TRUE) - 3 * IQR(ATM2$Cash, na.rm = TRUE)
upper_fence <- quantile(ATM2$Cash, 0.75, na.rm = TRUE) + 3 * IQR(ATM2$Cash, na.rm = TRUE)

ATM2_outliers <- ATM2 %>%
  filter(Cash < lower_fence | Cash > upper_fence)

ATM2_outliers
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Negative or zero values

Two Cash values were zero. This precludes Box-Cox transformation as well as multiplicative ETS methods; however, a simple workaround is to shift the values by adding a constant (this will be done during model development).

filter(ATM2, Cash <= 0)
## # A tsibble: 2 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2009-10-25     0
## 2 2010-03-30     0

Visualize time series

The plot below shows the ATM2 time series with missing observations and boundaries for outlier detection.

ATM2_plot <- ggtime::autoplot(ATM2, Cash) + 
  # Highlight outlier boundaries
  geom_hline(yintercept = lower_fence, linetype = 'dashed', color = 'goldenrod') +
  geom_hline(yintercept = upper_fence, linetype = 'dashed', color = 'goldenrod') +  
  labs(
    y = 'Amount (USD, x100)',
    title = 'ATM2 cash withdrawals, May 2009 to Apr 2010',
    subtitle = 'Raw data has 2 missing observations and no outliers'
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(face = 'bold'),
    plot.title = element_text(face = 'bold')
  )

# Add vertical lines to highlight missing data points
for (i in 1:nrow(ATM2_missing)) {
  ATM2_plot <- ATM2_plot + 
    geom_vline(xintercept = ATM2_missing[[i, 'Date']], 
               linetype = 'dashed', color = 'steelblue')
}

ATM2_plot

Data interpolation

I replaced the missing data with estimates interpolated with an ARIMA model.

ATM2_interpolated <- ATM2 %>%
  # Fit ARIMA model to data containing missing values
  model(ARIMA(Cash)) %>%
  # Estimate Cash for all time points
  interpolate(ATM2)

STL decomposition

STL decomposition of the interpolated time series showed that cash withdrawals from ATM2 do not appear to have a trend, but does have seasonality with a weekly period. Together, this indicates that the time series is not stationary, and that an ETS or ARIMA model could be suitable.

Of note, the remainder component shows some residual structure (ie, not white noise) after mid-February, which suggests that the trend and seasonal components do not fully account for changes in the data after this time point.

ATM2_decomp <- ATM2_interpolated %>%
  model(STL(Cash ~ season(window = 'periodic'), robust = TRUE)) %>%
  components()

ATM2_decomp %>%
  ggtime::autoplot() +
    # Uncomment next line to view a more limited period  
    # coord_cartesian(xlim = as.Date(c('2010-01-01', '2010-04-30'))) +
    ggtitle('STL decomposition of interpolated ATM2 withdrawals')

Variance and autocorrelation

  • The time series plot (top) shows non-constant variance and the (shifted) cash withdrawal amounts are all positive, so Box-Cox transformation could be helpful.

  • The ACF plot (lower left) has several significant spikes that slowly decrease in magnitude as lag increases, which is consistent with seasonality and non-stationarity.

  • The partial ACF plot (lower right) has several significant spikes up to lag=14, which suggests that a higher-order seasonal ARIMA model may be needed.

ATM2_interpolated %>%
  ggtime::gg_tsdisplay(Cash, plot_type = 'partial') +
    # I plotted a limited time period to make it easier to see changes in variance
    # To view the entire time series, comment out the next line
    coord_cartesian(xlim = as.Date(c('2009-10-01', '2010-01-31'))) +
    ggtitle('ATM2 withdrawals (interpolated)')

Box-Cox transformation

The Guerrero method indicated that a Box-Cox transformation with \(\lambda = 0.62\) would help stabilize the variance.

# Determine optimal value of lambda for Box-Cox transformation
ATM2_lambda <- ATM2_interpolated %>%
  features(Cash, features = guerrero) %>%
  pull(lambda_guerrero)

ATM2_lambda
## [1] 0.6746523

Differencing

The KPSS unit root test indicated that first-order seasonal differencing of the transformed data is sufficient to achieve stationarity.

ATM2_interpolated %>%
  features(box_cox(Cash, ATM2_lambda), unitroot_nsdiffs) %>%
  pull(nsdiffs)
## [1] 1

Models

Development

The data characteristics described above suggest that ETS and seasonal ARIMA models can be fit to the transformed data. I compared the performance of these models with the performance of a seasonal naive model as baseline.

# Create training and test datasets to develop models and assess performance
ATM2_train <- ATM2_interpolated %>%
  filter(year(Date) == 2009)

ATM2_test <- ATM2_interpolated %>%
  filter(year(Date) == 2010)
# ETS
ATM2_train_ets_fit <- ATM2_train %>%
  # Automatically determine the best method from the specified options
  # A constant is added to the data to ensure that Box-Cox transformation is
  # performed on positive values
  model(ETS = ETS(box_cox((Cash + 1), ATM2_lambda) ~ 
                    error(method = c('A', 'M')) +
                    trend(method = c('A', 'Ad', 'N')) +
                    season('A')))

# ARIMA
ATM2_train_arima_fit <- ATM2_train %>%
  # specify first-order difference (D), automatic determination of P and Q
  model(ARIMA = ARIMA(box_cox((Cash + 1), ATM2_lambda) ~ PDQ(D = 1)))   

# Baseline seasonal naive model
ATM2_train_snaive_fit = ATM2_train %>%
  model(SNAIVE = SNAIVE(box_cox((Cash + 1), ATM2_lambda)))
model_list <- list(ATM2_train_ets_fit, 
                   ATM2_train_arima_fit,
                   ATM2_train_snaive_fit)

compare_forecast_accuracy(model_list, ATM2_test) %>%
  arrange(RMSE)
## # A tibble: 3 × 5
##   model   RMSE   MAE   MPE  MAPE
##   <chr>  <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA   48.9  39.2  -Inf   Inf
## 2 ETS     49.2  39.3  -Inf   Inf
## 3 SNAIVE  51.6  41.1  -Inf   Inf

Comparison

MPE and MAPE are not usable metrics, most likely because the original data had two zero values. The model with the lowest RMSE and MAE is the \(ARIMA(0,0,0)(1,1,1)_7 + drift\) model. Due to the different model classes, comparisons using AIC or BIC are not valid.

Details of each model are shown below.

report(ATM2_train_arima_fit)
## Series: Cash 
## Model: ARIMA(0,0,0)(1,1,1)[7] w/ drift 
## Transformation: box_cox((Cash + 1), ATM2_lambda) 
## 
## Coefficients:
##          sar1     sma1  constant
##       -0.0731  -0.8956   -0.1824
## s.e.   0.0729   0.0380    0.0551
## 
## sigma^2 estimated as 34.42:  log likelihood=-763.42
## AIC=1534.84   AICc=1535.01   BIC=1548.73
report(ATM2_train_ets_fit)
## Series: Cash 
## Model: ETS(A,A,A) 
## Transformation: box_cox((Cash + 1), ATM2_lambda) 
##   Smoothing parameters:
##     alpha = 0.0001000339 
##     beta  = 0.0001000003 
##     gamma = 0.0001000695 
## 
##   Initial states:
##     l[0]        b[0]      s[0]     s[-1]    s[-2]    s[-3]   s[-4]    s[-5]
##  24.6899 -0.02201866 -17.38408 -8.869069 7.778726 1.990136 1.43907 4.962735
##     s[-6]
##  10.08249
## 
##   sigma^2:  33.847
## 
##      AIC     AICc      BIC 
## 2223.407 2224.752 2265.422
report(ATM2_train_snaive_fit)
## Series: Cash 
## Model: SNAIVE 
## Transformation: box_cox((Cash + 1), ATM2_lambda) 
## 
## sigma^2: 65.0693

Residual diagnostics

The residual plots of the \(ARIMA(0,0,0)(1,1,1)_7 + drift\) model shows several negative spikes, but they do not appear to have a clear pattern. All of the ACF spikes are non-significant. The distribution of residuals is close to normal (although there are a few outliers), and the mean is approximately zero. In addition, the Ljung-Box test P-value is not significant, so the null hypothesis that the residuals are white noise is not rejected. Collectively, these results suggest that this model is a good fit to the data.

ggtime::gg_tsresiduals(ATM2_train_arima_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a weekly period, so l = 2m = 2*7 = 14.
# Degrees of freedom = number of parameters in model, so dof = 2
ATM2_train_arima_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 14, dof = 2)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     8.39     0.754

Forecast future values

On the basis of the above results, I selected the seasonal \(ARIMA(0,0,0)(1,1,1)_7 + drift\) model as the final model. To forecast future values, I refit this model to the entire dataset (ie, train + test), and then used it to forecast cash withdrawals 1 month into the future.

# Refit the training model to the full dataset
ATM2_final_arima_fit <- ATM2_train_arima_fit %>%
  fabletools::refit(ATM2_interpolated)
# Forecast future values
ATM2_forecast <- ATM2_final_arima_fit %>%
  forecast(h = '1 month')

The forecast looks reasonable.

# Plot forecast
ATM2_forecast %>%
  ggtime::autoplot(ATM2_interpolated) +
    labs(
      y = 'Amount (USD, x100)',
      title = 'ATM2 cash withdrawals, May 2009 to Apr 2010 + Forecast',
      subtitle = latex2exp::TeX('ARIMA(0,0,0)(1,1,1)$_{7}$ model')    
    ) +
    theme_classic() +
    theme(
      axis.title = element_text(face = 'bold'),
      plot.title = element_text(face = 'bold')
    )

Save forecast

ATM2_forecast <- ATM2_forecast %>%
  as.data.frame() %>%
  select(Date, Cash_forecast = `.mean`) %>%
  # Round forecasts since ATMs don't dispense coins
  mutate(
    Cash_forecast = floor(Cash_forecast)
  )

# Workaround because Excel doesn't read date field correctly
ATM2_forecast$Date <- as.character(ATM2_forecast$Date)

ATM3

Data checks

Missing data

There is no missing data point in this time series.

filter(ATM3, is.na(Cash))
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Models

Mean

As noted earlier and shown below, this ATM only had 3 days of cash withdrawal activity. Given that the withdrawal amounts are similar, forecasting with a simple mean is sensible.

filter(ATM3, Cash > 0)
## # A tsibble: 3 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2010-04-28    96
## 2 2010-04-29    82
## 3 2010-04-30    85
# Mean model
ATM3_mean_fit <- ATM3 %>%
  filter(Cash > 0) %>%
  model(MEAN(Cash))
# Forecast 1 month into future
ATM3_mean_forecast <- ATM3_mean_fit %>%
  forecast(h = '1 month')
# Visualize forecast
ATM3_mean_forecast %>%
  ggtime::autoplot(ATM3) +
    # Only plot subset of time period for visibility
    coord_cartesian(xlim = as.Date(c('2010-04-20', '2010-05-30'))) +  
    labs(
      y = 'Amount (USD, x100)',
      title = 'ATM3 cash withdrawals, Apr 28-30 + Forecast',
      subtitle = 'MEAN model'
    ) +
  theme_classic() +
  theme(
    plot.title = element_text(face = 'bold'),
    axis.title = element_text(face = 'bold')
  )

Naive

Alternatively, a naive model, which uses the value of the last observation, could be used. This approach is also reasonable because the amount of the last withdrawal is close to the mean of the three withdrawals. However, the prediction interval is much wider than the forecast using the mean model.

# last withdrawal amount
slice_tail(ATM3, n = 1)
## # A tsibble: 1 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2010-04-30    85
# mean withdrawal amount
ATM3 %>%
  as.data.frame() %>%
  filter(Cash > 0) %>%
  summarise(
    mean = round(mean(Cash), 1)
  )
##   mean
## 1 87.7
# Naive model
ATM3_naive_fit <- ATM3 %>%
  filter(Cash > 0) %>%
  model(NAIVE(Cash))
# Forecast 1 month into future
ATM3_naive_forecast <- ATM3_naive_fit %>%
  forecast(h = '1 month')
# Visualize forecast
ATM3_naive_forecast %>%
  ggtime::autoplot(ATM3) +
    # Only plot subset of data for visibility
    coord_cartesian(xlim = as.Date(c('2010-04-20', '2010-05-30'))) +  
    labs(
      y = 'Amount (USD, x100)',
      title = 'ATM3 cash withdrawals, Apr 28-30 + Forecast',
      subtitle = 'NAIVE model'
    ) +
  theme_classic() +
  theme(
    plot.title = element_text(face = 'bold'),
    axis.title = element_text(face = 'bold')
  )

Model selection

With only three observations, neither model would likely be very reliable for long-term forecasts. However, for short-term forecasts, I would use the mean forecast.

Save forecast

ATM3_forecast <- ATM3_mean_forecast %>%
  as.data.frame() %>%
  select(Date, Cash_forecast = `.mean`) %>%
  # Round forecasts since ATMs don't dispense coins
  mutate(
    Cash_forecast = floor(Cash_forecast)
  )

# Workaround because Excel doesn't read date field correctly
ATM3_forecast$Date <- as.character(ATM3_forecast$Date)

ATM4

Data checks

Missing data

There is no missing data in this time series.

(ATM4_missing <- filter(ATM4, is.na(Cash)))
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Outliers

Cash withdrawals from ATM4 had one outlier that exceeded 3 interquartile ranges.

# Boundaries for outlier detection
lower_fence <- quantile(ATM4$Cash, 0.25, na.rm = TRUE) - 3 * IQR(ATM4$Cash, na.rm = TRUE)
upper_fence <- quantile(ATM4$Cash, 0.75, na.rm = TRUE) + 3 * IQR(ATM4$Cash, na.rm = TRUE)

ATM4_outliers <- ATM4 %>%
  filter(Cash < lower_fence | Cash > upper_fence)

ATM4_outliers
## # A tsibble: 1 x 2 [1D]
##   Date        Cash
##   <date>     <dbl>
## 1 2010-02-09 10919

Negative or zero values

All Cash values are positive, which satisfies the requirement for Box-Cox transformation if it is needed.

filter(ATM4, Cash <= 0)
## # A tsibble: 0 x 2 [?]
## # ℹ 2 variables: Date <date>, Cash <dbl>

Visualize time series

ATM4_plot <- ggtime::autoplot(ATM4, Cash) + 
  # Highlight outlier boundaries
  geom_hline(yintercept = lower_fence, linetype = 'dashed', color = 'goldenrod') +
  geom_hline(yintercept = upper_fence, linetype = 'dashed', color = 'goldenrod') +  
  labs(
    y = 'Amount (USD, x100)',
    title = 'ATM4 cash withdrawals, May 2009 to Apr 2010',
    subtitle = 'Raw data has no missing observations and 1 outlier'
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(face = 'bold'),
    plot.title = element_text(face = 'bold')
  )

ATM4_plot

Data interpolation

I replaced the outlier with an estimate interpolated with an ARIMA model.

ATM4_fill <- ATM4 %>%
  # Remove outlier
  anti_join(ATM4_outliers) %>%
  # Replace with NA
  fill_gaps()
## Joining with `by = join_by(Date, Cash)`
ATM4_interpolated <- ATM4_fill %>%
  # Fit ARIMA model to data containing missing values
  model(ARIMA(Cash)) %>%
  # Estimate Cash for all time points
  interpolate(ATM4_fill)

STL decomposition

STL decomposition of the interpolated time series showed that cash withdrawals from ATM4 do not have a trend, but there is seasonality with a weekly period. Together, this indicates that the time series is not stationary, and that an additive ETS or ARIMA model could be suitable. The remainder component appears to be white noise.

ATM4_decomp <- ATM4_interpolated %>%
  model(STL(Cash ~ season(window = 'periodic'), robust = TRUE)) %>%
  components()

ATM4_decomp %>%
  ggtime::autoplot() +
    # Uncomment next line to view a more limited period  
    # coord_cartesian(xlim = as.Date(c('2009-07-01', '2009-07-31'))) +  
    ggtitle('STL decomposition of interpolated ATM4 withdrawals')

Variance and autocorrelation

  • The time series plot (top) shows non-constant variance and the cash withdrawal amounts are all positive, so Box-Cox transformation could be helpful.

  • The ACF plot (lower left) appears to be somewhat sinusoidal and has three significant spikes that decrease in magnitude as lag increases, which is consistent with seasonality and non-stationarity.

  • The partial ACF plot (lower right) has several significant spikes up to lag=19, which suggests that a higher-order seasonal ARIMA model may be needed.

ATM4_interpolated %>%
  ggtime::gg_tsdisplay(Cash, plot_type = 'partial') +
    # I plotted a limited time period to make it easier to see changes in variance
    # To view the entire time series, comment out the next line
    coord_cartesian(xlim = as.Date(c('2009-10-01', '2010-01-31'))) +
    ggtitle('ATM4 withdrawals (interpolated)')

Box-Cox transformation

The Guerrero method indicated that a Box-Cox transformation with \(\lambda = 0.45\) would help stabilize the variance, which I confirmed visually.

# Determine optimal value of lambda for Box-Cox transformation
ATM4_lambda <- ATM4_interpolated %>%
  features(Cash, features = guerrero) %>%
  pull(lambda_guerrero)

ATM4_lambda
## [1] 0.4504502

Differencing

The KPSS unit root test indicated that no seasonal differencing of the transformed data is needed to achieve stationarity. I found this result surprising given the clear seasonality in the STL decomposition. Because of this discrepancy, I tested seasonal ARIMA models with or without differencing.

ATM4_interpolated %>%
  features(box_cox(Cash, ATM4_lambda), unitroot_nsdiffs) %>%
  pull(nsdiffs)
## [1] 0

Models

Development

The data characteristics described above suggest that ETS and seasonal ARIMA models can be fit to the transformed data. I compared the performance of these models with the performance of a seasonal naive model as baseline.

# Create training and test datasets to develop models and assess performance
ATM4_train <- ATM4_interpolated %>%
  filter(year(Date) == 2009)

ATM4_test <- ATM4_interpolated %>%
  filter(year(Date) == 2010)
# ETS
ATM4_train_ets_fit <- ATM4_train %>%
  # Automatically determine the best method from the specified options
  model(ETS = ETS(box_cox(Cash, ATM4_lambda) ~ 
                    error(method = c('A', 'M')) +
                    trend(method = c('A', 'Ad', 'N')) +
                    season('A')))

# ARIMA
ATM4_train_arima_fit <- ATM4_train %>%
  # Due to unclear need for differencing, I let D = 0 or 1;
  # PDQ parameters are determined automatically
  model(ARIMA = ARIMA(box_cox(Cash, ATM4_lambda) ~ PDQ(D = 0:1)))   

# Baseline seasonal naive model
ATM4_train_snaive_fit = ATM4_train %>%
  model(SNAIVE = SNAIVE(box_cox(Cash, ATM4_lambda)))
model_list <- list(ATM4_train_ets_fit, 
                   ATM4_train_arima_fit,
                   ATM4_train_snaive_fit)

compare_forecast_accuracy(model_list, ATM4_test) %>%
  arrange(RMSE)
## # A tibble: 3 × 5
##   model   RMSE   MAE    MPE  MAPE
##   <chr>  <dbl> <dbl>  <dbl> <dbl>
## 1 ARIMA   328.  269.  -613.  643.
## 2 ETS     347.  284.  -617.  652.
## 3 SNAIVE 1537. 1277. -2902. 2907.

Comparison

The \(ARIMA(2,0,0)(2,0,0)_7 + mean\) model has the lowest RMSE, MAE, MPE, and MAPE. Due to the different model classes, comparisons using AIC or BIC are not valid.

Details of each model are shown below.

report(ATM4_train_arima_fit)
## Series: Cash 
## Model: ARIMA(2,0,0)(2,0,0)[7] w/ mean 
## Transformation: box_cox(Cash, ATM4_lambda) 
## 
## Coefficients:
##          ar1     ar2    sar1   sar2  constant
##       0.0886  0.0335  0.2020  0.177   15.6554
## s.e.  0.0650  0.0642  0.0638  0.065    0.8530
## 
## sigma^2 estimated as 191.5:  log likelihood=-989.3
## AIC=1990.61   AICc=1990.96   BIC=2011.61
report(ATM4_train_ets_fit)
## Series: Cash 
## Model: ETS(A,N,A) 
## Transformation: box_cox(Cash, ATM4_lambda) 
##   Smoothing parameters:
##     alpha = 0.0001000831 
##     gamma = 0.0001026479 
## 
##   Initial states:
##      l[0]      s[0]      s[-1]    s[-2]    s[-3]    s[-4]    s[-5]    s[-6]
##  28.79997 -17.78852 -0.3709192 1.746823 3.157084 4.371296 3.884638 4.999602
## 
##   sigma^2:  156.3077
## 
##      AIC     AICc      BIC 
## 2596.336 2597.276 2631.349
report(ATM4_train_snaive_fit)
## Series: Cash 
## Model: SNAIVE 
## Transformation: box_cox(Cash, ATM4_lambda) 
## 
## sigma^2: 313.0914

Residual diagnostics

The residual plot of the \(ARIMA(2,0,0)(2,0,0)_7 + mean\) model does not appear to have a pattern. All of the ACF spikes are non-significant. The distribution of residuals is not quite normal and has a mean greater than zero. However, the Ljung-Box test P-value is not significant, so the null hypothesis that the residuals are white noise is not rejected. Collectively, these results suggest that this model is a good fit to the data.

ggtime::gg_tsresiduals(ATM4_train_arima_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a weekly period, so l = 2m = 2*7 = 14.
# Degrees of freedom = number of parameters in model, so dof = 2
ATM4_train_arima_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 14, dof = 2)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     12.1     0.436

Forecast future values

On the basis of the above results, I selected the seasonal \(ARIMA(2,0,0)(2,0,0)_7 + mean\) model as the final model. To forecast future values, I refit this model to the entire dataset (ie, train + test), and then used it to forecast cash withdrawals 1 month into the future.

# Refit the training model to the full dataset
ATM4_final_arima_fit <- ATM4_train_arima_fit %>%
  fabletools::refit(ATM4_interpolated)
# Forecast future values
ATM4_forecast <- ATM4_final_arima_fit %>%
  forecast(h = '1 month')

The forecast looks reasonable.

# Plot forecast
ATM4_forecast %>%
  ggtime::autoplot(ATM4_interpolated) +
    labs(
      y = 'Amount (USD, x100)',
      title = 'ATM4 cash withdrawals, May 2009 to Apr 2010 + Forecast',
      subtitle = latex2exp::TeX('ARIMA(2,0,0)(2,0,0)$_{7}$ + mean model')    
    ) +
    theme_classic() +
    theme(
      axis.title = element_text(face = 'bold'),
      plot.title = element_text(face = 'bold')
    )

Save forecast

ATM4_forecast <- ATM4_forecast %>%
  as.data.frame() %>%
  select(Date, Cash_forecast = `.mean`) %>%
  # Round forecasts since ATMs don't dispense coins
  mutate(
    Cash_forecast = floor(Cash_forecast)
  )

# Workaround because Excel doesn't read date field correctly
ATM4_forecast$Date <- as.character(ATM4_forecast$Date)

B. Power forecast

Data input and wrangling

filename <- 'ResidentialCustomerForecastLoad-624.xlsx'
url <- paste0(domain, path, filename)

power_df <- read.xlsx(url, sheet = 1)

The raw dataset is a tibble (dataframe) with 3 variables and 192 observations. For time series analysis, the dataframe needs to be converted to a tsibble, which in turn, requires the time variable to be converted a year-month data type and defined as the index variable.

class(power_df)
## [1] "data.frame"
glimpse(power_df)
## Rows: 192
## Columns: 3
## $ CaseSequence <dbl> 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 74…
## $ `YYYY-MMM`   <chr> "1998-Jan", "1998-Feb", "1998-Mar", "1998-Apr", "1998-May…
## $ KWH          <dbl> 6862583, 5838198, 5420658, 5010364, 4665377, 6467147, 891…
power_ts <- power_df %>%
  # Rename date column to be more meaningful
  rename('Date' = 'YYYY-MMM') %>%
  # Coerce date as year-month data type
  mutate(
    Date = yearmonth(Date)
  ) %>%
  # convert to tsibble
  as_tsibble(index = Date)

Data checks

Missing data

There is one missing observation for September 2008.

(power_missing <- filter(power_ts, is.na(KWH)))
## # A tsibble: 1 x 3 [1M]
##   CaseSequence     Date   KWH
##          <dbl>    <mth> <dbl>
## 1          861 2008 Sep    NA

Outliers

There were no outliers that exceeded 3 interquartile ranges; however, one exceeded 1.5 IQRs. Based on the time series plot below, which shows that this data point is clearly lower than other data points, I decided to use the less stringent outlier criteria and conclude that there is one outlier.

# Boundaries for outlier detection
lower_fence <- quantile(power_ts$KWH, 0.25, na.rm = TRUE) - 
                        3 * IQR(power_ts$KWH, na.rm = TRUE)
upper_fence <- quantile(power_ts$KWH, 0.75, na.rm = TRUE) + 
                        3 * IQR(power_ts$KWH, na.rm = TRUE)

power_outliers <- power_ts %>%
  filter(KWH < lower_fence | KWH > upper_fence)

power_outliers
## # A tsibble: 0 x 3 [?]
## # ℹ 3 variables: CaseSequence <dbl>, Date <mth>, KWH <dbl>
# Less stringent boundaries for outlier detection
lower_fence <- quantile(power_ts$KWH, 0.25, na.rm = TRUE) - 
                        1.5 * IQR(power_ts$KWH, na.rm = TRUE)
upper_fence <- quantile(power_ts$KWH, 0.75, na.rm = TRUE) + 
                        1.5 * IQR(power_ts$KWH, na.rm = TRUE)

power_outliers <- power_ts %>%
  filter(KWH < lower_fence | KWH > upper_fence)

power_outliers
## # A tsibble: 1 x 3 [1M]
##   CaseSequence     Date    KWH
##          <dbl>    <mth>  <dbl>
## 1          883 2010 Jul 770523

Negative or zero values

All KWH values were positive, which satisfies requirement for Box-Cox transformation if it is needed.

filter(power_ts, KWH <= 0)
## # A tsibble: 0 x 3 [?]
## # ℹ 3 variables: CaseSequence <dbl>, Date <mth>, KWH <dbl>

Visualize time series

The overall time series plot shows the presence of seasonality (monthly period), one missing observation, and one outlier.

power_plot <- ggtime::autoplot(power_ts, KWH) + 
  # Highlight outlier boundaries
  geom_hline(yintercept = lower_fence, linetype = 'dashed', color = 'goldenrod') +
  geom_hline(yintercept = upper_fence, linetype = 'dashed', color = 'goldenrod') +  
  labs(
    y = 'Power usage (kWH)',
    title = 'Residential power usage, 1998-2013',
    subtitle = 'Raw data has one missing observation and one outlier'
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(face = 'bold'),
    plot.title = element_text(face = 'bold'),
  )

# Add vertical lines to highlight missing data points
for (i in 1:nrow(power_missing)) {
  power_plot <- power_plot + 
    geom_vline(xintercept = power_missing[[i, 'Date']], 
               linetype = 'dashed', color = 'steelblue')
}

power_plot

Data interpolation

I replaced the missing observation and outlier with estimates interpolated with an ARIMA model.

power_fill <- power_ts %>%
  # Remove outlier
  anti_join(power_outliers) %>%
  # Replace with NA
  fill_gaps()
## Joining with `by = join_by(CaseSequence, Date, KWH)`
power_interpolated <- power_fill %>%
  # Fit ARIMA model to data containing missing values
  model(ARIMA(KWH)) %>%
  # Estimate KWH for all time points
  interpolate(power_fill)

STL decomposition

STL decomposition of the interpolated time series showed that power usage has a upward trend and seasonality, which indicate that the time series is not stationary. A seasonal plot further revealed that power usage peaked during winter and summer (ie, 6-month period). The remainder component appears to be white noise.

Of note, the trend component is relatively flat before 2009 and then subsequently trends upward. This suggests that either an additive or multiplicative ETS model could be suitable. In addition, the trend flattens out at the end of the time series, which may indicate a need for dampening in an ETS model.

power_interpolated %>%
  # Uncomment next line to view seasonality for a more limited period
  # filter(year(Date) >= 2005 & year(Date) <= 2010) %>%
  model(STL(KWH ~ season(window = 'periodic'), robust = TRUE)) %>%
  components() %>%
  ggtime::autoplot() +
    ggtitle('STL decomposition of interpolated power usage data')

ggtime::gg_season(power_interpolated, KWH) +
  labs(
    x = 'Month',
    y = 'Power usage (kWH)',
    title = 'Seasonality of residential power usage',
    subtitle = 'Interpolated data'
  ) +
  theme_classic() +
  theme(
    axis.title = element_text(face = 'bold'),
    plot.title = element_text(face = 'bold')
  )

Variance and autocorrelation

  • The time series plot (top) shows some evidence of non-constant variance and the power usage data are all positive, so Box-Cox transformation could be helpful.

  • The ACF plot (lower left) has a sinusoidal pattern of significant spikes, which is consistent with seasonality and non-stationarity.

  • The partial ACF plot (lower right) shows an exponential decay of spikes after lag=12,, which suggests that a higher-order seasonal ARIMA model may be needed.

power_interpolated %>%
  ggtime::gg_tsdisplay(KWH, plot_type = 'partial') +
  ggtitle('Residential power usage')

Box-Cox transformation

The Guerrero method indicated that a Box-Cox transformation with \(\lambda = -0.21\) would help stabilize the variance, which I confirmed visually.

# Determine optimal value of lambda for Box-Cox transformation
power_lambda <- power_interpolated %>%
  features(KWH, features = guerrero) %>%
  pull(lambda_guerrero)

power_lambda
## [1] -0.2057366

Differencing

The KPSS unit root test indicated that first-order seasonal differencing of the transformed data is sufficient to achieve stationarity.

power_interpolated %>%
  features(box_cox(KWH, power_lambda), unitroot_nsdiffs) %>%
  pull(nsdiffs)
## [1] 1

Models

Approach 1: Transformed data

Development

The data characteristics described above suggest that ETS and seasonal ARIMA models can be fit to the transformed data. I compared the performance of these models with the performance of a seasonal naive model as baseline.

# Create training and test datasets to develop models and assess performance
power_train <- power_interpolated %>%
  filter(year(Date) < 2011)

power_test <- power_interpolated %>%
  filter(year(Date) == 2010)
# ETS
power_train_ets_fit <- power_train %>%
  # Automatically determine the best method from the specified options
  model(ETS = ETS(box_cox(KWH, power_lambda) ~ 
                    error(method = c('A', 'M')) +
                    trend(method = c('A', 'Ad', 'M')) +
                    season(method = c('A', 'M'))))

# ARIMA
power_train_arima_fit <- power_train %>%
  model(ARIMA = ARIMA(box_cox(KWH, power_lambda) ~ PDQ(D = 1)))   

# Baseline seasonal naive model
power_train_snaive_fit = power_train %>%
  model(SNAIVE = SNAIVE(box_cox(KWH, power_lambda)))
# This block takes a few seconds to run
model_list <- list(power_train_ets_fit, 
                   power_train_arima_fit,
                   power_train_snaive_fit)

compare_forecast_accuracy(model_list, power_test) %>%
  arrange(RMSE)
## # A tibble: 3 × 5
##   model     RMSE     MAE    MPE  MAPE
##   <chr>    <dbl>   <dbl>  <dbl> <dbl>
## 1 SNAIVE  58042.  55941. -0.797 0.797
## 2 ARIMA  487284. 381902.  1.35  5.30 
## 3 ETS    652414. 457557.  3.45  6.01

Comparison

The seasonal naive model had the lowest RMSE, MAE, MPE, and MAPE. Due to the different model classes, comparisons using AIC or BIC are not valid.

Details of each model are shown below.

report(power_train_arima_fit)
## Series: KWH 
## Model: ARIMA(0,0,1)(2,1,0)[12] w/ drift 
## Transformation: box_cox(KWH, power_lambda) 
## 
## Coefficients:
##          ma1     sar1     sar2  constant
##       0.1699  -0.7387  -0.3878     8e-04
## s.e.  0.0962   0.0842   0.0840     4e-04
## 
## sigma^2 estimated as 1.398e-05:  log likelihood=608.86
## AIC=-1207.73   AICc=-1207.29   BIC=-1192.88
report(power_train_ets_fit)
## Series: KWH 
## Model: ETS(M,M,M) 
## Transformation: box_cox(KWH, power_lambda) 
##   Smoothing parameters:
##     alpha = 0.0001461708 
##     beta  = 0.0001251217 
##     gamma = 0.0003705203 
## 
##   Initial states:
##      l[0]     b[0]      s[0]     s[-1]     s[-2]   s[-3]    s[-4]    s[-5]
##  4.663809 1.000006 0.9994214 0.9973231 0.9990265 1.00168 1.002111 1.001629
##   s[-6]     s[-7]     s[-8]     s[-9]   s[-10]   s[-11]
##  0.9998 0.9981181 0.9986163 0.9995482 1.000767 1.001961
## 
##   sigma^2:  0
## 
##       AIC      AICc       BIC 
## -960.0912 -955.6564 -908.2436
report(power_train_snaive_fit)
## Series: KWH 
## Model: SNAIVE 
## Transformation: box_cox(KWH, power_lambda) 
## 
## sigma^2: 0

Residual diagnostics

As explained below, none of the three models had residuals consistent with white noise, which indicated that they are not a good fit to the data.

The residual plot of the seasonal naive model does not appear to have a pattern, and the distribution of residuals appears to be close to normal with a mean near zero. However, the ACF plot has a sinusoidal pattern and there is a significant spike at lag=12, which is consistent with the presence of residual seasonality. In addition, the Ljung-Box test P-value was significant, so the null hypothesis that the residuals are white noise was rejected. Collectively, these results suggest that this model is not a good fit to the data.

ggtime::gg_tsresiduals(power_train_snaive_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 0
power_train_snaive_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 SNAIVE    68.9  5.08e-10

The residual plot of the \(ARIMA(0,0,1)(2,1,0)_{12} + drift\) model does not appear to have a pattern. All but one of the ACF spikes are non-significant. The distribution of residuals is not quite normal and has a mean greater than zero. In addition, the Ljung-Box test P-value was significant, so the null hypothesis that the residuals are white noise was rejected. Collectively, these results suggest that this model is not a good fit to the data.

ggtime::gg_tsresiduals(power_train_arima_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 2
power_train_arima_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 2)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     30.1  0.000825

The residual plot of the ETS model does not appear to have a pattern. However, the ACF plot shows a decaying sinusoidal pattern with several significant spikes, which suggests residual seasonality. The distribution of residuals appears to be close to normal, although the mean appears to be less than zero. The Ljung-Box test P-value was significant, so the null hypothesis that the residuals are white noise was rejected. Collectively, these results suggest that this model is not a good fit to the data.

ggtime::gg_tsresiduals(power_train_ets_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 6 + (m-1) = 11
power_train_ets_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 11)
## # A tibble: 1 × 3
##   .model lb_stat     lb_pvalue
##   <chr>    <dbl>         <dbl>
## 1 ETS       36.0 0.00000000197

Approach 2: Un-transformed data

Next, I decided to try modeling without Box-Cox transformation. My rationale for this was to undo changes to the original data to determine whether this improves the fit.

Differencing

The KPSS unit root test indicated that first-order seasonal differencing of the un-transformed data is sufficient to achieve stationarity.

power_interpolated %>%
  features(KWH, unitroot_nsdiffs) %>%
  pull(nsdiffs)
## [1] 1

Development

# Create training and test datasets to develop models and assess performance
power_train <- power_interpolated %>%
  filter(year(Date) < 2011)

power_test <- power_interpolated %>%
  filter(year(Date) == 2010)
# ETS
power_train_ets_fit <- power_train %>%
  # Automatically determine the best method from the specified options
  model(ETS = ETS(KWH ~ error(method = c('A', 'M')) +
                              trend(method = c('A', 'Ad', 'M')) +
                              season(method = c('A', 'M'))))

# ARIMA
power_train_arima_fit <- power_train %>%
  model(ARIMA = ARIMA(KWH ~ PDQ(D = 1)))   

# Baseline seasonal naive model
power_train_snaive_fit = power_train %>%
  model(SNAIVE = SNAIVE(KWH))
# This block takes a few seconds to run
model_list <- list(power_train_ets_fit, 
                   power_train_arima_fit,
                   power_train_snaive_fit)

compare_forecast_accuracy(model_list, power_test) %>%
  arrange(RMSE)
## # A tibble: 3 × 5
##   model     RMSE     MAE   MPE  MAPE
##   <chr>    <dbl>   <dbl> <dbl> <dbl>
## 1 SNAIVE      0       0   0     0   
## 2 ARIMA  474883. 366900.  1.85  5.07
## 3 ETS    689525. 474174.  3.62  6.22

Comparison

I was not sure why all the accuracy metrics for the seasonal naive model are zero. I will examine its residuals below before jumping to the conclusion that this model is a perfect fit.

Compared with the \(ETS(M,M,M)\) model, the \(ARIMA(1,0,1)(2,1,0)_{12} + drift\) model has lower RMSE, MAE, MPE, and MAPE.

Due to the different model classes, comparisons using AIC or BIC are not valid. Details of each model are shown below.

report(power_train_arima_fit)
## Series: KWH 
## Model: ARIMA(1,0,1)(2,1,0)[12] w/ drift 
## 
## Coefficients:
##           ar1     ma1     sar1     sar2   constant
##       -0.4992  0.7717  -0.7303  -0.4233  187101.25
## s.e.   0.1332  0.0923   0.0845   0.0829   83347.93
## 
## sigma^2 estimated as 2.901e+11:  log likelihood=-2106.45
## AIC=4224.89   AICc=4225.51   BIC=4242.71
report(power_train_ets_fit)
## Series: KWH 
## Model: ETS(M,M,M) 
##   Smoothing parameters:
##     alpha = 0.06281518 
##     beta  = 0.0001000195 
##     gamma = 0.0001017283 
## 
##   Initial states:
##     l[0]     b[0]      s[0]     s[-1]     s[-2]    s[-3]   s[-4]    s[-5]
##  6218746 1.000604 0.9190375 0.7458877 0.8983215 1.200753 1.24708 1.205941
##      s[-6]     s[-7]     s[-8]     s[-9]   s[-10]   s[-11]
##  0.9757765 0.7731204 0.8219729 0.9322194 1.058363 1.221527
## 
##   sigma^2:  0.0077
## 
##      AIC     AICc      BIC 
## 4926.469 4930.904 4978.316
report(power_train_snaive_fit)
## Series: KWH 
## Model: SNAIVE 
## 
## sigma^2: 512316892917.58

Residual diagnostics

The residual plot of the seasonal naive model does not appear to have a pattern. Several spikes in the ACF plot exceed the critical values. The distribution of residuals is not normal and the mean is less than zero. In addition, the Ljung-Box test P-value was significant, so the null hypothesis that the residuals are white noise was rejected. Collectively, these results suggest that this model is not a good fit to the data.

ggtime::gg_tsresiduals(power_train_snaive_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 0
power_train_snaive_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 0)
## # A tibble: 1 × 3
##   .model lb_stat     lb_pvalue
##   <chr>    <dbl>         <dbl>
## 1 SNAIVE    62.1 0.00000000939

The residual plot of the \(ARIMA(1,0,1)(2,1,0)_{12} + drift\) model does not appear to have a pattern. All but one of the ACF spikes are non-significant (and the other only slightly exceeds the upper critical value). The distribution of residuals is close to normal, although the mean is greater than zero. The Ljung-Box test P-value was not significant, so the null hypothesis that the residuals are white noise was not rejected. Collectively, these results suggest that this model is a good fit to the data.

ggtime::gg_tsresiduals(power_train_arima_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 3
power_train_arima_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 3)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     14.1     0.118

The residual plot of the ETS model does not appear to have a pattern. However, the ACF plot has several significant spikes. The distribution of residuals is not quite normal, and the mean is less than zero. The Ljung-Box test P-value was significant, so the null hypothesis that the residuals are white noise was rejected. Collectively, these results suggest that this model is not a good fit to the data.

ggtime::gg_tsresiduals(power_train_ets_fit)

# Ljung-Box test of residuals
# The time series is seasonal with a semi-annual period, so l = 2m = 2*6 = 12
# Degrees of freedom = number of parameters in model, so dof = 6 + (m-1) = 11
power_train_ets_fit %>%
  residuals() %>%
  features(.resid, ljung_box, lag = 12, dof = 11)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ETS       37.7  8.04e-10

Forecast future values

On the basis of the above results, I selected the seasonal \(ARIMA(1,0,1)(2,1,0)_{12} + drift\) model of the un-transformed data as the final model. To forecast future values, I refit this model to the entire dataset (ie, train + test), and then used it to forecast power usage 1 year into the future.

# Refit the training model to the full dataset
power_final_arima_fit <- power_train_arima_fit %>%
  fabletools::refit(power_interpolated)
# Forecast future values
power_forecast <- power_final_arima_fit %>%
  forecast(h = '1 year')

The forecast looks reasonable.

# Plot forecast
power_forecast %>%
  ggtime::autoplot(power_interpolated) +
    labs(
      y = 'Power usage (kWH)',
      title = 'Residential power usage, 1998-2013 + Forecast',
      subtitle = latex2exp::TeX('ARIMA(1,0,1)(2,1,0)$_{12}$ + drift model')       
    ) +
    theme_classic() +
    theme(
      axis.title = element_text(face = 'bold'),
      plot.title = element_text(face = 'bold')
    )

Save forecast

power_forecast <- power_forecast %>%
  as.data.frame() %>%
  select(Date, KWH_forecast = `.mean`)

# Workaround because Excel doesn't read date field correctly
power_forecast$Date <- as.character(power_forecast$Date)

Export all forecasts

# List of forecast dataframes
forecasts_list <- list('ATM1' = ATM1_forecast, 
                    'ATM2' = ATM2_forecast, 
                    'ATM3' = ATM3_forecast, 
                    'ATM4' = ATM4_forecast,
                    'power' = power_forecast)

# Save to a single Excel file (1 forecast per sheet)
write.xlsx(forecasts_list, file = 'project1_forecasts.xlsx')

Session Details

pander::pander(sessionInfo())

R version 4.5.2 (2025-10-31)

Platform: aarch64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pander(v.0.6.6), latex2exp(v.0.9.8), forecast(v.9.0.0), fable(v.0.5.0), feasts(v.0.4.2), fabletools(v.0.6.1), ggtime(v.0.2.0), tsibble(v.1.2.0), openxlsx(v.4.2.8.1), lubridate(v.1.9.5), forcats(v.1.0.1), stringr(v.1.6.0), dplyr(v.1.2.0), purrr(v.1.2.1), readr(v.2.1.6), tidyr(v.1.3.2), tibble(v.3.3.1), ggplot2(v.4.0.2) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): gtable(v.0.3.6), anytime(v.0.3.12), xfun(v.0.56), bslib(v.0.10.0), lattice(v.0.22-7), numDeriv(v.2016.8-1.1), tzdb(v.0.5.0), quadprog(v.1.5-8), vctrs(v.0.7.2), tools(v.4.5.2), generics(v.0.1.4), curl(v.7.0.0), parallel(v.4.5.2), xts(v.0.14.1), pkgconfig(v.2.0.3), RColorBrewer(v.1.1-3), S7(v.0.2.1), distributional(v.0.7.0), lifecycle(v.1.0.5), compiler(v.4.5.2), farver(v.2.1.2), htmltools(v.0.5.9), sass(v.0.4.10), yaml(v.2.3.12), pillar(v.1.11.1), jquerylib(v.0.1.4), cachem(v.1.1.0), nlme(v.3.1-168), fracdiff(v.1.5-3), tidyselect(v.1.2.1), zip(v.2.3.3), digest(v.0.6.39), stringi(v.1.8.7), labeling(v.0.4.3), tseries(v.0.10-60), fastmap(v.1.2.0), grid(v.4.5.2), colorspace(v.2.1-2), cli(v.3.6.5), magrittr(v.2.0.4), utf8(v.1.2.6), withr(v.3.0.2), scales(v.1.4.0), timechange(v.0.4.0), TTR(v.0.24.4), rmarkdown(v.2.30), quantmod(v.0.4.28), otel(v.0.2.0), nnet(v.7.3-20), timeDate(v.4052.112), progressr(v.0.18.0), zoo(v.1.8-15), hms(v.1.1.4), urca(v.1.3-4), evaluate(v.1.0.5), knitr(v.1.51), ggdist(v.3.3.3), lmtest(v.0.9-40), rlang(v.1.1.7), Rcpp(v.1.1.1), glue(v.1.8.0), rstudioapi(v.0.18.0), jsonlite(v.2.0.0) and R6(v.2.6.1)

  1. Hyndman Chapter 9, section 9.10.↩︎