1 Introduction

Pada kali ini kita akan melakukan analisi time series dan forecasting dari data yang telah diunduh dari https://www.kaggle.com/shenba/time-series-datasets?select=Electric_Production.csv. Data tersebut berisi data produksi listrik tiap bulan (pencatatan di setiap tanggal 1 awal bulan) dari Januari 1985 sampai dengan Januari 2018.

2 Library & Setup

Adapun library yang kita gunakan untuk analisis kali ini adalah sebagai berikut :

# import library
library(tidyverse)
library(lubridate)
library(forecast)
library(TTR)
library(fpp)
library(tseries)
library(TSstudio)
library(padr)

3 Data Preparation

3.1 Read Data

Data yang akan digunakan untuk analisis kali ini adalah data Electric_Production.csv.

# read data
electric <- read.csv("data_input/Electric_Production.csv")
head(electric)
tail(electric)
# periksa data 
head(electric)

Cek dimensi data

# cek dimensi
dim(electric)
## [1] 397   2

Dari data diatas terdapat 397 baris dan 2 kolom

Deskripsi Data :

  1. DATE : Tanggal, Bulan dan Tahun produksi
  2. IPG2211A2N : Banyaknya Industrial Production perbulan (pencatatan dilakukan di awal bulan)

3.2 Karakteristik Time Series

1. Data harus terurut berdasarkan periode waktunya

# mengubah tipe data dan mengurutkan data dikolom date
electric <- electric %>% 
  mutate(DATE = mdy(DATE),
         MONTH = month(DATE, label = T)) %>% 
  arrange(DATE)
electric

Dari data yang di tampilkan diatas, data electric telah diurutkan berdasarkan tanggal awal di tiap bulan yang ada dikolom DATE.

2. Tidak boleh ada data periode/waktu yang terlewat

# periksa data
electric

Dari data diatas tidak terdapat data yang terlewat, maka bisa kita bisa lanjutkan kelangkah selanjutnya.

3. Data tidak boleh ada yang missing

# cek missing value
electric %>%
  is.na() %>% # cek NA
  colSums()
##       DATE IPG2211A2N      MONTH 
##          0          0          0

Dari data diatas tidak terdapat data yang missing value, maka bisa kita bisa lanjutkan kelangkah selanjutnya.

4 Time Series

Langkah selanjutnya adalah membuat objek time series, dengan frequency = ‘12’

# create ts object
electric_ts <- ts(data = electric$IPG2211A2N, start = ymd("1985-01-01"), frequency = 12)

Buat visualisasi

# visualize the data
electric_ts %>% autoplot()

ggplot(electric, aes(x = DATE, y = IPG2211A2N)) +
  geom_point(data = electric %>% filter(MONTH %in% c("Jan", "Apr", "Jul")),
             mapping = aes(x = DATE, y = IPG2211A2N, color = MONTH)) +
  scale_color_manual(values = c("red", "blue", "green")) +
  geom_line()

Berdasarkan plot di atas produksi barang electric IPG2211A2N terus mengalami peningkatan dari tahun ke tahun, terlihat juga terjadi peningkatan (setiap bulan Juli) dan penurunan (setiap bulan April) yang berulang di setiap tahunnya. Sehingga dapat disimpulkan bahwa terdapat efek seasonal tahunan pada data kita.

5 Decomposition

Kita gunakan decomposition dari fungsi decompose untuk menguraikan komponen-komponen dalam time series data.

# decompose
electric_dc <- decompose(x = electric_ts, type = "additive")

# visualize decomposition result
electric_dc %>%  autoplot()

# Analisis Seasonal

electric %>% 
  mutate(seasonal = electric_dc$seasonal) %>% 
  tail(12) 
electric %>% 
  mutate(seasonal = electric_dc$seasonal) %>% 
  tail(12) %>% 
  ggplot(aes(x = MONTH, y = seasonal)) +
  geom_col()

Dari plot di atas diketahui bahwa seasonal yang terjadi adalah setiap 12 bulan sekali, dimana pada bulan April dan Mei seasonal sangat negatif yang menandakan terjadinya penurunan. Sementara, pada bulan Juli dan Agustus seasonal sangat positif yang menandakan terjadinya peningkatan.

6 Build Model

6.1 Cross Validation

# untuk test (data yang paling baru)
test_electric <- tail(electric_ts, 12)
# untuk train (data yang lebih lama)
train_electric <- head(electric_ts, length(electric_ts) - length(test_electric))

6.2 Model Time Series

Setelah Anda memisahkan data electric_ts menjadi data train dan test, silahkan analisis pola tren dan seasonal pada data train_electric.

# your code here
decompose(x = train_electric, type = "additive") %>%  autoplot()

Berdasarkan grafik dekomposisi diatas, kita menggunakan model Holt-Winters, karena grafik tersebut mengandung pola trend dan seasonality.

Setelah kita menganalisis hasil dekomposisi pada data train_electric, kita sudah siap untuk membangun model. kita akan membuat model pertama menggunakan algoritma dari Holt-Winters. Kita dapat menggunakan fungsi HoltWinters().

# your code here
model_hw <- HoltWinters(x = train_electric)

Kita juga akan membuat model untuk meramalkan data train_electric menggunakan algoritme ARIMA. Kita akan buat sebuah model ARIMA menggunakan fungsi stlm().

# your code here
model_arima <- stlm(y = train_electric, s.window = 12, method = "arima")

7 Forecasting

Kita telah membuat model menggunakan Holt-Winters dan ARIMA, sekarang kita akan meramalkan banyaknya produksi barang IPG2211A2N dalam 12 bulan kedepan menggunakan fungsi forecast().

# your code here
hw_forecast <- forecast(model_hw, 12)
# your code here
arima_forecast <- forecast(model_arima, h = 12)

8 Evaluasi Model (Error)

Holt-Winters

# model evaluation
accuracy(hw_forecast, test_electric)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.07063815 2.483562 1.958386 0.01704223 2.150216 0.6986856
## Test set     1.28248169 4.948274 3.445381 1.02927682 3.175984 1.2291947
##                   ACF1 Theil's U
## Training set 0.1454255        NA
## Test set     0.1871180 0.4228492
accuracy(hw_forecast$mean, test_electric)
##                ME     RMSE      MAE      MPE     MAPE     ACF1 Theil's U
## Test set 1.282482 4.948274 3.445381 1.029277 3.175984 0.187118 0.4228492

ARIMA

# model evaluation
accuracy(arima_forecast, test_electric)
##                        ME     RMSE      MAE          MPE     MAPE      MASE
## Training set -0.003738892 1.942161 1.484658 -0.005357833 1.638681 0.5296755
## Test set     -0.214046896 4.134059 2.872826 -0.423369133 2.639645 1.0249266
##                      ACF1 Theil's U
## Training set -0.003086411        NA
## Test set      0.241600170  0.336749
accuracy(arima_forecast$mean, test_electric)
##                  ME     RMSE      MAE        MPE     MAPE      ACF1 Theil's U
## Test set -0.2140469 4.134059 2.872826 -0.4233691 2.639645 0.2416002  0.336749

Dari perbandingan nilai MAPE kedua model di atas, dapat disimpulkan bahwa model ARIMA lebih baik dibandingkan Holt-Winters, karena nilai MAPE lebih kecil.

9 Cek Asumsi

No-Autocorrelation

Menggunakan tes Ljung-Box

# your code here
Box.test(arima_forecast$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  arima_forecast$residuals
## X-squared = 0.0036961, df = 1, p-value = 0.9515

Hipotesis:

  • H0: error has no-autocorrelation
  • H1: error has autocorrelation

p-value = 0.9515, artinya > 0.05 -> Terima H0