We derive explicit expressions for the coefficients of certain quadratic
forms arising in the singular-spectrum analysis, a novel technique of
time series analysis. We also derive expressions for the first two moments of these quadratic forms under the assumption that the time series is a sum of a linear trend and white noise

V. G. Moskvina

Hrčak - Portal of scientific journals of Croatia

Mathematical Communications 5(2000), 161-171 161Distribution of random quadratic forms arising insingular-spectrum analysisValentinaG.Moskvina∗Abstract. We derive explicit expressions for the coefficients of cer-tain quadratic forms arising in the singular-spectrum analysis, a noveltechnique of time series analysis. We also derive expressions for thefirst two moments of these quadratic forms under the assumption thatthe time series is a sum of a linear trend and white noise.Key words: time series analysis, random quadratic forms, singular-spectrum analysis, trajectory matrix, centringAMS subject classifications: 62M10Received January 8, 2001 Accepted January 25, 20011. IntroductionOne of the methods of time series analysis, that has recently attracted a lot ofattention, is the so-called singular-spectrum analysis (SSA), see Elsher and Tsonis(1996) for an introduction and Golyandina, Nekrutkin and Zhigljavsky (2001) fora detailed description of the methodology and theory.Let y1, y2, . . . be a time series. SSA is based on the singular-value decompositionof the so-called “trajectory matrix”X = (xij)M,Kij=1 =y1 y2 y3 . . . yKy2 y3 y4 . . . yK+1.......... . ....yM yM+1 yM+2 . . . yN (1)where N is the number of observations, M is an integer called “lag parameter”(M ≤ N/2), K = N −M + 1 and xij = yi+j−1 for all i = 1 . . . ,M, j = 1 . . . ,K.In many cases performing a singe or double centring of the trajectory matrixis recommended, see Golyandina, Nekrutkin and Zhigljavsky (2001). Denote by Zthe matrix obtained from X after a centring operation is performed.Single centring means that the row averages are subtracted from each elementof the matrix X; that is, the elements of Z arezij = xij − 1MM∑i=1xij . (2)∗School of Mathematics, Cardiff University, P.O. Box 926, Cardiff CF24 4YH, United Kingdom,e-mail: MoskvinaV@Cardiff.ac.uk162 V.G.MoskvinaIf we subtract the row and column averages from xij and add the general mean,then we obtain double centring of the trajectory matrix. The resultedM×K-matrixZ = (zij)M,Kij=1 has the elementszij = xij − 1MM∑i=1xij − 1KK∑j=1xij +1MKM∑i=1K∑j=1xij . (3)Obviously, we could use the formulae (2) and (3) to compute the values zij forj > K; this corresponds to the use of yt with t > N .The following sum of squaresSSp,q =M∑i=1q∑j=p+1z2ij (4)is of high importance in the theory of SSA in relation to the problems of de-tection of structural heterogeneities in time series, see Chapter 3 in Golyandina,Nekrutkin and Zhigljavsky (2001), Moskvina and Zhigljavsky (2000) and the Web-site http://www.cf.ac.uk/maths/stats/changepoint/ . Here p and q are arbitrary in-tegers, 0 ≤ p < q.The case when p = 0 and q = K is of a particular interest. Another importantcase is when p ≥ K; in this case the sum of squares is computed for the subseriesthat contain points lying outside the interval [1, N ].We will demonstrate that the sum of squares SSp,q can be represented as aquadratic formSSp,q = eTBe, (5)where e = (e1, . . . , eq+M−1)T is some vector that can often be associated with theobservation noise in the original series y1, y2, . . . and B is a (q+M−1)× (q+M−1)matrix.In Section 2 we consider the most difficult case of double centring. We assumethat yt = a+ bt+ et (t = 1, 2, . . .), where a and b are unknown parameters and etis noise.If b is known (and we thus can assume b = 0), then there is no need to do doublecentring of the trajectory matrix X to achieve independence of a, the subtractionof the row averages would suffice. This case is considered in Section 3.In the simplest case, when both a and b are known (we thus can assume a =b = 0), we do not need to do centring at all to achieve independence of a and b.We thus have zij = xij = yi+j−1 in the formula for SSp,q; the case of no centringis considered in Section 4.In all three cases we shall compute the matrix B and the first two moments underthe assumption that et are independent normal random variables with Eet = 0 andvar(et) = σ2 = const.To compute the moments, we will use the results on distributions of quadraticforms of random variables, see e.g. Section 3.2 in Mathai and Provost (1992) andSection 2.5 in Searle (1971). These results imply thatE(eTBe) = σ2trB , (6)Distribution of random quadratic forms 163var(eTBe) = 2σ4trB2 (7)Note that the proof of (6) does not require normality of et and that both properties(6) and (7) can be easily generalized to the case when the components of e aredependent (e.g. to the case when et is an autoregressive process).2. Double centring2.1. NotationAssume that yt = a+ bt+et (t = 1, 2, . . .), where a and b are unknown parametersand et is noise. Consider the case of double centring, when the elements of thematrix Z are defined by (3).An important observation about the matrix Z is that its elements zij do notdepend on the actual values of the parameters a and b. Indeed, we can writexij = a+b(i+j−1)+εij, where εij = ei+j−1, and therefore we obtainzij = a+ b(i+ j − 1) + εij + 1KMM∑i=1K∑j=1(a+ b(i+ j − 1) + εij)− 1MM∑i=1(a+ b(i+ j − 1) + εij)− 1KK∑j=1(a+ b(i+ j − 1) + εij)= εij − ε.j − εi. + ε.. ,whereεi. =1KK∑k=1εik, ε.j =1MM∑i=1εij , ε.. =1MKM∑i=1K∑k=1εik (8)Denote Q = q−p. To formulate the results, we shall need the function wM,p,q(t)which is defined as follows (see also Figure 1):if M ≤ Q,wM,p,q(t) =t− p, for p < t ≤ p+MM, for p+M < t ≤ qq +M − t, for q < t ≤ q +M0, otherwise(9)and, if 0 < Q ≤M ,wM,p,q(t) =t− p, for p < t ≤ qQ, for q < t ≤M + pq +M − t, for M + p < t ≤M + q0, otherwise.(10)164 V.G.MoskvinaFigure 1. Function wM,p,q(t)We shall also need six matrices B1, . . . , B6 that are defined as follows.Matrix B1 is the diagonal matrix of size (M+Q−1)× (M+Q−1); its elements areb(1)ij ={wM,0,Q(i), if i = j0, otherwise. (11)Matrix B2 is presented in Figure 2. It is an (M+Q−1)× (M+Q−1)-matrix withA=min{M,Q}. The elements of B2 areb(2)ij =wi,0,A(j), 1 ≤ i < AwA,0,A(j), A ≤ i < max{M,Q}wM+Q−i,i−A,i(j), max{M,Q} ≤ i ≤M+Q−1.(12)Matrix B3 is the matrix V of Figure 2; its size is (M+K−1)× (M+K−1) and theparameters are A=M, v1=v2= K and a =M−1. The elements of B3 areb(3)ij =wi,0,K(j), if 1 ≤ i < MwM,0,K(j), if M ≤ i < KwM+K−i,i−K,i(j), if K ≤ i ≤M +K − 1.(13)Matrix B4 is the matrix U of Figure 2 with A = B = M and u1=u2=K−M+1.Its size is (M+K−1)× (M+K−1) and its elements areb(4)ij = wM,0,K(i)× wM,0,K(j), (14)where 1 ≤ i, j ≤ K+M−1.Distribution of random quadratic forms 1651 2 … B … B … 2 1 1 1 … 1 12 4 … 2B … 2B … 4 2 1 2 … 2 1 a……… … …… … …=A 2A … AB … AB … 2A A 1 2 … A … A … 2 1U =… … … …  … … … … … … …v1u1A 2A … AB … AB 2A A 1 2 … A … A … 2 1…… … ……… …… … …A-12 4 … 2B … 2B … 4 2 a 1 2 … 2 11 2 … B … B … 2 1 1 1 … 1 1u2 a v21 1 … 11 2 .. 2 1B-1… …… … …B 2 =1 2 … A … 2 1… … ……… …1 2 … A … 2 1A… … …… …1 2 … 2 11 … 1 1A0000Figure 2. Matrices U , V and B2.Matrix B5 is an (M+K−1) × (M+Q−1)-matrix V (see Figure 2) with A =min{M,Q}, v1==K, v1==Q, and a =M−1. The elements areB5 = b(5)ij =wi,0,Q(j), if 1 ≤ i ≤MwM,0,Q(j), if M ≤ i ≤ KwM+K−i,i−K,Q+i−K (j), if K ≤ i ≤M +K − 1.Matrix B6 is the matrix U with A = M, B = min{M,Q}, u1 = K−M+1 andu2 =| Q−M | +1. Its size is (M +K−1) × (M +Q−1) and the elements areb(6)ij = wM,0,K(i)× wM,0,Q(j) for 1≤ i≤K+M−1, 1≤j≤Q+M−1.Finally, as an addition to (8), we defineε˜i. =1Qq∑j=p+1εij , ε˜.. =1MQM∑i=1q∑j=p+1εij (15)2.2. The quadratic formTheorem 1. Let yt = a+ bt+ et for all t = 1, 2, . . . and some a and b. Then thesum of squares SSp,q can be represented as a quadratic formSSp,q = eTBe ,166 V.G.Moskvinawhere e = (e1, . . . , eq+M−1)T andB =(0 00 B1)− 1M(0 00 B2)+QK2(B3 00 0)− QK2M(B4 00 0)− 1K[(0 B50 0)+(0 0BT5 0)]+1KM[(0 B60 0)+(0 0BT6 0)]The structure of the matrix B is represented in Figure 3, where zero is codedby white and different patterns in right bottom, left top, and left bottom cornerscorrespond to the matrices B1 − 1MB2, QK2B3 − QK2MB4 and − 1KB5 + 1KMB6,respectively.Q+M-1M-1M-1MM-1M-1 MM-1 MMK-M+1p-TT Q+M-1Q-M-1M-1M-1Tp-TM-1Q-M+1Q-M+1K-M+1Q-M-1Figure 3. Matrix BProof of Theorem 1. According to the notations of (8) and (15) and the definitionof the sum of squares (4), we haveSSp,q =M∑i=1q∑j=p+1(εij − ε.j − εi. + ε..)2 .Distribution of random quadratic forms 167Expanding the square and simplifying separately each term in the sum, we obtainSSp,q = S1 −MS2 +QS3 −MQS4 − 2QS5 + 2MQS6 , (16)whereS1 =M∑i=1q∑j=p+1ε2ij , S2 =q∑j=p+1ε2.j , S3 =M∑i=1ε2i. ,S4 = ε2.. , S5 =M∑i=1ε˜i.εi. , S6 = ε˜..ε..To see how the simplification is made, consider the three most difficult terms:M∑i=1q∑j=p+1εi.ε.j=M∑i=1εi.q∑j=p+1ε.j= M∑i=11KK∑j=1εij q∑j=p+11MM∑i=1εij=Mε..Qε˜.. ,M∑i=1q∑j=p+1εijε.j =Mq∑j=p+1ε.j[1MM∑i=1εij]=Mq∑j=p+1ε2.j ,M∑i=1q∑j=p+1εijεi. = QM∑i=1εi. 1Qq∑j=p+1εij = Q M∑i=1εi.ε˜i. .Consider separately each sum in (16).It is easy to see from (11) and (16) thatS1 = eT(0 00 B1)e .Consider S2 and replace ε.j by the corresponding expression from (8). Using thenotation εij = ei+j−1 and e = (ep+1 . . . , eq+M−1)T we haveS2 =q∑j=p+1(1MM∑i=1εij)2 =1M2q∑j=p+1M∑i=1M∑i′=1ei+j−1ei′+j−1 .We now change the indices by the formula{t = i+ j − 1t′ = i′ + j − 1Since i, i′ = 1, . . . ,M , j = p+1, . . . , q we have t, t′ = p+1, . . . , q+M−1. Therefore,S2 =1M2p+M∑t=p+1p+M∑t′=p+1etet′ +p+M+1∑t=p+2p+M+1∑t′=p+2etet′ + . . .+q+M−1∑t=qq+M−1∑t′=qetet′︸ ︷︷ ︸q − p168 V.G.Moskvina=1M2e( 1M 0M,Q−10Q−1,M 0Q−1,Q−1)+ 01,1 01,M 01,Q−20M,1 1M 0M,Q−20Q−2,1 0Q−2,M 0Q−2,Q−2+ . . .+ 0Q−2,Q−2 0Q−2,M 0Q−2,10M,Q−2 1M 0M,101,Q−2 01,M 01,1+ ( 0Q−1,Q−1 0Q−1,M0M,Q−1 1M) eTwhere 1M is the M × M -matrix of ones and 0u,v is the u × v-matrix of zeros.Summing the matrices above we obtain the matrix B2, see (12), andMS2 =1MeT(0 00 B2)eMatrix B3, see (13), is built from the sum of squares S3 analogously to theconstruction of B2 from S2. In this case we use the substitution{t = i+ k − 1t′ = i+ k′ − 1where i = 1, . . . ,M, k, k′ = 1, . . . ,K and t, t′ = 1, . . . ,K+M−1. This givesQS3 =QK2eT(B3 00 0)e .The sum of squares S4 gives us the matrix 1K2M2B4; this can be easily seen ifwe consider B4 as a square of the (M+K−1)× (M+K−1)-matrix with elementswM,0,K(i), if i = j, and 0, otherwise. Thus,MQS4 =QK2MeT(B4 00 0)eWe can obtain the matrix B5 similarly to the case of B3 using the similaritybetween the sums of squares S3 and S5. In view of the differences in the numberof terms in the sums ε˜i. and εi., see (8) and (15), we obtain different sizes andmultipliers. We thus have2QS5 =1KeT[(0 B50 0)+(0 0BT5 0)]e .Finally, B6 is an (M+K−1)×(M+Q−1)-matrix which can be obtained similarlyto B4.2MQS6 =1KMeT[(0 B60 0)+(0 0BT6 0)]e .This completes the proof. ✷2.3. Computation of trB and trB2As a corollary of Theorem 1, we can derive the expressions for trB and trB2 inseveral particular cases, depending on the values of parameters p and q.Distribution of random quadratic forms 169Case 1: p ≥ N . We havetrB =Q(M − 1)(3K2 +M + 1)3K2trB2 =(2K2+Q2)(M2−1)(2M2+7)45K4+Q(M−1)(M2−M−1)3M− (M−1)(2M2−M−3)6.Case 2: p = 0, q = K. In this caseSSp,q =M∑i=1K∑k=1(εik − ε.k − εi. + ε..)2 ==M∑i=1K∑k=1ε2ik −KM∑i=1ε2i. −MK∑k=1ε2.k +MKε2.. = eTBeimplying Q = K, B3 = B5 = BT5 and B4 = B6 = BT6 ; therefore, the matrix B hasthe formB = B1 − 1MB2 +1KB3 − 1MKB4 .For i ≤ j and i+ j ≤ K +M we havebij =i(i+MK−M−K)MK , if i = j, 1 ≤ i ≤M − 1M − 1, if i = j,M ≤ i ≤ M+K2j−iM − 1, if M ≤ j, j −M + 1 ≤ i ≤ j− i(M+K−j)MK , if 1 ≤ i < j ≤M − 1i(M+K−j)MK , if K + 1 = K +M − 1, 1 ≤ i ≤ j −K− (M−i)(j−K)MK , if K + 1 ≤ j = K + i, 1 ≤ i ≤M − 10, otherwiseand bij = bji, bM+K−i,M+K−j = bij for 1 ≤ i ≤ M +K − 1, 1 ≤ j ≤M +K − 1.This yieldstrB =(M − 1)(3K2 −M − 1)3K,trB2 =(M2−1)(2M2+7)45K2+K(M−1)(3M2−M−1)3M− (M2−1)(2M−3)6+(M2−1)(2−3M2)15KM.3. Single centringConsider the case of single centring. Using the notation εij = ei+j−1 and (8) and(15) we obtainSSp,q =M∑i=1q∑j=p+1(εij − εi.)2 =M∑i=1q∑j=p+1ε2ij +QM∑i=1ε2i. − 2QM∑i=1ε˜i.εi.This yields that the sum of squares SSp,q can again be represented as a quadraticform eTBe, where the matrix B has the formB =(0 00 B1)+QK2(B3 00 0)− 1K[(0 B50 0)+(0 0BT5 0)].170 V.G.MoskvinaHere the matrices B1, B3, B5 are described as above and 0 is the zero-matrix ofsuitable size.Let us compute trB and trB2 in the two particular cases as in Section 2.Case 1: p ≥ N . We havetrB =MQ(K + 1)K;trB2 = M2Q− M(M2 − 1)3+Q2K4[(K −M + 1)2M2 − M(M − 1)6(5M2 − 8KM − 7M + 4K)]+2K2[(K −M + 1)(Q−M + 1)M2−M(M − 1)6(5M2 − 4M(K + Q)− 7M + 2(K + Q))]Case 2: p = 0, q = K. We haveeTBe = SSp,q =M∑i=1K∑j=1(εij − εi.)2 =M∑i=1K∑j=1ε2ij −KM∑i=1ε2i. ,implying Q = K, B3 = B5 = BT5 . Therefore, the matrix B has the formB = B1 − 1KB3;its elements arebij =K−1K wM,O,K(j), if i = j− 1Kwi, 0, K(j), if i = j, 1 ≤ i ≤M− 1KwM, 0, K(j), if i = j, M ≤ i ≤ K− 1KwM+K−i, i−K, i(j), if i = j K ≤ i ≤M +K − 1.The expressions for trB and trB2 aretrB =M(K − 1), trB2 = M2K − M(M2 − 1)3−M2 + M2(M2 − 1)6K2.4. No centringThe case without averaging is the easiest. In this case we haveSSp,q =M∑i=1q∑j=p+1ε2ij = eTBe, where B =(0 00 B1),trB =∑twM,p,q(t) =∑twM,0,Q(t) =MQDistribution of random quadratic forms 171andtrB2 =∑tw2M,p,q(t) ={13M(3MQ+ 1−M2) , if Q ≥M13Q(3MQ+ 1−Q2) , if Q ≤M .In particular, for p = 0 and q = K we have B = B1,trB =MK and trB2 =M2K − M(M2 − 1)3.References[1] J. Elsner, A.Tsonis, Singular Spectrum Analysis: A New Tool in Time Se-ries Analysis, Plenum Press, New York, 1996.[2] N.Golyandina, V.Nekrutkin, A. Zhigljavsky, Analysis of Time SeriesStructure: SSA and Related Techniques, Chapman and Hall/CRC, 2001.[3] A.M.Mathai, S. B. Provost, Quadratic Forms in Random Variables: The-ory and Applications, Marcel Dekker, New York, 1992.[4] V.Moskvina, A. Zhigljavsky, Application of singular-spectrum analysis tochange-point detection in time series, Journal of Time Series Analysis, submit-ted.[5] S.R. Searle, Linear models, Wiley, New York, 1971.

Distribution of random quadratic forms arising in singular-spectrum analysis

Abstract

Similar works

Full text

Available Versions

Hrčak - Portal of scientific journals of Croatia