DATA 624 Project 1: Forecasting ATM Cash, Residential Power and Waterflow

library(fpp3)
library(readxl)
library(writexl)

Part A: ATM Forecast

1. Reading and preparing the data

The goal of Part A is to forecast how much cash will be withdrawn from four ATMs during May 2010. The variable Cash is in hundreds of dollars. The following is going to be a process to read and clean thd data then deal with missing values and outliers and look at the patterns of each machine. Later compare some candidate models with time series cross validation and check the residuals of the winers.

atm_raw <- read_excel("data/ATM624Data.xlsx")
glimpse(atm_raw)

## 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, …

Excel always stores dates as a count of days. Using the origin 1899-12-30 turns it back into a calendar date, for example 39934 becomes 2009-05-01.

atm_raw |>
  summarise(rows = n(),
            no_atm_label = sum(is.na(ATM)),
            first_empty_date = as.Date(min(DATE[is.na(ATM)]), origin = "1899-12-30"),
            last_empty_date = as.Date(max(DATE[is.na(ATM)]), origin = "1899-12-30"))

atm <- atm_raw |>
  filter(!is.na(ATM)) |>
  mutate(DATE = as.Date(DATE, origin = "1899-12-30")) |>
  as_tsibble(key = ATM, index = DATE)

atm |> as_tibble() |> group_by(ATM) |>
  summarise(first = min(DATE), last = max(DATE), days = n(),
            missing = sum(is.na(Cash)), zeros = sum(Cash == 0, na.rm = TRUE),
            median = median(Cash, na.rm = TRUE), max = max(Cash, na.rm = TRUE))

The file has 14 rows from 2010-05-01 to 2010-05-14 that have a date but no atm label and no cash value. This are empty placeholder rows, not observations, so I drop them. That leaves 365 days for each ATM from 2009-05-01 to 2010-04-30 with no gaps in the calendar.

atm |>
  autoplot(Cash) +
  facet_grid(ATM ~ ., scales = "free_y") +
  labs(title = "Daily cash withdrawn by ATM", y = "Hundreds of dollars", x = NULL) +
  theme(legend.position = "none")

ATM1 and ATM2 move between 0 and ~150 with a repeating up and down pattern. ATM4 sits mostly below 1500 but one spike near 10000 flattens everything else on its panel. ATM3 is a flat line at 0 that only moves at the very end. Each machine needs its own treatment, so I handle them one at a time.

2. ATM3 has almost no history

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

ATM3 recorded cash on only the last three days (96, 82 and 85). Before that it is zero every day, which looks like a machine that was not working rather than filly working. 3 observations cannot support a trend, a weekly pattern or any model that has parameters to estimate, so ATM3 is forecast separately with a simple level later on. The remaining work is on ATM1, ATM2 and ATM4.

3. Missing values

atm_main <- atm |> filter(ATM != "ATM3")

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

There are 5 missing Cash values, all in June 2009 (3 for ATM1 and 2 for ATM2). Deleting those rows would break the daily sequence that the seasonal models need. Filling them with the overall mean would ignore the weekly cycle, and Thursdays are very low, so a mean would be badly wrong on some days. Instead I’m going to use the ARIMA interpolation from the fpp3 book. It fits a seasonal ARIMA to the series and fills each gap with the value that best agrees with the surrounding days and the weekly pattern.

atm_i <- atm_main |>
  model(ARIMA(Cash)) |>
  interpolate(atm_main)

weekday_median <- atm_main |> as_tibble() |>
  mutate(weekday = wday(DATE, label = TRUE, week_start = 1)) |>
  group_by(ATM, weekday) |>
  summarise(weekday_median = median(Cash, na.rm = TRUE), .groups = "drop")

atm_main |> filter(is.na(Cash)) |> as_tibble() |> select(ATM, DATE) |>
  left_join(atm_i |> as_tibble() |> rename(imputed = Cash), by = c("ATM", "DATE")) |>
  mutate(weekday = wday(DATE, label = TRUE, week_start = 1)) |>
  left_join(weekday_median, by = c("ATM", "weekday"))

Each imputed value is close to the median for the same weekday. For example, ATM1 on Saturday 06/13 was filled with 111 against a Saturday median of 96, and ATM2 on Thursday 06/18 was filled with 9.4 against a Thursday median of 8. This tells me the interpolation respects the weekly cycle.

4. Outliers

stl_parts <- atm_i |>
  model(STL(Cash ~ season(7), robust = TRUE)) |>
  components() |>
  as_tibble() |>
  group_by(ATM) |>
  mutate(iqr = IQR(remainder),
         q1 = quantile(remainder, 0.25),
         q3 = quantile(remainder, 0.75)) |>
  ungroup()

flag_count <- function(k) {
  stl_parts |>
    filter(remainder < q1 - k * iqr | remainder > q3 + k * iqr) |>
    count(ATM, name = paste0("beyond_", k, "_IQR"))
}

flag_count(3) |>
  full_join(flag_count(10), by = "ATM") |>
  full_join(flag_count(20), by = "ATM") |>
  mutate(across(-ATM, ~ replace_na(.x, 0L)))

For outliers I decomposed each series with a robust STL using a 7 day season. The remainder is what is left after the trend and the weekly pattern are removed, so an outlier is a day that is unusual for its weekday. I flag points whose remainder is far outside the interquartile range (IQR). At the common cutoff of 3 IQR, ATM1 has 42 flagged days and ATM2 has 33. That is too many to be errors. Most are unusually low days that appear throughout the year, which tells me these two machines simply have heavy tailed daily noise. Even at 10 IQR there are still 5 flagged days for ATM1. Only one point survives a 20 IQR cutoff, and it belongs to ATM4.

extreme <- stl_parts |>
  filter(remainder < q1 - 20 * iqr | remainder > q3 + 20 * iqr) |>
  select(ATM, DATE, Cash, remainder)
extreme

atm_i |>
  filter(ATM == "ATM4") |>
  autoplot(Cash) +
  geom_point(data = extreme, aes(DATE, Cash), colour = "red", size = 3) +
  labs(title = "ATM4",
       y = "Hundreds of dollars", x = NULL)

The flagged value is ATM4 at ~10000. That is roughly 27 times the ATM4 median of 404, it sits 35 IQR from its expected value, and the next largest value in the whole series is ~1800. A single day this extreme is almost certainly a recording error or a one time event that will not repeat in May, and if I keep it every variance estimate and prediction interval for ATM4 would be inflated. If I replace only this point. I deliberately keep the low days in ATM1 and ATM2, because they happen repeatedly and are part of how those machines behave.

atm_na <- atm_i |>
  mutate(Cash = if_else(paste(ATM, DATE) %in% paste(extreme$ATM, extreme$DATE),
                        NA_real_, Cash))

atm_clean <- atm_na |>
  model(ARIMA(Cash)) |>
  interpolate(atm_na)

atm_clean |>
  filter(paste(ATM, DATE) %in% paste(extreme$ATM, extreme$DATE)) |>
  as_tibble() |>
  mutate(weekday = wday(DATE, label = TRUE, week_start = 1)) |>
  left_join(weekday_median, by = c("ATM", "weekday"))

atm_clean |>
  autoplot(Cash) +
  facet_grid(ATM ~ ., scales = "free_y") +
  labs(title = "Cleaned daily cash", y = "Hundreds of dollars", x = NULL) +
  theme(legend.position = "none")

The interpolated value for ATM4 on 2010-02-09 is about 214, a little below the Tuesday median of 307 but well within the normal range of the machine. After cleaning, ATM4 is on a readable scale, with values from near 0 to about 1700.

5. Patterns in the cleaned data

atm_clean |>
  gg_season(Cash, period = "week") +
  facet_grid(ATM ~ ., scales = "free_y") +
  labs(title = "Weekly pattern by ATM", y = "Hundreds of dollars")

All three machines have a strong weekly pattern. Thursday is the lowest day for every ATM, and it is close to zero for ATM1 and ATM2. In some weeks, mostly earlier in the year, Monday or Tuesday was also near zero for ATM1 and ATM2, so the weekly shape is not perfectly stable. ATM4 is much noisier from week to week, but its Thursday dip is still visible.

atm_clean |>
  ACF(Cash, lag_max = 28) |>
  autoplot() +
  labs(title = "ACF of cleaned daily cash")

The ACF confirms a 7 day cycle. ATM1 and ATM2 have large spikes at lags 7, 14, 21 and 28, around 0.6 to 0.7, and ATM2 also has a repeating negative pattern between them. ATM4 has the same shape but a much weaker spike (about 0.15 at lag 7)

6. Choosing a transformation

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

The Guerrero method suggests lambda values of about 0.26 for ATM1, 0.67 for ATM2 and 0.45 for ATM4. These are all between a log (0) and no transformation, so the variance does grow somewhat with the level but not dramatically. Because the evidence is not overwhelming I’m not going to force a transormation. I’ll include both a Box Cox version and an untransformed version of the ETS and ARIMA models and let cross validation decide.

7. Candidate models and cross validation

In the following steps, I’ll compare six candidate models for each of ATM1, ATM2 and ATM4. - SNAIVE repeats the last observed week and is the benchmark any real model should beat. - ETS and ARIMA are fit on the raw data and again on the Box Cox scale (the _bc versions). - STL_ETS decomposes the Box Cox series with a weekly STL and forecasts the seasonally adjusted part with ETS.

To judge them I’ll use time series cross validation, training windows start at 200 days and grow by 28 days, and each window forecasts the next 31 days, the same length as May. The models will be ranked by RMSE. MAPE is not applicable because ATM2 has days with zero cash, which makes MAPE infinite, and AICc cannot compare transformed and untransformed fits.

fit_set <- function(data, lambda) {
  data |>
    model(
      SNAIVE   = SNAIVE(Cash),
      ETS      = ETS(Cash),
      ETS_bc   = ETS(box_cox(Cash, lambda)),
      ARIMA    = ARIMA(Cash),
      ARIMA_bc = ARIMA(box_cox(Cash, lambda)),
      STL_ETS  = decomposition_model(
        STL(box_cox(Cash, lambda) ~ season(7), robust = TRUE),
        ETS(season_adjust ~ season("N")))
    )
}

atm_names <- c("ATM1", "ATM2", "ATM4")
cv_acc <- list()

for (a in atm_names) {
  d <- atm_clean |> filter(ATM == a)
  lam <- lambdas |> filter(ATM == a) |> pull(lambda_guerrero)
  cv_acc[[a]] <- d |>
    stretch_tsibble(.init = 200, .step = 28) |>
    fit_set(lam) |>
    forecast(h = 31) |>
    accuracy(d) |>
    mutate(ATM = a)
}

cv_table <- bind_rows(cv_acc) |>
  select(ATM, .model, RMSE, MAE) |>
  arrange(ATM, RMSE)
cv_table

The Box Cox ETS has the lowest RMSE for ATM1 and ATM2, and the Box Cox ARIMA is best for ATM4. The best model beats SNAIVE, STL_ETS is the weakest model for ATM1 and ATM2 and does worse than SNAIVE there. Since RMSE alone does not show whether a model captured all the structure, I’ll run a Ljung Box test on the residuals of each model fit to the full data. A p value <0.05 means leftover autocorrelation, so the intervals cannot be trusted. The selection rule is the lowest cross validated RMSE among the models that pass.

lb_list <- list()
for (a in atm_names) {
  d <- atm_clean |> filter(ATM == a)
  lam <- lambdas |> filter(ATM == a) |> pull(lambda_guerrero)
  lb_list[[a]] <- fit_set(d, lam) |>
    augment() |>
    features(.innov, ljung_box, lag = 14, dof = 0) |>
    mutate(ATM = a)
}

cv_lb <- cv_table |>
  left_join(bind_rows(lb_list) |> select(ATM, .model, lb_pvalue), by = c("ATM", ".model"))
cv_lb

best <- cv_lb |>
  filter(lb_pvalue > 0.05) |>
  group_by(ATM) |>
  slice_min(RMSE, n = 1) |>
  ungroup()
best

The Ljung Box filter changed one decision. For ATM2 the ETS_bc had the lowest RMSE but failed the test (p of ~0.003), so I decided to use the ARIMA_bc, which is only about 4% worse on RMSE and passes comfortably.

8. Final models, diagnostics and forecasts

final_fc <- list()
final_fit <- list()

for (a in atm_names) {
  d <- atm_clean |> filter(ATM == a)
  lam <- lambdas |> filter(ATM == a) |> pull(lambda_guerrero)
  winner <- best |> filter(ATM == a) |> pull(.model)
  fit_all <- fit_set(d, lam)
  final_fit[[a]] <- fit_all |> select(all_of(winner))
  final_fc[[a]] <- final_fit[[a]] |>
    forecast(h = 31) |>
    as_tibble() |>
    mutate(ATM = a)
}

final_fit[["ATM1"]] |> gg_tsresiduals(lag_max = 28)

final_fit[["ATM1"]] |> augment() |> features(.innov, ljung_box, lag = 14, dof = 0)

No significant autocorrelation is left (p ~0.15). The large negative residuals are the low days I decided to keep, which make the tails heavier than normal, so the ATM1 intervals should be read as approximate.

final_fit[["ATM2"]] |> gg_tsresiduals(lag_max = 28)

final_fit[["ATM2"]] |> augment() |> features(.innov, ljung_box, lag = 14, dof = 0)

The residuals behave like white noise (p ~0.76), with only a small spike near lag 5.

final_fit[["ATM4"]] |> gg_tsresiduals(lag_max = 28)

final_fit[["ATM4"]] |> augment() |> features(.innov, ljung_box, lag = 14, dof = 0)

The residuals pass the test (p ~0.27). The wide histogram reflects how noisy daily cash at ATM4 is rather than a missing pattern.

fc_all <- bind_rows(final_fc) |>
  as_fable(index = DATE, key = c(ATM, .model), distribution = Cash, response = "Cash")

fc_all |>
  autoplot(atm_clean |> filter(ATM != "ATM3", DATE >= as.Date("2010-02-01")),
           level = c(80, 95)) +
  facet_grid(ATM ~ ., scales = "free_y") +
  labs(title = "May 2010 forecasts", y = "Hundreds of dollars", x = NULL)

ATM1 and ATM2 keep their weekly shape, including the low Thursdays. ATM4 is nearly flat because its weekly signal is weak compared with its daily noise. Some lower bounds fall below zero, which is impossible for cash, so I floored all values at zero in the Excel file.

atm3_level <- atm |> filter(ATM == "ATM3", Cash > 0) |> pull(Cash) |> mean()
atm3_level

## [1] 87.66667

Since ATM3 has no usable history, I decided to forecast it with a flat level equal to the mean of its three active days, which is ~87.7 per day. This is a placeholder

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

fc_main <- fc_all |>
  hilo(level = 80) |>
  as_tibble() |>
  transmute(DATE, ATM,
            Cash = pmax(0, round(.mean, 1)),
            Lower80 = pmax(0, round(`80%`$lower, 1)),
            Upper80 = pmax(0, round(`80%`$upper, 1)))

fc_atm3 <- tibble(DATE = may_dates, ATM = "ATM3",
                  Cash = round(atm3_level, 1), Lower80 = NA_real_, Upper80 = NA_real_)

atm_forecast <- bind_rows(fc_main, fc_atm3) |>
  mutate(DATE = format(DATE, "%Y-%m-%d")) |>
  arrange(ATM, DATE)

write_xlsx(
  list(Forecast = atm_forecast |> select(DATE, ATM, Cash),
       Forecast_with_intervals = atm_forecast),
  "ATM_Forecast_May2010.xlsx")

atm_forecast |> group_by(ATM) |>
  summarise(days = n(), May_total = sum(Cash), daily_mean = mean(Cash))

The workbook ATM_Forecast_May2010.xlsx has the daily point forecast on the first sheet and the 80% interval on the second. ATM4 accounts for about 64% of all cash forecast for May.

Techniques used and not used. I used interpolation for missing values, a robust STL screen for outliers, Box Cox transformations, benchmark, ETS, ARIMA and decomposition models, time series cross validation and Ljung Box tests. I did not use deletion or mean imputation because they ignore the weekly cycle. Also I did not use annual seasonality because the data covers only one year and I did not use regression with outside predictors such as paydays or holidays because no such data was provided, though they would be worth adding in a real business setting.

Part B: Residential Power Forecast

1. Reading and preparing the data

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

## 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 <- power_raw |>
  mutate(Month = yearmonth(lubridate::ym(`YYYY-MMM`))) |>
  select(Month, KWH) |>
  as_tsibble(index = Month)

power |> summarise(first = min(Month), last = max(Month), months = n(),
                   missing = sum(is.na(KWH)))

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

The month is stored as text such as 1998-Jan, so I converted it into a monthly index with lubridate::ym(). There is one missing value, in September 2008.

power |>
  autoplot(KWH) +
  labs(title = "Monthly residential power use", y = "KWH", x = NULL)

The series has a yearly cycle with peaks in summer and winter, and the swings get larger after 2010. Two values need attention. September 2008 is missing and July 2010 drops to about 770 thousand KWH, roughly a tenth of a normal July.

2. Missing value

power_i <- power |>
  model(ARIMA(KWH)) |>
  interpolate(power)

power_i |> filter(Month == yearmonth("2008 Sep"))

power |> as_tibble() |>
  filter(month(Month) == 9) |>
  select(Month, KWH)

The ARIMA interpolation filled September 2008 with about 7.6 million KWH, which is within the range of the other Septembers.

3. Outliers

power_stl <- power_i |>
  model(STL(KWH ~ season(window = "periodic"), robust = TRUE)) |>
  components()

b_iqr <- IQR(power_stl$remainder)
b_q <- quantile(power_stl$remainder, c(0.25, 0.75))

power_flag <- power_stl |>
  filter(remainder < b_q[1] - 3 * b_iqr | remainder > b_q[2] + 3 * b_iqr) |>
  select(Month, KWH, remainder)
power_flag

power_i |> as_tibble() |>
  filter(month(Month) == 7) |>
  select(Month, KWH)

power_i |> as_tibble() |>
  filter(month(Month) == 12) |>
  select(Month, KWH)

The robust STL screen flagged July 2010 and December 2013. July 2010 is clearly a data error, so I decided to replace it. December 2013 is the highest December on record, but winter peaks have been rising, so I decided to keep it as the most recent evidence of the upward drift.

power_na <- power_i |>
  mutate(KWH = if_else(Month == yearmonth("2010 Jul"), NA_real_, KWH))

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

power_clean |> filter(Month == yearmonth("2010 Jul"))

power_clean |>
  autoplot(KWH) +
  labs(title = "Cleaned monthly residential power use", y = "KWH", x = NULL)

The replacement for July 2010 is about 7.7 million, in line with the other Julys.

4. Patterns

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

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

power_clean |> ACF(KWH, lag_max = 36) |> autoplot()

Power use peaks twice a year, in January and in July and August, with lows in spring and October. The ACF peaks at lag 12 and repeats at 24 and 36, confirming yearly seasonality.

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

## [1] -0.2049062

The Guerrero lambda of about negative 0.20 is close to a log, which reflects seasonal swings growing with the level. As in Part A, I’ll compare Box Cox and untransformed versions.

5. Candidate models and cross validation

fit_power <- function(data, lambda) {
  data |>
    model(
      SNAIVE   = SNAIVE(KWH),
      ETS      = ETS(KWH),
      ETS_bc   = ETS(box_cox(KWH, lambda)),
      ARIMA    = ARIMA(KWH),
      ARIMA_bc = ARIMA(box_cox(KWH, lambda)),
      STL_ETS  = decomposition_model(
        STL(box_cox(KWH, lambda) ~ season(window = "periodic"), robust = TRUE),
        ETS(season_adjust ~ season("N")))
    )
}

power_cv <- power_clean |>
  stretch_tsibble(.init = 120, .step = 12) |>
  fit_power(power_lambda) |>
  forecast(h = 12) |>
  accuracy(power_clean) |>
  select(.model, RMSE, MAE, MAPE) |>
  arrange(RMSE)
power_cv

The same six candidates from Part A were compared with a 12 month season. Training windows start at 120 months, grow by 12, and forecast the next 12 months to match the 2014 task. The Box Cox ARIMA has the lowest RMSE and improves on SNAIVE by about 18%, with the untransformed ARIMA close behind.

6. Final model and forecast for 2014

power_best <- power_cv |> slice_min(RMSE, n = 1) |> pull(.model)
power_best

## [1] "ARIMA_bc"

power_fit <- fit_power(power_clean, power_lambda) |> select(all_of(power_best))
power_fit |> report()

## Series: KWH 
## Model: ARIMA(0,0,1)(2,1,0)[12] w/ drift 
## Transformation: box_cox(KWH, lambda) 
## 
## Coefficients:
##          ma1     sar1     sar2  constant
##       0.2462  -0.7202  -0.3546    0.0013
## s.e.  0.0833   0.0751   0.0769    0.0004
## 
## sigma^2 estimated as 1.505e-05:  log likelihood=751.74
## AIC=-1493.47   AICc=-1493.13   BIC=-1477.51

power_fit |> gg_tsresiduals(lag_max = 36)

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

The selected model is an ARIMA(0,0,1)(2,1,0)[12] with drift on the Box Cox scale. The Ljung Box p value of ~0.14 shows no leftover autocorrelation overall. One spike near lag 3 sits just outside the limits, but a single lag is not a concern when the overall test passes. The flat first year of residuals comes from the seasonal difference.

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

power_fc |>
  autoplot(power_clean |> filter(year(Month) >= 2008), level = c(80, 95)) +
  labs(title = "Residential power forecast for 2014", y = "KWH", x = NULL)

power_out <- power_fc |>
  hilo(level = 80) |>
  as_tibble() |>
  transmute(Month = format(Month, "%Y-%b"),
            KWH = round(.mean),
            Lower80 = round(`80%`$lower),
            Upper80 = round(`80%`$upper))

write_xlsx(list(Forecast = power_out |> select(Month, KWH),
                Forecast_with_intervals = power_out),
           "Residential_Power_Forecast_2014.xlsx")
power_out

The 2014 forecast keeps the two peak shape, with highs in January and August around 10.1 to 10.3 million KWH. The wide intervals are expected since weather drives power use and no weather data was available. The forecast is saved in Residential_Power_Forecast_2014.xlsx.

Part C (Bonus): Waterflow

1. Time sequencing and hourly aggregation

w1_raw <- read_excel("data/Waterflow_Pipe1.xlsx")
w2_raw <- read_excel("data/Waterflow_Pipe2.xlsx")

to_time <- function(x) as.POSIXct(round(x * 86400), origin = "1899-12-30", tz = "UTC")

w1 <- w1_raw |> mutate(Time = to_time(`Date Time`)) |> arrange(Time)
w2 <- w2_raw |> mutate(Time = to_time(`Date Time`)) |> arrange(Time)

bind_rows(Pipe1 = w1, Pipe2 = w2, .id = "Pipe") |>
  group_by(Pipe) |>
  summarise(readings = n(), start = min(Time), end = max(Time),
            readings_per_hour = n() / (as.numeric(difftime(max(Time), min(Time), units = "hours")) + 1))

Both files store the timestamp as an Excel day count, so I converted it to a date time first. Pipe 1 has irregular readings, about 4 per hour, while Pipe 2 already has one reading per hour. Since the two pipes cover different time spans, I decided to model them separately.

hourly <- function(df) {
  df |>
    mutate(Hour = floor_date(Time, "hour")) |>
    group_by(Hour) |>
    summarise(WaterFlow = mean(WaterFlow), .groups = "drop") |>
    as_tsibble(index = Hour)
}

p1_raw <- hourly(w1)
p2_raw <- hourly(w2)

c(Pipe1_hours = nrow(p1_raw), Pipe1_missing_hours = sum(count_gaps(p1_raw)$.n),
  Pipe2_hours = nrow(p2_raw), Pipe2_missing_hours = sum(count_gaps(p2_raw)$.n))

##         Pipe1_hours Pipe1_missing_hours         Pipe2_hours Pipe2_missing_hours 
##                 236                   4                1000                   0

fill_hours <- function(d) {
  d |>
    fill_gaps() |>
    mutate(WaterFlow = zoo::na.approx(WaterFlow))
}

p1 <- fill_hours(p1_raw)
p2 <- fill_hours(p2_raw)

c(Pipe1_hours = nrow(p1), Pipe2_hours = nrow(p2))

## Pipe1_hours Pipe2_hours 
##         240        1000

After taking the hourly mean, Pipe 2 has no gaps and Pipe 1 has 4 isolated missing hours, which I filled by linear interpolation. The ARIMA interpolation used in Parts A and B was unstable on this short series and filled the gaps with negative flows, which is physically impossible.

p1 |> autoplot(WaterFlow) + labs(title = "Pipe 1 hourly mean waterflow", x = NULL)

p2 |> autoplot(WaterFlow) + labs(title = "Pipe 2 hourly mean waterflow", x = NULL)

Both pipes fluctuate randomly around a constant level, no trend or daily cycle.

2. Is the data stationary?

stat_tests <- bind_rows(
  p1 |> features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs)) |> mutate(Pipe = "Pipe1"),
  p2 |> features(WaterFlow, list(unitroot_kpss, unitroot_ndiffs)) |> mutate(Pipe = "Pipe2")
)
stat_tests

p1 |> ACF(WaterFlow, lag_max = 48) |> autoplot() + labs(title = "Pipe 1 ACF")

p2 |> ACF(WaterFlow, lag_max = 48) |> autoplot() + labs(title = "Pipe 2 ACF")

Both series pass the KPSS test (p of 0.1, the largest value the test reports) and need no differencing. The ACFs agree, with almost every spike inside the limits and no slow decay. Both series are stationary and can be forecast, though there is very little structure to exploit.

3. Models and one week forecast

fit_pipe <- function(d) {
  d |> model(
    MEAN  = MEAN(WaterFlow),
    NAIVE = NAIVE(WaterFlow),
    ETS   = ETS(WaterFlow),
    ARIMA = ARIMA(WaterFlow)
  )
}

pipe_cv <- function(d) {
  d |> stretch_tsibble(.init = 96, .step = 24) |>
    fit_pipe() |> forecast(h = 24) |> accuracy(d) |>
    select(.model, RMSE, MAE) |> arrange(RMSE)
}

cv1 <- pipe_cv(p1); cv2 <- pipe_cv(p2)
cv1

cv2

fit_pipe(p1) |> select(ARIMA) |> report()

## Series: WaterFlow 
## Model: ARIMA(0,0,0) w/ mean 
## 
## Coefficients:
##       constant
##        19.8674
## s.e.    0.2720
## 
## sigma^2 estimated as 17.83:  log likelihood=-685.78
## AIC=1375.56   AICc=1375.61   BIC=1382.52

fit_pipe(p2) |> select(ARIMA) |> report()

## Series: WaterFlow 
## Model: ARIMA(0,0,0)(0,0,1)[24] w/ mean 
## 
## Coefficients:
##         sma1  constant
##       0.0846   39.5732
## s.e.  0.0314    0.5472
## 
## sigma^2 estimated as 256:  log likelihood=-4190.54
## AIC=8387.08   AICc=8387.1   BIC=8401.8

MEAN, ETS and ARIMA are tied for both pipes, while NAIVE is about 38% worse. The automatic ARIMA for Pipe 1 is ARIMA(0,0,0) with a mean, which is the MEAN model, and for Pipe 2 it only adds a weak seasonal term. Because the models tie within 1%, I’ll go with the simplest, the MEAN model.

pick_model <- function(cv) {
  close <- cv$.model[cv$RMSE <= min(cv$RMSE) * 1.01]
  if ("MEAN" %in% close) "MEAN" else cv$.model[1]
}

best1 <- pick_model(cv1)
best2 <- pick_model(cv2)
c(Pipe1 = best1, Pipe2 = best2)

##  Pipe1  Pipe2 
## "MEAN" "MEAN"

fit1 <- fit_pipe(p1) |> select(all_of(best1))
fit2 <- fit_pipe(p2) |> select(all_of(best2))

fit1 |> gg_tsresiduals(lag_max = 48)

fit1 |> augment() |> features(.innov, ljung_box, lag = 24, dof = 0)

fit2 |> gg_tsresiduals(lag_max = 48)

fit2 |> augment() |> features(.innov, ljung_box, lag = 24, dof = 0)

fc1 <- fit1 |> forecast(h = 168)
fc2 <- fit2 |> forecast(h = 168)

fc1 |> autoplot(p1 |> tail(24 * 7), level = c(80, 95)) + labs(title = "Pipe 1 one week forecast", x = NULL)

fc2 |> autoplot(p2 |> tail(24 * 7), level = c(80, 95)) + labs(title = "Pipe 2 one week forecast", x = NULL)

The residuals of both models look like white noise (p of ~0.58 for Pipe 1 and 0.12 for Pipe 2). The one week forecast is a flat line at the historical mean, about 19.9 for Pipe 1 and 39.6 for Pipe 2. The data are stationary and can be forecast, but the best forecast is the average level with an interval rather than a pattern.

mk <- function(fc) fc |> hilo(level = 80) |> as_tibble() |>
  transmute(DateTime = format(Hour, "%Y-%m-%d %H:%M"),
            WaterFlow = round(.mean, 3),
            Lower80 = round(`80%`$lower, 3), Upper80 = round(`80%`$upper, 3))

write_xlsx(list(Pipe1 = mk(fc1) |> select(DateTime, WaterFlow),
                Pipe2 = mk(fc2) |> select(DateTime, WaterFlow),
                Pipe1_with_intervals = mk(fc1),
                Pipe2_with_intervals = mk(fc2)),
           "Waterflow_Forecast.xlsx")

head(mk(fc1)); head(mk(fc2))

The forecast was saved in Waterflow_Forecast.xlsx, with 168 hourly values per pipe and their 80% intervals.