DATA 624 - Project 1: ATM Cash, Residential Power and Water Flow Forecasts

This project forecasts cash withdrawals from four ATMs for May 2010 (Part A) and monthly residential power use for 2014 (Part B), and, as the bonus, one week of hourly water flow in two pipes (Part C). The same workflow is used for each part: look at the data, repair only what is clearly wrong and say why, fit a set of benchmark, decomposition, exponential smoothing and ARIMA models, compare them on data held out from the end of the series, check the residuals of the winner, and refit it on all the data to produce the forecast. The comparison uses a single holdout for the ATMs and the water flow, and a holdout plus cross-validation for the power series; the reasons are given in each part.

Part A - ATM cash forecast

The file holds daily cash withdrawals in hundreds of dollars for four ATMs. The goal is a daily forecast for every ATM for all 31 days of May 2010, one ATM at a time, because a bank manager has to stock each machine separately.

A1. Reading the data

atm_raw <- read_excel("ATM624Data.xlsx") |>
  mutate(DATE = as.Date(DATE, origin = "1899-12-30"))   # Excel stores dates as day numbers

# the file ends with 14 rows (1-14 May 2010) that have a date but no ATM and no Cash
atm_raw |> filter(is.na(ATM)) |> summarise(rows = n(), first = min(DATE), last = max(DATE))

## # A tibble: 1 × 3
##    rows first      last      
##   <int> <date>     <date>    
## 1    14 2010-05-01 2010-05-14

atm <- atm_raw |> filter(!is.na(ATM)) |> as_tsibble(key = ATM, index = DATE)
# for each ATM: rows in the data against the number of days from the first to the last date
atm |>
  as_tibble() |>
  group_by(ATM) |>
  summarise(first = min(DATE), last = max(DATE), rows = n(),
            days_first_to_last = as.numeric(last - first) + 1) |>
  knitr::kable()

ATM	first	last	rows	days_first_to_last
ATM1	2009-05-01	2010-04-30	365	365
ATM2	2009-05-01	2010-04-30	365	365
ATM3	2009-05-01	2010-04-30	365	365
ATM4	2009-05-01	2010-04-30	365	365

The 14 blank rows at the end are placeholders and not observations, so they are dropped. After that every ATM has all 365 days from 1 May 2009 to 30 April 2010 and the number of rows equals the number of days from the first to the last date, so no date is missing and the daily sequence is intact. Keeping the complete daily index preserves the spacing of the series, so a missing value is filled in and not dropped.

A2. Exploratory data analysis

atm |>
  as_tibble() |>   # summary table of the raw data, not a forecasting step
  group_by(ATM) |>
  summarise(days = n(), missing = sum(is.na(Cash)), zeros = sum(Cash == 0, na.rm = TRUE),
            min = min(Cash, na.rm = TRUE), median = median(Cash, na.rm = TRUE),
            mean = mean(Cash, na.rm = TRUE), sd = sd(Cash, na.rm = TRUE),
            max = max(Cash, na.rm = TRUE)) |>
  knitr::kable(digits = 1)

ATM	days	missing	zeros	min	median	mean	sd	max
ATM1	365	3	0	1.0	91.0	83.9	36.7	180.0
ATM2	365	2	2	0.0	67.0	62.6	38.9	147.0
ATM3	365	0	362	0.0	0.0	0.7	7.9	96.0
ATM4	365	0	0	1.6	403.8	474.0	650.9	10919.8

# the only days on which ATM3 has a withdrawal
atm |> filter(ATM == "ATM3", Cash > 0) |> as_tibble() |> knitr::kable()

DATE	ATM	Cash
2010-04-28	ATM3	96
2010-04-29	ATM3	82
2010-04-30	ATM3	85

# the three largest ATM4 values, to show how far the outlier sits from the rest
atm |> filter(ATM == "ATM4") |> as_tibble() |> slice_max(Cash, n = 3) |> knitr::kable(digits = 1)

DATE	ATM	Cash
2010-02-09	ATM4	10919.8
2009-09-22	ATM4	1712.1
2010-01-29	ATM4	1574.8

atm |> autoplot(Cash) +
  facet_grid(vars(ATM), scales = "free_y") +
  labs(title = "Daily cash withdrawals by ATM", y = "Cash (hundreds of dollars)", x = "Date") +
  theme(legend.position = "none")

Three things stand out.

• ATM3 has a withdrawal of zero on 362 of its 365 days. The only non-zero values are the three days in the ATM3 table above. There is no history to model for this machine; it looks like it only started operating at the very end of the data. That is an assumption, not something the data proves.
• ATM4 has one value of 10,920 on 9 February 2010, while the next largest is 1,712 (largest-values table above). In dollars that is roughly 1.1 million in a day for one ATM, far outside anything else the machine did. It is most likely an error and would distort every model fitted to the series.
• ATM1 and ATM2 have a few missing values (3 and 2 in the summary table), and ATM2 has two days with zero withdrawals.

Level of each machine

The three working ATMs are on very different scales: the median withdrawal is 403.8 for ATM4, 91.0 for ATM1 and 67.0 for ATM2 (the ATM4 standard deviation is inflated by the outlier, so the median is the better guide). That is why the ATMs are modeled one at a time.

A3. Checking the weekly pattern

This section asks whether the weekly pattern of each working ATM is the same all year. It goes in five steps: spot a possible change, locate it, compare the weekday profiles, look at how deep the low day is, and describe an adjustment I considered and did not use. The recorded dates are not changed.

Step 1: spotting it. The heatmap has one cell per day, with weeks across and weekdays down. Dark blue cells are days with almost no withdrawals compared with that machine’s median.

# heatmap of cash by weekday and week; dark cells are days with almost no withdrawals (display only)
weekday_heatmap <- function(data, title, subtitle) {
  data |>
    filter(ATM != "ATM3") |>
    mutate(Cash = if_else(ATM == "ATM4" & Cash > 5000, NA_real_, Cash),   # leave out the ATM4 outlier
           weekday = factor(as.character(wday(DATE, label = TRUE, week_start = 1)),
                            levels = rev(c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"))),
           week_start = floor_date(DATE, "week", week_start = 1)) |>
    group_by(ATM) |>
    mutate(relative = pmin(Cash / median(Cash, na.rm = TRUE), 2)) |>   # base R pmin() caps the ratio at 2 so the color scale stays readable (display only)
    ungroup() |>
    ggplot(aes(week_start, weekday, fill = relative)) +
    geom_tile(colour = "white", linewidth = 0.2) +
    facet_wrap(vars(ATM), ncol = 1) +
    scale_fill_gradient2(low = "navy", mid = "whitesmoke", high = "darkorange", midpoint = 1,
                         na.value = "gray70", name = "Cash / median") +
    labs(title = title, subtitle = subtitle, x = "Week starting", y = "Weekday")
}

atm |> as_tibble() |>
  weekday_heatmap("Cash by weekday and week, as recorded",
                  "Dark = almost no withdrawals. Gray = missing value or the ATM4 outlier.")

In every machine the dark row is Thursday until mid-February and Tuesday afterwards. ATM2 has two dark rows, Wednesday and Thursday before the change and Monday and Tuesday after it, so the change involves more than the lowest day.

Step 2: locating it. For each week from January to April 2010 (weeks start on Monday), the table gives the weekday and date with the lowest withdrawal.

# weekday of the lowest withdrawal in each week, leaving out the missing values and the ATM4 outlier
atm |>
  filter(ATM != "ATM3", DATE >= as.Date("2010-01-04"), !is.na(Cash), !(ATM == "ATM4" & Cash > 5000)) |>
  as_tibble() |>
  mutate(week_start = floor_date(DATE, "week", week_start = 1)) |>
  group_by(ATM, week_start) |>
  slice_min(Cash, n = 1, with_ties = FALSE) |>
  ungroup() |>
  mutate(lowest = paste(wday(DATE, label = TRUE), format(DATE, "%d %b"))) |>
  select(week_start, ATM, lowest) |>
  pivot_wider(names_from = ATM, values_from = lowest) |>
  knitr::kable()

week_start	ATM1	ATM2	ATM4
2010-01-04	Thu 07 Jan	Thu 07 Jan	Sat 09 Jan
2010-01-11	Thu 14 Jan	Thu 14 Jan	Thu 14 Jan
2010-01-18	Thu 21 Jan	Thu 21 Jan	Wed 20 Jan
2010-01-25	Thu 28 Jan	Thu 28 Jan	Thu 28 Jan
2010-02-01	Thu 04 Feb	Thu 04 Feb	Sun 07 Feb
2010-02-08	Sun 14 Feb	Thu 11 Feb	Wed 10 Feb
2010-02-15	Mon 15 Feb	Mon 15 Feb	Sat 20 Feb
2010-02-22	Tue 23 Feb	Tue 23 Feb	Tue 23 Feb
2010-03-01	Tue 02 Mar	Tue 02 Mar	Sun 07 Mar
2010-03-08	Tue 09 Mar	Tue 09 Mar	Tue 09 Mar
2010-03-15	Tue 16 Mar	Tue 16 Mar	Tue 16 Mar
2010-03-22	Tue 23 Mar	Tue 23 Mar	Tue 23 Mar
2010-03-29	Tue 30 Mar	Tue 30 Mar	Tue 30 Mar
2010-04-05	Tue 06 Apr	Tue 06 Apr	Tue 06 Apr
2010-04-12	Tue 13 Apr	Tue 13 Apr	Tue 13 Apr
2010-04-19	Tue 20 Apr	Tue 20 Apr	Sat 24 Apr
2010-04-26	Tue 27 Apr	Tue 27 Apr	Tue 27 Apr

Down the ATM1 and ATM2 columns the lowest day is a Thursday in every week through the week of 1 February (ATM2 is still on Thursday 11 February in the week of 8 February). The weeks of 8 and 15 February are mixed (Sunday, Monday, Thursday), and from the week of 22 February the lowest day is a Tuesday in every week through the end of April. ATM4 makes the same move from Thursday to Tuesday, but its lowest day is less regular because the machine is noisier: in the five weeks before the move it is a Thursday in two (the others are Saturday, Wednesday and Sunday), and in the ten weeks after it a Tuesday in eight, with a Sunday (week of 1 March) and a Saturday (week of 19 April) in the other two.

Step 3: comparing the weekday profiles. The next graph gives the average cash for each weekday before the change (1 May 2009 to 8 February 2010) and after it (22 February to 30 March 2010, the last day used for fitting the models). The weeks in between are left out.

# average cash by weekday and period, leaving out the missing values and the ATM4 outlier (display only)
weekday_profile <- function(data) {
  data |>
    filter(ATM != "ATM3", !is.na(Cash), !(ATM == "ATM4" & Cash > 5000)) |>
    mutate(period = case_when(DATE < as.Date("2010-02-09") ~ "before",
                              DATE >= as.Date("2010-02-22") & DATE <= as.Date("2010-03-30") ~ "after"),
           weekday = wday(DATE, label = TRUE, week_start = 1)) |>
    filter(!is.na(period)) |>
    group_by(ATM, period, weekday) |>
    summarise(mean_cash = mean(Cash), .groups = "drop")
}

atm |> as_tibble() |> weekday_profile() |>
  mutate(period = if_else(period == "before", "before (1 May 2009 to 8 Feb 2010)",
                          "after (22 Feb to 30 Mar 2010)")) |>
  ggplot(aes(weekday, mean_cash, colour = period, group = period)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  facet_wrap(vars(ATM), scales = "free_y", ncol = 1) +
  scale_colour_manual(values = c("firebrick", "steelblue")) +
  labs(title = "Average cash by weekday, before and after the change",
       x = "Weekday", y = "Average cash (hundreds of dollars)", colour = NULL) +
  theme(legend.position = "bottom")

In every working machine the lowest weekday average is Thursday before the change and Tuesday after it. ATM2 also has a second low day, Wednesday before and Monday after. The other weekdays change too, and the two weekday profiles resemble a two-day shift.

Step 4: how deep is the low day? The table gives the withdrawals on the low day of each working machine: the Thursdays before 9 February and the Tuesdays from 22 February to 30 March.

# withdrawals on the low day before and after the change (display only)
atm |>
  filter(ATM != "ATM3", !is.na(Cash), !(ATM == "ATM4" & Cash > 5000)) |>
  as_tibble() |>
  mutate(weekday = wday(DATE, label = TRUE),
         low_day = case_when(DATE < as.Date("2010-02-09") & weekday == "Thu" ~ "Thursdays before 9 Feb",
                             DATE >= as.Date("2010-02-22") & DATE <= as.Date("2010-03-30") & weekday == "Tue" ~
                               "Tuesdays 22 Feb to 30 Mar")) |>
  filter(!is.na(low_day)) |>
  group_by(ATM, low_day) |>
  summarise(days = n(), mean = mean(Cash), median = median(Cash), sd = sd(Cash), .groups = "drop") |>
  knitr::kable(digits = 1)

ATM	low_day	days	mean	median	sd
ATM1	Thursdays before 9 Feb	40	15.3	14.0	7.0
ATM1	Tuesdays 22 Feb to 30 Mar	6	5.8	4.0	6.6
ATM2	Thursdays before 9 Feb	39	8.1	7.0	6.3
ATM2	Tuesdays 22 Feb to 30 Mar	6	1.5	1.5	1.0
ATM4	Thursdays before 9 Feb	40	62.3	48.6	54.6
ATM4	Tuesdays 22 Feb to 30 Mar	6	22.2	24.9	10.2

The low day is deeper after the change as well as on a different weekday. The median withdrawal on the low day falls from 14.0 to 4.0 for ATM1, from 7.0 to 1.5 for ATM2 and from 48.6 to 24.9 for ATM4. For ATM2 and ATM4 the low-day withdrawals are also less spread out after the change (standard deviation 1.0 against 6.3 and 10.2 against 54.6); for ATM1 the spread is about the same (6.6 against 7.0).

Step 5: an adjustment I considered and did not use. One option is to move the earlier observations two days earlier, so that their weekly pattern matches the later one. The graph shows what that would do: it overlays the average weekday profile before the change, moved two days earlier, on the profile after the change.

# what moving the earlier days two days earlier would do to the weekday profile (display only)
weekday_names <- c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")
profile <- atm |> as_tibble() |> weekday_profile()

bind_rows(
  profile |> filter(period == "after") |>
    mutate(weekday = factor(as.character(weekday), levels = weekday_names),
           period = "after (22 Feb to 30 Mar 2010)"),
  profile |> filter(period == "before") |>
    mutate(weekday = factor(weekday_names[(as.integer(weekday) - 1 - 2) %% 7 + 1], levels = weekday_names),
           period = "before, moved two days earlier")) |>
  ggplot(aes(weekday, mean_cash, colour = period, group = period)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  facet_wrap(vars(ATM), scales = "free_y", ncol = 1) +
  scale_colour_manual(values = c("firebrick", "steelblue")) +
  labs(title = "The earlier weekday profile moved two days earlier, against the later profile",
       x = "Weekday", y = "Average cash (hundreds of dollars)", colour = NULL) +
  theme(legend.position = "bottom")

The move brings the low days into the same weekday positions, but it would leave two dates empty at the transition, and those would then have to be filled. The table in step 4 also shows that withdrawals on the low day became lower after the change, so the difference involves both the timing and the level of the withdrawals. The data cannot show whether the change came from wrong dates or from a change in how the machines are operated. Moving the dates would impose an unverified explanation and create dates that are not in the original records, so I kept the recorded dates and treat the February change as a limitation: models fitted to the whole year are compared with models fitted only to the days after the change (section A7), and the May forecast assumes that the later pattern continues.

A4. Missing values and the outlier

The usual rule for a time series is not to delete rows: a bad or missing value is replaced by the mean or median of comparable days, for example the same weekday, and the replacement is documented. Kuhn and Johnson (Applied Predictive Modeling, sections 3.3 and 3.4) add two rules that apply here: an unusual value should first be checked for validity and recording errors, and should not be removed or changed hastily, and an imputation should be built from the training data only. Here every value that needs replacing is replaced by the median of the same weekday for that machine over the days before 9 February 2010, the days before the weekly pattern changed; this uses no day of the 31-day holdout. A median is used because it is not pulled by extreme values. The table lists every replaced value.

change_start <- as.Date("2010-02-09")   # first day of the weeks in which the weekly pattern changes

atm_flagged <- atm |>
  as_tibble() |>
  mutate(replace = is.na(Cash) | (ATM == "ATM4" & Cash > 5000),   # missing, or the 10,920 value
         weekday = wday(DATE, label = TRUE))

# median Cash for each ATM and weekday, from the days before the change that are not replaced
weekday_median <- atm_flagged |>
  filter(!replace, DATE < change_start) |>
  group_by(ATM, weekday) |>
  summarise(typical = median(Cash), .groups = "drop")

atm_clean <- atm_flagged |>
  left_join(weekday_median, by = c("ATM", "weekday")) |>
  mutate(Cash_original = Cash,
         Cash = if_else(replace, typical, Cash)) |>   # swap in the weekday median only where flagged
  as_tsibble(key = ATM, index = DATE)

atm_clean |>
  filter(replace) |>
  as_tibble() |>
  select(ATM, DATE, weekday, original = Cash_original, replacement = Cash) |>
  knitr::kable(digits = 1)

ATM	DATE	weekday	original	replacement
ATM1	2009-06-13	Sat	NA	96.5
ATM1	2009-06-16	Tue	NA	108.0
ATM1	2009-06-22	Mon	NA	90.0
ATM2	2009-06-18	Thu	NA	7.0
ATM2	2009-06-24	Wed	NA	26.0
ATM4	2010-02-09	Tue	10919.8	483.4

The table shows six values: the five values missing in the file (June 2009) and the ATM4 outlier on 9 February 2010. The ATM4 mean in the A2 table is above its median, partly because of that outlier, which is one reason to prefer the median. The zeros of ATM2 are not replaced; the next table shows where they fall and what is typical for that weekday in the same period.

# the ATM2 zeros, with the typical value for their weekday before or after the change
atm_clean |>
  as_tibble() |>
  filter(ATM == "ATM2") |>
  mutate(period = case_when(DATE < change_start ~ "before", DATE >= as.Date("2010-02-22") ~ "after")) |>
  group_by(period, weekday) |>
  mutate(typical_same_period = median(Cash[!replace])) |>
  ungroup() |>
  filter(Cash_original == 0) |>
  select(DATE, weekday, period, original = Cash_original, typical_same_period) |>
  knitr::kable(digits = 1)

DATE	weekday	period	original	typical_same_period
2009-10-25	Sun	before	0	62.0
2010-03-30	Tue	after	0	1.5

The zero on 30 March 2010 is a Tuesday, the low day after the change, where the typical ATM2 value is 1.5, so it is ordinary. The zero on 25 October 2009 is a Sunday, where the typical value before the change is 62.0, so it is unusual. A zero withdrawal can happen and nothing in the data shows it is an error, so both zeros stay in.

A5. Transformations and the training and test split

train_end <- as.Date("2010-03-30")   # last day used for fitting; the last 31 days are kept for testing
atm_train <- atm_clean |> filter(ATM != "ATM3", DATE <= train_end)   # the test set is 31 Mar to 30 Apr (31 days)

# Guerrero's lambda for each ATM, from training data only
atm_train |> features(Cash, guerrero) |> knitr::kable(digits = 3)

ATM	lambda_guerrero
ATM1	0.231
ATM2	0.524
ATM4	0.308

The Guerrero lambdas are 0.231 for ATM1, 0.524 for ATM2 and 0.308 for ATM4; each machine uses its own value. As Kuhn and Johnson do for the Box-Cox transformation (section 3.2), the lambdas are estimated from the training days only. Box-Cox versions of the ETS and ARIMA models are included in the comparison to see whether they help. The last 31 days (31 March to 30 April 2010) are held out because May is also 31 days. The training period includes the days from 22 February to 30 March, which follow the later weekly pattern, so the holdout tests whether a model has picked that pattern up.

A6. Choosing candidate ARIMA models

The first table gives the strength of the trend and of the weekly pattern for each machine (between 0 and 1), and the graph gives the autocorrelation of the cleaned series.

# strength of the trend and of the weekly pattern (between 0 and 1)
atm_train |> features(Cash, feat_stl) |>
  select(ATM, trend_strength, seasonal_strength_week) |>
  knitr::kable(digits = 2)

ATM	trend_strength	seasonal_strength_week
ATM1	0.34	0.77
ATM2	0.30	0.78
ATM4	0.27	0.35

# autocorrelation of the cleaned series: spikes at lags 7, 14, 21 mean a weekly pattern
atm_train |> ACF(Cash, lag_max = 28) |> autoplot() +
  labs(title = "ACF of daily cash by ATM")

The weekly pattern is strong for ATM1 and ATM2 and weaker for ATM4, whose withdrawals are noisier. The ACF of ATM1 and ATM2 has large spikes at lags 7, 14, 21 and 28, which confirms a weekly pattern, so the seasonal period is 7; ATM4 has a much smaller spike at lag 7, in line with its lower seasonal strength.

The next table gives the KPSS statistic and the number of seasonal and ordinary differences suggested for each machine.

atm_train |> features(Cash, list(unitroot_kpss, unitroot_nsdiffs, unitroot_ndiffs)) |>
  knitr::kable(digits = 3)

ATM	kpss_stat	kpss_pvalue	nsdiffs	ndiffs
ATM1	0.376	0.087	1	0
ATM2	2.242	0.010	1	1
ATM4	0.087	0.100	0	0

ATM1 and ATM2 need one seasonal difference, and ATM4 none.

The graph shows ATM1 after one seasonal difference, with its ACF and PACF.

atm_train |> filter(ATM == "ATM1") |>
  gg_tsdisplay(difference(Cash, 7), plot_type = "partial") +
  labs(title = "ATM1 after one seasonal difference (lag 7)")

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

After the seasonal difference the ACF has a spike at lag 1 and a spike at lag 7, and the PACF decays at the seasonal lags (the warnings about removed rows are the first 7 values lost to the seasonal difference). That suggests ARIMA(0,0,1)(0,1,1)[7] as a manual candidate, and ARIMA(1,0,0)(0,1,1)[7] as a close alternative. These are fitted next to the automatic ARIMA() and the full search (stepwise = FALSE, approx = FALSE), so that several ARIMA models are compared and not one automatic fit.

What the exploration led to

What the exploration shows	Where it is shown	What it led to
ATM3 is zero except for the last three days	A2 summary and ATM3 tables	ATM3 gets a simple forecast and no model comparison
One suspected erroneous ATM4 value; a few missing values in ATM1 and ATM2	A2 tables and time plot, A4	Replace them with same-weekday medians, keep the row order
A weekly pattern in every working ATM (strong for ATM1 and ATM2, weaker for ATM4)	A3 heatmap, A6 strength table and ACF	Seasonal period 7 in every seasonal model
The weekly pattern moves by two days between 9 and 22 February 2010, and the low day becomes deeper	A3 steps 1 to 4	Keep the recorded dates, take the replacement medians from before the change, compare models fitted to the whole year with models fitted to the later days only (A7), and assume the later pattern continues in May

A7. Fitting and comparing the models

Each ATM gets the same set of models fitted to the days of the whole year up to 30 March: the four benchmarks (mean, naive, seasonal naive, drift), automatic ETS, ETS on Box-Cox data, STL with ETS on the seasonally adjusted data, automatic ARIMA, the full ARIMA search, two manual seasonal ARIMA models and a Box-Cox ARIMA. Because the weekly pattern changed in February (section A3), three more models are fitted only to the days from 22 February: seasonal naive, ETS and automatic ARIMA. They show whether the older days hurt, at the cost of a very short training series; the first table gives the number of training days of each history.

recent_start <- as.Date("2010-02-22")   # first week with the later weekly pattern

# model definitions for the whole-year history; lambda is the Box-Cox parameter of the series
atm_specs <- function(lambda) {
  list(
    mean = MEAN(Cash),
    naive = NAIVE(Cash),
    snaive = SNAIVE(Cash),
    drift = RW(Cash ~ drift()),
    ets = ETS(Cash),
    ets_boxcox = ETS(box_cox(Cash, lambda)),   # Box-Cox inside ETS() is only mentioned in exercise 8.13 of the book
    stl_ets = decomposition_model(STL(Cash ~ season(window = 7), robust = TRUE),
                                  ETS(season_adjust ~ season("N"))),
    arima_auto = ARIMA(Cash),
    arima_search = ARIMA(Cash, stepwise = FALSE, approx = FALSE),
    arima_001_011 = ARIMA(Cash ~ 0 + pdq(0, 0, 1) + PDQ(0, 1, 1)),
    arima_100_011 = ARIMA(Cash ~ 0 + pdq(1, 0, 0) + PDQ(0, 1, 1)),
    arima_boxcox = ARIMA(box_cox(Cash, lambda))
  )
}

# fit all the whole-year models, or only those named in `which`
fit_atm <- function(series, which = NULL) {
  lambda <- series |> features(Cash, guerrero) |> pull(lambda_guerrero)
  specs <- atm_specs(lambda)
  if (!is.null(which)) specs <- specs[which]
  do.call(model, c(list(series), specs))
}

# the three models for the days from 22 February only
fit_atm_recent <- function(series, which = NULL) {
  specs <- list(snaive_recent = SNAIVE(Cash), ets_recent = ETS(Cash), arima_auto_recent = ARIMA(Cash))
  if (!is.null(which)) specs <- specs[which]
  do.call(model, c(list(series), specs))
}

working_atms <- c("ATM1", "ATM2", "ATM4")
atm_fits <- map(set_names(working_atms), \(k) fit_atm(atm_train |> filter(ATM == k)))
atm_fits_recent <- map(set_names(working_atms),
                       \(k) fit_atm_recent(atm_train |> filter(ATM == k, DATE >= recent_start)))

# training days used by each history
atm_train |>
  as_tibble() |>
  group_by(ATM) |>
  summarise(whole_year = n(), from_22_Feb = sum(DATE >= recent_start)) |>
  knitr::kable()

ATM	whole_year	from_22_Feb
ATM1	334	37
ATM2	334	37
ATM4	334	37

A helper computes the Ljung-Box test with the book’s rules: lag 2m = 14 for a weekly pattern and degrees of freedom p + q + P + Q for ARIMA models. Candidate models are then compared by RMSE and MAE on the holdout, together with the residual test.

# Ljung-Box p-value on innovation residuals, the same test as ljung_box(lag, dof) in the book;
# dof differs by model (number of ARMA terms for ARIMA models, 0 otherwise), so it is applied here with pchisq()
ljung_table <- function(fit, lag) {
  stats <- fit |> augment() |> features(.innov, ljung_box, lag = lag) |>
    select(.model, lb_stat)
  arima_terms <- fit |> select(starts_with("arima")) |> tidy() |>
    filter(term != "constant") |> count(.model, name = "dof")
  stats |> left_join(arima_terms, by = ".model") |>
    mutate(dof = replace_na(dof, 0),
           lb_pvalue = pchisq(lb_stat, df = pmax(lag - dof, 1), lower.tail = FALSE)) |>
    select(.model, lb_pvalue)
}

# holdout accuracy and residual test of one set of fits
score_fits <- function(fits, test, horizon, lag, history) {
  accuracy(forecast(fits, h = horizon), test) |>
    as_tibble() |> select(.model, RMSE, MAE) |>
    left_join(ljung_table(fits, lag), by = ".model") |>
    mutate(history = history)
}

atm_accuracy <- map(working_atms, \(k) {
  test <- atm_clean |> filter(ATM == k)
  bind_rows(score_fits(atm_fits[[k]], test, 31, 14, "whole year"),
            score_fits(atm_fits_recent[[k]], test, 31, 14, "from 22 Feb")) |>
    mutate(ATM = k)
}) |> list_rbind()

atm_accuracy |>
  arrange(ATM, RMSE) |>
  select(ATM, model = .model, history, RMSE, MAE, ljung_box_p = lb_pvalue) |>
  knitr::kable(digits = 3)

ATM	model	history	RMSE	MAE	ljung_box_p
ATM1	ets	whole year	12.748	9.878	0.012
ATM1	arima_001_011	whole year	12.848	10.011	0.088
ATM1	arima_100_011	whole year	12.861	10.013	0.051
ATM1	arima_auto_recent	from 22 Feb	13.195	10.574	0.126
ATM1	stl_ets	whole year	13.788	10.388	0.014
ATM1	arima_auto	whole year	13.846	10.782	0.116
ATM1	arima_boxcox	whole year	14.340	11.675	0.544
ATM1	ets_boxcox	whole year	14.921	11.960	0.180
ATM1	arima_search	whole year	15.584	11.815	0.158
ATM1	snaive	whole year	15.862	13.161	0.000
ATM1	snaive_recent	from 22 Feb	15.862	13.161	0.031
ATM1	ets_recent	from 22 Feb	16.021	12.122	0.217
ATM1	mean	whole year	32.237	20.174	0.000
ATM1	naive	whole year	67.198	63.290	0.000
ATM1	drift	whole year	70.341	65.946	0.000
ATM2	ets_boxcox	whole year	19.292	14.121	0.008
ATM2	ets	whole year	19.414	14.695	0.001
ATM2	arima_100_011	whole year	19.469	14.696	0.001
ATM2	arima_001_011	whole year	19.471	14.698	0.001
ATM2	arima_auto_recent	from 22 Feb	20.224	14.393	0.287
ATM2	ets_recent	from 22 Feb	20.485	13.708	0.664
ATM2	stl_ets	whole year	20.485	13.565	0.819
ATM2	arima_boxcox	whole year	20.747	15.711	0.291
ATM2	arima_auto	whole year	20.848	17.006	0.485
ATM2	arima_search	whole year	23.209	18.668	0.421
ATM2	snaive	whole year	26.821	19.839	0.000
ATM2	snaive_recent	from 22 Feb	26.821	19.839	0.057
ATM2	mean	whole year	37.194	32.386	0.000
ATM2	naive	whole year	72.909	62.710	0.000
ATM2	drift	whole year	77.406	67.851	0.000
ATM4	ets_boxcox	whole year	245.985	203.635	0.530
ATM4	stl_ets	whole year	251.690	207.075	0.051
ATM4	arima_boxcox	whole year	252.852	216.871	0.212
ATM4	ets_recent	from 22 Feb	259.887	207.516	0.577
ATM4	arima_100_011	whole year	268.278	226.576	0.408
ATM4	arima_001_011	whole year	268.508	226.698	0.401
ATM4	arima_search	whole year	268.727	226.967	0.326
ATM4	arima_auto	whole year	280.540	237.567	0.413
ATM4	arima_auto_recent	from 22 Feb	281.337	235.661	0.687
ATM4	mean	whole year	281.377	235.697	0.024
ATM4	ets	whole year	304.354	255.172	0.414
ATM4	snaive	whole year	329.584	260.177	0.000
ATM4	snaive_recent	from 22 Feb	329.584	260.177	0.215
ATM4	naive	whole year	466.056	378.190	0.000
ATM4	drift	whole year	493.200	411.604	0.000

RMSE and MAE are in hundreds of dollars. The naive and drift forecasts, which have no weekly pattern, are far behind for every ATM, and the plain mean is far behind for ATM1 and ATM2 but not for ATM4. The seasonal naive forecast is the right simple benchmark and is beaten by the seasonal ETS and ARIMA models; the seasonal naive forecast fitted from 22 February is identical to the whole-year one, because it only repeats the last week. RMSE gives more weight to large misses than MAE does; for the top models of each ATM the two measures usually give a similar order, with some differences (for ATM2 the STL model has the lowest MAE but not the lowest RMSE). The three models fitted only to the days from 22 February are not clearly better than the whole-year models. For ATM1 and ATM4 the best of them (RMSE 13.195 and 259.887) is behind the best whole-year model (12.748 and 245.985). For ATM2 the ARIMA fitted from 22 February (20.224) comes first among the models that pass the residual test, because the whole-year ETS and ARIMA models with a lower RMSE leave autocorrelation in the residuals (Ljung-Box p-values of 0.008 or less).

Selection rule. For each ATM I take the model with the lowest holdout RMSE among those whose residuals pass the Ljung-Box test at the 5% level, whichever history it was fitted to. Accuracy decides, and the residual test checks whether the residuals still show autocorrelation; passing it does not prove that there is none.

atm_selected <- atm_accuracy |>
  filter(lb_pvalue >= 0.05) |>   # keep models that pass the residual test
  group_by(ATM) |>
  slice_min(RMSE, n = 1) |>   # lowest holdout RMSE for each ATM
  ungroup()

atm_selected |> select(ATM, model = .model, history, RMSE, MAE, ljung_box_p = lb_pvalue) |>
  knitr::kable(digits = 3)

ATM	model	history	RMSE	MAE	ljung_box_p
ATM1	arima_001_011	whole year	12.848	10.011	0.088
ATM2	arima_auto_recent	from 22 Feb	20.224	14.393	0.287
ATM4	ets_boxcox	whole year	245.985	203.635	0.530

Kuhn and Johnson (section 4.6) advise favoring a simpler model whose performance is within a small tolerance of the best. The table lists, for each ATM, the models that pass the residual test and are within 1% of the lowest RMSE.

# passing models within 1% of the best RMSE, for each ATM
atm_accuracy |>
  filter(lb_pvalue >= 0.05) |>
  group_by(ATM) |>
  mutate(percent_above_best = 100 * (RMSE / min(RMSE) - 1)) |>
  filter(percent_above_best <= 1) |>
  ungroup() |>
  arrange(ATM, RMSE) |>
  select(ATM, model = .model, history, RMSE, percent_above_best) |>
  knitr::kable(digits = 3)

ATM	model	history	RMSE	percent_above_best
ATM1	arima_001_011	whole year	12.848	0.000
ATM1	arima_100_011	whole year	12.861	0.102
ATM2	arima_auto_recent	from 22 Feb	20.224	0.000
ATM4	ets_boxcox	whole year	245.985	0.000

For ATM2 and ATM4 the best model is alone in that range. For ATM1 the range holds two seasonal ARIMA models that differ only in a non-seasonal MA(1) against AR(1) term and by 0.1% in RMSE, and neither is clearly simpler, so I keep the lowest-RMSE choice. The rule therefore changes no choice for the ATMs.

# residual check for each selected model
for (k in atm_selected$ATM) {
  sel <- atm_selected |> filter(ATM == k)
  fits <- if (sel$history == "whole year") atm_fits[[k]] else atm_fits_recent[[k]]
  print(fits |> select(all_of(sel$.model)) |> gg_tsresiduals() +
          labs(title = paste(k, "-", sel$.model)))
}

The rule picks a different kind of model for each ATM. For ATM1 the automatic ETS has the lowest RMSE in the accuracy table (12.748), but its residuals fail the Ljung-Box test (p = 0.012), so the rule picks the ARIMA(0,0,1)(0,1,1)[7] model (12.848, p = 0.088). For ATM2 it picks the automatic ARIMA fitted from 22 February (20.224, p = 0.287). For ATM4 it picks the ETS on Box-Cox data (245.985, p = 0.530), which is not far below the plain mean (281.377). The residual plots support the choices. The autocorrelations of ATM1 stay inside the dashed lines except for one lag that is just outside, and those of ATM4 stay inside. The ATM1 residuals are bunched around zero with a few large values, mostly negative (days with much less cash than the model expects), so their tails are heavy and the prediction intervals are approximate; the ATM4 residuals have a few large positive values. The ATM2 model rests on only 37 training days, so the first seven residuals are zero and the autocorrelation bounds are wide; its choice is less certain than the others.

Check on the Tuesdays

In section A3 the low day is much lower after the change than before it, so a model fitted to the whole year may not have caught up with that. The next tables compare the forecast with the actual values for the four Tuesdays of the holdout (6, 13, 20 and 27 April), using the selected model of each ATM.

# forecast and actual for every holdout day, using the selected model of each ATM
holdout_check <- map(working_atms, \(k) {
  sel <- atm_selected |> filter(ATM == k)
  fits <- if (sel$history == "whole year") atm_fits[[k]] else atm_fits_recent[[k]]
  fits |> select(all_of(sel$.model)) |> forecast(h = 31) |>
    as_tibble() |>
    select(DATE, forecast = .mean) |>
    left_join(atm_clean |> filter(ATM == k) |> as_tibble() |> select(DATE, actual = Cash), by = "DATE") |>
    mutate(ATM = k, model = sel$.model, weekday = wday(DATE, label = TRUE))
}) |> list_rbind()

holdout_check |>
  filter(weekday == "Tue") |>
  transmute(ATM, model, DATE, actual, forecast, forecast_minus_actual = forecast - actual) |>
  knitr::kable(digits = 1)

ATM	model	DATE	actual	forecast	forecast_minus_actual
ATM1	arima_001_011	2010-04-06	4.0	16.4	12.4
ATM1	arima_001_011	2010-04-13	3.0	16.4	13.4
ATM1	arima_001_011	2010-04-20	3.0	16.4	13.4
ATM1	arima_001_011	2010-04-27	4.0	16.4	12.4
ATM2	arima_auto_recent	2010-04-06	2.0	0.6	-1.4
ATM2	arima_auto_recent	2010-04-13	1.0	0.2	-0.8
ATM2	arima_auto_recent	2010-04-20	1.0	0.5	-0.5
ATM2	arima_auto_recent	2010-04-27	2.0	0.3	-1.7
ATM4	ets_boxcox	2010-04-06	2.5	142.1	139.6
ATM4	ets_boxcox	2010-04-13	27.3	144.9	117.6
ATM4	ets_boxcox	2010-04-20	29.9	147.8	117.9
ATM4	ets_boxcox	2010-04-27	13.7	150.6	136.9

# how much of each machine's squared holdout error comes from the four Tuesdays
holdout_check |>
  group_by(ATM) |>
  summarise(mean_error_on_Tuesdays = mean((forecast - actual)[weekday == "Tue"]),
            percent_of_squared_error_from_Tuesdays =
              100 * sum(((forecast - actual)^2)[weekday == "Tue"]) / sum((forecast - actual)^2)) |>
  knitr::kable(digits = 2)

ATM	mean_error_on_Tuesdays	percent_of_squared_error_from_Tuesdays
ATM1	12.88	12.99
ATM2	-1.09	0.04
ATM4	128.01	3.52

For ATM2 the forecast is close to the actual Tuesdays (0.2 to 0.6 against 1 to 2). For ATM1 and ATM4 the selected models forecast Tuesdays that are far too high: 16.4 against actual values of 3 to 4 for ATM1, and 142.1 to 150.6 against 2.5 to 29.9 for ATM4. On average that is 12.9 and 128.0 hundred dollars too high. A likely reason is that these two models are fitted to the whole year and at the end of March had seen only six of the later Tuesdays (the table in section A3), so they had not yet picked up the deeper low day. I keep the selection rule unchanged and record this as a limitation. For the selected models, Tuesdays are 13.0% of ATM1’s squared holdout error, 0.04% of ATM2’s and 3.5% of ATM4’s, so the ranking of the models is also driven by the other days.

A8. ATM3

ATM3 has only three active days (96, 82 and 85 hundreds of dollars on 28 to 30 April 2010). With three points there is no way to estimate a weekly pattern or validate a model, so the candidates are the two benchmarks that make sense: the mean of the active days and the naive forecast (the last value).

atm3_active <- atm_clean |> filter(ATM == "ATM3", Cash > 0)
atm3_active |> select(DATE, Cash)   # the only active days

## # A tsibble: 3 x 2 [1D]
##   DATE        Cash
##   <date>     <dbl>
## 1 2010-04-28    96
## 2 2010-04-29    82
## 3 2010-04-30    85

tibble(method = c("mean of the three active days", "naive (last value)"),
       cash_hundreds = c(mean(atm3_active$Cash), last(atm3_active$Cash))) |>
  mutate(dollars = cash_hundreds * 100) |>
  knitr::kable(digits = 1)

method	cash_hundreds	dollars
mean of the three active days	87.7	8766.7
naive (last value)	85.0	8500.0

I use the mean of the three active days as the ATM3 forecast, about 8,767 dollars a day. The mean of the whole series would be close to zero and would describe a machine that was off, not one that has started operating. This forecast rests on the assumption that the last three days show how the machine will run in May; if ATM3 was only being tested, the forecast is too high, and I would want to confirm that with the bank. A flat forecast cannot show a weekly pattern, and no prediction interval is given because three observations cannot support one.

A9. Final models and the May 2010 forecast

Each selected model is refitted on all the days of its history through 30 April 2010 (cleaned as above) and used to forecast the 31 days of May. The forecast assumes that the later weekly pattern, with the low day on Tuesday, continues in May.

atm_fc <- map(working_atms, \(k) {
  sel <- atm_selected |> filter(ATM == k)
  m <- sel$.model
  # refit only the selected model, on all the days of its history through 30 April
  final_fit <- if (sel$history == "whole year") {
    fit_atm(atm_clean |> filter(ATM == k), which = m)
  } else {
    fit_atm_recent(atm_clean |> filter(ATM == k, DATE >= recent_start), which = m)
  }
  final_fit |> forecast(h = 31) |>
    hilo(level = 95) |>
    as_tibble() |>
    transmute(ATM = k, DATE, model = m, model_value = .mean, forecast = pmax(.mean, 0),   # cash cannot be negative
              lower95 = pmax(`95%`$lower, 0), upper95 = `95%`$upper)
}) |> list_rbind()

atm3_fc <- tibble(ATM = "ATM3", DATE = seq(as.Date("2010-05-01"), as.Date("2010-05-31"), by = "day"),
                  model = "mean of last 3 active days", forecast = mean(atm3_active$Cash),
                  lower95 = NA_real_, upper95 = NA_real_)

may_forecast <- bind_rows(atm_fc, atm3_fc) |>
  arrange(ATM, DATE) |>
  mutate(weekday = wday(DATE, label = TRUE), forecast_dollars = forecast * 100)

recent <- atm_clean |> filter(DATE >= as.Date("2010-03-01")) |> as_tibble()

ggplot() +
  geom_line(data = recent, aes(DATE, Cash), colour = "gray40") +
  geom_ribbon(data = may_forecast |> filter(!is.na(lower95)),
              aes(DATE, ymin = lower95, ymax = upper95), fill = "steelblue", alpha = 0.25) +
  geom_line(data = may_forecast, aes(DATE, forecast), colour = "steelblue") +
  facet_wrap(vars(ATM), ncol = 1, scales = "free_y") +
  labs(title = "May 2010 forecast by ATM (history from March 2010)",
       y = "Cash (hundreds of dollars)", x = "Date")

The table lists the forecast for the four Tuesdays of May (4, 11, 18 and 25).

may_forecast |>
  filter(ATM != "ATM3", weekday == "Tue") |>
  select(ATM, DATE, model_value, forecast, lower95, upper95) |>
  knitr::kable(digits = 1)

ATM	DATE	model_value	forecast	lower95	upper95
ATM1	2010-05-04	5.7	5.7	0.0	52.8
ATM1	2010-05-11	5.7	5.7	0.0	55.8
ATM1	2010-05-18	5.7	5.7	0.0	58.6
ATM1	2010-05-25	5.7	5.7	0.0	61.2
ATM2	2010-05-04	1.4	1.4	0.0	36.6
ATM2	2010-05-11	1.8	1.8	0.0	39.4
ATM2	2010-05-18	1.5	1.5	0.0	47.8
ATM2	2010-05-25	1.7	1.7	0.0	51.5
ATM4	2010-05-04	47.1	47.1	1.8	135.9
ATM4	2010-05-11	51.4	51.4	0.1	169.7
ATM4	2010-05-18	55.7	55.7	0.0	201.2
ATM4	2010-05-25	60.0	60.0	0.0	231.4

These are the lowest forecast days of each working ATM in the plot above, so the forecast keeps the Tuesday low day. The May Tuesday values are much lower than the April holdout forecasts in section A7 (for example 47.1 to 60.0 against 142.1 to 150.6 for ATM4), presumably because the final models also see April; the ATM4 values are still above the Tuesday values of 2.5 to 29.9 in the holdout. The full daily forecast for every ATM is in the Excel file Project1_Forecasts.xlsx, sheet ATM_May2010. Point forecasts and lower interval limits are floored at zero because cash cannot be negative; the floor shows in the lower95 column of the Tuesday table above for ATM1, ATM2 and the last two ATM4 Tuesdays.

A10. Techniques used and not used

Used: time plots, a weekday heatmap and profile plots, and ACF/PACF for the patterns; models fitted to the whole year compared with models fitted only to the days after the weekly pattern changed; same-weekday median imputation for the missing values and the outlier; Guerrero’s Box-Cox lambda; the four benchmarks; ETS with and without a Box-Cox transformation; STL with ETS; automatic, full-search and manual seasonal ARIMA; holdout accuracy (RMSE and MAE) and the Ljung-Box test to choose between them.

Tried and not kept:

• Mean, naive and drift forecasts ignore the weekly pattern. Naive and drift are far behind everything else, and the mean is far behind for ATM1 and ATM2, but they are the natural benchmarks and show how much the seasonal models add.
• Box-Cox transformations are used only where they pass the selection rule: the selected ATM4 model uses one, while the selected ATM1 and ATM2 models do not.
• Moving the earlier days two days earlier was considered and not used (section A3). It lines up the low days, but it would impose an unverified explanation for the February change and create dates that are not in the original records.
• Deleting the bad rows was rejected because the series needs a value for every day to keep its regular spacing.
• The mean of the whole year was rejected for ATM3 because it would describe a machine that was off.

Not used: regression models with weekday dummy variables could describe the weekly pattern directly, but this report compares the time-series methods listed above and does not add regression models. Replacing the weekday medians with a model-based interpolation was not done for the same reason. KNN imputation, which Kuhn and Johnson describe in section 3.4, was not used: the same-weekday median already summarizes the most comparable days, and only 6 values are replaced. Rolling cross-validation was not used for the ATMs; the comparison uses one holdout of exactly 31 days, the length of the May forecast.

Limitations

The forecast assumes that the later weekly pattern (Tuesday low day) continues through May. The data cannot show why the pattern changed in February, and a further change would break the forecast. In the April holdout the selected ATM1 and ATM4 models forecast the Tuesday low day far too high (section A7), and the May forecast for ATM4 may still do so; if the later pattern continues, the risk is over-supplying those machines on Tuesdays. The prediction intervals assume uncorrelated, roughly normal residuals, and the residuals for ATM1 and ATM2 are heavy-tailed, so the intervals should be read as approximate. Kuhn and Johnson (section 4.3) warn that a single test set is a weak basis for judging a model. The ATM models are chosen and reported on the same 31-day holdout, so the holdout errors are probably a little optimistic.

Part B - Residential power forecast

The file holds monthly residential power use in kilowatt hours (KWH) from January 1998 to December 2013. The goal is a monthly forecast for 2014.

B1. The data

power_raw <- read_excel("ResidentialCustomerForecastLoad-624.xlsx") |>
  mutate(month = yearmonth(`YYYY-MMM`)) |>
  select(month, KWH) |>
  as_tsibble(index = month)

# rows in the data against the number of months from the first to the last month
tibble(first = min(power_raw$month), last = max(power_raw$month), rows = nrow(power_raw),
       months_first_to_last = as.numeric(max(power_raw$month) - min(power_raw$month)) + 1) |>
  knitr::kable()

first	last	rows	months_first_to_last
1998 Jan	2013 Dec	192	192

power_raw |> filter(is.na(KWH) | KWH < 3e6)

## # A tsibble: 2 x 2 [1M]
##      month    KWH
##      <mth>  <dbl>
## 1 2008 Sep     NA
## 2 2010 Jul 770523

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

The first table shows 192 months from the first to the last month, with no month missing. The second table shows two bad values: September 2008 is missing, and July 2010 is 770,523 KWH, a visible spike downward in the time plot. The series has a strong annual pattern with a summer peak and a winter peak, and in the later years the level and the size of the seasonal swings are larger.

power_raw |> gg_season(KWH, labels = "both") +   # labels = "both" is not shown in the book (display only)
  labs(title = "Seasonal plot of residential power use", y = "KWH")

power_raw |> gg_subseries(KWH) +
  labs(title = "Seasonal subseries plot", y = "KWH")

The seasonal plot shows two peaks a year, a winter peak in January and a summer peak in July and August, with the lowest use in April and May and again in October and November. In the subseries plot, where the blue lines are the monthly means, August (about 8.3 million) and January (about 8.0 million) are the highest. The July panel also shows the outlier as one point far below everything else; the other Julys range from about 6 to 10 million, so 770,523 is most likely an error. The recent increase in level is not spread evenly across the year: since about 2010 the lines for January, February, June, August and December climb the most, April, May and November rise somewhat, and March and October stay roughly flat. The extra use is largest in the winter and summer months, which makes the seasonal swings look larger in the raw plot; the data has no temperature or other information that would explain why.

B2. Missing value and outlier

As in Part A, nothing is deleted, and the very low July value is treated as a possible recording error before anything else, as Kuhn and Johnson advise (section 3.3). Each bad month is replaced by the value of the same calendar month in the previous year. That keeps the seasonal pattern, and it uses only earlier data, so the replacement is available at every forecast origin used later in the cross-validation and nothing from after the month is used. The first table below gives the replaced values, and the second gives the previous and the following year of each month for comparison.

sep_2008 <- power_raw$KWH[power_raw$month == yearmonth("2007 Sep")]   # the same month one year earlier
jul_2010 <- power_raw$KWH[power_raw$month == yearmonth("2009 Jul")]

# replace the two bad months with the values computed above (data preparation)
power <- power_raw |>
  mutate(KWH_original = KWH,
         KWH = case_when(month == yearmonth("2008 Sep") ~ sep_2008,
                         month == yearmonth("2010 Jul") ~ jul_2010,
                         TRUE ~ KWH))

# the same series with the two bad months set to missing, used to score the forecasts on observed values only
power_observed <- power_raw |>
  mutate(KWH = if_else(month == yearmonth("2010 Jul"), NA_real_, KWH))

# the two replaced months, then the previous and following years behind each replacement
power |> filter(month %in% c(yearmonth("2008 Sep"), yearmonth("2010 Jul"))) |>
  as_tibble() |>
  knitr::kable(digits = 0, format.args = list(big.mark = ","))

month	KWH	KWH_original
2008 Sep	7,666,970	NA
2010 Jul	7,713,260	770,523

power_raw |> filter(month %in% yearmonth(c("2007 Sep", "2009 Sep", "2009 Jul", "2011 Jul"))) |>
  as_tibble() |>
  knitr::kable(digits = 0, format.args = list(big.mark = ","))

month	KWH
2007 Sep	7,666,970
2009 Jul	7,713,260
2009 Sep	7,583,146
2011 Jul	10,093,343

September 2007 and September 2009 are close, so the September replacement is reasonable. For July the following year (2011) is much higher than the previous year (2009), so the July replacement is probably too low and is only an estimate. The 770,523 value is a suspected error, not a confirmed typo, so it is replaced with a seasonal estimate and not with a corrected version of itself.

B3. Decomposition and autocorrelation

The cleaned series is decomposed with STL, and its ACF is shown, to see how strong the trend and the seasonal pattern are and how much differencing is needed.

power |> model(STL(KWH ~ season(window = 11), robust = TRUE)) |>
  components() |> autoplot() +
  labs(title = "STL decomposition of residential KWH")

power |> features(KWH, feat_stl) |> select(trend_strength, seasonal_strength_year) |>
  knitr::kable(digits = 2)

trend_strength	seasonal_strength_year
0.56	0.87

power |> ACF(KWH, lag_max = 48) |> autoplot() +
  labs(title = "ACF of residential KWH")

The seasonal component is strong (0.87) and the trend is moderate (0.56). The trend stays between about 5.9 and 6.5 million until 2009, climbs to about 7.4 million in 2011 and 2012, and ends near 7.3 million. The seasonal component keeps about the same size over the whole period in the STL plot, so the larger swings in the raw plot after 2010 appear to come mostly from the trend and from larger remainders, not from a growing seasonal pattern. The evidence for multiplicative seasonality is modest, so additive and multiplicative ETS forms and a Box-Cox version are compared on accuracy and not assumed. The remainder is largest at the end of the series, in December 2013 (about 9.6 million in the seasonal plot in B1, while the earlier Decembers are all below about 8.3 million). Nothing shows that this value is an error, so it is kept as a real observation. The ACF has a scalloped shape with large positive spikes at lags 12, 24, 36 and 48 that shrink only slowly (from about 0.75 to 0.5), and positive spikes of about 0.5 to 0.6 at lags 6 and 18 that reflect the two peaks a year (the spikes at lags 3, 9 and 15 are strongly negative, about -0.45, for the same reason). That slow decay at the seasonal lags is what a trending, seasonal series looks like, so a seasonal difference, and possibly a first difference, will be needed.

B4. Transformation, stationarity and the training and test split

The last 12 months (2013) are held out, so the models are trained on 1998 to 2012.

power_train <- power |> filter(year(month) <= 2012)
power_test <- power |> filter(year(month) == 2013)

power_train |> features(KWH, guerrero) |> knitr::kable(digits = 3)

lambda_guerrero
-0.169

power_train |> features(KWH, list(unitroot_kpss, unitroot_nsdiffs, unitroot_ndiffs)) |>
  knitr::kable(digits = 3)

kpss_stat	kpss_pvalue	nsdiffs	ndiffs
1.17	0.01	1	1

The Guerrero lambda is negative, which would be a strong transformation, so Box-Cox models are tested rather than assumed. The series needs one seasonal difference. After that I check whether another difference is needed and look at the ACF and PACF to choose manual ARIMA candidates.

power_train |> features(difference(KWH, 12), unitroot_kpss) |> knitr::kable(digits = 3)

kpss_stat	kpss_pvalue
0.101	0.1

power_train |> gg_tsdisplay(difference(KWH, 12), plot_type = "partial", lag_max = 36) +
  labs(title = "Residential KWH after one seasonal difference (lag 12)")

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

After the seasonal difference the KPSS test no longer rejects stationarity, so no further difference is needed. The ACF has a large negative spike at lag 12, and the PACF has negative spikes at lags 12 and 24 that shrink, which is the pattern of a seasonal moving average term, so (0,1,1) in the seasonal part. The warnings about removed rows are the first 12 values lost to the seasonal difference. Both plots also show a positive spike at lag 1, which a non-seasonal AR(1) or MA(1) term can capture, and a smaller one at lag 3. The manual candidates are ARIMA(1,0,0)(0,1,1)[12], ARIMA(0,0,1)(0,1,1)[12] and ARIMA(2,0,0)(0,1,1)[12]. ARIMA(1,0,0)(1,1,0)[12] is added as a contrast, a seasonal AR term in place of the seasonal MA term. They are fitted next to the automatic choices.

B5. Fitting and comparing the models

What the exploration led to

What the exploration shows	Where it is shown	What it led to
A strong annual pattern with two peaks a year	B1 seasonal and subseries plots, B3 seasonal strength	Seasonal naive as the benchmark, and seasonal ETS and seasonal ARIMA models
The level rises after 2009 and the swings look larger in the later years	B1 time plot, B3 STL plot	Box-Cox ETS compared with plain ETS, and drift terms in the ARIMA models
One seasonal difference is needed, and the ACF and PACF after it show a seasonal moving average pattern	B4 tests and plots	Manual ARIMA candidates with (0,1,1) in the seasonal part, next to the automatic and full-search ARIMA
One missing month and one suspected recording error	B1 tables, B2	Replaced with the same month of the previous year, and not scored as observed values

lambda_power <- power_train |> features(KWH, guerrero) |> pull(lambda_guerrero)

power_fit <- power_train |> model(
  mean = MEAN(KWH),
  naive = NAIVE(KWH),
  snaive = SNAIVE(KWH),
  drift = RW(KWH ~ drift()),
  ets = ETS(KWH),
  ets_boxcox = ETS(box_cox(KWH, lambda_power)),   # Box-Cox inside ETS() is only mentioned in exercise 8.13 of the book
  stl_ets = decomposition_model(STL(KWH ~ season(window = 11), robust = TRUE),
                                ETS(season_adjust ~ season("N"))),
  arima_auto = ARIMA(KWH),
  arima_search = ARIMA(KWH, stepwise = FALSE, approx = FALSE),
  arima_100_011 = ARIMA(KWH ~ pdq(1, 0, 0) + PDQ(0, 1, 1)),
  arima_001_011 = ARIMA(KWH ~ pdq(0, 0, 1) + PDQ(0, 1, 1)),
  arima_200_011 = ARIMA(KWH ~ pdq(2, 0, 0) + PDQ(0, 1, 1)),
  arima_100_110 = ARIMA(KWH ~ pdq(1, 0, 0) + PDQ(1, 1, 0))
)

# list the ARIMA and ETS specifications as a table (display only)
power_fit |> select(starts_with("arima"), ets, ets_boxcox) |>
  pivot_longer(everything(), names_to = "model", values_to = "specification") |>
  knitr::kable()

model	specification
arima_auto	<ARIMA(0,0,1)(2,1,0)[12] w/ drift>
arima_search	<ARIMA(3,0,1)(2,1,0)[12] w/ drift>
arima_100_011	<ARIMA(1,0,0)(0,1,1)[12] w/ drift>
arima_001_011	<ARIMA(0,0,1)(0,1,1)[12] w/ drift>
arima_200_011	<ARIMA(2,0,0)(0,1,1)[12] w/ drift>
arima_100_110	<ARIMA(1,0,0)(1,1,0)[12] w/ drift>
ets	<ETS(M,N,M)>
ets_boxcox	<ETS(M,N,A)>

power_accuracy <- accuracy(forecast(power_fit, h = 12), power) |>
  as_tibble() |> select(.model, RMSE, MAE, MASE) |>
  left_join(ljung_table(power_fit, 24), by = ".model")

power_accuracy |>
  arrange(RMSE) |>
  select(model = .model, RMSE, MAE, MASE, ljung_box_p = lb_pvalue) |>
  knitr::kable(digits = 3, format.args = list(big.mark = ","))

model	RMSE	MAE	MASE	ljung_box_p
stl_ets	947,371.5	583,774.9	0.961	0.046
arima_search	949,335.1	637,187.3	1.049	0.994
arima_auto	968,992.5	655,243.8	1.078	0.179
arima_100_011	980,554.1	600,155.0	0.988	0.432
arima_200_011	982,271.4	601,409.0	0.990	0.346
arima_001_011	994,084.6	613,661.1	1.010	0.465
arima_100_110	996,948.2	719,593.0	1.184	0.002
snaive	1,035,537.6	618,605.6	1.018	0.000
ets_boxcox	1,049,814.2	686,014.3	1.129	0.014
ets	1,111,268.4	722,113.2	1.189	0.029
naive	1,861,178.6	1,501,420.0	2.471	0.000
drift	1,865,709.5	1,502,832.8	2.474	0.000
mean	1,938,591.6	1,531,580.0	2.521	0.000

The models without a seasonal component (mean, naive, drift) are far behind. RMSE and MAE do not rank every model the same way: the ARIMA(1,0,0)(0,1,1) model has a lower MAE (600,155) than the full-search ARIMA (637,187), although its RMSE is higher (980,554 against 949,335). RMSE gives more weight to large misses, so a model with a few large errors ranks lower on RMSE than on MAE. The seasonal naive model, which just repeats last year, is strong for this series, and the best models improve on it by only about 5 to 9 percent in RMSE. A single test year is a fragile basis for choosing among models that close, as Kuhn and Johnson also say of a single test set (section 4.3), and 2013 ends with an unusually high December (section B3). The holdout is therefore repeated across several years with time series cross-validation: the training data grows by 12 months at a time (starting with 120 months) and each fit forecasts the next 12 months. The Box-Cox ETS is left out of the cross-validation because its lambda would have to be estimated again in every fold; it is compared on the 2013 holdout only, and the two bad months are not scored.

power_cv <- power |>
  stretch_tsibble(.init = 120, .step = 12) |>
  model(snaive = SNAIVE(KWH),
        ets = ETS(KWH),
        stl_ets = decomposition_model(STL(KWH ~ season(window = 11), robust = TRUE),
                                      ETS(season_adjust ~ season("N"))),
        arima_auto = ARIMA(KWH),
        arima_search = ARIMA(KWH, stepwise = FALSE, approx = FALSE),
        arima_100_011 = ARIMA(KWH ~ pdq(1, 0, 0) + PDQ(0, 1, 1)),
        arima_001_011 = ARIMA(KWH ~ pdq(0, 0, 1) + PDQ(0, 1, 1)),
        arima_200_011 = ARIMA(KWH ~ pdq(2, 0, 0) + PDQ(0, 1, 1)),
        arima_100_110 = ARIMA(KWH ~ pdq(1, 0, 0) + PDQ(1, 1, 0))) |>
  forecast(h = 12) |>
  accuracy(power_observed) |>   # the two bad months are not scored
  as_tibble() |> select(.model, cv_RMSE = RMSE, cv_MAE = MAE)

## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 12 observations are missing between 2014 Jan and 2014 Dec

power_cv |>
  mutate(percent_above_best = 100 * (cv_RMSE / min(cv_RMSE) - 1)) |>
  arrange(cv_RMSE) |>
  select(model = .model, cv_RMSE, cv_MAE, percent_above_best) |>
  knitr::kable(digits = 1, format.args = list(big.mark = ","))

model	cv_RMSE	cv_MAE	percent_above_best
arima_auto	779,327.5	553,316.6	0.0
arima_100_011	804,664.4	561,626.5	3.3
arima_search	805,174.7	577,103.6	3.3
arima_001_011	805,464.8	556,452.8	3.4
arima_200_011	807,585.6	563,291.5	3.6
arima_100_110	808,341.6	600,242.9	3.7
stl_ets	848,696.5	605,120.6	8.9
ets	890,273.1	626,391.0	14.2
snaive	939,756.7	695,725.9	20.6

The last stretch of data ends in December 2013 and has no future data to test on, so it is ignored; that is what the warning above says.

Selection rule. For Part B the choice rests on the cross-validated RMSE, because it uses several forecast origins instead of one. Among the models whose residuals also pass the Ljung-Box test in the 2013 fit, I take the one with the lowest cross-validated RMSE. The next-best model in the table above is about 3% higher, so the 1% tolerance step used in Part A would not change this choice.

power_selected <- power_accuracy |>
  left_join(power_cv, by = ".model") |>
  filter(lb_pvalue >= 0.05, !is.na(cv_RMSE)) |>   # keep models that pass the residual test
  slice_min(cv_RMSE, n = 1)   # lowest cross-validated RMSE

power_selected |> select(model = .model, cv_RMSE, holdout_RMSE = RMSE, ljung_box_p = lb_pvalue) |>
  knitr::kable(digits = 3, format.args = list(big.mark = ","))

model	cv_RMSE	holdout_RMSE	ljung_box_p
arima_auto	779,327.5	968,992.5	0.179

power_fit |> select(all_of(power_selected$.model)) |> gg_tsresiduals() +
  labs(title = paste("Residuals of", power_selected$.model))

The residuals of the selected model, fitted on the training years, have mean near zero and no obvious trend, and the histogram is close to bell-shaped. The ACF has one spike clearly outside the dashed lines, at lag 3, which I treat as a small remaining weakness; the Ljung-Box test over 24 lags does not reject at the 5% level (p = 0.179). The spread of the residuals is larger from 2010 on, so the prediction intervals are probably a little too narrow for the most recent years.

The single 2013 holdout and the cross-validation do not agree on which model is best: STL with ETS has the lowest holdout RMSE but its residuals narrowly fail the Ljung-Box test (p = 0.046), the full-search ARIMA is next, and the automatic ARIMA wins the cross-validation. The cross-validation covers several forecast years and gives a more reliable comparison than a single split, so it decides, and the other ARIMA models are within about 4% of the automatic one in the cross-validation table. That is not a clear difference.

B6. The 2014 forecast

The selected model is refitted on all data through December 2013 and forecasts the twelve months of 2014. The model is the automatic ARIMA() procedure, so its order is chosen again on the longer data; the specification it picks is printed first.

# the selected procedure only; its order is chosen again on the longer data
stopifnot(power_selected$.model == "arima_auto")
power_final <- power |> model(arima_auto = ARIMA(KWH))

power_final   # the specification the selected model chooses on all the data

## # A mable: 1 x 1
##                           arima_auto
##                              <model>
## 1 <ARIMA(0,0,4)(2,1,0)[12] w/ drift>

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

power_forecast <- power_fc |>
  hilo(level = 95) |>
  as_tibble() |>
  transmute(month, model = power_selected$.model, forecast_KWH = .mean,
            lower95 = `95%`$lower, upper95 = `95%`$upper)

power_forecast |> knitr::kable(digits = 0, format.args = list(big.mark = ","))

month	model	forecast_KWH	lower95	upper95
2014 Jan	arima_auto	10,291,438	9,084,135	11,498,741
2014 Feb	arima_auto	8,607,148	7,335,158	9,879,137
2014 Mar	arima_auto	7,159,784	5,887,058	8,432,509
2014 Apr	arima_auto	5,833,040	4,533,397	7,132,683
2014 May	arima_auto	5,942,376	4,639,971	7,244,781
2014 Jun	arima_auto	8,198,346	6,895,941	9,500,751
2014 Jul	arima_auto	9,474,812	8,172,407	10,777,217
2014 Aug	arima_auto	9,969,550	8,667,145	11,271,955
2014 Sep	arima_auto	8,485,269	7,182,864	9,787,674
2014 Oct	arima_auto	5,870,010	4,567,605	7,172,415
2014 Nov	arima_auto	6,145,967	4,843,562	7,448,372
2014 Dec	arima_auto	8,427,924	7,125,519	9,730,329

power_fc |> autoplot(power |> filter(year(month) >= 2008), level = 95) +
  labs(title = "Residential power use: 2014 forecast", y = "KWH", x = "Month")

# average monthly use in the last observed year and in the forecast year
tibble(year = c("2013 (actual)", "2014 (forecast)"),
       mean_KWH = c(mean(power_test$KWH), mean(power_forecast$forecast_KWH))) |>
  knitr::kable(digits = 0, format.args = list(big.mark = ","))

year	mean_KWH
2013 (actual)	7,618,335
2014 (forecast)	7,867,139

# residual check of the final fit on all the data (the earlier residual plot is from the fit on the training years)
power_final |> gg_tsresiduals() + labs(title = "Residuals of the final fit on all the data")

ljung_table(power_final, 24) |> knitr::kable(digits = 3)

.model	lb_pvalue
arima_auto	0.711

The forecast keeps the shape of the earlier years: a peak in January, another in July and August, and the lowest use in April, May, October and November. The 2014 average is higher than the 2013 average in the last table, which is consistent with the drift term of the model carrying the rise in level into 2014. The 95% intervals are about 2.4 to 2.6 million KWH wide, for example January from about 9.1 to 11.5 million.

Limitations of the 2014 forecast. The replacements for September 2008 and July 2010 are estimates, and December 2013 is an unusually high last value that carries weight in the refit. The order chosen on all the data, ARIMA(0,0,4)(2,1,0)[12] with drift, is not the one validated on the 2013 holdout, ARIMA(0,0,1)(2,1,0)[12], because the automatic procedure chooses the order again on the longer data. The residuals of the final fit in the plot above show no spike outside the dashed lines and the Ljung-Box test does not reject (p = 0.711), but their spread is larger from 2010 on and the last residual, for December 2013, is the largest, about 2.5 million, so the intervals are probably too narrow. The model uses only the history of the series, so it cannot anticipate weather or economic change.

B7. Techniques used and not used

Used: seasonal and subseries plots, seasonal differencing with KPSS and ACF/PACF to pick manual candidates, replacement of the missing month and the outlier with the same month of the previous year, Guerrero’s lambda, the benchmarks, ETS (with and without Box-Cox), STL with ETS, automatic, full-search and manual seasonal ARIMA, a one-year holdout and time series cross-validation, and the Ljung-Box test.

Tried and not kept: the Box-Cox ETS model (the Guerrero lambda is negative (-0.169), which would be a strong transformation; in the 2013 holdout table the Box-Cox ETS has a lower RMSE than the plain ETS, 1,049,814 against 1,111,268, but its residuals fail the Ljung-Box test (p = 0.014) and it stays behind the seasonal ARIMA models); the benchmarks other than seasonal naive, which cannot follow the seasonal pattern.

Not used: model-based interpolation for the two bad months, and regression models for the upward drift in level after 2010. The report compares the time-series methods listed above and does not add either. The ARIMA models include a drift term when the data call for it, and that is how the rise in level is carried into 2014.

Part C - Water flow (bonus)

Two files hold water-flow readings from two pipes. The goal is to put each on an hourly sequence, taking the mean of the readings within an hour, to decide whether the hourly series is stationary and can be forecast, and if so to forecast one week (168 hours) ahead.

C1. Reading the data

read_pipe <- function(file) {
  read_excel(file) |>
    rename(date_time = `Date Time`, flow = WaterFlow) |>
    mutate(date_time = as.POSIXct(round(date_time * 86400), origin = "1899-12-30", tz = "UTC"))   # Excel day numbers to date-times
}

pipes_raw <- bind_rows(pipe1 = read_pipe("Waterflow_Pipe1.xlsx"),
                       pipe2 = read_pipe("Waterflow_Pipe2.xlsx"), .id = "pipe")

# size, date range, missing values and spacing (in minutes) of the readings for each pipe
pipes_raw |>
  group_by(pipe) |>
  mutate(gap_minutes = as.numeric(difftime(date_time, lag(date_time), units = "mins"))) |>
  summarise(rows = n(), first = format(min(date_time), "%Y-%m-%d %H:%M"), last = format(max(date_time), "%Y-%m-%d %H:%M"),
            missing = sum(is.na(flow)), dup_times = sum(duplicated(date_time)),
            min_gap = min(gap_minutes, na.rm = TRUE), med_gap = median(gap_minutes, na.rm = TRUE),
            max_gap = max(gap_minutes, na.rm = TRUE)) |>
  knitr::kable(digits = 1)

pipe	rows	first	last	missing	dup_times	min_gap	med_gap	max_gap
pipe1	1000	2015-10-23 00:24	2015-11-01 23:35	0	0	0	10	82.1
pipe2	1000	2015-10-23 01:00	2015-12-03 16:00	0	0	60	60	60.0

Both files have 1,000 readings, with no missing values and no repeated times. The readings of Pipe 1 are spaced irregularly, from almost no time apart to 82.1 minutes, with a median of 10 minutes, and they cover about ten days (23 October to 1 November 2015). The readings of Pipe 2 are exactly 60 minutes apart from 23 October to 3 December 2015, so it is already an hourly series. The two pipes therefore have different time stamps and different lengths of history.

C2. The raw readings

pipes_raw |>
  group_by(pipe) |>
  summarise(mean = mean(flow), sd = sd(flow), min = min(flow), median = median(flow), max = max(flow)) |>
  knitr::kable(digits = 1)

pipe	mean	sd	min	median	max
pipe1	19.9	8.2	1.1	19.9	38.9
pipe2	39.6	16.0	1.9	39.7	78.3

ggplot(pipes_raw, aes(date_time, flow)) +
  geom_line(colour = "gray60", linewidth = 0.3) +
  geom_point(size = 0.6, colour = "steelblue") +
  facet_wrap(vars(pipe), scales = "free_x", ncol = 1) +
  labs(title = "Raw readings of each pipe, as recorded", x = "Date and time", y = "Water flow")

The average flow of Pipe 1 is 19.9 and that of Pipe 2 is 39.6, so Pipe 2 carries about twice as much, and it also varies more (standard deviation 16.0 against 8.2). All readings are positive, and none is extreme: the largest is 38.9 for Pipe 1 and 78.3 for Pipe 2. In the graph the readings of both pipes scatter around a steady level, with no visible trend or cycle.

C3. Putting the readings on an hourly sequence

The hour of each reading is found by rounding its time down to the hour, and the readings in the same hour are averaged. For Pipe 1 the table counts how many hours have 1, 2, 3 and so on readings, and lists the hours of its span that have no reading.

pipe1_hourly_raw <- pipes_raw |>
  filter(pipe == "pipe1") |>
  mutate(hour = floor_date(date_time, "hour")) |>
  group_by(hour) |>
  summarise(flow = mean(flow), readings = n(), .groups = "drop")

pipe1_hourly_raw |> count(readings, name = "hours") |> knitr::kable()

readings	hours
1	18
2	31
3	40
4	46
5	47
6	26
7	11
8	8
9	6
10	3

# every hour from the first to the last reading, and the hours without a reading
pipe1_hours <- tibble(hour = seq(min(pipe1_hourly_raw$hour), max(pipe1_hourly_raw$hour), by = "hour"))
pipe1_hours |> anti_join(pipe1_hourly_raw, by = "hour") |> knitr::kable()

hour
2015-10-27 17:00:00
2015-10-28 00:00:00
2015-11-01 05:00:00
2015-11-01 09:00:00

Most hours of Pipe 1 have between 3 and 6 readings, and the number runs from 1 to 10 (first table). Four hours in its span (23 October to 1 November) have no reading (second table), and none of them is next to another missing hour.

# an isolated missing hour gets the average of the hour before and the hour after
pipe1_filled <- pipe1_hours |>
  left_join(pipe1_hourly_raw, by = "hour") |>
  mutate(before = lag(flow), after = lead(flow),
         filled = is.na(flow),
         flow = if_else(filled, (before + after) / 2, flow))

pipe1_filled |> filter(filled) |> select(hour, before, after, flow) |> knitr::kable(digits = 1)

hour	before	after	flow
2015-10-27 17:00:00	25.9	19.1	22.5
2015-10-28 00:00:00	13.5	24.1	18.8
2015-11-01 05:00:00	16.9	16.8	16.8
2015-11-01 09:00:00	15.1	15.7	15.4

pipe1_ts <- pipe1_filled |> select(hour, flow) |> as_tsibble(index = hour)
# the same series with the filled hours set to missing, so they are not scored as observed values
pipe1_scoring <- pipe1_filled |>
  mutate(flow = if_else(filled, NA_real_, flow)) |>
  select(hour, flow) |>
  as_tsibble(index = hour)
pipe2_ts <- pipes_raw |>
  filter(pipe == "pipe2") |>
  transmute(hour = floor_date(date_time, "hour"), flow) |>
  as_tsibble(index = hour)

# hours in the span against hours with a value, for each pipe
tibble(pipe = c("pipe1", "pipe2"),
       hours_in_span = c(length(seq(min(pipe1_ts$hour), max(pipe1_ts$hour), by = "hour")),
                         length(seq(min(pipe2_ts$hour), max(pipe2_ts$hour), by = "hour"))),
       hours_with_value = c(nrow(pipe1_ts), nrow(pipe2_ts))) |>
  knitr::kable()

pipe	hours_in_span	hours_with_value
pipe1	240	240
pipe2	1000	1000

Each missing hour is filled with the average of the hour before and the hour after; the first table shows the three values for each. After the filling, Pipe 1 has a value in every one of its 240 hours, and Pipe 2 has a value in every one of its 1,000 hours (second table), so both series are complete and regular. The hourly means of Pipe 1 rest on different numbers of readings, so an hour with one or two readings may give a less precise hourly average than one with six; I keep this in mind but do not adjust for it.

C4. Is the hourly series stationary?

The graphs show the two hourly series with their averages, and the autocorrelation of each (48 hours for Pipe 1, a week for Pipe 2).

hourly_both <- bind_rows(pipe1 = as_tibble(pipe1_ts), pipe2 = as_tibble(pipe2_ts), .id = "pipe")

ggplot(hourly_both, aes(hour, flow)) +
  geom_line(colour = "gray50", linewidth = 0.4) +
  geom_hline(data = hourly_both |> group_by(pipe) |> summarise(average = mean(flow)),   # not in the book; draws each series' average (display only)
             aes(yintercept = average), colour = "firebrick", linetype = "dashed") +
  facet_wrap(vars(pipe), scales = "free", ncol = 1) +
  labs(title = "Hourly water flow with each series' average (red dashed)", x = "Hour", y = "Water flow")

pipe1_ts |> ACF(flow, lag_max = 48) |> autoplot() + labs(title = "Pipe 1: ACF, 48 hours")

pipe2_ts |> ACF(flow, lag_max = 168) |> autoplot() + labs(title = "Pipe 2: ACF, 168 hours")

Both hourly series move around a constant average with no trend, and their spread looks the same from start to end. The ACF of Pipe 1 stays at or inside the dashed lines at every lag up to 48 hours. The ACF of Pipe 2 has a few spikes slightly outside the lines, not at a regular spacing; with 168 lags, about 1 in 20 can fall outside by chance.

The tables give the average and standard deviation of the hourly values, the KPSS statistic and the number of differences suggested, the Ljung-Box test for white noise at lags 10 and 24, and the strength of a 24-hour pattern.

pipes_ts <- list(pipe1 = pipe1_ts, pipe2 = pipe2_ts)

# average and standard deviation of the hourly values of each pipe
map(pipes_ts, \(x) x |> as_tibble() |> summarise(mean = mean(flow), sd = sd(flow))) |>
  list_rbind(names_to = "pipe") |> knitr::kable(digits = 1)

pipe	mean	sd
pipe1	19.9	4.2
pipe2	39.6	16.0

map(pipes_ts, \(x) x |> features(flow, list(unitroot_kpss, unitroot_ndiffs))) |>
  list_rbind(names_to = "pipe") |> knitr::kable(digits = 3)

pipe	kpss_stat	kpss_pvalue	ndiffs
pipe1	0.193	0.1	0
pipe2	0.105	0.1	0

map(c(10, 24), \(lag) {
  map(pipes_ts, \(x) x |> features(flow, ljung_box, lag = lag)) |> list_rbind(names_to = "pipe") |>
    mutate(lag = lag)
}) |> list_rbind() |> select(pipe, lag, lb_stat, lb_pvalue) |> knitr::kable(digits = 3)

pipe	lag	lb_stat	lb_pvalue
pipe1	10	7.755	0.653
pipe2	10	6.458	0.775
pipe1	24	22.034	0.577
pipe2	24	32.474	0.116

map(pipes_ts, \(x) x |> features(flow, feat_stl)) |>
  list_rbind(names_to = "pipe") |>
  select(pipe, trend_strength, seasonal_strength_day) |> knitr::kable(digits = 2)

pipe	trend_strength	seasonal_strength_day
pipe1	0.04	0.21
pipe2	0.08	0.23

The KPSS test does not reject stationarity for either pipe (p-values of 0.1 or more), and no differencing is suggested. The Ljung-Box tests do not reject white noise at lag 10 or 24 for either pipe (all p-values above 0.1), so the tests did not detect autocorrelation at those lags. The trend strength is close to zero (0.04 and 0.08) and the strength of a 24-hour pattern is weak (0.21 and 0.23). The diagnostics are consistent with stationarity, and both series look like white noise around a constant mean: there is no trend or cycle to forecast, only the level. A series like this can be forecast, but only as a level: the forecast is the average.

C5. Fitting and comparing models

What the exploration led to

What the exploration shows	Where it is shown	What it led to
Pipe 1 readings are irregular and several share an hour; Pipe 2 is already hourly	C1 to C3 tables	Hourly means for both pipes
Four isolated missing hours in Pipe 1	C3 tables	Filled from the neighboring hours, and not scored as observed values
Little trend, weak 24-hour pattern, little detected autocorrelation, and KPSS does not reject stationarity	C4 tables and plots	The mean and small ARIMA models as the main candidates, next to the benchmarks and ETS

The same kinds of models as before are compared: the four benchmarks, ETS, the automatic ARIMA, and three small ARIMA models without seasonal terms. ARIMA(0,0,0) is white noise around a constant, and ARIMA(1,0,0) and ARIMA(0,0,1) add one autocorrelation term. The holdout is the last 48 hours for Pipe 1, which has only 240 hours of history, and the last 168 hours (the length of the forecast) for Pipe 2. The models are compared by RMSE and MAE on the holdout and by the Ljung-Box test at lag 24 on the residuals, and the selection rule is the Part A rule (models that pass the residual test, lowest holdout RMSE) with one addition from Kuhn and Johnson (section 4.6), who advise choosing a simpler model whose performance is within a small tolerance of the best. Here the tolerance is 1%, a threshold I chose (their example accepts a 4% loss), and the rule is applied across model families, not to the tuning parameters of one model. Among the models that pass the residual test, those within 1% of the lowest RMSE are kept, and the simplest of them is chosen; the order from simplest to most complex is the plain mean, ARIMA(0,0,0), ETS, the ARIMA models with one term, and the automatic ARIMA.

holdout_hours <- c(pipe1 = 48, pipe2 = 168)

fit_pipe <- function(series, which = NULL) {
  specs <- list(
    mean = MEAN(flow),
    naive = NAIVE(flow),
    snaive = SNAIVE(flow),
    drift = RW(flow ~ drift()),
    ets = ETS(flow),
    arima_auto = ARIMA(flow),
    arima_000 = ARIMA(flow ~ pdq(0, 0, 0) + PDQ(0, 0, 0)),
    arima_100 = ARIMA(flow ~ pdq(1, 0, 0) + PDQ(0, 0, 0)),
    arima_001 = ARIMA(flow ~ pdq(0, 0, 1) + PDQ(0, 0, 0))
  )
  if (is.null(which)) which <- names(specs)
  do.call(model, c(list(series), specs[which]))   # fit only the models named in `which`
}

pipe_scoring <- list(pipe1 = pipe1_scoring, pipe2 = pipes_ts$pipe2)   # holdout scored on observed hours only

pipe_train <- map(names(pipes_ts), \(k) {
  x <- pipes_ts[[k]]
  x |> filter(hour <= max(hour) - hours(holdout_hours[[k]]))
}) |> set_names(names(pipes_ts))
pipe_fits <- map(pipe_train, fit_pipe)

pipe_accuracy <- map(names(pipes_ts), \(k) {
  accuracy(forecast(pipe_fits[[k]], h = holdout_hours[[k]]), pipe_scoring[[k]]) |>
    as_tibble() |> select(.model, RMSE, MAE) |>
    left_join(ljung_table(pipe_fits[[k]], 24), by = ".model") |>
    mutate(pipe = k)
}) |> list_rbind()

pipe_accuracy |> arrange(pipe, RMSE) |>
  select(pipe, model = .model, RMSE, MAE, ljung_box_p = lb_pvalue) |> knitr::kable(digits = 3)

pipe	model	RMSE	MAE	ljung_box_p
pipe1	arima_100	3.947	2.994	0.454
pipe1	arima_auto	3.947	2.994	0.454
pipe1	arima_001	3.948	2.994	0.447
pipe1	arima_000	3.951	2.999	0.592
pipe1	mean	3.951	2.999	0.592
pipe1	ets	3.951	2.999	0.592
pipe1	naive	6.162	5.112	0.000
pipe1	snaive	6.362	5.145	0.000
pipe1	drift	7.196	6.314	0.000
pipe2	arima_auto	16.638	13.676	0.235
pipe2	ets	16.733	13.749	0.079
pipe2	arima_100	16.733	13.750	0.090
pipe2	arima_000	16.733	13.750	0.079
pipe2	mean	16.733	13.750	0.079
pipe2	arima_001	16.733	13.750	0.090
pipe2	naive	16.736	13.761	0.000
pipe2	drift	16.736	13.892	0.000
pipe2	snaive	21.570	17.605	0.000

# the specifications chosen by the automatic ARIMA and ETS
map(pipe_fits, \(f) f |> select(ets, arima_auto) |> pivot_longer(everything(), names_to = "model", values_to = "specification")) |>
  list_rbind(names_to = "pipe") |> knitr::kable()

pipe	model	specification
pipe1	ets	<ETS(M,N,N)>
pipe1	arima_auto	<ARIMA(1,0,0) w/ mean>
pipe2	ets	<ETS(M,N,N)>
pipe2	arima_auto	<ARIMA(0,0,0)(0,0,1)[24] w/ mean>

For both pipes the naive, seasonal naive and drift forecasts are no better than the mean, and for Pipe 1 they are clearly worse (RMSE 6.16, 6.36 and 7.20 against 3.95). They also fail the residual test (p about 0), because following the last value or the last day adds noise to a series with little detected autocorrelation (section C4). The mean, ARIMA(0,0,0) and ETS give the same RMSE for each pipe (3.95 for Pipe 1 and 16.73 for Pipe 2); the ETS(M,N,N) model in the second table is in effect the same flat forecast. The other ARIMA models differ from them by less than 1% (3.947 against 3.951 for Pipe 1, and 16.638 against 16.733 for Pipe 2). The automatic ARIMA chose ARIMA(1,0,0) with a mean for Pipe 1 and ARIMA(0,0,0)(0,0,1)[24] with a mean for Pipe 2. Two of the four filled hours of Pipe 1 (1 November at 05:00 and 09:00) fall in its 48-hour holdout; they are not scored, because they are filled values and not observations.

# the candidate models from the simplest to the most complex
simplicity_order <- c("mean", "arima_000", "ets", "arima_100", "arima_001", "arima_auto",
                      "drift", "naive", "snaive")

pipe_candidates <- pipe_accuracy |>
  filter(lb_pvalue >= 0.05) |>   # keep models that pass the residual test
  group_by(pipe) |>
  mutate(percent_above_best = 100 * (RMSE / min(RMSE) - 1)) |>
  ungroup() |>
  filter(percent_above_best <= 1)   # within 1% of the lowest holdout RMSE

pipe_candidates |>
  arrange(pipe, RMSE) |>
  select(pipe, model = .model, RMSE, percent_above_best, ljung_box_p = lb_pvalue) |>
  knitr::kable(digits = 3)

pipe	model	RMSE	percent_above_best	ljung_box_p
pipe1	arima_100	3.947	0.000	0.454
pipe1	arima_auto	3.947	0.000	0.454
pipe1	arima_001	3.948	0.001	0.447
pipe1	arima_000	3.951	0.086	0.592
pipe1	mean	3.951	0.086	0.592
pipe1	ets	3.951	0.088	0.592
pipe2	arima_auto	16.638	0.000	0.235
pipe2	ets	16.733	0.569	0.079
pipe2	arima_100	16.733	0.569	0.090
pipe2	arima_000	16.733	0.569	0.079
pipe2	mean	16.733	0.569	0.079
pipe2	arima_001	16.733	0.570	0.090

pipe_selected <- pipe_candidates |>
  group_by(pipe) |>
  slice_min(match(.model, simplicity_order), n = 1, with_ties = FALSE) |>   # the simplest of them
  ungroup()

pipe_selected |> select(pipe, model = .model, RMSE, MAE, ljung_box_p = lb_pvalue) |> knitr::kable(digits = 3)

pipe	model	RMSE	MAE	ljung_box_p
pipe1	mean	3.951	2.999	0.592
pipe2	mean	16.733	13.750	0.079

for (k in pipe_selected$pipe) {
  m <- pipe_selected |> filter(pipe == k) |> pull(.model)
  print(pipe_fits[[k]] |> select(all_of(m)) |> gg_tsresiduals() + labs(title = paste(k, "-", m)))
}

The first selection table lists the models that pass the residual test and are within 1% of the lowest RMSE; the naive, seasonal naive and drift models fail the residual test in the accuracy table above and are not in it. The lowest RMSE belongs to ARIMA(1,0,0) for Pipe 1 (3.947) and to the automatic ARIMA for Pipe 2 (16.638), but the plain mean is less than 1% above them (3.951 and 16.733), so it is the simplest model in the range and is selected for both pipes (second selection table). The ARIMA models’ advantage is within the 1% tolerance, so the simpler mean is preferred. In the residual plots the autocorrelations stay inside the dashed lines except for a few lags that are just outside, which is about what chance gives, and the histograms are close to bell-shaped. The Ljung-Box p-value of the mean model for Pipe 2 (0.079) is above the 5% limit but not by much. For a white-noise series the best forecast is the level of the series, and the plain mean estimates that level; the NAIVE() forecast, which repeats the last value, is no better here, and clearly worse for Pipe 1 (accuracy table above), because it carries the noise of one hour into every later hour.

C6. The one-week forecast

The selected model is refitted on all the hourly data of each pipe and forecasts the 168 hours that follow the last hour of the series. The table gives the first, middle and last hour of each forecast, and the graph shows the last seven days of history with the forecast and its 95% interval.

pipe_forecast <- map(names(pipes_ts), \(k) {
  m <- pipe_selected |> filter(pipe == k) |> pull(.model)
  fit_pipe(pipes_ts[[k]], which = m) |> forecast(h = 168) |>
    hilo(level = 95) |>
    as_tibble() |>
    transmute(pipe = k, hour, model = m, forecast_flow = .mean, lower95 = `95%`$lower, upper95 = `95%`$upper)
}) |> list_rbind()

# first, middle and last hour of each week-ahead forecast
pipe_forecast |> group_by(pipe) |> slice(c(1, 84, 168)) |> ungroup() |> knitr::kable(digits = 1)

pipe	hour	model	forecast_flow	lower95	upper95
pipe1	2015-11-02 00:00:00	mean	19.9	11.6	28.2
pipe1	2015-11-05 11:00:00	mean	19.9	11.6	28.2
pipe1	2015-11-08 23:00:00	mean	19.9	11.6	28.2
pipe2	2015-12-03 17:00:00	mean	39.6	8.1	71.0
pipe2	2015-12-07 04:00:00	mean	39.6	8.1	71.0
pipe2	2015-12-10 16:00:00	mean	39.6	8.1	71.0

recent_hours <- hourly_both |> group_by(pipe) |> filter(hour > max(hour) - days(7)) |> ungroup()

ggplot() +
  geom_line(data = recent_hours, aes(hour, flow), colour = "gray40") +
  geom_ribbon(data = pipe_forecast, aes(hour, ymin = lower95, ymax = upper95), fill = "steelblue", alpha = 0.25) +
  geom_line(data = pipe_forecast, aes(hour, forecast_flow), colour = "steelblue") +
  facet_wrap(vars(pipe), scales = "free", ncol = 1) +
  labs(title = "One-week forecast for each pipe (last 7 days of history)", x = "Hour", y = "Water flow")

For both pipes the forecast is a flat line at the average of the series (19.9 for Pipe 1 and 39.6 for Pipe 2, the averages in the first table of section C4), with a 95% interval of about 11.6 to 28.2 for Pipe 1 and about 8.1 to 71.0 for Pipe 2. The interval is about 1.96 standard deviations of the hourly values either side of the average (the standard deviations are in the first table of section C4); for Pipe 1 it is narrower than the spread of the raw readings in C2 because each hourly value is an average of several readings. This is what white noise implies: the forecast says where the level is, and nothing about which hour will be high or low. The forecasts of both pipes are in the Excel file, sheets Water_Pipe1_Week and Water_Pipe2_Week.

Limitations. The forecast assumes that the level stays where it was during the history, about ten days for Pipe 1 and about six weeks for Pipe 2; a change in demand or in the system would not be seen. Pipe 1 has only 240 hours of history, about ten days. The models are chosen and reported on one holdout each, which Kuhn and Johnson (section 4.3) describe as a weak basis, so the holdout errors are probably a little optimistic; the models that pass the residual test are within 1% of each other, so the choice among them rests on the preference for the simplest model and not on a measured difference.

C7. Techniques used and not used

Used: time plots and the ACF to look at the hourly series; the KPSS test and the number of differences for stationarity; the Ljung-Box test and the ACF to judge white noise; the strength of a trend and of a 24-hour pattern; the four benchmarks, ETS, and automatic and manual ARIMA models; holdout RMSE and MAE, with the simplest model within a 1% tolerance of the best chosen (Kuhn and Johnson, section 4.6); the mean of the readings within each hour; and the average of the neighboring hours to fill the four missing hours.

Not used: a Box-Cox transformation, because the spread of both series looks constant over time (section C4). A decomposition for forecasting, because the 24-hour pattern is weak. Rolling cross-validation, to keep the comparison simple. KNN or model-based interpolation for the four missing hours, because they are isolated and the series has no pattern to follow.

Forecast files

write_xlsx(list(
  ATM_May2010 = may_forecast |>
    transmute(DATE, ATM, forecast_cash_hundreds = forecast, forecast_cash_dollars = forecast_dollars,
              lower95_hundreds = lower95, upper95_hundreds = upper95, model),
  Power_2014 = power_forecast |> mutate(month = as.character(month)),
  Water_Pipe1_Week = pipe_forecast |> filter(pipe == "pipe1") |>
    transmute(hour = format(hour, "%Y-%m-%d %H:%M"), forecast_flow, lower95, upper95, model),
  Water_Pipe2_Week = pipe_forecast |> filter(pipe == "pipe2") |>
    transmute(hour = format(hour, "%Y-%m-%d %H:%M"), forecast_flow, lower95, upper95, model)
), "Project1_Forecasts.xlsx")

The forecasts are written to Project1_Forecasts.xlsx: ATM_May2010 (daily, per ATM, in hundreds and in dollars), Power_2014 (monthly, in KWH), and for the bonus Water_Pipe1_Week and Water_Pipe2_Week (hourly, 168 hours each). The empty interval cells for ATM3 are intentional (section A8).