Melanjutkan pemahaman mengenai analisis deret waktu, pada dokumen ini akan dilakukan simulasi pembangkitan data menggunakan model ARIMA(1,1,1) dengan parameter baru. Tujuan dari latihan ini adalah untuk mereplikasi siklus penuh pemodelan Box-Jenkins: mulai dari identifikasi secara visual (ACF/PACF), pengujian formal kestasioneran, penentuan kandidat, hingga estimasi model terbaik. Di akhir tahapan, kita akan membandingkan apakah model yang teridentifikasi oleh data sampel (berdasarkan nilai AIC terkecil dan auto.arima) akan persis sama dengan model Data Generating Process (DGP) awal yang kita gunakan untuk membangkitkan data.
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.2
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(TSA)
## Warning: package 'TSA' was built under R version 4.5.3
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.2
## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.Arima TSA
## plot.Arima TSA
membangkitkan data time series dengan model ARIMA(1,1,1). Pada contoh sebelumnya digunakan AR = 0.7 dan MA = -0.5. Kali ini, kita akan menggunakan nilai acak yang berbeda, misalnya AR = 0.4 dan MA = 0.6.
set.seed(456)
# Panjang data
n <- 200
# Parameter baru ARIMA(p=1, d=1, q=1)
ar_baru <- 0.4 # AR(1)
ma_baru <- 0.6 # MA(1)
# Simulasi data
ts_simulasi <- arima.sim(model = list(order = c(1,1,1), ar = ar_baru, ma = ma_baru), n = n)
# Plot data deret waktu asli
ts.plot(ts_simulasi, main = "Simulasi Data Baru ARIMA(1,1,1)", ylab = "Nilai", xlab = "Waktu")
Penjelasan:
Terlihat dari plot di atas bahwa data memiliki tren dan rata-rata yang berubah seiring waktu, yang merupakan indikasi awal bahwa data tidak stasioner.
par(mfrow=c(1,2))
acf(ts_simulasi, main="ACF Data Asli")
pacf(ts_simulasi, main="PACF Data Asli")
Penjelasan:
Plot ACF menurun secara sangat lambat (tails off slowly), yang merupakan karakteristik utama dari data yang tidak stasioner dalam rata-rata.
adf.test(ts_simulasi)
##
## Augmented Dickey-Fuller Test
##
## data: ts_simulasi
## Dickey-Fuller = -2.3894, Lag order = 5, p-value = 0.413
## alternative hypothesis: stationary
Penjelasan:
(Asumsi output: p-value > 0.05). Karena nilai p-value lebih besar dari taraf signifikansi 0.05, kita tolak H0 dan disimpulkan bahwa data tidak stasioner. Harus dilakukan proses differencing.
diff_ts <- diff(ts_simulasi)
head(diff_ts)
## Time Series:
## Start = 2
## End = 7
## Frequency = 1
## [1] 1.7418301 0.1248623 0.8115561 2.1000072 3.0871674 0.7864189
par(mfrow=c(1,2))
acf(diff_ts, main="ACF Data Diff-1")
pacf(diff_ts, main="PACF Data Diff-1")
Penjelasan:
Setelah di-differencing, plot ACF dan PACF sudah meluruh dengan cepat masuk ke dalam batas signifikansi (garis putus-putus biru).
adf.test(diff_ts)
## Warning in adf.test(diff_ts): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff_ts
## Dickey-Fuller = -5.0485, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
Karena nilai p-value sudah lebih kecil dari 0.05 (biasanya bernilai 0.01 pada data simulasi d=1), maka data hasil differencing sudah stasioner.
data_ts <- ts(diff_ts)
head(data_ts)
## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] 1.7418301 0.1248623 0.8115561 2.1000072 3.0871674 0.7864189
eacf(data_ts)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o x x x
## 1 x x o x o o o o o o o o o o
## 2 x x x o o o o o o o o o o o
## 3 x x x x o o o o o o o o o o
## 4 x x x x o o o o o o o o o o
## 5 o x x x x o o o o o o o o o
## 6 o x x x o x o o o o o o o o
## 7 o x x o o o o o o o o o o o
Berdasarkan plot ACF dan PACF (yang mengekor/terpotong pada lag-lag awal) serta melihat pola segitiga “O” pada output EACF, kita dapat menyusun beberapa kandidat model, antara lain:
ARIMA(1,1,1) (Model asli)
ARIMA(0,1,1)
ARIMA(1,1,0)
ARIMA(1,1,2)
Kita gunakan fungsi auto.arima dari package forecast untuk melihat rekomendasi sistem algoritma.
# Menggunakan data asli karena auto.arima akan otomatis mendeteksi d=1
model_auto <- auto.arima(ts_simulasi, trace = TRUE)
##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : 567.7085
## ARIMA(0,1,0) with drift : 731.2399
## ARIMA(1,1,0) with drift : 595.5844
## ARIMA(0,1,1) with drift : 597.0389
## ARIMA(0,1,0) : 729.2141
## ARIMA(1,1,2) with drift : 567.468
## ARIMA(0,1,2) with drift : 578.6167
## ARIMA(1,1,1) with drift : 565.3696
## ARIMA(2,1,1) with drift : 566.673
## ARIMA(2,1,0) with drift : 571.6485
## ARIMA(1,1,1) : 563.2927
## ARIMA(0,1,1) : 595.0887
## ARIMA(1,1,0) : 593.5284
## ARIMA(2,1,1) : 564.5705
## ARIMA(1,1,2) : 565.3698
## ARIMA(0,1,2) : 576.6733
## ARIMA(2,1,0) : 569.5711
## ARIMA(2,1,2) : 565.5826
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(1,1,1) : 566.4682
##
## Best model: ARIMA(1,1,1)
summary(model_auto)
## Series: ts_simulasi
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.4593 0.5656
## s.e. 0.0801 0.0804
##
## sigma^2 = 0.969: log likelihood = -280.17
## AIC=566.35 AICc=566.47 BIC=576.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.003235597 0.9769904 0.7961331 0.09290334 8.508666 0.6515119
## ACF1
## Training set -0.01395989
Untuk membandingkan nilai AIC (Akaike Information Criterion) secara akurat, kita harus melakukan fitting model menggunakan data asli (ts_simulasi) dengan parameter differencing d=1 di dalam argumen order=c(p,1,q).
cat("\nAIC Kandidat 1: ARIMA(1,1,1) =", AIC(arima(ts_simulasi, order=c(1,1,1), method="ML")))
##
## AIC Kandidat 1: ARIMA(1,1,1) = 566.3457
cat("\nAIC Kandidat 2: ARIMA(0,1,1) =", AIC(arima(ts_simulasi, order=c(0,1,1), method="ML")))
##
## AIC Kandidat 2: ARIMA(0,1,1) = 591.4524
cat("\nAIC Kandidat 3: ARIMA(1,1,0) =", AIC(arima(ts_simulasi, order=c(1,1,0), method="ML")))
##
## AIC Kandidat 3: ARIMA(1,1,0) = 595.3207
cat("\nAIC Kandidat 4: ARIMA(1,1,2) =", AIC(arima(ts_simulasi, order=c(1,1,2), method="ML")))
##
## AIC Kandidat 4: ARIMA(1,1,2) = 568.2587
Model terbaik adalah model yang memiliki nilai AIC terkecil.
Apakah model hasil pembangkitan sama dengan model terbaik dari proses pemodelan?
Seringkali dalam praktiknya, model terbaik yang direkomendasikan oleh ACF/PACF, EACF, maupun auto.arima BISA BERBEDA dengan model asli (ARIMA(1,1,1)) yang kita gunakan untuk membangkitkan data. Misalnya, algoritma mungkin memilih ARIMA(0,1,1) atau ARIMA(1,1,2) sebagai model dengan AIC terkecil meskipun data murni dibangkitkan dari (1,1,1).
Mengapa hal ini bisa terjadi?
Prinsip Parsimoni (Kriteria Penalti AIC), Nilai AIC dirancang untuk memberikan penalti pada model yang kompleks (memiliki banyak parameter). Jika koefisien AR(0.4) dan MA(0.6) menghasilkan pola yang bisa diaproksimasi dengan sangat baik oleh model yang lebih sederhana (seperti ARIMA(0,1,1)), maka rumus AIC akan memenangkan model yang lebih sederhana tersebut karena dianggap lebih efisien (parsimoni).
Pengaruh White Noise dan Variabilitas Sampel, Data time series dibangkitkan dari persamaan matematis yang ditambahkan dengan unsur error (acak/white noise). Pada sampel berukuran 200, variasi acak dari error ini terkadang mengaburkan “sinyal” asli dari parameter AR atau MA, sehingga pola yang tertangkap di sampel bergeser sedikit dari struktur aslinya.
Optimasi Maximum Likelihood (ML), Proses estimasi parameter di software (seperti fungsi arima) menggunakan iterasi numerik untuk mencari nilai Maximum Likelihood. Pada ruang sampel terbatas, kombinasi parameter lain kadang menghasilkan keakuratan (likelihood) yang secara kebetulan sedikit lebih tinggi pada data sampel spesifik tersebut dibandingkan parameter aslinya.