这将取决于异方差的性质。如果您想要一个预测区间，您通常需要一个参数规范，例如：

$$y_i \sim N(\mathbf{x}_i'\beta,\sigma_i(\mathbf{x}_i,\mathbf{z}_i ))$$ IE$y_i$正态分布，均值$\mathbf{x}_i'\beta$, 和标准差$\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$，其中标准差是某个已知函数$\mathbf{x}_i$或者可能是其他一些变量集$\mathbf{z}_i$, 这样你就可以估计每个的标准差$i^{th}$观察。

可能的功能示例包括：$\sigma^2_i(\mathbf{x}_i)=\sigma^2x_{i,k}$ （企业利润研究，Greene 的“Econometric Analysis”第 7 版 CH 9 中的一个示例），其中$x_{i,k}$是个$i^{th}$的观察$k^{th}$因变量，或者，如果使用时间序列数据，GARCH 和/或随机波动率规范。

您可以使用估计$\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$如果您愿意，可以作为预测区间的标准误差。我将在这里放弃正式的处理，因为考虑到估计错误 $\hat \sigma_i(\mathbf{x}_i,\mathbf{z}_i )$可能很复杂，但是对于足够大的样本，忽略估计误差不会对预测间隔产生太大影响。简而言之，这里没有必要打开那罐虫子。有关所有这些和更多示例的更详细解释，请参阅 Wooldridge 的书“Introductory Econometrics: A Modern Approach”，第 8 章。

问题在于，当人们提到异方差或“稳健”回归时，他们通常指的是异方差的精确性质（函数$\sigma_i(\mathbf{x}_i,\mathbf{z}_i )$) 是未知的，在这种情况下使用White或两步估计器。这些提供了一致的估计$var(\hat \beta)$但不是为了$\sigma_i$，因此您自然无法估计预测区间。 我认为无论如何，预测区间在这种情况下都没有意义。这些三明治类型估计器背后的想法是一致地估计系数的标准误差，$\hat \beta$，而无需为每个单独的观察提供准确的预测区间，从而使估计更加“稳健”。

编辑：

为了清楚起见，上面只考虑了最小二乘回归。其他形式的非参数回归，例如分位数回归，可以提供在没有残差标准误差参数说明的情况下获得预测区间的方法。

非参数分位数回归提供了一种非常通用的方法，它允许异方差和非线性。见第 9 节：http: //www.econ.uiuc.edu/~roger/research/rq/vig.pdf

更新： 90% 预测区间的合理近似值是 5th-percentile 回归曲线和 95th-percentile 回归曲线之间的空间。（根据曲线估计技术的细节和数据的稀疏性，您可能希望使用更像第 4 和第 96 个百分位数的东西来“保守”）。这种类型的非参数预测区间的直觉在维基百科上。

这个答案只是一个起点。对分位数回归预测区间进行了大量工作。或者只是制作非参数回归预测区间。

如果您的响应对解释变量的回归是一条直线，并且您的方差随着解释变量的增加而增加，则需要使用或（如果您的非常数方差更极端）作为权重的加权回归模型。这会通过您的 x 值对您的方差进行加权，因此存在比例关系。

这是包含在模型和预测中的权重的代码。请注意，您需要将权重添加到原始数据集和新数据集。

感谢@PopcornKing 提供 来自异方差数据计算预测区间的原始代码。

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