Skripsi Hubungan Simultan Utang Luar Negeri & Pertumbuhan Ekonomi

library(readxl)
library(dplyr)
library(AER)
library(lmtest)
library(sandwich)
library(ivreg2r)
library(tseries)
library(car)
library(scales)
library(ggplot2)
library(urca)
library(dplyr)
library(knitr)
library(broom)

# Import data
data1 <- read_excel(
  "D:/STIS/Semester 7/Skripsi/Data/data skripsi.xlsx",
  sheet = 6
)

# Data asli
data_asli <- data.frame(
  LPDB = data1$LGDP,
  RULN = data1$ULN_PERCENTOFGDP,
  FDI = data1$FDI,
  LDEF = data1$LDEF,
  LRES = data1$LRES,
  LEMPLOYMENT = data1$LEMPLOYMENT,
  TOT = data1$TOT,
  LEXCHANGE = data1$LEXCHANGE
)

# Membuat first difference
data <- data.frame(
  d_LPDB = diff(data1$LGDP),
  d_RULN = diff(data1$ULN_PERCENTOFGDP),
  d_FDI = diff(data1$FDI),
  d_LDEF = diff(data1$LDEF),
  d_LRES = diff(data1$LRES),
  d_LEMPLOYMENT = diff(data1$LEMPLOYMENT),
  d_TOT = diff(data1$TOT),
  d_LEXCHANGE = diff(data1$LEXCHANGE)
)

# ============================================================
# UJI STASIONERITAS AUGMENTED DICKEY-FULLER (ADF)
# Spesifikasi: Intercept + SIC/BIC
# Max Lag: 3
# Format Output: Menyerupai EViews
# ============================================================

# ------------------------------------------------------------
# 1. DATA YANG DIGUNAKAN
# ------------------------------------------------------------

data_level <- data_asli
data_diff  <- data


# ------------------------------------------------------------
# 2. FUNGSI P-VALUE MACKINNON (seperti EViews)
# ------------------------------------------------------------

pvalue_mackinnon <- function(stat, n) {

  # Koefisien respons permukaan MacKinnon (1996)
  # Model dengan intercept (drift), tanpa trend
  df_quantile <- function(p, n) {
    if (p == 0.01) {
      b0 <- -3.43035; b1 <- -6.5393; b2 <- -16.786
    } else if (p == 0.05) {
      b0 <- -2.86154; b1 <- -2.8903; b2 <- -4.234
    } else if (p == 0.10) {
      b0 <- -2.56677; b1 <- -1.5384; b2 <- -2.809
    } else {
      return(NA_real_)
    }
    return(b0 + b1 / n + b2 / n^2)
  }

  q01 <- df_quantile(0.01, n)
  q05 <- df_quantile(0.05, n)
  q10 <- df_quantile(0.10, n)

  if (stat <= q01) return(0.0000)
  if (stat <= q05) return(0.01 * (stat - q01) / (q05 - q01))
  if (stat <= q10) return(0.05 * (stat - q05) / (q10 - q05) + 0.05)
  return(0.10 + 0.90 * (stat - q10) / (abs(q10) * 2))
}


# ------------------------------------------------------------
# 3. FUNGSI UJI ADF (FORMAT EVIEWS)
# ------------------------------------------------------------

uji_adf_eviews <- function(dataset, nama_dataset) {

  variabel_numerik <- names(dataset)[
    sapply(dataset, is.numeric)
  ]

  hasil <- lapply(variabel_numerik, function(nama) {

    x <- dataset[[nama]]
    x <- x[is.finite(x)]
    n <- length(x)

    if (n < 10 || length(unique(x)) < 2) {
      return(data.frame(
        Dataset = nama_dataset,
        Variabel = nama,
        N = n,
        `Statistik ADF` = NA_real_,
        `P-value` = NA_real_,
        Kesimpulan = "Tidak dapat diuji",
        check.names = FALSE
      ))
    }

    # ------------------------------------------------------
    # MAX LAG = 3 (fixed, sesuai permintaan)
    # ------------------------------------------------------
    max_lag <- 3
    max_lag <- min(max_lag, floor(n / 2) - 2)
    max_lag <- max(0, max_lag)

    hasil_adf <- tryCatch(
      urca::ur.df(y = x, type = "drift",
                  lags = max_lag, selectlags = "BIC"),
      error = function(e) NULL
    )

    if (is.null(hasil_adf)) {
      return(data.frame(
        Dataset = nama_dataset,
        Variabel = nama,
        N = n,
        `Statistik ADF` = NA_real_,
        `P-value` = NA_real_,
        Kesimpulan = "Uji gagal",
        check.names = FALSE
      ))
    }

    stat <- as.numeric(hasil_adf@teststat[1, "tau2"])
    cval <- hasil_adf@cval["tau2", ]
    pval <- pvalue_mackinnon(stat, n)
    stasioner <- stat < cval["5pct"]

    data.frame(
      Dataset = nama_dataset,
      Variabel = nama,
      N = n,
      `Statistik ADF` = stat,
      `P-value` = pval,
      Kesimpulan = ifelse(stasioner, "Stasioner", "Tidak stasioner"),
      check.names = FALSE
    )
  })

  bind_rows(hasil)
}


# ------------------------------------------------------------
# 4. FUNGSI FORMAT TABEL GAYA EVIEWS
# ------------------------------------------------------------

tabel_eviews <- function(hasil, judul) {

  df_tampil <- hasil %>%
    mutate(
      `P-value` = ifelse(is.na(`P-value`), "NA",
                         formatC(`P-value`, format = "f", digits = 4)),
      `Statistik ADF` = formatC(`Statistik ADF`, format = "f", digits = 6)
    )

  knitr::kable(
    df_tampil,
    col.names = c("Dataset", "Variabel", "Obs",
                  "t-Statistic", "Prob.*", "Kesimpulan"),
    caption = judul,
    align = c("l", "l", "r", "r", "r", "l")
  )
}


# ------------------------------------------------------------
# 5. UJI ADF DATA LEVEL
# ------------------------------------------------------------

hasil_level <- uji_adf_eviews(
  dataset = data_level,
  nama_dataset = "Level"
)

tabel_eviews(
  hasil_level,
  "Augmented Dickey-Fuller Unit Root Test - Level (Intercept)"
)

# ------------------------------------------------------------
# 6. UJI ADF DATA FIRST DIFFERENCE
# ------------------------------------------------------------

hasil_diff <- uji_adf_eviews(
  dataset = data_diff,
  nama_dataset = "First Difference"
)

tabel_eviews(
  hasil_diff,
  "Augmented Dickey-Fuller Unit Root Test - First Difference (Intercept)"
)

lag_d_RULN <- lag(data$d_RULN)
lag_d_LPDB <- lag(data$d_LPDB)

eq_gdp <- ivreg2(
  d_LPDB ~
    d_FDI +
    d_LEMPLOYMENT +
    d_LRES
  |
    d_RULN
  |
    d_LEXCHANGE +
    d_TOT +
    d_LDEF,
  data = data,
  vcov = "iid"
)

summary(eq_gdp)

lag_d_LPDB <- lag(data$d_LPDB)
lag_d_RULN <- lag(data$d_RULN)

eq_uln <- ivreg2(
  d_RULN ~
    lag_d_RULN +
    d_LEXCHANGE +
    d_LDEF +
    d_TOT +
    d_LRES
  |
    d_LPDB
  |
    d_LEMPLOYMENT +
    d_FDI +
    lag_d_LPDB,
  data = data,
  vcov = "iid"
)

summary(eq_uln)

dataa <- read_excel(
  "D:/STIS/Semester 7/Skripsi/Data/data skripsi.xlsx",
  sheet = 3
)

ggplot(dataa, aes(x = LEXCHANGE, y = ULN_PERCENTOFGDP)) +
  geom_point(
    shape = 21, fill = "#5DA5DA", color = "#1F4E79",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#1F4E79",
    fill = "#A8C8E8", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  labs(
    x = "Nilai Tukar (Rupiah/US$)",
    y = "Utang Luar Negeri (% terhadap PDB)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = TOT, y = ULN_PERCENTOFGDP)) +
  geom_point(
    shape = 21, fill = "#5DA5DA", color = "#1F4E79",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#1F4E79",
    fill = "#A8C8E8", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  labs(
    x = "Terms of Trade (2000=100)",
    y = "Utang Luar Negeri (% terhadap PDB)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = LDEF, y = ULN_PERCENTOFGDP)) +
  geom_point(
    shape = 21, fill = "#5DA5DA", color = "#1F4E79",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#1F4E79",
    fill = "#A8C8E8", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(labels = scales::comma) +
  labs(
    x = "Defisit Anggaran (Milyar Rupiah)",
    y = "Utang Luar Negeri (% terhadap PDB)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = LRES, y = ULN_PERCENTOFGDP)) +
  geom_point(
    shape = 21, fill = "#5DA5DA", color = "#1F4E79",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#1F4E79",
    fill = "#A8C8E8", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_continuous(
    labels = label_number(scale = 1e-9, accuracy = 1),
    name = "Cadangan Devisa (Milyar US$)"
  ) +
  labs(y = "Utang Luar Negeri (% terhadap PDB)") +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = FDI, y = LGDP)) +
  geom_point(
    shape = 21, fill = "#7BC96F", color = "#1B5E20",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#2E7D32",
    fill = "#A5D6A7", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(
    labels = label_number(scale = 1e-9, suffix = "", accuracy = 1),
    name = "Pertumbuhan Ekonomi (Milyar US$)"
  ) +
  scale_x_continuous(labels = scales::comma) +
  labs(x = "Foreign Direct Investment (% terhadap PDB)") +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = LEMPLOYMENT, y = LGDP)) +
  geom_point(
    shape = 21, fill = "#7BC96F", color = "#1B5E20",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#2E7D32",
    fill = "#A5D6A7", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(
    labels = label_number(scale = 1e-9, suffix = "", accuracy = 1),
    name = "Pertumbuhan Ekonomi (Milyar US$)"
  ) +
  scale_x_continuous(
    labels = label_number(scale = 1e-6, accuracy = 1),
    name = "Jumlah Angkatan Kerja (Juta Jiwa)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

ggplot(dataa, aes(x = LRES, y = LGDP)) +
  geom_point(
    shape = 21, fill = "#7BC96F", color = "#1B5E20",
    size = 3, alpha = 0.8, stroke = 0.8
  ) +
  geom_smooth(
    method = "lm", se = TRUE, color = "#2E7D32",
    fill = "#A5D6A7", alpha = 0.25, linewidth = 1
  ) +
  scale_y_continuous(
    labels = label_number(scale = 1e-9, suffix = "", accuracy = 1),
    name = "Pertumbuhan Ekonomi (Milyar US$)"
  ) +
  scale_x_continuous(
    labels = label_number(scale = 1e-9, accuracy = 1),
    name = "Cadangan Devisa (Milyar US$)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.title = element_text(face = "bold", color = "#1A1A1A"),
    axis.text = element_text(color = "grey30"),
    panel.grid.major = element_line(color = "grey92", linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "white", color = NA)
  )

# ============================================================
# SUMMARY MODEL IV/2SLS DAN UJI ASUMSI
# ============================================================



# ------------------------------------------------------------
# 1. FUNGSI KOEFISIEN REGRESI
# ------------------------------------------------------------

ambil_koefisien <- function(model, nama_persamaan) {
  hasil <- broom::tidy(model)

  data.frame(
    Persamaan = nama_persamaan,
    Variabel = hasil$term,
    Koefisien = hasil$estimate,
    `Standar Error` = hasil$std.error,
    Prob = hasil$p.value,
    check.names = FALSE
  )
}

koef_gdp <- ambil_koefisien(eq_gdp, "Persamaan GDP")
koef_uln <- ambil_koefisien(eq_uln, "Persamaan ULN")

# ------------------------------------------------------------
# 2. FUNGSI MERAPIKAN DIAGNOSTIK IV
# ------------------------------------------------------------

rapikan_diagnostik <- function(data, nama_persamaan) {

  # Sesuaikan nama kolom
  nama_kolom <- names(data)

  kolom_uji <- intersect(
    c("test", "diagnostic", "name", "Variabel"),
    nama_kolom
  )[1]

  kolom_stat <- intersect(
    c("statistic", "stat", "Koefisien"),
    nama_kolom
  )[1]

  kolom_prob <- intersect(
    c("p_value", "p.value", "pvalue", "Prob"),
    nama_kolom
  )[1]

  if (is.na(kolom_uji) || is.na(kolom_stat)) {
    stop("Periksa struktur output diagnostics() pada ivreg2r.")
  }

  hasil <- data.frame(
    Persamaan = nama_persamaan,
    Variabel = as.character(data[[kolom_uji]]),
    Koefisien = as.numeric(data[[kolom_stat]]),
    `Standar Error` = NA_real_,
    Prob = if (!is.na(kolom_prob)) {
      as.numeric(data[[kolom_prob]])
    } else {
      NA_real_
    },
    check.names = FALSE
  )

  nama_baru <- c(
    underid = "Underidentification test",
    weak_id = "Weak identification test",
    overid = "Overidentification test (Sargan/Hansen)",
    endogeneity = "Endogeneity test",
    anderson_rubin_f = "Anderson-Rubin F test",
    anderson_rubin_chi2 = "Anderson-Rubin Chi-square test",
    stock_wright = "Stock-Wright test"
  )

  hasil %>%
    filter(Variabel %in% names(nama_baru)) %>%
    mutate(Variabel = unname(nama_baru[Variabel]))
}

diag_gdp <- rapikan_diagnostik(
  as.data.frame(ivreg2r::diagnostics(eq_gdp)),
  "Persamaan GDP"
)

diag_uln <- rapikan_diagnostik(
  as.data.frame(ivreg2r::diagnostics(eq_uln)),
  "Persamaan ULN"
)

# ------------------------------------------------------------
# 3. FUNGSI UJI ASUMSI
# ------------------------------------------------------------

baris_uji <- function(persamaan, nama_uji, statistik, prob) {
  data.frame(
    Persamaan = persamaan,
    Variabel = nama_uji,
    Koefisien = as.numeric(statistik),
    `Standar Error` = NA_real_,
    Prob = as.numeric(prob),
    check.names = FALSE
  )
}

# ============================================================
# 4. PERSAMAAN GDP
# ============================================================

mf_gdp <- model.frame(eq_gdp)
res_gdp <- as.numeric(residuals(eq_gdp))

# Normalitas
jb_gdp <- tseries::jarque.bera.test(res_gdp)

# Heteroskedastisitas: Breusch-Pagan
bp_data_gdp <- data.frame(
  res2 = res_gdp^2,
  d_RULN = mf_gdp$d_RULN,
  d_FDI = mf_gdp$d_FDI,
  d_LEMPLOYMENT = mf_gdp$d_LEMPLOYMENT,
  d_LRES = mf_gdp$d_LRES
)

aux_bp_gdp <- lm(
  res2 ~ d_RULN + d_FDI + d_LEMPLOYMENT + d_LRES,
  data = bp_data_gdp
)

LM_BP_gdp <- nobs(aux_bp_gdp) * summary(aux_bp_gdp)$r.squared
p_BP_gdp <- pchisq(LM_BP_gdp, df = 4, lower.tail = FALSE)

# Autokorelasi: Breusch-Godfrey orde 1
bg_gdp <- lmtest::bgtest(eq_gdp, order = 1, type = "Chisq")

# Multikolinearitas: VIF dari OLS bantu
aux_gdp <- lm(
  d_LPDB ~ d_RULN + d_FDI + d_LEMPLOYMENT + d_LRES,
  data = mf_gdp
)

vif_gdp <- car::vif(aux_gdp)

uji_asumsi_gdp <- bind_rows(
  baris_uji(
    "Persamaan GDP", "Normalitas (Jarque-Bera)",
    jb_gdp$statistic, jb_gdp$p.value
  ),
  baris_uji(
    "Persamaan GDP", "Heteroskedastisitas (Breusch-Pagan)",
    LM_BP_gdp, p_BP_gdp
  ),
  baris_uji(
    "Persamaan GDP", "Autokorelasi (Breusch-Godfrey)",
    bg_gdp$statistic, bg_gdp$p.value
  ),
  data.frame(
    Persamaan = "Persamaan GDP",
    Variabel = paste0("VIF: ", names(vif_gdp)),
    Koefisien = as.numeric(vif_gdp),
    `Standar Error` = NA_real_,
    Prob = NA_real_,
    check.names = FALSE
  )
)

# ============================================================
# 5. PERSAMAAN ULN
# ============================================================

mf_uln <- model.frame(eq_uln)
res_uln <- as.numeric(residuals(eq_uln))

# Normalitas
jb_uln <- tseries::jarque.bera.test(res_uln)

# Heteroskedastisitas: Breusch-Pagan
bp_data_uln <- data.frame(
  res2 = res_uln^2,
  lag_d_RULN = mf_uln$lag_d_RULN,
  d_LEXCHANGE = mf_uln$d_LEXCHANGE,
  d_LDEF = mf_uln$d_LDEF,
  d_TOT = mf_uln$d_TOT,
  d_LRES = mf_uln$d_LRES,
  d_LPDB = mf_uln$d_LPDB
)

aux_bp_uln <- lm(
  res2 ~ lag_d_RULN + d_LEXCHANGE + d_LDEF +
    d_TOT + d_LRES + d_LPDB,
  data = bp_data_uln
)

LM_BP_uln <- nobs(aux_bp_uln) * summary(aux_bp_uln)$r.squared
p_BP_uln <- pchisq(LM_BP_uln, df = 6, lower.tail = FALSE)

# Autokorelasi: Breusch-Godfrey orde 1
bg_uln <- lmtest::bgtest(eq_uln, order = 1, type = "Chisq")

# Multikolinearitas: VIF dari OLS bantu
aux_uln <- lm(
  d_RULN ~ d_LPDB + lag_d_RULN + d_LEXCHANGE +
    d_LDEF + d_TOT + d_LRES,
  data = mf_uln
)

vif_uln <- car::vif(aux_uln)

uji_asumsi_uln <- bind_rows(
  baris_uji(
    "Persamaan ULN", "Normalitas (Jarque-Bera)",
    jb_uln$statistic, jb_uln$p.value
  ),
  baris_uji(
    "Persamaan ULN", "Heteroskedastisitas (Breusch-Pagan)",
    LM_BP_uln, p_BP_uln
  ),
  baris_uji(
    "Persamaan ULN", "Autokorelasi (Breusch-Godfrey)",
    bg_uln$statistic, bg_uln$p.value
  ),
  data.frame(
    Persamaan = "Persamaan ULN",
    Variabel = paste0("VIF: ", names(vif_uln)),
    Koefisien = as.numeric(vif_uln),
    `Standar Error` = NA_real_,
    Prob = NA_real_,
    check.names = FALSE
  )
)

# ------------------------------------------------------------
# 6. GABUNGKAN HASIL PERSAMAAN GDP
# ------------------------------------------------------------

summary_gdp <- bind_rows(
  koef_gdp,
  diag_gdp,
  uji_asumsi_gdp
) %>%
  mutate(
    Koefisien = round(Koefisien, 4),
    `Standar Error` = round(`Standar Error`, 4),
    Prob = round(Prob, 4)
  )

# ------------------------------------------------------------
# 7. GABUNGKAN HASIL PERSAMAAN ULN
# ------------------------------------------------------------

summary_uln <- bind_rows(
  koef_uln,
  diag_uln,
  uji_asumsi_uln
) %>%
  mutate(
    Koefisien = round(Koefisien, 4),
    `Standar Error` = round(`Standar Error`, 4),
    Prob = round(Prob, 4)
  )

# ------------------------------------------------------------
# 8. TABEL RINGKASAN PERSAMAAN GDP
# ------------------------------------------------------------

knitr::kable(
  summary_gdp,
  caption = "Ringkasan Model IV/2SLS dan Uji Asumsi Persamaan GDP",
  col.names = c(
    "Persamaan",
    "Variabel/Pengujian",
    "Koefisien/Statistik",
    "Standar Error",
    "Prob."
  ),
  digits = 4,
  na = "-"
)

# ------------------------------------------------------------
# 9. TABEL RINGKASAN PERSAMAAN ULN
# ------------------------------------------------------------

knitr::kable(
  summary_uln,
  caption = "Ringkasan Model IV/2SLS dan Uji Asumsi Persamaan ULN",
  col.names = c(
    "Persamaan",
    "Variabel/Pengujian",
    "Koefisien/Statistik",
    "Standar Error",
    "Prob."
  ),
  digits = 4,
  na = "-"
)