Objective

This project has three parts. Part A forecasts daily cash withdrawals from four ATM machines for May 2010. Part B forecasts monthly residential power consumption for 2014. Part C, the bonus, takes two water flow data sets recorded on different time stamps, aggregates them to hourly means, tests whether they are stationary, and forecasts a week forward.

Each part ends with a forecast written to an Excel readable file, and the whole report is published to RPubs with the .Rmd submitted alongside it.

Pre-Requisite

Unlike the homework assignments, the data here does not come from a package. Four Excel files are required and they must sit in a data/ folder next to this .Rmd:

  1. ATM624Data.xlsx
  2. ResidentialCustomerForecastLoad-624.xlsx
  3. Waterflow_Pipe1.xlsx
  4. Waterflow_Pipe2.xlsx

Three packages are used directly. fpp3 provides the tsibble, feasts and fable functions that do all of the modelling. readxl reads the Excel files in. writexl writes the forecasts back out to Excel without needing a Java installation.

A fourth package, urca, must be installed even though it is never attached. feasts calls it for the unit root tests, and ARIMA() uses those tests internally when it selects orders automatically. Without urca present every ARIMA() call returns a NULL model with only a warning rather than an error, which is easy to mistake for a modelling problem. install.packages("urca") fixes it.

A note on method

Before starting I want to state the approach I follow in all three parts, because it is the same one each time and I would rather explain it once.

I plot the series first and let the plot tell me what structure is there. Then I deal with missing values and outliers explicitly rather than letting a model silently absorb them, because an outlier left in place will widen every prediction interval that follows it. Then I fit several candidate models rather than one, and I choose between them on a test set that the models never saw during fitting. Finally I check the residuals of the winner, because a model that leaves visible structure in its residuals is still leaving information on the table.

The one thing I deliberately do not do is pick a model on training accuracy. A model always fits the data it was fitted to better than it fits new data, so a training number measures memory, not forecasting skill.

Part A: ATM Cash Withdrawals

Reading and reshaping the data

The file holds a single long table: one row per ATM per day, with Cash in hundreds of dollars. The dates arrive as Excel serial numbers so they need converting, and Excel’s day zero is 1899-12-30.

atm_raw <- read_excel("data/ATM624Data.xlsx")
atm_raw$DATE <- as.Date(atm_raw$DATE, origin = "1899-12-30")

glimpse(atm_raw)
## Rows: 1,474
## Columns: 3
## $ DATE <date> 2009-05-01, 2009-05-01, 2009-05-02, 2009-05-02, 2009-05-03, 2009…
## $ 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, …

Two different kinds of missingness are present and they mean different things.

# Rows with no ATM label at all
atm_raw |> filter(is.na(ATM)) |> summarise(n = n(),
                                           from = min(DATE), to = max(DATE))
## # A tibble: 1 × 3
##       n from       to        
##   <int> <date>     <date>    
## 1    14 2010-05-01 2010-05-14
# Rows that name an ATM but have no Cash value
atm_raw |> filter(!is.na(ATM), is.na(Cash))
## # A tibble: 5 × 3
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2009-06-13 ATM1     NA
## 2 2009-06-16 ATM1     NA
## 3 2009-06-18 ATM2     NA
## 4 2009-06-22 ATM1     NA
## 5 2009-06-24 ATM2     NA

The fourteen rows with no ATM label are all dated May 2010. They are placeholder rows for the period I am being asked to forecast, so they carry no information and I drop them. The five rows that name an ATM but have no Cash are genuine gaps in the recorded history and I will fill them.

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

atm |> as_tibble() |> group_by(ATM) |>
  summarise(n = n(), first = min(DATE), last = max(DATE),
            missing = sum(is.na(Cash)),
            mean = round(mean(Cash, na.rm = TRUE), 1),
            median = round(median(Cash, na.rm = TRUE), 1),
            max = round(max(Cash, na.rm = TRUE), 1))
## # A tibble: 4 × 8
##   ATM       n first      last       missing  mean median    max
##   <chr> <int> <date>     <date>       <int> <dbl>  <dbl>  <dbl>
## 1 ATM1    365 2009-05-01 2010-04-30       3  83.9    91    180 
## 2 ATM2    365 2009-05-01 2010-04-30       2  62.6    67    147 
## 3 ATM3    365 2009-05-01 2010-04-30       0   0.7     0     96 
## 4 ATM4    365 2009-05-01 2010-04-30       0 474     404. 10920.

Each machine has exactly 365 daily observations running from 1 May 2009 to 30 April 2010, so the series are complete in time and I am forecasting the 31 days of May 2010. The summary already shows that the four machines are not alike. ATM4 averages 474 hundred dollars against ATM1’s 84, and its maximum of 10,920 is more than twenty times its own median. ATM3 averages 0.72, which is close enough to zero to be worth a separate look.

Visualising the four series

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

The plot settles most of the questions.

ATM1 and ATM2 are well behaved. Both fluctuate around a stable level with a dense high frequency wobble that is almost certainly a weekly cycle.

ATM3 is a flat line at zero until the very end, where three non-zero values appear.

ATM4 has one spike so large that it flattens everything else in its panel into a band near the bottom.

ATM3 is a different problem

atm |> filter(ATM == "ATM3", Cash > 0)
## # A tsibble: 3 x 3 [1D]
## # Key:       ATM [1]
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2010-04-28 ATM3     96
## 2 2010-04-29 ATM3     82
## 3 2010-04-30 ATM3     85

ATM3 has three non-zero days in its entire history, and they are the last three days of the record: 28, 29 and 30 April 2010. This is not a machine with missing data or a broken sensor. This is a machine that was switched on at the end of April 2010 and has been running for three days.

That distinction drives the whole treatment. If I feed 362 zeros and 3 real values into any seasonal model, the model will conclude the series is essentially zero and forecast close to zero for May. That forecast would be confidently wrong, because the machine is demonstrably now dispensing cash at a rate comparable to ATM1. The three real observations are the only relevant data, and the honest forecast is their mean, carried forward with a wide interval that admits how little I know. I return to this after handling the other three machines.

Handling ATM4’s outlier

atm |> filter(ATM == "ATM4") |> as_tibble() |>
  slice_max(Cash, n = 5) |> arrange(desc(Cash))
## # A tibble: 5 × 3
##   DATE       ATM     Cash
##   <date>     <chr>  <dbl>
## 1 2010-02-09 ATM4  10920.
## 2 2009-09-22 ATM4   1712.
## 3 2010-01-29 ATM4   1575.
## 4 2009-09-04 ATM4   1495.
## 5 2009-06-09 ATM4   1484.

A single day, 9 February 2010, records 10,919.76 hundred dollars, which is roughly $1.09 million out of one ATM in one day. The next largest day in the whole year is 1,712, and the median is 404. This is not a busy day, it is a data error or a one-off event such as a bulk reload, and either way it will not repeat in May.

Leaving it in would drag the level of any fitted model upward and inflate the residual variance, which widens every prediction interval. I replace it with the series median, which is a conservative choice because it does not assume the day was anything special.

atm <- atm |>
  mutate(Cash = if_else(ATM == "ATM4" & Cash > 5000,
                        median(atm$Cash[atm$ATM == "ATM4"], na.rm = TRUE),
                        Cash))

Filling the five missing values

The five gaps sit in ATM1 and ATM2 in June 2009. Because I have already seen a weekly rhythm, replacing a missing Wednesday with the overall series mean would distort it. I fill each gap with the median of that machine on that same day of the week, which respects the weekly structure.

atm <- atm |>
  mutate(dow = wday(DATE, label = TRUE)) |>
  group_by(ATM, dow) |>
  mutate(Cash = if_else(is.na(Cash), median(Cash, na.rm = TRUE), Cash)) |>
  ungroup() |>
  select(-dow)

sum(is.na(atm$Cash))
## [1] 0

No missing values remain.

Confirming the weekly season

atm |>
  filter(ATM != "ATM3") |>
  gg_season(Cash, period = "week") +
  facet_wrap(~ATM, scales = "free_y", ncol = 1) +
  labs(title = "Weekly seasonal pattern", y = "Cash (hundreds)")

atm |> as_tibble() |>
  filter(ATM != "ATM3") |>
  mutate(dow = wday(DATE, label = TRUE)) |>
  group_by(ATM, dow) |>
  summarise(mean = round(mean(Cash), 1), .groups = "drop") |>
  pivot_wider(names_from = ATM, values_from = mean)
## # A tibble: 7 × 4
##   dow    ATM1  ATM2  ATM4
##   <ord> <dbl> <dbl> <dbl>
## 1 Sun   103.   67.2  539.
## 2 Mon    86    58.7  481.
## 3 Tue    89.8  73.2  445 
## 4 Wed    82.2  43.4  414.
## 5 Thu    31.7  25.2  170 
## 6 Fri    98.6  92    574.
## 7 Sat    96.6  76    492.

The weekly pattern is real and it is the same shape on all three machines: averaged over the whole year, Thursday is far lower than every other day, at roughly a third of the level of the surrounding days. Whatever causes it is a property of the location rather than the machine. A model for these series has to carry a period-7 season.

A year-long average can hide a pattern that changes, so before relying on it I check which day is lowest month by month, and compare the most recent eight weeks with the rest of the year. This matters for forecasting: a model that adapts its seasonal shape will follow the recent pattern, while one that averages over the whole year will not.

low_by_month <- atm |> as_tibble() |>
  filter(ATM == "ATM1") |>
  mutate(month = yearmonth(DATE), dow = wday(DATE, label = TRUE, abbr = FALSE)) |>
  group_by(month, dow) |>
  summarise(mean = mean(Cash), .groups = "drop") |>
  group_by(month) |>
  slice_min(mean, n = 1, with_ties = FALSE) |>
  ungroup()

low_by_month
## # A tibble: 12 × 3
##       month dow       mean
##       <mth> <ord>    <dbl>
##  1 2009 May Thursday  11  
##  2 2009 Jun Thursday  16  
##  3 2009 Jul Thursday  13.6
##  4 2009 Aug Thursday  21  
##  5 2009 Sep Thursday  14.2
##  6 2009 Oct Thursday  10.8
##  7 2009 Nov Thursday  16  
##  8 2009 Dec Thursday  12.4
##  9 2010 Jan Thursday  19.2
## 10 2010 Feb Monday    74.2
## 11 2010 Mar Tuesday    6.8
## 12 2010 Apr Tuesday    3.5
low_recent_vs_earlier <- atm |> as_tibble() |>
  filter(ATM != "ATM3") |>
  mutate(dow = wday(DATE, label = TRUE, abbr = FALSE),
         period = if_else(DATE > max(DATE) - 56, "Last 8 weeks", "Earlier")) |>
  group_by(ATM, period, dow) |>
  summarise(mean = round(mean(Cash), 1), .groups = "drop") |>
  group_by(ATM, period) |>
  slice_min(mean, n = 1, with_ties = FALSE) |>
  ungroup()

low_recent_vs_earlier
## # A tibble: 6 × 4
##   ATM   period       dow       mean
##   <chr> <chr>        <ord>    <dbl>
## 1 ATM1  Earlier      Thursday  22  
## 2 ATM1  Last 8 weeks Tuesday    5.5
## 3 ATM2  Earlier      Thursday  13.2
## 4 ATM2  Last 8 weeks Tuesday    1.2
## 5 ATM4  Earlier      Thursday 120. 
## 6 ATM4  Last 8 weeks Tuesday   21.3
low_earlier <- low_recent_vs_earlier |>
  filter(ATM == "ATM1", period == "Earlier") |> pull(dow) |> as.character()
low_recent <- low_recent_vs_earlier |>
  filter(ATM == "ATM1", period == "Last 8 weeks") |> pull(dow) |> as.character()

For ATM1, the lowest day over the earlier part of the year is Thursday and the lowest day over the last eight weeks is Tuesday, and the same switch shows up on ATM2 and ATM4. The month by month table dates the change: Thursday is the low day in every month from May 2009 to January 2010, February 2010 has no clearly low day, and from March 2010 onward the dip falls on Tuesday. The low day has moved. That is an important finding for the forecast: May will follow the recent pattern, so I want models that let the seasonal shape adapt over time (ETS with a seasonal smoothing parameter, or STL with a finite seasonal window) rather than ones that fix a single shape from the whole year.

Do the series need transforming or differencing?

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

The Guerrero values come out at about 0.25 for ATM1, 0.67 for ATM2 and 0.45 for ATM4. All three sit between 0 and 1, which says the variance does grow somewhat with the level, but none is close enough to 0 to make a log transformation the obvious move.

I choose not to transform these series, for two reasons. ATM2 contains an exact zero, which rules out a log and makes small lambda Box-Cox values awkward. More importantly, ETS can handle growing variance directly through a multiplicative error term, so the structure the Guerrero value is pointing at can be captured inside the model rather than by reshaping the data first. That keeps the forecasts on the original scale, in hundreds of dollars, which matters for a figure a business reader has to act on.

atm |> filter(ATM != "ATM3") |>
  features(Cash, unitroot_nsdiffs)
## # A tibble: 3 × 2
##   ATM   nsdiffs
##   <chr>   <int>
## 1 ATM1        1
## 2 ATM2        1
## 3 ATM4        0
atm |> filter(ATM != "ATM3") |>
  features(Cash, unitroot_ndiffs)
## # A tibble: 3 × 2
##   ATM   ndiffs
##   <chr>  <int>
## 1 ATM1       1
## 2 ATM2       1
## 3 ATM4       0

The two tests split the machines. ATM1 and ATM2 each want one seasonal difference and one first difference, while ATM4 wants neither.

That difference is informative rather than awkward. ATM1 and ATM2 have enough persistence in their weekly pattern that differencing helps; ATM4, once its outlier is gone, is already stationary around a fixed level. It is another sign that ATM4 is the least structured of the three, which the accuracy table below confirms. I do not difference by hand: ARIMA applies whatever differencing it needs internally, and these tests simply tell me what to expect it to choose.

Fitting and comparing models

I hold back the last 28 days, four complete weeks, as a test set. Four weeks is long enough to expose a model that has failed to learn the weekly cycle and short enough to leave nearly a year for training.

I fit four candidates. Seasonal naive is the benchmark, and it is a serious one here because a series with a strong weekly cycle and no trend is exactly what it is built for, and because it simply repeats the most recent week it automatically picks up the move of the low day to Tuesday. ETS and ARIMA are chosen automatically. STL decomposition with an ETS model on the seasonally adjusted series is included because its default seasonal window lets the weekly shape change slowly over the year, which suits a pattern that has shifted.

atm_model <- atm |> filter(ATM != "ATM3")
train <- atm_model |> filter(DATE <= as.Date("2010-04-02"))

fit <- train |>
  model(
    snaive = SNAIVE(Cash),
    ets    = ETS(Cash),
    arima  = ARIMA(Cash),
    stl    = decomposition_model(STL(Cash ~ season(period = 7)),
                                 ETS(season_adjust ~ season("N")))
  )

fc_test <- fit |> forecast(h = 28)

accuracy(fc_test, atm_model) |>
  select(ATM, .model, RMSE, MAE, MAPE, MASE) |>
  arrange(ATM, RMSE)
## # A tibble: 12 × 6
##    ATM   .model  RMSE   MAE   MAPE  MASE
##    <chr> <chr>  <dbl> <dbl>  <dbl> <dbl>
##  1 ATM1  ets     12.4  10.0   70.8 0.544
##  2 ATM1  stl     12.6  10.1   82.1 0.550
##  3 ATM1  arima   13.0  10.5   93.1 0.570
##  4 ATM1  snaive  16.9  14.7   79.2 0.799
##  5 ATM2  stl     18.9  13.7  102.  0.660
##  6 ATM2  ets     19.8  14.2   56.5 0.681
##  7 ATM2  arima   21.7  16.9   51.2 0.812
##  8 ATM2  snaive  28.6  18.6   46.8 0.896
##  9 ATM4  stl    245.  204.   595.  0.579
## 10 ATM4  snaive 272.  196.   364.  0.556
## 11 ATM4  arima  281.  240.  1115.  0.680
## 12 ATM4  ets    313.  263.  1349.  0.747

Which measure I use, and which I ignore. I select on RMSE, because it penalises large misses more heavily than small ones, and a day when the machine runs out of cash is the costly mistake for a bank. I read MAE and MASE alongside it as checks. I deliberately ignore MAPE here. MAPE divides each error by the actual value, and these series have days with withdrawals close to zero, so a miss of a few hundred dollars on a near-zero day becomes an enormous percentage. That is why the MAPE column runs from about 50 percent to over 1,300 percent: it is measuring how close the actual values get to zero, not how good the forecasts are.

I also want to be precise about MASE, because it is easy to misread. MASE divides each model’s test error by the average error that a seasonal naive forecast made in-sample, one step ahead. A MASE below 1 therefore means the model’s errors over the four week test are smaller than seasonal naive’s one-step errors on the training data. It does not by itself mean the model beat seasonal naive on the test set. To answer that question I compare against the snaive row in the same table, which has its own MASE (0.80 on ATM1, for example).

ATM1. ETS is best on RMSE (12.4), MAE and MASE, with STL a close second and ARIMA third. All three beat seasonal naive (RMSE 16.9) clearly. ETS is the choice.

ATM2. STL is best on RMSE (18.9), MAE and MASE, with ETS second. Again all three beat seasonal naive (28.6) clearly. STL is the choice.

ATM4. This is the machine where the measures disagree, and the choice needs arguing rather than reading off. STL has the lowest RMSE (245 against 272 for seasonal naive), but seasonal naive has the lowest MAE (196 against 204) and the lowest MASE (0.556 against 0.579). ETS, which won on ATM1, comes last here on every measure, and ARIMA is third.

Read together, these numbers say STL and seasonal naive are about equally good on a typical day, and STL is better at avoiding the occasional large miss, which is exactly what RMSE rewards and MAE does not. I select STL for ATM4 on that basis, consistent with the RMSE rule I set out above, but the honest summary is that ATM4 carries little exploitable structure beyond its weekly cycle, and a forecast that simply repeats the last week would be nearly as good. The RMSE values for ATM4 are an order of magnitude larger than for the other machines only because the series itself is an order of magnitude larger.

I select the best model for each machine separately rather than forcing one model across all three, because the machines are genuinely different series and the table shows no single specification winning everywhere: ETS is first on ATM1 and last on ATM4.

A limitation of this comparison. The ranking rests on one four week window. A different four weeks could reorder close pairs such as ETS and STL on ATM1. Time series cross-validation, which repeats the comparison over many rolling windows, would give a more stable ranking; I did not do it here because it multiplies the fitting time by the number of windows, and on ATM1 and ATM2 the winners beat the benchmark by a wide margin that one window is enough to establish.

best <- accuracy(fc_test, atm_model) |>
  group_by(ATM) |>
  slice_min(RMSE, n = 1) |>
  ungroup() |>
  select(ATM, .model)

best
## # A tibble: 3 × 2
##   ATM   .model
##   <chr> <chr> 
## 1 ATM1  ets   
## 2 ATM2  stl   
## 3 ATM4  stl

Residual checks on the selected models

final_fit <- atm_model |>
  model(
    snaive = SNAIVE(Cash),
    ets    = ETS(Cash),
    arima  = ARIMA(Cash),
    stl    = decomposition_model(STL(Cash ~ season(period = 7)),
                                 ETS(season_adjust ~ season("N")))
  )

for (m in unique(best$ATM)) {
  chosen <- best$.model[best$ATM == m]
  print(final_fit |> filter(ATM == m) |> select(all_of(chosen)) |>
          gg_tsresiduals() + labs(title = paste(m, "-", chosen)))
}

final_fit |>
  augment() |>
  features(.innov, ljung_box, lag = 14, dof = 0) |>
  inner_join(best, by = c("ATM", ".model"))
## # A tibble: 3 × 4
##   ATM   .model lb_stat lb_pvalue
##   <chr> <chr>    <dbl>     <dbl>
## 1 ATM1  ets       28.8 0.0112   
## 2 ATM2  stl       22.5 0.0683   
## 3 ATM4  stl       46.4 0.0000245

The residual plots show series centred on zero with no remaining trend or seasonal shape, and the autocorrelation functions are mostly inside the significance bounds.

The Ljung-Box results are more mixed and I would rather report them plainly than gloss over them. ATM2 passes at the 5 percent level with a p-value near 0.07. ATM1 fails marginally at about 0.01. ATM4 fails decisively, with a p-value near 0.00002.

This means the residuals are not pure white noise, so some structure remains unmodelled. The practical consequence is specific and worth stating precisely: the point forecasts are still usable, because the test set comparison already showed these models matching or beating the benchmark on data they had not seen, but the prediction intervals are too narrow. Real uncertainty is wider than the intervals imply, especially for ATM4.

I chose not to chase the remaining structure with a more elaborate model. Daily cash withdrawal is driven by paydays, weekends and local events that no purely univariate model can see, so some autocorrelation in the residuals is expected here rather than a sign of a fixable error. The right response for a business forecast one month ahead is to treat the intervals as a lower bound on uncertainty, which is what I recommend in the conclusion.

Forecasting May 2010

fc_atm <- final_fit |>
  forecast(h = 31) |>
  inner_join(best, by = c("ATM", ".model"))

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

fc_low <- fc_atm |>
  as_tibble() |>
  mutate(dow = wday(DATE, label = TRUE, abbr = FALSE)) |>
  group_by(ATM, dow) |>
  summarise(mean_forecast = round(mean(.mean), 1), .groups = "drop") |>
  group_by(ATM) |>
  slice_min(mean_forecast, n = 1, with_ties = FALSE) |>
  ungroup()

fc_low
## # A tibble: 3 × 3
##   ATM   dow     mean_forecast
##   <chr> <ord>           <dbl>
## 1 ATM1  Tuesday           5.7
## 2 ATM2  Tuesday           1.8
## 3 ATM4  Tuesday          64.4
fc_low_day <- fc_low |> filter(ATM == "ATM1") |> pull(dow) |> as.character()

The forecasts carry the weekly cycle forward, and the intervals widen only gently because the series have no trend to extrapolate. The low day in the May forecast is Tuesday. That matches the recent pattern found earlier rather than the year-long average, which is what an adaptive model should do: the low day moved from Thursday to Tuesday in the final weeks of the data, and the selected models have followed it. A seasonal naive forecast would also have repeated the recent week, but a model that fixed one seasonal shape from the whole year would have put the dip on the wrong day for all of May.

Now ATM3, handled separately for the reasons given earlier.

atm3_vals <- atm |> filter(ATM == "ATM3", Cash > 0) |> pull(Cash)
atm3_mean <- mean(atm3_vals)
atm3_sd   <- sd(atm3_vals)

n3 <- length(atm3_vals)

# Option 1: interval from ATM3's own three values. A t multiplier with
# n - 1 = 2 degrees of freedom replaces 1.28, and sqrt(1 + 1/n) adds the
# uncertainty in the mean itself to the day-to-day variation.
half_own <- qt(0.90, df = n3 - 1) * atm3_sd * sqrt(1 + 1 / n3)

# Option 2: borrow day-to-day variability from ATM1, a machine at a similar
# level, by applying its coefficient of variation to ATM3's mean.
atm1_cash <- atm |> filter(ATM == "ATM1") |> pull(Cash)
atm1_cv   <- sd(atm1_cash) / mean(atm1_cash)
half_borrow <- qnorm(0.90) * atm1_cv * atm3_mean

half_used <- max(half_own, half_borrow)

c(mean = atm3_mean, sd_3_days = atm3_sd,
  half_width_own = half_own, half_width_borrowed = half_borrow)
##                mean           sd_3_days      half_width_own half_width_borrowed 
##           87.666667            7.371115           16.049307           48.864929
atm3_fc <- tibble(
  DATE = seq(as.Date("2010-05-01"), as.Date("2010-05-31"), by = "day"),
  ATM = "ATM3",
  .mean = atm3_mean,
  lower80 = pmax(0, atm3_mean - half_used),
  upper80 = atm3_mean + half_used
)

head(atm3_fc, 3)
## # A tibble: 3 × 5
##   DATE       ATM   .mean lower80 upper80
##   <date>     <chr> <dbl>   <dbl>   <dbl>
## 1 2010-05-01 ATM3   87.7    38.8    137.
## 2 2010-05-02 ATM3   87.7    38.8    137.
## 3 2010-05-03 ATM3   87.7    38.8    137.

With three observations I cannot estimate a weekly season or a trend, so the forecast is a flat line at their mean of about 87.7 hundred dollars.

The interval needs more thought than the point forecast. The three values happen to be close together (96, 82 and 85), so their standard deviation is small, about 7.4. Using 1.28 × sd directly would give an interval of about ±11 percent, narrower than ATM1’s, which would wrongly suggest I know more about a three day old machine than about one with a full year of history. A standard deviation estimated from three numbers is itself very uncertain, and three days cannot show the weekly swing that every other machine has.

So I compute the half-width two ways and use the wider. The first uses ATM3’s own values with a t multiplier on two degrees of freedom, which allows for the uncertainty in that small-sample standard deviation, giving about ±16. The second borrows the day-to-day variability of ATM1, a machine operating at a similar level, by applying ATM1’s coefficient of variation (about 0.43) to ATM3’s mean, giving about ±48.9. The 80 percent interval I report is 38.8 to 136.5 hundred dollars a day. That width is the correct signal to send: this number rests on three days of evidence and should be revisited as soon as more arrive.

One thing I did not do is give ATM3 a weekly pattern by borrowing ATM1’s seasonal shape. That would be a reasonable refinement, but three days give no evidence that ATM3 shares ATM1’s pattern. Its three days fall on a Wednesday, Thursday and Friday, so they include no Tuesday, the day the other machines now dip, and cannot tell me whether ATM3 dips too. I keep the flat forecast and let the wide interval carry that uncertainty, which is another reason the borrowed, wider interval is the right one to use.

Part A forecast output

atm_out <- fc_atm |>
  hilo(80) |>
  as_tibble() |>
  mutate(lower80 = pmax(0, `80%`$lower), upper80 = `80%`$upper) |>
  select(DATE, ATM, Forecast = .mean, lower80, upper80) |>
  bind_rows(atm3_fc |> select(DATE, ATM, Forecast = .mean, lower80, upper80)) |>
  arrange(ATM, DATE) |>
  mutate(DATE = as.character(DATE),
         across(c(Forecast, lower80, upper80), \(x) round(x, 2)))

write_xlsx(atm_out, "ATM_Forecast_May2010.xlsx")

atm_out |> group_by(ATM) |>
  summarise(days = n(),
            total_hundreds = round(sum(Forecast), 0),
            total_dollars = scales::comma(round(sum(Forecast) * 100, 0)))
## # A tibble: 4 × 4
##   ATM    days total_hundreds total_dollars
##   <chr> <int>          <dbl> <chr>        
## 1 ATM1     31           2388 238,808      
## 2 ATM2     31           1773 177,346      
## 3 ATM3     31           2718 271,777      
## 4 ATM4     31          13744 1,374,442
head(atm_out, 8)
## # A tibble: 8 × 5
##   DATE       ATM   Forecast lower80 upper80
##   <chr>      <chr>    <dbl>   <dbl>   <dbl>
## 1 2010-05-01 ATM1     87.0     56.3   118. 
## 2 2010-05-02 ATM1    101.      70.0   132. 
## 3 2010-05-03 ATM1     73.1     42.3   104. 
## 4 2010-05-04 ATM1      5.72     0      36.5
## 5 2010-05-05 ATM1    100.      69.3   131. 
## 6 2010-05-06 ATM1     79.4     48.6   110. 
## 7 2010-05-07 ATM1     85.6     53.0   118. 
## 8 2010-05-08 ATM1     87.0     54.5   120.

The file ATM_Forecast_May2010.xlsx contains all 124 rows, 31 days for each of the four machines, with an 80 percent interval alongside each point forecast. In business terms the four machines together are forecast to dispense on the order of 2.06 million dollars over May 2010, with ATM4 accounting for about two thirds of it on its own.

Part B: Residential Power Consumption

Reading the data

pwr_raw <- read_excel("data/ResidentialCustomerForecastLoad-624.xlsx")
glimpse(pwr_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…
pwr <- pwr_raw |>
  mutate(Month = yearmonth(`YYYY-MMM`)) |>
  select(Month, KWH) |>
  as_tsibble(index = Month)

pwr |> as_tibble() |>
  summarise(n = n(), from = min(Month), to = max(Month),
            missing = sum(is.na(KWH)))
## # A tibble: 1 × 4
##       n     from       to missing
##   <int>    <mth>    <mth>   <int>
## 1   192 1998 Jan 2013 Dec       1

192 monthly observations from January 1998 to December 2013, exactly the sixteen years described, with one missing value.

Plotting

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

Three features stand out. There is a slow upward trend across the sixteen years. There is a strong annual season with two peaks a year, a tall summer peak and a smaller winter one, which is the signature of air conditioning and heating load. And there is one dramatic downward spike near 2010 that breaks the pattern completely.

The gap and the spike

pwr |> filter(is.na(KWH))
## # A tsibble: 1 x 2 [1M]
##      Month   KWH
##      <mth> <dbl>
## 1 2008 Sep    NA
pwr |> as_tibble() |> slice_min(KWH, n = 3)
## # A tibble: 3 × 2
##      Month     KWH
##      <mth>   <dbl>
## 1 2010 Jul  770523
## 2 2005 Nov 4313019
## 3 1999 Dec 4419229

September 2008 is missing. July 2010 records 770,523 KWH when neighbouring Julys run between six and ten million. A summer month cannot plausibly be the lowest month in sixteen years of record; July is normally the annual peak. I treat it as a recording error rather than a real collapse in demand.

Both are the same kind of problem, a month whose value I do not trust, so I treat them the same way. I set the July 2010 value to missing and then interpolate both gaps from an ARIMA model fitted to the rest of the series, which fills them consistently with the trend and season around them instead of with a flat average.

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

# Orders are stated explicitly here rather than searched: automatic selection
# needs complete data to run its unit root tests, and the whole point of this
# step is that the data is not yet complete. A seasonal AR term with one
# seasonal difference is the right shape for a series like this.
pwr_clean <- pwr_flagged |>
  model(ARIMA(KWH ~ pdq(1, 0, 0) + PDQ(1, 1, 0))) |>
  interpolate(pwr_flagged)

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

Both months are now filled. The July 2010 figure of about 8.1 million sits naturally among its neighbouring Julys, which run 7.6 million in 2008 and 7.7 million in 2009, confirming the interpolation is sensible rather than merely non-missing.

pwr_clean |> filter(month(Month) == 7, year(Month) %in% 2008:2012)
## # A tsibble: 5 x 2 [1M]
##      Month       KWH
##      <mth>     <dbl>
## 1 2008 Jul  7643987 
## 2 2009 Jul  7713260 
## 3 2010 Jul  8114472.
## 4 2011 Jul 10093343 
## 5 2012 Jul  8886851
pwr_clean |> autoplot(KWH) +
  labs(title = "Monthly residential power consumption, cleaned",
       y = "KWH", x = NULL)

Decomposition

pwr_clean |>
  model(STL(KWH ~ trend(window = 21) + season(window = "periodic"))) |>
  components() |>
  autoplot() +
  labs(title = "STL decomposition of power consumption")

The decomposition separates the three pieces cleanly. The trend rises steadily with a flattening after about 2010. The seasonal component is stable in shape and large in magnitude relative to the remainder, confirming that the annual cycle is the dominant structure. The remainder is small and mostly patternless, which says trend and season together explain most of what happens in this series.

Transformation

pwr_clean |> features(KWH, features = guerrero)
## # A tibble: 1 × 1
##   lambda_guerrero
##             <dbl>
## 1          -0.212

The Guerrero method suggests a lambda of about -0.21. Looking back at the series plot, the seasonal swings do grow slightly as the level rises over the sixteen years, so a variance stabilising transformation is justified here in a way it was not for the ATM data.

I use a log transformation, which is Box-Cox with lambda = 0, rather than the suggested -0.21. The two values are close, so they stabilise the variance by almost the same amount, and the Guerrero estimate is itself uncertain at this sample size, so I do not treat -0.21 as more precise than “close to zero”. Given that, I prefer the transformation that is easier to explain: a log models proportional change, so the seasonal swing becomes a percentage of the level rather than a fixed number of kilowatt hours, which is how a business reader would naturally describe it.

fable back transforms the forecasts automatically. It also bias-adjusts them: the .mean it reports for a log model is the mean of the back-transformed forecast distribution, not simply exp() of the forecast on the log scale (which would be the median and would sit slightly low). The prediction intervals are percentiles of that same distribution, so no further adjustment is needed.

pwr_clean |> features(log(KWH), unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
pwr_clean |> features(log(KWH), unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

Both tests return 1: one seasonal difference to remove the annual cycle and one first difference to remove the trend. That matches what the plot and the decomposition already showed, a series with both a season and a trend.

It is worth noting that the ARIMA model selected later does not follow this exactly. It takes the seasonal difference but handles the trend with a drift term instead of a first difference, which is an equivalent way of accounting for steady growth and often a better behaved one for forecasting. These tests guide the search rather than dictate its result.

Fitting and comparing models

I hold back the last 24 months, two full years, as a test set. Two years is the right length for a monthly series with an annual cycle: it tests the seasonal pattern twice rather than once.

pwr_train <- pwr_clean |> filter(Month <= yearmonth("2011 Dec"))

pwr_fit <- pwr_train |>
  model(
    snaive = SNAIVE(KWH),
    ets    = ETS(log(KWH)),
    arima  = ARIMA(log(KWH)),
    stl    = decomposition_model(STL(log(KWH)),
                                 ETS(season_adjust ~ season("N")))
  )

pwr_fit
## # A mable: 1 x 4
##     snaive          ets                              arima
##    <model>      <model>                            <model>
## 1 <SNAIVE> <ETS(A,N,A)> <ARIMA(1,0,1)(2,1,0)[12] w/ drift>
## # ℹ 1 more variable: stl <model>
pwr_fc_test <- pwr_fit |> forecast(h = 24)

accuracy(pwr_fc_test, pwr_clean) |>
  select(.model, RMSE, MAE, MAPE, MASE) |>
  arrange(RMSE)
## # A tibble: 4 × 5
##   .model    RMSE     MAE  MAPE  MASE
##   <chr>    <dbl>   <dbl> <dbl> <dbl>
## 1 arima  785721. 626777.  8.15  1.03
## 2 ets    853626. 648947.  8.65  1.07
## 3 stl    882970. 735436.  9.90  1.21
## 4 snaive 965567. 742063.  9.30  1.22

The training-period models chosen automatically were ETS(A,N,A) on the log scale, which is additive level and season with no trend, and ARIMA(1,0,1)(2,1,0)[12] with drift.

ARIMA is best on every measure: RMSE (786 thousand KWH), MAE, MAPE (8.2 percent) and MASE (1.03). ETS is second on every measure. STL is third on RMSE but slightly worse than seasonal naive on MAPE (9.9 against 9.3 percent), so it is not a clear improvement over the benchmark. Seasonal naive is last on RMSE and MASE.

Why I chose ARIMA and not the others. ARIMA wins on all four measures, so the choice does not depend on which measure I prefer. ETS is the runner-up and a reasonable alternative; its weakness here is that it chose no trend term, so it cannot follow the slow rise in consumption, while ARIMA’s drift term can. STL does not reliably beat the benchmark, and seasonal naive, by construction, only repeats last year and so ignores the trend entirely.

Reading MASE correctly. Every MASE here is above 1, including the winner’s 1.03. That does not mean ARIMA failed to beat seasonal naive. MASE divides the test error by seasonal naive’s average in-sample, one-step error, and errors two years ahead are naturally larger than errors one step ahead. The fair comparison is with seasonal naive’s own MASE on the same test window, which is 1.22. ARIMA’s MAE is about 16 percent lower than seasonal naive’s, so it did beat the benchmark.

pwr_best <- accuracy(pwr_fc_test, pwr_clean) |>
  slice_min(RMSE, n = 1) |> pull(.model)
pwr_best
## [1] "arima"
pwr_final <- pwr_clean |>
  model(
    snaive = SNAIVE(KWH),
    ets    = ETS(log(KWH)),
    arima  = ARIMA(log(KWH)),
    stl    = decomposition_model(STL(log(KWH)),
                                 ETS(season_adjust ~ season("N")))
  )

pwr_final |> select(all_of(pwr_best)) |> report()
## Series: KWH 
## Model: ARIMA(0,0,1)(2,1,0)[12] w/ drift 
## Transformation: log(KWH) 
## 
## Coefficients:
##          ma1     sar1     sar2  constant
##       0.2602  -0.7115  -0.3561    0.0322
## s.e.  0.0846   0.0757   0.0776    0.0089
## 
## sigma^2 estimated as 0.008457:  log likelihood=172.75
## AIC=-335.49   AICc=-335.15   BIC=-319.53

The final model is refitted on all sixteen years, and when ARIMA() searches again on the full data it settles on a slightly different order from the one chosen on the training years: ARIMA(0,0,1)(2,1,0)[12] with drift rather than (1,0,1)(2,1,0)[12]. The seasonal part, the seasonal difference and the drift are unchanged, and these carry almost all of the forecast; only the short-run non-seasonal terms differ. I allow the re-search rather than freezing the training order, because the extra two years are real information and the test already established that the ARIMA family is the right choice for this series.

Residual check

pwr_final |> select(all_of(pwr_best)) |> gg_tsresiduals() +
  labs(title = paste("Residuals,", pwr_best))

pwr_final |> select(all_of(pwr_best)) |> augment() |>
  features(.innov, ljung_box, lag = 24, dof = 0)
## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 arima     25.7     0.366

The residuals are centred on zero, the histogram is roughly symmetric, and the autocorrelations sit inside the bounds at almost every lag. The Ljung-Box p-value of about 0.37 is comfortably above 0.05, so I fail to reject the hypothesis that the residuals are white noise.

This is a noticeably cleaner result than the ATM models produced, and it means the prediction intervals here can be taken at face value rather than treated as a lower bound. The model has taken the trend and the annual season out of the series and left behind something close to noise, which is exactly the outcome I want before forecasting.

Forecasting 2014

fc_pwr <- pwr_final |> select(all_of(pwr_best)) |> forecast(h = 12)

fc_pwr |>
  autoplot(pwr_clean |> filter(Month >= yearmonth("2008 Jan"))) +
  labs(title = "2014 monthly power consumption forecast",
       y = "KWH", x = NULL)

The forecast reproduces the familiar double peaked annual shape. The two highest months are January at about 10.3 million KWH and August at about 10.1 million, the winter heating peak and the summer cooling peak respectively, with the shoulder months of May and October forming the troughs near 5.9 million. The ratio between peak and trough is close to 1.75, so this is a series where the season matters more than the trend for any single month’s planning.

The 80 percent intervals stay at roughly ±12 percent of the forecast through the whole year. They widen in absolute terms in the high-consumption months simply because the level is higher, but in relative terms they barely grow, because most of the uncertainty a year ahead comes from the month-to-month noise rather than from uncertainty about the trend. Because the model was fitted on a log scale they are asymmetric on the original scale, slightly wider above the forecast than below it. That is the correct behaviour for a quantity that cannot go negative.

pwr_out <- fc_pwr |>
  hilo(80) |>
  as_tibble() |>
  mutate(Month = as.character(Month),
         Forecast_KWH = round(.mean, 0),
         lower80 = round(`80%`$lower, 0),
         upper80 = round(`80%`$upper, 0)) |>
  select(Month, Forecast_KWH, lower80, upper80)

write_xlsx(pwr_out, "Power_Forecast_2014.xlsx")
pwr_out
## # A tibble: 12 × 4
##    Month    Forecast_KWH lower80  upper80
##    <chr>           <dbl>   <dbl>    <dbl>
##  1 2014 Jan     10288648 9106320 11526799
##  2 2014 Feb      8500789 7492306  9558514
##  3 2014 Mar      6642050 5854077  7468498
##  4 2014 Apr      5976476 5267462  6720108
##  5 2014 May      5941318 5236475  6680576
##  6 2014 Jun      8261934 7281787  9289938
##  7 2014 Jul      9495766 8369245 10677292
##  8 2014 Aug     10056074 8863081 11307318
##  9 2014 Sep      8453019 7450203  9504799
## 10 2014 Oct      5849270 5155347  6577075
## 11 2014 Nov      6109819 5384986  6870043
## 12 2014 Dec      8273192 7291709  9302597

The total forecast consumption for 2014 is 93,848,355 KWH, and the monthly figures are written to Power_Forecast_2014.xlsx.

Part C (Bonus): Water Flow

Reading two differently stamped files

read_pipe <- function(path) {
  d <- read_excel(path, col_names = FALSE, skip = 1,
                  col_types = c("numeric", "numeric"))
  names(d) <- c("serial", "WaterFlow")
  d |> mutate(DateTime = as.POSIXct(serial * 86400,
                                    origin = "1899-12-30", tz = "UTC")) |>
    select(DateTime, WaterFlow)
}

pipe1 <- read_pipe("data/Waterflow_Pipe1.xlsx")
pipe2 <- read_pipe("data/Waterflow_Pipe2.xlsx")

bind_rows(
  pipe1 |> summarise(pipe = "Pipe1", n = n(), from = min(DateTime),
                     to = max(DateTime)),
  pipe2 |> summarise(pipe = "Pipe2", n = n(), from = min(DateTime),
                     to = max(DateTime))
)
## # A tibble: 2 × 4
##   pipe      n from                to                 
##   <chr> <int> <dttm>              <dttm>             
## 1 Pipe1  1000 2015-10-23 00:24:06 2015-11-01 23:35:43
## 2 Pipe2  1000 2015-10-23 01:00:00 2015-12-03 16:00:00

Both files hold 1000 readings, but over very different spans. Pipe 1’s thousand readings cover about ten days, so it records several times per hour. Pipe 2’s thousand readings cover about 41 days, which works out at roughly one reading per hour. This is the difference in time stamps the assignment refers to, and it is the reason the two series cannot simply be added together as they stand.

bind_rows(
  tibble(pipe = "Pipe1", gap_mins = as.numeric(diff(pipe1$DateTime), units = "mins")),
  tibble(pipe = "Pipe2", gap_mins = as.numeric(diff(pipe2$DateTime), units = "mins"))
) |>
  group_by(pipe) |>
  summarise(median_gap = round(median(gap_mins), 1),
            mean_gap = round(mean(gap_mins), 1))
## # A tibble: 2 × 3
##   pipe  median_gap mean_gap
##   <chr>      <dbl>    <dbl>
## 1 Pipe1         10     14.4
## 2 Pipe2         60     60

The gaps confirm it. Pipe 2 records exactly once an hour, on the hour. Pipe 1 is irregular: its median gap is 10 minutes and its mean gap about 14 minutes, so it averages roughly four readings an hour but with no fixed schedule. Aggregating both to an hourly mean puts them on a common footing, which is what the assignment asks for.

Time-base sequencing and hourly aggregation

For multiple recordings within an hour I take the mean, as instructed. The timestamp is floored to the hour so that every reading between 09:00:00 and 09:59:59 is assigned to the 09:00 bucket.

hourly <- function(d) {
  d |>
    mutate(Hour = floor_date(DateTime, "hour")) |>
    group_by(Hour) |>
    summarise(WaterFlow = mean(WaterFlow), n_readings = n(), .groups = "drop")
}

p1_hr <- hourly(pipe1)
p2_hr <- hourly(pipe2)

head(p1_hr, 10)
## # A tibble: 10 × 3
##    Hour                WaterFlow n_readings
##    <dttm>                  <dbl>      <int>
##  1 2015-10-23 00:00:00      26.1          4
##  2 2015-10-23 01:00:00      18.9          5
##  3 2015-10-23 02:00:00      15.2          2
##  4 2015-10-23 03:00:00      23.1          5
##  5 2015-10-23 04:00:00      15.5          2
##  6 2015-10-23 05:00:00      22.7          9
##  7 2015-10-23 06:00:00      20.6          8
##  8 2015-10-23 07:00:00      18.4          3
##  9 2015-10-23 08:00:00      20.8          3
## 10 2015-10-23 09:00:00      21.2          3
p1_hr |> summarise(min = min(n_readings), median = median(n_readings),
                   mean = round(mean(n_readings), 1), max = max(n_readings))
## # A tibble: 1 × 4
##     min median  mean   max
##   <int>  <dbl> <dbl> <int>
## 1     1      4   4.2    10

The n_readings column makes the aggregation visible: Pipe 1’s hours combine anywhere from a couple of readings to nearly ten into one mean, about four on average, while Pipe 2’s hours contain one reading each, so for Pipe 2 the “aggregation” changes nothing.

Why I keep the two pipes separate. The assignment can be read as asking for one combined series, for example total flow through both pipes. I analyse them separately for two reasons. First, they only overlap for about ten days, because Pipe 1 stops on 1 November while Pipe 2 runs to 3 December, so a combined series would either be ten days long or would have to invent Pipe 1’s values for the remaining month. Second, the two pipes turn out to behave differently (one is pure noise, the other has a weak daily cycle), and adding them together would blur exactly the finding that makes this part interesting.

p1_ts <- p1_hr |> select(Hour, WaterFlow) |> as_tsibble(index = Hour) |>
  fill_gaps() |>
  mutate(WaterFlow = if_else(is.na(WaterFlow),
                             mean(WaterFlow, na.rm = TRUE), WaterFlow))

p2_ts <- p2_hr |> select(Hour, WaterFlow) |> as_tsibble(index = Hour) |>
  fill_gaps() |>
  mutate(WaterFlow = if_else(is.na(WaterFlow),
                             mean(WaterFlow, na.rm = TRUE), WaterFlow))

c(pipe1_hours = nrow(p1_ts), pipe2_hours = nrow(p2_ts))
## pipe1_hours pipe2_hours 
##         240        1000

fill_gaps() inserts any hour that had no reading at all so that the series is regularly spaced, which the modelling functions require.

bind_rows(p1_ts |> mutate(pipe = "Pipe 1") |> as_tibble(),
          p2_ts |> mutate(pipe = "Pipe 2") |> as_tibble()) |>
  ggplot(aes(Hour, WaterFlow)) +
  geom_line() +
  facet_wrap(~pipe, scales = "free", ncol = 1) +
  labs(title = "Hourly mean water flow", y = "Water flow", x = NULL)

Pipe 1 looks like noise around a constant level near 20 with no visible trend. Pipe 2 is noisier and sits at a higher level near 40, again without an obvious trend or repeating daily shape.

Are the series stationary?

This is the question the assignment asks directly, so I test it rather than judge it by eye. unitroot_kpss tests the null hypothesis that the series is stationary, so a p-value above 0.05 means I fail to reject stationarity.

bind_rows(
  p1_ts |> features(WaterFlow, unitroot_kpss) |> mutate(pipe = "Pipe 1"),
  p2_ts |> features(WaterFlow, unitroot_kpss) |> mutate(pipe = "Pipe 2")
) |> select(pipe, kpss_stat, kpss_pvalue)
## # A tibble: 2 × 3
##   pipe   kpss_stat kpss_pvalue
##   <chr>      <dbl>       <dbl>
## 1 Pipe 1     0.244      0.1   
## 2 Pipe 2     0.438      0.0606
bind_rows(
  p1_ts |> features(WaterFlow, unitroot_ndiffs) |> mutate(pipe = "Pipe 1"),
  p2_ts |> features(WaterFlow, unitroot_ndiffs) |> mutate(pipe = "Pipe 2")
) |> select(pipe, ndiffs)
## # A tibble: 2 × 2
##   pipe   ndiffs
##   <chr>   <int>
## 1 Pipe 1      0
## 2 Pipe 2      0

Pipe 1 returns a KPSS p-value of 0.1, which is the ceiling the test reports, so there is no evidence at all against stationarity. Pipe 2 returns about 0.06, which clears the 0.05 threshold but only just, and I would rather flag that than round it up into a cleaner story than the data supports.

The unitroot_ndiffs result settles it. Both series need zero differences, and that function applies the KPSS test repeatedly rather than relying on one borderline p-value. Both series are stationary as they stand, with Pipe 2 the more marginal of the two.

That answers the first half of the assignment’s question. It also points at something important for the second half: a stationary series has no trend to extrapolate, so whatever forecasting value exists here has to come from autocorrelation instead.

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

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

acf_summary <- function(ts) {
  a <- ts |> ACF(WaterFlow, lag_max = 48)
  bound <- 1.96 / sqrt(nrow(ts))
  tibble(spikes_outside = sum(abs(a$acf) > bound),
         bound = round(bound, 3),
         lag24 = round(a$acf[24], 3))
}

acf_p1 <- acf_summary(p1_ts)
acf_p2 <- acf_summary(p2_ts)
bind_rows(acf_p1 |> mutate(pipe = "Pipe 1"), acf_p2 |> mutate(pipe = "Pipe 2"))
## # A tibble: 2 × 4
##   spikes_outside bound lag24 pipe  
##            <int> <dbl> <dbl> <chr> 
## 1              1 0.127 0.032 Pipe 1
## 2              7 0.062 0.118 Pipe 2

The autocorrelation plots are the crux of it, and they separate the two pipes. With 48 lags each tested at the 5 percent level, about two or three spikes are expected to cross the bound by chance alone even in pure noise.

For Pipe 1, 1 of 48 spikes cross the bound, which is no more than chance would produce, so it is evidence of nothing at all. Lag 24, where a daily cycle would appear, is 0.032, against a bound of ±0.127.

For Pipe 2, 7 of 48 spikes cross the bound, and lag 24 is 0.118 against a bound of ±0.062. That is more than chance would produce, so there is a weak daily signal here. It is small in absolute terms, but it is real.

I confirm this with a formal test rather than leaving it to the eye. The Ljung-Box test asks whether the autocorrelations up to a given lag are jointly different from zero, which is exactly the question the plots raise.

bind_rows(
  p1_ts |> features(WaterFlow, ljung_box, lag = 24) |> mutate(pipe = "Pipe 1"),
  p2_ts |> features(WaterFlow, ljung_box, lag = 24) |> mutate(pipe = "Pipe 2")
) |> select(pipe, lb_stat, lb_pvalue)
## # A tibble: 2 × 3
##   pipe   lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 Pipe 1    20.7   0.655  
## 2 Pipe 2    48.7   0.00205

The test agrees with the plots and makes the distinction sharp. Pipe 1 returns a p-value near 0.66, so I cannot reject white noise: this series is indistinguishable from random variation about a constant. Pipe 2 returns about 0.002, so I can reject it: there is genuine autocorrelation structure, however weak.

Can they be forecast?

The answer differs by pipe, which is the most interesting finding in this part.

A forecast is only informative when the past tells you something about the future. Both series are stationary, so neither has a trend to extrapolate. That leaves autocorrelation as the only possible source of forecasting power, and the tests above show Pipe 1 has none while Pipe 2 has a little.

On the diagnostics alone, my expectation was that Pipe 1 cannot be forecast in any useful sense, since it is white noise and the optimal forecast is its mean, while Pipe 2 might be forecast a little better than its mean, since a weak daily cycle exists that a model could try to exploit.

These are claims about forecasting, so I test them on a held out test set rather than asserting them from the diagnostics. The test window has to be long enough to test a daily cycle more than once. For Pipe 2, which has about six weeks of data, I hold back a full week (168 hours), the same length as the forecast I am asked to produce. Pipe 1 has only ten days, so a week would leave too little to train on; I hold back two days (48 hours) there.

I compare four candidates. The mean is the natural benchmark for a stationary series. Naive (repeat the last value) is the benchmark for a random walk, and is included to confirm that the series is not one. ARIMA and ETS are chosen automatically. I do not include seasonal naive or STL with a daily season, since neither pipe shows a daily cycle strong enough to justify a model built around one.

fit_pipe <- function(ts, h) {
  n <- nrow(ts)
  tr <- ts |> slice(1:(n - h))
  f <- tr |> model(mean = MEAN(WaterFlow),
                   naive = NAIVE(WaterFlow),
                   arima = ARIMA(WaterFlow),
                   ets = ETS(WaterFlow))
  accuracy(f |> forecast(h = h), ts) |> select(.model, RMSE, MAE, MASE) |>
    arrange(RMSE)
}

acc_p1 <- fit_pipe(p1_ts, h = 48)
acc_p2 <- fit_pipe(p2_ts, h = 168)
acc_p1
## # A tibble: 4 × 4
##   .model  RMSE   MAE  MASE
##   <chr>  <dbl> <dbl> <dbl>
## 1 ets     3.87  2.88 0.607
## 2 arima   3.87  2.88 0.608
## 3 mean    3.87  2.88 0.608
## 4 naive   6.09  5.07 1.07
acc_p2
## # A tibble: 4 × 4
##   .model  RMSE   MAE  MASE
##   <chr>  <dbl> <dbl> <dbl>
## 1 mean    11.3  7.38 0.757
## 2 ets     11.9  8.68 0.890
## 3 naive   16.9 14.6  1.49 
## 4 arima   22.6 20.1  2.06
p2_get <- function(m, col) acc_p2 |> filter(.model == m) |> pull({{ col }})
p2_arima_rmse <- p2_get("arima", RMSE); p2_mean_rmse <- p2_get("mean", RMSE)
p2_arima_mae  <- p2_get("arima", MAE);  p2_mean_mae  <- p2_get("mean", MAE)
p2_rmse_gain  <- 100 * (p2_mean_rmse - p2_arima_rmse) / p2_mean_rmse
p1_get <- function(m) acc_p1 |> filter(.model == m) |> pull(RMSE)
p1_best <- acc_p1$.model[1]

On Pipe 1 the mean has a test RMSE of 3.87, ARIMA 3.87 and naive 6.09. The mean and ARIMA are identical because ARIMA has selected the mean model (shown below). Naive is worse than the mean, confirming the series is not a random walk. Nothing beats the mean.

On Pipe 2, ARIMA’s test RMSE is 22.55 against 11.33 for the mean, about 99.1 percent higher, and its MAE is 20.05 against 7.38. This is the opposite of what I expected, and it is the most important result in Part C. The mean is the best model on every measure, ETS is a close second, and ARIMA is not just worse but the worst of the four, behind even naive.

To understand why, I look at which model ARIMA() actually chose on the training window, since that is the model that produced the test forecasts.

p2_train <- p2_ts |> slice(1:(nrow(p2_ts) - 168))
p2_train_fit <- p2_train |> model(arima = ARIMA(WaterFlow))
report(p2_train_fit)
## Series: WaterFlow 
## Model: ARIMA(1,1,0)(1,0,0)[24] w/ drift 
## 
## Coefficients:
##           ar1    sar1  constant
##       -0.5187  0.1598   -0.0171
## s.e.   0.0296  0.0344    0.4789
## 
## sigma^2 estimated as 193.2:  log likelihood=-3365.16
## AIC=6738.33   AICc=6738.38   BIC=6757.22
model_name <- function(fit) {
  sub("Model: ", "", grep("^Model:", capture.output(report(fit)), value = TRUE))
}
p2_train_model <- model_name(p2_train_fit)
p2_full_model  <- model_name(p2_ts |> model(arima = ARIMA(WaterFlow)))

p2_train_fit |>
  forecast(h = 168) |>
  autoplot(p2_ts |> slice((nrow(p2_ts) - 400):nrow(p2_ts)), level = 80) +
  labs(title = "Pipe 2: ARIMA chosen on the training window, forecast over the test week",
       y = "Water flow", x = NULL)

On the first five weeks of data, the automatic search picks ARIMA(1,1,0)(1,0,0)[24] w/ drift , which is not the model it picks on the full six weeks (ARIMA(0,0,1)(0,0,1)[24] w/ mean, shown below). The lesson is about how fragile automatic order selection is on a noisy series with only weak structure. The search compares models whose AICc values differ by very little, so small changes in the data can flip it to a different specification, and one that extrapolates something that is not really there produces forecasts far worse than the plain mean. A single week of extra data was enough to change the answer.

The decision for Pipe 2 is therefore the mean, not ARIMA. The mean is best on the test week on every measure, it needs no order selection that could go wrong, and its forecast is the same flat line that the full-data ARIMA would converge to after its first day anyway. The weak daily cycle is statistically real, which the Ljung-Box test and the seasonal MA term both confirm, but it is too weak to improve a week-ahead forecast, and trying to exploit it adds risk without adding accuracy.

Neither result is a failure of the modelling. Getting a flat forecast from a series with no usable structure is the correct answer, and recognising that is more useful than forcing a more elaborate model onto data that cannot support one. Had I chosen Pipe 2’s model from the diagnostics alone, I would have picked ARIMA; the test set is what caught the mistake.

p1_ts |> model(ARIMA(WaterFlow)) |> report()
## Series: WaterFlow 
## Model: ARIMA(0,0,0) w/ mean 
## 
## Coefficients:
##       constant
##        19.8927
## s.e.    0.2708
## 
## sigma^2 estimated as 17.68:  log likelihood=-684.72
## AIC=1373.44   AICc=1373.49   BIC=1380.4
p2_ts |> model(ARIMA(WaterFlow)) |> report()
## Series: WaterFlow 
## Model: ARIMA(0,0,1)(0,0,1)[24] w/ mean 
## 
## Coefficients:
##           ma1    sma1  constant
##       -0.0464  0.1172   39.7950
## s.e.   0.0311  0.0309    0.3764
## 
## sigma^2 estimated as 125.8:  log likelihood=-3835.06
## AIC=7678.13   AICc=7678.17   BIC=7697.76

The selected ARIMA orders tell the story precisely.

Pipe 1 gives ARIMA(0,0,0) w/ mean. Every order is zero, which means the search examined the data and concluded there was nothing to model: no AR terms, no MA terms, no differencing. The model is a constant plus white noise, and the constant is the series mean of about 19.9. This is the algorithm independently reaching the same conclusion the ACF plot and the Ljung-Box test did.

Pipe 2 gives ARIMA(0,0,1)(0,0,1)[24] w/ mean. The seasonal MA term at period 24 is the daily structure the ACF picked up, and it is the one piece of genuine signal in either pipe. Its coefficient is small, around 0.12, which is why it would barely move a forecast away from the mean. I report this model to show what structure exists, not to forecast with it, for the reasons given above. The non-seasonal MA term, about -0.05 with a standard error of 0.03, is not significantly different from zero; the search kept it because it lowered the AICc slightly, not because it captures anything real.

Week ahead forecast

For both pipes the chosen model is the mean. For Pipe 1 that is also what ARIMA selects, and for Pipe 2 it is what the test week supports. I fit it on the full series and forecast 168 hours, one week.

p1_final <- p1_ts |> model(mean = MEAN(WaterFlow))
p2_final <- p2_ts |> model(mean = MEAN(WaterFlow))

fc1 <- p1_final |> forecast(h = 168)
fc2 <- p2_final |> forecast(h = 168)

fc1 |> autoplot(p1_ts) + labs(title = "Pipe 1, one week ahead", y = "Water flow")

fc2 |> autoplot(p2_ts) + labs(title = "Pipe 2, one week ahead", y = "Water flow")

Both forecasts are flat lines at their series means, about 19.9 for Pipe 1 and 39.8 for Pipe 2, with constant width intervals.

That is exactly what the evidence supports. A flat forecast is not a lazy one here: for a stationary series with little or no autocorrelation, the mean genuinely is the best available prediction, and the prediction interval rather than the point forecast is where the useful information lives. For these pipes the practical takeaway is the expected operating range given by the 80 percent prediction interval, roughly 14.5 to 25.3 for Pipe 1 and 25 to 55 for Pipe 2, rather than any particular hour’s number.

Part C forecast output

water_out <- bind_rows(
  fc1 |> hilo(80) |> as_tibble() |>
    transmute(Pipe = "Pipe1", DateTime = format(Hour, "%Y-%m-%d %H:%M:%S"),
              Forecast = round(.mean, 3),
              lower80 = round(`80%`$lower, 3),
              upper80 = round(`80%`$upper, 3)),
  fc2 |> hilo(80) |> as_tibble() |>
    transmute(Pipe = "Pipe2", DateTime = format(Hour, "%Y-%m-%d %H:%M:%S"),
              Forecast = round(.mean, 3),
              lower80 = round(`80%`$lower, 3),
              upper80 = round(`80%`$upper, 3))
)

write_xlsx(water_out, "Waterflow_Forecast.xlsx")

water_out |> group_by(Pipe) |> summarise(hours = n(),
                                         mean_forecast = round(mean(Forecast), 2))
## # A tibble: 2 × 3
##   Pipe  hours mean_forecast
##   <chr> <int>         <dbl>
## 1 Pipe1   168          19.9
## 2 Pipe2   168          39.8
head(water_out, 5)
## # A tibble: 5 × 5
##   Pipe  DateTime            Forecast lower80 upper80
##   <chr> <chr>                  <dbl>   <dbl>   <dbl>
## 1 Pipe1 2015-11-02 00:00:00     19.9    14.5    25.3
## 2 Pipe1 2015-11-02 01:00:00     19.9    14.5    25.3
## 3 Pipe1 2015-11-02 02:00:00     19.9    14.5    25.3
## 4 Pipe1 2015-11-02 03:00:00     19.9    14.5    25.3
## 5 Pipe1 2015-11-02 04:00:00     19.9    14.5    25.3

Waterflow_Forecast.xlsx holds 336 rows, 168 hourly forecasts for each pipe.

Conclusion

Part A. The four ATMs turned out to be four different problems rather than one problem repeated. ATM1, ATM2 and ATM4 share a strong weekly cycle with one day far below the rest (Thursday from May 2009 to January 2010, Tuesday from March 2010 onward, which the adaptive models carried into the May forecast), and are well served by seasonal models chosen on a four week test set, with ETS winning on ATM1 and STL on ATM2 and ATM4. On ATM4 the choice of STL over seasonal naive rests on RMSE alone, since seasonal naive was slightly better on MAE, so ATM4’s forecast is the least certain improvement over simply repeating last week. ATM4 also needed a single extreme outlier removed first. ATM3 was not a forecasting problem at all in the usual sense: with three days of history it required a judgement call and a deliberately wide interval rather than a model.

The month’s forecast totals about $2.06 million across the four machines. One caveat travels with it: the Ljung-Box tests on the selected models were not clean, decisively so for ATM4, so the prediction intervals in Part A should be read as a lower bound on the true uncertainty rather than a calibrated range.

Part B. Residential power has a clear upward trend and a strong twice yearly season. One missing month and one implausible July reading were interpolated before modelling. A log transformation was justified here because the seasonal swings grow with the level. ARIMA with drift was best on every accuracy measure on the two year test set, beating seasonal naive’s MAE by about 16 percent, and it left residuals indistinguishable from white noise, so its intervals can be taken at face value. The 2014 forecast totals 93,848,355 KWH, with peaks in January and August.

Part C. The two pipes recorded on different time bases, Pipe 1 irregularly at about four readings an hour and Pipe 2 exactly once an hour, and were aggregated to a common hourly mean. Both proved stationary with zero differences required, so the answer to “is it stationary” is yes for both, though Pipe 2’s KPSS p-value of 0.06 is marginal.

The answer to “can it be forecast” turned out to differ by pipe, which was the most interesting result in the project. Pipe 1 is statistically indistinguishable from white noise: the Ljung-Box p-value is 0.66 and automatic ARIMA selection returns a model with every order set to zero. It cannot meaningfully be forecast, and the mean is the best available answer. Pipe 2 does carry a weak daily cycle, significant at lag 24, which ARIMA picks up as a seasonal MA term. But on a week-long test, the ARIMA chosen automatically had roughly twice the error of the plain mean, and the mean was best on every measure. So the answer for Pipe 2 is that its structure is statistically real but not usable: it can be forecast only in the sense of giving a reliable operating range, not in the sense of predicting particular hours. The mean is the forecast for both pipes. This was the clearest case in the project of the test set overruling what the diagnostics suggested.

The recurring lesson across all three parts is that the modelling itself was the quick bit. What actually determined the quality of each forecast was the work before it: recognising that ATM3’s zeros meant a new machine rather than missing data, that ATM4’s spike and July 2010’s trough were errors rather than events, and that the water flow series were white noise rather than something to model harder.

Model decisions at a glance

Every model I fitted, whether I used it, and why.

Series Models compared Chosen Why chosen Why the others were not
ATM1 SNAIVE, ETS, ARIMA, STL+ETS ETS Lowest test RMSE, MAE and MASE STL a close second; ARIMA third; SNAIVE clearly worst (RMSE 16.9 vs 12.4)
ATM2 SNAIVE, ETS, ARIMA, STL+ETS STL+ETS Lowest test RMSE, MAE and MASE ETS second; ARIMA third; SNAIVE clearly worst (28.6 vs 18.9)
ATM4 SNAIVE, ETS, ARIMA, STL+ETS STL+ETS Lowest test RMSE (245) SNAIVE close and slightly better on MAE and MASE; ARIMA third; ETS worst on every measure
ATM3 none (3 days of data) Mean of the 3 days Only relevant data; machine started 28 April Any model fitted to 362 zeros would forecast near zero
Power SNAIVE, ETS (log), ARIMA (log), STL+ETS (log) ARIMA (log) with drift Best on all four measures; clean residuals ETS has no trend term; STL no better than SNAIVE on MAPE; SNAIVE ignores the trend
Pipe 1 Mean, Naive, ARIMA, ETS Mean (ARIMA selects it) White noise: Ljung-Box p = 0.66, ARIMA(0,0,0); tied best on the two-day test Nothing to model; Naive much worse (RMSE 6.1 vs 3.9)
Pipe 2 Mean, Naive, ARIMA, ETS Mean Best on every measure on the one-week test (RMSE 11.3) ARIMA worst of the four (RMSE 22.6) despite the weak daily cycle; ETS close second but no better; Naive assumes a random walk the data does not show

What I did and did not do

Data decisions I made, and why

  • Dropped the 14 unlabelled May 2010 ATM rows: they are placeholders for the forecast period, not data.
  • Replaced ATM4’s 9 February 2010 value (about $1.09 million in one day) with the series median: it is more than six times the next largest day and will not recur.
  • Filled the five missing ATM1 and ATM2 days with the same-weekday median, so the weekly pattern is not distorted.
  • Treated ATM3 as a new machine rather than a broken series, and forecast it from its three real days only.
  • Set the July 2010 power reading to missing and interpolated it, together with the missing September 2008, from a seasonal ARIMA, because a July cannot plausibly be the lowest month in sixteen years.
  • Used a log transformation for power but no transformation for the ATMs, for the reasons given in each part.
  • Kept the two water pipes separate rather than combining them, because they overlap for only ten days and behave differently.

Things I deliberately did not do, and why

  • Did not choose models on training accuracy. Every choice was made on a test set the models had not seen.
  • Did not use MAPE to choose ATM models. Near-zero days make it meaningless for these series.
  • Did not run time series cross-validation. Each part uses one holdout window. Cross-validation would give a steadier ranking where models are close (ETS and STL on ATM1, STL and SNAIVE on ATM4), at the cost of much longer fitting time. Where the winners beat the benchmark by a wide margin, one window is enough.
  • Did not model paydays, holidays or other calendar effects for the ATMs. The data has no labels for them, and they are the most likely cause of the leftover autocorrelation in the ATM residuals. A regression with day-of-week and holiday dummies, or ARIMA with Fourier terms, would be the next step.
  • Did not give ATM3 a weekly pattern. Three days give no evidence it shares the other machines’ pattern.
  • Did not chase the remaining residual autocorrelation in the ATM models. Instead I report that the ATM intervals are too narrow and should be read as a lower bound on uncertainty.
  • Did not use ARIMA for Pipe 2, even though the diagnostics found a daily cycle. The one-week test showed the automatically chosen ARIMA doubling the mean’s error. For a stationary series whose structure is this weak, the mean is the correct forecast, and recognising that is the finding.
  • Did not combine the two pipes into one series. They overlap for only ten days and behave differently.

Output files

tibble(
  Part = c("A", "B", "C"),
  File = c("ATM_Forecast_May2010.xlsx", "Power_Forecast_2014.xlsx",
           "Waterflow_Forecast.xlsx"),
  Rows = c(nrow(atm_out), nrow(pwr_out), nrow(water_out)),
  Contents = c("31 daily forecasts x 4 ATMs, May 2010",
               "12 monthly KWH forecasts, 2014",
               "168 hourly forecasts x 2 pipes, one week")
)
## # A tibble: 3 × 4
##   Part  File                       Rows Contents                                
##   <chr> <chr>                     <int> <chr>                                   
## 1 A     ATM_Forecast_May2010.xlsx   124 31 daily forecasts x 4 ATMs, May 2010   
## 2 B     Power_Forecast_2014.xlsx     12 12 monthly KWH forecasts, 2014          
## 3 C     Waterflow_Forecast.xlsx     336 168 hourly forecasts x 2 pipes, one week