Overview

All three parts follow the same steps from Forecasting: Principles and Practice (fpp3):

  1. Tidy which is loading the Excel file into a tsibble, fixing the dates and dealing with missing values and outliers.
  2. Visualise with a time plot, a seasonal plot and the ACF to see what’s in there.
  3. Model by fitting a few candidates like benchmarks, ETS and ARIMA and comparing them on a hold-out.
  4. Forecast by refitting the winner on all the data and forecasting the horizon that was asked for.

The forecasts also get written to Project1_Forecasts.xlsx with one sheet per part at the end.

Part A – ATM cash withdrawals (May 2010)

Load and clean

# Read the raw file. DATE comes in as an Excel serial number, so convert it.
atm_raw <- read_excel("data/ATM624Data.xlsx")

atm <- atm_raw |>
  mutate(DATE = if (is.numeric(DATE)) as.Date(DATE, origin = "1899-12-30") else as.Date(DATE))

# How many rows, what date range and where are the NAs?
summary(atm)
##       DATE                   ATM            Cash        
##  Min.   :2009-05-01   Length   :1474   Min.   :    0.0  
##  1st Qu.:2009-08-01   N.unique :   4   1st Qu.:    0.5  
##  Median :2009-11-01   N.blank  :   0   Median :   73.0  
##  Mean   :2009-10-31   Min.nchar:   4   Mean   :  155.6  
##  3rd Qu.:2010-02-01   Max.nchar:   4   3rd Qu.:  114.0  
##  Max.   :2010-05-14   NAs      :  14   Max.   :10919.8  
##                                        NAs    :19

Problems found in the raw data:

  • 14 rows with no ATM and no Cash dates 2010‑05‑01 to 2010‑05‑14 are empty placeholders for the month we are asked to forecast so they are dropped.
  • 5 missing Cash values in June 2009 3 for ATM1 and 2 for ATM2.
  • ATM3 is all zeros except the last three days which are Apr 28–30 in 2010.
  • ATM4 has a wild value which is 10,919.76 on 2010‑02‑09, while the next highest day is about 1,712.
# Drop the placeholder rows for May 2010 and build a daily tsibble with one key per ATM.
atm_ts <- atm |>
  filter(!is.na(ATM)) |>
  as_tsibble(index = DATE, key = ATM)

# Treat the  huge ATM4 reading as a data error by  setting it to NA and fill it below.
atm_ts <- atm_ts |>
  mutate(Cash = if_else(ATM == "ATM4" & Cash > 5000, NA_real_, Cash))

# Fill the NAs which are the 5 original and the ATM4 outlier by fitting an ARIMA to each ATM. ATM3 is excluded here because it is nearly all zero.
atm_filled <- atm_ts |>
  filter(ATM != "ATM3") |>
  model(ARIMA(Cash)) |>
  interpolate(atm_ts |> filter(ATM != "ATM3"))

# Put ATM3 back so the full tsibble is complete again.
atm_ts <- bind_rows(atm_filled, atm_ts |> filter(ATM == "ATM3")) |>
  as_tsibble(index = DATE, key = ATM)

# Check that there's no missing values and 365 days per ATM.
atm_ts |> as_tibble() |> group_by(ATM) |> summarise(n = n(), missing = sum(is.na(Cash)))
## # A tibble: 4 × 3
##   ATM       n missing
##   <chr> <int>   <int>
## 1 ATM1    365       0
## 2 ATM2    365       0
## 3 ATM3    365       0
## 4 ATM4    365       0

Visualise

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

# Day-of-week pattern. ATM1 and ATM2 have a clear weekly cycle and notice Thursdays are low.
atm_ts |> filter(ATM %in% c("ATM1", "ATM2", "ATM4")) |>
  gg_season(Cash, period = "week") + facet_wrap(~ATM, scales = "free_y", ncol = 1) +
  labs(title = "Weekly seasonal plot", y = "Cash (x $100)")

atm_ts |> filter(ATM %in% c("ATM1", "ATM2", "ATM4")) |> ACF(Cash, lag_max = 28) |> autoplot() +
  labs(title = "ACF.... the spikes at lags 7, 14, 21 confirm weekly seasonality")

# Averaging cash by weekday for each month,  does the weekly pattern stay the same all year?
atm_ts |> filter(ATM %in% c("ATM1", "ATM2", "ATM4")) |> as_tibble() |>
  mutate(Month = yearmonth(DATE), Weekday = wday(DATE, label = TRUE, week_start = 1)) |>
  group_by(ATM, Month, Weekday) |> summarise(Cash = mean(Cash), .groups = "drop") |>
  # Express each cell relative to that ATM's overall average day so the three machines can share one scale.
  group_by(ATM) |> mutate(rel = Cash / mean(Cash)) |> ungroup() |>
  ggplot(aes(x = Month, y = Weekday, fill = rel)) + geom_tile() +
  facet_wrap(~ATM, ncol = 1) + scale_fill_viridis_c(name = "Cash /\navg day") +
  labs(title = "Cash by weekday and month, relative to each ATM's average day - the low day moves in Feb 2010",
       x = NULL, y = NULL)

What is the plot saying ?

  • ATM1 and ATM2 – Levels are stable with no trend and a strong seasonality weekly(lag‑7 spikes in the ACF).
  • ATM4 – same weekly pattern but much mores noise and more larger amounts . The spread grows with the level, so it would be worth trying a Box-Cox transformation.
  • The low day moved in Feb 2010. Before that the dead day was Thursday and after that it’s Tuesday, and it happened on all three machines at the same time. That looks more like a scheduling change like the refill day moving than customers changing their habits. So the models need to lean on the more recent weeks, and the hold out on April is after the switch so it tests for exactly that.
  • ATM3 is basically empty, there’s only three real days and those values are 96, 82 and 85 which they match ATM1 exactly. Looks like a new machine that is copying ATM1 so far so there’s not enough history to model it on its own.

Models

For ATM1, ATM2 and ATM4 I try three models, and for ATM4 I also try the Box-Cox versions since its variance grows with the level.

Model Why
SNAIVE weekly the benchmark, just repeat last week
ETS auto exponential smoothing with a weekly season
ARIMA auto seasonal ARIMA with period 7
ETS and ARIMA on Box-Cox Cash, ATM4 only same models but on a series with a stable variance

I hold out the last 4 weeks which are Apr 3 to 30 of 2010 to compare how accurate the forecasts are.

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

train <- atm_model |> filter(DATE <= as.Date("2010-04-02"))
test  <- atm_model |> filter(DATE >  as.Date("2010-04-02"))

# The three models for every ATM.
atm_fit <- train |>
  model(
    snaive = SNAIVE(Cash ~ lag("week")),
    ets    = ETS(Cash),
    arima  = ARIMA(Cash)
  )

# Box-Cox lambda for ATM4 using the Guerrero method and then the two transformed models.
lambda4 <- atm_model |> filter(ATM == "ATM4") |> features(Cash, guerrero) |> pull(lambda_guerrero)
lambda4
## [1] 0.4495623
atm4_fit <- train |> filter(ATM == "ATM4") |>
  model(
    ets_bc   = ETS(box_cox(Cash, lambda4)),
    arima_bc = ARIMA(box_cox(Cash, lambda4))
  )

# Accuracy on the 4 week hold out,the forecasts will get back transformed automatically.
atm_acc <- bind_rows(
  atm_fit  |> forecast(new_data = test) |> accuracy(atm_model),
  atm4_fit |> forecast(new_data = test |> filter(ATM == "ATM4")) |> accuracy(atm_model)
) |> select(ATM, .model, RMSE, MAE, MAPE) |> arrange(ATM, RMSE)

atm_acc |> kable(digits = 2, caption = "Hold-out accuracy (Apr 3-30, 2010), lowest RMSE first")
Hold-out accuracy (Apr 3-30, 2010), lowest RMSE first
ATM .model RMSE MAE MAPE
ATM1 ets 12.43 10.04 71.16
ATM1 arima 12.99 10.52 93.68
ATM1 snaive 16.93 14.71 79.16
ATM2 ets 19.86 14.17 55.44
ATM2 arima 21.73 16.89 51.26
ATM2 snaive 28.63 18.64 46.78
ATM4 ets_bc 256.27 215.19 698.92
ATM4 arima_bc 257.29 220.67 920.76
ATM4 snaive 271.89 195.88 364.03
ATM4 arima 281.30 239.82 1115.47
ATM4 ets 310.87 261.43 1329.30 
# Lowest hold out RMSE for each ATM.
atm_acc |> group_by(ATM) |> slice_min(RMSE, n = 1) |> select(ATM, .model) |> kable()
ATM .model
ATM1 ets
ATM2 ets
ATM4 ets_bc

ETS has the lowest hold out error for ATM1 and ATM2 with ARIMA being right behind it, and the Box-Cox ETS wins by a lot for ATM4. Before picking I refited the close ones on all the data and checked the residuals, because if the residuals still have autocorrelation then the prediction intervals come out too narrow.

atm12_fit  <- atm_model |> filter(ATM != "ATM4") |> model(ets = ETS(Cash), arima = ARIMA(Cash))
atm4_final <- atm_model |> filter(ATM == "ATM4") |> model(ets_bc = ETS(box_cox(Cash, lambda4)))

# Ljung-Box on the residuals. p > 0.05 means they look like white noise.
bind_rows(augment(atm12_fit), augment(atm4_final)) |>
  features(.innov, ljung_box, lag = 14) |>
  kable(digits = 3, caption = "Ljung-Box test (lag 14) on residuals, full data fits")
Ljung-Box test (lag 14) on residuals, full data fits
ATM .model lb_stat lb_pvalue
ATM1 arima 15.475 0.346
ATM1 ets 29.435 0.009
ATM2 arima 9.065 0.827
ATM2 ets 35.853 0.001
ATM4 ets_bc 16.127 0.306

The ETS residuals for ATM1 and ATM2 fail the test with p < 0.01 but the ARIMA residuals pass, and ARIMA was only about $100 a day worse on the hold out so I go with ARIMA for ATM1 and ATM2 since the accuracy is almost the same and the intervals can be trusted, and Box-Cox ETS for ATM4.

atm12_final <- atm12_fit |> select(arima)

atm12_final |> filter(ATM == "ATM1") |> report()
## Series: Cash 
## Model: ARIMA(0,0,1)(0,1,2)[7] 
## 
## Coefficients:
##          ma1     sma1     sma2
##       0.1896  -0.5778  -0.1067
## s.e.  0.0547   0.0506   0.0495
## 
## sigma^2 estimated as 553.8:  log likelihood=-1639.17
## AIC=3286.34   AICc=3286.46   BIC=3301.87
atm12_final |> filter(ATM == "ATM2") |> report()
## Series: Cash 
## Model: ARIMA(2,0,2)(0,1,1)[7] 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2     sma1
##       -0.4243  -0.9127  0.4603  0.8005  -0.7488
## s.e.   0.0547   0.0420  0.0849  0.0572   0.0388
## 
## sigma^2 estimated as 586:  log likelihood=-1648.62
## AIC=3309.24   AICc=3309.48   BIC=3332.52
atm4_final  |> report()
## Series: Cash 
## Model: ETS(M,N,A) 
## Transformation: box_cox(Cash, lambda4) 
##   Smoothing parameters:
##     alpha = 0.01429494 
##     gamma = 0.1751993 
## 
##   Initial states:
##      l[0]      s[0]     s[-1]    s[-2]    s[-3]    s[-4]    s[-5]    s[-6]
##  26.57888 -18.86539 -6.412597 2.822399 4.028224 9.683004 4.507669 4.236694
## 
##   sigma^2:  0.2216
## 
##      AIC     AICc      BIC 
## 4020.852 4021.473 4059.851
# Residual plots for two of the final models. Looking for no pattern and the ACF inside the bounds.
atm12_final |> filter(ATM == "ATM1") |> gg_tsresiduals(lag_max = 28)

atm4_final |> gg_tsresiduals(lag_max = 28)

atm_fc <- bind_rows(atm12_final |> forecast(h = 31),
                    atm4_final  |> forecast(h = 31))

atm12_final |> forecast(h = 31) |>
  autoplot(atm_model |> filter(DATE >= as.Date("2010-02-01")), level = 80) +
  facet_wrap(~ATM, scales = "free_y", ncol = 1) +
  labs(title = "May 2010 forecasts - ATM1 and ATM2 (hundreds of dollars)", y = "Cash (x $100)")

atm4_final |> forecast(h = 31) |>
  autoplot(atm_model |> filter(DATE >= as.Date("2010-02-01")), level = 80) +
  labs(title = "May 2010 forecast - ATM4 (hundreds of dollars)", y = "Cash (x $100)")

ATM3

With three days of real history there’s nothing to really fit, since those three days match ATM1 exactly the most reasonable assumption is that ATM3 keeps acting just like ATM1, so I use the ATM1 forecast as a proxy for ATM3 and flag it as low confidence. Once it has a few more weeks of its own data thene it should be reforecasted again

atm3_fc <- atm_fc |> filter(ATM == "ATM1") |> as_tibble() |>
  mutate(ATM = "ATM3", .model = "proxy_ATM1") |>
  select(ATM, DATE, .mean)

Final forecast (Part A)

atm_out <- atm_fc |> as_tibble() |> select(ATM, DATE, .mean) |>
  bind_rows(atm3_fc) |>
  mutate(Cash_hundreds = round(pmax(.mean, 0), 1),
         Cash_dollars  = round(Cash_hundreds * 100)) |>
  select(ATM, DATE, Cash_hundreds, Cash_dollars) |>
  arrange(ATM, DATE)

# Monthly totals which is the number a business person actually wants.
atm_out |> group_by(ATM) |>
  summarise(May_2010_total_dollars = sum(Cash_dollars),
            avg_per_day_dollars    = round(mean(Cash_dollars))) |>
  kable(format.args = list(big.mark = ","), caption = "May 2010 forecast totals")
May 2010 forecast totals
ATM May_2010_total_dollars avg_per_day_dollars
ATM1 242,640 7,827
ATM2 179,070 5,776
ATM3 242,640 7,827
ATM4 1,301,870 41,996 
atm_out |> pivot_wider(names_from = ATM, values_from = Cash_dollars, id_cols = DATE) |>
  kable(format.args = list(big.mark = ","), caption = "Daily forecast, dollars")
Daily forecast, dollars
DATE ATM1 ATM2 ATM3 ATM4
2010-05-01 8,720 6,640 8,720 33,720
2010-05-02 10,410 7,270 10,410 45,630
2010-05-03 7,310 1,360 7,310 46,930
2010-05-04 810 500 810 4,490
2010-05-05 10,000 9,780 10,000 58,250
2010-05-06 8,090 8,870 8,090 38,130
2010-05-07 8,670 6,530 8,670 61,910
2010-05-08 8,800 6,640 8,800 34,130
2010-05-09 10,280 7,270 10,280 46,070
2010-05-10 7,340 1,360 7,340 46,990
2010-05-11 860 500 860 4,910
2010-05-12 10,060 9,780 10,060 58,620
2010-05-13 8,040 8,870 8,040 38,680
2010-05-14 8,690 6,530 8,690 62,250
2010-05-15 8,800 6,640 8,800 34,550
2010-05-16 10,280 7,270 10,280 46,510
2010-05-17 7,340 1,360 7,340 47,050
2010-05-18 860 500 860 5,340
2010-05-19 10,060 9,780 10,060 59,000
2010-05-20 8,040 8,870 8,040 39,240
2010-05-21 8,690 6,530 8,690 62,590
2010-05-22 8,800 6,640 8,800 34,970
2010-05-23 10,280 7,270 10,280 46,950
2010-05-24 7,340 1,360 7,340 47,110
2010-05-25 860 500 860 5,770
2010-05-26 10,060 9,780 10,060 59,380
2010-05-27 8,040 8,870 8,040 39,800
2010-05-28 8,690 6,530 8,690 62,930
2010-05-29 8,800 6,640 8,800 35,400
2010-05-30 10,280 7,270 10,280 47,390
2010-05-31 7,340 1,360 7,340 47,180

Discussion (Part A)

I used ARIMA interpolation for the few missing daysthen seasonal naive, ETS and seasonal ARIMA for each ATM including the Box-Cox versions for ATM4. A 4 week hold out to shortlist and Ljung-Box on the residuals to break the tie between ETS and ARIMA, which is the reason ARIMA got picked for ATM1 and ATM2. I didn’t use regression with calendar or holiday dummies since there’s no holiday info and the seasonal models already catch the weekly pattern. Also no Neural nets or Prophet which is an overkill for 365 points and not in scope. One combined model across ATMs since the machines behave so differently.Some judgement calls had to be made, the ATM4 spike on Feb 9 is treated as a data error and not a real withdrawal day. Leaving it in would inflate every ATM4 forecast and its intervals. ATM3 is a proxy and not a model. The forecasts put the low day on Tuesdays which are May 4, 11, 18 and 25, not Thursdays.That’s on purpose. The weekly pattern moved in February 2010 and the seasonal ARIMA and ETS follow the recent pattern. If operations knows the refill day is about to move again then the forecast should be ran again. ATM4 is the one you to keep an eye out for. Its daily amounts are noisy so the point forecast is really just a weekly pattern average and the intervals are wide. When trying to decide how much cash to load the weekly total is more reliable and safer than any single day.

Part B – Residential power usage (2014 forecast)

Load and clean

power_raw <- read_excel("data/ResidentialCustomerForecastLoad-624.xlsx")

power <- power_raw |>
  mutate(Month = yearmonth(ym(`YYYY-MMM`))) |>
  select(Month, KWH) |>
  as_tsibble(index = Month)

summary(power$KWH)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      NAs 
##   770523  5429912  6283324  6502475  7620524 10655730        1
power |> filter(is.na(KWH) | KWH < 2e6)
## # A tsibble: 2 x 2 [1M]
##      Month    KWH
##      <mth>  <dbl>
## 1 2008 Sep     NA
## 2 2010 Jul 770523

There’s two issues we encountered here where one of them is a missing month which is Sep 2008 and the other one is reading that makes no sense which is Jul 2010 at 770,523 kWh, which is about about a tenth of a normal July. I treat the bad July as missing too and fill both with ARIMA interpolation.

power <- power |> mutate(KWH = if_else(KWH < 2e6, NA_real_, KWH))
power <- power |> model(ARIMA(KWH)) |> interpolate(power)
sum(is.na(power$KWH))
## [1] 0

Visualise

power |> autoplot(KWH / 1e6) + labs(title = "Monthly residential power usage", y = "GWh")

power |> gg_season(KWH / 1e6) + labs(title = "Seasonal plot - summer and winter peaks", y = "GWh")

power |> gg_subseries(KWH / 1e6) + labs(y = "GWh")

power |> ACF(KWH, lag_max = 36) |> autoplot() + labs(title = "ACF - strong annual seasonality")

The plot says there’s clear annual seasonality with peaks in Jan/Feb and Jul/Aug and during dips in spring and fall. The level is mostly flat from 1998 to 2010 and then steps up in 2011 by around 15 to 20%. The seasonal swings also got bigger after 2011 so trying the Box-Cox transformation is worth giving it a try.

Models

The candidates are seasonal naive, ETS, ARIMA and ARIMA with Box-Cox. The last 24 months which is 2012 and 2013 are the hold-out.

lambda_p <- power |> features(KWH, features = guerrero) |> pull(lambda_guerrero)
lambda_p
## [1] -0.2057366
p_train <- power |> filter(Month <= yearmonth("2011 Dec"))
p_test  <- power |> filter(Month >  yearmonth("2011 Dec"))

power_fit <- p_train |>
  model(
    snaive   = SNAIVE(KWH),
    ets      = ETS(KWH),
    arima    = ARIMA(KWH),
    arima_bc = ARIMA(box_cox(KWH, lambda_p))
  )

power_fit |> forecast(new_data = p_test) |> accuracy(power) |>
  select(.model, RMSE, MAE, MAPE) |> arrange(RMSE) |>
  kable(digits = 2, caption = "Hold-out accuracy (2012-2013)")
Hold-out accuracy (2012-2013)
.model RMSE MAE MAPE
arima 721232.5 532508.3 6.81
arima_bc 773435.0 617048.8 7.97
ets 850432.8 622085.5 8.41
snaive 965567.1 742063.5 9.30 
best_p <- power_fit |> forecast(new_data = p_test) |> accuracy(power) |>
  slice_min(RMSE, n = 1) |> pull(.model)
best_p
## [1] "arima"
# Refit on all the data with the winner.
power_full <- power |>
  model(
    snaive   = SNAIVE(KWH),
    ets      = ETS(KWH),
    arima    = ARIMA(KWH),
    arima_bc = ARIMA(box_cox(KWH, lambda_p))
  ) |> select(all_of(best_p))

report(power_full)
## 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_full |> gg_tsresiduals(lag_max = 36)

augment(power_full) |> features(.innov, ljung_box, lag = 24)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 arima     14.2     0.943

Forecast (Part B)

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

power_fc |> autoplot(power |> filter(Month >= yearmonth("2008 Jan")), level = c(80, 95)) +
  labs(title = "2014 monthly forecast", y = "KWH")

power_out <- power_fc |> hilo(95) |> unpack_hilo(`95%`) |> as_tibble() |>
  transmute(Month = as.character(Month),
            KWH_forecast = round(.mean),
            lower_95 = round(`95%_lower`), upper_95 = round(`95%_upper`))

power_out |> kable(format.args = list(big.mark = ","), caption = "2014 forecast with 95% interval")
2014 forecast with 95% interval
Month KWH_forecast lower_95 upper_95
2014 Jan 10,295,036 9,088,564 11,501,508
2014 Feb 8,617,878 7,345,438 9,890,318
2014 Mar 7,164,745 5,891,422 8,438,067
2014 Apr 5,831,199 4,530,511 7,131,888
2014 May 5,942,117 4,638,637 7,245,598
2014 Jun 8,199,229 6,895,748 9,502,710
2014 Jul 9,478,936 8,175,455 10,782,416
2014 Aug 9,973,049 8,669,569 11,276,530
2014 Sep 8,486,866 7,183,385 9,790,346
2014 Oct 5,870,156 4,566,676 7,173,637
2014 Nov 6,147,210 4,843,730 7,450,691
2014 Dec 8,421,706 7,118,225 9,725,186

Discussion (Part B)

What was used here was ARIMA interpolation for the missing month and the bad July 2010 reading, then SNAIVE, ETS, ARIMA and Box-Cox ARIMA with a 24 month hold out and residual checks. I didn’t use dynamic regression on temperature or customer counts since none of that was given. That would be the obvious next step since the usage is driven by weather. STL decomposition models since the seasonal ARIMA already handles the seasonality and the seasonal difference takes care of the 2011 step. The winner here is a seasonal ARIMA with drift and its residuals pass Ljung-Box. The 2011 jump is the main risk because the model assumes the post 2011 level sticks around, which is fine for a one year horizon.

Part C (bonus) – Water flow, two pipes

Load, time-sequence and aggregate by hour

read_pipe <- function(path, name) {
  read_excel(path) |>
    rename(dt = `Date Time`) |>
    mutate(dt = if (is.numeric(dt)) as.POSIXct(dt * 86400, origin = "1899-12-30", tz = "UTC")
                else as.POSIXct(dt, tz = "UTC"),
           Pipe = name)
}
pipe1_raw <- read_pipe("data/Waterflow_Pipe1.xlsx", "Pipe1")
pipe2_raw <- read_pipe("data/Waterflow_Pipe2.xlsx", "Pipe2")

# Pipe 1 has several readings per hour at random times so floor to the hour.
# Pipe 2 has exactly one reading per hour but the Excel serials have non rounded numbers
#         like 01:59:59.9999 instead of 02:00:00, so round to the nearest hour.
pipe1_raw |> summarise(from = min(dt), to = max(dt), n = n(),
                       median_gap_min = median(as.numeric(diff(dt), units = "mins")))
## # A tibble: 1 × 4
##   from                to                      n median_gap_min
##   <dttm>              <dttm>              <int>          <dbl>
## 1 2015-10-23 00:24:06 2015-11-01 23:35:43  1000           10.0
pipe2_raw |> summarise(from = min(dt), to = max(dt), n = n(),
                       median_gap_min = median(as.numeric(diff(dt), units = "mins")))
## # A tibble: 1 × 4
##   from                to                      n median_gap_min
##   <dttm>              <dttm>              <int>          <dbl>
## 1 2015-10-23 01:00:00 2015-12-03 16:00:00  1000           60.0
pipe1_h <- pipe1_raw |> mutate(hour = floor_date(dt, "hour")) |>
  group_by(hour) |> summarise(WaterFlow = mean(WaterFlow), n_readings = n()) |>
  mutate(Pipe = "Pipe1")

pipe2_h <- pipe2_raw |> mutate(hour = round_date(dt, "hour")) |>
  group_by(hour) |> summarise(WaterFlow = mean(WaterFlow), n_readings = n()) |>
  mutate(Pipe = "Pipe2")

# Both pipes lined up by the hour, which was the  layout the prompt wanted.
pipes_wide <- full_join(pipe1_h |> select(hour, Pipe1 = WaterFlow),
                        pipe2_h |> select(hour, Pipe2 = WaterFlow), by = "hour") |>
  arrange(hour)
head(pipes_wide, 10) |> kable(digits = 2, caption = "Hourly means, both pipes (first 10 hours)")
Hourly means, both pipes (first 10 hours)
hour Pipe1 Pipe2
2015-10-23 00:00:00 26.10 NA
2015-10-23 01:00:00 18.85 18.81
2015-10-23 02:00:00 15.16 43.09
2015-10-23 03:00:00 23.08 37.99
2015-10-23 04:00:00 15.48 36.12
2015-10-23 05:00:00 22.73 31.85
2015-10-23 06:00:00 20.59 28.24
2015-10-23 07:00:00 18.37 9.86
2015-10-23 08:00:00 20.79 26.68
2015-10-23 09:00:00 21.22 55.77 
# One hourly tsibble keyed by pipe. Fill the few empty hours in Pipe 1 by interpolation.
pipes <- bind_rows(pipe1_h, pipe2_h) |>
  select(Pipe, hour, WaterFlow) |>
  as_tsibble(index = hour, key = Pipe) |>
  fill_gaps()

pipes |> as_tibble() |> group_by(Pipe) |>
  summarise(hours = n(), missing_hours = sum(is.na(WaterFlow)),
            from = min(hour), to = max(hour))
## # A tibble: 2 × 5
##   Pipe  hours missing_hours from                to                 
##   <chr> <int>         <int> <dttm>              <dttm>             
## 1 Pipe1   240             4 2015-10-23 00:00:00 2015-11-01 23:00:00
## 2 Pipe2  1000             0 2015-10-23 01:00:00 2015-12-03 16:00:00
pipes <- pipes |> model(ARIMA(WaterFlow)) |> interpolate(pipes)

Pipe 1 covers about 10 days or 240 hours and Pipe 2 about 41 days or 1,000 hours so they only overlap for the first 10 days. Adding them into one total flow series would make a fake drop when Pipe 1 ends so each pipe gets its own model.

Is it stationary?

pipes |> autoplot(WaterFlow) + facet_wrap(~Pipe, scales = "free", ncol = 1) +
  labs(title = "Hourly mean water flow", y = "Flow")

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

# KPSS test. with p-value < 0.05  the means aren't stationary and it needs differencing.
pipes |> features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs, unitroot_nsdiffs)) |>
  kable(digits = 3, caption = "KPSS test and number of differences suggested")
KPSS test and number of differences suggested
Pipe kpss_stat kpss_pvalue ndiffs nsdiffs
Pipe1 0.079 0.1 0 0
Pipe2 0.105 0.1 0 0

Pipe 1 is flat around 20 with no trend and the ACF is inside the bounds, so it’s stationary and pretty much white noise while pipe 2 is bigger and bouncier around 40 but also no trend. The KPSS p-value is 0.1 for both pipes and ndiffs is 0 so no differencing needed and both can be forecast. There’s just very little structure to work with so the forecast is mostly the level with honest intervals.

Models and one week forecast (168 hours)

# Hold out the last 48 hours of each pipe to compare a few simple models.
pipe_train <- pipes |> group_by_key() |> slice_head(n = -48) |> ungroup()
pipe_test  <- pipes |> group_by_key() |> slice_tail(n = 48)  |> ungroup()

pipe_fit <- pipe_train |>
  model(
    mean  = MEAN(WaterFlow),
    naive = NAIVE(WaterFlow),
    ets   = ETS(WaterFlow),
    arima = ARIMA(WaterFlow)
  )

pipe_fit |> forecast(new_data = pipe_test) |> accuracy(pipes) |>
  select(Pipe, .model, RMSE, MAE) |> arrange(Pipe, RMSE) |>
  kable(digits = 2, caption = "48-hour hold-out accuracy")
48-hour hold-out accuracy
Pipe .model RMSE MAE
Pipe1 ets 6.41 3.97
Pipe1 arima 6.41 3.97
Pipe1 mean 6.41 3.97
Pipe1 naive 7.41 5.78
Pipe2 ets 17.10 13.71
Pipe2 mean 17.10 13.71
Pipe2 arima 17.30 13.85
Pipe2 naive 19.99 16.89 
best_pipe <- pipe_fit |> forecast(new_data = pipe_test) |> accuracy(pipes) |>
  group_by(Pipe) |> slice_min(RMSE, n = 1) |> select(Pipe, .model)
best_pipe
## # A tibble: 2 × 2
## # Groups:   Pipe [2]
##   Pipe  .model
##   <chr> <chr> 
## 1 Pipe1 ets   
## 2 Pipe2 ets
pipe_full <- pipes |>
  model(
    mean  = MEAN(WaterFlow),
    naive = NAIVE(WaterFlow),
    ets   = ETS(WaterFlow),
    arima = ARIMA(WaterFlow)
  ) |>
  pivot_longer(-Pipe, names_to = ".model", values_to = "fit") |>
  semi_join(best_pipe, by = c("Pipe", ".model"))

pipe_full |> select(-.model) |> report()
## # A tibble: 2 × 10
##   Pipe  .model sigma2 log_lik    AIC   AICc    BIC   MSE  AMSE   MAE
##   <chr> <chr>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 Pipe1 fit    28.5    -1059.  2123.  2123.  2134.  28.2  28.2 3.70 
## 2 Pipe2 fit     0.165  -6229. 12464. 12464. 12479. 257.  257.  0.331
pipe_fc <- pipe_full |> select(-.model) |> forecast(h = 168)

pipe_fc |> autoplot(pipes, level = 80) + facet_wrap(~Pipe, scales = "free", ncol = 1) +
  labs(title = "One-week-ahead hourly forecast", y = "Flow")

pipe_out <- pipe_fc |> hilo(80) |> unpack_hilo(`80%`) |> as_tibble() |>
  transmute(Pipe, hour, WaterFlow_forecast = round(.mean, 2),
            lower_80 = round(`80%_lower`, 2), upper_80 = round(`80%_upper`, 2))
head(pipe_out, 12) |> kable(caption = "First 12 hours of the Pipe 1 forecast")
First 12 hours of the Pipe 1 forecast
Pipe hour WaterFlow_forecast lower_80 upper_80
Pipe1 2015-11-02 00:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 01:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 02:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 03:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 04:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 05:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 06:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 07:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 08:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 09:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 10:00:00 19.47 12.63 26.3
Pipe1 2015-11-02 11:00:00 19.47 12.63 26.3

Discussion (Part C)

Both series are stationary since KPSS doesn’t reject and no differencing is suggested, so they could be forecast. There’s just not much to grab onto and no daily cycle in the ACF. The winner for both pipes is a simple exponential smoothing ETS(A,N,N) with a tiny smoothing parameter which is basically the mean, and it ties the mean on the hold out. So the week ahead forecast is a flat line around 19.5 for Pipe 1 and 39.5 for Pipe 2 with intervals that show the hourly noise.

Excel output

write_xlsx(
  list(
    ATM_May2010_daily   = atm_out,
    ATM_May2010_totals  = atm_out |> group_by(ATM) |>
                            summarise(May_2010_total_dollars = sum(Cash_dollars)),
    Power_2014_monthly  = power_out,
    Pipe1_week_hourly   = pipe_out |> filter(Pipe == "Pipe1") |> mutate(hour = format(hour, "%Y-%m-%d %H:%M")),
    Pipe2_week_hourly   = pipe_out |> filter(Pipe == "Pipe2") |> mutate(hour = format(hour, "%Y-%m-%d %H:%M"))
  ),
  "Project1_Forecasts.xlsx"
)

The forecast file Project1_Forecasts.xlsx goes in with this report.