Peramalan jumlah sesi akses LMS per bulan dengan exponential smoothing

# Contoh kode R: Exponential Smoothing (ETS)
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(forecast)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis: deret bulanan 8 tahun dengan tren dan musiman ----
bulan <- seq(as.Date("2017-01-01"), by = "month", length.out = 96)
tt <- seq_along(bulan)
tumbuh <- runif(1, 0.15, 0.5)                                # pertumbuhan total selama 8 tahun
level <- 1 + tumbuh * tt / 96 + cumsum(rnorm(96, 0, 0.004))   # level dengan sedikit pergeseran acak
pola <- c(-1, -0.8, -0.2, 0.1, 0.3, 0.6, 0.9, 0.7, 0.2, -0.1, 0.3, 1.2)   # pola musiman dasar
geser <- sample(0:11, 1)
musim <- 1 + runif(1, 0.06, 0.14) * pola[((tt - 1 + geser) %% 12) + 1] / sd(pola)
galat <- exp(rnorm(96, 0, runif(1, 0.025, 0.045)))
d <- data.frame(bulan = bulan, monthly_lms_sessions = round(max(100, 1) * level * musim * galat, 3))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: ETS otomatis (pemilihan error-tren-musiman berbasis AICc) ----
y <- ts(d$monthly_lms_sessions, start = c(2017, 1), frequency = 12)
latih <- window(y, end = c(2023, 12))
uji <- window(y, start = c(2024, 1))
fit <- ets(latih)                                  # model ETS(E,T,S) dipilih otomatis
print(summary(fit))
fc <- forecast(fit, h = length(uji), level = c(80, 95))
ak <- accuracy(fc, uji)
print(ak)
# pembanding: naive musiman (nilai bulan yang sama tahun lalu)
sn <- snaive(latih, h = length(uji))
ak_sn <- accuracy(sn, uji)
print(ak_sn)
n_par <- sum(names(fit$par) %in% c("alpha", "beta", "gamma", "phi"))
lb <- Box.test(residuals(fit), lag = 24, type = "Ljung-Box", fitdf = n_par)
print(lb)
fit_all <- ets(y, model = fit, use.initial.values = FALSE)    # struktur model yang sama, data lengkap
fc12 <- forecast(fit_all, h = 12)
alpha_hat <- unname(fit$par["alpha"])


# ---- Grafik ----
# 1) ramalan ETS dengan interval 80% dan 95% vs data uji
waktu <- function(x) as.numeric(time(x))
dh <- data.frame(t = waktu(latih), y = as.numeric(latih))
dh <- dh[dh$t >= 2021, ]                           # tampilkan 3 tahun terakhir data latih
du <- data.frame(t = waktu(uji), y = as.numeric(uji))
df <- data.frame(t = waktu(fc$mean), m = as.numeric(fc$mean),
                 lo80 = as.numeric(fc$lower[, 1]), hi80 = as.numeric(fc$upper[, 1]),
                 lo95 = as.numeric(fc$lower[, 2]), hi95 = as.numeric(fc$upper[, 2]))
p1 <- ggplot() +
  geom_ribbon(data = df, aes(t, ymin = lo95, ymax = hi95), fill = pal[2], alpha = 0.15) +
  geom_ribbon(data = df, aes(t, ymin = lo80, ymax = hi80), fill = pal[2], alpha = 0.3) +
  geom_line(data = dh, aes(t, y), colour = pal[1], linewidth = 0.7) +
  geom_line(data = df, aes(t, m), colour = pal[2], linewidth = 1) +
  geom_point(data = du, aes(t, y), colour = warna$teks, size = 1.8) +
  scale_x_continuous(breaks = 2021:2025) +
  labs(title = paste(fit$method, "forecast, 80/95% PI (dots = test)"), x = "Year", y = "monthly_lms_sessions")
# 2) komponen state ETS: level, slope, musiman
st <- fit$states[-1, , drop = FALSE]
kom <- data.frame(t = waktu(latih), Level = st[, "l"])
if ("b" %in% colnames(st)) kom$Slope <- st[, "b"]
if ("s1" %in% colnames(st)) kom$Season <- st[, "s1"]
kom$Error <- as.numeric(residuals(fit))
dk <- do.call(rbind, lapply(names(kom)[-1], function(k) data.frame(t = kom$t, komponen = k, nilai = kom[[k]])))
dk$komponen <- factor(dk$komponen, levels = names(kom)[-1])
p2 <- ggplot(dk, aes(t, nilai)) +
  geom_line(colour = pal[3], linewidth = 0.7) +
  facet_wrap(~komponen, ncol = 1, scales = "free_y", strip.position = "right") +
  theme(axis.text.y = element_text(size = rel(0.65))) +
  labs(title = "ETS components", x = "Year", y = NULL)
# 3) galat absolut per bulan uji: ETS vs naive musiman
dg <- rbind(data.frame(t = format(time(uji)), bln = month.abb[cycle(uji)], g = abs(as.numeric(uji - fc$mean)), model = "ETS"),
            data.frame(t = format(time(uji)), bln = month.abb[cycle(uji)], g = abs(as.numeric(uji - sn$mean)), model = "Seasonal naive"))
dg$bln <- factor(dg$bln, levels = month.abb[cycle(uji)])
p3 <- ggplot(dg, aes(bln, g, fill = model)) +
  geom_col(position = position_dodge(0.8), width = 0.75) +
  scale_fill_manual(values = c(ETS = pal[2], "Seasonal naive" = pal[4]), name = NULL) +
  theme(axis.text.x = element_text(size = rel(0.7))) +
  labs(title = "Absolute error by test month", x = "2024", y = "|Error|")
p_hasil <- if (tata == 1) ((p1 / p3) | p2) + plot_layout(widths = c(1.5, 1)) else (p2 | (p1 / p3)) + plot_layout(widths = c(1, 1.5))

print(p_hasil)