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)
| 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%)
| 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)
| 4 |
4 |
1 |
4 |
kbl(trace_tab, "Tabla D. Test de traza de Johansen (5%)")
Tabla D. Test de traza de Johansen (5%)
| 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)
| 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
| 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>