$\displaystyle\left|\frac{g^{\prime}(\widetilde{\xi}\,)}{g^{\prime}(\widetilde{\eta})}-1\right|=\frac{|g^{\prime}(\widetilde{\xi}\,)-g^{\prime}(\widetilde{\eta})|}{|g^{\prime}(\widetilde{\eta})|}=\frac{|g^{\prime\prime}(\widetilde{\chi})|}{|g^{\prime}(\widetilde{\eta})|}|\widetilde{\xi}-\widetilde{\eta}|$

for some $\displaystyle\widetilde{\chi}$ between
$\displaystyle\widetilde{\xi}$ and $\displaystyle\widetilde{\eta}$ (and
similarly, for $\displaystyle\sim$ replaced with $\displaystyle\wedge$)
because $\displaystyle g$ extends to a function in $\displaystyle
C^{2}[z_{i},z_{i+1}]$. So, there exists $\displaystyle\widetilde{c}>1$
(independent of $\displaystyle z_{i}$s) such that

$\displaystyle
1-\widetilde{c}\cdot\theta<\frac{1}{2}\left(\frac{g^{\prime}(\widetilde{\xi}\,)}{g^{\prime}(\widetilde{\eta})}+\frac{g^{\prime}(\widehat{\xi}\,)}{g^{\prime}(\widehat{\eta})}\right)<1+\widetilde{c}\cdot\theta.$

Therefore,

$\displaystyle
1-\widetilde{c}\cdot\theta<\frac{d(\widehat{\zeta})-d(\widetilde{\zeta})}{\widehat{\zeta}-\widetilde{\zeta}}<1+\widetilde{c}\cdot\theta$

where $\displaystyle\widetilde{c}>0$ works for all removed intervals in
$\displaystyle[0,1]$.

We have to maintain that

$\displaystyle(1+c\sqrt{\theta})^{-1}<d^{\prime}(x)<(1-c\sqrt{\theta})^{-1},\,\,x\in[\widetilde{\zeta},\widehat{\zeta}].$

in order to guarantee the second inequality in
$\displaystyle\eqref{eq:DiffeoBds}$. By choosing $\displaystyle d^{\prime}$ to
be close to the constant value $\displaystyle(1+c\sqrt{\theta})^{-1}$ or
$\displaystyle(1-c\sqrt{\theta})^{-1}$ on most of
$\displaystyle[\widetilde{\zeta},\widehat{\zeta}]$, we can make sure the
function

$\displaystyle
d^{\prime}\mapsto\text{AVG}(d^{\prime})\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle\displaystyle:$}}}{\vbox{\hbox{$\displaystyle\textstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptscriptstyle:$}}}{=}}\frac{1}{\widehat{\zeta}-\widetilde{\zeta}}\int_{\widetilde{\zeta}}^{\widehat{\zeta}}d^{\prime}(x)\,dx$

takes values as close as we want to $\displaystyle(1+c\sqrt{\theta})^{-1}$ or
$\displaystyle(1-c\sqrt{\theta})^{-1}$, respectively. Also, by choosing
$\displaystyle c$ sufficiently large,

$\displaystyle(1+c\sqrt{\theta})^{-1}<1-\widetilde{c}\cdot\theta<1+\widetilde{c}\cdot\theta<(1-c\sqrt{\theta})^{-1}.$

Since,
$\displaystyle\text{AVG}(\cdot):C[\widetilde{\zeta},\widehat{\zeta}]\to\mathbb{R}$
is continuous, the intermediate value theorem [69, Theorem 4.22] implies that
there is a continuous function $\displaystyle d^{\prime}$ such that (5.4) is
satisfied. As the final step, we have to check whether the first inequality in
(5.3) holds. To this end, note that

$\displaystyle\displaystyle
d(x)-x=\int_{\widetilde{\zeta}}^{x}d^{\prime}(y)\,dy+(d(\widetilde{\zeta})-\widetilde{\zeta})+(\widetilde{\zeta}-x),$

and hence,

$\displaystyle\displaystyle|d(x)-x|$
$\displaystyle\displaystyle\leq(1-c\sqrt{\theta})^{-1}|\widetilde{\zeta}-x|+|d(\widetilde{\zeta})-\widetilde{\zeta})|+|\widetilde{\zeta}-x|$
$\displaystyle\displaystyle\leq(1-c\sqrt{\theta})^{-1}\widetilde{\theta}^{2}+\theta+\widetilde{\theta}^{2}<c\sqrt{\theta}$

for $\displaystyle\theta$ sufficiently small. Therefore, $\displaystyle
d^{\prime}$ extends to a continuous and strictly increasing function on
$\displaystyle[0,1]$ in such a way that $\displaystyle d$ is continuous and
satisfies (5.3). Since $\displaystyle d$ is a diffeomorphism such that

$\displaystyle g|_{B_{\theta}}=\widetilde{g}\circ
d|_{B_{\theta}},\,m(B_{\theta})>1-c\sqrt{\theta},\,|d(x)-x|<c\sqrt{\theta},\,\Big{|}\frac{1}{d^{\prime}(x)}-1\Big{|}<c\sqrt{\theta},$

we have $\displaystyle\mathfrak{d}(g,\mathfrak{L}_{\theta})\leq
c\cdot\sqrt{\theta}.$ Hence, taking $\displaystyle\mathcal{L}_{\theta}$ and
$\displaystyle\mathcal{L}_{0}$ to be the transfer operators of
$\displaystyle\mathfrak{L}_{\theta}$ and $\displaystyle g$, respectively, we
have from (5.2),

(5.5)
$\|\mathcal{L}_{\theta}-\mathcal{L}_{0}\|_{\mathcal{B},\widetilde{\mathcal{B}}}\leq
12c\cdot\sqrt{\theta}.$

As a result, we have the continuity at $\displaystyle(0,0)$:

$\displaystyle\displaystyle\|\mathcal{L}_{\theta,is}-\mathcal{L}_{0}\|_{{\operatorname{BV}},L^{1}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta,is}-\mathcal{L}_{\theta}\|_{{\operatorname{BV}},L^{1}}+\|\mathcal{L}_{\theta}-\mathcal{L}_{0}\|_{{\operatorname{BV}},L^{1}}\to
0,\,\,\,\text{as}\,\,\,(\theta,s)\to(0,0).$

In fact, we have continuity in a neighbourhood of $\displaystyle(0,0)$:

$\displaystyle\displaystyle\|$
$\displaystyle\displaystyle\mathcal{L}_{\theta,is}\psi-\mathcal{L}_{\bar{\theta},i\bar{s}}\psi\|_{L^{1}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta,is}\psi-\mathcal{L}_{\bar{\theta},is}\psi\|_{L^{1}}+\|\mathcal{L}_{\bar{\theta},is}\psi-\mathcal{L}_{\bar{\theta},i\bar{s}}\psi\|_{L^{1}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta}(e^{ish_{\theta}}\psi)-\mathcal{L}_{\bar{\theta}}(e^{ish_{\bar{\theta}}}\psi)\|_{L^{1}}+\|\mathcal{L}_{\bar{\theta},is}\psi-\mathcal{L}_{\bar{\theta},i\bar{s}}\psi\|_{{\operatorname{BV}}}$
$\displaystyle\displaystyle\leq\|(\mathcal{L}_{\theta}-\mathcal{L}_{\bar{\theta}})(e^{ish_{\theta}}\psi)\|_{L^{1}}+\|\mathcal{L}_{\bar{\theta}}((e^{ish_{\theta}}-e^{ish_{\bar{\theta}}})\psi)\|_{L^{1}}+\|(\mathcal{L}_{\bar{\theta},is}-\mathcal{L}_{\bar{\theta},i\bar{s}})\psi\|_{{\operatorname{BV}}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta}-\mathcal{L}_{\bar{\theta}}\|_{{\operatorname{BV}},L^{1}}\|e^{ish_{\theta}}\psi\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta}}\|_{L^{1},L^{1}}\|\|(e^{ish_{\theta}}-e^{ish_{\bar{\theta}}})\psi\|_{L^{1}}+\|\mathcal{L}_{\bar{\theta},is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}$
$\displaystyle\displaystyle\leq\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})|s|\|h_{\theta}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta}}\|_{L^{1},L^{1}}|s|\|(h_{\theta}-h_{\bar{\theta}})\psi\|_{L^{1}}+\|\mathcal{L}_{\bar{\theta},is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}$
$\displaystyle\displaystyle\leq\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})|s|\|h_{\theta}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta}}\|_{L^{1},L^{1}}|s|\|h_{\theta}-h_{\bar{\theta}}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta},is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}}}\|\psi\|_{{\operatorname{BV}}}$

Therefore,

$\displaystyle\displaystyle\|\mathcal{L}_{\theta,is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}},L^{1}}$
$\displaystyle\displaystyle=\sup_{\|\psi\|_{BV}\leq
1}\|\mathcal{L}_{\theta,is}\psi-\mathcal{L}_{\bar{\theta},i\bar{s}}\psi\|_{L^{1}}$
(5.6)
$\displaystyle\displaystyle\leq\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})|s|\|h_{\theta}\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta}}\|_{L^{1},L^{1}}|s|\|h_{\theta}-h_{\bar{\theta}}\|_{{\operatorname{BV}}}+\|\mathcal{L}_{\bar{\theta},is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}}}$

Assume $\displaystyle\theta>\bar{\theta}$. Then,
$\displaystyle\mathfrak{L}_{\bar{\theta}}$ is a piece-wise expanding map and
$\displaystyle\mathfrak{L}_{\theta}$ is a spline approximation of
$\displaystyle\mathfrak{L}_{\bar{\theta}}$. So, from the general analysis we
did, we know that
$\displaystyle\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})\to
0$ as $\displaystyle\theta\to\bar{\theta}^{+}$. When
$\displaystyle\bar{\theta}>\theta$, $\displaystyle\mathfrak{L}_{\bar{\theta}}$
is a spline approximation of the piecewise expanding map
$\displaystyle\mathfrak{L}_{\theta}$. So, when $\displaystyle\bar{\theta}$ and
$\displaystyle\theta$ are close,
$\displaystyle\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})$
is close to $\displaystyle 0$ as well. So,
$\displaystyle\mathfrak{d}(\mathfrak{L}_{\theta},\mathfrak{L}_{\bar{\theta}})\to
0$ as $\displaystyle\theta\to\bar{\theta}^{-}$. So, from (5.2), we have,

$\displaystyle\lim_{(\theta,s)\to(\bar{\theta},\bar{s})}\|\mathcal{L}_{\theta,is}-\mathcal{L}_{\bar{\theta},i\bar{s}}\|_{{\operatorname{BV}},L^{1}}=0$

(B) with
$\displaystyle\mathcal{X}_{0}=\mathcal{X}^{(+)}_{3}={\operatorname{BV}}([0,1])$
and any $\displaystyle p_{0}\geq 1$, and (C) follows as in the previous
example from the results in [10, 22]. We essentially use the non-arithmeticity
of $\displaystyle h_{\theta}$ here. However,
$\displaystyle\mathcal{L}_{\theta}$ need not be a covering for [10, Lemma 1]
to hold.

###  5.3. Markov models

In addition to dynamical systems, our continuous Edgeworth expansion result in
4.2 is also applicable in the Markovian setting. Several ideas in [23, 63, 33,
22] give us conditions to establish first-order Edgeworth expansions under the
optimal moment conditions of the iid case. In the discussion below, using
Markov integral operators in place of transfer operators, we show that first-
order continuous Edgeworth expansions hold for $\displaystyle V-$geometrically
ergodic Markov processes, and hence, the bootstrap has the desired asymptotic
accuracy in that setting.

Let $\displaystyle(X,\mathcal{F},m)$ be a Borel probability space and
$\displaystyle J\subset\mathbb{R}$ be a neighbourhood of $\displaystyle 0$.
Recall that $\displaystyle\mathcal{M}_{1}(X)$ is the set of Borel probability
measures on $\displaystyle X$. For each $\displaystyle\theta\in J$, let
$\displaystyle\\{x^{\theta}_{n}\\}_{n\geq 0}$ be a time homogeneous Markov
chain on $\displaystyle X$ with initial measure
$\displaystyle\mu_{\theta}\in\mathcal{M}_{1}(X)$ and the Markov transition
operator

$\displaystyle\mathcal{L}_{\theta}(\varphi)(x)=\int_{X}\varphi(y)d\mathbf{P}_{\theta}(x,dy)={\mathbb{E}}(\varphi|x_{0}=x),\,\,\,\varphi\in
L^{1}(X).$

####  5.3.1. $\displaystyle V-$geometrically ergodic Markov chains

Let $\displaystyle V:X\to[1,\infty)$ be a measurable function with
$\displaystyle{\mathbb{E}}_{m}(V)<\infty$. Let
$\displaystyle\|f\|_{L^{\infty}_{V}}=\sup_{x\in X}|f(x)|/V(x)$ for
$\displaystyle f:X\to\mathbb{C}$, let

$\displaystyle
L^{\infty}_{V}\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle\displaystyle:$}}}{\vbox{\hbox{$\displaystyle\textstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptscriptstyle:$}}}{=}}\left\\{f:X\to\mathbb{C}\Big{|}\|f\|_{L^{\infty}_{V}}<\infty\right\\},$

and let $\displaystyle\|\cdot\|_{L^{\infty}_{V},L^{\infty}_{V}}$ be the
operator norm. Now suppose $\displaystyle\\{x^{0}_{n}\\}$ is irreducible and
aperiodic. In addition, suppose $\displaystyle\\{x^{0}_{n}\\}$ is
$\displaystyle V-$geometrically ergodic, i.e., there exist
$\displaystyle\nu\in\mathcal{M}_{1}(X)$, $\displaystyle C>0$ and
$\displaystyle\kappa\in(0,1)$ such that
$\displaystyle{\mathbb{E}}_{\nu}(V)<\infty$ and

(5.7) $\|\mathcal{L}^{n}-{\mathbb{E}}_{\nu}[\,\cdot\,]{\bf
1}_{X}\|_{L^{\infty}_{V},L^{\infty}_{V}}\leq C\kappa^{n}.$

As in [23], we assume that

> (i) there exist $\displaystyle N,L>0$ and $\displaystyle\kappa\in(0,1)$ such
> that $\displaystyle\mathcal{L}^{N}_{\theta}V\leq\kappa^{N}V+L{\bf
> 1}_{X},\,\,\theta\in J$, (ii) $\displaystyle\lim_{\theta\to
> 0}\|\mathcal{L}_{\theta}-\mathcal{L}\|_{L^{\infty},L^{\infty}_{V}}=0.$

Let $\displaystyle h_{\theta}\in L^{\infty}_{V},\,\theta\in J$ be such that
$\displaystyle h_{\theta}\to h_{0}$ in $\displaystyle L^{\infty}_{V}$ as
$\displaystyle\theta\to 0$ (often, this is vacuous in practice because we take
$\displaystyle h_{\theta}=h_{0}$ for all $\displaystyle\theta$) and
$\displaystyle h_{\theta}\in L^{3}$ for all $\displaystyle\theta$. We further
assume as in [63, Section 7] that

(5.8)
$\sup_{a=1,2,3}\sup_{j=0,\dots,a}\|\mathcal{L}_{\theta}(h^{a-j}_{\theta}V^{j/3})\|_{L^{\infty}_{V^{a/3}}}<\infty.$

We note from [63, Remark 7.6] that (5.8) is true as soon as
$\displaystyle\|h^{3}_{\theta}\|_{L^{\infty}_{V}}<\infty$ when the Markov
chains are stationary. Finally, we assume the following non-lattice condition
in [33, p. 435]:

> It is not the case that there exist $\displaystyle a\in\mathbb{R}$ and
> $\displaystyle b>0$, a $\displaystyle\nu-$full
> $\displaystyle\mathcal{L}-$absorbing set $\displaystyle U\in\mathcal{F}$ and
> $\displaystyle\zeta\in L^{\infty}$ such that for all $\displaystyle x\in U,$
> $\displaystyle h(y)+\zeta(y)-\zeta(x)\in a+b\mathbb{Z}.$

######  Remark 5.2:

Assuming a weak form of convergence in the Assumption (ii) as opposed to
convergence in $\displaystyle\|\cdot\|_{L^{\infty}_{V},L^{\infty}_{V}}$,
allows us to for many well-studied family of Markov chains, for example, the
auto-regressive processes in 2 below; see [23] for details.

######  Example 2:

Let $\displaystyle X=\mathbb{R}$ and $\displaystyle\\{x^{\theta}_{n}\\}_{n\geq
0}$ be defined by $\displaystyle x^{\theta}_{n}=\theta\cdot
x_{n-1}+\xi_{n},\,\,\,n\geq 1$ where $\displaystyle x_{0}$ is a
$\displaystyle\mathbb{R}$-valued random variable,
$\displaystyle\theta\in(-1,1)$ and $\displaystyle\xi_{n}$ is a sequence of
absolutely continuous iid random variables that are independent of
$\displaystyle x_{0}$. Then, one can show that
$\displaystyle\\{x^{\theta}_{n}\\}$ is $\displaystyle V-$geometrically ergodic
with $\displaystyle V(x)=1+|x|$, and
$\displaystyle\|\mathcal{L}_{\theta}-\mathcal{L}_{\theta_{0}}\|_{L^{\infty}_{V},L^{\infty}_{V}}\not\to
0$ but
$\displaystyle\|\mathcal{L}_{\theta}-\mathcal{L}_{\theta_{0}}\|_{L^{\infty},L^{\infty}_{V}}\to
0$ as $\displaystyle\theta\to\theta_{0}$.

Now, we verify our spectral assumptions. To this end, let
$\displaystyle\mathcal{B}=\mathcal{X}_{0}=\mathbb{C}\cdot{\bf 1}_{X}$,
$\displaystyle\mathcal{X}^{(+)}_{a}=\mathcal{X}_{a+1}=L^{\infty}_{V^{(a+1)/3}}$
for $\displaystyle a\in\\{0,1,2\\}$ and
$\displaystyle\mathcal{X}^{(+)}_{3}=\mathcal{X}_{4}=\widetilde{\mathcal{B}}=L^{1}$
as in [63, Section 7] with $\displaystyle r=2$ there. Then, it is easy to see
that $\displaystyle\mathcal{X}_{a},a=0,1,2,3,4$ forms a chain of Banach spaces
of the form (3.3) (here we use essentially the assumptions $\displaystyle
V\geq 1$, $\displaystyle{\mathbb{E}}_{m}(V)<\infty$ and
$\displaystyle{\mathbb{E}}_{\nu}(V)<\infty$). Recall that
$\displaystyle\mathcal{L}_{\theta,is}$ given by (2.8). It is straightforward
that $\displaystyle\mathcal{L}_{\theta,is}$ are linear operators on
$\displaystyle\mathcal{B}$. Applying (5.8) with $\displaystyle a=j$, we have
for all $\displaystyle\varphi\in L^{\infty}_{V^{a/3}}$,

$\displaystyle\displaystyle\|\mathcal{L}_{\theta,is}(\varphi)\|_{L^{\infty}_{V^{a/3}}}=\sup_{x}|(V(x))^{-a/3}\mathcal{L}_{\theta,is}(\varphi)(x)|\leq\|\varphi\|_{L^{\infty}_{V^{a/3}}}\sup_{x}|(V(x))^{-a/3}\mathcal{L}_{\theta}(V^{a/3})(x)|\lesssim\|\varphi\|_{L^{\infty}_{V^{a/3}}}.$

So, $\displaystyle\mathcal{L}_{\theta,is}$ are bounded linear operators on
$\displaystyle L^{\infty}_{V^{a/3}}$ for all $\displaystyle a=1,2,3$.

Now, under (i) and (ii), we have the following (see [23, Theorem 1]) with
$\displaystyle\nu_{0}\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle\displaystyle:$}}}{\vbox{\hbox{$\displaystyle\textstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptstyle:$}}}{\vbox{\hbox{$\displaystyle\scriptscriptstyle:$}}}{=}}\nu$:

> (1) For each $\displaystyle\bar{\kappa}\in(\kappa,1)$, there exists
> $\displaystyle{\varepsilon}$ such that for all
> $\displaystyle|\theta|<{\varepsilon}$, $\displaystyle\\{x^{\theta}_{n}\\}$
> has a unique invariant probability measure, $\displaystyle\nu_{\theta}$,
> such that $\displaystyle\nu_{\theta}(V)<\infty$ and
> $\sup_{|\theta|<{\varepsilon}}\|\mathcal{L}^{n}_{\theta}-{\mathbb{E}}_{\nu_{\theta}}[\,\cdot\,]{\bf
> 1}_{X}\|_{L^{\infty}_{V},L^{\infty}_{V}}\leq C\bar{\kappa}^{n},$ (2)
> $\displaystyle\lim_{\theta\to 0}\sup_{\|\varphi\|_{L^{\infty}}\leq
> 1}|\nu_{\theta}(\varphi)-\nu_{0}(\varphi)|=0.$

So, for all $\displaystyle\theta$ sufficiently close to $\displaystyle 0$,
$\displaystyle\\{x^{\theta}_{n}\\}$ is $\displaystyle V-$geometrically ergodic
with the same rate $\displaystyle\bar{\kappa}$. Now, we fix
$\displaystyle\bar{\kappa}$, and this fixes $\displaystyle{\varepsilon}$, and
we reduce the $\displaystyle\theta$ range from $\displaystyle J$ to
$\displaystyle J\cap(-{\varepsilon},{\varepsilon})$.

Next, for all $\displaystyle\varphi\in L^{\infty}$,

$\displaystyle\displaystyle\|\mathcal{L}_{\theta,is}\varphi-\mathcal{L}_{0,i\bar{s}}\varphi\|_{L^{\infty}_{V}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta}-\mathcal{L}\|_{L^{\infty},L^{\infty}_{V}}\|e^{ish_{\theta}}\varphi\|_{L^{\infty}}+\|\mathcal{L}_{0}\|_{L^{\infty}_{V},L^{\infty}_{V}}\|(e^{ish_{\theta}}-e^{i\bar{s}h_{0}})\varphi\|_{L^{\infty}_{V}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}_{\theta}-\mathcal{L}\|_{L^{\infty},L^{\infty}_{V}}\|\varphi\|_{L^{\infty}}+\|\mathcal{L}_{0}\|_{L^{\infty}_{V},L^{\infty}_{V}}\left(\|h_{\theta}-h_{0}\|_{L^{\infty}_{V}}+|s-\bar{s}|\|h_{0}\|_{L^{\infty}_{V}}\right)\|\varphi\|_{L^{\infty}},$

and hence,

$\displaystyle\lim_{(\theta,s)\to(0,\bar{s})}\|\mathcal{L}_{\theta,is}-\mathcal{L}_{0,i\bar{s}}\|_{L^{\infty},L^{\infty}_{V}}=0.$

A similar argument gives that, for all $\displaystyle|\theta|<{\varepsilon}$,

$\displaystyle\lim_{s\to\bar{s}}\|\mathcal{L}_{\theta,is}-\mathcal{L}_{\theta,i\bar{s}}\|_{L^{\infty},L^{\infty}_{V}}=0.$

Due to 3.1, 3.7 and 3.11, this weaker assumption is a sufficient replacement
for (A)(1) where $\displaystyle\mathcal{B}=\mathbb{C}\cdot{\bf 1}_{X}$
equipped with $\displaystyle\|\cdot\|_{L^{\infty}}$ and
$\displaystyle\widetilde{\mathcal{B}}=L^{1}$. (A)(2) follows from (1) because

$\displaystyle\displaystyle\|\mathcal{L}^{n}_{\theta,is}(c\cdot{\bf
1}_{X})\|_{L^{\infty}}\leq|c|\|\mathcal{L}^{n}_{\theta}({\bf
1}_{X})\|_{L^{\infty}}\leq|c|(C\bar{\kappa}^{n}+\|{\bf
1}_{X}\|_{L^{\infty}_{V}})$
$\displaystyle\displaystyle=C\bar{\kappa}^{n}\|c\cdot{\bf
1}_{X}\|_{L^{\infty}}+\|c\cdot{\bf 1}_{X}\|_{L^{\infty}_{V}}$
$\displaystyle\displaystyle\leq C\bar{\kappa}^{n}\|c\cdot{\bf
1}_{X}\|_{L^{\infty}}+{\mathbb{E}}_{m}(V)\|c\cdot{\bf 1}_{X}\|_{L^{1}}$

for all $\displaystyle c\in\mathbb{C}$.

Due to [63, Theorem 7.5], for $\displaystyle a=1,2,3$,

$\displaystyle\|\mathcal{L}_{\theta,is}(\varphi)\|_{L^{\infty}_{V^{a/3}}}\leq
C\kappa^{n}\|\varphi\|_{L^{\infty}_{V^{a/3}}}+\|\varphi\|_{L^{1}},$

for all $\displaystyle\varphi\in L^{\infty}_{V^{a/3}}$. So, we have (B)(3).
The same theorem gives the regularity required in (B)(1,2). From

$\displaystyle\|(h_{\theta})^{j}\varphi\|_{L^{(a+j)/3}}=\|(h_{\theta})^{j}\|_{L^{\infty}_{V^{j/3}}}\|\varphi\|_{L^{\infty}_{V^{a/3}}},$

we note that
$\displaystyle\mathcal{L}_{\theta,0}^{(j)}(\cdot)=\mathcal{L}_{\theta}((ih_{\theta})^{j}\,\cdot)$
is a linear operator from $\displaystyle L^{\infty}_{V^{a/3}}$ to
$\displaystyle L^{\infty}_{V^{(a+j)/3}}$. So, due to (5.8), for all
$\displaystyle\theta$, for all $\displaystyle a=0,1,2$ and for all
$\displaystyle\varphi\in L^{\infty}_{V^{a/3}}$, the Taylor expansion

$\displaystyle\mathcal{L}_{\theta,is}(\varphi)=\sum_{k=0}^{j}\frac{\mathcal{L}_{\theta}((ih_{\theta})^{k}\varphi)}{k!}s^{k}+\|\varphi\|_{L^{\infty}_{V^{a/3}}}o(s^{j}),\,\,\,s\in\mathbb{R}$

holds in $\displaystyle L^{\infty}_{V^{(a+j)/3}}$ for $\displaystyle
j=1,\dots,3-a$ with $\displaystyle L^{\infty}_{V^{0}}$ is understood to be
$\displaystyle\mathbb{C}\cdot{\bf 1_{X}}$.

We recall from [33, Lemma 10.1] that (5.7) and (1) above implies that for all
$\displaystyle\eta\in(0,1]$,

$\|\mathcal{L}_{\theta}^{n}-{\mathbb{E}}_{\nu}[\,\cdot\,]{\bf
1}_{X}\|_{L^{\infty}_{V^{\eta}},L^{\infty}_{V^{\eta}}}\leq C\kappa^{n},$

and we have (C)(1,3). (C)(4) follows from the non-lattice assumption; see [33,
Section 5.2]. (C)(2) is a given because
$\displaystyle\mathbf{P}_{\theta}(x,\cdot)$ are probability measures. In fact,
we have the stronger conclusion that $\displaystyle\mathcal{L}_{\theta}$ has a
spectral gap of $\displaystyle 1-\bar{\kappa}$ on $\displaystyle
L^{\infty}_{V^{\eta}_{\infty}}$ for all $\displaystyle\eta\in(0,1]$ and
$\displaystyle|\theta|<{\varepsilon}$.

######  Remark 5.3:

In the special case of $\displaystyle V\equiv 1$, we obtain the case of
uniformly ergodic Markov chains. We consider a single space
$\displaystyle\mathcal{B}=\widetilde{\mathcal{B}}=L^{\infty}$ and most of our
assumptions at the beginning of this section become vacuous, and we recover
the results of [13] for bounded observables.

######  Remark 5.4:

We could have considered the more general case of
$\displaystyle(Y,\mathcal{G},\widetilde{m})$ being another Borel probability
space, $\displaystyle\\{y_{n}\\}_{n\geq 0}$ being a sequence of iid sequence
of random variables on $\displaystyle Y$ with common distribution
$\displaystyle\widetilde{m}$ and independent of
$\displaystyle\\{x^{\theta}_{n}\\}$ for all $\displaystyle\theta\in J$,
$\displaystyle h_{\theta}:X\times X\times Y\to\mathbb{C}$, $\displaystyle
X^{\theta}_{n}=h_{\theta}(x^{\theta}_{n-1},x^{\theta}_{n},y_{n})$.

Part II $\displaystyle-$ Hall

, we have

$\displaystyle\sup_{z\in\mathbb{R}}\Big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=o_{\text{a.s.}}(n^{-1/2}).$

Our LABEL:thm:PivotedBootstrap and LABEL:cor:PivotedBootstrap generalize this
result where the existence of the first three finite moments follow from (B)
and the distribution being non-lattice follows from (C)(4). However, if we
have stronger assumptions on the common distribution, we can say more. In
particular, if, in the iid setting, the common distribution has four finite
moments and satisfies Cramér’s condition,

$\displaystyle\limsup_{|t|\to\infty}|{\mathbb{E}}(e^{isX})|<1\,,$

then the second order Edgeworth expansion exists, and it follows that,

$\displaystyle\sup_{z\in\mathbb{R}}\Big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=\mathcal{O}_{\text{a.s.}}(n^{-1}).$

See, for example, [28, Section 3.3]. So, in order to guarantee to improve the
asymptotic accuracy of our bootstrap, we need stronger spectral assumptions
that guarantee existence of at least four asymptotic moments and better
control of the characteristic function away from the origin. We plan to
explore this in a separate work.

###  7.2. Nonpivoted bootstrap

Suppose the asymptotic variance $\displaystyle\sigma^{2}$ is unknown and
cannot be estimated easily. Then write

(7.4) $T_{n}=\frac{1}{\sqrt{n}}\left(S_{n}(h)-nA\right).$

Using Algorithm LABEL:alg:nonpivot, we generate the bootstrap samples
$\displaystyle\\{x^{*}_{k}\\}_{0\leq k\leq n}$. The bootstrap estimator of
$\displaystyle T_{n}$ is given by

(7.5) $T^{*}_{n}=\frac{1}{\sqrt{n}}\left(S_{n}^{*}(h)-nA^{*}\right).$

######  Theorem 7.2:

Suppose that, almost surely, there exists $\displaystyle N$ such that
$\displaystyle\\{\widehat{g}\,|\,n\geq N\\}\cup\\{g\\}$ satisfy the
Assumptions (A), (B), (C), (D), and (E). Let $\displaystyle\mathcal{C}$ is a
space in (3.3) such that
$\displaystyle\mathcal{C}\hookrightarrow\mathcal{X}_{3}\hookrightarrow\mathcal{X}^{(+)}_{3}$
and assume that $\displaystyle\|\rho_{\mu^{*}}-\rho_{\mu}\|_{\mathcal{C}}\to
0$ as $\displaystyle n\to\infty$. Then,

$\displaystyle\sup_{z\in\mathbb{R}}\big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=\mathcal{O}_{\text{a.s.}}(|\sigma^{*}-\sigma|)+\mathcal{O}_{\text{a.s.}}(n^{-1/2}).$

###### Proof.

Even though the Birkhoff sum $\displaystyle S_{n}^{*}(h)$ is not normalized by
$\displaystyle\sigma^{*}$, due to uniformity in $\displaystyle x$ of we can
still use the conclusion of 4.2 to obtain

$\displaystyle\displaystyle\sup_{z\in\mathbb{R}}\bigg{|}{\mathbb{P}}(T^{*}_{n}\leq
z\,$
$\displaystyle\displaystyle|\,x_{0},\dots,x_{n-1})-\mathcal{E}_{n}\left(\frac{z}{\sigma^{*}}\right)\bigg{|}$
$\displaystyle\displaystyle=\sup_{z\in\mathbb{R}}\left|{\mathbb{P}}\left(\frac{S_{n}^{*}(h)-nA^{*}}{\sigma^{*}\sqrt{n}}\leq\frac{z}{\sigma^{*}}\,\Big{|}\,x_{0},\dots,x_{n-1}\right)-\mathcal{E}_{n}\left(\frac{z}{\sigma^{*}}\right)\right|$
$\displaystyle\displaystyle\leq\sup_{x\in\mathbb{R}}\left|{\mathbb{P}}\left(\frac{S^{*}_{n}(h)-nA^{*}}{\sigma^{*}\sqrt{n}}\leq
x\,\Big{|}\,\,x_{0},\dots,x_{n-1}\right)-\mathcal{E}_{n}\left(x\right)\right|=o_{\text{a.s.}}(n^{-1/2}).$

Similarly, we have

$\displaystyle\sup_{z\in\mathbb{R}}\left|{\mathbb{P}}_{\mu_{0}}(T_{n}\leq
z)-\mathcal{E}_{n}\left(\frac{z}{\sigma}\right)\right|=o(n^{-1/2}).$

As a result,

$\displaystyle\displaystyle\sup_{z\in\mathbb{R}}\Big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}$ $\displaystyle\displaystyle(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle\sup_{z\in\mathbb{R}}\left|\mathcal{E}_{n}\left(\frac{z}{\sigma}\right)-\mathcal{E}_{n}\left(\frac{z}{\sigma^{*}}\right)\right|+o_{\text{a.s.}}(n^{-1/2})$

Note that

$\displaystyle\displaystyle\mathcal{E}_{n}\left(\frac{z}{\sigma}\right)-\mathcal{E}_{n}\left(\frac{z}{\sigma^{*}}\right)$
$\displaystyle\displaystyle=\frac{1}{\sqrt{2\pi}\sigma^{*}}\int_{-\infty}^{z}e^{-\frac{y^{2}}{2(\sigma^{*})^{2}}}\,dy-\frac{1}{\sqrt{2\pi}\sigma}\int_{-\infty}^{z}e^{-\frac{y^{2}}{2\sigma^{2}}}\,dy$
$\displaystyle\displaystyle\phantom{aaaaaaaaaaaaaaa}+\frac{1}{\sqrt{n}}\left(P\left(\frac{x}{\sigma^{*}}\right)\mathfrak{n}\left(\frac{x}{\sigma^{*}}\right)-P\left(\frac{x}{\sigma}\right)\mathfrak{n}\left(\frac{x}{\sigma}\right)\right).$

Since the first term above is
$\displaystyle\mathcal{O}_{\text{a.s.}}(|\sigma^{*}-\sigma|)$ and the second
is $\displaystyle\mathcal{O}_{\text{a.s.}}(n^{-1/2})$, together we have

$\displaystyle\mathcal{O}_{\text{a.s.}}(|\sigma^{*}-\sigma|)+\mathcal{O}_{\text{a.s.}}(n^{-1/2}),$

which dominates the $\displaystyle o_{\text{a.s.}}(n^{-1/2})$ term.
$\displaystyle\square$

######  Remark 7.4:

We recall that the non-pivoted bootstrap does not depend on the knowledge of
$\displaystyle\sigma$ while the Gaussian approximation does (in fact, it is
impossible without $\displaystyle\sigma$). However, by 7.2, provided that
$\displaystyle|\sigma^{*}-\sigma|\lesssim n^{-1/2}$, the non-pivoted bootstrap
achieves an asymptotic accuracy of
$\displaystyle\mathcal{O}_{\text{a.s.}}(n^{-1/2})$, which is equal to that of
the Gaussian approximation. Hence, in a somewhat oracular way, the non-pivoted
bootstrap “knows” $\displaystyle\sigma$ without paying any price.

In general, we only know that
$\displaystyle|\sigma^{*}-\sigma|=o_{\text{a.s}}(1)$. In this case, by 3.12,
we have the following corollary.

######  Corollary 7.2:

Suppose that, almost surely, the Assumptions (A), (B), (C), and (D) hold with
$\displaystyle p_{0}\geq 2$, and $\displaystyle h$ is not $\displaystyle
g-$cohomologous to a constant in $\displaystyle L^{2}(m)$. Let
$\displaystyle\mathcal{C}$ is a space in (3.3) such that
$\displaystyle\mathcal{C}\hookrightarrow\mathcal{X}_{3}\hookrightarrow\mathcal{X}^{(+)}_{3}$
and assume that $\displaystyle\|\rho_{\mu^{*}}-\rho_{\mu}\|_{\mathcal{C}}\to
0$ as $\displaystyle n\to\infty$. Then,

$\displaystyle\sup_{z\in\mathbb{R}}\big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=o_{\text{a.s.}}(1).$

##  8\. Application of the Bootstrap to our Examples

In this section, we describe some of the standard techniques available to
apply the Algorithms LABEL:alg:pivot and LABEL:alg:nonpivot in the context of
dynamical systems.

###  8.1. Choice of $\displaystyle\widehat{g}$ and $\displaystyle\mu^{*}$

To approximate $\displaystyle g$, there are a two key methods: Either one can
use (8), (9), (12), and (13) of [59], or one can use standard spline
approximations; which method is more appropriate depends on the context.

For the choice of $\displaystyle\mu^{*}$, one may use the kernel estimator of
$\displaystyle\nu_{0}$ on [30, page 22], defined by

$\displaystyle\widehat{\rho}_{\nu_{0}}(x)=\frac{1}{n\mathfrak{b}}\sum_{k=1}^{n}K\Big{(}\frac{x-x_{k}}{\mathfrak{b}}\Big{)},$

where $\displaystyle K:\mathbb{R}\to\mathbb{R}$ is the density function for a
standard Gaussian random variable and $\displaystyle\mathfrak{b}$ is a
bandwidth. Alternatively, one can let $\displaystyle\mu^{*}$ (having density
$\displaystyle\rho_{\mu^{*}}$) be a fixed distribution that depends on neither
$\displaystyle x_{0},\dots,x_{n-1}$ nor $\displaystyle n$.

###  8.2. Suitability of $\displaystyle\widehat{g}$ and $\displaystyle\mu^{*}$

Now, we discuss the suitability of the choices of $\displaystyle\widehat{g}$
and $\displaystyle\mu^{*}$ to align with our abstract setting in
LABEL:part:EdgeExp, and hence, ensure the asymptotic accuracy of the
bootstrap. We require the following two conditions to be satisfied.

  * –

Almost surely, there is $\displaystyle N$ such that
$\displaystyle\\{\widehat{g}\,|\,n\geq N\\}\cup\\{g\\}$ satisfy the
Assumptions (A), (B), and (C).

  * –

There exists a space $\displaystyle\mathcal{C}$ in (3.3) such that
$\displaystyle\mathcal{C}\hookrightarrow\mathcal{X}_{3}\hookrightarrow\mathcal{X}^{(+)}_{3}$
and $\displaystyle\|\rho_{\mu^{*}}-\rho_{\mu}\|_{\mathcal{C}}\to 0$ as
$\displaystyle n\to\infty$.

In both smooth expanding and peicewise expanding maps, the (C)(1) holds. So,
$\displaystyle g$ possesses a unique acip $\displaystyle\nu_{0}$ (with density
$\displaystyle\rho_{\nu_{0}}$); moreover, $\displaystyle(g,\nu_{0})$ is strong
mixing, and hence, ergodic. As a result, $\displaystyle\nu_{0}-$almost every
orbit of $\displaystyle g$ is dense in the support of $\displaystyle\nu_{0}$;
cf. [77, p. 29]. So, $\displaystyle\nu_{0}-$almost surely, $\displaystyle
x_{0},\dots,x_{n-1},\dots$ is a dense orbit. Therefore, we can construct
spline approximations $\displaystyle\widehat{g}$. We have already shown that
spline approximations of expanding maps satisfy the Assumptions (A), (B), and
(C); see Sections 5.1 and 5.2, and hence, they are an ideal choice for our
simulations.

When a kernel estimators that changes with the sample size $\displaystyle n$,
the above suitability condition for stationary processes can be checked using
ideas in [80, Section 3] where the convergence properties of of
$\displaystyle\widehat{\rho}_{\nu_{0}}$ to $\displaystyle\rho_{\nu_{0}}$ are
discussed. In particular, if $\displaystyle\\{X_{n}\\}$ is stationary and
$\displaystyle\beta-$mixing, and $\displaystyle\rho_{\nu_{0}}$ is sufficiently
regular, then $\displaystyle\widehat{\rho}_{\nu_{0}}$ and
$\displaystyle\widehat{\rho}^{\prime}_{\nu_{0}}$ converge to
$\displaystyle\rho_{\nu_{0}}$ and $\displaystyle\rho^{\prime}_{\nu_{0}}$,
respectively and the convergence is uniform. In particular, for $\displaystyle
V-$geometrically ergodic Markov chains, which by definition are
$\displaystyle\beta-$mixing, we have
$\displaystyle\|\widehat{\rho}_{\nu_{0}}-\rho_{\nu_{0}}\|_{L^{\infty}}\to 0$.

Since the kernel estimators in dynamical systems and
$\displaystyle\beta-$mixing setting share many common properties [30], we
conjecture that, in the dynamical systems examples in Section 5, we also have
the uniform convergence of $\displaystyle\widehat{\rho}_{\nu_{0}}$ and
$\displaystyle\widehat{\rho}^{\prime}_{\nu_{0}}$ to
$\displaystyle\rho_{\nu_{0}}$ and $\displaystyle\rho^{\prime}_{\nu_{0}}$,
which implies
$\displaystyle\|\widehat{\rho}_{\nu_{0}}-\rho_{\nu_{0}}\|_{{\operatorname{BV}}}\to
0$. However, this is still an open problem. In general, it is interesting to
study the convergence of derivatives of kernel estimators of acips of
dynamical systems but this has not been pursued in the previous literature.

When the bootstrap initial measure, $\displaystyle\mu^{*}$, is fixed, the only
condition it should satisfy is being absolutely continuous, and we opt for
this option.

###  8.3. Improved accuracy of the non-pivoted bootstrap

Now, we explain why we our non-pivoted bootstrap algorithm exhibits better
asymptotic accuracy than $\displaystyle o_{\text{a.s.}}(1)$ when applied to
our examples in Section 5.1 and Section 5.2.

######  Example 3 (Piecewise $\displaystyle C^{2}$ expanding maps, a
continuation of Section 5.2):

Since we assume $\displaystyle x_{0},\dots,x_{n-1}$ is a sufficiently dense
partial orbit, we may assume that the mesh size, $\displaystyle\theta$, of the
spline approximation is $\displaystyle\mathcal{O}(n^{-1})$. So,

$\displaystyle\|\rho_{\nu_{\theta}}-\rho_{\nu_{0}}\|_{L^{1}}\lesssim\|\Pi_{\theta}-\Pi_{0}\|_{{\operatorname{BV}},L^{1}}\lesssim\|\mathcal{L}_{\theta}-\mathcal{L}_{0}\|_{{\operatorname{BV}},L^{1}}\lesssim\sqrt{\theta}=n^{-1/2}$

where $\displaystyle\nu_{\theta}$ is the unique acip of $\displaystyle
g_{\theta}$ and the implied constants are independent of
$\displaystyle\theta$. The estimates (from left to right) follow from the
proof of 3.2, the Cauchy integral representation of
$\displaystyle\Pi_{\theta}$ given in Section 3.2 and (5.5), respectively.

We recall from Section 3.2 that for all $\displaystyle k>0$,

$\frac{1}{k}{\mathbb{E}}_{\mu}(S_{\theta,k}(h))=A_{\theta}+\frac{1}{ik}m(\Pi^{(1)}_{\theta}\rho_{\mu})+\mathcal{O}(\kappa^{k})$

Since the convergence rate of $\displaystyle A_{\theta}\to A_{0}$ is
independent of the choice of the initial measure, we write,

$\displaystyle\displaystyle A_{\theta}-A_{0}$
$\displaystyle\displaystyle=\frac{1}{k}\left({\mathbb{E}}_{\nu_{\theta}}(S_{k,\theta}(h))-{\mathbb{E}}_{\nu_{0}}(S_{k,0}(h))\right)+o(k^{-1})$
$\displaystyle\displaystyle=\frac{1}{k}\sum_{j=1}^{k-1}\left[{\mathbb{E}}_{\nu_{\theta}}(h\circ
g_{\theta}^{j})-{\mathbb{E}}_{\nu_{0}}(h\circ g_{0}^{j})\right]+o(k^{-1})$
$\displaystyle\displaystyle=\frac{1}{k}\sum_{j=1}^{k-1}\left[{\mathbb{E}}_{\nu_{\theta}}(h)-{\mathbb{E}}_{\nu_{0}}(h)\right]+o(k^{-1})$
$\displaystyle\displaystyle|A_{\theta}-A_{0}|$
$\displaystyle\displaystyle\leq\frac{1}{k}\sum_{j=1}^{k-1}\left|\int
h(\rho_{\nu_{\theta}}-\rho_{\nu_{0}})\,dm\right|+o(k^{-1})$
$\displaystyle\displaystyle\leq\|h\|_{L^{\infty}}\|\rho_{\nu_{\theta}}-\rho_{\nu_{0}}\|_{L^{1}}+o(k^{-1})\lesssim
n^{-1/2}+o(k^{-1})$

Above, we use that $\displaystyle\nu_{\theta}$ is $\displaystyle
g_{\theta}-$invariant. Next, choosing $\displaystyle k=\mathcal{O}(n^{1/2})$,
we can obtain the best possible rate of convergence of the asymptotic means:

$\displaystyle|A^{*}-A|=|A_{\theta}-A_{0}|\lesssim n^{-1/2}.$

Now, we focus on the rate of convergence of the variance. We recall that for
all $\displaystyle k>0$,

$\frac{1}{k}{\mathbb{E}}_{\mu}([S_{\theta,k}(h_{\theta})-kA_{\theta}]^{2})=\sigma^{2}_{\theta}-\frac{1}{k}m(\Pi^{(2)}_{\theta}\rho_{\mu})+\mathcal{O}(\kappa^{n})$

Then,

$\displaystyle\displaystyle\sigma^{2}_{\theta}-\sigma^{2}_{0}$
$\displaystyle\displaystyle=\frac{1}{k}\left({\mathbb{E}}_{\nu_{\theta}}([S_{k,\theta}(h)-kA_{\theta}]^{2})-{\mathbb{E}}_{\nu_{0}}([S_{k,0}(h)-kA_{0}]^{2})\right)+o(k^{-1})$
$\displaystyle\displaystyle=\frac{1}{k}\left({\mathbb{E}}_{\nu_{\theta}}([S_{k,\theta}(h)]^{2})-{\mathbb{E}}_{\nu_{0}}([S_{k,0}(h)]^{2}\right)-k[A^{2}_{\theta}-A^{2}_{0}]+o(k^{-1})$

Since,
$\displaystyle|A^{2}_{\theta}-A^{2}_{0}|\lesssim|A_{\theta}-A_{0}|\lesssim
n^{-1/2}$ and the second term is $\displaystyle\mathcal{O}(kn^{-1/2})$.

To estimate the first term, note that,

$\displaystyle\displaystyle|{\mathbb{E}}_{\nu_{\theta}}([S_{k,\theta}(h)]^{2})-{\mathbb{E}}_{\nu_{0}}([S_{k,0}(h)]^{2})|$
$\displaystyle\displaystyle=\left|\sum_{j=0}^{k-1}\sum_{l=0}^{k-1}{\mathbb{E}}_{\nu_{\theta}}(h\circ
g_{\theta}^{j}\cdot h\circ g_{\theta}^{l})-{\mathbb{E}}_{\nu_{0}}(h\circ
g_{0}^{j}\cdot h\circ g_{0}^{l})\right|$
$\displaystyle\displaystyle=2\left|\sum_{j=1}^{k-1}(k-j)[{\mathbb{E}}_{\nu_{\theta}}(h\circ
g_{\theta}^{j}\cdot h)-{\mathbb{E}}_{\nu_{0}}(h\circ g_{0}^{j}\cdot
h)]\right|$ $\displaystyle\displaystyle\leq
2\sum_{j=1}^{k-1}(k-j)\left|{\mathbb{E}}_{m}(h\cdot\mathcal{L}^{j}_{\theta}(h\rho_{\nu_{\theta}}))-{\mathbb{E}}_{m}(h\cdot\mathcal{L}^{j}_{0}(h\rho_{\nu_{0}}))\right|$
$\displaystyle\displaystyle\leq
2\|h\|_{L^{\infty}}\sum_{j=1}^{k-1}(k-j)\|\mathcal{L}^{j}_{\theta}(h\rho_{\nu_{\theta}})-\mathcal{L}^{j}_{0}(h\rho_{\nu_{0}})\|_{L^{1}},$
$\displaystyle\displaystyle\|\mathcal{L}^{j}_{\theta}(h\rho_{\nu_{\theta}})-\mathcal{L}^{j}_{0}(h\rho_{\nu_{0}})\|_{L^{1}}$
$\displaystyle\displaystyle\leq\|\mathcal{L}^{j}_{\theta}(h(\rho_{\nu_{\theta}}-\rho_{\nu_{0}}))\|_{L^{1}}+\|(\mathcal{L}^{j}_{\theta}-\mathcal{L}^{j}_{0})(h\rho_{\nu_{0}})\|_{L^{1}}$
$\displaystyle\displaystyle\leq\sup_{\theta,r}\|\mathcal{L}^{r}_{\theta}\|_{L^{1},L^{1}}\cdot\|h\|_{L^{\infty}}\|\rho_{\nu_{\theta}}-\rho_{\nu_{0}}\|_{L^{1}}+\|\mathcal{L}^{j}_{\theta}-\mathcal{L}^{j}_{0}\|_{{\operatorname{BV}},L^{1}}\|h\rho_{\nu_{0}}\|_{{\operatorname{BV}}}$
$\displaystyle\displaystyle\lesssim
n^{-1/2}+\|\mathcal{L}^{j}_{\theta}-\mathcal{L}^{j}_{0}\|_{{\operatorname{BV}},L^{1}},$

and

$\displaystyle\displaystyle\|\mathcal{L}^{j}_{\theta}-\mathcal{L}^{j}_{0}\|_{{\operatorname{BV}},L^{1}}=\sum_{r=0}^{j-1}\|\mathcal{L}^{r}_{\theta}(\mathcal{L}_{\theta}-\mathcal{L}_{0})\mathcal{L}^{j-1-r}_{0}\|_{{\operatorname{BV}},L^{1}}$
$\displaystyle\displaystyle\lesssim\sup_{\theta,r}\|\mathcal{L}^{r}_{\theta}\|_{L^{1},L^{1}}\|\mathcal{L}_{\theta}-\mathcal{L}_{0}\|_{{\operatorname{BV}},L^{1}}\sum_{r=0}^{j-1}\|\mathcal{L}^{r}_{0}\|_{{\operatorname{BV}},{\operatorname{BV}}}$
$\displaystyle\displaystyle\lesssim jn^{-1/2},$

Therefore, we have

$\displaystyle|{\mathbb{E}}_{\nu_{\theta}}([S_{k,\theta}(h)]^{2})-{\mathbb{E}}_{\nu_{0}}([S_{k,0}(h)]^{2})|\lesssim
k^{3}n^{-1/2}.$

Finally,

$\displaystyle\displaystyle\sigma^{2}_{\theta}-\sigma^{2}_{0}\lesssim
k^{2}n^{-1/2}+kn^{-1/2}+o(k^{-1}),$

and choosing $\displaystyle k=\mathcal{O}(n^{1/6})$, we can obtain the best
possible rate of convergence of the standard deviations,

$\displaystyle|\sigma^{*}-\sigma|\lesssim|\sigma^{2}_{\theta}-\sigma^{2}_{0}|\lesssim
n^{-1/6}.$

So, in the case of piecewise expanding maps, we have the following improved
asymptotic accuracy of the non-pivoted bootstrap.

$\displaystyle\sup_{z\in\mathbb{R}}\Big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=\mathcal{O}_{\text{a.s.}}(n^{-1/6}).$

######  Example 4 (Smooth expanding maps, a continuation of Section 5.1.1):

As in the previous example, we assume that the mesh-size,
$\displaystyle\theta$, of the spline approximation,
$\displaystyle\mathfrak{s}$, obtained by $\displaystyle x_{0},\dots,x_{n-1}$
is $\displaystyle\mathcal{O}(n^{-1})$. Since

$\displaystyle\|g-\mathfrak{s}\|_{L^{\infty}}+\|g^{\prime}-\mathfrak{s}^{\prime}\|_{L^{\infty}}=\mathcal{O}(\theta^{2})+\mathcal{O}(\theta)=\mathcal{O}(n^{-1}),$

from Section 5.1, we have that

$\displaystyle\|\mathcal{L}_{\theta}-\mathcal{L}_{0}\|_{{\operatorname{BV}},L^{1}}\lesssim
n^{-1}.$

Arguing as in the previous example, for all $\displaystyle k>0$,

$\displaystyle|A^{*}-A|=|A_{\theta}-A_{0}|\lesssim
n^{-1}+\mathcal{O}(k^{-1}).$

It is easy to see that choosing $\displaystyle k=\mathcal{O}(n)$ gives the
best possible rate of convergence $\displaystyle\mathcal{O}(n^{-1})$. So,
arguing as in the previous example,.

$\displaystyle|\sigma^{*}-\sigma|\lesssim|\sigma^{2}_{\theta}-\sigma^{2}_{0}|\lesssim
k^{2}n^{-1}+kn^{-1}+o(k^{-1}),$

and choosing $\displaystyle k=\mathcal{O}(n^{1/3})$, the best possible rate of
convergence is $\displaystyle\mathcal{O}(n^{-1/3})$.

Therefore, for smooth expanding maps of $\displaystyle\mathbb{T}$, we have we
have the following improved asymptotic accuracy of the non-pivoted bootstrap.

$\displaystyle\sup_{z\in\mathbb{R}}\Big{|}{\mathbb{P}}_{\mu}(T_{n}\leq
z)-{\mathbb{P}}(T^{*}_{n}\leq
z\,|\,x_{0},\dots,x_{n-1})\Big{|}=\mathcal{O}_{\text{a.s.}}(n^{-1/3}).$

##  9\. Computer Simulation of the Bootstrap

###  9.1. Data generating processes

We let the choices of sample size be $\displaystyle n=25,50,100$ and take the
doubling map, the logistic map, and the drill map below as the choices for the
transformation function $\displaystyle g$.

####  9.1.1. Doubling map

The doubling map also called the dyadic transformation is given by

$\displaystyle g(x)=2x\mod 1,\,\,\,\,\,x\in[0,1].$

In general, $\displaystyle r$-adic transformations, including dyadic
transformations, allow us to study the statistical properties of the digits of
real numbers. Because of its simple yet chaotic nature, von Neumann proposed
the doubling map as a random number generator [76], and later, Rényi studied
its statistical properties [1]. In fact, it is an expanding (therefore,
hyperbolic) map, exhibits sensitive dependence on initial conditions, has a
unique exponentially mixing acip – the Lebesgue measure $\displaystyle[0,1]$ –
and also, is a Bernoulli map. Further, it is a $\displaystyle C^{\infty}-$map
of $\displaystyle\mathbb{T}$ with constant first derivative
$\displaystyle=2>1$, and hence, falls into the class of examples considered in
Section 5.1.

One could also consider the doubling map as the following map of the binary
expansions of data

(9.1) $\sum_{j=1}^{\infty}w_{j}2^{-j}\mapsto\sum_{j=1}^{\infty}w_{j+1}2^{-j}.$

In other words, the doubling map neglects the first binary digit but shift all
the other digit to the left by one unit. Unfortunately, in computer software,
the initial state $\displaystyle x_{0}$ only has a binary expansion up to a
finite order $\displaystyle K$, e.g., $\displaystyle K=32$ when $\displaystyle
x_{0}$ is a uniform random variable generated in R software; as a result of
this round-off error, by (9.1), after $\displaystyle K$ iterations the orbit
of the simulated doubling map will always end up at zero [7]. To overcome this
problem, we let $\displaystyle x_{i}=g(x_{i-1})+2^{-20}{\varepsilon}_{i}$,
where $\displaystyle\\{{\varepsilon}_{i}\\}$ are independent
Bernoulli$\displaystyle(1/2)$ random variables conditional on $\displaystyle
x_{i}\in[0,1]$. By the shadowing lemma, this perturbed orbit of
$\displaystyle\\{x_{i}\\}$ stays uniformly close to an unperturbed orbit; see
[60, p. 18].

####  9.1.2. Drill map

The motion of a rotary drill induces a piecewise expanding map of the
interval, $\displaystyle g:[0,1]\to[0,1]$, defined as follows:

$\displaystyle\displaystyle\alpha$
$\displaystyle\displaystyle=\frac{\Lambda}{(\Lambda-1)},\,\,\,q_{[\Lambda]-k}=\max\left\\{0,\frac{1}{2}\cdot\frac{\Lambda-1-k}{\Lambda-1}\right\\},\
\ k=1,\dots,[\Lambda],$ $\displaystyle\displaystyle d_{\Lambda}(x)$
$\displaystyle\displaystyle=\alpha\left(k-\sqrt{k^{2}-\frac{k}{\alpha}(k+1-2x)}\right),\
\ q_{[\Lambda]-k}<x\leq q_{[\Lambda]-k+1},\ k=1,2,\dots,[\Lambda],$
$\displaystyle\displaystyle g^{1/2}(x)$
$\displaystyle\displaystyle=x+d_{\Lambda}(x)\ \ \text{(mod $\displaystyle
1$)},\,\,\,\text{and}\,\,\,\,g(x)=g^{1/2}(g^{1/2}(x)),$

where $\displaystyle[\Lambda]$ denote the integer part of
$\displaystyle\Lambda$ and $\displaystyle\Lambda$ is a parameter indicating
the influence of gravity on fluid motion. We set $\displaystyle\Lambda=3$ in
our simulation and include the corresponding transformation function
$\displaystyle g$ in Figure 4.

Figure 4. Drill map when $\displaystyle\Lambda=3$.

This map was first considered in [43] to model the movement of an oil drill
with real world engineering applications in mind. This is a particular example
of transformations discussed in Section 5.2; see also [6, Section 1.2, Section
13.3].

####  9.1.3. Logistic map

The logistic map is useful as a discrete-time population model of various
biological species; see [50] for an in-depth discussion. The logistic map with
parameter $\displaystyle r\in(0,4]$ is defined by

$\displaystyle g(x)=rx(1-x),\,\,\,\,x\in[0,1].$

Here, we focus on the case $\displaystyle r=4$. It is well-known that
$\displaystyle g$ with $\displaystyle r=4$ exhibits sensitive dependency to
initial conditions and has a unique acip (despite not being hyperbolic); see
p. 34–35 of [60].

Even though $\displaystyle g$ does not belong to any of the examples in
Section 5, its simulations results are includes here because they were
comparable to those of the previous two examples. This is an indication that
asymptotic accuracy of the bootstrap may hold even in the case of mostly
hyperbolic maps considered in [79]. In fact, we believe that we can establish
continuous Edgeworth expansions for such maps based on ideas in [22]. This
will be pursued in the future.

Similar to the doubling map, the logistic map suffers from the round-off error
of computer software. As a remedy, we first transform the logistic map to the
tent map with $\displaystyle\mathcal{T}:[0,1]\to[0,1]$ defined by
$\displaystyle\mathcal{T}(x)=2\arcsin(\sqrt{x})/\pi$ and then add perturbation
to the tent map; see p. 33 of [60]. Specifically, we let

$\displaystyle\displaystyle y_{i-1}$
$\displaystyle\displaystyle=\mathcal{T}(x_{i-1}),$ $\displaystyle\displaystyle
y_{i}$
$\displaystyle\displaystyle=\begin{cases}2y_{i-1}+2^{-20}{\varepsilon}_{i},&\text{for}\
x<1/2,\\\ 2(1-y_{i-1})+2^{-20}(1-{\varepsilon}_{i}),&\text{for}\ x\geq
1/2,\end{cases}$ $\displaystyle\displaystyle x_{i}$
$\displaystyle\displaystyle=\mathcal{T}^{-1}(y_{i})$

where $\displaystyle\\{{\varepsilon}_{i}\\}$ are independent
Bernoulli$\displaystyle(1/2)$ random variables conditional on $\displaystyle
y_{i}\in[0,1]$. By the shadowing lemma, as in the doubling map context, the
perturbed orbit $\displaystyle\\{y_{i}\\}$ stays uniformly close to some
unperturbed orbit of the tent map, and in turn, $\displaystyle\\{x_{i}\\}$
shadows an orbit of the logistic map.

###  9.2. Quantity of interest

We aim to construct two-sided, upper-bounded, and lower-bounded 95% confidence
intervals for spatial average $\displaystyle A$. When constructing
$\displaystyle A$, we let $\displaystyle h(x)=x$, $\displaystyle h(x)=x^{2}$,
and $\displaystyle h(x)=x^{4}$; these choices of $\displaystyle h$ are closely
related to the mean, variance, and kurtosis, which have wide applications in a
variety of fields.

###  9.3. Bootstrap methods

To construct the confidence intervals mentioned in section 9.2, we apply the
pivoted bootstrap in Algorithm LABEL:alg:pivot and the non-pivoted bootstrap
in Algorithm LABEL:alg:nonpivot with bootstrap iteration $\displaystyle
B=1000$. We include more details of these algorithms below.

####  9.3.1. Estimation of the transformation $\displaystyle g$

Throughout our simulation, we assume that the support and the location of
discontinuities of transformation $\displaystyle g$ are known. Given this, we
approximate the transformation function $\displaystyle g$ by piecewise cubic
spline. When making extrapolations, we apply “FMM” [24] cubic splines in non-
pivoted bootstrap and “natural” cubic splines in pivoted bootstrap. When the
fitted value lies outside the given support, we take it to be the nearest
boundary point of the support, and then, move it inside the support by adding
or subtracting another small perturbation of $\displaystyle 2^{-20}$. In our
preliminary simulation, we used piecewise linear splines but found their
performance to be inferior to that of the piecewise cubic splines.

####  9.3.2.  Generation of initial bootstrap data $\displaystyle x_{0}^{*}$

In our simulation, we implement Line LABEL:lin:bootInit_pivoted of Algorithm
LABEL:alg:pivot and Line LABEL:lin:bootInit_non_pivoted of Algorithm
LABEL:alg:nonpivot, the generation of the initial bootstrap data
$\displaystyle x_{0}^{*}$, as follows:

1

Input: Data $\displaystyle x_{0},\dots,x_{n-1}$

Output: Initial bootstrap state $\displaystyle x_{0}^{*}$

2

3while  _$\displaystyle x_{0}^{*}\notin\text{support of}\ g$_ do

4 $\displaystyle ucv\leftarrow$ unbiased cross-validated bandwidth [70]

5 $\displaystyle\mathfrak{b}\leftarrow ucv/4$

6 Sample $\displaystyle{\varepsilon}_{0}\sim\text{Normal}(0,\mathfrak{b}^{2})$

7 Sample $\displaystyle x_{0}^{\prime}$ from
$\displaystyle\\{x_{0},\dots,x_{n-1}\\}$, independently of
$\displaystyle{\varepsilon}_{0}$

8 $\displaystyle x_{0}^{*}\leftarrow x_{0}^{\prime}+{\varepsilon}_{0}$

9

10 end while

return _$\displaystyle x_{0}^{*}$_

Algorithm 9.1 Generation of initial bootstrap data

######  Remark 9.1:

Generating $\displaystyle x_{0}^{*}$ with Algorithm 9.1 is equivalent to
generating $\displaystyle x_{0}^{*}$ from a kernel-estimated density function,
where the kernel is Gaussian and the bandwidth is given by
$\displaystyle\mathfrak{b}$; see [72, p. 471]. On the other hand, Algorithm
9.1 is computationally more economical than directly generating data from the
kernel-estimated density function; see the footnote on p. 1068 of [34].

######  Remark 9.2:

When the bandwidth is large, e.g., $\displaystyle\mathfrak{b}=ucv$, where
$\displaystyle ucv$ is the unbiased cross-validated bandwidth as in [70], the
bootstrap initial state could take a value that the true initial state rarely
takes, e.g., 0.99. When starting with such a extreme bootstrap initial state,
the bootstrap trajectory of $\displaystyle\\{x_{i}^{*}\\}$ could also be
extreme, and consequently the bootstrap distribution may contain too many
outliers. As a remedy, we let the bandwidth be rather small, i.e.,
$\displaystyle\mathfrak{b}=ucv/4$.

####  9.3.3. Estimation of long-run variance $\displaystyle\sigma^{2}$

Since we are not aware of any theoretically-justified estimator for
$\displaystyle\sigma$ in the dynamical system setting, we first obtain the
true $\displaystyle\sigma$ with Monte-Carlo simulation and then replace
$\displaystyle\widehat{\sigma}$ by this simulated true $\displaystyle\sigma$.
Similarly, we replace $\displaystyle\widehat{\sigma}^{*,b}$ by
$\displaystyle\tilde{\sigma}^{*}$, where

$\displaystyle\tilde{\sigma}^{*}\coloneqq\sqrt{\frac{1}{B}\sum_{b=1}^{B}\Big{(}\frac{S_{n}^{*,b}(h)-\bar{S}_{n}^{*}(h)}{\sqrt{n}}\Big{)}^{2}}.$

Indeed, in a similar way to Remark LABEL:rem:truebootexp, when both
$\displaystyle B$ and $\displaystyle n$ are large, recalling
$\displaystyle\sigma^{*}$ defined in (LABEL:eq:bootQuantities),

$\displaystyle\tilde{\sigma}^{*}\approx\sqrt{{\mathbb{E}}^{*}\bigg{(}\Big{(}\frac{S_{n}^{*}(h)-\bar{S}_{n}^{*}(h)}{\sqrt{n}}\Big{)}^{2}\bigg{)}}\approx\sqrt{{\mathbb{E}}^{*}\bigg{(}\Big{(}\frac{S_{n}^{*}(h)-nA^{*}}{\sqrt{n}}\Big{)}^{2}\bigg{)}}\approx\sigma^{*}.$

###  9.4. Non-bootstrap methods

####  9.4.1. Gaussian approximation

To construct the confidence intervals in section 9.2, alternatively we can
apply the Gaussian approximation, namely

$\displaystyle\frac{S_{n}(h)-nA}{\sqrt{n}\widehat{\sigma}}\Rightarrow N(0,1).$

As in Section 9.3.3, since we do not know any theoretically-justified
estimator for $\displaystyle\sigma$ in the dynamical system setting, we
replace $\displaystyle\widehat{\sigma}$ by the simulated true
$\displaystyle\sigma$ throughout our simulation.

####  9.4.2. $\displaystyle t$-approximation

To construct the confidence intervals, we can also apply a $\displaystyle
t$-approximation, which relies on the presumption that

$\displaystyle\frac{S_{n}(h)-nA}{\sqrt{n}s}\Rightarrow t_{n-1},$

where $\displaystyle s$ is sample standard deviation of
$\displaystyle\\{x_{i},i=0,\dots,n-1\\}$ and $\displaystyle t_{n-1}$ is a
$\displaystyle t$-distribution with degrees of freedom $\displaystyle n-1$.

###  9.5. Results

To obtain the empirical coverage of the confidence intervals, we run 700
iterations. The results obtained are included in Tables 1, 3, 2, 4, 6, 5, 7, 9
and 8. When reporting the data, we use the following abbreviations

  * –

t ($\displaystyle t$-approximation),

  * –

npboot (non-pivoted bootstrap),

  * –

Gaussian (Gaussian approxiamtion), and

  * –

pboot (pivoted bootstrap)

to differentiate how the $\displaystyle 95\%$ confidence intervals were
generated.

First, recall that the $\displaystyle t$-approximation and non-pivoted
bootstrap do not require any prior knowledge of the long-run variance
$\displaystyle\sigma^{2}$, while the Gaussian approximation and pivoted-
bootstrap heavily depend on $\displaystyle\sigma^{2}$. From Tables 1, 3, 2, 4,
6, 5, 7, 9 and 8, we see that, in general, the Gaussian approximation and
pivoted-bootstrap prevail over the $\displaystyle t$-approximation and non-
pivoted bootstrap, and hence, as expected, the knowledge of
$\displaystyle\sigma^{2}$ improves the accuracy.

Second, the non-pivoted bootstrap significantly outperforms the $\displaystyle
t$-approximation. In particular, the non-pivoted bootstrap does not suffer
from under-coverage in case of the doubling map and over-coverage in case of
the logistic map whereas the $\displaystyle t$-approximation suffers from
both. Indeed, the non-pivoted bootstrap has a performance that almost matches
the Gaussian approximation and the pivoted bootstrap. Consequently, when
$\displaystyle\sigma$ is unknown and cannot be easily estimated, the non-
pivoted bootstrap may be preferable.

Third, the Gaussian approximation and pivoted-bootstrap give comparable
results; none of them uniformly dominates the other. Since Gaussian
approximation fails to capture the asymmetry of the finite-sample
distribution, it performs slightly inferior to the pivoted-bootstrap on one-
sided confidence intervals. However, since the under-coverage on one tail of
the distribution is likely compensated by the over-coverage on the other tail,
the Gaussian approximation has a slight advantage when two-sided confidence
intervals are computed. In summary, when $\displaystyle\sigma$ is given, both
the Gaussian approximation and pivoted-bootstrap can be considered.

####  9.5.1. Two-sided confidence interval

| $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.753 | 0.740 | 0.749 | 0.729 | 0.709 | 0.760 | 0.717 | 0.729 | 0.747
npboot | 0.946 | 0.960 | 0.954 | 0.944 | 0.946 | 0.954 | 0.941 | 0.966 | 0.947
Gaussian | 0.957 | 0.957 | 0.947 | 0.956 | 0.957 | 0.944 | 0.960 | 0.957 | 0.953
pboot | 0.940 | 0.940 | 0.941 | 0.964 | 0.959 | 0.949 | 0.940 | 0.960 | 0.936
Table 1. Empirical coverage of two-sided $\displaystyle 95\%$ confidence intervals generated when data is generated by the doubling map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.943 | 0.934 | 0.923 | 0.913 | 0.926 | 0.924 | 0.856 | 0.907 | 0.910
npboot | 0.937 | 0.951 | 0.947 | 0.940 | 0.941 | 0.946 | 0.893 | 0.924 | 0.941
Gaussian | 0.956 | 0.954 | 0.953 | 0.947 | 0.959 | 0.960 | 0.951 | 0.963 | 0.953
pboot | 0.904 | 0.924 | 0.933 | 0.945 | 0.950 | 0.963 | 0.960 | 0.944 | 0.964
Table 2. Empirical coverage of two-sided $\displaystyle 95\%$ confidence intervals when data is generated by the drill map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.946 | 0.940 | 0.957 | 0.990 | 0.996 | 0.994 | 1.000 | 1.000 | 1.000
npboot | 0.923 | 0.930 | 0.949 | 0.934 | 0.950 | 0.949 | 0.939 | 0.953 | 0.929
Gaussian | 0.957 | 0.941 | 0.953 | 0.959 | 0.957 | 0.947 | 0.957 | 0.950 | 0.949
pboot | 0.960 | 0.939 | 0.957 | 0.949 | 0.940 | 0.943 | 0.963 | 0.946 | 0.964
Table 3. Empirical coverage of two-sided $\displaystyle 95\%$ confidence
intervals when data is generated by the logistic map.

####  9.5.2. Upper-bounded confidence interval

| $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.850 | 0.814 | 0.834 | 0.784 | 0.781 | 0.829 | 0.741 | 0.791 | 0.793
npboot | 0.936 | 0.947 | 0.957 | 0.944 | 0.947 | 0.961 | 0.930 | 0.966 | 0.944
Gaussian | 0.949 | 0.963 | 0.953 | 0.970 | 0.969 | 0.949 | 0.980 | 0.976 | 0.966
pboot | 0.944 | 0.941 | 0.937 | 0.959 | 0.956 | 0.944 | 0.946 | 0.960 | 0.941
Table 4. Empirical coverage of upper-bounded $\displaystyle 95\%$ confidence intervals when data is generated by the doubling map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.939 | 0.930 | 0.924 | 0.903 | 0.903 | 0.910 | 0.820 | 0.871 | 0.871
npboot | 0.923 | 0.944 | 0.947 | 0.930 | 0.939 | 0.956 | 0.884 | 0.930 | 0.934
Gaussian | 0.963 | 0.961 | 0.963 | 0.956 | 0.963 | 0.961 | 0.987 | 0.974 | 0.974
pboot | 0.931 | 0.937 | 0.930 | 0.958 | 0.936 | 0.940 | 0.956 | 0.954 | 0.966
Table 5. Empirical coverage of upper-bounded $\displaystyle 95\%$ confidence intervals when data is generated by the drill map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.967 | 0.953 | 0.954 | 0.989 | 0.984 | 0.989 | 0.997 | 0.999 | 0.997
npboot | 0.939 | 0.926 | 0.966 | 0.941 | 0.963 | 0.950 | 0.953 | 0.954 | 0.941
Gaussian | 0.946 | 0.933 | 0.944 | 0.946 | 0.934 | 0.950 | 0.929 | 0.924 | 0.936
pboot | 0.946 | 0.944 | 0.943 | 0.937 | 0.946 | 0.950 | 0.936 | 0.937 | 0.960
Table 6. Empirical coverage of upper-bounded $\displaystyle 95\%$ confidence
intervals when data is generated by the logistic map.

####  9.5.3. Lower-bounded confidence interval

| $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.813 | 0.826 | 0.830 | 0.861 | 0.843 | 0.840 | 0.899 | 0.857 | 0.874
npboot | 0.947 | 0.954 | 0.959 | 0.934 | 0.947 | 0.947 | 0.950 | 0.961 | 0.963
Gaussian | 0.963 | 0.950 | 0.954 | 0.939 | 0.949 | 0.943 | 0.936 | 0.933 | 0.944
pboot | 0.951 | 0.941 | 0.954 | 0.957 | 0.959 | 0.957 | 0.951 | 0.954 | 0.946
Table 7. Empirical coverage of lower-bounded $\displaystyle 95\%$ confidence intervals when data is generated by the doubling map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.959 | 0.946 | 0.937 | 0.966 | 0.964 | 0.941 | 0.986 | 0.977 | 0.971
npboot | 0.984 | 0.960 | 0.947 | 0.960 | 0.944 | 0.946 | 0.957 | 0.944 | 0.947
Gaussian | 0.951 | 0.957 | 0.939 | 0.936 | 0.947 | 0.946 | 0.931 | 0.943 | 0.924
pboot | 0.914 | 0.940 | 0.951 | 0.921 | 0.940 | 0.950 | 0.964 | 0.930 | 0.957
Table 8. Empirical coverage of lower-bounded $\displaystyle 95\%$ confidence intervals when the data is generated by the drill map. | $\displaystyle h(x)=x$ | $\displaystyle h(x)=x^{2}$ | $\displaystyle h(x)=x^{4}$
---|---|---|---
| $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$ | $\displaystyle n=25$ | $\displaystyle n=50$ | $\displaystyle n=100$
t | 0.941 | 0.939 | 0.949 | 0.987 | 0.996 | 0.990 | 1.000 | 1.000 | 1.000
npboot | 0.934 | 0.951 | 0.939 | 0.954 | 0.950 | 0.953 | 0.933 | 0.950 | 0.941
Gaussian | 0.961 | 0.949 | 0.956 | 0.971 | 0.970 | 0.960 | 0.973 | 0.973 | 0.969
pboot | 0.946 | 0.944 | 0.956 | 0.959 | 0.950 | 0.947 | 0.971 | 0.950 | 0.956
Table 9. Empirical coverage of lower-bounded $\displaystyle 95\%$ confidence
intervals when data is generated by the logistic map.

## References

  * [1] Alfréd Rényi  “Representations for real numbers and their ergodic properties”  In _Acta Math Acad Sci Hungary_ 8, 1957, pp. 477–493
  * [2] J H Ahlberg, E N Nilson and J L Walsh  “The theory of splines and their applications”  Academic Press, 1967
  * [3] Eli Barkai  “Aging in subdiffusion generated by a deterministic dynamical system”  In _Physical review letters_ 90.10, 2003, pp. 104101
  * [4] L Mark Berliner  “Statistics, probability and chaos”  In _Statistical Science_ , 1992, pp. 69–90
  * [5] Patrice Bertail and Stéphan Clémençon  “Regenerative block bootstrap for Markov chains”  In _Bernoulli_ 12.4, 2006, pp. 689–712
  * [6] Abraham Boyarsky and Pawel Góra  “Laws of Chaos”  Birkhäuser Basel, 1997
  * [7] Abraham Boyarsky and Pawel Góra  “Why computers like Lebesgue measure”  In _Computers & Mathematics with Applications_ 16.4, 1988, pp. 321–329
  * [8] Henk Broer and Floris Takens  “Dynamical Systems and Chaos”  Springer-Verlag New York, 2011
  * [9] Steven L Brunton and J Nathan Kutz  “Data-driven science and engineering: Machine learning, dynamical systems, and control”  Cambridge University Press, 2019
  * [10] Oliver Butterley and Peyman Eslami  “Exponential mixing for skew products with discontinuities”  In _Trans. Amer. Math. Soc._ 369, 2010, pp. 783–803
  * [11] Sangit Chatterjee and Mustafa R Yilmaz  “Chaos, fractals and statistics”  In _Statistical Science_ 7.1, 1992, pp. 49–68
  * [12] George Contopoulos  “Order and chaos in dynamical astronomy”  Springer Science & Business Media, 2004
  * [13] Somnath Datta and William P McCormick  “Some continuous Edgeworth expansions for Markov chains with applications to bootstrap”  In _Journal of Multivariate Analysis_ 52.1, 1995, pp. 83–106
  * [14] J. De Simoi and C. Liverani  “Limit theorems for fast–slow partially hyperbolic systems”  In _Invent. math._ , 2018, pp. 811–1016
  * [15] R.. Devaney  “An Introduction to Chaotic Dynamical Systems”  Addison-Wesley, 1987
  * [16] Thomas Duriez, Steven L Brunton and Bernd R Noack  “Machine learning control-taming nonlinear dynamics and turbulence”  Springer, 2017
  * [17] B. Efron  “Bootstrap methods: another look at the jackknife”  In _Ann. Statist._ 7.1, 1979, pp. 1–26
  * [18] Carl-Gustav Esséen  “Fourier analysis of distribution functions. A mathematical study of the Laplace-Gaussian law”  In _Acta Math._ 77, 1945, pp. 1–125
  * [19] William Feller  “An Introduction to Probability Theory and Its Applications, Vol. 2”  Wiley, 1982
  * [20] K. Fernando and C. Liverani  “Edgeworth expansions for weakly dependent random variables”  In _Ann. Inst. H. Poincaré Probab. Statist_ , 2021
  * [21] Kasun Fernando and Pratima Hebbar  “Higher Order Asymptotics for Large Deviations – Part I”  In _Annales de l’I.H.P. Probabilités et statistiques_ 121.3–4, 2021, pp. 219–257
  * [22] Kasun Fernando and Françoise Pène  “Expansions in the local and the central limit theorems for dynamical systems”  In _preprint_ , 2020
  * [23] Déborah Ferré, Loïc Hervé and James Ledoux  “Regular perturbatons of V–geometrically ergodic Markov chains”  In _J. Appl. Prob_ 50.1, 2013, pp. 184–194
  * [24] George Elmer Forsythe  “Computer methods for mathematical computations”  In _Prentice-Hall series in automatic computation_ 259, 1977
  * [25] Jürgen Franke, Jens-Peter Kreiss and Enno Mammen  “Bootstrap of kernel smoothing in nonlinear time series”  In _Bernoulli_ 8.1, 2002, pp. 1–37
  * [26] Fabian Fröhlich, Fabian J Theis and Jan Hasenauer  “Uncertainty analysis for non-identifiable dynamical systems: Profile likelihoods, bootstrapping and more”  In _International Conference on Computational Methods in Systems Biology_ , 2014, pp. 61–72  Springer
  * [27] Y. Guivarc’h and J. Hardy  “Théorèmes limites pour une classe de chaînes de Markov et applications aux difféomorphismes d’Anosov”  In _Annales de l’I.H.P. Probabilités et statistiques_ 24.1, 1988, pp. 73–98
  * [28] Peter Hall  “The Bootstrap and Edgeworth Expansion”  Springer, 1992
  * [29] Hanyuan Hang and Ingo Steinwart  “A Bernstein-type inequality for some mixing processes and dynamical systems with an application to learning”  In _Annals of Statistics_ 45.2, 2017, pp. 708–743
  * [30] Hanyuan Hang, Ingo Steinwart, Yunlong Feng and Johan AK Suykens  “Kernel density estimation for dynamical systems”  In _The Journal of Machine Learning Research_ 19.1, 2018, pp. 1260–1308
  * [31] Daisuke Haraki, Tomoya Suzuki, Hiroki Hashiguchi and Tohru Ikeguchi  “Bootstrap nonlinear prediction”  In _Physical Review E_ 75.5, 2007, pp. 056212
  * [32] Hubert Hennion and Loïc Hervé  “Limit Theorems for Markov Chains and Stochastic Properties of Dynamical Systems by Quasi-Compactness”  Springer-Verlag, 2001
  * [33] Loïc Hervé and Françoise Pène  “The Nagaev-Guivarc’h method via the Keller-Liverani theorem”  In _Bull. Soc. Math. France_ 138.3, 2010, pp. 415–489
  * [34] Joel L Horowitz  “Bootstrap methods for Markov processes”  In _Econometrica_ 71.4, 2003, pp. 1049–1082
  * [35] Valerie Isham  “Statistical aspects of chaos: a review”  In _Networks and chaos—statistical and probabilistic aspects_ 50, Monogr. Statist. Appl. Probab.  Chapman & Hall, London, 1993, pp. 124–200
  * [36] Jens Ledet Jensen  “Chaotic dynamical systems with a view towards statistics: a review”  In _Networks and chaos-statistical and probabilistic aspects_ 50, Monogr. Statist. Appl. Probab.  Chapman & Hall, London, 1993, pp. 201–250
  * [37] Moritz Jirak, Wei Biao Wu and Ou Zhao  “Sharp connections between Berry-Esseen characteristics and Edgeworth expansions for stationary processes”  In _Trans. Am. Math. Soc._ 374.6, 2021, pp. 4129–4183
  * [38] Tosio Kato  “Perturbation Theory for Linear Operators”  Springer-Verlag, 1995
  * [39] Gerhard Keller  “Stochastic stability in some chaotic dynamical systems”  In _Monatshefte für Mathematik_ 94.4, 1982, pp. 313–333
  * [40] Gerhard Keller and Carlangelo Liverani  “Stability of the spectrum for transfer operators”  In _Ann. Scuola Norm. Sup. Pisa_ 28.4, 1999, pp. 141–152
  * [41] Joseph B Keller  “The probability of heads”  In _The American Mathematical Monthly_ 93.3, 1986, pp. 191–197
  * [42] Hans R Kunsch  “The jackknife and the bootstrap for general stationary observations”  In _Annals of Statistics_ , 1989, pp. 1217–1241
  * [43] A Lasota and P Rusek  “An application of ergodic theory to the determination of the efficiency of cogged drilling bits”  In _Archiwum Górnictwa_ 3, 1974, pp. 281–295
  * [44] A. Lasota and M.C. Mackey  “Chaos, Fractals, and Noise, Stochastic Aspects of Dynamics.”  Springer, 1994
  * [45] Regina Y Liu and Kesar Singh  “Moving blocks jackknife and bootstrap capture weak dependence”  In _Exploring the limits of bootstrap_ 225, 1992
  * [46] Carlangelo Liverani  “Decay of correlations for piecewise expanding maps”  In _Journal of Statistical Physics_ 78.3–4, 1995, pp. 1111–11297
  * [47] Carlangelo Liverani  “Invariant measures and their properties. A functional analytic point of view”  In _Dynamical Systems. Part II: Topological Geometrical and Ergodic Properties of Dynamics_ Pisa: Centro di Ricerca Matematica Ennio De Giorgi, Scuola Normale Superiore, 2003, pp. 185–237
  * [48] Huma Lodhi and David Gilbert  “Bootstrapping parameter estimation in dynamic systems”  In _International Conference on Discovery Science_ , 2011, pp. 194–208  Springer
  * [49] AM Lyapunov  “The general problem of the stability of motion”  In _Comm. Soc. Math. Kharkow_ , 1892
  * [50] Robert M. May  “Stability and Complexity in Model Ecosystems”, Princeton Landmarks in Biology  Princeton University Press, 1973
  * [51] Kevin McGoff, Xin Guo and Anastasia Deckard  “The Local Edge Machine: inference of dynamic models of gene regulation”  In _Genome biology_ 17.1, 2016, pp. 1–13
  * [52] Kevin McGoff, Sayan Mukherjee, Andrew Nobel and Natesh Pillai  “Consistency of maximum likelihood estimation for some dynamical systems”  In _Annals of Statistics_ 43.1, 2015, pp. 1–29
  * [53] Kevin McGoff, Sayan Mukherjee and Natesh Pillai  “Statistical inference for dynamical systems: A review”  In _Statistics Surveys_ 9, 2015, pp. 209–252
  * [54] Kevin McGoff and Andrew B Nobel  “Empirical risk minimization and complexity of dynamical models”  In _Annals of Statistics_ 48.4, 2020, pp. 2031–2054
  * [55] Kevin McGoff and Andrew B Nobel  “Empirical risk minimization for dynamical systems and stationary processes”  In _Information and Inference: A Journal of the IMA_ , 2021
  * [56] Valérie Monbet and Pierre-François Marteau  “Non parametric resampling for stationary Markov processes: The local grid bootstrap approach”  In _Journal of statistical planning and inference_ 136.10, 2006, pp. 3319–3338
  * [57] Sergey V. Nagaev  “Some Limit Theorems for Stationary Markov Chains”  In _Theory of Probability and Applications_ 2.4, 1959, pp. 378–406
  * [58] Raymundo Navarrete and Divakar Viswanath  “Prediction of dynamical time series using kernel based regression and smooth splines”  In _Electronic Journal of Statistics_ 12.2, 2018, pp. 2217–2237
  * [59] Andrew Nobel  “Consistent estimation of a dynamical map”  In _Nonlinear dynamics and statistics_ Springer, 2001, pp. 267–280
  * [60] Edward Ott  “Chaos in Dynamical Systems”  Cambridge University Press, 2002
  * [61] Efstathios Paparoditis and Dimitris N Politis  “A Markovian local resampling scheme for nonparametric estimators in time series analysis”  In _Econometric Theory_ , 2001, pp. 540–566
  * [62] Efstathios Paparoditis and Dimitris N Politis  “The local bootstrap for Markov processes”  In _Journal of Statistical Planning and Inference_ 108.1-2, 2002, pp. 301–328
  * [63] Françoise Pène  “Probabilistic limit theorems via the operator perturbation method under optimal moment assumptions” https://hal-cnrs.archives-ouvertes.fr/hal-03241597, 2021
  * [64] Henri Poincaré  “Sur les propriétés des fonctions définies par les équations aux différences partielles”  Gauthier-Villars, 1879
  * [65] Robert F Port and Timothy Van Gelder  “Mind as motion: Explorations in the dynamics of cognition”  MIT press, 1995
  * [66] Charles Pugh and Michael Shub  “Stable Ergodicity”  In _Bull. Am. Math. Soc._ 41.1, 2003, pp. 1–41
  * [67] MB Rajarshi  “Bootstrap in Markov-sequences based on estimates of transition density”  In _Annals of the Institute of Statistical Mathematics_ 42.2, 1990, pp. 253–268
  * [68] Murray Rosenblatt  “A central limit theorem and a strong mixing condition”  In _Proceedings of the National Academy of Sciences of the United States of America_ 42.1 National Academy of Sciences, 1956, pp. 43
  * [69] Walter Rudin  “Principles of mathematical analysis”, International series in pure and applied mathematics  McGraw-Hill Inc., 1976
  * [70] David W Scott and George R Terrell  “Biased and unbiased cross-validation in density estimation”  In _Journal of the American Statistical Association_ 82.400, 1987, pp. 1131–1146
  * [71] Krishna V Shenoy, Maneesh Sahani and Mark M Churchland  “Cortical control of arm movements: a dynamical systems perspective”  In _Annual review of neuroscience_ 36, 2013, pp. 337–359
  * [72] BW Silverman and GA Young  “The bootstrap: to smooth or not to smooth?”  In _Biometrika_ 74.3, 1987, pp. 469–479
  * [73] Ingo Steinwart and Marian Anghel  “Consistency of support vector machines for forecasting the evolution of an unknown ergodic dynamical system from observations with unknown noise”  In _Annals of Statistics_ 37.2, 2009, pp. 841–875
  * [74] Peter Turchin  “Complex population dynamics”  Princeton university press, 2013
  * [75] Hal R Varian  “Dynamical systems with applications to economics”  In _Handbook of mathematical economics_ 1, 1981, pp. 93–110
  * [76] John Von Neumann  “Various techniques used in connection with random digits”  In _John von Neumann, Collected Works_ 5, 1963, pp. 768–770
  * [77] Peter Walters  “An introduction to ergodic theory”  Springer Science & Business Media, 2000
  * [78] E Weinan  “A proposal on machine learning via dynamical systems”  In _Communications in Mathematics and Statistics_ 5.1, 2017, pp. 1–11
  * [79] Lai-Sang Young  “Statistical properties of dynamical systems with some hyperbolicity”  In _Ann. Math._ 147.3, 1998, pp. 585–650
  * [80] Bin Yu  “Density estimation in the $\displaystyle L^{\infty}$ norm for dependent data with applications to the Gibbs sampler”  In _Annals of Statistics_ 21.2, 1993, pp. 711–735
  * [81] Anna Zincenko, Sergei Petrovskii, Vitaly Volpert and Malay Banerjee  “Turing instability in an economic–demographic dynamical system may lead to pattern formation on a geographical scale”  In _Journal of the Royal Society Interface_ 18.177, 2021