ARIMA
Leonardo
2/22/2021
Pengolahan Data untuk Model ARIMA
Library untuk ARIMA
# Untuk membaca data dengan format Excel
library(readxl)
# Untuk mengubah data menjadi Time Series, Acf, Pacf, Model Arima dan ADF Test
library(tseries)
# Untuk melihat signifikansi dari Koefesien
library(lmtest)
# Untuk memprediksi data nantinya
library(forecast)
# Untuk Plot Data
library(ggplot2) Karakteristik Data
Pada bagian ini kita akan melihat pengolahan data Total Barang Dalam Negri yang dimuat di Tanjung Priok. Data diperoleh dari bps.go.id pada bulan Juni tanggal 13 tahun 2020. Pertama kita akan melihat datanya terlebih dahulu.
data <- read_excel("Total Barang Dalam Negeri yang Dimuat di 5 Pelabuhan Utama, 2006-2020 (Ton).xlsx",
sheet = "Tanjung Priok", col_types = c("date", "numeric"))
# Statistika Deskriptif
summary(data) ## Waktu Total Barang (Ton)
## Min. :2006-01-01 00:00:00 Min. : 383471
## 1st Qu.:2009-07-24 06:00:00 1st Qu.: 707690
## Median :2013-02-15 00:00:00 Median : 989665
## Mean :2013-02-14 13:23:43 Mean : 960547
## 3rd Qu.:2016-09-08 12:00:00 3rd Qu.:1175159
## Max. :2020-04-01 00:00:00 Max. :1894263# Mengubah data menjadi bentuk Time Series
tpriok <- ts(data$`Total Barang (Ton)`)
# Membuat data train
tpriok_train <- ts(data$`Total Barang (Ton)`[1:138])
# Membuat Grafik Data
par(mfrow=c(2,1)) # Membuat plot Data menjadi 1 bagian
plot(tpriok,lwd = 2, main = 'Plot Data Tanjung Priok')
abline(h=mean(tpriok),lwd=2,lty = 2, col ='red') # Membuat garis rataan
Acf(tpriok_train, main ='Grafik ACF Data Train Tanjung Priok', lag.max = 36)
adf.test(tpriok_train)##
## Augmented Dickey-Fuller Test
##
## data: tpriok_train
## Dickey-Fuller = -1.4635, Lag order = 5, p-value = 0.7997
## alternative hypothesis: stationaryPersiapan Pengolahan Data
Dapat dilihat dari grafik bahwa data tidak stasioner dalam rataan dan dari ADF test kita memperoleh nilai \(p_{value} < \alpha = 0.05\) sehingga dapat dikatakan bahwa data Kemudian akan coba kita diferensiasi sebanyak 2 kali kemudian akan kita lihat grafik data dan ACF data yang sudah didiferensiasi dan juga nilai \(p_{value}\) dari ADF test masing-masing
tpriok_train_diff = diff(tpriok_train)
par(mfrow=c(2,1))
plot(tpriok_train_diff,lwd = 2, main = 'Plot Data Train
Tanjung Priok 1 Kali Diferensiasi')
abline(h=mean(tpriok_train_diff),lwd=2,lty = 2, col ='red')
Acf(tpriok_train_diff, main ='Grafik ACF Data Train Tanjung Priok 1 Kali
Diferensiasi ', lag.max = 36)
Plot Data dan ACF dari Data Prediksi yang Sudah Didiferensiasikan 1 Kali
adf.test(tpriok_train_diff)##
## Augmented Dickey-Fuller Test
##
## data: tpriok_train_diff
## Dickey-Fuller = -6.8938, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarytpriok_train_diff_2 = diff(tpriok_train_diff)
par(mfrow=c(2,1))
plot(tpriok_train_diff_2,lwd = 2,
main = 'Plot Data Train Tanjung Priok 2 Kali Diferensi')
abline(h=mean(tpriok_train_diff_2),lwd=2,lty = 2, col ='red')
Acf(tpriok_train_diff_2, main ='ACF Data Train Tanjung Priok 2 Kali
Diferensi', lag.max = 36)
Plot Data dan ACF dari Data Prediksi yang Sudah Didiferensiasikan 2 Kali
adf.test(tpriok_train_diff_2)##
## Augmented Dickey-Fuller Test
##
## data: tpriok_train_diff_2
## Dickey-Fuller = -10.611, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryDapat dilihat pada 2 buah data yang sudah didiferensiasi bahwa data tersebut sudah stasioner baik ketika didiferensiasi 1 kali dan 2 kali. Meskipun sudah stasioner pada diferensiasi yang pertama, kita dapat mendiferensiasikan lagi untuk mencari model yang terbaik.
Penentuan Model
Pada bagian ini kita akan menentukan Model dari data prediksi yang sudah didiferensiasi sebanyak 1 kali dan 2 kali. Pada bagian ini kita akan melihat terlebih dahulu grafik ACF dan PACFnya baru kemudian kita akan menentukan model-model mana saja yang mungkin.
par(mfrow = c(2,2))
Acf(tpriok_train_diff, main ='ACF Data Train Tanjung Priok 1 Kali
Diferensiasi', lag.max = 36)
Pacf(tpriok_train_diff, main ='PACF Data Train Tanjung Priok 1 Kali
Diferensiasi', lag.max = 36)
Acf(tpriok_train_diff_2, main ='ACF Data Train Tanjung Priok 2 Kali
Diferensiasi', lag.max = 36)
Pacf(tpriok_train_diff_2, main ='PACF Data Train Tanjung Priok 2 Kali
Diferensiasi', lag.max = 36)
Grafik ACF dan PACF untuk Penentuan Model
Untuk menentukan model dapat kita lihat pada lag ke berapa masih signifikan (melewati batas signifikansi). Jika sudah pernah melewati maka akan kita abaikan. Perlu diingat juga bahwa kita akan menentukan model berdasarkan prinsip Parsimoni.
Dapat dilihat pada grafik ACF dan PACF pada baris pertama, bahwa ACF melewati batas signifikansi pada lag-1 dan PACF melewati batas signifikansi pada lag-3. Kemudian dapat kita lihat kembali pada grafik ACF dan PACF pada baris pertama, bahwa ACF melewati batas signifikansi pada lag-8 dan PACF melewati batas signifikansi pada lag-8 Karena kita menggunakan Prinsip Parsimoni maka dapat kita ambil model yang parameter p dan q untuk model ARIMA nya tidak lebih dari 4. Sehingga dapat kita pilih model yang masih dalam range lag yang dilewati. Tapi untuk pengolahan data kali ini kita akan mencoba menggunakan model ARIMA sebagai berikut :
- \(ARIMA (3,1,1)\)
- \(ARIMA (0,1,0)\)
- \(ARIMA (0,2,0)\)
- \(ARIMA (3,2,3)\)
Model tidak terbatas dari 4 model ini saja. Tapi lebih dari ini, Untuk penentuan model dapat kita ambil nilai parameter p dan q = 0 jika pada lag pertama ACF/PACFnya nilainya negatif. Karena yang diharapkan nilai korelasinya positif. Tapi jika ingin diambil juga tida ada masalah.
Pada bagian ini kita akan melihat nilai AIC dari model yang kita miliki. Dan akan kita pilih model yang memiliki AIC yang paling rendah.
mod_1 = arima(tpriok_train,order=c(3,1,1))
mod_1##
## Call:
## arima(x = tpriok_train, order = c(3, 1, 1))
##
## Coefficients:
## ar1 ar2 ar3 ma1
## -0.0766 -0.1057 -0.2465 -0.5646
## s.e. 0.1446 0.1090 0.0966 0.1347
##
## sigma^2 estimated as 2.256e+10: log likelihood = -1827.79, aic = 3665.59mod_2 = arima(tpriok_train,order=c(0,1,0))
mod_2##
## Call:
## arima(x = tpriok_train, order = c(0, 1, 0))
##
##
## sigma^2 estimated as 3.356e+10: log likelihood = -1854.6, aic = 3711.2mod_3 = arima(tpriok_train,order=c(0,2,0))
mod_3##
## Call:
## arima(x = tpriok_train, order = c(0, 2, 0))
##
##
## sigma^2 estimated as 9.614e+10: log likelihood = -1912.64, aic = 3827.27mod_4 = arima(tpriok_train,order=c(3,2,3))
mod_4##
## Call:
## arima(x = tpriok_train, order = c(3, 2, 3))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## -0.7645 0.0122 -0.1828 -0.8130 -0.7969 0.6137
## s.e. 0.1407 0.1750 0.1079 0.1328 0.1658 0.1213
##
## sigma^2 estimated as 2.179e+10: log likelihood = -1815.56, aic = 3645.13mod_auto = auto.arima(tpriok_train, max.p=4,max.q=4,
seasonal =FALSE, stationary = TRUE)
mod_auto## Series: tpriok_train
## ARIMA(4,0,3) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 mean
## 0.7565 0.5499 -0.7161 0.3856 -0.3510 -0.6125 0.4553 831424.8
## s.e. 0.2085 0.1644 0.2017 0.1288 0.2125 0.1502 0.1554 195267.6
##
## sigma^2 estimated as 2.254e+10: log likelihood=-1837.87
## AIC=3693.75 AICc=3695.15 BIC=3720.09Dapat dilihat bahwa model \(ARIMA(3,2,3)\) memiliki nilai AIC yang paling rendah. Tapi karna prinsip parsimoni maka yang akan digunakan adalah model yang paling sederhana yakni \(ARIMA(3,1,1)\) Pada bagian ini juga dikenalkan syntax auto.arima yang dapat digunakan untuk menghemat waktu/referensi tapi tidak disarankan untuk langsung digunakan. Tidak ada jaminan bagi model dari auto.arima merupakan model yang terbaik.
Signifikansi dari Koefisien Data dan Pembuatan Model
coeftest(mod_1) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.076618 0.144612 -0.5298 0.59624
## ar2 -0.105661 0.108985 -0.9695 0.33230
## ar3 -0.246494 0.096608 -2.5515 0.01073 *
## ma1 -0.564618 0.134660 -4.1929 2.754e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Dapat dilihat bahwa \(p_{value}\) dari koefisien ar1, ar2 dan ar3 > \(\alpha\) maka dapat kita simpulkan bahwa hanya koefisien ar1 dan ar2 yang tidak siginifikan pada model yang kita miliki. Sehingga model kita hanya memiliki parameter ma 1 Sehingga kita memiliki model sebagai berikut \[\phi_3(B) (1-B)^1 Y_t = \theta_1(B) \epsilon_t\] \[(1-\phi_1B-\phi_2B^2-\phi_3B^3)(1-B)Y_t=(1-\theta_1B)\epsilon_t\] \[(1-\phi_1B-\phi_2B^2-\phi_3B^3-B+\phi_1B^2+\phi_2B^3+\phi_3B^4)Y_t=(1-\theta_1B)\epsilon_t\] \[(1-(1+\phi_1)B-(\phi_2-\phi_1)B^2-(\phi_3-\phi_2)B^3+\phi_3B^4)Y_t = \epsilon_t - \theta_1\epsilon_{t-1}\] \[Y_t-(1+\phi_1)Y_{t-1}-(\phi_2-\phi_1)Y_{t-2}-(\phi_3-\phi_2)Y_{t-3}+\phi_3Y_{t-4}=\epsilon_t - \theta_1\epsilon_{t-1}\] \[Y_t=(1+\phi_1)Y_{t-1}+(\phi_2-\phi_1)Y_{t-2}+(\phi_3-\phi_2)Y_{t-3}-\phi_3Y_{t-4}+\epsilon_t - \theta_1\epsilon_{t-1}\] Karena kita sudah mendapat nilai koefisien dari \(\phi_1,\phi_2, \phi_3\) dan \(\theta_1\) maka akan kita masukkan ke persamaan yang kita miliki, tapi kita hanya akan memasukkan parameter-parameter yang signifikan saja. Sehingga persamaan kita menjadi
\[Y_t=\epsilon_t+0.564618\epsilon_{t-1}\] Persamaan ini mirip dengan persamaan \(MA(1)\) tapi nilai koefisiennya pasti akan berbeda. Sehingga tidak sama dengan memasukkan model \(MA(1)\) langsung.
Diagnostik Model / Model Diagnostic
Pada bagian ini kita akan melihat residual dari model yang kita buat.
checkresiduals(mod_1)
Model Diagnostik dari Data Prediksi
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,1,1)
## Q* = 14.281, df = 6, p-value = 0.02665
##
## Model df: 4. Total lags used: 10 Dapat dilihat pada grafik diatas bahwa residu dari model kita bisa dikatakan stasioner, sedikit berkorelasi dan mengikut distribusi normal. Sehingga dapat dikatakan bahwa residu dari model ini sudah mengikuti white noise.
## Pemeriksaan Data pada data Test Pada bagian ini akan kita bangun secara manual data sisa 20% yang dijadikan data test. Kita akan mengambil nilai error dari model. Sehingga error saat t=1 adalah 0. Kemduian akan kita cari nilai MAPEnya
tpriok_pred <- tpriok_train - residuals(mod_1)
ts.plot(tpriok_pred,tpriok_train, xlab = 'Waktu', ylab = 'Total Barang',
col = c('red', 'blue'), main = 'Perbandingan Antara Data Asli dan dari Model')
accuracy(mod_1)## ME RMSE MAE MPE MAPE MASE
## Training set 15254.52 149641.3 109133.6 -0.04054429 12.54541 0.8626635
## ACF1
## Training set -0.01647759Dapat dilihat bahwa nilai MAPEnya adalah 23% sehingga dapat dikatakan bahwa model kita sangat baik. Berikutnya akan kita lihat hasil prediksi kita dalam bentuk grafik.
Forecast
Untuk memprediksi data di depan dapat digunakan syntax forecast. Perlu diperhatikan data apa yang akan di prediksi dan modelnya. Jangan sampai menggunakan model dari data asli ataupun menggunakan data untuk membangun model
fc = forecast(tpriok, model=mod_1,h=1)
summary(fc)##
## Forecast method: ARIMA(3,1,1)
##
## Model Information:
## Series: object
## ARIMA(3,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ma1
## -0.0766 -0.1057 -0.2465 -0.5646
## s.e. 0.0000 0.0000 0.0000 0.0000
##
## sigma^2 estimated as 2.256e+10: log likelihood=-2311.27
## AIC=4624.54 AICc=4624.56 BIC=4627.68
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 11935.14 178427.7 128932.6 -0.9497287 13.88677 0.8753136
## ACF1
## Training set 0.02047757
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 173 1015983 823511.7 1208455 721623.4 1310343plot(fc)
Grafik Prediksi Ke Depannya
Kesimpulan
Model dari pengolahan data Total Barang di Pelabuhan Tanjung Priok diperoleh bahwa model yang kita miliki adalah \(ARIMA (3,1,1)\) dan memiliki MAPE sebesar \(\approx 12\%\) sehingga model ini dapat digunakan untuk memprediksi data ke depannya tapi akan lebih baik jika kita memiliki data yang lebih banyak.