Time Series Modelling stationary ARIMA step by step



ARIMA (Autoregressive integrated moving average)
Penjeasan tentang ARIMA selengkapnya bisa lihat di wikipedia kilik link ini.

Tahapan ARIMA modeling:
  1. Exploratory data analysis
  2. Dekomposisi data
  3. Test stationarity
  4. Membuat Model ARIMA
  5. Forecasts

Library:
library(ggfortify)
library(tseries)
library(forecast)

Data
  • Disini kita pergunakan data Air Passanger ada dalam package R
  • Sekarang kita periksa apakah dataset AirPassangers sudah format time series?
data(AirPassengers)
AP <- AirPassengers
class(AP)
  • Hasilnya dataset AirPassangers sudah format time series
## [1] "ts"

STEP-1  EXPLORATORY DATA ANALYSIS (EDA)

  • Mari kita lihat dulu sepertia apa datasetnya
AP
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
  • Pemeriksaan nilai yang hilang (missing value)
sum(is.na(AP))
  • Hasilnya tidak ada missing value
## [1] 0
  • Pemeriksaan Frequency dengan fungsi frequency()
frequency(AP)
  • Hasilnya ada frequency=12
## [1] 12
  • Pemeriksaan cycle dengan fungsi cycle()
cycle(AP)
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949   1   2   3   4   5   6   7   8   9  10  11  12
## 1950   1   2   3   4   5   6   7   8   9  10  11  12
## 1951   1   2   3   4   5   6   7   8   9  10  11  12
## 1952   1   2   3   4   5   6   7   8   9  10  11  12
## 1953   1   2   3   4   5   6   7   8   9  10  11  12
## 1954   1   2   3   4   5   6   7   8   9  10  11  12
## 1955   1   2   3   4   5   6   7   8   9  10  11  12
## 1956   1   2   3   4   5   6   7   8   9  10  11  12
## 1957   1   2   3   4   5   6   7   8   9  10  11  12
## 1958   1   2   3   4   5   6   7   8   9  10  11  12
## 1959   1   2   3   4   5   6   7   8   9  10  11  12
## 1960   1   2   3   4   5   6   7   8   9  10  11  12
  • Summary data AirPassangers 
summary(AP)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   104.0   180.0   265.5   280.3   360.5   622.0
  • Plot data AirPassangers
  • Sumbu horisontal = "Periode"
  • Sumbu vertikal = "Penumpang (ribuan)"
  • Judul = "Air Passenger tahun 1949 sd 1961"
plot(AP,xlab="Periode", ylab = "Penumpang (ribuan)",main="Air Passenger tahun 1949 sd 1961")

  • Pemeriksaan Seasonal dengan fungsi boxplot()
boxplot(AP~cycle(AP),xlab="Periode", 
        ylab = "Penumpang (ribuan)" ,
        main ="Air Passenger tahun 1949 sd 1961")

Hasil EDA:

  • Jumlah penumpang bertambah setiap tahunnya
  • Trend meningkat secara linear
  • Disetiap bulan ke 6 dan 9, mean lebih tinggi dan variance lebih lebar karena banyak orang berlibur di US saat musim  panas
  • Data sudah bersih karena tidak ada outliers dan missing value sehingga tidak perlu cleansing data



STEP-2 Dekomposisi


  • Proses dekomposisi dllakukan untuk memriksa adanya trend, seasonal, dan random component,
  • Persamaannya:
Y[t]=T[t]S[t]e[t]
  • dimana
    • Y(t) jumlah penumpang pada waktu t,
    • T(t) trend pada waktu t,
    • S(t) seasonal pada waktu t,
    • e(t) random error pada waktu t.
  • untuk memeriksa dekomposisi lakukan fungsi berikut
dekomposisiAP <- decompose(AP,"multiplicative")
autoplot(dekomposisiAP)


Disini kita dapat melihat plot data, pola seasonal, trend yang naik, dan komponen random pada bagian  “remainder”.


STEP-3:TEST STATIONARITY


  • Time series dikatakan stationary dengan syarat mean, variance dan covariance bukan merupakan fungsi waktu.
  • Syarat model ARIMA, time series harus stationary

1. Test stationarity menggunakan Augmented Dickey-Fuller (ADF) test
  • Test stationary dilakukan dengan fungsi adf.test()
  • setting hipotesis:
    • Ho: time series bukan stationary
    • Ha: time series adalah stationary
adf.test(AP) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  AP
## Dickey-Fuller = -7.3186, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
  • Bila p-value kurang 5%, maka Ho ditolak. 
  • Karena p-value = 0.01 kurang dari 0.05 maka Ho ditolak sehingga time series adalah stationary.

2. Test Stationary menggunakan Autocorrelation
  • Test stationary yang lain adalah test autocorrelation. 
  • Kita gunakan fungsi ggAcf() ada dalam R package. 
  • Fungsi ini menghitung korelasi antara series dan lagnya dengan 95% confident interval.
  • Bila autokorelasi memotong garis putus-putus biru berarti lag tsb berkorelasi dengan seriesnya secara significant.
ggAcf(AP,
         xlab="Periode", 
         ylab = "Penumpang (ribuan)",
         main="Air Passenger tahun 1949 sd 1961")
  • Maximum pada lag bulan 1 or 12  menandakan a relasi positif dengan cycle 12 bulan.
  • Plot 
ggAcf(decomposeAP$random[7:138],
      xlab="Periode", 
      ylab = "Penumpang (ribuan)",
      main="Air Passenger tahun 1949 sd 1961")



  • Residuals tersebar di sekitar 0.

STEP-4 Membuat Model ARIMA

  • Fungsi untuk menghitung model ARIMA kita gunakan auto.arima().
arima_AP <- auto.arima(AP)
arima_AP
## Series: AP 
## ARIMA(2,1,1)(0,1,0)[12]                    
## 
## Coefficients:
##          ar1     ar2      ma1
##       0.5960  0.2143  -0.9819
## s.e.  0.0888  0.0880   0.0292
## 
## sigma^2 estimated as 132.3:  log likelihood=-504.92
## AIC=1017.85   AICc=1018.17   BIC=1029.35
  • Parameter ARIMA(2,1,1)(0,1,0)[12] 
  • Model ARIMA:
  • dimana E adalah beberapa error.
  • sekarang kita periksa apakah model ARIMA sudah fit dengan memeriksa residualnya
#library(ggfortify) gunakan library ini 
ggtsdiag(arima_AP)

  • Residual terlihat tersebar di sekitar 0 sebagai noise, tanpa pola. 
  • Model Arima ini secara fair dikatakan good fit

Part 5: FORECAST

  • forecast() dengan 95% confidence interval untuk 3 tahun (36 bulan) kedepan.
forecastAP <- forecast(arima_AP, level = c(95), h = 36)
autoplot(forecastAP)







No comments:

Post a Comment