An initial look at a summary of the data in the datas set shows that there are four separate ATMs.
It shows an equal number of entries per ATM but also indicates some rows with missing data. The first task is to decide how to handle this missing data.
# A tibble: 5 × 6
atm sum mean max min n
<chr> <dbl> <dbl> <dbl> <dbl> <int>
1 ATM1 NA NA NA NA 365
2 ATM2 NA NA NA NA 365
3 ATM3 263 0.721 96 0 365
4 ATM4 173026. 474. 10920. 1.56 365
5 <NA> NA NA NA NA 14
Inspecting for NA values.
atm_raw |>summarize_all(~sum(is.na(.)))
# A tibble: 1 × 3
date atm cash_in_hundreds
<int> <int> <int>
1 0 14 19
# A tibble: 19 × 3
date atm cash_in_hundreds
<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
6 2010-05-01 <NA> NA
7 2010-05-02 <NA> NA
8 2010-05-03 <NA> NA
9 2010-05-04 <NA> NA
10 2010-05-05 <NA> NA
11 2010-05-06 <NA> NA
12 2010-05-07 <NA> NA
13 2010-05-08 <NA> NA
14 2010-05-09 <NA> NA
15 2010-05-10 <NA> NA
16 2010-05-11 <NA> NA
17 2010-05-12 <NA> NA
18 2010-05-13 <NA> NA
19 2010-05-14 <NA> NA
There are 19 rows that contain NA values of 1,474 rows. This is 1.28% of the total data set.
There are two types of missing data.
1. There are 14 rows with a date but no atm and no cash value. As the goal is to forecast May of 2010, these are likely forecasting placeholders, for the time being they will be dropped from the dataset during analysis.
There are 5 values that have an ATM but no cash value. These can be imputed based on the other values in the data set for the specific ATMs’ missing data. For this small number of missing data points, a potential approach is to use the median value from the same ATM and the same day of the week.
# A tibble: 4 × 7
atm n sum mean max min median
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 365 30654 84.0 180 1 91
2 ATM2 365 22750. 62.3 147 0 66
3 ATM3 365 263 0.721 96 0 0
4 ATM4 365 173026. 474. 10920. 1.56 404.
The summary data provides an interesting observation. The max value of ATM4 is over $1 million dollars. When all the data is plotted there are two immediate takeaways: - ATM4 has the highest daily total cash being withdrawn - The withdrawal day of $1,091,976 dollars is an outlier.
atm |>ggplot() +geom_point(aes(x = date, y = cash_imputed, color = atm))
That single value from ATM4 is so extreme that it has the potential to skew summarized data and forecasts. To proceed with the analysis, removing this outlier from the data set should be considered.
With the goal of forecasting how much will be taken out of each ATM each day for May 2010, this single outlier will be treated as an isolated anomaly and therefore will be removed from the dataset, prior to analysis. It will be replaced in the same manner as the missing data was earlier, with the median for the day of the week for that specific ATM.
# A tibble: 4 × 7
atm n sum mean max min median
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 365 30654 84.0 180 1 91
2 ATM2 365 22750. 62.3 147 0 66
3 ATM3 365 263 0.721 96 0 0
4 ATM4 365 162407. 445. 1712. 1.56 403.
With the outlier removed, the highest value for ATM 4 is now $171207.50.
atm |>ggplot() +geom_boxplot(aes(x = atm, y = cash_imputed, color = atm))
The box plot shows the notable difference in the level of the withdrawals by ATM.
atm |>ggplot() +geom_line(aes(x = date, y = cash_imputed, color = atm))
The line plot shows more of the distribution, however they are difficult to separate out as they each exist at a different level. Below the line plots are separated by ATM.
atm |>ggplot() +geom_line(aes(x = date, y = cash_imputed, color = atm)) +facet_grid(atm ~ ., scales ="free_y")
The line plots provide a new piece of insight. Almost all of ATM3’s values are $0. This indicates several possible scenarios for ATM3:
It is a new ATM and has just started being used.
It has been broken, empty or non-operational until recently.
The shape of the data for ATM3 will be relevant to the approach for forecasting.
Looking at the shape of the other three plots and they all have similarities:
There does not seem to be a trend for any of the ATMs.
There withdrawals all seem to be around the same level per ATM.
Each atm has a consistent level of variation in itself.
It is possible that there is a seasonal variation based on day of the week.
Forecast Experimentation
Converting the data to a tsibble and creating separate tsibbles for each atm.
Considering a Box-Cox transformation: Since the variation does not appear to change based on level of the data, a Box-Cox transformation is not necessary for this data.
Decomposition supports that there is no sustained upward or downward trend but the underlying level does change over time.
The weekly seasonality appears to have a low withdrawal day one day per week.
Decomposition notes for ATM3:
Consistent with the earlier assessment that ATM3 was not being withdrawn from until late in April of 2010.
Forecasting
The following code fits ATM1, ATM2, and ATM4 with several forecasting models through March 31st 2010 as the training data. The models are then tested on the withheld month of April in order to discover which model has the lowest root mean square error for each model.
The models being trained include:
seasonal naive
exponential smoothing models with seasonality in additive, multiplicative, and an auto calculated model.
ARIMA models with seasonality of a 7 day period in 001113, 102111, and an auto calculated model.
The root mean square error (RMSE) data for ATM1 indicate that the exponential smoothing additive model with seasonality has the lowest RMSE (the auto-calculated model has identical settings).
The residual plot indicates the errors are centered around zero and look largely like noise. The acf has no determinate shape, although a couple positions exceed the boundary limits. The histogram of residuals appears normal.
With a model selected, the full data, not just the training set is fit below, and the month of MAY 2010 is forecast.
The 95% prediction interval is extrapolated. Then two columns are created:
“atm_total_cash_to_stock” is the upper limit of the 95% prediction interval, or the 97.5th percentile rounded up and is a forecast of what each ATM should be stocked with, in hundreds of dollars, for the day, to attempt to ensure there is adequate cash to cover the days withdrawals.
“atm_forecasted_cash_withdrawal_per_day” is the forecast of what will be taken out each day in hundreds of dollars. It is the rounded up mean forecast by the model.
The root mean square error (RMSE) data for ATM2 indicate that the exponential smoothing additive model with seasonality has the lowest RMSE (the auto-calculated model has identical settings).
The residual plot indicates the errors are centered around zero and look largely like noise. The acf has no determinate shape, although a couple positions exceed the boundary limits. The histogram of residuals appears normal.
With a model selected, the full data, not just the training set is fit below, and the month of MAY 2010 is forecasted.
The 95% prediction interval is extrapolated. Then two columns are created:
“atm_total_cash_to_stock” is the upper limit of the 95% prediction interval, or the 97.5th percentile rounded up and is a forecast of what each ATM should be stocked with, in hundreds of dollars, for the day, to attempt to ensure there is adequate cash to cover the days withdrawals.
“atm_forecasted_cash_withdrawal_per_day” is the forecast of what will be taken out each day in hundreds of dollars. It is the rounded up mean forecast by the model.
Generally the same workflow has been used to analyze ATM1, ATM2, and ATM4, but ATM3 must be considered separately. As discussed earlier, it mostly contains no transactions until recently the last three days of April.
autoplot(atm3, cash_imputed)
With only three points with data above zero, and those being the three most recent, it is harder to determine a forecast. As a result of this small number of non-zero points, a naive forecast approach is the most reasoned choice. This method assumes that the non-zero values represent a new operational state and that going forward, the atm will be actively aused.
The residual plot indicates the errors are centered around zero and look largely like noise. The acf has no determinate shape and appears to be noise, and all points are within the boundary lines. The histogram of residuals approaches is right skewed.
With a model selected, the full data, not just the training set is fit below, and the month of MAY 2010 is forecast.
The 95% prediction interval is extrapolated. Then two columns are created:
“atm_total_cash_to_stock” is the upper limit of the 95% prediction interval, or the 97.5th percentile rounded up and is a forecast of what each ATM should be stocked with, in hundreds of dollars, for the day, to attempt to ensure there is adequate cash to cover the days withdrawals.
“atm_forecasted_cash_withdrawal_per_day” is the forecast of what will be taken out each day in hundr3eds of dollars. It is the rounded up mean forecast by the model.
power_raw |>autoplot(kwh) +labs(title ="Power usage in Kilowatt-Hours",x ="Year, Month",y ="Kilowatt-Hours")
The plot of power usage shows several things: - There is the possibility of a slight upward trend, particularly starting around 2009. - The data has fluctuations which appear to be seasonality. - There is missing data around 2008. - There is one very low value in 2010 which appears to be an outlier.
Handling the anomalies
To start, we will examine the outlier.
power_raw |>filter(date >=yearmonth("2010 Jan") & date <=yearmonth("2010 Dec"))
# A tsibble: 12 x 3 [1M]
case_sequence date kwh
<dbl> <mth> <dbl>
1 877 2010 Jan 9397357
2 878 2010 Feb 8390677
3 879 2010 Mar 7347915
4 880 2010 Apr 5776131
5 881 2010 May 4919289
6 882 2010 Jun 6696292
7 883 2010 Jul 770523
8 884 2010 Aug 7922701
9 885 2010 Sep 7819472
10 886 2010 Oct 5875917
11 887 2010 Nov 4800733
12 888 2010 Dec 6152583
The outlier value is 2010 July at 770523. This is a tenth of the other values nearby in the timeline. This strongly indicates that this is an anomalous value. Without more information, it is difficult to determine why this one value is so much lower than all other values. It is potentially an error in either collection of data or recording of data. While it is theoretically possible that an unseasonably cool July in 2010 meant less power usage, considering how much lower it is than all other data in the data set, the likelihood that the value is accurate is low. Therefore, it is assumed that 2010 July is an isolated anomaly.
In this instance the median value of prior July energy usage will be used to impute the 2010 July value.
This code calculates the median of July energy usage prior to 2010 July. The dataframe will be updated later on.
Warning in report.mdl_df(power_fits): Model reporting is only supported for
individual models, so a glance will be shown. To see the report for a specific
model, use `select()` and `filter()` to identify a single model.
# A tibble: 8 × 11
.model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 snaive 6.31e+11 NA NA NA NA NA NA NA
2 smooth_addit… 3.73e+11 -2858. 5746. 5749. 5794. 3.44e11 3.76e11 4.54e+5
3 smooth_multi… 8.79e- 3 -2853. 5735. 5738. 5783. 3.55e11 3.84e11 7.21e-2
4 smooth_auto 8.68e- 3 -2851. 5732. 5735. 5780. 3.48e11 3.68e11 6.97e-2
5 arima_auto 3.28e+11 -2467. 4950. 4951. 4975. NA NA NA
6 arima_auto_b… 1.43e- 5 708. -1400. -1399. -1375. NA NA NA
7 arima301113 3.36e+11 -2468. 4955. 4956. 4986. NA NA NA
8 arima102111 3.51e+11 -2473. 4959. 4960. 4981. NA NA NA
# ℹ 2 more variables: ar_roots <list>, ma_roots <list>
The table of the RMSE is below. The lowest RMSE value is for the arima_auto model: ARIM(3,0,1)(2,1,0)[12] with drift.
power_test_forecasts |>filter(.model =="arima_auto") |>autoplot(power_ts, level =NULL) +labs(title ="Power usage: actual and forecast ARIMA(3,0,1)(2,1,0) with drift",x ="Year Month", y ="kWh",colour ="Model" ) +theme_minimal()
A look at the plots of residuals looks like the selected model is performing well. The residuals look like noise adn are centered around zero. The ACF has no specific patter and all values are within the boundary range. The histogram of residuals is approximately normal.
The intended approach to this task is slightly unclear. I see two possible approaches combining and averaging the two sets of pipe data together, or analyzing them separately. Both are attmepted here.
Combine both pipe data sets and forecast
The following code imports the data, converts the date time column to a date time data type and combines the two dataframes.
The data does not appear to be stationary. While it does not appear to have a trend or seasonal variation, it has significant non-seasonal variation and a notable change in variation. Around November 1, the level and the variation change significantly.
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).
Differencing centers the data around zero but does little to stabilize the variation.
water_lambda <- avg_pipe |>features(avg_hourly_water_flow, features = guerrero) |>pull(lambda_guerrero)avg_pipe |>autoplot(box_cox(avg_hourly_water_flow, water_lambda))
Box-Cox transformation is also not beneficial in stabilizing the variation.
Ultimately it may not be beneficial to transform this data to make it stationary.
A better approach may be to consider only data since November 1st as that data, when considered without the earlier data is stationary. From that we can will create a forecast.
# A mable: 1 x 1
arima_auto
<model>
1 <ARIMA(0,0,0)(0,0,1)[24] w/ mean>
Utilizing the auto ARIMA function, a lag of 24 is selected as a seasonal term. Interestingly the auto arima does not select a differencing coefficient.
The ARIMA(0, 0, 0)(0, 0, 1)[24] model is plotted above. Forecasts vary during the first 24 hours and then settles on the estimated mean. It is the influence of Pipe 1’s shorter data set that creates the earlier variation change in the combined data set. Ultimately, an analysis of the two data sets separately seems the preferred method. This is performed below.
The following code imports the data, converts the date time column to a date time data type and combines the two dataframes while retaining an id column.
avg_pipe |>autoplot(avg_hourly_water_flow) +facet_grid(pipe ~ .) +labs(title ="Hourly average water flow by pipe",x ="Date",y ="Average water flow" ) +guides(colour ="none") +theme_minimal()
The data is possibly stationary. The two sets of data have very different time spans and size of variation, which explains the variation difference in the original analysis of this data. However, as was learned in the earlier analysis, it is possible that there is a 24 hour seasonality.
The results of the the comparison of the RMSE indicates that for Pipe 1 a simple mean model have the lowest RMSE, while for Pipe 2 the exponential smoothing model (A, N, N) minimizes the RMSE. Although, in both instances, these are marginally improved over the next closest model. Additionally, the smoothing model for Pipe 2 is not notably different than the mean model used for Pipe 1.
pipe_final_forecast |>autoplot(avg_pipe_imputed, level =95) +facet_wrap(~ pipe, ncol =1) +labs(title ="Forecast for one week water flow",x ="Date",y ="Hourly average water flow" ) +theme_minimal()
The plot of the final forecasts show that both are essentially a mean forecast moving forward for the week, while the prediction intervals cover much of the variation within each model.
