All forecasting uses the fpp3 framework (tsibble, fable, feasts). Each part follows the same process: clean the data, look at it, hold out the most recent period, compare candidate models on that holdout, then refit the winner on all of the data and forecast. All forecasts are saved to project1_forecasts.xlsx.

Part A — ATM cash withdrawals

Data and cleaning

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

atm_raw |> filter(is.na(ATM)) |> summarise(rows = n(), from = min(DATE), to = max(DATE))
## # A tibble: 1 × 3
##    rows from       to        
##   <int> <date>     <date>    
## 1    14 2010-05-01 2010-05-14
atm_raw |> filter(!is.na(ATM), is.na(Cash))
## # A tibble: 5 × 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

Excel stores dates as day counts since 1899-12-30, so DATE is converted first. The file has 1474 rows. The 14 rows with no ATM are empty placeholders for May 2010, the forecast period, so they are dropped. That leaves 365 days per ATM, from 2009-05-01 to 2010-04-30. Cash is in hundreds of dollars.

atm <- atm_raw |>
  filter(!is.na(ATM)) |>
  as_tsibble(index = DATE, key = ATM)

atm |>
  autoplot(Cash) +
  facet_wrap(~ ATM, ncol = 1, scales = "free_y") +
  theme(legend.position = "none") +
  labs(title = "Daily cash withdrawals by ATM (hundreds of dollars)", y = NULL)

atm |>
  as_tibble() |>
  mutate(weekday = wday(DATE, label = TRUE)) |>
  group_by(ATM, weekday) |>
  summarise(median_cash = median(Cash, na.rm = TRUE), .groups = "drop") |>
  pivot_wider(names_from = weekday, values_from = median_cash) |>
  knitr::kable(digits = 0, caption = "Median withdrawals by weekday")
Median withdrawals by weekday
ATM Sun Mon Tue Wed Thu Fri Sat
ATM1 108 85 99 78 16 98 96
ATM2 64 64 82 28 8 104 72
ATM3 0 0 0 0 0 0 0
ATM4 533 481 307 406 67 572 441

What the plots show:

  • ATM1 and ATM2 have no trend and a strong weekly pattern, with one day each week far below the rest. Over the full year that low day is Thursday, but the table below shows it moved to Tuesday from February 2010 on. Any model used here has to let the weekly pattern change over time.
  • ATM3 reads 0 on every day except the last 3 (April 28 to 30, 2010). It looks like a machine that was just installed or switched on.
  • ATM4 has the same weekly shape at a higher level, plus one extreme value of 10,920 on 2010-02-09. That is more than 6 times the next largest day, which points to a data entry error or a one-off event, not normal demand.
  • There are 5 missing days, all at ATM1 and ATM2 in June 2009.

Cleaning: the ATM4 outlier is set to missing. Then each missing value is filled with that ATM’s median for the same weekday, which keeps the weekly pattern intact.

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

atm_clean |> filter(ATM != "ATM3") |> features(Cash, guerrero)
## # A tibble: 3 × 2
##   ATM   lambda_guerrero
##   <chr>           <dbl>
## 1 ATM1            0.249
## 2 ATM2            0.673
## 3 ATM4            0.451

The Guerrero method gives a different Box-Cox \(\lambda\) for each machine, so each machine is modeled with its own transformation. ATM2 has two days with \(0\), which is fine because its \(\lambda\) is positive.

atm_clean |>
  filter(ATM != "ATM3") |>
  as_tibble() |>
  mutate(Month = yearmonth(DATE), weekday = wday(DATE, label = TRUE)) |>
  group_by(ATM, Month) |>
  slice_min(Cash, n = 1, with_ties = FALSE) |>
  select(ATM, Month, weekday) |>
  pivot_wider(names_from = ATM, values_from = weekday) |>
  knitr::kable(caption = "Lowest weekday in each month")
Lowest weekday in each month
Month ATM1 ATM2 ATM4
2009 May Thu Thu Sat
2009 Jun Fri Thu Fri
2009 Jul Thu Thu Thu
2009 Aug Thu Thu Tue
2009 Sep Fri Fri Fri
2009 Oct Thu Sun Sat
2009 Nov Mon Mon Sun
2009 Dec Wed Thu Thu
2010 Jan Thu Thu Thu
2010 Feb Tue Thu Tue
2010 Mar Tue Tue Tue
2010 Apr Tue Tue Tue

ATM1, ATM2 and ATM4: model comparison

The data run through April 30, 2010. The last 30 days (April 2010) are held out, and 4 models are fit on the data before April:

  • SNAIVE: each day equals the same weekday last week. This is the benchmark.
  • ETS on the Box-Cox scale, with the form chosen by AICc.
  • ARIMA on the Box-Cox scale, with the seasonal ARIMA order chosen automatically and a weekly period.
  • STL + ETS: STL removes the weekly pattern, ETS forecasts the seasonally adjusted series, and the weekly pattern is added back. The STL seasonal window is set to 7 weeks so the weekly pattern can follow the shift from Thursday to Tuesday; with the default window, the forecasts kept placing the low day on Thursday.

Models not used: a trend model, because none of the series has a trend; Holt-Winters with multiplicative seasonality, because the Box-Cox transform already handles the changing spread; and models with yearly seasonality, because one year of daily data is not enough to estimate a yearly pattern.

atm_models <- function(data, lam) {
  data |>
    model(
      SNAIVE  = SNAIVE(Cash),
      ETS     = ETS(box_cox(Cash, lam)),
      ARIMA   = ARIMA(box_cox(Cash, lam)),
      STL_ETS = decomposition_model(
        STL(box_cox(Cash, lam) ~ season(period = 7, window = 7), robust = TRUE),
        ETS(season_adjust ~ season("N"))
      )
    )
}

holdout_start <- as.Date("2010-04-01")
atm_ids <- c("ATM1", "ATM2", "ATM4")

atm_results <- lapply(atm_ids, function(id) {
  d <- atm_clean |> filter(ATM == id)
  lam <- d |> features(Cash, guerrero) |> pull(lambda_guerrero)
  fit <- atm_models(d |> filter(DATE < holdout_start), lam)
  acc <- fit |> forecast(h = 30) |> accuracy(d) |> select(ATM, .model, RMSE, MAE, MAPE)
  list(id = id, lam = lam, data = d, acc = acc)
})

atm_acc <- bind_rows(lapply(atm_results, `[[`, "acc"))
atm_acc |>
  group_by(ATM) |>
  arrange(RMSE, .by_group = TRUE) |>
  knitr::kable(digits = 2, caption = "April 2010 holdout accuracy")
April 2010 holdout accuracy
ATM .model RMSE MAE MAPE
ATM1 ARIMA 14.33 12.06 45.44
ATM1 ETS 14.65 12.24 45.75
ATM1 SNAIVE 16.04 13.53 73.40
ATM1 STL_ETS 20.20 15.35 26.14
ATM2 ETS 19.73 14.68 109.25
ATM2 STL_ETS 20.59 13.77 83.95
ATM2 ARIMA 21.05 16.18 108.85
ATM2 SNAIVE 26.16 17.80 45.14
ATM4 ETS 245.91 208.77 652.33
ATM4 ARIMA 259.36 223.85 860.48
ATM4 STL_ETS 267.86 218.62 476.83
ATM4 SNAIVE 300.85 227.26 375.05 
atm_best <- atm_acc |>
  group_by(ATM) |>
  slice_min(RMSE, n = 1) |>
  ungroup() |>
  select(ATM, best_model = .model, RMSE)
knitr::kable(atm_best, digits = 2)
ATM best_model RMSE
ATM1 ARIMA 14.33
ATM2 ETS 19.73
ATM4 ETS 245.91

ARIMA wins for ATM1 and ETS wins for ATM2 and ATM4, and every winner beats the seasonal naive benchmark on RMSE. The top 2 models are close at every machine (for example 14.3 vs 14.7 at ATM1), so either would be a reasonable choice, and the winner is taken on RMSE. MAPE is reported but not used to choose: on the low days actual withdrawals are near zero, so a small miss in dollars becomes a huge percentage and MAPE swings wildly.

ATM1, ATM2 and ATM4: final models and residuals

Each machine’s best model is refit on all 365 days and forecast for the 31 days of May 2010.

atm_final <- lapply(atm_results, function(r) {
  best <- atm_best$best_model[atm_best$ATM == r$id]
  fit <- atm_models(r$data, r$lam) |> select(ATM, all_of(best))
  list(id = r$id, best = best, fit = fit, fc = forecast(fit, h = 31))
})

for (f in atm_final) {
  print(f$fit |> report())
}
## Series: Cash 
## Model: ARIMA(0,0,2)(0,1,1)[7] 
## Transformation: box_cox(Cash, lam) 
## 
## Coefficients:
##          ma1      ma2     sma1
##       0.1018  -0.1124  -0.6448
## s.e.  0.0524   0.0526   0.0427
## 
## sigma^2 estimated as 1.623:  log likelihood=-595.08
## AIC=1198.17   AICc=1198.28   BIC=1213.69
## # A mable: 1 x 2
## # Key:     ATM [1]
##   ATM                      ARIMA
##   <chr>                  <model>
## 1 ATM1  <ARIMA(0,0,2)(0,1,1)[7]>
## Series: Cash 
## Model: ETS(A,N,A) 
## Transformation: box_cox(Cash, lam) 
##   Smoothing parameters:
##     alpha = 0.0001000088 
##     gamma = 0.4006754 
## 
##   Initial states:
##      l[0]      s[0]     s[-1]    s[-2]    s[-3]    s[-4]    s[-5]    s[-6]
##  23.55319 -18.68819 -6.483719 2.372326 2.796957 5.325976 6.958391 7.718261
## 
##   sigma^2:  47.1024
## 
##      AIC     AICc      BIC 
## 3570.448 3571.069 3609.447 
## # A mable: 1 x 2
## # Key:     ATM [1]
##   ATM            ETS
##   <chr>      <model>
## 1 ATM2  <ETS(A,N,A)>
## Series: Cash 
## Model: ETS(M,N,A) 
## Transformation: box_cox(Cash, lam) 
##   Smoothing parameters:
##     alpha = 0.01389926 
##     gamma = 0.1729172 
## 
##   Initial states:
##      l[0]      s[0]    s[-1]    s[-2]    s[-3]    s[-4]    s[-5]   s[-6]
##  26.48241 -18.66685 -6.20789 3.688922 4.049661 7.173808 5.019026 4.94332
## 
##   sigma^2:  0.222
## 
##      AIC     AICc      BIC 
## 4025.598 4026.219 4064.597 
## # A mable: 1 x 2
## # Key:     ATM [1]
##   ATM            ETS
##   <chr>      <model>
## 1 ATM4  <ETS(M,N,A)>
for (f in atm_final) {
  print(f$fit |> gg_tsresiduals() + labs(title = paste(f$id, "residuals:", f$best)))
}

bind_rows(lapply(atm_final, function(f) {
  augment(f$fit) |> features(.innov, ljung_box, lag = 14, dof = 0)
})) |>
  knitr::kable(digits = 3, caption = "Ljung-Box test on residuals (lag 14)")
Ljung-Box test on residuals (lag 14)
ATM .model lb_stat lb_pvalue
ATM1 ARIMA 10.520 0.723
ATM2 ETS 32.139 0.004
ATM4 ETS 15.576 0.340

ATM1’s ARIMA\((0,0,2)(0,1,1)_7\) and ATM4’s ETS(M,N,A) leave residuals that look like white noise (Ljung-Box p-values of 0.72 and 0.34). ATM2’s ETS(A,N,A) still has some autocorrelation left (p = 0.004), so its prediction intervals are likely somewhat too narrow. On its low days the lower 95% bound also drops below zero, which is impossible for cash; in practice the floor is 0. Its point forecasts are still the most accurate of the 4 candidates on the holdout, so it is kept.

ATM3

ATM3 cannot be modeled like the others. With only 3 non-zero days there is no weekly pattern, trend, or variance to estimate. Any seasonal model would treat the 362 zeros before them as real demand and forecast close to zero. The more realistic reading is that the machine came online on April 28. The best available estimate for May is the mean of those 3 days, 87.7 per day, a flat forecast. Once a few weeks of data come in, ATM3 should be modeled like ATM1, with the weekly pattern.

atm3_level <- atm_clean |> filter(ATM == "ATM3") |> tail(3) |> pull(Cash) |> mean()
may_days <- seq(as.Date("2010-05-01"), as.Date("2010-05-31"), by = "day")
atm3_fc <- tibble(ATM = "ATM3", DATE = may_days, Forecast = atm3_level,
                  Model = "Mean of 3 active days")

May 2010 forecast

atm_fc_plot <- lapply(atm_final, function(f) {
  f$fc |>
    autoplot(atm_clean |> filter(ATM == f$id, DATE >= as.Date("2010-02-01"))) +
    labs(title = paste(f$id, "May 2010 forecast:", f$best), y = "Cash (hundreds)", x = NULL)
})
for (p in atm_fc_plot) print(p)

atm_fc <- bind_rows(
  bind_rows(lapply(atm_final, function(f) {
    f$fc |>
      as_tibble() |>
      transmute(ATM, DATE, Forecast = .mean, Model = f$best)
  })),
  atm3_fc
) |>
  arrange(ATM, DATE)

atm_fc |>
  group_by(ATM, Model) |>
  summarise(daily_mean = mean(Forecast), may_total_hundreds = sum(Forecast),
            may_total_dollars = sum(Forecast) * 100, .groups = "drop") |>
  knitr::kable(digits = 0, format.args = list(big.mark = ","),
               caption = "May 2010 forecast summary")
May 2010 forecast summary
ATM Model daily_mean may_total_hundreds may_total_dollars
ATM1 ARIMA 83 2,584 258,404
ATM2 ETS 62 1,927 192,718
ATM3 Mean of 3 active days 88 2,718 271,767
ATM4 ETS 420 13,032 1,303,165

The May forecasts keep the recent weekly shape: Tuesday is the low day at ATM1, ATM2 and ATM4, and ATM2 also has a low Monday, which matches its last 2 months of data. In total, the 4 ATMs are forecast to dispense about $2,026,054 in May 2010. ATM4 accounts for more than half of that. For cash planning, the daily values in the Excel file can be summed by restocking period, and the 80% and 95% intervals in the plots show how much buffer to keep. The ATM3 number is the least certain because it rests on 3 days of history.

Part B — Residential power usage

Data and cleaning

power_raw <- read_excel("ResidentialCustomerForecastLoad-624.xlsx") |>
  transmute(Month = yearmonth(`YYYY-MMM`), KWH) |>
  as_tsibble(index = Month)

power_raw |> filter(is.na(KWH) | KWH < 2e6)
## # A tsibble: 2 x 2 [1M]
##      Month    KWH
##      <mth>  <dbl>
## 1 2008 Sep     NA
## 2 2010 Jul 770523
power_raw |> autoplot(KWH) + labs(title = "Monthly residential power usage (KWH)", y = "KWH")

The series runs from January 1998 to December 2013. There are 2 data problems: September 2008 is missing, and July 2010 shows 770,523 KWH. That is about a tenth of a normal July and well below every other month in 16 years, so it is treated as an error. Both values are set to missing and filled in by an ARIMA model fit to the rest of the series (interpolate()), which uses the seasonal pattern and the nearby months.

power_na <- power_raw |> mutate(KWH = if_else(KWH < 2e6, NA_real_, KWH))
power <- power_na |>
  model(ARIMA(KWH)) |>
  interpolate(power_na)

power |> filter(Month %in% yearmonth(c("2008 Sep", "2010 Jul")))
## # A tsibble: 2 x 2 [1M]
##      Month      KWH
##      <mth>    <dbl>
## 1 2008 Sep 7477571.
## 2 2010 Jul 7747775.

Exploration

power |> gg_season(KWH, labels = "both") + labs(title = "Seasonal plot", y = "KWH")

power |>
  model(STL(KWH ~ season(window = "periodic"), robust = TRUE)) |>
  components() |>
  autoplot()

lambda_p <- power |> features(KWH, guerrero) |> pull(lambda_guerrero)
lambda_p
## [1] -0.2057365

Usage has 2 peaks a year: winter (January and February) and summer (July and August), with lows in spring (April and May) and fall (October and November). This pattern fits heating and air conditioning. The level is flat at about 74 to 77 million KWH a year through 2008, then rises to about 87 to 92 million KWH in 2011 to 2013, and the swings get larger as the level rises. The Guerrero \(\lambda\) is -0.21, close to a log transformation, which evens out that growing variance.

Model comparison

The last 2 years (2012 and 2013) are held out, and 5 models are compared:

power_train <- power |> filter(year(Month) < 2012)

power_models <- function(data) {
  data |>
    model(
      SNAIVE  = SNAIVE(KWH),
      ETS     = ETS(box_cox(KWH, lambda_p)),
      ARIMA   = ARIMA(box_cox(KWH, lambda_p)),
      STL_ETS = decomposition_model(
        STL(box_cox(KWH, lambda_p) ~ season(window = "periodic"), robust = TRUE),
        ETS(season_adjust ~ season("N"))
      ),
      TSLM    = TSLM(box_cox(KWH, lambda_p) ~ trend() + season())
    )
}

power_fit <- power_models(power_train)
power_acc <- power_fit |>
  forecast(h = 24) |>
  accuracy(power) |>
  select(.model, RMSE, MAE, MAPE) |>
  arrange(RMSE)
knitr::kable(power_acc, digits = 2, format.args = list(big.mark = ","),
             caption = "2012 to 2013 holdout accuracy")
2012 to 2013 holdout accuracy
.model RMSE MAE MAPE
ARIMA 773,435.0 617,048.8 7.97
STL_ETS 854,712.4 651,531.8 8.67
TSLM 885,627.0 598,796.6 7.60
SNAIVE 965,567.1 742,063.5 9.30
ETS 1,079,993.0 895,192.1 12.16 
power_fit |>
  select(-TSLM) |>
  forecast(h = 24) |>
  autoplot(power |> filter(year(Month) >= 2009), level = NULL) +
  labs(title = "Holdout forecasts vs actual, 2012 to 2013", y = "KWH")

ARIMA has the lowest RMSE on the 2012 to 2013 holdout, at about 8% MAPE. The regression model (TSLM) has a slightly lower MAPE but a higher RMSE, because a straight-line trend misses the step up in usage after 2009. ETS does worst, worse even than the seasonal naive benchmark. ARIMA is selected.

Final model and 2014 forecast

power_best <- power_acc$.model[1]
power_final <- power_models(power) |> select(all_of(power_best))
report(power_final)
## Series: KWH 
## Model: ARIMA(0,0,3)(2,1,0)[12] w/ drift 
## Transformation: box_cox(KWH, lambda_p) 
## 
## Coefficients:
##          ma1     ma2     ma3     sar1     sar2  constant
##       0.2674  0.0620  0.2196  -0.7230  -0.3805    0.0013
## s.e.  0.0759  0.0849  0.0710   0.0745   0.0765    0.0005
## 
## sigma^2 estimated as 1.414e-05:  log likelihood=758.83
## AIC=-1503.66   AICc=-1503.01   BIC=-1481.31
power_final |> gg_tsresiduals()

augment(power_final) |> features(.innov, ljung_box, lag = 24, dof = 0)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 ARIMA     16.3     0.879
power_fc <- power_final |> forecast(h = 12)
power_fc |>
  autoplot(power |> filter(year(Month) >= 2008)) +
  labs(title = paste("2014 monthly forecast:", power_best), y = "KWH")

power_fc_tbl <- power_fc |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  transmute(Month = format(Month, "%Y-%b"), Forecast_KWH = .mean,
            Lower_95 = `95%_lower`, Upper_95 = `95%_upper`)
knitr::kable(power_fc_tbl, digits = 0, format.args = list(big.mark = ","))
Month Forecast_KWH Lower_95 Upper_95
2014-Jan 10,414,868 8,468,970 12,750,740
2014-Feb 8,780,044 7,140,489 10,747,965
2014-Mar 7,079,855 5,807,942 8,594,948
2014-Apr 5,982,856 4,920,405 7,245,646
2014-May 5,941,894 4,888,079 7,194,114
2014-Jun 8,305,110 6,736,841 10,191,724
2014-Jul 9,585,835 7,726,695 11,834,786
2014-Aug 10,146,851 8,157,979 12,558,135
2014-Sep 8,521,497 6,904,626 10,468,488
2014-Oct 5,856,902 4,820,969 7,087,245
2014-Nov 6,131,467 5,037,595 7,432,740
2014-Dec 8,304,252 6,736,175 10,190,627

The final model is ARIMA\((0,0,3)(2,1,0)_{12}\) with drift on the Box-Cox scale. The residuals pass the Ljung-Box test (p = 0.88), so the model has captured the autocorrelation. The 2014 forecast repeats the 2 seasonal peaks: about 10.4 million KWH in January and 10.1 million in August, with lows of about 5.9 million in May and October. That adds up to about 95.1 million KWH for the year, a little above 2013, following the drift. The 95% intervals are about plus or minus 20% each month.

Part C — Waterflow (bonus)

Time-base and aggregate to hourly

read_pipe <- function(path, pipe) {
  read_excel(path) |>
    rename(DateTime = `Date Time`) |>
    mutate(DateTime = as.POSIXct(round(DateTime * 86400), origin = "1899-12-30", tz = "UTC"),
           Pipe = pipe)
}

water_raw <- bind_rows(read_pipe("Waterflow_Pipe1.xlsx", "Pipe1"),
                       read_pipe("Waterflow_Pipe2.xlsx", "Pipe2"))

water_raw |>
  group_by(Pipe) |>
  summarise(readings = n(), first = min(DateTime), last = max(DateTime))
## # A tibble: 2 × 4
##   Pipe  readings first               last               
##   <chr>    <int> <dttm>              <dttm>             
## 1 Pipe1     1000 2015-10-23 00:24:06 2015-11-01 23:35:43
## 2 Pipe2     1000 2015-10-23 01:00:00 2015-12-03 16:00:00

The timestamps are Excel day fractions, so they are converted to date-times and rounded to the second. Pipe1 has several readings at irregular times within each hour over about 10 days. Pipe2 has exactly 1 reading per hour over about 41 days. Each reading is assigned to the hour it falls in (floor_date), and readings in the same hour are averaged:

water <- water_raw |>
  mutate(Hour = floor_date(DateTime, "hour")) |>
  group_by(Pipe, Hour) |>
  summarise(WaterFlow = mean(WaterFlow), readings = n(), .groups = "drop") |>
  as_tsibble(index = Hour, key = Pipe)

water |> as_tibble() |> group_by(Pipe) |>
  summarise(hours = n(), mean_readings_per_hour = mean(readings))
## # A tibble: 2 × 3
##   Pipe  hours mean_readings_per_hour
##   <chr> <int>                  <dbl>
## 1 Pipe1   236                   4.24
## 2 Pipe2  1000                   1
water |> count_gaps()
## # A tibble: 4 × 4
##   Pipe  .from               .to                    .n
##   <chr> <dttm>              <dttm>              <int>
## 1 Pipe1 2015-10-27 17:00:00 2015-10-27 17:00:00     1
## 2 Pipe1 2015-10-28 00:00:00 2015-10-28 00:00:00     1
## 3 Pipe1 2015-11-01 05:00:00 2015-11-01 05:00:00     1
## 4 Pipe1 2015-11-01 09:00:00 2015-11-01 09:00:00     1
water_filled <- water |>
  fill_gaps() |>
  group_by_key() |>
  mutate(WaterFlow = approx(seq_along(WaterFlow), WaterFlow, seq_along(WaterFlow))$y) |>
  ungroup()

water_filled |> filter(Pipe == "Pipe1", Hour %in% (count_gaps(water) |> pull(.from)))
## # A tsibble: 4 x 4 [1h] <UTC>
## # Key:       Pipe [1]
##   Pipe  Hour                WaterFlow readings
##   <chr> <dttm>                  <dbl>    <int>
## 1 Pipe1 2015-10-27 17:00:00      22.5       NA
## 2 Pipe1 2015-10-28 00:00:00      18.8       NA
## 3 Pipe1 2015-11-01 05:00:00      16.8       NA
## 4 Pipe1 2015-11-01 09:00:00      15.4       NA

After aggregation Pipe1 has 236 hours, with an average of about 4 readings averaged into each hour. Pipe2 keeps its 1,000 hourly readings unchanged. Pipe1 is missing 4 single hours, which are filled by linear interpolation between the neighboring hours so that the series is regular for modeling. ARIMA interpolation was tried first, but for this mean-only series it returned values around \(-5.5\), below any real reading, so it was not used. The 2 pipes are kept separate because they cover different periods (10 days vs 41 days) and their levels differ by about 2 times.

water_filled |>
  autoplot(WaterFlow) +
  facet_wrap(~ Pipe, ncol = 1, scales = "free_x") +
  theme(legend.position = "none") +
  labs(title = "Hourly mean waterflow", y = "WaterFlow", x = NULL)

Stationarity

water_filled |> filter(Pipe == "Pipe1") |> ACF(WaterFlow, lag_max = 48) |> autoplot() +
  labs(title = "Pipe1 ACF")

water_filled |> filter(Pipe == "Pipe2") |> ACF(WaterFlow, lag_max = 48) |> autoplot() +
  labs(title = "Pipe2 ACF")

water_filled |>
  features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs, unitroot_nsdiffs)) |>
  knitr::kable(digits = 3)
Pipe kpss_stat kpss_pvalue ndiffs nsdiffs
Pipe1 0.193 0.1 0 0
Pipe2 0.105 0.1 0 0

Both series are stationary. Neither plot shows a trend or a change in spread. The KPSS test does not reject stationarity for either pipe (p-value of 0.1 or above, which is the top of the table), and unitroot_ndiffs and unitroot_nsdiffs both say no differencing is needed. Pipe1’s ACF has no spikes past the bounds, so its hourly flow is close to white noise around a mean of about 19.9. Pipe2 has a few small spikes, all below 0.1, including one at lag 24 that hints at a weak daily effect.

One-week forecast

water_fit <- water_filled |> model(ARIMA = ARIMA(WaterFlow))
water_fit
## # A mable: 2 x 2
## # Key:     Pipe [2]
##   Pipe                              ARIMA
##   <chr>                           <model>
## 1 Pipe1            <ARIMA(0,0,0) w/ mean>
## 2 Pipe2 <ARIMA(0,0,0)(0,0,1)[24] w/ mean>
water_fit |>
  augment() |>
  features(.innov, ljung_box, lag = 24, dof = 0)
## # A tibble: 2 × 4
##   Pipe  .model lb_stat lb_pvalue
##   <chr> <chr>    <dbl>     <dbl>
## 1 Pipe1 ARIMA     22.0     0.577
## 2 Pipe2 ARIMA     25.6     0.373
water_fc <- water_fit |> forecast(h = 24 * 7)
water_fc |>
  autoplot(water_filled) +
  facet_wrap(~ Pipe, ncol = 1, scales = "free_x") +
  labs(title = "One-week (168 hour) forecast", y = "WaterFlow", x = NULL)

Both series can be forecast, but there is very little structure to use, so the forecasts stay close to each pipe’s mean:

  • Pipe1: ARIMA picks \((0,0,0)\) with a mean, which is white noise. The forecast is a flat line at about 19.9 for the whole week, and the intervals show the hour-to-hour spread.
  • Pipe2: ARIMA picks a seasonal MA(1) with a 24-hour period. The first day of the forecast has a small daily ripple (about 38 to 43) that carries over from the last 24 hours. After that the forecast settles at the mean of about 39.6.

The residuals of both models pass the Ljung-Box test. Since these series are stationary and close to white noise, the honest one-week forecast is the mean level with wide intervals. Any model that seems to predict more detail than that would be fitting noise.

Saving the forecasts

water_fc_tbl <- water_fc |>
  hilo(level = 95) |>
  unpack_hilo(`95%`) |>
  as_tibble() |>
  transmute(Pipe, Hour = format(Hour, "%Y-%m-%d %H:%M"), Forecast = .mean,
            Lower_95 = `95%_lower`, Upper_95 = `95%_upper`)

write_xlsx(
  list(
    ATM_May2010 = atm_fc |>
      mutate(Forecast_hundreds = round(Forecast, 2),
             Forecast_dollars = round(Forecast * 100, 0)) |>
      select(ATM, DATE, Model, Forecast_hundreds, Forecast_dollars),
    Power_2014 = power_fc_tbl |> mutate(across(where(is.numeric), round)),
    Waterflow_Week = water_fc_tbl |> mutate(across(where(is.numeric), ~ round(.x, 3)))
  ),
  "project1_forecasts.xlsx"
)

The workbook project1_forecasts.xlsx has 3 sheets: ATM_May2010 (31 days for each of the 4 ATMs, in hundreds of dollars and in dollars), Power_2014 (12 months with 95% intervals), and Waterflow_Week (168 hours for each pipe with 95% intervals).