Replicated Lugsail Batch Means Estimator

• Replicated lugsail batch means estimator is a technique that combines lugsail lag-window adjustments with batch means and replication to improve covariance estimation in time series.
• It modifies standard weighted batch means through an inflation factor and weight to effectively counteract negative bias, even under strong positive autocorrelation.
• Replication across independent runs reduces variance by a factor of 1/R, making the method robust for simulation studies and MCMC output analysis.

The replicated lugsail batch means estimator is a methodology for the estimation of time-average covariance matrices of stationary stochastic processes, including those encountered in Markov chain Monte Carlo (@@@@2@@@@) output analysis and simulation studies. This framework integrates a lugsail transformation into traditional lag windows, weighted batch means (BM) estimators, and replication across independent runs to achieve superior finite-sample bias and variance properties, especially in the presence of strong positive autocorrelation (Vats et al., 2018).

1. Lugsail Lag-Window Construction

The lugsail lag window modifies any symmetric base kernel or lag window $g:\mathbb{R} \to \mathbb{R}$ with $g(0) = 1$, producing a new family of windows that address negative bias in spectral variance and covariance estimation. Given an inflation factor $\beta \geq 1$ and a weight $c \in [0, 1)$ (potentially sequence-dependent as $c_n$ with $c_n \to c$), the lugsail window is defined as:

$$g_L(x) = \frac{1}{1-c} g(x) - \frac{c}{1-c} g(\beta x).$$

Key properties:

• $g_L(0) = 1$.
• For $c > 1/\beta^q$ (where $q$ is the bias order for the base kernel) and original BM bias $O(b^{-q}) < 0$, $g_L$ induces a positive $O(b^{-q})$ first-order bias to offset negative bias.
• The window can exceed 1 initially, yielding a "lugsail" shape.

2. Weighted Batch Means Estimation with Lugsail Windows

Let $\{Y_t\}$ be a (possibly $p$-dimensional) stationary time series of length $n$. Choosing a maximal batch size $b$ ($1 \le b < n$), the construction proceeds as follows:

• Let $a_s = \lfloor n / s \rfloor$ be the number of non-overlapping batches for batch size $s$.
• For batch index $l = 0, ..., a_s-1$, the $l$-th batch mean is

$$\bar{Y}_l(s) = \frac{1}{s} \sum_{t=ls+1}^{(l+1)s} Y_t$$

• The corresponding second-difference weights:

$$\Delta_2^{(L)}(s) = g_L\left( \frac{s-1}{b} \right) - 2 g_L\left( \frac{s}{b} \right) + g_L\left( \frac{s+1}{b} \right)$$

The lugsail weighted BM estimator is then

$$\hat{\Sigma}_L = \sum_{s=1}^b \frac{s^2 \Delta_2^{(L)}(s)}{a_s - 1} \sum_{l=0}^{a_s - 1} (\bar{Y}_l(s) - \bar{Y})(\bar{Y}_l(s) - \bar{Y})^T$$

A jackknife-style equivalent is:

$$\hat{\Sigma}_L = \frac{1}{1-c} \hat{\Sigma}_b - \frac{c}{1-c} \hat{\Sigma}_{b/\beta}$$

where $\hat{\Sigma}_s$ is the standard BM estimator for batch size $s$.

3. Replicated Estimation Across Independent Runs

For multiple independent realizations (e.g., $R$ independent MCMC chains, each of length $n$), let $\hat{\Sigma}_L^{(r)}$ denote the lugsail BM estimate from the $r$-th run. The replicated lugsail estimator is:

$$\hat{\Sigma}_{RL} = \frac{1}{R} \sum_{r=1}^R \hat{\Sigma}_L^{(r)}$$

Empirical variance of the replicated estimator across $R$ chains is:

$$\widehat{\text{Var}}(\hat{\Sigma}_{RL}) = \frac{1}{R(R-1)} \sum_{r=1}^R \left( \hat{\Sigma}_L^{(r)} - \hat{\Sigma}_{RL} \right)\left( \hat{\Sigma}_L^{(r)} - \hat{\Sigma}_{RL} \right)^T$$

Replication substantially reduces estimator variance by a factor of $1/R$, with negligible increase in bias for large $R$.

4. Bias and Variance Under $\alpha$-Mixing Conditions

Under $\alpha$-mixing (strong mixing) with bounded fourth moments (Assumption A.1 of Vats & Flegal (2019)), if the ordinary BM exhibits first-order bias $O(b^{-1})$ and one chooses $b \to \infty$, $b/n \to 0$, the bias and variance for lugsail BM are:

• Bias:

$$\mathbb{E}\left[ \hat{\Sigma}_L \right] - \Sigma = \frac{\Gamma}{b} \frac{1 - \beta c}{1 - c} + o(b^{-1})$$

where $\Gamma = -\sum_{s=1}^\infty s \{ R(s) + R(s)^T \}$, with $R(s)$ the autocovariance function.

• Variance (for entry $(i,j)$):

$$\text{Var}(\hat{\Sigma}_L^{ij}) = \left[ \frac{1}{\beta} + \frac{\beta - 1}{\beta (1 - c)^2} \right] (\Sigma_{ii}\Sigma_{jj} + \Sigma_{ij}^2)\frac{b}{n} + o(b/n)$$

For the replicated estimator:

• Bias remains as above.
• Variance is reduced by $1/R$.

5. Recommended Parameter Choices and Practical Guidance

Parameter selection is guided by the underlying process's autocorrelation structure. Vats & Flegal (2019) prescribe rules-of-thumb:

• Moderate correlation ($\rho \in (0,0.7)$): Zero-lugsail with $\beta = 2$, $c = \beta^{-q}$ (zero first-order bias).
• Moderate to high correlation ($\rho \in (0.7,0.95)$): Adaptive lugsail with $\beta = 2$, $$c_n = \frac{\log n - \log b + 1}{\beta^q(\log n - \log b)+1} \to \beta^{-q}$$.
• High to extreme correlation ($\rho \in (0.95,1)$): Over-lugsail with $\beta = 3$, $c = \frac{2}{1 + 3^q}$, introducing a small positive $O(b^{-q})$ bias to counteract large negative higher-order terms.
• Batch size recommendation: $b = \lfloor n^{1/2} \rfloor$.

Replication ($R > 1$) is strongly recommended when computational resources allow, as it reduces estimator variance by $1/R$ with no substantial impact on bias.

6. Context, Applications, and Extension

The replicated lugsail batch means estimator is designed for broad applicability in MCMC, steady-state simulation, and time series analysis, particularly in regimes with substantial positive serial dependence. Its ability to convert any existing lag window to a lugsail form without new assumptions, and its compatibility with ordinary batch means estimators, makes it particularly practical for large-scale simulation studies (Vats et al., 2018). The approach is supported by theoretical results on bias and variance, substantially weakening mixing requirements for consistency.

The lugsail methodology has demonstrated improved finite-sample reliability over conventional kernels in vector autoregressive models and Bayesian logistic regression, where conventional methods suffer from severe negative bias.

7. Summary Table: Parameter Choices by Correlation Regime

Correlation Type	Lugsail Parameters	Bias Effect
Moderate ($\rho<0.7$)	$\beta=2$, $c=2^{-q}$	Zero first-order bias
Moderate–High ($0.7<\rho<0.95$)	$\beta=2$, $c_n$ adaptive	Bias $\rightarrow$ zero as $n\uparrow$
High–Extreme ($\rho>0.95$)	$\beta=3$, $c=2/(1+3^q)$	Small positive first-order bias

This estimator has become a preferred tool for practitioners requiring accurate time-average covariance estimation from highly correlated simulation or MCMC data, especially when computational efficiency and consistency under weak mixing are essential (Vats et al., 2018).