Data 624: Project One Part A
2022-03-21
pacman::p_load(tidyverse,DataExplorer,fpp3,tsibble,imputeTS,cowplot,xts)ATM_Data <- readxl::read_excel("Data/ATM624Data.xlsx",
col_types = c("date", "text", "numeric"))1 Introduction
For this assignment I have been tasked with
forecasting how much cash is taken out of 4 different ATM machines for
May 2010. The variable ‘Cash’ is provided in hundreds of dollars. I will
explore the data, looking for any irregular values or missing data. From
there we will fix and model the data and use those models to create
forecasts into the future.
I will also be displaying the power of the
purr map function. Due to the nature of this data, the Cash values for
each ATM need to be handled separately. This can be easily accomplished
by making a list of separate data frames split by ATM. From there, the
map function can be used to repeat the same code on each
dataframe.
2 Exploration
ATM_Data There are an irregular number of 0 values for
ATM 3. Not only that, but upon further inspection, the only values that
aren’t 0 match perfectly with another ATM. It seems apparent that none
of the cash values for ATM 3 can be trusted. I would advise that the
data be recollected or fixed before doing any modeling on it. For this
assignment I will drop it as a variable. The most I could do is use
other ATM values to fill in the blanks, but at that point those
predictions can hardly even be trusted to properly represent ATM
3.
ATM_Data %>% filter(Cash ==0)ATM_Data %>% filter(ATM=="ATM3") %>% count(Cash)ATM_Data %>%
filter(ATM %in% c("ATM1","ATM3") & DATE >= date("2010-04-28") ) The month we want to forecast for is found in
the data, full of NA values. It will be removed here and then re-added
back when forecasting later.
ATM_Data %>% filter(DATE >= as.Date( "2010-05-01") )ATM_Data <- ATM_Data %>%
filter(DATE < as.Date( "2010-05-01") & ATM !="ATM3" ) A small portion of the data for both ATM 1 and
2 appear to be missing.
ATM_Data %>% plot_missing()
ATM_Data %>% filter(!complete.cases(Cash)) A tsibble time series for the data as well as
a split version by ATM is created below. Splitting makes it so that each
ATM can be easily accessed separately when needed.
ATM_TS <- ATM_Data %>% mutate(DATE = as_date(DATE)) %>%
as_tsibble(index = DATE,key = ATM) By_ATM <- split(ATM_TS,ATM_TS$ATM) The initial time series plot is shown below.
It is clear that ATM 4 has way more cash taken out than the rest, with
one suspicious looking value that appears to be many times greater than
all other values.
autoplot(ATM_TS)
The patterns for each ATM can be seen clearly
by looking at each on their own scale. ATM 1 and 2 appear to be almost
identical in scale, and the patterns seem similar as well. Laying them
on top of each other further confirms this.
autoplot(ATM_TS) + facet_grid(ATM ~ ., scales = "free_y")
autoplot(bind_rows(By_ATM[1:2]))
The distribution of values for ATM 1 and 2
appear to be healthy. That one value for ATM 4 though is skewing it
entirely.
ggplot(ATM_TS,aes(ATM,Cash,group=ATM )) +
geom_boxplot(aes(fill=Cash)) +
facet_grid(ATM ~ ., scales = "free_y") 
3 Erroneous Value Handling
It seems that a single value of ATM 4 is more
than two times as large as the second largest value. The other values of
that month come no where even close to that value. Regardless if this is
an actual value or not, it does not give us any meaningful information,
as it would appear this only occurs one time and isn’t something we need
to concern our selves with for future predictions. With this in mind, we
chose to remove this value and treat it as missing. We will impute it
back in the following section.
By_ATM$ATM4 %>%
arrange(desc(Cash))By_ATM$ATM4 %>%
filter( month(DATE)==2) %>% autoplot()
ATM_TS <- ATM_TS %>% mutate(Cash=ifelse(Cash<2000,Cash,NA) )
By_ATM <- split(ATM_TS,ATM_TS$ATM) 4 Missing Data Handling
Exploring the missing values in more detail,
the bulk of them come from ATM 1 and 2 and appear to occur only in June
of 2009. The missing values exhibit a clear pattern.
bind_rows(By_ATM[1:2]) %>%
filter(month(DATE) == 6) %>% autoplot()
The surrounding months can be viewed to gain
a better understanding off the nearby values woithout missing
data.
bind_rows(By_ATM[1:2]) %>% mutate(Month = month(DATE) ) %>% filter( Month >=5 & Month <=8 ) %>%
autoplot() + facet_grid(Month~ ., scales = "free_y")
ts <- By_ATM %>%
map (
~xts( .$Cash ,order.by = .$DATE)
)The function below visualizes the distribution
of missing values within our data.
ts %>%
map2(c('/2009-8','/2009-8','2010'),
~ggplot_na_distribution(.x[.y] )
)
Comparing imputation techniques below, it
would appear that na_kalman makes more sense for this data as it more
accurately matches the existing values. It will be will used to impute
the missing values.
ts %>%
map2(c('2009','2009','2010'),
~ggplot_na_imputations(.x[.y],na_kalman(.x[.y]) )
)
ts %>%
map2(c('2009','2009','2010'),
~ggplot_na_imputations(.x[.y],na_interpolation(.x[.y]) )
)
complete_ts <- ATM_TS %>%group_by(ATM) %>% mutate(Cash=na_kalman(Cash) )%>% ungroup()
By_ATM <- split(complete_ts,complete_ts$ATM) All three of our series look natural with ATM 4
no longer skewed by that one outlier value.
autoplot(complete_ts,Cash) + facet_grid(ATM ~ ., scales = "free_y")
plot_boxplot(complete_ts,by="ATM")
5 Modeling
5.1 Data Handling
Iterating over the split data frames allows me
to view original data as well as the Box-Cox transformations for each
ATM. Doing so levels out the cash value for each ATM. Therefore it would
appear appropriate to apply this transformation for all our ATM
values.
By_ATM %>%
map2( c("ATM 1","ATM 2","ATM 4"),
~plot_grid(autoplot(.x, box_cox(Cash, features(.x, features = guerrero)$lambda_guerrero))+labs(y= "Lambda Cash"
,title = paste(.y ,"Lambda Transform"))
,autoplot(.),ncol=1 )
)
By_ATM <- By_ATM %>%
map( ~.x %>%
mutate(`Lam Cash`=box_cox(Cash, features(.,Cash, features = guerrero)$lambda_guerrero))
)5.2 Exponential Smoothing Models
5.2.1 Automatic Models
To start off with, the automatic ETS function is
used to find the best model for each ATM as a baseline.
auto_ETS <- By_ATM %>%
map( ~.x %>%
model( ETS(`Lam Cash`))
) The results of the automatic process below.
Both ATM 1 and 2 received the same model, ETS(A,N,A). ATM 3 got a
different model with ETS(M,N,A). These choice indicates that the
modeling algorithm did not find a noticeable trend for any of our ATM
values.
auto_ETS %>% map( ~.x %>%
pivot_longer(!ATM,names_to = "Model name",
values_to = "Orders")
)$ATM1
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM1 ETS(`Lam Cash`) <ETS(A,N,A)>
$ATM2
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM2 ETS(`Lam Cash`) <ETS(A,N,A)>
$ATM4
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM4 ETS(`Lam Cash`) <ETS(M,N,A)>5.2.2 Additional Modeling
A few few other model combinations are used for
comparison with the automatically generated ones.
compare_ETS <- By_ATM %>%
map( ~.x %>%
model(
ANN = ETS(`Lam Cash` ~ error("A") + trend("N") + season("N")),
AAN = ETS(`Lam Cash` ~ error("A") + trend("A") + season("N")),
MNN = ETS(`Lam Cash` ~ error("M") + trend("N") + season("N")),
Multiplicative = ETS(`Lam Cash` ~ error("M") + trend("A") +
season("M")),
`Multiplicative No Trend` = ETS(`Lam Cash` ~ error("M") + trend("N") +
season("M")),
`Damped Multiplicative` = ETS(`Lam Cash` ~ error("M") +
trend("Ad") + season("M"))
)
) compare_ETS %>% map2(auto_ETS,
~bind_rows(accuracy(.x),accuracy(.y)) %>%
arrange(RMSE)
)$ATM1
# A tibble: 7 x 11
ATM .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 Multiplicat~ Trai~ 0.0561 1.55 0.942 -9.59 19.1 0.821 0.739 0.109
2 ATM1 Multiplicat~ Trai~ 0.00993 1.57 0.969 -10.2 19.5 0.845 0.745 0.127
3 ATM1 Damped Mult~ Trai~ 0.0135 1.57 0.969 -10.1 19.5 0.845 0.746 0.131
4 ATM1 ETS(`Lam Ca~ Trai~ -0.0225 1.78 1.01 -Inf Inf 0.878 0.848 0.0934
5 ATM1 MNN Trai~ 0.0527 2.83 2.01 -26.3 43.7 1.75 1.35 0.0172
6 ATM1 ANN Trai~ -0.00196 2.92 2.06 -Inf Inf 1.79 1.39 0.0202
7 ATM1 AAN Trai~ -0.00115 2.95 2.10 -Inf Inf 1.83 1.40 0.0187
$ATM2
# A tibble: 7 x 11
ATM .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM2 ETS(`Lam Cash`) Trai~ -0.220 7.42 5.20 -Inf Inf 0.862 0.842 0.0310
2 ATM2 Multiplicative Trai~ 0.694 9.98 8.08 -99.1 146. 1.34 1.13 0.0538
3 ATM2 Damped Multipl~ Trai~ 0.754 10.1 8.29 -100. 146. 1.37 1.15 0.0278
4 ATM2 Multiplicative~ Trai~ -0.624 10.1 8.18 -119. 162. 1.36 1.15 0.0791
5 ATM2 MNN Trai~ 0.353 11.5 9.79 -126. 175. 1.62 1.31 0.0378
6 ATM2 AAN Trai~ 0.201 11.6 9.88 -Inf Inf 1.64 1.32 0.0145
7 ATM2 ANN Trai~ -0.796 11.7 9.89 -Inf Inf 1.64 1.33 0.0154
$ATM4
# A tibble: 7 x 11
ATM .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM4 ETS(`Lam Cash~ Trai~ -0.00734 9.93 7.84 -56.3 81.1 0.763 0.760 0.0822
2 ATM4 Multiplicativ~ Trai~ -0.194 9.94 7.80 -57.8 81.8 0.759 0.760 0.0901
3 ATM4 Damped Multip~ Trai~ -0.391 9.99 7.80 -58.7 82.1 0.759 0.765 0.0927
4 ATM4 Multiplicative Trai~ -0.733 10.0 7.83 -61.6 83.7 0.762 0.766 0.0876
5 ATM4 MNN Trai~ -0.00740 10.8 9.12 -70.5 97.2 0.888 0.828 0.0499
6 ATM4 ANN Trai~ 0.00654 10.8 9.12 -70.4 97.1 0.888 0.828 0.0499
7 ATM4 AAN Trai~ -0.178 11.0 9.23 -71.9 98.4 0.898 0.839 0.0489 compare_ETS %>%
map2(auto_ETS,
~ bind_rows(report(.x),glance(.y)) %>%
arrange(AICc)
)$ATM1
# A tibble: 7 x 10
ATM .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 ETS(`Lam Cash`) 3.26 -1288. 2595. 2596. 2634. 3.18 3.19 1.01
2 ATM1 Multiplicative No Tr~ 0.0445 -1327. 2674. 2675. 2713. 3.25 3.26 0.117
3 ATM1 Multiplicative 0.0445 -1328. 2679. 2680. 2726. 3.34 3.35 0.118
4 ATM1 Damped Multiplicative 0.0445 -1327. 2680. 2681. 2730. 3.31 3.32 0.118
5 ATM1 MNN 0.0860 -1468. 2941. 2941. 2953. 8.51 8.53 0.206
6 ATM1 ANN 8.56 -1468. 2941. 2941. 2953. 8.51 8.53 2.06
7 ATM1 AAN 8.78 -1471. 2953. 2953. 2972. 8.69 8.68 2.10
$ATM2
# A tibble: 7 x 10
ATM .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM2 ETS(`Lam Cash`) 56.5 -1808. 3637. 3637. 3676. 55.1 55.3 5.20
2 ATM2 Multiplicative 0.270 -1950. 3925. 3926. 3971. 107. 107. 0.406
3 ATM2 Multiplicative No T~ 0.238 -1954. 3928. 3929. 3967. 110. 111. 0.382
4 ATM2 Damped Multiplicati~ 0.280 -1957. 3940. 3941. 3990. 107. 107. 0.414
5 ATM2 AAN 136. -1971. 3953. 3953. 3972. 135. 135. 9.88
6 ATM2 ANN 138. -1975. 3956. 3956. 3968. 137. 137. 9.89
7 ATM2 MNN 0.280 -1977. 3961. 3961. 3972. 139. 139. 0.448
$ATM4
# A tibble: 7 x 10
ATM .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM4 ETS(`Lam Cash`) 0.194 -1907. 3834. 3835. 3873. 98.6 98.9 0.346
2 ATM4 Multiplicative No T~ 0.192 -1908. 3837. 3837. 3876. 98.7 99.0 0.341
3 ATM4 Damped Multiplicati~ 0.193 -1910. 3846. 3847. 3897. 99.9 100. 0.339
4 ATM4 Multiplicative 0.189 -1914. 3852. 3852. 3898. 100. 101. 0.337
5 ATM4 MNN 0.212 -1946. 3898. 3898. 3910. 117. 117. 0.387
6 ATM4 ANN 118. -1946. 3898. 3898. 3910. 117. 117. 9.12
7 ATM4 AAN 122. -1951. 3912. 3912. 3931. 120. 121. 9.23 5.2.3 Exponential Model Selection
With these results in mind, it would appear
that the automatically generated models by ETS to work best. Lowest AICc
scores and for the most part the lowest RMSE. The only exception is on
ATM1 with RMSE, but only by a tiny margin that is completely
negligible.
auto_ETS %>% map( ~.x %>%
pivot_longer(!ATM,names_to = "Model name",
values_to = "Orders")
)$ATM1
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM1 ETS(`Lam Cash`) <ETS(A,N,A)>
$ATM2
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM2 ETS(`Lam Cash`) <ETS(A,N,A)>
$ATM4
# A mable: 1 x 3
# Key: ATM, Model name [1]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM4 ETS(`Lam Cash`) <ETS(M,N,A)>5.3 ARIMA Modeling
5.3.1 Differencing
The ARIMA modeling process begins by figuring
out if the data is stationary or not. The ACF and PACF plots for ATM 1
and 2 do not appear to be stationary. On the other hand, ATM 4 shows
very few spikes that go outside the boundaries. It would seem that ATM 4
is already stationary without any differencing.
gg_tsdisplay(By_ATM$ATM1,`Lam Cash`,plot_type='partial')
gg_tsdisplay(By_ATM$ATM2,`Lam Cash`,plot_type='partial')
gg_tsdisplay(By_ATM$ATM4,`Lam Cash`,plot_type='partial')
Below I confirm whether ATM 4 is stationary
with the unitroot_kpss test. It would appear that it is
stationary.
By_ATM$ATM4 %>%
features(`Lam Cash`, unitroot_kpss) Meanwhile, ATM 1 and 2 appear to require one
seasonal difference. After applying that difference, the cash values for
each ATM appear to be stationary.
By_ATM$ATM1 %>%
features(`Lam Cash`, unitroot_nsdiffs)By_ATM$ATM2 %>%
features(`Lam Cash`, unitroot_nsdiffs)By_ATM$ATM1 %>%
mutate(dif = difference(`Lam Cash`,7)) %>%
features(dif, unitroot_ndiffs)By_ATM$ATM2 %>%
mutate(dif = difference(`Lam Cash`,7)) %>%
features(dif, unitroot_ndiffs)5.3.2 Auomatic Models
The ARIMA model building process is begun
the same way as ETS, by letting the algorithm choose the best model for
each, giving me a baseline to compare with.The automatically generated
models for each are shown below..
auto_ARIMA <- By_ATM %>%
map( ~ .x %>%
model(
stepwise = ARIMA(`Lam Cash`),
search = ARIMA(`Lam Cash`,stepwise = F))
)auto_ARIMA %>%
map( ~.x %>%
pivot_longer(!ATM,names_to = "Model name",
values_to = "Orders")
)$ATM1
# A mable: 2 x 3
# Key: ATM, Model name [2]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM1 stepwise <ARIMA(0,0,2)(0,1,1)[7]>
2 ATM1 search <ARIMA(0,0,2)(0,1,1)[7]>
$ATM2
# A mable: 2 x 3
# Key: ATM, Model name [2]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM2 stepwise <ARIMA(2,0,2)(0,1,1)[7]>
2 ATM2 search <ARIMA(5,0,0)(0,1,1)[7]>
$ATM4
# A mable: 2 x 3
# Key: ATM, Model name [2]
ATM `Model name` Orders
<chr> <chr> <model>
1 ATM4 stepwise <ARIMA(0,0,1)(2,0,0)[7] w/ mean>
2 ATM4 search <ARIMA(0,0,0)(2,0,0)[7] w/ mean>auto_ARIMA %>%
map( ~.x %>%
glance() %>%
arrange(AICc) %>%
select(ATM:BIC)
)$ATM1
# A tibble: 2 x 7
ATM .model sigma2 log_lik AIC AICc BIC
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 stepwise 3.17 -715. 1437. 1438. 1453.
2 ATM1 search 3.17 -715. 1437. 1438. 1453.
$ATM2
# A tibble: 2 x 7
ATM .model sigma2 log_lik AIC AICc BIC
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM2 stepwise 53.1 -1218. 2449. 2449. 2472.
2 ATM2 search 53.1 -1218. 2450. 2451. 2477.
$ATM4
# A tibble: 2 x 7
ATM .model sigma2 log_lik AIC AICc BIC
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM4 search 106. -1367. 2742. 2742. 2758.
2 ATM4 stepwise 105. -1366. 2742. 2742. 2762.5.3.3 Alternate Models
We test out a few other models for the
difference data and non difference data as shown below.
difference_models <- By_ATM[1:2] %>%
map( ~ .x %>%
model(
arima_000_110 = ARIMA(`Lam Cash` ~ pdq(0,0,0) + PDQ(1,1,0)),
arima_001_110 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,1,0)),
arima_001_111 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,1,1)),
arima_000_210 = ARIMA(`Lam Cash` ~ pdq(0,0,0) + PDQ(2,1,0)),
arima_100_012 = ARIMA(`Lam Cash` ~ pdq(1,0,0) + PDQ(0,1,2)),
arima_001_012 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,1,2)),
arima_002_210 = ARIMA(`Lam Cash` ~ pdq(0,0,2) + PDQ(2,1,0)),
arima_002_212 = ARIMA(`Lam Cash` ~ pdq(0,0,2) + PDQ(2,1,2)),
arima_005_011 = ARIMA(`Lam Cash` ~ pdq(0,0,5) + PDQ(0,1,1))
))difference_models %>%
map( ~.x %>%
glance() %>%
arrange(AICc) %>%
select(ATM:BIC)
)$ATM1
# A tibble: 9 x 7
ATM .model sigma2 log_lik AIC AICc BIC
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 arima_005_011 3.16 -713. 1440. 1440. 1467.
2 ATM1 arima_001_111 3.20 -716. 1441. 1441. 1456.
3 ATM1 arima_100_012 3.21 -717. 1442. 1442. 1458.
4 ATM1 arima_002_212 3.18 -714. 1442. 1443. 1470.
5 ATM1 arima_001_012 3.21 -716. 1443. 1443. 1462.
6 ATM1 arima_002_210 3.38 -725. 1460. 1460. 1479.
7 ATM1 arima_000_210 3.42 -728. 1462. 1462. 1474.
8 ATM1 arima_001_110 3.58 -736. 1478. 1478. 1490.
9 ATM1 arima_000_110 3.63 -739. 1482. 1482. 1489.
$ATM2
# A tibble: 9 x 7
ATM .model sigma2 log_lik AIC AICc BIC
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM2 arima_005_011 53.5 -1219. 2453. 2453. 2480.
2 ATM2 arima_001_111 55.9 -1228. 2464. 2464. 2480.
3 ATM2 arima_100_012 55.9 -1228. 2464. 2465. 2480.
4 ATM2 arima_002_212 55.4 -1225. 2465. 2465. 2492.
5 ATM2 arima_001_012 55.8 -1228. 2465. 2465. 2484.
6 ATM2 arima_000_210 57.8 -1234. 2475. 2475. 2486.
7 ATM2 arima_002_210 57.8 -1233. 2476. 2477. 2496.
8 ATM2 arima_000_110 59.6 -1240. 2484. 2484. 2492.
9 ATM2 arima_001_110 59.7 -1240. 2485. 2486. 2497.non_dif <- By_ATM$ATM4 %>%
model(
arima_000_100 = ARIMA(`Lam Cash` ~ pdq(0,0,0) + PDQ(1,0,0)),
arima_001_100 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,0,0)),
arima_001_101 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,0,1)),
arima_100_002 = ARIMA(`Lam Cash` ~ pdq(1,0,0) + PDQ(0,0,2)),
arima_001_002 = ARIMA(`Lam Cash` ~ pdq(0,0,1) + PDQ(1,0,2)),
arima_002_200 = ARIMA(`Lam Cash` ~ pdq(0,0,2) + PDQ(2,0,0)),
arima_002_202 = ARIMA(`Lam Cash` ~ pdq(0,0,2) + PDQ(2,0,2)),
arima_005_001 = ARIMA(`Lam Cash` ~ pdq(0,0,5) + PDQ(0,0,1))
)non_dif %>%
glance() %>%
arrange(AICc) %>%
select(ATM:BIC)5.3.4 ARIMA Model Selection
The automatic models are shown to be the best.
Both the stepwise and search models are almost exactly the same, so
either can be used for our ARIMA model for each ATM.
5.4 Comparison
In this section I will create a training data
set with one month from the data hidden. This will be used in order to
determine the final models used for each ATM.
train_ETS <-By_ATM %>%
map( ~.x %>%
filter_index(.~("2010-03-31")) %>%
model(ETS(`Lam Cash`))
)
train_ARIMA <- By_ATM %>%
map( ~.x %>%
filter_index(.~("2010-03-31")) %>%
model(ARIMA(`Lam Cash`))
)gg_tsresiduals(train_ETS$ATM1)
gg_tsresiduals(train_ETS$ATM2)
gg_tsresiduals(train_ETS$ATM4)
| The ETS Models for ATM 1 and 2 to have very low p-values. It might not be appropriate to use ETS for these ATM’s. ATM 4 appears to perform well enough as the p value indicates that the residuals are white noise.
train_ETS %>%
map(
~augment(.) %>%
features(.innov, ljung_box, lag=14,dof=4)
)$ATM1
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM1 ETS(`Lam Cash`) 21.3 0.0193
$ATM2
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM2 ETS(`Lam Cash`) 31.7 0.000450
$ATM4
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM4 ETS(`Lam Cash`) 12.6 0.249gg_tsresiduals(train_ARIMA$ATM1)
gg_tsresiduals(train_ARIMA$ATM2)
gg_tsresiduals(train_ARIMA$ATM4)
For ARIMA every models passes the test
therefore, all residuals appear to be white noise. I do note however,
that the value for ATM 4 is lower than the ETS model, possibly
indicating a weaker fit.
train_ARIMA %>%
map(
~augment(.) %>%
features(.innov, ljung_box, lag=14,dof=4)
)$ATM1
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM1 ARIMA(`Lam Cash`) 11.4 0.330
$ATM2
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM2 ARIMA(`Lam Cash`) 9.67 0.470
$ATM4
# A tibble: 1 x 4
ATM .model lb_stat lb_pvalue
<chr> <chr> <dbl> <dbl>
1 ATM4 ARIMA(`Lam Cash`) 15.6 0.110I use the accuracy function to test how well
each models fits the data. Both the ETS models and ARIMA models appear
to be closely matching, with ETS doing slightly better for ATM 4. With
the residual diagnostics from above in mind, it is clear which model
will be choose for each ATM.
pmap(list(train_ARIMA,train_ETS,By_ATM),
~bind_rows(
accuracy(..1),
accuracy(..2),
forecast( ..1,h = 31) %>% accuracy(..3),
forecast(..2,h = 31) %>% accuracy(..3)) %>%
select(-ME, -MPE, -ACF1)
) $ATM1
# A tibble: 4 x 8
ATM .model .type RMSE MAE MAPE MASE RMSSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM1 ARIMA(`Lam Cash`) Training 1.82 1.04 Inf 0.868 0.834
2 ATM1 ETS(`Lam Cash`) Training 1.84 1.04 Inf 0.870 0.846
3 ATM1 ARIMA(`Lam Cash`) Test 0.894 0.623 20.8 0.521 0.410
4 ATM1 ETS(`Lam Cash`) Test 0.895 0.624 20.9 0.521 0.411
$ATM2
# A tibble: 4 x 8
ATM .model .type RMSE MAE MAPE MASE RMSSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM2 ARIMA(`Lam Cash`) Training 7.31 5.11 Inf 0.825 0.812
2 ATM2 ETS(`Lam Cash`) Training 7.58 5.30 Inf 0.856 0.842
3 ATM2 ARIMA(`Lam Cash`) Test 6.14 4.47 Inf 0.722 0.682
4 ATM2 ETS(`Lam Cash`) Test 5.75 4.21 Inf 0.681 0.638
$ATM4
# A tibble: 4 x 8
ATM .model .type RMSE MAE MAPE MASE RMSSE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ATM4 ARIMA(`Lam Cash`) Training 10.3 8.61 86.4 0.826 0.779
2 ATM4 ETS(`Lam Cash`) Training 10.0 7.96 79.5 0.763 0.756
3 ATM4 ARIMA(`Lam Cash`) Test 9.36 7.47 124. 0.716 0.705
4 ATM4 ETS(`Lam Cash`) Test 8.36 6.52 108. 0.625 0.6305.5 Forecasting
5.5.1 Plotting Forecasts
Using results from the previous section I have
decided that we will use ARIMA for the first two ATM’s and will use an
ETS model for our fourth ATM. Shown below are the forecast generated
from the choose models for each ATM. There isn’t d anything unusual, or
clearly abnormal with any of these forecasts.
ATM_1_Forecast <- By_ATM$ATM1 %>%
model(ARIMA(`Lam Cash`)) %>%
forecast(h=31)
ATM_1_Forecast %>%
autoplot( By_ATM$ATM1,level=NULL) +
labs(title="Forecast for ATM 1")
ATM_2_Forecast <- By_ATM$ATM2 %>%
model(ARIMA(`Lam Cash`)) %>%
forecast(h=31)
ATM_2_Forecast %>%
autoplot( By_ATM$ATM2,level=NULL) +
labs(title="Forecast for ATM 2")
ATM_4_Forecast <- By_ATM$ATM4 %>%
model(ETS(`Lam Cash`)) %>%
forecast(h=31)
ATM_4_Forecast %>%
autoplot( By_ATM$ATM4,level=NULL) +
labs(title="Forecast for ATM 4")
5.5.2 Generating Predictions
In order to make more easily understood
predictions, I undo the Box-Cox Transformation using the mean value
forecast values. This generates the final May forecasts for this
assignment.
forecast_predictions <- bind_rows(
map2(By_ATM,list(ATM_1_Forecast,ATM_2_Forecast,ATM_4_Forecast),
~ .y %>%
mutate(Cash =inv_box_cox( .mean , features( .x, Cash, features = guerrero)$lambda_guerrero)) %>%
tibble()
)
)
forecast_predictions6 Conclusion
I was able to find appropriate models for each
of the ATM’s that seem to produce reasonable predictions. Using the map
function, I was able to easily produce results without having to repeat
much code. It seems that being mindful of the differences of each ATM
was very important in the end. This decision allowed me to discover that
different models worked best for each ATM. Using a single model would
not have generated proper results. This is something very important to
keep in mind when working with multivariate data. It would seem that
treating each as separate data frames is the most appropriate choice
unless you have reason to think otherwise.