Tentang

Report ini merupakan penyempurnaan Dive Deeper Time Series yang dapat dilihat di sini. Objek time series saya persempit ke >= tahun 1951 agar tidak mempertimbangkan downtrend di beberapa tahun sebelumnya.

Memuat Library R

library(tidyverse)
library(forecast)
library(MLmetrics)
library(tseries)
library(lubridate)
theme_set(theme_minimal())

Read Data

birth <- read.csv("data_input/nybirth.csv")
head(birth,100)
range(birth$date)
#> [1] "1946-01-01" "1959-12-01"

Membuat Object TS

birth_ts1 <- ts(data = birth$births,
                  start = c(1949,1),
                  frequency = 12)

Visualisasi Object TS

birth_ts1 %>% 
  autoplot()

Membuang data <= 1949-12-01

# library(lubridate)
new_dat <- birth %>% 
  mutate(date = ymd(date)) %>% 
  filter(date > "1950-12-01")
birth_ts2 <- ts(data = new_dat$births,
                  start = c(1951,1),
                  frequency = 12)
birth_ts2 %>% 
  autoplot()

class(birth_ts2)
#> [1] "ts"
birth_ts2 %>% 
  decompose() %>%  
  autoplot()

Splitting train dan test

Splitting data dua tahun terakhir sebagai data test, dan sisanya data train

birth_test <- tail(birth_ts2, 24) #hanya mengambil 24 bulan terakhir

birth_train <- head(birth_ts2, -24) #tidak mengambil 24 bulan terakhir
birth_train %>% 
  decompose() %>% 
  autoplot()

Model Fitting

Dengan mengggunakan metode smoothing auto
Dengan nilai smoothing mendefinisikan nilai alpha beta dan gamma (manual)

birth_hw <- HoltWinters(x = birth_train)

birth_hw2 <- HoltWinters(x = birth_train, alpha = 0.8, beta = 0 , gamma = 0.98)
#gamma = seasonal
#alpha = error (residual)
#beta = trend
birth_forecast_auto <- forecast(birth_hw, h = 24)
birth_forecast_manual <- forecast(birth_hw2, h = 24)

Visualisasi data dan model

birth_ts2 %>% 
  autoplot()+
  autolayer(birth_hw$fitted[,1], series =  "Auto", lwd = 1)+
  autolayer(birth_hw2$fitted[,1], series =  "Alpha 0.8; beta 0; gamma 0.98", lwd = 1)

birth_ts2 %>% 
  autoplot()+
  autolayer(birth_forecast_auto)

  # autolayer(birth_forecast_manual, col = "red")
birth_ts2 %>% 
  autoplot()+
  autolayer(birth_forecast_manual, col = "red")

Model Evaluation

MAPE(y_pred =  birth_forecast_manual$mean, y_true = birth_test) * 100 #$ mean adalah poin forecastnya  
#> [1] 3.100369
MAPE(y_pred =  birth_forecast_auto$mean, y_true = birth_test) * 100 #dikali 100(persen) supaya keluar nya bukan satu komaan
#> [1] 3.079305

nilai MAPE antara model auto dan manual terpaut 0.1 persen, rata-rata prosentase error dari model manual lebih besar

Seasonal ARIMA (SARIMA)

Memisahkan data untuk keperluan differencing

birth_train_a <- birth_train 
birth_test_a <- birth_test

Memeriksa stasionaritas data

birth_train_a%>% 
  adf.test() #P-value 0.01, data sebetulnya sudah stasioner 
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -5.0455, Lag order = 4, p-value = 0.01
#> alternative hypothesis: stationary

Data secara statistik sudah stasioner (p-value < 0.05),
namun setelah dilakukan pemeriksaan plot masih perlu dilakukan differencing

Seasonality Analysis

birth_train_a %>% 
  autoplot()

differed <- birth_train_a %>% 
  diff(lag = 12) %>% 
  diff(lag = 1)
library(TSstudio)

ts_seasonal(birth_train_a, type = "all")
ts_heatmap(birth_train_a)
ts_seasonal(differed, type = "all")
ts_heatmap(differed)

Berdasarkan hasil pemeriksaan seasonality plot dan heatmap, dapat kita lihat bahwa model memiliki autokorelasi yang lebih kecil dengan model sebelum differencing.

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

birth_train_a %>% 
  diff(lag = 12) %>% 
  diff(lag = 1 ) %>%
  tsdisplay(lag.max = 36)

Seasonal AR

*P : 1

*D : 1

*Q : 1

Keseluruhan Data

*p : 2, 3, 5

*d : 1

*q : 2

Kombinasi model yang mungkin terbentuk adalah:

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

  • ARIMA(2,1,2)(1,1,1) [12]

  • ARIMA(3,1,2)(1,1,1) [12]

  • ARIMA(5,1,2)(1,1,1) [12]

birth_arima212_111 <- arima(x= birth_train_a, order = c(2,1,2), seasonal = c(1,1,1))
birth_arima_312_111 <- arima(x= birth_train_a, order = c(3,1,2), seasonal = c(1,1,1))
birth_arima_512_111 <- arima(x= birth_train_a, order = c(5,1,2), seasonal = c(1,1,1))
birth_auto_arima <- auto.arima(birth_train_a)
birth_auto_arima
#> Series: birth_train_a 
#> ARIMA(2,0,0)(2,1,0)[12] with drift 
#> 
#> Coefficients:
#>          ar1      ar2     sar1     sar2   drift
#>       0.5143  -0.2139  -0.9442  -0.6365  0.0495
#> s.e.  0.1261   0.1221   0.1155   0.0978  0.0035
#> 
#> sigma^2 estimated as 0.3065:  log likelihood=-65.79
#> AIC=143.58   AICc=144.88   BIC=157.24

Memeriksa besaran informasi yang hilang (AIC)

# birth_arima212_111$aic
# birth_arima313_211$aic
birth_arima212_111$aic
#> [1] 154.3821
birth_arima_312_111$aic
#> [1] 156.1445
birth_arima_512_111$aic
#> [1] 158.0756
birth_auto_arima$aic
#> [1] 143.5841

Dari hasil penghitungan AIC dapat dilihat bahwa model SARIMA (3,1,3)(2,1,1)[12] memiliki lebih sedikit nilai AIC

Memeriksa besaran prosentase error

accuracy(birth_arima212_111)
#>                      ME      RMSE       MAE       MPE     MAPE      MASE
#> Training set 0.04696387 0.4870454 0.3775986 0.1403305 1.445971 0.3158515
#>                     ACF1
#> Training set -0.02000742
accuracy(birth_arima_312_111)
#>                      ME      RMSE       MAE      MPE     MAPE      MASE
#> Training set 0.05015363 0.4848575 0.3735564 0.153109 1.431428 0.3124703
#>                     ACF1
#> Training set -0.00885539
accuracy(birth_arima_512_111)
#>                      ME      RMSE       MAE       MPE    MAPE      MASE
#> Training set 0.04370168 0.4807244 0.3680617 0.1295703 1.41316 0.3078741
#>                        ACF1
#> Training set -0.00001872018
accuracy(birth_auto_arima)
#>                       ME      RMSE       MAE        MPE     MAPE      MASE
#> Training set -0.03226246 0.4944081 0.3789365 -0.1819979 1.457433 0.4482171
#>                    ACF1
#> Training set 0.01087832

Dari hasil penghitungan accuracy dapat dilihat bahwa model SARIMA (3,1,3)(2,1,1)[12] memiliki nilai MAPE yang lebih kecil

Model Forecasting

birth_a_forecast <- forecast(birth_arima_512_111, h = 24)
birth_au_forecast <- forecast(birth_auto_arima, h = 24)
birth_ts2 %>% 
  autoplot()+
  autolayer(birth_a_forecast, lwd = 1, series = "ARIMA (5,1,2)(1,1,1)", col = "cyan", alpha = 0.5)

  # autolayer(birth_au_forecast, lwd = 1, series = "ARIMA auto", col = "orange", alpha = 0.5)
  # autolayer(birth_forecast_manual, lwd = 1, series = "HoltWinters Manual", col = "red", alpha = 0.5)
birth_ts2 %>% 
  autoplot()+
  autolayer(birth_au_forecast, lwd = 1, series = "ARIMA auto", col = "green", alpha = 0.35)

birth_ts2 %>% 
  autoplot()+
  autolayer(birth_forecast_manual, lwd = 1, series = "Manual Holtwinters", col = "navyblue", alpha = 0.5)

Evaluasi model dengan data test

# seasonal ARIMA manual (5,1,2)(1,1,1)
MAPE(y_pred = birth_a_forecast$mean, birth_test_a)*100
#> [1] 3.463261
# auto ARIMA
MAPE(y_pred = birth_au_forecast$mean, birth_test_a)*100 
#> [1] 4.619402
# manual HoltWinters

MAPE(y_pred = birth_forecast_manual$mean, birth_test_a)*100
#> [1] 3.100369

Kesimpulan

Untuk workflow yang saya buat saat ini, model-model dengan tuning manual lebih unggul dibandingkan dengan yang auto generated. Namun memang, antara model Seasonal Arima (5,1,2)(1,1,1) dan Model Holtwinters (alpha = 0.8; beta = 0.0; gamma = 0.98) masih terpaut selisih 1 persen, dengan model HoltWinters manual lebih unggul (MAPE = 2.4 persen). Hal ini mungkin saja disebabkan karena pemotongan objek time series dari tahun 1951. Analisis saya, sepertinya pada tahun 1951 pun masih mengandung downtrend dari tahun-tahun sebelumnya. Apabila saya persempit lagi dengan mempertimbangkan uptrend yang lebih bersih (yang dapat dilihat konstan sejak 1952, kemungkinan semua model bisa menghasilkan error yang lebih kecil).