Import Library

# Library for data Processing
library(tidyverse)
# Library to forecast
library(forecast)
# Library for Processing date data type
library(lubridate)
# Library Simple Moving Average
library(TTR)
# Library MAPE
library(MLmetrics)
# Library for Processing date data type
library(zoo)
# Data visualisation
library(plotly)
library(xts)
library(TSstudio)
library(tseries)

Read Data

# install.packages("devtools")
devtools::install_github("FinYang/tsdl")
library(tsdl)
tsdl
## Time Series Data Library: 648 time series  
## 
##                        Frequency
## Subject                 0.1 0.25   1   4   5   6  12  13  52 365 Total
##   Agriculture             0    0  37   0   0   0   3   0   0   0    40
##   Chemistry               0    0   8   0   0   0   0   0   0   0     8
##   Computing               0    0   6   0   0   0   0   0   0   0     6
##   Crime                   0    0   1   0   0   0   2   1   0   0     4
##   Demography              1    0   9   2   0   0   3   0   0   2    17
##   Ecology                 0    0  23   0   0   0   0   0   0   0    23
##   Finance                 0    0  23   5   0   0  20   0   2   1    51
##   Health                  0    0   8   0   0   0   6   0   1   0    15
##   Hydrology               0    0  42   0   0   0  78   1   0   6   127
##   Industry                0    0   9   0   0   0   2   0   1   0    12
##   Labour market           0    0   3   4   0   0  17   0   0   0    24
##   Macroeconomic           0    0  18  33   0   0   5   0   0   0    56
##   Meteorology             0    0  18   0   0   0  17   0   0  12    47
##   Microeconomic           0    0  27   1   0   0   7   0   1   0    36
##   Miscellaneous           0    0   4   0   1   1   3   0   1   0    10
##   Physics                 0    0  12   0   0   0   4   0   0   0    16
##   Production              0    0   4  14   0   0  28   1   1   0    48
##   Sales                   0    0  10   3   0   0  24   0   9   0    46
##   Sport                   0    1   1   0   0   0   0   0   0   0     2
##   Transport and tourism   0    0   1   1   0   0  12   0   0   0    14
##   Tree-rings              0    0  34   0   0   0   1   0   0   0    35
##   Utilities               0    0   2   1   0   0   8   0   0   0    11
##   Total                   1    1 300  64   1   1 240   3  16  21   648

Dalam kasus kali ini saya akan menggunakan data sales. karena data pada tsdl sudah berbendtuk timeseries maka kita tidak perlu mengubahnya ke timeseries.

sales <- tsdl[[4]]
head(sales)
##         Jan    Feb    Mar    Apr    May    Jun
## 1960  87695  86890  96442  98133 113615 123924

Data Wrangling

Pertama kita lihat dulu jumlah data timeseries kita

length(sales)
## [1] 192
# NA check
anyNA(sales)
## [1] FALSE

Kita bisa melihat panjang total data 192 dan didapat FALSE, artinya data tidak memiliki NA. Tidak hanya NA yang harus kita periksa, kita juga harus mencoba mengetahui bagaimana distribusi dari data timeseries kami. kita dapat menggunakan sintaks summary untuk mengetahui distribusi data kita.

summary(sales)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   86890  128426  157459  162064  193556  255918

Minimal pada data time series kita adalah 86890, berarti data kita tidak memiliki angka 0.

Exploratory Data Analysis

Pada sub bagian ini, saya akan mencoba mengeksplorasi objek time series kita. Kita akan mencoba untuk lebih memahami objek kami, apakah kita dapat menemukan trend, seasonality, atau interpretasi atau wawasan lain berdasarkan visualisasi.

Decompose

Pada bagian ini kita akan mencoba mendekomposisi data, yang artinya kita akan mencoba memisahkan 3 komponen: Trend, Seasonality dan Error.

sales_deco <- decompose(sales, type = "multiplicative")

Visualize trend data

sales_deco$trend %>% autoplot()

Kami menemukan bahwa trend time series sudah smooth dan tidak ada seasonality lainnya di garis tersebut, itu berarti setup seasonality dan frequency di objek ts telah sesuai dengan baik.

selanjutnya, mari kita visualisasikan seasonality-nya

sales_deco$seasonal %>% autoplot()

bisa kita bilang bahwa object decompose kita bisa menangkap seasonality dengan baik karena dapat dilihat dalam grafik di atas terdapat pola berulang dengan periode yang tetap/sama.

untuk memastikannya, kita dapat melakukan decompose

autoplot(sales_deco)

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

dari adf.test terlihat bahwa nilai error sangatlah kecil. ini artinya adalah objek time series membentuk multiplicative yang berarti kenaikannya bertambah secara konstan.

sales %>% 
  diff() %>% 
  plot()

dari hasil tersebut terlihat bahwa objek time series memiliki komponen tren dengan pola multiplicative. kemudian dari hasil test didapat nilai p-value < 0.01. jika asumsi Ho data tidak stasioner dengan p-value > 0.05, maka tolak Ho, dimana data tersebut sudah stasioner. namun untuk keperluan permodelan time series perlu dilakukan differencing supaya tidak membentuk tren, dari plot diatas hasil 1x differencing data yang stasioner yang tidak membentuk tren

Cross Validation

untuk dapat membuat model yang baik, saya akan split data menjadi dua yaitu data train (data 13 tahun dari awal) dan data validation (data 2 tahun terakhir)

sales_train <- head(sales, length(sales) - 24)
sales_test <- tail(sales, 24)
autoplot(sales_train)+geom_line(data = sales_test, col = "tomato3")

Forecasting

setelah memiliki data train, maka saya akan membuat mdoel dengan menggunakan mode auto untuk mencari nilai p, d dan q

auto_sales <- auto.arima(sales_train, seasonal = F)
auto_sales
## Series: sales_train 
## ARIMA(5,1,1) with drift 
## 
## Coefficients:
##          ar1     ar2     ar3      ar4      ar5      ma1     drift
##       0.3313  0.3530  0.0740  -0.4880  -0.1942  -0.9136  676.7706
## s.e.  0.0791  0.0703  0.0751   0.0715   0.0810   0.0288   74.4922
## 
## sigma^2 estimated as 92839723:  log likelihood=-1767.39
## AIC=3550.78   AICc=3551.69   BIC=3575.73

dai hasil auto tersebut, didap bahwa model forecastnya adalah ARIMA(5,1,1). selanjutnya saya akan tuning pada time series.

tsdisplay(diff(sales_train))

terlihat bahwa plot ACF dan PACF adalah cut off lag, dimana dapat disimpulkan hasilnya untuk nilai masing-masing p,d,q

auto.arima (5,1,1) ACF = 2,4,6,8 PACF = 4,6,7,8

auto_sales_1 <- arima(sales_train, order = c(5,1,1)) # model ARIMA 5,1,1
auto_sales_2 <- arima(sales_train, order = c(2,1,1)) # model ARIMA 2,1,1 
auto_sales_3 <- arima(sales_train, order = c(4,1,1)) # model ARIMA 4,1,1
auto_sales_4 <- arima(sales_train, order = c(6,1,1)) # model ARIMA 6,1,1
auto_sales_5 <- arima(sales_train, order = c(8,1,1)) # model ARIMA 8,1,1
tsdisplay(auto_sales_1$residuals)

lalu, saya akan mencoba lakukan forecast selama 2 tahun berikutnya (24 bulan)

forecast1 <- forecast(auto_sales_1, h = 24)
forecast2 <- forecast(auto_sales_2, h = 24)
forecast3 <- forecast(auto_sales_3, h = 24)
forecast4 <- forecast(auto_sales_4, h = 24)
forecast5 <- forecast(auto_sales_5, h = 24)

lalu, saya akan membandingkan dari masing-masing forecasting untuk dicari nilai MAPE terhadap data validation

accuracy(forecast1, sales_test)
##                    ME     RMSE       MAE      MPE     MAPE      MASE       ACF1
## Training set 2855.153 10163.27  8534.194 1.488787 5.655918 0.9946131 -0.1512301
## Test set     4548.858 18623.23 16002.706 1.278122 7.078067 1.8650267  0.4005099
##              Theil's U
## Training set        NA
## Test set     0.9222993
accuracy(forecast2, sales_test)
##                      ME     RMSE      MAE       MPE     MAPE     MASE
## Training set   624.3751 13154.40 10708.32 0.1252036 7.101198 1.247995
## Test set     16353.1407 27394.95 22711.55 6.4402026 9.728148 2.646905
##                     ACF1 Theil's U
## Training set -0.01890237        NA
## Test set      0.54755395  1.341255
accuracy(forecast3, sales_test)
##                    ME     RMSE       MAE      MPE    MAPE     MASE       ACF1
## Training set 2750.872 10273.66  8634.079 1.420179 5.73474 1.006254 -0.1535234
## Test set     4919.772 19978.78 17098.020 1.385770 7.56048 1.992680  0.4198538
##              Theil's U
## Training set        NA
## Test set     0.9834761
accuracy(forecast4, sales_test)
##                    ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 2937.021  9039.445  7519.754 1.6119968 4.883666 0.8763857
## Test set     2722.246 16441.765 14524.122 0.5746104 6.467719 1.6927059
##                    ACF1 Theil's U
## Training set -0.1488308        NA
## Test set      0.1497860 0.8210069
accuracy(forecast5, sales_test)
##                    ME      RMSE       MAE       MPE     MAPE      MASE
## Training set 3125.053  8709.606  7214.865 1.7602088 4.630672 0.8408525
## Test set     2402.380 14785.429 12858.134 0.5217345 5.732507 1.4985442
##                    ACF1 Theil's U
## Training set -0.1392017        NA
## Test set      0.1048853 0.7392242
autoplot(forecast5$x)+
  autolayer(forecast5$fitted)+
  autolayer(forecast5$mean)+geom_line(data = sales_test, col = "tomato3")

Asumption Check

Setelah kita mendapatkan hasil dari masing-masing model. Mari kita coba untuk uji asumsi.

Normality Residuals

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

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

Untuk mengecek normality residual pada hasil forecasting time series kita bisa melakukan uji normality (shapiro test) dengan menggunakan fungsi shapiro.test(residual model)

#model 1
shapiro.test(forecast1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast1$residuals
## W = 0.98524, p-value = 0.07282
#model 2
shapiro.test(forecast2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast2$residuals
## W = 0.99433, p-value = 0.7665
#model 3
shapiro.test(forecast3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast3$residuals
## W = 0.98658, p-value = 0.1078
#model 4
shapiro.test(forecast4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast4$residuals
## W = 0.98493, p-value = 0.06651
#model 5
shapiro.test(forecast5$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  forecast5$residuals
## W = 0.99042, p-value = 0.319

Dari 5 model yang kita cek, semua model memiliki p-value lebih dari 0.05, dapat kita artikan bahwa residual semua model menyebar secar normal.

No-autocorrelation residual

\(H_0\): residual has no-autocorrelation \(H_1\): residual has autocorrelation

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

Mari kita uji Ljung-box dengan menggunakan fungsi Box.test(residual model, type = "Ljung-Box)

#model1
Box.test(forecast1$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast1$residuals
## X-squared = 3.9113, df = 1, p-value = 0.04796
#model2
Box.test(forecast2$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast2$residuals
## X-squared = 0.061105, df = 1, p-value = 0.8048
#model3
Box.test(forecast3$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast3$residuals
## X-squared = 4.0308, df = 1, p-value = 0.04468
#model4
Box.test(forecast4$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast4$residuals
## X-squared = 3.7881, df = 1, p-value = 0.05162
#model5
Box.test(forecast5$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  forecast5$residuals
## X-squared = 3.3138, df = 1, p-value = 0.0687

Dari hasil test di atas, hanya model 2, 4, dan 5 lolos dari uji autokorelasi karena p-valuenya lebih dari 0.05.

Conclusion

dari hasil plot forecast pada model terbaik adalah model yang memiliki MAPE terkecil yaitu pada model forecast 5 dengan angka 5.73%. Model forecast 5 juga lolos dari semua uji asumsi, jadi bisa kita simpulkan bahwa model forecast 5 adalah model yang cocok untuk memprediksi dataset sales kita.