Kapitola 1 Box-Jenkinsova metodologie 1.1 Úvod V předchozím jsme podali výklad analýzy jednorozměrných časových řad založený na dekompo- zičním principu. Pokoušeli jsme se vyčlenit (a posléze modelovat) deterministickou část časové řady (deterministický trend, deterministickou sezónní složku) a stochastickou část jsme hlouběji neanalyzovali - chápali jsme ji jako část zbytkovou. Ve skutečnosti však také v reziduální (zbytkové) části existují jisté systematičnosti, jež jsou velmi významné, a je proto nutné se jimi zabývat. Box-Jenkinsova metodologie (zkráceně B-J metodologie) ztotožňuje systematickou část časové řady s částí deterministickou a je založena na myšlence, že časová řada může být chápána jako řada stochastického charakteru. Zásluha Boxe a Jenkinse nespočívá v objevení principů, které budeme popisovat, ale ve vytvoření konkrétního postupu, jak tyto principy prakticky využívat. Téměř každý kvalitní statistický programový paket tuto metodologii obsahuje (SAS, BMDP, SPSS, STATGRAPHICS, RATS,) Výhody B-J metodologie: ˇ je flexibilní a rychle se adaptuje na změnu v charakteru modelovaného procesu ˇ v mnoha případech dává nejlepší výsledky (vzhledem k MSE-střední čtvercové chybě) Nevýhody B-J metodologie: ˇ musí být dostatečně dlouhé realizace ˇ ztrácí se možnost jednoduché interpretace výsledných modelů (lze těžko přesvědčit za- davatele, že řadu lze modelovat pomocí náhodných šoků. Jediným argumentem jsou zde často jen kvalitní předpovědi získané pomocí těchto modelů). 1.2 Základní pojmy V dalším budeme uvažovat centrované stacionární náhodné posloupnosti {Yt, t Z}, kde Yt L2 (, A, P), což je Hilbertův prostor reálných náhodných veličin s konečnými druhými momenty, ve kterém dvě náhodné veličiny X a Y považujeme za ekvivalentní, pokud P(X = Y ) = 1. 1 1.2. Základní pojmy 2 1.2.1 Operátor zpětného posunutí Definice 1.2.1. Nechť {Yt, t Z} je posloupnost náhodných veličin. Operátor zpětného posunutí je definován pomocí výrazu BYt = Yt-1, přičemž jej lze aplikovat několikanásobně jako Bj Yt = Yt-j. 1.2.2 Lineární proces Než zavedeme pojem lineárního procesu, vyslovme větu, která zabezpečuje jeho korektnost. Věta 1.2.2. Nechť {t, t Z} WN(0, 2 ) je bílým šumem, dále mějme posloupnost reálných čísel {j} j=0 takovou, že j=0 2 j < . Pak řada j=0 jj konverguje podle kvadratického středu, tj. existuje náhodná veličina Y L2 (, A, P) a platí Y = l.i.m. N N j=0 jj. Důkaz. Víme, že bílý šum t L2 (, A, P). Pro libovolná přirozená čísla k, N N platí N+k j=0 jj - N t=0 tt 2 = E N+k j=0 jj - N t=0 tt 2 = E N+k j=N+1 jj 2 = E N+k j=N+1 jj N+k h=N+1 hh = N+k j=N+1 N+k h=N+1 jh E jh nekorel. = 2 N+k j=N+1 2 j --- N 0 Posloupnost částečných součtů je tedy cauchyovská, tj. existuje k ní limita Y =l.i.m. N N j=0 jj. 2 Definice 1.2.3. Mějme {t, t Z} WN(0, 2 ) a posloupnost reálných čísel {j} j=0 takovou, že j=0 2 j < , pak lineární proces je definován vztahem Yt = j=0 jt-j. 1.2. Základní pojmy 3 Počítejme postupně střední hodnotu, rozptyl a autokovarianční funkci lineárního procesu a přesvědčeme se, že lineární proces je stacionární. EYt = E j=0 jt-j = j=0 j Et-j =0 = 0 DYt = D j=0 jt-j nekorel. = j=0 2 j Dt-j =2 = 2 j=0 2 j = 2 Y (t) = C(Ys, Ys+t) = EYsYs+t = E j=0 js-j h=0 hs+t-h = j=0 h=0 jhE s-js+t-h (t+j-h) = 2 j=0 jj+t protože (k) = 2 k = 0, 0 jinak. Ze Schwarzovy nerovnosti dostaneme |(t)| = |C(Ys, Ys+t)| = 2 j=0 jj+t DYsDYs+t = (0) = 2 j=0 2 j < . Podmínka stacionarity je tedy podmínka, že j=0 2 j < . Pokud zavedeme funkci Y (z) = j=0 jzj , pak podmínka j=0 2 j < implikuje, že funkce Y (z) je holomorfní uvnitř kružnice |z| < 1. Takže podmínku stacionarity lze vyslovit i pomocí podmínky: Y (z) je holomorfní pro |z| < 1, přičemž j=0 2 j < . Oba požadavky budou splněny, pokud bude platit Y (z) je holomorfní uvnitř a na jednotkové kružnici . 1.2. Základní pojmy 4 Lineární proces lze ještě zobecnit takto: Definice 1.2.4. Mějme {t, t Z} WN(0, 2 ) a posloupnost reálných čísel {j} j=- takovou, že j=- 2 j < , pak zobecněný lineární proces je definován vztahem Yt = j=- jt-j. Pro takto definovaný zobecněný lineární proces dokážeme obdobným způsobem jak pro oby- čejný lineární proces spočítat EYt = 0, DYt = 2 j=- 2 j a (t) = 2 j=- jj+t. Na závěr tohoto odstavce počítejme ještě spektrální hustotu zobecněného lineárního pro- cesu. Nejprve odvodíme spektrální hustotu bílého šumu, a to pomocí jeho autokovarianční funkce (t) f() = 1 2 t=- e-it (t) = 1 2 t=- e-it 2 (t) = 2 2 -, , 0 jinak kde (t) = 1 t = 0, 0 jinak. 2 2 /2 - Pak pomocí autokovarianční funkce zobecněného lineárního procesu počítáme spektrální hus- totu pro -, fY () = 1 2 t=- e-it (t) = 1 2 t=- e-it 2 j=- jj+t = 2 2 f() j=- j t=- j+te-it = f() j=- jeij t=- j+te-i(j+t) = f() t=- je-it 2 = f() t=- jeit 2 neboť |z|2 = z z Pokud položíme Y (z) = j=- jzj , pak můžeme psát fY () = f() e-i 2 = 2 2 Y e-i 2 = 2 2 Y ei 2 . 1.2. Základní pojmy 5 1.2.3 Lineární filtry Věta 1.2.5. Nechť {Xt, t Z} je (centrovaná) stacionární náhodná posloupnost a {j} j=- je absolutně konvergentní posloupnost reálných čísel (tj. j=- |j| < ). Pak platí Yt = j=- jXt-j L2 (, A, P) tj. {Yt, t Z} je stacionární náhodná posloupnost. Důkaz. Je zřejmé, že stačí dokázat existenci náhodných veličin Y (2) t = -1 j=- jXt-j a Y (1) t = j=0 jXt-j, protože pak bude platit Yt = Y (1) t + Y (2) t . Označme X(h) = EXtXt+|h| a X(0) = 2 X > 0. Pak pro libovolná přirozená čísla k, N N platí N+k j=0 jXt-j - N h=0 tXt-h 2 = E N+k j=0 jXt-j - N h=0 tXt-h 2 = E N+k j=N+1 jXt-j 2 = E N+k j=N+1 jXt-j N+k h=N+1 hXt-h = N+k j=N+1 N+k h=N+1 jh EXt-jXt-h |(j-h)|X(0)=2 X 2 X N+k j=N+1 N+k h=N+1 |j||h| = 2 X < N+k j=N+1 |j| 0 2 --- N 0. Posloupnost částečných součtů je tedy cauchyovská (podle kvadratického středu), tj. existuje k ní limita Y (1) t = l.i.m. N N j=0 jXt-j, Y (1) t L2 (, A, P), tj. Y (1) t má nulovou střední hodnotu a konečný rozptyl a je tedy stacionární. Podobně se dokáže i existence stacionární náhodné posloupnosti Y (2) t . 2 Můžeme tedy korektně vyslovit následující definici. 1.2. Základní pojmy 6 Definice 1.2.6. Nechť {Xt, t Z} je stacionární náhodná posloupnost a {j} j=- je absolutně konvergentní posloupnost reálných čísel (tj. j=- |j| < ). Pak Yt = j=- jXt-j nazveme lineárním filtrem procesu {Xt, t Z}. Věta 1.2.7. Mějme centrovanou stacionární náhodnou posloupnost {Xt, t Z} se spektrální hus- toutou fX(). Nechť {j} j=- je absolutně konvergentní posloupnost reálných čísel (tj. j=- |j| < ). Pak náhodná posloupnost Yt = j=- jXt-j je stacionární se spekt- rální hustotou fY () = fX() Y e-i 2 , kde Y (z) = j=- jzj pro |z| 1 se nazývá generující funkce filtru a Y () = Y (e-i ) přenosová funkce filtru. Důkaz. Stacionaritu jsme dokázali v předchozí větě. Označme proto Y (z) = j=- jzj pro |z| 1. S ohledem na to, že paltí Y (t) = - eit fY ()d, počítejme autokovarianční funkci Y (t) = C(Ys, Ys+t) = C j=- jXs-j, h=- hXs+t-h = j=- h=- jhC(Xs-j, Xs+t-h) = j=- h=- jhX(t + j - h) = j=- h=- jh - ei(t+j-h) fX()d = - eit j=- jeij h=- he-ih fX()d = - eit j=- jeij 2 fX()d = - eit h=- je-ih 2 fX ()d = - eit fX() Y e-i 2 =fY () d = - eit fX () Y ei 2 =fY () d. 2 1.2. Základní pojmy 7 1.2.4 Autokovarianční a autokorelační generující funkce stacionár- ních procesů Při vyšetřování vlastností autokovarianční a autokorelační funkce starionárních procesů je uži- tečné zavést následující transformaci G(z) = k=- (k)zk resp. R(z) = k=- (k)zk , které se nazývají autokovarianční, resp. autokorelační, generující funkce. Všimněme si ještě vztahu mezi spektrální hustotou a autokovarianční generující funkcí f() = 1 2 t=- (t)e-it = 1 2 G e-i . Příklad 1 : Bílý šum Mějme bílý šum t WN(0, 2 ) s autokovarianční funkcí (k)= 2 k = 0, 0 jinak , takže G(z) = k=- (k)zk = 2 a f() = 1 2 G e-i = 2 2 -, , 0 jinak . Příklad 2 : Zobecněný lineární proces Mějme posloupnost reálných čísel {j} j=- takovou, že j=- 2 j < , pak zobecněný lineární proces je definován vztahem Xt = j=- jt-j a má autokovarianční funkci tvaru X(k) = 2 j=- jj+k , takže GX (z) = k=- X(k)zk = 2 k=- zk j=- jj+k = 2 k=- j=- jz-j j+kzj+k = subst. h = j + k = 2 =G(z) h=- hzh j=- jz-j = G(z)X (z)X z-1 fX() = 1 2 GX e-i = 1 2 G e-i f() X e-i X ei = f() X e-i 2 . 1.2. Základní pojmy 8 Příklad 3 : Lineární filtr Mějme posloupnost reálných čísel {j} j=- takovou, že j=- |j| < . Dále nechť {Xt, t Z} je centrovaná stacionární náhodná posloupnost. Pak lineární filtr je definován vztahem Yt = j=- jXt-j a má autokovarianční funkci Y (k) = EYsYs+k = E j=- jXs-j h=- hXs+k-h = j=- h=- jhEXs-jXs+t-h = j=- h=- jhX(k + j - h) , takže GY (z) = k=- (k)zk = k=- zk j=- h=- jhX(k + j - h) = k=- j=- h=- jz-j hzh X(k + j - h)zk+j-h = subst. s = k + j - h = h=- hzh j=- hz-j s=- X(s)zs =GX (z) = GX(z)Y (z)Y z-1 fY () = 1 2 GY e-i = 1 2 GX e-i fX () Y e-i Y ei = fX () Y e-i 2 . 1.2. Základní pojmy 9 1.2.5 Definice ARMA procesu Definice 1.2.8. ARMA proces řádu p, q je definován vztahem Yt - 1Yt-1 - - pYt-p = t + 1t-1 + + qt-q , kde t WN(0, 2 ), přičemž pomocí operátoru zpětného chodu lze psát Yt ARMA(p, q) : (B)Yt = (B)t, kde (B) = 1 - 1B - - pBp (0 1) a (B) = 1 + 1B + + qBq (0 1). Řekneme, že {Yt, t Z} je ARMA(p, q) se střední hodnotou , jestliže {Yt - } je ARMA(p, q) proces. Speciální případy ARMA procesů nazýváme: Autoregresní proces (AR proces): Yt AR(p) ARMA(p, 0), tj. q = 0 Proces klouzavých součtů (MA proces): Yt MA(q) ARMA(0, q), tj. p = 0 1.2.6 Kauzalita Dříve než zavedeme pojem kauzality, všimněme si blíže AR(1) procesu. Autoregresní proces prvního řádu. Pro autoregresní proces prvního řádu Yt - 1Yt-1 = t postupně v k krocích upravujme Yt = 1Yt-1 + t = 1 (1Yt-2 + t-1) + t = 2 1Yt-2 + 1t-1 + t = 2 1 (1Yt-3 + t-2) + 1t-1 + t = 3 1Yt-3 + 2 1t-2 + 1t-1 + t ... = k 1 (1Yt-k-1 + t-k) + k-1 j=0 j 1t-j = k+1 1 Yt-k-1 + k j=0 j 1t-j (1) Uvažujme nejprve případ, kdy |1| < 1 a {Yt, t Z} je stacionární, tj. Yt L2 (, A, P) a EY 2 t < , pak Yt - k j=0 j 1t-j 2 = k+1 1 Yt-k-1 2 = E k+1 1 Yt-k-1 2 = 2k+2 1 0 E |Yt-k-1|2 =2 Y < --- k 0 1.2. Základní pojmy 10 tj. j=0 j 1t-j konverguje podle kvadratického středu k Yt a můžeme psát Yt = j=0 j 1t-j , tj. jde o lineární proces, tedy i o lineární filtr s generující funkcí filtru (pro |z| 1) AR(1)(z) = j=0 jzj = j=0 (1z)j = 1 1 - 1z protože j = j 1 a |1| < 1. Pak dokážeme spočítat EYt = E j=0 j 1t-j = j=0 j 1Et-j = 0 DYt = D j=0 j 1t-j nekorel. = j=0 2j 1 Dt-j = 2 j=0 2j 1 = 2 1 - 2 1 (t) = C(Ys, Ys+|t|) = EYs, Ys+|t| = E j=0 j 1s-j h=0 h 1s+|t|-h = j=0 h=0 j 1h 1Es-js+|t|-h = j=0 h=0 j 1h 1 (h - |t| - j) = 2 j=0 j 1 j+|t| 1 = |t| 1 1 - 2 1 2 . (t) = (t) (0) = |t| 1 tzv. ACF funkce. Pomocí generující funkce filtru dokážeme snadno spočítat i spektrální hustotu fAR(1)() = 2 2 AR(1) e-i 2 = 2 2 1 |1 - 1e-i|2 = 2 2 1 |e-i(ei - 1)|2 = 2 2 1 |ei - 1|2 . (2) Dále řešme případ, kdy |1| > 1. Použijeme-li vztah Yt-1 = 1 1 Yt - 1 1 t a postupně v k krocích budeme upravovat Yt = 1 1 Yt+1 - 1 1 t+1 = 1 1 1 1 Yt+2 - 1 1 t+2 - 1 1 t+1 = 1 2 1 Yt+2 - 1 2 1 t+2 - 1 1 t+1 = 1 2 1 1 1 Yt+3 - 1 1 t+3 - 1 2 1 t+2 - 1 1 t+1 = 1 3 1 Yt+3 - 1 3 1 t+3 - 1 2 1 t+2 - 1 1 t+1 ... = 1 k+1 1 Yt+k+1 - k j=0 1 k+1 1 t-j . 1.2. Základní pojmy 11 Stejně jako v předchozím případu lze snadno ukázat, že - j=0 1 j 1 t-j konverguje podle kva- dratického středu k Yt. Avšak vidíme, že Yt zde vyjadřujeme pomocí budoucích hodnot {s, s > t}. Tím porušujeme přirozenou podmínku, že Yt je na budoucnosti nezávislá a říkáme, že není kauzální. (3) V případě, že platí |1| = 1, pak AR(1) není stacionární, jde o tzv. náhodnou procházku. Nyní již můžeme zavést pojem kauzality. Definice 1.2.9. ARMA proces Yt ARMA(p, q) se nazývá kauzální, jestliže existuje absolutně konver- gentní posloupnost reálných čísel = {j} j=0, (tj. j=0 |j| < ) tak, že Yt = j=0 jt-j, tj. zkráceně Yt MA() : Yt = (B)t. Poznámka 1.2.10. Protože platí j=0 2 j j=0 |j| 2 < , pak kauzální proces Yt = j=0 jt-j je lineárním procesem. Protože lineární proces je stacionárním procesem, je kauzální ARMA proces Yt ARMA(p, q), kde j=0 |j| < , také stacionárním procesem. Autoregresní proces p -tého řádu: AR(p) : (B)Yt = t Poznámka 1.2.11. Mějme polynom (z) = 0 - 1z - - pzp a nechť 1 j jsou jeho kořeny, tj. 1 j = 0. Pak platí 0 - 1z - - pzp = p j z - 1 j = 0 j (1 - jz) , v našem případě 0 = 1 a p = 0. Proveďme tedy rozklad polynomu (z) na součin kořenových činitelů (z) = (1 - 1z)p1 . . . (1 - kz)pk , kde z01 = 1 1 , . . . , z0k = 1 k jsou rozdílné (reálné či komplexní) kořeny polynomu (z), p1, . . . , pk je jejich násobnost (přičemž platí p1 + + pk = p). 1.2. Základní pojmy 12 Budeme hledat takovou absolutně konvergentní posloupnosti čísel = {j} j=0, tak aby Yt = j=0 jt-j byl kauzální proces. Pokud použijeme operátor zpětného chodu, můžeme psát: (B)Yt = t, přitom hledáme (B) takové, aby platilo (B)(B) = 1 nebo (z)(z) = 1 čili (z) = 1 (z) . Z věty o rozkladu na částečné zlomky dostáváme (pokud pro názornost předpokládáme, že všechny kořeny jsou jednoduché) 1 (z) = 1 (1 - 1z) . . . (1 - pz) = c1 1 - 1z + + cp 1 - pz pro vhodná c1, . . . , cp. Pokud pro k = 1, . . . , p platí |kz| < 1, můžeme psát ck 1 - kz = ck j=0 (kz)j a dokázali jsme najít konvergentní řadu (z) = j=0 jzj = p k=1 ck j=0 j kzj = j=0 c1j 1 + + cpj p zj , přičemž j = c1j 1 + + cpj p, neboť (z) = 1 (z) je holomorfní pro |z| 1 pouze když |1| < 1, . . . , |p| < 1 1 1 > 1 |z01|>1 , . . . , 1 p > 1 |z0p|>1 , tedy všechny kořeny polynomu (z) musí ležet vně jednotkové kružnice. Tím jsme ukázali, že existuje řešení Yt = j=0 jt-j tzv. stochastické diferenční rovnice Yt - 1Yt-1 - . . . - pYt-p = t t WN(0, 2 ) (1.1) a tímto řešením je kauzální autoregresní posloupnost řádu p. Protože Yt je lineární proces, je toto řešení stacionární. Podmínka týkající se kořenů polynomu (z) je podstatná. Lze ukázat, že v případě, kdy alespoň jeden kořen polynomu (z) leží uvnitř nebo na hranici jednotkové kružnice, neexistuje kauzální posloupnost {Yt, t Z} splňující stochastickou diferenční rovnici (1.1). Snadno se dá ukázat, že toto řešení je jediné. 1.2. Základní pojmy 13 Střední hodnota, rozptyl, autokovariance a autokorelace AR(p) Pro kauzální AR(p) procesy počítejme nejprve EYt = E j=0 jt-j = j=0 jEt-j = 0 Abychom mohli spočítat rozptyl kauzálního AR(p) procesu, nejprve rovnici Yt = 1Yt-1 + + pYt-p + t vynásobíme výrazem Yt a spočítáme střední hodnoty obou stran, tj. EY 2 t = 1EYt-1Yt + + pEYt-pYt + EtYt. (A1) Protože EYt = 0, pak autokovarianční funkce je rovna (j) = C(Yt, Yt-j) = E(Yt - EYt)(Yt-j - EYt-j) = EYtYt-j a rozptyl (0) = EY 2 t = DYt. Dále spočtěme EYtt = E j=0 jt-j t = j=0 jEt-jt = j=0 j(j) = 2 , neboť (j) = 2 j = 0, 0 jinak. a 0 = 1. Vraťme se k rovnici (A1), pak po dosazení EYtt = 2 a (0) = EY 2 t dostaneme (0) = 1(1) + + p(p) + 2 (A2) Podělme obě strany rovnice (A2) výrazem (0) > 0 a protože pro autokorelaci platí (k) = (k) (0) , dostaneme (0) =1 = 1(1) + + p(p) + 2 (0) a odtud již plyne, že DYt = (0) = 2 1 - 1(1) - - p(p) . Při výpočtu autokovariance (nebo autokorelace ACF) budeme předpokládat, že k > 0, neboť (0) = DYt již jsme spočítali. Rovnici Yt = 1Yt-1 + + pYt-p + t vynásobíme výrazem Yt-k a spočítáme střední hodnoty obou stran, tj. EYtYt-k = 1EYt-1Yt-k + + pEYt-pYt-k + EtYt-k. (A3) 1.2. Základní pojmy 14 Protože EYt = 0, je (k) = cov(Yt, Yt-k) = E(Yt - EYt)(Yt-k - EYt-k) = EYtYt-k. Spočtěme EYt-kt = E j=0 jt-j-k t = j=0 jEt-j-kt = j=0 j (j + k) =0 = 0. Vraťme se k rovnici (A3), pak po dosazení EYtt = 0 a (k) = EYtYt-k dostaneme (k) = 1(k - 1) + + p(k - p) (A4) Podělme obě strany rovnice (A4) výrazem (0) a protože (k) = (k) (0) , dostaneme tzv. Yuleovy- Walkerovy rovnice. (k) = 1(k - 1) + + p(k - p) k 1 (A5) Explicitní vyjádření autokorelační funkce procesu AR(p) Při explicitním vyjádření autokorelační funkce procesu vyjdeme z Yuleo-Walkerových rovnic (k) = 1(k - 1) + + p(k - p) k 1. Označme B(k) = (k - 1), přičemž (0) = 1 a (-j) = (j) a hledejme řešení tzv. homogenní diferenční rovnice (k) - 1(k - 1) - - p(k - p) = 0 k 1 tj. (B)(k) = 0 . Poznámka: ŘEŠENÍ HOMOGENNÍ DIFERENČNÍ ROVNICE Mějme polynom (z) = 0 - 1z - - pzp a nechť 1 j jsou jeho kořeny, tj. 1 j = 0. Pak platí 0 - 1z - - pzp = p j z - 1 j = 0 j (1 - jz), v našem případě 0 = 1 a p = 0. 1. Nechť 1 j je kořen polynomu (z), pak k j je řešením (B)(k) = 0. Důkaz: (B)k j = (1 - B - - pBp )k j = k j - 1k-1 j - - pk-p j = k j 1 - 1 1 j - - p 1 p j = 1 j k j = 0. nebo ekvivalentně: jestliže uvažujeme faktorizaci (B) = 0 i(1-iB), tak mezi faktory je i člen (1 - jB) a platí (1 - jB)k j = k j - jB(k j ) = k j - j k-1 j = 0. 1.2. Základní pojmy 15 2. Nechť 1 1 , . . . , 1 p jsou různé jednoduché kořeny, pak c1k 1 + + cpk p je řešením ho- mogenní diferenční rovnice a c1, . . . , cp jsou konstanty, které jsou určeny počátečními podmínkami. Důkaz: (B)(c1k 1 + + cpk p) = c1 (B)k 1 " 1 1 " =0 + + cp (B)k p ,, 1 p =0 = 0. 3. Je-li kořen 1 j dvojnásobný kořen, pak k j , kk j jsou řešeními (B)(k) = 0. Důkaz: Díky faktorizaci můžeme psát (B) = (1 - jB)2 k=j(1 - kB). Pak (1 - jB)2 k j = (1 - 2jB + 2 j B2 )k j = k j - 2jk-1 j + 2 j k-2 j = 0 (1 - jB)2 kk j = (1 - 2jB + 2 j B2 )kk j = kk j - 2j(k - 1)k-1 j + 2 j (k - 2)k-2 j = kk j - 2kk j + 2k j + kk j - 2k j = 0 4. Analogicky dostaneme: je-li kořen 1 j r-tého řádu, pak k j , kk j , . . . , kr-1 k j jsou řešeními (B)(k) = 0. Shrneme-li tedy předchozí, za předpokladu, že 1 1 , . . . , 1 k jsou různé kořeny s násobnostmi p1, . . . , pk, přičemž p = p1 + + pk, pak řešení homogenní diferenční rovnice (B)(k) = 0 je tvaru (k) = k j=1 pj-1 s=0 cjsks k j , kde cjs jsou konstanty, které jsou určeny počátečními podmínkami. Dále položme j = rjeij . Pak máme (k) = k j=1 pj-1 s=0 cjsks rk j eikj , Vzhledem k tomu, že platí |j| = rj < 1, dostáváme odtud, že a (k) klesá pro k exponenciálně k nule, tj. (k) --- k 0 , což je velmi důležitá identifikační vlastnost autoregresních AR(p) procesů. 1.2. Základní pojmy 16 1.2.7 Invertibilita Víme, že pokud autoregresní proces konečného řádu AR(p) je kauzální, pak jej lze vyjádřit i pomocí MA procesu nekonečného řádu, tj. AR(p) MA() . Zajímá nás, za jakých podmínek můžeme MA proces konečného řádu vyjádřit pomocí auto- regresního procesu nekonečného řádu, tj. MA(q) AR() . Tuto vlastnost nazveme invertibilitou. Přesnou definici vyslovíme později. Nejprve si všimneme jednoduchého případu, a to MA(1) procesu. MA proces prvního řádu: Yt MA(1) : Yt = t + 1t-1 t WN(0, 2 ) . Nejprve označme 1 = - a předpokládejme, že |1| = || < 1 . Využijeme-li vztahu Yt = t + 1t-1 t = Yt + t-1, můžeme postupně upravovat t = Yt + t-1 = Yt + (Yt-1 + t-2) = Yt + Yt-1 + 2 t-2 = = k j=0 j Yt-j + k+1 t-k-1 a t - k j=0 j Yt-j 2 = E t - k j=0 j Yt-j 2 = E k+1 t-k-1 2 = 2(k+1) 2 --- k 0, tedy t = j=0 j Yt-j pro || < 1 . Budeme-li předpokládat, že |1| = || > 1 , použijme vztah t-1 = 1 1 Yt - 1 1 t a označme 1 = -. Pak můžeme opět postupně upravovat t = -1 Yt+1 + 1 t+1 = -1 Yt+1 + 1 1 Yt+2 + 1 t+2 = = - k+1 j=1 1 j Yt+j + 1 k+1 t+k+1. Obdobně jako v předchozím případě, lze ukázat, že posloupnost - N j=1 1 j Yt+j konverguje pro N také k t. Tento rozvoj však nemá praktický smysl, neboť nepozorovatelná ve- ličina t je vyjádřena pomocí teprve budoucích pozorovatelných hodnot {Ys, s > t}. Dále si všimněme vícenásobné reprezentace MA(1) procesů . Demonstrovat to můžeme pomocí následujícího příkladu. Nechť platí |1| > 1 , a uvažujme dva procesy, první neinver- tibilní, druhý invertibilní: (1) Yt = t + 1t-1 t WN(0, 2 ) (2) Xt = t + 1 1 t-1 t WN(0, 2 12 ) . Ukážeme, že oba dva procesy mají stejné první a druhé momenty. EYt = E (t + 1t-1) = Et + 1Et-1 = 0 , EXt = E t + 1 1 t-1 = Et + 1 1 Et-1 = 0 , 1.2. Základní pojmy 17 Y (k) = EYtYt+k = E (t + 1t-1) (t+k + 1t+k-1) = Ett+k + 1Ett+k-1 + 1Et-1t+k + 2 1Et-1t+k-1 = (k) + 1(k - 1) + 1(k + 1)2 1(k) = (1 + 2 1)(k) + 1(k - 1) + 1(k + 1) = 2 (1 + 2 1) k = 0 12 k = 1 0 jinak X(k) = EXtXt+k = E t + 1 1 t-1 t+k + 1 1 t+k-1 = Ett+k + 1 1 Ett+k-1 + 1 1 Et-1t+k + 1 2 1 Et-1t+k-1 = (k) + 1 1 (k - 1) + 1 1 (k + 1) 1 2 1 (k) = 1 + 1 2 1 (k) + 1 1 (k - 1) + 1 1 (k + 1) = 2 (1 + 2 1) k = 0 12 k = 1 0 jinak I když obě invertibilní i neinvertibilní MA reprezentace generují procesy se stejnými momenty prvního a druhého řádu, z praktických důvodů dáváme přednost procesu invertibilnímu, neboť nepozorovatelné veličiny t můžeme odhadnout pomocí přítomných a minulých hodnot pozo- rovatelných veličin {Xs, s < t}, kdežto u neinvertibilních MA reprezentací nepozorovatelné veličiny t neodhadneme, neboť nemáme ještě k dispozici budoucí hodnoty {Ys, s > t}. Nyní již můžeme podat definici invertibility. Definice 1.2.12. ARMA proces Yt ARMA(p, q) se nazývá invertibilní, jestliže existuje absolutně kon- vergentní posloupnost reálných čísel = {j} j=0 (tj. j=0 |j| < ,) tak, že t = j=0 jYt-j, tj. zkráceně Yt AR() : t = (B)Yt. Nyní vyšetřeme, za jakých podmínek je MA proces řádu q invertibilní. MA proces řádu q: Yt MA(q) : Yt = t + 1t-1 + + qt-q t WN(0, 2 ) . Naprosto analogickým postupem jako v případě kauzálního AR(p) procesu, lze ukázat, že všechny kořeny (z) musí ležet vně jednotkového kruhu. Proveďme tedy nejprve rozklad polynomu (z) = 1 + 1z + + qzq na součin kořenových činitelů (z) = (1 - 1z)q1 . . . (1 - kz)qk , kde z01 = 1 1 , . . . , z0k = 1 k jsou rozdílné (reálné či komplexní) kořeny polynomu (z), q1, . . . , qk je jejich násobnost (přičemž platí q1 + + qk = q). Nyní budeme hledat taková absolutně konvergentní = {j} j=0 (tj. j=0 |j| < ), aby t = j=0 jYt-j byl invertibilní proces. 1.2. Základní pojmy 18 Pokud použijeme operátor zpětného chodu, můžeme psát: Yt = (B)t, přitom hledáme (B) takové, aby platilo (B)(B) = 1 nebo (z)(z) = 1 čili (z) = 1 (z) . Z věty o rozkladu na částečné zlomky dostáváme (pokud pro názornost předpokládáme, že všechny kořeny jsou jednoduché) 1 (z) = 1 (1 - 1z) . . . (1 - pz) = c1 1 - 1z + + cp 1 - pz pro vhodná c1, . . . , cp. Pokud pro k = 1, . . . , p platí |kz| < 1 , můžeme psát ck 1 - kz = ck j=0 (kz)j a dokázali jsme najít konvergentní řadu (z) = j=0 jzj = p k=1 ck j=0 j kzj = j=0 c1j 1 + + cpj p zj , přičemž j = c1j 1 + + cpj p, neboť (z) = 1 (z) je holomorfní pro |z| 1 právě když |1| < 1, . . . , |p| < 1 1 1 > 1 |z01|>1 , . . . , 1 p > 1 |z0p|>1 , tedy všechny kořeny polynomu (z) musí ležet vně jednotkového kruhu. Na závěr tohoto odstavce ještě spočítejme střední hodnotu, rozptyl, autokovarianční funkci a také spektrální hustotu MA(q) procesu. Protože MA(q) proces je lineárním procesem, je vždy slabě stacionární, proto můžeme počítat EYt = E(t + 1t-1 + + qt-q) = 0 DYt = D(t + 1t-1 + + qt-q) = 2 (1 + 2 1 + + 2 q ) (t) = C(Ys, Ys+t) = EYsYs+t = q j=0 q h=0 jhEs-js+t-h = q j=0 q h=0 jh(h - t - j) = 2 q j=0 jj+t Protože 0 = 1 a j = 0 pro j > q, 1.2. Základní pojmy 19 dostáváme (t) = 2 (1 + 2 1 + + 2 q-2 + 2 q-1 + 2 q ) pro t = 0 2 (1 + 12 + + q-2q-1 + q-1q) t = 1 2 (2 + 13 + + q-2q) t = 2 ... ... 2 (q-1 + 1q) t = q - 1 2 q t = q 0 jinak = 2 q-t j=0 jj+t t = 0, . . . , q, 0 jinak. . Autokorelační funkce je pak rovna (t) = 1 t = 0 1 1+2 1++2 q q-t j=0 jj+t 1 t q, 0 1 0 jinak tedy, pro t > q je autokorelační funkce nulová, což je velmi důležitá identifikační vlastnost MA(q) procesů. Díky tomu, že MA(q) proces je lineárním procesem, spektrální hustota je rovna fMA(q)() = 2 2 e-i 2 . 1.2.8 Vícenásobná reprezentace MA(q) procesů Mějme MA proces řádu q: Yt MA(q) : Yt = t + 1t-1 + + qt-q t WN(0, 2 ). Proveďme tedy rozklad polynomu (z) = 1 + 1z + + qzq na součin kořenových činitelů (z) = j (1 - jz), Pak (protože MA(q) proces je lineárním procesem) autokovarianční generující funkce je rovna GY (z) = (z) z-1 2 . Dále platí (1 - jz)(1 - jz-1 ) = 1 - jz - jz-1 + 2 j = 2 j (-2 j - -1 j z - -1 j z-1 + 1) = 2 j 1 - 1 j z 1 - 1 j z-1 1.3. Vlastnosti ARMA(p, q) procesů 20 Tudíž GY (z) = 2 (z) z-1 = 2 j (1 - jz) j (1 - jz-1 ) = 2 j 2 j 2 j 1 - 1 j z (z) j 1 - 1 j z-1 (z-1) = 2 (z)(z-1 ) Takže proces Y t MA(q) : Yt = t + 1 t-1 + + q t-q t WN(0, 2 ) má stejnou autokovarianční generující funkcí GY (z) = 2 (z)(z-1 ) a jsou proto z hlediska prvních dvou momentů nerozlišitelné. Obecně můžeme dostat 2q různých procesů s funkcí s(z) = q j=1 (1 - 1 j z) s = 1, . . . , 2q Mezi všemi těmito procesy pouze jediný je invertibilní, a to ten, pro kterého platí invert j = j |j| < 1, -1 j |j| 1. Takže podmínka invertibility zajišťuje identifikovatelnost MA(q) procesu z hlediska prvních dvou momentů. 1.3 Vlastnosti ARMA(p, q) procesů Dříve než uvedeme nutnou a postačující podmínku pro kauzalitu a invertibilitu ARMA(p, q) procesů, vyšetřeme problematiku společných kořenů (z) a (z). 1.3.1 Společné kořeny polynomů (z) a (z) Mějme Yt ARMA(p, q) : Yt -1Yt-1 - -pYt-p = t +1t-1 + + qt-q t WN(0, 2 ) a předpokládejme že (z) a (z) mají společný kořen 1 . Pak můžeme psát (z) = (1 - z)(1 - 1z - - p-1zp-1 ) = (1 - z) (z) (z) = (1 - z)(1 + 1z + + q-1zq-1 ) = (1 - z) (z) Pak můžeme psát (1 - B) (B)Yt = (1 - B) (B)t. Pokud obě strany rovnice vydělíme výrazem (1 - B), dostaneme Yt ARMA(p - 1, q - 1) : (B)Yt = (B)t. Takže podmínka, že (z) a (z) nemají společné kořeny zajišťuje, že řády ARMA procesů nelze již snižovat. 1.3. Vlastnosti ARMA(p, q) procesů 21 1.3.2 Nutná a postačující podmínka kauzality a invertibility. V předchozích odstavcích jsme ukázali, že platí (1) Yt AR(p) : (B)Yt = t (z) = 0 pro z C |z| 1 AR(p) je kauzální; (2) Yt MA(q) : Yt = (B)t (z) = 0 pro z C |z| 1 MA(q) je invertibilní. Naprosto analogickým způsobem lze dokázat obecnější tvrzení: Věta 1.3.1. Nechť (B) a (B) nemají společné kořeny. Pak Yt ARMA(p, q) : (B)Yt = (B)t je kauzální (z) = 0 pro z C |z| 1. Yt ARMA(p, q) : (B)Yt = (B)t je invertibilní (z) = 0 pro z C |z| 1. Znamená to tedy, že Yt ARMA(p, q) je kauzálním a invertibilním ARMA procesem, jestliže všechny kořeny polynomů (z) a (z) leží vně jednotkového kruhu a koeficienty j a j jsou určeny ze vztahů (z) = j=0 jzj = (z) (z) pro |z| 1 (z) = j=0 jBj = (z) (z) pro |z| 1. V dalším budeme uvažovat pouze takové Yt ARMA(p, q) : (B)Yt = (B)t procesy, které splňují následující podmínky (P1) (B) a (B) nemají společné kořeny. (P2) Yt ARMA(p, q) je kauzální. (P3) Yt ARMA(p, q) je invertibilní. 1.3.3 Střední hodnota, rozptyl, autokovarianční a autokorelační funkce procesů ARMA(p, q) Střední hodnota Vzhledem ke kauzalitě ARMA(p, q) procesu můžeme počítat EYt = E j=0 jt-j = j=0 jEt-j = 0 Rozptyl Při odvození rozptylu nejprve rovnici Yt = 1Yt-1 + + pYt-p + t + 1t-1 + + qt-q 1.3. Vlastnosti ARMA(p, q) procesů 22 vynásobme výrazem Yt a spočtěme střední hodnoty obou stran, tj. EY 2 t = 1EYt-1Yt + + pEYt-pYt + EtYt + 1Et-1Yt + + qEt-qYt. (A6) Spočtěme pro i = 0, 1, . . . , q Et-iYt = Et-i j=0 jt-j = j=0 jEt-it-j = j=0 j(i-j) = 2 i (přičemž 0 = 1). Po dosazení do rovnice (A6) dostaneme (0) - 1(1) - - p(p) = 2 (1 + 11 + . . . + qq). (A7) Podělme obě strany rovnice (A7) výrazem (0). Vzhledem k tomu, že (k) = (k) (0) , dostaneme 1 - 1(1) - - p(p) = 2 (1 + 11 + . . . + qq) (0) takže DYt = (0) = 2 (1 + 11 + . . . + qq) 1 - 1(1) - - p(p) . Autokovarianční a autokorelační funkce (ACF) Při výpočtu autokovariance rovnici Yt - 1Yt-1 - - pYt-p = t + 1t-1 + + qt-q vynásobíme výrazem Yt-k a spočítáme střední hodnoty obou stran, takže dostaneme (k) - 1(k - 1) - - p(k - p) = EYt-kt + 1EYt-kt-1 + + qEYt-kt-q. (A8) Nejprve je třeba si uvědomit, že pro k > 0 platí EtYt-k = E t j=0 jt-k-j = j=0 jEtt-k-j = j=0 j(k + j) = 0 . Spočtěme pro i = 1, . . . , q Et-iYt-k = Et-i j=0 jt-k-j = j=0 jEt-it-j-k = j=0 j(k + j - i) = 2 i-k . Dosadíme-li do výrazu (A8) předchozí dva výsledky, dostáváme (k) - 1(k - 1) - - p(k - p) = 2 (11-k + + k k-k =0=1 +k+11 + + qq-k) 1.3. Vlastnosti ARMA(p, q) procesů 23 Vzhledem k tomu, že platí j = 0 pro j < 0, uvažujme dva případy: (a) jestliže 0 k max(p, q + 1) , pak (k) - 1(k - 1) - - p(k - p) = 2 (k + k+11 + + qq-k) . (A9) (b) pokud k > max(p, q + 1) , pak (k) - 1(k - 1) - - p(k - p) = 0 . (A10) Podělme obě strany rovnice (A9), resp (A10), výrazem (0). Dostaneme (k) - 1(k - 1) - - p(k - p) = 2 (k + k+11 + + qq-k) (0) resp. (k) - 1(k - 1) - - p(k - p) = 0. Nechť např. q + 1 > p. Pak máme více rovnic pro určení počátečních p podmínek. V tomto případě prvních q - p + 1 autokovariančních koeficientů jsou určeny z prvních q - p + 1 podmínek. Obecné řešení homogenní diferenční rovnice (k) - 1(k - 1) - - p(k - p) = 0 tj. (B)(k) = 0 je tvaru (k) = k j=1 pj-1 s=0 cjsks k j , kde 1 1 , . . . , 1 k jsou různé kořeny s násobnostmi p1, . . . , pk, přičemž p = p1 + + pk a cjs je právě p konstant, které jsou určeny počátečními podmínkami. Dále položme j = rjeij . Pak máme (k) = k j=1 pj-1 s=0 cjsks rk j eikj , Protože uvažujeme pouze kauzální ARMA(p, q) procesy, pro které platí |j| = rj < 1, dostáváme odtud, že a (k) klesá pro k exponenciálně k nule, tj. (k) --- k 0 a vidíme, že pro k > max(p, q + 1) má autokorelační funkce ARMA(p, q) procesů podobný charakter jako autoregresní AR(p) procesy. 1.3. Vlastnosti ARMA(p, q) procesů 24 1.3.4 Spektrální hustota ARMA(p, q) procesů Věta 1.3.2 (Spektrální hustota ARMA(p, q) procesů). Nechť (B)Yt = (B)t je kauzální a invertibilní ARMA(p, q) proces, přičemž (z) a (z) nemají společné kořeny. Pak spektrální hustota ARMA(p, q) procesu je rovna fY () = fARMA(p,q)() = 2 2 | (e-i )| 2 | (e-i)|2 pro -, . Důkaz. Protože všechny kořeny (z) a (z) leží vně jednotkového kruhu, ARMA(p, q) proces je kauzální a invertibilní. Kauzalita značí, že existuje absolutně konvergentní posloupnost reálných čísel = {j} j=0 (tj. j=0 |j| < ) taková, že platí Yt = j=0 jt-j kde t WN(0, 2 ). Víme, že spektrální hustota bílého šumu je rovna f() = 2 2 kde -, . Protože Yt je lineárním procesem, víme, že má spektrální hustotu fY () = e-i 2 f() = e-i 2 2 2 kde -, . Také (B)t jakožto lineární proces má spektrální hustotu tvaru e-i 2 2 2 pro -, . Rovněž (B)Yt jakožto lineární filtr má také spektrální hustotu, a ta je rovna e-i 2 fY () pro -, . Protože platí (B)Yt = (B)t, musí také platit e-i 2 fY () = e-i 2 2 2 pro -, . Odtud již dostáváme tvrzení věty fY () = fARMA(p,q)() = 2 2 | (e-i )| 2 | (e-i)|2 pro -, . 1.4. Stacionární procesy a nejlepší lineární predikce 25 1.4 Stacionární procesy a nejlepší lineární predikce Nechť {Yt, t Z} L2 (, A, P) je stacionární proces se střední hodnotou Y a autokovarianční funkcí Y (t). Pak náhodný proces {Yt-Y , t Z} má nulovou střední hodnotu (tj. je centrován) a má stejnou autokovarianční funkci Y (t). Uvažujme nejlepší lineární predikci ^Yt pomocí Yt-1, . . . , Yt-n (n 1), která je ortogonální projekcí ^Yt = Psp{1,Yt-1,...,Yt-n}(Yt). Lze snadno ukázat,že platí ^Yt = Psp{1,Yt-1,...,Yt-n}(Yt) = Y + Psp{Yt-1,...,Yt-n}(Yt). Takže bez újmy na obecnosti můžeme dále uvažovat pouze centrované stacionární procesy {Yt, t Z}, pro které platí ^Yt = Psp{1,Yt-1,...,Yt-n}(Yt) = Psp{Yt-1,...,Yt-n}(Yt). Nejprve definujme jednokrokovou predikci. Definice 1.4.1. Nechť {Yt, t Z} L2 (, A, P) je centrovaný stacionární proces. Označme pro n 1 Mn = sp{Y1, . . . , Yn}. Pak jednokroková (lineární) predikce je definována vztahem ^Yn+1 = ^Yn+1|n = 0 (= Y ) n = 0, Psp{Y1,...,Yn}(Yn+1) = PMn (Yn+1) n 1. Protože pro n 1 ^Yn+1 Mn, pak platí ^Yn+1 = n,1Yn + + n,nY1 a n,1, . . . , n,n minimalizují Yn+1 - ^Yn+1 2 = E|Yn+1 - ^Yn+1|2 . Podle projekční věty pro každé X L2 (, A, P) a pro každé Y Mn platí X - ^X, Y = X, Y - ^X, Y = 0 X, Y = ^X, Y což je EXY = E ^XY , takže jestliže pro j = 1, . . . , n položíme X = Yn+1 a Y = Yn+1-j, pak musí platit EYn+1Yn+1-j = E ^Yn+1Yn+1-j (j) = E Yn+1-j n i=1 n,iYn+1-i = n i=1 n,iEYn+1-iYn+1-j = n i=1 n,i(i - j) 1.4. Stacionární procesy a nejlepší lineární predikce 26 což lze maticově zapsat takto (1) (2) ... (n) = (0) (1) (n - 1) (1) (0) (n - 2) ... ... ... ... (n - 1) (n - 2) (0) n,1 n,2 ... n,n tj. n = n n. Projekční věta zaručuje existenci právě jednoho řešení ^Yn+1 Mn pro nějaké n Rn (kterých obecně může být více, jejich výsledkem je však pouze jediné ^Yn+1). Jestliže n je regulární, máme právě jediné n Rn a platí n = -1 n n . Následující věta dává postačující podmínku k tomu, aby pro každé n N byla n regulární maticí. Věta 1.4.2. Jestliže platí (0) > 0 a (h) --- h 0, pak kovarianční matice n = ((i - j))n i,j=1 je regulární pro každé n N. Důkaz. Tento důkaz se provádí sporem, viz Brockwel, Davis (1987), str. 160-161. 2 Důsledek 1.4.3. Označme Yn = (Yn, . . . , Y1) . Jestliže platí (0) > 0 a (h) --- h 0, pak nejlepší lineární predikce ^Yn+1 náhodné veličiny Yn+1 je tvaru ^Yn+1 = n,1Yn + + n,nY1 tj. ^Yn+1 = nYn přičemž n = -1 n n . Střední kvadratická chyba je rovna vn = MSE( ^Yn+1) = E(Yn+1 - ^Yn+1)2 = (0) - n-1 n n . (1.2) Důkaz. Tvrzení týkající se tvaru nejlepší lineární predikce a vektoru n plynou z předchozích poznámek a předešlé věty. Zbývá vypočítat střední kvadratickou chybu. E(Yn+1 - ^Yn+1)2 = E(Yn+1 - nYn)2 = EY 2 n+1 - 2E (nYnYn+1) + E (nYn) 2 . 1.4. Stacionární procesy a nejlepší lineární predikce 27 Nejprve počítejme EYnYn+1 = (EYnYn+1, EYn-1Yn+1, . . . , EY1Yn+1) = ((1), (2), . . . , (n)) = n. Dále si všimněme, že lze psát (nYn) 2 = nYnYnn a počítejme EYnYn = E Yn Yn-1 ... Y1 (Yn, . . . , Y1) = EY 2 n EYnYn-1 EYnY1 EYn-1Yn EY 2 n-1 EYn-1Y1 ... ... ... ... EY1Yn EY1Yn-1 EY 2 1 = (0) (1) (n - 1) (1) (0) (n - 2) ... ... ... ... (n - 1) (n - 2) (0) = n Takže můžeme pokračovat ve výpočtu střední kvadratické chyby E(Yn+1 - ^Yn+1)2 = EY 2 n+1 - 2nEYnYn+1 + nEYnYnn = (0) - 2nn + nnn = (0) - 2n-1 n n + n-1 n n-1 n n = (0) - n-1 n n. 2 Nyní definujme h-krokovou (lineární) predikci. Definice 1.4.4. Nechť {Yt, t Z} L2 (, A, P) je centrovaný stacionární proces. Označme pro n 1 Mn = sp{Y1, . . . , Yn}. Pak h-kroková (lineární) predikce je definována vztahem ^Yn+h = ^Yn+h|n = 0 (= Y ) n, h = 0, Psp{Y1,...,Yn}(Yn+h) = PMn (Yn+h) n, h 1. Obdobným způsobem jako u jednokrokové predikce můžeme odvodit, že jestliže platí (0) > 0 a (h) --- h 0, pak nejlepší lineární h-kroková predikce ^Yn+h náhodné veličiny Yn+h je tvaru ^Yn+h = (h) n,1Yn + + (h) n,nY1 tj. ^Yn+h = (h) n Yn přičemž (h) n = -1 n (h) n a (h) n = ((h), (h + 1), . . . , (h + n - 1)) . Střední kvadratická chyba je rovna v(h) n = MSE( ^Yn+h) = E(Yn+1 - ^Yn+1)2 = (0) - (h) n -1 n (h) n . 1.5. Rekurentní metody pro výpočet nejlepší lineární predikce 28 1.5 Rekurentní metody pro výpočet nejlepší lineární pre- dikce Věta 1.5.1 (Durbin-Levinsův algoritmus). Nechť {Yt, t Z} L2 (, A, P) je centrovaný stacionární proces s autokovarianční funkcí (h) takovou, že (0) > 0 a (h) --- h 0. Jestliže ^Yn+1 = Psp{Y1,...,Yn}(Yn+1) = n,1Yn + + n,nY1 je nejlepší lineární predikce, pak pro koeficienty n = (n,1, . . . , n,n) a pro střední kvadratické chyby vn = E Yn+1 - ^Yn+1 2 platí následující vztahy 1,1 = (1) (0) = (1) v0 = (0) (1.3) n,n = (n) - n-1n-1 /vn-1 (1.4) (1) n = n-1 - n,n n-1 vn = vn-1 1 - 2 n,n (1.5) kde n-1 = (n-1,1, . . . , n-1,n-1) n-1 = (n-1,n-1, . . . , n-1,1) n = (n,1, . . . , n,n-1, n,n) (1) n = (n,1, . . . , n,n-1) Důkaz. Mějme Mn = sp{Y1, . . . , Yn}. Označme Mn-1 = sp{Y2, . . . , Yn} a M n-1 = sp{Y1 - PMn-1 (Y1)}. Vidíme, že M n-1 je ortogonální komplement Mn-1 v Mn. Takže pro libovolné Y L2 (, A, P), tedy i pro Yn+1 L2 (, A, P) musí platit ^Yn+1 = PMn (Yn+1) = PMn-1 (Yn+1) + PM n-1 (Yn+1). Označme e1 = Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) tj. e1 = 1. 1.5. Rekurentní metody pro výpočet nejlepší lineární predikce 29 Pak platí ^Yn+1 = PMn-1 (Yn+1) + Fourier.koef. Yn+1, e1 e1 = PMn-1 (Yn+1) + Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) = PMn-1 (Yn+1) + Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) 2 Y1 - PMn-1 (Y1) (1.6) Označme a = Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) 2 . Počítejme predikci PMn-1 (Y1) = n-1,1Y2 + + n-1,n-1Yn a označme n-1 = (n-1,1, . . . , n-1,n-1) , pak složky vektoru n-1 jsou řešením rovnic EY1Y1+k (k) = E Y1+kPMn-1 (Y1) = n-1 j=1 n-1,j EY1+jY1+k (j-k) k = 1, . . . , n - 1, což lze zapsat maticově (1) (2) ... (n - 1) = (0) (1) (n - 2) (1) (0) (n - 3) ... ... ... ... (n - 2) (n - 3) (0) n-1,1 n-1,2 ... n-1,n-1 tj. n-1 = n-1 n-1, Jestliže platí (0) > 0 a (h) --- h 0, pak kovarianční matice n-1 je regulární a platí n-1 = -1 n-1n-1. Obdobně PMn-1 (Yn+1) = n-1,1Yn + + n-1,n-1Y2 a označíme-li n-1 = (n-1,1, . . . , n-1,n-1) , pak složky vektoru n-1 jsou řešením rovnic EYn+1Yn+1-k (k) = E Yn+1-kPMn-1 (Yn+1) = n-1 j=1 n-1,j EYn+1-jYn+1-k (j-k) k = 1, . . . , n - 1, 1.5. Rekurentní metody pro výpočet nejlepší lineární predikce 30 což maticově zapíšeme (1) (2) ... (n - 1) = (0) (1) (n - 2) (1) (0) (n - 3) ... ... ... ... (n - 2) (n - 3) (0) n-1,1 n-1,2 ... n-1,n-1 tj. n-1 = n-1 n-1, tedy n-1 = -1 n-1n-1 = n-1. Celkově tedy, označíme-li Y n-1 = (Y2, . . . , Yn) Yn-1 = (Yn, . . . , Y2) tak dostaneme PMn-1 (Y1) = n-1Y n-1 PMn-1 (Yn+1) = n-1Yn-1 Využijme vzorce (1.2) z důsledku 1.4.3 a počítejme střední kvadratické chyby obou predikcí Y1 - PMn-1 (Y1) 2 = (0) - n-1-1 n-1n-1 Yn+1 - PMn-1 (Yn+1) 2 = (0) - n-1-1 n-1n-1 Tedy vn-1 = Y1 - PMn-1 (Y1) 2 = Yn+1 - PMn-1 (Yn+1) 2 . Vraťme se k rovnici (1.6), pak ^Yn+1 = PMn (Yn+1) = PMn-1 (Yn+1) + a Y1 - PMn-1 (Y1) = aY1 + n-1,1Yn + + n-1,n-1Y2 - a (n-1,1Y2 + + n-1,n-1Yn) = aY1 + n-1 j=1 (n-1,j - an-1,n-j) Yn+1-j = n,nY1 + n,n-1Y2 + + n,1Yn a odtud dostaneme n,n = a n,j = n-1,j - an-1,n-j pro j = 1, . . . , n - 1. Nyní se vrátíme ke konstantě a. Protože Yn+1, Y1 - PMn-1 (Y1) = Yn+1, Y1 - Yn+1, PMn-1 (Y1) = EY1Yn+1 - E Yn+1 n-1 j=1 n-1,jY1+j = (n) - n-1 j=1 n-1,jEY1+jYn+1 = (n) - n-1 j=1 n-1,j(n - 1) = (n) - n-1n-1, 1.5. Rekurentní metody pro výpočet nejlepší lineární predikce 31 pak a = n,n = Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) 2 = (n) - n-1 j=1 n-1,j(n - 1) vn-1 = v-1 n-1 (n) - n-1n-1 . Zbývá dokázat, že platí vn = vn-1(1 - 2 n,n). Počítejme proto střední kvadratickou chybu predikce PMn (Yn+1): vn = Yn+1 - PMn (Yn+1) 2 = Yn+1 - PMn-1 (Yn+1) -PM n-1 (Yn+1) 2 = Yn+1 - PMn-1 (Yn+1) -PM n-1 (Yn+1), Yn+1 - PMn-1 (Yn+1) -PM n-1 (Yn+1) = Yn+1 - PMn-1 (Yn+1) 2 - 2 Yn+1 - PMn-1 (Yn+1), PM n-1 (Yn+1) + PM n-1 (Yn+1) 2 Protože PM n-1 (Yn+1) = a Y1 - PMn-1 (Y1) , tedy PM n-1 (Yn+1) 2 = a2 Y1 - PMn-1 (Y1) 2 = a2 vn-1. Zbývá dopočítat Yn+1 - PMn-1 (Yn+1), PM n-1 (Yn+1) = Yn+1 - PMn-1 (Yn+1), a Y1 - PMn-1 (Y1) = a Yn+1, Y1 - PMn-1 (Y1) - a PMn-1 (Yn+1), Y1 - PMn-1 (Y1) =0(ortogonal) Dále víme, že a = Yn+1, Y1 - PMn-1 (Y1) vn-1 = n,n, tedy a Yn+1, Y1 - PMn-1 (Y1) = a2 vn-1 a celkově vn = vn-1 - 2a2 vn-1 + a2 vn-1 = vn-1(1 - a2 ) = vn-1(1 - 2 n,n) . Tím jsou všechna tvrzení věty dokázána. 2 1.5. Rekurentní metody pro výpočet nejlepší lineární predikce 32 Důsledek 1.5.2 (Důsledek Durbin-Levinsonova algoritmu). Nechť {Yt, t Z} L2 (, A, P) je centrovaný stacionární proces s autokovarianční funkcí (h), pro kterou platí (0) > 0 a (h) --- h 0. Označme Mn = sp{Y1, . . . , Yn} Mn-1 = sp{Y2, . . . , Yn} a nejlepší lineární predikci ^Yn+1 = PMn (Yn+1) = n,nY1 + n,n-1Y2 + + n,1Yn, pak platí n,n = R Yn+1 - PMn-1 (Yn+1), Y1 - PMn-1 (Y1) . (1.7) Důkaz. Označme M n-1 = sp{Y1 - PMn-1 (Y1)}. Tedy M n-1 je ortogonální komplement Mn-1 v Mn. Takže lze psát ^Yn+1 = PMn (Yn+1) = PMn-1 (Yn+1) + PM n-1 (Yn+1). Protože PMn-1 (Yn+1) a Y1 - PMn-1 (Y1) jsou ortogonální tj. PMn-1 (Yn+1), Y1 - PMn-1 (Y1) = 0. Tedy Yn+1, Y1 - PMn-1 (Y1) = Yn+1 - PMn-1 (Yn+1), Y1 - PMn-1 (Y1) . Navíc platí (viz důkaz Durbin-Levinsonova algoritmu) Y1 - PMn-1 (Y1) 2 = Yn+1 - PMn-1 (Yn+1) 2 = vn-1, takže lze psát n,n = a = Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) 2 = Yn+1 - PMn-1 (Yn+1), Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) Yn+1 - PMn-1 (Yn+1) = E Yn+1 - PMn-1 (Yn+1) Y1 - PMn-1 (Y1) E Yn+1 - PMn-1 (Yn+1) 2 E Y1 - PMn-1 (Y1) 2 = C Yn+1 - PMn-1 (Yn+1), Y1 - PMn-1 (Y1) D Yn+1 - PMn-1 (Yn+1) D Y1 - PMn-1 (Y1) = R Yn+1 - PMn-1 (Yn+1), Y1 - PMn-1 (Y1) . 2 1.6. Parciální autokorelační funkce (PACF) 33 1.6 Parciální autokorelační funkce (PACF) Definice 1.6.1. Nechť {Yt, t Z} L2 (, A, P) je stacionární proces. Pak parciální autokorelační funkce je definována vztahem (1) = R(Yt, Yt+1) (k) = R(Yt - ^Yt, Yt-k - ^Yt-k) pro |k| > 1 kde ^Yt, resp. ^Yt-k jsou nejlepší lineární predikce Yt (resp. Yt-k) pomocí Yt-k+1, . . . , Yt-1. Nejlepší lineární predikce ^Yt a ^Yt-k jsou projekce ^Yt = PMk-1 (Yt) a ^Yt-k = PMk-1 (Yt-k) , kde Mk-1 = sp{Yt-k+1, . . . , Yt-1}. Přitom existují taková k-1 = (k-1,1, . . . , k-1,k-1) , že platí ^Yt = k-1,1Yt-1 + + k-1,k-1Yt-k+1 a také taková k-1 = (k-1,1, . . . , k-1,k-1) , že platí ^Yt-k = k-1,1Yt-k+1 + + k-1,k-1Yt-1, která minimalizují E(Yt - ^Yt)2 resp. E(Yt-k - ^Yt-k)2 , přičemž (jak již víme z důkazu Durbin-Levinsonova algoritmu) platí k-1,1 = k-1,1, . . . , k-1,k-1 = k-1,k-1 tj. k-1 = k-1. Celkově tedy, označíme-li Y k-1 = (Yt-k+1, . . . , Yt-1) Yk-1 = (Yt-1, . . . , Yt-k+1) tak dostaneme PMn-1 (Yt-k) = k-1Y k-1 PMn-1 (Yt) = k-1Yk-1 Víme, že pokud pro autokovarianční funkci (h) platí (0) > 0 a (h) --- h 0, pak matice k-1 je regulární a neznámé složky vektoru k-1 jsou rovny k-1 = -1 k-1k-1. Avšak podle důsledku 1.5.2 Durbin-Levinsonova algoritmu není třeba počítat inverzní matici -1 k-1, odtud k-1, následně ^Yt-k = PMn-1 (Yt-k) a ^Yt = PMn-1 (Yt) a nakonec korelační koeficient (k) = R(Yt - ^Yt, Yt-k - ^Yt-k), neboť platí (k) = k,k = R Yt - PMk-1 (Yt), Yt-k - PMk-1 (Yt-k) . 1.6. Parciální autokorelační funkce (PACF) 34 Shrneme-li předchozí, pak podle definice Forward ---- Backward ----- Yt-k Mk-1=sp{Yt-k+1,...,Yt-1} Yt-k+1, . . . , Yt-1 Yt r (F ) t-k = Yt-k - PMk-1 (Yt-k) r (B) t = Yt - PMk-1 (Yt) (k) = R r (F ) t-k, r (B) t a podle důsledku Durbin-Levinsonova algoritmu Mk=sp{Yt-k,Yt-k+1,...,Yt-1} Yt-k , Yt-k+1, . . . , Yt-1 Yt (k) = k,k , kde ^Yt = PMk (Yt) = n,1Yt-1 + + k,k =(k) Yt-k 1.7. Inovační algoritmus 35 1.7 Inovační algoritmus Základní myšlenkou Durbin-Levinsonova algoritmu je rozdělení Mn = sp{Yn, . . . , Y1} na dva ortogonální podprostory Mn-1 = sp{Yn, . . . , Y2} a M n-1 = sp{Y1 - PMn-1 (Y1)}. Následující rekurentní algoritmus spočívá v dekompozici Mn na n ortogonálních Hilbertových podprostorů pomocí Gram-Schmidtova algoritmu. Rekurentní algoritmus lze aplikovat nejen na stacionární procesy, ale obecně na procesy s konečnými druhými momenty. Pro jednoduchost předpokládejme, že jsou centrované. Nejprve zaveďme následujicí značení: (i, j) = EXiXj. Stejně označme Mn = sp{Yn, . . . , Y1} vn = Yn+1 - ^Yn+1 2 . Pokud označíme ^Yn = 0 (= Y ) pro n = 1 PMn-1 (Yn) pro n = 2, 3, . . . pak zřejmě Mn = sp{Yn - ^Yn, . . . , Y1 - ^Y1} n 1. Definujme tzv. inovaci vztahem Un+1 = Yn+1 - ^Yn+1 = Yn+1 - n j=1 n,jYn+1-j. Označme Un = (U1, . . . , Un) Yn = (Y1, . . . , Yn) ^Yn = ( ^Y1, . . . , ^Yn) . Pak lze psát Un = AnYn, kde matice An je dolní trojúhelníkovou maticí An = 1 0 0 -1,1 1 0 0 -2,2 -2,1 1 0 0 ... ... ... ... ... ... -n-2,n-2 -n-2,n-3 -n-2,1 1 0 -n-1,n-1 -n-1,n-2 -n-1,1 1 . 1.7. Inovační algoritmus 36 Všimněme si, že determinant matice je roven 1, takže existuje inverzní matice Cn = A-1 n , která je také dolní trojúhelníkovou maticí. Upravujme postupně ^Yn = Yn - Un = A-1 n Un - Un = A-1 n Un - In Un = nUn, kde n = Cn - In = 0 0 0 1,1 0 0 0 2,2 2,1 0 0 0 ... ... ... ... ... ... n-2,n-2 n-2,n-3 n-2,1 0 0 n-1,n-1 n-1,n-2 n-1,1 0 . Protože ^Yn = nUn = n(Yn - ^Yn) a protože n je dolní trojúhelníkovou maticí, můžeme psát ^Yn+1 = 0 n = 0 n j=1 n,j(Yn+1-j - ^Yn+1-j) n = 1, 2, . . . . Věta 1.7.1 (Inovační algoritmus). Nechť {Yt, t Z} je centrovaný náhodný proces s konečnými druhými momenty, přičemž kovarianční matice (EYiYj)n i,j=1 = ((i, j))n i,j=1 je regulární pro každé n N. Pak pro jed- nokrokovou predikci platí následující rekurentní vztahy ^Yn+1 = 0 n = 0 n j=1 n,j Yn+1-j - ^Yn+1-j n = 1, 2, . . . (1.8) v0 = (1, 1) (1.9) n,n-k = v-1 k (n + 1, k + 1) - k-1 j=0 k,k-jn,n-jvj k = 0, . . . , n - 1 (1.10) vn = (n + 1, n + 1) - n-1 j=0 2 n,n-jvj (1.11) Důkaz. Tvrzení (1.8) jsem dokázali již v předchozím textu. Pro n 1 proveďme následující přeindexování: ^Yn+1 = n j=1 n,j Yn+1-j - ^Yn+1-j = n,1 Yn - ^Yn + + n,n Y1 - ^Y1 = n-1 k=0 n,n-k Yk+1 - ^Yk+1 1.7. Inovační algoritmus 37 Pro k = 0, 1, . . . , n - 1 vynásobme obě strany předchozí rovnice výrazem Yk+1 - ^Yk+1 a vypočítejme střední hodnoty obou stran rovnic. Dostaneme: E ^Yn+1 Yk+1 - ^Yk+1 = n-1 j=0 n,n-jE Yj+1 - ^Yj+1 Yk+1 - ^Yk+1 nebo ekvivalentně pomocí skalárních součinů: ^Yn+1, Yk+1 - ^Yk+1 = n-1 j=0 n,n-j Yj+1 - ^Yj+1, Yk+1 - ^Yk+1 = n,n-k Yk+1 - ^Yk+1, Yk+1 - ^Yk+1 = n,n-k Yk+1 - ^Yk+1 2 = n,n-kvk, s využitím vztahů Yj+1 - ^Yj+1, Yk+1 - ^Yk+1 = 0 pro j = k. Dále díky tomu, že pro n > k Yn+1 - ^Yn+1, Yk+1 - ^Yk+1 = 0, dostáváme ^Yn+1, Yk+1 - ^Yk+1 = Yn+1, Yk+1 - ^Yk+1 = n,n-kvk. Dále upravujme n,n-kvk = Yn+1, Yk+1 - ^Yk+1 = Yn+1, Yk+1 =(n+1,k+1) - Yn+1, ^Yk+1 = (n + 1, k + 1) - Yn+1, k-1 j=0 k,k-j Yj+1 - ^Yj+1 = (n + 1, k + 1) - k-1 j=0 k,k-j Yn+1, Yj+1 - ^Yj+1 =n,n-j vj = (n + 1, k + 1) - k-1 j=0 k,k-jn,n-jvj Odtud jednoduchou úpravou dostaneme tvrzení (1.10): n,n-k = v-1 k (n + 1, k + 1) - k-1 j=0 k,k-jn,n-jvj . Nakonec díky tomu, že ^Yn+1, Yn+1 - ^Yn+1 = 0, 1.7. Inovační algoritmus 38 pak s využitím Pythagorovy věty dostaneme Yn+1 2 =(n+1,n+1) = Yn+1 - ^Yn+1 + ^Yn+1 2 = Yn+1 - ^Yn+1 2 =vn + ^Yn+1 2 . Tedy vn = (n + 1, n + 1) - ^Yn+1, ^Yn+1 = (n + 1, n + 1) - ^Yn+1, n-1 j=0 n,n-j Yj+1 - ^Yj+1 = (n + 1, n + 1) - n-1 j=0 n,n-j ^Yn+1, Yj+1 - ^Yj+1 =n,n-j vj = (n + 1, n + 1) - n-1 j=0 2 n,n-jvj čímž jsme dokázali poslední tvrzení (1.11). 2 Poznámka 1.7.2. Zatímco Durbin-Levinsův algoritmus dává koeficienty n,j v reprezentaci ^Yn+1 = n j=1 n,jYn+1-j = n-1 j=0 n,n-jYj+1, inovační algoritmus dává koeficienty n,j v ortogonálním rozvoji ^Yn+1 = n j=1 n,j Yn+1-j - ^Yn+1-j = n-1 j=0 n,n-j Yj+1 - ^Yj+1 . Poznámka 1.7.3. Inovační algoritmus dává ,,inovační reprezentaci samotných Yn+1, neboť platí Yn = Yn - ^Yn + ^Yn = (Yn - ^Yn) + (Cn - In)(Yn - ^Yn) = Cn(Yn - ^Yn) . a položíme-li n,0 = 1, můžeme psát Yn+1 = n j=0 n,j Yn+1-j - ^Yn+1-j = n j=0 n,n-j Yj+1 - ^Yj+1 . Tyto vztahy využijeme později při odvozování maximálně věrohodných odhadů neznámých parametrů n,j. 1.8. Výstavba modelů v B-J metodologii 39 1.8 Výstavba modelů v B-J metodologii 1.8.1 Odhady v ARMA procesech Určení vhodného ARMA(p, q) modelu pro danou realizaci stacionárního procesu {Yt, t Z} v sobě zahrnuje celou řadu problémů 1. výběr řádu modelu p a q, tj. provést identifikaci modelu; 2. odhad parametrů 1, . . . , p, 1, . . . , q a 2 ; 3. ověření vhodnosti modelu. Identifikace modelu je první fází výstavby modelu a jejím úkolem je rozhodnout, jaký typ modelu vybrat (tj. zda AR, MA nebo ARMA) a explicitně určit řád modelu. Před vlastní identifikací se doporučuje provést některé z následujících přípravných operací: 1. Pořídit grafický záznam řady. 2. Pokud střední hodnota je nenulová, provést odhad střední hodnoty a následně provést centrování. VLASTNÍ IDENTIFIKACE PROCESU Identifikace procesu se provádí na základě zkoumání průběhu ˇ odhadnuté autokorelační funkce ACF r(k) = ^(k); ˇ a parciální autokorelační funkce PACF a(k) = ^(k). Snažíme se především zjistit existenci případného identifikačního bodu k0. Je však nutné mít na paměti, že pracujeme pouze s odhadnutými hodnotami r(k) a a(k), takže naše závěry mohou být někdy dost zkreslené. Doporučuje se proto netrvat na jednoznačné identifikaci určitého modelu, ale přijmout a přezkoušet několik alternativ. AR(p) MA(q) ARMA(p, q) ACF neexistuje ko ko = q neexistuje ko (k) ve tvaru ,,U po q - p hodnotách je ve tvaru ,,U PACF ko = p neexistuje ko neexistuje ko (k) omezena křivkou ,,U po p - q hodnotách je omezena křivkou ,,U Identifikační bod ko je index, pro nějž platí: k > k0 je (k) = 0 (resp. (k) = 0). Křivka ,,U je křivka ve tvaru lineární kombinace klesajících geometrických posloupností a sinusoid s geometricky klesající amplitudou. 1.8. Výstavba modelů v B-J metodologii 40 1.8.2 Yuleovy-Walkerovy rovnice a odhad parametrů v AR(p) Nechť {Yt, t Z} je centrovaný kauzální autoregresní proces AR(p) : (B)Yt = t t WN(0, 2 ). Vraťme se k Yuleovým-Walkerovým rovnicím (0) =1 = 1(1) + + p(p) + 2 (0) 2 = (0) [1 - 1(1) - - p(p)] (k) = 1(k - 1) + + p(k - p) k 1 Označíme-li ^Rp = (^(i - j))p i,j=1 ^p = (^(1), . . . , ^(p)) p = (1, . . . , p) ^p = ( ^1, . . . , ^p) a v Yuleových-Walkerových rovnicích nahradíme (k) odpovídajícími odhady ^(k), pak (pokud ^(0) > 0) dostaneme tzv. Yuleovy-Walkerovy odhady: ^p = ^R-1 p ^p ^2 = ^(0) 1 - ^p ^R-1 p ^p . Věta 1.8.1. Nechť {Yt, t Z} je centrovaný kauzální autoregresní proces AR(p) : (B)Yt = t, kde t IID(0, 2 ) a ^p je Yuleovův-Walkerův odhad p = (1, . . . , p) , pak platí n ^p - p A Np O, 2 -1 p , kde p = ((i - j))p i,j=1 . Kromě toho platí ^2 P - 2 . Důkaz. viz Brockwell, Davis (1987), str. 255-257. 2 Z předchozích tvrzení plyne, že odhady získané řešením Yuleových-Walkerových rovnic jsou asymptoticky nestranné a lze pro ně konstruovat asymptotické intervaly spolehlivosti. V praktických situacích však skutečný řád p autoregresního procesu neznáme. V tom případě se využijí tvrzení následující věty. Věta 1.8.2. Nechť {Yt, t Z} je centrovaný kauzální autoregresní proces AR(p) : (B)Yt = t, kde t IID(0, 2 ) a ^m = ^m,1, . . . , ^m,m = ^R-1 m ^m, m > p, pak platí n ^m - m A Nm O, 2 -1 m , kde m = ((i - j))m i,j=1, m jsou koeficienty nejlepší lineární predikce ^Ym+1 = Psp{Ym,...,Y1}(Ym+1) = mYm, kde Ym = (Ym, . . . , Y1) , tj. m = R-1 m m, kde Rm = ((i - j))m ij=1 . Specielně pro m > p platí n ^m,m A N(0, 1) . Důkaz. viz Brockwell, Davis (1987), str. 255-257. 2 1.8. Výstavba modelů v B-J metodologii 41 1.8.3 Předběžné odhady v AR(p) a Durbin-Levinsův algoritmus. Předpokládejme, že máme k dispozici pozorování y1, . . . , yn centrované stacionární posloupnosti {Yt, t Z} AR(m) : (B)Yt = t t WN(0, 2 ). Za předpokladu, že ^(0) > 0, pak můžeme odhadnout neznámé parametry autoregresního modelu řádu m < n pomocí Yuleových-Walkerových rovnic. Odhadnutý AR(m) proces je tvaru Yt - ^m,1Yt-1 - - ^m,mYt-m = t t WN(0, ^vm), kde ^m = ^m,1, . . . , ^m,m = ^R-1 m ^m ^vm = ^(0) 1 - ^m ^R-1 m ^m . Jestliže ^(0) > 0, pak R1, R2, . . . nejsou singulární a můžeme využít Durbin-Levinsův algorit- mus pro postupné odhady autoregresních koeficientů ^1, ^2 a odhady variability bílého šumu ^v1, ^v2, . . . . Věta 1.8.3 (Durbin-Levinsův algoritmus). Jestliže ^(0) > 0, pak parametry ^m,1, . . . , ^m,m a ^vm autoregresního modelu Yt - ^m,1Yt-1 - - ^m,mYt-m = t t WN(0, ^vm), pro m = 1, . . . , n - 1 lze získat rekurzivně ze vztahů ^1,1 = ^(1) ^(0) = ^(1) v0 = ^(0) (1.12) ^m,m = ^(m) - ^m-1 ^m-1 /^vm-1 (1.13) ^ (1) m = ^m-1 - ^m,m ^ m-1 ^vm = ^vm-1 1 - ^2 m,m (1.14) kde ^m-1 = (^m-1,1, . . . , ^m-1,m-1) ^ m-1 = (^m-1,m-1, . . . , ^m-1,1) ^m = (^m,1, . . . , ^m,m-1, ^m,m) ^ (1) m = (^m,1, . . . , ^m,m-1) 1.8. Výstavba modelů v B-J metodologii 42 1.8.4 Předběžné odhady v MA(q) a inovační algoritmus. Jestliže chceme na základě pozorování y1, . . . , yn centrované stacionární posloupnosti provést odhad MA(m) (m = 1, 2, . . . , n - 1) ve tvaru Yt = t + ^m,1t-1 + + ^m,mt-m t WN(0, ^vm), můžeme využít inovační algoritmus. Věta 1.8.4 (Odhady parametrů MA procesů pomocí inovačního algoritmu). Jestliže ^(0) > 0, pak odhady parametrů MA procesů lze provést pomocí následujících rekurentních vztahů ^v0 = ^(0) ^m,m-k = ^v-1 k ^(m - k) - k-1 j=0 ^m,k-j ^m,m-j ^vj k = 0, . . . , m - 1 ^vm = ^(0) - m-1 j=0 ^2 m,m-j ^vj Označme ^m = ^m,1, . . . , ^m,m . Uvedeme nyní větu, která platí obecně pro ARMA(p, q) proces. Připomeňme, že pro MA(q) proces jsou 1 = 1, . . . , q = q a j = 0 pro j > q. Věta 1.8.5 (Asymptotické vlastnosti inovačních odhadů ^m.). Nechť {Yt, t Z} je kauzální a invertibilní ARMA pro- ces (B)Yt = (B)t, t IID(0, 2 ), E4 t < a (z) = j=0 jzj = (z) (z) , |z| 1. Pak pro libovolnou posloupnost kladných ce- lých čísel {m(n)} n=1 takovou, že m < n, m a m = o n 1 3 , když n , pro každé k platí n ^m,1 - 1, ^m,2 - 2, . . . , ^m,k - k, A Nk(0, A), kde A = (aij)k i,j=1 a aij = min(i,j) r=1 i-rj-r. Kromě toho platí ^vm P - 2 . Důkaz. viz Brockwell, Davis (1987), str. 239. 2 I když rekurentní odhady koeficientů MA procesů pomocí inovačního algoritmu jsou analo- gické jako rekurentní odhady koeficientů AR procesů pomocí Durbin-Levinsonova algoritmu, je 1.8. Výstavba modelů v B-J metodologii 43 mezi nimi přece jen jistý rozdíl. Pro odhady ^p = (^p,1, . . . , ^p,p) pomocí Durbin-Levinsonova algoritmu platí ^p P - p, avšak odhady ^q = (^q,1, . . . , ^q,q) nekonvergují podle pravděpodobnosti k q. Ke konvergenci podle pravděpodobnosti je třeba odhad (^m,1, . . . , ^m,q) , kde posloupnost {m(n)} n=1 splňuje podmínky předchozí věty. Výběr m (maximálně až do n 4 ) pro výběr pevné délky se volí tak, aby odhady (^m,1, . . . , ^m,q) se stabilizovaly. 1.8.5 Předběžné odhady v ARMA(p, q) procesu. Nechť {Yt, t Z} je kauzální a invertibilní ARMA proces {Yt, t Z} ARMA(p, q) : (B)Yt = (B)t t WN(0, 2 ). Z kauzality vyplývá, že existuje posloupnost {j} j=0 taková, že j=0 |j| < a platí Yt = j=0 jt-j, tj. pro |z| 1 dostáváme (z) = (z) (z) (z)(z) = (z). Koeficienty {j} j=0 se určí ze vztahu (1 - 1z - 2z2 - - pzp )(0 + 1z + 2z2 + ) = (1 + 1z + 2z2 + + qzq ) porovnáním koeficientů u mocnin proměnné z , tj. z0 : 0 = 1 0 = 1 z1 : 1 - 1 = 1 1 = 1 + 1 z2 : 2 - 11 - 2 = 2 2 = 2 + 11 + 2 z3 : 3 - 21 - 12 - 3 = 3 3 = 3 + 21 + 12 + 3 ... Obecně, položíme-li j = 0 pro j > q dostaneme 0 = 1 j = j + min(j,p) i=1 ij-i j = 1, 2, . . . Za předběžné odhady 1, 2, . . . , p+q použijeme inovační odhady ^m,1, . . . , ^m,p+q, jejichž asymptotické vlastnosti dává předchozí věta. Takže dostáváme ^2 = ^vm a ^m,j = j + min(j,p) i=1 i ^m,j-i j = 1, 2, . . . , p + q. Nejprve uvažujeme rovnice pro j = q+1, . . . , p+q, ve kterých j = 0, které maticově vyjádříme takto ^m,q+1 ^m,q+2 ... ^m,q+p-1 ^m,q+p = ^m,q ^m,q-1 ^m,q+1-p ^m,q+1 ^m,q ^m,q-1 ^m,q+2-p ... ... ... ... ... ^m,q+p-2 ^m,q+1 ^m,q ^m,q-1 ^m,q+p-1 ^m,q+1 ^m,q 1 2 ... p-1 p . 1.8. Výstavba modelů v B-J metodologii 44 Řešením těchto rovnic dostaneme odhady ^1, . . . , ^p. Nakonec získáme odhady ^1, . . . , ^q ze vztahů ^j = ^m,j - min(j,p) i=1 ^i ^m,j-i j = 1, . . . , q. Poznámka 1.8.6. Pro MA(q) platí ^j = ^m,j , neboť p = 0. 1.8.6 Maximálně věrohodné odhady. Předpokládejme, že {Yt, t Z} je Gaussovský proces s nulovou střední hodnotou a kovari- anční funkcí (i, j) = EXiXj. Označme Yn = (Y1, . . . , Yn) a n = ((i, j))n i,j=1 . Věrohodnostní funkce náhodného vektoru Yn je tvaru L(n) = (2)- n 2 |n|- 1 2 exp{-1 2 Yn-1 n Yn}. Dále označme Mn = sp{Yn, . . . , Y1} a ^Yn = 0 (= Y ) pro n = 1 PMn-1 (Yn) pro n = 2, 3, . . . pak zřejmě Mn = sp{Yn - ^Yn, . . . , Y1 - ^Y1} n 1. Pro nejlepší lineární predikce použijme inovační algoritmus, podle kterého ^Yn+1 = 0 n = 0 n j=1 n,j Yn+1-j - ^Yn+1-j n = 1, 2, . . . přičemž střední kvadratickou chybu označme vn = Yn+1 - ^Yn+1 2 . Označíme-li ^Yn = ( ^Y1, . . . , ^Yn) a Cn = 1 0 0 1,1 1 0 0 2,2 2,1 1 0 0 ... ... ... ... ... ... n-2,n-2 n-2,n-3 n-2,1 1 0 n-1,n-1 n-1,n-2 n-1,1 1 , 1.8. Výstavba modelů v B-J metodologii 45 pak můžeme psát ^Yn = (Cn - In)(Yn - ^Yn) . Postupně upravujme Yn = Yn - ^Yn + ^Yn = (Yn - ^Yn) + (Cn - In)(Yn - ^Yn) = Cn(Yn - ^Yn) . Tento výsledek použijme při vyjádření varianční matice n = EYnYn = E Cn(Yn - ^Yn)(Yn - ^Yn) Cn = CnE (Yn - ^Yn)(Yn - ^Yn) Cn. Nyní počítejme E (Yn-^Yn)(Yn-^Yn) = E(Y1-^Y1)2 v0 E(Y1-^Y1)(Y2-^Y2) =0 E(Y1-^Y1)(Yn-^Yn) =0 E(Y2-^Y2)(Y1-^Y1) =0 E(Y2-^Y2)2 v1 E(Y2-^Y2)(Yn-^Yn) =0 ... ... ... ... E(Yn-^Yn)(Y1-^Y1) =0 E(Yn-^Yn)(Yn-1-^Yn-1) =0 E(Yn-^Yn)2 vn-1 = diag{v0, . . . , vn-1} = Dn. Takže n = CnDnCn . Počítejme dále Yn-1 n Yn = (Yn - ^Yn) Cn [CnDnCn] -1 Cn(Yn - ^Yn) = (Yn - ^Yn) Cn (Cn) -1 D-1 n C-1 n Cn(Yn - ^Yn) = (Yn - ^Yn) D-1 n (Yn - ^Yn) = n j=1 (Yj - ^Yj)2 vj-1 . Dále zřejmě platí |n| = |CnDnCn| = |Cn| =1 |Dn| |Cn| =1 = v0v1 vn-1 . Takže věrohodnostní funkce náhodného vektoru Yn je tvaru L(n) = (2)- n 2 (v0v1 vn-1)- 1 2 exp -1 2 n j=1 (Yj - ^Yj)2 vj-1 . Nechť {Yt, t Z} je kauzální a invertibilní ARMA proces {Yt, t Z} ARMA(p, q) : (B)Yt = (B)t t N(0, 2 ). 1.8. Výstavba modelů v B-J metodologii 46 Ukazuje se (viz Brockwell, Davis (1978), str. 168/170), že k velkému zjednodušení jednokrokové predikce dojde, pokud inovační algoritmus aplikujeme ne na Yt, ale na následující transformo- vaný proces Wt = -1 Yt t = 1, . . . , m; m = max(p, q) -1 (B)Yt t > m. Poznamenejme nejprve, že zřejmě sp{Y1, . . . , Yn} = sp{W1, . . . , Wn} n 1. Označme Wj+1 = 0 = ^Y1 j = 0, Psp{W1,...,Wj}Wj+1 j 1. Pak platí Wt = -1 ^Yt t = 1, . . . , m; m = max(p, q), -1 ^Yt - 1Yt-1 - - pYt-p t > m, takže Yt - ^Yt = (Wt - Wt). Při aplikaci inovačního algoritmu na Wt dostaneme n,j a střední kvadratické chyby, které označme rj. Pak z předchozích vztahů vyplývá, že platí ^Yn+1 = n j=1 n,j(Yn+1-j - ^Yn+1-j) 1 n < m, 1Yn + + pYn+1-p + q j=1 n,j(Yn+1-j - ^Yn+1-j) n m. a vn = E(Yn+1 - ^Yn+1)2 = 2 E(Wn+1 - Wn+1)2 = 2 rn . Takže věrohodnostní funkce náhodného vektoru Yn je tvaru L(, , 2 ) = (22 )- n 2 (r0r1 rn-1)- 1 2 exp - 1 22 n j=1 (Yj - ^Yj)2 rj-1 . Pokud položíme ln L 2 = 0, a budeme předpokládat, že ^Yj a rj jsou nezávislé na 2 , dostaneme ^2 = 1 n S(^, ^) , kde S(^, ^) = n j=1 (Yj - ^Yj)2 rj-1 a ^ a ^ jsou hodnoty, které minimalizují tzv. redukovaný logaritmus věrohodnostní funkce l(^, ^) = ln 1 n S(^, ^) + 1 n n j=1 ln rj-1 . 1.8. Výstavba modelů v B-J metodologii 47 Poznámka 1.8.7. Alternativou k maximalizaci L(, , 2 ) je minimalizace váženého součtu čtverců S(, ) = n j=1 (Yj - ^Yj)2 rj-1 , přičemž ~2 = 1 n - p - q S(~, ~) a platí S(~, ~) ~2 A 2 (n - p - q). (viz Brockwell, Davis (1987), §8.9). Takto získané odhady se nazývají odhady metodou nejmenších čtverců. což vede k sys- tému nelineárních rovnic. Pokud chceme zkoumat asymptotické vlastnosti maximálně věrohodných odhadů, musíme zesílit předpoklady: nechť {Yt, t Z} je kauzální a invertibilní ARMA proces {Yt, t Z} ARMA(p, q) : (B)Yt = (B)t t IID(0, 2 ) a nechť (z) a (z) nemají společné kořeny. Pak, označíme-li maximálně věrohodný odhad neznámých parametrů = ( , , 2 ) symbolem ^MLE = (^ , ^ , ^2 ) , platí n ^MLE - A Nn+p+1(0, V ()), kde V () = 2 EUtUt EUtVt EVtUt EVtVt -1 , přičemž Ut = (Ut, . . . , Ut+1-p) Vt = (Vt, . . . , Vt+1-q) a {Ut, t Z} i {Vt, t Z} jsou autoregresní procesy (B)Ut = t (B)Vt = t (viz Brockwell, Davis (1987), §8.9). 1.9. Výstavba modelů a predikce v ARIMA procesech 48 1.9 Výstavba modelů a predikce v ARIMA procesech Až dosud jsme uvažovali pouze o procesech (slabě) stacionárních. Stacionární procesy se však v realitě vyskytují zřídka. Rozlišují se dva druhy nestacionarity: ˇ ve střední hodnotě ˇ v rozptylu. 1.9.1 Procesy nestacionární ve střední hodnotě Nyní je třeba si vysvětlit a odlišit pojmy Deterministický trend , kdy nestacionaritu ve střední hodnotě chápeme jako funkci času, tj. k jeho mo- delování použijeme například polynomický trend f(t) = 0 + 1t + + dtd periodický trend f(t) = + p j=1(jcosjt + jsinjt) Stochastický trend . U procesů ARMA jsem požadovali, aby všechny kořeny (z) = 1 - 1z - - pzp ležely vně jednotkové kružnice, tj. proces bude kauzální. Pokud nějaký kořen leží - na jednotkové kružnici mluvíme o procesu nestacionárním se stochastickým trendem - vně jednotkové kružnice mluvíme o procesu nestacionárním explozivního typu. Nestacionární proces obsahující stochastický trend lze převést na stacionární diferencová- ním. Zaveďme proto tzv. diferenční operátor: Yt = Yt - Yt-1 = (1 - B)Yt 2 Yt = (Yt) = (Yt - Yt-1) = (Yt - Yt-1) - (Yt-1 - Yt-2) = Yt - 2Yt-1 + Yt-2 = (1 - B)2 Yt ... d Yt = (1 - B)d Yt 1.9. Výstavba modelů a predikce v ARIMA procesech 49 Nestacionární proces se stochastickým trendem nazýváme integrovaným smíšeným mo- delem ARIMA(p, d, q) . Formálně jej zapíšeme pomocí operátoru zpětného chodu takto: ARIMA(p, d, q) : (B)(1 - B)d Yt = (B)t a položíme-li Wt = (1 - B)d Yt, pak Wt je stacionární ARMA(p, q). Zvláštní případy ARIMA(p, d, q) p d q Zkratka Název 0 IMA(d, q) Integrovaný proces klouzavých součtů 0 0 MA(q) Proces klouzavých součtů 0 ARI(p, d) Integrovaný autoregresní proces 0 0 AR(p) Autoregresní proces 0 0 I(d) Integrovaný proces 0 1 0 I(1) Náhodná procházka (random walk) Náhodná procházka: Yt = Yt-1 + t t V W(0, 2 ) Proces náhodné procházky je limitním případem procesu AR(1), kde 1 = 1. Pro náhodnou procházku platí: ˇ Hodnoty ACF = (k) budou klesat velmi pomalu (lineárně). ˇ Hodnoty PACF = (k) jsou logicky velmi podobné procesu AR(1). V praxi se používá obecnější tvar: Yt = + Yt-1 + t R. Potom, pokud budeme postupně upravovat, dostaneme Yt = + Yt-1 + t = Yt-2 + 2 + t + t-1 = = Y0 + t deterministický lineární trend + t j=1 j . 1.9. Výstavba modelů a predikce v ARIMA procesech 50 Poznámka 1.9.1. Tvary ACF = (k) a PACF = (k) procesů ARIMA(p, d, q) a náhodné procházky I(1) jsou prakticky totožné. Přítomnost jednotkových kořenů způsobuje ,,zakrytí prakticky všech identifikačních detailů těchto funkcí. Poznámka 1.9.2. Operátor (B) = (B)(1 - B)d se někdy nazývá zobecněný autoregresní operátor. Pokud (B) chápeme jako polynom v proměnné B, pak vzhledem ke kauzalitě modelu (1 - B)d Wt = (B)t má (B) právě p kořenů ležících vně jednotkového kruhu a d kořenů rovných 1. V praxi se nejprve diferencováním časové řady získá stacionární řada Wt a pro ni se vybuduje proces ARMA(p, q). Pokud jsme původně měli Y1, . . . , Yn, po diferencování zůstanou Wd+1, . . . , Wn. Poznámka 1.9.3. ARIMA(p, d, q) nemá smysl centrovat, neboť platí: d (Yt - Y ) = d Yt. Poznámka 1.9.4. Kromě trendů vyžadujících stochastické modelování mohou ARIMA modely zachytit i čistě deterministické trendy, pokud provedeme takovéto zobecnění ARIMA(p, d, q) modelů: ARIMA(p, d, q) : (B)(1 - B)d Yt = + (B)t R; , Pak této definici vyhovují procesy tvaru: 0 + 1t + + dtd polynomický trend řádu d + Yt . s využitím poznatků o diferencování polynomů lze totiž psát: (B)(1 - B)d (0 + 1t + + dtd + Yt) = (B)(d!d) =(1-1--p)d!d + (B)(1 - B)d Yt = + (B)t. 1.9. Výstavba modelů a predikce v ARIMA procesech 51 1.9.2 PROCESY NESTACIONÁRNÍ V ROZPTYLU Není-li splněna podmínka neměnnosti rozptylu v čase, je proces nestacionární v roz- ptylu. Takovýto proces je ovšem třeba nejprve vhodně transformovat. Vysvětleme si stručně pojem transformace stabilizující rozptyl. Situace nestabilního rozptylu nastává především v případě, kdy náhodná veličina Yt má rozdělení, které závisí na jediném parametru t, který obecně nemusí mít pro všechna t stejnou hodnotu. Předpokládejme, že tento parametr je zvolen tak, aby platilo Et Yt = t. Ve většině případů (ne však u normálního rozdělení) na t závisí i rozptyl veličiny Yt, takže můžeme psát Dt Yt = 2 (t). Přitom (t) bývá obvykle hladká funkce proměnné t. Protože t může souviset s časem t, není splněna podmínka neměnnosti rozptylu v čase. Vzniká tedy otázka, zda lze najít netriviální funkci g tak, aby náhodná veličina Zt = g(Yt) měla rozptyl nezávisející na t. (Požadavkem netriviality se vylučují konstantní funkce g, které by vedly k veličinám s nulovým rozptylem). Uvedená úloha v obecném případě nemá řešení. Používá se však určitých aproximací, které se ukázaly velmi užitečné. Pokud se zabýváme jen dostatečně hladkými funkcemi g, z Taylorova rozvoje dostaneme aproximaci g(Yt) g(t) + g (t)(Yt - t). Potom střední hodnotu lze aproximovat takto Eg(Yt) E [g(t) + g (t)(Yt - t)] = g(t) a rozptyl Dt [g(Yt)] [g (Yt)] 2 Dt Yt = [g (t)] 2 2 (t). Chceme, aby po transformaci byl rozptyl konstantní a nezávisel na střední hodnotě, tj. c2 = Dt [g(Yt)] = [g (t)] 2 2 (t) g (t) = c 2(t) , kde c je nějaká konstanta. Odtud snadno dostaneme tvar transformace stabilizující rozptyl g(t) = c 1 (t) dt + K. Konstanty c a K se volí tak, aby funkce g vypočtená podle předchozího vzorce měla výhodný tvar. Ukázalo se, že funkce g vypočtená podle předchozího vzorce nejen výrazně stabilizuje roz- ptyl, takže rozptyl Dg(Yt) závisí na t jen velmi málo, ale zároveň také rozdělení náhodné veličiny Zt = g(Yt) bývá již velmi blízké normálnímu, i když třeba samotné rozdělení veličiny Yt je výrazně nenormální. 1.9. Výstavba modelů a predikce v ARIMA procesech 52 1.9.3 VOLBA ŘÁDU MODELU Vhodná volba řádu modelu je důležitou součástí výstavby modelů v Box-Jenkinsonově meto- dologii. K identifikaci typu modelu nám poslouží především průběh odhadnuté autokorelační funkce ACF = (k) a parciální autokorelační funkce PACF = (k). Formálnější přístup při rozhodování o výběru řádů p a q je však založen na jistých kriteri- álních funkcích. Všimněme si nejprve kriteriální funkce odvozené pro AR(p) modely. Odvození kriteria FPE (Final prediction error) Mějme dvě nezávislé realizace (X1, . . . , Xn) a (Y1, . . . , Yn) procesu AR(p) s koeficienty p = (1, . . . , p) . Pomocí prvních realizací odvodíme maximálně věrohodné odhady ^p = ( ^1, . . . , ^p) . Definujme Yp = (Yn, Yn-1, . . . , Yn+1-p) a uvažujme jednokrokovou predikci ^Yn+1 = ^1Yn + + ^pYn+1-p = ^p Yp. Počítejme střední kvadratickou chybu této predikce: ^vn = E(Yn+1 - ^Yn+1)2 = E(Yn+1 - ^1Yn - - ^pYn+1-p)2 = E ^Yn+1 - 1Yn - - pYn+1-p - ( ^1 - 1)Yn - - ( ^p - p)Yn+1-p 2 = E n+1 - (^p - p) Yp 2 = E2 n+1 =2 +(^p - p) EYpYp p (^p - p) neboť E(n+1Yn-j) = 0 pro j = 0, . . . , p - 1. Protože (n)(^p - p) A Np(0, 2 -1 p ) nebo-li (^p - p) A Np(0, 2 n -1 p ). Odtud dostaneme, že n 2 (^p - p) p(^p - p) = n 2 Zp A 2 (p). 1.9. Výstavba modelů a predikce v ARIMA procesech 53 Protože platí E n 2 Zp = p, ihned dostaneme EZp = 2 p n a odtud ^vn 2 1 + p n . Jestliže ^2 je maximálně věrohodným odhadem neznámého parametru 2 , pak pro n platí n ^2 2 A 2 (n - p) a protože E n ^2 2 = n - p, pak E ^2 = n-p n 2 2 = n n-p E ^2 . Pokud ve výpočtu ^vn nahradíme 2 jeho nestranným odhadem ^2 , dostaneme následující kri- térium: FPE = ^2 n n-p 1 + p n = ^2 n + p n - p . Mnohem obecnějším kritériem, které není odvozeno pro AR(p) proces, je tzv. AIC krité- rium: AIC = nl(^pq, ^pq) + 2(p + q), nebo jeho korigovaná verze (neboť AIC nadhodnocuje řád modelu) AICC = -2 ln L(^pq, ^pq, S(^pq, ^pq)) + 2(p + q + 1)n n - p - q - 2 . Dalším kritériem je tzv. BIC kritérium: BIC = (n - p - q) ln n ^2 n - p - q + (p + q) ln P Y 2 t -n ^2 p+q . 1.9. Výstavba modelů a predikce v ARIMA procesech 54 1.9.4 Verifikace modelu ­ analýza reziduí Po odhadu neznámých parametrů ARMA(p, q) modelu metodou maximální věrohodnosti označme rezidua: Wt = Yt - ^Yt(^, ^) rt-1(^, ^) t = 1, . . . , n. Poznámka 1.9.5. Existují i další rezídua, např. ^t = ^ -1 (B) ^(B)Yt t = 1, . . . , n. Velikost těchto reziduí však závisí na měrných jednotkách náhodné proměnné Yt, proto se doporučuje používat raději Wt, která jsou standardizovaná a nezávislá na měřítku. Protože rezidua Wt by měla být pro n (pokud byl model vhodně zvolen) přibližně bílým šumem WN(0, 2 ) (popř. IID(0, 2 )), ihned se nabízí myšlenka, analyzovat nejprve jeho výběrovou autokorelační funkci ^W (h) = n-h t=1 (Wt - W)(Wt+h - W) n t=1(Wt - W)2 , kde W = n t=1 Wt. Celá řada článků je věnována odvození asymptotických intervalů spolehlivosti pro tuto výbě- rovou autokorelační funkci. Mnohem jednodušší je použití tzv. Portmanteaovy statistiky QW = h j=1 ^2 W (j) A 2 (h - p - q) pro h N. Existují také jisté modifikace této statistiky, např. (Ljung, Box, 1978) ~QW = n(n + 2) h j=1 ^2 W (j) n-j A 2 (h - p - q) pro h N, nebo (Granger, Anderson, 1978) ~QW W = n(n + 2) h j=1 ^W W (j) n-j A 2 (h) pro h N. 1.10. Modelování sezónnosti pomocí SARIMA modelů 55 1.10 Modelování sezónnosti pomocí SARIMA modelů Sezónnost je v Box-Jenkinsonově metodologii stejně jako trend modelována stochasticky. Nejprve zaveďme sezónní diferenční operátor o délce L > 0: LYt = Yt - Yt-L = (1 - BL )Yt 2 LYt = L(LYt) = L(Yt-Yt-L) = (Yt -Yt-L)-(Yt-L -Yt-2L) = Yt -2Yt-L +Yt-2L = (1-BL )2 Yt ... D L Yt = (1 - BL )D Yt Při konstrukci se uvažuje způsobem, který budeme demonstrovat pomocí následujícího příkladu: Nechť časová řada {Yt} vykazuje sezónnost o délce L = 12. 1. Zkonstruujeme nejprve ARIMA(P1, D1, Q1) model pro řadu lednových měření, tj. pro {S1 t = B12 Yt} 1(B12 )D1 12 Yt = 1(B12 ) (1) t ARIMA(P1, D1, Q1) kde časový index t odpovídá lednovým obdobím a o t se budeme zajímat později. Přitom 1(B12 ) = 1 - 1,1B12 - - 1,P1 B12P1 je tzv. sezónní autoregresní operátor SAR(P1) 1(B12 ) = 1 + 1,1B12 + + 1,Q1 B12Q1 je tzv. sezónní operátor klouzavých součtů SMA(Q1) D1 12 = (1 - BL )D1 je tzv. sezónní diferenční operátor SI(D1) 2. Podobné modely zkonstruujeme pro ostatní měsíce: 2(B12 )D2 12 Yt = 2(B12 ) (2) t ARIMA(P2, D2, Q2) ... 12(B12 )D12 12 Yt = 12(B12 ) (12) t ARIMA(P12, D12, Q12) 3. Předpokládejme přitom, že tyto modely jsou pro jednotlivé měsíce přibližně stejné, tj. P1 P12 P Q1 Q12 Q D1 D12 D 1(B12 ) 12(B12 ) (B12 ) 1(B12 ) 12(B12 ) (B12 ) 1.10. Modelování sezónnosti pomocí SARIMA modelů 56 4. Náhodné veličiny (j) t (j = 1, . . . , 12) by však v těchto modelech měly být pro různé měsíce mezi sebou korelované, neboť by měl existovat např. vztah mezi lednovými a únorovými hodnotami. Předpokládejme proto, že také řada t je popsána modelem ARIMA(p, d, q) tvaru (B)d t = (B)t ARIMA(p, d, q) kde t WN(0, 2 ) je bílý šum. 5. Spojme předchozí dva modely do jediného tzv. multiplikativního sezónního modelu řádu(p, d, q) × (P, D, Q)L (B)(BL )d D L Yt = (B)(BL )t SARIMA(p, d, q) × (P, D, Q)L L = 12. Příklad 1.10.1: Model SARIMA(0, 1, 1) × (0, 1, 1)12 má tvar: 12Yt = (1 - B)(1 - B12 )Yt = (1 + 1B)(1 + 1B12 )t, nebo ekvivalentně Yt - Yt-1 - Yt-12 + Yt-13 = t + 1t-1 + 1t-12 + 11t-13. Poznámka 1.10.1. Existují také aditivní sezónní modely, které se však používají jen zřídka. Jako příklad lze uvést model Yt = t + 1t-1 + 12t-12 + 13t-13. Výstavba sezónních modelů Označme řád běžného diferencování na odstranění trendu jako d a D jako řád diferencování na odstranění sezónnosti. V praktických situacích většinou d = 0, 1, 2 a D = 0, 1. Dále nechť L je délka sezóny. Výstavba sezónních modelů probíhá ve třech stejných fázích jako pro modely ARIMA. Všimněme si pouze FÁZE IDENTIFIKACE MODELU, neboť ostatní dvě fáze jsou totožné. 1. Odhad parametrů d, D : (a) Provede se studium odhadnuté autokorelační funkce ACF = ^(k), neboť identi- fikuje přítomnost trendu. Doporučuje se prozkoumat 4L hodnot ^(k). ˇ Určení D : Má-li funkce ^(k) v bodech L, 2L, 3L . . . lokální maxima, pak (bez ohledu na její průběh mezi těmito časovými body) je nutné položit D = 1. To plyne z toho, že hodnoty ^(L), ^(2L), ^(3L), . . . představují odhadnuté hod- noty autokorelační funkce pro řady {Sj t = BL Yt}, j = 1, . . . , L modelu (BL )D L Yt = (BL )t ARIMA(P, D, Q), přičemž nestacionaritě tohoto ARIMA modelu odpovídá pomalý pokles autokorelační funkce ^(L), ^(2L), ^(3L), . . ., tj. tuto řadu je nutno diferenco- vat (s krokem L) a pokládáme D = 1. ˇ Určení d : Jestliže funkce rk klesá mezi body jL a (j + 1)L pouze přibližně lineárně, je třeba provést také běžné diferencování. 1.10. Modelování sezónnosti pomocí SARIMA modelů 57 (b) Čísla d, D se také někdy určují tak, že se hledá nejmenší číslo mezi odhadnutými hodnotami ^2 Y , ^2 Yt , ^LYt , ^2 LYt , . . . rozptylů dané řady a jejich diferencí. 2. Odhad parametrů p, P, q, Q : Po určení řádu d a D zkonstruujeme řadu Wt = d D L Yt, pro kterou je nutné identifikovat model tvaru (B)(BL )Wt = (B)(BL )t SARIMA(p, 0, q) × (P, 0, Q)L. Pro tento účel se použije odhadnutá ACF = ^(k) a PACF = ^(k) řady Wt. (a) MA­Homogenní modely ˇ Jestliže ACF funkce ^(k) je zhruba významně nenulová v bodech 1, . . . , q L - q, . . . , L + q 2L - q, . . . , 2L + q ... QL - q, . . . , QL + q přičemž mezi těmito body se neodlišují významně od nuly ˇ a funkce ^(k) v jednotlivých úsecích mezi body jL a (j +1)L vždy v absolutní hodnotě klesá (geometricky nebo po sinusoidě s geometricky klesající ampli- tudou) a zároveň klesá, když ji sledujeme v bodech L, 2L, 3L, . . . , pak položíme p = 0 a P = 0 , tj. budeme identifikovat odpovídající model pro řadu Wt jako Wt = (B)(BL )t SARIMA(0, 0, q) × (0, 0, Q)L a tedy model pro řadu Yt jako d D L Yt = (B)(BL )t SARIMA(0, d, q) × (0, D, Q)L. (b) AR­Homogenní modely ˇ Jestliže naopak funkce ^(k) klesá v absolutní hodnotě (geometricky nebo po si- nusoidě s geometricky klesající amplitudou) v úsecích mezi body jL a (j+1)L a zároveň klesá, když ji sledujeme v bodech L, 2L, 3L, . . . 1.10. Modelování sezónnosti pomocí SARIMA modelů 58 ˇ a funkce ^(k) je zhruba významně nenulová v bodech 1, . . . , p L, . . . , L + p 2L, . . . , 2L + p ... PL, . . . , PL + p přičemž mezi těmito body se neodlišují významně od nuly, pak položíme q = 0 a Q = 0 , tj. budeme identifikovat odpovídající model pro řadu Wt jako (B)(BL )Wt = t SARIMA(p, 0, 0) × (P, 0, 0)L a tedy model pro řadu Yt jako (B)(BL )d D L Yt = t SARIMA(p, d, 0) × (P, D, 0)L. (c) Nehomogenní modely typu SARIMA(p, d, 0) × (0, D, Q)L nebo SARIMA(0, d, q) × (P, D, 0)L se většinou nepoužívají, neboť obvykle vedou při srovnání s předchozími tzv. ho- mogenní modely SARIMA(0, d, q) × (0, D, Q)L nebo SARIMA(p, d, 0) × (P, D, 0)L k odhadu neúnosně velkého počtu parametrů. (d) Identifikace obecných modelů SARIMA(p, d, q) × (P, D, Q)L, v nichž čísla p, q, P a Q mohou být vesměs nenulová, je již dosti komplikovanou zále- žitostí a obvykle zde hodně záleží na zkušenostech statistika, který analýzu provádí.