Python中OLS回归HAC标准误差计算及预测置信区间相关问题

你调用np.sqrt(reg.mse_resid)得到的是残差的标准误，这个值仅由OLS拟合的残差平方和、自由度决定，和你选择的协方差估计方法（普通OLS/Newey-West HAC）完全无关，所以无论你选什么cov_type，这个输出都不会变，它不是你需要的预测相关的稳健标准误。

场景1：计算回归系数的HAC稳健标准误 #

你拟合时指定cov_type='HAC'后，直接调用reg.bse就能得到所有系数的HAC稳健标准误，这部分你已经实现了。

场景2：计算预测值的HAC稳健标准误（用于构造预测置信区间） #

预测值的稳健标准误需要结合HAC估计得到的系数协方差矩阵、预测点的解释变量取值计算，步骤如下：

1. 拟合HAC回归（你的原有代码即可）

import statsmodels.formula.api as smf
import numpy as np

reg = smf.ols('y ~ 1 + x', data=df).fit(cov_type='HAC', cov_kwds={'maxlags':1})

2. 生成预测点的设计矩阵（包含常数项）

# 若要对原数据集样本做预测
X_pred = reg.model.exog

# 若要对新数据做预测，假设新数据存在df_new中，包含x列
# X_pred = sm.add_constant(df_new['x']).values

3. 计算每个预测值的稳健标准误

# 利用HAC估计的系数协方差矩阵计算预测值方差
pred_var = np.diag(X_pred @ reg.cov_params() @ X_pred.T)
# 得到预测值的稳健标准误
pred_se = np.sqrt(pred_var)

4. 构造95%置信区间

pred_values = reg.predict(X_pred)
ci_lower = pred_values - 1.96 * pred_se
ci_upper = pred_values + 1.96 * pred_se

如果你需要的是包含扰动项波动的预测区间而非仅系数不确定性带来的预测值置信区间，则需要在上述方差基础上加残差均方，计算方式为：

# 预测区间的标准误
pred_interval_se = np.sqrt(pred_var + reg.mse_resid)

这种场景下，HAC调整仅作用于系数不确定性对应的方差部分，残差均方仍和普通OLS结果一致，最终得到的预测区间会比经典OLS的预测区间更贴合自相关数据的真实不确定性。

内容的提问来源于stack exchange，提问作者Henrik Dalriksson