install.packages(c(
  "readxl",
  "dplyr",
  "ggplot2",
  "knitr",
  "writexl"
))
## Installing packages into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(knitr)
library(writexl)
install.packages("tidyr")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
library(tidyr)

Part A

atm <- read_excel("ATM624Data.xlsx")

atm <- atm |>
  mutate(
    DATE = as.Date(as.numeric(DATE), origin = "1899-12-30"),
    Cash = as.numeric(Cash)
  ) |>
  filter(
    !is.na(ATM),
    !is.na(DATE),
    !is.na(Cash)
  ) |>
  arrange(ATM, DATE)

head(atm)
summary(atm)
##       DATE                   ATM            Cash        
##  Min.   :2009-05-01   Length   :1455   Min.   :    0.0  
##  1st Qu.:2009-08-01   N.unique :   4   1st Qu.:    0.5  
##  Median :2009-10-31   N.blank  :   0   Median :   73.0  
##  Mean   :2009-10-30   Min.nchar:   4   Mean   :  155.6  
##  3rd Qu.:2010-01-29   Max.nchar:   4   3rd Qu.:  114.0  
##  Max.   :2010-04-30                    Max.   :10919.8
atm_summary <- atm |>
  group_by(ATM) |>
  summarise(
    Days = n(),
    Zero_Days = sum(Cash == 0),
    Average = mean(Cash),
    Median = median(Cash),
    Maximum = max(Cash)
  )

atm_summary
ggplot(atm, aes(x = DATE, y = Cash)) +
  geom_line() +
  facet_wrap(~ATM, scales = "free_y") +
  labs(
    title = "ATM Cash Withdrawals",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

atm_weekday <- atm |>
  mutate(
    Weekday = weekdays(DATE),
    Weekday = factor(
      Weekday,
      levels = c(
        "Monday", "Tuesday", "Wednesday",
        "Thursday", "Friday", "Saturday", "Sunday"
      )
    )
  ) |>
  group_by(ATM, Weekday) |>
  summarise(
    Average = mean(Cash),
    Median = median(Cash),
    .groups = "drop"
  )

atm_weekday
ggplot(atm_weekday,
       aes(x = Weekday, y = Median)) +
  geom_col() +
  facet_wrap(~ATM, scales = "free_y") +
  labs(
    title = "Median ATM Withdrawals by Weekday",
    x = "Day of Week",
    y = "Cash (Hundreds of Dollars)"
  ) +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

may_days <- data.frame(
  DATE = seq(
    as.Date("2010-05-01"),
    as.Date("2010-05-31"),
    by = "day"
  )
)

recent_atm <- atm |>
  filter(
    DATE >= as.Date("2010-03-06"),
    DATE <= as.Date("2010-04-30")
  ) |>
  mutate(Weekday = weekdays(DATE)) |>
  group_by(ATM, Weekday) |>
  summarise(
    Forecast = median(Cash),
    .groups = "drop"
  )

recent_atm
atm_forecast <- expand.grid(
  DATE = may_days$DATE,
  ATM = c("ATM1", "ATM2", "ATM3", "ATM4")
) |>
  mutate(Weekday = weekdays(DATE)) |>
  left_join(
    recent_atm,
    by = c("ATM", "Weekday")
  ) |>
  arrange(ATM, DATE)

atm3_average <- atm |>
  filter(ATM == "ATM3") |>
  summarise(Average = mean(Cash, na.rm = TRUE)) |>
  pull(Average)

atm_forecast <- atm_forecast |>
  mutate(
    Forecast = if_else(
      ATM == "ATM3",
      atm3_average,
      Forecast
    ),
    Forecast_USD = Forecast * 100
  )

atm_totals <- atm_forecast |>
  group_by(ATM) |>
  summarise(
    May_Forecast_USD = sum(Forecast_USD),
    .groups = "drop"
  )

atm_totals
ggplot(
  atm_forecast,
  aes(x = DATE, y = Forecast_USD)
) +
  geom_line() +
  geom_point() +
  facet_wrap(~ATM, scales = "free_y") +
  labs(
    title = "ATM Cash Withdrawal Forecasts - May 2010",
    x = "Date",
    y = "Forecasted Cash ($)"
  )

I analyzed daily cash withdrawals from four ATM machines between May 2009 and April 2010 to forecast withdrawals for May 2010. Since the Cash variable was recorded in hundreds of dollars, I multiplied the forecasted values by 100 to convert them into actual dollars.

I first examined the historical data, calculated summary statistics, and compared withdrawal patterns by day of the week. ATM1, ATM2, and ATM4 showed noticeable weekly patterns, so I used the median withdrawal for each weekday from the most recent eight weeks. Using the median helped reduce the influence of unusually large withdrawals.

ATM3 was different because 362 of its 365 historical observations were zero. For this ATM, I used the overall daily average, including zero-withdrawal days, rather than relying on the few positive observations.

The total forecasted cash withdrawals for May 2010 are approximately $1,646,003. ATM4 has the highest expected withdrawals at approximately $1.22 million, followed by ATM1 at $242,800, ATM2 at $177,650, and ATM3 at $2,234.

I chose these relatively simple forecasting methods because the ATMs showed different usage patterns. I did not apply the same ARIMA model to every ATM, and I avoided using an overall mean for ATM4 because of its extreme historical withdrawal. These forecasts are useful for cash planning, although unusual transactions and changes in ATM activity could affect the actual results.

Testing how accurately the forecasting method predicts April 2010

atm_train <- atm |>
  filter(DATE < as.Date("2010-04-01"))

atm_test <- atm |>
  filter(
    DATE >= as.Date("2010-04-01"),
    DATE <= as.Date("2010-04-30")
  )

atm_test |>
  count(ATM)
april_weekday <- atm_train |>
  filter(DATE >= as.Date("2010-02-04")) |>
  filter(ATM != "ATM3") |>
  mutate(Weekday = weekdays(DATE)) |>
  group_by(ATM, Weekday) |>
  summarise(
    Forecast = median(Cash, na.rm = TRUE),
    .groups = "drop"
  )

april_atm3 <- atm_train |>
  filter(ATM == "ATM3") |>
  summarise(
    Forecast = mean(Cash, na.rm = TRUE)
  ) |>
  mutate(
    ATM = "ATM3",
    Weekday = NA_character_
  )

april_results <- atm_test |>
  mutate(Weekday = weekdays(DATE)) |>
  left_join(april_weekday, by = c("ATM", "Weekday")) |>
  left_join(
    april_atm3 |> select(ATM, ATM3_Forecast = Forecast),
    by = "ATM"
  ) |>
  mutate(
    Forecast = if_else(
      ATM == "ATM3",
      ATM3_Forecast,
      Forecast
    )
  ) |>
  select(DATE, ATM, Cash, Forecast)

head(april_results)
april_accuracy <- april_results |>
  group_by(ATM) |>
  summarise(
    MAE = mean(abs(Cash - Forecast), na.rm = TRUE),
    RMSE = sqrt(mean((Cash - Forecast)^2, na.rm = TRUE)),
    Actual_April_USD = sum(Cash, na.rm = TRUE) * 100,
    Forecast_April_USD = sum(Forecast, na.rm = TRUE) * 100,
    .groups = "drop"
  )

april_accuracy
april_results |>
  select(DATE, ATM, Cash, Forecast) |>
  pivot_longer(
    cols = c(Cash, Forecast),
    names_to = "Type",
    values_to = "Amount"
  ) |>
  ggplot(aes(x = DATE, y = Amount, color = Type)) +
  geom_line() +
  facet_wrap(~ATM, scales = "free_y") +
  labs(
    title = "April 2010: Actual vs Forecasted ATM Withdrawals",
    x = "Date",
    y = "Cash (Hundreds of Dollars)",
    color = "Series"
  )

To evaluate my forecasting methods, I used April 2010 as a holdout period. I created forecasts using data before April and compared them with the actual withdrawals during April. This allowed me to evaluate the methods on observations that were not included in the training data.

I calculated Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) for each ATM. ATM1 had the lowest errors, with an MAE of 7.67 and RMSE of 10.54, measured in hundreds of dollars. ATM2 had an MAE of 12.95 and RMSE of 19.29. Both forecasts captured much of the recurring weekly withdrawal pattern.

ATM3 had an MAE of 8.77 and RMSE of 27.79. Its forecast was zero because the ATM had no positive withdrawals in the training period, but actual withdrawals occurred near the end of April. This suggests that historical averages may not be reliable when an ATM’s activity changes suddenly.

ATM4 had the largest errors, with an MAE of 188.68 and RMSE of 246.36. Although the forecast reflected some weekly patterns, it did not capture several large daily fluctuations.

Overall, the validation showed that the weekday median approach provided a reasonable baseline for ATM1 and ATM2. ATM3 and ATM4 were more difficult to forecast because of changing activity and unusually large withdrawals. These results highlight the limitations of the selected methods and the importance of reviewing forecasts before using them for cash planning.

Part B — Residential Power Forecasting

install.packages("readxl")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
install.packages("fpp3")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
library(readxl)
library(fpp3)
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.3 ──
## ✔ tibble      3.3.1     ✔ ggtime      1.1.0
## ✔ lubridate   1.9.5     ✔ feasts      0.5.0
## ✔ tsibble     1.2.0     ✔ fable       0.5.0
## ✔ tsibbledata 0.4.1
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
power <- read_excel("ResidentialCustomerForecastLoad-624.xlsx")

names(power)
## [1] "CaseSequence" "YYYY-MMM"     "KWH"
head(power)
summary(power)
##   CaseSequence        YYYY-MMM        KWH          
##  Min.   :733.0   Length   :192   Min.   :  770523  
##  1st Qu.:780.8   N.unique :192   1st Qu.: 5429912  
##  Median :828.5   N.blank  :  0   Median : 6283324  
##  Mean   :828.5   Min.nchar:  8   Mean   : 6502475  
##  3rd Qu.:876.2   Max.nchar:  8   3rd Qu.: 7620524  
##  Max.   :924.0                   Max.   :10655730  
##                                  NAs    :1
power_ts <- power |>
  mutate(
    Month = yearmonth(
      as.Date(paste0(`YYYY-MMM`, "-01"), format = "%Y-%b-%d")
    )
  ) |>
  as_tsibble(index = Month)

power_ts
power_ts |>
  autoplot(KWH) +
  labs(
    title = "Residential Electricity Consumption (1998-2013)",
    x = "Year",
    y = "Electricity Consumption (KWH)"
  )

power_ts |>
  gg_season(KWH) +
  labs(
    title = "Monthly Seasonal Patterns in Electricity Consumption",
    y = "KWH"
  )

power_ts |>
  gg_subseries(KWH) +
  labs(
    title = "Electricity Consumption by Month",
    y = "KWH"
  )

The residential electricity dataset contains 192 monthly observations from January 1998 through December 2013. The time series plot shows a recurring seasonal pattern, with electricity consumption generally increasing during winter and summer and decreasing during spring and fall. Consumption also appears to increase toward the end of the dataset.

One unusual observation occurs in July 2010, when electricity consumption drops significantly below the typical summer level. Since the reason for this observation is unknown, I will keep it in the initial analysis rather than assume it is an error. I compared ETS, ARIMA, and Seasonal Naïve (SNAIVE) models using 2013 as the testing period to determine which method was most appropriate for forecasting electricity consumption in 2014.

power_train <- power_ts |>
  filter(year(Month) <= 2012)

power_test <- power_ts |>
  filter(year(Month) == 2013)

nrow(power_train)
## [1] 180
nrow(power_test)
## [1] 12
power_models <- power_train |>
  model(
    ETS = ETS(KWH),
    ARIMA = ARIMA(KWH),
    SNAIVE = SNAIVE(KWH)
  )

report(power_models)
## Warning in report.mdl_df(power_models): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
power_fc_test <- power_models |>
  forecast(h = 12)

power_accuracy <- power_fc_test |>
  accuracy(power_test) |>
  select(.model, RMSE, MAE, MAPE)

power_accuracy |>
  arrange(RMSE)
power_fc_test |>
  autoplot(power_ts, level = NULL) +
  facet_wrap(~.model, ncol = 1) +
  labs(
    title = "2013 Electricity Forecast Comparison",
    x = "Year",
    y = "KWH"
  )

To determine which forecasting method was most appropriate, I divided the data into a training period from 1998 through 2012 and a testing period covering 2013. I compared three models: ETS, ARIMA, and Seasonal Naïve (SNAIVE).

The Seasonal Naïve model performed best, with an RMSE of 1,035,538 KWH, an MAE of 618,606 KWH, and a MAPE of 7.06%. ARIMA had a higher MAPE of 11.82%, while ETS had the highest at 18.04%.

These results suggest that the recurring monthly seasonal pattern is particularly important for forecasting electricity consumption. Although ETS and ARIMA are more complex methods, neither outperformed SNAIVE on the 2013 testing data. Therefore, I selected SNAIVE as the primary model for forecasting monthly residential electricity consumption in 2014.

final_power_model <- power_ts |>
  model(SNAIVE = SNAIVE(KWH))

report(final_power_model)
## Series: KWH 
## Model: SNAIVE 
## 
## sigma^2: 1381992371971.75
set.seed(12345678)

power_forecast_2014 <- final_power_model |>
  forecast(h = 12, bootstrap = TRUE, times = 1000)
power_forecast_2014 |>
  autoplot(power_ts) +
  labs(
    title = "Residential Electricity Consumption Forecast for 2014",
    x = "Year",
    y = "Electricity Consumption (KWH)"
  )

power_final_table <- power_forecast_2014 |>
  as_tibble() |>
  select(Month, Forecast_KWH = .mean)
power_final_table
power_final_table |>
  summarise(
    Total_2014_KWH = sum(Forecast_KWH),
    Average_Monthly_KWH = mean(Forecast_KWH)
  )

I used residential electricity consumption data from January 1998 through December 2013 to forecast monthly electricity usage for 2014. The data showed a clear seasonal pattern, with higher electricity consumption during winter and summer and lower consumption during spring and fall.

I compared ETS, ARIMA, and Seasonal Naïve (SNAIVE) models using data from 1998 through 2012 for training and 2013 for testing. SNAIVE performed best, with an RMSE of 1,035,538 KWH, MAE of 618,606 KWH, and MAPE of 7.06%. Based on these results, I selected SNAIVE for the final forecast.

After refitting the model using all available data, the forecasted total electricity consumption for 2014 was approximately 91,420,273 KWH, with an average monthly consumption of 7,618,356 KWH. January had the highest forecast at approximately 10,658,340 KWH.

The forecasts follow the seasonal pattern observed in previous years. The graph includes 80% and 95% prediction intervals to show the uncertainty surrounding the estimates.

One limitation is the unusually low electricity consumption recorded in July 2010. I kept this observation because its cause could not be determined from the dataset. The SNAIVE model also assumes that monthly consumption in 2014 will follow the corresponding months of 2013, so it does not account directly for changes in weather, demand, or long-term trends.

Part C

install.packages(c("fable", "fabletools", "distributional"))
## Installing packages into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
data1 <- read_excel("Waterflow_Pipe1.xlsx")
data2 <- read_excel("Waterflow_Pipe2.xlsx")

head(data1)
head(data2)
str(data1)
## tibble [1,000 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Date Time: num [1:1000] 42300 42300 42300 42300 42300 ...
##  $ WaterFlow: num [1:1000] 23.4 28 23.1 30 6 ...
str(data2)
## tibble [1,000 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Date Time: num [1:1000] 42300 42300 42300 42300 42300 ...
##  $ WaterFlow: num [1:1000] 18.8 43.1 38 36.1 31.9 ...
data1 <- data1 |>
  mutate(
    DateTime = as.POSIXct(
      `Date Time` * 86400,
      origin = "1899-12-30",
      tz = "UTC"
    )
  )

data2 <- data2 |>
  mutate(
    DateTime = as.POSIXct(
      `Date Time` * 86400,
      origin = "1899-12-30",
      tz = "UTC"
    )
  )

head(data1)
head(data2)
hourly1 <- data1 |>
  mutate(Hour = floor_date(
    round_date(DateTime, "minute"), "hour"
  )) |>
  group_by(Hour) |>
  summarise(
    WaterFlow = mean(WaterFlow, na.rm = TRUE),
    .groups = "drop"
  ) |>
  as_tsibble(index = Hour) |>
  fill_gaps()

hourly2 <- data2 |>
  mutate(Hour = floor_date(
    round_date(DateTime, "minute"), "hour"
  )) |>
  group_by(Hour) |>
  summarise(
    WaterFlow = mean(WaterFlow, na.rm = TRUE),
    .groups = "drop"
  ) |>
  as_tsibble(index = Hour) |>
  fill_gaps()
hourly1 |>
  as_tibble() |>
  summarise(
    First_Hour = min(Hour),
    Last_Hour = max(Hour),
    Hours_Recorded = n()
  )
hourly2 |>
  as_tibble() |>
  summarise(
    First_Hour = min(Hour),
    Last_Hour = max(Hour),
    Hours_Recorded = n()
  )
hourly1 |>
  gg_tsdisplay(WaterFlow, plot_type = "partial")
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

hourly2 |>
  gg_tsdisplay(WaterFlow, plot_type = "partial")

pipe1_complete <- hourly1 |>
  mutate(
    WaterFlow = approx(
      x = as.numeric(Hour),
      y = WaterFlow,
      xout = as.numeric(Hour),
      rule = 2
    )$y
  )

pipe2_complete <- hourly2 |>
  mutate(
    WaterFlow = approx(
      x = as.numeric(Hour),
      y = WaterFlow,
      xout = as.numeric(Hour),
      rule = 2
    )$y
  )

pipe1_complete |>
  gg_tsdisplay(WaterFlow, plot_type = "partial")

pipe2_complete |>
  gg_tsdisplay(WaterFlow, plot_type = "partial")

fit_pipe1 <- pipe1_complete |>
  model(
    MEAN = MEAN(WaterFlow),
    ARIMA = ARIMA(WaterFlow),
    SNAIVE = SNAIVE(WaterFlow ~ lag("day"))
  )

fit_pipe2 <- pipe2_complete |>
  model(
    MEAN = MEAN(WaterFlow),
    ARIMA = ARIMA(WaterFlow),
    SNAIVE = SNAIVE(WaterFlow ~ lag("day"))
  )

report(fit_pipe1)
## Warning in report.mdl_df(fit_pipe1): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
report(fit_pipe2)
## Warning in report.mdl_df(fit_pipe2): Model reporting is only supported for
## individual models, so a glance will be shown. To see the report for a specific
## model, use `select()` and `filter()` to identify a single model.
accuracy(fit_pipe1) |>
  select(.model, RMSE, MAE, MAPE) |>
  arrange(RMSE)
accuracy(fit_pipe2) |>
  select(.model, RMSE, MAE, MAPE) |>
  arrange(RMSE)
train_pipe1 <- pipe1_complete |>
  filter(Hour < max(Hour) - hours(47))

train_pipe2 <- pipe2_complete |>
  filter(Hour < max(Hour) - hours(47))

test_fit1 <- train_pipe1 |>
  model(
    MEAN = MEAN(WaterFlow),
    ARIMA = ARIMA(WaterFlow),
    SNAIVE = SNAIVE(WaterFlow ~ lag("day"))
  )

test_fit2 <- train_pipe2 |>
  model(
    MEAN = MEAN(WaterFlow),
    ARIMA = ARIMA(WaterFlow),
    SNAIVE = SNAIVE(WaterFlow ~ lag("day"))
  )

test_fc1 <- test_fit1 |> forecast(h = 48)
test_fc2 <- test_fit2 |> forecast(h = 48)

accuracy(test_fc1, pipe1_complete) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)
accuracy(test_fc2, pipe2_complete) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)
final_pipe1 <- pipe1_complete |>
  model(ARIMA = ARIMA(WaterFlow))

final_pipe2 <- pipe2_complete |>
  model(MEAN = MEAN(WaterFlow))

report(final_pipe1)
## Series: WaterFlow 
## Model: ARIMA(0,0,0) w/ mean 
## 
## Coefficients:
##       constant
##        19.8380
## s.e.    0.2774
## 
## sigma^2 estimated as 18.55:  log likelihood=-690.5
## AIC=1385.01   AICc=1385.06   BIC=1391.97
report(final_pipe2)
## Series: WaterFlow 
## Model: MEAN 
## 
## Mean: 39.5555 
## sigma^2: 257.5792
forecast_pipe1 <- final_pipe1 |>
  forecast(h = 168)

forecast_pipe2 <- final_pipe2 |>
  forecast(h = 168)

forecast_pipe1 |>
  as_tibble() |>
  select(Hour, .mean)
forecast_pipe2 |>
  as_tibble() |>
  select(Hour, .mean)
forecast_pipe1 |>
  autoplot(pipe1_complete) +
  labs(
    title = "Pipe 1: Seven-Day Water Flow Forecast",
    x = "Date",
    y = "Water Flow"
  )

forecast_pipe2 |>
  autoplot(pipe2_complete) +
  labs(
    title = "Pipe 2: Seven-Day Water Flow Forecast",
    x = "Date",
    y = "Water Flow"
  )

water_forecasts <- bind_rows(
  forecast_pipe1 |>
    as_tibble() |>
    transmute(
      Pipe = "Pipe 1",
      Hour,
      Forecast = .mean
    ),
  forecast_pipe2 |>
    as_tibble() |>
    transmute(
      Pipe = "Pipe 2",
      Hour,
      Forecast = .mean
    )
)

water_forecasts
water_forecasts |>
  count(Pipe)

For Part C, I analyzed two water flow datasets to forecast hourly water flow for the next seven days. I converted the original timestamps into hourly observations and averaged measurements that occurred within the same hour. Pipe 1 contained 240 hourly observations, including four missing values, while Pipe 2 contained 1,000 hourly observations with no missing values. I used interpolation to estimate the missing values in Pipe 1 before fitting the forecasting models.

I examined the time series plots, ACF, and PACF to identify possible trends and patterns. Both series fluctuated around relatively stable levels, although some autocorrelation was present. I compared MEAN, ARIMA, and Seasonal Naïve forecasting methods and evaluated their accuracy using the final 48 hours of each dataset as a testing period.

For Pipe 1, ARIMA produced the lowest testing RMSE of approximately 4.07, slightly outperforming MEAN at 4.08. For Pipe 2, MEAN performed best, with an RMSE of approximately 17.1, compared with 17.3 for ARIMA. The Seasonal Naïve method performed worse for both pipes.

Based on these results, I selected ARIMA for Pipe 1 and MEAN for Pipe 2. I refitted the selected models using the complete datasets and generated 168 hourly forecasts for each pipe, representing seven days of predicted water flow. The forecast graphs include 80% and 95% prediction intervals to illustrate uncertainty.

The two datasets ended on different dates, so their seven-day forecasts cover different calendar periods. Pipe 1’s forecast begins November 2, 2015, while Pipe 2’s begins December 3, 2015. Although ARIMA had a slightly lower testing RMSE for Pipe 1 (4.07 compared with 4.08 for MEAN), the difference was very small. The final ARIMA(0,0,0) model with a mean produces a constant forecast of approximately 19.84, so its forecasts are effectively the same as those from a MEAN model. This suggests that more complex forecasting models did not provide a meaningful improvement for these data.

Overall, the results suggest that relatively simple forecasting methods are appropriate for these datasets. The final output contains 336 hourly forecasts across both pipes. One limitation is that the datasets cover relatively short periods, so future changes in water usage may not be captured by the models.

write_xlsx(
  water_forecasts,
  "Part_C_Water_Flow_Forecasts.xlsx"
)
file.exists("Part_C_Water_Flow_Forecasts.xlsx")
## [1] TRUE