# Multivariate moment least-squares estimators for reversible Markov chains

Hyebin Song and Stephen Berg11footnotemark: 1
Department of Statistics, Pennsylvania State University Both authors
contributed equally. The authors gratefully acknowledge NSF support
DMS-2311141.Corresponding author<EMAIL_ADDRESS>
###### Abstract

Markov chain Monte Carlo (MCMC) is a commonly used method for approximating
expectations with respect to probability distributions. Uncertainty assessment
for MCMC estimators is essential in practical applications. Moreover, for
multivariate functions of a Markov chain, it is important to estimate not only
the auto-correlation for each component but also to estimate cross-
correlations, in order to better assess sample quality, improve estimates of
effective sample size, and use more effective stopping rules. Berg and Song
(2022) introduced the moment least squares (momentLS) estimator, a shape-
constrained estimator for the autocovariance sequence from a reversible Markov
chain, for univariate functions of the Markov chain. Based on this sequence
estimator, they proposed an estimator of the asymptotic variance of the sample
mean from MCMC samples. In this study, we propose novel autocovariance
sequence and asymptotic variance estimators for Markov chain functions with
multiple components, based on the univariate momentLS estimators from Berg and
Song (2022). We demonstrate strong consistency of the proposed
auto(cross)-covariance sequence and asymptotic variance matrix estimators. We
conduct empirical comparisons of our method with other state-of-the-art
approaches on simulated and real-data examples, using popular samplers
including the random-walk Metropolis sampler and the No-U-Turn sampler from
STAN.

Keywords: Markov chain Monte Carlo, shape-constrained inference, cross-
covariance sequence estimation, asymptotic variance

##  1 Introduction

In modern statistical applications, it is often necessary to analyze
intractable probability distributions, for which computing key summary
statistics such as expectations or quantiles analytically is challenging.
Markov chain Monte Carlo (MCMC) methods offer a practical approach for drawing
samples from these intractable distributions. An MCMC experiment involves
generating a sequence of random variables $X=(X_{0},X_{1},\dots)$, such that
the distribution of $X_{t}$ converges to the intractable target distribution
$\pi$. This simulated chain $X$ can then be used to estimate the quantity of
interest, such as the expectation with regard to the target distribution. For
instance, suppose we want to approximate $\mu_{g}=\int g(x)\pi(dx)$ where
$\pi$ is intractable. One can generate a Markov chain $X$ of length $M$ with
invariant distribution $\pi$, and use the MCMC estimator
$\hat{\mu}_{gM}=M^{-1}\sum_{t=0}^{M-1}g(X_{t})$ to estimate $\mu_{g}$.

The accuracy of an MCMC estimator can be quantified by its asymptotic
variance. Suppose we have a function $g:\mathsf{X}\to\mathbb{R}^{d}$ where
$d\geq 1$. Under mild conditions, a Markov chain central limit theorem holds
(e.g., Meyn and Tweedie, 2009), yielding

$\displaystyle\sqrt{M}(\hat{\mu}_{gM}-\mu_{g})\overset{d}{\to}\mathcal{N}(0,\Sigma).$
(1)

Due to correlations in $(X_{t})_{t\geq 0}$, $\Sigma$ usually differs from
$\textrm{Var}_{\pi}(g(X_{0}))$. In fact, $\Sigma$ aggregates the covariances
across all lags, i.e.,

$\displaystyle\Sigma=\sum_{k=-\infty}^{\infty}\textrm{Cov}_{\pi}(g(X_{0}),g(X_{k})).$
(2)

where we define
$\textrm{Cov}_{\pi}(g(X_{0}),g(X_{k}))=\textrm{Cov}_{\pi}(g(X_{|k|}),g(X_{0}))$
for $k<0$.

Historically, the assessment of uncertainty in MCMC estimates has often been
conducted on a component-wise basis. This amounts to considering univariate
CLT results
$\sqrt{M}(\hat{\mu}_{gM,j}-\mu_{gj})\overset{d}{\to}\mathcal{N}(0,\Sigma_{jj})$
for $j=1,\dots,d$, and therefore focusing only on the diagonals of the
asymptotic variance $\Sigma$. In particular, the quantity
$\sqrt{\Sigma_{jj}/M}$ is called the MCMC standard error of
$\hat{\mu}_{gM,j}$. Estimates of these MCMC standard errors have been used to
determine effective sample sizes, construct confidence intervals, and
implement stopping criteria to achieve a predefined level of precision for
each component in the Markov chain (Heidelberger and Welch, 1981; Jones et
al., 2006).

However, by accounting for dependencies among the components of functions of a
Markov chain, we can gain a more comprehensive understanding of the overall
quality of the Markov chain estimates and formulate a more effective stopping
criterion (Dai and Jones, 2017; Vats et al., 2019). For instance, using an
estimate of the asymptotic variance matrix $\Sigma$, we can construct
confidence regions for $\mu_{g}$ with smaller volumes compared to Bonferroni-
corrected confidence intervals. Additionally, an effective sample size can be
computed to summarize the overall inflation or deflation of sample size
relative to an independent and identically distributed (iid) chain. This
approach allows for the formulation of stopping rules based on an overall
accuracy criterion and offers practical advantages over a component-wise
criterion that terminates based on the worst-performing component, which might
be overly conservative.

Moreover, the estimation of the entire asymptotic variance matrix $\Sigma$ for
a multivariate chain, as opposed to only estimating its diagonals, is
essential in quantifying uncertainties in linear combinations of
$\hat{\mu}_{gM,j},j=1,\dots,d$. This is especially useful in devising MCMC
estimators with reduced variance using control variate methods (see, e.g.
Dellaportas and Kontoyiannis, 2012; Berg et al., 2019; Belomestny et al.,
2020; South et al., 2022). Control variate methods employ estimators of the
form
$\tilde{\mu}_{gM,C}=M^{-1}\sum_{t=0}^{M-1}\\{g(X_{t})-C^{\top}h(X_{t})\\}$,
where $h:\mathsf{X}\to\mathbb{R}^{p}$ is a control variate function chosen to
satisfy $\int h(x)\pi(dx)=0$, and $C\in\mathbb{R}^{p\times d}$. Since $h$ has
$\pi$ expectation 0, the terms $C^{\top}h(X_{t})$ do not lead to bias under
stationarity. If $C$ is chosen such that the asymptotic variance of
$\tilde{\mu}_{gM,C}$ is smaller than that of $\hat{\mu}_{gM}$, then it is
preferable to use $\tilde{\mu}_{gM,C}$ over $\hat{\mu}_{gM}$ to estimate
$\mu_{g}$. For example, Belomestny et al. (2020) utilize a spectral variance
estimator and minimize the estimated asymptotic variance of the control
variate estimator, and find improvements over alternative control variate
approaches that neglect the auto(cross)-covariance structure of the Markov
chain.

In terms of estimating $\Sigma$, given (2), one might consider estimating each
$\gamma_{g}(k)=\textrm{Cov}_{\pi}(g(X_{0}),g(X_{k}))$ first, and then using a
plug-in estimator $\hat{\Sigma}=\sum_{k=-\infty}^{\infty}\hat{\gamma}_{g}(k)$
to estimate $\Sigma$. A natural estimate for $\gamma_{g}(k)$ is the empirical
covariance matrix $r_{{g},M}(k)\in\mathbb{R}^{d\times d}$, defined as

$\displaystyle
r_{{g},M}(k)=\begin{cases}\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}(X_{t})\tilde{g}(X_{t+{|k|}})^{\top}&0\leq
k<M\\\
\\{\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}(X_{t})\tilde{g}(X_{t+|k|})^{\top}\\}^{\top}&-M<k<0\\\
0&|k|\geq M.\end{cases}$ (3)

where we define
$\tilde{g}(X_{t})=g(X_{t})-\frac{1}{M}\sum_{t=0}^{M-1}g(X_{t})$. Under mild
conditions, a Strong Law of Large Numbers (SLLN) holds such that
$r_{{g},M}(k)$ converges to $\gamma_{g}(k)$ for each fixed $k$. However, this
convergence does not imply the convergence of
$\hat{\Sigma}=\sum_{k=-\infty}^{\infty}r_{{g},M}(k)$ to
$\Sigma=\sum_{k=-\infty}^{\infty}\gamma_{g}(k)$ (e.g., Anderson, 1971). As a
result, alternative methods for estimating $\Sigma$ have been investigated in
previous literature.

In univariate settings where $d=1$, methods for estimating
$\Sigma\in\mathbb{R}_{+}$ with better statistical properties have been
proposed. These methods include spectral variance estimators (Anderson, 1971),
estimators based on batch means (Priestley, 1981; Flegal and Jones, 2010;
Chakraborty et al., 2022), and a class of methods for reversible Markov chains
including initial sequence estimators (Geyer, 1992) and the moment least
squares estimator (Berg and Song, 2022). The spectral variance estimators
utilize windowed autocovariance sequences as an input, where a window function
is often chosen to downweight the contribution from $r_{{g},M}(k)$ for larger
$k$ to the asymptotic variance estimate. The batch means and overlapping batch
means methods utilize scaled sample variances calculated from the sample means
of either non-overlapping or overlapping batches. The batch means and
overlapping batch means estimators turn out to be closely related to the
spectral variance estimators (Damerdji, 1991; Flegal and Jones, 2010). The
batch means and spectral variance estimators have known consistency
properties. In particular, with an appropriate choice of batch or window size,
they are a.s. consistent for $\Sigma$ (Flegal and Jones, 2010).

Geyer (1992), on the other hand, proposed initial sequence estimators for
estimating $\Sigma$ when $d=1$ by imposing various shape constraints on the
estimates of certain sums of autocovariances. More specifically, these
estimators exploit positivity, monotonicity, and convexity constraints
satisfied for reversible Markov chains by the sequence
$\Gamma_{g}=\\{\Gamma_{g}(k)\\}_{k=0}^{\infty}$ defined by
$\Gamma_{g}(k):=\gamma_{g}(2k)+\gamma_{g}(2k+1)$ for $k=0,1,2,\dots$. Initial
sequence estimators have very strong empirical performance and do not require
the choice of a tuning parameter value, making them very useful in practice.
For example, the widely used Stan software (Stan Development Team, 2019)
employs the initial sequence estimators to estimate the effective sample size
of Markov chain simulations. However, the consistency of these initial
sequence estimators of $\Sigma$ is still unknown. To our knowledge, the only
consistency results for the initial sequence estimators are that the initial
sequence estimates asymptotically do not underestimate $\Sigma$.

berg2022efficient further explored the idea of estimating the autocovariance
$\gamma_{g}$ by enforcing shape constraints and proposed a moment least
squares (momentLS) estimator for estimating $\gamma_{g}$ from a reversible and
geometrically ergodic Markov chain. Using the observation that $\gamma_{g}$
from a reversible Markov chain has a moment representation (see Section 2) and
thus satisfies certain shape constraints, they considered the projection of an
input autocovariance sequence estimate onto the set of $[-1+\delta,1-\delta]$
moment sequences, which they called a momentLS sequence estimator. Based on
this autocovariance sequence estimator, they also proposed an asymptotic
variance estimator to estimate the asymptotic variance of the sample mean of
univariate functions of the Markov chain iterates. They showed that both the
autocovariance sequence and asymptotic variance estimator are strongly
consistent, provided that the hyperparameter $\delta$ is chosen sufficiently
small.

Estimating the asymptotic variance matrix in multivariate settings where $d>1$
has also been considered. Many of the estimators can be considered as
multivariate generalizations of univariate estimators. For example,
multivariate counterparts to the spectral variance estimator and batch means
estimator were considered and their theoretical properties were studied in
vats2018strong; vats2019multivariate. A strong consistency property for the
multivariate spectral variance estimator and the strong consistency of
multivariate batch means with an appropriate choice of window function and
batch size were established in vats2018strong and vats2019multivariate,
respectively.

Various generalizations of the initial sequence estimators in
geyer1992practical have been proposed in previous literature. For example,
kosorok2000error constructed an initial positive (matrix) sequence estimator
by truncating the cumulative autocovariance matrix at the first time point
where the smallest eigenvalue becomes negative. dai2017multivariate proposed a
different truncation point where they suggested truncating when the
generalized variance associated with the current cumulative sum fails to
increase. Similar to initial sequence estimators in the univariate setting, no
consistency results have been established for the multivariate initial
sequence estimators of the asymptotic variance, although they are known to be
asymptotically conservative.

In this work, we propose a multivariate sequence estimator for an auto- or
cross-covariance sequence associated with a multivariate function
$g:\mathsf{X}\to\mathbb{R}^{d}$ with $d\geq 2$ of a Markov chain based on the
work of berg2022efficient. We show that our proposed multivariate sequence
estimator is strongly consistent and that the resulting multivariate matrix
estimator for the asymptotic variance matrix is also strongly consistent,
provided an appropriate choice of the hyperparameter $\delta$. The paper is
organized as follows: In Section 2, we describe our problem set-up and provide
some background on univariate momentLS estimators. In Section 3, we introduce
the proposed sequence and asymptotic variance matrix estimators for
multivariate functions of Markov chain iterates. The statistical guarantees of
the proposed estimators are studied in Section 4. In Section 5, we compare the
proposed estimator with competing methods using simulated and real data
examples.

##  2 Problem setup and univariate MomentLS estimator

We consider a $\psi$-irreducible, aperiodic Markov chain
$X=(X_{0},X_{1},\dots)$, which is a sequence of $\mathsf{X}$-valued random
variables. We assume that the state space $\mathsf{X}$ is a complete separable
metric space and let $\mathscr{X}$ be the associated Borel $\sigma$-algebra.
We let $Q:\mathsf{X}\times\mathscr{X}\to[0,1]$ be the transition kernel for
$X$. An initial measure $\nu$ on $\mathscr{X}$ and a transition kernel $Q$
define a Markov chain probability measure $P_{\nu}$ for
$X=(X_{0},X_{1},X_{2},\dots)$ on the canonical sequence space
$(\Omega,\mathcal{F})$. We write $E_{\nu}$ to denote expectation with respect
to $P_{\nu}$.

For a probability measure $\pi$ on $(\mathsf{X},\mathscr{X})$, a probability
kernel $Q$ is said to be $\pi$-stationary if $\pi(A)=\int Q(x,A)\pi(dx)$ for
all $A\in\mathscr{X}$. In addition, we say $Q$ satisfies the reversibility
property with respect to $\pi$ if
$\int_{A}\pi(dx)Q(x,B)=\int_{B}\pi(dy)Q(y,A)$ for all $A,B\in\mathscr{X}$. We
note that reversibility with respect to $\pi$ is a sufficient condition for
$\pi$-stationarity of $Q$. For a function $f:\mathsf{X}\to\mathbb{C}$ and a
transition kernel $Q$, we define the Markov operator $Q$ by

$\displaystyle Qf(x)=\int Q(x,dy)f(y)$ (4)

We define $Q^{0}f(x)=f(x)$, $Q^{1}f(x)=Qf(x)$, and $Q^{t}f(x)=Q(Q^{t-1}f)(x)$
for $t>1$, and we define $Q^{t}(x,A)=Q^{t}1_{A}(x)$, where $1_{A}(\cdot)$ is
the indicator function for the set $A$.

We let $L^{2}(\pi)$ be the space of functions which are square integrable with
respect to $\pi$, i.e.,
$L^{2}(\pi)=\\{f:\mathsf{X}\to\mathbb{C};\int|f(x)|^{2}\pi(dx)<\infty\\}$. For
functions $f,g\in L^{2}(\pi)$, we define an inner product
$\braket{f,g}_{\pi}=\int f(x)\overline{g(x)}\pi(dx).$ where $\overline{g(x)}$
is the conjugate of $g(x)$. We note that $L^{2}(\pi)$ is a Hilbert space
equipped with the inner product $\langle\cdot,\cdot\rangle_{\pi}$. For $f\in
L^{2}(\pi)$, we define $\|f\|_{L^{2}(\pi)}=\sqrt{\braket{f,f}_{\pi}}$. Also,
for an operator $T$ on $L^{2}(\pi)$, we define
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|T\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}=\sup_{f;\|f\|_{L^{2}(\pi)}\leq
1}\|Tf\|_{L^{2}(\pi)}$ and we say $T$ is bounded if
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|T\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}<\infty$.
We say $T^{*}$ is the adjoint of $T$ if $\langle Tf,g\rangle_{\pi}=\langle
f,T^{*}g\rangle_{\pi}$ for any $f,g\in L^{2}(\pi)$, and $T$ is self-adjoint if
$T=T^{*}$.

For a Markov kernel $Q$, we note that $Q$ is a linear operator since
$Q\\{af+bg\\}(x)=\int Q(x,dy)\\{af+bg\\}(y)=a\int Q(x,dy)f(y)+b\int
Q(x,dy)g(y)=aQf(x)+bQg(x)$ for any $a,b\in\mathbb{C}$ and $f,g\in L^{2}(\pi)$.
Since $Q$ is a contraction map,
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Q\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}\leq
1$ and $Q$ is bounded. If the transition kernel $Q$ is reversible with respect
to $\pi$, the Markov operator $Q$ is self-adjoint since $\langle
Q1_{A},1_{B}\rangle_{\pi}=\langle 1_{A},Q1_{B}\rangle_{\pi}$ for any
$A,B\in\mathscr{X}$, and simple functions are dense in $L^{2}(\pi)$.

The spectrum of the operator $Q$ plays a key role in determining the mixing
properties of a Markov chain with transition kernel $Q$. For an operator $T$
on the Hilbert space $L^{2}(\pi)$, the spectrum of $T$ is defined as
$\sigma(T)=\\{\lambda\in\mathbb{C};(T-\lambda I)^{-1}\mbox{does not exist or
is unbounded}\\}.$ For Markov operators $Q$, we define the spectral gap
$\delta(Q)$ of $Q$ as
$\delta(Q)=1-\sup\\{|\lambda|;\lambda\in\sigma(Q_{0})\\}$ where $Q_{0}$ is
defined as $Q_{0}f=Qf-E_{\pi}[f(X_{0})]f_{0}$ and $f_{0}\in L^{2}(\pi)$ is the
constant function such that $f_{0}(x)=1$ for all $x\in\mathsf{X}$. It is easy
to check that $Q_{0}$ is also self-adjoint when $Q$ is self-adjoint. If $Q$ is
reversible, $Q$ has a positive spectral gap $(\delta(Q)>0)$ if and only if the
chain is geometrically ergodic (Roberts1997-qk; kontoyiannis2012geometric).

We let the function of interest $g:\mathsf{X}\to\mathbb{R}^{d}$ be an
$\mathbb{R}^{d}$-valued function such that
$g(x)=[g_{1}(x),\dots,g_{d}(x)]^{\top}$ for each $x\in\mathsf{X}$ and each
$g_{j}$ is square integrable with respect to $\pi$. For any real-valued
$f,g:\mathsf{X}\to\mathbb{R}$, we define the cross-covariance sequence
$\gamma_{fg}(k)$ by

$\displaystyle\gamma_{fg}(k)$ $\displaystyle={\rm
Cov}_{\pi}\\{f(X_{0}),g(X_{k})\\}\quad\mbox{for }k\in\mathbb{Z}.$ (5)

When $f=g$, $\gamma_{ff}(k)$ is the autocovariance sequence for
$\\{f(X_{t})\\}_{t\geq 0}$ and we simply write it as $\gamma_{f}(k)$. We use
$\gamma_{fg}=\\{\gamma_{fg}(k)\\}_{k\in\mathbb{Z}}$ to denote the
auto(cross)-covariance sequence on $\mathbb{Z}$.

In the remainder, we consider a discrete time Markov chain
$X=\\{X_{t}\\}_{t=0}^{\infty}$ with stationary distribution $\pi$ and
$\pi$-reversible transition kernel $Q$ with a positive spectral gap. We
summarize our assumptions on the Markov chain $X$ as follows for future
reference:

  1. (A.1)

(Harris ergodicity) $X$ is $\psi$-irreducible, aperiodic, and positive Harris
recurrent.

  2. (A.2)

(Reversibility) $Q$ is $\pi$-reversible for a probability measure $\pi$ on
$(\mathsf{X},\mathscr{X})$.

  3. (A.3)

(Geometric ergodicity) There exists a real number $\rho<1$ and a non-negative
function $M$ on the state space $\mathsf{X}$ such that
$\|Q^{n}(x,\cdot)-\pi(\cdot)\|_{\rm TV}\leq M(x)\rho^{n},\mbox{ for all
}x\in\mathsf{X},$ where $\|\cdot\|_{\rm TV}$ is the total variation norm.

Throughout the paper, we assume that each component of
$g:\mathsf{X}\to\mathbb{R}^{d}$ is in $L^{2}(\pi)$, i.e,

  1. (B.1)

(Square integrability)$\int g_{j}(x)^{2}\pi(dx)<\infty$, for $j=1,\dots,d$.

For the definitions of $\psi$-irreducibility, aperiodicity, and positive
Harris recurrence, see e.g., meyn2009markov. Reversibility and geometric
ergodicity are the two key assumptions for the class of Markov chains we
consider in this paper. While reversibility is not necessary to obtain a
$\pi$-stationary chain, many popular Markov chains possess this property,
including Metropolis-Hastings and random scan Gibbs samplers. Geometric
ergodicity implies exponential convergence of the Markov chain $X$ to its
target distribution $\pi$. When the state space $\mathsf{X}$ is finite, all
irreducible and aperiodic Markov chains are geometrically ergodic. While
geometric ergodicity is not guaranteed for infinite state spaces, geometric
ergodicity is still of practical and theoretical importance. For instance,
geometric ergodicity is desirable for Markov chain samplers because non-
geometrically ergodic chains often exhibit slow convergence
(roberts1998markov). In addition, geometric ergodicity provides one of the
simplest sufficient conditions under which a Markov chain CLT holds.
Specifically, for a reversible and geometrically ergodic chain, the existence
of a finite second moment of the function of interest $g$ is sufficient to
establish a CLT (e.g., jones2004central)

$\displaystyle\sqrt{M}(\hat{\mu}_{gM}-\mu_{g})\to\mathcal{N}(0,\Sigma),$ (6)

where $\hat{\mu}_{gM}=\frac{1}{M}\sum_{t=0}^{M-1}g(X_{t})$ and $\mu_{g}=\int
g(x)\pi(dx)$.

berg2022efficient proposed shape-constrained estimators for both
autocovariance sequence and asymptotic variance estimation for the univariate
case $g:\mathsf{X}\to\mathbb{R}^{d}$ for $d=1$. These estimators are based on
the fact that the autocovariance sequence $\gamma_{g}$ from a reversible
Markov chain is a $[-1,1]$-moment sequence, i.e., there exists a positive
measure $F_{g}$ supported on $[-1,1]$ such that

$\displaystyle\gamma_{g}(k)=\int\alpha^{|k|}\,F_{g}(d\alpha),\;\;\;\;\;\;k\in\mathbb{Z}.$
(7)

In other words, each $\gamma_{g}(k)$ is the $|k|$th moment of $F_{g}$. The
moment representation (7) of $\gamma_{g}$ imposes specific shape restrictions
on $\gamma_{g}$, and to exploit this, berg2022efficient proposed to project a
point-wise consistent initial autocovariance sequence estimator $r_{{g},M}$
onto the set of moment sequences satisfying a representation as in (7).

To be more concrete, let us first introduce some definitions and notation. For
$a,b\in\mathbb{R}$ such that $-\infty<a\leq b<\infty$, we say that a sequence
$m$ is an $[a,b]$-moment sequence if there exists a positive measure $\mu$
supported on $[a,b]$ such that the equation

$\displaystyle m(k)=\int x^{|k|}\mu(dx)$ (8)

holds for any $k\in\mathbb{Z}$. The positive measure $\mu$ is called a
representing measure for the sequence $m$. We also define the moment space
$\mathscr{M}_{\infty}([a,b])$ as the set of $[a,b]$-moment sequences. When
$b>0$ and $a=-b$, we simply write
$\mathscr{M}_{\infty}(b)=\mathscr{M}_{\infty}([-b,b])$. We use
$\ell_{2}(\mathbb{Z})$ to denote the set of square-summable real-valued
sequences on $\mathbb{Z}$, i.e.,
$\ell_{2}(\mathbb{Z})=\\{m\in\mathbb{R}^{\mathbb{Z}};\|m\|_{2}^{2}<\infty\\}$
where we denote $\|m\|_{2}^{2}=\sum_{k\in\mathbb{Z}}m(k)^{2}$.

The moment least squares sequence (momentLS) estimator
$\Pi_{\delta}(r_{{g},M})$ resulting from an initial autocovariance sequence
estimator $r_{{g},M}\in\ell_{2}(\mathbb{Z})$ is defined as the projection of
$r_{{g},M}$ onto the $[-1+\delta,1-\delta]$ moment sequence space, i.e.,

$\displaystyle\Pi_{\delta}(r_{{g},M})=\underset{m\in\mathscr{M}_{\infty}({\delta})\cap\ell_{2}(\mathbb{Z})}{\arg\min}\sum_{k\in\mathbb{Z}}\\{r_{{g},M}(k)-m(k)\\}^{2}.$
(9)

Here $\delta>0$ is a hyperparameter for the momentLS estimator
$\Pi_{\delta}(r_{{g},M})$. We require the hyperparameter $\delta$ to be chosen
sufficiently small such that the true autocovariance $\gamma_{g}$ is a
feasible point for the optimization problem (9). That is, we require ${\rm
Supp}(F_{g})\subseteq[-1+\delta,1-\delta]$, where $F_{g}$ is the representing
measure for $\gamma_{g}$. For a measure $F$ on supported on $[-1,1]$, with the
definition

$\displaystyle\Delta(F)=1-\max\\{|x|;x\in{\rm Supp}(F)\\},$ (10)

the condition ${\rm Supp}(F_{g})\subseteq[-1+\delta,1-\delta]$ on $\delta$ is
equivalent to $0<\delta\leq\Delta(F_{g})$. For a reversible, geometrically
ergodic Markov chain, $\Delta(F_{g})\geq\delta(Q)>0$, and therefore such a
choice of $\delta$ exists. Once the autocovariance sequence $\gamma_{g}$ is
estimated by $\Pi_{\delta}(r_{{g},M})$, berg2022efficient proposed to estimate
the asymptotic variance $\Sigma=\sum_{k\in\mathbb{Z}}\gamma_{g}(k)$ in (6) by
$\hat{\Sigma}=\sum_{k\in\mathbb{Z}}\Pi_{\delta}(r_{{g},M})(k)$.

berg2022efficient showed for univariate $g:\mathsf{X}\to\mathbb{R}$ that if a
valid initial autocovariance sequence estimator $r_{{g},M}$ is used, in the
sense that the initial estimator $r_{{g},M}$ satisfies the following
conditions

  1. (R.1)

(a.s. elementwise convergence)
$r_{{g},M}(k)\underset{M\to\infty}{\to}\gamma_{g}(k)$ for each
$k\in\mathbb{Z}$, $P_{x}$-almost surely, for any initial condition
$x\in\mathsf{X}$,

  2. (R.2)

(finite support) $r_{{g},M}(k)=0$ for $k\geq n(M)$ for some $n(M)<\infty$, and

  3. (R.3)

(even function with a peak at 0) $r_{{g},M}(k)=r_{{g},M}(-k)$ and
$r_{{g},M}(0)\geq|r_{{g},M}(k)|$ for each $k\in\mathbb{Z}$,

and a valid $\delta$ is used, in the sense that $0<\delta\leq\Delta(F_{g})$,
then the momentLS sequence and asymptotic variance estimators are strongly
consistent for $\gamma_{g}$ and $\Sigma$, respectively. For future reference,
we summarize the consistency results for univariate momentLS estimators as
follows:

######  Theorem 1.

Suppose $X_{0},X_{1},...,$ is a Markov chain with transition kernel $Q$
satisfying (A.1)-(A.3), and suppose $g:\mathsf{X}\to\mathbb{R}$ satisfies
(B.1). Let $\gamma_{g}$ denote the autocovariance sequence as defined in (5),
and let $F_{g}$ denote the representing measure for $\gamma_{g}$. Suppose
$\delta$ is chosen so that $0<\delta\leq\Delta(F_{g})$. Let $r_{{g},M}$ be an
initial autocovariance sequence estimator for $\gamma_{g}$ satisfying
conditions (R.1) \- (R.3). Then for each initial condition $x\in\mathsf{X}$,

  1. 1.

$\sum_{k\in\mathbb{Z}}\\{\gamma_{g}(k)-\Pi_{\delta}(r_{{g},M})(k)\\}^{2}\underset{M\to\infty}{\to}0,\,\,\,P_{x}\mbox{-a.s.}$

  2. 2.

$P_{x}(\hat{\mu}_{g,M}\to F_{g}\mbox{ vaguely, as }M\to\infty)=1$, where
$\hat{\mu}_{g,M}$ and $F_{g}$ are the representing measures for
$\Pi_{\delta}(r_{{g},M})$ and $\gamma_{g}$, and

  3. 3.

$\sigma^{2}(\Pi_{\delta}(r_{{g},M}))\to\sigma^{2}(\gamma_{g})\,\,\,P_{x}\mbox{-a.s.}$

where we define $\sigma^{2}(m)=\sum_{k\in\mathbb{Z}}m(k)$ for a sequence $m$
on $\mathbb{Z}$.

##  3 Multivariate MomentLS estimator

We now introduce the multivariate momentLS estimator for the cross-covariance
sequence $\gamma_{g_{i},g_{j}}$ and asymptotic variance $\Sigma$ in the Markov
chain central limit theorem (6). Recall that for a Markov chain
$X=(X_{0},X_{1},\dots)$ and an $\mathbb{R}^{d}$-valued function
$g:\mathsf{X}\to\mathbb{R}^{d}$, we defined the cross (auto)-covariance
sequences $\gamma_{g_{i},g_{j}}$ for $i,j=1,\dots,d$ as

$\displaystyle\gamma_{g_{i},g_{j}}(k)={\rm
Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}))$

for $i=1,\dots,d$ and $k\in\mathbb{Z}$, where we define ${\rm
Cov}_{\pi}\\{g_{i}(X_{0}),g_{j}(X_{k})\\}={\rm
Cov}_{\pi}\\{g_{i}(X_{|k|}),g_{j}(X_{0})\\}$ for $k<0$, and we write
$\gamma_{g_{i}}=\gamma_{g_{i},g_{i}}$.

We first note that $\gamma_{g_{i},g_{j}}(k)=\gamma_{g_{j},g_{i}}(-k)$ for
$k\in\mathbb{Z}$ by definition, as for $k\in\mathbb{N}$,
$\gamma_{g_{j},g_{i}}(-k)={\rm Cov}_{\pi}(g_{j}(X_{|k|}),g_{i}(X_{0}))={\rm
Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}),)=\gamma_{g_{i},g_{j}}(k)$. In addition,
for a reversible Markov chain,
$\gamma_{g_{i},g_{j}}(k)=\gamma_{g_{j},g_{i}}(k)$ for $k\in\mathbb{Z}$. This
is because for $k\geq 0$, we have from self-adjointness of $Q$ that
$\gamma_{g_{i},g_{j}}(k)={\rm
Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}))=\langle\check{g}_{i},Q^{k}\check{g}_{j}\rangle_{\pi}=\langle
Q^{k}\check{g}_{i},\check{g}_{j}\rangle_{\pi}={\rm
Cov}_{\pi}(g_{i}(X_{k}),g_{j}(X_{0}))=\gamma_{g_{j},g_{i}}(k),$ where we
define the true-mean centered version $\check{g}_{i}$ of $g_{i}$ by
$\check{g}_{i}(x)=g_{i}(x)-E_{\pi}[g_{i}(X_{0})]$. For $k<0$, the result
follows from $\gamma_{g_{i},g_{j}}(k)=\gamma_{g_{j},g_{i}}(-k)$.

When $i=j$, i.e., $\gamma_{g_{i},g_{i}}$ is an autocovariance sequence,
$\gamma_{g_{i},g_{i}}$ admits the following moment representation,

$\gamma_{g_{i},g_{i}}(k)=\int\lambda^{|k|}F_{g_{i},g_{i}}(d\lambda)$

for any $k\in\mathbb{Z}$, where $F_{g_{i},g_{i}}$ is a positive measure with
${\rm Supp}(F_{g_{i},g_{i}})\subseteq[-1+\delta(Q),1-\delta(Q)]$, where we
recall that $\delta(Q)$ is the spectral gap of $Q$ (see, e.g.,
geyer1992practical). In other words, $\gamma_{g_{i},g_{i}}$ is a
$[-1+\delta(Q),1-\delta(Q)]$ moment sequence. In fact, $\gamma_{g_{i},g_{j}}$
when $i\neq j$ admits a similar mixture representation, although the
corresponding representing measure $F_{g_{i},g_{j}}$ may not necessarily be
positive. Moreover, all of these (signed) measures, denoted as
$F_{g_{i},g_{j}}$, arise from a spectral resolution of the Markov kernel $Q$.
A notable implication is that all $F_{g_{i},g_{j}}$ are supported on the same
interval $[-1+\delta(Q),1-\delta(Q)]$.

To be more concrete, let us first define a spectral resolution for a Hilbert
space $\mathcal{H}$ equipped with an inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$. Let $\|\cdot\|_{\mathcal{H}}$ be
the norm induced by the inner product. A function
$E:\mathbb{R}\to\mathscr{B}(\mathcal{H})$ is called a spectral resolution on
$\mathcal{H}$ if the function $E$ satisfies the following (see e.g., Chapter 6
in stein2009real): 1. $\|E(\lambda)f\|_{\mathcal{H}}$ is an increasing
function of $\lambda$ for every $f\in\mathcal{H}$, 2. there is an interval
$[a,b]$ such that $E(\lambda)=0$ if $\lambda<a$, and $E(\lambda)=I$ if
$\lambda\geq b$ where $I$ denotes the identity operator on $\mathcal{H}$, and
3. $\lim_{\mu\downarrow\lambda}E(\mu)f=E(\lambda)f$ for any $f\in\mathcal{H}$.
It can be shown that the property 1 is equivalent to the condition that for
any $\lambda<\mu$, $E(\mu)-E(\lambda)$ is an orthogonal projection operator.

For a reversible chain with a Markov kernel $Q$, by the spectral
representation theorem for self-adjoint bounded operators in Hilbert space
(e.g., Theorem 6.1 in stein2009real), we can obtain a spectral resolution $E$
such that for any $f,g\in L_{2}(\pi)$

$\displaystyle\langle Q_{0}f,g\rangle_{\pi}=\int\alpha\,F_{f,g}(d\alpha)$ (11)

where we denote by $F_{f,g}$ the complex-valued measure on
$\mathscr{B}_{\mathbb{R}}$ uniquely defined by the complex-valued function
$\Psi_{f,g}:\alpha\to\langle E(\alpha)f,g\rangle_{\pi}$ through
$F_{f,g}((-\infty,\alpha])=\Psi_{f,g}(\alpha)$. When $f=g$, we let
$F_{f}=F_{f,f}$ for notational simplicity. The following Proposition 1
summarizes the properties of the mixing measure $F_{g_{i},g_{j}}$ for the
cross-covariance sequence $\gamma_{g_{i},g_{j}}$.

######  Proposition 1.

Suppose $X$ satisfies (A.2) and $g:\mathsf{X}\to\mathbb{R}^{d}$ satisfies
(B.1). There exists a finite signed measure $F_{g_{i},g_{j}}$ on
$(\mathbb{R},\mathscr{B}_{\mathbb{R}})$ with ${\rm
Supp}(F_{g_{i},g_{j}})\subseteq[-1+\delta(Q),1-\delta(Q)]$ such that

$\displaystyle\gamma_{g_{i},g_{j}}(k)=\int\lambda^{|k|}\,F_{g_{i},g_{j}}(d\lambda)$
(12)

for each $k\in\mathbb{Z}$ and $i,j\in\\{1,\dots,d\\}$. Moreover, the mixing
measure $F_{g_{i},g_{j}}$ in (12) is unique, and is given by
$F_{g_{i},g_{j}}((-\infty,\alpha])=\langle E(\alpha)g_{i},g_{j}\rangle_{\pi}$
for $\alpha\in\mathbb{R}$ where $E$ is the spectral resolution in (11).
Finally, when $i=j$, the measure $F_{g_{i}}=F_{g_{i},g_{i}}$ is positive.

The proof is deferred to Supplementary Material S1.1.We note in general
$\gamma_{g_{i},g_{j}}(k)$ is not a moment sequence, as $F_{g_{i},g_{j}}$ is
not necessarily a positive measure when $i\neq j$. However, we can represent
$F_{g_{i},g_{j}}$ as a difference of positive measures, and consequently
represent $\gamma_{g_{i},g_{j}}(k)$ as a difference of moment sequences.

Note that for any $f\in\mathcal{H}$, $\Psi_{f,f}:\alpha\to\langle
E(\alpha)f,f\rangle_{\pi}$ is a positive, increasing function of $\alpha$, as
$\Psi_{f,f}(\alpha)=\|E(\alpha)f\|_{L^{2}(\pi)}^{2}$ since each $E(\alpha)$ is
a projection operator and $\alpha\to\|E(\alpha)f\|_{L^{2}(\pi)}$ is an
increasing function of $\alpha$. By the polarization identity and using the
fact that $\langle E(\alpha)g_{i},g_{j}\rangle_{\pi}$ is real due to
Proposition 1, for any $a,b\neq 0$,

$\displaystyle\Psi_{g_{i},g_{j}}(\alpha)$ $\displaystyle=\langle
E(\alpha)g_{i},E(\alpha)g_{j}\rangle_{\pi}$
$\displaystyle=\frac{1}{4ab}\\{\|E(\alpha)(ag_{i}+bg_{j})\|_{\pi}^{2}-\|E(\alpha)(ag_{i}-bg_{j})\|_{\pi}^{2}\\}$
$\displaystyle=\frac{1}{4ab}\\{\Psi_{ag_{i}+bg_{j},ag_{i}+bg_{j}}(\alpha)-\Psi_{ag_{i}-bg_{j},ag_{i}-bg_{j}}(\alpha)\\},$
(13)

and therefore, the measure $F_{g_{i},g_{j}}$ can be written as the difference
of two positive measures
$F_{g_{i},g_{j}}=\frac{1}{4ab}\\{F_{ag_{i}+bg_{j}}-F_{ag_{i}-bg_{j}}\\}$ where
$F_{ag_{i}+bg_{j}}$ and $F_{ag_{i}-bg_{j}}$ are positive measures uniquely
defined by positive, increasing functions $\Psi_{ag_{i}+bg_{j},ag_{i}+bg_{j}}$
and $\Psi_{ag_{i}-bg_{j},ag_{i}-bg_{j}}$ respectively. Consequently,
$\gamma_{g_{i},g_{j}}=\\{\gamma_{g_{i},g_{j}}(k)\\}_{k\in\mathbb{Z}}$ can also
be represented as the difference of two moment sequences, since

$\displaystyle\gamma_{g_{i},g_{j}}(k)$
$\displaystyle=\frac{1}{4ab}\int\alpha^{|k|}F_{ag_{i}+bg_{j}}(d\alpha)-\frac{1}{4ab}\int\alpha^{|k|}F_{ag_{i}-bg_{j}}(d\alpha)$
(14)

and both
$\gamma_{ag_{i}+bg_{j}}(k)=\int\alpha^{|k|}F_{ag_{i}+bg_{j}}(d\alpha)$ and
$\gamma_{ag_{i}+bg_{j}}(k)=\int\alpha^{|k|}F_{ag_{i}-bg_{j}}(d\alpha)$ are
moment sequences belonging to
$\gamma_{ag_{i}+bg_{j}}\in\mathscr{M}_{\infty}({\Delta(F_{ag_{i}+bg_{j}})})\cap\ell_{2}(\mathbb{Z})$
and
$\gamma_{ag_{i}-bg_{j}}\in\mathscr{M}_{\infty}({\Delta(F_{ag_{i}-bg_{j}})})\cap\ell_{2}(\mathbb{Z})$.

In light of this observation, we can estimate the $\gamma_{g_{i},g_{j}}$
sequence by employing the momentLS sequence estimators for the autocovariance
sequences of the two univariate processes
$\\{ag_{i}(X_{t})+bg_{j}(X_{t})\\}_{t\geq 0}$ and
$\\{ag_{i}(X_{t})-bg_{j}(X_{t})\\}_{t\geq 0}$. Note that the choice of $a$ and
$b$ determines the weighting applied to $g_{i}$ and $g_{j}$ when constructing
these two univariate processes. By Theorem 1, as long as we use $\delta$ such
that
$0<\delta\leq\min\\{\Delta(F_{ag_{i}+bg_{j}}),\Delta(F_{ag_{i}-bg_{j}})\\}$,
both $\Pi_{\delta}(r_{{ag_{i}+bg_{j}},M})$ and
$\Pi_{\delta}(r_{{ag_{i}-bg_{j}},M})$ will be strongly consistent for
$\gamma_{ag_{i}+bg_{j}}$ and $\gamma_{ag_{i}-bg_{j}}$, and therefore
$\hat{\gamma}_{M,ij}=(4ab)^{-1}\\{\Pi_{\delta}(r_{{ag_{i}+bg_{j}},M})-\Pi_{\delta}(r_{{ag_{i}-bg_{j}},M})\\}$
will be strongly consistent for $\gamma_{g_{i},g_{j}}$. Furthermore, the
following Proposition shows that
$\min\\{\Delta(F_{ag_{i}+bg_{j}}),\Delta(F_{ag_{i}-bg_{j}})\\}$ can be further
lower-bounded by $\min\\{\Delta(F_{g_{i}}),\Delta(F_{g_{j}})\\}$, therefore
the permissible range of values for $\delta$ for estimating
$\gamma_{ag_{i}+bg_{j}}$ and $\gamma_{ag_{i}-bg_{j}}$ can be given based on
$\Delta(F_{g_{i}})$ and $\Delta(F_{g_{j}})$ for any values of $a$ and $b$.

######  Proposition 2.

Let $g:\mathsf{X}\to\mathbb{R}^{d}$ be an $\mathbb{R}^{d}$-valued function
such that each $g_{j}\in L^{2}(\pi)$, $j=1,\dots,d$. For any
$\mathbf{c}\in\mathbb{R}^{d}$, ${\rm
Supp}(F_{\mathbf{c}^{\top}g})\subseteq\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})$

The proof is deferred to Supplementary Material S1.2. Note that if
$\delta_{i}$ and $\delta_{j}$ are valid in the sense that
$0<\delta_{i}\leq\Delta(F_{g_{i}})$ and $0<\delta_{j}\leq\Delta(F_{g_{j}})$,
then $\delta_{ij}=\min\\{\delta_{i},\delta_{j}\\}$ is also a valid input for
estimating $\gamma_{ag_{i}+bg_{j}}$ and $\gamma_{ag_{i}-bg_{j}}$. This is
because by Proposition 2,

$\displaystyle\operatorname{Supp}(F_{ag_{i}+bg_{j}})$
$\displaystyle\subseteq\operatorname{Supp}(F_{g_{i}})\cup\operatorname{Supp}(F_{g_{j}})$
$\displaystyle\subseteq[-1+\min\\{\Delta(F_{g_{i}}),\Delta(F_{g_{j}})\\},1-\min\\{\Delta(F_{g_{i}}),\Delta(F_{g_{j}})\\}]$
$\displaystyle\subseteq[-1+\delta_{ij},1-\delta_{ij}].$

and similarly
$\operatorname{Supp}(F_{ag_{i}-bg_{j}})\subseteq[-1+\delta_{ij},1-\delta_{ij}]$
as well.

We now propose the following momentLS sequence estimator for
$\gamma_{g_{i},g_{j}}$ on $\mathbb{Z}$. For a given vector
$\boldsymbol{\delta}=[\delta_{1},\dots,\delta_{d}]$ and a weight vector
$\mathbf{w}_{M}=[w_{M1},\dots,w_{Md}]$, we define the momentLS estimator for
$\gamma_{g_{i},g_{j}}$, $i\leq j$, as

$\displaystyle\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot;\mathbf{w}_{M})=\begin{cases}\Pi_{\delta_{i}}(r_{M,g_{i}})(\cdot)&i=j\\\
\frac{1}{4w_{Mi}w_{Mj}}\\{\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})(\cdot)-\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}-w_{Mj}g_{j}})(\cdot)\\}&i<j\end{cases}$

where $\delta_{ij}=\min\\{\delta_{i},\delta_{j}\\}$. When $i>j$, we simply let
$\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot;\mathbf{w}_{M})=\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ji}(\cdot;\mathbf{w}_{M})$.

The weight vector $\mathbf{w}_{M}$ determines relative weights for each
$g_{i}$. We allow $\mathbf{w}_{M}$ to be random as long as each element
$w_{Mi}$ converges to a non-zero number $w_{i}\neq 0$ almost surely. In the
next section, we will show that the momentLS estimator with random weights
$\mathbf{w}_{M}$ will converge to the same quantity as the momentLS estimator
with non-random weights $\mathbf{w}$ if $\mathbf{w}_{M}$ converges to
$\mathbf{w}$ elementwise almost surely. In terms of the actual choice of
weights, we propose to scale each $g_{i}$ based on their size, specifically
using the norm of each centered $g_{i}$, i.e.,
$\|\check{g}_{i}\|_{L^{2}(\pi)}=\sqrt{\int\check{g}_{i}(x)^{2}\pi(dx)}=\sqrt{\gamma_{g_{i}}(0)}$.
The lag-0 autocovariance $\gamma_{g_{i}}(0)$ is the variance of $g_{i}(X)$ for
$X\sim\pi$, and thus $\sqrt{\gamma_{i}(0)}>0$ provided $g_{i}$ is not
$\pi$-a.e. constant. We propose to use $\mathbf{w}_{M}=\mathbf{\tilde{w}}_{M}$
where the $i$th entry $\tilde{w}_{Mi}$ of $\mathbf{\tilde{w}}_{M}$ is the
inverse of the estimate of $\sqrt{\gamma_{g_{i}}(0)}$ based on the lag-0
empirical autocovariance, i.e.,

$\displaystyle\tilde{w}_{Mi}=1/\sqrt{r_{{g_{i}},M}(0)},\,\,\,i=1,\dots,d.$
(15)

When we use the choice $\mathbf{w}_{M}=\mathbf{\tilde{w}}_{M}$, we will drop
the dependence on $\mathbf{w}_{M}$ in
$\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot;\mathbf{w}_{M})$ and define
$\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot)=\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot;\mathbf{\tilde{w}}_{M})$.
Note for $\tilde{w}_{Mi}$ in (15), we have $\tilde{w}_{Mi}\to
1/\sqrt{\gamma_{g_{i}}(0)}$ almost surely as $M\to\infty$ for $i=1,\dots,d$.

We also define the momentLS estimator for the asymptotic variance of $\Sigma$
in the multivariate Markov chain CLT. We first propose a plug-in estimator
that essentially estimates each $\Sigma_{ij}$ element-wise. We define
$\hat{\Sigma}_{pw}^{(\boldsymbol{\delta})}\in\mathbb{R}^{d\times d}$ by

$\displaystyle\hat{\Sigma}_{pw,ij}^{(\boldsymbol{\delta})}=\sigma^{2}(\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij})=\sum_{k\in\mathbb{Z}}\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(k)=\sum_{k\in\mathbb{Z}}\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(k;\mathbf{\tilde{w}}_{M}),$
(16)

based on the momentLS estimates of $\gamma_{g_{i},g_{j}}$, for each $1\leq
i,j\leq d$, where we recall the definition
$\sigma^{2}(m)=\sum_{k\in\mathbb{Z}}m(k)$ for sequences $m$ on $\mathbb{Z}$.

One drawback of the proposed estimator is that the resulting asymptotic
variance estimator $\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$ is not
necessarily positive semi-definite. In such cases, we propose a principled
refinement of the first estimator which results in a valid covariance matrix.
The key idea is to re-estimate the eigenvalues of the asymptotic variance
matrix using momentLS estimators based on the eigenvectors from the first
estimator $\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$. To illustrate the idea,
suppose $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{d}$ are eigenvalues
of $\Sigma$ and let $\mathbf{U}=[\mathbf{U}_{1},\dots,\mathbf{U}_{d}]$ where
$U_{1},...,U_{d}$ denote the corresponding orthogonal eigenvectors. For
expositional simplicity, suppose for now that all eigenvalues are distinct. If
we choose $\mathbf{c}=\mathbf{U}_{j}$, the asymptotic variance
$\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})$ will be $\lambda_{j}$, since
$\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})=\sum_{k\in\mathbb{Z}}\mathbf{c}^{\top}\operatorname{Cov}_{\pi}(g(X_{0}),g(X_{k}))\mathbf{c}=\mathbf{c}^{\top}\Sigma\mathbf{c}=\lambda_{j}$.
Then, for each $j=1,\dots,d$, if we have a consistent (up to sign) estimator
$\hat{\mathbf{U}}_{Mj}$ for $\mathbf{U}_{j}$, we expect that we can
consistently estimate the $j$th eigenvalue $\lambda_{j}$ using the momentLS
estimator $\Pi_{\delta}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M})$, provided we
choose a sufficiently small $\delta>0$. Moreover, considering Proposition 2,
$\delta=\min\\{\delta_{1},\dots,\delta_{d}\\}$ will be a valid choice of
$\delta$ if each $\delta_{i}$ is chosen to be valid for estimating
$\gamma_{g_{i}}$ for $i=1,\dots,d$.

Now we introduce a positive semi-definite refinement estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$ from the plug-in estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$ and the momentLS asymptotic
variance estimator $\hat{\Sigma}^{(\boldsymbol{\delta})}$ as a combination of
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$ and
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$. Since
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$ is a symmetric matrix by
construction, there exist real eigenvalues
$\hat{\lambda}_{1}^{pw}\geq\dots\geq\hat{\lambda}_{d}^{pw}$ and eigenvectors
$[\hat{\mathbf{U}}_{M1},\dots,\hat{\mathbf{U}}_{Md}]$ such that
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}\hat{\mathbf{U}}_{Mj}=\hat{\lambda}_{j}^{pw}\hat{\mathbf{U}}_{Mj}$,
for $j=1,\dots,d$. We then re-estimate each $\lambda_{j}$ by
$\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M}))$,
the momentLS asymptotic variance estimator for
$\\{\hat{\mathbf{U}}_{M}^{\top}g(X_{t})\\}_{t\geq 0}$, where
$\underline{\delta}=\min_{1\leq i\leq d}\delta_{i}$. We define
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$ by

$\displaystyle\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}=\hat{\mathbf{U}}_{M}\hat{\mathbf{\Lambda}}_{\boldsymbol{\delta}}\hat{\mathbf{U}}_{M}^{\top}$
(17)

where
$\hat{\mathbf{U}}_{M}=[\hat{\mathbf{U}}_{M1},\dots,\hat{\mathbf{U}}_{Md}]\in\mathbb{R}^{d\times
d}$ and $\hat{\mathbf{\Lambda}}_{\boldsymbol{\delta}}$ is a $d$ by $d$
diagonal matrix with $j$th diagonal entry
$\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M}))$.
Finally, we define our momentLS estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}$ for the asymptotic variance $\Sigma$ as

$\hat{\Sigma}^{(\boldsymbol{\delta})}=\begin{cases}\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}&\lambda_{\rm
min}(\hat{\Sigma}_{pw}^{(\boldsymbol{\delta})})\geq 0\\\
\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}&\lambda_{\rm
min}(\hat{\Sigma}_{pw}^{(\boldsymbol{\delta})})<0.\end{cases}$ (18)

In the next section, we will show that both the plug-in estimator and the
refined positive semi-definite matrix estimator are strongly consistent for
$\Sigma$, provided that the input vector $\boldsymbol{\delta}$ is valid in the
sense that $0<\delta_{i}\leq\Delta(F_{g_{i}})$ for all $i=1,\dots,d$. This
ensures the strong consistency of the moment LS estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}$ for $\Sigma$.

##  4 Statistical guarantees

First, we have the following Lemma which shows that the momentLS sequence and
asymptotic variance estimators based on random weights $\mathbf{c}_{M}$ are
still strongly consistent for $\gamma_{\mathbf{c}^{\top}g}$, if
$\mathbf{c}=\lim_{M\to\infty}\mathbf{c}_{M}$ $P_{x}$-almost surely for any
initial condition $x\in\mathsf{X}$.

######  Lemma 1.

Let $\mathbf{c}_{M}=[c_{M1},\dots,c_{Md}]$ denote a sequence of random vectors
satisfying $\lim_{M\to\infty}\mathbf{c}_{M}=\mathbf{c}$, $P_{x}$-almost surely
for any initial condition $x\in\mathsf{X}$. Suppose the Markov chain $X$
satisfies (A.1) \- (A.3). Let $g:\mathsf{X}\to\mathbb{R}^{d}$ be a function
satisfying (B.1). Suppose $\delta$ is chosen so that
$0<\delta\leq\Delta(F_{\mathbf{c}^{\top}g})$. We have 1\.
$\|\Pi_{\delta}(r_{{\mathbf{c}_{M}^{\top}g},M})-\gamma_{\mathbf{c}^{\top}g}\|_{2}\to
0$, and 2\.
$|\sigma^{2}(\Pi_{\delta}(r_{{\mathbf{c}_{M}^{\top}g},M}))-\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})|\to
0$, $P_{x}$-almost surely, for any initial condition $x\in\mathsf{X}$.

The proof of Lemma 1 is in Supplementary Material S2.1. Now we present the
following Theorem which guarantees the almost sure convergence of the plug-in
sequence and asymptotic variance estimators.

######  Theorem 2.

Suppose the Markov chain $X$ satisfies (A.1) \- (A.3). Let
$g:\mathsf{X}\to\mathbb{R}^{d}$ be a function satisfying (B.1). Let
$\mathbf{w}_{M}\in\mathbb{R}^{d}$ be a random vector such that
$\mathbf{w}_{M}\to\mathbf{w}=[w_{1},\dots,w_{d}]$ $P_{x}$-almost surely for
any $x\in\mathsf{X}$ and $w_{i}\neq 0$ for all $1\leq i\leq d$. Suppose
$\boldsymbol{\delta}$ is given so that $0<\delta_{i}\leq\Delta(F_{g_{i}})$.
For any $1\leq i,j\leq d$, we have

  1. 1.

$\sum_{k\in\mathbb{Z}}\\{\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(k;\mathbf{w}_{M})-\gamma_{g_{i},g_{j}}(k)\\}^{2}\to
0$

  2. 2.

$|\sum_{k\in\mathbb{Z}}\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(k;\mathbf{w}_{M})-\sum_{k\in\mathbb{Z}}\gamma_{g_{i},g_{j}}(k)|\to
0$

$P_{x}$-almost surely, as $M\to\infty$, for any initial condition
$x\in\mathsf{X}$.

The proof of Theorem 2 is in Supplementary Material S2.2. The second statement
of Theorem 2 implies as a Corollary the strong consistency of
$\hat{\Sigma}_{pw}^{\boldsymbol{\delta}}$, based on the weights
$\tilde{\mathbf{w}}_{M}$ in (15), for $\Sigma$:

######  Corollary 1.

Suppose the Markov chain $X$ satisfies (A.1) \- (A.3). Let
$g:\mathsf{X}\to\mathbb{R}^{d}$ be a function satisfying (B.1). Assume that
${\rm Var}_{\pi}(g_{j}(X_{0}))>0$ for all $j=1,\dots,d$. Suppose
$\boldsymbol{\delta}$ is given so that $0<\delta_{i}\leq\Delta(F_{g_{i}})$ for
$i=1,\dots,d$. Then the plug-in estimator
$\hat{\Sigma}_{pw}^{(\boldsymbol{\delta})}\in\mathbb{R}^{d\times d}$ is
$P_{x}$-strongly consistent for $\Sigma\in\mathbb{R}^{d\times d}$ for any
$x\in\mathsf{X}$.

In Corollary 1, the condition ${\rm Var}_{\pi}(g_{j}(X_{0}))>0$ implies that
the function $g_{j}$ is not $\pi$-a.e. constant, and is necessary in order to
ensure $\tilde{w}_{Mi}$ in (15) converge to the finite value
$1/\text{Var}_{\pi}(g_{i}(X_{0}))$ for $i=1,...,d$.

We now move onto the convergence results for the refined positive semi-
definite estimator $\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$. We need the
following Lemma, which ensures that the convergence of the momentLS estimator
for $\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})$ can be controlled uniformly over
$\mathbf{c}$. In particular, recall that
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$ estimator relies on estimating
the asymptotic variance based on $\hat{\mathbf{U}}_{Mj}^{\top}g(X_{t})$ for
$j=1,\dots,d$. Lemma 2 allows us to control the distance
$|\sigma^{2}(\Pi_{\delta}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M}))-\sigma^{2}(\gamma_{\hat{\mathbf{U}}_{Mj}^{\top}g})|$
based on the convergence in the worst-case direction.

######  Lemma 2.

Suppose the Markov chain $X$ satisfies (A.1) \- (A.3). Let
$g:\mathsf{X}\to\mathbb{R}^{d}$ be a function which satisfies (B.1). Suppose
$\delta$ is chosen so that $0<\delta\leq\min_{1\leq i\leq
d}\Delta(F_{g_{i}})$. We have
$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}|\sigma^{2}(\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M}))-\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})|\to
0$ $P_{x}$-almost surely, for any $x\in\mathsf{X}$.

The proof is deferred to Supplementary Material S2.3. The proof follows
similar lines as in the proof of the strong convergence of
$\sigma^{2}(\Pi_{\delta}(r_{\mathbf{c}^{\top}g}))$ for fixed $\mathbf{c}$, but
requires some additional arguments to control the cross-products terms
$\gamma_{g_{i},g_{j}}(\cdot)$ by diagonal terms $\gamma_{g_{i}}(\cdot)$ and
$\gamma_{g_{j}}(\cdot)$ and obtain a bound which is independent of
$\mathbf{c}$ using Holder’s inequality. Finally, we present the following
Theorem which shows that the refined estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$ is also strongly consistent for
the true asymptotic variance $\Sigma$. The proof is deferred to Supplementary
Material S2.4.

######  Theorem 3.

Assume the same conditions as in Corollary 1. The refined estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$ converges $P_{x}$-almost surely
to the true asymptotic variance $\Sigma$ for any initial condition
$x\in\mathsf{X}$.

Finally, we summarize the strong consistency results for the momentLS
auto(cross)-covariance sequence estimators
$\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}$ and asymptotic variance
estimator $\hat{\Sigma}^{(\boldsymbol{\delta})}$. This Corollary is a direct
consequence of Theorem 2, Corollary 1, and Theorem 3.

######  Corollary 2.

Assume the same conditions as in Corollary 1. We have

  1. 1.

$\max_{1\leq i,j\leq
d}\|\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}-\gamma_{g_{i},g_{j}}\|_{2}\to
0$, $P_{x}$-almost surely

  2. 2.

$\hat{\Sigma}^{(\boldsymbol{\delta})}$ is a positive semi-definite matrix, and
$\|\hat{\Sigma}^{(\boldsymbol{\delta})}-\Sigma\|_{F}\to 0$, $P_{x}$-almost
surely

for each initial condition $x\in\mathsf{X}$.

##  5 Empirical studies

In this section, we conduct an empirical comparison between the proposed
momentLS estimator and other state-of-the-art methods for estimating the
asymptotic variance $\Sigma$ in the multivariate CLT. We assess the relative
error $\|\Sigma^{-1/2}(\hat{\Sigma}-\Sigma)\Sigma^{-1/2}\|_{F}$, which
quantifies the relative distance between the true asymptotic variance matrix
$\Sigma$ and the estimated variance matrix $\hat{\Sigma}$ produced by various
methods. Additionally, we evaluate the empirical coverage probabilities of the
confidence regions $C_{\alpha}(X;\hat{\Sigma})$ based on different estimates
of $\hat{\Sigma}$. Here, for a chain $X=(X_{0},\dots,X_{M-1})$,
$C_{\alpha}(X;\hat{\Sigma})$ is defined as follows:
$C_{\alpha}(X;\hat{\Sigma})=\left\\{\mu\in\mathbb{R}^{d}:M\left(\hat{\mu}_{gM}-\mu\right)^{\top}\hat{\Sigma}^{-1}\left(\hat{\mu}_{gM}-\mu\right)<\chi^{2}_{1-\alpha,d}\right\\},$
where $\hat{\mu}_{gM}=M^{-1}\sum_{t=0}^{M-1}g(X_{t})$ is the estimated mean
computed from the chain $X$.

We evaluate the performance of the different methods through simulated and
real-data examples. For the simulated examples, we consider two setups:
multivariate functions of a Markov chain generated using a finite state
Metropolis-Hastings sampler, and a Vector AR1 process. In both cases, we can
analytically compute the true asymptotic variances $\Sigma$. For real-data
examples, we investigate a Bayesian logistic regression model using the BUPA’s
liver disorders dataset from the UCI machine learning repository. We
considered two popular samplers: the random walk Metropolis sampler and the No
U-Turn Sampler (hoffman2014no) from STAN (stan2019) which is a dynamic variant
of Hamiltonian Monte Carlo.

###  5.1 Settings for simulated chains

##### Metropolis-Hastings chain

We consider a Metropolis-Hastings chain on a finite discrete state space
$(\mathsf{X},2^{\mathsf{X}})$ where $\mathsf{X}=\\{1,2,...,s\\}$. We randomly
generated the target distribution $\pi:2^{\mathsf{X}}\to[0,1]$ and the
proposal distribution $P:\mathsf{X}\times 2^{\mathsf{X}}\to[0,1]$, where we
drew a length $s$ random vector $U=[U_{1},U_{2},...,U_{s}]^{\top}$ and let
$\pi(\\{i\\})=U_{i}/(\sum_{i^{\prime}=1}^{s}U_{i^{\prime}})$, and for each
$i\in\mathsf{X}$ we drew a length $s$ random vector
$V_{i}=[V_{i1},\dots,V_{is}]^{\top}$ and let
$P(i,\\{j\\})=V_{ij}/\sum_{j=1}^{s}V_{ij}$. Then we constructed the
Metropolis-Hastings transition kernel $Q$ based on $P$ and $\pi$. Finally, we
simulated $g:\mathsf{X}\to\mathbb{R}^{d}$ by $g(i)\sim N_{d}(0,\Sigma_{g})$
i.i.d. for $i\in\\{1,\dots,s\\}$ where
$\textrm{diag}(\Sigma_{g})=[1,\dots,d^{2}]$ and
$\Sigma_{g,ij}=\sqrt{\Sigma_{g,ii}\Sigma_{g,jj}}\rho^{|i-j|}$. With this
choice of $g$, $g_{j}(X_{t}),j=1,\dots,d$ will be correlated with each other,
and have different magnitudes.

With a slight abuse of notation, let $Q\in\mathbb{R}^{s\times s}$ denote the
transition matrix defined as $Q_{ij}=Q(i,\\{j\\})$. Since $s=100$ is
relatively small, we can compute the eigenvalues and eigenvectors for $Q$. Let
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{s}$ denote the eigenvalues
of $Q$ and $\phi_{1},\dots\phi_{s}$ be the corresponding eigenvectors. Suppose
the eigenvectors $\phi_{i}$ are normalized so that
$\braket{\phi_{i},\phi_{j}}_{\pi}=1[i=j]$. Note we have $\lambda_{1}=1$ and
$\phi_{1}=\mathbf{1}_{s}$ since $Q\mathbf{1}_{s}=\mathbf{1}_{s}$. For each
$i=1,\dots,s$, we can write $g_{i}$ and
$\check{g}_{i}=g_{i}-E_{\pi}[g_{i}(X_{0})]\mathbf{1}_{s}$ as
$g_{i}(k)=\sum_{l=1}^{s}\braket{g_{i},\phi_{l}}_{\pi}\phi_{l}(k)$ and
$\check{g}_{i}(k)=\sum_{l=2}^{s}\braket{g_{i},\phi_{l}}_{\pi}\phi_{l}(k)$
since $\braket{g_{i},\phi_{1}}_{\pi}=E_{\pi}[g_{i}(X_{0})]$. Then,
$\textrm{Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}))=\braket{\check{g}_{i},Q{k}{}\check{g}_{j}}{}{\pi}=\sum_{l=2}^{s}\braket{g_{i},\phi_{l}}_{\pi}\braket{g_{j},\phi_{l}}_{\pi}\lambda_{l}^{|k|}$
since $\langle\phi_{l},\phi_{m}\rangle_{\pi}=1[l=m]$. In particular, we have
$F_{g_{i},g_{j}}=\sum_{l=2}^{s}\braket{g_{i},\phi_{l}}_{\pi}\braket{g_{j},\phi_{l}}_{\pi}\delta_{\\{\lambda_{l}\\}}$
where $\delta_{\\{a\\}}$ denotes a unit point mass at $a$. Note
$F_{g_{i},g_{j}}$ is in general a signed measure on
$(\mathbb{R},\mathscr{B}_{\mathbb{R}})$, since the weight
$\braket{g_{i},\phi_{l}}_{\pi}\braket{g_{j},\phi_{l}}_{\pi}$ at each
$\lambda_{l}$ is real, but not necessarily positive, except the case when
$i=j$, where we have
$\braket{g_{i},\phi_{l}}_{\pi}\braket{g_{i},\phi_{l}}_{\pi}=\langle
g_{i},\phi_{l}\rangle_{\pi}^{2}\geq 0$. Finally, we have the asymptotic
variance $\Sigma$ where
$\Sigma_{ij}=\sum_{l=2}^{s}\braket{g_{i},\phi_{l}}_{\pi}\braket{g_{j},\phi_{l}}_{\pi}\frac{1+\lambda_{l}}{1-\lambda_{l}}.$

##### Vector autoregressive process

As a second example, we consider a vector autoregressive process of order $1$
(VAR(1)). Here we have a stationary Markov chain on a continuous state space
$(\mathsf{X},\mathscr{X})=(\mathbb{R}^{d},\mathscr{B}_{\mathbb{R}^{d}})$ where
for $t=0,1,2,\ldots$, $X_{t}\in\mathbb{R}^{d}$, the process at time $t$ is
defined by $X_{t}=AX_{t-1}+\epsilon_{t}$ where $A$ is a $d\times d$ matrix,
$\epsilon_{t}\stackrel{{\scriptstyle iid}}{{\sim}}N_{d}(0,\Sigma_{\epsilon})$,
and $\Sigma_{\epsilon}$ is a $d\times d$ positive definite matrix.
osawa1988reversibility has showed that the VAR(1) process defined above is
reversible if and only if $A\Sigma_{\epsilon}$ is symmetric. Moreover, if
$\lambda_{\rm max}(A)<1$, the stationary distribution $\pi$ for the $X_{t}$
chain is $N_{d}(0_{d},V)$ where $V=(I-A^{2})^{-1}\Sigma_{\epsilon}$. We let
$A=D_{\rho}+(0.1)^{2}(1_{d}1_{d}^{\top}-I_{d})$ where $D_{\rho}={\rm
diag}(0.9,0.9,-0.9,-0.9)$, and $\Sigma_{\epsilon}=I_{d}$. The simulated
process is reversible because $A$ is symmetric. Moreover, since
$\operatorname{Cov}_{\pi}(X_{k},X_{0})=A\operatorname{Cov}_{\pi}(X_{k-1},X_{0})=\dots=A^{k}V$,
the asymptotic variance $\Sigma$ is
$\Sigma=\sum_{k\in\mathbb{Z}}A^{k}V=2\sum_{k=0}^{\infty}A^{k}-V=2(I-A)^{-1}V-V.$

####  5.1.1 Descriptions of estimators

We compared the asymptotic variance estimation performance of the following
estimators: 1. (SV(Bartlett)) the spectral variance estimator with the
modified Bartlett window $w_{M}(k)=(1-|k|/b_{M})1\\{|k|<b_{M}\\}$
(Damerdji1991-oj; flegal2010batch; vats2018strong), 2. (BM, OBM) the batch
means and overlapping batch means estimator (vats2019multivariate), 3. (mtv-
Init) the multivariate initial sequence estimator from dai2017multivariate,
and 4. (mtv-MLSE) our multivariate moment least squares estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}$ defined in (18). All empirical
experiments that we present in this section were performed using R software
(rlang). We used the momentLS R package available on GitHub
111https://github.com/hsong1/momentLS for computing univariate momentLS
estimators. All other estimators were computed using the mcmcse R package
(mcmcse_R).

Hyperparameters are required for the spectral variance estimator, BM, OBM, and
momentLS estimators. A truncation point for SV(Bartlett) or a batch size for
BM and OBM needs to be specified in advance. The vector
$\boldsymbol{\delta}=[\delta_{1},\dots,\delta_{d}]$ needs to be specified for
momentLS estimators, which determines the moment spaces onto which each
empirical autocovariance sequence $r_{gM,j}$, $j=1,\dots,d$ is projected. For
SV(Bartlett), BM, and OBM, we used the truncation point (batch size) tuning
method implemented in the R package mcmcse (liu2021batch). For the
multivariate momentLS estimator, we utilized the tuning method proposed in
berg2022efficient to select $\delta_{j}$ for each univariate sequence
$\\{g_{j}(X_{t})\\}_{t\geq 0}$, defined as
$\tilde{\delta}_{Mj}=0.8\frac{1}{L}\sum_{l=1}^{L}\hat{\delta}_{M/L,j}^{(l)}.$
Here, $\hat{\delta}_{M/L,j}^{(l)}$ is an estimator for $\Delta(F_{g_{j}})$
computed from the $l$th batch of empirical autocovariances. We used $L=5$.
Subsequently we let
$\boldsymbol{\delta}=[\tilde{\delta}_{M1},\dots,\tilde{\delta}_{Md}]$.

###  5.2 Results for simulated chains

We compared the relative error and coverage rate from various methods
(SV(Bartlett), BM, OBM, Mtv-Init, and Mtv-mLSE) with varying chain lengths
$M\in\\{5000,10000,20000,40000,80000\\}$. In particular, at each sample size
$M$, $B=1000$ independent chains were generated, and the various estimators
were computed on each of these replicate chains. For each method and at each
sample size $M$, we estimated the relative error and coverage rate as
$\widehat{\textrm{relative
error}}_{{}_{M}}=B^{-1}\sum_{b=0}^{B}\|\Sigma^{-1/2}(\hat{\Sigma}^{(b)}_{M}-\Sigma)\Sigma^{-1/2}\|_{F}$
and $\widehat{\textrm{coverage
probability}}_{{}_{M}}=B^{-1}\sum_{b=0}^{B}1\\{\mu_{g}\in
C_{\alpha}(X_{M}^{(b)};\hat{\Sigma}^{(b)}_{M})\\}$ where
$X_{M}^{(b)}=(X_{0}^{(b)},\dots,X_{M-1}^{(b)})$ is the $b$th simulated chain
of length $M$ and $\hat{\Sigma}^{(b)}_{M}$ is the estimated asymptotic
variance using $X_{M}^{(b)}$. We used $\alpha=0.05$ in
$C_{\alpha}(X^{(b)};\hat{\Sigma}^{(b)})$.

(a) Metropolis-Hastings

(b) Vector AR1

Figure 1:  Estimated relative errors
$\|\Sigma^{-1/2}(\hat{\Sigma}-\Sigma)\Sigma^{-1/2}\|_{F}$ and coverage
probabilities $P_{\pi}(\mu_{gM}\in C_{\alpha}(X;\hat{\Sigma}))$ with
$\alpha=0.05$ are computed using different methods for the Metropolis-Hastings
and Vector AR1 example. Error bars indicate one standard error. For coverage
probability plots, the dotted horizontal lines indicate the $95$% nominal
coverage level.

Figure 1 displays the estimated relative errors and coverage probabilities
from the SV(Bartlett), BM, OBM, mtv-Init, and mtv-mLSE methods. The
multivariate momentLS estimator (mtv-mLSE) appears to perform the best in
terms of both relative error and coverage probability in both examples,
demonstrating the smallest relative errors and empirical coverage rates
closest to the 95% nominal level across all sample sizes. Also, both mtv-mLSE
and mtv-Init estimators seem to exhibit faster convergence rates when compared
to the other methods, which was also the case for their counterparts in the
univariate setting (berg2022efficient). While the multivariate initial
sequence estimator (mtv-Init) performs second best in the Metropolis-Hastings
example, the estimator seems to struggle more in the Vector AR1 case, showing
relatively large relative errors compared to the other methods and tending to
provide under-coverage of the the true parameter. Both the spectral variance
(SV(Bartlett)) and overlapping batch means (OBM) estimators outperform the
batch means (BM) method, with SV and OBM performing very similarly in both
Metropolis-Hastings and VAR1 examples.

###  5.3 Bayesian logistic regression

In this section, we fit a Bayesian logistic regression model to the Liver
Disorders dataset sourced from the UCI machine learning repository. The
analysis focuses on assessing the relationship between blood test panel
predictors (MCV, ALKPHOS, SGPT, SGOT, GAMMAGT) and an individual’s drinking
habits. The dataset consists of $n=341$ observations after removing
duplicates. Each row represents an individual’s blood test results and their
average daily alcohol intake. The goal is to predict an individual’s drinking
status, which is classified as “heavy” if their average daily drink count
exceeds 3, based on the $5$ blood test panel results represented as
$\mathbf{x}_{i}=(x_{i1},\dots,x_{i5})\in\mathbb{R}^{5}$. We denote the binary
response variable as $Y_{i}$, where $Y_{i}=1$ if the average daily drink count
of the $i$th individual exceeds 3 and 0 otherwise. We suppose
$Pr(Y_{i}=1)=\frac{\exp(\alpha+\sum_{j=1}^{5}\beta_{j}x_{ij})}{1+\exp(\alpha+\sum_{j=1}^{5}\beta_{j}x_{ij})}$
and assign a $N(0,5)$ prior on the intercept $\alpha$ and independent $N(0,1)$
priors on the coefficients $\boldsymbol{\beta}=(\beta_{1},\dots,\beta_{5})$.

We sample $\\{\boldsymbol{\beta}(t)\\}_{t=0}^{M-1}$ from the posterior
distribution $\boldsymbol{\beta}|\\{Y_{i}\\}_{i=1}^{n}\sim\pi(\cdot)$ using
two popular samplers: 1. a random walk Metropolis sampler and 2. the No U-Turn
sampler (hoffman2014no). We considered varying chain lengths
$M\in\\{2500,5000,10000,20000,40000\\}$ for each sampler. Similar to previous
examples, we generated $B=1000$ parallel chains from each sampler and
estimated relative error and coverage probability.

For the random walk Metropolis sampler, we used the MCMClogit function in the
MCMCpack R package (JSSv042i09) with a burn-in period of $5000$ iterations.
Regarding the proposal step, the proposal distribution at time $t$ is taken to
be $N(\boldsymbol{\beta}(t),V_{prop})$ where
$V_{prop}=T\left(B_{0}+C^{-1}\right)^{-1}T$, where $T=1.1I_{d}$ is a diagonal
positive definite matrix, $B_{0}=\textrm{diag}(0.2,1,1,1,1,1)$ is the prior
precision matrix, and $C$ is the asymptotic variance-covariance matrix of the
MLEs.

The No U-Turn sampler is a dynamic variant of the Hamiltonian Monte Carlo
(HMC) method, which itself is a Metropolis-Hastings algorithm where proposal
states are generated by integrating a system of Hamiltonian equations using
the leapfrog integrator. HMC methods require the specification of three
algorithm parameters: 1. the discretization time (step size) $\epsilon$, 2.
the mass matrix $M_{mass}$, and 3. the number of leapfrog steps $L$.
hoffman2014no introduced the No U-Turn sampler (NUTS), which automatically
selects the value of $L$ during the sampling process. The No U-Turn sampler
implemented in STAN employs a warmup and sampling phases. During the warmup
phase, the sampler optimizes $\epsilon$ and $M_{mass}$ to match an acceptance
rate target. In the subsequent sampling phase, the sampler generates Markov
chain samples based on the prescribed transition mechanism, which depends on
the HMC parameters $(\epsilon,M_{mass})$. Since we want to generate parallel
chains from the same transition mechanism, we first learned $\epsilon$ and
$M_{mass}$ using $100000$ warm-up samples, and subsequently generated $B$
parallel chains of length $M+5000$ (with $0$ warmup iteration) using the same
$(\epsilon,M_{mass})$. The first $5000$ iterations were treated as iterations
from a burn-in period, and excluded in the final chains. All sampling was
performed using the CmdStan, the command-line interface to the Stan
statistical modeling language (stan2019), as the current R implementation of
Stan does not allow direct specification of $\epsilon$ and $M_{mass}$ during
the sampling phase.

To compare estimated asymptotic variances and coverage probabilities from the
competing methods, we also need accurate reference estimates of the posterior
mean and asymptotic variance for each coefficient. Since both quantities are
unknown, for both the random walk Metropolis sampler and NUTS, we
independently generated $B_{2}=10000$ independent chains
$\\{\boldsymbol{\beta}_{\rm par}^{(b)}(t)\\}_{t=0}^{M_{2}-1}$ of length
$M_{2}=40000$ (after burn-in periods of length $5000$) to estimate true mean
and asymptotic variance. Specifically, let
$\bar{\boldsymbol{\beta}}_{\textrm{par}}^{(b)}$ refer to the sample mean value
of $\boldsymbol{\beta}$ from the $b$th chain for $b=1,\dots,B_{2}$. We use
$\boldsymbol{\beta}_{\textrm{orcl}}=\frac{1}{B_{2}}\sum_{b=1}^{B_{2}}\bar{\boldsymbol{\beta}}_{\textrm{par}}^{(b)}$
to estimate the posterior mean, and
$\Sigma_{\textrm{orcl}}=\frac{M_{2}}{B_{2}-1}\sum_{b=1}^{B_{2}}(\bar{\boldsymbol{\beta}}_{\textrm{par}}^{(b)}-\boldsymbol{\beta}_{\rm
orcl})(\bar{\boldsymbol{\beta}}_{\textrm{par}}^{(b)}-\boldsymbol{\beta}_{\rm
orcl})^{\top}$ to estimate the asymptotic variance matrix.

(a) Random Walk Metropolis sampler

(b) No U-Turn sampler

Figure 2:  Estimated relative errors
$\|\Sigma^{-1/2}(\hat{\Sigma}-\Sigma)\Sigma^{-1/2}\|_{F}$ and coverage
probabilities $P_{\pi}(\mu_{gM}\in C_{\alpha}(X;\hat{\Sigma}))$ with
$\alpha=0.05$ for the Bayesian logistic model example using the random walk
Metropolis sampler and the No U-Turn sampler. Error bars indicate one standard
error. For coverage probability plots, the dotted horizontal lines indicate
the $95$% nominal coverage level.

Figure 2 shows the estimated relative errors and coverage probabilities for
Markov chains from the random walk Metropolis sampler and the No U-Turn
sampler. Similarly to the results in the simulated examples, the multivariate
momentLS estimator (mtv-mLSE) continues to perform the best in terms of both
relative error and coverage probability for both samplers on average. It
consistently achieves the smallest relative errors and empirical coverage
rates closest to the 95% nominal level across all sample sizes except when
$M=2500$. As in the previous case, the SV(Bartlett) and OBM estimators
performed quite similarly in both examples, consistently outperforming the
batch means (BM) method. While the multivariate initial sequence estimator
(mtv-Init) shows the second best performance for the random walk Metropolis
sampler, it appears to show relatively worse relative error performance
compared to SV(Bartlett)and OBM estimator for the no U-Turn sampler for the
considered chain lengths ranging from $2500$ to $40000$.

##  6 Conclusion

In this work, we have introduced novel multivariate momentLS sequence and
asymptotic variance estimators based on univariate momentLS sequence and
asymptotic variance estimators proposed in berg2022efficient. Our approach is
based on the observation that the cross-covariance sequence admits a similar
mixture representation to the autocovariance sequence, with a signed measure
as the mixing measure, and can be expressed as the difference between two
moment sequences. We have shown that the supports of the representing measures
of the cross-covariance sequences can be controlled by the supports of the
representing measures of each autocovariance sequence. Therefore we can use
the hyperparameter $\boldsymbol{\delta}$ based on estimating
$\Delta(F_{g_{j}})$ from each univariate chain $g_{j}(X_{t})$, $j=1,\dots,d$.
Theoretically, we have established strong consistency results for both the
sequence and asymptotic variance MomentLS estimators, provided that each
$\boldsymbol{\delta}_{i}$ is chosen to be sufficiently small. Furthermore,
through simulated and real-data Bayesian logistic regression examples, we have
demonstrated that our estimator empirically outperforms the competing methods.

Supplement to “Multivariate moment least-squares estimators for reversible
Markov chains”
Hyebin Song and Stephen Berg
Department of Statistics, Pennsylvania State University

##  S1 Proofs for Section 3

###  S1.1 Proof for Proposition 1

###### Proof.

Let $\mathcal{H}=L^{2}(\mathsf{X},\mathscr{X},\pi)=L^{2}(\pi)$ be the space of
square integrable functions with respect to $\pi$, with inner product
$\langle\cdot,\cdot\rangle_{\pi}$. Let $\mathscr{B}(\mathcal{H})$ be the
Banach algebra of all bounded linear operators $T$ on $\mathcal{H}$, with norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|T\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}=\sup\\{\|Tf\|_{L^{2}(\pi)}:f\in\mathcal{H},\|f\|_{L^{2}(\pi)}\leq
1\\}$. We recall that $\gamma_{g_{i},g_{j}}(k)={\rm
Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}))=\langle
g_{i},Q_{0}^{k}g_{j}\rangle_{\pi}=\langle Q_{0}^{|k|}g_{i},g_{j}\rangle_{\pi}$
for $k\geq 0$ since $Q_{0}$ is self-adjoint. When $k<0$,
$\gamma_{g_{i},g_{j}}(k)={\rm Cov}_{\pi}(g_{i}(X_{|k|}),g_{j}(X_{0}))=\langle
Q_{0}^{|k|}g_{i},g_{j}\rangle_{\pi}$.

Since $Q_{0}$ is self-adjoint and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Q_{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}\leq
1$, there exists a spectral resolution $\\{E(\lambda)\\}$ of $Q_{0}$ (e.g.,
Theorem 6.1 in stein2009real), which is a function from $\mathbb{R}$ to
$\mathscr{B}(\mathcal{H})$, such that

$\displaystyle Q_{0}=\int_{[a,b]}\lambda dE(\lambda),$ (S-1)

for $a=\min_{f\in\mathcal{H};\|f\|_{L^{2}(\pi)}\leq 1}\langle
Q_{0}f,f\rangle_{\pi}$ and $b=\max_{f\in\mathcal{H};\|f\|_{L^{2}(\pi)}\leq
1}\langle Q_{0}f,f\rangle_{\pi}$ (note that both $a$ and $b$ are real values
due to self-adjointness of $Q_{0}$), which means that for any
$f,g\in\mathcal{H}$,

$\displaystyle\langle Q_{0}f,g\rangle_{\pi}=\int\lambda F_{f,g}(d\lambda)$

where $F_{f,g}$ is the complex-valued measure on $\mathscr{B}_{\mathbb{R}}$
uniquely defined by the function $\Psi_{f,g}:\lambda\to\langle
E(\lambda)f,g\rangle_{\pi}$. We note that $F_{f,g}$ is well-defined since
$\Psi_{f,g}$ is of bounded variation: by the polarization identity, for
$f,g\in\mathcal{H}$,

$\displaystyle\Psi_{f,g}(\lambda)=\langle E(\lambda)f,g\rangle_{\pi}=\langle
E(\lambda)f,E(\lambda)g\rangle_{\pi}=Re(\langle
E(\lambda)f,E(\lambda)g\rangle_{\pi})+iIm(\langle
E(\lambda)f,E(\lambda)g\rangle_{\pi})$

where

$\displaystyle Re(\langle
E(\lambda)f,E(\lambda)g\rangle_{\pi})=\frac{1}{4}\\{\|E(\lambda)(f+g)\|_{L^{2}(\pi)}-\|E(\lambda)(f-g)\|_{L^{2}(\pi)}\\}$
$\displaystyle Im(\langle
E(\lambda)f,E(\lambda)g\rangle_{\pi})=\frac{1}{4}\\{\|E(\lambda)(f+ig)\|_{L^{2}(\pi)}-\|E(\lambda)(f-ig)\|_{L^{2}(\pi)}\\},$
(S-2)

and $\lambda\to\|E(\lambda)f\|_{L^{2}(\pi)}$ is an increasing function of
$\lambda$ for any $f\in\mathcal{H}$ by the properties of spectral resolutions.
Moreover, the complex measure $F_{f,g}$ is finite, i.e.,
$|F_{f,g}|(\mathbb{R})<\infty$ where we use $|\mu|$ denote the total variation
of a complex measure $\mu$. Note that $|F_{f,g}|$ is the same as the positive
measure defined by the total variation function $T_{\Psi_{f,g}}$ of
$\Psi_{f,g}$ (e.g., Theorem 3.29 in folland1999real). By the polarization
inequality (S-2), for any $x\in\mathbb{R}$,

$\displaystyle|F_{f,g}|((-\infty,x])$
$\displaystyle\leq\frac{1}{4}\\{\Psi_{f+g,f+g}(x)+\Psi_{f-g,f-g}(x)+\Psi_{f+ig,f+ig}(x)+\Psi_{f-ig,f-ig}(x)\\}$
$\displaystyle\leq\frac{1}{4}\sum_{j=1}^{4}\|f+a_{j}g\|_{L^{2}(\pi)}^{2}\,\,\textrm{where}\,a_{j}=\pm
1\,\textrm{or}\,\pm i$
$\displaystyle\leq\\{\|f\|_{L^{2}(\pi)}^{2}+\|g\|_{L^{2}(\pi)}^{2}+2\|f\|_{L^{2}(\pi)}\|g\|_{L^{2}(\pi)}\\}<\infty.$

In particular,
$|F_{f,g}|(\mathbb{R})\leq\|f\|_{L^{2}(\pi)}^{2}+\|g\|_{L^{2}(\pi)}^{2}+2\|f\|_{L^{2}(\pi)}\|g\|_{L^{2}(\pi)}<\infty$.

Let $f=g_{i}$ and $g=g_{j}$, which are real-valued functions by assumption.
For a function $\Phi$ continuous on $[a,b]$, the operator $\Phi(Q_{0})$
defined as the unique limit of approximating polynomials has the following
spectral representation

$\displaystyle\langle\Phi(Q_{0})f,g\rangle_{\pi}=\int\Phi(\lambda)F_{f,g}(d\lambda)$
(S-3)

for any $f,g\in\mathcal{H}$ (stein2009real). Let
$\Phi(\lambda)=\lambda^{|k|}$. Then

$\displaystyle\gamma_{g_{i},g_{j}}(k)=\langle
Q_{0}^{|k|}g_{i},g_{j}\rangle_{\pi}=\int\lambda^{|k|}F_{g_{i},g_{j}}(d\lambda)$

where we recall that $F_{g_{i},g_{j}}$ is the complex-valued measure defined
by $\Psi_{g_{i},g_{j}}:\lambda\to\langle E(\lambda)g_{i},g_{j}\rangle_{\pi}$.

We show that $F_{g_{i},g_{j}}$ is a real-valued measure for real-valued
$g_{i},g_{j}\in L_{2}(\pi)$, and non-negative if $g_{i}=g_{j}$. For a real-
valued function $f$, $Q_{0}f(x)=\int Q_{0}(x,dy)f(y)$ is also real-valued.
This implies for the operator $\varphi(Q_{0})$ which is defined based on a
real-valued function $\varphi:[a,b]\to\mathbb{R}$, $\varphi(Q_{0})f$ is also a
real-valued function. In particular, for any $\lambda\in\mathbb{R}$,

$\displaystyle E(\lambda)g_{i}=\begin{cases}0&\lambda<a\\\
\varphi^{\lambda}(Q_{0})g_{j}&a\leq\lambda\leq b\\\
g_{j}&\lambda>b\end{cases}$

is real-valued where $\varphi^{\lambda}(t)=1\\{t\leq\lambda\\}$ (c.f. p310 in
(stein2009real)). Therefore,
$\Psi_{g_{i},g_{j}}(\lambda)=F_{g_{i},g_{j}}((-\infty,\lambda])=Re(F_{g_{i},g_{j}}((-\infty,\lambda]))$,
i.e., $F_{g_{i},g_{j}}$ is a signed measure. When $g_{i}=g_{j}$,
$\Psi_{g_{i},g_{i}}(\lambda)=\langle
E(\lambda)g_{i},g_{i}\rangle_{\pi}=\langle
E(\lambda)g_{i},E(\lambda)g_{i}\rangle_{\pi}=\|E(\lambda)g_{i}\|_{L^{2}(\pi)}^{2}\geq
0$, and $F_{g_{i},g_{i}}$ is a positive measure.

We show that the support of $F_{g_{i},g_{j}}$ is contained in
$[-1+\delta(Q),1-\delta(Q)]$. Here, the support of $F_{g_{i},g_{j}}$ is
defined as ${\rm Supp}(F_{g_{i},g_{j}})={\rm
Supp}(F_{g_{i},g_{j}}^{+})\cup{\rm Supp}(F_{g_{i},g_{j}}^{-})$ where
$F_{g_{i},g_{j}}^{+}$ and $F_{g_{i},g_{j}}^{-}$ denote the positive and
negative parts of $F_{g_{i},g_{j}}$ in the Jordan-Hahn decomposition of
$F_{g_{i},g_{j}}$. First, we note that ${\rm
Supp}(F_{g_{i},g_{j}})\subseteq[a,b]$. This is because the distribution
functions of $F_{g_{i},g_{j}}^{+}$ and $F_{g_{i},g_{j}}^{-}$ are
$\frac{1}{2}(T_{\Psi_{g_{i},g_{j}}}+\Psi_{g_{i},g_{j}})$ and
$\frac{1}{2}(T_{\Psi_{g_{i},g_{j}}}-\Psi_{g_{i},g_{j}})$ respectively, where
$T_{\Psi_{g_{i},g_{j}}}$ is the total variation function of
$\Psi_{g_{i},g_{j}}$, and both $\Psi_{g_{i},g_{j}}$ and
$T_{\Psi_{g_{i},g_{j}}}$ are constant on $(-\infty,a)$ and $(b,\infty)$. Now
we show $[a,b]\subseteq[-1+\delta(Q),1-\delta(Q)]$. Since $Q_{0}$ is self-
adjoint,
$\max\\{|\lambda|;\lambda\in\sigma(Q_{0})\\}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Q_{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}$,
and therefore the spectral gap
$\delta(Q)=1-{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Q_{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{L^{2}(\pi)}=1-\max\\{|a|,|b|\\}$
by Proposition 6.2 in stein2009real. Then we have
$[a,b]\subseteq[-\max\\{|a|,|b|\\},\max\\{|a|,|b|\\}]=[-1+\delta(Q),1-\delta(Q)]$.

Finally, we discuss the uniqueness of $F_{g_{i},g_{j}}$. Suppose there exists
another finite signed measure $\tilde{F}$ with
$\operatorname{Supp}(\tilde{F})\subseteq[-1+\delta(Q),1-\delta(Q)]$ such that
for all $k\in\mathbb{N}$

$\displaystyle\int\lambda^{k}F_{g_{i},g_{j}}(d\lambda)=\int\lambda^{k}\tilde{F}(d\lambda).$

This implies for any real-valued polynomial $p_{N}$,

$\displaystyle\int p_{N}(\lambda)F_{g_{i},g_{j}}(d\lambda)=\int
p_{N}(\lambda)\tilde{F}(d\lambda).$ (S-4)

We show that $F_{g_{i},g_{j}}((-\infty,t])=\tilde{F}((-\infty,t])$ for any
$t\in\mathbb{R}$. This implies that $F_{g_{i},g_{j}}=\tilde{F}$ (e.g., Theorem
3.29 in folland1999real).

Let $a_{0}=-1+\delta(Q)$ and $b_{0}=1-\delta(Q)$ for notational simplicity. We
have that $F_{g_{i},g_{j}}((-\infty,t])=\tilde{F}((-\infty,t])=0$ for
$t<a_{0}$ and $F_{g_{i},g_{j}}((b_{0},\infty))=\tilde{F}((b_{0},\infty))=0$.
Taking $p_{N}(\lambda)=1$ in (S-4), we have
$F_{g_{i},g_{j}}((-\infty,b_{0}])=\tilde{F}((-\infty,b_{0}])$, and therefore,
$F_{g_{i},g_{j}}((-\infty,t])=\tilde{F}((-\infty,t])$ for all $t\geq b_{0}$.
Let $t\in[a_{0},b_{0})$ and $\epsilon>0$ be given. Choose
$\epsilon_{\lambda}>0$ such that both
$|F_{g_{i},g_{j}}|((t,t+\epsilon_{\lambda}])\leq\epsilon/4$ and
$|\tilde{F}|((t,t+\epsilon_{\lambda}])\leq\epsilon/4$ (such
$\epsilon_{\lambda}>0$ exists from the continuity properties of measures).
Define a continuous and bounded function $h_{t}$ such that

$\displaystyle h_{t}(\lambda)=\begin{cases}0&\lambda\leq
a_{0}-\epsilon_{\lambda}\\\
\frac{1}{\epsilon_{\lambda}}(\lambda-(a_{0}-\epsilon_{\lambda}))&a_{0}-\epsilon_{\lambda}<\lambda\leq
a_{0}\\\ 1&a_{0}<\lambda\leq t\\\
-\frac{1}{\epsilon_{\lambda}}(\lambda-(t+\epsilon_{\lambda}))&t<\lambda\leq
t+\epsilon_{\lambda}\\\ 0&\lambda>t+\epsilon_{\lambda}.\end{cases}$

Note that $h_{t}$ and $1_{(-\infty,t])}$ only differ on $(-\infty,a_{0})$ and
$(t,t+\epsilon_{\lambda})$ by definition of $h_{t}$, and therefore

$\displaystyle\int|h_{t}(\lambda)-1_{(-\infty,t]}(\lambda)||F_{g_{i},g_{j}}|(d\lambda)=\int_{(t,t+\epsilon_{\lambda})}|h_{t}(\lambda)-1_{(-\infty,t]}(\lambda)||F_{g_{i},g_{j}}|(d\lambda)\leq|F_{g_{i},g_{j}}|((t,t+\epsilon_{\lambda}])\leq\epsilon/4.$
(S-5)

Similarly,

$\displaystyle\int|h_{t}(\lambda)-1_{(-\infty,t])}||\tilde{F}|(d\lambda)\leq\epsilon/4$
(S-6)

Let $C_{0}=\max\\{|F_{g_{i},g_{j}}|(\mathbb{R}),|\tilde{F}|(\mathbb{R})\\}$.
Note $C_{0}<\infty$ because
$|F_{g_{i},g_{j}}|(\mathbb{R})\leq\|g_{i}\|_{L^{2}(\pi)}^{2}+\|g_{j}\|_{L^{2}(\pi)}^{2}+2\|g_{i}\|_{L^{2}(\pi)}\|g_{j}\|_{L^{2}(\pi)}<\infty$
and $\tilde{F}$ is assumed to be a finite measure. By the Weierstrass
approximation theorem, we can find a polynomial
$p_{N}(\lambda)=\sum_{l=0}^{k}\alpha_{l}\lambda^{k}$ with
$\alpha_{1},\dots,\alpha_{N}\in\mathbb{R}$ on $[a,b]$ such that
$|h_{t}(\lambda)-p_{N}(\lambda)|\leq\epsilon/(2C_{0})$ for all
$\lambda\in[a,b]$. Let
$\Delta_{1}(\lambda)=|h_{t}(\lambda)-1_{(-\infty,t]}(\lambda)|$ and
$\Delta_{2}(\lambda)=|h_{t}(\lambda)-p_{N}(\lambda)|$. We have,

$\displaystyle|F_{g_{i},g_{j}}((-\infty,t])-\tilde{F}((-\infty,t])|$
$\displaystyle=|\int 1_{(-\infty,t]}(\lambda)F_{g_{i},g_{j}}(d\lambda)-\int
1_{(-\infty,t]}(\lambda)\tilde{F}(d\lambda)|$
$\displaystyle\leq|\int\\{1_{(-\infty,t]}(\lambda)-h_{t}(\lambda)\\}F_{g_{i},g_{j}}(d\lambda)-\int\\{1_{(-\infty,t]}(\lambda)-h_{t}(\lambda)\\}\tilde{F}(d\lambda)|$
$\displaystyle\qquad+|\int\\{h_{t}(\lambda)-p_{N}(\lambda)\\}F_{g_{i},g_{j}}(d\lambda)-\int\\{h_{t}(\lambda)-p_{N}(\lambda)\\}\tilde{F}(d\lambda)|$
$\displaystyle\qquad+|\int p_{N}(\lambda)F_{g_{i},g_{j}}(d\lambda)-\int
p_{N}(\lambda)\tilde{F}(d\lambda)|$
$\displaystyle\leq\int\Delta_{1}(\lambda)|F_{g_{i},g_{j}}|(d\lambda)+\int\Delta_{1}(\lambda)|\tilde{F}|(d\lambda)+\int\Delta_{2}(\lambda)|F_{g_{i},g_{j}}|(d\lambda)+\int\Delta_{2}(\lambda)|\tilde{F}|(d\lambda)$

We have
$\int\Delta_{1}(\lambda)|F_{g_{i},g_{j}}|(d\lambda)+\int\Delta_{1}(\lambda)|\tilde{F}|(d\lambda)\leq\epsilon/2$
by (S-5) and (S-6). Moreover, we have,

$\displaystyle\int\Delta_{2}(\lambda)|F_{g_{i},g_{j}}|(d\lambda)+\int\Delta_{2}(\lambda)|\tilde{F}|(d\lambda)\leq\epsilon/(2C_{0})\\{|F_{g_{i},g_{j}}|([a,b])+|\tilde{F}|([a,b])\\}=\epsilon/2$

since $\sup_{\lambda\in[a,b]}\Delta_{2}(\lambda)\leq\epsilon/(2C_{0})$.
Therefore, since $\epsilon$ is arbitrary, we have
$F_{g_{i},g_{j}}((-\infty,t])=\tilde{F}((-\infty,t])$ for all
$t\in\mathbb{R}$. ∎

###  S1.2 Proof for Proposition 2

###### Proof.

By definition $\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}$, and using the
fact that $E(\alpha)$ is a linear operator, we have,

$\displaystyle\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(\alpha)=\langle
E(\alpha)\mathbf{c}^{\top}g,\mathbf{c}^{\top}g\rangle_{\pi}=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\langle
E(\alpha)g_{i},g_{j}\rangle_{\pi}.$

In particular, for any $a,b\in\mathbb{R}$ such that $a\leq b$,

$\displaystyle\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(b)-\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(a)$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\\{\langle
E(b)g_{i},g_{j}\rangle_{\pi}-\langle E(a)g_{i},g_{j}\rangle_{\pi}\\}$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\langle(E(b)-E(a))g_{i},g_{j}\rangle_{\pi}$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\langle(E(b)-E(a))g_{i},(E(b)-E(a))g_{j}\rangle_{\pi}$
$\displaystyle\leq\sum_{i=1}^{d}\sum_{j=1}^{d}|\mathbf{c}_{i}\mathbf{c}_{j}|\|(E(b)-E(a))g_{i}\|_{L^{2}(\pi)}\|(E(b)-E(a))g_{j}\|_{L^{2}(\pi)}$

where for the third line we use the fact that $E(b)-E(a)$ is also a projection
for $b\geq a$, and we use Cauchy-Schwartz inequality for the last inequality.

We note that for $f\in\mathcal{H}$ and $b\geq a$,
$\|(E(b)-E(a))f\|_{L^{2}(\pi)}^{2}=\Psi_{f,f}(b)-\Psi_{f,f}(a)$ since

$\displaystyle\|(E(b)-E(a))f\|_{L^{2}(\pi)}^{2}=\langle(E(b)-E(a))f,f\rangle_{\pi}=\langle
E(b)f,f\rangle_{\pi}-\langle E(a)f,f\rangle_{\pi}.$

Also, since $\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}$ is an increasing
function, we have,

$\displaystyle 0$
$\displaystyle\leq\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(b)-\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(a)$
$\displaystyle\leq\sum_{i=1}^{d}\sum_{j=1}^{d}|\mathbf{c}_{i}\mathbf{c}_{j}|\sqrt{\Psi_{g_{i},g_{i}}(b)-\Psi_{g_{i},g_{i}}(a)}\sqrt{\Psi_{g_{j},g_{j}}(b)-\Psi_{g_{j},g_{j}}(a)}.$
(S-7)

This implies that for any $(a,b]\subseteq\mathbb{R}$, if
$F_{g_{i}}((a,b])=\Psi_{g_{i},g_{i}}(b)-\Psi_{g_{i},g_{i}}(a)=0$ for all $i$
such that $\mathbf{c}_{i}\neq 0$, then $F_{\mathbf{c}^{\top}g}((a,b])=0$.

Now suppose there exists $x$ such that $x\in{\rm
Supp}(F_{\mathbf{c}^{\top}g})\cap\\{\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})\\}^{c}$. Since $x\in\\{\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})\\}^{c}$ and $\\{\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})\\}^{c}$ is open, we can find $\epsilon>0$ such that
$(x-\epsilon,x+\epsilon)\subseteq\\{\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})\\}^{c}$. Let
$A_{\epsilon}(x)=(x-\epsilon,x+\frac{\epsilon}{2}]$. We have
$F_{g_{i}}(A_{\epsilon}(x))=0$ for all $i$ such that $\mathbf{c}_{i}\neq 0$
because for such $i$,

$\displaystyle F_{g_{i}}(A_{\epsilon}(x))$
$\displaystyle=F_{g_{i}}(A_{\epsilon}(x)\cap\\{\cup_{i;\mathbf{c}_{i}\neq
0}{\rm Supp}(F_{g_{i}})\\}^{c})$ $\displaystyle\leq
F_{g_{i}}(\cap_{i;\mathbf{c}_{i}\neq 0}{\rm Supp}(F_{g_{i}})^{c})$
$\displaystyle\leq F_{g_{i}}({\rm Supp}(F_{g_{i}})^{c})=0.$

In other words, for any $i$ such that $\mathbf{c}_{i}\neq 0$,

$\displaystyle
F_{g_{i}}(A_{\epsilon}(x))=\Psi_{g_{i},g_{i}}(x+\frac{\epsilon}{2})-\Psi_{g_{i},g_{i}}(x-\epsilon)=0.$

Now, since $x\in{\rm Supp}(F_{\mathbf{c}^{\top}g})$,
$F_{\mathbf{c}^{\top}g}(A_{\epsilon}(x))\geq
F_{\mathbf{c}^{\top}g}((x-\frac{\epsilon}{2},x+\frac{\epsilon}{2}))>0$. On the
other hand, using (S1.2), we obtain

$\displaystyle F_{\mathbf{c}^{\top}g}(A_{\epsilon}(x))$
$\displaystyle=\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(x+\frac{\epsilon}{2})-\Psi_{\mathbf{c}^{\top}g,\mathbf{c}^{\top}g}(x-\epsilon)$
$\displaystyle\leq\left\\{\sum_{i=1}^{d}|\mathbf{c}_{i}|\sqrt{\Psi_{g_{i},g_{i}}(x+\frac{\epsilon}{2})-\Psi_{g_{i},g_{i}}(x-\epsilon)}\right\\}^{2}=0$

which is a contradiction. Therefore, we conclude that ${\rm
Supp}(F_{\mathbf{c}^{\top}g})\cap\\{\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})\\}^{c}=\emptyset$, i.e., ${\rm
Supp}(F_{\mathbf{c}^{\top}g})\subseteq\cup_{i;\mathbf{c}_{i}\neq 0}{\rm
Supp}(F_{g_{i}})$. ∎

##  S2 Proofs for Section 4

###  S2.1 Proof for Lemma 1

###### Proof.

First we note that for any $\mathbf{c}\in\mathbb{R}^{d}$,
$\mathbf{c}^{\top}g=\sum_{i=1}^{d}\mathbf{c}_{i}g_{i}$ is in $L_{2}(\pi)$
because

$\displaystyle\int\\{\sum_{i=1}^{d}\mathbf{c}_{i}g_{i}(x)\\}^{2}\pi(dx)$
$\displaystyle\leq\sum_{i=1}^{d}\mathbf{c}_{i}^{2}\int
g_{i}^{2}(x)\pi(dx)+2\sum_{1\leq i<j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}|\sqrt{\int g_{i}^{2}(x)\pi(dx)\int
g_{j}^{2}(x)\pi(dx)}<\infty$

where the last inequality follows from (B.1). Therefore, from Theorem 1, it
suffices to show that $r_{{\mathbf{c}_{M}^{\top}g},M}$ satisfies conditions
(R.1) \- (R.3), for $\gamma=\gamma_{\mathbf{c}^{\top}g}$ in (R.1). Obviously,
$r_{{\mathbf{c}_{M}^{\top}g},M}$ satisfies (R.2) and (R.3).

In the below, we show that $r_{{\mathbf{c}_{M}^{\top}g},M}$ satisfies (R.1),
i.e., we show
$r_{{\mathbf{c}_{M}^{\top}g},M}(k)\to\gamma_{\mathbf{c}^{\top}g}(k)$
$P_{x}$-almost surely for each $k\in\mathbb{Z}$, for each $x\in\mathsf{X}$.
Let the initial condition $x\in\mathsf{X}$ be given. First, we have,

$\displaystyle\gamma_{\mathbf{c}^{\top}g}(k)={\rm
Cov}_{\pi}(\mathbf{c}^{\top}g(X_{0}),\mathbf{c}^{\top}g(X_{k}))=\mathbf{c}^{\top}{\rm
Cov}_{\pi}(g(X_{0}),g(X_{k}))\mathbf{c}=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\gamma_{g_{i},g_{j}}(k)$

where we recall the definition $\gamma_{g_{i},g_{j}}(k)={\rm
Cov}_{\pi}(g_{i}(X_{t}),g_{j}(X_{t+k}))$.

We have,

$\displaystyle r_{{\mathbf{c}_{M}^{\top}g},M}(k)$
$\displaystyle=\frac{1}{M}\sum_{t=0}^{M-1-|k|}\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{Mi}\mathbf{c}_{Mj}\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})$
$\displaystyle=\sum_{1\leq i<j\leq
d}\mathbf{c}_{Mi}\mathbf{c}_{Mj}\left\\{\frac{1}{M}\sum_{t=0}^{M-1-|k|}\\{\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})+\tilde{g}_{i}(X_{t+|k|})\tilde{g}_{j}(X_{t})\\}\right\\}$
$\displaystyle\qquad+\sum_{i=1}^{d}\mathbf{c}_{Mi}^{2}\left\\{\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}_{i}(X_{t})\tilde{g}_{i}(X_{t+|k|})\right\\}$

By Lemma 8 in berg2022efficient, for any $f\in L_{2}(\pi)$ we have

$r_{f,M}(k)=\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{f}(X_{t})\tilde{f}(X_{t+|k|})\to\gamma_{f}(k)$

$P_{x}$-almost surely. Applying this Lemma to $f=g_{i},f=g_{j}$, and
$f=g_{i}+g_{j}$, we have

$\displaystyle r_{g_{i}+g_{j},M}(k)-r_{g_{i},M}(k)-r_{g_{j},M}(k)$
$\displaystyle=\frac{1}{M}\sum_{t=0}^{M-1-|k|}\\{(\tilde{g}_{i}(X_{t})+\tilde{g}_{j}(X_{t}))(\tilde{g}_{i}(X_{t+|k|})+\tilde{g}_{j}(X_{t+|k|}))-\tilde{g}_{i}(X_{t})\tilde{g}_{i}(X_{t+|k|})-\tilde{g}_{j}(X_{t})\tilde{g}_{j}(X_{t+|k|})\\}$
$\displaystyle=\frac{1}{M}\sum_{t=0}^{M-1-|k|}\\{\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})+\tilde{g}_{i}(X_{t+|k|})\tilde{g}_{j}(X_{t})\\}\underset{M\to\infty}{\to}\gamma_{g_{i}+g_{j}}(k)-\gamma_{g_{i}}(k)-\gamma_{g_{j}}(k)$
(S-8)

$P_{x}$-almost surely. Note by definition of $\gamma_{g_{i},g_{j}}$,
$\gamma_{g_{i}+g_{j}}(k)-\gamma_{g_{i}}(k)-\gamma_{g_{j}}(k)=\gamma_{g_{i},g_{j}}(k)+\gamma_{g_{j},g_{i}}(k)$
for each $k\in\mathbb{Z}$. Therefore, w.p $1$,

$\displaystyle\lim_{M\to\infty}r_{{\mathbf{c}_{M}^{\top}g},M}(k)$
$\displaystyle=\sum_{1\leq i<j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}\left\\{\lim_{M\to\infty}\frac{1}{M}\sum_{t=0}^{M-1-|k|}\\{\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})+\tilde{g}_{i}(X_{t+|k|})\tilde{g}_{j}(X_{t})\\}\right\\}$
$\displaystyle\qquad+\sum_{i=1}^{d}\mathbf{c}_{i}^{2}\left\\{\lim_{M\to\infty}\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}_{i}(X_{t})\tilde{g}_{i}(X_{t+|k|})\right\\}$
$\displaystyle=\sum_{1\leq i<j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}\\{\gamma_{g_{i},g_{j}}(k)+\gamma_{g_{j},g_{i}}(k)\\}+\sum_{i=1}^{d}\mathbf{c}_{i}^{2}\gamma_{g_{i},g_{i}}(k)=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\gamma_{g_{i},g_{j}}(k)=\gamma_{\mathbf{c}^{\top}g}(k).$

Thus $r_{{\mathbf{c}_{M}^{\top}g},M}$ satisfies condition (R.1) for
$\gamma=\gamma_{\mathbf{c}^{\top}g}$. Since (R.2)–(R.3) also hold for
$r_{{\mathbf{c}_{M}^{\top}g},M}$, the Lemma follows from Theorem 1. ∎

###  S2.2 Proof for Theorem 2

###### Proof.

When $i=j$, this is a direct consequence of the $\ell_{2}$-convergence of the
momentLS autocovariance sequence estimator for univariate functions of a
chain. Since each $w_{Mi}\to w_{i}\neq 0$, for a sufficiently large $M$,
$w_{Mi}\neq 0$ for $i=1,\dots,d$. Therefore, without loss of generality,
assume $w_{Mi}\neq 0$. When $i\neq j$, by (14),
$\gamma_{g_{i},g_{j}}(k)=\frac{1}{4w_{i}w_{j}}\\{\gamma_{w_{i}g_{i}+w_{j}g_{j}}(k)-\gamma_{w_{i}g_{i}-w_{j}g_{j}}(k)\\}$,
and we have

$\displaystyle\|\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(\cdot;\mathbf{w}_{M})-\gamma_{g_{i},g_{j}}(\cdot)\|_{2}$
$\displaystyle=\|\frac{1}{4w_{Mi}w_{Mj}}\\{\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})-\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})\\}-\frac{1}{4w_{i}w_{j}}(\gamma_{w_{i}g_{i}+w_{j}g_{j}}-\gamma_{w_{i}g_{i}-w_{j}g_{j}})\|_{2}$
$\displaystyle\leq I+II$

where

$\displaystyle I$
$\displaystyle=\|\frac{1}{4w_{Mi}w_{Mj}}\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})-\frac{1}{4w_{i}w_{j}}\gamma_{w_{i}g_{i}+w_{j}g_{j}}\|_{2}$
$\displaystyle II$
$\displaystyle=\|\frac{1}{4w_{Mi}w_{Mj}}\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}-w_{Mj}g_{j}})-\frac{1}{4w_{i}w_{j}}\gamma_{w_{i}g_{i}-w_{j}g_{j}}\|_{2}$

By the inequality

$\displaystyle\|a_{M}b_{M}-ab\|_{2}\leq\|a_{M}(b_{M}-b)\|_{2}+\|(a_{M}-a)b\|_{2},$
(S-9)

for scalar $a_{m},a\in\mathbb{R}$ and $b_{M},b\in\ell_{2}(\mathbb{Z})$, we
have

$\displaystyle
I=\|\frac{1}{4w_{Mi}w_{Mj}}\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})-\frac{1}{4w_{i}w_{j}}\gamma_{w_{i}g_{i}+w_{j}g_{j}}\|_{2}$
(S-10)
$\displaystyle\leq|\frac{1}{4w_{Mi}w_{Mj}}|\|\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})-\gamma_{w_{i}g_{i}+w_{j}g_{j}}\|_{2}+|\frac{1}{4w_{Mi}w_{Mj}}-\frac{1}{4w_{i}w_{j}}|\|\gamma_{w_{i}g_{i}+w_{j}g_{j}}\|_{2}$
(S-11)

Note
$\delta_{ij}=\min\\{\delta_{i},\delta_{j}\\}\leq\Delta(F_{w_{i}g_{i}+w_{j}g_{j}})$
by Proposition 2. Then, by Lemma 1 with
$\mathbf{c}_{M}=w_{Mi}\mathbf{e}_{i}+w_{Mj}\mathbf{e}_{j}$ and
$\mathbf{c}=w_{i}\mathbf{e}_{i}+w_{j}\mathbf{e}_{j}$ where $\mathbf{e}_{j}$ is
the $j$th canonical basis vector, we have
$\|\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}})-\gamma_{w_{i}g_{i}+w_{j}g_{j}}\|_{2}\to
0$, $P_{x}$-a.s for any $x\in\mathsf{X}$. The second term in (S-11) converges
to 0 $P_{x}$-a.s. since $w_{Mi^{\prime}}\overset{a.s.}{\to}w_{i^{\prime}}\neq
0$ $P_{x}$ a.s. for each $i^{\prime}=1,...,d$. Similarly, Term II also
converges to $0$ $P_{x}$-a.s..

The second claim can be shown similarly as follows:

$\displaystyle|\sum_{k\in\mathbb{Z}}\hat{\gamma}^{(\boldsymbol{\delta})}_{M,ij}(k;\mathbf{w}_{M})-\sum_{k\in\mathbb{Z}}\gamma_{g_{i},g_{j}}(k)|$
$\displaystyle\leq|\frac{1}{4w_{Mi}w_{Mj}}\sigma^{2}(\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}+w_{Mj}g_{j}}))-\frac{1}{4w_{i}w_{j}}\sigma^{2}(\gamma_{w_{i}g_{i}+w_{j}g_{j}})|$
$\displaystyle\quad+|\frac{1}{4w_{Mi}w_{Mj}}\sigma^{2}(\Pi_{\delta_{ij}}(r_{M,w_{Mi}g_{i}-w_{Mj}g_{j}}))-\frac{1}{4w_{i}w_{j}}\sigma^{2}(\gamma_{w_{i}g_{i}-w_{j}g_{j}})|$

We see each term converges to $0$ by using a similar inequality as in (S-9)
and applying the asymptotic variance convergence results in Theorem 1 and
Lemma 1. ∎

###  S2.3 Proof for Lemma 2

We prove the following lemma, which encompasses Lemma 2 as a part of its
conclusion.

######  Lemma 3.

Suppose the Markov chain $X$ satisfies (A.1) \- (A.3). Let
$g:\mathsf{X}\to\mathbb{R}^{d}$ be a function which satisfies (B.1). Suppose
$\delta$ is chosen so that $0<\delta\leq\underline{\delta}=\min_{1\leq i\leq
d}\Delta(F_{g_{i}})$. We have

  1. 1.

$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\|\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})-\gamma_{\mathbf{c}^{\top}g}\|_{2}\to
0$, and

  2. 2.

$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}|\sigma^{2}(\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M}))-\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})|\to
0$

$P_{x}$-almost surely for any $x\in\mathsf{X}$.

###### Proof.

Let $\mathbf{c}\in\mathbb{R}^{d}$ be given such that $\|\mathbf{c}\|_{2}=1$.
Let $\hat{F}_{\mathbf{c}^{\top}g,\delta M}$ and $F_{\mathbf{c}^{\top}g}$ be
the representing measures for $\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})$ and
$\gamma_{\mathbf{c}^{\top}g}$ respectively. Note
$\gamma_{\mathbf{c}^{\top}g}\in\mathscr{M}_{\infty}({\delta})\cap\ell_{2}(\mathbb{Z})$
since ${\rm Supp}(F_{\mathbf{c}^{\top}g})\subseteq\cup_{i=1}^{d}{\rm
Supp}(F_{g_{i}})\subseteq[-1+\underline{\delta},1-\underline{\delta}]\subseteq[-1+\delta,1-\delta]$
by Proposition 2. Using Lemma 1 in berg2022efficient,

$\displaystyle\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\|\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})-\gamma_{\mathbf{c}^{\top}g}\|_{2}^{2}$
$\displaystyle\leq\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{-\int\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle
F_{\mathbf{c}^{\top}g}(d\alpha)+\int\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle\hat{F}_{\mathbf{c}^{\top}g,\delta
M}(d\alpha)\right\\}$
$\displaystyle\leq\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{\sup_{\alpha\in[-1+\delta,1-\delta]}|\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\rangle|\right\\}\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\\{F_{\mathbf{c}^{\top}g}([-1,1])+\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\\}$ (S-12)

We first show that

$\displaystyle\limsup_{M\to\infty}\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{\sup_{\alpha\in[-1+\delta,1-\delta]}|\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\rangle|\right\\}=0,\,\,\mbox{$P_{x}$-
almost surely}$

By definition of $r$,

$\displaystyle r_{{\mathbf{c}^{\top}g},M}(k)$
$\displaystyle=\frac{1}{M}\sum_{t=0}^{M-1-|k|}\left(\sum_{i=1}^{d}\mathbf{c}_{i}\tilde{g}_{i}(X_{t})\right)\left(\sum_{j=1}^{d}\mathbf{c}_{j}\tilde{g}_{j}(X_{t+|k|})\right)$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d}\mathbf{c}_{i}\mathbf{c}_{j}\tilde{r}_{ij,M}(k)$
(S-13)

where we define

$\displaystyle\tilde{r}_{ij,M}(k)=\begin{cases}\frac{1}{M}\sum_{t=0}^{M-1-|k|}\tilde{g}_{i}(X_{t})\tilde{g}_{i}(X_{t+|k|})&i=j\\\
\frac{1}{2M}\sum_{t=0}^{M-1-|k|}\\{\tilde{g}_{i}(X_{t})\tilde{g}_{j}(X_{t+|k|})+\tilde{g}_{j}(X_{t})\tilde{g}_{i}(X_{t+|k|})\\}&i\neq
j.\end{cases}$

We also note from (S2.1) that for $i\neq j$,

$\tilde{r}_{ij,M}(k)=\\{r_{g_{i}+g_{j},M}(k)-r_{g_{i},M}(k)-r_{g_{j},M}(k)\\}/2$

and in particular, $|\tilde{r}_{ij,M}(k)|\leq
2^{-1}\\{r_{g_{i}+g_{j},M}(0)+r_{g_{i},M}(0)+r_{g_{j},M}(0)\\}$ $P_{x}$-almost
surely for all $1\leq i,j\leq d$, since $r_{{g},M}(0)\geq|r_{{g},M}(k)|$ a.s.,
for all $k\geq 0$ for empirical autocovariance sequences of the form in (3)
(see, e.g., berg2022efficient, Proposition 7). We note

$\displaystyle\gamma_{\mathbf{c}^{\top}g}(k)=\operatorname{Cov}_{\pi}((\mathbf{c}^{\top}g)(X_{0}),(\mathbf{c}^{\top}g)(X_{k}))=\mathbf{c}^{\top}\operatorname{Cov}_{\pi}(g(X_{0}),g(X_{k}))\mathbf{c}.$

Additionally, for
$\gamma_{g_{i},g_{j}}(k)=\operatorname{Cov}_{\pi}(g_{i}(X_{0}),g_{j}(X_{k}))$,
we have

$\displaystyle|\gamma_{g_{i},g_{j}}(k)|\leq\sqrt{\operatorname{Var}_{\pi}(g_{i}(X_{0}))\operatorname{Var}_{\pi}(g_{j}(X_{k}))}=\sqrt{\gamma_{g_{i}}(0)\gamma_{g_{j}}(0)}.$

For future reference to these bounds, we define

$\displaystyle\check{r}_{ij,M}$
$\displaystyle=0.5\\{r_{g_{i}+g_{j},M}(0)+r_{g_{i},M}(0)+r_{g_{j},M}(0)\\}$
$\displaystyle\check{\gamma}_{ij}$
$\displaystyle=\sqrt{\gamma_{g_{i}}(0)\gamma_{g_{j}}(0)}.$

We note
$\underset{M\to\infty}{\lim}\,\check{r}_{ij,M}=0.5\\{\gamma_{g_{i}+g_{j}}(0)+\gamma_{g_{i}}(0)+\gamma_{g_{j}}(0)\\}$
a.s. due to Assumption (B.1) for all $1\leq i,j\leq d$. In particular,
$\underset{M}{\lim\sup}\,|\check{r}_{ij,M}|<\infty$ a.s. for $1\leq i,j\leq
d$.

Coming back to bounding $\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle$,
we have by (S-13)

$\displaystyle|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|=|\sum_{k\in\mathbb{Z}}\sum_{1\leq
i,j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}(\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k))\alpha^{|k|}|$

We first show that the sequence
$\\{\mathbf{c}_{i}\mathbf{c}_{j}(\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k))\alpha^{|k|}\\}$
is absolutely summable. Then we can exchange the order of summation by Fubini.
We have

$\displaystyle\sum_{k\in\mathbb{Z}}\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}(\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k))\alpha^{|k|}|$
$\displaystyle\leq\sum_{k\in\mathbb{Z}}\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}||\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k)||\alpha|^{|k|}$
$\displaystyle\leq\sum_{k\in\mathbb{Z}}\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}||\alpha|^{|k|}\\{\check{r}_{ij,M}+\check{\gamma}_{ij}\\}$
$\displaystyle=\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}|\\{\check{r}_{ij,M}+\check{\gamma}_{ij}\\}\sum_{k\in\mathbb{Z}}|\alpha|^{|k|}$
$\displaystyle\leq d\frac{1+|\alpha|}{1-|\alpha|}\max_{1\leq i,j\leq
d}\\{\check{r}_{ij,M}+\check{\gamma}_{ij}\\}<\infty.$

where for the final inequality we use the fact that
$\|\mathbf{c}\|_{1}^{2}\leq d\|\mathbf{c}\|_{2}^{2}=d$. Therefore, by
exchanging the summation order, we obtain

$\displaystyle|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|$
$\displaystyle=|\sum_{1\leq i,j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}\sum_{k\in\mathbb{Z}}(\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k))\alpha^{|k|}|$
$\displaystyle\leq\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}||\sum_{k\in\mathbb{Z}}(\tilde{r}_{ij,M}(k)-\gamma_{g_{i},g_{j}}(k))\alpha^{|k|}|$
$\displaystyle\leq\left(\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}|^{2}\right)^{1/2}\left(\sum_{1\leq i,j\leq
d}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|^{2}\right)^{1/2}$
$\displaystyle=\left(\sum_{1\leq i,j\leq
d}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|^{2}\right)^{1/2}$

where for the second to last line we use Holder’s inequality and for the last
equality we use the fact that $\|\mathbf{c}\|_{2}=1$. Now

$\displaystyle\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|\right\\}$
$\displaystyle=\sup_{\alpha\in[-1+\delta,1-\delta]}\left\\{\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|\right\\}$
$\displaystyle\leq\sup_{\alpha\in[-1+\delta,1-\delta]}\left(\sum_{1\leq
i,j\leq
d}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|^{2}\right)^{1/2}$
$\displaystyle\leq\sup_{\alpha\in[-1+\delta,1-\delta]}\left(d^{2}\max_{1\leq
i,j\leq
d}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|^{2}\right)^{1/2}$
$\displaystyle=d\sup_{\alpha\in[-1+\delta,1-\delta]}\max_{1\leq i,j\leq
d}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|$

From (S2.1), we have that $\tilde{r}_{ij,M}(k)\to\gamma_{g_{i},g_{j}}(k)$
almost surely, for each $k\in\mathbb{Z}$. Following the same lines as in the
proof of proposition 8 in berg2022efficient except using the bounds
$|\tilde{r}_{ij}(k)|\leq\check{r}_{ij,M}$ and
$|\gamma_{g_{i},g_{j}}(k)|\leq\check{\gamma}_{ij}$, we obtain

$\displaystyle\lim_{M\to\infty}\sup_{\alpha\in[-1+\delta,1-\delta]}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|=0$
(S-14)

almost surely for each $1\leq i,j\leq d$, and thus

$\displaystyle\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|\right\\}\leq
d\max_{1\leq i,j\leq
d}\left(\sup_{\alpha\in[-1+\delta,1-\delta]}|\langle\tilde{r}_{ij,M}-\gamma_{g_{i},g_{j}},x_{\alpha}\rangle|\right)\to
0$ (S-15)

almost surely, since $d$ is finite and the order of supremums can be
exchanged.

Now we show that the term

$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{F_{\mathbf{c}^{\top}g}([-1,1])+\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\right\\}$

stays bounded. By definition of the representing measure,

$F_{\mathbf{c}^{\top}g}([-1,1])=\gamma_{\mathbf{c}^{\top}g}(0)=\mathbf{c}^{\top}\operatorname{Cov}_{\pi}(g(X_{0}),g(X_{0}))\mathbf{c}=\sum_{1\leq
i,j\leq d}\mathbf{c}_{i}\mathbf{c}_{j}\gamma_{g_{i},g_{j}}(0)$

and

$\sum_{1\leq i,j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}\gamma_{g_{i},g_{j}}(0)\leq\max_{i,j}\check{\gamma}_{ij}\sum_{1\leq
i,j\leq
d}\mathbf{c}_{i}\mathbf{c}_{j}\leq\max_{i,j}\check{\gamma}_{ij}\|\mathbf{c}\|_{1}^{2}\leq
d\max_{i,j}\check{\gamma}_{ij}.$

In particular,
$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}F_{\mathbf{c}^{\top}g}([-1,1])\leq
d\max_{i,j}\check{\gamma}_{ij}<\infty$. Thus, using the conclusion of Lemma 4
below, we have

$\displaystyle\limsup_{M\to\infty}\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\\{F_{\mathbf{c}^{\top}g}([-1,1])+\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\\}\leq d(\max_{i,j}\check{\gamma}_{ij})+C<\infty.$ (S-16)

where $C$ is a constant defined in Lemma 4.

Combining (S-12), (S-15), and (S-16), we conclude
$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\|\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})-\gamma_{\mathbf{c}^{\top}g}\|_{2}\to
0$, $P_{x}$-almost surely.

Having established the sequence convergence result, the asymptotic variance
consistency result can be established following similar lines in the proof of
Lemma 7 in berg2022efficient, except we using the supremum of the total size
of the representing measures

$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\left\\{F_{\mathbf{c}^{\top}g}([-1,1])+\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\right\\}$

for $B$ in their proof and using the supremum $\ell_{2}$ sequence convergence
result, i.e.,

$\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\|\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})-\gamma_{\mathbf{c}^{\top}g}\|_{2}\to
0.$

∎

######  Lemma 4.

Suppose $\delta$ is chosen so that $0<\delta\leq\underline{\delta}=\min_{1\leq
i\leq d}\Delta(F_{g_{i}})$. There exists a constant
$C=\frac{d(2-\delta)/\delta}{\inf_{\alpha,\alpha^{\prime}\in[-1+\delta,1-\delta]}\left\langle
x_{\alpha},x_{\alpha^{\prime}}\right\rangle}\max_{1\leq i,j\leq
d}\check{\gamma}_{ij}<\infty$, depending only on $\delta$ and $g$, such that

$\limsup_{M\to\infty}\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\leq C,$

$P_{x}$\- almost surely, where $\hat{F}_{\mathbf{c}^{\top}g,\delta M}$ is the
representing measure for $\Pi_{\delta}(r_{{\mathbf{c}^{\top}g},M})$ and
$\check{\gamma}_{ij}=\sqrt{\gamma_{g_{i}}(0)\gamma_{g_{j}}(0)}$ for $1\leq
i,j\leq d$.

###### Proof.

Let $\mathbf{c}\in\mathbb{R}^{d}$ such that $\|\mathbf{c}\|_{2}=1$ be given.
We have $\operatorname{Supp}(\hat{F}_{\mathbf{c}^{\top}g,\delta
M})\subseteq[-1+\delta,1-\delta]$ by the definition of the momentLS estimator
$\Pi_{\delta}(\cdot)$. Following the same lines as in the proof of Proposition
9 in berg2022efficient, we obtain

$\displaystyle\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])=\sum_{\alpha\in\operatorname{Supp}(\hat{F}_{\mathbf{c}^{\top}g,\delta
M})}\hat{F}_{\mathbf{c}^{\top}g,\delta
M}(\\{\alpha\\})\leq\frac{\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}\right\rangle|}{\inf_{\alpha,\alpha^{\prime}\in[-1+\delta,1-\delta]}\left\langle
x_{\alpha},x_{\alpha^{\prime}}\right\rangle}$ (S-17)

The denominator is deterministic and does not depend on $\mathbf{c}$. Let
$C_{0}:=\inf_{\alpha,\alpha^{\prime}\in[-1+\delta,1-\delta]}\left\langle
x_{\alpha},x_{\alpha^{\prime}}\right\rangle$. Now we show that for
$\tilde{C}=\frac{d(2-\delta)}{\delta}\max_{1\leq i,j\leq
d}\check{\gamma}_{ij}$, we have

$\displaystyle\limsup_{M\to\infty}\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}\right\rangle|\leq\tilde{C}$ (S-18)

almost surely. We have

$\displaystyle\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}\right\rangle|$
$\displaystyle\leq\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|+\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},\gamma_{\mathbf{c}^{\top}g}\right\rangle|.$

By (S-15),

$\underset{M\to\infty}{\lim\sup}\;\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},r_{{\mathbf{c}^{\top}g},M}-\gamma_{\mathbf{c}^{\top}g}\right\rangle|=0$

almost surely. For the second term

$\displaystyle|\left\langle
x_{\alpha},\gamma_{\mathbf{c}^{\top}g}\right\rangle|$
$\displaystyle=|\sum_{k\in\mathbb{Z}}\alpha^{|k|}\mathbf{c}^{\top}\operatorname{Cov}_{\pi}(g(X_{0}),g(X_{k}))\mathbf{c}|$
$\displaystyle=|\sum_{k\in\mathbb{Z}}\sum_{1\leq i,j\leq
d}\alpha^{|k|}\mathbf{c}_{i}\mathbf{c}_{j}\gamma_{g_{i},g_{j}}(k)|$
$\displaystyle\leq\sum_{1\leq i,j\leq
d}|\mathbf{c}_{i}\mathbf{c}_{j}|\sum_{k\in\mathbb{Z}}|\alpha^{|k|}\gamma_{g_{i},g_{j}}(k)|$
$\displaystyle\leq d(\max_{1\leq i,j\leq
d}\check{\gamma}_{ij})\sum_{k\in\mathbb{Z}}|\alpha|^{|k|}$

where for the second last inequality we use the Fubini-Tonelli Theorem and for
the last inequality we use
$\|\mathbf{c}\|_{1}\leq\sqrt{d}\|\mathbf{c}\|_{2}=\sqrt{d}$. Therefore,

$\displaystyle\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\sup_{\alpha\in[-1+\delta,1-\delta]}|\left\langle
x_{\alpha},\gamma_{\mathbf{c}^{\top}g}\right\rangle|\leq\frac{d(2-\delta)}{\delta}\max_{1\leq
i,j\leq d}\check{\gamma}_{ij}$

Therefore the inequality (S-18) holds with
$\tilde{C}=\frac{d(2-\delta)}{\delta}\max_{1\leq i,j\leq
d}\check{\gamma}_{ij}$. Thus, from (S-17), we have

$\displaystyle\underset{M\to\infty}{\lim\sup}\;\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}\;\hat{F}_{\mathbf{c}^{\top}g,\delta
M}([-1,1])\leq\tilde{C}/C_{0}=\frac{d(2-\delta)/\delta}{\inf_{\alpha,\alpha^{\prime}\in[-1+\delta,1-\delta]}\left\langle
x_{\alpha},x_{\alpha^{\prime}}\right\rangle}\max_{1\leq i,j\leq
d}\check{\gamma}_{ij}$

almost surely. ∎

###  S2.4 Proof for Theorem 3

###### Proof.

Let $\lambda_{1}\geq\dots\geq\lambda_{d}$ be eigenvalues of $\Sigma$ with
corresponding eigenvectors $[\mathbf{U}_{1},\dots,\mathbf{U}_{d}]$ such that
$\Sigma\mathbf{U}_{j}=\lambda_{j}\mathbf{U}_{j}$ for $j=1,\dots,d$. Also,
recall the definition of
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}=\hat{\mathbf{U}}_{M}\hat{\mathbf{\Lambda}}_{\boldsymbol{\delta}}\hat{\mathbf{U}}_{M}^{\top}$,
where $\hat{\mathbf{\Lambda}}_{\boldsymbol{\delta}}$ is a diagonal matrix with
$j$th diagonal entry
$\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M}))$
and $\hat{\mathbf{U}}_{M}=[\hat{\mathbf{U}}_{M1},\dots,\hat{\mathbf{U}}_{Md}]$
is the collection of eigenvectors for the original estimator
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}$, such that
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}\hat{\mathbf{U}}_{Mj}=\hat{\lambda}^{pw}_{j}\hat{\mathbf{U}}_{Mj}$
for $j=1,\dots,d$ and
$\hat{\lambda}^{pw}_{1}\geq\dots\geq\hat{\lambda}^{pw}_{d}$. For notational
simplicity, we denote
$\hat{\lambda}_{j}=\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\hat{\mathbf{U}}_{Mj}^{\top}g},M}))$,
$j=1,\dots,d$ and $\tilde{\Sigma}=\hat{\Sigma}^{(\boldsymbol{\delta})}_{psd}$.
Since
$\tilde{\Sigma}\hat{\mathbf{U}}_{M}=\hat{\mathbf{U}}_{M}\hat{\mathbf{\Lambda}}_{\boldsymbol{\delta}}$,
each $(\hat{\lambda}_{j},\hat{\mathbf{U}}_{Mj})$ is an eigenpair of
$\tilde{\Sigma}$ for $j=1,\dots,d$.

Let $r$ be the number of distinct eigenvalues in
$(\lambda_{1},\dots,\lambda_{d})$ and $\lambda_{(1)}>\dots>\lambda_{(r)}\geq
0$ be the ordered distinct eigenvalues. Define
$I_{j}=\\{k;\lambda_{k}=\lambda_{(j)}\\}$ and let $\mathbf{U}_{I_{j}}$ and
$\hat{\mathbf{U}}_{MI_{j}}$ be sub-matrices consisting of
$\\{\mathbf{U}_{k}\\}_{k\in I_{j}}$ and $\\{\hat{\mathbf{U}}_{Mk}\\}_{k\in
I_{j}}$, respectively, for $j=1,\dots,r$.

We have,

$\displaystyle\|\tilde{\Sigma}-\Sigma\|_{F}$
$\displaystyle=\|\sum_{j=1}^{d}\hat{\lambda}_{j}\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top}-\sum_{j=1}^{d}\lambda_{j}\mathbf{U}_{j}\mathbf{U}_{j}^{\top}\|_{F}$
$\displaystyle=\|\sum_{j=1}^{d}(\hat{\lambda}_{j}-\lambda_{j})\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top}+\sum_{j=1}^{d}\lambda_{j}(\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top}-\mathbf{U}_{j}\mathbf{U}_{j}^{\top})\|_{F}$
$\displaystyle\leq\sum_{j=1}^{d}|\hat{\lambda}_{j}-\lambda_{j}|+\sum_{j=1}^{r}\lambda_{(j)}\|\hat{\mathbf{U}}_{MI_{j}}\hat{\mathbf{U}}_{MI_{j}}^{\top}-\mathbf{U}_{I_{j}}\mathbf{U}_{I_{j}}^{\top}\|_{F}$
(S-19)

where for the last inequality we use the fact
$\|\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top}\|_{F}=\sqrt{\textrm{tr}(\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top}\hat{\mathbf{U}}_{Mj}\hat{\mathbf{U}}_{Mj}^{\top})}=1$
for $j=1,\dots,d$.

For the second term in (S-19),

$\displaystyle\|\hat{\mathbf{U}}_{MI_{j}}\hat{\mathbf{U}}_{MI_{j}}^{\top}-\mathbf{U}_{I_{j}}\mathbf{U}_{I_{j}}^{\top}\|_{F}=\frac{1}{\sqrt{2}}\|\textrm{Sin}\,\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})\|_{F}^{2}.$

where $\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})$ denotes the
$|I_{j}|$ by $|I_{j}|$ diagonal matrix whose diagonal entries are the
principal angles of the column spaces generated by $\mathbf{U}_{I_{j}}$ and
$\hat{\mathbf{U}}_{MI_{j}}$ respectively, and
$\sin\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})$ is defined entry-
wise. More specifically,

$\sin\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})={\rm
diag}\\{\sin(\cos^{-1}(\sigma_{1})),\dots,\sin(\cos^{-1}(\sigma_{|I_{j}|}))\\}$

where $\sigma_{1}\geq\dots\geq\sigma_{|I_{j}|}$ are the singular values of
$\hat{\mathbf{U}}_{MI_{j}}^{\top}\mathbf{U}_{I_{j}}$.

By a variant of the Davis-Kahan theorem (Theorem 2 in yu2015useful), for each
$j=1,\dots,r$,

$\|\sin\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})\|_{F}\leq\frac{2\min\left(|I_{j}|^{1/2}\|\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}-\Sigma\|_{\mathrm{op}},\|\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}-\Sigma\|_{\mathrm{F}}\right)}{\min\left(\lambda_{(j-1)}-\lambda_{(j)},\lambda_{(j)}-\lambda_{(j+1)}\right)}.$

where we define $\lambda_{(0)}=\infty$ and $\lambda_{(r+1)}=-\infty$. In
particular, since $d$ is fixed and
$\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw,ij}\to\Sigma_{ij}$ for each $i,j$,
$P_{x}$-a.s., we have
$\|\sin\Theta(\mathbf{U}_{I_{j}},\hat{\mathbf{U}}_{MI_{j}})\|_{F}\to 0$ for
each $j$, $P_{x}$-a.s.

For the first term in (S-19), we have for each $i=1,\dots,d$,

$\displaystyle|\hat{\lambda}_{i}-\lambda_{i}|$
$\displaystyle\leq|\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\hat{\mathbf{U}}_{Mi}^{\top}g},M}))-\sigma^{2}(\gamma_{\hat{\mathbf{U}}_{Mi}^{\top}g})|+|\sigma^{2}(\gamma_{\hat{\mathbf{U}}_{Mi}^{\top}g})-\sigma^{2}(\gamma_{\mathbf{U}_{i}^{\top}g})|$
$\displaystyle\leq\sup_{\mathbf{c};\|\mathbf{c}\|_{2}=1}|\sigma^{2}(\Pi_{\underline{\delta}}(r_{{\mathbf{c}^{\top}g},M}))-\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})|+|\sigma^{2}(\gamma_{\hat{\mathbf{U}}_{Mi}^{\top}g})-\sigma^{2}(\gamma_{\mathbf{U}_{i}^{\top}g})|$
(S-20)

The first term in (S-20) converges to $0$ almost surely due to Lemma 2 since
$\underline{\delta}=\min_{1\leq i\leq d}\delta_{i}\leq\min_{1\leq i\leq
d}\Delta(F_{g_{i}})$ by the condition of
$\boldsymbol{\delta}=[\delta_{1},\dots,\delta_{d}]$ in Theorem 2. For the
second term in (S-20), for each $\mathbf{c}\in\mathbb{R}^{d}$, we have
$\sigma^{2}(\gamma_{\mathbf{c}^{\top}g})=\mathbf{c}^{\top}\Sigma\mathbf{c}$
and

$\displaystyle|\sigma^{2}(\gamma_{\hat{\mathbf{U}}_{Mi}^{\top}g})-\sigma^{2}(\gamma_{\mathbf{U}_{i}^{\top}g})|=|\hat{\mathbf{U}}_{Mi}^{\top}\Sigma\hat{\mathbf{U}}_{Mi}-\mathbf{U}_{i}^{\top}\Sigma\mathbf{U}_{i}|$

Let $j$ be the index such that $\lambda_{i}=\lambda_{(j)}$. That is $i\in
I_{j}=\\{k;\lambda_{k}=\lambda_{(j)}\\}$. We note that for any orthogonal
matrix $\mathbf{O}\in\mathbb{R}^{|I_{j}|\times|I_{j}|}$, a column of the
matrix $\mathbf{U}_{I_{j}}\mathbf{O}\in\mathbb{R}^{d\times|I_{j}|}$ is also an
eigenvector corresponding to $\lambda_{(j)}$ since for $l=1,\dots,|I_{j}|$

$\displaystyle\Sigma(\mathbf{U}_{I_{j}}\mathbf{O})_{l}=\Sigma\left(\sum_{k=1}^{|I_{j}|}\mathbf{U}_{I_{j},k}\mathbf{O}_{kl}\right)=\sum_{k=1}^{|I_{j}|}\Sigma\mathbf{U}_{I_{j},k}\mathbf{O}_{kl}=\lambda_{(j)}\sum_{k=1}^{|I_{j}|}\mathbf{U}_{I_{j},k}\mathbf{O}_{kl}=\lambda_{(j)}(\mathbf{U}_{I_{j}}\mathbf{O})_{l}$

Moreover, by Theorem 2 in yu2015useful, for each $j=1,\dots,r$ there exists an
orthogonal matrix $\hat{\mathbf{O}}_{Mj}\in\mathbb{R}^{|I_{j}|\times|I_{j}|}$
such that

$\displaystyle\|\mathbf{U}_{I_{j}}-\hat{\mathbf{U}}_{MI_{j}}\hat{\mathbf{O}}_{Mj}\|_{F}\leq\frac{2^{3/2}\min\left(|I_{j}|^{1/2}\|\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}-\Sigma\|_{\mathrm{op}},\|\hat{\Sigma}^{(\boldsymbol{\delta})}_{pw}-\Sigma\|_{\mathrm{F}}\right)}{\min\left(\lambda_{(j-1)}-\lambda_{(j)},\lambda_{(j)}-\lambda_{(j+1)}\right)}.$
(S-21)

Let us define $\tilde{\mathbf{O}}_{Mj}=\hat{\mathbf{O}}_{Mj}^{-1}$, and let
$l$ be the column index for $\hat{\mathbf{U}}_{Mi}$ in
$\hat{\mathbf{U}}_{MI_{j}}$, i.e.,
$\hat{\mathbf{U}}_{MI_{j}}\mathbf{e}_{l}=\hat{\mathbf{U}}_{Mi}$. Note that
$\tilde{\mathbf{O}}_{Mj}$ is also an orthogonal matrix, and therefore
$(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\mathbf{e}_{l}$ is an eigenvector
of $\Sigma$ whose corresponding eigenvalue is $\lambda_{(j)}=\lambda_{i}$. In
particular,
$\lambda_{i}=\mathbf{U}_{i}^{\top}\Sigma\mathbf{U}_{i}=\mathbf{e}_{l}^{\top}(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})^{\top}\Sigma(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\mathbf{e}_{l}$.
Therefore, we have,

$\displaystyle|\sigma^{2}(\gamma_{\tilde{\mathbf{U}}_{i}^{\top}g})-\sigma^{2}(\gamma_{\mathbf{U}_{i}^{\top}g})|$
$\displaystyle=|\hat{\mathbf{U}}_{Mi}^{\top}\Sigma\hat{\mathbf{U}}_{Mi}-\mathbf{U}_{i}^{\top}\Sigma\mathbf{U}_{i}|$
$\displaystyle=|\mathbf{e}_{l}^{\top}\hat{\mathbf{U}}_{MI_{j}}^{\top}\Sigma\hat{\mathbf{U}}_{MI_{j}}\mathbf{e}_{l}-\mathbf{e}_{l}^{\top}(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})^{\top}\Sigma(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\mathbf{e}_{l}|$
$\displaystyle\leq\|\hat{\mathbf{U}}_{MI_{j}}\mathbf{e}_{l}-(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\mathbf{e}_{l}\|_{2}\\{\|\Sigma\hat{\mathbf{U}}_{MI_{j}}\mathbf{e}_{l}\|_{2}+\|\Sigma(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\mathbf{e}_{l}\|_{2}\\}$
$\displaystyle\leq
2\|\hat{\mathbf{U}}_{MI_{j}}-(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\|_{F}\|\Sigma\|_{F}$

where the first inequality uses the fact that

$\displaystyle|\mathbf{a}^{\top}\Sigma\mathbf{a}-\mathbf{b}^{\top}\Sigma\mathbf{b}|\leq\|\mathbf{a}-\mathbf{b}\|_{2}\\{\|\Sigma\mathbf{a}\|_{2}+\|\Sigma\mathbf{b}\|_{2}\\}$

and the second inequality uses the fact that for any
$\mathbf{A}\in\mathbb{R}^{d\times r}$ and $\mathbf{b}\in\mathbb{R}^{r}$,
$\|\mathbf{A}\mathbf{b}\|_{2}\leq\|\mathbf{A}\|_{F}\|\mathbf{b}\|_{2}$.

Now,

$\displaystyle\|\hat{\mathbf{U}}_{MI_{j}}-(\mathbf{U}_{I_{j}}\tilde{\mathbf{O}}_{Mj})\|_{F}=\|\hat{\mathbf{U}}_{MI_{j}}\hat{\mathbf{O}}_{Mj}-(\mathbf{U}_{I_{j}})\|_{F}\underset{M\to\infty}{\to}0$

$P_{x}$-almost surely by (S-21), which shows that
$|\hat{\lambda}_{j}-\lambda_{j}|\underset{M\to\infty}{\to}0$ $P_{x}$-almost
surely for each $j=1,\dots,d$. This implies
$\sum_{j=1}^{d}|\hat{\lambda}_{j}-\lambda_{j}|\underset{M\to\infty}{\to}0$,
$P_{x}$-a.s., since $d$ is finite. ∎

## References

  * Berg and Song [2022] Stephen Berg and Hyebin Song.  Efficient shape-constrained inference for the autocovariance sequence from a reversible Markov chain.  _arXiv preprint arXiv:2207.12705_ , 2022.
  * Meyn and Tweedie [2009] Sean Meyn and Richard L. Tweedie.  _Markov Chains and Stochastic Stability_.  Cambridge University Press, Cambridge, second edition, 2009.
  * Heidelberger and Welch [1981] Philip Heidelberger and Peter D Welch.  A spectral method for confidence interval generation and run length control in simulations.  _Communications of the ACM_ , 24:233–245, 1981.
  * Jones et al. [2006] Galin L Jones, Murali Haran, Brian S Caffo, and Ronald Neath.  Fixed-width output analysis for Markov chain Monte Carlo.  _Journal of the American Statistical Association_ , 101:1537–1547, 2006.
  * Dai and Jones [2017] Ning Dai and Galin L. Jones.  Multivariate initial sequence estimators in Markov chain Monte Carlo.  _Journal of Multivariate Analysis_ , 159:184–199, 2017\.
  * Vats et al. [2019] Dootika Vats, James M Flegal, and Galin L Jones.  Multivariate output analysis for Markov chain Monte Carlo.  _Biometrika_ , 106:321–337, 2019.
  * Dellaportas and Kontoyiannis [2012] Petros Dellaportas and Ioannis Kontoyiannis.  Control variates for estimation based on reversible Markov chain Monte Carlo samplers.  _Journal of the Royal Statistical Society: Series B_ , 74:133–161, 2012.
  * Berg et al. [2019] Stephen Berg, Jun Zhu, and Murray K. Clayton.  Control variates and Rao-Blackwellization for deterministic sweep Markov chains.  _arXiv preprint arXiv:1912.06926_ , 2019.
  * Belomestny et al. [2020] Denis Belomestny, Leonid Iosipoi, Eric Moulines, Alexey Naumov, and Sergey Samsonov.  Variance reduction for Markov chains with application to MCMC.  _Statistics and Computing_ , 30:973–997, 2020.
  * South et al. [2022] Leah F South, Marina Riabiz, Onur Teymur, and Chris J Oates.  Postprocessing of MCMC.  _Annual Review of Statistics and Its Application_ , 9:529–555, 2022.
  * Anderson [1971] T. W. Anderson.  _The statistical analysis of time series_.  John Wiley & Sons, Inc., New York-London-Sydney, 1971.
  * Priestley [1981] M.B. Priestley.  _Spectral Analysis and Time Series_.  Number Bd. 1-2 in Probability and mathematical statistics : a series of monographs and textbooks. Academic Press, 1981.
  * Flegal and Jones [2010] James M. Flegal and Galin L. Jones.  Batch means and spectral variance estimators in Markov chain Monte Carlo.  _The Annals of Statistics_ , 38:1034 – 1070, 2010.
  * Chakraborty et al. [2022] Saptarshi Chakraborty, Suman K Bhattacharya, and Kshitij Khare.  Estimating accuracy of the MCMC variance estimator: Asymptotic normality for batch means estimators.  _Statistics & Probability Letters_, 183:109337, 2022.
  * Geyer [1992] Charles J. Geyer.  Practical Markov chain Monte Carlo.  _Statistical Science_ , 7:473–483, 1992.
  * Damerdji [1991] Halim Damerdji.  Strong consistency and other properties of the spectral variance estimator.  _Management Science_ , 37:1424–1440, 1991.
  * Stan Development Team [2019] Stan Development Team.  Stan reference manual, version 2.28, 2019.
  * Vats et al. [2018] Dootika Vats, James M. Flegal, and Galin L. Jones.  Strong consistency of multivariate spectral variance estimators in Markov chain Monte Carlo.  _Bernoulli_ , 24:1860 – 1909, 2018.
  * Kosorok [2000] Michael R. Kosorok.  Monte Carlo error estimation for multivariate Markov chains.  _Statistics & Probability Letters_, 46:85–93, 2000.
  * Roberts and Rosenthal [1997] Gareth Roberts and Jeffrey Rosenthal.  Geometric ergodicity and hybrid Markov chains.  _Electronic Communications in Probability_ , 2:13–25, 1997\.
  * Kontoyiannis and Meyn [2012] Ioannis Kontoyiannis and Sean P Meyn.  Geometric ergodicity and the spectral gap of non-reversible Markov chains.  _Probability Theory and Related Fields_ , 154:327–339, 2012\.
  * Roberts and Rosenthal [1998] Gareth O Roberts and Jeffrey S Rosenthal.  Markov-chain Monte Carlo: some practical implications of theoretical results.  _The Canadian Journal of Statistics/La Revue Canadienne de Statistique_ , 26:5–20, 1998.
  * Jones [2004] Galin L. Jones.  On the Markov chain central limit theorem.  _Probability Surveys_ , 1:299 – 320, 2004.
  * Stein and Shakarchi [2009] Elias M Stein and Rami Shakarchi.  _Real analysis: measure theory, integration, and Hilbert spaces_.  Princeton University Press, 2009.
  * Hoffman et al. [2014] Matthew D Hoffman, Andrew Gelman, et al.  The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15:1593–1623, 2014\.
  * Ōsawa [1988] Hideo Ōsawa.  Reversibility of first-order autoregressive processes.  _Stochastic processes and their applications_ , 28:61–69, 1988.
  * R Core Team [2020] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2020.
  * Flegal et al. [2021] James M. Flegal, John Hughes, Dootika Vats, Ning Dai, Kushagra Gupta, and Uttiya Maji.  _mcmcse: Monte Carlo Standard Errors for MCMC_ , 2021.  R package version 1.5-0.
  * Liu et al. [2021] Ying Liu, Dootika Vats, and James M Flegal.  Batch size selection for variance estimators in MCMC.  _Methodology and Computing in Applied Probability_ , pages 1–29, 2021\.
  * Martin et al. [2011] Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park.  Mcmcpack: Markov chain Monte Carlo in R.  _Journal of Statistical Software_ , 42:1–21, 2011.
  * Folland [1999] Gerald B Folland.  _Real analysis: modern techniques and their applications_ , volume 40.  John Wiley & Sons, 1999.
  * Yu et al. [2015] Yi Yu, Tengyao Wang, and Richard J Samworth.  A useful variant of the Davis–Kahan theorem for statisticians.  _Biometrika_ , 102:315–323, 2015.