# -----------------------------------------------------------
# Install paket yang diperlukan (jika belum terpasang)
# -----------------------------------------------------------
pkgs <- c("readxl","forecast","ggplot2","dplyr","tidyr",
          "zoo","knitr","kableExtra","scales","patchwork")
new_pkgs <- pkgs[!sapply(pkgs, requireNamespace, quietly=TRUE)]
if (length(new_pkgs)) install.packages(new_pkgs)

library(readxl);   library(forecast); library(ggplot2)
library(dplyr);    library(tidyr);    library(zoo)
library(knitr);    library(kableExtra); library(scales)
library(patchwork)

theme_set(theme_bw(base_size = 11))

1 Persiapan Data

# ================================================================
#  SESUAIKAN PATH FILE EXCEL DI BAWAH INI
# ================================================================
FILE_PATH <- "C:/Users/ADVAN/Downloads/harga_cpo (2).xlsx"   # letakkan file xlsx di working directory
SHEET_NUM <- 3                  # "Rata-rata Bulanan"
N_TEST    <- 20                 # jumlah bulan data uji

# ----------------------------------------------------------------
# Fungsi parse angka format Indonesia
# (titik = pemisah ribuan, koma = desimal)
# ----------------------------------------------------------------
parse_indo <- function(x) {
  if (is.character(x)) {
    x <- trimws(x)
    x <- gsub("\\.", "",  x)  # hapus pemisah ribuan
    x <- gsub(",",   ".", x)  # ganti koma desimal → titik
  }
  as.numeric(x)
}

# ----------------------------------------------------------------
# Fungsi tambah 1 bulan ke objek Date (tanpa lubridate)
# ----------------------------------------------------------------
add_month <- function(d, n = 1) {
  m <- as.integer(format(d, "%m")) + n
  y <- as.integer(format(d, "%Y"))
  while (m > 12) { m <- m - 12; y <- y + 1 }
  while (m < 1)  { m <- m + 12; y <- y - 1 }
  as.Date(sprintf("%04d-%02d-01", y, m))
}

# ----------------------------------------------------------------
# Muat data
# ----------------------------------------------------------------
df_raw <- read_excel(FILE_PATH, sheet = SHEET_NUM)
colnames(df_raw) <- c("Bulan", "Rotterdam", "MDEX", "ICDX")
df_raw <- df_raw %>% mutate(across(c(Rotterdam, MDEX, ICDX), parse_indo))

# Data di sheet urut terbaru → terlama, balik ke kronologis
df <- df_raw[nrow(df_raw):1, ]

# Parse nama bulan Indonesia → Date
bln_map <- c(
  Januari="01",Februari="02",Maret="03",April="04",
  Mei="05",Juni="06",Juli="07",Agustus="08",
  September="09",Oktober="10",November="11",Desember="12"
)
parse_bulan <- function(s) {
  parts <- strsplit(trimws(s), " ")[[1]]
  as.Date(sprintf("%s-%s-01", parts[2], bln_map[parts[1]]))
}
df$Tanggal <- as.Date(sapply(df$Bulan, parse_bulan))
df <- df[, c("Tanggal","Rotterdam","MDEX","ICDX")]

cat("Rentang data:", format(min(df$Tanggal),"%B %Y"),
    "hingga", format(max(df$Tanggal),"%B %Y"), "\n")

## Rentang data: December 2019 hingga May 2026

cat("Total observasi:", nrow(df), "bulan\n")

## Total observasi: 78 bulan

# ----------------------------------------------------------------
# Penanganan outlier Rotterdam (nilai ~6909 pada Jun 2024)
# ----------------------------------------------------------------
out_idx <- which(df$Rotterdam > 5000)
if (length(out_idx) > 0) {
  cat("\n⚠️  Outlier terdeteksi – Rotterdam",
      format(df$Tanggal[out_idx], "%B %Y"),
      "| Nilai asli:", round(df$Rotterdam[out_idx], 2),
      "→ diganti interpolasi linear\n")
  df$Rotterdam[out_idx] <- NA
  df$Rotterdam <- as.numeric(na.approx(df$Rotterdam, na.rm = FALSE))
}

## 
## ⚠️  Outlier terdeteksi – Rotterdam June 2024 | Nilai asli: 10154375 → diganti interpolasi linear

# ----------------------------------------------------------------
# Split train / test
# ----------------------------------------------------------------
n_total <- nrow(df)
n_train <- n_total - N_TEST
df_train <- df[1:n_train, ]
df_test  <- df[(n_train + 1):n_total, ]

cat("\nData latih:", format(min(df_train$Tanggal),"%B %Y"),
    "–", format(max(df_train$Tanggal),"%B %Y"),
    "(", n_train, "bulan)\n")

## 
## Data latih: December 2019 – September 2024 ( 58 bulan)

cat("Data uji :", format(min(df_test$Tanggal),"%B %Y"),
    "–", format(max(df_test$Tanggal),"%B %Y"),
    "(", N_TEST,  "bulan)\n")

## Data uji : October 2024 – May 2026 ( 20 bulan)

df %>% head(5) %>%
  kbl(caption = "5 Observasi Pertama (Desember 2019 – ...)") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

5 Observasi Pertama (Desember 2019 – …)
Tanggal	Rotterdam	MDEX	ICDX
2019-12-01	665.00	602.74	597.58
2020-01-01	735.88	652.47	640.86
2020-02-01	840.62	740.26	771.11
2020-03-01	761.47	676.66	717.23
2020-04-01	644.26	569.75	590.41

2 Eksplorasi Data

2.1 Plot Deret Waktu

df_long <- df %>%
  pivot_longer(c(Rotterdam,MDEX,ICDX), names_to="Pasar", values_to="Harga") %>%
  mutate(Pasar = factor(Pasar, levels=c("Rotterdam","MDEX","ICDX")))

label_pasar <- c(
  Rotterdam = "Rotterdam (EU CIF, USD/ton)",
  MDEX      = "MDEX Malaysia (USD/ton)",
  ICDX      = "ICDX Indonesia (USD/ton)"
)

ggplot(df_long, aes(x=Tanggal, y=Harga, color=Pasar)) +
  geom_line(linewidth = 0.85) +
  geom_point(size = 1.1) +
  geom_vline(xintercept = as.numeric(df_test$Tanggal[1]),
             linetype="dashed", color="firebrick", linewidth=0.7) +
  annotate("text", x=df_test$Tanggal[1], y=Inf,
           label=" Data Uji →", hjust=0, vjust=1.5, color="firebrick", size=3.5) +
  facet_wrap(~Pasar, ncol=1, scales="free_y",
             labeller=labeller(Pasar=label_pasar)) +
  scale_color_manual(values=c(Rotterdam="#1f77b4",MDEX="#d62728",ICDX="#2ca02c")) +
  labs(title="Harga Rata-rata Bulanan CPO – Tiga Pasar Internasional",
       x=NULL, y="Harga (USD/ton)") +
  theme(legend.position="none", strip.text=element_text(face="bold"))

2.2 Plot ACF – Data Asli

par(mfrow=c(3,1), mar=c(3,4,3,1))
for (nm in c("Rotterdam","MDEX","ICDX")) {
  acf(df_train[[nm]], lag.max=20, main=paste("ACF –", nm, "(Data Latih)"),
      col="steelblue", lwd=2)
}

par(mfrow=c(1,1))

3 Fungsi Bantu

# ----------------------------------------------------------------
# Metrik kesalahan
# ----------------------------------------------------------------
calc_metrics <- function(actual, predicted) {
  valid <- !is.na(actual) & !is.na(predicted)
  e     <- actual[valid] - predicted[valid]
  c(RMSE = sqrt(mean(e^2)),
    MSE  = mean(e^2),
    MAPE = mean(abs(e / actual[valid])) * 100)
}

# ----------------------------------------------------------------
# Uji Ljung-Box (H0: tidak ada autokorelasi)
# Lolos jika p-value > 0.05
# ----------------------------------------------------------------
uji_lb <- function(resid, model_name = "", lag = NULL) {
  r <- na.omit(resid)
  if (is.null(lag)) lag <- min(10, max(1, floor(length(r)/5)))
  if (length(r) <= lag) lag <- max(1, length(r) - 1)
  bt    <- Box.test(r, lag=lag, type="Ljung-Box")
  lolos <- bt$p.value > 0.05
  cat(sprintf("  %-35s | p = %.4f | %s\n",
              model_name, bt$p.value,
              ifelse(lolos, "✅ Tidak ada autokorelasi",
                            "❌ Ada autokorelasi (residual tidak acak)")))
  list(p_value = bt$p.value, lolos = lolos)
}

# ----------------------------------------------------------------
# Double Moving Average
# ----------------------------------------------------------------
dma_compute <- function(y, m) {
  ma1 <- as.numeric(stats::filter(y, rep(1/m, m), sides=1))
  ma2 <- as.numeric(stats::filter(ma1, rep(1/m, m), sides=1))
  list(ma1=ma1, ma2=ma2, m=m, y=y)
}

dma_forecast_h <- function(obj, h) {
  valid <- which(!is.na(obj$ma1) & !is.na(obj$ma2))
  if (!length(valid)) return(rep(NA, h))
  T_v <- tail(valid, 1)
  a   <- 2*obj$ma1[T_v] - obj$ma2[T_v]
  b   <- (2/(obj$m - 1)) * (obj$ma1[T_v] - obj$ma2[T_v])
  a + b * seq_len(h)
}

dma_residuals <- function(y, m) {
  obj <- dma_compute(y, m)
  n   <- length(y)
  # Fitted 1-step: a_{t-1} + b_{t-1}
  a_vec <- 2*obj$ma1 - obj$ma2
  b_vec <- (2/(m-1)) * (obj$ma1 - obj$ma2)
  fitted_vals <- c(NA, head(a_vec + b_vec, n - 1))
  y - fitted_vals
}

# ----------------------------------------------------------------
# Rolling one-step-ahead forecast (level)
# fn_pred(y_available) → angka prediksi satu langkah
# ----------------------------------------------------------------
rolling_os <- function(y_train, y_test, fn_pred) {
  n <- length(y_test)
  p <- numeric(n)
  for (i in seq_len(n)) {
    y_now <- c(y_train, if (i > 1) y_test[1:(i-1)] else NULL)
    p[i]  <- tryCatch(fn_pred(y_now), error=function(e) NA_real_)
  }
  p
}

# ----------------------------------------------------------------
# ts() dari vektor + tanggal awal
# ----------------------------------------------------------------
make_ts <- function(y, d_start) {
  ts(y,
     start     = c(as.integer(format(d_start,"%Y")),
                   as.integer(format(d_start,"%m"))),
     frequency = 12)
}

4 Pemodelan pada Data Asli

Model yang dibangun di bagian ini langsung pada level harga (tanpa diferensiasi):

Naive Simple : \(\hat{Y}_{t+h} = Y_t\)
Naive Trend : \(\hat{Y}_{t+1} = 2Y_t - Y_{t-1}\)
Naive Drift : \(\hat{Y}_{t+h} = Y_t + h \cdot \dfrac{Y_t - Y_1}{t-1}\)
Holt’s Method: grid search \(\alpha \in (0,1)\), \(\beta \in (0,1)\), pilih AIC terkecil
Double MA : grid search \(m \in \{2,\ldots,\lfloor n/4\rfloor\}\), pilih RMSE latih terkecil

MARKETS <- c("Rotterdam","MDEX","ICDX")
results <- list()   # kumpulkan semua hasil

for (mkt in MARKETS) {
  cat("\n\n## Pasar:", mkt, "{.tabset}\n\n")

  y_train  <- df_train[[mkt]]
  y_test   <- df_test[[mkt]]
  n_tr     <- length(y_train)
  ts_train <- make_ts(y_train, df_train$Tanggal[1])
  res_mkt  <- list()

  cat("### Uji Diagnostik\n\n")
  cat("**Uji Ljung-Box pada residual fitted values (p > 0,05 = lolos)**\n\n```\n")

  # ============================================================
  ## Naive Simple
  # ============================================================
  ns_fit   <- naive(ts_train, h=N_TEST)
  ns_resid <- residuals(ns_fit)
  ns_diag  <- uji_lb(ns_resid, "Naive Simple")
  ns_os    <- rolling_os(y_train, y_test, function(y) tail(y,1))

  res_mkt[["Naive Simple"]] <- list(
    lolos    = ns_diag$lolos,
    resid    = as.numeric(ns_resid),
    fcast_ms = as.numeric(ns_fit$mean),
    fcast_os = ns_os
  )

  # ============================================================
  ## Naive Trend
  # ============================================================
  nt_fitted <- c(NA, NA, y_train[1:(n_tr-2)] + diff(y_train, differences=1)[1:(n_tr-2)])
  # Koreksi: Y_hat_{t} = 2*Y_{t-1} - Y_{t-2}
  nt_fitted2 <- c(NA, NA, 2*y_train[2:(n_tr-1)] - y_train[1:(n_tr-2)])
  nt_resid   <- y_train - nt_fitted2
  nt_diag    <- uji_lb(nt_resid, "Naive Trend")
  slope_tr   <- y_train[n_tr] - y_train[n_tr-1]
  nt_fcast_ms <- y_train[n_tr] + slope_tr * seq_len(N_TEST)
  nt_os       <- rolling_os(y_train, y_test, function(y) {
    n <- length(y)
    if (n < 2) y[n] else 2*y[n] - y[n-1]
  })

  res_mkt[["Naive Trend"]] <- list(
    lolos    = nt_diag$lolos,
    resid    = nt_resid,
    fcast_ms = nt_fcast_ms,
    fcast_os = nt_os
  )

  # ============================================================
  ## Naive Drift
  # ============================================================
  nd_fit   <- rwf(ts_train, h=N_TEST, drift=TRUE)
  nd_resid <- residuals(nd_fit)
  nd_diag  <- uji_lb(nd_resid, "Naive Drift")
  nd_os    <- rolling_os(y_train, y_test, function(y) {
    n <- length(y)
    if (n < 2) y[n] else y[n] + (y[n]-y[1])/(n-1)
  })

  res_mkt[["Naive Drift"]] <- list(
    lolos    = nd_diag$lolos,
    resid    = as.numeric(nd_resid),
    fcast_ms = as.numeric(nd_fit$mean),
    fcast_os = nd_os
  )

  # ============================================================
  ## Holt's Method – Grid Search (alpha x beta, step 0.05)
  # ============================================================
  # ============================================================
  # Holt's Method – Grid Search (REVISI AMAN)
  # ============================================================
  alpha_grid <- seq(0.05, 0.95, by=0.05)
  beta_grid  <- seq(0.05, 0.95, by=0.05)
  best_holt  <- list(aic=Inf, alpha=NA, beta=NA, fit=NULL)

  for (al in alpha_grid) {
    for (bt in beta_grid) {
      # Gunakan tryCatch untuk menangkap error numerik
      tryCatch({
        hf      <- holt(ts_train, alpha=al, beta=bt, h=N_TEST, initial="simple")
        aic_val <- hf$model$aic
        
        # Penanganan NA: Pastikan aic_val ada dan tidak NA sebelum dibandingkan
        if (!is.null(aic_val) && !is.na(aic_val)) {
          if (aic_val < best_holt$aic) {
            best_holt <- list(aic=aic_val, alpha=al, beta=bt, fit=hf)
          }
        }
      }, error = function(e) { NULL }) # Abaikan kombinasi yang error
    }
  }

  # Fallback: Jika grid search gagal total, gunakan optimasi otomatis dari fungsi holt
  if (is.null(best_holt$fit)) {
    best_holt$fit   <- holt(ts_train, h=N_TEST)
    best_holt$alpha <- best_holt$fit$model$par["alpha"]
    best_holt$beta  <- best_holt$fit$model$par["beta"]
  }

  # Jalankan uji diagnostik hanya jika model berhasil terbentuk
  res_resid_holt <- residuals(best_holt$fit)
  uji_lb(res_resid_holt, 
         sprintf("Holt (α=%.2f, β=%.2f)", best_holt$alpha, best_holt$beta))
  holt_diag  <- uji_lb(residuals(best_holt$fit), "")   # simpan hasil
  # Ulang agar lebih clean: re-run diam-diam
  {
    r_holt <- na.omit(residuals(best_holt$fit))
    lag_h  <- min(10, max(1, floor(length(r_holt)/5)))
    bt_h   <- Box.test(r_holt, lag=lag_h, type="Ljung-Box")
    holt_diag <- list(p_value=bt_h$p.value, lolos=bt_h$p.value>0.05)
  }

  holt_os <- rolling_os(y_train, y_test, function(y) {
    ts_y <- make_ts(y, df$Tanggal[1])
    tryCatch(
      as.numeric(holt(ts_y, alpha=best_holt$alpha, beta=best_holt$beta,
                      h=1, initial="simple")$mean),
      error=function(e) tail(y,1)
    )
  })

  res_mkt[["Holt"]] <- list(
    lolos    = holt_diag$lolos,
    resid    = as.numeric(residuals(best_holt$fit)),
    fcast_ms = as.numeric(best_holt$fit$mean),
    fcast_os = holt_os,
    alpha    = best_holt$alpha,
    beta     = best_holt$beta
  )

  # ============================================================
  ## Double Moving Average – Grid Search m
  # ============================================================
  m_grid   <- seq(2, min(floor(n_tr/4), 12))
  best_dma <- list(rmse=Inf, m=NA)

  for (m in m_grid) {
    resid_d <- dma_residuals(y_train, m)
    rmse_d  <- sqrt(mean(resid_d^2, na.rm=TRUE))
    if (!is.na(rmse_d) && rmse_d < best_dma$rmse) {
      best_dma <- list(rmse=rmse_d, m=m)
    }
  }

  dma_res_best <- dma_residuals(y_train, best_dma$m)
  uji_lb(dma_res_best, sprintf("DMA (m=%d)", best_dma$m))
  {
    r_dma  <- na.omit(dma_res_best)
    lag_d  <- min(10, max(1, floor(length(r_dma)/5)))
    bt_d   <- Box.test(r_dma, lag=lag_d, type="Ljung-Box")
    dma_diag <- list(p_value=bt_d$p.value, lolos=bt_d$p.value>0.05)
  }

  dma_obj     <- dma_compute(y_train, best_dma$m)
  dma_fcast_ms <- dma_forecast_h(dma_obj, N_TEST)
  dma_os       <- rolling_os(y_train, y_test, function(y) {
    dma_forecast_h(dma_compute(y, best_dma$m), 1)[1]
  })

  res_mkt[["DMA"]] <- list(
    lolos    = dma_diag$lolos,
    resid    = dma_res_best,
    fcast_ms = dma_fcast_ms,
    fcast_os = dma_os,
    m        = best_dma$m
  )

  cat("```\n")
  results[[mkt]] <- res_mkt
}

4.1 Pasar: Rotterdam

4.1.1 Uji Diagnostik

Uji Ljung-Box pada residual fitted values (p > 0,05 = lolos)

  Naive Simple                        | p = 0.7185 | ✅ Tidak ada autokorelasi
  Naive Trend                         | p = 0.0013 | ❌ Ada autokorelasi (residual tidak acak)
  Naive Drift                         | p = 0.7185 | ✅ Tidak ada autokorelasi
  Holt (α=0.82, β=0.00)             | p = 0.8859 | ✅ Tidak ada autokorelasi
                                      | p = 0.8859 | ✅ Tidak ada autokorelasi
  DMA (m=3)                           | p = 0.0413 | ❌ Ada autokorelasi (residual tidak acak)

4.2 Pasar: MDEX

4.2.1 Uji Diagnostik

Uji Ljung-Box pada residual fitted values (p > 0,05 = lolos)

  Naive Simple                        | p = 0.3685 | ✅ Tidak ada autokorelasi
  Naive Trend                         | p = 0.0125 | ❌ Ada autokorelasi (residual tidak acak)
  Naive Drift                         | p = 0.3685 | ✅ Tidak ada autokorelasi
  Holt (α=0.99, β=0.00)             | p = 0.3261 | ✅ Tidak ada autokorelasi
                                      | p = 0.3261 | ✅ Tidak ada autokorelasi
  DMA (m=2)                           | p = 0.0204 | ❌ Ada autokorelasi (residual tidak acak)

4.3 Pasar: ICDX

4.3.1 Uji Diagnostik

Uji Ljung-Box pada residual fitted values (p > 0,05 = lolos)

  Naive Simple                        | p = 0.7123 | ✅ Tidak ada autokorelasi
  Naive Trend                         | p = 0.0118 | ❌ Ada autokorelasi (residual tidak acak)
  Naive Drift                         | p = 0.7123 | ✅ Tidak ada autokorelasi
  Holt (α=1.00, β=0.00)             | p = 0.7694 | ✅ Tidak ada autokorelasi
                                      | p = 0.7694 | ✅ Tidak ada autokorelasi
  DMA (m=2)                           | p = 0.1241 | ✅ Tidak ada autokorelasi

5 Diferensiasi dan Pemodelan Data Diferensial

Data didiferensikan orde-1 (\(\Delta Y_t = Y_t - Y_{t-1}\)) sebelum pemodelan.

Model yang dibangun:

Simple Average: \(\hat{\Delta Y} = \bar{\Delta Y}\) (rata-rata seluruh perbedaan)
Moving Average: grid search \(m \in \{2,\ldots,\lfloor n/3\rfloor\}\) pada data diferensial
SES: grid search \(\alpha \in (0.05, 0.95)\) pada data diferensial

Ramalan level diperoleh dengan rekonstruksi: \(\hat{Y}_{T+h} = Y_T + \sum_{i=1}^{h} \hat{\Delta Y}_{T+i}\)

5.1 Plot Diferensial dan ACF

par(mfrow=c(3,2), mar=c(3,4,3,1))
for (nm in c("Rotterdam","MDEX","ICDX")) {
  dy <- diff(df_train[[nm]])
  plot(dy, type="l", main=paste("Δ Harga –", nm), ylab="Δ (USD/ton)",
       col="steelblue", lwd=1.2)
  abline(h=0, col="red", lty=2)
  acf(dy, lag.max=20, main=paste("ACF Δ –", nm), col="steelblue", lwd=2)
}

par(mfrow=c(1,1))

5.2 Pemodelan Data Diferensial

for (mkt in MARKETS) {
  cat("\n\n## Pasar:", mkt, "(Diferensial) {.tabset}\n\n")
  cat("### Uji Diagnostik\n\n")
  cat("**Uji Ljung-Box (p > 0,05 = lolos)**\n\n```\n")

  y_train  <- df_train[[mkt]]
  y_test   <- df_test[[mkt]]
  n_tr     <- length(y_train)
  dy_train <- diff(y_train)
  n_dtr    <- length(dy_train)
  y_last   <- tail(y_train, 1)

  # ============================================================
  ## Simple Average (pada data diferensial)
  # ============================================================
  sa_mean  <- mean(dy_train)
  sa_resid <- dy_train - sa_mean          # residual = Δy - mean(Δy)
  sa_diag  <- uji_lb(sa_resid, "Simple Average (diff)")
  sa_fcast_ms <- y_last + cumsum(rep(sa_mean, N_TEST))
  sa_os <- rolling_os(y_train, y_test, function(y) {
    tail(y,1) + mean(diff(y))
  })

  results[[mkt]][["SA_diff"]] <- list(
    lolos    = sa_diag$lolos,
    resid    = sa_resid,
    fcast_ms = sa_fcast_ms,
    fcast_os = sa_os
  )

  # ============================================================
  ## Moving Average pada data diferensial – Grid Search m
  # ============================================================
  m_grid2   <- seq(2, min(floor(n_dtr/3), 12))
  best_ma_d <- list(rmse=Inf, m=NA)

  for (m in m_grid2) {
    # Fitted: rollmean, align=right → nilai t = rata-rata m pengamatan sebelumnya
    ma_v <- as.numeric(rollmean(dy_train, k=m, align="right", fill=NA))
    # 1-step ahead: fitted[t] = MA_{t-1} → shift kanan 1
    fit_v  <- c(NA, head(ma_v, n_dtr-1))
    res_v  <- dy_train - fit_v
    rmse_v <- sqrt(mean(res_v^2, na.rm=TRUE))
    if (!is.na(rmse_v) && rmse_v < best_ma_d$rmse) {
      best_ma_d <- list(rmse=rmse_v, m=m)
    }
  }

  # Residual dengan m terbaik
  ma_v_best  <- as.numeric(rollmean(dy_train, k=best_ma_d$m, align="right", fill=NA))
  ma_res_d   <- dy_train - c(NA, head(ma_v_best, n_dtr-1))
  ma_diag_d  <- uji_lb(ma_res_d,
                        sprintf("Moving Average diff (m=%d)", best_ma_d$m))
  {
    r <- na.omit(ma_res_d); lag <- min(10, max(1,floor(length(r)/5)))
    bt <- Box.test(r, lag=lag, type="Ljung-Box")
    ma_diag_d <- list(p_value=bt$p.value, lolos=bt$p.value>0.05)
  }

  last_ma_val  <- tail(na.omit(ma_v_best), 1)
  ma_fcast_ms  <- y_last + cumsum(rep(last_ma_val, N_TEST))
  ma_os        <- rolling_os(y_train, y_test, function(y) {
    dy <- diff(y)
    ma_now <- if (length(dy) >= best_ma_d$m) mean(tail(dy, best_ma_d$m)) else mean(dy)
    tail(y,1) + ma_now
  })

  results[[mkt]][["MA_diff"]] <- list(
    lolos    = ma_diag_d$lolos,
    resid    = ma_res_d,
    fcast_ms = ma_fcast_ms,
    fcast_os = ma_os,
    m        = best_ma_d$m
  )

  # ============================================================
  ## SES pada data diferensial – Grid Search alpha
  # ============================================================
  # ============================================================
  # SES pada data diferensial – Grid Search alpha (REVISI)
  # ============================================================
  alpha_grid_s <- seq(0.05, 0.95, by=0.05)
  best_ses_d   <- list(aic=Inf, alpha=NA, fit=NULL)
  
  # Pastikan dy_train valid
  dy_train_clean <- na.omit(dy_train)
  
  if (length(dy_train_clean) > 0) {
    ts_dy <- ts(dy_train_clean, frequency=1)

    for (al in alpha_grid_s) {
      tryCatch({
        sf      <- ses(ts_dy, alpha=al, h=N_TEST, initial="simple")
        aic_val <- sf$model$aic
        if (!is.null(aic_val) && !is.na(aic_val) && aic_val < best_ses_d$aic) {
          best_ses_d <- list(aic=aic_val, alpha=al, fit=sf)
        }
      }, error=function(e) NULL)
    }
  }

  # Fallback jika grid search gagal
  if (is.null(best_ses_d$fit) && length(dy_train_clean) > 0) {
    best_ses_d$fit   <- ses(ts_dy, h=N_TEST)
    best_ses_d$alpha <- best_ses_d$fit$model$par["alpha"]
  }

  # Diagnostik
  if (!is.null(best_ses_d$fit)) {
    uji_lb(residuals(best_ses_d$fit), sprintf("SES diff (α=%.2f)", best_ses_d$alpha))
    
    # Simpan hasil diagnostik
    r_ses <- na.omit(residuals(best_ses_d$fit))
    lag_s <- min(10, max(1, floor(length(r_ses)/5)))
    bt_s  <- Box.test(r_ses, lag=lag_s, type="Ljung-Box")
    ses_diag_d <- list(p_value=bt_s$p.value, lolos=bt_s$p.value > 0.05)

    ses_diff_fcast_ms <- as.numeric(best_ses_d$fit$mean)
    ses_fcast_ms_lv   <- y_last + cumsum(ses_diff_fcast_ms)
    
    # Perbaikan pada rolling_os untuk SES diff
    ses_os <- rolling_os(y_train, y_test, function(y) {
      dy_now <- diff(y)
      if (length(dy_now) < 1) return(tail(y, 1)) # Proteksi jika diff kosong
      ts_dy2 <- ts(dy_now, frequency=1)
      pred   <- tryCatch(
        as.numeric(ses(ts_dy2, alpha=best_ses_d$alpha, h=1, initial="simple")$mean),
        error=function(e) mean(dy_now) # Fallback ke rata-rata jika ses gagal
      )
      tail(y,1) + pred
    })

    results[[mkt]][["SES_diff"]] <- list(
      lolos    = ses_diag_d$lolos,
      resid    = residuals(best_ses_d$fit),
      fcast_ms = ses_fcast_ms_lv,
      fcast_os = ses_os,
      alpha    = best_ses_d$alpha
    )
  }

  cat("```\n")
}

5.3 Pasar: Rotterdam (Diferensial)

5.3.1 Uji Diagnostik

Uji Ljung-Box (p > 0,05 = lolos)

  Simple Average (diff)               | p = 0.7185 | ✅ Tidak ada autokorelasi
  Moving Average diff (m=11)          | p = 0.9697 | ✅ Tidak ada autokorelasi
  SES diff (α=0.00)                  | p = 0.7185 | ✅ Tidak ada autokorelasi

5.4 Pasar: MDEX (Diferensial)

5.4.1 Uji Diagnostik

Uji Ljung-Box (p > 0,05 = lolos)

  Simple Average (diff)               | p = 0.3685 | ✅ Tidak ada autokorelasi
  Moving Average diff (m=3)           | p = 0.1769 | ✅ Tidak ada autokorelasi
  SES diff (α=0.00)                  | p = 0.3686 | ✅ Tidak ada autokorelasi

5.5 Pasar: ICDX (Diferensial)

5.5.1 Uji Diagnostik

Uji Ljung-Box (p > 0,05 = lolos)

  Simple Average (diff)               | p = 0.7123 | ✅ Tidak ada autokorelasi
  Moving Average diff (m=10)          | p = 0.4692 | ✅ Tidak ada autokorelasi
  SES diff (α=0.00)                  | p = 0.7123 | ✅ Tidak ada autokorelasi

6 Ringkasan Uji Diagnostik

model_labels <- c(
  "Naive Simple" = "Naive Simple",
  "Naive Trend"  = "Naive Trend",
  "Naive Drift"  = "Naive Drift",
  "Holt"         = "Holt's Method",
  "DMA"          = "Double MA",
  "SA_diff"      = "Simple Avg (Diff)",
  "MA_diff"      = "Moving Avg (Diff)",
  "SES_diff"     = "SES (Diff)"
)

all_mdl_keys <- names(model_labels)

diag_rows <- lapply(MARKETS, function(mkt) {
  lapply(all_mdl_keys, function(mdl) {
    r <- results[[mkt]][[mdl]]
    if (is.null(r)) return(NULL)
    data.frame(
      Pasar = mkt, Model = model_labels[mdl],
      Lolos = ifelse(r$lolos, "✅ Lolos", "❌ Tidak Lolos"),
      stringsAsFactors = FALSE
    )
  })
})
diag_df <- do.call(rbind, Filter(Negate(is.null), unlist(diag_rows, recursive=FALSE)))

diag_wide <- diag_df %>%
  pivot_wider(names_from=Pasar, values_from=Lolos)

kbl(diag_wide,
    caption="Hasil Uji Diagnostik Ljung-Box (H₀: tidak ada autokorelasi residual)") %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered"), full_width=FALSE) %>%
  column_spec(1, bold=TRUE)

Hasil Uji Diagnostik Ljung-Box (H₀: tidak ada autokorelasi residual)
Model	Rotterdam	MDEX	ICDX
Naive Simple	✅ Lolos	✅ Lolos	✅ Lolos
Naive Trend	❌ Tidak Lolos	❌ Tidak Lolos	❌ Tidak Lolos
Naive Drift	✅ Lolos	✅ Lolos	✅ Lolos
Holt’s Method	✅ Lolos	✅ Lolos	✅ Lolos
Double MA	❌ Tidak Lolos	❌ Tidak Lolos	✅ Lolos
Simple Avg (Diff)	✅ Lolos	✅ Lolos	✅ Lolos
Moving Avg (Diff)	✅ Lolos	✅ Lolos	✅ Lolos
SES (Diff)	✅ Lolos	✅ Lolos	✅ Lolos

7 Evaluasi pada Data Uji

Metrik dihitung untuk semua model (dengan catatan status diagnostik). Model terbaik dipilih berdasarkan RMSE terendah.

eval_list <- list()

for (mkt in MARKETS) {
  y_test <- df_test[[mkt]]
  for (mdl in all_mdl_keys) {
    r <- results[[mkt]][[mdl]]
    if (is.null(r)) next
    if (is.null(r$fcast_ms) || is.null(r$fcast_os)) next
    ms <- calc_metrics(y_test, r$fcast_ms)
    os <- calc_metrics(y_test, r$fcast_os)
    eval_list[[length(eval_list)+1]] <- data.frame(
      Pasar    = mkt,
      Model    = model_labels[mdl],
      Diagnostik = ifelse(r$lolos, "✅","❌"),
      RMSE_ms  = round(ms["RMSE"],4),
      MSE_ms   = round(ms["MSE"],4),
      MAPE_ms  = round(ms["MAPE"],4),
      RMSE_os  = round(os["RMSE"],4),
      MSE_os   = round(os["MSE"],4),
      MAPE_os  = round(os["MAPE"],4),
      row.names = NULL,
      stringsAsFactors = FALSE
    )
  }
}

eval_df <- do.call(rbind, eval_list)

# Tampilkan per pasar
for (mkt in MARKETS) {
  cat("\n**Pasar:", mkt, "**\n\n")
  tbl_mkt <- eval_df %>% filter(Pasar == mkt) %>% select(-Pasar) %>%
    arrange(RMSE_ms)
  
  print(
    kbl(tbl_mkt,
        col.names = c("Model","Diag","RMSE","MSE","MAPE (%)","RMSE","MSE","MAPE (%)"),
        caption = paste("Evaluasi –", mkt)) %>%
      kable_styling(bootstrap_options=c("striped","hover"), full_width=FALSE) %>%
      add_header_above(c(" "=2, "Multi-step Ahead"=3, "One-step Ahead"=3)) %>%
      row_spec(which(tbl_mkt$RMSE_ms == min(tbl_mkt$RMSE_ms)), bold=TRUE,
               background="#d4edda", color="#155724")
  )
}

## 
## **Pasar: Rotterdam **
## 
## <table class="table table-striped table-hover" style="color: black; width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Evaluasi – Rotterdam</caption>
##  <thead>
## <tr>
## <th style="empty-cells: hide;border-bottom:hidden;" colspan="2"></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Multi-step Ahead</div></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">One-step Ahead</div></th>
## </tr>
##   <tr>
##    <th style="text-align:left;"> Model </th>
##    <th style="text-align:left;"> Diag </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> Naive Trend </td>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> ❌ </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 98.4948 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 9701.224 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5.9785 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 77.5572 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 6015.121 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5.4862 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Moving Avg (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 152.8808 </td>
##    <td style="text-align:right;"> 23372.545 </td>
##    <td style="text-align:right;"> 10.2187 </td>
##    <td style="text-align:right;"> 73.7927 </td>
##    <td style="text-align:right;"> 5445.360 </td>
##    <td style="text-align:right;"> 4.8189 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Drift </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 158.3522 </td>
##    <td style="text-align:right;"> 25075.418 </td>
##    <td style="text-align:right;"> 10.5997 </td>
##    <td style="text-align:right;"> 72.7591 </td>
##    <td style="text-align:right;"> 5293.880 </td>
##    <td style="text-align:right;"> 4.7708 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Simple Avg (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 158.3522 </td>
##    <td style="text-align:right;"> 25075.418 </td>
##    <td style="text-align:right;"> 10.5997 </td>
##    <td style="text-align:right;"> 72.7591 </td>
##    <td style="text-align:right;"> 5293.880 </td>
##    <td style="text-align:right;"> 4.7708 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> SES (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 158.3953 </td>
##    <td style="text-align:right;"> 25089.072 </td>
##    <td style="text-align:right;"> 10.6027 </td>
##    <td style="text-align:right;"> 84.2519 </td>
##    <td style="text-align:right;"> 7098.389 </td>
##    <td style="text-align:right;"> 5.4478 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Holt's Method </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 167.8618 </td>
##    <td style="text-align:right;"> 28177.579 </td>
##    <td style="text-align:right;"> 11.2783 </td>
##    <td style="text-align:right;"> 94.4996 </td>
##    <td style="text-align:right;"> 8930.171 </td>
##    <td style="text-align:right;"> 6.2909 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Simple </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 227.1504 </td>
##    <td style="text-align:right;"> 51597.308 </td>
##    <td style="text-align:right;"> 15.8750 </td>
##    <td style="text-align:right;"> 74.6941 </td>
##    <td style="text-align:right;"> 5579.211 </td>
##    <td style="text-align:right;"> 4.8662 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Double MA </td>
##    <td style="text-align:left;"> ❌ </td>
##    <td style="text-align:right;"> 371.8764 </td>
##    <td style="text-align:right;"> 138292.072 </td>
##    <td style="text-align:right;"> 26.8102 </td>
##    <td style="text-align:right;"> 110.3024 </td>
##    <td style="text-align:right;"> 12166.628 </td>
##    <td style="text-align:right;"> 7.6963 </td>
##   </tr>
## </tbody>
## </table>
## **Pasar: MDEX **
## 
## <table class="table table-striped table-hover" style="color: black; width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Evaluasi – MDEX</caption>
##  <thead>
## <tr>
## <th style="empty-cells: hide;border-bottom:hidden;" colspan="2"></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Multi-step Ahead</div></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">One-step Ahead</div></th>
## </tr>
##   <tr>
##    <th style="text-align:left;"> Model </th>
##    <th style="text-align:left;"> Diag </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> Moving Avg (Diff) </td>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> ✅ </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 95.2610 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 9074.657 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 7.0157 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 62.0851 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 3854.559 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5.1596 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Double MA </td>
##    <td style="text-align:left;"> ❌ </td>
##    <td style="text-align:right;"> 104.8396 </td>
##    <td style="text-align:right;"> 10991.351 </td>
##    <td style="text-align:right;"> 8.9151 </td>
##    <td style="text-align:right;"> 69.1219 </td>
##    <td style="text-align:right;"> 4777.839 </td>
##    <td style="text-align:right;"> 5.5949 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Holt's Method </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 106.1584 </td>
##    <td style="text-align:right;"> 11269.611 </td>
##    <td style="text-align:right;"> 7.7831 </td>
##    <td style="text-align:right;"> 64.1131 </td>
##    <td style="text-align:right;"> 4110.485 </td>
##    <td style="text-align:right;"> 4.8112 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Trend </td>
##    <td style="text-align:left;"> ❌ </td>
##    <td style="text-align:right;"> 117.6731 </td>
##    <td style="text-align:right;"> 13846.951 </td>
##    <td style="text-align:right;"> 10.3364 </td>
##    <td style="text-align:right;"> 53.2096 </td>
##    <td style="text-align:right;"> 2831.257 </td>
##    <td style="text-align:right;"> 4.3404 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> SES (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 123.6882 </td>
##    <td style="text-align:right;"> 15298.760 </td>
##    <td style="text-align:right;"> 9.7340 </td>
##    <td style="text-align:right;"> 63.6458 </td>
##    <td style="text-align:right;"> 4050.787 </td>
##    <td style="text-align:right;"> 4.7708 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Drift </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 123.6896 </td>
##    <td style="text-align:right;"> 15299.108 </td>
##    <td style="text-align:right;"> 9.7342 </td>
##    <td style="text-align:right;"> 53.4082 </td>
##    <td style="text-align:right;"> 2852.431 </td>
##    <td style="text-align:right;"> 4.2889 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Simple Avg (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 123.6896 </td>
##    <td style="text-align:right;"> 15299.108 </td>
##    <td style="text-align:right;"> 9.7342 </td>
##    <td style="text-align:right;"> 53.4082 </td>
##    <td style="text-align:right;"> 2852.431 </td>
##    <td style="text-align:right;"> 4.2889 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Simple </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 165.4448 </td>
##    <td style="text-align:right;"> 27371.987 </td>
##    <td style="text-align:right;"> 14.5378 </td>
##    <td style="text-align:right;"> 54.1408 </td>
##    <td style="text-align:right;"> 2931.230 </td>
##    <td style="text-align:right;"> 4.3666 </td>
##   </tr>
## </tbody>
## </table>
## **Pasar: ICDX **
## 
## <table class="table table-striped table-hover" style="color: black; width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Evaluasi – ICDX</caption>
##  <thead>
## <tr>
## <th style="empty-cells: hide;border-bottom:hidden;" colspan="2"></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">Multi-step Ahead</div></th>
## <th style="border-bottom:hidden;padding-bottom:0; padding-left:3px;padding-right:3px;text-align: center; " colspan="3"><div style="border-bottom: 1px solid #ddd; padding-bottom: 5px; ">One-step Ahead</div></th>
## </tr>
##   <tr>
##    <th style="text-align:left;"> Model </th>
##    <th style="text-align:left;"> Diag </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##    <th style="text-align:right;"> RMSE </th>
##    <th style="text-align:right;"> MSE </th>
##    <th style="text-align:right;"> MAPE (%) </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> Holt's Method </td>
##    <td style="text-align:left;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> ✅ </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 73.1041 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5344.208 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5.3428 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 62.7612 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 3938.971 </td>
##    <td style="text-align:right;font-weight: bold;color: rgba(21, 87, 36, 255) !important;background-color: rgba(212, 237, 218, 255) !important;"> 5.4686 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Drift </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 73.9778 </td>
##    <td style="text-align:right;"> 5472.718 </td>
##    <td style="text-align:right;"> 5.3995 </td>
##    <td style="text-align:right;"> 52.4639 </td>
##    <td style="text-align:right;"> 2752.456 </td>
##    <td style="text-align:right;"> 4.5278 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Simple Avg (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 73.9778 </td>
##    <td style="text-align:right;"> 5472.718 </td>
##    <td style="text-align:right;"> 5.3995 </td>
##    <td style="text-align:right;"> 52.4639 </td>
##    <td style="text-align:right;"> 2752.456 </td>
##    <td style="text-align:right;"> 4.5278 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> SES (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 73.9882 </td>
##    <td style="text-align:right;"> 5474.259 </td>
##    <td style="text-align:right;"> 5.4001 </td>
##    <td style="text-align:right;"> 62.7605 </td>
##    <td style="text-align:right;"> 3938.885 </td>
##    <td style="text-align:right;"> 5.4686 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Moving Avg (Diff) </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 82.6896 </td>
##    <td style="text-align:right;"> 6837.566 </td>
##    <td style="text-align:right;"> 7.8208 </td>
##    <td style="text-align:right;"> 54.1314 </td>
##    <td style="text-align:right;"> 2930.209 </td>
##    <td style="text-align:right;"> 4.6591 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Simple </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 96.3729 </td>
##    <td style="text-align:right;"> 9287.731 </td>
##    <td style="text-align:right;"> 8.4865 </td>
##    <td style="text-align:right;"> 52.4980 </td>
##    <td style="text-align:right;"> 2756.039 </td>
##    <td style="text-align:right;"> 4.4899 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Naive Trend </td>
##    <td style="text-align:left;"> ❌ </td>
##    <td style="text-align:right;"> 181.4598 </td>
##    <td style="text-align:right;"> 32927.666 </td>
##    <td style="text-align:right;"> 18.5431 </td>
##    <td style="text-align:right;"> 64.3121 </td>
##    <td style="text-align:right;"> 4136.043 </td>
##    <td style="text-align:right;"> 5.0660 </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> Double MA </td>
##    <td style="text-align:left;"> ✅ </td>
##    <td style="text-align:right;"> 196.5525 </td>
##    <td style="text-align:right;"> 38632.868 </td>
##    <td style="text-align:right;"> 19.9913 </td>
##    <td style="text-align:right;"> 68.1661 </td>
##    <td style="text-align:right;"> 4646.613 </td>
##    <td style="text-align:right;"> 5.6483 </td>
##   </tr>
## </tbody>
## </table>

7.1 Model Terbaik per Pasar

# Pilih berdasarkan RMSE terkecil (multi-step)
best_ms_df <- eval_df %>%
  group_by(Pasar) %>%
  slice_min(RMSE_ms, n=1) %>%
  select(Pasar, Model, Diagnostik, RMSE_ms, MSE_ms, MAPE_ms) %>%
  ungroup()

best_os_df <- eval_df %>%
  group_by(Pasar) %>%
  slice_min(RMSE_os, n=1) %>%
  select(Pasar, Model, Diagnostik, RMSE_os, MSE_os, MAPE_os) %>%
  ungroup()

kbl(best_ms_df,
    caption="Model Terbaik – Multi-step Ahead") %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered"), full_width=FALSE) %>%
  column_spec(1:2, bold=TRUE)

Model Terbaik – Multi-step Ahead
Pasar	Model	Diagnostik	RMSE_ms	MSE_ms	MAPE_ms
ICDX	Holt’s Method	✅	73.1041	5344.208	5.3428
MDEX	Moving Avg (Diff)	✅	95.2610	9074.657	7.0157
Rotterdam	Naive Trend	❌	98.4948	9701.224	5.9785

kbl(best_os_df,
    caption="Model Terbaik – One-step Ahead") %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered"), full_width=FALSE) %>%
  column_spec(1:2, bold=TRUE)

Model Terbaik – One-step Ahead
Pasar	Model	Diagnostik	RMSE_os	MSE_os	MAPE_os
ICDX	Naive Drift	✅	52.4639	2752.456	4.5278
ICDX	Simple Avg (Diff)	✅	52.4639	2752.456	4.5278
MDEX	Naive Trend	❌	53.2096	2831.257	4.3404
Rotterdam	Naive Drift	✅	72.7591	5293.880	4.7708
Rotterdam	Simple Avg (Diff)	✅	72.7591	5293.880	4.7708

7.2 Plot Perbandingan Forecast pada Data Uji

test_dates <- df_test$Tanggal

plot_list <- lapply(MARKETS, function(mkt) {
  y_test    <- df_test[[mkt]]
  best_mdl  <- best_ms_df$Model[best_ms_df$Pasar == mkt]
  # Temukan key dari label
  mdl_key   <- names(model_labels)[model_labels == best_mdl]
  fcast_ms  <- results[[mkt]][[mdl_key]]$fcast_ms
  fcast_os  <- results[[mkt]][[mdl_key]]$fcast_os

  df_p <- data.frame(
    Tanggal = rep(test_dates, 3),
    Nilai   = c(y_test, fcast_ms, fcast_os),
    Seri    = rep(c("Aktual","Multi-step","One-step"), each=N_TEST)
  )

  ggplot(df_p, aes(x=Tanggal, y=Nilai, color=Seri, linetype=Seri)) +
    geom_line(linewidth=0.9) +
    geom_point(size=1.5) +
    scale_color_manual(values=c(Aktual="black","Multi-step"="red","One-step"="blue")) +
    scale_linetype_manual(values=c(Aktual="solid","Multi-step"="dashed","One-step"="dotted")) +
    labs(title=paste(mkt, "–", best_mdl),
         subtitle=sprintf("RMSE multi-step=%.2f | RMSE one-step=%.2f",
                          results[[mkt]][[mdl_key]]$fcast_ms |> (\(f) calc_metrics(y_test,f)["RMSE"])(),
                          results[[mkt]][[mdl_key]]$fcast_os |> (\(f) calc_metrics(y_test,f)["RMSE"])()),
         x=NULL, y="USD/ton", color=NULL, linetype=NULL) +
    theme(legend.position="top")
})

(plot_list[[1]] / plot_list[[2]] / plot_list[[3]]) +
  plot_annotation(title="Perbandingan Ramalan pada Data Uji")

8 Peramalan hingga Januari 2028

Model terbaik (multi-step) dilatih ulang pada seluruh data (79 bulan), kemudian digunakan untuk meramal Juli 2026 – Januari 2028 (19 bulan ke depan).

last_date      <- max(df$Tanggal)
n_horizon      <- 19   # Jul 2026 – Jan 2028

# Buat tanggal forecast
fc_dates <- Reduce(function(d, i) c(d, add_month(tail(d,1))),
                    seq_len(n_horizon - 1),
                    init = add_month(last_date))

final_fcast <- list()

for (mkt in MARKETS) {
  y_all    <- df[[mkt]]
  n_all    <- length(y_all)
  y_last_all <- tail(y_all, 1)

  best_mdl  <- best_ms_df$Model[best_ms_df$Pasar == mkt]
  mdl_key   <- names(model_labels)[model_labels == best_mdl]

  cat(sprintf("\n%-12s → Model terbaik: %s\n", mkt, best_mdl))

  ts_all <- make_ts(y_all, df$Tanggal[1])

  fcast_vals <- tryCatch({
    switch(mdl_key,
      "Naive Simple" = rep(y_last_all, n_horizon),
      "Naive Trend"  = y_last_all + (y_all[n_all]-y_all[n_all-1]) * seq_len(n_horizon),
      "Naive Drift"  = {
        drift <- (y_last_all - y_all[1])/(n_all-1)
        y_last_all + drift * seq_len(n_horizon)
      },
      "Holt" = {
        al <- results[[mkt]][["Holt"]]$alpha
        bt <- results[[mkt]][["Holt"]]$beta
        as.numeric(holt(ts_all, alpha=al, beta=bt, h=n_horizon, initial="simple")$mean)
      },
      "DMA" = {
        m <- results[[mkt]][["DMA"]]$m
        dma_forecast_h(dma_compute(y_all, m), n_horizon)
      },
      "SA_diff"  = {
        mean_dy <- mean(diff(y_all))
        y_last_all + cumsum(rep(mean_dy, n_horizon))
      },
      "MA_diff"  = {
        m <- results[[mkt]][["MA_diff"]]$m
        last_ma <- mean(tail(diff(y_all), m))
        y_last_all + cumsum(rep(last_ma, n_horizon))
      },
      "SES_diff" = {
        al <- results[[mkt]][["SES_diff"]]$alpha
        ts_dy <- ts(diff(y_all), frequency=1)
        pd <- as.numeric(ses(ts_dy, alpha=al, h=n_horizon, initial="simple")$mean)
        y_last_all + cumsum(pd)
      },
      rep(y_last_all, n_horizon)  # fallback
    )
  }, error = function(e) {
    cat("  ⚠️ Error:", conditionMessage(e), "→ fallback Naive Simple\n")
    rep(y_last_all, n_horizon)
  })

  final_fcast[[mkt]] <- fcast_vals
}

## 
## Rotterdam    → Model terbaik: Naive Trend
## 
## MDEX         → Model terbaik: Moving Avg (Diff)
## 
## ICDX         → Model terbaik: Holt's Method

# ----------------------------------------------------------------
# Plot gabungan: historis + forecast
# ----------------------------------------------------------------
plot_fcast_list <- lapply(MARKETS, function(mkt) {
  fv <- final_fcast[[mkt]]

  df_h <- data.frame(Tanggal=df$Tanggal, Nilai=df[[mkt]], Tipe="Historis")
  df_f <- data.frame(Tanggal=fc_dates,   Nilai=fv,       Tipe="Ramalan")
  df_c <- rbind(df_h, df_f)

  best_mdl <- best_ms_df$Model[best_ms_df$Pasar == mkt]

  ggplot(df_c, aes(x=Tanggal, y=Nilai, color=Tipe)) +
    geom_line(linewidth=0.9) +
    geom_point(data=df_f, size=2.5, shape=16) +
    geom_vline(xintercept=as.numeric(last_date), linetype="dotted",
               color="grey40", linewidth=0.7) +
    scale_color_manual(values=c(Historis="steelblue", Ramalan="firebrick")) +
    labs(title=paste(mkt, "–", best_mdl),
         x=NULL, y="USD/ton", color=NULL) +
    theme(legend.position="top")
})

(plot_fcast_list[[1]] / plot_fcast_list[[2]] / plot_fcast_list[[3]]) +
  plot_annotation(
    title   = "Ramalan Harga CPO – Juli 2026 hingga Januari 2028",
    caption = "Garis putus-putus vertikal: batas data historis"
  )

8.1 Tabel Ramalan Akhir

tbl_final <- data.frame(
  Bulan     = format(fc_dates, "%B %Y"),
  Rotterdam = round(final_fcast[["Rotterdam"]], 2),
  MDEX      = round(final_fcast[["MDEX"]], 2),
  ICDX      = round(final_fcast[["ICDX"]], 2)
)

kbl(tbl_final,
    col.names = c("Bulan","Rotterdam (USD/ton)","MDEX (USD/ton)","ICDX (USD/ton)"),
    caption   = "Ramalan Harga Rata-rata Bulanan CPO (Juli 2026 – Januari 2028)") %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered"), full_width=FALSE) %>%
  column_spec(1, bold=TRUE) %>%
  row_spec(nrow(tbl_final), bold=TRUE, background="#fff3cd")

Ramalan Harga Rata-rata Bulanan CPO (Juli 2026 – Januari 2028)
Bulan	Rotterdam (USD/ton)	MDEX (USD/ton)	ICDX (USD/ton)
June 2026	1689.43	1190.16	999.04
July 2026	1842.75	1243.54	1042.02
August 2026	1996.07	1296.93	1085.00
September 2026	2149.39	1350.32	1127.98
October 2026	2302.71	1403.70	1170.95
November 2026	2456.03	1457.09	1213.93
December 2026	2609.35	1510.48	1256.91
January 2027	2762.67	1563.86	1299.89
February 2027	2915.99	1617.25	1342.87
March 2027	3069.31	1670.64	1385.85
April 2027	3222.63	1724.02	1428.83
May 2027	3375.95	1777.41	1471.81
June 2027	3529.27	1830.80	1514.79
July 2027	3682.59	1884.18	1557.77
August 2027	3835.91	1937.57	1600.75
September 2027	3989.23	1990.96	1643.73
October 2027	4142.55	2044.34	1686.71
November 2027	4295.87	2097.73	1729.69
December 2027	4449.19	2151.12	1772.67

9 Informasi Parameter Model Terbaik

param_list <- lapply(MARKETS, function(mkt) {
  mdl_key  <- names(model_labels)[model_labels == best_ms_df$Model[best_ms_df$Pasar==mkt]]
  r        <- results[[mkt]][[mdl_key]]
  params   <- switch(mdl_key,
    "Holt"     = sprintf("α=%.2f, β=%.2f", r$alpha, r$beta),
    "DMA"      = sprintf("m=%d", r$m),
    "MA_diff"  = sprintf("m=%d", r$m),
    "SES_diff" = sprintf("α=%.2f", r$alpha),
    "–"
  )
  data.frame(Pasar=mkt, `Model (Multi-step)`=best_ms_df$Model[best_ms_df$Pasar==mkt],
             Parameter=params,
             `RMSE Test`=round(best_ms_df$RMSE_ms[best_ms_df$Pasar==mkt],4),
             check.names=FALSE, stringsAsFactors=FALSE)
})

kbl(do.call(rbind, param_list),
    caption="Ringkasan Parameter Model Terbaik per Pasar") %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered"), full_width=FALSE) %>%
  column_spec(1:2, bold=TRUE)

Ringkasan Parameter Model Terbaik per Pasar
Pasar	Model (Multi-step)	Parameter	RMSE Test
Rotterdam	Naive Trend	–	98.4948
MDEX	Moving Avg (Diff)	m=3	95.2610
ICDX	Holt’s Method	α=1.00, β=0.00	73.1041

10 Informasi Sesi R

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8  LC_CTYPE=English_Indonesia.utf8   
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=English_Indonesia.utf8    
## 
## time zone: Asia/Jakarta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.3.2  scales_1.4.0     kableExtra_1.4.0 knitr_1.50      
##  [5] zoo_1.8-14       tidyr_1.3.2      dplyr_1.2.1      ggplot2_4.0.3   
##  [9] forecast_8.24.0  readxl_1.4.5    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     xml2_1.3.8         stringi_1.8.7     
##  [5] lattice_0.22-6     digest_0.6.37      magrittr_2.0.3     evaluate_1.0.4    
##  [9] grid_4.4.1         RColorBrewer_1.1-3 fastmap_1.2.0      cellranger_1.1.0  
## [13] jsonlite_2.0.0     nnet_7.3-19        purrr_1.0.4        viridisLite_0.4.3 
## [17] textshaping_1.0.1  jquerylib_0.1.4    cli_3.6.3          rlang_1.2.0       
## [21] withr_3.0.2        cachem_1.1.0       yaml_2.3.10        tools_4.4.1       
## [25] parallel_4.4.1     colorspace_2.1-1   curl_6.4.0         vctrs_0.7.1       
## [29] R6_2.6.1           lifecycle_1.0.5    stringr_1.6.0      tseries_0.10-61   
## [33] pkgconfig_2.0.3    urca_1.3-4         pillar_1.11.0      bslib_0.9.0       
## [37] gtable_0.3.6       glue_1.8.0         quantmod_0.4.28    Rcpp_1.0.13       
## [41] systemfonts_1.3.2  xfun_0.56          tibble_3.2.1       lmtest_0.9-40     
## [45] tidyselect_1.2.1   rstudioapi_0.17.1  farver_2.1.2       nlme_3.1-164      
## [49] htmltools_0.5.8.1  labeling_0.4.3     svglite_2.2.2      rmarkdown_2.29    
## [53] xts_0.14.1         timeDate_4041.110  fracdiff_1.5-3     compiler_4.4.1    
## [57] S7_0.2.1           quadprog_1.5-8     TTR_0.24.4

Peramalan Harga CPO: Studi Komparatif Metode Deret Waktu

Tiga Pasar: Rotterdam (EU), MDEX Malaysia, ICDX Indonesia

04 May 2026

1 Persiapan Data

2 Eksplorasi Data

2.1 Plot Deret Waktu

2.2 Plot ACF – Data Asli

3 Fungsi Bantu

4 Pemodelan pada Data Asli

4.1 Pasar: Rotterdam

4.1.1 Uji Diagnostik

4.2 Pasar: MDEX

4.2.1 Uji Diagnostik

4.3 Pasar: ICDX

4.3.1 Uji Diagnostik

5 Diferensiasi dan Pemodelan Data Diferensial

5.1 Plot Diferensial dan ACF

5.2 Pemodelan Data Diferensial

5.3 Pasar: Rotterdam (Diferensial)

5.3.1 Uji Diagnostik

5.4 Pasar: MDEX (Diferensial)

5.4.1 Uji Diagnostik

5.5 Pasar: ICDX (Diferensial)

5.5.1 Uji Diagnostik

6 Ringkasan Uji Diagnostik

7 Evaluasi pada Data Uji

7.1 Model Terbaik per Pasar

7.2 Plot Perbandingan Forecast pada Data Uji

8 Peramalan hingga Januari 2028

8.1 Tabel Ramalan Akhir

9 Informasi Parameter Model Terbaik

10 Informasi Sesi R