Logo

Generalized Least SquaresΒΆ

In [1]: import statsmodels.api as sm

In [2]: import numpy as np

In [3]: data = sm.datasets.longley.load()

In [4]: data.exog = sm.add_constant(data.exog)

The Longley dataset is a time series dataset Let’s assume that the data is heteroskedastic and that we know the nature of the heteroskedasticity. We can then define sigma and use it to give us a GLS model

First we will obtain the residuals from an OLS fit

In [5]: ols_resid = sm.OLS(data.endog, data.exog).fit().resid

Assume that the error terms follow an AR(1) process with a trend resid[i] = beta_0 + rho*resid[i-1] + e[i] where e ~ N(0,some_sigma**2) and that rho is simply the correlation of the residuals a consistent estimator for rho is to regress the residuals on the lagged residuals

In [6]: resid_fit = sm.OLS(ols_resid[1:], sm.add_constant(ols_resid[:-1])).fit()

In [7]: print resid_fit.tvalues[0]
-1.43902298398

In [8]: print resid_fit.pvalues[0]
0.173784447887

While we don’t have strong evidence that the errors follow an AR(1) process we continue

In [9]: rho = resid_fit.params[0]

As we know, an AR(1) process means that near-neighbors have a stronger relation so we can give this structure by using a toeplitz matrix

In [10]: from scipy.linalg import toeplitz

In [11]: toeplitz(range(5))
Out[11]: 
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]])
In [12]: order = toeplitz(range(len(ols_resid)))

so that our error covariance structure is actually rho**order which defines an autocorrelation structure

In [13]: sigma = rho**order

In [14]: gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)

In [15]: gls_results = gls_model.fit()

of course, the exact rho in this instance is not known so it it might make more sense to use feasible gls, which currently only has experimental support

We can use the GLSAR model with one lag, to get to a similar result

In [16]: glsar_model = sm.GLSAR(data.endog, data.exog, 1)

In [17]: glsar_results = glsar_model.iterative_fit(1)

comparing gls and glsar results, we see that there are some small differences in the parameter estimates and the resulting standard errors of the parameter estimate. This might be do to the numerical differences in the algorithm, e.g. the treatment of initial conditions, because of the small number of observations in the longley dataset.

In [18]: print gls_results.params
[     -12.76564544       -0.03800132       -2.18694871       -1.15177649
       -0.06805356     1993.95292851 -3797854.90154772]

In [19]: print glsar_results.params
[      34.55678462       -0.03434101       -1.96214395       -1.00197296
       -0.0978046      1823.1828867  -3467960.63253558]

In [20]: print gls_results.bse
[     69.43080733       0.02624768       0.38239315       0.16525269
       0.17642833     342.63462757  670688.69930775]

In [21]: print glsar_results.bse
[     84.73371452       0.03280324       0.48054486       0.21138387
       0.22477437     445.82874779  871584.05169558]

This Page