Analisis Runtun Waktu (ARIMA) menggunakan R

 

Jumlah Barang yang dimuat di Bandara Utama Polonia (ton)

Tahun 2008-2018

Ø  Data bulanan jumlah barang yang dimuat di Bandara Utama Polonia pada bulan Januari 2008 hingga Desember 2018.

Ø  Ø  Data in sampel :

Tahun	Bulan	Polonia
2008	Januari	260
2008	Februari	190
2008	Maret	216
2008	April	265
2008	Mei	253
2008	Juni	273
2008	Juli	260
2008	Agustus	240
2008	September	337
2008	Oktober	302
2008	November	445
2008	Desember	312
2009	Januari	950
2009	Februari	1050
2009	Maret	1144
2009	April	856
2009	Mei	811
2009	Juni	929
2009	Juli	1098
2009	Agustus	1755
2009	September	888
2009	Oktober	828
2009	November	759
2009	Desember	1028
2010	Januari	1012
2010	Februari	1417
2010	Maret	795
2010	April	765
2010	Mei	981
2010	Juni	1094
2010	Juli	1244
2010	Agustus	1298
2010	September	1017
2010	Oktober	1169
2010	November	1311
2010	Desember	1578
2011	Januari	1351
2011	Februari	1531
2011	Maret	1635
2011	April	1373
2011	Mei	1251
2011	Juni	2069
2011	Juli	1327
2011	Agustus	1262
2011	September	1079
2011	Oktober	1250
2011	November	1400
2011	Desember	1615
2012	Januari	2171
2012	Februari	1531
2012	Maret	1549
2012	April	1116
2012	Mei	1133
2012	Juni	1243
2012	Juli	1623
2012	Agustus	1263
2012	September	1317
2012	Oktober	1294
2012	November	1563
2012	Desember	1903
2013	Januari	2017
2013	Februari	1570
2013	Maret	1068
2013	April	1470
2013	Mei	1357
2013	Juni	1614
2013	Juli	1375
2013	Agustus	2036
2013	September	1315
2013	Oktober	1546
2013	November	1149
2013	Desember	1449
2014	Januari	1224
2014	Februari	1348
2014	Maret	1425
2014	April	1385
2014	Mei	1014
2014	Juni	1461
2014	Juli	1256
2014	Agustus	1545
2014	September	1646
2014	Oktober	1336
2014	November	1186
2014	Desember	1284
2015	Januari	1192
2015	Februari	926
2015	Maret	1024
2015	April	996
2015	Mei	1089
2015	Juni	1133
2015	Juli	1497
2015	Agustus	1764
2015	September	1373
2015	Oktober	1085
2015	November	1271
2015	Desember	1543
2016	Januari	1461
2016	Februari	1166
2016	Maret	1169
2016	April	1173
2016	Mei	1193
2016	Juni	1707
2016	Juli	1439
2016	Agustus	1811
2016	September	1883
2016	Oktober	1519
2016	November	1542
2016	Desember	1235
2017	Januari	2155
2017	Februari	2003
2017	Maret	2362
2017	April	1910
2017	Mei	1733
2017	Juni	1470
2017	Juli	1665
2017	Agustus	1960
2017	September	1429
2017	Oktober	1206
2017	November	1615
2017	Desember	2236
2018	Januari	2186
2018	Februari	1723
2018	Maret	1582
2018	April	1410
2018	Mei	1562
2018	Juni	1287
2018	Juli	1644
2018	Agustus	1689
2018	September	1622
2018	Oktober	1448
2018	November	1639
2018	Desember	1915

Ø  Ø  Data out sampel

Tahun	Bulan	Polonia
2019	Januari	1705
2019	Februari	1100
2019	Maret	1181
2019	April	1036
2019	Mei	1227
2019	Juni	938
2019	Juli	3412
2019	Agustus	1283
2019	September	1228
2019	Oktober	1261
2019	November	1271
2019	Desember	1435

Analisis menggunakan R studio

> library(forecast) > library(FitAR) > library(lmtest) > library(tseries) > bandara <- read.csv(file = "C:/Users/ME/Documents/Semester 4 /Analisis Runtun Waktu/bandara.csv", header=TRUE) > #data barang bandara utama > bandara <- read.csv(file = "C:/Users/ME/Documents/Semester 4 /Analisis Runtun Waktu/bandara.csv", header=TRUE) > bandara > str(bandara) 'data.frame':  132 obs. of  7 variables: $ Tahun         : int  2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ... $ Bulan         : Factor w/ 13 levels "Agustus","April" ,..: 5 4 9 2 10 7 6 1 13 12 ... $ Polonia       : int  260 190 216 265 253 273 260 240 337 302 ... $ Soekarno.Hatta: int  7913 7166 13060 10367 10172 6856 9615 10045 12056 9930 ... $ Juanda        : int  689 520 749 506 715 658 644 601 517 669 ... $ Ngurah.Rai    : int  2052 1627 3111 2761 2506 1867 1496 1700 2185 2353 ... $ Hasanudin     : int  NA NA NA NA NA NA NA NA NA NA ...

 

Ø  Plot asli

> ts.plot(bandara$Polonia,main="Jumlah barang Polonia" , ylab= "Ton")

Dari plot time series, dapat dijelaskan bahwa data tidak stasioner, maka dilakukan transformasi log dan differencing untuk mendapatkan data runtun waktu yang stasioner. Hasilnya adalah sebagai berikut.

Ø  Plot ACF

> acf(bandara$Polonia)

 

Data plot ACF juga dapat dilihat bahwa grafik yang menurun perlahan-lahan menandakan data belum stasioner.

Ø  Uji Unit Root Test dengan ADF test

> adf.test(bandara$Polonia)        Augmented Dickey-Fuller Test   data:  bandara$Polonia Dickey-Fuller = -3.0067, Lag order = 5, p-value = 0.1579 alternative hypothesis: stationary

·    Hipotesis

H₀ : Data memuat unit root (tidak stasioner)

H₁ : Data tidak memuat unit root (stasioner)

·     Taraf signifikansi :  α = 0.05

·     Statistik uji : ADF test

·     Kriteria keputusan : H₀ ditolak jika p-value < α = 0.05

·     Kesimpulan : Oleh karena p-value = 0.1579 > α = 0.05, maka H₀ diterima sehingga dapat disimpulkan bahwa data memuat unit root (tidak stasioner).

Ø  Plot setelah Differencing

> ts.plot(dbarang) > ts.plot(dlbarang)

Ø  Uji Unit Root Test dengan ADF test setelah Differencing

> dlbarang <- diff(log(bandara$Polonia)) > dlbarang > dbarang <- diff(bandara$Polonia) > dbarang > adf.test(dbarang)           Augmented Dickey-Fuller Test   data:  dbarang Dickey-Fuller = -7.2231, Lag order = 5, p-value = 0.01 alternative hypothesis: stationary   Warning message: In adf.test(dbarang) : p-value smaller than printed p-value > adf.test(dlbarang)           Augmented Dickey-Fuller Test   data:  dlbarang Dickey-Fuller = -5.6883, Lag order = 5, p-value = 0.01 alternative hypothesis: stationary   Warning message: In adf.test(dlbarang) : p-value smaller than printed p-value

      

·   Hipotesis

H₀ : Data memuat unit root (tidak stasioner)

H₁ : Data tidak memuat unit root (stasioner)

 ·   Taraf signifikansi : α = 0.05

 ·  Statistik uji : ADF test

 ·   Kriteria keputusan : H₀ ditolak jika p-value < α = 0.05

 ·  Kesimpulan : Oleh karena p-value = 0.01 < α = 0.05, maka H₀ ditolak sehingga dapat disimpulkan bahwa data tidak memuat unit root (stasioner)

Ø  Plot ACF dan PACF berdasarkan data stasioner

> acf(dbarang) > pacf(dbarang)

Dengan membandingkan grafik dan garis batas, terlihat bahwa lag ke-1 dan lag ke-3 pada plot ACF signifikan, adapun pada plot PACF lag ke-1 dan lag ke-3 bersifat signifikan.

Ø Dari langkah sebelumnya diduga model ARIMA (1,1,1). Oleh karena itu, dalam melakukan pendugaan atau estimasi parameter, lakukan overfitting ke model ARIMA (1,1,2) dan ARIMA (2,1,1).

> m1 <- arima (Polonia,c(1,1,1)) > summary(m1)   Call: arima(x = Polonia, order = c(1, 1, 1))   Coefficients:          ar1      ma1       0.3401  -0.8141 s.e.  0.1110   0.0615   sigma^2 estimated as 77337:  log likelihood = -923.42,  aic = 1852.85   Training set error measures:                  ME     RMSE      MAE       MPE     MAPE      MASE Training set 38.649 277.0391 211.6909 0.7674744 16.12313 0.8839292                     ACF1 Training set -0.01387533 > coeftest(m1)   z test of coefficients:        Estimate Std. Error  z value  Pr(>|z|)    ar1  0.340118   0.111002   3.0641  0.002184 ** ma1 -0.814128   0.061527 -13.2320 < 2.2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   > m2 <- arima (Polonia,c(2,1,1)) > summary(m2)   Call: arima(x = Polonia, order = c(2, 1, 1))   Coefficients:          ar1      ar2      ma1       0.3380  -0.0315  -0.8037 s.e.  0.1118   0.0978   0.0712   sigma^2 estimated as 77268:  log likelihood = -923.37,  aic = 1854.75   Training set error measures:                    ME     RMSE      MAE       MPE     MAPE      MASE Training set 38.43688 276.9161 211.1871 0.7577198 16.08928 0.8818256                     ACF1 Training set -0.02954234 > coeftest(m2)   z test of coefficients:        Estimate Std. Error  z value Pr(>|z|)    ar1  0.337988   0.111750   3.0245  0.00249 ** ar2 -0.031487   0.097763  -0.3221  0.74740    ma1 -0.803729   0.071155 -11.2954  < 2e-16 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   > m3 <- arima (Polonia,c(1,1,2)) > summary(m3)   Call: arima(x = Polonia, order = c(1, 1, 2))   Coefficients:          ar1      ma1      ma2       0.3059  -0.7766  -0.0276 s.e.  0.2104   0.2024   0.1407   sigma^2 estimated as 77312:  log likelihood = -923.41,  aic = 1854.81   Training set error measures:                   ME     RMSE      MAE       MPE     MAPE      MASE Training set 38.57883 276.9943 211.5023 0.7643518 16.10928 0.8831415                     ACF1 Training set -0.01994682 > coeftest(m3)   z test of coefficients:        Estimate Std. Error z value  Pr(>|z|)    ar1  0.305908   0.210427  1.4537 0.1460163    ma1 -0.776581   0.202358 -3.8377 0.0001242 *** ma2 -0.027609   0.140722 -0.1962 0.8444574    --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   > m4 <- arima (Polonia,c(3,1,2)) > summary(m4)   Call: arima(x = Polonia, order = c(3, 1, 2))   Coefficients:           ar1     ar2      ar3      ma1      ma2       -0.0335  0.0230  -0.3296  -0.4102  -0.1612 s.e.   0.2503  0.1506   0.0941   0.2635   0.2362   sigma^2 estimated as 70607:  log likelihood = -917.6,  aic = 1847.21   Training set error measures:                    ME     RMSE     MAE       MPE     MAPE      MASE Training set 34.20985 264.7117 203.666 0.7190253 15.63094 0.8504204                     ACF1 Training set -0.02049689 > coeftest(m4)   z test of coefficients:        Estimate Std. Error z value  Pr(>|z|)    ar1 -0.033474   0.250278 -0.1337 0.8936014    ar2  0.022996   0.150616  0.1527 0.8786524    ar3 -0.329627   0.094093 -3.5032 0.0004597 *** ma1 -0.410165   0.263507 -1.5566 0.1195745    ma2 -0.161164   0.236157 -0.6824 0.4949577    --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > m5 <- arima (Polonia,c(3,1,1)) > summary(m5)   Call: arima(x = Polonia, order = c(3, 1, 1))   Coefficients:          ar1      ar2      ar3      ma1       0.1145  -0.0538  -0.3405  -0.5704 s.e.  0.1496   0.0959   0.0931   0.1535   sigma^2 estimated as 70810:  log likelihood = -917.79,  aic = 1845.57   Training set error measures:                   ME     RMSE      MAE       MPE     MAPE      MASE Training set 32.7519 265.0917 203.2646 0.5891464 15.62125 0.8487447                      ACF1 Training set -0.008775029 > coeftest(m5)   z test of coefficients:        Estimate Std. Error z value  Pr(>|z|)    ar1  0.114483   0.149636  0.7651 0.4442255    ar2 -0.053762   0.095923 -0.5605 0.5751568    ar3 -0.340451   0.093111 -3.6564 0.0002558 *** ma1 -0.570434   0.153499 -3.7162 0.0002022 *** --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Ø Diagnostic Check

Uji Satistik Ljung-Box

> tsdiag(m1) > tsdiag(m2) > tsdiag(m3)   > tsdiag(m4) > tsdiag(m5)

Kesimpulan: Hasil di atas menunjukkan bahwa residual model ARIMA(3,1,2) dan model ARIMA (3,1,1) yang telah memenuhi syarat white noise. Hal ini ditunjukkan oleh p‐value dari uji Ljung‐Box yang semua lagnya berada di atas garis batas biru.

Oleh karena AIC lebih kecil, BIC/SBC lebih kecil, dan uji diagnostik terpenuhi, maka model terbaik yang dapat dipilih adalah ARIMA (3,1,1).

Ø Estimasi Model ARIMA

> auto.arima(Polonia) Series: Polonia ARIMA(3,1,1) with drift   Coefficients:          ar1      ar2      ar3      ma1    drift       0.1538  -0.0407  -0.3319  -0.6327  11.3727 s.e.  0.1421   0.0959   0.0947   0.1434   7.0803   sigma^2 estimated as 72290:  log likelihood=-916.63 AIC=1845.27   AICc=1845.95   BIC=1862.52

Kesimpulan: Berdasarkan output di atas, maka model ARIMA(3,1,1) yang diperoleh dapat ditulis secara matematis seperti berikut.

(1 + 0.1538 B – 0.0407 B² – 0.3319 B³ ) (1 – B) Y_t = (1 – 0.6327 B) a_t ,

atau
Y_t = 0.8462 Y_t-1 + 0.1945 Y_t-2 + 0.2912 Y_t-3 – 0.3319 Y_t-4 + a_t  – 0.6327 a_{t -1} ,

dengan Yt adalah data asli pada waktu ke‐t.

Ø  Forecasting

> forecast(m5)     Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95 133       1813.814 1472.791 2154.837 1292.2648 2335.364 134       1722.366 1334.140 2110.591 1128.6259 2316.106 135       1623.372 1207.393 2039.351  987.1863 2259.558 136       1651.404 1233.736 2069.072 1012.6360 2290.172 137       1691.069 1265.885 2116.253 1040.8067 2341.332 138       1727.806 1290.389 2165.222 1058.8345 2396.777 139       1720.335 1260.693 2179.977 1017.3736 2423.297 140       1704.001 1226.319 2181.683  973.4488 2434.553 141       1690.026 1197.788 2182.263  937.2137 2442.838 142       1691.847 1188.529 2195.165  922.0885 2461.606 > plot(forecast(m5))

 

Kesimpulan: Hasil dari forecasting untuk 10 bulan ke depan tahun 2019 adalah sebagai berikut.

Ø  Perbandingan dengan data out sampel

Jika dibandingkan dengan data out sampel tahun 2019 hasil peramalan tidak begitu mendekati dengan data.