Capítulo 7 Modelos de Datos Panel, Panel VAR y Modelos No Lineales
7.1 Modelos de Datos Panel y Panel VAR
7.1.1 Motivación
Los datos panel son una extensión natural del análisis de series de tiempo en el que observamos una misma variable a lo largo del tiempo, pero para múltiples individuos. Denotemos con \(Y_{it}\) el componente o individuo \(i\), \(i = 1, 2, \ldots, N\) en el tiempo \(t\), \(t = 1, 2, \ldots, T\).
Conviene distinguir dos situaciones que aparecen recurrentemente en este capítulo, porque los métodos apropiados difieren entre ellas. En los paneles micro el número de individuos es grande y el número de periodos pequeño (\(N\) grande, \(T\) pequeño): es el caso de encuestas de hogares o de bases de empresas. En los paneles macro –los más cercanos al análisis de series de tiempo– el número de individuos es moderado y el número de periodos es relativamente largo (\(N\) moderado, \(T\) grande): países, estados o ciudades observados durante varias décadas. Las pruebas de raíz unitaria y de cointegración en panel de este capítulo están diseñadas para el segundo caso, mientras que los estimadores de panel dinámico por GMM que veremos más adelante fueron concebidos para el primero.
En general, escribiremos: \[\begin{equation} Y_{it} = \alpha_i + \beta X_{it} + U_{it} \tag{7.1} \end{equation}\]
Donde \(\beta\) es el efecto –común a todos los individuos– de la variable explicativa y \(\alpha_i\) es un término específico de cada individuo que recoge toda la heterogeneidad no observada que es constante en el tiempo. Podemos pensar en \(\alpha_i\) como el resultado de un conjunto de características fijas del individuo, observables o no: \[\begin{equation} \alpha_i = \beta_0 + \gamma' \mathbf{Z}_i \end{equation}\]
Donde \(\mathbf{Z}_i\) es un vector de características invariantes en el tiempo. La utilidad de los datos panel radica precisamente en que permiten controlar por \(\alpha_i\) sin tener que observar \(\mathbf{Z}_i\).
7.1.2 Estimadores básicos: agrupado, efectos fijos y efectos aleatorios
Antes de pasar a las pruebas de raíz unitaria y a los modelos dinámicos, conviene fijar los tres estimadores estáticos que utilizaremos a lo largo del capítulo. Todos parten de la ecuación (7.1) y se distinguen por el tratamiento que dan al término individual \(\alpha_i\):
Estimador agrupado (pooling). Ignora la estructura de panel y estima la ecuación por MCO sobre las \(N \times T\) observaciones, suponiendo que \(\alpha_i = \alpha\) para todo \(i\). Es consistente sólo si esa restricción es cierta; en caso contrario, la heterogeneidad no observada queda dentro del error y, si está correlacionada con \(X_{it}\), el estimador está sesgado. Este es el problema clásico de variable omitida en datos panel.
Estimador de efectos fijos (within). Trata a \(\alpha_i\) como un parámetro a estimar (o, equivalentemente, lo elimina). La transformación within consiste en restar a cada variable su promedio individual: \[\begin{equation} (Y_{it} - \bar{Y}_i) = \beta (X_{it} - \bar{X}_i) + (U_{it} - \bar{U}_i) \tag{7.2} \end{equation}\] Donde \(\bar{Y}_i = \frac{1}{T}\sum_t Y_{it}\). Como \(\alpha_i\) es constante en el tiempo, desaparece de la ecuación (7.2), de modo que \(\beta\) se estima consistentemente aun cuando \(\alpha_i\) esté correlacionado con los regresores. El costo es doble: no permite identificar el efecto de variables que no varían en el tiempo (quedan absorbidas por \(\alpha_i\)) y utiliza únicamente la variación dentro de cada individuo, por lo que es menos eficiente si la variación relevante es entre individuos.
Estimador de efectos aleatorios. Supone que \(\alpha_i\) es una variable aleatoria con media \(\beta_0\) y varianza \(\sigma^2_\alpha\), no correlacionada con los regresores. Bajo ese supuesto, la estimación por mínimos cuadrados generalizados aprovecha tanto la variación within como la variación between y resulta más eficiente que efectos fijos; pero si el supuesto de no correlación falla, es inconsistente.
La elección entre efectos fijos y aleatorios se resuelve, en la práctica, con la prueba de Hausman, cuya hipótesis nula es que ambos estimadores son consistentes –es decir, que \(Cov(\alpha_i, X_{it}) = 0\)– y por lo tanto que la diferencia entre sus coeficientes se debe únicamente al azar. Si se rechaza, se prefiere efectos fijos. Adicionalmente, una prueba \(F\) de significancia conjunta de los efectos individuales permite contrastar el modelo de efectos fijos contra el agrupado.
Ejemplo. Estimemos con la base Grunfeld la ecuación de inversión \(inv_{it} = \alpha_i + \beta_1 value_{it} + \beta_2 capital_{it} + U_{it}\) con los tres estimadores. Los resultados se reportan en el Cuadro 7.1.
library(plm)
data("Grunfeld", package = "plm")
mod_pool <- plm(inv ~ value + capital, data = Grunfeld,
index = c("firm", "year"), model = "pooling")
mod_fe <- plm(inv ~ value + capital, data = Grunfeld,
index = c("firm", "year"), model = "within")
mod_re <- plm(inv ~ value + capital, data = Grunfeld,
index = c("firm", "year"), model = "random")
vars_com <- c("value", "capital")
tab_panel <- data.frame(
Variable = c("$value_{it}$", "$capital_{it}$"),
Agrupado = round(coef(mod_pool)[vars_com], 4),
EF_within = round(coef(mod_fe)[vars_com], 4),
EA_random = round(coef(mod_re)[vars_com], 4))
knitr::kable(
tab_panel,
col.names = c("Variable", "Agrupado (MCO)", "Efectos fijos", "Efectos aleatorios"),
caption = "Ecuación de inversión de Grunfeld estimada con los tres estimadores de panel.",
align = c("l", "r", "r", "r"),
booktabs = TRUE,
escape = FALSE,
row.names = FALSE
)| Variable | Agrupado (MCO) | Efectos fijos | Efectos aleatorios |
|---|---|---|---|
| \(value_{it}\) | 0.1156 | 0.1101 | 0.1098 |
| \(capital_{it}\) | 0.2307 | 0.3101 | 0.3081 |
Las diferencias entre columnas son informativas: el coeficiente del capital pasa de 0.231 en el estimador agrupado a 0.31 al controlar por los efectos individuales, lo que sugiere que el modelo agrupado estaba recogiendo diferencias permanentes de escala entre empresas. El Cuadro 7.2 reporta las dos pruebas de especificación correspondientes.
prueba_F <- pFtest(mod_fe, mod_pool)
prueba_H <- phtest(mod_fe, mod_re)
knitr::kable(
data.frame(
Prueba = c("$F$ de efectos individuales",
"Hausman"),
Estadistico = round(c(as.numeric(prueba_F$statistic),
as.numeric(prueba_H$statistic)), 4),
gl = c(paste(prueba_F$parameter, collapse = ", "),
as.character(prueba_H$parameter)),
pvalor = signif(c(prueba_F$p.value, prueba_H$p.value), 4)),
col.names = c("Prueba", "Estadístico", "Grados de libertad", "Valor $p$"),
caption = paste("Pruebas de especificación del modelo de panel para la",
"ecuación de inversión de Grunfeld. La prueba $F$",
"contrasta efectos fijos contra el modelo agrupado; la de",
"Hausman, efectos fijos contra efectos aleatorios."),
align = c("l", "r", "c", "r"),
booktabs = TRUE,
escape = FALSE,
row.names = FALSE
)| Prueba | Estadístico | Grados de libertad | Valor \(p\) |
|---|---|---|---|
| \(F\) de efectos individuales | 49.1766 | 9, 188 | 0.0000 |
| Hausman | 2.3304 | 2 | 0.3119 |
La prueba \(F\) rechaza contundentemente la hipótesis de ausencia de efectos individuales, por lo que el estimador agrupado queda descartado. La prueba de Hausman, en cambio, no rechaza su hipótesis nula (valor \(p\) = 0.312), de modo que en este ejemplo los efectos aleatorios serían admisibles y más eficientes. Vale la pena notar que en la práctica aplicada se suele preferir el estimador de efectos fijos incluso en esta situación, por ser robusto a la correlación entre la heterogeneidad no observada y los regresores. En lo que resta del capítulo utilizaremos efectos fijos, tanto en el ejemplo de diferencias en diferencias como en las regresiones de cointegración en panel.
7.1.3 Pruebas de Raíces Unitarias en Panel
Las pruebas de raíces unitarias para panel suelen ser usadas principalmente en casos macroeconómicos. De forma similar al caso de series univariadas, asumiremos una forma de AR(1) para una serie en panel: \[\begin{equation} \Delta Y_{it} = \mu_i + \rho_i Y_{i,t-1} + \sum_{j = 1}^{k_i} \varphi_{ij} \Delta Y_{i,t-j} + \varepsilon_{it} \tag{7.3} \end{equation}\]
Donde \(i = 1, \ldots, N\), \(t = 1, \ldots, T\) y \(\varepsilon_{it}\) es una v.a. iid que cumple con: \[\begin{eqnarray*} \mathbb{E}[\varepsilon_{it}] & = & 0 \\ \mathbb{E}[\varepsilon_{it}^2] & = & \sigma^2_i < \infty \\ \mathbb{E}[\varepsilon_{it}^4] & < & \infty \end{eqnarray*}\]
Al igual que en el caso univariado, en este tipo de pruebas buscamos identificar cuándo las series son I(1) y cuándo I(0). Para tal efecto, la prueba de raíz unitaria que utilizaremos está basada en una prueba Dickey-Fuller Aumentada en la cual la hipótesis nula (\(H_0\)) es que todas las series en el panel son no estacionarias, es decir, son I(1). Es decir, \[\begin{equation} H_0 : \rho_i = 0 \quad \text{para todo } i = 1, \ldots, N \end{equation}\]
Por su parte, en el caso de la hipótesis alternativa tendremos dos: - \(H_1^A :\) Todas las series son I(0) – caso homogéneo, o
- \(H_1^B :\) Al menos una de las series es I(0) – caso heterogéneo
7.1.4 Ejemplo: Pruebas de Raíces Unitarias en Panel
Dependencias y configuración:
Los ejemplos utilizan bases de datos incluidas en el paquete plm, por lo que basta cargarlas con data():
Descripción de los datos: Grunfeld (1958) contiene 20 observaciones anuales de inversión bruta real (inv), valor real de la empresa (value) y valor real del capital (capital) para 10 grandes empresas de Estados Unidos durante 1935–1954. Para aplicar las pruebas de raíz unitaria necesitamos la serie de inversión de cada empresa en una columna distinta (formato ancho), que es el que espera la función purtest():
# Una columna por empresa (formato ancho)
Invest <- data.frame(split(Grunfeld$inv, Grunfeld$firm))
names(Invest) <- paste0("Firm_", 1:10)
# Series en logaritmos y en primeras diferencias
ts_LInvest <- ts(log(Invest), start = 1935, end = 1954, freq = 1)
ts_DLInvest <- diff(ts_LInvest, lag = 1, differences = 1)
str(Invest)## 'data.frame': 20 obs. of 10 variables:
## $ Firm_1 : num 318 392 411 258 331 ...
## $ Firm_2 : num 210 355 470 262 230 ...
## $ Firm_3 : num 33.1 45 77.2 44.6 48.1 74.4 113 91.9 61.3 56.8 ...
## $ Firm_4 : num 40.3 72.8 66.3 51.6 52.4 ...
## $ Firm_5 : num 39.7 50.7 74.2 53.5 42.6 ...
## $ Firm_6 : num 20.4 26 25.9 27.5 24.6 ...
## $ Firm_7 : num 24.4 23.2 32.8 32.5 26.6 ...
## $ Firm_8 : num 12.9 25.9 35 22.9 18.8 ...
## $ Firm_9 : num 26.6 23.4 30.6 20.9 28.8 ...
## $ Firm_10: num 2.54 2 2.19 1.99 2.03 1.81 2.14 1.86 0.93 1.18 ...La Figura 7.1 muestra la evolución de la inversión de las diez empresas.
matplot(1935:1954, Invest, type = "l", lty = 1, lwd = 1.5,
col = rainbow(10, v = 0.85),
xlab = "Tiempo", ylab = "Inversión bruta real")
legend("topleft", legend = names(Invest), ncol = 2, cex = 0.7,
col = rainbow(10, v = 0.85), lty = 1, lwd = 1.5)
Figura 7.1: Evolución de la inversión bruta real por empresa, panel de Grunfeld (1935-1954)
Existen varias pruebas de raíz unitaria en panel (Levin et al. 2002; Im et al. 2003; Maddala and Wu 1999; Hadri 2000). Apliquemos primero la de Levin, Lin y Chu (2002), que impone un coeficiente autorregresivo común a todas las series (hipótesis alternativa homogénea, \(H_1^A\)), sobre el logaritmo de la inversión en niveles y en primeras diferencias:
## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: ts_LInvest
## z = -0.63214, p-value = 0.2636
## alternative hypothesis: stationarity## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: ts_DLInvest
## z = -8.5076, p-value < 0.00000000000000022
## alternative hypothesis: stationarityEn niveles el valor \(p\) es superior a 0.10, por lo que no se rechaza la hipótesis nula de raíz unitaria en el panel; en primeras diferencias se rechaza de manera contundente. Nótese la advertencia que emite purtest() sobre la longitud de las series (\(T = 20\)): en paneles tan cortos la selección automática de rezagos y las aproximaciones asintóticas de la prueba deben tomarse con reserva.
Consideremos ahora la prueba de Im, Pesaran y Shin (2003), que promedia las pruebas ADF individuales y por lo tanto admite que el coeficiente autorregresivo difiera entre empresas (hipótesis alternativa heterogénea, \(H_1^B\)):
##
## Im-Pesaran-Shin Unit-Root Test (ex. var.: Individual
## Intercepts)
##
## data: ts_LInvest
## Wtbar = 0.6833, p-value = 0.7528
## alternative hypothesis: stationarity##
## Im-Pesaran-Shin Unit-Root Test (ex. var.: Individual
## Intercepts)
##
## data: ts_DLInvest
## Wtbar = -11.14, p-value < 0.00000000000000022
## alternative hypothesis: stationarityLa conclusión coincide con la de la prueba anterior: el panel de inversión es \(I(1)\). Conviene subrayar la diferencia en la interpretación de un rechazo en cada prueba. En la prueba de Levin, Lin y Chu el rechazo significa que todas las series del panel son estacionarias, mientras que en la de Im, Pesaran y Shin sólo permite afirmar que al menos una lo es. Por ello, cuando se sospecha que la velocidad de reversión a la media difiere entre individuos, la prueba IPS es la especificación adecuada, aunque su conclusión sea menos informativa.
7.1.5 Ejemplo: Diferencias en Diferencias en el INPC de servicios de análisis clínicos
Este ejemplo modela el efecto que tiene que una empresa de laboratorios de análisis clínicos compre a sus competidores en diferentes ciudades del país. Para lo cual usaremos información del Índice Nacional de Precios al Consumidor de los servicios de análisis clínicos.
Dependencias y configuración:
Importamos Datos desde un archivo de Excel. Los datos están en formato panel e incluyen: - INPC (Índice Nacional de Precios al Consumidor (2QJul2018 = 100)) para un tipo de servicio. Fuente: https://www.inegi.org.mx/programas/inpc/2018/#Tabulados
Variable dicotómica que indica el momento en que inició el tratamiento
Variable dicotómica que identifica al grupo expuesto al tratamiento
Interacción entre el tiempo y el grupo tratado, además de una variable de tendencia
Otras variables de control
Partamos de la siguiente ecuación: \[\begin{equation} y_{it} = \beta_1 + \beta_2 T_t + \beta_3 D_i + \theta T_t \times D_i + \boldsymbol{\gamma} \mathbf{X}_{it} + \varepsilon_{it} \tag{7.4} \end{equation}\]
Donde \(y_{it}\) es la variable sobre la cual queremos evaluar el efecto de un tratamiento, \(\mathbf{X}_{it}\) es un conjunto de variables de control, que incluye a las variables de efectos fijos (temporales, de individuos o del tipo que se decida) y, en su caso, una tendencia, entre otras. Finalmente, las variables \(T_t\) y \(D_i\) son variables dicotómicas que identifican con el valor de 1 el momento a partir del cual se observa el tratamiento y con 0 en cualquier otro caso, y con 1 los individuos que fueron tratados y con 0 en cualquier otro caso, respectivamente.
De esta forma, el producto \(T_t \times D_i\) sería una variable dicotómica que indicaría con un 1 a aquellos individuos que fueron tratados a partir del momento en que se implementó el tratamiento y con 0 en cualquier otro caso. Bajo este escenario, \(\theta\) es un coeficiente que captura el efecto del tratamiento, condicional a que se ha controlado por una serie de factores.
Otra forma de ver lo que se indica en la ecuación (7.4) es como la comparación de la diferencia en diferencia descrita por: \[\begin{eqnarray} E & = & [ (\overline{y}_{exit} | treatment) - (\overline{y}_{baseline} | treatment) ] \nonumber \\ & & - [ (\overline{y}_{exit} | placebo) - (\overline{y}_{baseline} | placebo) ] \tag{7.5} \end{eqnarray}\]
La expresión (7.5) cuantifica el efecto que tiene el tratamiento como la diferencia de medias de dos grupos: un grupo tratado y otro no tratado. Dentro de los cuales se han comparado las medias respecto de una línea base. Podemos plantear esto desde esta otra perspectiva. Pensemos que solo tenemos dos periodos \(t = {1, 2}\), en el cual en 1 no ha ocurrido el tratamiento y en 2 ya ha ocurrido. Así, el cambio en la variable de respuesta para los individuos tratados será (suponiendo que omitimos la matriz \(\mathbf{X}_{it}\) en la ecuación (7.4)): \[\begin{eqnarray} \mathbb{E}[ y_{i2} | D_i = 1 ] - \mathbb{E}[ y_{i1} | D_i = 1 ] & = & ( \beta_1 + \beta_2 + \beta_3 + \theta ) - ( \beta_1 + \beta_3 ) \nonumber \\ & = & \beta_2 + \theta \tag{7.6} \end{eqnarray}\]
Ahora, sobre los no tratados: \[\begin{eqnarray} \mathbb{E}[ y_{i2} | D_i = 0 ] - \mathbb{E}[ y_{i1} | D_i = 0 ] & = & ( \beta_1 + \beta_2 ) - ( \beta_1 ) \nonumber \\ & = & \beta_2 \tag{7.7} \end{eqnarray}\]
Así, la diferencia en diferencia resulta de restar la ecuación (7.7) a la ecuación (7.6): \[\begin{eqnarray} & & ( \mathbb{E}[ y_{i2} | D_i = 1 ] - \mathbb{E}[ y_{i1} | D_i = 1 ] ) - ( \mathbb{E}[ y_{i2} | D_i = 0 ] - \mathbb{E}[ y_{i1} | D_i = 0 ] ) \nonumber \\ & & = \theta \tag{7.8} \end{eqnarray}\]
En nuestro caso, utilizaremos como tratamiento la entrada al mercado de un laboratorio de análisis clínicos en una ciudad determinada.
## [1] 13420 17## # A tibble: 6 × 17
## Fecha Mes Anio Ciudad INPC Dummy_FIRMA Marca_01 Marca_02
## <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl>
## 1 Jul 2002 7 2002 Área Metro… 68.3… 0 0 0
## 2 Ago 2002 8 2002 Área Metro… 68.3… 0 0 0
## 3 Sep 2002 9 2002 Área Metro… 68.3… 0 0 0
## 4 Oct 2002 10 2002 Área Metro… 68.3… 0 0 0
## 5 Nov 2002 11 2002 Área Metro… 68.3… 0 0 0
## 6 Dic 2002 12 2002 Área Metro… 68.3… 0 0 0
## # ℹ 9 more variables: Marca_03 <dbl>, Marca_04 <dbl>, Marca_05 <dbl>,
## # Marca_06 <dbl>, Marca_07 <dbl>, Marca_08 <dbl>, Marca_09 <dbl>,
## # Marca_10 <dbl>, Marca_11 <dbl>Selección de ciudades:
Data <- Data[ which( Data$Ciudad != 'Atlacomulco, Edo. de Méx.' &
Data$Ciudad != 'Cancún, Q.R.' &
Data$Ciudad != 'Coatzacoalcos, Ver.' &
Data$Ciudad != 'Esperanza, Son.' &
Data$Ciudad != 'Izúcar de Matamoros, Pue.' &
Data$Ciudad != 'Pachuca, Hgo.' &
Data$Ciudad != 'Tuxtla Gutiérrez, Chis.' &
Data$Ciudad != 'Zacatecas, Zac.' ), ]Transformaciones:
# Log de INPC
Data$LINPC <- log( as.numeric(Data$INPC) , base = exp(1))
# Volvemos factor (ordenado) a la fecha. Los niveles se construyen de
# forma programática --de julio de 2002 al último mes disponible-- para
# evitar errores de captura en una lista larga de etiquetas.
meses_abr <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun",
"Jul", "Ago", "Sep", "Oct", "Nov", "Dic")
niveles_periodo <- paste(rep(meses_abr, times = 21),
rep(2002:2022, each = 12))
niveles_periodo <- niveles_periodo[
seq(which(niveles_periodo == "Jul 2002"),
which(niveles_periodo == "Oct 2022"))]
# Verificación: todas las fechas de la base deben tener un nivel asignado
stopifnot(all(Data$Fecha %in% niveles_periodo))
Data$Periodo <- factor(Data$Fecha, order = TRUE,
levels = niveles_periodo)La Figura 7.2 muestra la evolución del logaritmo del INPC de análisis clínicos para las 47 ciudades del panel (líneas grises); la línea roja es el promedio simple entre ciudades. Se aprecia una tendencia común al alza y una dispersión entre ciudades que se reduce hacia el final de la muestra.
Data %>%
ggplot( aes(x = Periodo, y = LINPC, group = Ciudad) ) +
geom_line(color = "grey70", linewidth = 0.3) +
stat_summary(aes(group = 1), fun = mean, geom = "line",
color = "darkred", linewidth = 0.9) +
scale_x_discrete(breaks = levels(Data$Periodo)[
seq(1, nlevels(Data$Periodo), by = 12)]) +
labs(x = NULL, y = "log(INPC) de análisis clínicos") +
theme(axis.text.x = element_text( size = 8, angle = 90, vjust = 0.5))
Figura 7.2: Logaritmo del INPC de análisis clínicos por ciudad (líneas grises) y promedio entre ciudades (línea roja)
Hagamos una prueba de Raíces Unitarias.
Conversión a formato ancho:
## # A tibble: 6 × 48
## Fecha `Acapulco, Gro.` `Aguascalientes, Ags.` Área Metropolitana d…¹
## <chr> <dbl> <dbl> <dbl>
## 1 Abr … 4.39 4.35 4.26
## 2 Abr … 4.39 4.36 4.31
## 3 Abr … 4.44 4.36 4.36
## 4 Abr … 4.47 4.33 4.40
## 5 Abr … 4.54 4.38 4.42
## 6 Abr … 4.60 4.41 4.44
## # ℹ abbreviated name: ¹`Área Metropolitana de la Cd. de México`
## # ℹ 44 more variables: `Campeche, Camp.` <dbl>,
## # `Cd. Acuña, Coah.` <dbl>, `Cd. Jiménez, Chih.` <dbl>,
## # `Cd. Juárez, Chih.` <dbl>, `Chetumal, Q.R.` <dbl>,
## # `Chihuahua, Chih.` <dbl>, `Colima, Col.` <dbl>,
## # `Córdoba, Ver.` <dbl>, `Cortazar, Gto.` <dbl>,
## # `Cuernavaca, Mor.` <dbl>, `Culiacán, Sin.` <dbl>, …Series de tiempo:
ts_LINPC <- ts( LINPC[c( 2:48 )],
start = c(2002, 7),
freq = 12)
ts_DLINPC <- diff(ts( LINPC[c( 2:48 )],
start = c(2002, 7),
freq = 12),
lag = 1,
differences = 1)Prueba de raíz unitaria en panel:
Rezagos óptimos: \[\begin{eqnarray*} p & = & Int(4*(T/100)^{(1/4)}) \\ & = & Int(4*(244/100)^{(1/4)}) \\ & = & Int(4.99) \\ & = & 4 \end{eqnarray*}\]
Prueba de Levin et al. (2002):
purtest(ts_LINPC, test = "levinlin", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts
## and Trend)
##
## data: ts_LINPC
## z = -56.022, p-value < 2.2e-16
## alternative hypothesis: stationaritypurtest(ts_DLINPC, test = "levinlin", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts
## and Trend)
##
## data: ts_DLINPC
## z = -107.84, p-value < 2.2e-16
## alternative hypothesis: stationarityPrueba de Choi (2001):
purtest(ts_LINPC, test = "logit", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Choi's Logit Unit-Root Test (ex. var.: Individual Intercepts
## and Trend)
##
## data: ts_LINPC
## L* = -85.506, df = 239, p-value < 2.2e-16
## alternative hypothesis: stationaritypurtest(ts_DLINPC, test = "logit", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Choi's Logit Unit-Root Test (ex. var.: Individual Intercepts
## and Trend)
##
## data: ts_DLINPC
## L* = -278.56, df = 239, p-value < 2.2e-16
## alternative hypothesis: stationarityPrueba de Hadri (2000):
purtest(ts_LINPC, test = "hadri", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Hadri Test (ex. var.: Individual Intercepts and Trend)
## (Heterosked. Consistent)
##
## data: ts_LINPC
## z = -2.7, p-value = 0.9965
## alternative hypothesis: at least one series has a unit rootpurtest(ts_DLINPC, test = "hadri", exo = "trend", # exo = c("none", "intercept", "trend"),
lags = "AIC", pmax = 4)##
## Hadri Test (ex. var.: Individual Intercepts and Trend)
## (Heterosked. Consistent)
##
## data: ts_DLINPC
## z = -8.9351, p-value = 1
## alternative hypothesis: at least one series has a unit rootEstimación del efecto del tratamiento. Estimamos la ecuación
(7.4) por efectos fijos de ciudad, incluyendo una tendencia
temporal. La variable Dummy_FIRMA es la interacción
\(T_t \times D_i\): toma el valor de 1 en las ciudades y meses en los que
la empresa ya había adquirido a su competidor, por lo que su coeficiente
es el estimador de diferencias en diferencias, \(\hat{\theta}\). Se
estiman dos especificaciones: una que incluye controles por la presencia
de cada marca en la ciudad y otra que sólo incluye el tratamiento y la
tendencia.
# Especificación 1: con controles por marca presente en la ciudad
fixed_1 <- plm( LINPC ~ Dummy_FIRMA + Marca_01 + Marca_02 + Marca_03 + Marca_04 + Marca_05 +
Marca_06 + Marca_07 + Marca_08 + Marca_09 + Marca_10 + Marca_11 + as.numeric(Periodo),
data = Data,
index=c("Ciudad", "Periodo"),
model = "within")
# Especificación 2: sólo tratamiento y tendencia
fixed_2 <- plm( LINPC ~ Dummy_FIRMA + as.numeric(Periodo),
data = Data,
index=c("Ciudad", "Periodo"),
model = "within")
did_tab <- data.frame(
Especificacion = c("Con controles de marca", "Sólo tratamiento y tendencia"),
Theta = round(c(coef(fixed_1)["Dummy_FIRMA"], coef(fixed_2)["Dummy_FIRMA"]), 5),
SE = round(c(sqrt(diag(vcov(fixed_1)))["Dummy_FIRMA"],
sqrt(diag(vcov(fixed_2)))["Dummy_FIRMA"]), 5),
t = round(c(coef(summary(fixed_1))["Dummy_FIRMA", 3],
coef(summary(fixed_2))["Dummy_FIRMA", 3]), 3),
p = signif(c(coef(summary(fixed_1))["Dummy_FIRMA", 4],
coef(summary(fixed_2))["Dummy_FIRMA", 4]), 4))
knitr::kable(
did_tab,
col.names = c("Especificación", "$\\hat{\\theta}$", "Error Est.", "Estad. $t$", "Valor $p$"),
caption = "Estimador de diferencias en diferencias del efecto de la adquisición sobre el logaritmo del INPC de análisis clínicos (efectos fijos por ciudad y tendencia temporal).",
align = c("l", "r", "r", "r", "r"),
booktabs = TRUE,
escape = FALSE,
row.names = FALSE
)| Especificación | \(\hat{\theta}\) | Error Est. | Estad. \(t\) | Valor \(p\) |
|---|---|---|---|---|
| Con controles de marca | -0.06483 | 0.01851 | -3.502 | 0.0004629 |
| Sólo tratamiento y tendencia | -0.02255 | 0.00413 | -5.460 | 0.0000000 |
El coeficiente estimado es negativo y estadísticamente significativo en ambas especificaciones. Dado que la variable dependiente está en logaritmos, un \(\hat{\theta}\) de -0.0225 implica que, en las ciudades tratadas y a partir del momento de la adquisición, el índice de precios de los servicios de análisis clínicos fue aproximadamente 2.25% menor que lo que sugeriría la trayectoria del grupo de control. El signo y la significancia del efecto son robustos a la inclusión de los controles por marca, pero su magnitud no lo es: con controles el efecto estimado es de 6.48%, casi el triple. Esa sensibilidad indica que parte del efecto que la especificación sencilla atribuye a la adquisición se explica por la composición de marcas presentes en cada ciudad, por lo que la especificación con controles es la preferible.
Dos advertencias son necesarias para interpretar este ejercicio. Primera, la validez del estimador depende del supuesto de tendencias paralelas: en ausencia del tratamiento, las ciudades tratadas y no tratadas habrían seguido trayectorias de precios paralelas. Segunda, cuando el tratamiento se adopta en momentos distintos según la unidad –como ocurre aquí, pues las adquisiciones no fueron simultáneas– la estimación por efectos fijos con una sola variable de interacción recupera un promedio ponderado de efectos que puede diferir del efecto promedio del tratamiento; la literatura reciente sobre DiD con adopción escalonada propone estimadores alternativos para ese caso.
7.1.6 Panel VAR
En esta sección extenderemos el caso del modelo VAR(p) a uno en forma panel. En este caso asumimos que las variables exógenas son los \(p\) rezagos de las \(k\) variables endógenas. Consideremos el siguiente caso de un modelo panel VAR con efectos fijos –ciertamente es posible hacer estimaciones con efectos aleatorios, no obstante, requiere de supuestos adicionales que no contemplamos en este libro–, el cual es la forma más común de estimación: \[\begin{equation} \mathbf{Y}_{it} = \mu_i + \sum_{l = 1}^p \mathbf{A}_l \mathbf{Y}_{i t - l} + \mathbf{B} \mathbf{X}_{it} + \varepsilon_{it} \tag{7.9} \end{equation}\]
Donde \(\mathbf{Y}_{it}\) es un vector de variables endógenas estacionarias para la unidad de corte transversal \(i\) en el tiempo \(t\), \(\mathbf{X}_{it}\) es una matriz que contiene variables exógenas, y \(\varepsilon_{it}\) es un término de error que cumple con: \[\begin{eqnarray*} \mathbb{E}[\varepsilon_{it}] & = & 0 \\ Var[\varepsilon_{it}] & = & \Sigma_\varepsilon \end{eqnarray*}\]
Donde \(\Sigma_\varepsilon\) es una matriz definida positiva.
7.1.6.1 ¿Por qué no estimar por MCO? El sesgo de Nickell
A diferencia del modelo panel estático, en la ecuación (7.9) los regresores incluyen rezagos de la variable dependiente. Esto invalida a los estimadores tradicionales de panel: el estimador within (de efectos fijos) elimina \(\mu_i\) restando la media individual, pero la media de \(\mathbf{Y}_{it}\) contiene información de todos los periodos, incluidos los errores pasados, por lo que el regresor transformado queda correlacionado con el error transformado. El sesgo resultante, documentado por Nickell (1981), es de orden \(1/T\): resulta despreciable en paneles largos, pero considerable en los paneles macro y microeconómicos típicos, donde \(T\) es pequeño (\(T \leq 10\), por ejemplo) aunque \(N\) sea grande. Por ello, la literatura de paneles dinámicos estima el PVAR por el método generalizado de momentos (GMM) sobre una transformación del modelo que elimine los efectos fijos sin inducir esa correlación.
7.1.6.2 Transformaciones: primeras diferencias y desviaciones ortogonales
La transformación más directa es la de primeras diferencias (FD): \[\begin{equation} \Delta \mathbf{Y}_{it} = \sum_{l = 1}^p \mathbf{A}_l \Delta \mathbf{Y}_{i t - l} + \mathbf{B} \Delta \mathbf{X}_{it} + \Delta \varepsilon_{it} \tag{7.10} \end{equation}\]
que elimina \(\mu_i\), pero con dos costos: (i) el nuevo error \(\Delta \varepsilon_{it} = \varepsilon_{it} - \varepsilon_{i,t-1}\) sigue por construcción un proceso MA(1), de modo que está correlacionado con \(\Delta \mathbf{Y}_{i,t-1}\) y la estimación requiere instrumentos; y (ii) cada hueco en el panel (dato faltante) elimina dos observaciones transformadas.
Una alternativa es la transformación de desviaciones ortogonales hacia adelante (FOD, forward orthogonal deviations) de Arellano y Bover (1995): a cada observación se le resta el promedio de las observaciones futuras disponibles de la misma unidad, reescalando para preservar la varianza: \[\begin{equation} \tilde{\mathbf{Y}}_{it} = c_{it} \left[ \mathbf{Y}_{it} - \frac{1}{T - t} \left( \mathbf{Y}_{i,t+1} + \ldots + \mathbf{Y}_{iT} \right) \right], \quad c_{it} = \sqrt{\frac{T-t}{T-t+1}}. \tag{7.11} \end{equation}\]
Como solo usa información futura, la transformación no introduce autocorrelación en los errores transformados (si \(\varepsilon_{it}\) no está autocorrelacionado, \(\tilde{\varepsilon}_{it}\) tampoco lo está) y solo se pierde la última observación de cada unidad, lo que la hace preferible en paneles cortos o con huecos. En el paquete panelvar de R, el argumento transformation de la función pvargmm() permite elegir entre ambas: "fd" (primeras diferencias) o "fod" (desviaciones ortogonales).
7.1.6.3 Estimación GMM: en diferencias y en sistema
Una vez transformado el modelo, el estimador GMM en diferencias de Arellano y Bond (1991) explota como instrumentos los niveles rezagados de las variables: bajo el supuesto de que los errores no están autocorrelacionados, los niveles \(\mathbf{Y}_{i,t-s}\) con \(s \geq 2\) no están correlacionados con \(\Delta \varepsilon_{it}\), lo que genera las condiciones de momentos \[\begin{equation} \mathbb{E}\left[ \mathbf{Y}_{i,t-s} \, \Delta \varepsilon_{it}' \right] = 0, \quad s = 2, 3, \ldots, t-1. \tag{7.12} \end{equation}\]
Cuando las series son muy persistentes (cercanas a raíz unitaria), los niveles rezagados son instrumentos débiles de las diferencias y el GMM en diferencias se comporta mal. El estimador GMM en sistema de Blundell y Bond (1998) añade al sistema las ecuaciones en niveles, instrumentadas con las diferencias rezagadas, bajo un supuesto adicional sobre las condiciones iniciales del proceso (que las desviaciones iniciales respecto del estado estacionario no estén correlacionadas con los efectos fijos). En pvargmm(), esta variante se activa con system_instruments = TRUE.
Dos aspectos prácticos completan la especificación:
- Proliferación de instrumentos. El número de instrumentos de la ecuación (7.12) crece con \(T^2\) y con el número de variables. Demasiados instrumentos sobreajustan las variables endógenas y debilitan las pruebas de especificación. Roodman (2009) recomienda limitar los rezagos usados como instrumentos (argumentos
max_instr_dependent_varsymin_instr_dependent_vars) o colapsar la matriz de instrumentos (collapse = TRUE). - Una y dos etapas. El estimador de dos etapas (
steps = "twostep") utiliza la matriz de ponderación óptima estimada en la primera etapa y es más eficiente, aunque sus errores estándar deben corregirse en muestras pequeñas (Windmeijer 2005).
7.1.6.4 Diagnóstico y selección de rezagos
La validez del conjunto de instrumentos se evalúa con la prueba de sobreidentificación de Hansen (\(J\)), cuya hipótesis nula es que los instrumentos son válidos (no correlacionados con el error); valores \(p\) demasiado bajos indican instrumentos inválidos, pero valores \(p\) sospechosamente cercanos a 1 suelen ser síntoma de demasiados instrumentos. Adicionalmente, se examina la autocorrelación de los errores en diferencias: se espera autocorrelación de orden 1 (por construcción del MA(1)), pero no de orden 2 — si la hubiera, los instrumentos con \(s = 2\) dejarían de ser válidos. Finalmente, el orden \(p\) del PVAR se selecciona con los criterios de Andrews y Lu (2001) (MMSC-AIC, MMSC-BIC, MMSC-HQIC), que son análogos a los criterios de información tradicionales pero se construyen a partir del estadístico \(J\) de Hansen; los calcularemos con la función Andrews_Lu_MMSC() en los ejemplos.
7.2 Ejemplos: Panel VAR(p)
Dependencias y configuración:
7.2.1 Ejemplo 1
La estimación sigue la literatura de panel dinámico de Arellano y Bond (1991), Blundell y Bond (1998) y Roodman (2009). Este conjunto de datos describe el empleo, los salarios, el capital y la producción de 140 empresas en el Reino Unido durante 1976–1984. Se estima un modelo en el que el empleo (\(n\)) es explicado por sus valores pasados (2 rezagos), los valores actuales y el primer rezago de los salarios (\(w\)) y la producción (\(ys\)), y el valor actual del capital (\(k\)).
## 'data.frame': 1031 obs. of 30 variables:
## $ c1 : chr "1-1" "2-1" "3-1" "4-1" ...
## $ ind : num 7 7 7 7 7 7 7 7 7 7 ...
## $ year : Factor w/ 9 levels "1976","1977",..: 2 3 4 5 6 7 8 2 3 4 ...
## $ emp : num 5.04 5.6 5.01 4.72 4.09 ...
## $ wage : num 13.2 12.3 12.8 13.8 14.3 ...
## $ cap : num 0.589 0.632 0.677 0.617 0.508 ...
## $ indoutpt: num 95.7 97.4 99.6 100.6 99.6 ...
## $ n : num 1.62 1.72 1.61 1.55 1.41 ...
## [list output truncated]Estimación. Replicamos la especificación del Cuadro 4b de Arellano y
Bond (1991): transformación en primeras diferencias
(transformation = "fd"), GMM en dos etapas y los niveles rezagados de
\(n\) desde el segundo rezago como instrumentos.
#?pvargmm
Arellano_Bond_1991_table4b <- pvargmm( dependent_vars = c("n"),
lags = 2,
exog_vars = c("w", "wL1", "k", "ys", "ysL1", "yr1979", "yr1980", "yr1981", "yr1982",
"yr1983", "yr1984"),
transformation = "fd", data = abdata, panel_identifier = c("id", "year"),
steps = c("twostep"),
system_instruments = FALSE,
max_instr_dependent_vars = 99,
min_instr_dependent_vars = 2L,
collapse = FALSE,
# La barra de avance de la estimación
# imprime una línea por iteración: útil
# en la consola, ilegible en el libro.
progressbar = FALSE)## ---------------------------------------------------
## Dynamic Panel VAR estimation, two-step GMM
## ---------------------------------------------------
## Transformation: First-differences
## Group variable: id
## Time variable: year
## Number of observations = 611
## Number of groups = 140
## Obs per group: min = 4
## avg = 4.364286
## max = 6
## Number of instruments = 38
##
## ===================
## n
## -------------------
## lag1_n 0.4742 *
## (0.1854)
## lag2_n -0.0530
## (0.0517)
## w -0.5132 ***
## (0.1456)
## wL1 0.2246
## (0.1419)
## k 0.2927 ***
## (0.0626)
## ys 0.6098 ***
## (0.1563)
## ysL1 -0.4464 *
## (0.2173)
## yr1979 0.0105
## (0.0099)
## yr1980 0.0247
## (0.0158)
## yr1981 -0.0158
## (0.0267)
## yr1982 -0.0374
## (0.0300)
## yr1983 -0.0393
## (0.0347)
## yr1984 -0.0495
## (0.0349)
## ===================
## *** p < 0.001; ** p < 0.01; * p < 0.05
##
## ---------------------------------------------------
## Instruments for equation
## Standard
## FD.(w wL1 k ys ysL1 yr1979 yr1980 yr1981 yr1982 yr1983 yr1984)
## GMM-type
## Dependent vars: L(2, 6)
## Collapse = FALSE
## ---------------------------------------------------
##
## Hansen test of overid. restrictions: chi2(25) = 30.11 Prob > chi2 = 0.22
## (Robust, but weakened by many instruments.)La lectura de los resultados sigue la lógica del panel dinámico. El coeficiente del primer rezago del empleo es positivo y significativo, lo que confirma una persistencia elevada del empleo a nivel de empresa, mientras que el segundo rezago no resulta significativo. Los salarios tienen el efecto negativo que anticipa la teoría de la demanda de trabajo y la producción, un efecto positivo. Respecto del diagnóstico, la prueba de sobreidentificación de Hansen arroja un estadístico de 30.112 con un valor \(p\) de 0.22: no se rechaza la validez del conjunto de instrumentos, aunque conviene recordar la advertencia de Roodman (2009) sobre la proliferación de instrumentos –aquí se permiten hasta 99 rezagos como instrumentos–, que tiende a debilitar justamente esta prueba.
7.2.2 Ejemplo 2:
Se utiliza un conjunto de datos panel de 265 municipios suecos durante 9 años (1979–1987). Las variables incluyen el gasto total (expenditures), los ingresos propios (revenues) y las transferencias intergubernamentales recibidas por el municipio (grants). Fuente: Dahlberg y Johansson (2000). Las transferencias del gobierno central al gobierno local son de tres tipos: apoyo a municipios con baja capacidad fiscal, transferencias para financiar ciertas actividades del gobierno local, y transferencias destinadas a inversiones específicas.
## [1] "id" "year" "expenditures" "revenues"
## [5] "grants"## id year expenditures revenues grants
## 1 114 1979 0.0229736 0.0181770 0.0054429
## 2 114 1980 0.0266307 0.0209142 0.0057304
## 3 114 1981 0.0273253 0.0210836 0.0056647
## 4 114 1982 0.0288704 0.0234310 0.0058859
## 5 114 1983 0.0226474 0.0179979 0.0055908
## 6 114 1984 0.0215601 0.0179949 0.0047536Estimación:
ex1_dahlberg_data <- pvargmm(dependent_vars = c("expenditures", "revenues", "grants"),
lags = 1,
transformation = "fod",
data = Dahlberg,
panel_identifier=c("id", "year"),
steps = c("twostep"),
system_instruments = FALSE,
max_instr_dependent_vars = 99,
max_instr_predet_vars = 99,
min_instr_dependent_vars = 2L,
min_instr_predet_vars = 1L,
collapse = FALSE,
progressbar = FALSE
)## Warning in pvargmm(dependent_vars = c("expenditures", "revenues",
## "grants"), : The matrix D_e is singular, therefore the general
## inverse is used## ---------------------------------------------------
## Dynamic Panel VAR estimation, two-step GMM
## ---------------------------------------------------
## Transformation: Forward orthogonal deviations
## Group variable: id
## Time variable: year
## Number of observations = 1855
## Number of groups = 265
## Obs per group: min = 7
## avg = 7
## max = 7
## Number of instruments = 252
##
## =========================================================
## expenditures revenues grants
## ---------------------------------------------------------
## lag1_expenditures 0.2846 *** 0.2583 ** 0.0167
## (0.0664) (0.0795) (0.0172)
## lag1_revenues -0.0470 0.0588 -0.0405 **
## (0.0637) (0.0726) (0.0151)
## lag1_grants -1.6746 *** -2.2367 *** 0.3204 ***
## (0.2818) (0.2846) (0.0521)
## =========================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
##
## ---------------------------------------------------
## Instruments for equation
## Standard
##
## GMM-type
## Dependent vars: L(2, 7)
## Collapse = FALSE
## ---------------------------------------------------
##
## Hansen test of overid. restrictions: chi2(243) = 263.01 Prob > chi2 = 0.18
## (Robust, but weakened by many instruments.)Los resultados muestran un sistema en el que las tres variables se determinan conjuntamente, lo que justifica tratarlas como endógenas en lugar de estimar tres ecuaciones separadas. En particular, el rezago de las transferencias entra con signo negativo y muy significativo en las ecuaciones de gasto y de ingresos –un municipio que recibió más transferencias el año anterior gasta y recauda menos este año–, mientras que el rezago del gasto es significativo en la ecuación de ingresos. En cambio, el rezago de los ingresos no resulta significativo ni en su propia ecuación ni en la del gasto.
Selección del orden de rezagos. Aplicamos los criterios MMSC de Andrews y Lu (2001), que se construyen a partir del estadístico \(J\) de Hansen y penalizan el número de parámetros y de instrumentos. Compararemos el PVAR(1) estimado con un PVAR(2):
## $MMSC_BIC
## [1] -1610.877
##
## $MMSC_AIC
## [1] -234.9924
##
## $MMSC_HQIC
## [1] -792.3698ex2_dahlberg_data <- pvargmm(dependent_vars = c("expenditures", "revenues", "grants"),
lags = 2,
transformation = "fod",
data = Dahlberg,
panel_identifier=c("id", "year"),
steps = c("twostep"),
system_instruments = FALSE,
max_instr_dependent_vars = 99,
max_instr_predet_vars = 99,
min_instr_dependent_vars = 2L,
min_instr_predet_vars = 1L,
collapse = FALSE,
progressbar = FALSE)## Warning in pvargmm(dependent_vars = c("expenditures", "revenues",
## "grants"), : The matrix D_e is singular, therefore the general
## inverse is used## ---------------------------------------------------
## Dynamic Panel VAR estimation, two-step GMM
## ---------------------------------------------------
## Transformation: Forward orthogonal deviations
## Group variable: id
## Time variable: year
## Number of observations = 1590
## Number of groups = 265
## Obs per group: min = 6
## avg = 6
## max = 6
## Number of instruments = 243
##
## =======================================================
## expenditures revenues grants
## -------------------------------------------------------
## lag1_expenditures 0.2572 ** 0.2486 * 0.0135
## (0.0893) (0.1018) (0.0178)
## lag1_revenues -0.1219 -0.0634 -0.0293
## (0.0879) (0.0997) (0.0171)
## lag1_grants -3.0718 *** -3.5849 *** 0.1581 *
## (0.5941) (0.6168) (0.0636)
## lag2_expenditures -0.0247 0.0252 0.0178
## (0.0791) (0.0834) (0.0157)
## lag2_revenues -0.2584 ** -0.2446 ** -0.0237
## (0.0785) (0.0800) (0.0157)
## lag2_grants -1.5265 *** -1.6687 *** 0.0702
## (0.1873) (0.2113) (0.0586)
## =======================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
##
## ---------------------------------------------------
## Instruments for equation
## Standard
##
## GMM-type
## Dependent vars: L(2, 6)
## Collapse = FALSE
## ---------------------------------------------------
##
## Hansen test of overid. restrictions: chi2(225) = 260.11 Prob > chi2 = 0.054
## (Robust, but weakened by many instruments.)## $MMSC_BIC
## [1] -1486.935
##
## $MMSC_AIC
## [1] -213.8917
##
## $MMSC_HQIC
## [1] -734.1071La comparación de los criterios MMSC entre el PVAR(1) y el PVAR(2) se realiza como con cualquier criterio de información: se elige la especificación con el menor valor. Conviene reportar los tres criterios (MMSC-AIC, MMSC-BIC y MMSC-HQIC), ya que –al igual que en el caso univariado– el criterio tipo BIC penaliza más la inclusión de parámetros y suele seleccionar órdenes menores.
Al igual que en el VAR de series de tiempo (Capítulo 5), la estabilidad del Panel VAR se verifica con los eigenvalores de la matriz compañera del sistema: el proceso es estable si todos los módulos son estrictamente menores que 1, es decir, si todas las raíces se encuentran dentro del círculo unitario. La Figura 7.3 muestra las tres raíces del modelo estimado, todas al interior del círculo, por lo que el PVAR(1) es estable:
## Eigenvalue stability condition:
##
## Eigenvalue Modulus
## 1 0.5370657+0.00000000i 0.53706574
## 2 0.0633651+0.06442426i 0.09036383
## 3 0.0633651-0.06442426i 0.09036383
##
## All the eigenvalues lie inside the unit circle.
## PVAR satisfies stability condition.plot(stab_ex1_dahlberg_data) +
coord_fixed() +
scale_x_continuous(name = "Parte real") +
scale_y_continuous(name = "Parte imaginaria") +
labs(title = "Raíces de la matriz compañera")
Figura 7.3: Raíces de la matriz compañera del Panel VAR(1) estimado. El proceso es estable porque todos los eigenvalores están dentro del círculo unitario.
Funciones de impulso-respuesta (comentadas; descomentar para ejecutar):
# ex1_dahlberg_data_oirf <- oirf(ex1_dahlberg_data, n.ahead = 8)
#
# ex1_dahlberg_data_girf <- girf(ex1_dahlberg_data, n.ahead = 8, ma_approx_steps= 8)
#
# ex1_dahlberg_data_bs <- bootstrap_irf(ex1_dahlberg_data, typeof_irf = c("GIRF"),
# n.ahead = 8,
# nof_Nstar_draws = 500,
# confidence.band = 0.95)
#
# plot(ex1_dahlberg_data_girf, ex1_dahlberg_data_bs)7.3 Cointegración en Panel
La extensión de los conceptos de cointegración al contexto panel sigue la misma lógica que en el caso univariado: si las series \(Y_{it}\) son \(I(1)\) pero existe una combinación lineal estacionaria, se dice que cointegran. El procedimiento más utilizado es el propuesto por Pedroni (1999, 2004), que plantea siete estadísticos de prueba bajo la hipótesis nula de no cointegración, permitiendo tanto coeficientes de cointegración homogéneos como heterogéneos entre individuos.
7.3.1 Prueba de Pedroni
Para cada individuo \(i = 1, \ldots, N\) se estima la regresión de cointegración de largo plazo: \[\begin{equation} Y_{it} = \alpha_i + \delta_i t + \boldsymbol{\beta}_i^{\prime} \mathbf{X}_{it} + e_{it} \tag{7.13} \end{equation}\]
donde \(\mathbf{X}_{it}\) es un vector \(m \times 1\) de regresores \(I(1)\), \(\alpha_i\) es un efecto fijo individual y \(\delta_i t\) es una tendencia determinista individual (que puede omitirse). La hipótesis nula es que los residuales \(\hat{e}_{it}\) tienen raíz unitaria, es decir, no existe cointegración para ningún individuo del panel: \[\begin{equation} H_0 : \hat{e}_{it} = \hat{e}_{i,t-1} + v_{it} \quad \forall \, i \end{equation}\]
Para construir los estadísticos se estima, para cada \(i\), la regresión auxiliar aumentada sobre los residuales: \[\begin{equation} \Delta \hat{e}_{it} = \hat{\rho}_i \hat{e}_{i,t-1} + \sum_{j=1}^{k_i} \hat{\psi}_{ij} \Delta \hat{e}_{i,t-j} + \hat{v}_{it} \tag{7.14} \end{equation}\]
A partir de los residuales \(\hat{v}_{it}\) y de las varianzas de largo plazo \(\hat{L}_{11i}^2\) (que corrigen la correlación serial), Pedroni (1999) construye siete estadísticos divididos en dos grupos:
Estadísticos de dimensión interna (within-dimension). Imponen un coeficiente autorregresivo \(\rho_i = \rho\) común al hacer la suma sobre individuos. Se componen del estadístico Panel-\(\nu\) (cociente de varianzas, diverge hacia la derecha bajo \(H_1\)), el Panel-\(\rho\) (análogo al estadístico \(\rho\) de Phillips-Perron), el Panel-\(t\) no paramétrico y el Panel-ADF paramétrico.
Estadísticos de dimensión entre grupos (between-dimension). Promedian estadísticos individuales sobre los \(N\) individuos, permitiendo que \(\rho_i\) difiera entre individuos bajo la alternativa. Se componen del Group-\(\rho\), el Group-\(t\) no paramétrico y el Group-ADF paramétrico.
Bajo la hipótesis nula de no cointegración, todos los estadísticos estandarizados convergen en distribución a una normal estándar: \[\begin{equation} \frac{\text{Estadístico} - \mu_{NT}}{\sigma_{NT}} \xrightarrow{d} N(0, 1) \end{equation}\]
donde \(\mu_{NT}\) y \(\sigma_{NT}\) son momentos que dependen del número de regresores \(m\) y se tabulan en Pedroni (1999). Los estadísticos Panel-\(\nu\) divergen hacia la derecha bajo la alternativa; los demás seis divergen hacia la izquierda, de modo que valores muy negativos llevan a rechazar \(H_0\).
La ventaja de los estadísticos between-dimension sobre los within-dimension es que permiten vectores de cointegración heterogéneos entre individuos bajo la hipótesis alternativa, lo que los hace más adecuados en paneles con coeficientes de largo plazo que varían entre individuos.
7.3.2 Prueba de Kao
Kao (1999) propone pruebas más restrictivas que imponen el mismo vector de cointegración para todos los individuos y no incluyen tendencia determinista individual. Sus estadísticos también se basan en los residuales de la regresión de cointegración y convergen a \(N(0,1)\) bajo la hipótesis nula de no cointegración.
7.3.3 Prueba de Westerlund
Westerlund (2007) propone cuatro estadísticos basados en la corrección de errores: en lugar de probar raíz unitaria en los residuales, se estima un modelo de corrección de errores para cada individuo y se prueba si el coeficiente de ajuste \(\alpha_i\) es significativamente negativo. Este enfoque tiene mayor potencia que las pruebas residuales cuando la dinámica de corto plazo varía entre individuos, y permite controlar la dependencia transversal mediante remuestreo (bootstrap).
7.3.4 Ejemplo: inversión y valor de la empresa en el panel de Grunfeld
Ilustremos la mecánica de las pruebas anteriores con el panel de Grunfeld (1958) que hemos utilizado a lo largo del capítulo (10 empresas, 1935–1954). La teoría de la inversión sugiere una relación de largo plazo entre la inversión de la empresa y su valor de mercado —en el espíritu de la \(q\) de Tobin (1969)—, de modo que postulamos la regresión de cointegración de la ecuación (7.13) con \(Y_{it} = ln(inv_{it})\) y \(X_{it} = ln(value_{it})\).
Paso 1: verificar que ambas series son \(I(1)\). Aplicamos la prueba de Levin-Lin-Chu en niveles y en diferencias:
LINV <- ts(data.frame(split(log(Grunfeld$inv), Grunfeld$firm)),
start = 1935, freq = 1)
LVAL <- ts(data.frame(split(log(Grunfeld$value), Grunfeld$firm)),
start = 1935, freq = 1)
purtest(LINV, test = "levinlin", exo = "intercept",
lags = "AIC", pmax = 4)## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: LINV
## z = -0.63214, p-value = 0.2636
## alternative hypothesis: stationarity## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: LVAL
## z = -0.74679, p-value = 0.2276
## alternative hypothesis: stationarity## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: diff(LINV)
## z = -8.5076, p-value < 2.2e-16
## alternative hypothesis: stationarity## Warning in selectT(l, theTs): the time series is short##
## Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts)
##
## data: diff(LVAL)
## z = -11.944, p-value < 2.2e-16
## alternative hypothesis: stationarityEn niveles no se rechaza la hipótesis nula de raíz unitaria para ninguna de las dos variables, mientras que en diferencias se rechaza contundentemente: ambas series son \(I(1)\), el escenario donde tiene sentido preguntarse por cointegración. (En paneles tan cortos como éste, \(T = 20\), distintas pruebas pueden arrojar conclusiones encontradas —la prueba IPS con tendencia, por ejemplo, rechaza la raíz unitaria en niveles—, por lo que en la práctica conviene contrastar más de una especificación antes de concluir.)
Paso 2: estimar la relación de largo plazo. Estimamos la regresión de cointegración con efectos fijos por empresa:
eq_lp <- plm(log(inv) ~ log(value), data = Grunfeld,
index = c("firm", "year"), model = "within")
coef(summary(eq_lp))## Estimate Std. Error t-value Pr(>|t|)
## log(value) 0.9713448 0.0959524 10.12319 1.607934e-19La elasticidad estimada de largo plazo es cercana a uno: en el largo plazo, la inversión crece proporcionalmente con el valor de mercado de la empresa.
Paso 3: prueba residual (idea de Kao). Si existe cointegración, los residuales \(\hat{e}_{it}\) de la regresión anterior deben ser estacionarios. Siguiendo la lógica de Kao (1999), estimamos la regresión ADF agrupada sobre los residuales, imponiendo el mismo coeficiente \(\rho\) para todas las empresas:
## Estimate Std. Error t value Pr(>|t|)
## lag(ehat) -0.2319281 0.04987178 -4.650488 0.000006207847El coeficiente de \(\hat{e}_{i,t-1}\) es negativo y su estadístico \(t\) es grande en valor absoluto. Cabe advertir que, al calcularse sobre residuales estimados, este estadístico no sigue la distribución \(t\) estándar: la prueba formal de Kao le aplica correcciones que lo llevan a una \(N(0,1)\). Dichas correcciones están implementadas en software especializado (por ejemplo, xtcointtest en Stata o las pruebas de panel de EViews); en R, las implementaciones disponibles se encuentran en paquetes archivados, por lo que aquí mostramos la mecánica de la prueba de forma transparente. Con un estadístico de esta magnitud, la evidencia contra la hipótesis nula de no cointegración es clara.
Paso 4: ecuación de corrección de error (idea de Westerlund). Como verificación complementaria, estimamos el mecanismo de corrección de error en panel: si hay cointegración, las desviaciones del equilibrio de largo plazo (\(\hat{e}_{i,t-1}\)) deben corregirse, es decir, su coeficiente \(\alpha\) debe ser negativo y significativo:
pGrunfeld <- pdata.frame(Grunfeld, index = c("firm", "year"))
pGrunfeld$ehat <- ehat
mce <- plm(diff(log(inv)) ~ lag(ehat) + lag(diff(log(inv))) +
diff(log(value)),
data = pGrunfeld, model = "within")
coef(summary(mce))## Estimate Std. Error t-value Pr(>|t|)
## lag(ehat) -0.2631642 0.05031209 -5.230636 5.005279e-07
## lag(diff(log(inv))) 0.1455486 0.06613555 2.200762 2.912503e-02
## diff(log(value)) 0.6766472 0.09021746 7.500180 3.582884e-12El coeficiente de ajuste estimado es \(\hat{\alpha} \approx -0.26\) y altamente significativo: cada año se corrige alrededor de una cuarta parte de la desviación respecto del equilibrio de largo plazo, lo que implica una vida media del desequilibrio de poco más de dos años. Ambos enfoques —el residual y el de corrección de error— apuntan a la misma conclusión: la inversión y el valor de la empresa cointegran en este panel.
7.4 Modelos No Lineales
7.4.1 Modelos de cambio de régimen
La exposición de los modelos de umbral de esta sección –SETAR, STAR y su extensión a varios regímenes– sigue a Franses y van Dijk (2000, cap. 3).
En años recientes, los modelos de serie de tiempo han sido incorporados en análisis de la existencia de diferentes estados que son generados por procesos estocásticos subyacentes. En esta sección del curso revisaremos algunos modelos de cambio de régimen. Restringimos nuestra revisión a modelos que asuman que la dinámica de las series puede ser descrita por modelos del tipo AR y dejamos fuera procesos del tipo MA.
En general distinguimos que existen dos tipos de modelos: 1. Modelos caracterizados por una variable observable, por lo que tenemos certeza de los regímenes en cada uno de los momentos.
- Modelos en los que el régimen no puede ser observado por una variable, pero sí conocemos el proceso estocástico subyacente.
7.4.2 Regímenes determinados por información observable
En estos casos asumimos que el régimen ocurre en un momento \(t\) y puede ser determinado por una variable observable. Este modelo es conocido como el modelo autorregresivo con umbral (TAR, Threshold Autoregressive model). En este caso también diremos que cuando el régimen está determinado por la información de la misma serie será llamado Self-Exciting TAR (SETAR).
Veamos un caso particular. Supongamos que existe un umbral, \(c\), para el régimen que está determinado por \(q_t = y_{t-1}\) y que el estado de la naturaleza nos permite establecer dos estados o regímenes: \[\begin{equation} y_t = \begin{cases} \phi_{01} + \phi_{11} y_{t-1} + \varepsilon_t \text{ si } y_{t-1} \leq c \\ \phi_{02} + \phi_{12} y_{t-1} + \varepsilon_t \text{ si } y_{t-1} > c \end{cases} \end{equation}\]
Donde asumiremos que \(\varepsilon_t\) es i.i.d y que cumple con: \[\begin{equation*} \mathbb{E}[\varepsilon_t | \Omega_{t-1}] = 0 \end{equation*}\]
Donde \(\Omega_{t-1} = \{ y_{t-1}, y_{t-2}, \ldots \}\). Existe una variante de este modelo que suaviza la transición entre regímenes conocido como Smooth Transition AR (STAR) y puede ser especificada en su modalidad de dos regímenes como: \[\begin{equation} y_t = (\phi_{01} + \phi_{11} y_{t-1}) (1 - G(y_{t-1}; \gamma, c)) + (\phi_{02} + \phi_{12} y_{t-1}) G(y_{t-1}; \gamma, c) + \varepsilon_t \end{equation}\]
Donde \(G(y_{t-1}; \gamma, c)\) es una función continua que suaviza la transición entre regímenes, propuesta en este contexto por Teräsvirta (1994). La práctica común es suponer que tiene una forma logística: \[\begin{equation} G(y_{t-1}; \gamma, c) = \frac{1}{1 + e^{-\gamma (y_{t-1} - c)}} \end{equation}\]
Es posible hacer extensiones de lo anterior a modelos de orden superior; dando como resultado: \[\begin{equation} y_t = \begin{cases} \phi_{01} + \phi_{11} y_{t-1} + \phi_{21} y_{t-2} + \ldots + \phi_{p_1 1} y_{t-p_1} + \varepsilon_t \text{ si } y_{t-1} \leq c \\ \phi_{02} + \phi_{12} y_{t-1} + \phi_{22} y_{t-2} + \ldots + \phi_{p_2 2} y_{t-p_2} + \varepsilon_t \text{ si } y_{t-1} > c \end{cases} \end{equation}\]
En el segundo caso: \[\begin{eqnarray*} y_t & = & (\phi_{01} + \phi_{11} y_{t-1} + \phi_{21} y_{t-2} + \ldots + \phi_{p_1 1} y_{t-p_1}) (1 - G(y_{t-1}; \gamma, c)) \\ & & + (\phi_{02} + \phi_{12} y_{t-1} + \phi_{22} y_{t-2} + \ldots + \phi_{p_2 2} y_{t-p_2}) G(y_{t-1}; \gamma, c) \\ & & + \varepsilon_t \end{eqnarray*}\]
De igual forma que en el caso de los modelos ARIMA y VAR, el número de rezagos se determina con criterios de información. Para un modelo SETAR de dos regímenes, Tong (1990) propone sumar los criterios de Akaike de los modelos \(AR\) de cada régimen: \[\begin{equation} AIC(p_1, p_2) = n_1 ln(\hat{\sigma}^2_1) + n_2 ln(\hat{\sigma}^2_2) + 2(p_1 + 1) + 2(p_2 + 1) \end{equation}\]
Los modelos SETAR y STAR generan procesos estacionarios siempre que cumplan ciertas condiciones. En este libro nos enfocaremos únicamente en el modelo SETAR. Las condiciones exactas se deben a Chan, Petruccelli, Tong y Woolford (1985) y las reproducimos en la formulación de Franses y van Dijk (2000, 79): el modelo SETAR de primer orden es estacionario si y sólo si se cumple alguna de las siguientes condiciones: 1. \(\phi_{11} < 1\), \(\phi_{12} < 1\), \(\phi_{11} \cdot \phi_{12} < 1\)
\(\phi_{11} = 1\), \(\phi_{12} < 1\), \(\phi_{01} > 0\)
\(\phi_{11} < 1\), \(\phi_{12} = 1\), \(\phi_{02} < 0\)
\(\phi_{11} = 1\), \(\phi_{12} = 1\), \(\phi_{02} < 0 < \phi_{01}\)
\(\phi_{11} \cdot \phi_{12} = 1\), \(\phi_{11} < 0\), \(\phi_{02} + \phi_{12} \cdot \phi_{01} > 0\)
Finalmente, en ocasiones podemos estar interesados en modelos con más de dos regímenes bajo una especificación SETAR o STAR. En el caso de un modelo SETAR, \(m\) regímenes quedan delimitados por \(m + 1\) puntos de corte \(c_0, c_1, \ldots, c_m\), de los cuales sólo \(m - 1\) son umbrales que deben estimarse, pues los dos extremos son infinitos: \[\begin{equation*} -\infty = c_0 < c_1 < \ldots < c_{m-1} < c_m = \infty \end{equation*}\]
Así, tendríamos ecuaciones: \[\begin{equation} y_t = \phi_{0j} + \phi_{1j} y_{t-1} + \varepsilon_t \text{ si } c_{j-1} < y_{t-1} < c_j \end{equation}\]
Para \(j = 1, 2, \ldots, m\). De forma similar, podemos recomponer el modelo STAR.
7.4.3 Regímenes determinados por variables no observables
Este tipo de modelos asume que el régimen ocurre en el momento \(t\) y que no puede ser observado, ya que este es determinado por un proceso no observable, el cual denotamos como \(s_t\). En el caso de dos regímenes, \(s_t\) puede ser asumido como que toma 2 valores: 1 y 2, por ejemplo. Supongamos que el proceso subyacente tiene una forma del tipo AR(1) dado por: \[\begin{equation} y_t = \begin{cases} \phi_{01} + \phi_{11} y_{t-1} + \varepsilon_t \text{ si } s_t = 1 \\ \phi_{02} + \phi_{12} y_{t-1} + \varepsilon_t \text{ si } s_t = 2 \end{cases} \tag{7.15} \end{equation}\]
O en un formato más corto de notación: \[\begin{equation} y_t = \phi_{0 s_t} + \phi_{1 s_t} y_{t-1} + \varepsilon_t \end{equation}\]
Para complementar el modelo, las propiedades del proceso \(s_t\) necesitan ser especificadas. El modelo más popular dentro de esta familia es el propuesto por Hamilton (1989), el cual es conocido como Markov Switching Model (MSM), en el cual el proceso \(s_t\) se asume como un proceso de Markov de primer orden. Esto implica que el régimen actual \(s_t\) sólo depende del período \(s_{t-1}\).
Así, el modelo es completado mediante la definición de las probabilidades de transición para moverse de un estado a otro: \[\begin{eqnarray*} \mathbb{P}(s_t = 1 | s_{t-1} = 1) & = & p_{11} \\ \mathbb{P}(s_t = 2 | s_{t-1} = 1) & = & p_{12} \\ \mathbb{P}(s_t = 1 | s_{t-1} = 2) & = & p_{21} \\ \mathbb{P}(s_t = 2 | s_{t-1} = 2) & = & p_{22} \end{eqnarray*}\]
Así, \(p_{ij}\) es igual a la probabilidad de que la cadena de Markov pase del estado \(i\) en el momento \(t-1\) al estado \(j\) en el tiempo \(t\). En todo caso asumiremos que \(p_{ij} > 0\) y que: \[\begin{eqnarray*} p_{11} + p_{12} = 1 \\ p_{21} + p_{22} = 1 \end{eqnarray*}\]
Otro tipo de probabilidades a analizar son las probabilidades incondicionales de \(\mathbb{P}(s_t = i)\), \(i = 1, 2\). Usando la teoría ergódica de las cadenas de Markov, estas probabilidades están dadas por: \[\begin{eqnarray*} \mathbb{P}(s_t = 1) & = & \frac{1 - p_{22}}{2 - p_{11} - p_{22}} \\ \mathbb{P}(s_t = 2) & = & \frac{1 - p_{11}}{2 - p_{11} - p_{22}} \end{eqnarray*}\]
Un caso más general es el de múltiples regímenes en el cual \(s_t\) puede tomar cualquier valor \(m > 2\), \(m \in \mathbb{N}\). Este modelo se puede escribir como: \[\begin{equation} y_t = \phi_{0j} + \phi_{1j} y_{t-1} + \varepsilon_t \text{ si } s_t = j \text{, para } j = 1, 2, \ldots, m \end{equation}\]
Donde las probabilidades de transición estarán dadas por: \[\begin{equation} p_{ij} = \mathbb{P}(s_t = j | s_{t-1} = i) \text{ para } i , j = 1, 2, \ldots, m \end{equation}\]
Donde la ecuación anterior satisface que \(p_{ij} > 0\), \(\forall i, j = 1, 2, \ldots, m\) y que: \[\begin{equation*} \sum_{j=1}^m p_{ij} = 1 \text{ para } i = 1, 2, \ldots, m \end{equation*}\]
Finalmente, plantearemos el procedimiento empírico seguido para la estimación de este tipo de modelos: 1. Estimar un proceso AR(p)
Probar una hipótesis nula de no linealidad de acuerdo con alguno de los modelos SETAR, STAR o MSM
Estimar los parámetros
Evaluar los resultados del modelo
Ajustar, en su caso, la estimación o modelo
Pronosticar o realizar análisis impulso-respuesta
7.5 Ejemplo: Modelos de Cambio de Régimen (TAR)
Dependencias y configuración:
Datos: tasas mensuales de muertes por influenza y neumonía en Estados
Unidos (11 años). La serie es la que utilizan Shumway y Stoffer
(2017) en el paquete astsa:
## X.8113721
## 1 0.4458291
## 2 0.3415985
## 3 0.2774243
## 4 0.2484958
## 5 0.2525427
## 6 0.2466902Conversión a series de tiempo:
La Figura 7.4 muestra la serie en niveles y en primeras diferencias. En niveles se aprecia con claridad el patrón estacional de los brotes de influenza; en diferencias destaca la asimetría que motiva un modelo de umbral: los incrementos son abruptos –el inicio del brote– y las caídas, aunque también rápidas, siguen una dinámica distinta.
par(mfrow = c(2, 1))
plot(flu, type = "b", col = "darkred", ylab = "Muertes por 10,000 habitantes",
xlab = "Meses",
main = "Tasas mensuales de muertes por gripe en Estados Unidos")
plot(D_flu, type = "b", col = "darkblue", ylab = "Diferencia mensual",
xlab = "Meses",
main = "Diferencias de las tasas mensuales de muertes por gripe en EE.UU.")
Figura 7.4: Tasa mensual de muertes por influenza y neumonía en Estados Unidos (panel superior) y su primera diferencia (panel inferior)
El paquete tsDyn simplifica la estimación en pocos pasos. Estimamos
primero un SETAR con cuatro rezagos (\(m = 4\)), variable de transición
\(y_{t-1}\) (thDelay = 0) y un umbral fijado en \(c = 0.05\):
## Warning:
## With the threshold you gave (0.05) there is a regime with less than trim=15% observations (86.51%, 13.49%, )## Warning: Possible unit root in the high regime. Roots are: 0.6182
## 0.6244 0.6244 0.6182##
## Non linear autoregressive model
##
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
## const.L phiL.1 phiL.2 phiL.3 phiL.4
## 0.004432028 0.501179574 -0.206693004 0.120140800 -0.122254140
##
## High regime:
## const.H phiH.1 phiH.2 phiH.3 phiH.4
## 0.4079353 -0.7483325 -1.0323129 -2.0450407 -6.7117769
##
## Threshold:
## -Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
## -Value: 0.05 (fixed)
## Proportion of points in low regime: 86.51% High regime: 13.49%
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1319314 -0.0217073 0.0015801 0.0196093 0.2674245
##
## Fit:
## residuals variance = 0.002166, AIC = -778, MAPE = 412.4%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const.L 0.0044320 0.0051791 0.8558 0.3938357
## phiL.1 0.5011796 0.0833821 6.0106 2.043e-08 ***
## phiL.2 -0.2066930 0.0608839 -3.3949 0.0009313 ***
## phiL.3 0.1201408 0.0576505 2.0839 0.0392868 *
## phiL.4 -0.1222541 0.0522363 -2.3404 0.0209137 *
## const.H 0.4079353 0.0314121 12.9866 < 2.2e-16 ***
## phiH.1 -0.7483325 0.1118329 -6.6915 7.455e-10 ***
## phiH.2 -1.0323129 0.1420203 -7.2688 4.018e-11 ***
## phiH.3 -2.0450407 0.7055163 -2.8986 0.0044573 **
## phiH.4 -6.7117769 0.8367945 -8.0208 7.911e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)+ (0) X(t-3)
##
## Value: 0.05 (fixed)# ask = FALSE evita que tsDyn pida un Enter entre los paneles del
# diagnóstico cuando el libro se compila desde la consola.
plot(D_flu_tar4_05, ask = FALSE)
Figura 7.5: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.6: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.7: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.8: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.9: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.10: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.11: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.12: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.13: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Figura 7.14: Diagnóstico del modelo SETAR con umbral fijo \(c = 0.05\): serie observada y ajustada, y distribución de las observaciones entre regímenes
Si no se especifica el umbral (th), setar lo busca automáticamente
sobre una grilla de valores candidatos, eligiendo el que minimiza la suma
de residuales al cuadrado del modelo:
## Warning: Possible unit root in the high regime. Roots are: 0.6316
## 0.7222 0.7222 0.6316##
## Non linear autoregressive model
##
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
## const.L phiL.1 phiL.2 phiL.3
## -0.00006287604 0.44640751264 -0.23158878472 0.10701408961
## phiL.4
## -0.14406085874
##
## High regime:
## const.H phiH.1 phiH.2 phiH.3 phiH.4
## 0.3486150 -0.5903335 -1.0318488 -2.4053812 -4.8052422
##
## Threshold:
## -Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
## -Value: 0.03798
## Proportion of points in low regime: 84.92% High regime: 15.08%
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.277686 -0.020554 0.005116 0.020762 0.119627
##
## Fit:
## residuals variance = 0.002405, AIC = -762, MAPE = 372.5%
##
## Coefficient(s):
##
## Estimate Std. Error t value Pr(>|t|)
## const.L -0.000062876 0.005572290 -0.0113 0.9910158
## phiL.1 0.446407513 0.089159422 5.0068 1.921e-06 ***
## phiL.2 -0.231588785 0.064328877 -3.6001 0.0004636 ***
## phiL.3 0.107014090 0.061034049 1.7534 0.0820953 .
## phiL.4 -0.144060859 0.055198346 -2.6099 0.0102109 *
## const.H 0.348615027 0.027460007 12.6954 < 2.2e-16 ***
## phiH.1 -0.590333531 0.109881869 -5.3724 3.869e-07 ***
## phiH.2 -1.031848767 0.149570541 -6.8987 2.641e-10 ***
## phiH.3 -2.405381238 0.691344783 -3.4793 0.0007014 ***
## phiH.4 -4.805242210 0.828241480 -5.8017 5.451e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold
## Variable: Z(t) = + (1) X(t) + (0) X(t-1)+ (0) X(t-2)+ (0) X(t-3)
##
## Value: 0.03798El umbral estimado es \(\hat{c} =\) 0.038, muy cercano al valor de 0.05 que habíamos impuesto, lo que indica que la partición en regímenes es robusta. La interpretación económica –o, en este caso, epidemiológica– es directa: cuando el cambio mensual en la tasa de mortalidad del mes anterior supera el umbral, la serie se encuentra en fase de brote y su dinámica es distinta de la del régimen de baja mortalidad.
Figura 7.15: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.16: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.17: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.18: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.19: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.20: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.21: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.22: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.23: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
Figura 7.24: Diagnóstico del modelo SETAR con umbral estimado por búsqueda en grilla
7.6 Resumen del capítulo
Este capítulo amplió el análisis de series de tiempo a dos dimensiones adicionales: la dimensión transversal (datos panel) y la no linealidad (modelos de cambio de régimen).
Datos panel. Los datos panel combinan series de tiempo para múltiples individuos o unidades (\(N\) individuos, \(T\) periodos). El modelo básico incluye efectos fijos \(\alpha_i\) que capturan heterogeneidad no observada constante en el tiempo. Las pruebas de raíz unitaria en panel (Levin et al. 2002; Im et al. 2003; Maddala and Wu 1999; Hadri 2000) amplían el poder estadístico al explotar la dimensión transversal; plantean hipótesis nulas de no estacionariedad con hipótesis alternativas homogéneas o heterogéneas según la prueba.
Diferencia en diferencias (DiD). El estimador DiD identifica el efecto causal de un tratamiento comparando la variación en el tiempo del grupo tratado con la del grupo de control. Su validez descansa en el supuesto de tendencias paralelas: en ausencia del tratamiento, ambos grupos habrían seguido trayectorias paralelas.
Panel VAR. El Panel VAR (PVAR) combina las ventajas del VAR clásico con la dimensión transversal del panel. Dado que el estimador within es sesgado en paneles dinámicos con \(T\) pequeño (sesgo de Nickell), se estima mediante GMM con instrumentos basados en los rezagos de las propias variables, controlando por efectos fijos mediante transformación de primeras diferencias o desviaciones ortogonales; la validez de los instrumentos se evalúa con la prueba de Hansen y las pruebas de autocorrelación de orden 1 y 2. La selección del número de rezagos sigue el criterio de Andrews y Lu (2001).
Cointegración en panel. Cuando las series del panel son \(I(1)\), las pruebas de Pedroni (1999), Kao (1999) y Westerlund (2007) permiten contrastar la existencia de relaciones de largo plazo. Las pruebas residuales verifican la estacionariedad de los residuales de la regresión de cointegración, mientras que el enfoque de corrección de errores verifica que las desviaciones del equilibrio se corrijan (\(\alpha < 0\)); ambos enfoques se ilustraron con el panel de Grunfeld.
Modelos no lineales. Cuando la dinámica de una serie cambia de forma discreta según el estado del sistema, los modelos de cambio de régimen son más apropiados que los modelos lineales.
TAR (Threshold Autoregressive): el régimen depende de si una variable observable supera un umbral \(c\). En la variante SETAR, el umbral lo determina la propia serie rezagada.
STAR (Smooth Transition AR): la transición entre regímenes es gradual y está gobernada por una función logística paramétrica \(G(y_{t-1};\gamma,c)\).
Markov Switching (MSM): el régimen \(s_t\) es una variable latente que sigue una cadena de Markov de primer orden. Las probabilidades de transición \(p_{ij} = \mathbb{P}(s_t = j \mid s_{t-1} = i)\) son estimadas junto con los parámetros del modelo AR condicional al régimen.
En todos los casos, el procedimiento empírico incluye: estimar un AR(p) lineal, probar la hipótesis de linealidad, estimar el modelo no lineal seleccionado, evaluar resultados y, finalmente, generar pronósticos o funciones de impulso-respuesta.
7.7 Ejercicios
Los siguientes ejercicios combinan preguntas teóricas con aplicaciones en R. Los marcados como (Teórico) se resuelven analíticamente; los marcados como (Computacional) utilizan los paquetes plm y tsDyn empleados en este capítulo, así como las bases de datos incluidas en plm.
(Teórico) Considere el modelo de efectos fijos \(Y_{it} = \alpha_i + \beta X_{it} + U_{it}\), con \(i = 1, \ldots, N\) y \(t = 1, \ldots, T\).
- Demuestre que la transformación within —restar a cada variable su media por individuo, por ejemplo \(Y_{it} - \bar{Y}_i\) con \(\bar{Y}_i = \frac{1}{T}\sum_t Y_{it}\)— elimina el efecto fijo \(\alpha_i\).
- Suponga ahora un panel dinámico: \(Y_{it} = \alpha_i + \rho Y_{i,t-1} + U_{it}\). Explique por qué, tras la transformación within, el regresor transformado \(Y_{i,t-1} - \bar{Y}_{i,-1}\) queda mecánicamente correlacionado con el error transformado, lo que genera el sesgo de Nickell.
- ¿Qué ocurre con ese sesgo cuando \(T \to \infty\)? ¿Por qué es especialmente preocupante en paneles macroeconómicos con \(T\) pequeño?
- Describa la solución de Arellano y Bond: qué transformación se aplica, qué se utiliza como instrumento y qué condiciones deben cumplir los instrumentos para ser válidos.
(Teórico) Sobre las pruebas de raíces unitarias en panel estudiadas en este capítulo:
- Escriba la hipótesis nula y la alternativa de la prueba de Levin, Lin y Chu (caso homogéneo, \(H_1^A\)) y de la prueba de Im, Pesaran y Shin (caso heterogéneo, \(H_1^B\)), en términos de los coeficientes \(\rho_i\) de la ecuación (7.3).
- Explique con precisión qué se puede concluir al rechazar \(H_0\) en la prueba IPS. ¿Es correcto afirmar que “todas las series del panel son estacionarias”?
- ¿Cuál de las dos pruebas resulta más apropiada cuando se sospecha que la velocidad de reversión a la media difiere entre individuos? Justifique.
- Explique cuál es la ventaja de las pruebas en panel frente a aplicar una prueba ADF individual a cada una de las \(N\) series.
(Teórico) Considere el modelo de diferencia en diferencias \(Y_{it} = \beta_0 + \beta_1 D_i + \beta_2 P_t + \beta_3 (D_i \times P_t) + U_{it}\), donde \(D_i\) indica pertenencia al grupo de tratamiento y \(P_t\) el periodo posterior a la intervención.
- Exprese las cuatro medias condicionales \(\mathbb{E}[Y_{it} | D_i, P_t]\) en términos de los parámetros del modelo.
- Demuestre que \(\beta_3\) es igual a la doble diferencia de medias: la variación del grupo tratado antes-después menos la variación del grupo de control antes-después.
- Enuncie el supuesto de tendencias paralelas y explique por qué es una condición de identificación no verificable directamente. ¿Qué evidencia empírica previa a la intervención suele presentarse para hacerlo plausible?
(Computacional) Utilice la base
Grunfelddel paqueteplm(inversión, valor de mercado y capital de 10 empresas estadounidenses, 1935-1954) para estimar el modelo \(inv_{it} = \alpha_i + \beta_1 value_{it} + \beta_2 capital_{it} + U_{it}\).- Estime el modelo por pooling, por efectos fijos (within) y por efectos aleatorios con
plm(), y compare los coeficientes estimados. - Aplique la prueba de Hausman con
phtest()y concluya qué especificación es preferible. - Interprete económicamente el coeficiente de \(value_{it}\) en la especificación elegida.
- Extraiga los efectos fijos estimados con
plm::fixef()–conviene usar el prefijo del paquete, porque otros paquetes también exportan una funciónfixef()– y comente qué tan heterogéneas son las empresas del panel.
- Estime el modelo por pooling, por efectos fijos (within) y por efectos aleatorios con
(Computacional) Utilice la base
Producdel paqueteplm(producto estatal bruto y otras variables para 48 estados de Estados Unidos, 1970-1986, tomada de Munnell (1990)) y analice el orden de integración del logaritmo del producto estatal bruto, \(ln(gsp_{it})\).- Construya la matriz de series por estado, como se hizo en este capítulo con la base Grunfeld, y aplique la prueba de Levin, Lin y Chu con
purtest(..., test = "levinlin"), con especificaciones de intercepto y de tendencia. - Aplique la prueba de Im, Pesaran y Shin (
test = "ips") con las mismas especificaciones. - Repita ambas pruebas sobre las primeras diferencias de las series.
- Concluya el orden de integración del panel y comente si las conclusiones de ambas pruebas coinciden. En caso contrario, explique la discrepancia a la luz de las hipótesis alternativas de cada prueba.
- Construya la matriz de series por estado, como se hizo en este capítulo con la base Grunfeld, y aplique la prueba de Levin, Lin y Chu con
(Computacional) Continúe con la serie de diferencias del logaritmo de los casos de influenza (
D_flu) empleada en el ejemplo SETAR de este capítulo.- Aplique la prueba de linealidad de Hansen con
setarTest()del paquetetsDyny concluya si se justifica un modelo de umbral frente a un AR lineal. - Estime un modelo SETAR con
m = 2y umbral buscado automáticamente, y compárelo mediante el criterio AIC con el modelo SETAR dem = 4estimado en el capítulo y con un AR lineal del mismo orden (funciónlinear()detsDyn). - Interprete los regímenes del modelo seleccionado: ¿cuántas observaciones caen en cada régimen y qué dinámica exhibe cada uno?
- Explique por qué un modelo lineal resulta insuficiente para una serie epidemiológica como esta, en la que los brotes generan episodios de crecimiento explosivo seguidos de reversiones rápidas.
- Aplique la prueba de linealidad de Hansen con