We will first install the packages and load the libraries that would be usefull for our project. I will install the packages in the console so that they are not being reinstall every time that i run my file.
library(fpp3)
Warning: package 'fpp3' was built under R version 4.5.3
── 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(lubridate)library(tseries)
Warning: package 'tseries' was built under R version 4.5.3
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
library(openxlsx)
Warning: package 'openxlsx' was built under R version 4.5.3
library(writexl)
Warning: package 'writexl' was built under R version 4.5.3
2) IMPORT DATA
Let’s import our ATM excel dataset, check the structure of the dataset for a better understanding of the data types we work to facilitate the cleaning of the data and a general descriptive statistics to have an overall initial statistic of our data.
# Load the raw dataatm_raw <-read_excel("ATM624Data.xlsx")head(atm_raw)
DATE ATM Cash
Min. :2009-05-01 Length:1474 Min. : 0.0
1st Qu.:2009-08-01 Class :character 1st Qu.: 0.5
Median :2009-11-01 Mode :character Median : 73.0
Mean :2009-10-31 Mean : 155.6
3rd Qu.:2010-02-01 3rd Qu.: 114.0
Max. :2010-05-14 Max. :10919.8
NA's :19
The starting date in our dataset is 2009-05-01 and it ends in 2010-05-14. The cash columns contains 19 missing dates which we should further when cleaning our data.
Let’s verify, if those missing data are in every individual columns.To perform this operation, we will consider the ATM column as our reference.
The output confirms that all missing are part of the Cash column with 3 missing cash value from ATM1, 2 missing cash value from ATM2 and 14 from the ATM categories.
Let’s remove the rows with missing ATM values, convert our dataframe into a time series with DATE as our index and plot the original data to observe the pattern of each key variable
# removing 14 rows atm_raw_clean <- atm_raw |>filter(!is.na(ATM))# conversion to time series.atm_raw_clean_ts <- atm_raw_clean |>as_tsibble(index=DATE, key= ATM)# Plotatm_raw_clean_ts |>autoplot(Cash) +labs(title ="CASH WITHDRAWAL IN ATMs over time",y="Cash in hundreds of dollars")
Interpreting
This plot exhibit that cash withdrawal from ATM1, ATM2, ATM3 fall mostly withim the same range whereas cash withdrawal from ATM4 is in a higher range with a huge outlier in 2010.
We will perform a case by case study meaning forecasting cash withdrawal amount in each ATM.
3) Splitting each ATM time series
We will filter each individual ATM in the ATM column from the overall clean time series.
This individual plot shows the ATM3 cash withdrawal amount is majority equal to zero except from certain case. Let’s check the non zero value.
# number rows in the ATM3 time seriesnrow(atm3_ts)
[1] 365
# View all non-zero rows in ATM3atm3_non_zero <- atm3_ts |>filter(Cash >0)atm3_non_zero
# 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
# Calculate the average cash withdrawalmean(atm3_non_zero$Cash)
[1] 87.66667
There are three non zeros value out of the 365 rows in the ATM3 time series.
4) Testing and Evaluation of Models
# ATM1: model test Candidate Comparison## Convert to a clean univariate tsts_atm1 <- atm_raw %>%filter(ATM =="ATM1") %>%arrange(DATE) %>%pull(Cash) %>%# extracts the column as a vector, not a data framets(frequency =7)fit_atm1_ets <-ets(ts_atm1)fit_atm1_arima <-auto.arima(ts_atm1, seasonal =TRUE, stepwise =FALSE, approximation =FALSE)fc_atm1 <-forecast(fit_atm1_arima, h =31) # h=31 because each day represents 1 observation and May has 31 days # ATM2: model test Candidate Comparison ## Convert to a clean univariate tsts_atm2 <- atm_raw %>%filter(ATM =="ATM2") %>%arrange(DATE) %>%pull(Cash) %>%# extracts the column as a vector, not a data framets(frequency =7)fit_atm2_ets <-ets(ts_atm2)fit_atm2_arima <-auto.arima(ts_atm2, seasonal =TRUE, stepwise =FALSE, approximation =FALSE)fc_atm2 <-forecast(fit_atm2_arima, h =31)# ATM3: Structural Shift Modeling because of too many zeros valuesmean_active_atm3 <-mean(atm3_non_zero$Cash)# Active baseline forecast for 31 days of Mayfc_atm3_point <-rep(mean_active_atm3, 31)# ATM4:model test Candidate Comparison## Convert to a clean univariate tsts_atm4 <- atm_raw |>filter(ATM =="ATM4") |>arrange(DATE) |>pull(Cash) |># extracts the column as a vector, not a data framets(frequency =7)fit_atm4_ets <-ets(ts_atm4)fit_atm4_arima <-auto.arima(ts_atm4, seasonal =FALSE, stepwise =FALSE)fc_atm4 <-forecast(fit_atm4_arima, h =31)
5) Assemble Forecast output Table for the month of May 2010
write.xlsx(forecast_results, "ATM_Cash_Withdrawal_May_2010_Forecast.xlsx", overwrite =TRUE)print("Forecast ATM Cash withdrawal of the month of May completed successfully.")
[1] "Forecast ATM Cash withdrawal of the month of May completed successfully."
PART B: Forecasting Power, ResidentialCustomerForecastLoad-624.xlsx
1) Loading the data and checking the characteristics
# load data resid_power<-read_excel("ResidentialCustomerForecastLoad-624.xlsx")head(resid_power)
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
The dataset is 192 rows and 3 columns dataset and the summary statistic shows that there is one missing residential power usage value. It is a monthly frequenced dataset.
2) Let’s convert the data to the time series and plot it
# The frequency is monthly (12) and the series starts in Jan 1998ts_resid_power <-ts(resid_power$KWH, start =c(1998, 1), frequency =12)ts_resid_power_clean <-tsclean(ts_resid_power)head(ts_resid_power_clean)
Jan Feb Mar Apr May Jun
1998 6862583 5838198 5420658 5010364 4665377 6467147
tail(ts_resid_power_clean)
Jul Aug Sep Oct Nov Dec
2013 8415321 9080226 7968220 5759367 5769083 6987874
3) Testing and Evaluating models
# Let's test and evaluate Exponential Smoothing / Holt-Winters Modelsfit_ses <-ses(ts_resid_power_clean, h =12)fit_holt <-holt(ts_resid_power_clean, h =12)fit_ets_m <-ets(ts_resid_power_clean, model ="MAM")fit_ets_a <-ets(ts_resid_power_clean, model ="AAA")fit_ets_auto <-ets(ts_resid_power_clean)# Let's test and evaluate Seasonal ARIMA Modelinglambda_PU <-BoxCox.lambda(ts_resid_power_clean)cat("Optimal Box-Cox Lambda:", lambda_PU, "\n")
Optimal Box-Cox Lambda: -0.1454175
# let's fit auto.arima with and without Box-Cox transformationfit_arima_raw <-auto.arima(ts_resid_power_clean, stepwise =FALSE, approximation =FALSE)fit_arima_bc <-auto.arima(ts_resid_power_clean, lambda = lambda_PU, stepwise =FALSE, approximation =FALSE)
4} Comparative models evaluation
# Let's check the AICc parameters from the modelscat("\n--- Model Comparison Summary ---\n")
--- Model Comparison Summary ---
cat("ETS Model Selected:", fit_ets_auto$method, "AICc:", fit_ets_auto$aicc, "\n")
ME RMSE MAE MPE MAPE MASE ACF1
Training set -7298.943 583970 444014.7 -0.7599925 6.764759 0.7453874 0.3058204
cat("\nIn-Sample Accuracy - ARIMA:\n")
In-Sample Accuracy - ARIMA:
print(forecast::accuracy(fit_arima_bc))
ME RMSE MAE MPE MAPE MASE
Training set 27171.71 566612.1 433555.9 -0.08069716 6.637629 0.7278298
ACF1
Training set 0.05365242
# Let's check residuals diagnostics to choose the best modelcheckresiduals(fit_arima_bc)
Ljung-Box test
data: Residuals from ARIMA(0,0,3)(2,1,0)[12] with drift
Q* = 11.874, df = 19, p-value = 0.8909
Model df: 5. Total lags used: 24
Interpretation
The model residuals exhibit a bell shape curve that is centered around zero with a constant variance across time which confirms that the seasonal differencing and Box-Cox transformation successfully stabilized the trend and heteroskedastic seasonal swings. The ACF plot displays all spikes lying within the 95% confidence blue dashed bands, indicating that the remaining variations behave essentially like white noise.
5) Let’s Generate and visualize the 2014 residendtial power usage Forecast with a 12 months observations.
# reporting forecastfc_2014 <-forecast(fit_arima_bc, h =12)fc_2014
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 2014 9465457 8407975 10678093 7903327 11391560
Feb 2014 8581952 7604389 9706173 7138886 10369106
Mar 2014 6580433 5856220 7409092 5510046 7895815
Apr 2014 5945084 5285857 6700206 4971004 7144112
May 2014 5886529 5234673 6633063 4923294 7071853
Jun 2014 8229114 7275475 9328683 6822251 9978397
Jul 2014 9523516 8397839 10825426 7864078 11596536
Aug 2014 10056917 8859403 11443516 8292071 12265528
Sep 2014 8492852 7504437 9633267 7034915 10307459
Oct 2014 5813572 5170885 6549424 4863830 6981853
Nov 2014 6097026 5418628 6874494 5094741 7331715
Dec 2014 7601354 6729884 8604540 6315203 9196549
# Plotingplot(fc_2014, main ="Residential KWH Power Usage Forecast for 2014",ylab ="Power Consumption (KWH)", xlab ="Time")grid()
Result
This plot exhibit the 1998–2013 historical power usage alongside the 2014 monthly forecast, which proves that the model successfully reproduces the characteristic bimodal seasonal pattern of peak consumption during summer and winter because they are cooling and heating months.
6) Let’s Add the forecast power usage to the existing files
# Let's generate the 2014 forecast datamonths_2014 <-c("2014-Jan", "2014-Feb", "2014-Mar", "2014-Apr", "2014-May", "2014-Jun", "2014-Jul", "2014-Aug", "2014-Sep", "2014-Oct", "2014-Nov", "2014-Dec")df_forecast <-data.frame(CaseSequence =seq(from =max(resid_power$CaseSequence) +1, length.out =12),`YYYY-MMM`= months_2014,KWH =round(as.numeric(fc_2014$mean), 0),check.names =FALSE)# Let's combine original data and new generated forecast datadf_updated <-bind_rows(resid_power , df_forecast)# Let's Write to fcombined Excel file output_file <-"ResidentialCustomerForecastLoad-624_Forecasted.xlsx"write_xlsx(df_updated, path = output_file)
