rm(list = ls()) # bersihkan objek lama (mis. dari Nomor 1) agar tidak tercampur # install.packages(c(“vars”, “urca”, “readxl”, “tseries”)) # sekali saja library(vars) # VAR, VARselect, vec2var, irf, fevd, diagnostik library(urca) # uji akar unit & uji kointegrasi Johansen library(readxl) # membaca file Excel
FOLDER <- file.path(Sys.getenv(ifelse(.Platform\(OS.type == "windows", "USERPROFILE", "HOME")), "Downloads") # <-- GANTI jika folder lain FILE_DATA <- file.path(FOLDER, "data_uas_no_2_arw_rapi.xlsx") # <-- GANTI jika nama beda if (!file.exists(FILE_DATA)) { # cari otomatis nama yang mirip (mis. "... (1).xlsx") cand <- list.files(FOLDER, pattern = "no[ _]2[ _]arw.*rapi.*\\.xlsx?\)“, ignore.case = TRUE, full.names = TRUE) if (length(cand) > 0) FILE_DATA <- cand[1] } if (!file.exists(FILE_DATA)) { # masih tidak ketemu -> pilih manual lewat jendela message(”File tidak ditemukan di: “, FILE_DATA,”pilih file secara manual.”) FILE_DATA <- file.choose() } cat(“File data:”, FILE_DATA, “”) SHEET <- “data”
SAMPEL <- “15thn” # “15thn” = Sep 2011 - Agu 2026 (180 obs); “semua” = Jan 2011 - Agu 2026 LAG_MAX <- 12 # lag maksimum yang diperiksa oleh VARselect KRITERIA <- “HQ(n)” # kriteria pemilihan lag: “AIC(n)”, “HQ(n)”, “SC(n)”, “FPE(n)” ECDET <- “const” # uji Johansen: “const” (konstanta di relasi kointegrasi), “trend”, “none” TARIF_UJI <- “5pct” # tingkat signifikansi nilai kritis: “10pct”, “5pct”, “1pct”
DATA_INFLASI <- “otomatis”
URUTAN_DASAR <- c(“lkurs”, “bi_rate”, “INFL”) # <– GANTI bila ingin urutan lain MODEL_PAKSA <- NULL # NULL = otomatis; atau “VECM”, “VAR_diff”, “VAR_level” HORIZON <- 24 # panjang IRF/FEVD (bulan) RUNS_BOOT <- 500 # jumlah replikasi bootstrap untuk selang kepercayaan IRF HORIZON_RAMAL <- 12 # bulan peramalan (Langkah 11)
bintang <- function(p) ifelse(is.na(p), ““, ifelse(p < 0.001,”*“, ifelse(p < 0.01,””, ifelse(p < 0.05, “”, ifelse(p < 0.1, ”.”, ””))))) statdes <- function(x) { x <- as.numeric(x); n <- length(x); m <- mean(x); s0 <- sqrt(mean((x - m)^2)) S <- mean((x - m)^3) / s0^3; K <- mean((x - m)^4) / s0^4 c(Rata2 = m, Median = median(x), Min = min(x), Maks = max(x), SimpBaku = sd(x), Skewness = S, Kurtosis = K, Obs = n) } # Uji akar unit satu deret: ADF (H0: akar unit), PP (H0: akar unit), KPSS (H0: stasioner) uji_akar <- function(x, tipe = c(”drift”, ”trend”)) { tipe <- match.arg(tipe) adf <- ur.df(x, type = tipe, lags = 12, selectlags = ”AIC”) pp <- tryCatch(ur.pp(x, type = ”Z-tau”, model = ifelse(tipe == ”drift”, ”constant”, ”trend”), lags = ”long”), error = function(e) NULL) kp <- ur.kpss(x, type = ifelse(tipe == ”drift”, ”mu”, ”tau”), lags = ”long”) data.frame( ADF_stat = round(adf@teststat[1, 1], 3), ADF_krit = adf@cval[1, ”5pct”], PP_stat = if (is.null(pp)) NA else round(pp@teststat, 3), PP_krit = if (is.null(pp)) NA else pp@cval[1, ”5pct”], KPSS_stat = round(kp@teststat, 3), KPSS_krit = kp@cval[1, ”5pct”], Stasioner_ADF = adf@teststat[1, 1] < adf@cval[1, ”5pct”], Stasioner_KPSS = kp@teststat < kp@cval[1, ”5pct”]) } # Orde integrasi: 0 bila level stasioner (ADF menolak akar unit), 1 bila selisih pertama stasioner orde_int <- function(x, tipe_level) { l <- uji_akar(x, tipe_level) if (l\(Stasioner_ADF) return(0) d1 <- uji_akar(diff(x), "drift") if (d1\)Stasioner_ADF) return(1) 2 } # Tentukan rank kointegrasi secara berurutan (r = 0, r <= 1, …) pilih_rank <- function(jo, tingkat = ”5pct”) { ts_ <- jo@teststat; cv <- jo@cval[, tingkat]; n <- length(ts_); rank <- 0 for (i in n:1) { if (ts_[i] > cv[i]) rank <- rank + 1 else break } rank } # Akar-akar matriks companion dari VAR level (untuk cek stabilitas / akar satuan pada VECM) akar_companion <- function(A_list) { k <- ncol(A_list[[1]]); p <- length(A_list) comp <- matrix(0, k p, k * p) comp[1:k, ] <- do.call(cbind, A_list) if (p > 1) comp[(k + 1):(k * p), 1:(k * (p - 1))] <- diag(k * (p - 1)) sort(Mod(eigen(comp)\(values), decreasing = TRUE) } # Ambil tabel IRF (nilai, batas bawah, batas atas) untuk satu pasangan impuls-respons tabel_irf <- function(ir, impuls, respons, hor = c(1, 3, 6, 12, 24)) { hor <- hor[hor <= HORIZON] idx <- hor + 1 # baris 1 = periode 0 data.frame(Bulan = hor, IRF = round(ir\)irf[[impuls]][idx, respons], 5), Bawah = round(ir\(Lower[[impuls]][idx, respons], 5), Atas = round(ir\)Upper[[impuls]][idx, respons], 5), Signifikan = with(ir, (Lower[[impuls]][idx, respons] > 0) | (Upper[[impuls]][idx, respons] < 0))) } # Grafik IRF 3x3 (baris = respons, kolom = impuls) dengan selang kepercayaan plot_irf_grid <- function(ir, nama, judul) { k <- length(nama); h <- 0:HORIZON par(mfrow = c(k, k), mar = c(3, 3.5, 2.5, 1), oma = c(0, 0, 2, 0)) for (r in nama) for (i in nama) { y <- ir\(irf[[i]][, r]; lo <- ir\)Lower[[i]][, r]; up <- ir$Upper[[i]][, r] plot(h, y, type =”n”, ylim = range(c(lo, up, 0)), xlab = ““, ylab =”“, main = paste0(”Shock “, i,” -> “, r), cex.main = 0.85) polygon(c(h, rev(h)), c(lo, rev(up)), col =”grey85”, border = NA) abline(h = 0, col = “grey40”, lty = 2); lines(h, y, lwd = 2, col = “firebrick”) } mtext(judul, outer = TRUE, font = 2) par(mfrow = c(1, 1), oma = c(0, 0, 0, 0)) }
dat <- read_excel(FILE_DATA, sheet = SHEET) dat\(tanggal <- as.Date(dat\)tanggal) head(dat) # CEK: tanggal, tahun, bulan, inflasi_mtm, bi_rate, kurs_usd, sampel_15thn if (SAMPEL == “15thn”) dat <- subset(dat, sampel_15thn == 1) dat <- dat[order(dat\(tanggal), ] stopifnot(!anyNA(dat[, c("inflasi_mtm", "bi_rate", "kurs_usd")])) cat("Sampel:", nrow(dat), "bulan |", format(min(dat\)tanggal)), “s.d.”, format(max(dat$tanggal)), “”) # HASIL YANG SEHARUSNYA MUNCUL (SAMPEL = “15thn”): # Sampel: 180 bulan | 2011-09-01 s.d. 2026-08-01 print(round(sapply(dat[, c(“inflasi_mtm”, “bi_rate”, “kurs_usd”)], statdes), 3))
awal <- c(as.integer(format(min(dat\(tanggal), "%Y")), as.integer(format(min(dat\)tanggal), “%m”))) mk <- function(x) ts(x, start = awal, frequency = 12) lkurs <- mk(log(dat\(kurs_usd)) bi_rate <- mk(dat\)bi_rate) infl <- mk(dat\(inflasi_mtm) lihk <- mk(log(100) + cumsum(log(1 + dat\)inflasi_mtm / 100)))
par(mfrow = c(4, 1), mar = c(3, 4, 2, 1)) plot(infl, main = “Inflasi bulanan (m-to-m, %)”, ylab = “%”); abline(h = 0, col = “grey60”) plot(bi_rate, main = “BI Rate (%)”, ylab = “%”) plot(lkurs, main = “ln Kurs USD/IDR”, ylab = “ln Rp”) plot(lihk, main = “ln Indeks harga (akumulasi inflasi)”, ylab = “ln”) par(mfrow = c(1, 1))
seri <- list(inflasi = infl, bi_rate = bi_rate, lkurs = lkurs, lihk = lihk) tipe <- c(inflasi = “drift”, bi_rate = “drift”, lkurs = “trend”, lihk = “trend”)
tab_level <- do.call(rbind, lapply(names(seri), function(n) data.frame(Variabel = n, Model = tipe[[n]], uji_akar(seri[[n]], tipe[[n]]), row.names = NULL))) cat(“=== Uji akar unit pada LEVEL ===”); print(tab_level, row.names = FALSE) tab_diff <- do.call(rbind, lapply(names(seri), function(n) data.frame(Variabel = paste0(“d(”, n, “)”), Model = “drift”, uji_akar(diff(seri[[n]]), “drift”), row.names = NULL))) cat(“=== Uji akar unit pada SELISIH PERTAMA ===”); print(tab_diff, row.names = FALSE)
orde <- sapply(names(seri), function(n) orde_int(seri[[n]], tipe[[n]])) cat(“integrasi (0 = stasioner di level, 1 = stasioner di selisih pertama):”); print(orde) # BACA: statistik ADF/PP < nilai kritis -> tolak akar unit (stasioner). # KPSS < nilai kritis -> tidak tolak stasioner. Bila ADF & KPSS sejalan, kesimpulannya kuat. # Kointegrasi Johansen hanya bermakna bila SEMUA variabel dalam sistem I(1).
bentrok <- tab_level\(Variabel[tab_level\)Stasioner_ADF != tab_level$Stasioner_KPSS] if (length(bentrok) > 0) cat(“: ADF dan KPSS memberi kesimpulan BERBEDA untuk:”, paste(bentrok, collapse = “,”), “(bukti campuran; jelaskan di laporan, mis. akibat lonjakan inflasi/outlier).”)
pakai_indeks <- switch(DATA_INFLASI, otomatis = (orde[[“inflasi”]] == 0), laju = FALSE, indeks = TRUE) NAMA_INFL <- if (pakai_indeks) “lihk” else “inflasi” cat(“>>> Variabel inflasi dalam model:”, if (pakai_indeks) “ln INDEKS HARGA (akumulasi inflasi), karena laju inflasi stasioner I(0)” else “LAJU inflasi (I(1) atau dipaksa)”, “”) URUTAN <- sub(“INFL”, NAMA_INFL, URUTAN_DASAR) Y <- cbind(lkurs, bi_rate, get(NAMA_INFL)); colnames(Y) <- c(“lkurs”, “bi_rate”, NAMA_INFL) Y <- Y[, URUTAN] n_var <- ncol(Y) cat(“Urutan variabel (Cholesky):”, paste(URUTAN, collapse = ” -> “),”“) if (!all(orde[c(if (pakai_indeks) ”lihk” else ”inflasi”, ”bi_rate”, ”lkurs”)] == 1)) warning(”Tidak semua variabel I(1). Interpretasikan kointegrasi dengan hati-hati (lihat catatan orde integrasi).”)
sel <- VARselect(Y, lag.max = LAG_MAX, type = “const”) print(sel\(selection) print(round(sel\)criteria, 4)) p_opt <- unname(sel$selection[KRITERIA]) K <- max(2, p_opt) # ca.jo membutuhkan K >= 2 cat(“terpilih (”, KRITERIA, “) p =”, p_opt, “-> K untuk uji Johansen =”, K, “”) # BACA: pilih lag dengan nilai kriteria TERKECIL. HQ/SC cenderung memilih lag lebih hemat, # AIC lebih longgar. Bila diagnostik (Langkah 7) menunjukkan autokorelasi, naikkan lag.
jo_trace <- ca.jo(Y, type = “trace”, ecdet = ECDET, K = K, spec = “transitory”) jo_eigen <- ca.jo(Y, type = “eigen”, ecdet = ECDET, K = K, spec = “transitory”) cat(“=== Uji Johansen (trace) ===”); print(summary(jo_trace)) cat(“=== Uji Johansen (max-eigen) ===”); print(summary(jo_eigen))
tab_joh <- function(jo, nama) { n <- length(jo@teststat) data.frame(Uji = nama, H0 = paste0(“r”, ifelse(seq_len(n) == n, “= 0”, paste0(“<=”, n - seq_len(n)))), Statistik = round(jo@teststat, 3), Kritis_5pct = jo@cval[, “5pct”], Keputusan = ifelse(jo@teststat > jo@cval[, “5pct”], “Tolak H0”, “Gagal tolak H0”)) } print(rbind(tab_joh(jo_trace, “Trace”), tab_joh(jo_eigen, “Max-eigen”)), row.names = FALSE) r_trace <- pilih_rank(jo_trace, TARIF_UJI); r_eigen <- pilih_rank(jo_eigen, TARIF_UJI) cat(“kointegrasi: trace =”, r_trace, “| max-eigen =”, r_eigen, “”) r <- r_trace # keputusan utama memakai uji trace # BACA: r = 0 -> tidak ada kointegrasi (pakai VAR pada selisih pertama). # 0 < r < jumlah variabel -> ada kointegrasi (pakai VECM dengan r relasi jangka panjang). # r = jumlah variabel -> semua variabel stasioner (pakai VAR level). # Bila trace dan max-eigen berbeda, jelaskan di laporan; trace dipakai sebagai keputusan utama.
JENIS <- if (!is.null(MODEL_PAKSA)) MODEL_PAKSA else if (r == 0) “VAR_diff” else if (r >= n_var) “VAR_level” else “VECM” cat(“>>> JENIS MODEL:”, JENIS, “”)
if (JENIS == “VECM”) { r_pakai <- max(1, min(r, n_var - 1)) vecm <- cajorls(jo_trace, r = r_pakai) # VECM terestriksi (OLS dua tahap) model_akhir <- vec2var(jo_trace, r = r_pakai) # representasi VAR level (untuk IRF, FEVD, diagnostik)
# (a) Hubungan JANGKA PANJANG (vektor kointegrasi beta, dinormalkan pada variabel pertama) beta <- vecm\(beta cat("\n=== Vektor kointegrasi (beta), dinormalkan pada", URUTAN[1], "===\n"); print(round(beta, 5)) b <- beta[, 1]; nm_b <- sub("\\.l[0-9]+\)“,”“, names(b)) cat(”jangka panjang (ECT_t = 0):“) cat(sprintf(” %s = %s“, URUTAN[1], paste(sprintf(”%+.4f %s”, -b[-1][seq_len(n_var - 1)], nm_b[2:n_var]), collapse = ” “))) if (length(b) > n_var) cat(sprintf(” konstanta: %+.4f“, -b[n_var + 1])) # (b) Hubungan JANGKA PENDEK & kecepatan penyesuaian (alpha) pada tiap persamaan cat(”=== Persamaan VECM: d(Y) ~ ECT(-1) + d(Y)(-1..) ===“) ringkas <- summary(vecm\(rlm) for (nm in names(ringkas)) { cat("\n---", nm, "---\n"); print(round(ringkas[[nm]]\)coefficients, 5)) } alpha <- sapply(names(ringkas), function(nm) ringkas[[nm]]$coefficients[”ect1”, c(”Estimate”, ”t value”, ”Pr(>|t|)”)]) cat(”=== Koefisien koreksi kesalahan (alpha) ===“); print(round(t(alpha), 5))
} else if (JENIS == “VAR_diff”) { dY <- diff(Y) model_akhir <- VAR(dY, p = max(1, K - 1), type = “const”) cat(“=== VAR pada selisih pertama, p =”, max(1, K - 1), “===”); print(summary(model_akhir)) } else { model_akhir <- VAR(Y, p = p_opt, type = “const”) cat(“=== VAR level, p =”, p_opt, “===”); print(summary(model_akhir)) }
cat(“=== Autokorelasi residual (Portmanteau) ===”) print(serial.test(model_akhir, lags.pt = 16, type = “PT.asymptotic”)) cat(“=== Autokorelasi residual (Breusch-Godfrey, lag 6) ===”) print(serial.test(model_akhir, lags.bg = 6, type = “BG”)) cat(“=== Normalitas residual (Jarque-Bera multivariat) ===”) print(normality.test(model_akhir, multivariate.only = TRUE)) cat(“=== Efek ARCH pada residual (multivariat, lag 5) ===”) print(arch.test(model_akhir, lags.multi = 5, multivariate.only = TRUE))
cat(“=== Stabilitas: akar karakteristik ===”) if (inherits(model_akhir, “varest”)) { print(round(roots(model_akhir), 4)) cat(“Semua akar < 1 ->”, ifelse(all(roots(model_akhir) < 1), “VAR stabil”, “VAR TIDAK stabil”), “”) } else { ak <- akar_companion(model_akhir$A) print(round(ak, 4)) cat(“Pada VECM dengan rank r =”, r, “, harus ada tepat”, n_var - r, “akar satuan (= 1) dan sisanya < 1. Akar yang mendekati 1:”, sum(ak > 0.9999), “”) } # BACA: Portmanteau/BG p-value > 0,05 -> tidak ada autokorelasi (bagus). Bila < 0,05, tambah lag. # JB p-value < 0,05 -> residual tidak normal (umum pada data keuangan; catat sebagai # keterbatasan, IRF bootstrap tetap dapat dipakai). ARCH p-value < 0,05 -> ada heteroskedastisitas.
var_g <- VAR(diff(Y), p = max(1, K - 1), type = “const”) tab_g <- do.call(rbind, lapply(URUTAN, function(v) { g <- causality(var_g, cause = v)\(Granger data.frame(Penyebab = v, F_stat = round(unname(g\)statistic), 3), df = paste(g\(parameter, collapse = ","), p_value = round(g\)p.value, 4), Sig = bintang(g$p.value)) })) cat(“=== Uji Granger (VAR selisih pertama): apakah variabel menyebabkan variabel lain? ===”) print(tab_g, row.names = FALSE)
if (JENIS == “VECM”) { cat(“=== Uji eksogenitas lemah (H0: alpha_j = 0, variabel j tidak menyesuaikan) ===”) for (j in seq_len(n_var)) { A <- diag(n_var)[, -j, drop = FALSE] res <- tryCatch(alrtest(jo_trace, A = A, r = r_pakai), error = function(e) NULL) if (!is.null(res)) cat(sprintf(” %-8s : LR = %.3f, p-value = %.4f %s“, URUTAN[j], res@teststat[1], res@pval[1], bintang(res@pval[1]))) } }
set.seed(123) ir <- irf(model_akhir, n.ahead = HORIZON, ortho = TRUE, cumulative = FALSE, boot = TRUE, ci = 0.95, runs = RUNS_BOOT) jud <- paste0(“IRF ortogonal -”, JENIS, ” (selang 95% bootstrap)“) plot_irf_grid(ir, URUTAN, jud)
tab_ir_all <- do.call(rbind, lapply(URUTAN, function(i) do.call(rbind, lapply(URUTAN, function(rr) { if (i == rr) return(NULL) cbind(Impuls = i, Respons = rr, tabel_irf(ir, i, rr)) })))) print(tab_ir_all, row.names = FALSE)
if (JENIS == “VAR_diff”) { ir_cum <- irf(model_akhir, n.ahead = HORIZON, ortho = TRUE, cumulative = TRUE, boot = TRUE, ci = 0.95, runs = RUNS_BOOT) plot_irf_grid(ir_cum, URUTAN, “IRF kumulatif (respons tingkat) - VAR selisih pertama”) }
fv <- fevd(model_akhir, n.ahead = HORIZON) hor_f <- c(1, 6, 12, 24); hor_f <- hor_f[hor_f <= HORIZON] for (v in URUTAN) { cat(“variabel”, v, “(persen ragam galat ramalan yang dijelaskan oleh shock masing-masing):”) print(round(100 * fv[[v]][hor_f, , drop = FALSE], 2)) } .
fc <- predict(model_akhir, n.ahead = HORIZON_RAMAL, ci = 0.95) for (v in names(fc\(fcst)) { cat("\nRamalan", v, "\n"); print(round(head(fc\)fcst[[v]], HORIZON_RAMAL), 4)) } if (JENIS != “VAR_diff” && “lkurs” %in% names(fc\(fcst)) cat("\nRamalan kurs (Rp/USD), 3 bulan pertama:", round(exp(fc\)fcst$lkurs[1:3, “fcst”])), “”) # CATATAN: pada VAR selisih pertama, ramalan adalah perubahan (selisih), bukan tingkat.
cat(“================ RINGKASAN INTERPRETASI ================”) cat(“Sampel :”, format(min(dat\(tanggal)), "s.d.", format(max(dat\)tanggal)), “(”, nrow(dat), “bulan )”) cat(“Variabel :”, paste(URUTAN, collapse = “,”), “”) cat(“Orde integrasi :”, paste(names(orde), orde, sep = “=”, collapse = “,”), “”) cat(“Lag (”, KRITERIA, “) :”, p_opt, “| K Johansen =”, K, “”) cat(“Kointegrasi : rank trace =”, r_trace, “, max-eigen =”, r_eigen, ifelse(r_trace > 0 & r_trace < n_var, “-> ADA kointegrasi (hubungan jangka panjang)”, ifelse(r_trace == 0, “-> TIDAK ada kointegrasi”, “-> semua variabel stasioner”)), “”) cat(“Model :”, JENIS, “”) if (JENIS == “VECM”) { a1 <- t(alpha)[1, ] cat(“ECT (alpha) pers.”, URUTAN[1], “:”, round(a1[1], 4), “p =”, round(a1[3], 4), bintang(a1[3]), ifelse(a1[1] < 0 & a1[3] < 0.05, “-> menyesuaikan menuju keseimbangan”, “-> penyesuaian tidak nyata”), “”) } # IRF: respons signifikan (selang 95% tidak memuat nol) pada horizon mana pun sig <- subset(tab_ir_all, Signifikan) if (nrow(sig) > 0) { agr <- aggregate(cbind(Bulan_signifikan = Bulan) ~ Impuls + Respons, sig, function(z) paste(z, collapse = “,”)) cat(“IRF signifikan (pasangan impuls -> respons, bulan signifikan pada horizon yang ditampilkan):”) print(agr, row.names = FALSE) } else cat(“IRF: tidak ada respons signifikan pada horizon yang ditampilkan.”) cat(“========================================================”)
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.