1. Read Dataset

dautonsa <- read.csv("DAUTONSA.csv")
dautonsa %>% head()

glimpse(dautonsa)

#> Rows: 630
#> Columns: 2
#> $ DATE     <chr> "1967-01-01", "1967-02-01", "1967-03-01", "1967-04-01", "1967…
#> $ DAUTONSA <dbl> 564.1, 509.1, 670.4, 710.2, 744.8, 780.2, 627.2, 517.2, 547.3…

2. Subset Dataset

dautonsa <- dautonsa %>% 
  mutate(DATE = as.Date(DATE)) %>% 
  filter(DATE >= "2015-01-01")

3. Range Dataset

dautonsa$DATE %>% range()

#> [1] "2015-01-01" "2019-06-01"

4. Missing Value

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

#>     DATE DAUTONSA 
#>        0        0

5. Duplicate Data

anyDuplicated(dautonsa)

#> [1] 0

6. Checking Dataset Structure

glimpse(dautonsa)

#> Rows: 54
#> Columns: 2
#> $ DATE     <date> 2015-01-01, 2015-02-01, 2015-03-01, 2015-04-01, 2015-05-01, …
#> $ DAUTONSA <dbl> 386.010, 419.495, 525.027, 477.314, 563.176, 491.397, 481.665…

1. Creating Time Series Object

ts_dautonsa <- ts(data = dautonsa$DAUTONSA,
                  start = c(2015,1),
                  end = c(2019,6),
                  frequency = 12) 

ts_dautonsa %>% autoplot()

2. Decomposition

Decomposition adalah suatu tahapan dalam time series analisis yang digunakan untuk menguraikan beberapa komponen dalam time series data.

💡 Komponen dalam time series :

• Trend : pola data secara general, cenderung untuk naik atau turun. Jika ada trend masih terdapat pola artinya masih ada pola yang belum terurai dengan baik.
• Seasonal : pola musiman yang membentuk pola berulang pada periode waktu yang tetap
• Error/Reminder/Random : pola yang tidak dapat ditangkap dalam trend dan seasonal Sebelum melakukan modeling forecasting kita perlu mengamati objek timeseries dari hasil decompose. Ide utama dari decompose adalah untuk menguraikan ketiga komponen dari objek ts (trend, seasonal, residual).

decom_dautonse <- decompose(ts_dautonsa)

decom_dautonse %>% autoplot() %>% ggplotly() 

Insight :

• Trend data kita cenderung menurun.
• Data kita memiliki pola seasonal.
• Pola dataset kita dapat dikategorikan sebagai additive seasonal karena kenaikan yang relatif konstan, tidak multiplikasi.

3. Cross Validation

# test menggunakan `tail()`
test_dautonse <-  tail(ts_dautonsa, 18)

test_dautonse %>% autoplot()

Data test kita terdiri dari Januari 2018 hingga Juli 2019 (1,5 tahun atau 18 bulan).

# train menggunakan `head()`
train_dautonse <- head(ts_dautonsa, -length(test_dautonse))
train_dautonse %>% autoplot()

Data train kita terdiri dari Januari 2015 hingga Desember 2018 (3 tahun atau 36 bulan).

1. Model Building

model_sma <- SMA(x = train_dautonse, n = 5)

2. Forecasting

forecast_sma <- forecast(object = model_sma, h = 18)
forecast_sma

#>          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
#> Jan 2018       343.7615 338.6416 348.8815 335.9312 351.5919
#> Feb 2018       337.2671 332.2438 342.2905 329.5846 344.9497
#> Mar 2018       347.6975 342.5187 352.8762 339.7772 355.6177
#> Apr 2018       354.9979 349.7102 360.2855 346.9111 363.0846
#> May 2018       351.8024 346.5623 357.0425 343.7883 359.8165
#> Jun 2018       367.0229 361.5559 372.4899 358.6619 375.3839
#> Jul 2018       369.1836 363.6843 374.6828 360.7732 377.5940
#> Aug 2018       359.8292 354.4691 365.1892 351.6317 368.0267
#> Sep 2018       356.5469 351.2356 361.8582 348.4239 364.6698
#> Oct 2018       340.6335 335.5591 345.7079 332.8729 348.3941
#> Nov 2018       328.1200 323.2319 333.0081 320.6443 335.5957
#> Dec 2018       326.2385 321.3783 331.0987 318.8055 333.6716
#> Jan 2019       303.8766 298.7670 308.9862 296.0621 311.6911
#> Feb 2019       297.7536 292.7468 302.7605 290.0964 305.4109
#> Mar 2019       306.5603 301.4053 311.7154 298.6764 314.4443
#> Apr 2019       312.5788 307.3224 317.8352 304.5398 320.6178
#> May 2019       309.3424 304.1402 314.5445 301.3863 317.2984
#> Jun 2019       322.2758 316.8559 327.6957 313.9868 330.5648

• Point Forecast: nilai forecast untuk periode yang ingin kita forecast
• Lo 80 & Hi 80 : rentang tebakan dari hasil forecast untuk confidence level 80%
• Lo 95 & Hi 95 : rentang tebakan dari hasil forecast untuk confidence level 95%

3. Visualization

train_dautonse %>% 
  autoplot() +
  autolayer(test_dautonse, series = "Test") + # data test
    autolayer(forecast_sma$mean, series = "Forecast Data Test") + # hasil forecast model
    autolayer(model_sma, series = "Fitted Values") # fitted value data train

4. Evaluation

# model additive
accuracy(model_sma, train_dautonse)

#>                 ME     RMSE      MAE       MPE     MAPE       ACF1 Theil's U
#> Test set -7.131075 38.49832 30.63705 -2.526263 7.435305 0.02793552 0.6443543

# evaluasi data test
accuracy(forecast_sma$mean, test_dautonse)

#>                 ME     RMSE      MAE     MPE     MAPE       ACF1 Theil's U
#> Test set -2.429805 35.46845 29.25938 -1.6535 8.718635 -0.1693064 0.6429105

Hasil forecast model Simple Moving Average dengan ordo 5 belum dapat mengikuti pola data aktualnya. Cenderung hasil prediksi masih berada ditengah kumpulam data (belum baik) meskipun selisih MAPE antara data train dan test relatif rendah dan mendekati.

1. Model Building

model_holt <- HoltWinters(x = train_dautonse)
model_holt

#> Holt-Winters exponential smoothing with trend and additive seasonal component.
#> 
#> Call:
#> HoltWinters(x = train_dautonse)
#> 
#> Smoothing parameters:
#>  alpha: 0.2283017
#>  beta : 0.007234357
#>  gamma: 1
#> 
#> Coefficients:
#>            [,1]
#> a    363.171163
#> b     -3.785401
#> s1  -102.778660
#> s2   -35.887200
#> s3    47.149162
#> s4    12.797192
#> s5    37.324159
#> s6     9.540481
#> s7    -5.618603
#> s8    15.832135
#> s9    24.632602
#> s10  -33.029591
#> s11  -35.162948
#> s12   -1.112163

2. Forecasting

forecast_holt <- forecast(object = model_holt, h = 18)
forecast_holt

#>          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
#> Jan 2018       256.6071 233.2549 279.9593 220.8929 292.3213
#> Feb 2018       319.7132 295.7515 343.6749 283.0669 356.3594
#> Mar 2018       398.9641 374.3996 423.5286 361.3960 436.5323
#> Apr 2018       360.8268 335.6656 385.9879 322.3460 399.3075
#> May 2018       381.5683 355.8161 407.3205 342.1837 420.9530
#> Jun 2018       349.9992 323.6612 376.3373 309.7186 390.2798
#> Jul 2018       331.0548 304.1356 357.9739 289.8855 372.2240
#> Aug 2018       348.7201 321.2244 376.2158 306.6691 390.7711
#> Sep 2018       353.7352 325.6670 381.8033 310.8086 396.6617
#> Oct 2018       292.2876 263.6508 320.9244 248.4913 336.0838
#> Nov 2018       286.3688 257.1669 315.5707 241.7083 331.0293
#> Dec 2018       316.6342 286.8705 346.3979 271.1145 362.1539
#> Jan 2019       211.1823 173.0636 249.3010 152.8848 269.4798
#> Feb 2019       274.2884 235.7260 312.8507 215.3123 333.2644
#> Mar 2019       353.5393 314.5325 392.5461 293.8836 413.1950
#> Apr 2019       315.4019 275.9500 354.8539 255.0654 375.7385
#> May 2019       336.1435 296.2456 376.0414 275.1250 397.1620
#> Jun 2019       304.5744 264.2299 344.9190 242.8727 366.2761

• Point Forecast: nilai forecast untuk periode yang ingin kita forecast
• Lo 80 & Hi 80 : rentang tebakan dari hasil forecast untuk confidence level 80%
• Lo 95 & Hi 95 : rentang tebakan dari hasil forecast untuk confidence level 95%

3. Visualization

train_dautonse %>% 
  autoplot() +
  autolayer(test_dautonse, series = "test") + # data test
    autolayer(forecast_holt$mean, series = "forecast data test") + # hasil forecast model
    autolayer(model_holt$fitted[,1], series = "fitted values")  # fitted value data train

4. Evaluation

# model additive
accuracy(model_holt$fitted[,1], train_dautonse)

#>                 ME     RMSE      MAE        MPE     MAPE        ACF1 Theil's U
#> Test set -3.580492 18.19397 14.62881 -0.8975071 3.658825 -0.07230681 0.3786306

# evaluasi data test
accuracy(forecast_holt$mean, test_dautonse)

#>                ME     RMSE      MAE      MPE     MAPE       ACF1 Theil's U
#> Test set 11.11905 25.78153 22.05057 3.465379 6.868374 0.06610927 0.4726422

Hasil forecast model Holt’s Exponential Smoothing sudah dapat mengikuti pola data aktualnya. Cenderung hasil prediksi sudah mengikuti pola data meskipun masih terdapat beberapa titik yang belum mendekati data aktual. Selain itu, ukuran MAPE yang dihasilkan sudah baik, train (3.65%) dan test (5.88%) serta selisih keduanya sekitar 2% alias tidak underfitting dan overfitting.

1. Stationarity Check

adf.test(train_dautonse)

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  train_dautonse
#> Dickey-Fuller = -3.42, Lag order = 3, p-value = 0.07034
#> alternative hypothesis: stationary

karena p-value < alpha, dimana 0.040 < 0.05 artinya data sudah stasioner.

2. Model Building

membuat model SARIMA menggunakan auto.arima() SARIMA(p,d,q)(P,D,Q)[seasonal/freq]

model_auto_sarima <- auto.arima(train_dautonse, seasonal = T)
model_auto_sarima

#> Series: train_dautonse 
#> ARIMA(0,0,0)(0,1,0)[12] with drift 
#> 
#> Coefficients:
#>         drift
#>       -3.4797
#> s.e.   0.3872
#> 
#> sigma^2 = 540.9:  log likelihood = -109.06
#> AIC=222.12   AICc=222.69   BIC=224.48

model_auto_sarima$aic

#> [1] 222.1198

3. Forecasting

# forecast, h = 18
length(test_dautonse)

#> [1] 18

forecast_sarima <- forecast(model_auto_sarima, h=18)

4. Visualization

# menggunakan autoplot dan autolayer
train_dautonse %>% 
  autoplot()+
  autolayer(test_dautonse)+
  autolayer(forecast_sarima$mean)

5. Evaluation

menggunakan accuracy() dari package forecast dan melihat nilai MAPEnya

# accuracy
accuracy(model_auto_sarima$fitted, train_dautonse) # error train

#>                 ME     RMSE      MAE          MPE     MAPE      ACF1 Theil's U
#> Test set 0.1629594 18.58972 12.34508 -0.008536156 3.052588 0.0431331 0.3643444

accuracy(forecast_sarima$mean, test_dautonse) # error test

#>                ME     RMSE      MAE       MPE     MAPE       ACF1 Theil's U
#> Test set 1.853778 22.78772 17.60102 0.6124365 5.520045 0.04267141  0.409262

Hasil forecast model Seasonal ARIMA sudah dapat mengikuti pola data aktualnya. Cenderung hasil prediksi sudah mengikuti pola data meskipun masih terdapat beberapa titik yang belum mendekati data aktual. Selain itu, ukuran MAPE yang dihasilkan sudah lebih baik, train (3.05%) dan test (5.52%) serta selisih keduanya sekitar 2% alias tidak underfitting dan overfitting.

1. Model Building

Dengan menggunakan informasi data train sales_train dan test sales_test, kita akan coba untuk membuat model dengan menggunakan metode stlm() dengan mengatur parameter method = “ets”.

Parameter stlm() yang harus di atur:

• y : object time series
• s.window : pola seasonal yang ingin ditangkap
• method : metode forecast yang akan digunakan, tersedia ets dan arima
• modeling

model_stlm <- stlm(y = train_dautonse, s.window = 12, method = "ets")

2. Forecasting

Melakukan forecast untuk 18 bulan mendatang

forecast_stlm <- forecast(model_stlm, h = 18)

3. Visualization

ts_dautonsa %>% autoplot(series = "actual") +
  autolayer(test_dautonse, series = "test data") +
  autolayer(forecast_stlm$mean, series = "forecast") +
  scale_color_manual(values = c("black", "blue", "firebrick"))

4. Evaluation

accuracy(model_stlm$fitted, train_dautonse) # error train

#>                 ME    RMSE      MAE        MPE     MAPE        ACF1 Theil's U
#> Test set -0.854042 15.3851 11.42279 -0.2259621 2.695135 0.003772181 0.2841413

accuracy(forecast_stlm$mean, test_dautonse)

#>                 ME    RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -16.43207 26.5702 21.59139 -4.906387 6.676358 0.2433919 0.4982757

Hasil forecast model STLM (Seasonal Trend with Loess Model) sudah dapat mengikuti pola data aktualnya. Cenderung hasil prediksi sudah mengikuti pola data meskipun masih terdapat beberapa titik yang belum mendekati data aktual. Selain itu, ukuran MAPE yang dihasilkan sudah lebih baik, train (2.69%) dan test (6.67%) serta selisih keduanya sekitar 4% alias tidak underfitting dan overfitting.

1. No-autocorrelation Residual

Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series bisa menggunakan :

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

H0: residual has no-autocorrelation H1: residual has autocorrelation

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

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

#> 
#>  Box-Ljung test
#> 
#> data:  model_auto_sarima$residuals
#> X-squared = 0.072718, df = 1, p-value = 0.7874

p-value > 0.05 (alpha) artinya gagal tolak H0, artinya residual/error pada data tidak terdapat autocorrelation.

2. Normality of 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(model_auto_sarima$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_auto_sarima$residuals
#> W = 0.93241, p-value = 0.02964

hist(model_auto_sarima$residuals, xlab = "Residuals", main = "Residual Model SARIMA")

p-value < alpha artinya tolak H0, artinya residual/error pada data tidak menyebar normal.