Bias matrix relating raw residual autocovariance to the truth

Builds the linear map A satisfying E[gamma_raw] = A gamma_true induced by projecting a design out of the data. Autocovariance estimated from GLM residuals is biased low, increasingly so for richer designs, because E[ehat ehat'] = M Sigma M for the residual-forming projection M. Solving A gamma = gamma_raw removes that bias.

Usage

acvf_bias_matrix(design, runs = NULL, censor = NULL, max_lag = 25L)

Arguments

design: Numeric design matrix (timepoints x regressors), the one whose projection produced the residuals.
runs: Optional run labels, length nrow(design). Each label must occupy one contiguous block and may not be missing. Lag products never cross a run boundary.
censor: Optional indices of censored timepoints (1-based), excluded exactly as in fit_noise().
max_lag: Highest lag to correct.

Value

A named list of (max_lag + 1) square matrices, one per run, rows indexed by the lag of the raw estimate and columns by the lag of the truth.

Details

Ordinarily you do not need this: pass design to fit_noise() and it builds and applies the matrices itself. Use this directly to inspect the bias, or to cache the matrices when many datasets share one design and segmentation and feed them back as acvf_correction.

One matrix is returned per run, because fit_noise() estimates each run separately and the residual operator restricted to a run still involves the whole design. Under pooling = "parcel", which pools lag products across runs, the run matrices are combined in proportion to the lag pairs each run contributes.

Projecting out run means and slow drift leaves the map with near-null directions (roughly a constant offset across lags). When the correction is applied, those directions are fixed by the same short-memory assumption the truncated system already makes (the autocovariance is near zero at the end of the budget) rather than taken from the data, which keeps the correction stable under censoring.

The correction is exact for noise whose autocovariance dies within max_lag. Under long memory it is partial, since truncating the system aliases in structure past the budget. With high-pass filtering (DCT, 128 s) a budget near 20-25 lags is typically enough; without high-pass it exceeds 150 and the approach stops being practical.

Cost is O(max_lag^2 * n_valid * r) per run for a design of rank r.

Examples

n <- 120
X <- cbind(1, poly(seq_len(n), 3))
A <- acvf_bias_matrix(X, max_lag = 5)
round(A[[1]][1:3, 1:3], 3)
#>        [,1]   [,2]   [,3]
#> [1,]  0.967 -0.064 -0.062
#> [2,] -0.032  0.935 -0.063
#> [3,] -0.032 -0.063  0.937