DATA 624 Project 1

Part A: ATM cash withdrawals | Part B: residential power usage | Part C: water flow (bonus)

knitr::opts_chunk$set(message = FALSE, warning = FALSE)
options(cli.unicode = FALSE)
library(fpp3)
library(readxl)
library(writexl)
library(lubridate)

Introduction

This report covers all three parts of the project.

• Part A forecasts daily cash withdrawals at four ATMs for May 2010.
• Part B forecasts monthly residential power usage for 2014.
• Part C (bonus) aggregates two water flow data sets to hourly means, tests whether they are stationary and forecastable, and produces a one week forecast.

The same approach is used throughout: inspect and clean the data, look at the patterns, compare several candidate models on a hold-out test set, check the residuals of the chosen model, then refit on all the data to produce the final forecast. All forecasts are written to a single Excel file at the end.

Part A: ATM cash withdrawals

Reading and inspecting the data

atm_raw <- read_excel("ATM624Data.xlsx")
str(atm_raw)

## tibble [1,474 x 3] (S3: tbl_df/tbl/data.frame)
##  $ DATE: num [1:1474] 39934 39934 39935 39935 39936 ...
##  $ ATM : chr [1:1474] "ATM1" "ATM2" "ATM1" "ATM2" ...
##  $ Cash: num [1:1474] 96 107 82 89 85 90 90 55 99 79 ...

head(atm_raw)

## # A tibble: 6 x 3
##    DATE ATM    Cash
##   <dbl> <chr> <dbl>
## 1 39934 ATM1     96
## 2 39934 ATM2    107
## 3 39935 ATM1     82
## 4 39935 ATM2     89
## 5 39936 ATM1     85
## 6 39936 ATM2     90

The DATE column arrives as an Excel serial number, so it needs converting. Cash is in hundreds of dollars.

atm_raw <- atm_raw |> mutate(DATE = as.Date(DATE, origin = "1899-12-30"))

range(atm_raw$DATE)

## [1] "2009-05-01" "2010-05-14"

table(atm_raw$ATM, useNA = "ifany")

## 
## ATM1 ATM2 ATM3 ATM4 <NA> 
##  365  365  365  365   14

The file holds one year of daily data for each of four ATMs (1 May 2009 to 30 April 2010), plus 14 extra rows dated May 2010 where both ATM and Cash are missing. Those rows are placeholders for the period we are asked to forecast, so they are removed.

atm_raw |> filter(is.na(ATM)) |> head(3)

## # A tibble: 3 x 3
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2010-05-01 <NA>     NA
## 2 2010-05-02 <NA>     NA
## 3 2010-05-03 <NA>     NA

atm <- atm_raw |> filter(!is.na(ATM))

Data quality issues

atm |> filter(is.na(Cash))

## # A tibble: 5 x 3
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2009-06-13 ATM1     NA
## 2 2009-06-16 ATM1     NA
## 3 2009-06-18 ATM2     NA
## 4 2009-06-22 ATM1     NA
## 5 2009-06-24 ATM2     NA

Five days have a missing Cash value: three at ATM1 and two at ATM2, all in June 2009.

atm |>
  group_by(ATM) |>
  summarise(n = n(), missing = sum(is.na(Cash)), zeros = sum(Cash == 0, na.rm = TRUE),
            min = min(Cash, na.rm = TRUE), median = median(Cash, na.rm = TRUE),
            mean = round(mean(Cash, na.rm = TRUE), 1), max = max(Cash, na.rm = TRUE))

## # A tibble: 4 x 8
##   ATM       n missing zeros   min median  mean    max
##   <chr> <int>   <int> <int> <dbl>  <dbl> <dbl>  <dbl>
## 1 ATM1    365       3     0  1       91   83.9   180 
## 2 ATM2    365       2     2  0       67   62.6   147 
## 3 ATM3    365       0   362  0        0    0.7    96 
## 4 ATM4    365       0     0  1.56   404. 474   10920.

Two further problems stand out:

• ATM3 is almost entirely zeros. Only the last three days of the series have any activity, which suggests the machine was installed (or brought online) at the very end of April 2010.
• ATM4 has an enormous maximum of 10,920 against a median of about 404.

atm |> filter(ATM == "ATM3", Cash > 0)

## # A tibble: 3 x 3
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2010-04-28 ATM3     96
## 2 2010-04-29 ATM3     82
## 3 2010-04-30 ATM3     85

atm |> filter(ATM == "ATM4") |> slice_max(Cash, n = 4) |> select(DATE, Cash)

## # A tibble: 4 x 2
##   DATE         Cash
##   <date>      <dbl>
## 1 2010-02-09 10920.
## 2 2009-09-22  1712.
## 3 2010-01-29  1575.
## 4 2009-09-04  1495.

The ATM4 value of $10,920 hundred (about $1.09 million in one day) on 9 February 2010 is more than six times the next largest value in the series. No ATM holds that much cash, so this is a data error rather than a real withdrawal. I treat it as missing.

atm |> ggplot(aes(DATE, Cash)) +
  geom_line() +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  labs(title = "Daily cash withdrawals by ATM (raw data)",
       y = "Cash (hundreds of dollars)", x = NULL)

Cleaning

Missing values and the single outlier are replaced with the median for that ATM on that day of the week. Using the weekday median rather than an overall median preserves the weekly pattern described below, and with only six values involved out of 1,460 the choice has very little influence on the forecasts.

atm_clean <- atm |>
  mutate(Cash = if_else(ATM == "ATM4" & Cash > 5000, NA_real_, Cash)) |>
  mutate(weekday = wday(DATE)) |>
  group_by(ATM, weekday) |>
  mutate(Cash = if_else(is.na(Cash), median(Cash, na.rm = TRUE), Cash)) |>
  ungroup() |>
  select(-weekday)

atm_ts <- atm_clean |> as_tsibble(index = DATE, key = ATM)

sum(is.na(atm_ts$Cash))

## [1] 0

atm_ts |> autoplot(Cash) +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  labs(title = "Daily cash withdrawals by ATM (cleaned)",
       y = "Cash (hundreds of dollars)", x = NULL) +
  guides(colour = "none")

Exploring the patterns

atm_clean |>
  mutate(weekday = wday(DATE, label = TRUE)) |>
  group_by(ATM, weekday) |>
  summarise(mean_cash = round(mean(Cash), 1), .groups = "drop") |>
  tidyr::pivot_wider(names_from = weekday, values_from = mean_cash)

## # A tibble: 4 x 8
##   ATM     Sun   Mon   Tue   Wed   Thu   Fri   Sat
##   <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ATM1  103.   86    89.8  82.2  31.7  98.6  96.6
## 2 ATM2   67.2  58.7  73.2  43.4  25.2  92    76  
## 3 ATM3    0     0     0     1.8   1.6   1.6   0  
## 4 ATM4  539.  481.  443.  414.  170   574.  492.

Averaged over the whole year, Thursday looks like the quiet day at every machine. That turns out to be misleading, so it is worth splitting the year into periods:

atm_clean |>
  filter(ATM != "ATM3") |>
  mutate(weekday = wday(DATE, label = TRUE),
         period = case_when(DATE < as.Date("2009-11-01") ~ "1. May-Oct 2009",
                            DATE < as.Date("2010-02-01") ~ "2. Nov 2009-Jan 2010",
                            TRUE ~ "3. Feb-Apr 2010")) |>
  group_by(ATM, period, weekday) |>
  summarise(mean_cash = round(mean(Cash), 1), .groups = "drop") |>
  tidyr::pivot_wider(names_from = weekday, values_from = mean_cash) |>
  arrange(ATM, period)

## # A tibble: 9 x 9
##   ATM   period                 Sun   Mon   Tue   Wed   Thu   Fri   Sat
##   <chr> <chr>                <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ATM1  1. May-Oct 2009      102.   94.8 114.   82.3  14.3  99.2 100. 
## 2 ATM1  2. Nov 2009-Jan 2010 103.   82.2  99.5  62.8  15.6 102.   95  
## 3 ATM1  3. Feb-Apr 2010      103.   72.2  30.8 101.   82.6  94    90.2
## 4 ATM2  1. May-Oct 2009       74.2  68.8  98.7  29.9   9.3 106.   90.3
## 5 ATM2  2. Nov 2009-Jan 2010  52.5  77.4  79.5  23.9   5.5 101.   55.7
## 6 ATM2  3. Feb-Apr 2010       69.1  19.9  16.1  90.1  76.5  54.5  65.8
## 7 ATM4  1. May-Oct 2009      552.  506.  545.  417.   57.8 578.  532. 
## 8 ATM4  2. Nov 2009-Jan 2010 627.  476.  555   298.   65.7 598   468. 
## 9 ATM4  3. Feb-Apr 2010      407.  437.  127.  522.  498.  540.  427.

atm_ts |>
  filter(ATM != "ATM3") |>
  gg_season(Cash, period = "week") +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  labs(title = "Weekly pattern by ATM", y = "Cash (hundreds of dollars)")

atm_ts |> filter(ATM == "ATM1") |> ACF(Cash, lag_max = 28) |> autoplot() +
  labs(title = "ACF of ATM1")

atm_ts |> filter(ATM == "ATM2") |> ACF(Cash, lag_max = 28) |> autoplot() +
  labs(title = "ACF of ATM2")

atm_ts |> filter(ATM == "ATM4") |> ACF(Cash, lag_max = 28) |> autoplot() +
  labs(title = "ACF of ATM4")

The main feature of ATM1, ATM2 and ATM4 is a strong weekly cycle, which the ACF confirms with clear spikes at lags 7, 14, 21 and 28. None of the series shows a meaningful long-run trend, so the level over the year is fairly stable.

The period table reveals something more important: the shape of the weekly cycle changed around February 2010, and it changed in the same way at all three machines. Through to January, Thursday was by far the quietest day (ATM2 averaged 9.4 on Thursdays against 105.9 on Fridays) and Tuesday was one of the busiest. From February onwards the roles swap: Tuesday collapses (ATM2 averages 16.1, ATM4 drops from about 550 to 112) while Thursday recovers to a normal level. Wednesday also rises at every machine.

Whether this is a genuine change in customer behaviour, a change in the servicing or replenishment schedule, or a labelling problem in the source system, it matters a great deal for a one-month forecast: the last three months are a better guide to May than the earlier part of the year. This is the single most important feature of the data, and it shapes the model choice below.

ATM4 operates at a completely different scale (averaging about 445 against roughly 84 for ATM1), and is also much more volatile from day to day.

ATM3 is a different problem entirely, and I handle it separately below.

Model selection

I hold out April 2010 (the last 30 days) as a test set and compare five candidate models on each of ATM1, ATM2 and ATM4:

• SNAIVE: the seasonal naive benchmark, repeating the same day of the previous week.
• ETS: automatic exponential smoothing.
• ARIMA: automatic ARIMA with weekly seasonality.
• ARIMA with a Box-Cox transformation, since the series are right-skewed and bounded below by zero.
• STL + ETS: an STL decomposition with ETS applied to the seasonally adjusted series.

atm_model_data <- atm_ts |> filter(ATM != "ATM3")

train <- atm_model_data |> filter(DATE < as.Date("2010-04-01"))

lambdas <- train |> features(Cash, features = guerrero)
lambdas

## # A tibble: 3 x 2
##   ATM   lambda_guerrero
##   <chr>           <dbl>
## 1 ATM1            0.430
## 2 ATM2            0.649
## 3 ATM4            0.395

fit_candidates <- train |>
  model(snaive = SNAIVE(Cash),
        ets = ETS(Cash),
        arima = ARIMA(Cash),
        arima_boxcox = ARIMA(box_cox(Cash, 0.5)),
        stl_ets = decomposition_model(STL(Cash), ETS(season_adjust ~ season("N"))))

fc_candidates <- fit_candidates |> forecast(h = 30)

accuracy(fc_candidates, atm_model_data) |>
  select(ATM, .model, RMSE, MAE, MAPE, MASE) |>
  arrange(ATM, RMSE)

## # A tibble: 15 x 6
##    ATM   .model        RMSE    MAE   MAPE  MASE
##    <chr> <chr>        <dbl>  <dbl>  <dbl> <dbl>
##  1 ATM1  stl_ets       11.6   9.11   72.2 0.494
##  2 ATM1  ets           11.6   9.12   65.1 0.494
##  3 ATM1  arima_boxcox  12.0   9.98   61.7 0.541
##  4 ATM1  arima         12.4   9.69   86.3 0.525
##  5 ATM1  snaive        16.0  13.5    73.4 0.733
##  6 ATM2  stl_ets       19.0  13.9    98.4 0.668
##  7 ATM2  ets           19.4  14.3    54.3 0.689
##  8 ATM2  arima         21.0  16.3    68.7 0.782
##  9 ATM2  arima_boxcox  21.1  15.0    61.9 0.721
## 10 ATM2  snaive        26.2  17.8    45.1 0.856
## 11 ATM4  stl_ets      249.  208.    560.  0.591
## 12 ATM4  arima_boxcox 261.  225.    880.  0.639
## 13 ATM4  arima        285.  243.   1082.  0.688
## 14 ATM4  snaive       301.  227.    375.  0.644
## 15 ATM4  ets          311.  263.   1279.  0.745

STL + ETS has the lowest test RMSE for all three ATMs, so I use it for each of them. The pattern of results is consistent and easy to explain:

• The STL decomposition handles the weekly cycle directly, and the ETS part then models the remaining level, so the method matches the structure of these series. Crucially, STL lets the seasonal shape evolve over time, so it tracks the February change described above instead of averaging the old and new patterns together. That is almost certainly why it wins here.
• ETS alone is a close second. The automatic selection picks ETS(A,N,A), that is additive weekly seasonality with no trend, which describes the data well.
• The ARIMA models are a little worse, and seasonal naive is clearly the weakest, which is expected since it throws away everything except last week’s value.
• The Box-Cox transformation does not help. It gives the lowest MAPE for ATM1 (MAPE is a poor guide here because some days have very small withdrawals, which inflates percentage errors), but its RMSE is no better.

fc_candidates |>
  autoplot(atm_model_data |> filter(DATE >= as.Date("2010-02-01")), level = NULL) +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  labs(title = "Test set forecasts for April 2010", y = "Cash (hundreds of dollars)")

Residual diagnostics

final_fit <- atm_model_data |>
  model(stl_ets = decomposition_model(STL(Cash), ETS(season_adjust ~ season("N"))))

final_fit |> filter(ATM == "ATM1") |> gg_tsresiduals() + labs(title = "ATM1 residuals")

final_fit |> filter(ATM == "ATM2") |> gg_tsresiduals() + labs(title = "ATM2 residuals")

final_fit |> filter(ATM == "ATM4") |> gg_tsresiduals() + labs(title = "ATM4 residuals")

final_fit |> augment() |> features(.innov, ljung_box, lag = 14)

## # A tibble: 3 x 4
##   ATM   .model  lb_stat lb_pvalue
##   <chr> <chr>     <dbl>     <dbl>
## 1 ATM1  stl_ets    33.1 0.00282  
## 2 ATM2  stl_ets    22.5 0.0683   
## 3 ATM4  stl_ets    47.0 0.0000191

The residuals are centred on zero and show no obvious change in variance over time. There is some remaining autocorrelation, so the models have not extracted every last pattern, but the residual ACF spikes are small and the point forecasts are still the best of the candidates tested. The ATM4 residuals have the heaviest tails, which reflects how variable that machine is; its prediction intervals should be read as approximate.

ATM3

atm_ts |> filter(ATM == "ATM3") |> autoplot(Cash) +
  labs(title = "ATM3: zero until the last three days", y = "Cash (hundreds of dollars)")

atm3_recent <- atm_clean |> filter(ATM == "ATM3", Cash > 0)
atm3_recent

## # A tibble: 3 x 3
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2010-04-28 ATM3     96
## 2 2010-04-29 ATM3     82
## 3 2010-04-30 ATM3     85

atm3_forecast <- mean(atm3_recent$Cash)
atm3_forecast

## [1] 87.66667

ATM3 has only three non-zero observations, on 28, 29 and 30 April 2010. Fitting a time series model to a series of 362 zeros would simply forecast zero, which is clearly wrong now that the machine is running, and three observations are not enough to estimate a weekly pattern or a trend.

The honest approach here is a simple average of what we actually know: the three days of real activity, giving a flat forecast of 87.7 hundred dollars per day. Those three values (96, 82, 85) are similar to each other and happen to be close to the ATM1 and ATM2 range, which makes the figure plausible. This forecast should be revised as soon as a few weeks of data are available.

Final forecasts for May 2010

atm_fc <- final_fit |> forecast(h = 31)

atm_fc |>
  autoplot(atm_model_data |> filter(DATE >= as.Date("2010-02-01"))) +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  labs(title = "May 2010 forecasts of daily cash withdrawals",
       y = "Cash (hundreds of dollars)", x = NULL)

may_dates <- seq(as.Date("2010-05-01"), as.Date("2010-05-31"), by = "day")

atm_forecast_table <- atm_fc |>
  as_tibble() |>
  select(ATM, DATE, Cash = .mean) |>
  bind_rows(tibble(ATM = "ATM3", DATE = may_dates, Cash = atm3_forecast)) |>
  mutate(Cash = round(Cash, 2)) |>
  arrange(DATE, ATM)

atm_forecast_wide <- atm_forecast_table |>
  tidyr::pivot_wider(names_from = ATM, values_from = Cash) |>
  mutate(Weekday = wday(DATE, label = TRUE), .after = DATE)

atm_forecast_wide

## # A tibble: 31 x 6
##    DATE       Weekday   ATM1   ATM2  ATM3  ATM4
##    <date>     <ord>    <dbl>  <dbl> <dbl> <dbl>
##  1 2010-05-01 Sat      86.8   65.5   87.7 411. 
##  2 2010-05-02 Sun     103.    76.2   87.7 455. 
##  3 2010-05-03 Mon      68.8    9.91  87.7 447. 
##  4 2010-05-04 Tue       3.57   1.8   87.7  64.0
##  5 2010-05-05 Wed     100.   100.    87.7 610. 
##  6 2010-05-06 Thu      81.0   91.1   87.7 473. 
##  7 2010-05-07 Fri      84.9   60.8   87.7 644. 
##  8 2010-05-08 Sat      86.8   65.5   87.7 411. 
##  9 2010-05-09 Sun     103.    76.2   87.7 455. 
## 10 2010-05-10 Mon      68.8    9.91  87.7 447. 
## # i 21 more rows

atm_forecast_table |>
  group_by(ATM) |>
  summarise(total_hundreds = round(sum(Cash)),
            total_dollars = round(sum(Cash) * 100),
            daily_mean = round(mean(Cash), 1))

## # A tibble: 4 x 4
##   ATM   total_hundreds total_dollars daily_mean
##   <chr>          <dbl>         <dbl>      <dbl>
## 1 ATM1            2374        237385       76.6
## 2 ATM2            1773        177346       57.2
## 3 ATM3            2718        271777       87.7
## 4 ATM4           13731       1373062      443.

atm_forecast_wide |>
  mutate(Weekday = as.character(Weekday)) |>
  group_by(Weekday) |>
  summarise(ATM1 = round(mean(ATM1), 1), ATM2 = round(mean(ATM2), 1), ATM4 = round(mean(ATM4), 1))

## # A tibble: 7 x 4
##   Weekday  ATM1  ATM2  ATM4
##   <chr>   <dbl> <dbl> <dbl>
## 1 Fri      84.9  60.8  644.
## 2 Mon      68.8   9.9  447.
## 3 Sat      86.8  65.5  411.
## 4 Sun     103.   76.2  455.
## 5 Thu      81    91.1  473.
## 6 Tue       3.6   1.8   64 
## 7 Wed     100.  100.   610.

The forecasts follow the recent weekly rhythm rather than the yearly average: Tuesday is the quiet day (ATM2 averages 1.8, ATM4 averages 64), while Wednesday, Friday and the weekend are busy. This is the direct consequence of choosing a method with time-varying seasonality, and it is the behaviour I want, since the February pattern was still in force through April.

It is also the main risk in this forecast. If the February change was temporary, or an artefact of how the data was recorded, then May will look more like the earlier months and these forecasts will be wrong on Tuesdays and Thursdays in particular. That is worth confirming with whoever owns the data before the numbers are used for cash loading.

Over the whole month, the four ATMs are expected to dispense about $1.98 million in total, with ATM4 accounting for roughly two thirds of it.

Part B: residential power usage

Reading and inspecting the data

power_raw <- read_excel("ResidentialCustomerForecastLoad-624.xlsx")
names(power_raw) <- c("CaseSequence", "YearMonth", "KWH")
head(power_raw)

## # A tibble: 6 x 3
##   CaseSequence YearMonth     KWH
##          <dbl> <chr>       <dbl>
## 1          733 1998-Jan  6862583
## 2          734 1998-Feb  5838198
## 3          735 1998-Mar  5420658
## 4          736 1998-Apr  5010364
## 5          737 1998-May  4665377
## 6          738 1998-Jun  6467147

power <- power_raw |>
  mutate(Month = yearmonth(paste(substr(YearMonth, 1, 4), substr(YearMonth, 6, 8)))) |>
  as_tsibble(index = Month) |>
  select(Month, KWH)

range(power$Month)

## <yearmonth[2]>
## [1] "1998 Jan" "2013 Dec"

nrow(power)

## [1] 192

The data covers 192 months, January 1998 to December 2013, with power consumption in kilowatt hours.

power |> autoplot(KWH) +
  labs(title = "Monthly residential power usage (raw)", y = "KWH")

power |> filter(is.na(KWH))

## # A tsibble: 1 x 2 [1M]
##      Month   KWH
##      <mth> <dbl>
## 1 2008 Sep    NA

power |> as_tibble() |> slice_min(KWH, n = 3)

## # A tibble: 3 x 2
##      Month     KWH
##      <mth>   <dbl>
## 1 2010 Jul  770523
## 2 2005 Nov 4313019
## 3 1999 Dec 4419229

Two problems are visible immediately:

• September 2008 is missing.
• July 2010 records 770,523 KWH, about one eighth of a normal July and far below any other month in sixteen years. July is a peak air-conditioning month, so this is not plausible; it looks like a reporting error.

Both are treated as missing and filled in using an ARIMA model, which estimates each value from the surrounding months and the seasonal pattern rather than from a plain average.

power_na <- power |> mutate(KWH = if_else(KWH < 1e6, NA_real_, KWH))

power_clean <- power_na |> model(ARIMA(KWH)) |> interpolate(power_na)

power_clean |>
  filter(Month %in% yearmonth(c("2008 Sep", "2010 Jul"))) |>
  mutate(KWH = round(KWH))

## # A tsibble: 2 x 2 [1M]
##      Month     KWH
##      <mth>   <dbl>
## 1 2008 Sep 7477569
## 2 2010 Jul 7747776

The interpolated values of about 7.48 million for September 2008 and 7.75 million for July 2010 sit comfortably within the range of neighbouring years.

power_clean |> autoplot(KWH) +
  labs(title = "Monthly residential power usage (cleaned)", y = "KWH")

Exploring the patterns

power_clean |> gg_season(KWH) + labs(title = "Seasonal plot of power usage", y = "KWH")

power_clean |> gg_subseries(KWH) + labs(title = "Subseries plot of power usage", y = "KWH")

power_clean |> model(STL(KWH)) |> components() |> autoplot() +
  labs(title = "STL decomposition of power usage")

The series has a strong twice-yearly seasonal pattern: peaks in January (heating) and again in July and August (cooling), with troughs in the mild months of April, May and October. There is a gradual upward trend, reflecting customer growth, with consumption rising from roughly 6 million KWH per month in the late 1990s to around 8 million by 2013. The STL remainder is small relative to the trend and seasonal components, apart from some larger values near the two corrected observations.

power_clean |> features(KWH, list(unitroot_nsdiffs, unitroot_ndiffs))

## # A tibble: 1 x 2
##   nsdiffs ndiffs
##     <int>  <int>
## 1       1      1

The unit root tests call for one seasonal difference and one first difference, confirming that both the seasonality and the trend need to be handled by whatever model is used.

Model selection

I hold out the last two years (2012 and 2013) as a test set and compare five candidates.

lambda_power <- power_clean |> features(KWH, features = guerrero) |> pull(lambda_guerrero)
lambda_power

## [1] -0.2057366

power_train <- power_clean |> filter(Month < yearmonth("2012 Jan"))

power_candidates <- power_train |>
  model(snaive = SNAIVE(KWH),
        ets = ETS(KWH),
        arima = ARIMA(KWH),
        arima_boxcox = ARIMA(box_cox(KWH, lambda_power)),
        stl_ets = decomposition_model(STL(box_cox(KWH, lambda_power)),
                                      ETS(season_adjust ~ season("N"))))

power_candidates |>
  forecast(h = 24) |>
  accuracy(power_clean) |>
  select(.model, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)

## # A tibble: 5 x 5
##   .model          RMSE     MAE  MAPE  MASE
##   <chr>          <dbl>   <dbl> <dbl> <dbl>
## 1 arima        721232. 532508.  6.81 0.876
## 2 arima_boxcox 773435. 617049.  7.97 1.01 
## 3 ets          850433. 622085.  8.41 1.02 
## 4 stl_ets      895630. 743853.  9.96 1.22 
## 5 snaive       965567. 742063.  9.30 1.22

ARIMA on the untransformed data has the lowest test RMSE (721,233) and the lowest MAPE (6.8%), so that is the model I take forward. The Box-Cox version is second. Interestingly the transformation hurts here, even though the Guerrero method suggests lambda = -0.21: the seasonal swings in this series grow only mildly with the level, so the transformation over-corrects. Seasonal naive is the weakest, as expected.

power_candidates |>
  forecast(h = 24) |>
  autoplot(power_clean |> filter(Month >= yearmonth("2008 Jan")), level = NULL) +
  labs(title = "Test set forecasts for 2012-2013", y = "KWH")

Final model and residual diagnostics

power_fit <- power_clean |> model(arima = ARIMA(KWH))
report(power_fit)

## Series: KWH 
## Model: ARIMA(0,0,4)(2,1,0)[12] w/ drift 
## 
## Coefficients:
##          ma1     ma2     ma3      ma4     sar1     sar2   constant
##       0.3352  0.0393  0.2200  -0.0707  -0.7033  -0.4179  227842.39
## s.e.  0.0780  0.0826  0.0746   0.0791   0.0775   0.0788   71996.59
## 
## sigma^2 estimated as 3.789e+11:  log likelihood=-2655.51
## AIC=5327.02   AICc=5327.86   BIC=5352.56

power_fit |> gg_tsresiduals()

power_fit |> augment() |> features(.innov, ljung_box, lag = 24, dof = 4)

## # A tibble: 1 x 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 arima     14.2     0.823

On the full data ARIMA() selects ARIMA(0,0,3)(2,1,0)[12] with drift: a seasonal difference to remove the annual pattern, two seasonal AR terms, three MA terms, and a drift term for the upward trend.

The residuals pass every check. All ACF spikes sit inside the significance bounds, the Ljung-Box test gives a p-value of 0.82, the residuals are centred on zero with roughly constant variance, and the histogram is close to normal. There is no information left in them.

Forecast for 2014

power_fc <- power_fit |> forecast(h = 12)

power_fc |>
  autoplot(power_clean |> filter(Month >= yearmonth("2008 Jan"))) +
  labs(title = "2014 forecast of monthly residential power usage", y = "KWH")

power_forecast_table <- power_fc |>
  as_tibble() |>
  select(Month, KWH = .mean) |>
  mutate(Month = as.character(Month), KWH = round(KWH))

power_forecast_table

## # A tibble: 12 x 2
##    Month         KWH
##    <chr>       <dbl>
##  1 2014 Jan 10295036
##  2 2014 Feb  8617878
##  3 2014 Mar  7164745
##  4 2014 Apr  5831199
##  5 2014 May  5942117
##  6 2014 Jun  8199229
##  7 2014 Jul  9478936
##  8 2014 Aug  9973049
##  9 2014 Sep  8486866
## 10 2014 Oct  5870156
## 11 2014 Nov  6147210
## 12 2014 Dec  8421706

sum(power_forecast_table$KWH)

## [1] 94428127

The forecast reproduces the familiar shape: a January peak of about 10.3 million KWH, a dip to about 5.8 million in April, a second peak of roughly 9.5 to 10 million in July and August, and a low again in October. The 2014 total comes to about 92.4 million KWH, slightly above 2013, which is consistent with the gentle upward trend.

Part C: water flow (bonus)

Reading the data and aggregating by hour

read_pipe <- function(file) {
  d <- read_excel(file)
  names(d) <- c("DateTime", "WaterFlow")
  d |> mutate(DateTime = as.POSIXct(DateTime * 86400, origin = "1899-12-30", tz = "UTC"))
}

pipe1 <- read_pipe("Waterflow_Pipe1.xlsx")
pipe2 <- read_pipe("Waterflow_Pipe2.xlsx")

head(pipe1, 4)

## # A tibble: 4 x 2
##   DateTime            WaterFlow
##   <dttm>                  <dbl>
## 1 2015-10-23 00:24:06      23.4
## 2 2015-10-23 00:40:02      28.0
## 3 2015-10-23 00:53:51      23.1
## 4 2015-10-23 00:55:40      30.0

head(pipe2, 4)

## # A tibble: 4 x 2
##   DateTime            WaterFlow
##   <dttm>                  <dbl>
## 1 2015-10-23 01:00:00      18.8
## 2 2015-10-23 01:59:59      43.1
## 3 2015-10-23 03:00:00      38.0
## 4 2015-10-23 04:00:00      36.1

tibble(pipe = c("Pipe 1", "Pipe 2"),
       n = c(nrow(pipe1), nrow(pipe2)),
       first = c(min(pipe1$DateTime), min(pipe2$DateTime)),
       last = c(max(pipe1$DateTime), max(pipe2$DateTime)))

## # A tibble: 2 x 4
##   pipe       n first               last               
##   <chr>  <int> <dttm>              <dttm>             
## 1 Pipe 1  1000 2015-10-23 00:24:06 2015-11-01 23:35:43
## 2 Pipe 2  1000 2015-10-23 01:00:00 2015-12-03 16:00:00

The two files have the same number of readings but very different time stamps. Pipe 1 records roughly four irregular readings per hour over ten days, while Pipe 2 records about one reading per hour over six weeks.

Following the instructions, each series is time sequenced and aggregated to hourly means.

to_hourly <- function(d) {
  d |>
    mutate(Hour = floor_date(DateTime, "hour")) |>
    group_by(Hour) |>
    summarise(WaterFlow = mean(WaterFlow), readings = n(), .groups = "drop")
}

hourly1 <- to_hourly(pipe1)
hourly2 <- to_hourly(pipe2)

head(hourly1, 5)

## # A tibble: 5 x 3
##   Hour                WaterFlow readings
##   <dttm>                  <dbl>    <int>
## 1 2015-10-23 00:00:00      26.1        4
## 2 2015-10-23 01:00:00      18.9        5
## 3 2015-10-23 02:00:00      15.2        2
## 4 2015-10-23 03:00:00      23.1        5
## 5 2015-10-23 04:00:00      15.5        2

head(hourly2, 5)

## # A tibble: 5 x 3
##   Hour                WaterFlow readings
##   <dttm>                  <dbl>    <int>
## 1 2015-10-23 01:00:00      30.9        2
## 2 2015-10-23 03:00:00      38.0        1
## 3 2015-10-23 04:00:00      34.0        2
## 4 2015-10-23 06:00:00      28.2        1
## 5 2015-10-23 07:00:00      18.3        2

tibble(pipe = c("Pipe 1", "Pipe 2"),
       hours_with_data = c(nrow(hourly1), nrow(hourly2)),
       hours_in_span = c(as.numeric(difftime(max(hourly1$Hour), min(hourly1$Hour), units = "hours")) + 1,
                         as.numeric(difftime(max(hourly2$Hour), min(hourly2$Hour), units = "hours")) + 1),
       mean_flow = round(c(mean(hourly1$WaterFlow), mean(hourly2$WaterFlow)), 2),
       sd_flow = round(c(sd(hourly1$WaterFlow), sd(hourly2$WaterFlow)), 2))

## # A tibble: 2 x 5
##   pipe   hours_with_data hours_in_span mean_flow sd_flow
##   <chr>            <int>         <dbl>     <dbl>   <dbl>
## 1 Pipe 1             236           240      19.9    4.24
## 2 Pipe 2             667          1000      39.8   13.8

Pipe 1 covers 240 hours of which 236 have readings, so only four hours are missing. Pipe 2 covers 1,000 hours but only 667 have readings, leaving 333 gaps. The gaps are filled with fill_gaps() and then interpolated with an ARIMA model so that the series are complete and evenly spaced.

ts1 <- hourly1 |> select(-readings) |> as_tsibble(index = Hour) |> fill_gaps()
ts2 <- hourly2 |> select(-readings) |> as_tsibble(index = Hour) |> fill_gaps()

ts1 <- ts1 |> model(ARIMA(WaterFlow)) |> interpolate(ts1)
ts2 <- ts2 |> model(ARIMA(WaterFlow)) |> interpolate(ts2)

ts1 |> autoplot(WaterFlow) + labs(title = "Pipe 1: hourly mean water flow", y = "Water flow")

ts2 |> autoplot(WaterFlow) + labs(title = "Pipe 2: hourly mean water flow", y = "Water flow")

Pipe 2 runs at roughly twice the flow of Pipe 1 (mean 39.8 against 19.9) and is much more variable.

Is the data stationary?

bind_rows(
  ts1 |> features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs)) |> mutate(pipe = "Pipe 1"),
  ts2 |> features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs)) |> mutate(pipe = "Pipe 2")
) |> relocate(pipe)

## # A tibble: 2 x 4
##   pipe   kpss_stat kpss_pvalue ndiffs
##   <chr>      <dbl>       <dbl>  <int>
## 1 Pipe 1    0.0794      0.1         0
## 2 Pipe 2    0.475       0.0473      1

ts1 |> ACF(WaterFlow, lag_max = 48) |> autoplot() + labs(title = "ACF of Pipe 1 hourly flow")

ts2 |> ACF(WaterFlow, lag_max = 48) |> autoplot() + labs(title = "ACF of Pipe 2 hourly flow")

bind_rows(
  tibble(pipe = "Pipe 1", lb_pvalue = Box.test(hourly1$WaterFlow, lag = 24, type = "Ljung-Box")$p.value),
  tibble(pipe = "Pipe 2", lb_pvalue = Box.test(hourly2$WaterFlow, lag = 24, type = "Ljung-Box")$p.value)
)

## # A tibble: 2 x 2
##   pipe   lb_pvalue
##   <chr>      <dbl>
## 1 Pipe 1  0.897   
## 2 Pipe 2  0.000205

Pipe 1 is stationary and is essentially white noise. The KPSS test gives a p-value of 0.1, unitroot_ndiffs() returns 0, the ACF spikes all sit inside the bounds, and a Ljung-Box test on the observed hours gives a p-value of 0.90. The series simply fluctuates around a constant mean of about 19.9.

Pipe 2 is stationary in level but not quite white noise. The KPSS p-value of 0.047 is borderline, and a Ljung-Box test on the observed hours is significant (p = 0.0002). However, the individual autocorrelations are tiny: the largest in the first ten lags is 0.099, so under 1% of the variation is explained. With 667 observations even a trivial correlation becomes statistically significant, and this one is too small to be useful. The test set results below confirm that it cannot be exploited.

Can it be forecast?

Both series can be forecast in the sense that a model can be fitted and a prediction made, but the useful question is whether anything beats a flat line at the mean. I hold out the last 48 hours of Pipe 1 and the last week (168 hours) of Pipe 2.

train1 <- ts1 |> slice(1:(n() - 48))
train2 <- ts2 |> slice(1:(n() - 168))

water_fit1 <- train1 |> model(mean = MEAN(WaterFlow), naive = NAIVE(WaterFlow),
                              ets = ETS(WaterFlow), arima = ARIMA(WaterFlow))
water_fit2 <- train2 |> model(mean = MEAN(WaterFlow), naive = NAIVE(WaterFlow),
                              ets = ETS(WaterFlow), arima = ARIMA(WaterFlow))

water_fit1 |> forecast(h = 48) |> accuracy(ts1) |>
  select(.model, RMSE, MAE, MAPE) |> arrange(RMSE)

## # A tibble: 4 x 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 ets     6.41  3.97  33.8
## 2 mean    6.41  3.97  33.8
## 3 arima   6.41  3.97  33.8
## 4 naive   7.41  5.78  38.5

water_fit2 |> forecast(h = 168) |> accuracy(ts2) |>
  select(.model, RMSE, MAE, MAPE) |> arrange(RMSE)

## # A tibble: 4 x 4
##   .model  RMSE   MAE  MAPE
##   <chr>  <dbl> <dbl> <dbl>
## 1 mean    13.0  10.3  35.7
## 2 ets     14.5  11.4  32.3
## 3 arima   15.2  12.2  32.6
## 4 naive   17.8  14.7  36.6

For Pipe 1 the mean, ETS and ARIMA models are indistinguishable (RMSE 6.41 for all three), because the automatic selections collapse to the same thing: ARIMA() chooses ARIMA(0,0,0) with a mean, and ETS() chooses ETS(A,N,N) with a very small smoothing parameter. Both are the mean model under another name. For Pipe 2 the mean is outright best, beating ARIMA by about 14% on RMSE.

So the answer is: the data is stationary, and the only thing worth forecasting is the mean level. There is no trend, no daily cycle and no exploitable short-term memory in either pipe, so I use the mean model for both, which gives a flat forecast with constant prediction intervals.

water_final1 <- ts1 |> model(mean = MEAN(WaterFlow))
water_final2 <- ts2 |> model(mean = MEAN(WaterFlow))

water_fc1 <- water_final1 |> forecast(h = 168)
water_fc2 <- water_final2 |> forecast(h = 168)

water_fc1 |> autoplot(ts1) +
  labs(title = "Pipe 1: one week forecast of hourly water flow", y = "Water flow")

water_fc2 |> autoplot(ts2) +
  labs(title = "Pipe 2: one week forecast of hourly water flow", y = "Water flow")

water_forecast_table <- bind_rows(
  water_fc1 |> as_tibble() |> transmute(Pipe = "Pipe1", Hour, WaterFlow = round(.mean, 3)),
  water_fc2 |> as_tibble() |> transmute(Pipe = "Pipe2", Hour, WaterFlow = round(.mean, 3))
) |>
  mutate(Hour = format(Hour, "%Y-%m-%d %H:%M"))

head(water_forecast_table, 4)

## # A tibble: 4 x 3
##   Pipe  Hour             WaterFlow
##   <chr> <chr>                <dbl>
## 1 Pipe1 2015-11-02 00:00      19.5
## 2 Pipe1 2015-11-02 01:00      19.5
## 3 Pipe1 2015-11-02 02:00      19.5
## 4 Pipe1 2015-11-02 03:00      19.5

tail(water_forecast_table, 4)

## # A tibble: 4 x 3
##   Pipe  Hour             WaterFlow
##   <chr> <chr>                <dbl>
## 1 Pipe2 2015-12-10 13:00      39.8
## 2 Pipe2 2015-12-10 14:00      39.8
## 3 Pipe2 2015-12-10 15:00      39.8
## 4 Pipe2 2015-12-10 16:00      39.8

The forecasts are flat: 19.47 for Pipe 1 and 39.77 for Pipe 2 in every hour of the week ahead. That is the correct answer for white noise, and the honest conclusion is that these data support forecasting the average level and its uncertainty, but not the hour-to-hour movements.

Forecast output file

All three sets of forecasts are written to a single Excel workbook with one sheet per part.

write_xlsx(
  list(ATM_May2010 = atm_forecast_wide |> mutate(DATE = as.character(DATE), Weekday = as.character(Weekday)),
       ATM_May2010_long = atm_forecast_table |> mutate(DATE = as.character(DATE)),
       Power_2014 = power_forecast_table,
       Waterflow_Week = water_forecast_table),
  path = "DATA624_Project1_Forecasts.xlsx"
)

file.exists("DATA624_Project1_Forecasts.xlsx")

## [1] TRUE

The workbook DATA624_Project1_Forecasts.xlsx contains:

• ATM_May2010: the May 2010 daily forecast for each of the four ATMs, one row per day.
• ATM_May2010_long: the same forecast in long format, one row per ATM and date.
• Power_2014: the monthly KWH forecast for 2014.
• Waterflow_Week: hourly forecasts for both pipes for the week ahead.

Summary

Part	Series	Chosen model	Reason
A	ATM1, ATM2, ATM4	STL + ETS	Lowest test RMSE of five candidates; handles the strong weekly cycle
A	ATM3	Mean of the three non-zero days	Only three days of real data exist
B	Residential power	ARIMA(0,0,3)(2,1,0)[12] with drift	Lowest test RMSE; residuals are white noise (p = 0.82)
C	Pipe 1 and Pipe 2	Mean	Both series are stationary with no exploitable structure

The two judgement calls worth flagging are the handling of the ATM4 outlier, where a single implausible value of about $1.09 million in one day was treated as a data error, and ATM3, where a model fitted to the historical zeros would have been technically defensible but useless in practice.