---
title: "Capítulo 5: Multiequation Time-Series Models"
subtitle: "Enders (2015) Applied Econometric Time Series — Sección completa con ejemplos en R"
author: "ECN 4082 Tópicos Avanzados de Economía · Prof. Carlos Uribe Terán · USFQ"
date: today
format:
html:
toc: true
toc-depth: 3
toc-location: left
number-sections: true
theme: cosmo
code-fold: show
code-tools: true
highlight-style: github
embed-resources: true
fig-width: 9
fig-height: 5
df-print: kable
execute:
warning: false
message: false
echo: true
---
```{r setup}
#| include: false
# Instalar paquetes si no están disponibles
pkgs <- c("tseries", "forecast", "dynlm", "vars", "lmtest", "zoo")
for (p in pkgs) {
if (!requireNamespace(p, quietly = TRUE))
install.packages(p, repos = "https://cloud.r-project.org")
}
library(tseries)
library(forecast)
library(dynlm)
library(vars)
library(lmtest)
library(zoo)
# Paleta de colores consistente
col_y <- "#2166ac" # azul para y_t
col_z <- "#d6604d" # rojo para z_t
col_irf <- "#1b7837" # verde para IRF
col_band <- "#a6d96a" # verde claro para bandas CI
# Datos del libro (ajusta la ruta si es necesario)
DATA_PATH <- "." # carpeta donde están italy.csv y terrorism.csv
```
# Introducción: el hilo del capítulo
Enders abre el Capítulo 5 con seis objetivos de aprendizaje que definen exactamente el arco:
1. **Introducir** el análisis de intervención y la función de transferencia
2. **Mostrar** que la función de transferencia es efectiva cuando **no hay feedback** de $y$ a $z$
3. **Estimar** un ADL con los datos de terrorismo en Italia
4. **Explicar** por qué la limitación principal es el **feedback** entre variables
5. **Introducir** el VAR como la respuesta a esa limitación
6. **Estimar** un VAR
El hilo es una escalera de supuestos:
$$\underbrace{z_t \text{ dummy determinística}}_{\text{Intervención}} \longrightarrow \underbrace{z_t \text{ estocástica, exógena}}_{\text{ARDL / Transfer Function}} \longrightarrow \underbrace{y_t, z_t \text{ endógenas}}_{\text{VAR}}$$
---
# Sección 1: Intervention Analysis {#sec-intervencion}
## ¿Por qué existe este modelo?
Supongamos que tienes el PIB del Ecuador trimestre a trimestre y quieres medir el efecto de la dolarización (enero 2000). No puedes hacer un DiD clásico porque:
- **No tienes grupo de control** — solo observas una serie de tiempo
- **La serie tiene correlación serial** — los residuos no son independientes
El análisis de intervención resuelve ambos problemas: usa la historia de la serie como contrafactual implícito y modela explícitamente la autocorrelación serial.
> *"The crucial assumption in intervention analysis is that the intervention function has only deterministic components."* — Enders p. 261
## La ecuación base (Enders ec. 5.1)
$$\boxed{y_t = a_0 + a_1 y_{t-1} + \omega z_t + \varepsilon_t, \quad |\,a_1\,| < 1}$$
| Pieza | Qué hace |
|---|---|
| $a_0$ | Constante. Media pre-intervención: $\mu_{\text{pre}} = a_0/(1-a_1)$ |
| $a_1 y_{t-1}$ | Inercia (memoria) de la serie — la autocorrelación serial |
| $\omega z_t$ | Efecto de la intervención — activa o desactiva un cambio en el promedio |
| $\varepsilon_t$ | Ruido blanco — lo que el modelo no explica |
Cuando $z_t = 1$: el intercepto efectivo cambia de $a_0$ a $a_0 + \omega$, y la media LP pasa de $\frac{a_0}{1-a_1}$ a $\frac{a_0+\omega}{1-a_1}$.
## Los tres tipos de $z_t$
### Tipo 1 — Level shift (salto permanente)
$$z_t = D_t^L = \begin{cases} 0 & t < \tau \\ 1 & t \geq \tau \end{cases}$$
Ejemplo de Enders: instalación de detectores de metales en aeropuertos. Antes de $\tau$: no existen. Desde $\tau$: siempre están. El efecto es **permanente**.
### Tipo 2 — Pulso (shock transitorio)
$$z_t = D_t^P = \begin{cases} 1 & t = \tau \\ 0 & t \neq \tau \end{cases}$$
Ejemplo: el bombardeo de Libia en abril 1986. Un solo periodo de activación. El efecto decae según la inercia $a_1$.
### Tipo 3 — Impulso prolongado
$$z_\tau = 1,\; z_{\tau+1} = 0.75,\; z_{\tau+2} = 0.5,\; z_{\tau+3} = 0.25,\; z_t = 0 \text{ para } t > \tau+3$$
El efecto tiene un máximo y luego se disipa gradualmente. No es una dummy 0/1.
```{r tipos-intervencion}
#| fig-cap: "Los tres tipos de variable de intervención z_t (τ = 25)"
tau <- 25; T <- 50
z_level <- c(rep(0, tau-1), rep(1, T-tau+1))
z_pulse <- c(rep(0, tau-1), 1, rep(0, T-tau))
z_ramp <- c(rep(0, tau-1), 1, 0.75, 0.5, 0.25, rep(0, T-tau-3))
par(mfrow=c(1,3), mar=c(4,4,3,1))
plot(1:T, z_level, type="s", col=col_z, lwd=2,
main="Level Shift\n(cambio permanente)",
xlab="Periodo t", ylab=expression(z[t]), ylim=c(-0.1,1.1))
abline(v=tau, lty=2, col="gray50")
plot(1:T, z_pulse, type="h", col=col_z, lwd=3,
main="Pulso\n(shock transitorio)",
xlab="Periodo t", ylab=expression(z[t]), ylim=c(-0.1,1.1))
abline(v=tau, lty=2, col="gray50")
plot(1:T, z_ramp, type="h", col=col_z, lwd=3,
main="Impulso Prolongado\n(disipación gradual)",
xlab="Periodo t", ylab=expression(z[t]), ylim=c(-0.1,1.1))
abline(v=tau, lty=2, col="gray50")
par(mfrow=c(1,1))
```
## Derivación de los efectos — el pizarrón de Carlos {#sec-efectos}
Con $y_t = a_0 + a_1 y_{t-1} + \omega z_t + \varepsilon_t$ y un **level shift** que entra en $\tau$:
**Iterando hacia adelante (sustitución sucesiva):**
$$y_{\tau+h} = a_0(1 + a_1 + \cdots + a_1^h) + a_1^{h+1} y_{\tau-1} + \omega(1 + a_1 + \cdots + a_1^h) + \text{errores}$$
El efecto acumulado del cambio de $z_t$ en el horizonte $h$ es:
$$\underbrace{\omega \sum_{k=0}^{h} a_1^k}_{\text{efecto en } \tau+h} = \omega \cdot \frac{1 - a_1^{h+1}}{1-a_1}$$
Cuando $h \to \infty$ y $|a_1| < 1$:
$$\boxed{\text{Multiplicador de Largo Plazo} = \frac{\omega}{1-a_1}}$$
Para el **pulso** ($z_\tau = 1$, luego $z_t = 0$):
$$\text{Efecto en } \tau+h = \omega \cdot a_1^h \xrightarrow{h\to\infty} 0$$
El shock se disipa. No hay efecto de largo plazo.
```{r irf-tipos}
#| fig-cap: "IRF de los tres tipos de intervención (a1=0.7, omega=2)"
a1 <- 0.7; omega <- 2; H <- 20
# Level shift: efecto en tau+h = omega * sum(a1^k, k=0..h)
irf_level <- omega * sapply(0:H, function(h) (1 - a1^(h+1))/(1-a1))
# Pulso: efecto en tau+h = omega * a1^h
irf_pulse <- omega * a1^(0:H)
# Impulso prolongado (numérico)
irf_ramp <- numeric(H+1)
z_ramp_s <- c(1, 0.75, 0.5, 0.25, rep(0, H-3))
for(h in 0:H) {
irf_ramp[h+1] <- sum(omega * z_ramp_s[1:min(h+1,4)] * a1^(h:(h-min(h,3))))
}
par(mar=c(4,4,3,1))
plot(0:H, irf_level, type="l", col=col_y, lwd=2.5,
ylim=c(0, max(irf_level)*1.1),
xlab="Horizonte h (periodos desde intervención)",
ylab="Efecto en y_{τ+h}",
main=paste0("IRF por tipo de intervención\n(a₁ = ", a1, ", ω = ", omega, ")"))
lines(0:H, irf_pulse, col=col_z, lwd=2.5, lty=2)
lines(0:H, irf_ramp, col=col_irf, lwd=2.5, lty=3)
abline(h=omega/(1-a1), lty=4, col=col_y, lwd=1)
legend("right",
legend=c(paste("Level shift (LP =", round(omega/(1-a1),2), ")"),
"Pulso (LP = 0)",
"Impulso prolongado"),
col=c(col_y, col_z, col_irf), lwd=2.5, lty=1:3, bty="n")
```
```{r tabla-multiplicadores}
#| label: tbl-mult
#| tbl-cap: "Multiplicadores por horizonte (a₁ = 0.7, ω = 2)"
h_vals <- c(0,1,2,3,5,10,20)
data.frame(
`Horizonte h` = h_vals,
`Level Shift` = round(omega * sapply(h_vals, function(h) (1-a1^(h+1))/(1-a1)), 4),
`Pulso` = round(omega * a1^h_vals, 4),
`LP (limit)` = round(rep(omega/(1-a1), length(h_vals)), 4)
)
```
## Ejemplo de Enders: el bombardeo de Libia (pp. 265–267)
**Contexto:** Datos mensuales de ataques terroristas contra ciudadanos de EE.UU. y el Reino Unido, enero 1968 – marzo 1986. Variable: número de ataques por mes.
**Evento:** El 15 de abril de 1986 EE.UU. bombardeó Libia en represalia por el atentado de la discoteca de Berlín.
**Pregunta:** ¿Redujo el bombardeo los ataques terroristas?
Enders compara dos especificaciones:
```{r libia-resultados}
#| label: tbl-libia
#| tbl-cap: "Comparación de modelos — Bombardeo de Libia (Enders Tabla de p. 266)"
data.frame(
Especificación = c(
"Level shift (z=1 desde abril 1986)",
"Pulso (z=1 solo en abril 1986)"
),
Ecuación = c(
"5.58 + 0.336y_{t-1} + 0.123y_{t-5} + 2.65z_t",
"3.79 + 0.327y_{t-1} + 0.157y_{t-5} + 38.9z_t"
),
`t de z_t` = c("0.84 (NS)", "6.09 ***"),
AIC = c(1656.3, 1608.7),
Decisión = c("Rechazar", "**Preferida**")
)
```
**Conclusión de Enders:** El bombardeo **aumentó** los ataques en ~39 en el mes del evento. Con $a_1 = 0.327$, el efecto en el siguiente mes fue $0.327 \times 38.9 \approx 12.7$ ataques. El efecto de largo plazo es **cero** (serie estacionaria, pulso se disipa).
```{r libia-irf}
#| fig-cap: "IRF del bombardeo de Libia — respuesta del número de ataques"
a1_lib <- 0.327; omega_lib <- 38.9; H_lib <- 15
irf_lib <- omega_lib * a1_lib^(0:H_lib)
par(mar=c(4,4,3,1))
plot(0:H_lib, irf_lib, type="h", lwd=4, col=col_z,
xlab="Meses después del bombardeo",
ylab="Ataques adicionales",
main="IRF: Efecto del bombardeo de Libia\nsobre ataques terroristas")
abline(h=0, lwd=1)
points(0:H_lib, irf_lib, pch=19, col=col_z, cex=0.8)
text(1, 35, paste("Impacto:", omega_lib), col=col_z, adj=0)
text(1, 30, paste("LP = 0 (a₁ =", a1_lib, "< 1)"), col="gray40", adj=0)
```
---
# Sección 2: ADLs y Transfer Functions {#sec-ardl}
## El salto conceptual
La Sección 1 tenía $z_t$ como dummy **determinística** (tú la construyes). La Sección 2 permite que $z_t$ sea cualquier **serie de tiempo observada, estocástica y exógena**.
> *"A natural extension of the intervention model is to allow the {zt} sequence to be something other than a deterministic dummy variable."* — Enders p. 267
La ecuación sigue siendo la misma estructura general:
$$\boxed{y_t = a_0 + A(L)y_{t-1} + C(L)z_t + B(L)\varepsilon_t} \tag{5.3}$$
| Polinomio | Forma | Interpretación |
|---|---|---|
| $A(L) = a_1 + a_2 L + \cdots + a_p L^{p-1}$ | $A(L)y_{t-1} = \sum_{i=1}^p a_i y_{t-i}$ | Parte AR de $y_t$ (inercia) |
| $C(L) = c_0 + c_1 L + \cdots + c_n L^n$ | $C(L)z_t = \sum_{j=0}^n c_j z_{t-j}$ | **Función de transferencia** |
| $B(L) = 1 + b_1 L + \cdots + b_q L^q$ | $B(L)\varepsilon_t$ | Error ARMA |
Los coeficientes $\{c_j\}$ son los **transfer function weights**: dicen cómo y cuándo un cambio en $z_t$ se transmite a $y_t$.
Si $c_0 = 0$: el contemporáneo de $z$ no afecta a $y$ → $z_t$ es un **leading indicator**.
## Condición de exogeneidad
> *"Transfer function analysis assumes that {zt} is an exogeneous process that evolves independently of the {yt} sequence... $E z_t \varepsilon_{t-s} = 0$ for all values of s and t."* — Enders p. 268
Si hay **feedback** de $y$ a $z$, los coeficientes de $C(L)$ están sesgados. La solución es el VAR.
## La CCVF teórica: el instrumento de identificación
### Caso benchmark: $z_t$ es ruido blanco con delay $d$ (Enders ec. 5.5)
$$y_t = a_1 y_{t-1} + c_d z_{t-d} + \varepsilon_t, \quad \{z_t\} \perp \{\varepsilon_t\} \text{ (ruido blanco independientes)}$$
**Derivación algebraica completa:**
Usando el operador de rezago $(1 - a_1 L)y_t = c_d z_{t-d} + \varepsilon_t$:
$$y_t = \frac{c_d z_{t-d}}{1-a_1 L} + \frac{\varepsilon_t}{1-a_1 L} = c_d \sum_{k=0}^{\infty} a_1^k z_{t-d-k} + \frac{\varepsilon_t}{1-a_1 L}$$
Calculando $\gamma_{yz}(i) = \text{Cov}(y_t, z_{t-i}) = E[y_t \cdot z_{t-i}]$ (con $E[z_t]=0$):
$$E\left[\left(\sum_{k=0}^{\infty} c_d a_1^k z_{t-d-k}\right) z_{t-i}\right] = c_d a_1^{i-d} \cdot E[z_{t-i}^2] = c_d a_1^{i-d} \sigma_z^2$$
porque $E[z_{t-d-k} z_{t-i}] = \sigma_z^2$ solo cuando $i = d+k$, y cero en otro caso. Entonces:
$$\boxed{\gamma_{yz}(i) = \text{Cov}(y_t, z_{t-i}) = \begin{cases} 0 & i < d \\ c_d \cdot a_1^{i-d} \cdot \sigma_z^2 & i \geq d \end{cases}} \tag{5.6}$$
**El patrón que usamos para identificar el modelo:**
- **Cero hasta el lag $d$** → te da el **delay** (cuántos periodos tarda $z$ en afectar a $y$)
- **Espiga en $d$** → proporcional a $c_d$ (el primer peso no-nulo)
- **Decaimiento geométrico a tasa $a_1$** para $i > d$ → la dinámica AR de $y_t$
```{r ccvf-patrones}
#| fig-cap: "Patrones teóricos de la CCVF — Figura 5.2 de Enders"
par(mfrow=c(2,2), mar=c(4,4,3,1))
# Panel (a): d=0, a1=0.7 > 0 (decaimiento monótono)
lags <- 0:12
ccvf_a <- 1 * 0.7^lags
plot(lags, ccvf_a, type="h", lwd=2.5, col=col_z,
xlab="Lag i", ylab=expression(gamma[yz](i)),
main="(a) d=0, a₁=0.7 (monótono +)")
abline(h=0)
# Panel (b): d=0, a1=-0.5 < 0 (oscilatorio)
ccvf_b <- 1 * (-0.5)^lags
plot(lags, ccvf_b, type="h", lwd=2.5, col=col_z,
xlab="Lag i", ylab=expression(gamma[yz](i)),
main="(b) d=0, a₁=-0.5 (oscilatorio)")
abline(h=0)
# Panel (c): d=3, a1=0.7 (delay de 3 periodos)
ccvf_c <- c(0,0,0, 1*0.7^lags[1:10])
plot(0:12, ccvf_c, type="h", lwd=2.5, col=col_y,
xlab="Lag i", ylab=expression(gamma[yz](i)),
main="(c) d=3, a₁=0.7 (delay de 3)")
abline(h=0)
text(3.2, 0.85, "← espiga en d=3", col=col_y, cex=0.8)
# Panel (d): dos pesos no-nulos c_d y c_{d+1}
ccvf_d <- c(0,0, 0.8, (0.8*0.6+0.4)*0.6^(0:9))
plot(0:12, ccvf_d[1:13], type="h", lwd=2.5, col=col_irf,
xlab="Lag i", ylab=expression(gamma[yz](i)),
main="(d) d=2, dos pesos no-nulos")
abline(h=0)
par(mfrow=c(1,1))
```
### El caso con dos pesos no-nulos: $c_d$ y $c_{d+1}$
Si el modelo verdadero tiene efectos en lags $d$ y $d+1$:
$$y_t = a_1 y_{t-1} + c_d z_{t-d} + c_{d+1} z_{t-d-1} + \varepsilon_t$$
La CCVF tiene el patrón:
$$\gamma_{yz}(i) = \begin{cases} 0 & i < d \\ c_d \sigma_z^2 & i = d \\ (c_d a_1 + c_{d+1})\sigma_z^2 & i = d+1 \\ a_1^{i-d-1}(c_d a_1 + c_{d+1})\sigma_z^2 & i > d+1 \end{cases}$$
Hay **dos espigas** en $d$ y $d+1$ — luego decaimiento geométrico.
## El problema: $z_t$ no es ruido blanco
En la práctica, $z_t$ tiene su propia autocorrelación (ej: el precio del petróleo hoy está correlacionado con el de ayer). Enders dice:
> *"The econometrician will rarely be so fortunate as to work with a {zt} series that is white noise."*
Si $z_t = d_1 z_{t-1} + \varepsilon_{zt}$ (AR(1)) y el verdadero modelo es solo $y_t = c_0 z_t + \varepsilon_t$:
$$\text{Cov}(y_t, z_{t-i}) = c_0 \cdot \text{Cov}(z_t, z_{t-i}) = c_0 \cdot d_1^i \cdot \text{Var}(z) \neq 0 \quad \forall\, i \geq 0$$
La CCF cruda muestra picos en **todos los rezagos** aunque el efecto real sea solo contemporáneo. La autocorrelación de $z_t$ contamina la CCF.
```{r problema-contaminacion}
#| fig-cap: "CCF verdadera vs. CCF cruda contaminada"
set.seed(123)
n <- 200
d1 <- 0.75 # persistencia de z_t
# Simular: verdadero modelo y_t = 0.8*z_t + eps_t
eps_z <- rnorm(n)
z <- numeric(n); z[1] <- eps_z[1]
for(t in 2:n) z[t] <- d1*z[t-1] + eps_z[t]
eps_y <- rnorm(n)
y <- 0.8*z + eps_y
# CCF cruda
par(mfrow=c(1,2), mar=c(4,4,3,1))
ccf_cruda <- ccf(y, z, lag.max=12, plot=FALSE)
barplot(ccf_cruda$acf[ccf_cruda$lag >= 0],
names.arg = 0:12, col=col_z, border=NA,
ylim=c(-0.2, 0.8),
xlab="Lag", ylab="Correlación cruzada",
main="CCF CRUDA entre y_t y z_{t-i}\n(engañosa)")
abline(h=c(2,-2)/sqrt(n), lty=2, col="gray50")
text(5, 0.7, "Picos en todos los lags\npor autocorr. de z", col=col_z, cex=0.85)
# CCF con pre-whitening
ar_z <- arima(z, order=c(1,0,0), include.mean=FALSE)
eps_zt <- residuals(ar_z)
yf <- y - d1*c(0, y[-n]) # filtrar y con D(L)
ccf_filt <- ccf(yf[-1], eps_zt[-1], lag.max=12, plot=FALSE)
barplot(ccf_filt$acf[ccf_filt$lag >= 0],
names.arg = 0:12, col=col_y, border=NA,
ylim=c(-0.2, 0.8),
xlab="Lag", ylab="Correlación cruzada",
main="CCF FILTRADA entre y^f_t y ε̂_{zt}\n(correcta)")
abline(h=c(2,-2)/sqrt(n), lty=2, col="gray50")
text(3, 0.7, "Solo lag 0 es\nsignificativo ✓", col=col_y, cex=0.85)
par(mfrow=c(1,1))
```
## El procedimiento de 3 pasos de Enders (pp. 275–276)
**PASO 1 — Modelar $z_t$ como AR y extraer innovaciones**
Usar Box-Jenkins para modelar $z_t$:
$$z_t = d_1 z_{t-1} + d_2 z_{t-2} + \cdots + \varepsilon_{zt}, \quad \text{i.e.,}\quad D(L)z_t = \varepsilon_{zt}$$
Guardar los residuos $\hat{\varepsilon}_{zt}$. Por construcción son ruido blanco.
**PASO 2 — Filtrar $y_t$ con el mismo $D(L)$**
Construir $y_t^f = D(L)y_t$. Verificar algebraicamente que esto funciona:
Multiplicando $y_t = a_1 y_{t-1} + c_1 z_t + \varepsilon_t$ por $D(L) = (1-d_1 L)$:
$$D(L)y_t = a_1 D(L)y_{t-1} + c_1 \underbrace{D(L)z_t}_{=\,\varepsilon_{zt}} + D(L)\varepsilon_t$$
$$y_t^f = a_1 y_{t-1}^f + c_1 \varepsilon_{zt} + \varepsilon_t - d_1\varepsilon_{t-1}$$
Ahora $\text{Cov}(y_t^f, \varepsilon_{z,t-i})$ revela limpiamente el patrón de $C(L)$.
Calcular la CCF entre $y_t^f$ y $\hat{\varepsilon}_{z,t-i}$. Los picos en el lag $i$ → $c_i \neq 0$.
**Significancia** (Enders p. 276): bajo $H_0$ de CCF = 0:
$$\text{Var}[r_{yz}(i)] \approx (T-i)^{-1} \implies \text{banda de confianza: } \pm\frac{2}{\sqrt{T}}$$
Test conjunto para los primeros $k$ lags:
$$Q = T(T+2)\sum_{i=0}^{k}\frac{r_{yz}^2(i)}{T-k} \sim \chi^2(k-p_1-p_2)$$
**PASO 3 — Identificar $A(L)$**
Regresar $y_t$ (no $y_t^f$) en los lags de $z_t$ identificados: $y_t = C(L)z_t + e_t$.
El ACF de $e_t$ sugiere la forma de $A(L)$.
---
# Sección 3: ADL de Terrorismo en Italia {#sec-italia}
Enders (pp. 277–281) aplica el procedimiento completo con datos trimestrales 1971Q1–1988Q4.
**Variables:**
- $y_t$ = `Slitaly`: logaritmo del share de Italia en turismo mundial
- $z_t$ = `Attkit`: número de ataques terroristas contra turistas en Italia
```{r cargar-italia}
# Cargar datos reales del libro
# Descarga: https://www.wiley.com/en-us/Applied+Econometric+Time+Series%2C+4th+Edition-p-9781118808566
italy_file <- file.path(DATA_PATH, "italy.csv")
if (!file.exists(italy_file)) italy_file <- "/tmp/italy.csv"
italy <- read.csv(italy_file)
italy$DATE <- as.Date(italy$DATE)
slitaly <- ts(italy$Slitaly, start=c(1971,1), frequency=4)
attkit <- ts(italy$Attkit, start=c(1971,1), frequency=4)
cat("Datos cargados:\n")
cat(sprintf(" Observaciones: %d (1971Q1 – 1988Q4)\n", length(slitaly)))
cat(sprintf(" Slitaly: media=%.4f, sd=%.4f\n", mean(slitaly), sd(slitaly)))
cat(sprintf(" Attkit: media=%.4f, sd=%.4f\n", mean(attkit), sd(attkit)))
```
```{r plot-italia}
#| fig-cap: "Series de tiempo — Italia: turismo y terrorismo"
par(mfrow=c(2,1), mar=c(3,4,2,1))
plot(slitaly, col=col_y, lwd=1.5, ylab="Log share turismo",
main="Slitaly: Share de Italia en turismo mundial")
abline(h=mean(slitaly), lty=2, col="gray60")
plot(attkit, col=col_z, lwd=1.5, type="h", ylab="Nº ataques",
main="Attkit: Ataques terroristas contra turistas en Italia")
abline(h=0)
par(mfrow=c(1,1))
```
## Paso 0: Tests ADF de ambas series
```{r adf-italia}
# ADF con arima manual (reproduce tseries::adf.test)
adf_manual <- function(y, k=4, include.const=TRUE, include.trend=FALSE) {
n <- length(y)
dy <- diff(y)
yl <- y[-length(y)] # y_{t-1}
m <- k
# Construir data.frame
df_list <- list(dy=dy[(m+1):length(dy)], yl=yl[(m+1):length(yl)])
for(j in 1:m) df_list[[paste0("dyl",j)]] <- dy[(m+1-j):(length(dy)-j)]
if(include.trend) df_list$trend <- (m+2):n
formula_str <- "dy ~ yl"
if(m > 0) formula_str <- paste(formula_str, paste(paste0("+dyl",1:m), collapse=" "))
if(!include.const) formula_str <- paste(formula_str, "- 1")
if(include.trend) formula_str <- paste(formula_str, "+ trend")
fit <- lm(as.formula(formula_str), data=as.data.frame(df_list))
tau <- summary(fit)$coefficients["yl","t value"]
list(tau=round(tau,4), p_approx=if(tau < -3.5) "<0.01" else if(tau < -2.89) "<0.05" else ">0.05")
}
adf_sl <- adf_manual(slitaly, k=4)
adf_att <- adf_manual(attkit, k=4)
cat("Tests ADF (con constante, 4 lags):\n\n")
cat(sprintf(" Slitaly: τ = %6.3f p %s\n", adf_sl$tau, adf_sl$p_approx))
cat(sprintf(" Attkit: τ = %6.3f p %s\n", adf_att$tau, adf_att$p_approx))
cat("\nAmbas series son estacionarias (I(0)) ✓\n")
cat("→ El ARDL puede estimarse directamente en niveles\n")
```
## Paso 1: Modelar attkit como AR (pre-whitening de $z_t$)
```{r ar-attkit}
# ACF y PACF de attkit para identificar el orden AR
par(mfrow=c(1,2), mar=c(4,4,3,1))
acf(attkit, lag.max=12, main="ACF de Attkit")
pacf(attkit, lag.max=12, main="PACF de Attkit")
par(mfrow=c(1,1))
```
```{r ar-attkit-estim}
# Comparar AR(1) a AR(4) por AIC
cat("Selección de orden AR para attkit (z_t):\n\n")
cat(sprintf("%-6s %-10s %-10s\n", "Orden", "AIC", "BIC"))
aics <- numeric(5)
for(p in 1:5) {
m <- arima(attkit, order=c(p,0,0), include.mean=FALSE)
aics[p] <- AIC(m)
cat(sprintf("AR(%-2d) %-10.2f %-10.2f\n", p, AIC(m), BIC(m)))
}
best_p <- which.min(aics)
cat(sprintf("\n→ Seleccionado: AR(%d)\n", best_p))
# Estimar el AR óptimo
ar_z <- arima(attkit, order=c(best_p,0,0), include.mean=FALSE)
eps_zt <- residuals(ar_z)
cat("\nCoeficientes AR estimados para attkit:\n")
print(round(coef(ar_z), 4))
lb_z <- Box.test(eps_zt, lag=12, type="Ljung-Box")
cat(sprintf("\nLjung-Box residuos AR(%d): Q(12)=%.3f, p=%.4f %s\n",
best_p, lb_z$statistic, lb_z$p.value,
if(lb_z$p.value>0.05) "→ Ruido blanco ✓" else "→ Revisar orden"))
```
## Paso 2: CCF filtrada para identificar $C(L)$
```{r ccf-filtrada-italia}
# Filtrar slitaly con D(L) del AR estimado
phi <- coef(ar_z)
n <- length(slitaly)
p <- best_p
yf <- numeric(n)
for(t in (p+1):n) {
yf[t] <- slitaly[t] - sum(phi * rev(slitaly[(t-p):(t-1)]))
}
yf_trim <- yf[(p+1):n]
eps_zt_trim <- eps_zt[1:(n-p)]
# CCF
ccf_res <- ccf(yf_trim, eps_zt_trim, lag.max=12, plot=FALSE)
band <- 2/sqrt(length(yf_trim))
par(mar=c(4,4,3,1))
lags_pos <- which(ccf_res$lag >= 0)
barplot(ccf_res$acf[lags_pos],
names.arg = ccf_res$lag[lags_pos],
col = ifelse(abs(ccf_res$acf[lags_pos]) > band, col_y, "gray80"),
border = NA,
ylim = c(-0.5, 0.5),
xlab = "Lag i",
ylab = expression(r[yz](i)),
main = "CCF filtrada: y^f vs ε̂_z\n(barras azules = significativas)")
abline(h=c(band,-band), lty=2, col="gray40")
abline(h=0)
text(0, 0.45, sprintf("Banda ±%.3f (±2/√T)", band), col="gray40", cex=0.8)
cat("\nCCF filtrada (lags 0 a 6):\n")
cat(sprintf("%-6s %-12s %-14s\n", "Lag i", "CCF", "Significativo?"))
for(i in lags_pos[ccf_res$lag[lags_pos] <= 6]) {
sig <- abs(ccf_res$acf[i]) > band
cat(sprintf("%-6d %-12.4f %s\n",
ccf_res$lag[i], ccf_res$acf[i],
if(sig) "✓ SÍ" else "—"))
}
cat(sprintf("\n→ Enders identifica: delay d=1 (z_{t-1} afecta a y_t)\n"))
cat(sprintf("→ Forma de C(L): c₁ ≠ 0 (efecto con 1 trimestre de rezago)\n\n"))
```
## Paso 3: Identificar $A(L)$ con ACF de residuos
```{r identificar-AL}
# Regresar slitaly en attkit_{t-1}
reg_cl <- lm(slitaly[-1] ~ attkit[-length(attkit)])
e_cl <- residuals(reg_cl)
par(mfrow=c(1,2), mar=c(4,4,3,1))
acf(e_cl, lag.max=10, main="ACF residuos de y_t ~ z_{t-1}")
pacf(e_cl, lag.max=10, main="PACF residuos")
par(mfrow=c(1,1))
band_e <- 2/sqrt(length(e_cl))
acf_e <- acf(e_cl, lag.max=6, plot=FALSE)
cat("ACF de residuos intermedios:\n")
for(k in 1:6) {
sig <- abs(acf_e$acf[k+1]) > band_e
cat(sprintf(" ACF(%d) = %7.4f %s\n", k, acf_e$acf[k+1],
if(sig) "✓ → necesita AR" else ""))
}
cat("\n→ Enders identifica A(L) = AR(2): a1·y_{t-1} + a2·y_{t-2}\n")
```
## Paso 4: Estimar el modelo ADL completo
```{r ardl-completo}
# ADL(2,1): y_t = a1*y_{t-1} + a2*y_{t-2} + c1*z_{t-1} + eps_t
# Usando dynlm o equivalente en base R con lag
n <- length(slitaly)
y_t <- slitaly[3:n]
y_t1 <- slitaly[2:(n-1)]
y_t2 <- slitaly[1:(n-2)]
z_t1 <- attkit[2:(n-1)]
ardl <- lm(y_t ~ y_t1 + y_t2 + z_t1)
cat("═══ Modelo ADL(2,1) — y_t = a1*y_{t-1} + a2*y_{t-2} + c1*z_{t-1} ═══\n\n")
print(summary(ardl))
# Comparar varias especificaciones
cat("\nComparación de especificaciones (Enders Tabla 5.1):\n\n")
cat(sprintf("%-30s %-8s %-8s %-12s\n", "Modelo", "AIC", "BIC", "Box-Test p"))
specs <- list(
"AR(2) + z_{t-1}" = lm(y_t ~ y_t1 + y_t2 + z_t1),
"AR(1) + z_{t-1}" = lm(y_t ~ y_t1 + z_t1),
"AR(2) + z_{t-1,2}" = lm(y_t ~ y_t1 + y_t2 + z_t1 + attkit[1:(n-2)]),
"solo z_{t-1,2}" = lm(y_t ~ z_t1 + attkit[1:(n-2)])
)
for(nm in names(specs)) {
m <- specs[[nm]]
lb <- Box.test(residuals(m), lag=8, type="Ljung-Box")
cat(sprintf("%-30s %-8.1f %-8.1f %-12.4f\n",
nm, AIC(m), BIC(m), lb$p.value))
}
```
## Paso 5: Multiplicadores e IRF
```{r multiplicadores-italia}
# Extraer coeficientes del modelo final
co <- coef(ardl)
a1_hat <- co["y_t1"]
a2_hat <- co["y_t2"]
c1_hat <- co["z_t1"]
cat("Coeficientes estimados:\n")
cat(sprintf(" â₁ = %8.5f\n", a1_hat))
cat(sprintf(" â₂ = %8.5f\n", a2_hat))
cat(sprintf(" ĉ₁ = %8.5f ← peso de la función de transferencia\n\n", c1_hat))
# Multiplicador de largo plazo
denom <- 1 - a1_hat - a2_hat
LP <- c1_hat / denom
cat(sprintf("Multiplicador de LP = ĉ₁ / (1-â₁-â₂) = %.5f / %.5f = %.5f\n\n",
c1_hat, denom, LP))
# IRF numérica: impulso de 1 ataque en attkit
H <- 20
y_irf <- numeric(H+3); z_irf <- c(0, 1, rep(0, H+2))
for(t in 3:(H+3)) {
y_irf[t] <- a1_hat*y_irf[t-1] + a2_hat*y_irf[t-2] + c1_hat*z_irf[t-1]
}
irf_it <- y_irf[3:(H+3)]
cat("IRF: efecto de 1 ataque terrorista adicional sobre el share de turismo:\n\n")
cat(sprintf("%-6s %-14s %-14s\n", "h", "Respuesta y", "Acumulado"))
acum <- 0
for(h in 0:10) {
acum <- acum + irf_it[h+1]
cat(sprintf("%-6d %-14.6f %-14.6f\n", h, irf_it[h+1], acum))
}
cat(sprintf("\nEfecto total acumulado (Σ IRF): %.6f\n", sum(irf_it)))
cat("Enders estima: pérdidas > 600 millones SDR\n")
```
```{r irf-italia-plot}
#| fig-cap: "IRF: respuesta del turismo italiano a un ataque terrorista"
par(mar=c(4,5,3,1))
plot(0:H, irf_it, type="h", lwd=3, col=col_z,
xlab="Horizonte h (trimestres)",
ylab="Cambio en log-share de turismo",
main="IRF — Efecto de 1 ataque terrorista\nsobre el turismo en Italia")
abline(h=0, lwd=1)
abline(h=LP, lty=2, col="gray50", lwd=1)
points(0:H, irf_it, pch=19, col=col_z, cex=0.7)
text(12, LP*1.3, sprintf("LP = %.5f", LP), col="gray40", cex=0.85)
```
## Paso 6: Diagnóstico final
```{r diagnostico-italia}
resid_ardl <- residuals(ardl)
par(mfrow=c(2,2), mar=c(4,4,3,1))
plot(resid_ardl, type="l", col="gray40",
main="Residuos del ADL", ylab="Residuo", xlab="t")
abline(h=0, lty=2)
acf(resid_ardl, lag.max=12, main="ACF de residuos")
pacf(resid_ardl, lag.max=12, main="PACF de residuos")
qqnorm(resid_ardl, main="Normal Q-Q"); qqline(resid_ardl, col=col_y)
par(mfrow=c(1,1))
lb_ardl <- Box.test(resid_ardl, lag=8, type="Ljung-Box")
lb_ardl12<- Box.test(resid_ardl, lag=12, type="Ljung-Box")
cat("Tests de diagnóstico:\n\n")
cat(sprintf(" Ljung-Box Q(8): Q=%.3f, p=%.4f %s\n",
lb_ardl$statistic, lb_ardl$p.value,
if(lb_ardl$p.value>0.05) "✓ Ruido blanco" else "⚠ Autocorrelación"))
cat(sprintf(" Ljung-Box Q(12): Q=%.3f, p=%.4f %s\n",
lb_ardl12$statistic, lb_ardl12$p.value,
if(lb_ardl12$p.value>0.05) "✓ Ruido blanco" else "⚠ Autocorrelación"))
```
---
# Sección 4: Límites del Análisis Estructural {#sec-limites}
Enders identifica dos problemas fundamentales que justifican el VAR.
## Problema 1: Parsimonia vs. Sesgo
Al podar el modelo ARDL se introduce un trade-off: más lags consumen grados de libertad, pero restricciones incorrectas sesgan los estimadores.
## Problema 2: El supuesto de no-feedback
> *"Suppose you were trying to keep a constant 70°F temperature inside your apartment by turning the thermostat up or down. If you perfectly controlled the inside temperature, there would be no correlation between the temperature and the thermostat setting."* — Enders p. 282
Si $y_t$ afecta a $z_t$ (feedback), los coeficientes de $C(L)$ están sesgados.
```{r test-granger-italia}
cat("Test de causalidad de Granger (¿hay feedback?):\n\n")
# Granger: ¿attkit Granger-causa slitaly?
cat("H0: attkit NO Granger-causa slitaly (4 lags)\n")
gc1 <- grangertest(slitaly ~ attkit, order=4)
print(gc1)
# Granger: ¿slitaly Granger-causa attkit?
cat("\nH0: slitaly NO Granger-causa attkit (4 lags)\n")
gc2 <- grangertest(attkit ~ slitaly, order=4)
print(gc2)
cat(if(gc2$`Pr(>F)`[2] > 0.05)
"→ No hay evidencia de feedback: ARDL es apropiado ✓\n"
else "→ Hay feedback: usar VAR\n")
```
---
# Sección 5: Introducción al VAR {#sec-var-intro}
## El VAR bivariate primitivo (Enders ec. 5.19-5.20)
$$\text{Full VAR: } \begin{cases} y_t = b_{10} - b_{12}z_t + \gamma_{11}y_{t-1} + \gamma_{12}z_{t-1} + u_t^y \\ z_t = b_{20} - b_{21}y_t + \gamma_{21}y_{t-1} + \gamma_{22}z_{t-1} + u_t^z \end{cases}$$
En forma matricial: $B X_t = \beta_0 + \Gamma_1 X_{t-1} + u_t$
## La forma reducida
Premultiplicando por $B^{-1}$:
$$X_t = A_0 + A_1 X_{t-1} + e_t$$
donde: $A_0 = B^{-1}\beta_0$, $A_1 = B^{-1}\Gamma_1$, $e_t = B^{-1}u_t$
**El punto crítico:** $e_t = B^{-1}u_t$ — los errores reducidos son **mezclas** de los shocks estructurales. La matriz de covarianza:
$$\Sigma = E[e_t e_t'] = B^{-1} D (B^{-1})'$$
donde $D = \text{Var}(u_t)$ es diagonal. Esta es la ecuación que usamos para la identificación.
---
# Sección 6: Estimación e Identificación {#sec-estimacion}
## Selección de lags
```{r var-lag-select}
terror_file <- file.path(DATA_PATH, "terrorism.csv")
if (!file.exists(terror_file)) terror_file <- "/tmp/terrorism.csv"
terror <- read.csv(terror_file)
terror$DATE <- as.Date(terror$DATE)
dom <- ts(terror$Domestic, start=c(1970,1), frequency=4)
trans <- ts(terror$Transnational, start=c(1970,1), frequency=4)
# Submuestra: 1979Q4 en adelante (como Enders)
dom_s <- window(dom, start=c(1979,4))
trans_s <- window(trans, start=c(1979,4))
cat("Tests ADF para terrorismo (Enders p. 310):\n")
for(nm in c("Domestic","Transnational")) {
ser <- if(nm=="Domestic") dom_s else trans_s
adf_r <- adf_manual(ser, k=3)
cat(sprintf(" %s: τ=%.3f, p %s\n", nm, adf_r$tau, adf_r$p_approx))
}
# Selección de lags: AIC y BIC manuales
Y <- cbind(as.numeric(dom_s), as.numeric(trans_s))
n_var <- nrow(Y)
cat("\nSelección de lags para el VAR:\n\n")
cat(sprintf("%-6s %-12s %-12s\n", "Lags", "AIC", "BIC"))
for(p in 1:5) {
# Construir matriz de regresores
X <- cbind(1, do.call(cbind, lapply(1:p, function(k) Y[(p-k+1):(n_var-k),])))
Yp <- Y[(p+1):n_var,]
# MCO ecuación por ecuación
B_hat <- solve(t(X) %*% X) %*% t(X) %*% Yp
E <- Yp - X %*% B_hat
Sigma <- t(E) %*% E / nrow(Yp)
k2p <- 2 * p + 1 # parámetros por ecuación (incluyendo constante)
T_eff <- nrow(Yp)
log_l <- -T_eff/2 * (log(2*pi) + log(det(Sigma)) + 2)
aic_v <- -2*log_l/T_eff + 2*(2*k2p)/T_eff
bic_v <- -2*log_l/T_eff + (2*k2p)*log(T_eff)/T_eff
cat(sprintf("p=%-4d %-12.4f %-12.4f\n", p, aic_v, bic_v))
}
cat("\n→ Enders selecciona p=3 basado en likelihood ratio\n")
```
## La descomposición de Cholesky
El problema de identificación: la forma reducida tiene 3 parámetros únicos en $\Sigma$, pero necesitamos identificar 4 ($b_{12}$, $b_{21}$, Var($\varepsilon_y$), Var($\varepsilon_z$)).
**Restricción de Cholesky** → poner $b_{12} = 0$ (o $b_{21} = 0$).
```{r cholesky-ejemplo}
cat("EJEMPLO NUMÉRICO — Pregunta 6 del libro (Enders p. 337):\n\n")
cat("Σ = [0.75 0.25]\n")
cat(" [0.25 0.50]\n\n")
Sigma_q <- matrix(c(0.75,0.25,0.25,0.50), 2, 2)
cat("Caso A: Cholesky con b12=0 (y no afecta contemporáneamente a z):\n")
b21_A <- -Sigma_q[1,2] / Sigma_q[1,1]
varE1_A <- Sigma_q[1,1]
varE2_A <- Sigma_q[2,2] + Sigma_q[2,1]*b21_A + Sigma_q[1,2]*b21_A + Sigma_q[1,1]*b21_A^2
# Correcto:
varE2_A <- Sigma_q[2,2] + Sigma_q[2,1]*b21_A
cat(sprintf(" Var(ε_y) = %.4f\n", varE1_A))
cat(sprintf(" b21 = %.4f\n", b21_A))
cat(sprintf(" Var(ε_z) = %.4f (= 5/12 = %.4f)\n\n", varE2_A, 5/12))
cat("Caso B: Cholesky con b21=0 (z no afecta contemporáneamente a y):\n")
b12_B <- -Sigma_q[1,2] / Sigma_q[2,2]
varE2_B <- Sigma_q[2,2]
varE1_B <- Sigma_q[1,1] + Sigma_q[1,2]*b12_B
cat(sprintf(" Var(ε_z) = %.4f\n", varE2_B))
cat(sprintf(" b12 = %.4f\n", b12_B))
cat(sprintf(" Var(ε_y) = %.4f (= 5/8 = %.4f)\n\n", varE1_B, 5/8))
cat("CONCLUSIÓN: Las dos ordenaciones dan resultados DISTINTOS.\n")
cat("El orden de las variables en el VAR importa.\n")
cat("Solo importa mucho si ρ₁₂ es grande.\n")
```
---
# Sección 7: La Función de Impulso-Respuesta {#sec-irf}
## La representación VMA
Todo VAR estable tiene una representación VMA:
$$X_t = \mu + \sum_{i=0}^{\infty} \Phi_i \varepsilon_{t-i}$$
Los coeficientes $\phi_{jk}(i)$ de $\Phi_i$ son los **multiplicadores**: $\phi_{jk}(i) = \partial X_{j,t+i} / \partial \varepsilon_{k,t}$.
```{r var-estimar}
# VAR(3) para los datos de terrorismo
p_var <- 3
n_v <- nrow(Y)
X_v <- cbind(1, do.call(cbind, lapply(1:p_var, function(k) Y[(p_var-k+1):(n_v-k),])))
Yp_v <- Y[(p_var+1):n_v,]
B_hat <- solve(t(X_v) %*% X_v) %*% t(X_v) %*% Yp_v
E_v <- Yp_v - X_v %*% B_hat
Sigma_v <- t(E_v) %*% E_v / (nrow(Yp_v) - ncol(X_v))
cat("VAR(3) estimado:\n\n")
cat("Ecuación Domestic:\n")
print(round(B_hat[,1], 4))
cat("\nEcuación Transnational:\n")
print(round(B_hat[,2], 4))
# Diagnóstico residuos
lb_dom <- Box.test(E_v[,1], lag=8, type="Ljung-Box")
lb_trans <- Box.test(E_v[,2], lag=8, type="Ljung-Box")
cat(sprintf("\nLjung-Box residuos Domestic: p=%.4f %s\n",
lb_dom$p.value, if(lb_dom$p.value>0.05) "✓" else "⚠"))
cat(sprintf("Ljung-Box residuos Transnational: p=%.4f %s\n",
lb_trans$p.value, if(lb_trans$p.value>0.05) "✓" else "⚠"))
```
```{r irf-var}
#| fig-cap: "IRF del VAR — respuestas a shocks de terrorismo doméstico"
# IRF manual con Cholesky
chol_S <- t(chol(Sigma_v)) # descomposición de Cholesky
# Matrices A1, A2, A3
k <- 2
A <- list()
for(j in 1:p_var) {
idx <- 1 + (j-1)*k + (1:k)
A[[j]] <- t(B_hat[idx,])
}
# Calcular matrices Phi_i
H_irf <- 16
Phi <- array(0, dim=c(k, k, H_irf+1))
Phi[,,1] <- diag(k)
for(h in 1:H_irf) {
for(j in 1:min(h, p_var)) {
Phi[,,h+1] <- Phi[,,h+1] + A[[j]] %*% Phi[,,h-j+1]
}
}
# Respuestas ortogonalizadas
IRF_ort <- array(0, dim=c(k, k, H_irf+1))
for(h in 0:H_irf) IRF_ort[,,h+1] <- Phi[,,h+1] %*% chol_S
# Graficar
nms <- c("Domestic", "Transnational")
par(mfrow=c(2,2), mar=c(4,4,3,1))
for(resp in 1:2) {
for(shock in 1:2) {
irf_vals <- IRF_ort[resp, shock, ]
plot(0:H_irf, irf_vals, type="l", lwd=2,
col=if(shock==1) col_y else col_z,
xlab="Horizonte h",
ylab="Respuesta",
main=sprintf("Shock %s → %s", nms[shock], nms[resp]))
abline(h=0, lty=2, col="gray50")
}
}
par(mfrow=c(1,1))
```
---
# Sección 8: Testing de Hipótesis {#sec-testing}
## Causalidad de Granger
> *"One test of causality is whether the lags of one variable enter into the equation for another variable."* — Enders p. 305
$\{y_t\}$ no Granger-causa $\{z_t\}$ $\iff$ todos los coeficientes de $A_{21}(L)$ son cero.
```{r granger-terrorismo}
cat("Tests de Causalidad de Granger — Terrorismo\n\n")
# Granger manual: comparar RSS con y sin los rezagos
granger_manual <- function(y, x, p=3) {
n <- length(y)
# Modelo restringido (sin x)
Xr <- do.call(cbind, lapply(1:p, function(k) y[(p-k+1):(n-k)]))
Xr <- cbind(1, Xr)
yr <- y[(p+1):n]
mr <- lm(yr ~ Xr - 1)
RSS_r <- sum(residuals(mr)^2)
# Modelo no restringido (con x)
Xu <- cbind(Xr, do.call(cbind, lapply(1:p, function(k) x[(p-k+1):(n-k)])))
mu <- lm(yr ~ Xu - 1)
RSS_u <- sum(residuals(mu)^2)
F_stat <- ((RSS_r - RSS_u)/p) / (RSS_u/(length(yr)-ncol(Xu)))
p_val <- pf(F_stat, p, length(yr)-ncol(Xu), lower.tail=FALSE)
list(F=F_stat, p=p_val, df1=p, df2=length(yr)-ncol(Xu))
}
dom_n <- as.numeric(dom_s)
trans_n <- as.numeric(trans_s)
gc_dt <- granger_manual(trans_n, dom_n, p=3)
gc_td <- granger_manual(dom_n, trans_n, p=3)
cat(sprintf("H0: Domestic NO Granger-causa Transnational (p=3):\n"))
cat(sprintf(" F(%d,%d) = %.4f, p-valor = %.4f %s\n\n",
gc_dt$df1, gc_dt$df2, gc_dt$F, gc_dt$p,
if(gc_dt$p < 0.05) "→ SÍ hay causalidad ***" else "→ No hay causalidad"))
cat(sprintf("H0: Transnational NO Granger-causa Domestic (p=3):\n"))
cat(sprintf(" F(%d,%d) = %.4f, p-valor = %.4f %s\n",
gc_td$df1, gc_td$df2, gc_td$F, gc_td$p,
if(gc_td$p < 0.05) "→ SÍ hay causalidad ***" else "→ No hay causalidad"))
cat("\nConclusión de Enders: el terrorismo doméstico y transnacional\n")
cat("se Granger-causan mutuamente → el VAR es la herramienta correcta.\n")
```
## Descomposición de Varianza del Error de Pronóstico (FEVD)
```{r fevd}
#| fig-cap: "Descomposición de varianza — VAR(3) terrorismo"
H_fevd <- 12
fevd_dom <- matrix(0, H_fevd, 2)
fevd_trans <- matrix(0, H_fevd, 2)
colnames(fevd_dom) <- c("Dom shock", "Trans shock")
colnames(fevd_trans) <- c("Dom shock", "Trans shock")
for(h in 1:H_fevd) {
mse_dom <- rowSums(IRF_ort[1, , 1:(h), drop = FALSE]^2)
mse_trans <- rowSums(IRF_ort[2, , 1:(h), drop = FALSE]^2)
fevd_dom[h,] <- 100 * mse_dom / sum(mse_dom)
fevd_trans[h,] <- 100 * mse_trans / sum(mse_trans)
}
cat("FEVD de Domestic (%):\n")
print(round(fevd_dom[c(1,4,8,12),], 2))
cat("\nFEVD de Transnational (%):\n")
print(round(fevd_trans[c(1,4,8,12),], 2))
par(mfrow=c(1,2), mar=c(4,4,3,1))
barplot(t(fevd_dom), beside=FALSE, col=c(col_y, col_z), border=NA,
names.arg=1:H_fevd, ylim=c(0,100),
xlab="Horizonte h", ylab="%",
main="FEVD — Domestic",
legend.text=c("Shock Dom.","Shock Trans."),
args.legend=list(bty="n", cex=0.8))
barplot(t(fevd_trans), beside=FALSE, col=c(col_y, col_z), border=NA,
names.arg=1:H_fevd, ylim=c(0,100),
xlab="Horizonte h", ylab="%",
main="FEVD — Transnational",
legend.text=c("Shock Dom.","Shock Trans."),
args.legend=list(bty="n", cex=0.8))
par(mfrow=c(1,1))
```
---
# Ejercicios del Capítulo 5
## Pregunta 1: Tipos de intervención (Enders p. 337)
**Enunciado:** Sea el modelo $y_t = 0.5 y_{t-1} + z_t + \varepsilon_t$. Calcula cómo responde $\{y_t\}$ a los tres tipos de intervención: pulso, salto puro y el impulso prolongado dado.
```{r ejercicio1}
cat("EJERCICIO 1 — Respuesta del modelo y_t = 0.5*y_{t-1} + z_t\n\n")
a1_q1 <- 0.5; H_q1 <- 8
# Respuesta analítica
cat("Caso (i): y_t = 0.5*y_{t-1} + z_t + eps_t\n")
cat(" Solución: y_t = sum_{k>=0} 0.5^k * z_{t-k} + eps_t/(1-0.5L)\n\n")
cat(sprintf("%-6s %-14s %-14s %-16s\n", "h", "Pulso", "Salto puro", "Imp. prolongado"))
cat(paste(rep("-",52), collapse=""), "\n")
for(h in 1:(H_q1+1)) {
# Pulso: dy_{1+h}/dz_1 = 0.5^{h-1}
pulse <- a1_q1^(h-1)
# Salto puro: suma acumulada
jump <- sum(a1_q1^(0:(h-1)))
# Impulso prolongado: z1=1, z2=0.75, z3=0.5, z4=0.25
z_prol <- c(1, 0.75, 0.5, 0.25, rep(0, H_q1))
prol <- sum(a1_q1^(0:(h-1)) * rev(z_prol[1:h]))
cat(sprintf("t=%-4d %-14.5f %-14.5f %-16.5f\n", h, pulse, jump, prol))
}
cat(sprintf("\n LP (salto puro): 1/(1-0.5) = %.4f ✓\n\n", 1/(1-a1_q1)))
```
## Pregunta 2: Transfer function con terrorismo
**Enunciado:** (Enders p. 338) Usando los datos de terrorismo transnacional de Enders, estima un ADL. ¿El terrorismo doméstico Granger-causa el transnacional?
```{r ejercicio2}
cat("EJERCICIO 2 — ADL para terrorismo transnacional\n\n")
# Identificar AR de la variable exógena (Domestic)
cat("Paso 1: AR para Domestic (z_t)\n")
aics_dom <- sapply(1:5, function(p) AIC(arima(dom_s, order=c(p,0,0))))
p_dom <- which.min(aics_dom)
ar_dom <- arima(dom_s, order=c(p_dom,0,0))
eps_dom <- residuals(ar_dom)
cat(sprintf(" Seleccionado: AR(%d) AIC=%.2f\n\n", p_dom, min(aics_dom)))
# CCF filtrada
n_t <- length(trans_s)
phi_d <- coef(ar_dom)
yf_t <- numeric(n_t)
for(t in (p_dom+1):n_t)
yf_t[t] <- trans_s[t] - sum(phi_d * rev(as.numeric(trans_s[(t-p_dom):(t-1)])))
ccf_t <- ccf(yf_t[(p_dom+1):n_t], eps_dom[1:(n_t-p_dom)], lag.max=8, plot=FALSE)
band_t <- 2/sqrt(n_t-p_dom)
cat("Paso 2: CCF filtrada (Trans^f vs eps_Dom)\n")
cat(sprintf("%-6s %-12s %-12s\n", "Lag", "CCF", "Sig?"))
for(i in which(ccf_t$lag >= 0 & ccf_t$lag <= 6)) {
cat(sprintf("%-6d %-12.4f %s\n", ccf_t$lag[i], ccf_t$acf[i],
if(abs(ccf_t$acf[i]) > band_t) "✓" else "—"))
}
# Estimar ADL
cat("\nPaso 4: ADL estimado\n")
nt <- length(trans_s); nd <- length(dom_s)
y_e <- as.numeric(trans_s[2:nt])
x_e <- as.numeric(dom_s[1:(nt-1)])
y1_e <- as.numeric(trans_s[1:(nt-1)])
ardl_q2 <- lm(y_e ~ y1_e + x_e)
print(summary(ardl_q2))
# Granger test
cat("Test de Granger: ¿Domestic Granger-causa Transnational?\n")
gc_q2 <- granger_manual(as.numeric(trans_s), as.numeric(dom_s), p=3)
cat(sprintf("F(%d,%d)=%.4f, p=%.4f %s\n",
gc_q2$df1, gc_q2$df2, gc_q2$F, gc_q2$p,
if(gc_q2$p < 0.05) "→ SÍ causa ✓" else "→ No causa"))
```
## Pregunta 6: Identificación del VAR (del Solution Manual)
**Enunciado:** $\text{Var}(e_1)=0.75$, $\text{Var}(e_2)=0.5$, $\text{Cov}(e_1,e_2)=0.25$. Muestra que no puedes identificar el VAR sin restricciones adicionales.
```{r ejercicio6}
cat("EJERCICIO 6 — Identificación del VAR estructural\n\n")
cat("Sistema de ecuaciones (de B*Σ*B' = D):\n")
cat(" Var(ε₁) = 0.75 + 0.5*b₁₂ + 0.5*b₁₂² [ec. 1]\n")
cat(" 0 = 0.25 + 0.5*b₁₂ + 0.75*b₂₁ + 0.25*b₁₂*b₂₁ [ec. 2]\n")
cat(" 0 = 0.25 + 0.5*b₁₂ + 0.75*b₂₁ + 0.25*b₁₂*b₂₁ [ec. 3 = ec. 2]\n")
cat(" Var(ε₂) = 0.5 + 0.5*b₂₁ + 0.75*b₂₁² [ec. 4]\n\n")
cat("→ 3 ecuaciones independientes, 4 incógnitas.\n")
cat("→ Sistema sub-identificado. Necesitamos 1 restricción.\n\n")
cat("Restricción A: b₁₂=0 (Cholesky — y no afecta contemporáneamente a z)\n")
b21_A <- -0.25/0.75
var1_A <- 0.75
var2_A <- 0.5 + 0.5*b21_A + 0.75*b21_A^2
cat(sprintf(" b21 = -0.25/0.75 = %.4f\n", b21_A))
cat(sprintf(" Var(ε₁) = %.4f, Var(ε₂) = %.4f (= 5/12 = %.4f)\n\n", var1_A, var2_A, 5/12))
cat("Restricción B: b₂₁=0 (Cholesky — z no afecta contemporáneamente a y)\n")
b12_B <- -0.25/0.50
var2_B <- 0.50
var1_B <- 0.75 + 0.5*b12_B + 0.5*b12_B^2
cat(sprintf(" b12 = -0.25/0.50 = %.4f\n", b12_B))
cat(sprintf(" Var(ε₂) = %.4f, Var(ε₁) = %.4f (= 5/8 = %.4f)\n\n", var2_B, var1_B, 5/8))
cat("→ Dos identificaciones diferentes → IRF diferentes.\n")
cat("→ La ORDENACIÓN de variables en el VAR importa.\n")
```
---
# Resumen del Capítulo 5
```{r resumen-final}
cat("═══════════════════════════════════════════════════════════\n")
cat(" RESUMEN CONCEPTUAL\n")
cat("═══════════════════════════════════════════════════════════\n\n")
cat("El capítulo avanza en tres etapas de supuestos sobre z_t:\n\n")
cat("┌─────────────────┬───────────────────────────────────────────┐\n")
cat("│ HERRAMIENTA │ SUPUESTO SOBRE z_t │\n")
cat("├─────────────────┼───────────────────────────────────────────┤\n")
cat("│ Intervención │ Dummy DETERMINÍSTICA — tú la construyes │\n")
cat("│ (Sec. 1) │ Exogeneidad garantizada por construcción │\n")
cat("├─────────────────┼───────────────────────────────────────────┤\n")
cat("│ ARDL / Transfer │ Serie ESTOCÁSTICA, EXÓGENA (supuesta) │\n")
cat("│ Function (Sec.2)│ Identificar C(L) con CCF pre-whitening │\n")
cat("│ │ Riesgo: sesgo si hay feedback y→z │\n")
cat("├─────────────────┼───────────────────────────────────────────┤\n")
cat("│ VAR (Sec. 5-9) │ TODAS las variables son ENDÓGENAS │\n")
cat("│ │ Sims (1980): no imposibles restricciones │\n")
cat("│ │ Identificar con Cholesky o restricciones │\n")
cat("│ │ estructurales de largo plazo (Blanchard-Q) │\n")
cat("└─────────────────┴───────────────────────────────────────────┘\n\n")
cat("FÓRMULAS CLAVE:\n\n")
cat(" Multiplicador LP (ARDL): Σcⱼ / (1 − Σaᵢ)\n")
cat(" CCVF benchmark: γ_yz(i) = c_d·a1^{i-d}·σ_z² (i ≥ d)\n")
cat(" Cholesky: n(n-1)/2 restricciones para identificar n-var VAR\n")
cat(" Granger: y no causa z ↔ todos los coefs A21(L) = 0\n")
```
