library(MASS)        
library(coda)  
library(ggplot2)
library(HDInterval)  
library(bayesplot)

## This is bayesplot version 1.15.0

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

set.seed(345)

Simulasi Data Time Series AR(4)

n_data = 500
phi_asli = c(0.8, -0.64, 0.512, -0.4096)
sigma2_asli = 4

# Vektor y
y = numeric(n_data)

# 4 Data Pertama Pembangkitan Acak
y[1:4] = rnorm(4, 0, sqrt(sigma2_asli))

# Looping pembangkitan AR(4) 
for (t in 5:n_data) {
  mu_t = phi_asli[1]*y[t-1] + phi_asli[2]*y[t-2] + phi_asli[3]*y[t-3] + phi_asli[4]*y[t-4]
  y[t] = mu_t + rnorm(1, 0, sqrt(sigma2_asli))
}

Gibbs Sampling

Persiapan

n_iter = 15000
burn_in = 5000
n_chains = 3      # DITAMBAHKAN: Menentukan jumlah rantai (chain)
p = 4

# Variabel respon (Y) dan Matriks Desain (X)
y_resp = y[5:n_data]
X_mat = matrix(0, nrow = length(y_resp), ncol = p)

for (i in 1:length(y_resp)) {
  # Mengambil y[t-1], y[t-2], y[t-3], y[t-4]
  X_mat[i, ] = c(y[i+3], y[i+2], y[i+1], y[i]) 
}

# Parameter Prior
v0 = 0.3          # Shape untuk Inverse Gamma
w0 = 3.3333       # Scale untuk Inverse Gamma
var_prior = 1000  # Varians untuk Normal
Sigma_prior_inv = diag(1/var_prior, p)

Gibss Sampler Dengan Multiple Chain (Nested Loop)

list_chains = list()

for (c in 1:n_chains) {
  
  # Set seed berbeda untuk tiap rantai agar pergerakannya independen
  set.seed(c * 345) 
  
  simpan_phi = matrix(0, nrow = n_iter, ncol = p)
  simpan_sigma2 = numeric(n_iter)
  
  phi_curr = rep(0, p)
  sigma2_curr = 1
  
# Iterasi MCMC
  for (iter in 1:(n_iter + burn_in)) {
    
    # A. Update Parameter Phi
    Sigma_n = solve((t(X_mat) %*% X_mat) / sigma2_curr + Sigma_prior_inv)
    mu_n = Sigma_n %*% ((t(X_mat) %*% y_resp) / sigma2_curr)
    phi_curr = mvrnorm(1, mu = mu_n, Sigma = Sigma_n)
    
    # B. Update Parameter Sigma2
    resid = y_resp - (X_mat %*% phi_curr)
    SSR = sum(resid^2)
    shape_post = v0 + (length(y_resp) / 2)
    scale_post = w0 + (SSR / 2)
    sigma2_curr = 1 / rgamma(1, shape = shape_post, rate = scale_post)
    
    # C. Simpan hasil 
    if (iter > burn_in) {
      simpan_phi[iter - burn_in, ] = phi_curr
      simpan_sigma2[iter - burn_in] = sigma2_curr
    }
  }
  
  # Format hasil per chain menjadi objek mcmc
  res_mcmc = mcmc(cbind(simpan_phi, sigma2 = simpan_sigma2))
  colnames(res_mcmc) = c("phi_1", "phi_2", "phi_3", "phi_4", "sigma2")
  list_chains[[c]] = res_mcmc
}

# Gabungkan seluruh chain menjadi mcmc.list
mcmc_gabungan = as.mcmc.list(list_chains)
matriks_sampel = as.matrix(mcmc_gabungan)

Diagnostik Konvergensi & Ringkasan

# Ringkasan MCMC
summary(mcmc_gabungan)

## 
## Iterations = 1:15000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 15000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## phi_1   0.7695 0.04135 0.0001949      0.0001942
## phi_2  -0.6821 0.04958 0.0002337      0.0002328
## phi_3   0.4806 0.05020 0.0002367      0.0002364
## phi_4  -0.4062 0.04181 0.0001971      0.0001971
## sigma2  4.1160 0.26362 0.0012427      0.0012427
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%     75%   97.5%
## phi_1   0.6882  0.7417  0.7696  0.7971  0.8505
## phi_2  -0.7799 -0.7156 -0.6821 -0.6486 -0.5843
## phi_3   0.3821  0.4468  0.4805  0.5145  0.5791
## phi_4  -0.4874 -0.4343 -0.4062 -0.3780 -0.3243
## sigma2  3.6340  3.9317  4.1064  4.2854  4.6642

# Effective Size MCMC
print(effectiveSize(mcmc_gabungan))

##    phi_1    phi_2    phi_3    phi_4   sigma2 
## 45352.86 45380.38 45134.96 45000.00 45000.00

# Diagnostik Gelman-Rubin (PSRF harus < 1.1)
print(gelman.diag(mcmc_gabungan))

## Potential scale reduction factors:
## 
##        Point est. Upper C.I.
## phi_1           1          1
## phi_2           1          1
## phi_3           1          1
## phi_4           1          1
## sigma2          1          1
## 
## Multivariate psrf
## 
## 1

Algoritma Gibbs Sampler berhasil mengestimasi parameter model dengan tingkat akurasi dan presisi yang sangat tinggi.

Akurasi Nilai Tengah (Mean): Nilai Posterior Mean untuk seluruh parameter berada sangat dekat dengan nilai parameter asli (parameter simulasi).

\(\phi_1 = 0.7695\) (Sangat dekat dengan nilai asli 0.8)
\(\phi_2 = -0.6821\) (Sangat dekat dengan nilai asli -0.64)
\(\phi_3 = 0.4806\) (Sangat dekat dengan nilai asli 0.512)
\(\phi_4 = -0.4062\) (Sangat dekat dengan nilai asli -0.4096)
\(\sigma^2 = 4.1160\) (Sangat dekat dengan varians asli 4)

Presisi (Standard Error & HPD): Nilai Time-series SE sangat kecil (berada di kisaran 0.0001 hingga 0.001), yang menandakan rendahnya galat standar dari estimasi MCMC. Selain itu, nilai parameter asli seluruhnya masuk ke dalam rentang interval kredibel 95% (kuantil 2.5% hingga 97.5%). Hal ini membuktikan bahwa model tidak bias dan spesifikasi prior yang diberikan sudah sesuai.

Visualisasi Diagnostik

# Visualisasi Diagnostik MCMC dengan tema warna khas bayesplot
color_scheme_set("viridis")

# Plot 1: Trace Plot 
plot_trace = mcmc_trace(mcmc_gabungan, facet_args = list(nrow = 5)) + 
  labs(title = "Trace Plot Parameter AR(4)")
print(plot_trace)

# Plot 2: Density Plot 
plot_density = mcmc_areas(mcmc_gabungan, prob = 0.95, point_est = "mean") + 
  labs(title = "Posterior Density Plot (95% HPD Interval)")
print(plot_density)

Gelman-Rubin Diagnostic (PSRF): Nilai taksiran titik (Point est.) dan taksiran multivariat (Multivariate psrf) untuk seluruh parameter menunjukkan angka tepat 1. Karena nilainya jauh di bawah ambang batas kritis (< 1.1), ini memastikan bahwa ketiga rantai (chains) yang dijalankan dari titik awal berbeda telah berbaur sepenuhnya dan berhasil konvergen ke satu distribusi posterior yang sama.
Effective Sample Size (ESS): Nilai ESS yang dihasilkan sangat optimal, yakni berkisar di angka 45.000 untuk seluruh parameter. Karena total iterasi sampel yang disimpan adalah 45.000 (3 chains × 15.000 iterasi), nilai ESS ini mengindikasikan bahwa hampir tidak ada autokorelasi antar-sampel (perfect mixing). Sampel yang dihasilkan benar-benar independen dan sangat informatif.
Trace Plot (Visual): Mengonfirmasi hasil numerik ESS dan PSRF, visualisasi Trace Plot menampilkan pola “ulat berbulu” (hairy caterpillar) yang tebal dan stabil di sekitar nilai rata-rata. Tidak terlihat adanya tren naik/turun atau stak pada satu titik, menandakan proses MCMC telah mencapai kondisi stasioner.
Kerapatan Posterior (Density Plot): Distribusi peluang untuk masing-masing parameter membentuk kurva lonceng yang halus dan unimodal (satu puncak tunggal) yang berpusat pada nilai Posterior Mean-nya. Bentuk distribusi pita bayangan yang sempit dan simetris ini memperkuat kesimpulan bahwa tingkat ketidakpastian model sangat rendah.

Prediksi (T>400)

waktu_plot = 401:n_data
mean_pred = numeric(length(waktu_plot))
hpd_bawah = numeric(length(waktu_plot))
hpd_atas = numeric(length(waktu_plot))

# Looping sederhana untuk memprediksi tiap titik waktu
for (i in 1:length(waktu_plot)) {
  t_skrg = waktu_plot[i]
  y_lalu = c(y[t_skrg-1], y[t_skrg-2], y[t_skrg-3], y[t_skrg-4])
  
  # Hitung rata-rata prediksi untuk seluruh iterasi MCMC
  mu_sampel = simpan_phi %*% y_lalu
  
  # Tambahkan efek varians error (sigma2)
  y_pred_sampel = mu_sampel + rnorm(n_iter, 0, sqrt(simpan_sigma2))
  
  # Simpan rata-rata dan interval HPD
  mean_pred[i] = mean(y_pred_sampel)
  hpd_val = hdi(y_pred_sampel, credMass = 0.95)
  hpd_bawah[i] = hpd_val[1]
  hpd_atas[i] = hpd_val[2]
}

Visualisasi Hasil Forecast

df_forecast = data.frame(
  t = waktu_plot,
  y = y[waktu_plot],
  mean = mean_pred,
  hpdlower = hpd_bawah,
  hpdupper = hpd_atas
)

plot_forecast = ggplot(df_forecast, aes(x = t)) +
  geom_ribbon(aes(ymin = hpdlower, ymax = hpdupper), fill = "grey70", alpha = 0.5) +
  geom_line(aes(y = y, color = "Data Aktual (Y)"), size = 0.7) +
  geom_line(aes(y = mean, color = "Posterior Mean"), linetype = "dashed", size = 0.8) +
  scale_color_manual(name = "Legenda:", values = c("Data Aktual (Y)" = "red", "Posterior Mean" = "blue")) +
  labs(
    title = "Peramalan Model AR(4) dengan MCMC Bayesian", 
    x = "Waktu (t)", 
    y = "Nilai Observasi"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    legend.position = "bottom"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

print(plot_forecast)

Grafik peramalan pada periode evaluasi (\(t > 400\)) menunjukkan performa model yang sangat baik dalam mereplikasi data observasi.

Pergerakan Nilai Tengah: Garis putus-putus biru (Posterior Mean) bergerak konsisten dan sangat merekat dengan pola fluktuasi garis merah (Data Aktual). Model AR(4) terbukti mampu menangkap pola autoregresif dari empat lag sebelumnya dengan sangat sensitif.
Kuantifikasi Ketidakpastian (HPD Band): Pita abu-abu transparan mewakili batas bawah dan batas atas dari ketidakpastian prediksi (95% HPD). Dapat diamati bahwa hampir seluruh lonjakan dan penurunan data aktual berada dengan aman di dalam area pita abu-abu tersebut. Model berhasil memberikan rentang jaring pengaman (safety mechanism) yang akurat tanpa menjadi terlalu lebar atau terlalu sempit.

Timeseries Bayesian (tugas PSB)

Muthi’ah Iffa Karimah

2026-05-31