Libraries and Setup
Di era ini, kita sering membutuhkan analisis prediktif untuk membantu kita dalam membuat keputusan. Salah satu hal penting dalam prediksi adalah memprediksi untuk data-data di masa depan. Jenis prediksi ini sering juga disebut dengan peramalan.
Peramalan banyak dibutuhkan di berbagai situasi: menentukan apakah akan membangun pembangkit listrik lain dalam lima tahun ke depan membutuhkan prakiraan permintaan di masa depan; staf penjadwalan di pusat panggilan minggu depan membutuhkan prakiraan volume panggilan; persediaan persediaan membutuhkan prakiraan kebutuhan persediaan1.
Tujuan dibuatnya artikel ini adalah untuk memperkenalkan salah satu metode peramalan dengan melibatkan variabel prediktor, yaitu ARIMAX
. Secara khusus, artikel ini bertujuan untuk:
– Memperkenalkan peramalan yang melibatkan prediktor
– Memperkenalkan dan aplikasi dari ARIMAX
– Membandingkan hasil peramalan ARIMA dengan ARIMAX
# Import library
library(fpp3)
library(forecast)
library(lmtest)
library(padr)
library(tseries)
Tentang ARIMA
Auto Regressive Integrated Moving Average (ARIMA)(p,d,q) merupakan versi lanjutan dari model Auto Regressive (AR), Moving Average (MA), dan Auto Regressive Moving Average (ARMA)2. Model ARIMA merupakan model yang diaplikasikan pada permasalahan deret waktu/time series. ARIMA menggabungkan tiga jenis pemodelan ke dalam satu model3:
- I: Differencing dilambangkan dengan \(d\). I memberi tahu kita jumlah seri berbeda yang diubah antara pengamatan berturut-turut terhadap seri aslinya.
- AR: Auto Regressive dilambangkan dengan \(p\). AR memberi tahu kita orde
dari lags yang diperlukan untuk menyesuaikan proses AR dengan seri stasioner. ACF dan PACF membantu kami mengidentifikasi parameter terbaik untuk proses AR. - MA: Moving Average dilambangkan dengan \(q\). MA memberitahu kita
jumlah error terms dalam rangkaian yang akan diregresikan untuk mengurangi
perbedaan error proses AR ke white noise.
Tentang ARIMAX
ARIMAX atau Regresi ARIMA merupakan perpanjangan dari model ARIMA. Dalam peramalan, metode ini juga melibatkan variabel independen 4. Model ARIMAX merepresentasikan komposisi rangkaian waktu keluaran menjadi komponen-komponen berikut: autoregressive (AR), moving average (MA), terintegrasi (I), dan facktor eksternal (X) 5. Faktor eksternal (X) mencerminkan penggabungan tambahan dari nilai sekarang \(u_i(t)\) dan nilai masa lalu \(u_i(t-j)\) dari input faktor eksternal (variabel independen) ke dalam model ARIMAX6.
Rumus Multiple linear regression models:
\[Y = \beta_0 + \beta_1*x_1+…+\beta_i*x_i+\varepsilon\]
Di mana \(Y\) merupakan sebuah variabel dependen dari variabel prediktor \(x_i\) dan \(\varepsilon\) biasanya diasumsikan sebagai error/white noise. Kami akan mengganti \(\varepsilon\) dengan \(n_t\) pada persamaan. Error \(\phi_t\) diasumsikan mengikuti hasil dari model ARIMA. Sebagai contoh, jika \(n_t\) mengikuti model ARIMA (1,1,1), dapat kita tuliskan
\[Y = \beta_0 + \beta_1x_1+\beta_2x_2+…+\beta_ix_i+\eta_t\]
\[(1-\phi_1B)(1-B)\eta_t = (1+\phi_1B)\varepsilon_t\]
Di mana \(\varepsilon_t\), merupakan seri white noise. Model ARIMAX memiliki two error terms; the error dari model regresi yang dinotasikan dengan \(\phi_t\) dan error dari model ARIMA model yang dinotasikan dengan \(\varepsilon_t\).
Studi Kasus: Peramalan Konsumsi berdasarkan Tingkat Pendapatan, Produksi, Pengangguran, dan Dana Cadangan
Pada kasus ini, akan diramalkan persentase perubahan ekonomi di USA menggunakan data us_change
dari library fpp3
.
us_change
#> # A tsibble: 198 x 6 [1Q]
#> Quarter Consumption Income Production Savings Unemployment
#> <qtr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1970 Q1 0.619 1.04 -2.45 5.30 0.9
#> 2 1970 Q2 0.452 1.23 -0.551 7.79 0.5
#> 3 1970 Q3 0.873 1.59 -0.359 7.40 0.5
#> 4 1970 Q4 -0.272 -0.240 -2.19 1.17 0.700
#> 5 1971 Q1 1.90 1.98 1.91 3.54 -0.100
#> 6 1971 Q2 0.915 1.45 0.902 5.87 -0.100
#> 7 1971 Q3 0.794 0.521 0.308 -0.406 0.100
#> 8 1971 Q4 1.65 1.16 2.29 -1.49 0
#> 9 1972 Q1 1.31 0.457 4.15 -4.29 -0.200
#> 10 1972 Q2 1.89 1.03 1.89 -4.69 -0.100
#> # … with 188 more rows
Data di atas merupakan data sosial-ekonomi di United States pada quarter pertama tahun 1970 sampai quarter kedua tahun 2019 yang terdiri dari:
Quarter
: Quarter dan tahunConsumption
: Tingkat konsumsiIncome
: Tingkat pendapatanProduction
: Tingkat ProduksiSavings
: Dana cadanganUnemployment
: Tingkat pengangguran
Eksploratory Data Analysis (EDA)
Sebelum melakukan pemodelan, dilakukan EDA terlebih dahulu dengan membuat line plot dari setiap variabel baik dependen maupun independen untuk mengetahui pola dari setiap variabel apakah sudah stasioner atau belum.
us_change %>%
pivot_longer(-Quarter, names_to = "variable", values_to = "value") %>%
ggplot(aes(x = Quarter, y = value)) +
geom_line() +
facet_grid(variable ~ ., scales = "free_y") +
labs(title = "Perubahan Sosio-Ekonomi USA dari waktu ke waktu") +
theme_minimal()
Dari plot di atas, terlihat sepertinya kelima variabel (Consumption, Income, Production, Savings, dan Unemployment) sudah stationer. Namun, akan tetap dilakukan pengujian secara statistik untuk mengecek stationarity data supaya hasil yang diperoleh dapat bersifat objektif. Dalam hal ini, dapat menggunakan ADF test atau KPSS test (lebih baik jika mencoba keduanya).
ADF test
Untuk mendapatkan hasil yang pasti dan objektif, kita bisa melakukan Augmented Dickey-Fuller (ADF) test dengan menggunakan fungsi adf.test()
dari library tseries
.
H0: Punya unit root (tidak stationer)
H1: Tidak punya unit root (stationer)
p-value < 0.05 (alpha), data stationer
KPSS test
Disarankan pula melakukan uji lainnya (KPSS test) untuk mendapatkan kesimpulan yang relatif konstan dan pasti berdasarkan data historis dengan menggunakan fungsi kpss.test()
dari library tseries
H0: rata-rata dan variansi konstan (Data stationer)
H1: rata-rata dan variansi tidak konstan (Data tidak stationer)
df <- us_change[, -1]
stationary_test <- data.frame("ADF" = double(), "KPSS" = double())
for (i in 1:ncol(df)) {
stationary_test[i, "ADF"] <- adf.test(pull(df[, i]))$p.value
stationary_test[i, "KPSS"] <- kpss.test(pull(df[, i]))$p.value
}
stationary_test %>%
mutate(variable = colnames(df)) %>%
select(variable, ADF, KPSS)
#> variable ADF KPSS
#> 1 Consumption 0.01 0.1
#> 2 Income 0.01 0.1
#> 3 Production 0.01 0.1
#> 4 Savings 0.01 0.1
#> 5 Unemployment 0.01 0.1
Berdasarkan p-value pada ADF test, seluruh variabel mempunyai p-value 0.01 (Stasioner) dan KPSS test 0.1 (Stasioner), sehingga didapatkan kesimpulan bahwa kelima variabel tersebut signifikan stasioner.
Cross Validation
Data us_change
akan dibagi menjadi 2 subset data, yaitu sebanyak 4 tahun (2016 – 2019) untuk data test
dan 35 tahun (1970 – 2015) untuk data train
test <- us_change %>%
mutate(year = year(Quarter)) %>%
filter(year >= 2016)
train <- us_change %>%
mutate(year = year(Quarter)) %>%
filter(year < 2016)
Model Fitting dengan ARIMA
Kita akan mecoba melakukan fitting model menggunakan model ARIMA terlebih dahulu
fit_arima <- auto.arima(ts(train$Consumption, frequency = 4), seasonal = F)
summary(fit_arima)
#> Series: ts(train$Consumption, frequency = 4)
#> ARIMA(1,0,3) with non-zero mean
#>
#> Coefficients:
#> ar1 ma1 ma2 ma3 mean
#> 0.5747 -0.3581 0.093 0.1946 0.7418
#> s.e. 0.1526 0.1635 0.081 0.0857 0.0936
#>
#> sigma^2 estimated as 0.3533: log likelihood=-163.03
#> AIC=338.06 AICc=338.53 BIC=357.35
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.001107612 0.5862617 0.4378066 -35.61455 161.7278 0.672272
#> ACF1
#> Training set -0.002685885
Dari output di atas diperoleh model ARIMA(1,0,3) dengan nilai RMSE pada data training sebesar 0.58
Kita akan mencoba melakukan forecasting terhadap data test, kemudian mengkalkulasi error dari kedual model (ARIMA dan ARIMAX)
prediction_arima <- forecast(object = fit_arima, h = nrow(test))
Sebelum mengkalkulasi error yang diperoleh terlebih dahulu kita akan melihat visualisasi hasil forecast dari kedua model
prediction_arima %>%
autoplot() +
theme_minimal()
Model Fitting dengan ARIMAX
Kemudian kita akan mencoba melakukan fitting model ARIMAX dengan variabel dependen adalah tingkat konsumsi dan variabel independen adalah tingkat pendapatan, tingkat produksi, dana cadangan, dan tingkat pengangguran. Lalu kita akan mencoba membandingkan hasil dari model ARIMA dan ARIMAX.
Adapun apabila ingin forecasting, namun tidak memiliki nilai predictor di masa depan, bisa dilakukan forecasting terlebih dahulu ke prediktornya, kemudian dilakukan forecast terhadap variabel targetnya.
fit_arimax <- train %>%
model(regarima = ARIMA(Consumption ~ Income + Production + Savings + Unemployment))
report(fit_arimax)
#> Series: Consumption
#> Model: LM w/ ARIMA(0,1,2) errors
#>
#> Coefficients:
#> ma1 ma2 Income Production Savings Unemployment
#> -1.0853 0.1087 0.7446 0.0384 -0.0527 -0.2095
#> s.e. 0.0717 0.0698 0.0419 0.0244 0.0031 0.1049
#>
#> sigma^2 estimated as 0.1023: log likelihood=-49.6
#> AIC=113.2 AICc=113.84 BIC=135.66
Dari output di atas diperoleh model ARIMAX(0,1,2)
prediction_arimax <- forecast(object = fit_arimax, new_data = test)
prediction_arimax %>%
autoplot(train) +
theme_minimal()
Error
Berdasarkan kedua plot di atas model ARIMAX lebih bisa memprediksi pola Consumption dibandingkan model ARIMA. Untuk membuktikan hal tersebut kita akan menghutung error dari kedua model tersebut, kemudian membandingkan hasilnya
print(paste("RMSE model ARIMA:", round(accuracy(object = prediction_arima, data = test)[2], 2)))
#> [1] "RMSE model ARIMA: 0.59"
print(paste("RMSE model ARIMAX:", round(forecast::accuracy(object = prediction_arimax, data = us_change)$RMSE, 2)))
#> [1] "RMSE model ARIMAX: 0.12"
Dari output di atas diketahui bahwa model ARIMAX menghasilkan error yang lebih kecil dibandingkan model ARIMA. Sehingga, final model yang akan digunakan adalah model ARIMAX. Model ARIMAX tersebut harus memenuhi beberapa asumsi supaya hasil yang peramalan data di masa depan bersifat BLUE (best, linier, unbiased, estimation)
Asumsi
Dalam pemodelan Time Series, terdapat 2 asumsi yang harus terpenuhi, yakni normality of residual dan no autocorrelation. Pada
gg_tsresiduals(fit_arimax)
gg_tsresiduals(fit_arimax)
Kesimpulan
ARIMAX model adalah metode yang dapat dijadikan solusi dalam time series forecasting yang melibatkan exogenous factor. Hal ini karena tidak selalu suatu variabel time series dapat dilakukan peramalan hanya berdasarkan informasi variabel itu sendiri di masa lalu, ada kemungkinan variabel tersebut juga sangat berkorelasi erat dengan faktor-faktor eksternal, seperti data-data dalam case sosial-ekonomi. Namun, hal yang perlu digaris bawahi adalah model ARIMAX cukup sulit untuk diinterpretasikan tidak seperti halnya model regresi linier karena estimate koefisien yang dihasilkan juga bergantung pada lag dari target variabel (pola variabel target di masa lalu).
References
- Forecasting: Principles and Practice↩︎
- Epidemiology and ARIMA model of positive-rate of influenza viruses among children in Wuhan, China: A nine-year retrospective study↩︎
- Machine Learning using R↩︎
- Comparison of Prediction Accuracy of Multiple Linear
Regression, ARIMA and ARIMAX Model for Pest Incidence of Cotton with Weather Factors↩︎ - Container Throughput Forecasting Using Dynamic Factor Analysis and ARIMAX Model↩︎
- Forecasting: Principles and Practice↩︎