Kapitola 1 Úvod Pojem náhodného procesu je zobecněním pojmu náhodné veličiny. Zatímco náhodná veličina je reálná funkce jedné proměnné - elementárního jevu, je náhodný proces reálnou funkcí dvou proměnných - elementárního jevu a jedné reálné proměnné. Tou obvykle bývá čas. Je možno uvést bezpočet příkladů náhodných procesů. Uveďme namátkou tyto ˇ Ve fyzikálních a technických vědách: seismický záznam v geofyzice, řada nejvyšších denních teplot v meterologii, 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, změny v tloušťce drátu v průběhu jeho délky, změny v počtu výzev na určité telefonní lince, atd. ˇ V biologických vědách: sledování různých parametrů znečištění ovzduší, EEG, EKG záznamy v medicině, procesy množení (např. bakterií), apod. ˇ Ve společenských vědách: změny v počtu obyvatelstva, procesy mortality a invalidity obyvatelstva, aj. ˇ V ekonomice změny poptávky po určitém výrobku, analýza vývoje kursu akcií na burze, objem zemědělské produkce, počet čekajících v letecké dopravě, atd. Tyto procesy, napohled rozmanité, lze jednotně popsat matematickým pojmem náhodného (stochastického) procesu. Ta část matematické statistiky, která se zmíněnými procesy za- bývá, se také nazývá statistickou dynamikou. Cílem analýzy náhodných procesů je konstrukce odpovídajícího modelu, což umožní poro- zumět mechanismu, na jehož základě jsou generovány sledované údaje. Znalost modelu dále umožňuje předpovídat budoucí vývoj a je-li možné řídit a optimalizovat činnost příslušného systému (vhodnou volbou vstupních parametrů a počátečních podmínek). 1 1.1. Historie 2 1.1 Historie K nejstarším záznamům ve tvaru časových řad patří astronomická pozorování. Grafická zná- zornění časových řad v podobě, na kterou jsme zvyklí teď, se začala objevovat na počátku 19. století (např. záznamy zemědělské produkce - známá Beveridgeova řada popisující cenový index pšenice v západní Evropě v letech 1500-1869). V 19. století a na počátku 20. století pod vlivem bouřlivého rozvoje fyziky byla snaha pohlí- žet na časové řady, jako kdyby byly generovány deterministickým mechanismem (s případnými chybami při měření jednotlivých hodnot). Hodně se proto pracovalo na rozpoznání periodických složek v časové řadě a jako matematický aparát se uplatnila především FOURIEROVA ana- lýza. Schusterův periodogram (1898) položil základy pro moderní spektrální analýzu časových řad. Tento deterministický přístup přetrval ještě během první čtvrtiny 20. století, přestože byl kritizován pro neschopnost vyjádřit nepravidelnosti v amplitudách a ve vzdálenosti mezi body zvratu (tj. lokálními extrémy) časových řad. Velký přelom v rozvoji disciplíny proto znamenal přístup Yulea (1927) a Slutského (1937), kteří ukázali, že řada tvořena vzájemně závislými pozorováními může být úspěšně generována navzájem nezávislými "šoky" (přesně řečeno řadou nekolerovaných veličin s nulovou střední hodnotou a konstantním rozptylem, pro kterou se vžilo označení bílý šum). Bylo opravdu vel- kým překvapením, když se potvrdilo, že pomocí modelů využívajících bílý šum (např. autore- gresních modelů) lze rozumně popsat časové řady, které se často vyznačují cyklickým pohybem s proměnnou amplitudou a proměnnou vzdáleností mezi body zvratu. O další rozvoj disciplíny směrem k praktickým aplikacím se zasloužili Brown, Durbin, Kendall, Pawloski, Wold a další. Kapitola 2 Náhodné procesy 2.1 Definice náhodného procesu Definice 2.1.1. Nechť je dán pravděpodobnostní prostor (, A, P), indexová množina T R a reálná funkce X : × T R definovaná pro a t T. Jestliže je X(, t) pro t T borelovsky měřitelná vzhledem k A, tj. pro B B a pro t T X-1 (B) = { : X(, t) B} A, kde B je -algebra borelovských podmnožin, pak tuto funkci nazýváme (n-rozměrným) ná- hodným procesem. Náhodný proces X(, t) při pevném se nazývá realizace (trajektorie) procesu. Pravděpodobnostní míru PX(B) = P(X-1 (B)) B B nazýváme rozdělení pravděpodobností náhodného procesu X(, t). Poznámka 2.1.2. Obdobně jako u náhodných veličin kdy místo X(), píšeme pouze X , u náhodných procesů místo {X(, t), , t T} píšeme {Xt, t T} . Definice 2.1.3. Pokud T = Z = {0, 1, 2, . . .} nebo T Z, mluvíme o procesu s diskrétním časem nebo o časové řadě či náhodné posloup- nosti. Pokud T = t1, t2 , kde - t1 < t2 +, říkáme, že {Xt, t T} je náhodný proces se spojitým časem. 3 2.1. Definice náhodného procesu 4 Definice 2.1.4. Dvojice (S, S), kde S je množina hodnot náhodných veličin Xt a S je -algebra podmnožin S, se nazývá stavový prostor procesu {Xt, t T}. Pokud náhodné veličiny Xt nabývají pouze diskrétních hodnot, říkáme, že jde o proces s diskrétními stavy. Navývá-li hodnot z nějakého intervalu, mluvíme o procesu se spojitými stavy. Rozdělení pravděpodobností PX náhodného procesu {Xt, t T} jednoznačně definuje rozdělení každého konečněrozměrného náhodného vektoru X = (Xt1 , . . . , Xtn ) , kde t1, . . . , tn jsou libovolné body z množiny T. Definice 2.1.5. Nechť Tn je množina všech vektorů Tn = {t = (t1, . . . , tn) : t1 t2 tn; ti T; i = 1, . . . , n}. Pak (konečně dimenzionální) distribuční funkcí náhodného procesu rozumíme funkci Ft(x) = Ft1,...,tn (x1, . . . , xn) = P(Xt1 x1, . . . , Xtn xn) = PXt ((-, x1 >, . . . , (-, xn >) pro t = (t1, . . . , tn) Tn a x = (x1, . . . , xn) Rn . Pro různá n a pro různé hodnoty t1, . . . , tn dostáváme celý systém distribučních funkcí, označme jej F, a který nemůže být úplně libovolný, ale zřejmě musí splňovat tzv. Kolmogorovy podmínky konzistence (K1) Podmínka symetrie: pro libovolnou permutaci i1, . . . , in čísel 1, . . . , n platí Fti1 ,...,tin (xi1 , . . . , xin ) = Ft1,...,tn (x1, . . . , xn) (K2) Podmínka konzistence: Ft1,...,tn,tn+1 (x1, . . . , xn, ) = Ft1,...,tn (x1, . . . , xn) Každému náhodnému procesu lze tedy přiřadit konzistentní systém distribučních funkcí. K da- nému konzistentnímu systému distribučních funkcí existuje vždy takový náhodný proces, že jeho systém distribučních funkcí je totožný se zadaným systémem, což říká následující věta. Věta 2.1.6. Kolmogorova věta K systému distribučních funkcí, které splňují Kolmogorovy podmínky konzistence, existuje pravděpodobnostní prostor (, A, P) a náhodný proces {Xt, t T} tak, že F je jeho systé- mem distribučních funkcí. Důkaz. Protože jde v podstatě o problém z teorie míry, nebudeme důkaz uvádět, lze ho najít např. v Neubrunn, T., Riečan, B.: Miera a integrál, Bratislava. Veda 1981. 2 2.2. PŘÍKLADY NÁHODNÝCH PROCESŮ 5 2.2 PŘÍKLADY NÁHODNÝCH PROCESŮ Příklad 2.2.1 (Sinusoida s náhodnou fází a amplitudou). Nechť A a jsou nezávislé náhodné veličiny A 0, Rs(0, 2). Náhodný proces {Xt, t T} může být definován Xt = r-1 A cos(t + ) , 0, r > 0. Příklad 2.2.2 (Binární proces). Binární proces {Xt, t = 1, 2, . . .} je posloupnost nezávislých alternativních náhodných proměn- ných: Xt A 1 2 . Příklad 2.2.3 (Náhodná procházka). Náhodná procházka {Xt, t = 1, 2, . . .} je definován takto X0 = 0 , Xt = Xt-1 + t , kde t jsou nezávislé stejně rozdělené náhodné veličiny s nulovou střední hodnotou a konstant- ním rozptylem, tj. t L(0, 2 ). Příklad 2.2.4 (Markovovův proces). Markovovův proces {Xt, t = 1, 2, . . .} je definován takto P(a < Xt b; X1 = x1, . . . , Xn = xn) = P(a < Xt b; Xn = xn), tj. chování procesu v budoucnosti závisí na jeho stavu v posledním časovém okamžiku, ale ne na cestě, kterou do tohoto stavu dospěl. 2.3 Stochastické procesy druhého řádu Definice 2.3.1. Řekneme, že náhodný proces {Xt, t T} je striktně stacionární, jestliže pro t = (t1, . . . , tn) Tn a pro = (t1 + h, . . . , tn + h) Tn platí Ft(x) = Ft1,...,tn (x1, . . . , xn) = F1,...,n (x1, . . . , xn) = F (x). Rovnost lze interpretovat tak, že základní pravděpodobnostní charakteristiky procesu se nemění při posunutí v čase. Definice 2.3.2. Existuje-li pro t T střední hodnota EXt, pak nazýváme funkci t = EXt střední hodnotu náhodného procesu. 2.3. Stochastické procesy druhého řádu 6 Definice 2.3.3. Náhodný proces {Xt, t T} nazýváme procesem druhého řádu, jestliže pro t T platí EX2 t < a říkáme, že náhodný proces má konečné druhé momenty. Poznámka 2.3.4. Pokud EX2 t < , pak ze Schwarzovy nerovnosti plyne E|Xt| (E12 E|Xt|2 ) 1 2 = (E|Xt|2 ) 1 2 < , tj. existuje střední hodnota EXt = t a rozptyl DXt = EX2 t - (EXt)2 - 2 t . Definice 2.3.5. Náhodný proces {Xt, t T} nazýváme stacionární ve střední hodnotě, pokud pro t T je střední hodnota konstantní, tj. EXt = . Pokud EXt = 0, náhodný proces nazýváme centrovaným. Definice 2.3.6. Uvažujme náhodný proces {Xt, t T}, který má konečné druhé momenty. Pak funkci (s, t) = C(Xs, Xt) = E(Xs - EXs)(Xt - EXt) nazveme autokovarianční funkcí. Tato reálná funkce dvou proměnných dává informaci o lineárním vztahu mezi jakoukoliv dvojicí náhodných veličin Xs a Xt. Definice 2.3.7. Náhodný proces {Xt, t T} se nazývá kovariančně stacionární, pokud pro t, s T platí (s, t) = (0, |s - t|) , což budeme také psát ve formě (s, t) = (s - t) , tj. autokovarianční funkce závisí na svých argumentech pouze prostřednictvím jejich rozdílů. 2.3. Stochastické procesy druhého řádu 7 Poznámka 2.3.8. Protože C(Xs, Xt) = C(Xt, Xs), pak pro kovariančně stacionární procesy platí (-t) = (t) a všechny náhodné veličiny Xt mají tentýž konečný rozptyl DXt = C(Xt, Xt) = (t - t) = (0). Ze Schwarzovy nerovnosti dále plyne |(t)| = |C(X0, Xt)| DX0DXt = (0) . Definice 2.3.9. Náhodný proces {Xt, t T} se nazývá (slabě) stacionární, je-li kovariančně stacio- nární, tj. (s, t) = (s - t) , pro t, s T a navíc stacionární ve střední hodnotě, tj. EXt = pro t T. Poznámka 2.3.10. Přívlastek slabě se většinou vynechává. Lze snadno ukázat, že je-li proces striktně stacionární, je také stacionární. Opačná implikace však neplatí. Definice 2.3.11. Nechť náhodný proces {Xt, t T} je stacionární. Označme (0) = 2 a zaveďme funkci (t) = (t) 2 = (t) (0) . Tuto funkci nazveme autokorelační funkcí stacionárního náhodného procesu Definice 2.3.12. Řekneme, že náhodný proces {t, t T} je bílým šumem ( White Noise), jestliže t jsou nekorelované náhodné veličiny s nulovou střední hodnotou, tj. Et = 0 Dt = 2 C(t, s) = 0 (s = t), značíme t WN(0, 2 ). 2.3. Stochastické procesy druhého řádu 8 Pokud jsou navíc nejen nekolerované, ale i nezávislé, značíme je symbolem IID(independent identical defined), píšeme t IID(0, 2 ). Věta 2.3.13. Náhodné procesy t WN(0, 2 ) a t IID(0, 2 ) jsou stacionárními náhodnými procesy. Důkaz. Zřejmý. 2 Definice 2.3.14. Náhodný proces {Xt, t T} se nazývá gaussovským (normálním), jestliže pro každé při- rozené n a libovolná čísla tj T, j = 1, . . . , n, je jeho n-rozměrná distribuční funkce Ft1,...,tn (x1, . . . , xn) distribuční funkcí n-rozměrného normálního rozdělení. Věta 2.3.15. Gaussův náhodný proces {Xt, t T} je stacionární, právě když je striktně stacionární. Důkaz. Je triviální a plyne z vlastností normálního rozdělení. 2 Definice 2.3.16. Řekneme, že náhodný proces {Xt, t T} splňuje lineární regresní model, pokud pro jeho střední hodnotu platí t T : EXt = t = m j=0 jfj(t), kde f0, . . . , fm jsou známé funkce definované na T, = (0, . . . , m) je neznámý vektor regresních koeficientů. 2.3. Stochastické procesy druhého řádu 9 2.3.1 Příklady Příklad 2.3.17. Mějme náhodnou veličinu U N(0, 1) a označme (u) = P(U u); u R. Pro t T R zaveďme Xt = U t . Vypočítejme distribuční funkci Ft(x) = P(Xt x). Ft(x) = P(Xt x) = P(U t x) = P(U x t ) = (x t ) t > 0; x R 0 t = 0; x 0 1 t = 0; x > 0 P(U x t ) = 1 - (x t ) t < 0; x R . Příklad 2.3.18. Mějme posloupnost nezávislých náhodných veličin, pro něž platí: Xt N(1, 1) je-li t liché Ex(1) je-li t sudé . Proces je (slabě) stacionární, neboť EXt = 1, DXt = 1 a (s, t) = 0 pro s = t (jsou nezávislé). Tento proces však není strikně stacionární, neboť pro liché a sudé t je distribuční funkce rozdílná. Příklad 2.3.19. Nechť stacionární náhodný proces {Yt, t Z}. Definujme Xt = Yt je-li t liché Yt + 1 je-li t sudé . I když X(t + h, t) = C(Xt+h, Xt) = C(Yt+h, Yt) = Y (h), tj. proces je kovariančně stacionární, přesto není (slabě) stacionární, protože střední hodnota není konstatní. Příklad 2.3.20. Mějme náhodný proces {Xt, t Z}, který je definován takto Xt = A cos t + B sin t , kde A, B jsou nekorelované náhodné veličiny, pro něž EA = EB = 0, DA = DB = 1, -, . 2.3. Stochastické procesy druhého řádu 10 Zjistěme, zda je proces (slabě) stacionární. Protože náhodné veličiny A a B jsou nekolerované, pak C(A, B) = E(A B) = 0. Vypočtěme nejprve střední hodnotu procesu: EXt = E(A cos t + B sin t) = cos t EA =0 + sin t EB =0 = 0 . Autokovarianční funkce: (t + h, t) = C(Xt+h, Xt) = E(Xt+h Xt) = E{[A cos (t + h) + B sin (t + h)] [A cos t + B sin t]} = E{A2 cos (t + h) cos t + AB cos (t + h) sin t + AB sin (t + h) cos t + B2 sin (t + h) sin t} = cos (t + h) cos t EA2 =1 + cos (t + h) sin t EAB =0 + sin (t + h) cos t EAB =0 + sin (t + h) sin t EB2 =1 = cos (t + h) cos t + sin (t + h) sin t = cos t[cos t cos h - sin t sin h] + sin t[sin t cos h + cos t sin h] = cos h[cos2 t + sin t2 t] = cos h = (h) Tedy (t + h, t) nezávisí na t, proto proces je (slabě) stacionární. Příklad 2.3.21 (Náhodná procházka). Nechť X0 = 0 ; Xt = Xt-1 + t , kde t jsou nezávislé stejně rozdělené náhodné veličiny t L(0, 2 ). Zjistěme, zda je proces (slabě) stacionární. Vypočtěme nejprve střední hodnotu procesu: EXt = E t i=1 i = t i=1 Ei =0 = 0 tj. konstatní střední hodnota. Autokovarianční funkce: (t + h, t) = C(Xt+h, Xt) = E(Xt+h Xt) = E t+h i=1 i t i=1 i = t i=1 E2 i = t 2 , závisí na t, tedy proces není (slabě) stacionární. 2.3. Stochastické procesy druhého řádu 11 Příklad 2.3.22 (MA(1) proces). Mějme Xt = t + t-1 , kde t WN(0, 2 ). Zjistěme, zda je proces (slabě) stacionární. Vypočtěme nejprve střední hodnotu procesu: EXt = E(t + t-1) = Et =0 + Et-1 =0 = 0 . Autokovarianční funkce (t + h, t) = (t + h, t) = C(Xt+h, Xt) = E(Xt+h Xt) = E(t+h + t+h-1)(t + t-1) = Et+h t + Et+h-1 t + Et+h t-1 + 2 Et+h-1 t-1 = = 2 (1 + 2 ) h = 0 2 h = 1 0 jinak = (h) Tedy MA(1) proces je (slabě) stacionární. Autokorelační funkce: (h) = 1 h = 0 1+2 h = 1 0 jinak . Příklad 2.3.23 (Markovovův skokový proces vzniku a zániku). Mějme náhodný proces {Xt, t R}, pro který platí P(Xt = 0) = + P(Xt = 1) = + > 0, > 0 P(Xs = x|Xt = y) = pxy(t - s) - < s t < , přičemž p00(t) = 1 - p01(t) = + + + e-(+)t t 0 p11(t) = 1 - p10(t) = + + + e-(+)t t 0. Konstanty a se interpretují jako intenzity vzniku a zániku. Rozhodněme, zda je proces stacionární. Výpočet střední hodnoty EXt = x=0,1 xP(Xt = x) = + a proces je tedy stacionární ve střední hodnotě. Výpočet autokovarianční funkce: platí (s, t) = C(Xs, Xt) = E(XsXt) - (EXs)(EXt), proto pro - < s t < nejprve počítejme E(XsXt) = x=0,1 y=0,1 xyP(Xs = x, Xt = y) = 1 1 P(Xs = 1, Xt = 1) = P(Xt = 1) P(Xs = 1|Xt = 1) p11(t-s) = + + + + e-(+)(t-s) 2.3. Stochastické procesy druhého řádu 12 tudíž (s, t) = + 2 + (+)2 e-(+)(t-s) - + 2 = (+)2 e-(+)(t-s) = (t - s) . Protože pro každou kovarianční funkci platí (t) = (-t), dostáváme pro - < t < (t) = (+)2 e-(+)|t|. Proces je tedy kovariančně stacionární a celkově je slabě stacionární. Výpočet autokorelační funkce: (t) = (t) (0) = e-(+)|t| . Korelace mezi náhodnými veličinami tohoto procesu pro |t| exponenciálně klesá k nule. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Autocorelation function (t)=exp{-(+)|t|}; =2.5; =2.7 Obrázek 2.1: Markovovův skokový proces vzniku a zániku: autokorelační funkce. 2.3. Stochastické procesy druhého řádu 13 2.3.2 Vlastnosti autokovariační funkce Třebaže v praktických situacích máme co činit jen s reálnými náhodnými veličinami, v teorii bývá výhodné pracovat někdy s komplexními náhodnými veličinami. Komplexní veličinou rozumíme veličinu X = Y + iZ, kde Y a Z jsou reálné náhodné veličiny. Komplexním náhodným procesem nazveme systém komplexních náhodných veličin {Xt, t T}. Mnoho dalších úvah se bude týkat právě komplexních procesů. Slovo ,,komplexní se bude vynechávat, když bude zřejmé ze souvislosti. Existují-li střední hodnoty EY a EZ, definuje se střední hodnota komplexní náhodné veličiny X = Y + iZ EX = EY + iEZ . Budeme se nyní zabývat základními vlastnostmi autokovarianční funkce (s, t). Přitom se samozřejmě předpokládá, že jde o proces s konečnými druhými momenty. Jelikož autokovarianční funkce procesu zůstává stejná při změně střední hodnoty, budeme také pro jednoduchost předpokládat, že střední hodnota procesu je rovna nule, tj. že proces je centrován. Nejprve uvedeme definici tzv. pozitivně semidefinitní funkce. Definice 2.3.24. Nechť f(s, t) je funkce dvou proměnných definovaná na T × T. Říkáme, že f je pozi- tivně semidefinitní, platí-li pro jakékoli přirozené číslo n, pro libovolná komplexní čísla c1, . . . , cn a libovolné body t1, . . . , tn T vztah n j=1 n k=1 f(tj, tk)cjck 0 . (2.1) Funkce jedné proměnné g(t), t T se nazývá pozitivně semidefinitní, platí-li pro každné přirozené n, libovolná komplexní čísla c1, . . . , cn a libovolné body t1, . . . , tn T a tj -tk T pro j, k = 1, . . . , n vztah n j=1 n k=1 g(tj - tk)cj ck 0 . (2.2) Ze souvislosti bude vždy patrno, zda jde o definici ve smyslu (2.1) nebo (2.2). Věta 2.3.25. Autokovarianční funkce (s, t) je pozitivně semidefinitní. 2.3. Stochastické procesy druhého řádu 14 Důkaz. Nechť {Xt, t T} je centrovaný proces s autokovarianční funkcí (s, t). Pak zřejmě platí 0 E n j=1 cjXtj 2 = E n j=1 cjXtj n j=k ck Xtk = n j=1 n j=k cj ckE(Xtj Xtk ) = n j=1 n j=k cj ck(tj, tk). Věta 2.3.26. Autokovarianční funkce (s, t) je hermitovsky symetrická, tj. pro s, t T platí (s, t) = (t, s) Důkaz. Plyne ze základních vlastností integrálu. Platí (s, t) = E(Xs Xt) = E(Xt Xs) = (t, s). 2 Věta 2.3.27. Je-li funkce (s, t) pozitivně semidefinitní a hermitovsky symetrická, existuje takový ná- hodný proces (dokonce normální), že (s, t) je jeho autokovarianční funkcí. Důkaz. Viz Doob, J.L.: Stochastic processes. New York, Wiley 1953 - kap. 2, §3. 2 Věta 2.3.28. Pro autokovarianční funkci (s, t) platí nerovnosti (s, s) 0 a |(s, t)| (s, s) (t, t). Důkaz. První nerovnost plyne z definice autokovarianční funkce a druhá ze Schwarzovy ne- rovnosti. 2 Lemma 2.3.29. Součet dvou pozitivně semidefinitních hermitovsky symetrických funkcí je opět funkce pozi- tivně semidefinitní a hermitovsky symetrická. Důkaz. Nechť f1(s, t) a f2(s, t) jsou pozitivně semidefinitní. Položme f(s, t) = f1(s, t) + f2(s, t). Pro libovolná komplexní čísla c1, . . . , cn platí n j=1 n j=k cjckf(tj, tk) = n j=1 n j=k cj ckf1(tj, tk) + n j=1 n j=k cj ckf2(tj, tk). Každý z obou výrazů na pravé straně je nezáporný. Musí být tudíž nezáporný i výraz vlevo, čímž je zaručena pozitivní semidefinitnost funkce f. Jsou-li f1 a f2 hermitovsky symetrické, je i funkce f hermitovsky symetrická; to je ihned zřejmé. 2 2.3. Stochastické procesy druhého řádu 15 Věta 2.3.30. Součet dvou autokovariačních funkcí je opět autokovarianční funkcí. Důkaz. Důkaz plyne z lemmatu 2.3.29 a z vět 2.3.25, 2.3.26 a 2.3.27. 2 Věta 2.3.31. Reálná část autokovarianční funkce je též autokovarianční funkcí. Imaginární část je auto- kovarianční funkcí jen tehdy, je-li rovna identicky nule. Důkaz. Nechť {Zt, t T} je komplexní náhodný proces s autokovariační funkcí (s, t). Bez újmy na obecnosti budeme předpokládat, že náhodný proces má nulovou střední hodnotu. Položme Zt = Xt + iYt, kde Xt a Yt jsou reálná náhodné veličiny. Z předpokladu EZt = plyne EXt = EYt = 0. Máme (s, t) = EZs Zt = E(Xs + iYs)(Xt - iYt) = EXsXt + EYsYt + i(EYsXt - EXsYt) Reálná část (s, t) je rovna EXsXt + EYsYt. Je tedy rovna součtu autokovariační funkce procesu {Xt, t T} a autokovariační funkce procesu {Yt, t T}. Podle věty 2.3.30 je autokovariančí funkcí. Imaginární část (s, t) je rovna EYsXt - EXsYt. V bodech s = t je rovna nule. Má-li být autokovariační funkcí, musí splňovat druhou nerovnost z věty 2.3.28. To je možné jen tehdy, je-li stále rovna nule. (Na druhé straně funkce identicky rovná nule je autokovariační funkcí - např. procesu, který je stále roven nule). 2 2.3.3 Spojitost, derivace a integrál procesu Budeme nyní předpokládat, že T = (a, b) R. Tam, kde budeme pojednávat o spojitosti a o derivaci procesu, pripustíme i možnost a = - nebo b = . V celé této kapitole budeme předpokládat, že proces je centrován. Tento předpoklad je dosti podstatný, neboť jinak by různé nespojitosti apod. mohly být zaviněny málo hladkým průběhem střední hodnoty. Samozřejmě budeme předpokládat, že proces má konečné druhé momenty. Všechny úvahy se týkají komplexních procesů. 2.3. Stochastické procesy druhého řádu 16 SPOJITOST NÁHODNÉHO PROCESU Pokud se zajímáme o spojitost procesu {Xt, t T} v bodě t0 T, budeme studovat chování náhodných veličin Xt při t t0. Jestliže Xt konvergují v nějakém smyslu k Xt0 , je možno mluvit o spojitosti procesu Xt v bodě t0. Z různých typů konvergencí se ukazuje v tomto případě jako nejužitečnější konver- gence podle kvadratického středu. Definice 2.3.32. Řekneme, že náhodný proces {Xt, t T} je spojitý podle středu v bodě t0 T , jestliže při t t0 konvergují Xt k Xt0 podle kvadratického středu, tj. když E|Xt - Xt0 |2 0 pro t t0. V tom případě píšeme Xt0 = l.i.m. tt0 Xt (zkratka z anglického "limit in the mean"). Je-li proces {Xt, t T} spojitý v každém bodě množiny T , říkáme stručně, že je spojitý. Poznámka 2.3.33. Z teorie pravděpodobnosti je známo, že konvergence podle kvadratického středu implikuje kon- vergenci podle pravděpodobnosti. Věta 2.3.34 (kritérium spojitosti procesu). Proces {Xt, t T} je spojitý právě tehdy, když je jeho autokovarianční funkce (s, t) spojitá v bodech (s, t), pro něž s = t. Důkaz. Bez újmy na obecnosti můžeme předpokládat, že proces je centrovaný. Je-li proces {Xt, t T} spojitý, pak platí |(s, t) - (s0, t0)| = |EXs Xt - EXs0 Xt0 | = = | E(Xs - Xs0 )( Xt - Xt0 ) (1) + EXs0 ( Xt - Xt0 ) (2) + E(Xs - Xs0 ) Xt0 (3) | trojúhel.ner. |E(Xs - Xs0 )( Xt - Xt0 )| + |EXs0 ( Xt - Xt0 )| + |E(Xs - Xs0 ) Xt0 | Schwarz.ner. E|Xs-Xs0 |2 E| Xt- Xt0 |2 0 1 2 + E|Xs0 |2 E| Xt- Xt0 |2 0 1 2 + E|Xs-Xs0 |2 E| Xt0 |2 0 1 2 pro s s0, t t0 (využili jsme vlastnosti spojitosti skalárního součinu). Funkce (s, t, ) je tudíž spojitá všude, a tedy také na diagonále s = t. 2.3. Stochastické procesy druhého řádu 17 Předpokládejme nyní, že (s, t, ) je spojitá na diagonále s = t. Máme E|Xs - Xt|2 = E(Xs - Xt)( Xs - Xt) = EXs Xs - EXs Xt - EXt Xs + EXt Xt = (s, s) - (s, t) - (t, s) + (t, t) Při pevném t a při s t z našeho předpokladu vyplývá (s, s) (t, t), (s, t) (t, t), (t, s) (t, t), takže E|Xs - Xt|2 0, tj. konverguje podle kvadratického středu. 2 DERIVACE NÁHODNÉHO PROCESU Derivaci náhodného procesu budeme definovat obdobně, jako se definuje derivace funkce. Definice 2.3.35. Řekneme, že náhodný proces {Xt, t T} má v bodě t0 T derivaci Xt0 , jestliže platí l.i.m. h0 Xt0+h - Xt0 h = Xt0 pro t0 + h T. Má-li náhodný proces {Xt, t T} derivaci ve všech bodech t T, říkáme stručně, že má derivaci. Věty, které dávají nutnou a postačující podmínku pro existenci derivace náhodného procesu, lze najít v knize Anděl, J.: Statistická analýza časových řad. Praha. SNTL 1976 INTEGRÁL NÁHODNÉHO PROCESU Definice 2.3.36. Nechť T = (a, b) je konečný nedegenerovaný interval. Nechť T1, T2 T, kde T1 < T2. Pomocí dělení Dn intervalu T1, T2 takového, že T1 = t0 < t1 < < tn = T2, utvoříme součet In = n-1 k=0 Xtk (tk+1 - tk). Označme n = max 0kn-1 (tk+1 - tk) normu dělení Dn. Jestliže posloupnost In konverguje podle kvadratického středu k náhodné veličině I pro n při každé takové posloupnosti dělení Dn, že n 0, náhodná veličina I se nazývá Riemannův integrál procesu {Xt, t T} na intervalu T1, T2 . Píšeme pak I = T2 T1 Xt dt. 2.3. Stochastické procesy druhého řádu 18 Věta 2.3.37 (kritérium existence integrálu procesu). Má-li centrovaný náhodný proces {Xt, t T} s konečnými druhými momenty autokovarianční funkci (s, t), pak Rie- mannův integrál T2 T1 Xt dt existuje právě tehdy, když existuje Riemannův integrál T2 T1 T2 T1 (s, t) ds dt. Důkaz. Lze najít v knize Anděl, J.: Statistická analýza časových řad. Praha. SNTL 1976. 2 V praxi samozřejmě není nutno počítat T2 T1 T2 T1 (s, t) ds dt, abychom se přesvědčili, že exis- tuje T2 T1 Xt dt. Stačí např. zjistit, že (s, t) je spojitá, protože pak už Riemannův integrál T2 T1 T2 T1 (s, t) ds dt musí existovat vzhledem k omezenosti uzavřené množiny T1, T2 × T1, T2 . Poznámka 2.3.38. Někdy je třeba užívat integrál procesu v jiném smyslu. Chápeme-li proces X(, t) jako funkci dvou proměnných, můžeme pro každé definovat Lebesgueův integrál T2 T1 X(, t) dt. Musíme mít ovšem zaručeno, že tento integrál vůbec existuje. V definici náhodného procesu jsme totiž nečinili žádné předpoklady o chování funkcí X(, ), tj. o chování množiny funkcí proměnné t při pevném , takže ani nemáme zaručeno, že tyto funkce jsou vůbec měřitelné. Pokud se stane, že při každém je X(, ) spojitou funkcí proměnné t, pak je samozřejmě existence Lebesgueuva integrálu zaručena. 2.3.4 Spektrální rozklad autokovariančních funkcí stacionárních pro- cesů V celém tomto odstavci budeme předpokládat, že náhodný proces {Xt, t T} je stacionární, centrovaný a druhého řádu (tj. s konečnými druhými momenty). Významnou vlastností stacionárních náhodných procesů je vlastnost, že jeho autokovariační funkci lze vyjádřit jako (nespočetný) součet harmonických funkcí s různými frekvencemi a amplitudami. Věta 2.3.39 (Herglotzova věta). Je-li {Xt, t Z} stacionární posloupnost, pak se její autokovarianční funkce (t) dá vyjádřit ve tvaru (t) = - eit dF(), kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. 2.3. Stochastické procesy druhého řádu 19 Důkaz. Protože (t) je pozitivně semidefinitní, platí pro každé přirozené n, libovolná kom- plexní čísla c1, . . . , cn C a libovolné body t1, . . . , tn T a tj - tk T j, k = 1, . . . , n vztah n j=1 n k=1 (tj - tk)cj ck 0. (2.3) Položme tj = j, cj = e-ij , kde - . Definujme fn() = 1 2n n j=1 n k=1 (j - k)e-i(j-k) . Vzhledem k (2.3) platí fn() 0 pro - . Po substituci j - k = s dostáváme vyjádření fn() = 1 2n n-1 s=-n+1 (n - |s|) (s)e-is a celkově fn() = 1 2 n-1 s=-n+1 1 - |s| n (s)e-is |s| < n, 0 |s| n. (2.4) Uvědomíme-li si, že platí - eikx dx = 2 k = 0, 0 k Z, k = 0, (2.5) což lze ověřit třeba i rozkladem eikx = cos kx + i sin kx, počítejme - eit fn() d = - eit 1 2 n-1 s=-n+1 1 - |s| n (s)e-is d = 1 2 n-1 s=-n+1 1 - |s| n (s) - ei(t-s) d =2 pro t=s; =0 jinak = (t) 1 - |t| n |t| < n, 0 |t| n. Označme Fn() = - fn(x) dx. (2.6) 2.3. Stochastické procesy druhého řádu 20 Protože fn(x) je nezáporná, je Fn() neklesající funkce. Je jasné, že Fn(-) = 0. Dále s využitím vztahu (2.5) počítejme Fn() = - fn() d = - 1 2 n-1 s=-n+1 1 - |s| n (s)e-is d = 1 2 n-1 s=-n+1 1 - |s| n (s) - e-is d =2 pro s=0; =0 jinak = (0) . Položme Fn() = 0 pro < -, Fn() = (0) pro > . Podle 1. Hellyho věty lze z posloupnosti funkcí {Fn()} vybrat podposloupnost {Fnk ()}, která konverguje k nějaké neklesající zprava spojité funkci F() ve všech bodech spojitosti této limitní funkce (jde o tzv. konvergenci v podstatě). Pro F() musí zřejmě platit F(- - 0) = 0 a F( + 0) = (0); ovšem vzhledem ke spojitosti zprava musí dokonce platit F() = (0). Podle 2. Hellyho věty (vzhledem k ohraničenosti funkce eit ) platí pro každé celé t lim nk - eit dFnk () = - eit dF(). Dále platí (t) = lim nk (t) 1 - |t| nk = lim nk - eit fnk () d = lim nk - eit dFnk () = - eit dF() . 2 Má-li funkce F() skok v bodě = , pozměníme ji tak, že velikost tohoto skoku dáme odtud do bodu = -. Tím se hodnota integrálu nezmění a docílili jsme však toho, že F() = (0). 2.3. Stochastické procesy druhého řádu 21 Věta 2.3.40 (Bochnerova věta). Je-li {Xt, t R} stacionární proces spojitý podle středu, pak se jeho autokovarianční funkce (t) dá vyjádřit ve tvaru (t) = - eit dF(), kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. Důkaz. Opět vyjdeme z pozitivně semidefinitní funkce (t), tj. pro každé přirozené n, libovolná komplexní čísla c1, . . . , cn C a libovolné body t1, . . . , tn T a tj - tk T j, k = 1, . . . , n platí vztah n j=1 n k=1 (tj - tk)cj ck 0. (2.7) Dosadíme-li do (2.7) tj = j N , kde N N cj = e-ij , kde - , dostaneme, že funkce fN n () = 1 2n n j=1 n k=1 j - k N e-i(j-k) = 1 2n n-1 s=-n+1 (n - |s|) s N e-is = 1 2 n-1 s=-n+1 1 - |s| n s N e-is je pro každé -, nezáporná. Proto každá funkce FN n () = - fN n (x)dx je neklesající a zřejmě FN n (-) = 0. S využitím vztahu (2.5) platí FN n () = - fN n () d = - 1 2 n-1 s=-n+1 1 - |s| n s N e-is d = 1 2 n-1 s=-n+1 1 - |s| n s N - e-is d =2 pro s=0; =0 jinak = (0) . 2.3. Stochastické procesy druhého řádu 22 Podle 1. Hellyho věty lze z posloupnosti funkcí F N 1 (), FN 2 (), . . . vybrat podposloupnost {F N nk ()}, která konverguje k nějaké neklesající zprava spojité funkci F N (), pro kterou platí FN (- - 0) = 0 a FN ( + 0) = FN () = (0). Položme N (t) = N -N eit dFN N . (2.8) Z 2. Hellyho věty plyne, že pro každé celé k platí N k N = N -N ei k N dFN N = - eikx dFN (x) (2.9) = lim nj - eikx fN nj (x)dx = lim nj k N 1 - |k| nj = k N . Ke každému t R a každému přirozenému N existuje celé k (závislé na t a N) tak, že platí 0 t - k N < 1 N . (2.10) Při této volbě platí k N t pro N . Ze spojitosti procesu {Xt, t T} plyne spojistost funkce (t). Proto s přihlédnutím k (2.9) máme (t) = lim N k N = lim N N k N . (2.11) Ze zřejmé rovnosti N (t) = N k N + N (t) - N k N plyne vzhledem k (2.10) lim N N (t) = (t) + lim N N (t) - N k N (2.12) Dosadíme-li z (2.8) a z první části (2.9), máme N (t) - N k N = N -N ei k N ei(t- k N ) - 1 dFN N Nechť = t - k N . Z (2.10) plyne 0 N < 1. Ze Schwarzovy nerovnosti plyne N (t) - N k N 2 N -N ei k N 2 dFN N =F N ()=(0) N -N ei(t- k N ) - 1 2 dFN N 2.3. Stochastické procesy druhého řádu 23 První integrál na pravé straně je roven (0). Protože ei - 1 2 = ei - 1 e-i - 1 = 2 - 2 cos , máme N (t) - N k N 2 2(0) N -N (1 - cos ) dFN N = 2(0) - (1 - cos Nx) dFN (x). Z nerovnosti cos b cos pro - , 0 b < 1, plyne 1 - cos Nx 1 - cos x, takže N (t) - N k N 2 2(0) - (1 - cos Nx) dFN (x) = 2(0) (0) - reálná část 1 N . Vzhledem ke spojitosti (t) a k tomu, že (0) je nezáporné (a tudíž reálné) číslo, je lim N N (t) - N k N = 0. Podle (2.12) tedy platí (t) = lim N N (t). (2.13) 2 Položme FN () = 0 pro < -N FN N pro -N N (0) pro > N. Je-li (0) = 0, je tvrzení věty zřejmé (stačí vzít F() 0). Nechť tedy (0) > 0. Pak můžeme bez újmy na obecnosti předpokládat (0) = 1, neboť jinak bychom místo (t) mohli vyšetřovat (t) (0) . Vztah (2.8) můžeme napsat jako N (t) = - eit dFN (). 2.3. Stochastické procesy druhého řádu 24 Jelikož FN () má všechny vlastnosti distribuční funkce, je {N (t)} posloupnost charakteristic- kých funkcí, která ve všech bodech t podle (2.13) konverguje ke spojité funkci (t). Podle známé věty o vztahu mezi charakteristickými a distribučními funkcemi je (t) také charakteristická funkce, která odpovídá nějaké distribuční funkci F() takové, že FN F v podstatě. Tudíž (t) = - eit dF(), což je vztah, který bylo třeba dokázat. Jelikož F je distribuční funkce, platí F(-) = 0 a F() = 1 = (0). Vzorci (t) = - eit dF() resp. (t) = - eit dF() se říká spektrální rozklad kovarianční funkce. Funkce F() se nazývá spektrální distri- buční funkce. Je-li F() absolutně spojitá, pak existuje taková funkce f(), že pro náhodné stacionární posloupnosti, resp. pro stacionární náhodné procesy platí F() = - f(x)dx resp. F() = - f(x)dx. (2.14) Jelikož F() je neklesající, je f() skoro všude nezáporná. Je-li třeba, pozměníme ji na množině míry nula tak, aby byla všude nezáporná. Tím se integrál (2.14) nezmění. Funkce f() se nazývá spektrální hustota. Existuje-li spektrální hustota, pak můžeme psát (t) = - eit f()d (2.15) resp. (t) = - eit f()d . (2.16) Všimněme si ještě, zda a jak se dá na základě nějaké jednoduché vlastnosti kovarianční funkce (t) poznat, zda vůbec spektrální hustota existuje. Věta 2.3.41. K existenci spektrální hustoty stacionární náhodné posloupnosti stačí, aby pro její kovari- anční funkci platilo t=- |(t)| < K existenci spektrální hustoty spojitého stacionární náhodného procesu stačí, aby pro její kovarianční funkci platilo - |(t)|dt < . 2.3. Stochastické procesy druhého řádu 25 Důkaz. Viz Gichman, I.I., Skorochod, A.V.: Teorija slučajnych processov. Moskva. Nauka 1971. 2 V následujících dvou větách je zodpovězena otázka, jak vypočítat spektrální hustotu z ko- varianční funkce. Věta 2.3.42. Existuje-li spektrální hustota f() stacionární posloupnosti a má-li variaci konečnou na -, , pak platí f() = 1 2 t=- e-it (t) (2.17) ve všech bodech spojitosti funkce f(), což je skoro všude vzhledem k Lebesqueově míře. Důkaz. Ze vzorce (2.15) vidíme, že až na normující konstantu 1 2 jsou (t) Fourierovy ko- eficienty funkce f() vzhledem k ortogonálnímu systému funkcí {e-it }. Zbytek tvrzení plyne z faktu, že funkce s konečnou variací má nejvýše spočetně bodů nespojitosti. 2 Věta 2.3.43. Existuje-li spektrální hustota f() spojitého stacionárního procesu a je-li autokovarianční funkce absolutně integrovatelná, tj. - |(t)|dt < , pak f() = 1 2 - e-it (t) dt. (2.18) Důkaz. Ze vzorce (2.16) vidíme, že až na normující konstantu 1 2 jsou je mezi (t) a f() stejný vztah jako mezi charakteristickou funkcí a hustotou rozdělení. Proto lze přímo převzít vzorec pro výpočet hustoty z charakteristické funkce. 2 Věta 2.3.44. Spektrální hustota f() reálného spojitého stacionárního procesu nebo reálné stacionární posloupnosti je sudá funkce v tom smyslu, že pro ni platí f() = f(-) (2.19) skoro všude vzhledem k Lebesqueově míře. Důkaz. Nechť {Xt, t T} je spojitý stacionární proces. Jelikož je reálný, platí pro každé t T, že (t) = (-t). Proto vzhledem k (2.16) (t) = - eit f()d = - e-it f()d = (-t). 2.3. Stochastické procesy druhého řádu 26 Substitucí se snadno zjistí, že pravá strana je rovna - eit f(-)d takže - eit f()d = - eit f(-)d. (2.20) Je-li f() = 0 skoro všude, je tvrzení věty zřejmé. Předpokládejme tedy, že - f()d = C > 0. Bez újmy na obecnosti můžeme položit C = 1 (jinak stačí místo f() uvažovat f() C ). Pak vzo- rec (2.20) ukazuje, že charakteristické funkce příslušející hustotám f() a f(-) jsou totožné. Vzhledem k vzájemně jednoznačnému vztahu mezi rozdělením pravděpodobnosti a charakte- ristickou funkcí odtud vyplývá tvrzení věty. Pro stacionární posloupnosti je důkaz obdobný. 2 2.3. Stochastické procesy druhého řádu 27 2.3.5 Příklady Mějme náhodný proces Yt = r cos(t0 + ) , kde r R . . .amplituda, 0 [0, ] . . .frekvence Rs(-, ) . . . náhodná fáze s hustotou f(x) = 1 2 x [-, ] 0 jinak . Není těžké ukázat, že EYt = 0 a autokovarianční funkce (h) = C(Yt, Yt+h) = 1 2 r2 cos h0, tedy rozptyl je roven D(Yt) = 1 2 r2 . Ekvivalentně lze psát Yt = U cos t0 + V sin t0, kde U = r cos V = -r sin náhodné amplitudy s vlastnostmi EU = EV = 0, C(U, V ) = 0, DU = DV = 1 2 r2 = 2 . Zobecníme-li tento proces pro nějaké přírozené číslo k Yt = k j=1 rj cos(tj + j) = k j=1 Uj cos tj + Vj sin tj, kde Uj = rj cos j Vj = -rj sin j náhodné amplitudy s vlastnostmi EUj = EVj = 0, C(Ui, Vj) = 0, DUj = DVj = 1 2 r2 j = 2 j . pak dostaneme EYt = 0, (h) = C(Yt, Yt+h) = k j=1 1 2 r2 j cos hj, DYt = k j=1 1 2 r2 j = k j=1 2 j . Se zobecněním můžeme jít ještě dále a definovat (konvergují-li příslušné nekonečné řady) Yt = j=1 rj cos(tj + j) = j=1 Uj cos tj + Vj sin tj a bude platit EYt = 0, (h) = C(Yt, Yt+h) = j=1 1 2 r2 j cos hj, DYt = j=1 1 2 r2 j = j=1 2 j . 0.5236 1.0472 1.5708 0 2 4 6 8 10 12 14 p(j )=j 2 j=1,2,3 1 2 =0.5 2 2 =12.5 3 2 =4.5 0.5236 1.0472 1.5708 0 2 4 6 8 10 12 14 I()=j 2 /j j=1,2,3 1 2 =0.5 2 2 =12.5 3 2 =4.5 Example: Yt =1cos(t/6+1 ) + 5cos(t/3+2 ) + 3cos(t/2+3 ) j ~ Rs (-,) j=1,2,3 -15 -10 -5 0 5 10 15 -8 -6 -4 -2 0 2 4 6 8 10 (t)=cov(Y ,Y +t ) Obrázek 2.2: Harmonický proces: autokovarianční funkce a diskrétní spektrum definované vztahem p(j) = F(j)-F(j-1) = 2 j , j = 1, 2, . . . vyjádřené pomocí čárkového a sloupcového grafu. Na závěr tohoto odstavce ukážeme graf s několika příklady autokovariančních funkcí a jim 2.3. Stochastické procesy druhého řádu 28 odpovídajících spektrálních hustot. -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 1.2 2 =1 (t)=C(Y ,Y+t ) White noise t ~ WN(0,2 ) t=0, 1,2,... -6 -4 -2 0 2 4 6 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 2 /(2)=0.15915 spectral density f()= 2 /(2)I[- ] White noise t ~ WN(0,2 ) t=0, 1,2,... -4 -3 -2 -1 0 1 2 3 4 0 0.5 1 1.5 2 2.5 (t)=C(Y ,Y+t ) -15 -10 -5 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 spectral density f()=2(1-cost0 )/(t0 2 ) -6 -4 -2 0 2 4 6 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 (t)=C(Y ,Y+t )=(sin2 t-sin1 t)./[t(2 -1 )] Autocovariance and spectral density 2 3 4 5 6 7 0 0.5 1 1.5 spectral density f()=I[1 ,2 ] 1 = 2 =2 Obrázek 2.3: Příklady autokovarianční funkce a spektrální hustoty. 2.3. Stochastické procesy druhého řádu 29 2.3.6 Hilbertovy prostory spjaté s procesy druhého řádu Připomeňme si z funkcionální analýzy následující pojmy a vlastnosti: UNITÁRNÍ PROSTORY. Komplexní lineární prostor H se nazývá unitární, jestliže pro každé dva prvky x a y z H existuje komplexní číslo x, y , nazývané skalární či vnitřní součin, tak že pro každé x, y, z H a C platí (a) x, y = y, x ; (b) x + y, z = x, z + y, z ; (c) x, y = x, y ; (d) x, x 0; (e) x, x = 0 x = 0. NORMA. V unitárním prostoru H definujeme normu vztahem x = x, x . CAUCHY-SCHWARZOVA NEROVNOST. V unitárním prostoru platí | x, y | x y a | x, y | = x y x = x,y y,y y. ORTOGONALITA. Řekneme, že x a y z unitárního prostoru H jsou ortogonální, pokud platí x, y = 0 a značíme xy. ORTOGONÁLNÍ A ORTONORMÁLNÍ MNOŽINY. Řekneme, že množina M H je ortogonální, jestliže pro každé různé prvky x, y M platí xy. Jestliže navíc pro x M platí x = 1, pak množina M se nazývá ortonormální. Poznámka: Je-li M ortogonální množina, pak množina { x x : x M} je ortonormální. VLASTNOSTI NORMY. Mějme unitární prostor H s normou definovanou vztahem x = x, x . Pak pro každé x, y H a pro každé C platí (a) x + y 2 = x 2 + y 2 + x, y + y, x ; (b) x + y x + y (tzv. trojúhelníková nerovnost); (c) x = || x ; (d) x 0; (e) x = 0 x = 0; (f) x + y 2 + x - y 2 = 2 x 2 + 2 y 2 (tzv. rovnoběžníková rovnost); KONVERGENCE PODLE NORMY. Řekneme že posloupnost prvků {xn} z unitárního prostoru H konverguje podle normy k x H, jestliže xn -x 0 pro n . SPOJITOST SKALÁRNÍHO SOUČINU. Jestliže {xn} a {yn} jsou prvky z unitárního prostoru H takové, že xn - x 0 a yn - y 0 pro n , pak platí (a) xn x a (b) xn, yn x, y pro n . CAUCHYOVSKÁ POSLOUPNOST. Řekneme, že posloupnost prvků {xn} z unitárního prostoru H je cauchyovská, pokud xn - xm 0 pro n, m . 2.3. Stochastické procesy druhého řádu 30 HILBERTOVY PROSTORY. Hilbertův prostor je úplný unitární prostor, tj. takový, ve kterém každá cauchyovská posloupnost {xn} konverguje podle normy k nějakému prvku x H, tj. xn - xm ---- n,m 0 x H : xn - x --- n 0. UZAVŘENÝ PODPROSTOR. Řekneme, že lineární podprostor M Hilbertova prostoru H je uzavřeným podprostorem H, jestliže M obsahuje všechny limitní body, tj. jestliže platí, že xn - x 0, pak x M. ORTOGONÁLNÍ KOMPLEMENT. Ortogonální komplement množiny M je množina M všech prvků H, které jsou ortogonální ke každému prvku z M. Tedy M = {y H : x, y = 0, tj. xy, x M} . PROJEKČNÍ VĚTA. Jestliže M je uzavřený podprostor Hilbertova prostoru a x H, pak (a) existuje jediný prvek ^x M, takový, že x - ^x = inf yM x - y (b) ^x M a x - ^x = inf yM x - y ^x M a (x - ^x) M . Prvek ^x se nazývá ortogonální projekcí prvku x z H do M a značíme ^x = PM(x) a zobrazení PM : H M se nazývá projekcí H do M. VLASTNOSTI PROJEKCE. Nechť H je Hilbertův prostor a PM je projekcí H do M. Pak pro každé x, y, xn H a pro každé , C platí (a) Každý prvek x H má jedinou reprezentaci jako součet prvku z M a prvku z M , tj. x = PM(x) + (I - PM)(x), kde I značí identické zobrazení (b) PM(x + y) = PM(x) + PM(y) (c) x 2 = PM(x) 2 + (I - PM)(x) 2 (d) xn - x --- n 0 PM(xn) --- n PM(x) (e) x M PM(x) = x (f) x M PM(x) = 0 (g) jestliže M1 a M2 jsou dva podprostory H takové, že M1 M2, pak PM1 (PM2 (x)) = PM1 (x) UZÁVĚR. Nechť M je podprostor Hilbertova prostoru H. Uzávěrem M (také budeme zna- čit sp(M), anglicky ,,closed span ) množiny M nazveme nejmenší uzavřenou množinu obsahující M. Poznámka: Platí M = sp(M) = {x H : xn - x --- n 0, xn L(M)}, kde L(M) je množina všech lineárních kombinací prvků množiny H, tzv. lineární obal množiny H. PROJEKCE NA KONEČNÉ ORTONORMÁLNÍ MNOŽINĚ. Jestliže {e1, . . . , en} je ortonormální podmnožina Hilbertova prostoru H a M = sp{e1, . . . , en}, pak pro každé x H platí (a) PM(x) = n i=1 x, ei ei (b) PM(x) 2 = n i=1 | x, ei |2 2.3. Stochastické procesy druhého řádu 31 (c) x - n i=1 x, ei ei 2 x - n i=1 iei 2 pro 1, . . . , n C (d) x - n i=1 x, ei ei 2 = x - n i=1 iei 2 i = x, ei i = 1, . . . , n (e) n i=1 | x, ei |2 x (tzv. Besselova nerovnost) Poznámka: koeficienty i = x, ei se nazývají Fourierovy koeficienty vzhledem k množině {e1, . . . , en}. SEPARABILITA. Hilbertův prostor H nazveme separabilním, právě když H = sp{et, t T} , kde T je spočetná indexová množina. ORTONORMÁLNÍ REPREZENTACE V SEPARABILNÍM HILBERTOVĚ PROSTORU. Nechť H = sp{e1, e2, . . .} je separabilní Hilbertův prostor, kde {ei} i=1 je ortonormální množina. Pak pro každé x, y H platí (a) Množina všech konečných lineárních kombinací {e1, . . . , en} je hustá, tj. pro x H a > 0 k N a 1, . . . , n C, že platí x - n i=1 iei < (b) x = i=1 x, ei ei pro x H, tj. x - n i=1 x, ei ei --- n 0 (c) x 2 = i=1 | x, ei |2 (tzv. Parsevalova identita) (d) x, y = i=1 x, ei ei, y (e) x = 0 x, ei = 0 i = 1, 2, . . . Zaveďme následující prostory náhodných veličin: ˇ Označme L2 (, A, P) , resp. L2 C(, A, P) množinu všech reálných, resp. komplex- ních náhodných veličin definovaných nad týmž pravděpodobnostním prostorem (, A, P), které mají konečné druhé momenty, tj. platí EX2 < , resp. E|X|2 < . Do tohoto prostoru zahrnujeme také všechny konstanty z R, resp. z C, které považujeme za náhodné veličiny s nulovým rozptylem. V tomto prostoru vytvoříme třídy ekvivalentních náhodných veličin takto: řekneme, že dvě náhodné veličiny jsou ekvivalentní, pokud se liší jen na množině míry nula. Zřejmě X a Y jsou ekvivalentní právě tehdy, platí-li E|X - Y |2 = 0. V takto definovaném prostoru tříd ekvivalentních náhodných veličin definujeme pro každé X, Y L2 (, A, P), resp. X, Y L2 C(, A, P), skalární součin předpisem X, Y = E(XY ) resp. X, Y = E(X Y ) a odpovídající normu X = X, X = EX2 resp. X = X, X = E(X X) = E|X|2. Přechod ke třídám je nutný proto, abychom zaručili platnost požadavku x, x = 0 x = 0. 2.3. Stochastické procesy druhého řádu 32 Věta 2.3.45. Prostory L2 (, A, P) a L2 C(, A, P) jsou Hilbertovy prostory. Důkaz. Viz Brockwell, P.J. a Davis, R.A.: Theory and methods, 2-nd ed., Springer-Verlag, New York, 1991, §2.10. 2 Již dříve jsme definovali pojem spojitosti podle středu v bodě t0 T takto E|Xt - Xt0 |2 0 pro t t0. což jsme značili Xt0 = l.i.m. tt0 Xt (zkratka z anglického "limit in the mean") a je-li proces {Xt, t T} spojitý v každém bodě množiny T , říkali jsme stručně, že je spojitý. Tutéž spojitost můžeme definovat i pomocí výše uvedené normy takto Xt - Xt0 2 = E|Xt - Xt0 |2 0 pro t t0 a pro každý uzavřený podprostor M L2 C(, A, P) díky projekční větě můžeme definovat nejlepší střední kvadratickou predikci prvku Y L2 C(, A, P) pomocí M. Definice 2.3.46. Jestliže M je uzavřený podprostor H, kde H = L2 (, A, P), resp. H = L2 C(, A, P), pak nejlepší střední kvadratická predikce Y H v M je prvek ^Y M takový, že Y - ^Y 2 = inf ZM Y - Z 2 = inf ZM E|Y - Z|2 tj. ^Y = PM(Y ). Nyní definujme Definice 2.3.47. Jestliže M je uzavřený podprostor H, kde H = L2 (, A, P), resp. H = L2 C(, A, P), a X H, pak definujme podmíněnou střední hodnotu při dané M předpisem EMX = E(X|Y M) = PM(X). Dále definujme Definice 2.3.48. Nechť X, Z1, . . . , Zn H, kde H = L2 (, A, P), resp. H = L2 C(, A, P). Pak podmíněná střední hodnota X při daném náhodném vektoru Z = (Z1, . . . , Zn) je dána vztahem E(X|Z) = EM(Z)X = E(X|Y M(Z)), kde M(Z) je uzavřený podprostor všech náhodných veličin (Z) z H, které jsou borelovskou funkcí náhodného vektoru Z, tj. : Rn C, resp. : Cn R. 2.3. Stochastické procesy druhého řádu 33 Výpočet podmíněných středních hodnot bývá však velmi obtížný. Z těchto důvodů se ome- zujeme většinou na jednodušší podprostor, a to M = sp{1, Z1, . . . , Zn} = {1, Z1, . . . , Zn}. Připomeňme vlastnost predikce ^x M prvku x H. Platí ^x = PM(x) M ^x M (x - ^x) M tj. pro každé y M platí x - ^x, y = x, y - ^x, y = 0 a odtud dostaneme tzv. projekční rovnice ^x, y = x, y . Dále již uvažujme Hilbertův prostor H = L2 (, A, P) a jeho podprostor M = sp{1, Z1, . . . , Zn} = {1, Z1, . . . , Zn}, kde Z1, . . . , Zn L2 (, A, P). Pak projekce je dána vztahem ^X = PM(X) = EMX = arg inf Y M X - Y 2 = arg inf Y M E(X - Y )2 a projekční rovnice jsou tvaru E (Y EMX) = E(Y X) . Pro každý prvek z M (tedy i pro 1, Z1, . . . , Zn) platí tyto rovnice, tj. pro Y = 1 máme E(1 EMX) = E(1 X) E(1 n i=0 iZi) = EX n i=0 iEZi = EX a pro Y = Zj, j = 1, . . . , n dostaneme E(Zj EMX) = E(Zj X) E(Zj n i=0 iZi) = E(ZjX) n i=0 iE(ZiZj) = E(ZjX) Celkem dostáváme systém n + 1 rovnic. Definujme proto nyní nejlepší lineární predikci pomocí obecnějších systémů náhodných ve- ličin druhého řádu {Zt, t T}. 2.3. Stochastické procesy druhého řádu 34 Definice 2.3.49. Nechť X H a pro každé t T také Zt H, kde T je indexová množina, H = L2 (, A, P), resp. H = L2 C(, A, P). Pak nejlepší lineární predikci náhodné veličiny X pomocí {Zt, t T} rozumíme Psp{Zt,tT }(X) . Vidíme, že při hledání nejlepší lineární predikce vystačíme se znalostí kovarianční funkce a není třeba znát ani momenty vyšších řádů. PREDIKCE U NORMÁLNĚ ROZDĚLENÝCH NÁHODNÝCH VELIČIN. Je-li sdružené rozdělení náhodných veličin X, Z1, . . . , Zn normální, tj. (X, Z1, . . . , Zn) Nn+1(, ) , kde = X Z1 Z2 ... Zn = X Z a = 2 X | XZ1 XZn XZ1 | 2 Z1 Z1Z2 Z1Zn XZ2 | Z1Z2 2 Z2 Z2Zn ... | ... ... ... ... XZn | 2 Zn = 2 X XZ XZ ZZ Pak rozdělení náhodné veličiny X při daném Z má opět normální rozdělení X|Z N X|Z, 2 X|Z , kde X|Z = X + XZ-1 ZZ(Z - Z) a 2 X|Z = 2 X + XZ-1 ZZZX Odtud vidíme, že podmíněná střední hodnota je lineární funkcí náhodného vektoru Z = (Z1, . . . , Zn) . To znamená, že v tomto případě je nejlepší lineární predikce totožná s optimální pre- dikcí (ve smyslu minimální střední kvadratické odchylky) založenou na podmíněných středních hodnotách. 2.4. Odhady středních hodnot a autokovariancí 35 2.4 Odhady středních hodnot a autokovariancí Stochastický proces je matematickým modelem reálného děje náhodného charakteru, který probíhá nepřetržitě v čase. Můžeme jej však pozorovat jen v konečných časových intervalech a na základě těchto pozorování určit odhady hodnot charakteristik tohoto procesu - střední hodnoty, rozptylu, kovarianční funkce, atd. Jestliže máme k dispozici dostatečný počet pozorování realizací náhodného procesu, můžeme 1. Přibližně určit charakteristiky každé realizace náhodného procesu. 2. Přibližné celkové charakteristiky lze získat zprůměrováním předchozích. Tato metoda zpracování je však poměrně složitá a vzniká otázka, či by nebylo možné pro stacionární náhodný proces zaměnit tento složitý přístup za mnohem jednodušší, který se zakládá na předpokladu, že střední hodnota nezávisí na čase a korelační funkce na začátku výpočtu. Kromě toho vzniká otázka, zda při zpracování pozorování stacionárního náhodného procesu je třeba disponovat několika jejich realizacemi. Protože náhodný proces je stacionární a homogenní v čase, je přirozené předpokládat, že jedna jediná realizace s dostatečnou délkou je postačujícím materiálem na získání charakteristik náhodného procesu. Při podrobnějším zkoumání této otázky se ukázalo, že existuje takováto možnost, ale ne pro všechny stacionární náhodné procesy. Tedy jestliže jediná realizace náhodného procesu pozorovaná v dostatečně dlouhém čase může být považovaná za určitého reprezentanta všech možných realizací, říkáme, že takovéto stacionární stochastické procesy mají ergodickou vlastnost. Jestliže určitý náhodný proces nemá tuto vlastnost ergodičnosti, i když je stacionární, potom jeho různé realizace, které se vyskytují s určitými pravděpodobnostmi, mají různý charakter průběhů. V tomto duchu, jako by šlo o realizace různých jednodušších stacionárních procesů, které mají svoje individuální charakteristiky. V některých případech na neergodičnost stacionárního procesu může působit už jen výskyt jediného náhodného sčítance (tj. náhodné proměnné nezávislé na čase). Příklad 2.4.1. Nechť {Y (t) = X(t) + Z, t R} je náhodný proces, kde {X(t), t R} je ergodický stacionární proces definovaný na pravděpo- dobnostním prostoru (, A, P) a Z náhodná veličina definovaná na témže pravděpodobnostním prostoru se střední hodnotou Z, rozptylem 2 Z a pro niž pro každé t R platí C(X(t), Z) = 0. Potom Y (t) = X + Z Y (t) = C(Y (s), Y (s + t)) = C(X(s) + Z, X(s + t) + Z) = = C(X(s), X(s + t)) X (t) + C(X(s + t), Z) =0 + C(Z, X(s + t)) =0 + C(Z, Z) 2 Z = X(t) + 2 Z. 2.4. Odhady středních hodnot a autokovariancí 36 Tedy náhodný proces {Y (t), t R} je stacionární proces, ale nemůžeme ho považovat za ergo- dický, neboť se dá očekávat, že každá jeho realizace se bude charakterem svého průběhu lišit od jiných - v závislosti od toho jakou hodnotu při dané realizaci nabyla náhodná veličina Z. Kovarianční funkce stacionárního procesu Y (t), t R se od kovarianční funkce stacionárního ergodického procesu {X(t), t R} liší o kladnou složku 2 Z. Takže pro t se hodnoty Y (t) nezmenšují k nule, ale od určitého času tm zůstávají konstantní (= 2 Z). Nyní budeme definovat ergodičnost stacionárních procesů přesněji matematicky v souvislosti s konstrukcí odhadů některých charakteristik stacionárních procesů. 2.4.1 Odhady střední hodnoty Nechť {Y (t), t R} je stochastický proces 2. řádu, který pozorujeme v časovém intervalu [0, T]. Nechť jeho konstantní střední hodnota je neznámá a je třeba ji odhadnout. Definice 2.4.2. Odhad střední hodnoty ^ stacionárního náhodného procesu {Y (t), t [0, T]} pomocí me- tody nejmenších čtverců (MNČ)je definován vztahem: T 0 (Y (t) - ^)2 dt = min T 0 (Y (t) - )2 dt. Poznámka 2.4.3. Stále budeme předpokládat, že integrály vystupující v jednotlivých vztazích existují a dají se v nich zaměnit pořadí integrování a střední hodnoty E. Snadno lze odvodit, že odhad střední hodnoty pomocí MNČ je roven ^ = 1 T T 0 Y (t)dt (2.21) neboť 0 = d d T 0 Y (t)2 - 2Y (t) + 2 dt = -2 T 0 Y (t)dt + 2 T 0 dt =T = 2T - 2 T 0 Y (t) dt . Věta 2.4.4. Odhad střední hodnoty pomocí metody nejmenších čtverců je nestranný a jeho střední kvadratická chyba je rovna MSE(^) = 2 T T 0 1 - u T Y (u) du. (2.22) Důkaz. Nestrannost: E^ = E 1 T T 0 Y (t)dt = 1 T T 0 EY (t) =(stac.) dt = 1 T T 0 dt =T = . 2.4. Odhady středních hodnot a autokovariancí 37 Střední kvadratická chyba v případě nestranného odhadu je rozptylem tohoto odhadu MSE(^) = E (^ - )2 = E (^ - E^)2 = D(^). Počítejme MSE(^) = E (^ - )2 = E 1 T T 0 Y (t) dt - 2 = E 1 T T 0 (Y (t) - ) dt 2 = 1 T2 E T 0 T 0 (Y (s) - )(Y (t) - ) ds dt = 1 T2 T 0 T 0 E [(Y (s) - )(Y (t) - )] Y (t-s)(stac.) ds dt = 1 T2 T 0 T 0 Y (t - s) ds dt Uvažujme transformaci u = t - s v = t s Jakobiánem |J| = 1. Protože s, t [0, T], pak platí -T u T 0 v = t T a tudíž u v = s + u T = u, tedy max{0, u} v min{T, T + u}. Tak dostaneme MSE(^) = 1 T2 T -T min{T,T +u} max{0,u} Y (u) dv du = 1 T2 0 -T Y (u) T +u 0 dv du + T 0 Y (u) T u dv du = 1 T2 0 -T Y (u)(T + u) du + T 0 Y (u)(T - u) du = 1 T2 T -T Y (u)(T - |u|) du = 2 T2 T 0 (T - u)Y (u) du = 2 T T 0 1 - u T Y (u) du = D(^) = D 1 T T 0 Y (t) dt , 2 Pro další studium ergodických procesů je vhodné vyslovit následující definici: Definice 2.4.5. Řekneme, že stacionární proces {Y (t), t R} 2. řádu je ergodický ve střední hodnotě, pokud platí lim T D 1 T T 0 Y (t) dt = 0. (2.23) 2.4. Odhady středních hodnot a autokovariancí 38 Věta 2.4.6. Nechť pro stacionární proces {Y (t), t R} 2. řádu s kovarianční funkcí Y (t) platí lim t 1 T T 0 1 - u T |Y (u)| du = 0. Potom je náhodný proces {Y (t), t R} ergodický ve střední hodnotě. Důkaz. Tvrzení věty plyne ze vztahů (2.22), (2.23) a nerovnosti T 0 1 - u T Y (u)du T 0 1 - u T |Y (u)| du. Důsledek 2.4.7. Nechť lim t Y (t) = 0. Pak stacionární proces s kovarianční funkcí Y (t) je ergodický ve střední hodnotě. Důkaz. Jestliže lim t Y (t) = 0, pak také lim t |Y (t)| = 0. Pak pro libovolně malé > 0 existují dostatečné velká T, T0 R (T0 < T) taková, že pro každé t > T0, platí |Y (t)| < . Pak lim T D 1 T T 0 Y (t) dt = lim T 2 T T 0 1 - u T Y (u) du lim T 2 T T 0 1 - u T |Y (u)| du lim T 2 T T 0 |Y (u)| du = lim T 2 T T0 0 |Y (u)| Y (0) du + T T0 |Y (u)| < du lim T 2 T0 T Y (0) + 1 - T0 T = 0 ergodicita ve střední hodnotě. 2 Poznamenejme, že jestliže platí lim t Y (t) = 0, pak také pro autokorelační funkci platí lim t Y (t) = lim t Y (t) Y (0) = 0, což znamená, že síla lineárních vazeb mezi jednotlivými náhodnými veličinami, které tvoří daný stacionární proces {Y (t), t R}, jakmile se tyto od sebe neustále vzdalují, postupně slábne, tj. jejich korelační koeficient 0. 2.4. Odhady středních hodnot a autokovariancí 39 DISKRÉTNÍ NÁHODNÉ PROCESY Při pozorování stacionárních procesů {Y (t), t R} druhého řadu se spojitým časem nejčastěji pozorujeme jen určitou jejich konečnou diskrétní část, tj. pro n N v diskrétních časových okamžicích t1, . . . , tn R pozorujeme jen náhodný vektor Y = (Yt1 , . . . , Ytn ) = (Y1, . . . , Yn) , který nazýváme diskrétním pozorováním náhodného procesu {Y (t), t R} (anebo dis- kretizací náhodného procesu {Y (t), t R} se spojitým časem), kde jsme položili ti = i, i = 1, . . . , n. Pak lze snadno ukázat, že obdobným diskrétním ekvivalentem odhadu střední hodnoty je odhad Y = 1 T n i=1 Yti T n R ti+t/2 ti-t/2 Y (t)dt = 1 n n t=1 Yt kde t = T n . 2.4.2 Odhady autokovarianční a autokorelační funkce Odhad kovarianční funkce lze analogicky jako v případě střední hodnoty nalézt ve tvaru ^Y () = 1 T - T - 0 [(Y (t) - ^) (Y (t + ) - ^)] dt. Podobně jak jsme výše definovali ergodičnost ve střední hodnotě pro stacionární proces {Y (t), t R}, můžeme definovat i jeho ergodičnost v rozptylu, pokud platí lim T D 1 T T 0 (Y (t) - )2 dt = 0 a jeho ergodičnost v kovarianční funkci, jestliže platí lim T D 1 T T - 0 (Y ( + t) - ) (Y (t) - ) dt = 0. Snadno lze ukázat, že obdobnými diskrétními ekvivalenty jsou následující odhady Odhad autokovarianční funkce ck = 1 n - k n-k t=1 (Yt - Y )(Yt+k - Y ) pro k = 0, 1, . . . , n - 1. Odhad autokorelační funkce ACF rk = ck c0 pro k = 0, 1, . . . , n - 1. Aby tyto odhady měly praktický význam, požaduje se obvykle n > 50 a k < n 4 , neboť odhady {ck}n-1 k=0, resp. {rk}n-1 k=0, nejsou lineárně nezávislé a s rostoucím k roste i jejich rozptyl. Kapitola 3 Analýza časových řad V dalších budeme u náhodného procesu {Xt, t T} indexovou množinu T interpretovat jako čas a pokud T = Z, budeme tento proces nazývat časovou řadou. 3.1 Základní přístupy k analýze časových řad V analýze časových řad se setkáváme s následujícími základními přístupy: 1. V časové doméně: a) klasická dekompozice časových řad je založena na regresní analýze; b) neoklasická dekompozice časových řad, tzv. Box-Jenkinsonova metodologie, je zalo- žena na korelační analýze; 2. Ve spektrální doméně: je spektrální analýza časových řad založena na Fourierově ana- lýze. 3.1.1 Klasická dekompozice Klasická dekompozice časových řad vychází z předpokladu, že náhodný proces, který generuje časovou řadu, je závislý pouze na čase. Snaží se pak časovou řadu rozdělit na ˇ deterministikou a ˇ náhodnou složku. Deterministickou složku dále rozkládá na trend a sezónní složku. Jednotlivé složky ˇ Trt trend odráží dlouhodobé změny v chování časových řad; ˇ Szt sezónní složka popisuje periodické změny (kratšího rázu) např. vliv střídání ročních období; ˇ t náhodné fluktuace modelují drobné a v jednotlivostech nepostižitelné příčiny ko- lísání časových řad. 40 3.2. Lineární dynamické modely 41 Při klasické dekompozici časových řad se používají především ˇ aditivní modely: Yt = Trt + Szt + t ˇ multipikativní modely Yt = Trt Szt t, které se transformují logaritmováním na aditivní modely. Klíčovým nástrojem klasické dekompozice časových řad je regresní analýza, kde jednotlivá pozorování se obvykle berou jako navzájem nekorelovaná, tj. {t, t Z} se chápe jako bílý šum. 3.1.2 Box-Jenkinsonova metodologie Box-Jenkinsonova metodologie na rozdíl od klasické dekompozice předpokládá, že všechny složky časové řady, tj. trend i cyklická složka, mají náhodný charakter. Těžištěm je pak korelační analýza a cílem je vyšetřit vzájemnou závislost jednotlivých prvků řady s různým zpožděním a závislost na různě zpožděném (náhodném) vstupu. 3.1.3 Spektrální analýza časových řad Spektrální analýza časových řad (Fourierovská analýza) považuje časovou řadu za "směs" si- nusovek a kosinovek o rozličných amplitudách a frekvencích. Tato koncepce pak umožní provést explicitní popis periodického chování časové řady a pře- devším vystopovat ty významné složky periodicity, které se podílejí na věcných vlastnostech zkoumaného procesu. V této koncepci není stěžejním faktorem časová proměnná, ale přávě faktor frekvenční. Spektrální analýza je také vhodná při srovnávání chování několika řad, neboť umožňuje nalézt případné časové zpoždění mezi dvěma řadami, ale také dovolí porovnat řady v rámci jednotlivých frekvencí. 3.2 Lineární dynamické modely V praxi se často setkáváme s tím, že hodnoty určité časové řady evidentně nejsou jen funkcí času, či předchozích pozorování, ale jsou vysvětlovány pomocí dalších časových řad, kterým říkáme faktorové časové řady a mluvíme o tzv. příčinných (kauzálních, faktorových) modelech. Modely se konstruují na základě teoretických předpokladů. Při jejich studiu se používají ˇ odhadovací metody, které jsou většinou modifikacemi základní metody nejmenších čtverců a lineární regresní analýzy; ˇ Box-Jenkinsonovy modely, v nichž kromě popisované řady a bílého šumu vystupují další vysvětlující časové řady. Příklad: Ct = + Ct-1 + Xt + Pt + t kde {Ct, t Z} časová řada výdajů obyvatelstva na nákup {Xt, t Z} časová řada disponibilních peněžních příjmů {Pt, t Z} časová řada cenový index spotřebního zboží {t, t Z} bílý šum , , , neznámé parametry 3.3. Klasická dekompozice časových řad 42 3.3 Klasická dekompozice časových řad Dekompozicí časové řady rozumíme rozklad časové řady na deterministickou a náhodnou složku, která má v případě aditivního modelu tvar Yt = Trt + Szt + t multiplikativního modelu Yt = Trt Szt t Jednotlivé složky Trt, Szt trend a sezónní složka mají deterministický charakter t náhodné fluktuace mají stochastický charakter, přičemž {t, t Z} je bílý šum, pro který platí Et = 0, C(t, s) = Ets = 0 s = t, Dt = 2 s = t. někdy se navíc předpokládá jeho normalita, tj. t N(0, 2 ). 3.3.1 Obecné lineární regresní modely Mějme regresní model plné hodnosti: Y = X + h(X) = h(X X) = p + 1 = k n > p + 1 Ln(0, 2 In) kde vektor závisle proměnných Y = (Y1, . . . , Yn) matice plánu X = (xij) i = 1, . . . , n; j = 0, . . . , p vektor chyb = (1, . . . , n) E = 0; D = 2 In; Tento model se také nazývá regresní model plné hodnosti s pevným plánem, neboť regresory xij (i = 1, . . . , n, j = 1, . . . , k) jsou nenáhodné, tj. pevně dané. Podmínka D = 2 In znamená, že náhodné veličiny Y1, . . . , Yn mají různé střední hod- noty (které jsou známou funkcí regresorů) a stejné rozptyly - mluvíme o homogenitě roz- ptylu. Odhad neznámých parametrů provedený metodou nejmenších čtverců je řešením normálních rovnic X X = X Y a platí: = (X X) -1 X Y Označme ˇ Y = X^ = X (X X) -1 X H Y = HY ˇ ^ = Y - Y = (I - H M )Y = MY = M(X + )) = MX =0 + M = (I - H) ˇ s2 = SSE n-p-1 = 1 n-p-1 (Y -Y) (Y -Y) = 1 n-p-1 ^ ^ = 1 n-p-1 Y (I-H)Y = 1 n-p-1 (I-H) 3.3. Klasická dekompozice časových řad 43 Pak platí: ˇ E = , ˇ Es2 = E(SSE) n-p-1 = 2 , tj. s2 je nestranným odhadem rozptylu, ˇ D = 2 (X X)-1 . Platí-li navíc Nn(0, 2 In) , pak ˇ Y Nn(X, 2 In) ˇ Nn(O, 2 (I - H)) ˇ Np+1(, 2 (X X)-1 ) ˇ SSE 2 2 (n - p - 1) ˇ a s2 jsou stochasticky nezávislé ˇ Tj = ^j-j s2vjj t(n - p - 1), kde (X X)-1 = (vij)i,j=0,...,p ˇ F = 1 qs2 (^2 - 2) W-1 (^2 - 2) F(q, n - p - 1), kde (X X)-1 = V U U W , = 1 2 ; ^ = ^1 ^2 a h(W) = q ˇ T = c ^-c s2c (X X)-1c t(n - p - 1), kde c = (c0, c1, . . . , cp) a E(c ) = c ˇ Označme i-tý řádek matice plánu X jako xi = (xi0, . . . , xip), pak Yi = xi + i N(xi, 2 ), ^Yi = xi N(xi, 2 xi(X X)-1 xi) Yi - ^Yi = xi( - ) + i N(0, 2 (1 + xi(X X)-1 xi)). Intervaly spolehlivosti pro parametry j j = 0, . . . , p j - t1- 2 (n-p-1)s vjj, j + t1- 2 (n-p-1)s vjj pro střední hodnotu predikce " xi b - t1- 2 (n-p-1)s p xi(X X)-1xi, xi b + t1- 2 (n-p-1)s p xi(X X)-1xi " E ^Yi =Exi b =xi pro predikci ^Yi = xi b " xi b - t1- 2 (n-p-1)s p 1+xi(X X)-1xi, xi b + t1- 2 (n-p-1)s p 1+xi(X X)-1xi " i = 1, . . . , n kde t1- 2 (n-p-1) je 1 - 2 kvantil Studentova rozdělení o n-p-1 stupních volnosti 3.3. Klasická dekompozice časových řad 44 3.3.2 Rozšířený lineární regresní model a vážená metoda nejmenších čtverců Následující věta ukazuje, jakým způsobem lze lineární regresní model rozšířit i na případ, kdy rozptyl není homogenní. Věta 3.3.1. Y = X + Ln(0, 2 V) V > 0 h(X) = k ^ = (X V-1 X)-1 X V-1 Y. Důkaz. Jelikož jsme předpokládali, že V > 0, tj. V je pozitivně definitní, takže existuje V- 1 2 , která je symetrická a regulární. Proto h(V- 1 2 X) = h(X) = k = h(X V-1 X) = h(X V- 1 2 V- 1 2 X) takže X V-1 X je regulární. Položme Z = V- 1 2 Y, F = V- 1 2 X, = V- 1 2 . Pak z Y = X + plyne, že V- 1 2 Y = V- 1 2 X + V- 1 2 , tj. Z = F + . Pak E = EV- 1 2 = V- 1 2 E =0 = 0 a D = D(V- 1 2 ) = 2 V- 1 2 VV- 1 2 = 2 V- 1 2 V 1 2 V 1 2 V- 1 2 = 2 In a tento model již splňuje předpoklady klasického regresního modelu, ve kterém odhad vektoru neznámých parametrů metodou nejmenších čtverců je roven ^ = (F F)-1 F Z = (X V- 1 2 V- 1 2 X)-1 X V- 1 2 V- 1 2 Y = (X V-1 X)-1 X V-1 Y. Poznámka 3.3.2. Nejčastěji se matice V uvažuje ve tvaru V = diag{v1, . . . , vn}, tj. jde o diagonální matici. Položíme-li W = V-1 = diag{ 1 v1 , . . . , 1 vn } = diag{w1, . . . , wn}, přičemž prvky w1, . . . , wn se nazývají váhami (tedy čím je rozptyl větší, tím je váha pozorování menší). Pak odhad neznámých parametrů metodou nejmenších čtverců: ^ = (X WX)-1 X WY se nazývá vážená metoda nejmenších čtverců. 3.3. Klasická dekompozice časových řad 45 3.3.3 Modelování trendu Klasická dekompozice časových řad deterministickou složku trend modeluje pomocí regres- ních modelů. Do trendu se obvykle zahrnují i cyklické složky s dlouhou periodou. Ke stanovení trendů lze v závislosti na jejich typu a charakteru náhodných fluktuací v dané časové řadě použít různé přístupy: ˇ parametrický přístup, který předpokládá určitý typ rozdělení, obvykle normální roz- dělení bílého šumu; ˇ neparametrický přístup, kam patří různé metody založené na jádrových odhadech, vyhlazovacích splajnech, waveletech apod. Regresní modely můžeme dále rozdělit na ˇ modely globálního trendu a ˇ modely postupného nebo lokálního trendu. MODELY GLOBÁLNÍHO TRENDU Mějme časovou řadu {Yt, t T} a n pozorování této řady v časových okamžicích t1 < t2 < ... < tn. Označme: Yi = Yti . Předpokládejme, že pozorování Yi vyhovují modelu: Yi = f(ti) + i i = 1, . . . , n (3.1) kde f(t) je neznámá trendová funkce uvažované řady vybraná z předem dané třídy funkcí i jsou náhodné fluktuace nevykazující systematickou odchylku, tj. Ei = 0, s konstantním rozptylem Di = 2 a nekolerované, t.j. C(i, j) = 0 i = j Regresní model trendu Parametrický přístup k odhadu trendové funkce vede na regresní modely trendu. Předpoklá- dáme, že trendová funkce f(t) patří do třídy funkcí, které jsou určeny konečným počtem para- metrů, přičemž f(t) = 0 + 11(t) + + pp(t) kde 0, 1, . . . , p neznámé parametry 0, 1, . . . , p známé funkce Nejprve uvažujme lineární regresní model: Y = X + h(X) = h(X X) = k = p + 1 n > p + 1 Nn(0, 2 In) (3.2) přičemž vektor závisle proměnných Y = (Y1, . . . , Yn) matice plánu X = (xij) = (j(ti)) i = 1, . . . , n; j = 0, . . . , p; 0(ti) 1 vektor chyb = (1, . . . , n) E = 0; D = 2 In; 3.3. Klasická dekompozice časových řad 46 Odhad neznámých parametrů provedený metodou nejmenších čtverců je řešením normálních rovnic X X = X Y a platí: ^ = (X X) -1 X Y (3.3) Označme reziduální rozptyl s2 = s2 k = SSE n - k = 1 n - p - 1 (Y - X) (Y - X) kde k = p + 1 (3.4) Z velkého okruhu trendových funkcí, které vedou k lineárnímu regresnímu modelu, se zamě- říme na ˇ polynomický trend: f(t) = 0 + 1t + + ptp ˇ periodický trend: f(t) = + p j=1(jcosjt + jsinjt) V případě polynomického trendu, matice plánu je tvaru X = 1 t1 t2 1 tp 1 ... ... ... ... ... 1 tn t2 n tp n . Kromě neznámých parametrů = (0, . . . , p) zbývá určit vhodný stupeň polynomu p . Pro odhad stupně polynomu se nabízí 2 intuitivní metody (1) ,,od nejnižšího stupně k nejvyššímu : začneme se stupněm p = 0 , postupně stupeň zvyšujeme a testujeme hypotézu H0 : p = 0 proti alternativě H1 : p = 0 pomocí statistiky (viz Anděl) Tp = ^p - p s2 kvpp t(n - p - 1), kde (X X) -1 = (vij)p i,j=0. Jestliže H0 zamítneme zvyšujeme stupeň polynomu. (2) ,,od maximálního stupně dolů : zvolme p = pmax . Testujeme opět H0 : p = 0 proti alternativě H1 : p = 0 pomocí Tp. Jestliže H0 nezamítneme snižujeme stupeň polynomu. Obě metody nedávají uspokojivé výsledky (viz Anderson(1971)). Penalizační metoda odhadu počtu regresních koeficientů Předpokládejme, že k0 je skutečný počet regresních parametrů. Lze ukázat, že platí E(s2 k) > 2 pro k < k0 E(s2 k) = 2 k k0 Zůstává problém, jak z grafu hodnot s2 k určit právě tu hodnotu k0, od níž počínaje již graf dostává vodorovný charakter. Tento problém se řeší zavedením tzv. penalizační funkce a např. Anděl navrhuje místo hodnot s2 k použít její modifikaci Ak = s2 k(1 + kwn) . 3.3. Klasická dekompozice časových řad 47 Penalizační funkce wn ˇ nesmí být příliš velká - aby nezkreslila klesající charakter s2 k pro k < k0; ˇ nesmí být příliš malá - aby z hodnot s2 k oscilujících kolem 2 vytvořila pro k 0 rostoucí posloupnost; Za odhad ^k se bere hodnota k {0, 1, . . . , kmax}, pro kterou Ak nabývá svého minima. Konstanta kmax je maximální počet parametrů, které jsme ochotni uvažovat a o němž jsme si jisti, že splňuje podmínku k0 kmax. Za dosti obecných podmínek týkajících se rozumné volby hodnot ti lze ukázat (Geweke a Meese(1981), Anděl a kol.(1981)), že pokud wn > 0 wn --- n 0 nwn --- n ^k k0 podle pravděpodobnosti. V praxi se osvědčilo volit wn = 1 4 n , tj. Ak = s2 k 1 + k 4 n . Další kriteria pro určení počtu regresních koeficientů Akaikeovo informační kritérium (1972) AICk = ln s2 k + 2k n nadhodnocuje k0 Swarz (1978) a Rissanen (1978) SRk = ln s2 k + k ln n n Hannan a Quinn (1979) HQk = ln s2 k + 2kc ln ln n n c > 1; obvykle c = 2 nebo 3. 50 100 150 200 250 300 7 7.5 8 8.5 9 9.5 10 10.5 BOX JENKINS SERIES D : CHEMICAL PROCESS VISCOSITY READINGS dgr=2 dgr=4 Obrázek 3.1: Vstupní data spolu s polynomickým trendem různých řádů. 3.3. Klasická dekompozice časových řad 48 Obrázek 3.2: Intervaly spolehlivosti pro pa- rametry i, . . . , k pro jednotlivé stupně poly- nomů. (Tečkovaná čára značí polohu nuly) 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 r2 =0.1947 p=0.0000 r 2 =0.3780 p=0.0000 r 2 =0.3870 p=0.0000 r2 =0.3945 p=0.0000 r2 =0.3947 p=0.0000 r 2 =0.3952 p=0.0000 r2 =0.3988 p=0.0000 r2 =0.4016 p=0.0000 coef. i degreeofpolynom SERIES C CHEMICAL PROCESS TEMPERATURE READINGS: EVERY MINUTE Obrázek 3.3: Penalizační kritéria pro odhad stupně polynomu. 1 2 3 4 5 6 7 8 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.3 Sk (Mean Square Error) 1 2 3 4 5 6 7 8 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Ak 1 2 3 4 5 6 7 8 -1.5 -1.45 -1.4 -1.35 -1.3 -1.25 -1.2 -1.15 AICk 1 2 3 4 5 6 7 8 -1.5 -1.45 -1.4 -1.35 -1.3 -1.25 -1.2 SR k 1 2 3 4 5 6 7 8 -1.5 -1.45 -1.4 -1.35 -1.3 -1.25 -1.2 HQ k (c=2) SERIES C CHEMICAL PROCESS TEMPERATURE READINGS: EVERY MINUTE 1 2 3 4 5 6 7 8 -1.5 -1.45 -1.4 -1.35 -1.3 -1.25 -1.2 -1.15 HQ k (c=3) 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=1; period=24.50 12.25 16.33 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=2; period=24.50 12.25 16.33 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=3; period=24.50 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=4; period=12.25 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=5; period=12.25 Polynomial and trigonometric trends 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895 0 5 10 15 20 25 30 35 40 Lepus-Olomouc dgr of pol.trend=6; period=12.25 Obrázek 3.4: Ukázka různých typů regresních modelů. 3.3. Klasická dekompozice časových řad 49 PERIODICKÝ TREND Je-li f(t) periodická funkce s periodou T, pak frekvencí rozumíme veličinu = 2 T . Uvažujme model: Yi = f(ti) + i Ei = 0; Di = 2 ; C(i, j) = 0; i = j; i, j = 1, . . . , n kde (a) f(ti) = + p j=1(jcosjti + jsinjti) nebo (b) f(ti) = + p j=1 jcos(jti + j) j = 2 j + 2 j j = arctan j j Jde o nelineární regresní model vzhledem k (3p + 1) neznámých parametrů: (a) 1, . . . , p 1, . . . , p 1, . . . , p (b) 1, . . . , p 1, . . . , p 1, . . . , p Odhad vektoru neznámých parametrů pomocí metody nejmenších čtverců minimalizuje vý- raz (a) S(, 1, . . . , p, 1, . . . , p, 1, . . . , p) = n i=1 (Yi - f(ti))2 (b) S(, 1, . . . , p, 1, . . . , p, 1, . . . , p) = n i=1 (Yi - f(ti))2 Numericky lze systém nelineárních rovnic řešit např. pomocí Gauss-Newtonovy metody. Lineární model pro známé frekvence Situace se zjednoduší, pokud frekvence 1, . . . , p jsou známé. Pak model (a) je lineární a matice plánu je tvaru: Xn×(2p+1) = 1 cos 1t1 sin 1t1 cos pt1 sin pt1 ... ... ... ... ... ... 1 cos 1tn sin 1tn cos ptn sin ptn . a X X(2p+1)×(2p+1) = n 0 0 0 n 2 ... ... ... 0 0 0 n 2 ^ = 1 n n i=1 Yi ^j = 2 n n i=1 Yi cos jti j = 1, . . . , p ^j = 2 n n i=1 Yi sin jti Neznámé parametry modelu (b) získáme ze vztahů ^j = ^2 j + ^2 j j = 1, . . . , p ^j = arctan ^j ^j Pokud časová řada vykazuje (po odečtení např. lineárního trendu) přibližně periodické cho- vání, je třeba rozhodnout, které frekvence se na tvorbě periodického trendu výrazně uplatňují. Pro nalezení významných period je výhodné užít metod spektrální analýzy časových řad. 3.3. Klasická dekompozice časových řad 50 MODELY LOKÁLNÍHO POSTUPNÉHO TRENDU Hlavní myšlenka lokální (vážené) metody nejmenších čtverců spočívá v tom, že prove- deme odhad trendu Trt polynomem na lokálním intervalu [t - s, t + s] na rozdíl od klasické (vážené) metody nejmenších čtverců, kdy trend odhadujeme polynomem na celém intervalu možných hodnot parametru t, který označíme [T1, T2]. Parametr s > 0 se nazývá šířka vyhlazovacího okénka. Interval [t-s, t+s] vyhlazovací okénko. I když vyhlazovací funkce, se kterou pracujeme, není polynomická funkce, může být za před- pokladu, že je lokálně hladká (tj. existují její spojité derivace až do nějakého vhodně zvoleného řádu), lokálně rozvedena do Taylorovy řady kolem bodu t. Proto může být dobře aproxi- mována lokálním polynomem, což lze provést metodou nejmenších čtverců, případně váženou metodou nejmenších čtverců. Popsaná lokální (vážená) metoda nejmenších čtverců se někdy též nazývá klouzavá po- lynomická metoda, protože kolem bodu t, v němž má být trend odhadnut, je umístěno vyhlazovací okénko [t - s, t + s] a odhad trendu Trt se ,,pohybuje spolu s t. Zvolme tento přístup: uvnitř vyhlazovacího okénka [t-s, t+s] aproximujme neznámý trend polynomem stupně m (x) = m j=0 j(t) (x - t)j . Koeficienty j(t) (j = 0, 1, . . . , m) uvádíme jakožto funkci bodu t, (který je středem okénka [t - s, t + s]), abychom zdůraznili, že tyto koeficienty budou pro každé t jiné. Neznámé koeficienty j(t) polynomu (x) odhadneme (váženou) metodou nejmen- ších čtverců, kde matice plánu X je tvořena prvky xij = (ti - t)j , přičemž j = 0, 1, . . . , m a index i nabývá pouze těch hodnot, pro které platí |ti - t| s. Je zřejmé, že platí (t) = 0(t) . Volba šířky vyhlazovacího okénka S rostoucím s pracujeme s větším počtem pozorování ve vyhlazovacím okénku [t - s, t + s], proto bude klesat rozptyl odhadu trendu, což však bude mít za následek nárůst jeho vychýlení od skutečné hodnoty. Vychýlení odhadu záleží na derivaci trendové funkce a projevuje se tak, že odhad Trt má tendenci podhodnocovat velikost lokálních extrémů trendové funkce, mluvíme o přehlazení. Pokud naopak budeme použivat úzké vyhlazovací okénko, odhad trendu bude méně vychýlen, ale na úkor velké variability odhadu. V tomto případě mluvíme o podhlazení trendové funkce. 3.3. Klasická dekompozice časových řad 51 V dalším budeme předpokládat, že posloupnost časových okamžiků t1, t2, . . . , tn - t1 t2 - ti ti+1 - tn-1 tn je ekvidistantní, tj. položíme-li pro i = 1, . . . , n - 1 = ti+1 - ti, pak ti = t1 + (i - 1) i = ti - t1 + 1 a položíme-li pro i = 1, . . . , n t i = ti - t1 + 1, můžeme bez újmy na obecnosti uvažovat pouze o časových řadách, pro něž platí ti = i i = 1, . . . , n. KLOUZAVÉ PRŮMĚRY Pokud zvolíme symetrické vyhlazovací okénko kolem bodu t - r = 2s + 1 členů t - s t t + s tak, že obsahuje úsek časové řady s lichým počtem členů r a položíme s = r - 1 2 , pak (pro jednoduchost místo j(t) pišme j) (x) = m j=0 j(x - t)j pro x [t - s, t + s]. Zaveďme substituci = x - t. Potom (x) = (t + ) = m j=0 jj [-s, s] (t) = 0 . 3.3. Klasická dekompozice časových řad 52 Máme tedy model Yt+ = m j=0 jj + t+ [-s, s]; Et+ = 0; Dt+ = 2 ; C(t+ , t+ ) = 0; = . Máme-li k dispozici časovou řadu délky n, pak pro každé t = s + 1, . . . , n - s zvlášť do- stáváme klasický lineární regresní model Y(t) = X + ; E = 0, D = 2 I2s+1 tj. Y(t) = Yt-s ... Yt ... Yt+s = 1 -s (-s)2 . . . (-s)m ... ... ... . . . ... 1 0 02 . . . 0m ... ... ... . . . ... 1 s s2 . . . sm 0 ... m + t-s ... t ... t+s = x1 ... xs+1 ... x2s+1 0 ... m + t-s ... t ... t+s Odhad neznámých parametrů provedený metodou nejmenších čtverců je řešením normálních rovnic X X = X Y(t) a platí ^(t) = (X X) -1 X Y(t) Y(t) = X (X X) -1 X projekční matice H Y(t) = HY(t), přičemž v tomto případě X X = 2s + 1 0 s =-s 2 0 s =-s 4 s =-s m 0 s =-s 2 0 s =-s 4 0 s =-s 6 s =-s 2 0 s =-s 4 0 s =-s 6 0 0 s =-s 4 0 s =-s 6 0 s =-s 8 s =-s 4 0 s =-s 6 0 s =-s 8 0 ... ... ... ... ... ... ... s =-s m 0 s =-s m+2 0 s =-s 2m , X Y(t) = s =-s Yt+ s =-s Yt+ ... s =-s m Yt+ a tudíž odhady ^Yt = xs+1(t) = (1, 0, . . . , 0) ^0 ^1 ... ^m = ^0 . Označme H = h1 ... hs+1 ... h2s+1 . Prvky projekční matice se nazývají váhy. Pak odhad ^Yt = hs+1Y(t) . 3.3. Klasická dekompozice časových řad 53 Pro prvních s hodnot, tj. pro t = 1, . . . , s vytváříme společný lineární regresní model: Y(F ) = Y1 ... Ys Ys+1 Ys+2 ... Y2s+1 = x1 ... xs xs+1 xs+2 ... x2s+1 0 1 ... m + 1 ... s s+1 s+2 ... 2s+1 Y(F ) = X (X X)-1 X Y(F ) = HY(F ) ^Y1 = h1Y(F ) ... ^Ys = hsY(F ) Obdobně pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n vytváříme společný lineární regresní model: Y(L) = Yn-2s ... Yn-s-1 Yn-s Yn-s+1 ... Yn = x1 ... xs xs+1 xs+2 ... x2s+1 0 1 ... m + n-2s ... n-s-1 n-s n-s+1 ... n Y(L) = X (X X)-1 X Y(L) = HY(L) ^Yn-s+1 = hs+2Y(L) ... ^Yn = h2s+1Y(L) Pak předchozí úvahy můžeme shrnout takto: pro prvních s hodnot, tj. pro t = 1, . . . , s dostáváme ^Yt = htY(F ) , pro tzv. "středové" hodnoty, tj. pro t = s + 1, . . . , n - s ^Yt = hs+1Y(t) a pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n ^Yt = ht-n+2s+1Y(L) . Pro ilustraci uvažujme příklad 1: m = 2, s = 2, r = 2s + 1 = 5 Matice plánu X = x1 x2 x3 x4 x5 = 1 -2 4 1 -1 1 1 0 0 1 1 1 1 2 4 , informační matice X X = 5 0 10 0 10 0 10 0 34 a projekční matice H = h1 h2 h3 h4 h5 = 31/35 9/35 -3/35 -1/7 3/35 9/35 13/35 12/35 6/35 -1/7 -3/35 12/35 17/35 12/35 -3/35 -1/7 6/35 12/35 13/35 9/35 3/35 -1/7 -3/35 9/35 31/35 váhy pro první bod (asymetrické váhy) váhy pro druhý bod (asymetrické váhy) váhy pro "středové" body (symetrické váhy) váhy pro předposlední bod (asymetrické váhy) váhy pro poslední bod (asymetrické váhy) 3.3. Klasická dekompozice časových řad 54 Ukažme ještě další příklad 2: m = 3, s = 2, r = 2s + 1 = 5 Matice plánu X = x1 x2 x3 x4 x5 = 1 -2 4 -8 1 -1 1 -1 1 0 0 0 1 1 1 1 1 2 4 8 , informační matice X X= 5 0 10 0 0 10 0 34 10 0 34 0 0 34 0 130 a projekční matice H = h1 h2 h3 h4 h5 = 69/70 2/35 -3/35 2/35 -1/70 2/35 27/35 12/35 -8/35 2/35 -3/35 12/35 17/35 12/35 -3/35 2/35 -8/35 12/35 27/35 2/35 -1/70 2/35 -3/35 2/35 69/70 váhy pro první bod (asymetrické váhy) váhy pro druhý bod (asymetrické váhy) váhy pro "středové" body (symetrické váhy) váhy pro předposlední bod (asymetrické váhy) váhy pro poslední bod (asymetrické váhy) Vidíme, že ˇ Součet vah v jednom řádku je roven jedné (jde o prvky projekční matice s jednotkovou normou). ˇ Středové váhy jsou symetrické kolem prostřední hodnoty. ˇ Je-li m sudé, pak "středové" váhy řádu m a m + 1 pro stejnou délku r = 2s + 1 jsou totožné. Jednoduché klouzavé pruměry Uvažujme nejprve: m = 0, r = 2s + 1 Spočtěme postupně X2s+1 = (1, . . . , 1) , X X = 2s + 1, X Y(t) = s =-s Yt+ , normální rovnice: (2s + 1)^0 = s =-s Yt+ ^0 = 1 2s + 1 s =-s Yt+ . Pro t = s + 1, . . . , n - s máme ^Yt = ^(t) = ^0 = 1 2s + 1 s =-s Yt+ Pro prvních s hodnot, tj. pro t = 1, . . . , s dostáváme ^Yt = 1 2s + 1 2s+1 =1 Y 3.3. Klasická dekompozice časových řad 55 a pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n jsou odhady ve tvary ^Yt = 1 2s + 1 2s+1 =1 Yn+1- . Uvažujme příklad: m = 1, s = 2, r = 2s + 1 = 5 Matice plánu X = x1 x2 x3 x4 x5 = 1 1 1 1 1 , informační matice X X = 5 a projekční matice H = h1 h2 h3 h4 h5 = 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 1/5 váhy pro první bod váhy pro druhý bod váhy pro "středové" body váhy pro předposlední bod váhy pro poslední bod Dále uvažujme: m = 1, r = 2s + 1 Spočtěme pomocný vzorec: s =-s 2 = 2 s =1 2 = 2 1 6 s(s + 1)(2s + 1) = s(s + 1)(2s + 1) 3 . Normální rovnice 2s + 1 0 0 s =-s 2 0 1 = s =-s Yt+ s =-s Yt+ ^0 = 1 2s+1 s =-s Yt+ ^1 = 3 s(s+1)(2s+1) s =-s Yt+ Pro t = s + 1, . . . , n - s máme ^Yt = ^(t) = ^0 = 1 2s + 1 s =-s Yt+ což je stejné jako v pro m = 0. Pro prvních s hodnot, tj. pro t = 1, . . . , s dostáváme ^F 0 = 1 2s + 1 2s+1 =1 Y , ^F 1 = 3 s(s + 1)(2s + 1) s =-s Ys++1 a ^Y1 = ^F 0 - s^F 1 , ^Y2 = ^F 0 - (s - 1)^F 1 , . . . , ^Ys = ^F 0 - ^F 1 . Pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n jsou odhady ve tvary ^L 0 = 1 2s + 1 2s+1 =1 Yn-+1, ^L 1 = 3 s(s + 1)(2s + 1) s =-s Yn-2s+-1 3.3. Klasická dekompozice časových řad 56 a ^Yn-s+1 = ^L 0 + ^L 1 , ^Yn-s+2 = ^L 0 + 2^L 1 , . . . , ^Yn = ^L 0 + s^L 1 . Uvažujme příklad: m = 1, s = 2, r = 2s + 1 = 5 Matice plánu X = x1 x2 x3 x4 x5 = 1 -2 1 -1 1 0 1 1 1 2 , informační matice X X = 5 0 0 10 a projekční matice H = h1 h2 h3 h4 h5 = 3/5 2/5 1/5 0 -1/5 2/5 3/10 1/5 1/10 0 1/5 1/5 1/5 1/5 1/5 0 1/10 1/5 3/10 2/5 -1/5 0 1/5 2/5 3/5 váhy pro první bod (asymetrické váhy) váhy pro druhý bod (asymetrické váhy) váhy pro "středové" body (symetrické váhy) váhy pro předposlední bod (asymetrické váhy) váhy pro poslední bod (asymetrické váhy) Centrované klouzavé průměry Často dochází k situaci, kdy chceme průměrovat hodnoty přes sudý počet období (např. pro vyrovnávání sezónních fluktuací), např. o délce 12 při naměřených měsíčních pozorováních nebo o délce 4 při čtvrtletních pozorováních. Klouzavý průměr délky 12 by sice vyrovnal z velké části sezónní fluktuace v řadě, přitom však např. aritmetický průměr lednové až prosincové hodnoty by padl mezi body červen a červenec. Proto se zavedl jiný postup: t - 6 t - 5 t - 4 t - 3 t - 2 t - 1 t t + 1 t + 2 t + 3 t + 4 t + 5 t + 6 - ? p1 = 1 12 (Yt-6 + + Yt+5) - 6 p2 = 1 12 (Yt-5 + + Yt+6) ˇ V 1. kroku se vypočtou p1 = 1 12 (Yt-6 + + Yt+5) a p2 = 1 12 (Yt-5 + + Yt+6) ˇ V 2. kroku p3 = 1 2 (p1 + p2) = 1 212 (Yt-6 + 2Yt-5 + 2Yt+5 + Yt+6) Obecně centrovaný jednoduchý klouzavý průměr délky 2s+1 je tvaru ^Yt = 1 4s (Yt-s + 2Yt-s+1 + 2Yt+s-1 + Yt+s) . Volba řádu klouzavých průměrů Budeme se snažit najít nějaké objektivní kritérium. Předpokládejme, že Yt = Trt + t = 0 + 1t + + mtm + t t WN(0, 2 ). 3.3. Klasická dekompozice časových řad 57 Zaveďme obecně diference: yt = yt - yt-1 2 yt = yt - yt-1 = yt - 2yt-1 + yt-2 ... k yt = yt - k 1 yt-1 + k 2 yt-2 - + - 1)k yt-k Polynom 0 + 1t + + mtm bude při každé následující diferenci snižovat svůj řád o jedničku. Konečně při diferenci řádu m + 1 se tento polynom vynuluje. Bílý šum vytváří při k -té diferenci k t = t - k 1 t-1 + k 2 t-2 - + (-1)k t-k Spočtěme střední hodnotu a rozptyl této k -té diference: Ek t = E t - k 1 t-1 + k 2 t-2 - + (-1)k t-k = 0 Dk t = D t - k 1 t-1 + k 2 t-2 - + (-1)k t-k = 2 1 + k 1 2 + k 2 2 + + 1 Výraz ve složených závorkách je zřejmě koeficient členu xk ve výrazu (1 + x)k (1 + x)k a je tedy s využitím binomické věty roven 2k k . Takže Dk t = 2k k 2 . Proveďme standardizaci k -té diference bílého šumu Uk,t = k t - 0 2k k . Položme S2 k = n t=k+1 U2 k,t. Pak ES2 k = (n - k)2 . Zavedeme-li pro k m + 1 veličinu Vk = n t=k+1 k Yt 2 2k k (n - k) , pak Vk je pro k m + 1 odhadem rozptylu bílého šumu. Počítejme V1, V2, . . . dokud nezazna- menáme, že tyto hodnoty začínají konvergovat k nějaké konstantě. 3.3. Klasická dekompozice časových řad 58 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Simulace sinusovky s bilym sumem 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=0 delky r=2s+1=7 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=0 delky r=2s+1=21 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=0 delky r=2s+1=31 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Simulace sinusovky s bilym sumem 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=1 delky r=2s+1=7 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=1 delky r=2s+1=21 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=1 delky r=2s+1=31 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Simulace sinusovky s bilym sumem 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=2 delky r=2s+1=7 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=2 delky r=2s+1=21 0 20 40 60 80 100 0.5 1 1.5 2 2.5 3 3.5 Klouzave prumery radu m=2 delky r=2s+1=31 Obrázek 3.5: Různá volba řádu a délky klouzavých průměrů pro simulovaná data 3.3. Klasická dekompozice časových řad 59 0 5 10 15 20 25 30 35 40 35 40 45 50 55 60 65 70 Daily morning temperature of a cow. V(k) rad diference 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0 delky r=2s+1=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0 delky r=2s+1=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0 delky r=2s+1=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1 delky r=2s+1=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1 delky r=2s+1=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1 delky r=2s+1=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2 delky r=2s+1=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2 delky r=2s+1=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2 delky r=2s+1=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3 delky r=2s+1=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3 delky r=2s+1=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3 delky r=2s+1=21 Obrázek 3.6: Různá volba řádu a délky klouzavých průměrů pro reálná data 3.3. Klasická dekompozice časových řad 60 EXPONENCIÁLNÍ VYROVNÁVÁNÍ Základy položili Holt (1957) a Brown (1961). Na rozdíl od klouzavých průměrů vychází z poly- nomiální lokální vážené metody nejmenších čtverců, kde váhy jednotlivých čtverců se směrem do minulosti exponenciálně snižují = odtud název metody. Uvažujeme-li vyhlazovací okno pouze směrem do minulosti, pak pro každé t máme následu- jící regresní model, přičemž = 0, 1, . . . Yt- = m j=0 (-)j j + t- , Et- = 0; Eqs = 0, q = s Dt- = - 2 ; (0, 1) . Použijeme-li maticový zápis, dostaneme Y = X + E = 0 D = 2 W-1 (0, 1) kde Y = Yt-0 Yt-1 Yt-2 ... X = 1 0 0 0 1 (-1)1 (-1)2 (-1)m 1 (-2)1 (-2)2 (-2)m ... ... ... ... a W = 0 0 0 0 1 0 0 0 2 ... ... ... ... = t t-1 t-2 ... . Odhad parametrů metodou nejmenších vážených čtverců (neboť rozptyly nejsou kon- stantní) je dán vzorcem: ^ = (X WX)-1 X WY kde X WX = =0 =0(-)1 =0(-)m =0(-)1 =0(-)2 =0(-)m+1 ... ... ... ... =0(-)m =0(-)m+1 =0(-)2m a X WY = =0 Yt- =0(-)1 Yt- ... =0(-)m Yt- . 3.3. Klasická dekompozice časových řad 61 JEDNODUCHÉ EXPONENCIÁLNÍ VYROVNÁVÁNÍ Vyjádříme-li explicitně normální rovnice X WX = X WY pro m = 0, přičemž použijeme označení ^0 = b0, dostaneme bo =0 = =0 Yt- Pomocné vztahy 3.3.3. Protože platí (0, 1), pak =0 = 1 1- . Protože ^Yt = b0, dostáváme ^Yt = b0 = (1 - ) =0 Yt- = (1 - )Yt + (1 - ) =1 Yt- subst.k=-1 = = (1 - )Yt + (1 - ) k=0 k+1 Yt-1-k = = (1 - )Yt + (1 - ) k=0 k Yt-1-k ^Yt-1 = (1 - )Yt + ^Yt-1 Dostáváme tedy jednoduchý rekurentní vzorec: ^Yt = (1 - )Yt + ^Yt-1 . Lze snadno řídit adaptivnost metody: 1. při malém metoda rychle reaguje na změny v charakteru dat, neboť se prosadí vliv prvního sčítance 2. při větším zesílí se vyrovnávací schopnost. Předpovídání (predikce) Chceme-li použít metodu jednoduchého exponenciálního vyrovnávání pro předpovídání, pak klademe ^Yt+ (t) = ^Yt pro libovolné > 0. Přitom ^Yt+ (t) označuje předpověď hodnoty Yt+ konstruované v čase t na základě hodnot Yt, Yt-1, . . .. Speciálně v koncovém bodě ^Yn+ (n) = ^Yn. 3.3. Klasická dekompozice časových řad 62 Vysvětlení rekurentního vzorce jiným způsobem Upravujme postupně ^Yt = (1 - )Yt + ^Yt-1 = (1 - )Yt + ^Yt-1 + ^Yt-1 - ^Yt-1 = ^Yt-1 + (1 - )(Yt - ^Yt-1) Podívejme se na odhad ^Yt-1, kterou hodnotu předpovídá, považujeme-li t-1 za poslední známou hodnotu, tj. ^Yt-1 = ^Yt(t - 1) a odtud ^Yt-1 je předpověď pro čas t zkonstruovaná v časovém okamžiku t - 1 (tj. z dostupných hodnot t - 1, t - 2, . . .). Pak ^Yt = ^Yt-1 + (1 - ) (Yt - ^Yt-1) chyba predikce = ^Yt-1 + (1 - ) Yt - ^Yt(t - 1) Jinými slovy: ˇ pro opravu předchozí vyrovnané hodnoty použijeme (jakmile dostaneme pozorování Yt) vhodně diskontovanou chybu předpovědi ^Yt(t - 1) konstruované v čase t - 1. Předpovědní intervaly na hladině významnosti Podle Bowerman,B.L., OConnell, R.T.: Time series and forecasting. North Scituate, Massa- chusetts, Duxbury Press (1979) mají tvar ^Yn+ (n) - u1- 2 d (n) , ^Yn+ (n) + u1- 2 d (n) kde (n) = 1 n n t=1 |Yt - ^Yt(t - 1)| a d = 1.25 pro libovolné > 0. Výhodou předpovědních intervalů konstruovaných v rámci exponenciálního vyrovnávání je jejich snadné přizpůsobení při dodání dalších pozorování časové řady. 3.3. Klasická dekompozice časových řad 63 Hledání optimálního parametru Na základě bohatých praktických zkušeností s danou metodou je možné se pro většinu apli- kací omezit při výběru na interval [0.7, 1). Výběr z tohoto intervalu lze například provést pomocí simulací: volíme postupně rovno 0.70, 0.72, 0.74, . . ., 0, 98. Pak vybereme tu hodnotu , která v dané řadě poskytuje nejlepší předpovědi. Poznámka 3.3.4. Někdy se však neklade požadavek, aby bylo z intervalu [0.7, 1). Jakmile však optimální nalezená hodnota leží vně tohoto intervalu, je to indikací toho, že s použitou metodou není něco v pořádku, tj. např. vhodnější by byla metoda dvojitého exponenciálního vyrovnávání. Při řešení konkrétního příkladu se setkáváme s těmito problémy: 1. chceme-li použít rekurentní vzorec, potřebujeme ^Y0. 2. potřebujeme určit . Proto se algoritmus rozpadá do dvou fází: 1. volba konstanty . (a) spočteme ^Y0 jako aritmetický průměr prvních n1 členů (doporučuje se volit n1 = 6 < n nebo n1 = n 2 ). Nebereme n1 = n z toho důvodu, že nejprve chceme simulačně určit hodnotu pro 0.70, 0.72, 0.74, . . ., 0, 98 a pokud bychom vzali arit- metický průměr ze všech hodnot, nemohli bychom posoudit kvalitu vyrovnávání v historických datech, neboť by se všechny účastnily na odhadu počátečního ^Y0. (b) Pro i [0.7, 1) se spočtou odhady ^Y1 = (1 - i)Y1 + i ^Y0 ... ^Yn = (1 - i)Yn + i ^Yn-1 předpovědi ^Yt(t - 1) = ^Yt-1 a parametr se určí ze vztahu = arg min i SSE = arg min i n t=1(Yt - ^Yt-1)2 2. Nejprve se určí ^Y0 = 1 n2 n2 i=1 Yi. Většinou se bere n2 = n. Pro zvolené v 1. fázi se spočtou odhady ^Y1, . . . , ^Yn, pak předpovědi ^Yn+1(n) = ^Yn = (1 - )Yn + i ^Yn-1. Pokud chceme vypočítat ^Yn+1, musíme si počkat, až dostaneme hodnotu Yn+1, jinak musíme vystačit se vzorcem ^Yn+ (n) = ^Yn pro > 0. 3.3. Klasická dekompozice časových řad 64 Adaptivní řídící proces - tj. změna během výpočtu Adaptivní řídící proces během výpočtu pomocí jistého indikátoru poruchy signalizuje, že vy- rovnávání přestává být pro příslušnou vyrovnávací konstantu adekvátní (někdy dokonce tuto hodnotu podle potřeby automaticky mění). Jako indikátor poruchy se často používá veličina I(, t) = Y (, t) D(, t) , kde označíme-li ej() = Yj - ^Yj(j - 1), pak Y (, t) = t j=1 ej() a D(, t) = 1 t t j=1 |ej()| . Zvedne-li se I(, t) nad jistou kontrolní mez K, je to signál ke změně hodnoty nebo dokonce přestává vyhovovat uvažovaný typ trendu (a je např. nutné přejít na dvojité exponenciální vyrovnávání). Mez K se obvykle volí mezi 4 až 6. Další metoda: Souběžně se provádí 3 procedury pro - 0.05 + 0.05. Pokud ˇ například D(, t) min (D( - 0.05, t), D( + 0.05, t)) , pokračuje se, ˇ pokud D( + 0.05, t) < D(, t) se změní na + 0.05 a pokračuje se pro + 0.05 + 0.1. 3.3. Klasická dekompozice časových řad 65 DVOJITÉ EXPONENCIÁLNÍ VYROVNÁVÁNÍ Vyjádříme-li explicitně normální rovnice X WX = X WY pro m = 1 , dostaneme =0 - =0 - =0 =0 2 0 1 = =0 Yt- - =0 Yt- . Pomocné vztahy 3.3.5. Pro (0, 1) platí =0 = 1 1 - =0 = (1 - )2 =0 2 = (1 + ) (1 - )3 Důkaz. Protože (0, 1), pak =0 = 1 1 - d d =0 = =0 -1 = 1 (1 - )2 d d =0 -1 = =0 ( - 1)-2 = 2 (1 - )3 Díky předchozím vztahům dostaneme =0 = =0 -1 = (1 - )2 =0 2 = =0 ( - 1) + =0 = 2 =0 ( - 1)-2 + (1 - )2 = 22 (1 - )3 + (1 - )2 = (1 + ) (1 - )3 . 2 3.3. Klasická dekompozice časových řad 66 Pokračujme v řešení systému normálních rovnic 1 1- - (1-)2 - (1-)2 (1+) (1-)3 0 1 = =0 Yt- - =0 Yt- . Pokud postupně první a druhou rovnici vynásobíme výrazem (1-), resp. -(1-)2 , dostaneme: 1 - 1- -(1+) 1- 0 1 = (1 - ) =0 Yt- (1 - )2 =0 Yt- . Pomocné vztahy 3.3.6. Jestliže definujeme tzv. vyrovnávací statistiky 1. a 2. řádu St = (1 - ) =0 Yt- a S [2] t = (1 - ) =0 St- , pak platí St = (1 - )Yt + St-1 S [2] t = (1 - )St + S [2] t-1 a také (1 - )2 =0 Yt- = S [2] t - (1 - )St. Důkaz. Upravujme postupně St = (1 - ) =0 Yt- = (1 - )Yt + (1 - ) =1 Yt- subst.k=-1 = = (1 - )Yt + (1 - ) k=0 k Yt-1-k St-1 = (1 - )Yt + St-1 S [2] t = (1 - ) =0 St- = (1 - )St + (1 - ) =1 St- subst.:k=-1 = = (1 - )St + (1 - ) k=0 k+1 St-1-k = (1 - )St + (1 - ) k=0 k St-1-k S [2] t-1 = (1 - )St + S [2] t-1 . 3.3. Klasická dekompozice časových řad 67 Všimněme si opět S [2] t . Protože St-k = (1 - ) =0 Yt-k- , pak S [2] t = (1 - ) k=0 k St-k = (1 - )2 k=0 k =0 Yt-k- = (1 - )2 k=0 k 0 Yt-k + 1 Yt-k-1 + 2 Yt-k-2 + = (1 - )2 k=0 k Yt-k první řada =(1-)St +(1 - )2 { 1 Yt-1 + 2 Yt-2 + 3 Yt-3 + + 2 Yt-2 + 3 Yt-3 + 4 Yt-4 + + 3 Yt-3 + 4 Yt-4 + 5 Yt-5 + + } = (1 - )St + (1 - )2 {Yt-1 + 22 Yt-2 + 33 Yt-3 + } = (1 - )St + (1 - )2 =0 Yt- A odtud plyne (1 - )2 =0 Yt- = S [2] t - (1 - )St. 2 Pokračujme opět v řešení systému normálních rovnic, přičemž použijeme označení b0 = ^0 a b1 = ^1. b0 - 1 - b1 = St b0 = St + 1 - b1 b0 - (1 + ) 1 - b1 = S [2] t - (1 - )St Dosaďme za b0 do druhé rovnice: St + 2 1 - b1 - (1 + ) 1 - b1 = S [2] t - St - St - 1 - b1 = S [2] t - St b1 = 1 - (St - S [2] t ) Dopočítejme b0 b0 = St + 1 - b1 = St + 1 - 1 - (St - S [2] t ) = 2St - S [2] t . Protože ^Yt = b0, dostáváme celkově ^Yt = 2St - S [2] t St = (1 - )Yt + St-1 S [2] t = (1 - )St + S [2] t-1 3.3. Klasická dekompozice časových řad 68 5 10 15 20 25 100 150 200 250 300 350 400 Jednoduche exponencialni vyrovnani =0.70 5 10 15 20 25 100 150 200 250 300 350 400 Dvojite exponencialni vyrovnani =0.84 5 10 15 20 25 100 150 200 250 300 350 400 Jednoduche exponencialni vyrovnani adaptivni procedura Mesicni exportni odbyt (v tisicich Kc) jisteho strojirenskeho podniku. 5 10 15 20 25 100 150 200 250 300 350 400 Dvojite exponencialni vyrovnani adaptivni procedura 2 4 6 8 10 12 14 16 0 2 4 6 8 x 10 4 Jednoduche exponencialni vyrovnani =0.70 2 4 6 8 10 12 14 16 0 2 4 6 8 x 10 4 Dvojite exponencialni vyrovnani =0.70 2 4 6 8 10 12 14 16 0 2 4 6 8 x 10 4 Jednoduche exponencialni vyrovnani adaptivni procedura Rocni vyvozy zahr. obchodu (v milionech Kcs) v CSSR v letech 1965-1980 2 4 6 8 10 12 14 16 0 2 4 6 8 x 10 4 Dvojite exponencialni vyrovnani adaptivni procedura 3.3. Klasická dekompozice časových řad 69 5 10 15 20 25 30 35 -200 0 200 400 600 800 1000 1200 Jednoduche exponencialni vyrovnani =0.70 5 10 15 20 25 30 35 -200 0 200 400 600 800 1000 1200 Dvojite exponencialni vyrovnani =0.82 5 10 15 20 25 30 35 -200 0 200 400 600 800 1000 1200 Jednoduche exponencialni vyrovnani adaptivni procedura Zaznamy okr.spravy o mesicnich poctech deti ve skolach v prirode(1980-1982). 5 10 15 20 25 30 35 -200 0 200 400 600 800 1000 1200 Dvojite exponencialni vyrovnani adaptivni procedura 5 10 15 20 25 30 26 27 28 29 30 31 32 Jednoduche exponencialni vyrovnani =0.94 5 10 15 20 25 30 26 27 28 29 30 31 32 Dvojite exponencialni vyrovnani =0.84 5 10 15 20 25 30 26 27 28 29 30 31 32 Jednoduche exponencialni vyrovnani adaptivni procedura Tydenni objemy zakazek (v tisicich Kc) prijatych jistou zakazkovou prodejnou. 5 10 15 20 25 30 26 27 28 29 30 31 32 Dvojite exponencialni vyrovnani adaptivni procedura 3.3. Klasická dekompozice časových řad 70 3.3.4 Sezónnost Sezónní složka je periodická složka, jejíž délka periody je kratší. Délku periody lze zpravidla rozpoznat už na základě intuitivní úvahy, neboť pravidelnost oscilací v sezónních řadách se velmi často připisuje působení sluneční soustavy na průběh těchto dějů. Vliv oběhu Země kolem Slunce se promítá buď přímo (klimatické podmínky), anebo ve zprostředkovaných, nepřímých souvislostech (např. ukazatele cestovního ruchu). Metoda malého trendu Uvažujme regresní model ve tvaru: Yt = Trt + Szt + t t WN(0, 2 ). Přeindexujme Y1, . . . , Yn na Yjk, j = 1, . . . , r r . . . počet sezón 1, . . . , n jk, k = 1, . . . , d d . . . délka sezóny. Předpokládejme, že trend je konstantní pro j-tou sezónu, tj. Trj = mj a rovněž sezónní hodnota je konstantní pro k-tou sezónní složku, tj. Szk = sk. Regresní model můžeme napsat ve tvaru Yjk = mj + sk + jk t WN(0, 2 ), j = 1, . . . , r; k = 1, . . . , d. Maticově, lze tento model rozepsat takto Y11 ... ... Y1d Y21 ... ... Y2d ... ... ... Yr1 ... ... Yrd = 1 0 0 | 1 0 0 ... ... ... ... ... ... ... | 0 1 ... ... ... ... ... ... ... ... ... | ... ... ... 0 1 0 0 | 0 0 1 0 1 0 0 | 1 0 0 ... ... ... ... ... ... ... | 0 1 ... ... ... ... ... ... ... ... ... | ... ... ... 0 0 1 0 0 | 0 0 1 0 0 0 0 1 0 | 1 0 0 ... ... ... ... ... ... ... | 0 1 ... ... ... ... ... ... ... ... ... | ... ... ... 0 0 0 0 1 0 | 0 0 1 0 0 0 0 1 | 1 0 0 ... ... ... ... ... ... ... | 0 1 ... ... ... ... ... ... ... ... ... | ... ... ... 0 0 0 0 0 1 | 0 0 1 m1 ... mr s1 ... sd + 11 ... ... 1d 21 ... ... 2d ... ... ... r1 ... ... rd Matice plánu však není plné hodnosti, neboť když sečteme prvních r sloupců, dostaneme vektor samých jedniček, což je rovno také součtu posledních d sloupců. Proto přidejme ještě jednu podmínku, a to d k=1 sk = 0. 3.3. Klasická dekompozice časových řad 71 Potom X X = d 0 0 | 1 1 0 d ... ... | ... ... ... ... ... ... ... 0 | ... ... ... ... 0 0 d | 1 1 1 1 | r 0 0 ... ... ... ... | 0 r ... ... ... ... ... ... | ... ... ... 0 1 1 | 0 0 r a X Y = Y1 ... Yr Y1 ... Yd , kde využíváme tzv. tečkové notace Yj = d i=1 Yji Yk = r i=1 Yik. Normální rovnice X X = X Y můžeme přepsat pro j =1, . . . , r k =1, . . . , d do tvaru dmj + d i=1 si =0 = Yj mj = 1 d Yj = Yj r i=1 mi + rsk = Yk rsk = Yk - r i=1 mi = r i=1 (Yik - mi) sk = 1 r r i=1 (Yik - mi) 3.3. Klasická dekompozice časových řad 72 Model: polynomický trend za celé období spolu se sezónností Uvažujme regresní model ve tvaru: Yt = Trt + Szt + t t WN(0, 2 ). Přeindexujme Y1, . . . , Yn na Yjk, j = 1, . . . , r r . . . počet sezón 1, . . . , n jk, k = 1, . . . , d d . . . délka sezóny. Předpokládejme, že trend je polynomický za celé období, tj. Trt = 0 + 1t + + mtm a sezónní hodnota je konstantní pro k-tou sezónní složku, tj. Szk = sk. Regresní model můžeme napsat ve tvaru: Yjk = 0 + 1t + + mtm + sk + jk j = 1, . . . , r k = 1, . . . , d t = (j - 1)d + k Matice plánu je pak tvaru X = 1 1 1m | 1 0 0 0 1 2 2m | 0 1 0 0 ... ... ... ... | 0 0 ... 0 ... ... ... ... | ... ... ... ... 0 1 d dm | 0 0 0 1 1 d + 1 (d + 1)m | 1 0 0 0 1 d + 2 (d + 2)m | 0 1 0 0 ... ... ... ... | 0 0 ... 0 ... ... ... ... | ... ... ... ... 0 1 2d (2d)m | 0 0 0 1 ... ... ... ... | 1 0 0 0 ... ... ... ... | 0 1 0 0 ... ... ... ... | 0 0 ... 0 ... ... ... ... | ... ... ... ... 0 ... ... ... ... | 0 0 0 1 1 (r - 1)d + 1 [(r - 1)d + 1]m | 1 0 0 0 1 (r - 1)d + 2 [(r - 1)d + 2]m | 0 1 0 0 ... ... ... ... | 0 0 ... 0 ... ... ... ... | ... ... ... ... 0 1 rd (rd)m | 0 0 0 1 Matice plánu však není plné hodnosti, neboť když sečteme posledních d sloupců, dostaneme vektor samých jedniček. Proto použijeme tvz. metodu horního rohu a položíme první sezónu rovnu nule s1 = 0. 3.3. Klasická dekompozice časových řad 73 Za těchto předpokladů lze model maticově napsat ve tvaru Y11 ... ... ... Y1d Y21 ... ... ... Y2d ... ... ... ... Yr1 ... ... ... Yrd = 1 1 1m | 0 0 0 1 2 2m | 1 0 0 ... ... ... ... | 0 ... ... 0 ... ... ... ... | ... ... ... 0 1 d dm | 0 0 1 1 d + 1 (d + 1)m | 0 0 0 1 d + 2 (d + 2)m | 1 0 0 ... ... ... ... | 0 ... ... 0 ... ... ... ... | ... ... ... 0 1 2d (2d)m | 0 0 1 ... ... ... ... | 0 0 0 ... ... ... ... | 1 0 0 ... ... ... ... | ... ... ... 0 ... ... ... ... | 0 ... ... 0 ... ... ... ... | 0 0 1 1 (r - 1)d + 1 [(r - 1)d + 1]m | 0 0 0 1 (r - 1)d + 2 [(r - 1)d + 2]m | 1 0 0 ... ... ... ... | 0 ... ... 0 ... ... ... ... | ... ... ... 0 1 rd (rd)m | 0 0 1 0 1 ... m s2 ... sd + 11 ... ... ... 1d 21 ... ... ... 2d ... ... ... ... r1 ... ... ... rd Odhad vektoru neznámých parametrů metodou nejmenších čtverců je dán vztahem ^ = ^0 ^1 ... ^m s2 ... sd = (X X)-1 X Y Srovnání obou modelů na konkrétních příkladech Závěrem proveďme jednoduché srovnání obou předchozích modelů. Jak je vidět z následujících obrázků i velice jednoduchý model s konstantním trendem po celou dobu sezóny dává u obou časových řad lepší výsledky než komplikovanější a výpočetně náročnější druhý model, viz např. hodnota SSE = n i=1 (Yi - Yi)2 . 3.3. Klasická dekompozice časových řad 74 5 10 15 20 25 30 35 100 200 300 400 500 600 700 800 900 1000 1100 Zaznamy okr.spravy o mesicnich poctech deti ve skolach v prirode(1980-1982). Metoda konstantniho trendu v jednotlivych sezonach-aditivni model Y jk =m j +s k + jk j=1,...,3 k=1,...,12 s 1 =-74.9444 s 2 =6.7222 s 3 =-16.9444 s 4 =-55.6111 s 5 =-180.2778 s 6 =-182.9444 s 7 =-45.9444 s 8 =97.7222 s 9 =183.3889 s 10 =137.7222 s 11 =-43.9444 s 12 =175.0556 SSE=52446.611111 m 1 =330.5833 m 2 =567.6667 m 3 =801.5833 1 2 3 4 5 6 7 8 9 10 11 12 1 -200 -150 -100 -50 0 50 100 150 200 sezonnislozkysk (k=1,...,12) 5 10 15 20 25 30 35 100 200 300 400 500 600 700 800 900 1000 1100 Metoda polynomickeho trendu za cele obdobi+sezonnost-aditivni model Zaznamy okr.spravy o mesicnich poctech deti ve skolach v prirode(1980-1982). Y jk = 0 + 1 t+s 2 +...+s d + jk j=1,...,3 k=1,...,12 t=(j-1)d+k s 2 = 62.04167 s 3 = 18.75 s 4 =-39.54167 s 5 =-183.8333 s 6 = -206.125 s 7 = -88.75 s 8 = 35.29167 s 9 = 101.3333 s 10 = 36.04167 s 11 = -165.25 s 12 = 34.125 SSE=52466.666667 0 =236.5417 1 =19.6250 1 2 3 4 5 6 7 8 9 10 11 12 1 -250 -200 -150 -100 -50 0 50 100 150 sezonnislozkysk (k=1,...,12) 3.3. Klasická dekompozice časových řad 75 10 20 30 40 50 60 70 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 Monthly accidental deaths in the USA, 1973 - 1978 Metoda konstantniho trendu v jednotlivych sezonach-aditivni model Y jk =m j +s k + jk j=1,...,6 k=1,...,12 s 1 =-743.7361 s2 =-1503.9028 s3 =-723.9028 s4 =-522.9028 s 5 =338.4306 s6 =807.5972 s7 =1665.0972 s 8 =961.4306 s 9 =-87.4028 s10 =196.9306 s11 =-320.5694 s 12 =-67.0694 SSE=3955160.263889 m 1 =9651.7500 m2 =8718.5000 m3 =8585.8333 m4 =8396.7500 m 5 =8576.8333 m6 =8796.7500 1 2 3 4 5 6 7 8 9 10 11 12 1 -2000 -1500 -1000 -500 0 500 1000 1500 2000 sezonnislozkys k (k=1,...,12) 10 20 30 40 50 60 70 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 Metoda polynomickeho trendu za cele obdobi+sezonnost-aditivni model Monthly accidental deaths in the USA, 1973 - 1978 Y jk = 0 + 1 t+ 2 t2 +s 2 +...+s d + jk j=1,...,6 k=1,...,12 t=(j-1)d+k s 2 =-740.26057 s3 = 57.992433 s4 = 275.59236 s5 = 1151.8725 s 6 = 1634.333 s 7 = 2503.4736 s 8 = 1809.7946 s 9 = 769.29573 s 10 = 1060.3105 s 11 = 547.83883 s 12 = 804.71409 SSE=4400371.750159 0 =9133.8707 1 =-71.9782 2 =0.8265 1 2 3 4 5 6 7 8 9 10 11 12 1 -1000 -500 0 500 1000 1500 2000 2500 3000 sezonnislozkys k (k=1,...,12) 3.3. Klasická dekompozice časových řad 76 10 20 30 40 50 60 70 80 60 80 100 120 140 160 180 200 Monthly housing starts, construction contracts and average new home mortgate rates. Jan 83 - Oct 89. Metoda konstantniho trendu v jednotlivych sezonach-aditivni model Yjk =mj +sk +jk j=1,...,7 k=1,...,12 s 1 =-14.5972 s2 =-33.4972 s 3 =-37.4790 s 4 =-37.0076 s 5 =-1.5790 s 6 =21.0638 s 7 =25.0638 s 8 =27.6210 s 9 =16.1924 s 10 =12.0495 s 11 =4.3352 s 12 =10.9638 SSE=6168.987177 m 1 =145.8700 m2 =147.3833 m 3 =143.6250 m 4 =151.7833 m 5 =137.0833 m 6 =123.7500 m 7 =117.6583 1 2 3 4 5 6 7 8 9 10 11 12 1 -40 -30 -20 -10 0 10 20 30 sezonnislozkys k (k=1,...,12) 10 20 30 40 50 60 70 80 60 80 100 120 140 160 180 200 Metoda polynomickeho trendu za cele obdobi+sezonnost-aditivni model Monthly housing starts, construction contracts and average new home mortgate rates. Jan 83 - Oct 89. Y jk = 0 + 1 t+ 2 t2 +s 2 +...+s d + jk j=1,...,7 k=1,...,12 t=(j-1)d+k s 2 =0.746834 s 3 = 36.4787 s4 = 59.4528 s 5 = 63.812 s 6 = 66.7562 s 7 = 55.7426 s 8 = 52.0426 s 9 = 44.7991 s10 = 51.9264 s 11 = 20.6894 s 12 = 2.17645 SSE=6407.767450 0 =99.2911 1 =0.7714 2 =-0.0140 1 2 3 4 5 6 7 8 9 10 11 12 1 0 10 20 30 40 50 60 70 sezonnislozkysk (k=1,...,12) 3.3. Klasická dekompozice časových řad 77 10 20 30 40 50 60 70 80 10 12 14 16 18 20 22 24 26 28 Monthly housing starts, construction contracts and average new home mortgate rates. Jan 83 - Oct 89. Metoda konstantniho trendu v jednotlivych sezonach-aditivni model Yjk =mj +sk +jk j=1,...,7 k=1,...,12 s1 =-2.4461 s2 =-3.4113 s3 =-4.8158 s 4 =-4.7619 s 5 =0.3358 s6 =1.1881 s7 =2.7885 s8 =3.4471 s 9 =1.9625 s 10 =2.4702 s11 =1.1171 s12 =1.2892 SSE=73.522167 m1 =16.4329 m2 =17.4090 m3 =19.0226 m 4 =20.4023 m 5 =21.3367 m6 =21.5169 m7 =21.1252 1 2 3 4 5 6 7 8 9 10 11 12 1 -5 -4 -3 -2 -1 0 1 2 3 4 sezonnislozkysk (k=1,...,12) 10 20 30 40 50 60 70 80 10 12 14 16 18 20 22 24 26 28 Metoda polynomickeho trendu za cele obdobi+sezonnost-aditivni model Monthly housing starts, construction contracts and average new home mortgate rates. Jan 83 - Oct 89. Y jk = 0 + 1 t+ 2 t2 +s 2 +...+s d + jk j=1,...,7 k=1,...,12 t=(j-1)d+k s 2 =-0.033788 s 3 = 4.979 s4 = 5.749 s 5 = 7.2698 s 6 = 7.8515 s 7 = 6.2928 s 8 = 6.729 s 9 = 5.307 s10 = 5.413 s 11 = 2.3544 s 12 = 1.3124 SSE=80.231739 0 =10.4319 1 =0.1885 2 =-0.0013 1 2 3 4 5 6 7 8 9 10 11 12 1 -1 0 1 2 3 4 5 6 7 8 sezonnislozkysk (k=1,...,12) 3.3. Klasická dekompozice časových řad 78 3.3.5 Analýza reziduální (náhodné) složky Někdy se stane, že časová řada předložená k analýze nevykazuje při zběžné prohlídce nebo grafickém znázornění výskyt žádné systematické složky, takže se zdá, že je tvořena pouze bílým šumem. Pro jistotu však je vhodné provést nějaký objektivní statistický test, který by tuto hypotézu potvrdil. Testy náhodnosti se také někdy používají v situaci, kdy je nutné po provedení dekompozice určité řady ověřit, zda jsme skutečně všechny systematické složky z řady eliminovali tak, že zbytek po eliminaci jsou opravdu už jen náhodné fluktuace ve tvaru bílého šumu. Takové použití testů z této kapitoly však není z čistě teoretického hlediska korektní, neboť testované časové řady jsou pak již zatíženy určitými odhadovými chybami a často ani nejsou tvořeny nekorelovanými veličinami. 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 vzestupna tendence 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 sestupna tendence 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 odpuzovani sousedu 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 zmena podminek 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 odlehle pozorovani 0 1 2 3 4 5 6 7 8 9 10 11 0 2 4 6 8 10 periodicita Obrázek 3.7: Typy porušení náhodnosti. TESTY NÁHODNOSTI Testy náhodnosti jsou testy statistické hypotézy, že Y1, . . . , Yn jsou náhodným výběrem z nějakého rozdělení s distribuční funkcí F. Jinými slovy test náhodnosti je test nulové hypotézy, která tvrdí, že sdružená distribuční funkce je součinem marginálních, tj. platí P(Y1 y1, . . . , Yn yn) = n i=1 F(yi). Ve všech testech budeme testovat tuto nulovou hypotézu: H0 : časová řada je realizací vzájemně nezávislých stejně rozdělených náhodných veličin Jestliže časová řada má nenulovou střední hodnotu , budeme vyšetřovat centrovanou časo- vou řadu Zt = Yt - . 3.3. Klasická dekompozice časových řad 79 TEST ZALOŽENÝ NA ZNAMÉNKÁCH PRVNÍCH DIFERENCÍ (The Difference-Sign Test) Označme S počet kladných prvních diferencí posloupnosti Y1, . . . , Yn. Jsou-li některé sousední hodnoty stejné, vyškrtněme je kromě té první. Věta 3.3.7. Za platnosti H0 má statistika S tyto vlastnosti: ES = n-1 2 DS = n+1 12 U = S - ES DS A N(0, 1) Důkaz. Tento test poprvé zavedli Moore, G.H., Wallis, W.A.: Time Series Significance Tests Based on Signs of Differences, Journal of the American Statistical Association, Vol. 38, Issue 222, (Jun., 1943), 153-164. Položme Vt = 1 Yt+1 > Yt 0 Yt+1 < Yt Vt A 1 2 EVt = 1 2 DVt = (1 - 1 2 )1 2 = 1 4 . Pomocí Vt vyjádřeme S = n-1 t=1 Vt ES = (n - 1)1 2 a DS = n-1 t=1 DVt +2 n-1 t,s=1,t Yt+1 a Yt je dolní bod zvratu, pokud Yt-1 > Yt < Yt+1 jsou-li některé sousední hodnoty stejné vyškrtnout až na první. Věta 3.3.8. Za platnosti H0 má statistika SZ tyto vlastnosti: ESZ = 2(n - 2) 3 DSZ = 16n - 29 90 UZ = SZ - ESZ DSZ A N(0, 1) Důkaz. viz Moore, G.H., Wallis, W.A.: A Significance Test for Time Series Analysis, Journal of the American Statistical Association, Vol. 36, Issue 215, (Sep., 1941), 401-409 2 Hypotézu H0 zamítáme, pokud počet bodů zvratu je příliš velký nebo malý, tj. pokud |UZ| u1- 2 . TEST ZALOŽENÝ NA KENDALOVĚ KOEFICIENTU Označme S . . . počet takových dvojic Yt a Ys, že Yt < Ys pro t < s. jsou-li některé hodnoty stejné vyškrtnou se až na první. Protože dvojic (Yt, Ys), kde t = s, je n 2 = n(n-1) 2 a možnosti Yt < Ys a Yt > Ys mají za platnosti nulové hypotézy stejnou pravděpodobnost, pak ES = n(n-1) 2 1 2 = n(n-1) 4 , proto častěji se používá tzv. Kendalův koeficient pořadové korelace definovaný = 4S n(n - 1) - 1 . Věta 3.3.9. Za platnosti H0 má statistika tyto vlastnosti: E = 0 D = 2(2n + 5) 9n(n - 1) U = - E D A N(0, 1) Důkaz. viz Mann, H.B.: Non-parametric tests against trend, Econometrica, 13, (1945), 245- 259. 2 Hypotézu H0 zamítáme, pokud počet definovaných dvojic je příliš velký nebo malý, tj. pokud |U | u1- 2 . 3.3. Klasická dekompozice časových řad 81 TEST ZALOŽENÝ NA SPEARMANOVĚ KOEFICIENTU Nechť náhodné veličiny R1, . . . , Rn označují pořadí hodnot náhodných veličin ve výběru Y1, . . . , Yn. Spearmanův koeficient pořadové korelace je definovaný takto: = 1 - 6 n(n2 - 1) n i=1 (i - Ri)2 . Věta 3.3.10. Za platnosti H0 má statistika tyto vlastnosti: E = 0 D = 1 (n - 1) U = - E D A N(0, 1) Důkaz. viz Daniels,H.E.: Rank correlation and population models, J. R. Statist. Soc., B, 12 (1950), 171-181. 2 Hypotézu H0 zamítáme, pokud |U| u1- 2 TEST ZALOŽENÝ NA POČTU ITERACÍ DVOU PRVKŮ Mějme posloupnost n prvků, kde je n1 prvků prvního druhu a n2 prvků druhého druhu: n = n1 + n2. Vybírejme náhodně z tohoto souboru postupně prvky (bez vracení) až do úplného vybrání. Všech takto vybíraných n-tic je n n1 = n n2 a všechny mají stejnou pravděpodobnost 1 ( n n1 ) = 1 ( n n2 ) . Iterací (sérií) rozumíme každou částečnou posloupnost této posloupnosti obsa- hující pouze prvky jednoho druhu a ohraničenou z každé strany prvkem druhého druhu, nebo začátkem, nebo koncem posloupnosti. Příklad iterace (série): aa 1 bb 2 a 3 bbb 4 aa 5 bb 6 aaaa 7 Nechť SI je počet iterací. Pak platí ESI = 2n1n2 n + 1 a DSI = 2n1n2(2n1n2 - n) n2(n - 1) . Pokud n1 , n2 a n1 n2 c (0, 1), pak UI = SI - ESI DSI A N(0, 1). Iterace se používají právě k ověřování hypotézy, že daná posloupnost prvků dvojího druhu je uspořádána náhodně, tj. že při jejím vzniku mělo každé z možných uspořádání stejnou prav- děpodobnost. Jako alternativní hypotézy přicházejí v úvahu 3 možnosti: 3.3. Klasická dekompozice časových řad 82 H1 prvky stejného druhu mají tendenci sdružovat se ve skupiny, takže posloupnosti s malým počtem iterací mají větší pravděpodobnost, než ostatní. Hypotézu H0 zamítáme proti alternativě H1, jestliže UI u . H2 prvky stejného druhu tvoří zřídkakdy větší skupiny, takže posloupnosti s větším počtem iterací mají větší pravděpodobnost, než ostatní. Hypotézu H0 zamítáme proti alter- nativě H2, jestliže UI u1- . H3 spojuje obě předchozí. Hypotézu H0 zamítáme proti alternativě H3, jestliže |UI| u1- 2 . MEDIÁNOVÝ TEST Mediánový test pro testování hypotézy, že Y1, . . . , Yn je posloupnost nezávislých pozorování, volí výběrový medián ~Y = Y" n+1 2 " n liché 1 2 Y(n 2 ) + Y(n 2 +1) n sudé , kde Y(1), . . . , Y(n) je uspořádáním Y1, . . . , Yn. Je-li Yi > ~Y nahradíme jej prvkem a Yi < ~Y nahradíme jej prvkem b , pro i = 1, . . . , n. Medián rozdělí n prvků přesně na dvě skupiny n = 2m. Označme Sm počet iterací (sérií) v posloupnosti vytvořené podle předchozího pravidla. Věta 3.3.11. Za platnosti H0 má statistika Sm tyto vlastnosti: ESm = m + 1 DSm = m(m-1) (2m-1) USm = Sm - ESm DSm A N(0, 1). Alternativa H1 je citlivá vůči vzestupné či sestupné tendenci, dlouhým periodám, náhlé změně podmínek. Alternativa H2 je citlivá vůči krátké periodě, odpuzování sousedů. 3.3. Klasická dekompozice časových řad 83 0 20 40 60 80 100 120 140 160 180 200 -1 0 1 2 3 4 5 Y t = + t ; =2; t WN(0,1) H0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 1 0.122169444 100 Test based on runs up and down: 1 1.0108213 138 Test based on Kendall coefficient: 1 1.20890794 10522 Test based on Spearman coefficient: 1 1.204696 0.085398635 Median test: 1 0.141778031 102 0 20 40 60 80 100 120 140 160 180 200 -10 -5 0 5 10 Y t = + a sin(4 t/n) + t ; n=200; =2; a=5; t WN(0,1) H0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 1 1.3438639 105 Test based on runs up and down: 1 0.84235108 127 Test based on Kendall coefficient: 0 5.499263 7348 Test based on Spearman coefficient: 0 5.6226341 -0.39857796 Median test: 0 13.185357 8 Obrázek 3.8: Demonstrace citlivosti testů náhodnosti na simulovaných datech. 3.3. Klasická dekompozice časových řad 84 0 10 20 30 40 50 60 600 800 1000 1200 1400 1600 1800 Ex 3.5: plastics.dat Monthly sales of product A for a plastics manufacturer H 0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 0 5.5441595 42 Test based on runs up and down: 0 4.5601364 24 Test based on Kendall coefficient: 0 4.2221884 1216 Test based on Spearman coefficient: 0 3.8574334 0.50219505 Median test: 0 4.6874741 13 0 20 40 60 80 100 120 20 40 60 80 100 120 140 160 180 Ex 7.5: condmilk.dat Stocks of evaporated and sweetened condensed milk. H 0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 0 2.3618875 52 Test based on runs up and down: 0 11.92609 24 Test based on Kendall coefficient: 1 0.31755367 3500 Test based on Spearman coefficient: 1 0.24887228 -0.022814084 Median test: 0 6.9671546 23 0 20 40 60 80 100 120 20 30 40 50 60 70 80 90 100 Ex 8.8: schizo.dat Daily perceptual speed scores for a schizophrenic patient. H0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 1 0.95668921 55 Test based on runs up and down: 1 0.29464381 78 Test based on Kendall coefficient: 0 7.1177388 2001 Test based on Spearman coefficient: 0 7.0656999 -0.64771165 Median test: 0 7.211366 21 0 20 40 60 80 100 120 200 300 400 500 600 700 800 900 1000 1100 writing.dat Industry sales for printing and writing paper Jan 1963 - Dec 1972(in Thousands of french francs). H 0 : Test hypothesis of data being white noise (1=not rejected; 0=rejected): Test based on sign differences: 1 1.7320508 65 Test based on runs up and down: 0 2.7633623 66 Test based on Kendall coefficient: 0 9.0275973 5560 Test based on Spearman coefficient: 0 7.814135 0.71632058 Median test: 0 6.7838084 24 Obrázek 3.9: Demonstrace citlivosti testů náhodnosti na reálných datech. Kapitola 4 Spektrální analýza jednorozměrných časových řad 4.1 Úvod Pojem spektra se vyskytuje nejen v teorii náhodných procesů, ale také v matematice, fyzice a technice. Jestliže nějaký proces vlnění je součtem harmonických vlnění (tzv. harmonik), tak spektrum procesu vlnění se nazývá funkce, která popisuje rozdělení amplitud podle jednotlivých frekvencí. Spektrum ukazuje, která vlnění převládají v daném procesu a jaká je jeho vnitřní struktura. Spektrum v případě stacionárního náhodného procesu dává rozdělení rozptylů ná- hodných amplitud podle různých frekvencí vlnění. V celém tomto odstavci proto budeme předpokládat, že náhodný proces {Yt, t T} je stacionární, centrovaný a druhého řádu (tj. s konečnými druhými momenty). Významnou vlastností stacionárních náhodných procesů je vlastnost, že jeho autokovariační funkci lze vyjádřit jako (nespočetný) součet harmonických funkcí s různými frekvencemi a amplitudami. ˇ Podle Herglotzovy věty platí: Je-li {Yt, t Z} stacionární posloupnost, pak se její autokovarianční funkce (t) dá vyjá- dřit ve tvaru (t) = - eit dF(), kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. ˇ Podle Bochnerovy věty platí: Je-li {Yt, t R} stacionární proces spojitý podle středu, pak se jeho autokovarianční funkce (t) dá vyjádřit ve tvaru (t) = - eit dF(), kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. 85 4.2. PERIODOGRAM 86 Vzorci (t) = - eit dF() resp. (t) = - eit dF() se říká spektrální rozklad kovarianční funkce. Funkce F() se nazývá spektrální distri- buční funkce. ˇ Je-li F() absolutně spojitá, pak existuje taková funkce f(), že pro náhodné stacionární posloupnosti, resp. pro stacionární náhodné procesy platí F() = - f(x)dx resp. F() = - f(x)dx. (4.1) ˇ K existenci spektrální hustoty stacionární náhodné posloupnosti, resp. spojitého stacio- nární náhodného procesu, stačí, aby pro její kovarianční funkci platilo t=- |(t)| < resp. - |(t)|dt < . ˇ Existuje-li spektrální hustota f() stacionární posloupnosti a má-li variaci konečnou na -, , pak platí f() = 1 2 t=- e-it (t) (4.2) ve všech bodech spojitosti funkce f(), což je skoro všude vzhledem k Lebesqueově míře. ˇ Existuje-li spektrální hustota f() spojitého stacionárního procesu a je-li autokovarianční funkce absolutně integrovatelná, tj. - |(t)|dt < , pak f() = 1 2 - e-it (t) dt. (4.3) ˇ Spektrální hustota f() reálného spojitého stacionárního procesu nebo reálné stacionární posloupnosti je sudá funkce v tom smyslu, že pro ni platí f() = f(-) (4.4) skoro všude vzhledem k Lebesqueově míře. 4.2 PERIODOGRAM V dalším budeme předpokládat, že {Yt, t Z} je centrovaná stacionární náhodná po- sloupnost. 4.2. PERIODOGRAM 87 Definice 4.2.1. Nechť Y1, . . . , Yn jsou pozorování náhodné posloupnosti {Yt, t Z}. Pak periodogram definujeme vztahem In() = 1 2n n t=1 Yte-it 2 [-, ]. Lemma 4.2.2. Položme An() = 2 n n t=1 Yt cos t Bn() = 2 n n t=1 Yt sin t, pak platí In() = 1 4 A2 n() + B2 n() Důkaz. In() = 1 2n n t=1 Yte-it 2 = 1 2n n t=1 Yt cos t - i n t=1 Yt sin t 2 = = 1 2n n t=1 Yt cos t 2 + n t=1 Yt sin t 2 = 1 4 A2 n() + B2 n() . 2 Poznámka 4.2.3. Někteří autoři definují periodogram poněkud jinak: I n() = 2 n n t=1 Yte-it 2 = A2 n() + B2 n() = 4In(). Lemma 4.2.4. Pokud označíme pro k = 0, 1, . . . , n - 1 Ck = 1 n - k n-k t=1 YtYt+k C k = 1 n n-k t=1 YtYt+k pak platí In() = 1 2 C0 + 2 n-1 k=1 1 - k n Ck cos k = 1 2 C 0 + 2 n-1 k=1 C k cos k 4.2. PERIODOGRAM 88 Důkaz. In() = 1 2n n t=1 Yt cos t 2 + n t=1 Yt sin t 2 = 1 2n n t=1 Yt cos t n s=1 Ys cos s + n t=1 Yt sin t n s=1 Ys sin s = 1 2n n t=1 n s=1 YtYs (cos t cos s + sin t sin s) = 1 2n n t=1 n s=1 YtYs cos (s - t) Zavedeme-li dále substituci k = s - t , pak -n + 1 k n - 1 a 1 t n 1 s = t + k n 1 - k t n - k týká se kladných k max(1, 1 - k) t min(n, n - k). týká se záporných k a pak platí In() = 1 2n n-1 k=-n+1 cos k min(n,n-k) t=max(1,1-k) YtYt+k. Nyní vezměme zvlášť případy, kdy k = 0 a ostatní, přičemž využijme faktu, že funkce cos je sudou funkcí. Dostaneme proto In() = 1 2 1 n n t=1 Y 2 t C0 + 1 2 -1 k=-n+1 n - |k| n cos k 1 n - |k| n t=1-k YtYt+k C-k=Ck + 1 2 n-1 k=1 n - k n cos k 1 n - k n-k t=1 YtYt+k Ck = = 1 2 n-1 k=-(n-1) 1 - |k| n Ck cos k = 1 2 C0 + 2 n-1 k=1 1 - k n Ck cos k In() = 1 2 1 n n t=1 Y 2 t C 0 + 1 2 -1 k=-n+1 cos k 1 n n t=1-k YtYt+k C -k=C k + 1 2 n-1 k=1 cos k 1 n n-k t=1 YtYt+k C k = 1 2 C 0 + 2 n-1 k=1 C k cos k 2 4.2. PERIODOGRAM 89 Poznámka 4.2.5. K numerickému výpočtu hodnot periodogramu se často používají právě předchozí vzorce. Poznámka 4.2.6. Pro teoretické účely bývá výhodnější tato varianta In() = 1 2 n-1 k=-(n-1) 1 - |k| n Ck cos k = 1 2 n-1 k=-(n-1) C k cos k. Pro náhodnou posloupnost {Yt, t T Z} platí f() = 1 2 t=- (t) cos t. Veličiny 1 - k n Ck, (resp. C k) můžeme považovat za jakýsi odhad (k) a periodogram se tudíž dá považovat za empirický odhad spektrální hustoty. Vlastnosti tohoto odhadu udává následující věta. Věta 4.2.7. Jestliže {Yt, t T Z} je stacionární náhodná posloupnost s nulovou střední hodnotou a se spojitou spektrální hustotou f(), pak má periodogram In() následující vlastnosti: lim n EIn() = f() -, . lim n DIn() = f2 () = 0, (-, ), 2f2 () = 0, . Důkaz. Nejprve připomenme následující vztahy: (1) n t=1 eit = ei ein - 1 ei - 1 = ei ei n 2 1 2i ei n 2 - e-i n 2 ei 1 2 1 2i ei 1 2 - e-i 1 2 = ei n+1 2 sin n 2 sin 1 2 pro = 2k n t=1 eit = n pro = 2k (2) n t=1 eit 2 = n t=1 eit n s=1 e-is = n t=1 n s=1 ei(t-s) = ei n+1 2 sin n 2 sin 1 2 2 = sin2 n 2 sin2 1 2 . . . tzv. Fejérovo jádro pro = 2k n t=1 eit 2 = n2 pro = 2k 4.2. PERIODOGRAM 90 (3) Nechť funkce g() je integrovatelná na -, , pak pro každý bod -, , v němž je g() spojitá, platí: lim n 1 2n - sin2 n 2 ( - ) sin2 1 2 ( - ) g() d = g() Jde o známou vlastnost Fejérova jádra (viz např. Anděl, J.: Statistická analýza časových řad, Praha 1976, str. 97). Nyní počítejme střední hodnotu EIn() = E 1 2n n t=1 Yte-it 2 = E 1 2n n t=1 n s=1 YtYse-i(t-s) = 1 2n n t=1 n s=1 e-i(t-s) E(YtYs) (t,s)=(t-s) = 1 2n n t=1 n s=1 e-i(t-s) - ei(t-s) f() d =(t-s) = 1 2n - n t=1 n s=1 e-i(t-s) ei(t-s) f() d = 1 2n - n t=1 n s=1 eit(-) e-is(-) =| Pn t=1 eit(-) | 2 f() d = 1 2n - sin2 n 2 ( - ) sin2 1 2 ( - ) f() d viz (3) --- n f() Tím jsme dokázali první vlastnost. Důkaz druhé vlastnosti provedeme pouze pro Yt IID(0, 2 ), přičemž budeme předpokládat, že pro t Z navíc platí EY 4 t = 4 < . Více lze najít např. v knize Anděl. J.: Statistická analýza časových řad, Praha 1976, str. 100-110. Pro bílý šum platí E(YtYs) = (t - s) = 2 s = t, 0 jinak tj. (h) = 2 h = 0, 0 jinak. Proto EIn() = 1 2n n t=1 n s=1 e-i(t-s) E(YtYs) = 1 2n n t=1 EY 2 t = 1 2n n2 = 2 2 . Dále pomocí autokovarianční funkce počítejme spektrální hustotu f() = 1 2 t=- e-it (t) = 2 2 -, , 0 jinak 2 2 /2 - Než budeme pokračovat dále, poznamenejme, že pro Yt IID(0, 2 ) platí E (YtYsYuYv) = 4 t = s = u = v, 4 t = s = u = v, 4 t = u = s = v, 4 t = v = s = u, 0 jinak. 4.2. PERIODOGRAM 91 Pak pro 1, 2 -, počítejme E [In(1)In(2)] = E 1 2n n t=1 Yte-it1 2 1 2n n t=1 Yte-it2 2 = 1 42n2 E n t=1 n s=1 e-it1 eis1 YtYs n u=1 n v=1 e-iu2 eiv2 YuYv = 1 42n2 n t=1 n s=1 n u=1 n v=1 ei(s-t)1 ei(v-u)2 E (YtYsYuYv) = 1 42n2 n t=1 EY 4 t + n t=1 n u=1,t=u EY 2 t Y 2 u + n t=1 n s=1,s=t ei(s-t)1 ei(s-t)2 EY 2 t Y 2 s + + n t=1 n u=1,u=t ei(u-t)1 ei(t-u)2 EY 2 t Y 2 u = 1 42n2 n4 + n(n-1)4 + 4 n t=1 n s=1 ei(s-t)(1+2) - n4 + 4 n t=1 n s=1 ei(s-t)(1-2) - n4 = 1 42 4 + 4 - 34 n + 4 n2 n t=1 eit(1+2) 2 + n t=1 eit(1-2) 2 Pro 1 = 2 a 1 2 = 2k (k = 0, 1) platí E [In(1)In(2)] = 1 42 4 + 4 - 34 n + 4 n2 sin2 n 2 (1 + 2) sin2 1 2 (1 + 2) + sin2 n 2 (1 - 2) sin2 1 2 (1 - 2) Pro 1 = 2 = = 2k (k = 0, 1) platí EI2 n() = 1 42 4 + 4 - 34 n + 4 n2 sin2 n sin2 + n2 = 1 42 24 + 4 - 34 n + 4 n2 sin2 n sin2 Pro 1 = 2 = = 2k (k = 0, 1) platí EI2 n() = 1 42 4 + 4 - 34 n + 4 n2 n2 + n2 = 1 42 34 + 4 - 34 n Odtud díky vztahu DIn() = EI2 n() - (EIn())2 , kde EIn() = 2 2 dostaneme DIn() = 1 42 4 + 4-34 n + 4 n2 sin2 n sin2 --- n 4 42 = f2 () pro = 0, (-, ), 1 42 24 + 4-34 n --- n 2 4 42 = 2f2 () pro = 0, . Poznámka 4.2.8. Z předchozí věty vyplývá 1. Periodogram In() je asymptoticky nestranným odhadem spektrální hustoty. 2. Periodogram In() není konzistentním odhadem spektrální hustoty, neboť jeho rozptyl nekonverguje k nule, vzrůstá-li neomezeně délka posloupnosti n. 4.2. PERIODOGRAM 92 0 50 100 150 200 -5 -4 -3 -2 -1 0 1 2 3 4 5 y t = sin(250t /n) + sin(2125t /n) + t ; t ~ WN(0,1); n=1001 250 300 350 400 450 -4 -3 -2 -1 0 1 2 3 4 5 500 550 600 650 700 -6 -4 -2 0 2 4 6 750 800 850 900 950 -6 -4 -2 0 2 4 6 Obrázek 4.1: Simulovaná časová řada. 4.2. PERIODOGRAM 93 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=125 n=251 log(I n ( j )) Periodogram I n ( j ) j =j (2 /m) j=1,...m 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=250 n=501 log(I n ( j )) 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=375 n=751 log(I n ( j )) 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=500 n=1001 log(I n ( j )) Obrázek 4.2: Ukázka nekonzistentnosti periodogramu, neboť jeho variabilita se ne- snižuje se vzrůstajícím počtem pozorování. 4.2. PERIODOGRAM 94 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=125 n=250 log(I n ( j )) j =j(2/m)j=1,...m Periodogram z prvni casti dat 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=125 n=250 log(I n ( j )) j =j(2/m)j=1,...m Redukce rozptylu zprumerovanim 4 periodogramu bez prekryvani dat 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=125 n=250 log(I n ( j )) j =j(2/m)j=1,...m Redukce rozptylu zprumerovanim 7 periodogramu s prekryvani dat delky 125 0 0.5 1 1.5 2 2.5 3 -10 -5 0 5 m=125 n=250 log(I n ( j )) j =j(2/m)j=1,...m Redukce rozptylu zprumerovanim 12 periodogramu s prekryvani dat delky 187 Obrázek 4.3: Na obrázku je demonstrován tzv. Welchův přístup, který vycházel ze zmenšo- vání variability průměrováním periodogramu přes nepřekrývající, resp. překrývající podúseky dat. Periodogram lze ještě více vyhladit použitím různých vah (viz. Neparametrické odhady spektrální hustoty) 4.2. PERIODOGRAM 95 Příklad 4.2.9. Metoda skrytých period. Uvažujme nyní takové časové řady, které můžeme rozložit na součet harmonických frek- vencí, jejichž délky period lze vyjádřit jako podíl Tk = n k , kde n je počet naměřených hodnot a 0 < k n. Označme dále fk = 1 Tk = k n k-tá frekvence k = 2fk = 2 k n k-tá úhlová frekvence Maximální délka periody, kterou jsme schopni určit, je rovna počtu pozorování, tj. Tmax = n , tedy k = 1 a minimální frekvence má velikost min = 2 n . Nejkratší zjistitelná perioda je Tmin = 2 . Této délce odpovídá frekvence max = , tzv. Nyquistova frekvence. Z předchozích úvah vyplývá, že k může nabývat hodnot: k = 1, 2, . . . , n 2 pro n sudé 1, 2, . . . , n-1 2 pro n liché. Model časové řady pak můžeme zapsat ve tvaru Yt = st + t , kde t označuje ekvidistatní časové okamžiky měření (pro jednoduchost předpokládáme, že intervaly mají jednotkovou velikost), n je počet naměřených hodnot, t je bílý šum: t WN(0, 2 ), tj. Et = 0; Dt = 2 ; C(t, s) = 0; t = s, st je periodická funkce tvaru st = p j=1(j cos tj + j sin tj) pro n liché p j=1(j cos tj + j sin tj) + n 2 (-1)t pro n sudé kde j, j R, p N jsou daná čísla a nazýváme je parametry modelu. Pokud n je sudé číslo, může být mezi vybranými periodami i perioda délky Tmin = 2 , což odpovídá frekvenci max = . Vzhledem k tomu, že platí sin n = 0 a cos n = (-1)n , proto zapisujeme koeficient odpovídající této frekvenci zvlášť. Ekvivalentně můžeme st napsat ve tvaru: st = p j=1 j cos(tj - j) pro n liché p j=1 j cos(tj - j) + n 2 (-1)t pro n sudé kde k = 2 k + 2 k je amplituda k-té harmonické složky k = arctan k k k > 0 arctan k k + k < 0 je fázový posun k-té harmon. složky Pokud zahrneme do počtu frekvencí i frekvenci nulovou (kdy k = 0), pak přibude člen, který označíme 0 2 , přitom 0 = 0. 4.2. PERIODOGRAM 96 V dalším pro jednoduchost předpokládejme, že n je liché číslo. Použijeme-li trigonometrické vyjádření komplexního čísla: a + bi = r(cos + i sin ) a vztahů cos kx = 1 2 eikx + e-ikx sin kx = 1 2i eikx - e-ikx pak dostaneme Fourierovu řadu konečné délky st = 0 2 + p k=1 (k cos tk + k sin tk) = 0 2 + p k=1 1 2 eitk + e-itk k + 1 2i eitk - e-itk k = 0 2 c0 + p k=1 1 2 (k - ik) ck eitk + 1 2 (k + ik) ck e-itk = p k=-p ckeitk , přičemž pro k = 1, . . . , p platí k = ck + ck = 2Re(ck) k = i(ck - ck) = -2Im(ck). Vypočítejme střední hodnotu tohoto modelu EYt = st, autokovarianční funkce je tvaru (k) = C(Yt, Yt+k) = 0 k > 0, t R 2 k = 0, t R a spektrální hustota f() = 1 2 t=- e-it (t) = 2 2 -, , 0 jinak 2 2 /2 - Všimněme si periodogramu této náhodné posloupnosti In() = 1 2n n t=1 p j=-p cjeitj + t e-it 2 = 1 2 1 n n t=1 p j=-p cjeit(j -) I (1) n () + 1 n n t=1 te-it I (2) n () 2 4.2. PERIODOGRAM 97 Upravujme dále I(1) n () = 1 n p j=-p cj n t=1 eit(j -) součet geom.řady = 1 n p j=-p cjei(j -) ein(j -) - 1 ei(j -) - 1 gn(j -) Je-li různé od všech 1, . . . , p, pak platí zřejmě lim n I(1) n () = 0. Existuje-li takové k, že platí = k, pak I(1) n (k) = n ck n + 1 n p j=-p,j=k gj(j - ) = nck + 1 n p j=-p,j=k gn(j - ) 0 a pro n vzrůstá jeho absolutní hodnota nade všechny meze. Druhý člen I (2) n () je náhodná veličina s nulovou střední hodnotou a s rozptylem DI(2) n () = E I(2) n () 2 = E 1 n n t=1 te-it 2 = 1 n E n t=1 te-it n s=1 teis = 1 n n t=1 n s=1 e-i(t-s) Ets nekorel. = 1 n n t=1 E2 t = 1 n n2 = 2 Pokud bude platit, že náhodná posloupnost splňuje model Yt = st + t = 0 2 + p j=1 (j cos tj + j sin tj) + t t WN(0, 2 ), bude mít periodogram pro velká n v bodech = 1, . . . , p výrazně velké hodnoty (řádu n), jinde jeho hodnoty budou relativně malé, budou kolísat kolem hodnoty 2 2 . Periodogram je tedy dobrým ukazatelem periodicit. Z uvedeného příkladu vyplývá, že významná lokální maxima v průběhu periodogramu by měla identifikovat periodickou strukturu uvažovaného modelu tak, že vyznačí neznámé frek- vence 1, . . . , p. Nějaký vhodný statistický test by pak měl rozhodnout, jaké hodnoty periodogramu můžeme opravdu považovat za významně velké ve srovnání s hodnotami ostatními. 4.2. PERIODOGRAM 98 4.2.1 Test R. A. Fishera R. A. Fisher odvodil test, kterým se dá zjistit významnost nejvyšších hodnot periodogramu. Uvažujme posloupnost nezávislých náhodných veličin Y1, . . . , Yn. Budeme testovat hypotézu: H0 : Yt = t t WN(0, 2 ) Předpokládejme, že n = 2m + 1 je liché číslo. (Je-li n sudé, obvykle se vynechá první člen jakožto časově nejvzdálenější od současnosti). Uvažujme hodnoty periodogramu v bodech k = 2 n k k = 1, . . . , m. Vypočítejme hodnoty periodogramu In(k) k = 1, . . . , m, srovnejme je podle velikosti a označme V1, . . . , Vm. Položme (tzv. Fisherova statistika) W = V1 V1 + + Vm , která nabývá hodnot mezi nulou a jedničkou. Budou-li všechny veličiny téměř stejné, bude hodnota W blízká číslu 1 m . Bude-li naopak veličina V1 nabývat velmi vysokých hodnot ve srovnání s ostatními veličinami V2, . . . , Vm, bude hodnota W blízká číslu 1 . Je tedy vidět, že velké hodnoty (které jsou blízké 1) budou tvořit kritický obor naší hypotézy proti alternativě H1 : Yt = p j=1 j cos(tj - j) + t t WN(0, 2 ) R. A. Fisher odvodil (viz Anděl, J.: Statistická analýza časových řad, Praha 1976, str. 79- 86) distribuční funkci statistiky W za platnosti hypotézy H0 : (za předpokladu, že uvažujeme gaussovský bílý šum) 1 - FW |H0 (x) = P(W > x|H0) = m(1 - x)m-1 - m 2 (1 - 2x)m-1 + , kde 0 < x < 1 a na pravé straně sčítáme tak dlouho, dokud jsou členy (1 - kx) kladné, což lze také zapsat takto P(W > x|H0) = k=1,2,... m k [max(0, 1 - kx)]m-1 = k=1,2,... m k [(1 - kx)+]m-1 . 4.2. PERIODOGRAM 99 Pak hypotézu H0 zamítáme (na hladině významnosti ), pokud 1 - FW |H0 (w) = P(W > w|H0) = W |H0 , kde w je skutečná hodnota Fisherovy statistiky při daných hodnotách konkrétní časové řady a W |H0 je tzv p-value. V případě, že pomocí Fisherova testu zjistíme signifikantní periodicitu určité frekvence, je na místě otázka, jak statisticky testovat případné další periodicity o jiných frekvencích. Whittle doporučil, aby se v případě významnosti největší hodnoty periodogramu V1 tato hodnota vynechala. Dále pak na základě veličin V2, . . . , Vm položíme W(2) = V2 V2 + + Vm a stanovíme P(W(2) > w(2) ) podle stejného vzorce, kde místo m dosadíme m - 1. Vyjde-li i tato druhá největší hodnota významná, opět se vynechá a m se zmenší o další jedničku. Když takto stanovíme frekvence 1, . . . , p, získáme model se známými frekvencemi a zbylé neznámé koeficienty odhadneme metodou nejmenších čtverců. 4.2.2 Siegelův test V praxi se ukázalo, že při platnosti alternativní hypotézy s p > 1 (tzn., že se uplatňuje více než jedna periodicita) nezamítá Fisherův test nulovou hypotézu s příliš velkou pravděpodobností. Jinými slovy lze říci, že Fisherův test nemá při alternativní hypotéze pro p > 1 takovou sílu jako v případě, kdy p = 1. (Tehdy je jeho síla přijatelná a dokonce v jistém smyslu optimální). Proto byly hledány modifikace, které tento nedostatek odstraňují. Uvedeme tzv. Siegelův test. Zde se místo statistiky W používá statistika T tvaru T = m j=1 In(j) m i=1 In(i) - gF + , kde (x)+ = max(x, 0), gF je 100(1 - )% kritická hodnota Fisherova testu (tj. P(W > gF |H0) = ) a 0 < < 1 je předem zvolená konstanta (doporučuje se volit = 0.6). Nulovou hypotézu pak zamítáme, pokud T > t, kde t je kritická hodnota tohoto testu. Kritické hodnoty bývají tabelovány pro různá a m. Jako významné jsou uznány pouze ty periodicity, jejichž odpovídající sčítance přispěly do celkové hodnoty testovací statistiky T. 4.2. PERIODOGRAM 100 20 40 60 80 100 120 -10 0 10 20 30 40 50 60 70 80 yt = a + bt + Ccos(t-()) + t ; =2/T; ()=2(+/T)=2/T(T+); t ~ WN(0,2 ) (1): a=30; b=0; C=20; T=40; =0; = 0; =10; WN = N(0, 2 ) 0 20 40 60 80 100 120 -5 0 5 10 15 Realne Fourierovy koeficienty ak 0 20 40 60 80 100 120 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Imaginarni Fourierovy koeficienty bk 0 20 40 60 80 100 120 0 5 10 15 20 25 AMPLITUDOVE SPEKTRUM A=2|z k |; z k =a k +i.b k stredova symetrie 0 20 40 60 80 100 120 -200 -150 -100 -50 0 50 100 150 200 FAZOVE SPEKTRUM k =arctg(b k /a k ) stredova antisymetrie stupen 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 x 10 4 PERIODOGRAM Ik =(2/n)|zk | 2 ; zk =ak +i.bk Fisher's test 40.0 Obrázek 4.4: Příklad periodické časové řady. 4.2. PERIODOGRAM 101 20 40 60 80 100 120 0 20 40 60 80 100 120 yt = a + bt + Ccos(t-()) + t ; =2/T; ()=2(+/T)=2/T(T+); t ~ WN(0,2 ) (2): a=30; b=0.5; C=20; T=40; =0; = 0; =10; WN = N(0, 2 ) 0 20 40 60 80 100 120 -5 0 5 10 Realne Fourierovy koeficienty ak 0 20 40 60 80 100 120 -10 -8 -6 -4 -2 0 2 4 6 8 10 Imaginarni Fourierovy koeficienty bk 0 20 40 60 80 100 120 0 2 4 6 8 10 12 14 16 18 20 AMPLITUDOVE SPEKTRUM A=2|z k |; z k =a k +i.b k stredova symetrie 0 20 40 60 80 100 120 -200 -150 -100 -50 0 50 100 150 200 FAZOVE SPEKTRUM k =arctg(b k /a k ) stredova antisymetrie stupen 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 x 10 4 PERIODOGRAM Ik =(2/n)|zk | 2 ; zk =ak +i.bk Fisher's test 40.0 120.0 60.0 Obrázek 4.5: Příklad periodické časové řady spolu s trendem. Na obrázku periodogramu je vidět, jak přítomnost trendu ovlivní velikost nízkých frekvencí. Proto se doporučuje u nesta- cionárních časových řad trend z řady nejprve odstranit. 4.3. Odhady spektrální hustoty 102 4.3 Odhady spektrální hustoty 4.3.1 Neparametrické odhady spektrální hustoty (Window Spectral Estimation) Neparametrické odhady spektrální hustoty centrované stacionární náhodné posloupnosti {Yt, t Z} jsou založeny na zlepšení vlastností periodogramu. Periodogram je empirickým odhadem spektrální hustoty, který je asymptoticky nestranný, avšak nekonzistentní. Připomeňme, že platí (viz lemma 4.2.4) In() = 1 2 C0 + 2 n-1 k=-(n-1) C k cos k = 1 2 n-1 k=-(n-1) C k cos k, neboť pro k = 0, 1, 2, . . . , (n - 1) platí C k = C -k, přičemž C k = 1 n n-k t=1 YtYt+k k = 0, . . . , n - 1. Vzhledem k tomu, že platí cos k = 1 2 eik + e-ik , upravujme In() = 1 2 n-1 k=-(n-1) C k 1 2 eik + e-ik = 1 2 C 0 + 1 2 -1 k=-(n-1) C k 1 2 eik + -1 k=-(n-1) C k 1 2 e-ik + n-1 k=1 C k 1 2 eik + n-1 k=1 C k 1 2 e-ik = 1 2 C 0 + 1 2 n-1 s=1 C -s 1 2 e-is + -1 k=-(n-1) C k 1 2 e-ik + -1 s=-(n-1) C -s 1 2 e-is + n-1 k=1 C k 1 2 e-ik = 1 2 C 0 + 1 2 -1 k=-(n-1) 1 2 (C -k + C k) 2C k =2C -k e-ik + n-1 k=1 1 2 (C -k + C k) 2C k e-ik = 1 2 n-1 k=-(n-1) C ke-ik . 4.3. Odhady spektrální hustoty 103 Uvědomíme-li si, že periodogram jakožto odhad spektrální hustoty, je založen na všech možných odhadech autokovariační funkce v bodech k = 0, 1, 2, . . . , (n - 1), tj. C 0 = 1 n Y 2 1 + + Y 2 n n členů C 1 = C -1 = 1 n (Y1Y2 + + Yn-1Yn + Y3Yn) n-1 členů ... C n-3 = C -(n-3) = 1 n (Y1Yn-2 + Y2Yn-1 + Y3Yn) 3 členy C n-2 = C -(n-2) = 1 n (Y1Yn-1 + Y2Yn) 2 členy C n-1 = C -(n-1) = 1 n Y1Yn 1 člen a tedy je založen i na velmi málo kvalitních odhadech. K určitému zlepšení jistě dojde, pokud budeme používat jen m n nejkvalitnějších odhadů. Mluvíme pak o prostém useknutém periodogramu ^fn() = 1 2 m k=-m C k cos k = 1 2 m k=-m C k e-ik , což lze také zapsat takto ^fn() = 1 2 n-1 k=-(n-1) w(k)C k cos k = 1 2 n-1 k=-(n-1) w(k)C ke-ik , kde w(k) = 1 |k| m 0 |k| > m . Označme Fourierovu transformaci funkce w(k) W() = 1 2 k=- w(k)e-ik = 1 2 m k=-m e-ik a položme s = k + m + 1, pak k = s - m - 1 a W() = 1 2 2m+1 s=1 e-i(s-m-1) = 1 2 ei(m+1) 2m+1 s=1 e-is = 1 2 eim 1 - e-i(2m+1) 1 - e-i = 1 2 eim e-i 2m+1 2 ei 2m+1 2 - e-i 2m+1 2 e-i 1 2 ei 1 2 - e-i 1 2 = 1 2 sin(m + 1 2 ) sin 1 2 = Dm(), kde Dm() je tzv. Dirichletovo jádro. 4.3. Odhady spektrální hustoty 104 Vzhledem k tomu, že lze psát In() = 1 2 n-1 k=-(n-1) C ke-ik , pak pomocí inverzní Fourierovy transformace platí C k = - In()e-ik d . Počítejme postupně ^fn() = 1 2 n-1 k=-(n-1) w(k)C ke-ik = 1 2 n-1 k=-(n-1) w(k) - In()eik d e-ik = - In() 1 2 n-1 k=-(n-1) w(k)e-ik(-) W (-) d = - In()W( - )d . Jde o tzv. vyhlazený periodogram (smoothed periodogram). Funkce W() se nazývá spek- trální okénko (spectral window). Tato funkce má do jisté míry aproximovat Diracovu funkci a platí pro ni - W()d = 1. Takto počítat odhad spektrální hustoty by však bylo (vzhledem k málo hladkému průběhu periodogramu) nepohodlné, proto se obvykle odhad počítá podle vzorce ^fn() = 1 2 n-1 k=-(n-1) w(k)C ke-ik , přičemž inverzní Fourierovy transformace w(k) = - W()eik d k = 0, 1, 2, . . . (n - 1) se nazývá korelační okénko (covariance lag window, nebo time-domaing window). Typickými korelačními okénky jsou tzv. useknutá okénka, pro která existuje takové přiro- zené číslo m (bod useknutí, truncation point) tak, že w(k) = 0 pro |k| > m (m se obvykle volí v rozmezí od n 6 do n 5 ). 4.3. Odhady spektrální hustoty 105 Příklady korelačních a spektrálních okének Prostý useknutý odhad w(k) = 1 0 < |k| m 0 |k| > m W() = 1 2 sin(m + 1 2 ) sin 1 2 -6 -4 -2 0 2 4 6 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Lag window w(k) -5 0 5 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Spectral window W()-Dirichlet kernel -6/7 -4/7 -2/7 0 2/7 4/7 6/7 Obrázek 4.6: Korelační a spektrální okénko v případě prostého useknutého periodogramu. Bartletovo okénko w(k) = 1 - |k| m 0 < |k| m 0 |k| > m W() = 1 2m sin2 m 2 sin2 2 = Fm() W() je v tomto případě Fejérovým jádrem. Parzenovo okénko w(k) = 1 - 6 k m 2 + 6 |k| m 3 |k| < m 2 2 1 - |k| m 3 m 2 < |k| m 0 |k| > m W() = 3 8m3 sin m 4 1 2 sin 2 4 1 - 2 3 sin2 2 kde m je nějaké sudé číslo. Obecné Tukeovo okénko w(k) = 1 - 2a + 2a cos k m |k| m 0 |k| > m W() = aDm - m + (1 - 2a)Dm() + aDm + m kde a (0, 1 4 . Tukey-Hanningovo okénko w(k) = 1 2 1 + cos k m |k| m 0 |k| > m W() = 1 4 Dm - m + 1 2 Dm() + 1 4 Dm + m 4.3. Odhady spektrální hustoty 106 Tukey-Hammingovo okénko w(k) = 0.54 + 0.46 cos k m |k| m 0 |k| > m W() = 0.23Dm - m + 0.54Dm() + 0.23Dm + m Daniellovo okénko . Na závěr ještě uvedeme jedno neuseknuté korelační okénko. Mějme pro (0, ) následující spektrální okénko W() = 1 2 || < 0 || > , které je vlastně hustotou náhodné veličiny s rovnoměrně spojitým rozdělením na intervalu (-, ). Počítejme nejprve pro k = 1, 2, . . . (n - 1) odpovídající korelační okénko: w(k) = - W()eik d = - 1 2 eik d = 1 2 eik ik - = 1 k 1 2i eik - e-ik sin k = sin k k . Pro k = 0 máme w(0) = - W()d = - 1 2 d = 1 . A celkově dostaneme wk = 1 k = 0 sin k k k = 1, 2, . . . .