1 Pendahuluan

1.1 Latar Belakang

Analisis deret waktu (time series) merupakan salah satu bidang statistika yang berfokus pada pemodelan data yang diobservasi secara berurutan dalam waktu. Salah satu kelas model yang paling fundamental adalah Model Autoregresif (AR), di mana nilai saat ini dari suatu variabel dijelaskan oleh nilai-nilai sebelumnya (lag) dari variabel yang sama.

Pada analisis ini, kita menggunakan Model AR orde 4 (AR(4)), artinya nilai \(y_t\) dipengaruhi oleh empat observasi sebelumnya:

\[y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma^2)\]

di mana:

  • \(\phi_1, \phi_2, \phi_3, \phi_4\) adalah koefisien autoregresif yang mengukur pengaruh masing-masing lag
  • \(\varepsilon_t\) adalah galat/noise yang diasumsikan berdistribusi normal dengan rata-rata 0 dan variansi \(\sigma^2\)

1.2 Mengapa Pendekatan Bayesian?

Berbeda dengan pendekatan frekuentis yang menghasilkan estimasi titik (point estimate), inferensi Bayesian menghasilkan distribusi posterior lengkap untuk setiap parameter. Hal ini memberikan beberapa keunggulan:

  1. Kuantifikasi ketidakpastian yang lebih kaya — kita mendapatkan interval kredibel, bukan sekadar confidence interval
  2. Integrasi pengetahuan prior — informasi sebelumnya dapat dimasukkan secara formal melalui distribusi prior
  3. Propagasi ketidakpastian ke prediksi — forecast yang dihasilkan merupakan distribusi prediktif, bukan hanya nilai titik

1.3 Tujuan Analisis

Analisis ini bertujuan untuk:

  1. Membangkitkan data simulasi dari proses AR(4) dengan parameter yang diketahui
  2. Mengestimasi parameter model AR(4) menggunakan metode Bayesian via MCMC (Markov Chain Monte Carlo) dengan bantuan Stan melalui paket rstan
  3. Mengevaluasi kualitas estimasi dengan membandingkan posterior mean terhadap nilai parameter sebenarnya
  4. Melakukan peramalan (forecasting) 100 langkah ke depan menggunakan distribusi posterior predictive

1.4 Parameter Simulasi

Sesuai dengan ilustrasi pada slide kuliah, parameter yang digunakan adalah:

\[\phi_1 = 0.8, \quad \phi_2 = -0.64, \quad \phi_3 = 0.512, \quad \phi_4 = -0.4096, \quad \sigma^2 = 4\]

2 Persiapan

2.1 Pemuatan Paket

library(MCMCpack)
library(bayesplot)
library(ggplot2)
library(dplyr)

3 Simulasi Data AR(4)

3.1 Pembangkitan Data

Kita membangkitkan 500 observasi dari proses AR(4) dengan parameter yang sudah ditentukan. Seed dipatok agar hasil dapat direproduksi.

set.seed(1061)

# --- Parameter true yang akan diestimasi ---
n_obs    <- 500
par_phi1 <- 0.8
par_phi2 <- -0.64
par_phi3 <- 0.512
par_phi4 <- -0.4096
par_sigma <- sqrt(4)   # standar deviasi, bukan variansi

# --- Bangkitkan deret AR(4) ---
seri <- numeric(n_obs)

for (t in 5:n_obs) {
  seri[t] <- par_phi1 * seri[t-1] +
             par_phi2 * seri[t-2] +
             par_phi3 * seri[t-3] +
             par_phi4 * seri[t-4] +
             rnorm(1, mean = 0, sd = par_sigma)
}

cat("Ringkasan statistik deret waktu simulasi:\n")
## Ringkasan statistik deret waktu simulasi:
summary(seri)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -7.937146 -1.996888 -0.003371 -0.057965  1.823308  8.788618

Berdasarkan hasil pembangkitan data simulasi untuk proses Autoregresif orde 4 (AR(4)) dengan total 500 observasi, ringkasan statistik menunjukkan bahwa data yang terbentuk memiliki karakteristik yang stabil dan stasioner di sekitar nilai nol. Hal ini ditunjukkan oleh nilai rata-rata (mean) sebesar -0.057965 dan median sebesar -0.003371, di mana kedua nilai tersebut sangat dekat dengan angka nol sebagai pusat distribusinya. Fluktuasi data bergerak secara simetris dalam rentang yang wajar, dengan nilai minimum berada di -7.937146 dan nilai maksimum mencapai 8.788618. Selain itu, sebaran kuartil pertama (1st Qu.) sebesar -1.996888 dan kuartil ketiga (3rd Qu.) sebesar 1.823308 menegaskan bahwa sebagian besar observasi mengelompok rapat di sekitar nilai tengah. Karakteristik statis yang berpusat di nol ini mengonfirmasi secara empiris bahwa proses AR(4) yang dibangkitkan memenuhi asumsi kestabilan deret waktu, sehingga siap dan layak untuk digunakan dalam tahapan pemodelan estimasi parameter selanjutnya.

3.2 Visualisasi Data Simulasi

df_seri <- data.frame(
  waktu = 1:n_obs,
  nilai = seri
)

ggplot(df_seri, aes(x = waktu, y = nilai)) +
  geom_line(colour = "#2c7bb6", linewidth = 0.6, alpha = 0.85) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey50") +
  labs(
    title    = "Data Simulasi Proses AR(4)",
    subtitle = paste0("n = ", n_obs, " | φ₁=0.8, φ₂=-0.64, φ₃=0.512, φ₄=-0.4096, σ²=4"),
    x        = "Indeks Waktu (t)",
    y        = "y(t)"
  ) +
  theme_bw(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold"),
    plot.subtitle = element_text(colour = "grey40")
  )

Visualisasi plot deret waktu menunjukkan perilaku data simulasi dari proses Autoregresif orde 4 (AR(4)) sepanjang 500 indeks waktu. Grafik ini mengonfirmasi bahwa data yang dibangkitkan bersifat stasioner, ditandai dengan fluktuasi nilai \(y(t)\) yang bergerak secara acak dan konsisten di sekitar garis rata-rata nol (ditunjukkan oleh garis putus-putus abu-abu). Selain itu, tidak terdapat tren naik ataupun turun yang monoton, serta varians data terlihat stabil di sepanjang waktu pengamatan. Pola naik-turun yang rapat dan berulang ini secara visual mencerminkan adanya ketergantungan yang kuat antara nilai saat ini dengan nilai-nilai pada empat lag waktu sebelumnya, sesuai dengan karakteristik teoritis dari model AR(4) yang diaplikasikan.

4 Pemodelan Bayesian AR(4)

4.1 Pembagian Data: Training dan Forecasting

Kita gunakan 400 observasi pertama sebagai data latih (training), sementara 100 observasi sisanya menjadi periode peramalan yang akan dibandingkan dengan prediksi model.

n_train    <- 400
n_forecast <- n_obs - n_train   # = 100

data_train <- seri[1:n_train]

# Data yang dikirim ke Stan
input_stan <- list(
  N = n_train,
  y = data_train
)

cat(sprintf("Data latih  : t = 1 s/d %d\n", n_train))
## Data latih  : t = 1 s/d 400
cat(sprintf("Data tes    : t = %d s/d %d\n", n_train + 1, n_obs))
## Data tes    : t = 401 s/d 500
cat(sprintf("Horizon forecast: %d langkah\n", n_forecast))
## Horizon forecast: 100 langkah

Sebelum masuk ke tahap estimasi parameter, data deret waktu yang berjumlah 500 observasi terlebih dahulu dibagi menjadi dua bagian secara berurutan (time-based split). Bagian pertama adalah data latih (training data) yang mencakup observasi dari \(t = 1\) sampai \(t = 400\), yang berfungsi sebagai basis data untuk mengestimasi distribusi posterior dari parameter model AR(4). Bagian kedua adalah data tes (testing data) dari \(t = 401\) sampai \(t = 500\), yang sengaja disimpan dan tidak diikutkan dalam pemodelan. Data tes sepanjang 100 langkah waktu ini nantinya akan digunakan sebagai standar evaluasi (ground truth) untuk diuji kesesuaiannya dengan hasil peramalan (horizon forecast) yang diproduksi oleh model Bayesian.

4.2 Penyusunan Data Matriks Lag (Desain Matriks)

Model AR(4) pada dasarnya bisa diestimasi menggunakan fungsi regresi linear berganda jika kita memecah data menjadi variabel dependen \(y_t\) dan variabel independen berupa lag 1 sampai lag 4.

# Membuat data frame untuk regresi linear lag
N <- length(data_train)
y_target <- data_train[5:N]
lag1     <- data_train[4:(N-1)]
lag2     <- data_train[3:(N-2)]
lag3     <- data_train[2:(N-3)]
lag4     <- data_train[1:(N-4)]

df_reg_train <- data.frame(
  y = y_target,
  lag1 = lag1,
  lag2 = lag2,
  lag3 = lag3,
  lag4 = lag4
)

4.3 Sampling MCMC dengan MCMCpack

# Jalankan MCMC tanpa Stan compiler
hasil_fit <- MCMCregress(
  y ~ lag1 + lag2 + lag3 + lag4 - 1, # -1 artinya tanpa intersep (sesuai rumus AR)
  data   = df_reg_train,
  burnin = 8000,
  mcmc   = 8000,
  thin   = 1,
  seed   = 42
)

# Mengubah nama kolom agar sesuai parameter
colnames(hasil_fit) <- c("koef1", "koef2", "koef3", "koef4", "var_galat")

5 Hasil Estimasi Posterior

5.1 Ringkasan Posterior

# Ambil summary dari objek MCMCpack
ringkasan_mcmc <- summary(hasil_fit)

# Hitung kuantil manual untuk tabel yang rapi
stats_raw <- ringkasan_mcmc$statistics
quant_raw <- ringkasan_mcmc$quantiles

tabel_posterior <- data.frame(
  Parameter   = colnames(hasil_fit),
  Nilai_True  = c(par_phi1, par_phi2, par_phi3, par_phi4, par_sigma^2),
  Post_Mean   = round(stats_raw[, "Mean"], 4),
  Post_Median = round(quant_raw[, "50%"], 4),
  CI_2.5      = round(quant_raw[, "2.5%"], 4),
  CI_97.5     = round(quant_raw[, "97.5%"], 4),
  Geweke_p    = round(2 * (1 - pnorm(abs(geweke.diag(hasil_fit)$z))), 4) # Alternatif Rhat untuk 1 rantai
)

tabel_posterior

Interpretasi: Hasil estimasi Bayesian menunjukkan bahwa model berhasil menduga seluruh parameter dengan akurat, di mana nilai rata-rata posterior (Post_Mean) berada sangat dekat dengan nilai parameter sejati (Nilai_True). Secara rinci, koefisien autoregresif terestimasi sebesar \(\phi_1 = 0.8161\), \(\phi_2 = -0.6306\), \(\phi_3 = 0.3978\), dan \(\phi_4 = -0.3225\), dengan variansi galat (\(\sigma^2\)) sebesar 4.4554. Ketepatan estimasi ini diperkuat oleh nilai sejati masing-masing parameter yang seluruhnya berhasil masuk ke dalam rentang Interval Kredibel 95% (antara CI_2.5 hingga CI_97.5), menandakan ketidakpastian model yang sangat terukur. Selain itu, berdasarkan nilai p-value dari uji diagnostik Geweke (Geweke_p) yang seluruhnya jauh di atas taraf signifikansi 0.05 (p > 0.05), dapat disimpulkan secara formal bahwa rantai algoritma MCMC telah mencapai kondisi stasioner dan konvergen dengan sangat baik.

5.2 Perbandingan Estimasi vs Nilai Sejati

selisih_df <- data.frame(
  Parameter  = c("phi1", "phi2", "phi3", "phi4", "sigma2"),
  Nilai_True = c(par_phi1, par_phi2, par_phi3, par_phi4, par_sigma^2),
  Post_Mean  = round(tabel_posterior$Post_Mean, 4),
  Selisih    = round(abs(tabel_posterior$Post_Mean -
                         c(par_phi1, par_phi2, par_phi3, par_phi4, par_sigma^2)), 4)
)

selisih_df

Tabel perbandingan antara nilai sejati (Nilai_True) dan rata-rata posterior (Post_Mean) menunjukkan bahwa model Bayesian AR(4) memiliki tingkat akurasi pendugaan yang sangat tinggi. Kedekatan hasil estimasi ini ditunjukkan oleh nilai selisih absolut (Selisih) yang sangat kecil untuk seluruh parameter. Koefisien \(\phi_1\) dan \(\phi_2\) memiliki akurasi paling tinggi dengan selisih masing-masing hanya sebesar 0.0161 dan 0.0094 dari nilai sejatinya. Sementara itu, untuk koefisien lag yang lebih dalam yaitu \(\phi_3\) dan \(\phi_4\), selisihnya sedikit meningkat menjadi 0.1142 dan 0.0871, yang merupakan hal wajar dalam estimasi deret waktu karena informasi dari lag yang lebih jauh cenderung memiliki ketidakpastian yang lebih besar. Terakhir, variansi galat (\(\sigma^2\)) terestimasi sebesar 4.4554 dengan selisih 0.4554 dari nilai sejati 4.0000. Secara keseluruhan, nilai selisih yang berada di bawah angka 0.5 ini membuktikan bahwa algoritma MCMC berbasis regresi lag mampu memulihkan (recover) parameter struktural dari data simulasi AR(4) secara optimal.

6 Diagnostik Konvergensi MCMC

6.1 Trace Plot

Trace plot yang baik menunjukkan mixing yang baik — rantai bergerak bebas di sekitar nilai stasioner tanpa tren atau stuck.

mcmc_trace(
  as.array(hasil_fit)
) +
  ggtitle("Trace Plot — AR(4) Bayesian (MCMCpack)") +
  theme_bw()

Grafik trace plot di atas menunjukkan visualisasi pergerakan sampel simulasi MCMC sepanjang 8.000 iterasi untuk kelima parameter model AR(4). Secara keseluruhan, kelima plot menunjukkan pola bauran yang sangat ideal, menyerupai bentuk ulat bulu (fuzzy caterpillar plot). Rantai bergerak secara konsisten dan acak di sekitar nilai stasionernya tanpa memperlihatkan adanya tren naik-turun yang struktural, pola siklis, ataupun gejala macet (stuck). Karakteristik visual yang stabil dan berfluktuasi cepat ini mengonfirmasi secara empiris bahwa algoritma Gibbs Sampler dari MCMCpack telah menjelajahi ruang parameter secara efisien dan mencapai konvergensi yang sempurna sejak awal iterasi pasca-burn-in.

6.2 Distribusi Posterior

library(MCMCpack)
library(coda)

# 1. Buat list kosong untuk menampung rantai MCMC
rantai_list <- list()

# 2. Jalankan 3 rantai dengan seed berbeda
for(i in 1:3) {
  fit_rantai <- MCMCregress(
    y ~ lag1 + lag2 + lag3 + lag4 - 1,
    data   = df_reg_train,
    burnin = 8000,
    mcmc   = 8000,
    seed   = 42 + i
  )
  # Mengubah nama kolom di setiap rantai secara konsisten sebelum digabung
  colnames(fit_rantai) <- c("koef1", "koef2", "koef3", "koef4", "var_galat")
  rantai_list[[i]] <- fit_rantai
}

# 3. Gabungkan menjadi objek mcmc.list
hasil_fit <- mcmc.list(rantai_list)

# 4. DEFINISIKAN ARRAY AMAN DI SINI
n_iter  <- nrow(hasil_fit[[1]])
n_chain <- length(hasil_fit)
n_param <- ncol(hasil_fit[[1]])

array_aman <- array(NA, dim = c(n_iter, n_chain, n_param))
for(i in 1:n_chain) {
  array_aman[, i, ] <- as.matrix(hasil_fit[[i]])
}
dimnames(array_aman) <- list(NULL, NULL, colnames(hasil_fit[[1]]))
library(bayesplot)

# Mengubah ke matriks agar nama parameter aman
matriks_fit <- as.matrix(hasil_fit)

# Menggunakan mcmc_dens (bukan dens_overlay) agar aman untuk data matriks
mcmc_dens(
  matriks_fit
) +
  ggtitle("Distribusi Posterior — Setiap Parameter AR(4)") +
  theme_bw()

Grafik plot densitas di atas menggambarkan bentuk distribusi posterior penuh untuk kelima parameter model AR(4) hasil simulasi MCMC. Seluruh kurva densitas, baik untuk parameter koefisien autoregresif (koef1 hingga koef4) maupun variansi galat (var_galat), secara visual membentuk kurva unipolar tunggal yang simetris menyerupai distribusi normal (lonceng terbalik). Bentuk kurva yang mulus, kokoh, dan tidak memiliki pola multimodal (banyak puncak) ini mengonfirmasi secara empiris bahwa proses penaksiran parameter berjalan sangat stabil dan tidak mengalami bias. Selain itu, puncak kurva yang tajam (narrow band) mencerminkan tingkat presisi estimasi yang tinggi, di mana area konsentrasi massa probabilitas terbesar terkunci tepat di sekitar nilai rata-rata posterior masing-masing parameter.

6.3 Autocorrelation Function (ACF) Sampel MCMC

ACF yang turun cepat menuju nol menandakan sampel yang hampir tidak berkorelasi — ideal untuk inferensi.

library(bayesplot)
mcmc_acf(
  array_aman,
  pars = c("koef1", "var_galat"),
  lags = 40
) +
  theme_bw()

Grafik di atas menampilkan plot autokorelasi (ACF) dari sampel MCMC untuk parameter koef1 dan var_galat di ketiga rantai simulasi hingga lag ke-40. Seluruh panel grafik menunjukkan pola penurun nilai autokorelasi yang sangat ideal, di mana nilai korelasi langsung anjlok tajam mendekati angka nol tepat pada lag pertama (\(t=1\)) dan terus bertahan stabil di sekitar garis nol pada lag-lag berikutnya. Pola peluruhan yang sangat cepat (die out) ini mengindikasikan bahwa sampel-sampel parameter yang dihasilkan oleh algoritma Gibbs Sampler bersifat independen (hampir tidak memiliki korelasi antar-iterasi). Karakteristik rantai yang bersih dari autokorelasi ini membuktikan secara empiris bahwa penjelajahan ruang parameter berjalan sangat efisien, rantai tidak mengalami gejala lambat bergerak (sluggish mixing), dan informasi distribusi posterior yang diperoleh memiliki kualitas yang sangat tepercaya untuk inferensi model.

7 Peramalan (Forecasting) Bayesian

7.1 Ekstraksi Sampel Posterior

# Mengubah objek mcmc.list menjadi matriks datar agar mudah diekstrak kolomnya
mat_sampel <- as.matrix(hasil_fit)

# Mengekstrak sampel parameter berdasarkan nama kolom matriks
vec_koef1    <- mat_sampel[, "koef1"]
vec_koef2    <- mat_sampel[, "koef2"]
vec_koef3    <- mat_sampel[, "koef3"]
vec_koef4    <- mat_sampel[, "koef4"]
vec_vargalat <- mat_sampel[, "var_galat"]

n_sampel <- length(vec_koef1)
cat("Jumlah sampel posterior yang diekstrak:", n_sampel, "\n")
## Jumlah sampel posterior yang diekstrak: 24000

7.2 Posterior Predictive Forecasting

Untuk setiap sampel dari distribusi posterior, kita simulasikan 100 langkah ke depan menggunakan rekursi AR(4). Hasilnya adalah distribusi prediktif (posterior predictive distribution).

# Matriks untuk menyimpan forecast: baris = sampel, kolom = horizon
mat_forecast <- matrix(NA, nrow = n_sampel, ncol = n_forecast)

# 4 nilai terakhir dari training sebagai inisialisasi rolling window
init_window <- seri[(n_train - 3):n_train]  # y[397], y[398], y[399], y[400]

for (s in 1:n_sampel) {
  window <- init_window   # reset rolling window untuk setiap sampel

  for (h in 1:n_forecast) {
    y_baru <- vec_koef1[s] * window[4] +
              vec_koef2[s] * window[3] +
              vec_koef3[s] * window[2] +
              vec_koef4[s] * window[1] +
              rnorm(1, 0, sqrt(vec_vargalat[s]))

    mat_forecast[s, h] <- y_baru

    # Geser window: buang yang paling lama, masukkan nilai baru
    window <- c(window[2], window[3], window[4], y_baru)
  }
}

7.3 Ringkasan Distribusi Prediktif

pred_mean  <- apply(mat_forecast, 2, mean)
pred_lower <- apply(mat_forecast, 2, quantile, probs = 0.025)
pred_upper <- apply(mat_forecast, 2, quantile, probs = 0.975)
pred_med   <- apply(mat_forecast, 2, quantile, probs = 0.50)

df_forecast <- data.frame(
  t          = (n_train + 1):n_obs,
  pred_mean  = pred_mean,
  pred_med   = pred_med,
  pred_lower = pred_lower,
  pred_upper = pred_upper
)

head(df_forecast, 10)

Tabel di atas menampilkan hasil ringkasan dari Distribusi Prediktif Posterior (Posterior Predictive Distribution) untuk 10 langkah pertama pada periode peramalan (horizon forecast), yaitu dari indeks waktu \(t = 401\) hingga \(t = 410\).

7.4 Visualisasi Forecast

df_aktual <- data.frame(
  t     = 1:n_obs,
  nilai = seri,
  tipe  = ifelse(1:n_obs <= n_train, "Training", "Aktual (tes)")
)

ggplot() +
  # Seluruh data aktual
  geom_line(
    data    = df_aktual,
    aes(x = t, y = nilai, colour = tipe),
    linewidth = 0.55
  ) +
  # Interval kredibel 95%
  geom_ribbon(
    data  = df_forecast,
    aes(x = t, ymin = pred_lower, ymax = pred_upper),
    fill  = "#fdae61",
    alpha = 0.40
  ) +
  # Garis forecast (mean posterior)
  geom_line(
    data      = df_forecast,
    aes(x = t, y = pred_mean, colour = "Forecast (mean post.)"),
    linewidth = 1.0,
    linetype  = "dashed"
  ) +
  # Garis pemisah train/forecast
  geom_vline(
    xintercept = n_train,
    linetype   = "dotted",
    colour     = "grey30",
    linewidth  = 0.8
  ) +
  annotate("text", x = n_train - 10, y = max(seri) * 0.9,
           label = "Train", hjust = 1, colour = "grey30", size = 3.5) +
  annotate("text", x = n_train + 10, y = max(seri) * 0.9,
           label = "Forecast", hjust = 0, colour = "grey30", size = 3.5) +
  scale_colour_manual(
    values = c(
      "Training"             = "#4575b4",
      "Aktual (tes)"         = "#d73027",
      "Forecast (mean post.)" = "#1a9850"
    )
  ) +
  labs(
    title    = "Peramalan Bayesian AR(4) — 100 Langkah ke Depan",
    subtitle = "Pita kuning = Interval Kredibel 95% | Garis hijau putus = Posterior Predictive Mean",
    x        = "Indeks Waktu (t)",
    y        = "y(t)",
    colour   = NULL
  ) +
  theme_bw(base_size = 12) +
  theme(
    legend.position   = "bottom",
    plot.title        = element_text(face = "bold"),
    plot.subtitle     = element_text(colour = "grey40", size = 10)
  )

Grafik di atas memvisualisasikan performa peramalan Bayesian model AR(4) untuk 100 langkah ke depan (\(t = 401\) s/d \(t = 500\)) berbatasan dengan data latihnya. Garis hijau putus-putus menggambarkan nilai rata-rata prediktif posterior (Posterior Predictive Mean), yang terlihat bergerak stabil menuju nilai rata-rata jangka panjang deret waktu (yaitu nol) seiring bertambahnya horizon ramalan, mencerminkan sifat stasioneritas yang melekat pada model AR. Keunggulan utama dari pendekatan Bayesian ini ditunjukkan oleh pita kuning transparan yang membentuk Interval Kredibel 95%. Pita ini terbukti mampu menangkap dan membungkus mayoritas fluktuasi data aktual periode tes (garis merah) secara proporsional. Pola pita yang lebar secara konstan ini mengonfirmasi bahwa model tidak hanya memberikan tebakan titik, melainkan berhasil menyebarkan ketidakpastian parameter secara jujur untuk memayungi risiko pergerakan data acak di masa depan.

8 Evaluasi Akurasi Forecast

aktual_tes <- seri[(n_train + 1):n_obs]

mae  <- mean(abs(pred_mean - aktual_tes))
rmse <- sqrt(mean((pred_mean - aktual_tes)^2))
mape <- mean(abs((pred_mean - aktual_tes) / aktual_tes)) * 100

# Proporsi observasi yang masuk ke dalam interval kredibel 95%
cakupan_95 <- mean(aktual_tes >= pred_lower & aktual_tes <= pred_upper) * 100

data.frame(
  Metrik = c("MAE", "RMSE", "MAPE (%)", "Cakupan CI-95% (%)"),
  Nilai  = round(c(mae, rmse, mape, cakupan_95), 4)
)

Tabel di atas menyajikan metrik evaluasi akurasi peramalan untuk menguji seberapa dekat hasil prediksi model dengan data aktual periode tes. Berdasarkan nilai rata-rata kesalahan absolut, model menghasilkan nilai MAE sebesar 2.3165 dan RMSE sebesar 2.9599. Sementara itu, persentase rata-rata kesalahan absolutnya atau MAPE menunjukkan angka 98.6596%. Metrik evaluasi paling krusial dalam pendekatan Bayesian ditunjukkan oleh Cakupan CI-95% yang mencapai 94.0000%, yang berarti dari 100 langkah data aktual di masa depan, ada 94 data yang berhasil masuk ke dalam koridor interval kredibel yang dipetakan oleh model.

9 Kesimpulan

Analisis Bayesian pada Model Autoregresif orde 4 (AR(4)) menggunakan pendekatan praktis tanpa Stan (MCMCpack) ini menghasilkan beberapa temuan utama:

1. Kualitas Estimasi Parameter

Secara keseluruhan, model Bayesian berhasil menduga seluruh parameter struktural dari data simulasi dengan tingkat akurasi yang sangat tinggi. Hal ini dibuktikan oleh nilai rata-rata posterior (Posterior Mean) yang posisinya sangat dekat dengan nilai parameter sejati (True Value), dengan nilai selisih absolut yang berada di bawah angka 0.5.

Parameter Nilai True Posterior Mean Selisih Absolut
\(\phi_1\) 0.8000 0.8157 0.0157
\(\phi_2\) −0.6400 -0.6295 0.0105
\(\phi_3\) 0.5120 0.3973 0.1147
\(\phi_4\) −0.4096 -0.322 0.0876
\(\sigma^2\) 4.0000 4.4489 0.4489

Ketepatan estimasi ini diperkuat oleh fakta bahwa seluruh nilai sejati parameter sukses masuk ke dalam rentang Interval Kredibel 95%. Hal ini menunjukkan bahwa pemodelan MCMC berbasis regresi lag mampu memulihkan parameter asli data secara optimal.

2. Konvergensi dan Kualitas Rantai MCMC

Evaluasi diagnostik secara formal dan visual mengonfirmasi bahwa algoritma Gibbs Sampler telah mencapai kondisi stasioner dan konvergen yang sempurna: * Uji Geweke: Nilai p-value uji Geweke yang berada jauh di atas taraf signifikansi 0.05 (\(p > 0.05\)) untuk seluruh parameter membuktikan secara statistik bahwa rantai MCMC sudah stasioner. * Trace Plot & Densitas: Grafik trace plot menunjukkan bauran rantai yang ideal menyerupai “ulat bulu” tanpa tren struktural, didukung oleh plot densitas yang membentuk kurva unipolar tunggal simetris (menyerupai lonceng normal). * Autokorelasi (ACF): Grafik ACF yang langsung anjlok tajam mendekati nol sejak lag pertama menandakan sampel posterior yang dihasilkan bersifat independen dan bebas dari gejala sluggish mixing.

3. Evaluasi Kemampuan Peramalan (Forecasting)

Model Bayesian AR(4) menunjukkan performa peramalan 100 langkah ke depan yang sangat tangguh dalam mengonfirmasi ketidakpastian masa depan. Metrik akurasi menunjukkan nilai MAE sebesar 2.3125 dan RMSE sebesar 2.9578. Meskipun nilai MAPE tinggi (98.2897%) akibat sifat matematis data yang berfluktuasi di sekitar angka nol, model ini terbukti sangat objektif secara Bayesian. Hal ini divalidasi oleh tingkat Cakupan Interval Kredibel 95% yang mencapai 94% (sangat dekat dengan nilai teoritis 95%), yang berarti pita prediktif posterior mampu memayungi risiko pergerakan data aktual secara jujur dan proporsional.

4. Keunggulan Pendekatan

Dengan beralih menggunakan paket MCMCpack, analisis inferensi Bayesian berbasis MCMC terbukti dapat diselesaikan secara penuh, valid, dan akurat. Pendekatan ini menawarkan solusi operasional yang jauh lebih praktis dan ramah pengguna karena terbebas dari kendala teknis error kompilasi C++ eksternal (Rtools), tanpa mengorbankan kekayaan informasi distribusi posterior yang menjadi keunggulan utama mazhab Bayesian.