Intervalos de predicción con heterocedasticidad

Estoy usando R para realizar una regresión lineal. He visto formas de calcular los intervalos de predicción, pero estos dependen de datos homoscedastic. ¿Hay alguna forma de calcular los intervalos de predicción con datos heteroscedásticos?

Dependería de la naturaleza de la heterocedasticidad. Si desea un intervalo de predicción, generalmente necesita una especificación paramétrica como:

$$y_i \sim N(\mathbf{x}_i'\beta,\sigma_i(\mathbf{x}_i,\mathbf{z}_i ))$$ es decir $y_i$ normalmente se distribuye con media $\mathbf{x}_i'\beta$y desviación estándar $\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$, donde la desviación estándar es alguna función conocida de la $\mathbf{x}_i$ o quizás algún otro conjunto de variables $\mathbf{z}_i$ , de esa manera puede estimar la desviación estándar para cada $i^{th}$ observación.

Los ejemplos de posibles funciones incluyen; $\sigma^2_i(\mathbf{x}_i)=\sigma^2x_{i,k}$ (Estudios de beneficios empresariales, un ejemplo del "Análisis econométrico" de Greene, 7ª edición, CH 9), donde $x_{i,k}$ es el $i^{th}$ observación de la $k^{th}$ variable dependiente o, si se trabaja con datos de series temporales, GARCH y / o especificaciones de volatilidad estocástica.

Puedes usar las estimaciones $\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$como los errores estándar para sus intervalos de predicción si lo desea. Voy a renunciar a un tratamiento formal aquí porque los errores de estimación en cuenta $\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$puede ser complicado pero, con una muestra suficientemente grande, ignorar el error de estimación no afecta tanto el intervalo de predicción. En resumen, no es necesario abrir esa lata de gusanos aquí. Para obtener una explicación más detallada de todo esto y más ejemplos, consulte el libro de Wooldridge "Introducción a la econometría: un enfoque moderno" , Capítulo 8.

El problema es que cuando las personas se refieren a la regresión heteroscedastica o "robusta", generalmente se refieren a la situación en la que la naturaleza precisa de la heteroscedasticidad (la función $\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$) no se conoce, en cuyo caso se utiliza un estimador White o de dos pasos . Estos ofrecen estimaciones consistentes para$var(\hat \beta)$ pero no para el $\sigma_i$y, por lo tanto, no tiene una forma natural de estimar los intervalos de predicción. Yo diría que los intervalos de predicción no son significativos en este contexto de todos modos. La idea detrás de estos estimadores tipo sándwich es estimar consistentemente el error estándar de los coeficientes,$\hat \beta$, sin la carga de ofrecer intervalos de predicción precisos para cada observación individual, lo que hace que las estimaciones sean más "robustas".

Editar:

Para ser claros, lo anterior solo considera la regresión de mínimos cuadrados. Otras formas de regresión no paramétrica, como la regresión cuantil, pueden ofrecer medios para obtener un intervalo de predicción sin especificación paramétrica del error estándar residual.

La regresión cuantil no paramétrica ofrece un enfoque muy general que permite tanto la heterocedasticidad como la no linealidad. Consulte la sección 9: http://www.econ.uiuc.edu/~roger/research/rq/vig.pdf

ACTUALIZACIÓN: Una aproximación razonable para un intervalo de predicción del 90% es el espacio entre la curva de regresión del percentil 5 y la curva de regresión del percentil 95. (Dependiendo de los detalles de la técnica de estimación de curvas y la escasez de datos, es posible que desee usar algo más como los percentiles 4 y 96 para ser "conservador"). La intuición para este tipo de intervalo de predicción no paramétrica está aquí en wikipedia .

Esta respuesta es solo un punto de partida. Se ha realizado una cantidad significativa de trabajo en intervalos de predicción de regresión cuantil . O simplemente haga intervalos de predicción de regresión no paramétricos .

Si la regresión de su respuesta en su variable explicativa es una línea recta y su varianza aumenta con la variable explicativa, se necesita un modelo de regresión ponderado con o (si su varianza no constante es más extrema) como su peso. Esto pondera su varianza por su valor x, de modo que hay una relación proporcional.

Aquí hay un código con los pesos incluidos en el modelo y la predicción. Tenga en cuenta que debe agregar los pesos tanto a su conjunto de datos original como a su nuevo conjunto de datos.

Gracias a @PopcornKing por su código original de Cálculo de intervalos de predicción a partir de datos heterocedásticos .

library(ggplot2)
dummySamples <- function(n, slope, intercept, slopeVar){
  x = runif(n)
  y = slope*x+intercept+rnorm(n, mean=0, sd=slopeVar*x)
  return(data.frame(x=x,y=y))
}

myDF <- dummySamples(20000,3,0,5)
plot(myDF$x, myDF$y)
w = 1/myDF$x**2
t = lm(y~x, data=myDF, weights=w) 
summary(t)

newdata = data.frame(x=seq(0,1,0.01))
w = 1/newdata$x**2
p1 = predict.lm(t, newdata, interval = 'prediction', weights=w)
a <- ggplot()
a <- a + geom_point(data=myDF, aes(x=x,y=y), shape=1)
a <- a + geom_abline(intercept=t$coefficients[1], slope=t$coefficients[2])         
a <- a + geom_abline(intercept=t$coefficients[1],   slope=t$coefficients[2], color='blue')  
a <- ggplot()
a <- a + geom_point(data=myDF, aes(x=x,y=y), shape=1)
a <- a + geom_abline(intercept=t$coefficients[1], slope=t$coefficients[2],  color='blue')
newdata$lwr = p1[,c("lwr")]
newdata$upr = p1[,c("upr")]
a <- a + geom_ribbon(data=newdata, aes(x=x,ymin=lwr, ymax=upr),   fill='yellow', alpha=0.3)
a