Paquetes y Utilidades

library(readxl)
library(zoo)
library(urca)
library(vars)
library(lmtest)
library(sandwich)
library(knitr)

kbl <- function(x, cap=""){ knitr::kable(x, caption = cap, align = "c") }

1 Carga y armado de series

DATA_FILE <- "DATA NUEVA TESIS.xlsx"     
SHEET     <- "Hoja2"                      

stopifnot(file.exists(DATA_FILE))
df <- read_excel(DATA_FILE, sheet = SHEET)

cols_req <- c("Año/Trimestre","PIB","Exportaciones","Importaciones","Remesas","IED")
stopifnot(all(cols_req %in% names(df)))

df <- df[order(df$`Año/Trimestre`), ]

get_q <- function(qchar){ list(y=as.integer(substr(qchar,1,4)),
                               q=as.integer(substr(qchar,6,6))) }
st <- get_q(df$`Año/Trimestre`[1])

y_pib <- ts(df$PIB,           start=c(st$y, st$q), frequency=4)
y_x   <- ts(df$Exportaciones, start=c(st$y, st$q), frequency=4)
y_m   <- ts(df$Importaciones, start=c(st$y, st$q), frequency=4)
y_rem <- ts(df$Remesas,       start=c(st$y, st$q), frequency=4)
y_ied <- ts(df$IED,           start=c(st$y, st$q), frequency=4)

ln_pib <- log(y_pib); ln_x <- log(y_x); ln_m <- log(y_m)
ln_rem <- log(y_rem); ln_ied <- log(y_ied)

dln_pib <- diff(ln_pib); dln_x <- diff(ln_x); dln_m <- diff(ln_m)
dln_rem <- diff(ln_rem); dln_ied <- diff(ln_ied)

Y  <- cbind(ln_pib, ln_ied, ln_rem, ln_x, ln_m)
colnames(Y) <- c("ln_pib","ln_ied","ln_rem","ln_x","ln_m")
DY <- cbind(dln_pib, dln_ied, dln_rem, dln_x, dln_m)
colnames(DY) <- paste0("d", colnames(Y))

1.1 Gráfica de trayectorias (log)

par(mfrow=c(2,3), mar=c(3,3,2,1))
plot(ln_pib, main="ln(PIB)")
plot(ln_ied, main="ln(IED)")
plot(ln_rem, main="ln(Remesas)")
plot(ln_x,   main="ln(Exportaciones)")
plot(ln_m,   main="ln(Importaciones)")
par(mfrow=c(1,1))

2 Pruebas de raíz unitaria

one_ur <- function(x, nombre){
  adf_L <- tryCatch(summary(ur.df(x, type="trend", selectlags="AIC"))@teststat[1], error=function(e) NA)
  pp_L  <- tryCatch(summary(ur.pp(x, type="Z-tau", model="trend"))@teststat, error=function(e) NA)
  kpssL <- tryCatch(kpss.test(as.numeric(na.omit(x)), null="Level")$p.value, error=function(e) NA)
  kpssT <- tryCatch(kpss.test(as.numeric(na.omit(x)), null="Trend")$p.value, error=function(e) NA)

  dx <- diff(x)
  adf_D <- tryCatch(summary(ur.df(dx, type="drift", selectlags="AIC"))@teststat[1], error=function(e) NA)
  pp_D  <- tryCatch(summary(ur.pp(dx, type="Z-tau", model="constant"))@teststat, error=function(e) NA)

  data.frame(Variable=nombre,
             ADF_nivel=round(adf_L,3), PP_nivel=round(pp_L,3),
             KPSS_Level_p=round(kpssL,3), KPSS_Trend_p=round(kpssT,3),
             ADF_dif=round(adf_D,3), PP_dif=round(pp_D,3))
}

UR_TAB <- rbind(
  one_ur(ln_pib,"ln_pib"),
  one_ur(ln_ied,"ln_ied"),
  one_ur(ln_rem,"ln_rem"),
  one_ur(ln_x,  "ln_x"),
  one_ur(ln_m,  "ln_m")
)

kbl(UR_TAB, "Tabla B. Pruebas de raíz unitaria (estadísticos y p-valores KPSS)")
Tabla B. Pruebas de raíz unitaria (estadísticos y p-valores KPSS)
Variable ADF_nivel PP_nivel KPSS_Level_p KPSS_Trend_p ADF_dif PP_dif
ln_pib -3.931 -3.213 NA NA -14.313 -9.832
ln_ied -4.073 -5.768 NA NA -9.369 -15.595
ln_rem -3.669 -3.111 NA NA -9.230 -14.035
ln_x -4.330 -3.810 NA NA -11.537 -9.456
ln_m -2.940 -2.564 NA NA -5.525 -7.988

3 Selección de rezagos y cointegración

sel <- VARselect(Y, lag.max = 4, type = "const")$selection
K <- if (!is.na(sel["SC(n)"])) sel["SC(n)"] else if (!is.na(sel["HQ(n)"])) sel["HQ(n)"] else 2L
K <- as.integer(K)
K <- max(2L, min(4L, K))

jo <- try(ca.jo(Y, type = "trace", ecdet = "trend", K = K, spec = "transitory"), silent = TRUE)

# Fallback
if (inherits(jo, "try-error")) {
  jo <- try(ca.jo(Y, type = "trace", ecdet = "const", K = K, spec = "transitory"), silent = TRUE)
}

# Fallback 2
if (inherits(jo, "try-error")) {
  ln_open <- log((y_x + y_m)/y_pib)
  Y_alt <- cbind(ln_pib, ln_ied, ln_rem, ln_open)
  colnames(Y_alt) <- c("ln_pib","ln_ied","ln_rem","ln_open")
  sel_alt <- VARselect(Y_alt, lag.max = 4, type = "const")$selection
  K <- if (!is.na(sel_alt["SC(n)"])) sel_alt["SC(n)"] else if (!is.na(sel_alt["HQ(n)"])) sel_alt["HQ(n)"] else 2L
  K <- max(2L, min(4L, as.integer(K)))
  jo <- ca.jo(Y_alt, type = "trace", ecdet = "trend", K = K, spec = "transitory")
  Y <- Y_alt  
}

# Rango r (5%)
pick_r <- function(j) sum(j@teststat > j@cval[,"5pct"])
r <- min(pick_r(jo), ncol(Y) - 1L)

# Tabla de traza (USAR names(), no rownames())
ts_vals <- as.numeric(jo@teststat)        # vector de estadísticos
cv_vals <- as.numeric(jo@cval[,"5pct"])   # críticos al 5%

L <- min(length(ts_vals), length(cv_vals))
ts_vals <- ts_vals[seq_len(L)]
cv_vals <- cv_vals[seq_len(L)]


hyp_rown <- rownames(jo@cval)
if (is.null(hyp_rown) || length(hyp_rown) < L || all(hyp_rown == "")) {
  # Orden típico del test de traza: r ≤ m-1, r ≤ m-2, ..., r = 0
  hyp_rown <- paste0("r \u2264 ", (L-1):0)  # \u2264 = "≤"
} else {
  hyp_rown <- hyp_rown[seq_len(L)]
}

trace_tab <- data.frame(
  Hipotesis = hyp_rown,
  Trace     = round(ts_vals, 2),
  CV_5pct   = round(cv_vals, 2),
  row.names = NULL
)

kbl(trace_tab, "Tabla D. Test de traza de Johansen (5%)")
Tabla D. Test de traza de Johansen (5%)
Hipotesis Trace CV_5pct
r <= 4 | 9.89 12.25
r <= 3 | 26.27 25.32
r <= 2 | 45.44 42.44
r <= 1 | 78.90 62.99
r = 0 | 115.54 87.31
kbl(as.data.frame(t(sel)), "Tabla C. Selección de rezagos (AIC/HQ/SC/FPE)")
Tabla C. Selección de rezagos (AIC/HQ/SC/FPE)
AIC(n) HQ(n) SC(n) FPE(n)
4 4 1 4
kbl(trace_tab, "Tabla D. Test de traza de Johansen (5%)")
Tabla D. Test de traza de Johansen (5%)
Hipotesis Trace CV_5pct
r <= 4 | 9.89 12.25
r <= 3 | 26.27 25.32
r <= 2 | 45.44 42.44
r <= 1 | 78.90 62.99
r = 0 | 115.54 87.31
cat(sprintf("\n**Decisión de modelo**: con K=%d, r=%d ⇒ %s\n",
            K, r, ifelse(r >= 1, "VECM", "VAR en diferencias")))
## 
## **Decisión de modelo**: con K=2, r=4 ⇒ VECM

Resultados del modelo principal

is_vecm <- r >= 1L

VECM (si r ≥ 1)

vecm_lr  <- cajorls(jo, r = r)     # largo + corto
vecm_var <- vec2var(jo, r = r)     # IRF/FEVD

# Tabla 1: Vectores de cointegración (elasticidades LP)
TABLA1 <- round(vecm_lr$beta, 4)
kbl(TABLA1, "Tabla 1. Cointegración (elasticidades de largo plazo)")
Tabla 1. Cointegración (elasticidades de largo plazo)
ect1 ect2 ect3 ect4
ln_pib.l1 1.0000 0.0000 0.0000 0.0000
ln_ied.l1 0.0000 1.0000 0.0000 0.0000
ln_rem.l1 0.0000 0.0000 1.0000 0.0000
ln_x.l1 0.0000 0.0000 0.0000 1.0000
ln_m.l1 -0.2107 -3.3029 0.6550 -0.7986
trend.l1 -0.0133 0.0399 -0.0409 0.0017
# Tabla 2: Corto plazo (Δln_pib) - robusta (coeficientes)
coefs_mat <- vecm_lr$rlm$coefficients
icol <- which(grepl("pib", colnames(coefs_mat), ignore.case=TRUE))[1]
TABLA2 <- data.frame(Parametro = rownames(coefs_mat),
                     Estimate  = round(coefs_mat[,icol],4),
                     row.names = NULL)
kbl(TABLA2, "Tabla 2. Ecuación de corto plazo (Δln_pib) - VECM")
Tabla 2. Ecuación de corto plazo (Δln_pib) - VECM
Parametro Estimate
ect1 -0.5501
ect2 -0.0103
ect3 -0.0305
ect4 -0.0341
constant 5.2615
ln_pib.dl1 0.3018
ln_ied.dl1 -0.0012
ln_rem.dl1 0.0435
ln_x.dl1 -0.1474
ln_m.dl1 -0.0182
# IRFs (respuesta de ln_pib a shocks externos)
shocks <- intersect(c("ln_ied","ln_rem","ln_x","ln_m","ln_open"), colnames(Y))
if (length(shocks) > 0) {
  par(mfrow=c(2,2), mar=c(3,3,2,1))
  for (imp in head(shocks, 4)) {
    plot(irf(vecm_var, impulse=imp, response="ln_pib", n.ahead=20, ortho=FALSE, boot=FALSE),
         main=paste("IRF:", gsub("ln_","",imp), "→ PIB"))
  }
  par(mfrow=c(1,1))
}

# FEVD de ln_pib (h=4,12,20)
fe <- fevd(vecm_var, n.ahead=20)
if ("ln_pib" %in% names(fe)) {
  H <- c(4,12,20)
  TABLA3 <- sapply(H, function(h) fe$ln_pib[h, ])
  colnames(TABLA3) <- paste0("h=", H)
  kbl(round(TABLA3,3), "Tabla 3. FEVD de ln_pib")
  plot(fe)
}

Causalidad de Granger y diagnósticos

jo <- ca.jo(DY, ecdet="const", type="trace", K=2)
vecm_var <- vec2var(jo, r=1)
# Causalidad de Granger hacia (Δ)ln_pib
suppressWarnings({
  if (!exists("var_d") && !exists("vecm_var")) {
    stop("No hay objeto de modelo disponible (ni var_d ni vecm_var).")
  }
  modelo <- if (exists("var_d")) var_d else vecm_var

  objetivo <- if ("dln_pib" %in% names(modelo$y)) "dln_pib" else "ln_pib"

  shocks_posibles <- c("dln_ied","dln_rem","dln_x","dln_m","ln_ied","ln_rem","ln_x","ln_m")
  shocks <- intersect(shocks_posibles, colnames(modelo$y))
  shocks <- setdiff(shocks, objetivo)  # no testear objetivo causando a sí mismo
  if (length(shocks) == 0) stop("No se encontraron shocks válidos en el sistema.")

  get_p <- function(causa) {
    res <- try(vars::causality(modelo, cause = causa), silent = TRUE)
    if (inherits(res, "try-error") || is.null(res$Granger$p.value)) return(NA_real_)
    as.numeric(res$Granger$p.value)
  }

  pvals <- sapply(shocks, get_p)
  gr_tab <- data.frame(
    Shock   = shocks,
    `p-value` = round(pvals, 4),
    `Signif.` = cut(pvals,
                    breaks = c(-Inf, .01, .05, .10, Inf),
                    labels = c("***","**","*","")),
    check.names = FALSE
  )

  # Tabla en Rmd y consola
  if (requireNamespace("kableExtra", quietly = TRUE)) {
    print(
      knitr::kable(gr_tab,
                   caption = "Tabla 4. Causalidad de Granger hacia el (Δ)PIB",
                   align = "lcc") |>
        kableExtra::kable_styling(full_width = FALSE)
    )
  } else {
    print(gr_tab)
  }

  try(utils::write.csv(gr_tab, "Tabla4_Granger.csv", row.names = FALSE), silent = TRUE)
})
##           Shock p-value Signif.
## dln_ied dln_ied      NA    <NA>
## dln_rem dln_rem      NA    <NA>
## dln_x     dln_x      NA    <NA>
## dln_m     dln_m      NA    <NA>