On variance estimation under shifts in the mean
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Axt, Ieva; Fried, Roland Article — Published Version On variance estimation under shifts in the mean AStA Advances in Statistical Analysis Provided in Cooperation with: Springer Nature Suggested Citation: Axt, Ieva; Fried, Roland (2020) : On variance estimation under shifts in the mean, AStA Advances in Statistical Analysis, ISSN 1863-818X, Springer, Berlin, Heidelberg, Vol. 104, Iss. 3, pp. 417-457, https://doi.org/10.1007/s10182-020-00366-5 This Version is available at: https://hdl.handle.net/10419/288278 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/
Vol.:(0123456789) AStA Advances in Statistical Analysis (2020) 104:417–457 https://doi.org/10.1007/s10182-020-00366-5 1 3 ORIGINAL PAPER On variance estimation undershifts inthemean IevaAxt1· RolandFried1 Received: 18 April 2019 / Accepted: 18 March 2020 / Published online: 1 April 2020 © The Author(s) 2020 Abstract In many situations, it is crucial to estimate the variance properly. Ordinary variance estimators perform poorly in the presence of shifts in the mean. We investigate an approach based on non-overlapping blocks, which yields good results in changepoint scenarios. We show the strong consistency and the asymptotic normality of such blocks-estimators of the variance under independence. Weak consistency is shown for short-range dependent strictly stationary data. We provide recommendations on the appropriate choice of the block size and compare this blocks-approach with difference-based estimators. If level shifts occur frequently and are rather large, the best results can be obtained by adaptive trimming of the blocks. Keywords Blockwise estimation· Change-point· Trimmed mean 1 Introduction We consider a sequence of random variables Y1,…,YN generated by the model Most of the time we assume that X1,…,XN are i.i.d. random variables with E( X t) = 𝜇 and Var( X t) =𝜎 2 , but this will be relaxed occasionally to allow for a short-range dependent strictly stationary sequence. The observed data y1,…,yN are affected by an unknown number K of level shifts of possibly different heights h1,…,hK at different time points t1,…,tK . Our goal is the estimation of the variance 𝜎2 . Without loss of generality, we will set 𝜇=0 in the following. (1) Y t=Xt+ K ∑ k=1 hkIt≥tk . * Ieva Axt [email protected] Roland Fried [email protected] 1 TU Dortmund University, Dortmund, Germany
418 I.Axt, R.Fried 1 3 In Sect.2 we analyse estimators of 𝜎2 from the sequence of observations (Yt)t≥1 by combining estimates obtained from splitting the data into several blocks. Without the need of explicit distributional assumptions the mean of the blockwise estimates turns out to be consistent if the size and the number of blocks increases, and the number of jumps increases slower than the number of blocks. If many jumps in the mean are expected to occur, an adaptively trimmed mean of the blockwise estimates can be used, see Sect.3. In Sect.4 a simulation study is conducted to assess the performance of the proposed approaches. In Sect.5 the estimation procedures are applied to real data sets, while Sect.6 summarizes the results of this paper. 2 Estimation ofthevariance byaveraging When dealing with independent identically distributed data the sample variance is the common choice for estimation of 𝜎2 . However, if we are aware of a possible presence of level shifts at unknown locations, it is reasonable to divide the sample Y1,…,YN into m non-overlapping blocks of size n=⌊N∕m⌋ and to calculate the average of the m sample variances derived from the different blocks. A similar approach has been used in Dai etal. (2015) in the context of repeated measurements data and in Rooch etal. (2019) for estimation of the Hurst parameter. The blocks-estimator 𝜎 2 Mean of the variance investigated here is defined as where S 2 j = 1 n−1∑n t=1 (Yj,t−Yj) 2 , Y j= 1 n∑n t=1 Yj, t and Yj,1,…,Yj,n are the observations in the jth block. We are interested in finding the block size n which yields a low mean squared error (MSE) under certain assumptions. In what follows, we will concentrate on the situation where all jump heights are positive. This is a worse scenario than having both, positive and negative jumps, since the data are more spread in the former case resulting in a larger positive bias of most scale estimators. 2.1 Asymptotic properties We will use some algebraic rules for derivation of the expectation and the variance of quadratic forms in order to calculate the MSE of 𝜎 2 Mean , see Seber and Lee (2012). Let B be the number of blocks with jumps in the mean and K≥B the total number of jumps. The expected value and the variance of 𝜎 2 Mean are given as follows: (2) 𝜎 2Mean = 1 m m ∑ j=1 S2 j ,
419 1 3 On variance estimation undershifts inthemean where 𝜈4 =E ( X4 1) ,A=𝕀 n − 1 n 1 n 1 T n , 𝕀n is the unit matrix, 1n=(1, …,1)T and 𝜇j contains the expected values of the random variables in the perturbed block j=1, …,B , i.e., 𝜇j =(𝜇 j,1 ,…,𝜇 j,n ) T =(E(Y j,1 ),…,E(Y j,n )) T . The term 𝜇T j A𝜇j∕(n−1 ) is the empirical variance of the expected values E(Yj,1),…,E(Yj,n) in block j. In a jumpfree block, we have 𝜇T j A𝜇j= 0 , since all expected values and therefore the elements of 𝜇j are equal. The blocks-estimator (2) estimates the variance consistently if the number of blocks grows sufficiently fast as is shown in Theorem1. Theorem1 Let Y1,…,YN with Yt =X t + ∑K k=1 h k I t ≥ tk from Model (1) be segregated into m blocks of size n, where t1,…,tK are the time points of the jumps of size h1,…,hK , respectively. Let B out of m blocks be contaminated by K1,…, KB jumps, respectively, with ∑B j=1 Kj=K=K(N ) . Moreover, let E(|X1|4)<∞ , K�∑ K k=1hk �2 =o(m ) and m→∞ , whereas the block size n can be fixed or increasing as N→∞ . Then 𝜎 2 Mean = 1 m∑ m j =1 S2 j →𝜎 2 almost surely. Proof Without loss of generality assume that the first B out of m blocks are contaminated by K1,…, KB jumps, respectively. Let the term S2 j,0 denote the empirical variance of the uncontaminated data in block j, while S2 j,h is the empirical variance when Kj level shifts are present. Moreover, Yj,1,…,Yj,n are the observations in the jth block, 𝜇j,t =E ( Y j,t) and 𝜇 j= 1 n∑n t=1 E � Yj,t � . Then we have For the second term in the last Eq.(3), we have almost surely E ( 𝜎2Mean ) =𝜎2+1 m(n−1) B ∑ j=1 𝜇T jA𝜇j, Var( 𝜎2Mean ) =1 m ( 𝜈4 n−𝜎4(n−3) n(n−1) ) +4𝜎2 m2(n−1)2 B ∑ j=1 𝜇T jA𝜇j , (3) 𝜎 2Mean =1 m m ∑ j=1 S2 j=1 m m ∑ j=B+1 S2 j,0 +1 m B ∑ j=1 S2 j,h =1 m m ∑ j=B+1 S2 j,0 +1 m B ∑ j=1 1 n−1 n ∑ t=1(Xj,t+𝜇j,t−Xj−𝜇j) 2 =1 m m ∑ j=1 S2 j,0 +1 m B ∑ j=1 2 n−1 n ∑ t=1 (Xj,t−Xj)(𝜇j,t−𝜇j) +1 m B ∑ j=1 1 n−1 n ∑ t=1( 𝜇j,t−𝜇j ) 2.
420 I.Axt, R.Fried 1 3 The term 1 B∑ B j=1 1 n∑ n t=1 � � � (Xj,t−Xj) � � � in (4) is a random variable with finite moments if n and B are fixed. This random variable converges to the term E� 1 n∑ n t=1 ��� (Xj,t−Xj) ���� almost surely if B→∞ . In the case of B→∞ and n→∞ this term converges to E(| X 1|) almost surely due to Theorem2 of Hu etal. (1989) and the condition E(|X1|4)<∞ , since S 2 j −E ( S2 j ) are uniformly bounded with P ( | S2 j−E ( S2 j )| >t)→0∀ t due to Chebyshev’s inequality and Var (S 2 j )→ 0 . Moreover, we used the fact that B� � �∑ K k=1hk � � � ≤K � � �∑ K k=1hk � � � =o(m ) . The following is valid for the third term in (3): The first term of the last equation in (3) converges almost surely to 𝜎2 due to the results on triangular arrays in Theorem2 of Hu etal. (1989), assuming that the condition E(|X1|4)<∞ holds, since S 2 j −E ( S2 j ) are uniformly bounded with P ( | S2 j−E ( S2 j )| >t)→0∀ t due to Chebyshev’s inequality and Var (S 2 j )→ 0 . Application of Slutsky’s Theorem proves the result. ◻ Remark 2 1. If the jump heights are bounded by a constant h≥hk,k=1, …,K, the strongest restriction arises if all heights equal this upper bound resulting in the constraint K�∑ K k=1hk �2 =K3h2=o(m ) . Consistency is thus guaranteed if the number of blocks grows faster than K3 . (4) | | | | | | 2 m(n−1) B ∑ j=1 n ∑ t=1 (Xj,t−Xj)(𝜇j,t−𝜇j) | | | | | | ≤2 m(n−1) B ∑ j=1 n ∑ t=1| | |(Xj,t−Xj)| | || | |(𝜇j,t−𝜇j)| | | ≤2 m(n−1) B ∑ j=1 n ∑ t=1| | |(Xj,t−Xj)| | || | | | | | K ∑ k=1 hk| | | | | | =B | | | | | | K ∑ k=1 hk | | | | | | 2 m n n−1 1 B B ∑ j=1 1 n n ∑ t=1 | | | (Xj,t−Xj) | | | ⟶ 0. 1 m B ∑ j=1 1 n−1 n ∑ t=1(𝜇j,t−𝜇j)2≤1 m B ∑ j=1 1 n−1 n ∑ t=1 ( K ∑ k=1 hk )2 =B m n n−1 ( K ∑ k=1 hk ) 2 ⟶0.
421 1 3 On variance estimation undershifts inthemean 2. By the Central Limit Theorem, the estimator 𝜎 2 Mean is asymptotically normal if no level shifts are present and the block size n is fixed. Its asymptotic efficiency relative to the ordinary sample variance is in case of i.i.d data with finite fourth moments, where (see Angelova (2012) for the variance Var (S 2 1) of the sample variance). E.g., under normality the efficiency of the blocks estimater with fixed n is (n−1)n−1<1. 3. The asymptotic efficiency of the blocks-estimator relative to the sample variance is 1 if n→∞ . The next Theorem shows that 𝜎 2 Mean is asymptotically not only normal but even fully efficient in case of a growing block size. Theorem3 Assume that the i.i.d. random variables Y1,…,YN are segregated into m blocks of size n, with m,n→∞ such that m=o(n),n=o(N) . Moreover, assume that 𝜈4 =E ( X 4 1) < ∞ . Then we have Proof Rewriting the estimator 𝜎 2 Mean we get For the second term of the numerator in (5), we have that Var(S2) Var( 𝜎2 Mean ) = 𝜈 4 N−𝜎 4 (N−3) N(N−1) 1 m( 𝜈4 n −𝜎4(n−3) n(n−1)) = 𝜈4−𝜎 4 (N−3) (N−1) 𝜈4−𝜎4(n−3) (n−1) N→∞ ⟶ 𝜈4−𝜎4 𝜈4−𝜎4(n−3) (n−1) Var( 𝜎2 Mean ) =Var ( 1 m m ∑ j =1 S2 j ) =1 m Var ( S2 1 ) =1 m ( 𝜈4 n−𝜎4(n−3) n(n−1) ) √ N � 𝜎2Mean −𝜎2 �d ⟶N � 0, 𝜈4−𝜎4 �. (5) 𝜎2Mean −𝜎2 � Var� 𝜎2Mean� = 1 m(n−1) ∑ m j=1 ∑ n t=1 � Xj,t−Xj �2 −𝜎2 �1 m�𝜈4 n−𝜎4(n−3) n(n−1)� = √ N 1 m(n−1)∑m j=1∑n t=1X2 j,t−n m(n−1)∑m j=1X2 j−𝜎2 � 𝜈4−𝜎4(n−3) n−1 .
422 I.Axt, R.Fried 1 3 since m=o(n) . Convergence of the term (6) in mean implies convergence in probability to zero. Application of the Central Limit Theorem to the remaining terms of (5) yields the desired result. ◻ Remark 4 In the proof of Theorem3, we have assumed that m=o(n) , i.e., the block size grows faster than the number of blocks. This condition can be dropped using the Lyapunov condition under the assumption of finite eighth moments, as will be shown in the following. We set with E(Ti,j)=0 and ∑m i j=1 E(T 2 i,j )=1∀i , where i denotes the ith row of the triangular array. The Lyapunov condition [see corollary 1.9.3 in Serfling (1980)] is the following: With 𝛿=2 and existing moments 𝜈1,…,𝜈8 of X1 we get where g is a function of the existing moments 𝜈1,…,𝜈8 with g(𝜈1,…,𝜈8)=O(1) . See Angelova (2012) for the fourth central moment of the sample variance. Therefore, the condition m=o(n) can be dropped. (6) E�� ���� � √Nn m(n−1) m � j=1 X2 j � ���� �� = √ N m n n−1 m � j=1 E�X2 j�=√Nn n−1 𝜎 2 n = √ mn n n−1 𝜎2 n = � m n n n−1 𝜎2→0, T i,j= S 2 i,j− 𝜎2 √ mi ni ( 𝜈4−𝜎4(ni−3) ni−1 ) , ∃ 𝛿>0∶lim i→∞ m i ∑ j=1 E (| Ti,j | 2+𝛿 ) = 0 m i ∑ j=1 E(|Ti,j|4)= m i ∑ j=1 1 m2 i n2 i(𝜈4−𝜎4(ni−3) ni−1)2 ⋅E (( S2 i,j−𝜎2 ) 4 ) =mi m2 i n2 i(𝜈4−𝜎4(ni−3) ni−1)2 ⋅O(1 n2 i)⋅g(𝜈1,…,𝜈8 ) =o ( 1 m i) →0,
423 1 3 On variance estimation undershifts inthemean 2.2 Choice oftheblock size When choosing blocks of length n=2 , the estimator 𝜎 2 Mean results in a differencebased estimator which considers ⌊N∕2⌋ consecutive non-overlapping differences: Difference-based estimators have been considered in many papers, see Von Neumann etal. (1941), Rice (1984), Gasser etal. (1986), Hall etal. (1990), Dette etal. (1998), Munk etal. (2005), Tong etal. (2013), among many others. Dai and Tong (2014) discussed estimation approaches based on differences in nonparametric regression context, Wang et al. (2017) considered an estimation technique which involves differences of second order, while Tecuapetla-Gómez and Munk (2017) proposed a difference-based estimator for m-dependent data. An ordinary differencebased estimator of first order, which considers all N−1 consecutive differences, is [see, e.g., Von Neumann etal. (1941)]: where A= AT A and A such that AY = ( Y 2 −Y 1 ,…,Y N −Y N−1)T. Theorems 1.5 and 1.6 in Seber and Lee (2012) are used to calculate the expectation and the variance of the estimator 𝜎 2 Diff : where 𝜇=E(Y) , 𝜈i =E ( X i 1) and a is a vector of the diagonal elements of A. When no changes in the mean are present both, the difference-based and the averaging estimators, are unbiased. For the variance of the estimators we have For example, for N=100 and a block size n=10, we get Var( 𝜎2 Mean ) = 0.0222 , while Var( 𝜎2 Diff ) = 0.0302 . Therefore, when no changes in the mean are present the block-estimator can have smaller variance than the difference-based estimator. In Sect. 4 we compare the performance of the proposed estimation procedures with that of the difference-based estimator (7) in different change-point scenarios. 𝜎 2Mean,n=2=1 2 ⌊ N∕2 ⌋ ⌊N∕2⌋ � j=1 � Y2j−Y2j−1 � 2 . (7) 𝜎 2Diff =1 2(N−1) N−1 ∑ j=1 ( Yj+1−Yj ) 2=1 2(N−1)YTAY , E ( 𝜎2Diff ) =𝜎2+ 1 2(N−1)𝜇TA𝜇, Var( 𝜎2Diff ) =1 4(N−1) 2 ( 𝜈4(4N−6)+2𝜎4+4𝜎2𝜇TA2𝜇+4𝜈3𝜇TAa ), Var ( 𝜎2Diff ) =𝜈4 4N−6 4(N−1)2+𝜎 4 2(N−1)2 , Var( 𝜎2Mean ) = 𝜈4 N −𝜎4(n−3) N(n−1) .
424 I.Axt, R.Fried 1 3 In the following, we investigate the proper choice of the block size for the estimator 𝜎 2 Mean . All calculations in this paper have been performed with the statistical software +R+, version 3.5.2, R Core Team (2018). For known jump positions the MSE of the blocks-variance estimator 𝜎 2 Mean can be determined analytically. The position of the jump is relevant for the performance of this approach. Therefore, it is reasonable to consider different positions of the K jumps to get an overall assessment of the performance of the blocks-estimator. For every K∈{1, 3, 5} , we generate K jumps of equal heights h=𝛿⋅𝜎 , with 𝛿∈{0, 0.1, 0.2, …, 4.9, 5} , at positions sampled randomly from a uniform distribution on the values maxn(N−⌊N∕n⌋n)+1, …,N−maxn(N−⌊N∕n⌋n) without replacement, and calculate the MSE for every reasonable block size n∈{2, 3, 4, …,⌊N∕2⌋} . This is repeated 1000 times, leading to 1000 MSE values for every h and n based on different jump positions. The average of these MSE values is taken for each h and n. Data are generated from the standard normal or the t5 -distribution. Panel (a) of Fig.1 shows the block size nopt which yields the least theoretical MSE value of the estimator 𝜎 2 Mean depending on the jump height h=𝛿⋅𝜎 with K∈{1, 3, 5} jumps for N=1000 observations and normal distribution. We observe that nopt decreases for 𝜎 2 Mean as the jump height grows. Blocks of size 2 (resulting in a non-overlapping difference-based estimator) are preferred when h≈4𝜎 and K=5 , while larger blocks lead to better results in case of smaller or less jumps. 012345 02040608 0 δ MSE−optimal block length (a) 012345 0.0020.004 0.006 δ MSE (b) 012345 0.0020.004 0.006 δ MSE (c) Fig. 1 a MSE-optimal block length nopt of 𝜎 2 Mean , b MSE regarding nopt of 𝜎 2 Mean and c MSE of 𝜎 2 Mean when choosing n = √ N K+1 for K=1 (solid line), K=3 (dashed line) and K=5 (dotted line) with N=1000 , Y t=Xt+ ∑K k=1 hIt≥t k , where Xt∼N(0, 1) and h=𝛿 ⋅ 𝜎,𝛿∈{0, 0.1, …,5} .
431 1 3 On variance estimation undershifts inthemean due to absolute summability of the autocovariance function 𝛾 and the condition K�∑ K k=1hk �2 =o(mlog(N)−1 ) . Therefore, the term A2,N converges to zero in rth mean with r=2 , which implies convergence in probability. The argument in the probability A3,N is deterministic. Set 𝜖N=log(N)−1. For large N we have E( 1 m B ∑ j=1 2 n−1 n ∑ t=1 Xj,t(𝜇j,t−𝜇j) ) =1 m B ∑ j=1 2 n−1 n ∑ t=1 E(Xj,t)(𝜇j,t−𝜇j) =1 m B ∑ j=1 2𝜇 n−1 n ∑ t=1 (𝜇j,t−𝜇j)=0 Var(1 m B ∑ j=1 2 n−1 n ∑ t=1 Xj,t(𝜇j,t−𝜇j)) =4 m2(n−1)2Var(B ∑ j=1 n ∑ t=1 Xj,t(𝜇j,t−𝜇j)) =4 m2(n−1)2Cov(B ∑ j=1 n ∑ t=1 Xj,t(𝜇j,t−𝜇j), B ∑ l=1 n ∑ s=1 Xl,s(𝜇l,s−𝜇l) ) ≤4 m2(n−1)2 B ∑ j=1 n ∑ t=1 B ∑ l=1 n ∑ s=1|𝜇j,t−𝜇j||𝜇l,s−𝜇l||Cov(Xj,t,Xl,s) | ≤4 m2(n−1)2(K ∑ k=1 hk)2m ∑ j=1 n ∑ t=1 m ∑ l=1 n ∑ s=1|Cov(Xj,t,Xl,s)| =4 m2(n−1)2(K ∑ k=1 hk)2(NVar(X1,1)+2 N−1 ∑ u=1|𝛾(u)|(N−u)) ≤4N m2(n−1)2(K ∑ k=1 hk)2(Var(X1,1)+2 ∞ ∑ u=1 | 𝛾(u) | ) ⟶0,
432 I.Axt, R.Fried 1 3 since K�∑ K k=1hk �2 =o(mlog(N)−1 ) holds. Altogether we get the result 𝜎 2 Mean p →𝜎 2 . ◻ Remark 7 If the jump heights are bounded by a constant h ≥ hk,k=1, …,K, the strongest restriction arises if all heights equal this upper bound resulting in the constraint K3 h 2 =o ( mlog(N) −1) . 3 Trimmed estimator ofthevariance So far we have considered cases where the number of changes in the mean K is rather small with respect to the number of the blocks m and thus to the number of observations N. However, there might be situations in which level shifts occur frequently. Asymptotically, the blocks-estimator 𝜎 2 Mean [see (2)] is still a good choice for the estimation of the variance as long as the number of level shifts grows slowly, see e.g. Theorem1. However, if many jumps are present in a finite sample the blocks-estimator is no longer a good choice and will become strongly biased. We propose an asymmetric trimmed mean of the blockwise estimates instead of their ordinary average, i.e., large estimates are removed and the average value of the remaining ones is calculated. We do not consider a symmetric trimmed mean, since the sample variance is positively biased in the presence of level shifts, so that estimates from blocks containing a level shift are expected to show up as upper outliers. Moreover, we suggest using rather many small blocks to account for potentially many level shifts. In the next Sects.3.1 and 3.2, the choice of the trimming fraction is discussed. 3.1 Estimation withafixed trimming fraction The trimmed blocks-estimator is given as A 3,N≤P ( 1 m B ∑ j=1 1 n−1 n ∑ t=1(𝜇j,t−𝜇j)2> 𝜖N 3 ) =P(log(N)1 m B ∑ j=1 1 n−1 n ∑ t=1(𝜇j,t−𝜇j)2>1 3 ) ≤P(log(N)B m n n−1(K ∑ k=1 hk)2 >1 3) ≤P ( log(N)K m n n−1 ( K ∑ k=1 hk ) 2 >1 3 ) =0,
433 1 3 On variance estimation undershifts inthemean where S2 (1) ≤⋯≤S 2 (m) are the ordered blockwise estimates, m is the number of blocks and CN,Tr,𝛼 is a sample and distribution dependent correction factor to ensure unbiasedness in the absence of level shifts. In practice this constant can be simulated under the assumption of observing data from a known location-scale family. For example, for standard normal distribution, 𝛼=0.2 , N=1000 and n=20 ( m=50 ) we generate 1000 samples of length N=1000 and calculate the average of the uncorrected trimmed variance estimates. The reciprocal of this average value yields C1000,Tr,0.2 =1.198 . As an example, we generate 1000 time series of length N∈{1000, 2500} from normal and t5 -distribution. We add K=N ⋅ p jumps to the generated data at randomly chosen positions, as was done in Sect. 2.2, with p∈{0, 2∕1000, 4∕1000, 6∕1000, 10∕1000} and height h∈{0, 2, 3, 5} . We choose n=20 to ensure that the number of jump-contaminated blocks is sufficiently smaller than the total number of blocks. Table 2 shows the simulated MSE of the trimmed estimator (12) for 𝛼∈{0.1, 0.3, 0.5} . Clearly, the performance of the trimmed estimator depends on the number of jumps in the mean and the trimming parameter 𝛼 . Larger values of 𝛼 are required when dealing with many jumps but lead to an increased MSE if there are only a few jumps. Therefore, it is reasonable to choose 𝛼 adaptively, as will be described in the next Sect.3.2. 3.2 Adaptive choice ofthetrimming fraction Instead of using a fixed trimming fraction, we can choose 𝛼 adaptively, yielding the adaptive trimmed estimator 𝜎 2 Tr,ad with (12) 𝜎 2Tr,𝛼=CN,Tr,𝛼 1 m− ⌊ 𝛼m ⌋ m−⌊𝛼m⌋ � j=1 S2 (j) , Table 2 Simulated MSE ⋅102 of 𝜎 2 Tr,𝛼 for normally and t5 -distributed data with N=1000 and different 𝛼 , h and K K h N(0,1) t5 𝛼=0.1 𝛼=0.3 𝛼=0.5 𝛼=0.1 𝛼=0.3 𝛼=0.5 0 0 0.24 0.28 0.35 2.39 4.56 7.05 2 2 0.28 0.30 0.38 1.89 3.96 6.39 8 0.26 0.28 0.35 1.50 3.51 6.01 4 2 0.35 0.34 0.39 1.49 3.12 5.50 8 0.53 0.43 0.46 1.61 2.76 5.13 6 2 0.54 0.46 0.48 1.22 2.47 4.81 8 2.80 0.54 0.51 8.74 1.99 4.32 10 2 1.14 0.76 0.67 2.03 1.58 3.56 8 69.38 1.04 0.77 199.96 1.43 2.79
434 I.Axt, R.Fried 1 3 where S2 (1) ≤⋯≤S 2 (m) are the ordered blockwise sample variances and 𝛼adapt is the adaptively chosen percentage of the blocks-estimates which will be removed. 3.2.1 Estimation undernormality We use the approach for outlier detection discussed in Davies and Gather (1993) to determine 𝛼adapt , assuming that the underlying distribution is normal. In this case, the distribution of the sample variance is well known, i.e., in block j we have that Since the true variance 𝜎2 is not known, we propose to replace 𝜎2 by an appropriate initial estimate, such as the median of the blocks-estimates, i.e., where CN,Med is a finite sample correction factor. Subsequently, we remove those values which exceed q 𝜒2 n−1 ,1−𝛽 m , the (1−𝛽m) -quantile of the 𝜒2 n−1 -distribution, with We will refer to the adaptively trimmed estimator based on the approach of Davies and Gather (1993) under normality as 𝜎 2 norm Tr,ad . The choice of 𝛽m results in a probability of 1−𝛽 that no observation (block in our case) is trimmed if T1,…,Tm are i.i.d. 𝜒2 n−1 -distributed, i.e., Furthermore, we expect that roughly m ⋅ 𝛽m blocks are trimmed on average in the absence of level shifts. The following simulations suggest that the adaptive trimming fraction is slightly larger than 𝛽m which can be explained by the fact that we need to use an estimate such as 𝜎 2 Med instead of the unknown 𝜎2 . We generated 10,000 sequences of observations of size N∈{1000, 2500} for 𝛽∈{0.05, 0.1} . Table3 shows the average number of trimmed blocks in the absence of level shifts. (13) 𝜎 2Tr,ad =CN,Tr,ad 1 m− ⌊ 𝛼adaptm ⌋ m−⌊𝛼 adapt m⌋ � j=1 S2 (j) , T j∶= (n−1)S 2 j 𝜎 2∼𝜒2 n−1 . (14) 𝜎 2 Med =C N,Med ⋅med{S2 1 ,…,S2 m } , (15) T j= (n−1)S 2 j 𝜎 2 Med , (16) 𝛽m =1−(1−𝛽) 1∕m and 𝛽∈(0, 1) . P (T1≤q𝜒2 n−1,1−𝛽m,…,Tm≤q𝜒2 n−1,1−𝛽m)= ( P(T1≤q𝜒2 n−1,1−𝛽m) )m = ( 1−𝛽 m) m= ( (1−𝛽)1∕m ) m=1−𝛽 .
435 1 3 On variance estimation undershifts inthemean We suggest choosing a small block size, e.g. n=20 , to cope with a possibly large number of change-points. In this way, it is ensured that the number of uncontaminated blocks is much larger than the number of perturbed blocks. Moreover, we choose 𝛽=0.05 . The correction factor in (13) needs to be simulated taking into account that the percentage of the omitted block-estimates is no longer fixed. Therefore, for given N and 𝛽 we generate 1000 sequences of observations. In each simulation run we calculate the block-estimates S2 j ,j=1, …,m , and the initial estimate of the variance 𝜎 2 Med . Subsequently, we remove the values T j=(n−1)S2 j ∕ 𝜎2 Med which exceed the quantile q 𝜒2 n−1 ,1−𝛽 m . Then the average value of the remaining block-estimates is computed. The procedure yields 1000 estimates. The correction factor is the reciprocal of the average of these values. For N=1000 and 𝛽=0.05 the simulated correction factor is C0.05 1000,Tr,ad = 1.0020 , while for N=2500 we have C0.05 2500,Tr,ad = 1.0009 , so both are nearly 1 and could be neglected with little loss. 3.2.2 Estimation underunknown distribution When no distributional assumptions are made and the block size n is large, one can use that √ n � S2 j−𝜎2 � ∕ √ 𝜈4−𝜎 4 is approximately standard normal if the fourth moment exists. The fourth central moment 𝜈 4=E (( X1−E(X1) ) 4 ) needs to be estimated properly in the presence of level shifts. We can estimate it in blocks and then compute the median of the blocks-estimates 𝜇 4,Med , as was done in (14). Then, values which exceed the (1−𝛽m) -quantile of the standard normal distribution are removed. The corresponding adaptively trimmed estimator is denoted as 𝜎 2 other Tr,ad , where, again, 𝛽=0.05 [see (16)] will be used in what follows. Table 4 shows the average number of trimmed blocks in the absence of level shifts for normally distributed data, analogously to Table 3. We observe that the average number of trimmed block-estimates is much larger than the values for the trimming procedure which is based on the normality assumption. Remark 8 We do not use distribution dependent correction factors for the estimators 𝜇 4,Med and 𝜎 2 Med , since the underlying distribution is not known. The asymptotic (17) T j= √ n � S2 j− 𝜎2Med � ∕ � 𝜇 4,Med − 𝜎22 Med , Table 3 Average number of trimmed blocks in the absence of level shifts for normally distributed data and different N and 𝛽 for the estimator 𝜎 2 norm Tr,ad 𝛽 N 1000 2500 0.05 0.0647 0.0583 0.1 0.1268 0.1208
436 I.Axt, R.Fried 1 3 distribution of the blockwise empirical variances is normal and therefore, for large block sizes, the distribution is expected to be roughly symmetric. A correction factor can then be neglected with little loss, since the sample median of symmetrically distributed random variables estimates their expected value. 3.3 Choice ofthetrimming fraction based ontheMOSUM procedure In Sect.2.4 we have used the MOSUM procedure proposed by Eichinger and Kirch (2018) to estimate the unknown number of jumps K. The corresponding estimate K can also be used to determine the trimming fraction 𝛼 of the trimmed estimator in (12), i.e., we can set 𝛼= K∕m . We will denote the corresponding estimator as 𝜎 2 mos Tr,ad . We compare the average number of trimmed blocks for the three estimators 𝜎 2 mos Tr,ad , 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad in a simulation study. We generate 1000 sequences of observations of length N=1000 for different values of K and h. To every sequence we add K jumps of height h at randomly chosen positions, as was done in Sect.2.2. Table 5 shows the average number of trimmed blocks for data sampled from the standard normal distribution and the corresponding MSE in different scenarios. Table 4 Average number of trimmed blocks in the absence of level shifts for normally distributed data and different values of N and 𝛽 for the estimator 𝜎 2 other Tr,ad 𝛽 N 1000 2500 n = 20 0.05 1.2327 2.1748 0.1 1.6007 2.8159 n=40 0.05 0.4018 0.6036 0.1 0.5676 0.8829 Table 5 Average number of trimmed blocks ( K ) and resulting MSE for normally distributed data, N=1000 and different values of K and h K h K MSE 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 Mean 0 0 0.09 0.11 0.41 0.23 0.23 0.26 0.22 2 2 2.06 0.45 1.38 0.25 0.26 0.29 0.25 2 5 2.06 1.82 2.13 0.24 0.23 0.26 1.20 2 8 2.07 1.96 2.20 0.23 0.23 0.25 6.71 6 2 5.52 1.10 3.05 0.30 0.53 0.47 0.33 6 5 5.63 5.16 5.26 0.29 0.26 0.29 1.91 6 8 5.68 5.51 5.43 0.43 0.23 0.27 11.20 10 2 8.40 1.66 4.20 0.34 1.11 1.11 0.46 10 5 8.72 8.16 7.89 0.37 0.33 0.39 2.18 10 8 8.78 8.72 8.22 0.84 0.27 0.32 11.65
437 1 3 On variance estimation undershifts inthemean Moreover, the MSE of the averaging estimator 𝜎 2 Mean is displayed in the table for comparison. We observe that 𝜎 2 mos Tr,ad trims more blocks on average and yields a smaller MSE than 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad if the jumps are small. The MSE of the averaging estimator 𝜎 2 Mean is similar to that of 𝜎 2 mos Tr,ad in that case. As opposed to this, the three trimmed estimators trim similarly many blocks on average if the jumps are large, with 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad leading to smaller MSE values than 𝜎 2 mos Tr,ad . This can be explained by the fact that the former methods base the decision whether to trim a block or not directly on the criterion whether the empirical variance in a block is unusual or not. The averaging estimator 𝜎 2 Mean is outperformed by the three trimmed estimators when dealing with large jumps. For the heavy-tailed t5 -distribution, the results are found in Table6. The MOSUMbased procedure yields slightly better results for large jump heights, i.e., h>2 , since the adaptively trimmed estimators overestimate the number of contaminated blocks in many cases. For small jumps, the adaptively trimmed procedures yield smaller MSE values. The three trimmed estimators perform better than 𝜎 2 Mean in every case. The results for the data generated by the AR(1) model with 𝜙=0.5 are displayed in Table7. The MOSUM-based procedure trims considerably more blocks than the adaptively trimmed estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad yielding higher MSE values. Moreover, the averaging estimator 𝜎 2 Mean is outperformed by the adaptively trimmed estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad in every case. We conclude that the estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad perform better when dealing with large jumps and normally distributed data. For small jumps less contaminated blocks are trimmed by the adaptively trimmed estimators resulting in a higher MSE values in many scenarios. In this case the MOSUM-based estimator 𝜎 2 mos Tr,ad yields the best results which are similar to those of 𝜎 2 Mean . For heavy-tailed data, the MOSUM-based procedure is preferable under large jumps but is inferior to the other Table 6 Average number of trimmed blocks ( K ) and resulting MSE for t5 -distributed data, N=1000 and different values of K and h K h K MSE 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 Mean 0 0 0.07 2.04 1.56 2.17 2.98 2.81 2.11 2 2 2.06 2.26 1.85 2.16 2.54 2.33 2.29 2 5 2.05 3.55 2.99 1.63 2.72 2.46 4.01 2 8 2.04 3.73 3.06 1.74 2.90 2.61 8.80 6 2 5.45 2.38 2.26 3.81 1.91 1.91 1.85 6 5 5.57 6.12 5.43 1.98 2.58 2.51 3.93 6 8 5.65 6.86 5.87 2.19 2.78 2.56 15.48 10 2 8.32 2.60 2.79 4.97 1.61 2.94 2.87 10 5 8.64 8.48 7.41 2.04 2.43 3.74 4.57 10 8 8.74 9.63 8.18 2.49 2.53 3.06 13.99
438 I.Axt, R.Fried 1 3 trimmed estimates if the jumps are small. When dealing with positively correlated processes the adaptively trimmed estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad perform best. Therefore, the choice of the most appropriate estimator depends on the knowledge about the underlying distribution and the height of the jumps. In the simulation study in Sect.4, we concentrate on the adaptively trimmed estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad since they yield very good results in many scenarios, and on the averaging estimator 𝜎 2 Mean since it performs well when dealing with a few small jumps and independent data. 4 Simulations In this section we compare the estimators 𝜎 2 Mean , 𝜎 2 Diff , 𝜎 2 Tr,0.5 with the block size n=20 , 𝜎 2 norm Tr,ad with the block size n=20 and 𝛽=0.05 , 𝜎 2 other Tr,ad with the block size n=40 and 𝛽=0.05 , 𝜎 2 mosum Mean and 𝜎 2 mosum W in different scenarios. We generate 1000 sequences of observations of length N∈{200, 1000, 2500} from the standard normal and the t5 distribution. We add K jumps of heights h∈{0, 2, 3, 5, 8} to the data at randomly chosen positions as was done in Sect.2.2. K is chosen dependent on the number of observations, i.e., K=p ⋅ N with p∈{0, 2∕1000, 4∕1000, 6∕1000, 10∕1000} . Table 8 shows the simulated MSE for the normal distribution. The estimators 𝜎 2 Mean and 𝜎 2 mosum Mean yield similar results. We conclude that the estimation of the number of jumps K [required in the rule (8)] does not have a large effect on the estimator. The estimators 𝜎 2 Mean and 𝜎 2 mosum Mean yield, the best results if the jump heights are not very large, i.e., h≤2 . However, the MSE of taking the ordinary average is much larger than that of the other estimators if the jump heights are large. Large jumps result in large blockwise estimates, which have a strong impact on the ordinary average. The trimmed estimator 𝜎 2 norm Tr,ad yields the best results among all methods considered here for normally distributed data if the jumps are rather high. When many Table 7 Average number of trimmed blocks ( K ) and resulting MSE for data generated from the AR(1) model with 𝜙=0.5 and N=1000 for different values of K and h K h K MSE 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 mos Tr,ad 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 Mean 0 0 5.08 0.66 1.27 6.01 2.69 1.84 1.17 2 2 6.79 1.04 1.84 6.77 2.51 1.64 5.46 2 5 6.47 2.29 2.75 5.94 2.69 1.80 2.65 2 8 6.37 2.48 2.83 5.72 2.72 1.72 1.32 6 2 9.41 1.54 2.82 7.26 1.74 1.16 20.85 6 5 9.00 5.33 5.55 5.43 2.50 1.64 12.41 6 8 8.96 5.80 5.71 5.26 2.58 1.54 3.29 10 2 11.53 1.99 3.65 7.37 1.22 1.04 42.02 10 5 11.26 8.02 7.87 4.72 2.06 1.37 29.70 10 8 11.30 8.92 8.27 4.52 2.40 1.35 13.75
439 1 3 On variance estimation undershifts inthemean Table 8 Simulated MSE ⋅102 of 𝜎 2 Mean , 𝜎 2 Diff , 𝜎 2 Tr,0.5 , 𝜎 2 norm Tr,ad , 𝜎 2 other Tr,ad , 𝜎 2 mosum Mean and 𝜎 2 mosum W for normally distributed data and different sample sizes N, jump heights h⋅𝜎 and number of jumps K=p ⋅ N with p∈{0, 2∕1000, 4∕1000, 6∕1000, 10∕1000} K h 𝜎 2 Mean 𝜎 2 Diff 𝜎 2 Tr,0.5 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 mosum Mean 𝜎 2 mosum W N=200 0 0 1.10 1.51 1.64 1.10 1.38 1.19 0.98 1 2 1.34 1.54 1.77 1.54 1.79 1.29 2.92 3 1.81 1.60 1.87 1.43 1.50 1.73 1.01 5 5.23 2.03 1.94 1.32 1.46 5.00 0.95 8 26.31 4.41 2.01 1.30 1.48 24.27 0.99 2 2 1.56 1.57 2.28 2.15 3.11 1.59 2.13 3 2.25 1.76 2.45 2.13 2.37 2.40 1.11 5 7.15 3.21 2.52 1.55 1.94 7.27 2.40 8 37.09 12.17 2.38 1.34 1.77 40.25 5.95 N=1000 0 0 0.21 0.30 0.35 0.23 0.27 0.22 0.18 2 2 0.25 0.30 0.38 0.29 0.28 0.27 0.21 3 0.36 0.31 0.35 0.25 0.25 0.40 0.21 5 1.23 0.37 0.36 0.23 0.26 1.28 0.30 8 6.47 0.72 0.35 0.20 0.24 6.69 1.04 4 2 0.29 0.31 0.39 0.39 0.34 0.29 0.23 3 0.46 0.33 0.41 0.38 0.26 0.47 0.30 5 1.79 0.56 0.42 0.25 0.27 1.96 0.79 8 10.40 1.95 0.46 0.24 0.27 10.57 4.48 6 2 0.32 0.32 0.48 0.59 0.49 0.33 0.33 3 0.51 0.38 0.49 0.49 0.35 0.55 0.63 5 2.03 0.87 0.47 0.26 0.27 2.29 2.31 8 11.68 4.01 0.51 0.25 0.28 13.25 15.74 10 2 0.46 0.34 0.67 1.16 1.16 0.40 0.55 3 0.66 0.50 0.71 0.88 0.58 0.71 1.79 5 2.28 1.87 0.72 0.32 0.38 3.07 11.23 8 12.36 10.57 0.77 0.27 0.32 18.10 69.41 N=2500 0 0 0.08 0.12 0.13 0.08 0.10 0.09 0.08 5 2 0.11 0.12 0.15 0.14 0.11 0.12 0.09 3 0.17 0.13 0.14 0.11 0.10 0.18 0.14 5 0.70 0.18 0.15 0.10 0.10 0.76 0.36 8 4.08 0.53 0.15 0.09 0.10 4.40 2.14 10 2 0.13 0.13 0.17 0.27 0.18 0.13 0.16 3 0.21 0.15 0.19 0.20 0.12 0.24 0.44 5 0.87 0.37 0.21 0.11 0.11 1.02 2.39 8 4.99 1.76 0.18 0.09 0.10 6.43 12.63
440 I.Axt, R.Fried 1 3 small level shifts are present 𝜎 2 other Tr,ad outperforms 𝜎 2 norm Tr,ad , although the latter makes use of the exact normality assumption. The estimator 𝜎 2 other Tr,ad tends to remove more block-estimates than 𝜎 2 norm Tr,ad in the absence of level shifts, see Tables 3 and 4. Therefore, we also expect that more blocks are trimmed away by 𝜎 2 other Tr,ad if level shifts are present, reducing the risk of including perturbed blocks in the trimmed estimator 𝜎 2 other Tr,ad . The trimmed estimator 𝜎 2 Tr,0.5 with a fixed trimming fraction also yields good results. However, this estimation procedure requires the knowledge of the underlying distribution to compute the finite sample correction factor, see Sect. 3.1. The difference-based estimator 𝜎 2 Diff performs well as long as the jumps are moderately high. Table 9 shows the results for the t5 distribution. The estimation procedures 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad yield the best results in this scenario. In Table10 the simulated MSE is presented when the data is generated from the autoregressive (AR) model with 𝜙=0.5 , i.e., the data is positively correlated. The performance of the difference-based estimator worsens considerably then. This is due to the fact that this estimation procedure makes explicit use of the assumption of uncorrelatedness. While 𝜎 2 Diff underestimates the true variance drastically (resulting in a high MSE value) when no changes in the mean are present, the performance seems to improve slightly when dealing with many high jumps. This can be explained by the fact that the positive bias, which arises from the jumps, compensates for the negative bias which arises from the (incorrect) assumption of uncorrelatedness. The blocks-estimator 𝜎 2 Mean exhibits a similar behaviour, since the block size is small when the number of jumps is high, while correlated data require large block sizes to ensure satisfying results. For dependent data the best results are obtained when using the adaptively trimmed estimators. Since the variance 𝜎2 is underestimated when the data are dependent, the values Tj in (15) and (17) get larger resulting in a higher trimming parameter 𝛼adapt . Therefore, more blocks are trimmed away ensuring that the perturbed ones are not involved in the calculation of the overall estimate. Table 8 (continued) K h 𝜎 2 Mean 𝜎 2 Diff 𝜎 2 Tr,0.5 𝜎 2 norm Tr,ad 𝜎 2 other Tr,ad 𝜎 2 mosum Mean 𝜎 2 mosum W 15 2 0.15 0.13 0.23 0.47 0.33 0.15 0.38 3 0.26 0.19 0.27 0.37 0.15 0.27 1.27 5 1.19 0.68 0.28 0.13 0.11 1.27 8.19 8 7.03 3.81 0.24 0.08 0.11 7.54 53.93 25 2 0.21 0.16 0.47 1.22 0.96 0.19 2.00 3 0.39 0.32 0.51 0.94 0.32 0.38 7.90 5 1.86 1.68 0.54 0.17 0.17 1.82 52.63 8 11.14 10.37 0.48 0.11 0.13 11.16 356.78
447 1 3 On variance estimation undershifts inthemean of 101,070 values of the ordinary sample variance for time series which correspond to pixels without virus adhesion. Since changes in the mean are not expected there, we use these data to get some insight into the typical value range of the variance. The sample variance of the contaminated data (upper panel) is 1.59 ×10−4 which is not within the typical range of values, since it exceeds the upper whisker of the boxplot. The other estimation procedures discussed in this paper yield values within the interval [1.1 ×10−4, 1.2 ×10−4] which are well within the interquartile range. We conclude that these approaches yield reasonable estimates for these data. 5.3 PAMONO data withtrend Again, we consider a PAMONO dataset, see Sect.5.2. Panel (a) of Fig.6 shows a time series corresponding to a pixel, which seems to exhibit a virus adhesion as well as a linear trend. Panel (b) of Fig.6 shows the differenced data, i.e., Yt−Yt−1,t=2, …, 388 . The differences of first order appear to be independent and scattered around a fixed mean. Few large differences can be observed which presumably originate from the jumps in the mean at the corresponding time points. The existence of the trend could be explained by the fact that the surface, on which the fluid for virus adhesion is placed, was heated up over time. N=388 observations are available. We apply the estimation procedures 𝜎 2 Mean [using K∈{1, …,5} in the formula (8)], 𝜎 2 Diff , 𝜎 2 Tr,0.5 , 𝜎 2 norm Tr,ad , 𝜎 2 other Tr,ad , 𝜎 2 mosum Mean and 𝜎 2 mosum W to the data and get estimated 0040030020010 0.5500.555 0.5600.565 0.5700.575 Time Intensity (a) 0040030020010 −0.004 0.002 Time (b) Fig. 6 a Intensity over time for one pixel and b corresponding differenced series
448 I.Axt, R.Fried 1 3 values for the variance, which range from 𝜎 2 Mean =0.95 ×10− 6 (with K=5 ) to 𝜎 2 mosum W =1.66 ×10− 6 . The empirical variance of the observations has the value 26.48 ×10−6 , which is much larger than the other estimates. According to our experience the PAMONO data can be assumed to be uncorrelated after differencing. The sample variance of differenced data is 1.93 ×10−6 , which is an estimator of 2𝜎2 , yielding the value 0.97 ×10−6 as an estimate for 𝜎2 , which is near the estimated value of 𝜎 2 Mean . We conclude that the proposed procedures yield reasonable results even in this situation, where a linear trend is present. 6 Conclusion In the presence of level shifts, ordinary variance estimators like the empirical variance perform poorly. In this paper, we considered several estimation procedures in order to account for possible changes in the mean. Estimation of 𝜎2 based on pairwise differences is popular in nonparametric regression and works well in the presence of level shifts and an unknown error distribution if the data are independent and the fraction of shifts is asymptotically negligible. However, we have identified scenarios where estimation based on longer blocks is to be preferred. If only a few small level shifts are expected in a long sequence of observations our recommendation is to use the mean of the blocks-variances 𝜎 2 Mean . This estimation procedure does not require knowledge of the underlying distribution, performs well in the aforementioned situation and is asymptotically even as efficient as the ordinary sample variance if there are no level shifts. If many or large level shifts are expected to occur we recommend using the adaptive trimmed estimators 𝜎 2 norm Tr,ad and 𝜎 2 other Tr,ad . These procedures are constructed for independent data and use either the exact 𝜒2 -distribution or the asymptotic normal distribution of the blockwise estimates, where the second and the fourth moments need to be estimated. We have found these trimming approaches to work reasonably well even under moderate autocorrelations, although many blocks are trimmed away then, presumably due to the underestimation of the unknown variance in the formula (17). Therefore, when no changes in the mean are present the trimmed estimators suffer efficiency loss. On the other hand, we expect that many perturbed blocks are trimmed away in the presence of level shifts reducing the bias of the estimator. The trimming approach could be extended to dependent data in future work. In many applications we rather wish to estimate the standard deviation 𝜎 , e.g. for standardization. If only few jumps of moderate heights are expected to occur, either the average value of the blockwise standard deviations or the square root of the blocks-variance estimator 𝜎 2 Mean can be used. Otherwise, the square root of the trimmed estimator 𝜎 2 other Tr,ad can be recommended. For a large sample size N the finite sample correction factors can be neglected with little loss, see “Appendix”. An interesting extension will be to consider situations where not only the level but also the variability of the data can change. Suitable approaches for such scenarios
449 1 3 On variance estimation undershifts inthemean might be constructed by combining the ideas discussed here with those presented by Wornowizki etal. (2017), where tests for changes in variability have been investigated using blockwise approaches, assuming a constant mean. This will be an issue for future work. Acknowledgements Open Access funding provided by Projekt DEAL. This work has been supported by the Collaborative Research Center “Statistical modelling of nonlinear dynamic processes” (SFB 823) of the German Research Foundation (DFG), which is gratefully acknowledged. The authors would like to thank Sermad Abbas for providing the R+code to extract the PAMONO time series and the referees for their constructive comments which lead to substantial improvements. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/. Appendix Blockwise estimation ofthestandard deviation In many applications, we do not wish to estimate the variance 𝜎2 but rather the standard deviation 𝜎 , e.g. for standardization. Estimation bytheblockwise average We will consider the two blocks-estimators where CN,1 and CN,2 are sample dependent correction factors to ensure unbiasedness when no changes in the mean are present. The block size n can be chosen according to the rule (8). If the number of change-points is not known in practice, it can be estimated as is done in Sect.2.4. For normally distributed data, the correction factors CN,1 and CN,2 can be determined analytically. To derive the correction factor CN,1 for the estimator 𝜎 corr Mean,1, we will first consider the exact distribution of the empirical variance when dealing with jumps in the mean in order to derive the distribution of 𝜎 Mean,1 in (18). Given (18) 𝜎 corr Mean,1 =CN,1𝜎 Mean,1 =CN,1 1 m m ∑ j=1 Sj and (19) 𝜎 corr Mean,2 =CN,2𝜎 Mean,2 =CN,2 √ √ √ √ 1 m m ∑ j=1 S2 j =CN,2 √ 𝜎2Mean ,
450 I.Axt, R.Fried 1 3 independent identically normally distributed data Xj,1,…,Xj,n it is well known that ( n−1)S 2 j 𝜎 2∼𝜒2 n−1 in a jth jump-free block with n observations. The situation is different in the presence of jumps. Without loss of generality the following lemma is expressed in terms of the first block consisting of the observation times t=1, …,n and containing K1 ≤K jumps. Lemma 9 Assume that X1,…,Xn∼ N (0, 𝜎 2) and Y t=Xt+ ∑ K1 k=1 hkIt≥t k for t=1, …,n . Then we have for S 2 1 = 1 n−1∑n t=1 (Yt−Y 1 ) 2 that where 𝜆 1=1 𝜎 2 ∑ n t=1� 𝜇1,t−𝜇1 � 2,𝜇1,t= ∑ K1 k=1 hkIt≥t k and 𝜇1 = 1 n∑ n t=1∑ K1 k=1 hkIt≥t k . Proof Y1=X1+𝜇1 and Yt −Y 1 =X t −X 1 −𝜇 1 +𝜇 1,t ,t=1, …,n , are independent, since X1 and Xt −X 1 are independent and the remaining terms are deterministic constants. Hence, S2 1 and Y1 are independent. Furthermore, with ∑ n t=1 � Yt−𝜇1 𝜎�2 ∼𝜒2 n,𝜆 1 , since Yt∀t are independent and n 𝜎2 X 2 1 ∼𝜒 2 1 . The moment-generating function at z∈ℝ of both sides and the independence of S2 1 and Y1 yield: In the following, we assume that B≤K blocks are contaminated by K1,…, KB jumps, respectively, with ∑B k=1 K k = K . Without loss of generality assume that the n−1 𝜎 2S2 1∼𝜒2 n−1,𝜆1 (the non-central chi-squared distribution) , n ∑ t=1 ( Yt−𝜇1 𝜎 ) 2 = n ∑ t=1 ( Yt−Y1+Y1−𝜇1 𝜎 )2 = n ∑ t=1(Yt−Y1 𝜎)2 + n ∑ t=1(Y1−𝜇1 𝜎)2 +2(Y1−𝜇1 𝜎)n ∑ t=1(Yt−Y1 𝜎) = n−1 𝜎 2S2 1+ n 𝜎 2 ( Y1−𝜇1 ) 2 +0= n−1 𝜎 2S2 1+ n 𝜎 2X 2 1 ( 1−2⋅z)−n∕2exp (𝜆 1z 1−2z ) =M𝜒2 n,𝜆1 (z)=Mn−1 𝜎2S2 1 (z)⋅M𝜒2 1(z) =Mn−1 𝜎2S2 1 (z)⋅(1−2⋅z)−1∕2 ⇔Mn−1 𝜎2S2 1 (z)=(1−2⋅z)−(n−1)∕ 2 exp (𝜆1z 1−2z)=M𝜒2 n−1,𝜆1 (z ) ⇒n−1 𝜎 2S2 1∼𝜒2 n−1,𝜆1.
451 1 3 On variance estimation undershifts inthemean jumps are contained in the first B blocks, while the last m−B>0 blocks do not contain any jumps. The square root of a 𝜒2 n−1,𝜆 j -distributed random variable ( n−1)S 2 j ∕ 𝜎2 is 𝜒 -distributed with n−1 degrees of freedom and non-centrality parameter √𝜆j , see e.g. Miller (1964). We hence have √ n−1Sj∕𝜎∼𝜒n−1, √ 𝜆j , j=1, …,m , where 𝜆j=0 for the last blocks j=B+1, …,m , i.e., √ n−1S j ∕𝜎∼𝜒 n−1 . The expected value of Sj is given as where F1,1(a,b,z) represents the generalized hypergeometric function, see Olver etal. (2010) for more details. When no changes in the mean are present we have that 𝜆j=0∀j and therefore F1,1(− 0.5, 0.5(n−1),−0.5𝜆j)=1 . The exact finite sample correction factor is given as which is the reciprocal of the term C n,𝜆 j in (20) when no level shifts are present, since F1,1(− 0.5, 0.5(n−1),−0.5𝜆j)=1 , j=1, …,m , in this case. For the second estimator (19), we have the following statements on its expectation and a suitable finite sample correction factor: We used the fact that 𝜎 Mean,2 follows a scaled 𝜒u,v distribution with u=m(n−1) degrees of freedom and the non-centrality parameter v = �∑ B j=1𝜆j , since n−1 𝜎 2 ∑m j=B+1 S2 j ∼𝜒2 (m−B)(n−1) and n−1 𝜎2 ∑B j=1S2 j∼𝜒2 B(n−1), ∑ B j=1 𝜆 j . Using this information we can determine the expectation of the estimator straightforwardly, see Miller (1964) and Olver etal. (2010). The correction factor is the reciprocal of Dn,𝜆1,…,𝜆B, where we have F 1,1 � −0.5, 0.5m(n−1),−0.5 ∑ B j=1𝜆j � = 1 in the absence of level shifts. (20) E (Sj)=𝜎 √ 2 √n−1 Γ(0.5n) Γ(0.5(n−1)) F1,1(− 0.5, 0.5(n−1),−0.5𝜆j) =∶ 𝜎Cn,𝜆j , C N,1 = √ n−1 √2 Γ(0.5(n−1)) Γ(0.5n) , 𝜎 corr Mean,2 =CN,2 𝜎 √m(n−1) � n−1 𝜎2 m � j=B+1 S2 j+n−1 𝜎2 B � j=1 S2 j �1∕2 , E �𝜎 corr Mean,2�=CN,2 𝜎√2 √m(n−1) Γ(0.5(m(n−1)+1)) Γ(0.5m(n−1)) ×F1,1�−0.5, 0.5m(n−1),−0.5 B � j=1 𝜆j� ∶= CN,2𝜎Dn,𝜆1,…,𝜆B, CN,2 = √ m(n−1) √2 Γ(0.5m(n−1)) Γ(0.5(m(n−1)+1)).
452 I.Axt, R.Fried 1 3 The following consistency statements are valid for the two introduced uncorrected estimators 𝜎 Mean,1 = 1 m∑m j=1 S j [as defined in (18)] and 𝜎 Mean,2 = � 1 m ∑ m j=1S2 j [as defined in (19)] of 𝜎 : Corollary 10 Under the conditions of Theorem1 the estimators 𝜎 Mean,1 and 𝜎 Mean,2 converge almost surely to 𝜎 , as N→∞ . Proof The strong consistency of 𝜎 Mean,2 follows immediately from the Continuous Mapping Theorem. For 𝜎 Mean,1 , we have due to Theorem2 of Hu etal. (1989) that since Sj −E ( S j) are uniformly bounded with P ( | S j −E ( S j)| >t)→0∀ t due to Chebyshev’s inequality and Var(Sj) → 0 . Let Sj,h be the sample standard deviation in the perturbed block while Sj,0 is the estimate in the uncontaminated block. We have that i.e., it suffices to show 1 m�∑ m j=B+1E � Sj,0 � + ∑ B j=1E � Sj,h �� → 𝜎 . For the first of these two terms, we have since the consistency and the decreasing variance of S1,0 implies convergence of the expectation, see Lemma 1.4A in Serfling (1980). Using Jensen’s inequality we get for the second term 1 m m ∑ j=1 ( Sj−E ( Sj )) → 0 almost surely, 1 m m ∑ j=1 ( Sj−E ( Sj )) =𝜎 Mean,1 −1 m (m ∑ j=B+1 E ( Sj,0 ) + B ∑ j=1 E ( Sj,h )), 1 m m ∑ j=B+1 E ( Sj,0 ) =m−B mE ( S1,0 ) ⟶𝜎as N⟶∞ , 1 m B ∑ j=1 E(Sj,h)=1 m B ∑ j=1 E (√ S2 j,h ) ≤1 m B ∑ j=1√ √ √ √E(S2 j,0 + n ∑ t=1 2(Xj,t−Xj)(𝜇j,t−𝜇j)+(𝜇j,t−𝜇j)2 n−1) =1 m B ∑ j=1√ √ √ √E(S2 j,0)+ n ∑ t=1(𝜇j,t−𝜇j)2 n−1=1 m B ∑ j=1√ √ √ √𝜎2+ n ∑ t=1(𝜇j,t−𝜇j)2 n−1 ≤B m√ √ √ √ √ 𝜎2+n n−1(K ∑ k=1 hk)2 ⟶0,
453 1 3 On variance estimation undershifts inthemean where 𝜇j,t and 𝜇j are defined in the proof of Theorem1. Remark 11 The correction factors CN,1 and CN,2 from (18) and (19) satisfy where CN,1 = 𝜎 ∕E( 𝜎 Mean,1) and CN,2 = 𝜎 ∕E( 𝜎 Mean,2) in the absence of level shifts. This can be shown with Lemma 1.4A in Serfling (1980), since 𝜎 Mean,1 and 𝜎 Mean,2 are consistent estimators and their variances tend to zero which implies convergence of the means and thus the above statement. Therefore, for large N and n we can neglect the correction factors and use the estimators 𝜎 Mean,1 and 𝜎 Mean,2 instead of 𝜎 corr Mean,1 and 𝜎 corr Mean,1 with block sizes n→∞ . Trimmed estimation When dealing with a large number of level shifts, as is discussed in Sect. 3, the square root of the variance estimator 𝜎 2 Tr,ad from (13) can be used to estimate the standard deviation 𝜎 . For large N and n, a correction factor to ensure unbiasedness when no changes in the mean are present can be neglected. Table12 shows the simulated finite sample correction factors for normally and t5 -distributed data as well as for the stationary AR(1)-process with normal errors and parameter 𝜙∈{0.3, 0.6} . We observe that the correction factors are nearly one except for strongly correlated data, i.e., AR-process with parameter 𝜙=0.6 . See Figs.7, 8, 9, 10, 11 and 12. CN,1 → 1 and CN,2 → 1 as N → ∞, Table 12 Simulated finite sample correction factors for the adaptively trimmed estimation procedures for normally and t5 -distributed data as well as for the stationary AR(1)-process with normal errors and parameter 𝜙∈{0.3, 0.6} , denoted by AR(0.3) and AR(0.6) N(0,1) t5 AR(0.3) AR(0.6) √ 𝜎2norm Tr,ad N = 1000, n = 50 1.0025 1.0660 1.0213 1.0939 N = 5000, n = 100 0.9999 1.0413 1.0108 1.0450 √ 𝜎2other Tr,ad N = 1000, n = 50 1.0082 1.0666 1.0310 1.1179 N = 5000, n= 100 1.0014 1.0334 1.0135 1.0550
454 I.Axt, R.Fried 1 3 012345 02040608 0 δ MSE−optimal block length (a) 012345 0.022 0.0240.026 0.028 δ MSE (b) 012345 0.022 0.0240.026 0.028 δ MSE (c) Fig. 7 a MSE-optimal block length nopt of 𝜎 2 Mean , b MSE regarding nopt of 𝜎 2 Mean and c MSE of 𝜎 2 Mean when choosing n = √ N K+1 for K=1 (—), K=3 (- - -) and K=5 ( ⋅⋅⋅ ) with N=1000 , Y t=Xt+ ∑K k =1 hIt≥t k , where Xt∼t5 and h=𝛿 ⋅ 𝜎,𝛿∈{0, 0.1, …,5} 012345 0.0000.005 0.0100.015 0.0200.025 0.030 δ MSE (a) 012345 0.0200.025 0.0300.035 0.0400.045 0.050 δ MSE (b) Fig. 8 MSE of 𝜎 2 Mean when choosing n = √ N K+1 for true K=5 (—), K=0 ( ) K=1 ( --- ) K=2 ( ···· ) K=3 ( -·- ) K=4 ( ––– ) and K=6 ( –-– ) with N=1000 and h=𝛿 ⋅ 𝜎,𝛿∈{0, 0.1, …,5} , Y t=Xt+ ∑5 k =1 hIt≥t k , where a Xt∼N(0, 1) and b Xt∼t3
455 1 3 On variance estimation undershifts inthemean 012345 02040608 0 δ MSE−optimal block length (a) 012345 0.0008 0.0012 0.0016 δ MSE (b) 012345 0.0008 0.0012 0.0016 δ MSE (c) Fig. 9 a MSE-optimal block length nopt of 𝜎 2 Mean , b MSE regarding nopt of 𝜎 2 Mean and c MSE of 𝜎 2 Mean when choosing n = √ N K+1 for K=1 (—), K=3 (- - -) and K=5 ( ⋅⋅⋅ ) with N=2500 , Y t=Xt+ ∑K k =1 hIt≥t k , where Xt∼N(0, 1) and h=𝛿 ⋅ 𝜎,𝛿∈{0, 0.1, …,5} 012345 02040608 0 δ MSE−optimal block length (a) 012345 0.0090.011 0.013 δ MSE (b) 012345 0.0090.011 0.013 δ MSE (c) Fig. 10 a MSE-optimal block length nopt of 𝜎 2 Mean , b MSE regarding nopt of 𝜎 2 Mean and c MSE of 𝜎 2 Mean when choosing n = √ N K+1 for K=1 (—), K=3 (- - -) and K=5 ( ⋅⋅⋅ ) with N=2500 , Y t=Xt+ ∑K k =1 hIt≥t k , where Xt∼t5 and h=𝛿 ⋅ 𝜎,𝛿∈{0, 0.1, …,5}
456 I.Axt, R.Fried 1 3 References Angelova, J.A.: On moments of sample mean and variance. Int. J. Pure Appl. Math. 79(1), 67–85 (2012) Dai, W., Tong, T.: Variance estimation in nonparametric regression with jump discontinuities. J. Appl. Stat. 41(3), 530–545 (2014) Dai, W., Ma, Y., Tong, T., Zhu, L.: Difference-based variance estimation in nonparametric regression with repeated measurement data. J. Stat. Plan. Inference 163, 1–20 (2015) Davies, L., Gather, U.: The identification of multiple outliers. J. Am. Stat. Assoc. 88(423), 782–792 (1993) −2 −1 012 5000 6000 7000 8000 9000 10000 11000 Theoretical Quantiles Sample Quantiles Fig. 11 Q–Q plot of the maximal monthly discharge of the Nile river in the period 1871–1984 −3 −2 −1 0123 0.28 0.29 0.30 0.31 0.32 0.33 Theoretical Quantiles Sample Quantiles Fig. 12 Q–Q plot of a pixel with a virus adhesion from the PAMONO data