ENH: Add the Andrews (1991) AR(1) plug-in bandwidth to the kernel estimators - #964

Target 9.1: merge after v9.0.0 is tagged. Closes #887.

Adds bandwidth_method="newey-west" | "andrews" to every kernel covariance estimator
(including the Andrews, Gallant, NeweyWest and ZeroLag aliases), to
PreWhitenedRecolored, which passes it to its kernel, and to the fit methods of
DynamicOLS, FullyModifiedOLS and CanonicalCointegratingReg. It also makes the
automatic estimates of every kernel the estimators of Andrews (1991) and Newey and West
(1994).

bandwidth_method:

• It accepts the spellings that kernel names accept: case, hyphens and underscores are
  ignored, so "Newey_West" and "neweywest" are "newey-west". It is exposed as a
  read-only bandwidth_method property.
• "newey-west" is the default.
• "andrews" implements Andrews (1991, Section 6, eqs. 6.3-6.4).
  • An AR(1) is fitted to each column: OLS with an intercept, and sigma^2 = RSS/(T-1),
    as R's ar(method="ols") does.
  • The optimal bandwidth is c_k (alpha(q) T)^(1/(2q+1)). The kernel constants are
    Andrews' constants (Bartlett 1.1447, Parzen 2.6614, QS 1.3221, Tukey-Hanning
    1.7462), and the derivation reproduces them.
  • Weights are respected. Constant and zero-weight columns are skipped.
  • Complex data are handled as separate real and imaginary columns.
  • An explicit bandwidth overrides the method.

Kernel weights (changes results). The kernels that are 0 beyond 1 (all except QS)
weight lag h by k(h/(H+1)) for a bandwidth H, so that H is the number of lags, as in
Newey and West (1987): Bartlett(x, bandwidth=L) is the estimator with L lags, as
sandwich::NeweyWest(lag=L), Stata's newey and statsmodels. The optimal bandwidth S of
Andrews (1991) and Newey and West (1994) is that of k(h/S), and it was used as H, so
the automatic estimates used the weights of S + 1.

• The automatic bandwidth is now H = max(S - 1, 0), whose weights are those of S. QS,
  which weights every lag by k(h/H), keeps H = S. opt_bandwidth is still S.
• With force_int=True, H is rounded up, which is the integer part of S, the
  m = [gamma T^(1/3)] of Newey and West (1994, eq. 2-2).
• A fractional bandwidth now uses every lag with a positive weight. The lag
  floor(H) + 1 was dropped.
• An explicit integer bandwidth is unchanged, so the EViews-validated cointegration
  results, which pass integer bandwidths, are unchanged.
• The docstrings described k(h/H) and now describe the weights that are used.
• The automatic long-run covariances of the bounded kernels fall, by 0.6% to 4.9% in
  AR(1) examples.

Pilot lag. This branch includes the commit of #983 (9.0), which uses the integer
part of 4 (T/100)^r for the pilot lag of "newey-west", as Newey and West (1994),
auto_bandwidth and sandwich::bwNeweyWest do. It drops out when this branch is rebased
on main after #983 is merged.

Other changes:

• Extreme scales. Both estimates of alpha(q) are ratios that do not depend on the
  scale of the data. The data are divided by the power of 2 nearest their largest
  absolute value first. This is exact, so results are unchanged, and data scaled by
  1e-200 to 1e200 no longer overflow or underflow. All columns share one scale, since a
  separate scale for each column would change their weights in Andrews' estimate.
• Tukey-Parzen constant. It was 1.8576. The formula gives 1.85874, so it is now
  1.8587, and its optimal bandwidth is 0.059% larger. Every constant is now the formula
  rounded to 4 decimals.

Validation: R 4.6.1 with sandwich 3.1.3, 6 seeded datasets (univariate and
multivariate, T = 50 to 500, nonzero means, non-uniform and zero weights). The tests
store the R values as constants and use rtol 1e-10; the largest difference is 2.2e-15.

• "andrews": opt_bandwidth equals bwAndrews(x, kernel=..., approx="AR(1)", weights=w, prewhite=FALSE), and the long-run covariance equals
  kernHAC(lm(x ~ 1), bw=bwAndrews, kernel=..., prewhite=FALSE, adjust=FALSE), for the
  Bartlett, Parzen, QS and Tukey-Hanning kernels.
• "newey-west": opt_bandwidth equals bwNeweyWest, the long-run covariance equals
  kernHAC with that bandwidth (Bartlett, Parzen, QS), and Bartlett with
  force_int=True equals NeweyWest(lm(x ~ 1), lag=floor(bwNeweyWest)).
• A fractional bandwidth equals NeweyWest(lag=5.4) and kernHAC(bw=6.4).
• For every kernel, the automatic estimate equals an independent sum of k(h/S) times
  the autocovariances.
• The kernel constants are checked against the formula by quadrature.

Tests:

• arch/tests/covariance/test_andrews_bandwidth.py (242 tests), including the R
  comparisons above, the accepted spellings, and the scale invariance from 1e-200 to
  1e200 (exact for powers of 2, without warnings).
• DOLS, FMOLS and CCR: the default equals "newey-west", "andrews" gives the bandwidth
  of the kernel estimator on the same residuals or scores, an explicit bandwidth
  overrides the method, and the spellings are accepted (60 tests). Validation rows for
  every estimator.
• Fast suite 27,547 passed, and 27,601 after a rebase on main with BUG: Use the integer part for the pilot lag of the automatic bandwidth #983. flake8, ruff,
  black, isort, numpydoc lint and mypy are clean on the changed files.

For review:

• PhillipsPerron and phillips_ouliaris could also take bandwidth_method. That belongs
  in a separate issue.

Thanks to @ZacKienzle2 for the proposal.

🤖 Generated with Claude Code