0.1 R Markdown

library(dplyr) # data wrangling
library(lubridate) # date manipulation
library(forecast) # time series library
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
library(tseries) # adf.test
library(fpp) # usconsumtion
library(TSstudio) # mempercantik visualisasi timeseries
library(ggplot2)
library(tidyr)

1 Apa Itu Metode Peramalan Deret Waktu?

Metode Peramalan Deret Waktu atau Time Series adalah satu metode dalam statistika yang dapat melakukan pendugaan nilai suatu peubah di masa depan berdasarkan data-data historis sebagai dasar pemodelannya. Berbeda dengan analisis regresi, metode ini membutuhkan peubah Y dalam bentuk rentetan waktu secara terurut dalam suatu nilai numerik dengan peubah Y dipengaruhi oleh nilai Y itu sendiri pada masa lampau (\(Y_{t-1}\)). Adapun Proses untuk menduga nilai yang akan terjadi di masa depan disebut forecasting.

Analisa ini dapat dijumpai dalam melakukan peramalan curah hujan di suatu daerah dalam satu tahun, harga daging kambing bulanan, jumlah penumpang pesawat per semester dan berbagai kasus lainnya.

2 Persiapan Data

Data yang akan digunakan pada metode peramalan deret waktu ini adalah data iklim di kota Delhi, India yang dikumpulkan mulai 1 Januari 2013 hingga 24 April 2017.Data tersebut bisa diakses pada https://www.kaggle.com/sumanthvrao/daily-climate-time-series-data.

2.1 Membaca Data

Pertama, pembacaan data yang akan digunakan sebagai pemodelan deret waktu.

iklim <- read.csv("DailyDelhiClimateTrain.csv")


glimpse(iklim)
#> Rows: 1,462
#> Columns: 5
#> $ date         <chr> "2013-01-01", "2013-01-02", "2013-01-03", "2013-01-04"...
#> $ meantemp     <dbl> 10.000000, 7.400000, 7.166667, 8.666667, 6.000000, 7.0...
#> $ humidity     <dbl> 84.50000, 92.00000, 87.00000, 71.33333, 86.83333, 82.8...
#> $ wind_speed   <dbl> 0.0000000, 2.9800000, 4.6333333, 1.2333333, 3.7000000,...
#> $ meanpressure <dbl> 1015.667, 1017.800, 1018.667, 1017.167, 1016.500, 1018...

2.2 Deskripsi Data

  • Data tersebut meliputi
    • date
      • tanggal data tersebut direkam
    • meantemp
      • temperatur rata-rata pada hari tersebut
    • humidity
      • kelembapan udara
    • wind_speed
      • Kecepatan Angin
    • meanpressure
      • tekanan rata-rata

3 Pra-proses Data

Terdapat 4 peubah yang dapat digunakan untuk melakukan peramalan deret waktu pada set data tersebut. Peubah humidity selanjutnya dipilih untuk dilakukan peramalan kelembapan udara yang akan terjadi di masa selanjutnya. Adapun peubah date yang bertipe chr akan diubah menjadi peubah dengan tipe date untuk diketahui rentang waktu pada sub-set data tersebut.

iklim <- iklim %>% 
  mutate(date = ymd(date))


range(iklim$date)
#> [1] "2013-01-01" "2017-01-01"

Rentang tahun yang digunakan pada data train tersebut terdiri dari 4 tahun dari 1 Januari 2013 hingga 1 Januari 2017.

Adapun prasyarat yang harus dilakukan pada peramalan data deret waktu terdiri dari:

  • Data tidak boleh ada yang NA
  • Data harus terurut berdasarkan waktu
  • Tidak diperkenankannya data dengan tahun yang hilang
colSums(is.na(iklim))
#>         date     meantemp     humidity   wind_speed meanpressure 
#>            0            0            0            0            0

Pengecekan tersebut menunjukan data tak memiliki NA, data pada sub-set data ini sudah terurut dalam bentuk harian per bulan dan per tahun, dan tidak terdapat data dengan tahun yang hilang.

4 Membangun Model

Data yang akan digunakan pada analisa harus diubah menjadi data dengan tipe time series. Proses konversi data tersebut menjadi data bertipe time series adalah sebagai berikut:

iklim_ts <- ts(data = iklim$humidity, start = c(2013,1,1), frequency = 365)

Karena data yang terekam berupa data harian, maka frequency yang digunakan adalah 365 untuk kita lihat polanya dalam bentuk tahunan.

iklim_ts %>% ts_plot()

##Validasi Silang

Data akan dibagi menjadi 2 yakni data train yang akan digunakan untuk pembentukan model dan data test yang akan digunakan sebagai pengujian akurasi terhadap model data trainnya.

Model yang akan dibentuk diambil dari 3 tahun pertama dan 1 tahun selanjutnya sebagi data test.

iklim_train <- head(iklim_ts, 3*365)
iklim_test <- tail(iklim_ts, length(iklim_ts)-length(iklim_train))

Selanjutnya akan dilakukan pengecekan menggunakan decompose untuk mencari tahu faktor additive atau multiplikatif yang terdapat dalam data. Pengecekan dilakukan terhadap data train

iklim_train%>% decompose() %>% autoplot()

Data tersebut memiliki pola additive dan pola musiman. Metode peramalan deret waktu yang akan digunakan berdasarkan pola data tersebut adalah Holt’s Winters Exponential (Triple Exponential Smoothing) dan Seasonal Arima. Selain dengan kedua pemodelan tersebut, akan pula diuji menggunakan pemodelan STLM (Seasonal Trend with Loess Model).

4.1 Pemodelan dengan Holt’s Winters Exponential (Triple Exponential Smoothing)

merupakan metode forecasting yang tepat digunakan untuk data yang memiliki efek trend dan seasonal

iklim_hoW <- HoltWinters(iklim_train)

4.2 Pemodelan dengan Seasonal ARIMA (SARIMA)

Sebelum dilakukan pemodelan menggunakan seasoanl Arima, akan dilakukan pengujian untuk mengetahui kestasioneran data. Adapun hipotesis yang diberikan adalah :

H0 : data tidak stasioner H1 : data stasioner

adf.test(x = iklim_train)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  iklim_train
#> Dickey-Fuller = -3.5499, Lag order = 10, p-value = 0.0375
#> alternative hypothesis: stationary

Karena p-value yang didapatkan < 0.05 kita dapat menyimpulkan bahwa tersebut telah stasioner.

tsdisplay(iklim_train, lag.max = 36)

Fungsi tsdisplay dapat kita gunakan untuk menentukan kemungkinan ordo ARIMA yang dapat dibentuk dengan deskripsi sebagai berikut:

  • tails off: garis saat lag ke-1 mengalami penurunan yang lambat.
  • cuts off lag: garis saat lag ke-1 mengalami penurunan yang cepat.

Selain menentukan ordo ARIMA secara manual, R menyediakan fitur * ‘auto.arima’ yang dapat digunakan untuk menentukan ordo model secara otomatis. Parameter yang dibutuhkan adalah:

  • y: peubah yang digunakan untk pemodelan.
  • seasonal: tipe pemodelan yang digunakan. Seasonal = TRUE untuk seasonal ARIMA
iklim_model_auto <- auto.arima(y = iklim_train, seasonal = TRUE)
iklim_model_auto
#> Series: iklim_train 
#> ARIMA(5,1,1)(0,1,0)[365] 
#> 
#> Coefficients:
#>          ar1      ar2      ar3     ar4     ar5      ma1
#>       0.7803  -0.0786  -0.0562  0.0953  0.0267  -0.9876
#> s.e.  0.0377   0.0469   0.0470  0.0474  0.0379   0.0074
#> 
#> sigma^2 estimated as 113.5:  log likelihood=-2760.45
#> AIC=5534.89   AICc=5535.05   BIC=5567.03

Model yang terbentuk berdasarkan auto.arima adalah ARIMA(5,1,1)(0,1,0)[365]

4.3 Pemodelan dengan STLM (Seasonal Trend with Loess Model)

STL secara konsep akan melakukan smoothing terhadap data tetangga setiap masing-masing observasi dengan memberikan bobot yang lebih berat terhadap data yang dekat dengan observed data. Kekurangan dari STL hanya bisa melakukan decompose pada additive data, apabila terdapat multiplicative data dapat menggunakan transformasi log()[^12].

Untuk memodelkan hasil STL, kita bisa menerapkan STLM(Seasonal Trend with Loess Model) dimana kita bisa menerapkan metode exponential smoothing (ETS) dan ARIMA. Selain itu, STLM dapat digunakan sebagai alternative cara untuk menangkap seasonal yang belum bisa ditangkap oleh metode ETS dan ARIMA biasa.

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
iklim_stlm <- stlm(y = iklim_train, s.window = 365, method = "ets")

5 Evaluasi Model

Setelah mendapatkan kedua model tersebut, langkah selanjutnya adalah mengevaluasi model terbaik dalam peramalan data. Proses yang akan dilakukan adalah menggunakan forecast masing-masing model dan diuji menggunakan data test.

5.1 Model Holt’s Winters Exponential

iklim_how_forecast <- forecast(object = iklim_hoW, h = 365)
iklim_ts %>% autoplot(series = "train data")+
  autolayer(iklim_test, series = "test data") +
  autolayer(iklim_how_forecast$mean, series = "forecast")

accuracy(iklim_how_forecast$mean, x = iklim_test)
#>                ME     RMSE      MAE     MPE     MAPE      ACF1 Theil's U
#> Test set 4.735448 14.04932 11.04442 5.39748 19.78077 0.6992682  1.873943

5.2 Model Seasonal ARIMA (SARIMA)

iklim_auto_forecast <- forecast(object = iklim_model_auto, h = 365)
plot(iklim_auto_forecast)

iklim_ts %>% autoplot(series = "train data")+
  autolayer(iklim_test, series = "test data") +
  autolayer(iklim_auto_forecast$mean, series = "forecast")

accuracy(iklim_auto_forecast$mean, x = iklim_test)
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -5.299036 16.67471 13.23241 -15.67275 27.89623 0.7263831  2.839585

5.3 Model STLM (Seasonal Trend with Loess Model)

iklim_stlm_forecast <- forecast(iklim_stlm, h = 365)
iklim_ts %>% autoplot(series = "train data")+
  autolayer(iklim_test, series = "test data") +
  autolayer(iklim_stlm_forecast$mean, series = "forecast")

accuracy(iklim_stlm_forecast$mean, x = iklim_test)
#>                ME    RMSE      MAE      MPE     MAPE      ACF1 Theil's U
#> Test set 1.031508 12.1086 9.561429 -2.16612 18.15752 0.6923795  1.763546

Berdasarkan MAPE yang terdapat pada ketiga model, dapat disimpulkan pemodelan menggunakan STLM adalah yang terbaik dengan indikator paramater MAPE terkecil yakni 18.15%.

5.4 Uji Asumsi

Langkah selanjutnya adalah melakukan pengujian asumsi yang terdiri dari uji normalitas dan ketiadaan autokorelasi terhadap peramalan yang telah dibuat

5.4.1 Uji Normalitas

  1. Normality: Shapiro.test

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

shapiro.test(iklim_stlm_forecast$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  iklim_stlm_forecast$residuals
#> W = 0.9928, p-value = 3.671e-05
hist(iklim_stlm_forecast$residuals, breaks = 100)

plot(iklim_stlm_forecast$residuals)

5.4.2 Uji Autokorelasi

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

Box.test(iklim_stlm_forecast$residuals, type = "Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  iklim_stlm_forecast$residuals
#> X-squared = 1.8305, df = 1, p-value = 0.1761

Berdasarkan pengujian dua asumsi tersebut terdapat pemenuhan asumsi untuk ketiadaan autokorelasi namun masih terdapat asumsi normalitas yang dilanggar. Pelanggaran normalitas bisa terlihat dari frekuensi histogram yang tidak tersebar pada nilai rata-rata.

Salah satu teknik yang bisa dilakukan untuk mengatasi masalah normalitas tersebut adalah dengan melakukan seasonal adjustment. Teknik tersebut berusaha untuk menghilangkan efek musiman pada data.

6 Saran

Pemodelan selanjutnya dapat mempertimbangkan teknik pemodelan yang lain seperti SARIMA dengan ordo model manual, melakukan seasonal adjustment yakni dengan menghilangkan faktor musimam sehingga murni melihat perilaku data berdasarkan trendnya. Terdapat pula 3 faktor lain yang dapat diteliti yakni suhu rata-rata, kecepatan angin, atau tekanan udara.