Stochastické modely časových řad (9. cvičení) 1 M5201 - 9. CVIČENÍ: Exponenciální vyrovnávání 1 Exponenciální vyrovnávání Hlavní myšlenka lokální (vážené) metody nejmenších čtverců spočívá v tom, že provedeme odhad trendu polynomem na lokálním intervalu (t — s, t + s) na rozdíl od klasické (vážené) metody nejmenších čtverců, kdy trend odhadujeme polynomem na celém intervalu možných hodnot parametru t. Exponenciální vyrovnávání aproximuje neznámý trend polynomem stupně m, který se odhaduje na základě pozorování z asymetrického vyhlazovacího okénka směrem do minulosti. Exponenciální vyrovnávání vychází z polynomiální lokální vážené metody nejmenších čtverců, kde váhy jednotlivých čtverců se směrem do minulosti exponenciálně snižují - odtud název metody. Uvažujeme-li vyhlazovací okno pouze směrem do minulosti, pak pro každé í,r = 0,1,... dostaneme regresní model tvaru Yt_T=j2(-Tyi33(t)+st- 3=0 Est-T=0; Esq£s = 0, q ^ s; Ľ£í_r = a Ta2; a £ (0,1). tj. matice vah je rovna W = diag{wi,..., wn,...} = diag{a°, a1,..., aT ...}. Odhad parametrů f3 provedeme metodou nejmenších vážených čtverců (neboť rozptyly nejsou konstantní) která je dána vzorcem: -lv'i P = (X'WX)^X'WY kde X'WX: , oo oo OO oo T=0 T=0 V=0 T=0 oo , E(-^r«T \ T = 0 oo T = 0 E(-^)2m«T X'WY = f EX>w \ T = 0 OO T = 0 >=0 Tento přístup založený na vážené polynomiální metodě nejmenších čtverců se nazývá Brownův přístup. Značení: Pro dobrou srozumitelnost zavedeme následující značení. Nechť {Yt,t £ Z} je náhodná posloupnost, její realizace v časových okamžicích ti,t2, ■ ■ ■ ,tn označme yi, 7/2 > •••) Vn-Symbolem k včetně. Vt\k označme odhad hodnoty Yt v čase t na základě hodnot do časového okamžiku Jestliže k < t, pak yt|fc nazýváme predikcí. 2 Stochastické modely časových řad (9. cvičení) pokud k = t, ýt\t nazýváme filtrací je-li k = n > t, pak ýt\n nazýváme vyrovnáním (smoothing). 2 Jednoduché exponenciální vyrovnávání Exponenciální vyrovnávání pro m = 0 se nazývá jednoduché exponenciální vyrovnávání. i l-a' Použijeme-li označení /3q(í) = &o(í) a uvážíme—li, že pro a £ (0,1) je aT = t^—, dostaneme r=0 b0(t) J2 «r = E «ryí-- 6o(í) = £ = (1 - a) X>Tyt_T subst. k = r-l T=0 T=0 T=0 Abychom získali rekurentní vztah, upravujme Ýt = (l-a)E%0aTYt-T = (l-a)Yt + (l-a)E?=iaTYt = (l-a)Yt + (l-a)j:r=oak+1yt-i-k oo = (1 - a)Yt + a(l-a)Y, ^t-i-k = (1 - ot)Yt + ati k=0 ýt-! Protože predikce o r (r > 0) kroků dopředu pro jednoduché exponenciální vyrovnávání je rovna Yt+T\t = Yf = bo(t), můžeme předchozí rekurentní vztah přepsat pro realizace a dále upravovat Vt+i\t = (1 - ot)yt + až/í|í-i = (1 - a)yt + aytii-i + 2/ŕ|ŕ—i - ž/í|í-i = ýt\t-i + (1 - a) (yt - ž/í|í-i) chyba predikce ££|£_i a o rekurentním vzorci s chybou predikce Ét\t-i se říká, že je ve formě korekce chyby předpovědi (error correction form). 3 Ad hoc přístupy Holta a Winterse Pokud chceme na základě pozorování yi, ■ ■ ■ ,yt sestrojit předpověď budoucí hodnoty yt+i v čase t + 1, označme ji yt+i\t, Pak nejjednodušším odhadem může být obyčejmý průměr. Tato předpověď je vhodná, pokud hodnoty časové řady náhodně kolísají kolem střední hodnoty, která se v čase nemění. Jako rozumější se však jeví použít pro predikci budoucí hodnoty ve větší míře pozorování, která jsou časově nejbliže. Pak se nabízejí vážené průměry t-i yt+i\t = ^2wjityt-j, (1) j=0 kde součet vah je roven jedné, tj. Y^j=q wj,t = 1- Exponenciální vyrovnávání je založeno na myšlence použití vah, které do minulosti klesají exponenciálně. Stochastické modely časových řad (9. cvičení) 3 S využitím vztahu t± ■ 1 - a* £V =--, pro a G (0,1), (2) n 1 - a 3=0 chceme-li, aby součet vah, které exponenciálně klesají, byl roven jedné, položíme wht = T--joP . (3) 1 — cr Protože pro t —> oo konvergují váhy wjj —> Wj = (1 — ol)cx? , můžeme uvažovat jednokrokovou předpověď ze všech minulých pozorování ve tvaru oo Vt+i\t = (1 - oí) alyt-j, Pro Ol £ (0,1) • (4) j=o Analogicky jako u Brownova přístupu odvodíme rekurentní vztahy oo Vt+i\t = (1 - oí)yt + (l-a)J2 alyt-j 3=1 oo . k. 1-fc = (1 - a)yt + a(l - a) ^2 akyt k=o = (1 - a)yt + aytii-i Obdobně získáme i tvar využívající korekci chyby předpovědi vt+i\t = (i - ")ž/í + "ýí|í-i + ýí|í-i - vt\t-i = yt\t-i + (i - - yt|i-i) = Vt\t-i + (1 - a)éí|í-i Na tomto ad-hoc přístupu se nám podařilo ukázat, že se v podstatě jedná o jednoduché exponenciální vyrovnávání, které přepokládá model Yt =A)(r) +£í s lokální hladinou f3a(t). Použijeme-li značení obvyklá pro tento přístup, kdy váhy mají tvar Wj=l3(l-I3y, (5) tj. a = 1 — f3, místo f3a(t), píšeme Lt (level). Odvozené vztahy v novém značení: Vt+i\t = Pyt + (1 - P)Vt\t-i = Vt\t-i + Pet\t-i (6) Lt+1 = Lt + /3eí+i|i (7) 4 Stochastické modely časových řad (9. cvičení) 4 Holtovo exponenciální vyrovnávání Oproti jednoduchému exponenciálnímu vyrovnávání Holtova metoda předpokládá lokálně lineární trend, jehož koeficienty f3o(t) i /?i(í) se v čase mění. Hodnota časové řady v okamžiku t je určena jednak její úrovní /3o(t), jednak směrnicí /?i(í). V Holtově metodě se úroveň v čase t značí symbolem Lt (zkratka pro level) a směrnice jako T (zkratka pro trend). Úroveň Lt je zároveň vyrovnanou hodnotou realizace yt v okamžiku t. Směrnice lokálně lineárního trendu T (někdy mluvíme krátce o trendu) vyjadřuje očekávanou změnu úrovně časové řady při jednotkové časové změně. Pokud chceme pomocí Holtovy metody přepovídat hodnotu časové řady o h > 0 jednotek dopředu, položíme yt+h\t = Lt + Tth. (8) Takže, je-li h = 1, dostaneme jednokrokovou předpověď jako yt+1\t = Lt + Tt . (9) Protože by přibližně mělo platit, že realizace yt+i ~ Lt+i, pak se jeví vhodné získat Lt+i, jako konvexní lineární kombinaci hodnot (Lt+T) a yt+i- V Holtově metodologii bývá zvykem místo a G (0,1) používat f3 = 1 — a, takže konvexní lineární kombinace bude mít tvar Lt+1 = (l-P)(Lt + Tt)+pyt+1 . (10) Hodnota f3 se nazývá vyrovnávací konstanta pro úroveň řady. Analogickou úvahu použijeme i pro směrnici trendu Tt. Z přepokladu, že řada má lokálně lineární trend vyplývá, že by přibližně mělo platit Tt+i « Tt, ale zároveň má také smysl očekávat, že směrnice trendu je přibližně rozdílem sousedních úrovní, tj. Tt+i ~ Lt+i — Lt . Novou hodnotu směrnice Tt+i budeme uvažovat jako konvexní lineární kombinaci Tt+i = (l-7)Tt + 7(Lt+i-Lt), kde 7 G (0,1) (11) 7 je tzv. vyrovnávací konstanta pro lineární růst (pro směrnici). Na závěr odstavce ještě ukážeme přepsání předchozích rekurentních vztahů do chybového Stochastické modely časových řad (9. cvičení) 5 tvaru. Lt+1 = (1 - P)(Lt + Tt) + pyt+i = (1 - P)(Lt + Tt) + /3yí+i + /3yt+1,t - f3yt+1{t = /3(yt+1 - /3)(Lí + Tt) + /3yt+1|t £t+i|t Lt+Tt = /3et+l|t + Lt+Tt Tt+i = (1 - 7)Tt + 7(Lí+i - Lt) = Tt - 7Tt + 7 Lt+i -7Lt lt+tt+/3ít+1|t = Tt - 7Tt + 7(Lt + Tt + /3et+1|t) - 7^í = Tt + 7/3eí+i|i • 5 Holtovo-Wintersovo exponenciální vyrovnávání V případě, kdy časová řada má sezónní charakter, nevystačíme se žádnou z předchozích metod. Rozšíření Holtovy metody na sezónní časové řady je známo jako Holtova Wintersova metoda. Autorem je Holtův student Peter R. Winters. Holtova-Wintersova metoda je založena na třech vyrovnávacích konstantách. Jedna je pro hladinu, druhá pro trend a třetí pro sezónnost. Dle charakteru dat využívá aditivní nebo multiplikativní notaci. Uvažujme časovou řadu s lokálně lineárním trendem a sezónností s periodou p > 2. Stejně jako u Holtovy metody označme symbolem Lt úroveň v čase t, symbolem Tt směrnici lokálně lineárního trendu a symbolem St sezónní výkyv čase t. Součet úrovně Lt s hodnotou sezónního výkyvu St představuje v okamžiku t vyrovnanou hodnotu realizace yt. Předpověď hodnoty časové řady o h > 0 jednotek dopředu je pak dána vztahem yt+h\t = Lt + St-p+h + Tth, (12) takže v případě jednokrokové predikce platí Vt+i\t = Lt + St+i-p + Tt (13) Protože by mělo přibližně platit yt+i ~ Lt+i + St+i-p a Lt+i ~ Lt + Tt, má smysl získat úroveň Lt+i jako konvexní lineární kombinaci hodnot (Lt + St) a (yt+i - St+i-p), tj. Lt+1 = (1 - P){Lt + Tt) + p{yt+i - St+i-P). (14) 6 Stochastické modely časových řad (9. cvičení) Protože řada má lokálně lineární trend, mělo by přibližně platit Tt+i ~ Tt, ale zároveň lze směrnici lokálně lineárního trendu vyjádřit pomocí rozdílu sousedních hladin Tt+i ~ Lt+i — Lt. Oba předchozí vztahy využijeme při konstrukci směrnice lokálně linárního trendu díky konvexní linární kombinaci Tt+i = (l-7)Tt + 7(Lt+i-Lt), kde 7 £ (0,1) se nazývá vyrovnávací konstanta pro směrnici trendu. Pro sezónní výkyvy musí platit vztah St+i ~ St+i-p > a také St+i ~ yt+i — Lt+i Tedy označíme-li symbolem 5 £ (0,1) vyrovnávací konstantu pro sezónní výkyvy, pak Tt-\) = Tt-\ + aßet Optimální model je vybrán na základě tzv. AIC kritéria AIC = -2log(Likelihood) + 2m, kde m je počet parametrů. Při výstupu procedura ets() používá následující notaci Trendová Sezónní komponenta komponenta N A M (None) (Additive) (Multiplicative) N (None) N, N N, A N, M A (Additive) A, N A, A A, M Ad (Additive damped) Ad, N Ad, A Ad, M M (Multiplicative) M, N M, A M, M Md (Multiplicative damped) Md,N Md,A Md,M Označení optimálního modelu tvoří trojici ETS(E,T,S), kde E error možné hodnoty A, M T trend N, A, Ad, M, Md S seasonal N, A, M 8 Stochastické modely časových řad (9. cvičení) příklad 1 Metodu exponenciálního vyrovnávání aplikujeme na časovou řadu v souboru wool.txt. Jedná se o čtvrtletní údaje o množství vlny (v tunách) vyprodukované v Austrálii v letech 1965-1994. Datový soubor obsahuje pouze jeden sloupec s daty, proto bude nejjednodušší k načtení použít funkci scan: > íileDat <- paste(data.library, "wool.txt", sep = "") > wool <- scan(íileDat) Z načtených dat vytvoříme časovou řadu a vykreslíme ji. > woolTS <- ts(wool, start = 1965, frequency = 4) > str(woolTS) Time-Series [1:119] from 1965 to 1994: 6172 6709 6633 6660 6786 ... > par (mar = c(2, 2, 0, 0) + 0.5) > plot (woolTS) 1965 1970 1975 1980 1985 1990 1995 Obrázek 1: Produkce vlny v Austrálii v letech 1965-1994- Na načtená data nejprve použijeme jednoduché exponenciální vyrovnávání (lokálně konstantní trend). Potřebujeme určit jedinou vyrovnávací konstantu a. Model bude ve tvaru Yt+h = Lt, kde Lt je dáno vztahem Lt = aYt + (1 - a)Lt-i Stochastické modely časových řad (9. cvičení) 9 K tomu můžeme použít funkci HoltWinters, která je sice určená pro obecnější Holt— Wintersovo exponenciální vyrovnávání, kde počítáme nejen s úrovní (jako v případě jednoduchého exponenciálního vyrovnávání), ale i s lineárním trendem a sezónní složkou. Abychom odhadli pouze model jednoduchého exponenciálního vyrovnávání s jedinou vyrovnávací konstantou a, je třeba zadat, že nás nezajímají vyrovnávací konstanty f3 a 7, a toho dosáhneme zadáním argumentu beta=FALSE, gamma=FALSE. Odhadnutou konstantu Lf pro úroveň v čase t odpovídající poslednímu pozorování zobrazíme pomocí funkce coef f icients. V R je tato konstanta označena jako a. > x <- woolTS > modell <- HoltWinters (x, beta = FALŠE, gamma = FALŠE) > summary(modell) Lenj ;th Class Mode fitted 236 mts numeric X 119 ts numeric alpha 1 -none- numeric beta 1 -none- logical gamma 1 -none- logical coefficients 1 -none- numeric seasonal 1 -none- character SSE 1 -none- numeric call 4 -none- call > modell$alpha [1] 0.376407 > coefficients(modeliJ a 5711.117 Na základě tohoto modelu provedeme predikci na dva roky dopředu. Ta by měla vypadat takto Vt+h\t = Lt = a Výsledky jednoduchého exponenciálního vyrovnávání i predikci vykreslíme do grafu. > předl <- predict (modell, n. ahead = 12, prediction, interval = TRUE) > summary (předl) f it upr lwr Min. 5711 Min. 6951 Min. 3727 1st Qu. 5711 1st Qu. 7173 1st Qu. 3885 Median 5711 Median 7365 Median 4057 Mean 5711 Mean 7349 Mean 4073 3rd Qu. 5711 3rd Qu. 7537 3rd Qu. 4250 Max. 5711 Max. 7695 Max. 4471 10 Stochastické modely časových řad (9. cvičení) > par (mar = c(2, 2, 1, 0) + 0.5) > plot(modeli, predicted.values = předl) Holt-Winters filtering 1965 1970 1975 1980 1985 1990 1995 Obrázek 2: Jednoduché exponenciální vyrovnávání pomocí funkce HoltWinters () pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Z grafu je jasně patrné, že predikované hodnoty mají stejnou konstantní úroveň a ta se rovná odhadnuté úrovni pro poslední pozorování Lf - v R jí odpovídá konstanta a. Zkusíme použít poněkud složitější model - Holtovo exponenciální vyrovnávání, které nevyužívá lokálně konstantní trend jako v předchozím případě, ale lokálně lineární trend. Tento typ exponenciálního vyrovnávání vyžaduje určení dvou vyrovnávacích konstant a a f3. Model exponenciálního vyrovnávání má tvar Yt+h = Lt + Tth, kde L t a Tt jsou dány Lt = ayt + (l-a)(Lt_i+Tt_i) Tt = l3(Lt - Lt_i) + (1 - /3)Tt_i Opět použijeme funkci HoltWinters, tentokrát použijeme argument gamma=FALSE, protože sezónní složka nás zatím nezajímá. > model2 <- HoltWinters (x, gamma = FALŠE) > summary (model2) Length Class Mode fitted 351 mts numeric x 119 ts numeric alpha 1 -none- numeric Stochastické modely časových řad (9. cvičení) 11 beta 1 -none- numeric gamma 1 -none- logical coefficients 2 -none- numeric seasonal 1 -none- character SSE 1 -none- numeric call 3 -none- call > model2$alpha alpha 0.4180789 > model2$beta beta 0.1901893 > coefficients (model2) a b 5888.1247 205.5479 Na základě druhého modelu znovu provedeme predikci na dva roky dopředu. Ta bude ve tvaru Ýt+h\t = a + bí = 5888.1247 + 205.5479í. Výsledky jednoduchého exponenciálního vyrovnávání i predikci vykreslíme do grafu. > pred2 <- predict(model2, n.ahead = 12, prediction.interval = TRUE) > summary (pred2) f it upr lwr Min. 6094 Min. 7432 Min. 4016 1st Qu. 6659 1st Qu. 8503 1st Qu. 4406 Median 7224 Median 9767 Median 4682 Mean 7224 Mean 9880 Mean 4568 3rd Qu. 7789 3rd Qu. 11173 3rd Qu. 4784 Max. 8355 Max. 12694 Max. 4822 > par (mar = c(2, 2, 1, 0) + 0.5) > plot(model2, predicted.values = pred2) 12 Stochastické modely časových řad (9. cvičení) Holt-Winters filtering n-1-1-1-1-1-1— 1965 1970 1975 1980 1985 1990 1995 Obrázek 3: Holtova metoda pomocí funkce HoltWinters () pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Porovnáme-li tento graf s grafem pro jednoduché vyrovnávání, můžeme si všimnout rozdílu v predikci. U jednoduchého exponenciálního vyrovnávání je predikce konstantní, zatímco u Holtovy metody je to rostoucí lineární funkce. To souvisí s tím, že u jednoduchého exponenciálního vyrovnávání pracujeme s lokálním polynomem nultého stupně, u Holtovy metody s lokálním polynomem prvního stupně. Protože data, která máme k dispozici, jsou čtvrtletní, je možné (a při pohledu „od oka" to vypadá hodně pravděpodobně), že bychom v nich mohli objevit sezónní chování. Proto vyzkoušíme ještě třetí metodu - Holtovo-Wintersovo vyrovnávání, které bere v úvahu i sezón-nost. Model bude ve tvaru Ýt+h = Lt + Tth + St_p+h+, kde Lf, Tf a St jsou dány Lt = a(Yt - st-p) + (1 - a){Lt-! + Tt_i) Tt = f3(Lt - Lt-!) + (1 - p)Tt-i St = 7(yt - Lt) + (1 - 7)St_p Holtovo-Wintersovo vyrovnávání vyžaduje tři vyrovnávací konstanty a, f3 a 7. K jejich určení opět použijeme funkci HoltWinters, tentokrát již bez argumentů udávajících vynechané konstanty. > model3 <- HoltWinters (x) > suwmaxy (model3) Stochastické modely časových řad (9. cvičení) 13 Lenj ;th Class Mode fitted 460 mts numeric X 119 ts numeric alpha 1 -none- numeric beta 1 -none- numeric gamma 1 -none- numeric coefficients 6 -none- numeric seasonal 1 -none- character SSE 1 -none- numeric call 2 -none- call > model3$alpha alpha 0.6521521 > model3$beta beta 0.007690488 > model3$gamma gamma 0.6478487 > coefficients (model3) 6222.72882 b si s2 s3 s4 18.95975 -220.20379 -899.45453 -147.85037 168.78680 Hodnoty sezónních složek vykreslíme do grafu. > parGnar = c(2, 2, 0, 0) + 0.5) > plot (1:4, coefficients(model3) [3:6], type = "o") > for (k in 1:4) abline(v = k, col = "gray", lty = 1) 14 Stochastické modely časových řad (9. cvičení) Obrázek 4: Hodnoty sezónních složek z Holtova-Wintersova exponenciálního vyrovnávání pomocí funkce HoltWinters() pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Na základě tohoto modelu znovu provedeme predikci na dva roky dopředu. Výsledky Holtova-Wintersova exponenciálního vyrovnání i predikci vykreslíme do jednoho grafu. Predikci dostaneme ze vztahu Yt+h\t = a + b/i + St+p+i+(h-i) mod p, > pred3 <- predict(model3, n.ahead = 12, prediction.interval = T) > summary (pred3) f it upr lwr Min. 5361 Min. 6397 Min. 3419 1st Qu. 5894 1st Qu. 7256 1st Qu. 4152 Median 6152 Median 7706 Median 4410 Mean 6071 Mean 7713 Mean 4430 3rd Qu. 6329 3rd Qu. 8229 3rd Qu. 4764 Max. 6619 Max. 8883 Max. 5156 > parGnar = c(2, 2, 1, 0) + 0.5) > plot(model3, predicted.values = pred3) Holt-Winters filtering 1965 1970 1975 1980 1985 1990 1995 Obrázek 5: Holtovo-Wintersovo exponenciální vyrovnávání pomocí funkce HoltWinters() pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Dále se podíváme, jak dopadnou výsledky v případě použití funkce ets() z balíčku forecast. > library(forecastJ Stochastické modely časových řad (9. cvičení) 15 Funkce ets() umožňuje vyrovnávání jak v aditivních, tak v multiplikativních modelech (to umí i HoltWinters ()), a navíc umí odhadnout i exponenciální vyrovnávání s tlumícím faktorem. My zatím nebudeme uvažovat tak komplikované modely, provedeme obyčejné Holtovo-Wintersovo exponenciální vyrovnávání bez tlumícího faktoru a budeme chtít pouze aditivní model, to znamená, že uděláme přibližně totéž co v předchozím modelu (model3). Odhad vyrovnávacích konstant bez tlumícího faktoru provedeme v případě funkce ets() přidáním argumentu damped=FALSE, argumentem additive. only=TRUE vyloučíme možnost multiplikativních modelů. > x <- woolTS > modelETSl <- ets(x, damped = FALŠE, additive. only = TRUE) > summary (modelETSl) ETS(A,N,A) Call: ets(y = x, damped = FALŠE, additive.only = TRUE) Smoothing parameters: alpha = 0.7488 gamma = le-04 Initial states: 1 = 6644.9513 s=31.5458 408.2945 129.0941 -568.9344 sigma: 436.0035 AIC AICc BIC 2027.196 2027.946 2043.871 In-sample error measures: ME RMSE MAE -7.8706447 436.0035123 340.1038656 MASE 0.6055871 Prohlédneme si odhadnuté parametry. > modelETSl$par alpha gamma 1 sO sl 7.487666e-01 1.000492e-04 6.644951e+03 3.154583e+01 4.082945e+02 s2 1.290941e+02 Na základě AIC kritéria funkce ets() vybrala pro naše data jednoduché exponenciální vyrovnávání se sezónní složkou - tj. lokálně konstantní trend s odchylkami odpovídajícím čtvrtletím (značení ETS(A, N,A)). Všechny komponenty jsou aditivní. Odhadnutý model se liší od modelu 3 v tom, že se v něm nevyskytuje lineární trend. Výsledky exponenciálního vyhlazení opět vykreslíme. Pokud použijeme funkci plot na model, který je výstupem funkce ets, dostaneme graf, v němž jsou vykreslena zvlášť původní data, odhadnuté úrovně Lf, trendová komponenta Tt a cyklická komponenta St- Pokud zadáme argument plot .type=single, vše se zobrazí do jediného grafu. MPE MAPE -0.5567476 6.0128422 16 Stochastické modely časových řad (9. cvičení) > par (mar = c(2, 2, 1, 0) + 0.5) > plot(modelETSl) Decomposition by ETS(A,N,A) method 1965 1970 1975 1980 1985 1990 1995 Time Obrázek 6: Dekompozice pomocí funkce ets() pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 > par (mar = c(2, 2, 1, 0) + 0.5) > plot(modelETSl, plot.type = "single", col = 1:3, yláb = "") Decomposition by ETS(A,N,A) method o o o 00 o o o CD O o o o o o OJ o 1965 1970 1975 1980 1985 1990 1995 Obrázek 7: Holtovo-Wintersovo exponenciální vyrovnávání pomocí funkce ets() pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Stochastické modely časových řad (9. cvičení) Opět se podíváme na predikované hodnoty pro dva následující roky. > predETSl <- forecast(modelETSl) > summary (predETSl) Forecast method: ETS(A,N,A) Model Information: ETS(A,N,A) Call: ets(y = x, damped = FALSE, additive.only = TRUE) Smoothing parameters: alpha = 0.7488 gamma = le-04 Initial states: 1 = 6644.9513 s=31.5458 408.2945 129.0941 -568.9344 sigma: 436.0035 AIC AICc BIC 2027.196 2027.946 2043.871 In-sample error measures: ME RMSE MAE -7.8706447 436.0035123 340.1038656 MASE 0.6055871 Forecasts: Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 1994 Q4 5975. 187 5416. 426 6533. 948 5120. 636 6829. 738 1995 Ql 5374. 703 4676. 665 6072. 741 4307. 146 6442. 260 1995 Q2 6072. 705 5258. 858 6886. 552 4828. 034 7317. 377 1995 Q3 6351. 917 5436. 827 7267. 007 4952. 407 7751. 427 1995 Q4 5975. 187 4968. 989 6981. 384 4436. 341 7514. 033 1996 Ql 5374. 703 4284. 989 6464. 417 3708. 129 7041. 277 1996 Q2 6072. 705 4905. 415 7239. 995 4287. 489 7857. 921 1996 Q3 6351. 917 5111. 913 7591. 921 4455. 495 8248. 339 MPE MAPE -0.5567476 6.0128422 > par (mar = c(2, 2, 1, 0) + 0.5) > plot (predETSl) 18 Stochastické modely časových řad (9. cvičení) Forecasts from ETS(A,N,A) o o o oo o o o o o o co o o o LO o o o ^- 1965 1970 1975 1980 1985 1990 1995 Obrázek 8: Predikce pomocí funkce forecast () pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 Pro zajímavost si můžeme ještě na závěr vyzkoušet, jak by dopadlo exponenciální vyrovnání naší časové řady v případě, že bychom funkci ets() povolili použít všechny možnosti, které umí, včetně tlumícího faktoru a multiplikativních modelů. > x <- woolTS > modelETS2 <- ets(x) > summary (modelETS2) ETS(M,N,A) Call: ets(y = x) Smoothing parameters: alpha = 0.941 gamma = le-04 Initial states: 1 = 6650.654 s=33.2949 461.4852 135.7959 -630.576 sigma: 0.0748 AIC AICc BIC 2016.480 2017.230 2033.155 In-sample error measures: ME RMSE MAE MPE -6.3777596 443.0083548 349.5425906 -0.4032783 MASE MAPE 6.0751771 0.6223936 Stochastické modely časových řad (9. cvičení) 19 > par (mar = c(2, 2, 1, 0) + 0.5) > plot (modelETS2) Decomposition by ETS(M,N,A) method —i-1-1-i-1-1-r1 1965 1970 1975 1980 1985 1990 1995 Time Obrázek 9: Dekompozice pomocí funkce ets(): Produkce vlny v Austrálii v letech 1965-1994 > par(mar = c(2, 2, 1, 0) + 0.5) > plot(modelETS2, plot.type = "single", col = 1:3, yláb = "") Decomposition by ETS(M,N,A) method "~i-1-1-1-1-1-r 1965 1970 1975 1980 1985 1990 1995 "O CD > čB co o CD > CD O co CO CD co Obrázek 10: Holtovo-Wintersovo exponenciální vyrovnávání pomocí funkce ets() pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 20 Stochastické modely časových řad (9. cvičení) > predETS2 <- forecast(modelETS2) > summary (predETS2) Forecast method: ETS(M,N,A) Model Information: ETS(M,N,A) Call: ets(y = x) Smoothing parameters: alpha = 0.941 gamma = le-04 Initial states: 1 = 6650.654 s=33.2949 461.4852 135.7959 -630.576 sigma: 0.0748 AIC AICc BIC 2016.480 2017.230 2033.155 In-sample error measures: ME RMSE MAE MPE MAPE -6.3777596 443.0083548 349.5425906 -0.4032783 6.0751771 MASE 0.6223936 Forecasts: Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 1994 Q4 5969. 934 5397. 324 6542. 543 5094. 203 6845. 664 1995 Ql 5306. 109 4563. 842 6048. 377 4170. 909 6441. 310 1995 Q2 6072. 065 5142. 938 7001. 192 4651. 088 7493. 042 1995 Q3 6397. 812 5299. 973 7495. 651 4718. 813 8076. 812 1995 Q4 5969. 934 4746. 830 7193. 037 4099. 358 7840. 509 1996 Ql 5306. 109 3992. 769 6619. 450 3297. 529 7314. 690 1996 Q2 6072. 065 4642. 732 7501. 397 3886. 090 8258. 040 1996 Q3 6397. 812 4851. 590 7944. 034 4033. 069 8762. 555 > par (mar = c(2, 2, 1, 0) + 0.5) > plot (predETS2) Stochastické modely časových řad (9. cvičení) 21 o o o oo o o o o o o to o o o LO o o o ^- 1965 Forecasts from ETS(M,N,A) 1970 1975 1980 1985 1990 1995 Obrázek 11: Predikce pomocí funkce forecastO pro časovou řadu: Produkce vlny v Austrálii v letech 1965-1994 7 Ukol: Aplikujte exponenciální vyhlazování na časovou řadu s měsíčními údaji o počtu smrtelných úrazů v USA v letech 1973-1978 v souboru deaths.dat.