Electricity
Before electricity generation began over 100 years ago, houses were lit with kerosene lamps, food was cooled in iceboxes, and rooms were warmed by wood-burning or coal-burning stoves. Beginning with Benjamin Franklin’s experiment with a kite one stormy night in Philadelphia, the principles of electricity gradually became understood. Electricity is a basic part of nature and it is one of our most widely used forms of energy. Many cities and towns were built alongside waterfalls (a primary source ofmechanical energy) that turned water wheels to perform work.
Dataset
This dataset contain hourly power consumption reported by AEP. AEP a major investor-owned electric utility in the United States, delivering electricity to more than five million customers in 11 states. With this dataset, we are going to predict power consumption for the next four weeks with Time Series.
power <- read.csv("AEP_hourly.csv")
str(power)## 'data.frame': 121273 obs. of 2 variables:
## $ Datetime: chr "2004-12-31 01:00:00" "2004-12-31 02:00:00" "2004-12-31 03:00:00" "2004-12-31 04:00:00" ...
## $ AEP_MW : num 13478 12865 12577 12517 12670 ...About the columns :
Datetime: Date and time of the dataAEP_MW: Reported Power consumption by AEP, in Megawatt(MW)
Data Wrangling
Working with Time Series, there are conditions to be met :
Period has to be sorted from oldest to newest
Same interval between each data
No missing interval
No missing values
power <- power %>%
mutate(Datetime = ymd_hms(Datetime)) %>%
arrange(Datetime) #arrange the dataset by its DatetimeFinding the range time of the dataset
power$Datetime %>%
range()## [1] "2004-10-01 01:00:00 UTC" "2018-08-03 00:00:00 UTC"Create a sequence list of date based on the range, and compare it to the dataset to find missing date. This dataset is missing 27 rows.
date_range <- seq(from = as.POSIXct("2004-10-01 01:00:00", tz="UTC"), to = as.POSIXct("2018-08-03 00:00:00", tz="UTC"), by = "hour")
date_range[!date_range %in% power$Datetime]## [1] "2004-10-31 02:00:00 UTC" "2005-04-03 03:00:00 UTC"
## [3] "2005-10-30 02:00:00 UTC" "2006-04-02 03:00:00 UTC"
## [5] "2006-10-29 02:00:00 UTC" "2007-03-11 03:00:00 UTC"
## [7] "2007-11-04 02:00:00 UTC" "2008-03-09 03:00:00 UTC"
## [9] "2008-11-02 02:00:00 UTC" "2009-03-08 03:00:00 UTC"
## [11] "2009-11-01 02:00:00 UTC" "2010-03-14 03:00:00 UTC"
## [13] "2010-11-07 02:00:00 UTC" "2010-12-10 00:00:00 UTC"
## [15] "2011-03-13 03:00:00 UTC" "2011-11-06 02:00:00 UTC"
## [17] "2012-03-11 03:00:00 UTC" "2012-11-04 02:00:00 UTC"
## [19] "2012-12-06 04:00:00 UTC" "2013-03-10 03:00:00 UTC"
## [21] "2013-11-03 02:00:00 UTC" "2014-03-09 03:00:00 UTC"
## [23] "2014-03-11 14:00:00 UTC" "2015-03-08 03:00:00 UTC"
## [25] "2016-03-13 03:00:00 UTC" "2017-03-12 03:00:00 UTC"
## [27] "2018-03-11 03:00:00 UTC"Fill the missing row and fill the power consumption with the average between previous and next rows.
power <- power %>%
pad(start_val = min(power$Datetime), end_val = max(power$Datetime)) %>%
mutate(AEP_MW = na.approx(AEP_MW))## pad applied on the interval: hourRecheck missing date.
date_range[!date_range %in% power$Datetime]## POSIXct of length 0Exploratory Data Analysis
The least power consumption is 9581 MW, on Sunday, 2 October 2016 05:00. And the most power consumption is 25.695 MW, on Monday, 20 October 14:00. Both date are in the real data(not in the missing date). Both mean and median of power consumption are at 15.000 MW.
summary(power)## Datetime AEP_MW
## Min. :2004-10-01 01:00:00 Min. : 9581
## 1st Qu.:2008-03-17 13:45:00 1st Qu.:13629
## Median :2011-09-02 02:30:00 Median :15309
## Mean :2011-09-02 01:51:02 Mean :15499
## 3rd Qu.:2015-02-16 14:15:00 3rd Qu.:17200
## Max. :2018-08-03 00:00:00 Max. :25695power[power$AEP_MW %in% c(9581,25695),]## Datetime AEP_MW
## 35534 2008-10-20 14:00:00 25695
## 105223 2016-10-02 05:00:00 9581Time Series
Explore the dataset by creating a timeseries object, and find the trend and seasonal
power$AEP_MW %>%
head(24*7*4) %>%
ts(start = 2004, frequency = 24) %>%
decompose() %>%
autoplot()
From the plot, we can see trend is showing ups and downs, indicating this is a multiseasonal timeseries dataset. Hence, we have to reanalyze this dataset with multi seasonal timeseries.
power$AEP_MW %>%
head(24*7*4) %>%
msts(seasonal.periods = c(24, 24*7)) %>%
mstl() %>%
autoplot()
Based on the multiseasonal timeseries plot above, we can see the trend, and there are daily and weekly seasonal. And from the data plot, we can assume, this is a multiplicative timeseries dataset.
#assign the multiseasonal timeseries
power_msts <- power$AEP_MW %>%
msts(seasonal.periods = c(24, 24*7))Cross Validation
Split the dataset into train and test dataset. Since this is a timeseries dataset, splitting should be done sequentially, not randomly. For test dataset, we will forecast the power consumption for the next 4 weeks. And the rest of the data, goes to train dataset. And since this is a big and grain dataset, we will only use 30 months of data for the train dataset.
power_test <- tail(power_msts, 24*7*4)
power_train <- head(power_msts, -length(power_test))
power_train <- tail(power_train, 24*7*4*30)Model & Forecast
We will create three model with Holt-Winters, Seasonal Arima(SARIMA) and TBATS.
Holt-Winters
pmodel_hw <- HoltWinters(x = power_train, seasonal = "multiplicative")
ppred_hw <- forecast(object = pmodel_hw, h = 24*7*4)SARIMA
For Seasonal ARIMA, data has to be stationary. From Augmented Dickey-Fuller test shows p-value 0.01 < 0.05 (alpha), indicating our timeseries already stationary, and no differencing is needed.
options(warn = - 1)
adf.test(power_train)##
## Augmented Dickey-Fuller Test
##
## data: power_train
## Dickey-Fuller = -11.943, Lag order = 27, p-value = 0.01
## alternative hypothesis: stationarypmodel_stlm_sarima <- stlm(power_train, method = "arima", lambda = 0)
ppred_stlm_sarima <- forecast(pmodel_stlm_sarima, h = 24*7*4)TBATS
# pmodel_tbats <- tbats(power_train, use.box.cox = FALSE,
# use.trend = TRUE,
# use.damped.trend = TRUE)
#
# saveRDS(pmodel_tbats, file = "pmodel_tbats.RDS")
pmodel_tbats <- readRDS("pmodel_tbats.RDS")
ppred_tbats <- forecast(pmodel_tbats, h = 24*7*4)Model Evaluation
TBATS is the best model with lowest MAPE(Mean Absolute Percentage Error) at 8.6%. And SARIMA at 9.4%. Both model MAPE below 10%. We can also see the forecast from the plot below.
accuracy_result <- rbind(accuracy(ppred_hw$mean, power_test),
accuracy(ppred_stlm_sarima$mean, power_test),
accuracy(ppred_tbats$mean, power_test))
rownames(accuracy_result) <- c("Holt-Winters","SARIMA","TBATS")
accuracy_result## ME RMSE MAE MPE MAPE ACF1
## Holt-Winters -1122.5396 3041.602 2403.169 -9.793437 16.863292 0.9559026
## SARIMA -1345.2072 1800.833 1441.590 -8.839780 9.435308 0.9737576
## TBATS 566.4253 1790.770 1400.362 2.532889 8.601969 0.9776165
## Theil's U
## Holt-Winters 5.451440
## SARIMA 2.778478
## TBATS 2.462018Holt-Winters
power_train %>%
autoplot() +
autolayer(ppred_hw, series = "hw") +
autolayer(power_test, series = "actual")
SARIMA
power_train %>%
autoplot() +
autolayer(ppred_stlm_sarima, series = "sarima") +
autolayer(power_test, series = "actual")
TBATS
power_train %>%
autoplot() +
autolayer(ppred_tbats, series = "tbats") +
autolayer(power_test, series = "actual")
Assumption
There are two assumption in time series model :
- Normality of Residual
- No-Autocorrelation Residual
We will test two of our model, SARIMA and TBATS, on these assumption.
Normality of Residual
\(H_0\): residual distributed normally
\(H_1\): residual not distributed normally
We can use Shapiro-Wilk normality test and histogram. But since Shapiro-Wilk test can only handle 5000 data, we will add Anderson-Darling normality test. Even tough the histogram for SARIMA and TBATS shows bell-curve, but both Shapiro-Wilk and Anderson-Darling tests for SARIMA and TBATS shows p-value < 0.05, meaning : Reject H0, accept H1,Residual are not distributed normally. Our model did not met the Normality of Residual Assumption.
Shapiro-Wilk Test
SARIMA
set.seed(310)
x <- as.data.frame(ppred_stlm_sarima$residuals) %>%
sample_n(5000)
shapiro.test(ts(x))##
## Shapiro-Wilk normality test
##
## data: ts(x)
## W = 0.9691, p-value < 2.2e-16TBATS
set.seed(310)
x <- as.data.frame(ppred_tbats$residuals) %>%
sample_n(5000)
shapiro.test(ts(x))##
## Shapiro-Wilk normality test
##
## data: ts(x)
## W = 0.98937, p-value < 2.2e-16Anderson-Darling Test
SARIMA
nortest::ad.test(ppred_stlm_sarima$residuals)##
## Anderson-Darling normality test
##
## data: ppred_stlm_sarima$residuals
## A = 54.472, p-value < 2.2e-16TBATS
nortest::ad.test(ppred_tbats$residuals)##
## Anderson-Darling normality test
##
## data: ppred_tbats$residuals
## A = 29.241, p-value < 2.2e-16Histogram
SARIMA
hist(ppred_stlm_sarima$residuals)
plot(density(ppred_stlm_sarima$residuals))
TBATS
hist(ppred_tbats$residuals)
plot(density(ppred_tbats$residuals))
No-Autocorrelation Residual
\(H_0\): Residual has No-Autocorrelation
\(H_1\): Residual has Autocorrelation
We can use acf plot and Ljung-Box test to test No-Autocorrelation. A good acf plot shows Lag not crossing the limit line. Our model shows some lags are cross the limit line, this could indicate Autocorrelation. But the better way to test Autocorrelation is by using the Ljung-Box test.
SARIMA model, with p-value 0.9565 > 0.05, meaning : Fail to reject H0, residual has No-Autocorrelation. TBATS mode, with p-value < 0.05, meaning : reject H0, residual has Autocorrelation.
acf plot
SARIMA
acf(pmodel_stlm_sarima$residuals)
TBATS
acf(ppred_tbats$residuals)
Ljung-Box
SARIMA
Box.test(pmodel_stlm_sarima$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: pmodel_stlm_sarima$residuals
## X-squared = 0.028168, df = 1, p-value = 0.8667TBATS
Box.test(ppred_tbats$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: ppred_tbats$residuals
## X-squared = 135.1, df = 1, p-value < 2.2e-16Conclusion
Both model did not met the Normality Assumption. While doing a better job on forecasting, TBATS also did not met the No-Autocorrelation Assumption. On the other hand, SARIMA is slightly below TBATS, but it met the No-Autocorrelation Assumption. Personally, I will choose the SARIMA model.