STL with Multiple Seasonal Periods

In this case, we’re going to predict the value of hourly power consumption data comes from PJM’s website and are in megawatts(Mw). PJM Interconnection LLC (PJM) is a regional transmission organization (RTO) in the United States. It is part of the Eastern Interconnection grid operating an electric transmission system serving all or parts of Delaware, Illinois, Indiana, Kentucky, Maryland, Michigan, New Jersey, North Carolina, Ohio, Pennsylvania, Tennessee, Virginia, West Virginia, and the District of Columbia. The data contains the hourly power consumption estimated by each electricity provider. dataset.

#load the package
library(lubridate)
library(dplyr)
library(forecast)
library(ggplot2)
library(scales)

#import data
pjm<- read.csv("data_input/pjm.csv") %>% filter(provider == "AEP")
str(pjm)

#> 'data.frame':    45336 obs. of  3 variables:
#>  $ datetime: Factor w/ 45336 levels "2013-06-01T00:00:00Z",..: 1 2 3 4 5 6 7 8 9 10 ...
#>  $ provider: Factor w/ 8 levels "AEP","COMED",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ cons    : num  13477 12699 12274 11904 11862 ...

After you have read the time series data into R, the next step is to store the data in a time series object in R, so that you can use R many functions for analysing time series data. To store the data in a time series object, we can use the ts() function in R. Sometimes the time series data set that you have may have been collected at regular intervals that were less than one year, for example, monthly or quarterly. In this case, you can specify the number of times that data was collected per hour by using the frequency parameter in the ts() function. Because our each row representing a data within hourly interval, you can set frequency=24, and we will only used AEP provider.

pjm$datetime <- ymd_hms(pjm$datetime)
ts_train<-pjm$cons %>% ts(freq= 24)

To estimate the trend component and seasonal component of a seasonal time series, we can use the decompose() function in R. This function estimates the trend, seasonal, and irregular components of a time series. Now let’s inspect the seasonal and trend component.

ts_train %>% 
  tail(24*7*4) %>% 
  decompose() %>% 
  autoplot()

Decomposing a time series means separating it into it’s constituent components.

1. The first panel from top is the original, observed time series. Note that a series with multiplicative effects can often be transformed into one with additive effect through a log transformation.

2. The second panel plots the seasonal component. with the figure being computed by taking the average for each time unit over all periods and then centering it around the mean.

3. The third panel plots the trend component,We see that the estimated trend component shows a pattern. This pattern in trend might be sourced from uncaptured extra seasonality from higher natural period in this case,so it can be considered as multi-seasonal data. To solve this complex seasonality, we need to convert the data into msts() object which accept multiple frequency setting.

4. The bottom-most panel the error component, which is determined by removing the trend and seasonal figure.
   To deal with such series, we will use the msts class which handles multiple seasonality time series. This allows you to specify all of the frequencies that might be relevant. It is also flexible enough to handle non-integer frequencies.
   

msts_cons<-pjm$cons %>% msts( seasonal.periods = c(24, 24*7))
msts_cons  %>% head(  24 *7 *4 ) %>% mstl() %>% autoplot()    

Now we can see a clearer trend and could confirm the daily and weekly seasonality for the data. Before we’re going to modelling our data, we’re going to split our data become data train and data test.

msts_train <- head(msts_cons, length(msts_cons) - 24*7)
msts_test <- tail(msts_cons,  24*7)

#subset to more recent period
msts_train <- tail(msts_train, 24*7*4*3)
autoplot(msts_train)

We will make prediction for 7 days and plot them. In modelling the log transformation, we can use lambda = 0 in stlm() setting.

stlm_model <- msts_train %>%
  stlm(lambda = 0) %>% 
  forecast(h = 24*7) 
plot(stlm_model)