Pendahuluan

Forecasting sangatlah berharga bagi bisnis karena memberikan kemampuan untuk membuat keputusan bisnis yang terinformasi dan mengembangkan strategi berbasis data. Perusahaan pun dapat membuat keputusan keuangan dan operasional berdasarkan kondisi pasar, sehingga prediksi masa depan juga bisa terlihat.

Hal ini karena forecasting dilakukan dengan mengumpulkan data di masa lalu dan dianalisis untuk menemukan pola, yang pada akhirnya digunakan untuk memprediksi tren dan perubahan di masa depan.

> Maksud dan Tujuan

Penelitian ini bertujuan untuk membuat model machine learning yang dapat memprediksi total quantity harian dari product yang terjual pada studi kasus FMCG dalam kurun waktu 1 tahun.

> Benefit

Hasil forecasting yang dilakukan dapat menjadi salah satu acuan dalam membuat kebijakan mengenai manajemen stok product.

> Batasan

Pada penelitian ini masih terdapat permasalahan outlier dan memiliki pola seasonal yang cukup rumit.

Data Wrangling

1. Import Library

Sebelum melakukan analisis kita perlu memuat library-library yang dibutuhkan dalam analisa ini.

Install terlebih dahulu jika belum pernah menggunakan sebelumnya.

# Import library untuk data wrangling
library(tidyverse)
library(dplyr)
library(DescTools)
library(padr)

# Import library untuk visualisasi
library(plotly) # Interactive Plotting
library(ggplot2)
library(TSstudio)

# Machine Learning - Time Series 
library(forecast)
library(tseries) 
library(MLmetrics)

2. Read Data

Data pada penelitian ini berisikan informasi mengenai Usia, Gender, Status, Income, dan Riwayat Transaksi, Quantity, serta Amount pada tahun 2022.

data_prod <- read.csv("Data/dataset-ts.csv", sep = ";") %>%
  janitor::clean_names()

# Cek struktur data
glimpse(data_prod)

## Rows: 365
## Columns: 2
## $ date           <chr> "2022-01-01", "2022-01-02", "2022-01-03", "2022-01-04",…
## $ total_quantity <int> 49, 59, 63, 45, 30, 71, 37, 47, 47, 60, 75, 35, 28, 31,…

Definisi Variabel:

date: Tanggal transaksi penjualan
total_quantity: Jumlah Quantity produk yang terjual

3. Data Cleansing

Dibagian ini variabel data akan disesuaikan tipe datanya.

> Mengubah tipe variabel

data_prod <- data_prod %>% 
  # Ubah type data yang belum sesuai
  mutate(date = as.Date(date)
         )

glimpse(data_prod)

## Rows: 365
## Columns: 2
## $ date           <date> 2022-01-01, 2022-01-02, 2022-01-03, 2022-01-04, 2022-0…
## $ total_quantity <int> 49, 59, 63, 45, 30, 71, 37, 47, 47, 60, 75, 35, 28, 31,…

> Karakteristik data time series

Data harus terurut berdasarkan periode waktunya

data_prod <- data_prod %>%
  arrange(date)

Tidak boleh ada waktu atau periode yang terlewat/bolong

data_prod <- data_prod %>% 
  pad()

## pad applied on the interval: day

range(data_prod$date)

## [1] "2022-01-01" "2022-12-31"

nrow(data_prod)

## [1] 365

Data yang digunakan ada sebanyak 365 hari, dimulai dari tanggal 1 Januari 2022 hingga 31 Desember 2022.

Data tidak boleh ada yang missing

data_prod %>% 
  is.na() %>% 
  colSums()

##           date total_quantity 
##              0              0

Data yang digunakan sudah tidak memiliki missing value, sehingga dapat dilanjutkan ke tahap selanjutnya.

Eksplorasi Data

Ekplorasi data dilakukan agar kita menjadi lebih paham dengan gambaran awal dari data yang akan digunakan.

> Trend data

data_prod %>% 
  ts_plot(title = "Trend Quantity Penjualan Produk Selama Tahun 2022",
          Xtitle = "Date", 
          Ytitle = "Total Quantity",
          line.mode = "lines+markers")

Insight :

Seperti yang terlihat dari hasil plot bahwa: 1. Total Quantity Produk sangat berfluktuatif, dan secara sekilas mayoritas bergerak diantara 30 hingga 80. 2. Terdapat data yang jauh dari pergerakan data lainnya, seperti pada 8 Maret total quantity hanya 16 dan 2 Maret total quantity mencapai 119.

> Dekomposisi

Tujuan dari melakukan dekomposisi adalah untuk melihat apakah dalam data terdapat trend dan seasonal. Hal ini untuk memastikan pola trend dan seasonal sudah ditangkap dengan baik dan frequency yang kita gunakan sudah tepat.

Decompose data akan memecah data ts menjadi 3 komponen (trend, seasonality, error). * Trend: Pola data secara umum, apakah mengalami kenaikan atau penurunan * Seasonality: Pola berulang dengan periode waktu yang tetap * Remainder/Error: Komponen lain yang tidak tertangkap oleh Trend dan Seasonality.

data_prod %>% 
  .$total_quantity %>% 
  ts(frequency = 7) %>% 
  decompose() %>% 
  autoplot() +
  labs(title = "Dekomposisi Mingguan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Pada bagian ini data masih terlihat berfluktuasi sehingga trend nya belum terlihat jelas.
Seasonal: Terlihat bahwa data memiliki pola mingguan, karena terlihat pola berulang pada seasonal.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

data_prod %>% 
  .$total_quantity %>% 
  ts(frequency = 7*4) %>% 
  decompose() %>% 
  autoplot() +
  labs(title = "Dekomposisi Bulanan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Terlihat data memiliki pola menurun, namun masih terlihat berfluktuasi.
Seasonal: Terlihat bahwa data memiliki pola mingguan, karena terlihat pola berulang pada seasonal.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

data_prod %>% 
  .$total_quantity %>% 
  msts(seasonal.periods = c(7, 7*4)) %>% 
  mstl() %>% 
  autoplot() +
  labs(title = "Dekomposisi Mingguan + Bulanan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Terlihat data memiliki pola menurun dan fluktuasi menjadi lebih smooth.
Seasonal: Terlihat bahwa data memiliki pola bulanan, karena terlihat pola berulang pada seasonal28.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

Modelling: Forecasting menggunakan Metode ARIMA

Dalam pembangunan model akan digunakan beberapa metode, diantaranya ARIMA dan Seasonal ARIMA (SARIMA). Model tersebut dipilih karena pada data yang digunakan terdapat trend dan seasonal.

> Cross Validation Dataset

Dilakukan cross validation untuk membagi dataset menjadi data train/data latih dan data test. Data test bertujuan untuk mengetahui seberapa baik model dalam memprediksi data yang belum pernah dilihat dalam pelatihan (generalisasi).

Dalam kasus ini data test diambil sebanyak 85 sampel berhubung keseluruhan data hanya berjumlah 365.

n_train <- 10*4*7
n_test <- nrow(data_prod) - n_train

data_prod_ts1 <- ts(data_prod[,2], start = 1, frequency = 1)
data_prod_ts7 <- ts(data_prod[,2], start = 1, frequency = 7)
data_prod_ts28 <- ts(data_prod[,2], start = 1, frequency = 7*4)

data_train_s1 <- head(x = data_prod_ts1, n = n_train)
data_test_s1 <- tail(x = data_prod_ts1, n = n_test)

data_train_s7 <- head(x = data_prod_ts7, n = n_train)
data_test_s7 <- tail(x = data_prod_ts7, n = n_test)

data_train_s28 <- head(x = data_prod_ts28, n = n_train)
data_test_s28 <- tail(x = data_prod_ts28, n = n_test)

data_prod_ms <- msts(data = data_prod[,2], seasonal.periods = c(7, 7*4, 7*4*2, 7*4*3))

data_train_ms <- head(x = data_prod_ms, n = n_train)
data_test_ms <- tail(x = data_prod_ms, n = n_test)

> Stasionarity

Stasionarity time series memiliki arti bahwa pada data time series yang kita miliki tidak memiliki trend maupun seasonal dan memiliki variansi konstan. Melakukan pengujian menggunakan adf.test()

# Augmented Dickey-Fuller Test
data_train_s1 %>% 
  adf.test()

## Warning in adf.test(.): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -5.8487, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

Menguji stasioneritas data H0: data tidak stasioner H1: data stasioner p-value < alpha (0.05), artinya data stasioner.

Apabila data yang kita miliki belum stasioner, maka kita bisa melakukan differencing.

> Differencing

Jadi untuk membuat datanya stasioner, cara yang paling umum digunakan adalah dengan melakukan differencing, yaitu mengurangi data saat ini dengan data sebelumnya. Terkadang, tergantung pada kompleksitas data, jumlah differencing bisa lebih dari 1 kali.

Namun data yang dimiliki sudah stasioner, sehingga tidak perlu dilakukan differencing. Selanjutnya akan dilakukan pemodelan.

1. ARIMA (Autoregressive Integrated Moving Average)

ARIMA adalah gabungan antara dua metode, yaitu Auto Regressive (AR) dan Moving Average (MA). I nya menjelaskan Integrated. Tujuan utama dari ARIMA adalah melakukan autocorrelation pada data.

# create model
auto_arima <- auto.arima(y = data_train_s1)

auto_arima_f <- forecast(auto_arima, h = 85)

# Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_prod_ts1, # data aktual full
              forecast.obj = auto_arima_f, # hasil forecast
              test = data_test_s1) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA tanpa seasonal adalah ARIMA(0,0,0).
Terlihat bahwa hasil forecasting membentuk garis lurus di rata-rata data.
Model belum bisa memprediksi fluktuasi data.

2. SARIMA (Seasonal ARIMA)

Pada metode ARIMA dapat pula diikut sertakan pola seasonal yang dimiliki oleh data, misalnya ingin melihat pengaruh pola mingguan atau bulanan, dst.

auto_sarima7 <- stlm(y = data_train_s7, method = "arima", s.window = 7)

auto_sarima7_f <- forecast(auto_sarima7, h = 85)

# Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_prod_ts7, # data aktual full
              forecast.obj = auto_sarima7_f, # hasil forecast
              test = data_test_s7) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA adalah ARIMA(1,0,1), dimana seasonal yang digunakan adalah pola mingguan.
Terlihat bahwa hasil forecasting sudah bisa mengikuti fluktuasi dari data meskipun masih cukup jauh dari data aktualnya.

auto_sarima28 <- stlm(y = data_train_s28, method = "arima", s.window = 28)

auto_sarima28_f <- forecast(auto_sarima28, h = 85)

# Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_prod_ts28, # data aktual full
              forecast.obj = auto_sarima28_f, # hasil forecast
              test = data_test_s28) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA adalah ARIMA(0,0,0), dimana seasonal yang digunakan adalah pola bulanan.
Terlihat bahwa hasil forecasting bisa mengikuti fluktuasi dari data meskipun masih cukup jauh dari data aktualnya.

3. Multiple Seasonal

Selain menggunakan 1 pola seasonal yang dimiliki oleh data, misalnya ingin melihat pengaruh pola mingguan atau bulanan, dst. Dapat pula menyertakan beberapa pola seasonal secara sekaligus untuk melihat pengaruh pola seasonalnya secara bersamaan.

stlm_ms <- stlm(y = data_train_ms, method = "arima", allow.multiplicative.trend = TRUE, s.window = c(7, 7*4, 7*4*2, 7*4*3))

stlm_ms_f <- forecast(stlm_ms, h = 85)

# Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_prod_ms, # data aktual full
              forecast.obj = stlm_ms_f, # hasil forecast
              test = data_test_ms) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA adalah ARIMA(3,0,5), dimana seasonal yang digunakan adalah pola mingguan, 1 bulanan, 2 bulan, dan 3 bulan secara bersamaan.
Terlihat bahwa hasil forecasting sudah bisa mengikuti fluktuasi dari data dan sudah mendekati data aktualnya.

> Evaluasi Model

eval <-
  rbind(
    accuracy(auto_arima$fitted, data_train_s1),
    accuracy(auto_sarima7$fitted, data_train_s7),
    accuracy(auto_sarima28$fitted, data_train_s28),
    accuracy(stlm_ms$fitted, data_train_ms)
    )

rownames(eval) <- c("ARIMA","SARIMA(Seasonal 7)","SARIMA(Seasonal 28)", "Multiple Seasonal ARIMA")
round(eval,2)

##                            ME  RMSE   MAE    MPE  MAPE  ACF1 Theil's U
## ARIMA                    0.00 17.27 13.44 -13.86 32.10 -0.02      0.61
## SARIMA(Seasonal 7)      -0.01 14.70 11.48 -10.59 26.88  0.06      0.53
## SARIMA(Seasonal 28)      0.00 16.39 12.91 -12.46 30.41 -0.04      0.59
## Multiple Seasonal ARIMA -0.04  9.53  7.67  -4.02 16.98 -0.03      0.33

Insight :

Dari hasil evaluasi di atas, dapat disimpulkan bahwa model telah dapat memprediksi total quantity produk. Adapun model terbaik adalah metode ARIMA dengan Multiple Seasonal dengan nilai error terkecil. Dimana nilai MAE-nya adalah 7.67, berarti rentang kesalahan data prediksi terhadap data aktual senilai (+/-) 7.67.

Improving Model

Dalam kasus ini, Model ARIMA dengan Multiple Seasonal menjadi model terbaik. Maka dari itu dicoba dilakukan treatment dalam meningkatkan performa model. Adapun treatment yang akan dilakukan yaitu dengan menangani Outlier.

Pola persebaran data total quantity:

hist(data_train_ms, breaks = 100, col = "navy", border = "darkred", main = "Sebaran Data Total Quantity", xlab = "Total Quantity")

Insight :

Data Total Quantity paling banyak tersebar di rentang 20 hingga 80.
Terdapat data yang cukup jauh dari pola persebaran datanya (>80).

summary(data_train_ms)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.00   40.00   50.00   51.47   62.00  119.00

Insight :

Data Total Quantity paling rendah adalah 15 dan paling tinggi adalah 119.
Dengan Q3 adalah 62, maka dapat disimpulkan bahwa 75% dari data kurang dari 62, sehingga data 119 bisa dikatakan eksteme.

Maka akan dilalakukan penanganan untuk outlier dengan teknik Winsorize, yaitu memotong data ekstrim dan menggantinya dengan angka yang tidak ektrem. Dalam kasus ini, saya memotong 7.5% data total quantity dari data train dan menggantinya dengan angka pada data ke 7.5%, serta memotong 7.5% tertinggi dari data train dan menggantinya dengan angka pada data ke 92.5%.

data_train_ms_out <- Winsorize(x = data_train_ms, probs = c(0.075, 0.925))

stlm_ms_out <- stlm(y = data_train_ms_out, method = "arima", s.window = c(7, 7*4, 7*4*2, 7*4*3))

stlm_ms_out_f <- forecast(stlm_ms_out, h = 85)

> Evaluasi Model

eval <-
  rbind(
    accuracy(stlm_ms$fitted, data_train_ms),
    accuracy(stlm_ms_out$fitted, data_train_ms_out)
    )

rownames(eval) <- c("Multiple Seasonal ARIMA", "Multiple Seasonal ARIMA without Outlier")
round(eval,2)

##                                            ME RMSE  MAE   MPE  MAPE  ACF1
## Multiple Seasonal ARIMA                 -0.04 9.53 7.67 -4.02 16.98 -0.03
## Multiple Seasonal ARIMA without Outlier -0.04 8.22 6.59 -2.89 14.19 -0.05
##                                         Theil's U
## Multiple Seasonal ARIMA                      0.33
## Multiple Seasonal ARIMA without Outlier      0.36

Insight :

Dari hasil evaluasi di atas, dapat disimpulkan bahwa model ARIMA dengan Multiple Seasonal tanpa Outlier memberikan nilai error terkecil. Dimana nilai MAE-nya adalah 6.59, berarti rentang kesalahan data prediksi terhadap data aktual senilai (+/-) 6.59.

# Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_prod_ms, # data aktual full
              forecast.obj = stlm_ms_out_f, # hasil forecast
              test = data_test_ms) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA adalah ARIMA(3,0,5), dimana seasonal yang digunakan adalah pola mingguan, 1 bulanan, 2 bulan, dan 3 bulan secara bersamaan.
Terlihat bahwa hasil forecasting sudah bisa mengikuti fluktuasi dari data dan sudah mendekati data aktualnya.

> Evaluasi hasil forecasting

eval_test <-
  rbind(
    accuracy(stlm_ms_out_f$mean, data_test_ms)
    )

rownames(eval_test) <- c("Multiple Seasonal ARIMA without Outlier")
round(eval_test,2)

##                                            ME  RMSE   MAE    MPE  MAPE  ACF1
## Multiple Seasonal ARIMA without Outlier -5.36 17.94 14.34 -22.12 36.52 -0.09
##                                         Theil's U
## Multiple Seasonal ARIMA without Outlier      0.88

Insight :

Dari hasil evaluasi forecasting di atas diketahui nilai MAE-nya adalah 14.34, berarti rentang kesalahan data prediksi terhadap data aktual senilai (+/-) 14.34.

Assumption

Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena dengan menggunakan residual data, kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model. Metode forecasting yang baik menghasilkan nilai residual berikut ini:

Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
Residual memiliki rata-rata 0.

Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

> No-autocorrelation residual

H0: residual has no-autocorrelation H1: residual has autocorrelation yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series bisa menggunakan beberapa cara, yaitu: melakukan uji Ljung-box dengan menggunakan fungsi Box.test(residual model, type = “Ljung-Box)

# menggunakan Ljung-Box test
Box.test(stlm_ms_out_f$residuals, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  stlm_ms_out_f$residuals
## X-squared = 0.64771, df = 1, p-value = 0.4209

Insight :

Nilai p-value = 0.4209 > alpha = 0.05, artinya H0 diterima, residual tidak ada autokorelasi pada residual.

> Normality residual

H0: residual menyebar normal H1: 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).

shapiro.test(stlm_ms_out_f$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  stlm_ms_out_f$residuals
## W = 0.9963, p-value = 0.7578

Insight :

Nilai p-value = 0.7578 > alpha = 0.05, artinya H0 diterima, residual menyebar normal.

Kesimpulan

Dari hasil analisis dan pembahasan yang telah diuraikan sebelumnya, didapatkan beberapa kesimpulan sebagai berikut:

Forecasting menggunakan multiple seasonality menghasilkan prediksi yang lebih baik.
Model terbaik dalam memprediksi data total quantity harian adalah model ARIMA(3,0,5) dengan multiple seasonality. Dan dengan melakukan outlier treatment, yaitu memotong data yang ekstrim dengan metode Winsorize dapat meningkatkan akurasi dari prediksi.
Model ARIMA(3,0,5) dengan multiple seasonality tanpa Outlier juga memenuhi uji asumsi terhadap distribusi residual yang normal dan juga tidak adanya autokorelasi terhadap residual.

Rekomendasi

Model yang dibuat sudah mampu membaca pola seasonal dari data, namun tidak semua seasonal dapat ditangkap dengan baik. Maka untuk membuat model yang lebih baik bisa menggunakan metode yang lebih kompleks atau menggunakan deep learning.
Pola total quantity data terlihat inkonsisten hal ini menyebabkan seasonal sulit ditangkap, hal ini mungkin ada pengaruh dari faktor luar yang tidak dimasukan kedalam model.
Data total quantity berasal dari penjumlahan penjualan quantity semua produk yang ada yaitu Cashew,Cheese Stick,Choco Bar,Coffee Candy,Crackers,Ginger Candy,Oat,Potato Chip,Thai Tea, dan Yoghurt dengan pola penjualan yang mungkin berbeda. Sebaiknya bisa memprediksi untuk masing-masing produk agar forecasting quantity harian untuk tiap produk tepat dan pola penjualan tiap produk bisa teramati.

Additional: Forecasting Produk dengan penjualan total quantity terbanyak

1. Read Data

data_ts2 <- read.csv("Data/dataset-ts2.csv", sep = ";") %>%
  janitor::clean_names()

# Cek struktur data
glimpse(data_ts2)

## Rows: 2,666
## Columns: 3
## $ product_name   <chr> "Cashew", "Cashew", "Cashew", "Cashew", "Cashew", "Cash…
## $ date           <chr> "2022-01-02", "2022-01-03", "2022-01-04", "2022-01-06",…
## $ total_quantity <int> 11, 6, 8, 3, 1, 7, 4, 3, 3, 2, 7, 1, 6, 1, 1, 2, 5, 3, …

Cek produk yang memiliki total quantity terbanyak:

data_ts2 %>% 
  group_by(product_name) %>% 
  summarise(Quantity = sum(total_quantity)) %>% 
  arrange(desc(Quantity))

Didapat bahwa produk Thai Tea memiliki total quantity terbanyak. Maka sebagai contoh akan menggunakan akan menggunakan data produk Thai Tea.

# Only Use 1 Produk
data_tea <- data_ts2 %>% 
  filter(product_name ==  "Thai Tea") %>% 
  select(date, total_quantity)

2. Data Cleansing

> Mengubah tipe variabel

Dibagian ini variabel data akan disesuaikan tipe datanya.

data_tea <- data_tea %>% 
  # Ubah type data yang belum sesuai
  mutate(date = as.Date(date)
         )

glimpse(data_tea)

## Rows: 330
## Columns: 2
## $ date           <date> 2022-01-01, 2022-01-02, 2022-01-03, 2022-01-04, 2022-0…
## $ total_quantity <int> 3, 2, 9, 16, 7, 19, 3, 6, 13, 13, 4, 8, 20, 22, 5, 2, 8…

> Karakteristik data time series

Data harus terurut berdasarkan periode waktunya

data_tea <- data_tea %>%
  arrange(date)

Tidak boleh ada waktu atau periode yang terlewat/bolong

range(data_tea$date)

## [1] "2022-01-01" "2022-12-30"

nrow(data_tea)

## [1] 330

Insight :

Rentang data berawal dari 1 Januari 2022 hingga 30 Desember 2022, harusnya hingga 31 Desember 2022.
Jumlah baris pada data ada 330, artinya ada hari yang hilang/tidak ada transaksi.

Dari permasalahan diatas maka perlu ditambahkan data/rentang waktu yang hilang.

Des_31 <- c("2022-12-31", 0)
data_tea <- rbind(data_tea, Des_31)

data_tea <- data_tea %>% 
  pad()

## pad applied on the interval: day

range(data_tea$date)

## [1] "2022-01-01" "2022-12-31"

nrow(data_tea)

## [1] 365

Insight :

Rentang data sudah berawal dari 1 Januari 2022 hingga 31 Desember 2022.
Jumlah baris pada data ada 365, artinya tepat 1 Tahun.

Data tidak boleh ada yang missing

data_tea %>% 
  is.na() %>% 
  colSums()

##           date total_quantity 
##              0             34

Insight :

Terdapat 34 missing value pada kolom total_quantity.

Oleh karena itu maka akan dilakukan penanganan dengan cara inputasi nilai, nilai yang dimasukan adalah 0 (tidak ada transaksi).

data_tea <- data_tea %>% 
  mutate(total_quantity = as.numeric(ifelse(is.na(total_quantity), 0, total_quantity)))

data_tea %>% 
  is.na() %>% 
  colSums()

##           date total_quantity 
##              0              0

Insight :

Sudah tidak terdapat missing value.

3. Eksplorasi Data

Ekplorasi data dilakukan agar kita menjadi lebih paham dengan gambaran awal dari data yang akan digunakan.

> Trend data

data_tea %>% 
  ts_plot(title = "Trend quantity penjualan produk selama 2022",
          Xtitle = "Date", 
          Ytitle = "Total Quantity",
          line.mode = "lines+markers")

Insight :

Seperti yang terlihat dari hasil plot bahwa: 1. Total Quantity Produk Thai Tea sangat berfluktuatif, dan secara sekilas mayoritas bergerak diantara 5 hingga 20. 2. Terdapat data yang jauh dari pergerakan data lainnya, seperti 2 Maret total quantity mencapai 28. 3. Terdapat cukup banyak data bernilai 0, ini dikarenakan proses inputasi missing value.

> Dekomposisi

Tujuan dari melakukan dekomposisi adalah untuk melihat apakah dalam data terdapat trend dan seasonal.

data_tea %>% 
  .$total_quantity %>% 
  ts(frequency = 7) %>% 
  decompose() %>% 
  autoplot() +
  labs(title = "Dekomposisi Mingguan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Pada bagian ini data masih terlihat berfluktuasi sehingga trend nya belum terlihat jelas.
Seasonal: Terlihat bahwa data memiliki pola mingguan, karena terlihat pola berulang pada seasonal.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

data_tea %>% 
  .$total_quantity %>% 
  ts(frequency = 7*4) %>% 
  decompose() %>% 
  autoplot() +
  labs(title = "Dekomposisi Bulanan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Terlihat data memiliki pola menurun, namun masih terlihat berfluktuasi.
Seasonal: Terlihat bahwa data memiliki pola mingguan, karena terlihat pola berulang pada seasonal.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

data_tea %>% 
  .$total_quantity %>% 
  msts(seasonal.periods = c(7, 7*4)) %>% 
  mstl() %>% 
  autoplot() +
  labs(title = "Dekomposisi Mingguan + Bulanan") +
  theme(plot.title = element_text(hjust = 0.5))

Insight :

Trend: Terlihat data memiliki pola menurun dan fluktuasi menjadi lebih smooth.
Seasonal: Terlihat bahwa data memiliki pola bulanan, karena terlihat pola berulang pada seasonal28.
Error: Terlihat error masih memiliki pola tertentu, artinya masih ada pola yang belum tertangkap oleh model.

> Cross Validation Dataset

Dalam kasus ini data test diambil sebanyak 85 sampel berhubung keseluruhan data hanya berjumlah 365.

n_train <- 10*4*7
n_test <- nrow(data_tea) - n_train

data_tea_ts1 <- ts(data_tea[,2], start = 1, frequency = 1)
data_tea_ts7 <- ts(data_tea[,2], start = 1, frequency = 7)
data_tea_ts28 <- ts(data_tea[,2], start = 1, frequency = 7*4)

tea_train_s1 <- head(x = data_tea_ts1, n = n_train)
tea_test_s1 <- tail(x = data_tea_ts1, n = n_test)

tea_train_s7 <- head(x = data_tea_ts7, n = n_train)
tea_test_s7 <- tail(x = data_tea_ts7, n = n_test)

tea_train_s28 <- head(x = data_tea_ts28, n = n_train)
tea_test_s28 <- tail(x = data_tea_ts28, n = n_test)

data_tea_ms <- msts(data = data_tea[,2], seasonal.periods = c(7, 7*4, 7*4*2, 7*4*3))

tea_train_ms <- head(x = data_tea_ms, n = n_train)
tea_test_ms <- tail(x = data_tea_ms, n = n_test)

> Stasionarity

Stasionarity time series memiliki arti bahwa pada data time series yang kita miliki tidak memiliki trend maupun seasonal dan memiliki variansi konstan. Melakukan pengujian menggunakan adf.test()

# Augmented Dickey-Fuller Test
tea_train_s1 %>% 
  adf.test()

## Warning in adf.test(.): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -6.5343, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

H0: data tidak stasioner H1: data stasioner p-value < alpha (0.05), artinya data stasioner.

Didapatkan bahwa data yang dimiliki sudah stasioner, maka tidak perlu melakukan differencing.

Modelling: Forecasting menggunakan Metode ARIMA

Sama seperti sebelumnya, dalam pembangunan model akan digunakan beberapa metode, diantaranya ARIMA dan Seasonal ARIMA (SARIMA). Model tersebut dipilih karena pada data yang digunakan terdapat trend dan seasonal.

tea_auto_arima <- auto.arima(y = tea_train_s1)
tea_auto_sarima7 <- stlm(y = tea_train_s7, method = "arima", s.window = 7)
tea_auto_sarima28 <- stlm(y = tea_train_s28, method = "arima", s.window = 28)
tea_stlm_ms <- stlm(y = tea_train_ms, method = "arima", allow.multiplicative.trend = TRUE, s.window = c(7, 7*4, 7*4*2, 7*4*3))

> Evaluasi Model

eval2 <-
  rbind(
    accuracy(tea_auto_arima$fitted, tea_train_s1),
    accuracy(tea_auto_sarima7$fitted, tea_train_s7),
    accuracy(tea_auto_sarima28$fitted, tea_train_s28),
    accuracy(tea_stlm_ms$fitted, tea_train_ms)
    )

rownames(eval2) <- c("ARIMA","SARIMA(Seasonal 7)","SARIMA(Seasonal 28)", "Multiple Seasonal ARIMA")
round(eval2,2)

##                           ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## ARIMA                   0.00 5.37 4.18 -Inf  Inf  0.00         0
## SARIMA(Seasonal 7)      0.01 4.40 3.49 -Inf  Inf  0.00         0
## SARIMA(Seasonal 28)     0.00 5.05 4.00 -Inf  Inf  0.00         0
## Multiple Seasonal ARIMA 0.00 2.80 2.26  NaN  Inf -0.01         0

Insight :

Dari hasil evaluasi di atas, dapat disimpulkan bahwa model telah dapat memprediksi total quantity produk Thai Tea. Adapun model terbaik adalah metode ARIMA dengan Multiple Seasonal dengan nilai error terkecil. Dimana nilai MAE-nya adalah 2.26, berarti rentang kesalahan data prediksi terhadap data aktual senilai (+/-) 2.26.

tea_stlm_ms_f <- forecast(tea_stlm_ms, h = 85)

#Visualisasi perbandingan prediksi dan aktual
test_forecast(actual = data_tea_ms, # data aktual full
              forecast.obj = tea_stlm_ms_f, # hasil forecast
              test = tea_test_ms) # data test

Insight :

Model terbaik yang terbentuk dengan metode ARIMA adalah ARIMA(4,0,2), dimana seasonal yang digunakan adalah pola mingguan, 1 bulanan, 2 bulan, dan 3 bulan secara bersamaan.
Terlihat bahwa hasil forecasting sudah bisa mengikuti fluktuasi dari data dan sudah mendekati data aktualnya.

> Evaluasi Hasil Forecasting

eval_test2 <-
  rbind(
    accuracy(tea_stlm_ms_f$mean, tea_test_ms)
    )

rownames(eval_test2) <- c("Multiple Seasonal ARIMA")
round(eval_test2,2)

##                            ME RMSE  MAE  MPE MAPE  ACF1 Theil's U
## Multiple Seasonal ARIMA -1.13 6.91 5.73 -Inf  Inf -0.12         0

Insight :

Dari hasil evaluasi forecasting di atas diketahui nilai MAE-nya adalah 5.73, berarti rentang kesalahan data prediksi terhadap data aktual senilai (+/-) 5.73.

Asumsi

Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
Residual memiliki rata-rata 0.

Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

> No-autocorrelation residual

H0: residual has no-autocorrelation H1: residual has autocorrelation yang diinginkan p-value > 0.05 (alpha), no-autocorrelation

# menggunakan Ljung-Box test
Box.test(tea_stlm_ms_f$residuals, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  tea_stlm_ms_f$residuals
## X-squared = 0.050082, df = 1, p-value = 0.8229

Insight :

Nilai p-value = 0.1192 > alpha = 0.05, artinya H0 diterima, residual tidak ada autokorelasi pada residual.

> Normality residual

H0: residual menyebar normal H1: 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).

shapiro.test(tea_stlm_ms_f$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  tea_stlm_ms_f$residuals
## W = 0.99791, p-value = 0.9793

Insight :

Nilai p-value = 0.9793 > alpha = 0.05, artinya H0 diterima, residual menyebar normal.

Kesimpulan

Dari hasil analisis dan pembahasan yang telah diuraikan sebelumnya, didapatkan beberapa kesimpulan sebagai berikut:

Forecasting menggunakan multiple seasonality menghasilkan prediksi yang lebih baik.
Model terbaik dalam memprediksi data total quantity harian produk Thai Tea adalah model ARIMA(4,0,2) dengan multiple seasonality.
Model ARIMA(4,0,2) dengan multiple seasonality juga memenuhi uji asumsi terhadap distribusi residual yang normal dan juga tidak adanya autokorelasi terhadap residual.

Prediksi Total Quantity harian dari Produk Terjual

Machine Learning Regression (Time Series)

Anugrah Hidayat

July 15, 2023

Pendahuluan

> Maksud dan Tujuan

> Benefit

> Batasan

Data Wrangling

1. Import Library

2. Read Data

3. Data Cleansing

> Mengubah tipe variabel

> Karakteristik data time series

Eksplorasi Data

> Trend data

> Dekomposisi

Modelling: Forecasting menggunakan Metode ARIMA

> Cross Validation Dataset

> Stasionarity

> Differencing

1. ARIMA (Autoregressive Integrated Moving Average)

2. SARIMA (Seasonal ARIMA)

3. Multiple Seasonal

> Evaluasi Model

Improving Model

> Evaluasi Model

> Evaluasi hasil forecasting

Assumption

> No-autocorrelation residual

> Normality residual

Kesimpulan

Rekomendasi

Additional: Forecasting Produk dengan penjualan total quantity terbanyak

1. Read Data

2. Data Cleansing

> Mengubah tipe variabel

> Karakteristik data time series

3. Eksplorasi Data

> Trend data

> Dekomposisi

> Cross Validation Dataset

> Stasionarity

Modelling: Forecasting menggunakan Metode ARIMA

> Evaluasi Model

> Evaluasi Hasil Forecasting

Asumsi

> No-autocorrelation residual

> Normality residual

Kesimpulan