pyfgls — Feasible GLS for Python

Install

pyfgls is on PyPI:

pip install pyfgls                    # core
pip install "pyfgls[formula,plot]"    # + formula interface and diagnostic plots
pip install --upgrade pyfgls          # update to the latest release

Requires Python 3.10 or newer; tested on 3.10–3.13 on Linux, macOS and Windows. Core dependencies: NumPy, SciPy, pandas, scikit-learn, joblib. Optional: formulaic (formula extra), matplotlib (plot extra), statsmodels (only used in the test suite).

A virtual environment keeps these dependencies apart from other projects:

python3 -m venv .venv
source .venv/bin/activate            # Windows: .venv\Scripts\activate
python -m pip install --upgrade pip
python -m pip install pyfgls

Check the installation; this prints the installed version:

python -c "import pyfgls; print(pyfgls.__version__)"

pip says "No matching distribution found for pyfgls"? Your pip runs on Python 3.9 or older (on macOS, the built-in python3 is 3.9), and older versions of pip do not say that the Python version is the reason. Check with python3 --version, install Python 3.10+ (from python.org, or brew install python@3.12) and use a virtual environment:

shell

python3.12 -m venv .venv && source .venv/bin/activate
python -m pip install pyfgls

Development version from GitHub:

pip install "git+https://github.com/emiliogcs/pyfgls.git"   # latest main
git clone https://github.com/emiliogcs/pyfgls.git && cd pyfgls
pip install -e ".[dev]"                                      # editable, with test tools

See CONTRIBUTING.md for the development workflow.

Quickstart

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.913171       0.876773
      std err   0.131491       0.195064
      CI width  0.515434       0.764637
x1    coef      1.412530       0.992690
      std err   0.067605       0.025818
      CI width  0.265006       0.101206

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

import pyfgls as pf
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

import pyfgls as pf
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.849550              0.716812
      std err       0.364408              0.299758
      CI width      1.428453              1.175031
x1    coef          1.001988              1.044079
      std err       0.168275              0.133527
      CI width      0.659625              0.523416

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.

4. Robust versus naive intervals

When a regressor responds to past shocks, the textbook GLS variance ([PG26] eq. 27) ignores that the AR filter was estimated and its intervals are too short. The default variance (eq. 29) accounts for it.

import pyfgls as pf
from pyfgls.simulate import serial_dgp

# AR(1) errors (rho = 0.5) and a regressor x_t = 0.5 x_{t-1} + v_t + e_{t-1}
# that responds to last period's shock: pre-determined, not exogenous
sim = serial_dgp(200, ar=(0.5,), rho_x=0.5, gamma=1.0, burn=100, random_state=3)
model = pf.SerialFGLS(sim.y, sim.X)
robust = model.fit()  # default: [PG26] eq. (29)
naive = model.fit(cov_type="naive")  # eq. (27): valid for exogenous regressors only
print(f"nonexogeneity_share = {robust.nonexogeneity_share:.3f}")
print(pf.compare([naive, robust], names=["naive", "robust"]))

nonexogeneity_share = 0.342
                   naive    robust
const coef      1.219315  1.219315
      std err   0.144750  0.146042
      CI width  0.567409  0.572474
x1    coef      0.987909  0.987909
      std err   0.052561  0.067858
      CI width  0.206037  0.266000

The estimates are identical; only the variance differs. For the slope the robust standard error is about 30% larger: the naive interval would under-cover. With exogenous regressors the share is close to 0 and the two coincide.

5. Formula interface

With the formula extra, models can be written as formulas on a pandas DataFrame; lag() creates lags, and the names and the (time) index are carried into the results.

import pandas as pd
import pyfgls as pf
from pyfgls.simulate import serial_dgp

sim = serial_dgp(240, burn=100, random_state=4)
df = pd.DataFrame(
    {"y": sim.y, "x": sim.X[:, 1]},
    index=pd.period_range("2005-01", periods=240, freq="M"),
)
res = pf.SerialFGLS.from_formula("y ~ x + lag(x)", data=df, missing="drop").fit()
print(res.summary())

==============================================================================
                           FGLS Regression Results
==============================================================================
Dep. variable:                       y  Inference:                      normal
Method:                  FGLS (serial)  AR order k*:        1  [BIC (kmax=12)]
Observations:                      239  AR coefs rho:                   0.8244
Effective obs.:                    238  Non-exog. share:                 0.064
Covariance:                     robust
------------------------------------------------------------------------------
                 coef    std err          z      P>|z|      0.025      0.975
------------------------------------------------------------------------------
const          1.6427     0.3616      4.543      0.000     0.9340     2.3514
x              0.9894     0.0413     23.971      0.000     0.9085     1.0703
lag(x)        -0.0072     0.0411     -0.175      0.861    -0.0878     0.0734
==============================================================================
Notes:
[1] Covariance: [PG26] eq. (29), robust to non-exogenous (pre-determined) regressors.
[2] Collinear lagged columns dropped from the Durbin regression: L1.const, L1.x

missing="drop" drops the first month, which lag(x) leaves empty. Note [2]: the Durbin regression already contains the lags of every regressor, so the lag of the constant and L1.x, which duplicates lag(x), are dropped there.

When to use which estimator

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).

Compared with the usual alternatives:

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.

User guide

Data requirements

• 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".

Inference options (cov_type)

Use use_t=True for Student-t critical values with n_eff − q degrees of freedom.

Reading summary()

print(res.summary()) prints three blocks: a header that describes the fit, the coefficient table, and numbered notes. res.summary(alpha=0.10) gives 90% intervals. Quick example 1 (serial correlation) and Quick example 2 (heteroskedasticity) show typical output.

Header

Coefficient table

res.to_frame() returns the whole table as a DataFrame, and res.to_dict() and res.to_json() serialise the fit.

Notes

The notes say which variance formula was used, with its equation in the paper, and report what the estimator did on its own: columns of the Durbin regression dropped because they are collinear (the lag of the constant, of trends, of seasonal dummies, or a lag that duplicates another regressor), or a declared assumption that the data contradict (exog="exogenous" with a large nonexogeneity_share). Problems found while fitting are raised as warnings at fit time (see the FAQ).

Comparing fits

Each table has one row block per coefficient with its estimate, standard error and confidence-interval width, and one column per fit.

Interpreting nonexogeneity_share

The robust variance of [PG26] eq. (29) is the naive variance (eq. (27)) plus a term that accounts for the estimation of the AR filter. That term is proportional to \(g_j = E[\tilde x_t u_{t-j}]\), which is zero when the regressors are exogenous.

res.nonexogeneity_share reports the average, over coefficients, of the correction relative to the naive variance, and res.serial.nonexogeneity_by_coef gives it per coefficient:

• ≈ 0 — the regressors behave as exogenous; naive and robust intervals coincide.
• > 0.1 — the first-stage correction matters; the naive intervals would under-cover.

The robust variance is valid in both cases, so it is always the default. Declaring exog="exogenous" does not change the estimates; it only adds a note to the summary when the share is large.

Lagged dependent variables

When lags of \(y\) are regressors (e.g. a dynamic model \(y_t = \alpha y_{t-1} + \gamma w_t + u_t\) with AR errors), OLS is inconsistent and the Durbin-regression coefficients on lagged \(y\) mix \(\alpha\) with the AR coefficients of the errors. Declare the lagged dependent variables so pyfgls can recover the AR coefficients from the lagged stochastic regressors (the Wallis ratio, [PG26] Section 4.2):

X = pd.DataFrame({"const": 1.0, "ylag": y.shift(1), "w": w}).iloc[1:]
res = pf.SerialFGLS(y.iloc[1:], X, lagged_dep_vars=["ylag"]).fit()
res.serial.wallis_regressors  # regressor(s) used for the ratio

• wallis_regressor="auto" (default) uses the regressor whose contemporaneous Durbin coefficient has the largest \(|t|\); "median" takes the median of the ratios over all eligible regressors; or name one.
• At least one stochastic regressor that is not a lag of \(y\) is required.
• If you forget to declare a regressor that equals a lag of \(y\), pyfgls warns (LaggedDependentWarning) and refuses to run the standard algorithm.

This path is less tested than the standard one; the test suite checks it by simulation (bias and coverage of \(\hat\alpha\) and \(\hat\gamma\)).

Many candidate covariates

The adaptive Lasso is designed for large candidate sets:

Z = pd.concat([df[cols], np.log(df[cols]), df[cols] ** 2], axis=1)  # your dictionary
res = pf.LassoFGLS(y, X, z=Z, random_state=0, n_jobs=-1).fit()
res.selected_z

Tips:

• Candidates are standardised internally; their units do not matter.
• z may have more columns than observations; cross-validation then becomes noisier and slower (a SampleSizeWarning tells you).
• Use pyfgls.hetero.make_dictionary(X, names, transforms=("poly2", "log", "trig")) to build the default dictionary and extend it.
• With dependent data use fold_scheme="blocked".
• lambda_rule="1se" gives sparser fits; the default "min" follows [GP24].
• Extreme weights trigger a WeightWarning; weight_clip=(0.01, 0.99) winsorises the fitted log-variance (off by default to match the paper).

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))

Reproducibility

Some steps of the estimators are random. Fix the seeds below and a fit gives the same numbers every time it is run.

• Serial correlation only (SerialFGLS, or FGLS(..., hetero=None)): the fit is deterministic. There is nothing to seed except the bootstrap.

• Heteroskedasticity (LassoFGLS, FGLS): the cross-validation that tunes the adaptive Lasso (and SVR) splits the sample into random folds. Without a seed the selected covariates, and with them the estimates, can change from run to run. Pass random_state= (an int or a numpy.random.Generator) to the model:

  python

  res = pf.LassoFGLS(y, X, random_state=0).fit()

  fold_scheme="blocked" uses contiguous blocks instead of random folds and needs no seed.

• Bootstrap (fit(cov_type="bootstrap")): pass boot_random_state=; by default it takes the model's random_state. Each replication draws from its own child seed (numpy.random.SeedSequence.spawn), so the result does not depend on n_jobs.

• Simulated data: every design in pyfgls.simulate takes random_state.

n_jobs only spreads the folds and the bootstrap replications over cores; with a fixed seed the results are the same for any value.

When you report results, record the versions as well, since numerical libraries can change the last digits:

import numpy
import pandas
import scipy
import sklearn

import pyfgls

for module in (pyfgls, numpy, scipy, pandas, sklearn):
    print(module.__name__, module.__version__)

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.

Theory

Serial correlation ([PG26])

Model: \(y_t = x_t'\beta + u_t\) with \(u_t = C(L)e_t\) invertible, so that \(A(L)u_t = e_t\) with \(A(L) = 1 - \sum_{j\ge1}\rho_j L^j\) and \(e_t\) a martingale difference sequence.

Exogenous versus pre-determined regressors

The regressors are pre-determined when \(x_t\) is uncorrelated with current and future innovations, and exogenous when, in addition, it is uncorrelated with past innovations. Many time-series regressors are only pre-determined: they react to past shocks.

Under pre-determinedness alone, OLS is generally inconsistent — its error \(u_t\) contains past innovations that are correlated with \(x_t\). GLS, which works with the innovations \(e_t\) themselves, remains consistent. This reverses the usual advice that OLS is the "robust" choice ([PG26] Section 3).

Estimating the filter

The feasible version needs an estimate of \(A(L)\) that is itself consistent under pre-determinedness. Estimating it from OLS residuals (Cochrane–Orcutt) fails because the residuals inherit the OLS inconsistency. Instead [PG26] use the Durbin regression (eq. (14))

whose regressors are all uncorrelated with \(e_t\). The order \(k\) is chosen by BIC, with every candidate order estimated on the same sample \(t = k_{\max}+1,\dots,T\) so that the criteria are comparable:

with \(\hat\sigma_k^2\) the mean squared residual of the order-\(k\) Durbin regression. Lagged deterministic terms (constant, trends, seasonal dummies) are collinear with their contemporaneous values and are dropped automatically.

FGLS and its variance

The data are quasi-differenced with \(\hat\rho_1,\dots,\hat\rho_{k^*}\) and OLS is applied (eq. (15)):

Because the filter is estimated, the variance of \(\hat\beta\) has two parts (Theorem 2): the usual GLS term and a correction for the first-stage estimation of \(\rho\), which involves \(g_j = E[\tilde x_t u_{t-j}]\). Under exogeneity \(g_j = 0\) and the correction vanishes (the naive variance, eq. (27)); otherwise it does not, and the robust variance of eqs. (28)–(29) — the default in pyfgls — must be used. Its relative size is reported as nonexogeneity_share.

Lagged dependent variables

If lags of \(y\) are among the regressors, the coefficients on lagged \(y\) in the Durbin regression no longer estimate \(\rho\). Declare them (lagged_dep_vars=[...]) and pyfgls recovers \(\rho_j\) from the ratio of the coefficients of a lagged and a contemporaneous stochastic regressor (Wallis, 1967; [PG26] Section 4.2).

Heteroskedasticity of unknown form ([GP24])

Model: \(y = X\beta + u\), \(E[u\mid X] = 0\), \(\operatorname{Var}(u_t) = v_t\) with an unknown variance function.

Estimating the skedastic function

GLS weights each observation by \(1/v_t\). [GP24] estimate \(\log v_t\) by regressing \(\log\max(\hat u_t^2,\delta^2)\) — the log of the squared OLS residuals, floored to avoid \(\log 0\) — on a set of candidate covariates \(z_t\) with an adaptive Lasso:

1. a ridge regression gives preliminary coefficients;
2. each candidate is penalised in inverse proportion to \(|\hat\gamma^{\text{ridge}}_j|^{\psi}\);
3. the penalty level \(\lambda\) and the exponent \(\psi\) are chosen jointly by cross-validation.

With \(\ell_t = \log\max(\hat u_t^2,\delta^2)\),

and the estimator is weighted least squares with weights \(1/\hat v_t\): \(\hat\beta = (X'\hat V^{-1}X)^{-1}X'\hat V^{-1}y\) with \(\hat V = \operatorname{diag}(\hat v_t)\).

The Lasso can handle many candidates, even more than observations, and selects a parsimonious model. If nothing is selected the weights are constant and FGLS reduces to OLS, so nothing is lost under homoskedasticity.

Inference

Because the estimated variance function is only an approximation, pyfgls uses a heteroskedasticity-robust sandwich for the weighted estimator (eq. (6)), with a finite-sample bias correction that inflates the squared residuals by the leverage of the weighted regression and by the degrees of freedom of the skedastic fit (eq. (8)). The resulting intervals have close to nominal coverage and are much shorter than OLS intervals under heteroskedasticity.

The parametric alternatives of Romano and Wolf (2017) ("wls_s1", "wls_s2") and a support-vector regression ("svr") are also available.

Both problems ([PG26] Supplement S-4)

With serially correlated errors whose innovations are heteroskedastic, the two corrections are applied in sequence:

1. F-C: the serial-correlation FGLS of [PG26];
2. the adaptive-Lasso skedastic fit of [GP24] on the F-C residuals, with candidate covariates built from the unfiltered regressors;
3. F-CH: weighted least squares on the quasi-differenced data.

The variance is the bias-corrected sandwich of [GP24] applied to the filtered regressors. pyfgls optionally adds a first-stage term for the estimated AR filter, the analogue of [PG26] eq. (29) for heteroskedastic innovations (first_stage_correction, on by default for pre-determined regressors). Under homoskedastic innovations F-CH is essentially equivalent to F-C.

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.

API reference

Generated from the docstrings of pyfgls 0.1.0.

Estimators

pyfgls.FGLS(
    y,
    X,
    *,
    serial='auto',
    hetero='lasso',
    z='auto',
    ar_order=None,
    ic='bic',
    kmax=None,
    min_lags=None,
    ic_penalty='order',
    lagged_dep_vars=None,
    wallis_regressor='auto',
    exog='predetermined',
    first_stage_correction=None,
    variance=None,
    delta=0.1,
    delta_scale=False,
    cv_folds=10,
    psi_grid=PSI_GRID,
    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',
    allow_gaps=False,
    check_stationarity=False,
)

>>> from pyfgls.simulate import combined_dgp
>>> sim = combined_dgp(300, ar=(0.5,), skedastic="x2", random_state=1)
>>> res = FGLS(sim.y, sim.X, random_state=0).fit()
>>> res.method
'F-CH (serial+hetero)'
>>> res.ar_order
1
>>> list(res.compare_stages().columns)
['F-C (serial)', 'F-CH (serial+hetero)']

from_formula(formula, data, **kwargs)

>>> import pandas as pd
>>> from pyfgls.simulate import serial_dgp
>>> sim = serial_dgp(300, random_state=0)
>>> df = pd.DataFrame({"y": sim.y, "x": sim.X[:, 1]})
>>> SerialFGLS.from_formula("y ~ x", data=df).fit().params.index.tolist()
['const', 'x']

fit(
    cov_type='robust',
    *,
    use_t=False,
    n_boot=999,
    bootstrap='wild',
    wild_dist='rademacher',
    boot_random_state=None,
)

pyfgls.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,
)

>>> from pyfgls.simulate import serial_dgp
>>> sim = serial_dgp(500, random_state=0)
>>> res = SerialFGLS(sim.y, sim.X).fit()
>>> res.ar_order
1
>>> round(float(res.params["x1"]), 2)
1.0

pyfgls.LassoFGLS(
    y,
    X,
    *,
    z='auto',
    skedastic='lasso',
    variance=None,
    delta=0.1,
    delta_scale=False,
    cv_folds=10,
    psi_grid=PSI_GRID,
    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',
)

>>> from pyfgls.simulate import hetero_dgp
>>> sim = hetero_dgp(400, skedastic="x2", random_state=0)
>>> res = LassoFGLS(sim.y, sim.X, random_state=0).fit()
>>> res.hetero.skedastic
'lasso'
>>> ols = LassoFGLS(sim.y, sim.X, skedastic="none").fit()
>>> bool(res.bse["x1"] < ols.bse["x1"])
True

Results

pyfgls.results.FGLSResults(
    params,
    cov,
    method,
    cov_type,
    nobs,
    n_eff,
    use_t,
    resid,
    resid_final,
    data,
    serial=None,
    hetero=None,
    stages=dict(),
    options=dict(),
    notes=(),
)

conf_int(alpha=0.05)

cov_params()

predict(X_new=None)

wald_test(restrictions, value=None)

compare_ols(*, cov_type=None, alpha=0.05)

compare_stages(alpha=0.05)

to_frame(alpha=0.05)

to_dict()

to_json(**kwargs)

summary(alpha=0.05)

plot_diagnostics(**kwargs)

pyfgls.results.compare(results, *, alpha=0.05, names=None)

Benchmarks

pyfgls.benchmarks.ols(
    y,
    X=None,
    *,
    cov_type='HAC',
    kernel=None,
    bandwidth=None,
    small_sample=False,
    use_t=False,
)

>>> from pyfgls.simulate import serial_dgp
>>> sim = serial_dgp(300, random_state=0)
>>> ols(sim.y, sim.X).method
'OLS'

pyfgls.benchmarks.cochrane_orcutt(
    y,
    X=None,
    *,
    ar_order=1,
    iterate=False,
    max_iter=100,
    tol=1e-08,
    use_t=False,
)

pyfgls.benchmarks.durbin(
    y,
    X=None,
    *,
    ar_order=None,
    kmax=None,
    use_t=False,
)

pyfgls.benchmarks.infeasible_gls(
    y,
    X=None,
    *,
    rho=None,
    omega=None,
    use_t=False,
)

pyfgls.benchmarks.hybrid_gls(
    y,
    X=None,
    *,
    alpha=0.05,
    lags=1,
    random_state=None,
    use_t=False,
)

Diagnostics

pyfgls.diagnostics.PredeterminedTest(
    table,
    ar_order,
    regressors,
    alpha,
    advice,
)

pyfgls.diagnostics.test_predetermined(
    y,
    X,
    max_lag=4,
    *,
    ar_order=None,
    kmax=None,
    regressors=None,
    alpha=0.05,
    scale='innovations',
)

LM_l = n (SSR_0 - SSR_1) / SSR_0  ~  chi^2(l m)   (m regressors tested),

pyfgls.diagnostics.ADFResult(statistic, lags, nobs, critical_values)

pyfgls.diagnostics.adf_test(x, *, max_lags=None)

pyfgls.diagnostics.check_stationarity(data)

pyfgls.diagnostics.weight_diagnostics(res)

pyfgls.diagnostics.ar_roots(res)

pyfgls.diagnostics.nonexogeneity_share(res)

scikit-learn wrapper

pyfgls.sklearn.FGLSRegressor(
    serial='auto',
    hetero='lasso',
    *,
    fit_intercept=True,
    ar_order=None,
    kmax=None,
    ic='bic',
    z='auto',
    cv_folds=10,
    random_state=None,
    n_jobs=None,
    cov_type='robust',
)

>>> from pyfgls.simulate import hetero_dgp
>>> sim = hetero_dgp(200, random_state=0)
>>> reg = FGLSRegressor(serial=None, random_state=0).fit(sim.X[:, 1:], sim.y)
>>> reg.coef_.shape
(1,)

fit(X, y)

predict(X)

Plots

pyfgls.plotting.acf(x, nlags=20)

pyfgls.plotting.plot_intervals(
    results,
    *,
    params=None,
    alpha=0.05,
    labels=None,
    ax=None,
)

pyfgls.plotting.plot_diagnostics(res, *, nlags=20, figsize=(10.0, 7.0))

Simulation designs

pyfgls.simulate.Simulation(y, X, beta, u, e, info=dict())

pyfgls.simulate.arma_errors(e, ar=(), ma=())

pyfgls.simulate.serial_dgp(
    T,
    *,
    ar=(0.8,),
    ma=(),
    rho_x=0.0,
    gamma=1.0,
    beta=(1.0, 1.0),
    burn=0,
    random_state=None,
)

>>> sim = serial_dgp(500, random_state=0)
>>> sim.X.shape
(500, 2)

pyfgls.simulate.hetero_dgp(
    T,
    *,
    skedastic='x2',
    beta=(1.0, 1.0),
    n_extra=0,
    extra_corr=0.0,
    random_state=None,
)

pyfgls.simulate.combined_dgp(
    T,
    *,
    ar=(0.5,),
    ma=(),
    skedastic='x2',
    beta=(1.0, 1.0),
    driver='x',
    phi_prob=0.75,
    gamma=0.0,
    burn=100,
    random_state=None,
)

pyfgls.simulate.predictive_dgp(
    T,
    *,
    theta=0.5,
    rho_x=0.0,
    gamma=0.0,
    beta=(1.0, 1.0),
    burn=100,
    random_state=None,
)

Exceptions and warnings

import warnings
from pyfgls.exceptions import PyFGLSWarning

warnings.simplefilter("ignore", PyFGLSWarning)

pyfgls.exceptions.PyFGLSWarning

pyfgls.exceptions.NonStationaryWarning

pyfgls.exceptions.CollinearityWarning

pyfgls.exceptions.ConvergenceWarning

pyfgls.exceptions.WeightWarning

pyfgls.exceptions.SampleSizeWarning

pyfgls.exceptions.SpecificationWarning

pyfgls.exceptions.LaggedDependentWarning

pyfgls.exceptions.ExperimentalWarning

pyfgls.exceptions.PyFGLSError

pyfgls.exceptions.InputError

pyfgls.exceptions.MissingDataError

pyfgls.exceptions.TimeIndexError

pyfgls.exceptions.CollinearityError

pyfgls.exceptions.InsufficientDataError

Replication

The tests/statistical suite reproduces the papers' simulations with fewer replications (pytest --runslow tests/statistical, run nightly in CI). The scripts in benchmarks/ produce full tables.

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

Runtimes on one core of an Apple M1 Pro: about 1.5 minutes for Table 1 with 10,000 replications, 3–4 minutes for the general case with 2,000, and about 20 minutes for the [GP24] tables with 2,000 replications at T = 100, 200, 400. The scripts use every core by default (--jobs -1), which divides these times roughly by the number of cores; --seed fixes the draws.

[PG26] Table 1 (T = 500, AR(1) 0.8, non-exogenous regressor)

[PG26] general case (T = 200, \(\rho_x\) = 0.8)

AR(2) errors (1.34, −0.42), exogenous regressor: MSE×100 of OLS 10.7 and FGLS 0.405 (paper: 11.5 and 0.41). Robust 95% intervals cover between 0.93 and 0.95 for AR(1) and AR(2) errors with \(\gamma \in \{0, 0.25, 0.5\}\); the naive intervals under-cover when \(\gamma \ne 0\).

[GP24] Table 1 (slope eMSE of Lasso-FGLS relative to OLS)

pyfgls (2,000 replications per cell) / paper (10,000 replications):

At T = 400 every ratio matches the paper; at T = 100–200 the pyfgls ratios are 0.03–0.06 higher. The tuning-free benchmarks agree with the paper at T = 100 (WLS-S1: x 0.92 vs 0.91, x² 0.71 vs 0.68, (log x)² 0.52 vs 0.54; infeasible GLS within 0.02), and so do the interval lengths of Table 2 (Lasso: x² 1.12 vs 1.10, (log x)² 0.32 vs 0.32, exp 1.35 vs 1.34), so the small-sample difference comes from the tuning of the adaptive Lasso (an open item). The 95% intervals from eqs. (6)+(8) cover 0.94–0.96 throughout.

[PG26] Supplement S-4: serial correlation and heteroskedasticity (T = 200)

Slope MSE × 100 and coverage of 95% intervals, 1,000 replications, candidate covariates \(W = (w, w^2, \cos w, \cos 2w)\).

Non-exogenous regressor (\(x_t + \gamma e_{t-1}\)): coverage of the F-CH interval with and without the first-stage correction is 0.942 / 0.777 (AR(1) 0.5, γ = 0.5) and 0.923 / 0.860 (AR(1) 0.8, γ = 1). The correction is therefore on by default for pre-determined regressors.

FAQ and troubleshooting

pip install pyfgls says "No matching distribution found". Why? pyfgls requires Python 3.10 or newer, and pip skips releases that do not support the Python it runs on. Older versions of pip do not say that this is the reason. On macOS, the built-in python3 is 3.9. Check with python3 --version, install Python 3.10+ (python.org or brew install python@3.12) and use a virtual environment: python3.12 -m venv .venv && source .venv/bin/activate && python -m pip install pyfgls.

Why not Cochrane–Orcutt? It estimates the AR filter from OLS residuals. When the regressors are only pre-determined, OLS is inconsistent, and so is the filter: in the [PG26] Table 1 design the Cochrane–Orcutt estimate of \(\rho = 0.8\) is biased by about −0.17 and the slope by +0.03, while FGLS based on the Durbin regression is essentially unbiased. pyfgls.benchmarks.cochrane_orcutt is available for comparison.

Why a robust variance for FGLS? Estimating the filter adds variability that the textbook GLS variance ignores whenever the regressors are correlated with past errors ([PG26] eq. (29)). For the heteroskedastic estimator, the fitted variance function is only an approximation, so a sandwich with a bias correction is used ([GP24] eq. (8)).

Isn't OLS with HAC standard errors safer? Only if the regressors are exogenous — otherwise OLS is inconsistent. Even then, GLS can be far more precise and HAC intervals tend to under-cover with persistent errors.

What if there is no serial correlation or heteroskedasticity? BIC selects \(k^* = 0\) and/or the Lasso selects no covariate; you get OLS with valid standard errors.

Can I use panel data? Not yet: a MultiIndex raises a TimeIndexError.

Which warnings and errors can I get, and what do I do about them?

How do I silence the warnings? All warnings derive from pyfgls.exceptions.PyFGLSWarning: warnings.simplefilter("ignore", pf.exceptions.PyFGLSWarning).

Citation

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 to cite the software itself. License: MIT © Emilio González-Coya. Questions and bug reports: GitHub Issues. Changes between versions: CHANGELOG.md.