Volatilitas perubahan harian harga minyak sawit mentah dan kaitannya dengan harga minyak mentah

# Contoh kode R: GARCH Volatility Modelling
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(rugarch)
})
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: imbal hasil harian (%) dengan volatilitas berkelompok ----
n <- 1000
hari <- seq(as.Date("2021-01-04"), by = "day", length.out = 1500)
hari <- hari[!format(hari, "%u") %in% c("6", "7")][1:n]      # hari kerja saja
simulasi_garch <- function(n, omega, a, b, nu, burn = 300) {
  z <- rt(n + burn, df = nu) / sqrt(nu / (nu - 2))             # galat t terstandar (varians 1)
  h <- numeric(n + burn); e <- numeric(n + burn)
  h[1] <- omega / (1 - a - b)
  for (i in 2:(n + burn)) {
    h[i] <- omega + a * e[i - 1]^2 + b * h[i - 1]
    e[i] <- sqrt(h[i]) * z[i]
  }
  e[-seq_len(burn)]
}
a1 <- runif(1, 0.07, 0.14)
b1 <- runif(1, 0.8, 0.975 - a1)                              # persistensi < 0.975
nu <- runif(1, 5, 9)
s_y <- runif(1, 0.9, 2)                                       # simpangan baku tak bersyarat (%)
s_x <- runif(1, 0.6, 1.1)                                     # SD deret terkait (mis. indeks pasar)
x1 <- simulasi_garch(n, s_x^2 * (1 - 0.07 - 0.88), 0.07, 0.88, 8)
d <- data.frame(tanggal = hari,
                crude_oil_price_change = round(x1, 4),
                exchange_rate_change = round(rnorm(n, 0, runif(1, 0.3, 0.6)), 4))
gam <- runif(1, 0.3, 0.8)                                     # kepekaan Y terhadap X1
eps <- simulasi_garch(n, s_y^2 * (1 - a1 - b1), a1, b1, nu)
d$cpo_price_change <- round(runif(1, 0.01, 0.06) + gam * d$crude_oil_price_change + 0.1 * d$exchange_rate_change + eps, 4)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Uji efek ARCH (LM Engle) pada residual model rata-rata ----
r0 <- residuals(lm(cpo_price_change ~ crude_oil_price_change, data = d))
e2 <- r0^2
L <- embed(e2, 6)                                             # e2_t dan 5 lag
arch_lm <- lm(L[, 1] ~ L[, -1])
stat_lm <- nrow(L) * summary(arch_lm)$r.squared
p_arch <- pchisq(stat_lm, df = 5, lower.tail = FALSE)
cat(sprintf("ARCH-LM(5): statistik = %.2f, p = %.4g\n", stat_lm, p_arch))

# ---- Analisis: sGARCH(1,1) dengan galat Student-t, X1 di persamaan rata-rata ----
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                   mean.model = list(armaOrder = c(0, 0), include.mean = TRUE,
                                     external.regressors = matrix(d$crude_oil_price_change)),
                   distribution.model = "std")
fit <- tryCatch(ugarchfit(spec, data = d$cpo_price_change, solver = "solnp"), error = function(e) NULL)
if (is.null(fit) || convergence(fit) != 0) fit <- ugarchfit(spec, data = d$cpo_price_change, solver = "hybrid")   # cadangan
print(fit)
cf <- coef(fit)
pers <- unname(cf["alpha1"] + cf["beta1"])                     # persistensi volatilitas
paruh <- log(0.5) / log(pers)                                 # waktu paruh guncangan (hari)
sig <- as.numeric(sigma(fit))
zres <- as.numeric(residuals(fit, standardize = TRUE))
lb2 <- Box.test(zres^2, lag = 10, type = "Ljung-Box")         # sisa efek ARCH setelah GARCH
print(lb2)
# VaR 95% satu hari ke depan (X1 diasumsikan pada rata-ratanya)
fc <- ugarchforecast(fit, n.ahead = 1, external.forecast = list(mregfor = matrix(mean(d$crude_oil_price_change))))
var95 <- -(as.numeric(fitted(fc)) + as.numeric(sigma(fc)) * qdist("std", 0.05, mu = 0, sigma = 1, shape = cf["shape"]))
cat(sprintf("Persistensi = %.3f, waktu paruh = %.1f hari, VaR 95%% 1 hari = %.3f%%\n", pers, paruh, var95))


# ---- Grafik ----
mu_t <- as.numeric(fitted(fit))
dg <- data.frame(tgl = d$tanggal, r = d$cpo_price_change, lo = mu_t - 2 * sig, hi = mu_t + 2 * sig, sig = sig)
# 1) imbal hasil harian dengan pita +-2 sigma bersyarat
p1 <- ggplot(dg, aes(tgl)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = pal[2], alpha = 0.22) +
  geom_line(aes(y = r), colour = pal[1], linewidth = 0.3) +
  labs(title = "Daily returns with +/-2 conditional SD band", x = NULL, y = "cpo_price_change")
# 2) volatilitas bersyarat
s_tb <- sqrt(cf["omega"] / (1 - pers))
p2 <- ggplot(dg, aes(tgl, sig)) +
  geom_hline(yintercept = s_tb, linetype = "dashed", colour = warna$redup) +
  geom_area(fill = pal[3], alpha = 0.2) +
  geom_line(colour = pal[3], linewidth = 0.6) +
  labs(title = "Conditional volatility (dashed = long-run)", x = NULL, y = "Sigma (%)")
# 3) QQ residual terstandar terhadap distribusi t terstandar
dq <- data.frame(teori = qdist("std", ppoints(length(zres)), mu = 0, sigma = 1, shape = cf["shape"]),
                 sampel = sort(zres))
p3 <- ggplot(dq, aes(teori, sampel)) +
  geom_abline(slope = 1, intercept = 0, colour = pal[2], linewidth = 0.9) +
  geom_point(colour = pal[4], alpha = 0.5, size = 1.2) +
  labs(title = sprintf("Q-Q vs Student-t (df = %.1f)", cf["shape"]), x = "Theoretical quantile", y = "Standardized residual")
p_hasil <- if (tata == 1) p1 / (p2 | p3) else (p1 / p2) | p3

print(p_hasil)