Seasonal Autoregressive Integrated Moving Average (SARIMA)dengan R [pengunjung museum Avila Adobe]

SARIMA?

Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan pengembangan dari model Autoregressive Integrated Moving Average (ARIMA) pada data runtun waktu yang memiliki pola musiman. Metode ini dipopulerkan oleh George Box dan Gwilym Jenskins sekitar tahun 1970-an, model ini telah banyak dipelajari secara luas dan mengadopsi salah satu model yaitu ARIMA model. [jika ingin mengenal analisis ARIMA klik disini]

Adapun Notasi SARIMA adalah sebagai berikut:

dimana (p, d, q) merupakan bagian yang tidak musiman dari model, (P, D, Q) merupakan bagian musiman dari model dan s merupakan jumlah periode musiman (periodecity).

Bentuk umum dari model SARIMA adalah sebagai berikut:

Model SARIMA melibatkan langkah-langkah berikut:

1. Mengidentifikasi struktur SARIMA (p, d, q) (P, D, Q)^s
2. Mengestimasi parameter yang tidak diketahui
3. Melakukan uji Goodness of Fit terhadap estimasi residual
4. Melakukan peramalan dari data yang diketahui

DATA:

Data yang digunakan adalah data jumlah pengunjung museun Avila Adobe di Los Angeles dari Januari 2014 — Maret 2020 yang diperoleh dari https://www.kaggle.com/

Data dapat di download di link berikut:

https://www.kaggle.com/cityofLA/los-angeles-museum-visitors

Goals:

Melakukan Prediksi Untuk 5 bulan kedepan, tepatnya bulan april — agustus 2020

Kali ini peramalan menggunakan bantuan program R. Maka pastikan program sudah terinstall pada laptop kalian !

Langkah 1: Statistika Deskriptif

Pertamma masukkan data yang telah di download pada program r, disini peneliti menggunakan RStudio

pengunjung=read.csv("D:\\KULIAH\\semester 6\\P data mining\\UAS medium\\sarima pengunjung museum//museum-visitors.csv",sep = ";")
View(pengunjung)

Kemudian ubah data ke bentuk data time series dan disimpan pada objek “pengunjung”.

pengunjung=ts(pengunjung$Avila.Adobe, start=c(2014,1), freq=12)

Gambarkan plot dari jumlah pengunjung museum Avila Adobe Los Angeles.

ts.plot(pengunjung, col="blue", main="TS:pengunjung museum Avila Adobe Los Angeles", lwd=2)

Dari gambar diatas dapat diketahui bahwa jumlah pengunjung dalam enam tahun terakhir cukup fluktuatif dan membentuk pola musiman, setiap diawal tahun (tepatnya dibulan ferbruari) terjadi penurunan. Dan juga diketauhi hampir disetiap bulan kelima (mei) menjadi puncak pengunjung.

Melihat Statistika desktiprif dari data secara ringkas.

summary(pengunjung)

Dalam kurun waktu 6 tahun tersebut, rata-rata jumlah pengunjung sebanyak 22.559 orang.

Langkah 2: Melakukan uji stasioneritas

Dibutuhkan packages berikut.

library(forecast)
library(tseries)

Berdasarkan plot data sebelumnya, kita ketahui data mengandung musiman. Hal tersebut bisa menjadi petunjuk bahwa data tidak stasioner. Namun untuk membuktikannya lebih lanjut dapat menggunakan uji stasioneritas dari plot ACF (Autocorrelation Function) , data belum stasioner ketikan grafik terlihat meluruh secara melambat menuju 0. sepeti kasus berikut (contoh).

Kemudian stasioneritas data juga dapat dilihat dari uji ADF (Augmented Dickey-Fuller). Dengan H0: data mengandung unit root (tidak stasioner). daerah kritis tolak H0 jika p-value < alfa.

Pada praktik kali ini didapatkan hasil sebagai berikut, uji ADF (Augmented Dickey-Fuller).

#uji stasioneritas data
adf.test(pengunjung) #H0: data mengandung unit root (tidak stasioner)

Hipotesis

H0= Data mengandung unit root (yang mengakibatkan data tidak stasioner)

H1= Data tidak mengandung unit root (yang mengakibatkan data stasioner

Keputusan: Tolak H0 karena p-value <α yaitu 0,01< 0.05

Kesimpulan: Dengan menggunakan tingkat kepercayaan sebesar 95%, data yang ada mendukung hipotesis nol, yang artinya data mengandung unit root yang mengakibatkan data tidak stasioner.

Dilihat dari plot ACF

Acf(pengunjung,lag.max = 36)

Plot tidak meluruh lambat menuju 0. Sehingga dari uji ACF dan plot ACF data sudah stasioner (tidak perlu dilakukan differencing)

Langkah 3: Estimasi Model

Selanjutnya dilakukan estimasi model dapat dlihat dari plot ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function) sebagai berikut.

#identifikasi model
par(mfrow=c(2,1))
Pacf(pengunjung, lag.max=36)
Acf(pengunjung,lag.max = 36)

Grafik ACF dan PACF pengunjung musium

Plot di atas digunakan untuk menentukan order dalam model SARIMA (p, d, q), (P,D,Q). Pertama menentukan (p,d,q) nilai p atau AR diperoleh dari jumlah lag 4 data awal grafik PACF yang melewati rentang atau batasan rata-rata atau pada grafik digambarkan dengan garis putus-putus berwarna biru. diketahui pada plot ACF menunjukan bahwa terlihat data meluruh menuju 0 secara eksponensial sehingga q=0. Pada plot PACF terdapat satu lag yang keluar makan p=1.

Kedua menentukan (P,D,Q) ini menunjukan musiman suatu data. Nilai P atau SAR dilihat dari lag ke 12, lag 24, lag 36 grafik PACF yang melewati rentang atau batasan rata-rata atau pada grafik digambarkan dengan garis putus-putus berwarna biru. Begitu juga dengan melihat nilai P pada grafik ACF yaitu lag ke 12, 24, 36. Maka diperoleh Q=1 karena lag ke 12 keluar. Pada plot PACF tidak terdapat lag keluar sehingga P=0.

note: nilai diffrencing = 0 karen diawal kita tidak melakukan diffrencing untuk men-stasioner kan data.

Sehingga didapatkan model awal SARIMA (1,0,0) (0,0,1). Namun untuk medapatkan model terbaik dilakukan overfitting model disekitar model utama. Berikut hasil uji kecocokan model secara keseluruhan.

Peneliti hanya menggunakan 3 model.

printstatarima <- function (x, digits = 4,se=TRUE){
  if (length(x$coef) > 0) {
    cat("\nCoefficients:\n")
    coef <- round(x$coef, digits = digits)
    if (se && nrow(x$var.coef)) {
      ses <- rep(0, length(coef))
      ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
      coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
      coef <- rbind(coef, s.e. = ses)
      statt <- coef[1,]/ses
      pval  <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
      coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
      coef <- t(coef)
    }
    print.default(coef, print.gap = 2)
  }
}

fitting model

#fitting model
model1=Arima(log(pengunjung), order=c(0,0,1), seasonal = list(order=c(0,0,1), period=12),
             include.mean = FALSE)
summary(model1)
printstatarima(model1)
model2=Arima(log(pengunjung),order=c(1,0,0), seasonal = list(order=c(1,0,0), periode=12),
             include.mean = FALSE)
summary(model2)
printstatarima(model2)
model3=Arima(log(pengunjung),order=c(1,0,1), seasonal = list(order=c(1,0,0), periode=12),
             include.mean = FALSE)
summary(model3)
printstatarima(model3)

dari model 1, diketahui yang signifikan yaitu ar dan sma. namun peneliti coba mencari model selanjutnya dengan memasukan sar=1 dan menurunkan sma=0

pada model 2 terlihat semua koefisiennya signifikan. Kemudian dicoba lagi dengan menambahkan ar=1.

pada model 3 terlihat ma1 tidak signifikan.

Dari hasil identifikasi model, maka didapatkan ringkasan sebagai berikut.

Berdasarkan tabel diatas model signifikan pada p-value < α =0,05. Maka, didapatkan ada 2 model yang signifikan semua koefisien modelnya.

Langkah 4: Uji Diagnostik

Kemudian dilakukan uji diagnostik residual data, terdapat tiga uji diagnostik residual data yaitu uji homoskedastisitas, uji autokorelasi, dan uji normalias,. Berikut hasil dari uji diagnostik residual model SARIMA (1,0,0) (0,0,1) sebagai berikut:

#uji diagnostik 
par(mfrow=c(1,1))
tsdiag(model1)

Berdasarkan plot gambar diatas dapat dilihat pada plot ACF residual data merupakan model white noise ditandai dengan tidak adanya lag yang keluar dari garis batas interval. Sedangkan dilihat dari plot p-value for Ljung-Box statistic, nilai p-value ada yang lebih dari α=0,05 yang artinya tidak terdapat autokorelasi pada data sehingga asumsi no-autokorelasi terpenuhi.

Maka model cocok digunakan untuk peramalan.

Langkah 5: Peramalan

Sebagai catatan: Jika model yang signifikan ada lebih dari satu, maka gunakan model dengan nilai AIC terkecil.

Selanjutnya melakukan permalan pengunjung untuk 5 hari ke depan, gunakan sintaks berikut.

#peramalan 
pengunjung.pred=predict(model1,n.ahead=5)
#membalikkan ke data awal
pengunjung.pred2=exp(pengunjung.pred$pred)
pengunjung.pred2

Jangan lupa karena tadi kita menggunakan data tranformasi “log”, maka kita harus mengembalikannya ke data semula menggunakan tranformasi “exp”.

maka didapatkan hasil sebagai berikut.

Berdasarkan tabel diatas menunjukkan hasil peramalan jumlah pengunjung 5 hari ke depan yaitu mulai bulan april 2020 — agustus 2020. Terlihat hasil peramalan cukup berfluktuasi, ada pengurangan jumlah pengunjung pada bulan juni dan agustus.

hasil prediksi

#nilai prediksi
pengunjung.fitted2=exp(fitted(model1))
pengunjung.fitted2

Kemudian data aktual, fitting, dan peramalan ditampilan dalam bentuk grafik.

### Grafik data asli dan data prediksi
ts.plot(pengunjung, col="blue", xlim=c(2014,2021), main="prediction of passengers for the 5 periods")
lines(pengunjung.fitted2, col="red")
lines(pengunjung.pred2,col="green")

garis biru adalah data aktual, garis merah data prediksi, dan hijau data ramalan.

Langkah 6: Akurasi Peramalan

Dapat melihat nilai keakuratan peramalan dengan menggunakan metode Mean Absobute Percentage Error (MAPE).

mape= mean(abs(pengunjung-(pengunjung.fitted2))/(pengunjung))*100
mape

didapatkan nilai akurasi sebesar 80,78%. Berdasarkan nilai MAPE, peramalan metode Seasonal Autoregressive Moving Average (SARIMA) secara umum, menunjukan bahwa hasil peramalan yang cukup baik/akurat, karena memililki akurasi > 80%.