Marie Forbelská Stochastické modelování jednorozměrných časových řad I 1 12 24 36 1 2 3 4 x 10 4 1999 2000 2001 Punkevní jeskyne 1 12 24 36 5000 10000 15000 20000 1999 2000 2001 Sloupsko-Sosuvske jeskyne 1 12 24 36 5000 10000 15000 20000 1999 2000 2001 Katerinska jeskyne 1 12 24 36 5000 10000 15000 20000 1999 2000 2001 Balcarka Návštěvnost v Moravském krasu od ledna 1999 do prosince 2001 Brno 2OO7 Příro dovědecká fa kulta Masa rykova Unive rzita Ústav matematiky a statistiky iii Obsah Seznam použitých symbolů 1 1 Náhodné procesy 3 1.1 Definice náhodného procesu . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Příklady náhodných procesů . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Stochastické procesy druhého řádu . . . . . . . . . . . . . . . . . . 6 1.3.1 Striktní a slabá stacionarita . . . . . . . . . . . . . . . . . . 6 1.3.2 Příklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4 Kanonické rozklady náhodných procesů . . . . . . . . . . . . . . . 14 1.4.1 Elementární procesy a kanonické rozklady . . . . . . . . . . 14 1.4.2 Příklady kanonických rozkladů . . . . . . . . . . . . . . . . 16 1.5 Vlastnosti autokovariační funkce . . . . . . . . . . . . . . . . . . . 19 1.6 Spojitost, derivace a integrál procesu . . . . . . . . . . . . . . . . . 22 1.6.1 Spojitost náhodného procesu . . . . . . . . . . . . . . . . . 22 1.6.2 Derivace náhodného procesu . . . . . . . . . . . . . . . . . 24 1.6.3 Integrál náhodného procesu . . . . . . . . . . . . . . . . . . 24 1.7 Spektrální rozklad autokovariančních funkcí stacionárních procesů 26 1.7.1 Herglotzova a Bochnerova věta . . . . . . . . . . . . . . . . 26 1.7.2 Příklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 1.8 Hilbertovy prostory spjaté s procesy druhého řádu . . . . . . . . . 37 1.8.1 Základní metrické a topologické pojmy . . . . . . . . . . . . 37 1.8.2 Hilbertův prostor náhodných veličin druhého řádu . . . . . 40 1.8.3 Predikce v případě normálně rozdělených náhodných veličin. 43 1.9 Odhady středních hodnot a autokovariancí . . . . . . . . . . . . . . 44 1.9.1 Odhady střední hodnoty . . . . . . . . . . . . . . . . . . . . 45 1.9.2 Odhady autokovarianční a autokorelační funkce . . . . . . . 50 2 Analýza časových řad 51 2.1 Základní přístupy k analýze časových řad . . . . . . . . . . . . . . 51 2.1.1 Klasická dekompozice . . . . . . . . . . . . . . . . . . . . . 51 2.1.2 Box-Jenkinsonova metodologie . . . . . . . . . . . . . . . . 52 2.1.3 Spektrální analýza časových řad . . . . . . . . . . . . . . . 52 2.1.4 Lineární dynamické modely . . . . . . . . . . . . . . . . . . 52 iv 2.2 Klasická dekompozice časových řad . . . . . . . . . . . . . . . . . . 53 2.2.1 Obecné lineární regresní modely a metoda nejmenších čtverců 53 2.2.2 Rozšířený lineární regresní model a vážená metoda nejmenších čtverců . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.2.3 Modelování trendu . . . . . . . . . . . . . . . . . . . . . . . 57 2.2.4 Sezónnost . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 2.2.5 Analýza reziduální (náhodné) složky . . . . . . . . . . . . . 97 3 Spektrální analýza jednorozměrných časových řad 105 3.1 Úvod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.2 Periodogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.3 Metoda skrytých period. . . . . . . . . . . . . . . . . . . . . . . . . 116 3.3.1 Regresní model s periodickým trendem . . . . . . . . . . . . 116 3.3.2 Test R. A. Fishera . . . . . . . . . . . . . . . . . . . . . . . 118 3.3.3 Siegelův test . . . . . . . . . . . . . . . . . . . . . . . . . . 120 3.3.4 Odhady neznámých parametrů modelu a periodogram . . . 123 3.4 Odhady spektrální hustoty . . . . . . . . . . . . . . . . . . . . . . . 126 3.4.1 Neparametrické odhady spektrální hustoty (Window Spectral Estimation) . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Literatura 133 1 Seznam použitých symbolů {x : } množina prvků x, které splňují podmínku N množina přirozených čísel R množina reálných čísel R+ množina kladných reálných čísel R+ 0 množina nezáporných reálných čísel Rm reálný m-rozměrný euklidovský prostor Sm povrch jednotkové koule v Rm a, b, . . . reálná čísla a, b, . . . vektory reálných čísel (sloupcové) a euklidovská norma vektoru a, tj. a = aa A, B, . . . matice a A norma vektoru a vzhledem k dané pozitivně definitní matici A, tj. a A = aAa In jednotková matice typu (n × n) 0 nulový vektor nebo nulová matice 1n n-rozměrný sloupcový vektor samých jedniček rk(A) hodnost matice tr(A) stopa matice A transponovaná matice A-1 inverzní matice Apseudoinverzní matice |A|, det(A) determinant matice A B Kroneckerův součin matic A a B vec(A) operátor, který ze sloupců matice A vytváří vektor A B C D bloková matice f(x)dx Lebesgueův integrál f(x)dF(x) Lebesgue-Stieltjesův integrál (x) gamma funkce B(a, b) beta funkce ln x přirozený logaritmus se základem e prostor elementárních jevů A -algebra 2 Bm borelovská -algebra v m-rozměrném euklidovském prostoru f : (Rm , Bm ) R funkce f : Rm R je Bm -měřitelná IA indikátor množiny A A B sjednocení množin A B průnik množin A - B rozdíl množin prázdná množina P(A) pravděpodobnost jevu A (, A, P) pravděpodobnostní prostor X, Y, . . . náhodné veličiny x, y, . . . realizace náhodných veličin X, Y, . . . náhodný vektor x, y, . . . realizace náhodných vektorů f(x), g(x), . . . hustoty pravděpodobnosti p(x), q(x), . . . pravděpodobnostní funkce F(x), G(x), . . . distribuční funkce F-1 (x), G-1 (x), . . . kvantilové funkce (inverzní distribuční funkce) E střední hodnota C(X, Y ) kovariance náhodných veličin X a Y DX rozptyl náhodné veličiny X C(X, Y) kovarianční matice náhodných vektorů X a Y DX varianční matice náhodného vektoru X X N(, 2 ) náhodná veličina X se řídí normálním rozdělením se střední hodnotou a rozptylem 2 (x) hustota pravděpodobnosti standardizovaného normálního rozdělení N(0, 1) (x) distribuční funkce standardizovaného normálního rozdělení N(0, 1) Rs(a, b) rozdělení rovnoměrně spojité na intervalu (a, b) (a < b) Rs(Sm ) rovnoměrné rozdělení na povrchu m-rozměrné jednotkové koule (a) Gamma rozdělení s parametrem a (a, b) Beta rozdělení s parametry a a b 2 (n) 2 -rozdělení o n stupních volnosti Nm(, ) m-rozměrné normální rozdělení se střední hodnotou a kovarianční maticí 3 Kapitola 1 Náhodné procesy Ú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, 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á 4 se zmíněnými procesy zabývá, se také nazývá statistickou dynamikou. Matematicko-statistické metody zpracování pozorování těchto náhodných procesů mají velmi široké uplatnění v praxi. Jsou založené na teorii pravděpodobnosti náhodných procesů, kterou se nyní budeme zabývat. 1.1 Definice náhodného procesu Definice 1.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 pro t T je X(, t) borelovsky měřitelná funkce vzhledem k A (tj. pro B B a pro t T platí 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)) nazýváme rozdělení pravděpodobností náhodného procesu X(, t). Poznámka 1.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 1.1.3. Pokud indexová množina T = Z = {0, 1, 2, . . .} nebo T Z, mluvíme o procesu s diskrétním časem nebo o časové řadě či náhodné posloupnosti. Pokud indexová množina T = t1, t2 , kde - t1 < t2 +, říkáme, že {Xt, t T } je náhodný proces se spojitým časem. 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. Nabý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 n rozměrného náhodného vektoru X = (Xt1 , . . . , Xtn ) , kde t1, . . . , tn jsou libovolné body z množiny T . Definice 1.1.4. Nechť T n je množina všech vektorů T n = {t = (t1, . . . , tn) : t1 t2 tn; ti T ; i = 1, . . . , n}. 5 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) T n 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, 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 dané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 1.1.5. 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. 6 1.2 Příklady náhodných procesů Příklad 1.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 1.2.2. Binární proces {Xt, t = 1, 2, . . .} je posloupnost nezávislých alternativních náhodných proměnných: Xt A 1 2 . Příklad 1.2.3. Náhodná procházka {Xt, t = 1, 2, . . .} je definována 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 konstantním rozptylem, tj. t L(0, 2 ). Příklad 1.2.4. Markovovův proces {Xt, t = 1, 2, . . .} je definován takto P(a < Xt+1 b|X1 = x1, . . . , Xt = xt) = P(a < Xt+1 b|Xt = xt), 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. 1.3 Stochastické procesy druhého řádu 1.3.1 Striktní a slabá stacionarita Definice 1.3.1. Řekneme, že náhodný proces {Xt, t T } je striktně stacionární, jestliže pro t = (t1, . . . , tn) T n a pro = (t1 + h, . . . , tn + h) T n 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. 7 Definice 1.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. Definice 1.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 1.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 1.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 1.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 1.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ů. 8 Poznámka 1.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 1.3.9. Náhodný proces {Xt, t T } se nazývá (slabě) stacionární, je-li kovariančně stacioná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 1.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 1.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 pro- cesu Definice 1.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), 9 značíme t WN(0, 2 ). 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 1.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ý. Definice 1.3.14. Náhodný proces {Xt, t T } se nazývá gaussovským (normálním), jestliže pro každé přirozené 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 1.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í. Definice 1.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ů. 10 1.3.2 Příklady Příklad 1.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), x R. 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 1.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 1.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í konstantní. 11 Příklad 1.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, -, . 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čítejme 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 + h) - t] = cos h = (h) Tedy (t + h, t) nezávisí na t, proto proces je (slabě) stacionární. Náhodná procházka. Příklad 1.3.21. 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 12 tj. proces je konstantní ve střední hodnotě. Dále počítejme autokovarianční funkci: (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í. MA(1) proces. Příklad 1.3.22. 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 . Markovovův skokový proces vzniku a zániku. Příklad 1.3.23. Mějme náhodný proces {Xt, t R}, pro který platí P(Xt = 0) = + P(Xt = 1) = + > 0, > 0 P(Xt = x|Xs = y) = pyx(t - s) - < s t < , 13 přičemž p00(t) = 1 - p01(h) = + + + e-(+)h h 0 p11(t) = 1 - p10(h) = + + + e-(+)h h 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(Xs = 1) P(Xt = 1|Xs = 1) p11(t-s) = + + + + e-(+)(t-s) 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 Autokorelační funkce 14 1.4 Kanonické rozklady náhodných procesů 1.4.1 Elementární procesy a kanonické rozklady Definice 1.4.1. Nechť je dán pravděpodobnostní prostor (, A, P) a indexová množina T R. Náhodný proces {Xt, t T } nazveme elementárním, pokud lze vyjádřit ve tvaru Xt = Uf(t) , kde U je náhodná veličina definovaná na (, A, P) a f(t) je reálná funkce definovaná na daném T R. Příklady: * Xt = U sin t, * Xt = Ut2 . Určeme charakteristiky elementárních náhodných procesů s konečnými druhými momenty (tj. EU2 < ) a označme EU = U DU = 2 U . Postupně počítejme EXt = f(t)EU = f(t)U (t, t + h) = C(Xt, Xt+h) = E(Xt - EXt)(Xt+h - EXt+h) = E(XtXt+h) - EXtEXt+h = E(U2 f(t)f(t + h)) - (f(t)EU)(f(t + h)EU) = f(t)f(t + h)EU2 - f(t)f(t + h)(EU)2 = f(t)f(t + h) (EU2 - (EU)2 ) DU = f(t)f(t + h)DU = f(t)f(t + h)2 U 15 Definice 1.4.2. Nechť je dán pravděpodobnostní prostor (, A, P), indexová množina T R, spočetně mnoho náhodných veličin U1, U2, . . . definovaných na (, A, P) a spočetně mnoho reálných funkcí f1(t), f2(t), . . . definovaných na daném T R. Kanonickým rozkladem náhodného procesu {Xt, t T } nazveme pak roz- klad Xt = X(t) + i=1 Uifi(t) , kde U1, U2, . . . jsou vzájemně nezávislé centrované náhodné veličiny s konečnými druhými momenty a X(t) je reálná funkce definovaná na daném T R. Označíme-li pro i = 1, 2, . . . DUi = 2 Ui , pak lze snadno spočítat: EXt = E(X(t) + i=1 Uifi(t)) = X (t) + i=1 fi(t) EUi =0 = X(t) (t, t + h) = C(Xt, Xt+h) = E(Xt - EXt)(Xt+h - EXt+h) = E(XtXt+h) - EXtEXt+h = i=1 E(U2 i fi(t)fi(t + h)) - i=1 (fi(t)EUi)(fi(t + h)EUi) = i=1 [fi(t)fi(t + h)EU2 i - f(t)f(t + h)(EUi)2 ] = i=1 fi(t)fi(t + h) (EU2 i - (EUi)2 ) DUi=2 Ui = i=1 fi(t)fi(t + h)2 Ui 16 1.4.2 Příklady kanonických rozkladů Příklad 1.4.3. Mějme nezávislé náhodné veličiny Ui N(0, 2 ), i = 1, 2 a pro libovolné reálné číslo mějme Xt = U1 cos t + U2 sin t . Jaké vlastnosti má tento náhodný proces? EXt = E(U1 cos t + U2 sin t) = cos tEU1 + sin tEU1 = 0 (t, t + h) = C(Xt, Xt+h) = E(Xt - EXt)(Xt+h - EXt+h) = E(XtXt+h) = E(U1 cos t + U2 sin t)(U1 cos (t + h) + U2 sin (t + h)) = cos t cos (t + h)EU1 + sin t sin (t + h)EU2 = 2 [cos t cos (t + h) + sin t sin (t + h)] cos (t+h-t) = 2 cos h = (h) DXt = (0) = 2 (h) = (h) (0) = cos h Proces je slabě stacionární. Rozšířený Poissonův proces. Příklad 1.4.4. Mějme náhodný proces {Xt, t R}, který splňuje tyto vlastnosti (i) X0 = 0 (ii) pro s, t T, s t mají náhodné veličiny Xt - Xs Poissonovo rozdělení s parametrem (t - s), kde > 0, tj. Xt - Xs Po((t - s)) . (iii) a Xt2 -Xt1 , . . . , Xtn -Xtn-1 jsou pro každé přirozené n a t1 t2 . . . tn sdruženě nezávislé náhodné veličiny. Z vlastností (i) a (ii) vyplývá, že pro t T, t 0 má náhodná veličina Xt Po(t) a pro t T, t < 0 platí, že náhodná veličina Xt Po(-t), 17 celkově Xt Po(|t|). Podle známých vlastností Poissonova rozdělení je střední hodnota tohoto náhodného procesu rovna EXt = |t| , tedy tento proces není stacionární ve střední hodnotě. Protože rozptyl náhodné veličiny s Poissonovým rozdělením je roven střední hodnotě, tedy tento proces má konečné druhé momenty a platí DXt = |t| . Pro výpočet autokovariační funkce nejprve předpokládejme, že 0 s t. Pak platí (s, t) = C(Xs, Xt) = C( Xs - 0 =Xs-X0 , (Xt - Xs) + Xs) = C(Xs - X0, Xt - Xs) =0(nekorel.) + C(Xs, Xs) =DXs = s Za předpokladu, že s < 0, t > 0, pak z vlastností (i) a (ii) vyplývá, že (s, t) = C(Xs, Xt) = 0. Ostatní případy dostaneme obdobně. Celkově tedy (s, t) = min(|s|, |t|) st > 0 0 st < 0 . Tento proces druhého řádu není stacionárním procesem. V následujícím příkladu však ukážeme, že stochastický proces, který vznikne z jeho jednokrokových přírůstků, je už stacionárním procesem druhého řádu. Příklad 1.4.5. Proces jednokrokových přírůstků rozšířeného Poissonova procesu. Nechť {Xt, t T } je rozšířený Poissonův proces z předchozího příkladu a nechť Yt = Xt+1 - Xt t R. Určeme střední hodnotu a autokovariační funkci náhodného procesu {Yt, t R}. Protože EXt = |t| , EYt = E(Xt+1 - Xt) = (t + 1) - t = , 18 proces je stacionární ve střední hodnotě. Pro výpočet autokovariační funkce nejprve předpokládejme, že |t-s| 1. Pak náhodné veličiny Ys = Xs+1 - Xs a Yt = Xt+1 - Xt jsou podle (iii) stochasticky nezávislé, proto platí Y (s, t) = C(Ys, Yt) = 0 pro |t - s| 1. Dále předpokládejme, že s t s + 1. Potom postupně dostaneme (s, t) = C(Ys, Yt) = C(Xs+1 - Xs, Xt+1 - Xt) = C((Xs+1 - Xt) + (Xt - Xs), (Xt+1 - Xs+1) + (Xs+1 - Xt)) = C(Xs+1 - Xt, Xt+1 - Xs+1) =0 + C(Xt - Xs, Xt+1 - Xs+1) =0 + C(Xs+1 - Xt, Xs+1 - Xt) =0 +C(Xt - Xs, Xs+1 - Xt) = (s + 1 - t). Použitím symetrie autokovarianční funkce nakonec dostaneme Y (s, t) = (1 - |t - s|) |t - s| < 1 0 |t - s| 1 . Tento stochastický proces je procesem druhého řádu se střední hodnotou EYt = a autokorelační funkcí (h) = 1 - |h| |h| < 1 0 |h| 1 . 19 1.5 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 1.5.1. Nechť f(s, t) je funkce dvou proměnných definovaná na T ×T . Říkáme, že f je pozitivně 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)cj ck 0 . (1.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 . (1.2) 20 Ze souvislosti bude vždy patrno, zda jde o definici ve smyslu (1.1) nebo (1.2). Věta 1.5.2. Autokovarianční funkce (s, t) je pozitivně semidefinitní. 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 k=1 ck Xtk = n j=1 n k=1 cj ckE(Xtj Xtk ) = n j=1 n k=1 cjck(tj, tk). Věta 1.5.3. Autokovarianční funkce (s, t) je hermitovsky symetrická, tj. pro s, t T platí (s, t) = (t, s) Důkaz. Snadno dokážeme, protože (s, t) = E(Xs Xt) = E(Xt Xs) = (t, s). Věta 1.5.4. 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. Věta 1.5.5. 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á je důsledkem Schwarzovy nerovnosti. 21 Lemma 1.5.6. Součet dvou pozitivně semidefinitních hermitovsky symetrických funkcí je opět funkce pozitivně 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 k=1 cj ckf(tj, tk) = n j=1 n k=1 cj ckf1(tj, tk) + n j=1 n k=1 cjckf2(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é. Věta 1.5.7. Součet dvou autokovariačních funkcí je opět autokovarianční funkcí. Důkaz. Důkaz plyne z lemmatu 1.5.6 a z vět 1.5.2, 1.5.3 a 1.5.4. Věta 1.5.8. Reálná část autokovarianční funkce je též autokovarianční funkcí. Imaginární část je autokovarianč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) = C(Zs, Zt) = E(Zs - EZs)(Zt - EZt). Bez újmy na obecnosti budeme předpokládat, že náhodný proces má nulovou střední hodnotu, tj. 0 = EZt = E(Xt + iYt) = EXt + EYt, což implikuje, že EXt = EYt = 0. Počítejme Z (s, t) = EZs Zt = E(Xs + iYs)(Xt - iYt) = EXsXt + EYsYt + i(EYsXt - EXsYt) 22 Reálná část Z(s, t) je rovna Re(Z(s, t)) = EXsXt + EYsYt = X(s, t) + Y (s, t). Je tedy rovna součtu autokovariační funkce procesu {Xt, t T } a autokovariační funkce procesu {Yt, t T }. Podle věty 1.5.7 je autokovariančí funkcí. Imaginární část Z(s, t) je rovna Im(Z (s, t)) = EYsXt - EXsYt. Připomeňme, že pro libovolnou autokovarianční funkce (s, t) musí platit: (i) (s, s) 0 (ii) 0 |(s, t)| (s, s) (t, t). V bodech s = t dostaneme Im(Z(s, s)) = EYsXs - EXsYs = 0. Druhá nerovnost však je splněna 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. 1.6 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ů. 1.6.1 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ší konvergence podle kvadratického středu. 23 Definice 1.6.1. Ř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 1.6.2. Z teorie pravděpodobnosti je známo, že konvergence podle kvadratického středu implikuje konvergenci podle pravděpodobnosti. Věta 1.6.3 (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í pro s, t, s0, t0 T 0 |(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. 24 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 pro s t, tj. konverguje podle kvadratického středu. 1.6.2 Derivace náhodného procesu Derivaci náhodného procesu budeme definovat obdobně, jako se definuje derivace funkce. Definice 1.6.4. Řekneme, že náhodný proces {Xt, t T } má v bodě t0 T derivaci X t0 , jestliže platí l.i.m. h0 Xt0+h - Xt0 h = X t0 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 1.6.3 Integrál náhodného procesu Definice 1.6.5. Nechť T = (a, b) je konečný nedegenerovaný interval. Nechť [T 1, T 2] T , kde T1 < T2. Pomocí dělení Dn intervalu [T 1, T 2] takového, že T1 = t0 < t1 < < tn = T2, utvoříme součet In = n-1 k=0 Xtk (tk+1 - tk). 25 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 [T 1, T 2]. Píšeme pak I = T2 T1 Xt dt. Věta 1.6.6 (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 Riemannů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. V praxi samozřejmě není nutno počítat T2 T1 T2 T1 (s, t) ds dt, abychom se přesvědčili, že existuje 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 [T 1, T 2] × [T 1, T 2]. Poznámka 1.6.7. 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. 26 1.7 Spektrální rozklad autokovariančních funkcí stacionárních procesů 1.7.1 Herglotzova a Bochnerova věta 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 1.7.1 (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á. Důkaz. Protože (t) je pozitivně semidefinitní funkce, platí 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 vztah n j=1 n k=1 (tj - tk)cj ck 0. (1.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 (1.3) platí fn() 0 pro - . 1 2 3 4 5 1 2 3 4 5 (n - 1) k rát(n - 1) k rát s = 1 s = - 1 (n - 2) k rát (n - 2) k rát s = 2 s = - 2 (n - 3) k rát (n - 3) k rát s = 3 s = - 3 n krát s = 0 k j 27 Po substituci j - k = s, kde s {-n + 1, . . . , n - 1} dostáváme vyjádření (viz obrázek) 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 pro - , 0 jinak. (1.4) Uvědomíme-li si, že platí - eikx dx = 2 k = 0, 0 k Z, k = 0, (1.5) což lze ověřit třeba i rozkladem eikx = cos kx + i sin kx, počítejme pro t Z - 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. (1.6) Protože fn(x) je nezáporná, je Fn() neklesající funkce. Je jasné, že Fn(-) = 0. Dále s využitím vztahu (1.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) . 28 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() . 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). Věta 1.7.2 (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á. 29 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. (1.7) Dosadíme-li do (1.7) tj = j N , kde N N cj = e-ij , kde - , dostaneme (opět s využitím substituce j - k = s, viz předchozí věta), ž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 (1.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) . Podle 1. Hellyho věty lze z posloupnosti funkcí FN 1 (), FN 2 (), . . . vybrat podposloupnost {FN nk ()}, která konverguje k nějaké neklesající zprava spojité funkci FN (), pro kterou platí FN (- - 0) = 0 a FN ( + 0) = FN () = (0). 30 Položme N (t) = N -N eit dFN N . (1.8) Z 2. Hellyho věty plyne, že pro každé celé k platí N k N = N -N ei k N dFN N = subst. x = N = - eikx dFN (x) = lim nj - eikx dFN nj (x) = lim nj - eikx fN nj (x)dx (1.9) = 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 . (1.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 (1.9) máme (t) = lim N k N = lim N N k N . (1.11) Ze zřejmé rovnosti N (t) = N k N + N (t) - N k N plyne vzhledem k (1.10) lim N N (t) = (t) + lim N N (t) - N k N (1.12) Dosadíme-li z (1.8) a z první části (1.9), máme N (t) - N k N = N -N ei k N ei(t- k N ) - 1 dFN N Nechť = t - k N < 1 N . Z (1.10) plyne 0 N < 1. 31 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 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 = subst. x = N = 2(0) (1 - cos Nx) dFN (x). Z nerovnosti cos b cos pro - , 0 b < 1, plyne 1 - cos N <1 x 1 - cos x, takže N (t) - N k N 2 2(0) (1 - cos Nx) dFN (x) 2(0) (1 - cos x) dFN (x) = 2(0) (0) - reálná část 1 N , neboť 1 N = N 1 N = - eix dFN (x) = (cos x + i sin x)dFN (x). 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 (1.12) tedy platí (t) = lim N N (t). (1.13) 32 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 (1.8) můžeme napsat jako N (t) = - eit dFN (). Jelikož FN () má všechny vlastnosti distribuční funkce, je {N (t)} posloupnost charakteristických funkcí, která ve všech bodech t podle (1.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í distribuč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. (1.14) 33 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 (1.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 resp. (t) = - eit f()d . (1.15) 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 1.7.3. K existenci spektrální hustoty stacionární náhodné posloupnosti stačí, aby pro její kovarianč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 < . Důkaz. Viz Gichman, I.I., Skorochod, A.V.: Teorija slučajnych processov. Moskva. Nauka 1971. V následujících dvou větách je zodpovězena otázka, jak vypočítat spektrální hustotu z kovarianční funkce. Věta 1.7.4. 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) (1.16) ve všech bodech spojitosti funkce f(), což je skoro všude vzhledem k Lebesqueově míře. Důkaz. Ze vzorce (1.15) vidíme, že až na normující konstantu 1 2 jsou (t) Fourierovy koeficienty 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ě 34 bodů nespojitosti (Připomeňme, že variace je difinována takto b a (f) = sup Dn n k=1 |f(xk) - f(xk-1)|, kde Dn = {a = x0 < x1 < < xn = b} je dělení intervalu a, b .) Věta 1.7.5. 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. (1.17) Důkaz. Ze vzorce (1.15) vidíme, že až na normující konstantu 1 2 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. Věta 1.7.6. 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(-) (1.18) 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 (1.15) (t) = - eit f()d = - e-it f()d = (-t). Substitucí se snadno zjistí, že pravá strana je rovna - eit f(-)d takže - eit f()d = - eit f(-)d. (1.19) 35 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 vzorec (1.19) 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 charakteristickou funkcí odtud vyplývá tvrzení věty. Pro stacionární posloupnosti je důkaz obdobný. 1.7.2 Příklady Příklad 1.7.7. 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 . Příklad 1.7.8. Zobecníme-li předcházející náhodný 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 . Na závěr tohoto odstavce ukážeme graf s několika příklady autokovariančních funkcí a jim odpovídajících spektrálních hustot. 36 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: Y t =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 1.1: 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. -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 1.2: Příklady autokovarianční funkce a spektrální hustoty. 37 1.8 Hilbertovy prostory spjaté s procesy druhého řádu 1.8.1 Základní metrické a topologické pojmy 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); 38 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 (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 . 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 pro- storu 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) 39 (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 M, tzv. lineární obal množiny M. 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 (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 n N a 1, . . . , n C taková, ž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, . . . 40 1.8.2 Hilbertův prostor náhodných veličin druhého řádu 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. komplexní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. Věta 1.8.1. 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., SpringerVerlag, New York, 1991, §2.10. 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ý. 41 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 1.8.2. 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 1.8.3. 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 1.8.4. 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. Výpočet podmíněných středních hodnot bývá však velmi obtížný. Z těchto důvodů se omezujeme většinou na jednodušší podprostor, a to M = sp{1, Z1, . . . , Zn} = {1, Z1, . . . , Zn}. 42 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) 43 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 veličin druhého řádu {Zt, t T }. Definice 1.8.5. 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ů. 1.8.3 Predikce v případě 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) 44 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í predikcí (ve smyslu minimální střední kvadratické odchylky) založenou na podmíněných středních hodnotách. 1.9 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, autokovarianč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í charak- teristiky. 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). 45 Příklad 1.9.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ěpodobnostní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 . Tedy náhodný proces {Y (t), t R} je stacionární proces, ale nemůžeme ho považovat za ergodický, 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. Autokovarianční funkce stacionárního procesu Y (t), t R se od autokovarianč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ů. 1.9.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 1.9.2. Odhad střední hodnoty ^ stacionárního náhodného procesu {Y (t), t [0, T ]} pomocí metody nejmenších čtverců (MNČ) je definován vztahem: ^ = arg min R T 0 (Y (t) - ) 2 dt. Poznámka 1.9.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. 46 Snadno lze odvodit, že odhad střední hodnoty pomocí MNČ je roven ^ = 1 T T 0 Y (t)dt (1.20) 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 1.9.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. (1.21) 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 = . Střední kvadratická chyba v případě nestranného odhadu je rozptylem tohoto odhadu MSE(^) = E (^ - )2 = E (^ - E^)2 = D(^). 47 Počítejme MSE(^) = E (^ - )2 = E 1 T T 0 Y (t) dt - 2 = E 1 T T 0 (Y (t) - ) dt 2 = 1 T 2 E T 0 T 0 (Y (s) - )(Y (t) - ) ds dt = 1 T 2 T 0 T 0 E [(Y (s) - )(Y (t) - )] Y (t-s)(stac.) ds = 1 T 2 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 T 2 T -T min{T,T +u} max{0,u} Y (u) dv du = 1 T 2 0 -T Y (u) T +u 0 dv du + T 0 Y (u) T u dv du = 1 T 2 0 -T Y (u)(T + u) du + T 0 Y (u)(T - u) du = 1 T 2 T -T Y (u)(T - |u|) du = 2 T 2 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 . 48 Pro další studium ergodických procesů je vhodné vyslovit následující definici: Definice 1.9.5. Řekneme, že stacionární proces {Y (t), t R} je ergodický ve střední hodnotě, pokud platí lim T D 1 T T 0 Y (t) dt = 0. (1.22) Věta 1.9.6. Nechť pro stacionární proces {Y (t), t R} s autokovarianč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ů (1.21), (1.22) a nerovnosti T 0 1 - u T Y (u)du T 0 1 - u T |Y (u)| du. Důsledek 1.9.7. Nechť lim t Y (t) = 0. Pak stacionární proces s autokovarianč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)| < . 49 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ě. 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. 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 diskretizací 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 . 50 1.9.2 Odhady autokovarianční a autokorelační funkce Odhad autokovarianč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 autokovarianč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í od- hady 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. 51 Kapitola 2 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. 2.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ě analýze. 2.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 * T rt trend odráží dlouhodobé změny v chování časových řad; 52 * 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 kolísání časových řad. Při klasické dekompozici časových řad se používají především * aditivní modely: Yt = T rt + Szt + t * multipikativní modely Yt = T rt 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. 2.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. 2.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" sinusovek 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ředevší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í. 2.1.4 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ů. 53 2.2 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 = T rt + Szt + t, multiplikativního modelu Yt = T rt Szt t. Jednotlivé složky T rt, 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 s nulovou střední hodnotou Et = 0, který je nekorelovaný, tj. C(t, s) = Ets = 0 s = t, Dt = 2 s = t. Značíme Ln(0, 2 In) , kde pro n N, t Z je n-rozměrný vektor tvaru = (t, t+1, . . . , t+n-1) , přičemž vždy E=0=(0, . . ., 0) a varianční matice D = (C(i, j)) n i,j=1 =2 In, kde In je jednotková matice. Pokud navíc budeme předpokládat normalitu, budeme značit t N(0, 2 ) , popř. Nn(0, 2 In) . 2.2.1 Obecné lineární regresní modely a metoda nejmenších čtverců Mějme regresní model plné hodnosti: Y = X + h(X) = h(X X) = p + 1 = k n > p + 1 Ln(0, 2 In) s vektorem závisle proměnných Y = (Y1, . . . , Yn) maticí plánu X = (xij) i = 1, . . . , n; j = 0, . . . , p vektorem 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í hodnoty (které jsou známou funkcí regresorů) a stejné rozptyly - mluvíme o homogenitě rozptylu. 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. 54 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) Pak platí (viz např. Zvára, K.: Regresní analýza. Praha. Academia. 1989): * 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(XX)-1c t(n - p - 1), kde c = (c0, c1, . . . , cp) a E(c ) = c * Označme i-tý řádek matice plánu X jako x i = (xi0, . . . , xip), pak Yi = x i + i N(x i, 2 ), Yi = x i N(x i, 2 x i(X X)-1 xi) Yi - Yi = x i( - ) + i N(0, 2 (1 + x i(X X)-1 xi)). 55 V následující tabulce uvádíme horní a dolní meze příslušných intervalů spo- lehlivosti: Intervaly spolehlivosti pro parametry j dolní mez j - t1- 2 (n-p-1)s vjj (j = 0, . . . , p) horní mez j + t1- 2 (n-p-1)s vjj pro střední hodnotu predikce dolní mez x i b - t1- 2 (n-p-1)s p x i(XX)-1xi E bYi =Ex i b =x i (i = 1, . . . , n) horní mez x i b + t1- 2 (n-p-1)s p x i(XX)-1xi pro predikci dolní mez x i b - t1- 2 (n-p-1)s p 1+x i(XX)-1xi bYi = x i b (i = 1, . . . , n) horní mez x i b + t1- 2 (n-p-1)s p 1+x i(XX)-1xi kde t1- 2 (n-p-1) je 1- 2 kvantil Studentova rozdělení o n-p-1 stupních volnosti Příklad 2.2.1. Regresní přímka v klasickém lineárním regresním modelu -1 -0.5 0 0.5 1 1.5 -4 -2 0 2 4 6 X Y Obrázek 2.1: Ukázka klasického regresního modelu s homogenním rozptylem. Klasickým speciálním případem lineárního modelu je jednoduchá lineární regrese, kdy předpokládáme, že nezávislé náhodné veličiny Yi (i = 1, . . . , n) mají normální rozdělení Yi N(i = 0 + 1xi, 2 ) , kde xi jsou dané konstanty, které nejsou všechny stejné. Rozptyly Yi jsou stejné, kdežto střední hodnoty lze vyjádřit jako lineární funkci známých konstant xi pomocí neznámých parametrů 0, 1. V tomto případě Y = Y1 ... Yn , matice plánu: X = 1 x1 ... ... 1 xn , = 1 ... n Nn(0, 2 In). 2.2.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 2.2.2. Mějme regresní model, ve kterém Y = X + , Ln(0, 2 V), V > 0, a hodnost matice h(X) = k (tj. V je pozitivně definitní), pak odhad pomocí metody nejmenších čtverců je roven = (X V-1 X)-1 X V-1 Y . 56 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 2.2.3. 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ů. 57 2.2.3 Modelování trendu Klasická dekompozice časových řad deterministickou složku trend modeluje pomocí regresní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í rozdě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 (2.1) kde f(t) je neznámá trendová funkce 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 pro 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 parametrů, 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) (2.2) přičemž 58 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; 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 (2.3) Označme reziduální rozptyl s2 = s2 k = SSE n - k = 1 n - p - 1 (Y - X) (Y - X) kde k = p + 1 (2.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)). 59 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) . 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 infor. 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. 60 61 Příklad 2.2.4. Počty zemřelých osob po dopravním úrazu v ČR v letech 1980 - 2005 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 Obrázek 2.2: Vstupní data spolu s polynomickým trendem 6. řádu. 1 2 3 4 5 6 7 8 0 0.5 1 1.5 2 2.5 3 3.5 x 10 4S k (Mean Square Error) 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 x 10 4 A k 1 2 3 4 5 6 7 8 8.5 9 9.5 10 10.5 11 AIC k 1 2 3 4 5 6 7 8 8.5 9 9.5 10 10.5 11 SRk 1 2 3 4 5 6 7 8 9.2 9.4 9.6 9.8 10 10.2 10.4 10.6 10.8 HQk (c=2) 1 2 3 4 5 6 7 8 9.8 10 10.2 10.4 10.6 10.8 11 HQk (c=3) Obrázek 2.3: Penalizační kritéria pro odhad stupně polynomu. 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 r 2 =0.4086 p=0.0004 r 2 =0.4673 p=0.0007 r 2 =0.7479 p=0.0000 r 2 =0.8792 p=0.0000 r 2 =0.8883 p=0.0000 r2 =0.9520 p=0.0000 r 2 =0.9576 p=0.0000 r 2 =0.9627 p=0.0000 coef. i degreeofpolynom Obrázek 2.4: Intervaly spolehlivosti pro parametry i (i = 1 . . . , k) pro jednotlivé stupně polynomů. (Tečkovaná čára značí polohu nuly) 62 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) = + Pp j=1(j cosjti + j sinj ti) nebo (b) f(ti) = + Pp j=1 jcos(j ti + j ) j = q 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) = Pn i=1 (Yi - f(ti))2 (b) S(, 1, . . . , p, 1, . . . , p, 1, . . . , p) = Pn 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) = 0 B @ 1 c11 s12 cp1 sp1 ... ... ... ... ... ... 1 c1n s1n cpn spn 1 C A , kde cji = cos jti sji = sin jti pro j = 1, . . . , p i = 1, . . . , n Pokud n = 2m + 1 ti = i j = 2j n pro některá j {1, . . . , m} počítejme postupně (1) Pro k = 1, 2, . . . platí nX t=1 eikjt = nX t=1 (cos kj t + i sin kj t) = eikj (1 - eikjn ) 1 - eikj = eik 2j n 1 - eik 2j n (1 - ei2kj | {z } =1 ) = 0 Protože tedy nX t=1 cos kjt + i nX t=1 sin kjt = 0, 63 pak platí nX t=1 cos kj t = nX t=1 cos k 2j n t = nX t=1 sin kjt = nX t=1 sin k 2j n t = 0 (k = 1, 2, . . .) (2.5) (2) S využitím vztahu sin cos = 1 2 sin 2 a (2.5) dostaneme nX t=1 cjtsjt = nX t=1 cos jt sin j t = 1 2 nX t=1 sin 2j t = 0 . (3) Protože cos2 = 1 2 (1 + cos 2) , pak nX t=1 c2 jt = nX t=1 cos2 jt = 1 2 nX t=1 (1 + cos 2) = n 2 . (4) Obdobně, protože sin2 = 1 2 (1 - sin 2) , pak nX t=1 s2 jt = nX t=1 sin2 jt = 1 2 nX t=1 (1 - sin 2) = n 2 . (5) Použijeme-li vztah 1 2 (cos 2 + cos 2) = cos( + ) cos( - ) , pak pro nX t=1 cjtcht = nX t=1 cos 2j n t cos 2h n t nejprve vypočteme a ze vztahů + = 2j n t - = 2h n t 2 = 2(j+h) n t (sečtením rovnic) 2 = 2(j-h) n t (odečtením rovnic) takže pro j = h platí nX t=1 cjtcht = 1 2 nX t=1 cos 2(j+h) n t | {z } =0 +1 2 nX t=1 cos 2(j-h) n t | {z } =0 = 0 . (6) Protože 1 2 (cos 2 - cos 2) = sin( + ) sin( - ) , pak pro nX t=1 sjtsht = nX t=1 sin 2j n t sin 2h n t 64 nejprve vypočteme a ze vztahů + = 2j n t - = 2h n t 2 = 2(j+h) n t (sečtením rovnic) 2 = 2(j-h) n t (odečtením rovnic) takže pro j = h platí nX t=1 sjtsht = 1 2 nX t=1 cos 2(j-h) n t | {z } =0 -1 2 nX t=1 cos 2(j+h) n t | {z } =0 = 0 . (7) Analogicky, protože 1 2 (sin 2 + sin 2) = sin( + ) cos( - ) , pak pro j = h platí nX t=1 sjtcht = 1 2 nX t=1 sin 2(j+h) n t | {z } =0 +1 2 nX t=1 sin 2(j-h) n t | {z } =0 = 0 . Nyní, jestliže využijeme předchozích vztahů, můžeme spočítat matici X X(2p+1)×(2p+1) = 0 B B B B @ n 0 0 0 n 2 ... ... ... 0 0 0 n 2 1 C C C C A a odtud velmi snadno z normálních rovnic dostaneme odhady neznámých parametrů ve tvaru ^ = 1 n nX t=1 Yt ^j = 2 n nX t=1 Yt cos jt j = 1, . . . , p. ^j = 2 n nX t=1 Yt sin jt Neznámé parametry modelu (b) získáme ze vztahů ^j = q ^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é chová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. 65 Příklad 2.2.5. Průměrné roční průtoky vody v řece Nigeru v Coulicouro (Mali) v letech 1907 až 1957 (převzato z knihy Anděl, J.: Statistická analýza časových řad, Praha SNTL 1976) 1907 1917 1927 1937 1947 1957 20 30 40 50 60 70 80 90 Obrázek 2.5: Vstupní data spolu s lineárním a trigonomickým trendem (s periodou délky 25.5 roků). Data jsou uvedena v kubických stopách za sekundu (krát 10-3). Pro známou frekvenci (získanou pomocí metody skrytých period, viz odstavec 3.3) = 2 T = 2 25.5 = 0.2464 budeme uvažovat regresní model tvaru Yt = a + bt + cos(t) + sin(t) + t , t N(0, 2 ) nebo ekvivalentní model Yt = a + bt + cos(t + ) + t , t N(0, 2 ) s maticí plánu X = 0 B @ 1 t1 cos(t1) sin(t1) ... ... ... ... 1 tn cos(tn) sin(tn) 1 C A a vektorem neznámých parametrů = 0 B B @ a b 1 C C A . Pomocí metody nejmenších čtverců obdržíme odhady ba = 54.2645 b = 9.0107 b = 2.4084 bb = 0.2101 b = -8.1678 b = -1.277 , přitom první pozorování konané v roce 1907 odpovídá t = 1. 66 MODELY LOKÁLNÍHO POSTUPNÉHO TRENDU Hlavní myšlenka lokální (vážené) metody nejmenších čtverců spočívá v tom, že provedeme 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ředpokladu, ž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 aproximová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á polynomická 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) = mX 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í b(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 cTrt 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. 67 V dalším budeme předpokládat, že posloupnost časových okamžiků t1, t2, . . . , 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 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) = mX j=0 j (x - t)j pro x [t - s, t + s]. Zaveďme substituci = x - t. Potom (x) = (t + ) = mX j=0 j j [-s, s] b(t) = c0 . Máme tedy model Yt+ = mX j=0 j j +t+ , kde [-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ášť dostáváme klasický lineární regresní model Y(t) = X + ; E = 0, D = 2 I2s+1 tj. Y(t) = 0 B B B B B B B @ Yt-s ... Yt ... Yt+s 1 C C C C C C C A = 0 B B B B B B B @ 1 -s (-s)2 . . . (-s)m ... ... ... . . . ... 1 0 02 . . . 0m .. . .. . .. . . . . .. . 1 s s2 . . . sm 1 C C C C C C C A 0 B @ 0 ... m 1 C A+ 0 B B B B B B B @ t-s ... t ... t+s 1 C C C C C C C A = 0 B B B B B B B @ x 1 ... x s+1 ... x 2s+1 1 C C C C C C C A 0 B @ 0 ... m 1 C A+ 0 B B B B B B B @ t-s ... t ... t+s 1 C C C C C C C A 68 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) bY(t) = X ` X X ´-1 X | {z } projekční matice H Y(t) = HY(t), přičemž v tomto případě X X = 0 B B B B B B B B B B B B B B B B B B B B B @ 2s + 1 0 sP =-s 2 0 sP =-s 4 sP =-s m 0 sP =-s 2 0 sP =-s 4 0 sP =-s 6 sP =-s 2 0 sP =-s 4 0 sP =-s 6 0 0 sP =-s 4 0 sP =-s 6 0 sP =-s 8 sP =-s 4 0 sP =-s 6 0 sP =-s 8 0 ... ... ... ... ... ... ... sP =-s m 0 sP =-s m+2 0 sP =-s 2m 1 C C C C C C C C C C C C C C C C C C C C C A , X Y(t) = 0 B B B B B B B B B @ sP =-s Yt+ sP =-s Yt+ ... sP =-s m Yt+ 1 C C C C C C C C C A a tudíž odhady ^Yt = x s+1 b(t) = (1, 0, . . . , 0) 0 B B B B @ ^0 ^1 ... ^m 1 C C C C A = ^0 . Označme H = 0 B B B B B B @ h 1 ... h s+1 ... h 2s+1 1 C C C C C C A . Prvky projekční matice se nazývají váhy. Pak odhad ^Yt = h s+1Y(t) . Pro prvních s hodnot, tj. pro t = 1, . . . , s vytváříme společný lineární regresní model: Y(F ) = 0 B B B B B B B B B B B @ Y1 ... Ys Ys+1 Ys+2 ... Y2s+1 1 C C C C C C C C C C C A = 0 B B B B B B B B B B B @ x 1 ... x s x s+1 x s+2 ... x 2s+1 1 C C C C C C C C C C C A 0 B B B @ 0 1 ... m 1 C C C A + 0 B B B B B B B B B B B @ 1 ... s s+1 s+2 ... 2s+1 1 C C C C C C C C C C C A bY(F ) =X (X X) -1 X Y(F ) = HY(F ) ^Y1 = h 1Y(F ) ... ^Ys = h sY(F ) 69 Obdobně pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n vytváříme společný lineární regresní model: Y(L) = 0 B B B B B B B B B B B @ Yn-2s ... Yn-s-1 Yn-s Yn-s+1 ... Yn 1 C C C C C C C C C C C A = 0 B B B B B B B B B B B @ x 1 ... x s x s+1 x s+2 ... x 2s+1 1 C C C C C C C C C C C A 0 B B B @ 0 1 ... m 1 C C C A + 0 B B B B B B B B B B B @ n-2s ... n-s-1 n-s n-s+1 ... n 1 C C C C C C C C C C C A bY(L) =X (X X) -1 X Y(L) =HY(L) ^Yn-s+1 =h s+2Y(L) ... ^Yn =h 2s+1Y(L) Pak předchozí úvahy můžeme shrnout takto: pro prvních s hodnot, tj. pro t = 1, . . . , s dostáváme ^Yt =h tY(F ) , pro tzv. "středové" hodnoty, tj. pro t = s + 1, . . . , n - s ^Yt =h s+1Y(t) a pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n ^Yt =h t-n+2s+1Y(L) . Příklad 2.2.6. Uvažujme nejprve m = 2, s = 2, r = 2s + 1 = 5 , pak matice plánu X = 0 B B B B @ x 1 x 2 x 3 x 4 x 5 1 C C C C A = 0 B B B B @ 1 -2 4 1 -1 1 1 0 0 1 1 1 1 2 4 1 C C C C A , informační matice X X = 0 @ 5 0 10 0 10 0 10 0 34 1 A a projekční matice H= 0 B B B B B B B @ h 1 h 2 h 3 h 4 h 5 1 C C C C C C C A = 0 B B B B B B B @ 31 35 9 35 -3 35 -5 35 3 35 9 35 13 35 12 35 6 35 -5 35 -3 35 12 35 17 35 12 35 -3 35 -5 35 6 35 12 35 13 35 9 35 3 35 -5 35 -3 35 9 35 3 35 1 C C C C C C C A 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) Příklad 2.2.7. Uvažujme ještě m = 3, s = 2, r = 2s + 1 = 5 , pak matice plánu X = 0 B B B B @ x 1 x 2 x 3 x 4 x 5 1 C C C C A = 0 B B B B @ 1 -2 4 -8 1 -1 1 -1 1 0 0 0 1 1 1 1 1 2 4 8 1 C C C C A , informační matice X X = 0 B B @ 5 0 10 0 0 10 0 34 10 0 34 0 0 34 0 130 1 C C A a 70 projekční matice H= 0 B B B B B B B @ h 1 h 2 h 3 h 4 h 5 1 C C C C C C C A = 0 B B B B B B B @ 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 1 C C C C C C C A 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é průměry I. Uvažujme nejprve: m = 0, r = 2s + 1 Spočtěme postupně X2s+1 = (1, . . . , 1) , X X = 2s + 1, X Y(t) = sX =-s Yt+ , normální rovnice: (2s + 1)^0 = sX =-s Yt+ ^0 = 1 2s + 1 sX =-s Yt+ . Pro t = s + 1, . . . , n - s máme ^Yt = ^(t) = ^0 = 1 2s + 1 sX =-s Yt+ . Pro prvních s hodnot, tj. pro t = 1, . . . , s dostáváme ^Yt = 1 2s + 1 2s+1X =1 Y a pro posledních s hodnot, tj. pro t = n - s + 1, . . . , n jsou odhady ve tvary ^Yt = 1 2s + 1 2s+1X =1 Yn+1- . 71 Příklad 2.2.8. Mějme m = 1, s = 2, r = 2s + 1 = 5 , pak matice plánu X = 0 B B B B @ x 1 x 2 x 3 x 4 x 5 1 C C C C A = 0 B B B B @ 1 1 1 1 1 1 C C C C A , informační matice X X = ` 5 ´ a projekční matice H = 0 B B B B @ h 1 h 2 h 3 h 4 h 5 1 C C C C A = 0 B B B B @ 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 1 C C C C A 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 II. Dále uvažujme: m = 1, r = 2s + 1 Spočtěme pomocný vzorec: sX =-s 2 = 2 sX =1 2 = 2 1 6 s(s + 1)(2s + 1) = s(s + 1)(2s + 1) 3 . Normální rovnice 0 @ 2s + 1 0 0 sP =-s 2 1 A ,, 0 1 = 0 B B @ sP =-s Yt+ sP =-s Yt+ 1 C C A ^0 = 1 2s+1 sP =-s Yt+ ^1 = 3 s(s+1)(2s+1) sP =-s Yt+ Pro t = s + 1, . . . , n - s máme ^Yt = ^(t) = ^0 = 1 2s + 1 sX =-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+1X =1 Y , ^F 1 = 3 s(s + 1)(2s + 1) sX =-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+1X =1 Yn-+1, ^L 1 = 3 s(s + 1)(2s + 1) sX =-s Yn-2s+-1 a ^Yn-s+1 = ^L 0 + ^L 1 , ^Yn-s+2 = ^L 0 + 2^L 1 , . . . , ^Yn = ^L 0 + s ^L 1 . 72 Příklad 2.2.9. Nechť m = 1, s = 2, r = 2s + 1 = 5 , pak matice plánu X= 0 B B B B @ x 1 x 2 x 3 x 4 x 5 1 C C C C A = 0 B B B B @ 1 -2 1 -1 1 0 1 1 1 2 1 C C C C A , informační matice X X= ,, 5 0 0 10 a projekční matice H= 0 B B B B B B B @ h 1 h 2 h 3 h 4 h 5 1 C C C C C C C A = 0 B B B B B B B @ 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 1 C C C C C C C A 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: * 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 W N(0, 2 ). 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 73 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 n 1 + `k 1 ´2 + `k 2 ´2 + + 1 o 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 q`2k k ´ . Položme S2 k = Pn t=k+1 U2 k,t. Pak ES2 k = (n - k)2 . Zavedeme-li pro k m + 1 veličinu Vk = Pn 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 nezaznamenáme, že tyto hodnoty začínají konvergovat k nějaké konstantě. 74 Příklad 2.2.10. Klouzavé průměry různých řádů a délek 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0, r=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0, r=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=0, r=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1, r=3 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1, r=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=1, r=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2, r=5 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2, r=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=2, r=21 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3, r=5 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3, r=11 10 20 30 40 50 60 70 30 40 50 60 70 80 90 100 m=3, r=21 Obrázek 2.6: Ilustrace vlivu řádu a délky klouzavých průměrů na vyhlazování časové řady pro data ,,Daily morning temperature of a cow . (Makridakis, Wheelwright and Hyndman (1998) Forecasting: methods and applications, John Wiley & Sons: New York. Exercises 2.3 and 2.4.) 75 Příklad 2.2.11. Klouzavé průměry různých řádů a délek 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=0, r=3 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=0, r=13 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=0, r=25 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=1, r=3 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=1, r=11 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=1, r=21 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=2, r=5 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=2, r=11 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=2, r=21 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=3, r=5 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=3, r=11 2001 2002 2003 2004 2005 2006 2007 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 x 10 6 m=3, r=21 Obrázek 2.7: Ilustrace vlivu řádu a délky klouzavých průměrů na vyhlazování časové řady pro data ,,Návštěvnost v hromadných ubytovacích zařízeních v ČR - počet hostů . Údaje o návštěvnosti v hromadných ubytovacích zařízeních vycházejí z pravidelného šetření organizovaného ČSÚ. Host v ubytovacím zařízení je každá osoba (nezahrnuje se personál a majitelé ubytovacího zařízení, kteří v ubytovacím zařízení bydlí), která použila služeb zařízení k přechodnému ubytování včetně dětí. 76 Příklad 2.2.12. Klouzavé průměry různých řádů a délek 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=0, r=3 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=0, r=7 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=0, r=11 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=1, r=3 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=1, r=7 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=1, r=11 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=2, r=5 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=2, r=11 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=2, r=15 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=3, r=5 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=3, r=11 1980 1985 1990 1995 2000 2005 800 900 1000 1100 1200 1300 1400 1500 1600 1700 m=3, r=15 Obrázek 2.8: Ilustrace vlivu řádu a délky klouzavých průměrů na vyhlazování časové řady pro data ,,Počty zemřelých osob po dopravním úrazu v ČR v letech 1980 - 2005 (viz příklad 2.2.4). 77 Příklad 2.2.13. Klouzavé průměry různých řádů a délek 10 20 30 40 50 20 30 40 50 60 70 80 90 m=0, r=3 10 20 30 40 50 20 30 40 50 60 70 80 90 m=0, r=11 10 20 30 40 50 20 30 40 50 60 70 80 90 m=0, r=21 10 20 30 40 50 20 30 40 50 60 70 80 90 m=1, r=3 10 20 30 40 50 20 30 40 50 60 70 80 90 m=1, r=11 10 20 30 40 50 20 30 40 50 60 70 80 90 m=1, r=21 10 20 30 40 50 20 30 40 50 60 70 80 90 m=2, r=5 10 20 30 40 50 20 30 40 50 60 70 80 90 m=2, r=11 10 20 30 40 50 20 30 40 50 60 70 80 90 m=2, r=21 10 20 30 40 50 20 30 40 50 60 70 80 90 m=3, r=5 10 20 30 40 50 20 30 40 50 60 70 80 90 m=3, r=11 10 20 30 40 50 20 30 40 50 60 70 80 90 m=3, r=21 Obrázek 2.9: Ilustrace vlivu řádu a délky klouzavých průměrů na vyhlazování časové řady pro data Průměrné roční průtoky vody v řece Nigeru v Coulicouro (Mali) v letech 1907 až 1957 (viz příklad 2.2.5). 78 EXPONENCIÁLNÍ VYROVNÁVÁNÍ Základy položili Holt (1957), Winters (1960) a Brown (1961). Na rozdíl od klouzavých průměrů vychází z polynomiá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, = 0, 1, . . . máme následující regresní model (pro lepší přehlednost místo j (t) opět budeme psát j ) Yt- = mX 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 = 0 B B B @ Yt-0 Yt-1 Yt-2 ... 1 C C C A X = 0 B B B @ 1 0 0 0 1 (-1)1 (-1)2 (-1)m 1 (-2)1 (-2)2 (-2)m ... ... ... ... 1 C C C A a W = 0 B B B @ 0 0 0 0 1 0 0 0 2 ... ... ... ... 1 C C C A = 0 B B B @ t t-1 t-2 ... 1 C C C A . Odhad parametrů metodou nejmenších vážených čtverců (neboť rozptyly nejsou konstantní) je dán vzorcem: b = (X WX)-1 X WY kde X WX = 0 B B B B B B B B B @ P =0 P =0 (-)1 P =0 (-)m P =0 (-)1 P =0 (-)2 P =0 (-)m+1 ... ... ... ... P =0 (-)m P =0 (-)m+1 P =0 (-)2m 1 C C C C C C C C C A a X WY = 0 B B B B B B B B B @ P =0 Yt- P =0 (-)1 Yt- ... P =0 (-)m Yt- 1 C C C C C C C C C A . 79 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 X =0 = X =0 YtPomocné vztahy 2.2.14. Protože platí (0, 1), pak P =0 = 1 1- . Protože ^Yt = b0, dostáváme ^Yt = b0 = (1 - ) P =0 Yt- = (1 - )Yt + (1 - ) P =1 Yt- subst.k=-1 = = (1 - )Yt + (1 - ) P k=0 k+1 Yt-1-k = = (1 - )Yt + (1 - ) X k=0 k Yt-1-k | {z } ^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. 80 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) | {z } 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, Massachusetts, Duxbury Press (1979) mají tvar ^Yn+ (n) - u1- 2 d (n) , ^Yn+ (n) + u1- 2 d (n) kde (n) = 1 n Pn 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. 81 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 aplikací 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 2.2.15. 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 vy- rovná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 aritmetický 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 nP t=1 (Yt - ^Yt-1)2 2. Nejprve se určí ^Y0 = 1 n2 Pn2 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. 82 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 vyrovná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) = tX j=1 ej() a D(, t) = 1 t tX 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. 83 DVOJITÉ EXPONENCIÁLNÍ VYROVNÁVÁNÍ Vyjádříme-li explicitně normální rovnice X WX = X WY pro m = 1 , dostaneme 0 B @ P =0 - P =0 - P =0 P =0 2 1 C A 0 1 ! = 0 B @ P =0 Yt- - P =0 Yt- 1 C A . Pomocné vztahy 2.2.16. Pro (0, 1) platí X =0 = 1 1 - X =0 = (1 - )2 X =0 2 = (1 + ) (1 - )3 Důkaz. Protože (0, 1), pak X =0 = 1 1 - d d X =0 ! = X =0 -1 = 1 (1 - )2 d d X =0 -1 ! = X =0 ( - 1)-2 = 2 (1 - )3 Díky předchozím vztahům dostaneme X =0 = X =0 -1 = (1 - )2 X =0 2 = X =0 ( - 1) + X =0 = 2 X =0 ( - 1)-2 + (1 - )2 = 22 (1 - )3 + (1 - )2 = (1 + ) (1 - )3 . 84 Pokračujme v řešení systému normálních rovnic 1 1- (1-)2 - (1-)2 (1+) (1-)3 ! 0 1 ! = 0 B @ P =0 Yt- - P =0 Yt- 1 C A . Pokud postupně první a druhou rovnici vynásobíme výrazem (1 - ), resp. -(1 - )2 , dostaneme: 1 - 1- -(1+) 1- ! 0 1 ! = 0 B @ (1 - ) P =0 Yt(1 - )2 P =0 Yt- 1 C A . Pomocné vztahy 2.2.17. Jestliže definujeme tzv. vyrovnávací statistiky 1. a 2. řádu St = (1 - ) X =0 Yt- a S [2] t = (1 - ) X =0 St- , pak platí St = (1 - )Yt + St-1 S [2] t = (1 - )St + S [2] t-1 a také (1 - )2 X =0 Yt- = S [2] t - (1 - )St. Důkaz. Upravujme postupně St = (1 - ) X =0 Yt- = (1 - )Yt + (1 - ) X =1 Yt- subst.k=-1 = = (1 - )Yt + (1 - ) X k=0 k Yt-1-k | {z } St-1 = (1 - )Yt + St-1 85 S [2] t = (1 - ) X =0 St- = (1 - )St + (1 - ) X =1 St- subst.:k=-1 = = (1 - )St + (1 - ) X k=0 k+1 St-1-k = (1 - )St + (1 - ) X k=0 k St-1-k | {z } S [2] t-1 = (1 - )St + S [2] t-1 . Všimněme si opět S [2] t . Protože St-k = (1 - ) X =0 Yt-k- , pak S [2] t = (1 - ) X k=0 k St-k = (1 - )2 X k=0 k X =0 Yt-k= (1 - )2 X k=0 k ^ 0 Yt-k + 1 Yt-k-1 + 2 Yt-k-2 + ~ = (1 - )2 X k=0 k Yt-k | {z } 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 X =0 YtA odtud plyne (1 - )2 X =0 Yt- = S [2] t - (1 - )St. 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 86 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 Příklad 2.2.18. Exponenciální vyrovnávání různých řádů: adaptivní a neadaptivní přístup 5 10 15 20 25 600 800 1000 1200 1400 1600 1800 2000 adaptivní procedura 5 10 15 20 25 600 800 1000 1200 1400 1600 1800 2000 =0.70 5 10 15 20 25 600 800 1000 1200 1400 1600 1800 2000 adaptivní procedura 5 10 15 20 25 600 800 1000 1200 1400 1600 1800 2000 =0.70 Obrázek 2.10: Ilustrace vlivu řádu a způsobu výběru parametru (adaptivní a neadaptivní) při exponenciálním vyrovnání časové řady pro data ,,Počty zemřelých osob po dopravním úrazu v ČR v letech 1980 - 2005 (viz příklad 2.2.4). 87 Příklad 2.2.19. Exponenciální vyrovnávání různých řádů: adaptivní a neadaptivní přístup 10 20 30 40 50 10 20 30 40 50 60 70 80 90 100 110 adaptivní procedura 10 20 30 40 50 10 20 30 40 50 60 70 80 90 100 110 =0.70 10 20 30 40 50 10 20 30 40 50 60 70 80 90 100 110 adaptivní procedura 10 20 30 40 50 10 20 30 40 50 60 70 80 90 100 110 =0.74 Obrázek 2.11: Ilustrace vlivu řádu a způsobu výběru parametru (adaptivní a neadaptivní) při exponenciálním vyrovnání časové řady pro data ,,Průměrné roční průtoky vody v řece Nigeru v Coulicouro (Mali) v letech 1907 až 1957 . 88 Příklad 2.2.20. Exponenciální vyrovnávání různých řádů: adaptivní a neadaptivní přístup 10 20 30 40 50 60 70 20 30 40 50 60 70 80 90 100 adaptivní procedura 10 20 30 40 50 60 70 20 30 40 50 60 70 80 90 100 =0.80 10 20 30 40 50 60 70 20 30 40 50 60 70 80 90 100 adaptivní procedura 10 20 30 40 50 60 70 20 30 40 50 60 70 80 90 100 =0.86 Obrázek 2.12: Ilustrace vlivu řádu a způsobu výběru parametru (adaptivní a neadaptivní) při exponenciálním vyrovnání časové řady pro data ,,Daily morning temperature of a cow . 89 2.2.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 W N(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 jk W N(0, 2 ), j = 1, . . . , r; k = 1, . . . , d. Maticově, lze tento model rozepsat takto 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ Y11 ... ... Y1d Y21 ... ... Y2d ... ... ... Yr1 ... ... Yrd 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A = 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 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 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A 0 B B B B B B B B @ m1 ... mr s1 ... sd 1 C C C C C C C C A + 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 11 ... ... 1d 21 ... ... 2d ... ... ... r1 ... ... rd 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A 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 dX k=1 sk = 0. 90 Potom X X = 0 B B B B B B B B B B B B B B B @ 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 1 C C C C C C C C C C C C C C C A a X Y = 0 B B B B B B B B @ Y1 ... Yr Y1 ... Yd 1 C C C C C C C C A , kde využíváme tzv. tečkové notace Yj = Pd i=1 Yji, Yk = Pr i=1 Yik. Normální rovnice X X = X Y můžeme přepsat (pro j = 1, . . . , r, k = 1, . . . , d) do tvaru dmj + dX i=1 si | {z } =0 = Yj mj = 1 d Yj = Yj rX i=1 mi + rsk = Yk rsk = Yk - rX i=1 mi = rX i=1 (Yik - mi) sk = 1 r rX i=1 (Yik - mi) Model: polynomický trend za celé období spolu se sezónností Uvažujme regresní model ve tvaru: Yt = Trt + Szt + t t W N(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 91 Matice plánu je pak tvaru X = 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 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 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A 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. 92 Za těchto předpokladů lze model maticově napsat ve tvaru 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ Y11 ... ... .. . Y1d Y21 ... ... .. . Y2d .. . ... ... ... Yr1 ... ... .. . Yrd 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A = 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 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 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A 0 B B B B B B B B B B B @ 0 1 . .. m s2 . .. sd 1 C C C C C C C C C C C A + 0 B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B @ 11 ... . .. ... 1d 21 ... ... .. . 2d . .. ... ... .. . r1 .. . ... ... rd 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A Odhad vektoru neznámých parametrů metodou nejmenších čtverců je dán vztahem ^ = 0 B B B B B B B B B B B @ ^0 ^1 ... ^m bs2 ... bsd 1 C C C C C C C C C C C A = (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 = nX i=1 (Yi - bYi)2 . 93 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 s3 =-16.9444 s4 =-55.6111 s5 =-180.2778 s6 =-182.9444 s7 =-45.9444 s 8 =97.7222 s 9 =183.3889 s 10 =137.7222 s 11 =-43.9444 s12 =175.0556 SSE=52446.611111 m 1 =330.5833 m 2 =567.6667 m3 =801.5833 1 2 3 4 5 6 7 8 9 10 11 12 1 -200 -150 -100 -50 0 50 100 150 200 sezonnislozkys k (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 s4 =-39.54167 s5 =-183.8333 s6 = -206.125 s7 = -88.75 s8 = 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 sezonnislozkys k (k=1,...,12) 94 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 s 2 =-1503.9028 s 3 =-723.9028 s 4 =-522.9028 s5 =338.4306 s 6 =807.5972 s 7 =1665.0972 s 8 =961.4306 s 9 =-87.4028 s10 =196.9306 s 11 =-320.5694 s 12 =-67.0694 SSE=3955160.263889 m 1 =9651.7500 m 2 =8718.5000 m 3 =8585.8333 m 4 =8396.7500 m5 =8576.8333 m 6 =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 s 3 = 57.992433 s 4 = 275.59236 s 5 = 1151.8725 s6 = 1634.333 s 7 = 2503.4736 s 8 = 1809.7946 s 9 = 769.29573 s 10 = 1060.3105 s11 = 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) 95 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 s1 =-14.5972 s2 =-33.4972 s3 =-37.4790 s4 =-37.0076 s5 =-1.5790 s6 =21.0638 s7 =25.0638 s8 =27.6210 s9 =16.1924 s10 =12.0495 s11 =4.3352 s12 =10.9638 SSE=6168.987177 m1 =145.8700 m2 =147.3833 m3 =143.6250 m4 =151.7833 m5 =137.0833 m6 =123.7500 m7 =117.6583 1 2 3 4 5 6 7 8 9 10 11 12 1 -40 -30 -20 -10 0 10 20 30 sezonnislozkysk (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. Yjk =0 +1 t+2 t2 +s2 +...+sd +jk j=1,...,7 k=1,...,12 t=(j-1)d+k s2 =0.746834 s3 = 36.4787 s4 = 59.4528 s5 = 63.812 s6 = 66.7562 s7 = 55.7426 s8 = 52.0426 s9 = 44.7991 s10 = 51.9264 s11 = 20.6894 s12 = 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) 96 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 Y jk =m j +s k + jk j=1,...,7 k=1,...,12 s 1 =-2.4461 s 2 =-3.4113 s 3 =-4.8158 s 4 =-4.7619 s 5 =0.3358 s 6 =1.1881 s 7 =2.7885 s 8 =3.4471 s 9 =1.9625 s 10 =2.4702 s 11 =1.1171 s 12 =1.2892 SSE=73.522167 m 1 =16.4329 m 2 =17.4090 m 3 =19.0226 m 4 =20.4023 m 5 =21.3367 m 6 =21.5169 m 7 =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 sezonnislozkys k (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. Yjk =0 +1 t+2 t2 +s2 +...+sd +jk j=1,...,7 k=1,...,12 t=(j-1)d+k s2 =-0.033788 s3 = 4.979 s4 = 5.749 s5 = 7.2698 s6 = 7.8515 s7 = 6.2928 s8 = 6.729 s9 = 5.307 s10 = 5.413 s11 = 2.3544 s12 = 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) 97 2.2.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 2.13: 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) = nY 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 98 Jestliže časová řada má nenulovou střední hodnotu , budeme vyšetřovat centrovanou časovou řadu Zt = Yt - . 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 2.2.21. 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 = Pn-1 t=1 Vt ES = (n - 1)1 2 DS = Pn-1 t=1 DVt + 2 Pn-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 2.2.22. 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 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 . 100 Věta 2.2.23. 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. Hypotézu H0 zamítáme, pokud počet definovaných dvojic je příliš velký nebo malý, tj. pokud |U | u1- 2 . 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) nX i=1 (i - Ri)2 . Věta 2.2.24. 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. 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 obsahují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): 0 @ aa|{z} 1 bb|{z} 2 a|{z} 3 bbb|{z} 4 aa|{z} 5 bb|{z} 6 aaaa| {z } 7 1 A 101 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 pravděpodobnost. Jako alternativní hypotézy přicházejí v úvahu 3 možnosti: 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 alternativě 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 = 8 < : 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. 102 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 2.2.25. 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ů. 103 0 20 40 60 80 100 120 140 160 180 200 -1 0 1 2 3 4 5 Yt = + 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 Yt = + 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 2.14: Demonstrace citlivosti testů náhodnosti na simulovaných datech. 104 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 2.15: Demonstrace citlivosti testů náhodnosti na reálných datech. 105 Kapitola 3 Spektrální analýza jednorozměrných časových řad 3.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. (a) 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) = R - eit dF() , kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. (b) 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) = R - eit dF() , kde F() je taková neklesající, zprava spojitá funkce, že F(-) = 0 a F() = (0). Přitom F() je jediná. Vzorci (t) = Z - eit dF() resp. (t) = Z - eit dF() se říká spektrální rozklad autokovarianční funkce. Funkce F() se nazývá spektrální distribuční funkce. 106 Připomeňme: 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() = Z f(x)dx resp. F() = Z f(x)dx. (3.1) K existenci spektrální hustoty stacionární náhodné posloupnosti, resp. spojitého stacionární náhodného procesu, stačí, aby pro její kovarianční funkci platilo X t=|(t)| < resp. Z |(t)|dt < . Existuje-li spektrální hustota f() stacionární posloupnosti a má-li variaci konečnou na -, , pak platí f() = 1 2 X t=- e-it (t) (3.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. Z |(t)|dt < , pak f() = 1 2 Z - e-it (t) dt. (3.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(-) (3.4) skoro všude vzhledem k Lebesqueově míře. 3.2 Periodogram V dalším budeme předpokládat, že {Yt, t Z} je centrovaná stacionární náhodná posloupnost. Definice 3.2.1. Nechť Y1, . . . , Yn jsou pozorování náhodné posloupnosti {Yt, t Z}. Pak periodogram definujeme vztahem In() = 1 2n ˛ ˛ ˛ ˛ nP t=1 Yte-it ˛ ˛ ˛ ˛ 2 [-, ]. 107 Lemma 3.2.2. Položme An() = r 2 n nX t=1 Yt cos t Bn() = r 2 n nX t=1 Yt sin t, pak platí In() = 1 4 ^ A2 n() + B2 n() ~ . Důkaz. In() = 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it ˛ ˛ ˛ ˛ ˛ 2 = 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yt cos t - i nX t=1 Yt sin t ˛ ˛ ˛ ˛ ˛ 2 = = 1 2n " nX t=1 Yt cos t !2 + nX t=1 Yt sin t !2# = 1 4 ^ A2 n() + B2 n() ~ . Poznámka 3.2.3. Někteří autoři definují periodogram poněkud jinak: I n() = 2 n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it ˛ ˛ ˛ ˛ ˛ 2 = ^ A2 n() + B2 n() ~ = 4In(). Lemma 3.2.4. Pokud označíme pro k = 0, 1, . . . , n - 1 Ck = 1 n - k n-kX t=1 YtYt+k C k = 1 n n-kX t=1 YtYt+k pak platí In() = 1 2 " C0 + 2 n-1X k=1 ` 1 - k n ´ Ck cos k # = 1 2 " C 0 + 2 n-1X k=1 C k cos k # . 108 Důkaz. In() = 1 2n " nX t=1 Yt cos t !2 + nX t=1 Yt sin t !2# = 1 2n " nX t=1 Yt cos t ! nX s=1 Ys cos s ! + nX t=1 Yt sin t ! nX s=1 Ys sin s !# = 1 2n nX t=1 nX s=1 YtYs (cos t cos s + sin t sin s) = 1 2n nX t=1 nX 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 |{z} ) t min(n, z }| { n-k). týká se záporných k a pak platí In() = 1 2n n-1X k=-n+1 cos k min(n,n-k) X 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 nX t=1 Y 2 t | {z } C0 + 1 2 -1X k=-n+1 n - |k| n cos k 1 n - |k| nX t=1-k YtYt+k | {z } C-k=Ck + 1 2 n-1X k=1 n - k n cos k 1 n - k n-kX t=1 YtYt+k | {z } Ck = = 1 2 n-1X k=-(n-1) ,, 1 - |k| n Ck cos k = 1 2 " C0 + 2 n-1X k=1 ,, 1 - k n Ck cos k # In() = 1 2 1 n nX t=1 Y 2 t | {z } C 0 + 1 2 -1X k=-n+1 cos k 1 n nX t=1-k YtYt+k | {z } C -k =C k + 1 2 n-1X k=1 cos k 1 n n-kX t=1 YtYt+k | {z } C k = 1 2 C 0 + 2 n-1X k=1 C k cos k ! . 109 Poznámka 3.2.5. K numerickému výpočtu hodnot periodogramu se často používají právě předchozí vzorce. Poznámka 3.2.6. Pro teoretické účely bývá výhodnější tato varianta In() = 1 2 n-1X k=-(n-1) " 1 - |k| n " Ck cos k = 1 2 n-1X k=-(n-1) C k cos k. Pro náhodnou posloupnost {Yt, t T Z} platí f() = 1 2 X 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 3.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) nX t=1 eit = ei ein - 1 ei - 1 = ei e i n 2 1 2i " e i 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 nX t=1 eit = n pro = 2k (2) ˛ ˛ ˛ ˛ ˛ nX t=1 eit ˛ ˛ ˛ ˛ ˛ 2 = nX t=1 eit ! nX s=1 e-is ! = nX t=1 nX 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 ˛ ˛ ˛ ˛ ˛ nX t=1 eit ˛ ˛ ˛ ˛ ˛ 2 = n2 pro = 2k 110 (3) Nechť funkce g() je integrovatelná na -, , pak pro každý bod -, , v němž je g() spojitá, platí: lim n 1 2n Z 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 ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it ˛ ˛ ˛ ˛ ˛ 2 = E 1 2n nX t=1 nX s=1 YtYse-i(t-s) = 1 2n nX t=1 nX s=1 e-i(t-s) E(YtYs) | {z } (t,s)=(t-s) = 1 2n nX t=1 nX s=1 e-i(t-s) Z - ei(t-s) f() d | {z } =(t-s) = 1 2n Z - nX t=1 nX s=1 e-i(t-s) ei(t-s) f() d = 1 2n Z - nX t=1 nX s=1 eit(-) e-is(-) | {z } =| Pn t=1 eit(-) |2 f() d = 1 2n Z 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 nX t=1 nX s=1 e-i(t-s) E(YtYs) = 1 2n nX t=1 EY 2 t = 1 2n n2 = 2 2 . Dále pomocí autokovarianční funkce počítejme spektrální hustotu f() = 1 2 X t=- e-it (t) = ( 2 2 -, , 0 jinak 2 2 /2 - 111 Než budeme pokračovat dále, poznamenejme, že pro Yt IID(0, 2 ) platí E (YtYsYuYv) = 8 >>>>>>< >>>>>>: 4 t = s = u = v, 4 t = s = u = v, 4 t = u = s = v, 4 t = v = s = u, 0 jinak. Pak pro 1, 2 -, počítejme E [In(1)In(2)] = E " 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it1 ˛ ˛ ˛ ˛ ˛ 2 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it2 ˛ ˛ ˛ ˛ ˛ 2# = 1 42n2 E ( nX t=1 nX s=1 e-it1 eis1 YtYs nX u=1 nX v=1 e-iu2 eiv2 YuYv ) = 1 42n2 nX t=1 nX s=1 nX u=1 nX v=1 ei(s-t)1 ei(v-u)2 E (YtYsYuYv) = 1 42n2 8 < : nX t=1 EY 4 t + nX t=1 nX u=1,t=u EY 2 t Y 2 u + nX t=1 nX s=1,s=t ei(s-t)1 ei(s-t)2 EY 2 t Y 2 s + + nX t=1 nX u=1,u=t ei(u-t)1 ei(t-u)2 EY 2 t Y 2 u 9 = ; = 1 42n2 ( n4 + n(n-1)4 + 4 nX t=1 nX s=1 ei(s-t)(1+2) - n4 = +4 nX t=1 nX s=1 ei(s-t)(1-2) - n4 ) = 1 42 ( 4 + 4 - 34 n + 4 n2 "˛ ˛ ˛ ˛ ˛ nX t=1 eit(1+2) ˛ ˛ ˛ ˛ ˛ 2 + ˛ ˛ ˛ ˛ ˛ nX 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 ­ff = 1 42 24 + 4 - 34 n + 4 n2 sin2 n sin2 ff . Pro 1 = 2 = = 2k (k = 0, 1) platí EI2 n() = 1 42 4 + 4 - 34 n + 4 n2 ^ n2 + n2~ ff = 1 42 34 + 4 - 34 n ff 112 Odtud díky vztahu DIn() = EI2 n() - (EIn())2 , kde EIn() = 2 2 dostaneme DIn() = 8 < : 1 42 n 4 + 4-34 n + 4 n2 sin2 n sin2 o ---- n 4 42 = f2() pro = 0, (-, ), 1 42 n 24 + 4-34 n o ---- n 2 4 42 = 2f2() pro = 0, . Poznámka 3.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. 113 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 3.1: Simulovaná časová řada. 114 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 3.2: Ukázka nekonzistentnosti periodogramu, neboť jeho variabilita se nesnižuje se vzrůstajícím počtem pozorování. 115 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 3.3: Na obrázku je demonstrován tzv. Welchův přístup, který vycházel ze zmenšová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) 116 3.3 Metoda skrytých period. 3.3.1 Regresní model s periodickým trendem Uvažujme nyní takové časové řady, které můžeme rozložit na součet harmonických frekvencí, 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 W N(0, 2 ), tj. Et = 0; Dt = 2 ; C(t, s) = 0; t = s, st je periodická funkce tvaru st = ( Pp j=1(j cos tj + j sin tj) pro n liché Pp 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 = ( Pp j=1 j cos(tj - j ) pro n liché Pp j=1 j cos(tj - j ) + n 2 (-1)t pro n sudé kde k = p 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é harm. složky 117 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. 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 X k=1 (k cos tk + k sin tk) = 0 2 + p X k=1 h 1 2 " eitk + e-itk " k + 1 2i " eitk - e-itk " k i = 0 2|{z} c0 + p X k=1 2 6 4 1 2 (k - ik) | {z } ck eitk + 1 2 (k + ik) | {z } ck e-itk 3 7 5 = p X 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 X t=- e-it (t) = ( 2 2 -, , 0 jinak 2 2 /2 - Všimněme si periodogramu této náhodné posloupnosti In() = 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 p X j=-p cjeitj + t ! e-it ˛ ˛ ˛ ˛ ˛ 2 = 1 2 ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ 1 n nX t=1 p X j=-p cj eit(j -) | {z } I (1) n () + 1 n nX t=1 te-it | {z } I (2) n () ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ ˛ 2 118 Upravujme dále I(1) n () = 1 n p X j=-p cj nX t=1 eit(j-) | {z } součet geom.řady = 1 n p X j=-p cj ei(j-) ein(j -) - 1 ei(j -) - 1 | {z } 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 X j=-p,j=k gn(j - ) = nck | {z } + 1 n p X j=-p,j=k gn(j - ) | {z } 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 nX t=1 te-it ˛ ˛ ˛ ˛ ˛ 2 = 1 n E nX t=1 te-it ! nX s=1 teis ! = 1 n nX t=1 nX s=1 e-i(t-s) Ets | {z } nekorel. = 1 n nX t=1 E2 t = 1 n n2 = 2 Pokud bude platit, že náhodná posloupnost splňuje model Yt = st + t = 0 2 + p X j=1 (j cos tj + j sin tj ) + t t W N(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é frekvence 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. 3.3.2 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 W N(0, 2 ) 119 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 = Pp j=1 j cos(tj - j) + t t W N(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) = X k=1,2,... m k ! [max(0, 1 - kx)]m-1 = X k=1,2,... m k ! [(1 - kx)+]m-1 . Pak hypotézu H0 zamítáme (na hladině významnosti ), pokud 1 - FW |H0 (w) = P(W > w|H0) = W |H0 , 120 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ů. 3.3.3 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 = mX j=1 0 B B @ In(j) mP i=1 In(i) - gF 1 C C A + , 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. 121 20 40 60 80 100 120 -10 0 10 20 30 40 50 60 70 80 y t = 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|zk |; zk =ak +i.bk stredova symetrie 0 20 40 60 80 100 120 -200 -150 -100 -50 0 50 100 150 200 FAZOVE SPEKTRUMk =arctg(bk /ak ) stredova antisymetrie stupen 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 x 10 4 PERIODOGRAM I k =(2/n)|z k | 2 ; z k =a k +i.b k Fisher's test 40.0 Obrázek 3.4: Příklad periodické časové řady. 122 20 40 60 80 100 120 0 20 40 60 80 100 120 y t = 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|zk |; zk =ak +i.bk stredova symetrie 0 20 40 60 80 100 120 -200 -150 -100 -50 0 50 100 150 200 FAZOVE SPEKTRUMk =arctg(bk /ak ) stredova antisymetrie stupen 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 x 10 4 PERIODOGRAM I k =(2/n)|z k | 2 ; z k =a k +i.b k Fisher's test 40.0 120.0 60.0 Obrázek 3.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 nestacionárních časových řad trend z řady nejprve odstranit. 123 3.3.4 Odhady neznámých parametrů modelu a periodogram Jestliže jsme pomocí předchozích testů nalezli významné frekvence 1, . . . , p, můžeme uvažovat regresní model Yt = st + t , t W N(0, 2 ) (např. gaussovský bílý šum) kde st = + p X j=1 (j cos tj + j sin tj) = + p X j=1 j cos(tj + j ) a j = q 2 j + 2 j j = arctan j j . Uvažujeme-li regresní model v maticovém zápisu Y = X + , ve kterém Y = 0 B @ Y1 ... Yn 1 C A , = 0 B @ 1 ... n 1 C A , = 0 B B B B B B B @ 1 1 ... p p 1 C C C C C C C A . s maticí plánu X = 0 B @ 1 c11 s12 cp1 sp1 ... ... ... ... ... ... 1 c1n s1n cpn spn 1 C A , kde cji = cos jti sji = sin jti pro j = 1, . . . , p i = 1, . . . , n , pak, jak již bylo ukázáno na straně 62 X X = 0 B B B B @ n 0 0 0 n 2 ... ... ... 0 0 0 n 2 1 C C C C A . Z normálních rovnic X X = X Y dostaneme odhady neznámých parametrů ve tvaru ^ = 1 n nX t=1 Yt ^j = 2 n nX t=1 Yt cos jt j = 1, . . . , p. ^j = 2 n nX t=1 Yt sin jt 124 Připomeňme, že In() = 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it ˛ ˛ ˛ ˛ ˛ 2 = 1 4 ^ A2 n() + B2 n() ~ kde An() = r 2 n nX t=1 Yt cos t Bn() = r 2 n nX t=1 Yt sin t. Z předchozího dostáváme ^j = q 2 n An(j) ^j = q 2 n Bn(j) ^j = q ^2 j + ^2 j = q 2 n A2 n(j) + 2 n B2 n(j) = q 2 n p A2 n(j) + B2 n(j) = q 2 n p 4In(j) = 2 q 2 n In(j) ^j = arctan Bn(j) An(j) , takže, pokud spočítáme periodogram, jeho hodnoty pro významné frekvence j můžeme použít i pro odhady neznámých parametrů. 125 Příklad 3.3.1. Průměrné roční průtoky vody v řece Nigeru v Coulicouro (Mali) v letech 1907 až 1957 (převzato z knihy Anděl, J.: Statistická analýza časových řad, Praha SNTL 1976) 1907 1917 1927 1937 1947 1957 20 30 40 50 60 70 80 90 Obrázek 3.6: Vstupní data spolu s lineárním a trigonomickým trendem (s periodou délky 25.5 roků). Data jsou uvedena v kubických stopách za sekundu (krát 10-3). 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 250 300 350 25.50 Významná hodnota periodogramu byla pomocí Fisherova testu určena pro frekvenci bF = 2 bTF = 2 25.5 = 0.2464 Pomocí iterativní Damsleth-Spotvoll metody získáme bD = 2 bTD = 2 27.9 = 0.2252 přitom uvažujeme regresní model tvaru Yt = a + bt + cos(t) + sin(t) + t , kde t N(0, 2 ). Pomocí metody nejmenších čtverců obdržíme odhady baF = 54.2645 bF = 9.0107 bbF = 0.2101 bF = -8.1678 a baD = 54.2645 bD = 3.7386 bbD = 0.2101 bD = -11.9439 . 126 3.4 Odhady spektrální hustoty 3.4.1 Neparametrické odhady spektrální hustoty (Window Spectral Estimation) Neparametrické odhady spektrální hustoty centrované stacionární náhodné po- sloupnosti {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 3.2.4) In() = 1 2n ˛ ˛ ˛ ˛ ˛ nX t=1 Yte-it ˛ ˛ ˛ ˛ ˛ 2 = 1 2 " C 0 + 2 n-1X k=1 C k cos k # . Využijme dále vztahů C k = C -k, kde C k = 1 n n-kX t=1 YtYt+k pro k = 0, 1, 2, . . . , (n - 1) a cos k = 1 2 " eik + e-ik " . Upravujme postupně In() = 1 2 " C 0 + n-1X k=1 C k eik + n-1X k=1 C k e-ik # = 1 2 2 4C 0 + -1X s=-(n-1) C -se-is + n-1X k=1 C k e-ik 3 5 = 1 2 n-1X k=-(n-1) C k e-ik . 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. 127 C 0 = 1 n ` Y 2 1 + + Y 2 n ´ | {z } n členů C 1 = C -1 = 1 n (Y1Y2 + + Yn-1Yn + Y3Yn) | {z } n-1 členů ... C n-3 = C -(n-3) = 1 n (Y1Yn-2 + Y2Yn-1 + Y3Yn) | {z } 3 členy C n-2 = C -(n-2) = 1 n (Y1Yn-1 + Y2Yn) | {z } 2 členy C n-1 = C -(n-1) = 1 n Y1Yn | {z } 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 mX k=-m C k cos k = 1 2 mX k=-m C k e-ik , což lze také zapsat takto ^fn() = 1 2 n-1X k=-(n-1) w(k)C k cos k = 1 2 n-1X k=-(n-1) w(k)C k e-ik , kde w(k) = 1 |k| m 0 |k| > m . Označme Fourierovu transformaci funkce w(k) W () = 1 2 X k=- w(k)e-ik = 1 2 mX k=-m e-ik a řadu přeindexujeme tak, aby indexy šly od 1 do 2m + 1, tj. položme s = k + m + 1, pak k = s - m - 1 a 128 (i) pro = 2k W () = 1 2 2m+1X s=1 e-i(s-m-1) = 1 2 ei(m+1) 2m+1X 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 a (ii) pro = 2k W () = 2m + 1. Vzhledem k tomu, že lze psát In() = 1 2 n-1X k=-(n-1) C k e-ik , vidíme, že In() je Fourierovou transformací C k , takže naopak lze pomocí inverzní Fourierovy transformace psát C k = Z - In()eik d . Počítejme postupně ^fn() = 1 2 n-1X k=-(n-1) w(k)C k e-ik = 1 2 n-1X k=-(n-1) w(k) Z - In()eik d e-ik = Z - In() 1 2 n-1X k=-(n-1) w(k)e-ik(-) | {z } W (-) d = Z In()W ( - )d . 129 Jde o tzv. vyhlazený periodogram (smoothed periodogram). Funkce W () se nazývá spektrální okénko (spectral window). Tato funkce má do jisté míry aproximovat Diracovu funkci a platí pro ni Z 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-1X k=-(n-1) w(k)C k e-ik , přičemž inverzní Fourierova transformace w(k) = Z 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řirozené čí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 ). 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 Lag window w(k) -2 0 2 -0.2 0 0.2 0.4 0.6 0.8 1 Spectral window W()-Dirichlet kernel -6/7 -4/7 -2/7 0 2/7 4/7 6/7 Obrázek 3.7: Korelační a spektrální okénko v případě prostého useknutého perio- dogramu. 130 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. -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 Lag window w(k) -4 -2 0 2 4 0 0.2 0.4 0.6 0.8 Spectral window W()-Fejer kernel Obrázek 3.8: Bartletovo korelační a spektrální okénko. Parzenovo okénko w(k) = 8 >>>< >>>: 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. -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 Lag window w(k) -4 -2 0 2 4 0 2 4 6 8 10 12 Spectral window W() Obrázek 3.9: Parzenovo korelační a spektrální okénko. 131 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 . Pokud a = 1 4 , pak se nazývá Tukey-Hanningovo okénko. -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 Lag window w(k) -4 -2 0 2 4 0 0.2 0.4 0.6 0.8 Spectral window W() Obrázek 3.10: Tukey-Hanningovo korelační a spektrální okénko. 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 ´ -10 -5 0 5 10 0 0.2 0.4 0.6 0.8 1 Lag window w(k) -4 -2 0 2 4 0 0.2 0.4 0.6 0.8 Spectral window W() Obrázek 3.11: Tukey-Hammingovo korelační a spektrální okénko. 132 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 (-, ). Pro k = 1, 2, . . . (n - 1) počítejme nejprve odpovídající korelační okénko: w(k) = Z W ()eik d = Z - 1 2 eik d = 1 2 eik ik ­ - = 1 k 1 2i " eik - e-ik " | {z } sin k = sin k k . Pro k = 0 je zřejmě rovno jedné, celkově tedy wk = ( 1 k = 0 sin k k k = 1, 2, . . . . -10 -5 0 5 10 -0.2 0 0.2 0.4 0.6 0.8 1 Lag window w(k) -4 -2 0 2 4 0 0.05 0.1 0.15 0.2 0.25 Spectral window W() Obrázek 3.12: Daniellovo korelační a spektrální okénko. 133 Literatura [1] ANDĚL, J. Matematická statistika. SNTL/ALFA Praha, 1978. [2] ANDĚL, J. Statistické metody. Matfyzpress Praha, 1993. [3] ANDĚL, J. Statistická analýza časových řad. Praha. SNTL 1976. [4] BEKLOVÁ, M., NĚMCOVÁ, M., PIKULA, J. Longsterm trends in fluctuation of the population level of the chosen game species in the ČSSR. Proceedings of a XVI. International Congres of game biologists 25.9.-2.10.1983, High Tatras, ČSSR, 1983 [5] BOWERMAN, B.L., OCONNELL, R.T. Time series and forecasting. North Scituate, Massachusetts, Duxbury Press. 1979. [6] BOX, G., JENKINS, G. Time series analysis - forecasting and control. Holden-Day 1976. [7] BROCKWELL, P.J., DAVIS, R.A. Time Series: Theory and Methods. SpringerVerlag, New York, 2-nd edition, 1991. [8] BROWN, R.G. Statistical forecasting for inventory control, New York. McGrawHill. 1959. [9] CIPRA, T. Analýza časových řad s aplikacemi v ekonomii. SNTL, Praha, 1986. [10] DAMSLETH, E., SPJTVOLL, E. Estimation of Trigonometric Components in Time Series, J.Amer. Statistics, Assoc. 77. 1982, pp. 382-387. [11] DANIELS, H.E. Rank correlation and population models, J. R. Statist. Soc., B, 12 1950. 171-181. [12] DOOB, J.L. Stochastic processes. New York, Wiley 1953. [13] FORBELSKÁ, M. Detekce periodicity v hydrologických datech. In XIII. letní škola bometriky, Biometrické metody a modely v současné vědě a výzkumu. 1. vyd. Brno: ÚKZÚZ Brno, 1998. s. 173-178. [14] GICHMAN, I.I., SKOROCHOD, A.V. Teorija slučajnych processov. Moskva. Nauka 1971. [15] HOLT, C.C. Forecasting seasonal and trends by exponentially weighted moving averages. Office of Naval Research, Research Memorandum No. 52. 1957. [16] HAMILTON, J.D. Time Series Analysis. Princeton University Press. 1994. [17] KALMAN, R. A new approach to linear filtering and prediction problems. Trans. ASME J. Basic Eng. D 82 (1960), 34-45. [18] Mann, H.B.: Non-parametric tests against trend, Econometrica, 13, 1945. 245-259. 134 [19] 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. [20] 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. [21] NEUBRUNN, T., RIEČAN, B. Miera a integrál. Bratislava. Veda 1981. [22] PRIESTLEY, M. Spectral analysis and time series. Academic Press 1989. [23] RAO, R.C. Lineární metody statistické indukce a jejich aplikace. ACADEMIA Praha, 1978. [24] STUART, A. The Power of Two Difference-Sign Tests, Vol. 47, Issue 259, Sep. 1952, 416-424. [25] WINTERS, P.R. (1960) Forecasting sales by exponentially weighted moving averages Management Science. 6. 324342. [26] ZVÁRA, K. Regresní analýza Praha. Academia. 1989.