Estimating a change point in a sequence of very high-dimensional covariance matrices
H. Dette G. M. Pan Q. Yang

=16.5cm
=23.0cm
=-1cm
=0.1cm
=7.0mm
























equationsection

    equationsectionplainthmTheorem[section]
                
                
                propProposition[section]
                lemLemma[section]
                coroCorollary[section]
                rmkRemark[section]
                egExample[section]
                clmClaim[section]
                deffDefinition[section]
                assuAssumption[section]
                condCondition[section]
                
                
                
                
                
                
                
                
                Estimating a change point in a sequence of very high-dimensional covariance matrices
                    H. Dette, G. M. Pan and Q. Yang
                    
                ====================================================================================
                
                
                This paper considers the problem of estimating a change point in the covariance matrix in a sequence of high-dimensional vectors, where the dimension is substantially larger than the sample size.
                
                A two-stage approach is proposed to efficiently estimate the location of the change point. The first step consists of a   reduction of the  dimension to identify elements of the covariance matrices corresponding to significant changes. In a second step we use   the components after dimension reduction to determine the position of  the change point. Theoretical properties are developed for both steps and numerical studies are conducted to support the new methodology.
                
                Keywords: High-dimensional covariance matrices; change point analysis; dimension reduction.
                
                
                
                
                
                § INTRODUCTION
                
                
                Change point detection has a long history having it origins in quality control [see <cit.> or <cit.> for early references] and
                it has been an active field of research until today since the phenomena of sudden changes arise in various areas, such as financial data
                (house market, stock), signal processing, genetic engineering, seismology, machine learning.
                 In the last decades numerous authors have worked on this problem from several perspectives including the construction of tests for   the hypotheses of the existence of change points
                and the estimation of their locations. We refer to <cit.> and <cit.> for some recent reviews on this subject.
                
                An important  problem in the detection of structural breaks in multivariate data
                 is the detection of changes in a sequence of  means.  <cit.>, <cit.>, <cit.> and <cit.> investigated  this problem   using different
                variants of CUSUM statistics. More recently, the   high-dimensional case
                 (dimension larger than the sample size) has been discussed by several authors as well. <cit.> considered a maximum of statistics across panels coordinate-wise
                 to test the hypothesis of at least one change point in a sequence of high dimensional mean vectors
                 [see also <cit.> who studied   relevant changes in this context]. <cit.>  suggested sparsified binary segmentation for this problem, while <cit.>
                 investigated a double CUSUM approach  transferring - roughly speaking -  the high-dimensional data to a univariate CUSUM statistic.
                 We also mention the work of <cit.>, who looked at the problem under sparse alternatives and   <cit.>, who
                 considered the situation, where  at  certain  time  points,  the
                mean structure changes in a sparse subset of the coordinates.
                
                
                
                While substantial effort has been spent on change point analysis for the multivariate mean, the problem of detecting structural breaks
                in the covariance matrix  has not been studied so intensively in the literature. For a fixed dimension, say p,  <cit.> developed  nonparametric change point analysis based on the well known CUSUM approach.
                 <cit.> proposed a general approach to detect relevant change points in a parameter of a time series. In an online supplement to this paper a test for a relevant change in the covariance matrix is proposed, where the dimension is also fixed.
                Recently, <cit.> considered the case where the dimension is increasing with the sample size and
                demonstrated by means of a simulation study that tests of stability of the whole covariance matrix have severe size distortions. As an alternative they proposed and investigated change point analysis
                based on PCA. In an unpublished preprint  <cit.>  also looked at the high-dimensional setting and suggested  a multiscale approach
                under sparsity assumptions, while <cit.> considered the problem of detecting multiple change points in the situation p=O(n/log n) (here n is the sample size) and
                investigated optimality properties of the   binary segmentation  [see <cit.>] and the wild binary segmentation  algorithm [see <cit.>] for localising  multiple changes points in a sequence
                of high dimensional covariance matrices.
                
                The purpose of the present  paper is to propose an alternative  estimator of the change point  in  a sequence of  very high dimensional  covariance matrices and to investigate its theoretical and empirical properties,  where we do not impose any sparsity assumptions on the matrices.
                 To be precise, suppose that
                 _1,⋯,_k_0, _k_0+1,⋯,_n
                 are p-dimensional   observations with common mean vector  and existing covariance matrices.
                 The parameter k_0 defines  the true change point in the structure of the covariance matrices. That is,
                the first k_0 observations have covariance matrix Σ_1=(σ_1(a,b))_a,b=1,…,p∈^p× p,
                while the last n-k_0 observations have covariance matrices Σ_2=(σ_2(a,b))_a,b=1,…,p∈^p× p and Σ_1≠Σ_2. We are interested in estimating the point k_0. One difficulty in dealing with changes in the covariance matrix is the dimensionality since there are p(p+1)/2 positions needed to be compared.  This brings in much noisy information when there are many equal components in the matrices Σ_1 and Σ_2, leading to a loss of accuracy in detection of the change point. Taking this consideration into account, we propose to proceed in two steps to identify the change point.
                First, we  apply a dimension reduction technique reducing the dimension from p(p+1)/2 in the original problem to a substantially
                 smaller value, say m.   Roughly speaking, we   only keep the components
                 in the analysis for which a weighted mean of the squared differences of the covariance estimators from the samples
                 _1,⋯,_k and  _k+1,⋯,_n exceeds a given threshold (the mean is calculated summing with respect to  the different values of the potential change points of k).
                 Therefore our approach is vaguely related to the estimation of sparse covariance matrices, which has found considerable attention in the literature
                 [see <cit.>,  <cit.>  or <cit.> among many others]. However, in contrast to this work, we do not assume a sparse structure of the covariance matrix, but identify important components by  thresholding  a weighted sum of the (squared) differences corresponding to all potential samples before and after a postulated change point.
                In a second step after dimension reduction, we use a CUSUM type statistic based on the reduced components to  locate the change point.
                
                An outline of the paper is given as follows. In Section <ref>, we introduce our main methodology – both the dimension reduction step and the detection step.
                In particular, a  bootstrap method is suggested to select the threshold used for the dimension reduction  (see the discussion in Section <ref>).
                Theoretical results are developed in Section <ref>, where we prove that (asymptotically) we identify all relevant components correctly and that we estimate the location of the change k_0 consistently.     In Section <ref> we investigate the finite sample properties of the new method and demonstrate that it yields precise estimates of the change point in situations, where the dimension is substantially larger than the sample size. We also provide a comparison with two alternative methods which are most similar in spirit to our approach and have recently been proposed in the literature.
                Finally, all proofs and technical details are deferred to an appendix in Section  <ref>.
                
                
                
                
                § METHODOLOGY
                
                Let  _1,⋯,_k_0,_k_0+1,⋯,_n denote a sample of independent p-dimensional random vectors   with common mean  and existing covariance matrices. The position k_0 is the “true” change point of the covariance matrices, i.e.
                the first k_0 random variables _1,⋯,_k_0  have covariance matrix Σ_1=(σ_1(a,b))_a,b=1,…,p∈^p× p,
                while the last (n-k_0) random variables _k_0+1,⋯,_n  have covariance matrices Σ_2=(σ_2(a,b))_a,b=1,…,p∈^p× p and Σ_1≠Σ_2.  Our aim is to estimate  the location k_0 of the change.
                For this purpose  we proceed in two steps.
                
                  *  Step 1 consists of a dimension reduction. If
                
                    Σ̂^k_1 = (σ̂^k_1 (a,b ) )^p_a,b=1
                
                   and Σ̂^n_k+1 = ( σ̂_k+1^n  (a,b) )^p_a,b=1 denote the respective estimators of the covariance matrices from the data _1,…,_k and _k+1,…, _n,
                    we - roughly speaking - only keep components in the change point analysis for which the quantity
                 
                    ∑^n-2_k=2 k(n-k) ( σ̂^k_1(a,b)  - σ̂^n_k+1 (a,b) )^2
                
                 is sufficiently large.
                
                  *  Step 2 consists of the  detection of a change point in the data obtained after dimension reduction. For this purpose let σ̃^k_1 and σ̃^n_k+1 denote the vectors  containing all  elements of the matrices  Σ̂^k_1 and Σ̂^n_k+1 corresponding to components  which have been identified in the first step of the procedure. Then - roughly speaking - we propose to estimate the change point by maximizing the statistic
                    
                    Ũ_n(k) =  (n-k)k (σ̃^k_1 - σ̃^n_k+1 ) ^2_2
                
                    where ·_2 denotes the Euclidean norm.
                
                We will give a detailed explanation of these two steps in the following subsections, where the statistics under consideration will be slightly modified.
                The proposed methodology depends on a regularisation parameter, say τ,  determining the amount of dimension reduction for Step 1,  and in Section <ref> we introduce a data-driven
                method for choosing this  threshold.
                
                
                
                 §.§ Dimension reduction  
                
                For i=1,…,n denote  by _i=(X_i1,⋯,X_ip)^T the ith observation, let =1/n∑_i=1^n _i=(X̅_1,⋯,X̅_p) be the sample mean and define
                
                    _i=_i- =  ( Ẋ_i1, … , Ẋ_ip )^T  =    (X_i1-X̅_1 , …  , X_ip-X̅_p  )^T
                
                as the vector of centered observations. We introduce the following  statistic
                
                    V_k=   ( V_k (a,b)  ) _1 ≤ a ≤ b ≤ p   =   1/k(k-1)∑∑_i≠ j≤ k(_i_i^T)∘(_j_j^T) 
                          +1/(n-k)(n-k-1)∑∑_i≠ j>k(_i_i^T)∘(_j_j^T)      -2/k(n-k)∑_i≤ k∑_j>k(_i_i^T)∘(_j_j^T),
                 where “vech(H)” indicates the half-vectorization p(p+1)/2 vector by vectorizing only the lower triangular part of the symmetric matrix H=(H(a,b))^p_a,b=1 and “∘” is the Hadamard product (or entrywise product) of the vectors  and . Obviously,
                 V_k is a p(p+1)/2-dimensional vector.
                
                
                We first give an intuitive illustration of the motivation behind the construction of the statistic V_k defined in (<ref>), which
                is in fact  motivated by being an  approximation of the statistic
                
                    Ṽ_k   =   1/k^2∑_i, j=1^k(_i_i^T)∘(_j_j^T)
                    +1/(n-k)^2∑_i, j=k+1^n(_i_i^T)∘(_j_j^T)       -2/k(n-k)∑_i=1^k∑_j=k+1^n(_i_i^T)∘(_j_j^T) 
                       =    ( 1/k∑_i=1^k(_i_i^T)  - 1/n-k∑_i=k+1^n(_i_i^T)  ) ∘ ( 1/k∑_i=1^k(_i_i^T)  - 1/n-k∑_i=k+1^n(_i_i^T)  ).
                
                Note that   the vector Ṽ_k coincides with the vector
                
                    (  ( σ̂^k_1(a,b) - σ̂_k+1^n(a,b)  )^2    )_1 ≤  a ≤ b ≤ p
                
                of the squared  (componentwise) differences of the elements of the covariance estimators
                Σ̂^k_1 = 1/k∑^k_i=1_i ^T_i and
                Σ̂^n_k+1 = 1/n-k∑^n_i=k+1_i ^T_i. Consequently,
                at the “true”  change point   k=k_0, one can verify that Ṽ_k_0 is an
                estimator of vech(Σ_1-Σ_2)^2, which will be used to measure the difference between the two population covariance matrices.
                The difference between  Ṽ_k and V_k consists  in the fact that  in the statistic V_k we omit the terms corresponding to i=j to eliminate the influence of  the
                covariances of the random variables (_i_i^T).
                
                However, in the change point problem we actually do not know the location of k_0, and we have to consider all the positions k
                as long as the statistic is well defined. In particular, we obtain for the expectation of the  component V_k(a,b) of the vector V_k
                    corresponding to the   position (a,b) in the matrices Σ_1 and Σ_2
                    V_k(a,b) = {[                               (1-2/n)^2(σ_1(a,b)-σ_2(a,b))^2                                                       k=k_0,;             k_0(k_0-1)/k(k-1)(1-2/n)^2(σ_1(a,b)-σ_2(a,b) )^2                                                       k>k_0,; (n-k_0)(n-k_0-1)/(n-k)(n-k-1)(1-2/n)^2(σ_1(a,b)-σ_2(a,b) )^2                                                       k<k_0. ].
                
                 Note also, that the
                  mean of V_k(a,b) always achieves the largest value 
                    (1-2/n)^2(σ_1(a,b)-σ_2(a,b) )^2
                
                  at the true change point position k_0 because the coefficients k_0(k_0-1)/k(k-1) and (n-k_0)(n-k_0-1)/(n-k)(n-k-1) before (1-2/n)^2(σ_1(a,b)-σ_2(a,b))^2 are smaller than 1  when k≠ k_0. Moreover, these coefficients  are only related to
                k and are not influenced by the position (a,b). Consequently, for any fixed k, larger values of V_k(a,b) indicate a larger difference between
                σ_1(a,b) and σ_2(a,b), thus implying a significant component.
                Additionally, instead of investigating each value of k separately, we suggest a weighted sum
                
                    D=1/n-3∑_k=2^n-2k(n-k)/nV_k,
                
                  to identify   the largest   components among  the p(p+1)/2 entries.
                 The weights k(n-k)/n are introduced to address the different sizes of   the variance of V_k for different values of k. By selecting the largest entries in the vector D, we are able to identify the components with the largest changes.
                
                
                In view of this discussion, we conduct the dimension reduction as follows.
                Let D(a,b) denote the elements of the vector D corresponding to the position (a,b) in the matrices Σ_1 and Σ_2.
                We determine   all  components which are larger than a critical value τ, which will be   specified in Section <ref>,   define
                
                    D_τ={(a,b): D(a,b) > τ, 1≤ a≤ b ≤ p}
                
                as the set of all corresponding components and denote by m=# D_τ  its cardinality.
                In this way, we reduce the p(p+1)/2-dimensional vector to
                a vector of   dimension m.    In the next step  we will simply work with the m-dimensional vectors corresponding to the components  identified by the set D_τ. In Theorems <ref> and <ref>, it will be shown that after dimension reduction with an appropriate
                threshold all  entries with no difference are discarded, while all the entries with a sufficiently large difference   are kept.
                
                
                
                
                
                
                
                
                 §.§ Change point detection after dimension reduction
                
                For the estimation of the change point we propose the test statistic
                
                
                
                
                
                    U_n(k)=1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1
                    (_i-_j)^T(_t-_l),
                
                where _i is an m-dimensional subvector of (_i_i^T), only keeping the m components of the index set
                D_τ defined in (<ref>). Then the estimator of the  change point k_0 is defined by
                
                    k̂=1≤ k≤ nmax  U_n(k).
                
                The motivation behind the construction of U_n(k) stems from the fact the statistic U_n(k) is related to a CUSUM type
                statistic which is frequently used in change point analysis. To be precise, consider the CUSUM statistic
                
                    Ũ_n(k)    =   1/n^4∑^k_i,t=1∑^n_j,l=k+1 (_i - _j)^T (_t - _l) 
                       =   1/n^4{ (n-k)^2 (∑^k_i=1_i )^T (∑^k_j=1_j ) - 2k(n-k) (∑^k_i=1_i )^T (∑^n_j=k+1_j ) 
                               + k^2 ( ∑^n_i=k+1_i )^T ( ∑^n_j=k+1_j ) }
                       =   1/n^4 (n-k) ∑^k_i=1_i - k ∑^n_j=k+1_j  ^2_2
                     = 1/n^4 (n-k)k (σ̃^k_1 - σ̃^n_k+1 ) ^2_2
                
                where
                ·_2 denotes the usual Euclidean norm   and
                σ̃^k_1  and σ̃^n_k+1 respectively denote the vectors containing the elements of the covariance estimators
                 Σ̂^k_1   and Σ̂_k+1^n corresponding to positions identified in the first step.
                 Observing the definition of _i and noting that the difference
                
                    Ũ_n(k) - U_n(k)    =   1/n^4∑^k_i=1∑^n_j=k+1_i - _j ^2_2 +
                     1/n^4∑^k_i=1∑^n_(j≠ l)=k+1 ( _i - _j )^T  ( _i - _l ) 
                               + 1/n^4∑^k_(i≠ t)=1∑^n_j=k+1( _i - _j )^T  ( _t - _j )
                is of smaller order than U_n(k) when k is far from 1 and n, we see that the statistic U_n(k) is a CUSUM type statistic obtained from the components identified in the first step. It is therefore related to the statistic in equation (2.12) in  <cit.>, who proposed an estimator of the change point based on a quadratic form using  ALL elements of the difference Σ̂^k_1 - Σ̂_k+1^n.
                Note that in the definition of U_n we eliminate the influence of the covariances cov(^T_i_i) by omitting terms corresponding to i=t and j=l in Ũ_n(k). As a consequence, we avoid the estimation of such   higher order moments.
                
                
                We will show in Theorem <ref> that - under appropriate regularity conditions -  the statistic
                k̂ in (<ref>) is in fact a consistent estimator   of the change point  if the dimension and sample size converge to infinity
                and the threshold is chosen appropriately. More precisely,  we can choose τ=C·max(log p,log n) with a sufficiently large constant C
                and  the dimension p can be of polynomial order of the sample size n.
                
                
                
                 §.§ Selecting the threshold τ via resampling 
                
                For a data driven choice of the threshold  τ we propose  a bootstrap   approach, which mimics the distributional properties in the case of no
                change point.
                 To be precise define
                
                    =  (_1,⋯, _n )=  ((_1_1^T),(_2_2^T),⋯,(_n_n^T) )
                    ∈^p(p+1)/2× n,      = 1/√(2) (_2-_1,  _4-_3, … , _2⌊n/2⌋ -_2⌊n/2⌋-1 )  ∈^p(p+1)/2×⌊n/2⌋,
                
                and denote by  _(i)^T  the i-th row of the matrix ,  i.e.
                
                    =   (_(1),⋯, _(p(p+1)/2))^T.
                
                For each ℓ =1,⋯, p(p+1)/2, the ⌊n/2⌋-dimensional vector _(ℓ ) can be considered as a combination of ⌊n/2⌋
                 observations and   we  denote  the empirical  standard deviation of these  ⌊n/2⌋ observations by o_ℓℓ.
                  The construction of the matrix  ensures  that  the means of its columns are zero  except for at most one position
                  (note that the means of the columns of  have one change point at k_0), which does not have a substantial effect on the standard deviation
                 provided that the sample size is not too small.
                If the variance of the random variables in the ℓth  row of the matrix  is constant, it is easy to see that o_ℓℓ^2 estimates
                this variance. If this assumption is not satisfied, generally speaking, o_ℓℓ^2 always estimates n×(T̅_(l)), where T̅_(l) is the average
                of the random variables in the ℓth  row of the matrix . Note that
                 the  factor 1/√(2)  in (<ref>) reflects the fact that the matrix  is formed from differences of two consecutive columns of
                 the matrix .
                
                In order to  estimate  the threshold τ   let Υ denote  the
                 p(p+1)/2×p(p+1)/2  diagonal matrix, with entries o_11, … , o_p(p+1)/2, p(p+1)/2.
                We generate a new data matrix =(y_ij) ∈^p(p+1)/2 × n  with independent  standard normal distributed entries and define
                
                    ^*=Υ.
                
                In other words, in the bootstrap we replace  the quantities (_1_1^T) , … , (_n_n^T)
                 by the n columns of ^∗. As a consequence  the terms  Ẋ_1aẊ_1b , … , Ẋ_naẊ_nb
                are replaced by  o_ℓℓy_ℓ 1,…, o_ℓℓ y_ℓ n where the random variables y_ℓ1,…, y_ℓ n
                are independent standard normal distributed and the index ℓ corresponds to the position (a,b).
                
                
                Next we calculate for the matrix ^*  the quantities V^*_k and D^* defined in (<ref>) and (<ref>) respectively,  and obtain the vector
                
                    D^*=(D^* (a,b) : 1 ≤ a ≤ b ≤ p).
                
                The  threshold  τ is finally defined  as the largest entry of D^*, i.e.
                
                    τ=max_1 ≤ a ≤ b ≤ p D^*(a,b).
                
                
                § ASYMPTOTIC PROPERTIES
                
                
                In this section we discuss the theoretical properties of our approach. For this purpose we need several assumptions, which will be stated first, beginning with conditions on the distribution of the random vectors _i.
                
                
                Denote _i=(X_i1,⋯,X_ip), i=1,⋯,n. For any 1≤ a≤ p, X_ia is a  sub-Gaussian random variable, i.e.
                there are positive constants C_1, C_2 (independent of the indices i and a) such that for every t>0,
                
                    (|X_ia|>t)≤ C_1e^-C_2t^2.
                
                Moreover, the covariance matrices before and after the change point
                satisfy Σ_ν_op≤ M  (ν=1,2) for some  positive constant M, where ·_op denotes the spectral norm.
                Our next assumption specifies the size of the change, which can be detected using the threshold τ. Note that the dimension p is increasing with the sample size
                and a difference between the matrices Σ_1 and Σ_2 might vanish asymptotically if p,n →∞
                although it is visible for any  fixed p (for example if Σ_1 - Σ_2 = e_pe^T_p/p where e^T_p = (0,…, 0,1)).
                
                 The smallest nonzero entry of the matrix Σ_1-Σ_2 satisfies
                    | σ_1(a,b)-σ_2(a,b) | >
                    λ≥ C√(τ/n)max{n^2/(n-k_0)^2, n^2/k_0^2}.
                
                Note that condition (<ref>)
                implies that
                
                    λ≥ C√(τ/n)max{n/n-k_0, n/k_0, √(n/k_0),  √(n/n-k_0), n/k_0√(n-k_0/k_0),n/n-k_0√(k_0/n-k_0)}.
                 For some small positive constant c we have
                    p^2n   =   o(e^cτ) 
                    
                    p^2n   =   o(e^cn^1/4√(τ)) 
                    
                     p^2n^2    =   o(e^c√(τ)).
                
                Recall the definition of the set 𝒟_τ in (<ref>) and define
                
                    N = {(a,b)  1 ≤ a ≤ b ≤ p;  σ_1(a,b) =σ_2(a,b) }
                
                as the set of indices corresponding to equal elements in the matrices Σ_1 and Σ_2. Then    under Assumption <ref>
                    {⋃_(a,b)∈ N{D(a,b)>τ}} =  (  N∩ D_τ≠∅ )
                    ≤ c_1p^2n[e^-c_2τ+e^-c_2n^1/4√(τ)+ne^-c_2√(τ)],
                
                where c_1 and c_2 are some constants. In particular, if Assumption    <ref> is also satsified
                
                    {⋃_(a,b)∈ N{D(a,b)>τ}}
                     → 0,
                
                i.e, after dimension reduction, all the entries with no difference are discarded.
                
                Recall the definition of the set 𝒟_τ in (<ref>) and define
                
                    P = { (a,b) : 1 ≤ a ≤ b ≤ p; |σ_1(a,b)-σ_2(a,b)| > λ}
                
                as the set of components which differ by more than λ. Then under Assumption <ref> and   <ref> we have
                
                    {⋂_(a,b) ∈ P{D(a,b)>τ}} = ℙ (𝒫⊂𝒟_τ) ≥
                    1-c_3p^2n[e^-c_4τ+e^-c_4n^1/4√(τ)+ne^-c_4√(τ)],
                
                where c_3 and c_4 are  constants. In particular, if Assumption    <ref> is also satisfied
                
                    {⋂_(a,b) ∈ P{D(a,b)>τ}}
                    → 1,
                
                i.e, after dimension reduction, all components corresponding to a difference larger than λ are kept.
                
                
                 In our next result we establish the asymptotic consistency of the estimator k̂. Here and throughout this paper the symbol i.p.⟶ denotes convergence in probability.
                
                
                If Assumptions <ref> and   <ref> are satisfied, we have
                
                    {|k̂k_0-1|≥ϵ}≤
                    c_5p^2n[e^-c_6τ+ne^-c_6√(τ)],
                
                where c_5 and c_6 are  constants.
                 In particular, if Assumption    <ref> is also satisfied it follows
                
                    k̂/k_0i.p.⟶ 1.
                 Theorem <ref>, Theorem <ref> and Theorem <ref> are still true if
                Assumption <ref> is replaced by the following Assumption <ref>.
                
                 Assume that k_0>n^ϵ for some 0<ϵ<1 and that there exists a positive constant M<∞ such that
                pn^-M→ 0 and  pn e^-cτ→ 0 for some sufficiently small positive constant c. Note that M could be any large positive constant.
                 Note that we can choose τ=C·max(log p,log n) in  Assumption <ref>, where
                C is a sufficiently large constant. Then the only requirements are  k_0>n^ϵ for the location of the change and
                 that the dimension p cannot exceed a polynomial  order of the sample size n (but the degree of the polynomial can be  arbitrary).
                 Moreover, the inequality (<ref>)  also qualitatively describes a relation  between the location of the change and the
                size of the differences between the elements of the covariance matrices before and after the change.
                For example, if  p>n,  we have τ=C·max(log p,log n)=C·log p  and, if
                 k_0 is proportional to n, this means that the smallest non  zero element of the matrix Σ_1 - Σ_2 should satisfy
                
                    |σ_1(a,b)-σ_2(a,b) |  >  C√(log p/n).
                
                 This is a well-known order to  distinguish signal from noise in  covariance matrix estimation, see for example <cit.>,
                 who considered covariance estimators based on thresholding.
                
                
                On the other hand, the choice τ=C·max(log p,log n) is not  possible in Assumption <ref>.  However,
                if we choose τ=C·max((log p)^2,(log n)^2) with a sufficiently large constant C, there is no restriction
                on the dimension p and n.
                 Now, if k_0 is proportional to n, this means that the smallest non  zero element of the matrix Σ_1 - Σ_2 has  to satisfy
                
                    |σ_1(a,b)-σ_2(a,b)  |  >  C max(log p,log n)/√(n)
                
                for some constant C.  This means that the procedure estimates k_0 consistently  even if the differences between
                the elements of the two matrices are very small.
                
                
                § FINITE SAMPLE PROPERTIES
                
                In this section we investigate the finite sample properties of the new change point estimator by means of a simulation study
                and compare our approach with two alternative methods proposed by <cit.> and <cit.>, which
                are most similar in spirit as the procedure proposed in the present paper.
                
                To be precise let r_0=k_0/n be the “true”  change point fraction and  let  r̂=k̂/n, where  k̂
                is the new change point estimator  defined   in (<ref>). All the numerical results below are calculated from 200
                replications  and we  obtain the simulated estimates  r̂_1,…, r̂_K  of r_0.
                In the following discussion we present
                the mean
                
                    r̅=1/K∑_i=1^K r̂_i ,
                
                the standard deviation
                
                    std (r̂) =√(1/K-1∑_i=1^K (r̂_i-r̅)^2)
                
                and the corresponding mean squared error
                
                    MSE=1/K∑_i=1^K (r̂_i-r_0)^2 =K-1/K· std^2  (r̂) +(r̅-r_0)^2.
                
                Throughout this section we denote  by   blk(,)
                 a block-diagonal matrix composed by matrices  and  of appropriate dimension.
                 Σ_1=_p is always  the identity matrix and we consider  four different choices for the matrix
                 Σ_2 to investigate the performance of the new estimator under the following  alternatives
                
                  *  case 1:Σ_2=1.5∗_p;             case 2:Σ_2=2∗_p;
                
                  *  case 3:Σ_2= blk(4,_p-1);         case 4:Σ_2= blk(8,_p-1).
                
                Cases 1 and 2 indicate that there are many (p positions) small disturbances between Σ_1 and Σ_2, with a
                 magnitude increasing  from 0.5 to 1. On the other hand
                  there is only one disturbance between the two population covariance matrices in cases 3 and 4, but the magnitude is more significant (3 and 7 respectively).
                
                The true change point fraction is  chosen as  r_0=k_0/n=0.5 and  the first k_0 and  the last (n-k_0)   samples
                are generated from a  multivariate normal distribution N _p(0, Σ_1) and  a  N_p(0, Σ_2) distribution, respectively.
                
                
                
                
                 §.§ Performance of the new  estimator 
                
                
                
                In order to investigate the finite sample properties of the new estimator  we choose two sample sizes
                n=100 and n=200 and consider different  dimensions p ranging from 5 to 500. For each  pair (n,p), the mean change point fraction,
                standard deviation (std) and mean squared error (MSE) are recorded for all four cases under consideration,
                and the results are summarized in Table <ref> (n=100) and Table <ref> (n=200).
                 The numerical results from the two tables can be summarized as follows:
                
                 
                  (1)
                When the sample size n increases, the  performance of the estimator  is better.
                
                  (2) The dimension p of the data  does not have a significant influence on the performance of the estimators. In particular
                the mean squared error is remarkably stable with respect to the dimension in all four cases under consideration.
                
                  (3)   When the magnitude of the disturbance between Σ_1 and Σ_2 increases,
                 the estimator  performs   better (compare the results from  case 1 with case 2 or from case 3 with  case 4).
                
                
                Next we investigate the influence of   the dimension reduction step on performance
                of  the estimator. To this end, we consider  case 1 and case 3 with sample sizes  n=200
                and present   in Table <ref> the corresponding results without dimension reduction.
                In other words we apply the estimator (<ref>) based on all components.
                We note that the computation time without dimension reduction  is substantially larger
                because we work with p(p+1)/2-dimensional vectors.
                
                Comparing Table <ref> with the corresponding  results  in Table <ref>, we observe the
                following.
                 
                  (1)  In case 1, the  differences in the bias of  r̂ are  negligible (in both cases the mean is very close to 0.5).
                 On the other hand the standard deviations and as consequence the MSE in Table <ref>
                are smaller, which means that the   estimator without dimension reduction is more stable.
                Note that there are many small disturbances between two population matrices and thus keeping all positions
                promotes a safer estimation.
                On the other hand  the MSE results in Table <ref>  from the estimator using dimension reduction
                are already satisfactory.
                 
                  (2) In  case 3 when there is only one significant disturbance the situation is different. Although
                 the bias  of r̂ in Table <ref> is smaller, its
                 standard deviation and MSE  increase very  fast to an unacceptable level with increasing dimension.
                 This means that the estimator without dimension reduction is not reliable if the dimension is large.
                 
                  (3) As an interesting phenomenon we note  that in contrast to  Table <ref>
                 the standard deviation and mean squared error  in Table <ref> show a downward trend in case 1
                 but  an upward tendency in case 3. This observation can be explained by
                 the fact that without dimension reduction the gap between Σ_1 and Σ_2 increases
                 with the dimension p in case 1 but decreases in case 3.
                
                
                 §.§ Comparison with an estimator based on a quadratic form
                
                We first compare the new method with the estimator (2.12) suggested in <cit.>.
                Note that this statistic involves an inverse matrix Σ̂_n^-1, where Σ̂_n is an estimator of the long-run covariance satisfying their condition (2.6). Under our setting, this long-run covariance reduces to the population covariance matrix of a p(p+1)/2-dimensional random vector.
                As consequence the dimension p has to be substantially  smaller than the sample size n  to estimate the inverse
                of the covariance matrix precisely. In order to get a larger range for the dimension  p, we use  the sample size  n=400 in this subsection and let
                the dimension p vary from 5 to 50 (for larger values of p the method of <cit.> shows some instabilties).
                The location of the change point is assumed to be k_0=200.
                 In Table <ref> we display the results of the estimator  in <cit.> and the estimator  (<ref>)  proposed  in this paper, where we restrict to the case
                 2 and case  3 for the sake of brevity (the cases 1 and 4 show a very similar picture).
                We observe that the new estimator always performs better. While this superiority is only minor for small dimension, it becomes substantial for p=40, 50, in this case
                the mean squared error of the estimator  in  (2.12)  suggested by   <cit.> is  very large  (compared to the cases p ≤ 20), while
                the new estimate shows a remarkable stability with respect to different dimensions. The differences are also visualized  in Figure <ref>   for the
                cases   2 and 3, respectively, where we show the histograms of both estimates obtained from the different simulation runs.
                 The sample size is n=400, the  dimension is  p=25 and the  change point is located at  k_0=200 (red line).
                
                
                
                
                 §.§ Comparison with a multiscale estimator 
                
                
                We conclude this section with a brief comparison with the procedure in <cit.> who proposed a multiscale
                approach for the localization of the   change point.
                For the sake of comparison  we use  the same design as in  Section 5.1 of this paper.  We also  performed a comparison under scenarios considered
                in Section <ref>. Here the method proposed by these authors does not yield reliable estimates of the change point and
                the  results are not displayed for the sake of brevity.
                
                
                 To be precise  we summarize the setting  in <cit.>   here briefly.  The covariance matrix before  the change point is  Σ_1=_p
                   and the  matrix  Σ_2  after the change point is  generated as follows.  First a Poisson distributed random variable K∼ Poiss(3)
                   is generated. Then  the matrix Σ_2 is composed as a block-diagonal matrix
                of K (symmetric) matrices of size 2× 2 with ones on their diagonals and their off-diagonal element drawn uniformly from the set  [-0.6;-0.3]∪[0.3;0.6].
                The remaining  (p-2k)×(p-2k) diagonal block  of  the  matrix  Σ_2 is the identity matrix and all other elements of Σ_2 are 0.
                We consider a  sample of   n=1000 observations with  dimension p=50, where the “true”  change point
                is given by   k_0=n/2=500. The procedure of  <cit.> also requires the specification of a  set  ℐ_s  corresponding to observations without change points
                and we use ℐ_s=[1,2,⋯,100] (as suggested in their paper) in our simulation.  Note also that according to (2.5) in Section 2.3 of <cit.>
                a change point estimator is  only defined if there exists a narrowest window detecting a change-point. This was in 199 of  the 200 replications the
                 case.
                
                 In  Table <ref> we show the simulated mean, standard deviation and mean squared error of both estimates for 200 simulation runs.
                 We observe that the estimator proposed in this paper shows a substantially better performance than the multiscale estimator introduced by
                  <cit.>.  Histograms of the simulated change point positions for both methods in Figure <ref> point to the same conclusion.
                
                
                 Acknowledgements.
                The authors would like to thank
                 M. Stein who typed parts of this manuscript with considerable technical expertise.
                 We are also grateful to V.  Avanesov and N. Buzun for making the R-code
                 of their procedure available to us.
                The work of H. Dette   was partially supported by the Deutsche
                Forschungsgemeinschaft (DFG Research Unit 1735). Guangming Pan was supported in part by by a MOE Tier 2 grant 2014-T2-2-060 and by a MOE
                Tier 1 Grant RG25/14 at the Nanyang Technological University, Singapore.
                
                
                
                
                chicago
                
                
                § PROOF OF MAIN RESULTS
                
                
                We first introduce some auxiliary results that are   frequently used in the proofs. Lemma <ref> is a direct conclusion of Lemma 2.7.7 in <cit.>.
                
                For any 1≤ a,b≤ p, X_iaX_ib is a sub-exponential random variable, i.e. there exist positive constants C_1, C_2 > 0 (which do not depend on the indices i and a) such that for every t>0,
                
                    (|X_iaX_ib|>t)≤ C_1e^-C_2t.
                
                
                Let W_1,⋯,W_N be independent, mean zero, sub-exponential random variables. Then, there exist positive constants C_1, C_2 > 0 such that for every t>0,
                
                    {|1/N∑_i=1^N W_i| >t}≤ C_1e^-C_2N·min (t^2,t).
                
                Let W_1,⋯,W_N be independent, mean zero, sub-gaussian random variables. Then, there exist positive constants C_1, C_2 such that for every t>0,
                
                    {|1/N∑_i=1^N W_i| >t}≤ C_1e^-C_2Nt^2.
                
                
                 §.§ Proof of Theorem <ref>
                
                
                Observing   the construction of the statistic V_k in (<ref>) and (<ref>), we may assume without loss of generality that =0.  The components of the vector V_k and D corresponding to the entry in the position (a,b) of the matrices Σ_1,Σ_2, 1≤ a,b≤ p are given by
                
                    V_k(a,b)   =   1/k(k-1)∑∑_i≠ j≤ k(Ẋ_iaẊ_ib)(Ẋ_jaẊ_jb)
                    +1/(n-k)(n-k-1)∑∑_i≠ j>k(Ẋ_iaẊ_ib)(Ẋ_jaẊ_jb)      -2/k(n-k)∑_i≤ k∑_j>k(Ẋ_iaẊ_ib)(Ẋ_jaẊ_jb),
                    D(a,b)   =   1/n-3∑_k=2^n-2k(n-k)/nV_k(a,b) = D^(1)(a,b)+D^(2)(a,b)+D^(3)(a,b),
                
                where Ẋ_ia=X_ia-X̅_a=X_ia-1/n∑_j=1^n X_ja by (<ref>) and the terms D^(ℓ)(a,b) are defined by
                
                    D^(1)(a,b)   =   1/n-3∑_k=2^⌊√(n)⌋k(n-k)/nV_k(a,b), 
                    
                    D^(2)(a,b)   =   1/n-3∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1k(n-k)/nV_k(a,b),
                    
                    D^(3)(a,b)   =   1/n-3∑_k=n-⌊√(n)⌋^n-2k(n-k)/nV_k(a,b).
                
                 The reason
                for this decomposition of  D(a,b) is that all the k's or (n-k)'s  in D^(2)(a,b) are sufficiently large, while  both D^(1)(a,b) and D^(3)(a,b) only involve (⌊√(n)⌋-1) terms and thus the coefficient 1/n-3=⌊√(n)⌋-1/n-3·1/⌊√(n)⌋-1  gives us an extra factor of order 1/√(n) in the calculations.
                
                To be precise, let
                
                    X̅_k(a,b)=1/k∑_i=1^k X_iaX_ib,   Y̅_k(a,b)=1/n-k∑_i=k+1^n X_iaX_ib,       X̅_ka=1/k∑_i=1^k X_ia,      Y̅_ka=1/n-k∑_i=k+1^n X_ia,      Ẋ̅̇_k(a,b)=1/k∑_i=1^k Ẋ_iaẊ_ib,   Ẏ̅̇_k(a,b)=1/n-k∑_i=k+1^n Ẋ_iaẊ_ib
                
                and
                
                    V_k1   =k(n-k)n[Ẋ̅̇_k(a,b)-Ẏ̅̇_k(a,b)]^2,   V_k2=k(n-k)n[1k-1Ẋ̅̇_k(a,b)^2+1n-k-1Ẏ̅̇_k(a,b)^2],
                    
                    
                    V_k3   =k(n-k)n [1k(k-1)∑_i=1^k
                    (Ẋ_iaẊ_ib)^2+1(n-k)(n-k-1)∑_i=k+1^n
                    (Ẋ_iaẊ_ib)^2 ],
                
                then a straightforward but tedious calculation yields
                
                    Ẋ̅̇_k(a,b)   =X̅_k(a,b)-X̅_kaX̅_b-X̅_kbX̅_a+X̅_aX̅_b, 
                    Ẏ̅̇_k(a,b)   =Y̅_k(a,b)-Y̅_kaX̅_b-Y̅_kbX̅_a+X̅_aX̅_b
                
                and
                
                    k(n-k)/nV_k(a,b)=V_k1+V_k2-V_k3.
                
                With these notations we decompose the quantities D^(i)(a,b) as follows:
                
                    D^(i)(a,b)≜ A^(i)+B^(i)-C^(i),
                
                 where
                
                    A^(1)=1/n-3∑_k=2^⌊√(n)⌋V_k1,     B^(1)=1/n-3∑_k=2^⌊√(n)⌋V_k2,     C^(1)=1/n-3∑_k=2^⌊√(n)⌋V_k3,
                       A^(2)=1/n-3∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1V_k1,   B^(2)=1/n-3∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1V_k2,   C^(2)=1/n-3∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1V_k3, 
                       A^(3)=1/n-3∑_k=n-⌊√(n)⌋^n-2V_k1,   B^(3)=1/n-3∑_k=n-⌊√(n)⌋^n-2V_k2,   C^(3)=1/n-3∑_k=n-⌊√(n)⌋^n-2V_k3.
                
                
                
                  Without loss of generality, we can assume that σ_1(a,b)=σ_2(a,b) =0.
                 Observing the decomposition (<ref>) the assertion (<ref>)
                follows from
                
                    {⋃_(a,b)∈𝒩{D^(i)(a,b)>cτ}}→ 0,    i=1,2,3.
                
                When we prove these results we derive exponential inequalities for all three probabilities, which directly yield the estimate in (<ref>).
                
                In (<ref>) and hereafter in the proof, c and c_i(i=1,2,⋯) indicate some positive constants that may change from line to line.
                According to the decomposition (<ref>), it is sufficient to derive exponential inequalities, which will be used to
                verify the following results:
                
                    p^2·{A^(i)>cτ}→ 0,   p^2·{B^(i)>cτ}→ 0,   p^2·{C^(i)>cτ}→ 0,   i=1,2,3.
                 The case i=1,3. For the index i=1 and i=3 the arguments are very similar and for the sake of brevity we only
                 consider the case   i=1 in (<ref>). For the statistic A^(1) we find that
                
                    {A^(1)>cτ}
                       =   {1√(n)-1∑_k=2^⌊√(n)⌋√(n)-1n-3k(n-k)n[Ẋ̅̇_k(a,b)-Ẏ̅̇_k(a,b)]^2>cτ}   ≤   ∑_k=2^⌊√(n)⌋{√(n)-1n-3k(n-k)n[Ẋ̅̇_k(a,b)-Ẏ̅̇_k(a,b)]^2>cτ}   ≤   ∑_k=2^⌊√(n)⌋{|Ẋ̅̇_k(a,b)-Ẏ̅̇_k(a,b)|>√(cτ M_n,k,α)},
                
                where we use the notation M_n,k,α = √(n)· n/k(n-k) for the sake of a transparent notation.
                Considering the components that constitute Ẋ̅̇_k(a,b) and Ẏ̅̇_k(a,b) in (<ref>), we use  Lemma <ref> and <ref> to calculate the following probabilities to get an upper bound of (<ref>)
                
                    {|X̅_k(a,b)|>√(cτ M_n,k,α)}   ≤    c_1e^-c_2kmin(τ M_n,k,α,√(τ M_n,k,α) )   ≤    c_1e^-c_2min(τ√(n),√(τ√(n))),{|Y̅_k(a,b)|>√(cτ M_n,k,α)}   ≤    c_1e^-c_2(n-k)min(τ M_n,k,α,√(τ M_n,k,α) )   ≤   c_1 e^-c_2min(τ n, √(τ n(n-√(n))))
                
                (k=2,…,  ⌊√(n)⌋). Similarly, it follows from Lemma <ref>
                    {|X̅_kaX̅_b|>√(cτ M_n,k,α)}
                       ≤   {|X̅_ka|> c_1(τ M_n,k,α)^1/4}+{|X̅_b|> c_2(τ M_n,k,α)^1/4}
                       ≤   c_3[e^-c_4k√(τ M_n,k,α)+e^-c_5n√(τ M_n,k,α)]≤ c_3[e^-c_4√(τ√(n))+e^-c_5√(τ)n],
                
                for k=2,…,  ⌊√(n)⌋, and using similar arguments we obtain for
                
                    {|Y̅_kaX̅_b|>√(cτ M_n,k,α)}   ≤   
                     c_3  [e^-c_4√(τ n(n-√(n)))+e^-c_5√(τ)n ], 
                    {|X̅_aX̅_b|>√(cτ M_n,k,α)}   ≤    c_3e^-c_4n√(τ M_n,k,α)≤  c_3e^-c_4√(τ)n,
                k=2,…⌊√(n)⌋. The same estimates (with different constants c,c_1,c_2,…) can be derived for the remaining terms involving X̅_kbX̅_a and Y̅_kbX̅_a.
                Combining (<ref>), (<ref>)-(<ref>) we obtain the upper bound
                
                    {A^(1)>cτ}   ≤    c_1√(n)(e^-c_2min(τ√(n),√(τ√(n)))+e^-c_2min(τ n, √(τ n(n-√(n))))      +e^-c_2√(τ√(n))+e^-c_2√(τ n (n-√(n)))+e^-c_2√(τ)n).
                
                Since the smallest absolute value among the exponents in (<ref>) is √(τ√(n))
                we have
                
                    {A^(1)>cτ}≤
                     c_3 ·√(n)e^-c√(τ)n^1/4
                
                 for some small positive constant c. Consequently, using the assumption (<ref>)   it follows that p^2·{A^(1)>cτ}→ 0.
                
                
                Using similar arguments, we can investigate the other terms.
                To be precise consider the statistic  B^(1), for which we obtain the estimate
                
                    {B^(1)>cτ}   =   {1n-3∑_k=2^⌊√(n)⌋k(n-k)n[1k-1Ẋ̅̇_k(a,b)^2+1n-k-1Ẏ̅̇_k(a,b)^2]>cτ}
                       ≤   ∑_k=2^⌊√(n)⌋{√(n)-1n-3k(n-k)n[1k-1Ẋ̅̇_k(a,b)^2+1n-k-1Ẏ̅̇_k(a,b)^2]>cτ}
                    
                       ≤   ∑_k=2^⌊√(n)⌋ [
                    { |Ẋ̅̇_k(a,b) |>√(cτ M^(1)_n,k,α)}+
                    { |Ẏ̅̇_k(a,b) |>√(cτ M^(2)_n,k,α)} ] ,
                
                where the quantities M^(1)_n,k,α and M^(2)_n,k,α are defined by
                
                    M^(1)_n,k,α = √(n)· n/(n-k)    M^(2)_n,k,α=  √(n)· n/k,
                
                respectively. Again we  consider  the components that constitute Ẋ̅̇_k(a,b) and Ẏ̅̇_k(a,b) in (<ref>)
                separately and obtain by an application of   Lemma <ref> and <ref> the following estimates
                (for k=2 , … ,   ⌊√(n)⌋)
                
                    {|X̅_k(a,b)|>√(cτ M^(1)_n,k,α)}   ≤   
                    c_1e^-c_2kmin(τ M^(1)_n,k,α,√(τ M^(1)_n,k,α) )≤  c_1 e^-c_2min(τ√(n),√(τ√(n))),{|Y̅_k(a,b)|>√(cτ M^(2)_n,k,α)}   ≤   
                    c_1e^-c_2(n-k)min(τ M^(2)_n,k,α,√(τ M^(2)_n,k,α) )   ≤    c_1e^-c_2min(τ n(n-√(n)), √(τ n(n-√(n))^2)).
                
                Similarly,  Lemma <ref> gives for   k=2 , … , ⌊√(n)⌋
                    {|X̅_kaX̅_b|>√(cτ M^(1)_n,k,α)}   ≤     c_3[e^-c_4k√(τ M^(1)_n,k,α)+e^-c_5n√(τ M^(1)_n,k,α)]
                       ≤    c_3   [e^-c_4√(τ√(n))+e^-c_5n√(τ√(n))],{|Y̅_kaX̅_b|>√(cτ M^(2)_n,k,α)}   ≤   
                    c_3[e^-c_4(n-k)√(τ M^(2)_n,k,α)+e^-c_5n√(τ M^(2)_n,k,α)] 
                       ≤   
                    c_3[e^-c_4(n-√(n))√(τ n)+e^-c_5n√(τ n)],{|X̅_aX̅_b|>√(cτ M^(1)_n,k,α)}   ≤    c_3e^-c_4n√(τ M^(1)_n,k,α)≤ c_3e^-c_4n√(τ√(n)),{|X̅_aX̅_b|>√(cτ M^(2)_n,k,α)}   ≤    c_3e^-c_4n√(τ M^(2)_n,k,α)≤ c_3e^-c_4n√(τ n).
                
                Thus, observing that the   smallest absolute value among the exponents in these estimates  is  given by √(τ√(n)),
                 an upper bound for the probability in (<ref>) is obtained as
                
                    {B^(1)>cτ}≤ c_1  √(n) e^-c√(τ) n^1/4
                
                 for some small positive constant c. Consequently, observing Assumption <ref> we have p^2{B^(1)>cτ}→ 0  as n,p →∞.
                
                Finally the term  C^(1) is estimated as follows
                
                    {C^(1)>cτ}   ≤   ∑_k=2^⌊√(n)⌋{√(n)-1n-3k(n-k)n [1k(k-1)∑_i=1^k
                    (Ẋ_iaẊ_ib)^2+1(n-k)(n-k-1)∑_i=k+1^n
                    (Ẋ_iaẊ_ib)^2 ]>cτ}   ≤   ∑_k=2^⌊√(n)⌋[
                    {1k∑_i=1^k(Ẋ_iaẊ_ib)^2>c_1τn√(n)n-k}+
                     {1n-k∑_i=k+1^n(Ẋ_iaẊ_ib)^2>c_2τ n√(n)k}]    ≤   ∑_k=2^⌊√(n)⌋[
                    k{|Ẋ_iaẊ_ib|>√(c_1τ√(n))}+ (n-k)
                    {|Ẋ_iaẊ_ib|>√(c_2τ n)}] .
                
                Observing that
                Ẋ_iaẊ_ib=X_iaX_ib-X_iaX̅_b-X_ibX̅_a+X̅_aX̅_b,
                we obtain from  Lemma <ref> and  <ref> the estimates
                
                    k{|X_iaX_ib|>√(c_1τ√(n))}   ≤   
                    c √(n)e^-c_4√(τ√(n)),
                     (n-k){|X_iaX_ib|>√(c_2τ n)}   ≤   
                    cn e^-c_6√(τ n), 
                    k{|X_iaX̅_b|>√(c_1τ√(n))}   ≤   
                     c√(n)[e^-c_7√(τ√(n))+e^-c_8n√(τ√(n))],
                     (n-k){|X_iaX̅_b|>√(c_2τ n)}   ≤     cn[e^-c_7√(τ n)+e^-c_8n√(τ n)], 
                      k{|X̅_aX̅_b|>√(c_1τ√(n))}   ≤    c √(n)e^-c_7n√(τ√(n)),
                     (n-k){|X̅_aX̅_b|>√(c_2τ n)}   ≤    c n e^-c_7n√(τ n),
                
                whenever k=2, …  , ⌊√(n)⌋.
                Summarizing we have
                
                    {C^(1)>cτ}   ≤    c_1 [n(e^-c_2√(τ√(n))+e^-c_2n√(τ√(n)))+
                    n√(n)(e^-c_2√(τ n)+e^-c_2n√(τ n))]   ≤   c_3ne^-c√(τ√(n))
                
                for some small positive constant c. Now assumption (<ref>)  implies   p^2{C^(1)>cτ}→ 0, if
                
                    p^2ne^-c √(τ)n^1/4→ 0,
                
                 and therefore   the proof
                of (<ref>)  in the case i=1 is completed (the case i=3 follows by exactly the same arguments).
                
                 The case i=2.
                For the term A^(2), we get that
                
                    {A^(2)>cτ}   ≤   ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1 [
                    {|Ẋ̅̇_k(a,b)|>√(c_1τ M_n,k)}+ { |Ẏ̅̇_k(a,b) |>√(c_2τ M_n,k)}] ,
                
                where we use the notation  M_n,k= n/k(n-k). Similar calculations as given  for the term  A^(1) give for the  summands  in
                the representation (<ref>)  of  Ẋ̅̇_k(a,b) and Ẏ̅̇_k(a,b) the estimates (here we use the fact that for i=2 we have k=⌊√(n)⌋ +1, …, n-⌊√(n)⌋-1)
                
                    {|X̅_k(a,b)|>√(cτ M_n,k)}   ≤      c_1e^-c_2kmin(τ M_n,k,√(τ M_n,k) )≤
                    c_1  e^-c_2min(τ ,√(τ√(n))),{|Y̅_k(a,b)|>√(cτ M_n,k)}   ≤     c_1e^-c_2(n-k)min(τ M_n,k,√(τ M_n,k) )≤
                    c_1  e^-c_2min(τ , √(τ√(n))) ,{|X̅_kaX̅_b|>√(cτ M_n,k)}   ≤    c_3[e^-c_4k√(τ M_n,k)+e^-c_5n√(τ M_n,k)]≤ c_3[e^-c_4√(τ√(n))+e^-c_5√(τ n)], {|Y̅_kaX̅_b|>√(cτ M_n,k)}   ≤    c_3[e^-c_4(n-k)√(τ M_n,k)+e^-c_5n√(τ M_n,k)]≤ c_3[e^-c_4√(τ√(n))+e^-c_5√(τ n)],{|X̅_aX̅_b|>√(cτ M_n,k)}   ≤    c_3e^-c_4n√(τ M_n,k)≤  c_3 e^-c_4√(τ n),
                
                and consequently we obtain
                
                    p^2·{A^(2)>cτ}≤  c_1  n p^2 max{ e^-cτ ,  e^-c√(τ√(n))},
                
                 for some small positive constant c, which converges to 0 under the stated assumptions (<ref>) and (<ref>).
                
                For the term B^(2) we have
                
                    {B^(2)>cτ}   ≤   ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1 [ {  |Ẋ̅̇_k(a,b) |>√(c_1τnn-k)}+
                    { |Ẏ̅̇_k(a,b) |>√(c_2τnk)}] ,
                
                where the two probabilities can be bounded    taking into account the representation  (<ref>) and the
                estimates
                
                    {  |X̅_k(a,b) |>√(c_1τ nn-k)}≤
                    c_1e^-c_2kmin(τ n/n-k,√(τ n/n-k))≤ c_1 e^-c_2min(τ√(n),√(τ n)),      {|Y̅_k(a,b)|>√(cτ nk)}≤
                    c_1e^-c_2(n-k)min(τ n/k,√(τ n/k))≤ c_1e^-c_2min(τ√(n),√(τ n)),       {|X̅_kaX̅_b|>√(cτ nn-k)}≤
                    c_3[e^-c_4k√(τ n/n-k)+e^-c_5n√(τ n/n-k)]≤ c_3
                    [e^-c_4√(τ n)+e^-c_5n√(τ)],      {|Y̅_kaX̅_b|>√(cτ nk)}≤ c_3[e^-c_4(n-k)√(τ n/k)+e^-c_5n√(τ n/k)]≤ c_3 [e^-c_4√(τ n)+e^-c_5n√(τ)],      {|X̅_aX̅_b|>√(cτ nn-k)}≤  c_3e^-c_4n√(τ n/n-k)≤ c_3  e^-c_4n√(τ),      {|X̅_aX̅_b|>√(cτ nk)}≤ c_3 e^-c_4n√(τ)
                
                for k=⌊√(n)⌋+1 , ⌊√(n)⌋+2 , … , n-⌊√(n)⌋-1.
                Therefore, we obtain as an  upper bound for the probability in  (<ref>)
                
                    p^2·{B^(2)>cτ}   ≤    c_1p^2n(e^-c_2min(τ√(n),√(τ n))+e^-c_2√(τ n)+e^-c_2n√(τ))   =   c_1p^2n(e^-c_2√(n)√(τ)+e^-c_2√(τ n)+e^-c_2n√(τ))= o(1),
                
                 by assumption (<ref>).
                
                
                
                
                For the term C^(2), we observe (note that n/(n-k)>1 and n/k>1)
                
                    {C^(2)>cτ}   ≤   ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1 [ {1k∑_i=1^k(Ẋ_iaẊ_ib)^2>c_1τnn-k}+ {1n-k∑_i=k+1^n(Ẋ_iaẊ_ib)^2>c_2τnk} ]   ≤   ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1k{|Ẋ_iaẊ_ib|>√(c_1τ)}
                        +    ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1(n-k){|Ẋ_iaẊ_ib|>√(c_2τ)},
                
                where the two probabilities can be bounded using the representation
                Ẋ_iaẊ_ib=X_iaX_ib-X_iaX̅_b-X_ibX̅_a+X̅_aX̅_b. This gives for k= ⌊√(n)⌋+1,  … , n-⌊√(n)⌋-1 the estimates
                
                    {|X_iaX_ib|>c√(τ)}   ≤    c_1 e^-c_2√(τ),{|X_iaX̅_b|>c√(τ)}   ≤   [{|X_ia|>c_1τ^1/4}+{|X̅_b|>c_2τ^1/4}]
                    ≤ c_3 [e^-c_4√(τ)+e^-c_5n√(τ)],{|X̅_aX̅_b|>c√(τ)}   ≤   [{|X̅_a|>c_1τ^1/4}+{|X̅_b|>c_2τ^1/4}]≤ c_3 e^-c_4n√(τ).
                
                Observing that ∑^n-⌊√(n)⌋-1_k=⌊√(n)⌋+1 k ≤ n^2, ∑_k=⌊√(n)⌋+1^n-⌊√(n)⌋-1(n-k)≤ n^2  we can bound  the probability on the left hand side of (<ref>) by
                
                    {C^(2)>cτ}≤ c_1 n^2(e^-c_2√(τ)+e^-c_3n√(τ)),
                
                which vanishes asymptotically by assumption (<ref>).
                
                Therefore, we have established (<ref>) for all indices i = 1,2,3 which implies (<ref>), and the assertion of Theorem <ref> follows. A careful inspection of our arguments shows that we have also established the estimate (<ref>) in Theorem <ref>.
                
                
                
                 §.§ Proof of Theorem <ref>
                
                
                Recall that  the set 𝒫 is the set of indices corresponding to elements with
                |σ_1(a,b)-σ_2(a,b) |>λ, where   λ satisfies (<ref>), and  note that the assertion (<ref>) is equivalent to
                
                    {⋃_(a,b)∈𝒫{D(a,b)≤τ}}→ 0.
                
                For a proof of this statement we derive several  exponential inequalities which directly yield the estimate in (<ref>).
                For this purpose we  introduce the decomposition
                
                    D(a,b)=1/n-3∑_k=2^n-2[V_k1+V_k2-V_k3] =  A+B-C,
                
                 where the quantities A, B and C are given by
                
                    A=1/n-3∑_k=2^n-2V_k1,   B=1/n-3∑_k=2^n-2V_k2,   C=1/n-3∑_k=2^n-2V_k3,
                
                respectively, and the statistics V_k1, V_k2 and V_k_3 have been defined in (<ref>).
                Observing the inclusion
                 
                    {D(a,b)≤τ} ⊂{A≤3τ}∪{|B-C|≥ 2τ}⊂{A≤3τ}∪{B≥τ}∪{C≥τ}
                
                  the assertion of Theorem <ref> follows  by deriving exponential inequalities for the probabilities for these three events, which are used to prove
                
                    p^2·{A≤3τ}→ 0,   p^2·{B≥τ}→ 0,   p^2·{C≥τ}→ 0.
                
                We now investigate the three probabilities in (<ref>) separately. First, for the last two terms B and C,
                from (<ref>), it is easy to see that B=B^(1)+B^(2)+B^(3) and C=C^(1)+C^(2)+C^(3), where - as in the proof of Theorem <ref> - the indices (1), (2) and (3) correspond to the summation with respect to the sets { 2,…,⌊√(n)⌋}, {⌊√(n)⌋ +1, …, n-⌊√(n)⌋ -1} and {n-⌊√(n)⌋,…,n-2  }, respectively.
                Moreover,  using similar arguments as given in the proof of Theorem <ref>, the following results can be established.
                
                
                
                
                    p^2·{B^(i)>cτ}   ≤   c_1  p^2√(n) e^-c√(τ) n^1/4, i=1,3,   p^2·{B^(2)>cτ}≤ c_1p^2ne^-c_2√(τ n),
                     p^2·{C^(i)>cτ}   ≤   c_3p^2ne^-c√(τ)n^1/4, i=1,3,  
                    p^2·{C^(2)>cτ}≤ c_1p^2 n^2e^-c_2√(τ).
                
                The only difference when analyzing the above probabilities lies in the expectation of X_iaX_ib, (X_iaX_ib), which is not necessarily zero now. But this does not affect the proof  of (<ref>) since (X_iaX_ib)  is always bounded by Assumption <ref>. Hence
                
                    p^2·{B≥τ}→ 0,   p^2·{C≥τ}→ 0.
                
                In order to show the remaining exponential equation and prove the  assertion  p^2·{A≤3τ}→ 0, we note that σ_1(a,b)≠σ_2(a,b), and a straightforward calculation gives
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                    (X̅_k(a,b)-Y̅_k(a,b))={[ n-k_0/n-k(σ_1(a,b)-σ_2(a,b))                    : k≤ k_0,;    k_0/k(σ_1(a,b)-σ_2(a,b) )                     : k>k_0. ].
                
                Let
                
                    A_k=Ẋ̅̇_k(a,b)-Ẏ̅̇_k(a,b)-(X̅_k(a,b)-Y̅_k(a,b)),   B_k=(X̅_k(a,b)-Y̅_k(a,b)),
                
                and define b_n,k≜k(n-k)/n. Then one can observe that
                
                    |B_k|≥min{n-k_0n,k_0n}·|σ_1(a,b)-σ_2(a,b)|
                
                and
                
                    1n-3∑_k=2^n-2b_n,kB_k^2    ≥   min{(n-k_0)^2n^2,k_0^2n^2}·  n(σ_1(a,b)-σ_2(a,b) )^2·1n-3∑_k=2^n-2k(n-k)n^2   ≥   min{(n-k_0)^2n^2,k_0^2n^2}· n(σ_1(a,b)-σ_2(a,b) )^26≥min{(n-k_0)^2n^2,k_0^2n^2}· nλ^26.
                
                Consequently, by assumption (<ref>) (with a sufficiently large constant C) we  obtain the estimate
                
                    3τ≤1/2(n-3)∑_k=2^n-2b_n,kB_k^2.
                
                On the other hand, we have by the definition of A, A_k and B_k
                    {A≤3τ}
                       =   {1n-3∑_k=2^n-2b_n,k[A_k+B_k]^2≤ 3τ}
                    
                       ≤   {2n-3∑_k=2^n-2b_n,kA_kB_k≤ 3τ-1n-3∑_k=2^n-2b_n,kB_k^2 }   ≤   {2n-3∑_k=2^n-2b_n,kA_kB_k≤ -12(n-3)∑_k=2^n-2b_n,kB_k^2 }   ≤   {2n-3|∑_k=2^n-2b_n,kA_kB_k|≥12(n-3)|∑_k=2^n-2b_n,kB_k^2| }   ≤    n·{|b_n,kA_k|≥c|B_k|min{(n-k_0)^2n^2,k_0^2n^2}· n(σ_1(a,b)-σ_2(a,b))^2 }   =    n·{|A_k|≥ cϕ},
                
                where ϕ=1/|B_k|1/b_n,kmin{(n-k_0)^2/n^2,k_0^2/n^2}· n(σ_1(a,b)-σ_2(a,b) )^2. In order to investigate the probability {|A_k|≥ cϕ} in (<ref>)   note that 1/|B_k|≥1/|σ_1(a,b)-σ_2(a,b) | by (<ref>). If
                
                    |σ_1(a,b)-σ_2(a,b) |≥λ≥ C√(τn)max{n^2(n-k_0)^2,n^2k_0^2},
                
                 it is therefore easy to see that
                
                    ϕ≥nk(n-k)min{(n-k_0)^2n^2,k_0^2n^2}· n|σ_1(a,b)-σ_2(a,b) |≥
                    C nk(n-k)√(nτ).
                
                Observing the decomposition (<ref>), the term A_k can be written as
                
                    A_k
                       =   [X̅_k(a,b)-Y̅_k(a,b)-(X̅_k(a,b)-Y̅_k(a,b))]-X̅_kaX̅_b-X̅_kbX̅_a+Y̅_kaX̅_b+Y̅_kbX̅_a,
                
                and by Lemma <ref> and Lemma <ref> we obtain
                
                    {|A_k|≥ cϕ}   ≤   {|X̅_k(a,b)-X̅_k(a,b)|≥ c_1ϕ}
                    +{|Y̅_k(a,b)-Y̅_k(a,b)|≥ c_2ϕ}      +2{|X̅_ka|≥ c_3√(ϕ)}+2{|X̅_b|≥ c_4√(ϕ)}+2{|Y̅_ka|≥ c_3√(ϕ)}+2{|X̅_b|≥ c_4√(ϕ)}   ≤   c_5[e^-c_6kmin{ϕ^2,ϕ}+e^-c_6(n-k)min{ϕ^2,ϕ}+e^-c_6kϕ+e^-c_6(n-k)ϕ+e^-c_6nϕ]   ≤   c_5[e^-c_6min{τ,√(nτ)}+e^-c_6min{τ,√(nτ)}+e^-c_6√(nτ)+e^-c_6√(nτ)+e^-c_6√(nτ)] = o(1),
                
                where the last estimate follows from (<ref>) and (<ref>).
                Combining (<ref>) with (<ref>), it follows that p^2·{A≤3τ}→ 0, which completes the proof of Theorem <ref>.
                
                
                
                
                
                 §.§ Proof of Theorem <ref>
                
                Recall the definition of the statistic U_n(k) in (<ref>)
                  and the definition of the change point estimator k̂ in (<ref>).
                
                
                
                Let _i indicate the  m-dimensional subvector of (_i_i^T) corresponding to the components   in the set D_τ in (<ref>). Obviously,
                
                    {|k̂k_0-1|≥ϵ}≤{k̂≥(1+ϵ)k_0}+{k̂≤ (1-ϵ)k_0},
                
                and we will derive exponential bounds for the  two terms on the right-hand side
                to prove that the probability vanishes asymptotically. We only consider the first term because the
                second term can be handled similarly. It is sufficient to show that
                
                    P{⋃_k≥(1+ϵ)k_0(U_n(k)≥ U_n(k_0))}≤∑_k≥(1+ϵ)k_0^n P{U_n(k)≥ U_n(k_0)}≤ c_1p^2n[e^-c_2τ+ne^-c_2√(τ)],
                
                 which follows if the estimate
                
                    P{U_n(k)≥ U_n(k_0)}≤ c_1p^2[e^-c_2τ+ne^-c_2√(τ) ]
                
                holds. By Assumption <ref> this term is of order o(1/n)
                uniformly with respect to   (1 + ε)k_0 ≤ k ≤ n.
                For a proof of this statement define the vectors _i=(W_i1,⋯,W_im) by
                
                    _i=_i+_i,
                
                and denote by _1 and _2 the m-dimensional vectors containing the elements   of the matrices Σ_1 and Σ_2, respectively, which correspond to positions (a,b) ∈𝒟_τ identified in Step 1 of the procedure.
                We will make use of the decomposition
                
                    U_n(k)   =   1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1[(_i-_j)+(_i-_j)]^T
                    [(_t-_l)+(_t-_l)]
                    
                       =   A(k)+B(k)+C(k)+E(k),
                
                where
                
                    A(k)   =   1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1(_i-_j)^T(_t-_l), 
                    
                    B(k)   =   1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1(_t-_l)^T(_i-_j), 
                    
                    C(k)   =   1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1(_i-_j)^T(_t-_l), 
                    
                    E(k)   =   1/n^4∑^k∑^k_(i≠ t)=1∑^n∑^n_(j≠ l)=k+1(_i-_j)^T(_t-_l),
                
                and begin investigating the constant terms E(k) and E(k_0). For this purpose we make use of the notation
                
                    a_n,k=  (n-k)(n-k-1)
                
                and obtain by a direct calculation
                
                    E =E(k_0)-E(k)=1/n^4[k_0(k_0-1)a_n,k_0-k_0(k_0-1)a_n,k]_1-_2^2   =   1/n^4k_0(k_0-1)[(n-k)(k-k_0)+(n-k-1)(k-k_0)+(k-k_0)^2]_1-_2^2.
                
                Observing the inclusion
                
                    { U_n(k) ≥ U_n(k_0)}⊂{ A(k) + B(k) + C(k) -( A(k_0) + B(k_0) + C(k_0) )≥ E }
                
                we now
                  investigate the other terms A(k) - A(k_0), B(k) - B(k_0), C(k) - C(k_0). A straightforward but tedious calculation yields the decomposition
                
                    B(k)-B(k_0)
                    
                    
                    
                    
                    = -B_1+B_2+B_3,
                
                where
                
                    B_1=1/n^4(k_0-1)[a_n,k_0-a_n,k](_1-_2)^T∑_i=1^k_0_i,      B_2=1/n^4k_0(k_0-1)(k-k_0)(_1-_2)^T∑_j=k_0+1^n_j,      B_3=1/n^4k_0(k_0-1)(n-k-1)(_1-_2)^T∑_j=k_0+1^k_j.
                
                Using (<ref>)  the term B_1 can be handled as follows
                
                    {|B_1|≥ cE}   =   {|(_1-_2)^T∑_i=1^k_0_i|≥ c k_0_1-_2^2}   ≤   
                    
                    {∑_i=1^k_0_i ^2 
                    ≥ (c k_0_1-_2)^2  }   ≤   
                     ∑_j=1^m{|∑_i=1^k_0W_ij|≥c k_0_1-_2√(m)}   =   ∑_j=1^m{|1k_0∑_i=1^k_0W_ij|≥ c_1-_2√(m)}.
                
                Note  that for each j=1,⋯, m, there exists a position (a,b) such that W_ij can be written as
                
                    W_ij=Ẋ_iaẊ_ib-(X_iaX_ib)=X_iaX_ib-(X_iaX_ib)-X_iaX̅_b-X_ibX̅_a+X̅_aX̅_b,
                
                and thus we obtain
                
                    1/k_0∑_i=1^k_0W_ij=1/k_0∑_i=1^k_0[X_iaX_ib-(X_iaX_ib)]-X̅_k_0aX̅_b-X̅_k_0bX̅_a+X̅_aX̅_b.
                
                With these notations the probability in (<ref>) can be further bounded using Lemmas <ref> and <ref>, that is
                
                    {|B_1|≥ cE}   ≤    m[{|1k_0∑_i=1^k_0[X_iaX_ib-(X_iaX_ib)]|≥ c_1_1-_2√(m)}+2{|X̅_k_0a|≥ c_2√(_1-_2√(m))}      +3{|X̅_b|≥ c_3√(_1-_2√(m))}+{|X̅_a|≥ c_4√(_1-_2√(m))}]   ≤   c_6m[e^-c_7min{k_0_1-_2^2/m, k_0_1-_2/√(m)}+
                    e^-c_7k_0_1-_2/√(m)+e^-c_7n_1-_2/√(m)]   ≤   c_6p^2[e^-c_7min{τ,√(nτ)}+e^-c_7√(nτ)+e^-c_7√(nτ)]≤ c_8p^2e^-c_7min{τ,√(nτ)},
                
                where the last inequality follows from the fact
                  that k_0_1-_2^2/m≥ Cτ and k_0_1-_2/√(m)≥ C√(nτ) if the smallest nonzero entry of (Σ_1-Σ_2) satisfies |σ_1(a,b)-σ_2(a,b) |≥λ≥ C√(τ/n)n/k_0. Note that Assumption <ref> gives p^2n=o(e^cτ), p^2n=o(e^c√(nτ)), and therefore it follows that
                
                    {|B_1|≥ cE}=o(1/n).
                
                Similarly, each component of the vector B_2 can be represented as
                
                    1/n-k_0∑_i=k_0+1^nW_ij=1/n-k_0∑_i=k_0+1^n[X_iaX_ib-(X_iaX_ib)]-Y̅_k_0aX̅_b-Y̅_k_0bX̅_a+X̅_aX̅_b
                
                for some (a,b) ∈𝒫.
                We can find that
                
                    {|B_2|≥ cE}   =   {|(_1-_2)^T∑_i=k_0+1^n_i|≥ c(n-k_0+n-k-1)_1-_2^2}   ≤   {_1-_2∑_i=k_0+1^n_i≥ c (n-k_0+n-k-1)_1-_2^2  }   ≤   m·{|∑_i=k_0+1^nW_ij|≥c(n-k_0+n-k-1)_1-_2√(m)}   ≤   m·{|1/n-k_0∑_i=k_0+1^nW_ij|≥c_1-_2√(m)}   ≤    c_1m [e^-c_2min{(n-k_0)_1-_2^2/m,(n-k_0)_1-_2/√(m)}+e^-c_2(n-k_0)_1-_2/√(m)+e^-c_2n_1-_2/√(m)]   ≤     c_1p^2[e^-c_3min{τ,√(nτ)}+e^-c_3√(nτ)+e^-c_3√(nτ)]≤ c_4p^2e^-c_3min{τ,√(nτ)},
                
                where the last inequality follows from the fact that    |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)n/n-k_0 for all (a,b) ∈𝒫, which implies
                (n-k_0)_1-_2^2/m≥ Cτ and (n-k_0)_1-_2/√(m)≥ C√(nτ).
                So we get
                
                    {|B_2|≥ cE}=o(1/n).
                
                In order to get a similar result for the term  B_3 we note that
                
                    {|B_3|≥ cE}
                       ≤   ℙ{_1-_2^2
                    
                    ∑_j=1^m(∑_i=k_0+1^kW_ij)^2≥ c [2(k-k_0)+(k-k_0)^2/n-k-1]^2_1-_2^4  }   ≤   m ℙ{|∑_i=k_0+1^kW_ij|≥c[2(k-k_0)+(k-k_0)^2n-k-1]_1-_2 /√(m)}   ≤   m·ℙ{|1k-k_0∑_i=k_0+1^kW_ij|≥2c_1-_2√(m)}   ≤    c_1m [e^-c_2min{(k-k_0)_1-_2^2/m,(k-k_0)_1-_2/√(m)}+e^-c_2(k-k_0)_1-_2/√(m)+e^-c_2n_1-_2/√(m)]   ≤    c_1m [e^-c_2ϵmin{k_0_1-_2^2/m,k_0_1-_2/√(m)}+e^-c_2ϵ k_0_1-_2/√(m)+e^-c_2n_1-_2/√(m)]    ≤   
                    c_3p^2e^-c_4min{τ,√(nτ)}.
                
                 By condition (<ref>) in Assumption <ref> and Assumption <ref>, we get
                
                    {|B_3|≥ cE}=o(1/n).
                
                
                
                Combining (<ref>)  and (<ref>), we can conclude that
                
                    {|B(k)-B(k_0)|≥ cE}≤  c_1p^2e^-c_2min{τ,√(nτ)}  =o(1/n)
                
                uniformly with respect to k ≥ k_0(1+ ε).
                Similarly, one can see that the estimate
                
                    {|C(k)-C(k_0)|≥ cE}≤  c_1p^2e^-c_2min{τ,√(nτ)} =o(1/n).
                
                Finally, we investigate the terms A(k) and A(k_0) introducing the decomposition
                
                    A(k)=A_1(k)-2A_2(k)+A_3(k),
                
                where
                
                    A_1(k)=1/n^4(n-k)(n-k-1)∑^k∑^k_(i≠ t)=1_i^T_t,  
                    A_3(k)=1/n^4k(k-1)∑^n∑^n_(j≠ l)=k+1_j^T_l,      A_2(k)=1/n^4(k-1)(n-k-1)∑_i=1^k∑_l=k+1^n_i^T_l.
                
                Then
                
                    A(k)-A(k_0)=[A_1(k)-A_1(k_0)]-2[A_2(k)-A_2(k_0)]+[A_3(k)-A_3(k_0)].
                
                In order to get that {|A(k)-A(k_0)|≥ cE}=o(1/n), it is sufficient to show that {|A_i(k)-A_i(k_0)|≥ cE}=o(1/n), i=1,2,3.
                In the following we only show that
                
                    {|A_1(k)-A_1(k_0)|≥ cE}≤ c_1p^2[e^-c_2τ+ne^-c_2√(τ) ]
                    =o(1/n)
                 uniformly with respect to k≥ (1 + ε)k_0. The other two terms can be treated similarly.
                Define
                
                    H_1=a_n,k_0-a_n,k/n^4∑^k_0∑^k_0_i≠ t_i^T_t,  
                    H_2=a_n,k/n^4∑^k∑^k_i≠ t=(k_0+1)_i^T_t,   H_3=a_n,k/n^4∑_i=1^k_0∑_t=k_0+1^k_i^T_t,
                
                and note that
                
                    A_1(k)-A_1(k_0)=-H_1+H_2+2H_3.
                
                First, we obtain for the term H_1 using (<ref>)
                
                    {|H_1|≥ cE}={|∑^k_0∑^k_0_i≠ t_i^T_t|≥ c k_0(k_0-1) _1-_2^2 }   =   {|∑_j=1^m∑^k_0∑^k_0_i≠ tW_ijW_tj|≥ c k_0(k_0-1) _1-_2^2 }   ≤   m·{|∑^k_0∑^k_0_i≠ tW_ijW_tj|≥ c k_0(k_0-1) _1-_2^2m}
                       ≤   m·{|1k_0∑_i=1^k_0W_ij|≥ c_1 _1-_2√(m)}+m·{|1k_0∑_i=1^k_0W_ij^2|≥ c(k_0-1) _1-_2^2m}   ≤     c_2p^2e^-c_3min{τ,√(nτ)}+p^2{|1k_0∑_i=1^k_0W_ij^2|≥ c(k_0-1) _1-_2^2m},
                
                where the last inequality follows by similar arguments as used in the derivation of   (<ref>) and (<ref>). For the second term, we use the decomposition (<ref>)
                
                    p^2{|1k_0∑_i=1^k_0W_ij^2|≥c2(k_0-1) _1-_2^2m}   ≤    p^2 k_0·{|W_ij|≥√(c2(k_0-1))_1-_2√(m)}   ≤    p^2 k_0·{ |X_iaX_ib-(X_iaX_ib)|≥ c_4√(k_0)_1-_2√(m)}      +2p^2 k_0·{ |X_iaX̅_b|≥ c_4√(k_0)_1-_2√(m)}+p^2 k_0·{ |X̅_aX̅_b|≥ c_4√(k_0)_1-_2√(m)}   ≤    c_5p^2 k_0e^-c_6√(k_0)_1-_2/√(m)≤ c_5p^2 ne^-c_7√(τ)=o(1n),
                
                where the last line uses Lemma <ref> for the sub-exponential random variable X_iaX_ib and the probability of this sub-exponential term is also the leading one among the remaining three terms. Moreover, because |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)√(n/k_0) for all (a,b)∈𝒫 we have
                √(k_0)_1-_2/√(m)≥ C√(τ). Then the order o(1/n)  comes from the assumptions (<ref>) - (<ref>).
                Combining this estimate with (<ref>) gives
                
                    ℙ (|H_1| ≥ cE) ≤ c_2p^2e^-c_3min{τ,√(nτ)}+ c_5p^2 ne^-c_7√(τ)= o(1/n).
                
                Next for the term H_2, similarly  we can calculate
                
                    {|H_2|≥ cE}   =   ℙ{|∑^k∑^k_i≠ t=(k_0+1)_i^T_t|≥ c k_0(k_0-1)n-k[n-kn-k-1(k-k_0)+(k-k_0)+(k-k_0)^2n-k-1]_1-_2^2 }   ≤   m·ℙ{|∑^k∑^k_i≠ t=(k_0+1)W_ijW_tj|≥ c k_0(k_0-1)n-k[2(k-k_0)+(k-k_0)^2n-k-1]_1-_2^2m}   ≤   
                    
                    
                    m·ℙ{|1/k-k_0∑_i=k_0+1^k W_ij|≥√(c2k_0(k_0-1)(n-k)[2(k-k_0)+1n-k-1])_1-_2√(m)}      +m·ℙ{|1k-k_0∑_i=k_0+1^k W_ij^2|≥c2k_0(k_0-1)n-k[2+(k-k_0)n-k-1]_1-_2^2m}.
                
                Using similar arguments as in the discussion of the term (<ref>), the above probability can be further bounded by
                
                    {|H_2|≥ cE}   ≤    c_1m e^-c_2 (k-k_0)min{k_0^2/(n-k_0)(k-k_0)_1-_2^2/m,√(k_0^2/(n-k_0)(k-k_0))_1-_2/√(m)}+c_2mne^-c_3√(k_0^2/n-k_0)_1-_2/√(m)   ≤     c_1p^2e^-c_4min{τ,√(nτ)}+c_2p^2ne^-c_5√(τ)=o(1/n),
                
                where the last line is due to the observation that if the smallest nonzero entry of (Σ_1-Σ_2) satisfies |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)n/k_0√(n-k_0/k_0) for some large C in (<ref>), then
                
                    k_0^2/(n-k_0)_1-_2^2/m≥ Cτ√(k_0^2(k-k_0)/n-k_0)_1-_2/√(m)≥ C√(nτ).
                
                 So together with p^2n=o(e^cτ), p^2n=o(e^c√(nτ))  and p^2n^2=o(e^c√(τ)) in Assumption <ref>, we can find the probability to be of order o(1/n).
                For the term H_3, according to (<ref>) in Assumption <ref>, we have
                
                    {|H_3|≥ cE}≤ m·ℙ{|∑_i=1^k_0∑_t=k_0+1^kW_ijW_tj|≥ c k_0(k_0-1)n-k[2(k-k_0)]_1-_2^2m}   ≤   m·ℙ{|∑_i=1^k_0W_ijk_0|≥√((k_0-1)n-k_0)_1-_2√(m)}
                    +m·ℙ{|∑_t=k_0+1^kW_tjk-k_0|≥ 2c√((k_0-1)n-k_0)_1-_2√(m)}   ≤   c_1m[e^-c_2k_0min{k_0/n-k_0_1-_2^2/m,√(k_0/n-k_0)_1-_2/√(m)}+ e^-c_3(k-k_0)min{k_0/n-k_0_1-_2^2/m,√(k_0/n-k_0)_1-_2/√(m)}]   ≤   c_1me^-c_4min{τ,√(nτ)}=o(1/n)
                
                when the smallest nonzero entry of (Σ_1-Σ_2) satisfies |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)n/k_0√(n-k_0/k_0) for some large C and p^2n=o(e^cτ), p^2n=o(e^c√(nτ)).
                Combining these arguments gives
                
                    ℙ (|A_ℓ (k) - A_ℓ(k_0)| ≥ cE) ≤ c_1p^2[e^-c_2τ+ne^-c_2√(τ) ]= o(1/n),     ℓ=1,2,3,
                
                where we note once again that the cases ℓ =2,3 follow by similar arguments as given for ℓ =1.
                From (<ref>), (<ref>), (<ref>) we therefore obtain
                (<ref>), which proves
                
                    ℙ{k̂≥ (1+ ϵ)k_0 }≤ c_5p^2n[e^-c_6τ+ne^-c_6√(τ)]
                    → 0.
                
                By the discussion at the beginning of the proof and (<ref>)
                the assertion of Theorem <ref> follows.
                
                
                
                
                
                 §.§ Proof of Corollary <ref>
                
                The difference  in proving Theorem <ref>, <ref> and   <ref> under Assumption <ref> and Assumption <ref> consists only in a different treatment of   the terms C^(2) in (<ref>), H_1 in (<ref>) (and H_2 in (<ref>)), for which we need to make use of the following Proposition <ref>. The proof of this result is postponed to Section <ref>.
                
                
                Suppose y_1,⋯, y_k (k≥ n^ϵ for some 0<ϵ<1) are independent sub-exponential random variables.  Let Δ>max_i [y_i^2]. Then
                for any positive constants c>0, M >0 there exists a constant n_0= n_0(c,M) ∈ℕ, such that for all n ≥ n_0.
                
                    {1/k∑_i=1^k y_i^2>Δ}<cn^-M.
                 First, we discuss the differences in the proof of Theorem <ref>   and look at the term C^(2) in (<ref>)
                recalling the representation  Ẋ_iaẊ_ib=X_iaX_ib-X_iaX̅_b-X_ibX̅_a+X̅_aX̅_b. Proposition <ref> gives for the sum corresponding to the first term
                
                    {1/k∑_i=1^k(X_iaX_ib)^2>cτ}≤ c· n^-M,  ∀ M>0.
                
                Moreover, for k=⌊ n^1/2⌋ +1,…, n- ⌊ n^1/2⌋ -1 we have
                
                    {1k∑_i=1^k(X_iaX̅_b)^2>cτ}   ≤     k{|X_iaX̅_b|>c√(τ)}   ≤    n  [{|X_iaX̅_b|>c√(τ),|X̅_b|>1}+{|X_iaX̅_b|>c√(τ),|X̅_b|≤ 1}]   ≤   n  [{|X̅_b|>1}+{|X_ia|>c√(τ)}]≤ c_1n [e^-c_2n+e^-c_3τ],{1k∑_i=1^k(X̅_aX̅_b)^2>cτ}   ≤    k{|X̅_aX̅_b|>c√(τ)}   ≤    n [{|X̅_a|>c_1τ^1/4}+{|X̅_b|>c_2τ^1/4}]≤ c_3p^2n^2e^-c_4n√(τ),
                
                and the probability in  (<ref>)   can be bounded by
                
                    {C^(2)>cτ}≤ c_1[ n^-M+ n^2e^-c_2n+ n^2e^-c_3τ+ n^2e^-c_4n√(τ)].
                
                Consequently, if
                
                    pn^-M→ 0,    p^2n^2e^-cn→ 0,    p^2n^2 e^-cτ→ 0
                
                for some small positive constant c, it follows that p^2·{C^(2)>cτ}→ 0. Here M>0 could be any large positive constant. These estimates show
                 that (<ref>) holds for the case i=2 as long as
                
                
                
                
                
                
                
                
                
                    pn^-M→ 0,    p^2n^2e^-cn→ 0,    p^2n^2 e^-cτ→ 0,   p^2n e^-cn^1/4√(τ)→ 0.
                
                for some small positive constant c.
                Note that  these  conditions contain  (<ref>)  and that
                pn^-M→ 0 implies p^2n^2e^-cn→ 0 and p^2n e^-cn^1/4√(τ)→ 0.
                Consequently,   (<ref>) holds  if
                
                    pn^-M→ 0,    p^2n^2 e^-cτ→ 0,
                
                
                
                
                 where c is some small positive constant and M>0 could be any large constant.
                
                
                
                Next, we discuss the differences in the proof of Theorem <ref> and look exemplarily at the term H_1.
                For the second term in (<ref>), recall that
                 
                    W_ij=X_iaX_ib-(X_iaX_ib)-X_iaX̅_b-X_ibX̅_a+X̅_aX̅_b
                
                 in equation (<ref>).  Proposition <ref> gives
                 
                    {|1k_0∑_i=1^k_0(X_iaX_ib-(X_iaX_ib))^2|≥c2(k_0-1) _1-_2^2m}≤ c n^-M.
                
                 In addition,
                 
                    {|1k_0∑_i=1^k_0(X_iaX̅_b)^2|≥c2(k_0-1) _1-_2^2m}≤   k_0·{ |X_iaX̅_b|≥ c_4√(k_0)_1-_2√(m)}   ≤    n [{|X_iaX̅_b|>c_4√(k_0)_1-_2√(m),|X̅_b|>1}+{|X_iaX̅_b|>c_4√(k_0)_1-_2√(m),|X̅_b|≤ 1}]   ≤    n [{|X̅_b|>1}+{|X_ia|>c_4√(k_0)_1-_2√(m)}]≤ c_1 n[e^-c_2n+e^-c_3k_0_1-_2^2/m],
                
                and
                
                    {|1k_0∑_i=1^k_0(X̅_aX̅_b)^2|≥c2(k_0-1) _1-_2^2m}   ≤    p^2 k_0·{ |X̅_aX̅_b|≥ c_4√(k_0)_1-_2√(m)}   ≤   c_1p^2ne^-c_2n√(k_0)_1-_2√(m).
                
                So, when the smallest nonzero entry of (Σ_1-Σ_2) satisfies |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)√(n/k_0) for some large C,
                √(k_0)_1-_2/√(m)≥ C√(τ), the second term can be bounded by
                
                    {|1k_0∑_i=1^k_0W_ij^2|≥c(k_0-1)2_1-_2^2m}   ≤    c[n^-M+ne^-cn+ne^-ck_0_1-_2^2/m+ne^-cn√(k_0)_1-_2/√(m)]   ≤   c[n^-M+ne^-cn+ne^-cτ+ ne^-cn√(τ)],
                
                and we obtain
                
                    {|H_1|≥ cE}≤  c[p^2e^-cmin{τ,√(nτ)}+p^2n^-M+p^2ne^-cn+p^2ne^-cτ+ p^2ne^-cn√(τ)]=o(1/n),
                
                where  the last estimate follows from  pn^-M→ 0 and  pn e^-cτ→ 0 in Assumption <ref>.
                
                Finally,  the term H_2 in (<ref>)
                can be treated using the same derivation as in (<ref>) and  we obtain
                
                    {|H_2|≥ cE}≤  c_1p^2e^-c_4min{τ,√(nτ)}+c_2p^2[n^-M+ne^-cn+ne^-cτ+ ne^-cn√(τ)]=o(1/n),
                
                where we use the fact  that
                
                    k_0^2/(n-k_0)_1-_2^2/m≥ Cτ
                
                if the smallest nonzero entry of (Σ_1-Σ_2) satisfies |σ_1(a,b)-σ_2(a,b) |≥ C√(τ/n)n/k_0√(n-k_0/k_0) for some large C in (<ref>).
                
                Together with  the conditions pn^-M→ 0 and  pn e^-cτ→ 0  from Assumption <ref> it follows that  the probability is
                of order  o(1/n).
                
                Adjusting the above three terms in the proof of Theorem <ref>, Theorem <ref> and Theorem <ref>, we complete the proof of Corollary <ref>.
                
                
                
                
                 §.§ Proof of Proposition <ref>
                
                Denote C≜Δ-1/k∑_i=1^k  y_i^2, C>0. By   Theorem 4.1 in <cit.> and the following remark, we have for any h>2
                    {1k∑_i=1^k y_i^2>Δ}   =   {1k∑_i=1^k (y_i^2- y_i^2)>C}≤1C^hk^h[∑_i=1^k (y_i^2- y_i^2)]^h
                       ≤   C_1^hC^hk^h(hlog h)^h·max{[|∑_i=1^k (y_i^2- y_i^2)|^2]^h/2, ∑_i=1^k|y_i^2- y_i^2|^h},
                
                where C_1 is an absolute constant. We calculate the upper bounds for the two terms using the fact that the random variables y_i  are sub-exponential. For i=1,…,k we have
                
                    |y_i^2- y_i^2|^h
                       ≤    2^h+1 y_i^2h
                    =2^h+1 [ y_i^2hI{|y_i|≤ h(log h)^2}+ y_i^2hI{|y_i|>h(log h)^2}]   ≤   2^h+1[h^2h(log h)^4h+c_1]≤ c_2· 2^h h^2h(log h)^4h,
                
                and thus for any positive integer M>0, we can find h>2 such that
                
                    C_1^h/C^hk^h(h/log h)^h·∑_i=1^k|y_i^2- y_i^2|^h≤ c_3^h k·[(hlog h)^3/k]^h≤ cn^-M,
                
                where the last inequality is due to the fact that k>n^ϵ.
                Moreover, since the  random variables y_i are independent, we obtain
                
                    [|∑_i=1^k (y_i^2- y_i^2)|^2]^h/2=[∑_i=1^k(y_i^2- y_i^2)^2]^h/2≤ c_4^h k^h/2.
                
                Therefore, for any positive integer M>0, there exists a constant h, such that
                
                    C_1^h/C^hk^h(h/log h)^h·[|∑_i=1^k (y_i^2- y_i^2)|^2]^h/2≤ c_5^h(h/√(k)log h)^h≤ cn^-M,
                
                where the last inequality is also based on the fact that k>n^ϵ. Combining (<ref>), (<ref>) and (<ref>) the assertion of Proposition <ref> follows.
