为相关误差构造 GLS 模型，并与 GLSAR 对照

原作：statsmodels developers。本文依据官方教程 Generalized Least Squares 完整翻译整理，保留原示例的步骤、代码与输出，并对容易误读的统计假设补充说明。

版本与复现范围：2026 年 10 月 5 日核验的 stable 页面显示 statsmodels 0.15.0。下文数值均来自原页面已有输出，并非本次执行结果；原表中的 2026 年 8 月 27 日是该次输出记录的日期，不能据此认定文章发表日期。

普通最小二乘回归得到残差后，如果怀疑相邻时点的误差存在相关性，可以把这种结构写成矩阵，再交给广义最小二乘（GLS）。本例使用 Longley 时间序列数据：先拟合 OLS，再从滞后残差估计一阶相关系数，最后比较显式指定协方差结构的 GLS 与 GLSAR。

读取数据并加入常数项

import numpy as np

import statsmodels.api as sm

Longley 是一个时间序列数据集。通过 statsmodels 内置数据接口读取后，先给自变量矩阵添加常数列，再查看前五行。

data = sm.datasets.longley.load()
data.exog = sm.add_constant(data.exog)
print(data.exog.head())

   const  GNPDEFL       GNP   UNEMP   ARMED       POP    YEAR
0    1.0     83.0  234289.0  2356.0  1590.0  107608.0  1947.0
1    1.0     88.5  259426.0  2325.0  1456.0  108632.0  1948.0
2    1.0     88.2  258054.0  3682.0  1616.0  109773.0  1949.0
3    1.0     89.5  284599.0  3351.0  1650.0  110929.0  1950.0
4    1.0     96.2  328975.0  2099.0  3099.0  112075.0  1951.0

原文从“假设存在异方差，而且知道异方差的形式”引入 sigma。需要区分的是，下面真正构造的 rho ** order 描述的是误差随时间距离变化的自相关结构：它的对角线全为 1，并未给每个时点设置不同的方差。这并不影响用 GLS 演示相关误差的目的，但不能将其解释为已经估计出一套时变方差模型。

从 OLS 残差估计一阶关系

首先拟合普通最小二乘回归，提取残差：

ols_resid = sm.OLS(data.endog, data.exog).fit().resid

原文接着假定误差满足如下带常数项的一阶自回归关系：

ε_i = β₀ + ρ ε_i−1 + η_i，其中创新项原文记为 η ∼ N(0, Σ²)。

把当前残差对上一期残差回归，就能估计这里的 ρ。原文将这称为带 trend 的 AR(1)；代码实际添加的是截距，并没有加入随时间变化的趋势列。np.asarray(...)[1:] 取第 2 个至最后一个残差，[:-1] 取第 1 个至倒数第 2 个残差，两组观测按相邻期对齐。

resid_fit = sm.OLS(
    np.asarray(ols_resid)[1:], sm.add_constant(np.asarray(ols_resid)[:-1])
).fit()
print(resid_fit.tvalues[1])
print(resid_fit.pvalues[1])

原页面给出的斜率 t 值和 p 值分别为：

-1.4390229839613828
0.17378444789154043

原作者明确指出，这里没有强证据表明误差遵循 AR(1) 过程。后续计算仍继续，是为了演示如何把假设转成 GLS 的输入，而不是确认经济数据的真实误差机制。取得回归斜率：

rho = resid_fit.params[1]

用 Toeplitz 矩阵表示时间距离

对于平稳的一阶自回归相关结构，越近的时点具有越强的相关性幅度；若 ρ 为负，相关系数的符号会随间隔奇偶交替。Toeplitz 矩阵可以方便地表达时点之间的距离。先看五个时点的距离矩阵：

from scipy.linalg import toeplitz

toeplitz(range(5))

array([[0, 1, 2, 3, 4],
       [1, 0, 1, 2, 3],
       [2, 1, 0, 1, 2],
       [3, 2, 1, 0, 1],
       [4, 3, 2, 1, 0]])

第 i,j 个元素就是时点距离 |i−j|。将矩阵扩展到全部残差的长度：

order = toeplitz(range(len(ols_resid)))

此时 rho ** order 是逐元素乘方，得到 sigma_ij = ρ^|i−j|。主对角线为 1，相邻时点为 ρ，隔两期为 ρ²。这个矩阵给出协方差的相对结构；在此同方差 AR(1) 假设下，整体方差尺度另行估计。原文以 sigma 这个接口参数名将其传入 GLS。

sigma = rho**order
gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)
gls_results = gls_model.fit()

编辑补注：这种平稳 AR(1) 相关矩阵要求 |ρ| < 1。把教程换成其他数据时，不能直接假定残差回归得到的系数一定满足条件，也应检查数据的时间顺序、缺失值和矩阵的数值性质。

与一阶 GLSAR 示例对照

真实的 ρ 并不已知，因此原文提出，可考虑使用可行广义最小二乘。原页面说明这部分支持仍具有实验性质，并给出一阶 GLSAR 示例：

glsar_model = sm.GLSAR(data.endog, data.exog, 1)
glsar_results = glsar_model.iterative_fit(1)
print(glsar_results.summary())

务必保留这里的参数含义：GLSAR(..., 1) 中的整数 1 表示 AR 的阶数，并不是把相关系数设成 1；iterative_fit(1) 中的 1 是最大迭代次数。不能把这一行改写成“迭代直到收敛”。

进一步静态核对 同版 GLSAR 源码 可见：整数阶数会把 AR 系数初始化为零；更新系数的循环是 range(maxiter - 1)。因此这里 maxiter=1 不会进入该更新循环。白化操作还会舍弃开头的一个观测。下方输出不能用来证明“已估计并收敛的一阶相关误差模型与显式 GLS 相同”。

以下保留原页面完整回归摘要，日期、时间与所有数值均为来源页面记录：

                           GLSAR Regression Results
==============================================================================
Dep. Variable:                 TOTEMP   R-squared:                       0.996
Model:                          GLSAR   Adj. R-squared:                  0.992
Method:                 Least Squares   F-statistic:                     295.2
Date:                Thu, 27 Aug 2026   Prob (F-statistic):           6.09e-09
Time:                        07:05:20   Log-Likelihood:                -102.04
No. Observations:                  15   AIC:                             218.1
Df Residuals:                       8   BIC:                             223.0
Df Model:                           6
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const      -3.468e+06   8.72e+05     -3.979      0.004   -5.48e+06   -1.46e+06
GNPDEFL       34.5568     84.734      0.408      0.694    -160.840     229.953
GNP          -0.0343      0.033     -1.047      0.326      -0.110       0.041
UNEMP        -1.9621      0.481     -4.083      0.004      -3.070      -0.854
ARMED        -1.0020      0.211     -4.740      0.001      -1.489      -0.515
POP          -0.0978      0.225     -0.435      0.675      -0.616       0.421
YEAR       1823.1829    445.829      4.089      0.003     795.100    2851.266
==============================================================================
Omnibus:                        1.960   Durbin-Watson:                   2.554
Prob(Omnibus):                  0.375   Jarque-Bera (JB):                1.423
Skew:                           0.713   Prob(JB):                        0.491
Kurtosis:                       2.508   Cond. No.                     4.80e+09
==============================================================================

原摘要还附有两项限制：标准误以误差协方差矩阵设定正确为前提；条件数约为 4.8e+09，非常大，可能意味着严重多重共线性或其他数值问题。高 R² 并不会消除这些限制。

逐项比较参数和标准误

原文最后分别输出两种拟合的系数和标准误：

print(gls_results.params)
print(glsar_results.params)
print(gls_results.bse)
print(glsar_results.bse)

为便于逐行比较，下面把原页面的四列结果合成一张表，数值保留来源显示的精度：

原作者指出，两组参数和标准误存在差别，可能涉及算法的数值差异，例如初始条件的处理；Longley 的样本量又很小。结合前面对 iterative_fit(1) 的核对，还必须考虑两段代码实际使用的相关系数不同。这里的结果适合学习接口与计算路径，不足以裁定哪一种模型对这些经济变量的解释更可靠。

复用示例前的静态检查

本次只阅读代码、核对官方文档及源码，没有安装 statsmodels，也没有运行任何回归。示例使用内置数据和数值运算，没有发现外部命令执行、硬编码凭证或拼接不可信输入的路径；这不等于已经完成包级安全审计或证明没有漏洞。

原代码没有验证估计的 ρ 是否有限、是否处于平稳区间。下面是编辑补充的防御性检查，应放在构造 sigma 之前。它没有出现在原教程中，也未经执行：

if not np.isfinite(rho) or abs(rho) >= 1:
    raise ValueError("rho must be finite and satisfy abs(rho) < 1")

如果另做迭代 GLSAR 实验，可以增加迭代上限并检查 converged、history 与最终 rho，但那将是另一个实验，不能继续沿用本页的数值输出。本稿保留原代码，不凭空提供所谓“修复后结果”。实际分析还需独立诊断残差、确认观测次序，并评估小样本与病态设计矩阵造成的不确定性。

来源、归属与许可

原文：statsmodels — Generalized Least Squares。补充核对：GLSAR 实现。网页版权标识：© Copyright 2009–2025, Josef Perktold, Skipper Seabold, Jonathan Taylor, statsmodels-developers。中文翻译、编辑说明及示意图：未完纪整理。

代码同时保留 statsmodels 的 BSD 三条款许可声明；完整声明保留在本文下方。软件许可与本次文章授权分别记录，不将开源软件许可自动推定为任何第三方材料的许可。

Copyright (C) 2006, Jonathan E. Taylor
All rights reserved.

Copyright (c) 2006-2008 Scipy Developers.
All rights reserved.

Copyright (c) 2009-2018 statsmodels Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
  a. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
  b. Redistributions in binary form must reproduce the above copyright
     notice, this list of conditions and the following disclaimer in the
     documentation and/or other materials provided with the distribution.
  c. Neither the name of statsmodels nor the names of its contributors
     may be used to endorse or promote products derived from this software
     without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.