Errores autocorrelacionados

2023

Errores autocorrelacionados

Considere el caso en el que el supuesto (S4) no se cumple, pero el supuesto (S3) (y los otros supuestos restantes (S1), (S2), (S5) y (S6)) siguen siendo válidos. Entonces podemos escribir el siguiente modelo como:

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \quad \mathbb{E} \left( \boldsymbol{\varepsilon} | \mathbf{X}\right) = \boldsymbol{0},\quad \mathbb{V}{\rm ar}\left( \boldsymbol{\varepsilon} | \mathbf{X} \right) = \mathbb{E} \left( \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top \right)= \mathbf{\Sigma} \neq \sigma_\epsilon^2 \mathbf{I} \]

El caso en el que la matriz de varianza-covarianza del error no sea diagonal, sino con elementos diagonales iguales, expresados como:

Errores autocorrelacionados

\[ \mathbf{\Sigma} = \begin{bmatrix} \sigma^2 & \sigma_{1,2} & \sigma_{1,3} & ... & \sigma_{1,N}\\ \sigma_{2,1} & \sigma^2 & \sigma_{2,3} & ... & \sigma_{2,N}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \sigma_{N,1} & \sigma_{N,2} & \sigma_{N3,} & ... & \sigma^2 \end{bmatrix}, \quad \sigma_{i,j} = \sigma_{j,i} = \mathbb{C}{\rm ov}(\epsilon_i, \epsilon_j) \neq 0, \quad i,j = 1,...,N \]

se conoce como autocorrelación (o correlación en serie).

Definición de Autocorrelación

Autocorrelación significa la correlación de los valores de una misma variable ordenados en el tiempo (con datos de series temporales) o en el espacio (con datos espaciales). En el contexto de la regresión, decimos que hay autocorrelación en los términos de error cuando:

\[Cov(\varepsilon_{t}, \varepsilon_{t+s}) = E(\varepsilon_{t}, \varepsilon_{t-s}) \neq 0\]

Significa afirmar que los valores pasados, presentes o futuros de los errores están correlacionados. Perciba que las notaciones $t-s$ y $t+s$ son equivalentes en el conceptos de covarianza. Sin embargo en el análisis de regresión es más adecuado hablar de la correlación de los errores presentes con sus valores pasados.

Definición de Autocorrelación

El tipo mas común de autocorrelación es aquella dada por un proceso autoregresivo de primer orden, AR(1):

\[\varepsilon_{t} = \rho\varepsilon_{t-1} + u_{t}\] $\rho$ es el coeficiente de autocorrelación de los errores $(-1<\rho<1)$ y $u_{t}$ son errores IID.

Causas de la autocorrelación

Cabe destacar que la autocorrelación suele estar presente en los datos de series de tiempo. Para datos transversales, los errores pueden estar correlacionados en términos de distancia social o geográfica. Por ejemplo, la distancia entre ciudades, pueblos, barrios, etc.

Causas de la autocorrelación

Inercia: las series temporales económicas suelen presentar ciclos, es decir, períodos de crecimiento o decaimiento. Cuando este comportamiento se refleja en los factores no observados ($\varepsilon_{t}$), es común que los cambios en la tendencia ocurrieron lentamente. \[Y_{t} = \beta_{0} + \beta_{1}X_{t} + \varepsilon_{t} {\; \rightarrow\;} \varepsilon_{t} = \rho\varepsilon_{t-1} + u_{t}\]
Fallas de especificación: la autocorrelación puede deberse a la ausencia de un regresor importante o transformación de las variables existentes. Los errores expresan así un patrón sistemático debido a la ausencia de dicha información.

\[Y_{t} = \beta_{0} + \beta_{1}X_{t} + \varepsilon_{t} {\; \rightarrow\;} Y_{t} = \beta_{0} + \beta_{1}X_{t} + \beta_{2}X^{2}_{t} + \varepsilon_{t}\]

Causas de la autocorrelación

Desfases: las decisiones económicas en un período $t$ dependen, muchas veces, de información desfasada del período ${t-1}$. Desconsiderar este tipo de relación sujeta los errores a la correlación serial.

\[\begin{eqnarray*} Y_{t} = \beta_{0} + \beta_{1}X_{t} + \varepsilon_{t} {\; \rightarrow\;} Y_{t} = \beta_{0} + \beta_{1}X_{t} + \beta_{2}X_{t-1} + \varepsilon_{t}\\ o Y_{t} = \beta_{0} + \delta X_{t} + \beta_{2}X_{t-1} + \varepsilon_{t} \end{eqnarray*}\]

Consecuencias de la autocorrelación

Al igual que antes con la heterocedasticidad:

Los estimadores de MCO permanecen insesgados y consistentes;
Los estimadores de MCO ya no son eficientes;
El estimador de varianza de las estimaciones de MCO es sesgado e inconsistente;
Estadística-$t$ de las estimaciones de MCO no son válidas.

Consecuencias de la autocorrelación

Identificación

Tratamiento de la autocorrelación

Asumamos que $\mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}$ es el verdadero modelo.
Realice un análisis visual
Cree un modelo para los residuales de MCO. Por ejemplo, la regresión auxiliar residual en la prueba de Breusch-Godfrey incluye las variables explicativas del modelo original, así como los residuos de cierto orden:

\[ \widehat{\epsilon}_i = \alpha_0 + \alpha_1 X_{1, i} + ... + \alpha_k X_{k,i} + \rho_1 \widehat{\epsilon}_{i - 1} + \rho_2 \widehat{\epsilon}_{i - 2} + ... + \rho_p \widehat{\epsilon}_{i - p} + u_t \]

Pruebe la hipótesis nula de que los residuos están autocorrelacionados utilizando la regresión auxiliar residual de MCO: \[ H_0: \rho_1 = \rho_2 = ... = \rho_p = 0 \]

Tratamiento de la autocorrelación

Si no rechazamos la hipótesis nula, podemos utilizar los estimadores de MCO habituales.
Si rechazamos la hipótesis nula, podemos proceder de dos maneras:

Utilice los estimadores de MCO, pero corrija sus estimadores de varianza (es decir, hágalos consistentes);
En lugar de MCO, use FGLS (y sus variaciones).
Intente especificar un modelo diferente

Ejemplo con R

Simularemos el siguiente modelo:

\[ \begin{cases} Y_{i} &= \beta_0 + \beta_1 X_{1,i} + \beta_2 X_{2,i} + \epsilon_i \\ \epsilon_i &= \rho \epsilon_{i-1} + u_i,\ |\rho| < 1,\ u_i \sim \mathcal{N}(0, \sigma^2),\quad \epsilon_0 = 0 \end{cases} \]

set.seed(123)

N        <- 200
beta_vec <- c(10, 5, -3)
rho      <- 0.8

x1 <- seq(from = 0, to = 5, length.out = N)
x2 <- sample(seq(from = 3, to = 17, length.out = 80), size = N, replace = TRUE)

# residuos autocorrelacionados:
ee <- rnorm(mean = 0, sd = 3, n = N)
for(i in 2:N){
  ee[i] <- rho * ee[i-1] + ee[i]
}

x_mat    <- cbind(1, x1, x2)
y        <- x_mat %*% beta_vec + ee
data_mat <- data.frame(y, x1, x2)

Prueba para la autocorrelación

Podemos examinar la presencia de autocorrelación a partir de los gráficos de residuos, así como realizar una serie de pruebas formales. Comenzaremos estimando nuestro modelo a través de MCO, como lo haríamos habitualmente.

mdl_1 <- lm(y ~ x1 + x2, data = data_mat)
print(round(coef(summary(mdl_1)), 5))

##             Estimate Std. Error   t value Pr(>|t|)
## (Intercept) 10.61601    0.97983  10.83457        0
## x1           4.89018    0.20457  23.90497        0
## x2          -2.93686    0.07545 -38.92241        0

Correlograma residual

Comenzamos calculando los residuos del modelo de MCO:

\[ \widehat{\boldsymbol{\varepsilon}} = \mathbf{Y} - \mathbf{X} \widehat{\boldsymbol{\beta}}_{\rm MCO} \]

Usaremos estos residuales como estimaciones para el término de error verdadero (no observado) $\boldsymbol{\varepsilon}$ y examinaremos sus autocorrelaciones. Tenga en cuenta que a partir de la estructura del MCO es cierto que $\mathbb{E} \left( \widehat{\boldsymbol{\varepsilon}} | \mathbf{X}\right) = \boldsymbol{0}$.

Correlograma residual

Entonces, suponga que queremos calcular la correlación muestral entre $\widehat{\epsilon}_i$ y $\widehat{\epsilon}_{i - k}$. Para calcular la autocorrelación de orden $k$ tenemos:

\[ \widehat{\rho}(k) = \widehat{\mathbb{C}{\rm orr}}(\widehat{\epsilon}_i, \widehat{\epsilon}_{i-k}) = \dfrac{\widehat{\mathbb{C}{\rm ov}}(\widehat{\epsilon}_i, \widehat{\epsilon}_{i-k})}{\sqrt{\widehat{\mathbb{V}{\rm ar}}(\widehat{\epsilon}_i)}\sqrt{\widehat{\mathbb{V}{\rm ar}}(\widehat{\epsilon}_{i-k})}} = \dfrac{\sum_{i = k + 1}^N \widehat{\epsilon}_i \widehat{\epsilon}_{i-k}}{\sum_{i = 1}^N \widehat{\epsilon}_i^2} \]

Correlograma residual

Además, suponga que estamos interesados en probar si la autocorrelación de orden k $\rho(k)$ es significativamente diferente de cero:

\[ \begin{aligned} H_0: \rho(k) = 0\\ H_1: \rho(k) \neq 0 \end{aligned} \]

Bajo la hipótesis nula, es cierto que $\widehat{\rho}(k) \stackrel{a}{\sim} \mathcal{N}(0, 1/N)$, donde $\stackrel{a}{\sim}$ indica distribución asintótica: la distribución cuando el tamaño de la muestra $N \rightarrow \infty$.

Correlograma residual

En este caso, significa que para muestras grandes la distribución es aproximadamente normal. En consecuencia, una estadística adecuada se puede construir como:

\[ Z = \dfrac{\widehat{\rho}(k) - 0}{\sqrt{1/N}} = \widehat{\rho}(k) \cdot \sqrt{N} \stackrel{a}{\sim} \mathcal{N}(0, 1) \]

Correlograma residual

Entonces, con un nivel de signficancia de $5\%$, el valor crítico es $Z_c \approx 1.96$. En este caso, se rechazaría la hipótesis nula cuando:

\[ \widehat{\rho}(k) \cdot \sqrt{N} \geq 1.96, \text{ o } \widehat{\rho}(k) \cdot \sqrt{N} \leq - 1.96 \]

o alternativamente, si queremos tener la notación similar a “$\text{estimación} \pm \text{valor crítico} \cdot \text{se}(\text{estimación})$”, lo podemos escribir como: \[ \widehat{\rho}(k) \pm 1.96 / \sqrt{N} \]

Ejemplo con R

Como regla general, a veces podemos tomar $2$ en vez de $1.96$. También podemos calcular esto a través del R:

z_c <- qnorm(p = 1 - 0.05 / 2, mean = 0, sd = 1)
print(z_c)

## [1] 1.959964

Ejemplo con R

Comenzaremos calculando manualmente las autocorrelaciones y sus límites de confianza para $k = 1,...,20$:

rho <- NULL
e <- mdl_1$residuals
for(k in 1:20){
  r <- sum(e[-c(1:k)] * e[1:(N-k)]) / sum(e^2)
  rho <- c(rho, r)
}
rho_lower <- - z_c / sqrt(N)
rho_upper <- + z_c / sqrt(N)

Ejemplo con R

y podemos graficarlos:

plot(1:k, rho, ylim = c(min(rho_lower), 1), pch = 21, bg = "red")
segments(x0 = 1:k, y0 = rho, x1 = 1:k, y1 = 0, lwd = 2)
abline(h = 0)
# Graficando los límites de confianza:
abline(h = rho_upper, lty = 2, col = "blue")

abline(h = rho_lower, lty = 2, col = "blue")

Ejemplo con R

Alternativamente, podemos usar las funciones disponibles en el R:

forecast::Acf(e)

Ejemplo con R

Podemos ver en los gráficos que hay algunas autocorrelaciones de la muestra $\widehat{\rho}(k)$, que son significativamente diferentes de cero (es decir, sus valores están por encima de la línea horizontal azul)

Ejemplo con R

Al igual que con la mayoría de las herramientas de diagnóstico visual, es posible que no siempre sea una respuesta clara solo a partir de los gráficos. Por tanto, podemos aplicar una serie de pruebas estadísticas diferentes para comprobar si existen correlaciones muestrales estadísticamente significativas.

Pruebas de autocorrelación

Prueba de Durbin-Watson (DW)

La prueba DW se utiliza para probar la hipótesis de que los residuales están autocorrelacionados de orden 1, es decir, que en el siguiente modelo:

\[ \epsilon_i = \rho \epsilon_{i-1} + v_i \]

la hipótesis que se está probando es:

\[ \begin{cases} H_0: \rho = 0\qquad \text{(no hay autocorrelación)}\\ H_1: \rho \neq 0\qquad \text{(autocorrelación de orden 1)} \end{cases} \]

Prueba de Durbin-Watson (DW)

Como no observamos los verdaderos términos del error, usamos los residuales de MCO, osea $\widehat{\epsilon}_i$ y calculamos la estadística de Durbin-Watson como:

\[ DW = \dfrac{\sum_{i = 2}^N (\widehat{\epsilon}_i - \widehat{\epsilon}_{i-1})^2}{\sum_{i = 1}^N \widehat{\epsilon}_i^2} \]

Como regla rápida, si el estadístico DW está cerca de $2$, entonces rechazamos la hipótesis nula de que no hay autocorrelación.

Prueba de Durbin-Watson (DW)

Si asumimos que $\sum_{i = 2}^N \widehat{\epsilon}_i^2 \approx \sum_{i = 2}^N \widehat{\epsilon}_{i-1}^2$, entonces podemos reescribir la estadística DW como:

\[ DW = \dfrac{\sum_{i = 2}^N \widehat{\epsilon}_i^2 - 2 \sum_{i = 2}^N \widehat{\epsilon}_i\widehat{\epsilon}_{i-1} + \sum_{i = 2}^N \widehat{\epsilon}_{i-1}^2}{\sum_{i = 1}^N \widehat{\epsilon}_i^2} \approx \dfrac{2\left[ \sum_{i = 2}^N \widehat{\epsilon}_i^2 - \sum_{i = 2}^N \widehat{\epsilon}_i\widehat{\epsilon}_{i-1} \right]}{\sum_{i = 1}^N \widehat{\epsilon}_i^2} = 2 \left[ 1 - \widehat{\rho}(1)\right] = \begin{cases} 4,\ \text{ if } \widehat{\rho}(1) = -1\\ 2,\ \text{ if } \widehat{\rho}(1) = 0\\ 0,\ \text{ if } \widehat{\rho}(1) = 1 \end{cases} \]

lo que ayuda a comprender por qué esperamos que el estadístico DW esté cerca de $2$ bajo la hipótesis nula.

Ejemplo con R

print(lmtest::dwtest(mdl_1, alternative = "two.sided"))

## 
##  Durbin-Watson test
## 
## data:  mdl_1
## DW = 0.56637, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is not 0

porque $p-valor < 0.05$, rechazamosla hipótesis nula a una nivel de significancia de $5\%$ y concluimos que los residuos están autocorrelacionados en orden 1.

Prueba de Durbin-Watson (DW)

Como no observamos los verdaderos términos del error, usamos los residuales de MCO, osea $\widehat{\epsilon}_i$ y calculamos la estadística de Durbin-Watson como:

\[ DW = \dfrac{\sum_{i = 2}^N (\widehat{\epsilon}_i - \widehat{\epsilon}_{i-1})^2}{\sum_{i = 1}^N \widehat{\epsilon}_i^2} \]

Como regla rápida, si el estadístico DW está cerca de $2$, entonces rechazamos la hipótesis nula de que no hay autocorrelación.

Prueba Q de Ljung-Box

La prueba Q de Ljung-Box (alguna veces llamada de Prueba de Portmanteau) es usada para probar si las observaciones en el tiempo son aleatorias e independiente. En particular, para un $k$ dado, se prueba lo siguiente:

\[ \begin{align*} H_{0}& :\text{ las autocorrelaciones hasta el retraso $k$ son todas 0}\\ H_{A}& :\text{ las autocorrelaciones de $k$ o más retrasos difieren de 0} \end{align*} \]

Prueba Q de Ljung-Box

La estadística de prueba es calculada como:

\[Q_{k} = n(n+2)\sum^{k}_{j=1}\frac{\widehat{\rho}^{2}_{j}}{n-j}\]

Que es aproximadamente distribuida como una $\chi^{2}_{k}$.

Ejemplo con R

print(Box.test(residuals(mdl_1), lag = 1, type = "Ljung"))

## 
##  Box-Ljung test
## 
## data:  residuals(mdl_1)
## X-squared = 101.11, df = 1, p-value < 2.2e-16

porque $p-valor < 0.05$, rechazamosla hipótesis nula a una nivel de significancia de $5\%$ y concluimos que los residuos están autocorrelacionados en orden 1.

Prueba de Breusch-Godfrey (BG)

La prueba de BG puede aplicarse al modelo, con una variable de respuesta con un lag (retardo) en el lado derecho, por ejemplo: $Y_i = \beta_0 + \beta_1 X_{1,i} + ... + \beta_k X_{k,i} + \beta_{k+1} Y_{i - 1} + \epsilon_i$.

En general, estimamos los parámetros del modelo a través de MCO y calculamos los residuos: \[ \widehat{\boldsymbol{\varepsilon}} = \mathbf{Y} - \mathbf{X} \widehat{\boldsymbol{\beta}} \] donde $\mathbf{X}$ que contiene $X_{1,i},...,X_{k,i}, Y_{i - 1}$ para $i = 1,...,N$.

Prueba de Breusch-Godfrey (BG o LM)

Entonces, estimamos una regresión auxiliar en $\widehat{\boldsymbol{\varepsilon}}$ como: \[ \widehat{\epsilon}_i = \alpha_0 + \alpha_1 X_{1,i} + ... + \alpha_k X_{k,i} + \rho_i \widehat{\epsilon}_{i-1} + ... + \rho_p \widehat{\epsilon}_{i-p} + u_i \] (Nota: si incluimos $Y_{i - 1}$ en la regresión inicial, entonces necesitamos incluirlos también en la regresión auxiliar).

Queremos probar la hipótesis nula de que los coeficientes residuales rezagados son insignificantes:

\[ \begin{cases} H_0: \rho_1 = ... = \rho_p = 0\\ H_1: \rho_j \neq 0, \text{ para algún } j \end{cases} \]

Prueba de Breusch-Godfrey (BG o LM)

Podemos realizar la prueba-$F$, o la prueba-Chi-cuadrado: \[ LM = (N-p)R_{\widehat{\epsilon}}^2 \sim \chi^2_p \] donde $R_{\widehat{\epsilon}}^2$ es el $R$-cuadrado de la regresión auxiliar en $\widehat{\boldsymbol{\varepsilon}}$.

La prueba BG es una prueba más general en comparación con DW.

Ejemplo con R

Observando las gráficas del correlograma. Probemos la hipótesis nula de que al menos uno de los primeros tres rezagos de los residuos no es cero, es decir $p = 3$:

bg_t <- lmtest::bgtest(mdl_1, order = 3)
print(bg_t)

## 
##  Breusch-Godfrey test for serial correlation of order up to 3
## 
## data:  mdl_1
## LM test = 100.72, df = 3, p-value < 2.2e-16

Desde que $p-valor < 0.05$, rechazamos la hipótesis nula y concluimos que existe una autocorrelación en los residuos.

Pruebas de autocorrelación: Comparación

Característica	Durbin-Watson	Ljung-Box	Breusch-Godfrey
Propósito	Autocorrelación de primer orden	Autocorrelación en varios rezagos	Autocorrelación de orden superior
Orden de Autocorrelación	Primer orden	Múltiples (conjunto)	Superior (específico)
Variables Rezagadas	No	Sí	Sí
Distribución	Normal	$\chi^2$ con $k$ g.l.	$\chi^2$ con $p$ g.l.
Ventajas	Simple, directa	Múltiples rezagos	Orden superior, modelos con rezagados
Desventajas	Solo primer orden, no rezagados	No específica, menos intuitiva	Regresión adicional, más compleja

Heterocedasticidad y Autocorrelación-Errores estándar consistentes (HAC)

Hasta ahora, hemos asumido que los elementos diagonales de $\mathbf{\Sigma}$ son constantes, es decir, que los residuos son autocorrelacionados pero homocedásticos. Podemos generalizar aún más esto para el caso de errores estándar heterocedásticos y autocorrelacionados:

\[ \mathbf{\Sigma} = \begin{bmatrix} \sigma_1^2 & \sigma_{1,2} & \sigma_{1,3} & ... & \sigma_{1,N}\\ \sigma_{2,1} & \sigma_2^2 & \sigma_{2,3} & ... & \sigma_{2,N}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \sigma_{N,1} & \sigma_{N,2} & \sigma_{N3,} & ... & \sigma_N^2 \end{bmatrix}, \quad \sigma_{i,j} = \sigma_{j,i} = \mathbb{C}{\rm ov}(\epsilon_i, \epsilon_j) \neq 0, \quad i,j = 1,...,N \]

Heterocedasticidad y Autocorrelación-Errores estándar consistentes (HAC)

podemos usar el estimador MCO y corregir los errores estándar. En este caso, los errores estándar corregidos se conocen como errores estándar HAC (heterocedasticidad y autocorrelación consistente) o errores estándar de Newey-West.

Heterocedasticidad y Autocorrelación-Errores estándar consistentes (HAC)

La matriz de covarianza de White asume residuos no autocorrelacionados. Por otro lado, Newey-West propuso un estimador de covarianza más general, que es robusto tanto a la heterocedasticidad como a la autocorrelación, de los residuos con covarianza desconocida. El estimador de covarianza del coeficiente HAC maneja la autocorrelación con rezagos de hasta orden $p$; se supone que los rezagos mayores que $p$ son no significativos y, por lo tanto, pueden ignorarse. Se define como:

Heterocedasticidad y Autocorrelación-Errores estándar consistentes (HAC)

\[ HAC_{\widehat{\boldsymbol{\beta}}_{MCO}} = \left( \mathbf{X}^\top \mathbf{X}\right)^{-1} N \widehat{\mathbf{\Sigma}} \left( \mathbf{X}^\top \mathbf{X}\right)^{-1} \]

donde:

\[ \widehat{\mathbf{\Sigma}} = \widehat{\boldsymbol{\Gamma}}(0) + \sum_{j = 1}^p \left[ 1 - \dfrac{j}{p + 1}\right] \left[ \widehat{\boldsymbol{\Gamma}}(j) + \widehat{\boldsymbol{\Gamma}}(-j)\right] \] \[ \widehat{\boldsymbol{\Gamma}}(j) = \dfrac{1}{N} \left( \sum_{i = 1}^N\widehat{\epsilon}_{i}\widehat{\epsilon}_{i-j} \mathbf{X}_i\mathbf{X}_{i-j}^\top\right),\quad \mathbf{X}_i = \left[1,\ X_{1,i},\ X_{2,i},\ ...,\ X_{k,i} \right]^\top \]

Heterocedasticidad y Autocorrelación-Errores estándar consistentes (HAC)

En ausencia de correlación serial tendríamos que: \[ \widehat{\mathbf{\Sigma}} = \widehat{\boldsymbol{\Gamma}}(0) \] que es equivalente a los estimadores de White (HCE).

Nota: HAC no solo corrige la autocorrelación, sino también la heterocedasticidad.

No se alarme si ve errores estándar de HAC ligeramente diferentes en diferentes programas estadísticos; hay una serie de variaciones diferentes de $\widehat{\mathbf{\Sigma}}$.

Ejemplo con R

Usando las funciones del R tenemos que:

#V_HAC <- sandwich::vcovHAC(mdl_1)
V_HAC <- sandwich::NeweyWest(mdl_1, lag = 1)
print(V_HAC)

##             (Intercept)           x1           x2
## (Intercept)  1.43560709 -0.303808295 -0.043572462
## x1          -0.30380830  0.183608167 -0.002018108
## x2          -0.04357246 -0.002018108  0.004421496

Ejemplo con R

entonces los coeficientes, sus errores estándar y $p$-valores se pueden resumir como:

# print(lmtest::coeftest(mdl_1, sandwich::vcovHAC(mdl_1)))
print(lmtest::coeftest(mdl_1, sandwich::NeweyWest(mdl_1, lag = 1))[, ])

##              Estimate Std. Error    t value      Pr(>|t|)
## (Intercept) 10.616005 1.19816822   8.860196  4.719858e-16
## x1           4.890185 0.42849524  11.412460  1.741811e-23
## x2          -2.936860 0.06649433 -44.167068 3.799730e-104

Note que usamos [, ] para que podamos tratar la salida como un data.frame si necesitáramos redondear los resultados.

Ejemplo con R

Comparamos los residuales sesgados en la salida de MCO:

print(coef(summary(mdl_1)))

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 10.616005 0.97982714  10.83457 9.311699e-22
## x1           4.890185 0.20456771  23.90497 3.887306e-60
## x2          -2.936860 0.07545421 -38.92241 1.933624e-94

vemos que los errores estándar de HAC son algo mayores, en comparación con los que obtenemos con MCO.

GLS factible - Procedimiento Cochrane-Orcutt (CORC)

Alternativamente al HAC, podemos hacer algunas suposiciones sobre la naturaleza de la autocorrelación y emplear un estimador GLS más eficiente.

Para el caso, cuando los residuos están autocorrelacionados en el lap 1, pero no heteroscedásticos:

\[ \epsilon_i = \rho \epsilon_{i-1} + u_i,\ |\rho| < 1,\ u_i \sim \mathcal{N}(0, \sigma^2), \quad i \in \{0, \pm 1, \pm 2,... \} \]

GLS factible - Procedimiento Cochrane-Orcutt (CORC)

notamos que:

\[ \begin{aligned} \mathbb{V}{\rm ar}(\epsilon_i) &= \mathbb{V}{\rm ar}(\rho\epsilon_{i-1} + u_i)\\ &= \mathbb{V}{\rm ar}(\rho(\rho\epsilon_{i-2} + u_{i-1}) + u_i) \\ &= \mathbb{V}{\rm ar}(\rho^2(\rho\epsilon_{i-3} + u_{i-2}) + u_i + \rho u_{i-1})\\ &= ... \\ &= \mathbb{V}{\rm ar} \left(\sum_{j = 0}^\infty \rho^j u_{i-j}\right)\\ &=\sum_{j = 0}^\infty \rho^{2j} \cdot \mathbb{V}{\rm ar} \left(u_{i-j}\right)\\ &= \dfrac{\sigma^2}{1 - \rho^2} \end{aligned} \]

GLS factible - Procedimiento Cochrane-Orcutt (CORC)

\[ \begin{aligned} \mathbb{C}{\rm ov}(\epsilon_i,\ \epsilon_{i - k}) &= \mathbb{C}{\rm ov}(\rho(\rho\epsilon_{i-2} + u_{i-1}) + u_i,\ \epsilon_{i - k}) \\ &=\mathbb{C}{\rm ov}(\rho^2(\rho\epsilon_{i-3} + u_{i-2}) + \rho u_{i-1},\ \epsilon_{i - k}) \\ &= ... \\ &=\mathbb{C}{\rm ov}(\rho^k\epsilon_{i-k},\ \epsilon_{i - k}) \\ &= \rho^k \mathbb{C}{\rm ov}(\epsilon_{i-k},\ \epsilon_{i - k}) \\ &= \rho^k \sigma^2, \quad \text{ dado que } \mathbb{C}{\rm ov}(u_i, \epsilon_j ) = 0, i \neq j \end{aligned} \]

GLS factible - Procedimiento Cochrane-Orcutt (CORC)

En consecuencia, podemos reescribir la matriz de covarianza como:

\[ \mathbf{\Sigma} = \dfrac{\sigma^2}{1 - \rho^2} \begin{bmatrix} 1 & \rho & \rho^2 & ... & \rho^{N-1}\\ \rho & 1 & \rho & ... & \rho^{N-2}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ \rho^{N-1} & \rho^{N-2} & \rho^{N-3} & ... & 1 \end{bmatrix} \]

En este caso, el conocimiento del parámetro $\rho$ nos permite aplicar empíricamente el GLS - como FGLS.

Estimador Cochrane–Orcutt (CORC):

Estimamos $\boldsymbol{\beta}$ via MCO de:

\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]

y calculamos los residuos $\widehat{\boldsymbol{\varepsilon}}^{(1)} = \mathbf{Y} - \mathbf{X} \widehat{\boldsymbol{\beta}}$, donde ${(1)}$ denota la primera iteración.

Estime la siguiente regresión residual a través de MCO: \[ \widehat{\epsilon}^{(1)}_i = \rho \widehat{\epsilon}^{(1)}_{i-1} + \widehat{u}_i \] y obtenemos $\widehat{\rho}^{(1)}$.
Calcule las siguientes variables transformadas:

\[ Y_i^* = Y_i - \widehat{\rho}^{(1)} Y_{i-1}, \qquad \boldsymbol{X}_i^* = \boldsymbol{X}_i - \widehat{\rho}^{(1)} \boldsymbol{X}_{i-1},\text{ donde }\mathbf{X}_i = \left[1,\ X_{1,i},\ X_{2,i},\ ...,\ X_{k,i} \right]^\top \]

Estimador Cochrane–Orcutt (CORC):

Aplicar MCO al siguiente modelo:

\[ (Y_i - \widehat{\rho}^{(1)} Y_{i-1}) = \beta_0 (1 - \widehat{\rho}^{(1)}) + \beta_1 (X_{1,i} - \widehat{\rho}^{(1)}{X}_{i-1}) + ... + \beta_k (X_{k,i} - \widehat{\rho}^{(1)}{X}_{k-1}) + u_i \]

o, en forma más compacta: \[ \mathbf{Y}^* = \mathbf{X}^* \boldsymbol{\beta} + \boldsymbol{u} \] para obtener las estimaciones de MCO de $\widetilde{\boldsymbol\beta}$.

Estimador Cochrane–Orcutt (CORC):

Nota: si decidimos usar una columna de unos para la constante $\beta_0$, entonces realmente habremos estimado $\widetilde{\beta}_0^* = \widetilde{\beta}_0 (1 - \widehat{\rho}^{(1)})$ y necesitará transformar el término de intersección como $\widetilde{\beta}_0 = \widehat{\beta}_0^* / (1 - \widehat{\rho}^{(1)})$. Por otro lado, si transformamos la columna de intersección para tener $1 - \widehat{\rho}^{(1)}$ en cambio, en la matriz de diseño, entonces no necesitaremos transformar el coeficiente de intersección.

Estimador Cochrane–Orcutt (CORC):

Habiendo estimado el modelo, calcule los residuos sobre los datos no transformados: $\widehat{\boldsymbol{\varepsilon}}^{(2)} = \mathbf{Y} - \mathbf{X} \widetilde{\boldsymbol{\beta}}$ y vaya al paso (2). Repita este procedimiento hasta $\widehat{\rho}$ que converja: es decir si el cambio en $\widehat{\rho}^{(K)}$, en comparación con la iteración anterior $\widehat{\rho}^{(K-1)}$ no es más que $\Delta$ (por ejemplo $\Delta = 0.001$) detenga el procedimiento.

Estimador Cochrane–Orcutt (CORC):

El valor final $\widehat{\rho}^{(K)}$ luego se utiliza para obtener las estimaciones de FGLS (CORC) de $\widehat{\boldsymbol\beta}$. Nuevamente, dependiendo de su especificación, es posible que deba transformar el coeficiente de intersección:$\widehat{\beta}_0 = \widehat{\beta}_0^* / (1 - \widehat{\rho}^{(K)})$.

Estos estimadores de FGLS no son insesgados pero son consistentes y asintóticamente eficientes.

Ejemplo con R

El procedimiento se puede realizar con las funciones del R de la siguiente manera:

mdl_1_CORC <- orcutt::cochrane.orcutt(mdl_1)
#
print(coef(summary(mdl_1_CORC)))

##              Estimate Std. Error t value      Pr(>|t|)
## (Intercept) 10.402800 1.50824063   6.897  7.124332e-11
## x1           5.082365 0.48877934  10.398  1.899433e-20
## x2          -2.978490 0.04202203 -70.879 1.185668e-141

Ejemplo con R

la estimación de $\rho$ es:

print(mdl_1_CORC$rho)

## [1] 0.7088817

que está cerca del valor real de $\rho$ usado en la generación de datos.

Elegiendo entre HAC y FGLS

Se ha vuelto más popular estimar modelos por MCO y corregir los errores estándar para formas bastante arbitrarias de autocorrelación y heterocedasticidad.

Elegiendo entre HAC y FGLS

Con respecto a HAC:

El cálculo de errores estándar robustos generalmente funciona bien en grandes conjuntos de datos;
Con el aumento en el poder computacional, no solo ha sido posible estimar (rápidamente) modelos en grandes conjuntos de datos, sino también calcular sus estimadores de covarianza robusta (HAC);
Si bien FGLS ofrece una eficiencia teórica, implica hacer suposiciones adicionales sobre la matriz de covarianza de los errores, que pueden no ser fáciles de probar/verificar, lo que puede amenazar la consistencia del estimador.

Elegiendo entre HAC y FGLS

Respecto a FGLS:

si incluimos $Y_{i-1}$ en el lado derecho de la ecuación - el FGLS no solo es ineficiente, sino que también es inconsistente;
en la mayoría de las aplicaciones de FGLS, se supone que los errores siguen un proceso autorregresivo de primer orden. Puede ser mejor evaluar las estimaciones de MCO y utilizar una corrección sólida en sus errores estándar para formas más generales de autocorrelación;
Además de imponer una suposición de la estructura de covarianza residual con respecto a la autocorrelación, GLS también requiere la suposición S3 para mantenerse, a diferencia de HAC.

Simulación Monte Carlo: MCO vs FGLS

Para ilustrar los efectos de la autocorrelación sobre los errores estándar de las estimaciones y la eficiencia entre MCO y FGLS, realizaremos una simulación de Monte Carlo. Simularemos el siguiente modelo:

\[ \begin{cases} Y_{i}^{(m)} &= \beta_0 + \beta_1 X_{1,i}^{(m)} + \beta_2 X_{2,i}^{(m)} + \epsilon_i^{(m)} \\ \epsilon_i^{(m)} &= \rho \epsilon_{i-1}^{(m)} + u_i^{(m)} ,\ |\rho| < 1,\ u_i^{(m)} \sim \mathcal{N}(0, \sigma^2),\quad \epsilon_0^{(m)} = 0,\quad i = 1,...,N,\quad m = 1,..., MC \end{cases} \]

Nosotros simularemos un total de $MC = 1000$ muestras de este modelo con coeficientes específicos y estimaremos los parámetros a través de los métodos estudiados. Lo haremos con el siguiente código:

Ejemplo con R

set.seed(123)
# Número de simulaciones:
MC <- 1000
# Fijamos los parámetros
N <- 100
beta_vec <- c(1, 0.2, -0.3)
rho <- 0.8
# matriz de parámetros para cada muestra:
beta_est_ols  <- NULL
beta_pval_ols <- NULL
beta_pval_hac <- NULL
#
beta_est_fgls  <- NULL
beta_pval_fgls <- NULL
#
beta_est_gls  <- NULL
beta_pval_gls <- NULL
#
for(i in 1:MC){
  # simulamos los datos
  x1 <- seq(from = 0, to = 5, length.out = N)
  x2 <- sample(seq(from = 3, to = 17, length.out = 80), size = N, replace = TRUE)
  ee <- rnorm(mean = 0, sd = 3, n = N)
  for(i in 2:N){
    ee[i] <- rho * ee[i-1] + ee[i]
  }
  y  <- cbind(1, x1, x2) %*% beta_vec + ee
  data_mat <- data.frame(y, x1, x2)
  # estimamos via MCO
  mdl_0 <- lm(y ~ x1 + x2, data = data_mat)
  # corregimos el error estándar de MCO:
  # mdl_hac <- lmtest::coeftest(mdl_0, vcov. = sandwich::vcovHAC(mdl_0))
  mdl_hac <- lmtest::coeftest(mdl_0, vcov. = sandwich::NeweyWest(mdl_1, lag = 1))
  # estimamos via CORC:
  suppressWarnings({
    mdl_fgls <- orcutt::cochrane.orcutt(mdl_0)
  })
  # Guarde las estimaciones de cada muestra:
  beta_est_ols  <- rbind(beta_est_ols, coef(summary(mdl_0))[, 1])
  beta_est_fgls  <- rbind(beta_est_fgls, coef(summary(mdl_fgls))[, 1])
  # Guarde los p-valores de cada muestra:
  beta_pval_ols <- rbind(beta_pval_ols, coef(summary(mdl_0))[, 4])
  beta_pval_hac <- rbind(beta_pval_hac, mdl_hac[, 4])
  beta_pval_fgls <- rbind(beta_pval_fgls, coef(summary(mdl_fgls))[, 4])
}

Ejemplo con R

En cuanto a la tasa de rechazo de la hipótesis nula de que un parámetro no significativo, tenemos que:

alpha = 0.05
a1 <- colSums(beta_pval_ols < alpha) / MC
a2 <- colSums(beta_pval_hac < alpha) / MC
a3 <- colSums(beta_pval_fgls < alpha) / MC
#
a <- t(data.frame(a1, a2, a3))
colnames(a) <- names(a1)
rownames(a) <- c("OLS: H0 rejection rate", "HAC: H0 rejection rate", "FGLS: H0 rejection rate")
print(a)

##                         (Intercept)    x1    x2
## OLS: H0 rejection rate        0.373 0.505 0.757
## HAC: H0 rejection rate        0.457 0.365 0.919
## FGLS: H0 rejection rate       0.135 0.122 0.999

Ejemplo con R

Entonces, el FGLS parece ser mejor en algunos casos de parámetros, mientras que HAC puede ser similar al MCO. Con todo, si solo tenemos una autocorrelación del problema de orden 1, es posible que las correcciones no tengan un gran impacto en nuestras conclusiones en algunos casos.
Podemos mirar los histogramas de estimación de coeficientes:

Ejemplo con R

# histogramas de las estimaciones:
layout(matrix(c(1, 2, 3, 3), nrow = 2, byrow = TRUE))
for(j in ncol(beta_est_ols):1){
  p1 <- hist(beta_est_ols[, j], plot = FALSE, breaks = 25)
  p2 <- hist(beta_est_fgls[, j], plot = FALSE, breaks = 25)
  #
  plot(p1, col = rgb(0,0,1,1/4), main = colnames(beta_est_ols)[j], ylim = c(0, max(p1$counts, p2$counts)))
  plot(p2, col = rgb(1,0,0,1/4), add = T)
  abline(v = beta_vec[j], lwd = 3, lty = 2, col = "red")
  #
  legend("topleft", legend = c("True Value", "OLS", "FGLS"), 
         lty = c(3, NA, NA), lwd = c(3, 0, 0), pch = c(NA, 22, 22),
         col = c("red", "black", "black"), pt.cex = 2,
         pt.bg = c(NA, rgb(0,0,1,1/4), rgb(1,0,0,1/4))
         )
}

Ejemplo con R

Para una estructura de correlación de orden superior y más compleja de los residuales, este puede no ser siempre el caso, por lo tanto, si detectamos una autocorrelación, deberíamos considerar ello de alguna manera.

Referencias

Robert Grant - Data Visualization - Charts, Maps, and Interactive Graphics (2018)
https://little-book-of-r-for-multivariate-analysis.readthedocs.io/en/latest/src/multivariateanalysis.html
https://cran.r-project.org/web/packages/naniar/vignettes/naniar-visualisation.html