---
title: "Getting started with fgls"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with fgls}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

```{r setup}
library(fgls)
```

**fgls** estimates linear regressions by *feasible generalized least squares*
(FGLS). It implements two recent methods:

* **Serial correlation** (Perron and González-Coya, 2026): the errors are
  approximated by an autoregression whose order is chosen by BIC and whose
  coefficients come from a *Durbin regression*. The resulting estimator is
  consistent even when the regressors are only *pre-determined*, a situation
  where OLS and Cochrane-Orcutt are not. A variance estimator accounts for the
  estimated filter.
* **Heteroskedasticity of unknown form** (González-Coya and Perron, 2024): the
  log of the squared residuals is modeled by an *adaptive Lasso* over a
  dictionary of candidate covariates, and the weighted regression is paired with
  a bias-corrected sandwich variance that stays valid when that model is wrong.

The two corrections can be combined (Supplement S-4 of the 2026 paper). The R
package is a port of the Python package
[pyfgls](https://pypi.org/project/pyfgls/), and the two agree to machine
precision on deterministic quantities (see `vignette("validation")`).

## Three examples

### Serial correlation

`simulate_serial()` generates the design of Table 1 of the 2026 paper:
AR(1) errors with coefficient 0.8, and a regressor that responds to past
innovations, \(x_t = v_t + e_{t-1}\), so it is pre-determined but not exogenous.

```{r serial}
sim <- simulate_serial(500, seed = 1)
fit <- serial_fgls(y ~ x, data = sim)
summary(fit)
```

The true slope is 1. OLS is badly biased in this design; FGLS is not:

```{r serial-compare}
compare_fgls(OLS = lm(y ~ x, data = sim), FGLS = fit)
```

`fit$ar_order` is the order chosen by BIC, `fit$rho` the AR coefficients and
`fit$nonexogeneity_share` the relative size of the correction that makes the
variance robust to non-exogenous regressors.

### Heteroskedasticity of unknown form

```{r hetero}
het <- simulate_hetero(200, skedastic = "log2", seed = 2)
fit_h <- lasso_fgls(y ~ x, data = het, seed = 1)
fit_h$skedastic$selected
compare_fgls(OLS = lm(y ~ x, data = het), FGLS = fit_h)
```

The candidate covariates default to \(x, x^2, \log x, (\log x)^2, \cos x,
\cos 2x\) for every non-constant regressor (`make_dictionary()`); supply your own
with `z`. The penalty and the exponent of the adaptive weights are chosen by
10-fold cross-validation, so set `seed` for reproducible results.

### Both

```{r combined}
both <- simulate_combined(200, ar = 0.5, skedastic = "x2", seed = 3)
fit_b <- fgls(y ~ x, data = both, seed = 1)
fit_b
compare_fgls(`F-C` = fit_b$stages$serial, `F-CH` = fit_b)
```

`fgls()` is the general function: `serial = "none"` gives the
heteroskedasticity-only estimator and `hetero = "none"` the serial-only one.

## Working with fits

Fits behave like `lm` objects:

```{r methods}
coef(fit)
confint(fit, level = 0.9)
predict(fit, newdata = data.frame(x = c(-1, 0, 1)), interval = "confidence")
vcov(fit, type = "naive") # PG26 eq. (27), valid only for exogenous regressors
```

They work with `lmtest::coeftest()`, `lmtest::coefci()` and the broom generics
`tidy()`, `glance()` and `augment()`. `plot(fit)` draws diagnostics: the
autocorrelation function of the innovations (residuals of the final filtered
regression) should look like white noise.

## Time series inputs

Rows must be in time order. Besides data frames, `fgls()` accepts `ts`, `zoo`
and `xts` objects; a data frame can carry a time index through `time =`, which is
checked for sorting and gaps:

```{r ts}
dates <- seq(as.Date("2000-01-01"), by = "month", length.out = 500)
df <- data.frame(date = dates, sim)
serial_fgls(y ~ x, data = df, time = "date")$ar_order
```

## Choosing options

* With several stochastic regressors, the BIC that penalizes only the AR order
  can select the largest order; use `control = fgls_control(ic_penalty =
  "params")`.
* Declare lagged dependent variables among the regressors in `lagged_dep_vars`.
* For errors known to be MA(1), such as multi-step forecast errors, use
  `serial = "ma1"`.
* `exog = "exogenous"` never changes the estimates. In the serial-only
  estimator it only changes the notes, and the robust variance is valid either
  way. In the combined estimator it also drops the first-stage correction from
  the variance by default (`first_stage_correction = FALSE`); that correction is
  needed when the regressors are only pre-determined, so keep the default
  `exog = "predetermined"` unless they are known to be exogenous.

See `vignette("serial-correlation")`, `vignette("heteroskedasticity")` and
`vignette("combined")` for the methods in more detail.

## References

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>