ÚVOD DO ČASOVÝCH ŘAD1 VÍTĚZSLAV VESELÝ Katedra aplikované matematiky a informatiky Ekonomicko-správní fakulta MU v Brně Email: vesely@econ.muni.cz OBSAH 1 Úvod 2 Základní pojmy 2.1 Definice časové řady 2.2 Systém distribučních funkcí 2.3 Momentové charakteristiky (střední hodnota a autokovarianční funkce) 2.4 Striktní a slabá stacionarita 2.5 Lineární transformace časových řad 2.6 Odhady momentových funkcí 3. Predikce v časových řadách 3.1 Princip ortogonální projekce 3.2 Durbin-Levinsonův a inovační algoritmus. 4 Modelování a odhad parametrů 4.1 Předzpracování * Ošetření specifických efektů * Stabilizace rozptylu * Identifikace periodických komponent 4.2 Dekompoziční model * Užití lineárního regresního modelu Date: 24.listopad 2003, korekce 9.2.2009. 1sborník semináře ANALÝZA DAT, listopad 2003 (TriloByte s r.o.) ˇ Filtrační techniky 4.3 Box-Jenkinsovy modely * ARMA modely pro stacionární časové řady * ARIMA modely pro kovariančně stacionární časové řady s trendem * SARIMA modely pro kovariančně stacionární časové řady s trendem a sezónní komponentou. 5 Demonstrační ukázky pomocí knihovny TSA-M 2 1. Úvod Časová řada: = soubor pozorování náhodných veličin {xt, t T}, kde t je zpravidla čas a T R. 1.1. UKÁZKY ČASOVÝCH ŘAD. Obr. 1.1: Měření proudu procházejícího odporem r v obvodu se střídavým napětím, v(t) = a cos(vt + ), tj. xt = a r cos(vt + ), kde a a jsou zatíženy náhodnou chybou. Obr. 1.2: Růst populace v USA v létech 1790­1980 sledovaný v desetiletých intervalech. Obr. 1.3: Počet stávek v USA v létech 1951­1980. Obr. 1.4: Měsíčný počty tisíců cestujících v mezinárodní letecké dopravě v létech 1949­1960. Obr. 1.5: Počty slunečních skvrn v létech 1770­1869. Obr. 1.6: Počet úmrtí při nehodách v USA v létech 1973­1978. 1.2. OBLASTI UPLATNĚNÍ METOD ANALÝZY ČASOVÝCH ŘAD. fyzika, technika: * seismický záznam v geofyzice. * řada nejvyšších denních teplot v meteorologii. * průběh výstupního signálu určitého elektrického přístroje. * tenzometrické měření povrchového napětí v provozu namáhané strojní součástky. biologie, ekologie: sledování různých parametrů znečištění ovzduší. 3 medicína: záznam EKG nebo EEG. společenské vědy: změny v počtu a složení obyvatelstva. sociologie: vývoj rozvodovosti. ekonomie: teorie časových řad = jedna z nejdůležitějších kvantitativních metod pro analýzu ekonomických dat, např.: * analýza poptávky po určitém výrobku * analýza objemu zemědělské produkce * analýza počtu cestujících v letecké dopravě * analýza vývoje kurzu akcií na burze 0 10 20 30 40 50 60 70 80 90 100 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Current through a resistor BD Example 1.1.1 Obr. 1.1 0 2 4 6 8 10 12 14 16 18 20 0 0.5 1 1.5 2 2.5 x 10 8 BD Example 1.1.2. Population in USA at ten-year intervals, 1790 - 1980 Obr. 1.2 0 5 10 15 20 25 30 3000 3500 4000 4500 5000 5500 6000 6500 BD Example 1.1.3. Strikes in the USA, 1951 - 1980 Obr. 1.3 0 50 100 150 100 200 300 400 500 600 700 BD Example 1.1.6. Airline passenger monthly totals (in thousands), Jan.49-Dec.60 Obr. 1.4 4 0 10 20 30 40 50 60 70 80 90 100 0 20 40 60 80 100 120 140 160 BD Example 1.1.5. The Wolfer sunspot numbers, 1770 - 1869 Obr. 1.5 0 10 20 30 40 50 60 70 80 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 BD Example 1.1.6. Monthly accidental deaths in the USA, 1973 - 1978 Obr. 1.6 1.3. CÍL ANALÝZY. Porozumění mechanismu, jímž se generují sledované údaje. Znalost modelu tohoto mechanismu znalost algoritmu, jímž můžeme chování tohoto mechanismu simulovat na počítači schopnost popsat s jistou přesností jeho chování: * mezi časovými okamžiky měření (interpolace) * v budoucnosti (extrapolace, prognóza) * s cílem řídit a optimalizovat činnost určitého systému vhodnou volbou vstupních a počátečních podmínek (regulace), např. regulace složitých technologických procesů. 5 2. Základní pojmy 2.1. DEFINICE ČASOVÉ ŘADY. Definice 2.1. Náhodný (stochastický) proces X je neprázdný systém (T = ) náhodných veličin definovaných na stejném pravděpodobnostním prostoru (, A, P). Píšeme X := {Xt | t T} nebo stručně {Xt}. Speciální případy: T R . . . proces se spojitým časem nebo náhodná funkce. T Z . . . proces s diskrétním časem, náhodná posloupnost nebo časová řada. Poznámka 2.2. Indexová množina T je obvykle uspořádaná a interpretuje se jako spojitý nebo diskrétní časový interval. Může ale být také neuspořádaná, například souřadnice bodů v rovině (v meteorologii) nebo v 3-D prostoru (geofyzika). Definice 2.3. Pro každý náhodný případ dostáváme funkci x : T R jako výsledek náhodného experimentu: x(t) := Xt(). Tato funkce se nazývá pozorování (trajektorie, realizace) procesu X. Poznámka 2.4. Obrázky 1.1-1.6 ilustrují trajektorie různých náhodných procesů (časových řad). Příklad 2.5 (Příklady časových řad). (1) Sinusovka s náhodnou amplitudou a fází (Obr. 1.1). (2) Binární proces házení mincí. (3) Náhodná procházka. (4) Proces větvení. 2.2. SYSTÉM DISTRIBUČNÍCH FUNKCÍ. Definice 2.6 (Konzistentní systém distribučních funkcí procesu X). Označme T := {t | t = [t1, t2, . . . , tn] Tn , ti = tj, pro i = j, n N}. Pro každé t T libovolné délky n N nechť Ft(x) značí adruženou distribuční funkci marginálního náhodného vektoru Xt vybraného ze stochastického procesu X = {Xt}tT v časových okamžicích t1, t2, . . . , tn. Systém {Ft}tT úplně popisuje stochastické chování 6 procesu X a nazývá se konzistentním systémem distribučních funkcí procesu X (viz následující tvrzení). Věta 2.7. Systém distribučních funkcí {Ft}tT z definice 2.6 se nazývá konzistentní protože má pro každé x = [x1, x2, . . . , xn]T Rn a n N následující dvě vlastnosti konzistence: (i) Ftp (xp) = Ft(x) pro každou permutaci p indexů {1, 2, . . . , n}. Zde značí tp, resp. xp vektory t, resp. x se složkami permutovanými podle permutace p. (ii) limxi Ft(x) = Ft(x1, . . . , xi-1,, xi+1, . . . , xn) =: =: Ft(i)(x(i)) pro každé i {1, 2, . . . , n}, kde t(i), resp. x(i)) značí vektor t, resp. x s vynechanou i-tou složkou. Definice 2.8. Stochastický proces se nazývá normální nebo gaussovský, když každá distribuční funkce Ft jeho konzistentního systému (t T) je sdruženou distribuční funkcí normálně rozloženého marginálního náhodného vektoru Xt. Definice 2.9. Jestliže X = {Xt} je časová řada, kde všechny náhodné veličiny Xt, t T jsou vzájemně nezávislé a stejně rozložené se střední hodnotou a rozptylem 2 , budeme psát X IID(, 2 ) 2.3. MOMENTOVÉ CHARAKTERISTIKY. Nyní zavedeme momentové funkce jako analogie střední hodnoty a kovarianční matice náhodného vektoru, který můžeme pokládat za speciální případ konečné časové řady (T = {1, 2, . . . , n}). Definice 2.10. Jsou-li dány stochastické procesy X = {Xt}tT a Y = {Yt}tT , oba definované na témž pravděpodobnostním prostoru, pak definujeme momentové funkce 1. a 2. řádu následovně: (1) střední hodnota X: X : T R, kde X (t) := EXt, pokud střední hodnoty existují pro všechna t T. 7 (2) autokovarianční funkce X: X : T × T R, kde X (r, s) := cov(Xr, Xs), pokud kovariance existují pro všechna r, s T. (3) rozptyl X: 2 X : T R+ , kde 2 X (t) := var(Xt) = cov(Xt, Xt) = X (t, t), pokud rozptyly existují pro všechna t T. (4) autokorelační funkce X: X : T × T [-1, 1], kde X (r, s) := 8 < : X (r,s) X (r,r) X (s,s) pro p X (r, r) p X (s, s) = 0 0 jinak, pokud korelace existují pro všechna r, s T. (5) vzájemná (křížová) kovarianční funkce X a Y : XY : T × T R, kde XY (r, s) := cov(Xr, Ys), pokud kovariance existují pro všechna r, s T. (6) vzájemná (křížová) korelační funkce X a Y : XY : T × T [-1, 1], kde XY (r, s) := 8 < : XY (r,s) X (r,r) Y (s,s) pro p X (r, r) p Y (s, s) = 0 0 jinak, pokud korelace existují pro všechna r, s T. 2.4. STRIKTNÍ A SLABÁ STACIONARITA. Definice 2.11. Časová řada X := {Xt | t Z} se nazývá striktně stacionární, jestliže každá distribuční funkce jeho konzistentního systému {Ft}tT, nezávisí na posunutí: Ft() Ft+h() pro každé t T a h Z. Definice 2.12. Časová řada X := {Xt | t Z} se nazývá (slabě) stationární, jestliže jsou splněny následující podmínky: (1) X má konečné druhé momenty: 2 X (t) < pro každé t Z. (2) X (r, s) = X (r + h, s + h) každé r, s, h Z. (3) X () X je konstantní funkce. Jestliže jsou splněny pouze podmínky (1) a (2), pak X se nazývá kovariančně stacionární. 8 Poznámka 2.13. (1) Zřejmě ze (2) vyplývá při r = s, že rozptyl stacionární časové řady je rovněž konstantní funkce: 2 x() 2 X . (2) Jestliže platí (3), pak z X (r, s) = EXrXs - (EXr)(EXs) = EXrXs - 2 x plyne, že (2) je ekvivalentní (a může tak být nahrazena) podmínkou: EXrXs = EXr+hXs+h pro každé r, s, h Z. Celkem vidíme, že všechny první a druhé momenty stacionární časové řady jsou invariantní k posunutí. Proto je slabá stacionarita také někdy označována jako stacionarita druhého řádu. Poznámka 2.14. Zřejmě podmínka (2) z definice 2.12 může být také ekvivalentně nahrazena modifikovanou podmínkou: (2') x(r, s) závisí pouze na rozdílu svých argumentů r - s. Můžeme proto autokovarianční i autokorelační funkci stacionární časové řady zavést jako funkci jen jedné proměnné: X (h) := X (t + h, t) X (h) := X (t + h, t) = X (t + h, t) X X = X (h) X (0) 2 X = X (t, t) = X (0), (2.1) kde t, h Z jsou libovolné. Věta 2.15. Každá striktně stacionární časová řada s konečnými druhými momenty je stacionární. Poznámka 2.16. Obecně stationarita nemá za následek striktní stacionaritu. Věta 2.17. Každá gaussovská stacionární časová řada je striktně stacionární. Definice 2.18. Časová řada X = {Xt} se nazývá bílý šum se střední hodnotou a rozptylem 2 , jestliže X (t) a 9 X (r, s) = ( 2 pro r = s 0 jinak . Píšeme X WN (, 2 ). Stacionární časová řada, která není bílým šumem, se někdy také nazývá barevný šum. Poznámka 2.19. Je snadné ověřit následující implikace: X IID(, 2 ) X WN (, 2 ) X je stacionární. Poznamenejme také, že žádná z opačných implikací neplatí. (viz též 2.16). 2.5. LINEÁRNÍ TRANSFORMACE ČASOVÝCH ŘAD. Definice 2.20. Řekneme, že časová řada X := {Xt} má ohraničené druhé momenty, jestliže existuje konstanta C R s vlastností E|Xt|2 C < pro každé t. Věta 2.21. Časová řada X := {Xt} má ohraničené druhé momenty právě když existují reálné konstanty X < a X < takové, že |X (t)| X a |X (t)| X pro každé t, tj. právě když X má ohraničenou střední hodnotu i rozptyl. Věta 2.22 (Hlavní věta o konvergenci). Jestliže := {k} 1 (tj.P k|k| < ) a U := {Ut} je časová řada s ohraničenými druhými momenty, pak nekonečná řada Xt : 2 = P k kUt-k konverguje podle kvadratického středu1 pro každé t a lineárně transformovaná časová řada X := {Xt} má rovněž ohraničené druhé momenty. 1tj. částečné součty konvergují dle kvadratického středu: E|Xt - PN k=-N kUt-k|2 0. Součet v tomto smyslu označujeme symbolem 2 = 10 Poznámka 2.23. Lineární transformace s koeficienty := {k}, která každou časovou řadu s ohraničenými druhými momenty transformuje opět na řadu s ohraničenými druhými momenty se nazývá stabilní. Z předchozí věty tedy vyplývá, že postačující podmínkou pro stabilitu lineární transformace je absolutní sumovatelnost posloupnosti jejích koeficientů. Proto se tato podmínka také někdy nazývá podmínkou stability. Důsledek 2.24. Jestliže 1 a U = {Ut} je stacionární se střední hodnotou U a autokovarianční funkcí U , pak platí (1) Xt 2 = P k kUt-k konverguje pro každé t (2) X := {Xt} je stacionární se střední hodnotou a autokovarianční funkcí: X = U X j j (2.2a) X (h) = X j X k U (h - j + k), h 0. (2.2b) Důsledek 2.25. 2 X = X (0) = X j X k U (k - j) (2.3a) 2 X 2 U ( X j |j|)2 (2.3b) Definice 2.26. Časová řada X, která je lineární transformací řady U jako ve větě 2.22 se také někdy nazývá časovou řadou (lineárně) generovanou řadou U. Jestliže je k = 0 pro každé k < 0, pak X se nazývá kauzální (kauzálně generovanou), neboť její hodnoty Xt nezávisí na budoucích hodnotách U , > t generující řady (jsou podmíněny jen hodnotami U , t). V takovém případě částečné součty při sumaci dle kvadratického středu nabudou tvaru PN k=0 kUt-k a samotnou řadu tak můžeme psát ve tvaru Xt 2 = P k=0 kUt-k. 11 Poznámka 2.27 (Operátor zpětného posuvu). Značíme pro k N: BUt := Ut-1, atd. rekurentně Bk Ut = B(Bk-1 Ut) := Ut-k B-1 Ut := Ut+1, atd. rekurentně B-k Ut := B-1 (B-(k-1) Ut) := Ut+k B0 Ut := Ut. Pak formálně místo Xt 2 = P k=0 kUt-k píšeme stručně Xt = (B)Ut, kde (B) := P k=0 kBk , neboť P k=0 kUt-k = P k=0 kBk Ut = ( P k=0 kBk )Ut je operátor získaný formálním dosazením operátoru B za proměnnou v řadě (polynomu) (z) = P k=0 kzk , která(ý) hraje roli přenosové charakteristiky lineární transformace (ta je totiž konvolučního typu). 2.6. ODHADY MOMENTOVÝCH FUNKCÍ. Definice 2.28. Nechť x = [x1, . . . , xn] je n pozorování (xt = Xt() pro t = 1, . . . , n) stacionární časové řady se střední hodnotou , rozptylem 2 , autokovarianční funkcí () a autokorelační funkcí (). Pak jejich odhady spočteme takto: b := 1 n Pn j=1 xj . . . odhad ; b(h) := 1 n Pn-h j=1 (xj+h - b)(xj - b), 0 h n - 1, b(h) := b(-h), -(n - 1) h < 0 . . . odhad (h); b2 := b(0) . . . odhad rozptylu; b(h) := b(h) b(0) , -(n - 1) h n - 1 (viz rovnici (2.1)) . . . odhad autokorelační funkce v případě b(0) = 0, jinak b(h) := 0. Věta 2.29. Nechť X := [X1, . . . , Xn] marginální vektor v X korespondující s vektorem pozorování x. Potom jak matice bn := [b(i - j)]i,j = 2 6 6 4 b(0) b(1) . . . b(n - 1) b(1) b(0) . . . b(n - 2) . . . . . . . . . . . . . . . . . . . . . . . . . . b(n - 1) b(n - 2) . . . b(0) 3 7 7 5 , která je odhadem varianční matice varX, tak i matice bRn := bn b(0) , která je odhadem korelační matice (X), jsou symetrické a pozitivně semidefinitní. 12 0 100 200 300 400 500 600 700 800 900 1000 -3 -2 -1 0 1 2 3 White noise WN(0,1) Gaussovský bílý šum WN(0,1) 0 5 10 15 20 25 30 35 40 45 50 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Autocorrelation function of white noise WN(0,1) Autokorelační funkce gaussovského bílého šumu WN(0,1) Poznámka 2.30. (1) Odhad je spolehlivý pouze pro n > 50 a h < n 4 . 13 (2) Odhad b(h) je vychýlený (Eb(h) = (h)), protože součet pozorování je dělen n a nikoliv počtem stupňů volnosti n-1-h. Poznamenejme, že věta 2.29 neplatí pro nevychýlený odhad (matice bn ztratí očekávanou vlastnost pozitivní semidefinitnosti). V každém případě je ale odhad b(h) asymptoticky nevychýlený v tom smyslu, že Eb(h) (h) pro n . Navíc je také konzistentní podle kvadratického středu v tom smyslu, že E|b(h) - (h)|2 0 pro n , kde konvergence je dokonce rychlejší než v případě nevychýleného odhadu. (3) Z algebraického hlediska b(h) lze psát ve tvaru skalárního součinu b(h) = 1 n x0, xh , kde xh := [0, . . . , 0 | {z } h , x1 - b, . . . , xn - b, 0, . . . , 0 | {z } n-1-h ]. Vektor x0 tak reprezentuje the původní vektor pozorování (doplněný na konci n-1 nulami), zatímco xh je jeho kopie posunutá (zpožděná) o h. Zřejmě x0 2 = xh 2 = Pn j=1|xj - b|2 . Ze Schwarzovy nerovnosti dostáváme | x0, xh | x0 2 , což vede na |b(h)| 1 n x0 2 = 1 n x0, x0 = b(0). Odtud vidíme, že odhad autokorelační funkce zachovává její přirozenou vlastnost |b(h)| 1. (4) Vzhledem ke (3) je možno odhad b(h) interpretovat geometricky jako kosimus úhlu mezi původním a posunutým vektorem pozorování x0 a xh, což lze interpretovat jako míru jejich lineární závislosti (podobnosti): nula znamená ortogonalitu (úplnou lineární nezávislost=žádná korelace mezi oběma vektory), 1 odpovídá lineární závislosti (úplná korelace: jeden z nich je skalárním násobkem druhého). (5) Trend je charakterizován korelacemi na velkou vzdálenost, což se projevuje pomalým poklesem (h) při h . Periodická složka se projeví oscilatorickým chováním b(h) s dominantní periodou této složky, případně směsí takových složek. 14 3. Predikce v časových řadách 3.1. PRINCIP ORTOGONÁLNÍ PROJEKCE. Definice 3.1 (Prostor L2(, A, P)). L2 := L2(, A, P) definujeme jako množinu všech (komplexních) náhodných veličin nad týmž pravděpodobnostním prostorem (, A, P), které mají konečné druhé momenty (resp. rozptyly - viz dále 3.5), tj. L2(, A, P) := {X | X náhodná veličina nad (, A, P), E|X|2 < }. Poznamenejme, že do tohoto prostoru zahrnujeme také všechny konstanty z C, které považujeme za náhodné veličiny s nulovým rozpty- lem. Věta 3.2. L2(, A, P) je Hilbertův prostor se skalárním součinem X, Y = EXY a normou X 2 = p X, X = p E|X|2. Důkaz. L2(, A, P) je obdobou funkcionálního prostoru L2(I) tvořeného funkcemi absolutně integrovatelnými v kvadrátu na intervalu I R. Totiž E|X|2 = R |X()|2 dP(), takže namísto s Lebesgueovým integrálem pracujeme s obecnějším pojetím integrálu, kde Lebesgueova míra je nahrazena pravděpodobnostní mírou P: * Skalární součin X, Y existuje a je konečný pro každé X, Y L2(, A, P), jak snadno nahlédneme z nerovnosti 4|XY | = (|X| + |Y |)2 - (|X| - |Y |)2 (|X| + |Y |)2 + (|X| - |Y |)2 = = 2(|X|2 + |Y |2 ), odkud užitím |Y | = |Y | dostáváme |XY | 1 2 (|X|2 + |Y |2 ), takže |EXY | Z |X()Y ()| dP() 1 2 Z |X()|2 dP() + Z |Y ()|2 dP() < . 15 L2(, A, P) je vektorovým prostorem. Je uzavřený na násobení skaláry c C, neboť E|cX|2 = |c|2 E|X|2 < . Uzavřenost vzhledem ke sčítání plyne z: |X + Y |2 (|X| + |Y |)2 = |X|2 + 2|XY | + |Y |2 E|X + Y |2 E|X|2 + 2E|XY | + E|Y |2 < . * Ověření, že L2(, A, P) je úplný, neboli Hilbertův prostor, je složitější, ale provádí se opět zcela analogicky jako v případě funkcionálního prostoru L2(I). Podrobnosti lze nalézt například v monografii [BD93, §2.10]. Důsledek 3.3 (Schwarzova nerovnost). |EXY |, |EXY | X 2 Y 2 = p E|X|2 p E|Y |2, X, Y L2(, A, P). Důsledek 3.4. X L2(, A, P) EX existuje. Důkaz. |EX| = |E(1.X)| p E|1|2 | {z } 1 p E|X|2 = p E|X|2 < . Důsledek 3.5. X, Y L2(, A, P) X - EX, Y - EY L2(, A, P) a X - EX, Y - EY = E(X - EX)(Y - EY ) = cov(X, Y ) existuje a splňuje Schwarzovu nerovnost |cov(X, Y )| p E|X - EX|2 p E|Y - EY |2 = X Y . 16 Důsledek 3.6. (X, Y ) = ( cov(X,Y ) X Y pro xy = 0 0 pro xy = 0 je tzv. korelační koeficient náhodných veličin X a Y , pro nějž platí |(X, Y )| 1 a speciálně -1 (X, Y ) 1 v případě reálných náhodných veličin X a Y . Věta 3.7 (Věta o ortogonální projekci: [BD93], §2.3). Je-li L uzavřený lineární podprostor Hilbertova prostoru H (například H = L2) a x H, pak (i) existuje jediný prvek bx L takový, že x - bx = inf yL x - y (*) a (ii) nějaký prvek má minimální vlastnost (*) právě když x - bx L, tj. x - bx y (neboli x - bx, y = 0) pro každé y L. Píšeme bx = PLx, kde PL je tzv. operátor ortogonální projekce na podprostor L. Definice 3.8. Nechť X0 := 1, X1, . . . , Xn L2 a L := L(X0, X1, . . . , Xn) je lineární podprostor v L2 generovaný náhodnými veličinami X0, X1, . . . , Xn. Je-li X L2 nějaká další náhodná veličina, pak bX = PLX se nazývá nejlepší2 lineární predikce náhodné veličiny X pomocí náhodných veličin X0, X1, . . . , Xn. Poznámka 3.9. Operátor PL je korektně definován, neboť L má konečnou dimenzi (dim L n + 1) a je tak automaticky uzavřený. Protože bX L(X0, X1, . . . , Xn) existuje vektor 0 := [0, 1, . . . , n]T Rn+1 takový, že bX = 0 X0 |{z} 1 +1 + + nXn. (3.1a) 2rozumí se ve smyslu minimální střední kvadratické chyby. 17 I když bX je jediný, nemusí být 0 obecně jednoznačně určen. Je tomu tak jen když X0, X1, . . . , Xn jsou lineárně nezávislé, tj. tvoří bázi prostoru L. Věta 3.10. Vektor 0 ze vztahu (3.1a) se spočte jako řešení následujícího systému n + 1 lineárních rovnic (tzv. normální systém rovnic), kde i = 0, 1, . . . , n: X0, Xi 0 + X1, Xi 1 + + Xn, Xi n = X, Xi (3.1b) neboli dle věty 3.2 E(X0Xi)0 + E(X1Xi)1 + + E(Xn, Xi)n = E(XXi). (3.1c) Důkaz. Dle 3.7: bX = PLX X - bX Y pro každé Y L X - bX Xi pro každé i = 0, 1, . . . , n X - bX, Xi = 0 pro každé i = 0, 1, . . . , n bX, Xi = X, Xi pro každé i = 0, 1, . . . , n Pn j=0 jXj, Xi = X, Xi pro každé i = 0, 1, . . . , n Pn j=0 j Xj, Xi = X, Xi pro každé i = 0, 1, . . . , n, což je právě (3.1b). Příklad 3.11 (n = 0). L = L(X0) = L(1) = R je lineární podprostor všech (reálných) konstant v L2. Pak (3.1b) se redukuje jen na jednu rovnici: E(X0X0 | {z } 1 )0 = E(X X0 |{z} 1 ). Odtud máme 0 = EX a tedy bX = 0X0 = EX. Tedy EX je nejlepší lineární predikce náhodné veličiny X pomocí konstanty. Dosažená minimální střední kvadratická chyba je E(X - EX)2 = varX. Důsledek 3.12. Jestliže EX = EX1 = . . . EXn = 0, pak bX = PLX = PL(X1,...,Xn)X, 0 = 0 a (3.1c) přejde v systém n rovnic, kde i = 1, . . . , n: E(X1Xi)1 + + E(Xn, Xi)n = E(XXi). (3.1d) 18 Položíme-li X := [X1, . . . , Xn]T , můžeme jej přepsat do maticového tvaru (varX) = cov(X, X)T , (3.1e) kde := [1, . . . , n]T Rn . Důkaz. První rovnice (i = 0) v (3.1c) bude tvaru: E(1.1) | {z } 1 0 + E(X1.1) | {z } 0 1 + + E(Xn.1) | {z } 0 n = E(X.1) | {z } 0 . Odtud dostáváme 0 = 0, takže ve zbývajících rovnicích pro i = 1, . . . , n bude E(X0X1)0 = 0, což je právě systém normálních rovnic pro PL(X1,...,Xn)X . Důsledek 3.13 (Nejlepší lin. predikce pro stacionární časové řady). Nechť X := {Xt} je stacionární časová řada s nulovou střední hodnotou X = 0 a autokovarianční funkcí X . Pro dané k N nechť bX (k) t+k značí nejlepší (k-krokovou) lineární predikci náhodné veličiny Xt+k pomocí n předchozích veličin Xt, Xt-1, . . . , Xt+1-n ve tvaru bX (k) t+k = Pn j=1 (k) n,jXt+1-j. Pak vektor (k) n := [ (k) n1 , . . . , (k) nn ]T nalezneme řešením tzv. Yule-Walkerova systému lineárních rov- nic: 2 6 6 6 4 (0) (1) (n - 1) (1) (0) (n - 2) ... ... ... (n - 1) (n - 2) (0) 3 7 7 7 5 | {z } n 2 6 6 6 6 4 (k) n1 (k) n2 ... (k) nn 3 7 7 7 7 5 | {z } (k) n = 2 6 6 6 4 (k) (k + 1) ... (k + n - 1) 3 7 7 7 5 | {z } (k) n (3.2) Důkaz. Plyne z 3.12 po záměnách: X ; Xt+k, Xj ; Xt+1-j, j ; (k) nj a tedy pak varX = [cov(XiXj)] ; [cov(Xt+1-i, Xt+1-j)] = [(t + 1 - i - (t + 1 - j))] = [(j - i)] = n a podobně 19 cov(X, X)T ; [E(Xt+k, Xt+1-j)] = = [(t + k - (t + 1 - j))] = [(k + j - 1)] = (k) n . Věta 3.14 (Střední kvadratická chyba). Střední kvadratická chyba v (k) n := E(Xt+k -X (k) t+k)2 nejlepší k-krokové lineární predikce se spočte dle vztahu: v(k) n = (0) - nX j=1 (k) nj (k + j - 1) = (0) - T n (k) n . (3.3) Důkaz. v (k) n = E(Xt+k - X (k) t+k)2 = Xt+k - X (k) t+k 2 = Xt+k - X (k) t+k, Xt+k - X (k) t+k = Xt+k - X (k) t+k, Xt+k Xt+k - X (k) t+k, X (k) t+k | {z } 0 = Xt+k, Xt+k - X (k) t+k, Xt+k = Xt+k 2 - Pn j=1 (k) nj Xt+1-j, Xt+k = Xt+k 2 - Pn j=1 (k) nj Xt+1-j, Xt+k = E(X2 t+k) - Pn j=1 (k) nj E(Xt+1-jXt+k) = (0)- Pn j=1 (k) nj (t+1-j -(t+k)) = (0)- Pn j=1 (k) nj (1-j -k) = (0) - Pn j=1 (k) nj (k + j - 1). Definice 3.15. Nechť X, Y, X1, . . . , Xn L2 a bX a bY jsou nejlepší lineární predikce náhodných veličin X a Y pomocí X1, . . . , Xn. Pak (X, Y | X1, . . . , Xn) := (X - bX, Y - bY ) se nazývá parciální korelační koeficient náhodných veličin X a Y při daných X1, . . . , Xn. Interpretace: parciální korelace mezi X a Y při daných X1, . . . , Xn je korelace mezi X a Y očištěná o tu svoji část, která je zprostředkována náhodnými veličinami X1, . . . , Xn. Definice 3.16. Nechť X = {Xt | t Z} je stacionární časová řada s autokorelační funkcí X (). Pak parciální autokorelační funkci 20 (pacf) X () časové řady X definujeme následovně: X (0) = X (0) = 1, X (1) = X (1) = (Xt+1, Xt) X (h) = (Xt+h, Xt | Xt+1, . . . , Xt+h-1) for h 2. Věta 3.17. Jestliže X = {Xt | t Z} je stacionární časová řada se střední hodnotou X = 0 a parciální autokorelační funkcí X (), pak X (n) = n,n pro n 1, kde n = [n,1, . . . , n,n]T je řešením problému nejlepší lineární predikce. Důkaz. See [BD93, §5.2]. X (h) lze počítat rekurzivně užitím mezivýsledných vztahů (3.4b) dále uvedeného Durbin-Levinsonova algoritmu 3.19. Jiný postup počítá (k) n,n z (3.2) pomocí Cramerova pravidla podle následujícího důsledku. Bohužel tato metoda není výpočetně příliš efektivní pro velká n. Důsledek 3.18. Je-li matice n := [x(i - j)]n i,j=1 regulární, pak X (n) = n,n = det n detn , kde n := [(:, 1), . . . , (:, n - 1), n] a n = [X (1), . . . , X (n)]T . 21 3.2. DURBIN-LEVINSONŮV A INOVAČNÍ ALGORITMUS. Věta 3.19 (Durbin-Levinsonův algoritmus pro rekur. výpočet n). Nechť X = {Xt | t Z} je stacionární časová řada s nulovou střední hodnotou X = 0 a autokovarianční funkcí X (h) 0 pro h . Jestliže bXn+1 = n,1Xn + + n,nX1 je nejlepší lineární predikce, pak pro koeficienty n,j a střední kvadratickou chybu vn = E|Xn+1 - bXn+1|2 platí následující rekurentní vztahy. Počáteční krok pro n = 0: v0 = X (0). (3.4a) Rekurentní krok pro n > 0: n,n = X (n) =0 pro n=1 z }| { n-1X j=1 n-1,j X (n - j) v-1 n-1, (3.4b) vn = vn-1(1 - |n,n|2 ), (3.4c) 2 6 6 6 4 n,1 n,2 ... n,n-1 3 7 7 7 5 = 2 6 6 6 4 n-1,1 n-1,2 ... n-1,n-1 3 7 7 7 5 - n,n 2 6 6 6 4 n-1,n-1 n-1,n-2 ... n-1,1 3 7 7 7 5 pro n > 1. (3.4d) Důkaz. Viz [BD93, str. 162]. Představme ještě jiný algoritmus pro výpočet jednokrokové nejlepší lineární predikce. Jedná se o tzv inovační algoritmus (IA), který může být na rozdíl od Durbin-Levinsonova algoritmu (DLA) aplikován na libovolnou i nestacionární časovou řadu s nulovou střední hodnotou X = 0. Tento algoritmus není nic jiného než souřadnicový přepis Gram-Schmidtova ortogonalizačního procesu. 22 Věta 3.20 (Inovační algoritmus pro rekurentní výpočet n). Nechť X = {Xt | t Z} je časová řada s nulovou střední hodnotou X = 0 s maticí smíšených momentů [i,j]n i,j=1, i,j := E(XiXj), která je regulární pro každé n. Pak jednokrokové nejlepší lineární predikce bXn+1 := bX (1) n+1, n 0, a jejich střední kvadratické chyby vn := v (1) n spočteme rekurentně pomocí následujících vztahů: bXn+1 = ( 0 pro n = 0 Pn j=1 n,j(Xn+1-j - bXn+1-j) pro n 1 , (3.5a) užitím v0 = (1, 1) (3.5b) n,n-k = v-1 k ((n + 1, k + 1) - k-1X j=0 k,k-jn,n-jvj) (3.5c) vn = (n + 1, n + 1) - n-1X j=0 2 n,n-jvj. (3.5d) postupně v pořadí: v0; 1,1, v1; 2,2, 2,1, v2; . . . 23 4. Modelování a odhad parametrů 4.1. PŘEDZPRACOVÁNÍ. 4.1.1. Ošetření specifických efektů. a) Volba okamžiků pozorování: * jsou přímo diskrétní svou povahou, např. úroda obilí za jednotlivé roky. * vznikají diskretizací spojité časové řady, např. teplota v danou denní dobu na daném místě. * vznikají akumulací (agregací) hodnot za určité období, např. denní množství srážek, roční výroba závodu. Místo akumulace se někdy provádí průměrování. Je-li dána možnost volby, je třeba jí věnovat pozornost: * málo bodů unikne charakteristický rys řady. * mnoho bodů zvýší se výpočetní náročnost. * ekvidistantní diskretizace zpravidla usnadní numerické zpracování, ale neumožňuje adaptivně měnit hustotu diskretizace v závislosti na lokálním charakteru řady. * při agregaci se mohou porušit vlastnosti původní řady. b) Problémy s kalendářem: * různá délka kalendářních měsíců. * 4 nebo 5 víkendů v měsíci. * různý počet pracovních dnů v měsíci. * pohyblivé svátky: např. svátek na začátku měsíce sníží prodej potravin za tento měsíc, ale zvýší jej za předchozí v důsledku efektu předzásobení. Příklad `očištění' časové řady např. od proměnlivé délky měsíce: zavedeme `standardní' měsíc o délce 30 dnů a pak údaj třebas o produkci za leden přenásobíme korekčním faktorem 30 31 . 24 c) Problémy s nekompatibilitou jednotlivých měření: Příklad: hodnota nějakého ukazatele se jeden rok týká např. 85 podniků, další rok jen 82 apod. d) Problémy s délkou časových řad: * zvětšení počtu měření (např. půlením časových intervalů mezi body) nemusí vždy znamenat zvětšení množství informace. * někdy se mohou objevit protichůdné tendence: metoda zpracování vyžaduje delší řadu, ale na druhé straně řada vzniklá dlouhodobým sledováním může měnit charakteristiky svého modelu v čase obtíže s konstrukcí mo- delu. 4.1.2. Stabilizace rozptylu. Většina běžně užívaných modelů předpokládá alespoň kovariančně stacionární chování analyzované řady X := {Xt}, tj. zejména konstantní rozptyl 2 X (tzv. homoskedasticita). Pokud je tato podmínka porušena (tzv. heteroskedasticita), je třeba buď užít vhodný model a nebo řadu transformovat na řadu s konstantním rozptylem. Takové transformaci říkáme stabilizace rozptylu. Obvykle předpokládáme exponenciální model závislosti směrodatné odchylky na střední hodnotě: X (t) = 0X (t) . V praxi se nejčastěji vyskytují 2 případy: * = 0: řada je již homoskedastická a není třeba ji transformovat; * = 1: směrodatná odchylka řady závisí lineárně na střední hod- notě. Existují postupy, které umožňují odhadnout parametr z pozorované řady. Heteroskedasticitu pak odstraníme jednou z níže uvedených transformací, kde zvolíme parametr = 1 - . 25 Box-Coxova transformace pro řadu Xt > 0: Yt := ( X t -1 pro = 0 ln(Xt) pro = 0 (lineární nárůst při = 1) . Mocninná transformace pro řadu Xt > 0: Yt := ( X t pro = 0 ln(Xt) pro = 0 (lineární nárůst při = 1) . Poznámka 4.1. (1) Jestliže není splněna podmínka Xt > 0 nebo pozorované hodnoty jsou blízké nule, nelze uvedené transformace přímo použít. V takovém případě je možno řadu vhodně posunout do kladných hodnot. Tento posuv však představuje natolik umělý zásah do stochastického chování originální řady, že může vést k jejímu znehodnocení a nevěrohodnosti výsledného modelování. (2) Mocninná transformace není spojitá pro 0, takže je nutno se vyhýbat malým nenulovým hodnotám parametru . (3) Postup odhadování parametru , resp. lze nalézt např. v [Cip86, s.144-145]. 4.1.3. Identifikace periodických komponent. Identifikace je založena na rozkladu T-periodické funkce x(t) do její Fourierovy řady: x(t) = X k=- cke i2kt T = a0 2 + X k=1 Ak cos 2kt T - k , 26 kde se pak snažíme identifikovat energeticky nejsilnější harmonické komponenty, tj. takové, které mají velkou hodnotu kvadrátu své amplitudy A2 k, což je veličina úměrná energii (nebo střednímu výkonu) této složky na intervalu délky T. K tomu slouží tzv. periodogram: Periodogram = odhad energetické spektrální hustoty reálné T-periodické funkce x(t), kde za ck, k = 1, . . . , m, m := N-1 2 , dosadíme odhad bck spočtený pomocí DFTN (operátor Diskrétní Fourierovy transformace -- viz dále) po ekvidistantní diskretizaci jedné periody x(t): T = Nt, xn = x(nt), n = 0, 1, . . . , N, x := [x0, . . . , xN-1]. Periodogram je tedy posloupnost hodnot {Ik}m k=1, kde Ik = I(k) = 2T|bck|2 = 2Nt 1 N Xk 2 = 2 N t|Xk|2 , k = 1, . . . , m (4.1) a X := [X1, . . . , XN ] je transformací vektoru x pomocí DFTN : Xk = N-1X t=0 xte-i 2kt N = NX t=1 xte-i 2kt N . Druhá rovnost je zde důsledkem N-periodicity x0 = xN . Z hlediska kvalitativních závěrů nezáleží na multiplikativní konstantě, proto bez újmy na obecnosti můžeme předpokládat t = 1 a dostáváme tak výsledný vztah pro hodnoty Ik periodogramu ve tvaru Ik = I(k) = 2 N |Xk|2 = (4.2) = 2 N 8 < : NX t=1 xt cos 2kt N !2 + NX t=1 xt sin 2kt N !2 9 = ; , k = 1, . . . , m. (4.3) 27 Velká hodnota Ik indikuje velkou energii k-té harmonické komponenty, tj. silné zastoupení složky xk(t) = ckeikt + c-ke-ikt = Ak cos(kt - k) v rozvoji x(t) do Fourierovy řady (k := 2k T ). Věta 4.2. Ik = 2 X |h| x) = mX j=1 1-jx>0 (-1)j-1 m j (1 - jx)m-1 , 0 < x < 1, přičemž P(W > x) m(1 - x)m-1 je dobrá aproximace pro m 50. Nechť gF (1 - ) značí (1 - )-kvantil rozdělení statistiky W (tabelováno), tj. P(W gF (1 - )) = 1 - . Pak H0 zamítáme s rizikem , pokud W > gF (1 - ). Je-li v takovém případě Ik0 = maxk=1,...,m Ik, pak k0-tou harmonickou komponentu považujeme za významnou. Další významnou harmonickou komponentu zjistíme z periodogramu délky m - 1, kde vynecháme Ik0 , a tak pokračujeme dále, dokud není hypotéza H0 přijata. 4.1.5. Siegelův test periodicity. Tento test je vhodnější než Fisherův v případě většího množství harmonických komponent. Siegelova statistika je tvaru T = mX k=1 (Yk- gF (1-))+ , 0 < < 1 (doporučená hodnota je = 0, 6). 29 H0 zamítáme s rizikem , jestliže T > t(1 - ), kde t(1 - ) je (1 - )-kvantil rozdělení T (tabelováno). Poznámka 4.1 (Praktická doporučení). (1) Vzhledem k povaze nulové hypotézy ve výše uvedených testech by v v testovaných datech neměl být obsažen trend. Doporučuje se protom alespoň hrubé předběžné odstranění dlouhodobého trendu. Jeho přítomnost totiž snižuje citlivost testů. (2) Vzhledem k povaze rozvoje do Fourierovy řady musí být N dělitelné délkami period všech harmonických složek, které se snažíme detekovat (obvykle stačí, aby N bylo dělitelné délkoou dominantní periody). V případě potřeby je proto třeba od začátku řady ubrat co nejmenší počet pozorování a snížit tak celkovou délku n na N = n - r, které vyhovuje. Činíme tak ovšem jen pro účel identifikace periodických komponent. (3) Požadavek (2) je obtížné splnit, jestliže délka dominantní periody není předem známa. V takovém případě postupně zkoušíme r = 0, 1, . . . a vybereme N = n - r, pro nějž jsou v periodogramu nejvyšší a nejostřejší špičky (lokální maxima), z nichž je co nejvíce vůči sobě navzájem harmonických (perioda jedné komponenty dělí periodu jiné komponenty). 4.2. DEKOMPOZIČNÍ MODEL. Aditivní model (AM): Xt = Trt + Pt z }| { Szt + Ct | {z } Dt +Et, Multiplikativní model (MM): Xt = Trt. Pt z }| { Szt.Ct | {z } Dt .Et, kde značí Trt . . . dlouhodobý trend 30 Szt . . . sezónní komponenta (perioda je předem známa) Ct . . . cyklická komponenta (perioda není předem známa) Pt . . . celková periodická komponenta Dt = X (t) . . . deterministická komponenta (střední hodnota Et . . . stacionární náhodná komponenta (obvykle bílý šum nebo IID): E = 0 pro AM, resp. E = 1 pro MM. AM lze použít jen pro homoskedastická data, neboť X (t) = (Dt + Et) = (Et) = E je konstantní vzhledem ke stacionaritě Et. Naopak MM lze použít jen pro lineárně heteroskedastická data, neboť X (t) = (Dt.Et) = Dt(Et) = X (t)E závisí lineárně na střední hodnotě X (t). V takovém případě lze MM převést na AM logaritmickou transformací (srov. s 4.1.2). 4.2.1. Parametrizované modely. Označení . t = (t1, . . . , tn)T , ti = tj pro i = j . . . vektor "času", X := Xt = (Xt1 , . . . , Xtn )T . . . náhodný vektor populace vybrané z časové řady, x := xt = (x1, . . . , xn)T , xi = Xti () . . . pozorování vektoru populace X. A. Parametrické metody (P). Tyto metody používáme v případě, že je znám analytický tvar Dt = D(t; 1, . . . , p) v aditivním modelu Xt = Dt + Et (multiplikativní model převedeme na aditivní logaritmováním) s neznámými parametry = (1, . . . , p)T , které odhadneme minimalizací výrazu N(1, . . . , p) = D(t; 1, . . . , p) - x 2 ve vhodné normě . Zpravidla se používá euklidovská norma a klasická metoda nejmenších čtverců: N(1, . . . , p) = nX i=1 |D(ti, 1, . . . , p) - xi|2 . 31 Lokální extrém hledáme řešením systému rovnic N(1, . . . , p) j = 0 pro j = 1, . . . , p s ověřením podmínky pro lokální minimum (Hessova matice je pozitivně definitní): J = 2 N(1, . . . , p) i j i,j > 0. Pokud D závisí lineárně na 1, . . . , p, pak tento systém je systémem p lineárních rovnic pro p neznámých 1, . . . , p a dostáváme klasický lineární regresní model. V případě nelineární regrese použijeme pro minimalizaci N(1, . . . , p) vhodnou optimalizační metodu, například Optimization Toolbox v prostředí systému MATLAB. Lineární regresní model: D(t; 1, . . . , p) = 1(t)1 + + p(t)p, kde i jsou dané. X = D + E, E = (Et1 , . . . , Etn )T IID(0, 2 ), kde D = [1(t), . . . , p(t)], j(t) = (j(t1), . . . , j(tn))T je matice s plnou sloupcovou hodností p. Pro vektor pozorování přejde výše uvedený model v x = D + e, e = (e1, . . . , en)T , kde ei = Eti (). Klasické dekompoziční metody * Metoda malého trendu předpokládá v každé sezónní periodě přibližně konstantní trend, který se odhadne jako průměr. Po jeho odečtení odhadneme sezónní složku sezónním průměrováním (viz [BD93, str. 21]). 32 ˇ Metoda sezónních klouzavých průměrů využívá pro odhad trendu obyčejných klouzavých průměrů s délkou rovnou sezónní periodě. Dále postupujeme jako v předchozím případě (viz [BD93, str. 22]). * lineární regresní model (ozn. polfs), v němž trend je modelován polynomem (pol) a periodická složka zkrácenou Fourierovou řadou (fs=Fourier Series) v goniometrickém tvaru nazývaný též metoda skrytých period, do níž zařadíme pouze významné harmonické složky s periodou Tj dle perio- dogramu: Trt = 1tp + 2tp-1 + + p+1, 0 p (4.4a) Pt = q X j=1 aj cos 2t Tj + bj sin 2t Tj . (4.4b) V modelu tedy odhadujeme (p+1)+2q neznámých parametrů i, i = 1, 2, . . . , p + 1 a aj, bj, j = 1, 2, . . . , q. * regresní model (ozn. polper), v němž trend je opět modelován polynomem (4.4a) a u periodické složky Pt (per) jsou její hodnoty v jedné její základní (sezónní) periodě přímo odhadovanými parametry: sk := Pk pro k = 2, . . . , s a s1 := P1 určíme z dodatečné podmínky 1 s Ps k=1 sk = 0: Celkem tedy v modelu odhadujeme p + 1 + s - 1 = p + s parametrů a pak klademe Pt = sk, t = (j -1)s+k (tj., je-li Pt k-tou hodnotou v j-té základní periodě). 4.2.2. Filtrační techniky. Za vhodných předpokladů konstruujeme lineární nebo nelineární operátor S (smoother = vyhlazovač) takový, že S({Xt}) =: { bDt} {Dt} je dobrý odhad Dt v tom smyslu, že { bDt} {Dt} ve vhodné konvergenci při n , tj. odhad je asymptoticky přesný. 33 Běžné vyhlazovací techniky * Metoda klouzavých průměrů (lineární konvoluční filtr): Yt = q X =-q w Xt+ = q X =-q w Dt+ | {z } bDt:=S({Dt}) + q X =-q w Et+ | {z } bEt:=S({Et}) , kde požadujeme: (1) Co nejmenší zkreslení Dt, tj. Dt Pq =-q w Dt+ . Standardní požadavek je, aby bez zkreslení byly filtrem zachovány polynomy dostatečně vysokého stupně. Je-li Dt po částech polynomického tvaru, bude pak aproximace bDt velmi dobrá. Odpověď dává následující věta: Věta 4.2. Platí Pt = Pq =-q Pt+ pro každý polynom Pt = c0 +c1t+ +cktk nejvýše k-ho stupně právě když váhy splňují následující 2 podmínky q X =-q w = 1 q X =-q r w = 0 pro r = 1, . . . , k (4.5) (2) Co největší redukci rozptylu náhodné složky Et, tj. 0 Pq =-q w Et+ . Mírou pro tento účel je tzv. redukční intenzita klouzavých průměrů definovaná jako podíl rozptylů 2 Y /2 X za předpokladu, že Et WN (0, 2 ). Spočtěme nejprve rozptyl 2 Y : 2 Y = var bEt = var( Pq =-q w Et+ ) = Pq =-q w2 var(Et+ ) = 2 Pq =-q w2 = 2 w 2 . Odtud 2 Y /2 X = var bEt/varEt = 2 w 2 /2 = w 2 . Požadujeme tedy w 2 < 1 co nejmenší. Dá se ukázat, 34 že z tohoto pohledu vyhovují nejlépe obyčejné klouzavé průměry w = 1 2q+1 , kde w 2 = Pq =-q 1 (2q+1)2 = 2q+1 (2q+1)2 = 1 2q+1 . Ty však zase nejsou příliš optimální z hlediska požadavku (1). Zachovávají totiž pouze polynomy nejvýše stupně 1, tj. (po částech) lineární průběh Dt a nehodí se tedy na řady s vyšší dynamikou ve střední hodnotě. Klasickou technikou hledání vah je postupná polynomiální regrese, kdy postupně zleva prokládáme segmenty dat délky 2q + 1 regresní polynom zvoleného stupně. Jsou-li data ekvidistantní, pak takto získané váhy nezávisí na poloze segmentu a jsou tedy korektně definovány. Hledání optimálních vah zachovávajících polynomy do zadaného stupně při co nejlepší redukční intenzitě vedou na úlohu kvadratického programování (minimalizujeme w 2 =Pq =-q w2 ) při lineárních omezeních daných rovnicemi (4.5). Příkladem kvalitně navržených vah jsou tzv. Spencerovy 15-ti bodové váhy: [w0, . . . , w7] = 1 320 [74, 67, 46, 21, 3, -5, -6, -3], kde w- = w pro || 7. Tyto váhy zachovávají polynomy až do stupně 3 při redukční intenzitě 0.1926, což je velmi dobrý výsledek ve srovnání s dolní mezí 1 15 = 0.0667 danou obyčejnými klouzavými průměry délky 15. * Exponenciální vyrovnávání * Jádrové vyhlazování * Waveletové vyhlazování 35 4.3. BOX-JENKINSOVY MODELY. 4.3.1. ARMA modely pro stacionární časové řady. Definice 4.3 (ARMA proces). Stochastický proces X = {Xt | t Z} se nazývá ARMA procesem řádu p, q (0 p, q < ), píšeme X ARMA(p, q), jestliže Xt = 1Xt-1 + 2Xt-2 + + pXt-p | {z } Autoregresní část (AR) + Zt + 1Zt-1 + 2Zt-2 + + qZt-q | {z } Část klouzavých součtů (MA=Moving Average) , (4.6a) kde Z := {Zt | t Z} WN (0, 2 ) a p = 0, q = 0, = 0. Přepíšeme-li (4.6a) do ekvivalentního tvaru Xt - 1Xt-1 - 2Xt-2 - - pXt-p = Zt + 1Zt-1 + 2Zt-2 + + qZt-q, (4.6b) lze použít též zkrácený zápis dle poznámky 2.27: (B)Xt = (B)Zt (4.6c) kde při 0 = 0 = 1 je (z) = 1 - 1z - - pzp , (z) = 1 + 1z + + qzq , z C. Poznámka 4.4. V předchozí definici bez újmy na obecnosti předpokládáme 0 = 1, neboť v opačném případě stačí nahradit bílý šum {Zt} WN (0, 2 ) šumem {0Zt} WN (0, (0)2 ). 36 Poznámka 4.5 (Speciální případy). a) Autoregresní proces (AR proces): X AR(p) = ARMA(p, 0) : (B)Xt = Zt, neboť (z) 1. Tedy Xt = 1Xt-1 + 2Xt-2 + + pXt-p + Zt. (4.6d) V tomto případě připouštíme i p = za předpokladu, že := {j} j=1 1. b) Proces klouzavých součtů (MA proces): X MA(q) = ARMA(0, q) : Xt = (B)Zt, neboť (z) 1. Tedy Xt = Zt + 1Zt-1 + 2Zt-2 + + qZt-q. (4.6e) V tomto případě připouštíme i q = za předpokladu, že := {j} j=1 1. c) Bílý šum (jediný proces, který je současně AR i MA): X ARMA(0, 0) = AR(0) = MA(0) = WN (0, 2 ) : Xt = Zt. d) Smíšený ARMA proces: X ARMA(p, q), 0 < p, q < : Netriviální směs autoregrese a klouzavých součtů. Definice 4.6. Buď X = {Xt | t Z}, X ARMA(p, q). X se nazývá kauzální ARMA proces, jestliže existuje = {j} j=0, P j=0|j| < (tj. 37 1) tak, že Xt = X j=0 jZt-j (zkráceně Xt = (B)Zt), t Z. (4.7a) X se nazývá invertibilní ARMA proces, jestliže existuje = {j} j=0, P j=0|j| < (tj. 1) tak, že X j=0 jXt-j = Zt (zkráceně (B)Xt = Zt), t Z. (4.7b) Poznámka 4.7. Tedy kauzalita znamená, že ARMA proces X je kauzálním lineárním procesem, který lze generovat bílým šumem {Zt}, neboli X MA(). Podobně invertibilita znamená, že bílý šum {Zt} je kauzálním lineárním procesem, který je generován ARMA procesem X, což je ekvivalentní s X AR(). Věta 4.8. Nechť {Xt} ARMA(p, q) : (B)Xt = (B)Zt takový, že polynomy (z) a (z) nemají společné kořeny. Pak platí (i) {Xt} je kauzální právě když (z) = 0 pro všechna z C, |z| 1 (všechny kořeny (z) musí být vně uzavřeného jednotkového kruhu). V takovém případě j jsou určeny Taylorovým rozvojem racionální funkce lomené (z) := (z) (z) = P j=0 jzj , přičemž posloupnost := {j} 1 a je určena jednoznačně. (ii) {Xt} je invertibilní právě když (z) = 0 pro všechna z C, |z| 1 (všechny kořeny (z) musí být vně uzavřeného jednotkového kruhu). V takovém případě j jsou určeny Taylorovým rozvojem racionální funkce lomené (z) := (z) (z) = P j=0 jzj , přičemž posloupnost := {j} 1 a je určena jednoznačně. 38 Věta 4.9 (Autokovarianční funkce MA procesu). {Xt} MA(q), q je stacionární náhodný proces s nulovou střední hodnotou X = 0 a autokovarianční funkcí X (h) = 2 q X k=0 h+kk pro h 0. (4.8a) Speciálně pro q < je X (h) = ( 2 Pq-h k=0 h+kk pro 0 h q 0 pro h > q (4.8b) a zejména X (q) = 2 q = 0, neboť 0 = 1. Důkaz. Vztahy jsou důsledkem 2.24. Důsledek 4.10. 2 X = X (0) = 2 q X k=0 |k|2 . 2 X = 2 (1 + |1|2 + |2|2 + + +|q|2 ) pro q < . (4.9) X (h) = Pq-h k=0 h+kk Pq k=0|k|2 pro h 0. (4.10) Věta 4.11 (Parciální autokorelační funkce kauzálního AR procesu). Nechť {Xt} AR(p), p < je kauzální AR proces. Pak {Xt} je stacionární s nulovou střední hodnotou X = 0 a parciální autokorelační funkcí X , pro niž X (p) = p = 0 a X (h) = 0 pro h > p. Na obrázcích 2, 3 a 4 jsou typické průběhy odhadnuté autokorelační a parciální autokorelační funkce simulovaných procesů po řadě AR(2), MA(2) a ARMA(2, 2). Čárkovaný pás udává interval spolehlivosti 95% pro nulovou hodnotu. Vidíme, že v případě procesu AR(2) na obr. 2, resp. MA(2) na obr. 3 skutečně X (h) 0, resp. X (h) 0 pro h > 2 v souladu s větami 4.11, resp. 4.9. V ostatních případech obálka X (h), resp. X (h) (v případě ARMA(2, 2) na obr. 4 dokonce 39 obou) vykazuje exponenciálně klesající, případně navíc i oscilatorický, průběh. Dá se totiž ukázat, že X i X jsou v těchto případech lineární kombinací klesajících geometrických posloupností a kosinusovek s geometricky klesající amplitudou. 40 Pozorování procesu, n = 500 0 50 100 150 200 250 300 350 400 450 500 -8 -6 -4 -2 0 2 4 6 8 10 12 Data with the mean estimate Autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Autocorrelation function Parciální autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Partial autocorrelation function Obrázek 2. AR(2) : = [0.5, 0.2], 2 = 2.25. 41 Pozorování procesu, n = 500 0 50 100 150 200 250 300 350 400 450 500 -8 -6 -4 -2 0 2 4 6 8 Data with the mean estimate Autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Autocorrelation function Parciální autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Partial autocorrelation function Obrázek 3. MA(2) : = [-0.5, -0.2], 2 = 2.25. 42 Pozorování procesu, n = 500 0 50 100 150 200 250 300 350 400 450 500 -10 -8 -6 -4 -2 0 2 4 6 8 10 Data with the mean estimate Autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Autocorrelation function Parciální autokorelační funkce 0 5 10 15 20 25 30 35 40 45 50 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Partial autocorrelation function Obrázek 4. ARMA(2, 2) : = [0.5, 0.2], = [-0.6, 0.3], 2 = 2.25. 43 4.3.2. (S)ARIMA modely pro kovariančně stacionární časové řady. Do popisu těchto modelů vstupuje operátor diferencování, který lze také ekvivalentně vyjádřit pomocí operátoru zpětného posuvu B popsaném v poznámce 2.27: Diferencování časové řady na vzdálenost s N: Yt := sXt := Xt -Xt-s = Xt -Bs Xt = (1-Bs )Xt, atd. rekurentně pro libovolný řád diferencování d N: Yt := d sXt := s(d-1 s Xt) = (1 - Bs )d Xt. Pro s > 1 se d s = (1 - Bs )d také nazývá operátorem sezónního diferencování řádu d, neboť s hraje roli délky sezónní periody. Pro s = 1 značíme d := d 1 = (1 - B)d a operátor nazýváme operátorem nesezónního diferencování řádu d. Konstrukce modelů: {Xt} SARIMA(p, d, q, P, D, Q, s), jestliže pro Vt := D s Xt platí: Vt = e1Vt-s + e2Vt-2s + + eP Vt-P s | {z } Sezónní autoregresní část +Et+ + e1Et-s + e2Et-2s + + eQEt-Qs | {z } Sezónní část klouzavých součtů , (4.11a) kde {d Et} ARMA(p, q) s parametry jako v (4.11b), přičemž P, Q, d, D 0 a s > 1 jsou opět celá čísla a eP = 0, eQ = 0. {Xt} ARIMA(p, d, q) = SARIMA(p, d, q, 0, 0, 0, 1) neboli Xt = Vt = Et, tj. pro Wt := d Xt ARMA(p, q) platí: Wt = 1Wt-1 + 2Wt-2 + + pWt-p + Zt+ + 1Zt-1 + 2Zt-2 + + qZt-q. (4.11b) Ve všech výše uvedených modelech navíc předpokládáme, že popisují kauzální časovou řadu, neboli současně platí {Xt} MA(), 44 přičemž posloupnost MA-koeficientů je absolutně sumovatelná a příslušná řada konverguje dle kvadratického středu. Použití a interpretace modelu ARIMA: Jestliže aplikujeme operátor na (alespoň lokálně) polynomiální průběh, dojde ke snížení stupně polynomu (nejméně) o 1 (snadno se ověří). Tedy d sníží stupeň nejméně o d. Model ARIMA(p, d, q) se tedy hodí na stacinárně kovarianční časovou řadu X := {X}t, kde stacionarita je porušena jen ve střední hodnotě X (t), která není konstantní a má po částech polynomiální charakter stupně nejvýše d. Nesezónním diferencováním řádu d dosáhneme buď konstantní nebo dokonce nulové střední hodnoty u diferencované řady. Je-li takto získaná řada stacionární (v což obvykle doufáme), pak lze použít model ARIMA. Použití a interpretace modelu SARIMA: Zde předpokládáme, že všechny částečné mezisezónní řady {Xk+s} pro k = 0, 1, . . . , s - 1 se chovají jako ARIMA(P, D, Q) s týmiž parametry ei a ej (tj. nezávisí na k). Reziduální řada {Et} získaná z tohoto modelu je ovšem dekorelována jen mezisezónně a zatížena zbytkovým kolísáním ve střední hodnotě E(t) a tudíž předpokládáme pro ni opět model ARIMA, avšak s jinými řády p, d, q a jinými parametry i a j. Obvykle stačí volit D = 1 (odstraní mezisezónně konstantní sezónní složku). Výjimečně volíme D = 2, jestliže je tendence k mezisezónnímu lineárnímu nárůstu nebo poklesu. Pozorované mezisezónní řady bývají totiž většinou velmi krátké a tudíž takové chování je z nich obtížně ověřitelné. 45 5. Demonstrační ukázky pomocí knihovny TSA-M V práci [Ves00] lze nalézt podrobný popis knihovny funkcí TSA-M vytvořené autorem tohoto příspěvku jako toolbox pro numerické výpočetní prostředí MATLAB. Některé možnosti knihovny jsou zde ilustrovány na reálném datovém souboru zachycujícím měsíční počty úmrtí při automobilových nehodách v USA v letech 1973­1978 převzatém z [BD93, Example 1.1.6] (viz Obr. 1.6 z odst. 1.1). V ukázkách je modelována predikce o jednu sezónu (12 měsíců) jednak pomocí klasického lineárního regresního modelu a jednak užitím modelu SARIMA(0, 1, 1, 0, 1, 1, 12), který se po podrobnější analýze ukázal pro zvolená data jako optimální. Ukázky současně názorně demonstrují posloupnost systematických kroků nutných při analýze časové řady: Předzpracování Volba modelu Identifikace modelu Předběžný odhad jeho parametrů Finální odhad parametrů Verifikace mo- delu. Literatura [And76] Jiří Anděl, Statistická analýza časových řad, SNTL, Praha, 1976. [BD93] Peter J. Brockwell and Richard A. Davis, Time series: Theory and methods, 2-nd ed., Springer-Verlag, New York, 1991 (corrected 2-nd printing 1993). [BD94] , ITSM for Windows. A user's guide to time series modelling and forecasting, Springer-Verlag, New York, 1994, 2 diskettes included. [BD02] , Introduction to time series and forecasting, 2-nd ed., Springer-Verlag, New York, 2002. [Cip86] Tomáš Cipra, Analýza časových řad s aplikacemi v ekonomii, SNTL, Praha, 1986. [Ham94] James D. Hamilton, Time series analysis, Princeton University Press, Princeton, NJ 08540, 1994. [Pri89] M. B. Priestley, Spectral analysis and time series, vol. I­II, Academic Press, London, 1989. [Ves00] V. Veselý, Knihovna program°u TSA-M pro analýzu časových řad, Zborník zo XIV. letnej školy biometriky: Biometrické metódy a modely v pôdohospodárskej vede, výskume a výuke, Račkova dolina, 2. - 6. októbra 2000 (P. Fľak, ed.), VUZVN (Research Inst. of Animal Production Nitra, Slovakia), ASAPV (Agency of Slovak Academy for Agricultural Science), 2000, pp. 239­248. 46 Doc. RNDr. Vítězslav Veselý, CSc., Katedra aplikované matematiky a informatiky, Ekonomicko-správní fakulta MU v Brně, Lipová 41a, 602 00 Brno tel.: +420 5 49498330, fax: +420 5 49491720 E-mail address: vesely@econ.muni.cz, URL: http://www.math.muni.cz/~vesely 47