# Prediction Intervals with Heteroscedasticity

I am using R to perform linear regression. I have seen ways to calculate prediction intervals, but these depend on homoscedastic data. Is there a way to calculate prediction intervals with heteroscedastic data?

It would depend on the nature of the heteroskedasticity. If you wanted a prediction interval, you usually need a parametric specification like:

$$y_i \sim N(\mathbf{x}_i'\beta,\sigma_i(\mathbf{x}_i,\mathbf{z}_i ))$$ i.e. $y_i$ is normally distributed with mean $\mathbf{x}_i'\beta$, and standard deviation $\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$, where the standard deviation is some known function of the $\mathbf{x}_i$ or perhaps some other set of variables $\mathbf{z}_i$ , that way you can estimate the standard deviation for each $i^{th}$ observation.

Examples of possible functions include; $\sigma^2_i(\mathbf{x}_i)=\sigma^2x_{i,k}$ (Studies of firm profits, an example from Greene's "Econometric Analysis" 7th edition CH 9), where $x_{i,k}$ is the $i^{th}$ observation of the $k^{th}$ dependent variable, or, if working with time series data, GARCH and/or stochastic volatility specifications.

You can use the estimates $\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$ as the standard errors for your prediction intervals if you like. I will forgo a formal treatment here because accounting for estimation errors in $\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$ can be complicated but, with a sufficiently large sample, ignoring the estimation error does not effect the prediction interval that much. In short, it is not necessary to open that can of worms here. For a more detailed explanation of all this and more examples, see Wooldridge's book "Introductory Econometrics: A Modern Approach", Ch 8.

The problem is that when people refer to heteroskedastic or "robust" regression, they are usually referring to the situation in which the precise nature of the heteroskedasticity (the function $\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$) is not known, in which case a White or two-step estimator is used. These offer consistent estimates for $var(\hat \beta)$ but not for the $\sigma_i$, and so you have no naturally way to estimate prediction intervals. I would argue that prediction intervals are not meaningful in this context anyway. The idea behind these sandwich type estimators is to consistently estimate the standard error of the coefficients, $\hat \beta$, without the burden of offering accurate prediction intervals for each individual observation, thus making the estimates more "robust".

## Edit:

Just to be clear, the above only considers least squares regression. Other forms of non-parametric regression, such as quantile regression, may offer means of obtaining a prediction interval without parametric specification of residual standard error.

Nonparametric quantile regression gives a very general approach that allows for both heteroscedasticity and nonlinearity. See section 9: http://www.econ.uiuc.edu/~roger/research/rq/vig.pdf

UPDATE: A reasonable approximation for a 90% prediction interval is the space between the 5th-percentile regression curve and the 95th-percentile regression curve. (Depending on the details of the curve estimation technique and the sparsity of the data, you might want to use something more like the 4th and 96th percentiles to be "conservative"). Intuition for this type of nonparametric prediction interval is here on wikipedia.

This answer is just a starting point. A significant amount of work has been done on quantile regression prediction intervals. Or just make nonparametric regression prediction intervals.

• True, but how do you get prediction intervals in quantile regression? – Zachary Blumenfeld Oct 2 '15 at 9:36

If the regression of your response on your explanatory variable is a straight line and your variance increases with the explanatory variable, a weighted regression model is needed with $w&space;=&space;1/x_{i}$ or $w&space;=&space;1/x_{i}^{2}$ (if your nonconstant variance is more extreme) as your weight. This weights your variance by your x value, so that there's a proportional relationship.

Here's code with the weights included in the model and prediction. Notice that you need to add the weights to both your original dataset and your new dataset.

Thanks to @PopcornKing for his original code from Calculating prediction intervals from heteroscedastic data.

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