C Mínimos cuadrados generalizados en R

Uno de los modelos más utilizados en estadística para el caso de datos no homogéneos es el conocido modelo clásico de regresión lineal: \[Y =\sum\limits_{j=0}^{p}X_{j}\beta_{j} +\varepsilon,\] donde \(\boldsymbol{\beta }=(\beta_{0}, \ldots,\beta_{p})^{\top}\in \mathbb{R}^{p+1}\) es un vector desconocido, \(\left\{ X_{j}:j=0, \ldots,p\right\}\) un conjunto de variables explicativas (por comodidad asumiremos que \(X_0=1\)) y \(\varepsilon\) un error aleatorio (independiente) de media cero con \(Var(\varepsilon)=\sigma^{2}\).

Supongamos que el objetivo es la estimación eficiente de los parámetros de \(\boldsymbol{\beta}\) (variación de gran escala), a partir de una muestra aleatoria simple (los datos observados) \[\left\{ \left( {X_1}_i, \ldots, {X_p}_i, Y_{i} \right) : i = 1, \ldots, n \right\}.\] Bajo las hipótesis anteriores: \[\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon},\] siendo \(\mathbf{Y}=\left( Y_1, \ldots, Y_n \right)^{\top}\), \(\mathbf{X}\) una matriz \(n\times (p+1)\) con \(\mathbf{X}_{ij} = {X_{(j-1)}}_i\) y \(\boldsymbol{\varepsilon} = \left( \varepsilon_1, \ldots,\varepsilon_n \right)^{\top}\), con \[Var\left( \mathbf{Y} \right) =Var\left( \boldsymbol{\varepsilon} \right) = \sigma^{2} \mathbf{I}_{n},\] siendo \(\mathbf{I}_{n}\) la matriz identidad \(n\times n\). En estas condiciones el estimador lineal insesgado de \(\boldsymbol{\beta}\) más eficiente (el óptimo bajo normalidad) resulta ser el estimador de mínimos cuadrados ordinarios (OLS, ordinary least squares): \[\begin{equation} \hat{\boldsymbol{\beta}}_{ols}=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Y}, \tag{C.1} \end{equation}\] que se obtiene al minimizar la correspondiente suma de cuadrados \[\hat{\boldsymbol{\beta}}_{ols} = \arg \min_{\mathbf{\beta }}\left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)^{\top }\left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)\] o al maximizar la verosimilitud asumiendo normalidad. Además, puede verse fácilmente que \[Var(\hat{\boldsymbol{\beta}}_{ols})=\sigma^{2}(\mathbf{X}^{\top}\mathbf{X})^{-1}.\]

Sin embargo la suposición de que los errores son independientes e idénticamente distribuidos influye crucialmente en la inferencia. En el modelo anterior, en lugar de errores incorrelados, si suponemos que: \[Var\left( \boldsymbol{\varepsilon} \right) =\boldsymbol{\Sigma},\] obtenemos el modelo lineal de regresión generalizado y en este caso el estimador lineal óptimo de \(\boldsymbol{\beta}\) es el estimador de mínimos cuadrados generalizados (GLS, generalized least squares; WLS, weighted least squares, como caso particular si la matriz de covarianzas es diagonal): \[\begin{equation} \hat{\boldsymbol{\beta}}_{gls} =(\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{Y}. \tag{C.2} \end{equation}\]

Podemos ver cómo se obtiene este estimador desde varios puntos de vista. Uno de ellos es aplicar una transformación lineal a los datos de forma que sean incorrelados y tenga sentido aplicar el estimador OLS. Para ello se emplea una factorización de la matriz de covarianzas (ver p.e. la Sección 7.3 de Fernández-Casal y Cao, 2021), por ejemplo la factorización de Cholesky: \[\boldsymbol{\Sigma} = L L^\top,\] donde \(L\) es una matriz triangular inferior (fácilmente invertible), por lo que: \[\begin{aligned} Var\left( L^{-1}\mathbf{Y} \right) & =Var\left( L^{-1}\boldsymbol{\varepsilon} \right) = L^{-1} \boldsymbol{\Sigma} L^{{-1}^\top} \\ &= L^{-1} LL^\top L^{\top^{-1}} = \mathbf{I}_{n}. \end{aligned}\] Aplicando esta transformación al modelo anterior: \[\begin{aligned} L^{-1}\mathbf{Y} &= L^{-1}\mathbf{X}\boldsymbol{\beta} + L^{-1}\boldsymbol{\varepsilon}\\ \tilde{\mathbf{Y}} &= \tilde{\mathbf{X}}\boldsymbol{\beta} + \tilde{\boldsymbol{\varepsilon}}, \end{aligned}\] como \(Var(\tilde{\boldsymbol{\varepsilon}}) = Var\left( L^{-1}\boldsymbol{\varepsilon} \right) = \mathbf{I}_{n}\), al minimizar la suma de cuadrados \[\left( \tilde{\mathbf{Y}}-\tilde{\mathbf{X}}\mathbf{\beta }\right)^{\top }\left( \tilde{\mathbf{Y}}-\tilde{\mathbf{X}}\mathbf{\beta }\right) \\ = \left( L^{-1}\mathbf{Y}-L^{-1}\mathbf{X}\mathbf{\beta }\right)^{\top }\left( L^{-1}\mathbf{Y}-L^{-1}\mathbf{X}\mathbf{\beta }\right) \\ = \left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)^{\top } L^{{-1}^\top}L^{-1} \left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right) \\ = \left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)^{\top } \boldsymbol{\Sigma}^{-1}\left( \mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)\] obtenemos el estimador lineal óptimo: \[\begin{aligned} \hat{\boldsymbol{\beta}}_{gls} &=(\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}})^{-1}\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{Y}} \\ &=\left(\left( L^{-1}\mathbf{X}\right)^{\top}\left( L^{-1}\mathbf{X}\right)\right)^{-1}\left( L^{-1}\mathbf{X}\right)^{\top}L^{-1}\mathbf{Y} \\ &= \left(\mathbf{X}^{\top}L^{{-1}^\top}L^{-1} \mathbf{X}\right)^{-1} \mathbf{X}^{\top}L^{{-1}^\top}L^{-1} \mathbf{Y} \\ &= (\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{Y}. \end{aligned}\]

Si \(\boldsymbol{\Sigma}=\sigma^{2} \mathbf{I}_{n}\), los estimadores (C.1) y (C.2) coinciden; pero en caso contrario las estimaciones basadas en el modelo anterior pueden llegar a ser altamente ineficientes. Puede verse fácilmente que en el caso general: \[Var\left( \hat{\boldsymbol{\beta}}_{gls} \right)=(\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1}, \\ Var\left( \hat{\boldsymbol{\beta}}_{ols} \right) =(\mathbf{X}^{\top}\mathbf{X})^{-1} (\mathbf{X}^{\top}\boldsymbol{\Sigma}\mathbf{X})(\mathbf{X}^{\top}\mathbf{X})^{-1},\] resultando además que \(Var( \hat{\boldsymbol{\beta}}_{ols}) - Var( \hat{\boldsymbol{\beta}}_{gls} )\) es una matriz semidefinida positiva (e.g. Searle, 1971, Sección 3.3).

En la práctica la matriz de covarianzas suele ser desconocida y habría que estimarla. El procedimiento habitual, para evitar la estimación de \(n(n+1)/2\) parámetros adicionales a partir del conjunto de \(n\) observaciones, suele ser asumir un modelo paramétrico adecuado (ver comentarios al final de la Sección 1.2.1).

Para más detalles, incluyendo el enfoque de máxima verosimilitud y su aplicación en la práctica para la regresión con series de tiempo, ver por ejemplo el apéndice Time-Series Regression and Generalized Least Squares in R de Fox (2019). La teoría para la estimación por máxima verosimilitud en el caso temporal (unidimensional) es totalmente análoga a la del caso espacial (multidimensional), descrita en la Sección 3.3.3 (bajo la hipótesis habitual de normalidad).