pyfgls — Feasible GLS for Python

pyfgls implements the feasible generalized least squares (FGLS) estimators of

• [PG26] Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics 258, 106346 — serially correlated errors;
• [GP24] González-Coya, E. and Perron, P. (2024). Estimation in the presence of heteroskedasticity of unknown form: a Lasso-based approach. Journal of Econometric Methods 13(1), 29–48 — heteroskedasticity of unknown form;

and their combination for time series with both problems ([PG26] Supplement S-4).

The message of both papers is the same: GLS is not a fragile luxury. Done the way they propose, it is consistent in situations where OLS is not, it is (essentially) never less precise than OLS, it is often far more precise, and its confidence intervals have accurate coverage and are shorter. pyfgls packages those procedures with statsmodels-style results, pandas support and the papers' simulation designs.

import pyfgls as pf

res = pf.FGLS(y, X).fit()  # time series: serial correlation + heteroskedasticity
print(res.summary())

Contents

1. Install
2. Which estimator should I use?
3. Guidance for practitioners
4. Quick examples
5. Full functionality and usage
6. What the estimators do
7. Replication of the papers
8. Implementation choices and open items
9. Citing

Install

pip install pyfgls                         # once released on PyPI
pip install "pyfgls[formula,plot]"         # + formula interface and diagnostic plots
pip install -e ".[dev]"                    # development install from a clone

Python ≥ 3.10. Core dependencies: NumPy, SciPy, pandas, scikit-learn, joblib. Optional: formulaic (formula extra), matplotlib (plot extra), statsmodels (only used in the test suite).

Which estimator should I use?

FGLS is the general class: serial=None turns off the serial-correlation correction and hetero=None the heteroskedasticity correction, so FGLS(y, X, hetero=None) is the same estimator as SerialFGLS(y, X) and FGLS(y, X, serial=None) the same as LassoFGLS(y, X).

Guidance for practitioners

Time-series regressions with serially correlated errors ([PG26], Section 6)

The concluding remarks of [PG26] translate into the following recommendations.

1. Use FGLS as the default. In a linear time-series regression with stationary errors, estimate by FGLS with the AR order chosen by BIC and build confidence intervals with the exogeneity-robust variance of eq. (29) — SerialFGLS(y, X).fit() does exactly this. If the regressors happen to be exogenous, the robust correction vanishes asymptotically, so nothing is lost; if they are not, it is what makes the intervals reliable. OLS, by contrast, is generally inconsistent when the regressors react to past shocks, and Cochrane–Orcutt inherits that inconsistency.
2. Check that the regressors are pre-determined. FGLS needs the regressors to be uncorrelated with current and future innovations (not with past ones). Use the variable-addition LM test of Remark 9, pf.diagnostics.test_predetermined(y, X). If it rejects, add lags of the relevant regressors to the model and refit; if the problem comes from omitted variables that are not observed, an instrumental-variables procedure is needed instead (outside the scope of this package).
3. Large negative moving-average components need care. When the errors may contain an MA term with a large negative coefficient (typical of multi-step forecast errors), a finite autoregression approximates them poorly and BIC may stop too early, which shows up as intervals that under-cover. If you know the MA structure, use a tailored transformation such as the GMA estimator (SerialFGLS(..., ma1=True)).
4. Lagged dependent variables. With lags of y among the regressors, OLS is inconsistent under serial correlation; declare them (lagged_dep_vars=["ylag"]) so the AR filter is recovered from the lagged regressors (Section 4.2).
5. Multi-step predictive regressions. Under rational expectations the h-step forecast errors are MA(h−1) and both OLS and FGLS are consistent; FGLS is more efficient as long as the implied MA process is invertible. If it may be non-invertible, FGLS is only consistent with exogenous regressors — be cautious.
6. Know your maintained assumptions. The advantages of GLS over OLS with HAC standard errors rest on stationarity, invertibility of the error process and pre-determined regressors. When these are plausible — as they usually are — the gains in precision and in the reliability of inference are substantial. Use check_stationarity=True to screen for unit roots, and watch for NonStationaryWarning (an estimated AR filter with a root near the unit circle).
7. Heteroskedastic innovations ([PG26] Remark 12). Either keep FGLS for the serial correlation and make the variance heteroskedasticity-robust (cov_type="HC3"), or — recommended — add the Lasso-based weighting of [GP24]: FGLS(y, X) (F-CH). With homoskedastic innovations F-CH is equivalent to correcting for serial correlation only, so there is no cost.

How much does nonexogeneity_share matter? It reports how large the eq. (29) correction is relative to the textbook GLS variance: near 0 the regressors behave as exogenous; large values mean the naive interval would under-cover.

Cross sections with heteroskedasticity of unknown form ([GP24], Section 6)

1. Lasso-FGLS is a safe replacement for OLS + robust standard errors. In all the designs of [GP24] it was never noticeably worse than OLS (at most a few percent of MSE when there is no heteroskedasticity, because the Lasso then selects nothing and the weights are constant), and it is markedly better when heteroskedasticity is present.
2. You do not need to know the variance function. It is enough that the candidate covariates z include at least one variable correlated with what drives the variance, plus some non-linear transforms. The Lasso picks what fits the variance pattern in the sample; it need not pick the "true" form, and it does not need to.
3. Prefer many variables and few transforms. Adding many candidate variables works well; adding many functional forms of the same variable creates near-collinearity and helps little. Trigonometric terms (cos w, cos 2w, i.e. Chebyshev-type terms) are a good flexible choice: LassoFGLS(y, X, transforms=("linear", "trig")) or build your own z with pyfgls.hetero.make_dictionary. The default dictionary {x, x², log x, (log x)², cos x, cos 2x} per regressor follows the paper's simulations.
4. Always use the bias-corrected robust variance (the default, eq. (8)). The fitted variance function is an approximation, so model-based standard errors are not reliable; the corrected sandwich gives close-to-nominal coverage with intervals that are much shorter than OLS-HC intervals under heteroskedasticity.
5. Avoid unrestricted parametric WLS when the variance driver is uncertain. The Romano–Wolf WLS estimators (skedastic="wls_s1", "wls_s2") can be much worse than OLS when the true driver is not among their covariates; SVR ("svr") is competitive in MSE but tends to give less accurate coverage. The adaptive Lasso dominates both.

Quick examples

1. Serial correlation with a non-exogenous regressor

The [PG26] Table 1 design: AR(1) errors with ρ = 0.8 and a regressor w_t = v_t + e_{t−1} that responds to last period's shock.

import pyfgls as pf
from pyfgls.simulate import serial_dgp

sim = serial_dgp(500, random_state=0)
res = pf.SerialFGLS(sim.y, sim.X).fit()
print(res.summary())
print(res.compare_ols())

==============================================================================
                           FGLS Regression Results
==============================================================================
Dep. variable:                       y  Inference:                      normal
Method:                  FGLS (serial)  AR order k*:        1  [BIC (kmax=15)]
Observations:                      500  AR coefs rho:                   0.7668
Effective obs.:                    499  Non-exog. share:                 0.015
Covariance:                     robust
------------------------------------------------------------------------------
                 coef    std err          z      P>|z|      0.025      0.975
------------------------------------------------------------------------------
const          0.8768     0.1951      4.495      0.000     0.4945     1.2591
x1             0.9927     0.0258     38.449      0.000     0.9421     1.0433
==============================================================================
Notes:
[1] Covariance: [PG26] eq. (29), robust to non-exogenous (pre-determined) regressors.
[2] Collinear lagged columns dropped from the Durbin regression: L1.const

                   OLS  FGLS (serial)
const coef      0.9132         0.8768
      std err   0.1315         0.1951
      CI width  0.5154         0.7646
x1    coef      1.4125         0.9927
      std err   0.0676         0.0258
      CI width  0.2650         0.1012

The true slope is 1: OLS is off by 0.41 (it is inconsistent here), FGLS is on target with an interval less than half as wide.

2. Heteroskedasticity of unknown form

from pyfgls.simulate import hetero_dgp

sim = hetero_dgp(400, skedastic="x2", random_state=1)  # Var(u | x) = x²
res = pf.LassoFGLS(sim.y, sim.X, random_state=0).fit()
print(res.summary())

==============================================================================
                           FGLS Regression Results
==============================================================================
Dep. variable:                       y  Inference:                      normal
Method:                  FGLS (hetero)  Skedastic fit:                   lasso
Observations:                      400  Selected z:              x1^2, log(x1)
Effective obs.:                    400  df_hat:                              3
Covariance:                     robust  Var. ratio max/min:                20.1
------------------------------------------------------------------------------
                 coef    std err          z      P>|z|      0.025      0.975
------------------------------------------------------------------------------
const          1.3531     0.2825      4.790      0.000     0.7995     1.9068
x1             0.7736     0.1455      5.317      0.000     0.4884     1.0588
==============================================================================
Notes:
[1] Covariance: [GP24] eqs. (6)+(8), sandwich with bias correction.

3. Both problems

from pyfgls.simulate import combined_dgp

sim = combined_dgp(300, ar=(0.5,), skedastic="x2", random_state=2)
res = pf.FGLS(sim.y, sim.X, random_state=0).fit()
print(res.compare_stages())  # F-C: serial correction only; F-CH: both

                F-C (serial)  F-CH (serial+hetero)
const coef            0.8496                0.7168
      std err         0.3644                0.2998
      CI width        1.4285                1.1750
x1    coef            1.0020                1.0441
      std err         0.1683                0.1335
      CI width        0.6596                0.5234

Persistent AR(2) errors, T = 200: the OLS-HAC interval is centred far from the truth and is about three times as wide as the FGLS interval.

Full functionality and usage

Input data

• y, X may be NumPy arrays, lists, pandas Series/DataFrame; names and the index are carried into the results. Include a column of ones for an intercept (a warning is issued if there is none; implicit intercepts such as a full set of dummies are recognised).
• Time order matters for serial-correlation corrections. A pandas index must be sorted and unique; a DatetimeIndex/PeriodIndex with gaps raises TimeIndexError unless allow_gaps=True (lags are then positional, e.g. trading days). Panels (MultiIndex) are not supported.
• Missing values raise MissingDataError by default; missing="drop" drops leading/trailing rows for time series (interior gaps would corrupt the lags) and any incomplete row for cross sections. Infinite values always raise.
• Rank-deficient X raises CollinearityError naming the offending columns. Lagged constants, trends and seasonal dummies in the Durbin regression are dropped automatically.
• Formulas (needs formulaic): pf.FGLS.from_formula("y ~ x1 + np.log(x2) + lag(x1)", data=df, missing="drop"); the intercept is named const.

SerialFGLS — serial correlation ([PG26])

pf.SerialFGLS(
    y,
    X,
    ar_order=None,
    ic="bic",
    kmax=None,
    min_lags=None,
    ic_penalty="order",
    lagged_dep_vars=None,
    wallis_regressor="auto",
    exog="predetermined",
    ma1=False,
    missing="raise",
    allow_gaps=False,
    check_stationarity=False,
).fit(cov_type="robust", use_t=False)

LassoFGLS — heteroskedasticity ([GP24])

pf.LassoFGLS(
    y,
    X,
    z="auto",
    skedastic="lasso",
    variance=None,
    delta=0.1,
    delta_scale=False,
    cv_folds=10,
    psi_grid=(0, 0.25, 0.5, 0.75, 1, 2),
    lambda_rule="min",
    fold_scheme="random",
    cv_weights="fold",
    weight_clip=None,
    transforms=("linear", "poly2", "log", "trig"),
    random_state=None,
    n_jobs=None,
    missing="raise",
).fit(cov_type="robust")

FGLS — both, and the general interface

pf.FGLS(y, X, serial="auto", hetero="lasso", z="auto", first_stage_correction=None, ...)

serial ∈ {"auto", "fixed", "ma1", None}, hetero ∈ {"lasso", "svr", "wls_s1", "wls_s2", "known", None}; all SerialFGLS and LassoFGLS options are accepted. The combined fit (F-CH) applies the [GP24] weights to the quasi-differenced regression, with candidate covariates built from the unfiltered regressors. first_stage_correction (default on for pre-determined regressors) adds the eq. (29)-type term for the estimated AR filter to the sandwich variance; in our simulations with non-exogenous regressors it raises the coverage of 95% intervals from 0.78 to 0.94. cov_type: "robust", "HC0", "HC3", "bootstrap".

Results (FGLSResults)

Diagnostics

from pyfgls import diagnostics

diagnostics.test_predetermined(y, X, max_lag=4)  # [PG26] Remark 9 (experimental)
diagnostics.check_stationarity(res.data)  # ADF screening (warnings)
diagnostics.weight_diagnostics(res)  # concentration of the FGLS weights

Benchmarks, bootstrap, simulation designs

pf.benchmarks.ols(y, X, cov_type="HAC")  # Andrews QS kernel, automatic bandwidth
pf.benchmarks.ols(y, X, cov_type="NW")  # Newey–West
pf.benchmarks.cochrane_orcutt(y, X)  # two-step CO (invalid if not exogenous)
pf.benchmarks.durbin(y, X)  # beta^D from the Durbin regression
pf.benchmarks.infeasible_gls(y, X, rho=[0.8])  # known AR coefficients (or omega=)
pf.benchmarks.hybrid_gls(y, X)  # [PG26] Suppl. S-1, exogenous only

res = pf.SerialFGLS(y, X).fit(cov_type="bootstrap", n_boot=999, boot_random_state=1)

from pyfgls.simulate import serial_dgp, hetero_dgp, combined_dgp, predictive_dgp

scikit-learn

from pyfgls.sklearn import FGLSRegressor
from sklearn.model_selection import TimeSeriesSplit, cross_val_score

cross_val_score(FGLSRegressor(hetero=None), X_features, y, cv=TimeSeriesSplit(5))

Warnings and errors

All warnings derive from pyfgls.exceptions.PyFGLSWarning (NonStationaryWarning, CollinearityWarning, SampleSizeWarning, SpecificationWarning, LaggedDependentWarning, WeightWarning, ConvergenceWarning, ExperimentalWarning) and all errors from PyFGLSError (InputError, MissingDataError, TimeIndexError, CollinearityError, InsufficientDataError).

Performance

python benchmarks/performance.py (T = 10,000, 10 regressors, 100 candidate covariates, one core): serial FGLS ≈ 2 s, Lasso-FGLS ≈ 3 s, combined ≈ 5 s. n_jobs parallelises the cross-validation folds and the bootstrap.

What the estimators do

Serial correlation ([PG26]). The errors are approximated by an AR(k) process. Its coefficients come from the Durbin regression of y_t on its own lags, x_t and the lags of x_t (eq. 14) — every regressor there is uncorrelated with the current innovation, so the estimates are consistent whether or not the regressors are exogenous (unlike Cochrane–Orcutt, which uses OLS residuals). The order k* minimises BIC with all orders estimated on a common sample. The data are then quasi-differenced and OLS is applied (eq. 15). The estimator is √T-consistent and asymptotically normal under pre-determinedness; with non-exogenous regressors the first stage inflates its variance by a term that eq. (29) estimates from by-products of the two steps, with no bandwidth or other tuning parameter.

Heteroskedasticity ([GP24]). log(max(û², δ²)) from OLS is regressed on the candidate covariates by an adaptive Lasso: ridge first step, penalty weights |γ̂_ridge|^(−ψ), and the pair (λ, ψ) chosen jointly by 10-fold cross-validation with common folds. The exponentiated fit gives the GLS weights. Because the variance function is only approximated, the variance of the estimator is a heteroskedasticity-robust sandwich with the Miller–Startz/Rothenberg finite-sample correction (eqs. 6 and 8), which accounts for the leverage of the weighted regression and the degrees of freedom of the Lasso fit.

Both ([PG26] Supplement S-4). F-C (the serial-correlation FGLS), then the Lasso skedastic fit on the F-C residuals, then weighted least squares on the quasi-differenced data (F-CH), with the bias-corrected sandwich on the filtered regressors.

More detail: docs/theory.

Replication of the papers

The slow tests (pytest --runslow tests/statistical, nightly in CI) re-run the papers' simulations with fewer replications. Highlights (full tables in docs/replication.md):

python benchmarks/replicate_pg26_tables.py --table 1 --reps 10000
python benchmarks/replicate_pg26_tables.py --table general --reps 2000
python benchmarks/replicate_gp24_tables.py --reps 2000 --T 100 200 400

Implementation choices and open items

Where the papers leave room, pyfgls makes the following choices (all overridable):

• kmax = int(12·(T/200)^¼); the [PG26] simulations use 12 for T = 200 and 500.
• Variances are scaled by the effective sample T − k*; normal critical values unless use_t=True.
• The Lasso degrees of freedom count the intercept; the Romano–Wolf variance uses HC3 OLS residuals (eq. 7).
• The ridge weights of the adaptive Lasso are re-estimated within each CV fold so the CV error is not flattered by weights that have seen the test fold (cv_weights="full" computes them once).
• The Lasso path is computed exactly (LARS); coordinate descent did not converge on the highly correlated default dictionary.
• The first k* observations are dropped (no Prais–Winsten).

Open items for the authors:

• Combined variance. The Supplement does not write the first-stage term for F-CH; pyfgls adds a heteroskedasticity-robust analogue of eq. (29) by default (coverage 0.94 with vs 0.78 without it for non-exogenous regressors).
• Lag penalty with many regressors. The eq. (14) BIC penalises only the AR order; with several stochastic regressors it tends to select kmax. ic_penalty="params" is available; the default follows the paper.
• Small-sample Lasso efficiency. At T = 100–200 our Lasso-FGLS eMSE ratios are 0.03–0.06 above [GP24] Table 1 (e.g. x²: 0.75–0.77 vs 0.70), although WLS-S1/S2, the infeasible GLS and the interval lengths agree with the paper and the T = 400 ratios match; the difference must come from details of the glmnet-based tuning.
• Multi-regressor Wallis ratio (wallis_regressor): "auto" uses the regressor with the largest |t|; [PG26] reports insensitivity to the choice.
• Remark 9 test. Read literally (n·R² of the Durbin residuals on X and lagged x), the test rejects correctly specified models whose regressors are pre-determined but not exogenous, because the Durbin residuals are not orthogonal to X and still contain past innovations. test_predetermined therefore uses the variable-addition form n(SSR₀ − SSR₁)/SSR₀ on the innovation scale (residuals and regressors filtered with the Durbin AR filter); simulated size 0.04, power 1.00 against an omitted x_{t−1}. scale="errors" keeps the literal reading.
• GMA. θ̂ is the exact Gaussian MLE from the residuals of the Durbin regression with k_T = int(4(T/100)^{2/9}) lags (§5), and all T observations are transformed.

Citing

If you use pyfgls, please cite the papers that the methods come from:

• Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics 258, 106346. https://doi.org/10.1016/j.jeconom.2026.106346
• González-Coya, E. and Perron, P. (2024). Estimation in the presence of heteroskedasticity of unknown form: a Lasso-based approach. Journal of Econometric Methods 13(1), 29–48. https://doi.org/10.1515/jem-2023-0007

@article{PerronGonzalezCoya2026,
  author  = {Perron, Pierre and Gonz{\'a}lez-Coya, Emilio},
  title   = {Feasible {GLS} for time series regression},
  journal = {Journal of Econometrics},
  volume  = {258},
  pages   = {106346},
  year    = {2026},
  doi     = {10.1016/j.jeconom.2026.106346}
}

@article{GonzalezCoyaPerron2024,
  author  = {Gonz{\'a}lez-Coya, Emilio and Perron, Pierre},
  title   = {Estimation in the presence of heteroskedasticity of unknown form:
             a {L}asso-based approach},
  journal = {Journal of Econometric Methods},
  volume  = {13},
  number  = {1},
  pages   = {29--48},
  year    = {2024},
  doi     = {10.1515/jem-2023-0007}
}

See also CITATION.cff. License: MIT © Emilio González-Coya.