1 Library

library(dplyr) # for data wrangling
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate) # to dea with date
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(padr) # for padding
library(forecast) # time series library
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries) # for adf.test
library(MLmetrics) # calculate error
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall

2 Objektif

Pada kali ini saya akan melakukan time series dengan dataset dari electric production. Dataset yang saya gunakan ini berasal dari https://www.kaggle.com/kandij/electric-production

3 Read Data

data_elektronik <- read.csv("dataset/Electric_Production2.csv")
data_elektronik

4 Data Preprocessing

Selanjutnya dengan melakukan preprocessing data dari data yang saya gunakan, maka saya akan sesuaikan dengan aturan format tanggalnya. Untuk mengubah format tanggal, saya menggunakan library(lubridate) dengan format tanggal mdy atau month-date-year.

data_elektronik_clean <- data_elektronik %>% 
  mutate(DATE = mdy(DATE))

data_elektronik_clean

5 Missing Value

Langkah selanjutnya yaitu dengan mengecek apakah interval waktu ada yang terlewat atau tidak dan apakah data ada missing atau tidak.

data_elektronik_clean %>% 
  pad() %>% 
  is.na() %>% 
  colSums()
## pad applied on the interval: month
##  DATE Value 
##     0     0

Dari data yang saya gunakan, tidak ada interval waktu yang terlewat dan tidak ada data yang hilang.

6 Membuat Object Time Series

Selanjutnya, membuat object time series. Pada tahap ini membuat dengan menggunakan fungsi ts(). Dari data yang saya punya, data tersebut memiliki nilai frequency yaitu 12, yang berarti data saya yaitu data bulanan per tahun.

elektronik_ts <- ts(data = data_elektronik_clean$Value, start = c(1985,1),frequency = 12)

6.1 Visualisasi Hasil Objek Time Series

elektronik_ts %>% 
  autoplot()

Dari hasil visualisasi yang saya dapatkan, maka saya simpulkan dari data yang miliki dengan type = addictive. Karena jarak antar garis tidak begitu signifikan.

7 Splitting

Tahap selanjutnya yaitu dengan melakukan proses splitting dengan membagi kedalam data train dan data test. Untuk data test saya menggunakan 3 tahun dan data train yaitu sisa dari data test.

# data test
elektronik_test <- tail(elektronik_ts, 3*12)
# data train
elektronik_train <- head(elektronik_ts, length(elektronik_ts)-length(elektronik_test))

8 Decompose

Tahap selanjutnya yaitu melakukan decompose() terhadap data train diatas. Decompose ini digunakan untuk menentukkan addictive atau multiplicative.

elektronik_train %>% 
  decompose(type = "additive") %>% 
  autoplot()

9 Build Model

Tahap selanjutnya, model yang sudah didapatkan diatas akan dilakukan build model dengan menggunakan beberapa metode. Beberapa metode yang dibawah ini:

9.1 Holt Winter

Metode pertama dengan menggunakan HoltWinter. Pada metode ini saya melakukan uji coba dengan dua kondisi yaitu: 1. Kondisi pertama model dengan nilai smoothing nya auto. 2. Model kedua dengan nilai smoothing dengan mendefinisikan nilai alpha, beta dan gamma secara manual.

# 1
elektronik_holt <- HoltWinters(x = elektronik_train, seasonal = "additive")

# 2
elektronik_holt_manual <- HoltWinters(elektronik_train, alpha = 0.8, beta = 0.01, gamma = 1)

9.1.1 Forecasting

Tahap selanjutnya melakukan forecasting dari model holt winter tersebut. Forecasting ini untuk 3 tahun kedepan.

elektronik_forecast_holt <- forecast(object = elektronik_holt, h = 36)

elektronik_forecast_holt_manual <- forecast(object = elektronik_holt_manual, h = 36)

9.1.2 Visualisasi

Dari hasil forecast tersebut, langkah selanjutnya yaitu melakukan visualisasi dari hasil data train, data test dan hasil forecast dari holt winter.

elektronik_train %>% 
  autoplot() +
  autolayer(elektronik_test, series = "Data Test") + 
  autolayer(elektronik_forecast_holt$mean, series = "Model Forecast Auto") +
  autolayer(elektronik_forecast_holt_manual$mean, series = "Model Forecast Manual") +
  autolayer(elektronik_holt_manual$fitted[,1], series = "Holt Winter Manual")

Dari hasil visualisasi diatas, model sudah mendekat atau meniru dari data test.

9.1.3 Model Evaluation

Tahap selanjutnya, melakukan model evaluation dari model holt winter.

accuracy(elektronik_forecast_holt_manual$mean, elektronik_test)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -3.899256 5.434653 4.730196 -3.993342 4.687219 0.3765252  0.537377

Dari model evaluation yang saya dapatkan MAPE dari metode holtwinter sebesar 4.68% yang berarti model memiliki kecenderungan untuk salah prediksi sebesar 4.68%.

9.2 ARIMA

Pada metode kedua yaitu ARIMA. Menggunakan metode ARIMA ini melakukan pengujian menggunakan adf.test().

adf.test(elektronik_train)
## Warning in adf.test(elektronik_train): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  elektronik_train
## Dickey-Fuller = -5.7829, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

Menguji stasioneritas data H0: data tidak stasioner H1: data stasioner

pvalue < alpha artinya tolak h0 artinya dengan sudah stasioner.

9.2.1 Model Fitting

Tahap selanjutnya yaitu melakukan model fitting dengan melakukan tunning model untuk p,d dan q nya menggunakan acf dan pacf. ACF untuk melihat korelasi suatu data dengan mempertimbangkan data tengahnya. PACF untuk melihat korelasi data tanpa mempertimbangkan data tengahnya.

elektronik_train %>% 
  tsdisplay()

Kombinasi model ARIMA yang terbentuk :

ARIMA (p,d,q)

  • p : 1, 2, 3, 4, 5, 7, 8
  • d : 0
  • q : 0
  1. ARIMA (1,0,0)
  2. ARIMA (2,0,0)
  3. ARIMA (3,0,0)
  4. ARIMA (4,0,0)
  5. ARIMA (5,0,0)
  6. ARIMA (7,0,0)
  7. ARIMA (8,0,0)
elektronik_arima1 <- Arima(y = elektronik_ts,order = c(1,0,0))

elektronik_arima2 <- Arima(y = elektronik_ts,order = c(2,0,0))

elektronik_arima3 <- Arima(y = elektronik_ts,order = c(3,0,0))

elektronik_arima4 <- Arima(y = elektronik_ts,order = c(4,0,0))

elektronik_arima5 <- Arima(y = elektronik_ts,order = c(5,0,0))

elektronik_arima7 <- Arima(y = elektronik_ts,order = c(7,0,0))

elektronik_arima8 <- Arima(y = elektronik_ts,order = c(8,0,0))

9.2.2 Cek AIC

Setelah itu, cek AIC terlebih dahulu untuk masing-masing model.

elektronik_arima1$aic
## [1] 2737.192
elektronik_arima2$aic
## [1] 2641.803
elektronik_arima3$aic
## [1] 2461.876
elektronik_arima4$aic
## [1] 2309.696
elektronik_arima5$aic
## [1] 2226.276
elektronik_arima7$aic
## [1] 2225.666
elektronik_arima8$aic
## [1] 2219.857

Selanjutnya untuk membuat model arima secara otomatis menggunakan auto.arima()

elektronik_auto <- auto.arima(elektronik_ts)

elektronik_auto
## Series: elektronik_ts 
## ARIMA(2,1,1)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1
##       0.5503  -0.0683  -0.9477  -0.7635
## s.e.  0.0544   0.0549   0.0193   0.0331
## 
## sigma^2 estimated as 5.838:  log likelihood=-888.05
## AIC=1786.11   AICc=1786.27   BIC=1805.86

Model dengan menggunakan ARIMA yang saya dapatkan adalah elektronik.auto dengan AIC dengan masing-masing model yang paling kecil sebesar 1786.11.

9.2.3 Model Forecasting

Tahap selanjutnya yaitu dengan melakukan forecasting untuk 3 tahun ke depan

elektronik_forecast_arima <- forecast(object = elektronik_auto, h = 3*12)

9.2.4 Visualisasi

Tahap selanjutnya yaitu dengan melakukan visualisasi.

elektronik_train %>% 
  autoplot() +
  autolayer(elektronik_test, series = "Data Test") +
  autolayer(object = elektronik_auto$fitted, series = "Fitted") +
  autolayer(elektronik_forecast_arima$mean, series = "Arima Forecast")

Dari tahap diatas, maka data telah mampu mempelajari data testnya.

9.2.5 Model Evaluasi

Tahap selanjutnya yaitu model evaluasi dengan memeriksa nilai MAPE untuk nilai aic terkecil yaitu elektronik_auto.

accuracy(elektronik_auto)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.1474299 2.363849 1.791855 -0.2115265 1.942956 0.6278541
##                      ACF1
## Training set -0.002715831

Dari model evaluasi diatas, saya mendapatkan nilai MAPE sebesar 1.942% yang berarti model memiliki kecenderungan untuk salah prediksi sebesar 1.942%.

9.3 Seasonal ARIMA (SARIMA(p,d,q)(P,D,Q)[seasonal])

Dari data yang saya gunakan, hasil yang saya dapatkan adalah p-value < alpha dan tidak melakukan differencing karena data juga telah stationer.

Selanjutnya menentukkan nilai difference untuk index utama dan seasonal.

d = 0 D = 0

Untuk seasonal, AR murni * P : 1 # lag = 12 signifikan * D : 0 * Q : 0

Untuk keseluruhan data, MA murni

  • p : 1, 2, 3, 4, 5, 7, 8
  • d : 0
  • q : 0

Kombinasi model yang mungkin terbentuk adalah :

SARIMA : (p,d,q)(P,D,Q)[m]

  • SARIMA (1,0,0)(1,0,0)[12]
  • SARIMA (2,0,0)(1,0,0)[12]
  • SARIMA (3,0,0)(1,0,0)[12]
  • SARIMA (4,0,0)(1,0,0)[12]
  • SARIMA (5,0,0)(1,0,0)[12]
  • SARIMA (7,0,0)(1,0,0)[12]
  • SARIMA (8,0,0)(1,0,0)[12]

9.3.1 Model Fitting

Tahap selanjutnya melakukan model fitting dengan 7 model yang berbeda dan 1 model auto dengan menggunakan auto.arima().

elektronik_sarima1 <- Arima(y = elektronik_train, order = c(1,0,0), seasonal = c(1,0,0))

elektronik_sarima2 <- Arima(y = elektronik_train, order = c(2,0,0), seasonal = c(1,0,0))

elektronik_sarima3 <- Arima(y = elektronik_train, order = c(3,0,0), seasonal = c(1,0,0))

elektronik_sarima4 <- Arima(y = elektronik_train, order = c(4,0,0), seasonal = c(1,0,0))

elektronik_sarima5 <- Arima(y = elektronik_train, order = c(5,0,0), seasonal = c(1,0,0))

elektronik_sarima8 <- Arima(y = elektronik_train, order = c(8,0,0), seasonal = c(1,0,0))

elektronik_sarima_auto <- auto.arima(elektronik_train)

9.3.2 Cek AIC

Setelah itu, cek AIC terlebih dahulu untuk masing-masing model.

elektronik_sarima1$aic
## [1] 1779.861
elektronik_sarima2$aic
## [1] 1781.834
elektronik_sarima3$aic
## [1] 1770.673
elektronik_sarima4$aic
## [1] 1772.574
elektronik_sarima5$aic
## [1] 1769.704
elektronik_sarima8$aic
## [1] 1770.839
elektronik_sarima_auto$aic
## [1] 1569.832

Model dengan menggunakan SARIMA yang saya dapatkan adalah elektronik_sarima_auto dengan AIC dengan masing-masing model yang paling kecil sebesar 1569.832.

9.3.3 Cek Residual

Selanjutnya cek residual untuk memastikan model forecast tidak menghasilkan residual yang berkorelasi. Jika menghasilkan residual yang berkorelasi artinya masih terdapat informasi yang tertinggal.

elektronik_sarima_auto$residuals %>% 
  tsdisplay()

Dari hasil cek residual data diatas, data yang kurang dari 15 tidak ada residual yang saling berkorelasi.

9.3.4 Model Forecasting

Tahap selanjutnya yaitu dengan melakukan forecasting untuk 3 tahun ke depan

elektronik_forecast_sarima <- forecast(object = elektronik_sarima_auto, h = 3*12)

9.3.5 Visualisasi

Tahap selanjutnya yaitu dengan melakukan visualisasi.

elektronik_train %>% 
  autoplot() +
  autolayer(elektronik_test, series = "Data Test") +
  autolayer(object = elektronik_sarima_auto$fitted, series = "Fitted") +
  autolayer(elektronik_forecast_sarima$mean, series = "Sarima Forecast")

Dari tahap diatas, maka data telah mampu mempelajari data testnya.

9.3.6 Model Evaluasi

Tahap selanjutnya, dengan menggunakan model evaluasi dengan menggunakan accuracy() dan melihat nilai MAPE.

accuracy(elektronik_forecast_sarima$mean, elektronik_test)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -1.767665 4.112563 3.202789 -1.901711 3.127623 0.4039153 0.3886494

Dari model evaluasi diatas, saya mendapatkan nilai MAPE sebesar 3.127% yang berarti model memiliki kecenderungan untuk salah prediksi sebesar 3.127%.

9.4 Cek Asumsi

Tahap selanjutnya yaitu dengan cek asumsi pada time series di data yang saya gunakan untuk mengukur apakah residual sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria. Dibawah ini adalah asumsi untuk mengujinya.

9.4.1 No-autocorrelation residual

\(H_0\): residual tidak memiliki autocorrelation \(H_1\): residual mempunyai autocorrelation

yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

acf(elektronik_sarima_auto$residuals)

Box.test(elektronik_sarima_auto$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  elektronik_sarima_auto$residuals
## X-squared = 0.0064716, df = 1, p-value = 0.9359

\(H_0\): residual tidak memiliki autocorrelation \(H_1\): residual memiliki autocorrelation

Karena nilai pvalue > alpha, artinya residual tidak ada autokorelasi

9.4.2 Normality residual

\(H_0\): residual menyebar normal \(H_1\): residual tidak menyebar normal

yang diinginkan p-value > 0.05 (alpha), residual menyebar normal.

shapiro.test(elektronik_sarima_auto$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  elektronik_sarima_auto$residuals
## W = 0.98678, p-value = 0.002241

Dari hasil diatas, maka yang saya simpulkan data yang saya dapatkan residual tidak menyebar normal

hist(elektronik_sarima_auto$residuals)

10 Kesimpulan

Dari ketiga model yang saya buat, maka saya simpulkan hasil yang menghasilkan MAPE yang paling rendah yaitu dengan menggunakan metode ARIMA(). Dengan menggunakan ARIMA() menghasilkan nilai MAPE 1.942% yang berarti model memiliki kecenderungan untuk salah prediksi sebesar 1.942%