广义最小二乘法¶
[1]:
import numpy as np
import statsmodels.api as sm
Longley 数据集是一个时间序列数据集
[2]:
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 并使用它来获得 GLS 模型
首先，将从 OLS 拟合中获得残差
[3]:
ols_resid = sm.OLS(data.endog, data.exog).fit().resid
假设误差项遵循具有趋势的 AR(1) 过程
\(\epsilon_i = \beta_0 + \rho\epsilon_{i-1} + \eta_i\)
其中 \(\eta \sim N(0,\Sigma^2)\)
并且 \(\rho\) 只是残差的相关性，rho 的一致估计量是对滞后残差回归残差
[4]:
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])
-1.4390229839613828
0.17378444789154032
虽然没有强有力的证据表明误差遵循 AR(1) 过程，但我们将继续进行
[5]:
rho = resid_fit.params[1]
众所周知，AR(1) 过程意味着近邻具有更强的关系，因此可以通过使用 Toeplitz 矩阵来赋予这种结构
[6]:
from scipy.linalg import toeplitz
toeplitz(range(5))
[6]:
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]])
[7]:
order = toeplitz(range(len(ols_resid)))
因此，误差协方差结构实际上是 rho**order，它定义了自相关结构
[8]:
sigma = rho ** order
gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)
gls_results = gls_model.fit()
当然，在这种情况下，确切的 rho 是未知的，因此使用可行 GLS 可能更有意义，但目前仅有实验性支持。
可以使用具有一个滞后的 GLSAR 模型获得类似的结果
[9]:
glsar_model = sm.GLSAR(data.endog, data.exog, 1)
glsar_results = glsar_model.iterative_fit(1)
print(glsar_results.summary())
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, 03 Oct 2024 Prob (F-statistic): 6.09e-09
Time: 15:44:52 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
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.8e+09. This might indicate that there are
strong multicollinearity or other numerical problems.
/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/scipy/stats/_axis_nan_policy.py:418: UserWarning: `kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=15 observations were given.
return hypotest_fun_in(*args, **kwds)
比较 gls 和 glsar 结果，可以发现参数估计值和参数估计值的标准误差之间存在一些微小的差异。 这可能是由于算法中的数值差异，例如，初始条件的处理方式，因为 Longley 数据集中的观测值数量很少。
[10]:
print(gls_results.params)
print(glsar_results.params)
print(gls_results.bse)
print(glsar_results.bse)
const -3.797855e+06
GNPDEFL -1.276565e+01
GNP -3.800132e-02
UNEMP -2.186949e+00
ARMED -1.151776e+00
POP -6.805356e-02
YEAR 1.993953e+03
dtype: float64
const -3.467961e+06
GNPDEFL 3.455678e+01
GNP -3.434101e-02
UNEMP -1.962144e+00
ARMED -1.001973e+00
POP -9.780460e-02
YEAR 1.823183e+03
dtype: float64
const 670688.699310
GNPDEFL 69.430807
GNP 0.026248
UNEMP 0.382393
ARMED 0.165253
POP 0.176428
YEAR 342.634628
dtype: float64
const 871584.051696
GNPDEFL 84.733715
GNP 0.032803
UNEMP 0.480545
ARMED 0.211384
POP 0.224774
YEAR 445.828748
dtype: float64
最后更新时间：2024 年 10 月 3 日