ČASOVÉ ŘADY VÍTĚZSLAV VESELÝ PODPORA VÝUKY NA POČÍTAČOVÉ SÍTI http://www.math.muni.cz/~veselý Algoritmy pro časové řady v MATLABu: d = vesely('cas .rady') ... připojení a výpis obsahu knihovny d ... výstup=úplná cesta do adresáře knihovny. 1. ÚVOD 1.1. Ukázky časových řad. Časová řada = soubor pozorování nějaké veličiny {xt, ťST}, kde t je zpravidla čas aTCl. Obr. 1.1: Měření proudu procházejícího odporem r v obvodu se střídavým napětím, v(t) — acos(vt + 0), tj. xt — fcos(vt + 0). Obr. 1.2 (MATLAB: getdata( 'brockwell.dat/uspop.dat')): Růst populace v USA v létech 1790-1980 sledovaný v desetiletých intervalech. Obr. 1.3 (MATLAB: getdata( 'brockwell .dat/strikes .dat')): Počet stávek v USA v létech 1951-1980. Obr. 1.4: Výsledky Národní a Americké ligy v baseballu v létech 1933-1980. Obr. 1.5 (MATLAB: getdata( 'brockwell.dat/sunspots.dat')): Počty slunečních skvrn v létech 1770-1869. Obr. 1.6 (MATLAB: getdata( 'brockwell .dat/deaths . dat')): Počet úmrtí při nehodách v USA v létech 1973-1978. 1.2. Oblasti uplatnění metod analýzy časových řad. fyzika, technika: • seismický záznam v geofyzice. • řada nej vyšších denních teplot v meteorologii. • průběh výstupního signálu určitého elektrického přístroje. • tenzometrické měření povrchového napětí v provozu namáhané strojní součástky, biologie, ekologie: sledování různých parametrů znečištění ovzduší. medicína: záznam EKG nebo EEG. společenské vědy: změny v počtu a složení obyvatelstva, sociologie: vývoj rozvodovosti. EKONOMIE: Teorie časových řad = jedna z nej důležitějších kvantitativních metod pro analýzu ekonomických dat, např.: • analýza poptávky po určitém výrobku • analýza objemu zemědělské produkce • analýza počtu cestujících v letecké dopravě • analýza vývoje kurzu akcií na burze Date: 27. května 1999. 1 2 vítězslav veselý 1.3. Cíl analýzy. Porozumění mechanismu, jímž se generují sledované údaje. Znalost modelu tohoto mechanismu =>• znalost algoritmu, jímž můžeme chování tohoto mechanismu simulovat na počítači =>• schopnost popsat s jistou přesností jeho chování: • mezi časovými okamžiky měření (interpolace) • v budoucnosti (extrapolace, prognóza) • s cílem řídit a optimalizovat činnost určitého systému vhodnou volbou vstupních a počátečních podmínek (regulace), např. regulace složitých technologických procesů. 1.4. Některé specifické problémy analýzy časových řad. a) Volba okamžiků pozorování: • jsou přímo diskrétní svou povahou, např. úroda obilí za jednotlivé roky. • vznikají diskretizací spojité časové řady, např. teplota v danou denní dobu na daném místě. • vznikají akumulací (agregací) hodnot za určité období, např. denní množství srážek, roční výroba závodu. Místo akumulace se někdy provádí průměrování. Je-li dána možnost volby, je třeba jí věnovat pozornost: • málo bodů =>• unikne charakteristický rys řady. • mnoho bodů =>• zvýší se výpočetní náročnost. • ekvidistantní diskretizace zpravidla usnadní numerické zpracování, ale neumožňuje adaptivně měnit hustotu diskretizace v závislosti na lokálním charakteru řady. • při agregaci se mohou porušit vlastnosti původní řady. b) Problémy s kalendářem: • různá délka kalendářních měsíců. • 4 nebo 5 víkendů v měsíci. • různý počet pracovních dnů v měsíci. • pohyblivé svátky: např. svátek na začátku měsíce sníží prodej potravin za tento měsíc, ale zvýší jej za předchozí v důsledku efektu predzásobení. Příklad 'očištění' časové řady např. od proměnlivé délky měsíce: zavedeme 'standardní' měsíc o délce 30 dnů a pak údaj třebas o produkci za leden přenáso-bíme korekčním faktorem c) Problémy s nekompatibilitou jednotlivých měření: Příklad: hodnota nějakého ukazatele se jeden rok týká např. 85 podniků, další rok jen 82 apod. d) Problémy s délkou časových řad: • zvětšení počtu měření (např. půlením časových intervalů mezi body) nemusí vždy znamenat zvětšení množství informace. • někdy se mohou objevit protichůdné tendence: metoda zpracování vyžaduje delší řadu, ale na druhé straně řada vzniklá dlouhodobým sledováním může měnit charakteristiky svého modelu v čase =>• obtíže s konstrukcí modelu. 1.5. Označení a základy matematické statistiky. s :— v, resp. v —: s . .. označení výrazu v symbolem s. Číselné obory: N .. . množina všech přirozených čísel N0 — N U {0} .. . množina všech nezáporných celých čísel "E .. . množina všech celých čísel M .. . množina všech reálných čísel M+ .. . množina všech nezáporných reálných čísel C .. . množina všech komplexních čísel (•)+ : M —¥ M+ . .. zobrazení definované předpisem (x)+ — max(0, x) (a,b) ... otevřený interval časové rady 3 [a, b] ... uzavřený interval J(a, b) — {x | min(a, b) < x < max(a, 6)} J[a, 6] — | min(a, 6) < x < max(a, 6)}. Vektory: a; — (a?i,. .., Kn)T ěC ... zpravidla sloupcové x + h — [x\ + h,. .. ,xn + h)T, IiěC * = (*i, • • -,tk)T e Mk, ti G {1,. ..,«} pro i - 1,.. ., k, k < n =>• xt := («tl, . • ., «tfc)T G C* 1 < i < n =>• as (i) - (aľi,. .., a?t-_i, a?f-+i, .. ., xn)T f(x) =/(a;i,.. .,£„), cřa; - daľi .. . cfen 0, 0nxi ... nulový sloupcový vektor délky n. Matice: A, Amxn — [ciij] — [A(i, j)] . .. matice rozměru m x n 3l(A) — {y\y — Ax] .. . obor hodnot operátoru A IN (A) — {x | Ax — 0} ... jádro operátoru A AT — [a j i] . .. transponovaná matice A* — [ä j i] ... hermitovsky sdružená matice /, /„ — Inxn — [Síj] .. . jednotková matice řádu n \A\ ... determinant čtvercové matice A 0, 0mXn — [0] • • • nulová matice rozměru m x n a?i 0 ... 0 0 x2 ... 0 diag(a?) diagonální matice 0 0 ... x„ A(i, :) — (au, .. ., cLin) —: r,- ... i-tý řádek matice A ve stylu MATLAB A(:, j) — (aij, .. ., amj)T —: sj ... j-tý sloupec matice A ve stylu MATLAB A — [ri; .. .; rm] — [si,. .., s„] ... blokový zápis matice A ve stylu MATLAB A > 0 (resp. A > 0) ... pozitivně (semi)definitní matice. Norma a skalární součin vektorů: (x, y) - Ya=i xi Vi - V*x> speciálně: y - A*xA$_yi_= (x,A(:,i)) pro i = 1,2,. . .,n ||a;|| — \J(x,x) — \/£ľ=i \xi\2 ■ ■ ■ Euklidovská norma. Schwarzova nerovnost: |(a;,y)| < ||as||||y||, kde rovnost nastane právě když vektory x a y jsou lineárně závislé. Pravděpodobnostní prostor (íí, A, P): Í2 .. . základní prostor elementárních jevů A C 2n ... tr-algebra náhodných jevů P .. . pravděpodobnostní míra na A Všechny náhodné veličiny budeme vždy uvažovat nad týmž pravděpodobnostním prostorem. Komplexní náhodnou veličinou budeme rozumět veličinu X — Xr + iXi, kde Xr a Xi jsou reálné náhodné veličiny představující po řadě reálnou a imaginární část X X — (Xi,. .., X„)T .. . (komplexní) náhodný vektor tvořený (komplexními) náhod, veličinami X,. Střední hodnota: (i — fix — EX . .. střední hodnota náhodné veličiny X fj, — — EX — (EX\, • • •, EXn)T . .. střední hodnota náhodného vektoru X. Rozptyl a kovariance náhodných veličin: er2 = a\ = varX = E|X - EX|2 = E|X|2 - |EX|2 > 0 . .. rozptyl X (TXY = cov(X, Y) = F,(X - F,X)(Y - EY) = VXY - (EX)(ĚY) . .. kovariance laľ cov(X, X) = varX, cov(y, X) = cov(X,Y) cov(Er- Xr, ^2S Ya) — ^2r ^2S cov(Xr, Ys) a odtud speciálně: var(X + Y) = varX + cov(X, Y) + cov(V, X) + varY = varX + 2Re cov(X, Y) + varY. Varianční a kovarianční matice náhodných vektorů: Ex = varX = [covpf;, Xj)] = E(X - EX)(X- EX)* = EXX* - (EX)(EX)* . .. varianční matice X Exy = cov(X, Y) = [coy(Xi,Yj)] = E(X - EX)(Y- EY)* = EXY* - (EX)(EY)* ... .. . kovarianční matice X a Y 4 vítězslav veselý cov(X,X) = varX, cov(Y,X) = cov(X, Y)* =>• varX = (varX)* .. . ... varianční matice X je hermitovská. Pro konstantní komplexní vektory a a c a matice B a D odpovídajících rozměrů platí: cov(a + BX,c+ DY) = cov(BX, DY) = B cov(X, Y) D*. JJ. X = Y var(a + BX) = cov(a + BX, a + BX) = cov(BX, BX) = B var(X) B*. \yb* = B 0 < var(6*X) — 6*varX6=> varX > 0 .. . varianční matice je pozitivně semidefinitní. Ex je tedy celkem hermitovská a pozitivně semidefinitní a má proto reálná nezáporná vlastní čísla A,-. 1 1 11 Zřejmě existuje matice E|-, jejíž vlastní čísla jsou A? a přitom platí: Ex — E|- E^. Dále platí cov(£r- Es ^'») = Er- Es cov(^r-, Ys) a odtud speciálně: var(X + Y) = varX + cov(X, Y) + cov(Y, X) + varY = varX + 2Re cov(X, Y) + varY. 1.6. Normální rozdělení a rozdělení z něj odvozená. 1. Normální rozdělení: X ~ N(fi, cr2), fi = EX, cr2 = varX . . . . .. reálná náhodná veličina s normálním (gaussovským) rozdělením; ,— (»-M)2 f(x) = (v2ir• a + bX ~ N(a + bfi,b2(t2) pro a,b G M; X~Nn(fj,,V) =>• a + 5X~ Nm(a + Bfj,,BVBT) pro a G Mm a matici 5 = Bmxn nad M. 2. Rozdělení x2(n): Nechť č/,- ~ JV(0,1) pro i — 1, .. ., n jsou stochasticky nezávislé, pak náhodná veličina C — Eľ=i U? ~ X2(ra) má Pearsonovo "chí kvadrát" rozdělení o n stupních volnosti; Xa(n) ■ ■ ■ a-kvantil pro C. Platí Ci ~ X2(rai) Pro ' — 1; • • •; m stochasticky nezávislé =>• C — YľiĹi ~ X2(ral+ra2+• • -t-"m)- 3. Studentovo t rozdělení: Nechť U ~ N(0, 1) a C ~ X2(^) Jsou stochasticky nezávislé, pak náhodná veličina T — Y ~ ť(Ar) má Studentovo ť rozdělení o k stupních volnosti; 4. Fisher-Snedecorovo F rozdělení: Nechť Ci ~ X2(rai) a C2 ~ X2(ra2) jsou stochasticky nezávislé, pak náhodná veličina F — ~ F{ni,n2) ma Fisher-Snedecorovo i*1 rozdělení s n\ a «2 stupni volnosti; ^„(rai, «2) • • • a-kvantil pro F. časové rady 5 1.7. Prostor L2(Q,A,P). L2(Q,A,P) definujeme jako množinu všech (komplexních) náhodných veličin nad týmž pravděpodobnostním prostorem (£l,A, P), které mají konečné druhé momenty (resp. rozptyly - viz dále 1.11), tj. L2(Q,A, P) :={x\x náhodná veličina nad (£l,A, P), E|X|2 < oo}. Poznamenejme, že do tohoto prostoru zahrnujeme také všechny konstanty z (C, které považujeme za náhodné veličiny s nulovým rozptylem. Věta 1.8. L2(Q,A,P) je Hilbertův prostor se skalárním součinem (x,y) — EXY a normou \\X\\2 = ^X) = ^Ě\Xf. Důkaz. L2(Q, A, P) je obdobou funkcionálního prostoru L2(d) tvořeného funkcemi absolutně integro-vatelnými v kvadrátu na intervalu 3 C M. Totiž E|X|2 — fn |X(w)|2 dP(u>), takže namísto s Lebes-gueovým integrálem pracujeme s obecnějším pojetím integrálu, kde Lebesgueova míra je nahrazena pravděpodobnostní mírou P: • Skalární součin (X,y) existuje a je konečný pro každé X,y £ L2(Q,A,P), jak snadno nahlédneme z nerovnosti 4\XY\ = (\X\ + |Y|)2 - (\X\ - |Y|)2 < (\X\ + |Y|)2 + (\X\ - |Y|)2 = 2(\x\ + |Y|)2, odkud užitím \y\ — \y\ dostáváme l^l<^(l^|2 + l^|2), takže \vxy\< f \X(lo)yJJ)\ dP(Lo) < \ \ f \X(lo)\2 dP(Lo)+ f \Y(lo)\2 dP(u) < oo. Jet * Uct Jet • L2(Q, A, P) je vektorovým prostorem. Je uzavřený na násobení skaláry c G C, neboť E|cX|2 — |c|2E|X|2 < oo. Uzavřenost vzhledem ke sčítání plyne z: \X + Y|2 < (\x\ + \y\Ý = \X\2 + 2\xy\ + |y|2 =>• F,\X + Y|2 < E|X|2 + 2F,\xy\ + E|Y|2 < oo. • Ověření, že L2(Q, A, P) je úplný, neboli Hilbertův prostor, je složitější, ale provádí se opět zcela analogicky jako v případě funkcionálního prostoru L2(d). Podrobnosti lze nalézt například v monografii [1, §2.10]. □ Důsledek 1.9 (Schwarzova nerovnost). \vxy\, \-exy\ < ||X||2||Y||2 - ^ĚW^Ělr]2, x,y G L2(Q,A,P). Důsledek 1.10. X £ L2(Q,A,P) =>• EX existuje. Důkaz. \vx\ = \e(1.x)\ < VW? VWF = VW< <». □ Důsledek 1.11. X,Y G L2(Sl,A,P) =>• X-EX, Y - EY G L2(Sl,A, P) (X - EX, Y - EY) = E(X - EX)(Y - EY) = cov(X, Y) existuje a splňuje Schwarzovu nerovnost \cov(X, Y)| < y/E\X - EX|2^/E|Y - EY|2 = • X, Y nekorelované, avšak nikoliv naopak: X,Y nekorelované X,Y stochasticky nezávislé. časové rady 7 2.4. IDENTIFIKACE PERIODICKÝCH KOMPONENT. PERIOD O GRAM. 2.4.1. Fourierova řada (FŘ). Komplexní harmonické kmity: E = {ek{t)}k£Z, ek(t) = -^e^ = ^^kt = -^=(cosujkt + ismujkt), kde pro k ^ 0 značí fk — M ... frekvenci |Ar|-ho harmonického kmitu s periodou yj^-, uk — "ŽTfk ■ ■ ■ Jeno úhlovou rychlost. E je ortonormální systém v Ĺ2([0,T]): (ej,ek) = Jo e*^1 dt = SJik = { J ^ ^.jj _B je dokonce úplný (báze) v Ĺ2([0,T])(bez důkazu): oo x(t) G L2([0,T]) libovolná a T-periodická =>• a;(ť) — xkek(t), k = — oo kde řada konverguje (bez ohledu na pořadí) dle normy 11 _ 112 v prostoru L^^O, T]) a souřadnice xk jsou jednoznačně určeny vztahem 00 00 (x: ek} — ^ *y ^ Xjej, ek} — ^ ] £j £k) — %k• J--00 j=-00 " T" ' Komplexní tvar rozvoje do Fourierovy řady: 00 1 00 k = — oo kde (2.4.1) Cfe :=<*(*) - ^=(x,ek) = -^= J x{t)-j=e-i^tdt=^ x^e^-* dt. ck ... k-tý komplexní Fourierův koeficient funkce x, co ... střední hodnota (stejnosměrná složka) funkce x, {ck]k£i ■■■ komplexní fourierovské spektrum funkce x, xk(t) := ckeiUkt + c_ke-iUkt, k > 1 fc-tá harmonická komponenta funkce x o frekvenci fk. Tvrzení 2.4.1.1 (Dirichletova věta). Jestliže T-periodická funkce x(t) má konečnou variaci na [0, T], pak její FŘ konverguje bodové k x (t) v každém bodě t, kde je x(t) spojitá a k \(x(t + 0) + x(t — 0)) v každém bodě nespojitosti (u funkce s konečnou variací je to vždy pouze konečný skok). Je-li x(t) reálná funkce, lze (2.4.1) upravit na další dva ekvivalentní tvary: Zřejmě c_k — cit- Označíme-li Ck - 7^{ak - ibk); ak, bk G M, pak cq — %■ G M, 60 — 0 a dostáváme vítězslav veselý Goniometrický tvar Fourierovy řady: £(ť) = c0 + J2(c^ÍWkt + c-ke-iwkt) k-1 cke k=i - tt + y^(afe cosUkt + bk sin ujkt), k=i kde ak — 2Re(cfc) — |r f^ x(t) cos u>kt dt bk - -2lm(ck) — |r f^ x(t) sin cjf-t dt. Položíme-li — Ak cos cpk a bk — Ak sin ipk, dostaneme Amplitudově-fázový tvar Fourierovy řady: (2.4.2) x(t) - y + ^2 Ak(cosípk cosuikt + simpk sinuikt) = -y + ^Ak cos(uikt - ipk) k = l k = l "f" ' kde (2.4.3) fk Ak fk c0 = f {Ak}k£j$ {fk}k£N ^\ = y/4+bl bk — arctan ak amplituda k-té harmonické složky funkce x, fázový posuv k-té harmonické složky funkce x (—7r < ipk < ir), střední hodnota (stejnosměrná složka) funkce x, amplitudové spektrum funkce x, fázové spektrum funkce x. 2.4.2. Parsevalova identita (PI). Nechť x £ Í2([0, T]) (tj. s konečnou energií na intervalu [0, T]) je T-periodická, potom Odtud pak dostáváme: Komplexní tvar Parsevalovy identity: 1"M2= 77T / \x{t)\2dt (FŘ je konečná: fmax < 4ry udává konečný frekvenční obsah x(t)), {ck}k£% je Ař-periodická posloupnost (ck — Ck+mN, m G Z). Po úpravě: -i [ir] JV-i £n = X! Cfce8'2^ í^c^t = X cfcei2*\ (2.4.8) k-0 k-0 kde první ze sum jsme upravili záměnou sčítacího indexu r — k + N na tvar: -1 JV-l JV-l s využitím vztahů - [^fi] + JV = [f ] + 1 a cr = cr-N pro r = [f ] + 1,.. ., N - 1 (Ař-periodicita). Označíme-li x — (aľo, a?i, .. ., kjv-i)t vzorky a?(í) na její jedné periodě a c— (co, ?i,. .. ,Cjv_i)T odhady Fourierových koeficientů, pak (2.4.8) určuje tzv. operátor diskrétní Fourierovy transformace (a; — DFTjy-(c)) dle následující definice. 10 vítězslav veselý Definice 2.4.3.1 (Diskrétní Fourierova transformace (DFT)). DFT^ :CN^>-CNje lineární operátor X - W%x určený maticí N x N: W^- — [W^-n] 0. = J2 WknWňn' = J2 wnr~s) = £ «"> 1 = wns- n = 0 n = 0 n = 0 Odtud N pro r — s q_1 — 0 pro r ^ s neboť = W^(r_,) = lag^lpror^sv důsledku 0 < \r - s\ < N - 1. □ Poznámka 2.4-3.3. V systému MATLAB jsou operátory DFT^- realizovány pomocí algoritmu tzv. rychlé Fourierovy transformace implementovaných v procedurách f f t a if f t takto: X_ = W^x = DFT^(a;) = fft(x) a x = (W^)-1^ = lw+X = ^DFT+0*) = ifft(X). Důsledek 2.4.3.4. Podle (2.4.8), věty 2.4.3.2 a poznámky 2.4.3.3 platí x = W^a=DFT^(c) = jVifft(g), c = (W+)-1* - -W],x = —DFT^-(a?) = -ttt(x) 2.4.4. Periodogram. Periodogram — odhad energetické spektrální hustoty (2.4.6) reálné T-periodické funkce x(t), kde za Ck, k — 1,... ,m, m := [ ], dosadíme odhad Ck z (2.4.7) spočtený pomocí DFT^- dle 2.4.3.4 po ekvidistantní diskretizaci jedné periody x(t): T - NAt, xn - x(nAt), n = 0, 1,. .., N. Periodogram je tedy posloupnost hodnot kde h_ - I(uk) - 2T|cfc |2 - 2NAt a JV-l JV Xk = 2^ xte * " = 2_^xte 1 ™ t=o t=i Druhá rovnost je zde důsledkem jV-periodicity xq — xjy. N Xk = ^At\Xk\", k=í, (2.4.9) časové Rady 11 Z hlediska kvalitativních závěrů nezáleží na multiplikativní konstantě, proto bez újmy na obecnosti můžeme předpokládat At — 1 a dostáváme tak výsledný vztah pro hodnoty Ik periodogramu ve tvaru 2, 2 I /A 2nkt\ . 2Tvkt\\ , Ik = I(uk) =-\Xk\2 = - l í^cos—1 +1^8!!!—1 \,k = l,...,m. (2.4.10) Velká hodnota Ik indikuje velkou energii k-té harmonické komponenty, tj. silné zastoupení složky xk(t) — cketWkt + c_ke~tWkt — Ak cos(w;-ť — cpk) v rozvoji x(t) do Fourierovy řady. Věta 2.4.4.1. N - 1 4-2 J2 pro k -1,2,..., \h\.e-^j I Protože dle důkazu věty 2.4.3.2 platí N N ^2 e~t3^~ = et3^~ = 0 pro k ^ O(modAř), t-i lze psát \s = l / \t=l — }_^}_^{xs-x){xt-x)e N s-lt-1 N-h 2 -J7^2(xt+h-x)(xt-x)e- N \h\ 0 pro j - 1, 2, .. ., m, je 0 < W < 1. Nechť Xt =Pt + Et, Et ~ WN(0, er2) gaussovský. Platí-li nulová hypotéza Ho : Xt — Et, pak testová statistika W má na intervalu [0,1] rozdělení, pro nějž platí m / "v p(w>a;)= (-iy-1(J)(i-jxr-1,o<*0 přičemž P{W > x) m(l — :c)m_1 je dobrá aproximace pro m < 50. Nechť (?_f(1 — a) značí (1 — a)-kvantil rozdělení statistiky W (tabelováno), tj. P(W < <7_f(1 — a)) — í — a. Pak ířo zamítáme s rizikem a, pokud > <7_f(1 — a). Je-li v takovém případě Iko — maxfc-i) m Ik, pak fco-tou harmonickou komponentu považujeme za významnou. Další významnou harmonickou komponentu zjistíme z periodogramu délky m — l, kde vynecháme Iko, a tak pokračujeme dále, dokud není hypotéza Ho přijata. 12 vítězslav veselý 2.4.6. Siegelův test periodicity. Tento test je vhodnější než Fisherův v případě většího množství harmonických komponent. Siegelova statistika je tvaru m T\ — ^^(Yk — A (?_f(1 — «)) + , 0 < A < 1 (doporučená hodnota je A — 0, 6). k-l Ho zamítáme s rizikem a, jestliže T\ > t\(í — a), kde t\(í — a) je (1 — a)-kvantil rozdělení T\ (tabelováno). časové Rady 13 2.5. SEPARACE JEDNOTLIVÝCH SLOŽEK. Označení. t — (t\,. .., t„)T, t{ ^ t j pro i ^ j ... vektor "času" , X :— Xt — (Xtl,.. ., ^t„)T • • • náhodný vektor populace vybrané z časové řady, x .— xt — (a?i, .. ., Kn)T, #i — . .. pozorování vektoru populace X. A. Parametrické metody (P). Tyto metody používáme v případě, že je znám analytický tvar Dt — D(t; fli,.. ., /3p) v aditivním modelu Xt — Dt + Et (multiplikativní model převedeme na aditivní logaritmováním) s neznámými parametry f3 — (j3i,. .., j3p)T, které odhadneme minimalizací výrazu N(/3i, .. ., j3p) — \\D(t; f)\, .. ., j3p) — x\\2 ve vhodné normě || • ||. Zpravidla se používá euklidovská norma a klasická metoda nej menších čtverců: N(l31,...,l3p) = J2\D(U,l3i,...,/3p)-xi ! = 1 Lokální extrém hledáme řešením systému rovnic 0 pro j = 1,.. .,p s ověřením podmínky pro lokální minimum (jakobián je pozitivně definitní matice): \d2N (ä,..., ppy j dfc dj3j > 0. i, 3 Pokud D závisí lineárně na . .. ,/3p, pak tento systém je systémem p lineárních rovnic pro p neznámých f)\,. .., f)p a dostáváme klasický lineární regresní model. V případě nelineární regrese použijeme pro minimalizaci N(fli,. .., /3P) vhodnou optimalizační metodu, například Optimization Toolbox v prostředí systému MATLAB. Lineární regresní model: D(t; fli, .. ., /3p) — ípi(t)/3i + • • • + ípp(t)/3p, kde ip{ jsou dané. X = D/3 + E, E= (Etl,...,Eu)T ~IID(0,). B. Neparametrické metody (NP). Za vhodných předpokladů konstruujeme lineární nebo nelineární operátor s (smoother — vyhlazo-vač) takový, že s({xt})=: je dobrý odhad Dt v tom smyslu, že {Dt} —¥ {Dt} ve vhodné konvergenci při n —¥ oc, tj. odhad je asymptoticky přesný. Značení metod: AM-P ... parametrická 1 , . , , , _ _ metoda v aditivním modelu. AM-N P ... neparametrická J MM-P ... parametrická 1 , . ,,. ... ,. . . „t, • n / ř metoda v multiplikativním modelu MM-NP ... neparametrická J 14 vítězslav veselý 2.5.2. Metody pro odhad sezónní složky. x — (xi, .. ., xn)T . .. ekvidistantní pozorování náhodného procesu X — {Xt\t £ 7L} e — (ei,. .., en)T . .. pozorované chyby í G{l,2,...,n} Předpokládáme model: - aditivní: Xt = Trt + Szt + Ct + Et (AM) nebo - multiplikativní: Xt = TrtSztCtEt (MM). Nechť d je perioda Szt a n — md, pak zavedeme matici m x d X = X\ . . . Xk . . . X d Xd+l ■ ■ ■ Xd+k ■ ■ ■ X2d = YxiM kde Xjk — £(j_i)d+fc (např. j =rok, fc=mésíc). Označíme-li vec operátor, který skládá sloupce matice pod sebe do jednoho sloupcového vektoru, pak můžeme psát x — vec(XT). V MATLABu: X = reshape(a;, d, m).' Řádek X(j, :) představuje pozorování hodnot j-té periody Szt spolu s Trt, Ct a Et. Sloupec X{\, k) představuje pozorování k-té hodnoty v každé periodě Szt spolu s Trt, Ct a Et. Nechť E — [ejtk] je analogicky zavedená matice chyb. platí pro s •— Sí+---+Sd Hledaná složka Szt je určena hodnotami s :— (si,.. ., Sd) :— (Szi, .. ., Szd) jedné své periody, přičemž T - . J v AM .„ _ A. s = i i (2.5.4) 1 1 v MM v ; Budeme hledat odhady % splňující (2.5.4). 2.5.2.1. Metoda malého trendu pro odhad Szt (AM-NP, MM-NP). Předpoklad: Trt + Ct, resp. TrtCt je přibližně konstantní v každé periodě Szt: v AM: Xjk-Sk-ejk - Trt + Ct & m j \ . ,,,, // \ m n > v i-te periodě, ť = (i — lid + Ar. v MM: xjik/(sk ejtk) = TrtCt Rí rrij j j v ' yj ' Algoritmus: (1) fhj - i J2k=i xiM ■ ■ ■ odhad mj Pro 3 = 1, (2) pro Ar = 1, .. ., d spočteme: vAM: sk = ^EjLi (xi* v MM: sk = ^E™=i^W%, Odhady zřejmě splňují (2.5.4): časové Rady 15 v AM: v MM: ^ a i \ k=l \ j=l 1 \ ^ I l_ \ ^ Xjtk k=i \ j=i 3 ) 1A1 i - — 1^— jL xiM - —m-l. j = l J fc = l Dále pokračujeme kroky v 2.5.2.3. 2.5.2.2. Metoda klouzavého průměru pro odhad Szt (AM-NP, MM-NP). Zavedeme vektor vah w — (w\,.. ., u>2q+i)T pro výpočet obyčejného průměru takto: w — ^(1, 1, „ ., ť)T pro periodu liché délky d — 2q + 1, (2g+l)x w — ^(1/2, 1, ..., 1,1/2)T pro periodu sudé délky d — 2q. (2g+l)x Spočtěme klouzavý průměr mt J29r=-q wTxt+T pro q + 1 < t < n - q, xt pro 1 < t < q nebo n — q + 1 0, t = 1, 2,. .., n. Tento model je vhodný pro popis trendu s konstantním koeficientem růstu Trt+2 - Trt+i = apt+2 - a/?ŕ+1 Trt+Í - Trt af?+1 - a/3* 1 ' kde při a > 0 dostáváme pro j) > 1 růst a pro 0 < j) < 1 pokles (obr. 2.5.1). V případě a < 0 je tomu naopak. Algoritmus: Zlogaritmováním provedeme transformaci na aditivní model s lineárním trendem: InTrt = lna + ťln/?, /3 = (ln/?, Ina)T. Odhad /3 spočteme metodou lineární regrese dle odst. 2.5.1.1. Odlogaritmováním hodnot získaných podle vztahů (2.5.1)-(2.5.3) pak dostaneme výsledky pro původní model. Kontrolní kritérium adekvátnosti tohoto modelu: t+1 ftá /?, resp. A ln xt — ln xt+i — ln xt ftá ln j3. xt Kontrolujeme tedy, zda tato veličina osciluje náhodně kolem konstantní úrovně. 2.5.3.8. Modifikovaný exponenciální trend (AM-P). Model: Trt = ^ + a/3t, a, /?, 7 G M, f} > 0, t= 1, 2,. .., n. Jedná se o tříparametrické zobecnění předchozího modelu. Používá se především v případě, kdy trend má konstantní koeficient růstu, resp. poklesu a přitom je shora, resp. zdola asymptoticky omezen. Tehdy volíme 0 < /? < 1, přičemž 7 pak odpovídá asymptotické úrovni, ke které trend roste (a < 0), resp. klesá (a > 0). Viz obr. 2.5.2. 18 vítězslav veselý a) Rostoucí: a — —4, /? — 0, 95, 7 — 5 b) Klesající: a — 4, /? — 0, 95, 7=1 Obrázek 2.5.2. Modifikovaný exponenciální trend Algoritmus: • Soubor pozorování x rozdělíme na stejně velké třetiny délky m (pokud n ^ 3m, vynecháme ze začátku řady ramod3 pozorovaných hodnot). • Pozorování v jednotlivých třetinách sečteme. Užitím vzorce pro součet m členů geometrické řady S :— YltLi — PÍP™ ~ 1)/(P ~ 1) Pak obdržíme tři rovnice o třech neznámých a, /?, 7: £1 := Lt=i R* Er=irrt = m7 + aá1 £2 R* Ľt2ľm+iTrt = mj + al3mS £3 •— Z^t=2m+1 R* = mi + a/32mS. • Systém rovnic vyřešíme postupnou eliminací jednotlivých neznámých: £3 — £2 £3 - £2 = a/?raS(/?ra - 1) £2 - £1 = aS{pm - 1) Pak £2 - £1 = aS{fT - 1) =>• «=^db;-ty a £1 = m7 + 7 «5 S{j3m - 1) m • Obdrželi jsme tak přibližné vztahy pro výpočet neznámých parametrů a, /?, 7: 5 £2 — £l /?m - 1 £2 — £l S(j3m - 1) £1 -aS Kontrolní kritérium adekvátnosti tohoto modelu: ATrt+1 _ 7 + a/3t+2 - (7 + a/?t+1) _ a/?t+1(/? - 1) ATrt 7 + a/3t+1 - (7 + a/?* = f). Kontrolujeme tedy, zda veličina Aií+i/A3)t = (xt+2 — xt+i)/(xt+i — xt) osciluje náhodně kolem konstantní úrovně. časové Rady 19 2.5.3.9. Logistický trend (AM-P) Model: 7 a,/3, f G M, /? > 0, 7 > 0, t = 1,2, 1 + a/?*' Tento model dále rozšiřuje možnosti modifikovaného exponenciálního trendu, neboť má inflexní bod v t — —lna/ ln/?, který zachycuje změnu akcelerace růstu (obr. 2.5.3). Toto tvrzení nyní ověříme: . 7a(ln p)Pt „ ^ a/?* ln /? kde jsme použili 7 — Trt = 7 — 1+^gt = i+^t • Oba členy Trt a 7 — Trt působí proti sobě a příslušný lokální extrém odpovídající inflexnímu bodu nalezneme řešením rovnice Tr" = 0: 0 = Trf - -í^[Trí(7 - Trt) - TrtTr[] = -^^(7 - 2Trt) <^> 7 = 2Trt = Tq|L7 l + aPť = 2 & aP* = 1 lna + *ln/? = 0 í =-£f • Algoritmus: • Transformace časové řady: «t — ^±, kde Aij — aľt+i — aľ*. Označíme-li yt := Trt, pak Aí/t m Tr{ --^(7 ~ Vt) a tedy ^ m -ln/? + í^jft. Protože a;t = yt + et ^ — y, + a AíCt — Ai/t + Aet, dostáváme (i) -\ríp+^{xt-et) + e(2) „(2) • Odhad parametrů P a 7: Parametry odhadneme v lineárním modelu — 1 a?i kde e. (3) _ „(2) lnj ■et. + ef], kde /?i In/?, /?2 = -= —• 7 7 Pak P = e^ a7 = f-. • Odhad parametru a: Odhad a nalezneme z a;t»á 7/(1+ «/?*), odkud 7 7 a/? «J--1 ^ ln a + ť ln /? ln ---1 xt \xt Chybu této aproximace potlačíme zprůměrováním pro t — 1,2,...,«: i £(ln a + t In /?) = ln a + I (ln /?) = I £ ln( J - 1). t-i t-i * Odtud po dosazení odhadů /?i — ln/? a 7 — dostáváme výsledný odhad pro lna: ln a — — (n + l)/?i tzv. Rhodesův vztah. Tedy pak a _ „ln 8 Kontrolní kritérium adekvátnosti tohoto modelu: Křivka 1. diferencí xt+\ — xt ftá Tr£ osciluje kolem křivky podobné normální hustotě (obr. 2.5.3b)). Navíc Í/Trt se zřejmě chová jako modifikovaný exponenciální trend, takže kontrolujeme, zda í 1 ] Xt + l I Xt+l xt osciluje kolem konstantní hodnoty. 20 vítězslav veselý a) Logistická křivka b) Derivace logistické křivky Obrázek 2.5.3. Logistický trend: a = 4, /? = 0, 95, 7 = 5 a) Gompertzova křivka b) Derivace Gompertzovy křivky Obrázek 2.5.4. Logistický trend: a = - ln 160, /? = 0, 95, 7 = ln 160 2.5.3.10. Gompertzova křivka (MM-P). Model: Trt = e^+api, a, /?, 7 G M, /? > 0, í = 1, 2, .. ., n. Má podobné vlastnosti jako logistický trend, tj. opět je možno modelovat inflexi — tentokrát pro a < — 1 v bodě t — — ln(—a)/ ln /?. V tomto případě však 1. derivace Tr't není symetrická kolem bodu inflexe, aleje zešikmena doleva (obr. 2.5.4). Algoritmus: Zlogaritmováním převedeme úlohu na modifikovaný exponenciálni trend. Kontrolní kritérium adekvátnosti tohoto modelu: lnTrt se chová jako modifikovaný exponenciálni trend, takže stačí kontrolovat, zda Qn.xt+2 - lnaľt+i)/(lnaľt+i -lnaľt) osciluje kolem konstantní hodnoty. časové Rady 21 Věta 4.18 (Durbin-Levinsonův algoritmus pro rekurentní výpočet 4>„). Nechť X — {Xt \t G "Z,} je stacionární náhodný proces s nulovou střední hodnotou fix — 0 a autoko-varianční funkcí Jx(h) —± 0 pro h —¥ oo. Jestliže Xn+i — <&n^Xn + ■ ■ ■ + <^„t„Xi je nejlepší lineární predikce jako ve 4-10, pak pro koeficienty <&n,j a střední kvadratickou chybu v„ — E|Xn+i — Xn+i|2 platí následující rekurentní vztahy. Počáteční krok pro n — 0: v o - 7x(0). Rekurentní krok pro n > 0: (4.4a) = 0 pro n = l Vn - (1 - |$n,n|2), " $n,l ' $n-l,l $n,2 - $n-l,2 $n-l,n-2 . $n,n-l . . $n-l,l pro ra > 1. (4.4b) (4.4c) (4.4d) □ Důkaz. Viz [1, str. 162]. Definice 4.19 (ARMA proces). Stochastický proces X — {Xt \ t G 2Z} se nazývá ARMA procesem řádu p, q (0 < p, q < oo), píšeme X ~ ARMA(p, q), jestliže xt = Si*t_i + $2^t-2 + • • • + $P*t-P + + ©i^t-i + ®iZt-2 + ■■■ + QqZt_q, (4.5a) Autoregresní část (AR) „2\ Část klouzavých součtů (MA=Moving Average) kde Z := {Zt \ t G TL) ~ WN{0, er2) a $p ^ 0, Qq ^ 0, cr ^ 0. Přepíšeme-li (4.5a) do ekvivalentního tvaru ■ $pXt_p = Zt + 9i^_i + Q2Zt-2 + ■■■ + QgZt t — q, (4.5b) (4.5c) Xt -$2^t-2- • lze použít též zkrácený zápis $(5)Xt - Q(B)Zt nebo $(«)X(«) - 0(«)^(«) kde při $o = ©o = 1 $(2) = 1 - -----$pžP, 9(2) = l + 0i2 + • • • + 9g29, « G C. Poznámka 4-20. V předchozí definici bez újmy na obecnosti předpokládáme 0o — 1, neboť v opačném případě stačí nahradit bílý šum {Zt} ~ WN(0,(t2) šumem {@oZt} ~ WN(0, (0o2Xt-2 + ■■■ + %Xt-p + Zt. (4.5d) V tomto případě připouštíme i p — oo za předpokladu, že $ := {<^j}JL1 G i\. b) Proces klouzavých součtů (MA proces): X ~ MA(q) = ARMA(0, q) : Xt = @{B)Zt, neboť $(2) = 1. Tedy xt = Zt+ BiZt-i + QiZt-i + ■■■ + QqZt-q. (4.5e) V tomto případě připouštíme i q — oo za předpokladu, že 0 := {Qj}'jL1 G i\. 22 vítězslav veselý c) Bílý šum (jediný proces, který je současně AR i MA procesem): X ~ ARMA(0, 0) = AR(0) = MA{0) = WN{0, er2) : Xt = Zt d) Smíšený ARMA proces: X ~ ARMA(p, q), 0 < p, q < oo: Netriviální směs autoregrese a klouzavých součtů. Definice 4.22. Buď X = {Xt \ t G 7L}, X ~ ARMA(p, q). X se nazývá kauzální ARMA proces, jestliže existuje ip — {^jlJLoi YľjLol^jl ^ 00 ^ ^ ^i) ^e oo Xt = Í>jZt-j (zkráceně Xt = ip(B)Zt), t G X. (4.6a) X se nazývá invertibilní ARMA proces, jestliže existuje 7r — {^j}JL0, ^JLol^jl 00 OJ- n ^ ^i) tak, že oo ^TTjXt-j = Zt (zkráceně n(B)Xt = Zt), t G X. (4.6b) j=o Poznámka 4-23. Tedy kauzalita znamená, že ARMA proces X je kauzálním lineárním procesem, který lze generovat bílým šumem {Zt}, neboli X ~ MA(oo). Podobně invertibilita znamená, že bílý šum {Zt} je kauzálním lineárním procesem, který je generován ARMA procesem X, což je ekvivalentní s X ~ AR(oo). Výše jsme použili také skutečnost, že ipo — tto — 1 (viz dále odst. 4.32, který se týká výpočtu kauzální, resp. invertibilní reprezentace). Věta 4.24 (Autokovarianční funkce MA procesu). {Xt} ~ MA(q), q < oo je stacionární náhodný proces s nulovou střední hodnotou fix — 0 a autokovarianční funkcí i 7x(/i) - 0. (4.7a) k=o Speciálně pro q < oo je lx{h) = {(ľ2^@h+k®k P™U q a zejména fx{l) — ""29g ^ 0, neboť 0O = 1. Důkaz. Podle 4.4 je {Xt} stacionární a platí i fJ-x — fJ-z 9j — 0, neboť fiz — 0. j = 0 lx{h)= J2QjQkMh-j+k),káelz{h) = r wohhZl- V sumě zůstanou nenulové členy jen pro h — j + k — 0, tj. pro j — h + k, kdy tedy dostáváme (4.7a) i -tx{h) = Y,®h+k®k-fz{0). Pokud je q < oo, pak @h+k — 0 pro h + k > q, tj. pro k > q — h, takže (4.7a) lze psát ve tvaru (4.7b). □ časové Rady 23 Důsledek 4.25. <4=7x(0) = ,|2. tťo (4-8a) <4 ^ <^2(l + lQi|2 + l©2|2 + • • • + +|6g|2) pro q < oo. PxW-Eff 9*+*f* profe>0. (4.8b) Věta 4.26 (Parciální autokorelační funkce kauzálního AR procesu). Nechť {Xt} ~ AR(p), p < oo je kauzální AR proces. Pak {Xt} je stacionární s nulovou střední hodnotou fix — O a parciální autokorelační funkcí c*x, pro niž ctx(p) — $p ^ O a c*x(h) — O pro h> p. Důkaz. S uvážením 4.23 je v důsledku kauzality {Xt} ~ MA(oo) a tedy stacionární s nulovou střední hodnotou podle 4.24. Podle (4.6a) je Xt — E^lo ^j^t-j a tedy Xt G £>(Zt, Zt-i,. ..) —: Ht (uzávěr lineárního podprostoru v L^i^l, A, P) generovaného náhodnými veličinami Zu, u < t). Položíme-li $j — O pro j > p, můžeme dle (4.5d) pro každé n > p psát Veličiny Zt jsou nekorelované, takže Zt ± Zu pro t ^ u. Zejména Zt _L Zu pro u < t a tedy Zt _L £>t-i s využitím spojitosti skalárního součinu v L^^ljA, P). Podle věty o ortogonální projekci je Xt — Yľj-i ®jXt-j jednoznačně určená nejlepší lineární predikce Xt pomocí veličin Xt-\, .. .,Xt-n pro každé n > p. Podle věty 4.13 pak a(p) — $p a a(n) — $n — O pro n > p. □ Na obrázcích 4.1, 4.2 a 4.3 jsou typické průběhy odhadnuté autokorelační a parciální autokorelační funkce simulovaných procesů po řadě AR(2), MA(2) a ARMA(2, 2). Čárkovaný pás udává interval spolehlivosti 95% pro nulovou hodnotu. Vidíme, že v případě procesu AR(2) na obr. 4.1, resp. MA(2) na obr. 4.2 skutečně ax(h) Pá O, resp. px(h) Pá O pro h > 2 v souladu s větami 4.26, resp. 4.24. V ostatních případech obálka px(h), resp. ctx(h) (v případě ARMA(2, 2) na obr. 4.3 dokonce obou) vykazuje exponenciálně klesající, případně navíc i oscilatorický, průběh. Dá se totiž ukázat, že p x i ctx jsou v těchto případech lineární kombinací klesajících geometrických posloupností a kosinusovek s geometricky klesající amplitudou. vítězslav veselý Pozorování procesu, n — 500 50 XOO X50 200 250 300 350 400 450 500 Autokorelační funkce XO 15 20 25 30 35 40 45 50 Parciální autokorelační funkce 5 XO 15 20 25 30 35 40 45 50 Obrázek 4.1. Proces Aff(2) : * = [0.5, 0.2], er2 = 2.25 časové Rady Pozorování procesu, n — 500 50 XOO X50 200 250 300 350 400 450 500 Autokorelační funkce XO 15 20 25 30 35 40 45 50 Parciální autokorelační funkce XO 15 20 25 30 35 40 45 50 Obrázek 4.2. Proces MA(2) : ® = [-0.5, -0.2], er2 = 2.25 vítězslav veselý Pozorování procesu, n — 500 50 ÍOO 150 200 250 300 350 400 4 50 500 Autokorelační funkce T XO 15 20 25 30 35 40 45 50 Parciální autokorelační funkce r O 5 XO 15 20 25 30 35 40 45 50 Obrázek 4.3. Proces ARMA(2, 2) : $ = [0.5,0.2], 0 = [-0.6,0.3], er2 = 2.25 časové Rady 27 5.2.3. Postup hledání SARIMA modelu pro pozorovanou časovou řadu. 1. Odstranění případné heteroskedasticity Viz 5.1.3/1. 2. Určení sezónní periody s Za s zvolíme společnou periodu nejvíce dominantních periodických složek, neboli nejmenší společný násobek jejich period. Dominantní periodické složky nejsnáze identifikujeme z periodogramu. Pokud společná perioda vychází příliš dlouhá, volíme za s periodu dominantní složky s nej větším počtem jejích harmonických komponent indikovaných v periodogramu jako významných, tj. jejich frekvence jsou (alespoň přibližně) celočíselnými násobky frekvence této dominantní složky. Je přitom třeba dávat pozor, aby byl periodogram spočten z vhodně zkrácené časové řady, tak aby v ní byl obsažen celý počet period délky s. Periodu s výrazně dominantní periodické složky lze zpravidla dosti přesně odečíst i jako periodu oscilací v odhadnuté autokorelační funkci px{)- 3. Určení řádů diferencování d a D (a) Analogicky jako v 5.1.3/2(a) se snažíme zvolit nejmenší d a D tak, aby diferencovaná řada AdAfXt se jevila jako stacionární. V praxi je zpravidla 0 < d, D < 1. Zejména případ D > 1 je velmi málo pravděpodobný, neboť odpovídá proměnlivému kolísání sezónních amplitud, které však bývá většinou spojeno s obdobně proměnným rozptylem. Jeho odstranění vhodnou transformací ad 1 pak obvykle odstraní i kolísání v amplitudách (viz soubor airpass.m). (b) Zkoumáme chování odhadnuté autokorelační funkce pw{) diferencované řady. Pokud ještě vykazuje oscilace s periodou s (lze samozřejmě citlivěji ověřit opět z periodogramu diferencované řady), pak je třeba řád D ještě zvýšit (oscilace s jinou periodou svědčí pro použití zobecněného SARIMA modelu popsaného v 5.2.5). Pokud pw(h) klesá pro h — 1,2,. .., | pomalu (lineárně) a nikoliv exponenciálně, pak je třeba ještě zvětšit řád d. (c) Analogicky jako v 5.1.3(c) hledáme diferencovanou řadu s nej menším rozptylem. 4. Určení řádů P, Q, p, q (a) rozborem autokorelační a parciální autokorelační funkce pw{) a ®w{) Dle 4.34 musí pro h — 1,2,3,... průběh pw (hs) a aw (hs) korespondovat s modelem ARMA(P, Q) a pro h — 1, 2, .. ., s — 1 s modelem ARMA(p, q). Zpravidla se používají modely typu SARIMA(0, d, q, 0, D, Q, s) nebo SARIMA(p, d, 0, P, D, 0, s). Naopak modely SARIMA(Q, d, q, P, D, 0, s) nebo SARIMA(p, d, 0, 0, D, Q, s) se nedoporučují, neboť obvykle vedou k vysokým řádům a tedy k velkému počtu parametrů. Identifikace obecných modelů SARIMA(p, d,q,P,D,Q,s) bývá většinou dosti obtížná. (b) užitím ztrátové funkce (viz 4.36/4) 1° Nejprve určíme P, Q : (5.2b) W ~ SARIMA(p, 0, q, P, 0, Q, s) (54a) ~ ARMA(P, Q), $(B)Wt(k) = ®{B)E[k' pro k — 1, 2,. .., s. Pak P, Q nalezneme dle minima ztrátové funkce pro j Ylk=i ^t*^ případně můžeme spočíst optimální Pf., Q\. pro každý parciální proces a za P, Q vybrat dvojici, která se nejčastěji vyskytuje mezi všemi dvojicemi P/,, Qk-2° Pak určíme p, q : Minimalizujeme ztrátovou funkci pro W ~ SARIMA(p,Q,q, P,Q,Q,s), kde P, Q jsou již pevné a p, q vybíráme z dvojic ve vhodném rozmezí p — 0, 1, .. .,pmax a g — 0, 1, .. ., (?max-Pokud je P — 0 (resp. Q — 0), volíme většinou nejprve 0 — p — pmax (resp. 0 — q — (?max) pevně. 5. Odhad parametrů modelu a jeho verifikace V zásadě postupujeme jako ve 4.36/5-7. 28 vítězslav veselý 5.2.4. Konstrukce predikcí v SARIMA modelu. Přímý postup zpětného integrování analogický k 5.1.4 je algoritmicky komplikovaný vzhledem k sezónnímu diferencování As. Užijeme proto jiný přístup: V zásadě lze totiž postupovat jako ve 4.36/7(4), neboť z hlediska predikcí se {Xt} řídí ARMA modelem užít 4> . Závěrem poznamenejme, že pro výpočet 4> , ®* a c lze využít procedury ARMApar: [#*, &*, c, #*] = ARMApar ([p, d, q, P, D, Q, «],{#, <ř, 0,&}, fiW). Vlastní predikce ; pak spočteme jako c+armasim(íc, 4> , @ , ícO, 2O). časové Rady 29 LITERATURA [1] Peter J. Brockwell and Richard A. Davis, Time series: Theory and methods, 2-nd ed., Springer-Verlag, New York, 1991 (corrected 2-nd printing 1993). [2] Tomas Cipra, Analýza časových řad s aplikacemi v ekonomii, SNTL, Praha, 1986. Katedra aplikované matematiky, Masarykova univerzita v Brně, Janáčkovo nám. 2a, 662 95 Brno E-mail address: vesely@@math .muni . cz