# Strong Invariance Principles for Ergodic Markov Processes

Ardjen Pengel  Delft Institute of Applied Mathematics, Delft University of
Technology Joris Bierkens  Delft Institute of Applied Mathematics, Delft
University of Technology

###### Abstract

Strong invariance principles describe the error term of a Brownian
approximation of the partial sums of a stochastic process. While these strong
approximation results have many applications, the results for continuous-time
settings have been limited. In this paper, we obtain strong invariance
principles for a broad class of ergodic Markov processes. The main results
rely on ergodicity requirements and an application of Nummelin splitting for
continuous-time processes. Strong invariance principles provide a unified
framework for analysing commonly used estimators of the asymptotic variance in
settings with a dependence structure. We demonstrate how this can be used to
analyse the batch means method for simulation output of Piecewise
Deterministic Monte Carlo samplers. We also derive a fluctuation result for
additive functionals of ergodic diffusions using our strong approximation
results.

Keywords: Strong invariance priciple, piecewise deterministic Markov
processes, diffusions, asymptotic variance estimation.

##  1 Introduction

Let $X=(X_{k})_{k\in\mathbb{N}}$ be a stochastic sequence defined on a common
probability space and consider the partial sum process $S_{n}$, given by
$S_{n}=\sum_{k=1}^{n}X_{k}$. We say that a strong invariance principle holds
for $X$ if there exist a probability space $(\Omega,\mathcal{F},\mathbb{P})$
on which we can construct a sequence of random variables
$X^{\prime}=(X^{\prime}_{k})_{k\in\mathbb{N}}$ and a Brownian motion
$W=(W(t))_{t\geq 0}$, such that $X$ and $X^{\prime}$ are equal in law and

$\absolutevalue{S_{n}^{\prime}-\mu n-\sigma
W(n)}=\mathscr{O}(\psi_{n})\quad\textrm{a.s.,}$

where $S_{n}^{\prime}$ denotes the partial sum process of $X^{\prime}$, $\mu$
and $\sigma$ are finite constants determined by the law of the process,
$\mathscr{O}$ describes the asymptotic regime, and $\psi_{n}$ the
corresponding approximation error. More specifically, if $S=(S_{t})_{t\geq 0}$
denotes a stochastic process and $\psi=(\psi_{t})_{t\geq 0}$ is some positive
sequence, we write

$S_{T}=o(\psi_{T})\ \textrm{a.s.}\ \quad\textrm{and}\quad S_{T}=O(\psi_{T})\
\textrm{a.s.}$

to denote

$\mathbb{P}\left(\lim_{T\rightarrow\infty}{S_{T}}\big{/}{\psi_{T}}=0\right)=1\
\quad\textrm{and}\quad\mathbb{P}\left(\limsup_{T\rightarrow\infty}\
{\absolutevalue{S_{T}}}\big{/}{\psi_{T}}<\infty\right)=1\ $

respectively. For notational convenience we will usually make no distinction
between $X$ and $X^{\prime}$.

For a sequence of independent, identically distributed random variables with
mean zero and unit variance, the Komlós-Major-Tusnády approximation [46; 47]
asserts that if $E\absolutevalue{X_{1}}^{p}<\infty$ for some $p>2$, then on a
suitably enriched probability space, we can construct a Brownian motion
$W=\\{W(t),t\geq 0\\}$ such that

$S_{n}=W(n)+o(n^{1/p})\quad\textrm{a.s.}$ (1)

If we can additionally assume that the moment-generating function exists in an
area around zero, i.e., $\mathbb{E}e^{t\absolutevalue{X}}<\infty$ for some
$t>0$, then

$S_{n}=W(n)+O(\log n)\quad\textrm{a.s.}$ (2)

Furthermore, if only existence of the the second moment is assumed, Major [54]
showed that there exists a sequence $t_{n}\sim n$ such that

$S_{n}=W(t_{n})+o(n^{1/2})\quad\textrm{a.s.}$ (3)

The error terms appearing in the strong invariance principles (1), (2), and
(3) are optimal. The approximation error appearing in the strong invariance
principle also quantifies the convergence rate in the functional central limit
theorem, as shown in Csörgö and Horváth [17; Theorem 1.16 and Theorem 1.17].
These strong approximation results are powerful tools used to obtain numerous
results in both probability and statistics as seen in, e.g., Csörgö and Hall
[20], Csörgö and Révész [19], Parzen [65], and Shorack and Wellner [74].

Naturally, it is of great interest to extend these results beyond the i.i.d.
setting. Berkes et al. [5] gives an extensive overview of invariance
principles for dependent sequences. Moreover, they show that for stationary
sequences, which can be represented as a function of an i.i.d sequence, the
Komlós-Major-Tusnády bound (1) can be attained given the appropriate decay of
coupling coefficients. Employing this result, Cuny et al. [21] show that
partial sums of functions of random iterates also inherit the optimal bound.
For most settings with dependent variables attaining the optimal rates still
remains an open problem. Berkes et al. [4] considers the class of stationary
and weakly $m$-dependent processes and derive a strong invariance principle
attaining a rate close to the Komlós-Major-Tusnády bound.

In Markovian settings, strong approximation results were obtained by Cuny et
al. [21], Csáki and Csörgő [16], Vats et al. [80], and Merlevède et al. [55],
among others. The strong invariance principle of Merlevède et al. [55] attains
the Komlós-Major-Tusnády bound given in (2). The results of Csáki and Csörgő
[16] and Merlevède et al. [55] are established through an application of
Nummelin splitting, introduced in the seminal papers of Athreya and Ney [2]
and Nummelin [60]. Provided that the transition operator of the chain
satisfies a one-step minorisation condition, a bivariate process can be
constructed, on an enlarged probability space, such that this process
possesses a recurrent atom and the first coordinate of the constructed process
is equal in law to the original Markov chain. Since this constructed process
has a recurrent atom, it inherits a regenerative structure, i.e., if $X$
denotes the considered Markov chain and $R$ denotes the hitting time of the
recurrent atom, then

  1. i.

$\left(X_{k+R}\right)_{k\in\mathbb{N}}$ is independent of $(X_{k})_{0\leq
k\leq R}$ and $R$

  2. ii.

$\left(X_{k+R}\right)_{k\in\mathbb{N}}$ is equal in law to
$(X_{k})_{k\in\mathbb{N}}$.

Consequently, the chain can be divided into independent identically
distributed cycles, and the Komlós-Major-Tusnády approximation can be applied
to obtain a strong invariance principle.

Strong invariance principles for Markov chains are useful tools for analysing
estimators of the asymptotic variance of Markov Chain Monte Carlo (MCMC)
sampling algorithms. Damerdji [22; 23], Flegal et al. [33], and Vats et al.
[80] show strong consistency of the batch means and spectral variance
estimators for MCMC simulation output using the appropriate strong
approximation results.

Recently, there has been growing interest in Monte Carlo algorithms based on
Piecewise Deterministic Markov Processes (PDMPs), see for example Bierkens et
al. [8] and Bouchard-Côté et al. [10]. The main appeal of these processes is
their non-reversible nature. Sampling algorithms based on non-reversible
processes can offer advantageous properties to their reversible MCMC
counterparts in terms of both asymptotic variance and convergence to
equilibrium. In order for many useful results regarding estimation of the
asymptotic variance of Markov chain simulation output to carry over to PDMP-
based methods, it is required that a strong invariance principle holds for the
underlying continuous-time process.

In this paper, we obtain strong approximation results for a broad class of
(continuous-time) ergodic Markov processes. Firstly, we show that the strong
invariance principle given in Theorem 3.1 can be obtained directly through
ergodicity and moment conditions. However, the resulting error rate is not
explicit and therefore less convenient to work with. A natural approach for
obtaining a more refined strong invariance principle would be through
regenerative properties of the process. If the process we consider is
regenerative then the approximation results of Csáki and Csörgő [16] would
hold. However, it is in general not possible to show that the transition
semigroup satisfies a minorisation condition such that a regenerative
structure can be obtained. The resolvent chain, on the other hand, does
satisfy a one-step minorisation condition. Utilising this result, Löcherbach
and Loukianova [52] extend the concept of Nummelin splitting to Harris
recurrent Markov processes. Hence we can redefine the process such that it is
embedded in a richer process which is endowed with a recurrent atom. Although
the resulting cycles are not independent and we therefore do not have
regeneration in the classic sense, we do obtain short range dependence.
Therefore we can utilise the approximation results of Berkes et al. [4] to
obtain a strong invariance principle attaining a convergence rate of order
$O(T^{1/4}\log(T))$. This result is formulated in Theorem 3.5. Although the
nearly optimal bound of Berkes et al. [4] does not carry over, to the best of
our knowledge, there are currently no approaches established that lead to
superior rates for the class of processes we consider.

Furthermore, we show that the required conditions for our obtained strong
invariance principle can be simplified to more natural ergodicity and moment
conditions. Finally, we discuss some applications of our obtained strong
invariance principles. We demonstrate how strong approximation results can be
utilised for analysing the batch means estimator of the asymptotic variance of
continuous-time Monte Carlo samplers and we show that for diffusion processes
the magnitude of increments can be described.

This article is organised as follows. In Section 2, we describe the setting
and review Nummelin splitting in continuous time as introduced by Löcherbach
and Loukianova [52] and discuss other applicable results. In Section 3, the
main results of the paper are given. In Section 4, we give a brief
introduction of Piecewise Deterministic Markov processes and discuss the
estimation of the asymptotic variance. Section 5 illustrates the applicability
of our results to diffusion processes. In Section 6, the proofs of the main
results are given.

##  2 Preliminaries on Markov processes and Nummelin splitting in continuous
time

Let $X=(X_{t})_{t\geq 0}$ be a stochastic process defined on a filtered
probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq
0},\mathbb{P}_{x})$, with Polish state space $(E,\mathscr{E})$ and initial
value $X_{0}=x.$ We consider the case where $X$ is a positive Harris recurrent
strong Markov process with transition semigroup given by $(P_{t})_{t\geq 0}$
with finite invariant measure $\pi$. By definition of positive Harris
recurrence, $\pi$ can be normalised to be a probability measure and we have
that

$\pi(A)>0\implies\mathbbm{P}_{x}\left(\int_{0}^{\infty}\mathbbm{1}_{\\{X_{s}\in
A\\}}ds=\infty\right)=1,\quad x\in E.$ (4)

Throughout this paper we will additionally require ergodicity of the
considered processes. We say that a Markov process $X$ is _ergodic_ with
convergence rate $\Psi$ if

$\norm{P_{t}(x,\cdot)-\pi}_{TV}\leq V(x)\Psi(t),\quad\textrm{for all}\ x\in E\
\textrm{and}\ t\geq 0,$ (5)

where $V$ is some positive $\pi$-integrable function and $\Psi$ some positive
rate tending to zero. Furthermore, a process is called polynomially or
exponentially ergodic if $\Psi$ decays accordingly. For a more thorough
discussion of these definitions we refer to Meyn and Tweedie [56].

The resolvent chain $\bar{X}=(\bar{X}_{n})_{n\geq 0}$ is obtained by observing
the process at independent exponential times, i.e., $\bar{X}_{n}:=X_{T_{n}}$
for $n\geq 0$. Here $(T_{n})_{n\geq 0}$ denote the sampling times at which we
observe the process $X$, which are defined as $T_{0}:=0$ and
$T_{n}:=\sum_{k=1}^{n}\sigma_{k}$, where $(\sigma_{k})_{k\geq 1}$ denote a
sequence of i.i.d. exponential random variables. The resolvent chain will
inherit positive Harris recurrence from the original process. The transition
kernel of the process $\bar{X}=(\bar{X}_{n})_{n\in\mathbb{N}}$ is given by

$U(x,A)=\int_{0}^{\infty}P_{t}(x,A)e^{-t}dt,$ (6)

and satisfies the one-step minorisation condition, see for example Höpfner and
Löcherbach [42] or Revuz [68],

$U(x,A)\geq h\otimes\nu(x,A),$ (7)

where $h\otimes\nu(x,A)=h(x)\nu(A)$, with $h(x)=\alpha\mathbbm{1}_{C}(x)$ for
some $\alpha\in(0,1)$, a measurable set $C$ with $\pi(C)>0$, and $\nu(\cdot)$
a probability measure equivalent to $\pi(\cdot\cap C)$.

The minorisation condition of the resolvent chain motivates the introduction
of the kernel $K((x,u),dy):E\times[0,1]\rightarrow E$ given by

$K((x,u),dy)=\left\\{\begin{array}[]{ll}\nu(dy)&\mbox{for}\ (x,u)\in
C\times[0,\alpha]\vskip 12.0pt plus 4.0pt minus 4.0pt\\\ W(x,dy)&\mbox{for}\
(x,u)\in C\times(\alpha,1]\vskip 12.0pt plus 4.0pt minus 4.0pt\\\
U(x,dy)&\mbox{for}\ x\notin C\end{array}\right.,$ (8)

where the residual kernel $W(x,dy)$ is defined as

$W(x,dy)=\frac{U(x,dy)-\alpha\nu(dy)}{1-\alpha}.$ (9)

Since the resolvent chain is also positive Harris recurrent, it will hit $C$
infinitely often. Given that the resolvent chain has hit $C$, with probability
$\alpha$ the chain will move independently of its past according to the small
measure $\nu$ and with probability $(1-\alpha)$ it will move according to the
residual kernel $W$. By the Borel-Cantelli lemma the residual chain will
generate according to $\nu$ infinitely often. Let $R_{k}$ denote the $k$-th
time that the resolvent chain moves according to $\nu$. The randomised
stopping times $(R_{k})_{k}$ serve as regeneration epochs for the resolvent;
for every $k$, $\bar{X}_{R_{k}}$ has law $\nu$ and is independent of both its
past and of $R_{k}$. The implied regenerative properties that the process $X$
obtains through its resolvent are made explicit with the approach of
Löcherbach and Loukianova [52]. Their framework requires the following
regularity conditions on the transition semigroup of the process $X$:

######  Assumption 1.

  * (i)

The semigroup $(P_{t})_{t\geq 0}$ is Feller, i.e., for every
$A\in\mathscr{E}$, the mapping $x\longmapsto P_{t}(x,A)$ is bounded.

  * (ii)

There exists a $\sigma$-finite measure $\Lambda$ on $(E,\mathscr{E})$ such
that for every $t>0$, $P_{t}(x,dy)=p_{t}(x,dy)\Lambda(dy)$, with
$(t,x,y)\longmapsto p_{t}(x,y)$ jointly measurable.

At the so-called sampling times of the process $X$, we can apply the Nummelin
splitting technique to the resolvent chain. We then fill in the original
process between the sampling times. Following this procedure, Löcherbach and
Loukianova [52] construct on an extended probability space a process $Z$ with
state space $E\times[0,1]\times E$, that admits a recurrent atom. The first
coordinate of $Z$ has the same law as the original process $X$, the second
coordinate denotes the auxiliary variables employed in order to generate draws
from the resolvent chain via the splitting procedure, and the third coordinate
corresponds to the subsequent values of the resolvent chain.

The process $Z=(Z^{1}_{t},Z^{2}_{t},Z^{3}_{t})_{t\geq 0}$ can be constructed
according to the following procedure. Firstly, let $Z^{1}_{0}=X_{0}=x$.
Independently from $Z_{1}$ generate $Z_{0}^{2}\sim U[0,1]$, where $U[0,1]$
denotes the uniform distribution on the unit interval. Given
$\\{Z_{0}^{2}=u\\}$, draw $Z_{0}^{3}$ according to $K((x,u),dx^{\prime})$.
Then inductively for $n\geq 1$, on $Z_{0}=(x,u,x^{\prime})$:

  1. I.

Choose $\sigma_{n+1}$ according to

$\left(\frac{p_{t}(x,x^{\prime})}{u(x,x^{\prime})}\mathbbm{1}_{\\{0<u(x,x^{\prime})<\infty\\}}+\mathbbm{1}_{\\{u(x,x^{\prime})\in\\{0,\infty\\}\\}}\right)e^{-t}dt\
\ \textrm{on}\ \ \mathbb{R}_{+}.$ (10)

The next sampling time $T_{n+1}$ is given by $T_{n}+\sigma_{n+1}$.

  2. II.

On $\\{\sigma_{n+1}=t\\}$, put $Z^{2}_{T_{n}+s}:=u$ and
$Z^{3}_{T_{n}+s}:=x^{\prime}$ for all $0\leq s<t$.

  3. III.

Draw a bridge of $Z^{1}$ conditioned on its endpoints $Z^{1}_{T_{n}}$ and
$Z^{1}_{T_{n+1}}$, so that for every $0\leq s<t$ we obtain

$Z^{1}_{T_{n}+s}\sim\frac{p_{s}(x,y)p_{t-s}(y,x^{\prime})}{p_{t}(x,x^{\prime})}\mathbbm{1}_{\\{p_{t}(x,x^{\prime})>0\\}}\Lambda(dy)$
(11)

Let $Z^{1}_{T_{n}+s}:=x_{0}$ for some fixed $x_{0}\in E$ on
${\\{p_{t}(x,x^{\prime})=0\\}}$. Moreover, given $Z^{1}_{T_{n}+s}=y$ on
$s+u<t$ we have that

$Z^{1}_{T_{n}+s+u}\sim\frac{p_{u}(y,y^{\prime})p_{t-s-u}(y^{\prime},x^{\prime})}{p_{t-s}(y,x^{\prime})}\mathbbm{1}_{\\{p_{t-s}(y,x^{\prime})>0\\}}\Lambda(dy)$
(12)

Again, on ${\\{p_{t-s}(y,x^{\prime})=0\\}}$, let $Z^{1}_{T_{n}+s}=x_{0}.$

  4. IV.

At jump time $T_{n+1}$ we have $Z^{1}_{T_{n+1}}:=Z^{3}_{T_{n}}=x^{\prime}$.
Draw $Z^{2}_{T_{n+1}}$ independently of $Z_{s},s<T_{n+1}$, uniformly on the
unit interval. Given $\\{Z^{2}_{T_{n+1}}=u^{\prime}\\}$, generate

$Z^{3}_{T_{n+1}}\sim K((x^{\prime},u^{\prime}),dx^{\prime\prime}).$

Note that in the construction of $Z$ the inter-sampling times
$(\sigma_{n})_{n\geq 1}$ are drawn according to (10), their conditional
distribution given the starting and endpoint of the sampled chain. Equation
(11) and (12), describe the distributions of points in a bridge of the process
$X$. The first coordinate of $Z$ consists of bridges drawn according to the
law of the original process $X$, between realisations of the resolvent chain.
The results of Löcherbach and Loukianova [52; 53] that we work with are given
in the following propositions. Firstly, we see that the first coordinate of
$Z$ has the desired distribution.

######  Proposition 2.1.

The constructed process $Z$ is a Markov process with respect to its natural
filtration $\mathbb{F}$. Moreover, the first coordinate $Z^{1}$ is equal in
law to our process $X$, namely,

$\mathcal{L}((X_{t})_{t\geq 0}|X_{0}=x)=\mathcal{L}((Z^{1}_{t})_{t\geq
0}|Z^{1}_{0}=x).$

Moreover, $(T_{n}-T_{n-1})_{n\geq 1}$ are i.i.d exponential random variables
and are independent of $Z^{1}$; therefore, we also have that

$\mathcal{L}((X_{T_{n}})_{n\geq 0}|X_{0}=x)=\mathcal{L}((Z^{1}_{T_{n}})_{n\geq
0}|Z^{1}_{0}=x).$

Moreover, the process $X$ is embedded in a richer process $Z$, which admits a
recurrent atom $A:=C\times[0,\alpha]\times E$ in the sense of the following
proposition.

######  Proposition 2.2.

Let $(S_{n},R_{n})$ be a sequence of stopping times defined as
$S_{0}=R_{0}:=0$ and

$S_{n+1}:=\inf\\{T_{m}>R_{n}:Z_{T_{m}}\in A\\}\ \ \textrm{and}\ \
R_{n+1}:=\inf\\{T_{m}:T_{m}>S_{n+1}\\}.$

Then $Z_{R_{n}^{+}}$ is independent of $\mathcal{F}_{R_{n-1}}$ for all $n\geq
1$ and $(Z_{R_{n}})_{n\geq 1}$ is an i.i.d sequence with

$Z_{R_{n}}\sim\nu(dx)\lambda(du)K((x,u),dx^{\prime})\ \ \textrm{for all}\
n\geq 1.$

The stopping times $\\{S_{n}\\}_{n}$ thus denote the hitting times of the
recurrent atom $A$ for the jump process $(Z_{T_{n}})_{n}$, and
$\\{R_{n}\\}_{n}$ denote the implied regeneration epochs of the process $Z.$
As a direct consequence, we obtain the following regenerative structure for
the original process.

######  Proposition 2.3.

Let $f$ be a measurable $\pi$-integrable function, then we can construct a
sequence of increasing stopping times $\\{R_{n}\\}_{n}$ with $R_{0}=0$ and

$\xi_{n}:=\int_{R_{n-1}}^{R_{n}}f(X_{s})\ ds,\ \ n\geq 1,$

such that the sequence $\\{\xi_{n}\\}_{n}$ is a stationary sequence under
$\mathbb{P}_{\nu}.$ Moreover, for $n\geq 2$, $\xi_{n}$ is independent of
$\mathcal{F}_{R_{n-2}}.$

The regenerative structure given in Proposition 2.3 was also noted by Sigman
[75]. They define a process $X$ to be _one-dependent regenerative_ if there
exists, on a possibly enlarged probability space, a sequence of randomised
stopping times $R_{n}$ with corresponding cycle lengths
$\rho_{n}=R_{n+1}-R_{n}$ such that $\\{\left(X_{R_{n}+t}\right)_{t\geq
0},(\rho_{n+k})_{k}\\}$ has the same distribution for each $n\geq 1$ and are
independent of $\\{(\rho_{n})_{n=1}^{k-1},(X_{t})_{t<R_{n-1}}\\}$ for $n\geq
2.$ Note that according to this definition the initial cycle is allowed to
have a different distribution. Löcherbach and Loukianova [52] give a
constructive approach towards this result, in which they explicitly define the
corresponding stopping times and the recurrent atom. By the implied
regenerative structure of $X$, we obtain the following characterisation of the
stationary measure.

######  Proposition 2.4 (Sigman [75; Theorem 2]).

Let $X$ be a positive recurrent one-dependent regenerative process, then we
can characterise its stationary measure as follows

$\pi(A)=\dfrac{1}{\varrho}\
\mathbb{E}_{\nu}\int_{0}^{R_{1}}\mathbbm{1}_{\\{X_{s}\in A\\}}\ ds,$ (13)

where $\varrho$ is defined as $\mathbb{E}_{\nu}R_{1}.$ Moreover, we have the
following erdogic law of large numbers

$\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(X_{s})ds=\dfrac{1}{\varrho}\
\mathbb{E}_{\nu}\int_{0}^{R_{1}}f(X_{s})ds\ \ \textrm{a.s.},$ (14)

for all $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ with
$\pi(\absolutevalue{f})<\infty$. Note that the normalisation constant of $\pi$
in Proposition 2.4 is finite and non-zero due to the positive Harris
recurrence of the process.

######  Remark 2.5.

The framework of Löcherbach and Loukianova [52; 53] does not require
ergodicity. Moreover, it is important to note that contrary to the classically
regenerative case, Proposition 2.4 does not imply convergence in total
variation to the stationary measure. For a counterexample see Sigman [75;
Remark 3.2]. For our applications we will require ergodicity and hence we must
additionally impose this as stated in (5). ∎

Furthermore, we will rely on the following properties of the resolvent chain.
Granted that the process $X$ is aperiodic and positive Harris recurrent, then
also the resolvent $\bar{X}$ will inherit these properties, as seen in Meyn
and Tweedie [57; Propostion 5.4.5] and Tweedie [79; Theorem 3.1] respectively.
Moreover, by Down et al. [29; Theorem 5.3], exponential convergence to
stationarity is equivalent for $X$ and $\bar{X}$. The split chain of the
resolvent in turn obtains aperiodicity and positive Harris recurrence from
$\bar{X}$, as seen in for example Nummelin [61]. Following a co-de-
initialising argument of Roberts and Rosenthal [71], we see that the split
chain inherits the rate of convergence of the resolvent chain. To conclude, we
see that the split chain inherits aperiodicity, positive Harris recurrence,
and the rate of ergodicity from the process $X$.

##  3 Main Theorems

The most straightforward approach for obtaining a strong approximation result
for Markov processes would be through ergodicity requirements. Kuelbs and
Philipp [48] show that a multivariate strong invariance principle holds for
sums of random vectors satisfying a strong mixing condition; see also Theorem
6.1. This mixing condition can easily be satisfied by guaranteeing an
appropriate rate of ergodicity of the process. All proofs are provided in
Section 6.

######  Theorem 3.1.

Let $X$ be polynomially ergodic of order
$\beta\geq(1+\varepsilon)(1+2/\delta)$ for some $\varepsilon,\delta>0$. Then
for every initial distribution and for all $f:E\rightarrow\mathbb{R}^{p}$ with
$\pi(\norm{f}^{2+\delta})<\infty$, we can construct a process that is equal in
law to $X$ together with a standard $p$-dimensional Brownian motion
$W=\\{W(t)\\}_{t\geq 0}$ on some probability space such that

$\norm{\int_{0}^{T}f(X_{t})\
dt-T\pi(f)-\Sigma_{f}^{1/2}W(T)}=O(\psi_{T})\quad\textrm{a.s.}$ (15)

with

$\psi_{T}=T^{1/2-\min(\delta/(2\delta+4),\ \lambda)}\ \textrm{for some}\
\lambda\in(0,1/2),$

and positive semi-definite $p\times p$ covariance matrix $\Sigma_{f}$ given by

$\Sigma_{f}=\int_{0}^{\infty}\operatorname*{Cov_{\pi}}(f(X_{0}),f(X_{s}))\
ds+\int_{0}^{\infty}\operatorname*{Cov_{\pi}}(f(X_{s}),f(X_{0}))\ ds,$ (16)

with all entries converging absolutely and integration of matrices defined
element-wise.

######  Remark 3.2.

The asymptotic covariance matrix $\Sigma_{f}$ given in Theorem 3.1 cannot be
simplified. Only for the univariate case $(p=1)$ and for reversible processes
do we obtain that

$\Sigma_{f}=2\int_{0}^{\infty}\operatorname*{Cov_{\pi}}(f(X_{0}),f(X_{s}))\
ds.$ (17)

As a result of the reversibility, the cross-covariance matrices in (16) will
be symmetric and thus the asymptotic covariance can be expressed as (17). ∎

The rate $\psi_{T}$ appearing in Theorem 3.1 will depend on the dependence and
moment structure of the considered process. For processes admitting higher
order moments and having faster decaying levels of dependence the
approximation bound $\psi_{T}$ will tend to infinity at a slower rate. This
can be interpreted as the magnitude of the difference between the centred
additive functional of the process and the approximating Brownian motion being
smaller. Although result (15) has useful applications for arbitrary
$\lambda\in(0,1/2)$, many refined limit theorems require an explicit remainder
term, where more insight is given regarding the impact of the moment and
dependence structure on the approximation error. In order to derive a more
refined strong invariance principle we will make us of splitting arguments.
Following the continuous time Nummelin splitting technique, as introduced by
Löcherbach and Loukianova [52], it follows that the process can be embedded in
a richer process, which admits a recurrent atom. Hence the process can be
redefined such that it can be split in identically distributed blocks of
strongly mixing random variables. Therefore we can utilise the approximation
results for weakly $m$-dependent sequences of Berkes et al. [4] to obtain a
strong invariance principle; see also Theorem 6.4.

######  Theorem 3.3.

Let $X=(X_{t})_{t\geq 0}$ be an aperiodic, positive Harris recurrent Markov
process for which Assumption 1 is satisfied. Let $f:E\rightarrow\mathbb{R},$
be a given $\pi$-integrable function. Define the sequence of random times
$\\{R_{n}\\}_{n=1}^{\infty}$ and $\\{\xi_{n}\\}_{n=1}^{\infty}$ as given in
Proposition 2.2 and 2.3. Moreover, assume that

$\mathbb{E}_{\nu}[R_{1}^{q}]<\infty\quad\textrm{for some}\ \ q\geq 1,$ (18)
$\mathbb{E}_{\nu}\left[\left(\int_{0}^{R_{1}}f(X_{s})ds\right)^{2+\delta}\right]<\infty\quad\textrm{for
some}\ \ \delta>0.$ (19)

Then for every initial distribution we can construct a process, on an enriched
probability space, that is equal in law to $X$ together with two standard
Brownian motions $W_{1}$ and $W_{2}$ such that

$\displaystyle\absolutevalue{\int_{0}^{T}f(X_{s})ds-T\pi(f)-W_{1}(s_{T}^{2})-W_{2}(t_{T}^{2})}=O(\psi_{T})\
\textrm{a.s.,}$ (20)

where $\\{\sigma^{2}_{T}\\}$ and $\\{\tau^{2}_{T}\\}$ are non-decreasing
sequences with $\sigma^{2}_{T}=\frac{\sigma^{2}_{\xi}}{\varrho}T+O(T/\log T)$,
$\tau_{T}^{2}=O(T/\log T)$.

as $T\rightarrow\infty$, and $\psi_{T},\pi(f),\varrho,\textrm{and
}\sigma_{\xi}$ are defined in equations (22) to (25).

In Theorem 3.3 we obtain an explicit approximation error. In alignment with
expectations, we see that faster convergence to the stationary measure and the
existence of higher order moments will result in an improved approximation
error. However, the required moment conditions for Theorem 3.3 stated in (18)
and (19) are impractical and would be burdensome, if not impossible, to verify
directly for most applications. For classically regenerative Markov chains
this problem also arises, see the analogous requirements of regenerative
simulation given in Mykland et al. [59] and the strong approximation result of
Csáki and Csörgő [16]. Hobert et al. [41] were the first to simplify moment
conditions of this form and give practical sufficient conditions for
regenerative simulation. More specifically, in their main result they show
that polynomial or geometric ergodicity and moment conditions with respect to
the stationary measure are sufficient to guarantee finiteness of the variance
of a cycle. This result was generalised to higher order cycle moments by Jones
et al. [44] and Bednorz and Łatuszyński [3]; hence simplifying the required
conditions of Csáki and Csörgő [16]. However, the aforementioned approaches
are all for Markov chains satisfying a one-step minorisation condition, i.e.,
for the classically regenerative setting. Since our setting involves a more
complicated reconstruction of the process of interest, the results do not
immediately carry over. In Theorem 19, we show that the cycle moment
conditions (18) and (19) required for Theorem 3.3 can also be guaranteed with
more easily verifiable ergodicity and moment conditions.

######  Theorem 3.4.

Let $X=(X_{t})_{t\geq 0}$ be an aperiodic, positive Harris recurrent Markov
process for which Assumption 1 is satisfied. Moreover, let $X$ be ergodic with
rate $\Psi$, which decays at least polynomially of order $\beta\geq q+1$, then

$\mathbb{E}_{\nu}[R_{1}^{q}]<\infty\quad\textrm{for all}\ \ q\geq 1.$

Moreover, for all measurable $f:E\rightarrow\mathbb{R}$ with
$\pi(\absolutevalue{f}^{p+\varepsilon})<\infty$ where $p\geq 1$ and
$q>(p+\varepsilon)/\varepsilon$, we have that

$\mathbb{E}_{\nu}\left[\left(\int_{0}^{R_{1}}f(X_{s})ds\right)^{p+\varepsilon}\
\right]<\infty.$

By combining Theorem 3.3 and Theorem 3.4 we obtain the desired strong
invariance principle.

######  Theorem 3.5.

Let $X=(X_{t})_{t\geq 0}$ be an aperiodic, positive Harris recurrent Markov
process for which Assumption 1 is satisfied. Moreover, let $X$ be ergodic with
rate $\Psi$, which decays at least polynomially of order $\beta\geq q+1$ with
$q>(2+\delta+\varepsilon)/\varepsilon$ for some $\varepsilon,\delta>0$. Then
for every initial distribution and for all measurable
$f:E\rightarrow\mathbb{R}$ with $\pi(|f|^{2+\delta+\epsilon})<\infty$ we can,
on an enriched probability space, define a process that is equal in law to $X$
and two standard Brownian motions $W_{1}$ and $W_{2}$ such that

$\absolutevalue{\int_{0}^{T}f(X_{s})ds-T\pi(f)-W_{1}(\sigma_{T}^{2})-W_{2}(\tau_{T}^{2})}=O(\psi_{T})\
\textrm{a.s.,}$ (21)

where $\\{\sigma^{2}_{T}\\}$ and $\\{\tau^{2}_{T}\\}$ are non-decreasing
sequences with $\sigma^{2}_{T}=\frac{\sigma^{2}_{\xi}}{\varrho}T+O(T/\log T)$,
$\tau_{T}^{2}=O(T/\log T)$, and

$\displaystyle\psi_{T}$ $\displaystyle=\max\left\\{T^{1/2q}\log
T,T^{1/(2+\delta)}\log^{2}(T)\right\\},$ (22) $\displaystyle\pi(f)$
$\displaystyle=\frac{1}{\varrho}\ \mathbb{E}_{\nu}\int_{0}^{R_{1}}f(X_{s})\
ds,$ (23) $\displaystyle\varrho$
$\displaystyle=\mathbb{E}_{\nu}[R_{1}],\textrm{ and}$ (24)
$\displaystyle\sigma_{\xi}$
$\displaystyle=\sqrt{\operatorname*{Var_{\nu}}(\xi_{1})+2\operatorname*{Cov_{\nu}}(\xi_{1},\xi_{2})}\
.$ (25)

###### Proof.

The assertion follows immediately from Theorem 3.3 and Theorem 3.4. ∎

The appearance of the second Brownian motion in Theorem 3.5 is inherited from
the strong invariance principle of Berkes et al. [4]. Although we obtain
different time perturbations of the Brownian motions, all desired properties
carry over. The second Brownian motion appearing in (21) is of a smaller
magnitude, and will therefore be asymptotically negligible in typical
applications. Furthermore, even though the two Brownian motions are not
independent, their correlation decays over time

$\displaystyle\operatorname*{Corr}\left(W_{1}(\sigma_{t}^{2}),W_{2}(\tau_{s}^{2})\right)\rightarrow
0,\quad\textrm{as}\ t,s\rightarrow\infty.$ (26)

Note that the nearly optimal convergence rate of Berkes et al. [4] does not
carry over. Instead we obtain an approximation error that cannot be improved
beyond $O(T^{1/4}\log(T))$. Obtaining an approximation error superior to
$O(T^{1/4}\log(T))$ remains an open problem for the class of processes we
consider. A possible approach for attaining a better convergence rate would be
to extend to results of Berkes et al. [4] to a multivariate setting and follow
the approach of Merlevède et al. [55].

######  Remark 3.6.

Note that if the considered process would be regenerative then the strong
invariance principle of Csáki and Csörgő [16] would also be applicable. ∎

######  Remark 3.7.

The aforementioned results require a certain degree of polynomial ergodicity,
but can mutatis mutandis be seen to hold assuming exponential ergodicity. ∎

##  4 Estimation of the asymptotic variance of Piecewise Deterministic Monte
Carlo samplers

Recently, there has been growing interest in Piecewise Deterministic Monte
Carlo (PDMC) algorithms, which are based on Piecewise Deterministic Markov
Processes (PDMPs). It is well known that non-reversibility (a property
possessed by PDMPs) can significantly improve performance of sampling methods,
in terms of both convergence rate to equilibrium and asymptotic variance, see
for example, the results of Hwang et al. [43] and Lelievre et al. [49]
regarding convergence to stationarity and Duncan et al. [30] and Rey-Bellet
and Spiliopoulos [69] regarding the asymptotic variance. Furthermore, PDMPs
have piecewise deterministic paths and can therefore be simulated without
discretisation error, in contrast to for example Langevin and Hamiltonian
dynamics. The primary sampling algorithms belonging to this class are the Zig-
Zag sampler and the Bouncy Particle sampler, introduced by Bierkens and
Roberts [7] and Bouchard-Côté et al. [10] respectively. Moreover, since these
processes maintain the correct target distribution if sub-sampling is
employed, they enjoy advantageous scaling properties to large datasets as
explained in Bierkens et al. [8].

The aim of such a sampling method is to explore some probability distribution
$\pi$ on $\mathfrak{X}\subseteq\mathbb{R}^{d}$, such that we can infer
important characteristics of this distribution. Typically, the objective is to
compute expectations with respect to this distribution, in other words, for
some $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$, we are interested interested
in $\pi(f)$, the expectation with respect to our target distribution. Existing
PDMC samplers consist of a position and a velocity component. We will consider
processes $Z=(Z_{t})_{t\geq 0}$ with $Z_{t}=(X_{t},V_{t})$, where $X_{t}$ and
$V_{t}$ denote the position and velocity component respectively. Our process
takes values in $E=\mathfrak{X}\times\mathcal{V}$, where $\mathfrak{X}$
denotes the state-space of the position component and $\mathcal{V}$ denotes
the space of attainable velocities. Piecewise Deterministic Markov processes
are characterised by their deterministic dynamics between random event times
along with a Markov kernel that describes the transitions at events. More
specifically, their deterministic dynamics are described by some ordinary
differential equation. Both the Zig-Zag process and the Bouncy Particle
sampler have piecewise linear trajectories characterised by

$\frac{dX_{t}}{dt}=V_{t}\ \quad\textrm{and}\ \quad\frac{dV_{t}}{dt}=0.$

Thus the velocity describes the rate of change of the position, and does not
change along the deterministic dynamics. Changes in the velocity occur
according to some inhomogeneous Poisson rate $\lambda(Z_{t})$. The Poisson
events consist of changes in the velocity component of our process. The
fundamental idea behind these sampling methods is to choose the event rate and
the changes in velocity such that the position component explores the state-
space according to the target distribution $\pi$. The event rate increases as
the position is moving towards regions of lower probability mass. Since the
exact specification of these processes are in principle not relevant for our
results, we refer to Bierkens et al. [8] and Bouchard-Côté et al. [10] for the
simulation scheme of the Zig-Zag process and Bouncy Particle Sampler
respectively. It can be shown that under very mild regularity conditions both
sampling processes admit a stationary distribution given by

$\displaystyle\mu(dx,dv)=\pi(dx)\nu(dv),$ (27)

where the target distribution $\pi$ is the marginal distribution of the
position component and $\nu$ is the marginal distribution of the velocity
component. Moreover, an ergodic law of large numbers holds

$\displaystyle\lim_{T\rightarrow\infty}\hat{\mu}_{T}(f):=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}f(X_{s},V_{s})\
ds=\int_{E}f(x,v)\mu(dx,dv)=:\mu(f),$

for all $\mu$-integrable $f$. Given the independence of position and velocity
at equilibrium, as seen in (27), the time average of the position component is
a natural estimator for the space average $\pi(f)$. For obtaining insight
regarding the accuracy of this estimator we require a central limit theorem to
hold.

######  Corollary 4.1.

Let $Z=(X_{t},V_{t})_{t\geq 0}$ be polynomially ergodic of order
$\beta\geq(1+\varepsilon)(1+2/\delta)$ for some $\varepsilon,\delta>0$ with
stationary measure $\mu$ satisfying (27). Then for all
$f:E\rightarrow\mathbb{R}^{p}$ with $\pi(\norm{f}^{2+\delta})<\infty$, a
central limit theorem (CLT) holds:

$\frac{1}{\sqrt{T}}\int_{0}^{T}(f(X_{s},V_{s})-\mu(f))\
ds\overset{d}{\longrightarrow}\mathcal{N}_{p}(0,\Sigma_{f}).$ (28)

Additionally, also a functional central limit theorem (FCLT) holds:

$\left(\frac{1}{\sqrt{n}}\int_{0}^{nt}(f(X_{s},V_{s})-\mu(f))\
ds\right)_{t\in[0,1]}\overset{d}{\longrightarrow}\Sigma_{f}^{1/2}W\
\textrm{as}\ n\rightarrow\infty,$ (29)

where

$\Sigma_{f}=\int_{0}^{\infty}\operatorname*{Cov_{\mu}}(f(X_{0},V_{0}),f(X_{s},V_{s}))\
ds+\int_{0}^{\infty}\operatorname*{Cov_{\mu}}(f(X_{s},V_{s}),f(X_{0},V_{0}))\
ds,$ (30)

$W=(W_{t})_{t\in[0,1]}$ denotes a standard $p$-dimensional Brownian motion,
and the weak convergence is with respect to the Skorohod topology on $D[0,1]$,
the space of real-valued càdlàg functions with domain $[0,1]$.

###### Proof.

By Philipp and Stout [66; Theorem 1.E], the FCLT immediately follows from the
strong invariance principle formulated in Theorem 3.1. Similarly, by Damerdji
[23; Proposition 2.1] the CLT follows. ∎

With regard to sufficient conditions for the required ergodicity, we remark
that in Bierkens et al. [9] exponential ergodicity of the Zig-Zag process is
shown for target distributions that have a non-degenerate local maximum and
appropriately decaying tails. In Deligiannidis et al. [27] and Durmus et al.
[31] conditions for exponential ergodicity of the Bouncy Particle Sampler are
given. Utilising hypocoercivity techniques, Andrieu et al. [1] establish
polynomial rates of convergence for PDMPs with heavy-tailed stationary
distributions.

In order to assess the accuracy of the ergodic average resulting from our
sampling method, we need to estimate the asymptotic variance appearing in
(30). Moreover, the asymptotic variance is also useful for determining the
efficiency of the sampling algorithm via measures such as the effective sample
size, see for example Gong and Flegal [38] or Vats et al. [81]. The estimation
of the asymptotic variance is also required for the implementation of stopping
rules, which consists of justifiable criteria for termination of the
simulation. In order to validate stopping rules that guarantee a desired level
of precision, Glynn and Whitt [36] show that the estimator of the asymptotic
covariance matrix must be strongly consistent. Damerdji [22; 23] gives a
general framework for the strong consistency of batch means, spectral variance
and standardised time series methods. The central assumption for their results
is the formulation of a strong invariance principle for the underlying
sampling process. We see that the interplay between the approximation error in
the strong invariance principle and the various tuning parameters of the
considered estimation methods plays a crucial role. Although consistency can
usually be attained under more lenient assumptions (such as for example only
an FCLT), these results provide insufficient guidance for the practical
implementation of these methods, since they provide no explicit choice for the
tuning parameters of the respective methods. The tuning parameters greatly
affect the quality of the estimation procedure, hence choosing these according
to more refined criteria usually results in better performance.

For simplicity, we will mainly consider the one-dimensional case, i.e. our
quantity of interest is given by $\pi(f)$, with
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ a given $\pi$-integrable function. Let
the simulation output, which in our case consists of the position component of
a PDMP, be given by $(X_{t})_{t\in[0,T]}$. Note that from Corollary 4.1 also a
(functional) central limit theorem follows for the position component $X$. We
are interested in estimating the asymptotic variance (30); which we will
denote with $\sigma^{2}_{f}$, when we are not considering the multivariate
setting.

###  4.1 Batch means

The batch means method divides the obtained sample trajectory of our process
into non-overlapping parts. The sample variance of the means of the obtained
batches gives rise to a natural estimator for the asymptotic variance. More
specifically, we divide our simulation output in $k_{T}$ batches of length
$\ell_{T}$ such that $k_{T}=\lfloor T/\ell_{T}\rfloor$. We proceed by
computing the sample average of each obtained batch;

$\bar{Z}_{i}(\ell_{T}):=\frac{1}{\ell_{T}}\int_{(i-1)\ell_{T}}^{i\ell_{T}}f(X_{s})ds,\quad
i=1,\dots,k_{T}.$

If a functional central limit theorem holds for our process, it follows that
the computed means $\bar{Z}_{i}(\ell_{T})$ are asymptotically independent and
identically distributed for all fixed amount of batches. Hence, we can
heuristically reason that the sample variance of
$(\bar{Z}_{i}(\ell_{T})_{i=1}^{k_{T}}$ will be close to
$\operatorname*{Var}(\bar{Z}_{i}(\ell_{T})).$ Moreover, since each
$\bar{Z}_{i}(\ell_{T})$ is also an empirical mean, it is reasonable to expect
their variance to be approximately $\sigma_{f}^{2}/\ell_{T}.$ The batch means
estimator of the asymptotic variance is defined by correcting the sample
variance of the batch means $(\bar{Z}_{i}(\ell_{T})_{i=1}^{k_{T}}$ by a factor
$\ell_{T}$, namely

$\hat{\sigma}_{T}^{2}=\frac{\ell_{T}}{k_{T}-1}\sum_{i=1}^{k_{T}}\left(\bar{Z}_{i}(\ell_{T})-\frac{1}{k_{T}}\sum_{i=1}^{k_{T}}\bar{Z}_{i}(\ell_{T})\right)^{2}.$
(31)

Following the framework of Damerdji [23], we impose the following conditions
on the amount of batches and their length.

######  Assumption 2.

Let the amount of batches $k_{T}$ and their lengths $\ell_{T}$ be such that

  1. 1.

$k_{T}\rightarrow\infty$ , $\ell_{T}\rightarrow\infty$, and
$\ell_{T}/T\rightarrow 0$ as $T\rightarrow\infty$,

  2. 2.

$\ell_{T}$ and $\ell_{T}/T$ are both monotonically increasing,

  3. 3.

there exists a constant $c\geq 1$ such that
$\sum_{n=1}^{\infty}k_{n}^{-c}<\infty$.

The first requirement of Assumption 2 is a necessary condition for consistency
as seen from the results of Glynn and Whitt [35]. The second requirement is
solely for technical reasons and the third requirement ensures that the amount
of the batches grows fast enough; if we choose $\ell_{T}=T^{\alpha}$ the
requirement holds for all $\alpha\in(0,1)$, since we can choose
$c>1/(1-\alpha)$.

######  Theorem 4.2.

Let $Z$ be polynomially ergodic of order
$\beta\geq(1+\varepsilon)(1+2/\delta)$ for some $\varepsilon,\delta>0$ with
stationary measure $\mu$ satisfying (27) and let
$\pi(\absolutevalue{f}^{2+\delta})<\infty$. Assume that Assumption 2 holds and
that

$\frac{\psi_{T}^{2}}{\ell_{T}}\log(T)\rightarrow 0,\ \textrm{as}\
T\rightarrow\infty,$ (32)

where

$\psi_{T}=T^{1/2-\lambda}\ \textrm{with}\ \lambda=\min(\delta/(2\delta+4),\
\lambda^{\prime})\ \textrm{for some}\ \lambda^{\prime}\in(0,1/2),$ (33)

then for every initial distribution
$\hat{\sigma}^{2}_{T}\rightarrow\sigma^{2}_{f}$ as $T\rightarrow\infty$ with
probability 1.

###### Proof.

The result follows from Theorem 3.1 and Damerdji [23; Theorem 3.3]. ∎

We see from the required assumption (32) that a larger approximation error in
the strong invariance principle, which corresponds to higher orders of
dependence, results in a larger required batch size $\ell_{T}$. This is in
agreement with the idea behind batching methods; every batch should give a
proper representation of the dependence structure of the process. Otherwise a
structural bias will be introduced in the estimation procedure. On the other
hand, choosing the batch size larger than necessary, will result in a lower
amount of batches $k_{T}$ leading to a higher variance for the estimator.
Strong approximations can also be used to characterise the mean squared error
and obtain a central limit theorem for the batch means estimator.

######  Theorem 4.3.

Let $Z$ be polynomially ergodic of order
$\beta\geq(1+\varepsilon)(1+2/\delta)$ for some $\varepsilon,\delta>0$ with
stationary measure $\mu$ satisfying (27) and let
$\pi(\absolutevalue{f}^{2+\delta})<\infty$. Let the initial distribution be
given by $\mu$ and assume that Assumption 2 holds and
$\mathbb{E}_{\mu}C^{2}<\infty$, where $C$ is defined in (54). Then we have
that

$\mathbb{E}_{\mu}\absolutevalue{\hat{\sigma}^{2}_{T}-\sigma_{f}^{2}}^{2}=2\sigma_{f}^{4}\frac{\ell_{T}}{T}+O\left(\frac{\psi_{T}}{\sqrt{T}}\log^{\frac{1}{2}}T\right)+O\left(\ell_{T}^{-1}\psi_{T}^{2}\log
T\right),$ (34)

with $\psi_{T}$ given in (33). Moreover, if $\ell_{T}^{-1}\psi_{T}(T\log
T)^{1/2}\rightarrow 0$ as $T\rightarrow\infty$, then we obtain a CLT for the
batch means estimator

$\sqrt{k_{T}}(\hat{\sigma}^{2}_{T}-\sigma^{2}_{f})\xrightarrow{d}\mathcal{N}(0,2\sigma_{f}^{4})\
\textrm{as}\ T\rightarrow\infty.$ (35)

###### Proof.

By the imposed conditions of the process, the strong invariance principle
formulated in Theorem 3.1 holds. The first claim then follows by Damerdji [24;
Theorem 1 and Lemma 3] and the second by Damerdji [24; Proposition 2]. ∎

The first and second term in (34) describe the variance, whereas the third
term represents the bias. Note that the second term does not depend on
$\ell_{T}$ and tends to zero. The obtained bounds for the variance are sharp,
whereas, the bounds for bias have room for improvement.

####  4.1.1 Batch size selection for PDMC

Glynn and Whitt [35] show that there exists no consistent estimator of
$\sigma_{f}^{2}$ with a fixed amount of batches. Hence the amount of batches
should explicitly depend on the length of the simulation $T$. For the standard
choice $\ell_{T}=T^{\alpha}$ we see that for $\alpha>1-2\lambda^{\prime}$ we
obtain both strong consistency and $L^{2}$-convergence of the batch means
estimator. Theorem 4.3 suggests that $\alpha^{*}=1-\lambda$ would be optimal
in the mean squared error sense, for $T$ of increasing magnitude, where
$\lambda$ is defined in (33). The well-known results of Chien et al. [14],
Goldsman et al. [37], and Song and Schmeiser [76] obtain a bound for the bias
of order $O(\ell_{T}^{-1})$, which implies an optimal (in the MSE sense) batch
size of ${\ell}^{\diamond}_{T}\asymp{T}^{1/3}.$ However, the aforementioned
results require the sampling process to be stationary, uniformly ergodic, and
satisfy moment condition $\pi(f^{12})<\infty$. Obtaining the bias term of
order $O(\ell_{T}^{-1})$ for batch means under milder conditions remains an
unaddressed problem. Theorem 4.2 and 4.3 only require a strong invariance
principle, which we have shown holds under polynomial ergodicity; a very
reasonable assumption for simulation output. Moreover, these results do not
require stationarity and thus hold for every initial distribution. Theorem 4.3
imposes more demanding conditions on $\ell_{T}$ than aforementioned
frameworks, however, it is quite reasonable to let the batch size depend on
the dependence structure of the process through $\psi_{T}$, instead of only on
through the constant $\int s\gamma(s)ds$, as is the case in the aforementioned
results. Moreover, in practice the performance of batch means methods with
batch size ${\ell}^{\diamond}_{T}$ are often found to be sub-optimal whereas
larger batch sizes see better finite sample performance, as noted by for
example Flegal et al. [33].

An alternative approach for determining the optimal batch size was given by
Chien [15], who obtains an optimal batch size of $\tilde{\ell}_{T}\asymp
T^{1/2}$ by minimising the distance between the cumulants of the studentised
ergodic average and a standard Gaussian, which suggest that the resulting
confidence intervals enjoy improved finite-sample properties.

####  4.1.2 Asymptotic normality of the batch means estimator

We see that given polynomial ergodicity, also the central limit theorem for
the batch-means estimator carries over to the PDMC setting. The results of
Sherman and Goldsman [73] require uniform ergodicity and the moment condition
$\pi(f^{12})<\infty$ in order to obtain asymptotic normality of the batch
means estimator. Since uniform ergodicity is not attainable for most practical
problems, less stringent conditions on the rate of ergodicity are desired.
Theorem 4.3 places more restricting conditions on the batch size and excludes
the choice $\ell^{\diamond}_{T}\asymp T^{1/3}.$ Chakraborty et al. [13] obtain
a CLT for the batch means estimator assuming reversibility, stationarity,
geometric ergodicity and moment condition $\pi(f^{8})<\infty.$ Moreover, the
required batch size must be such that $k_{T}=o(\ell_{T}^{2})$. Hence their
result is also unable to guarantee asymptotic normality for batch size
$\ell^{\diamond}_{T}.$ We see that Theorem 4.3 gives more practical conditions
for guaranteeing asymptotic normality of the batch-means estimator, in
particular, the results are applicable to non-reversible processes.

####  4.1.3 Spectral variance and overlapping batch means estimators for the
PDMC standard error

Analogous to the batch means method, given the strong invariance principle
formulated in Theorem 3.1, many results for other estimators of the asymptotic
variance also carry over. Flegal et al. [33] gives more convenient
alternatives for some of the requirements of the framework given in Damerdji
[22]. The results of Flegal et al. [33] regarding spectral variance and
overlapping batch means estimators for MCMC output are thus also applicable
for PDMC, with minor adjustments to their assumptions. Note that the assumed
minorisation condition and geometric ergodicity of the Markov chain in Flegal
et al. [33] are only imposed such that the strong invariance principle of
Csáki and Csörgő [16] holds. Although implementation of spectral variance
estimators for continuous-time output might be impractical, these estimators
are still of theoretical interest. Numerous estimation methods, such as
overlapping batch means and certain standardised time series methods, with
feasible implementation for PDMC output, can be show to be (asymptotically)
equivalent to spectral estimators. Furthermore, we expect that results of Vats
et al. [81] regarding the multivariate extensions of batch means also carry
over. Similarly, we expect the results of Vats et al. [80] and Liu et al. [51]
regarding spectral variance and generalised overlapping batch means estimators
respectively to remain valid in the continuous-time setting. Hence also the
implications for the optimal values of the tuning parameters of these
estimation methods for the asymptotic variance remains valid. Lastly, note
that our results hold for all sampling algorithms that produce continuous-time
output, and are not restricted to the PDMP setting.

##  5 Increments of Additive Functionals of Ergodic Markov processes

Strong approximation results enable various asymptotic properties of Brownian
motion to carry over to other stochastic processes. In this section, we show
that the strong invariance principle given in Theorem 3.5 can be used to show
that the increments of additive functionals of Markov processes are of the
same magnitude as Brownian increments, provided we have sufficient decay of
the approximation error. The following theorem describes the magnitude of the
fluctuations of Brownian increments over subintervals of length $a_{T}.$

######  Theorem 5.1 (Csörgö and Révész [18; Theorem 1]).

Let $W=(W_{t})_{t\geq 0}$ denote a Brownian motion, and let $a_{T}$ be a
positive non-decreasing function of $T$ such that $0<a_{T}\leq T$ and
$T/a_{T}$ is non-decreasing. Then

$\limsup_{T\rightarrow\infty}\sup_{0\leq t\leq T-a_{T}}\sup_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{W_{t+u}-W_{t}}=1\quad\textrm{a.s.},$ (36)

where

$\beta_{T}=\left(2a_{T}\left[\log\frac{T}{a_{T}}+\log\log
T\right]\right)^{-1/2}.$

Taking $a_{T}=T$ gives the law of iterated logarithm, and for $a_{T}=c\log T$
with $c>0$, the Erdös-Rényi law of large numbers for Brownian motion is
obtained. This fluctuation result has been extended to other processes such as
integrated Brownian motion, fractional Brownian motion, and non-stationary
Gaussian processes, see Li [50], El-Nouty [32] and Ortega [64] respectively.
While these fluctuation results are of independent interest, they are also
used as building blocks in applications, such as proving convergence
properties of kernel density estimators, see for example Révész [67] and
Deheuvels [26]. These fluctuation results are also used for proving almost
sure convergence of various estimators of the asymptotic variance in
simulation output settings, see the references given in Section 4.1. By the
Komlós-Major-Tusnády approximation the fluctuation result immediately carries
over for i.i.d. sequences satisfying appropriate moment conditions, as seen in
Csörgö and Révész [19; Theorem 3.1.1 and 3.2.1].

In order to describe the fluctuations of additive functionals over an interval
of a specified length $a_{T}$, we require an explicit remainder term for the
Brownian approximation, as given in Theorem 3.5. However, due to the
appearance of the second Brownian motion in this invariance principle and the
perturbed time sequences, it is not immediate that the Brownian fluctuation
result carries over. Berkes et al. [4] show that the magnitude of the
increments of partial sums of weakly $m$-dependent sequences are indeed given
by Theorem 5.1, due to the smaller scaling of the second Brownian motion.
However, in our case the perturbed time sequences are random since they depend
on the amount of one-dependent regenerative cycles of the process, hence the
desired result does not follow directly by Berkes et al. [4; Theorem 4].

######  Theorem 5.2.

Let $X=(X_{t})_{t\geq 0}$ be an aperiodic, positive Harris recurrent Markov
process for which Assumption 1 is satisfied. Moreover, let $X$ be ergodic with
rate $\Psi$, which decays at least polynomially of order $\beta\geq q+1$ with
$q>(2+\delta+\varepsilon)/\varepsilon$ for some $\varepsilon,\delta>0$.
Consider a function $f:E\rightarrow\mathbb{R}$ with $\pi(f)=0$ and
$\pi(|f|^{2+\delta+\epsilon})<\infty$. Let $a_{T}$ be a given positive non-
decreasing function of $T$ such that

  1. 1.

$\displaystyle 0<a_{T}\leq T$,

  2. 2.

$\displaystyle T/a_{T}$ is non-decreasing,

  3. 3.

$\displaystyle a_{T}$ is regularly varying at $\infty$ with index
$\zeta\in(0,1].$

Suppose that $\beta_{T}\psi_{T}=o(1)$, where

$\beta_{T}=\left(2a_{T}\left[\log\frac{T}{a_{T}}+\log\log
T\right]\right)^{-1/2},$

and

$\psi_{T}=\max\left\\{T^{1/2q}\log T,T^{1/(2+\delta)}\log^{2}(T)\right\\}.$

Then we have that

$\limsup_{T\rightarrow\infty}\sup_{0\leq t\leq T-a_{T}}\sup_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{\int_{t}^{t+u}f(X_{s})ds}\leq\frac{\sigma^{2}_{\xi}}{\varrho}\quad\textrm{a.s.}$
(37)

As noted by Berkes et al. [4], the spit invariance principle also implies the
distributional version of Theorem 5.2; with similar adaptations to their
argument this would also hold in our case. Since the approximation error
$\psi_{T}$ of Theorem 3.5 cannot be guaranteed to be smaller than
$O(T^{1/4}\log(T))$, the fluctuation result given in Theorem 5.2 cannot
describe the magnitude of increments over slowly growing time intervals
$a_{T}$.

###  5.1 Application to diffusion processes

Diffusions are an important class of processes for which the strong
approximation given in Theorem 3.5 and the related fluctuation result given in
Theorem 5.2 are applicable. Let $X=(X_{t})_{t\geq 0}$ denote a one-dimensional
diffusion process that is defined as the solution of the following time-
homogeneous stochastic differential equation (SDE)

$\vspace{-0.15cm}\left\\{\begin{array}[]{ll}dX_{t}=b(X_{t})dt+\sigma(X_{t})dW_{t}\\\
X_{0}\sim\mu,\end{array}\right.$ (38)

where $\mu$ is the initial distribution of the process,
$\mathfrak{X}\subseteq\mathbb{R}$ denotes the state-space,
$b:\mathfrak{X}\xrightarrow{}\mathbb{R}$ and
$\sigma:\mathfrak{X}\xrightarrow{}\mathbb{R}$ denote the drift and volatility
function respectively, and the process $W$ is a Brownian motion. We assume
that all required regularity conditions hold such that the existence and
uniqueness of a strong solution of the SDE is guaranteed. For example, we can
impose Lipschitz conditions on the drift and volatility of the SDE. For a more
detailed explanation, we refer to Rogers and Williams [72].

For diffusion processes to admit the desired ergodic properties we must impose
additional regularity conditions. The scale function of a one-dimensional
diffusion process is given by

$s(u)=\int_{x}^{u}\exp\left[-2\int_{x}^{z}\frac{b(y)}{\sigma^{2}(y)}dy\right]dz\
\textrm{and must satisfy}\lim_{u\rightarrow\pm\infty}s(u)=\pm\infty.$ (39)

If condition (39) holds it follows that the diffusion is recurrent, the time
for the process to return to any bounded subset of its state space is a.s.
finite. The speed density of the diffusion process
$m:\mathfrak{X}\rightarrow\mathbb{R}^{+}$, given by
$m(u)=\left(s(u)\sigma^{2}(u)\right)^{-1},$ must be Lebesque integrable for
the diffusion to be positive Harris recurrent. For higher dimensional
diffusion process Bhattacharya [6] gives conditions that guarantee positive
Harris recurrence. In order for the obtained strong invariance principles
given in Theorem 3.1 and Theorem 3.5 to hold, we require polynomial or
exponential convergence to stationarity. These ergodicity requirements are
usually established through Foster–Lyapunov drift conditions; see Down et al.
[29] and Fort and Roberts [34] for exponential and polynomial ergodicity
respectively. These results have been applied to various classes of diffusion
processes, see for example Cattiaux et al. [12; Theorem 8.3 and 8.4] and
Stramer and Tweedie [77; Theorem 3.1 and 4.1].

In order for the strong approximation result in Theorem 3.5 and the related
fluctuation result of Theorem 5.2 to hold, the Nummelin splitting scheme of
Löcherbach and Loukianova [52] must be applicable. Therefore we must impose
regularity conditions such that Assumption 1 is satisfied, i.e., the
transition semigroup of the diffusion must be Feller and admit densities with
respect to some dominating measure. Under appropriate growth and continuity
conditions on the drift and volatility, diffusion processes are Feller, see
for example Williams [82; Theorem 2.2]. Moreover, if the volatility function
$\sigma$ is strictly positive (positive-definite in the multivariate case),
the diffusion is elliptic and therefore admits transition densities; Stroock
and Varadhan [78; Theorem 3.2.1]. Hence, Assumption 1 is satisfied.
Alternatively, for multivariate diffusions, we can impose the parabolic
Hörmander condition which ensures that the propagation of the noise through
the different coordinates is sufficient, such that the transition density
exists, see for example Rogers and Williams [72; Theorem 38.16].

####  5.1.1 Discussion and suggestions for further research

We see that Theorem 3.1 and Theorem 3.5 are applicable to a broad class of
diffusions and extend the current results on strong approximations for
diffusion processes. Heunis [40] and Mihalache [58] obtain strong invariance
principles for diffusions and a complementary fluctuation result and change
point test respectively. The results of Mihalache [58] yield an explicit
approximation error comparable to that of Theorem 3.5, but are only applicable
to stochastic integrals with respect to Brownian, i.e., diffusion processes
with no drift. The results of Heunis [40] give an implicit approximation error
and hold for singular diffusions. The strong invariance principle of Heunis
[40] is not covered by our results since singular diffusions generally do not
satisfy the mixing properties required for our framework.

The obtained strong invariance principles offer various applications for
diffusion processes, see for example Csörgö and Hall [20] and their given
references.

Following the approach of Berkes et al. [4; Proposition 2], Theorem 3.5 can be
used to obtain a change-point test for diffusions. If the diffusion process we
consider has a drift that enforces mean-reversion, we could construct a test
for the existence of a deterministic linear trend over specified time periods.
This approach would require continuous-time output of a diffusion process, and
is therefore more of theoretical interest. However, it is plausible that the
asymptotic behaviour of the change-point test should carry over to the high-
frequency setting, where the diffusion is observed discretely and it is
assumed the inter-observation times tend to zero.

##  6 Proofs

###  6.1 Theorem 3.1

Kuelbs and Philipp [48] give a strong invariance principle for mixing random
variables. In order to state their result, we first briefly introduce mixing
coefficients. Let $\mathscr{A}$ and $\mathscr{B}$ denote two sub
$\sigma$-algebras of our probability space. The $\alpha$-mixing coefficients
of two $\sigma$-algebras quantify their dependence as follows

$\alpha(\mathscr{A},\mathscr{B})=\sup\\{\Pr(F\cap
G)-\Pr(F)\Pr(G):F\in\mathscr{A},\ G\in\mathscr{B}\\}.$

The mixing coefficients of a stochastic process $X$, endowed with its natural
filtration, are defined as
$\alpha_{X}(s):=\sup_{t}\alpha(\mathcal{F}_{-\infty}^{t},\mathcal{F}_{t+s}^{\infty})\
\textrm{for }s>0,$ with $\mathcal{F}_{-\infty}^{t}=\sigma(X_{u}:u\leq t)$ and
$\mathcal{F}_{t+s}^{\infty}=\sigma(X_{u}:u\geq t+s)$. The mixing coefficients
of a process measure the dependence between events in terms of units of time
that they are apart. For a stationary Markov process the mixing coefficients
simplify as $\alpha(s)=\alpha(\sigma(X_{0}),\sigma(X_{s}))$, as shown in for
example Bradley [11; page 118].

######  Theorem 6.1 (Kuelbs and Philipp [48; Theorem 4]).

Let $\xi=(\xi_{k})_{k=1}^{\infty}$ be a weakly stationary sequence taking
values in $\mathbb{R}^{p}$ with mean zero and
$\sup_{k}\mathbb{E}\norm{\xi_{k}}^{2+\delta}\leq 1$, for some
$\delta\in(0,1]$. Moreover, let $\alpha_{\xi}$ the $\alpha$-mixing
coefficients of $\xi$ decay polynomially with rate
$n^{-(1+\varepsilon)(1+2/\delta)}$ for some $\varepsilon>0$. Then we can
redefine $\xi$ on a new probability space on which we can also construct a
$p$-dimensional Brownian motion $W$ with covariance matrix $\Sigma_{\xi}$,
with absolutely converging entries

$(\Sigma_{\xi})_{ij}=\mathbb{E}[\xi_{i1}\xi_{j1}]+\sum_{k=2}^{\infty}\mathbb{E}[\xi_{i1}\xi_{jk}]+\sum_{k=2}^{\infty}\mathbb{E}[\xi_{ik}\xi_{j1}],\
\textrm{for}\ 1\leq i,j\leq p,$

such that

$\norm{\sum_{k=1}^{n}\xi_{k}-W(n)}=O(n^{1/2-\lambda_{\xi}})\ \textrm{a.s.}$

for some $\lambda_{\xi}\in(0,1/2)$ depending only on $\varepsilon,\delta$ and
$p$.

The following lemmata are useful in the proof of Theorem 3.1.

######  Lemma 6.2 (Douc et al. [28; Theorem F.3.3]).

Let $X$ be a stationary ergodic Markov process with rate of convergence to
stationarity given by $\psi$, then $\alpha_{X}(s)$, the $\alpha$-mixing
coefficients of the process $X$, decay according to $\psi$, i.e., for all
$s\geq 0$ we have that

$\alpha_{X}(s)\leq\pi(V)\Psi(s),$

where $\Psi$ and $V$ are as stated in (5).

######  Lemma 6.3 (Davydov [25] and Rio [70]).

Let $(\Omega,\mathscr{F},\Pr)$ be a probability space and $\mathscr{A}$ and
$\mathscr{B}$ be two sub $\sigma$-algebras and consider random variables $X$
and $Y$ that are measurable with respect to these $\sigma$-algebras
respectively. Moreover, assume that $X\in L^{p}(\Pr)$ and $Y\in L^{q}(\Pr)$,
for some $p,q\geq 1$. Then we can bound their covariance in terms of the
$\alpha$-mixing coefficients as follows

$\lvert\operatorname*{Cov}(X,Y)\rvert\leq
8\alpha\left(\mathscr{A},\mathscr{B}\right)^{1/r}\norm{X}_{p}\norm{Y}_{q},\
\textrm{where}\ p,q,r\in[1,\infty]\ \textrm{and}\
\frac{1}{p}+\frac{1}{q}+\frac{1}{r}=1.$

####  6.1.1 Proof of Theorem 3.1

###### Proof.

For notational convenience introduce $Y=(Y_{t})_{t\geq 0}$, where
$Y_{t}:=f(X_{t})-\pi(f)$ for $t\geq 0$, and $\xi=(\xi_{k})_{k=1}^{n}$, with
$n:=n_{T}:=\lfloor T\rfloor$ and $\xi_{k}:=\int_{k-1}^{k}Y_{t}dt$ for
$k=1,\cdots,n$. Note that $Y_{t}$ is a $p$-dimensional vector, i.e.,
$Y_{t}=(Y_{1t},\cdots,Y_{pt})^{\top}$ and therefore also each $\xi_{k}$ is a
$p$-dimensional vector, $\xi_{k}=(\xi_{1k},\cdots,\xi_{pk})^{\top}.$
Furthermore, by definition $n$ is a function of the sample size $T$, however,
for notational convenience we suppress this. Since we are in the setting of
Lemma 6.2, $X$ has polynomially decaying $\alpha$-mixing coefficients, which
we will denote with $(\alpha_{X}(s))_{s\geq 0}$. Consequently, we have that
$Y$ and $\xi$ are both stationary and strongly mixing processes with
polynomially decaying $\alpha$-mixing coefficients $(\alpha_{Y}(s))_{s\geq 0}$
and $(\alpha_{\xi}(h))_{h\in\mathbb{N}}$ respectively. This can easily be seen
by observing that
$\sigma\left(f(X_{t})\right)\subseteq\sigma\left(X_{t}\right)$ and
$\sigma\left(\xi_{k}\right)\subseteq\sigma\left(X_{s}:s\leq k\right)$. In
order to show that a strong invariance holds for $Y$, we will show that it
holds for $\xi$ and determine the growth rate of the corresponding remainder
terms. Moment conditions for $\xi$ are directly inherited by the assumed
moment conditions for $X$. By an application of Jensen’s inequality we see
that

$\displaystyle\pi\left(\norm{\xi_{k}}^{2+\delta}\right)$
$\displaystyle=\mathbb{E}_{\pi}\norm{\int_{k-1}^{k}Y_{s}ds}^{2+\delta}\leq\
\mathbb{E}_{\pi}\int_{k-1}^{k}\norm{Y_{s}}^{2+\delta}ds=\pi\left(\norm{f-\pi(f)}^{2+\delta}\right)<\infty.$

Therefore, by Theorem 6.1, we can redefine $\xi$ on a new probability space on
which we can also construct a $p$-dimensional Brownian motion $W$ with
covariance matrix $\Sigma_{\xi}$, with absolutely converging entries

$(\Sigma_{\xi})_{ij}=\mathbb{E}[\xi_{i1}\xi_{j1}]+\sum_{k=2}^{\infty}\mathbb{E}[\xi_{i1}\xi_{jk}]+\sum_{k=2}^{\infty}\mathbb{E}[\xi_{ik}\xi_{j1}],\
\textrm{for}\ 1\leq i,j\leq p,$

such that

$\norm{\sum_{k=1}^{n}\xi_{k}-W(n)}=O(n^{1/2-\lambda_{\xi}})\ \textrm{a.s.}$

for some $\lambda_{\xi}\in(0,1/2)$ depending only on $\varepsilon,\delta$ and
$p$.

The claim follows if we show that

$\displaystyle\norm{\sum_{k=1}^{n}\xi_{k}-\int_{0}^{T}Y_{t}dt}$
$\displaystyle=O(T^{1/(2+\delta)})\ \textrm{a.s. for}\ T\rightarrow\infty,$
(40) $\displaystyle\norm{W_{T}-W_{n}}$ $\displaystyle=o(T^{1/(2+\delta)})\
\textrm{a.s. for}\ T\rightarrow\infty,\ \textrm{and}$ (41)
$\displaystyle\Sigma_{f}$ $\displaystyle=\Sigma_{\xi}.$ (42)

In order to show that (40) holds, we note that

$\norm{\int_{0}^{T}Y_{s}ds-\sum_{k=1}^{n}\xi_{k}}=\norm{\int_{n}^{T}Y_{s}ds}\leq\int_{n}^{n+1}\norm{Y_{s}}ds.$
(43)

By a Borel-Cantelli argument it will follow that

$\norm{\int_{n}^{T}f(Z_{s})\ ds}=O(T^{1/(2+\delta)})\quad\textrm{a.s. for}\
T\rightarrow\infty.$ (44)

Indeed, let $\varepsilon>0$ be given and introduce the event

$A_{n,\varepsilon}=\left\\{\int_{n}^{n+1}\norm{Y_{s}}ds>n^{(1+\varepsilon)/(2+\delta)}\right\\}.$

By Markov’s inequality it follows that the introduced sequence of events
satisfies

$\sum_{n=1}^{\infty}\Pr\left(A_{n,\varepsilon}\right)\leq\sum_{n=1}^{\infty}\Pr\left(\int_{n}^{n+1}\norm{Y_{s}}^{2+\delta}ds>n^{1+\varepsilon}\right)\leq\pi(\norm{f-\pi(f)})^{2+\delta}\sum_{n=1}^{\infty}\frac{1}{n^{1+\varepsilon}}<\infty.$

The Borel-Cantelli lemma states that $\mathbb{P}(\limsup
A_{n,\varepsilon})=0$. In order for (44) to hold, it is sufficient to note
that from a certain point in time the remainder term will always satisfy the
asserted bound for almost every sample path, i.e.,

$\displaystyle\Pr\left(\bigcup_{s\geq 0}\bigcap_{T\geq
s}\left\\{\norm{\int_{n}^{T}Y_{u}du}\leq T^{1/(2+\delta)}\right\\}\right)$
$\displaystyle\geq\Pr\left(\bigcap_{m=1}^{\infty}\liminf_{n\rightarrow\infty}A_{n,1/m}^{c}\right)=\lim_{m\rightarrow\infty}\Pr\left(\liminf_{n\rightarrow\infty}A_{n,1/m}^{c}\right)=1,$

where the first inequality follows due to monotonicity and the second-last
equality by the continuity of measures. A similar Borel-Cantelli argument also
shows that (41) holds. Introduce sequence of events

$B_{n,\varepsilon}=\left\\{\norm{W(T)-W(n)}>n^{(1+\varepsilon)/\kappa}\right\\},$

for given $\varepsilon>0$ and some $\kappa>2+\delta$. Since all moments of
$\norm{W(T)-W(n)}$ are finite, we have by Markov’s inequality that the
introduced sequence of events satisfies

$\sum_{n=1}^{\infty}\Pr\left(B_{n,\varepsilon}\right)\leq\sum_{n=1}^{\infty}\Pr\left(\norm{B_{T}-B_{n}}^{\kappa}>n^{1+\varepsilon}\right)\leq\mathbb{E}_{\Pr}\norm{W(T)-W(n)}^{\kappa}\sum_{n=1}^{\infty}\frac{1}{n^{1+\varepsilon}}<\infty.$

By completely analogous reasoning as the previous Borel-Cantelli argument we
see that

$\norm{W(T)-W(n)}=O(n^{1/\kappa})=o(T^{1/(2+\delta)})\ \textrm{a.s. for}\
T\rightarrow\infty.$

Therefore the term (41) will be asymptotically negligible. Finally, we see
that by Lemma 6.3 it immediately follows that the asserted asymptotic variance
$\Sigma_{f}$ is finite, i.e., all entries

$(\Sigma_{f})_{ij}=\int_{0}^{\infty}\operatorname*{Cov_{\pi}}(f_{i}(X_{0}),f_{j}(X_{s}))\
ds+\int_{0}^{\infty}\operatorname*{Cov_{\pi}}(f_{i}(X_{s}),f_{j}(X_{0}))\ ds\
\ \textrm{for}\ 1\leq i,j\leq p$ (45)

converge absolutely. Indeed,

$\displaystyle\int_{0}^{\infty}\lvert\operatorname*{Cov_{\pi}}(f_{i}(X_{0}),f_{j}(X_{s}))\rvert
ds$ $\displaystyle\leq
8\int_{0}^{\infty}\alpha_{Y}(s)^{\delta/(2+\delta)}\pi(\lvert
Y_{i0}\rvert^{2+\delta})^{1/(2+\delta)}\pi(\lvert
Y_{js}\rvert^{2+\delta})^{1/(2+\delta)}ds$ $\displaystyle\leq
8\pi(V)^{\delta/(2+\delta)}\pi(\lvert
Y_{i0}\rvert^{2+\delta})^{1/(2+\delta)}\pi(\lvert
Y_{j0}\rvert^{2+\delta})^{1/(2+\delta)}\int_{0}^{\infty}\alpha_{X}(s)^{\delta/(2+\delta)}ds<\infty.$

By definition $\alpha$-mixing sequences are monotonically decreasing and
bounded by $1/4$. The second term of (45) is treated similarly. In order to
show that $\Sigma_{f}=\Sigma_{\xi}$, we will show that all entries are equal.

Firstly, we decompose the asymptotic covariance matrix as follows

$\displaystyle\lim_{T\rightarrow\infty}\operatorname*{Var_{\pi}}\left(\frac{1}{\sqrt{T}}\int_{0}^{T}Y_{t}dt\right)=\lim_{T\rightarrow\infty}\operatorname*{Var_{\pi}}\left(\frac{1}{\sqrt{T}}\left(\int_{0}^{n}Y_{t}dt+\int_{n}^{T}Y_{t}dt\right)\right)=$
$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\operatorname*{Var_{\pi}}\left(\int_{0}^{n}Y_{t}dt\right)+\lim_{T\rightarrow\infty}\frac{1}{T}\operatorname*{Var_{\pi}}\left(\int_{n}^{t}Y_{t}dt\right)+\lim_{T\rightarrow\infty}\frac{1}{T}\operatorname*{Cov_{\pi}}\left(\int_{0}^{n}Y_{t}dt,\int_{n}^{t}Y_{t}dt\right)+$
$\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\operatorname*{Cov_{\pi}}\left(\int_{n}^{t}Y_{t}dt,\int_{0}^{n}Y_{t}dt,\right).$
(46)

Let $\Sigma_{T1},\Sigma_{T2},\Sigma_{T3}$ and $\Sigma_{T4}$ denote the four
terms in (6.1.1). We will show that the entry-wise convergence gives us the
desired result. For $1\leq i,j\leq p$ we obtain the following expressions for
the elements of the matrices in (6.1.1):

$\displaystyle(\Sigma_{T1})_{ij}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{n}\int_{0}^{n}\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})\
dtds$ (47) $\displaystyle(\Sigma_{T2})_{ij}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{n}^{T}\int_{n}^{T}\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})\
dtds$ (48) $\displaystyle(\Sigma_{T3})_{ij}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{n}\int_{n}^{T}\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})\
dtds$ (49) $\displaystyle(\Sigma_{T4})_{ij}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{n}\int_{n}^{T}\operatorname*{Cov_{\pi}}(Y_{is},Y_{jt})\
dtds$ (50)

We see that $(\Sigma_{T1})_{ij}$ tends to the asymptotic variance
$(\Sigma_{f})_{ij}$, since

$\left(\operatorname*{Var_{\pi}}\left(\frac{1}{\sqrt{n}}\sum_{k=1}^{n}\xi_{k}\right)\right)_{ij}=\frac{T}{n}\cdot\frac{1}{T}\int_{0}^{n}\int_{0}^{n}\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})\
dtds.$

Finally, another application of Lemma 6.3 gives us that
$(\Sigma_{T2})_{ij},(\Sigma_{T3})_{ij}$ and $(\Sigma_{T4})_{ij}$ tend to zero
since

$\displaystyle\frac{1}{T}\int_{n}^{T}\int_{n}^{T}\absolutevalue{\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})}\
dtds$
$\displaystyle\leq\frac{1}{T}C_{f,V}\int_{n}^{T}s^{-\beta\delta/(2+\delta)}ds=o(1),$
$\displaystyle\frac{1}{T}\int_{0}^{n}\int_{n}^{T}\absolutevalue{\operatorname*{Cov_{\pi}}(Y_{it},Y_{js})}\
dtds$
$\displaystyle\leq\frac{n}{T}C_{f,V}\int_{n}^{T}s^{-\beta\delta/(2+\delta)}\
ds=o(1),\ $
$\displaystyle\frac{1}{T}\int_{0}^{n}\int_{n}^{T}\absolutevalue{\operatorname*{Cov_{\pi}}(Y_{is},Y_{jt})}\
dtds$
$\displaystyle\leq\frac{n}{T}C_{f,V}\int_{n}^{T}s^{-\beta\delta/(2+\delta)}\
ds=o(1),$

with $C_{f,V}=8\pi(V)^{\delta/(2+\delta)}\pi(\lvert
Y_{i0}\rvert^{2+\delta})^{1/(2+\delta)}\pi(\lvert
Y_{j0}\rvert^{2+\delta})^{1/(2+\delta)}<\infty.$

Note that we have now shown our result assuming stationarity, i.e., with
initial distribution $\pi$. However, by a standard argument this can be
extended to every initial distribution. By Kallenberg [45; Theorem 20.10] we
have that for ergodic Markov processes every bounded harmonic function is
constant. By following the argument of Meyn and Tweedie [57; Proposition
17.1.6] it follows that the strong invariance principle holds for every
initial distribution. ∎

###  6.2 Theorem 3.3

Berkes et al. [4] give a strong invariance principle for weakly $m$-dependent
processes, which are defined as processes that can be approximated by
$m$-dependent processes in the $L^{p}$-sense, with a sufficiently decaying
approximation error (rate function in terminology of Berkes et al. [4]). Their
strong invariance principle, stated in Theorem 6.4, is obtained through a
classical blocking argument for $m$-dependent random variables. By dividing an
$m$-dependent sequence into non-overlapping long and short blocks, two
sequences of independent random variables are obtained; these can both be
approximated by a Brownian motion. Trivially, stationary $m$-dependent
processes satisfying appropriate moment conditions fall into their framework.
For more details we refer to Berkes et al. [4].

######  Theorem 6.4 (Berkes et al. [4; Theorem 2]).

Let $\xi=(\xi_{k})_{k=1}^{\infty}$ be a centered stationary sequence with
$\sup_{k}\mathbb{E}\absolutevalue{\xi_{k}}^{2+\delta}<\infty$, for some
$\delta>0$. Moreover, let $\xi$ be weakly $m$-dependent in $L^{p}$ with an
exponentially decaying rate function $\kappa$, i.e.,

$\kappa(m)\ll\exp(-cm),\quad\textrm{for some }c>0.$

Then the series

$\sigma_{\xi}^{2}=\sum_{k=0}^{\infty}\mathbb{E}\xi_{0}\xi_{k}$

is absolutely convergent, and we can redefine $\xi$ on a new probability space
on which we can also construct two standard Brownian motions $W_{1}$ and
$W_{2}$ such that

$\absolutevalue{\sum_{k=1}^{n}\xi_{k}-W_{1}(s_{n}^{2})-W_{2}(t_{n}^{2})}=O(n^{1/p}\log^{2}n)\
\textrm{a.s.,}$

where $\\{s_{n}\\}$ and $\\{t_{n}\\}$ are non-decreasing deterministic
sequences with $s_{n}^{2}=\sigma_{\xi}^{2}n+O(n/\log n)$ and
$t_{n}^{2}=O(n/\log n)$.

As noted by Berkes et al. [4] the perturbed time sequences $\\{s_{n}\\}$ and
$\\{t_{n}\\}$ are deterministic and can be explicitly calculated.

####  6.2.1 Proof of Theorem 3.3

###### Proof.

By Proposition 2.1 we see that we can redefine our process such that it is
embedded in a richer process $Z$. We will identify $X$ as the first coordinate
of the process $Z.$ Following Proposition 2.2, we introduce the sequence of
stopping times $(S_{n},R_{n})$ defined as $S_{0}=R_{0}:=0$ and

$S_{n+1}:=\inf\\{T_{m}>R_{n}:Z_{T_{m}}\in A\\}\ \ \textrm{and}\ \
R_{n+1}:=\inf\\{T_{m}:T_{m}>S_{n+1}\\}.$

Then $Z_{R_{n}^{+}}$ is independent of $\mathcal{F}_{R_{n-1}}$ for all $n\geq
1$ and $(Z_{R_{n}})_{n\geq 1}$ is an i.i.d sequence with

$Z_{R_{n}}\sim\nu(dx)\lambda(du)K((x,u),dx^{\prime})\ \ \textrm{for all}\
n\geq 1.$

As a direct consequence, the sequence $\\{\xi_{n}\\}_{n}$ defined as

$\xi_{n}:=\int_{R_{n-1}}^{R_{n}}\\{f(X_{s})-\pi(f)\\}\ ds,\ \ n\geq 1,$

is stationary under $\mathbb{P}_{\nu}.$ Moreover, by Proposition 2.3 for
$n\geq 2$, $\xi_{n}$ is independent of $\mathcal{F}_{R_{n-2}}.$

Let $N(T)$ denote the amount of regenerations of the resolvent chain up to
time $T$, namely

$N(T)=\max\\{k:R_{k}\leq T\\}.$

It immediately follows that

$\int_{0}^{T}\\{f(X_{s})-\pi(f)\\}\
ds=\sum_{k=1}^{N(T)}\xi_{k}+\int_{N(T)}^{T}\\{f(X_{s})-\pi(f)\\}\ ds.$

Consequently, we have that

$\Bigg{|}\int_{0}^{T}\\{f(X_{s})-\pi(f)\\}\
ds-\sum_{k=1}^{N(T)}\xi_{k}\Bigg{|}\leq\int_{N(T)}^{T}\absolutevalue{f(X_{s})-\pi(f)}ds.$
(51)

By an analogous Borel-Cantelli argument as given in Theorem 3.1 for the
remainder term defined in (43), we have that

$\int_{N(T)}^{T}\absolutevalue{f(X_{s})-\pi(f)}ds=O(T^{1/(2+\delta)})\quad\textrm{a.s.
for}\ T\rightarrow\infty.$ (52)

Furthermore, by Proposition 2.3, the sequence $\\{\xi_{k}\\}_{k=1}^{\infty}$
is a stationary $m$-dependent sequence. By the imposed moment conditions and
stationarity, we have by the reasoning given in Berkes et al. [4; Section 3.1]
that $\\{\xi_{k}\\}_{k=1}^{\infty}$ is also a weakly $m$-dependent process
with a rate function $\kappa(m)$ equal to zero for $m\geq 1.$ Hence by Theorem
6.4, we can redefine $(\xi_{k})_{k}$ on a new probability space on which we
can also construct two standard Brownian motions $W_{1}$ and $W_{2}$ such that

$\absolutevalue{\sum_{k=1}^{n}\xi_{k}-n\mathbb{E}_{\nu}\xi_{1}-W_{1}(s_{n}^{2})-W_{2}(t_{n}^{2})}=O(n^{1/p}\log^{2}n)\
\textrm{a.s.,}$ (53)

where $\\{s_{n}\\}$ and $\\{t_{n}\\}$ are increasing deterministic sequences
with $s_{n}^{2}=\sigma_{\xi}^{2}n+O(n/\log n)$ and $t_{n}^{2}=O(n/\log n)$.
Note that by Proposition 2.4 we have that

$\pi(f)=\frac{1}{\varrho}\ \mathbb{E}_{\nu}\int_{0}^{R_{1}}f(X_{s})ds.$

Hence

$\mathbb{E}_{\nu}\xi_{1}=\mathbb{E}_{\nu}\int_{0}^{R_{1}}\\{f(X_{s})-\pi(f)\\}\
ds=\varrho\cdot\pi(f-\pi(f))=0.$

Furthermore, by definition of big $O$ in (53), there exists an almost surely
finite random variable $C$ such that for almost all sample paths $\omega$ we
have that $\textrm{for all }n\geq N_{0}\equiv N_{0}(\omega)$ we have that

$\left.\frac{1}{n^{1/p}\log^{2}n}\absolutevalue{\sum_{k=1}^{n}\xi_{k}(\omega)-W_{1}(s_{n}^{2},\omega)-W_{2}(t_{n}^{2},\omega)}\right.<C(\omega)$
(54)

From standard results in renewal theory, the Marcinkiewicz-Zygmund law for
first passage times and the corresponding central limit theorem, see for
example Gut [39; Chapter 3: Theorem 4.4 and 5.1], it follows that the amount
of regeneration times of the resolvent displays the following asymptotic
behaviour

$N(T)=\frac{T}{\varrho}+O(T^{1/q})\quad\textrm{a.s. for
}T\rightarrow\infty\textrm{ and }\ 1<q\leq 2.$ (55)

Since $N(T)$ is almost surely increasing and tends to infinity, we have that
for almost every sample path $\omega$ there exists a $T_{0}\equiv
T_{0}(\omega)$ such that $N(T)(\omega)\geq N_{0}$ for all $T\geq T_{0}$. Hence
we obtain from (54) that

$\left.\limsup_{T\rightarrow\infty}\frac{\left\lvert\sum_{k=1}^{N(T)}\xi_{k}-W_{1}(s_{N(T)}^{2})-W_{2}(t_{N(T)}^{2})\right\lvert}{{N(T)}^{1/p}\log^{2}(N(T))}\right.<C\quad\textrm{a.s.
},$ (56)

where $s_{N(T)}^{2}$ and $t_{N(T)}^{2}$ are almost surely increasing
sequences, which given $N(T)$ are deterministic with

$\displaystyle s_{N(T)}^{2}$ $\displaystyle=\sigma_{\xi}^{2}N(T)+O(N(T)/\log
N(T))$ $\displaystyle t_{N(T)}^{2}$ $\displaystyle=O(N(T)/\log N(T)).$

We see that (56) can be reformulated as

$\displaystyle\left\lvert\sum_{k=1}^{N(T)}\xi_{k}-W_{1}(s^{2}_{N(T)})-W_{2}(t_{N(T)}^{2})\right\lvert$
$\displaystyle=O({N(T)}^{1/p}\log^{2}N(T)))\quad\textrm{a.s. }$ (57)
$\displaystyle=O({T}^{1/p}\log^{2}T))\quad\textrm{a.s. }$ (58)

Here the second equality follows by (55). Furthermore, the asymptotic
behaviour of $N(T)$ motivates the introduction of
$\sigma^{2}_{T},\tau^{2}_{T}$ and $\beta_{T}$ defined as

$\displaystyle\sigma^{2}_{T}$
$\displaystyle=s^{2}_{N(T)}-\sigma^{2}_{\xi}\beta_{T}$
$\displaystyle\tau^{2}_{T}$ $\displaystyle=t^{2}_{N(T)}-\beta_{T}$
$\displaystyle\beta_{T}$ $\displaystyle=N(T)-T/\varrho$

From (55) we see that $\beta_{T}=O(T^{1/q})$, i.e.,
$\limsup\beta_{T}/T^{1/q}<\infty$ a.s. Hence by Theorem 5.1 (see also Theorem
1.2.1 of Csörgö and Révész [19]) we see that

$\big{\lvert}W_{1}(s^{2}_{N(T)})-W_{1}(\sigma^{2}_{T})\big{\lvert}=O(T^{1/(2q)}\log
T)\ \textrm{and}\
\big{\lvert}W_{2}(t^{2}_{N(T)})-W_{1}(\tau^{2}_{T})\big{\lvert}=O(T^{1/(2q)}\log
T),$ (59)

with

$\sigma^{2}_{T}=\frac{\sigma^{2}_{\xi}}{\varrho}T+O(T/\log T)\ \textrm{and}\
\tau_{T}^{2}=O(T/\log T).$

Combining results (52), (57), and (59) concludes the proof. By the same
arguments given in the proof of Theorem 3.1, the strong invariance principle
holds for every initial distribution. ∎

Furthermore, we have by Berkes et al. [4; Proposition 1] that

$\operatorname*{Corr}(W_{1}(s_{n}),W_{2}(t_{m}))\rightarrow 0\quad\textrm{as}\
m,n\rightarrow\infty.$

Hence (26) immediately follows.

###  6.3 Theorem 3.4

Note that by Proposition 2.1 $(Z^{1}_{T_{n}},Z^{2}_{T_{n}})_{n}$, the jump
chain of the first two coordinates of $Z$, has the same distribution as the
split chain of the resolvent. From (7) and (8) we see that
$(Z^{1}_{T_{n}},Z^{2}_{T_{n}})_{n}$ is a Markov chain taking values in
$E^{\prime}:=E\times[0,1]$ that moves according to the kernel

$U^{\prime}((x,u),(dx,du))=\nu(dx)\lambda(du)\mathbbm{1}_{\\{u\leq\alpha\mathbbm{1}_{C}(x)\\}}+W(x,dx)\lambda(du)\mathbbm{1}_{\\{u>\alpha\mathbbm{1}_{C}(x)\\}},$
(60)

where $\lambda$ denotes Lebesgue measure on the unit interval. Observe that
the kernel of the split chain of the resolvent also satisfies a one-step
minorisation condition $U^{\prime}\geq s\otimes\nu\otimes\lambda,$ i.e.,

$U^{\prime}((x,u),(dx,du))\geq s(x,u)\nu(dx)\lambda(du),$ (61)

where

$s(x,u)=\mathbbm{1}_{\\{u\leq\alpha\mathbbm{1}_{C}(x)\\}}.$

Moreover, the split chain of the resolvent is aperiodic, positive Harris
recurrent and inherits the rate of convergence to stationarity from $X$.

######  Lemma 6.5 (Hobert et al. [41; Lemma 1]).

Let $(X_{t})_{t\geq 0}$ be a positive Harris recurrent Markov process with
invariant distribution $\pi$. Then for any $\pi$-integrable function
$g:E^{[0,\infty)}\rightarrow\mathbb{R}$ we have the following inequality holds

$\mathbb{E}_{\pi}\absolutevalue{g}\geq c\ \mathbb{E}_{\nu}\absolutevalue{g},$
(62)

where $c=\alpha\pi(C)$.

###### Proof.

Since the resolvent chain has the same stationary distribution as the process
$X$, i.e., $\pi=\pi U$, the claim follows with the identical argument of
Hobert et al. [41; Lemma 1]. ∎

####  6.3.1 Proof of Theorem 3.4

###### Proof.

Firstly, by the construction of the randomised stopping times $(S_{n})_{n}$
and $(R_{n})_{n}$ we see that $R_{n}=S_{n}+\sigma_{n+1}$, where $\sigma_{n+1}$
has a standard exponential distribution. Hence, by the triangle inequality in
$L^{q}(\pi)$ we only need to show that $\mathbb{E}_{\pi}[{S_{1}}^{q}]<\infty$,
with

$S_{1}=\inf\\{T_{n}:Z_{T_{n}}\in C\times[0,\alpha]\times E\\}$

Let $\bar{Z}=(\bar{Z}_{n})_{n}$ denote the jump chain of the process $Z$,
i.e., $\bar{Z}_{n}=Z_{T_{n}}$, where the $(T_{n})_{n}$ denote the jump times.
Let $\bar{X}=(\bar{X}_{n})_{n\geq 0}$ again denote the resolvent chain. Let
$N_{t}$ denote the amount of jumps up to time $t$. Let $\bar{\tau}_{A}$ denote
the hitting time of the recurrent atom for jump chain $\bar{Z}$, i.e.,

$\displaystyle\bar{\tau}_{A}:$ $\displaystyle=\inf\\{n\geq 0:\bar{Z}_{n}\in
A\\}=\inf\\{n\geq 0:\bar{Z}_{n}\in C\times[0,\alpha]\times E\\}.$

Employing the relationship between the expectation of positive random
variables and tail probabilities we can express the expectation of interest as
follows

$\displaystyle\mathbb{E}_{\pi}S_{1}^{q}$
$\displaystyle=\int_{0}^{\infty}t^{q}\mathbb{P}_{\pi}\left(S_{1}>t\right)dt$
$\displaystyle=\int_{0}^{\infty}t^{q}\sum_{m=0}^{\infty}\mathbb{P}_{\pi}\left(\bar{\tau}_{A}>m;N_{t}=m\right)dt$
$\displaystyle=\int_{0}^{\infty}t^{q}\sum_{m=0}^{\infty}\underbrace{\mathbb{P}_{\pi}\left(\bar{\tau}_{A}>m;N_{t}=m;\bar{Z}_{0}\in
C\right)}_{=0}+\
\mathbb{P}_{\pi}\left(\bar{\tau}_{A}>m;N_{t}=m;\bar{Z}_{0}\notin C\right)dt$
$\displaystyle=\int_{0}^{\infty}t^{q}\int_{E^{\prime}}\int_{E^{\prime}}\sum_{m=0}^{\infty}\frac{t^{m-1}}{(m-1)!}e^{-t}\sum_{k=m+1}^{\infty}\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dz)\mathbbm{1}_{A}(z)\pi(dx)dt$
$\displaystyle=\int_{E^{\prime}}\int_{E^{\prime}}\sum_{m=0}^{\infty}\int_{0}^{\infty}\frac{t^{m+q-1}}{(m-1)!}e^{-t}\sum_{k=m+1}^{\infty}\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dz)\mathbbm{1}_{A}(z)dt\pi(dx)$
$\displaystyle=\int_{E^{\prime}}\int_{E^{\prime}}\sum_{k=1}^{\infty}\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dz)\mathbbm{1}_{A}(z)\sum_{m=0}^{k-1}\frac{\Gamma(m+q)}{\Gamma(m)}\pi(dx)$
$\displaystyle=\int_{E^{\prime}}\int_{E^{\prime}}\sum_{k=1}^{\infty}\frac{\Gamma(k+q)}{(q+1)\Gamma(k-1)}\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dz)s(z)\pi(dx).$

Note that ${\Gamma(k+q)}/{\Gamma(k-1)}$ can be dominated by some polynomial
$\psi(k)$ with a leading term of order $k^{q+1}$. By Nummelin and Tuominen
[63; Proposition 1.6] we have that

$\int_{E^{\prime}}\int_{E^{\prime}}\sum_{k=0}^{\infty}\psi(k)\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dz)s(z)\pi(dx)<\infty.$

It follows that $\mathbb{E}_{\pi}S_{1}^{q}<\infty.$

The second statement of Theorem 3.4 follows according to the argument of
Bednorz and Łatuszyński [3; Lemma 2] with some minor adaptations. We give the
proof for completion.

$\displaystyle\left[\mathbb{E}_{\pi}\xi_{1}^{p}\right]^{1/p}$
$\displaystyle\leq\left[\mathbb{E}_{\pi}{\left(\int_{0}^{R_{1}}\absolutevalue{f(X_{s})}ds\right)}^{p}\right]^{1/p}$
$\displaystyle=\left[\mathbb{E}_{\pi}\left(\int_{0}^{\infty}\absolutevalue{f(X_{s})}\mathbbm{1}_{\\{R_{1}\geq
s\\}}ds\right)^{p}\right]^{1/p}$
$\displaystyle\leq\int_{0}^{\infty}\left[\mathbb{E}_{\pi}\left(\absolutevalue{f(X_{s})}^{p}\mathbbm{1}_{\\{R_{1}\geq
s\\}}ds\right)\right]^{1/p}$
$\displaystyle\leq\int_{0}^{\infty}\left[\mathbb{E}_{\pi}\absolutevalue{f(X_{s})}^{p+\varepsilon}\right]^{p/(p+\varepsilon)}\left[\mathbb{E}_{\pi}\mathbbm{1}_{\\{R_{1}\geq
s\\}}\right]^{\varepsilon/(p+\varepsilon)}ds$
$\displaystyle\leq\pi\left(\absolutevalue{f}^{p+\varepsilon}\right)^{p/(p+\varepsilon)}\int_{0}^{\infty}\left[\mathbb{P}_{\pi}(R_{1}\geq
s)\right]^{\varepsilon/(p+\varepsilon)}ds$
$\displaystyle\leq\pi\left(\absolutevalue{f}^{p+\varepsilon}\right)^{p/(p+\varepsilon)}(1+\pi(R_{1}^{q}))\int_{1}^{\infty}s^{-\varepsilon
q/(p+\varepsilon)}ds<\infty.$

Here the inequalities follow by Minkowski’s integral inequality, Hölder’s
inequality, stationarity, and Markov’s inequality. Note that the integral on
the last line is finite due to the imposed condition on the rate of polynomial
ergodicity $q$. An application of Lemma 6.5 concludes the proof.

∎

######  Remark 6.6.

For the exponentially erogdic case we would make use of Nummelin and Tuominen
[62; Lemma 2.8] which states that for an exponentially ergodic Markov chain
there exists an $r>1$ such that

$\int_{E^{\prime}}\int_{E^{\prime}}\sum_{k=0}^{\infty}r^{k}\left(U^{\prime}-\nu\otimes\lambda\otimes
s\right)^{k}(x,dy)\mathbbm{1}_{C}(y)\pi(dx)<\infty.$

∎

###  6.4 Proof of Theorem 5.2

###### Proof.

Firstly, by Theorem 19 there exist two standard Brownian motions $W_{1}$ and
$W_{2}$ such that

$\absolutevalue{\int_{0}^{T}f(X_{s})ds-
W_{1}(\sigma^{2}_{T})-W_{2}(\tau^{2}_{T})}=O(\psi_{T})\ \textrm{a.s.,}$

where $\\{\sigma^{2}_{T}\\}$ and $\\{\tau^{2}_{T}\\}$ are non-decreasing
sequences with

$\sigma^{2}_{T}=\frac{\sigma^{2}_{\xi}}{\varrho}T+O(T/\log T)\ \textrm{and}\
\tau_{T}^{2}=O(T/\log T).$

as $T\rightarrow\infty$, where $\sigma^{2}_{\xi}\textrm{ and}\ \varrho$ are
defined in Theorem 19. An application of our strong invariance principle gives
the following

$\displaystyle\limsup_{T\rightarrow\infty}\max_{0\leq t\leq
T-a_{T}}\max_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{\int_{0}^{t+u}f(X_{u})du-\int_{0}^{t}f(X_{u})du}$
$\displaystyle\leq\limsup_{T\rightarrow\infty}\max_{0\leq t\leq
T-a_{T}}\max_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{W_{1}(\sigma^{2}_{t+u})-W_{1}(\sigma^{2}_{t})}$
$\displaystyle+\limsup_{T\rightarrow\infty}\max_{0\leq t\leq
T-a_{T}}\max_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{W_{2}(\tau^{2}_{t+u})-W_{2}(\tau^{2}_{t})}$
$\displaystyle+\beta_{T}O(\psi_{T})$ $\displaystyle=:A_{1}+A_{2}+A_{3}.$

Since $\beta_{T}\psi_{T}=o(1)$, it immediately follows that
$\limsup_{T}A_{3}=0$ almost surely. In order to use the arguments of Berkes et
al. [4; Theorem 4] for the terms $A_{1}$ ad $A_{2}$, we require the following
properties of the sequence $\sigma^{2}_{T}$; for any $\varepsilon>0$ there
exists some $T_{0}$ such that for all $T\geq T_{0}$

$\sigma^{2}_{T}\leq\left(\frac{\sigma^{2}_{\xi}}{\varrho}+\varepsilon\right)T\quad\textrm{and}\quad\sup_{u\geq
0}\\{\sigma^{2}_{u+a_{T}}-\sigma^{2}_{u}\\}\leq\left(\frac{\sigma^{2}_{\xi}}{\varrho}+\varepsilon\right)a_{T}.$
(63)

From the proof of Theorem 3.3 we see that the process
$\\{\xi_{n}\\}_{n\in\mathbb{N}}$ with
$\xi_{n}:=\int_{R_{k-1}}^{R_{k}}f(X_{s})\ ds$, where $\\{R_{k}\\}_{k}$ denotes
the regeneration times of the resolvent, satisfies the following strong
approximation result

$\absolutevalue{\sum_{k=1}^{n}\xi_{k}-W_{1}(s_{n}^{2})-W_{2}(t_{n}^{2})}=O(n^{1/p}\log^{2}n)\
\textrm{a.s.,}$

where $\\{s_{n}\\}$ and $\\{t_{n}\\}$ are increasing deterministic sequences
with $s_{n}^{2}=\sigma_{\xi}^{2}n+O(n/\log n)$ and
$t_{n}^{2}\sim\sigma_{\xi}^{2}n/\log n$. By the properties of the sequences
$s^{2}_{n}$ and $t^{2}_{n}$ described in Theorem 6.4 we have that for any
$\varepsilon>0$ there exists some $n_{0}$ such that for all $n\geq n_{0}$

$s^{2}_{n}\leq(\sigma^{2}_{\xi}+\varepsilon)n\ \quad\textrm{and}\
\quad\sup_{k}\\{s^{2}_{k+{a_{n}}}-s^{2}_{k}\\}\leq(\sigma^{2}_{\xi}+\varepsilon)a_{n}$

Since $\sigma^{2}_{T}=\frac{\sigma^{2}_{\xi}}{\mu}T+O(T/\log(T))$, the first
required property described in (63) can be easily seen to hold. Since
$\limsup_{k\rightarrow\infty}(s^{2}_{k+1}-s^{2}_{k})=\sigma^{2}_{\xi}$ and
$N(u)$ tends to infinity almost surely as $u\rightarrow\infty$, it immediately
follows that

$\limsup_{u\rightarrow\infty}(s^{2}_{N(u)+1}-s^{2}_{N(u)})=\sigma^{2}_{\xi}.$

Hence for $u$ sufficiently large we have that

$s^{2}_{N(u+a_{T})}-s^{2}_{N(u)}=\sum_{j=1}^{N(u+a_{T})-N(u)}s^{2}_{N(u)+j}-s^{2}_{N(u)+j-1}\leq(\sigma^{2}_{\xi}+\varepsilon)(N(u+a_{T})-N(u))$

Since $\sigma^{2}_{T}$ was defined as
$\sigma^{2}_{T}:=s^{2}_{N(T)}-\sigma^{2}_{\xi}\beta_{T}$ with
$\beta_{T}:=N(T)-T/\varrho$ (see proof of Theorem 3.3),

we have that for all $\varepsilon>0$ and $u$ sufficiently large

$\displaystyle\sigma^{2}_{u+a_{T}}-\sigma^{2}_{u}$
$\displaystyle=s^{2}_{N(u+a_{T})}-s^{2}_{N(u)}-\sigma^{2}_{\xi}(\beta_{u+a_{T}}-\beta_{u})$
$\displaystyle\leq(\sigma^{2}_{\xi}+\varepsilon)(N(u+a_{T})-N(u))-\sigma^{2}_{\xi}(N(u+a_{T})-N(u))+\frac{\sigma^{2}_{\xi}}{\varrho}a_{T}$
$\displaystyle=\varepsilon(N(u+a_{T})-N(u))+\frac{\sigma^{2}_{\xi}}{\varrho}a_{T}.$

Since $a_{T}\rightarrow\infty$ as $T\rightarrow\infty$ by assumption, and we
have that the amount of regenerations of the resolvent chain in an interval of
length $a_{T}$ will be equal to $a_{T}/\varrho+o(a_{T})$, it follows that for
all $\varepsilon>0$ there exists some $T_{0}$ such that for all $T\geq T_{0}$

$\displaystyle\sup_{u\geq
0}\\{\sigma^{2}_{u+a_{T}}-\sigma^{2}_{u}\\}\leq\left(\frac{\sigma^{2}_{\xi}}{\varrho}+2\varepsilon\right)a_{T}.$

Consequently, we have also shown that the required properties given in (63)
hold. Hence for $T\geq T_{0}$ we obtain

$\displaystyle\max_{0\leq t\leq T-a_{T}}\max_{0\leq u\leq
a_{T}}\beta_{T}\absolutevalue{W_{1}(\sigma^{2}_{t+u})-W_{1}(\sigma^{2}_{t})}$
$\displaystyle\leq\sup_{0\leq t\leq\sigma^{2}_{T-a_{T}}}\sup_{0\leq
u\leq(\sigma^{2}_{\xi}/\varrho+\varepsilon)a_{T}}\beta_{T}\absolutevalue{W_{1}(t+u)-W_{1}(t)}$
$\displaystyle\leq\sup_{0\leq
t\leq(\sigma^{2}_{\xi}/\varrho+\varepsilon)(T-a_{T})}\sup_{\ 0\leq
u\leq(\sigma^{2}_{\xi}/\varrho+\varepsilon)a_{T}}\beta_{T}\absolutevalue{W_{1}(t+u)-W_{1}(t)}$
$\displaystyle=\sup_{0\leq
t\leq\tilde{T}_{\varepsilon}-\tilde{a}_{T,\varepsilon}}\sup_{\ 0\leq
u\leq\tilde{a}_{T,\varepsilon}}\beta_{T}\absolutevalue{W_{1}(t+u)-W_{1}(t)},$

where $\tilde{T}_{\varepsilon}$ and $\tilde{a}_{T,\varepsilon}$ are defined as
$(\sigma^{2}_{\xi}/\varrho+\varepsilon)T$ and
$(\sigma^{2}_{\xi}/\varrho+\varepsilon)a_{T}$ respectively. Introduce

$\tilde{\beta}_{T,\varepsilon}:=\left(2\tilde{a}_{T,\varepsilon}\left[\log\frac{\tilde{T}_{\varepsilon}}{\tilde{a}_{T,\varepsilon}}+\log\log\tilde{T}_{\varepsilon}\right]\right)^{-1/2},$

then by Theorem 5.1 we have that

$\limsup_{T\rightarrow\infty}\sup_{0\leq
t\leq\tilde{T}_{\varepsilon}-a_{T,\varepsilon}}\sup_{0\leq u\leq
a_{T,\varepsilon}}\tilde{\beta}_{T,\varepsilon}\absolutevalue{W(t+u)-W_{t}}={\sigma^{2}_{\xi}}/{\varrho}\quad\textrm{a.s.}$

Similarly, it can be shown that $\limsup A_{2}=0$ almost surely, which
completes the proof.
∎

## Acknowledgements

This work is part of the research programme ‘Zigzagging through computational
barriers’ with project number 016.Vidi.189.043, which is financed by the Dutch
Research Council (NWO).

## References

  * Andrieu et al. [2021] C. Andrieu, P. Dobson, and A. Q. Wang.  Subgeometric hypocoercivity for piecewise-deterministic markov process monte carlo methods.  _Electronic Journal of Probability_ , 26:1–26, 2021.
  * Athreya and Ney [1978] K. B. Athreya and P. Ney.  A new approach to the limit theory of recurrent markov chains.  _Transactions of the American Mathematical Society_ , 245:493–501, 1978.
  * Bednorz and Łatuszyński [2007] W. Bednorz and K. Łatuszyński.  A few remarks on “fixed-width output analysis for markov chain monte carlo” by jones et al.  _Journal of the American Statistical Association_ , 102(480):1485–1486, 2007.
  * Berkes et al. [2011] I. Berkes, S. Hörmann, J. Schauer, et al.  Split invariance principles for stationary processes.  _The Annals of Probability_ , 39(6):2441–2473, 2011.
  * Berkes et al. [2014] I. Berkes, W. Liu, and W. B. Wu.  Komlós–major–tusnády approximation under dependence.  _The Annals of Probability_ , 42(2):794–817, 2014\.
  * Bhattacharya [1978] R. Bhattacharya.  Criteria for recurrence and existence of invariant measures for multidimensional diffusions.  _The Annals of Probability_ , pages 541–553, 1978.
  * Bierkens and Roberts [2017] J. Bierkens and G. Roberts.  A piecewise deterministic scaling limit of lifted metropolis–hastings in the curie-weiss model.  _The Annals of Applied Probability_ , 27(2):846–882, 2017.
  * Bierkens et al. [2019a] J. Bierkens, P. Fearnhead, G. Roberts, et al.  The zig-zag process and super-efficient sampling for bayesian analysis of big data.  _The Annals of Statistics_ , 47(3):1288–1320, 2019a.
  * Bierkens et al. [2019b] J. Bierkens, G. O. Roberts, and P.-A. Zitt.  Ergodicity of the zigzag process.  _Ann. Appl. Probab._ , 29(4):2266–2301, 08 2019b.  doi: 10.1214/18-AAP1453.  URL https://doi.org/10.1214/18-AAP1453.
  * Bouchard-Côté et al. [2018] A. Bouchard-Côté, S. J. Vollmer, and A. Doucet.  The bouncy particle sampler: A nonreversible rejection-free markov chain monte carlo method.  _Journal of the American Statistical Association_ , 113(522):855–867, 2018.
  * Bradley [2005] R. C. Bradley.  Basic properties of strong mixing conditions. a survey and some open questions.  _Probability Surveys_ , 2:107–144, 2005.
  * Cattiaux et al. [2011] P. Cattiaux, D. Chafai, and A. Guillin.  Central limit theorems for additive functionals of ergodic markov diffusions processes.  _arXiv preprint arXiv:1104.2198_ , 2011.
  * Chakraborty et al. [2019] S. Chakraborty, S. K. Bhattacharya, and K. Khare.  Estimating accuracy of the mcmc variance estimator: a central limit theorem for batch means estimators.  _arXiv preprint arXiv:1911.00915_ , 2019.
  * Chien et al. [1997] C. Chien, D. Goldsman, and B. Melamed.  Large-sample results for batch means.  _Management Science_ , 43(9):1288–1295, 1997\.
  * Chien [1988] C.-H. Chien.  Small-sample theory for steady state confidence intervals.  In _Proceedings of the 20th conference on Winter simulation_ , pages 408–413, 1988.
  * Csáki and Csörgő [1995] E. Csáki and M. Csörgő.  On additive functionals of markov chains.  _Journal of Theoretical Probability_ , 8(4):905–919, 1995.
  * Csörgö and Horváth [1993] M. Csörgö and L. Horváth.  _Weighted approximations in probability and statistics_.  J. Wiley & Sons, 1993.
  * Csörgö and Révész [1979] M. Csörgö and P. Révész.  How big are the increments of a wiener process?  _The Annals of Probability_ , pages 731–737, 1979.
  * Csörgö and Révész [2014] M. Csörgö and P. Révész.  _Strong approximations in probability and statistics_.  Academic press, 2014.
  * Csörgö and Hall [1984] S. Csörgö and P. Hall.  The komlós-major-tusnády approximations and their applications.  _Australian Journal of Statistics_ , 26(2):189–218, 1984.
  * Cuny et al. [2018] C. Cuny, J. Dedecker, and F. Merlevède.  On the komlós, major and tusnády strong approximation for some classes of random iterates.  _Stochastic Processes and their Applications_ , 128(4):1347–1385, 2018.
  * Damerdji [1991] H. Damerdji.  Strong consistency and other properties of the spectral variance estimator.  _Management Science_ , 37(11):1424–1440, 1991\.
  * Damerdji [1994] H. Damerdji.  Strong consistency of the variance estimator in steady-state simulation output analysis.  _Mathematics of Operations Research_ , 19(2):494–512, 1994.
  * Damerdji [1995] H. Damerdji.  Mean-square consistency of the variance estimator in steady-state simulation output analysis.  _Operations Research_ , 43(2):282–291, 1995.
  * Davydov [1968] Y. A. Davydov.  Convergence of distributions generated by stationary stochastic processes.  _Theory of Probability & Its Applications_, 13(4):691–696, 1968.
  * Deheuvels [2000] P. Deheuvels.  Uniform limit laws for kernel density estimators on possibly unbounded intervals.  In _Recent advances in reliability theory_ , pages 477–492. Springer, 2000.
  * Deligiannidis et al. [2019] G. Deligiannidis, A. Bouchard-Côté, A. Doucet, et al.  Exponential ergodicity of the bouncy particle sampler.  _The Annals of Statistics_ , 47(3):1268–1287, 2019.
  * Douc et al. [2018] R. Douc, E. Moulines, P. Priouret, and P. Soulier.  _Markov chains_.  Springer, 2018.
  * Down et al. [1995] D. Down, S. P. Meyn, and R. L. Tweedie.  Exponential and uniform ergodicity of markov processes.  _The Annals of Probability_ , pages 1671–1691, 1995.
  * Duncan et al. [2016] A. B. Duncan, T. Lelievre, and G. Pavliotis.  Variance reduction using nonreversible langevin samplers.  _Journal of statistical physics_ , 163(3):457–491, 2016.
  * Durmus et al. [2020] A. Durmus, A. Guillin, P. Monmarché, et al.  Geometric ergodicity of the bouncy particle sampler.  _Annals of Applied Probability_ , 30(5):2069–2098, 2020.
  * El-Nouty [1999] C. El-Nouty.  On the large increments of fractional brownian motion.  _Statistics & probability letters_, 41(2):169–178, 1999.
  * Flegal et al. [2010] J. M. Flegal, G. L. Jones, et al.  Batch means and spectral variance estimators in markov chain monte carlo.  _The Annals of Statistics_ , 38(2):1034–1070, 2010.
  * Fort and Roberts [2005] G. Fort and G. O. Roberts.  Subgeometric ergodicity of strong markov processes.  _The Annals of Applied Probability_ , 15(2):1565–1589, 2005.
  * Glynn and Whitt [1991] P. W. Glynn and W. Whitt.  Estimating the asymptotic variance with batch means.  _Operations Research Letters_ , 10(8):431–435, 1991.
  * Glynn and Whitt [1992] P. W. Glynn and W. Whitt.  The asymptotic validity of sequential stopping rules for stochastic simulations.  _The Annals of Applied Probability_ , 2(1):180–198, 1992.
  * Goldsman et al. [1990] D. Goldsman, M. Meketon, and L. Schruben.  Properties of standardized time series weighted area variance estimators.  _Management Science_ , 36(5):602–612, 1990.
  * Gong and Flegal [2016] L. Gong and J. M. Flegal.  A practical sequential stopping rule for high-dimensional markov chain monte carlo.  _Journal of Computational and Graphical Statistics_ , 25(3):684–700, 2016.
  * Gut [2009] A. Gut.  _Stopped random walks_.  Springer, 2009.
  * Heunis [2003] A. J. Heunis.  Strong invariance principle for singular diffusions.  _Stochastic processes and their applications_ , 104(1):57–80, 2003.
  * Hobert et al. [2002] J. P. Hobert, G. L. Jones, B. Presnell, and J. S. Rosenthal.  On the applicability of regenerative simulation in markov chain monte carlo.  _Biometrika_ , 89(4):731–743, 2002.
  * Höpfner and Löcherbach [2003] R. Höpfner and E. Löcherbach.  _Limit theorems for null recurrent Markov processes_.  Number 768. American Mathematical Soc., 2003.
  * Hwang et al. [1993] C.-R. Hwang, S.-Y. Hwang-Ma, and S.-J. Sheu.  Accelerating gaussian diffusions.  _The Annals of Applied Probability_ , pages 897–913, 1993.
  * Jones et al. [2006] G. L. Jones, M. Haran, B. S. Caffo, and R. Neath.  Fixed-width output analysis for markov chain monte carlo.  _Journal of the American Statistical Association_ , 101(476):1537–1547, 2006.
  * Kallenberg [1997] O. Kallenberg.  _Foundations of modern probability_ , volume 2.  Springer, 1997.
  * Komlós et al. [1976] J. Komlós, P. Major, and G. Tusnády.  An approximation of partial sums of independent rv’s, and the sample df. ii.  _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , 34(1):33–58, 1976.
  * [47] M. Komlos.  Tusnady (1975) an approximation of partial sums of rv’s, and the sample df, i, z.  _Wahrscheinlichkeitstheorie und Verw. Gebiete_ , 32:111–131.
  * Kuelbs and Philipp [1980] J. Kuelbs and W. Philipp.  Almost sure invariance principles for partial sums of mixing b-valued random variables.  _The Annals of Probability_ , pages 1003–1036, 1980.
  * Lelievre et al. [2013] T. Lelievre, F. Nier, and G. A. Pavliotis.  Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion.  _Journal of Statistical Physics_ , 152(2):237–274, 2013.
  * Li [1992] W. V. Li.  Limit theorems for the square integral of brownian motion and its increments.  _Stochastic processes and their applications_ , 41(2):223–239, 1992.
  * Liu et al. [2021] Y. Liu, D. Vats, and J. M. Flegal.  Batch size selection for variance estimators in mcmc.  _Methodology and Computing in Applied Probability_ , pages 1–29, 2021\.
  * Löcherbach and Loukianova [2008] E. Löcherbach and D. Loukianova.  On nummelin splitting for continuous time harris recurrent markov processes and application to kernel estimation for multi-dimensional diffusions.  _Stochastic Processes and their Applications_ , 118(8):1301–1321, 2008.
  * Löcherbach and Loukianova [2009] E. Löcherbach and D. Loukianova.  The law of iterated logarithm for additive functionals and martingale additive functionals of harris recurrent markov processes.  _Stochastic processes and their applications_ , 119(7):2312–2335, 2009.
  * Major [1979] P. Major.  An improvement of strassen’s invariance principle.  _The Annals of Probability_ , 7(1):55–61, 1979\.
  * Merlevède et al. [2015] F. Merlevède, E. Rio, et al.  Strong approximation for additive functionals of geometrically ergodic markov chains.  _Electronic Journal of Probability_ , 20, 2015.
  * Meyn and Tweedie [1993] S. P. Meyn and R. L. Tweedie.  Stability of markovian processes ii: Continuous-time processes and sampled chains.  _Advances in Applied Probability_ , pages 487–517, 1993.
  * Meyn and Tweedie [2012] S. P. Meyn and R. L. Tweedie.  _Markov chains and stochastic stability_.  Springer Science & Business Media, 2012.
  * Mihalache [2012] S. Mihalache.  Strong approximations and sequential change-point analysis for diffusion processes.  _Statistics & Probability Letters_, 82(3):464–472, 2012.
  * Mykland et al. [1995] P. Mykland, L. Tierney, and B. Yu.  Regeneration in markov chain samplers.  _Journal of the American Statistical Association_ , 90(429):233–241, 1995.
  * Nummelin [1978] E. Nummelin.  A splitting technique for harris recurrent markov chains.  _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , 43(4):309–318, 1978.
  * Nummelin [2004] E. Nummelin.  _General irreducible Markov chains and non-negative operators_.  Number 83. Cambridge University Press, 2004.
  * Nummelin and Tuominen [1982] E. Nummelin and P. Tuominen.  Geometric ergodicity of harris recurrent marcov chains with applications to renewal theory.  _Stochastic Processes and Their Applications_ , 12(2):187–202, 1982.
  * Nummelin and Tuominen [1983] E. Nummelin and P. Tuominen.  The rate of convergence in orey’s theorem for harris recurrent markov chains with applications to renewal theory.  _Stochastic Processes and Their Applications_ , 15(3):295–311, 1983.
  * Ortega [1984] J. Ortega.  On the size of the increments of nonstationary gaussian processes.  _Stochastic processes and their applications_ , 18(1):47–56, 1984.
  * Parzen [1979] E. Parzen.  Nonparametric statistical data modeling.  _Journal of the American statistical association_ , 74(365):105–121, 1979.
  * Philipp and Stout [1975] W. Philipp and W. Stout.  _Almost sure invariance principles for partial sums of weakly dependent random variables_ , volume 161.  American Mathematical Soc., 1975.
  * Révész [1982] P. Révész.  On the increments of wiener and related processes.  _The Annals of Probability_ , pages 613–622, 1982.
  * Revuz [2008] D. Revuz.  _Markov chains_.  Elsevier, 2008.
  * Rey-Bellet and Spiliopoulos [2015] L. Rey-Bellet and K. Spiliopoulos.  Irreversible langevin samplers and variance reduction: a large deviations approach.  _Nonlinearity_ , 28(7):2081, 2015.
  * Rio [1993] E. Rio.  Covariance inequalities for strongly mixing processes.  In _Annales de l’IHP Probabilités et statistiques_ , volume 29, pages 587–597, 1993.
  * Roberts and Rosenthal [2001] G. O. Roberts and J. S. Rosenthal.  Markov chains and de-initializing processes.  _Scandinavian Journal of Statistics_ , 28(3):489–504, 2001.
  * Rogers and Williams [2000] L. C. G. Rogers and D. Williams.  _Diffusions, Markov processes and martingales: Volume 2, Itô calculus_ , volume 2.  Cambridge university press, 2000.
  * Sherman and Goldsman [2002] M. Sherman and D. Goldsman.  Large-sample normality of the batch-means variance estimator.  _Operations Research Letters_ , 30(5):319–326, 2002.
  * Shorack and Wellner [2009] G. R. Shorack and J. A. Wellner.  _Empirical processes with applications to statistics_.  SIAM, 2009.
  * Sigman [1990] K. Sigman.  One-dependent regenerative processes and queues in continuous time.  _Mathematics of Operations Research_ , 15(1):175–189, 1990.
  * Song and Schmeiser [1995] W. T. Song and B. W. Schmeiser.  Optimal mean-squared-error batch sizes.  _Management Science_ , 41(1):110–123, 1995.
  * Stramer and Tweedie [1999] O. Stramer and R. Tweedie.  Langevin-type models i: Diffusions with given stationary distributions and their discretizations.  _Methodology and Computing in Applied Probability_ , 1(3):283–306, 1999.
  * Stroock and Varadhan [1997] D. W. Stroock and S. S. Varadhan.  _Multidimensional diffusion processes_ , volume 233.  Springer Science & Business Media, 1997.
  * Tweedie [1993] S. Tweedie.  Generalized resolvents and harris recurrence.  _Doeblin and modern probability_ , 149, 1993.
  * Vats et al. [2018] D. Vats, J. M. Flegal, G. L. Jones, et al.  Strong consistency of multivariate spectral variance estimators in markov chain monte carlo.  _Bernoulli_ , 24(3):1860–1909, 2018.
  * Vats et al. [2019] D. Vats, J. M. Flegal, and G. L. Jones.  Multivariate output analysis for markov chain monte carlo.  _Biometrika_ , 106(2):321–337, 2019.
  * Williams [2006] D. Williams.  _Stochastic integrals: proceedings of the LMS Durham Symposium, July 7-17, 1980_ , volume 851.  Springer, 2006.