Marie Forbelská Stochastické modelování jednorozměrných časových řad II 0 50 100 150 200 250 300 -4 -2 0 2 4 6 0 10 20 30 40 50 -0.5 0 0.5 1 0 10 20 30 40 50 -0.5 0 0.5 1 ARMA(2, 2) : Yt =0.5Yt-1 +0.2Yt-2 +t -0.4t-1 +0.3t-2 AR roots=3.81 1.31; MA roots=1.83 1.83; t W N(0, 1) ACF PACF Brno 2OO7 Příro dovědecká fa kulta Masa rykova Unive rzita Ústav matematiky a statistiky iii Obsah Seznam použitých symbolů 1 1 Box-Jenkinsova metodologie 3 1.1 Úvod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Základní pojmy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Operátor zpětného posunutí . . . . . . . . . . . . . . . . . . 4 1.2.2 Lineární proces . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3 Lineární filtry . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.4 Autokovarianční a autokorelační generující funkce stacionárních procesů . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.5 Definice ARMA procesu . . . . . . . . . . . . . . . . . . . . 12 1.2.6 Kauzalita . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.7 Invertibilita . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.2.8 Vícenásobná reprezentace MA(q) procesů . . . . . . . . . . 28 1.2.9 Společné kořeny polynomů (z) a (z) . . . . . . . . . . . 29 1.2.10 Nutná a postačující podmínka kauzality a invertibility ARMA procesu. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.2.11 Střední hodnota, rozptyl, autokovarianční a autokorelační funkce procesů ARMA(p, q) . . . . . . . . . . . . . . . . . . 30 1.2.12 Spektrální hustota ARMA(p, q) procesů . . . . . . . . . . . 32 1.3 Stacionární procesy a nejlepší lineární predikce . . . . . . . . . . . 35 1.4 Rekurentní metody pro výpočet nejlepší lineární predikce . . . . . 38 1.4.1 Důsledek Durbin-Levinsonova algoritmu . . . . . . . . . . . 43 1.5 Parciální autokorelační funkce (PACF) . . . . . . . . . . . . . . . . 44 1.6 Inovační algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 1.6.1 Jednokroková nejlepší lineární predikce v AR(p) . . . . . . 50 1.6.2 Vícekroková nejlepší lineární predikce v AR(p) . . . . . . . 51 1.6.3 PACF pro AR(p), MA(q) a ARMA(p, q) . . . . . . . . . . 52 1.6.4 Jednokroková nejlepší lineární predikce v MA(q) . . . . . . 53 1.6.5 Nejlepší lineární predikce v ARMA(p, q) . . . . . . . . . . . 54 1.7 Výstavba modelů v B-J metodologii . . . . . . . . . . . . . . . . . 57 1.7.1 Odhady v ARMA procesech . . . . . . . . . . . . . . . . . . 57 1.7.2 Yuleovy-Walkerovy rovnice a odhad parametrů v AR(p) . . 61 iv 1.7.3 Předběžné odhady v AR(p) a Durbin-Levinsův algoritmus. 62 1.7.4 Předběžné odhady v MA(q) a inovační algoritmus. . . . . . 63 1.7.5 Předběžné odhady v ARMA(p, q) procesu. . . . . . . . . . 64 1.7.6 Maximálně věrohodné odhady. . . . . . . . . . . . . . . . . 66 1.8 Výstavba modelů a predikce v ARIMA procesech . . . . . . . . . . 71 1.8.1 Procesy nestacionární ve střední hodnotě . . . . . . . . . . 71 1.8.2 Procesy nestacionární v rozptylu . . . . . . . . . . . . . . . 73 1.8.3 Volba řádu modelu . . . . . . . . . . . . . . . . . . . . . . . 81 1.8.4 Verifikace modelu ­ analýza reziduí . . . . . . . . . . . . . . 83 1.9 Modelování sezónnosti pomocí SARIMA modelů . . . . . . . . . . 84 2 Dynamické lineární modely 89 2.1 Motivační příklad . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 2.2 Stavově-prostorové modely . . . . . . . . . . . . . . . . . . . . . . . 91 2.3 Stacionární stavově-prostorové modely . . . . . . . . . . . . . . . . 94 2.4 Nejlepší lineární predikce pomocí projekce náhodných vektorů druhého řádu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 2.5 Kalmanův iterační proces . . . . . . . . . . . . . . . . . . . . . . . 105 Literatura 107 1 Seznam použitých symbolů N množina přirozených čísel Z množina celých čísel R množina reálných čísel Rm reálný m-rozměrný euklidovský prostor a, b, . . . reálná čísla a, b, . . . vektory reálných čísel (sloupcové) a euklidovská norma vektoru a, tj. a = aa A, B, . . . matice In jednotková matice typu (n × n) 0 nulový vektor nebo nulová matice 1n n-rozměrný sloupcový vektor samých jedniček h(A) hodnost matice A transponovaná matice A-1 inverzní matice Apseudoinverzní matice, tj. A= AA- A |A|, det(A) determinant matice A B C D bloková matice ln x přirozený logaritmus se základem e prostor elementárních jevů A -algebra Bm borelovská -algebra v m-rozměrném euklidovském prostoru IA indikátor množiny A A B sjednocení množin A B průnik množin A - B rozdíl množin prázdná množina x, y skalární součin P(A) pravděpodobnost jevu A (, A, P) pravděpodobnostní prostor X, Y, . . . náhodné veličiny x, y, . . . realizace náhodných veličin X, Y, . . . náhodný vektor x, y, . . . realizace náhodných vektorů 2 EX střední hodnota náhodné veličiny X C(X, Y ) kovariance náhodných veličin X a Y DX rozptyl náhodné veličiny X C(X, Y) kovarianční matice náhodných vektorů X a Y DX varianční matice náhodného vektoru X X L(, 2 ) náhodná veličina X s rozdělením se střední hodnotou a rozptylem 2 X N(, 2 ) náhodná veličina X se řídí normálním rozdělením se střední hodnotou a rozptylem 2 Rs(a, b) rozdělení rovnoměrně spojité na intervalu (a, b), a, b R, a 0, L N, tj. LXt = Xt - Xt-L fX() spektrální hustota náhodného procesu {Xt, t T } (t) autokovarianční funkce stacionárního procesu (t) autokorelační funkce stacionárního procesu (t) parciální autokorelační funkce stacionárního procesu ARMA(p, q) ARMA proces řádu p a q (p, q N) AR(p) AR (autokorelační) proces řádu p (p N) MA(q) MA proces řádu q (q N) ARIMA(p, d, q) smíšený proces řádu p, d a q (p, d, q N) 3 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 dekompozič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 pro- cesu * 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 zadavatele, ž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ů). 4 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.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 5 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. 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 . 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 nekorel. = s - j = s + t - h h = j + t = 2 j=0 jj+t . 6 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 j=0 2 j < . Pokud zavedeme funkci (z) = j=0 jzj , pak podmínka j=0 2 j < implikuje, že funkce (z) je holomorfní uvnitř kružnice |z| < 1. Takže podmínku stacionarity lze vyslovit i pomocí podmínky (z) je holomorfní pro |z| < 1, přičemž j=0 2 j < . Oba požadavky budou splněny, pokud bude platit (z) je holomorfní uvnitř a na jednotkové kružnici . 7 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 procesu. 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í hustotu 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 (z) = j=- jzj , pak můžeme psát fY () = f() e-i 2 = 2 2 e-i 2 = 2 2 ei 2 . 8 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 (1) t = -1 j=jXt-j a Y (2) t = j=0 jXt-j, protože pak bude platit Yt = Y (1) t + Y (2) t . Označme X(h) = EXtEt+|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 . 9 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í hustoutou 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 spektrální hustotou fY () = fX() e-i 2 , kde (z) = j=- jzj |z| 1 se nazývá generující funkce filtru a () = e-i přenosová funkce filtru. Důkaz. Stacionaritu jsme dokázali v předchozí větě. Nyní 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. Označme (z) = j=- jzj pro |z| 1. 10 Pak, protože platí Y (t) = - eit fY ()d, dostaneme fY () = fX() e-i 2 = fX() ei 2 . 1.2.4 Autokovarianční a autokorelační generující funkce stacionárních procesů Při vyšetřování vlastností autokovarianční a autokorelační funkce starionárních procesů je užiteč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 . Bílý šum Příklad 1.2.8. 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 R(z) = k=- (k)zk = 1. Zobecněný lineární proces Příklad 1.2.9. 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. 11 Odtud autokovarianční generujíci funkce a spektrální hustota jsou rovny GX (z) = k=- (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=- hz-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 , kde X(z) = j=- jzj pro |z| 1 se nazývá generující funkce filtru a X() = X e-i přenosová funkce filtru (zobecněný lineární proces je také filtrem). Lineární filtr Příklad 1.2.10. 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+k-h = j=- h=jhX(k + j - h) , 12 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 , kde Y (z) = j=- jzj pro |z| 1 se nazývá generující funkce filtru a Y () = Y e-i přenosová funkce filtru. 1.2.5 Definice ARMA procesu Definice 1.2.11. 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 13 1.2.6 Kauzalita Dříve než zavedeme pojem kauzality, všimněme si blíže AR(1) procesu. Autoregresní procesy prvního řádu. Yt = 0.5Yt-1 + t, t N(0, 1) 0 50 100 150 200 250 300 -4 -3 -2 -1 0 1 2 3 4 Yt = -0.5Yt-1 + t, t N(0, 1) 0 50 100 150 200 250 300 -4 -3 -2 -1 0 1 2 3 4 Yt = 0.85Yt-1 + t, t N(0, 1) 0 50 100 150 200 250 300 -6 -4 -2 0 2 4 6 Yt = -0.25Yt-1 + t, t N(0, 1) 0 50 100 150 200 250 300 -3 -2 -1 0 1 2 3 4 Obrázek 1.1: Ukázky autoregresních procesů 1. řádu 14 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 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 . 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|) = E(Ys Ys+|t|) = E j=0 j 1s-j h=0 h 1 s+|t|-h = j=0 h=0 jhE s-js+|t|-h nekorel. = s - j = s + |t| - h h = j + |t| = 2 j=0 j 1 j+|t| 1 = |t| 1 1 - 2 1 2 . 15 Autokorelační funkce (ACF) je pak tvaru (t) = (t) (0) = |t| 1 Pomocí generující funkce filtru AR(1)(z) = j=0 j 1yj = 1 1 - 1z pro |z| < 1 a |1| < 1 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 . ACF pro Yt = 0.5Yt-1 + t -10 -8 -6 -4 -2 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 ACF pro Yt = -0.5Yt-1 + t -10 -8 -6 -4 -2 0 2 4 6 8 10 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Obrázek 1.2: Ukázky autokorelačních funkcí fAR(1)() pro Yt = 0.5Yt-1 + t -3 -2 -1 0 1 2 3 0.1 0.2 0.3 0.4 0.5 0.6 fAR(1)() pro Yt = -0.5Yt-1 + t -3 -2 -1 0 1 2 3 0.1 0.2 0.3 0.4 0.5 0.6 Obrázek 1.3: Ukázky spektrálních hustot (2) Dále řešme případ, kdy |1| > 1. Použijeme-li vztah Yt-1 = 1 1 Yt - 1 1 t 16 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 , stejně jako v předchozím případu - j=0 1 j 1 t-j konverguje podle kvadratické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.12. ARMA proces Yt ARMA(p, q) se nazývá kauzální, jestliže existuje absolutně konvergentní 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.13. 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. 17 Autoregresní proces p -tého řádu: AR(p) : (B)Yt = t 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). Budeme hledat takovou absolutně konvergentní posloupnosti čísel = {j} j=0, tak aby Yt = j=0 jt-j = (B)t byl kauzální proces. Takže postupně odvo- zujme t = (B) Yt =(B)t = (B)(B)t, tj. (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 18 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é. 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) 19 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 j2 (j) = 2 , kde (j) = 1 j = 0, 0 jinak. 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) Připomeňme, že s využitím vztahu EYt = 0, je (k) = C(Yt, Yt-k) = E(Yt - EYt)(Yt-k - EYt-k) = EYtYt-k. 20 Spočtěme EYt-kt = E j=0 jt-j-k t = j=0 jEt-j-kt = j=0 j2 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. 21 2. Nechť 1 1 , . . . , 1 p jsou různé jednoduché kořeny, pak c1k 1 + + cpk p jsou řešením homogenní diferenční rovnice a c1, . . . , cp jsou konstanty, které jsou určeny počátečními podmínkami. 3. Je-li kořen 1 j dvojnásobný kořen, pak k j a 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 - 2tk j + kk j + 2k j - 2k j t = 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 m jsou různé kořeny s násobnostmi p1, . . . , pm, přičemž p = p1 + + pm, pak řešení homogenní diferenční rovnice (B)(k) = 0 je tvaru (k) = m 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) = m 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ů. 22 1.2.7 Invertibilita Víme, že kauzální autoregresní proces konečného řádu p lze vyjádřit 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í autoregresního procesu nekonečného řádu, tj. MA(q) AR() . Nejprve si všimneme jednoduchého případu, a to MA(1) procesu. Yt = t + 0.5t-1, t N(0, 1) 0 50 100 150 200 250 300 -3 -2 -1 0 1 2 3 Yt = t - 0.5t-1, t N(0, 1) 0 50 100 150 200 250 300 -4 -3 -2 -1 0 1 2 3 Yt = t + 0.85t-1, t N(0, 1) 0 50 100 150 200 250 300 -4 -3 -2 -1 0 1 2 3 Yt = t - 0.25t-1, t N(0, 1) 0 50 100 150 200 250 300 -3 -2 -1 0 1 2 3 4 Obrázek 1.4: Ukázky MA procesů prvního řádu 23 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 tYt + 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 . Dále využijme vztah t-1 = 1 1 Yt - 1 1 t a označme 1 = -. Za předpokládejme, že platí |1| = || > 1 , 24 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. I když posloupnost N j=1 1 j Yt+j konverguje pro N také k t, tento rozvoj nemá praktický smysl, neboť t je vyjadřena pomocí budoucích hodnot {Ys, s > t}. Dále si všimněme dalšího důležitého faktu, že pokud platí |1| > 1 , a uvažujeme-li dva procesy (1) Yt = t + 1t-1 t WN(0, 2 ), (2) Xt = t + 1 1 t-1 t WN(0, 2 12 ), pak oba dva procesy mají stejné první a druhé momenty, neboť EYt = E (t + 1t-1) = Et + 1Et-1 = 0 , EXt = E t + 1 1 t-1 = Et + 1 1 Et-1 = 0 , 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 = pokud: t = t + k t = t + k - 1 t - 1 = t + k t - 1 = t + k - 1 pak: k = 0 k = 1 k = -1 k = 0 = 2 + 2 12 = 2 (1 + 2 1) k = 0 12 k = 1 0 jinak 25 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 = 2 12 + 1 2 1 2 12 = 2 (1 + 2 1) k = 0 1 1 2 12 = 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 pozorovatelný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.14. ARMA proces Yt ARMA(p, q) se nazývá invertibilní, jestliže existuje absolutně konvergentní 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 invertibilní MA proces řádu q. 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)r1 . . . (1 - kz)rk , kde z01 = 1 1 , . . . , z0k = 1 k jsou rozdílné (reálné či komplexní) kořeny polynomu (z), r1, . . . , rk je jejich násobnost (přičemž platí r1 + + rk = q). Nyní budeme hledat taková absolutně konvergentní = {j} j=0(tj. j=0 |j| < ), aby t = j=0 jYt-j 26 byl invertibilní 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 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 jsEs-js+t-h = s - j = s + t - h h = j + t = 2 q j=0 jj+t 27 Protože 0 = 1 a j = 0 pro j > q, dostáváme (t) = 8 >>>>>>>>>>>>>< >>>>>>>>>>>>>: 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 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 fY () = 2 2 e-i 2 . Autokorelační funkce (t) pro Yt = t - 0.4t-1 + 0.2t-2 - 0.3t-3 -5 -4 -3 -2 -1 0 1 2 3 4 5 -0.5 0 0.5 1 Spektrální hustota fMA(3)() pro Yt = t - 0.4t-1 + 0.2t-2 - 0.3t-3 -3 -2 -1 0 1 2 3 0.1 0.2 0.3 0.4 0.5 Autokorelační funkce (t) pro Yt = t - 0.2t-1 + 0.1t-2 + 0.3t-3 -5 -4 -3 -2 -1 0 1 2 3 4 5 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Spektrální hustota fMA(3)() pro Yt = t - 0.2t-1 + 0.1t-2 + 0.3t-3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 0.3 Obrázek 1.5: Ukázky autokorelačních funkcí a spektrálních hustot pro MA(3) procesy. 28 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 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ů. 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). 29 1.2.9 Společné kořeny polynomů (z) a (z) Mějme Yt ARMA(p, q) : Yt - 1Yt-1 - - pYt-p = t + 1t-1 + + qt-q, kde 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) tj. (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.2.10 Nutná a postačující podmínka kauzality a invertibility ARMA procesu. 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.2.15. Nechť (B) a (B) nemají společné kořeny. Pak (i) Yt ARMA(p, q) : (B)Yt = (B)t je kauzální (z) = 0 pro z C |z| 1. (ii) 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 30 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.2.11 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 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 = i2 (přičemž 0 = 1). Po dosazení do rovnice (A6) dostaneme (0) - 1(1) - - p(p) = 2 (1 + 11 + . . . + qq) (A7). 31 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 dosta- neme (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 s 0 platí EtYt-s = E t j=0 jt-s-j = j=0 jEtt-s-j = 0. Spočtěme pro i = 0, 1, . . . , q Et-iYt-k = Et-i j=0 jt-j-k = j=0 jEt-it-j-k = t - i = t - j - k j = i - k = 2 i-k k i přičemž 0 = 1 k q 0 k > i neboť j = 0 pro j < 0 Uvážíme-li, že j = 0 pro j < 0, potom pro 0 k max(p, q + 1) platí (k) - 1(k - 1) - - p(k - p) = 2 (k + k+11 + + qq-k) (A9) a pro k > max(p, q + 1) platí (k) - 1(k - 1) - - p(k - p) = 0 . 32 Podělme obě strany rovnice (A9) 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) = m j=1 pj -1 s=0 cjsks k j , kde 1 1 , . . . , 1 m jsou různé kořeny s násobnostmi p1, . . . , pm, přičemž p = p1 + + pm a cjs je právě p konstant, které jsou určeny počátečními podmínkami. 1.2.12 Spektrální hustota ARMA(p, q) procesů Věta 1.2.16 (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 () = 2 2 e-i 2 | (e-i)| 2 pro -, . 33 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 () = 2 2 e-i 2 | (e-i)| 2 pro -, . 34 Na závěr tohoto odstavce jsou vykresleny příklady tří realizací AR, MA a ARMA procesů spolu s jejich teoretickými spektrálními hustotami. AR(2) : Yt = 0.5Yt-1 + 0.2Yt-2 + t, t N(0, 1) 0 50 100 150 200 250 300 -4 -2 0 2 4 MA(2) : Yt = t - 0.5t-1 - 0.2t-1, t N(0, 1) 0 50 100 150 200 250 300 -4 -2 0 2 4 ARMA(2, 2) : Yt = 0.5Yt-1 + 0.2Yt-2 + t - 0.4t-1 + 0.3t-1, t N(0, 1) 0 50 100 150 200 250 300 -4 -2 0 2 4 6 Obrázek 1.6: Ukázky realizací AR, MA a ARMA procesů fAR() = 2 2 1 |(e-i)|2 -3 -2 -1 0 1 2 3 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 fMA() = 2 2 e-i 2 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 0.3 fARMA() = 2 2 |(e-i )|2 |(e-i)|2 -3 -2 -1 0 1 2 3 0.2 0.4 0.6 0.8 1 1.2 1.4 pro -, . Obrázek 1.7: Ukázky spektrálních hustot AR, MA a ARMA procesů. 35 1.3 Stacionární procesy a nejlepší lineární pre- dikce 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.3.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 = EXY , 36 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 = EYn+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) 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.3.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. Důsledek 1.3.3. Označme Yn = (Yn, . . . , Y1) . 37 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 . 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 = nYnY nn a počítejme EYnY n = E 2 6 6 6 4 0 B B B @ Yn Yn-1 ... Y1 1 C C C A (Yn, . . . , Y1) 3 7 7 7 5 = 0 B B B @ EY 2 n EYnYn-1 EYnY1 EYn-1Yn EY 2 n-1 EYn-1Y1 ... ... ... ... EY1Yn EY1Yn-1 EY 2 1 1 C C C A = 0 B B B @ (0) (1) (n - 1) (1) (0) (n - 2) ... ... ... ... (n - 1) (n - 2) (0) 1 C C C A = 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 - 2 nEYnYn+1 + nEYnY nn = (0) - 2 nn + nnn = (0) - 2 n-1 n n + n-1 n n-1 n n = (0) - n-1 n n. Nyní definujme h-krokovou (lineární) predikci. 38 Definice 1.3.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á 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.4 Rekurentní metody pro výpočet nejlepší lineární predikce Věta 1.4.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,j (j = 1, . . . , n) a střední kvadratické chyby 39 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. 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 . 40 Počítejme postupně PMn-1 (Y1) = n-1,1Y2 + + n-1,n-1Yn a označíme-li n-1 = (n-1,1, . . . , n-1,n-1) , pak složky vektoru n-1 jsou pro k = 1, . . . , n - 1 řešením rovnic EY1Y1+k = E Y1PMn-1 (Y1) = n-1 j=1 n-1,jEY1+jY1+k, 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,1Y2 + + n-1,n-1Yn a označíme-li n-1 = (n-1,1, . . . , n-1,n-1) , pak složky vektoru n-1 jsou pro k = 1, . . . , n - 1 řešením rovnic EYn+1Yn+1-k = E Yn+1PMn-1 (Yn+1-k) = n-1 j=1 n-1,jEYn+1-jYn+1-k, 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. 41 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.3.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 - j), 42 pak a = n,n = Yn+1, Y1 - PMn-1 (Y1) Y1 - PMn-1 (Y1) 2 = (n) - n-1 j=1 n-1,j(n - j) vn-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) E = ˙ 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) ¸ | {z } =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). 43 1.4.1 Důsledek Durbin-Levinsonova algoritmu Důsledek 1.4.2. 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, 44 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) . 1.5 Parciální autokorelační funkce (PACF) Definice 1.5.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 , 45 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.4.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 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 . 46 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 . 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 . 47 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.6.1. 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 jednokrokovou 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.6.4) 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 Obě strany předchozí rovnice vynásobme výrazem Yk+1 - Yk+1 a vypočítejme střední hodnoty. Dostaneme: E Yn+1 Yk+1 - Yk+1 = n-1 j=0 n,n-jE Yj+1 - Yj+1 Yk+1 - Yk+1 48 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 > s 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-jvj = (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, 49 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-jvj = (n + 1, n + 1) - n-1 j=0 2 n,n-jvj čímž jsme dokázali poslední tvrzení (1.11). Poznámka 1.6.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.6.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. 50 1.6.1 Jednokroková nejlepší lineární predikce v AR(p) Nejprve si všimněme, jaké vlastnosti má predikce v případě autoregresních procesů řádu p. Věta 1.6.4. Nechť {Yt, t Z} L2 (, A, P) je centrovaný nedegenerovaný kauzální AR(p) proces Yt = 1Yt-1 + +pYt-p +t. Pak pro nejlepší lineární predikci platí Yn+1 = 0 n = 0 min(n,p) j=1 nYn+1-j n = 1, 2, . . . . Důkaz. Vzhledem k tomu, že autokovarianční funkce (h) exponenciálně klesá k nule, stačí předpokládat, že proces není degenerovaný, tj. rozptyl (0) > 0. Nejlepší lineární predikce podle definice je rovna Yn+1 = Yn+1|n = 0 (= Y ) n = 0, Psp{Y1,...,Yn}(Yn+1) = PMn (Yn+1) n 1. Předpokládejme tedy, že n p a postupně upravujme Yn+1 = Psp{Y1,...,Yn}(Yn+1) = Psp{Y1,...,Yn}(1Yn + + pYn+1-p + n+1) = p j=1 jPsp{Y1,...,Yn}(Yn+1-j) + Psp{Y1,...,Yn}(n+1). Připomeňme, že pro projekci v případě j = 1, . . . , n platí Psp{Y1,...,Yn}(Yj) = Yj, neboť Yj sp{Y1, . . . , Yn}. Dále počítejme pro j = 1, . . . , n skalární součin n+1, Yj kauzal. = n+1, k=1 jj-k = k=1 j E(n+1j-k) =0 = 0, tj. n+1 Y1, . . . , Yn, a Psp{Y1,...,Yn}(n+1) = 0, takže Yn+1 = 1Yn + + pYn+1-p jestliže n p, čímž dostáváme tvrzení věty. Tedy v případě AR(p) procesu Yt = 1Yt-1 + + pYt-p + t 51 jsou koeficienty n,1, . . . , n,1 nejlepší lineární predikce pro n p rovny n,1 = 1 ... n,p = p n,p+1 = 0 ... n,n = 0 1.6.2 Vícekroková nejlepší lineární predikce v AR(p) Podle definice h-kroková 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. Počítejme postupně Yn+2|n = Psp{Y1,...,Yn}(Yn+2) = Psp{Y1,...,Yn}(1Yn+1 + + pYn+2-p + n+2) = p j=1 jPsp{Y1,...,Yn}(Yn+2-j) + Psp{Y1,...,Yn}(n+2) =0(viz předchozí důkaz) = 1Psp{Y1,...,Yn}(Yn+1) + p j=2 j Psp{Y1,...,Yn}(Yn+2-j) =Yn+2-j = 1Yn+1|n + 2Yn + + pYn+2-p ... Yn+p|n = Psp{Y1,...,Yn}(Yn+p) = Psp{Y1,...,Yn}(1Yn+p-1 + + pYn + n+p) = 1Yn+p-1|n + + p-1Yn+1|n + pYn pro s > p Yn+p|n = Psp{Y1,...,Yn}(Yn+p) = Psp{Y1,...,Yn}(1Yn+p-1 + + pYn + n+p) = 1Yn+s-1|n + + pYn+s-p|n 52 1.6.3 PACF pro AR(p), MA(q) a ARMA(p, q) Věta 1.6.5. Nechť {Yt, t Z} L2 (, A, P) je centrovaný nedegenerovaný kauzální AR(p) proces Yt = 1Yt-1 + + pYt-p + t. Pak platí (1) (p) = p (2) (k) = 0 pro k > p. Důkaz. Již dříve jsme ukázali, že v případě AR(p) procesu Yt = 1Yt-1 + + pYt-p + t jsou koeficienty n,1, . . . , n,1 nejlepší lineární predikce pro n p rovny n,1 = 1 ... n,p = p n,p+1 = 0 ... n,n = 0 Tedy pokud přímo n = p, tak podle důsledku Durbin­Lewinsonova algoritmu platí (p) = p,p = p. Jestliže k > p , pak je parciální autokorelační funkce nulová (k) = 0 , což je velmi důležitá identifikační vlastnost AR(p) procesů. Poznámka 1.6.6. Parciální autokorelační koeficienty (1), . . . , (p-1) lze určit jako 1,1, . . . , p-1,p-1 z Durbin­Levinsonova algoritmu. Důsledek 1.6.7. Nechť {Yt, t Z} L2 (, A, P) je centrovaný nedegenerovaný invertibilní MA(q) (resp. ARMA(p, q)) proces. Pak neexistuje takové k0 N, že pro k > k0 platí (k) = 0. Důkaz. Využijeme toho, že proces MA(q) (resp. ARMA(p, q)) je invertibilní. Pak existuje absolutně konvergentní posloupnost reálných čísel = {j} j=0 (tj. j=0 |j| < ) taková, že t = j=0 jYt-j, tj. zkráceně Yt AR() : t = (B)Yt, 53 tj. p = , takže podle předchozí věty nenajdeme k0 N takové, že pro k > k0 platí (k) = 0. 1.6.4 Jednokroková nejlepší lineární predikce v MA(q) Pro jednokrokovou predikci v případě MA(q) procesů je velmi užitečný inovační algoritmus. Nejprve uvedeme podrobně rekurentní vzorce pro stacionární procesy, pro které platí (i, j) = (i - j). Predikci pomocí inovací lze vypočítat pomocí následujícího vzorce Yn+1 = 0 n = 0 n,1 Yn -Yn + n,n Y1 -Y1 n = 1, 2, . . . , přičemž pro n = 0 v0 = (0), dále n,n = (n) v0 n,n-1 = (n-1) v1 - 1,1n,n v0 v1 n,n-2 = (n-2) v2 - 2,2n,n v0 v2 - 2,1n,n-1 v1 v2 ... n,2 = (2) vn-2 - n-2,n-2n,n v0 vn-2 - - n-2,1n,3 vn-3 vn-2 (n-2) členů n,1 = (1) vn-1 - n-1,n-1n,n v0 vn-1 - n-1,n-2n,n-1 v1 vn-1 - - n-1,1n,2 vn-2 vn-1 (n-1) členů a nakonec vn = (0) - 2 n,nv0 - - -2 n,1vn-1. Vzhledem k tomu, že MA(q) proces má autokovarianční funkci (k) = 0 pro k > q, pak v případě, že n > q jsou koeficienty n,n = 0, . . . , n,q+1 = 0 a teprve n,q = 0, . . . , n,1 = 0, takže nejlepší lineární predikce je tvaru Yn+1 = 0 n = 0 min(n,q) j=1 n,j Yn+1-j - Yn+1-j n = 1, 2, . . . 54 1.6.5 Nejlepší lineární predikce v ARMA(p, q) 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. Yt MA(), takže pro |z| 1 dostáváme (z) = (z) (z) (z)(z) = (z). Koeficienty {j} j=0 se určí ze vztahu (1-1z-1z2 - -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 j = 0 pro j > q j > p a označíme-li m = max(p, q) , dosta- neme 0 = 1 j = j + min(j,p) i=1 ij-i pro 1 j m j = p i=1 ij-i pro j > m a vidíme, že pro j > m se koeficienty k neprosadí. Pokud bychom použili predikci pomoci inovací, bude vždy pro n > m platit Yn+1 = n j=1 n,j Yn+1-j - Yn+1-j takže použijeme vždy n předchozích inovací a ztrácíme tak výhodu, která byla u MA procesu konečného řádu. 55 Nechceme-li o tuto výhodu přijít, ukázalo se, že díky jednoduché transformaci vuyžijeme jednak možnosti použít maximálně q předchozích inovací a také toho, že díky AR části je proces lineární kombinací předchozích p hodnot. Položme nejprve m = max(p, q) a definujme Wt = -1 Yt 1 t m -1 (B)Yt = -1 (Yt - 1Yt-1 - - pYt-p) t > m tedy pro 1 t m jde o ARMA(p, q) proces s jednotkovým rozptylem a pro t > m jde o MA(q) proces (opět s jednotkovým rozptylem). Zkoumejme jednokrokovou predikci tohoto transformovaného procesu. Zřejmě Mn = sp{Y1, . . . , Yn} = sp{W1, . . . , Wn}, takže položíme-li W1 = 0 = Y = W pro n = 1, pak pro 1 t m Wt = Psp{Y1,...,Yt-1}(-1 Yt) = -1 Yt a pro t > m Wt = Psp{Y1,...,Yt-1}(-1 (B)Yt) = -1 Psp{Y1,...,Yt-1} (Yt - 1Yt-1 - - pYt-p) = -1 Yt - 1Yt-1 - - pYt-p . Z předchozího také plyne, že Wt - Wt = -1 (Yt - Yt). Použijeme-li inovační algoritmus na proces Wt, dostaneme Wn+1 = n j=1 n,j Wn+1-j - Wn+1-j 1 n m - 1, q j=1 n,j Wn+1-j - Wn+1-j n m. Koeficienty n,j se určí pomocí autokovarianční funkce procesu Wt (viz inovační algoritmus). Zpětnou transformací k původnímu procesu bude nejlepší lineární 56 predikce o jeden krok dopředu rovna Yn+1 = 0 n = 1 n j=1 n,j Yn+1-j - Yn+1-j 1 n m - 1, p j=1 jYn+1-j + q j=1 n,j Yn+1-j - Yn+1-j n m. Při odvození predikce o h > 1 kroků dopředu opět vyjdeme z transformovaného procesu Wt Wn+h|n = Psp{Y1,...,Yn}(-1 Yn+h) n + h m, Psp{Y1,...,Yn} -1 (Yn+h - p j=1 jYn+h-j) n + h > m. 57 1.7 Výstavba modelů v B-J metodologii 1.7.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, tj r(k) = (k) * a parciální autokorelační funkce PACF, tj. 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 neexistuje ko ko = q po q - p hodnotách (k) ve tvaru ,,U je ve tvaru ,,U PACF neexistuje ko neexistuje ko ko = p po p - q hodnotách (k) omezena křivkou ,,U 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. 58 Na následujícím obrázku vidíme simulovaná data kauzálního AR(2) procesu Yt = 0.5Yt-1 + 0.2Yt-2 + t, t N(0, 1), 0 50 100 150 200 250 300 -4 -2 0 2 4 Obrázek 1.8: Simulovaná data pro kauzální AR(2) proces. jehož kořeny polynomu (z) = 1 - 1z - 2z2 = 1 - 0.5z - 0.2z2 jsou rovny z01 = 3.81 a z02 = 1.31. Na následujících dvou grafech jsou vykresleny odhady autokorelační ACF a parciální autokorelační funkce PACF této konkrétní realizace. Odhad ACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Odhad PACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Obrázek 1.9: Odhady ACF a PACF pro simulovaná data z AR(2) procesu. Vidíme, že odhad autokorelační funkce (t) exponenciálně klesá k nule pro t a odhady parciální autokorelační funkce (t) pro t > k0 = 2 se výrazně neliší od nuly. 59 Naproti tomu, pro simulovaná data z invertibilního MA(2) procesu Yt = t - 0.5t-1 - 0.2t-2, t N(0, 1), 0 50 100 150 200 250 300 -4 -2 0 2 4 Obrázek 1.10: Simulovaná data pro invertibilní MA(2) proces. s kořeny polynomu (z) = 1 + 1z + 2z2 = 1 - 0.5z - 0.2z2 rovnými hodnotám z01 = 3.81 a z02 = 1.31, vidíme na odhadech ACF a PACF zobrazených na následujících dvou obrázcích, že naopak odhad parciální autokorelační funkce (t) exponenciálně klesá k nule pro t a odhady autokorelační funkce (t) pro t > k0 = 2 se výrazně neliší od nuly. Odhad ACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Odhad PACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Obrázek 1.11: Odhady ACF a PACF pro simulovaná data z MA(2) procesu. 60 Na závěr si všimneme simulovaných dat z kauzálního a invertibilního ARMA(2, 2) procesu Yt = 0.5Yt-1 + 0.2Yt-2 + t - 0.4t-1 + 0.3t-2, t N(0, 1), 0 50 100 150 200 250 300 -4 -2 0 2 4 6 Obrázek 1.12: Simulovaná data pro invertibilní ARMA(2, 2) proces. jehož všechny kořeny polynomů (z) = 1-1z-2z2 = 1-0.5z-0.2z2 a (z) = 1+1z+2z2 = 1-0.4z+0.3z2 leží vně jednotkové kružnice, neboť AR : z01 = 3.81 a z02 = 1.31 MA : z01 = 1.83 a z02 = 1.83. Podíváme-li se na grafy odhadů autokorelační a parciální autokorelační funkce, vidíme že obě funkce exponenciálně klesají k nule pro t . V tomto případě nelze určit identifikační bod k0, od kterého by se některá z korelačních funkcí již výrazně nelišila od nuly. Odhad ACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Odhad PACF : (t) 0 10 20 30 40 50 -0.5 0 0.5 1 Obrázek 1.13: Odhady ACF a PACF pro simulovaná data z ARMA(2, 2) procesu. 61 1.7.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 - pR-1 p p . Věta 1.7.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 (1991), str. 255-257. 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. 62 Věta 1.7.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 mYm = Psp{Ym,...,Y1}Ym+1, přičemž Ym = (Ym, . . . , Y1) , tj. m = R-1 m m, přičemž Rm = ((i - j))m ij=1. Speciálně pro m > p platí n m,m A N(0, 1). Důkaz. viz Brockwell, Davis (1991), str. 255-257. 1.7.3 Předběžné odhady v AR(p) a Durbin-Levinsův algo- ritmus. 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 - mR-1 m m . Jestliže (0) > 0, pak R1, R2, . . . nejsou singulární a můžeme využít DurbinLevinsův algoritmus pro postupné odhady autoregresních koeficientů 1, 2 a odhady variability bílého šumu v1, v2, . . . . 63 Věta 1.7.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 = b(1) b(0) = (1) v0 = (0) (1.12) m,m = (m) - m-1m-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.7.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.7.4. 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-jm,m-jvj k = 0, . . . , m - 1 vm = (0) - m-1 j=0 2 m,m-jvj 64 Označme m = m,1, . . . , m,m . Pak platí věta Věta 1.7.5. Nechť {Yt, t Z} je kauzální a invertibilní ARMA proces (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 celý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. I když rekurentní odhady koeficientů MA procesů pomocí inovačního algoritmu jsou analogické jako rekurentní odhady koeficientů AR procesů pomocí Durbin-Levinsonova algoritmu, je 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 stabili- zovaly. 1.7.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 65 (z) = (z) (z) (z)(z) = (z). Koeficienty {j} j=0 se určí ze vztahu (1-1z-1z2 - -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 j = 0 j > p 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 do- stáváme 2 = vm a m,j = j + min(j,p) i=1 im,j-i j = 1, 2, . . ., p + q. Nejprve uvažujeme rovnice pro j = q + 1, . . . , p + q 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 . Ř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 im,j-i j = 1, . . . , q. Poznámka 1.7.6. Pro MA(q) platí j = m,j , neboť p = 0. 66 1.7.6 Maximálně věrohodné odhady. Předpokládejme, že {Yt, t Z} je Gaussovský proces s nulovou střední hodnotou a kovarianč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 Y n-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 , pak můžeme psát Yn = (Cn - In)(Yn - Yn) . 67 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 = EYnY n = E Cn(Yn - Yn)(Yn - Yn) C n = CnE (Yn - Yn)(Yn - Yn) C n Nyní počítejme E h (Yn-bYn)(Yn-bYn) i = 0 B B B B B B B B B B B B @ E(Y1-bY1)2 | {z } v0 E(Y1-bY1)(Y2-bY2) | {z } =0 E(Y1-bY1)(Yn-bYn) | {z } =0 E(Y2-bY2)(Y1-bY1) | {z } =0 E(Y2-bY2)2 | {z } v1 E(Y2-bY2)(Yn-bYn) | {z } =0 . .. ... ... . .. E(Yn-bYn)(Y1-bY1) | {z } =0 E(Yn-bYn)(Yn-1-bYn-1) | {z } =0 E(Yn-bYn)2 | {z } vn-1 1 C C C C C C C C C C C C A = diag{v0, . . . , vn-1} = Dn. Takže n = CnDnC n . Počítejme dále Y n-1 n Yn = (Yn - Yn) C n [CnDnC n] -1 Cn(Yn - Yn) = (Yn - Yn) C n (C n) -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| = |CnDnC n| = |Cn| =1 |Dn| |C n| =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 . 68 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 ). Ukazuje se (viz Brockwell, Davis (1991), 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í transformovaný proces Wt = -1 Yt t = 1, . . . , m; m = max(p, q) (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, 69 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 . Poznámka 1.7.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 (1991), §8.9). Takto získané odhady se nazývají odhady metodou nejmenších čtverců. což vede k systé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 ) , 70 platí n MLE - A Nn+p+1(0, V ()), kde V () = 2 EUtU t EUtV t EVtU t 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 (1991), §8.9). 71 1.8 Výstavba modelů a predikce v ARIMA pro- cesech Až dosud jsme uvažovali pouze o procesech (slabě ) stacionárních. V reálných situacích se však se stacionárními procesy setkáváme pouze zřídka. Obecně rozlišujeme dva druhy nestacionarity: * ve střední hodnotě * v rozptylu. 1.8.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 modelová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 72 Nestacionární proces se stochastickým trendem nazýváme integrovaným smíšeným modelem 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, takže (a) hodnoty ACF = (k) budou klesat velmi pomalu (lineárně), (b) hodnoty PACF = (k) jsou logicky velmi podobné procesu AR(1). V praxi se používá následující modifikace: 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 . 73 Poznámka 1.8.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í téměř všech identifikačních detailů těchto funkcí. Poznámka 1.8.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.8.3. ARIMA(p, d, q) nemá smysl centrovat, neboť platí: d (Yt - Y ) = d Yt. Poznámka 1.8.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.8.2 Procesy nestacionární v rozptylu Není-li splněna podmínka neměnnosti rozptylu v čase, je proces nestacionární v rozptylu. Takovýto proces je ovšem třeba nejprve vhodně transformovat. Vysvětleme si stručně pojem transformace stabilizující rozptyl. 74 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 Et g(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 (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 rozptyl, takže rozptyl Dt g(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í. 75 Mocninné transformace Pro přehlednost vynechejme index t a uvažujme kladnou náhodnou veličinu X z rozdělení, které závisí na parametru se střední hodnotou EX = (pokud tomu tak není, provede se vhodná reparametrizace) a rozptylem DX = 2 () = ( )2 , , R, tj. X L(, 2 2 ). Podle obecného vzorce se transformace stabilizující rozptyl vypočítá takto: g() = cd () + K = c d + K = c ln || + K = 1, c 1- 1+ K = 1. . Položme v dalším = 1 - a tento parametr nazvěme transformačním parametrem pro mocninnou trans- formaci. Různou volbou c a K dostaneme následující často užívané transformace * Box-Coxova mocninná transformace pro kladné náhodné veličiny při volbě c = a K = 0 = 0 = 1, - 1 = - 1 1- = 0 = 1, a odtud g(X) = X() = ln X = 0 ( = 1), X -1 = 0 ( = 1). . * Box-Coxova mocninná transformace s posunutím se použije v případě, že hodnoty náhodné veličiny nejsou kladné. Nalezneme proto takové reálné číslo a tak, aby pro všechny realizace platilo x+a > 0 a transformace bude mít tvar: g(X + a) = (X + a)() = ln(X + a) = 0 ( = 1), (X+a) -1 = 0 ( = 1). . * Mocninná transformace se znaménkem lze opět použít v případě, že náhodné veličiny nejsou kladné: g(X) = sign(X)|X|() = sign(X) ln |X| = 0 ( = 1), sign(X)|X| -1 = 0 ( = 1). 76 Odhad transformačního parametru mocninné transformace * Parametrický přístup pomocí metody maximální věrohodnosti. Mějme nezávislé realizace náhodné veličiny X L(X, 2 X 2 X ). Předpokládejme, že existuje takové = 1 - , že transformovaný náhodný vektor Y = (Y1 = g(X1), . . . , Yn = g(Xn)) je výběr z normálního rozdělení se střední hodnotou a rozptylem 2 . Označme y = (y1, . . . , yn) realizaci náhodného výběru. Hledejme maximum věrohodnostní funkce pro = (, 2 ) , tj. pro funkci L(, 2 ) = n i=1 - 1 22 exp 1 2 yi - 2 = (22 )- n 2 exp - 1 2 n i=1 yi - 2 , což je stejná úloha jako hledat maximum logaritmu věrohodnostní funkce l(, 2 ) = ln L(, 2 ) = - n 2 ln(2) - n 2 ln(2 ) - 1 2 n i=1 yi - 2 . Maxima nalezneme, položíme-li l = 0 a l 2 = 0. 0 = l = 2 22 n i=1 (yi - ) 0 = l 2 = - n 22 + 1 24 n i=1 (yi - )2 a odtud pak dostaneme ^ = 1 n n i=1 yi = y, 2 = 1 n n i=1 (yi - y)2 = s2 . Upravme nyní logaritmus věrohodnostní funkce takto: l(, 2 ) = - n 2 ln(2) - n 2 ln(2 ) - 1 22 n i=1 [(yi - y) + (y - )] 2 = - n 2 ln(2) - n 2 ln(2 ) - 1 22 n i=1 (yi - y)2 + n(y - )2 = - n 2 ln(2) - n 2 ln(2 ) - 1 22 ns2 + n(y - )2 77 Nyní dokažme, že funkce l(, 2 ) nabývá v bodě (^, ^2 ) = (y, s2 ) svého maxima. Platí l(y, s2 ) = - n 2 ln(2) - n 2 ln(s2 ) - n 2 , Ověřme, zda platí nerovnost l(, 2 ) ? l(y, s2 ) -n 2 ln(2) - n 2 ln(2 ) - ns2 +n(y-)2 22 ? -n 2 ln(2) - n 2 ln(s2 ) - n 2 -1 2 ln(2 ) - s2 22 - (y-)2 22 ? -n 2 ln(s2 ) - 1 2 0 ? s2 22 - 1 2 -ln s 1. člen + (y-)2 22 0 Protože pro všechna kladná x = s > 0 platí ln x < x2 -1 2 , je první i druhý člen nezáporný a nerovnost platí. Celkově jsme tedy dostali, že max ,2 l(, 2 ) = l(y, s2 ) = - n 2 ln(2) - n 2 ln(s2 ) - n 2 a max ,2 L(, 2 ) = L(y, s2 ) = 2s2 - n 2 e- n 2 . Nyní toto maximum vyjádřeme v původních proměnných xi, kdy yi = g(xi) = ln xi = 0, x i -1 = 0. Nejprve vypočtěme jakobián této transformace: |J| = n i=1 dyi dxi = n i=1 x-1 i = n i=1 x-1 i . Pak max ,2 L(, 2 , ) = 2s2 () - n 2 e- n 2 |J| = 2s2 () - n 2 e- n 2 n i=1 x-1 i = 2s2 () - n 2 e- n 2 n i=1 e(-1) ln xi max ,2 l(, 2 , ) = - n 2 ln(2) - n 2 ln(s2 ()) - n 2 + ( - 1) n i=1 ln xi. 78 Nyní hledejme maximum funkce l(^, ^2 , ) = l(y, s2 , ) pro parametr . Protože maximum vzhledem k nezávisí na konstantách, budeme maximalizovat funkci l () = - n 2 ln(s2 ()) + ( - 1) n i=1 ln xi. Teoretickým odvozením maximálně věrohodného odhadu parametru se zde již dále nebudeme zabývat, ale ukážeme si jednodušší přístup, který pro ekvidistantní hodnoty 1 < 2 < < m (pro dostatečně velké m) ze vhodně zvoleného intervalu ( 1, 2), (kde 1, 2 R, 1 < 2) vypočítá hodnoty l () a hledá argument ^ maxima těchto hodnot. Ve své práci Box a Cox (1964) odvodili asymptotické rozdělení statistiky K = -2 l () - l (^) A 2 (1), takže můžeme zkonstruovat jednostranný asymptotický interval spolehlivosti pro parametr 1 - = P K < 2 1-(1) = P -2 l () - l (^) < 2 1-(1) = P l (^) - 1 2 2 1-(1) =D l () , tj. všechna splňující nerovnost l () D leží v intervalu spolehlivosti a jsou tedy přijatelná. Testování hypotéz typu H0 : = 0 proti alternativě H1 : > 0: 1. Budeme testovat hypotézu H1 0 : = 1. Pokud hypotézu nezamítneme, tj. l (1) D, nemusíme data transformovat. 2. Pokud předchozí hypotézu zamítneme, můžeme testovat další hypotézu H2 0 : = 0. Pokud H2 0 nezamítneme, tj. l (0) D l (1) < D, transformace bude tvaru yi = ln xi. Pokud však se l (0) < D l (1) < D, provedeme transformaci yi = x ^ i - 1 ^ . 79 * Jednoduchý algoritmus v praktických úlohách 1. Algoritmus nejprve zkontroluje vstupní data tak, aby byla nezáporná, tj. případně přičte kladnou konstantu. 2. Upravený vektor dat rozdělí na krátké úseky o délce 4 až 12 údajů. 3. V každém úseku dat se provede pokud možno robustní odhad střední hodnoty ^i (průměr, medián) a robustní odhad variability ^2 i (např. max-min, interkvartilové rozpětí). 4. Protože předpokládáme, že platí () = pak logaritmovanáním dostaneme vztah ln(()) = ln a + ln(), takže neznámé můžeme odhadnout pomocí metody nejmenších čtverců díky hodnotám zi = ln(i) a ui = ln ^i v regresním modelu zi = a + ui + i i WN(0, 2 ). 5. Pro odhad = 1 - ^ pomocí t-statistiky zkonstruujeme interval spolehlivosti I(). ­ Pokud tento interval bude obsahovat nulu, tj. 0 I() data se nebudou transformovat. ­ Pokud 0 / I() 1 I(), volí se logaritmická transformace yi = ln xi. ­ Jinak se volí mocninná transformace yi = x ^ i - 1 ^ . 80 International airline passenger monthly totals (in thousands), Jan.49-Dec.60 20 40 60 80 100 120 140 150 200 250 300 350 400 450 500 550 600 Observed data 20 40 60 80 100 120 140 7 7.5 8 8.5 9 9.5 10 10.5 Data transformed by means of Box-Cox 0 100 200 300 400 500 600 0 2 4 6 8 10 12 14 16 18 20 Histogram of observed data 6 7 8 9 10 11 0 2 4 6 8 10 12 14 16 18 Histogram of transformed data (Box-Cox) -3 -2 -1 0 1 2 3 -800 -780 -760 -740 -720 -700 -680 Loglikelihood function max = 0.15 max L = -0.24 U = 0.53 10 2 log(median) log(IQR) 1 =1.215; =1-1 =-0.2152 0.4925 Obrázek 1.14: Ukázka mocninné transformace. Protože 0 leží uvnitř intervalu spolehlivosti pro parametr , lze doporučit logaritmickou transformaci. 81 1.8.3 Volba řádu modelu Vhodná volba řádu modelu je důležitou součástí výstavby modelů v Box-Jenkinsonově metodologii. 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) EYpY p 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 ). 82 Odtud dostaneme, že n 2 (p - p) p(p - p) = n 2 Zp A 2 (p). 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í nc2 2 A 2 (n - p) a protože E nc2 2 = n - p, pak E2 = n-p n 2 2 = n n-p E2 . Pokud ve výpočtu vn nahradíme 2 jeho nestranným odhadem 2 , dostaneme následující krité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 n2 n - p - q + (p + q) ln P Y 2 t -nc2 p+q . 83 1.8.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.8.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 b2 W (j) n-j A 2 (h - p - q) pro h N, nebo (Granger, Anderson, 1978) ~QWW = n(n + 2) h j=1 bW W (j) n-j A 2 (h) pro h N. 84 1.9 Modelování sezónnosti pomocí SARIMA mo- delů Sezónnost je v Box-Jenkinsonově metodologii stejně jako trend modelována sto- chasticky. 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) 85 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 ) 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.9.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.9.2. 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é. 86 1. Odhad parametrů d, D : (a) Provede se studium odhadnuté autokorelační funkce ACF = (k), neboť identifikuje 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é hodnoty 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 diferencovat (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í. (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 bo- dech 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 87 * 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í amplitudou) 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 sinusoidě 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, . . . * 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. 88 (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. homogenní 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í. 89 Kapitola 2 Dynamické lineární modely 2.1 Motivační příklad (1) Představme si, že jsme se ztratili na moři, popřípadě v podmínkách vnitrozemského státu v hlubokém lese. Nacházíme se v nějaké pozici x , která objektivně existuje, ale kterou neznáme a snažíme se ji určit. Skutečná pozice x je neznámým stavem. (2) Pomocí hvězd odhadneme naši pozici jako y0 s určitou nepřesností, neurčitostí 2 y0 . Na základě dosavadních znalostí známe hustotu pozice x , kterou označíme p1|0(x) = p(x|y0) , se střední hodnotou ^x1|0 = y0 a rozptylem 2 1|0 = 2 y0 . (3) Mraky se rozestoupí a jasněji vidíme hvězdy. Proto upřesníme náš odhad jako y1 s větší jistotou 2 y1 < 2 y0 . Na základě širších znalostí dostaneme novou hustotu pozice x jako p1|1(x) = p(x|y0, y1) se střední hodnotou ^x1|1 a rozptylem 2 1|1. Střední hodnotu vypočteme jako vážený průměr pozorování y0 a y1, kde váha je o to větší, oč je pozorování přesnější, tj. má menší rozptyl (přitom součet vah je roven 1) ^x1|1 = 1 2 y0 1 2 y0 + 1 2 y1 y0 + 1 2 y1 1 2 y0 + 1 2 y1 y1 = 2 y1 2 y0 +2 y1 =(1-K) y0 + 2 y0 2 y0 +2 y1 =K y1 = (1-K)y0 +Ky1 Za předpokladu, že y0 a y1 jsou nezávislá pozorování, pak rozptyl váženého 90 průměru je roven 2 1|1 = 2 y1 2 y0 +2 y1 2 2 y0 + 2 y0 2 y0 +2 y1 2 2 y1 = 2 y0 2 y0 +2 y1 2 y1 = K2 y1 = (1 - K)2 y0 tj. 2 1|1 < 2 y1 Celkově dostaneme tyto rekurentní vztahy ^x1|1 = ^x1|0 + K(y1 - ^x1|0) 2 1|1 = (1 - K)2 1|0, tj. využili jsme informací z prvního i druhého měření. | | | ^x1|0 = y0 ^x1|1 y1 (4) Nyní zdynamizujeme (rozpohybujeme) náš příklad (a) Zatímco děláme svoje měření, nestojíme na místě, ale pohybujeme se, tj. existuje jednoduchý dynamický model dx dt = v konstantní posun + w náhodná složka . Předtím, než provedeme další měření, uděláme predikci ^x2|1 = F1 ^x1|1 na základě předchozího stavu ^x1|1 a dynamického modelu (nějaké funkce přechodu F1) s určitou dávkou neurčitosti 2 2|1 (b) Provedeme pomocí hvězd další měření polohy y2 s neurčitostí, nepřesností 2 y2 . (5) Všechny předchozí informace shrneme do odhadu polohy ^x2|2 = ^x2|1 + K(y2 - ^x2|1) s vahou K = 2 2|1 2 2|1 + 2 y2 (též tzv. Kalmanův zisk) 2 2|2 = (1 - K)2 2|1. 91 | | | ^x1|0 = y0 ^x1|1 y1 | | | ^x2|1 = F1 ^x1|1^x2|2 y2 pohyb 2.2 Stavově-prostorové modely Místo jednorozměrné náhodné posloupnosti {Yt, t Z} uvažujme posloupnost w-rozměrných náhodných vektorů {Yt, t Z}, Yt Rw , které splňují tzv. Datové a stavové rovnice DATOVÁ ROVNICE Yt = GtXt + Wt t = 1, 2, 3, . . . STAVOVÁ ROVNICE Xt+1 = FtXt + Vt t = 1, 2, 3, . . . přičemž Xt . . . je tzv. stavový v-rozměrný náhodný vektor Wt . . . je šum měření Vt . . . je šum procesu Gt . . . je posloupnost matic typu w × v (popisují vztah pozorování ke stavu) Ft . . . je posloupnost matic typu v × v (modelují dynamiku - tzv. matice přechodu) Dále platí EVt = 0 EWt = 0 D Wt Vt = Rt St S t Qt tj. EWtW t = Rt EVtV t = Qt EWtV t = St C(Xt, (W t, V t) ) = 0, tj. jsou nekorelované Všechny náhodné vektory mají konečné druhé momenty. 92 Příklad 2.2.1. NÁHODNÁ PROCHÁZKA S TRENDEM Mějme R, šum procesu Vt WN(0, 2 v), náhodné veličiny T rt, přičemž T r0 = 0 = 0. Dále nechť pro t = 1, 2, . . . platí C(T rt, Vt) = 0 tj. T rt a Vt jsou nekorelované, což značíme T rt Vt. Definujme T rt+1 = T rt + + Vt = T rt-1 + + Vt-1 + + Vt = T rt-1 + 2 + Vt + Vt-1 = po t krocích = T r0 =0=0 +t + t j=1 Vj Položme Xt = T rt Vt = Vt 0 Ft = 1 1 0 1 . Pak Xt+1 = T rt+1 = 1 1 0 1 T rt + Vt 0 = FtXt + Vt t = 1, 2, . . . . Označme šum měření Wt WN(0, 2 w) a položme Yt = 1 0 =Gt T rt + Wt = GtXt + Wt t = 1, 2, . . . . Jestliže X1 = T r1 , V1, W1, V2, W2, . . . jsou nekorelované, dostáváme stavověprostorovou reprezentaci náhodné procházky, pro kterou platí EVt = 0 DVt = EVtV t = 2 v 0 0 0 = Qt = Q EWt = EWt = 0 DWt = EWtW t = EW2 t = 2 w = Rt = R EVtW t = 0 0 = St = S. 93 Příklad 2.2.2. SEZÓNNÍ ŘADA SE ŠUMEM Uvažujme sezónu délky d a sezónní komponenty s1, . . . , sd přičemž platí st+d = st a s1 + + sd = 0. Takže dostaneme st+1 = st+1-d ... st+1 + st + st-1 + + st+1-d+1 = 0 Odtud získáme deterministickou rovnici st+1 = -st - st-1 - - st+2-d Přidejme šum procesu Vt WN(0, 2 v) a dostaneme po přeznačení stochastickou rovnici Yt+1 = -Yt - Yt-1 - - Yt+2-d + Vt . Položme Xt+1 = Yt+1 Yt Yt-1 ... Yt+3-d = -1 -1 -1 -1 1 0 0 0 0 ... ... ... ... ... ... ... 0 0 0 0 1 0 =Ft Yt Yt-1 Yt-2 ... Yt+2-d =Xt + Vt 0 0 ... 0 =Vt tj. stavově-prostorový model sezónní řady se šumem je roven Xt+1 = FXt + Vt Yt = 1 0 0 =Gt Xt 94 2.3 Stacionární stavově-prostorové modely DATOVÁ ROVNICE Yt = GXt + Wt t = 1, 2, 3, . . . STAVOVÁ ROVNICE Xt+1 = FXt + Vt t = 1, 2, 3, . . . přičemž Xt . . . je tzv. stavový v-rozměrný náhodný vektor Wt . . . je šum měření Vt . . . je šum procesu G . . . je matice typu w × v (popisují vztah pozorování ke stavu) F . . . je matice typu v × v tzv. matice přechodu Dále platí EVt = 0 EWt = 0 D Wt Vt = Rt St S t Qt tj. EWtW t = Rt EVtV t = Qt EWtV t = St C(Xt, (W t, V t) = 0, tj. jsou nekorelované Všechny náhodné vektory mají konečné druhé momenty. Stavová rovnice se nazývá stabilní (také kauzální), právě když všechna vlastní čísla matice F leží uvnitř jednotkové kružnice, tj. det(I - Fz) = 0 pro |z| < 1. Pokud je systém stabilní (kauzální), pak Xt+1 = j=0 Fj Vt-j Yt = Wt + G j=0 Fj Vt-1-j 95 Příklad 2.3.1. AUTOREGRESNÍ MODEL ŘÁDU p AR(p) : Yt = 1Yt-1 + pYt-p + t , kde t WN(0, 2 ). Yt+1 Yt Yt-1 ... Yt+1-p =Xt+1 = 1 2 p-1 p 1 0 0 0 0 ... ... ... ... ... ... ... 0 0 0 0 1 0 =F Yt Yt-1 Yt-2 ... Yt-p =Xt + t 0 0 ... 0 =Vt Yt = 1 0 0 =G Xt Příklad 2.3.2. MA PROCES ŘÁDU q MA(q) : Yt = t + 1t-1 + + qt-q , kde t WN(0, 2 ). t t-1 t-2 ... t+1-q =Xt+1 = 0 0 0 0 1 0 0 0 0 ... ... ... ... ... ... ... 0 0 0 0 1 0 =F t-1 t-2 t-3 ... t-q =Xt + t 0 0 ... 0 =Vt Yt = 1 q =G Xt + t =Wt 96 Příklad 2.3.3. ARMA PROCES ŘÁDU p, q ARMA(p, q) : Yt = 1Yt-1 + pYt-p + t + 1t-1 + + qt-q , kde t WN(0, 2 ). Yt Yt-1 Yt-2 ... ... Yt+1-p t+1 t ... t+1-q =Xt+1 = 1 p-1 p 1 1 q-1 q 1 0 0 0 ... ... ... ... ... ... ... ... ... ... 1 ... ... ... ... 0 ... ... ... ... 1 ... ... ... ... ... ... ... 0 0 1 0 =F Yt-1 Yt-2 Yt-3 ... ... Yt-p t t-1 ... t-q =Xt + 0 0 0 ... ... 0 t+1 0 ... 0 =Vt Yt = 1 p 1 1 q =G Xt 2.4 Nejlepší lineární predikce pomocí projekce náhodných vektorů druhého řádu Mějme pravděpodobnostní prostor (, A, P). Pro pevně zvolené v N označme Lv 2 = {X = (X1, . . . , Xv) : X1 L2 (, A, P), . . . , Xv L2 (, A, P)} a označme L 2 = v=1 Lv 2. Pak lze nad tímto prostorem definovat skalární součin pro X Lv 2 a Y Lw 2 (v, w N) předpisem X, Y = EXY za předpokladu, že existuje sdružené rozdělení náhodného vektoru Z = X Y a Z Lv+w 2 . 97 Označme pro Y0, . . . , Yt Lw 2 Mt = sp{Y0, . . . , Yt} uzavřený podprostor generovaný všemi možnými lineárními kombinacemi typu C0Y0 + + CtYt, kde C0, . . . , Ct jsou reálné matice. Pak uvažujme nad L 2 projekci X Lv 2 do Mt PMt (X) = (PMt (X1), . . . , PMt (Xv)) , kterou budeme značit různými způsoby, a to PMt (X) = ^X = Pt(X) = Pt(X|Y0, . . . , Yt). Připomeňme vlastnosti predikce, které v následujících důkazech využijeme (a) vždy existuje jediný vektor Pt(X) takový, že pro Y Mt = sp{Y0, . . . , Yt} platí X - ^X, Y = 0 X, Y = ^X, Y EXY = E ^XY . Protože ^X Mt, pak EX ^X = E ^X ^X . (b) Jestliže X, Y1, . . . , Yt mají sdružené normální rozdělení, pak (pokud Y0 = 1 = (1, . . . , 1) ) platí ^X = Pt(X) = E(X|Y1, . . . , Yt) t 1. (c) Predikce ^X = Pt(X) je lineární v tom smyslu, že pro libovolnou matici A Rk+v a X, Z Lv 2 platí: Pt(AX) = APt(X) Pt(X + Z) = Pt(X) + Pt(Z) (d) Pokud X Lv 2 a Y Lw 2 ,pak Pt(X|Y) = MY, kde M Rv+w , pro níž platí M = EXY [EYY ]a Aznačí pseudoinverzní matici k matici A, což je taková matice, pro níž A= AA- A Každá matice má alespoň jednu pseudoinverzní matici. Pokud matice A je regulární, pak A= A-1 . 98 Připomeňme opět definici stavového modelu DATOVÁ ROVNICE Yt = GtXt + Wt t = 1, 2, 3, . . . STAVOVÁ ROVNICE Xt+1 = FtXt + Vt t = 1, 2, 3, . . . přičemž Xt . . . je tzv. stavový v-rozměrný náhodný vektor Wt . . . je šum měření Vt . . . je šum procesu Gt . . . je posloupnost matic typu w × v Ft . . . je posloupnost matic typu v × v Dále platí EVt = 0 EWt = 0 D Wt Vt = E Wt Vt W t V t = EWtW t EWtV t EVtW t EVtV t = Rt St S t Qt C(Xt, (W t, V t) = 0, tj. jsou nekorelované: Xt (W t, V t) Označme ^Xt|k = Psp{Y0,...,Yk}(Xt) = P(Xt|Y0, . . . , Yk) t|k = E(Xt - ^Xt|k)(Xt - ^Xt|k) Pokud k = t - 1 . . . jde o tzv. problém (jednokrokové) predikce k = t . . . jde o tzv. problém filtrace k = n > t. . . jde o tzv. problém vyhlazení Přidejme předpoklady pro t Wt {Y0, . . . , Yt-1}, tj. jsou nekorelované Vt {Y0, . . . , Yt-1} St = 0 (tj. šumy procesu Vt a měření Wt jsou nekorelované) 99 Věta 2.4.1. Jednokroková Kalmanova predikce Xt = Xt|t-1 = Pt-1(Xt) = P(Xt|Y0, . . . , Yt-1) = Psp{Y0,...,Yt-1}(Xt) a chybová predikční kovarianční matice t|t-1 = E(Xt - Xt)(Xt - Xt) = E(Xt - Xt|t-1)(Xt - Xt|t-1) jsou jednoznačně určeny (1) počátečními podmínkami: X1 = X1|0 = P(X1|Y0) = Psp{Y0}(X1) 1|0 = 1 = X1X1 - bX1 bX1 kde X1X1 = EX1X 1 bX1 bX1 = EX1X 1 (2) a platí pro ně následující rekurentní vztahy: Xt+1 = Xt+1|t = FtXt|t-1 + Kt+1|t(Yt - GtXt|t-1), kde Kt+1|t je tzv. predikční Kalmanův zisk, pro nějž platí: Kt+1|t = Xt+1It - ItIt přičemž Xt+1It = Ftt|t-1G t ItIt = Gtt|t-1G t + Rt a t+1|t = Xt+1Xt+1 - bXt+1 bXt+1 přičemž Xt+1Xt+1 = FtXtXt F t + Qt bXt+1 bXt+1 = FtbXt bXt F t + Kt+1|tItIt K t+1|t a kde It jsou inovace pro Yt, tj. It = Yt - Yt = Yt - Psp{Y0,...,Yt-1}(Yt) 100 Důkaz. Nejprve definujme inovaci pro Yt I0 = Y0 It = Yt - Yt = Yt - Psp{Y0,...,Yt-1}(Yt) = Yt - Pt-1(Yt) = Yt - Pt-1(GtXt + Wt) = Yt - GtPt-1(Xt) - Pt-1(Wt) Díky nezávislosti náhodných vektorů Wt {Y0, . . . , Yt-1} platí Pt-1(Wt) = Psp{Y0,...,Yt-1}(Wt) = 0, takže dostaneme It = Yt - GtXt = GtXt + Wt - GtXt = Gt Xt - Xt + Wt . Je třeba si uvědomit, že inovace jsou ortogonální (tj. nekorelované) I0 I1 I2 . . . It, takže pro libovolné X platí Pt(X) = P(X|Y0, . . . , Yt) = P(X|I0, . . . , It) = P(X|I0, . . . , It-1) + P(X|It) = Pt-1(X) + P(X|It) = Pt-1(X) + MIt, kde M = EXI t[EItI t]- . Takže Xt+1 = Xt+1|t = Pt(Xt+1) = Pt-1(Xt+1) + P(Xt+1|It) = Pt-1(FtXt + Vt) + EXt+1I t[EItI t]- It a označíme-li Xt+1It = EXt+1I t ItIt = EItI t, 101 pak Xt+1 = Ft Pt-1(Xt) = bXt|t-1 + Pt-1(Vt) =0 +Xt+1It - ItIt (Yt - GtXt) Vyjádřeme nyní Xt+1It = EXt+1I t = E(FtXt + Vt) Gt(Xt - Xt) + Wt = E Ft(Xt -Xt)+FtXt +Vt Gt(Xt -Xt)+Wt = Ft E(Xt - Xt)(Xt - Xt) =t|t-1 G t + Ft E(Xt -Xt)W t =0(nekorel.) +Ft EXt(Xt - Xt) =0(nekorel.) + Ft EXtW t =0(nekorel.) + EVt(Xt - Xt) =0(nekorel.) G t + EVtW t =St=0 = Ftt|t-1G t Dále počítejme ItIt = EItI t = E Gt(Xt - Xt) + Wt Gt(Xt - Xt) + Wt = GtE(Xt - Xt)(Xt - Xt) G t + Gt E(Xt - Xt)W t =0(nekorel.) + EWt(Xt - Xt) =0(nekorel.) G t + EWtW t = Gtt|t-1G t + Rt Tedy celkově máme Xt+1 = Xt+1|t = FtXt|t-1 + Xt+1It - ItIt (Yt - GtXt) =It a Kt+1|t = Xt+1It - ItIt je tzv. Kalmanův predikční zisk a můžeme tedy psát Xt+1 = FtXt|t-1 + Kt+1|t(Yt - GtXt|t-1) Zbývá najít rekurentní vztah pro t+1|t. Přitom využijeme důležitý vztah, který vychází z vlastností ortogonální projekce, tj. že pro Y Psp{Y0,...,Yt-1} platí Xt - Xt, Y = 0 Xt, Y = Xt, Y EXtY = EXtY 102 a protože Xt Psp{Y0,...,Yt-1}, dostaneme EXtXt = EXtXt . Proto počítejme 1|0 = 1 = E(X1 - X1)(X1 - X1) = EX1X 1 - EX1X 1 =E bX1 bX 1 - EX1X1 =E bX1 bX 1 +EX1X 1 = EX1X 1 - EX1X 1 = X1X1 - bX1 bX1 Úplnou matematickou indukcí obdobně dokážeme, že pokud budeme předpokládat, že platí t|t-1 = XtXt - bXt bXt , pak t+1|t = E(Xt+1 - Xt+1)(Xt+1 - Xt+1) = EXt+1X t+1 - EXt+1X t+1 =E bXt+1 bX t+1 - EXt+1X t+1 =E bXt+1 bX t+1 +EXt+1X t+1 = EXt+1X t+1 - EXt+1X t+1 = Xt+1Xt+1 - bXt+1 bXt+1 Xt+1Xt+1 = EXt+1X t+1 = E(FtXt) + Wt)(FtXt) + Wt) = Ft EXtXt =XtXt F t + Ft EXtW t =0(nekorel.) + EWtX t =0(nekorel.) F tFtEVtV t = FtXtXt F t + Qt bXt+1 bXt+1 = EXt+1X t+1 = E FtXt + Kt+1|t(Yt - GtXt) FtXt + Kt+1|t(Yt - GtXt) = Ft EXtX t bXt bXt F t + Ft EXt =It (Yt - GtXt) =0(nekorel.) K t+1|t + Kt+1|t E(Yt -GtXt)Xt =0(nekorel.) F t + Kt+1|t E(Yt -GtXt)(Yt -GtXt) =ItIt K t+1|t = FtXtXt F t + Kt+1|tItIt K t+1|t 103 Věta 2.4.2. Pro Kalmanovu filtraci Xt|t = P(Xt|Y0, . . . , Yt) a filtrovací chybovou kovarianční matici t|t = E(Xt - Xt|t)(Xt - Xt|t) pro t 1 platí následující rekurentní vztahy Xt|t = Xt|t-1 + Kt|t(Yt - GtXt|t-1) , kde Kt|t je tzv. filtrační Kalmanův zisk, pro nějž platí Kt|t = XtIt - ItIt , přičemž XtIt = t|t-1G t ItIt = Gtt|t-1G t + Rt a t|t = (I - Kt|tGt)t|t-1 , kde I je jednotková matice řádu v × v. Důkaz. Využijme opět inovací I0 = Y0 It = Yt - Yt = Yt - Psp{Y0,...,Yt-1}(Yt) = Yt - Pt-1(Yt) = Yt - Pt-1(GtXt + Wt) = Yt - GtPt-1(Xt) - Pt-1(Wt = Yt - GtXt = Gt(Xt - Xt) + Wt které jsou navzájem kolmé (tj. nekolerované). 104 Počítejme Xt|t = P(Xt|Y0, . . . , Yt) = P(Xt|I0, . . . , It) = P(Xt|I0, . . . , It-1) + P(Xt|It) = Xt + MIt = Xt + EXtI t [EXtI t] - It = Xt + XtIt - ItIt (Yt - GtXt) přičemž XtIt = EXtI t = EXt Gt(Xt - Xt) + Wt = E (Xt - Xt) + Xt Gt(Xt - Xt) + Wt = E(Xt - Xt)(Xt - Xt) Gt + E(Xt - Xt)W t =0(nekorel.) + EXt(Xt - Xt) =0(nekorel.) G t + EXtWt =0(nekorel.) = t|t-1G t Takže celkově dostaneme Xt|t = Xt + XtIt - ItIt ozn.Kt|t (Yt - GtXt) =It , odtud Xt|t - Xt = Kt|tI. Zbývá dopočítat t|t. Víme, že t|t-1 = E(Xt - Xt)(Xt - Xt) = E (Xt - Xt|t) + (Xt|t - Xt|t-1) =Kt|tIt (Xt - Xt|t) + (Xt|t - Xt|t-1) =Kt|tIt = E(Xt - Xt|t)(Xt - Xt|t) t|t + E(Xt - Xt|t)I t =0(nekorel.) K t|t + Kt|t EIt(Xt - Xt|t) =0(nekorel.) +Kt|t EItI t K t|t = t|t + Kt|tItIt K t|t 105 Protože Kt|t = XtIt - ItIt = t|t-1G t- ItIt dostáváme t|t-1 = t|t + XtIt - ItIt ItIt - ItIt XtIt = t|t + XtIt - ItIt XtIt = t|t + t|t-1G t- ItIt =Kt|t Gtt|t-1 Odtud t|t = t|t-1 - Kt|tGtt|t-1 = (I - Kt|tGt)t|t-1. 2.5 Kalmanův iterační proces Shrňme předchozí výsledky Kalmanovy predikce a filtrace takto 1) Protože Kt+1|t = Xt+1It - ItIt = Ftt|t-1G t- ItIt Kt|t = XtIt - ItIt = t|t-1G t- ItIt Kt+1|t = FtKt|t . 2) Dále platí FtXt|t = FtXt|t-1 + FtKt|t(Yt - GtXt) = FtXt|t-1 + Kt+1|t(Yt - GtXt) = Xt+1|t . 3) Budeme se snažit nově vyjádřit t+1|t . Protože Xt+1 - Xt+1|t = FtXt + Vt - FtXt|t = Ft(Xt - Xt|t) + Vt a vzhledem k tomu, že Vt (Xt - Xt|t), dostáváme t+1|t = E(Xt+1 - Xt+1|t)(Xt+1 - Xt+1|t) = E Ft(Xt - Xt|t) + Vt Ft(Xt - Xt|t) + Vt = FtE(Xt - Xt|t)(Xt - Xt|t) F t + EVtV t = Ftt|tF t + Qt . Všechny předchozí mezivýsledky použijeme pro odvození velmi jednoduchého Kalmanova iteračního procesu, který je spojením filtrace a predikce. 106 Kalmanův iterační proces (I) Počáteční podmínky X1|0 = X1 = EX1 při Y0 = 1 1|0 = E(X1 - EX1)(X1 - EX1) = DX1 (II) Datový (filtrační) krok Kalmanova filtru Nejprve se spočítá tzv. Kalmanův zisk (nebo též Kalmanovo zesílení) Kt|t = t|t-1G t(Gtt|t-1G t + Rt)- , pak Xt|t = Xt|t-1 + Kt|t(Yt - GtXt|t-1) a filtrační chybovou kovarianční matici t|t = (I - Kt|tGt)t|t-1. (III) Časový (predikční) krok Kalmanova filtru Xt+1|t = FtXt|t t+1|t = Ftt|tF t + Qt. t - 2 t - 1 t t + 1 Xt-2|t-2 Ft-2 matice přechodu Xt-1|t-2 (predikce) filtrace Xt-1|t-1 Ft-1 matice přechodu Xt|t-1 (predikce) filtrace ^Xt|t Ft matice přechodu Xt+1|t Yt-1 Yt 107 Literatura [1] ANDĚL, J. Matematická statistika. SNTL/ALFA Praha, 1978. [2] ANDĚL, J. Statistické metody. Matfyzpress Praha, 1993. [3] ANDĚL, J. Statistická analýza časových řad. Praha. SNTL 1976. [4] BOX, G.E.P, COX, D.R.: Analysis of Transformations. Journals of the Royal Statistical Society, Biometrika 26, 1964, 211-252. [5] BOX, G., JENKINS, G. Time series analysis - forecasting and control. Holden-Day 1976. [6] BROCKWELL, P.J., DAVIS, R.A. Time Series: Theory and Methods. Springer-Verlag, New York, 1991. [7] BROCKWELL, P.J., DAVIS, R.A. Introduction to time series and forecasting. Springer-Verlag, New York, 2002. [8] CIPRA, T. Analýza časových řad s aplikacemi v ekonomii. SNTL, Praha, 1986. [9] GICHMAN, I.I., SKOROCHOD, A.V. Teorija slučajnych processov. Moskva. Nauka 1971. [10] GRANGER, C.W.J., ANDERSEN, A. An introduction to bilinear time series models . Vandenhoeck and Ruprecht, Göttingen 1978. [11] HAMILTON, J.D. Time Series Analysis. Princeton University Press. 1994. [12] KALMAN, R. A new approach to linear filtering and prediction problems. Trans. ASME J. Basic Eng. D 82 (1960), 34-45. [13] KUBÁČKOVÁ, L., KUBÁČEK, L., KUKUČA, J. Pravdepodobnostť a štatistika v geodézii a geofyzike. Veda, Bratislava (1982). [14] LJUNG, G. M., BOX, G. E. P. On a measure of lack of fit in time series models. Biometrika 65. 1978, 553­564. [15] NEUBRUNN, T., RIEČAN, B. Miera a integrál. Bratislava. Veda 1981. 108 [16] PRIESTLEY, M. Spectral analysis and time series. Academic Press 1989. [17] RAO, R.C. Lineární metody statistické indukce a jejich aplikace. ACADEMIA Praha, 1978.