Serial correlation and heteroskedasticity together

library(fgls)

Section S-4 of the Supplement to Perron and González-Coya (2026) combines the two corrections:

  • F-C corrects for serial correlation: Durbin regression, BIC lag selection, quasi-differencing;
  • F-CH then models the variance of the filtered residuals with the adaptive Lasso of González-Coya and Perron (2024). The candidate covariates are built from the unfiltered regressors (or supplied in z), and the quasi-differenced regression is re-estimated by weighted least squares.

The variance of F-CH is the bias-corrected sandwich of the 2024 paper applied to the filtered regressors (Supplement eqs. A.4–A.5). F-C intervals are built the same way with unit weights.

sim <- simulate_combined(200, ar = c(1.34, -0.42), skedastic = "exp", seed = 5)
fit <- fgls(y ~ x, data = sim, seed = 1)
summary(fit)
#> 
#> Call:
#> fgls(formula = y ~ x, data = sim, seed = 1)
#> 
#> Method:                 F-CH (serial+hetero)
#> Covariance:             robust
#> Observations:           200 (effective 198)
#> Inference:              normal
#> AR order k*:            2 [BIC, kmax = 12]
#> AR coefficients:        1.3981, -0.4621
#> Non-exogeneity share:   0.011
#> Skedastic fit:          lasso
#> Selected z:             log(x), log(x)^2
#> Lasso psi, lambda:      2, 0.06188 (glmnet scale; 12.25 in GP24 eq. 4)
#> df_hat:                 3
#> Variance ratio max/min: 42.7
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  -0.5773     2.5323  -0.228     0.82    
#> x             1.0163     0.1035   9.821   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Notes:
#> [1] Covariance: PG26 Supplement (A.4)-(A.5) sandwich (robust), with the
#>     first-stage (Durbin) correction.
compare_fgls(OLS = lm(y ~ x, data = sim), `F-C` = fit$stages$serial, `F-CH` = fit)
#>          term statistic        OLS        F-C       F-CH
#> 1 (Intercept)  estimate  3.0197767  3.4310357 -0.5772954
#> 2 (Intercept) std.error  3.1225187  3.7986467  2.5323070
#> 3 (Intercept)  ci.width 12.2400485 14.8904213  9.9264609
#> 4           x  estimate  0.7746907  0.9815785  1.0162912
#> 5           x std.error  0.8656833  0.2010900  0.1034812
#> 6           x  ci.width  3.3934163  0.7882582  0.4056387

fgls() is the general function; this combined estimator is its default.

The first-stage correction

The Supplement does not write out the correction for the estimated AR filter that eq. (29) of the paper provides in the serial-only case. fgls adds the heteroskedasticity-robust analogue of that term by default when the regressors are declared pre-determined (first_stage_correction = TRUE, the default for exog = "predetermined"). With exogenous regressors it is negligible. With non-exogenous ones it is needed for correct coverage. The package’s slow Monte Carlo acceptance tests (tests/testthat/test-statistical-s4.R in the source package, run with FGLS_SLOW_TESTS=true) check this with AR(1) errors with coefficient 0.8, \(v(x) = x^2\), a regressor that responds to past innovations (\(\gamma = 0.5\) and 1) and \(T = 200\): the corrected 95% interval must cover between 91% and 97% of the time, and at least three percentage points more often than the uncorrected one. The replication notes of pyfgls report coverage of 0.923 with the correction and 0.860 without it for \(\gamma = 1\) (1,000 replications).

sim_ne <- simulate_combined(300, ar = 0.8, skedastic = "x2", gamma = 1, seed = 6)
with_fs <- fgls(y ~ x, data = sim_ne, seed = 1)
without_fs <- fgls(y ~ x, data = sim_ne, seed = 1, first_stage_correction = FALSE)
c(with = sqrt(vcov(with_fs)[2, 2]), without = sqrt(vcov(without_fs)[2, 2]))
#>       with    without 
#> 0.04187342 0.03825273

No loss under homoskedasticity

When the innovations are homoskedastic the Lasso often selects nothing, and F-CH then equals F-C, as in this example:

hom <- simulate_combined(200, ar = 0.5, skedastic = "homoskedastic", seed = 4)
fit_h <- fgls(y ~ x, data = hom, seed = 1)
fit_h$skedastic$selected
#> character(0)
all.equal(coef(fit_h), coef(fit_h$stages$serial))
#> [1] TRUE

In other samples it selects a few covariates by chance, and F-CH then differs somewhat from F-C. On average nothing is lost: in the homoskedastic row of Table S.9 (vignette("validation")) the two have essentially the same mean squared error and both cover close to 95%.

Bootstrap cross-check

vcov = "bootstrap" re-runs the whole estimator on residual-based sieve bootstrap samples, using the estimated AR filter with wild innovations. It is an independent check of the analytic variance. The regressors are held fixed, so with non-exogenous regressors it does not reproduce their feedback.

fb <- fgls(y ~ x, data = sim, seed = 1, vcov = "bootstrap", control = fgls_control(n_boot = 19))
cbind(analytic = sqrt(diag(vcov(fit))), bootstrap = sqrt(diag(vcov(fb))))
#>              analytic bootstrap
#> (Intercept) 2.5323070 2.7106793
#> x           0.1034812 0.1042649