arXiv:1209.0703v3 [math.ST] 4 May 2016

Bernoulli 22(3), 2016, 1808–1838
DOI: 10.3150/15-BEJ712

Markov Chain Monte Carlo
confidence intervals
YVES F. ATCHADÉ
1

University of Michigan, 1085 South University, Ann Arbor, 48109, MI, United States.
E-mail: yvesa@umich.edu

For a reversible and ergodic Markov chain {Xn , n ≥ 0} with invariant distribution π, we show
that a valid confidence interval for π(h) can be constructed whenever the asymptotic variance
σP2 (h) is finite and positive. We do not impose any additional condition on the convergence rate
of the Markov chain. The confidence interval is derived using the so-called fixed-b lag-window
estimator of σP2 (h). We also derive a result that suggests that the proposed confidence interval
procedure converges faster than classical confidence interval procedures based on the Gaussian
distribution and standard central limit theorems for Markov chains.
Keywords: Berry–Esseen bounds; confidence interval; lag-window estimators; martingale
approximation; MCMC; reversible Markov chains

1. Introduction
Confidence intervals play an important role in Monte Carlo simulation (Robert and
Casella [26], Asmussen and Glynn [1]). In Markov Chain Monte Carlo (MCMC), the
existing literature requires the Markov chain to be geometrically ergodic for the validity
of confidence interval procedures (Jones et al. [15], Flegal and Jones [8], Atchadé [3]). The
main objective of this work is to simplify some of these assumptions. We show that for a
reversible ergodic Markov chain, a valid confidence interval can be constructed whenever
the asymptotic variance itself is finite. No additional convergence rate assumption on the
Markov chain is required.
Let {Xn , n ≥ 0} be a reversible stationary Markov chain with invariant distribution
π. For h ∈ L2 (π), the asymptotic variance of h is denoted σP2 (h) (see (2) below for the
definition). A remarkable result by C. Kipnis and S. R. Varadhan (Kipnis and Varadhan
Pn
1 √
[19]) says that if 0 < σP2 (h) < ∞, then σP (h)
i=1 (h(Xi ) − π(h)) converges weakly to
n
R
def
N(0, 1) where π(h) = h(z)π(dz). In order to turn this result into a confidence interval
for π(h), an estimator σn of σP (h) is needed. A common practice consists in choosing
σn as a consistent estimator of σP (h). However, consistent estimation of σP (h) typically

This is an electronic reprint of the original article published by the ISI/BS in Bernoulli,
2016, Vol. 22, No. 3, 1808–1838. This reprint differs from the original in pagination and
typographic detail.
1350-7265

c

2016 ISI/BS

2

Y.F. Atchadé

requires further assumptions on the convergence rate of the Markov chain (typically geometric ergodicity), and on the function h. Instead of insisting on consistency, we consider
the so-called fixed-b approach developed by Kiefer, Vogelsang and Bunzel [18], Kiefer and
Vogelsang [17], where the proposed estimator σn is known to be inconsistent. Using this
inconsistent estimator we show in Theorem 2.2 that a Studentized analog of the Kipnis–
P
def
Varadhan’s theorem holds: if 0 < σP2 (h) < ∞, then Tn = σn1√n ni=1 (h(Xi ) − π(h)) converges weakly to a (non-Gaussian) distribution. The theorem extends to nonstationary
Markov chains that satisfy a very mild ergodicity assumption. To a certain extent, the
result is a generalization of Atchadé and Cattaneo [4] which establishes the same limit
theorem for geometrically ergodic (but not necessarily reversible) Markov chains. The result is particularly relevant for Markov chains with sub-geometric convergence rates. For
such Markov chains, the author is not aware of any result that guarantees the asymptotic
validity of confidence intervals. However, it is important to point out that the finiteness
of σP2 (h) carries some implications in terms of convergence rate of P , and is not always
easy to check. But the main point of this work is that the finiteness of σP2 (h) is all that
is needed for consistent confidence interval.
that there exists a pair of random
As we shall see, Theorem 2.2 comes from the fact P
n
variables (N, D), say, such that the joint process ( √1n i=1 (h(Xi ) − π(h)), σn2 ) converges
2
weakly to (σP (h)N, σP (h)D). As a result, σP (h) cancels out in the limiting distribution
of Tn . This approach to confidence intervals is closely related to the standardized time
series method of Schruben [29] (see also Glynn and Iglehart [9]), well known in operations
research. Indeed in its simplest form, the standardized time series method is the analog
of the fixed-b procedure using the batch-mean estimator with a fixed number of batches.
Despite this close connection, this paper focuses only on the fixed-b confidence interval.
We also compare the fixed-b lag-window estimators with the more commonly used
lag-window estimators. We limit this comparison to the case of geometrically ergodic
Markov chains. We prove in Theorem√ 2.6 that the convergence rate of the fixed-b lagwindow estimator is of order log(n)/ n, better than the fastest rate achievable by the
more commonly used lag-window estimator. Similar comparisons based on the convergence of Tn has been reported elsewhere in the literature. Jansson [13] studied stationary
Gaussian moving average models and established that the rate of convergence of Tn is
n−1 log(n). Sun, Phillips and Jin [30] obtained the rate n−1 , under the main assumption that the underlying process is Gaussian and stationary. It seems unlikely that the
convergence rate n−1 will hold without
√ the Gaussian assumption. However, it is unclear
whether the convergence rate log(n)/ n obtained in Theorem 2.6 is tight.
We organize the paper as follows. Section 2 contains the main results, including the
rate of convergence of the fixed-b lag-window estimator in Section 2.4. We present a
simulation example to illustrate the finite sample properties of the confidence intervals
in Section 2.5. All the main proofs are postponed to Section 3 and the Appendix.

1.1. Notation
Throughout the paper (X, B) denotes a measure space with a countably generated sigmaalgebra B with a probability measure of interest π. We denote L2 (π) the usual space of

Markov Chain Monte Carlo confidence intervals

3

L2 -integrable functions with respect to π, with norm k · k and associated inner product
h·i, and we denote L20 (π) the subspace of L2 (π) of functions orthogonal to the constants:
R
def
L20 (π) = {f ∈ L2 (π): f (x)π(dx) = 0}.
For a measurable function f : X → R, a probability measure ν on (X, B) and a
def R
def
Markov kernel Q on X, we use the notation: ν(f ) = f (x)ν(dx), f¯ = f − π(f ),
R
def
def
Qf (x) = f (y)Q(x, dy), and Qj f (x) = Q{Qj−1 f }(x), with Q0 f (x) = f (x). For V : X →
[0, ∞), we define LV as the space of all measurable real-valued functions f : X → R
def

s.t. |f |V = supx∈X |f (x)|/V (x) < ∞. For two probability measures ν1 , ν2 , we denotes
def
kν1 − ν2 ktv = sup|f |≤1 |ν1 (f ) − ν2 (f )|, the total variation distance between ν1 and ν2 ,
def

and kν1 − ν2 kV = sup{f,|f |V ≤1} |ν1 (f ) − ν2 (f )|, its V -norm generalization.
For sequences {an , bn } of real nonnegative numbers, the notation an . bn means that
an ≤ cbn for all n, and for some constant c that does not depend on n. For a random
sequence {Xn }, we write Xn = Op (an ) if the sequence |Xn |/an is bounded in probability.
We say that Xn = op (an ) if Xn /an converges in probability to zero as n → ∞.

2. Monte Carlo confidence intervals for reversible
Markov chains
Throughout the paper, P denotes a Markov kernel on (X, B) that is reversible with respect
to π. This means that for any pair f, g ∈ L2 (π), hf, P gi = hg, P f i. We assume that P
satisfies the following.
A1 For π-almost all x ∈ X,
lim kP n (x, ·) − πktv = 0.

n→∞

(1)

Remark 1. Assumption A1 is very basic. For instance, if P is φ-irreducible, and aperiodic (in addition to being reversible with respect to π), then A1 holds. If in addition P is
Harris recurrent, then (1) holds for all x ∈ X. If P is a Metropolis–Hastings kernel, Harris
recurrence typically follows from π-irreducibility. All these statements can be found, for
instance, in Tierney [31].
Throughout the section, unless stated otherwise, {Xn , n ≥ 0} is a (nonstationary)
Markov chain on (X, B) with transition kernel P and started at some arbitrary (but
fixed) point x ∈ X for which (1) holds. The Markov kernel P induces in the usual way
a self-adjoint operator (also denoted P ) on the Hilbert space L20 (π) that maps h 7→ P h.
This operator P admits a spectral measure E on [−1, 1], and for h ∈ L20 (π) we will
def
write µh (·) = hh, E(·)hi for the associated nonnegative Borel measure on [−1, 1]. Assumption A1 implies that µh does not charge 1 or −1, that is µh ({−1, 1}) = 0. This is
Lemma 5 of Tierney [5].

4

Y.F. Atchadé

2.1. Confidence interval for π(h)
Let h ∈ L20 (π). We define
def

σP2 (h) =

Z

1

−1

1+λ
µh (dλ),
1−λ

(2)

that we call the asymptotic variance of h. The terminology comes from the fact that if the
Markov chain is assumed stationary, a calculation (see, e.g., Häggström and Rosenthal
[11], Theorem 4) using the properties of the spectral measure µh gives
"

lim nE

n→∞

n

−1

n
X

!2 #

h(Xk )

k=1

= σP2 (h).

(3)

For nonstationary Markov chains, such as the one considered in this paper, it is unclear
whether (3) continues to hold in complete generality. The estimation of σP2 (h) is often
2
of interest because when (3) holds,
mean squared error of
Pn σP (h)/n approximates the
2
−1
h(X
).
An
estimate
of
σ
(h)
is often also sought
the Monte Carlo estimate n
k
P
k=1
in order to exploit the Kipnis–Varadhan theorem for confidence interval purposes. It is
known (Häggström and Rosenthal [11], Theorem 4) that σP2 (h) can also be written as
+∞
X

σP2 (h) =

ℓ=−∞

γ|ℓ| (h),

(4)

def

where for ℓ ≥ 0, γℓ (h) = hh, P ℓ hi. This suggests the so-called lag-window estimator of
σP2 (h)
def
σb2n =

n−1
X

ℓ=−n+1


ℓ
γn,|ℓ| ,
w
bn

def



where γn,ℓ = n−1

n−ℓ
X
j=1

(5)

(h(Xj ) − π̂n (h))(h(Xj+ℓ ) − π̂n (h)).

Pn
In the above display, π̂n (h) = n−1 k=1 h(Xk ), 1 ≤ bn ≤ n is an integer such that bn → ∞,
as n → ∞, and w : R → R is an even function (w(−x) = w(x)) with support [−1, 1], that
is, w(x) 6= 0 on (−1, 1) and w(x) = 0 for |x| ≥ 1. Since w has support [−1, 1], the actual
range for ℓ in the summation defining σb2n is −bn + 1 ≤ ℓ ≤ bn − 1.
The lag-window estimator σb2n can be applied more broadly in time series and the
method has a long history. Some of the earlier work go back to the 1950s (Grenander and
Rosenblatt [10], Parzen [24]). Convergence results specific to nonstationary Markov chains
have been established recently (see, e.g., Damerdji [6], Flegal and Jones [8], Atchadé [3]
and the references therein); however, under assumptions that are much stronger than

Markov Chain Monte Carlo confidence intervals

5

A1. It remains an open problem whether σb2n can be shown to converge to σP2 (h) assuming only A1. In particular, the author is not aware of any result that establishes the
consistency of σb2n without assuming that P is geometrically ergodic.
However, if the goal is to construct a confidence interval for π(h), we will now see that
it is enough to assume A1 and σP2 (h) < ∞. Consider the lag-window estimator obtained
by setting bn = n. This writes
def

σn2 =

n−1
X

w

ℓ=−n+1

 
ℓ
γn,|ℓ| .
n

(6)

This estimator is well known to be inconsistent for estimating σP2 (h), but has recently
attracted a lot of interest in the Econometrics literature under the name of fixed-b asymptotics (Kiefer, Vogelsang and Bunzel [18], Kiefer and Vogelsang [17], Sun, Phillips and
Jin [30], see also Neave [22] for some pioneer work). This paper takes inspiration from
this literature. However, unlike these works, we exploit the Markov structure and we do
def R 1
not impose any stationary assumption. We introduce the function v(t) = 0 w(t − u) du,
t ∈ [0, 1], and the kernel φ : [0, 1] × [0, 1] → R, where
Z 1
φ(s, t) = w(s − t) − v(s) − v(t) +
v(t) dt,
s, t ∈ [0, 1].
(7)
0

We say that a kernel k : [0, 1] × [0,P1] →P
R is positive definite if for all n ≥ 1, all
n
n
a1 , . . . , an ∈ R, and t1 , . . . , tn ∈ [0, 1], i=1 j=1 ai aj k(ti , tj ) ≥ 0. We will assume that
the weight function w in (6) is such that the following holds.

A2 The function w : R → R is an even function, with support [−1, 1], and of class C 2 on
(−1, 1). Furthermore, the kernel φ defined in (7) is positive definite, and not identically
zero.

Example 1. Assumption A2 holds for the function w given by w(u) = (1 − u2 )1(−1,1) (u).
Indeed in this case, a simple calculation gives that φ(s, t) = 2(s − 0.5)(t − 0.5), which
(by its multiplicative form) is clearly positive definite. In this particular case, solving
R1
R1
φ(s, t)u(t) dt = αu(s) yields the unique eigenvalue α = 2 0 (t − 0.5)2 dt = 1/6.
0
A general approach to guarantee that φ as in (7) is positive definite is to start with a
positive definite function w, as the next lemma shows.

Lemma 2.1. Suppose that the kernel [0, 1] × [0, 1] → R defined by (s, t) 7→ w(s − t) is
continuous and positive definite. Then φ as in (7) is also positive definite.
Proof. By Mercer’s theorem (see Theorem A.1), there exist nonnegative numbers
R1
{λj , j ≥ 0}, orthonormal functions ξj : [0, 1] → R such that 0 w(t − s)ξj (s) ds = λj ξj (t),
and
X
w(t − s) =
λj ξj (t)ξj (s),
j≥0

6

Y.F. Atchadé

and the series converges uniformly and absolutely. It is easy to show that one can inR1
R1
P
terchange integral and sum and write v(t) = 0 w(t − s) ds = j≥0 λj ξj (t) 0 ξj (s) ds,
R1
R1
R1R1
P
2
j≥0 λj ( 0 ξj (t) dt) , and then we get
0 v(t) dt = 0 0 w(t − s) ds dt =
φ(s, t) =

X
j≥0


Z
λj ξj (t) −

1

ξj (t) dt

0


Z
ξj (s) −

0

1


ξj (s) ds .

This expression of φ easily shows that it is positive definite.



The usual approach for showing that the kernel (s, t) 7→ w(s − t) is positive definite
is by showing that the weight function t 7→ w(t) is a characteristic function (or more
generally the Fourier transform of a positive measure) and applying Bochner’s theorem.
This approach shows that A2 holds for the Bartlett function w(x) = (1 − |x|)1(−1,1) (x),
the Parzen function

2
3

if |x| ≤ 12 ,
 1 − 6x + 6|x| ,
w(x) = 2(1 − |x|)3 ,
if 21 ≤ |x| ≤ 1,


0,
if |x| > 1,

and for a number of others weight functions (see, e.g., Hannan [12], pages 278–279 for
details). In the case of the Bartlett function, the kernel φ is given by
φ(s, t) =

2
3

− s(1 − s) − t(1 − t) − |s − t|.

For the Parzen function, we have
v(s) =

3
3
4
+ s ∧ (1 − s) − 2(s ∧ (1 − s)) + (s ∧ (1 − s))
8

and

Z

1

v(t) dt =

0

23
,
40

def

where a ∧ b = min(a, b).
Assumption A2 implies that φ, considered as a linear operator on L2 [0, 1] (φf (s) =
R1
φ(s, t)f (t) dt) is self-adjoint, compact and positive. Therefore, it has only nonnegative
0
eigenvalues, and a countable number of positive eigenvalues. We denote {αj , j ∈ I} the
set of positive eigenvalues of φ (each repeated according to its multiplicity). The index
set I ⊆ {1, 2, . . .} is either finite or I = {1, 2, . . .}. We introduce the random variable Tw
defined as
def

Tw = qP

Z0

2
i∈I αi Zi

Here is the main result.

i.i.d.

where {Z0 , Zi , i ∈ I} ∼ N(0, 1).

Markov Chain Monte Carlo confidence intervals

7

Theorem 2.2. Assume A1–A2, and h ∈ L2 (π). If 0 < σP2 (h) < ∞, then as n → ∞,
w

σn2 → σP2 (h)

X

αi Zi2

and

def

Tn =

n

1 X
w
√
(h(Xk ) − π(h)) → Tw ,
σn n
k=1

i∈I

i.i.d.

where {Zi , i ∈ I} ∼ N(0, 1).
Proof. See Section 3.1.



The theorem implies that the confidence interval
r
σn2
π̂n (h) ± t1−α/2
,
n

(8)

is an asymptotically valid Monte Carlo confidence interval for π(h), where t1−α/2 is the
(1 − α/2)-quantile of the distribution of Tw . These quantiles are intractable in general
but can be easily approximated by Monte Carlo simulation (see Section 2.3).
The assumption that σP2 (h) is finite can be difficult to check. When P is known to
satisfy a drift condition, one can find whole class of functions for which the asymptotic
variance is finite, as the following proposition shows. The proposition uses Markov chain
concepts that have not been defined above, and we refer the reader to Meyn and Tweedie
[21] for details.
Proposition 2.3. Suppose that P is φ-irreducible and aperiodic, with invariant distribution π. Suppose also that there exist measurable functions V, f : X → [1, ∞), constant
b < ∞, and some petite set C ∈ B such that
P V (x) ≤ V (x) − f (x) + b1C (x),

x ∈ X.

(9)

If π(f V ) < ∞, then for all h ∈ Lf , σP2 (h) < ∞.
Proof. This is a well-known result. We give the proof only for completeness.
PWithout any
loss of generality, suppose that π(h) = 0. We recall that σP2 (h) = π(h2 ) + 2 j≥1 hh, P j hi.
R
Since |hh, P j hi| ≤ |h(x)||P j h(x)|π(dx), we obtain
X

Z
X
|hh, P j hi| ≤ |h|f |h(x)|
kP j (x, ·) − π(·)kf π(dx).
j≥0

j≥0

Since P is φ-irreducible and aperiodic, and under the drift condition (9), Meyn and
Tweedie [21], Theorem 14.0.1 implies that there exists a finite constant B such that
P
j
j≥0 kP (x, ·) − π(·)kf ≤ BV (x), x ∈ X. We conclude that
σP2 (h)

≤ 2B|h|f

Z

|h(x)|V (x)π(dx) ≤

2B|h|2f

Z

f (x)V (x)π(dx) < ∞.



8

Y.F. Atchadé

Remark 2. Proposition 2.3 has a number of well-known special cases. The most common
case is when f = λV for some λ ∈ (0, 1), in which case P is geometrically ergodic and
σP2 (h) < ∞ for all h ∈ LV 1/2 . Another important special case is f = V α , for some α ∈
[0, 1). Such drift condition implies that the Markov chain converges at a polynomial rate.
If α ≥ 0.5, then Proposition 2.3 implies that σP2 (h) < ∞ for all h ∈ LV α−0.5 . To see this,
notice that (9) with f = V α , and Jarner and Roberts [14], Lemma 3.5 imply that P V 1/2 ≤
V 1/2 − cV α−1/2 + b1 1C . Since π(V α ) < ∞, the claim follows from Proposition 2.3.

2.2. Example: Metropolis Adjusted Langevin Algorithm for
smooth densities
We give another example where it is possible to check that σP2 (h) < ∞ without geometric
ergodicity. Take X = Rd equipped with the usual Euclidean inner product h·, ·i2 , norm
| · |, and the Lebesgue measure denoted dx. We consider a probability measure π that
has a density with respect to the Lebesgue measure, and in a slight abuse of notation we
use the same symbol to represent π and its density: π(x) = e−u(x) /Z, for some function
u : X → R that we assume is differentiable, with gradient ∇u.
2
Let qσ (x, ·) denotes the density of the Gaussian distribution N(x − σ2 ρ(x)∇u(x), σ 2 Id ),
2
where the term ρ(x) ≥ 0 is used to modulate the drift − σ2 ∇u(x), and σ > 0 is a scaling
constant. We consider the Metropolis–Hastings algorithm that generates a Markov chain
{Xn , n ≥ 0} with invariant distribution π as follows. Given Xn = x, we propose Y ∼
qσ (x, ·). We either “accept” Y and set Xn+1 = Y with probability α(x, Y ), or we “reject”
Y and set Xn+1 = x, where


π(y) qσ (y, x)
α(x, y) = min 1,
.
π(x) qσ (x, y)
def

When ρ(x) = 0, we get the Random Walk Metropolis (RWM), and when ρ(x) = 1, we get
the Metropolis Adjusted Langevin Algorithm (MaLa). However, we are mainly interested
in the case where
τ
def
,
x∈X
(10)
ρ(x) =
max(τ, |∇u(x)|)
for some given constant τ > 0, which corresponds to the truncated MaLa proposed by
Roberts and Tweedie [28]. The truncated MaLa combines the stability of the RWM and
the mixing of the MaLa. It is known to be geometrically ergodic whenever RWM is
geometrically ergodic (Atchadé [2]). However, checking in practice that the truncated
MaLa is geometrically ergodic can be difficult, as this involves checking conditions on
the curvature of the log-density. We show in the next result that if the gradient of the
log-density u is Lipschitz and unbounded then P satisfies a drift condition of the type
(9), and σP2 (h) is guaranteed to be finite for certain functions.

Markov Chain Monte Carlo confidence intervals

9

B1 Suppose that u is bounded from below, continuously differentiable, and ∇u is
Lipschitz, and
lim sup|∇u(x)| = +∞.
|x|→∞

def

Theorem 2.4. Assume B1 and (10). Set V (x) = a + u(x), where a ∈ R is chosen such
that V ≥ 1. Then there exist b, r ∈ (0, ∞) such that
P V (x) ≤ V (x) −

σ2
2
ρ(x)|∇u(x)| + b1{|x|≤r}(x),
4

x ∈ X.

(11)

R
In particular, if u(x)|∇u(x)|e−u(x) dx < ∞, then σP2 (h) < ∞ for all h ∈ Lf , where
f (x) = ρ(x)|∇u(x)|2 .
Proof. See Section 3.2.



Remark 3. This result can be useful in contexts where the log-density u is known to
have a Lipschitz gradient, but is too complicated to allow an easy verification of the
geometric ergodicity conditions.

2.3. On the distribution of the random variable Tw
It is clear that the limiting distribution Tw used for constructing the confidence interval
(8) depends on the choice of w. More research is needed to explain how to best choose w
in this regard. But from the limited simulations done in this paper, we found that weight
functions w with large characteristic exponents lead to heavy-tailed limiting distributions
Tw , and wider confidence intervals. The characteristic exponent of a weight function
w is the largest number r > 0 such that limu→0 |u|−r (1 − w(u)) ∈ (0, ∞). Overall, we
recommend the use of the Bartlett weight function w(u) = (1 − |u|)1(−1,1) (u), which has
characteristic exponent 1, and has behaved very well in the simulations conducted.
Another issue is how to compute the quantiles of Tw . As defined, the distribution of
Tw is intractable in general, as it requires knowing the eigenvalues of φ. But the next
result gives a straightforward method for approximate simulation from Tw .
Proposition 2.5. Let {Zj , 1 ≤ j ≤ N } be i.i.d. standard normal random variables. Then
) def
=
T(N
w

PN

qP P
N
N
i=1

j=1 Zj

i−1 j−1
j=1 φ( N , N )Zi Zj

w

→ Tw

as N → ∞.

Remark 4. As pointed out by a referee, one can also approximately sample from Tw by
i.i.d.
generating X1:N ∼ N(0, 1), and compute TN , with h(x) = x. The approach in Propo2
2
sition 2.5 is similar, but replaces σN
by σ̌N
as defined in (17). By Lemma 3.4, the two
approaches are essentially equivalent.

10

Y.F. Atchadé
Table 1. Approximations of t such that P(Tw >
t) = α/2

w(u) = (1 − u2 )+
Parzen
Bartlett

α = 10%

α = 5%

15.49 (0.06)
4.11 (0.01)
3.77 (0.005)

31.21 (0.19)
5.64
4.78 (0.01)

Proof of Proposition 2.5. Let {0, αj , j ∈ I} be the eigenvalues of φ, with associated
eigenfunctions {Ψ0 , Ψj , j ∈ I} (Ψ0 ≡ 1). By Mercer’s theorem (see Theorem 14 in the
Appendix),


N X
N
X
X
i−1 k−1
Zi Zk = N
αj
,
φ
N
N
i=1
j∈I

k=1



N
1 X
i−1
√
Zi
Ψj
N
N i=1

!2

.

Hence,
)
T(N
w

√ PN
1/ N i=1 Ψ0 ((i − 1)/N )Zi
= qP
.
√ PN
2
N
Ψ
((i
−
1)/N
)Z
)
α
(1/
j
i
j
i=1
j∈I

P
i−1
It is an application of Lemma 3.3 that as N → ∞, { √1N N
i=1 Ψ0 ( N )Zi ,
P
N
i−1
√1
i=1 Ψj ( N )Zi , j ∈ I} converges weakly to {Z0 , Zj , j ∈ I}. The result then follows
N
from the continuous mapping theorem.

We use Proposition 2.5 to approximately simulate Tw for the function w(u) = (1 −
u2 )1(−1,1) (u), and for the Bartlett and Parzen functions. Table 1 reports the 95% and
(N )

97.5% quantiles, computed based on 10 000 independent samples of Tw , with N = 3000.
We replicate these estimates 50 times to evaluate the Monte Carlo errors reported in
parenthesis.
√
As explained in Example 1, in the case w(u) = (1 − u2 )1(−1,1) (u), Tw = 6T1 , where
Tν denotes the student’s distribution with ν degree of freedom; thus, is this case we can
compute accurately the quantiles. In particular, the 95% and 97.5% quantiles are 15.465
and 31.123, respectively.

2.4. Rate of convergence of σn2
An interesting question is understanding how the lag-window estimators σn2 and σb2n
compare. On one hand, the asymptotic behavior of σb2n is better understood. In the
stationary case, the best rate of convergence of σb2n towards σP2 (h) is n−q/(1+2q) (see,
e.g., Parzen [24], Theorem 5A–B), where q is the largest number q ∈ (0, r] such that
P
q
j
j≥1 j γj (h) < ∞, where γj (h) = hh, P hi, and r is the characteristic exponent of w.

Markov Chain Monte Carlo confidence intervals

11

This optimal rate is achieved by choosing bn ∝ n1/(1+2q) . Hence, the optimal rate in the
case of a geometrically ergodic Markov chain is n−r/(1+2r) . However, it is well documented
(see, e.g., Newey and West [23]) that the finite sample properties of σb2n are very sensitive
to the actual constant in bn ∝ n1/(1+2q) , and some tuning is often required in practice.
On the other hand, the fixed-b framework has the advantage that it requires no tuning,
since bn = n. Furthermore, we establish in this section that σn2 has a better convergence
rate. Reversibility plays no role in this discussion. We further simplify the analysis by
assuming that P satisfies a geometric ergodicity assumption:
(G) There exists a measurable function V : X → [1, ∞) such that π(V ) < ∞, and for
all β ∈ (0, 1],
kP n (x, ·) − π(·)kV β ≤ Cρn V β (x),

n ≥ 0, x ∈ X.

(12)

Denote Lip1 (R) the set of all bounded Lipschitz functions f : R → R such that
def

|f (x) − f (y)|
≤ 1.
|x − y|
x6=y

|f |Lip = sup

For P, Q two probability measures on R, we define
def

d1 (P, Q) =

sup
f ∈Lip1 (R)

Z

f dP −

Z

f dQ .

d1 (P, Q) is the Wasserstein metric between P, Q. An upper bound on d1 (Pn , P ) gives a
Berry–Esseen-type bound on the rate of weak convergence of Pn to P . In a slight abuse
of notation, if X, Y are random variables, and X ∼ P and Y ∼ Q, we shall also write
d1 (X, Y ) to mean d1 (P, Q).
Theorem 2.6. Suppose that A2 and (G) hold. Suppose also that I is finite. For δ ∈
[0, 1/4), let h ∈ LV δ be such that π(h) = 0, and σP2 (h) = 1. Then
log(n)
d1 (σn2 , χ2 ) . √
n

as n → ∞,

(13)

P
where χ2 = i∈I αi Zi2 , {Zi , i ∈ I} are i.i.d. N(0, 1), and {αi , i ∈ I} is the set of positive
eigenvalues of φ.
Proof. See Section 3.3.



Remark 5. The assumption that I is finite is mostly technical and it seems plausible
that this result continues to hold without that assumption. For example, I is finite for
the kernel w(u) = (1 − u2 )1(−1,1) (u).

12

Y.F. Atchadé

2.5. A simulation example
This section illustrates the finite sample behavior of the fixed-b confidence interval procedure. We will compare the fixed-b procedure and the standard confidence interval
procedure based on σb2n (using a Gaussian limit). As example, we consider the posterior distribution of a logistic regression model, and use the Random Walk Metropolis
algorithm (Robert and Casella [26]).
Let X = Θ = Rd equipped with its Borel sigma-algebra, and π be absolutely continuous
w.r.t. the Lebesgue measure dθ with density still denoted by π. We write |θ| for the
Euclidean norm of θ. Let qΣ denotes the density of the normal distribution N(0, Σ) on
Θ with covariance matrix Σ. The Random Walk Metropolis algorithm (RWMA) is a
popular MCMC algorithm that generates a Markov chain with invariant distribution π
and transition kernel given by
Z
PΣ (θ, A) = 1A (θ) + α(θ, θ + z)(1A (θ + z) − 1A (θ))qΣ (z) dz,
θ ∈ Θ, A ∈ B(Θ),
X

def

where 1A denotes the indicator function, and α(θ, ϑ) = min(1, π(ϑ)
π(θ) ) is the acceptance
probability.
We assume that π is the posterior distribution from a logistic regression model. More
precisely, we assume that we have binary responses yi ∈ {0, 1}, where
yi ∼ B(p(x′i θ)),

i = 1, . . . , n,

and xi ∈ Rd is a vector of covariate, and θ ∈ Rd is the vector of parameter.
B(p) denotes
ex
the Bernoulli distribution with parameter p ∈ (0, 1), and p(x) = 1+e
x is the cdf of the
logistic distribution. Let X ∈ Rn×d denote the matrix with ith row x′i . Let ℓ(θ|X) denotes
the log-likelihood function of the model. We assume a Gaussian prior N(0, s2 Id ) for θ,
with s = 20. The posterior distribution of θ then becomes
2

π(θ|X) ∝ eℓ(θ|X)e−1/(2s

)|θ|2

.

It is known that for this target distribution the RWM is geometrically ergodic (see,
e.g., Atchadé [3], Section 5.2). Therefore, for all polynomial functions Theorem 2.2 holds.
It is also known that with an appropriate choice of bn , σb2n converges in probability to
σP2 (h) (see, e.g., Atchadé [3], Theorem 4.1, and Corollary 4.1). So we will compare the
fixed-b confidence intervals and the classical confidence intervals based on σb2n .
We simulate a Gaussian dataset with n = 250, d = 15, and simulate the components of
the true value of β from a U(−10, 10). We first run the adaptive chain for 106 iterations
and take the sample posterior mean of β as the “true” posterior mean. We focus on the
coefficient β1 . Each sampler is run for 30 000 iterations, with no burn-in period. For the
RMW, we use a covariance matrix Σ = cI15 , where c is chosen such that the acceptance
probability in stationarity is about 30%, obtained from a preliminary run.
From each sampler, we compute the fixed-b 95% confidence interval, and a classical
95% confidence interval. To explore the range of behavior of the classical procedure,

Markov Chain Monte Carlo confidence intervals

13

Table 2. Coverage probability and half-length for
fixed-b confidence intervals

w(u) = (1 − u2 )+
Parzen
Bartlett

Coverage

Half-length

0.945 ± 0.03
0.94 ± 0.03
0.955 ± 0.03

0.10 ± 0.01
0.03 ± 0.002
0.02 ± 0.001

we use bn = nδ for different values of δ ∈ (0, 1). To estimate coverage probability and
half-length of these confidence intervals, K = 200 replications are performed. The result
is summarized in Table 2 for the fixed-b procedure, and in Figure 1 for the classical
procedure.
We see from the results that using bn = n gives very good coverage, except for the
choice w(u) = (1 − u2 )+ , which generates significantly wider intervals. This is somewhat
expected given the very heavy tail of the limiting distribution. The result also shows that
the confidence interval procedure based on σb2n works equally well when bn is carefully
chosen, but can perform poorly otherwise.
We also test the conclusion of Theorem 2.6 by comparing the finite sample convergence
rate of the two confidence interval procedures. Here, we use only the Bartlett function.
For the standard procedure, we use the best choice of δ (δ ≈ 0.66), as given by the
previous simulation. We compute the confidence intervals after MCMC runs of length
n, where n ∈ {100, . . . , 104 }. Each run is repeated 30 times to approximate the coverage
probabilities and interval lengths. The result is plotted on Figure 2, and is consistent
with Theorem 2.6 that the fixed-b procedure has faster convergence. The price to pay is
a (slightly) wider interval length as seen on Figure 2.

3. Proofs
3.1. Proof of Theorem 2.2
Let φ as in (7). Assumption A2 and Mercer’s theorem implies that the kernel φ has
a countable number of positive eigenvalues {αi , i ∈ I} with associated eigenfunctions
{Ψj , j ∈ I} such that
X
φ(s, t) =
αj Ψj (s)Ψj (t),
(s, t) ∈ [0, 1] × [0, 1],
(14)
j∈I

R1
where the convergence of the series is uniform on [0, 1] × [0, 1]. Since 0 φ(s, t) dt = 0,
0 is also an eigenvalue of φ with eigenfunction Ψ0 (x) ≡ 1. Hence, we define Ī = {0} ∪ I,
α = {αj , j ∈ Ī}, with α0 = 0, and ℓ2 (α) the associated Hilbert space of real
qnumbers
P 2
P
2
sequences {xj , j ∈ Ī} such that j xj < ∞, equipped with the norm kxkα =
j αj xj
def P
and the inner product hx, yiα = j αj xj yj . We will need the differentiability of the

14

Y.F. Atchadé

Figure 1. Coverage probability and confidence interval half-length for parameter β1 for different
values of δ using σb2n , and bn = nδ . The dashed line is the 95% confidence band estimated from
200 replications.

Figure 2. Coverage probability and confidence interval half-length for parameter β1 as function
of number of MCMC iterations. The square-line corresponds to using σn2 .

Markov Chain Monte Carlo confidence intervals

15

eigenfunction Ψj . This is given by Kadota’s theorem (Kadota [16]). Under the assumption
that w is continuously twice differentiable, the eigenfunctions Ψj , j ∈ I are continuously
differentiable (with derivative Ψ′ ) and
X
∂2
φ(s, t) =
αj Ψ′j (s)Ψ′j (t),
∂s∂t

(s, t) ∈ [0, 1] × [0, 1],

j∈I

(15)

where again the convergence of the series is uniform on [0, 1] × [0, 1]. The expansions (14)
and (15) easily imply that
X
j∈I

αj < ∞,

sup

X

t∈[0,1] j∈I

sup

X

t∈[0,1] j∈I

2

αj |Ψj (t)| < ∞ and
(16)

2

αj |Ψ′j (t)| < ∞.

It is easy to check that σn2 can also be written as
σn2



n
n
i−j
1 XX
w
(h̄(Xi ) − πn (h̄))(h̄(Xj ) − πn (h̄))
=
n i=1 j=1
n
=



n
n  
i−j
1 XX
w
− vn,i − vn,j + un h̄(Xi )h̄(Xj ),
n i=1 j=1
n

Pn

Pn Pn

w( i−j
n ). Notice that vn,i is a RieR 1=1
mann sum approximation of v(i/n), where v(t) = 0 w(t − u) du, and un approximates
R 1R 1
R1
0 0 w(t − u) du dt = 0 v(t) dt. In view of this, we introduce
where vn,i = n−1

def

σ̌n2 =

i−ℓ
ℓ=1 w( n ),

and un = n−2

i=1
def





 Z 1

n
n  
1 XX
i−j
i−1
j−1
w
−v
−v
+
v(t) dt h̄(Xi )h̄(Xj )
n i=1 j=1
n
n
n
0

!2 (17)




n X
n
n
X
X
X
1
i−1
1
i−1 j −1
=
φ
Ψℓ
h̄(Xi )h̄(Xj ) =
h̄(Xi ) .
,
αℓ √
n i=1 j=1
n
n
n i=1
n
ℓ∈I

The last equality uses the Mercer’s expansion for φ as given in (14). This implies that
Tn =

Pn

i=1 h̄(Xi )

√
σn n

= qP

√ Pn
1/(σP (h) n) i=1 Ψ0 ((i − 1)/n)h̄(Xi )
.
√ Pn
2
2
2
2
i=1 Ψℓ ((i − 1)/n)h̄(Xi )) + (σn − σ̌n )/σP (h)
ℓ∈I αℓ (1/(σP (h) n)

16

Y.F. Atchadé

Hence, the proof of the theorem boils down to the limiting behavior of the ℓ2 (α)-valued
process
)
(


n
X
i−1
1
√
h̄(Xi ), j ∈ Ī ,
Ψj
n
σP (h) n i=1
i−1
1 √ Pn
and the remainder (σn2 − σ̌n2 ). In Lemma 3.5, we show that { σP (h)
i=1 Ψℓ ( n )h̄(Xi ),
n
ℓ ∈ Ī} converges weakly to {Zℓ , ℓ ∈ Ī}, and that σn2 − σ̌n2 converges in probability to zero.
This is done first in the stationary case in Lemmas 3.3–3.4, and in the nonstationary case
in Lemma 3.5. Hence, the theorem follows by applying Slutszy’s theorem and the continuous mapping theorem. Everything rely on a refinement of the martingale approximation
of Kipnis and Varadhan [19] that we establish first in Lemma 3.2.
3.1.1. Martingale approximation for Markov chains
Throughout this section, unless stated otherwise, {Xn , n ≥ 0} denotes a stationary reversible Markov chain with invariant distribution π and transition kernel P , and we
def
fix h ∈ L20 (π). We denote Fn = σ(X0 , . . . , Xn ). We introduce the probability measure
π̄(dx, dy) = π(dx)P (x, dy) on X × X, and we denote L2 (π̄) the associated L2 -space with
def RR
norm |||f |||2 =
|f (x, y)|2 × π(dx)P (x, dy). For ε > 0, define
def

Uε (x) =

X
j≥0

1
P j h(x),
(1 + ε)j+1

def

Gε (x, y) = Uε (y) − P Uε (x).

Since P is a contraction of L20 (π), it is clear that Uε ∈ L2 (π), and Gε ∈ L2 (π̄). Furthermore, for all ε > 0,
kUε k ≤ ε−1 khk and |||Gε ||| ≤ 2kUε k.

(18)

When σP2 (h) < ∞ a stronger conclusion is possible, and this is the key observation
made by Kipnis and Varadhan [19], Theorem 1.3. We summarize their result as follows.
Lemma 3.1 (Kipnis and Varadhan [19]). Suppose that h ∈ L20 (π), and σP2 (h) < ∞.
Then for any sequence {εn , n ≥ 0} of positive numbers such that limn εn = 0,
lim

n→∞

√
εn kUεn k = 0.

Furthermore, there exists G ∈ L2 (π̄), with
σP2 (h) = |||G|||2 , and limn |||Gεn − G||| = 0.

R

P (x, dz)G(x, z) = 0 (π-a.e.) such that

For n ≥ 1, define the process
⌊nt⌋
X
1
√
Bn (t) =
G(Xi , Xi−1 ),
σP (h) n i=1

0 ≤ t ≤ 1,

Markov Chain Monte Carlo confidence intervals

17

and let {B(t), 0 ≤ t ≤ 1} denotes the standard Brownian motion. It is an easy consequence of Lemma 3.1 that {G(Xi , Xi−1 ), 1 ≤ i ≤ n} is a stationary martingale difference sequence with finite variance. Therefore, by the weak invariance principle for staw
tionary martingales, Bn → B in D[0, 1] equipped with the Skorohod metric. In Corollary 1.5, [19], it is shown that the Markov chain {Xn , n ≥ 0} inherits this weak invariance principle. For the purpose of this paper, we need some refinements of this result.
def
Let {an,k , 0 ≤ k ≤ n} be a sequence of real numbers. Set |an |∞ = sup0≤k≤n |an,k |, and
def Pn
|an |tv = k=1 |an,k − an,k−1 |.
Lemma 3.2. Let h ∈ L20 (π) be such that σP2 (h) < ∞.
(1) If |an |∞ + |an |tv is bounded in n, then
n
X

an,i−1 h(Xi ) =

n
X

an,i−1 G(Xi , Xi−1 ) + Rn ,

(19)

i=1

i=1

where n−1 E(|Rn |2 ) → 0 as n → ∞.
(2) If f : [0, 1] → R is a continuously differentiable function, then
R1
h(Xi ) converges weakly to 0 f (t) dB(t), as n → ∞.

1 √
σP (h) n

Pn

i−1
i=1 f ( n ) ×

def Pn
Proof. Set Sn = i=1 an,i−1 h(Xi ). The function Uε satisfies (1 + ε)Uε (x) − P Uε (x) =
h(x), π-a.e. x ∈ X. This is used to write

an,k−1 h(Xk ) = an,k−1 (εUε (Xk ) + Uε (Xk ) − P Uε (Xk ))
= an,k−1 εUε (Xk ) + an,k−1 (Uε (Xk ) − P Uε (Xk−1 ))
+ (an,k−1 P Uε (Xk−1 ) − an,k P Uε (Xk )) + (an,k − an,k−1 )P Uε (Xk ).
It follows that
Sn = ε

n
X

an,k−1 Uε (Xk ) +

k=1

+

n
X

k=1

n
X

an,k−1 G(Xk , Xk−1 )

k=1

an,k−1 (Gε (Xk , Xk−1 ) − G(Xk , Xk−1 ))

+ (an,0 P Uε (X0 ) − an,n P Uε (Xn )) +

n
X

k=1

(an,k − an,k−1 )P Uε (Xk ),

which is valid for any ε > 0. In particular with ε = εn = 1/n, we have
Sn =

n
X

k=1

an,k−1 G(Xk , Xk−1 )+

n
X

k=1

an,k−1 (Gεn (Xk , Xk−1 )−G(Xk , Xk−1 ))+Rn(1) +Rn(2) +Rn(3) ,

18

Y.F. Atchadé

where
def

Rn(1) = εn

n
X

an,k−1 Uεn (Xk ),

k=1
def

Rn(3) =

n
X

k=1

def

Rn(2) = (an,0 P Uεn (X0 ) − an,n P Uεn (Xn ))

and

(an,k − an,k−1 )P Uεn (Xk ).

By stationarity and the martingale property,
" n
!2 #
n
X
1
1X 2
E
an,k−1 → 0,
an,k−1 (Gε (Xk , Xk−1 ) − G(Xk , Xk−1 ))
= |||Gεn − G|||2
n
n
k=1

k=1

using Lemma 3.1, and the assumption on an . The other remainders are also easily dealt
with.
n

√
1 X
1
2
√ E1/2 (|Rn(3) |2 ) ≤ √
|an,k − an,k−1 |E1/2 (|P Uεn (Xk )| ) = εn kUεn k|an |tv → 0,
n
n
k=1

using Lemma 3.1 and the assumption on an . Similarly,
√
1
2
√ E1/2 (|Rn(2) | ) ≤ 2|an |∞ εn kUεn k → 0 and
n
n

√
1X
1
√ E1/2 (|Rn(1) |2 ) ≤ εn kUεn k
|an,k−1 | → 0.
n
n
k=1

This proves part (1) of the lemma. For part (2), we use part (1) with an,i = f (i/n) to
conclude that


n
X
i−1
1
√
f
h(Xi )
n
σP (h) n i=1


n
X
i−1
1
√
f
G(Xi , Xi−1 ) + op (1)
σP (h) n i=1
n
Z 1
=
f (t) dBn (t) + op (1),

=

0

where An = op (1) means that An converges in probability to zero as n → ∞. To conclude
R1
R1
the proof, it suffices to show that 0 f (t) dBn (t) converges weakly to 0 f (t) dB(t). This
follows from the weak convergence continuous mapping theorem by noticing that B has
R1
continuous sample path (almost surely), and the map D[0, 1] → R, x 7→ 0 f (t) dx(t) is
R1
continuous at all points x0 ∈ C[0, 1], where the integral 0 f (t) dx(t) is understood as

Markov Chain Monte Carlo confidence intervals

19

a Riemann–Stietjes integral. To see the continuity, take {xn } a sequence of elements
in D[0, 1] that converges to x0 ∈ C[0, 1] in the Skorohod metric. Since x0 ∈ C[0, 1], the
R1
sequence {xn } converges to x0 in C[0, 1] as well. By integration by part, 0 f (t) dxn (t) =
R1
f (1)xn (1) − f (0)xn (0) − 0 xn (t)f ′ (t) dt, and
Z

1

0

f (t) dxn (t) −

Z

0



1

f (t) dx0 (t) ≤ |xn − x0 |∞ 2|f |∞ +

Z

0

1

′



|f (t)| dt → 0,

as n → ∞.



Lemma 3.3. Let h ∈ L20 (π) be such that σP2 (h) < ∞. Define
)
(


Z 1

n
X
i
−
1
1
def
def
(n)
√
Ψj
h̄(Xi ), j ∈ Ī
and Z =
Z
=
Ψj (t) dB(t), j ∈ Ī .
σP (h) n i=1
n
0
Then as n → ∞, Z (n) converges weakly to Z in ℓ2 (α).
w

Proof. We need to show that for all u ∈ ℓ2 (α), hZ (n) , uiα → hZ, uiα , and that {Z (n) } is
tight.
P
i−1
1 √ Pn
For u ∈ ℓ2 (α), hZ (n) , uiα = σP (h)
j αj uj Ψj (t).
i=1 fu ( n )h̄(Xi ), where fu (t) =
n
From basic results in calculus, it follows from Kadota’s theorem that fuR is continuously
w
1
differentiable on [0, 1]. Hence, by Lemma 3.2, part (2), hZ (n) , uiα → 0 fu (t) dB(t) =
hu, Ziα . To show that {Z (n) } is tight, it suffices to show that
!
∞
X
2
(n)
lim sup E
(20)
hZ , ej iα = 0.
N →∞ n≥1

j=N

We have
2

E(hZ (n) , ej iα ) =

 


n
n
αj X X
k−1
i−1
Ψ
π(hP |i−k| h)
Ψ
j
j
σP2 (h)n i=1
n
n

αj
= 2
σP (h)n

k=1

Z

1

n X
n
X

−1 i=1 k=1

Ψj



 

k − 1 |i−k|
i−1
Ψj
λ
µh (dλ).
n
n

By Fubini’s theorem, for N ≥ 1,
!
 


Z 1X
∞
n X
n X
∞
X
1
k − 1 |i−k|
i−1
2
(n)
E
hZ , ej iα = 2
Ψj
λ
µh (dλ).
αj Ψj
σP (h)n −1 i=1
n
n
j=N

k=1 j=N

P
Let ε > 0. By uniform convergence of the series j αj Ψj (s)Ψj (t), we can find N0 such
P
that for any N ≥ N0 and for all s, t ∈ [0, 1], | ℓ≥N αℓ Ψℓ (t)Ψℓ (s)| ≤ ε. So that for all

20

Y.F. Atchadé

n ≥ 1,
E

∞
X

ℓ=N

hZ

(n)

2
, eℓ iα

!

ε
≤ 2
σP (h)n

Z

1

n X
n
X

|i−j|

λ

−1 i=1 j=1

ε
µh (dλ) ≤ 2
σP (h)

Z

1

−1

1+λ
µh (dλ) = ε,
1−λ

since ε > 0 is arbitrary, this proves (20).



Lemma 3.4. Let h ∈ L20 (π) be such that σP2 (h) < ∞. Then as n → ∞, E(|σn2 − σ̌n2 |) =
O(1/n). Hence σn2 − σ̌n2 converges in probability to 0, as n → ∞.
Proof. Comparing the expression of σn2 and σ̌n2 , we see that
!2
n
1 X
√
h̄(Xi )
v(t) dt
n i=1
0
!
!


n
n 
i−1
1 X
1 X
√
vn,i − v
h̄(Xi )
−2 √
h̄(Xi ) .
n
n i=1
n i=1


Z
2
2
σn − σ̌n = un −

1



(21)

Pn
Since the sequence E[( √1n i=1 h̄(Xi ))2 ] converges to the finite limit σ 2 (h) by assumption,
it is bounded, and there exists a finite constant c1 such that
E(|σn2 − σ̌n2 |)
≤

c21

un −

Z

1

0

2c1 1/2
v(t) dt +
E
n

"

!2 #



n
1 X
i−1
√
n vn,i − v
h(Xi )
.
n i=1
n

def

−1
Set an,0 = 0, an,i = n(vn,i − v( i−1
n )). We recall that vn,i = n
R1
0 w(t − u) du, and write

an,i = n

=n

n Z
X
ℓ=1

i−ℓ
ℓ=1 w( n ),

and v(t) =




 
i−1
i−1 ℓ−1
−w
−
− u du
w
n
n
n
(ℓ−1)/n

n Z
X
ℓ=1

Pn

ℓ/n

ℓ/n

(ℓ−1)/n



ℓ−1
−u
n

Z

0

1

w

′





ℓ−1
i−1 ℓ−1
dt du.
−
−t u−
n
n
n

Using this expression, it is easy to show that |an |∞ ≤ |w′ |∞ /2. And since w is of class
2
C
on w′ using the above expression shows that |an |tv = |an,1 | +
Pn, a mean-value theorem
′
′′
i=2 |an,i − an,i−1 | ≤ (|w |∞ + |w |∞ )/2. We are then in position to apply Lemma 3.2(1)
to obtain
"
!2 #
n
1 X
√
E
= O(1).
an,i h̄(Xi )
n i=1

Markov Chain Monte Carlo confidence intervals

21

P P
By similar arguments as above, and since un = n−2 ni=1 n=1 w( i−j
n ) is a Riemann sum
R1
R1
1
approximation of 0 v(t) dt, we obtain that |un − 0 v(t) dt| = O( n ). In conclusion,
E(|σn2

− σ̌n2 |)

 
1
.
=O
n

(22)


Lemma 3.5. Assume A1. Suppose that the Markov chain {Xn , n ≥ 0} starts at X0 = x
for x ∈ X such that (1) holds. Let h ∈ L20 (π) be such that σP2 (h) < ∞. Then as n → ∞,
w
σn2 − σ̌n2 converges in probability to zero, and Z (n) → Z in ℓ2 (α).
Proof. Ergodicity is equivalent to the existence of a successful coupling of the Markov
chain and its stationary copy. More precisely, we can construct a process {(Xn , X̃n ), n ≥
0} such that {Xn , n ≥ 0} is a Markov chain with initial distribution δx and transition
kernel P , {X̃n , n ≥ 0} is a Markov chain with initial distribution π and transition kernel
P , and there exists a finite (coupling) time τ such that Xn = X̃n for all n ≥ τ . For a
proof of this result, see for instance Lindvall [20], Theorem 14.10; see also Roberts and
Rosenthal [27], Proposition 28. We use a wide “tilde” to denote quantities computed
from the stationary chain {X̃n , n ≥ 0}.
Since Xn = X̃n for all n ≥ τ , and in view of the expression of σn2 − σ̌n2 given in (21),
2
2
it is straightforward to check that σn2 − σ̌n2 − (σ^
n − σ̌n ) converges to zero in probability.
g
(n)
(n)
The convergence of kZ − Z k is handled similarly.
α

(n) k2
kZ (n) − Zg
α

=

X

αℓ

X

αℓ

ℓ∈I

=

k=1

ℓ∈I

τ
≤
n

!2
 
n
1 X
k
√
(h(Xk ) − h(X̃k ))
Ψℓ
n
n

!2
 
τ −1
1 X
k
√
(h(Xk ) − h(X̃k ))
Ψℓ
n
n
k=1



sup

X

t∈[0,1] ℓ∈I

αℓ |Ψℓ (t)|

2

 X
τ

(h(Xk ) − h(X̃k ))

k=1

2

!

,

which converges almost surely to zero, given (16), and since τ is finite almost surely. 

3.2. Proof of Theorem 2.4
def

Since u is bounded from below, we can choose a = 1 − inf x∈X u(x) such that V (x) =
2
a + u(x) ≥ 1. Let qσ (x, y) be the density of the proposal N(x − σ2 ρ(x)∇u(x), σ 2 Id ), and

22

Y.F. Atchadé
def

define R(x) = {y ∈ Rp : α(x, y) < 1}. We have
Z
P V (x) − V (x) = α(x, y)(V (y) − V (x))qσ (x, y) dy
=

Z

R(x)

+

Z

[α(x, y) − 1](V (y) − V (x))qσ (x, y) dy

(23)

(V (y) − V (x))qσ (x, y) dy.

Since ∇u is Lipschitz, with Lipschitz constant L, say, we have by Taylor expansion
V (y) − V (x) ≤ h∇u(x), y − xi2 +

L
|y − x|2 .
2

Integrating both sides, and using the fact that ρ(x)|∇u(x)| ≤ τ , we get


Z
L σ4
σ2
2
2
ρ(x)2 |∇u(x)| + dσ 2
(V (y) − V (x))qσ (x, y) dy ≤ − ρ(x)|∇u(x)| +
2
2 4
(24)


σ2
L τ 2 σ4
2
2
≤ − ρ(x)|∇u(x)| +
+ dσ .
2
2
4
We also have
π(y) qσ (y, x)
π(x) qσ (x, y)

2
2
σ2
σ2
1
1
.
= exp V (x) − V (y) − 2 x − y + ρ(y)∇u(y) + 2 y − x − ρ(x)∇u(x)
2σ
2
2σ
2
If y ∈ R(x), we necessarily have

π(y) qσ (y,x)
π(x) qσ (x,y)

< 1, which translates to
2

V (y) − V (x) > −

2

1
σ2
σ2
1
x − y + ρ(y)∇u(y) + 2 y − x − ρ(x)∇u(x) .
2
2σ
2
2σ
2

Hence, if y ∈ R(x),
[α(x, y) − 1](V (y) − V (x))

2
2
σ2
σ2
1
1
≤ [α(x, y) − 1] − 2 x − y + ρ(y)∇u(y) + 2 y − x − ρ(x)∇u(x)
2σ
2
2σ
2

2
σ
2
2
ρ2 (y)|∇u(y)| − ρ2 (x)|∇u(x)|
= [1 − α(x, y)]
8

2
− 2 hy − x, ρ(x)∇u(x) + ρ(y)∇u(y)i
σ

Markov Chain Monte Carlo confidence intervals


σ 2 2 4τ
τ + 2 |y − x| .
≤
8
σ

23

Hence,
Z

σ2 τ 2 τ
[α(x, y) − 1](V (y) − V (x))qσ (x, y) dy ≤
+
8
2
R(x)

r

dσ 2 +

σ2 τ 2
.
2

(25)

We combine (23)–(25) to conclude that
P V (x) − V (x) ≤ −

σ2
2
ρ(x)|∇u(x)| + K,
2

q
2 4
2 2
2 2
def 2
where K = L2 ( τ 4σ + dσ 2 ) + σ 8τ + τ2 dσ 2 + σ 2τ . Since f (x) = σ2 ρ(x)|∇u(x)|2 is continuous and f (x) → ∞, as kxk → ∞ by assumption, the results follow readily.

3.3. Proof Theorem 2.6
We follow Dedecker and Rio [7], Theorem
2.1. With the geometric ergodicity assumption,
Pn
the martingale approximation to i=1 h(Xi ) can be constructed more explicitly than in
Lemmas 3.1 and 3.2. Define
X
g(x) =
P j h̄(x),
x ∈ X.
j≥0

By the geometric ergodicity assumption, g is well-defined and belongs to LV δ . Then we
def
define D0 = 0, and Dk = g(Xk ) − P g(Xk−1 ), k ≥ 1. It is easy to see that {Dk , k ≥ 0}
is a martingale-difference sequence with respect to the natural filtration of {Xn , n ≥ 0}.
Using this martingale, we define
def
σ̄n2 =

X

αℓ

ℓ∈I

def

and we recall that σ̌n2 =
σn2

=

X
ℓ∈I

αℓ

P

ℓ∈I αℓ (

√1
n


 !2
n
1 X
i−1
√
Ψℓ
Di ,
n i=1
n
Pn

i−1
2
i=1 Ψℓ ( n )h(Xi )) .

Hence,


 !2
n
i−1
1 X
√
Ψℓ
Di + (σn2 − σ̌n2 ) + (σ̌n2 − σ̄n2 ).
n i=1
n

Although the martingales are constructed differently, the argument in Lemma 3.4 carries
through and shows that E(|σn2 − σ̌n2 |) = O(1/n). The proof is similar to the proof of

24

Y.F. Atchadé

√
Lemma 3.4 and is omitted. Also E(|σ̌n2 − σ̄n2 |) = O(1/ n). To see this, use the Cauchy–
Schwarz inequalities for sequences in ℓ2 (α) and for random variables to write
E(|σ̌n2 − σ̄n2 |)
"
!
!#




n
n
X
1 X
1 X
i−1
i−1
√
=E
(h(Xi ) − Di )
(h(Xi ) + Di )
αℓ √
Ψℓ
Ψℓ
n
n
n i=1
n i=1
ℓ∈I

"(

≤E

X

≤

αℓ

ℓ∈I

×
(

X

(

X

×

αℓ

ℓ∈I

"

αℓ E

ℓ∈I

(

!2 )1/2


n
i−1
1 X
√
(h(Xi ) − Di )
Ψℓ
n
n i=1

X

!2 #)1/2


n
i−1
1 X
√
Ψℓ
(h(Xi ) − Di )
n
n i=1

"

αℓ E

ℓ∈I

!2 )1/2 #


n
i−1
1 X
√
(h(Xi ) + Di )
Ψℓ
n
n i=1

!2 #)1/2


n
i−1
1 X
√
Ψℓ
(h(Xi ) + Di )
.
n
n i=1

By the martingale approximation, we have
n
X

Ψℓ

i=1






n−1
i−1
(h(Xi ) − Di ) = Ψℓ (0)P g(X0 ) − Ψ
P g(Xn )
n
n





n 
X
i−1
i−2
Ψℓ
+
− Ψℓ
P g(Xi−1 ).
n
n
i=2

The details of these calculations can be found for instance in [4], Proposition A1. It is
then easy to show that
X
ℓ∈I

"

αℓ E

n
X

Ψℓ

i=1



!2 #

i−1
(h(Xi ) − Di )
n



X
X
2
2
′
αℓ |Ψℓ (t)| |h|2V δ .
αℓ |Ψℓ (t)| + 3 sup
≤ 6 sup
0≤t≤1

ℓ∈I

0≤t≤1

ℓ∈I

For the second term, notice that
n
X
i=1

Ψℓ








n
n
X
X
i−1
i−1
i−1
Ψℓ
Ψℓ
(h(Xi ) + Di ) = 2
Di +
(h(Xi ) − Di ).
n
n
n
i=1
i=1

Markov Chain Monte Carlo confidence intervals

25

Hence, with similar calculations, we obtain
!2 #
" n
X
X i − 1
(h(Xi ) + Di )
Ψℓ
αℓ E
n
i=1
ℓ∈I
X
2
≤ 2|h|2V δ n sup
αℓ |Ψℓ (t)|
0≤t≤1

ℓ∈I



X
X
2
2
+ 6 2 sup
αℓ |Ψℓ (t)| + sup
αℓ |Ψ′ℓ (t)| |h|2V δ .
0≤t≤1

0≤t≤1

ℓ∈I

ℓ∈I

√
Given (16), these calculations show that E(|σ̌n2 − σ̄n2 |) = O(1/ n). We conclude that


 !2

n
X
1
i−1
1 X
2
Ψℓ
Di + O p √ ,
σn =
αj √
n i=1
n
n
ℓ∈I

which implies that
1
d1 (σn2 , χ2 ) . d1 (σ̄n2 , χ2 ) + √ .
n

(26)

Therefore, we only need to focus on the term d1 (σ̄n2 , χ2 ).
P
P
On the Euclidean space RI , we P
define the norms kxk2α = Pi∈I αi x2i , kxk2 = i∈I x2i
and the inner-products hx, yiα = i∈I αi xi yi , and hx, yi = i∈I xi yi . For a sequence
(a1 , a2 , . . .), we use the notation ai:k = (ai , . . . , ak ) (and ai:k is the empty set if i > k). We
introduce new random variables {Zi,j , i ∈ I, 1 ≤ j ≤ n} which are i.i.d. N(0, 1), and set
Pk
def Pk
Sℓ:k = ( j=ℓ Z1j , . . . , j=ℓ ZIj )T ∈ RI , so that
χ

2 dist.

=

X
i∈I

n

αi

1 X
√
Zi,j
n j=1

!2

1
= √ S1:n
n

2

.
α

For 1 ≤ ℓ ≤ k ≤ n, and omitting the dependence on n, we set Bℓ:k as the RI×(k−ℓ+1)
matrix
 
j
,
i ∈ I, ℓ ≤ j ≤ k.
Bℓ:k (i, j) = Ψi
n
By the Mercer’s expansion for φ, we have
σ̄n2 =

X
i∈I

αi

  !2
n
1 X
1
k
√
= √ B1:n D1:n
Dk
Ψi
n
n
n
k=1

2

.
α

For f ∈ Lip1 (R), we introduce the function fα : R|I| → R, defined as fα (x) = f (kxk2α ). As
a matter of telescoping the sums, we have
E[f (σ̄n2 ) − f (χ2 )]

26

Y.F. Atchadé



 
1
1
= E fα √ B1:n D1:n − fα √ S1:n
n
n





n
X
1
1
1
1
E fα √ B1:ℓ D1:ℓ + √ Sℓ+1:n − fα √ B1:ℓ−1 D1:ℓ−1 + √ Sℓ:n
=
n
n
n
n
ℓ=1




n
X
1
1
√
√
B1:ℓ−1 D1:ℓ−1 +
Bℓ:ℓ Dℓ
=
E fα,n,ℓ+1
n
n
ℓ=1


1
1
− fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √ Sℓ:ℓ ,
n
n
where we define
 

1
fα,n,ℓ (x) = E fα x + √ Sℓ:n
n
def

and set fα,n,n+1 (x) = fα (x).

First, we claim that fα,n,ℓ is differentiable everywhere on RI . To prove this, it suffices
to obtain the almost everywhere differentiability of z ∈ RI 7→ fα (x + z) for any x ∈ RI . By
Rademacher’s theorem, f as a Lipschitz function is differentiable almost everywhere on
R. If E is the set of points where f is not differentiable, we conclude that fα is differentiable at all points z ∈
/ {z ∈ RI : kx + zk2α ∈ E}. Now by Ponomarëv [25], Theorem 2, the
Lebesgue measure of the set {z ∈ RI : kx + zk2α ∈ E} is zero, which proves the claim.
As a result, the function x 7→ fα,n,ℓ (x) is differentiable with derivative

 
 
1
1
′
x + √ Sℓ:n , h
.
∇fα,n,ℓ (x) · h = 2E fα x + √ Sℓ:n
n
n
α

By writing this expectation wrt the distribution of x +
∇fα,n,ℓ (x) · h = 2

Z


fα′ (z)hz, hiα exp −

√1 Sℓ:n ,
n

we get


n
(kxk2 − 2hx, zi) µn,ℓ (dz),
2(n − ℓ + 1)

where µn,ℓ is the distribution of √1n Sℓ:n . This implies that fα,n,ℓ is infinitely differentiable
with second derivatives given by
∇(2) fα,n,ℓ (x) · (h1 , h2 )


n
= −2
n−ℓ+1

Z
′
× fα (z)hz, h1 iα hx − z, h2 i exp −


n
2
(kxk − 2hx, zi) µn,ℓ (dz)
2(n − ℓ + 1)
 

 
 r
1
S
Sℓ:n
n
√ ℓ:n
= 2E fα′ x + √ Sℓ:n
, h1
, h2 ,
x
+√
n
n−ℓ+1
n−ℓ+1
n−ℓ+1
α

Markov Chain Monte Carlo confidence intervals

27

which implies after some easy calculations that for h ∈ RI ,
r


n
(2)
2
kxkα .
|∇ fα,n,ℓ (x) · (h, h)| . khk 1 +
n−ℓ+1
Similarly for h ∈ RI ,
|∇(3) fα,n,ℓ (x) · (h, h, h)| .

r

r


n
n
khk3 1 +
kxkα .
n−ℓ+1
n−ℓ+1

(27)

(28)

Now, by Taylor expansion we have




1
1
1
1
fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √ Bℓ:ℓ Dℓ − fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √ Sℓ:ℓ
n
n
n
n


1
1
= √ ∇fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 · (Bℓ:ℓ Dℓ − Sℓ:ℓ )
n
n


1 (2)
1
(3)
+
∇ fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 · [(Bℓ:ℓ Dℓ , Bℓ:ℓ Dℓ ) − (Sℓ:ℓ , Sℓ:ℓ )] + ̺n,ℓ ,
2n
n
where, using (28),
r
r

n
ℓ−1
B1:ℓ−1 D1:ℓ−1
(3)
−3/2
√
|̺n,ℓ | .
1+
n
n−ℓ+1
n−ℓ+1
ℓ−1

α


(kBℓ:ℓ Dℓ k3α + kSℓ:ℓ k3α ).

It follows that
n−1
X
ℓ=1

(3)

E(|̺n,ℓ |) . n−1

n
n
X
X
1
1
√ + n−1/2
. n−1/2 log(n).
ℓ
ℓ
ℓ=1
ℓ=1

By first conditioning on Fℓ−1 , we have




1
E ∇fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 · (Bℓ:ℓ Dℓ − Sℓ:ℓ ) = 0.
n
def 1 (2)
fα,n,ℓ ( √1n B1:ℓ−1 D1:ℓ−1 ),
2∇

Writing Kn,ℓ =

we have


1
√
B1:ℓ−1 D1:ℓ−1 · [(Bℓ:ℓ Dℓ , Bℓ:ℓ Dℓ ) − (Sℓ:ℓ , Sℓ:ℓ )]
∇ fα,n,ℓ+1
n
X ℓ ℓ
X ℓ ℓ
Ψi
Ψj
Kn,ℓ (i, j) −
Ψj
Kn,ℓ (i, j)Zi,ℓ Zjℓ .
Ψi
= Dℓ2
n
n
n
n
i,j
i,j
(2)



Therefore,




1
E ∇(2) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 · [(Bℓ:ℓ Dℓ , Bℓ:ℓ Dℓ ) − (Sℓ:ℓ , Sℓ:ℓ )]|Fℓ−1
n

(29)

28

Y.F. Atchadé
=

X
i,j

Ψi

   
ℓ
ℓ
Ψj
Kn,ℓ+1 (i, j)[E(Dℓ2 |Fℓ−1 ) − δij ],
n
n

where δij = 1 if i = j and zero otherwise. We claim that the proof will be finished if we
show that for all i, j ∈ I, and 1 ≤ ℓ ≤ n,
√

n
.
(30)
n−ℓ+1
Pn
To prove this claim, it suffice to use (30) to show that |n−1 ℓ=1 Ψi ( nℓ )Ψj ( nℓ )E(Kn,ℓ+1 (i, j))| .
P
n
n−1/2 log(n) for i 6= j, and |n−1 ℓ=1 Ψi ( nℓ )Ψj ( nℓ )E(Kn,ℓ+1 (i, j)[E(Dℓ2 |Fℓ−1 ) − 1])| .
n−1/2 log(n) for all i, j ∈ I. To show this, write
1/2

E

2

[(Kn,ℓ (i, j) − Kn,ℓ+1 (i, j)) ] .

   
n−1
1X
ℓ
ℓ
Ψj
E(Kn,ℓ+1 (i, j))
Ψi
n
n
n
ℓ=1
( n−1    )
ℓ
1X
ℓ
=
Ψj
E(Kn,n (i, j))
Ψi
n
n
n
ℓ=1
" ℓ−1    #
n−1
ℓ
ℓ
1X 1X
Ψj
[E(Kn,ℓ (i, j) − Kn,ℓ+1 (i, j))].
+
Ψi
n
n
n
n
ℓ=1

k=1

By the convergence of Riemann sums, | n1
and (30), this implies that

Pn−1
ℓ=1

Ψi ( nℓ )Ψj ( nℓ )| . n−1 . Combined with (27)

!
   
n
n
√ X
ℓ
1 √
log(n)
ℓ
1X
1
Ψj
E(Kn,ℓ+1 (i, j)) ≤
. √ .
Ψi
n+ n
n
n
n
n
k
n
ℓ=1

k=1

For the second term, notice from the definition of Dℓ at the beginning of the proof
that E(Dℓ2 |Fℓ−1 ) − 1 = G(Xℓ−1 ) − π(G), where G(x) = P g 2 (x) − (P g(x))2 . Since h ∈ LV δ
by geometric ergodicity, the solution
for δ < 1/4, G ∈ LV 2δ , and 2δ < 1/2. Therefore,P
of the Poisson equation for G defined as U (x) = j≥0 P j (G(x) − π(G)) is well-defined,
U ∈ LV 2δ , and we have almost surely
U (Xℓ−1 ) − P U (Xℓ−1 ) = E(Dℓ2 |Fℓ−1 ) − 1.
Notice that, since 2δ < 1/2, for any p ≥ 2 such that 2pδ ≤ 1, the geometric ergodicity
assumption (G) implies that supk≥1 E(|U (Xk )|p ) < ∞. Now we use the usual martingale
approximation trick (see, e.g., Atchadé and Cattaneo [4], Proposition A1) to write
   
n−1
ℓ
ℓ
1X
Ψj
E(Kn,ℓ+1 (i, j)[E(Dℓ2 |Fℓ−1 ) − 1])
Ψi
n
n
n
ℓ=1

Markov Chain Monte Carlo confidence intervals
   
1
1
1
Ψj
E(Kn,2 (i, j)U (X0 ))
= Ψi
n
n
n

 

1
1
1
Ψj 1 −
E(Kn,n (i, j)U (Xn−1 ))
− Ψi 1 −
n
n
n
   
n−1 
1X
ℓ
ℓ
+
Ψj
Kn,ℓ+1 (i, j)
E Ψi
n
n
n
ℓ=1

 



ℓ−1
ℓ−1
− Ψi
Ψj
Kn,ℓ (i, j) U (Xℓ−1 ) .
n
n

29

We now use the fact that Ψi Ψj is of class C1 (see Theorem A.1(ii)), (27), and (30) to
conclude that
   
n−1
ℓ
ℓ
1X
Ψj
E(Kn,ℓ+1 (i, j)[E(Dℓ2 |Fℓ−1 ) − 1])
Ψi
n
n
n
ℓ=1

n−1
1
log(n)
1 X 1/2
2
.√ +
E (|Kn,ℓ+1 (i, j) − Kn,ℓ+2 (i, j)| ) . √ .
n n
n
ℓ=1

This proves the claim. It remains to establish (30). Write Eℓ to denote the expectation
operator wrt n−1/2 Sℓ:n . We then have for any h1 , h2 ∈ RI ,
2Kn,ℓ · (h1 , h2 )


1
= ∇(2) fα,n,ℓ √ B1:ℓ−1 D1:ℓ−1 · (h1 , h2 )
n


n
=2
n−ℓ+1
 


 
Sℓ:n
Sℓ:n
1
Sℓ:n
1
′
√ B1:ℓ−1 D1:ℓ−1 + √ , h1
√ , h2
× Eℓ fα √ B1:ℓ−1 D1:ℓ−1 + √
n
n
n
n
n
α




1
n−ℓ
Sℓ
=
∇(2) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √
· (h1 , h2 )
n−ℓ+1
n
n

 

n
1
+
O √ .
n−ℓ+1
n
Therefore,
2(Kn,ℓ − Kn,ℓ+1 ) · (h1 , h2 )


Sℓ
1
· (h1 , h2 )
= ∇(2) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √
n
n

30

Y.F. Atchadé


Bℓ D ℓ
1
· (h1 , h2 )
− ∇(2) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √
n
n



 

1
n
1
Sℓ
1
(2)
· (h1 , h2 ) +
O √
∇ fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √
−
n−ℓ+1
n
n
n−ℓ+1
n

 

1
Sℓ
Sℓ
Bℓ D ℓ
Bℓ D ℓ
= ∇(3) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + t √ + (1 − t) √
· h1 , h2 , √ − √
n
n
n
n
n



 

1
1
n
1
Sℓ
−
· (h1 , h2 ) +
O √ ,
∇(2) fα,n,ℓ+1 √ B1:ℓ−1 D1:ℓ−1 + √
n−ℓ+1
n
n
n−ℓ+1
n

for some t ∈ (0, 1). Using (27) and (28), (30) follows from the above.

Appendix: Mercer’s theorem
We recall Mercer’s theorem below. Part (i) is the standard Mercer’s theorem, and part
(ii) is a special case of a result due to T. Kadota (Kadota [16]).
Theorem A.1 (Mercer’s theorem). (i) Let k : [0, 1] × [0, 1] → R be a continuous positive semidefinite kernel. Then there exist nonnegative numbers {λj , j ≥ 0}, and orthonorR1
mal functions {φj , j ≥ 0}, φj ∈ L2 ([0, 1]), such that 0 k(x, y)φj (y) dy = λj φj (x) for all
x ∈ [0, 1], j ≥ 0, and
lim

sup

n→∞ x,y∈[0,1]

k(x, y) −

n
X

λj φj (x)φj (y) = 0.

j=0

Furthermore, if λj 6= 0, φj is continuous.
(ii) Let k as above. If in addition k is of class C2 on [0, 1] × [0, 1], then for λj 6= 0, φj
is of class C1 on [0, 1] and
lim

sup

n→∞ x,y∈[0,1]

n
X
∂2
λj φ′j (x)φ′j (y) = 0.
k(x, y) −
∂x∂y
j=0

By setting x = y, in both expansions, it follows that
X
2
sup
λj |φj (x)| ≤ sup k(x, x) < ∞
0≤x≤1

(A.1)

0≤x≤1

j≥0

and
sup

X

0≤x≤1 j≥0

2

λj |φ′j (x)| ≤ sup

0≤x≤1

∂2
k(u, v)|u=x,v=x < ∞.
∂u∂v

(A.2)

Markov Chain Monte Carlo confidence intervals

31

Acknowledgments
I am grateful to Matias Cattaneo and Shane Henderson for very helpful discussions.
This work is partly supported by NSF Grants DMS-09-06631, NSF-SES 1229261 and
DMS-12-28164.

References
[1] Asmussen, S. and Glynn, P.W. (2007). Stochastic Simulation: Algorithms and Analysis.
Stochastic Modelling and Applied Probability 57. New York: Springer. MR2331321
[2] Atchadé, Y.F. (2006). An adaptive version for the Metropolis adjusted Langevin algorithm
with a truncated drift. Methodol. Comput. Appl. Probab. 8 235–254. MR2324873
[3] Atchadé, Y.F. (2011). Kernel estimators of asymptotic variance for adaptive Markov chain
Monte Carlo. Ann. Statist. 39 990–1011. MR2816345
[4] Atchadé, Y.F. and Cattaneo, M.D. (2014). A martingale decomposition for quadratic
forms of Markov chains (with applications). Stochastic Process. Appl. 124 646–677.
MR3131309
[5] Chan, K.S. and Geyer, C. (1994). Comment on “Markov chains for exploring posterior
distributions.” Ann. Statist. 22 1747–1758.
[6] Damerdji, H. (1995). Mean-square consistency of the variance estimator in steady-state
simulation output analysis. Oper. Res. 43 282–291. MR1327416
[7] Dedecker, J. and Rio, E. (2008). On mean central limit theorems for stationary sequences.
Ann. Inst. Henri Poincaré Probab. Stat. 44 693–726. MR2446294
[8] Flegal, J.M. and Jones, G.L. (2010). Batch means and spectral variance estimators in
Markov chain Monte Carlo. Ann. Statist. 38 1034–1070. MR2604704
[9] Glynn, P.W. and Iglehart, D.L. (1990). Simulation output analysis using standardized
time series. Math. Oper. Res. 15 1–16. MR1038232
[10] Grenander, U. and Rosenblatt, M. (1953). Statistical spectral analysis of time series arising from stationary stochastic processes. Ann. Math. Statist. 24 537–558.
MR0058901
[11] Häggström, O. and Rosenthal, J.S. (2007). On variance conditions for Markov chain
CLTs. Electron. Commun. Probab. 12 454–464 (electronic). MR2365647
[12] Hannan, E.J. (1970). Multiple Time Series. New York: Wiley. MR0279952
[13] Jansson, M. (2004). The error in rejection probability of simple autocorrelation robust
tests. Econometrica 72 937–946. MR2051441
[14] Jarner, S.F. and Roberts, G.O. (2002). Polynomial convergence rates of Markov chains.
Ann. Appl. Probab. 12 224–247. MR1890063
[15] Jones, G.L., Haran, M., Caffo, B.S. and Neath, R. (2006). Fixed-width output analysis
for Markov chain Monte Carlo. J. Amer. Statist. Assoc. 101 1537–1547. MR2279478
[16] Kadota, T.T. (1967). Term-by-term differentiability of Mercer’s expansion. Proc. Amer.
Math. Soc. 18 69–72. MR0203397
[17] Kiefer, N.M. and Vogelsang, T.J. (2005). A new asymptotic theory for
heteroskedasticity-autocorrelation robust tests. Econometric Theory 21 1130–1164.
MR2200988
[18] Kiefer, N.M., Vogelsang, T.J. and Bunzel, H. (2000). Simple robust testing of regression hypotheses. Econometrica 68 695–714. MR1769382

32

Y.F. Atchadé

[19] Kipnis, C. and Varadhan, S.R.S. (1986). Central limit theorem for additive functionals
of reversible Markov processes and applications to simple exclusions. Comm. Math.
Phys. 104 1–19. MR0834478
[20] Lindvall, T. (1992). Lectures on the Coupling Method. Wiley Series in Probability and
Mathematical Statistics: Probability and Mathematical Statistics. New York: Wiley.
MR1180522
[21] Meyn, S.P. and Tweedie, R.L. (1993). Markov Chains and Stochastic Stability. Communications and Control Engineering Series. London: Springer. MR1287609
[22] Neave, H.R. (1970). An improved formula for the asymptotic variance of spectrum estimates. Ann. Math. Statist. 41 70–77. MR0251876
[23] Newey, W.K. and West, K.D. (1994). Automatic lag selection in covariance matrix
estimation. Rev. Econ. Stud. 61 631–653. MR1299308
[24] Parzen, E. (1957). On consistent estimates of the spectrum of a stationary time series.
Ann. Math. Statist. 28 329–348. MR0088833
[25] Ponomarëv, S.P. (1987). Submersions and pre-images of sets of measure zero. Sibirsk.
Mat. Zh. 28 199–210. MR0886871
[26] Robert, C.P. and Casella, G. (2004). Monte Carlo Statistical Methods, 2nd ed. Springer
Texts in Statistics. New York: Springer. MR2080278
[27] Roberts, G.O. and Rosenthal, J.S. (2004). General state space Markov chains and
MCMC algorithms. Probab. Surv. 1 20–71. MR2095565
[28] Roberts, G.O. and Tweedie, R.L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli 2 341–363. MR1440273
[29] Schruben, L. (1983). Confidence interval estimation using standardized time series. Oper.
Res. 31 1090–1108.
[30] Sun, Y., Phillips, P.C.B. and Jin, S. (2008). Optimal bandwidth selection
in heteroskedasticity-autocorrelation robust testing. Econometrica 76 175–194.
MR2374985
[31] Tierney, L. (1994). Markov chains for exploring posterior distributions. Ann. Statist. 22
1701–1762. MR1329166
Received February 2014 and revised October 2014

