R Notebook

UTS PEMODELAN STATISTIKA DAN SIMULASI

ALICE DELICA PERMATA

2404220017

library(forecast)

## Warning: package 'forecast' was built under R version 4.4.3

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

set.seed(2026)

# Parameter Awal
AR <- 0.7
MA <- 0.4
SAR <- 0.8
SMA <- 0.5
sigma_error <- 0.1
diff_non_s <- 1
diff_s <- 1

# Fungsi Simulasi
run_simulation <- function(n_samples, missing_pct = 0, n_replications = 100) {
  results <- data.frame(AR_est = numeric(), MA_est = numeric(),
                        SAR_est = numeric(), SMA_est = numeric())
  
  for (i in 1:n_replications) {
    # Membangkitkan Data (burn-in period digunakan agar data lebih stabil)
    burn_in <- 100
    data_sim <- arima.sim(n = n_samples + burn_in,
                          model = list(ar = AR, ma = MA),
                          sd = sigma_error)
    data_sim <- data_sim[(burn_in + 1):(burn_in + n_samples)]
    data_ts <- ts(data_sim, frequency = 12)
    
    # Skenario Missing Value
    if (missing_pct > 0) {
      indices <- sample(1:length(data_ts), size = floor(missing_pct * length(data_ts)))
      data_ts[indices] <- NA
      data_ts <- na.interp(data_ts) 
    }
    
    # Estimasi Parameter
    fit <- tryCatch({
      Arima(data_ts, order = c(1, diff_non_s, 1),
            seasonal = list(order = c(1, diff_s, 1), period = 12),
            method = "ML")
    }, error = function(e) { return(NULL) })
    
    if (is.null(fit)) { next }
    
    # Ekstrak nilai koefisien
    ar1_val  <- if("ar1" %in% names(fit$coef)) fit$coef["ar1"] else NA
    ma1_val  <- if("ma1" %in% names(fit$coef)) fit$coef["ma1"] else NA
    sar1_val <- if("sar1" %in% names(fit$coef)) fit$coef["sar1"] else NA
    sma1_val <- if("sma1" %in% names(fit$coef)) fit$coef["sma1"] else NA
    
    results <- rbind(results, data.frame(
      AR_est  = ar1_val,
      MA_est  = ma1_val,
      SAR_est = sar1_val,
      SMA_est = sma1_val
    ))
  }
  return(results)
}

res_sken1 <- run_simulation(n_samples = 36,  missing_pct = 0.0)  # n36, MV 0%
res_sken2 <- run_simulation(n_samples = 120, missing_pct = 0.0)  # n120, MV 0%
res_sken3 <- run_simulation(n_samples = 36,  missing_pct = 0.1)  # n36, MV 10%
res_sken4 <- run_simulation(n_samples = 120, missing_pct = 0.1)  # n120, MV 10%

Menghitung Bias

calc_bias <- function(df, p_true) {
  means <- colMeans(df, na.rm = TRUE)
  return(means - p_true)
}

true_vals <- c(AR, MA, SAR, SMA)
bias_s1 <- calc_bias(res_sken1, true_vals)
bias_s2 <- calc_bias(res_sken2, true_vals)
bias_s3 <- calc_bias(res_sken3, true_vals)
bias_s4 <- calc_bias(res_sken4, true_vals)

summary_table <- data.frame(
  Parameter = c("AR(1)", "MA(1)", "SAR(1)", "SMA(1)"),
  True_Value = true_vals,
  Skenario_1_n36_MV0  = bias_s1,
  Skenario_2_n120_MV0 = bias_s2,
  Skenario_3_n36_MV10 = bias_s3,
  Skenario_4_n120_MV10= bias_s4
)
print(summary_table)

##         Parameter True_Value Skenario_1_n36_MV0 Skenario_2_n120_MV0
## AR_est      AR(1)        0.7        -0.77979389          -0.9986498
## MA_est      MA(1)        0.4        -0.06772435           0.1456429
## SAR_est    SAR(1)        0.8        -1.23916841          -0.8253322
## SMA_est    SMA(1)        0.5        -0.53838954          -1.4174841
##         Skenario_3_n36_MV10 Skenario_4_n120_MV10
## AR_est           -0.7282071          -0.95013020
## MA_est           -0.1768758           0.09202833
## SAR_est          -1.2064704          -0.84985045
## SMA_est          -0.6014135          -1.40152039

plot_data <- rbind(
  data.frame(Value = res_sken1$AR_est, Skenario = "S1: n=36, MV=0%"),
  data.frame(Value = res_sken2$AR_est, Skenario = "S2: n=120, MV=0%"),
  data.frame(Value = res_sken3$AR_est, Skenario = "S3: n=36, MV=10%"),
  data.frame(Value = res_sken4$AR_est, Skenario = "S4: n=120, MV=10%")
)

ggplot(plot_data, aes(x = Skenario, y = Value, fill = Skenario)) +
  geom_boxplot() +
  geom_hline(yintercept = AR, linetype = "dashed", color = "red", linewidth = 1) +
  labs(title = "Distribusi Estimasi Parameter AR(1) lintas 4 Skenario Kombinasi",
       subtitle = "Garis merah putus-putus adalah Parameter Sebenarnya (True Value = 0.7)",
       y = "Nilai Estimasi", x = "Skenario") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))