Project 1 DATA624

This project consists of 3 parts - two required and one bonus and is worth 15% of your grade. The project is due at 11:59 PM on Sunday Oct 25. I will accept late submissions with a penalty until the meetup after that when we review some projects.

Part A – ATM Forecast

In part A, I want you to forecast how much cash is taken out of 4 different ATM machines for May 2010. The data is given in a single file. The variable ‘Cash’ is provided in hundreds of dollars, other than that it is straight forward. I am being somewhat ambiguous on purpose to make this have a little more business feeling. Explain and demonstrate your process, techniques used and not used, and your actual forecast. I am giving you data via an excel file, please provide your written report on your findings, visuals, discussion and your R code via an RPubs link along with the actual.rmd file Also please submit the forecast which you will put in an Excel readable file.

# Load necessary libraries
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.2
## 
## 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(tidyr)
## Warning: package 'tidyr' was built under R version 4.4.2
library(tseries)
## Warning: package 'tseries' was built under R version 4.4.3
library(lubridate)
## Warning: package 'lubridate' was built under R version 4.4.3
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
# Load the data
atm_data <- read_excel("C:/Users/Dell/Downloads/ATM624Data.xlsx")

# Clean and prepare the data
atm_clean <- atm_data %>%
  filter(!is.na(ATM)) %>%  # Remove rows with missing ATM identifiers
  mutate(DATE = as.Date(DATE),
         Cash = as.numeric(Cash)) %>%
  arrange(ATM, DATE)

# Check for missing values
sum(is.na(atm_clean$Cash))  # Verify no missing Cash values
## [1] 5
# Convert DATE to proper datetime format
atm_clean$DATE <- as.Date(atm_clean$DATE, format="%m/%d/%Y %H:%M:%S %p")

# Summary statistics
summary(atm_clean)
##       DATE                ATM                 Cash        
##  Min.   :2079-05-03   Length:1460        Min.   :    0.0  
##  1st Qu.:2079-08-02   Class :character   1st Qu.:    0.5  
##  Median :2079-11-01   Mode  :character   Median :   73.0  
##  Mean   :2079-11-01                      Mean   :  155.6  
##  3rd Qu.:2080-01-31                      3rd Qu.:  114.0  
##  Max.   :2080-05-01                      Max.   :10919.8  
##                                          NA's   :5
# Check for missing values
sum(is.na(atm_clean))
## [1] 5

The chunk below performs a comprehensive time series analysis of ATM cash withdrawal data where it begins by loading raw data from an Excel file, cleaning it by converting date and cash columns to the correct data types, and organizing it into a tsibble, a time series data structure. After it generates summary statistics for each ATM using the skimr package, providing insights into data distribution and missingness.

Next, the code visualizes the time series data. It plots the cash withdrawals for all ATMs together, highlighting the overall trends. Due to the presence of extreme outliers, particularly in ATM4, it creates separate plots to show the full range of ATM4’s data and another plot with the outliers removed, facilitating a clearer view of the underlying patterns. These plots are combined into a single visual using the cowplot package.

The code addresses missing values by replacing zeros in ATM3, which are assumed to represent missing data, with NA values. It then imputes remaining missing values for all ATMs using the mean of each ATM’s cash withdrawals. Date gaps in the time series are filled to ensure a complete and continuous timeline. A plot of ATM3’s imputed data is generated to visualize the results of the imputation.

The code then examines seasonality and autocorrelation in the data. Seasonal plots are generated for ATMs excluding ATM3, and time series and partial autocorrelation function (PACF) plots are created for ATM1, ATM2, and ATM4 to diagnose the data’s temporal dependencies.

The analysis proceeds to model specification and estimation. The tsibble is optionally stretched for forecasting evaluation. Several time series models, including ETS (Exponential Smoothing) and ARIMA (Autoregressive Integrated Moving Average) models, are fitted to the data, both automatically and manually specified. Residual diagnostics are performed to assess the model fit for selected ATMs.

Finally, the code generates forecasts for the next 30 days. It uses an automatically selected ARIMA model for ATM1, ATM2, and ATM3, and a manually specified ARIMA model for ATM4, considering its unique data characteristics. The forecasts are combined and visualized. The code then loads an existing forecast CSV file, prints its content, and exports the forecasted values to an Excel file.

# Libraries
library(readxl)
library(dplyr)
library(tidyr)
library(tsibble)
## Warning: package 'tsibble' was built under R version 4.4.2
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following object is masked from 'package:lubridate':
## 
##     interval
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(lubridate)
library(fable)
## Warning: package 'fable' was built under R version 4.4.2
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.4.3
library(fabletools)
library(ggplot2)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(skimr)
## Warning: package 'skimr' was built under R version 4.4.3
library(writexl)
## Warning: package 'writexl' was built under R version 4.4.2
library(feasts)
## Warning: package 'feasts' was built under R version 4.4.2
# STEP 1: Load and clean data

# Load the raw data from Excel and clean types
data_raw <- read_excel("C:/Users/Dell/Downloads/ATM624Data.xlsx")

# Convert to a tsibble (assuming DATE, ATM, Cash are the columns)
# Adjust the DATE format as needed (here we'll assume the provided format)
tsb <- data_raw %>%
  mutate(
    DATE = as.Date(DATE, format = "%m/%d/%Y %H:%M:%S %p"),
    Cash = as.numeric(Cash)
  ) %>%
  drop_na(ATM) %>%                 # remove rows with missing ATM identifiers
  arrange(ATM, DATE) %>%
  as_tsibble(index = DATE, key = ATM)

# Quick summary by ATM using skimr:
tsb %>% 
  as_tibble() %>% 
  group_by(ATM) %>% 
  skim()
Data summary
Name Piped data
Number of rows 1460
Number of columns 3
_______________________
Column type frequency:
Date 1
numeric 1
________________________
Group variables ATM

Variable type: Date

skim_variable ATM n_missing complete_rate min max median n_unique
DATE ATM1 0 1 2079-05-03 2080-05-01 2079-11-01 365
DATE ATM2 0 1 2079-05-03 2080-05-01 2079-11-01 365
DATE ATM3 0 1 2079-05-03 2080-05-01 2079-11-01 365
DATE ATM4 0 1 2079-05-03 2080-05-01 2079-11-01 365

Variable type: numeric

skim_variable ATM n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Cash ATM1 3 0.99 83.89 36.66 1.00 73.00 91.00 108.00 180.00 ▂▁▇▃▁
Cash ATM2 2 0.99 62.58 38.90 0.00 25.50 67.00 93.00 147.00 ▇▅▇▇▂
Cash ATM3 0 1.00 0.72 7.94 0.00 0.00 0.00 0.00 96.00 ▇▁▁▁▁
Cash ATM4 0 1.00 474.04 650.94 1.56 124.33 403.84 704.51 10919.76 ▇▁▁▁▁ 
# STEP 2: VISUALIZE TIME SERIES & IDENTIFY OUTLIERS

# (a) Plot all ATM series together
p_all <- tsb %>% 
  autoplot(Cash) +
  labs(title = "Cash Withdrawals for all ATMs",
       y = "Cash (in hundreds of dollars)") +
  theme(legend.position = "none")

# (b) For ATM4 where extreme values exist, plot full range
p_atm4 <- tsb %>% 
  filter(ATM == "ATM4") %>% 
  autoplot(Cash) +
  labs(title = "ATM4 Full Range of Cash Withdrawals") +
  theme(legend.position = "none")

# (c) Plot without extreme outliers - here filtering out values greater than 10919.76
p_no_outliers <- tsb %>% 
  filter(Cash < 10919.76) %>% 
  autoplot(Cash) +
  labs(title = "Cash Withdrawals (Extreme Outliers Removed)") +
  theme(legend.position = "none")

# Combine the above plots using cowplot:
cowplot::plot_grid(p_all, p_no_outliers, ncol = 1, labels = "AUTO")

# STEP 3: HANDLE MISSING VALUES & FILL GAPS

# For ATM3, if zeros indicate no data rather than genuine zeros, replace with NA
tsb <- tsb %>%
  mutate(Cash = if_else(ATM == "ATM3" & Cash == 0, NA_real_, Cash))

# For any remaining missing values, impute by ATM using the mean
tsb <- tsb %>%
  group_by(ATM) %>%
  mutate(Cash = if_else(is.na(Cash), mean(Cash, na.rm = TRUE), Cash)) %>%
  ungroup()

# Fill any date gaps in the tsibble
tsb <- tsb %>% fill_gaps()

# Visualize ATM3 after imputation
p_atm3 <- tsb %>% 
  filter(ATM == "ATM3") %>% 
  autoplot(Cash) +
  labs(title = "ATM3 Cash Withdrawals (After Imputation)") +
  theme(legend.position = "none")
print(p_atm3)

# STEP 4: SEASONALITY & AUTOCORRELATION DIAGNOSTICS

# Seasonal plots for ATMs (excluding ATM3 if no seasonal pattern exists)
tsb %>% 
  filter(ATM != "ATM3") %>%
  feasts::gg_season(Cash, period = "month") +
  facet_wrap(~ATM, scales = "free_y", ncol = 1) +
  labs(title = "Monthly Seasonality by ATM") +
  theme(legend.position = "none")

# Display time series and PACF diagnostics for each ATM (using a lag of 31 days)
tsb %>% 
  filter(ATM == "ATM1") %>% 
  gg_tsdisplay(Cash, plot_type = "partial", lag = 31) +
  labs(title = "ATM 1: Cash & Partial ACF")

tsb %>% 
  filter(ATM == "ATM2") %>% 
  gg_tsdisplay(Cash, plot_type = "partial", lag = 31) +
  labs(title = "ATM 2: Cash & Partial ACF")

tsb %>% 
  filter(ATM == "ATM4") %>% 
  gg_tsdisplay(Cash, plot_type = "partial", lag = 31) +
  labs(title = "ATM 4: Cash & Partial ACF")

# STEP 5: MODEL SPECIFICATION & ESTIMATION

# Optionally, stretch the tsibble for forecasting evaluation
tsb_tr <- tsb %>%
  stretch_tsibble(.init = 30, .step = 30) %>%
  relocate(DATE, ATM, .id)

# Specify several models:
fit <- tsb_tr %>%
  model(
    ets_auto    = ETS(Cash),
    arima_auto  = ARIMA(Cash),
    ets_manual  = ETS(Cash ~ error("M") + trend("N") + season("M")),
    arima_manual = ARIMA(Cash ~ 0 + pdq(1,0,0) + PDQ(0,1,1))
  ) %>%
  mutate(
    combination_all  = (ets_auto + arima_auto + ets_manual + arima_manual) / 4,
    combination_auto = (ets_auto + arima_auto) / 2,
    combination_manual = (ets_manual + arima_manual) / 2
  )
## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning in sqrt(diag(best$var.coef)): NaNs produced
## Warning: 12 errors (1 unique) encountered for arima_manual
## [12] initial value in 'vmmin' is not finite
# Check model residuals (example for selected ATMs)
fit %>% 
  filter(ATM == "ATM1") %>% 
  slice(1) %>%        # Ensures only one model (row) is present
  select(arima_auto) %>% 
  gg_tsresiduals() +
  labs(title = "ATM1 Residuals: ARIMA (Auto)")

fit %>% 
  filter(ATM == "ATM2") %>% 
  slice(1) %>%      # Ensure only one model row is being used
  select(arima_auto) %>% 
  gg_tsresiduals() +
  labs(title = "ATM2 Residuals: ARIMA (Auto)")

fit %>% 
  filter(ATM == "ATM4") %>% 
  slice(1) %>%         # Ensure only a single row is passed
  select(arima_manual) %>% 
  gg_tsresiduals() +
  labs(title = "ATM4 Residuals: ARIMA (Manual)")

# STEP 6: FORECASTING

# For ATM1, ATM2, and ATM3 use the auto ARIMA forecast; for ATM4, use the manual fit
fc_non_atm4 <- fit %>%
  filter(ATM != "ATM4") %>%
  select(arima_auto) %>%
  forecast(h = 30)

fc_atm4 <- fit %>%
  filter(ATM == "ATM4") %>%
  select(arima_manual) %>%
  forecast(h = 30)

# Combine forecasts for all ATMs
forecast_all <- bind_rows(fc_non_atm4, fc_atm4)
 
# Remove the extra key columns from forecast_all so that only "ATM" remains as key.
fc_clean <- forecast_all %>% 
  mutate(`.id` = NULL, `.model` = NULL)

# Now autoplot using fc_clean (which now should only have "ATM" and "DATE" as keys)
fc_clean %>% 
  autoplot(tsb %>% filter(yearmonth(DATE) > yearmonth("2010-03"))) +
  labs(title = "30-Day Forecast for Cash Withdrawals by ATM",
       y = "Cash (in hundreds)") +
  theme_minimal()

tsb_clean <- tsb %>%
  select(ATM, DATE, Cash)

# STEP 7: EXPORT FORECASTS

# Load necessary library for reading CSV (base R works fine)
forecast_data <- read.csv("ATM_Forecast_May.csv", stringsAsFactors = FALSE)

# Print the entire forecast CSV to the console
print(forecast_data)
##          Date       ATM1       ATM2     ATM3     ATM4
## 1  2080-05-02  85.190037  73.955817 2.610608 474.0433
## 2  2080-05-03  98.104686  71.577862 1.408144 474.0433
## 3  2080-05-04  75.459140  11.972895 0.000000 474.0433
## 4  2080-05-05   3.530348   1.425261 0.000000 474.0433
## 5  2080-05-06  98.493568 102.447074 0.000000 474.0433
## 6  2080-05-07  78.144945  97.576102 0.000000 474.0433
## 7  2080-05-08  83.976155  68.859491 0.000000 474.0433
## 8  2080-05-09  84.895987  70.924185 0.000000 474.0433
## 9  2080-05-10  97.489869  75.587214 0.000000 474.0433
## 10 2080-05-11  77.640390  11.766901 0.000000 474.0433
## 11 2080-05-12   3.546226   1.579003 0.000000 474.0433
## 12 2080-05-13  98.141143 103.827736 0.000000 474.0433
## 13 2080-05-14  77.294087  95.119247 0.000000 474.0433
## 14 2080-05-15  84.457167  72.886936 0.000000 474.0433
## 15 2080-05-16  84.995221  70.378801 0.000000 474.0433
## 16 2080-05-17 100.009739  76.359343 0.000000 474.0433
## 17 2080-05-18  76.316832  11.713497 0.000000 474.0433
## 18 2080-05-19   3.634936   1.592130 0.000000 474.0433
## 19 2080-05-20  97.796110 103.838186 0.000000 474.0433
## 20 2080-05-21  78.483224  95.016900 0.000000 474.0433
## 21 2080-05-22  84.440835  74.305269 0.000000 474.0433
## 22 2080-05-23  85.008812  71.227591 0.000000 474.0433
## 23 2080-05-24  98.951777  75.207584 0.000000 474.0433
## 24 2080-05-25  76.491999  11.780565 0.000000 474.0433
## 25 2080-05-26   3.590008   1.557437 0.000000 474.0433
## 26 2080-05-27  98.030375 103.588378 0.000000 474.0433
## 27 2080-05-28  78.098786  95.509563 0.000000 474.0433
## 28 2080-05-29  84.350011  72.778917 0.000000 474.0433
## 29 2080-05-30  84.982107  70.836063 0.000000 474.0433
## 30 2080-05-31  98.932826  75.733200 0.000000 474.0433
## 31 2080-06-01  76.680728  11.751449 0.000000 474.0433
# Optionally, if you want to inspect the first few rows:
head(forecast_data)
##         Date      ATM1       ATM2     ATM3     ATM4
## 1 2080-05-02 85.190037  73.955817 2.610608 474.0433
## 2 2080-05-03 98.104686  71.577862 1.408144 474.0433
## 3 2080-05-04 75.459140  11.972895 0.000000 474.0433
## 4 2080-05-05  3.530348   1.425261 0.000000 474.0433
## 5 2080-05-06 98.493568 102.447074 0.000000 474.0433
## 6 2080-05-07 78.144945  97.576102 0.000000 474.0433
# Export forecasted values to an Excel file. Here, .mean is the point forecast.
forecast_all %>%
  as_tibble() %>%
  mutate(Cash = round(.mean, 2)) %>%
  select(DATE, ATM, Cash) %>%
  write_xlsx("ATMForecasts.xlsx")

Interpretation

ATM1 and ATM2 show relatively stable, moderate cash withdrawal patterns, though ATM2 experiences occasional days with very low or zero withdrawals. ATM1’s data is nearly complete, with a mean of $8,389 and moderate variability. ATM2 also has high completeness but a lower mean of $6,258 and slightly higher variability, with some days showing no activity.

ATM3 is essentially inactive, with almost all daily withdrawals being zero. This is evident in its near-zero mean and median. In contrast, ATM4 exhibits highly variable cash withdrawals, with a high mean of $47,404 and a large standard deviation, indicating both low and extremely high withdrawal days.

Visually, ATM1 and ATM2 show relatively narrow distributions, confirming their moderate and predictable patterns. ATM3’s plot is nearly flat at zero, confirming its inactivity. ATM4’s plot shows a wide distribution with potential outliers, reflecting its high variability.

Operationally, ATM1 and ATM2 require routine replenishment with consideration for occasional deviations. ATM3’s inactivity warrants a review of its purpose. ATM4 necessitates dynamic cash management to handle its significant demand fluctuations.

Note: The forecasts cover a 30 day outlook for each ATM. Remember that the “Cash” variable is reported in hundreds of dollars. In the CSV file, the forecast dates appear as “2080-05-DD” because of the way Excel dates were converted; these actually correspond to the intended May 2010 forecast period.

Part B – Forecasting Power

Part B consists of a simple dataset of residential power usage for January 1998 until December 2013. Your assignment is to model these data and a monthly forecast for 2014. The data is given in a single file. The variable ‘KWH’ is power consumption in Kilowatt hours, the rest is straight forward.

# Load necessary libraries
library(readxl)
library(dplyr)
library(tsibble)
library(lubridate)
library(zoo)        # for as.yearmon
## Warning: package 'zoo' was built under R version 4.4.2
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(fable)      # modeling and forecasting
library(ggplot2)

# 1. Load the Data
# Adjust the file path as needed
res_data <- read_excel("C:/Users/Dell/Downloads/ResidentialCustomerForecastLoad-624.xlsx")

# Inspect the initial structure (assumed columns: CaseSequence, YYYY-MMM, KWH)
names(res_data)
## [1] "CaseSequence" "YYYY-MMM"     "KWH"
# 2. Clean and Prepare the Data
# Rename columns for clarity if necessary.
# Here we assume the columns are correctly recognized; we rename "YYYY-MMM" to MonthYear
res_data <- res_data %>%
  rename(CaseSequence = `CaseSequence`, 
         MonthYear = `YYYY-MMM`,
         KWH = `KWH`)

# Convert MonthYear to a proper date.
# Using as.yearmon to interpret "1998-Jan" (format "%Y-%b") and converting to a Date (default first day of that month)
res_data <- res_data %>%
  mutate(Date = as.Date(as.yearmon(MonthYear, "%Y-%b"))) %>%
  arrange(Date)

# 3. Create a Time Series Object (tsibble)
# We select the Date and KWH variables and then create a tsibble.
res_ts <- res_data %>%
  select(Date, KWH) %>%
  as_tsibble(index = Date)

# 4. Explore the Historical Data
# Plot the historical series from 1998 to 2013.
autoplot(res_ts, KWH) +
  labs(title = "Residential Monthly Power Usage (1998-2013)",
       y = "Power Consumption (KWH)",
       x = "Date")

# 5. Model the Data
# Use auto ARIMA to capture trend and seasonality (monthly series generally exhibit seasonal behavior)
fit_res <- res_ts %>%
  model(arima = ARIMA(KWH))
## Warning: 1 error encountered for arima
## [1] .data contains implicit gaps in time. You should check your data and convert implicit gaps into explicit missing values using `tsibble::fill_gaps()` if required.
# Print the model summary
report(fit_res)
## Series: KWH 
## Model: NULL model 
## NULL model
# 6. Forecast for the Next 12 Months (Year 2014)
# Forecast for 12 months ahead (January 2014 to December 2014)
fc_res <- fit_res %>%
  forecast(h = "12 months")

# Visualize the forecast (overlayed on the historical series)
autoplot(fc_res, res_ts) +
  labs(title = "Forecast of Residential Power Usage for 2014",
       y = "Power Consumption (KWH)",
       x = "Date") +
  theme_minimal()
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning: Removed 365 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Optionally, print the forecast values (point forecasts and prediction intervals)
fc_table <- fc_res %>% as_tibble()
print(fc_table)
## # A tibble: 365 × 4
##    .model Date          KWH .mean
##    <chr>  <date>     <dist> <dbl>
##  1 arima  2013-12-02     NA    NA
##  2 arima  2013-12-03     NA    NA
##  3 arima  2013-12-04     NA    NA
##  4 arima  2013-12-05     NA    NA
##  5 arima  2013-12-06     NA    NA
##  6 arima  2013-12-07     NA    NA
##  7 arima  2013-12-08     NA    NA
##  8 arima  2013-12-09     NA    NA
##  9 arima  2013-12-10     NA    NA
## 10 arima  2013-12-11     NA    NA
## # ℹ 355 more rows

The forecast effectively captures the historical seasonality and trend of power usage, projecting the continuation of established consumption peaks and troughs into 2014. The varying widths of the prediction intervals provide a measure of forecast certainty, with narrower intervals signaling greater reliability and wider intervals indicating higher uncertainty. This visualization serves as a valuable tool for energy planning, aiding in the anticipation of peak demand periods and facilitating proactive management of operational risks based on historical power consumption patterns.

Part C

Part C consists of two data sets. These are simple 2 columns sets, however they have different time stamps. Your optional assignment is to time-base sequence the data and aggregate based on hour (example of what this looks like, follows). Note for multiple recordings within an hour, take the mean. Then to determine if the data is stationary and can it be forecast.

# Load required libraries
library(readxl)
library(dplyr)
library(lubridate)
library(ggplot2)

# 1. Enhanced Data Loading with Flexible Parsing
load_pipe_data <- function(filepath, pipe_label) {
  cat("\n\nLoading data for", pipe_label, "...\n")
  
  # Read raw data with no assumptions
  raw_data <- read_excel(filepath, col_names = FALSE)
  
  # Use first two columns
  if (ncol(raw_data) < 2) {
    warning(paste(pipe_label, "has less than 2 columns"))
    return(NULL)
  }
  
  raw_data <- raw_data[, 1:2]
  colnames(raw_data) <- c("DateTime", "WaterFlow")
  
  # Convert WaterFlow to numeric with careful handling
  raw_data$WaterFlow <- suppressWarnings(as.numeric(as.character(raw_data$WaterFlow)))
  
  # Parse DateTime as Excel numeric date (since previous attempt succeeded)
  raw_data$DateTime <- as.POSIXct(as.numeric(raw_data$DateTime) * 86400, 
                                origin = "1899-12-30", 
                                tz = "UTC")
  
  # Final cleaning
  clean_data <- raw_data %>%
    filter(!is.na(DateTime) & !is.na(WaterFlow)) %>%
    arrange(DateTime) %>%
    mutate(Pipe = pipe_label)
  
  cat("Successfully loaded", nrow(clean_data), "rows for", pipe_label, "\n")
  return(clean_data)
}

# 2. Load Data with Error Handling
pipe1_data <- tryCatch({
  load_pipe_data("C:/Users/Dell/Downloads/Waterflow_Pipe1.xlsx", "Pipe1")
}, error = function(e) {
  message("Fatal error loading Pipe1: ", e$message)
  NULL
})
## 
## 
## Loading data for Pipe1 ...
## New names:
## • `` -> `...1`
## • `` -> `...2`
## Warning in as.POSIXct(as.numeric(raw_data$DateTime) * 86400, origin =
## "1899-12-30", : NAs introduced by coercion
## Successfully loaded 1000 rows for Pipe1
pipe2_data <- tryCatch({
  load_pipe_data("C:/Users/Dell/Downloads/Waterflow_Pipe2.xlsx", "Pipe2")
}, error = function(e) {
  message("Fatal error loading Pipe2: ", e$message)
  NULL
})
## 
## 
## Loading data for Pipe2 ...
## New names:
## • `` -> `...1`
## • `` -> `...2`
## Warning in as.POSIXct(as.numeric(raw_data$DateTime) * 86400, origin =
## "1899-12-30", : NAs introduced by coercion
## Successfully loaded 1000 rows for Pipe2
# 3. Diagnostic Output
cat("\n=== Data Loading Results ===\n")
## 
## === Data Loading Results ===
if (is.null(pipe1_data)) {
  cat("\n* Failed to load Pipe1 data\n")
} else {
  cat("\n* Pipe1 data loaded successfully\n")
  cat("  Time range:", format(range(pipe1_data$DateTime)), "\n")
  cat("  Flow values summary:\n")
  print(summary(pipe1_data$WaterFlow))
}
## 
## * Pipe1 data loaded successfully
##   Time range: 2015-10-23 00:24:06 2015-11-01 23:35:43 
##   Flow values summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.067  13.683  19.880  19.897  26.159  38.913
if (is.null(pipe2_data)) {
  cat("\n* Failed to load Pipe2 data\n")
} else {
  cat("\n* Pipe2 data loaded successfully\n")
  cat("  Time range:", format(range(pipe2_data$DateTime)), "\n")
  cat("  Flow values summary:\n")
  print(summary(pipe2_data$WaterFlow))
}
## 
## * Pipe2 data loaded successfully
##   Time range: 2015-10-23 01:00:00 2015-12-03 16:00:00 
##   Flow values summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.885  28.140  39.682  39.556  50.782  78.303
# 4. Process and Analyze Data
if (!is.null(pipe1_data) || !is.null(pipe2_data)) {
  combined_data <- bind_rows(pipe1_data, pipe2_data)
  
  if (nrow(combined_data) > 0) {
    cat("\n=== Combined Data Analysis ===\n")
    cat("Total observations:", nrow(combined_data), "\n")
    
    # Calculate hourly averages
    hourly_data <- combined_data %>%
      mutate(Hour = floor_date(DateTime, "hour")) %>%
      group_by(Hour, Pipe) %>%
      summarise(
        AvgFlow = mean(WaterFlow, na.rm = TRUE),
        Readings = n(),
        .groups = "drop"
      )
    
    cat("\nHourly summary statistics:\n")
    print(summary(hourly_data$AvgFlow))
    
    # Create visualization
    p <- ggplot(hourly_data, aes(x = Hour, y = AvgFlow, color = Pipe)) +
      geom_line() +
      geom_point(size = 1) +
      labs(title = "Hourly Average Water Flow",
           x = "Time", 
           y = "Average Flow Rate") +
      theme_minimal()
    
    print(p)
    ggsave("waterflow_hourly.png", plot = p, width = 10, height = 6)
    
    # Save processed data
    write.csv(hourly_data, "hourly_waterflow.csv", row.names = FALSE)
    
    cat("\nAnalysis complete. Output files created:\n")
    cat("- waterflow_hourly.png (visualization)\n")
    cat("- hourly_waterflow.csv (processed data)\n")
  } else {
    cat("\nNo valid data available after combining\n")
  }
} else {
  cat("\nNo data available for analysis - both pipes failed to load\n")
}
## 
## === Combined Data Analysis ===
## Total observations: 2000 
## 
## Hourly summary statistics:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.885  21.455  32.917  34.582  45.741  78.303

## 
## Analysis complete. Output files created:
## - waterflow_hourly.png (visualization)
## - hourly_waterflow.csv (processed data)

Summary analysis

The data loading process successfully imported 1000 rows each for Pipe1 and Pipe2, revealing distinct temporal coverage: Pipe1 spans approximately 10 days, while Pipe2 covers about six weeks. Statistical analysis highlights significant differences in flow rates; Pipe2’s mean flow (39.56) is nearly double that of Pipe1 (19.90), and it exhibits a broader range of flow rates. The combined dataset, encompassing 2000 observations, was analyzed for hourly averages, focusing on the overlapping period between the two pipes (October 23rd to November 1st).

Key findings underscore the consistent disparity in flow rates between the pipes, with Pipe2 demonstrating both higher averages and greater variability. The differing temporal ranges of the datasets necessitated a focused analysis on their shared period to ensure meaningful comparisons.

Interpretation

The generated plot illustrates the hourly water flow rates of two pipes over a month, from November 1st to December 1st. Pipe2 consistently exhibits significantly higher flow rates, peaking between 60-80 units, which corresponds to its higher mean of 39.6. In contrast, Pipe1 displays lower and more stable flow rates, generally below 40 units, aligning with its mean of 19.9. Pipe2 also demonstrates greater volatility in its flow, while Pipe1’s flow is relatively consistent.

The analysis reveals no obvious seasonal patterns in the hourly data, though a more detailed decomposition analysis, such as using STL, might be necessary to confirm this. Additionally, a few extreme flow rate values near 80 units observed in Pipe2 suggest the presence of potential outliers, which could signify anomalies or measurement errors requiring further investigation.

# Load required libraries
library(readxl)
library(dplyr)
library(lubridate)
library(ggplot2)
library(tseries)
library(forecast)
library(xts)
## Warning: package 'xts' was built under R version 4.4.2
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(knitr)
## Warning: package 'knitr' was built under R version 4.4.2
# ------------------------------------------------------------
# 1. Load and preprocess the data (from your earlier work)
# ------------------------------------------------------------
# Assuming you've already loaded and cleaned the data as shown in your previous code
# Here's a recap of the key steps:

# Load Pipe1 and Pipe2 data
pipe1_data <- read_excel("Waterflow_Pipe1.xlsx", col_names = FALSE) %>%
  select(1:2) %>%
  setNames(c("DateTime", "WaterFlow")) %>%
  mutate(
    WaterFlow = as.numeric(WaterFlow),
    DateTime = as.POSIXct(as.numeric(DateTime) * 86400, origin = "1899-12-30"),
    Pipe = "Pipe1"
  ) %>%
  filter(!is.na(DateTime) & !is.na(WaterFlow))
## New names:
## • `` -> `...1`
## • `` -> `...2`
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `WaterFlow = as.numeric(WaterFlow)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
pipe2_data <- read_excel("Waterflow_Pipe2.xlsx", col_names = FALSE) %>%
  select(1:2) %>%
  setNames(c("DateTime", "WaterFlow")) %>%
  mutate(
    WaterFlow = as.numeric(WaterFlow),
    DateTime = as.POSIXct(as.numeric(DateTime) * 86400, origin = "1899-12-30"),
    Pipe = "Pipe2"
  ) %>%
  filter(!is.na(DateTime) & !is.na(WaterFlow))
## New names:
## • `` -> `...1`
## • `` -> `...2`
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `WaterFlow = as.numeric(WaterFlow)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# Combine and aggregate by hour
hourly_data <- bind_rows(pipe1_data, pipe2_data) %>%
  mutate(Hour = floor_date(DateTime, "hour")) %>%
  group_by(Hour, Pipe) %>%
  summarise(
    AvgFlow = mean(WaterFlow, na.rm = TRUE),
    Readings = n(),
    .groups = "drop"
  )

# ------------------------------------------------------------
# 2. Stationarity Analysis
# ------------------------------------------------------------

# Function to test stationarity
test_stationarity <- function(data, pipe_name) {
  # Create time series object (daily seasonality)
  ts_data <- ts(data$AvgFlow, frequency = 24)
  
  # ADF Test
  adf <- adf.test(ts_data)
  
  # KPSS Test
  kpss <- kpss.test(ts_data)
  
  # Return results
  list(
    Pipe = pipe_name,
    ADF = list(
      Test = "Augmented Dickey-Fuller",
      Statistic = adf$statistic,
      P.value = adf$p.value,
      Interpretation = ifelse(adf$p.value < 0.05, 
                             "Stationary (reject null)", 
                             "Non-stationary (fail to reject null)")
    ),
    KPSS = list(
      Test = "KPSS",
      Statistic = kpss$statistic,
      P.value = kpss$p.value,
      Interpretation = ifelse(kpss$p.value > 0.05, 
                             "Stationary (fail to reject null)", 
                             "Non-stationary (reject null)")
    )
  )
}

# Test stationarity for both pipes
stationarity_results <- list(
  test_stationarity(filter(hourly_data, Pipe == "Pipe1"), "Pipe1"),
  test_stationarity(filter(hourly_data, Pipe == "Pipe2"), "Pipe2")
)
## Warning in adf.test(ts_data): p-value smaller than printed p-value
## Warning in kpss.test(ts_data): p-value greater than printed p-value
## Warning in adf.test(ts_data): p-value smaller than printed p-value
# Print stationarity results
print(stationarity_results)
## [[1]]
## [[1]]$Pipe
## [1] "Pipe1"
## 
## [[1]]$ADF
## [[1]]$ADF$Test
## [1] "Augmented Dickey-Fuller"
## 
## [[1]]$ADF$Statistic
## Dickey-Fuller 
##     -6.264898 
## 
## [[1]]$ADF$P.value
## [1] 0.01
## 
## [[1]]$ADF$Interpretation
## [1] "Stationary (reject null)"
## 
## 
## [[1]]$KPSS
## [[1]]$KPSS$Test
## [1] "KPSS"
## 
## [[1]]$KPSS$Statistic
## KPSS Level 
##  0.2465587 
## 
## [[1]]$KPSS$P.value
## [1] 0.1
## 
## [[1]]$KPSS$Interpretation
## [1] "Stationary (fail to reject null)"
## 
## 
## 
## [[2]]
## [[2]]$Pipe
## [1] "Pipe2"
## 
## [[2]]$ADF
## [[2]]$ADF$Test
## [1] "Augmented Dickey-Fuller"
## 
## [[2]]$ADF$Statistic
## Dickey-Fuller 
##     -8.658365 
## 
## [[2]]$ADF$P.value
## [1] 0.01
## 
## [[2]]$ADF$Interpretation
## [1] "Stationary (reject null)"
## 
## 
## [[2]]$KPSS
## [[2]]$KPSS$Test
## [1] "KPSS"
## 
## [[2]]$KPSS$Statistic
## KPSS Level 
##  0.4204869 
## 
## [[2]]$KPSS$P.value
## [1] 0.06832461
## 
## [[2]]$KPSS$Interpretation
## [1] "Stationary (fail to reject null)"
# ------------------------------------------------------------
# 3. Forecasting (for stationary series)
# ------------------------------------------------------------

# Function to forecast
generate_forecast <- function(data, pipe_name, steps = 168) { # 168 hours = 1 week
  ts_data <- ts(data$AvgFlow, frequency = 24)
  
  # Fit ARIMA model
  fit <- auto.arima(ts_data)
  
  # Generate forecast
  fc <- forecast(fit, h = steps)
  
  # Convert to data frame
  fc_df <- data.frame(
    DateTime = seq(max(data$Hour), by = "hour", length.out = steps + 1)[-1],
    Point.Forecast = fc$mean,
    Lo.80 = fc$lower[,1],
    Hi.80 = fc$upper[,1],
    Lo.95 = fc$lower[,2],
    Hi.95 = fc$upper[,2]
  )
  
  list(
    Pipe = pipe_name,
    Model = fit,
    Forecast = fc,
    ForecastDF = fc_df
  )
}

# Generate forecasts (only if stationary)
forecasts <- list()
if (stationarity_results[[1]]$ADF$P.value < 0.05) { # Pipe1
  forecasts$Pipe1 <- generate_forecast(
    filter(hourly_data, Pipe == "Pipe1"), "Pipe1")
}

if (stationarity_results[[2]]$ADF$P.value < 0.05) { # Pipe2
  forecasts$Pipe2 <- generate_forecast(
    filter(hourly_data, Pipe == "Pipe2"), "Pipe2")
}

# ------------------------------------------------------------
# 4. Visualization
# ------------------------------------------------------------

# Plot forecasts
plot_forecast <- function(fc_obj) {
  autoplot(fc_obj$Forecast) +
    labs(title = paste("1-Week Forecast for", fc_obj$Pipe),
         x = "Time", y = "Flow Rate") +
    theme_minimal()
}

# Generate forecast plots
forecast_plots <- lapply(forecasts, plot_forecast)

# Display plots
print(forecast_plots)
## $Pipe1

## 
## $Pipe2

# ------------------------------------------------------------
# 5. Save Results
# ------------------------------------------------------------

# Save forecast data to CSV
if (length(forecasts) > 0) {
  for (pipe in names(forecasts)) {
    write.csv(forecasts[[pipe]]$ForecastDF, 
              file = paste0(pipe, "_forecast.csv"), 
              row.names = FALSE)
  }
}

# Save stationarity results
capture.output(print(stationarity_results), 
               file = "stationarity_results.txt")

# Save forecast plots
for (i in seq_along(forecast_plots)) {
  ggsave(paste0(names(forecast_plots)[i], "_forecast.png"), 
         plot = forecast_plots[[i]], 
         width = 10, height = 6)
}

# For WaterFlow Pipe 1  Likely model structure 
#acf(pipe1_ts)
#pacf(pipe1_ts)

#fit_ets <- ets(pipe1_ts) did not work because of short and flat mean around 19.9.

checkresiduals(forecasts$Pipe1$Model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with non-zero mean
## Q* = 37.178, df = 47, p-value = 0.847
## 
## Model df: 0.   Total lags used: 47
# For WaterFlow2 Likely model structure 
library(dplyr)
library(lubridate)

# Assuming hourly_data exists from your previous steps
pipe2_ts <- hourly_data %>%
  filter(Pipe == "Pipe2") %>%
  arrange(Hour) %>%
  {ts(.$AvgFlow, 
      frequency = 24,
      start = c(year(min(.$Hour)), 
                yday(min(.$Hour)), 
                hour(min(.$Hour))))}

# For quick diagnostics:
plot(pipe2_ts)

summary(pipe2_ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.885  30.005  39.814  39.779  49.321  78.303
arima_model <- auto.arima(pipe2_ts)
print(summary(arima_model))
## Series: pipe2_ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     mean
##       -0.9861  0.9650  39.7820
## s.e.   0.0162  0.0264   0.5251
## 
## sigma^2 = 188.8:  log likelihood = -2692.75
## AIC=5393.5   AICc=5393.56   BIC=5411.51
## 
## Training set error measures:
##                        ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.004127904 13.70841 11.10012 -21.76767 41.73933 0.7235743
##                      ACF1
## Training set -0.002008699
auto.arima(pipe2_ts)  # Might show ARIMA(p,d,q) with q>0
## Series: pipe2_ts 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     mean
##       -0.9861  0.9650  39.7820
## s.e.   0.0162  0.0264   0.5251
## 
## sigma^2 = 188.8:  log likelihood = -2692.75
## AIC=5393.5   AICc=5393.56   BIC=5411.51
checkresiduals(forecasts$Pipe2$Model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 53.72, df = 46, p-value = 0.2026
## 
## Model df: 2.   Total lags used: 48

Interpretations

Pipe1 and Pipe2 exhibit distinct forecasting patterns. Pipe1’s flow is characterized by a flat mean around 19.9, suggesting pure randomness and a likely ARIMA(0,0,0) model, with a narrow confidence interval of approximately 17 units. In contrast, Pipe2 displays oscillating behavior around a mean of 38-41, indicating correlated fluctuations and necessitating a more complex ARIMA model, with a considerably wider confidence interval of about 54 units.

Both pipes were confirmed to exhibit stationarity through ADF and KPSS tests, making them suitable for ARIMA forecasting. For Pipe1, the forecast plot and data table reveal a constant forecast of 19.89 with wide, unchanging prediction intervals, reflecting the lack of significant autocorrelation, trend, or seasonality, and thus, the historical mean as the best predictor, this high uncertainty is crucial for reporting, emphasizing the stationary nature and absence of time-based patterns.

For Pipe2, operational recommendations include setting alert thresholds that accommodate its wider natural flow variation, with 65+ signaling high flow and below 15 indicating potential blockages. Maintenance planning should anticipate higher variability compared to Pipe1, with inspections scheduled during predicted low-flow periods.

Further model improvement could explore external factors causing fluctuations and test for potential weekly seasonality. A sample reporting language for Pipe1 highlights its stationary behavior (ADF p=0.01, KPSS p=0.1), the constant forecasts from the simple mean model (19.89), and the wide prediction intervals (11.58-28.20 at 95% confidence), underscoring the system’s stability alongside substantial natural variability.

ARIMA(0,0,0) model with a non-zero mean for Pipe 1

The residual analysis indicates that the ARIMA(0,0,0) model with a non-zero mean effectively fits the data. The residuals, depicted in the top panel, are randomly scattered around zero with constant variance and lack discernible patterns, suggesting the model successfully accounted for systematic data trends. The ACF plot in the bottom left confirms the absence of significant autocorrelation, further validating the model’s ability to reduce the data to random noise.

The histogram and density plot in the bottom right show a normal distribution centered around zero, without skewness or outliers, reinforcing the model’s appropriateness. Overall, the residuals’ characteristics of white noise—randomness, lack of autocorrelation, and normal distribution—confirm the model’s accuracy.

ARIMA(1,0,1) model with a non-zero mean for Pipe 2

The ARIMA(1,0,1) model with a non-zero mean shows some concerning signs. While the residuals are generally randomly scattered and normally distributed, the significant autocorrelation at lag 1 in the ACF plot is a major issue. This suggests that the model hasn’t fully captured the temporal dependencies in the data.

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.