Hubungan jangka panjang produktivitas tenaga kerja, upah riil, dan belanja pendidikan: VECM

# Contoh kode R: Johansen Cointegration and VECM
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(urca)
  library(vars)
})
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 I(1) yang berkointegrasi (tren stokastik bersama) ----
nt <- 120
w1 <- cumsum(rnorm(nt, 0, 0.6))                                  # tren stokastik (random walk) pertama
w2 <- cumsum(rnorm(nt, 0, 0.6))                                  # tren stokastik kedua
b1 <- runif(1, 0.6, 1.1); b2 <- runif(1, 0.3, 0.7) * sample(c(-1, 1), 1)
u <- as.numeric(arima.sim(list(ar = runif(1, 0.2, 0.5)), n = nt, sd = runif(1, 0.45, 0.8)))  # galat keseimbangan stasioner
zx1 <- w1 + rnorm(nt, 0, 0.1)
zx2 <- w2 + rnorm(nt, 0, 0.1)
zy <- b1 * zx1 + b2 * zx2 + u                                   # hubungan jangka panjang Y = b1 X1 + b2 X2 + u
sk <- function(v) (v - mean(v)) / sd(v)
d <- data.frame(periode = seq(as.Date("2015-01-01"), by = "month", length.out = nt),
                labour_productivity = round(pmax(100 + 15 * sk(zy), 0), 3),
                real_wage = round(pmax(50 + 10 * sk(zx1), 0), 3),
                education_spending = round(pmax(30 + 6 * sk(zx2), 0), 3))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Uji akar unit ADF: level vs diferensi pertama ----
dv <- d[, c("labour_productivity", "real_wage", "education_spending")]
adf_tau <- function(v, tipe) { u_ <- ur.df(v, type = tipe, lags = 4, selectlags = "AIC"); c(u_@teststat[1], u_@cval[1, 2]) }
adf <- t(sapply(dv, function(v) c(adf_tau(v, "drift"), adf_tau(diff(v), "drift"))))
colnames(adf) <- c("tau_level", "kritis5_level", "tau_diff", "kritis5_diff")
print(round(adf, 3))                                            # tolak H0 akar unit bila tau < nilai kritis
i1 <- sum(adf[, "tau_level"] > adf[, "kritis5_level"] & adf[, "tau_diff"] < adf[, "kritis5_diff"])

# ---- Uji kointegrasi Johansen (trace) ----
cj <- ca.jo(dv, type = "trace", ecdet = "const", K = 2, spec = "transitory")
print(summary(cj))
stat <- rev(cj@teststat); krit <- rev(cj@cval[, 2])             # urutan r = 0, r <= 1, r <= 2
rank_uji <- 0
while (rank_uji < 3 && stat[rank_uji + 1] > krit[rank_uji + 1]) rank_uji <- rank_uji + 1
rank_pakai <- min(max(rank_uji, 1), 2)                          # cadangan: VECM butuh 1 <= r <= 2

# ---- VECM: vektor kointegrasi (beta) dan koefisien penyesuaian (alpha) ----
vecm <- cajorls(cj, r = 1)                                      # vektor pertama, dinormalisasi pada Y
print(vecm$rlm); print(vecm$beta)
beta <- vecm$beta[, 1]
alpha <- coef(vecm$rlm)["ect1", ]
lr_x1 <- -unname(beta[2]); lr_x2 <- -unname(beta[3])            # Y = lr_x1 X1 + lr_x2 X2 + konstanta
ect <- as.numeric(as.matrix(cbind(dv, 1)) %*% beta)             # galat koreksi (ECT)
p_alpha <- summary(vecm$rlm)[[1]]$coefficients["ect1", 4]       # uji t alpha pada persamaan dY

# ---- IRF dari representasi VAR level ----
vv <- vec2var(cj, r = rank_pakai)
ir <- irf(vv, impulse = c("real_wage", "education_spending"), response = "labour_productivity", n.ahead = 20, boot = TRUE, runs = 100)


# ---- Grafik ----
nm3 <- make.unique(c("labour_productivity", "real_wage", "education_spending"), sep = " ")
dl <- data.frame(periode = rep(d$periode, 3), nilai = c(d$labour_productivity, d$real_wage, d$education_spending),
                 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 = "Levels of the three series (I(1))", x = NULL, y = NULL)
de <- data.frame(periode = d$periode, ect = ect - mean(ect))
p2 <- ggplot(de, aes(periode, ect)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_hline(yintercept = c(-2, 2) * sd(de$ect), linetype = "dashed", colour = pal[2]) +
  geom_area(fill = pal[4], alpha = 0.25) +
  geom_line(colour = pal[4], linewidth = 0.7) +
  labs(title = "Error-correction term (deviation from LR)", x = NULL, y = "ECT (centred)")
kej <- c("real_wage", "education_spending")
lab_s <- make.unique(c("Shock: real_wage", "Shock: education_spending"), sep = " ")
di <- do.call(rbind, lapply(1:2, function(i) data.frame(
  h = 0:20, 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)
p3 <- ggplot(di, aes(h, irf)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_ribbon(aes(ymin = lo, ymax = hi, fill = guncangan), alpha = 0.2, 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]) +
  labs(title = "IRF of labour_productivity from VECM (vec2var)", x = "Horizon", y = "Response")
p_hasil <- if (tata == 1) p1 | (p2 / p3) else p1 / (p2 | p3)

print(p_hasil)