INTRODUCCIÓN

Modelo Lineal Clásico

El modelo lineal clásico es definido como

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

donde :

• \(\mathbf{Y}\) es un vector de datos observables (variable respuesta);

• \(\boldsymbol{\beta}\) es un vector desconocido de parámetros;

• \(\mathbf{X}\) es la matriz de diseño;

• \(\boldsymbol{\varepsilon}\) es el vector de errores aleatorios y \(\boldsymbol{\varepsilon}\sim N(\mathbf{0}, \sigma^2 \mathbf{I})\).

Modelo Lineal Clásico

• Hay muchas situaciones reales en las que estas suposiciones no son apropiadas. En algunos casos, el sistema de medición utilizado podría ser una fuente de variabilidad, y el error es proporcional a la cantidad medida. Otras veces esto ocurre cuando los errores están correlacionados.

• Además, cuando la distribución subyacente es continua, pero asimétrica, como lognormal, gamma, etc., la varianza no es constante y, en muchos casos, la varianza es una función de la media.

Modelo Lineal Clásico

• Un punto importante es que la varianza constante está vinculada a la suposición de distribución normal para la respuesta.

• En la presencia de heterocedasticidad la variancia de los errores será diferente para cada valor de \(X\) y los estimadores obtenidos usando mínimos cuadrados ordinarios aunque permanecen insesgados y consistentes dejan de ser eficientes (no presentan variancia mínima).

REGRESIÓN LINEAL MÚLTIPLE GENERAL (RLMG)

Mínimos cuadrados generalizados

Regresión lineal múltiple general (RLMG)

Regresión lineal múltiple general (RLMG)

Regresión lineal múltiple general (RLMG)

Consecuencias cuando se usa MCO en una RLMG

Consecuencias cuando se usa MCO en una RLMG

Consecuencias cuando se usa MCO en una RLMG

Consecuencias cuando se usa MCO en una RLMG

Consecuencias cuando se usa MCO en una RLMG

Consecuencias cuando se usa MCO en una RLMG

Mínimos cuadrados generalizados (MCG)

• Luego, premultiplicar nuestro modelo inicial a través de esta matriz produce: \[\begin{equation} \mathbf{\Psi}^\top \mathbf{Y} = \mathbf{\Psi}^\top\mathbf{X} \boldsymbol{\beta} + \mathbf{\Psi}^\top\boldsymbol{\varepsilon} \tag{4.1} \end{equation}\]
• Porque \(\mathbf{\Omega}\) es no singular, así como \(\mathbf{\Psi}\), por lo tanto, siempre podemos multiplicar la expresión anterior por \(\left(\mathbf{\Psi}^\top\right)^{-1}\) para volver al modelo inicial. En otras palabras, simplemente hemos escalado ambos lados de la ecuación inicial.

Mínimos cuadrados generalizados (MCG)

• Entonces, las siguientes propiedades son válidas para los residuos: \[ \begin{aligned} \mathbb{E} \left( \mathbf{\Psi}^\top\boldsymbol{\varepsilon} \right) &= \mathbf{\Psi}^\top \mathbb{E} \left( \boldsymbol{\varepsilon} \right) = \boldsymbol{0}\\\\ \mathbb{V}{\rm ar} \left( \mathbf{\Psi}^\top\boldsymbol{\varepsilon} \right) &= \mathbb{V}{\rm ar} \left( \mathbf{\Psi}^\top\boldsymbol{\varepsilon} \right)\\ &= \mathbf{\Psi}^\top \mathbb{V}{\rm ar} \left( \boldsymbol{\varepsilon} \right) \mathbf{\Psi} \\ &= \sigma^2 \mathbf{\Psi}^\top \mathbf{\Omega} \mathbf{\Psi} \\ &= \sigma^2 \mathbf{\Psi}^\top \left(\mathbf{\Omega}^{-1} \right)^{-1} \mathbf{\Psi} \\ &= \sigma^2 \mathbf{\Psi}^\top \left(\mathbf{\Psi} \mathbf{\Psi}^\top \right)^{-1} \mathbf{\Psi} \\ &= \sigma^2 \mathbf{\Psi}^\top \left(\mathbf{\Psi}^\top \right)^{-1} \mathbf{\Psi} ^{-1} \mathbf{\Psi} \\ &= \sigma^2 \mathbf{I} \end{aligned} \] lo que significa que ahora las suposiciones (S3) - (S4) se cumplen para el modelo, definido en (4.1).

Mínimos cuadrados generalizados (MCG)

• En consecuencia, el estimador de MCO de \(\boldsymbol{\beta}\) de la regresión en (4.1) es: \[ \begin{aligned} \widehat{\boldsymbol{\beta}} &= \left( \mathbf{X}^\top \mathbf{\Psi} \mathbf{\Psi}^\top \mathbf{X}\right)^{-1} \mathbf{X}^\top \mathbf{\Psi} \mathbf{\Psi}^\top \mathbf{Y} = \left( \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{X}\right)^{-1} \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{Y} \end{aligned} \] Este estimador se llama estimador de mínimos cuadrados generalizados (MCG) de \(\boldsymbol{\beta}\). Entonces, el estimador de MCG se puede formular de la siguiente manera:

Mínimos cuadrados generalizados (MCG)

• Sea \(\mathbf{\Omega}^{-1} = \mathbf{\Psi} \mathbf{\Psi}^\top\) y \(\mathbf{Y}^* = \mathbf{X}^* \boldsymbol{\beta} + \boldsymbol{\varepsilon}^*\), donde \(\mathbf{Y}^* = \mathbf{\Psi}^\top \mathbf{Y}\), \(\mathbf{X}^* = \mathbf{\Psi}^\top\mathbf{X}\) y \(\boldsymbol{\varepsilon}^* =\mathbf{\Psi}^\top\boldsymbol{\varepsilon}\).

• Entonces el estimador de MCG de \(\boldsymbol{\beta}\) es: \[ \widehat{\boldsymbol{\beta}}_{MCG} = \left( \mathbf{X}^{*\top} \mathbf{X}^* \right)^{-1} \mathbf{X}^{*\top} \mathbf{Y}^* = \left( \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{X}\right)^{-1} \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{Y} \]

Mínimos cuadrados generalizados (MCG)

• Alternativamente, también se puede obtener como una solución al minimizar la función de criterio de MCG: \[ \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right)^\top \mathbf{\Omega}^{-1} \left( \mathbf{Y} - \mathbf{X} \boldsymbol{\beta} \right) \rightarrow \min_{\boldsymbol{\beta}} \]
• Esta función de criterio puede considerarse como una generalización de la función (\(RSS\)), que se minimiza en el caso del MCO. El efecto de tal ponderación es claro cuando \(\mathbf{\Omega}\) es diagonal.
• A cada observación se le asigna simplemente un peso proporcional a la inversa de la varianza del error.

Mínimos cuadrados generalizados (MCG)

• A continuación, necesitamos estimar la varianza del error: \[ \begin{aligned} \widehat{\sigma}^2 &= \dfrac{1}{N-(k+1)} \left( \mathbf{Y}^* - \mathbf{X}^* \widehat{\boldsymbol{\beta}}_{MCG}\right)^\top \left( \mathbf{Y}^* - \mathbf{X}^* \widehat{\boldsymbol{\beta}}_{MCG}\right) \\ &= \dfrac{1}{N-(k+1)} \left( \mathbf{\Psi}^\top \mathbf{Y} - \mathbf{\Psi}^\top\mathbf{X} \widehat{\boldsymbol{\beta}}_{MCG}\right)^\top \left( \mathbf{\Psi}^\top \mathbf{Y} - \mathbf{\Psi}^\top\mathbf{X}\widehat{\boldsymbol{\beta}}_{MCG}\right) \\ &= \dfrac{1}{N-(k+1)} \left( \mathbf{Y} - \mathbf{X} \widehat{\boldsymbol{\beta}}_{MCG}\right)^\top \mathbf{\Psi} \mathbf{\Psi}^\top \left( \mathbf{Y} - \mathbf{X}\widehat{\boldsymbol{\beta}}_{MCG}\right) \end{aligned} \] que conduce al siguiente estimador insesgado y consistente de \(\sigma^2\):

\[ \widehat{\sigma}^2 = \dfrac{1}{N-(k+1)} \left( \mathbf{Y} - \mathbf{X} \widehat{\boldsymbol{\beta}}_{GLS}\right)^\top \mathbf{\Omega}^{-1} \left( \mathbf{Y} - \mathbf{X}\widehat{\boldsymbol{\beta}}_{GLS}\right) \]

Mínimos cuadrados generalizados (MCG)

• Debido a nuestro modelo RLMG convenientemente escrito \(\mathbf{Y}^* = \mathbf{X}^* \boldsymbol{\beta} + \boldsymbol{\varepsilon}^*\). Las siguientes propiedades del MCG pueden derivarse siguiendo los mismos principios que en el caso de MCO:

• El MCG es un estimador insesgado: \[ \mathbb{E} \left( \widehat{\boldsymbol{\beta}}_{MCG} \right) = \boldsymbol{\beta} + \left( \mathbf{X}^{*\top} \mathbf{X}^* \right)^{-1} \mathbf{X}^{*\top} \mathbb{E} \left(\boldsymbol{\varepsilon}^* \right) = \boldsymbol{\beta} \]

• La matriz de varianza-covarianza del MCG es: \[ \begin{aligned} \mathbb{V}{\rm ar} \left( \widehat{\boldsymbol{\beta}}_{MCG} \right) = \sigma^2 \left( \mathbf{X}^{*\top} \mathbf{X}^* \right)^{-1} = \sigma^2 \left( \mathbf{X}^\top \mathbf{\Psi} \mathbf{\Psi}^\top\mathbf{X} \right)^{-1} = \sigma^2 \left( \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{X} \right)^{-1} \end{aligned} \]

Mínimos cuadrados generalizados (MCG)

• Si los errores se distribuyen normalmente, entonces: \[ \widehat{\boldsymbol{\beta}}_{MCG} | \mathbf{X} \sim \mathcal{N} \left(\boldsymbol{\beta},\ \sigma^2 \left( \mathbf{X}^\top \mathbf{\Omega}^{-1} \mathbf{X} \right)^{-1}\right) \] \(\widehat{\sigma}^2\) es insesgado y consistente.

• El estimador MCG es BLUE para el RLMG. En consecuencia, para el RLMG, el estimador de MCO es ineficiente.

Mínimos cuadrados ponderados (MCP)

• Es particularmente sencillo obtener estimaciones de MCG, cuando los errores son heterocedásticos pero no correlacionados, ya que esto implica que la matriz \(\mathbf{\Omega}\) sea diagonal, con elementos diagonales no constantes: \[ \mathbf{\Omega} = \begin{bmatrix} \omega_1^2 & 0 & 0 &... & 0\\ 0 & \omega_2^2 & 0 & ... &0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & ... & \omega_N^2 \end{bmatrix} \iff \mathbf{\Omega}^{-1} = \begin{bmatrix} \omega_1^{-2} & 0 & 0 &... & 0\\ 0 & \omega_2^{-2} & 0 & ... &0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & ... & \omega_N^{-2} \end{bmatrix} \]

• Entonces, si seleccionamos matriz \(\mathbf{\Psi}\) como: \[ \mathbf{\Psi} = \begin{bmatrix} \omega_1^{-1} & 0 & 0 &... & 0\\ 0 & \omega_2^{-1} & 0 & ... &0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & ... & \omega_N^{-1} \end{bmatrix} \]

Mínimos cuadrados ponderados (MCP)

• Entonces, nuestra típica expresión de regresión \(\mathbf{\Psi}^\top \mathbf{Y} = \mathbf{\Psi}^\top\mathbf{X} \boldsymbol{\beta} + \mathbf{\Psi}^\top\boldsymbol{\varepsilon}\) se puede escribir como:

• En otras palabras, todas las variables, incluido el término constante, se multiplican por los mismos pesos \(\omega_i^{-1}\).

Mínimos cuadrados ponderados (MCP)

• Hay varias formas de determinar los pesos utilizados en la estimación de mínimos cuadrados ponderados. En el caso más simple, podemos suponer que \(\mathbb{E}(\epsilon_i^2)\) es proporcional a \(Z_i^2\), donde \(Z_i\) es alguna variable observada.
• Por ejemplo \(Z_i\) puede ser la población o los ingresos.
• Alternativamente, a veces \(\mathbb{E}(\epsilon_i^2)\) puede ser proporcional al tamaño de la muestra utilizada para obtener un valor promedio (o total) de observación \(i\), que se guarda en el conjunto de datos.

Mínimos cuadrados ponderados (MCP)

• Si la observación \(Y_i\) es un promedio de \(N_i\) observaciones igualmente variables (no correlacionadas), entonces \(\mathbb{V}{\rm ar} \left( Y_i \right) = \sigma^2 / N_i\) y \(\omega_i = 1/\sqrt{N_i}\).
• Si la varianza de \(Y_i\) es proporcional a algún predictor \(Z_i\), entonces \(\mathbb{V}{\rm ar} \left( Y_i \right) = Z_i^2\sigma^2\) y \(\omega_i = Z_i\).

Mínimos cuadrados generalizados factibles (MCGF)

• En la práctica, la verdadera forma de \(\mathbf{\Omega}\) es desconocido.

• Incluso en el caso más simple de \(\mathbf{\Omega} = \text{diag} \left( \omega_1^2, ..., \omega_N^2 \right)\), necesitaríamos estimar un total de \(k+1 + N\) parámetros (\(\beta_0,...,\beta_k\) y \(\omega_1^2, ..., \omega_N^2\)), lo cual es imposible ya que solo tenemos \(N\) observaciones (y una regla empírica dice que necesitaríamos un mínimo de 5 a 20 observaciones por parámetro para obtener estimaciones satisfactorias de \(\boldsymbol{\beta}\)).

Mínimos cuadrados generalizados factibles (MCGF)

• Por lo tanto, debemos asumir algunas condiciones simplificadoras para \(\mathbf{\Omega}\). Por ejemplo, si \(\mathbf{\Omega}\) es diagonal, supongamos que \(\omega_i^2 = \alpha_0 + \alpha_1 X_{m,i}\) para algunos \(m = 1,...,k\), ahora solo necesitamos estimar dos parámetros adicionales.

• El caso, donde usamos una matriz estimada \(\widehat{\mathbf{\Omega}}\), se conoce como mínimos cuadrados generalizados factibles (o estimables) (MCGF). Los estimadores tienen buenas propiedades en muestras grandes.

Mínimos cuadrados generalizados factibles (MCGF)

• En consecuencia, utilizando la matriz estimada de \(\mathbf{\Omega}\) da como resultado el siguiente estimador MCGF:

• El estimador MCGF se define como: \[ \widehat{\boldsymbol{\beta}}_{MCGF} = \left( \mathbf{X}^\top \widehat{\mathbf{\Omega}}^{-1} \mathbf{X}\right)^{-1} \mathbf{X}^\top \widehat{\mathbf{\Omega}}^{-1} \mathbf{Y} \] donde \(\widehat{\mathbf{\Omega}} = \widehat{\mathbf{\Omega}}(\boldsymbol{\theta})\) es una estimación paramétrica (con vector de parámetro \(\boldsymbol{\theta}\)) de la verdadera matriz desconocida \(\mathbf{\Omega}\).

Una nota sobre interpretación de coeficientes

Una nota sobre interpretación de coeficientes

Una nota sobre interpretación de coeficientes

• Las estimaciones MCG (así como MCP y MCGF) de \(\boldsymbol{\beta}\) conservan su interpretación de coeficientes en la que una unidad aumenta en una variable explicativa\(X_j\) afecta la variable dependiente \(Y\).

Errores heterocedásticos

Errores heterocedásticos

• Los estimadores de MCO permanecerán insesgados;

• El estimador de la varianza de MCO es sesgado e inconsistente;

• La estadística-\(t\) habitual de las estimaciones de MCO no tienen la distribución \(t\) de Student, incluso en muestras grandes.

Errores heterocedásticos

Errores heterocedásticos

Ejemplo

set.seed(123)

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

x1 <- seq(from = 0, to = 5, length.out = N)
x2 <- sample(seq(from = 3, to = 17, length.out = 80), size = N, replace = TRUE)
e  <- rnorm(mean = 0, sd = 1:N, n = N)

x_mat <- cbind(1, x1, x2)

y <- x_mat %*% beta_vec + e

data_mat <- data.frame(y, x1, x2)

Prueba de heterocedasticidad

• Podemos examinar la presencia de heterocedasticidad 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) -0.89066   28.42373 -0.03134  0.97503
## x1           7.21929    5.93429  1.21654  0.22524
## x2          -1.83213    2.18885 -0.83703  0.40359

Diagnóstico gráfico del residuo

• Una forma de investigar la existencia de heterocedasticidad es examinar visualmente los residuos del modelo MCO.

• Si son homocedásticos, no debe haber patrón en los residuos. Si los errores son heterocedásticos, mostrarían una variación creciente o decreciente de alguna manera sistemática. Por ejemplo, la variación puede aumentar con valores mayores de \(\widehat {Y}\), o con valores mayores de \(X_ {j}\).

Diagnóstico gráfico del residuo

par(mfrow = c(2, 2))

plot(mdl_1$residuals, main = "Run-Sequence Plot")


plot(mdl_1$fitted.values, mdl_1$residuals, main = "Residuals vs Fitted")

plot(x1, mdl_1$residuals, main = bquote("Residuals vs"~X[1]))

plot(x2, mdl_1$residuals, main = bquote("Residuals vs"~X[2]))

Diagnóstico gráfico del residuo

• Observamos que los residuos parecen tener una varianza diferente, dependiendo del valor de \(\widehat {Y}\) y \(X_1\).

Pruebas de heterocedasticidad

La prueba Goldfeld-Quandt

La prueba Goldfeld-Quandt

La prueba Goldfeld-Quandt

La prueba Goldfeld-Quandt

La prueba Goldfeld-Quandt

• En nuestro modelo de ejemplo, vemos que los residuos se pueden ordenar por los valores ajustados o por \(X_1\). Si especificamos ordenar por \(X_1\):

GQ_t <- lmtest::gqtest(mdl_1, alternative = "two.sided", order.by = ~ x1)
print(GQ_t)

## 
##  Goldfeld-Quandt test
## 
## data:  mdl_1
## GQ = 10.217, df1 = 97, df2 = 97, p-value < 2.2e-16
## alternative hypothesis: variance changes from segment 1 to 2

La prueba Goldfeld-Quandt

GQ_t <- lmtest::gqtest(mdl_1, alternative = "two.sided", order.by = order(mdl_1$fitted.values))
print(GQ_t)

## 
##  Goldfeld-Quandt test
## 
## data:  mdl_1
## GQ = 5.7826, df1 = 97, df2 = 97, p-value = 3.736e-16
## alternative hypothesis: variance changes from segment 1 to 2

La prueba Goldfeld-Quandt

• Por otro lado, si hubiéramos especificado ordenar por \(X_2\):

GQ_t <- lmtest::gqtest(mdl_1, alternative = "two.sided", order.by = ~ x2)
print(GQ_t)

## 
##  Goldfeld-Quandt test
## 
## data:  mdl_1
## GQ = 1.1747, df1 = 97, df2 = 97, p-value = 0.4292
## alternative hypothesis: variance changes from segment 1 to 2

La prueba Goldfeld-Quandt

Una prueba de heterocedasticidad general

• La siguiente prueba se utiliza para la heterocedasticidad condicional, que está relacionada con las variables explicativas. Esta prueba es una generalización de la prueba de Breusch-Pagan.

Una prueba de heterocedasticidad general

Una prueba de heterocedasticidad general

Una prueba de heterocedasticidad general

Una prueba de heterocedasticidad general

Una prueba de heterocedasticidad general

• Para nuestros datos, los resultados de la prueba son:

BP_t <- lmtest::bptest(mdl_1)
print(BP_t)

## 
##  studentized Breusch-Pagan test
## 
## data:  mdl_1
## BP = 35.896, df = 2, p-value = 1.604e-08

• El \(p\)-valor es menor que el nivel de significancia de 0.05 elegido, por lo que los residuos son heterocedásticos.

La prueba White

La prueba White

La prueba White en R

• La prueba también se realiza de manera similar a como fue para la regresión univariante con una variable explicativa:

#
W_t <- lmtest::bptest(mdl_1, ~ x1 + I(x1^2) + x2 + I(x2^2) + x1:x2)
print(W_t)

## 
##  studentized Breusch-Pagan test
## 
## data:  mdl_1
## BP = 43.195, df = 5, p-value = 3.373e-08

Errores estándar consistentes con heterocedasticidad (HCE)

Errores estándar consistentes con heterocedasticidad (HCE)

Errores estándar consistentes con heterocedasticidad (HCE)

Errores estándar consistentes con heterocedasticidad (HCE)

Ejemplo en R

• En nuestro caso simulado, la corrección de White se puede estimar manualmente:

xtx_inv <- solve(t(x_mat) %*% x_mat)

V_HC <- xtx_inv %*% t(x_mat) %*% diag(mdl_1$residuals^2) %*% x_mat %*% xtx_inv

print(V_HC)

##                      x1         x2
##    693.25117 -91.740409 -52.244062
## x1 -91.74041  46.570134   2.341761
## x2 -52.24406   2.341761   4.749053

Ejemplo en R

V_HC_1 <- sandwich::vcovHC(mdl_1, type = "HC0")
print(V_HC_1)

##             (Intercept)         x1         x2
## (Intercept)   693.25117 -91.740409 -52.244062
## x1            -91.74041  46.570134   2.341761
## x2            -52.24406   2.341761   4.749053

Ejemplo en R

print(lmtest::coeftest(mdl_1, vcov. = V_HC_1))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.32966 -0.0338   0.9730
## x1           7.21929    6.82423  1.0579   0.2914
## x2          -1.83213    2.17923 -0.8407   0.4015

en comparación con los errores estándar sesgados:

print(coef(summary(mdl_1)))

##               Estimate Std. Error     t value  Pr(>|t|)
## (Intercept) -0.8906595  28.423734 -0.03133506 0.9750341
## x1           7.2192910   5.934290  1.21653830 0.2252355
## x2          -1.8321252   2.188846 -0.83702804 0.4035910

Ejemplo en R

Ejemplo en R

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

• Un estimador de HC alternativo ajusta los grados de libertad de \(\text{HC0}_{\widehat{\boldsymbol{\beta}}_{MCO}}\):

\[\text{HC1}_{\widehat{\boldsymbol{\beta}}_{MCO}} = \dfrac{N}{N-(k+1)}\left( \mathbf{X}^\top \mathbf{X}\right)^{-1} \mathbf{X}^\top \text{diag}(\widehat{\epsilon}_1^2,...,\widehat{\epsilon}_N^2)\mathbf{X} \left( \mathbf{X}^\top \mathbf{X}\right)^{-1} \]

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Especificaciones alternativas de HCE

Ejemplo en R

print(lmtest::coeftest(mdl_1, vcov. = sandwich::vcovHC(mdl_1, type = "HC0")))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.32966 -0.0338   0.9730
## x1           7.21929    6.82423  1.0579   0.2914
## x2          -1.83213    2.17923 -0.8407   0.4015

print(lmtest::coeftest(mdl_1, vcov. = sandwich::vcovHC(mdl_1, type = "HC1")))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.52939 -0.0336   0.9733
## x1           7.21929    6.87600  1.0499   0.2950
## x2          -1.83213    2.19576 -0.8344   0.4051

Ejemplo en R

print(lmtest::coeftest(mdl_1, vcov. = sandwich::vcovHC(mdl_1, type = "HC2")))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.59922 -0.0335   0.9733
## x1           7.21929    6.89407  1.0472   0.2963
## x2          -1.83213    2.20143 -0.8322   0.4063

print(lmtest::coeftest(mdl_1, vcov. = sandwich::vcovHC(mdl_1, type = "HC3")))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.87208 -0.0331   0.9736
## x1           7.21929    6.96475  1.0365   0.3012
## x2          -1.83213    2.22390 -0.8238   0.4110

Ejemplo en R

print(lmtest::coeftest(mdl_1, vcov. = sandwich::vcovHC(mdl_1, type = "HC4")))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.89066   26.72908 -0.0333   0.9735
## x1           7.21929    6.92585  1.0424   0.2985
## x2          -1.83213    2.21168 -0.8284   0.4085

Mínimos cuadrados ponderados (MCP)

• Después de detectar heterocedasticidad en los errores, es posible que deseemos imponer una estructura en la matriz de covarianza residual y estimar los coeficientes mediante MCG.

• Si sabemos que no existe una autocorrelación en los errores, entonces la covarianza es diagonal. Esto conduce a un caso específico de mínimos cuadrados ponderados por MCG (MCP).

Mínimos cuadrados ponderados (MCP)

• En realidad, no conocemos la verdadera forma de heterocedasticidad. Entonces, incluso sabiendo que la matriz de covarianza es diagonal todavía no dice nada sobre los elementos diagonales, podrían ser cualquiera de los siguientes:

\[ \begin{aligned} \mathbf{\Sigma} &= \sigma^2 \cdot \text{diag}\left(X_{j,1}, X_{j,2},..., X_{j,N} \right) \quad \text{es decir, la varianza es proporcional a } X_j\\ \mathbf{\Sigma} &= \sigma^2 \cdot \text{diag}\left(X_{j,1}^2, X_{j,2}^2,..., X_{j,N}^2 \right)\\ \mathbf{\Sigma} &= \sigma^2 \cdot \text{diag}\left(\sqrt{X_{j,1}}, \sqrt{X_{j,2}},..., \sqrt{X_{j,N}} \right), \text{ si } X_{j,i} \geq 0,\ \forall i = 1,...,N \end{aligned} \]

Mínimos cuadrados ponderados (MCP)

• Además, en un entorno de regresión múltiple, el patrón de heterocedasticidad puede depender de más de una variable explicativa; incluso podría estar relacionado con variables, no incluidas en el modelo. Entonces, ¿cómo seleccionamos la forma más probable de heterocedasticidad?

Mínimos cuadrados ponderados (MCP)

• En general, una especificación de heterocedasticidad, que funciona bastante bien, es:

Mínimos cuadrados ponderados (MCP)

Procedimiento MCG factible

Procedimiento MCG factible:

Procedimiento MCG factible:

Ejemplo con R

• Comenzamos transformando manualmente los datos y estimando MCO en las variables transformadas:

resid_data <- data.frame(log_e2 = log(mdl_1$residuals^2), x1, x2)
#
resid_mdl  <- lm(log_e2 ~ x1 + x2, data = resid_data)
h_est <- exp(resid_mdl$fitted.values)
#
data_weighted = data.frame(data_mat / sqrt(h_est), 
                           weighted_intercept = 1 / sqrt(h_est))
mdl_w_ols <- lm(y ~ -1 + weighted_intercept + x1 + x2, data = data_weighted)
print(round(coef(summary(mdl_w_ols)), 5))

##                    Estimate Std. Error  t value Pr(>|t|)
## weighted_intercept -1.39237   10.19482 -0.13658  0.89151
## x1                  7.97645    4.00325  1.99249  0.04770
## x2                 -2.02459    0.88157 -2.29656  0.02270

Ejemplo con R

• A continuación, realizamos MCP (donde \(\mathbf{\Omega}^{-1} = \text{diag} \left(\widehat{h}_1^{-1},...,\widehat{h}_N^{-1} \right)\)).

• Manualmente:

beta_wls <- solve(t(x_mat) %*% diag(1 / h_est) %*% x_mat) %*% t(x_mat) %*% diag(1 / h_est) %*% y
#
resid_wls   <- y - x_mat %*% beta_wls
sigma2_wls  <- (t(resid_wls) %*% diag(1 / h_est) %*% resid_wls) / (N - length(beta_wls))
beta_wls_se <- c(sigma2_wls) * solve(t(x_mat) %*% diag(1 / h_est) %*% x_mat)
#
print(data.frame(est = beta_wls, se = sqrt(diag(beta_wls_se))))

##          est         se
##    -1.392372 10.1948202
## x1  7.976445  4.0032528
## x2 -2.024586  0.8815747

Ejemplo con R

• Nota: en la mayoría de los programas, es necesario indicar los pesos, que son inversamente proporcionales a las varianzas.

• En otras palabras, necesita suministrar \(\text{weights}= 1 / \sigma_i^2 = 1 / \omega_i^2\) que son los elementos diagonales de \(\mathbf{\Omega}^{-1}\), no \(\mathbf{\Psi}\) (independientemente de su especificación de \(\mathbf{\Sigma}\)).

• Luego, el software toma automáticamente la raíz cuadrada de los valores especificados al calcular.

Ejemplo con R

mdl_wls <- lm(y ~ x1 + x2, data = data_mat, weights =  1 / h_est)
print(round(coef(summary(mdl_wls)), 5))

##             Estimate Std. Error  t value Pr(>|t|)
## (Intercept) -1.39237   10.19482 -0.13658  0.89151
## x1           7.97645    4.00325  1.99249  0.04770
## x2          -2.02459    0.88157 -2.29656  0.02270

Ejemplo con R

Ejemplo con R

e_star <- 1 / sqrt(h_est) * mdl_wls$residuals

par(mfrow = c(2, 2))

plot(mdl_1$fitted.values, mdl_1$residuals, main = "OLS Residuals vs Fitted")
plot(x1, mdl_1$residuals, main = bquote("OLS Residuals vs"~X[1]))
plot(mdl_wls$fitted.values, e_star, col = "blue", main = "WLS Residuals vs Fitted")
plot(x1, e_star, col = "blue", main = bquote("WLS Residuals vs"~X[1]))

Ejemplo con R

Ejemplo con R

Ejemplo con R

par(mfrow = c(2, 2))

plot(mdl_1$fitted.values, mdl_1$residuals, main = "OLS Residuals vs Fitted")
plot(x1, mdl_1$residuals, main = bquote("OLS Residuals vs"~X[1]))
plot(mdl_wls$fitted.values, mdl_wls$residuals, col = "blue", main = "WLS Residuals (of ORIGINAL data) vs Fitted")
plot(x1, mdl_wls$residuals, col = "blue", main = bquote("WLS Residuals (of ORIGINAL data) vs"~X[1]))

Ejemplo con R

Ejemplo con R

• Al observar los datos originales, sin transformar, parecería que no tomamos en cuenta la heterocedasticidad, pero este no es el caso. Como se mencionó, hemos creado un modelo sobre los datos transformados, por lo tanto, debemos analizar los residuos de los datos transformados.

Ejemplo con R

• Como se mencionó anteriormente, hemos intentado calcular aproximadamente los pesos, utilizando los logaritmos de los residuos al cuadrado. Por lo tanto, es posible que no siempre podamos capturar toda la heterocedasticidad. Sin embargo, es mucho mejor de lo que era inicialmente.

Ejemplo con R

• Por otro lado, si usáramos \(1 / i\) como los pesos, en lugar de \(1 / \sqrt{\widehat{h}_{i}}\):

mdl_wls_2 <- lm(y ~ x1 + x2, data = data_mat, weights =  1 / (1:N)^2)
print(round(coef(summary(mdl_wls_2)), 5))

##             Estimate Std. Error   t value Pr(>|t|)
## (Intercept) 15.34254    2.37257   6.46663  0.00000
## x1           5.33401    2.90696   1.83491  0.06803
## x2          -3.32553    0.22286 -14.92183  0.00000

Ejemplo con R

par(mfrow = c(1, 2))


plot(mdl_wls_2$fitted.values, 1 / 1:N * mdl_wls_2$residuals, col = "red", main = "WLS Residuals vs Fitted")
plot(x1, 1 / 1:N * mdl_wls_2$residuals, col = "red", main = bquote("WLS Residuals vs"~X[1]))

Ejemplo con R

• Los errores resultantes del modelo transformado ya no son heterocedásticos.

Ejemplo con R

• Esta es una conclusión importante con respecto a MCP (y MCGF en general): si intentamos aproximar los resultados, entonces, sin importar la aproximación, todavía existe la posibilidad de que los residuos de MCP todavía contengan algún tipo de heterocedasticidad (débil).

• Como tal, debemos examinar diferentes gráficos de residuos para determinar qué tipo de ponderaciones son apropiadas.

Ejemplo con R

• En este ejemplo, la varianza verdadera (desconocida) del término de error es \(\mathbb{V}{\rm ar}\left( \boldsymbol{\varepsilon} | \mathbf{X} \right) = \mathbf{\Sigma} = \sigma_\epsilon^2 \cdot \text{diag}(1^2, 2^2,...,N^2)\).

Respecto a R^2

• En varios programas, el coeficiente de determinación de MCP utiliza una media ponderada al calcular \(\text{TSS}\) (suma total de cuadrados). Como resultado, el reportado \(R^2\) ahora mide la proporción de variación total ponderada de \(Y\), \(Y^*\) explicado por la ponderación de X, \(X^*\).

• Para obtener una expresión más convencional de \(R^2\) se usa la expresión general (o pseudo): \[ R^2_g = \mathbb{C}{\rm orr}(Y, \widehat{Y})^2 \] donde \(\widehat{Y}\) son los valores ajustados de la variable dependiente original (es decir, no ponderada).

Ejemplo en R

• Tenga en cuenta que la media ponderada se calcula como: \[ \overline{Y}^{(w)} = \dfrac{\sum_{i = 1}^N Y_i/ h_i}{\sum_{i = 1}^N 1/ h_i} \]. Podemos verificar esto examinando los resultados manuales con el resultado:

w     <- 1 / h_est
RSS_W <- sum(w * mdl_wls$residuals^2)
TSS_W <- sum(w * (y - weighted.mean(y, w))^2)
print(1 - RSS_W / TSS_W)

## [1] 0.04576885

Ejemplo en R

print(summary(mdl_wls)$r.squared)

## [1] 0.04576885

Respecto a R^2

r2g_ols <- cor(y, mdl_1$fitted)^2
r2g_wls <- cor(y, mdl_wls$fitted)^2
print(paste0("OLS pseudo-R^2 = ", r2g_ols))

## [1] "OLS pseudo-R^2 = 0.011880153644293"

Respecto a R^2

print(paste0("WLS pseudo-R^2 = ", r2g_wls))

## [1] "WLS pseudo-R^2 = 0.0118801535929629"

Respecto a R^2

• Como se mencionó anteriormente con respecto a \(R^2\), no siempre es una buena medida de la adecuación general del modelo. Incluso si el \(R^2\) es mayor el MCO, si los residuos no se ajustan a los supuestos (S2) al (S6), entonces el modelo y los resultados de la prueba de hipótesis no son válidos.

Elegiendo entre HCE y WLS

• El estimador de mínimos cuadrados generalizados requiere que conozcamos la forma subyacente de la matriz de varianza-covarianza.

Elegiendo entre HCE y WLS

Elegiendo entre HCE y WLS

Elegiendo entre HCE y WLS

Elegiendo entre HCE y MCG

Elegiendo entre HCE y MCG

Simulación de Monte Carlo: MCO vs MCGF

Simulación de Monte Carlo: MCO vs MCGF

set.seed(123)
# Number of simulations:
MC <- 1000
# Fixed parameters
N <- 100
beta_vec <- c(10, 5, -3)
# matrix of parameter estimates for each sample:
beta_est_ols  <- NULL
beta_pval_ols <- NULL
beta_pval_hce <- NULL
#
beta_est_wls  <- NULL
beta_pval_wls <- NULL
#
beta_est_gls  <- NULL
beta_pval_gls <- NULL
#
for(i in 1:MC){
  # simulate the data:
  x1 <- seq(from = 0, to = 5, length.out = N)
  x2 <- sample(seq(from = 3, to = 17, length.out = 80), size = N, replace = TRUE)
  e  <- rnorm(mean = 0, sd = 1:N, n = N)
  y  <- cbind(1, x1, x2) %*% beta_vec + e
  data_mat <- data.frame(y, x1, x2)
  # estimate via OLS
  mdl_0 <- lm(y ~ x1 + x2, data = data_mat)
  # correct OLS se's:
  mdl_hce <- lmtest::coeftest(mdl_0, vcov. = sandwich::vcovHC(mdl_0, type = "HC0"))
  # estimate via WLS
  resid_data <- data.frame(log_e2 = log(mdl_0$residuals^2), x1, x2)
  h_est   <- exp(lm(log_e2 ~ x1 + x2, data = resid_data)$fitted.values)
  mdl_wls <- lm(y ~ x1 + x2, data = data_mat, weights =  1 / h_est)
  # estimate via GLS (by knowing the true covariance matrix)
  mdl_gls <- lm(y ~ x1 + x2, data = data_mat, weights =  1 / (1:N)^2)
  # Save the estimates from each sample:
  beta_est_ols  <- rbind(beta_est_ols, coef(summary(mdl_0))[, 1])
  beta_est_wls  <- rbind(beta_est_wls, coef(summary(mdl_wls))[, 1])
  beta_est_gls  <- rbind(beta_est_gls, coef(summary(mdl_gls))[, 1])
  # Save the coefficient p-values from each sample:
  beta_pval_ols <- rbind(beta_pval_ols, coef(summary(mdl_0))[, 4])
  beta_pval_hce <- rbind(beta_pval_hce, mdl_hce[, 4])
  beta_pval_wls <- rbind(beta_pval_wls, coef(summary(mdl_wls))[, 4])
  beta_pval_gls <- rbind(beta_pval_gls, coef(summary(mdl_gls))[, 4])
}

Simulación de Monte Carlo: MCO vs MCGF

• En primer lugar, es interesante ver cuántas veces habríamos rechazado la hipótesis nula de que un parámetro es significativo con un nivel de significancia \(\alpha = 0.05\). Dividiremos el número de veces que hubiéramos rechazado la hipótesis nula por el número total de muestras Monte Carlo para obtener la tasa de rechazo:

alpha = 0.05
a1 <- colSums(beta_pval_ols < alpha) / MC
a2 <- colSums(beta_pval_hce < alpha) / MC
a3 <- colSums(beta_pval_wls < alpha) / MC
a4 <- colSums(beta_pval_gls < alpha) / MC
#
a <- t(data.frame(a1, a2, a3, a4))
colnames(a) <- names(a1)
rownames(a) <- c("OLS: H0 rejection rate", "HCE: H0 rejection rate", "WLS: H0 rejection rate", "GLS: H0 rejection rate")
print(a)

Simulación de Monte Carlo: MCO vs MCGF

##                        (Intercept)    x1    x2
## OLS: H0 rejection rate       0.062 0.254 0.569
## HCE: H0 rejection rate       0.109 0.216 0.584
## WLS: H0 rejection rate       0.195 0.369 0.982
## GLS: H0 rejection rate       0.856 0.605 1.000

Simulación de Monte Carlo: MCO vs MCGF

• Comentaremos los resultados con R. Vemos que hubiéramos rechazado la hipótesis nula de que \(X_2\) es significativo alrededor del \(56.9\%\) de las muestras simuladas con MCO. Si tuviéramos que corregir la heterocedasticidad, esto habría aumentado a \(58.4\%\).

• Por otro lado, si tuviéramos que volver a estimar el modelo con MCP, habríamos rechazado la hipótesis nula de que \(X_2\) es significativo alrededor de \(98.2\%\) de las muestras simuladas.

Simulación de Monte Carlo: MCO vs MCGF

• Una posible explicación del rendimiento relativamente bajo de HCE es el hecho de que usa \(\widehat{\mathbf{\Sigma}} = \text{diag}(\widehat{\epsilon}_1 ^ 2, ..., \widehat{\epsilon}_N^2)\) para la matriz de covarianza. En este caso, es claramente inferior a la especificación de la matriz de covarianza utilizada en MCP. Por otro lado, como ya hemos mencionado, los HCE solo son útiles en muestras grandes.

• Debido a la gran varianza de \(\epsilon_i\), vemos que a menudo no rechazaríamos la hipótesis nula de que \(X_1\) es significativo Por otro lado, como podemos ver, usar WLS reduce este riesgo de manera significativa.

Simulación de Monte Carlo: MCO vs MCGF

• Desafortunadamente, en aplicaciones empíricas nunca conoceremos la verdadera estructura de covarianza, por lo que los resultados de GLS solo se presentan como el mejor de los casos, lo que esperamos lograr con FGLS.

Simulación de Monte Carlo: MCO vs MCGF

• Podemos mirar los histogramas de estimación de coeficientes:

# plot the estimate histograms:
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_wls[, j], plot = FALSE, breaks = 25)
  p3 <- hist(beta_est_gls[, 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, p3$counts)))
  plot(p2, col = rgb(1,0,0,1/4), add = T)
  plot(p3, col = "green", add = T)
  abline(v = beta_vec[j], lwd = 3, lty = 2, col = "red")
  #
  legend("topleft", legend = c("True Value", "OLS", "WLS", "GLS"), 
         lty = c(3, NA, NA, NA), lwd = c(3, 0, 0, 0), pch = c(NA, 22, 22, 22),
         col = c("red", "black", "black", "black"), pt.cex = 2,
         pt.bg = c(NA, rgb(0,0,1,1/4), rgb(1,0,0,1/4), "green")
         )
}

Simulación de Monte Carlo: MCO vs MCGF

Simulación de Monte Carlo: MCO vs MCGF

• Las estimaciones MCO son menos eficientes que las WLS en términos de la varianza estimada. Por otro lado, ambas estimaciones son insesgadas: su promedio es igual al valor real del parámetro. En consecuencia, dado que los estimadores de MCO son menos eficientes, necesitaríamos una muestra más grande para lograr un nivel de precisión similar al de la WLS.

• Como se mencionó anteriormente, GLS es el mejor de los casos, cuando conocemos exactamente la verdadera estructura de covarianza y no necesitamos estimarla.