Interaksi belanja iklan, lalu lintas situs, dan penjualan daring kuartalan: VAR dan uji Granger

# Contoh kode R: VAR and Granger Causality
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(vars)
  library(car)
})
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: tiga deret bulanan dari proses VAR(2) yang stasioner ----
nt <- 120
# matriks koefisien lag 1 dan lag 2; diulang sampai prosesnya stabil (semua akar < 0.9)
repeat {
  A1 <- diag(runif(3, 0.3, 0.55))
  A1[1, 2] <- runif(1, 0.25, 0.45)                              # X1 -> Y
  A1[1, 3] <- runif(1, 0.1, 0.35) * sample(c(-1, 1), 1)        # X2 -> Y
  A1[2, 3] <- runif(1, -0.15, 0.15)
  A1[3, 2] <- runif(1, -0.1, 0.2)
  A2 <- diag(runif(3, -0.15, 0.15))
  A2[1, 2] <- runif(1, 0.05, 0.2)
  komp <- rbind(cbind(A1, A2), cbind(diag(3), matrix(0, 3, 3)))  # matriks pendamping (companion)
  if (max(Mod(eigen(komp)$values)) < 0.9) break
}
z <- matrix(0, nt + 50, 3)
for (i in 3:(nt + 50)) z[i, ] <- A1 %*% z[i - 1, ] + A2 %*% z[i - 2, ] + rnorm(3, 0, c(0.6, 1, 1))
z <- scale(z[-(1:50), ])                                        # buang periode pemanasan (burn-in)
d <- data.frame(periode = seq(as.Date("2015-01-01"), by = "month", length.out = nt),
                online_sales = round(pmax(100 + 15 * z[, 1], 0), 3),
                advertising_spend = round(pmax(50 + 10 * z[, 2], 0), 3),
                site_traffic = round(pmax(30 + 6 * z[, 3], 0), 3))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: Vector Autoregression (VAR) ----
dv <- d[, c("online_sales", "advertising_spend", "site_traffic")]
sel <- VARselect(dv, lag.max = 6, type = "const")               # pemilihan panjang lag
print(sel$selection)
p_lag <- max(1, as.integer(sel$selection[["AIC(n)"]]))
fit <- do.call(VAR, list(y = dv, p = p_lag, type = "const"))  # do.call: nilai lag tersimpan di objek
print(summary(fit))

# ---- Diagnostik: stabilitas dan autokorelasi residual ----
akar <- roots(fit)                                              # modulus akar; < 1 berarti VAR stabil
st <- serial.test(fit, lags.pt = 12, type = "PT.asymptotic")    # uji Portmanteau
print(akar); print(st)

# ---- Kausalitas Granger ----
# uji Wald pada persamaan online_sales: apakah semua lag X1 (atau X2) bersama-sama nol?
eq_y <- fit$varresult$online_sales
gr_x1 <- car::linearHypothesis(eq_y, paste0("advertising_spend.l", 1:p_lag, " = 0"))
gr_x2 <- car::linearHypothesis(eq_y, paste0("site_traffic.l", 1:p_lag, " = 0"))
p_x1 <- gr_x1[2, "Pr(>F)"]; p_x2 <- gr_x2[2, "Pr(>F)"]
print(gr_x1); print(gr_x2)
print(causality(fit, cause = "advertising_spend")$Granger)                 # versi blok (X1 -> semua variabel lain)

# ---- Impulse response function (IRF) dan dekomposisi varians (FEVD) ----
ir <- irf(fit, impulse = c("advertising_spend", "site_traffic"), response = "online_sales", n.ahead = 10,
          ortho = TRUE, boot = TRUE, runs = 200)
fv <- fevd(fit, n.ahead = 10)
print(round(fv$online_sales, 3))


# ---- Grafik ----
nm3 <- make.unique(c("online_sales", "advertising_spend", "site_traffic"), sep = " ")
dl <- data.frame(periode = rep(d$periode, 3), nilai = c(d$online_sales, d$advertising_spend, d$site_traffic),
                 seri = factor(rep(nm3, each = nt), levels = nm3))
p1 <- ggplot(dl, aes(periode, nilai, colour = seri)) +
  geom_line(linewidth = 0.7, show.legend = FALSE) +
  facet_wrap(~seri, ncol = if (tata == 1) 1 else 3, scales = "free_y") +
  scale_colour_manual(values = pal[1:3]) +
  labs(title = "The three series", x = NULL, y = NULL)
kej <- c("advertising_spend", "site_traffic")
lab_s <- make.unique(c("Shock: advertising_spend", "Shock: site_traffic"), sep = " ")
di <- do.call(rbind, lapply(1:2, function(i) data.frame(
  h = 0:10, irf = ir$irf[[kej[i]]][, 1], lo = ir$Lower[[kej[i]]][, 1], hi = ir$Upper[[kej[i]]][, 1],
  guncangan = lab_s[i])))
di$guncangan <- factor(di$guncangan, levels = lab_s)
p2 <- ggplot(di, aes(h, irf)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_ribbon(aes(ymin = lo, ymax = hi, fill = guncangan), alpha = 0.22, show.legend = FALSE) +
  geom_line(aes(colour = guncangan), linewidth = 1, show.legend = FALSE) +
  facet_wrap(~guncangan, ncol = 2) +
  scale_colour_manual(values = pal[2:3]) + scale_fill_manual(values = pal[2:3]) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(title = "Orthogonal IRF of online_sales (95% bootstrap)", x = "Horizon", y = "Response")
nm_s <- make.unique(c("online_sales", "advertising_spend", "site_traffic"), sep = " ")
dfv <- data.frame(h = rep(1:10, 3), porsi = 100 * c(fv$online_sales),
                  sumber = factor(rep(nm_s, each = 10), levels = nm_s))
p3 <- ggplot(dfv, aes(factor(h), porsi, fill = sumber)) +
  geom_col(width = 0.75, colour = warna$panel, linewidth = 0.3) +
  scale_fill_manual(values = pal[1:3], name = "Shock") +
  labs(title = "FEVD of online_sales", x = "Horizon", y = "Share of variance (%)")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)