1.Intro

Data berikut merupakan data Produksi Industri: listrik dan gas dari tahun 2000-2019. Dari data ini akan dilakukan prediksi untuk melihat berapa energi yang dihasilkan pada pabrik yang terletak di US.

1.1 Library

Setup library

library(tidyverse)
library(lubridate)
library(padr)
library(zoo)
library(fpp)
library(TSstudio)
library(forecast)
library(TTR)
library(tseries)

theme_set(theme_minimal())

1.2 Read Data

Membaca data dari dataset IPG2211A2N

indu <- read.csv("data/IPG2211A2N.csv") %>% 
  tail(228) ## Untuk memotong jumlah data yang tadinya dari 1939 menjadi mulai dari 2000
str(indu)

## 'data.frame':    228 obs. of  2 variables:
##  $ DATE      : chr  "2000-06-01" "2000-07-01" "2000-08-01" "2000-09-01" ...
##  $ IPG2211A2N: num  91.3 96.3 99.7 91 83.4 ...

1.3 Data Preprocessing

Merubah tipe data DATE menjadi tipe data datetime menggunakan libarary lubridate

indu$DATE <- ymd(indu$DATE)

Cek Range atau periode waktu

range(indu$DATE)

## [1] "2000-06-01" "2019-05-01"

Cek apakah data sudah memenuhi syarat data time series

indu %>% 
  arrange(DATE) %>% 
  pad() %>% 
  anyNA()

## pad applied on the interval: month

## [1] FALSE

1.4 Eksplorasi Data

ggplot(data = indu, aes(x = DATE, y = IPG2211A2N)) +
  geom_point() +
  geom_line()

2.Time Series

2.1 Membuat object time series

Melakukan perubahan data menjadi timeseries (ts) menggunakan library ts.

indu_ts <- ts(data = indu$IPG2211A2N, start = c(2000,1), frequency = 12)

2.3 Visualisasi object time series

indu_ts %>% 
  autoplot()

2.4 Decomposition

Menggunakan fungsi decompose() menjadi beberapa plot yaitu: 1. Trend 2. Seasonal 3. Random/Error

indu_dc <- decompose(indu_ts)

2.5 Visualisasi hasil decompose

plot(indu_dc)

2.6 Seasonal analysis

autoplot(indu_dc$x - indu_dc$trend)

2.7 Trend Analysis

autoplot(indu_dc$x - indu_dc$seasonal)

3. Cross Validation

Sebelum membuat model, pisahkan data menjadi data train dan test sebanyak 36 bulan atau 3 tahun

indu_test <- indu_ts %>% 
  tail(36)
indu_train <- indu_ts %>% 
  head(-36)

4. Model Building

4.1 Holtwinters

model_hol_indu <- HoltWinters(indu_train, seasonal = "multiplicative")
forecast_model_hol_indu <- forecast(model_hol_indu, h = 36)

Cek akurasi dari hasil prediksi Holtwinters

accuracy(forecast_model_hol_indu, indu_test)

##                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 0.2323497 2.897883 2.204945 0.1856393 2.190787 0.7317331 0.1263927
## Test set     5.3113959 6.794215 5.657476 4.9389534 5.286636 1.8774901 0.4104344
##              Theil's U
## Training set        NA
## Test set     0.6421105

4.2 STLM (Sarima)

Metode SARIMA digunakan karena pada saat dilakukan decompose() terdapat seasonal. Untuk membuat model SARIMA dengan fungsi stlm().

model_st_indu <- stlm(y = indu_train, method = "arima")
summary(model_st_indu$model)

## Series: x 
## ARIMA(2,1,1) 
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5093  -0.054  -0.8833
## s.e.  0.0872   0.080   0.0519
## 
## sigma^2 estimated as 5.494:  log likelihood=-432.57
## AIC=873.14   AICc=873.36   BIC=886.15
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.2053264 2.319345 1.828789 0.1640977 1.84274 0.6098118
##                      ACF1
## Training set -0.005504033

Membuat prediksi untuk model SARIMA

forecast_model_st_indu <- forecast(model_st_indu, h = 36)

forecast_model_st_indu %>% 
  autoplot() +
  autolayer(indu_test)

5. Evaluation

5.1 Komparasi model

Membandingkan performa masing-masing model dengan data test.

compare_forecast <- function(x, test){
  lapply(x, function(x) forecast::accuracy(x, test)["Test set", ]) %>%
    lapply(., function(x) x %>% t() %>% as.data.frame) %>% 
    bind_rows() %>% 
    mutate(model = names(x)) %>%
    select(model, everything())}

forecast_list <- list(
  "Holt_winters" = forecast_model_hol_indu,
  "SARIMA" = forecast_model_st_indu)

compare_forecast(forecast_list, indu_test)

##          model       ME     RMSE      MAE      MPE     MAPE    MASE      ACF1
## 1 Holt_winters 5.311396 6.794215 5.657476 4.938953 5.286636 1.87749 0.4104344
## 2       SARIMA 3.279125 5.043834 4.003223 3.001917 3.733724 1.32851 0.3463216
##   Theil's U
## 1 0.6421105
## 2 0.4721886

INSIGHT: Dari hasil evaluasi didapatkan MAPE pada model SARIMA lebih kecil 2% yang membuat model SARIMA dapat lebih tepat digunakan untuk melakukan prediksi keluaran energi dari proses produksi pabrik.

5.2 Visualisasi

test_forecast() digunakan untuk melihat keakuratan model SARIMA untuk prediksi data test

test_forecast(actual = indu_ts,
              forecast.obj = forecast_model_st_indu,
              train = indu_train,
              test = indu_test)

6. Assumption

6.1 No auto-correlation residual

Harapan dilakukan test ini adalah p-values > alpha (0.05)

Box.test(model_st_indu$residuals, type = 'Ljung-Box')

## 
##  Box-Ljung test
## 
## data:  model_st_indu$residuals
## X-squared = 0.0059079, df = 1, p-value = 0.9387

6.2 Normality residual