# 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
# d. Skenario 4: KOMBINASI Outlier & Missing Value (Faktorial)
data_kombinasi <- data_asli
idx_out_komb <- sample(1:n, size = floor(0.05 * n))
idx_miss_komb <- sample(1:n, size = floor(0.05 * n))
# Injeksi independen
data_kombinasi[idx_out_komb] <- data_kombinasi[idx_out_komb] + rnorm(length(idx_out_komb), mean = 10, sd = 2)
data_kombinasi[idx_miss_komb] <- 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"]
# Model 4: Data Kombinasi (Outlier + Missing)
fit_kombinasi <- arima(data_kombinasi, order = c(1,0,0), include.mean = FALSE)
est_kombinasi <- coef(fit_kombinasi)["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,
Est_Kombinasi = est_kombinasi
)
hasil_simulasi <- rbind(hasil_simulasi, temp_df)
}
}
head(hasil_simulasi)
## Replikasi Ukuran_Sampel Est_Normal Est_Outlier Est_Missing Est_Kombinasi
## ar1 1 50 0.7811526 0.1640142 0.7810390 0.3373695
## ar11 2 50 0.7544650 0.3073940 0.7460712 0.5527449
## ar12 3 50 0.6653868 -0.1239402 0.6932031 0.3923494
## ar13 4 50 0.5097080 0.1356213 0.5036944 0.1825201
## ar14 5 50 0.7213433 0.4180318 0.7464278 0.3712395
## ar15 6 50 0.5674159 0.3093908 0.5715335 0.2459184
# 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
hasil_simulasi$Bias_Kombinasi <- hasil_simulasi$Est_Kombinasi - 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, Est_Kombinasi, Bias_Kombinasi) ~ Ukuran_Sampel,
data = hasil_simulasi,
FUN = mean
)
# Merapikan Tabel Output
kable(ringkasan_hasil,
digits = 3,
col.names = c("Ukuran Sampel (n)",
"Estimasi", "Bias",
"Estimasi", "Bias",
"Estimasi", "Bias",
"Estimasi", "Bias"),
caption = "<b>Tabel 1. Ringkasan Rata-rata Estimasi dan Bias Parameter AR(1) dari 100 Replikasi</b>",
align = "c") %>%
# Menambahkan tema styling dari kableExtra
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE,
position = "center",
font_size = 12) %>%
# Mengelompokkan header atas (Header bertingkat)
add_header_above(c(" " = 1, "Normal" = 2, "Outlier" = 2, "Missing Value" = 2, "Kombinasi" = 2)) %>%
# Memberikan penekanan pada ukuran sampel
column_spec(1, bold = TRUE)
Tabel 1. Ringkasan Rata-rata Estimasi dan Bias Parameter AR(1) dari
100 Replikasi
|
|
Normal
|
Outlier
|
Missing Value
|
Kombinasi
|
|
Ukuran Sampel (n)
|
Estimasi
|
Bias
|
Estimasi
|
Bias
|
Estimasi
|
Bias
|
Estimasi
|
Bias
|
|
50
|
0.670
|
-0.030
|
0.223
|
-0.477
|
0.669
|
-0.031
|
0.241
|
-0.459
|
|
100
|
0.691
|
-0.009
|
0.206
|
-0.494
|
0.690
|
-0.010
|
0.223
|
-0.477
|
|
200
|
0.696
|
-0.004
|
0.225
|
-0.475
|
0.695
|
-0.005
|
0.221
|
-0.479
|
# Visualisasi sebaran estimasi menggunakan Boxplot (Format 2x2)
par(mfrow=c(2,2))
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)
boxplot(hasil_simulasi$Est_Kombinasi ~ hasil_simulasi$Ukuran_Sampel,
main = "Estimasi Kombinasi", xlab = "N", ylab = "Estimasi phi_1", col="mediumpurple")
abline(h = true_ar1, col = "red", lwd = 2, lty = 2)
