ÚVOD DO ČASOVÝCH ŘAD VÍTĚZSLAV VESELÝ 1. ÚVOD Časová řada: = soubor pozorování náhodných veličin {xt,t £ r}, kde t je zpravidla čas a T C IR. 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, vit) = acos(vt + 0), tj. xt = - cos(vt + 0), kde a a 0 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 UPLATNENÍ 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ší. 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 Date: 12.1eden 2004. BD Example 1.1.1 BD Example 1.1.2 90 100 Obr. 1.1 Obr. 1.2 BD Example 1.1.3. BD Example 1.1.6. Strikes in the USA, 1951 - 1980 Obr. 1.3 Obr. 1.4 BD Example 1.1.5. BD Example 1.1.6. Tlie Wolf er Lspot numbers, 177 0 Monthly accidental deaths in the USA, 1973 Obr. 1.5 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ů. 2. ZÁKLADNÍ POJMY 2.1. DEFINICE ČASOVÉ ŘADY. Definice 2.1. Náhodný (stochastický) proces X je neprázdný systém (T ^ 0) náhodných veličin definovaných na stejném pravděpodobnostním prostoru (0,^4, P). Píšeme X := {Xt 1t G T} nebo stručně {^t}. Speciální případy: TC1... proces se spojitým časem nebo náhodná funkce. rcZ... 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 u G 0 dostáváme funkci x : T —>■ IR jako výsledek náhodného experimentu: xit) := Xt(u). 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 7 := {t\t = [ti, t2,. . ., tn] G Tn} ti ^ tj, pro i / j, n £ N}. Pro každé t G 7 libovolné délky n G N nechť Ft{x) značí sdruženou distribuční funkci marginálního náhodného vektoru Xt vybraného ze stochastického procesu X = {Xt}teT v časových okamžicích ti,t2,. . . ,tn. Systém {Ft}te7 úplně popisuje stochastické chování 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}te7 z definice 2.6 se nazývá konzistentní protože má pro každé x = [x±, x2,. . ., xn]T G IRn a n G N následující dvě vlastnosti konzistence: (i) Ft (xp) = Ft{x) pro každou permutaci p indexů {l,2,...,n}. Zde značí tp, resp. Xp vektory t, resp. x se složkami permutovanými podle permutace p. (ii) lim^oo Ft(x) = Ft(xi,..., Xi_i, oo, xi+1, ...,xn)=: = : Ft(i)(x(i)) pro každé i G {l,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ý marginální náhodný vektor Xt (t G 7) je normálně rozložený neboli když každá distribuční funkce Ft konzistentního systému tohoto procesu popisuje sdružené normální rozložení. Definice 2.9. Jestliže X = {Xt} je časová řada, kde všechny náhodné veličiny Xt, t G T, jsou vzájemně stochasticky nezávislé a stejně rozložené se střední hodnotou /i a rozptylem R, kde "f x (r, s) := cov(Xr}Xs)} pokud kovariance existují pro všechna r,s£ľ. (3) rozptyl X: a\ : T -^ R+, kde o~x(t) := vaľ(Xt) = cov(Xt}Xt) = "fX(t,t), pokud rozptyly existují pro všechna t G T. (4) autokorelační funkce X: px '■ T x T —> [—1,1], kde / x í / tXÍCf í \ pro V/^x{r,r)y/-fx{s,s) ^ 0 pX(r,s) := < V^x(r,r)^lx(s,s) [O jinak, pokud korelace existují pro všechna r, s G T. (5) vzájemná (křížová) kovarianční funkce X a Y: ^fXY '■ T x T —y R, kde 7xy(r, s) := cov(Xr, Y,), pokud kovariance existují pro všechna r, s G T. (6) vzájemná (křížová) korelační funkce X a Y: PXY :T x T -► [-1,1], kde / x í / ]xÝf f , Pro V^x{r,r)y/-fY{s,s) ŕ 0 pXY(r,s) := < V^(r'r)V^(s>s) [O jinak, pokud korelace existují pro všechna r, s G T. 2.4. STRIKTNÍ A SLABÁ STACIONARITA. Definice 2.11. Časová řada X := {Xt 11 G Z} se nazývá striktně stacionární, jestliže každá distribuční funkce jejího konzistentního systému {Ft}teg-, nezávisí na posunutí: Ft(-) = Ft+h(-) pro každé t G 7 a /i G Z. Definice 2.12. Časová řada X := {Xt 11 G Z} se nazývá (slabě) stationární, jestliže jsou splněny následující podmínky: (1) X má konečné druhé momenty: &x(t) < oo pro každé t G Z. (2) "fX(r, 5) = 7x(r + /i, s + /i) každé r,s,li£ Z. (3) px(') = f*i je konstantní funkce. Jestliže jsou splněny pouze podmínky (1) a (2), pak X se nazývá kovariančně stacionární. 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: (Tx(-) = o\. (2) Jestliže platí (3), pak z 7X(r, s) = EXrXs - (EXr)(EXs) = EXrXs - p\ plyne, že (2) je ekvivalentní (a může tak být nahrazena) podmínkou: EXrXs = EXr+hXs+h pro každé r,s,h G 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') 7x(y, 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é: 7x(h) :=jx(t + h,t) px{h):=px{t + h,t) =----------------= ——- (2.1) vx°x 7x(U) u x 7x(M) = 7x(0), kde í,/i£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 /i a rozptylem x(t) = /i a I cr2 pro r = s I 0 jinak Píšeme X ~ WN(ß, a2 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(/i,a2) => X ~ WN(/i,a2) => X je stacionární. Poznamenejme také, že žádná z opačných implikací neplatí. 2.5. LINEÁRNÍ TRANSFORMACE ČASOVÝCH RAD. Definice 2.20. Řekneme, že časová řada X := {Xt} má ohraničené druhé momenty, jestliže existuje konstanta C G IR s vlastností E|Xt|2 < C < oo pro každé t. Věta 2.21. Časová řada X := {Xt} má ohraničené druhé momenty právě když existují reálné konstanty fix < oo a a x < oo takové, že \/J,x(t)\ < fíj a \crx(t)\ < o 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 $ := {$fc} G í\ (tj. J^J^fcl < ooj a U := {Ut} je časová řada s ohraničenými druhými momenty, pak nekonečná řada Xt := ^2k <&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. tj. částečné součty konvergují dle kvadratického středu: E\Xt — ^2k=_N <5>kUt-k\2 —* 0. Součet v tomto smyslu označujeme symbolem =. V dalším textu bude sumace řad vždy uvažována v tomto smyslu, i když symbol 2 nad rovnítkem bude vynechán. Poznámka 2.23. Lineární transformace s koeficienty $ := {$fc}, 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 $ £ í\ a U = {Ut} je stacionární se střední hodnotou pu a autokovarianční funkcí ^jj, pak platí 2 (1) Xt = ^2k ^kUt-k konverguje pro každé t (2) X := {Xt} je stacionární se střední hodnotou a autokovarianční funkcí: Px = Pu Y^ $j (2-2a) i lx(h) = YJYJ ®&klu(h -j + k), h> 0. (2.2b) Důsledek 2.25. a x 7x(o)=x;e^$^(a;-^ (2-3a) <4<^(I>'-l2) (2-3b) i 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 $/; = 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 Ut,t > t generující řady (jsou podmíněny jen hodnotami Ut,t < t). V takovém případě částečné součty při sumaci dle kvadratického středu nabudou tvaru ^2k=0 ^kUt-k a samotnou řadu tak můžeme psát ve tvaru Xt = Y<ľ=o ®kUt-k- Poznámka 2.27 (Operátor zpětného posuvu). Značíme pro HN: BUt := Ut-u atd. rekurentně BkUt = B(Bk-lUt) := Ut-k B~lUt := Ut+1, atd. rekurentně B~kUt := B^iB'^-^Ut) := Ut+k B°Ut := Ut. Pak formálně místo Xt = £~ 0 Wt-* píšeme stručně Xt = $(B)Ut, kde $(B):=ZT=o®kB\nehoť Y,ľ=o ^kUt-k = Sfclo $kBkUt = (YlkLo $kBk)Ut je operátor získaný formálním dosazením operátoru B za proměnnou v řadě (polynomu) n o£ wfck. d_ fc e noise WIST C O , 3. ) :~*«°A*°<^«^^^ 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 < -|. (2) Odhad ^f(h) je vychýlený (ß/y(h) ^ "f(h)), protože součet pozorování je dělen n a nikoliv počtem stupňů volnosti n — l — h. Poznamenejme, že věta 2.29 neplatí pro nevychýlený odhad (matice Yn ztratí očekávanou vlastnost pozitivní semide-finitnosti). V každém případě je ale odhad j(h) asymptoticky nevychýlený v tom smyslu, že F/y(h) —y "f (h) pro n —y oo. Navíc je také konzistentní podle kvadratického středu v tom smyslu, že E|íy(/i) — 7(/^)|2 —> 0 pro n —y oo, kde konvergence je dokonce rychlejší než v případě nevychýleného odhadu. (3) Z algebraického hlediska ^f(h) lze psát ve tvaru skalárního součinu 7(/i) = ^(xo, Xh), kde Xh '■= [0,. . ., 0, X\ — ju,. . . , xn — /i, 0,. . . , Oj. Vektor x0 tak repre- h n—l—h zentuje 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ě ||íKo||2 = ll^/ill2 = XľLil^j ~~ A* P'• Ze Schwarzovy nerovnosti dostáváme \{x0,xh)\ < ||a?o||2, což vede na \j(h)\ < ^\\x0\\2 = ±(x0,x0) = 7(0). Odtud vidíme, že odhad autokorelační funkce zachovává její přirozenou vlastnost IW|<1. (4) Vzhledem ke (3) je možno odhad p(h) interpretovat geometricky jako kosinus ú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á orto-gonalitu (ú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 ^f(h) při h —y 00. Periodická složka se projeví oscilatorickým chováním ^f(h) s dominantní periodou této složky, případně směsí takových složek. 3. Predikce v časových řadách 3.1. PRINCIP ORTOGONÁLNÍ PROJEKCE. Definice 3.1 (Prostor L2(íl,A,P)). L2 := L2(íl}A} P) definujeme jako množinu všech (reálných nebo komplexních) náhodných veličin nad týmž pravděpodobnostním prostorem (fž,./l, P), které mají konečné druhé momenty (resp. rozptyly - viz dále 3.5), tj. L2(tt,A, P) := {X I X náhodná veličina nad (0,^1, P), E|X|2 < 00}. Poznamenejme, že do tohoto prostoru zahrnujeme také všechny konstanty z C, které považujeme za náhodné veličiny s nulovým rozptylem. Věta 3.2 ([BD93, s.46-48, §2.10*]). L2(íl,A, P) je Hubertův prostor2 se skalárním součinem (X,Y) = EXY a normou \\X\\2 = \/(X, X) = yE|X|2_ Tedy konvergence indukovaná touto normou je právě konvergence podle kvadratického středu. Důsledek 3.3 (Schwarzova nerovnost). \EXY\ < ||X||2||V||2 = VE|X|VE|1T, X,Y e L2(Ü,A,P) Důsledek 3.4. X G L2(íl,A,P) =*► EX existuje. Důkaz. \EX\ = \E(1.X)\ < VEN2 VEI^P = VEI^P < 00. 1 Důsledek 3.5. X,Y eL2(íl,A,P) =>• X-EX, Y -EV G L2(íl,A,P) a (X -EX, V - EV) = E(X - EX)(Y - EV) = cov(X, V) existuje a splňuje Schwarzovu nerovnost \cov(X,Y)\ < y/E\X - EX\2\fE\Y - EV|2 = oxoY. D Hubertův prostor představuje zobecnění euklidovského prostoru připouštějící i nekonečnou dimenzi v případě, že prvky prostoru jsou funkce, náhodné veličiny apod. Důsledek 3.6. cov(x,y) ( covI,ľ / n p(X,Y) = \ -^ Pr° ax(J^0 [ O pro oxoy = O je tzv. korelační koeficient náhodných veličin X &Y, pro nějž platí \p(X,Y)\ < 1 a speciálně — 1 < p(X} Y) < 1 v případě reálných náhodných veličin X aY. 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 G H, pak následující tvrzení jsou ekvivalentní: (i) existuje jediný prvek x G L takový, že \\x — x\\ = inf llx — y\\ (*) yeL a (ii) nějaký prvek má minimální vlastnost (*) právě když x — x _l_ L, tj. x — x _l_ y (neboli [x — x,y) = 0) pro každé y G L. Píšeme x = Plx, kde Pl je tzv. operátor ortogonální projekce na podprostor L. Definice 3.8. Nechť X0 := 1, X\,. . ., Xn G L2 a L := C(X0} Xí}. . ., Xn) je lineárni podprostor v L2 generovaný náhodnými veličinami X0}Xi}. . . }Xn. Je-li X G L2 nějaká další náhodná veličina, pak X = P^X se nazývá nejlepší3 lineární predikce náhodné veličiny X pomocí náhodných veličin X0} Xí}. . ., Xn. Poznámka 3.9. Operátor Pl je korektně definován, neboť L má konečnou dimenzi (dimL < n + 1) a je tak automaticky uzavřený. Protože X G C(X0} Xí}. . . , Xn), existuje vektor ß0 := [ß0, ßi,. . . , ßn]T G IRn+1 takový, že X = ß0 X0 +ß1X1 + --- + ßnXn. (3.1a) i I když X je jediný, nemusí být ß0 obecně jednoznačně určen. Je tomu tak jen když X0} X\,. . ., 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 rovnicj, kde i = 0,1,. . . , n: (X0,Xi)ß0 + {XuX^ßy + ■■■ + (Xn,Xi)ßn = (X,Xi) (3.1b) neboli dle věty 3.2 E(X0Xi)ß0 + E(X1Xí)/3i + • • • + E(Xn,Xi)ßn = E(XXt). (3.1c) Důkaz. Dle 3.7: X = PLX ^ X - X ± Y pro každé YeL <£> X - X ± X% pro každé i = 0,1,..., n <£> (X - X,Xi) = 0 pro každé i = 0,1,..., n <£> (X, Xt) = (X, Xt) pro každé i = 0,1,. . . , n <^ (^?=o ßjX17 Xi) = (X} Xi) pro každé i = 0,1,. . . , n <^ Sj=o ßAXji Xi) = (X} Xi) pro každé i = 0,1,. . . , n, což je právě (3.1b). D Příklad 3.11 (n = 0). L = jC(Xo) = jC(1) = IR je lineární podprostor všech (reálných) konstant v L2. Pak (3.1b) se redukuje jen na jednu rovnici: E(X0X0)ß0 = E(X X0 ). i i rozumí se ve smyslu minimální střední kvadratické chyby. Odtud máme ß0 = EX a tedy X = ßoX0 = 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. = 0; pak X = PLX 1,. . . , n; P C(X! (3.1d) (3.1e) Důsledek 3.12. Jestliže EX = EX: = ... EX ßo = 0 a (3.1c) přejde v systém n rovnic, kde i E^X^i + ■■■ + E(Xn,Xi)ßn = E(XXi). Položíme-li X := [Xi,. . . ,Xn]T, můžeme jej přepsat do maticového tvaru (varX)/3 = cov(X, X)T, kde ß :=[/3i,..., ßn]T G Rn. Důkaz. První rovnice (i = 0) v (3.1c) bude tvaru: E(l.l) ßo + E(XlI) ß, + • • • + E(Xn.l) ßn = E(X.l). 10 0 0 Odtud dostáváme ß0 = 0, takže ve zbývajících rovnicích pro i = 1, E(X0Xi)/30 = 0, což je právě systém normálních rovnic pro Pc(xu...,xn)X. Důsledek 3.13 (Nejlepší lineární predikce pro stacionární časové řady). NechťX := {Xt} je stacionární časová řada s nulovou střední hodnotou fix = 0 a auto- kovarianční funkcí ~fx ■ Pro dané k £ N nechť X\\ značí nejlepší (k-krokovou) lineární n bude D í-i? X t+i-n ve predikci náhodné veličiny Xt+k pomocí n předchozích veličin Xt,X- tvaru Xí+fc = Ej=i ®n jXt+i-j ■ Pak vektor $^ ' := [<&nl ,. . . , $L]T nalezneme řešením tzv. Yule-Walkerova systému lineárních rovnic; 7(0) 7(1) 7(1) 7(0) 7(n 7(n 1) 2) Pj{n — 1) 7(n — 2) 7(0) ní $ (k) n2 ' l(k) ■ 7(fc + l) _"f(k + n — 1), (3.2) * « „(*) (k) Důkaz. Plyne z 3.12 po záměnách: X -^» Xt+fc, Xj -^» Xí+i_j, /3j -^» $^ a tedy pak varX = [cov^-X,)] ^ [cov(Xí+1_8,Xí+1_,)] = [y(t + l-i-(t + l-j))] = [y(j-i)] = Tr. a podobně 'w. ' D cov(X,X)^ ^ [E(Xí+fc,Xí+1_,)] = ••• = h(t+k-(t+l-j))] = h(k+j-l)] Věta 3.14 (Střední kvadratická chyba). Střední kvadratická chyba vn := ElXt+fc — Xí+fc|2 nejlepší k-krokové lineární predikce se spočte dle vztahu: v. Důk az. v. (k) nk) = 7(0) - Y, *nh(k +J~1)= 7(0) - Tnlík). 3 = 1 T7IV v(fc) |2 — \\Y v(k) i|2 _ lv hj\At+k - At+k\ - \\Xt+k - At+k\\ - {At+k (3.3) {At+k - At+k,At+k/ \At+k — A+ , ,„, A xt+k - t+ki y^~t+kl {Xt+k,X-, t+ki ^t+k/ y{k) Y ^t+k,^t+k ,(k) y(k)\ I Y^Ki Y ^ \At+k,At+ki IX M o >(fc)/ t+k\\~ - \Ej = l ®nj Xt+1_^Xt+k) - ||Xí+fc|| - Yľj = l ®nj (Xt+l-j,Xt+k) 3=í -nlnxt+l-3Xi+k) = 7(o) - e;=i *S}-S)7(i-i-fc) = 7(o)-E;=i*S)- E(^2+fc) - ZU WhiA*) = 7(0) - E;=i ^7(í + 1 - 3 - (t + k)) 7(0) - £?=i $g}7(l - J - A:) = 7(0) - EL «^í* + J " !)■ D Definice 3.15. Nechť X} Y, X\,. . . , Xn G L2 a X a Y jsou nejlepší lineární predikce náhodných veličin X a Y pomocí X\,..., Xn. Pak /?(X, K | X\,..., Xn) := p(X — X} Y — Y) se nazývá parciální korelační koeficient náhodných veličin XaY při daných Xi,... ,Xn. Interpretace: parciální korelace mezi XaY při daných X\,. . . , 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 X\,. . . , Xn. Definice 3.16. Nechť X = {Xt\t G Z} je stacionární časová řada s autokorelační funkcí px(-)- Pak parciální autokorelační funkci (pacf) ax(-) časové řady X definujeme následovně: ax(0)=px(0) = l, ax(l)=Px(l) = p(Xt+uXt) ax{h) = p(Xt+h,Xt \Xt+1,.. .,Xí+/l_i) pro h>2. Věta 3.17 ([BD93, §5.2]). Jestliže X = {Xt \t G Z} je stacionární časová řada se střední hodnotou px = 0 a parciální autokorelační funkcí ctx(-), pak ctx(n) = §n,n Vro n ^ 1; kde <&„ = [$n,i5 • • • ? ^n,™]7 je řešením problému nejlepší lineární predikce. ctx(h) lze počítat rekurzivně užitím mezivýsledných vztahů (3.4b) dále uvedeného Dur-bin-Levinsonova algoritmu 3.19. Jiný postup počítá $„,1 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 Tn := [^x{i — Í)]™?=i regulární, pak kde r* := [r(:, 1),. . ., T(:, n - l),fn] «7« = [lx(l), ■ ■ ■ , lx{n)]T. 3.2. DURBIN-LEVINSONÜV A INOVAČNÍ ALGORITMUS. Věta 3.19 (Durbin-Levinsonův algoritmus pro rekur. výpočet <&„ [BD93, Prop.5.2.1]). Nechť X = {Xt \t G Z} je stacionární časová řada s nulovou střední hodnotou px = 0 a autokovarianční funkcí ^fx(h) —>■ 0 pro h —y oo. Jestliže Xn+\ = $rJilXrj + - ' mJr*&n,nXi je nejlepší lineární predikce, pak pro koeficienty $nj a střední kvadratickou chybu vn = E|Xn+i — Xn+i\2 platí následující rekurentní vztahy. Počáteční krok pro n = 0: i>o = 7*(0). (3.4a) Rekurentní krok pro n > 0: =0 pro n=l $* [7x1 n) n-l -i,jlx\n j)]'- n-li $n,2 $ n.n — 1 -i(l-|$ $n-l,l $n-l,2 n,n I y; $ n —l,n —1 $, ^n-l,n-l &n-l,n-2 $ ra-1,1 pro n > 1. (3.4b) (3.4c) (3.4d) 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 fix = 0. Tento algoritmus není nic jiného než souřadnicový přepis Gram-Schmidtova ortogonalizačního procesu. Věta 3.20 (Inovační algoritmus pro rekurentní výpočet <&„ [BD93, Prop.5.2.2]). Nechť X = {Xt \t £ Z} je časová řada s nulovou střední hodnotou fix = 0 s maticí 3ii,j = l> "-«j •" smíšených momentů [/€;, nokrokové nejlepší lineární predikce Xn+Í := X^1} n > 0; a jejich střední kvadratické chyby vn := vn spočteme rekurentně pomocí následujících vztahů: Fj(XiXj), která je regulární pro každé n. Pak jed- X, n+í 0 pro n = 0 YJj = l ®n,j(Xn+1-j - -^rx+l-j) pTO U > l, užitím Vq k(1A) 0 n.n —k -lU(n + l,k + l) k-í / j ^k,k—j^n,n—j^j^ j=0 Qk.k-iQr. vn = K{n + 1, n + 1 postupně v pořadí: v0; 0i,i, v\\ 02,2, 02,i, v2] ■ ■ ■ n-l j=o n,n-jVj- (3.5a) (3.5b) (3.5c) (3.5d) 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 predzá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' 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í modelu. 4.1.2. Stabilizace rozptylu. Většina běžně užívaných modelů předpokládá alespoň kovariančně stacionární chocání analyzované řady X := {^t}, tj. zejména konstantní rozptyl o\ (tzv. homoskedasti-cita). 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ě: ax(t) = (Jo^x(t)@- V praxi se nejčastěji vyskytují 2 případy: • 0 = 0: řada je již homoskedastická a není třeba ji transformovat; • 0 = 1: směrodatná odchylka řady závisí lineárně na střední hodnotě. Existují postupy, které umožňují odhadnout parametr 0 z pozorované řady. He-teroskedasticitu pak odstraníme jednou z níže uvedených transformací, kde zvolíme parametr A = 1 — 0. Box-Coxova transformace pro řadu Xt > 0: Yt ■= Í^Ť1 pro A ^ 0 I ln(Xt) pro A = 0 (lineární nárůst při 0 = 1) Mocninná transformace pro řadu Xt > 0: Yt .= { X? pro A ^ 0 ln(Xt) pro A = 0 (lineární nárůst při 0 = 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 A —y 0, takže je nutno se vyhýbat malým nenulovým hodnotám parametru A. (3) Postup odphadování parametru 0, resp A 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 xit) do její Fourierovy řady: oo oo 0 , ^-~\ i2Trkt CLq ^-~\ /Z7T/CÍ \ a II x(t) = 2_^ cke T = — + 2_^Akcos{—------,k = l,...,m. (4.3) Velká hodnota Ik indikuje velkou energii k-té harmonické komponenty, tj. silné zastoupení složky xkit) = ckeiUJkt + c_ke~iUJkt = Ak cos(ukt — ^/, ;2nhk I\ — 1 h = 2 2^ l(h)e 1k pro k = 1,2,..., —— \h\ 0 pro j = 1, 2,. . . , m, je 0 < W < 1. Nechť Vt = Pt + Pt, £t ~ WW(0, a2) gaussovský. Platí-li nulová hypotéza H0 : Xt = Et, pak testová statistika W má na intervalu [0,1] rozdělení, pro nějž platí p(w > x) = J2 i-1)1'1 (7) í1 - ^r~\ o < x < 1, l-jaľ>0 přičemž PiW > x) ss ra(l — x)m_1 je dobrá aproximace pro m < 50. Nechť <7f(1 — a) značí (1 — a)-kvantil rozdělení statistiky W (tabelováno), tj. P(W < (/f(1 — o)) = 1 — a. Pak H0 zamítáme s rizikem a, pokud W > (/f(1 — «)• Je-li v takovém případě Ik0 = maxjt=i,...,m Ik, Pak A;0-tou harmonickou komponentu považujeme za významnou. Další významnou harmonickou komponentu zjistíme z pe-riodogramu délky m — 1, kde vynecháme 7fc0, 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 m gF(l — a))+, 0 < A < 1 (doporučená hodnota je A = 0, 6). fc=i H0 zamítáme s rizikem a, jestliže Tx > t\(l — a), kde t\(l — a) je (1 — a)-kvantil rozdělení T x (tabelováno). Poznámka 4.3 (Praktická doporučení). (1) Vzhledem k povaze nulové hypotézy ve výše uvedených testech by v testovaných datech neměl být obsažen trend. Doporučuje se proto 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élkou 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. Pt Aditivní model (AM): Xt = Trt + 1Šž?+CÍ+Et, Dt Pt Multiplikativní model (MM): Xt = Trt.'Sz^Cl.Et, Dt kde značí Trt ■ ■ ■ dlouhodobý trend Szt ■ ■ . sezónní komponenta (perioda je předem známa) Ct . . . cyklická komponenta (perioda není předem známa) Pt ■ ■ ■ celková periodická komponenta Dt = Hx{t) . . . deterministická komponenta (střední hodnota Et ■ ■ ■ stacionární náhodná komponenta (obvykle bílý šum nebo IID): \íe = 0 pro AM, resp. he = 1 pro MM. AM lze použít jen pro homoskedastická data, neboť er x if) = a(Dt-\-Et) = cr(Et) = oe je konstantní vzhledem ke stacionaritě Et. Naopak MM lze použít jen pro lineárně heteroskedastická data, neboť (Tx(t) = a(Dt.Et) ~-Dtcr(Et) = Hx{t)(TE závisí lineárně na střední hodnotě jix{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 = X: = (ti,... = xt = t )T ti f= t j ■ ■ 7 Xtn pro i ^ J x : = = xt = (xi,.. x )T ) 0. j * ,j Pokud D závisí lineárně na ßi}. . . ,ßp, pak tento systém je systémem p lineárních rovnic pro p neznámých ßi}. . ., ßp a dostáváme klasický lineární regresní model. V případě nelineární regrese použijeme pro minimalizaci N(ß\,. . . , ßp) vhodnou optimalizační metodu, například Optimization Toolbox v prostředí systému MATLAB. Lineární regresní model: Dit] ßi}. . . , ßp) = (fi(t)ßi + • • • + (fp(t)ßp} kde (f i jsou dané. X = Dß + E, E = (Etl,..., Etnf ~ IID(0, a2), kde D = [■ {Dt} ve vhodné konvergenci při n —y oo, tj. odhad je asymptoticky přesný. Běžné vyhlazovací techniky: • Metoda klouzavých průměrů (lineární konvoluční filtr): Q Q Q Yt = ^2 WTXt+T = ^2 WrDt+T + ^2 WrEt+T, Dt:=S({Dt}) kde požadujeme: (1) Co nejmenší zkreslení Dt, tj. Dt ~ St=- wrDt+T- 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 Dt velmi dobrá. Odpověď dává následující věta: Věta 4.4. Pro každý polynom Pt = c0 + c\t + • • • + c^t nejvýše k-ho stupně platí Pt = ^29T=-q wrPt+T právě když váhy splňují následující 2 podmínky ^2wT = l, T = -q y TrwT = 0 pro r = 1,. . . , k. (4.5) T = -q (2) Co největší redukci rozptylu náhodné složky Et} tj. 0 ~ St=- wrEt+T-Mírou pro tento účel je tzv. redukční intenzita klouzavých průměrů definovaná j ako podíl rozptylů o\ja\ za předpokladu, že Et ~ WN(0, a2). Spočtěme nejprve rozptyl oY: o\ = vaiEt = var(^'__ wTEt+T) = V9 wlvzr(Et+T) = a2 V9 w2 = o-2\\w\\2. / 'i— — q t V l\> ) / ;>->— — q r W W Odtud o-y/o-j, = var Éŕ/var Et = ct2||«;||2/(T2 = ||i«||2. Požadujeme tedy \\w\\2 < 1 co nejmenší. Dá se ukázat, že z tohoto pohledu vyhovují nejlépe obyčejné klouzavé průměry wT = j-tj-, kde II™ IP = Er = -g (2g + l)2 = (2q + í)2 = 2^+T ■ ^ VŠ&k Z&Se neJSOU PřílÍŠ 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 = ^29 _ 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,..., wr] = ^[74, 67,46, 21, 3, —5, —6, —3], kde w_T = wT pro \t\ < 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í j^ = 0.0667 danou obyčejnými klouzavými průměry délky 15. • Exponenciální vyrovnávání: viz [Cip86, kap.7]. • Jádrové vyhlazování: viz [Ves94], [Ves96a]. • Waveletové vyhlazování: viz [Ves96b], [Ves97b], [Ves97a] a případně zobecnění těchto technik pro neortogonální a přeurčené bázové systémy [CDS98], [Ves02], [ZVH04]. 4.3. BOX-JENKINSOVY MODELY. 4.3.1. ARMA modely pro stacionární časové řady. Definice 4.5 (ARMA proces). Stochastický proces X = {Xt \t G Z} se nazývá ARMA procesem řádu p, q (0 < p}q < oo), píšeme X ~ ARMA(p} q), jestliže Xt = $!**_! + $2Xt_2 + • • • + %Xt_p + Zt + Q1Zt_1 + 02Zt_2 + • • • + 0gZt_g, Autoregresní část (AR) Část klouzavých součtů (MA=Moving Average) (4.6a) kde Z := {Zt \t G Z} ~ WN(0,a2) a $p ^ 0, 0g ^ 0, d/0. Přepíšeme-li (4.6a) do ekvivalentního tvaru Xt - ^Xt-! - $2Xt_2----------%Xt_p = Zt + Q1Zt_1 + 02Zt_2 + • • • + 0gZt_g, (4.6b) lze použít též zkrácený zápis dle poznámky 2.27: ${B)Xt = 6{B)Zt kde při $o = 0o = 1 je $(z) = i-$lZ---------pzp, 6(z) = i + elZ +■■■ + eqzq, zeC. Poznámka 4.6. V předchozí definici bez újmy na obecnosti předpokládáme 0O = 1, neboť v opačném případě stačí nahradit bílý šum {Zt} ~ WN(0, a2) šumem {Q0Zt} ~ MV(O,(0oa)2). 4.6c Poznámka 4.7 (Speciální případy). a) Autoregresní proces (AR proces): X ~ AR(p) = ARMA(p, 0) : = {^j}^l0, X^ol^'l ^ °° W- ^ ^ ^0 ^ak, ^e oo Aľt = J^ V'j^í-j (zkráceně Xt = ý(B)Zt), t G Z. (4.7a) j=o X se nazývá invertibilní ARMA proces, jestliže existuje TT = {^j)f=o, EjlolV'jl < oo (tj. vr G £i) tak, že OO ^TTjXt-j = Zt (zkráceně Tľ(B)Xt = Zt), t G Z. (4.7b) j=o Poznámka 4.9. 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(oo). 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í sl~ AR(oo). Věta 4.10. Nechť {Xt} ~ ARMA(p,q) : $(B)Xt = 0{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 G C; \z\ < 1 (všechny kořeny &(z) musí být vně uzavřeného jednotkového kruhu). V takovém případě i/jj jsou určeny Taylorovým rozvojem racionální funkce lomené i^(z) := ^ffr = Y^TLo^jZ3> přičemž posloupnost ip := {ipj} G í\ a je určena jednoznačně. (ii) {Xt} je invertibilní právě když ©(z) ^ O pro všechna z G C; \z\ < 1 (všechny kořeny ©(z) musí být vně uzavřeného jednotkového kruhu). V takovém případě iTj jsou určeny Taylorovým rozvojem racionální funkce lomené tt(z) := ^W- = Sj^o ^jZ1 ' přičemž posloupnost it := {ttj} G i\ a je určena jednoznačně. Věta 4.11 (Autokovarianční funkce M A procesu). {Xt} ~ MA(q), q < oo je stacionární náhodný proces s nulovou střední hodnotou fix = 0 a autokovarianční funkcí lx{h) = 0. (4.8a) fc=o Speciálně pro q < oo je t=o ö/^+fcöfc pro 0 < h < q lx{h) = { ^fc=ü '^ * r ^ ^ ^ (4.8b) n pro h > q a zejména "fx (q) = c2©g 7^ 0; neboť 0O = 1. Důkaz. Vztahy jsou důsledkem 2.24. D Důsledek 4.12. 9 2 4 = 7x(0) = a2]T|O fei fc=o (4.9) <4 = a\l + |0!^ + \Q2\ +■■■ + +|0g|) pro g < oo PxW=E^=°0fen+ffcpro/.>O. (4.10) Věta 4.13 (Parciální autokorelační funkce kauzálního AR procesu). Nechť{Xt} ~ AR(p)} p < oo je kauzální AR proces. Pak {Xt} je stacionární s nulovou střední hodnotou fix = 0 a parciální autokorelační funkcí «j, pro niž a x (p) = $p 7^ 0 a ctx(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ě ctx(h) ss 0, resp. p x (h) ~ 0 pro h > 2 v souladu s větami 4.13, resp. 4.11. V ostatních případech obálka px(h)} resp. a x (h) (v případě ARMA(2,2) na obr. 4 dokonce obou) vykazuje exponenciálně klesající, případně navíc i oscilatorický, průběh. Dá se totiž ukázat, že px i «x jsou v těchto případech lineární kombinací klesajících geometrických posloupností a kosinu-sovek s geometricky klesající amplitudou. 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 G N: Yt := AsXt := Xt - Xt_s = Xt - BsXt = (1 - Bs)Xt, atd. rekurentně pro libovolný t := AJAf-1^) = (1 - Ba)dXt. řád diferencování d G N: K := AdXt := AJA^X) = (1 - Bs)dXt. Pozorování procesu, n = 500 Autokorelační funkce t Ťm_........ „ _ „ m "P "Ť* T «P =P CT, «P *P CT, «=, «P -H=- CT, «P "P < ......Aiilii"11*.. ................................... Parciální autokorelační funkce Obrázek 2. Aiž(2) : $ = [0.5,0.2], a2 = 2.25. Pro s > 1 se Af = (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 Ad := Af = (l-B)d a operátor nazýváme operátorem nesezonniho diferencování řádu d. Konstrukce modelů: {Xt} ~ SARIMA(p, d, q, P, D, Q, s), jestliže pro Vt := Af Xt platí: K = $iK_s + $2Vť_2s + • • • + $pK-ps +Pt + 0iEť_s + 02Et-2s + ■■■ + QnE, t — Qsf Sezónní autoregresní část Sezónní část klouzavých součtů (4.11a) kde {AdEt} ~ 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 $p ^ 0, 0g 7^ 0. {Xt} ~ ARIMA(p, d, q) = SARIMA(p, d, q, 0, 0, 0,1) neboli Xt = Vt = Pt, tj. pro Wt := AdXt ~ ARMA(p, q) platí: Wť = $iWť_i + $2Wť_2 + • • • + $P^Í-P + Zť + 0!^.! + 02Zt_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} ~ AP4(oo), přičemž posloupnost MA-koeficientů je absolutně sumovatelná a příslušná řada konverguje dle kvadratického středu. Pozorování procesu, n = 500 Autokorelační funkce Parciální autokorelační funkce ±-H- ___ «=P ^p <=+=> ■=■ => T* T •+» *+» *¥* ť ^«LJ^^i___________-t _ í-élJ^^ _J>_ Obrázek 3. MA(2) : 0 = [-0.5, -0.2], a2 = 2.25. Použití a interpretace modelu ARIMA: Jestliže aplikujeme operátor A na (alespoň lokálně) polynomiální průběh, dojde ke snížení stupně polynomu (nejméně) o 1 (snadno se ověří). Tedy Ad sníží stupeň nejméně o d. Model ARIMA(p} d} q) se tedy hodí na stacionárně kovarianční časovou řadu X : = {Xt}, kde stacionarita je porušena jen ve střední hodnotě /j,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 {Xfc+TS}T pro k = 0,1,. . . , s — 1 se chovají jako ARIMA(P} D} Q) s týmiž parametry $; a 0j (tj. nezávisí na k). Reziduálni řada {Et} získaná z tohoto modeluje ovšem dekorelována jen mezisezónně a zatížena zbytkovým kolísáním ve střední hodnotě HeÍ^) 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 $; a Qr 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é. Pozorovaní procesu, n = 500 Autokorelační funkce :TÍ ■TfT T°'±Ai^=-""i,''í"p-=T. ,T,T, =.i^.-^J,-í'.>ť?>J* Parciálni autokorelační funkce Obrázek 4. ARMA(2,2) : $ = [0.5,0.2], 0 = [-0.6,0.3], a2 = 2.25. 5. Demonstrační ukázky pomocí knihovny TSA-M V práci [VesOO] 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í MÁTL AB. 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í —y Volba modelu —y Identifikace modelu —y Předběžný odhad jeho parametrů —y Finální odhad parametrů —y Verifikace modelu. Pro odhad parametrů se standardně používají techniky založené na minimalizaci střední kvadratické chyby užitím ortogonální projekce. V případě (S)AR(I)MA modelů se však obvykle užívají pouze v roli předběžného odhadu, z něhož startujeme iterace pro nalezení finálních, zpravidla tzv. maximálně věrohodných odhadů, které maximalizují příslušnou věrohodnostní funkci. Při verifikaci modelu hrají důležitou roli rezidua (rozdíl mezi daty a modelem). Je-li model adekvátní pro data, pak reziduálni řada by neměla vykazovat významné korelace. V případě (S)AR(I)MA modelů by dokonce měla být pozorováním bílého šumu {Zt}. Pro tento účel se používají nejrůznější testy, z nichž nejznámější je tzv. Portmanteau test: podrobněji viz např. [Cip86, kap.10, odst.12.3] a [BD93, §9.4]. Literatura [And76] Jiří Anděl, Statistická analýza časových řad, SNTL, Praha, 1976. [BD93] Peter J. Brockwell and Richard A. Davis, Time senes: 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. [CDS98] S. S. Chen, D. L. Donoho, and M. A. Saunders, Atomic decomposition by basis pursuit, SIAM J. Sei. Comput. 20 (1998), no. 1, 33-61, reprinted in SIAM Review, 43 (2001), no. 1, pp.129-159. [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. [Ves94] V. Veselý, Jádrové odhady — implementace, ilustrativní příklady, aplikace, sborník zimní školy ROBUST'94, Malenovice (J. Antoch and G. Dohnal, eds.), JCMF (Jedn. čes. matematiků a fyziků), leden 1994, pp. 160-171. [Ves96a] _______, Kernel and wavelet smoothing: Basic theory and examples of practical data analysis, Proceedings and abstracts of ESES'96 — Environmental Statistics and Earth Science (Brno) (I. Horová, J. Jurečková, and J. Vosmanský, eds.), Masaryk University of Brno, Czech Rep., August 1996, pp. 163-174. [Ves96b] _______, Waveletová analýza dat, sborník celost, konf. ANALÝZA DAT'96, hotel Technik, Lázně Bohdaneč (K. Kupka, ed.), TriloByte, Ltd. (Pardubice, Czech Rep.), listopad 1996, pp. 1-13. [Ves97a] _______, Wavelety — úvod do teorie a aplikace, sborník celost, konf. ANALÝZA DAT'97, hotel Technik, Lázně Bohdaneč, 4.11-7.11.1997 (K. Kupka, ed.), TriloByte, Ltd. (Pardubice, Czech Rep.), listopad 1997, pp. 83-97. [Ves97b] _______, Wavelety a jejich použití při filtraci dat, sborník letní školy ROBUST'96, Lednice, září 1996 (J. Antoch and G. Dohnal, eds.), JCMF (Jedn. čes. matematiků a fyziků), 1997, pp. 241-272. [VesOO] _______, Knihovna programů TS A-M pro analýzu časových řad, Zborník zo XIV. letnej školy biometriky: Biometrické metody a modely v pôdohospodárskej vede, výskume a výuke, Rackova 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. [Ves02] _______, Hubert-space techniques for spectral representation in terms of overcomplete bases, Proceedings of the summer school DATASTAT'2001, Cihák near Zamberk (I. Horová, ed.), Folia Fac. Sei. Nat. Univ. Masaryk. Brunensis, Mathematica, vol. 11, Dept. of Appl. Math., Masaryk University of Brno, Czech Rep., 2002, pp. 259-273. [ZVH04] J. Zelinka, V. Veselý, and I. Horová, Comparative study of two kernel smoothing techniques, Proceedings of the summer school DATASTAT'2003, Svratka (I. Horová, ed.), Folia Fac. Sei. Nat. Univ. Masaryk. Brunensis, Mathematica, Dept. of Appl. Math., Masaryk University of Brno, Czech Rep., 2004, to appear. 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: veselySecon.muni.cz, URL: http://www.math.muni.cz/~vesely