Resultados Empiricos

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)
}

##No lo utilizamos por eliminar nuestro factor de equilibrio el ECT (término de corrección de error) ### 4.2 VAR en diferencias (si r = 0)

sel_d <- VARselect(DY, lag.max=4, type="const")$selection
p <- if (!is.na(sel_d["SC(n)"])) sel_d["SC(n)"] else if (!is.na(sel_d["HQ(n)"])) sel_d["HQ(n)"] else 2L
p <- max(1L, as.integer(p))
var_d <- VAR(DY, p=p, type="const")

# Tabla 2: ecuación Δln_pib (robusta a 2 o 4 columnas en coeftest)
ct_raw <- coeftest(var_d$varresult$dln_pib)
mat    <- as.matrix(ct_raw)
rn     <- rownames(mat)
kcols  <- ncol(mat)
dfres  <- tryCatch(df.residual(var_d$varresult$dln_pib), error = function(e) NA)

if (kcols == 4) {
  # Formato clásico: Estimate, Std. Error, t/z value, Pr(>|t|)
  TABLA2 <- data.frame(
    Parametro = rn,
    Estimate  = mat[,1],
    `Std.Error` = mat[,2],
    `t value` = mat[,3],
    `Pr(>|t|)` = mat[,4],
    row.names = NULL
  )
} else if (kcols == 2) {
  # Estimate y Std.Error -> calculamos t y p con df residuales
  tval <- mat[,1] / mat[,2]
  # t-Student con dfres, si no, con normal
  pval <- if (is.finite(dfres)) 2*pt(abs(tval), df = dfres, lower.tail = FALSE) else 2*pnorm(abs(tval), lower.tail = FALSE)
  TABLA2 <- data.frame(
    Parametro = rn,
    Estimate  = mat[,1],
    `Std.Error` = mat[,2],
    `t value` = tval,
    `Pr(>|t|)` = pval,
    row.names = NULL
  )
} else {
  
  TABLA2 <- data.frame(Parametro = rn, row.names = NULL)
  nm <- colnames(mat); if (is.null(nm)) nm <- paste0("Col", seq_len(kcols))
  for (j in seq_len(kcols)) TABLA2[[nm[j]]] <- mat[,j]
  # si existen al menos 2 columnas, intenta derivar t y p
  if (!("t value" %in% names(TABLA2)) && kcols >= 2) {
    TABLA2[["t value"]] <- TABLA2[[2]] != 0
    suppressWarnings({
      TABLA2[["t value"]] <- TABLA2[[1]] / TABLA2[[2]]
    })
  }
  if (!("Pr(>|t|)" %in% names(TABLA2)) && "t value" %in% names(TABLA2)) {
    tv <- TABLA2[["t value"]]
    TABLA2[["Pr(>|t|)"]] <- if (is.finite(dfres)) 2*pt(abs(tv), df = dfres, lower.tail = FALSE) else 2*pnorm(abs(tv), lower.tail = FALSE)
  }
}

# Redondeo amable
TABLA2[] <- lapply(TABLA2, function(x) if (is.numeric(x)) round(x, 4) else x)

# Tabla
kbl(TABLA2, "Tabla 2. Corto plazo (Δln_pib) - VAR(Δ)")

Tabla 2. Corto plazo (Δln_pib) - VAR(Δ)

Parametro	Estimate	Std.Error	t.value	Pr…t..
dln_pib.l1	0.1106	0.1229	0.9003	0.3716
dln_ied.l1	-0.0078	0.0073	-1.0751	0.2866
dln_rem.l1	0.0237	0.0444	0.5343	0.5951
dln_x.l1	-0.2222	0.0342	-6.4955	0.0000
dln_m.l1	0.1138	0.0609	1.8676	0.0667
const	0.0156	0.0038	4.0627	0.0001

# IRFs (respuesta de Δln_pib)
shocks <- intersect(c("dln_ied","dln_rem","dln_x","dln_m"), colnames(DY))
if (length(shocks) > 0) {
  par(mfrow=c(2,2), mar=c(3,3,2,1))
  for (imp in head(shocks, 4)) {
    plot(irf(var_d, impulse=imp, response="dln_pib", n.ahead=20, boot=FALSE),
         main=paste("IRF:", gsub("dln_","Δ",imp), "→ ΔPIB"))
  }
  par(mfrow=c(1,1))
}

# FEVD de Δln_pib (h=4,12,20)
fe <- fevd(var_d, n.ahead=20)
if ("dln_pib" %in% names(fe)) {
  H <- c(4,12,20)
  TABLA3 <- sapply(H, function(h) fe$dln_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)
})

## <table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
## <caption>Tabla 4. Causalidad de Granger hacia el (Δ)PIB</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;">  </th>
##    <th style="text-align:left;"> Shock </th>
##    <th style="text-align:center;"> p-value </th>
##    <th style="text-align:center;"> Signif. </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> dln_ied </td>
##    <td style="text-align:left;"> dln_ied </td>
##    <td style="text-align:center;"> 0.7787 </td>
##    <td style="text-align:center;">  </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> dln_rem </td>
##    <td style="text-align:left;"> dln_rem </td>
##    <td style="text-align:center;"> 0.4574 </td>
##    <td style="text-align:center;">  </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> dln_x </td>
##    <td style="text-align:left;"> dln_x </td>
##    <td style="text-align:center;"> 0.0000 </td>
##    <td style="text-align:center;"> *** </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> dln_m </td>
##    <td style="text-align:left;"> dln_m </td>
##    <td style="text-align:center;"> 0.0461 </td>
##    <td style="text-align:center;"> ** </td>
##   </tr>
## </tbody>
## </table>