Forecast of Electricity Demand in Malta

library(tidyverse)
library(readxl)
library(ggplot2)
library(lubridate)
library(writexl)
library(magrittr)

1, Tidy Data

Here is the history data about monthly electricity demand satisfied by non-renewable energy technologies from 2013 to 2017, which are downloaded from Malta National Statistics Office.

head(ti)

ti %>%
  mutate(month = month(.$date)) %>%
  mutate(year = as.character(year(.$date))) %>%
  ggplot() +
    geom_line(mapping = aes(x = month, y = demand, color = year))

ti %>%
  ggplot() +
    geom_line(mapping = aes(x = date, y = demand)) +
    labs(
      title = "Monthly Electricity Demand from Non-Renewables",
      subtitle = "in Malta from 2013 to 2017",
      y = "Electricity Demand from Non-Renewables (MWh)",
      x = "Index (Month)"
      )

2, Decomposition

ts <-
  ts(ti$demand, frequency = 12)

li_decomp <-
  ts %>%
  decompose(type = c("multiplicative"), filter = NULL)

The seasonal compoents is the most important part, and it’s shown in the following figure.

ts <-
  ts(ti$demand, frequency = 12)
li_decomp <-
  ts %>%
  decompose(type = c("multiplicative"), filter = NULL)
ti %>%
  mutate(coef_decomp = as.numeric(li_decomp$seasonal)) %>%
  ggplot() +
    geom_line(mapping = aes(x = date, y = coef_decomp)) +
    labs(
      title = "Seasonal Coefficients of Decomposition",
      subtitle = "for Monthly Non-Re E-Demand in Malta from 2013 to 2017",
      y = "Seasonal Coefficients",
      x = "Index (Month)"
      )

The random component is checked using ACF and PACF.

acf(li_decomp$random, na.action = na.pass)

pacf(li_decomp$random, na.action = na.pass)

We can say that the random component can be treated as white noise.

3, Model Trending Component using Linear Regression

acf(li_decomp$trend, na.action = na.pass)

pacf(li_decomp$trend, na.action = na.pass)

arima(as.ts(li_decomp$trend))

Call:
arima(x = as.ts(li_decomp$trend))

Coefficients:
        intercept
      178153.1495
s.e.     337.1318

sigma^2 estimated as 5455566:  log likelihood = -440.4,  aic = 884.8

Though there is some correlation in ACF of the trending component, the “arima” function shows that no benefit, which is the same as expected becuase the length of the sampling is one month. So linear regression is used.

df_trend <- data.frame(x = 1:length(li_decomp$trend), y = li_decomp$trend)
lm_trend <- lm(y ~ x, df_trend)
par(mfrow = c(2, 2))
plot(lm_trend)

The illustration of linear regression anad forecast of trending component is shown in the following figure.

par(mfrow = c(1, 1))
plot(
  predict(
    lm_trend, newdata = data.frame(x = c(1:(length(df_trend$x) + 12 * 10)))
    )
  )
lines(df_trend$x, df_trend$y)

4, Forecast

trend_new <- predict(
  lm_trend, newdata = data.frame(
    x = c((length(df_trend$x) + 1):(length(df_trend$x) + 12 * 20))
    )
  )
seasonal_new <- rep(li_decomp$seasonal[1:12], 20)
ti_new <-
  tibble(
    index = 1:(12 * 20), demand = trend_new * seasonal_new  # multiplicative
    ) %>%
  mutate(date = rep(ti$date[60], length(.$index)))
for (i in 1:length(ti_new$index)) {
  month(ti_new$date[i]) = i + month(ti$date[60])
}
ti_new %<>%
  select(date, demand) %>%
  mutate(observation = rep(FALSE, length(.$date)))
ti_fore <- ti %>%
  bind_rows(ti_new)
ti_fore %>%
  ggplot() +
    geom_line(mapping = aes(x = date, y = demand, color = observation)) +
    labs(title = "Forecast of Monthly Electricity Demand from Non-Renewables",
         subtitle = "in Malta from 2018 to 2037 using Data from 2013 to 2017",
         y = "Electricity Demand from Non-Renewables (MWh)",
         x = "Index (Month)")