Project 1

Part A: ATM Forecast

The objective of this analysis is to forecast daily cash withdrawals for four ATM machines during May 2010. The dataset contains historical daily cash withdrawal records from May 1, 2009, through April 30, 2010. Cash withdrawals are measured in hundreds of dollars. The analysis begins with data preparation and exploratory data analysis to identify trends, seasonal patterns, missing values, and unusual observations. Different forecasting methods will then be evaluated to determine which approaches provide reliable predictions for each ATM. The final forecasts will be used to estimate the expected cash requirements for May 2010 and provide information that may support ATM cash management and replenishment planning.

Import and data preparation

load_ATM_data <- read_excel("ATM624Data.xlsx") |>
  mutate(DATE = as.Date(DATE, origin = "1899-12-30"))

tsibble_ATM_data <- load_ATM_data |>
  filter(!is.na(ATM)) |>
  as_tsibble(key = ATM, index = DATE)

glimpse(tsibble_ATM_data)

## Rows: 1,460
## Columns: 3
## Key: ATM [4]
## $ DATE <date> 2009-05-01, 2009-05-02, 2009-05-03, 2009-05-04, 2009-05-05, 2009…
## $ ATM  <chr> "ATM1", "ATM1", "ATM1", "ATM1", "ATM1", "ATM1", "ATM1", "ATM1", "…
## $ Cash <dbl> 96, 82, 85, 90, 99, 88, 8, 104, 87, 93, 86, 111, 75, 6, 102, 73, …

range(tsibble_ATM_data$DATE)

## [1] "2009-05-01" "2010-04-30"

Interpretation:

After importing and preparing the dataset, we obtained 1,460 observations and three variables: DATE, ATM, and Cash. The dataset contains four ATM machines, each with 365 daily observations from May 1, 2009, through April 30, 2010. The DATE variable was converted from Excel’s numeric date format into R’s Date format. The dataset was then converted into a tsibble, using ATM as the key and DATE as the time index.

Exploratory data analysis

# Plot the four ATMs
tsibble_ATM_data |>
  autoplot(Cash) +
  facet_wrap(vars(ATM), scales = "free_y") +
  labs(
    title = "Daily ATM Cash Withdrawals",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

# Descriptive statistics
ATM_summary <- tsibble_ATM_data |>
  as_tibble() |>
  group_by(ATM) |>
  summarise(
    Start = min(DATE),
    End = max(DATE),
    Observations = n(),
    Missing = sum(is.na(Cash)),
    Zeros = sum(Cash == 0, na.rm = TRUE),
    Maximum = max(Cash, na.rm = TRUE),
    .groups = "drop"
  )

ATM_summary

## # A tibble: 4 × 7
##   ATM   Start      End        Observations Missing Zeros Maximum
##   <chr> <date>     <date>            <int>   <int> <int>   <dbl>
## 1 ATM1  2009-05-01 2010-04-30          365       3     0    180 
## 2 ATM2  2009-05-01 2010-04-30          365       2     2    147 
## 3 ATM3  2009-05-01 2010-04-30          365       0   362     96 
## 4 ATM4  2009-05-01 2010-04-30          365       0     0  10920.

# Check missing dates
has_gaps(tsibble_ATM_data)

## # A tibble: 4 × 2
##   ATM   .gaps
##   <chr> <lgl>
## 1 ATM1  FALSE
## 2 ATM2  FALSE
## 3 ATM3  FALSE
## 4 ATM4  FALSE

scan_gaps(tsibble_ATM_data)

## # A tsibble: 0 x 2 [?]
## # Key:       ATM [0]
## # ℹ 2 variables: ATM <chr>, DATE <date>

Interpretation:

The time-series plots reveal different withdrawal patterns across the four ATM machines. ATM1 and ATM2 exhibit frequent fluctuations, with recurring sharp decreases that may indicate weekly seasonality. ATM3 remains at zero for most of the observation period, followed by a sudden increase in withdrawals at the end of April 2010. ATM4 has considerably larger withdrawal amounts than the other ATMs and contains an extreme observation around March 2010. These differences suggest that a single forecasting method may not be appropriate for all four machines.

tsibble_ATM_data |>
  gg_season(Cash, period = "week")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

tsibble_ATM_data |>
  ACF(Cash) |>
  autoplot()

## Interpretation: The ACF plots suggest serial dependence and possible weekly seasonality for ATM1 and ATM2. ATM3 exhibits unusual behavior because most observations are zero, followed by a sudden increase. ATM4 shows weaker autocorrelation patterns, but the presence of an extreme observation makes it difficult to determine whether the series behaves like white noise.

Investigate missing values and outliers

# Missing Cash values
tsibble_ATM_data |>
  filter(is.na(Cash))

## # A tsibble: 5 x 3 [1D]
## # Key:       ATM [2]
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2009-06-13 ATM1     NA
## 2 2009-06-16 ATM1     NA
## 3 2009-06-22 ATM1     NA
## 4 2009-06-18 ATM2     NA
## 5 2009-06-24 ATM2     NA

# ATM3 nonzero withdrawals
tsibble_ATM_data |>
  filter(ATM == "ATM3", Cash > 0)

## # A tsibble: 3 x 3 [1D]
## # Key:       ATM [1]
##   DATE       ATM    Cash
##   <date>     <chr> <dbl>
## 1 2010-04-28 ATM3     96
## 2 2010-04-29 ATM3     82
## 3 2010-04-30 ATM3     85

# ATM4 largest withdrawals
tsibble_ATM_data |>
  filter(ATM == "ATM4") |>
  arrange(desc(Cash)) |>
  slice_head(n = 10)

## Warning: Current temporal ordering may yield unexpected results.
## ℹ Suggest to sort by `ATM`, `DATE` first.

## # A tsibble: 10 x 3 [1D]
## # Key:       ATM [1]
##    DATE       ATM     Cash
##    <date>     <chr>  <dbl>
##  1 2010-02-09 ATM4  10920.
##  2 2009-09-22 ATM4   1712.
##  3 2010-01-29 ATM4   1575.
##  4 2009-09-04 ATM4   1495.
##  5 2009-06-09 ATM4   1484.
##  6 2009-09-05 ATM4   1301.
##  7 2010-02-28 ATM4   1276.
##  8 2009-08-23 ATM4   1246.
##  9 2009-06-14 ATM4   1221.
## 10 2009-09-29 ATM4   1195.

Split the data into training data and test data

# Training period: May 2009 - March 2010
train_ATM <- tsibble_ATM_data |>
  filter(DATE < as.Date("2010-04-01"))

# Testing period: April 2010
test_ATM <- tsibble_ATM_data |>
  filter(
    DATE >= as.Date("2010-04-01"),
    DATE <= as.Date("2010-04-30")
  )

# Impute missing values in training data only
train_ATM <- train_ATM |>
  group_by(ATM) |>
  mutate(
    Cash = approx(
      x = which(!is.na(Cash)),
      y = Cash[!is.na(Cash)],
      xout = seq_along(Cash),
      rule = 2
    )$y
  ) |>
  ungroup()

# Check training and testing periods
range(train_ATM$DATE)

## [1] "2009-05-01" "2010-03-31"

range(test_ATM$DATE)

## [1] "2010-04-01" "2010-04-30"

Fit five forecasting models

fit_train <- train_ATM |>
  model(
    NAIVE = NAIVE(Cash),
    SNAIVE = SNAIVE(Cash ~ lag("week")),
    DRIFT = RW(Cash ~ drift()),
    ETS = ETS(Cash),
    ARIMA = ARIMA(Cash)
  )
fit_train

## # A mable: 4 x 6
## # Key:     ATM [4]
##   ATM     NAIVE   SNAIVE         DRIFT          ETS                    ARIMA
##   <chr> <model>  <model>       <model>      <model>                  <model>
## 1 ATM1  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,A)> <ARIMA(0,0,1)(0,1,2)[7]>
## 2 ATM2  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,A)> <ARIMA(2,0,2)(0,1,1)[7]>
## 3 ATM3  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,N)>   <ARIMA(0,0,0) w/ mean>
## 4 ATM4  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(M,N,A)>   <ARIMA(0,0,0) w/ mean>

# Forecast April 2010
fc_test <- fit_train |>
  forecast(new_data = test_ATM)
fc_test

## # A fable: 600 x 5 [1D]
## # Key:     ATM, .model [20]
##    ATM   .model DATE      
##    <chr> <chr>  <date>    
##  1 ATM1  NAIVE  2010-04-01
##  2 ATM1  NAIVE  2010-04-02
##  3 ATM1  NAIVE  2010-04-03
##  4 ATM1  NAIVE  2010-04-04
##  5 ATM1  NAIVE  2010-04-05
##  6 ATM1  NAIVE  2010-04-06
##  7 ATM1  NAIVE  2010-04-07
##  8 ATM1  NAIVE  2010-04-08
##  9 ATM1  NAIVE  2010-04-09
## 10 ATM1  NAIVE  2010-04-10
## # ℹ 590 more rows
## # ℹ 2 more variables: Cash <dist>, .mean <dbl>

Plot models for each ATM

# Convert forecasts to a regular tibble for plotting
forecast_plot_data <- fc_test |>
  as_tibble()

# Actual April observations
actual_plot_data <- test_ATM |>
  as_tibble()

ggplot() +
  geom_line(
    data = actual_plot_data,
    aes(x = DATE, y = Cash),
    color = "black",
    linewidth = 1.2
  ) +
  geom_line(
    data = forecast_plot_data,
    aes(
      x = DATE,
      y = .mean,
      color = .model
    ),
    linewidth = 0.8
  ) +
  facet_wrap(~ ATM, scales = "free_y", ncol = 2) +
  labs(
    title = "ATM Forecast Model Comparison - April 2010",
    subtitle = "Black = Actual Withdrawals",
    x = "Date",
    y = "Cash (Hundreds of Dollars)",
    color = "Model"
  ) +
  theme_minimal()

## Interpretation: Forecast Model Comparison

The five forecasting models were evaluated using the April 2010 test data. The results indicate that ETS provides the lowest RMSE for ATM1 and ATM2, suggesting that it produces more accurate forecasts than the other models during the testing period. Both ATMs exhibit recurring weekly fluctuations, which are captured more effectively by the seasonal forecasting methods than by the Naive and Random Walk with Drift models. For ATM3, the forecasting models struggle to predict the sudden increase in withdrawals toward the end of April because most of the training observations are zero. This suggests that the forecasting results should be interpreted cautiously. For ATM4, the Random Walk with Drift model achieves the lowest RMSE, although the forecasting errors remain relatively large. The extreme historical observation may also affect model performance. Overall, the results demonstrate that different ATM machines require different forecasting approaches because their withdrawal patterns are not identical.

Compare forecasting accuracy

accuracy_table <- accuracy(fc_test, test_ATM) |>
  select(ATM, .model, RMSE, MAE, MAPE)

# Rank models by RMSE for each ATM
accuracy_table |>
  as_tibble() |>
  group_by(ATM) |>
  arrange(RMSE, .by_group = TRUE)

## # A tibble: 20 × 5
## # Groups:   ATM [4]
##    ATM   .model  RMSE    MAE   MAPE
##    <chr> <chr>  <dbl>  <dbl>  <dbl>
##  1 ATM1  ETS     11.6   9.11   64.7
##  2 ATM1  ARIMA   12.4   9.71   86.4
##  3 ATM1  SNAIVE  16.0  13.5    73.4
##  4 ATM1  NAIVE   59.0  50.3   522. 
##  5 ATM1  DRIFT   60.3  51.8   529. 
##  6 ATM2  ETS     19.4  14.3    57.4
##  7 ATM2  ARIMA   21.0  16.3    79.5
##  8 ATM2  SNAIVE  26.2  17.8    45.1
##  9 ATM2  DRIFT   54.8  41.2  1140. 
## 10 ATM2  NAIVE   55.0  41.3  1143. 
## 11 ATM3  ARIMA   27.8   8.77  100  
## 12 ATM3  DRIFT   27.8   8.77  100  
## 13 ATM3  ETS     27.8   8.77  100  
## 14 ATM3  NAIVE   27.8   8.77  100  
## 15 ATM3  SNAIVE  27.8   8.77  100  
## 16 ATM4  DRIFT  291.  248.   1240. 
## 17 ATM4  ARIMA  293.  249.   1240. 
## 18 ATM4  NAIVE  297.  253.   1271. 
## 19 ATM4  SNAIVE 301.  227.    375. 
## 20 ATM4  ETS    386.  302.   1791.

Interpretation:

The best forecasting model for each ATM was selected primarily based on the lowest RMSE during the April 2010 testing period. MAE and MAPE were also examined as supplementary accuracy measures, although MAPE is less reliable when actual withdrawals are zero or close to zero.

Visualize model accuracy

accuracy_table |>
  ggplot(aes(
    x = .model,
    y = RMSE,
    fill = .model
  )) +
  geom_col() +
  facet_wrap(~ ATM, scales = "free_y") +
  labs(
    title = "Forecast Accuracy by ATM",
    x = "Forecasting Model",
    y = "RMSE"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

## Interpretation: We can Visualize those result base on the bar plot we have, it indicate ETS to be best model for ATM1, and ETS to be best for ATM2 and ATM3 all model might be good to work with, and for ATM4 the best model is RW/Drift.

best_models <- accuracy_table |>
  as_tibble() |>
  group_by(ATM) |>
  slice_min(RMSE, n = 1, with_ties = TRUE) |>
  ungroup()

best_models

## # A tibble: 8 × 5
##   ATM   .model  RMSE    MAE   MAPE
##   <chr> <chr>  <dbl>  <dbl>  <dbl>
## 1 ATM1  ETS     11.6   9.11   64.7
## 2 ATM2  ETS     19.4  14.3    57.4
## 3 ATM3  ARIMA   27.8   8.77  100  
## 4 ATM3  DRIFT   27.8   8.77  100  
## 5 ATM3  ETS     27.8   8.77  100  
## 6 ATM3  NAIVE   27.8   8.77  100  
## 7 ATM3  SNAIVE  27.8   8.77  100  
## 8 ATM4  DRIFT  291.  248.   1240.

Interpretation:

Base on the observation of the 4 ATMs the best model base on the RMSE, MAE and MAPE result we can conlude for ATM1 the best model would be ETS Model, and for ATM2 the best model would be ETS model, and for ATM3 the best model would be ARIMA model, and for ATM4 the model woukd be RW/Drift.

forecasting May 2010

# Prepare complete historical data
final_ATM_data <- tsibble_ATM_data |>
  group_by(ATM) |>
  mutate(
    Cash = approx(
      x = which(!is.na(Cash)),
      y = Cash[!is.na(Cash)],
      xout = seq_along(Cash),
      rule = 2
    )$y
  ) |>
  ungroup()
final_ATM_data

## # A tsibble: 1,460 x 3 [1D]
## # Key:       ATM [4]
##    DATE       ATM    Cash
##    <date>     <chr> <dbl>
##  1 2009-05-01 ATM1     96
##  2 2009-05-02 ATM1     82
##  3 2009-05-03 ATM1     85
##  4 2009-05-04 ATM1     90
##  5 2009-05-05 ATM1     99
##  6 2009-05-06 ATM1     88
##  7 2009-05-07 ATM1      8
##  8 2009-05-08 ATM1    104
##  9 2009-05-09 ATM1     87
## 10 2009-05-10 ATM1     93
## # ℹ 1,450 more rows

# Fit all five models using the full history
final_fit <- final_ATM_data |>
  model(
    NAIVE = NAIVE(Cash),
    SNAIVE = SNAIVE(Cash ~ lag("week")),
    DRIFT = RW(Cash ~ drift()),
    ETS = ETS(Cash),
    ARIMA = ARIMA(Cash)
  )
final_fit

## # A mable: 4 x 6
## # Key:     ATM [4]
##   ATM     NAIVE   SNAIVE         DRIFT          ETS                    ARIMA
##   <chr> <model>  <model>       <model>      <model>                  <model>
## 1 ATM1  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,A)> <ARIMA(0,0,1)(0,1,2)[7]>
## 2 ATM2  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,A)> <ARIMA(2,0,2)(0,1,1)[7]>
## 3 ATM3  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(A,N,N)>           <ARIMA(0,0,2)>
## 4 ATM4  <NAIVE> <SNAIVE> <RW w/ drift> <ETS(M,N,A)>   <ARIMA(0,0,0) w/ mean>

# Forecast May 2010
final_fc <- final_fit |>
  forecast(h = 31)
final_fc

## # A fable: 620 x 5 [1D]
## # Key:     ATM, .model [20]
##    ATM   .model DATE      
##    <chr> <chr>  <date>    
##  1 ATM1  NAIVE  2010-05-01
##  2 ATM1  NAIVE  2010-05-02
##  3 ATM1  NAIVE  2010-05-03
##  4 ATM1  NAIVE  2010-05-04
##  5 ATM1  NAIVE  2010-05-05
##  6 ATM1  NAIVE  2010-05-06
##  7 ATM1  NAIVE  2010-05-07
##  8 ATM1  NAIVE  2010-05-08
##  9 ATM1  NAIVE  2010-05-09
## 10 ATM1  NAIVE  2010-05-10
## # ℹ 610 more rows
## # ℹ 2 more variables: Cash <dist>, .mean <dbl>

Interpretation: May 2010 Forecasts

The May 2010 forecasts reveal different expected withdrawal patterns across the four ATM machines. ATM1 and ATM2 display recurring weekly fluctuations, suggesting that the seasonal patterns observed in the historical data may continue during the forecast period. ATM3 produces a relatively constant forecast, although the sudden increase in withdrawals near the end of April makes this prediction uncertain. For ATM4, the Random Walk with Drift model predicts a gradual decline in daily withdrawals throughout May 2010. However, the large historical outlier and forecasting errors suggest considerable uncertainty in these predictions. Overall, the results indicate that ATM1 and ATM2 exhibit more consistent weekly patterns, while ATM3 and ATM4 require additional investigation before their forecasts are used for cash replenishment decisions.

selected_models <- tibble(
  ATM = c("ATM1", "ATM2", "ATM3", "ATM4"),
  .model = c("ETS", "ETS", "ARIMA", "DRIFT")
)

# Keep only the selected forecasts
selected_fc <- final_fc |>
  as_tibble() |>
  inner_join(
    selected_models,
    by = c("ATM", ".model")
  )

# Plot May forecasts
ggplot(selected_fc,
       aes(x = DATE, y = .mean)) +
  geom_line(color = "steelblue", linewidth = 1) +
  facet_wrap(~ ATM, scales = "free_y") +
  labs(
    title = "Selected ATM Forecasts for May 2010",
    x = "Date",
    y = "Forecast Cash (Hundreds of Dollars)"
  ) +
  theme_minimal()

Interpretation:

The forecast plots indicate recurring weekly fluctuations for ATM1 and ATM2. ATM3 exhibits relatively stable predicted withdrawals as most its data show no withdrawals befor the sudden increase, while ATM4 shows a gradual downward trend. These patterns reflect the behavior of the selected forecasting models during May 2010.

# Historical data
historical_data <- final_ATM_data |>
  as_tibble()

# Selected May forecasts
forecast_data <- selected_fc

ggplot() +
  geom_line(
    data = historical_data,
    aes(x = DATE, y = Cash),
    color = "gray50",
    linewidth = 0.5
  ) +
  geom_line(
    data = forecast_data,
    aes(x = DATE, y = .mean),
    color = "blue",
    linewidth = 1
  ) +
  facet_wrap(~ ATM, scales = "free_y") +
  labs(
    title = "ATM Cash Withdrawal Forecasts for May 2010",
    subtitle = "Gray = Historical Data, Blue = Forecast",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  ) +
  theme_minimal()

selected_fc_dist <- final_fc |>
  semi_join(
    selected_models,
    by = c("ATM", ".model")
  )

selected_fc_dist |>
  autoplot(final_ATM_data, level = c(80, 95)) +
  facet_wrap(vars(ATM), scales = "free_y") +
  labs(
    title = "ATM Cash Withdrawal Forecasts for May 2010",
    x = "Date",
    y = "Cash (Hundreds of Dollars)"
  )

Interpretation: May 2010 Forecasts

The selected forecasting models produce different patterns for the four ATM machines during May 2010. ATM1 and ATM2 exhibit recurring weekly fluctuations, with relatively high withdrawals on certain days and sharp decreases on others. These patterns are consistent with the seasonal behavior observed in the historical data.For ATM3, the selected ARIMA model produces forecasts based on the historical withdrawal patterns. However, because ATM3 recorded zero withdrawals for most of the historical period before experiencing a sudden increase near the end of April, these forecasts should be interpreted cautiously. For ATM4, the Random Walk with Drift model predicts a gradual decrease in daily cash withdrawals throughout May. Nevertheless, the extreme observation in the historical data and the relatively high forecasting errors indicate substantial uncertainty. Overall, the forecasts suggest that ATM1 and ATM2 have more predictable weekly withdrawal patterns, while ATM3 and ATM4 require greater caution when estimating future cash requirements.

Final Stage

#install.packages("writexl")
library(writexl)

## Warning: package 'writexl' was built under R version 4.5.2

# Total expected withdrawals for May
monthly_forecast <- selected_fc |>
  group_by(ATM) |>
  summarise(
    May_Total_USD = sum(.mean) * 100,
    .groups = "drop"
  )

monthly_forecast

## # A tibble: 4 × 2
##   ATM   May_Total_USD
##   <chr>         <dbl>
## 1 ATM1        238361.
## 2 ATM2        182757.
## 3 ATM3           402.
## 4 ATM4       1454932.

# Export daily forecasts
daily_forecast <- selected_fc |>
  transmute(
    ATM,
    DATE,
    Model = .model,
    Forecast_USD = .mean * 100
  )

write_xlsx(
  list(
    Daily_Forecasts = daily_forecast,
    Monthly_Totals = monthly_forecast
  ),
  "ATM_May_2010_Forecasts.xlsx"
)

Conclusion For Part A:

The objective of this project was to forecast daily cash withdrawals for four ATM machines during May 2010 using historical observations from May 2009 through April 2010. Five forecasting methods were evaluated using April 2010 as the testing period: Naive, Seasonal Naive, Random Walk with Drift, ETS, and ARIMA. The results indicate that ETS provides the best forecasting accuracy for ATM1 and ATM2 based on RMSE, while Random Walk with Drift performs best for ATM4. ATM3 presents a particular forecasting challenge because withdrawals remained at zero for most of the historical period before increasing sharply near the end of April. The analysis demonstrates that forecasting performance depends on the characteristics of each ATM. Weekly seasonality, unusual observations, and changes in withdrawal activity can affect forecasting accuracy. Although the selected models provide estimates of future cash requirements, prediction uncertainty should be considered when planning ATM cash replenishment. Further improvements could include investigating unusual observations, incorporating additional historical data, and evaluating the forecasting models over multiple testing periods.

Part B: Forcating Power, Residential Customer Forecast Load

Load the data

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

## # A tibble: 192 × 3
##    CaseSequence `YYYY-MMM`     KWH
##           <dbl> <chr>        <dbl>
##  1          733 1998-Jan   6862583
##  2          734 1998-Feb   5838198
##  3          735 1998-Mar   5420658
##  4          736 1998-Apr   5010364
##  5          737 1998-May   4665377
##  6          738 1998-Jun   6467147
##  7          739 1998-Jul   8914755
##  8          740 1998-Aug   8607428
##  9          741 1998-Sep   6989888
## 10          742 1998-Oct   6345620
## # ℹ 182 more rows

Understanding the data

glimpse(load_Power_data)

## Rows: 192
## Columns: 3
## $ CaseSequence <dbl> 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 74…
## $ `YYYY-MMM`   <chr> "1998-Jan", "1998-Feb", "1998-Mar", "1998-Apr", "1998-May…
## $ KWH          <dbl> 6862583, 5838198, 5420658, 5010364, 4665377, 6467147, 891…

summary(load_Power_data)

##   CaseSequence     YYYY-MMM              KWH          
##  Min.   :733.0   Length:192         Min.   :  770523  
##  1st Qu.:780.8   Class :character   1st Qu.: 5429912  
##  Median :828.5   Mode  :character   Median : 6283324  
##  Mean   :828.5                      Mean   : 6502475  
##  3rd Qu.:876.2                      3rd Qu.: 7620524  
##  Max.   :924.0                      Max.   :10655730  
##                                     NA's   :1

colSums(is.na(load_Power_data))

## CaseSequence     YYYY-MMM          KWH 
##            0            0            1

sum(duplicated(load_Power_data))

## [1] 0

Preparing the data

is_tsibble(load_Power_data)

## [1] FALSE

# Convert to tsible
power_data <- load_Power_data |>
  mutate(Date = yearmonth(`YYYY-MMM`)) |>
  select(Date, KWH) |>
  as_tsibble(
    index = Date
  )
power_data

## # A tsibble: 192 x 2 [1M]
##        Date     KWH
##       <mth>   <dbl>
##  1 1998 Jan 6862583
##  2 1998 Feb 5838198
##  3 1998 Mar 5420658
##  4 1998 Apr 5010364
##  5 1998 May 4665377
##  6 1998 Jun 6467147
##  7 1998 Jul 8914755
##  8 1998 Aug 8607428
##  9 1998 Sep 6989888
## 10 1998 Oct 6345620
## # ℹ 182 more rows

is_tsibble(power_data)

## [1] TRUE

has_gaps(power_data)

## # A tibble: 1 × 1
##   .gaps
##   <lgl>
## 1 FALSE

scan_gaps(power_data)

## # A tsibble: 0 x 1 [?]
## # ℹ 1 variable: Date <mth>

count_gaps(power_data)

## # A tibble: 0 × 3
## # ℹ 3 variables: .from <mth>, .to <mth>, .n <int>

range(power_data$Date)

## <yearmonth[2]>
## [1] "1998 Jan" "2013 Dec"

power_data <- power_data |>
  fill_gaps(KWH = NA_real_)

has_gaps(power_data)

## # A tibble: 1 × 1
##   .gaps
##   <lgl>
## 1 FALSE

Interpretation: Data Preparation

The residential electricity consumption dataset contains 192 monthly observations, with one missing value in the KWH variable. The YYYY-MMM column was converted into a monthly time index using yearmonth(), and the dataset was transformed into a tsibble for time series analysis. The dataset was also examined for missing observations and gaps in the monthly sequence. Any missing time periods were made explicit using fill_gaps(). The missing electricity consumption values will be investigated and treated appropriately before fitting forecasting models.

Exploratory Analysis

power_data |>
  autoplot(KWH) + labs(title = "Monthly Power Customer Consomption in KWH",
                       x = "Date",
                       y = "Electricity Consumption (KWH)")

# Seasonal patterns
power_data |>
  gg_season(KWH) + labs(title = "Monthly Electricity Consumption by Season",
                        y = "Electricity Consumption (KWH)",
                        x = "Month")

## Interpretation: Exploratory Analysis

The time series plot reveals a recurring seasonal pattern in monthly residential electricity consumption. Electricity demand generally increases during winter and summer and decreases during spring and autumn. There also appears to be some increase in the magnitude of peak consumption toward the end of the observation period, although the overall trend is not consistently upward.

The seasonal plot provides further evidence of annual seasonality. Electricity consumption is generally higher during January and the summer months, particularly July and August, while lower consumption is observed during March, April, May, October, and November. An unusually low observation appears around July 2010, which may represent an outlier or an unusual event that requires further investigation.

power_data |>
  filter(
    Date >= yearmonth("2010 Jun"),
    Date <= yearmonth("2010 Aug")
  )

## # A tsibble: 3 x 2 [1M]
##       Date     KWH
##      <mth>   <dbl>
## 1 2010 Jun 6696292
## 2 2010 Jul  770523
## 3 2010 Aug 7922701

power_data |>
  ACF(KWH) |> autoplot()

## Interpretation: Unusual Observation in July 2010

Further investigation of the residential electricity consumption dataset identified an unusually low observation in July 2010, when consumption was recorded as 770,523 KWH. In comparison, electricity consumption was 6,696,292 KWH in June 2010 and 7,922,701 KWH in August 2010. The July observation was approximately 89.5% lower than the average of the two neighboring months, representing a substantial deviation from the surrounding consumption levels.

This unusual observation may indicate a data-entry error, incomplete reporting, or an actual disruption in electricity consumption. However, without additional information about the data collection process, the cause cannot be determined. Therefore, the observation was retained in the original analysis to preserve the integrity of the dataset. Its potential influence on model estimation and forecasting accuracy should be considered when interpreting the results.

Interpretation: Autocorrelation Analysis

The ACF plot reveals a strong seasonal pattern in monthly electricity consumption. Significant positive autocorrelations occur at several lags, particularly lag 12, indicating that electricity consumption in a given month is strongly associated with consumption during the same month in the previous year. Negative autocorrelations at several intermediate lags suggest differences in electricity demand between seasons. These results provide additional evidence of annual seasonality and suggest that seasonal forecasting models, such as Seasonal ARIMA or seasonal ETS, may be appropriate.

Transformation

# Using box-cox, log transformation
lambda <- power_data |>
  features(KWH, features = guerrero) |>
  pull(lambda_guerrero)
lambda

## [1] 0.1132906

# Box_cox transformation
box_cox_power <- power_data |>
  mutate(box_cox_pw = box_cox(KWH, lambda = lambda))
box_cox_power

## # A tsibble: 192 x 3 [1M]
##        Date     KWH box_cox_pw
##       <mth>   <dbl>      <dbl>
##  1 1998 Jan 6862583       43.7
##  2 1998 Feb 5838198       42.7
##  3 1998 Mar 5420658       42.3
##  4 1998 Apr 5010364       41.9
##  5 1998 May 4665377       41.4
##  6 1998 Jun 6467147       43.3
##  7 1998 Jul 8914755       45.3
##  8 1998 Aug 8607428       45.1
##  9 1998 Sep 6989888       43.8
## 10 1998 Oct 6345620       43.2
## # ℹ 182 more rows

# Visualize the box-cox transformation

box_cox_power |>
  autoplot(box_cox_pw) + labs(title = "Box Cox Transformation of Power Customer Data")

# Log Transformation
log_power_tr <- power_data |>
  mutate(log_trans = log(KWH))

# Visualize

log_power_tr |> autoplot(log_trans) + labs(title = "Log Transformation of Power Customer Data")

Interpretation:

Interpretation: Data Transformation

A Box-Cox transformation was applied to investigate whether stabilizing the variance could improve the modeling of residential electricity consumption. A logarithmic transformation was also examined for comparison. Both transformations changed the scale of the observations while preserving the underlying seasonal structure.

The transformed series continued to exhibit recurring seasonal fluctuations. Therefore, transformation alone did not eliminate seasonality, and forecasting models capable of handling seasonal patterns were still required. The effectiveness of the Box-Cox transformation was subsequently evaluated through out-of-sample forecast accuracy.

Train/test split

# Training period: yearMonth 2014
train_Power <- power_data |>
  filter(Date <= yearmonth("2011 Dec"))

# Testing period: 2012 - 2013
test_Power <- power_data |>
  filter(
    Date >= yearmonth("2012 Jan"),
    Date <= yearmonth("2013 Dec")
  )

# Impute missing values in training data only
train_Power_data <- train_Power |>
  arrange(Date) |>
  mutate(
    KWH = approx(
      x = which(!is.na(KWH)),
      y = KWH[!is.na(KWH)],
      xout = seq_along(KWH),
      rule = 2
    )$y
  ) |>
  ungroup()

# Check training and testing periods
range(train_Power_data$Date)

## <yearmonth[2]>
## [1] "1998 Jan" "2011 Dec"

range(test_Power$Date)

## <yearmonth[2]>
## [1] "2012 Jan" "2013 Dec"

# Check missing values
sum(is.na(train_Power$KWH))

## [1] 1

sum(is.na(train_Power_data$KWH))

## [1] 0

# Visualize the training data 

train_Power_data |>
  autoplot(KWH) + labs(title = "Residential Electricity Consumption After Imputation")

## Interpretation: To better have a training data that can capture the dataset we used training data from 1998- 2011 and for testing the data we used 2012-2013 that will gives us better training and we can then forecast with accuraccy.

Model fitting

 # Box-Cox lambda from training data
lambda_train <- train_Power_data |>
  features(KWH, guerrero) |>
  pull(lambda_guerrero)

# Fit models
power_models <- train_Power_data |>
  model(SNAIVE = SNAIVE(KWH),
        ETS = ETS(KWH),
        ARIMA = ARIMA(KWH),
        ARIMA_BoxCox = ARIMA(box_cox(KWH, lambda_train)))

glance(power_models)

## # A tibble: 4 × 11
##   .model     sigma2 log_lik   AIC  AICc   BIC      MSE     AMSE     MAE ar_roots
##   <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>   <dbl> <list>  
## 1 SNAIVE   1.46e+12     NA    NA    NA    NA  NA       NA       NA      <NULL>  
## 2 ETS      1.45e- 2  -2696. 5421. 5424. 5468.  6.99e11  7.55e11  0.0774 <NULL>  
## 3 ARIMA    6.82e+11  -2352. 4714. 4715. 4729. NA       NA       NA      <cpl>   
## 4 ARIMA_B… 5.29e- 3    195. -380. -380. -365. NA       NA       NA      <cpl>   
## # ℹ 1 more variable: ma_roots <list>

# Forecast 24 months
power_test_forecasts <- power_models |>
  forecast(h = 24)

# Compare models
power_accuracy <- power_test_forecasts |>
  accuracy(test_Power)

power_accuracy |>
  select(.model, RMSE, MAE, MAPE) |>
  arrange(RMSE)

## # A tibble: 4 × 4
##   .model           RMSE      MAE  MAPE
##   <chr>           <dbl>    <dbl> <dbl>
## 1 ETS           861658.  646527.  8.56
## 2 SNAIVE        965567.  742063.  9.30
## 3 ARIMA        1317472.  843866. 10.4 
## 4 ARIMA_BoxCox 1982572. 1066206. 12.8

Interpretation:

Four forecasting models were fitted to the residential electricity consumption data: Seasonal Naïve, ETS, ARIMA, and ARIMA with a Box-Cox transformation. The models were evaluated using the 2012–2013 testing period. The ETS model achieved the lowest RMSE of approximately 861,658 KWH, MAE of 646,527 KWH, and MAPE of 8.56%, outperforming the other three models. The Seasonal Naïve model was the second-best performer, with an RMSE of approximately 965,567 KWH. ARIMA and ARIMA_BoxCox produced higher forecasting errors, with ARIMA_BoxCox performing the worst. These results suggest that ETS is the most suitable model among the four candidates for forecasting residential electricity consumption in 2014.

Model diagnostics

# Residual ACF, Ljung - Box
power_residuals <- power_models |>
  augment()
power_residuals

## # A tsibble: 672 x 6 [1M]
## # Key:       .model [4]
##    .model     Date     KWH .fitted .resid .innov
##    <chr>     <mth>   <dbl>   <dbl>  <dbl>  <dbl>
##  1 SNAIVE 1998 Jan 6862583      NA     NA     NA
##  2 SNAIVE 1998 Feb 5838198      NA     NA     NA
##  3 SNAIVE 1998 Mar 5420658      NA     NA     NA
##  4 SNAIVE 1998 Apr 5010364      NA     NA     NA
##  5 SNAIVE 1998 May 4665377      NA     NA     NA
##  6 SNAIVE 1998 Jun 6467147      NA     NA     NA
##  7 SNAIVE 1998 Jul 8914755      NA     NA     NA
##  8 SNAIVE 1998 Aug 8607428      NA     NA     NA
##  9 SNAIVE 1998 Sep 6989888      NA     NA     NA
## 10 SNAIVE 1998 Oct 6345620      NA     NA     NA
## # ℹ 662 more rows

power_models |>
  select(ARIMA) |>
  gg_tsresiduals()

power_models |>
  select(ETS) |>
  gg_tsresiduals()

power_models |>
  select(ARIMA_BoxCox) |>
  gg_tsresiduals()

power_models |>
  select(SNAIVE) |>
  gg_tsresiduals()

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

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

## Warning: Removed 12 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_rug()`).

# Plot ACF of All four models
power_residuals |>
  ACF(.innov, lag_max = 36) |>
  autoplot() +
  facet_wrap(vars(.model), scales = "free_y") +
  labs(
    title = "Residual ACF Comparison",
    x = "Lag",
    y = "ACF"
  )

power_ljungbox <- power_models |>
  augment() |>
  features(
    .innov,
    ljung_box,
    lag = 24,
    dof = 0
  )

power_ljungbox

## # A tibble: 4 × 3
##   .model       lb_stat lb_pvalue
##   <chr>          <dbl>     <dbl>
## 1 ARIMA          17.8  0.814    
## 2 ARIMA_BoxCox    7.42 1.000    
## 3 ETS            34.0  0.0841   
## 4 SNAIVE         63.5  0.0000204

Interpretation:

Residual diagnostics were performed using residual time plots, ACF plots, and the Ljung-Box test. The residual ACF plots show that the ARIMA and ARIMA_BoxCox models have relatively weak remaining autocorrelation, while the Seasonal Naïve model exhibits significant autocorrelation at certain lags. Using 24 lags and no degrees-of-freedom adjustment, the Ljung-Box test produced p-values of 0.8140 for ARIMA, 0.9995 for ARIMA_BoxCox, and 0.0841 for ETS. These results do not provide sufficient evidence to reject the null hypothesis of uncorrelated residuals at the 5% significance level. In contrast, the Seasonal Naïve model produced a p-value of approximately 0.00002, indicating significant remaining autocorrelation. Although ETS does not have the highest Ljung-Box p-value, it achieved the lowest forecasting errors on the testing dataset. The results should be interpreted as preliminary until model-specific degrees-of-freedom adjustments are considered. # Forecast evaluation

# Comparing Diagnostics results
model_comparison <- power_accuracy |>
  select(.model, RMSE, MAE, MAPE) |>
  left_join(
    power_ljungbox,
    by = ".model"
  ) |>
  arrange(RMSE)

model_comparison

## # A tibble: 4 × 6
##   .model           RMSE      MAE  MAPE lb_stat lb_pvalue
##   <chr>           <dbl>    <dbl> <dbl>   <dbl>     <dbl>
## 1 ETS           861658.  646527.  8.56   34.0  0.0841   
## 2 SNAIVE        965567.  742063.  9.30   63.5  0.0000204
## 3 ARIMA        1317472.  843866. 10.4    17.8  0.814    
## 4 ARIMA_BoxCox 1982572. 1066206. 12.8     7.42 1.000

Interpretation:

Interpretation: Forecast Evaluation

The four forecasting models were compared using RMSE, MAE, and MAPE calculated from the 2012–2013 testing dataset. The ETS model achieved the lowest RMSE of approximately 861,658 KWH and MAPE of 8.56%, followed by the Seasonal Naïve model with an RMSE of approximately 965,567 KWH and MAPE of 9.30%. The ARIMA and ARIMA_BoxCox models produced higher forecasting errors, with the transformed ARIMA model performing the worst.

Although the ARIMA and ARIMA_BoxCox models showed less evidence of residual autocorrelation based on the Ljung–Box test, ETS achieved the best out-of-sample forecasting accuracy. Therefore, ETS was selected as the preferred forecasting model for residential electricity consumption in 2014. # Final Forecast

# Impute missing values in the full historical data
power_data_final <- power_data |>
  arrange(Date) |>
  mutate(
    KWH = approx(
      x = which(!is.na(KWH)),
      y = KWH[!is.na(KWH)],
      xout = seq_along(KWH),
      rule = 2
    )$y
  )

# Verify
sum(is.na(power_data_final$KWH))

## [1] 0

# Fit the selected model
final_model <- power_data_final |>
  model(ETS = ETS(KWH))

report(final_model)

## Series: KWH 
## Model: ETS(M,N,M) 
##   Smoothing parameters:
##     alpha = 0.1137195 
##     gamma = 0.0001001086 
## 
##   Initial states:
##     l[0]      s[0]     s[-1]     s[-2]    s[-3]    s[-4]    s[-5]     s[-6]
##  6134773 0.9547463 0.7533319 0.8602685 1.191176 1.255577 1.199206 0.9982427
##      s[-7]     s[-8]     s[-9]   s[-10]   s[-11]
##  0.7683657 0.8160779 0.9174797 1.061846 1.223683
## 
##   sigma^2:  0.0144
## 
##      AIC     AICc      BIC 
## 6224.057 6226.784 6272.919

# Forecast 12 months of 2014
forecast_2014 <- final_model |>
  forecast(h = 12)

forecast_2014 |>
  autoplot(power_data_final) +
  labs(
    title = "Residential Electricity Consumption Forecast for 2014",
    x = "Year",
    y = "KWH"
  )

# Display 2014 monthly point forecasts
forecast_2014 |>
  as_tibble() |>
  select(Date, .mean)

## # A tibble: 12 × 2
##        Date    .mean
##       <mth>    <dbl>
##  1 2014 Jan 9304198.
##  2 2014 Feb 8073653.
##  3 2014 Mar 6975902.
##  4 2014 Apr 6204893.
##  5 2014 May 5842200.
##  6 2014 Jun 7589986.
##  7 2014 Jul 9117199.
##  8 2014 Aug 9546786.
##  9 2014 Sep 9056676.
## 10 2014 Oct 6541029.
## 11 2014 Nov 5727865.
## 12 2014 Dec 7259239.

# Display 80% and 95% prediction intervals
forecast_2014 |>
  hilo(level = c(80, 95))

## # A tsibble: 12 x 6 [1M]
## # Key:       .model [1]
##    .model     Date
##    <chr>     <mth>
##  1 ETS    2014 Jan
##  2 ETS    2014 Feb
##  3 ETS    2014 Mar
##  4 ETS    2014 Apr
##  5 ETS    2014 May
##  6 ETS    2014 Jun
##  7 ETS    2014 Jul
##  8 ETS    2014 Aug
##  9 ETS    2014 Sep
## 10 ETS    2014 Oct
## 11 ETS    2014 Nov
## 12 ETS    2014 Dec
## # ℹ 4 more variables: KWH <dist>, .mean <dbl>, `80%` <hilo>, `95%` <hilo>

Interpretation: Final Forecast for 2014

The ETS model was selected to forecast monthly residential electricity consumption for January through December 2014 because it achieved the lowest forecasting errors during the testing period. After refitting the model using the complete historical dataset from January 1998 through December 2013, the selected specification was ETS(M,N,M), representing multiplicative errors, no trend, and multiplicative seasonality.

The forecast exhibits recurring seasonal fluctuations, with electricity consumption expected to increase during periods of higher seasonal demand and decrease during periods of lower demand. The absence of a trend component suggests that the underlying level of electricity consumption is expected to remain relatively stable throughout 2014, while seasonal variations continue.

The estimated level smoothing parameter is approximately 0.1137, indicating gradual adjustments to the consumption level. The seasonal smoothing parameter is approximately 0.0001, suggesting that the seasonal pattern remains relatively stable over time.

The 80% and 95% prediction intervals represent the uncertainty surrounding the monthly forecasts. The wider 95% intervals provide a broader range of plausible future consumption values.

Overall, the ETS(M,N,M) model captures the recurring seasonal behavior observed in the historical data and provides monthly electricity consumption forecasts for 2014.

Final stage

# Final stage: Export 2014 forecasts
library(dplyr)
library(writexl)

# Extract monthly forecasts
monthly_forecast <- forecast_2014 |>
  mutate(
    Lower_80 = quantile(KWH, 0.10),
    Upper_80 = quantile(KWH, 0.90),
    Lower_95 = quantile(KWH, 0.025),
    Upper_95 = quantile(KWH, 0.975)
  ) |>
  as_tibble() |>
  transmute(
    Date = as.character(Date),
    Model = .model,
    Forecast_KWH = .mean,
    Lower_80,
    Upper_80,
    Lower_95,
    Upper_95
  )

# Calculate total forecast consumption for 2014
yearly_forecast <- monthly_forecast |>
  summarise(
    Year = 2014,
    Total_Forecast_KWH = sum(Forecast_KWH),
    Average_Monthly_KWH = mean(Forecast_KWH),
    .groups = "drop"
  )

yearly_forecast

## # A tibble: 1 × 3
##    Year Total_Forecast_KWH Average_Monthly_KWH
##   <dbl>              <dbl>               <dbl>
## 1  2014          91239626.            7603302.

# Export to Excel
write_xlsx(
  list(
    Monthly_Forecasts = monthly_forecast,
    Annual_Summary = yearly_forecast
  ),
  "KWH_2014_Forecasts.xlsx"
)

Conclusion: Part B

The objective of Part B was to analyze historical residential electricity consumption and develop monthly forecasts for 2014. The analysis included data preparation, exploratory analysis, transformations, model fitting, diagnostic testing, and forecast evaluation.

The exploratory analysis revealed a strong annual seasonal pattern in electricity consumption, with higher demand during winter and summer and lower consumption during other months. An unusually low observation was identified in July 2010 and retained because its cause could not be verified.

Four forecasting models were evaluated using the 2012–2013 testing period: Seasonal Naïve, ETS, ARIMA, and ARIMA with a Box-Cox transformation. Among these models, ETS achieved the lowest RMSE of approximately 861,658 KWH and MAPE of 8.56%, indicating the strongest out-of-sample forecasting performance.

The ETS model was subsequently refitted using the complete historical dataset. The selected ETS(M,N,M) specification captured the recurring seasonal fluctuations without including a trend component. The resulting forecasts indicate that residential electricity consumption is expected to continue exhibiting seasonal variation throughout 2014.

Although the ETS model provided the best forecasting accuracy among the models considered, the unusual July 2010 observation and uncertainty in future electricity demand remain important limitations. Overall, the analysis demonstrates the importance of seasonal modeling, out-of-sample evaluation, and residual diagnostics when forecasting residential electricity consumption.