# Memuat library yang dibutuhkan
library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library (knitr)
library (kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.3
# 1. PERSIAPAN PARAMETER & SKENARIO
set.seed(2026) 
n_replikasi <- 100
true_ar1 <- 0.7

# Skenario 1: Variasi Ukuran Sampel
ukuran_sampel <- c(50, 100, 200) 
# Kecil, Sedang, Besar

# Tempat untuk menyimpan hasil
hasil_simulasi <- data.frame()
# 2. PROSES SIMULASI MONTE CARLO & ESTIMASI
for (n in ukuran_sampel) {
  for (i in 1:n_replikasi) {
    
    # a. Membangkitkan Data Normal (AR(1) / ARIMA(1,0,0))
    data_asli <- arima.sim(model = list(ar = true_ar1), n = n)
    
    # b. Skenario 2: Menyisipkan Outlier (5%)
    data_outlier <- data_asli
    idx_outlier <- sample(1:n, size = floor(0.05 * n))
    data_outlier[idx_outlier] <- data_outlier[idx_outlier] + rnorm(length(idx_outlier), mean = 10, sd = 2)
    
    # c. Skenario 3: Menyisipkan Missing Value (5% NA)
    data_missing <- data_asli
    idx_missing <- sample(1:n, size = floor(0.05 * n))
    data_missing[idx_missing] <- NA
    
    # --- PROSES ESTIMASI PARAMETER ---
    # Model 1: Data Normal
    fit_normal <- arima(data_asli, order = c(1,0,0), include.mean = FALSE)
    est_normal <- coef(fit_normal)["ar1"]
    
    # Model 2: Data dengan Outlier
    fit_outlier <- arima(data_outlier, order = c(1,0,0), include.mean = FALSE)
    est_outlier <- coef(fit_outlier)["ar1"]
    
    # Model 3: Data dengan Missing Value (Kalman Filter internal)
    fit_missing <- arima(data_missing, order = c(1,0,0), include.mean = FALSE)
    est_missing <- coef(fit_missing)["ar1"]
    
    # --- MENYIMPAN HASIL KE DATA FRAME ---
    temp_df <- data.frame(
      Replikasi = i,
      Ukuran_Sampel = n,
      Est_Normal = est_normal,
      Est_Outlier = est_outlier,
      Est_Missing = est_missing
    )
    hasil_simulasi <- rbind(hasil_simulasi, temp_df)
  }
}
head(hasil_simulasi)
##      Replikasi Ukuran_Sampel Est_Normal Est_Outlier Est_Missing
## ar1          1            50  0.7811526  0.16401425   0.7810390
## ar11         2            50  0.7529569  0.53659186   0.7397839
## ar12         3            50  0.7225699  0.16361118   0.7096838
## ar13         4            50  0.6103371  0.40125144   0.5862750
## ar14         5            50  0.7294486  0.46792616   0.7227662
## ar15         6            50  0.5157346  0.08734972   0.5224046
# 3. EVALUASI DAN ANALISIS BIAS

# Menghitung Bias (Estimasi - True Parameter)
hasil_simulasi$Bias_Normal <- hasil_simulasi$Est_Normal - true_ar1
hasil_simulasi$Bias_Outlier <- hasil_simulasi$Est_Outlier - true_ar1
hasil_simulasi$Bias_Missing <- hasil_simulasi$Est_Missing - true_ar1

# Menghitung Rata-rata Estimasi dan Bias per Skenario Ukuran Sampel
ringkasan_hasil <- aggregate(
  cbind(Est_Normal, Bias_Normal, Est_Outlier, Bias_Outlier, Est_Missing, Bias_Missing) ~ Ukuran_Sampel, 
  data = hasil_simulasi, 
  FUN = mean
)

# Merapikan Tabel Output
kable(ringkasan_hasil, 
      digits = 3, # Membulatkan angka menjadi 3 di belakang koma
      col.names = c("Ukuran Sampel (n)", 
                    "Estimasi Normal", "Bias Normal", 
                    "Estimasi Outlier", "Bias Outlier", 
                    "Estimasi Missing", "Bias Missing"),
      caption = "<b>Tabel 1. Ringkasan Rata-rata Estimasi dan Bias Parameter AR(1) dari 100 Replikasi</b>",
      align = "c") %>% # Membuat semua teks rata tengah (center)

# Menambahkan tema styling dari kableExtra
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
              full_width = FALSE, 
              position = "center") %>%

# Mengelompokkan header atas (Header bertingkat) agar lebih mudah dibaca
add_header_above(c(" " = 1, "Skenario Normal" = 2, "Skenario Outlier" = 2, "Skenario Missing Value" = 2)) %>%

# Memberikan warna khusus pada baris Skenario Outlier (opsional, untuk penekanan)
  row_spec(0, bold = TRUE, color = "white", background = "#2C3E50") %>% # Warna header
  column_spec(1, bold = TRUE) # Menebalkan kolom ukuran sampel
Tabel 1. Ringkasan Rata-rata Estimasi dan Bias Parameter AR(1) dari 100 Replikasi
Skenario Normal
Skenario Outlier
Skenario Missing Value
Ukuran Sampel (n) Estimasi Normal Bias Normal Estimasi Outlier Bias Outlier Estimasi Missing Bias Missing
50 0.690 -0.010 0.248 -0.452 0.689 -0.011
100 0.679 -0.021 0.220 -0.480 0.679 -0.021
200 0.691 -0.009 0.220 -0.480 0.691 -0.009
# Visualisasi sebaran estimasi menggunakan Boxplot
par(mfrow=c(1,3)) # Menyiapkan jendela plot 1 baris 3 kolom
boxplot(hasil_simulasi$Est_Normal ~ hasil_simulasi$Ukuran_Sampel, 
        main = "Estimasi Normal", xlab = "N", ylab = "Estimasi phi_1", col="lightblue")
abline(h = true_ar1, col = "red", lwd = 2, lty = 2)

boxplot(hasil_simulasi$Est_Outlier ~ hasil_simulasi$Ukuran_Sampel, 
        main = "Estimasi Outlier", xlab = "N", ylab = "Estimasi phi_1", col="salmon")
abline(h = true_ar1, col = "red", lwd = 2, lty = 2)

boxplot(hasil_simulasi$Est_Missing ~ hasil_simulasi$Ukuran_Sampel, 
        main = "Estimasi Missing", xlab = "N", ylab = "Estimasi phi_1", col="lightgreen")
abline(h = true_ar1, col = "red", lwd = 2, lty = 2)