A Method of Trend Extraction Using Singular Spectrum Analysis

A METHOD OF TREND EXTRA CTION USING SINGULAR SPECTR UM ANAL YSIS Theo dore Alexandro v Cen ter for Industrial Mathematics, Univ ersity of Bremen, German y (theo dore@math.uni-bremen.de) Octob er 22, 2018 Abstract The pap er presents a new metho d of trend extraction in the frame- work of the Sing ular Sp ectrum Analysis (SSA) approach. This metho d is easy to us e , do es not need sp ecification of models of time s eries a nd trend, allows to extract trend in the presence o f noise and o scillations and ha s o nly tw o parameter s (b esides basic SSA parameter called win- dow length). One par ameter manages sca le of the extr acted trend and another is a metho d s pecific threshold v alue. W e pr op ose pro cedures for the choice of the par ameters. The presented metho d is ev aluated on a simulated time ser ies with a p olyno mial trend and an o s cillating com- po nent with unkno wn per io d and on the seasonally adjusted mont hly data of unemployment level in Alask a for the p erio d 19 76/01 -2006 /09. Keyw ords: Time ser ies; trend extraction; Singular Sp ectrum An alysis. 1 INTR ODUCTION T rend e xtr action is an imp ortan t task in app lied time series analysis, in particular in economics and engineering. W e present a new metho d of trend extraction in the framewo rk of the Singular S p ectrum Analysis approac h. T rend is usually defi ned as a smo oth add itiv e co mp onent conta in in g information ab out time series global c h ange. This defin ition is r ather v ague (whic h t yp e of smo othness is used? whic h kind of in formation is cont ained in the trend?). I t may sound strange, bu t there is no m ore p recise defin ition of the trend accepted by the m a j ority of researc hers and practitioners. Eac h approac h to trend extraction defines trend with resp ect to the m athematica l to ols used (e .g. using F ourier t r ansformation or deriv ativ es). Thus in the 1 corresp onding literature one can fin d v arious sp ecific definitions of the trend . F or furth er discussion on trend issues we refer to [2]. Singular S p ectrum Analysis (SSA) is a g ener al approac h to time series analysis and forecast. Algo r ithm of SSA is similar to that of P r incipal Com- p onents Analysis (PCA) of m ultiv ariate data. I n con trast to P CA which is applied to a matrix, SSA is applied to a time series and pro vides a rep re- sen tation of the giv en time serie s in terms of eig env alues and eigen v ectors of a matrix made of t h e time series. The basic idea of SSA h as b een p ro- p osed by [5] for dimension calculation and r econstruction of attrac tors of dynamical systems, see historica l reviews in [10] and in [11]. In this pap er w e mostly follo w the notations of [11 ]. SSA can b e used for a wide r ange of tasks: trend or quasi-p erio dic com- p onent detection and extractio n , d enoising, forecasting, change- p oin t detec- tion. T he pr esen t b ib liograph y on SSA in cludes tw o m onographes, sev eral b o ok c hapters, and o ver a hundred pap ers. F or more details see references at the website SS Aw iki: ht tp://www.math.u ni-bremen.de/ ∼ theo dore/ssa wiki. The metho d p resen ted in this pap er has b een fir st p rop osed in [3] and is studied in detail in the author’s unpu blished Ph.D. thesis [1] a v ailable only in Russ ian at h ttp://www.p dmi.ras.ru/ ∼ theo/autossa. The prop osed method is easy to u se (h as only t w o parameters), d o es not need s p ecification of mo dels of time series and trend , allo w s one to sp ecify desired trend sca le, and extracts trend in the presence of n oise and oscillati ons . The outline of this p ap er is as follo ws. Section 2 introd u ces SSA, form u- lates p rop erties of trend s in SSA and p resen ts the already existing metho ds of trend extrac tion in SS A. Section 3 prop oses our metho d of trend extrac- tion. In section 4 we discus s the frequen cy prop erties of add itiv e comp onen ts of a time series and p resen t ou r pro cedure for the c h oice of first parameter of the metho d, a lo w-fr equency b oundary . Section 5 starts with in vesti gation of the role of the second metho d parameter, the lo w -fr equency con tribu tion, based on a simulat ion example. Then we p r op ose a h euristic strategy for the c hoice of th is parameter. In section 6, applications of the prop osed m etho d to a sim ulated time series with a p olynomial trend and oscil lations and on the unemplo yment lev el in Alask a are considered. Finally , section 7 offers conclusions. 2 2 SINGULAR SPECTR UM ANAL YSIS Let us ha v e a time series F = ( f 0 , . . . , f N − 1 ), f n ∈ R , of length N , and w e are lo oking for some sp ecific additiv e comp on ent of F (e.g. a trend). The cen tral idea of SSA is to em b ed F in to h igh-dimensional euclidean space, then find a sub s pace corresp ond ing to the sough t-for comp onent and, fi nally , reconstruct a time s eries comp onent c orr esp onding to this su bspace. The c hoice of the su bspace is a crucial question in SSA. The basic SSA algorithm consists of decomp osition of a time series and r econstruction of a desired additiv e comp onent. These tw o steps are summarized b elo w; for a detaile d description, s ee page 16 of [11]. Decomp osition. Th e decomp osition take s a time s er ies of length N and comes up w ith an L × K matrix. This stage starts b y defin ing a parameter L (1 < L < N ), called the window length, a n d constructing the so- called tra jectory matrix X ∈ R L × K , K = N − L + 1, with stepwise tak en p ortions of the orig in al time series F as columns: F = ( f 0 , . . . , f N − 1 ) → X = [ X 1 : . . . : X K ] , X j = ( f j − 1 , . . . , f j + L − 2 ) T . ( 1) Note that X is a Hankel matrix and (1) defi n es one-to -one corresp ondence b et ween series of length N and Hank el matrices of size L × K . Then Sin gular V alue Decomp osition (SVD) of X is app lied, w here j -th comp onen t of SVD is sp ecified b y j -th eigen v alue λ j and eigen vect or U j of XX T : X = d X j =1 p λ j U j V j T , V j = X T U j . p λ j , d = m ax { j : λ j > 0 } . Since the matrix XX T is p ositiv e-definite, their eigen v alues λ j are p ositiv e. The SVD comp onents are num b ered in the d ecreasing order of eigen v alues λ j . W e define j -th Empirical Orthogonal F unction (EO F) as th e sequence of elements of the j -th eigen v ector U j . T he triple ( p λ j , U j , V j ) is called j -th eigen triple, p λ j is called the j -th sin gu lar v alue, U j is the j -th left singular v ector and V j is the j -th righ t sin gular v ector. Reconstruction. Reconstru ction g o es fr om an L × K matrix in to a time series of le n gth N . This stage co mbines (i) selec tion of a subgroup J ⊂ { 1 , . . . , L } of SVD comp onents; (ii) hank elization (a v eraging along en- tries w ith indices i + j = const ) of the L × K matrix from th e selected J comp onen ts of the SVD; (iii) reconstruction of a time series comp onen t of length N from the Hank el ma tr ix b y the ment ioned one-to-one corr esp on- dence (like in (1) bu t in the reve r s e dir ection, s ee b elo w the exact formulae). 3 The resu lt of the r econstruction stage is a time s er ies additiv e co mp onent : X J = X j ∈J p λ j U j V j T → G = ( g 0 , . . . , g N − 1 ) . F or the sak e of brevit y , let us d escrib e the hank elization of the matrix X J and the s u bsequent reconstruction of a time series comp onen t G as b eing app lied to a matrix Y =  y ij  i = L,j = K i,j =1 as it is introdu ced in [11]. First w e in tro duce L ∗ = min { L, K } , K ∗ = max { L, K } and defin e an L ∗ × K ∗ matrix Y ∗ as giv en by Y ∗ = Y if L 6 K and Y ∗ = Y T if L > K . Th en the elemen ts of the time series G = ( g 0 , . . . , g N − 1 ) formed from the matrix Y are calculated b y a v eraging al ong cross-diagonals of matrix Y ∗ as g n =                  1 n +1 n +1 P m =1 y ∗ m,n − m +2 , 0 6 n < L ∗ − 1 , 1 L ∗ L ∗ P m =1 y ∗ m,n − m +2 , L ∗ − 1 6 n < K ∗ , 1 N − n N − K ∗ +1 P m = n − K ∗ +2 y ∗ m,n − m +2 , K ∗ 6 n < N . (2) Changing the wind o w length p arameter and, what is more imp ortan t, the subgrou p J of SVD comp onen ts used for reconstruction, one can c h ange the outp u t time series G . In the problem of trend extraction, we are lo oking for G approximat in g a tr en d of a time series. Thus, the t r end extractio n problem in SSA is reduced to (i) the c hoice of a windo w length L used for decomp osition and (ii) the selection of a su bgroup J of SVD comp onen ts used for reconstruction. The fi rst pr oblem is thoroughly discussed in section 1.6 of [1 1 ]. In this paper , w e p rop ose a solution for the second problem. Note that f or the reconstru ction of a time series comp onent, SSA consid- ers the wh ole time series, as its algorithm uses S VD of th e tra jectory m atrix built fr om all parts of the time series. Therefore, SS A is not a lo cal metho d in contrast to a lin ear fi ltering or wa v elet metho ds . On the other h an d , this prop erty makes SSA r obust to outliers, see [11] for more details. An essentia l disad v anta ge of SSA is its computational complexity for th e calculatio n of S VD. Th is shortcoming can b e reduced by using mo d ern [9] and parallel algorithms for SVD. Moreo v er, for trend revision in case of receiving n ew data p oin ts, a computationally attractiv e algorithm of [12 ] for up d ating SVD can b e used. It is worth to men tion here that the similar id eas of using SVD of the tra- jectory matrix ha ve b een prop osed in other areas, e.g. in signal extraction in o ceanology [8] and estimatio n of parameters of damp ed complex exp onen tial signals [13]. 4 2.1 T rend in SSA SSA is a non p arametric approac h whic h do es n ot need a pr iori sp ecification of mo dels of time series and trend, neither d eterministic nor sto c hastic ones. The classes of trends and residuals which can b e successfully s eparated by SSA are c haracterized as follo ws. First, since we extract an y trend b y select in g a sub group of all d SVD comp onen ts, this tr en d should generate less than d SVD comp onent s. F or an infi n ite time series, a class of su c h trend s coincides with the class of time series go verned b y fin ite difference equations [11]. This class can b e d escrib ed explicitly as linear co mbinations o f pro du cts of p olynomials, exp onen tials and sin es [6 ]. An elemen t of th is class suits w ell for repr esen tation of a smo oth and slo w v arying trend. Second, a residual should b elong to a cl ass of time series w hic h ca n b e separated from a trend. The separabilit y theory du e to [14] helps us deter- mine this class. In [14] it w as pro ved that (i) an y deterministic fun ction can b e asymptotic ally separated from an y ergo dic sto c hastic noise as the time series length and window length tend to infinit y; (ii) un der some conditions an y trend can b e separated from an y quasi-p erio dic comp onen t, see also [11]. These pr op erties of SSA m ak e this approac h feasible for trend extraction in the p resence of n oise and quasi-p erio dic oscillating components. Finally , as trend is a smo oth and s lo w v arying time series comp onen t, it generates SVD comp onent s with smo oth and slo w v arying EOFs. Eigen ve c- tors repr esent an orthonormal basis of a tra jectory v ector sp ace spann ed on the column s of tra jectory matrix. T h us e ach EO F is a linea r com bination of p ortions of the corresp onding time series and in h erits its global smo oth- ness prop erties. This idea is considered in detail in [11] for the cases of p olynomial and exp onen tial trends. 2.2 Existing metho ds of t rend extraction in SSA A naiv e app roac h to trend extrac tion in SSA is to reconstruct a tr end from sev eral fir st S VD co mp onent s. Despite its simplicit y , this approac h works in man y real-life cases for the follo wing reason. An eigen v alue repr esen ts a con tribution of the corresp onding SVD comp onent int o the form of the time series, see sectio n 1 .6 o f [11]. Since a trend usually c haracterizes the shap e of a time series, its eigenv alues are larger th an the other ones, that implies small order n umbers of the trend SVD comp onent s. Ho we ver, the selecti on pro cedure fails when the v alues of a trend are s m all enough as c omp ared with a r esidual, or w hen a trend has a complicated structure (e. g. a high- 5 order p olynomial) and is characte r ized by man y (not only by the first on es) SVD comp onen ts. A smarter w ay of selecting tr en d SVD comp onen ts is to c h o ose the com- p onents with smooth and slo w v arying EOFs (w e ha ve explained this f act ab o ve) . At presen t, there exist only one parametric method of [15] whic h follo ws this approac h. In [15] it w as prop osed using the Kendall co r rela- tion co efficien t for testing for m onotonic gro wth of an EOF. Un fortunately , this met h o d is far from per f ect since it is not p ossible to establish wh ich kinds of tr en d ca n b e extrac ted by it s means. This metho d seems to be aimed at extraction of monotonic tren d s b ecause their EOFs are usually monotonic. Ho we ver, ev en a monotonic trend can p ro duce non-monotonic EOF, esp ecially in case of noisy obs er v ations. An example could b e a linear trend which generate s a linear and a constan t EO Fs. If there is a n oise or another time series comp onen t added, then t h is comp onent is often mixed with trend comp onen ts corrup ting its EOFs. Then , ev en in case of v ery small corruption, the constant EOF can b e highly non monotonic. Nat- urally , the metho d using the Ken d all co rr elation coefficient does not suit for non monotonic tr ends p ro ducing n on monotonic EOFs. F or example, a p olynomial of lo w order w hic h is often used f or trend mod elling u sually pro du ces non monotonic EOFs, f or details see e.g. [11]. 3 PR OPOSED METHOD F OR TREND EXTRA C- TION In this s ection, we presen t our metho d of tren d extractio n . First, follo w - ing [11], w e introd uce the p erio dogram of a time series. Let u s consid er the F ourier representa tion of the elemen ts of a time s er ies X of lengt h N , X = ( x 0 , . . . , x N − 1 ), see e.g . section 7.3 of [7 ]: x n = c 0 + X 1 6 k 6 N − 1 2  c k cos(2 πn k / N ) + s k sin(2 π nk / N )  + ( − 1) n c N/ 2 , where k ∈ N , 0 6 n 6 N − 1, and c N/ 2 = 0 if N is an odd n um b er. T hen the p erio dogram of X at th e frequencies ω ∈ { k / N } ⌊ N/ 2 ⌋ k =0 is d efined as I N X ( k / N ) = N 2      2 c 2 0 , k = 0 , c 2 k + s 2 k , 0 < k < N / 2 , 2 c 2 N/ 2 , if N an ev en n umb er and k = N / 2 . (3) 6 Note that this per io d ogram is different from the perio d ogram usu ally used in sp ectral analysis, see e.g. [4] or [7 ]. T o sh o w this difference, let us denote the k -th elemen t o f the discrete F ourier transform of X as F k ( X ) = N − 1 X n =0 e − i 2 π nk / N x n , then th e p erio dogram I N X ( ω ) at the frequencies ω ∈ { k / N } ⌊ N/ 2 ⌋ k =0 is calculated as I N X ( k / N ) = 1 N ( 2 |F k ( X ) | 2 , if 0 < k < N / 2 , |F k ( X ) | 2 , if k = 0 or N is ev en and k = N / 2 . One ca n see that in addition to the n ormalizatio n differen t fr om that in [4] and [7] , the v alues for f requencies in the interv al (0 , 0 . 5) a r e m ultiplied b y t wo . This is done to ensu re the follo win g prop ert y: || X || 2 2 = N − 1 X n =0 x 2 n = ⌊ N/ 2 ⌋ X k =0 I N X ( k / N ) . (4) Let us introd uce th e cum u lativ e co ntribution of the f requencies [0 , ω ] as π N X ( ω ) = P k :0 ≤ k / N ≤ ω I N X ( k / N ), ω ∈ [0 , 0 . 5]. T hen, for a give n ω 0 ∈ (0 , 0 . 5), w e defi ne the cont r ibution of lo w frequ encies from the in terv al [0 , ω 0 ] to X ∈ R N as C ( X, ω 0 ) = π N X ( ω 0 ) /π N X (0 . 5) . (5) Then, giv en parameters ω 0 ∈ (0 , 0 . 5) and C 0 ∈ [0 , 1], w e prop ose to select those SVD comp onen ts w hose eigen v ectors sa tisfy the follo wing criterion: C ( U j , ω 0 ) > C 0 , (6) where U j is th e corresp onding j -th eigen v ector. One ma y in terp r et this metho d as sele ction of S VD comp onents with EOFs most ly c haracterized b y low-frequency fluctuations. It is w orth noting here that w h en w e apply C , π or I (defi ned ab o ve for a time series) to a vecto r , they are simply applied to a series of elemen ts of th e v ector. Ha ving the trend SVD comp onents selec ted using (6), one reconstructs the trend according to secti on 2. Th e question is ho w to s elect ω 0 and how to define the threshold C 0 . These issues are discussed in sections 4 an d 5, resp ectiv ely . 7 4 THE LOW-FREQUE NCY BOUND AR Y ω 0 The lo w-frequency b oundary ω 0 manages the scale of the extrac ted trend: the lo wer is ω 0 , the slo wer v aries th e extracted trend. Selectio n of ω 0 can b e done a priori based on ad d itional information ab out the data th us pr esp ec- ifying th e desired scale of the trend. F or example, if we assu me to ha v e a quasi-p erio dic co mp onent with kno wn p erio d T , then we should select ω 0 < 1 /T in order n ot to include this comp onen t in the trend. F or extraction of a trend of m on thly data with p ossible seasonal oscillati ons o f p erio d 12, we suggest to select ω 0 < 1 / 12 , e.g. ω 0 = 0 . 075. In this pap er w e also prop ose a metho d of selecti on of ω 0 considering a time series p erio dogram. Since a trend is a slow v arying comp onen t, its p erio dogram has large v alues close to zero frequency and small v alues for other frequencies. The problem of selec ting ω 0 is the problem of finding suc h a lo w-frequency v alue that the frequencies corresp onding to the large trend p er io d ogram v alues are insid e the interv al [0 , ω 0 ]. A t the same time, ω 0 cannot b e to o large b ecause then an oscillating comp onent with a fre- quency less than ω 0 can b e included in the trend p ro duced. Considering the p erio dogram of a tren d, w e could find the pr op er v alue of ω 0 but for a giv en time series its trend is u n kno wn . What w e pr op ose is to c ho ose ω 0 based on the p erio dogram of the original time series. The foll owing pr op osition subs tantiat es this app roac h. Prop osition 4.1. L et us have two time series G = ( g 0 , . . . , g N − 1 ) and H = ( h 0 , . . . , h N − 1 ) of length N , then for e ach k : 0 ≤ k ≤ ⌊ N / 2 ⌋ the fol lowing ine quality holds:   I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )   6 2 q I N G ( k / N ) I N H ( k / N ) . (7) Pr o of of Pr op osition 4.1. Let us first consider the case when 0 < k < N / 2. W e denote as c k ,X and s k ,X the co efficien ts of F ourier r epresen tation of a time series X used in the p erio d ogram definition (3). Then, b y this defin ition, I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N ) = = N 2  c 2 k ,G + H + s 2 k ,G + H − c 2 k ,G − s 2 k ,G − c 2 k ,H − s 2 k ,H  . Since c k ,G + H = 2 N ℜF k ( G + H ) = c k ,G + c k ,H (where ℜ z denotes an imaginary part of a complex n umb er z ) and, analog ous ly , s k ,G + H = s k ,G + s k ,H , w e ha ve I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N ) = N ( c k ,G c k ,H + s k ,H s k ,H ) . (8) 8 Let us consider the per io d ograms m ultiplication used in t h e r ight part of (7): I N G ( k / N ) I N H ( k / N ) = N 2 4  c 2 k ,G + s 2 k ,G   c 2 k ,H + s 2 k ,H  . (9) Since for all real a, b, c, and d it holds that ( a 2 + b 2 )( c 2 + d 2 ) = ( | ac | + | bd | ) 2 + ( | ad | − | bc | ) 2 , then I N G ( k / N ) I N H ( k / N ) = N 2 4 ( | c k ,G c k ,H | + | s k ,G s k ,H | ) 2 +( | c k ,G s k ,H | − | c k ,H s k ,G | ) 2 . (10) Finally , taking the s quare of (8), dividin g it b y four and taking int o ac- coun t (10), w e ha v e 1 4  I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )  2 = = N 2 4 ( c k ,G c k ,H + s k ,G s k ,H ) 2 6 N 2 4 ( | c k ,G c k ,H | + | s k ,G s k ,H | ) 2 6 6 N 2 4 ( | c k ,G c k ,H | + | s k ,G s k ,H | ) 2 + ( | c k ,G s k ,H | − | c k ,H s k ,G | ) 2 = = I N G ( k / N ) I N H ( k / N ) and the inequalit y in (7) holds 0 < k < N / 2. Second, w e consider the case when k = 0 o r k = N / 2. Again, b y the definition of the p erio dogram 2 q I N G ( k / N ) I N H ( k / N ) = 2 q N 2 c 2 k ,G c 2 k ,H = 2 N | c k ,G c k ,H | . A t the same time,   I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )   = N   c 2 k ,G + H − c 2 k ,G − c 2 k ,H   = N | 2 c k ,G c k ,H | whic h leads for k = 0 or k = N/ 2 to   I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )   = 2 q I N G ( k / N ) I N H ( k / N ) and the result in (7) holds with equalit y . Corollary 4 .2. L et us define f or a time series F of length N the fr e quency supp ort of the p erio do gr am I N F as a subset of fr e quencies { k / N } ⌊ N/ 2 ⌋ k =0 such that I N F ( k ′ / N ) > 0 for k ′ / N fr om this subset. If the f r e quency supp orts of two time series G and H ar e disjoint then I N G + H ( k / N ) = I N G ( k / N ) + I N H ( k / N ) 9 Let us demonstrate that when supp orts of p erio dograms o f time series G and H a r e nearly disjoint, th e p erio dogram of the su m G + H is clo s e to the su m of their p erio dograms. The fact that the p erio dograms of G and H are v ery d ifferen t at k/ N can b e expressed as I N G ( k / N ) /I N H ( k / N ) = d ≫ 1 , since withou t loss of generalit y w e can assume I N G ( k / N ) > I N H ( k / N ). Then using Prop osition 4.1 we ha ve that   I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )   6 6 2 q I N G ( k / N ) I N H ( k / N ) = 2 √ d I N G ( k / N ) ≪ I N G ( k / N ) , that means that the difference   I N G + H ( k / N ) − I N G ( k / N ) − I N H ( k / N )   is sig- nifican tly smaller than the v alue of th e largest p erio dogram (of I N G , I N H ) at the p oin t k / N . In man y applications, the giv en time series ca n b e modelled as made of a trend with large p erio dogram v alues at lo w-frequency in terv al [0 , ω 0 ], oscillati ons with p erio d s smaller than 1 /ω 0 , and noise whose frequ ency con- tribution sp reads o ve r all the fr equ encies [0 , 0 . 5] bu t is relativ ely small. In this case the p erio dogram supp orts of the trend and the residual can b e considered as nearly disj oint. Therefore, from Corollary 4.2, w e conclud e that the p erio d ogram of the time series is approxi m ately equal to the sum of the per io d ograms of th e trend, oscillations and noise. F or a time series X of length N , w e prop ose to select the v alue of the parameter ω 0 according to the follo wing rule: ω 0 = max k / N , 0 6 k 6 N/ 2  k / N : I N X (0) , . . . , I N X ( k / N ) < M N X  , (11) where M N X is the median of the v alues of p erio d ogram of X . The mo d - elling of a time series as a sum of a trend , oscillati ons and a noise (let us supp ose to ha v e a normal noise) motiv ates this rule as follo ws. Sin ce the frequency su pp orts of the trend and oscillat in g co mp onent s do not o ve r lap, only the v alues of the noise p erio dogram can mix with the v alues of the trend p erio dogram. First, the v alues of the noise p erio dogram for neighbor- ing ordinates are asymptotic ally indep endent (see e.g. sec tion 7. 3.2 of [7]). Second, supp osing a relat ively long time series an d narrow frequency sup- p orts of trend and oscillating comp onen ts, the median of v alues of the time series p er io d ogram gi ves an estimation of the median of the v alues of the 10 noise p erio dogram. S in ce a trend is sup p osed to ha v e large con tribution to the shap e of the time se r ies (i.e. a l arge L 2 -norm) compared to the noise and its frequency supp ort i s quite narro w compared to the whole in terv al [0 , 0 . 5], its p erio dogram v alues are r elativ ely larger than th e median of the noise p erio dogram v alues due to (4). Therefore, the condition used in (11 ) is fulfilled only for such a frequ ency ω 0 that the trend p erio dogram v alues is close to zero (outside the trend fr equency int erv al). Large n oise p erio d ogram v alues in this frequency region can lead to selecting larger than necessary ω 0 . But remem b er th at we co mp are the p er io d ogram v alues with their m edian and the noise p erio dogram v alues are indep end en t (asymptotically). Hence, with p r obabilit y app ro ximately equal to 1 − 0 . 5 m (e.g. this v alue is equal to 0.9375 for m = 4) we select the m -th p oint (of the grid { k / N } ) lo cated to the right side of the tr en d frequency interv al (where the trend p eridogram v alues are larger then the noised p erio dogram median). Note that the lengths N of the time series and L of e igenv ector are differen t ( L < N ) whic h causes differen t resolution of their p erio dograms. Ha ving estimated ω 0 after co n sideration of the p er io d ogram of the original time series, one should select ω ′ 0 = ⌈ Lω 0 ⌉ /L. (12) Dep endence of ω 0 on the time series resolution. Let us defin e the resolution ρ of the original time series as ρ = ( τ n +1 − τ n ) − 1 , where τ n is the time of n -th measur emen t. If one ha v e estimate d ω 0 for the data with resolution ρ and there co mes th e same data but measured with higher res- olution ρ ′ = mρ ( m ∈ R ) thus increasing the data length in m times, then in order to extract the same trend , one sh ould tak e the new thr eshold v alue ω ′ 0 = ω 0 /m . In a similar manner, after decima tion of the data reducing the resolution in m times, th e v alue ω ′ 0 = mω 0 should b e tak en. Example 4.3 (Th e c hoice of ω 0 for a noised exp onent ial trend) . L et us c onsider an examp le of sele ction of the tr eshold ω 0 for an exp onential tr end and a white Gaussian noise which also demo nstr ates Pr op osition 4.1. L et the time series F = G + H b e of length N = 1 20 , wher e the c omp onents G and H ar e define d as g n = Ae 0 . 01 n , h n = B ε n , ε n ∼ iid N (0 , 1) and A, B ar e sele cte d so that || G || 2 = || H || 2 = P N − 1 n =0 g n = P N − 1 n =0 h n = 1 . The normaliza tion is done to ensur e that P 60 k =0 I N G ( k / N ) = P 60 k =0 I N H ( k / N ) = 1 . Figur e 1 sho ws a) the simulate d time series F , b) its c omp onents, c) th e p erio do gr ams of the c omp onents, d) the p erio do gr ams zo ome d to gether with a line c orr esp onding to the me dian of the noise p erio do gr am values e qual 11 to 0 . 0126, e) the p erio do gr am I N F of F and a kind o f “c onfidenc e” interval of its estimat ion I N G + I N H c alculate d ac c or ding Pr op osition 4.1 and a line c orr esp onding to the me dian M 120 F of t he time series p erio do gr am v alues (use d for estimating ω 0 ), and f ) the discr ep ancy, the differ enc e b etwe en I N F and I N G + I N H to gether with the values of this differ e nc e estimate d in the right side of (7). Note tha the me dian o f the p erio do gr am values of F is e qual to 0.0141 , which is clo se to the me dian o f the noise p erio do gr am val u es e qual to 0.0126. The value of ω 0 estimate d ac c or ding to the pr op ose d rule (11) is e qual to 6 / 12 0 = 0 . 05 . a) Origi nal time series F = G + H 0 20 40 60 80 100 120 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 n F = G + H b) Time series comp onents G and H 0 20 40 60 80 100 120 −0.2 −0.1 0 0.1 0.2 0.3 n G H c) Periodog rams of comp onents G and H 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 ω I 120 G ( ω ) I 120 H ( ω ) d) Per i o dograms of G and H (zoomed) 0 0.1 0.2 0.3 0.4 0.5 0 0.02 0.04 0.06 0.08 0.1 ω I 120 G ( ω ) I 120 H ( ω ) this line corresponds to the median of the noise periodogram values, equal to 0.0126 e) Periodogram o f F and its estimates 0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 ω I G + H ( ω ) I G ( ω ) + I H ( ω ) + C ( ω ) I G ( ω ) + I H ( ω ) − C ( ω ) this line corresponds to M 120 F , the median of the time series periodogram values, equal to 0.0141 estimated value of ω 0 f ) Discrepancy 0 0.1 0.2 0.3 0.4 0.5 −0.1 −0.05 0 0.05 0.1 ω I G + H ( ω ) − I G ( ω ) − I H ( ω ) + C ( ω ) − C ( ω ) Figure 1: T he c hoice of ω 0 for an exp onential trend and Gaussian noise; The v alue C ( ω ) used in the legends is equal to 2 q I N G ( ω ) I N H ( ω ). 5 THE LO W-F REQUENCY CONTRIBUTION C 0 Before suggesting a pro cedure for selection of the second paramete r of the prop osed metho d, the lo w-frequency threshold C 0 , w e inv estiga te the effect of the c hoice of C 0 on the qualit y of the trend extracted. F or th is aim, we con- 12 sider a time series mo del with a trend that generates SVD comp onen ts with kno wn n u m b ers. Then , for a s u fficien t n umb er of sim ulated time series, we compare our trend extractio n procedu re with a SSA-based pr o cedure which simply reconstru cts the trend using the kno wn trend SVD comp onents. 5.1 A sim ulation example: an exp onen tial trend plus a Gaus- sian noise The mo d el considered is the s ame as in examp le ab o ve. Let the time series F = ( f 0 , . . . , f N − 1 ) consist of an exp onen tial trend t n plus a Gaussian white noise r n : f n = t n + r n , t n = e αn , r n = σ e αn ε n , ε n ∼ iidN (0 , 1) . (13) According to [11], for suc h a time series with mo derate noise the first SVD comp onen t corresp onds to the trend. W e c onsid ered only the noise lev- els wh en this is tru e (empirically c hec k ed). Note that the noise r n has a m ultiplicativ e m o del as its standard deviation is prop ortional to the trend . In the follo win g, we consider the follo wing prop er ties. Fi r st, w e calculate the difference b et w een the trend ˆ t n ( C 0 ) resulted from our metho d with C 0 used and the reconstru ction ˜ t n of the first S VD c omp onen t exploiting the w eight ed mean square error (MSE ) b ecause this measure is more rele v ant for a model with a m ultiplicativ e noise th an a sim p le MSE: D ( ˆ t n ( C 0 ) , ˜ t n ) = 1 N N − 1 X n =0 e − 2 αn  ˆ t n ( C 0 ) − ˜ t n  2 . (14) This m easur e compares our trend and the ideal SS A trend. Second, w e calculate the weigh ted mean square errors b etw een ˆ t n ( C 0 ), ˜ t n and the tru e trend t n separately: D ( ˆ t n ( C 0 )) = 1 N N − 1 X n =0 e − 2 αn  t n − ˆ t n ( C 0 )  2 , D ( ˜ t n ) = 1 N N − 1 X n =0 e − 2 αn  t n − ˜ t n  2 . (15) 5.1.1 Sc heme of estimation of the errors using simulation The errors (14), (15) are estimated using the follo wing scheme. W e sim ulate S r ealizat ions of the time series F acco r d ing to the mo d el (13) and calculate the mean of D ( ˆ t n ( C 0 ) , ˜ t n ) for all v alues of C 0 from the large grid 0:0 . 01:1: D ( ˆ t n ( C 0 ) , ˜ t n ) = 1 S S X s =1 D ( ˆ t ( s ) n ( C 0 ) , ˜ t ( s ) n ) , (16) 13 where ˆ t ( s ) n ( C 0 ) and ˜ t ( s ) n denote trends of th e s -th sim ulated time series. T h e mean errors D ( ˆ t n ( C 0 )), D ( ˜ t n ) b etw een the true trend t n and the extrac ted trends ˆ t n ( C 0 ) and ˜ t n , resp ectiv ely , are calculated similarly . Let us also denote the min imal v alues of the mean errors as D min ( ˆ t n , ˜ t n ) = m in C 0 D ( ˆ t n ( C 0 ) , ˜ t n ) , D min ( ˆ t n ) = min C 0 D ( ˆ t n ( C 0 )) (17) and the v alue of C 0 pro vidin g the minimal mean error b et w een the extracted trend and the ideal S SA trend as C opt 0 = arg min C 0 D ( ˆ t n ( C 0 ) , ˜ t n ) , so that D min ( ˆ t n , ˜ t n ) = D ( ˆ t n ( C opt 0 ) , ˜ t n ). The simulate d time series are of length N = 47. In order to ac hiev e the b est separabilit y [11] we ha ve selected the SSA window length L = ⌈ N / 2 ⌉ = 24. The estimates of the mean errors are calculated on S = 10 4 realizatio n s of the time s eries. W e consider different v alues of the mo del parameters α and σ . The v alues of α are 0 (co r resp ondin g to a constan t trend), 0.01 and 0.02 whic h corresp ond to the increase of trend v alues (from t 0 to t N − 1 ) in 1, 1.6 and 2.5 times, resp ectiv ely . The lev els of noise are 0 . 2 6 σ 6 1 . 6. It was empirically c hec ked that for su c h le vels of noise the fi rst SVD co mp onent corresp onds to the trend. Moreo v er, w e estimated the probabilit y of the t yp e I error of not selecting the first SVD comp onen t as the ratio of times when the fi rst comp on ent is not iden tified as a trend comp onen t by our p r o cedure to the num b er of rep etitions S . Choice of ω 0 . In order to select the lo w-frequ en cy threshold ω 0 , w e con- sidered sev eral sim ulated time series with differen t α and the maximal noise σ = 1 . 6. Tw o examples of their p erio dograms for α = 0 and α = 0 . 02 are dep icted in Figure 2. The median v alues for th e p erio dograms depicted in Figure 2 are 2.936 a n d 2.924 whic h leads to ω 0 = 0 for α = 0 and ω 0 = ⌈ 1 / N · L ⌉ /L = 1 / 24 ≅ 0 . 042 for α = 0 . 0 2 estimated using (12). W e decided to take the same ω 0 = 0 . 042 (the largest one) for all α co n sidered. 5.1.2 Sim ulation results Figure 3 sho ws the evo lu tion of the square ro ots of th e mean errors and C opt 0 as a fun ction of σ . The v alues α = 0 and α = 0 . 02 are u sed. The square ro ots 14 0 0.1 0.2 0.3 0.4 0.5 0 50 100 ω Periodograms of the simulated time series I F ( ω ) , α = 0 I F ( ω ) , α = 0 . 0 2 this line corresponds to the medians of the time series periodogram values, equal to 2.936 for α =0 and 2.924 for α =0.02 Figure 2: The p erio dograms of t wo time ser ies of the m o del (13) with σ = 1 . 6 and α = 0 , 0 . 02. of the mean errors (i.e. sta n dard deviations) are tak en for b etter comparison with σ whic h is the standard deviation m ultiplier of the noise. The plots of the m in imal mean error D min ( ˆ t n , ˜ t n ) and the optimal C opt 0 for α = 0 . 02 are depicted in F igur e 3, where the v alues for α = 0 are also sho wn in gra y colo r . The estimates for α = 0 . 01 are not rep orted h ere. The in terpr etation of the pro duced results is as follo w s . First, the trend extracted with the optimal C 0 is v ery similar to the id eal S SA trend , r econ- structed b y the first SVD comp onen t since D min ( ˆ t n , ˜ t n ) ≪ D min ( ˆ t n ) (the error b etw een our trend and the id eal trend is muc h s m aller th an the error of the ideal trend itself ), especially when σ 6 0 . 8. Moreo v er, the estimated probabilit y of th e t yp e I error (i.e . th e pr obabilit y of not selec ting the first SVD comp onen t) is less than 0.05 for σ 6 1 . 4. All this allo ws us to conclude that in case of an exp onen tial trend and a wh ite Gaussian noise the prop osed metho d of tr en d extraction with an optimal C 0 with high probabilit y selects the r equ ired first S VD comp onent corresp onding to the tr end. The trend ˆ t n ( C opt 0 ) extracted with an optimal C 0 estimates the tru e trend quite go o d when comparing the deviation q D min ( ˆ t n ) with the noise stan- dard deviation σ . F or example, for σ = 1 . 6 th e v alue of q D min ( ˆ t n ) is appro ximately equal to 0 . 5. Note that for d ifferent α the mean errors D min ( ˆ t n ) are very similar though the used optimal v alues of C 0 are qu ite differen t (F igur e 3). This sho ws that the metho d adapts to the c han ge of the mo d el parameter α . Let us consider the dep enden ce of inaccuracy of the prop osed trend ex- 15 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 σ Squ a r e r o o t s o f t he me a n D -e rr or s D min ( ˆ t n , ˜ t n ) , α = 0 D min ( ˆ t n , ˜ t n ) , α = 0 . 02 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0.8 0.85 0.9 0.95 1 σ C 0 providing m inimal D ( ˆ t n ( C 0 ) , ˜ t n ) C opt 0 , α = 0 C opt 0 , α = 0 . 02 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0.4 0.5 σ Squ a r e r o o t s o f t he me a n D -e rr or s D min ( ˆ t n ), α = 0 D min ( ˆ t n ), α = 0 . 0 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0.4 0.5 σ Square ro o t s o f t he me a n D -e rr or s D ( ˜ t n ), α = 0 D ( ˜ t n ), α = 0 . 0 2 Figure 3: Th e square ro ots of the mean errors D min ( ˆ t n , ˜ t n ) (top left) D min ( ˆ t n ) (b ottom left) and D ( ˜ t n ) (b ottom righ t) as well as the optimal C 0 v alue providing a minimal mean err or D min ( ˆ t n , ˜ t n ) b et we en the extracted trend and the ideal S SA trend (top righ t); all for α = 0 and α = 0 . 02. traction metho d on the v alue of C 0 . As ab o v e, the inaccuracy is measured with the minimal mean err or D min ( ˆ t n , ˜ t n ) b etw een the extract ed trend and the ideal SS A trend. Figure 4 sho ws the graphs of this error as a fu nction of C 0 for d ifferen t exp onen tials α and n oise lev els σ . One can see that it is cru cial not to select to o large C 0 since in this case the trend comp onen t can b e not in cluded in the reconstruction (that is also confirmed by the esti m ated probabilit y of the type I error whic h is not rep orted here). A t the same time without significan t loss of acc u racy one can choose C 0 smaller than C opt 0 (corresp onding to the b est accuracy). T his is true du e to th e small con tribution of eac h of n oise comp on ents wh ich can b e erron eously included for C 0 < C opt 0 . 16 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 C 0 D ( ˆ t n ( C 0 ) , ˜ t n ) α = 0 , σ = 0 . 8 α = 0 , σ = 1 . 4 α = 0 . 02 , σ = 0 . 8 α = 0 . 02 , σ = 1 . 4 Figure 4: The error D ( ˆ t n ( C 0 ) , ˜ t n ) as a fu nction of C 0 for differen t com bin a- tions of α = 0 , 0 . 02 and σ = 0 . 8 , 1 . 4. 5.2 Heuristic pro cedure for t he c hoice of C 0 Based on the obser v ations of section 5.1, w e prop ose the follo win g heur istic pro cedure for c ho osing th e v alue of the metho d lo w-fr equ ency threshold C 0 . As discuss ed, trend EOFs v ary slo w. First we sho w that this prop er ty is inherited by the trend elemen tary reconstructed comp onents, the time series comp onen ts eac h r econstructed fr om one tren d SVD comp onen t. Prop osition 5.1. L et ( √ λ, U, V ) b e an eigentriple o f SSA de c omp osition of a time series F , U = ( u 1 , . . . , u L ) T , V = ( v 1 , . . . , v L ) T , and G b e a time series r e c onstructe d by this eigentriple. If it is true that ∃ δ 1 , δ 2 ∈ R : ∀ k , 1 6 k 6 L − 1 : | u k +1 − u k | < δ 1 , | v k +1 − v k | < δ 2 , then for the elements of G = ( g 0 , . . . , g N − 1 ) the fol lowing holds: ∃ ǫ ( δ 1 , δ 2 ) : ∀ n, L ∗ − 1 6 n < K ∗ :    g n +1 − g n    < ǫ ( δ 1 , δ 2 ) , wher e L ∗ = min { L, K } , K ∗ = max { L, K } . Pr o of of Pr op osition 5.1. O ne can easil y pro v e this p r op osition taking in to accoun t ho w the elemen tary r econstructed comp onent G is constructed from its eigent r ip le ( √ λ, U, V ), see section 2. First, the matrix Y = √ λU V T is constructed. Second, th e hanke lization of Y is p erformed. 17 Let us sho w how to ca lculate ǫ u sing (2) for δ 1 , δ 2 when L 6 K . F or other cases ǫ ( δ 1 , δ 2 ) is ca lculated similarly .    g n +1 − g n    = √ λ L     L X m =1  u m v n − m +3 − u m v n − m +2      < < √ λ L L X m =1 | u m || v n − m +3 − v n − m +2 | < < √ λ L δ 2 L X m =1 | u m | < δ 2 √ λ L  u 1 + ( L − 1) δ 1  . Let us ha ve a time s er ies F and denote its trend extracted w ith the metho d with parameters ω 0 , C 0 as T ( ω 0 , C 0 ). In order to prop ose the pro cedure selecting C 0 , we first define the norm alized con tribution of lo w- frequency oscillations in the resid ual F − T ( ω 0 , C 0 ) as: R F ,ω 0 ( C 0 ) = C  F − T ( ω 0 , C 0 ) , ω 0  C ( F , ω 0 ) − 1 , where C is defin ed in equation (5). Based on Prop osition 5.1 , w e exp ect that the ele mentary reconstructed comp onen ts corresp onding to a tren d hav e large cont r ibution of low fre- quencies. Thus, the maximal v alues of C 0 whic h lead to selection of tr en d- corresp onding SVD comp onen ts should generate ju mps of R F ,ω 0 ( C 0 ). Exploiting this id ea, we p rop ose the f ollo w ing w a y of c ho osing C 0 : C R 0 = min  C 0 ∈ [0 , 1] : R F ,ω 0 ( C 0 + ∆ C ) − R F ,ω 0 ( C 0 ) ≥ ∆ R  , (18) where ∆ C is a searc h step and ∆ R is the giv en threshold. On one h and, this strategy is heur istic and requires selection of ∆ R , b ut on the other h and, the sim ulation results and app lication to d ifferen t time series sho we d its abilit y to choose reasonable C 0 in man y cases. Based on this empirical exp erience, w e suggest us in g 0 . 05 ≤ ∆ R ≤ 0 . 1. The step ∆ C is to b e chosen as small as p ossible to discrimin ate id entificatio n s o ccurring at differen t v alues of C 0 . T o redu ce computational time, we commonly take ∆ C ≥ 0 . 01 and suggest a default v alue of ∆ C = 0 . 01. 6 EXAMPLES Sim ulated example with p olynomial trend. The firs t example illus- trates the c h oice of parameters ω 0 and C 0 . W e simulated a time series of 18 0 50 100 150 200 250 300 − 20 0 20 Original time series 0 50 100 150 200 250 300 − 20 0 20 Original time series Original trend Extracted trend 0 0.05 0.1 0.15 0.2 0.25 0 2000 4000 6000 frequency Time series periodogram 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 C 0 threshold ∆ R =0.05 R F,0.02 ( C 0 + ∆ C )− R F,0.02 ( C 0 ), ∆ C =0.01 C 0 R =0.53 obtained Figure 5: Simulated example with a p olynomial trend: original time series (top left); th e orig in al trend and an extracted one with L = 180, ∆ C = 0 . 01, and ∆ R = 0 . 05 (top righ t); zo omed ti m e series p er io d ogram inside ω ∈ [0 , 0 . 25] (b ottom left); th e v alues o f R F ,ω 0 ( C 0 + ∆ C ) − R F ,ω 0 ( C 0 ) u sed for the c hoice of C 0 resulted in a v alue C R 0 = 0 . 53 (bottom righ t). length N = 300 , sho wn in Figure 5, c ontaining a p olynomial trend, an exp onent ially-mo du lated sine w a ve, and a wh ite Gaussian noise, whose n -th elemen t is exp ressed as f n = 10 − 11 ( n − 10)( n − 70)( n − 160) 2 ( n − 290) 2 + exp(0 . 01 n ) sin(2 π n/ 12) + ε n , ε n is iidN (0 , 5 2 ). The p erio d of the sine wa v e is assu med to b e unknown. W e ha v e c hosen the wind o w length L = N/ 2 = 150 for ac hieving b etter separabilit y of trend and residu al. The v alue ω 0 = 6 / N = 0 . 02 w as selec ted using (11), w here the calculated median v alue is M N X ≅ 37 . 06. The searc h for C 0 using (18) has b een done with ste p ∆ C = 0 . 01 and ∆ R = 0 . 05. As sho wn in Figure 5, d espite of the strong noise and oscillatio ns , the extracte d trend appro ximates the original one v ery w ell. The ac h iev ed m ean squ are error is 0 . 79. F or example, the ideal lo w p ass fi lter with the cutoff frequency 0 . 02 pro duced the error of 3 . 14. This su p eriorit y is ac h iev ed mostly due to b etter appr o ximation at the first and last 50 p oin ts of the time series. All the calculations w ere p er f ormed using our Matlab-based soft wa re AutoSSA a v ailable at http://www.p dmi.ras.ru/ ∼ theo/autossa. 19 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 0 0.02 0.04 0.06 0.08 0.1 0 0.5 1 1.5 2 x 10 9 frequency Time series periodogram Figure 6: Unemp lo ymen t lev el in Alask a: original d ata (left-hand side panel), zo omed p erio dogram (right -h an d side pan el). T rends of the unemplo yment lev el. Let u s demonstrate extraction of tren d s of differen t scale. W e to ok the data of the un emplo yment lev el (unemplo yed p ersons) in Ala s k a for the p er io d 1976/0 1-2006/09 (mon thly data, seasonally adjus ted), provided b y the Bureau of Lab or Statistics at h ttp://www.bls.go v under the identifier LASST02000 004 (Figure 6). This time series is t ypical for economical app lications, where data con tain rela- tiv ely li ttle noise and are sub ject to abrupt c hanges. Eco n omists are often in terested in the “short” term tr end wh ic h includes cyclical fluctuations and is r eferred to as trend-cycle. The length of the data is N = 369. F or ac hieving b etter s ep arabilit y of trend and residual w e selected L close to N/ 2 bu t divisible b y the p erio d T = 12 of p robable seasonal oscil lations: L = 12 ⌊ N/ 24 ⌋ = 180. W e extracted trends of differen t scale s using the follo wing v alues of ω 0 : 0.01, 0 .02, 0.05, 0.0 75, and 0.095, see Figure 7 for t h e results. The v alue 0 . 095 ≅ ⌈ 33 / 369 · 180 ⌉ / 1 80 was s elected according to (12), where M N X ≅ 5 . 19 · 10 5 . The v alue 0.075 is th e default v alue for mon thly data (section 4 ). Other v alues (0.0 1, 0.02 a n d 0.05) w ere considered for b etter illustration of ho w the v alue of ω 0 influences th e scale of the extracted trend. Th e search for C 0 w as p erformed as describ ed in section 5 in the in terv al [0 . 5 , 1] with the step ∆ C = 0 . 01 and ∆ R = 0 . 05. 20 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 Original time series Trend with ω 0 =0.01 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 Original time series Trend with ω 0 =0.02 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 Original time series Trend with ω 0 =0.05 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 Original time series Trend with ω 0 =0.075 Jan80 Jan85 Jan90 Jan95 Jan00 Jan05 1 1.5 2 2.5 3 x 10 4 Original time series Trend with ω 0 =0.095 Figure 7: Un emplo yment lev el in Alask a: extracted trends of different scales with ω 0 = 0 . 01 , 0 . 02 , 0 . 05 , 0 . 075, and 0 . 095 ( L = 180, ∆ C = 0 . 01, and ∆ R = 0 . 05). 7 CONCLUSIONS SSA is an attractiv e approac h to trend extraction b ecause it: (i) requires no mo del sp ecification of time series and trend, (ii) extracts trend of n oisy time series con taining oscillatio n s of unknown p erio d . In th is pap er, we presen ted a metho d whic h inherits these prop erties and is easy to use since it requires selection of only t w o parameters. A cknow le dgments. Th e author warmly thanks his Ph.D. thesis advisor Nina Gol yandina for her sup ervision and guidance that help ed h im to pro- duce the r esults p resen ted in this pap er. The author greatly appreciate s an anon ymous review er for his v aluable and constructiv e comments. 21 References [1] Ale xandro v, T. (2 006). Softwar e p ackage for auto matic extr action and for e c ast of additive c omp onents of time series in the fr amework of the Caterpil lar-SSA appr o ach . PhD thesis, St.P etersbur g State Un iv er- sit y . In Russian, a v ailable at h ttp://www.p dmi.ras.ru / ∼ theo/autossa. [2] Ale xandro v, T., Biancon cini, S ., D a gum, E. B., Maass, P., and McElro y, T. S . (2008 ). A review of some mod ern approac hes to the problem of trend extraction. US C ensus Bureau T ec hRep ort RRS 2008 / 0 3. [3] Ale xandro v, T., a n d Gol y an dina, N. (2005) . Automatic extrac- tion and forecast of time s er ies cyclic comp onents within the f ramew ork of SSA. In “Pro c. of th e 5t h St.Pete r sburg W orkshop on Sim ulation”, 45–50 . [4] Brill inger , D. R. (2001) . Time Series: Data A nalysis and The- ory . So ciet y for I n dustrial and App lied M athematics, Philadelphia, P A, USA. [5] Broomhead, D. S. and King, G. P. (1986). Extracting qualitativ e dynamics fr om exp er im ental data. Physic a D 20 , 217–236 . [6] Buc h st a ber, V. M. (1995 ). Time series analysis and grassmannians. In “Amer. Math. So c. T rans”, v ol. 162. AMS, 1–17. [7] Cha tfield, C. (2003). The Analy sis of Time Series: A n Intr o duction . 6th ed., Chapman & Hall/CR C. [8] Cole b rook, J. M. (1978). Con tin uous plankton records — zo oplank- ton and e vir on m en t, northeast A tlant ic and North Sea, 1948–197 5. Oc e anol. A cta. 1 , 9–23. [9] Drma ˇ c, Z., and Veseli ´ c, K. (2 005). New fast and accurate jacobi svd algorithm: I, I I. T ec h . Rep. LAP A CK W orking Note 169, Dep. of Mathematics, Un iv ersit y of Zagreb, Croatia. [10] G hil, M.; Allen , R. M.; Det t inger, M. D.; Ide, K.; K on - drashov, D.; M a nn, M. E.; R ober tson, A.; S aunders, A.; Tian, Y.; V aradi, F.; and Yiou, P. (2002) . Adv anced sp ectral metho ds for climatic time series. R ev. Ge ophys. 40 , 1 , 1– 41. 22 [11] G ol y a ndina, N. E.; Nekrutkin, V. V., and Zhiglj a vs ky, A. A. (2001 ). Analysis of Time Series Structure: SS A and Related T ec h- niques. Boca Raton, FL: Chapman&Hall/CR C. [12] G u, M., and Eisenst a t, S. C. (1993 ). A stable and fast al- gorithm for u p dating the singular v alue d ecomp osition. T ec h. Rep. Y ALEU/DCS/RR-96 6, Dep. of Computer Science, Y al e Univ ers it y . [13] Kumaresan, R., and Tufts, D. W. ( 1980). Data-adaptiv e p rin- cipal co mp onent signal pr o cessing. In “Proc. of IEEE Conference On Decision and Con trol”, Albuquerque, 949–9 54. [14] Nekr utkin, V. (1996). Theoretica l proper ties of the “Caterpillar” metho d of time serie s analysis. In “Pro c. 8t h IEEE Signal Pro cessing W orkshop on Statistical S ignal and Ar r a y Pro cessing”, IEEE Comp u ter So ciet y, 395– 397. [15] V a ut ard, R.; Yiou, P., and Ghil, M. (1992). Singular-sp ectrum analysis: A toolkit for short, noisy c haotic signals. Phy sic a D 5 8 , 95– 126. 23