Padahal kesempatan ini, kita akan melakukan prediksi penjualan. Data yang akan digunakan adalah data penjualan yang terdapat di dalam package tsdl
Import Library
library(tidyverse) ## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.2 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(padr)## Warning: package 'padr' was built under R version 4.1.1library(forecast) ## Warning: package 'forecast' was built under R version 4.1.1## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(tseries) ## Warning: package 'tseries' was built under R version 4.1.1library(MLmetrics) ## Warning: package 'MLmetrics' was built under R version 4.1.1##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:base':
##
## Recalllibrary(imputeTS)## Warning: package 'imputeTS' was built under R version 4.1.1##
## Attaching package: 'imputeTS'## The following object is masked from 'package:tseries':
##
## na.removelibrary(tsdl)Memuat Data (Load Data)
Pada kesempatan kali ini, kita akan menggunakan sales dataset pertama yang ada di dalam package tsdl
devtools::install_github("FinYang/tsdl")## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install Rtools 4.0 from https://cran.r-project.org/bin/windows/Rtools/.## Skipping install of 'tsdl' from a github remote, the SHA1 (56e09154) has not changed since last install.
## Use `force = TRUE` to force installationsales <- subset(tsdl,"Sales")[[1]]
head(sales, 5*12)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1960 87695 86890 96442 98133 113615 123924 128924 134775 117357 114626
## 1961 92188 88591 98683 99207 125485 124677 132543 140735 124008 121194
## 1962 101007 94228 104255 106922 130621 125251 140318 146174 122318 128770
## 1963 108497 100482 106140 118581 132371 132042 151938 150997 130931 137018
## 1964 109894 106061 112539 125745 136251 140892 158390 148314 144148 140138
## Nov Dec
## 1960 107677 108087
## 1961 111634 111565
## 1962 117518 115492
## 1963 121271 123548
## 1964 124075 136485Eksplorasi Data (Explanatory Data Analysis)
Dataset yang disediakan oleh tsdl adalah dataset yang sudah dalam bentuk time series object (ts), sehingga kita tidak perlu lagi mengkonversinya ke dalam kelas ts
class(sales)## [1] "ts"Dalam analisis data time series, ketentuan yang perlu kita perhatikan dalam pembuatan model antara lain:
- Data harus urut berdasarkan periode waktunya
- Tidak boleh ada periode waktu yang terlewat
- Tidak boleh ada data yang hilang/kosong (missing value/NA)
Dataset yang disediakan oleh package tsdl ini adalah dataset yang telah bersih dan siap digunakan, sehingga pada ekslorasi data kali ini kita bisa langsung melihat bagaimana pola datanya.
sales %>%
autoplot()
Berdasarkan grafik tersebut, kita tahu bahwa dataset yang kita gunakan memiliki trend dan seasonal yang additive karena data bergerak secara konstan disekitar rata-rata-nya. Agar lebih terlihat bagaimana error, trend dan seasonal-nya, kita akan lakukan decompose pada dataset kita
sales %>%
decompose() %>%
autoplot()
Dari grafik yang ditampilkan, dapat kita lihat bahwa dataset kita memiliki trend yang meningkat dan pola seasonal tahunan.
Pembuatan Model
Cross Validation
Sebelum membuat model, kita akan membagi dataset kita menjadi data train dan data test. Karena kita hanya akan memprediksi penjualan satu tahun ke depan, maka data test yang akan kita gunakan juga data penjualan selama satu tahun terakhir. Sedangkan sisa dari data test akan digunakan sebagai data train
sales_test <- tail(sales, 1*12)
sales_train <- head(sales, length(sales) - length(sales_test))Model Holt’s Winter
Seperti yang telah kita ketahui bahwa data kita memiliki trend dan seasonal, oleh karena itu model yang akan kita gunakan disini Holt’s Winter atau triple smoothing exponential. Sebagai pembanding, kita akan membuat juga model Holt’s Winter yang otomatis ditentukan oleh R.
hw_manual <- ets(sales_train, model = "AAA")
hw_auto <- ets(sales_train, model = "ZZZ")
hw_manual## ETS(A,A,A)
##
## Call:
## ets(y = sales_train, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.1151
## beta = 0.0047
## gamma = 1e-04
##
## Initial states:
## l = 108151.6497
## b = 410.1542
## s = -3828.341 -4532.028 7198.531 6026.774 25452.78 23065.42
## 10417.93 9748.918 -11127.96 -15004.95 -26057.41 -21359.66
##
## sigma: 6058.866
##
## AIC AICc BIC
## 4087.316 4091.094 4141.596hw_auto## ETS(M,A,A)
##
## Call:
## ets(y = sales_train, model = "ZZZ")
##
## Smoothing parameters:
## alpha = 0.0949
## beta = 0.0047
## gamma = 1e-04
##
## Initial states:
## l = 108151.8419
## b = 430.8957
## s = -3828.366 -4531.961 7198.939 6026.738 25452.68 23065.48
## 10418.04 9749.063 -11128.3 -15005.03 -25862.1 -21555.2
##
## sigma: 0.0351
##
## AIC AICc BIC
## 4043.173 4046.951 4097.453Dari model hw_manual dan hw_auto nilai betha dan gamma yang dihasilkan sama, namun untuk nilai alpha kedua model berbeda, dimana model hw_manual mempunyai alpha sebesar 0.1151 dan model hw_auto mempunyai alpha sebesar 0.0949.
Berdasarkan perbedaan tersebut, AIC yang dihasilkan oleh kedua model juga berbeda. Model hw_auto menghasilkan AIC yang lebih kecil, yaitu 4043.173 dibanding AIC hw_manual yang AIC-nya sebesar 4087.316. Dengan demikian dapat dikatakan model hw_auto lebih baik dibanding model hw_manual.
Model SARIMA
Selain model Holt’s Winter kita juga akan membuat model dengan menggunakan SARIMA. Untuk membuat model ARIMA terlebih dahulu kita cek apakah data kita sudah stasioner atau belum
sales_train %>%
adf.test()## Warning in adf.test(.): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: .
## Dickey-Fuller = -9.4641, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryDengan menggunakan Dickey-Fuller test, p-value yang diperoleh adalah 0.01, dimana ini lebih kecil dari alpha (0.05), yang berarti data yang akan kita gunakan sudah stasioner, sehingga tidak perlu dilakukan differencing. Karena tidak ada differencing, maka nilai d dan D kita adalah 0
Selanjutnya akan kita lihat grafik ACF dan PACF kita untuk menentukan nilai p, q, P, dan Q
sales_train %>%
tsdisplay()
Untuk keseluruhan data, grafik ACF kita adalah dies down sedangkan grafik PACF kita adalah cut off pada lag ke-1, ke-3, ke-4, ke-5, ke-7, ke-9, ke-11, ke-13, ke-14 dan ke-29. Namun untuk pembuatan model kita hanya akan coba untuk menggunakan lag ke-1, lag ke-3, dan lag ke-4.
Dengan demikian, model AR kita adalah sebagai berikut:
- p : 1, 3, 4
- d : 0
- q : 0
Sementara untuk data seasonal kita, baik grafik ACF maupun PACF keduanya memiliki pola yang dies down.
Dengan demikian, model MA kita adalah sebagai berikut:
- P : 0
- D : 0
- Q : 0
Sama seperti Holt’s Winter, kita juga akan membuat model SARIMA yang otomatis ditentukan oleh R sebagai pembanding
sarima_1 <- Arima(y = sales_train, order = c(1,0,0), seasonal = c(0,0,0))
sarima_2 <- Arima(y = sales_train, order = c(3,0,0), seasonal = c(0,0,0))
sarima_3 <- Arima(y = sales_train, order = c(4,0,0), seasonal = c(0,0,0))
sarima_auto <- auto.arima(sales_train, seasonal = T)
sarima_1## Series: sales_train
## ARIMA(1,0,0) with non-zero mean
##
## Coefficients:
## ar1 mean
## 0.9390 157074.06
## s.e. 0.0264 16171.64
##
## sigma^2 estimated as 207603624: log likelihood=-1979.07
## AIC=3964.15 AICc=3964.28 BIC=3973.72sarima_2## Series: sales_train
## ARIMA(3,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 mean
## 0.8245 0.2822 -0.1682 156751.15
## s.e. 0.0731 0.0936 0.0746 15715.38
##
## sigma^2 estimated as 199508772: log likelihood=-1974.54
## AIC=3959.07 AICc=3959.42 BIC=3975.04sarima_3## Series: sales_train
## ARIMA(4,0,0) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 mean
## 0.7840 0.3487 0.0347 -0.2482 157061.99
## s.e. 0.0717 0.0930 0.0940 0.0734 12012.95
##
## sigma^2 estimated as 188473871: log likelihood=-1969.01
## AIC=3950.02 AICc=3950.5 BIC=3969.17sarima_auto## Series: sales_train
## ARIMA(2,1,0)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 sma1 sma2
## -0.9922 -0.6617 -0.6097 -0.1424
## s.e. 0.0577 0.0595 0.0932 0.0875
##
## sigma^2 estimated as 30226168: log likelihood=-1678.43
## AIC=3366.85 AICc=3367.22 BIC=3382.44Pemilihan Model
Dari semua model yang ada akan kita seleksi untuk kita gunakan sebagai model untuk memprediksi penjualan satu tahun ke depan. Hal yang menjadi acuan untuk pemilihan model adalah nilai AIC dari setiap model. Model dengan nilai terkecil lah yang akan kita pilih.
print(paste("AIC hw_manual: ", hw_manual$aic))## [1] "AIC hw_manual: 4087.31600835051"print(paste("AIC hw_auto: ", hw_auto$aic))## [1] "AIC hw_auto: 4043.17302289997"print(paste("AIC sarima1: ", sarima_1$aic))## [1] "AIC sarima1: 3964.14601258481"print(paste("AIC sarima2: ", sarima_2$aic))## [1] "AIC sarima2: 3959.07166889384"print(paste("AIC sarima2: ", sarima_3$aic))## [1] "AIC sarima2: 3950.01636930519"print(paste("AIC sarima_auto: ", sarima_auto$aic))## [1] "AIC sarima_auto: 3366.85087666638"Berdasarkan nilai AIC-nya, model sarima auto adalah model yang memiliki nilai AIC terkecil. Dengan demikian model sarima_auto lah yang akan kita gunakan untuk memprediksi penjualan setahun ke depan.
Peramalan Penjualan
Selanjutnya kita akan melakukan peramalan penjualan selama satu tahun ke depan dengan menggunakan model sarima_auto
sales_forecast <- forecast(sarima_auto, h = 1*12)
sales_forecast## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 1975 199369.6 192322.6 206416.7 188592.1 210147.2
## Feb 1975 190906.7 183859.4 197954.0 180128.8 201684.6
## Mar 1975 205527.0 198104.7 212949.2 194175.6 216878.3
## Apr 1975 212611.3 203826.6 221396.1 199176.2 226046.4
## May 1975 238528.7 229703.6 247353.8 225031.8 252025.5
## Jun 1975 233450.3 224096.3 242804.3 219144.6 247756.0
## Jul 1975 250963.6 241006.6 260920.5 235735.7 266191.4
## Aug 1975 254216.0 244129.9 264302.2 238790.6 269641.5
## Sep 1975 223614.0 213034.8 234193.2 207434.5 239793.5
## Oct 1975 235537.7 224590.5 246485.0 218795.3 252280.1
## Nov 1975 221058.7 209903.9 232213.5 203998.9 238118.5
## Dec 1975 223192.7 211626.1 234759.3 205503.1 240882.3sales_train %>%
autoplot() +
autolayer(sales_test, series = "Data Tes") +
autolayer(sales_forecast$mean, series = "Peramalan") +
autolayer(sarima_auto$fitted, series = "Model")
Jika kita perhatikan, model yang telah kita buat sudah mampu menangkap pola data dengan sangat baik, baik itu di data train maupun di data test.
Evaluasi Model
Untuk evaluasi model, kita akan coba cek berapa MAPE yang dihasilkan oleh model.
accuracy(sales_forecast$mean, sales_test)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 1239.729 6936.803 5683.172 0.4176778 2.474218 -0.1439487 0.4704209Dapat kita lihat model menghasilkan MAPE yang cukup kecil, yaitu 2.47%.
Selain dengan menggunakan nilai MAPE, kita juga perlu mengecek autokorelasi dan normalitas model.
Untuk menguji normalitas model, kita akan menggunakan Shapiro test
shapiro.test(sarima_auto$residuals)##
## Shapiro-Wilk normality test
##
## data: sarima_auto$residuals
## W = 0.98363, p-value = 0.03337Nilai p-value model sebesar 0.03337, lebih kecil dari alpha (0.05), yang berarti error (residual) dari model telah terdistribusi secara normal. Adapun distribusi dari error model adalah sebagai berikut
hist(sarima_auto$residuals, breaks = 20)
Sementara untuk autocorrelation, kita akan mengujinya dengan menggunakan uji L-Jung
Box.test(sarima_auto$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: sarima_auto$residuals
## X-squared = 0.12817, df = 1, p-value = 0.7203Berdasarkan uji L-Jung, model mempunyai p-value yang lebih besar dari aplha (0.05), yaitu 0.7203, yang berarti residual (error) dari modal tidak terjadi autokorelasi.
Kesimpulan
Dari hasil evaluasi model yang telah dilakukan, dapat disimpulkan bahwa model Auto Sarima (sarima_auto) merupakan model yang terbaik dengan nilai AIC sebesar 3366.85 dan MAPE sebesar 2.47%. Model juga sudah baik karena residual sudah terdistribusi secara normal dan tidak mempunyai autokorelasi