Peramalan kunjungan portal pembelajaran daring universitas dengan ARIMA

# Contoh kode R: ARIMA (Box-Jenkins)
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 ----
bulan <- seq(as.Date("2017-01-01"), by = "month", length.out = 96)
tt <- seq_along(bulan)
musim <- sin(2 * pi * tt / 12 + runif(1, 0, 2 * pi)) * runif(1, 0.3, 0.8)
galat <- as.numeric(arima.sim(list(ar = runif(1, 0.3, 0.7)), n = 96, sd = 0.35))
z <- runif(1, 0.008, 0.02) * tt + musim + galat
d <- data.frame(bulan = bulan, monthly_portal_visits = round(100 + 15 * 0.6 * (z - mean(z)), 3))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: ARIMA Box-Jenkins (identifikasi otomatis + evaluasi out-of-sample) ----
y <- ts(d$monthly_portal_visits, start = c(2017, 1), frequency = 12)
latih <- window(y, end = c(2023, 12))
uji <- window(y, start = c(2024, 1))
fit <- auto.arima(latih, stepwise = TRUE, approximation = FALSE)
print(summary(fit))
fc <- forecast(fit, h = length(uji))
ak <- accuracy(fc, uji)
print(ak)
lb <- Box.test(residuals(fit), lag = 12, type = "Ljung-Box", fitdf = length(coef(fit)))
print(lb)
ord <- arimaorder(fit)
label_model <- sprintf("ARIMA(%d,%d,%d)%s", ord[1], ord[2], ord[3],
                       if (length(ord) > 3) sprintf("(%d,%d,%d)[12]", ord[4], ord[5], ord[6]) else "")
fit_all <- Arima(y, model = fit)          # model yang sama untuk peramalan 12 bulan ke depan
fc12 <- forecast(fit_all, h = 12)


# ---- Grafik ----
dd <- data.frame(t = as.numeric(time(y)), y = as.numeric(y), bagian = ifelse(time(y) >= 2024, "Test", "Training"))
df <- data.frame(t = as.numeric(time(fc$mean)), m = as.numeric(fc$mean),
                 lo = as.numeric(fc$lower[, 2]), hi = as.numeric(fc$upper[, 2]))
p1 <- ggplot() +
  geom_ribbon(data = df, aes(t, ymin = lo, ymax = hi), fill = pal[2], alpha = 0.2) +
  geom_line(data = dd, aes(t, y, colour = bagian), linewidth = 0.8) +
  geom_line(data = df, aes(t, m, colour = "Forecast"), linewidth = 1, linetype = "dashed") +
  scale_colour_manual(values = c(Training = pal[1], Test = warna$teks, Forecast = pal[2]), name = NULL) +
  labs(title = paste(label_model, "forecast vs test data (95% interval)"), x = "Year", y = "monthly_portal_visits")
ac <- acf(residuals(fit), lag.max = 24, plot = FALSE)
da <- data.frame(lag = as.numeric(ac$lag[-1]) * 12, acf = as.numeric(ac$acf[-1]))
batas <- 1.96 / sqrt(length(residuals(fit)))
p2 <- ggplot(da, aes(lag, acf)) +
  geom_hline(yintercept = c(-batas, batas), linetype = "dashed", colour = pal[2]) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_segment(aes(xend = lag, yend = 0), colour = pal[1], linewidth = 1) +
  labs(title = "Residual ACF", x = "Lag (months)", y = "ACF")
p3 <- ggplot(data.frame(r = as.numeric(residuals(fit))), aes(r)) +
  geom_histogram(aes(y = after_stat(density)), bins = 18, fill = pal[3], colour = warna$panel, alpha = 0.85) +
  geom_density(colour = pal[2], linewidth = 0.9) +
  labs(title = "Residual distribution", x = "Residual", y = "Density")
p_hasil <- if (tata == 1) p1 / (p2 | p3) + plot_layout(heights = c(1.3, 1)) else p1 / (p2 | p3)

print(p_hasil)