Stochastické modely časových řad (6. cvičení) 1 M5201 - 6. CVIČENÍ: Generování náhodných procesů (búý šum, IID proces; AR, M A, ARMA procesy), odhad jejich ACF a PACF 1 Bílý šum a IID proces Náhodný proces {st,t G T} je bílý šum (značíme Et ~ WN(0,a2), jestliže Et jsou nekorelované náhodné veličiny s nulovou střední hodnotou a konstantním rozptylem, tj. platí Est = 0, Det = a2, C{et,e8) = 0, s^t. Bílý šum, v němž jsou náhodné veličiny Et nejen nekorelované, ale také nezávislé, se nazývá IID proces a značí se Et ~ IID(0, a2). Autokorelační funkce stacionárního náhodného procesu {Xt, t G T} je definována vztahem g{) = wr kde 7(í) je autokovarianční funkce a 7(0) = D{Xt). Parciální autokorelační funkce stacionárního náhodného procesu {Xt, t G T} G L (u, A, P) je definována vztahem a(l)=R(Xt,Xt+1), a(k) = R(Xt - Xt, Xt-k - Xt-k) pro \k\ > 1, kde Xt, resp. Xt-k jsou nejlepší lineární predikce Xt, resp. Xt-k pomoci Xt-k+i, ■ ■ ■ ,Xt-±. PŘÍKLAD 1 Mějme náhodný proces {et,t G T}, kde všechny náhodné veličiny Et mají normální rozdělení N(0,a2) a jsou navzájem nezávislé. Je zřejmé, že takový proces představuje bílý šum (jedná se dokonce o IID proces). Hodnotu rotpylu zvolme a2 = 1 a vygenerujeme 300 hodnot tohoto procesu. Vygenerovanou řadu vykreslíme. > x <- morm(300, mean = 0, sd = 1) > plot (x, type = "1", xláb = "t") 2 Stochastické modely časových řad (6. cvičení) 50 —I— 200 —I— 250 —r~ 300 100 150 Obrázek 1: Simulovaná data bílého šumu st ~ N(0,1) Nyní prozkoumáme autokorelační a parciální autokorelační funkci bílého šumu. Vyjdeme při tom ze znalosti autokovarianční funkce pro bílý šum 7(fc) Odtud ihned vidíme, že g(k) _ íl k = 0 [O jinak a(k)=0 k = 1,2, Podívejme se, jak tomuto faktu odpovídají odhady autokorelační funkce a parciální autokorelační funkce, které získáme na základě vygenerovaných dat. K výpočtu a vykreslení výběrové autokorelační funkce použijeme v R funkci acf () a pro výběrovou parciální autokorelační funkci existuje analogická funkce pacf (). > acfWN <- acf(x, main = "Výberová ACF", xlab = "k") Výberová ACF n _L_L i i i Obrázek 2: Výběrová autokorelační funkce bílého šumu st ~ N(0,1) V obrázku odhadnuté autokorelační funkce jsou dvě modré přerušované čáry. Ty představují 95% pás spolehlivosti pro autokorelační funkci bílého šumu. To znamená, že v případě Stochastické modely časových řad (6. cvičení) 3 bílého šumu by hodnoty výběrové autokorelační funkce pro k > 0 měly ležet s pravděpodobností 95 % uvnitř tohoto pásu. > pelcíWN <- pacf(x, main = "Výběrová PACF", xláb = "k", ylim = c(-0.2, 1)) Výberová PACF I I I I Obrázek 3: Výběrová parciální autokorelační funkce bílého šumu st ~ N(0,1) Na obrázku opět vidíme, že hodnoty výběrové parciální autokorelační funkce leží uvnitř 95% pásu spolehlivosti zkonstruovaného pro případ bílého šumu, což souhlasí s tím, že naše zkoumaná řada byla simulována jako realizace bílého šumu. Pro testování nulové hypotézy Hq: časová řada představuje bílý šum byl odvozen Box-Piercův a Box-Ljungův test. Vyzkoušejme oba tyto testy na simulovaná data bílého šumu. > Box.test00 Box-Pierce test data: x X-squared = 2.2612, df = 1, p-value = 0.1327 > Box. test (x, type = "Ljung-Box") Box-Ljinig test data: x X-squared = 2.2839, df = 1, p-value = 0.1307 V obou případech jsme dostali p-hodnoty (p-value) větší než 0.05, a proto testy nezamítají hypotézu, že testovaná data jsou bílým šumem (což, jak víme, skutečně platí). 4 Stochastické modely časových řad (6. cvičení) 2 Kauzální AR procesy Uvažujme autoregresní model p-tého řádu. Tento model je definován vztahem: Yt - (fiYt-i-----(z) = 1 — ip±z — (f2Z2 — ... — (fpzp ležely vně jednotkové kružnice. Pro autokorelační a parciální autokorelační funkci kauzálního AR(p) procesu platí následující: • g{k) klesá pro k —> oo exponenciálně k nule • Je-li {Yt,t £ Z} centrovaný nedegenerovaný kauzální AR(p) proces, pak platí: (1) a(p) = (fp (2) a(p) = 0 pro k > p (důležitá identifikační vlastnost AR(p) procesu) PŘÍKLAD 2 Dále budeme pracovat s AR procesem druhého řádu s konkrétními hodnotami tpi = 0.9 a (f2 = —0.25, tj. proces zadaný rovnicí Yt = 0.9Yt-i -0.25Yt_2+Et. U takto zadaného modelu nevíme, zda se jedná o kauzální nebo nekauzální proces. Abychom to byli schopni zjistit, potřebujeme znát kořeny polynomu 1 — if\z — ■ ■ ■ — (ppzp . K tomu použijeme funkci polyroot (), která předpokládá, že koeficienty polynomu řádu p jsou zadány v pořadí od nejnižší mocniny k nej vyšší, tj. aQ,a±,... ,ap . Až budeme znát kořeny polynomu Q (z), pomocí funkce Mod() zjistíme jejich absolutní hodnoty. > ARphi <- c(0.9, -0.25) > ARroots <- polyroot(c(1, -ARphi)) > Mod(ARroots) [1] 2 2 Pro větší názornost si ještě vykresleme obrázek, v němž bude znázorněna poloha kořenů polynomu $>(z) vůči jednotkové kružnici. Stochastické modely časových řad (6. cvičení) 5 > x <- seq(-l, 1, length = 1000) > y <- sqrt(l - x'2) > plot(Re(ARroots), Im(ARroots), type = "p", pch = 19, xlab = "Reálna cast", ylab = "Imaginární cast", main = "Jednotková kružnice", ylim = c(-2.5, 2.5), xlim = c(-2.5, 2.5)) > lines(c(x, rev(x)), c(y, -rev(y))) > ábline(h = 0) > abline(v = 0) > legend(c(-2.2), legend = "Kořeny procesu AR(2)", pch = 19, bty = "n") Jednotková kružnice CM - CjJ _ • Kořeny procesu AR(2) -2-10 1 Reálna cast Obrázek 4: Poloha kořenů polynomu §{z) vzhledem k jednotkové kružnici Vidíme, že kořeny leží vně jednotkové kružnice, takže jde o kauzální AR proces. Nyní pomocí R vygenerujeme 300 hodnot takto zadaného AR procesu. Použijeme k tomu funkci arima. sim(), která je určena ke generování ARIMA procesů. Parametry těchto procesů se zadávají do argumentu model a musí se jednat o datový objekt typu seznam (list), který má položky s názvy ar, ma a order. V případě AR procesu stačí zadat seznam s jedinou položkou ar obsahující vektor s koeficienty tpi,... ,(pp. Po vygenerování získanou časovou řadu vykreslíme. > ar.sim <- arima. sim(model = list(ar = ARphi), n = 300) > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot(ar.sim) 6 Stochastické modely časových řad (6. cvičení) ■* -ľ OJ I i t "1-1- 0 50 100 150 Obrázek 5: Simulovaná data AR procesu Yt —I- 200 = 0.9Yt- —i-r 250 300 - 0.25yt-2 + et Pomocí funkce ARMAacf () můžeme vypočítat teoretickou autokorelační funkci p{k) (ACF) a teoretickou parciální autokorelační funkci a{k) (PACF) AR procesu Yt = 0.9Yí_i — 0.25Yí_2 + Et- Výběrové protějšky těchto funkcí vypočtené na základě vygenerovaných hodnot procesu získáme opět pomocí funkcí acf (), popř. pacf (). > tARacf <- ARMAacf(ar = ARphi, lag.max = 20) > par(mírow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > plot(0: (length(tARacf) - í), tARacf, type = "h", main = "Teoretická ACF") > áblineCh = 0) Teoretická ACF ~l-1-1-1- 0 5 10 15 Obrázek 6: Teoretická ACF pro AR proces Yt = 0.9Yí-i —r- 20 0.25yt-2 + Et > vARací <- ací(ar.sim, lag.max = 20, plot = FALSE) > par (mírow = c(l, 1), mar = c (2, 2, 3, 0) + 0.05) > plot(vARacf, main = "Výberová ACF") Stochastické modely časových řad (6. cvičení) 7 Výberová ACF J_L T" t J_L ~i-1-1-1-r 0 5 10 15 20 Obrázek 7: Výběrová ACF pro simulovaná data AR procesu Yt = 0.9Yí_i — 0.25Yi_2 + £t Z grafů je dobře vidět že jak teoretická, tak výběrová autokorelační funkce klesá exponenciálně k nule. > tARpací <- ARMAací (ar = ARphi, lag.max = 20, pací = TRUE) > par(mírow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > plot(l:length(tARpacf), tARpací, type = "h", main = "Teoretická PACF") > áblineCh = 0) Teoretická PACF Obrázek -1-1-1- 5 10 15 Teoretická PACF pro AR proces Yt = 0.9Yt-! 20 0.25yt-2 + et > vARpací <- pací (ar. sim, lag.max = 20, plot = FALSE) > par (mírow = c(l, 1), mar = c (2, 2, 3, 0) + 0.05) > plot (vARpací, main = "Výberová PACF") 8 Stochastické modely časových řad (6. cvičení) Výberová PACF CD O O C\J O CM CĎ 5 10 15 20 Obrázek 9: Výběrová PACF pro simulovaná data AR procesu Yt = Q.QYt-i — 0.25YÍ-2 + £t Teoretická parciální autokorelační funkce nabývá nenulových hodnot pouze pro k = 1 a k = 2 a v souladu s tím je také její výběrový protějšek, kde pro k > 2 jsou hodnoty výběrové parciální autokorelační funkce blízké nule a leží uvnitř 95%ního pásu spolehlivosti pro bílý šum. Pomocí funkce arima() můžeme provést odhad parametrů AR(2) modelu. Při použití této funkce je nutné do argumentu order zadat řád odhadovaného ARIMA modelu. AR proces v našem příkladu je speciálním případem ARIMA procesu s řádem c(2,0,0). > ar.fit <- arimaCar. sim, order = c(2, 0, 0)) > print(ar.fit) Series: ar.sim ARIMA(2,0,0) with non-zero mean Call: arima(x = ar.sim, order = c(2, 0, 0)) Coefficients: arl ar2 intercept 0.8714 -0.2373 -0.1692 s.e. 0.0560 0.0561 0.1466 sigma~2 estimated as 0.8694: log likelihood = -405.08 AIC = 818.17 AICc = 818.3 BIC = 832.98 Jak je vidět z výsledků, odhady se příliš neliší od skutečných hodnot parametrů. Kvalitu vyrovnání dat pomocí právě odhadnutého modelu můžeme posoudit na základě reziduí. Ve statistické knihovně k tomu máme k dispozici funkci tsdiagO, která se aplikuje na odhadnutý AR model. > tsdiagCar.fit) Stochastické modely časových řad (6. cvičení) 9 Standardized Residuals U A L 1.1,1,1, JÍnlil 1111, hli lll .11 1.1 h .11 i 1 ll lllllj, ill J H 1,1.1 J - |'V'TI|I } i1 n Til' l i■ r h j f T 1 n III í'111'ľ jllp TI f ľ ľ II 'It 'l'íp ~i-1-1-1-1-1-r- O 50 100 150 200 250 300 Time p values for Ljung-Box statistic ° -1-1-1-1-1- 2 4 6 8 10 lag Obrázek 10: Diagnostické grafy pro odhadnutý AR(2) model pro simulovaná data AR procesu Yt = 0.9Yt-i ~ 0.25Yí_2 + et Jestliže byl model zvolen a odhadnut správně, rezidua by měla představovat bílý šum. První graf zobrazuje samotná rezidua. Ta by neměla vykazovat žádnou systematičnost. V druhém grafu je vykreslena výběrová autokorelační funkce reziduí, která by měla mít významně nenulovou hodnotu pouze v bodě k = 0. Třetí graf ukazuje p-hodnoty Box-Ljungova testu, zda rezidua (pro různá zpoždění k) představují bílý šum. Je-li model určen správně, měly by všechny hodnoty ležet nad hladinou významnosti testu, která je v grafu znázorněna přerušovanou čarou. Všechny grafy nasvědčují tomu, že AR proces 2. řádu představuje vhodný model, protože po jeho odstranění z časové řady jsme získali bílý šum. Vzhledem k tomu, že jsme data vygenerovali právě pomocí tohoto AR procesu, je výsledek v souladu s očekáváním. Nyní, když jsme pro data nalezli vhodný model, jsme schopni predikovat hodnoty procesu pro několik kroků dopředu. Pro tento účel je v souboru FunkceM5201 .R vytvořena funkce PlotPredictARIMA. Abychom ji mohli používat, musíme tento soubor načíst příkazem source(). Pokud by nás zajímal zdrojový kód funkce, stačí jako příkaz zadat do R název funkce PlotPredictARIMA a zdrojový kód se vypíše na obrazovku. Funkci použijeme pro predikci o 10 kroků dopředu. Jako argumenty zadáme zkoumanou časovou řadu (ar.sim), odhadnutý model pro tuto řadu (ar.fit) a počet predikovaných hodnot (n.ahead=10). 10 Stochastické modely časových řad (6. cvičení) > fileSkript <- paste(data.library, "FunkceM5201.R", sep = "") > source(fileSkript) > par (mlrow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > PlotPredictARIMA(ar.sim, ar.fit, n.ahead = 10) 4 -r 0 50 100 150 200 250 300 Obrázek 11: Predikce na základě odhadnutého AR(2) modelu pro simulovaná data AR procesu Yt = 0.9yt_i - 0.25yt-2 + et Když jsme odhadovali parametry AR procesu pomocí funkce arima, museli jsme dopředu určit řád odhadovaného modelu. V knihovně f orecast máme k dispozici funkci auto. arima O, kterou také můžeme použít k odhadu modelu pro naše simulovaná data. Tato funkce se od funkce arima liší předpokladem, že vstupní data jsou generována obecným ARIMA procesem a tím, že v průběhu výpočtu zároveň zjišťuje vhodný řád modelu a zároveň odhaduje parametry pro tento model. Zadáme-li volbu trace=TRUE, na obrazovku se nám vypíší všechny typy ARIMA modelů, které funkce aut o. arima O brala při výpočtech v úvahu. Vyzkoušejme tedy, jaký model odhadneme pomocí této funkce. > library(forecast) > ar.íitF <- auto.arimaCar.sim, trace = TRUE) ARIMA(2,0,2) with non-zero mean 820 6625 ARIMA(0,0,0) with non-zero mean 1038.305 ARIMA(1,0,0) with non-zero mean 834 5473 ARIMA(0,0,1) with non-zero mean 872 167 ARIMA(1,0,2) with non-zero mean 821 9867 ARIMA(3,0,2) with non-zero mean 822 1195 ARIMA(2,0,1) with non-zero mean 818 8195 ARIMA(2,0,1) with zero mean 818 033 ARIMA(1,0,1) with zero mean 820 6264 ARIMA(3,0,1) with zero mean 819 195 ARIMA(2,0,0) with zero mean 817 0522 ARIMA(2,0,0) with non-zero mean 817 9045 ARIMA(1,0,0) with zero mean 833 37 ARIMA(3,0,0) with zero mean 817 1965 Best model: ARIMA(2,0,0) with zero mean > priu t Car. f i tF) Stochastické modely časových řad (6. cvičení) 11 Series: ar.sim ARIMA(2,0,0) with zero mean Call: auto.arima(x = ar.sim, trace = TRUE) Coefficients: arl ar2 0.8750 -0.2343 s.e. 0.0561 0.0561 sigma~2 estimated as 0.8731: log likelihood = -405.74 AIC = 817.47 AICc = 817.56 BIC = 828.59 V knihovně forecast je ještě k dispozici funkce forecast (), kterou použijeme při predikci pro právě odhadnutý model. > par(mfrow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > plot(forecast(ar.fitF), n = 10) Forecasts from ARIMA(2,0,0) with zero mean 0 50 100 150 200 250 300 Obrázek 12: Predikce na základě odhadnutého ARIMA(2,0,0) modelu pro simulovaná data AR procesu Yt = 0.9Yt_i - 0.25^-2 + et 3 Invertibilní MA procesy Uvažujme nyní MA model g-tého řádu: Yt = Et + OiEt-i H-----h 0q£t-q ■ MA proces Yt ~ MA(q) se nazývá invertibilní, jestliže je možné ho zapsat jako proces AR(oo), tj. Et = Y^jLo^jYt-j- MA(q) proces je vždy slabě stacionární. Nutnou a postačující podmínkou invertibility MA procesu je, aby všechny kořeny polynomu Q (z) = 1 + 0\z + 02Z2 + • • • + 9qzq ležely vně jednotkové kružnice. Pro autokorelační a parciální autokorelační funkci invertibilního MA(q) procesu platí následující: • r?(fc) Je nulová pro všechna k > q (důležitá identifikační vlastnost MA(q) procesů) 12 Stochastické modely časových řad (6. cvičení) • Neexistuje takové k$ £ N, že by pro k > k$ platilo a{k) = 0. PŘÍKLAD 3 V příkladu budeme pracovat s MA procesem druhého řádu s konkrétními hodnotami 0\ = 0.75 a 62 = 0.2, tj. s procesem tvaru Yt = et + 0.75£í_i + 0.2£í_2 U takto přímo zadaného modelu nevíme, zda se jedná o invertibilní či neinvertibilní proces. Abychom byli schopni rozhodnout o invertibilitě zadaného procesu, budou nás zajímat kořeny polynomu 1 + 0\z + • • • + 9qzq . K tomu opět použijeme funkci polyrootO, která předpokládá koeficienty polynomu v pořadí aQ,a±,... ,aq . Jestliže kořeny leží vně jednotkové kružnice, MA proces je invertibilní. > MAtheta <- c(0.75, 0.2) > print(MAroots <- polyroot(c(l, MAtheta))) [1] -1.875+1.218349i -1.875-1.218349i > Mod(MAroots) [1] 2.236068 2.236068 Kořeny opět znázorníme jako body v komplexní rovině, abychom viděli jejich polohu vůči jednotkové kružnici. > x <- seq(-l, 1, length = 1000) > y <- sqrttt - x'2) > plot (Re (MAroots) , Im(MAroots) , type = "p", pch = 19, xlab = "Reálna cast", ylab = "Imaginárni cast", main = "Jednotková kružnice", ylim = c (-3, 3), xlim = c (-3, 3)) > lines(c(x, rev(x)), c(y, -rev(y))) > ábline(h = 0) > abline(v = 0) > legend("bottomright", legend = "Kořeny procesu MA(2)", pch = 19) Stochastické modely časových řad (6. cvičení) 13 Jednotková kružnice • • • Kořeny procesu MA(2) -3-2-10 1 2 3 Reálna cast Obrázek 13: Poloha kořenů polynomu Q(z) vzhledem k jednotkové kružnici Vidíme, že kořeny leží vně jednotkové kružnice, takže jde o invertibilní MA proces. Pomocí funkce arima. sim() provedeme simulaci pro 300 hodnot. Jak víme, funkce arima. sim() simuluje data obecně pro zadaný ARIMA model. V našem případě, kdy generujeme časovou řadu MA procesu 2. řádu (což je speciální připad ARIMA procesu), zadáme parametry modelu do argumentu model jako seznam s jedinou složkou ma = MAtheta. Vygenerovaná data vykreslíme. > jna.sim <- arima. sim(model = list(ma = MAtheta.), n = 300) > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot (ma. sim) 14 Stochastické modely časových řad (6. cvičení) co — co _ I O 50 100 150 200 250 300 Obrázek 14: Simulovaná data MA procesu Yt = £t + 0.75et_i + 0.2et-2 Pomocí funkce ARMAacf () můžeme vypočítat teoretickou autokorelační funkci p{k) (ACF) a teoretickou parciální autokorelační funkci a{k) (PACF) MA procesu Yt = £t + 0.75et_i + 0.2et-2- Výběrové protějšky těchto funkcí získáme pomocí funkcí acf (), resp. pacf (). > tMAací <- ARMAacf (ma = MAtheta, lag.max = 20) > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot(0: (length(tMAacf) - 1), tMAací, type = "h", main = "Teoretická ACF") > ábline(h = 0) Teoretická ACF ~i-1-1-1-r 0 5 10 15 20 Obrázek 15: Teoretická ACF pro MA proces Yt = et + 0.75^-1 + 0.2et-2 > vMAací <- ací(ma.sim, lag.max = 20, plot = FALŠE) > par(mírow = c(l, 1), mar = c(2, 2, 3, 0) + 0.05) > plot(vMAacf, main = "Výberová ACF") Stochastické modely časových řad (6. cvičení) 15 Výberová ACF n-1- t ~i-1-1-1- 0 5 10 15 Obrázek 16: Výběrová ACF pro simulovaná data MA procesu Yt J_L 20 et + 0.75£í_i + 0.2£í_2 Vidíme, že teoretická autokorelační funkce je nenulová pouze pro k = 0,1, 2. Tomu odpovídá i výběrová autokorelační funkce, která je významně nenulová pouze v těchto bodech a všude jinde její hodnoty leží uvnitř 95%ního pásu spolehlivosti pro bílý šum. > tMApací <- ARMAacf (ma = MAtheta, lag.max = 20, pací = TRUE) > par(mírow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > plot(l:length(tMApacf), tMApací, type = "h", main = "Teoretická PACF") > áblineCh = 0) Teoretická PACF -1-1-1-r 5 10 15 20 Obrázek 17: Teoretická PACF pro MA proces Yt = et + 0.75^-1 + 0.2^-2 > vMApací <- pacíCma.sim, lag.max = 20, plot = FALSE) > par (mírow = c(l, 1), mar = c (2, 2, 3, 0) + 0.05) > plot (vMApací, main = "Výberová PACF") 16 Stochastické modely časových řad (6. cvičení) Výberová PACF "1-r -1-1-1-r 5 10 15 20 Obrázek 18: Výběrová PACF pro simulovaná data MA procesu Yt = £t + 0.75et-i + 0.2et-2 Vidíme, že teoretická parciální autokorelační funkce klesá v absolutní hodnotě exponenciálně k nule a tomu odpovídá i průběh výběrové parciální autokorelační funkce. Pomocí funkce arima() provedeme odhad parametrů MA(2) modelu na základě simulovaných dat. > ma.fit <- arima(ma. sim, order = c(0, 0, 2)) > print (ma. f i t) Series: ma.sim ARIMA(0,0,2) with non-zero mean Call: arima(x = ma.sim, order = c(0, 0, 2)) Coefficients: mal ma2 intercept 0.7108 0.1472 -0.0703 s.e. 0.0561 0.0575 0.1019 sigma~2 estimated as 0.9057: log likelihood = -411.09 AIC = 830.19 AICc = 830.32 BIC = 845 I v tomto případě jsou odhady blízké skutečným hodnotám parametrů. Vhodnost použitého modelu posoudíme podle reziduí. Diagnostické grafy získáme opět díky funkci tsdiagO, kterou použijeme na odhadnutý MA model. > tsdiagCma.fit) Stochastické modely časových řad (6. cvičení) 17 Standardized Residuals MA ľ |ľ|ľ| f' Tlľ ľiľllľ iJi,,,,, iJii......U J j. i ii. ,.m i.ii Jiud jLí .i i .i .n iJjI mi t™ ľ......II llll'ľ ľ] i |"i| i|p ij'i ni'iijf'iiii' rnr 150 Time ACF of Residuals -i-'-I-r- p values for Ljung-Box statistic Obrázek 19: Diagnostické grafy pro odhadnutý MA(2) model pro simulovaná data MA procesu Yt = st + 0.75£í_i + 0.2£í_2 Žádný z diagnostických grafů nenaznačuje, že by náš model byl nevhodný. Rezidua nevykazují žádnou systematičnost, výběrová autokorelační funkce kromě k = 0 leží v 95%ním pásu spolehlivosti pro bílý šum a ani Ljung-Boxův test nezamítá hypotézu o tom, že rezidua představují bílý šum. Odhadnutý model použijeme k predikci o 10 kroků dopředu pomocí funkce PlotPredictARIMA. > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > PlotPredictARIMA (ma., sim, ma. fit, n. ahead = 10) 18 Stochastické modely časových řad (6. cvičení) 3 - -3 - O 50 100 150 200 250 300 Obrázek 20: Predikce na základě odhadnutého M A (2) modelu pro simulovaná data MA procesu Yt = Et + 0.75£í_i + 0.2et-2 Ještě provedeme obecnější odhad modelu, kdy použijeme funkci auto.arima(). Již víme, že tato funkce odhaduje zároveň řád ARIMA modelu i jeho parametry > ma.fitF <- auto. arima(ma. sim, trace = TRUE) ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA (2,0,2 (0,0,0 (1,0,0 (0,0 (1,0 (3,0 (2,0 (2,0 (1,0 (3,0 (2,0,0 (2,0,2 (1,0,0 (3,0,2 with non-zero mean 829 2493 with non-zero mean 955 3027 with non-zero mean 853 376 with non-zero mean 834 5677 with non-zero mean 834 5871 with non-zero mean 831 7525 with non-zero mean 827 8052 with zero mean 826 4774 with zero mean 831 0727 with zero mean 828 398 with zero mean 829 5102 with zero mean 827 8355 with zero mean 851 7538 with zero mean 830 3901 Best model: ARIMA(2,0,1) with zero mean > priu t (ma. f i tF) Series: ma.sim ARIMA(2,0,1) with zero mean Call: auto.arima(x = ma.sim, trace = TRUE) Coefficients: arl ar2 mal 0.3151 -0.0629 0.4022 s.e. 0.2033 0.1343 0.1993 sigma~2 estimated as 0.9063: log likelihood = -411.2 AIC = 830.39 AICc = 830.53 BIC = 845.21 Stochastické modely časových řad (6. cvičení) 19 Druhý odhadnutý modelu použijeme při volání funkce forecastO, pomocí níž budeme chtít zjistit predikci o 10 kroků dopředu. > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot(forecast(ma.fitF), n = 10) Forecasts from ARIMA(2,0,1) with zero mean 0 50 100 ŕ 150 200 250 300 Obrázek 21: Predikce na základě odhadnutého APJMA(2,0,0) modelu pro simulovaná data MA procesu Yt = s + 0.75£t_i + 0.25£t_2 4 Kauzální a invertibilní ARMA procesy Uvažujme proces ARMA(p,q). Tento model je definován vztahem: Yt - WiYt-i-----^fpYt-p = st + Oist-i H-----h OgSt-q ■ Kauzalita a invertibilita jsou definovány obdobně jako v případě AR a MA procesu: • ARMA proces Yt ~ ARMA(p, q) se nazývá kauzální, jestliže je možné ho zapsat jako proces MA(oo), tj. Yt = YfjLoýjtt-j- • ARMA proces Yt ~ ARMA(p, q) se nazývá invertibilní, jestliže je možné ho zapsat jako proces AR(oo), tj. st = Ylj^o^jYt-j- Analogicky platí, že nutnou a postačující podmínkou kauzality ARMA(p,q) procesu je, aby všechny kořeny polynomu $>(z) = 1 — k$ platilo g{k) = 0. Od k = q — p klesá autokorelační funkce exponenciálně k nule. • Neexistuje takové ko, že by pro všechna k > ko platilo a{k) = 0. Od k = p — q klesá parciální autokorelační funkce exponenciálně k nule. 20 Stochastické modely časových řad (6. cvičení) PŘÍKLAD 4 Dále budeme pracovat s ARMA(2,2) procesem s konkrétními hodnotami tpi = 0.9, ARphi <- c(0.9, -0.5) > ARroots <- polyroot(c(l, -ARphi)) > Mod(ARroots) [1] 1.414214 1.414214 Vidíme, že kořeny leží vně jednotkové kružnice, takže jde o kauzální ARMA proces. Teď najdeme kořeny polynomu 1 + 0.75z + 0.25z2 a zjistíme jejich absolutní hodnoty. > MAtheta <- c(0.75, 0.25) > MAroots <- polyroot(c(1, MAtheta)) > Mod(MAroots) [1] 2 2 Vidíme, že kořeny leží vně jednotkové kružnice, takže jde o invertibilní ARMA proces. Ještě vykresleme polohu kořenů pro oba polynomy vzhledem k jednotkové kružnici v komplexní rovině. > par (m írow = c(l, 2)) > plot (Re (ARroots) , Im(ARroots) , type = "p", pch = 19, xlab = "Reálna cast", ylab = "Imaginárni cast", main = "Kořeny AR časti", ylim = c (-2, 2), xlim = c (-2, 2)) > lines(c(x, rev(x)), c(y, -rev(y))) > ábline(h = 0) > abline(v = 0) > plot (Re (MAroots) , Im(MAroots) , type = "p", pch = 19, xlab = "Reálna cast", ylab = "Imaginárni cast", main = "Kořeny MA časti", ylim = c (-2, 2), xlim = c (-2, 2)) > lines(c(x, rev(x)), c(y, -rev(y))) > ábline(h = 0) > abline(v = 0) Stochastické modely časových řad (6. cvičení) 21 Kořeny AR časti Kořeny MA časti T" -r -1 o 1 Reálna cast -r -1 o 1 Reálna cast Obrázek 22: Poloha kořenů polynomů $>(z) a Q(z) vzhledem k jednotkové kružnici Pomocí funkce arima.simO provedeme simulaci 300 hodnot procesu. Jako argument model zadáme seznam se dvěma položkami ar (parametry popisující AR část ARMA procesu) a ma (parametry popisující MA část ARMA procesu). > par (mírow = c(l, 1)) > arma.sim <- arima. s im (model = list (ar = ARphi, ma = MAtheta) , n = 300) > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot (arma.sim) 0 50 100 150 200 250 300 Obrázek 23: Simulovaná data ARMA procesu Yt - 0.9Yí_i + 0.5Yt_2 = et + 0.75^-1 + 0.25£í_2 Pomocí funkce ARMAacf () vypočítáme teoretickou autokorelační funkci p(k) (ACF) a teoretickou parciální autokorelační funkci a(k) (PACF). Výběrové protějšky těchto funkcí získáme opět pomocí funkcí acf (), popř. pacf (). > tARMAací <- ARMAacf(ar = ARphi, ma = MAtheta, lag.max = 20) > par(mfrow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > plot(0: (length(tARMAací) - í), tARMAací, type = "h", main = "Teoretická ACF") > ábline(h = 0) 22 Stochastické modely časových řad (6. cvičení) Teoretická ACF 0 5 10 Obrázek 24: Teoretická ACF pro ARMA proces Yt- 15 20 -0.9yt_i+0.5yt_2 = et+0.75et_i+0.25et_2 > vARMAací <- ací(arma.sim, lag.max = 20, plot = FALSE) > par (mírow = c(l, 1), mar = c (2, 2, 3, 0) + 0.05) > plot (vARMAací, main = "Výberová ACF") Výberová ACF .0 0.2 0.4 0.6 0.8 1.0 ...... ... 1.......................................................................... \ ...... 1 1 -0.2 0. ... ..................................... ....... .................. —r 0 5 10 15 20 Obrázek 25: Výběrová ACF pro simulovaná data ARMA procesu Yf — 0.9Yí_i + 0.5Yt_2 = et + 0.75£í_i + 0.25£i_2 > tARMApaci <- ARMAaci(ar = ARphi, ma = MAtheta, lag.max = 20, paci = TRUE) > par (mirow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot(l:length(tARMApacf), tARMApaci, type = "h", main = "Teoretická PACF") > abline(h = 0) Stochastické modely časových řad (6. cvičení) 23 Teoretická PACF -1-1-1-r 5 10 15 20 Obrázek 26: Teoretická PACF pro ARMA proces Yt - 0.9Yt_i + 0.5^-2 = et + 0.75et_i + 0.25£í_2 > vARMApacf <- pacf (arma. si m, lag.max = 20, plot = FALSE) > par (mirow = c(l, 1), mar = c (2, 2, 3, 0) + 0.05) > plot (vARMApacf, main = "Výberová PACF") Výberová PACF -1-1-1-r 5 10 15 20 Obrázek 27: Výběrová PACF pro simulovaná data ARMA procesu Yt — 0.9Yí-i + 0.5YÍ-2 = Et + 0.75£í_i + 0.25£í_2 Pomocí funkce arima() provedeme odhad parametrů ARMA(2, 2) modelu. Do argumentu order musíme zadat řád odhadovaného ARIMA modelu. ARMA proces v našem příkladu je speciálním případem ARIMA procesu s řádem c(2,0,2). > arma.fit <- arima(arma.sim, order = c(2, 0, 2)) > print (arma. f it) Series: arma.sim ARIMA(2,0,2) with non-zero mean Call: arima(x = arma.sim, order = c(2, 0, 2)) Coefficients: 24 Stochastické modely časových řad (6. cvičení) arl ar2 mal ma2 intercept 0.9072 -0.5185 0.6644 0.2077 -0.1461 s.e. 0.0944 0.0723 0.1045 0.0923 0.1710 sigma~2 estimated as 0.9377: log likelihood = -417.41 AIC = 846.82 AICc = 847.11 BIC = 869.04 Kvalitu vyrovnání dat pomocí právě odhadnutého modelu posoudíme na základě reziduí s využitím funkce tsdiagO, kterou aplikujeme na odhadnutý ARMA(2,2) model. > tsdiagCarma.fitJ Standardized Residuals ihJli.nl li Li Jj J li il.li ll II hli ll .d i hLi ii .ii lili Mi. ,i i .u i w II1 ľ pi" 111 .....i i|'|ľľ| 1 ||||ľ'i 'l|l ľ1 "'l'ľ'l'í i'| "ľ ľ ji1 mi pH m in íl ľ T-1-1-1-1-1-1- 0 50 100 150 200 250 300 Time ACF of Residuals - 1 1 ' JI ' T-1-1-1-T 0 5 10 15 20 Lag p values for Ljung-Box statistic ° -1-1-1-1-1- 2 4 6 8 10 lag Obrázek 28: Diagnostické grafy pro odhadnutý ARMA(2,2) model pro simulovaná data ARMA procesuj - 0.9^-1 + 0.5yt_2 = et + 0.75et_i + 0.25et_2 Z grafů opět vidíme, že po odstranění systematické složky z časové řady dostaneme pouze bílý šum. Na základě odhadnutého modelu provedeme predikci pro 10 budoucích hodnot. K tomu opět použijeme funkci PlotPredictARIMA. Jako argumenty zadáme zkoumanou časovou řadu (arma.sim), odhadnutý model pro tuto řadu (arma.fit) a počet predikovaných hodnot (n.ahead > par (mfrow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > PlotPredictARIMA(ajrma.sim, arma.fit, n.ahead = 10) Stochastické modely časových řad (6. cvičení) 25 6 - O 50 100 150 200 250 300 Obrázek 29: Predikce na základě odhadnutého ARMA(2,2) modelu pro simulovaná data ARMA procesu Yt - 0.9Yt_i + 0.5yt_2 = et + 0.75et_i + 0.25et_2 Ještě se podívejme, jak odhadne řád modelu a parametry procesu funkce auto.arima. > arma.fitF <- auto. arima(arma. sim, tra.ce = TRUE) ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA ARIMA (2,0,2 (0,0,0 (1,0,0 (0,0 (1,0 (3,0 (2,0 (2,0 (1,0 (3,0 (2,0 (1,0 (3,0 (2,0 (2,0 (1,0 (3,0 with with with with with with with with with with with with with with with with with non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean non-zero mean zero mean zero mean zero mean zero mean zero mean zero mean zero mean 841. 6855 1332.781 1109.331 1045.133 883. 1501 843. 9904 845. 0074 843. 5533 929. 7609 845. 8741 840. 1418 881. 8383 842. 563 843. 601 841. 9994 928. 0157 844. 4304 Best model: ARIMA(2,0,2) with zero mean > priu t(arma.fi tF) Series: arma.sim ARIMA(2,0,2) with zero mean Call: auto.arima(x = arma.sim, trace = TRUE) Coefficients: arl ar2 mal ma2 0.9057 -0.5157 0.6681 0.2104 s.e. 0.0945 0.0726 0.1044 0.0921 sigma~2 estimated as 0.94: log likelihood = -417.77 AIC = 845.55 AICc = 845.75 BIC = 864.06 26 Stochastické modely časových řad (6. cvičení) Na úplný závěr provedeme predikci pro 10 budoucích hodnot na základě posledního modelu, přičemž použijeme funkci forecast(). > par (mírow = c(l, 1), mar = c (2, 2, 1, 0) + 0.05) > plot (forecast (arma. íitF) , n = 10) _Forecasts from ARIMA(2,0,2) with zero mean CD - 0 50 100 150 200 250 300 Obrázek 30: Predikce na základě odhadnutého APJMA(2,0,3) modelu pro simulovaná data ARMA procesu Yt - 0.9Yt_i + 0.5yt_2 = et + 0.75£t_i + 0.25^-2 5 Úkol: 1. Rozhodněte o kauzalitě, resp. invertibilitě následujících procesů: (a) Yt = 0.75Yí_i + 0.5yt_2 + et (b) Yt = et + 0.75£í_i + 0.75£í_2 (c) Yt = 0.8Yí_i - 0.2yt_2 + et + 0.6£í_i - 0.2£í_2 + 0.1et_3 2. Pro procesy z 1. úkolu, u nichž jste ověřili, že jsou kauzální (u AR procesu), invertibilní (u MA procesu) nebo zároveň kauzální i invertibilní (u ARMA procesu) vygenerujte 300 hodnot časové řady, vykreslete autokorelační a parciální autokorelační funkci a odhadněte parametry modelu např. s využitím funkce arima.