Comparing inference methods on one estimate

The package computes inference for a fixed-dimensional time-series parameter by six methods. All of them start from the same estimate and the same observation-level influence contributions, so a difference between their answers comes from the normalizer and its reference law, not from a difference in the setup. This vignette shows how to run them, how to read the comparison table, and how to change the tuning values that four of them require. The data are synthetic throughout.

The default method is the affine-equivariant adjusted-range increment hull, which is the subject of the package’s main reference. The other five are provided so that it can be compared with methods already in use.

The six methods

`method`	Normalizer	Tuning	Reference law
`"hull"`	increment hull of the centered influence path	none	simulated Brownian gauge law
`"ldl"`	componentwise adjusted ranges after lag-zero prewhitening	coordinate order	simulated independent-component law
`"shao"`	integrated outer product of the path	integration rule	simulated Brownian quadratic law
`"hac"`	kernel long-run covariance estimate	kernel, bandwidth	chi-squared
`"fixedb"`	Bartlett estimate with bandwidth fraction b	b	simulated fixed-b law
`"ewc"`	equal-weighted cosine (EWC) estimate	number of terms	scaled F

The first three are self-normalized: the normalizer has a nondegenerate random limit, the unknown long-run covariance cancels, and no smoothing parameter is chosen. "hac" estimates the long-run covariance matrix consistently and uses the usual chi-squared reference. "fixedb" and "ewc" use reference laws derived by holding their smoothing parameter fixed as the sample grows, retaining the randomness of the covariance estimate. For EWC, the default number of terms increases with sample size; the scaled F reference uses the selected number at each sample size. HAC stands for heteroskedasticity and autocorrelation consistent; HAR, as in the EWC literature, stands for heteroskedasticity and autocorrelation robust.

One estimate, six answers

library(aersn)
set.seed(2026)
n <- 400
e <- matrix(rnorm(2 * n), n, 2)
Y <- e
for (t in 2:n) Y[t, ] <- 0.5 * Y[t - 1, ] + e[t, ]
Y <- sweep(Y, 2, c(0.12, -0.05), `+`)

fit <- aersn_mean(Y, names = c("m1", "m2"))
fit
#> Affine-equivariant adjusted-range self-normalization (increment hull)
#>   model: sample mean (psi_t = Y_t - Ybar) 
#>   n = 400 observations; q = 2 parameters
#>   centering: calendar time, tau(r) = r 
#>   estimate:
#>      m1      m2 
#> 0.14530 0.03354 
#>   hull diagnostics: numerical rank 2 of 2  ; condition 1.953 ; min projected range 2.158 ; spread 1.510 
#> Use aersn_test(), confint(), aersn_contrast(), aersn_region(), plot().

aersn_compare() runs the methods on the same null value and level:

cmp <- aersn_compare(fit, null = c(0, 0), draws = 2000, seed = 1)
cmp
#> Comparison of inference methods on one estimate
#>   n = 400   q = 2   level = 0.95
#>   null value: m1 = 0, m2 = 0
#>   estimate:   m1 = 0.1453, m2 = 0.03354
#>  method statistic scale critical p     reject half_width
#>  hull   1.374     gauge 2.398    0.335 no     0.3177    
#>  ldl    1.286     wald  4.915    0.412 no     0.2936    
#>  shao   24.24     wald  96.15    0.348 no     0.3213    
#>  hac    3.315     wald  5.991    0.191 no     0.2133    
#>  fixedb 4.679     wald  25.32    0.435 no     0.3901    
#>  ewc    2.748     wald  7.335    0.292 no     0.2596    
#> 
#> Statistics are on different scales and are not comparable as numbers: 
#>   the hull statistic is a gauge, the others are Wald statistics, 
#>   and each is referred to its own law. The comparable columns are 
#>   p, reject and half_width (for the first coordinate).
#> 
#> Tuning actually used:
#>   hull   none
#>   ldl    order=1, 2; lag_zero_divisor=n; factor=unit lower triangular L
#>   shao   integration=calendar
#>   hac    kernel=Bartlett; bandwidth_rule=short (floor(4 (n/100)^(2/9))); bandwidth=6; lag_truncation=5; prewhite=FALSE; adjust=FALSE; center=TRUE
#>   fixedb b=0.5; b_grid=0.5; m=200; kernel=Bartlett
#>   ewc    nu=21; rule=floor(0.4 n^(2/3)); center=TRUE; basis=type II cosine

Read the table by column. The statistic and critical_value columns are not comparable across rows: the hull statistic is a gauge, homogeneous of degree one in the estimation error, and the other five are Wald statistics, homogeneous of degree two. Each is referred to its own law, so their magnitudes have different units. The columns that can be compared are p_value, reject and half_width, the last being the half-width of the interval for the first coordinate. The tuning column records the values actually used, including any that were selected from the data or rounded.

The same information is available one method at a time:

aersn_test(fit, null = c(0, 0), method = "shao", draws = 2000, seed = 1)
#> 
#>  Quadratic self-normalization (Shao) test, q = 2
#> 
#> data:  fit
#> T = 24.24, q = 2, n = 400
#> alternative hypothesis: true parameter is not equal to the null value
#> null value: m1 = 0, m2 = 0 
#> estimate:   m1 = 0.1453, m2 = 0.03354 
#> critical value at level 0.95: 96.15  (do not reject)
#> tuning actually used: integration=calendar
#> p-value = 0.348 (Monte Carlo s.e. 0.011, resolution 0.00050)
#> reference: matched-grid Monte Carlo law for quadratic self-normalization: q = 2, uniform grid with n = 400 intervals, 2000 draws, seed 1, calendar-time integration
aersn_test(fit, null = c(0, 0), method = "hac", kernel = "Parzen",
           bandwidth = "andrews")
#> 
#>  HAC long-run covariance test, q = 2
#> 
#> data:  fit
#> T = 2.814, q = 2, n = 400
#> alternative hypothesis: true parameter is not equal to the null value
#> null value: m1 = 0, m2 = 0 
#> estimate:   m1 = 0.1453, m2 = 0.03354 
#> critical value at level 0.95: 5.991  (do not reject)
#> tuning actually used: kernel=Parzen; bandwidth_rule=Andrews (1991) plug-in; bandwidth=15.6197; lag_truncation=15; prewhite=FALSE; adjust=FALSE; center=TRUE
#> p-value = 0.2449 (closed-form law)
#> reference: chi-squared law with 2 degrees of freedom

Intervals, contrasts and regions

Every method supports the same inference objects. The three interval types keep their meaning: "simultaneous" projects the joint region and holds for all linear contrasts at once, "joint" treats the selected contrasts as a lower-dimensional problem, and "marginal" builds a separate one-dimensional interval per coordinate without simultaneous coverage.

confint(fit, method = "ewc")
#> Simultaneous (joint-region projection) Equal-weighted cosine (EWC) confidence intervals, level 0.95
#>      lower  upper
#> m1 -0.1143 0.4049
#> m2 -0.2943 0.3613
#> critical value 7.335 (reference dimension 2); scaled F law: 2 nu/(nu - 2 + 1) times F(2, 20) with nu = 21 cosine terms
aersn_contrast(fit, rbind("m1 - m2" = c(1, -1)), method = "fixedb", b = 0.5,
               draws = 2000, seed = 1)
#> Simultaneous (joint-region projection) Bartlett fixed-b intervals for linear contrasts, level 0.95
#>         estimate   lower  upper
#> m1 - m2   0.1118 -0.6603 0.8838
#> critical value 25.32 (reference dimension 2); matched-grid Monte Carlo law for Bartlett fixed-b: q = 2, uniform grid with n = 400 intervals, 2000 draws, seed 1, b = 0.500000 (bandwidth in grid intervals 200)

Inference on a target is carried out by rebuilding the method on the target’s influence contributions. For the hull, and for any method whose normalizer is a linear functional of outer products, this agrees with transforming the full-dimensional normalizer. It differs where the method itself is not linear in that sense: the LDL factorization is recomputed for the target, and an automatic HAC bandwidth is selected again from the target’s contributions.

Confidence regions keep the geometry of their method. The hull region is a convex polygon in two dimensions; the other five are ellipsoids.

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(aersn_region(fit, method = "hull", draws = 2000, seed = 1),
     null = c(0, 0), main = "Increment hull")
plot(aersn_region(fit, method = "shao", draws = 2000, seed = 1),
     null = c(0, 0), main = "Quadratic self-normalization")

par(op)

Changing the tuning values

Four methods need a tuning value. The package reports the value it used rather than only the value requested, which matters when a rule selects from the data or when a fraction is rounded to an integer lag.

HAC kernel and bandwidth

Three kernels are available, with five ways to set the bandwidth.

grid <- expand.grid(kernel = c("Bartlett", "Parzen", "Quadratic Spectral"),
                    rule = c("short", "long", "andrews", "newey-west"),
                    stringsAsFactors = FALSE)
for (i in seq_len(nrow(grid))) {
  out <- tryCatch({
    nz <- aersn_hac_lrv(fit, kernel = grid$kernel[i],
                        bandwidth = grid$rule[i])
    sprintf("bandwidth %6.3f, lag %s", nz$tuning$bandwidth,
            nz$tuning$lag_truncation)
  }, error = function(e) "not supported")
  cat(sprintf("%-20s %-12s %s\n", grid$kernel[i], grid$rule[i], out))
}
#> Bartlett             short        bandwidth  6.000, lag 5
#> Parzen               short        bandwidth  6.000, lag 5
#> Quadratic Spectral   short        bandwidth  6.000, lag 5
#> Bartlett             long         bandwidth 27.000, lag 26
#> Parzen               long         bandwidth 27.000, lag 26
#> Quadratic Spectral   long         bandwidth 27.000, lag 26
#> Bartlett             andrews      bandwidth 10.346, lag 10
#> Parzen               andrews      bandwidth 15.620, lag 15
#> Quadratic Spectral   andrews      bandwidth  7.759, lag NA
#> Bartlett             newey-west   bandwidth 10.000, lag 9
#> Parzen               newey-west   not supported
#> Quadratic Spectral   newey-west   not supported

The Newey-West rule follows sandwich::NeweyWest(), which applies Bartlett weights with bandwidth floor(bw) + 1; it is therefore offered for the Bartlett kernel only, and the other kernels report that rather than quietly substituting a different rule. A numeric bandwidth, or a lag truncation for the two compactly supported kernels, can be given directly:

aersn_hac_lrv(fit, kernel = "Bartlett", lag = 8)$tuning[c("bandwidth",
                                                          "lag_truncation")]
#> $bandwidth
#> [1] 9
#> 
#> $lag_truncation
#> [1] 8

Three quantities are easy to confuse and are kept distinct: the real-valued bandwidth h in the weight w(lag/h); the lag truncation L, which for the Bartlett kernel corresponds to h = L + 1; and the fixed-b fraction, which is a third parameterisation.

The fixed-b fraction

The bandwidth in grid intervals must be a whole number, so the fraction actually used is round(b n) / n:

for (b in c(0.2, 0.5, 0.7, 0.9, 1)) {
  nz <- aersn_fixed_b_normalizer(fit, b = b)
  cat(sprintf("requested b = %.2f -> m = %3d, realized b = %.4f\n",
              b, nz$tuning$m, nz$tuning$b_grid))
}
#> requested b = 0.20 -> m =  80, realized b = 0.2000
#> requested b = 0.50 -> m = 200, realized b = 0.5000
#> requested b = 0.70 -> m = 280, realized b = 0.7000
#> requested b = 0.90 -> m = 360, realized b = 0.9000
#> requested b = 1.00 -> m = 400, realized b = 1.0000

Each fraction has its own reference law, keyed by the realized value, so a law simulated for one b cannot be used with another. At b = 1 the estimate is exactly twice the quadratic self-normalizer, the statistic is exactly half the quadratic statistic, and the two tests agree:

max(abs(aersn_fixed_b_normalizer(fit, b = 1)$matrix -
          2 * aersn_shao_normalizer(fit)$matrix))
#> [1] 1.776357e-15

The number of cosine terms

aersn_ewc_lrv(fit)$tuning[c("nu", "rule")]
#> $nu
#> [1] 21
#> 
#> $rule
#> [1] "floor(0.4 n^(2/3))"
aersn_test(fit, method = "ewc", nu = 30)$critical.value
#> [1] 6.884802

The default is floor(0.4 n^(2/3)). Admissible values run from the parameter dimension to n - 1. Raising the number of terms lowers the critical value, because the reference law moves towards chi-squared, and raises the bias of the estimate under strong dependence.

nu has its own argument in aersn_test() and the other inference functions. Without it, R’s partial matching would send nu = 30 to the null argument, since nu is a prefix of null.

What each method assumes

All six need the influence contributions to satisfy a functional central limit theorem with a nonsingular long-run covariance matrix, and the estimator to be asymptotically linear in them. Beyond that:

The hull needs no further condition and its statistic is unchanged by any nonsingular reparameterization of the parameters.
The componentwise method uses a reference law that treats its q ratios as independent. Diagonalizing the sample lag-zero covariance does not diagonalize the long-run covariance, so that law is correct only when the transformed long-run covariance is diagonal in the limit. The statistic also depends on the order of the coordinates and is not affine equivariant.
Quadratic self-normalization integrates in calendar time by default. Under nonuniform variance accumulation the integration should follow the variance-accumulation profile, which integration = "profile" does.
HAC inference relies on the bandwidth conditions under which the estimate is consistent; a particular finite-sample bandwidth does not verify them.
Fixed-b and EWC reference laws are fixed-smoothing asymptotic laws for dependent data, not exact finite-sample distributions.

The package checks what it can check, such as the dimension, the grid, the conditioning of a normalizer and the match between a statistic and its reference law. It cannot check the dependence conditions, and none of the methods removes the size distortion that all of them show under strong persistence; the manuscript’s simulations report it for each.

aersn_normalizer(fit, "ldl")$notes
#> [1] "Diagonalizing the sample lag-zero covariance does not imply a diagonal long-run covariance; the independent-component reference law needs that additional condition."
#> [2] "Not affine equivariant: the value depends on coordinate order and on nonsingular reparameterization."