DATA_624_Project_1

Part A

In part A, I want you to 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. I am being somewhat ambiguous on purpose to make this have a little more business feeling. Explain and demonstrate your process, techniques used and not used, and your actual forecast. I am giving you data via an excel file, please provide your written report on your findings, visuals, discussion and your R code via an RPubs link along with the actual.rmd file Also please submit the forecast which you will put in an Excel readable file.

Reading in the data

library(fpp3)

── Attaching packages ──────────────────────────────────────────── fpp3 1.0.3 ──

✔ tibble      3.3.1     ✔ tsibble     1.2.0
✔ dplyr       1.2.0     ✔ tsibbledata 0.4.1
✔ tidyr       1.3.2     ✔ ggtime      1.0.0
✔ lubridate   1.9.5     ✔ feasts      0.5.0
✔ ggplot2     4.0.2     ✔ fable       0.5.0

── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
✖ lubridate::date()    masks base::date()
✖ dplyr::filter()      masks stats::filter()
✖ tsibble::intersect() masks base::intersect()
✖ tsibble::interval()  masks lubridate::interval()
✖ dplyr::lag()         masks stats::lag()
✖ tsibble::setdiff()   masks base::setdiff()
✖ tsibble::union()     masks base::union()

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.1     ✔ readr   2.2.0
✔ purrr   1.2.1     ✔ stringr 1.6.0

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter()     masks stats::filter()
✖ tsibble::interval() masks lubridate::interval()
✖ dplyr::lag()        masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readxl)
library(forecast)
library(lubridate)
library(writexl)
library(gridExtra)

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

library(imputeTS)

raw_data <- read_excel("ATM624Data.xlsx")

ATM_df <- raw_data |>
  mutate(DATE = as.Date(DATE, origin = "1899-12-30"))

head(ATM_df)

# A tibble: 6 × 3
  DATE       ATM    Cash
  <date>     <chr> <dbl>
1 2009-05-01 ATM1     96
2 2009-05-01 ATM2    107
3 2009-05-02 ATM1     82
4 2009-05-02 ATM2     89
5 2009-05-03 ATM1     85
6 2009-05-03 ATM2     90

Cleaning up the data

ATM_df |>
  filter(is.na(ATM)) |>
  arrange(desc(Cash)) |>
  select(DATE, ATM, Cash)

# A tibble: 14 × 3
   DATE       ATM    Cash
   <date>     <chr> <dbl>
 1 2010-05-01 <NA>     NA
 2 2010-05-02 <NA>     NA
 3 2010-05-03 <NA>     NA
 4 2010-05-04 <NA>     NA
 5 2010-05-05 <NA>     NA
 6 2010-05-06 <NA>     NA
 7 2010-05-07 <NA>     NA
 8 2010-05-08 <NA>     NA
 9 2010-05-09 <NA>     NA
10 2010-05-10 <NA>     NA
11 2010-05-11 <NA>     NA
12 2010-05-12 <NA>     NA
13 2010-05-13 <NA>     NA
14 2010-05-14 <NA>     NA

Looking at the data, it appears there are 14 rows with no ATM specified so I will filter those data points because they’re irrelevant for the forecasts.

ATM_df <- ATM_df |>
  filter(!is.na(ATM))

Plotting the data

ggplot(ATM_df, aes(x = DATE, y = Cash, color = ATM)) +
  geom_line() +
  facet_wrap(~ ATM, scales = "free_y", ncol = 2, axes = "all_x") +
  labs(title = "Daily Cash Withdrawals by ATM",
       x = "Date",
       y = "Cash (hundreds of dollars)") +
  theme_minimal() +
  theme(legend.position = "none")

From my initial observations of the plots for the 4 ATMs, a couple things stand out. ATM3 appears to have no withdrawals while the other ATMs have regular activity. ATM4 appears to have one spike in withdrawal where a lot of cash was taken out around February 2010. ATM1 and ATM2 both appear to be seasonal with ATM1 being more steady while ATM2 appears to have higher withdrawals in 2009 that declined in 2010. For the ATM forecasts, I will do them each separately.

ATM 1

ATM1_data <- ATM_df |>
  filter(ATM == "ATM1") |>
  arrange(DATE)

ggplot(ATM1_data, aes(x = DATE, y = Cash)) +
  geom_line(color = "steelblue") +
  labs(
    title = "Daily Cash Withdrawals for ATM 1",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

summary(ATM1_data)

      DATE                ATM                 Cash       
 Min.   :2009-05-01   Length:365         Min.   :  1.00  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 73.00  
 Median :2009-10-30   Mode  :character   Median : 91.00  
 Mean   :2009-10-30                      Mean   : 83.89  
 3rd Qu.:2010-01-29                      3rd Qu.:108.00  
 Max.   :2010-04-30                      Max.   :180.00  
                                         NA's   :3       

Looking at the plot for ATM 1, it appears to have weekly seasonality with no apparent negative or positive trend.

Grabbing the summary statistics for ATM 1, we can see that there are 3 missing values and 365 data points. In order to handle the 3 missing values, I’ll impute the data using the mean value because it’s only missing a few of the values and there are no extreme values.

ATM1_clean <- ATM1_data |>
  mutate(Cash = if_else(is.na(Cash), mean(Cash, na.rm = TRUE), Cash))
summary(ATM1_clean)

      DATE                ATM                 Cash       
 Min.   :2009-05-01   Length:365         Min.   :  1.00  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 73.00  
 Median :2009-10-30   Mode  :character   Median : 90.00  
 Mean   :2009-10-30                      Mean   : 83.89  
 3rd Qu.:2010-01-29                      3rd Qu.:108.00  
 Max.   :2010-04-30                      Max.   :180.00  

I’ll calculate the lambda to see if a box-cox transformation is necessary for stabilizing the variance.

ATM1_ts<- ATM1_clean |>
  as_tsibble(index = DATE)

ATM1_lambda <- ATM1_ts |>
  features(Cash, features = guerrero) |>
  pull(lambda_guerrero)
ATM1_lambda

[1] 0.2365721

Based on the lambda value, I will use a box-cox transformation to stabilize the variance. The data is not stationary since there is clear seasonality so I will use unit root tests to calculate the differences required to make the data stationary.

ATM1_ts |>
  mutate(bc_Cash = box_cox(Cash, ATM1_lambda)) |>
  features(bc_Cash, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       1      0

The result of the unit root tests showed that 1 seasonal diff is needed and since the data is weekly seasonal, I’ll use a lag of 7.

ATM1_ts |>
  mutate(d_bc_Cash = difference(box_cox(Cash, ATM1_lambda), lag = 7)) |>
  features(d_bc_Cash, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       0      0

Recalculating it after the seasonal difference, it appears like no additional difference is needed. Next I’ll check the plots for the data using gg_tsdisplay.

ATM1_ts |>
  mutate(d_bc_Cash = difference(box_cox(Cash, ATM1_lambda), lag = 7)) |>
  gg_tsdisplay(d_bc_Cash, plot_type = "partial")

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

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

Looking at the plots, the acf plot shows a spike at lag 7 and the pacf plot shows spikes at lag 7, 14, and 21. This confirms a strong weekly seasonality.

ATM1_ARIMA_fit <- ATM1_ts |>
  model(ARIMA(box_cox(Cash, ATM1_lambda)))
report(ATM1_ARIMA_fit)

Series: Cash 
Model: ARIMA(0,0,2)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM1_lambda) 

Coefficients:
         ma1      ma2     sma1
      0.0952  -0.1156  -0.6473
s.e.  0.0525   0.0529   0.0424

sigma^2 estimated as 1.496:  log likelihood=-580.54
AIC=1169.08   AICc=1169.2   BIC=1184.61

ATM1_ARIMA_fit |>
  gg_tsresiduals()

Using the Auto ARIMA function, I got an ARIMA(0,0,2)(0,1,1) model. Next I’ll try an ARIMA(0,0,0)(0,1,1) model to see how the plots compare.

ATM1_ARIMA_000_011_fit <- ATM1_ts |>
  model(ARIMA(box_cox(Cash, ATM1_lambda) ~ pdq(0,0,0) + PDQ(0,1,1)))
report(ATM1_ARIMA_000_011_fit)

Series: Cash 
Model: ARIMA(0,0,0)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM1_lambda) 

Coefficients:
         sma1
      -0.6372
s.e.   0.0425

sigma^2 estimated as 1.521:  log likelihood=-584.4
AIC=1172.81   AICc=1172.84   BIC=1180.57

ATM1_ARIMA_000_011_fit |>
  gg_tsresiduals()

The Auto ARIMA had a better AICc score and a better acf plot looks better than that of the ARIMA(0,0,0)(0,1,1) model. Next, I’ll try an auto ETS model and compare it just to check.

ATM1_ETS_fit <- ATM1_ts |>
  model(ETS(box_cox(Cash, ATM1_lambda)))
report(ATM1_ETS_fit)

Series: Cash 
Model: ETS(A,N,A) 
Transformation: box_cox(Cash, ATM1_lambda) 
  Smoothing parameters:
    alpha = 0.0001000375 
    gamma = 0.3516381 

  Initial states:
     l[0]      s[0]     s[-1]     s[-2]    s[-3]    s[-4]     s[-5]    s[-6]
 7.353716 -3.855517 0.4027751 0.3108518 1.002135 1.029524 0.6801793 0.430052

  sigma^2:  1.5304

     AIC     AICc      BIC 
2319.657 2320.278 2358.656 

ATM1_ETS_fit |>
  gg_tsresiduals()

ATM1_models <- ATM1_ts |>
  model(
    auto_ARIMA = ARIMA(box_cox(Cash, ATM1_lambda)),
    manual_ARIMA = ARIMA(box_cox(Cash, ATM1_lambda) ~ pdq(0,0,0) + PDQ(0,1,1)),
    auto_ETS = ETS(box_cox(Cash, ATM1_lambda))
  )

accuracy(ATM1_models) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)

# A tibble: 3 × 3
  .model        RMSE   MAE
  <chr>        <dbl> <dbl>
1 manual_ARIMA  25.0  15.9
2 auto_ARIMA    25.1  16.2
3 auto_ETS      25.1  16.2

This comparison of the RMSE and MAE values shows that the manual ARIMA better fits the historical data compared to the Auto ARIMA and Auto ETS. Even though the manual ARIMA has lower values, I think using the auto ARIMA to generate the forecasts is a better choice because it has a lower AICc value and its acf plot was within the blue dashed lines compared to the acf plot of the manual ARIMA.

ATM1_ARIMA_fit |>
  forecast(h = 31) |>
  autoplot(ATM1_ts) +
  labs(title = "ATM1 Cash forecasts")

ATM1_forecast <- ATM1_ARIMA_fit |>
  forecast(h = 31)

ATM1_export <- ATM1_forecast |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(cols = c(`80%`, `95%`)) |>
  as_tibble() |>
  select(DATE, .mean, `80%_lower`, `80%_upper`, `95%_lower`, `95%_upper`) |>
  rename(ATM1_Cash = .mean)

write_xlsx(ATM1_export, "ATM1_Cash_Forecast.xlsx")

ATM 2

ATM2_data <- ATM_df |>
  filter(ATM == "ATM2") |>
  arrange(DATE)

ggplot(ATM2_data, aes(x = DATE, y = Cash)) +
  geom_line(color = "steelblue") +
  labs(
    title = "Daily Cash Withdrawals for ATM 2",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

summary(ATM2_data)

      DATE                ATM                 Cash       
 Min.   :2009-05-01   Length:365         Min.   :  0.00  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 25.50  
 Median :2009-10-30   Mode  :character   Median : 67.00  
 Mean   :2009-10-30                      Mean   : 62.58  
 3rd Qu.:2010-01-29                      3rd Qu.: 93.00  
 Max.   :2010-04-30                      Max.   :147.00  
                                         NA's   :2       

Looking at the plot for ATM 2, it also appears to have weekly seasonality that has slightly declined over time.

The summary statistics for ATM 2 show that there are 2 missing values and 365 data points. In order to handle the 2 missing values, I will impute the data using the mean value because the data is only missing a few values with no extreme outliers.

ATM2_data <- ATM_df |>
  filter(ATM == "ATM2") |>
  mutate(Cash = if_else(is.na(Cash), mean(Cash, na.rm = TRUE), Cash))
summary(ATM2_data)

      DATE                ATM                 Cash       
 Min.   :2009-05-01   Length:365         Min.   :  0.00  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 26.00  
 Median :2009-10-30   Mode  :character   Median : 66.00  
 Mean   :2009-10-30                      Mean   : 62.58  
 3rd Qu.:2010-01-29                      3rd Qu.: 93.00  
 Max.   :2010-04-30                      Max.   :147.00  

ATM2_ts<- ATM2_data |>
  as_tsibble(index = DATE)

ATM2_lambda <- ATM2_ts |>
  features(Cash, features = guerrero) |>
  pull(lambda_guerrero)
ATM2_lambda

[1] 0.7240522

Based on the lambda value generated from the Guerrero method, I’ll a box-cox transformation to stabilize the variance.

Next I’ll use unit root tests to figure out how many differences is needed to make the data stationary. Data being stationary is a requirement for ARIMA models.

ATM2_ts |>
  mutate(bc_Cash = box_cox(Cash, ATM2_lambda)) |>
  features(bc_Cash, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       1      1

Based on the result, I will use 1 difference and 1 seasonal difference with lag 7 based on weekly seasonality.

ATM2_ts |>
  mutate(diff_Cash = difference(difference(box_cox(Cash, ATM2_lambda), lag = 7), differences = 1)) |>
  features(diff_Cash, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       0      0

ATM2_ts |>
  mutate(diff_Cash = difference(difference(box_cox(Cash, ATM2_lambda), lag = 7), differences = 1)) |>
  gg_tsdisplay(diff_Cash, plot_type = "partial")

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

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

The plot for the differenced and seasonal differenced ATM2 data flutuates around 0. The acf and pacf still have spikes that go past the blue dashed lines though. the notable spikes at lag 7 in both plots show that this is weekly seasonality.

ATM2_ARIMA_fit <- ATM2_ts |>
  model(ARIMA(box_cox(Cash, ATM2_lambda)))
report(ATM2_ARIMA_fit)

Series: Cash 
Model: ARIMA(2,0,2)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM2_lambda) 

Coefficients:
          ar1      ar2     ma1     ma2     sma1
      -0.4245  -0.9067  0.4682  0.7993  -0.7251
s.e.   0.0587   0.0440  0.0889  0.0586   0.0409

sigma^2 estimated as 67.59:  log likelihood=-1261.81
AIC=2535.62   AICc=2535.86   BIC=2558.91

ATM2_ARIMA_fit |>
  gg_tsresiduals()

Using the auto ARIMA function, it used a ARIMA(2,0,2)(0,1,1)[7] model. Next I’ll try out some other ARIMA models.

ATM2_ARIMA_000_011_fit <- ATM2_ts |>
  model(ARIMA(box_cox(Cash, ATM2_lambda) ~ pdq(0,0,0) + PDQ(0,1,1)))
report(ATM2_ARIMA_000_011_fit)

Series: Cash 
Model: ARIMA(0,0,0)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM2_lambda) 

Coefficients:
         sma1
      -0.6041
s.e.   0.0425

sigma^2 estimated as 71.08:  log likelihood=-1272.28
AIC=2548.56   AICc=2548.59   BIC=2556.32

ATM2_ARIMA_000_011_fit |>
  gg_tsresiduals()

ATM2_ARIMA_111_011_fit <- ATM2_ts |>
  model(ARIMA(box_cox(Cash, ATM2_lambda) ~ pdq(1,1,1) + PDQ(0,1,1)))
report(ATM2_ARIMA_111_011_fit)

Series: Cash 
Model: ARIMA(1,1,1)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM2_lambda) 

Coefficients:
         ar1      ma1     sma1
      0.0340  -1.0000  -0.6004
s.e.  0.0531   0.0087   0.0432

sigma^2 estimated as 71.58:  log likelihood=-1272.75
AIC=2553.51   AICc=2553.62   BIC=2569.02

ATM2_ARIMA_111_011_fit |>
  gg_tsresiduals()

ATM2_ARIMA_002_011_fit <- ATM2_ts |>
  model(ARIMA(box_cox(Cash, ATM2_lambda) ~ pdq(0,0,2) + PDQ(0,1,1)))
report(ATM2_ARIMA_002_011_fit)

Series: Cash 
Model: ARIMA(0,0,2)(0,1,1)[7] 
Transformation: box_cox(Cash, ATM2_lambda) 

Coefficients:
         ma1      ma2     sma1
      0.0299  -0.0922  -0.6271
s.e.  0.0526   0.0511   0.0434

sigma^2 estimated as 70.67:  log likelihood=-1270.41
AIC=2548.83   AICc=2548.94   BIC=2564.35

ATM2_ARIMA_002_011_fit |>
  gg_tsresiduals()

ATM2_models <- ATM2_ts |>
  model(
    auto_ARIMA = ARIMA(box_cox(Cash, ATM2_lambda)),
    auto_ETS = ETS(box_cox(Cash, ATM2_lambda))
  )

accuracy(ATM2_models) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)

# A tibble: 2 × 3
  .model      RMSE   MAE
  <chr>      <dbl> <dbl>
1 auto_ARIMA  24.4  17.1
2 auto_ETS    25.3  17.8

Since the Auto ARIMA had a better AICc score compared to my manual ARIMA models and a lower RMSE than the Auto ETS model, I’ll use the Auto ARIMA model to generate the forecast for ATM2.

ATM2_ARIMA_fit |>
  gg_tsresiduals()

Looking at the plots for the auto ARIMA model above, it does resemble white noise so I’ll use the ARIMA model to generate the forecasts for ATM 1.

ATM2_ARIMA_fit |>
  forecast(h = 31) |>
  autoplot(ATM2_ts) +
  labs(title = "ATM 2 Cash forecasts")

ATM2_forecast <- ATM2_ARIMA_fit |>
  forecast(h = 31)

ATM2_export <- ATM2_forecast |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(cols = c(`80%`, `95%`)) |>
  as_tibble() |>
  select(DATE, .mean, `80%_lower`, `80%_upper`, `95%_lower`, `95%_upper`) |>
  rename(ATM2_Cash = .mean)

write_xlsx(ATM2_export, "ATM2_Cash_Forecast.xlsx")

ATM 3

ATM3_data <- ATM_df |>
  filter(ATM == "ATM3") |>
  arrange(DATE)

ggplot(ATM3_data, aes(x = DATE, y = Cash)) +
  geom_line(color = "steelblue") +
  labs(
    title = "Daily Cash Withdrawals for ATM 3",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

summary(ATM3_data)

      DATE                ATM                 Cash        
 Min.   :2009-05-01   Length:365         Min.   : 0.0000  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 0.0000  
 Median :2009-10-30   Mode  :character   Median : 0.0000  
 Mean   :2009-10-30                      Mean   : 0.7205  
 3rd Qu.:2010-01-29                      3rd Qu.: 0.0000  
 Max.   :2010-04-30                      Max.   :96.0000  

Looking at the plot for ATM 3, it appears to be mostly 0 values with only a few data points at the very end. In order to work with the data, I’ll filter the 0s out first because they will bring down the forecast.

ATM3_nonzero <- ATM3_data |>
  filter(Cash > 0)

nrow(ATM3_nonzero)

[1] 3

ATM3_nonzero_ts <- ATM3_nonzero |>
  as_tsibble(index = DATE)

For this data set, I believe using either a simple mean model or a naive model will be the best way forward since there are only 3 data points.

ATM3_mean_fit <- ATM3_nonzero_ts |>
  model(mean_model = MEAN(Cash))
glance(ATM3_mean_fit)

# A tibble: 1 × 2
  .model     sigma2
  <chr>       <dbl>
1 mean_model   54.3

ATM3_naive_fit <- ATM3_nonzero_ts |>
  model(naive_model = NAIVE(Cash))
glance(ATM3_naive_fit)

# A tibble: 1 × 2
  .model      sigma2
  <chr>        <dbl>
1 naive_model   144.

Looking at the two sigma2 values above, I believe the simple mean model will work better since it has a lower sigma2 value.

ATM3_forecast <- ATM3_mean_fit |>
  forecast(h = 31)

autoplot(ATM3_forecast, ATM3_nonzero_ts)

ATM3_forecast <- ATM3_mean_fit |>
  forecast(h = 31)

ATM3_export <- ATM3_forecast |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(cols = c(`80%`, `95%`)) |>
  as_tibble() |>
  select(DATE, .mean, `80%_lower`, `80%_upper`, `95%_lower`, `95%_upper`) |>
  rename(ATM3_Cash = .mean)

write_xlsx(ATM3_export, "ATM3_Cash_Forecast.xlsx")

ATM 4

ATM4_data <- ATM_df |>
  filter(ATM == "ATM4") |>
  arrange(DATE)

ggplot(ATM4_data, aes(x = DATE, y = Cash)) +
  geom_line(color = "steelblue") +
  labs(
    title = "Daily Cash Withdrawals for ATM 4",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

summary(ATM4_data)

      DATE                ATM                 Cash          
 Min.   :2009-05-01   Length:365         Min.   :    1.563  
 1st Qu.:2009-07-31   Class :character   1st Qu.:  124.334  
 Median :2009-10-30   Mode  :character   Median :  403.839  
 Mean   :2009-10-30                      Mean   :  474.043  
 3rd Qu.:2010-01-29                      3rd Qu.:  704.507  
 Max.   :2010-04-30                      Max.   :10919.762  

Looking at the data plot for ATM 4, there is one very big anomaly so I’ll impute the data with the median value because of that extreme outlier value.

ATM4_data <- ATM_df |>
  filter(ATM == "ATM4") |>
  mutate(Cash = if_else(Cash > 5000, median(Cash), Cash))
summary(ATM4_data)

      DATE                ATM                 Cash         
 Min.   :2009-05-01   Length:365         Min.   :   1.563  
 1st Qu.:2009-07-31   Class :character   1st Qu.: 124.334  
 Median :2009-10-30   Mode  :character   Median : 403.839  
 Mean   :2009-10-30                      Mean   : 445.233  
 3rd Qu.:2010-01-29                      3rd Qu.: 704.192  
 Max.   :2010-04-30                      Max.   :1712.075  

ggplot(ATM4_data, aes(x = DATE, y = Cash)) +
  geom_line(color = "steelblue") +
  labs(
    title = "Daily Cash Withdrawals for ATM 4",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

ATM4_ts<- ATM4_data |>
  as_tsibble(index = DATE)

Next I’ll use the Guerrero method to calculate the lambda to see if a box-cox transformation is needed.

ATM4_lambda <- ATM4_ts |>
  features(Cash, features = guerrero) |>
  pull(lambda_guerrero)
ATM4_lambda

[1] 0.4525693

Based on the lambda value, I’ll use a box-cox transformation to stabilize the variance and use the unit root tests to check if a difference or seasonal difference is needed to make the data stationary.

ATM4_ts |>
  mutate(bc_Cash = box_cox(Cash, ATM4_lambda)) |>
  features(bc_Cash, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       0      0

Based on the unit root tests, neither difference is needed.

ATM4_ts |>
  mutate(bc_Cash = box_cox(Cash, ATM4_lambda)) |>
  gg_tsdisplay(bc_Cash, plot_type = "partial")

ATM4_ARIMA_fit <- ATM4_ts |>
  model(ARIMA(box_cox(Cash, ATM4_lambda)))
report(ATM4_ARIMA_fit)

Series: Cash 
Model: ARIMA(0,0,1)(2,0,0)[7] w/ mean 
Transformation: box_cox(Cash, ATM4_lambda) 

Coefficients:
         ma1    sar1    sar2  constant
      0.0769  0.2088  0.1996   17.1612
s.e.  0.0528  0.0516  0.0525    0.7422

sigma^2 estimated as 182.2:  log likelihood=-1466.35
AIC=2942.7   AICc=2942.87   BIC=2962.2

ATM4_ARIMA_fit |>
  gg_tsresiduals()

ATM4_ARIMA_000_200_fit <- ATM4_ts |>
  model(ARIMA(box_cox(Cash, ATM4_lambda) ~ pdq(0,0,0) + PDQ(2,0,0)))
report(ATM4_ARIMA_000_200_fit)

Series: Cash 
Model: ARIMA(0,0,0)(2,0,0)[7] w/ mean 
Transformation: box_cox(Cash, ATM4_lambda) 

Coefficients:
        sar1    sar2  constant
      0.2124  0.1892   17.3553
s.e.  0.0517  0.0521    0.6919

sigma^2 estimated as 182.8:  log likelihood=-1467.41
AIC=2942.82   AICc=2942.93   BIC=2958.42

Since the Auto ARIMA has basically the same AICc value and is simpler, I’ll compare the auto ARIMA to an auto ETS model.

ATM4_models <- ATM4_ts |>
  model(
    auto_arima = ARIMA(box_cox(Cash, ATM4_lambda)),
    auto_ets = ETS(box_cox(Cash, ATM4_lambda))
  )

accuracy(ATM4_models) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)

# A tibble: 2 × 3
  .model      RMSE   MAE
  <chr>      <dbl> <dbl>
1 auto_ets    344.  262.
2 auto_arima  352.  274.

Since the auto ETS had a better RMSE and MAE score than the auto ARIMA, I’ll use that to generate the forecast after I verify that the residuals resemble white noise.

ATM4_ETS_fit <- ATM4_ts |>
  model(ETS(box_cox(Cash, ATM4_lambda)))
glance(ATM4_ETS_fit)

# A tibble: 1 × 9
  .model                      sigma2 log_lik   AIC  AICc   BIC   MSE  AMSE   MAE
  <chr>                        <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ETS(box_cox(Cash, ATM4_lam…  0.223  -2006. 4033. 4033. 4072.  171.  171. 0.370

ATM4_ETS_fit |>
  gg_tsresiduals()

Looking at the plots above, it does resemble white noise so I’ll use the auto ETS model to generate the forecasts for ATM 4.

ATM4_ETS_fit |>
  forecast(h = 31) |>
  autoplot(ATM4_ts) +
  labs(title = "ATM 4 Cash forecasts")

ATM4_forecast <- ATM4_ETS_fit |>
  forecast(h = 31)

ATM4_export <- ATM4_forecast |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(cols = c(`80%`, `95%`)) |>
  as_tibble() |>
  select(DATE, .mean, `80%_lower`, `80%_upper`, `95%_lower`, `95%_upper`) |>
  rename(ATM4_Cash = .mean)

write_xlsx(ATM4_export, "ATM4_Cash_Forecast.xlsx")

Part B

Part B consists of a simple dataset of residential power usage for January 1998 until December 2013. Your assignment is to model these data and a monthly forecast for 2014. The data is given in a single file. The variable ‘KWH’ is power consumption in Kilowatt hours, the rest is straight forward. Add this to your existing files above.

raw_data <- read_excel("ResidentialCustomerForecastLoad-624.xlsx")

power_df <- raw_data |>
  mutate(`YYYY-MMM` = yearmonth(`YYYY-MMM`))
summary(power_df)

  CaseSequence      YYYY-MMM             KWH          
 Min.   :733.0   Min.   :1998 Jan   Min.   :  770523  
 1st Qu.:780.8   1st Qu.:2001 Dec   1st Qu.: 5429912  
 Median :828.5   Median :2005 Dec   Median : 6283324  
 Mean   :828.5   Mean   :2005 Dec   Mean   : 6502475  
 3rd Qu.:876.2   3rd Qu.:2009 Dec   3rd Qu.: 7620524  
 Max.   :924.0   Max.   :2013 Dec   Max.   :10655730  
                                    NA's   :1         

power_ts <- power_df |>
  mutate(Month = yearmonth(`YYYY-MMM`)) |>
  as_tsibble(index = Month)

autoplot(power_ts, KWH) +
  labs(
    title = "Monthly Residential Power Usage (1998 - 2013)",
    y = "Power Consumption (KWH)",
    x = "Time"
  )

There is one missing value so I’ll impute the data with the median because there is one extreme value.

power_ts_imputed <- power_ts |>
  mutate(KWH = replace_na(KWH, median(KWH, na.rm = TRUE)))
summary(power_ts_imputed)

  CaseSequence      YYYY-MMM             KWH               Month         
 Min.   :733.0   Min.   :1998 Jan   Min.   :  770523   Min.   :1998 Jan  
 1st Qu.:780.8   1st Qu.:2001 Dec   1st Qu.: 5434539   1st Qu.:2001 Dec  
 Median :828.5   Median :2005 Dec   Median : 6283324   Median :2005 Dec  
 Mean   :828.5   Mean   :2005 Dec   Mean   : 6501333   Mean   :2005 Dec  
 3rd Qu.:876.2   3rd Qu.:2009 Dec   3rd Qu.: 7608792   3rd Qu.:2009 Dec  
 Max.   :924.0   Max.   :2013 Dec   Max.   :10655730   Max.   :2013 Dec  

autoplot(power_ts_imputed, KWH) +
  labs(
    title = "Monthly Residential Power Usage (1998 - 2013)",
    y = "Power Consumption (KWH)",
    x = "Time"
  )

For the outlier value, I’ll impute it with the median value.

power_ts_median_imputed <- power_ts_imputed |>
  mutate(
    KWH = ifelse(
      KWH < 2000000, 
      median(KWH[KWH >= 2000000], na.rm = TRUE), 
      KWH
    )
  )

autoplot(power_ts_median_imputed, KWH) +
  labs(
    title = "Monthly Residential Power Usage (1998 - 2013)",
    y = "Power Consumption (KWH)",
    x = "Time"
  )

I’ll calculate the lambda to see if a box-cox transformation is necessary for stabilizing the variance.

Power_lambda <- power_ts_median_imputed |>
  features(KWH, features = guerrero) |>
  pull(lambda_guerrero)
Power_lambda

[1] -0.2134793

Based on the resulting lambda value, I will use a box-cox transformation to stabilize the variance. The data is not stationary since there is clear seasonality so I will use unit root tests to calculate the differences required to make the data stationary.

power_ts_median_imputed |>
  mutate(bc_KWH = box_cox(KWH, Power_lambda)) |>
  features(bc_KWH, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       1      1

Based on the unit root test result, I will use 1 difference and 1 seasonal difference with lag 12 based on annual seasonality.

power_ts_median_imputed |>
  mutate(diff_KWH = difference(difference(box_cox(KWH, Power_lambda), lag = 12), differences = 1)) |>
  features(diff_KWH, list(unitroot_nsdiffs, unitroot_ndiffs))

# A tibble: 1 × 2
  nsdiffs ndiffs
    <int>  <int>
1       0      0

power_ts_median_imputed |>
  mutate(diff_KWH = difference(difference(box_cox(KWH, Power_lambda), lag = 12), differences = 1)) |>
  gg_tsdisplay(diff_KWH, plot_type = "partial")

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

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

Looking at the acf plot, the strong spike at lag 12 confirms that there is annual seasonality.

Power_ARIMA_fit <- power_ts_median_imputed |>
  model(ARIMA(box_cox(KWH, Power_lambda)))
report(Power_ARIMA_fit)

Series: KWH 
Model: ARIMA(1,0,0)(2,1,0)[12] w/ drift 
Transformation: box_cox(KWH, Power_lambda) 

Coefficients:
         ar1     sar1     sar2  constant
      0.2569  -0.7427  -0.3819     9e-04
s.e.  0.0748   0.0733   0.0747     3e-04

sigma^2 estimated as 1.234e-05:  log likelihood=770.62
AIC=-1531.25   AICc=-1530.9   BIC=-1515.28

Power_ARIMA_fit |>
  gg_tsresiduals()

The acf plot for the auto ARIMA has a couple spikes at 3 and 6 that go past the blue dashed lines but the rest are within the blue dashed lines. Next I’ll check an auto ETS model.

Power_ETS_fit <- power_ts_median_imputed |>
  model(ETS(box_cox(KWH, Power_lambda)))
report(Power_ETS_fit)

Series: KWH 
Model: ETS(M,N,M) 
Transformation: box_cox(KWH, Power_lambda) 
  Smoothing parameters:
    alpha = 0.1056055 
    gamma = 0.2065322 

  Initial states:
     l[0]      s[0]     s[-1]     s[-2]    s[-3]    s[-4]    s[-5]    s[-6]
 4.516702 0.9991198 0.9980019 0.9996028 1.001731 1.002178 1.001732 1.000175
     s[-7]     s[-8]     s[-9]   s[-10]   s[-11]
 0.9980366 0.9985391 0.9991438 1.000211 1.001529

  sigma^2:  0

      AIC      AICc       BIC 
-1156.797 -1154.070 -1107.935 

Power_ETS_fit |>
  gg_tsresiduals()

The auto ETS model also a couple spikes that go past the blue dashed lines at lag 1 and slightly at lag 4.

Power_models <- power_ts_median_imputed |>
  model(
    auto_ARIMA = ARIMA(box_cox(KWH, Power_lambda)),
    auto_ETS = ETS(box_cox(KWH, Power_lambda))
  )

accuracy(Power_models) |>
  select(.model, RMSE, MAE) |>
  arrange(RMSE)

# A tibble: 2 × 3
  .model        RMSE     MAE
  <chr>        <dbl>   <dbl>
1 auto_ETS   651744. 490988.
2 auto_ARIMA 652389. 498499.

Comparing the Auto ARIMA and Auto ETS models, the Auto ETS model fits the historical data better hence the slightly lower RMSE and MAE values. Therefore I will use the Auto ETS model to generate the forecasts for 2014.

Power_ETS_fit |>
  forecast(h = 12) |>
  autoplot(power_ts_median_imputed) +
  labs(title = "Power KWH forecasts")

power_forecast <- Power_ETS_fit |>
  forecast(h = 12)

power_export <- power_forecast |>
  hilo(level = c(80, 95)) |>
  unpack_hilo(cols = c(`80%`, `95%`)) |>
  as_tibble() |>
  select(Month, .mean, `80%_lower`, `80%_upper`, `95%_lower`, `95%_upper`) |>
  rename(Forecast_KWH = .mean)

write_xlsx(power_export, "2014_Power_Forecast.xlsx")