HANDS ON PROGRES UAS TDS 1 NOVA

# ============================================================
#   BAB 4 — HASIL DAN PEMBAHASAN
#   Pemodelan ARIMA dan SARIMA pada Data IDX Kompas 100
# ============================================================
rm(list = ls())
set.seed(42)

library(tseries)   # Uji ADF, Jarque-Bera

library(forecast)  # ARIMA, SARIMA, auto.arima, forecast()

library(TSA)       # EACF

library(tseries)

# 4.1 Membaca Data
# ============================================================
data_idx <- read.csv(
  "C:/Users/Msi user/Downloads/Tugas Dalam Statistika I/Data Historis IDX Kompas 100 (1).csv",
  sep          = ",",
  header       = TRUE,
  fileEncoding = "UTF-8-BOM"
)

# Bersihkan nama kolom
names(data_idx) <- c("Tanggal", "Terakhir", "Pembukaan",
                     "Tertinggi", "Terendah", "Volume")

# Konversi format tanggal
data_idx$Tanggal <- as.Date(data_idx$Tanggal, format = "%m/%d/%Y")

# Bersihkan kolom numerik — hilangkan titik/koma pemisah ribuan
clean_num <- function(x) {
  x <- gsub("\\.", "", as.character(x))
  x <- gsub(",", ".", x)
  as.numeric(x)
}

data_idx$Terakhir   <- clean_num(data_idx$Terakhir)
data_idx$Pembukaan  <- clean_num(data_idx$Pembukaan)
data_idx$Tertinggi  <- clean_num(data_idx$Tertinggi)
data_idx$Terendah   <- clean_num(data_idx$Terendah)

# Urutkan dari terlama ke terbaru
data_idx <- data_idx[order(data_idx$Tanggal), ]

# Lihat struktur data
str(data_idx)

head(data_idx)

tail(data_idx)

# 4.2 Statistik Deskriptif Harga Penutupan IDX Kompas 100
# ============================================================
cat("======================================================\n")

cat("   STATISTIK DESKRIPTIF HARGA PENUTUPAN (TERAKHIR)\n")

cat("   IDX KOMPAS 100\n")

cat("======================================================\n")

summary(data_idx$Terakhir)

cat(sprintf("\nJumlah observasi  : %d\n",   nrow(data_idx)))

cat(sprintf("Standar deviasi   : %.2f\n",  sd(data_idx$Terakhir, na.rm = TRUE)))

cat(sprintf("Koefisien variasi : %.2f%%\n",
            sd(data_idx$Terakhir, na.rm = TRUE) /
              mean(data_idx$Terakhir, na.rm = TRUE) * 100))

cat(sprintf("Missing value     : %d\n",    sum(is.na(data_idx$Terakhir))))

cat(sprintf("Periode           : %s s.d. %s\n",
            format(min(data_idx$Tanggal), "%d %B %Y"),
            format(max(data_idx$Tanggal), "%d %B %Y")))

# 4.3 Membentuk Objek Time Series
# ============================================================
# frequency = 5 (hari kerja per minggu)
data.ts <- ts(
  data_idx$Terakhir,
  frequency = 5
)

cat("\nPanjang data :", length(data.ts), "observasi\n")

cat("Frekuensi    :", frequency(data.ts), "hari kerja per minggu\n")

# 4.4 Visualisasi Data Time Series
# ============================================================
ts.plot(
  data.ts,
  type = "l",
  ylab = "Indeks IDX Kompas 100",
  xlab = "Waktu (Minggu ke-)",
  col  = "#1565C0",
  lwd  = 1.2
)
title(main = "Time Series Plot Indeks IDX Kompas 100\n(Data Harian Hari Kerja)")
abline(h   = mean(data.ts),
       col = "firebrick",
       lty = 2,
       lwd = 1.5)
legend("topleft",
       legend = c("Indeks IDX Kompas 100", "Rata-rata"),
       col    = c("#1565C0", "firebrick"),
       lty    = c(1, 2), lwd = 2, cex = 0.85)

# 4.5 Visualisasi Seasonal Plot
# ============================================================
seasonplot(
  data.ts,
  year.labels = TRUE,
  main        = "Seasonal Plot Indeks IDX Kompas 100\n(Siklus Mingguan — 5 Hari Kerja)",
  ylab        = "Indeks IDX Kompas 100",
  col         = rainbow(ceiling(length(data.ts) / 5)),
  cex         = 0.7
)

# 4.6 Monthplot (Pola Per Hari Kerja dalam Siklus Mingguan)
# ============================================================
monthplot(
  data.ts,
  ylab   = "Indeks IDX Kompas 100",
  col    = "#1565C0",
  main   = "Monthplot Indeks IDX Kompas 100\nper Hari Kerja dalam Siklus Mingguan",
  xlab   = "Hari ke- (1=Senin, 5=Jumat)"
)

# 4.7 Uji Kehomogenan Ragam (Fligner-Killeen)
# ============================================================
# H0: Ragam data homogen
# H1: Ragam data tidak homogen
tahun <- cycle(data.ts)
fligner.test(data_idx$Terakhir ~ tahun)

# 4.8 Splitting Data (Training 80% — Testing 20%)
# ============================================================
n <- length(data.ts)

train.ts <- subset(data.ts, start = 1,            end = floor(0.8 * n))
test.ts  <- subset(data.ts, start = floor(0.8 * n) + 1)

cat("\nPanjang data training :", length(train.ts), "observasi\n")

cat("Panjang data testing  :", length(test.ts),  "observasi\n")

# 4.9 Visualisasi Data Training dan Testing
# ============================================================
ts.plot(
  train.ts,
  col  = "#1565C0",
  ylab = "Indeks IDX Kompas 100",
  xlab = "Waktu (Minggu ke-)"
)
title(main = "Training Data — Indeks IDX Kompas 100 (80%)")
points(train.ts, pch = 20, col = "#1565C0", cex = 0.3)

ts.plot(
  test.ts,
  col  = "firebrick",
  ylab = "Indeks IDX Kompas 100",
  xlab = "Waktu (Minggu ke-)"
)
title(main = "Testing Data — Indeks IDX Kompas 100 (20%)")
points(test.ts, pch = 20, col = "firebrick", cex = 0.3)

# BAB 4.10 — ARIMA NON-SEASONAL
# ============================================================

# 4.10.1 Plot ACF dan PACF Data Training
acf0 <- acf(
  train.ts,
  main    = "ACF Data Training — IDX Kompas 100",
  lag.max = 40,
  xaxt    = "n",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(train.ts), labels = 0:40)

pacf0 <- pacf(
  train.ts,
  main    = "PACF Data Training — IDX Kompas 100",
  lag.max = 40,
  xaxt    = "n",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(train.ts), labels = 0:40)

# 4.10.2 Differencing Non-Seasonal (d = 1)
# ============================================================
d1 <- diff(train.ts)

ts.plot(
  d1,
  type = "l",
  ylab = "Δ Indeks IDX Kompas 100",
  xlab = "Waktu (Minggu ke-)",
  col  = "#1565C0"
)
title(main = "Plot Differencing Non-Seasonal (d = 1)")

# 4.10.3 Uji ADF — Data Asli dan Setelah Differencing
# ============================================================
cat("======================================================\n")

cat("  UJI ADF — DATA ASLI (TRAINING)\n")

cat("  H0: Data tidak stasioner (ada unit root)\n")

cat("  H1: Data stasioner\n")

cat("======================================================\n")

print(adf.test(train.ts))

cat("\n======================================================\n")

cat("  UJI ADF — SETELAH DIFFERENCING d = 1\n")

cat("======================================================\n")

print(adf.test(d1))

# 4.10.4 Plot ACF dan PACF Setelah Differencing d = 1
# ============================================================
acf1 <- acf(
  d1,
  lag.max = 40,
  xaxt    = "n",
  main    = "ACF Setelah Differencing d = 1",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(d1), labels = 0:40)

pacf1 <- pacf(
  d1,
  main    = "PACF Setelah Differencing d = 1",
  lag.max = 40,
  xaxt    = "n",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(d1), labels = 0:40)

# EACF untuk identifikasi kandidat model ARIMA
eacf(d1)

# 4.10.5 Kandidat dan Estimasi Model ARIMA
# ============================================================
# Berdasarkan ACF/PACF/EACF, kandidat model:
KandidatARIMA <- c(
  "ARIMA(0,1,0)",
  "ARIMA(0,1,1)",
  "ARIMA(1,1,0)",
  "ARIMA(1,1,1)",
  "ARIMA(0,1,2)"
)
KandidatARIMA

# Estimasi
model_a1 <- Arima(train.ts, order = c(0,1,0))
model_a2 <- Arima(train.ts, order = c(0,1,1))
model_a3 <- Arima(train.ts, order = c(1,1,0))
model_a4 <- Arima(train.ts, order = c(1,1,1))
model_a5 <- Arima(train.ts, order = c(0,1,2))

# Fungsi AICc
calc_AICc <- function(model, n) {
  k   <- length(model$coef)
  aic <- AIC(model)
  aic + (2 * k * (k + 1)) / (n - k - 1)
}

n_train <- length(train.ts)

compARIMA <- data.frame(
  "Kandidat Model" = KandidatARIMA,
  "Nilai AIC"  = round(c(model_a1$aic, model_a2$aic, model_a3$aic,
                          model_a4$aic, model_a5$aic), 2),
  "Nilai AICc" = round(c(calc_AICc(model_a1, n_train),
                          calc_AICc(model_a2, n_train),
                          calc_AICc(model_a3, n_train),
                          calc_AICc(model_a4, n_train),
                          calc_AICc(model_a5, n_train)), 2),
  "Nilai BIC"  = round(c(model_a1$bic, model_a2$bic, model_a3$bic,
                          model_a4$bic, model_a5$bic), 2),
  check.names = FALSE
)

cat("\n======================================================\n")

cat("  TABEL PERBANDINGAN KANDIDAT MODEL ARIMA\n")

cat("======================================================\n")

print(compARIMA)

# Model ARIMA terbaik (AIC terkecil)
daftar_model_arima <- list(model_a1, model_a2, model_a3, model_a4, model_a5)
model_arima_best   <- daftar_model_arima[[which.min(compARIMA$`Nilai AIC`)]]

cat("\n======================================================\n")

cat("  MODEL ARIMA TERBAIK\n")

cat("======================================================\n")

summary(model_arima_best)

# 4.10.6 Diagnostik Residual ARIMA — Uji Ljung-Box
# ============================================================
# H0: Tidak terdapat autokorelasi pada residual
# H1: Terdapat autokorelasi pada residual
lags <- c(5, 10, 15, 20, 25, 30)

hasil_ljungbox_arima <- data.frame(
  Lag       = lags,
  Statistik = NA,
  P_Value   = NA
)

for (i in seq_along(lags)) {
  test <- Box.test(
    residuals(model_arima_best),
    lag   = lags[i],
    type  = "Ljung-Box",
    fitdf = length(model_arima_best$coef)
  )
  hasil_ljungbox_arima$Statistik[i] <- round(test$statistic, 4)
  hasil_ljungbox_arima$P_Value[i]   <- round(test$p.value,   4)
}

cat("\n======================================================\n")

cat("  UJI LJUNG-BOX RESIDUAL ARIMA\n")

cat("======================================================\n")

print(hasil_ljungbox_arima)

# BAB 4.11 — ARIMA SEASONAL (SARIMA)
# ============================================================

# 4.11.1 Seasonal Differencing s = 5 (siklus mingguan 5 hari kerja)
D1 <- diff(train.ts, 5)

ts.plot(
  D1,
  type = "l",
  ylab = "Seasonal Differencing",
  xlab = "Waktu (Minggu ke-)",
  col  = "#1565C0"
)
title(main = "Plot Seasonal Differencing (D=1, s=5)")

# 4.11.2 Plot ACF Setelah Seasonal Differencing
# ============================================================
acf2 <- acf(
  D1,
  lag.max = 40,
  xaxt    = "n",
  main    = "ACF Setelah Seasonal Differencing (D=1, s=5)",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(D1), labels = 0:40)

# 4.11.3 Differencing Gabungan (d=1, D=1, s=5)
# ============================================================
d1D1 <- diff(D1)

ts.plot(
  d1D1,
  type = "l",
  ylab = "Differencing Gabungan",
  xlab = "Waktu (Minggu ke-)",
  col  = "#1565C0"
)
title(main = "Plot Differencing Gabungan (d=1, D=1, s=5)")

cat("\n======================================================\n")

cat("  UJI ADF — DIFFERENCING GABUNGAN (d=1, D=1, s=5)\n")

cat("======================================================\n")

print(adf.test(d1D1))

# 4.11.4 Identifikasi Model SARIMA — ACF & PACF Differencing Gabungan
# ============================================================
acf3 <- acf(
  d1D1,
  lag.max = 40,
  xaxt    = "n",
  main    = "ACF Differencing Gabungan (d=1, D=1, s=5)",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(d1D1), labels = 0:40)

pacf3 <- pacf(
  d1D1,
  lag.max = 40,
  xaxt    = "n",
  main    = "PACF Differencing Gabungan (d=1, D=1, s=5)",
  col     = "#1565C0"
)
axis(1, at = 0:40 / frequency(d1D1), labels = 0:40)

# 4.11.5 Identifikasi Komponen Seasonal (Lag Kelipatan s=5)
# ============================================================
acf3_fresh      <- acf(d1D1, lag.max = 40, plot = FALSE)
acf3_fresh$lag  <- acf3_fresh$lag * frequency(d1D1)
acf3.1          <- as.data.frame(cbind(acf3_fresh$acf, acf3_fresh$lag))
names(acf3.1)   <- c("ACF", "Lag")
acf3.2          <- acf3.1[which(acf3.1$Lag %% 5 == 0 & acf3.1$Lag > 0), ]

barplot(
  height    = acf3.2$ACF,
  names.arg = acf3.2$Lag,
  ylab      = "ACF",
  xlab      = "Lag (kelipatan s=5)",
  col       = "#1565C0",
  main      = "ACF pada Lag Musiman (kelipatan s=5)"
)
abline(h = c( 1.96 / sqrt(length(d1D1)),
             -1.96 / sqrt(length(d1D1))),
       col = "red", lty = 2)

pacf3_fresh      <- pacf(d1D1, lag.max = 40, plot = FALSE)
pacf3_fresh$lag  <- pacf3_fresh$lag * frequency(d1D1)
pacf3.1          <- as.data.frame(cbind(pacf3_fresh$acf, pacf3_fresh$lag))
names(pacf3.1)   <- c("PACF", "Lag")
pacf3.2          <- pacf3.1[which(pacf3.1$Lag %% 5 == 0 & pacf3.1$Lag > 0), ]

barplot(
  height    = pacf3.2$PACF,
  names.arg = pacf3.2$Lag,
  ylab      = "PACF",
  xlab      = "Lag (kelipatan s=5)",
  col       = "#1565C0",
  main      = "PACF pada Lag Musiman (kelipatan s=5)"
)
abline(h = c( 1.96 / sqrt(length(d1D1)),
             -1.96 / sqrt(length(d1D1))),
       col = "red", lty = 2)

# EACF
eacf(d1D1)

# 4.12 Kandidat Model SARIMA(p,1,q)(P,1,Q)[5]
# ============================================================
# Berdasarkan ACF, PACF, dan EACF:
#   SARIMA(0,1,0)(0,1,1)[5]
#   SARIMA(0,1,0)(0,1,2)[5]
#   SARIMA(0,1,0)(1,1,1)[5]
#   SARIMA(0,1,1)(0,1,1)[5]
#   SARIMA(0,1,2)(0,1,1)[5]

# 4.13 Estimasi Parameter Model SARIMA
# ============================================================
tmodel1 <- Arima(train.ts, order = c(0,1,0),
                 seasonal = list(order = c(0,1,1), period = 5))
summary(tmodel1)

tmodel2 <- Arima(train.ts, order = c(0,1,0),
                 seasonal = list(order = c(0,1,2), period = 5))
summary(tmodel2)

tmodel3 <- Arima(train.ts, order = c(0,1,0),
                 seasonal = list(order = c(1,1,1), period = 5))
summary(tmodel3)

tmodel4 <- Arima(train.ts, order = c(0,1,1),
                 seasonal = list(order = c(0,1,1), period = 5))
summary(tmodel4)

tmodel5 <- Arima(train.ts, order = c(0,1,2),
                 seasonal = list(order = c(0,1,1), period = 5))
summary(tmodel5)

# 4.14 Pemilihan Model SARIMA Terbaik
# ============================================================
KandidatModelSARIMA <- c(
  "SARIMA(0,1,0)(0,1,1)[5]",
  "SARIMA(0,1,0)(0,1,2)[5]",
  "SARIMA(0,1,0)(1,1,1)[5]",
  "SARIMA(0,1,1)(0,1,1)[5]",
  "SARIMA(0,1,2)(0,1,1)[5]"
)

compmodelSARIMA <- data.frame(
  "Kandidat Model" = KandidatModelSARIMA,
  "Nilai AIC"  = round(c(tmodel1$aic, tmodel2$aic, tmodel3$aic,
                          tmodel4$aic, tmodel5$aic), 2),
  "Nilai AICc" = round(c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc,
                          tmodel4$aicc, tmodel5$aicc), 2),
  "Nilai BIC"  = round(c(tmodel1$bic, tmodel2$bic, tmodel3$bic,
                          tmodel4$bic, tmodel5$bic), 2),
  check.names = FALSE
)

cat("\n======================================================\n")

cat("  TABEL PERBANDINGAN KANDIDAT MODEL SARIMA\n")

cat("======================================================\n")

print(compmodelSARIMA)

# Model terbaik (AIC terkecil)
daftar_sarima  <- list(tmodel1, tmodel2, tmodel3, tmodel4, tmodel5)
idx_terbaik    <- which.min(compmodelSARIMA$`Nilai AIC`)
model_sarima_best <- daftar_sarima[[idx_terbaik]]

cat("\nModel SARIMA terbaik :", KandidatModelSARIMA[idx_terbaik], "\n")

# 4.15 Perbandingan dengan auto.arima
# ============================================================
auto_model <- auto.arima(
  train.ts,
  seasonal      = TRUE,
  stepwise      = TRUE,
  approximation = FALSE,
  D             = 1,
  max.P         = 2,
  max.Q         = 2
)
auto_model

PerbandinganModel <- data.frame(
  Model = c(KandidatModelSARIMA[idx_terbaik], "SARIMA auto.arima"),
  AIC   = c(round(model_sarima_best$aic, 2), round(auto_model$aic, 2))
)
print(PerbandinganModel)

cat(sprintf("\nModel dengan AIC terkecil: %s\n",
            PerbandinganModel$Model[which.min(PerbandinganModel$AIC)]))

# 4.16 Pengujian Signifikansi Parameter Model SARIMA Terbaik
# ============================================================
printstatarima <- function(x, digits = 4, se = TRUE) {
  if (length(x$coef) > 0) {
    cat("\nCoefficients:\n")
    coef <- round(x$coef, digits = digits)
    if (se && nrow(x$var.coef)) {
      ses        <- rep(0, length(coef))
      ses[x$mask]<- round(sqrt(diag(x$var.coef)), digits = digits)
      coef       <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
      coef       <- rbind(coef, s.e. = ses)
      statt      <- coef[1,] / ses
      pval       <- 2 * pt(abs(statt),
                            df         = length(x$residuals) - 1,
                            lower.tail = FALSE)
      coef       <- rbind(coef,
                          t     = round(statt, digits = digits),
                          sign. = round(pval,  digits = digits))
      coef <- t(coef)
    }
    print.default(coef, print.gap = 2)
  }
}

printstatarima(model_sarima_best)

# 4.17 Diagnostik Residual Model SARIMA Terbaik
# ============================================================
png("diagnostik_sarima.png", width = 1200, height = 900)

tsdisplay(
  residuals(model_sarima_best),
  lag.max = 40,
  main = paste("Diagnostik Residual",
               KandidatModelSARIMA[idx_terbaik]),
  col = "#1565C0"
)

dev.off()

# --- Uji Ljung-Box ---
# H0: Tidak terdapat autokorelasi pada residual
# H1: Terdapat autokorelasi pada residual
hasil_ljungbox_sarima <- data.frame(
  Lag       = lags,
  Statistik = NA,
  P_Value   = NA
)
for (i in seq_along(lags)) {
  test <- Box.test(
    residuals(model_sarima_best),
    lag  = lags[i],
    type = "Ljung-Box"
  )
  hasil_ljungbox_sarima$Statistik[i] <- round(test$statistic, 4)
  hasil_ljungbox_sarima$P_Value[i]   <- round(test$p.value,   4)
}
cat("\n======================================================\n")

cat("  UJI LJUNG-BOX RESIDUAL SARIMA\n")

cat("======================================================\n")

print(hasil_ljungbox_sarima)

# --- Uji Normalitas Jarque-Bera ---
# H0: Residual berdistribusi normal
# H1: Residual tidak berdistribusi normal
cat("\n======================================================\n")

cat("  UJI JARQUE-BERA RESIDUAL SARIMA\n")

cat("======================================================\n")

print(jarque.bera.test(residuals(model_sarima_best)))

par(mfrow = c(1,2))
hist(residuals(model_sarima_best),
     main   = paste("Histogram Residual", KandidatModelSARIMA[idx_terbaik]),
     col    = "lightblue", border = "white",
     xlab   = "Residual", breaks = 30)
qqnorm(residuals(model_sarima_best),
       main = paste("QQ-Plot Residual", KandidatModelSARIMA[idx_terbaik]))
qqline(residuals(model_sarima_best), col = "red", lwd = 2)

par(mfrow = c(1,1))

# 4.18 Perbandingan Metrik Evaluasi ARIMA vs SARIMA
# ============================================================
mape_func <- function(actual, predicted) {
  mean(abs((actual - predicted) / actual)) * 100
}

fc_arima  <- forecast(model_arima_best,  h = length(test.ts))
fc_sarima <- forecast(model_sarima_best, h = length(test.ts))

actual <- as.numeric(test.ts)

rmse_a <- round(sqrt(mean((actual - fc_arima$mean)^2)), 4)
mae_a  <- round(mean(abs(actual - fc_arima$mean)), 4)
mape_a <- round(mape_func(actual, as.numeric(fc_arima$mean)), 4)

rmse_s <- round(sqrt(mean((actual - fc_sarima$mean)^2)), 4)
mae_s  <- round(mean(abs(actual - fc_sarima$mean)), 4)
mape_s <- round(mape_func(actual, as.numeric(fc_sarima$mean)), 4)

nama_arima_best  <- KandidatARIMA[which.min(compARIMA$`Nilai AIC`)]
nama_sarima_best <- KandidatModelSARIMA[idx_terbaik]

TabelPerbandingan <- data.frame(
  Model      = c(nama_arima_best, nama_sarima_best),
  RMSE       = c(rmse_a, rmse_s),
  MAE        = c(mae_a,  mae_s),
  "MAPE (%)" = c(mape_a, mape_s),
  AIC        = c(round(AIC(model_arima_best),  2),
                 round(AIC(model_sarima_best), 2)),
  BIC        = c(round(BIC(model_arima_best),  2),
                 round(BIC(model_sarima_best), 2)),
  check.names = FALSE
)

cat("\n======================================================\n")

cat("  PERBANDINGAN METRIK EVALUASI MODEL\n")

cat("======================================================\n")

print(TabelPerbandingan)

best_model <- ifelse(rmse_s < rmse_a, nama_sarima_best, nama_arima_best)
cat(sprintf("\n★ Model terbaik berdasarkan RMSE : %s\n", best_model))

cat("\n--- Akurasi ARIMA ---\n");  accuracy(fc_arima, test.ts)

cat("\n--- Akurasi SARIMA ---\n"); accuracy(fc_sarima, test.ts)

# 4.19 Visualisasi Forecasting — Validasi Data Testing
# ============================================================
# Plot ARIMA
plot(fc_arima, col = "#1565C0",
     main = paste("Forecasting", nama_arima_best, "— Validasi Data Testing"),
     xlab = "Waktu (Minggu ke-)", ylab = "Indeks IDX Kompas 100",
     shadecols = c("mistyrose", "lightyellow"))
lines(test.ts, col = "darkgreen", lwd = 2)
legend("topleft",
       legend = c("Training", paste("Forecast", nama_arima_best), "Aktual Testing"),
       col    = c("#1565C0", "blue", "darkgreen"),
       lty = 1, lwd = 2, cex = 0.8)

# Plot SARIMA
plot(fc_sarima, col = "#1565C0",
     main = paste("Forecasting", nama_sarima_best, "— Validasi Data Testing"),
     xlab = "Waktu (Minggu ke-)", ylab = "Indeks IDX Kompas 100",
     shadecols = c("thistle1", "lavender"))
lines(test.ts, col = "darkgreen", lwd = 2)
legend("topleft",
       legend = c("Training", paste("Forecast", nama_sarima_best), "Aktual Testing"),
       col    = c("#1565C0", "purple", "darkgreen"),
       lty = 1, lwd = 2, cex = 0.8)

# Overlay Aktual vs ARIMA vs SARIMA
ts.plot(
  test.ts, fc_arima$mean, fc_sarima$mean,
  gpars = list(
    col  = c("darkgreen", "red", "purple"),
    lty  = c(1, 2, 3), lwd = c(2, 1.5, 1.5),
    main = "Overlay: Aktual vs ARIMA vs SARIMA (Periode Testing)\nIndeks IDX Kompas 100",
    ylab = "Indeks IDX Kompas 100", xlab = "Waktu (Minggu ke-)"
  )
)
legend("topleft",
       legend = c("Aktual",
                  sprintf("%s (MAPE = %.2f%%)", nama_arima_best,  mape_a),
                  sprintf("%s (MAPE = %.2f%%)", nama_sarima_best, mape_s)),
       col = c("darkgreen", "red", "purple"),
       lty = c(1, 2, 3), lwd = 2, cex = 0.8)

# 4.20 Tabel dan Visualisasi Hasil Forecasting Model Terbaik
# ============================================================
forecast_terbaik <- cbind(fc_sarima$mean, fc_sarima$lower, fc_sarima$upper)
forecast_terbaik <- as.data.frame(forecast_terbaik)
colnames(forecast_terbaik) <- c("Forecast", "Lo 80", "Lo 95", "Hi 80", "Hi 95")

cat("\nTabel Hasil Forecasting", nama_sarima_best, "(20 baris pertama):\n")

print(head(round(forecast_terbaik, 2), 20))

# Visualisasi
plot(fc_sarima$mean, type = "l", col = "#1565C0", lwd = 2,
     xlab = "Periode", ylab = "Indeks IDX Kompas 100",
     main = paste("Hasil Forecasting", nama_sarima_best,
                  "\npada Data Testing"))
lines(fc_sarima$lower[,2], col = "firebrick", lwd = 1.5, lty = 2)
lines(fc_sarima$upper[,2], col = "darkgreen", lwd = 1.5, lty = 2)
lines(as.numeric(test.ts),  col = "black",    lwd = 1.5)
legend("topleft",
       legend = c("Forecast", "Lower 95%", "Upper 95%", "Aktual"),
       col    = c("#1565C0", "firebrick", "darkgreen", "black"),
       lty    = c(1,2,2,1), lwd = 2, cex = 0.8)

# 4.21 Forecast 1 Minggu ke Depan (5 Hari Kerja)
# ============================================================
fc_future <- forecast(model_sarima_best, h = 5)

plot(
  fc_future,
  col       = "#1565C0",
  fcol      = "purple",
  main      = paste("Peramalan Indeks IDX Kompas 100\n5 Hari Kerja ke Depan —",
                    nama_sarima_best),
  xlab      = "Waktu (Minggu ke-)",
  ylab      = "Indeks IDX Kompas 100",
  shadecols = c("thistle1", "lavender")
)

cat("\nForecast 5 Hari Kerja ke Depan:\n")

print(round(as.data.frame(fc_future), 2))

# 4.22 KAJIAN SIMULASI
# ============================================================

# 4.22.1 Simulasi Data ARIMA(1,1,1)
cat("\n======================================================\n")

cat("  KAJIAN SIMULASI — ARIMA(1,1,1)\n")

cat("======================================================\n")

n      <- 200
phi1   <-  0.5; phi2   <- -0.7
theta1 <- -0.3; theta2 <- -0.5
sd_e   <- sqrt(1.20)

e <- rnorm(n, mean = 0, sd = sd_e)
x <- numeric(n)
for (t in 3:n) {
  x[t] <- phi1*x[t-1] + phi2*x[t-2] +
          e[t] + theta1*e[t-1] + theta2*e[t-2]
}
x_d1     <- cumsum(x)
data_sim <- cumsum(x_d1)

# Hapus 50 data awal (burn-in)
data_trim  <- data_sim[51:200]
n_data     <- length(data_trim)
train_size <- floor(0.8 * n_data)

train_sim <- ts(data_trim[1:train_size])
test_sim  <- ts(data_trim[(train_size+1):n_data], start = train_size + 1)

ts.plot(train_sim, main = "Data Training — Simulasi ARIMA(1,1,1)",
        ylab = "Nilai", xlab = "Waktu", col = "#1565C0")

# ADF
cat("ADF Data Asli Simulasi:\n");   print(adf.test(train_sim))

train_d1 <- diff(train_sim);        cat("ADF d=1:\n"); print(adf.test(train_d1))

train_d2 <- diff(train_d1);         cat("ADF d=2:\n"); print(adf.test(train_d2))

par(mfrow = c(1,2))
acf(train_d2,  main = "ACF Differencing ke-2 — Simulasi ARIMA")
pacf(train_d2, main = "PACF Differencing ke-2 — Simulasi ARIMA")

par(mfrow = c(1,1))
eacf(as.matrix(train_d2))

# Fitting
model_sim <- forecast::Arima(train_sim, order = c(1,1,1))
cat("\nKoefisien ARIMA(1,1,1):\n"); print(model_sim$coef)

result_sim <- data.frame(
  Model = "ARIMA(1,1,1)",
  AIC   = round(AIC(model_sim), 4),
  BIC   = round(BIC(model_sim), 4),
  AICc  = round(calc_AICc(model_sim, length(train_sim)), 4)
)
print(result_sim)

# Diagnostik
ts.plot(residuals(model_sim), main = "Residual ARIMA(1,1,1) — Simulasi",
        ylab = "Residual", col = "purple")
abline(h = 0, col = "red", lty = 2)

acf(residuals(model_sim), main = "ACF Residual — Simulasi ARIMA")

cat("\nLjung-Box:\n");   print(Box.test(residuals(model_sim), lag=20, type="Ljung-Box"))

cat("\nShapiro-Wilk:\n"); print(shapiro.test(residuals(model_sim)))

par(mfrow = c(1,2))
hist(residuals(model_sim), main = "Histogram Residual Simulasi ARIMA",
     col = "lightblue", border = "black")
qqnorm(residuals(model_sim), main = "QQ Plot Residual Simulasi ARIMA")
qqline(residuals(model_sim), col = "red", lwd = 2)

par(mfrow = c(1,1))

# Perbandingan parameter asli vs estimasi
tabel_param_arima <- data.frame(
  Parameter = c("phi1", "phi2", "theta1", "theta2"),
  Asli      = c(phi1, phi2, theta1, theta2),
  Estimasi  = round(as.numeric(model_sim$coef[1:4]), 4),
  Selisih   = round(abs(c(phi1,phi2,theta1,theta2) -
                        as.numeric(model_sim$coef[1:4])), 4)
)
cat("\n======================================================\n")

cat("  PARAMETER ASLI vs ESTIMASI — SIMULASI ARIMA\n")

cat("======================================================\n")

print(tabel_param_arima)

# ============================================================
# 4.22.2 Simulasi Data SARIMA(1,1,1)(1,1,1)[5]
# ============================================================

sim_sarima_data <- arima.sim(
  model = list(
    order = c(1,1,1),
    seasonal = list(order = c(1,1,1), period = 5),
    ar  = 0.5,
    ma  = -0.3,
    sar = 0.4,
    sma = -0.5
  ),
  n = 300
)

# Hapus 50 data awal (burn-in)
sim_sar_trim <- sim_sarima_data[51:300]

# Training-testing
n_sar <- length(sim_sar_trim)
train_size_sar <- floor(0.8 * n_sar)

train_sar <- ts(
  sim_sar_trim[1:train_size_sar],
  frequency = 5
)

test_sar <- ts(
  sim_sar_trim[(train_size_sar + 1):n_sar],
  frequency = 5
)

# Fitting model SARIMA simulasi
model_sar_sim <- Arima(
  train_sar,
  order = c(1,1,1),
  seasonal = list(order = c(1,1,1), period = 5)
)

summary(model_sar_sim)

# Parameter asli
phi1_s   <- 0.5
theta1_s <- -0.3
Phi1_s   <- 0.4
Theta1_s <- -0.5

# Differencing gabungan
D1_sar   <- diff(train_sar, lag = 5)
d1D1_sar <- diff(D1_sar)

ts.plot(d1D1_sar,
        main = "Differencing Gabungan (d=1,D=1) — Simulasi SARIMA",
        ylab = "Δ Nilai", xlab = "Waktu", col = "#1A237E")
abline(h = 0, col = "red", lty = 2)

cat("ADF Differencing Gabungan — Simulasi SARIMA:\n")

print(adf.test(d1D1_sar))

par(mfrow = c(1,2))
acf(d1D1_sar,  lag.max = 40, main = "ACF Diff. Gabungan — Simulasi SARIMA",
    col = "#1A237E", xaxt = "n")
axis(1, at = 0:40/5, labels = 0:40)
pacf(d1D1_sar, lag.max = 40, main = "PACF Diff. Gabungan — Simulasi SARIMA",
     col = "#1A237E", xaxt = "n")
axis(1, at = 0:40/5, labels = 0:40)

par(mfrow = c(1,1))

# Fitting
png("Diagnostik_Residual_SARIMA.png",
    width = 1600,
    height = 1200)

tsdisplay(
  residuals(model_sar_sim),
  lag.max = 40,
  main = "Diagnostik Residual SARIMA(1,1,1)(1,1,1)[5] - Simulasi",
  col = "#00796B"
)

dev.off()

# Uji Ljung-Box
cat("\nLjung-Box SARIMA Simulasi:\n")

print(Box.test(
  residuals(model_sar_sim),
  lag = 20,
  type = "Ljung-Box"
))

# Uji Jarque-Bera
cat("\nJarque-Bera SARIMA Simulasi:\n")

print(jarque.bera.test(residuals(model_sar_sim)))

tabel_param_sarima <- data.frame(
  Parameter = c("phi1", "theta1", "Phi1", "Theta1"),
  Asli      = c(phi1_s, theta1_s, Phi1_s, Theta1_s),
  Estimasi  = round(c(
    model_sar_sim$coef["ar1"],
    model_sar_sim$coef["ma1"],
    model_sar_sim$coef["sar1"],
    model_sar_sim$coef["sma1"]
  ), 4),
  Selisih = round(abs(
    c(phi1_s, theta1_s, Phi1_s, Theta1_s) -
    c(
      model_sar_sim$coef["ar1"],
      model_sar_sim$coef["ma1"],
      model_sar_sim$coef["sar1"],
      model_sar_sim$coef["sma1"]
    )
  ), 4)
)

#   SELESAI — BAB 4
# ============================================================
cat("\n============================================================\n")

cat("  BAB 4 SELESAI DIEKSEKUSI\n")

cat("  Model Terbaik ARIMA  :", nama_arima_best,  "\n")

cat("  Model Terbaik SARIMA :", nama_sarima_best, "\n")

cat("  MAPE ARIMA           :", mape_a, "%\n")

cat("  MAPE SARIMA          :", mape_s, "%\n")

cat("============================================================\n")