Bootstrap for change point detection
Nazar Buzun, Valeriy Avanesov
{buzun,avanesov}@wias-berlin.de

arXiv:1710.07285v1 [math.ST] 19 Oct 2017

October 23, 2017

1

Introduction

The problem of change point detection appears each time one needs to explore a set of random
data and make a decision about homogeneity of its structure. In other words, the problem can
be stated as two following questions: were there any structural changes in the nature of observed
data? At which moments, if so? The present work mainly focuses on the sequential or online
change point detection. In this case the data is aggregated from running random process. Formally
a time moment τ is a change point, if stochastic properties of the observed signal {Yt }nt=1 have
undergone changes in its distribution:
(
Yt v IP1 t < τ,
Yt v IP2 t ≥ τ.
The goal is to find such structural breaks as soon as possible. Such problem arises across many
scientific areas: quality control Lai (1995), cybersecurity Blazek and Kim (2001), Wang et al.
(2004), econometrics Spokoiny (2009), Mikosch and Starica (2004), geodesy e.t.c. Article Shiryaev
(1963) describes classical results in change point detection theory. Overview of the state-of-art
methods are presented in Polunchenko and Tartakovsky (2011) and Shiryaev (2010).
This research considers sequential hypothesis testing, in which each hypothesis (IP1 = IP2 )
monitors the presence of change point through Likelihood Ratio Test (LRT) using sliding window.
At each time step the procedure extracts a data slice, splits it in two parts of equal size and
executes LRT on it. High values of LRT indicate possible distribution difference in the window
parts (IP1 6= IP2 ). Procedures with LRT are rather popular in related literature. The work
Quandt (1960) proposes application of LRT for detection of breaks in linear regression model.
It was further developed by many authors, e.g. Haccou et al. (1987), Srivastava and Worsley
(1986). Papers Liu et al. (2008), Zou et al. (2007) investigate LRT for change point detection
for nonparametric case. Nonparametric approaches are easily adaptable for complex data but
in general they need more information for model building than their parametric alternatives.
Introduction of parametric assumption: IP1 , IP2 ∈ {IP (θ) : θ ∈ IRp } allows to reduce the suffisient
number of observations as soon as IP (θ) has less degrees of freedom than nontapametric model.
The state-of-the-art review of parametric models based on LRT and its application to economics
and bio-informatics are presented by Chen and Gupta (2012). The paper Gombay (2000) explores
how LRT can be used for sequential change point detection in case IP (θ) is exponential family.
The LRT statistic requires its quantiles or critical values to be set from the signal data {Yt }nt=1 .
Many works are dedicated to asymptotic behaviour of LRT, e.g. Jandhyala and Fotopoulos (1999)
obtains lower and upper bounds for distribution of asymptotic maximum likelihood estimator. The
work Kim (1994) provides a very detailed study of its asymptotic behaviour in linear regression
models. Similar results for change in mean of a Gaussian process are given in Fotopoulos et al.
1

(2010). In Biau et al. (2016) an approach with Wiener process and Donsker–Prohorov Theorem
describes relatively general method for LRT-like statistics distribution approximation.
Instead of asymptotic distribution for LRT one may find a benefit of resampling and bootstrap.
This technique is popular, e.g. Frick et al. (2014), Spokoiny (2009), since it provides a way to
simulate a complex distribution of LRT statistic (for wide family of IP (θ)) through empirical data
distribution. Using bootstrap one can generate LRT[ statistic multiple times in order to obtain
quantile distribution of the initial LRT. Both LRT and LRT[ statistics have (ref. Sections 5.2.
and 5.3.) approximation with the following forms with high probability
LRT ≈ kξ + ∆k,

LRT[ ≈ kξ [ + ∆[ k.

(Qf)

Larger ∆ values correspond to more pronounced hypothesis rejection (more apparent changes
in data sequence). Argument ξ could be treated as a noise component. For LRT critical value
calibration one requires data without change points and consequently with ∆ = 0. Section 2
contains description of a modified LRT which enable the calibration even if data contains change
points.
The cornerstone of the novel change point detection procedure is the concept of change-point
pattern. The geometry of a pattern depends on a type of transition region between two distributions that the data obeys before and after a change respectively. Three examples are presented
at the Fig. 1. The triangle (spades) pattern appears in case of an abrupt transition from IP (θ1 )
to IP (θ2 ). A smooth transition between two distributions entails trapezium change-point pattern.
And a horn pattern appears due to an abrupt change in variance. Processing of a change-point
pattern instead of a single LRT-value allows to reduce noise influence ξ(t) and false-alarm rate.
The presence of change-point patterns is the corollary of (Qf) representation.
60
25.0

55

22.5

50
45

20.0

40
35

15.0

30

LRT statistics

LRT statistics

17.5

12.5
10.0

25
20
15

7.5

10
5

5.0

0
2.5

-5

0.0

-10
-15

-2.5

-20
0

100

200

300

400

500

600

700

800

900

1,000

0

250

500

750

time

1,000

1,250

1,500

time

20.0
17.5
15.0
12.5

LRT statistics

10.0
7.5
5.0
2.5
0.0
-2.5
-5.0
-7.5
-10.0
-12.5
-15.0
-17.5
0

100

200

300

400

500

600

700

800

900

1,000

time

Figure 1: Types of change point and the geometry of change-point patterns: triangle pattern
– abrupt mean transition, trapezium pattern – smooth mean transition, horn pattern – abrupt
variance transition.
In case of a single change point one may find the pattern position by maximising convolution
with a pattern function Pτ (t) (ref. Section 2 for details):
X
argmax
Pτ (t)kξ(t) + ∆(t)k.
τ

t

P
In order to set critical value correctly quantiles of the statistic maxτ t Pτ (t)kξ(t)k should
P
P
be close in distribution to quantiles of maxτ t Pτ (t)kξ [ (t)k. Assuming that ξ(t) ≈ t+h
i=t−h ξ i
2

(independent random vectors sum), we have made probability measures comparison using technique from article Chernozhukov et al. (2013b). In Section 5.8. we prove Bootstrap approximation
illustrating on the way useful mathematical concepts such as Linderberg telescoping sums, anticoncentration of normal vector, Slepian bridge and empirical covariance matrix deviations.
Section
P
5.9. extends the statements for statistics of type maxt Qt (ξ) (in particular maxτ t Pτ (t)kξ(t)k,
Section 5.10.).
The last part of this paper (Section 5.11.) contains some specification for aforementioned results
with generalised linear models (GLM).

2

Procedure

This section provides description of the Change Point Detection algorithm which employs
Likelihood Ratio Test (LRT). Let (IP (θ), θ ∈ IRp , L(θ) = log(∂ n IP (θ)/∂Y )) be a parametric
assumption about the nature of data inside the window (Yt−h , . . . , Yt+h−1 ) with central point t
and size 2h. Here and further we assume, that the observations {Yi }ni=1 are independent, so
X
L(θ, Y) =
li (θ).
(L)
i

Denote argmax of the Likelihood function and the “real” model parameter value as follows
θ∗ = argmax IEL(θ, Y).

θb = argmax L(θ, Y),
θ

θ

The algorithm sequentially computes LRT statistic (Th (t)) for each t in the sliding window procedure. The LRT statistic itself corresponds to the gain from window split into two parts (Yl , Yr ):
b Y),
Th (t) = L(θbl , Yl ) + L(θbr ; Yr ) − L(θ,
Yl = (Yt−h , . . . , Yt−1 ),

Yr = (Yt , . . . , Yt+h−1 ),

θbl = argmax L(θ, Yl ),

θbr = argmax L(θ, Yr )

θ

(T)

θ

According to the Theorem 3, encountering change point, statistic 2Th (t) ≈ kξ(t) + ∆(t)k2 starts
growing according to change point pattern type (for example spades, trapezium, horn, ref. the
Figure 1). In order to match pattern positions, the procedure monitors 2h values of the LRT
simultaneously and convolves them with each of the predefined pattern functions Pτ (t):
TPh (τ ) =

X

p
Pτ (t) 2Th (t).

(TP)

t

√
High values of TPh (τ ) correspond to a sufficient correlation of 2Th and Pτ (similar to the
dependence on t). The algorithm marks a time moment τ at a scale h as a change point, if the
test statistic TPh (τ ) exceeds a calibrated (by bootstrap procedure) critical value zh :
{τ is a change point } ⇔ {∃h : TPh (τ ) > zh }.
The greater window size h is chosen, the more probably the algorithm will mark τ as a change
point. Again, small windows may mark τ faster.
Weighted bootstrap procedure enables resampling of the statistic max1≤τ ≤n TPh (τ ) and thus
calculation of the critical value zh for the window size 2h. It generates a sequence of weighted
likelihood functions, where each element is a convolution of independent likelihood components
and weight vector (u[1 , . . . , u[n ):
X
L[ (θ, Y) =
u[i li (θ),
(Lb)
i

3

where {u[i }ni=1 are i.i.d. and u[i ∈ N (1, 1). At each weights generation one gets a new value of L[ (θ)
b fluctuations.
and its optimal parameter θ[ and thus bootstrap procedure enables to estimate L(θ)
The corresponding bootstrap LRT statistic is
Th[ (t) = L[ (θl[ , Yl ) + L[ (θr[ , Yr ) − sup{L[ (θ, Yl ) + L[ (θ + θbr − θbl , Yr )},

(Tb)

θ

θ[ = argmax L[ (θ, Y).
θ

Parameter (θbr − θbl ) is required for condition Th[ ≈ ξ [ (ref. Theorem 6). In this case one can
estimate max1≤τ ≤n TP[h (τ ) quantiles under the null hypothesis (∆[ (t) ∝ θbr (t) − θbl (t)) instead of
the false assumption (∆[ (t) = 0).
Empirical bootstrap version generates subsamples of data {Yk } from the complete dataset with
random independent indexes of size n. In this case
L (θ, Y) =

X

lk(i) (θ),

(Le)

i

where {k(i)}ni=1 are i.i.d. and k(i) ∈ {1, . . . , n}. For all window positions θbr = θbl = θb and here
bias correction is not required. So the corresponding LRT statistic is like (T):
Th (t) = L (θl , Yl ) + L (θr , Yr ) − L (θ , Y),

(Te)

θ = argmax L (θ, Y).
θ

Empirical bootstrap works better in the application but less suitable for theoretical investigations
(the distribution is discontinuous).

3

Main results

Below we present the Theorems that describes difference between probabilistic measures of
TPh (τ ) and TP[h (τ ) (precision of the bootstrap calibration) and LRT sensitivity to parameter θ∗
transition at change point. In independent models each noise vector ξ lr (t) = ξ(t) ∈ IRp is a sum
of independent vectors (ref. Section 5.2. for ξ lr (t) definition)
ξ lr (t) =

t−1
X

ξi −

ξi ,

ξ i ∝ ∇li (θ∗ ).

i=t

i=t−h

Aggregate all ξ i into one vector

t+h−1
X

ξ T = (ξ T1 , . . . , ξ Tn ).

Theorem 1. Let dataset size be n and the window equal to h. Include conditions from lemmas
3, 6 and 9. Then for each fixed z




[
[
IP max TPh (τ ) > z − IP
max TPh (τ ) > z ≤ ∆T V
1≤τ ≤n

1≤τ ≤n

[
∆T V = C1 µZ + C2 k Var(ξ lr ) − Var [ (ξ lr[ )k1/2
∞ + 7CA (♦ + ♦ ),

where
p
√
k Var(ξ lr ) − Var [ (ξ lr[ )k∞ ≤ 10 log(np) 2hk Var(ξ)k∞ (3 + kbk) + kbk2 ,
2

kbk = max
t

t+2h
X
i=t

4

kIEξ i k2∞ ,

µ3Z ≤ 2hIE kξk3∞ ,
1/2

C1 = 5Cµ1/3 CA ,

C2 = 4CΣ CA .

Constants Cµ ∼ log2 (n), CΣ ∼ log(n) and CA ≈ p3/2 log(n) are described in Section 5.10..
The proof is a direct consequence of Theorems 3, 6 and 9.
Remark.

1. Parameters asymptotic
µZ ∼

log1/2 (n)
,
(2h)1/6

k Var(ξ lr ) − Var [ (ξ lr[ )k1/2
∞ ∼

log1/4 (n)
,
(2h)1/4

♦ + ♦[ ∼

p
h1/2

.

2. For quantile estimation of the statistic max1≤τ ≤n TPh (τ ) with quantile of max1≤τ ≤n TP[h (τ )
one has to show that


[
IP max TPh (τ ) > z (α) − α ≤ ∆T V ,
1≤τ ≤n

for z [ (α) defined by equation
IP

[


max

1≤τ ≤n

TP[h (τ )


> z (α) = α.
[

This statement is a consequence of the Theorem (1) but not a direct one since the argument
z [ (α) is random and depends on max1≤τ ≤n TPh (τ ). Involving sandwich Lemma 21 fulfills
this issue.
The next part of this Section evaluates the smallest parameter θ∗ transition that is sufficient
for change point detection in a fixed position τ and window size 2h. Let zh (α) be a quantile of
P
t Pτ (t)kξ lr (t)k such that
!
P
X
Var( t Pτ (t)kξ lr (t)k)
P
.
IP
Pτ (t)kξ lr (t)k > zh (α) = α ≤
(zh (α) − t Pτ (t)IEkξ lr (t)k)2
t
Section 5.17. provides upper bound for α and is summarized in following statement.
P
√
Theorem 2. Let t Pτ (t) = 0 and IEkξ lr (t)k = p. The sufficient condition for abrupt type
change point detection of size ∆ with probability 1 − e−x in position τ using triangle pattern (P) is
kDlr (θr∗ − θl∗ )(τ )k = ∆ > 5p1/4 (x + log(2h))1/4 ex/2 + 21♦,
where matrix Dlr is defined in Theorem 3.

4

Experiments

In order to substantiate patterns utility we compare procedure from Section 2 with the similar
one but without pattern (i.e. Pτ (t) = 1I[τ = t]). The experiment scenario is follows. The dataset
{Yi } consists of 500 normal random vectors from IR5 with one change point at position τ ∗ = 250.
Yi ∈ N (0, I5 ),
Yi ∈ N (0.25, I5 ),

0 ≤ i < 250
250 ≤ i < 500.

The procedure searches for change point location as τb = argmaxτ TPh (τ ). Then the quality of the
detection is measured by average error |b
τ − τ ∗ | (c.p. position error) and fraction of the detected
change points (POWER) (ref. Figure 2).
5

35.0

0.85

32.5

0.80

30.0

0.75
0.70
0.65

25.0

0.60

22.5

0.55

POWER

c.p. position error

27.5

20.0

0.50
0.45

17.5

0.40

15.0

0.35
0.30

12.5

0.25

10.0

0.20
0.15

7.5

0.10
80

90

100

110

120

130

140

150

160

170

180

190

200

210

220

230

240

80

90

100

110

120

130

140

window size
without pattern

150

160

170

180

190

200

210

220

230

240

window size

triangle pattern

without pattern

triangle pattern

Figure 2: Change point localisation test and power test for the procedure from Section 2. One
case with triangle pattern and the other case without pattern.
The second experiment describes bootstrap convergence depending on window size (2h). We set
bootstrap confidence level equal to 0.1 and compute p-value from real distribution with bootstrap
quantile z [ .


IP [

max TP[h (τ ) > z [

1≤τ ≤n


IP

max TPh (τ ) > z

[

1≤τ ≤n

= 0.1,




− 0.1 = O

1
hβ


.

From the plot below (ref. Figure 3) one can observe that
1
β> ,
2
which suppose better convergence in comparison with the theoretical study (ref. Theorem 1),
where β = 1/6.
-2.6
-2.7
-2.8
-2.9

log(P - Pb)

-3.0
-3.1
-3.2
-3.3
-3.4
-3.5
-3.6
-3.7
-3.8
3.0

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

3.9

4.0

4.1

4.2

4.3

4.4

4.5

4.6

log(2h)
homogeneous data

data with change point

Figure 3: Bootstrap convergence. Homogeneous data: Yi ∈ N (0, Ip ), 0 ≤ i < 6h, data with change
point: Yi ∈ N (0, Ip ), 0 ≤ i < 3h, Yi ∈ N (0.3, Ip ), 3h ≤ i < 6h. The parameters are p = 30,
h ∈ {10, 20, 30, 40, 50}.
The last experimental part presents results of the comparison of the proposed algorithm of
change point detection (LRTOffline) 2 with two other methods: Bayesian online changepoint
6

detection (BOCPD) from Adams and MacKay (2007) and cpt.meanvar(PELT,. . .) (RMeanVar)
from R package. The first method is constructed for online inference, but so far as it returns
CP location with each CP signal, it is also applicable for offline testing scenario. The idea of
this method is predictive filtering: its forecasts a new data point using only the information have
been observed already, where the distribution family is fixed (Normal for the tests in this paper).
Bayesian inference calculates the length of the observed data (from the last CP). The second
algorithm also uses preliminary specified model. Its design focuses into finding multiple changes
in mean and variance in Normally (another distributions also supported) distributed data. The
returned set of change points is the result of sequential testing H0 (existing number of change
points) against H1 (one extra change point) applying the likelihood ratio statistic of the whole
data coupled with the penalty for CP count. RMeanVar performs better than well known method
CUSUM due to synchronous changes in both data parameters mean and variance.
Quality of measurements uses Normalised Mutual Information (NMI). The next equation defines NMI measure of two partitions (X, Y ) of time range by change points
NMI(X, Y ) = 2

H(X) + H(Y ) − H(X, Y )
.
H(X) + H(Y )

1.0

1.0

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

NMI

NMI

H(X) and H(X, Y ) and entropy functions. Higher NMI values (they are in [0, 1]) correspond to
better quality.
Synthetic test data have been generated with different values of the distribution parameter
transition (∆). Each ∆ value corresponds to 10 sampled data sequences over which one compute
measure average. Each data sequence has two, one or none change points. The data has two
distributions: normal (N (θ(1), θ(2))) and Poisson (Po(θ)). Parametric assumption for all methods
is N (θ(1), θ(2)), so Poisson data corresponds to misspecification scenario.

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0

0.0
0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

0.2

0.3

0.4

0.5

0.6

0.7

delta
LRTOffline

BOCPD

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

delta
RMeanVar

LRTOffline

BOCPD

RMeanVar

Figure 4: First data: N (θ(1), θ(2)), second data: Po(θ), data size = 340, parametric assumption
for all methods is N (θ(1), θ(2)), NMI – Normalized Mutual Information between predicted and
reference partitions of time interval with change points, change points count per test {0, 1, 2}.
In the tests with normal data all the methods achieves similar NMI scores. In the tests with
Poisson data (misspecification) RMeanVar has relatively low quality and LRTOffline outperforms
slightly BOCPD method.
Change point detection algorithm (LRTOffline) implementation link:
https://github.com/nazarblch/cpd

References
R. P. Adams and D. J. MacKay. Bayesian online changepoint detection. 2007.

7

G. Biau, K. Bleakley, and D. M. Mason.
Electron. J. Statist., 10(2):2097–2123, 2016.
http://dx.doi.org/10.1214/16-EJS1164.

Long signal change-point detection.
doi:
10.1214/16-EJS1164.
URL

R. Blazek and H. Kim. A novel approach to detection of denial–of–service attacks via adaptive
sequential and batch–sequential change–point detection methods. In D. Kurlander, M. Brown,
and R. Rao, editors, Proc. of IEEE Workshop on Systems, Man, and Cybernetics Information
Assurance, pages 41–50. ACM Press, June 2001.
J. Chen and A. Gupta.
Parametric Statistical Change Point Analysis: With Applications to Genetics, Medicine, and Finance. Springer, 2012. ISBN 0817648003. URL
http://www.amazon.com/Parametric-Statistical-Change-Point-Analysis/dp/0817648003.
V. Chernozhukov, D. Chetverikov, and K. Kato. Comparison and anti-concentration bounds for
maxima of gaussian random vectors. CeMMAP working papers, 2013a.
V. Chernozhukov, D. Chetverikov, and K. Kato. Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors. Ann. Statist., 41(6):2786–2819,
12 2013b. doi: 10.1214/13-AOS1161. URL http://dx.doi.org/10.1214/13-AOS1161.
S. B. Fotopoulos, V. K. Jandhyala, and E. Khapalova. Exact asymptotic distribution of changepoint mle for change in the mean of gaussian sequences. The Annals of Applied Statistics, pages
1081–1104, 2010.
K. Frick, A. Munk, and H. Sieling. Multiscale change point inference. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 76(3):495–580, 2014.
E. Gombay. Sequential change-point detection with likelihood ratios. Statistics & probability
letters, 49(2):195–204, 2000.
P. Haccou, E. Meelis, and S. Van De Geer. The likelihood ratio test for the change point problem
for exponentially distributed random variables. Stochastic processes and their applications, 27:
121–139, 1987.
B. Jandhyala and S. B. Fotopoulos. Capturing the distributional behaviour of the maximum
likelihood estimator of a changepoint. Biometrika, 86(1):129–140, 1999.
H.-J. Kim. Tests for a change-point in linear regression. Lecture Notes-Monograph Series, pages
170–176, 1994.
V. Koltchinskii. A remark on low rank matrix recovery and noncommutative Bernstein
type inequalities, volume Volume 9 of Collections, pages 213–226. Institute of Mathematical Statistics, Beachwood, Ohio, USA, 2013. doi: 10.1214/12-IMSCOLL915. URL
http://dx.doi.org/10.1214/12-IMSCOLL915.
T. L. Lai. Sequential changepoint detection in quality control and dynamical systems. Journal of
the Royal Statistical Society. Series B (Methodological), pages 613–658, 1995.
B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection. Ann. Statist., 28(5):1302–1338, 10 2000. doi: 10.1214/aos/1015957395. URL
http://dx.doi.org/10.1214/aos/1015957395.
Y. Liu, C. Zou, and R. Zhang. Empirical likelihood ratio test for a change-point in linear regression
model. Communications in Statistics—Theory and Methods, 37(16):2551–2563, 2008.
T. Mikosch and C. Starica.
Changes of structure in financial time series and the garch model.
Econometrics 0412003, EconWPA, 2004.
URL
http://ideas.repec.org/p/wpa/wuwpem/0412003.html.
8

A. Polunchenko and A. Tartakovsky. State-of-the-art in sequential change-point detection.
Methodol. Comput. Appl. Probab., 14:649–684, 2011.
R. E. Quandt. Tests of the hypothesis that a linear regression system obeys two separate regimes.
Journal of the American statistical Association, 55(290):324–330, 1960.
A. Shiryaev. Quickest detection problems: Fifty years later. Sequential Anal.: Design Methods
and Applicat., 29:345–385, 2010.
A. N. Shiryaev. On optimum methods in quickest detection problems. Theory of Probability &
Its Applications, 8(1):22–46, 1963.
V. Spokoiny. Multiscale local change point detection with applications to value-at-risk. Ann. of
Stat., 2009.
V. Spokoiny. Penalized maximum likelihood estimation and effective dimension.
arXiv:1205.0498, 2012.

eprint

M. Srivastava and K. J. Worsley. Likelihood ratio tests for a change in the multivariate normal
mean. Journal of the American Statistical Association, 81(393):199–204, 1986.
H. Wang, D. Zhang, and K. G. Shin. Change-point monitoring for the detection of dos attacks.
Dependable and Secure Computing, IEEE Transactions on, 1(4):193–208, 2004.
C. Zou, Y. Liu, P. Qin, and Z. Wang. Empirical likelihood ratio test for the change-point problem.
Statistics & probability letters, 77(4):374–382, 2007.

5

Appendix

5.1.

Likelihood function restrictions

Assume that for each data subset Y = {Yi }t+h
i=t the likelihood function L(θ) = L(θ, Y) has
rather precise approximation by its quadratic Tailor expansion in a local region Θ(r) with central
point θ∗ :
1
L(θ) ≈ L(θ∗ ) + ξ T D(θ − θ∗ ) − kD(θ − θ∗ )k2 ,
2
where
ξ = D−1 ∇L(θ∗ ), D2 (θ) = −∇2 IEL(θ), D = D(θ∗ ).
The following functions characterize the error of quadratic approximation:
1
α(θ, θ0 ) = L(θ) − L(θ0 ) − (θ − θ0 )T ∇L(θ0 ) + kD(θ − θ0 )k2 ,
2
χ(θ, θ0 ) = D−1 ∇α(θ, θ0 ) = D−1 (∇L(θ) − ∇L(θ0 )) + D(θ − θ0 ).
Let with probability 1 − e−x their upper bounds in region Θ(r) satisfy conditions
|α(θ, θ∗ )|
≤ ♦(r, x),
kD(θ − θ∗ )k

kχ(θ, θ∗ )k ≤ ♦(r, x),

where
♦(r, x) = {δ(r) + 6ν0 ωz(x, p)}r,
p
√
12p
z(x, p) = 6p + 2x +
,
g
9

(A)

√
1
r
ω ∼ √ , δ(r) ∼ √ , g ∼ h.
h
h
The stochastic part of the likelihood for independent data has denotation
ζ(θ) = L(θ) − IEL(θ) =

n
X

ζi (θ)

(Zeta)

i=1

The next restriction for the Fisher matrix −D2 (θ) and components ζi (θ, Yi ) deviations ensure
condition (A).
kD−1 D2 (θ)D−1 − Ip k ≤ δ(r),
(dD)
For all kγ1 k = kγ2 k = 1 and |λ| < g and i


λ2 νi2
λ T −1 2
−1
log IEexp
γ1 D ∇ ζi (θ)D γ2 ≤
,
ω
2

X

νi2 = ν02 .

(ED2i)

i

Condition (dD) is responsible for the quadratic approximation of IEL(θ), so
k∇IEL(θ) − ∇IEL(θ∗ ) − D(θ − θ∗ )k ≤ δ(r)r.
In its turn (ED2i) enables linear approximation of ζ(θ).
Lemma 1 (Deviations of empirical process norm). Let condition ED2i is fulfilled, then in the
local region Θ(r) with probability 1 − e−x the next statement holds for all θ, θ0 ∈ Θ(r)
kD−1 (∇ζ(θ) − ∇ζ(θ0 ))k ≤ 6ν0 ωz(x, p)r.
Paper Spokoiny (2012) contains proof for this statement. If the considered point is MLE
b then its concentration in the region Θ(r) follows from condition (L) and Theorem 2.1
(θ = θ)
from Spokoiny (2012).
IEL(θ∗ ) − IEL(θr ) = Ω(r2 ),

5.2.

kD(θ∗ − θr )k = r.

(L)

LRT theorem

Further consider a fixed window position t and window size 2h. We are going to derive explicit
dependence between statistic Th (t) and parameter difference from left and right part of the window
(θr∗ − θl∗ ). Approximation of Th (t) by its quadratic form splits noise and deterministic parts, such
that 2Th (t) ≈ kD(θr∗ − θl∗ ) + ξ lr k2 . In the fixed window position the likelihood function has view
L(θ) = Ll (θ) + Lr (θ) = L(θ, Yl ) + L(θ, Yr ),
D2 = −∇2 IEL(θ∗ ),

Dk2 = −∇2 IELk (θk∗ ),

ξ k = Dk−1 Lk (θk∗ ),

i = {l, r}.

Assume that exist a local region where parameter θ concentrates
Θ(r) = {θ : kD(θ − θ∗ )k < r}.
b θbl , θbr ∈
From condition (A) for function L(θ) it holds with probability 1 − e−x in local region (θ,
Θ(r)) (ref. Theorem 3.2 in Spokoiny (2012) with Th (t) = L(θ, θ0 )) that
p
Dl (θb − θbl )
2Th (t) −
Dr (θb − θbr )

√
≤ 2♦( 2r, x).

b θbl , θbr using Theorem 3.1 from Spokoiny (2012) with notation ξ i (θ) =
Find relation between θ,
Dk−1 ∇Lk (θ)
D(θb − θ) ≤ D−1 {Dl ξ l (θ) + Dr ξ r (θ)} + 2♦(r, x).
10

√
D−1 Dl {ξ l (θ) − Dl (θ − θbl )}
2r, x).
≤
2♦(
D−1 Dr {ξ r (θ) − Dr (θ − θbr )}
Define vector θe that is close to θb

θe = argmin

2

Dl (θ − θbl )

2



+ Dr (θ − θbr )

,

θ

θe = (Dl2 + Dr2 )−1 (Dl2 θbl + Dr2 θbr ).
Dl (θb − θbl )
Dl (θe − θbl )
−
Dr (θb − θbr )
Dr (θe − θbr )

√
e ≤ 2♦(r, x) + 2♦( 2r, x).
≤ D(θb − θ)

Dl (θe − θbl )
= Dlr (θbr − θbl ) ,
Dr (θe − θbr )

Dlr = Dl D−1 Dr .

The temporary result is (with probability 1 − 3e−x )
p

√
≤ 4♦( 2r, x) + 2♦(r, x).

2Th (t) − Dlr (θbr − θbl )

Involve ξ l and ξ r by means of Fisher expansion (equation 3.7 in Spokoiny (2012)) for the model
with two independent components
√
Dr D−1 {Dl (θbl − θl∗ ) − ξ l }
≤
♦(
2r, x).
−1
∗
Dl D {Dr (θbr − θr ) − ξ r }
The final result is Theorem 3, which enables to describe Th (t) function depending on change point
type and subsequently choose appropriate pattern Ph (t) (ref. Section 2).
b θbl , θbr ∈ Θ(r) and the
Theorem 3. Assume that MLE parameters belong to the local region θ,
likelihood has a fit quadratic expansion (A), then with probability 1 − 4e−x for each t
p
√
2Th (t) − kDlr (θr∗ − θl∗ )(t) + ξ lr (t)k ≤ 7♦( 2r, x),
where
ξ lr (t) = Dlr {Dl−2 ∇L(θl∗ , Yl ) + Dr−2 ∇L(θr∗ , Yr )},

5.3.

Dlr = Dl D−1 Dr .

Bootstrap Wilks and Fisher expansions

As it was mentioned in Section 2 the bootstrap procedure allows to yield likelihood function
with two options: each likelihood component is multiplied by weight (weighted bootstrap) or new
data is resampled (empirical bootstrap). The likelihood function in weighted bootstrap case is a
zipped sum with i.i.d weights (u[1 , . . . , u[n ) and independent {li (θ)}ni=1 :
[

L (θ) =

n
X

u[i li (θ),

i=1

ζ [ (θ) = L[ (θ) − L(θ).
Each weight element has Var [ u[i = 1 and IE [ u[i = 1, which is made with a view
IE [ L[ (θ) = L(θ),
[

[

Var ∇L (θ) =

n
X
i=1

11

∇li (θ)∇li (θ)T .

It is expected that Var [ ∇L[ (θ) is close to Var ∇L(θ), which depends on IE∇li (θ) values. For
example in i.i.d models IE∇li (θ∗ ) = 0 and IE∇li (θ) is close to zero in Θ(r).
Here variable ξ has a bootstrap duplicate
X
b =
ξ [ = D−1 ∇ζ [ (θ)
ξ i εi[ , εi[ = u[i − 1.
(xib)
i

Let function α[ (θ, θ0 ) denotes quadratic approximation error for the weighted likelihood function.
Assume further then θ, θ0 ∈ Θ(r).
1
α[ (θ, θ0 ) = L[ (θ) − L[ (θ0 ) − (θ − θ0 )T ∇L[ (θ0 ) + kD(θ − θ0 )k2 .
2
The mean and deviation of the approximation error are (ref. Theorem 3.2 in Spokoiny (2012))
kD−1 ∇IE [ α[ (θ, θ0 )k = kD−1 ∇α(θ, θ0 )k ≤ ♦(r, x),
−1

[

[

[

[

S (θ, θ0 ) = D {∇α (θ, θ0 ) − IE ∇α (θ, θ0 )} =

n
X

D−1 {∇li (θ) − ∇li (θ0 )}εi[ .

(Sb)

i=1

Define function
o

S [ = S [ (θ, θ0 ) − IEY S [ (θ, θ0 )
=

n
X

D−1 {∇ζi (θ) − ∇ζi (θ0 )}εi[ .

i=1
o

The function S [ (taking into account Lemma 1) with probability 1 − e−x fulfils
v
u n
X
uX
o
kS [ k ≤ ♦(r, x)t (νi2 /ν02 )(εi[ )2 ,
νi2 /ν02 = 1.
i=1

i

Let vector ε [ = (ε1[ , . . . , εn[ )T has a restricted exponential moment:
[ γT ε [

log IE e

kγk2
,
≤
2

(Eu)

then by means of the quadratic form deviation (Lemma 16) with variable B 1/2 ε [ where B =
diag(ν12 /ν02 , . . . , νn2 /ν02 ) one get with probability 1 − 2e−x
v
r !
u n
p
uX
√
ν
2x
i
t (ν 2 /ν 2 )(ε [ )2 ≤ tr(B) + 2x max
.
=1+O
0
i
i
i
ν0
n
i=1
Consequently with high probability
o

kS [ k ≤ 2♦(r, x).
In order to find a bound for expectation IEY S(θ, θ0 ) one may employ the deviation bound for
matrix sub-Gaussian sums (ref. Lemma 13). Restrict the norms of IEY S(θ, θ0 ) components and
expand it by Tailor.
[

IEY S (θ, θ0 ) =

n
X

D−1 IEY {∇li (θ) − ∇li (θ0 )}εi[

i=1

=

( n
X

)
D−1 IEY {∇2 li (θ1 )}D−1 εi[

i=1

12

D(θ − θ0 ),

D−1 ∇2 IEli (θ1 )D−1

op

≤ Ci (r),

(dDi)

which leads with probability 1 − 2e−x to
s
X
kIEY S (θ, θ0 )k ≤ 2(x + log p)
Ci2 (r) r.
[

i

Theorem 4 (Weighted bootstrap Wilks). Under conditions (ED2i), (A), (dDi) with sub-Gaussian
[
bootstrap weights (condition Eu) in local region Θ(r) it holds with probability (1−2e−x )(1−4e−x )
(x relates to Y generation and x[ relates to bootstrap weights generation)
|α[ (θ, θ0 )| ≤ ♦[ (r, x)kD(θ − θ0 )k,
where

r
[

♦ (r, x) =

2+O

2x[
n

!!

(Ab)

s
X
♦(r, x) + 2(x[ + log p)
Ci2 (r) r.
i

A modification of Fisher expansion (Theorem 2.2 in Spokoiny (2012)) for the weighted likelihood could be proved using the following property
χ[ (θ, θ0 ) = D−1 (∇L[ (θ) − ∇L[ (θ0 )) + D(θ − θ0 ),
χ[ (θ, θ0 ) = D−1 ∇α[ (θ, θ0 ).
Theorem 5 (Weighted bootstrap Fisher). Under conditions from Theorem 4 it holds with prob[
ability (1 − 2e−x )(1 − 4e−x ) that
b = kD(θ[ − θ)
b − D−1 ∇L[ (θ)k
b ≤ ♦[ (r, x),
kχ[ (θ[ , θ)k
where θ[ , θb are MLE parameters of the weighted and non-weighted likelihood functions.
Theorem 4 enables to prove a theorem similar to Theorem 3 for the bootstrap LRT statistic Th[ .
The proof steps are the same as in Theorem 3.
Theorem 6 (Weighted bootstrap LRT). Assume that MLE parameters belong to the local region
θ[ , θl[ , θr[ ∈ Θ(r), the likelihood has a good quadratic expansion (conditions from Theorem 4), then
[
with probability 1 − 8e−x − 16e−x for each window position t
q
√
2Th[ (t) − Dlr (θbr − θbl )(t) + ξ [lr (t) ≤ 7♦[ ( 2r, x),
where
ξ [lr (t) = Dlr {Dl−2 ∇L[ (θbl , Yl ) + Dr−2 ∇L[ (θbr , Yr )}.
In empirical bootstrap case (Le) with sample size h and dataset size n the function corresponding to quadratic approximation error is


S (θ, θ0 ) =

h
X
i=1

r

n −1
D {∇lk(i) (θ) − ∇lk(i) (θ0 )},
h

where random indexes k(i) ∈ {1, . . . , n} and independent. Define
ui = (0, . . . , 1 , . . . , 0)T ,
k(i)

and
∇l(θ) = (∇l1 (θ), . . . , ∇ln (θ)).
13

Rewrite S  (θ, θ0 ) as
r



S (θ, θ0 ) =
o

o

Define function S  with u =

Ph

h

X
n −1
D {∇l(θ) − ∇l(θ0 )}
(ui − IEui ).
h
i=1

i=1 (ui

− IEui )

o

S  = S  (θ, θ0 ) − IEY S  (θ, θ0 )
r r
r
n −1 2
n
h
o
−1
=
D (∇ ζ(θ1 )u)D
D(θ − θ0 ).
h
h
n
With B = diag(ν12 /ν02 , . . . , νn2 /ν02 ) using Lemma 1 one get
q
o
o
o
o
kS  k ≤ ♦(r, x) (u)T B u ≤ ♦(r, x) kuk∞ ,
o

where kuk∞ is the maximal count of recurrences in bootstrap sample.


h
o
1
e
kuk∞
u11
1−
IEe
≤ hIEe = h
+
≤ h(1 + e).
n
n


Therefore by means of Chernoff bound with probability 1 − e−x
o

kS  k ≤ (log(4h) + x )♦(r, x).
In order to find a bound for expectation IEY S  (θ, θ0 ) one may employ the matrix Bernstein bound
(ref. Lemma 12). Restrict the norms of IEY S  (θ, θ0 ) components and expand it by Tailor.
h r
X
n −1

D IEY {∇lk(i) (θ) − ∇lk(i) (θ0 )}
IEY S (θ, θ0 ) =
h
i=1
( h
)r
Xn
h
D−1 IEY {∇2 lk(i) (θ1 )}D−1
D(θ − θ0 ),
=
h
n
i=1
n −1 2
D ∇ IEli (θ1 )D−1
≤ Ci (r),
h
op

which leads with probability 1 − 2e−x to

q 
2 
[
kIEY S (θ, θ0 )k ≤
Rx + v 5xp r,
3 p

(dDi)

where R = kC(r)k∞ and v = kC(r)k. Summarize bounds for S  (θ, θ0 ).
Theorem 7 (Empirical bootstrap Wilks). Under conditions (ED2i), (A), (dDi) in local region

Θ(r) it holds with probability (1 − e−x )(1 − 3e−x ) (x relates to Y generation and x relates to
bootstrap indexes {k(i)} generation)
|α (θ, θ0 )| ≤ ♦ (r, x)kD(θ − θ0 )k,
where

q 
2

♦ (r, x) = (log(4h) + x ) ♦(r, x) +
kC(r)k∞ xp + kC(r)k 5xp r.
3
Theorem 8 (Empirical bootstrap LRT). Assume that MLE parameters belong to the local region
θ , θl , θr ∈ Θ(r), the likelihood has a good quadratic expansion (conditions from Theorem 7), then

with probability 1 − 4e−x − 12e−x for each t
q
√
2Th (t) − kξ lr (t)k ≤ 7♦ ( 2r, x),






(Ae)

where
b Yl ) + D−2 ∇L (θ,
b Yr )}.
ξ lr (t) = Dlr {Dl−2 ∇L (θ,
r
14

5.4.

Bootstrap measure approximation scheme

Approximation of IP measure and variable f (X) (f is non-random) by corresponded bootstrap
(empirical) measure IP [ and variable f (X [ ) could be done in three steps: 1) approximate f (·) by
smooth function f4 (·) and 1I by smooth indicator; 2) make Gaussian approximation and Gaussian
e < x).
comparison for variables f4 (X) and f4 (X [ ); 3) find anti-concentration bound for IP (f (X)
[
|IP {f (X) < x} − IPP
{f (X [ ) < x}| ≤ ???
X [ = i Xi ε[i
w
e ∈ N (0, Σ )
X −→ X
l
w
e b ∈ N (0, Σ [ )
X [ −→ X

smoothing: 1I[f (X) < x − 4] ≤ f4 (X) ≤ 1I[f (X) < x + 4]
e < x} − IP {f (X)
e < x + 4}| ≤ CA 4
anti-concentration: |IP {f (X)
e ≤
|IEf4 (X) − IEf4 (X)|

Cµ
42

e − IE [ f4 (X
e [ )| ≤
|IEf4 (X)

5.5.

P

i

IEkXi k3

CΣ
kΣ
4

− Σ [ k∞

Anti-concentration

Anti-concentration property can be interpreted as an asymptotic of probability measure depending on event size, tending to zero. Denote by Aε \ A a region of size ε around event A.
Then anti-concentration is higher when probability of Aε \ A is lower. The next lemma deals with
anti-concentration in one dimensional case where the random variable is maximum of a Gaussian
vector.
e ∈ N (m, Σ) ∈ IRp ,
Lemma 2. Let X
p
σ1 ≤ Σii ≤ σ2 ,
p
ei − mi )/ Σii ,
ap = IE max(X
i

then ∀c
ei ∈ [c, c + ε]) ≤ εCA (σ1 , σ2 , ap ),
IP (max X
i


r

σ 
4 σ2
σ2
σ1
σ2
1
+2−
≈ 4 2 ap .
CA (σ1 , σ2 , ap ) =
ap +
−1
2 log
σ1 σ1
σ1
ε
σ2
σ1
Proof may be found by reference Chernozhukov et al. (2013a).
e has a low parameter σ1 then low-variance
Remark. Consider the case when m = 0. If vector X
e could be dropped from consideration. Find the bound σ for such low-variance
components of X
components from condition that with high probability
ei < max X
ei
max X
i: σi <σ

i

p
ei is σ( 2 log(p) + 3), since IE maxi: σ <σ X
ei ≤
The
upper
probabilistic
bound
for
max
X
i:
σ
<σ
i
i
p
σ 2 log(p) and maximum of a Gaussian vector deviated from its expectation such as one dimensional Gaussian variable:


ei > IE max X
ei + r ≤ e−r2 /σ2 .
IP max X
i: σi <σ

i: σi <σ

15

As for extension of this lemma one can employ it for a Composite maximum function of type
maxt Qt (X) in case
max Qt (X) = max{γtT X}, γt ∈ Υ.
t

γt ,t

Apply Lemma 2 with replacement Xi = γtT X, i = (γt , t).
Lemma 3. For parameters
σ1 ≤

q
γtT Σγt ≤ σ2 ,

γt ∈ Υ,

e − m)
e − m)
IE maxt Qt (X
γ T (X
≤
,
ap = IE max tp
t,γt
σ1
γtT Σγt
it holds that
e ∈ [c, c + ε]) ≤ εCA (σ1 , σ2 , ap ).
IP (max Qt (X)
t

5.6.

Gaussian approximation

The goal of Gaussian approximation is to measure the difference between a distribution of
independent random vectors sum and distribution of the closest Gaussian vector. It makes sequential replacement of random arguments in function IEf (X1 +. . .+Xn ) by its normal duplicates
with the same moments:
e1 + . . . + Xn ), . . . , IEf (X
e1 + . . . + X
en ),
IEf (X1 + . . . + Xn ), IEf (X
IEXi = mi ,

Var Xi = Si ,

ei ∈ N (mi , Si ).
X

Each replacement yields difference element erri according to Tailor expansion with t ∈ [0, 1]
ei )
IEf (X + Xi ) − IEf (X + X
o
n

1
1
ei ) X
ei ⊗ X
ei ⊗ X
ei .
= IE tr ∇3 f (X + tXi ) Xi ⊗ Xi ⊗ Xi − IE tr ∇3 f (X + tX
6
6
The first and the second elements of the expansion are substituted. Use property for a tensor T
and any vectors a, b, c
X
Tijk ai bj ck ≤ kT k1 kak∞ kbk∞ kck∞ .
tr{T a ⊗ b ⊗ c} =
i,j,k

Then

1
1
ei )k1 kX
ei k3∞ .
erri ≤ IE ∇3 f (X + tXi ) 1 kXi k3∞ + IEk∇3 f (X + tX
6
6
Consider the case when function f is a composition of a smooth indicator g4 and some other
differentiable function h (for example smooth max). The third gradient of the composition is
000
00
00
0
∇3 (g4 ◦ h) = g4
∇h ⊗ ∇h ⊗ ∇h + 2g4
∇2 h ⊗ ∇h + g4
∇h ⊗ ∇2 h + g4
∇3 h,

∇3 (g4 ◦ h)

1

000
00
| ∇2 h
≤ |g4
| k∇hk31 + 3|g4

1

0
k∇hk1 + |g4
| ∇3 h

1

.

0
Assume that g4 grows from 0 to 1 in interval [z, z + 4]. So g4
= 0 outsize [z, z + 4]. Furthermore
in this case
∇3 f = ∇3 (g4 ◦ h) = ∇3 (g4 ◦ h) 1I[z ≤ h(x) ≤ z + 4].

Move to approximation of a distribution function. The aim is to find upper bound for
!
!
!
!
X
X
ei ≤ z .
EG = IP H
Xi ≤ z − IP H
X
i

i

16

Let h be a smooth approximation for function H such that
H(x) − 4 ≤ h(x) ≤ H(x) + 4,
1I[H(x) ≥ z + 4] ≤ g4 (h(x)) ≤ 1I[H(x) ≥ z − 4].
P
Since IE 1I[H(X) ≤ z] = IP (H(X) ≤ z), one get with X = i Xi

(SmAprx)

IP (H(X) ≤ z − 4)
≤ IEf (X)
e +
≤ IEf (X)

X

erri

i

e ≤ z + 4) +
≤ IP (H(X)

X

erri .

i

Subsequently
e ≤ z ± 24) ≤
IP (H(X) ≤ z) − IP (H(X)

X

erri .

i

e (ref. Section 5.5.):
Suppose an anti-concentration property for random variable H(X)
e ∈ [z, z + 4]) ≤ CA 4.
IP (H(X)
Then
EG ≤ 2CA 4 +

X

erri

i

and
EX (1I[z ≤ h(X + tXi ) ≤ z + 4])
≤ EX (1I[z ≤ H(X + tXi ) ≤ z + 24])


e
≤ EXe 1I[z ≤ H(X + tXi ) ≤ z + 24] + EG
≤ CA (24 + kXi k∞ ) + EG .
Demote the sum of the third moments
µ3X

X
ei
=
IE kXi k3∞ + IE X
i

3
∞


.

Assume restriction for the third derivatives and define constant Cµ :
1
∇3 (g4 ◦ h)
6
Then
X
i

1

≤

1
Cµ .
43

!
X
1
erri ≤ 3 Cµ (24CA + EG ) µ3X + O
IE kXi k4∞ .
4
i

Furthermore, neglecting the sum of the fourth moments
EG .

1
Cµ (24CA + EG )µ3X + 24CA ,
43

e ≤ z) ≤ 5C 1/3 CA µX .
IP (H(X) ≤ z) − IP (H(X)
µ
17

(ErrG)

Consider linear forms (AX) with sparse-row matrices and random vector X = (X1 , . . . , Xn )
with independent elements (sub-vectors). In order to make approximation by Gaussian vector
e one have to group (X1 , . . . , Xn ) such that elements in one group has no common non-zero
AX
coefficients in each row of matrix A. This allows following representation
AX = Z = Z1 + . . . + Zh ,

Zi = AFi X,

vectors {Zi } are independent and Fi is a filter for the i-th group that sets matrix columns to
0 related to the other groups. In case each row of the matrix A has no more than h non-zero
elements then minimal groups count equals to h. The next statement confirms it.
Lemma 4. Let each element in a set of subsets {Ms } has size h. Then subsets {Z1 , . . . , Zh } exist
with properties
[
[
\
\
Zs =
Ms , Z k
Ms = 1, Zk
Zs = 0.
s

s

S
Proof. Build subsets {Z1 , . . . , Zh } constructively. Take one element from s Ms . Exclude this
element from all {Ms } and add
Z1 . Mark subsets which contain this element as {M0 s }.
S it to S
add it to Z1 . Repeat this procedure until
Take
from s Ms \ s M0 s and S
S another
S element
0
| s Ms \ s M s | = 0. Then do the same steps for s Ms \ Z1 and obtain Z2 with the required
properties by construction.
Summarise results of this Section. Apply equation (ErrG) to Zi instead of Xi .
Lemma 5. Let matrix A has at most h non-zero elements in each row and non-zero elements
correspond to independent elements in vector X = (X1 , . . . , Xn ). Then Gaussian approximation
e has following upper bound
with vector X


e ≤ z ≤ 5Cµ1/3 CA µZ ,
IP (H(AX) ≤ z) − IP H(AX)
(ErrGA)
where
µ3Z = IE kXk3∞

h
X

|||AFi |||31 ≤ IE kXk3∞ hkAk3∞ ,

i=1

and Cµ , CA correspond to
1
1
∇3 (g4 ◦ h(AX)) 1 ≤ 3 Cµ ,
6
4
e ∈ [z, z + 4]) ≤ CA 4.
IP (H(AX)

5.7.

Gaussian comparison

This Section deals with extension of the Gaussian approximation ( ref. Section5.6.) for the
e are slightly different. Let as previously
case P
when the second moments of Var X and Var X
n
X = i=1 Xi and
1
ei
Var Xi − Var X
∼ .
n
∞
e
Then after sequential replacement Xi → Xi the approximation bound term will be
n
X

erri ∼ 1,

i=1

while

1
∼√ .
n
∞
The next lemma resolves this problem. At first one is able to make Gaussian approximation with
e and then compare two Gaussian vectors with different variances.
equal variances (Var X = Var X),
e
Var X − Var X

18

Lemma 6. Let X and Y be two zero mean Gaussian vectors with ΣX = Var(X) and ΣY =
Var(Y ) . Let also f (X) be a smooth function. Then
 = IEf (X) − IEf (Y ) ≤

1
kΣX − ΣY k∞ kIE∇2 f k1 ,
2

(EErrGC)

where k∇2 f k1 = supX k∇2 f (X)k1 .
Proof. Without loss of generality assume that X and Y are given on the same probability space
and independent. For each t ∈ [0, 1] , define
√
√
Z(t) = t X + 1 − t Y,
√
√

Ψ (t) = IEf (Z(t)) = IEf t X + 1 − t Y .
Z1
 = |Ψ (1) − Ψ (0)| =

Ψ 0 (t)dt .

0

1 
Ψ 0 (t) = IE[∇f (Z(t))> Z 0 (t)] = IE t−1/2 X − (1 − t)−1/2 Y
2

>


∇f (Z(t)) .

To compute this expectation, we apply the Stein identity. Let W be a zero mean Gaussian vector.
Then for any C 1 vector function s it holds
IE[W s(W )> ] = Var(W ) IE[∇s(W )].




IE ∇f (Z(t))X > = t1/2 ΣX IE ∇2 f Z(t)




IE ∇f (Z(t))Y > = (1 − t)1/2 ΣY IE ∇2 f Z(t) ,
 

1 
tr ΣX − ΣY IE ∇2 f Z(t)
2


1
≤ kΣX − ΣY k∞ IE ∇2 f Z(t)
2

Ψ 0 (t) ≤

1

≤

1
kΣX − ΣY k∞ IEk∇2 f k1 .
2

Find an upper bound for distribution difference
EGC = |IP (H(X) ≤ z) − IP (H(Y ) ≤ z)|
Identically to Section 5.6. consider a composition of a smooth indicator g4 and a differentiable
smooth approximation h for function H. The approximation satisfies condition (SmAprx). As0
sume that g4 grows from 0 to 1 in interval [z, z + 4]. So g4
= 0 outsize [z, z + 4]. Assume
following restriction for the second derivative of the function f = g4 ◦ h and define CΣ :
1 2
1
k∇ g4 ◦ hk1 ≤ 2 CΣ .
2
4
Then by means of Lemma 6 and taking into account IP = IE 1I
|IP (H(X) ≤ z) − IP (H(Y ) ≤ z ± 24)|
≤ IEf (X) − IEf (Y )


2
≤ k∇ g4 ◦ hk1 IE 1I h(X) ∈ [z, z + 4] kΣX − ΣY k∞ ,


1
≤ 2 CΣ IE 1I H(X) ∈ [z, z + 24] kΣX − ΣY k∞ .
4
19

Following the logic from Section 5.6. the anti-concentration (ref. Section 5.5.) property allows
24-shift elimination and provides an upper bound for IE 1I[H(X) ∈ [z, z + 24]].
IE 1I[H(X) ∈ z ± 24] ≤ 24CA ,
EGC ≤

2
CΣ 4CA kΣX − ΣY k∞ + 24CA ,
42

Optimize over 4 values
1/2

1/2
EGC ≤ 4CΣ CA kΣX − ΣY k∞
.

Summarise the discussion related to EGC bound.
Lemma 7. Let function H(X) with Gaussian argument satisfies anti-concentration with CA
constant, X and Y be two zero mean Gaussian vectors with ΣX = Var(X) and ΣY = Var(Y ) .
Then under condition (SmAprx)
1/2

1/2
|IP (H(X) ≤ z) − IP (H(Y ) ≤ z)| ≤ 4CΣ CA kΣX − ΣY k∞
.

(ErrGC)

Consider a random vector X = (X1 , . . . , Xn ) with independent sub-vectors {Xi } and matrix
A which consists of blocks of size [p × p], p = dim(Xi ). Define also a bootstrap vector with
independent random weights {εi[ }, such that X [ = (ε1[ X1 , . . . , εn[ Xn ).


a11 Ip . . . a1n Ip
.
...
A=
an1 Ip . . . ann Ip
The next comparison of H(AX) and H(AX [ ) my means of (ErrGC) require upper bound for the
maximal element of the covariance matrices difference.
[
kΣAX − ΣAX
k∞ =

max

aT ,bT ∈rowsA

|aT (Σ − Σ [ )b|.

Let for a fixed rows a, b with probability 1 − e−x
|aT (Σ − Σ [ )b| =

X

ai bj Σij −

ij

X

ai bj Xi Xj ≤ (x).

ij

Note that elements in sum (ai bi Xi Xj ) are independent due to the specific block structure of matrix
A. Then the joint bound with probability 1 − e−x is
[
kΣAX − ΣAX
k∞ ≤ (x + 2 log(np)).

Involve the upper bound for covariance matrix deviations (ErrVD) with εi = Xi and Ui = ai /V


√
2
2
2
2
(x) = V
Rεε x + 2vεε 5x + δ kbk ,
3
p
where vεε = δ kΣk∞ (3 + kbk),
X
X
V δ = kAk∞ , V 2 =
ai bj Σij ≤ |||A|||2 kΣk∞ , kbk2 =
IEXi2 .
ij

i: ai >0

√
Finally under assumption 23 Rεε x < vεε 5x with probability 1 − 1/n
p
[
kΣAX − ΣAX
k∞ ≤ 10 log(np)|||A|||kAk∞ kΣk∞ (3 + kbk) + kAk2∞ kbk2 .
20

(dSA)

5.8.

Smooth-max properties

Denote by h(x) a smooth maximum function which converges to maxi xi when β → ∞.
X
h(x) = β −1 log u(x), u(x) =
eβxi , x = (x1 , x2 , . . .).
i

Characterize its derivatives which will be used further in Tailor expansions.
Lemma 8. All derivatives of h(x) of order m = {1, 2, 3} have following upper bounds ∀x
k∇(m) h(x)k1 ≤ β m−1 .
Proof.
∇u
,
u


∇ ⊗ ∇u ∇u ⊗ ∇u
2
−1
∇ h(x) = β
−
,
u
u2


∇ ⊗ ∇ ⊗ ∇u ∇ ⊗ ∇u ⊗ ∇u ∇ ⊗ (∇u ⊗ ∇u)
∇u ⊗ ∇u ⊗ ∇u
3
−1
−
∇ h(x) = β
−
+2
.
u
u2
u2
u3
P
(i) that satisfies condition i pi = 1. The first tensor norm equals to the convoluDefine pi = ∇u
u
tion maximum with vectors α, φ, γ under restriction kαk∞ = 1, kφk∞ = 1, kγk∞ = 1.
X

X
X
αT ∇2 h(x)γ = β
pi αi γi −
pi α i
pj γj
∇h(x) = β −1

oo

= β (IEαγ − IEαIEβγ) = βIE αγ ≤ βkαk∞ kγk∞ ,
X

∇3ijk h(x)αi φj γk = β 2 (IEαφγ − IEαIEφγ − IEαφIEγ − IEαγIEφ + 2IEαIEφIEγ)

ijk



ooo
= β IE αφγ ≤ β 2 kαk∞ kφk∞ kγk∞ .
2

Taking maximum provides the required restriction for L1 tensor norms.
The next property of h(x) with x ∈ IRp characterise the error of smooth maximum approximation
max(xi ) ≤ h(x) ≤ max(xi ) + β −1 log(p).
i

i

Applying indicator for both parts of this inequality yields statement similar to SmAprx which is
used in Gaussian approximation (ref. Section 5.6.) and Gaussian comparison (ref. Section 5.7.).
Lemma 9. For a smooth indicator function g4 (g4 grows from 0 to 1 inside interval [z, z + 4])
it holds with 4 = β −1 log(p) that




1I max xi > z + 4 ≤ g4 h(x) ≤ 1I max xi > z − 4 .
0≤i≤p

0≤i≤p

This result allows to use inequalities (ErrG) and (ErrGC) with following restrictions for the
derivatives. Let g4 (x) = g(x/4) and H(x) = max(x) then
Cµ =


max(g 0 , g 00 , g 000 )
1 + 3 log p + log2 p ,
6
CΣ =

max(g 0 , g 00 )
(1 + log p) .
2
21

5.9.

Composite smooth-max properties

Consider a composite maximum function with its smooth approximation h(q(X))
max Qt (x) ≈ h(q(x)),

q(x) = (q1 (x), . . . , qT (x)).

1≤t≤T

where h in the smooth approximation for maximum from Section 5.8. and qt is a smooth approximation for Qt with property:
∀t : Qt (x) ≤ qt (x) ≤ Qt (x) + βq−1 log(p).
Combination with maximum approximation property h(x) ≤ maxt (xt ) + β −1 log(p) leads to statement
max(Qt (x)) ≤ h(q(x)) ≤ max(Qt (x)) + (βq−1 + β −1 ) log(p).
t

t

This allows to extend Lemma 9.
Lemma 10. For a smooth indicator function g4 (g4 grows from 0 to 1 inside interval [z, z + 4])
it holds with 4 = (βq−1 + β −1 ) log(p) for all x and z that




1I max Qt (x) > z + 4 ≤ g4 (h(q(x)) ≤ 1I max Qt (x) > z − 4 .
1≤t≤T

1≤t≤T

Assume also restriction for derivatives of the functions qt :
∀t : k∇(m) qt (X)k1 ≤

Cqm−1
,
4m−1

from which follows
k∇h(q(X))k1 ≤ k∇hk1 k∇qt k1 ≤ 1,
k∇2 h(q(X))k1 ≤ k∇2 hk1 k∇qt k21 + k∇hk1 k∇2 qt k1 ≤

log T
Cq
+
,
4
4

k∇3 h(q(X))k1 ≤ k∇3 hk1 k∇qt k31 + 3k∇2 hk1 k∇2 qt k1 k∇qt k1 + k∇hk1 k∇3 qt k1
≤

1
(log2 T + 3Cq log T + Cq2 ).
42

So one can override the constants used in (ErrG) and (ErrGC) for the case H(X) = max1≤t≤T Qt (X):
Cµ =


max(g 0 , g 00 , g 000 )
1 + 3 log T + 3Cq + log2 T + 3Cq log T + Cq2 ,
6
CΣ =

max(g 0 , g 00 )
(1 + log T + Cq ) ,
2

where g(x/4) = g4 (x).

5.10.

LRT bootstrap precision

Consider statistic TPh (τ ) and the corresponded bootstrap analogy TP[h (τ ) defined in Section 2.
One could also repeat the further logic for empirical bootstrap TPh (τ ). Describe the bootstrap
approximation for the quadratic form of the statistic TPh (τ ) on the grounds of Theorems 3 and 6.
The quadratic form of TPh (τ ) is
(
)
X
max
Pτ (t)kξ lr (t)k .
(TPQ)
1≤τ ≤n

t

22

The corresponded bootstrap quadratic form is
)
(
X
[
max
Pτ (t)kξ lr (t)k .
1≤τ ≤n

t

Our aim is to show that these two forms are close by distribution. For simplification assume that
for all window positions the true model parameter is fixed θl∗ = θr∗ = θ∗ . Then ξ lr (t) doesn’t pay
attention on parameter changes and
ξ lr (t) =

t+h
X

ξi −

i=t

t+2h
X

ξi ,

i=t+h

where ξ i = D−1 ∇li (θ∗ ) and D−1 = Dlr Dl−2 = Dlr Dr−2 . Use smooth-max approximation for
the composite maximum function (ref. Section 5.9.) with argument ξ lr = Z = Aξ, where ξ =
(ξ 1 , . . . , ξ n ).
(
)
X
1
Pτ (t)kZ(t)k ≈ h(q(Z)),
√ max
p τ
t
where
qτ (Z) =

X
t

kWt Zk
Pτ (t) √ ,
p

Wt = diag(0, . . . , 1tp , . . . , 1tp+p , . . . , 0).

By means of Lemmas 5, 7 and 3 we are to prove the next Theorem.
Theorem 9 (Bootstrap Approx). Let window size be equal to h, the dataset size P
– n, the model
dimension – p, pattern functions Pτ (t) be independent from kξ lr (t)k and normalized t |Pτ (t)| = 1.
Then the distribution difference between TPh (τ ) and TP[h (τ ) is
!
!
X
X
IP max
Pτ (t)kξ lr (t)k > z − IP [ max
Pτ (t)kξ [lr (t)k > z
1≤τ ≤n

1≤τ ≤n

t

t

1/2

[
≤ 5Cµ1/3 CA µZ + 4CΣ CA kΣAξ − ΣAξ
k1/2
∞ ,

where the constants Cµ and CΣ correspond to composite maximum function (ref. Section 5.9.)
and
p
√
[
kΣAξ − ΣAξ
k∞ ≤ 10 log(np) 2hkΣξ k∞ (3 + kbk) + kbk2 ,
2

kbk = max
t

t+2h
X

kIEξ i k2∞ ,

i=t

µ3Z ≤ IE kξk3∞ hkAk3∞ = IE kξk3∞ h.
Proof. One should estimate distribution difference from replacement of the random argument in
[
[
T
e
e
statistic (TPQ):
P ξ → ξ → ξ . Note that ξ ∈ N (IEξ, Var(ξ)) and ξ ∈ N (0, diag(ξ i ξ i )). Taking
into account t |Pτ (t)| = 1 estimate q’s derivatives required for (ErrG) and (ErrGC).
k∇qτ (Z)k1 =

X
t

2

k∇ qτ (Z)k1 =

X
t

k∇3 qτ (Z)k1 ≤

X
t

kWt Zk1
|Pτ (t)| √
≤ 1,
pkWt Zk

√
2 p
kWt − ZWt Z T /kWt Zk2 k1
|Pτ (t)|
≤
,
√
pkWt Zk
kWt Zk

kWt ⊗ Wt Zk1 X
kZWt Z T ⊗ Wt Zk1
6p
3|Pτ (t)| √
+
3|P
(t)|
≤
.
√
τ
3
5
pkWt Zk
pkWt Zk
kWt Zk2
t
23

Then the constant Cq figuring in Section 5.9. has bound
√
6p4
Cq ≤
. 1,
kWt Zk
and therefore with T = n one get

max(g 0 , g 00 , g 000 )
Cµ =
5 + 6 log n + log2 n ,
6
max(g 0 , g 00 )
(2 + log n) .
2
Finally the bootstrap approximation error for statistic (TPQ) is
!
!
X
X
IP max
Pτ (t)kξ lr (t)k > z − IP [ max
Pτ (t)kξ lr[ (t)k > z
CΣ =

1≤τ ≤n

1≤τ ≤n

t

t

!
≤ IP

max

1≤τ ≤n

X

!
− IP [

Pτ (t)kξ lr (t)k > z

t

max

1≤τ ≤n

X

Pτ (t)ke
ξ lr (t)k > z

t

!

!
+ IP

max

1≤τ ≤n

X

− IP [

Pτ (t)ke
ξ lr (t)k > z

t

max

1≤τ ≤n

X

Pτ (t)kξ lr[ (t)k > z

t

1/2

[
1/2
≤ 5Cµ1/3 CA µZ + 4CΣ CA kΣAξ − ΣAξ
k∞
,

From (dSA) follows
p
[
kΣAξ − ΣAξ
k∞ ≤ 10 log(np)|||A|||kAk∞ kΣk∞ (3 + kbk) + kAk2∞ kbk2
p
√
= 10 log(np) 2hkΣξ k∞ (3 + kbk) + kbk2 .
e = Ae
For parameter CA from Lemma 3 one have to estimate σ1 , σ2 and ap for the case X
ξ
X
pσ12 ≤
Pτ (t)Pτ (s)γtT Wt ΣWs γs ≤ pσ22 ,
t,s

γ = (γτ , . . . , γτ +h ),

e
γt ∈ argmax γtT Wt X.
kγt k=1

The domain of γ can be restricted with high probability by condition
pσ12 ≤

X
t,s

Pτ (t)Pτ (s)

e
eT
X
X
Wt ΣWs
≤ pσ22 ,
e
e
kWt Xk
kWs Xk

P
e
e
Define random variables w12 , w22 with the same distribution as kWt XkkW
s Xk and Pτ =
t Pτ (t)Wt ,
then
e T Pτ ΣPτ X
e ≤ pσ 2 w2 .
pσ12 w12 ≤ X
2 2
In the middle is a chi-square random variable with degrees of freedom dependent from matrix
B = Σ 1/2 Pτ ΣPτ Σ 1/2 . With regarding chi-square deviations (Theorem 17) it holds with high
probability that
p
e T Pτ ΣPτ X
e ≥ tr(B) − tr(B 2 ) log n,
X
and correspondingly one may set bounds σ1 , σ2 as follows
p
tr(B) − tr(B 2 ) log n
2
σ1 = min
, σ22 = max kWt ΣWt kop /p.
2
t
t
p z (Wt ΣWt , log n)
24

P
For parameter ap one have to find bound for IE maxt Qt (Ae
ξ). From conditions t |Pτ (t)| = 1 and
kWt Zk ≤ z(Wt ΣWt , x) (Theorem 16) follows that
IE max Qt (Ae
ξ) ≤ z(Wt ΣWt , log n),
t

and furthermore

√ 2
pz (Wt ΣWt , log n)
IE maxt Qt (Ae
ξ)
p
≤
.
ap ≤
σ1
(tr(B) − tr(B 2 ) log n)1/2

The final estimation for CA is
CA ≈ 4

5.11.

pz 4 (Wt ΣWt , log n)
σ2
p
a
≤
max
kW
ΣW
k
.
p
t
t
op
t
σ12
(tr(B) − tr(B 2 ) log n)3/2

Generalized linear model case

Generalized linear models (GLM) are frequently used for modeling the data with special structure: categorical data, binary data, Poisson and exponential data, volatility models, etc. All these
examples can be treated in a unified way by a GLM approach. This section specifies the results
and conditions to this case. Let Y = (Y1 , . . . , Yn ) ∼ IP be a sample of independent r.v.’s. The
parametric GLM is given by Yi ∼ IPψiT θ , where ψi are given factors in IRp . Generalised linear
model may be presented in form
L(θ) = S T θ − A(θ)
where
A(θ) =

n
X

g(ψiT θ),

S=

i=1

n
X

Yi ψi .

i=1

This model has following properties (ref. Section 5.1.)
2

2

−∇ IEL(θ) = D (θ) =

n
X

g 00 (ψiT θ)ψi ψiT ,

i=1

δ(r) = ag δψ r,
where
ag = max
x

g 000 (x)
,
g 00 (x)

δψ = max D−1 ψi .
i

As g(·) is convex, it holds g 00 (u) ≥ 0 for any u and thus D2 (θ) ≥ 0 . An important feature of a
GLM is that the stochastic component ζ(θ) of L(θ) is linear in θ : with εi = Yi − IEYi
ζ(θ) =

n
X

ξ = D−1

εi ψi> θ,

i=1

n
X

εi ψi .

i=1

Linearity in θ of the stochastic component ζ(θ) and concavity of the deterministic part IEL(θ)
allow for a simple and straightforward proof of the result about localisation of the MLE θb in the
region Θ(r0 ) (ref. condition L and Theorem 2.1 in Spokoiny (2012)).
Theorem 10. If for some r0 > 0 , D(θ) fulfils dD condition with δ(r0 ) < 1 , and if ξ is a
sub-Gaussian vector (SG) with variance matrix B then

IP θb 6∈ Θ(r0 ) ≤ 2e−x
where lover bound for r0 is defined by z(B, x) from Lemma 16
r0 {1 − δ(r0 )} ≥ 2z(B, x).
25

Proof. The function L(θ) is concave in θ because
−∇2 L(θ) = D(θ) ≥ 0.
Therefore, it suffices to check that for each θ with kD(θ−θ∗ )k = r0 that with probability 1−2e−x
L(θ∗ ) − L(θ) > 0.
1
2
D(θ◦ )(θ − θ∗ )
2
1 − δ(r0 )
kD0 (θ − θ∗ )k2
≥ (S − IES)> (θ − θ∗ ) +
2
1 − δ(r0 ) 2
1 − δ(r0 ) 2
= ξ > D0 (θ − θ∗ ) +
r0 ≥ −kξk r0 +
r0 .
2
2

L(θ∗ ) − L(θ) = (θ − θ∗ )> ∇L(θ∗ ) +


If kξk ≤ r0 1 − δ(r0 ) /2 , then this implies L(θ∗ ) − L(θ) > 0 , and the result follows. Theorem
16 gives a probabilistic restriction for kξk:


IP kξk > z(B, x) ≤ 2e−x .

b
As a corollary, we obtain Fisher and Wilks expansions (condition A holds for the case θ = θ)
for generalized linear models with
♦(r0 , x) = δ(r0 )r0 = ag δψ r20 .
Theorem 11. Suppose the conditions of Theorem 10 for some r0 . Then with probability 1−2e−x

D θb − θ∗ − ξ ≤ r0 δ(r0 ),
b θ∗ ) − kξk2 ≤ 2r2 δ(r0 ) + r2 δ 2 (r0 ).
2L(θ,
0
0
q
b θ∗ ) − kξk ≤ 3r0 δ(r0 ).
2L(θ,
Proof.
α(θ, θ∗ ) = A(θ) − A(θ∗ ) − (θ − θ∗ )> ∇A(θ∗ ) − kD(θ − θ∗ )k2 /2
=


1
(θ − θ∗ )> D(θ∗ ) − D(θ◦ ) (θ − θ∗ ) ,
2

where θ◦ is a point on the interval between θ and θ∗ . Condition (dD) implies
α(θ, θ∗ ) ≤

δ(r0 )
D0 (θ − θ∗ )
2

2

≤

δ(r0 ) 2
r.
2 0

Assume that for independent εi = Yi − IEYi
 1 2 2
log IE exp λ s−1
i = 1, . . . , n,
i ε i ≤ ν0 λ ,
2
n
X
2
V0 =
s2i ψi ψi> .

|λ| ≤ g1 .

(SGeps)

i=1
2

The squared norm kξk is a quadratic form of εi and one can employ upper bound from
Section 5.14..
26

√
√
p + 2x according to Lemma 15, fix

Theorem 12. Suppose (SGeps), z(p, x) ≤

r0 = 4ν0 z(p, x),
Then the conditions of Theorem 10 are fulfilled with δ(r0 ) ≤ ag (r0 ) δΨ r0 and

ν2
log IE exp γ > ξ ≤ 0 kDV0−1 γk2 ,
2

kγk ≤ g =

g1
.
δΨ si

Proof. For each θ ∈ Θ(r0 ) and i ≤ n , it holds
D−1 ψi

ψi> θ − ψi> θ∗ =

∗

D(θ) − D(θ ) =

>

D(θ − θ∗ ) ≤ kD−1 ψi k r0 ≤ δΨ r0 .

n
X


g 00 (ψi> θ) − g 00 (ψi> θ∗ ) ψi ψi> .

i=1

g 00 (ψi> θ) − g 00 (ψi> θ∗ ) =
max
i≤n


g 000 (ψi> θ◦ ) >
ψi θ − ψi> θ g 00 (ψi> θ∗ ).
> ∗
00
g (ψi θ )


g 000 (ψi> θ◦ ) >
ψi θ − ψi> θ ≤ ag (r0 ) δΨ r0
> ∗
00
g (ψi θ )

The last statement of theorem follows from
n
X
 > −1

log IE exp γ V0 (S − IES) =
log IE exp λi s−1
i εi ,
i=1

|λi | = |γ > V0−1 ψi | si ≤ g kV0−1 ψi k si ≤ g1 .
n
n
 > −1

ν02 X 2 ν02 X > −1
ν2
λi =
γ V0 ψi ψi> s2i V0−1 γ = 0 kγk2 .
log IE exp γ V0 (S − IES) ≤
2 i=1
2 i=1
2

Finally, condition SGeps ensure all the required likelihood restrictions (condition A and θb ∈
Θ(r)) from Section 5.1..
5.11..1

Bootstrap Wilks and Fisher expansions

This Section projects statements from Section 5.3. with overriding some aspects for GLM. In
this case function Sb doesn’t depend on Y

[

[

S (θ, θ0 ) = IEY S (θ, θ0 ) =

n
X

D−1 {∇li (θ) − ∇li (θ0 )}ε [

i=1

=

n
X

D−1 {∇g(ψiT θ) − ∇g(ψiT θ0 )}ε [ .

i=1

The gradients difference has view with θ1 ∈ (θ, θ0 )
∇g(ψiT θ) − ∇g(ψiT θ0 ) = g 00 (ψiT θ1 )ψi ψiT D−1 D(θ − θ0 ).
Then S [ (θ, θ0 ) may be presented as a sum of independent matrices
!
n
X
S [ (θ, θ0 ) =
Si D(θ − θ0 ),
i=1

27

where
Si =

di dTi εi[ ,

q
di = g 00 (ψiT θ1 )D−1 ψi ,

kdi k2 ≤ (g 00 (ψiT θ∗ ) + δ(r0 )) δψ2 = δd2 .
From deviation bound for matrix Gaussian sums (Lemma 14) with probability 1 − 2e−x
n
X

√
≤ δd2 ν0 ksk 2x.

Si

i=1

op

Theorem 13. Theorems 5, 4 and 6 have in GLM case following upper bound
√
♦[ (r0 , x) = r0 δd2 ν0 ksk 2x + 2♦(r0 ).
b Conditionally on Y from Theorem 12
An important variable in Theorem 13 is ξ [ = D−1 L[ (θ).
n
n


1 X 2 1 X > −1
log IE [ exp γ > V0−1 (S [ − IES [ ) ≤
λi =
γ V0 ψi ψi> ε2i V0−1 γ
2 i=1
2 i=1

By definition of V0
IE

( n
X

)
V0−1 ψi ψi> εi V0−1

2

=I

i=1

Applying consequence of Theorem 12 it holds with probability 1 − e−x for ηi = ε2i − IEε2i
n
X


V0−1 ψi ψi> ηi V0−1

i=1

op

where

s
vψ ≤ δψ max
i

Theorem 14. Assume (2/3)Rψ xp + vψ

Eηi2
,
s2i

Rψ . δψ2 ksk∞ ν0 .

p
5xp ≤ 1 then with probability 1 − e−x


ν2
log IE [ exp γ > ξ [ ≤ [ kγk2 ,
2
5.11..2

p
2
≤ Rψ xp + vψ 5xp ,
3

ν[ =

√

2 DV0−1

op

,

kγk ≤ g.

Bootstrap likelihood argmax concentration

An important property of GLM likelihood function is convexity: −∇2 L = −∇2 IEL ≥ 0. This
property is useful for MLE concentration proof (Theorem 10). Bootstrap likelihood is convex with
high probability under an additional condition described in the following statement.
Theorem 15. Assume that ∀i
q
√
g 00 (ψiT θ) D−1 (θ)ψi
2x < 1,
then with probability 1 − e−x
−∇2 L[ (θ) > 0.
Proof.
−∇2 L[ (θ) =

n
X

!
g 00 (ψiT θ)ψi ψiT ui = D(θ) I +

i=1

X
i

28

di dTi ε[i

D(θ).

!
−∇2 L[ (θ) ≥ 0 ⇔ 1 − λmax −

X

di dTi ε[i

≥ 0,

i

p
P
T
where di = g 00 (ψiT θ)D−1 (θ)ψi ,
i di di = I. Use Lemma 14 to get matrix deviation bound
which states that with probability 1 − 2e−x
!
X
√
λmax −
di dTi ε[i ≤ max kdi k 2x.
i

i

Using previous result one is able to prove θb[ concentration analogically to Theorem 10.
Theorem 16. Under condition from lemma 15


b > r[ ≤ 4e−x ,
IP D(θb[ − θ)

r[0 = 8ν0[ z(p, x).

0

b on the border.
Proof. Detect sufficiently large radius r of local region Θ(r) with L[ (θ) < L[ (θ)
[
[
With ζ = L − L:
b = ζ [ (θ) − ζ [ (θ)
b + L(θ) − L(θ).
b
L[ (θ) − L[ (θ)
b ≤ − 1 − δ(2r) r2 .
L(θ) − L(θ)
2
Use variable Sb, from its definition
b − ∇ζ [ (θ)
b T (θ − θ)
b = S [ (θ, θ)
b T D(θ − θ).
b
ζ [ (θ) − ζ [ (θ)
According to Theorem 13

b ≤ ♦[ (r, x). Gradient ζ [ (θ)
b may be estimated by means of
S [ (θ, θ)

sub exponential property 14 and bound for ξ [ provided by Theorem 15.
b = ξ[
D−1 ∇ζ [ (θ)
and
IP

[



ξ

[

>

ν0[ z(p, x)



≤ 2e−x .

Finally, the proving statement requires condition (assume 1 − δ(2r) > 1/2 and ♦[ (r, x)/r < 1)
1 2
b + rν [ z(p, x) ≤ 2rν [ z(p, x).
r > r2 S [ (θ, θ)
0
0
4

5.12.

Matrix Bernstein inequality

Lemma 11 (MasterPbound). Assume that S1 , . . . , Sn are independent Hermitian matrices of the
same size and Z = ni=1 Si . Then
!
n
X
1
IEλmax (Z) ≤ inf log tr exp
log IEeθSi ,
(1)
θ>0 θ
i=1
!
n
X
IP {λmax (Z) ≥ z} ≤ inf e−θz tr exp
log IEeθSi .
(2)
θ>0

i=1

29

Proof. By the Markov inequality
IP {λmax (Z) ≥ z} ≤ inf e−θz IE exp(θλmax (Z)).
θ

Recall the spectral mapping theorem: for any function f : IR → IR and Hermitian matrix A
eigenvalues of f (A) are equal to eigenvalues of A . Thus

exp(θλmax (Z)) = exp(λmax (θZ)) = λmax exp(θZ) ≤ tr eθZ .
Therefore,
IP {λmax (Z) ≥ z} ≤ inf e−θz IE tr exp(θZ),

(3)

θ

and (??) follows.
To prove (??) fix θ . Using the spectral mapping theorem one can get that
IEλmax (Z) =

 1

1
1
IEλmax (θZ) = log IE exp λmax (θZ) = log IEλmax exp(θZ) .
θ
θ
θ

Thus we get
IEλmax (Z) ≤

1
log tr IE exp(θZ).
θ

(4)

P
The final step in proving the master inequalities is to bound from above IE tr exp ( ni=1 Si ) . To
do this we use Jensen’s inequality for the convex function tr exp(H + log(X)) (in matrix X ),
where H is deterministic Hermitian matrix. For a random Hermitian matrix X one can write
IE tr exp(H + X) = IE tr exp(H + log eX ) ≤ tr exp(H + log IEeX ).

(5)

Convexity of function (tr exp(H + log(X))) is followed from
tr exp(H + log(X)) = max[tr(Y H) − (D(Y ; X) − tr X)],
Y 0

where D(Y ; X) is relative entropy
D(Y ; X) = φ(X) − [φ(Y ) + h∇φ(Y ), X − Y i],

φ(X) = tr(X log X)

due to the partial maximum and D(Y ; X) are concave functions.
DenotePby IEi the conditional expectation with respect to random matrix Xi . To bound
IE tr exp ( ni=1 Si ) we use (??) for the sum of independent Hermitian matrices by taking the
conditional expectations with respect to i -th matrix:
!
!
n
n−1
X
X
IE tr exp
Si = IEIEn tr exp
Si + Sn
i=1

i=1

≤ IE tr exp

n−1
X

!
Si + log(IEn exp(Sn ))

i=1

≤ tr exp

n
X

!
log IEeθSi

.

i=1

To complete the prove of the Master’s theorem combine (??) and (??) with (??).
The next statement was taken from Koltchinskii (2013) with the proof sketch.
30

(6)

Lemma 12 (Bernstein inequality for moment restricted matrices). Suppose that


kSi kop
2
∀i : IEψ
≤ 1,
M
n
X

2

v =

IESi2

i=1

R = Mψ

−1



2 nM 2
δ v2

,
op


,

δ ∈ (0, 2/ψ(1)).

Then under condition zR ≤ (e − 1)(1 + δ)v2
z2
≥ z} ≤ 2p exp −
2(1 + δ)v2 + 2Rz/3


IP {kZkop
α


.

2

If ψ(u) = eu − 1 then R = M log1/α ( 2δ nM
+ 1).
v2
Proof. According to Master bound one have to estimate IEeθS for S in S1 , . . . , Sn . Denote a
function
eu − 1 − u
f (u) =
u2
Taylor expansion yields
IEeθS ≤ Ip + θ2 IES 2 f (θ kSk)
log IEeθS ≤ θ2 IES 2 f (θ kSk) ≤ θ2 f (θτ )IES 2 + Ip θ2 IE kSk2 f (θ kSk)I(kSk ≥ τ )


  τ −1
kSk   τ −1
2
2
2
IE kSk f (θ kSk)I(kSk ≥ τ ) ≤ M IEψ
ψ
≤ M2 ψ
M
M
M
  −1
δv2
R
2
M ψ
=
M
2n
δv2
log IEeθS ≤ θ2 f (θR)IES 2 + Ip θ2
2n
!
!


n
2
X
X
2
2
2 δv
tr exp
log IE exp(θSi ) ≤ tr exp θ f (θR)IE
Si exp θ
.
2
i=1
i
Consequence. In case ψ(u) = eu − 1 with probability 1 − 2e−x
p
2
z ≤ Rxp + v 5xp ,
3

R ≈ M.

Condition for function ψ(u) = eu − 1 follows from sub-Gaussian moment restriction


kSi kop
IE exp 2
≤ 2.
M
Lemma 13 (Deviation bound for matrix convolution with sub-Gaussian weights). Let a set of
symmetric [p × p] matrices (A1 , . . . , An ) satisfy
X

A2i ≤ v2 .

i

31

Let εi be independent sub-Gaussian, i = 1, . . . , n .
Z =

n
X

εi Ai

i=1

fulfills


p
2
IP kZkop ≥ 2xp v ≤ 2e−x ,
where xp = x + log p.
Proof. Apply the Master inequalities for the case
1
1
log IEeεi Ai = Ui Λ2i UiT = A2i ,
2
2
(
)


X
1 2 2
θεi Ai
tr exp
log IEe
≤ p exp
θ v .
2
i

Ai = Ui Λi UiT ,

Lemma 14 (Deviation bound for rank one matrix convolution with sub-Gaussian weights). Let
vectors u1 , . . . , un in IRp satisfy
kui k ≤ δ
for a fixed constant δ . Let εi be independent sub-Gaussian, i = 1, . . . , n . Then for each vector
b = (b1 , . . . , bn )> ∈ IRn , the matrix Z1 with
Z1 =

n
X

εi bi ui u>
i

i=1

fulfills

IP kZ1 kop

√


≥ δ kbk 2x ≤ 2e−x .
2

2 /2

Proof. As εi are i.i.d. standard sub-Gaussian and IEeaεi ≤ ea
the Master inequality and property
exp(εi uu> ) =

for |a| < 1/2 , it follows from

ui u>
i
exp(εi u2 )
u2

that


IP kZ1 kop ≥ z ≤ 2 inf e
θ>0

−θz

X

n
>
tr exp
log IE exp(θεi bi ui ui )
i=1

X

n
θ2 b2i kui k4 ui u>
i
−θz
≤ 2 inf e tr exp
.
θ>0
2
kui k2
i=1
2
Moreover, as kui k ≤ δ and Ui = ui u>
i /kui k is a rank-one projector with tr Ui = 1 , it holds

 2X

 2 4X

n
n
θ δ
θ
θ2 δ 4 kbk2
2
4
2
tr exp
bi kui k Ui ≤ exp tr
bi Ui = exp
.
2 i=1
2 i=1
2
32

√
This implies for z = δ 2 kbk 2x
IP kZ1 kop



1 2 4 2
≥ z ≤ 2 inf exp −θz + θ δ kbk = 2e−x
θ>0
2


and the assertion follows.

5.13.

Covariance matrices

Consider a sequence independent random variables {εi εTj }, cor(εi , εj ) = Σij , flatted into one
vector ε. The subject of interest is upper bound for operator norm of
U blockDiag (εεT )U T − Iq ,
X
U blockDiag (εεT )U T =
Ui εi εTj UjT , UΣU T = Iq ,

UiT Uj

op

≤ δ2.

ij

Analogically divide ε into mean and stochastic parts
ε = IEε + (ε − IEε) = B + ζ.
Then initial term includes three parts:
U blockDiag(BB > )U >
+ U blockDiag(ζB > )U > + U blockDiag(Bζ T )U >
+ U blockDiag(ζζ T − Σ)U > .
Estimate successively each component.
U blockDiag(BB > )U >

op

= sup
γ

≤ sup
γ

n
X

γ T Ui Bi BjT UjT γ

i=1
n
X

Bi2 UiT Uj

op

≤ δ 2 kBk2 .

i=1

For the second and the third component one may apply
P Master bound 11, in which one have to
estimate exponential moments of each element of the i Ui Ai UiT :
tr log IE exp{Ui Ai UiT }.
With condition IEAi = 0 an intuition hint is
1
T
3
log IEeUi Ai Ui ≤ IE(Ui Ai UiT )2 + O( Ui Ai UiT op ).
2
Consequently by means of Bernstein matrix inequality 12 one have to restrict the second moment
and tail by log IE exp{ Ui Aij UjT op }.
X

v2 =

IE(Ui Aij UjT )2

ij

2
vζζ
/δ 2 =

X

≤ max UiT Uj
i

Ui IEA2ij UjT

ij

,
op

Ui IE(ζi ζjT − Σij )2 UjT

≤ max kΣij kop

op

IE(ζiT Σij−1 ζj )ζi ζjT Σij−1


2

≤ max kΣij kop (λ − 1) +
ij

op

op

ij

ij

X

−I

IE(ζiT Σij−1 ζj )2
33

op

1I(ζiT Σij−1 ζj


> λ) ,

λ > 1.

Upper bound for IE kξk4 1I(kξk >
2
= 9δ 2 maxij kΣij kop p.
vζζ

√

λ) from SGI with ξ 2 = ζiT Σij−1 ζj leads to asymptotic

X

2
vζB
/δ 2 =

Ui IE(ζi BiT )2 UjT

ij

op

≤ max kΣij kop
ij

X

UjT Ui

op

Bi2

ij

≤ max kΣij kop δ 2 kBk2 .
ij

Es for exponential moments for tails restriction
(
log IE exp{ Ui Aij UjT

op

} = log IE exp

)
sup uTi Aij uj

,

kui k2 , kuj k2 ≤ δ 2 .

ui ,uj

Sub-Gaussian properties for exponential moments (Lemma 17) lead to ∃Mζζ , MζB :
(
)
 2

2 T
2δ
T
2
log IE exp sup
ui (ζi ζj − Σij )uj ≤ log IE exp
(kζi k − λmin (Σii ))
Mζζ
ui ,uj Mζζ
≤ SG{kξk2 , 2δ 2 kΣii kop /Mζζ } −
(

2 T
ui ζi BjT uj
log IE exp sup
M
ui ,uj
ζB

)
≤ SG{kξk, 2δ 2

where Mζζ = 3δ 2 p kΣii kop and MζB = 3δ 2 p
−x

12 with probability 1 − 2e

2δ 2
λmin (Σii ) ≤ log(2),
Mζζ

q

q

kΣii kop kBj k /MζB } ≤ log(2),

kΣii kop kBi k. Finally, as a consequence of Theorem

and xq = x + log(2q) and R∗∗ ≈ M∗∗

p
2
≤
R
x
+
2v
5xq ,
ζB
q
ζB
op
3
p
2
U blockDiag(ζζ > )U > op ≤ Rζζ xq + 2vζζ 5xq .
3
So, the summarized error with probability 1 − 2e−x is
p
2
U blockDiag (εεT )U T − Iq op ≤ Rεε xq + 2vεε 5xq + δ 2 kBk2 ,
3
where


q
2
Rεε ≈ MζB + Mζζ = 3δ p max kΣii kop + kΣii kop kBk∞
U blockDiag(ζB > )U >

(ErrVD)

i

and
vεε = vζB + vζζ

5.14.

q
√
= δ max kΣii kop (3 p + kBk).
i

Sub-Gaussian vectors

Consider vector ξ which has restricted exponential or sub-Gaussian moments: ∃ g > 0 :

log IE exp γ > ξ ≤ kγk2 /2,
γ ∈ IRp , kγk ≤ g.
(SG)
Define short notation for upper bounds of such characteristics
log IE exp λX) ≤ SG(X, λ).
For ease of presentation, assume below
large, namely, 0.3g ≥
√ that g is sufficiently
2
examples of an i.i.d. sample, g  n . Define xc = g /4.
34

√
p . In typical

√
p . Then for each x > 0

IP kξk ≥ z(p, x) ≤ 2e−x + 8.4e−xc 1I(x < xc ),

Lemma 15. Let (SG) hold and 0.3g ≥

(7)

where z(p, x) is defined by
1/2
√
p + 2 px + 2x
, x ≤ xc ,
−1
g + 2g (x − xc ),
x > xc .

(
def

z(p, x) =

Usually the second term in (??) can be simply ignored. Define

def
p = tr B ,

def

v2 = tr(B 2 ),


def
λ = λmax B .

def

zc2 = p + vg + λg2 /2,
p
p/λ + gv/λ + g2 /2
def
.
gc =
1 + v/(λg)
p
Lemma 16. Let (SG) hold and 0.3g ≥ p/λ . Then for each x > 0

IP kB 1/2 ξk ≥ z(B, x) ≤ 2e−x + 8.4e−xc 1I(x < xc ),
where z(B, x) is defined by
(p
p + 2vx1/2 + 2λx, x ≤ xc ,
z(B, x) =
zc + 2λ(x − xc )/gc , x > xc .
def

p
√
√
p + 2vx1/2 + 2λx can be upper bounded by p + 2λx :
(√
√
p + 2λx,
x ≤ xc ,
z(B, x) ≤
zc + 2λ(x − xc )/gc , x > xc .

The upper quantile z(B, x) =

Lemma 17. Suppose (SG). For any µ < 1 with g2 > pµ , it holds
 µkξk2  

p
IE exp
1I kξk ≤ g/µ − p/µ ≤ 2(1 − µ)−p/2 .
2
If g is sufficiently large than approximately
 µkξk2 
IE exp
≤ (1 − µ)−p/2 .
2
Lemma 18. Suppose (SG) and kBkop = 1 . For any µ < 1 with g2 /µ ≥ p , it holds
p


IE exp µkB 1/2 ξk2 /2 1I kBξk ≤ g/µ − p/µ ≤ 2det(IIp − µB)−1/2 .
−p/2
Proof. With cp (B) = 2π
det(B −1/2 )
Z


1
cp (B) exp γ > ξ −
kB −1/2 γk2 1I(kγk ≤ g)dγ
2µ
 µkB 1/2 ξk2  Z
 1
= cp (B) exp
exp − µ1/2 B 1/2 ξ − µ−1/2 B −1/2 γ
2
2
 µkB 1/2 ξk2 

= µp/2 exp
IPξ kµ−1/2 B 1/2 ε + B 1/2 ξk ≤ g/µ ,
2
35

2



1I(kγk ≤ g)dγ

(8)

where ε denotes a standard normal vector in IRp and IPξ means the conditional probability
 given
ξ . Moreover, for any u ∈ IRp and r ≥ p1/2 + kuk , it holds in view of IP kB 1/2 εk2 > p ≤ 1/2

√ 
IP kB 1/2 ε − uk ≤ r ≥ IP kB 1/2 εk ≤ p ≥ 1/2.
This implies


p

exp µkB 1/2 ξk2 /2 1I kBξk ≤ g/µ − p/µ
Z


1
−p/2
≤ 2µ
cp (B) exp γ > ξ −
kB −1/2 γk2 1I(kγk ≤ g)dγ.
2µ
Further, by (SG)
Z



1
kB −1/2 γk2 1I(kγk ≤ g)dγ
exp γ > ξ −
2µ
Z

 kγk2
1
−
kB −1/2 γk2 dγ
≤ cp (B) exp
2
2µ

cp (B)IE

≤ det(B −1/2 ) det(µ−1 B −1 − Ip )−1/2 = µp/2 det(IIp − µB)−1/2
and (??) follows.
The next object of interest in this topic is IEkξkr 1I(kξk > t). Rather useful form of it is
Z+∞
IEkξk 1I(kξk > t) = IP (kξk > t)t + r
IP (kξk > t)tr−1 dt.
r

r

t

With x0 = (t −

√

p)2 /2

Z+∞
Z+∞
√
√
√
IP (kξk > t)tr−1 dt =
2e−x ( 2x + p)r−1 d 2x ≤
t

2e−x0 tr−1
1 − (r − 1) log(x0 )/x0

x0

Consequently, if (r − 1) log(x0 )/x0 ≤ 1/2, than
√

IEkξkr 1I(kξk > t) ≤ 2e−(t−

p)2 /2


tr + 2rtr−1 .

(SGI)

Analogically one is able to restrict moment with exponent part and r log(x0 )/x0 + α ≤ 1/2
√ 2

(SGIexp)
IEkξkr eαkξk 1I(kξk > t) ≤ 4e−(t− p) /2+αt tr + rtr−1 .

5.15.

Gaussian quadratic forms

The next result explains the concentration effect of γ > Bγ for a standard Gaussian vector γ
and a symmetric matrix B . We use a version from Laurent and Massart (2000) with a complete
proof.
Theorem 17. Let γ be a standard normal Gaussian vector and B be symmetric positively
definite p × p -matrix. Then with p = tr(B) , v2 = tr(B 2 ) , and λ = kBkop , it holds for each
x≥0


>
2
IP γ Bγ > z (B, x) ≤ e−x ,
(9)
def

z(B, x) =
36

p
p + 2vx1/2 + 2λx.

In particular, it implies

IP kB 1/2 γk > p1/2 + (2λx)1/2 ≤ e−x .
Also

IP γ > Bγ < p − 2vx1/2 ≤ e−x .

(10)

If B is symmetric but non necessarily positive then

IP γ > Bγ − p > 2vx1/2 + 2λx ≤ 2e−x .
Proof. Normalisation by λ reduces the statement to the case with λ = 1 . Further, the standard
rotating arguments allow to reduce the Gaussian quadratic form kγk2 to the chi-squared form:
>

γ Bγ =

p
X

λj νj2

j=1

with independent standard normal r.v.’s νj . Here λj ∈ [0, 1] are eigenvalues of B , and p =
λ1 +. . .+λp , v2 = λ21 +. . .+λ2p . One can easily compute the exponential moment of (γ > Bγ −p)/2 :
for each positive µ < 1
p

1 X
log IE exp µ(γ > Bγ − p)/2 =
−µλj − log(1 − µλj ) .
2 j=1

Lemma 19. Let µλj < 1 and λj ≤ 1 . Then
p

µ2 v2
1 X
.
−µλj − log(1 − µλj ) ≤
2 j=1
4(1 − µ)
Proof. In view of µλj < 1 , it holds for every j
−µλj − log(1 − µλj ) =

∞
X
(µλj )k
k=2

k

∞
(µλj )2 X
(µλj )2
(µλj )2
≤
(µλj )k ≤
≤
,
2 k=0
2(1 − µλj )
2(1 − µ)

and thus
p
p
X
1 X
(µλj )2
µ2 v2
−µλj − log(1 − µλj ) ≤
≤
.
2 j=1
4(1 − µ)
4(1 − µ)
j=1

The next technical lemma is helpful.
Lemma 20. For each v > 0 and x > 0 , it holds


inf −µ vx1/2 + x +
µ>0

37

µ2 v2
4(1 − µ)


≤ −x.

(11)

Proof. Let pick up
µ = 1−

1
x1/2
=
,
2x1/2 /v + 1
x1/2 + v/2

so that µ/(1 − µ) = 2x1/2 /v . Then

−µ vx1/2 + x +

µ2 v2
4(1 − µ)


= −µ vx1/2 + x + v2 /4 +
=−

µv2
4(1 − µ)

2 2x1/2 v
x1/2
1/2
x
+
v/2
+
= −x
x1/2 + v/2
4

and the result follows.
Now we apply the Markov inequality


log IP γ > Bγ > p + 2vx1/2 + 2x = log IP (γ > Bγ − p)/2 > vx1/2 + x




1/2
>
≤ inf −µ vx + x + log IE exp µ(γ Bγ − p)/2
µ>0



1/2

≤ inf −µ vx
µ>0


+x +

µ2 v2
4(1 − µ)


≤ −x

and the first assertion (??) follows. The second statement follows from the first one by tr(B 2 ) ≤
kBkop tr(B) = λ p .
Similarly for any µ > 0

 µ
√ 
√ 
IP γ > Bγ − p < −2v x ≤ exp −µv x IE exp − (γ > Bγ − p) .
2
By (??)
p

1 X
log IE exp −µ(γ > Bγ − p)/2 =
µλj − log(1 + µλj ) .
2 j=1

and
p

p

∞

p

1 X
1 X X (−µλj )k X (µλj )2
µ2 v2
µλj − log(1 + µλj ) =
≤
=
.
2 j=1
2 j=1 k=2
k
4
4
j=1
√
Here the choice µ = 2 x/v yields (??).
One can put together the arguments used for obtaining the lower and the upper bound for
getting a bound for a general quadratic form γ > Bγ , where B is symmetric but not necessarily
positive.

5.16.

Sandwich lemma

Lemma 21. Let differentiable measure IP [ depends on r.v. from a continuous measure IP , z =
(z1 , . . . , zK ) is a multivariate quantile. Assume following error in total variation distance between
measures




[
IP z1 , . . . , zK − IP z1 , . . . , zK < δ.
38

Then each quantile zk[ (α), 1 ≤ k ≤ K from measure IP [ may be bounded by quantile from
measure IP :
zk (α + δ) ≤ zk[ (α) ≤ zk (α − δ),
where



IP zk (α) = 1 − α,

IP

[



zk[ (α)


= 1 − α.

And if q [ is the multiplicity correction parameter such that


[
[
[
IP z1 (q α), . . . , zK (q α) = 1 − α,
then

IP z [ (q [ α) − (1 − α) ≤ (2K + 1) δ.
Proof. Define two sets
Z+ (δ) = {z : IP (z) ≤ 1 − α + δ},
Z− (δ) = {z : IP (z) ≥ 1 − α − δ}.
For all points z from Z+ ∩Z− it holds that |IP (z) − (1 − α)| ≤ δ. If IP [ (z [ ) = 1 then z [ ∈ Z+ since
for all fixed z ∈ IRK \ Z+ : IP [ (z) > IP (z) − δ ≥ 1 − α. Analogically z [ ∈ Z− and z [ ∈ Z+ ∩ Z− .
In case K = 1 one can choose non-random quantiles in the border of Z+ ∩ Z− which will bound
z [ . So each component of z [ could be bounded in the same way:
zk (α + δ) ≤ z [ (α) ≤ zk (α − δ).
In case K > 1 the these bounds become random because of multiplicity correction which involves
random multiplier q [ . We have to bound zk (q [ α + δ) by a non-random quantile in order to use it
as an argument for measure IP .
z(α + δ) ≤ z [ (α) ≤ z(α − δ)
z [ (q [ α) ≤ z(q [ α − δ) ≤ z [ (q [ α − 2δ)




[
[
[
[ [
1 − α ≤ IP z(q α − δ) ≤ IP z (q α − 2δ)
Lemma 22. For differentiable measure IP (ξ < x) and event A:
IP (ξ < x, A)0x
≤ 1.
IP (ξ < x)0x



 X
[
[ [
[
[ [
IP z (q α − 2δ) = IP z (q α) +
(IP [ )0z[ (zk[ )0 2δ,
k

i

(IP

[ 0
)zk (zk[ )0

=

[
(IP [ (z1[ , . . . , zk[ , . . . , zK
))0z[

k

[ (z [ ))0
k z[

(IP
k


IP [ z [ (q [ α − 2δ) ≤ 1 − α + 2Kδ.

≤ 1,



[
1 − α ≤ IP z(q α − δ) ≤ 1 − α + 2Kδ,
[



[
1 − α − 2Kδ ≤ IP z(q α + δ) ≤ 1 − α.
[

39

According to the arguments from the beginning of the proof z(q [ α − δ) and z(q [ α − δ) belongs to
Z+ (2Kδ + δ) ∩ Z− (2Kδ + δ). Due to one dimensional parametrization if z(q [ α − δ) there exist two
fixed points on the border of Z+ (2Kδ + δ) ∩ Z− (2Kδ + δ) such that
z+ = max z(q [ α − δ),

z− = min z(q [ α + δ).

Finally,
z− ≤ z(q [ α + δ) ≤ z [ (q [ α) ≤ z(q [ α − δ) ≤ z+ ,
and subsequently
1 − α − (2K + 1)δ ≥ IP (z− ) ≤ IP (z [ (q [ α)) ≤ IP (z+ ) ≤ 1 − α + (2K + 1)δ.

5.17.

TP variance
√

P
p then t Pτ (t)IEkξlr (t)k = 0 and
v

!
u 2h
2h
uX
X
X
Var
Pτ (t)kξlr (t)k ≤
Pτ2 (t) Var t
kξlr (t)k2  .

Assume here that

P

t

Pτ (t) = 0 and IEkξlr (t)k =

t

t=1

t=1

 v

2
v
u 2h
u 2h
uX
uX
kξlr (t)k2  = 2hp − IE t
kξlr (t)k2  .
Var t
t=1

t=1

Since each kξlr (t)k2 is close to normal random vector norm, its lower bound may be taken from
Theorem 17. So with probability e−x for all 1 ≤ t ≤ 2h
√
kξlr (t)k2 > p − 2 p(x + log(2h))1/2 .
Then with probability e−x
v

u 2h
uX
√
kξlr (t)k2  ≤ 4h p(x + log(2h))1/2 .
Var t
t=1

Consider triangle pattern example

0,



(t − τ )/h + 1/2,
Pτ (t) =

(τ − t)/h + 1/2,



0,

t < τ − h,
τ − h ≤ t ≤ τ,
τ ≤ t ≤ τ + h,
t > τ + h.

In this case through integral sum one get
2h
X
t=1

1
Pτ2 (t) ≈ h.
6

Finally
!
Var

X

Pτ (t)kξlr (t)k

t

and

r
zh (x) ≤

2 √
≤ h2 p(x + log(2h))1/2 ,
3

2 1/4
hp (x + log(2h))1/4 ex/2 .
3
40

(P)

The abrupt type change point statistic without noise component has view
p
2Th (t) = (Pτ (t) + 1/2)∆ ± 7♦, τ − h ≤ t ≤ τ + h.
The sufficient condition for change point detection of size ∆ in position τ is
τ +h
X

p
Pτ (t) 2Th (t) > zh (x).

t=τ −h

Equivalently
7
1
∆h + h♦ >
6
2

r

2 1/4
hp (x + log(2h))1/4 ex/2 .
3

41

