Package {fgls}

Type:	Package
Title:	Feasible GLS for Serial Correlation and Heteroskedasticity
Version:	0.1.0
Description:	Feasible generalized least squares (FGLS) estimators for linear regressions whose errors are serially correlated, heteroskedastic of unknown form, or both. Serial correlation is handled by the autoregressive filter of Perron and González-Coya (2026) <doi:10.1016/j.jeconom.2026.106346>, estimated from a Durbin regression whose order is chosen by an information criterion, with a variance estimator that stays valid when the regressors are only pre-determined rather than strictly exogenous. Heteroskedasticity is handled by the adaptive Lasso estimate of the skedastic function of González-Coya and Perron (2024) <doi:10.1515/jem-2023-0007>, with a bias-corrected sandwich variance. Both corrections can be combined. Benchmark estimators (ordinary least squares with heteroskedasticity and autocorrelation consistent standard errors, Cochrane-Orcutt, infeasible GLS) and the simulation designs of both papers are included.
License:	MIT + file LICENSE
URL:	https://github.com/emiliogcs/fgls, https://emiliogcs.com/fgls/
BugReports:	https://github.com/emiliogcs/fgls/issues
Depends:	R (≥ 4.1.0)
Imports:	glmnet (≥ 4.1), graphics, grDevices, stats, utils
Suggests:	broom, e1071, generics, jsonlite, knitr, lmtest, parallel, rmarkdown, sandwich, spelling, testthat (≥ 3.0.0), tibble, xts, zoo
VignetteBuilder:	knitr, rmarkdown
Config/testthat/edition:	3
Encoding:	UTF-8
Language:	en-US
RoxygenNote:	7.3.1
NeedsCompilation:	no
Packaged:	2026-10-23 15:26:52 UTC; emiliogonzalez
Author:	Emilio González-Coya [aut, cre, cph]
Maintainer:	Emilio González-Coya <emiliogc@bu.edu>
Repository:	CRAN
Date/Publication:	2026-11-07 15:40:08 UTC

fgls: Feasible GLS for serial correlation and heteroskedasticity

Description

Feasible generalized least squares (FGLS) for linear regressions with serially correlated errors (Perron and González-Coya, 2026), heteroskedastic errors of unknown form (González-Coya and Perron, 2024), or both (Supplement S-4 of the 2026 paper).

Author(s)

Maintainer: Emilio González-Coya emiliogc@bu.edu (ORCID) [copyright holder]

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

See Also

Useful links:

• https://github.com/emiliogcs/fgls

• https://emiliogcs.com/fgls/

• Report bugs at https://github.com/emiliogcs/fgls/issues

Augmented Dickey-Fuller test (screening)

Description

The regression \Delta x_t = a + g x_{t-1} + \sum_j c_j \Delta x_{t-j} + e_t with the number of lags chosen by BIC (at most \lfloor 12 (T/100)^{1/4} \rfloor); returns the t statistic on g with MacKinnon's asymptotic critical values. Used by fgls(check_stationarity = TRUE) to screen the response and the regressors; it is a screening device, not a replacement for a careful unit-root analysis.

Usage

adf_test(x, max_lags = NULL)

## S3 method for class 'fgls_adf'
print(x, ...)

Arguments

Value

A list of class "fgls_adf" with statistic, lags, nobs, critical_values and unit_root_rejected (at 5%). The print method returns x invisibly.

See Also

Other diagnostics: test_predetermined(), weight_diagnostics()

Examples

set.seed(1)
adf_test(cumsum(rnorm(200)))$unit_root_rejected
adf_test(rnorm(200))$unit_root_rejected

Compare fits side by side

Description

A table of coefficients, standard errors and confidence-interval widths for several fits, e.g. FGLS against OLS with HAC standard errors (PG26) or OLS with HC3 standard errors (GP24), or the F-C and F-CH stages of a combined fit.

Usage

compare_fgls(..., level = 0.95, ols_vcov = NULL)

Arguments

Details

Plain stats::lm() fits are accepted and converted to OLS with HAC (Andrews quadratic-spectral) standard errors when any other fit corrects for serial correlation, and with HC3 standard errors otherwise; override with ols_vcov. Other models that inherit from "lm" (glm, multivariate mlm fits, weighted fits, fits with an offset, ...) are rejected, because the conversion would replace them by unweighted OLS; use ols_robust() for OLS with robust standard errors.

Value

A data frame with columns term, statistic ("estimate", "std.error", "ci.width") and one column per fit.

See Also

Other fgls: fgls(), fgls-methods, fgls_control(), lasso_fgls(), plot.fgls(), serial_fgls(), summary.fgls()

Examples

sim <- simulate_serial(300, seed = 1)
compare_fgls(
  OLS = stats::lm(y ~ x, data = sim),
  FGLS = serial_fgls(y ~ x, data = sim)
)

Feasible GLS for serial correlation and/or heteroskedasticity

Description

fgls() fits a linear regression by feasible generalized least squares, correcting for serially correlated errors (Perron and González-Coya, 2026), for heteroskedasticity of unknown form (González-Coya and Perron, 2024), or for both (the "F-CH" estimator of Supplement S-4 of the 2026 paper). serial_fgls() and lasso_fgls() are shortcuts for the two single corrections.

Usage

fgls(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  serial = c("auto", "fixed", "ma1", "none"),
  hetero = c("lasso", "svr", "wls_s1", "wls_s2", "known", "none"),
  vcov = c("robust", "naive", "nonrobust", "HC0", "HC1", "HC2", "HC3", "bootstrap"),
  z = "auto",
  ar_order = NULL,
  kmax = NULL,
  ic = c("bic", "aic"),
  lagged_dep_vars = NULL,
  wallis_regressor = "auto",
  exog = c("predetermined", "exogenous"),
  first_stage_correction = NULL,
  variance = NULL,
  delta = 0.1,
  psi_grid = c(0, 0.25, 0.5, 0.75, 1, 2),
  cv_folds = 10L,
  lambda_rule = c("min", "1se"),
  weight_clip = NULL,
  seed = NULL,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  check_stationarity = FALSE,
  control = fgls_control(),
  contrasts = NULL,
  model = TRUE,
  weights,
  ...
)

Arguments

Value

An object of class "fgls", a list with components including

coefficients, vcov

Estimates and their covariance matrix.

residuals, fitted.values

y_t - x_t'\hat\beta and x_t'\hat\beta on the full sample.

innovations

Residuals of the final (filtered and/or weighted) regression.

ar_order, rho, bic_path, nonexogeneity_share

Serial block (NULL without serial correction); more in serial.

skedastic, weights

Heteroskedasticity block: selected covariates, coefficients, \psi, \lambda, degrees of freedom, fitted variances; weights are 1/\hat v_t. For the Lasso, skedastic$lambda is the selected penalty on the scale of 'glmnet', whose objective is that of GP24 eq. (4) divided by the number of observations n of the skedastic regression, and skedastic$lambda_paper = n \lambda is the same penalty on the scale of eq. (4). Both refer to the standardized candidate covariates.

stages

Intermediate estimators (stages$serial is F-C). They are "fgls" objects with the same methods as the final fit; vcov(fit$stages$serial, type =) accepts "robust", "HC0" and "HC3".

nobs, n_eff

Number of observations and of observations in the final regression (T - k^*).

Use the methods summary(), coef(), vcov(), confint(), predict(), residuals(), plot() and compare_fgls().

Serial correlation (PG26)

The errors are approximated by an AR(k) process. Its coefficients are the coefficients on the lagged dependent variable in the Durbin regression of y_t on y_{t-1}, \dots, y_{t-k}, x_t, \dots, x_{t-k} (eq. 14), and k is chosen by BIC over 0, \dots, k_{max} on a common sample. The data are then quasi-differenced and the model is estimated by OLS on the filtered data (eq. 15), dropping the first k observations. Unlike Cochrane-Orcutt, this is consistent when the regressors are only pre-determined (uncorrelated with current and future innovations) rather than strictly exogenous. The default variance (vcov = "robust", eqs. 28–29) accounts for the estimation of the filter and is valid in both cases; vcov = "naive" (eq. 27) is valid only for exogenous regressors. The share by which the robust variance exceeds the naive one is reported as nonexogeneity_share. If the selected order is 0, FGLS equals OLS.

Heteroskedasticity (GP24)

The log of the squared OLS residuals (floored at delta^2) is regressed on candidate covariates z_t by an adaptive Lasso: a ridge first step gives weights |\hat\gamma_j|^{-\psi}, and the penalty \lambda and the exponent \psi are chosen jointly by 10-fold cross-validation with the same folds for every \psi. The exponentiated fit gives the weights of a weighted least-squares regression. Inference uses a sandwich variance with the bias correction of eq. (8), so it stays valid when the skedastic model is wrong. If the Lasso selects no covariate the weights are constant and FGLS equals OLS (no loss under homoskedasticity).

Both (Supplement S-4)

With serial != "none" and hetero != "none" the serial correction is applied first ("F-C"), the skedastic function is fitted to the residuals of the quasi-differenced regression with candidate covariates built from the unfiltered regressors, and the quasi-differenced regression is re-estimated by weighted least squares ("F-CH"). The F-C stage is kept in fit$stages$serial.

Choosing options

• Time series with possibly autocorrelated errors and no concern about heteroskedasticity: hetero = "none" (or serial_fgls()).

• Cross-sectional data: serial = "none" (or lasso_fgls()). Rows can be in any order.

• Lagged dependent variables among the regressors must be declared in lagged_dep_vars; the AR coefficients are then recovered by the Wallis ratio (PG26 Section 4.2).

• Errors known to be MA(1), e.g. multi-step forecast errors: serial = "ma1" (GMA estimator, PG26 Section 5).

• With several stochastic regressors the BIC that penalizes only the AR order can select kmax; use control = fgls_control(ic_penalty = "params").

• Results involving cross-validation depend on the random folds: set seed for reproducibility.

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

See Also

Other fgls: compare_fgls(), fgls-methods, fgls_control(), lasso_fgls(), plot.fgls(), serial_fgls(), summary.fgls()

Examples

# Serial correlation with a non-exogenous regressor (PG26 Table 1 design)
sim <- simulate_serial(300, seed = 1)
fit <- fgls(y ~ x, data = sim, hetero = "none")
summary(fit)

# Heteroskedasticity of unknown form (GP24)
het <- simulate_hetero(200, skedastic = "log2", seed = 2)
fit2 <- fgls(y ~ x, data = het, serial = "none", seed = 1)
fit2$skedastic$selected

# Both (Supplement S-4)
both <- simulate_combined(200, ar = 0.5, skedastic = "x2", seed = 3)
fit3 <- fgls(y ~ x, data = both, seed = 1)
compare_fgls(fit3$stages$serial, fit3)

Tidy, glance and augment fgls fits

Description

Methods for the 'broom' generics (from the 'generics' package), registered when 'generics' is loaded.

Usage

## S3 method for class 'fgls'
tidy(x, conf.int = FALSE, conf.level = 0.95, ...)

## S3 method for class 'fgls'
glance(x, ...)

## S3 method for class 'fgls'
augment(
  x,
  data = NULL,
  newdata = NULL,
  se_fit = FALSE,
  interval = c("none", "confidence"),
  conf.level = 0.95,
  ...
)

Arguments

Value

A data frame (a tibble if 'tibble' is installed): one row per coefficient for tidy(), one row for glance(), one row per observation for augment() (columns .fitted, optionally .lower, .upper and .se.fit, then .resid, .innovation and, with a heteroskedasticity correction, .weight). nonexogeneity.share in glance() is NA when the share was not computed.

Examples

if (requireNamespace("broom", quietly = TRUE)) {
  fit <- serial_fgls(y ~ x, data = simulate_serial(200, seed = 1))
  broom::tidy(fit, conf.int = TRUE)
  broom::glance(fit)
  head(broom::augment(fit))
}

Conditions signaled by fgls

Description

All warnings raised by the package inherit from class "fgls_warning" and all errors from "fgls_error", so they can be caught or muffled as a group with tryCatch(), withCallingHandlers() or suppressWarnings() (whose classes argument, available since R 4.0.0, filters by class).

Warning classes

• fgls_warning_nonstationary: the estimated autoregressive filter has a root on, inside or close to the unit circle.

• fgls_warning_collinearity: a matrix that should be positive definite was not numerically so (for example a nearly singular cross-product of the regressors, or a covariance matrix with a negative eigenvalue), and a pseudo-inverse or eigenvalue clipping was used instead. Collinear lag columns dropped from the Durbin regression (lags of the intercept, of a trend or of seasonal dummies) do not signal it: they are recorded in fit$serial$dropped_lags and, for serial-only fits, in the notes of summary.fgls().

• fgls_warning_convergence: an iterative routine did not fully converge.

• fgls_warning_weight: the estimated weights are extreme.

• fgls_warning_sample_size: the sample is small relative to the number of parameters, lags or folds.

• fgls_warning_specification: the specification looks unusual (e.g. no intercept, the selected order hits kmax).

• fgls_warning_lagged_dependent: a regressor looks like a lag of the dependent variable (also a specification warning).

• fgls_warning_experimental: the feature is experimental.

Error classes

fgls_error_input (invalid arguments) and its subclasses fgls_error_missing, fgls_error_time_index, fgls_error_collinearity and fgls_error_insufficient_data.

Examples

# an interior missing value in a time series is an error of class fgls_error_missing
df <- data.frame(y = c(1, 2, NA, 4, 5, 3, 2), x = c(2, 1, 3, 5, 4, 2, 1))
tryCatch(serial_fgls(y ~ x, data = df), fgls_error = function(e) class(e)[1])

# warnings can be caught by class too
withCallingHandlers(
  serial_fgls(y ~ x - 1, data = simulate_serial(100, seed = 1)),
  fgls_warning_specification = function(w) invokeRestart("muffleWarning")
)$ar_order

Methods for fgls fits

Description

Standard accessors for objects returned by fgls(), serial_fgls(), lasso_fgls() and the benchmark estimators. They make the fits work with generic tools such as lmtest::coeftest() and lmtest::coefci().

Usage

## S3 method for class 'fgls'
vcov(object, type = NULL, ...)

## S3 method for class 'fgls'
confint(object, parm, level = 0.95, type = NULL, ...)

## S3 method for class 'fgls'
residuals(object, type = c("response", "innovations"), ...)

## S3 method for class 'fgls'
fitted(object, ...)

## S3 method for class 'fgls'
predict(
  object,
  newdata,
  se.fit = FALSE,
  interval = c("none", "confidence"),
  level = 0.95,
  type = "response",
  na.action = stats::na.pass,
  ...
)

## S3 method for class 'fgls'
nobs(object, ...)

## S3 method for class 'fgls'
df.residual(object, ...)

## S3 method for class 'fgls'
formula(x, ...)

## S3 method for class 'fgls'
model.matrix(object, ...)

## S3 method for class 'fgls'
model.frame(formula, ...)

## S3 method for class 'fgls'
terms(x, ...)

## S3 method for class 'fgls'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Details

Inference is asymptotically normal (as in the papers) unless the fit was made with use_t = TRUE. Accordingly, df.residual() returns Inf for normal inference, so that lmtest::coeftest() reproduces summary.fgls(), and n_{eff} - q otherwise. nobs() returns the number of observations in the final regression, n_{eff} = T - k^*.

Value

The usual values: a coefficient covariance matrix (vcov), a matrix of confidence limits (confint), vectors of residuals or fitted values, predictions (a vector, or a list/matrix with se.fit/interval), or the corresponding model components. The print() method returns x invisibly.

Variances and the sandwich package

The variance estimators of the papers are available through vcov(fit, type = ...) and can be passed to other tools (for example lmtest::coeftest(fit, vcov. = vcov(fit, type = "naive"))). The package provides no sandwich::estfun() or sandwich::bread() methods, because the robust variance of PG26 eq. (29) and the bias-corrected sandwich of GP24 eq. (8) are not of the generic form those functions assume: sandwich::vcovHAC() and related functions would not reproduce them.

Variance types

vcov(fit, type = t) returns the stored covariance when t is the type the fit was made with (e.g. "robust", or "HAC" for ols_robust(), whose stored vcov_type label also shows the kernel and bandwidth). Other types are recomputed from the stored estimates only for fits made by fgls(), serial_fgls() and lasso_fgls(), and for the F-C stage of a combined fit (fit$stages$serial, types "robust", "HC0" and "HC3"). The benchmark estimators store only the covariance they were fitted with; refit them with another vcov instead.

See Also

Other fgls: compare_fgls(), fgls(), fgls_control(), lasso_fgls(), plot.fgls(), serial_fgls(), summary.fgls()

Examples

sim <- simulate_serial(200, seed = 1)
fit <- serial_fgls(y ~ x, data = sim)
vcov(fit)
vcov(fit, type = "naive")
confint(fit)
predict(fit, newdata = data.frame(x = c(-1, 0, 1)), interval = "confidence")
head(residuals(fit, type = "innovations"))
nobs(fit)

Advanced options for fgls

Description

Tuning options of fgls() that most users never change. The defaults follow the papers and the Python package 'pyfgls'.

Usage

fgls_control(
  ic_penalty = c("order", "params"),
  min_lags = NULL,
  fold_scheme = c("random", "blocked"),
  foldid = NULL,
  cv_weights = c("fold", "full"),
  ridge_scale = c("standardized", "original"),
  n_lambda = 100L,
  lambda_min_ratio = NULL,
  lasso_engine = c("lars", "glmnet"),
  glmnet_thresh = 1e-07,
  delta_scale = FALSE,
  transforms = c("linear", "poly2", "log", "trig"),
  svr_cost = c(0.1, 1, 10, 100),
  svr_epsilon = c(0.05, 0.25, 0.5, 1),
  n_boot = 999L,
  bootstrap = c("wild", "iid"),
  wild_dist = c("rademacher", "mammen"),
  boot_seed = NULL,
  cores = 1L
)

Arguments

Value

A list of class "fgls_control".

See Also

Other fgls: compare_fgls(), fgls(), fgls-methods, lasso_fgls(), plot.fgls(), serial_fgls(), summary.fgls()

Examples

fgls_control(ic_penalty = "params", n_boot = 199)

Feasible GLS under heteroskedasticity of unknown form (González-Coya and Perron, 2024)

Description

Shortcut for fgls() with serial = "none". The log of the squared OLS residuals is regressed on candidate covariates z by an adaptive Lasso tuned by cross-validation; the exponentiated fit gives the weights of a weighted least-squares regression. Inference uses the sandwich of eq. (6) with the bias correction of eq. (8), so it remains valid when the skedastic model is misspecified.

Usage

lasso_fgls(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  hetero = c("lasso", "svr", "wls_s1", "wls_s2", "known", "none"),
  vcov = c("robust", "nonrobust", "HC0", "HC1", "HC2", "HC3", "bootstrap"),
  z = "auto",
  seed = NULL,
  ...
)

Arguments

Value

An object of class "fgls"; see fgls().

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

See Also

Other fgls: compare_fgls(), fgls(), fgls-methods, fgls_control(), plot.fgls(), serial_fgls(), summary.fgls()

Other heteroskedasticity: make_dictionary()

Examples

het <- simulate_hetero(200, skedastic = "x2", seed = 1)
fit <- lasso_fgls(y ~ x, data = het, seed = 1)
fit
fit$skedastic$selected
ols <- lasso_fgls(y ~ x, data = het, hetero = "none")
compare_fgls(OLS = ols, FGLS = fit)

Candidate covariates for the skedastic function

Description

Builds the default dictionary of candidate covariates used by fgls() and lasso_fgls() with z = "auto": for every non-constant regressor x, the columns x, x^2, \log x and (\log x)^2 (only when x > 0), \cos x and \cos 2x. This contains the specification z = (1, x, [\log x]^2, x^2, \cos x, \cos 2x) of the simulations of González-Coya and Perron (2024); the intercept of the skedastic regression is always included and never penalized. Binary (dummy) columns are passed through untransformed, undefined or constant transforms are skipped and duplicated columns are removed.

Usage

make_dictionary(
  x,
  transforms = c("linear", "poly2", "log", "trig"),
  skip = NULL
)

Arguments

Value

A numeric matrix with one column per candidate covariate and informative column names.

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

See Also

Other heteroskedasticity: lasso_fgls()

Examples

X <- cbind(const = 1, x = seq(1, 4, length.out = 10), d = rep(0:1, 5))
colnames(make_dictionary(X))
colnames(make_dictionary(X, transforms = c("linear", "poly2")))

Benchmark estimators

Description

Estimators used as comparisons in the papers and in compare_fgls():

Usage

ols_robust(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  vcov = c("HAC", "NW", "HC0", "HC1", "HC2", "HC3", "nonrobust"),
  kernel = NULL,
  bandwidth = NULL,
  small_sample = FALSE,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  model = TRUE
)

cochrane_orcutt(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  ar_order = 1L,
  iterate = FALSE,
  max_iter = 100L,
  tol = 1e-08,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  model = TRUE
)

durbin_estimator(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  ar_order = NULL,
  kmax = NULL,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  model = TRUE
)

infeasible_gls(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  rho = NULL,
  omega = NULL,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  model = TRUE
)

hybrid_gls(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  alpha = 0.05,
  lags = 1L,
  seed = NULL,
  use_t = FALSE,
  time = NULL,
  allow_gaps = FALSE,
  model = TRUE
)

Arguments

Details

• ols_robust(): ordinary least squares with heteroskedasticity and autocorrelation consistent standard errors (vcov = "HAC": Andrews' 1991 quadratic-spectral kernel with automatic bandwidth, the benchmark of PG26; "NW": Newey-West Bartlett kernel with \lfloor 4 (T/100)^{2/9} \rfloor lags) or White-type standard errors ("HC0"–"HC3").

• cochrane_orcutt(): two-step (or iterated) Cochrane-Orcutt FGLS. The two-step version estimates the AR coefficients from OLS residuals and is not valid with non-exogenous regressors (PG26 Section 3); use serial_fgls() instead.

• durbin_estimator(): the Durbin regression (PG26 eq. 14) of y_t on y_{t-j}, x_t and x_{t-j}, j = 1, \dots, k, with the AR order k chosen by BIC unless given. The coefficients on the stochastic regressors are those on x_t. Lags of the constant, of trends and of periodic dummies are collinear with the regressors and dropped, so the Durbin coefficients c on these regressors estimate (I - \sum_j \rho_j P_j)\beta, where the dropped lags are x_{t-j} = x_t P_j (plus kept lags). They are converted back to \hat\beta = (I - \sum_j \hat\rho_j P_j)^{-1} c, e.g. c / (1 - \sum_j \hat\rho_j) for the intercept, with a delta-method covariance, and are NA (with a note) when that is not possible. Residuals, fitted values and predictions use these coefficients; innovations are the residuals of the Durbin regression. (The Python package 'pyfgls' reports the untransformed c.)

• infeasible_gls(): GLS with known AR coefficients (rho) or a known error covariance matrix (omega, up to scale); a simulation benchmark.

• hybrid_gls(): the hybrid OLS/GLS estimator of PG26 Supplement S-1 (exogenous regressors only): OLS when the Breusch-Godfrey test does not reject, otherwise quasi-differencing with a random AR(1) coefficient of the sign of the first-order residual autocorrelation.

Value

An object of class "fgls" (see fgls()). The non-exogeneity share of PG26 eq. (29) is not computed for these estimators and is NA.

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

Examples

sim <- simulate_serial(300, seed = 1)
ols_robust(y ~ x, data = sim)
cochrane_orcutt(y ~ x, data = sim)
durbin_estimator(y ~ x, data = sim)
infeasible_gls(y ~ x, data = sim, rho = 0.8)

Diagnostic plots for fgls fits

Description

Base-graphics diagnostics: (1) residuals against fitted values; (2) the autocorrelation function of the innovations (the residuals of the final filtered and/or weighted regression, which should look like white noise); (3) the lag-order criterion path with the selected order (serial correction); (4) the fitted variances against the fitted values (heteroskedasticity correction), or a normal quantile plot of the innovations otherwise.

Usage

## S3 method for class 'fgls'
plot(x, which = 1:4, ask = FALSE, lag.max = 20L, ...)

Arguments

Value

x, invisibly.

See Also

Other fgls: compare_fgls(), fgls(), fgls-methods, fgls_control(), lasso_fgls(), serial_fgls(), summary.fgls()

Examples

fit <- fgls(y ~ x, data = simulate_combined(200, seed = 1), seed = 1)
plot(fit)

Feasible GLS for serially correlated errors (Perron and González-Coya, 2026)

Description

Shortcut for fgls() with hetero = "none": the errors are approximated by an AR(k) process whose coefficients come from the Durbin regression and whose order is chosen by BIC (eq. 14), the data are quasi-differenced and OLS is applied (eq. 15). The default variance is the exogeneity-robust estimator of eqs. (28)–(29), valid when the regressors are only pre-determined.

Usage

serial_fgls(
  formula,
  data,
  subset,
  na.action = stats::na.fail,
  serial = c("auto", "fixed", "ma1"),
  vcov = c("robust", "naive", "HC0", "HC1", "HC2", "HC3", "bootstrap"),
  ar_order = NULL,
  kmax = NULL,
  ic = c("bic", "aic"),
  lagged_dep_vars = NULL,
  exog = c("predetermined", "exogenous"),
  ...
)

Arguments

Value

An object of class "fgls"; see fgls().

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

See Also

Other fgls: compare_fgls(), fgls(), fgls-methods, fgls_control(), lasso_fgls(), plot.fgls(), summary.fgls()

Examples

sim <- simulate_serial(500, seed = 42)
fit <- serial_fgls(y ~ x, data = sim)
fit
fit$ar_order
fit$nonexogeneity_share

# Lagged dependent variable among the regressors (Wallis recovery)
df <- data.frame(y = sim$y[-1], ylag = sim$y[-500], x = sim$x[-1])
serial_fgls(y ~ ylag + x, data = df, lagged_dep_vars = "ylag")

Simulate the designs of the papers

Description

Data-generating processes used in the Monte Carlo studies of Perron and González-Coya (2026) (simulate_serial(), simulate_combined(), simulate_predictive()) and González-Coya and Perron (2024) (simulate_hetero()). They are handy for examples, teaching and replication.

Usage

simulate_serial(
  n,
  ar = 0.8,
  ma = numeric(0),
  rho_x = 0,
  gamma = 1,
  beta = c(1, 1),
  burn = 0,
  seed = NULL
)

simulate_hetero(
  n,
  skedastic = "x2",
  beta = c(1, 1),
  n_extra = 0L,
  extra_corr = 0,
  seed = NULL
)

simulate_combined(
  n,
  ar = 0.5,
  ma = numeric(0),
  skedastic = "x2",
  beta = c(1, 1),
  driver = c("x", "mixture"),
  phi_prob = 0.75,
  gamma = 0,
  burn = 100L,
  seed = NULL
)

simulate_predictive(
  n,
  theta = 0.5,
  rho_x = 0,
  gamma = 0,
  beta = c(1, 1),
  burn = 100L,
  seed = NULL
)

Arguments

Details

• simulate_serial(): y_t = \beta_0 + \beta_1 x_t + u_t with ARMA errors driven by standard normal innovations e_t and a regressor x_t = \rho_x x_{t-1} + v_t + \gamma e_{t-1}, which is pre-determined and exogenous only when \gamma = 0. The defaults (\rho_x = 0, \gamma = 1, AR(1) coefficient 0.8, no burn-in, u_0 = 0) give the design of Table 1 of the 2026 paper.

• simulate_hetero(): y_t = \alpha + \beta x_t + \sqrt{v(x_t)} \varepsilon_t, x_t \sim U[1, 4] (eq. 9 of the 2024 paper), optionally with extra candidate covariates w_j = c x + (1 - c) q_j, q_j \sim U[1, 4].

• simulate_combined(): ARMA errors whose innovations are heteroskedastic, e_t = \sqrt{v(w_t)} \varepsilon_t, with w_t = x_t (driver = "x") or w_t = \phi_t x_t + (1 - \phi_t) q_t, \phi_t \sim Bernoulli(phi_prob) (driver = "mixture") (Supplement S-4).

• simulate_predictive(): the predictive regression y_{t+2} = \alpha + \beta x_t + u_{t+2} with MA(1) errors (Supplement S-3), returned aligned so that row t holds y_{t+2} and x_t.

Value

A data frame with columns y and x (and w or ⁠w1, w2, ...⁠ where relevant), with attributes beta (true coefficients), u (errors), e (innovations) and info (design parameters and, for the heteroskedastic designs, the true variances v).

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Examples

sim <- simulate_serial(200, seed = 1)
head(sim)
attr(sim, "beta")

het <- simulate_hetero(100, skedastic = "log2", seed = 2)
summary(attr(het, "info")$v)

Summarize an fgls fit

Description

The summary contains the coefficient table (estimates, standard errors, z or t statistics, p-values and confidence intervals) and the main by-products of the estimator: the selected AR order k^* and AR coefficients with the non-exogeneity share for serial correlation, and the selected candidate covariates, the tuning parameters and the degrees of freedom of the skedastic fit for heteroskedasticity.

Usage

## S3 method for class 'fgls'
summary(object, level = 0.95, ...)

## S3 method for class 'summary.fgls'
print(
  x,
  digits = max(3L, getOption("digits") - 3L),
  signif.stars = getOption("show.signif.stars"),
  ...
)

Arguments

Value

summary() returns an object of class "summary.fgls": a list with the call, the coefficient table (coefficients), the confidence intervals (conf.int), a named list of header items (info) and the notes. Its print() method returns its argument invisibly.

See Also

Other fgls: compare_fgls(), fgls(), fgls-methods, fgls_control(), lasso_fgls(), plot.fgls(), serial_fgls()

Examples

sim <- simulate_serial(200, seed = 1)
s <- summary(serial_fgls(y ~ x, data = sim))
s
s$coefficients

Test whether lags of the regressors are missing (pre-determinedness)

Description

Experimental. A variable-addition LM test in the spirit of Remark 9 of Perron and González-Coya (2026): FGLS needs the regressors to be pre-determined (uncorrelated with current and future innovations), which can fail when lags of the regressors belong in the model. The steps are (a) the Durbin regression (AR order by BIC) giving \hat\beta^D; (b) the residuals \hat u^D = y - X\hat\beta^D, demeaned; (c) for l = 1, \dots, max_lag, the statistic LM_l = n (SSR_0 - SSR_1)/SSR_0, asymptotically \chi^2 with as many degrees of freedom as added lags, from the auxiliary regressions without and with lags 1..l of the stochastic regressors.

Usage

test_predetermined(
  formula,
  data,
  max_lag = 4L,
  ar_order = NULL,
  kmax = NULL,
  regressors = NULL,
  alpha = 0.05,
  scale = c("innovations", "errors"),
  na.action = stats::na.fail,
  time = NULL,
  allow_gaps = FALSE
)

## S3 method for class 'fgls_predetermined'
print(x, ...)

Arguments

Details

With scale = "innovations" (default) the auxiliary regression is run on the innovation scale: \hat u^D and all regressors are filtered with the Durbin AR filter and lagged \hat u^D enter as nuisance regressors, which keeps the \chi^2 approximation valid with pre-determined but non-exogenous regressors. scale = "errors" uses \hat u^D and the raw regressors (the literal reading of Remark 9).

Value

An object of class "fgls_predetermined": a list with

• table: a data frame with columns lags, LM, df and p.value, one row per number of added lags;

• ar_order: the AR order of the Durbin regression;

• regressors: the names of the regressors whose lags were tested;

• alpha: the level used for the advice;

• rejected: TRUE if some test rejects at level alpha;

• advice: a sentence summarizing the result;

• scale: the scale used, "innovations" or "errors".

The print method returns x invisibly.

References

Perron, P. and González-Coya, E. (2026). Feasible GLS for time series regression. Journal of Econometrics, 258, 106346. doi: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. doi:10.1515/jem-2023-0007

Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003

See Also

Other diagnostics: adf_test(), weight_diagnostics()

Examples

sim <- simulate_serial(300, seed = 1)
test_predetermined(y ~ x, data = sim)

Diagnostics of an fgls fit

Description

• weight_diagnostics(): how concentrated the FGLS weights 1/\hat v_t are: the ratio of the largest to the smallest fitted variance, the Kish effective sample size (\sum w)^2 / \sum w^2 as a share of the sample, and the share of the total weight carried by the top 1% of observations.

• ar_roots(): roots of the estimated lag polynomial 1 - \hat\rho_1 z - \dots - \hat\rho_k z^k; the filter is invertible when all lie outside the unit circle.

• nonexogeneity_share(): the relative size of the first-stage correction of PG26 eq. (29), averaged over coefficients. Values near zero mean the regressors behave as exogenous; large values mean the naive variance (eq. 27) would understate the uncertainty. It is computed from the robust variance whatever vcov the fit uses, and is NA for the benchmark estimators (see ols_robust()) and for serial = "ma1", which do not compute it.

Usage

weight_diagnostics(fit)

ar_roots(fit)

nonexogeneity_share(fit)

Arguments

Value

• weight_diagnostics(): a named numeric vector with elements variance_ratio, effective_n_share and top1pct_weight_share.

• ar_roots(): a complex vector with the k^* roots (empty when k^* = 0).

• nonexogeneity_share(): a single number.

See Also

Other diagnostics: adf_test(), test_predetermined()

Examples

fit <- fgls(y ~ x, data = simulate_combined(200, seed = 1), seed = 1)
weight_diagnostics(fit)
Mod(ar_roots(fit))
nonexogeneity_share(fit)