Inferencia en un modelo lineal con heteroscedasticidad condicional

Question

Supongamos que observo vectores de variables independientes $\vec{x}$ y $\vec{z}$ y la variable dependiente $y$ . Me gustaría montar un modelo de la forma: $$y = \vec{x}^{\top}\vec{\beta_1} + \sigma g\left(\vec{z}^{\top} \vec{\beta_2}\right) \epsilon,$$ donde $g$ es una función de valor positivo dos veces diferenciable, $\sigma$ es un parámetro de escala desconocido, y $\epsilon$ es una variable aleatoria gaussiana de media cero y varianza unitaria (se supone que es independiente de $\vec{x}$ y $\vec{z}$ ). Esta es esencialmente la configuración de la prueba de heteroscedasticidad de Koenker (al menos hasta donde yo entiendo).

Tengo $n$ de observaciones de $\vec{x}, \vec{z}$ y $y$ y me gustaría estimar $\vec{\beta_1}$ y $\vec{\beta_2}$ . Sin embargo, tengo algunos problemas:

1. No estoy seguro de cómo plantear el problema de estimación como algo parecido a los mínimos cuadrados (supongo que hay un truco bien conocido). Mi primera conjetura sería algo así como $$min_{\vec{\beta_1}, \vec{\beta_2}} \left(\sum_{i=1}^n \frac{\left(y_i - \vec{x_i}^{\top}\vec{\beta_1}\right)^2}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)\left(\sum_{i=1}^n \frac{1}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)^{-1},$$ pero no estoy seguro de cómo resolverlo numéricamente (tal vez un método iterativo cuasi-Newton podría servir).
2. Suponiendo que pueda plantear el problema de una manera sana y encontrar algunas estimaciones $\hat{\beta}_1, \hat{\beta}_2$ me gustaría conocer la distribución de las estimaciones para que Por ejemplo Puedo realizar pruebas de hipótesis. No me importaría probar los dos vectores de coeficientes por separado, pero preferiría alguna forma de probarlos, Por ejemplo $H_0: \vec{w_1}^{\top} \vec{\beta_1} + \vec{w_2}^{\top} \vec{\beta_2} \le c$ para un determinado $\vec{w_1}, \vec{w_2}, c$ .

Answer 1

En un contexto algo más general con $Y$ un $n$ -vector de dimensiones de $y$ -observaciones (las respuestas, o variables dependientes), $X$ un $n \times p$ matriz de $x$ -observaciones (covariables, o variables dependientes) y $\theta = (\beta_1, \beta_2, \sigma)$ los parámetros tales que $Y \sim N(X\beta_1, \Sigma(\beta_2, \sigma))$ entonces la probabilidad logarítmica negativa es $$l(\beta_1, \beta_2, \sigma) = \frac{1}{2}(Y-X\beta_1)^T \Sigma(\beta_2, \sigma)^{-1} (Y-X\beta_1) + \frac{1}{2}\log |\Sigma(\beta_2, \sigma)|$$ En la pregunta del OP, $\Sigma(\beta_2, \sigma)$ es diagonal con $$\Sigma(\beta_2, \sigma)_{ii} = \sigma^2 g(z_i^T \beta_2)^2$$ por lo que el determinante se convierte en $\sigma^{2n} \prod_{i=1}^n g(z_i^T \beta_2)^2$ y la probabilidad logarítmica resultante se convierte en $$\frac{1}{2\sigma^2} \sum_{i=1}^n \frac{(y_i-x_i^T\beta_1)^2}{ g(z_i^T \beta_2)^2} + n \log \sigma + \sum_{i=1}^n \log g(z_i^T \beta_2)$$ Hay varias formas de abordar la minimización de esta función (suponiendo que los tres parámetros son independientes de la variación).

• Se puede intentar minimizar la función mediante un algoritmo de optimización estándar recordando la restricción de que $\sigma > 0$ .
• Se puede calcular el perfil menos la probabilidad logarítmica de $(\beta_1, \beta_2)$ minimizando el exceso de $\sigma$ por el hecho de ser fijo $(\beta_1, \beta_2)$ y luego introducir la función resultante en un algoritmo estándar de optimización sin restricciones.
• Puede alternar entre la optimización de cada uno de los tres parámetros por separado. Optimizar sobre $\sigma$ puede hacerse de forma analítica, optimizando sobre $\beta_1$ es un problema de regresión por mínimos cuadrados ponderados, y la optimización sobre $\beta_2$ equivale a ajustar un modelo lineal generalizado gamma con $g^2$ el enlace inverso.

La última sugerencia me atrae porque se basa en soluciones que ya conozco bien. Además, la primera iteración es algo que consideraría hacer de todos modos. Es decir, calcular primero una estimación inicial de $\beta_1$ por mínimos cuadrados ordinarios ignorando la posible heteroscedasticidad, y luego ajustar un glm gamma a los residuos al cuadrado para obtener una estimación inicial de $\beta_2$ $-$ sólo para comprobar si el modelo más complicado parece valer la pena. Las iteraciones que incorporan la heterocedasticidad en la solución de mínimos cuadrados como pesos podrían entonces mejorar la estimación.

En cuanto a la segunda parte de la pregunta, probablemente consideraría calcular un intervalo de confianza para la combinación lineal $w_1^T\beta_1 + w_2^T\beta_2$ ya sea utilizando la asintótica estándar de MLE (comprobando con simulaciones que la asintótica funciona) o mediante bootstrap.

Editar: Por asintótica estándar de MLE Me refiero a utilizar la aproximación normal multivariante a la distribución de la MLE con matriz de covarianza la información inversa de Fisher. La información de Fisher es por definición la matriz de covarianza del gradiente de $l$ . Depende en general de los parámetros. Si puede encontrar una expresión analítica para esta cantidad, puede intentar introducir la MLE. Como alternativa, puede estimar la información de Fisher mediante el método observado La información de Fisher, que es el hessiano de $l$ en el MLE. Su parámetro de interés es una combinación lineal de los parámetros en los dos $\beta$ -por lo que a partir de la normal multivariante aproximada de la MLE se puede encontrar una aproximación normal de la distribución de los estimadores como se describe aquí . De este modo se obtiene un error estándar aproximado y se pueden calcular los intervalos de confianza. Está bien descrito en muchos libros de estadística (matemática), pero una presentación razonablemente accesible que puedo recomendar es Con toda probabilidad por Yudi Pawitan. De todos modos, la derivación formal de la teoría asintótica es bastante complicada y se basa en una serie de condiciones de regularidad, y sólo da validez a asintótica distribuciones. Por ello, en caso de duda, siempre haría algunas simulaciones con un nuevo modelo para comprobar si puedo confiar en los resultados para parámetros y tamaños de muestra realistas. Un bootstrapping simple y no paramétrico en el que se muestrean las triplas $(y_i,x_i,z_i)$ del conjunto de datos observados con reemplazo puede ser una alternativa útil si el procedimiento de ajuste no requiere demasiado tiempo.