Forecasting with Machine Learning

library(tidyverse)
library(forecast)
library(gridExtra)
library(Metrics)
library(tvReg)
library(tsibble)
library(fable)
library(fable.prophet)
library(hms)
library(lubridate)

Load dataset

url_link = "https://gist.github.com/novrisuhermi/86d88b5d82cea37007979921c38fc207/raw/03ded2a94d48aa77dd18d37b11574801aaa62364/data2.csv"
data <- read_csv(url_link) %>% mutate(time=hms::as_hms(ts), day=day(ts))

Rows: 3600 Columns: 3
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl  (2): rate, occupancy
dttm (1): ts

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(data)

tail(data)

train    <- data %>% filter(day!=25)
test     <- data %>% filter(day==25)
data.ts  <- data %>% select(ts,occupancy,rate) %>% as_tsibble(index = ts)
train.ts <- train %>% select(ts,occupancy,rate) %>% as_tsibble(index = ts)
test.ts  <- test %>% select(ts,occupancy,rate) %>% as_tsibble(index = ts)

Modeling

Time series data

y <- ts(data = train$rate, frequency = 144)
x <- ts(data = train$occupancy, frequency = 144)
autoplot(y) + xlab("") + ylab("Energy consumption (kWh)")

Linear Regression

m <- list()
m[["LM"]] <- train.ts %>% model(lm=TSLM(rate~log(occupancy)))

m[["LM"]] %>% forecast(test.ts) %>% autoplot()

ARIMA

m[["ARIMA"]] <- auto.arima(y)
summary(m[["ARIMA"]])

Series: y 
ARIMA(1,0,2)(0,1,0)[144] 

Coefficients:
         ar1      ma1      ma2
      0.9662  -0.4559  -0.1216
s.e.  0.0053   0.0184   0.0183

sigma^2 = 151.4:  log likelihood = -13011.89
AIC=26031.78   AICc=26031.79   BIC=26056.2

Training set error measures:
                    ME     RMSE      MAE        MPE    MAPE      MASE         ACF1
Training set 0.1073758 12.04051 8.923401 -0.4944098 7.43582 0.5325004 -0.002891708

m[["ARIMA"]] %>% forecast::forecast(h=144) %>% autoplot()

ARIMAX

m[["ARIMAX"]] <- auto.arima(y, xreg = x)
summary(m[["ARIMAX"]])

Series: y 
Regression with ARIMA(1,0,3)(0,1,0)[144] errors 

Coefficients:
         ar1      ma1      ma2     ma3    xreg
      0.9346  -0.4514  -0.1321  0.0278  0.0364
s.e.  0.0097   0.0197   0.0198  0.0189  0.0037

sigma^2 = 148:  log likelihood = -12973.41
AIC=25958.82   AICc=25958.85   BIC=25995.45

Training set error measures:
                    ME     RMSE      MAE        MPE     MAPE      MASE          ACF1
Training set 0.0835479 11.90231 8.856013 -0.5202593 7.365541 0.5284791 -0.0002280971

m[["ARIMAX"]] %>% forecast::forecast(h=144, xreg=test$occupancy) %>% autoplot()

Prophet Forecast

m[["Prophet"]] <- train.ts %>% model(prophet = prophet(rate ~ growth("linear") + season("day", type = "multiplicative")))

components(m[["Prophet"]]) %>% autoplot()

m[["Prophet"]] %>% forecast(test.ts) %>% autoplot()

components(m[["Prophet"]])

Prophet Forecast with covariate

m[["Prophet2"]] <- train.ts %>% model(prophet = prophet(rate ~ growth("linear") + season("day", type = "multiplicative") + occupancy))

components(m[["Prophet2"]]) %>% autoplot()

m[["Prophet2"]] %>% forecast(test.ts) %>% autoplot()

components(m[["Prophet2"]])

NNETAR

m[["NNETAR"]] <- nnetar(y, size=20, lambda="auto")
m[["NNETAR"]]

Series: y 
Model:  NNAR(17,1,20)[144] 
Call:   nnetar(y = y, size = 20, lambda = "auto")

Average of 20 networks, each of which is
a 18-20-1 network with 401 weights
options were - linear output units 

sigma^2 estimated as 0.001969

m[["NNETAR"]] %>% forecast::forecast(h=144, PI=TRUE) %>% autoplot()

NNETAR with covariate

m[["NNETAR2"]] <- nnetar(y, size=20, xreg=x, lambda="auto")
m[["NNETAR2"]]

Series: y 
Model:  NNAR(17,1,20)[144] 
Call:   nnetar(y = y, size = 20, xreg = x, lambda = "auto")

Average of 20 networks, each of which is
a 19-20-1 network with 421 weights
options were - linear output units 

sigma^2 estimated as 0.001873

m[["NNETAR2"]] %>% forecast::forecast(xreg=test$occupancy, PI=TRUE) %>% autoplot()

Forecasting and Evaluation

Forecasting

prediction <- test %>% select(ts, rate) %>% rename(actual = rate, time = ts)
n_test <- dim(prediction)[1]
prediction['LM'] <- (m[["LM"]] %>% forecast(test.ts))$.mean
prediction['ARIMA'] <- (m[["ARIMA"]] %>% forecast::forecast(h=144))$mean
prediction['ARIMAX'] <- (m[["ARIMAX"]] %>% forecast::forecast(h=144, xreg=test$occupancy))$mean
prediction['Prophet'] <- (m[["Prophet"]] %>% forecast(test.ts))$.mean
prediction['Prophet2'] <- (m[["Prophet2"]] %>% forecast(test.ts))$.mean
prediction['NNETAR'] <- (m[["NNETAR"]] %>% forecast::forecast(h=144))$mean
prediction['NNETAR2'] <- (m[["NNETAR2"]] %>% forecast::forecast(xreg=test$occupancy))$mean

head(prediction)

plotForecast <- function(data){
  p1 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = LM)) + 
        xlab("") + ylab("") + ggtitle("Linear Model") + scale_x_datetime(date_labels = "%H:%M") 

  p2 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = ARIMA)) +
        xlab("") + ylab("") + ggtitle("ARIMA") + scale_x_datetime(date_labels = "%H:%M") 
  
  p3 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = ARIMAX)) +
        xlab("") + ylab("") + ggtitle("ARIMAX") + scale_x_datetime(date_labels = "%H:%M") 
  
  p4 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = Prophet)) + 
        xlab("") + ylab("") + ggtitle("Prophet") + scale_x_datetime(date_labels = "%H:%M") 
  
  p5 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = Prophet2)) + 
        xlab("") + ylab("") + ggtitle("Prophet with covariate") + scale_x_datetime(date_labels = "%H:%M")
  
  p6 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = NNETAR)) +
        xlab("") + ylab("") + ggtitle("NNAR") + scale_x_datetime(date_labels = "%H:%M") 
  
  p7 <- ggplot(data = data) + 
        geom_line(color = "red", aes(x = time, y = actual), size = 0.5, linetype = 2) +
        geom_line(color = "black", aes(x = time, y = NNETAR2)) +
        xlab("") + ylab("") + ggtitle("NNAR with covariate") + scale_x_datetime(date_labels = "%H:%M") 
   
        
  grid.arrange(p1,p2,p3,p4,p5,p6,p7,nrow=3)
}

plotForecast(prediction)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Evaluation

RMSE

sort(apply(X = prediction[,-c(1,2)], MARGIN = 2, FUN = rmse, actual = prediction$actual))

Prophet2   ARIMAX       LM    ARIMA  NNETAR2  Prophet   NNETAR 
10.58774 12.36074 12.78291 13.31192 13.82411 16.01882 21.33545 

MAPE

sort(apply(X = prediction[,-c(1,2)], MARGIN = 2, FUN = mape, actual = prediction$actual))

  Prophet2     ARIMAX      ARIMA         LM    Prophet    NNETAR2     NNETAR 
0.07253363 0.07927125 0.08764947 0.09318489 0.10055044 0.10085085 0.14022036 

sMAPE

sort(apply(X = prediction[,-c(1,2)], MARGIN = 2, FUN = smape, actual = prediction$actual))

  Prophet2     ARIMAX      ARIMA         LM    NNETAR2    Prophet     NNETAR 
0.07089585 0.07606659 0.08267374 0.08474103 0.09000829 0.10201565 0.13480095