 |
ARIMA Intervensi dengan R
|
Halo teman-teman, berjumpa lagi dengan blog sederhana ini. Pada beberapa unggahan sebelumnya, kita telah belajar dan berbagai bersama mengenai pemodelan Autoregressive Integrated and Moving Average (ARIMA). Kali ini kita akan coba ulas bersama secara singkat mengenai tetangganya model ARIMA. Bukan tetangga sebenarnya, tetapi lebih ke saudara dari model ARIMA. Karena model satu ini memang merupakan hasil pengembangan dari ARIMA.
Model yang satu ini kita namakan sebagai ARIMA Intervensi atau bisa juga disebut ARIMA dengan intervensi. ARIMA Intervensi ini adalah bentuk Intervention Analysis, namun bentuknya berbeda dengan model yang akhir-akhir ini berkembang di tengah masyarakat, yaitu Causal Impact Time Series Model. ARIMA Intervensi merupakan model ARIMA dengan variabel eksogen berupa error sebenarnya, bentuknya mirip ARIMAX. Kalau ARIMAX variabel eksogennya merupakan variabel lain yang tidak berasal dari dirinya sendiri.
Baik, sebelum menuju praktik pemodelan menggunakan R, kita akan ulas terlebih dahulu mengenai pengertian atau ARIMA Intervensi secara teoritis. ARIMA Intervensi ini merupakan alat analisis data runtun waktu (time series) di mana sebuah data runtun waktu dipengaruhi oleh kejadian eksternal yang menyebabkan adanya perubahan pada pola data selama periode waktu tertentu. Menurut Morokimban, dkk. (2021), ARIMA Intervensi merupakan bentuk intervention analysis untuk mengolah data deret waktu yang dipengaruhi oleh peristiwa di luar kendali atau yang disebut intervensi. Menurut Damayanti dan Siska (2021), ARIMA Intervensi merupakan model intervensi data time series yang digunakan untuk pemodelan sekaligus peramalan data yang mengandung intervensi baik dari faktor internal maupun eksternal. Dari beberapa pengertian berdasarkan kajian di atas, bisa kita tarik benang merah bahwa ARIMA Intervensi merupakan bentuk model analisis intervensi pada data runtun waktu (time series) dengan memasukkan intervensi baik sebagai faktor internal maupun sebagai faktor eksternal.
Intervensi yang dimaksud dalam ARIMA Intervensi atau analisis intervensi ini banyak kita temui di sekitar kita. Bisa berupa kebijakan pemerintah, bencana alam, atau kejadian lain yang membuat sebuah data runtun waktu secara nyata berubah pola, baik perubahan itu menjadi sebuah pola baru yang permanen, atau bersifat sementara atau singkat. Sebagai contoh yang nyata dan masih berlangsung adalah pandemi Covid-19 yang menerpa dunia sejak awal 2020 lalu. Dari data, efek dari pandemi ini terlihat begitu mengubah polanya. Terhadap jumlah penumpang kereta api harian misalkan, sebelum pandemi Covid-19, mungkin datanya terlihat normal tanpa fluktuasi yang berarti, tetapi sejak pandemi Covid-19 menerpa, pola data jumlah penumpang kereta api itu pastinya akan berubah drastis atau signifikan menurun. Ini kalau intervensi itu dipahami sebuah kejadian alam.
Berbeda bila intervensi yang disebabkan oleh sebuah kebijakan pemerintah, misalkan Pembatasan Sosial Berskala Besar (PSBB) yang sekitar kuartal atau triwulan II - 2020 mulai diterapkan di sejumlah wilayah di Indonesia khususnya. Tentu efeknya terhadap perubahan pola data runtun waktu akan berbeda, bila efek yang dimaksud adalah dampak secara alami pandemi Covid-19 dengan dampak yang ditimbulkan akibat kebijakan pemerintah dengan penerapan PSBB atau lockdown. Efek yang ditimbulkan oleh intervensi secara alami pandemi Covid-19 akan terlihat lebih smooth perubahan pola data runtun waktunya dibandingkan efek kebijakan PSBB atau lockdown yang kita saksikan bersama menimbulkan shock atau kejutan tajam terhadap data.
Demikian halnya ketika pemulihan pasca pandemi Covid-19, efek yang alami akibat pandemi, dengan asumsi tanpa kebijakan pemerintah, secara shock tidak begitu tajam, tetapi waktu untuk pemulihan akan semakin lama karena terjadi secara alami. Berbeda dengan efek kebijakan merespon pandemi dengan penerapan SPBB dan lockdown sejumlah wilayah, secara tajam dan tiba-tiba membuat data runtun waktu menjadi berubah drastis namun pemulihan data menuju normal semakin cepat seiring dengan penanganan serius terhadap pandemi Covid-19. Ini kita bicara sekilas mengenai intervensi yang terjadi secara alami dan intervensi buatan berupa kebijakan pemerintah.
Secara umum, ARIMA Intervensi itu sendiri bisa kita maknai sebagai model yang menggabungkan formulasi ARIMA dengan intervensi yang mengikuti sebuah sebaran fungsi tertentu. Untuk lebih tergambar, bisa kita lihat pada ilustrasi matematis berikut:
 |
Persamaan umum ARIMA dengan Intervensi fungsi Pulse dan Step
|
Dari ilustrasi di atas, selanjutnya kita akan berkenalan dengan dua jenis fungsi intervensi yang populer digunakan dalam penelitian dengan alat analisis ARIMA Intervensi, yaitu fungsi Pulse dan fungsi Step. Fungsi Pulse ditandai sebagai sebuah dummy variable yang bernilai 1 ketika sebuah intervensi terjadi pada waktu T, kemudian efek intervensi itu lenyap atau hilang dan data runtun waktu kembali normal seperti sedia kala. Beberapa peristiwa yang diduga dapat kita nyatakan mengikuti sebaran fungsi Pulse misalnya peristiwa Bom Bali I, Bom Bali II, Gunung Merapi meletus, Gunung Krakatau meletus, pemadaman listrik sehari, dan sejenisnya. Sedangkan fungsi Step ditandai sebagai sebuah dummy variable yang bernilai 1 saat kejadian, peristiwa, atau intervensi terjadi dan efeknya gradual naik atau turun secara perlahan sejak waktu T dan efeknya bisa sementara atau bisa juga bersifat permanen menjadikan data runtun waktu berpola baru. Beberapa identifikasi pola intervensi diilustrasikan pada gambar berikut:
 |
Pola residual penentuan order fungsi intervensi Pulse dan Step
|
Berikutnya, kita beranjak mengangkat sebuah topik yang mungkin menarik untuk dijadikan pemicu ide penelitian. Kali ini kita akan membahas mengenai penerapan ARIMA Intervensi atau model intervensi untuk peramalan laju pertumbuhan ekonomi Indonesia triwulan II hingga triwulan IV tahun 2022.
Meramal pertumbuhan ekonomi Indonesia triwulan II hingga IV tahun 2022 ini sebenarnya telah dilakukan berbagai pihak. Oleh karena itu, kita juga mau mencoba melakukan peramalan kira-kira berapa persen pertumbuhan ekonomi Indonesia pada triwulan II 2022 sampai triwulan IV 2022 tahun ini? Menurut Laporan World Economic Outlook (WEO) International Monetary Fund (IMF), dalam periode 2021-2023 ekonomi Indonesia diramalkan akan tumbuh kuat sebesar 3,3 persen, 5,6 persen, dan 6,0 persen. Menurut Asian Development Bank (ADB), ekonomi Indonesia pada 2022 diperkirakan akan tumbuh 5,2 persen. Sedangkan menurut Bank Indonesia (BI), ekonomi Indonesia pada tahun 2022 akan berada pada kisaran 4,7 persen hingga 5,5 persen.
Menanggapi beberapa peramalan atau prediksi pertumbuhan ekonomi Indonesia selama 2022 dan 2023 tersebut, kita juga akan mencoba melakukannya melalui pemodelan analisis intervensi ARIMA. Adapun efek intervensi yang terakomodir pada residual peramalan ARIMA sebelum intervensi memiliki pola yang mirip ilustrasi berikut sehingga kita akan tentukan order intervensi yang tentatif adalah b = 0, r = 2, dan s = 0.
 |
Pola efek residual Intervensi Step - 34 pertumbuhan ekonomi Indonesia
|
Intervensi dari pemodelan yang akan kita gunakan adalah kebijakan pemerintah memberlakukan Pembatasan Sosial Berskala Besar (PSBB) yang dilakukan sekitar triwulan II - 2020 lampu. Sedangkan data runtun waktu yang akan kita gunakan adalah data laju pertumbuhan ekonomi Indonesia triwulanan year on year (y on y) bersumber dari laman Badan Pusat Statistik Republik Indonesia (bps.go.id) periode 2012 hingga triwulan I - 2022. Datanya telah saya siapkan dan dapat teman-teman unduh pada tautan berikut. Untuk melakukan pemodelan, kita dapat menggunakan beberapa code sebagai beriku:
#Import data laju pertumbuhan ekonomi Indonesia Triwulan 1/2012 - Triwulan I/2022 year on year
library(readxl)
lpeind <- read_excel("C:/Users/Joko Ade/Downloads/lpe_ind.xlsx")
head(lpeind, 10)
## # A tibble: 10 x 1
## Y
## <dbl>
## 1 6.5
## 2 6.3
## 3 6.4
## 4 6.17
## 5 6.11
## 6 6.02
## 7 5.81
## 8 5.62
## 9 5.72
## 10 5.21
#Mengattach data
attach(lpeind)
## The following object is masked from lpeind (pos = 3):
##
## Y
## The following object is masked from lpeind (pos = 16):
##
## Y
#Melihat sekilas data
str(lpeind)
## tibble [41 x 1] (S3: tbl_df/tbl/data.frame)
## $ Y: num [1:41] 6.5 6.3 6.4 6.17 6.11 6.02 5.81 5.62 5.72 5.21 ...
#konversi ke data runtun waktu
#frekuensi = 4 karena triwulanan
lpe_y <- ts(Y, start = c(2012,1), frequency = 4)
str(lpe_y)
## Time-Series [1:41] from 2012 to 2022: 6.5 6.3 6.4 6.17 6.11 6.02 5.81 5.62 5.72 5.21 ...
lpe_y
## Qtr1 Qtr2 Qtr3 Qtr4
## 2012 6.50 6.30 6.40 6.17
## 2013 6.11 6.02 5.81 5.62
## 2014 5.72 5.21 5.12 5.01
## 2015 5.02 4.83 4.74 4.78
## 2016 5.15 4.94 5.21 5.03
## 2017 4.94 5.01 5.06 5.19
## 2018 5.06 5.27 5.17 5.18
## 2019 5.06 5.05 5.01 4.96
## 2020 2.97 -5.32 -3.49 -2.17
## 2021 -0.70 7.07 3.51 5.02
## 2022 5.01
#Visualisasi grafik garis
plot(lpe_y, col = "blue", lwd = 1, main = "Laju Pertumbuhan Ekonomi Indonesia \nTw I/2012 - Tw I/2022 year on year", xlab = "Periode",
ylab = "persen")
#Membubuhi grid pada grafik
grid()
#memberi garis vertikal warna putus-putus pada grafik
abline(v =2020.2, col = "red", lty=2)
#menambahkan teks keterangan pada garis vertikal merah
text(2019.2, 0, c("Efek Kebijakan PSBB"))
 |
Visualisasi 1
|
#Melihat Data Pertumbuhan yang Paling Ekstrem Bawah
which.min(lpe_y)
## [1] 34
#Data ke-34 merupakan Periode Intervensi
#Library yang dibutuhkan
library(lmtest)
library(forecast)
library(tseries)
#Periode Pra-Intervensi dari Q1 2012 - Q1 2020
pre <- window(lpe_y, start=c(2012), end=c(2020,1))
plot(pre)
 |
Visualisasi 2
|
#Uji Stasioneritas pada mean
#Uji Augmented Dickey Fuller (ADF)
#Uji Pada Level d= 0
adf.test(pre)[4]
## $p.value
## [1] 0.835727
#Uji pada Diff 1, d = 1
adf.test(diff(pre))[4]
## Warning in adf.test(diff(pre)): p-value greater than printed p-value
## $p.value
## [1] 0.99
#Uji pada Diff 2, d = 2
adf.test(diff(diff(pre))) #Stasioner di Diff 2, d = 2
##
## Augmented Dickey-Fuller Test
##
## data: diff(diff(pre))
## Dickey-Fuller = -1.7452, Lag order = 3, p-value = 0.6708
## alternative hypothesis: stationary
#Uji Philip - Perons (PP)
#uji pada level d = 0
pp.test(pre)[4]
## $p.value
## [1] 0.5120746
#uji pada Diff 1 d = 1
pp.test(diff(pre))[4] #======> Stasioner pada diff 1
## Warning in pp.test(diff(pre)): p-value smaller than printed p-value
## $p.value
## [1] 0.01
#Transformasi Box-Cox agar stasioner pada Ragam
#Mencari Lambda Opsional
require(car)
lamdafm1 <- boxCox(pre~1, family = "yjPower")
 |
Visualisasi 3
|
lamdaopt <- lamdafm1$x[which.max(lamdafm1$y)]
lamdaopt
## [1] 2
lpe_box <- forecast(as.numeric(lpeind$Y), lambda = lamdaopt)
#Data yang digunakan pemodelan arima pra Intervensi
lpe_bc <- lpe_box$fitted
#Plot pre_box
plot(lpe_bc)
 |
Visualisasi 4
|
#Uji Stasioneritas pada Mean
adf.test(lpe_bc)
##
## Augmented Dickey-Fuller Test
##
## data: lpe_bc
## Dickey-Fuller = -4.1786, Lag order = 3, p-value = 0.01249
## alternative hypothesis: stationary
pp.test(lpe_bc)
##
## Phillips-Perron Unit Root Test
##
## data: lpe_bc
## Dickey-Fuller Z(alpha) = -16.139, Truncation lag parameter = 3, p-value
## = 0.1115
## alternative hypothesis: stationary
#Cek ACF dan PACF model lpe_bc, d = 0
par(mfrow=c(1, 2))
acf(lpe_bc)
pacf(lpe_bc)
 |
Visualisasi 5
|
#dari hasil cek ACF dan PACF model yang mungkin adalah ARIMA dengan AR max 2, MA max 3, d = 0
#Karena stasioner pada level, include.mean = T
#Model ARIMA(1, 0, 0)
modpre1 <- arima(lpe_bc, order = c(1, 0, 0), include.mean = T)
coeftest(modpre1)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.766292 0.094528 8.1065 5.210e-16 ***
## intercept 4.660828 0.981590 4.7482 2.052e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre1)
 |
Visualisasi 6
|
jarque.bera.test(modpre1$residuals)
##
## Jarque Bera Test
##
## data: modpre1$residuals
## X-squared = 389.47, df = 2, p-value < 2.2e-16
ks.test(modpre1$residuals, ecdf(modpre1$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre1$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre1)
 |
Visualisasi 7
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 9.7998, df = 6, p-value = 0.1333
##
## Model df: 2. Total lags used: 8
#Model ARIMA(2, 0, 0)
modpre2 <- arima(lpe_bc, order = c(2, 0, 0), include.mean = T)
coeftest(modpre2)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.91351 0.15067 6.0630 1.336e-09 ***
## ar2 -0.18714 0.15078 -1.2411 0.2146
## intercept 4.61122 0.84587 5.4514 4.996e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre2)
 |
Visualisasi 8
|
jarque.bera.test(modpre2$residuals)
##
## Jarque Bera Test
##
## data: modpre2$residuals
## X-squared = 336.47, df = 2, p-value < 2.2e-16
ks.test(modpre2$residuals, ecdf(modpre2$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre2$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre2)
 |
Visualisasi 9
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,0) with non-zero mean
## Q* = 8.8084, df = 5, p-value = 0.117
##
## Model df: 3. Total lags used: 8
#Model ARIMA(0, 0, 1)
modpre3 <- arima(lpe_bc, order = c(0, 0, 1), include.mean = T)
coeftest(modpre3)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.956935 0.090862 10.5317 < 2.2e-16 ***
## intercept 4.559979 0.495293 9.2066 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre3)
 |
Visualisasi 10
|
jarque.bera.test(modpre3$residuals)
##
## Jarque Bera Test
##
## data: modpre3$residuals
## X-squared = 304.76, df = 2, p-value < 2.2e-16
ks.test(modpre3$residuals, ecdf(modpre3$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre3$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre3)
 |
Visualisasi 11
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1) with non-zero mean
## Q* = 13.844, df = 6, p-value = 0.03143
##
## Model df: 2. Total lags used: 8
#Model ARIMA(0, 0, 2)
modpre4 <- arima(lpe_bc, order = c(0, 0, 2), include.mean = T)
coeftest(modpre4)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.36946 0.11889 3.1075 0.001887 **
## ma2 0.82551 0.15069 5.4780 4.3e-08 ***
## intercept 4.54043 0.54446 8.3393 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre4)
 |
Visualisasi 12
|
jarque.bera.test(modpre4$residuals)
##
## Jarque Bera Test
##
## data: modpre4$residuals
## X-squared = 332.21, df = 2, p-value < 2.2e-16
ks.test(modpre4$residuals, ecdf(modpre4$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre4$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre4)
 |
Visualisasi 13
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,2) with non-zero mean
## Q* = 10.08, df = 5, p-value = 0.07301
##
## Model df: 3. Total lags used: 8
#Model ARIMA(1, 0, 1)
modpre5 <- arima(lpe_bc, order = c(1, 0, 1), include.mean = T)
coeftest(modpre5)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.69845 0.13327 5.2409 1.598e-07 ***
## ma1 0.17990 0.18170 0.9901 0.3221
## intercept 4.63510 0.90055 5.1469 2.648e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre5)
 |
Visualisasi 14
|
jarque.bera.test(modpre5$residuals)
##
## Jarque Bera Test
##
## data: modpre5$residuals
## X-squared = 353.17, df = 2, p-value < 2.2e-16
ks.test(modpre5$residuals, ecdf(modpre5$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre5$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre5)
 |
Visualisasi 15
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 8.8106, df = 5, p-value = 0.1169
##
## Model df: 3. Total lags used: 8
#Model ARIMA(1, 0, 2)
modpre6 <- arima(lpe_bc, order = c(1, 0, 2), include.mean = T)
coeftest(modpre6)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.54006 0.14612 3.6961 0.000219 ***
## ma1 0.15137 0.11074 1.3669 0.171643
## ma2 0.72187 0.14776 4.8853 1.033e-06 ***
## intercept 4.53515 0.87443 5.1864 2.144e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre6)
 |
Visualisasi 16
|
jarque.bera.test(modpre6$residuals)
##
## Jarque Bera Test
##
## data: modpre6$residuals
## X-squared = 565.6, df = 2, p-value < 2.2e-16
ks.test(modpre6$residuals, ecdf(modpre6$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre6$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre6)
 |
Visualisasi 17
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2) with non-zero mean
## Q* = 6.096, df = 4, p-value = 0.1921
##
## Model df: 4. Total lags used: 8
#Model ARIMA(2, 0, 2)
modpre7 <- arima(lpe_bc, order = c(2, 0, 2), include.mean = T)
coeftest(modpre7)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.830270 0.211014 3.9347 8.331e-05 ***
## ar2 -0.363821 0.200154 -1.8177 0.06911 .
## ma1 -0.032675 0.154397 -0.2116 0.83239
## ma2 0.761140 0.129992 5.8553 4.761e-09 ***
## intercept 4.493167 0.684919 6.5601 5.376e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre7)
 |
Visualisasi 18
|
jarque.bera.test(modpre7$residuals)
##
## Jarque Bera Test
##
## data: modpre7$residuals
## X-squared = 665.39, df = 2, p-value < 2.2e-16
ks.test(modpre7$residuals, ecdf(modpre7$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre7$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre7)
 |
Visualisasi 19
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2) with non-zero mean
## Q* = 3.7047, df = 3, p-value = 0.2952
##
## Model df: 5. Total lags used: 8
#Model ARIMA(0, 0, 3)
modpre8 <- arima(lpe_bc, order = c(0, 0, 3), include.mean = T)
coeftest(modpre8)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 1.10450 0.15507 7.1225 1.060e-12 ***
## ma2 0.96164 0.16201 5.9355 2.930e-09 ***
## ma3 0.85714 0.13826 6.1996 5.659e-10 ***
## intercept 4.56896 0.70982 6.4368 1.220e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre8)
 |
Visualisasi 20
|
jarque.bera.test(modpre8$residuals)
##
## Jarque Bera Test
##
## data: modpre8$residuals
## X-squared = 1507, df = 2, p-value < 2.2e-16
ks.test(modpre8$residuals, ecdf(modpre8$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre8$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre8)
 |
Visualisasi 21
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,3) with non-zero mean
## Q* = 2.0352, df = 4, p-value = 0.7293
##
## Model df: 4. Total lags used: 8
#Model ARIMA(1, 0, 3)
modpre9 <- arima(lpe_bc, order = c(1, 0, 3), include.mean = T)
coeftest(modpre9)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.018778 0.213337 0.0880 0.9299
## ma1 1.095300 0.190945 5.7362 9.682e-09 ***
## ma2 0.950676 0.205288 4.6309 3.640e-06 ***
## ma3 0.855374 0.142139 6.0179 1.767e-09 ***
## intercept 4.570049 0.719035 6.3558 2.073e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre9)
 |
Visualisasi 22
|
jarque.bera.test(modpre9$residuals)
##
## Jarque Bera Test
##
## data: modpre9$residuals
## X-squared = 1498, df = 2, p-value < 2.2e-16
ks.test(modpre9$residuals, ecdf(modpre9$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre9$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre9)
 |
Visualisasi 23
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,3) with non-zero mean
## Q* = 2.0581, df = 3, p-value = 0.5604
##
## Model df: 5. Total lags used: 8
#Model ARIMA(2, 0, 3)
modpre10 <- arima(lpe_bc, order = c(2, 0, 3), include.mean = T)
coeftest(modpre10)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.019015 0.205595 0.0925 0.9263
## ar2 0.037517 0.184237 0.2036 0.8386
## ma1 1.092562 0.178878 6.1078 1.010e-09 ***
## ma2 0.935785 0.208368 4.4910 7.088e-06 ***
## ma3 0.843218 0.146820 5.7432 9.290e-09 ***
## intercept 4.569282 0.741197 6.1647 7.060e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modpre10)
 |
Visualisasi 24
|
jarque.bera.test(modpre10$residuals)
##
## Jarque Bera Test
##
## data: modpre10$residuals
## X-squared = 1500.1, df = 2, p-value < 2.2e-16
ks.test(modpre10$residuals, ecdf(modpre10$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modpre10$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modpre10)
 |
Visualisasi 24
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,3) with non-zero mean
## Q* = 2.0412, df = 3, p-value = 0.5639
##
## Model df: 6. Total lags used: 9
#Perbandingan AIC
knitr::kable(cbind(c("ARIMA (1, 0, 0)", "ARIMA (2, 0, 0)", "ARIMA (0, 0, 1)",
"ARIMA (0, 0, 2)", "ARIMA (1, 0, 1)", "ARIMA (1, 0, 2)", "ARIMA (2, 0, 2)",
"ARIMA (0, 0, 3)", "ARIMA (1, 0, 3)", "ARIMA (2, 0, 3)"),
c(modpre1$aic, modpre2$aic, modpre3$aic, modpre4$aic,
modpre5$aic, modpre6$aic, modpre7$aic, modpre8$aic, modpre9$aic, modpre10$aic),
c(accuracy(modpre1)[5], accuracy(modpre2)[5], accuracy(modpre3)[5], accuracy(modpre4)[5],
accuracy(modpre5)[5], accuracy(modpre6)[5], accuracy(modpre7)[5], accuracy(modpre8)[5],
accuracy(modpre9)[5], accuracy(modpre10)[5])),
col.names = c("Model", "Nilai AIC", "MAPE"))
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
## Warning in trainingaccuracy(object, test, d, D): test elements must be within
## sample
| Model |
Nilai AIC |
MAPE |
| ARIMA (1, 0, 0) |
158.645853576448 |
NaN |
| ARIMA (2, 0, 0) |
159.13975976907 |
NaN |
| ARIMA (0, 0, 1) |
163.315744573826 |
NaN |
| ARIMA (0, 0, 2) |
164.129032776725 |
NaN |
| ARIMA (1, 0, 1) |
159.584637257102 |
NaN |
| ARIMA (1, 0, 2) |
156.623374786337 |
NaN |
| ARIMA (2, 0, 2) |
155.600581219104 |
NaN |
| ARIMA (0, 0, 3) |
144.460666107613 |
NaN |
| ARIMA (1, 0, 3) |
146.452892844622 |
NaN |
| ARIMA (2, 0, 3) |
148.411825121056 |
NaN |
#Dengan pertimbangan signifikansi variabel, model terbaik adalah model ARIMA(0, 0, 3)
#Identifikasi Order Intervensi
preforcast <- forecast(modpre8, h = 8)
## Error in ts(x): 'ts' object must have one or more observations
par(mfrow=c(1, 1))
#menampilkan Forecast Pra Intervensi
preforcast
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 42 5.339365 3.786503 6.892227 2.96446779 7.714263
## 43 4.242082 1.954650 6.529513 0.74375714 7.740406
## 44 4.919852 2.190542 7.649162 0.74573287 9.093972
## 45 4.568962 1.545709 7.592215 -0.05470469 9.192628
## 46 4.568962 1.545709 7.592215 -0.05470469 9.192628
## 47 4.568962 1.545709 7.592215 -0.05470469 9.192628
## 48 4.568962 1.545709 7.592215 -0.05470469 9.192628
## 49 4.568962 1.545709 7.592215 -0.05470469 9.192628
#Identifikasi intervensi order dengan plot residual
error_iden <- rep(0, 33)
error_iden[1:33] <- preforcast$residuals[1:33]
error_iden[34:41] <- Y[34:41] - preforcast$mean
plot(error_iden, type = "h", xlab = "Waktu", ylab = "residuals", xaxt = "n", ylim = c(-15, 15))
abline(h=c(-3*sd(preforcast$residuals), 3*sd(preforcast$residuals)), col = "blue", lty = 2)
abline(h=0, col = "blue", lty = 2)
abline(v=34, col = "red", lty = 3, lwd = 1.5)
text(34, 12, "T", cex = 1, pos = 1)
axis(1, at = c(0, 10, 20, 30, 40), labels = c("T-34", "T-24", "T-14", "T-8", "T+6"))
 |
Visualisasi 25
|
#Estimasi parameter intervensi Pulse
library(TSA)
#Model ARIMA Intervensi order b = 0, r = 2, s = 0
modintervensi <- arimax(lpe_y, order=c(0, 0, 3),#Order ARIMA
seasonal=list(order=c(0, 0, 0), frequency=4), #Tanpa Musiman frekuensi bulanan (12)
include.mean=TRUE, #karena with non-zero mean jadi TRUE
xtransf= data.frame(Qtr2=1 * (seq(lpe_y) == 34)), #cek bahwa waktu Pra-Intervensi 1 sampai data ke-33, maka seq(cds) == 34 atau mulai data 34
transfer=list(c(2, 0)), method = "ML") # list(c(r, s))
modintervensi
##
## Call:
## arimax(x = lpe_y, order = c(0, 0, 3), seasonal = list(order = c(0, 0, 0), frequency = 4),
## include.mean = TRUE, method = "ML",
## xtransf = data.frame(Qtr2 = 1 * (seq(lpe_y) == 34)), transfer = list(c(2,
## 0)))
##
## Coefficients:
## ma1 ma2 ma3 intercept Qtr2-AR1 Qtr2-AR2 Qtr2-MA0
## -0.5085 0.3336 0.5997 5.240 0.9212 -0.2585 -10.4511
## s.e. 0.1444 0.1559 0.1273 0.166 0.0349 0.0427 0.5215
##
## sigma^2 estimated as 0.4767: log likelihood = -45.75, aic = 105.51
coeftest(modintervensi)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.508485 0.144378 -3.5219 0.0004284 ***
## ma2 0.333602 0.155866 2.1403 0.0323299 *
## ma3 0.599703 0.127298 4.7110 2.465e-06 ***
## intercept 5.240047 0.166034 31.5602 < 2.2e-16 ***
## Qtr2-AR1 0.921226 0.034889 26.4046 < 2.2e-16 ***
## Qtr2-AR2 -0.258476 0.042742 -6.0473 1.473e-09 ***
## Qtr2-MA0 -10.451127 0.521535 -20.0392 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(modintervensi)
 |
Visualisasi 26
|
jarque.bera.test(modintervensi$residuals)
##
## Jarque Bera Test
##
## data: modintervensi$residuals
## X-squared = 9.2253, df = 2, p-value = 0.009926
ks.test(modintervensi$residuals, ecdf(modintervensi$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: modintervensi$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(modintervensi)
 |
Visualisasi 27
|
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,3) with non-zero mean
## Q* = 5.4273, df = 3, p-value = 0.1431
##
## Model df: 7. Total lags used: 10
AIC(modintervensi)
## [1] 107.5082
#Diagnosa Intervensi Step
step34 <- filter(1*(seq(lpe_y) == 34), filter = c(0.921226, -0.258476), method = "rec", sides = 1) * -10.451127 #b = 0, r = 1, s = 0
step34
## Time Series:
## Start = 1
## End = 41
## Frequency = 1
## [1] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [6] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [11] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [16] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [21] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [26] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## [31] 0.00000000 0.00000000 0.00000000 -10.45112700 -9.62784992
## [36] -6.16806017 -3.19360926 -1.34774037 -0.41610212 -0.03496555
## [41] 0.07534124
plot(step34)
 |
Visualisasi 28
|
#Memodelkan ARIMA Intervensi Step34
#untuk menguji signifikansi intervensi pada model
arima_intervensi <- Arima(lpe_y, order = c(0, 0, 3),seasonal = c(0, 0, 0), xreg = step34,
include.mean = T)
summary(arima_intervensi)
## Series: lpe_y
## Regression with ARIMA(0,0,3) errors
##
## Coefficients:
## ma1 ma2 ma3 intercept xreg
## -0.5098 0.3336 0.6003 5.2408 1.0002
## s.e. 0.1388 0.1546 0.1265 0.1506 0.0297
##
## sigma^2 estimated as 0.5424: log likelihood=-45.75
## AIC=103.51 AICc=105.98 BIC=113.79
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02915366 0.690108 0.483772 4.494498 16.5069 0.2462113
## ACF1
## Training set 0.06937737
coeftest(arima_intervensi) #Intervensinya Signifikan
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.509806 0.138849 -3.6716 0.000241 ***
## ma2 0.333592 0.154589 2.1579 0.030934 *
## ma3 0.600306 0.126497 4.7456 2.079e-06 ***
## intercept 5.240765 0.150592 34.8011 < 2.2e-16 ***
## xreg 1.000170 0.029673 33.7069 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tsdiag(arima_intervensi)
 |
Visualisasi 29
|
jarque.bera.test(arima_intervensi$residuals)
##
## Jarque Bera Test
##
## data: arima_intervensi$residuals
## X-squared = 9.1631, df = 2, p-value = 0.01024
ks.test(arima_intervensi$residuals, ecdf(arima_intervensi$residuals))
##
## One-sample Kolmogorov-Smirnov test
##
## data: arima_intervensi$residuals
## D = 0.02439, p-value = 1
## alternative hypothesis: two-sided
checkresiduals(arima_intervensi)
 |
Visualisasi 30
|
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(0,0,3) errors
## Q* = 5.0604, df = 3, p-value = 0.1674
##
## Model df: 5. Total lags used: 8
AIC(arima_intervensi)
## [1] 103.508
#Peramalan hasil model Intervensi
xreg_rob = forecast(auto.arima(step34), h = 3)$mean
#Melakukan peramalan 3 triwulan ke depan Q2 2022 - Q4 2022
fc <- forecast(arima_intervensi, xreg = xreg_rob)
fc
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2022 Q2 5.506978 4.541276 6.472679 4.030064 6.983891
## 2022 Q3 5.293738 4.221357 6.366119 3.653673 6.933803
## 2022 Q4 4.886250 3.774480 5.998020 3.185944 6.586556
#Melihat Akurasi Peramalan
accuracy(arima_intervensi, xreg = xreg_rob)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.02915366 0.690108 0.483772 4.494498 16.5069 0.2462113
## ACF1
## Training set 0.06937737
#Melihat R Square Model ARIMA (0, 0, 3) Tanpa Intervensi
library(caret)
arimasaja <- Arima(lpe_bc, order = c(0, 0, 3), include.mean = T)
pred_arimasaja <- arimasaja$fitted
R2(pred_arimasaja, as.numeric(lpe_y), formula = "traditional")
## [1] 0.1421703
#Melihat R Square Model ARIMA dengan Intervensi
pred <- arima_intervensi$fitted
R2(pred, lpe_y, formula = "traditional")
## [1] 0.9283124
#Menampilkan Prediksi dan nilai asli
library(ggfortify)
prediksi <- as.numeric(pred)
trueval <- as.numeric(lpe_y)
predtrueval <- data.frame(prediksi, trueval)
predtrueval.ts <- ts(predtrueval)
theme_set(theme_bw())
autoplot(predtrueval.ts,xlab = "Periode", ylab = "Laju Pertumbuhan (%)") +
ggtitle("Plot Runtun Waktu Prediksi dan Data Asli") +
theme(plot.title = element_text(hjust = .5))
 |
Visualisasi 31
|
pesimis <- as.data.frame(fc$lower)[2]
optimis <- as.data.frame(fc$upper)[2]
a <- plot(lpe_y, col = "blue", type = "b", xlim = c(2011.1, 2023.1), ylim = c(-10, 13), pch = 19,
ylab = "Pertumbuhan (persen)", xlab = "Periode",
main = "Peramalan ARIMA(0, 0, 3) Intervensi Step-34 \n Laju Pertumbuhan Ekonomi Indonesia Q2 - Q4 2022")
b <- lines(pred, col = "red", type = "b")
c <- lines(fc$mean, col = "yellow", type = "b", pch = 19)
d <- lines(ts(pesimis$`95%`, start = 2022.1, frequency = 4), col = "red", type = "b")
e <- lines(ts(optimis$`95%`, start = 2022.1, frequency = 4), col = "green", type ="b")
text(2022.2,10, c("Skenario\n Optimis"), cex = 0.8)
text(2022.2,-3, c("Skenario\n Pesimis"), cex = 0.8)
abline(v =2020.2, col = "red", lty=2)
abline(v =2022, col = "red", lty=2)
text(2019.1, 0.2, c("Kebijakan PSBB"), cex = 0.8, col = "blue")
legend(2011.4,-7, legend = c("Data Aktual", "Prediksi dari ARIMA Intervensi", "Forecast", "Upper Forecast 95% - Optimis",
"Lower Forecast 95% - Pesimis"), col = c("blue", "red", "yellow", "green", "red"),
lty = 1, cex = 0.6, box.col = "white", horiz = F)
 |
Visualisasi 32
|
Berdasarkan model intervensi atau ARIMA Intervensi tentatif yang kita peroleh, dengan tingkat keyakinan 95 persen, pertumbuhan ekonomi Indonesia triwulan II 2022 diperkirakan sebesar 4,03 persen hingga level optimis sebesar 6,98 persen. Sedangkan skenario moderat, petumbuhan ekonomi Indonesia pada triwulan II 2022 diperkirakan sebesar 5,51 persen. Kemudian, dengan tingkat keyakinan 95 persen, pertumbuhan ekonomi Indonesia
triwulan III 2022 diperkirakan sebesar 3,65 persen hingga level optimis
sebesar 6,93 persen. Sedangkan skenario moderat, petumbuhan ekonomi
Indonesia pada triwulan III 2022 diperkirakan sebesar 5,29 persen. Sementara itu, dengan tingkat keyakinan 95 persen, pertumbuhan ekonomi Indonesia
triwulan IV 2022 diperkirakan sebesar 3,19 persen hingga level optimis
sebesar 6,59 persen. Sedangkan untuk skenario moderat, petumbuhan ekonomi
Indonesia pada triwulan IV 2022 diperkirakan sebesar 4,88 persen. Beberapa asumsi yang digunakan dalam peramalan pertumbuhan ekonomi ini adalah tidak adanya perubahan drastis pola konsumsi masyarakat, pandemi Covid-19 tetap terkendali, demikian halnya inflasi yang juga terkendali pada rentang batas normal.
Dari hasil prediksi tersebut, pertumbuhan ekonomi Indonesia pada tahun 2022 berada pada rentang 3,95 persen hingga 6,39 persen secara cummulative to cummulative (c to c). Adapun menurut skenario moderat, pertumbuhan ekonomi Indonesia tahun 2022 adalah sebesar 5,17 persen (c to c). Pertumbuhan skenario moderat dari model kita ini agaknya tidak jauh berbeda dengan prediksi dari ADB yang memprediksi pertumbuhan ekonomi Indonesia tahun 2022 sebesar 5,2 persen.
Dari hasil pemodelan ARIMA (0, 0, 3) dengan intervensi Step - 34 di atas, kita kemudian dapat menyusun persamaan model tentatif sebagai berikut:
 |
Model ARIMA Intervensi tentatif
|
Demikian sedikit ulasan mengenai ARIMA dengan intervensi fungsi Pulse dan Step menggunakan R. Selamat memahami dan mempraktikkan!
Sangat inspiratif mas, jadi pengen mencoba juga, forecasting dengan memperhatikan efek pandemi...
BalasHapus