1 Intro

Data berikut merupakan data perubahan iklim di Delhi dari tahun 2013-2016, India. Kita akan mencoba memprediksi temperature di Delhi menggunakan data time series berikut.

2 Setup

Set up chunk di awal untuk mengatur format chunk pada markdown

options(scipen = 9999)
rm(list=ls())

Setup library yang akan digunakan

library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # imputasi missing value TS
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz

library(forecast) # algoritma forecasting
library(TTR) # SMA function
library(tseries) # adf.test

# supaya semua plot memiliki theme_minimal()
theme_set(theme_minimal())

3 Data Preprocessing

3.1 Input data dan Struktur data

Impor data dari dataset Daily delhi climate

delhi <- read.csv("DailyDelhiClimateTrain.csv")
str(delhi)
## 'data.frame':    1461 obs. of  5 variables:
##  $ date        : chr  "1/1/2013" "1/2/2013" "1/3/2013" "1/4/2013" ...
##  $ meantemp    : num  10 7.4 7.17 8.67 6 ...
##  $ humidity    : num  84.5 92 87 71.3 86.8 ...
##  $ wind_speed  : num  0 2.98 4.63 1.23 3.7 ...
##  $ meanpressure: num  1016 1018 1019 1017 1016 ...

3.2 Data Wrangling

Dalam LBB ini saya hanya akan menggunakan variabel date dan meantemp yang merupakan rata rata temperature udara. Lalu variabel date yang sebelumnya character, diubah menjadi tipe date.

delhi <- delhi %>% 
  select(date, meantemp) %>% 
  mutate(date = mdy(date)) %>% 
  arrange(date)

head(delhi)

Setelah itu kita cek apakah sudah memenuhi syarat time series yaitu :

  1. Data tidak boleh ada yang NA
  2. Data harus terurut berdasarkan waktu
  3. Data tidak boleh ada waktu yang bolong
delhi %>% 
  pad() %>%
  anyNA()
## pad applied on the interval: day
## [1] FALSE

Hasil cek menyatakan bahwa data diatas tidak memiliki data NA, dan data date sudah berurutan dan tidak ada yang hilang, langkah selanjutkan adalah mengubah tipe data menjadi tipe data ts.

4 Exploratory Data Analysis

4.1 Time Series

Sebelum melakukan EDA, ubah data menjadi time series (ts) dengan frekuensi tahunan

delhi_ts <- ts(data = delhi$meantemp, start = c(2013,1), frequency = 365)

4.2 Cross Validation

Sebelum melakukan pemodelan, kita akan memisahkan data menjadi data train dan data test dengan data test sebanyak 4 bulan atau 120 hari

delhi_test <-  delhi_ts %>% tail(120)
delhi_train <- delhi_ts %>% head(-120)

4.3 Decompose

Lalu kita lakukan EDA dengan menggunakan fungsi decompose() yang memiliki fungsi untuk memecah komponen ts menjadi beberapa komponen, yaitu :

  1. Trend adalah pergerakan mean secara global
  2. Seasonal adalah pola data per frequency
  3. Error adalah nilai yang tidak dapat oleh Trend dan Seasonal
delhi_dc <- delhi_train %>% 
  decompose()

autoplot(delhi_dc)

4.4 Seasonal Analysis

Apabila ingin melihat pola atau seasonal

autoplot(delhi_dc$x - delhi_dc$trend)

4.5 Trend Analysis

Apabila ingin melihat trend naik atau turun dengan menggunakan adjusted seasonal yaitu membuang efek seasonal data.

autoplot(delhi_dc$x - delhi_dc$seasonal)

Dari hasil plot diatas, dapat dikatakan bahwa data iklim delhi memiliki pola seasonal dan tidak memiliki trend. Selanjutnya dilanjutkan dengan membangun model

5 Model Building

Karena data memiliki pola seasonal dan trend maka saya akan mencoba model HoltWinters dan ARIMA, dan akan dibandingkan mana model yang lebih baik saat evaluasi.

5.1 Holt-Winters

Untuk model holtwinters, kita akan set beta menjadi False karena data tidak memiliki trend.

delhi_hw <- HoltWinters(delhi_train, beta = F)

5.2 SARIMA

Pada hasil decompose diatas, data kita terdapat seasonal, sehingga metoda ARIMA yang digunakan adalah SARIMA, Seasonal ARIMA adalah metode ARIMA di mana data time series memiliki pola seasonal. Kita akan menggunakan metode SARIMA dengan fungsi stlm()

delhi_arima <- stlm(y = delhi_train,
                   method = "arima")

summary(delhi_arima$model)
## Series: x 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.6873  -0.9829
## s.e.  0.0225   0.0068
## 
## sigma^2 estimated as 1.734:  log likelihood=-2270.23
## AIC=4546.47   AICc=4546.49   BIC=4562.07
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.03516891 1.315521 1.005676 -0.1237221 4.013427 0.3918013
##                    ACF1
## Training set 0.00754072

6 Forecast

Setelah kita selesai membuat model, kita akan melakukan forecast pada kedua model tersebut

6.1 Forecast Holt-Winters

hw_forecast <- forecast(object = delhi_hw, h = 120)
hw_forecast %>% 
  autoplot() +
  autolayer(delhi_test)

6.2 Forecast SARIMA

arima_forecast <- forecast(object = delhi_arima, h = 120)
arima_forecast %>% 
  autoplot() +
  autolayer(delhi_test)

7 Evaluasi

7.1 Model Forecast Comparison

Gunakan custom function compare_forecast() untuk membandingkan performa masing-masing model di data test

compare_forecast <- function(x, test){
  lapply(x, function(x) forecast::accuracy(x, test)["Test set", ]) %>%
    lapply(., function(x) x %>% t() %>% as.data.frame) %>% 
    bind_rows() %>% 
    mutate(model = names(x)) %>%
    select(model, everything())}

forecast_list <- list(
  "Holt_winters" = hw_forecast,
  "SARIMA" = arima_forecast)

compare_forecast(forecast_list, delhi_test)

Pada hasil evaluasi, MAPE pada model SARIMA lebih rendah hampir 10% yang berarti model SARIMA dapat memprediksi iklim lebih tepat dibanding model Holt-Winters

7.2 Visualisasi

Menggunakan fungsi test_forecast(), kita akan melihat bagaimana visualisasi keakuratan model SARIMA dalam memprediksi data test

test_forecast(actual = delhi_ts, 
              forecast.obj = arima_forecast, 
              train = delhi_train, 
              test = delhi_test)

8 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. Kita akan menggunakan tes pada residual data, mengapa menggunakan residual data? Karena kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model.

Metode forecasting yang baik menghasilkan nilai residual berikut ini:

  1. Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
  2. Residual berdistribusi normal di sekitaran nilai 0.

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

8.1 No auto-correlation residual

Untuk melakukan tes asumsi NO-Auto Correlation, lebih pasti jika melakukan uji Ljung-box dengan hipotesis:

  • \(H_0\): tidak terdapat autokorelasi pada residual
  • \(H_1\): terdapat autokorelasi pada residual

Harapan: p-value > alpha (0.05)

Box.test(delhi_arima$residuals, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  delhi_arima$residuals
## X-squared = 0.076423, df = 1, p-value = 0.7822

8.2 Normality residual

Untuk melakukan tes asumsi normality, lebih pasti jika melakukan uji shapiro test menggunakan fungsi shapiro.test(model$residuals) dengan hipotesis:

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

Harapan: p-value > alpha (0.05)

shapiro.test(delhi_arima$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  delhi_arima$residuals
## W = 0.98983, p-value = 0.00000004998

9 Conclusion

  • Dari hasil evaluasi model yang telah dilakukan, dapat disimpulkan bahwa model SARIMA merupakan model yang terbaik dibandingkan dengan model Holt-Winters
  • Dari hasil uji asumsi, tidak ada auto korelasi pada residual (p-value >0.05), tetapi residual forecast masih tidak terdistribusi secara normal.