arXiv:1509.05720v3 [] 19 Oct 2016

SUBSAMPLING FOR GENERAL STATISTICS UNDER LONG
RANGE DEPENDENCE WITH APPLICATION TO CHANGE
POINT ANALYSIS
ANNIKA BETKEN AND MARTIN WENDLER

Abstract. In the statistical inference for long range dependent time series
the shape of the limit distribution typically depends on unknown parameters.
Therefore, we propose to use subsampling. We show the validity of subsampling for general statistics and long range dependent subordinated Gaussian
processes which satisfy mild regularity conditions. We apply our method to
a self-normalized change-point test statistic so that we can test for structural
breaks in long range dependent time series without having to estimate any nuisance parameter. The finite sample properties are investigated in a simulation
study. We analyze three data sets and compare our results to the conclusions
of other authors.

1. Introduction
1.1. Long Range Dependence. While most statistical research is done for independent data or short memory time series, in many applications there are also time
series with long memory in the sense of slowly decaying correlations: in hydrology
(starting with the work of Hurst [31]), in finance (e.g. Lo [39]), in the analysis of
network traffic (e.g. Leland, Taqqu, Willinger and Wilson [37]) and in many other
fields of research.
As model of dependent time series we will consider subordinated Gaussian processes: Let (ξn )n∈N be a stationary sequence of centered Gaussian variables with
Var(ξn ) = 1 and covariance function γ satisfying
(1)

γ(k) := Cov(ξ1 , ξk+1 ) = k −D Lγ (k)

for D > 0 and a slowly varying function Lγ . If D < 1, the spectral density f of
(ξn )n∈N is not continuous, but has a pole at 0. The spectral density has the form
f (x) = |x|D−1 Lf (x)
for a function Lf which is slowly varying at the origin (see Proposition 1.1.14 in
Pipiras and Taqqu [44]).
Furthermore, let G : R → R be a measurable function such that E[G2 (ξ1 )] < ∞.
The stochastic process (Xn )n∈N given by
Xn := G(ξn )
2010 Mathematics Subject Classification. 60G15; 62G09; 60G22.
Key words and phrases. Subsampling; Gaussian Processes; Long Range Dependence; ChangePoint Test.
supported by the German National Academic Foundation, Collaborative Research Center SFB
823 Statistical modelling of nonlinear dynamic processes.
1

2

ANNIKA BETKEN AND MARTIN WENDLER

P∞
is called longPrange dependent if
n=0 |Cov(X1 , Xn+1 )| = ∞, and short range
∞
dependent if n=0 |Cov(X1 , Xn+1 )| < ∞.
Pn
In limit theorems for the partial sum Sn = i=1 Xi , the normalization and the
shape of the limit distribution not only depend on the decay of the covariances γ(k)
as k → ∞, but also on the function G. More precisely, Taqqu [53] and Dobrushin
and Major [25] independently proved that
n
X
1
(Xi − E[Xi ]) ⇒ C(r, H)gr Zr,H (1)
Lγ (n)r/2 nH i=1
1
1
if the Hurst parameter H := max{1 − rD
2 , 2 } is greater than 2 . Here, r denotes
the Hermite rank of the function G, C(r, H) is a constant, gr is the first non-zero
coefficient in the expansion of G as a sum of Hermite polynomials and Zr,H is a
Hermite process. For more details on Hermite polynomials and limit theorems for
subordinated Gaussian processes we recommend the book of Pipiras and Taqqu [44].
In this case (rD < 1), the process (Xn )n∈N is long range dependent as the covariances are not summable. Note that the limiting random variable C(r, H)Zr,H (1) is
Gaussian only if the Hermite rank r = 1.
If rD = 1, the process (Xn )n∈N might be short or long range dependent according
to the slowly varying function P
Lγ . If rD > 1, the process is short range dependent.
n
In this case, the partial sum i=1 (Xi − E[Xi ]) has (with proper normalization)
always a Gaussian limit.
There are other models for long memory processes: Fractionally integrated autoregressive moving average processes can show long range dependence, see Granger
and Joyeux [28]. General linear processes with slowly decaying coefficients were
studied by Surgailis [52].

1.2. Subsampling. For practical applications the parameters D, r and the slowly
varying function Lγ are unknown and thus the scaling needed in the limit theorems
and the shape of the asymptotic distribution are not known, either. That makes it
difficult to use the asymptotic distribution for statistical inference. The situation
gets even more complicated if one is not interested in partial sums, but in nonlinear
statistical functionals. For example, U -statistics can have a limit distribution which
is a linear combination of random variables related to different Hermite ranks, see
Beutner and Zähle [13]. Self-normalized statistics typically converge to quotients
of two random variables (e.g. McElroy and Politis [42]). The change-point test
proposed by Berkes, Horváth, Kokoszka and Shao [11] converges to the supremum
of a fractional Brownian bridge under the alternative hypothesis.
To overcome the problem of the unknown shape of the limit distribution and to
avoid the estimation of nuisance parameters, one would like to use nonparametric
methods. However, Lahiri [36] has shown that the popular moving block bootstrap might fail under long range dependence. Another nonparametric approach
is subsampling (also called sampling window method), first studied by Politis and
Romano [45], Hall and Jing [30], and Sherman and Carlstein [50]. The idea is the
following: Let Tn = Tn (X1 , . . . , Xn ) be a series of statistics converging in distribution to a random variable T . However, as we typically just have one sample,
we observe only one realization of Tn and therefore cannot estimate the distribution of Tn . If l = ln is a sequence with ln → ∞ and ln = o(n), then Tl also

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

3

converges in distribution to T and we have multiple (though dependent) realizations Tl (X1 , . . . , Xl ), Tl (X2 , . . . , Xl+1 ),. . ., Tl (Xn−l+1 , . . . , Xn ), which can be used
to calculate the empirical distribution function.
Note that we do not need to know the limit distribution. In our example (selfnormalized change point test statistic, see Section 3), the shape of the distribution
depends on two unknown parameters, but we can still apply subsampling. However,
for other statistics, one needs an unknown scaling to achieve convergence. If this is
the case, one has to estimate the scaling parameters before applying subsampling.
UnderPlong range dependence the validity of subsampling for the sample mean
n
X̄ = n1 i=1 Xi has been investigated in the literature starting with Hall, Jing
and Lahiri [29] for subordinated Gaussian processes. Nordman and Lahiri [43] and
Zhang, Ho, Wendler and Wu [58] studied linear processes with slowly decaying
coefficients. For the case of Gaussian processes an alternative proof can be found
in the book of Beran, Feng, Ghosh and Kulik [10].
It was noted by Fan [26] that the proof in [29] can be easily generalized to other
statistics than the sample mean. However, the assumptions on the Gaussian process
are restrictive (see also [42]). Their conditions imply that the sequence (ξn )n∈N is
completely regular, which might hold for some special cases (see Ibragimov and
Rozanov [32]), but excludes many examples:
Example 1 (Fractional Gaussian Noise). Let (BH (t))t∈[0,∞) be a fractional Brownian motion, i.e. a centered, self-similar Gaussian process with covariance function

1
|t|2H + |s|2H − |t − s|2H
E [BH (t)BH (s)] =
2
for some H ∈ ( 12 , 1). Then, (ξn )n∈N given by ξn = BH (n) − BH (n − 1) is called
fractional Gaussian noise. By self-similarity we have
X

n
3n
X
ξi ,
ξj = corr (BH (n), BH (3n) − BH (2n))
corr
i=1

j=2n+1

= corr (BH (1), BH (3) − BH (2)) .
As a result, the correlations of linear combinations of observations in the past and
future do not vanish if the gap between past and future grows. Thus, fractional
Gaussian noise is not completely regular.
Jach, McElroy and Politis [33] provided a more general result on the validity
of subsampling. They assume that the function G has Hermite rank 1, that G is
invertible and Lipschitz-continuous and that the process (ξn )n∈N has a causal representation as a functional of an independent sequence of random variables. These
assumptions are difficult to check in practice. Moreover, although not explicitly
stated in [33], the statistic Tn has to be Lipschitz-continuous (uniformly in n),
which is not satisfied by many robust estimators (see Section 3 for an example).
The main aim of this paper is to establish the validity of the subsampling method
for general statistics Tn without any assumptions on the continuity of the statistic,
on the function G and only mild assumptions on the Gaussian process (ξn )n∈N .
Independently of our research, similar theorems have been proved by Bai, Taqqu
and Zhang [6]. We will discuss their results after our main theorem in Section 2.
In Section 3 we will apply our theorem to a self-normalized, robust change-point
statistic. The finite sample properties of this test will be investigated in a simulation

4

ANNIKA BETKEN AND MARTIN WENDLER

study in Section 4. Finally, the proof of the main result and the lemmas needed
can be found in Section 5.
2. Main Results
2.1. Statement of the Theorem. For a statistic Tn = Tn (X1 , . . . , Xn ) the subsampling estimator F̂l,n of the distribution function FTn with FTn (t) = P (Tn ≤ t)
is defined in the following way: For t ∈ R let
F̂l,n (t) =

n−l+1
X
1
1{Tl (Xi ,...,Xi+l−1 )≤t} .
n − l + 1 i=1

Our first assumption guarantees the convergence of the distribution function FTn :
Assumption 1. (Xn )n∈N is a stochastic process and (Tn )n∈N is a sequence of
statistics such that Tn ⇒ T in distribution as n → ∞ for a random variable T with
distribution function FT .
This is a standard assumption for subsampling, see for example [45]. If the
distribution does not converge, we cannot expect the distribution of Tl to be close
to the distribution of Tn .
Next, we will formulate our conditions on the sequence of random variables
(Xn )n∈N :
Assumption 2. Xn = G(ξn ) for a measurable function G and a stationary, Gaussian process (ξn )n∈N with covariance function
γ(k) := Cov(ξ1 , ξ1+k ) = k −D Lγ (k)
such that the following conditions hold:
(1) D ∈ (0, 1] and Lγ is a slowly varying function with
max
k̃∈{k+1,...,k+2l0 −1}

Lγ (k) − Lγ (k̃) ≤ K

l0
min {Lγ (k), 1}
k

for a constant K < ∞ and all l0 ∈ {lk , . . . , k}.
(2) (ξn )n∈N has a spectral density f with f (x) = |x|D−1 Lf (x) for a slowly
varying function Lf which is bounded away from 0 on [0, π] such that
limx→0 Lf (x) ∈ (0, ∞] exists.
While we have some regularity conditions on the underlying Gaussian process
(ξn )n∈N , we do not impose any conditions on the function G: no finite moments
or continuity are required, so that our results are applicable for heavy-tailed random variables and robust test statistics. In the next subsection we will show that
Assumption 2 holds for some standard examples of long range dependent Gaussian
processes.
Furthermore, we need a restriction on the growth rate of the block length l:
Assumption 3. Let (ln )n∈N be a non-decreasing
sequence of integers such that

l = ln → ∞ as n → ∞ and ln = O n(1+D)/2− for some  > 0.
If the dependence of the underlying process (ξn )n∈N gets stronger, the range√of
possible values for l gets smaller. A popular choice for the block length is l ≈ C n
(see for example [29]), which is allowed for all D ∈ (0, 1]. Now, we can state our
main result:

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

5

Theorem 1. Under Assumptions 1, 2 and 3 we have
P

FTn (t) − F̂l,n (t) −→ 0
as n → ∞ for all points of continuity t of FT . If FT is continuous, then
P

sup FTn (t) − F̂l,n (t) −→ 0.
t∈R

As a result, we have a consistent estimator for the distribution function of Tn .
It is possible to build tests and confidence intervals based on this estimator.
If D > 1, the process (ξn )n∈N is strongly mixing due to Theorem 9.8 in the book
of Bradley [18]. The statements of Theorem 1 hold by Corollary 3.2 in [45] for any
block length l satisfying l → ∞ and l = o(n).
In a recent article, Bai et al. [6] have shown that subsampling is consistent
for long range dependent Gaussian processes without any extra assumptions on
the slowly varying function Lf , but with a stronger restriction on the block size l,
namely l = o(n2−2H Lγ (n)). In another article by Bai and Taqqu [5], the validity
of subsampling is shown under the mildest possible assumption on the block length
(l = o(n)). The condition on the spectral density is slightly stronger than our
condition, the case limx→0 Lf (x) = ∞ is not allowed.
2.2. Examples for our Assumptions. We will now give two examples of Gaussian processes satisfying Assumption 2:
Example 2 (Fractional Gaussian Noise). Fractional Gaussian Noise (ξn )n∈N with
Hurst parameter H as introduced in Example 1 has the covariance function


1
γ(k) =
|k − 1|2H − 2|k|2H + |k + 1|2H = H(2H − 1) k −D + h(k)k −D−1
2
for D = 2 − 2H and a function h bounded by a constant M < ∞. This can be
easily seen by means of a Taylor expansion. Hence, Lγ (k) = H(2H − 1)(1 + h(k)/k)
and for all k̃ ≥ k
M
h(k) h(k̃)
1
≤ H(2H − 1)
−
=: K .
k
k
k
k̃
This implies part 1 of Assumption 2. For the second part note that the spectral
density f corresponding to fractional Gaussian noise is given by
Lγ (k) − Lγ (k̃) ≤ H(2H − 1)

f (λ) = C(H)(1 − cos(λ))

∞
X

|λ + 2kπ|

D−3

k=−∞

P∞
D−3
1 − cos(λ) k=−∞ |λ + 2kπ|
=λ
C(H)
,
λ2
λD−3
see Sinai [51]. The slowly varying function
P∞
D−3
1 − cos(λ) k=−∞ |λ + 2kπ|
Lf (λ) = C(H)
λ2
λD−3
is bounded away from 0 because this holds for the first factor (1 − cos(λ))/λ2 and
since
P∞
D−3
D−3
|λ + 0π|
k=−∞ |λ + 2kπ|
≥
= 1.
λD−3
λD−3
D−1

6

ANNIKA BETKEN AND MARTIN WENDLER

Example 3 (Gaussian FARIMA processes). Let (εn )n∈Z be Gaussian white noise
with variance σ 2 = Var(ε0 ). Then, for d ∈ (0, 1/2), a FARIMA(0, d, 0) process
(ξn )n∈N is given by
∞
X
Γ(j + d)
ξn =
εn−j .
Γ(j + 1)Γ(d)
j=0
According to Pipiras and Taqqu [44], Section 1.3, it has the specral density

1−D
2
σ2
|λ|
−iλ −2d
D−1 σ
f (λ) =
|1 − e |
= |λ|
2π
2π |1 − e−iλ |
with D = 1 − 2d ∈ (0, 1). As |1 − e−iλ | ≤ λ, part 2 of Assumption 2 holds. For
part 1 we have by Corollary 1.3.4 of [44] that
Γ(k + d)
Γ(1 − 2d)
.
Γ(1 − d)Γ(d) Γ(k − d + 1)
1/2 x x

Recall that by the Stirling formula Γ(x) = 2π
1+O(x−1 ) . Consequently,
x
e
k−d+1 
Γ(1 − 2d) −2d+1 2d−1  k + d k+d 
k
1 
γ(k) = σ 2
e
k
1+O
.
Γ(1 − d)Γ(d)
k
k−d+1
k

Using a Taylor expansion of (k + d) log(k + d) − log(k) + (k − d + 1) log(k) −
log(k − d + 1) , it easily follows that
γ(k) = σ 2

γ(k) = k −D Lγ (k)
with Lγ (k) = C + O(1/k) for some constant C. Part 1 of Assumption 2 follows in
the same way as in Example 2.
It would be interesting to know, if the sampling window method is also consistent for long range dependent linear processes and general statistics without the
assumption of Gaussianity. However, this seems to be a very difficult problem and
is beyond the scope of this article.
3. Applications
3.1. Robust, Self-Normalized Change-Point Test. In this paper, the main
motivation for considering subsampling procedures in order to approximate the
distribution of test statistics consists in avoiding the choice of unknown parameters.
As an example we will consider a self-normalized test statistic that can be applied
to detect changes in the mean of long range dependent and heavy-tailed time series.
Given observations X1 , . . . , Xn with Xi = µi + G(ξi ) we are concerned with a
decision on the change-point problem
H : µ1 = . . . = µn
against
A : µ1 = . . . = µk 6= µk+1 = . . . = µn for some k ∈ {1, . . . , n − 1} .
Under the hypothesis H we assume that the data generating process (Xn )n∈N is stationary, while under the alternative A there is a change in location at an unknown
point in time. This problem has been widely studied: Csörgő and Horváth [21]
give an overview of parametric and non-parametric methods that can be applied in
order to detect change-points in independent data.

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

7

Many commonly used testing procedures are based on Cusum (cumulative sum)
test statistics, but when applied to data sets generated by long range dependent
processes, these change-point tests often falsely reject the hypothesis of no change
in the mean (see also Baek and Pipiras [4]). Furthermore, the performance of
Cusum-like change-point tests is sensitive to outliers in the data.
In contrast, testing procedures that are based on rank statistics have the advantage of not being sensitive to outliers in the data. Rank-based tests were introduced
by Antoch, Hušková, Janic and Ledwina [3] for detecting changes in the distribution function of independent random variables. Wilcoxon-type rank tests have
been studied by Wang [56] in the presence of linear long memory time series and
by Dehling, Rooch and Taqqu [22] for subordinated Gaussian sequences.
Note that the normalization of the Wilcoxon change-point test statistic as proposed in [22] depends on the slowly varying function Lγ , the LRD parameter D
and the Hermite rank r of the class of functions 1{Xi ≤x} − F (x), x ∈ R. Although
many authors assume r = 1 and while there are well-tried methods to estimate D,
estimating Lγ does not seem to be an easy task. For this reason, the Wilcoxon
change-point test does not seem to be suitable for applications to real data.
To avoid these issues, Betken [12] proposes an alternative normalization for the
Wilcoxon change-point test. This normalization approach has originally been established by Lobato [40] for decision on the hypothesis that a short range dependent
stochastic process is uncorrelated up to a lag of a certain order. In change-point
analysis, the normalization has recently been applied to several test statistics: Shao
and Zhang [49] define a self-normalized Kolmogorov-Smirnov test statistic that
serves to identify changes in the mean of short range dependent time series. Shao
[48] adopted the normalization so as to define an alternative normalization for a
Cusum test which detects changes in the mean of short range dependent as well as
long range dependent time series.
For the definition of the self-normalized
Wilcoxon test statistic, we introduce
Pn
the ranks Ri := rank(Xi ) = j=1 1{Xj ≤Xi } for i = 1, . . . , n. It seems natural to
transfer the normalization that has been used in [48] to the Cusum test statistic
of the ranks in order to establish a self-normalized version of the Wilcoxon test
statistic, which is robust to outliers in the data. Therefore, the corresponding
two-sample test statistic is defined by
Pk
Gn (k) := 

i=1

1
n

Pk

t=1

Ri −

St2 (1, k) +

1
n

k
n

Pn

Pn

i=1

t=k+1

Ri

St2 (k + 1, n)

1/2 ,

where
St (j, k) :=

t
X
h=j

k

Rh − R̄j,k



with R̄j,k :=

X
1
Rt .
k − j + 1 t=j

The self-normalized Wilcoxon change-point test rejects the hypothesis for large
values of maxk∈{bnτ1 c,...,bnτ2 c} |Gn (k)|, where 0 < τ1 < τ2 < 1. The proportion of
the data that is included in the calculation of the supremum is restricted by τ1 and
τ2 . A common choice for these parameters is τ1 = 1 − τ2 = 0.15; see Andrews [2].
For long range dependent subordinated Gaussian processes (Xn )n∈N , the asymptotic distribution of the test statistic under the hypothesis H can be derived by the

8

ANNIKA BETKEN AND MARTIN WENDLER

continuous mapping theorem (see Theorem 1 in [12]):

Tn (τ1 , τ2 ) :=
⇒

sup
τ1 ≤λ≤τ2

max

k∈{bnτ1 c,...,bnτ2 c}

Rλ
0

(Zr (t) −

|Gn (k)|
|Zr (λ) − λZr (1)|
R 1−λ
+ 0 (Zr? (t) −

t
2
λ Zr (λ)) dt

t
?
1−λ Zr (1

− λ))2 dt

1/2

.

Here, Zr is an r-th order Hermite process with Hurst parameter H := max{1 −
rD 1
?
2 , 2 } and Zt (r) = Zr (1) − Zr (1 − t). A comparison of Tn (τ1 , τ2 ) with the critical
values of its limit distribution still presupposes determination of these parameters.
We can bypass the estimation of D and r by applying the subsampling procedure
since, due to the convergence of Tn (τ1 , τ2 ), Assumption 1 holds.
Note that even under the alternative A (change in location), we have to find
the quantiles of the distribution under the hypothesis (stationarity). As the block
length l is much shorter than the sample size n, most blocks will not be contaminated by the change-point so that the distribution of the test statistic will not
change that much. The accuracy and the power of the test will be investigated by
a simulation study in Section 4.
If the distribution of Xi is notPcontinuous, there might be ties in the data and
n
consideration of the ranks Ri = j=1 1{Xj ≤Xi } may not be appropriate. We proPn
pose to use a modified statistic based on the modified ranks R̃i = j=1 (1{Xj <Xi } +
1
2 1{Xj =Xi } ) in this case. For the convergence of the corresponding self-normalized
change point test see Appendix A.
The test statistic Tn (τ1 , τ2 ) is designed for the detection of a single change-point.
An extension of the testing procedure that allows for multiple change-points is possible by adapting Shao’s testing procedure which takes this problem into consideration (see [48]). For convenience, we describe the construction of the modified test
statistic in the case of two change-points. The general idea consists in dividing the
sample given by X1 , . . . , Xn according to the pair (k1 , k2 ) of potential change-point
locations and to compute the original test statistic with respect to the subsamples
X1 , . . . , Xk2 and Xk1 +1 , . . . , Xn . We reject the hypothesis for large values of the
sum of the corresponding single statistics.
For ε ∈ (0, τ2 − τ1 ) define Tn (τ1 , τ2 , ε) := sup(k1 ,k2 )∈Ωn (τ1 ,τ2 ,ε) |Gn (k1 , k2 )|, where
Ωn (τ1 , τ2 , ε) := {(k1 , k2 ) : bnτ1 c ≤ k1 < k2 ≤ bnτ2 c, k2 − k1 ≥ bnεc} and

(1)

Pk1
Gn (k1 , k2 ) := 

Ri

i=1

1
n

Pk1 
t=1

2

(1)
St (1, k1 )

Pk2
+

+

1
n

Pk2

t=k1 +1



(2)
St (k1

k1
k2

i=1



t=k1 +1

Ri

(1)

Pk2

Pk2
(2)

i=k1 +1

1
n

−

−

Ri

(1)
St (k1

+ 1, k2 )

Pn

Ri

k2 −k1
n−k1

2
+ 1, k2 ) +

1
n

i=k1 +1

2 1/2

(2)


2 1/2
(2)
t=k2 +1 St (k1 + 1, n)

Pn

,

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

9

where
(1)

Ri

:=

k2
X

1{Xj ≤Xi } ,

(2)

Ri . =

j=1
(h)

St (j, k) :=

n
X

1{Xj ≤Xi } ,

j=k1 +1
t 
X

(h)

Ri

(h)

− R̄j,k



k

(h)

with R̄j,k :=

i=j

X (h)
1
R .
k − j + 1 t=j t

The distribution of the test statistic converges to a limit T (r, τ1 , τ2 , ε) (see Appendix B), so subsampling can be applied. The critical values corresponding to the
asymptotic distribution of the test statistic are reported in Table 1.
Table 1. Simulated critical values for the distribution of
T (1, τ1 , τ2 , ε) when [τ1 , τ2 ] = [0.15, 0.85] and ε = 0.15. The sample
size is 1000, the number of replications is 10, 000.

H
H
H
H
H
H

= 0.501
= 0.6
= 0.7
= 0.8
= 0.9
= 0.999

10%
5%
17.79 19.76
19.80 22.38
22.08 24.95
24.24 27.61
26.50 30.11
28.28 32.32

1%
24.13
27.68
30.46
34.04
37.78
41.24

3.2. Data Examples. We will revisit some data sets which have been analyzed
before in the literature. We will use the self-normalized Wilcoxon change-point test
combined with subsampling and compare our findings to the conclusions of other
authors.
The plot in Figure 1 depicts the annual volume of discharge from the Nile river
at Aswan in 108 m3 for the years 1871 to 1970. The data set has been analyzed for
the detection of a change-point by numerous authors under differing assumptions
concerning the data generating random process and by usage of diverse methods.
Amongst others, Cobb [19], MacNeill, Tang and Jandhyala [41], Wu and Zhao [57]
and Shao [48] provided statistically significant evidence for a decrease of the Nile’s
annual discharge towards the end of the 19th century. The construction of the
Aswan Low Dam between 1898 and 1902 serves as a popular explanation for an
abrupt change in the data.
The value of the self-normalized Wilcoxon test statistic computed with respect
to the data is given by Tn (τ1 , τ2 ) = 13.48729. For a level of significance of 5%, the
self-normalized Wilcoxon
change-point test rejects the hypothesis for every possible

value of H ∈ 12 , 1 . Furthermore, we approximate the distribution of the selfnormalized
Wilcoxon test statistic by the sampling window method with block size
√
l = b nc = 10. The subsampling-based test decision also indicates the existence
of a change-point in the mean of the data, even if we consider the 99%-quantile of
F̂l,n .
In particular, previous analysis of the Nile data done by Wu and Zhao [57] and
Balke [7] suggests that the change in the discharge volume occurred in 1899. We
applied the self-normalized Wilcoxon test and the sampling window method to the
corresponding pre-break and post-break samples. Neither of these two approaches

600

800

1000 1200 1400

ANNIKA BETKEN AND MARTIN WENDLER

annual discharge

10

1880

1900

1920

1940

1960

time
Figure 1. Measurements of the annual discharge of the river Nile
at Aswan in 108 m3 for the years 1871-1970. The dotted line indicates the location of the change-point; the dashed lines designate
the sample means for the pre-break and post-break samples.

leads to rejection of the hypothesis, so that it seems reasonable to consider both
samples as stationary. At this point, it is interesting to note that, based on the whole
sample, local Whittle estimation with bandwidth parameter m = bn2/3 c suggests
the existence of long range dependence characterized by an Hurst parameter Ĥ =
0.962, whereas the estimates for the pre-break and post-break samples given by
Ĥ1 = 0.517 and Ĥ2 = 0.5, respectively, should be considered as indication of
short range dependent data. In this regard, our findings support the conjecture
of spurious long memory caused by a change-point and therefore coincide with the
results of Shao [48].
The second data set consists of the seasonally adjusted monthly deviations of
the temperature (degrees C) for the northern hemisphere during the years 1854 to
1989 from the monthly averages over the period 1950 to 1979. The data results
from spatial averaging of temperatures measured over land and sea. At first sight,
the plot in Figure 2 may suggest an increasing trend as well as an abrupt change of
the temperature deviations. Statistical evidence for a positive deterministic trend
implies affirmation of the conjecture that there has been global warming during the
last decades.
In scientific discourse, the question of whether the Northern hemisphere temperature data acts as an indicator for global warming of the atmosphere is a controversial
issue. Deo and Hurvich [24] provided some indication for global warming by fitting

0.0

0.5

11

−1.5 −1.0 −0.5

Northern hemisphere temperature

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

1860

1900

1940

1980

time
Figure 2. Monthly temperature of the northern hemisphere for
the years 1854-1989 from the data base held at the Climate Research Unit of the University of East Anglia, Norwich, England.
The temperature anomalies (in degrees C) are calculated with respect to the reference period 1950-1979. The dotted line indicates
the location of the potential change-point; the dashed lines designate the sample means for the pre-break and post-break samples.
a linear trend to the data. Beran and Feng [8] considered a more general stochastic model by the assumption of so-called semiparametric fractional autoregressive
(SEMIFAR) processes. Their method did not deliver sufficient statistical evidence
for a deterministic trend. Wang [55] applied another method for the detection of
gradual change to the global temperature data and did not detect an increasing
trend , either. Nonetheless, he offers an alternative explanation for the occurrence
of a trend-like behavior by pointing out that it may have been generated by stationary long range dependent processes. In contrast, it is shown in Shao [48] that
the existence of a change-point in the mean yields yet another explanation for the
performance of the data.
The value of the self-normalized Wilcoxon test statistic computed with respect
to the data is given by Tn (τ1 , τ2 ) = 18.98636. Consequently, the self-normalized
Wilcoxon change-point test would reject the hypothesis for every possible value of
H ∈ 12 , 1 at a level of significance of 1%. In addition, an application of the sampling window method with respect to the self-normalized Wilcoxon test statistic
based on comparison of Tn (τ1 , τ2 ) with the 99%-quantile of the sampling distribution F̂l,n yields a test decision in favor of the alternative hypothesis for any choice
of the block length l ∈ {bnγ c| γ = 0.3, 0.4, . . . , 0.9} = {9, 19, 40, 84, 177, 371, 778}.

12

ANNIKA BETKEN AND MARTIN WENDLER

All in all, both testing procedures provide strong evidence for the existence of a
change in the mean.
According to Shao [48] the change-point is located around October 1924. Based
on the whole sample local Whittle estimation with bandwidth m = bn2/3 c provides
an estimator Ĥ = 0.811. The estimated Hurst parameters for the pre-break and
post-break sample are Ĥ1 = 0.597 and Ĥ2 = 0.88, respectively. Neither of both
testing procedures, i.e. subsampling with respect to the self-normalized Wilcoxon
test statistic and comparison of the value of Tn (τ1 , τ2 ) with the corresponding critical values of its limit distribution, provides evidence for another change-point in
the pre-break or post-break sample.
Moreover, computation of the test statistic that allows for two change-point
locations yields Tn (τ1 , τ2 , ε) = 17.88404 (for τ1 = 1 − τ2 = ε = 0.15), i.e. if
compared to the values in Table 1, the test statistic only surpasses the critical value
corresponding to H = 0.501 and a significance level of 10%, but does not exceed any
of the other values. Subsampling with respect to the test statistic Tn (τ1 , τ2 , ε) does
not support the conjecture of two changes, either. In fact, subsampling leads to a
rejection of the hypothesis when the block length equals l = bn0.7 c = 177 (based on
a comparison of Tn (τ1 , τ2 , ε) with the 95%-quantile of the corresponding sampling
distribution F̂l,n ), but yields a test decision in favor of the hypothesis for block
lengths l ∈ {bnγ c| γ = 0.5, 0.6, 0.8, 0.9} = {40, 84, 371, 778} and for comparison
with the 90%-quantile of F̂l,n .
Therefore, it seems safe to conclude that the appearance of long memory in
the post-break sample is not caused by another change-point in the mean. The
pronounced difference between the local Whittle estimators Ĥ1 and Ĥ2 suggests a
change in the dependence structure of the times series. Another explanation might
be a gradual change of the temperature in the post-break period. We conjecture
that our test has only low power in the case of a gradual change, because the denominator of our self-normalized test statistic is inflated as the ranks systematically
deviate from the mean rank of the first and second part. When using subsampling,
the trend also appears in subsamples so that we fail to approximate the distribution
under the hypothesis.
As pointed out by one of the referees, the Northern hemisphere temperature data
does not seem to be second-order stationary; the variance in the first part of the
time series seems to be higher. Such a change in variance should also result in a
loss of power. The reason is that the ranks in the part with the higher variance are
more extreme, so that the distance to the mean rank of this part is larger. This
leads to a higher value of the denominator of our self-normalized test statistic and
consequently to a lower value of the ratio.
The third data set consists of the arrival rate of Ethernet data (bytes per 10
milliseconds) from a local area network (LAN) measured at Bellcore Research and
Engineering Center in 1989. For more information on the LAN traffic monitoring we
refer to Leland and Wilson [38] and Beran [9]. Figure 3 reveals that the observations
are strongly right-skewed. As the self-normalized Wilcoxon test is based on ranks,
we do not expect that this will affect our analysis.
Coulon, Chabert and Swami [20] examined this data set for change-points before.
The method proposed in their paper is based on the assumption that a FARIMA
model holds for segments of the data. The number of different sections and the

13

8000
6000
4000
0

2000

bytes per 0.01 seconds

10000

12000

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

0

1000

2000

3000

4000

time

Figure 3. Ethernet traffic in bytes per 10 milliseconds from a
LAN measured at Bellcore Research Engineering Center.
location of the change-points are chosen by a model selection criterion. The algorithm proposed by Coulon et al. [20] detects multiple changes in the parameters of
the corresponding FARIMA time series.
In contrast, an application of the self-normalized Wilcoxon change-point test
does not provide evidence for a change-point in the mean: the value of the test
statistic is given by Tn (τ1 , τ2 ) = 3.270726, i.e. even for a level of significance
of 10%, the self-normalized Wilcoxon
change-point test does not reject the hy
pothesis for any value H ∈ 12 , 1 . Furthermore, subsampling with respect to the
self-normalized Wilcoxon test statistic does not lead to a rejection of the hypothesis , either (for any choice of block length l ∈ {bnγ c| γ = 0.3, 0.4, . . . , 0.9} =
{12, 27, 63, 144, 332, 761, 1745} and for comparison with the 90%-quantile of the
corresponding sampling distribution F̂l,n ).
Taking into consideration that the data set contains ties (the value 0 appears
several times), we also applied the self-normalized Wilcoxon test statistic based on
the modified ranks R̃i and used subsampling with respect to this statistic. Both
approaches did not lead to a rejection of the hypothesis.
An application of the test statistic constructed for the detection of two changes
yields a value of Tn (τ1 , τ2 , ε) = 15.24527 when ε = τ1 = 1 − τ2 = 0.15. Clearly, this
does not lead to a rejection of the hypothesis for any value of the parameter H. In

14

ANNIKA BETKEN AND MARTIN WENDLER

addition, subsampling based on comparison of Tn (τ1 , τ2 , ε) with the 90%-quantile
of the corresponding sampling distribution F̂l,n does not provide evidence for the
assertion of multiple changes for any block lenght l ∈ {bnγ c| γ = 0.5, 0.6, 0.7, 0.8} =
{63, 144, 332, 761} in the data, either.
These results do not coincide with the analysis of the previous authors. On the
one hand this may be due to the fact that the applied methods differ considerably
from the testing procedures applied before. On the other hand, the change-point
estimation algorithm proposed in Coulon, Chabert and Swami [20] is not robust
to skewness or heavy-tailed distributions and decisively relies on the assumption
of FARIMA time series. However, this seems to contradict observations made by
Bhansali and Kokoszka [14] as well as Taqqu and Teverovsky [54] who stress that
the model that fits the Ethernet traffic data is very unlikely to be FARIMA.
Estimation of the Hurst parameter by the local Whittle procedure with bandwidth parameter m = bn2/3 c yields an estimate of Ĥ = 0.845 and therefore indicates the existence of long range dependence. This is consistent with the results of
Leland et al. [37] and Taqqu and Teverovsky [54].
In the three data examples, we find that the results obtained by subsampling
and by parameter estimation are in good accordance with each other. The methods
take into account long range dependence or heavy tails, but still detect a change in
location in the first two examples. For the third data example our analysis supports
the hypothesis of stationarity.
4. Simulations
We will now investigate the finite sample performance of the subsampling procedure with respect to the self-normalized Wilcoxon test and with respect to the
classical Wilcoxon change-point test. Moreover, we will compare these results to the
performance of the tests when the test decision is based on critical values obtained
from the asymptotic distribution of the test statistic.
For this purpose, we consider subordinated Gaussian time series (Xn )n∈N , Xn =
G(ξn ), where (ξn )n∈N is fractional Gaussian noise (introduced in Examples 1 and
2) with Hurst parameter H ∈ {0.6, 0.7, 0.8, 0.9} and covariance function


D
−D
(1 − D) ,
γ(k) ∼ k
1−
2
where D = 2−2H. Initially, we take G(t) = t, so that (Xn )n∈N has normal marginal
distributions. We also consider the transformation

− 21 

1
βk 2
βk
−β
G(t) =
k(Φ(t))
−
(β − 1)2 (β − 2)
β−1
(with Φ denoting the standard normal distribution function) so as to generate
Pareto-distributed data with parameters k, β > 0 (referred to as Pareto(β, k)). In
both cases, the Hermite rank r of 1{G(ξi )≤x} − F (x), x ∈ R, equals r = 1 and
Z
1
J1 (x)dF (x) = √ ;
2 π
R
see [22].
Under the above conditions, the critical values of the asymptotic distribution of
the self-normalized Wilcoxon test statistic are reported in Table 2 in [12]. The limit

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

15

of the Wilcoxon change-point test statistic can be found in [22], the corresponding
critical values can be taken from Table 1 in [12].
The frequencies of rejections of both testing procedures are reported in Table
2 and Table 3 for the self-normalized Wilcoxon change-point test and in Table 4
and Table 5 for the classical Wilcoxon test (without self-normalization). The calculations are based on 5, 000 realizations of time series with sample size n = 300
and n = 500. We have chosen block lengths l = ln = bnγ c with γ ∈ {0.4, 0.5, 0.6}.
As level of significance we chose 5%, i.e. we compare the values of the test statistic with the corresponding critical values of its asymptotic distribution and the
corresponding quantile of the empirical distribution function F̂l,n , respectively.
For the usual testing procedures the estimation
of the Hermite rank r, the slowly
R
varying function Lγ and the integral J1 (x)dF (x) is neglected. Yet, for every
simulated time series we estimate the Hurst parameter H by the local Whittle
estimator Ĥ proposed in Künsch [35]. This estimator is based on an approximation
of the spectral density by the periodogram at the Fourier frequencies. It depends on
the spectral bandwidth parameter m = m(n) which denotes the number of Fourier
1
+ m
frequencies used for the estimation. If the bandwidth m satisfies m
n −→
0 as n −→ ∞, the local Whittle estimator is a consistent estimator for H; see
Robinson [47]. For convenience we always choose m = bn2/3 c in this article. The
critical values corresponding to the estimated values of H are determined by linear
interpolation.
Under the alternative A we analyze the power of the testing procedures (the
frequency of rejection) by considering different choices for the height of the level
shift (denoted by h) and the location [nτ ] of the change-point. In the tables the
columns that are superscribed by “h = 0” correspond to the frequency of a type 1
error, i.e. the rejection rate under the hypothesis H.
For the self-normalized Wilcoxon change-point test (based on the asymptotic
distribution), the empirical size almost equals the level of significance of 5% for
normally distributed data (see Table 2). The sampling window method yields rejection rates that slightly exceed this level. For Pareto(3, 1) time series both testing
procedures lead to similar results and tend to reject the hypothesis too often when
there is no change. With regard to the empirical power, it is notable that for
fractional Gaussian noise time series the sampling window method yields considerably better power than the test based on asymptotic critical values. If Pareto(3,
1)-distributed time series are considered, the empirical power of the subsampling
procedure is still better than the empirical power that results from using asymptotic critical values. However, in this case, the deviation of the rejection rates is
rather small. While the empirical size is not much affected by the Hurst parameter
H, the empirical power is lower for H = 0.8, 0.9.
Considering the classical Wilcoxon test (without self-normalization), it is notable
that for both procedures the empirical size is in most cases not close to the nominal
level of significance (5%), ranging from 1.1% to 20.8% using subsampling and from
2.6% to 36.0% using asymptotic critical values. In general, the sampling window
method becomes more conservative for higher values of the Hurst parameter H,
while the test based on the asymptotic distribution becomes more liberal. Under the
alternative, the usual application of the Wilcoxon test yields better power than the
sampling window method, especially for high values of H. It should be emphasized

16

ANNIKA BETKEN AND MARTIN WENDLER

that this comparison is problematic because the rejection frequencies under the
hypothesis differ.
We conclude that the self-normalized Wilcoxon change-point test is more reliable
than the classical change-point test. The reason might be that in the scaling of the
classical test, the estimator Ĥ of the Hurst parameter enters as a power of the
sample size n. Thus, a small error in this estimation might lead to a large error
in the value of the test statistic. By using the sampling window method for the
self-normalized version, we avoid the estimation of unknown parameters so that the
performance is similar to the performance of the classical testing procedure which
compares the values of the test statistic with the corresponding critical values.
Note that in most cases covered by our simulations the choice of the block length
for the subsampling procedure does not have a big impact on the frequency of a type
1 error. Considering the self-normalized Wilcoxon change-point test, an increase of
the block length tends to go along with a decrease in power, especially for big values
of the Hurst parameter H and Pareto-distributed random variables. For smaller
values of H the effect is not pronounced. We recommend using a block length bn0.4 c
or bn0.5 c for the self-normalized change-point test as the choice l = bn0.6 c implies
worse properties in most cases.
An application of the subsampling testing procedure to the classical (non-selfnormalized) Wilcoxon test for different choices of the block length shows the opposite effect on the rejection rate under the alternative: an increase of the block
length results in a higher frequency of rejections. Here, the block length bn0.6 c
leads to better results in many cases. However, we recommend to not use this test,
but to self-normalize the test statistic instead.
An alternative way of choosing the block length would be to apply the datadriven block selection rule proposed by Götze and Rac̆kauskas [27] and Bickel and
Sakov [15]. Although the algorithm had originally been implemented for applications of the m-out-of-n bootstrap to independent and identically distributed data, it
also lead to satisfactory simulation results in applications to long range dependent
time series (see [33]). Another general approach to the selection of the block size
in the context of hypothesis testing is given by Algorithm 9.4.2 in Politis, Romano
and Wolf [46].

Table 2. Rejection rates of the self-normalized Wilcoxon change-point test obtained by subsampling (left) with block
length l = bnγ c, γ ∈ {0.4, 0.5, 0.6}, and by comparison with asymptotic critical values (right) for fractional Gaussian
noise of length n with Hurst parameter H.
sampling window method

asymptotic distribution

fGn
n
H = 0.6 300

500

H = 0.7 300

500

H = 0.8 300

500

H = 0.9 300

h=0
0.041
0.064
0.070
0.053
0.060
0.069
0.057
0.070
0.077
0.056
0.059
0.065
0.070
0.067
0.073
0.066
0.068
0.069
0.093
0.074
0.073
0.079
0.067
0.063

h = 0.5
0.263
0.313
0.322
0.396
0.421
0.411
0.155
0.171
0.177
0.183
0.193
0.192
0.126
0.117
0.114
0.121
0.114
0.119
0.126
0.097
0.095
0.105
0.091
0.087

h=1
0.700
0.742
0.705
0.859
0.861
0.829
0.412
0.423
0.403
0.513
0.508
0.476
0.251
0.234
0.218
0.295
0.278
0.257
0.208
0.161
0.145
0.194
0.166
0.146

h = 0.5
0.502
0.570
0.555
0.697
0.720
0.697
0.291
0.313
0.314
0.382
0.382
0.387
0.223
0.208
0.201
0.217
0.210
0.205
0.209
0.169
0.165
0.185
0.162
0.152

h=1
0.952
0.964
0.943
0.994
0.995
0.991
0.759
0.763
0.737
0.856
0.854
0.819
0.526
0.494
0.466
0.591
0.567
0.532
0.462
0.397
0.367
0.461
0.416
0.391

τ = 0.25

τ = 0.5

h=0

h = 0.5

h=1

h = 0.5

h=1

0.044

0.209

0.521

0.424

0.861

0.049

0.303

0.687

0.577

0.958

0.053

0.108

0.268

0.228

0.611

0.048

0.133

0.359

0.302

0.730

0.048

0.081

0.144

0.141

0.362

0.053

0.085

0.198

0.163

0.462

0.057

0.065

0.106

0.125

0.308

0.051

0.068

0.120

0.128

0.350

17

500

l
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41

τ = 0.5

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

τ = 0.25

sampling window method

asymptotic distribution

τ = 0.25
n
300

500

H = 0.7

300

500

H = 0.8

300

500

H = 0.9

300

500

l
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41

h=0
0.041
0.067
0.070
0.055
0.066
0.071
0.057
0.064
0.077
0.066
0.068
0.072
0.070
0.068
0.073
0.064
0.063
0.064
0.093
0.071
0.073
0.073
0.064
0.065

h = 0.5
0.847
0.871
0.831
0.947
0.946
0.921
0.571
0.527
0.527
0.693
0.684
0.632
0.355
0.284
0.284
0.401
0.379
0.345
0.253
0.168
0.168
0.256
0.219
0.190

h=1
0.977
0.946
0.946
0.997
0.994
0.976
0.821
0.738
0.738
0.904
0.893
0.838
0.574
0.454
0.454
0.609
0.581
0.509
0.396
0.254
0.254
0.405
0.340
0.296

τ = 0.5
h = 0.5
0.990
0.990
0.979
0.999
0.999
0.996
0.990
0.876
0.842
0.949
0.942
0.921
0.703
0.666
0.633
0.738
0.714
0.688
0.597
0.532
0.482
0.585
0.547
0.503

h=1
1.000
1.000
1.000
1.000
1.000
1.000
0.994
0.990
0.975
0.999
0.998
0.994
0.931
0.905
0.857
0.948
0.933
0.903
0.832
0.772
0.729
0.839
0.802
0.762

τ = 0.25

τ = 0.5

h=0

h = 0.5

h=1

h = 0.5

h=1

0.056

0.820

0.912

0.984

0.999

0.061

0.920

0.970

0.996

1.000

0.070

0.529

0.702

0.856

0.982

0.076

0.663

0.820

0.940

0.995

0.072

0.297

0.428

0.640

0.875

0.069

0.369

0.510

0.715

0.920

0.073

0.165

0.236

0.499

0.738

0.068

0.199

0.296

0.529

0.782

ANNIKA BETKEN AND MARTIN WENDLER

Pareto(3, 1)
H = 0.6

18

Table 3. Rejection rates of the self-normalized Wilcoxon change-point test obtained by subsampling (left) with
block length l = bnγ c, γ ∈ {0.4, 0.5, 0.6}, and by comparison with asymptotic critical values (right) for Pareto(3,
1)-transformed fractional Gaussian noise of length n with Hurst parameter H.

Table 4. Rejection rates of the classical Wilcoxon change-point test obtained by subsampling (left) with block length
l = bnγ c, γ ∈ {0.4, 0.5, 0.6}, and by comparison with asymptotic critical values (right) for fractional Gaussian noise of
length n with Hurst parameter H.
sampling window method

asymptotic distribution

fGn
n
H = 0.6 300

500

H = 0.7 300

500

H = 0.8 300

500

H = 0.9 300

h=0
0.066
0.054
0.059
0.063
0.058
0.062
0.052
0.049
0.051
0.042
0.039
0.046
0.028
0.029
0.034
0.023
0.028
0.039
0.009
0.009
0.015
0.008
0.011
0.021

h = 0.5
0.20
0.223
0.264
0.285
0.345
0.397
0.080
0.095
0.120
0.104
0.131
0.160
0.030
0.038
0.057
0.031
0.044
0.071
0.010
0.014
0.029
0.006
0.009
0.021

h=1
0.232
0.411
0.529
0.436
0.663
0.789
0.088
0.158
0.227
0.153
0.267
0.373
0.031
0.048
0.088
0.036
0.070
0.129
0.006
0.009
0.028
0.003
0.011
0.032

h = 0.5
0.386
0.439
0.663
0.569
0.627
0.683
0.162
0.206
0.267
0.249
0.287
0.343
0.054
0.075
0.070
0.064
0.097
0.137
0.016
0.021
0.011
0.015
0.029
0.058

h=1
0.591
0.784
0.870
0.856
0.952
0.975
0.302
0.466
0.593
0.539
0.689
0.789
0.092
0.179
0.272
0.162
0.273
0.391
0.020
0.060
0.153
0.026
0.086
0.197

τ = 0.25

τ = 0.5

h=0

h = 0.5

h=1

h = 0.5

h=1

0.026

0.096

0.160

0.223

0.727

0.036

0.148

0.256

0.378

0.897

0.035

0.067

0.228

0.167

0.665

0.030

0.079

0.259

0.225

0.714

0.077

0.153

0.421

0.245

0.673

0.050

0.112

0.439

0.226

0.714

0.36

0.484

0.739

0.524

0.830

0.319

0.439

0.743

0.511

0.845

19

500

l
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41

τ = 0.5

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

τ = 0.25

sampling window method

asymptotic distribution

τ = 0.25
n
300

500

H = 0.7

300

500

H = 0.8

300

500

H = 0.9

300

500

l
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41
9
17
30
12
22
41

h=0
0.170
0.130
0.109
0.163
0.132
0.114
0.224
0.175
0.140
0.208
0.167
0.143
0.203
0.160
0.137
0.190
0.160
0.137
0.128
0.097
0.092
0.112
0.100
0.095

h = 0.5
0.949
0.963
0.962
0.991
0.997
0.997
0.785
0.802
0.789
0.921
0.931
0.925
0.508
0.496
0.484
0.639
0.649
0.626
0.150
0.128
0.125
0.159
0.161
0.170

h=1
0.742
0.861
0.871
0.916
0.976
0.989
0.568
0.680
0.708
0.763
0.862
0.891
0.326
0.347
0.364
0.445
0.513
0.556
0.077
0.071
0.077
0.089
0.101
0.106

τ = 0.5
h = 0.5
0.991
0.996
0.998
1.000
1.000
1.000
0.939
0.955
0.959
0.989
0.992
0.994
0.743
0.776
0.791
0.865
0.886
0.890
0.320
0.403
0.481
0.402
0.518
0.571

h=1
0.923
0.991
0.998
0.993
0.999
1.000
0.796
0.949
0.976
0.956
0.996
0.998
0.565
0.808
0.881
0.770
0.929
0.961
0.336
0.550
0.677
0.436
0.680
0.771

τ = 0.25

τ = 0.5

h=0

h = 0.5

h=1

h = 0.5

h=1

0.108

0.938

0.985

0.998

1.000

0.128

0.988

0.999

1.000

1.000

0.179

0.833

0.969

0.974

0.999

0.191

0.940

0.994

0.996

1.000

0.204

0.729

0.925

0.918

0.993

0.212

0.805

0.963

0.948

0.999

0.309

0.712

0.901

0.848

0.966

0.27

0.726

0.911

0.851

0.975

ANNIKA BETKEN AND MARTIN WENDLER

Pareto(3, 1)
H = 0.6

20

Table 5. Rejection rates of the classical Wilcoxon change-point test obtained by subsampling (left) with block length
l = bnγ c, γ ∈ {0.4, 0.5, 0.6}, and by comparison with asymptotic critical values (right) for Pareto(3, 1)-transformed
fractional Gaussian noise of length n with Hurst parameter H.

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

21

5. Proofs
5.1. Auxiliary Results.
Lemma 1. Under Assumption 2, there is a constant KD < ∞, such that for all
Pl
x1 , . . . , xl ∈ R with Var( i=1 xi ξi ) = 1
l
X

x2i ≤ KD .

i=1

Proof. Recall that we can rewrite the covariances as
Z π
eikλ f (λ)dλ
γ(k) =
−π

and that the spectral density f can be written as f (λ) = Lf (|λ|)|λ|D−1 . By our
assumptions Lf (x) ≥ Cmin for a constant Cmin > 0, so that we can conclude that
1 =Var

l
X


xi ξi =

i=1

=2

π

π

−π

X

X

ei(j−k)λ f (λ)dλ =

xj xk

1≤j,k≤l

Z

xj xk γ(j − k)

1≤j,k≤l

Z

X

=

X

1≤j,k≤l

0 1≤j,k≤l

≥2Cmin π D−1

Z

π

0
π

0

l
X

ei(j−k)λ Lf (|λ|)|λ|D−1 dλ

−π

xj xk ei(j−k)λ Lf (λ)λD−1 dλ = 2
Z

π

Z
xj xk

l
X

2

xj e−ijλ Lf (λ)λD−1 dλ

j=1

2

xj e−ijλ dλ.

j=1

We rewrite the integrand as
l
X

2

xj e−ijλ

X

=

j=1

xj xk e−ijλ eikλ =

j=1

1≤j,k≤l

=

l
X

x2j +

j=1

=

l
X

l
X

X

x2j +

X

xj xk e−i(j−k)λ

j6=k



xj xk e−i(j−k)λ + e−i(k−j)λ



j<k

X

x2j + 2

j=1

xj xk cos((k − j)λ) =

j<k

X
1≤j,k≤l

As a result, we have
Z
0

π

l
X

2

xj e−ijλ dλ =

π

Z
0

j=1

=

X

xj xk cos((k − j)λ)dλ

1≤j,k≤l

X
1≤j,k≤l

Z

π

cos((k − j)λ)dλ

xj xk
0

xj xk cos((k − j)λ).

22

ANNIKA BETKEN AND MARTIN WENDLER

=

l
X

x2j

cos(0)dλ +
0

j=1

=π

π

Z

l
X

X

Z

π

cos((k − j)λ)dλ

xj xk
0

j6=k

x2j .

j=1

All in all, this yields
1 = Var

l
X

!
xi ξi

≥ 2Cmin π

D−1

π

Z
0

i=1

l
X

2

xj e−ijλ dλ = 2Cmin π D

l
X

x2j .

j=1

j=1
D

Therefore, the statement of the lemma holds with KD = 1/(2Cmin π ).

0
< ∞ and l0 ∈ N such
Lemma 2. Under Assumption 2, there are constants KD
that
l
X
0 D/2
l
xi ≤ K D
i=1

for all l ≥ l0 and x1 , . . . , xl ∈ R with Var

P

l
i=1


xi ξi = 1.

Proof. The statement of the proof is equivalent to the existence of a constant C > 0,
Pl
such that for all x1 , . . . , xl ∈ R with i=1 xi = 1, we have
!
l
X
Var
xi ξi ≥ Cl−D .
i=1

P

l
?
?
x
=
1
be
the
values
that
minimize
Var
x
ξ
i=1 i
i=1 i i .
Pl
?
Then µ̂ξ (ξ1 , . . . , ξn ) :=
i=1 xi ξi is the best linear unbiased estimator for µ :=
E(ξ1 ). For a process (ζn )n∈N with spectral density
Let x?1 , . . . , x?l ∈ R with

Pl

fζ (x) =

1
1 − eix
2π

D−1

,

we have
Var (µ̂ζ (ζ1 , . . . , ζn )) ≥ C1 l−D
for a constant C1 > 0 by a Corollary of Adenstedt [1] (see p. 1101). We rewrite the
spectral density fζ of (ζn )n∈N with the help of the spectral density f of (ξn )n∈N as
D−1

fζ (x) = f (x)

1 − eix
.
2π|x|D−1 Lf (x)
D−1

|1−eix |
Note that the function g with g(x) = 2π|x|D−1 Lf (x) is bounded, as we assumed that
Lf is bounded away from 0. Hence, we have
Var (µ̂ξ (ξ1 , . . . , ξn )) ≥
for all l ≥ l0 by Lemma 4.4 in [1].

1
Var (µ̂ζ (ζ1 , . . . , ζn )) ≥ Cl−D
g(0)


SUBSAMPLING UNDER LONG RANGE DEPENDENCE

23

The next Lemma deals with the ρ-mixing coefficient, which is defined in the
following way: Let A, B be two σ-fields. Then
ρ(A, B) := sup corr(X, Y ),
where the supremum is taken over all A-measurable random variables X and all
B-measurable random variables Y . For details we recommend the book of Bradley
[18].
Lemma 3. Under Assumption 2, there are constants C1 , C2 < ∞ such that

ρ(k, l) := ρ σ(ξi , 1 ≤ i ≤ l), σ(ξj , k + l + 1 ≤ j ≤ k + 2l)
−D

≤ C1 (k/l)

Lγ (k) + C2 l2 k −D−1 max{Lγ (k), 1}

for all k ∈ N and all l ∈ {lk , . . . , k}.
Proof. Kolmogorov and Rozanov [34] proved that there exist real numbers a1 , a2 , . . . , al ,
b1 , b2 , . . . , bl such that
l
l
X

X

ai ξi ,
bj ξk+l+j
ρ σ(ξi , 1 ≤ i ≤ l), σ(ξj , k + l + 1 ≤ j ≤ k + 2l) = Cov
i=1

and Var

Cov

Pl



i=1

l
X

ai ξi = Var

ai ξi ,

l
X

Pl

j=1 bj ξk+l+j

bj ξk+l+j



j=1

= 1. The triangular inequality yields



j=1

i=1

≤

l
X
i=1

ai

l
X

bj |γ(k)| +

j=1

l X
l
X

|ai ||bj | |γ(k) − γ(k + l + j − i)| .

i=1 j=1

We will treat the two summands on the right hand side separately. For the first
term, it follows by Lemma 2 that
l
X
i=1

ai

l
X

bj |γ(k)| =

l
X
i=1

j=1

ai

l
X

bj |γ(k)| ≤ Kd02 lD Lγ (k)k −D .

j=1

Before we deal with the second summand, we observe that by Hölder’s inequality
and Lemma 1
v
v
u l
u l
l
l
u X
X
X
p
p
u X
√
√
t
2
|ai | ≤ l
ai ≤ KD l and
|bj | ≤ tl
b2j ≤ KD l.
i=1

i=1

j=1

j=1

Due to Assumption 2
sup
|k−k̃|≤2l−1

for some constant K.

Lγ (k) − Lγ (k̃) ≤ K

l
k

24

ANNIKA BETKEN AND MARTIN WENDLER

Consequently, for all k̃ ∈ {k + 1, . . . , k + 2l − 1}
γ(k) − γ(k̃) ≤ Lγ (k) k −D − k̃ −D + |Lγ (k) − Lγ (k̃)|k̃ −D

≤ Lγ (k) k −D − (k + 2l − 1)−D + |Lγ (k) − Lγ (k̃)|k −D
l
≤ Cd k −D−1 lLγ (k) + K k −D max{Lγ (k), 1}
k
≤ C3 k −D−1 l max{Lγ (k), 1}
Pl
Pl
for some constants Cd , C3 . Combining this with the bounds for i=1 |ai |, j=1 |bj |,
we finally arrive at
l
X
i=1

|ai |

l
X

|bj | |γ(k) − γ(k + l + j − i)| ≤ KD l

max

γ(k) − γ(k̃)

k̃∈{k+1,...,k+2l−1}

j=1

= KD C3 k −D−1 l2 max{Lγ (k), 1}.

5.2. Proof of the Main Result. Let t be a point of continuity of FT . In order
to simplify notation, we write N = n − l + 1 and Tl,i = Tl (Xi , . . . , Xi+l−1 ). The
triangular inequality yields
|F̂l,n (t) − FTn (t)| ≤ |F̂l,n (t) − FT (t)| + |FT (t) − FTn (t)|.
The second term on the right-hand side of the above inequality converges to zero
because of Assumption 1. As L2 -convergence implies stochastic convergence, it
suffices to show that


E |F̂l,n (t) − FT (t)|2 −→ 0
in order to prove that the first term converges to zero, as well. We have


E |F̂l,n (t) − FT (t)|2
 
2

2

2
2
(t) − E F̂l,n (t) + (FT (t)) − 2FT (t) E F̂l,n (t) + E F̂l,n (t)
= E F̂l,n
2

= Var(F̂l,n (t)) + E F̂l,n (t) − FT (t) .
Furthermore, stationarity of the process (Xn )n∈N and Assumption 1 imply
E F̂l,n (t) =

N

1 X
l→∞
E 1{Tl,i ≤t} = P (Tl,1 ≤ t) = FTl (t) −−−→ FT (t).
N i=1

It remains to show that Var(F̂l,n (t)) −→ 0. Again, it follows by stationarity of
(Xn )n∈N that
N




1
2 X
Var F̂l,n (t) = Var 1{Tl,1 ≤t} + 2
(N − i + 1)Cov 1{Tl,1 ≤t} , 1{Tl,i ≤t}
N
N i=2

≤

N

2 X
Cov 1{Tl,1 ≤t} , 1{Tl,i ≤t} .
N i=1

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

25

Recall that by Assumption 3, we have l ≤ Cl n(1+D)/2− for some constants Cl and
 > 0. For n large enough such that l < 12 bn1−/2 c, we split the sum of covariances
into two parts:
N

1 X
Cov 1{Tl,1 ≤t} , 1{Tl,i ≤t}
N i=1

1
=
N
≤

≤

bn1−/2 c

X

Cov 1{Tl,1 ≤t} , 1{Tl,i ≤t}



i=1

bn1−/2 c
1
+
N
N
bn1−/2 c
1
+
N
N

N
X

+

1
N

N
X

Cov 1{Tl,1 ≤t} , 1{Tl,i ≤t}



i=bn1−/2 c+1

ρ(σ(Xi , 1 ≤ i ≤ l), σ(Xj , k ≤ j ≤ k + l − 1))

k=bn1−/2 c+1
NX
−l−1

ρ(k, l),

k=bn1−/2 c−l

where

ρ(k, l) := ρ σ(Xi , 1 ≤ i ≤ l), σ(Xj , k + l + 1 ≤ j ≤ k + 2l) .
Obviously, the first summand converges to zero by Assumption 3. For the second
summand note that as a consequence of Potter’s Theorem (Theorem 1.5.6 in the
book of Bingham, Goldie and Teugels [17]), there is a constant CL such that Lγ (k) ≤
CL k D/2 for all k ∈ N. This together with Lemma 3 yields
1
N

NX
−l−1

≤ CL C1
≤

ρ(k, l)

k=bn1−/2 c−l

lD
N

NX
−l−1

k −D k D/2 + CL C2

k=bn1−/2 c/2

l2
N

NX
−l−1

k −D−1 k D/2

k=bn1−/2 c/2

CL C1 ClD 2D(1−/2) nD(((1+D)/2−)−(1−/2)+/2(1−/2))

+ CL C2 Cl2 21+D(1−/2) n((1+D−2)−(D+1)(1−/2)+(1−/2)D/2)


2
3
n→∞
≤ C n−D((1−D)/2+ /4) + n−( 2 −D+D/4) −−−−→ 0
for some constant C < ∞. Thus, we have proved that Var(F̂l,n (t)) → 0 as n → ∞
and that the first conjecture of Theorem 1 holds.
The second assertion of Theorem 1 follows from
P

FTn (t) − F̂l,n (t) −→ 0
by the usual Glivenko-Cantelli argument for the uniform convergence of empirical
distribution functions; see for example section 20 in the book of Billingsley [16]. 
Acknowledgements
We thank the referee for his careful reading of the article and his thoughtful
comments which lead to a significant improvement of the article. We also thank
Norman Lambot for reading the article, thereby helping to reduce the number of
misprints.

26

ANNIKA BETKEN AND MARTIN WENDLER

References
1. Rolf K. Adenstedt, On large-sample estimation for the mean of a stationary
random sequence, The Annals of Statistics (1974), 1095–1107.
2. Donald W. K. Andrews, Tests for parameter instability and structural change
with unknown change point, Econometrica 61 (1993), 821–856.
3. Jaromı́r Antoch, Marie Hušková, Alicja Janic, and Teresa Ledwina, Data driven
rank test for the change point problem, Metrika 68 (2008), 1–15.
4. Changryong Baek and Vladas Pipiras, On distinguishing multiple changes in
mean and long-range dependence using local Whittle estimation, Electronic
Journal of Statistics 8 (2014), no. 1, 931–964.
5. Shuyang Bai and Murad S. Taqqu, Canonical correlation between blocks
of long-memory time series and consistency of subsampling, arXiv preprint
arXiv:1512.00819 (2015).
6. Shuyang Bai, Murad S. Taqqu, and Ting Zhang, A unified approach to selfnormalized block sampling, Stochastic Processes and their Applications 126
(2016), no. 8, 2465–2493.
7. Nathan S. Balke, Detecting level shifts in time series, Journal of Business &
Economic Statistics 11 (1993), no. 1, 81–92.
8. J. Beran and Y. Feng, SEMIFAR models - a semiparametric framework for
modelling trends, long-range dependence and nonstationarity, Computational
Statistics & Data Analysis 40 (2002), no. 2, 393–419.
9. Jan Beran, Statistics for long-memory processes, Chapman & Hall, 1994.
10. Jan Beran, Yuanhua Feng, Sucharita Ghosh, and Rafal Kulik, Long-memory
processes, Springer-Verlag Berlin Heidelberg, 2013.
11. István Berkes, Lajos Horváth, Piotr Kokoszka, and Qi-Man Shao, On discriminating between long-range dependence and changes in mean, The Annals of
Statistics 34 (2006), 1140–1165.
12. Annika Betken, Testing for change-points in long-range dependent time series
by means of a self-normalized Wilcoxon test, Journal of Time Series Analysis
37 (2016), 185–809.
13. Eric Beutner and Henryk Zähle, Continuous mapping approach to the asymptotics of U- and V-statistics, Bernoulli 20 (2014), no. 2, 846–877.
14. Rajendra J. Bhansali and Piotr S. Kokoszka, Estimation of the long-memory
parameter: a review of recent developments and an extension, Lecture NotesMonograph Series (2001), 125–150.
15. Peter J. Bickel and Anat Sakov, On the choice of m in the m out of n bootstrap
and confidence bounds for extrema, Statistica Sinica (2008), 967–985.
16. Patrick Billingsley, Probability and measure, John Wiley & Sons, Inc., 1995.
17. N. H. Bingham, C. M. Goldie, and J. L. Teugels, Regular variation, Cambridge
University Press, 1987.
18. Richard C. Bradley, Introduction to strong mixing conditions, Kendrick press,
2007.
19. George W. Cobb, The problem of the Nile: conditional solution to a changepoint
problem, Biometrika 65 (1978), no. 2, 243 – 251.
20. Martial Coulon, Marie Chabert, and Ananthram Swami, Detection of multiple changes in fractional integrated ARMA processes, Signal Processing, IEEE
Transactions on 57 (2009), no. 1, 48–61.

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

27

21. Miklós Csörgő and Lajos Horváth, Limit theorems in change-point analysis,
Wiley Chichester; New York, 1997.
22. Herold Dehling, Aeneas Rooch, and Murad S. Taqqu, Non-parametric changepoint tests for long-range dependent data, Scandinavian Journal of Statistics 40
(2013), 153 – 173.
23. Herold Dehling, Aeneas Rooch, and Martin Wendler, Two-Sample U-Statistic
Processes for Long-Range Dependent Data, arXiv preprint arXiv:1404.0551
(2014).
24. Rohit S. Deo and Clifford M. Hurvich, Linear trend with fractionally integrated
errors, Journal of Time Series Analysis 19 (1998), no. 4, 379–397.
25. Roland L. Dobrushin and Peter Major, Non-central limit theorems for nonlinear functionals of Gaussian fields, Zeitschrift für Wahrscheinlichkeitstheorie
und verwandte Gebiete 50 (1979), no. 1, 27–52.
26. Zhenwen Fan, Statistical issues and developments in time series analysis and
educational measurement, BiblioBazaar, 2012.
27. Friedrich Götze and Alfredas Račkauskas, Adaptive choice of bootstrap sample
sizes, Lecture Notes-Monograph Series (2001), 286–309.
28. Clive W. J. Granger and Roselyne Joyeux, An introduction to long-memory
time series models and fractional differencing, Journal of Time Series Analysis
1 (1980), no. 1, 15–29.
29. Peter Hall, Bing-Yi Jing, and Soumendra Nath Lahiri, On the sampling window
method for long-range dependent data, Statistica Sinica 8 (1998), no. 4, 1189–
1204.
30. Peter Hall and Bingyi Jing, On sample reuse methods for dependent data, Journal of the Royal Statistical Society. Series B (Methodological) (1996), 727–737.
31. Harold Edwin Hurst, Methods of using long-term storage in reservoirs, ICE
Proceedings, vol. 5, Thomas Telford, 1956, pp. 519–543.
32. Ildar Abdulovich Ibragimov and Yurii Antolevich Rozanov, Gaussian random
processes, Springer, New York, 1978.
33. Agnieszka Jach, Tucker McElroy, and Dimitris N. Politis, Subsampling inference
for the mean of heavy-tailed long-memory time series, Journal of Time Series
Analysis 33 (2012), no. 1, 96–111.
34. A. N. Kolmogorov and Yu. A. Rozanov, On strong mixing conditions for stationary Gaussian processes, Theory of Probability & Its Applications 5 (1960),
no. 2, 204–208.
35. Hans R. Künsch, Statistical aspects of self-similar processes, Proceedings of
the first World Congress of the Bernoulli Society, vol. 1, VNU Science Press
Utrecht, The Netherlands, 1987, pp. 67–74.
36. Soumendra Nath Lahiri, On the moving block bootstrap under long range dependence, Statistics & Probability Letters 18 (1993), no. 5, 405–413.
37. Will E. Leland, Murad S. Taqqu, Walter Willinger, and Daniel V. Wilson,
On the self-similar nature of Ethernet traffic (extended version), Networking,
IEEE/ACM Transactions on 2 (1994), no. 1, 1–15.
38. Will E. Leland and Daniel V. Wilson, High time-resolution measurement and
analysis of LAN traffic: Implications for LAN interconnection, INFOCOM’91.
Proceedings. Tenth Annual Joint Conference of the IEEE Computer and Communications Societies. Networking in the 90s., IEEE, IEEE, 1991, pp. 1360–
1366.

28

ANNIKA BETKEN AND MARTIN WENDLER

39. Andrew W. Lo, Long-term memory in stock market prices, Tech. report, National Bureau of Economic Research, 1989.
40. Ignacio N. Lobato, Testing that a dependent process is uncorrelated, Journal of
the American Statistical Association 96 (2001), 1066–1076.
41. I. B. Macneill, S. M. Tang, and V. K. Jandhyala, A Search for the Source of
the Nile’s Change-Points , Environmetrics 2 (1991), no. 3, 341 – 375.
42. Tucker McElroy and Dimitris Politis, Self-normalization for heavy-tailed time
series with long memory, Statistica Sinica 17 (2007), no. 1, 199.
43. Daniel J. Nordman and Soumendra N. Lahiri, Validity of the sampling window method for long-range dependent linear processes, Econometric Theory 21
(2005), no. 06, 1087–1111.
44. Vladas Pipiras and Murad S. Taqqu, Long-range dependence and self-similarity,
Cambridge University Press, 2011.
45. Dimitris N. Politis and Joseph P. Romano, Large sample confidence regions
based on subsamples under minimal assumptions, The Annals of Statistics
(1994), 2031–2050.
46. Dimitris N. Politis, Joseph P. Romano, and Michael Wolf, Subsampling, 1999.
47. Peter M. Robinson, Gaussian semiparametric estimation of long range dependence, The Annals of Statistics (1995), 1630–1661.
48. Xiaofeng Shao, A simple test of changes in mean in the possible presence of
long-range dependence, Journal of Time Series Analysis 32 (2011), 598–606.
49. Xiaofeng Shao and Xianyang Zhang, Testing for change points in time series,
Journal of the American Statistical Association 105 (2010), 1228–1240.
50. Michael Sherman and Edward Carlstein, Replicate histograms, Journal of the
American Statistical Association 91 (1996), no. 434, 566–576.
51. Yu G. Sinai, Self-similar probability distributions, Theory of Probability & Its
Applications 21 (1976), no. 1, 64–80.
52. D. Surgailis, Zones of attraction of self-similar multiple integrals, Lithuanian
Mathematical Journal 22 (1982), no. 3, 327–340.
53. Murad S. Taqqu, Convergence of integrated processes of arbitrary Hermite rank,
Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 50 (1979),
53–83.
54. Murad S. Taqqu and Vadim Teverovsky, Robustness of Whittle-type estimators for time series with long-range dependence, Communications in statistics.
Stochastic models 13 (1997), no. 4, 723–757.
55. Lihong Wang, Gradual changes in long memory processes with applications,
Statistics 41 (2007), no. 3, 221–240.
, Change-point detection with rank statistics in long-memory time-series
56.
models, Australian & New Zealand Journal of Statistics 50 (2008), 241–256.
57. Wei Biao Wu and Zhibiao Zhao, Inference of trends in time series, Journal
of the Royal Statistical Society: Series B (Statistical Methodology) 69 (2007),
no. 3, 391 – 410.
58. Ting Zhang, Hwai-Chung Ho, Martin Wendler, and Wei Biao Wu, Block sampling under strong dependence, Stochastic Processes and their Applications 123
(2013), no. 6, 2323–2339.

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

29

Appendix A. A Modified Change Point Test for Data with Ties
If the distribution of Xi = G(ξi ) is not continuous, there is a positive probability
that Xi = Xj for some i 6= j, so there might be ties in the sample. We propose to
Pn
use the following test statistic based on the modified ranks R̃i = j=1 (1{Xj <Xi } +
1
2 1{Xj =Xi } ):
Pk
Pn
k
i=1 R̃i − n
i=1 R̃i
T̃n (τ1 , τ2 ) :=
max
n P
o1/2 ,
Pn
k∈{bnτ1 c,...,bnτ2 c}
k
1
2 (1, k) + 1
2 (k + 1, n)
S̃
S̃
t=1 t
t=k+1 t
n
n
where
S̃t (j, k) =

t 
X
R̃h −
h=j


k
X
1
R̃i .
k − j + 1 i=j

To be able to apply subsampling, we need T̃n to converge in distribution, which we
will show now:
Lemma 4. Let (ξn )n∈N be a stationary sequence of centered standard Gaussian
variables with covariance function γ(k) = k −D Lγ (k) for a D ∈ (0, 1) and a slowly
varying function Lγ . Let Xi = G(ξi ) for a function G, piecewise monotone on
finitely many pieces. Then T̃n (τ1 , τ2 ) ⇒ T for some random variable T .
Proof. Let h(x, y) = 1{G(x)<G(y)} + 12 1{G(x)=G(y)} − 12 . We define the modified
Wilcoxon process (W̃n (λ))λ∈[0,1] by
W̃n (λ) :=

[nλ]
1 X
ndn i=1

n
X

h(ξi , ξj )

j=[nλ]+1

p
Pn
with dn = Var( i=1 ξi ). From Theorem 2.2 in Dehling, Rooch, Wendler [23], we
have the weak convergence of this process W̃n to the limit process W with
W (λ)
Z
Z Z

= −(1 − λ)Z(λ) ϕ(x)dh̃(x) − λ(Z(1) − Z(λ))
ϕ(y)dh(x, y)(y) ϕ(x)dx.
Here, Z is a fractional Brownian motion, ϕ is the density function of the standard
normal distribution and h̃(x) = E[h(x, ξi )]. Following the proof of Theorem 1 in
[12], we can express T̃n (τ1 , τ2 ) as a function of W̃n :
Tn (τ1 , τ2 )
=

sup
τ1 ≤λ≤τ2

Rλ
0

W̃n (λ)
R1
(t)
(W̃n (t) − ccnn(λ)
W̃n (λ))2 dt + λ (W̃n (t) −

1/2
1−cn (t)
2
1−cn (λ) W̃n (λ)) dt

.

Note that cn (λ) converges to λ uniformly, so we have the asymptotic equivalence
Tn (τ1 , τ2 )
≈

sup
τ1 ≤λ≤τ2

Rλ
0

W̃n (λ)
R1
(W̃n (t) − λt W̃n (λ))2 dt + λ (W̃n (t) −

1/2
1−t
2
1−λ W̃n (λ)) dt

.

30

ANNIKA BETKEN AND MARTIN WENDLER

By the continuous mapping theorem, we get
Tn (τ1 , τ2 )
⇒

sup
τ1 ≤λ≤τ2

Rλ
0

(W (t) −

t
2
λ W (λ)) dt

|W (λ)|
R1
+ λ (W (t) −

1/2
1−t
2
1−λ W (λ)) dt

=: T.


Appendix B. A Test for Multiple Change Points
For testing the alternative hypothesis of two change points, we suggest to use
the test statistic Tn (τ1 , τ2 , ε) = sup(k1 ,k2 )∈Ωn (τ1 ,τ2 ,ε) |Gn (k1 , k2 )|. Some calculations
yield
Gn (k1 , k2 )
W̃n (λ1 , λ2 )
=

λ
R1
0

2
W̃n (r, λ2 )− λr1 W̃n (λ1 , λ2 ) dr

+
λ
R2
λ1

+

λ
R2
λ1

2
−r
W̃n (r, λ2 )− λλ22−λ
W̃n (λ1 , λ2 ) dr
1

 12

|Wn∗ (λ2 , λ1 )|
2
1
Wn∗ (r, λ1 )− λr−λ
Wn∗ (λ2 , λ1 ) dr
2−λ1

+

R1
λ2

2
1−r
Wn∗ (λ2 , λ1 ) dr
Wn∗ (r, λ1 )− 1−λ
2

 12

+ oP (1),
where
W̃n (λ, τ ) := Wn (λ, λ) − Wn (λ, τ ),

Wn∗ (λ, τ ) := Wn (λ, λ) − Wn (τ, λ)

with
bnλc

Wn (λ, τ ) =

X

n
X





1{Xi ≤Xj } −

1
2





i=1 j=bnτ c+1

, 0 ≤ λ ≤ τ ≤ 1.

Define
d2n := Var 

n
X

Hr (ξj ) ,

j=1

where Hr denotes the r-th order
 Hermite polynomial and r designates the Hermite
rank of the class of functions 1{G(ξi )≤x} − F (x), x ∈ R . It can be shown that
1
ndn Wn (λ, τ ) converges in distribution to
Z
1
{(1 − τ )Zr (λ) − λ(Zr (1) − Zr (τ ))}
Jr (x)dF (x), 0 ≤ λ ≤ τ ≤ 1,
r!
where Zr is an r-th order Hermite process with Hurst parameter H := max{1 −
rD 1
2 , 2 } and where

Jr (x) = E Hr (ξi )1{G(ξi )≤x} .

SUBSAMPLING UNDER LONG RANGE DEPENDENCE

31

As a result, under the hypothesis the limiting distribution of Tn (τ1 , τ2 , ε) is given
by T (r, τ1 , τ2 , ε) = supτ1 ≤λ1 <λ2 ≤τ2 , λ2 −λ1 ≥ε Gr (λ1 , λ2 ) with
Gr (λ1 , λ2 )
Zr (λ1 ) −
=

R λ1
0

Zr (t) −

2
t
λ1 Zr (λ1 ) dt

+

R λ2
λ1

Zr (t) −

Zr (λ2 ) −
+
λ
R2
λ1

Zr (t) −

λ2 −t
λ2 −λ1 Zr (λ1 )

−

λ1
λ2 Zr (λ2 )
t−λ1
λ2 −λ1 Zr (λ2 )

1−λ2
1−λ1 Zr (λ1 )

2
t−λ1
λ2 −λ1 Zr (λ2 ) dt

+

−

R1
λ2

Ruhr-Universität Bochum, Germany
E-mail address: annika.betken@rub.de
Ernst-Moritz-Arndt-Universität Greifswald, Germany
E-mail address: martin.wendler@uni-greifswald.de

 12

−

2
λ2 −t
λ2 −λ1 Zr (λ1 ) dt

λ2 −λ1
1−λ1 Zr (1)

Zr (t) −

1−t
1−λ2 Zr (λ2 )

 21 .

−

2
t−λ2
1−λ2 Zr (1) dt

