1. Intro

Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan model dalam analisis runtun waktu (time series). Metode SARIMA merupakan perkembangan dari metode ARIMA, dimana metode ini pertama kali diperkenalkan oleh BoxJenkins pada tahun 1970. Musiman adalah kecenderungan mengulangi pola tingkah gerak dalam periode musim, biasanya satu tahun untuk data bulanan. Model SARIMA merupakan model ARIMA yang digunakan untuk menyelesaikan time series musiman yang terdiri dari dua bagian, yaitu bagian tidak musiman (non-musiman) dan bagian musiman. Bagian non-musiman dari metode ini adalah model ARIMA. Secara umum bentuk modelnya adalah ARIMA musiman atau ARIMA (p, d, q)(P, Q, S) S.

Musiman didefinisikan sebagai suatu pola yang berulang-ulang dalam selang waktu yang tetap. Notasi umum SARIMA adalah :

\[SARIMA (p,d,q)(P,D,Q)_s\] dimana :

(p,d,q) = Bagian yang tidak musiman dari model

(P,D,Q) = Bagian musiman dari model

s = Jumlah periode per musim

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 yang digunakan dalam ilustrasi penerapan SARIMA adalah data jumlah penumpang penerbangan domestik dari suatu Bandara X. Berdasarkan dari data tersebut, dilakukan peramalan 12 periode kedepan.

2. Input Data

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)

penumpang=read.delim("clipboard")
penumpang

Dikonversi menjadi data runtun waktu.

penumpang=ts(penumpang,start=c(2008,1),frequency=12)
penumpang=penumpang[,-1]
penumpang
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2008 328446 282242 310829 275959 281806 281428 313235 298881 221823 340406
## 2009 331528 289383 319658 311808 337949 356683 391103 381083 308402 433271
## 2010 388819 349148 337764 405263 422810 441067 476358 367383 455251 477450
## 2011 457763 413489 436203 429024 448132 518562 529225 386770 548400 504393
## 2012 530692 494799 544229 522512 542413 528568 602625 528977 651697 586731
## 2013 624398 518487 585913 560613 590532 627701 512111 658784 682269 685386
## 2014 617838 479197 538497 510996 554213 609753 455747 743304 604342 636143
## 2015 565027 474994 501031 513301 568271 518583 608491 699259 545042 605010
## 2016 656208 571726 618357 596085 675702 523451 825715 792232 649375 718451
## 2017 672153 546280 613925 615053 618375 539747 889548 676016 696308 655005
##         Nov    Dec
## 2008 294919 309608
## 2009 427065 417994
## 2010 449044 474367
## 2011 516875 521433
## 2012 607904 608329
## 2013 593909 624290
## 2014 593137 644533
## 2015 587291 671396
## 2016 699285 693048
## 2017 666423
ts.plot(penumpang, col='blue', main="Jumlah Penumpang Penerbangan Domestik", lwd=2)

Dapat dilihat bahwa pola datanya mengandung tren naik dan memiliki variasi data musiman. Karena data mengandung musiman, maka hal tersebut menunjukkan data tidak stasioner. Untuk melihatnya, dapat dilakukan uji ADF (Augmented Dickey-Fuller Test) dan KPSS test.

3. Seasonal Autoregressive Integrated Moving Average (SARIMA) Method

3.1 Stationary Data

### Uji stasioneritas data ### 
adf.test(penumpang) 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  penumpang
## Dickey-Fuller = -3.2741, Lag order = 4, p-value = 0.07891
## alternative hypothesis: stationary
kpss.test(penumpang, null="Level")
## Warning in kpss.test(penumpang, null = "Level"): p-value smaller than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  penumpang
## KPSS Level = 2.1533, Truncation lag parameter = 4, p-value = 0.01

Dengan uji ADF, diperoleh p-value sebesar 0.07891, karena p-value = 0.07891 > α = 0.05, maka disimpulkan bahwa data yang ada gagal menolak H0 yang artinya data mengandung unit root sehingga mengakibatkan data tidak stasioner. Sedangkan dengan uji KPSS, diperoleh p-value sebesar 0.01, artinya data tidak mengandung unit root sehingga mengakibatkan data tidak stasioner.

Ketidakstasioneran data juga dapat dibuktikan oleh grafik ACF dan PACF.

### Melihat adanya musiman (selain dari plot) ###
par (mfrow=c(1,2))
Acf(penumpang, lag.max = 36)
Pacf(penumpang, lag.max = 36)

Berdasarkan plot di atas, dapat dilihat bahwa data grafik ACF menurun perlahan-lahan mendekati nilai 0 serta bentuk grafiknya bergelombang. Hal tersebut berarti data mengandung musiman dan tidak stasioner. Jika data tidak stasioner, maka perlu dilakukan penyesuaian untuk menghasilkan data yang stasioner. Salah satu cara yang umum dipakai adalah metode pembedaan (differencing).

3.2 Differensi Musiman Orde 1

Dilakukan transformasi musiman yaitu dengan diferensi musiman orde 1.

### Jika data mengandung musiman maka, dilakukan differencing musiman ###
penumpang.dslog = diff(log(penumpang), differences=1, lag=12)
ts.plot(penumpang.dslog, col="red", main="Penumpang Penerbangan Domestik dslog")

### Uji stationeritas data differencing musiman ###
adf.test(penumpang.dslog)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  penumpang.dslog
## Dickey-Fuller = -2.7533, Lag order = 4, p-value = 0.2642
## alternative hypothesis: stationary

Pada plot data differensi musiman orde 1, data tidak berfluktuasi disekitar rata-rata sehingga data diferensi musiman orde 1 tidak stasioner.

Dengan uji ADF, juga diperoleh kesimpulan bahwa data yang ada gagal menolak H0 yang artinya data mengandung unit root sehingga mengakibatkan data tidak stasioner karena p-value = 0.2642.

3.3 Differensi Non-Musiman Orde 1

Karena data masih belum stasioner, maka data harus didefferensi kembali yaitu dengan differensi non-musiman orde 1.

### Data belum stasioner maka, dilakukan differensing lagi yang non-musiman ###
penumpang.ddslog =  diff(penumpang.dslog, differences=1)
ts.plot(penumpang.ddslog, col="red", main="Penumpang Penerbangan Domestik ddslog")

### Uji stationeritas data differencing musiman dan non-musiman ###
adf.test(penumpang.ddslog)
## Warning in adf.test(penumpang.ddslog): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  penumpang.ddslog
## Dickey-Fuller = -6.2622, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Pada plot data differensi non-musiman orde 1, karena data telah berfluktuasi disekitar rata-rata, maka data diferensi non-musiman orde 1 telah stasioner. Uji ADF juga menunjukkan bahwa data yang ada menolak H0 yang artinya data tidak mengandung unit root sehingga mengakibatkan data stasioner.

3.4 Identifikasi Model SARIMA

Karena data sudah stasioner, maka selanjutnya mengidentifikasi model SARIMA dengan menggambar plot ACF dan PACF dari data hasil differensi non-musiman orde 1.

### Melakukan estimasi model dengan melihat plot ACF dan PACF ###
par(mfrow=c(1,2))
Acf(penumpang.ddslog, lag.max=37)
Pacf(penumpang.ddslog, lag.max=37)

Pada model SARIMA notasi yang digunakan dalam membuat model yaitu (p,d,q) sebagai notasi non musiman dan (P,D,Q) sebagai notasi musimannya. Output lag yang digunakan pada grafik ACF dan PACF hanya 4 lag pertama.

Berdasarkan plot tersebut, batang PACF keluar hingga lag ke-2 dan pada lag 12, 24, 36 tidak ada lag yang keluar maka menunjukkan AR yaitu p = 2 dan SAR yaitu P = 0. Pada batang ACF keluar hingga lag ke-1 dan pada lag 12, 24, 36 keluar hingga lag ke-3 maka menunjukkan order MA yaitu q = 1 dan Q = 3. Serta sebelumnya dilakukan differensi musiman orde 1 (d = 1) dan non-musiman orde 1 (D = 1).

Dengan begitu diperoleh model utama SARIMA (2,1,1)(0,1,3). Overfitting terhadap model dapat dipilih model dengan order lebih rendah atau kombinasi order pada model utama. Kali ini akan diuji 6 model, yaitu model SARIMA (2,1,1)(0,1,3), SARIMA (2,1,1)(0,1,2), SARIMA (0,1,1)(0,1,2), SARIMA (1,1,0)(0,1,2), SARIMA (1,1,1)(0,1,0), SARIMA (2,1,1)(0,1,0).

3.5 Estimasi Parameter Model SARIMA

Selanjutnya dilakukan estimasi parameter untuk keenam model untuk melihat signifikansi dari koefisien model.

### Fitting Model ###
model1 = Arima(log(penumpang), order = c(2,1,1), seasonal = list(order = c(0,1,3), period = 12), include.mean = FALSE)
summary(model1)
## Series: log(penumpang) 
## ARIMA(2,1,1)(0,1,3)[12] 
## 
## Coefficients:
##           ar1      ar2      ma1     sma1     sma2     sma3
##       -0.3473  -0.1252  -0.4463  -0.0484  -0.5287  -0.1783
## s.e.   0.2478   0.1784   0.2399   0.1795   0.2199   0.1603
## 
## sigma^2 = 0.008465:  log likelihood = 100.05
## AIC=-186.09   AICc=-184.95   BIC=-167.45
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE     MASE
## Training set -0.002509932 0.08434072 0.06415786 -0.01987401 0.4868681 0.513831
##                     ACF1
## Training set -0.00398209
model2 = Arima(log(penumpang), order = c(2,1,1), seasonal = list(order = c(0,1,2), period = 12), include.mean = FALSE)
summary(model2)
## Series: log(penumpang) 
## ARIMA(2,1,1)(0,1,2)[12] 
## 
## Coefficients:
##           ar1      ar2      ma1     sma1     sma2
##       -0.3807  -0.1475  -0.4454  -0.0177  -0.5783
## s.e.   0.2247   0.1675   0.2185   0.1397   0.1883
## 
## sigma^2 = 0.008577:  log likelihood = 99.26
## AIC=-186.51   AICc=-185.67   BIC=-170.53
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set -0.002797484 0.08532263 0.06429591 -0.02212178 0.4878298 0.5149366
##                      ACF1
## Training set -0.005544792
model3 = Arima(log(penumpang), order = c(0,1,1), seasonal = list(order = c(0,1,2), period = 12), include.mean = FALSE)
summary(model3)
## Series: log(penumpang) 
## ARIMA(0,1,1)(0,1,2)[12] 
## 
## Coefficients:
##           ma1     sma1     sma2
##       -0.7056  -0.0342  -0.6392
## s.e.   0.0651   0.1597   0.2094
## 
## sigma^2 = 0.008477:  log likelihood = 97.48
## AIC=-186.96   AICc=-186.57   BIC=-176.31
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set -0.004225837 0.08565799 0.06327716 -0.03282879 0.4799424 0.5067776
##                    ACF1
## Training set -0.1161501
model4 = Arima(log(penumpang), order = c(1,1,0), seasonal = list(order = c(0,1,2), period = 12), include.mean = FALSE)
summary(model4)
## Series: log(penumpang) 
## ARIMA(1,1,0)(0,1,2)[12] 
## 
## Coefficients:
##           ar1     sma1     sma2
##       -0.5631  -0.0076  -0.5115
## s.e.   0.0824   0.1523   0.2248
## 
## sigma^2 = 0.01018:  log likelihood = 90.52
## AIC=-173.03   AICc=-172.64   BIC=-162.38
## 
## Training set error measures:
##                        ME      RMSE        MAE         MPE      MAPE      MASE
## Training set -0.001624374 0.0938753 0.06829058 -0.01370462 0.5176222 0.5469294
##                    ACF1
## Training set -0.2084378
model5 = Arima(log(penumpang), order = c(1,1,1), seasonal = list(order = c(0,1,0), period = 12), include.mean = FALSE)
summary(model5)
## Series: log(penumpang) 
## ARIMA(1,1,1)(0,1,0)[12] 
## 
## Coefficients:
##           ar1      ma1
##       -0.3117  -0.5361
## s.e.   0.1263   0.1145
## 
## sigma^2 = 0.0103:  log likelihood = 92.81
## AIC=-179.62   AICc=-179.39   BIC=-171.63
## 
## Training set error measures:
##                        ME       RMSE      MAE         MPE     MAPE      MASE
## Training set -0.001713167 0.09489911 0.070025 -0.01370127 0.532239 0.5608201
##                     ACF1
## Training set -0.01095692
model6 = Arima(log(penumpang), order = c(2,1,1), seasonal = list(order = c(0,1,0), period = 12), include.mean = FALSE)
summary(model6)
## Series: log(penumpang) 
## ARIMA(2,1,1)(0,1,0)[12] 
## 
## Coefficients:
##           ar1      ar2      ma1
##       -0.3843  -0.0642  -0.4678
## s.e.   0.2588   0.1942   0.2482
## 
## sigma^2 = 0.01039:  log likelihood = 92.87
## AIC=-177.73   AICc=-177.34   BIC=-167.08
## 
## Training set error measures:
##                       ME       RMSE        MAE         MPE      MAPE      MASE
## Training set -0.00164896 0.09485056 0.07003067 -0.01323821 0.5322795 0.5608655
##                      ACF1
## Training set -0.005102761
printstatarima <- function (x, digits = 4,se=T,...){
  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)
  }
}

printstatarima(model1) 
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.3473  0.2478  -1.4015  0.1637
## ar2   -0.1252  0.1784  -0.7018  0.4842
## ma1   -0.4463  0.2399  -1.8604  0.0653
## sma1  -0.0484  0.1795  -0.2696  0.7879
## sma2  -0.5287  0.2199  -2.4043  0.0178
## sma3  -0.1783  0.1603  -1.1123  0.2683
printstatarima(model2) 
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.3807  0.2247  -1.6943  0.0929
## ar2   -0.1475  0.1675  -0.8806  0.3803
## ma1   -0.4454  0.2185  -2.0384  0.0437
## sma1  -0.0177  0.1397  -0.1267  0.8994
## sma2  -0.5783  0.1883  -3.0712  0.0026
printstatarima(model3)
## 
## Coefficients:
##                  s.e.         t   sign.
## ma1   -0.7056  0.0651  -10.8387  0.0000
## sma1  -0.0342  0.1597   -0.2142  0.8308
## sma2  -0.6392  0.2094   -3.0525  0.0028
printstatarima(model4)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.5631  0.0824  -6.8337  0.0000
## sma1  -0.0076  0.1523  -0.0499  0.9603
## sma2  -0.5115  0.2248  -2.2754  0.0247
printstatarima(model5)
## 
## Coefficients:
##                 s.e.        t  sign.
## ar1  -0.3117  0.1263  -2.4679  0.015
## ma1  -0.5361  0.1145  -4.6821  0.000
printstatarima(model6)
## 
## Coefficients:
##                 s.e.        t   sign.
## ar1  -0.3843  0.2588  -1.4849  0.1402
## ar2  -0.0642  0.1942  -0.3306  0.7415
## ma1  -0.4678  0.2482  -1.8848  0.0619

Berikut hasil estimasi parameter keenam model:

Dengan menggunakan tingkat kepercayaan 95%, dapat disimpulkan bahwa model SARIMA (1,1,1)(0,1,0) menolak H0 yang artinya model tersebut signifikan.

3.6 Uji Diagnostik

Dilakukan uji diagnostik untuk melihat ada autokorelasi atau tidak. Uji ini dilakukan pada model yang signifikan saja, yaitu pada model SARIMA (1,1,1)(0,1,0).

### Uji diagnostik ###
tsdiag(model5)

jarque.bera.test(model5$residuals)
## 
##  Jarque Bera Test
## 
## data:  model5$residuals
## X-squared = 2.6573, df = 2, p-value = 0.2648

Berdasarkan output di atas, residual model SARIMA (1,1,1)(0,1,0) bersifat white noise (WN) karena tidak ada lag ke ≥ 1 yang keluar dari batas interval. Selain itu, p-value pada L-jung Box juga di atas batas 5% yang artinya residual tidak mengandung korelasi. Sehingga uji diagnostik model SARIMA (1,1,1)(0,1,0) terpenuhi.

Sehingga dapat disimpulkan bahwa model SARIMA (1,1,1)(0,1,0) adalah model yang paling tepat digunakan untuk peramalan jumlah penumpang penerbangan domestik pada 12 periode ke depan.

3.7 Peramalan

### Peramalan ###
penumpang.pred=predict(model5, n.ahead = 12)
penumpang.pred
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2017                                                                        
## 2018 13.37392 13.16710 13.28368 13.28556 13.29093 13.15494 13.65456 13.38006
##           Sep      Oct      Nov      Dec
## 2017                            13.40626
## 2018 13.40963 13.34849 13.36577         
## 
## $se
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2017                                                                      
## 2018 0.1026820 0.1110449 0.1161091 0.1217178 0.1268482 0.1318479 0.1366440
##            Aug       Sep       Oct       Nov       Dec
## 2017                                         0.1015124
## 2018 0.1412836 0.1457738 0.1501303 0.1543637
### Membalikkan ke data awal ###
penumpang.pred2 = exp(penumpang.pred$pred)
penumpang.pred2
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2017                                                                        
## 2018 643011.8 522877.2 587525.7 588636.0 591805.7 516558.6 851330.3 646972.6
##           Sep      Oct      Nov      Dec
## 2017                            664146.0
## 2018 666392.7 626864.2 637791.7

Diketahui bahwa dengan menggunakan model SARIMA (1,1,1)(0,1,0), pada Desember 2017 diramalkan jumlah penumpang penerbangan domestik adalah sebanyak 664146 orang, Januari 2018 sebanyak 643012 orang, Februari 2018 sebanyak 522877 orang, Maret 2018 sebanyak 587526 orang, April 2018 sebanyak 588636 orang, Mei 2018 sebanyak 591806 orang, Juni 2018 sebanyak 516559 orang, Juli 2018 sebanyak 851330 orang, Agustus 2018 sebanyak 646973 orang, September 2018 sebanyak 666393 orang, Oktober 2018 sebanyak 626864 orang, dan November 2018 sebanyak 637792 orang.