---
title: "Validation: replicating the papers and agreement with pyfgls"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Validation: replicating the papers and agreement with pyfgls}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

The results below were **precomputed** (the Monte Carlo simulations take tens of
minutes) by `tools/replicate_tables.R` and `tools/compare_python_r.R` in the
source repository, and are stored in the package:

```{r load}
rep <- readRDS(system.file("extdata", "replication.rds", package = "fgls"))
xl <- readRDS(system.file("extdata", "cross_language.rds", package = "fgls"))
rep$replications
```

The papers use 10,000 replications; with 1,000 the Monte Carlo standard error
of a coverage rate near 0.95 is about 0.007. The same designs are run, with
tolerances, by the package's acceptance tests (`FGLS_SLOW_TESTS=true`).

## Perron and González-Coya (2026), Table 1

AR(1) errors (0.8), a regressor that responds to past innovations, \(T = 500\);
slope estimates (true value 1):

```{r t1}
knitr::kable(rep$pg26_table1, digits = 3)
rep$pg26_table1_median_k
```

OLS is inconsistent, Cochrane-Orcutt is biased (its estimate of \(\rho\) is
biased downwards), and FGLS matches the infeasible GLS.

## Perron and González-Coya (2026), Tables 2--5 (selected cells)

\(T = 200\), \(x_t = 0.8 x_{t-1} + v_t + \gamma e_{t-1}\); MSE \(\times 100\) and
coverage of 95% intervals with the robust (eq. 29) and naive (eq. 27)
variances:

```{r t25}
knitr::kable(rep$pg26_general, digits = 3)
```

The robust intervals have close to nominal coverage for all \(\gamma\); the naive
ones under-cover when the regressor is not exogenous (\(\gamma > 0\)).

## González-Coya and Perron (2024), Table 1

Empirical MSE of the adaptive-Lasso FGLS slope relative to OLS, \(T = 400\), with
the paper's values, and coverage of the eq. (6)+(8) intervals:

```{r gp24}
knitr::kable(rep$gp24_table1, digits = 3)
```

At \(T = 100\)--\(200\) the ratios of this implementation (and of pyfgls) are a few
hundredths above the paper's; the T = 400 column is matched.

## Supplement, Table S.9: F-C and F-CH

\(T = 200\); MSE \(\times 100\), coverage and average standard errors:

```{r s9}
knitr::kable(rep$s9, digits = 3)
```

F-CH improves on F-C under heteroskedasticity and loses nothing under
homoskedasticity. With 1,000 replications some F-CH coverage rates are
slightly above 0.95; the intervals are conservative, not too short.

## Agreement with pyfgls

`fgls` is a port of the Python package pyfgls. Both packages fitted the same
`r length(xl$designs[[1]]$seeds)` data sets per design with the same
cross-validation folds (pyfgls `r xl$designs[[1]]$pyfgls`):

```{r xl}
s <- xl$summary
s <- s[s$quantity %in% c("beta", "se", "kstar", "selected", "coverage_95"), ]
knitr::kable(s, digits = 15, row.names = FALSE)
```

Serial-correlation quantities agree to rounding error, and the selected AR
order is identical in every data set. In about 1% of the Lasso data sets the
degrees of freedom differ by one, because pyfgls 0.1.0 counts coefficients of
order \(10^{-17}\) (left by interpolating its Lasso path at a knot) as selected.
With the fixes merged into pyfgls after 0.1.0 (`emiliogcs/pyfgls#2`), everything
agrees and the selected sets are identical:

```{r xl-dev}
xl$pyfgls_pr2$version
knitr::kable(xl$pyfgls_pr2$summary, digits = 15, row.names = FALSE)
```

The deterministic parts are also checked against golden fixtures generated by
pyfgls (`tests/testthat/test-golden-*.R`), to relative precision \(10^{-8}\) or
better.