Data 624 Project 1
Part A – ATM Forecast
ATM624Data.xlsx
In part A, we will forecast how much cash is taken out of 4 different ATM machines for May 2010. The data is given in a single file. The variable ‘Cash’ is provided in hundreds of dollars, other than that it is straight forward.
ATMdataRaw<-read_excel("data/ATM624Data.xlsx") %>% mutate(DATE = excel_numeric_to_date(DATE)) %>%
filter(DATE<as_date("2010/05/01"))%>%filter(is.na(Cash)==F)
head(ATMdataRaw) %>%kable() %>%
kable_styling(
full_width = F, position="center", bootstrap_options = c("hover")) %>%
scroll_box() %>%
kable_classic_2()| DATE | ATM | Cash |
|---|---|---|
| 2009-05-01 | ATM1 | 96 |
| 2009-05-01 | ATM2 | 107 |
| 2009-05-02 | ATM1 | 82 |
| 2009-05-02 | ATM2 | 89 |
| 2009-05-03 | ATM1 | 85 |
| 2009-05-03 | ATM2 | 90 |
ATMdata <-ATMdataRaw %>%
pivot_wider(names_from = ATM, values_from = Cash) %>%
remove_empty(which = "cols")
range(ATMdata$DATE)[1] "2009-05-01" "2010-04-30"The format of the data is shown above. This dataframe is provided in a “long” format and provides the DATE, the ATM, and how much Cash was withdrawn from each ATM. The range of this dataset is also computed showing withdrawals from May 2009 to April of 2010 - A 1 year sample with a sample frequency of 1-day.
Using the plot_time_series function from timetk, We can plot each of the time series data and get a general gist of our time series. At first glance, we can see that ATM1 and ATM2 seem like fairly stable time series while ATM3 only has a few points that will likely make modeling this series uneccessary and predictions should probably simply be a mean of the data. ATM4 has an extreme outlier that will need to be adjusted for in our dataset.
interactive = F
ATMdataRaw %>%
plot_time_series(DATE,Cash,
.interactive=interactive,
.color_var = ATM,
.facet_vars = ATM,
.facet_ncol = 2,
.smooth = F)
#ATM_ts<-ts(select(ATMdata,-DATE), start = c(2009,05,01), frequency = 1)Looking at the lag diagnostics of all of our features show that there are some autocorrelations experienced in ATM1 and ATM2. We expect that the outlier in ATM4 messes with the ACF and PACF making it differ. again, ATM3 shows ACFs that are virtually 0 on its lags due to no variance in the data with the exception of the final few points. likely an ATM that recently came on-line and only has a few withdrawals recorded.
ATMdataRaw %>% group_by(ATM) %>%
plot_acf_diagnostics(DATE,Cash,
.interactive=interactive,
.show_white_noise_bars = T,
.lags = 25)
Models to be Evaluated
Along these datasets, we will attempt various models and see which performs best on our time series data. The models we will attempt are:
Model 1: ARIMA - Known as Auto Regressive Integrated Moving Average which is a class of models that explains a time series based on its own past values (lagged values), and uses these lags (both seasonal and non seasonal) to produce predictions. ARIMA models are very flexible and consider differencing, auto-regression and moving averages in its forecasting computation.
Model 2: Boosted ARIMA - This is an ARIMA model that also uses boosting to improve model errors on exogenous regressors, which can be engineered time features or other external features of dataset that are not explicitly time.
Model 3: ETS - Known as Exponential Smoothing based on a 3-character string that denotes the type of model being done e.g. [A,N,N].
R Documentationdefines this clearly stating “The first letter denotes the error type; the second letter denotes the trend type; and the third letter denotes the season type. In all cases,”N“=none,”A“=additive,”M“=multiplicative and”Z“=automatically selected. So, for example,”ANN" is simple exponential smoothing with additive errors, “MAM” is multiplicative Holt-Winters’ method with multiplicative errors, and so on."Model 4: Prophet - An open source tool developed by facebook used for making forecasts on univariate time series data. These models are additive regression models with non-linear trend fitting based on seasonality and holiday effects. They are robust to missing data, trend shifts and outlier handling.
Model 5: Naive - Naive forecasting is a simpler forecasting method that simply uses the last known value as the forecast for future predictions
#ATM1<-ATM_ts[,1]
#ATM2<-ATM_ts[,2]
#ATM3<-ATM_ts[,3]
#ATM4<-ATM_ts[,4]
ATM1<-ATMdataRaw %>% filter(ATM=="ATM1")
ATM2<-ATMdataRaw %>% filter(ATM=="ATM2")
ATM3<-ATMdataRaw %>% filter(ATM=="ATM3")
ATM4<-ATMdataRaw %>% filter(ATM=="ATM4")
#ggtsdisplay((ATM1))
#ATM 1
Exploration
A seasonal diagnostic of the data was performed to see if there are any trends that could potentially weigh in on the forecasts. the most notable ones seem to on the weeklies where Thursdays tend to have a lower median cash withdrawal than other days.
dataset<-ATM1
dataset %>%
plot_seasonal_diagnostics(DATE,Cash,
.interactive = interactive)
Additionally, STL decomposition can show the breakdown of the dataset by trend and seasonality (if applicable). In this case, we notice that the overall trend starts to decline but flatens out and remains stagnant.
dataset %>%
plot_stl_diagnostics(DATE,Cash,
.interactive = interactive,
.frequency = "auto",.trend = "auto",
.feature_set = c("observed", "season", "trend"))
Lastly, Anomoly diagnostics allow us to assess where there are potential outliers in this dataset. The ATM1 dataset does appear to have outliers but none that are expected to adversely effect modelling. nonetheless, the function step_clean_ts is used on all models to clean for outliers as well as linearly impute missing data.
dataset %>%
plot_anomaly_diagnostics(DATE,Cash,
.interactive = interactive,
.ribbon_alpha = 0,
.anom_size = 3)
Modeling
The dataset for ATM1 is split into a train/testing dataset of 90/10 as shown below.
dataset<-dataset %>% rename("date"="DATE", "value"="Cash")
set.seed(1)
splits <- initial_time_split(dataset, prop = 0.9)
bind_rows(
training(splits) %>% mutate(set = "train"),
testing(splits) %>% mutate(set = "test")
) %>% ggplot()+
geom_line(mapping = aes(x = date, y = value, color = set))+
defaulttheme
modelexport<-function(dataset, RENAME = T){
if(RENAME == T){
dataset<-dataset %>% rename("date"="DATE", "value"="Cash")
}
set.seed(1)
splits <- initial_time_split(dataset, prop = 0.9)
default_recipe<-recipe(value ~ date, data = training(splits)) %>%
step_ts_clean(value)
model_fit_arima <- arima_reg() %>%
set_engine(engine = "auto_arima")# %>%
wflw_fit_arima <- workflow() %>%
add_recipe(default_recipe) %>%
add_model(model_fit_arima) %>%
fit(training(splits))
model_fit_arima_boosted <- arima_boost(
min_n = 2,
learn_rate = 0.015,
mtry =3) %>%
set_engine(engine = "auto_arima_xgboost")
recipe2<- recipe(value ~ date, data = training(splits)) %>%
step_timeseries_signature(date)
wflw_fit_arima_boosted <- workflow() %>%
add_recipe(recipe2) %>%
add_model(model_fit_arima_boosted) %>%
fit(training(splits))
# Model 3: ets ----
model_fit_ets <- exp_smoothing() %>%
set_engine(engine = "ets")
wflw_fit_ets <- workflow() %>%
add_recipe(default_recipe) %>%
add_model(model_fit_ets) %>%
fit(training(splits))
# Model 4: prophet ----
model_fit_prophet <- prophet_reg() %>%
set_engine(engine = "prophet")
wflw_fit_prophet <- workflow() %>%
add_recipe(default_recipe) %>%
add_model(model_fit_prophet) %>%
fit(training(splits))
# Model 5: snaive ----
model_fit_naive <- naive_reg() %>%
set_engine(engine = "naive")
wflw_fit_naive <- workflow() %>%
add_recipe(default_recipe) %>%
add_model(model_fit_naive) %>%
fit(training(splits))
models_tbl <- modeltime_table(
wflw_fit_arima,
wflw_fit_arima_boosted,
wflw_fit_ets,
wflw_fit_prophet,
wflw_fit_naive
)
calibration_tbl <- models_tbl %>%
modeltime_calibrate(new_data = testing(splits))
calibration_tbl
}
calibration_tbl<-modelexport(ATM1)The Five different models mentioned above were created and tested on the data as shown in the forecast and the figures are shown below. the Naive model clearly does the poorest based on the trend evaluation. Prophet also does not seem to capture the magnitude of the frequencies as well either. Based on the table blow. The ARIMA models do the best at predicting with the lowest RMSE (root mean square error), and highest rsq (r-squared - coefficient of determination). The Boosted Arima provides very marginal improvement over the ARIMA, or ETS models. The preferred model of choice would be the ARIMA or ETS.
for (i in 1:nrow(calibration_tbl)){
assign(paste0("p",i),
calibration_tbl[i,]%>%
modeltime_forecast(
new_data = testing(splits),
actual_data = dataset
) %>%
plot_modeltime_forecast(
.legend_max_width = 25, # For mobile screens
.interactive = F,
.title = NULL,
.color_lab = NULL,
.conf_interval_show = F
)
)
}
gridExtra::grid.arrange(p1,p2,p3,p4,p5, ncol = 2)
calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = F,
.title = "ATM1 Model Accuracies"
)| ATM1 Model Accuracies | ||||||||
|---|---|---|---|---|---|---|---|---|
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
| 1 | ARIMA(0,0,1)(0,1,1)[7] | Test | 9.45 | 32.56 | 0.27 | 21.53 | 12.60 | 0.85 |
| 2 | ARIMA(0,0,1)(0,1,1)[7] W/ XGBOOST ERRORS | Test | 9.43 | 31.96 | 0.27 | 21.20 | 12.59 | 0.85 |
| 3 | ETS(A,N,A) | Test | 9.59 | 32.58 | 0.28 | 21.65 | 12.71 | 0.85 |
| 4 | PROPHET | Test | 31.97 | 320.01 | 0.92 | 54.42 | 43.43 | 0.01 |
| 5 | NAIVE | Test | 20.49 | 277.41 | 0.59 | 34.42 | 32.12 | NA |
ATM 2
Exploration
Similarly for ATM2, we wee that there are no noteable month to month, or quarterly trends within the dataset. However, like ATM1, there is a significant dip in withdrawals on Wednesdays and Thursdays. This may imply a frequency of 7 for a weekly seasonality based trend.
dataset<-ATM2
dataset %>%
plot_seasonal_diagnostics(DATE,Cash,
.interactive = interactive)
The diagnostic plot shows that the overall trend of withdrawals at this ATM is continuously decreasing while the seasonality does not vary much.
dataset %>%
plot_stl_diagnostics(DATE,Cash,
.interactive = interactive,
.frequency = "auto",.trend = "auto",
.feature_set = c("observed", "season", "trend"))
This ATM proves to be the most robust in terms of data completion and lack of anomalies. only 1 anomaly is detected shown by the red dot below and this is not expected to have any impact on the forecast predictions of any of the models.
dataset %>%
plot_anomaly_diagnostics(DATE,Cash,
.interactive = interactive,
.ribbon_alpha = 0,
.anom_size = 3)
Modeling
The following figure represents that 90/10 test/train split for the model evalation
dataset<-dataset %>% rename("date"="DATE", "value"="Cash")
set.seed(1)
splits <- initial_time_split(dataset, prop = 0.9)
bind_rows(
training(splits) %>% mutate(set = "train"),
testing(splits) %>% mutate(set = "test")
) %>% ggplot()+
geom_line(mapping = aes(x = date, y = value, color = set))+
defaulttheme
Passing our ATM2 data into the function created to model ATM1 yields the results shown below. Unlike ATM1, ATM2 obtains a more significant boost in performance when using the Boosted ARIMA model in lieu of the ARIMA model. ETS also performs better than the ARIMA model. both Prophet and NAIVE models do not perform well. The model of choice for ATM2 would be the Boosted ARIMA model.
calibration_tbl<-modelexport(ATM2)
for (i in 1:nrow(calibration_tbl)){
assign(paste0("p",i),
calibration_tbl[i,]%>%
modeltime_forecast(
new_data = testing(splits),
actual_data = dataset
) %>%
plot_modeltime_forecast(
.legend_max_width = 25, # For mobile screens
.interactive = F,
.title = NULL,
.color_lab = NULL,
.conf_interval_show = F
)
)
}
gridExtra::grid.arrange(p1,p2,p3,p4,p5, ncol = 2)
calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = F,
.title = "ATM3 Model Accuracies"
)| ATM3 Model Accuracies | ||||||||
|---|---|---|---|---|---|---|---|---|
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
| 1 | ARIMA(2,0,2)(0,1,2)[7] | Test | 16.56 | Inf | 0.44 | 59.08 | 24.09 | 0.65 |
| 2 | ARIMA(3,0,3)(0,1,1)[7] W/ XGBOOST ERRORS | Test | 13.92 | Inf | 0.37 | 48.08 | 20.43 | 0.73 |
| 3 | ETS(A,N,A) | Test | 15.42 | Inf | 0.41 | 44.21 | 21.77 | 0.70 |
| 4 | PROPHET | Test | 46.88 | Inf | 1.26 | 104.91 | 55.92 | 0.24 |
| 5 | NAIVE | Test | 33.92 | Inf | 0.91 | 65.80 | 44.79 | NA |
ATM 3
the following figures below show the timeseries data for ATM3. This ATM only has about 3 points that are not 0. It is likely that this ATM was only recently commissioned at the end April 2010, or that it was out of service for the periods prior to that.
plot_time_series(ATM3, DATE, Cash,
.smooth = F,
.interactive = F)
ATM3 %>% filter(Cash>0) %>%
plot_time_series(DATE,Cash,.smooth = F, .interactive = F)
Due to the lack of information associated with ATM3, it is inappropriate to spend any appreciable effort in forecasting any advanced models. the most appropriate forecast for this ATM is likely either NAIVE (using the last known withdrawal), or simply using the mean of the withdrawals, in this case $87.7.
ATM3 %>% filter(Cash>0) %>% group_by(ATM) %>%
summarise(mean = mean(Cash)) # A tibble: 1 x 2
ATM mean
<chr> <dbl>
1 ATM3 87.7ATM 4
Exploration
ATM4 shows similar trends to ATM1 and ATM2, however there is a major outlier showing a withdrawal greater than 9000. this is likely a false value as this exceeds typical ATM withdrawal limits and can be removed from this dataset.
dataset<-ATM4
dataset %>%
plot_anomaly_diagnostics(DATE,Cash,
.interactive = interactive,
.ribbon_alpha = 0,
.anom_size = 3)
The STL decomposition shows a potentially weekly seasonality as well as an overall decreasing trend which is consistent with the other ATMs.
dataset<-ATM4 %>% filter(Cash<9000)
dataset %>%
plot_stl_diagnostics(DATE,Cash,
.interactive = interactive,
.frequency = "auto",.trend = "auto",
.feature_set = c("observed", "season", "trend"))
Similarly, after removing the outlier, we notice that Thursdays seem to be a day where ATMs are used less frequently than other days.
dataset %>%
plot_seasonal_diagnostics(DATE,Cash,
.interactive = interactive)
The following represents the 90/10, test/train split for ATM4.
Modeling
dataset<-dataset %>% rename("date"="DATE", "value"="Cash")
set.seed(1)
splits <- initial_time_split(dataset, prop = 0.9)
bind_rows(
training(splits) %>% mutate(set = "train"),
testing(splits) %>% mutate(set = "test")
) %>% ggplot()+
geom_line(mapping = aes(x = date, y = value, color = set))+
defaulttheme
ATM4 is very similar to ATM1 where the ARIMA models perform best followed very closely by the ETS. Boosted ARIMA provided no significant benefit over ARIMA and thus ARIMA or ETS would be the models of choice for ATM4.
calibration_tbl<-modelexport(ATM4)
for (i in 1:nrow(calibration_tbl)){
assign(paste0("p",i),
calibration_tbl[i,]%>%
modeltime_forecast(
new_data = testing(splits),
actual_data = dataset
) %>%
plot_modeltime_forecast(
.legend_max_width = 25, # For mobile screens
.interactive = F,
.title = NULL,
.color_lab = NULL,
.conf_interval_show = F
)
)
}
gridExtra::grid.arrange(p1,p2,p3,p4,p5, ncol = 2)
calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = F,
.title = "ATM4 Model Accuracies"
)| ATM4 Model Accuracies | ||||||||
|---|---|---|---|---|---|---|---|---|
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
| 1 | ARIMA(0,0,0) WITH NON-ZERO MEAN | Test | 227.15 | 546.16 | 0.71 | 62.39 | 273.80 | NA |
| 2 | ARIMA(0,0,0) WITH NON-ZERO MEAN W/ XGBOOST ERRORS | Test | 235.23 | 1038.67 | 0.68 | 67.01 | 281.28 | 0.07 |
| 3 | ETS(M,N,A) | Test | 303.46 | 751.80 | 0.95 | 73.15 | 395.43 | 0.11 |
| 4 | PROPHET | Test | 285.39 | 737.99 | 0.89 | 70.18 | 362.49 | 0.04 |
| 5 | NAIVE | Test | 449.84 | 1038.28 | 1.40 | 81.38 | 523.53 | NA |
Part B – Forecasting Power
ResidentialCustomerForecastLoad-624.xlsx
Part B consists of a simple dataset of residential power usage for January 1998 until December 2013. the goal is to model these data and a monthly forecast for 2014.
res_data<-read_excel("Data/ResidentialCustomerForecastLoad-624.xlsx") %>%
mutate(`YYYY-MMM`= ym(`YYYY-MMM`)) %>%
rename(date = `YYYY-MMM`, value = KWH) %>%
filter(is.na(value)==F)
head(res_data)# A tibble: 6 x 3
CaseSequence date value
<dbl> <date> <dbl>
1 733 1998-01-01 6862583
2 734 1998-02-01 5838198
3 735 1998-03-01 5420658
4 736 1998-04-01 5010364
5 737 1998-05-01 4665377
6 738 1998-06-01 6467147Exploration
We’ll start out by looking at the dataset overall for its behavior. From the time series plot below, we see that there is a very slight increasing trend with some cyclicality. We also notice that there is an increase in the fluctuations over time and may need to perform transformations in order to assess some of the models.
res_data %>%
plot_time_series(date,value,
.interactive=interactive,
.facet_ncol = 2,
.smooth = F)
This dataset is a monthly dataset, however, we will look at seasonality in multiple ways just to cover all bases. From the seasonal diagnostic plot below, it is clear that the dataset has a very strong monthly seasonal component to it that can and should be incorporated into the model
res_data %>% plot_seasonal_diagnostics(date,value,
.interactive = F)
Looking at the decomposition plot allows us to isolate different portions of the time series. The trend component shows a continous increase (not clearly shown in the raw observed plot). and the remainder highlights the increase in frequency of our dataset making a stronger case for a transformation.
res_data %>%
plot_stl_diagnostics(date,value,
.interactive = F,
.feature_set = c("observed","season","trend","remainder")
)
The data was split into a 90/10 test/train split and evaluated amongst the model types used in Part A. This includes ARIMA, Boosted ARIMA, ETS, Prophet, and NAIVE. The results presented below are based on the non-transformed residential customer data and show that our optimal model based on the Prophet forecast developed by Facebook. Second to this is the Boosted ARIMA model with a drift component however this is not much better than the regular ARIMA model with a drift component. The pdq PQD Values for the ARIMA models are (0,0,1) (2,1,1) which implies one order of non-seasonal moving average, 2 orders of seasonal autoregression, 1 order of seasonal differencing and 1 order of seasoanl moving average.
dataset<-res_data
set.seed(1)
splits <- initial_time_split(dataset, prop = 0.9)
calibration_tbl<-modelexport(res_data, RENAME = F)
for (i in 1:nrow(calibration_tbl)){
assign(paste0("p",i),
calibration_tbl[i,]%>%
modeltime_forecast(
new_data = testing(splits),
actual_data = res_data
) %>%
plot_modeltime_forecast(
.legend_max_width = 25, # For mobile screens
.interactive = F,
.title = NULL,
.color_lab = NULL,
.conf_interval_show = T
)
)
}
gridExtra::grid.arrange(p1,p2,p3,p4,p5, ncol = 2)
Similarly, applying the same models to a log-transformed dataset is shown below. The log-transformed dataset performed better than the untransformed dataset with R-squared values reaching 0.79 on the ETS and 0.78 with the ARIMA modeling, both out-performing Prophet. The plot below shows the forecasted results of the models.
calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = F,
.title = "Residential Customer Forecast Accuracies"
)| Residential Customer Forecast Accuracies | ||||||||
|---|---|---|---|---|---|---|---|---|
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
| 1 | ARIMA(0,0,1)(2,1,1)[12] WITH DRIFT | Test | 919056.4 | 11.37 | 0.62 | 12.50 | 1326028.2 | 0.42 |
| 2 | ARIMA(0,0,1)(2,1,0)[12] WITH DRIFT W/ XGBOOST ERRORS | Test | 902228.1 | 11.10 | 0.61 | 12.21 | 1289179.9 | 0.46 |
| 3 | ETS(M,N,A) | Test | 1233567.4 | 15.86 | 0.84 | 16.48 | 1584978.7 | 0.10 |
| 4 | PROPHET | Test | 658686.8 | 8.13 | 0.45 | 8.49 | 945626.4 | 0.68 |
| 5 | NAIVE | Test | 1885544.2 | 22.18 | 1.28 | 26.65 | 2426191.7 | NA |
res_datatrans<-
res_data %>% mutate(value = log(value))
calibration_tbl<-modelexport(res_datatrans, RENAME = F)
calibration_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = F,
.title = "Log Transformed Residential Customer Forecast Accuracies"
)| Log Transformed Residential Customer Forecast Accuracies | ||||||||
|---|---|---|---|---|---|---|---|---|
| .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
| 1 | ARIMA(0,0,1)(0,1,1)[12] WITH DRIFT | Test | 0.08 | 0.48 | 0.40 | 0.49 | 0.10 | 0.78 |
| 2 | ARIMA(0,0,0)(2,0,0)[12] WITH NON-ZERO MEAN W/ XGBOOST ERRORS | Test | 0.37 | 2.34 | 1.91 | 2.33 | 0.42 | 0.04 |
| 3 | ETS(A,N,A) | Test | 0.07 | 0.44 | 0.36 | 0.44 | 0.09 | 0.79 |
| 4 | PROPHET | Test | 0.10 | 0.63 | 0.51 | 0.63 | 0.12 | 0.74 |
| 5 | NAIVE | Test | 0.27 | 1.70 | 1.40 | 1.72 | 0.34 | NA |
for (i in 1:nrow(calibration_tbl)){
assign(paste0("p",i),
calibration_tbl[i,]%>%
modeltime_forecast(
new_data = testing(splits),
actual_data = res_datatrans
) %>%
plot_modeltime_forecast(
.legend_max_width = 25, # For mobile screens
.interactive = F,
.title = NULL,
.color_lab = NULL,
.conf_interval_show = T
)
)
}
gridExtra::grid.arrange(p1,p2,p3,p4,p5, ncol = 2)
Based on the results, we will use the ETS model to make our final predictions. The model will be re-fit to the entire series and then predicted. The excel file transforms the logged data backed to its original KWH format.
predictions<-
calibration_tbl %>%
# Remove ARIMA model with low accuracy
filter(.model_id == 3) %>%
# Refit and Forecast Forward
modeltime_refit(res_datatrans) %>%
modeltime_forecast(h = "12 months", actual_data = res_datatrans)
predictions%>%
plot_modeltime_forecast(.interactive = FALSE)
#predictions %>% filter(.key=="prediction") %>%
# select(.index,.value) %>%
# mutate(.value = 10^(.value)) %>%
# rename(date = .index, predicted_KWH = .value) %>%
# xlsx::write.xlsx("data/res_data_predictions.xlsx")