Model ARIMA (2,2,2)
Model ARIMA(2,2,2), maka \(Y_t\) diperoleh dari penjabaran operator backshift sehingga untuk model ARIMA(2,2,2): \[ \begin{aligned} \phi_p(B)(1 - B)^d Y_t &= \mu + \theta_q(B) e_t \\ (1 - \phi_1 B - \phi_2 B^2)(1 - B)^2 Y_t &= \mu + (1 + \theta_1 B + \theta_2 B^2) e_t \\ (1 - \phi_1 B - \phi_2 B^2)(1 - 2B + B^2) Y_t &= \mu + (1 + \theta_1 B + \theta_2 B^2) e_t \\ (1 - 2B + B^2 - \phi_1 B + 2\phi_1 B^2 - \phi_1 B^3 - \phi_2 B^2 + 2\phi_2 B^3 - \phi_2 B^4) Y_t &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} \\ (1 - (2+\phi_1)B + (1+2\phi_1-\phi_2)B^2 - (\phi_1-2\phi_2)B^3 - \phi_2 B^4) Y_t &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} \\ Y_t - (2+\phi_1)Y_{t-1} + (1+2\phi_1-\phi_2)Y_{t-2} - (\phi_1-2\phi_2)Y_{t-3} - \phi_2 Y_{t-4} &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} \\ Y_t &= \mu + (2+\phi_1)Y_{t-1} - (1+2\phi_1-\phi_2)Y_{t-2} + (\phi_1-2\phi_2)Y_{t-3} + \phi_2 Y_{t-4} + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} \end{aligned} \]
Substitusi Parameter Ke Model ARIMA(2,2,2)
# Inisialisasi Parameter
phi1 <- -0.5
phi2 <- -0.7
theta1 <- -0.3
theta2 <- -0.5
mu <- 0
set.seed(123)
n_generate <- 200
sd_e <- sqrt(1.02)
# Generate error term (e_t)
e <- rnorm(n_generate, mean = 0, sd = sd_e)
# Inisialisasi vektor Y
Y <- numeric(n_generate)
# Memberikan nilai awal 0 untuk t=1 sampai t=4 karena lag terdalam adalah t-4
Y[1:4] <- 0
# 1. Membuat pengali (koefisien) AR dan MA secara terpisah di luar loop
coef_ar1 <- 2 + phi1
coef_ar2 <- -(1 + 2*phi1 - phi2)
coef_ar3 <- phi1 - 2*phi2
coef_ar4 <- phi2
coef_ma1 <- theta1
coef_ma2 <- theta2
# 2. Looping eksekusi persamaan Y_t
for (t in 5:n_generate) {
# Mendefinisikan suku AR secara terpisah
ar1_term <- coef_ar1 * Y[t-1]
ar2_term <- coef_ar2 * Y[t-2]
ar3_term <- coef_ar3 * Y[t-3]
ar4_term <- coef_ar4 * Y[t-4]
# Mendefinisikan suku MA secara terpisah
ma1_term <- coef_ma1 * e[t-1]
ma2_term <- coef_ma2 * e[t-2]
# Menyatukan semuanya menjadi Y_t sesuai rumus
Y[t] <- mu + ar1_term + ar2_term + ar3_term + ar4_term + e[t] + ma1_term + ma2_term
}
data_arima222 <- YGENERATE 200 DATA LALU HAPUS 50 DATA PERTAMA
data_bersih <- data_arima222[51:200]VISUALISASI DATA MENJADI TIME-SERIES OBJECT
data_ts <- ts(data_bersih)
plot(data_ts, main="Plot Data Time-Series ARIMA(2,2,2)",
ylab="Nilai", xlab="Waktu", col="blue", lwd=2)
SPLIT DATA (150 DATA)
80% Train = 120 data, 20% Test = 30 data
train_data <- window(data_ts, end=120)
test_data <- window(data_ts, start=121)VISUALISASI DATA TRAIN MENJADI TIME-SERIES OBJECT
plot(train_data, main="Train Data (120 Data)",
ylab="Nilai", xlab="Waktu", col="darkgreen", lwd=2)
Cek kestasioneran data
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:
\(𝐻_0:\)Data tidak stasioner (unit roots mendekati 1)
\(𝐻_1:\)Data stasioner (unit roots tidak mendekati 1)
cat("\n--- Uji ADF Data Train ---\n")##
## --- Uji ADF Data Train ---print(adf.test(train_data)) ##
## Augmented Dickey-Fuller Test
##
## data: train_data
## Dickey-Fuller = -0.42348, Lag order = 4, p-value = 0.9838
## alternative hypothesis: stationarytrain_diff1 <- diff(train_data, differences = 1)
cat("\n--- Uji ADF Differencing Orde 1 ---\n")##
## --- Uji ADF Differencing Orde 1 ---print(adf.test(train_diff1))## Warning in adf.test(train_diff1): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: train_diff1
## Dickey-Fuller = -4.0967, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationarytrain_diff2 <- diff(train_data, differences = 2)
cat("\n--- Uji ADF Differencing Orde 2 ---\n")##
## --- Uji ADF Differencing Orde 2 ---print(adf.test(train_diff2))## Warning in adf.test(train_diff2): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: train_diff2
## Dickey-Fuller = -7.4423, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5\%\) di atas dapat diketahui bersama bahwa \(p-value = 0.01 < \alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_0\) artinya data yang digunakan merupakan data yang stasioner.
Identifikasi kandidat model
Identifikasi kandidat model diperoleh berdasarkan nilai p dan q dimana nilai d=2 ACF, PACF, dan EACF (Data yang sudah stasioner)
acf(train_diff2, main="ACF Data Train", lag.max=20)
Berdasarkan plot ACF terlihat bahwa plot cuts off setelah lag ke-7
sehingga kandidat model yang diperoleh adalah ARIMA(0,2,7)
pacf(train_diff2, main="PACF Data Train", lag.max=20)
Berdasarkan plot PACF di atas kandidat model yang diperoleh adalah
ARIMA(6,2,0)
cat("\n--- Extended ACF (EACF) ---\n")##
## --- Extended ACF (EACF) ---eacf(train_diff2)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o x o x o o o o o o o
## 1 x x x x x o x o x o o o o o
## 2 x x o x o o o o o o o o o o
## 3 x x o x o o o o o o o o o o
## 4 o x o x o o o o o o o o o o
## 5 o x x x o o o o o o o o o o
## 6 x x x x o o o o o o o o o o
## 7 x x o x o x o o o o o o o oBerdasarkan hasil EACF di atas kandidat model yang diperoleh adalah ARIMA(2,2,4),ARIMA(3,2,4), dan ARIMA(4,2,4)
Kandidat Model
Diperoleh beberapa kandidat model yaitu: ARIMA(0,2,7),ARIMA(6,2,0),ARIMA(2,2,4),ARIMA(3,2,4), dan ARIMA(4,2,4).
# ARIMA(0,2,7)
arima027 <- arima(train_diff2, order = c(0,2,7), include.mean = TRUE, method = "ML")
arima027##
## Call:
## arima(x = train_diff2, order = c(0, 2, 7), include.mean = TRUE, method = "ML")
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ma1 ma2 ma3 ma4 ma5 ma6 ma7
## -2.6663 1.371 2.1484 -2.5871 0.2809 0.7559 -0.3028
## s.e. NaN 0.092 0.3751 0.4658 0.2782 0.0765 NaN
##
## sigma^2 estimated as 1.146: log likelihood = -184.7, aic = 383.39# ARIMA(6,2,0)
arima620 <- arima(train_diff2, order = c(6,2,0), include.mean = TRUE, method = "ML")## Warning in log(s2): NaNs produced## Warning in log(s2): NaNs producedarima620##
## Call:
## arima(x = train_diff2, order = c(6, 2, 0), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6
## -1.8300 -2.8139 -2.5231 -1.9137 -0.8698 -0.3466
## s.e. 0.0876 0.1732 0.2594 0.2598 0.1749 0.0887
##
## sigma^2 estimated as 2.435: log likelihood = -219.83, aic = 451.65# ARIMA(2,2,4)
arima224 <- arima(train_diff2, order = c(2,2,4), include.mean = TRUE, method = "ML")
arima224##
## Call:
## arima(x = train_diff2, order = c(2, 2, 4), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 ma4
## -0.6341 -0.7193 -2.1035 0.6523 1.0308 -0.5795
## s.e. 0.0761 0.0683 0.1252 0.3199 0.3184 0.1176
##
## sigma^2 estimated as 1.007: log likelihood = -176.49, aic = 364.99# ARIMA(3,2,4)
arima324 <- arima(train_diff2, order = c(3,2,4), include.mean = TRUE, method = "ML")
arima324##
## Call:
## arima(x = train_diff2, order = c(3, 2, 4), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 ma4
## -0.9939 -0.9269 -0.3348 -1.7826 -0.0906 1.5431 -0.6694
## s.e. 0.1339 0.1210 0.1275 0.1156 0.1711 0.1400 0.0915
##
## sigma^2 estimated as 0.9664: log likelihood = -174.91, aic = 363.82# ARIMA(4,2,4)
arima424 <- arima(train_diff2, order = c(4,2,4), include.mean = TRUE, method = "ML")
arima424##
## Call:
## arima(x = train_diff2, order = c(4, 2, 4), include.mean = TRUE, method = "ML")
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## -0.9242 -0.7895 -0.1987 0.1180 -1.8726 -0.0140 1.6495 -0.7627
## s.e. 0.1388 0.1968 0.1925 0.1332 NaN 0.1034 0.1039 NaN
##
## sigma^2 estimated as 0.9306: log likelihood = -174.26, aic = 364.52AICKandidatModel <- c(383.39, 451.65, 364.99, 363.82, 364.52)
KndidatModelARIMA <- c("ARIMA(0,2,7)","ARIMA(6,2,0)","ARIMA(2,2,4)","ARIMA(3,2,4)","ARIMA(4,2,4)")
compmodelARIMA <- cbind(KndidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA## Kandidat Model Nilai AIC
## 1 ARIMA(0,2,7) 383.39
## 2 ARIMA(6,2,0) 451.65
## 3 ARIMA(2,2,4) 364.99
## 4 ARIMA(3,2,4) 363.82
## 5 ARIMA(4,2,4) 364.52Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu ARIMA(3,2,4)
Model ARIMA (3,2,4)
Model ARIMA(3,2,4), maka \(Y_t\) diperoleh dari penjabaran operator backshift sehingga untuk model ARIMA(3,2,4): \[ \begin{aligned} \phi_p(B)(1 - B)^d Y_t &= \mu + \theta_q(B) e_t \\ (1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)(1 - B)^2 Y_t &= \mu + (1 + \theta_1 B + \theta_2 B^2 + \theta_3 B^3 + \theta_4 B^4) e_t \\ (1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)(1 - 2B + B^2) Y_t &= \mu + (1 + \theta_1 B + \theta_2 B^2 + \theta_3 B^3 + \theta_4 B^4) e_t \\ (1 - 2B + B^2 - \phi_1 B + 2\phi_1 B^2 - \phi_1 B^3 - \phi_2 B^2 + 2\phi_2 B^3 - \phi_2 B^4 - \phi_3 B^3 + 2\phi_3 B^4 - \phi_3 B^5) Y_t &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \theta_3 e_{t-3} + \theta_4 e_{t-4} \\ (1 - (2+\phi_1)B + (1+2\phi_1-\phi_2)B^2 - (\phi_1-2\phi_2+\phi_3)B^3 - (\phi_2-2\phi_3)B^4 - \phi_3 B^5) Y_t &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \theta_3 e_{t-3} + \theta_4 e_{t-4} \\ Y_t - (2+\phi_1)Y_{t-1} + (1+2\phi_1-\phi_2)Y_{t-2} - (\phi_1-2\phi_2+\phi_3)Y_{t-3} - (\phi_2-2\phi_3)Y_{t-4} - \phi_3 Y_{t-5} &= \mu + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \theta_3 e_{t-3} + \theta_4 e_{t-4} \\ Y_t &= \mu + (2+\phi_1)Y_{t-1} - (1+2\phi_1-\phi_2)Y_{t-2} + (\phi_1-2\phi_2+\phi_3)Y_{t-3} + (\phi_2-2\phi_3)Y_{t-4} + \phi_3 Y_{t-5} + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \theta_3 e_{t-3} + \theta_4 e_{t-4} \end{aligned} \] Dengan mensubstitusi nilai koefisien model (\(\phi_1 = -0.9939\), \(\phi_2 = -0.9269\), \(\phi_3 = -0.3348\), \(\theta_1 = -1.7826\), \(\theta_2 = -0.0906\), \(\theta_3 = 1.5431\), \(\theta_4 = -0.6694\)) dan \(\mu = 0\) ke dalam persamaan terakhir ARIMA(3,2,4) diperoleh \[ \begin{aligned} Y_t &= (2 + (-0.9939))Y_{t-1} - (1 + 2(-0.9939) - (-0.9269))Y_{t-2} \\ &\quad + ((-0.9939) - 2(-0.9269) + (-0.3348))Y_{t-3} + ((-0.9269) - 2(-0.3348))Y_{t-4} \\ &\quad + (-0.3348)Y_{t-5} + e_t + (-1.7826)e_{t-1} + (-0.0906)e_{t-2} \\ &\quad + (1.5431)e_{t-3} + (-0.6694)e_{t-4} \\ Y_t &= 1.0061 Y_{t-1} + 0.0609 Y_{t-2} + 0.5251 Y_{t-3} - 0.2573 Y_{t-4} - 0.3348 Y_{t-5} \\ &\quad + e_t - 1.7826 e_{t-1} - 0.0906 e_{t-2} + 1.5431 e_{t-3} - 0.6694 e_{t-4} \end{aligned} \]
Diagnostic model
ARIMA100diag <- stats::arima(train_diff2, order = c(3,2,4), method = "ML")
checkresiduals(ARIMA100diag)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,2,4)
## Q* = 6.0717, df = 3, p-value = 0.1082
##
## Model df: 7. Total lags used: 10Berdasarkan visualisasi checkresiduals, plot sisaan tampak berfluktuasi secara stasioner di sekitar angka nol tanpa membentuk pola tertentu. Pada plot ACF, tidak terlihat adanya nilai korelasi lag yang signifikan (keluar dari batas garis putus-putus biru). Selain itu, histogram sisaan menunjukkan bentuk kurva yang menyerupai lonceng simetris. Secara visual, sisaan terindikasi berdistribusi normal dan saling bebas (white noise).
sisaan <- ARIMA100diag$residuals
# Uji formal normalitas data
ks.test(sisaan,"pnorm")##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.074578, p-value = 0.5278
## alternative hypothesis: two-sidedUji Normalitas (Kolmogorov-Smirnov) Pengujian normalitas secara analitik menggunakan Kolmogorov-Smirnov Test memiliki rumusan hipotesis sebagai berikut:
\(H_0\): Sisaan berdistribusi normal
\(H_1\): Sisaan tidak berdistribusi normal
Berdasarkan hasil Kolmogorov-Smirnov Test di atas, didapatkan nilai p-value > \(\alpha = 0.05\), maka keputusan yang diambil adalah terima \(H_0\). Secara matematis, hal ini membuktikan bahwa sisaan dari model menyebar secara normal.
# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: sisaan
## t = -0.74324, df = 117, p-value = 0.4588
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2448570 0.1112232
## sample estimates:
## mean of x
## -0.06681686Uji Nilai Tengah Sisaan (T-Test) Pengujian ini menggunakan One Sample t-test untuk memastikan keacakan sisaan dengan hipotesis:
\(H_0 : \mu = 0\) (Nilai tengah sisaan sama dengan nol)
\(H_1 : \mu \neq 0\) (Nilai tengah sisaan tidak sama dengan nol)
Berdasarkan hasil uji-t nilai p-value > \(\alpha = 0.05\), maka terima \(H_0\). Ini menunjukkan bahwa secara statistik, nilai harapan (expected value) dari galat terbukti sama dengan nol.
# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")##
## Box-Ljung test
##
## data: sisaan
## X-squared = 13.535, df = 23, p-value = 0.9394Uji Asumsi White Noise (Ljung-Box) Pengujian kebebasan sisaan dilakukan melalui Ljung-Box Test dengan rumusan hipotesis:
\(H_0\): Sisaan saling bebas / tidak terdapat autokorelasi (white noise)
\(H_1\): Sisaan terautokorelasi
Hasil pengujian memberikan nilai p-value > \(\alpha = 0.05\), maka diambil keputusan terima \(H_0\). Hal ini menyimpulkan bahwa sisaan bersifat saling bebas / tidak terdapat autokorelasi (white noise)