Simulasi dan Pemodelan ARIMA (2, 2, 2)
Gracelinda Chrisageulis
2026-05-06
library(forecast)## Warning: package 'forecast' was built under R version 4.5.3library(tseries)## Warning: package 'tseries' was built under R version 4.5.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(TSA)## Warning: package 'TSA' was built under R version 4.5.3## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarset.seed(123)
n <- 200
phi1 <- 0.5
phi2 <- -0.7
teta1 <- -0.3
teta2 <- -0.5
ar <- numeric(4)
ar[1] <- phi1 + 2
ar[2] <- phi2 - 2*phi1 - 1
ar[3] <- phi1 - 2*phi2
ar[4] <- phi2
ma <- numeric(2)
ma[1] <- -teta1
ma[2] <- -teta2
et <- rnorm(n, mean = 0, sd = sqrt(1.05))
Yt <- numeric(n)
for (t in 5:n) {
Yt[t] <- ar[1] * Yt[t-1] +
ar[2] * Yt[t-2] +
ar[3] * Yt[t-3] +
ar[4] * Yt[t-4] +
et[t] +
ma[1] * et[t-1] +
ma[2] * et[t-2]
}
data <- ts(Yt[-(1:50)])
ts.plot(data, col = "pink", lwd = 1.5,
main = "Plot Data Yt - 150 Observasi (ARIMA 2,2,2)",
ylab = "Yt", xlab = "Waktu")
points(data, pch = 20, col = "pink", cex = 0.5)
# 2. Split Train Test
train.data <- data[1:120]
test.data <- data[121:150]Visualisasi data training digunakan untuk mencari model terbaik
ts.plot(train.data, col="pink")
title(main = "Plot data train")
points(train.data, pch = 20, col = "pink")
#3.Model ARIMA # Uji stasioneritas Mengidentifikasi kestasioneran data
dilakukan dengan uji ADF.
Uji Augmented Dickey-Fuller (ADF) adalah salah satu alat identifikasi untuk menguji kestasioneran data dan memiliki hipotesis sebagai berikut:
\(H_0\) : Data tidak stasioner (unit roots mendekati 1)
\(H_1\): Data stasioner (unit roots tidak mendekati 1)
adf.test(train.data, alternative = "stationary",k=trunc((length(train.data)-1)^(1/3)))##
## Augmented Dickey-Fuller Test
##
## data: train.data
## Dickey-Fuller = -2.3354, Lag order = 4, p-value = 0.4373
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF, diperoleh nilai p-value sebesar 0,4373 yang lebih besar dari 0,05. Hal ini menunjukkan bahwa H0 tidak dapat ditolak, sehingga data bersifat tidak stasioner. Oleh karena itu, perlu dilakukan proses differencing agar data menjadi stasioner sebelum dilakukan tahap identifikasi model selanjutnya.
Differencing
#Differencing Ordo 1
data.dif1 <- diff(train.data, difference=1)
plot.ts(data.dif1)
points(data.dif1)
adf.test(data.dif1, alternative = "stationary",k=trunc((length(train.data)-1)^(1/3)))##
## Augmented Dickey-Fuller Test
##
## data: data.dif1
## Dickey-Fuller = -1.7984, Lag order = 4, p-value = 0.6603
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF pada differencing pertama, diperoleh nilai p-value sebesar 0,6603 yang lebih besar dari 0,05. Hal ini menunjukkan bahwa H0 tidak dapat ditolak, sehingga data masih belum stasioner. Oleh karena itu, perlu dilakukan differencing tahap kedua agar data menjadi stasioner sebelum melanjutkan ke proses identifikasi model tentatif.
#Differencing Ordo 2
data.dif2 <- diff(train.data, difference=2)
plot.ts(data.dif2)
points(data.dif2)
adf.test(data.dif2, alternative = "stationary",k=trunc((length(train.data)-1)^(1/3)))## Warning in adf.test(data.dif2, alternative = "stationary", k =
## trunc((length(train.data) - : p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: data.dif2
## Dickey-Fuller = -4.2162, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF pada differencing kedua, diperoleh nilai p-value sebesar 0,01 yang lebih kecil dari 0,05. Hal ini menunjukkan bahwa H0 ditolak, sehingga dapat disimpulkan bahwa data telah stasioner setelah dilakukan differencing orde ke-2. Oleh karena itu, proses selanjutnya adalah melakukan identifikasi model tentatif.
#4. Identifikasi kandidat model Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=0\)
#Plot ACF
acf(data.dif2, lag.max = 20, col = "pink")
#Plot PACF
pacf(data.dif2, lag.max = 20, col = "pink")
#EACF
eacf(data.dif2)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x x o x o o o o o o o o
## 1 x x x x x x o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x o o o o o o o o o o o o oBerdasarkan hasil EACF di atas, kandidat model yang diperoleh adalah ARMA(2,2), karena pola segitiga dari simbol “o” mulai terbentuk secara jelas pada orde AR = 2 dan MA = 2 serta berlanjut secara konsisten ke arah kanan bawah. Oleh karena itu, model tersebut dapat dipertimbangkan sebagai model yang paling sesuai untuk merepresentasikan data.
#5. KANDIDAT MODEL
library(forecast)
# Model kandidat (berdasarkan EACF)
model_222 <- Arima(data, order = c(2,2,2))
# Ringkasan model
summary(model_222)## Series: data
## ARIMA(2,2,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 0.5305 -0.7244 0.1668 0.4470
## s.e. 0.1029 0.0778 0.1174 0.1184
##
## sigma^2 = 0.9761: log likelihood = -206.88
## AIC=423.75 AICc=424.18 BIC=438.74
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.03190303 0.9680213 0.7607377 0.007802163 0.2785288 0.1929199
## ACF1
## Training set 0.004254808Model yang diperoleh: \[ (1 - B)^2 Y_t = \phi_1 (1 - B)^2 Y_{t-1} + \phi_2 (1 - B)^2 Y_{t-2} + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} \\Y_t = 0.5305Y_{t-1} - 0.7244Y_{t-2} + \epsilon_t + 0.1668\epsilon_{t-1} + 0.4470\epsilon_{t-2} \] Model ARIMA(2,2,2) memiliki parameter yang stabil dan memberikan nilai AIC, AICc, serta BIC yang menunjukkan kecocokan model yang baik. Nilai error yang rendah menandakan akurasi model cukup tinggi, dan residual yang tidak berkorelasi (ACF1 mendekati nol) menunjukkan bahwa model sudah memenuhi asumsi white noise. Dengan demikian, model ini layak digunakan untuk analisis dan peramalan.
#6. PEMILIHAN MODEL TERBAIK (AIC, BIC, AICc)
# Karena hanya satu model, langsung tampilkan kriterianya
n <- length(data)
k <- length(model_222$coef)
AICc_manual <- AIC(model_222) + (2*k*(k+1))/(n-k-1)
kriteria <- data.frame(
Model = "ARIMA(2,2,2)",
AIC = AIC(model_222),
BIC = BIC(model_222),
AICc = AICc_manual
)
kriteria## Model AIC BIC AICc
## 1 ARIMA(2,2,2) 423.7546 438.7406 424.0304Model ARIMA(2,2,2) memiliki nilai AIC sebesar 423,7546, BIC sebesar 438,7406, dan AICc sebesar 424,0304. Nilai-nilai ini menunjukkan tingkat kecocokan model yang baik dalam menjelaskan data dengan mempertimbangkan kompleksitas model. Karena hanya satu model yang diuji, maka model ARIMA(2,2,2) dipilih sebagai model terbaik dan layak digunakan untuk analisis serta peramalan.
#7. DIAGNOSTIK MODEL (Normalitas Residual)
# Histogram residual
hist(residuals(model_222),
main = "Histogram Residual ARIMA(2,2,2)",
xlab = "Residual")
# QQ Plot
qqnorm(residuals(model_222))
qqline(residuals(model_222), col = "pink")
# Uji Shapiro-Wilk
shapiro.test(residuals(model_222))##
## Shapiro-Wilk normality test
##
## data: residuals(model_222)
## W = 0.98976, p-value = 0.3452Hasil uji Shapiro-Wilk menunjukkan nilai p-value sebesar 0,3452, yang lebih besar dari 0,05. Hal ini menunjukkan bahwa tidak terdapat cukup bukti untuk menolak hipotesis nol, sehingga residual dapat dianggap berdistribusi normal. Dengan demikian, asumsi normalitas pada model ARIMA(2,2,2) telah terpenuhi.