DATA 624 Project 1 (A)

Project Introduction

The goal of this project is to forecast cash withdrawals for four ATMs from May 3, 2009 to May 3, 2010, using time series data provided in a single file. The data includes the variable ‘Cash,’ which represents the amount withdrawn, measured in hundreds of dollars. Given the limited initial information about these ATMs and the withdrawals, this analysis will adopt a business-oriented approach, with flexibility in data interpretation to reflect real-world forecasting challenges.

ATM machines

Loading the data that will be used for the analysis.

atm_data <- read_xlsx("ATM624Data.xlsx")
# tail(atm_data)

Since Data column is in excel serial time, we need to transform the variable to a readable format:

atm_data_tr <- atm_data |> group_by(ATM) |> mutate(Date= as.Date(DATE,origin = "1900-01-01")) |> select(Date,ATM,Cash) |> as_tsibble(key = ATM, index = Date)

atm_data_tr |> head() |>  kable() |> kable_styling() |>  kable_classic()

Date	ATM	Cash
2009-05-03	ATM1	96
2009-05-04	ATM1	82
2009-05-05	ATM1	85
2009-05-06	ATM1	90
2009-05-07	ATM1	99
2009-05-08	ATM1	88

Time Series Exploration

In this section, we examine the presence and impact of missing values within the ATM withdrawal dataset.

missing1 <- plot_intro(atm_data_tr, title = 'Missing Information on the ATM Cash Dataset',
           ggtheme = theme_minimal())

# Plot missing volume in Column 
missing2 <- plot_missing(atm_data_tr,title = 'Information about Missing Values in the ATM Cash dataset',ggtheme = theme_minimal())

grid.arrange(missing1, missing2, nrow = 1,layout_matrix= rbind(1,2))

Missing data accounts for only 1.29% of the total observations, a relatively small proportion that suggests a minor impact on overall data integrity. Nonetheless, handling missing values is essential to ensure that our forecasting models are robust and not influenced by gaps in the dataset.

Visualizing ATM Cash Withdrawal

atm_data_tr |> filter(!is.na(Cash)) |> autoplot(Cash) + facet_wrap(~ATM, scales = "free",ncol=2) +xlab("Daily Cash Withdrawal") +ggtitle("ATM Data Distributions without Missing Value")

atm_1 <- atm_data_tr |> filter(ATM=="ATM1") |> ungroup() |> select(-ATM)
atm_2 <- atm_data_tr |> filter(ATM=="ATM2")|> ungroup() |> select(-ATM)
atm_3 <- atm_data_tr |> filter(ATM=="ATM3")|> ungroup() |> select(-ATM)
atm_4 <- atm_data_tr |> filter(ATM=="ATM4")|> ungroup() |> select(-ATM)

plotResiduals <- function(data){
plot1 <- data |> autoplot(.innov) + labs(title = "innovation Residuals")
plot2 <-data |> ACF(.innov)  |> autoplot() + labs(title = "Auto Correlation Function")
plot3 <-data  |> ggplot(aes(x = .innov)) + geom_histogram(binwidth = bins_cal) + labs(title = "Histogram of residuals")
  
return( grid.arrange(plot1, plot2, plot3, nrow = 2,layout_matrix= rbind(1,c(2,3))))
}

bins_cal <- function(data){
  num_bins <- nclass.FD(data)
bin_width <- (max(data) - min(data)) / num_bins
return(bin_width)
}

get_IQR <- function(var){
  Q1 <- quantile(var, 0.25)
  Q3 <- quantile(var, 0.75)
  IQR_value <- IQR(var)
  print(Q1)
    print(Q3)
    print(IQR_value)

    df <- data.frame(Q1,
                   Q3,
                   IQR_value,
                   )
  return(df)
}


atm_data_tr |> filter(!is.na(Cash)) |> ggplot(aes(x = Cash)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + facet_wrap(~ATM, scales = "free",ncol = 2)

atm_data_tr |> filter(!is.na(Cash)) |> ggplot(aes(x = Cash)) +   geom_histogram(binwidth = bins_cal) + facet_wrap(~ATM, scales = "free",ncol = 2)

atm_data_tr |> filter(!is.na(Cash)) |> ggplot(aes(x=Date, y = Cash)) +   geom_point() + facet_wrap(~ATM, scales = "free",ncol = 2)

ATM 1 exhibits a notable data pattern, with a high concentration of outliers in the first quartile and a histogram that suggests a bimodal distribution. The histogram displays two distinct peaks: a smaller, left-skewed distribution in the lower range of values and a significantly larger peak on the right. This structure implies two different patterns or clusters of cash withdrawals, possibly indicating varying withdrawal behaviors or conditions over time.

The outliers in the lower quantile may reflect unusual low-transaction days or system anomalies that deviate from typical usage patterns. Meanwhile, the dominant peak in the higher range suggests a more consistent, likely primary, withdrawal pattern.

In contrast to ATM 1, ATM 2 shows a different pattern. ATM 2’s withdrawal data reveals a right-skewed, bimodal distribution with two distinct peaks, where the right peak has a slightly higher volume than the left. This indicates two common withdrawal patterns, with the larger right-side mode suggesting a preference for higher withdrawal amounts. Despite the right skew, the tails are relatively short, meaning most transactions are contained within a central range without extreme high or low values.

The nearly balanced but right-leaning bimodal structure points to consistent withdrawal behaviors with a slight tendency toward larger transactions. In forecasting, capturing both peaks accurately will be essential to represent ATM 2’s dual demand structure, especially given the greater frequency of higher withdrawals.

ATM 3, on the other hand, stands out due to its minimal variation. ATM 3 presents a unique data pattern, consisting almost entirely of zero values, with only three non-zero observations scattered throughout the dataset. These few non-zero values appear as extreme outliers relative to the distribution, given the otherwise consistent trend of zero withdrawals. The presence of only three non-zero entries suggests that ATM 3 was likely inactive or experienced minimal use, rendering the non-zero values as anomalies within an otherwise dormant dataset.

Given this structure, ATM 3’s data lacks the variability necessary for reliable forecasting. Attempting to model this series would yield predictions that are not meaningful in a real-world context, as the limited non-zero values do not provide enough information to infer any pattern or trend. Therefore, it is sensible to exclude ATM 3 from the forecasting analysis, focusing instead on ATMs with richer, more representative datasets.

Finally, ATM 4 displays a distinct challenge due to an extreme outlier. ATM 4 displays an extreme outlier, the cause of which is currently unknown. However, the substantial deviation of this rightmost value warrants further investigation. To assess the degree of this deviation, we will calculate the z-score, which will indicate how many standard deviations this value is from the mean. This analysis will help determine if it is an anomaly or a data entry error that needs to be addressed.

atm4_outlier <- atm_4 |> filter(Cash == max(Cash))

z_scores_manual <- (atm4_outlier$Cash -mean(atm_4$Cash)) / sd(atm_4$Cash)

atm4_outlier["Z-score"] <- z_scores_manual
atm4_outlier  |>  kable() |> kable_styling() |>  kable_classic()

Date	Cash	Z-score
2010-02-11	10919.76	16.04723

# sum(is.na(atm_4$Cash))

On February 11, 2010, a cash withdrawal of 10,919.76 was recorded, with a z-score of 16. This indicates that the value is 16 standard deviations above the mean, giving it an almost zero probability of occurrence. Such an extreme outlier suggests a one-time anomaly or potential recording error. Given its unique and isolated nature, removing this value will likely improve data accuracy without compromising the dataset.

To address this outlier, a more representative value will be imputed for February 11, 2010, in the following sections.

atm_4[atm_4$Date ==  atm4_outlier$Date[1],"Cash"] <- NA

After removing the outlier, the distribution will look as follows:

atm_4_box <- atm_4 |> filter(!is.na(Cash)) |> ggplot(aes(x = Cash)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 4")
atm_4_point <- atm_4 |> filter(!is.na(Cash)) |> ggplot(aes(x=Date, y = Cash)) +   geom_point() + ggtitle("ATM 4")
atm_4_hist <- atm_4 |> filter(!is.na(Cash)) |> ggplot(aes(x = Cash)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 4")
arr4 <- grid.arrange(atm_4_point,atm_4_box, atm_4_hist, nrow = 2,layout_matrix= rbind(1,c(2,3)))

Imputing Data

In this analysis step, multiple imputation is conducted on missing data within the ATM datasets to ensure complete and accurate time series for subsequent analysis. Using the mice package, Predictive Mean Matching (PMM) is applied, a method that fills in missing values by selecting observed values from similar cases, thereby maintaining realistic values within the range of existing data. Five imputations are generated for each ATM dataset, and one completed version is selected for analysis. This process is applied separately to atm_1, atm_2, and atm_4 to address any gaps, reducing potential biases associated with incomplete data. Each imputed dataset is then converted to tsibble format, indexed by Date, to support time series analysis. This approach preserves the structure and integrity of the datasets, ensuring robust and reliable forecasting outcomes.

atm_1_imputation <- mice(atm_1, m = 5, method = 'pmm', maxit = 5, seed = 123, printFlag = FALSE)
atm_1_c <- complete(atm_1_imputation, 1) |> as_tsibble(index=Date)

atm_2_imputation <- mice(atm_2, m = 5, method = 'pmm', maxit = 5, seed = 123, printFlag = FALSE)
atm_2_c <- complete(atm_2_imputation, 1) |> as_tsibble(index=Date)

atm_4_imputation <- mice(atm_4, m = 5, method = 'pmm', maxit = 5, seed = 123, printFlag = FALSE)
atm_4_c <- complete(atm_4_imputation, 1) |> as_tsibble(index=Date)

Visualizations after Imputed Data

atm_1_box <- atm_1_c |> ggplot(aes(x = Cash)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 1")
atm_1_point <-atm_1_c |> ggplot(aes(x=Date, y = Cash)) +   geom_point() + ggtitle("ATM 1")
atm_1_hist <- atm_1_c |> ggplot(aes(x = Cash)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 1")

atm_2_box <- atm_2_c |> ggplot(aes(x = Cash)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 2")
atm_2_point <-atm_2_c |> ggplot(aes(x=Date, y = Cash)) +   geom_point() + ggtitle("ATM 2")
atm_2_hist <- atm_2_c |> ggplot(aes(x = Cash)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 2")

atm_4_box <- atm_4_c |> ggplot(aes(x = Cash)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 4")
atm_4_point <- atm_4_c |> ggplot(aes(x=Date, y = Cash)) +   geom_point() + ggtitle("ATM 4")
atm_4_hist <- atm_4_c |> ggplot(aes(x = Cash)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 4")


arr1 <- grid.arrange(atm_1_point,atm_1_box, atm_1_hist, nrow = 2,layout_matrix= rbind(1,c(2,3)))

arr4 <- grid.arrange(atm_4_point,atm_4_box, atm_4_hist, nrow = 2,layout_matrix= rbind(1,c(2,3)))

arr2 <- grid.arrange(atm_2_point,atm_2_box, atm_2_hist, nrow = 2,layout_matrix= rbind(1,c(2,3)))

Data Transformation

To address the impact of outliers in the ATM withdrawal data for ATMs 1, 2, and 4, a Box-Cox transformation will be applied. Each of these ATMs presents unique challenges with outliers: ATM 1 has a concentration of low-value outliers, ATM 2 shows a right-skewed, bimodal distribution, and ATM 4 contains an extreme high-value outlier. These outliers can distort model outcomes, increase forecast errors, and obscure meaningful patterns within the data.

The Box-Cox transformation is well-suited for stabilizing variance and reducing skewness, as it compresses extreme values without sacrificing the integrity of the data’s structure. By adjusting data distributions closer to normality, the Box-Cox transformation will help reduce the influence of outliers and provide a more balanced dataset for forecasting.

cash_1_lambda <- atm_1_c |> features(Cash,features = guerrero) |> pull(lambda_guerrero)
cash_2_lambda <- atm_2_c |> features(Cash,features = guerrero) |> pull(lambda_guerrero)
cash_4_lambda <- atm_4_c |> features(Cash,features = guerrero) |> pull(lambda_guerrero)

# #Boxcox
atm_1_c$Cash_TR <- box_cox(atm_1_c$Cash, cash_1_lambda)
atm_2_c$Cash_TR <- box_cox(atm_2_c$Cash, cash_2_lambda)
atm_4_c$Cash_TR <- box_cox(atm_4_c$Cash, cash_4_lambda)

Data Post-Transformation

atm_1_boxc <- atm_1_c |> ggplot(aes(x = Cash_TR)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 1")
atm_1_pointc <-atm_1_c |> ggplot(aes(x=Date, y = Cash_TR)) +   geom_point() + ggtitle("ATM 1")
atm_1_histc <- atm_1_c |> ggplot(aes(x = Cash_TR)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 1")

atm_2_boxc <- atm_2_c |> ggplot(aes(x = Cash_TR)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 2")
atm_2_pointc <-atm_2_c |> ggplot(aes(x=Date, y = Cash_TR)) +   geom_point() + ggtitle("ATM 2")
atm_2_histc <- atm_2_c |> ggplot(aes(x = Cash_TR)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 2")

atm_4_boxc <- atm_4_c |> ggplot(aes(x = Cash_TR)) +  geom_boxplot(outlier.colour="red", outlier.shape=8,outlier.size=4) + ggtitle("ATM 4")
atm_4_pointc <- atm_4_c |> ggplot(aes(x=Date, y = Cash_TR)) +   geom_point() + ggtitle("ATM 4")
atm_4_histc <- atm_4_c |> ggplot(aes(x = Cash_TR)) +  geom_histogram(binwidth = bins_cal) + ggtitle("ATM 4")


grid.arrange(atm_1_pointc,atm_1_boxc, atm_1_histc, nrow = 2,layout_matrix= rbind(1,c(2,3)))

grid.arrange(atm_4_pointc,atm_4_boxc, atm_4_histc, nrow = 2,layout_matrix= rbind(1,c(2,3)))

grid.arrange(atm_2_pointc,atm_2_boxc, atm_2_histc, nrow = 2,layout_matrix= rbind(1,c(2,3)))

atm1_comp <- atm_1_c |> model(stl = STL(Cash_TR)) |> components()
atm2_comp <- atm_2_c |> model(stl = STL(Cash_TR)) |> components() 
atm4_comp <- atm_4_c |> model(stl = STL(Cash_TR)) |> components()

ATM 1

To gain insights into the ATM 1 data, STL decomposition will be applied to break down the time series into its seasonal, trend, and residual components, providing a clearer understanding of the underlying patterns in the data.

atm1_comp |> autoplot()

As part of this process, visualizing the ACF and PACF will aid in identifying any autocorrelation patterns and seasonal structures within the time series, which are essential for informing model selection and parameter tuning.

atm_1_c |> gg_tsdisplay(Cash_TR,plot_type = "partial",lag_max = 120)

Both the autocorrelation function (ACF) and partial autocorrelation function (PACF) were examined to uncover underlying patterns and seasonality in the data. The ACF plot reveals a distinct weekly seasonal pattern, with positive autocorrelations appearing at every 7-day interval, indicating that observations separated by a week tend to be similar. The autocorrelation at the first 7-day lag is the strongest, while subsequent weekly lags, such as 14, 21, and 27, display progressively weaker correlations, showing a gradual decay in influence. This slow decay in autocorrelation at these intervals suggests a persistent structure in the data that could indicate a trend, even if this trend does not stand out clearly in the STL decomposition’s trend component. This behavior implies that, despite the seasonal patterns, there may be underlying trends that are either subtle or partially absorbed by seasonal or noise components in the STL.

Additionally, the ACF plot shows clusters of negative autocorrelations between these multiples of 7, where only a few values surpass the lower confidence interval, hinting at a minor but significant clustering effect. The PACF complements this, showing a strong positive correlation around 0.65 at lag 7, consistent with the weekly seasonality observed in the ACF, along with weaker, significant negative correlations around -0.25 at other lags. These insights are instrumental for model selection, particularly in guiding differenciation choices for ARIMA models, where the observed 7-day seasonality suggests a seasonal differenciation component, and in refining the seasonal window for ETS models to capture this weekly pattern effectively.

Model Preparations

To forecast the amount of cash withdrawn from this ATM, several models will be implemented, each designed to capture the weekly seasonality and trends identified in the data. The models under consideration include a Seasonal Naïve model (SNAIVE), an ETS model with additive error, damped trend, and additive seasonality over a 7-day period, a Seasonal ARIMA model incorporating weekly seasonal differenciation (with parameters informed by this analysis), and a dynamic ARIMA model that will automatically determine its optimal parameters, likely selecting seasonal components to align with the weekly structure as needed.

For the ARIMA models, it is essential to ensure the data is stationary to produce reliable forecasts. We will evaluate the data’s stationarity by applying unit root tests, which assess the need for seasonal and non-seasonal differenciation, guiding our preparation of the dataset for effective ARIMA modeling.

bind_cols(
atm_1_c |>
  features(Cash_TR, unitroot_ndiffs),

atm_1_c |>
  features(Cash_TR, unitroot_nsdiffs)
) |>  kable() |> kable_styling() |>  kable_classic()

ndiffs	nsdiffs
0	1

Unit root tests indicate that only one seasonal difference is required to make the data stationary for the ARIMA model. Ensuring stationarity enhances model accuracy and reliability in forecasting. While ARIMA will be a primary model, we will also explore other forecasting models to determine the most effective approach for our data.

atm_1_c$diff <- difference(atm_1_c$Cash_TR,lag = 7)

atm_1_c |>
  features(diff, unitroot_kpss)  |>  kable() |> kable_styling() |>  kable_classic()

kpss_stat	kpss_pvalue
0.0153101	0.1

The KPSS test was conducted to assess the stationarity of the data. The KPSS statistic of 0.0153 and a p-value of 0.1 suggest that the data is likely stationary, as the high p-value (greater than common significance levels like 0.05 or 0.01) does not provide evidence to reject the null hypothesis of stationarity.

atm_1_c |> gg_tsdisplay(diff,plot_type='partial',lag_max = 50)

In examining the differenced data, both the ACF and PACF plots reveal key insights into the autoregressive structure within the dataset. The ACF displays a single significant autocorrelation at lag 7, indicating a weekly pattern that persists even after differenciation, with this autocorrelation showing on the negative side, suggesting a potentially dampened weekly relationship. The PACF also reflects a significant correlation at lag 7, followed by a gradual decay at multiples of 7 (such as lags 14, 21, etc.), with all correlations on the negative side.

Since both the ACF and PACF indicate autocorrelation at lag 7, it suggests that an autoregressive component is likely important in capturing the weekly structure in the model. This alignment across both charts highlights that a weekly autoregressive term could enhance the model’s ability to capture the underlying seasonality and is essential for the parameterization of the autoregressive part in our forecasting models.

Model Training

In this section, we describe the approach used to train and evaluate several forecasting models on our dataset. For training, we used 80% of the data, allowing for a robust evaluation on the remaining test set. The models explored include a Seasonal Naïve model (SNAIVE), an ETS model with additive error, damped trend, and additive seasonality over a 7-day period, a Seasonal ARIMA model with weekly seasonal differenciation, and a standard ARIMA model (ARIMA). By applying these models, we aim to identify the best-performing approach based on accuracy and generalizability.

atm_count1 <- atm_1_c |> nrow() * 0.8
atm_count1 <- floor(atm_count1)
atm_1_train <- atm_1_c |> select(-Cash) |> slice(1:atm_count1)

atm1_forecast_len <- nrow(atm_1_c) - atm_count1

atm1_models <- atm_1_train |> model(snaive=SNAIVE(Cash_TR),
                                AAdA=ETS(Cash_TR ~ error("A")+trend("Ad")+season("A",period=7)),
                                sarima = ARIMA(Cash_TR ~ pdq(0, 0, 0) + PDQ(0, 1, 1, period = 7)),
                                arima = ARIMA(Cash_TR)
                                )

atm_1_fit_snaive <- atm1_models |> select(snaive)
atm_1_fit_aada <- atm1_models |> select(AAdA)
atm_1_fit_sarima <- atm1_models |> select(sarima)
atm_1_fit_sarima_auto <- atm1_models |> select(arima)

Forecast

This section presents the forecasted values produced by each model, evaluated against the actual values in the test set. By comparing these forecasts with the observed data, we can later assess each model’s accuracy and ability to generalize beyond the training period.

atm1_forecast1 <-  atm_1_fit_snaive |> forecast(h=atm1_forecast_len)
atm1_forecast2 <- atm_1_fit_aada |> forecast(h=atm1_forecast_len)
atm1_forecast3 <-  atm_1_fit_sarima |> forecast(h=atm1_forecast_len)
atm1_forecast4 <-  atm_1_fit_sarima_auto |> forecast(h=atm1_forecast_len)

atm1_forecast1 |> autoplot(level = NULL) + autolayer(atm_1_c,Cash_TR)+ggtitle("Snaive Forecast")

atm1_forecast2 |> autoplot(level = NULL) + autolayer(atm_1_c,Cash_TR)+ggtitle("ETS (A,Ad,A) Forecast")

atm1_forecast3 |> autoplot(level = NULL) + autolayer(atm_1_c,Cash_TR)+ggtitle("SARIMA(0,0,0)(0,1,1)[7]")

atm1_forecast4 |> autoplot(level = NULL) + autolayer(atm_1_c,Cash_TR)+ggtitle("SARIMA Forecast (Automatic Parameter Selection)")

Model Evaluation

In this section, we evaluate the performance of various forecasting models using key accuracy metrics. The models assessed include Seasonal Naïve (snaive), ETS(AAA), Seasonal ARIMA (sarima), and ARIMA with automatic parameter calculation, with each model’s performance measured on both training and test datasets. Metrics considered are Mean Error (ME), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Percentage Error (MPE), Mean Absolute Percentage Error (MAPE), and Mean Absolute Scaled Error (MASE). These metrics provide a comprehensive view of each model’s accuracy and reliability in both training and test phases, highlighting model robustness and the ability to generalize to new data.

bind_rows(
  atm_1_fit_snaive |> accuracy(),
atm_1_fit_aada |> accuracy(),
atm_1_fit_sarima |> accuracy(),
atm_1_fit_sarima_auto |> accuracy(),
atm1_forecast1 |> accuracy(atm_1_c),
atm1_forecast2 |> accuracy(atm_1_c),
atm1_forecast3 |> accuracy(atm_1_c),
atm1_forecast4 |> accuracy(atm_1_c)) |> select(-MAPE) |> select(-MPE) |>  kable() |> kable_styling() |>  kable_classic()

.model	.type	ME	RMSE	MAE	MASE	RMSSE	ACF1
snaive	Training	-0.0087658	1.563949	0.8135924	1.0000000	1.0000000	0.0969108
AAdA	Training	-0.0357342	1.145855	0.6216602	0.7640929	0.7326679	0.1105523
sarima	Training	-0.0039664	1.165200	0.6459359	0.7939306	0.7450370	0.1331017
arima	Training	-0.0044840	1.153682	0.6376787	0.7837816	0.7376721	-0.0021455
snaive	Test	0.3445666	3.585271	2.7188384	3.3417695	2.2924474	-0.1617226
AAdA	Test	-0.0533964	3.405268	2.1982535	2.7019100	2.1773524	0.0277751
sarima	Test	-0.0887833	3.367566	2.2474552	2.7623846	2.1532455	-0.0435956
arima	Test	-0.0928764	3.369542	2.2430224	2.7569361	2.1545091	-0.0403931

Seasonal Naive Residuals

atm1_models |> select(snaive) |> report()

## Series: Cash_TR 
## Model: SNAIVE 
## 
## sigma^2: 2.4545

atm1_snaive <- atm1_models |> select(snaive) |> augment() 

plotResiduals(atm1_snaive)

ETS(A,Ad,A) Residuals

atm1_models |> select(AAdA) |> report()

## Series: Cash_TR 
## Model: ETS(A,Ad,A) 
##   Smoothing parameters:
##     alpha = 0.02894509 
##     beta  = 0.0001055486 
##     gamma = 0.0001008201 
##     phi   = 0.9686722 
## 
##   Initial states:
##      l[0]        b[0]      s[0]      s[-1]    s[-2]    s[-3]     s[-4]
##  7.932775 0.006368925 -3.992474 -0.1460253 1.199331 0.511737 0.9288118
##      s[-5]     s[-6]
##  0.9143184 0.5843007
## 
##   sigma^2:  1.3693
## 
##      AIC     AICc      BIC 
## 1763.124 1764.434 1810.922

atm1_aaa <- atm1_models |> select(AAdA) |> augment()

plotResiduals(atm1_aaa)

Seasonal ARIMA (0, 0, 0),(0, 1, 1)m=7 Residuals

atm1_models |> select(sarima) |> report()

## Series: Cash_TR 
## Model: ARIMA(0,0,0)(0,1,1)[7] 
## 
## Coefficients:
##          sma1
##       -0.9033
## s.e.   0.0350
## 
## sigma^2 estimated as 1.396:  log likelihood=-457.35
## AIC=918.71   AICc=918.75   BIC=926.01

atm1_sarima <- atm1_models |> select(sarima) |> augment()
plotResiduals(atm1_sarima)

Seasonal ARIMA (0, 0, 2),(0, 1, 1)m=7 Residuals (Automatic Parameter Selection)

atm1_arima <- atm1_models |> select(arima) |> report()

## Series: Cash_TR 
## Model: ARIMA(0,0,1)(0,1,1)[7] 
## 
## Coefficients:
##          ma1     sma1
##       0.1385  -0.9084
## s.e.  0.0599   0.0341
## 
## sigma^2 estimated as 1.373:  log likelihood=-454.71
## AIC=915.42   AICc=915.51   BIC=926.38

atm1_arima <- atm1_models |> select(arima) |> augment()
plotResiduals(atm1_arima)

Based on the provided metrics the Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) provide a direct measure of average prediction error magnitude, with lower values indicating a more accurate model. In the training data, the RTS (A,Ad,A) model demonstrates the lowest RMSE (1.146) and MAE (0.622), followed closely by the ARIMA model with an RMSE of 1.154 and MAE of 0.638. This pattern continues in the test set, where the RTS (A,Ad,A) model again achieves the lowest RMSE (3.405) and MAE (2.198), while the ARIMA and SARIMA models have slightly higher error values. These results suggest that the RTS (A,Ad,A) model captures the underlying structure of the data more effectively than the others, resulting in consistently lower errors across both datasets.

Moving beyond simple error measures, Mean Absolute Scaled Error (MASE) and Root Mean Squared Scaled Error (RMSSE) compare each model to a naive baseline, highlighting relative improvements. In the training data, the RTS (A,Ad,A) model achieves a MASE of 0.764 and an RMSSE of 0.733, signaling its efficiency not only in surpassing the naive baseline but also in performing more effectively than the ARIMA model (MASE of 0.784 and RMSSE of 0.738). This trend is maintained in the test results, where the RTS (A,Ad,A) model yields the lowest MASE (2.702) and RMSSE (2.177). These metrics reinforce the RTS (A,Ad,A) model’s ability to generalize well; it not only consistently outperforms the naive model but also manages to keep forecast errors low in both training and test sets.

The Autocorrelation of residuals at Lag 1 (ACF1) further informs our understanding of each model’s fit by assessing serial correlation in residuals, where lower values indicate a model that captures data dependencies without leaving patterns in the errors. On the training set, the ARIMA model exhibits the lowest autocorrelation (ACF1 = -0.002), showing a strong fit to the data. However, on the test set, the RTS (A,Ad,A) model achieves the lowest ACF1 (0.028), indicating it is particularly successful at avoiding autocorrelation in predictions on unseen data, which is essential for accurate forecasting.

In conclusion, the RTS (A,Ad,A) model stands out as the most effective overall. It demonstrates consistently lower error values across both training and test datasets, minimizes autocorrelation, and excels in generalizing to new data. While the ARIMA model performs well, particularly on the training data, the RTS (A,Ad,A) model’s steady performance across all metrics makes it the optimal choice, offering both accuracy and reliability for future forecasting.

Selected Model Forecasts*

 atm1_final_forecast <- atm_1_fit_aada |> forecast(h=atm1_forecast_len*2)

atm1_final_forecast |> autoplot() + autolayer(atm_1_train,Cash_TR) + labs(title = "ETS (A,Ad,A)")

atm1_final_forecast |> autoplot() + autolayer(atm_1_c,Cash_TR) + labs(title = "ETS (A,Ad,A)")

Saving Predictions

Forecast predictions are saved in both the original and transformed scales to ensure clarity and interpretability. While transformations improve model stability and performance, reverting predictions to the original scale makes results more accessible. This dual-saving approach maintains analytic precision and allows for straightforward interpretation of forecasted values.

cash1_predicted_values <-  data.frame(Date=atm1_final_forecast$Date,
                         Cash_Yhat=atm1_final_forecast$.mean,
                         Cash_OG_Scale= round(inv_box_cox(atm1_final_forecast$.mean, cash_1_lambda),2) )

write_xlsx(cash1_predicted_values,"ATM1_predicted.xlsx")

ATM 2

To gain insights into the ATM 2 data, STL decomposition will be applied to break down the time series into its seasonal, trend, and residual components, providing a clearer understanding of the underlying patterns in the data.

atm2_comp |> autoplot()

atm_2_c |> gg_tsdisplay(Cash_TR,plot_type = "partial",lag_max = 120)

Model Preparations

bind_cols(
atm_2_c |>
  features(Cash_TR, unitroot_ndiffs),

atm_2_c |>
  features(Cash_TR, unitroot_nsdiffs)
)|>  kable() |> kable_styling() |>  kable_classic()

ndiffs	nsdiffs
1	1

Unit root tests indicate that one seasonal difference and a first order difference are required to make the data stationary for the ARIMA model. Ensuring stationarity enhances model accuracy and reliability in forecasting. While ARIMA will be a primary model, we will also explore other forecasting models to determine the most effective approach for our data.

atm_2_c$diff <- difference(atm_2_c$Cash_TR,lag = 7)

atm_2_c |>
  features(diff, unitroot_kpss)|>  kable() |> kable_styling() |>  kable_classic()

kpss_stat	kpss_pvalue
0.0157759	0.1

The KPSS test was conducted to assess the stationarity of the data. The KPSS statistic of 0.01577 and a p-value of 0.1 suggest that the data is likely stationary. Since seasonal differenciation alone achieved stationarity, first-order differenciation will not be applied, as unnecessary differenciation could introduce additional noise into the data.

atm_2_c |> gg_tsdisplay(diff,plot_type='partial',lag_max = 50)

The analysis of the differenced data from the second ATM highlights a clear weekly structure without significant non-seasonal trends. The ACF reveals a single, notable negative autocorrelation at lag 7, pointing to a weekly pattern with a dampened effect rather than a strong or alternating weekly relationship. The PACF also shows a prominent negative correlation at lag 7, with a quick decay at subsequent multiples of 7 (lags 14, 21, and 28), indicating that the weekly influence weakens over time but remains significant. Importantly, neither the ACF nor PACF displays significant non-seasonal autocorrelations, suggesting that first-order differenciation is unnecessary for achieving stationarity. This allows the model to concentrate on the seasonal AR component to capture the weekly cycle accurately, maintaining model simplicity without additional non-seasonal terms.

Based on the findings, the suggested model specification would indeed be:

ARIMA(0,0,0)(1,1,1)[7]

Model Training

In this section, we describe the approach used to train and evaluate several forecasting models on our dataset. For training, we used 80% of the data, allowing for a robust evaluation on the remaining test set. The models explored include a Seasonal Naïve model (SNAIVE), an ETS model with additive error, damped trend, and additive seasonality over a 7-day period, another ETS model with additive error, additive trend, and additive seasonality over a 7-day period, a Seasonal ARIMA model with weekly seasonal differenciation, and a ARIMA model with automatic parameter selection. By applying these models, we aim to identify the best-performing approach based on accuracy and generalizability.

atm_count2 <- atm_2_c |> nrow() * 0.8
atm_count2 <- floor(atm_count2)
atm_2_train <- atm_2_c |> select(-Cash) |> slice(1:atm_count2)

atm2_forecast_len <- nrow(atm_2_c) - atm_count2

atm2_models <- atm_2_train |> model(
  snaive=SNAIVE(Cash_TR),
  AAdA=ETS(Cash_TR ~ error("A")+trend("Ad")+season("A",period=7)),
  AAA=ETS(Cash_TR ~ error("A")+trend("A")+season("A",period=7)),
  sarima = ARIMA(Cash_TR ~ pdq(0, 0, 0) + PDQ(1, 1, 1, period = 7)),
  arima = ARIMA(Cash_TR))

atm_2_fit_snaive <- atm2_models |> select(snaive)
atm_2_fit_aada <- atm2_models |> select(AAdA)
atm_2_fit_aaa <- atm2_models |> select(AAA)
atm_2_fit_sarima <- atm2_models |> select(sarima)
atm_2_fit_sarima_auto <- atm2_models |> select(arima)

Forecast

atm2_forecast1 <-  atm_2_fit_snaive |> forecast(h=atm2_forecast_len)
atm2_forecast2 <- atm_2_fit_aada |> forecast(h=atm2_forecast_len)
atm2_forecast3 <- atm_2_fit_aaa |> forecast(h=atm2_forecast_len)
atm2_forecast4 <-  atm_2_fit_sarima |> forecast(h=atm2_forecast_len)
atm2_forecast5 <-  atm_2_fit_sarima_auto |> forecast(h=atm2_forecast_len)

atm2_forecast1 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR) +ggtitle("Snaive Forecast")

atm2_forecast2 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("ETS (A,Ad,A) Forecast")

atm2_forecast3 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("ETS (A,A,A) Forecast")

atm2_forecast4 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("SARIMA(0,0,0)(1,1,1)[7]")

atm2_forecast5 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("SARIMA Forecast (Automatic Parameter Selection)")

Model Evaluation

In this section, we evaluate the performance of various forecasting models using key accuracy metrics. The models assessed include Seasonal Naïve (snaive),ETS(A,Ad,A), ETS(A,A,A), Seasonal ARIMA (sarima), and ARIMA, with each model’s performance measured on both training and test datasets. Metrics considered are Mean Error (ME), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Percentage Error (MPE), Mean Absolute Percentage Error (MAPE), and Mean Absolute Scaled Error (MASE). These metrics provide a comprehensive view of each model’s accuracy and reliability in both training and test phases, highlighting model robustness and the ability to generalize to new data.

bind_rows(
  atm_2_fit_snaive |> accuracy(),
  atm_2_fit_aada |> accuracy(),
atm_2_fit_aaa |> accuracy(),
atm_2_fit_sarima |> accuracy(),
atm_2_fit_sarima_auto |> accuracy(),
atm2_forecast1 |> accuracy(atm_2_c),
atm2_forecast2 |> accuracy(atm_2_c),
atm2_forecast3 |> accuracy(atm_2_c),
atm2_forecast4 |> accuracy(atm_2_c),
atm2_forecast5 |> accuracy(atm_2_c),
) |> select(-MPE) |> select(-MAPE)  |>  kable() |> kable_styling() |>  kable_classic()

.model	.type	ME	RMSE	MAE	MASE	RMSSE	ACF1
snaive	Training	-0.2014018	8.724314	5.969562	1.0000000	1.0000000	0.0756357
AAdA	Training	-0.7512715	6.301713	4.316326	0.7230558	0.7223163	0.0434730
AAA	Training	-0.0395356	6.255091	4.265839	0.7145984	0.7169723	0.0409818
sarima	Training	0.2112957	6.246644	4.175999	0.6995486	0.7160041	0.0381978
arima	Training	0.2112957	6.246644	4.175999	0.6995486	0.7160041	0.0381978
snaive	Test	6.4119893	19.352448	17.221288	2.8848496	2.2182201	-0.0999983
AAdA	Test	2.2820044	20.090270	17.668523	2.9597688	2.3027909	0.0324322
AAA	Test	2.7861691	20.208498	17.774306	2.9774891	2.3163425	0.0378059
sarima	Test	3.1511564	20.256205	18.384055	3.0796321	2.3218107	0.0098848
arima	Test	3.1511564	20.256205	18.384055	3.0796321	2.3218107	0.0098848

Seasonal Naive Residuals

atm2_models |> select(snaive) |> report()

## Series: Cash_TR 
## Model: SNAIVE 
## 
## sigma^2: 76.3409

atm2_snaive <- atm2_models |> select(snaive) |> augment() 

plotResiduals(atm2_snaive)

ETS(A,Ad,A) Residuals

atm2_models |> select(AAdA) |> report()

## Series: Cash_TR 
## Model: ETS(A,Ad,A) 
##   Smoothing parameters:
##     alpha = 0.006985195 
##     beta  = 0.003644821 
##     gamma = 0.0001004065 
##     phi   = 0.89318 
## 
##   Initial states:
##      l[0]     b[0]     s[0]     s[-1]    s[-2]    s[-3]     s[-4]    s[-5]
##  23.84392 0.307956 -18.2396 -9.431998 7.973254 3.154046 0.8653569 5.206796
##     s[-6]
##  10.47215
## 
##   sigma^2:  41.4135
## 
##      AIC     AICc      BIC 
## 2758.652 2759.961 2806.450

atm2_aada <- atm2_models |> select(AAdA) |> augment()

plotResiduals(atm2_aada)

ETS(A,A,A) Residuals

atm2_models |> select(AAA) |> report()

## Series: Cash_TR 
## Model: ETS(A,A,A) 
##   Smoothing parameters:
##     alpha = 0.01356248 
##     beta  = 0.0001000013 
##     gamma = 0.0001000351 
## 
##   Initial states:
##      l[0]        b[0]      s[0]     s[-1]    s[-2]    s[-3]    s[-4]    s[-5]
##  24.56583 -0.02452592 -18.35748 -9.378274 8.205675 3.118232 1.170683 4.922015
##     s[-6]
##  10.31915
## 
##   sigma^2:  40.6578
## 
##      AIC     AICc      BIC 
## 2752.315 2753.434 2796.436

atm2_aaa <- atm2_models |> select(AAA) |> augment()

plotResiduals(atm2_aaa)

Seasonal ARIMA (0, 0, 0),(1, 1, 1)m=7 Residuals

atm2_models |> select(sarima) |> report()

## Series: Cash_TR 
## Model: ARIMA(0,0,0)(1,1,1)[7] w/ drift 
## 
## Coefficients:
##          sar1     sma1  constant
##       -0.0466  -0.9034   -0.1956
## s.e.   0.0647   0.0301    0.0487
## 
## sigma^2 estimated as 40.4:  log likelihood=-936.21
## AIC=1880.42   AICc=1880.57   BIC=1895.03

atm2_sarima <- atm2_models |> select(sarima) |> augment()
plotResiduals(atm2_sarima)

ARIMA(0,0,0)(1,1,1)[7] with drift Residuals

atm2_models |> select(arima) |> report()

## Series: Cash_TR 
## Model: ARIMA(0,0,0)(1,1,1)[7] w/ drift 
## 
## Coefficients:
##          sar1     sma1  constant
##       -0.0466  -0.9034   -0.1956
## s.e.   0.0647   0.0301    0.0487
## 
## sigma^2 estimated as 40.4:  log likelihood=-936.21
## AIC=1880.42   AICc=1880.57   BIC=1895.03

atm2_arima <- atm2_models |> select(arima) |> augment()
plotResiduals(atm2_arima)

To determine the best model for forecasting, we evaluated several key metrics across both training and test datasets, focusing on RMSE (Root Mean Square Error), MAE (Mean Absolute Error), MASE (Mean Absolute Scaled Error), and ACF1 (autocorrelation of residuals at lag 1). In the training data, the SARIMA and ARIMA models achieved the lowest RMSE (6.25), MAE (4.18), and MASE (0.70), as well as the lowest ACF1 (0.038), indicating minimal residual autocorrelation and suggesting a strong fit. However, when applied to the test set, these models exhibited higher RMSE and MASE values (RMSE of 20.26 and MASE of 3.08), hinting at overfitting and a reduced ability to generalize to new data. In contrast, the Seasonal Naïve (snaive) model performed best on the test set in terms of raw error (RMSE of 19.35 and MAE of 17.22), though its high MASE (2.88) indicated less stability when scaled against the benchmark. Balancing training and test performance, the AAdA model emerged as the most consistent, with slightly higher test RMSE than the snaive model (20.09 versus 19.35) but a more favorable MASE (2.96), suggesting improved generalizability. Thus, while SARIMA and ARIMA excelled in the training data, the AAdA model offers a more balanced performance across both datasets, making it the most reliable choice for forecasting in this analysis, with the snaive model serving as a simple alternative.

Selected Model Forecasts*

 atm2_final_forecast <- atm_2_fit_aada |> forecast(h=atm2_forecast_len*2)

atm2_final_forecast |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR) + labs(title = "ETS (A,Ad,A)")

Saving Predictions

cash2_predicted_values <-  data.frame(Date=atm2_final_forecast$Date,
                         Cash_Yhat=atm2_final_forecast$.mean,
                         Cash_OG_Scale= round(inv_box_cox(atm2_final_forecast$.mean, cash_2_lambda),2) )

write_xlsx(cash1_predicted_values,"ATM2_predicted.xlsx")

ATM 4

To gain insights into the ATM 4 data, STL decomposition will be applied to break down the time series into its seasonal, trend, and residual components, providing a clearer understanding of the underlying patterns in the data.

atm4_comp |> autoplot()

atm_4_c |> gg_tsdisplay(Cash_TR,plot_type = "partial",lag_max = 50)

Following decomposition, the data for ATM 4 reveals a distinct weekly seasonal pattern, with residuals exhibiting characteristics of white noise, suggesting that seasonality effectively captures the underlying structure. The autocorrelation function (ACF) indicates strong autocorrelation at multiples of 7 lags, forming a wave-like pattern within the positive side of the confidence interval. This pattern confirms a consistent weekly cycle, where each week’s values correlate strongly with those from prior weeks at similar intervals.

The partial autocorrelation function (PACF) further supports this seasonal pattern, showing significant positive correlations at lags 7 and 14, with a smaller correlation at lag 21 that just exceeds the confidence interval. This indicates that while the weekly influence remains pronounced, it diminishes slightly over successive weekly lags, confirming a tapering seasonal effect.

The white noise quality of the residuals implies that the seasonal decomposition successfully accounts for most of the structure in ATM 4’s data, leaving minimal unexplained variation.

Model Preparations

bind_cols(
atm_4_c |>
  features(Cash_TR, unitroot_ndiffs),

atm_4_c |>
  features(Cash_TR, unitroot_nsdiffs)
) |>  kable() |> kable_styling() |>  kable_classic()

ndiffs	nsdiffs
0	0

Unit root tests indicate that no seasonal difference or order difference are required to make the data stationary for the ARIMA model. Therefore, among the models, we can try ****ARIMA(0,0,0)(0,0,2)[7]***.

Model Training

In this section, we describe the approach used to train and evaluate several forecasting models on our dataset. For training, we used 80% of the data, allowing for a robust evaluation on the remaining test set. The models explored include a Seasonal Naïve model (SNAIVE), an ETS model with additive error, damped trend, and additive seasonality over a 7-day period, another ETS model with additive error, additive trend, and additive seasonality over a 7-day period, a Seasonal ARIMA model with weekly seasonal differenciation, and a ARIMA model with automatic parameter selection. By applying these models, we aim to identify the best-performing approach based on accuracy and generalizability.

atm_count4 <- atm_4_c |> nrow() * 0.8
atm_count4 <- floor(atm_count4)
atm_4_train <- atm_4_c |> select(-Cash) |> slice(1:atm_count4)

atm4_forecast_len <- nrow(atm_4_c) - atm_count4

atm4_models <- atm_4_train |> model(
  snaive=SNAIVE(Cash_TR),
  AAdA=ETS(Cash_TR ~ error("A")+trend("Ad")+season("A",period=7)),
  AAA=ETS(Cash_TR ~ error("A")+trend("A")+season("A",period=7)),
  sarima = ARIMA(Cash_TR ~ 1 + pdq(2, 0, 1) + PDQ(1,0,1, period = 7)),
  arima = ARIMA(Cash_TR))

atm_4_fit_snaive <- atm4_models |> select(snaive)
atm_4_fit_aada <- atm4_models |> select(AAdA)
atm_4_fit_aaa <- atm4_models |> select(AAA)
atm_4_fit_sarima <- atm4_models |> select(sarima)
atm_4_fit_sarima_auto <- atm4_models |> select(arima)

Forecast

atm4_forecast1 <-  atm_4_fit_snaive |> forecast(h=atm4_forecast_len)
atm4_forecast2 <- atm_4_fit_aada |> forecast(h=atm4_forecast_len)
atm4_forecast3 <- atm_4_fit_aaa |> forecast(h=atm4_forecast_len)
atm4_forecast4 <-  atm_4_fit_sarima |> forecast(h=atm4_forecast_len)
atm4_forecast5 <-  atm_4_fit_sarima_auto |> forecast(h=atm4_forecast_len)

atm2_forecast1 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR) +ggtitle("Snaive Forecast")

atm2_forecast2 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("ETS (A,Ad,A) Forecast")

atm2_forecast3 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("ETS (A,A,A) Forecast")

atm2_forecast4 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("SARIMA(2,0,1)(1,0,1)[7]")

atm2_forecast5 |> autoplot(level = NULL) + autolayer(atm_2_c,Cash_TR)+ggtitle("SARIMA Forecast (Automatic Parameter Selection)")

Model Evaluation

In this section, we evaluate the performance of various forecasting models using key accuracy metrics. The models assessed include Seasonal Naïve (snaive),ETS(A,Ad,A), ETS(A,A,A), Seasonal ARIMA (sarima), and ARIMA with automatic parameter selection, with each model’s performance measured on both training and test datasets. Metrics considered are Mean Error (ME), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), Mean Percentage Error (MPE), Mean Absolute Percentage Error (MAPE), and Mean Absolute Scaled Error (MASE). These metrics provide a comprehensive view of each model’s accuracy and reliability in both training and test phases, highlighting model robustness and the ability to generalize to new data.

bind_rows(
  atm_4_fit_snaive |> accuracy(),
  atm_4_fit_aada |> accuracy(),
atm_4_fit_aaa |> accuracy(),
atm_4_fit_sarima |> accuracy(),
atm_4_fit_sarima_auto |> accuracy(),
atm4_forecast1 |> accuracy(atm_4_c),
atm4_forecast2 |> accuracy(atm_4_c),
atm4_forecast3 |> accuracy(atm_4_c),
atm4_forecast4 |> accuracy(atm_4_c),
atm4_forecast5 |> accuracy(atm_4_c),
) |> select(-MPE) |> select(-MAPE) |>  kable() |> kable_styling() |>  kable_classic()

.model	.type	ME	RMSE	MAE	MASE	RMSSE	ACF1
snaive	Training	0.0128818	17.49901	13.791648	1.0000000	1.0000000	0.0321332
AAdA	Training	-0.0316413	12.44117	9.914755	0.7188956	0.7109643	0.0757388
AAA	Training	-0.1078960	12.46447	9.943767	0.7209992	0.7122959	0.0715012
sarima	Training	-0.0162114	12.45118	10.056435	0.7291685	0.7115359	0.0032610
arima	Training	-0.0012250	13.80178	11.603466	0.8413401	0.7887175	0.0800296
snaive	Test	0.1476361	23.02361	18.945796	1.3737152	1.3157093	-0.3312529
AAdA	Test	0.3081586	17.26549	14.011429	1.0159358	0.9866553	0.0461409
AAA	Test	0.6714368	17.19633	13.984603	1.0139907	0.9827029	0.0551939
sarima	Test	0.1256425	16.97507	13.815929	1.0017605	0.9700586	0.0386921
arima	Test	-0.0421462	14.49909	11.655598	0.8451200	0.8285664	-0.0814381

Seasonal Naive Residuals

atm4_models |> select(snaive) |> report()

## Series: Cash_TR 
## Model: SNAIVE 
## 
## sigma^2: 307.2935

atm4_snaive <- atm4_models |> select(snaive) |> augment() 

plotResiduals(atm4_snaive)

ETS(A,Ad,A) Residuals

atm4_models |> select(AAdA) |> report()

## Series: Cash_TR 
## Model: ETS(A,Ad,A) 
##   Smoothing parameters:
##     alpha = 0.05214333 
##     beta  = 0.0001000293 
##     gamma = 0.0001005893 
##     phi   = 0.9609351 
## 
##   Initial states:
##      l[0]       b[0]      s[0]     s[-1]   s[-2]    s[-3]    s[-4]    s[-5]
##  33.65347 -0.1588139 -17.89534 -1.522121 3.59086 2.668411 4.431615 3.700939
##     s[-6]
##  5.025638
## 
##   sigma^2:  161.4163
## 
##      AIC     AICc      BIC 
## 3155.883 3157.192 3203.681

atm4_aada <- atm4_models |> select(AAdA) |> augment()

plotResiduals(atm4_aada)

ETS(A,A,A) Residuals

atm4_models |> select(AAA) |> report()

## Series: Cash_TR 
## Model: ETS(A,A,A) 
##   Smoothing parameters:
##     alpha = 0.05970484 
##     beta  = 0.0001000007 
##     gamma = 0.000100081 
## 
##   Initial states:
##      l[0]        b[0]      s[0]     s[-1]  s[-2]    s[-3]    s[-4]    s[-5]
##  33.56418 -0.00400619 -17.91052 -1.627621 3.3051 3.387347 4.327095 3.490388
##     s[-6]
##  5.028212
## 
##   sigma^2:  161.445
## 
##      AIC     AICc      BIC 
## 3154.976 3156.094 3199.097

atm4_aaa <- atm4_models |> select(AAA) |> augment()

plotResiduals(atm4_aaa)

Seasonal ARIMA SARIMA(2,0,1)(1,0,1)[7] with mean Residuals

atm4_models |> select(sarima) |> report()

## Series: Cash_TR 
## Model: ARIMA(2,0,1)(1,0,1)[7] w/ mean 
## 
## Coefficients:
##          ar1      ar2      ma1    sar1     sma1  constant
##       0.7759  -0.0039  -0.6821  0.9997  -0.9860    0.0018
## s.e.  0.1084   0.0544   0.1142  0.0001   0.0044    0.0002
## 
## sigma^2 estimated as 158.3:  log likelihood=-1158.72
## AIC=2331.43   AICc=2331.83   BIC=2357.17

atm4_sarima <- atm4_models |> select(sarima) |> augment()
plotResiduals(atm4_sarima)

ARIMA(0,0,0)(2,0,0)[7] with Mean Residuals

atm4_models |> select(arima) |> report()

## Series: Cash_TR 
## Model: ARIMA(0,0,0)(2,0,0)[7] w/ mean 
## 
## Coefficients:
##         sar1    sar2  constant
##       0.2340  0.1417   18.3669
## s.e.  0.0581  0.0584    0.7921
## 
## sigma^2 estimated as 192.5:  log likelihood=-1181.18
## AIC=2370.37   AICc=2370.51   BIC=2385.07

atm4_arima <- atm4_models |> select(arima) |> augment()
plotResiduals(atm4_arima)

The evaluation of forecasting models indicates that the ARIMA model provides the best balance of accuracy and generalizability, making it the preferred choice for predicting future values in this dataset. With the lowest RMSE (14.50) and MAE (11.66) on the test set, as well as favorable MASE (0.845) and RMSSE (0.829) values, the ARIMA model outperforms all other models in terms of minimizing forecast errors on new data. This strong test performance highlights its ability to generalize effectively beyond the training data, establishing it as a reliable option for forecasting.

The ARIMA model’s residual autocorrelation further supports its selection. With an ACF1 value of -0.081 on the test set, ARIMA demonstrates minimal residual correlation, indicating that it has captured the data’s structure well, with little unmodeled pattern remaining. This low residual autocorrelation suggests that ARIMA has effectively accounted for most predictive signals, enhancing the reliability of its future forecasts by minimizing bias from residual structure.

While the ARIMA model’s training set errors are slightly higher than those of the ETS models AAdA (ETS(A, Ad, A)) and AAA (ETS(A, A, A)), its RMSE of 13.80 and MAE of 11.60 remain competitive, showing a strong fit without overfitting. Unlike other models, which display higher errors on the test set relative to training, ARIMA’s consistent performance across datasets suggests it is well-suited to adapting to new data. This stability is crucial for long-term forecasting, where changing data conditions are expected.

Selected Model Forecasts*

In summary, the Seasonal ARIMA with the automatic parameter selection model is the most effective choice, combining low test set errors, minimal residual autocorrelation, and robust generalizability. Its balanced performance across key metrics makes it an optimal choice for capturing the dataset’s dynamics and producing accurate forecasts with minimal risk of overfitting or residual bias.

 atm4_final_forecast <- atm_4_fit_sarima_auto |> forecast(h=atm4_forecast_len*2)

atm4_final_forecast |> autoplot() + autolayer(atm_4_train,Cash_TR) + labs(title = "ARIMA(0,0,0)(2,0,0)[7] w/ mean")

atm4_final_forecast |> autoplot() + autolayer(atm_4_c,Cash_TR) + labs(title = "ARIMA(0,0,0)(2,0,0)[7] w/ mean")

Saving Predictions

cash4_predicted_values <-  data.frame(Date=atm4_final_forecast$Date,
                         Cash_Yhat=atm4_final_forecast$.mean,
                         Cash_OG_Scale= round(inv_box_cox(atm4_final_forecast$.mean, cash_4_lambda),2) )

write_xlsx(cash1_predicted_values,"ATM4_predicted.xlsx")

Conclusion

This project successfully demonstrated a comprehensive approach to forecasting ATM cash withdrawals using three distinct datasets. By systematically addressing data preparation, including outlier management, missing value imputation, and transformations to reduce skewness, each dataset was optimized for effective modeling. STL decomposition revealed key seasonal patterns, especially weekly cycles, which guided the selection of models suited to capture these trends.

Various forecasting models, such as Seasonal Naïve, ETS, and SARIMA, were evaluated based on their performance across training and test datasets using metrics like RMSE, MAE, and ACF. The ETS model with additive error, damped trend, and additive seasonality (ETS(A,Ad,A)) consistently emerged as the most reliable for future predictions, balancing accuracy and robustness across datasets. SARIMA models also showed potential in capturing specific seasonal characteristics, confirming that multiple approaches may be suitable depending on dataset nuances.

DATA 624 Project 1 (A)

Lewris Mota Sanchez

2024-10-19

Project Introduction

ATM machines

Time Series Exploration

Imputing Data

Data Transformation

ATM 1

Model Preparations

Model Training

Forecast

Model Evaluation

Saving Predictions

ATM 2

Model Preparations

Model Training

Forecast

Model Evaluation

Saving Predictions

ATM 4

Model Preparations

Model Training

Forecast

Model Evaluation

Saving Predictions

Conclusion