Stochastické modely časových řad (11. cvičení) 1 M5201 - 11. CVIČENÍ: Ověřováni vhodnosti modelů pomocí analýzy reziduí a vhodných testů 1 Míry influence v regresní diagnostice V praxi se často můžeme setkat s jevem, že v souboru dat se vyskytují některé hodnoty výrazně se lišící od hodnot ostatních. V literatuře se v průběhu minulých let rozvinuly dva směry, které se svým způsobem snaží s jejich existencí vyrovnat. • metody robustní statistiky • metody regresní diagnostiky Regresní diagnostika jde cestou detekovat více či méně ojedinělá data a dát původci dat možnost rozhodnout se, jak s nimi v případě výskytu dále naložit, tj. zdaje v souboru ponechat či vyloučit, věnovat jim menší váhu při zpracování, popřípadě je vhodně transformovat apod. V rámci regresní diagnostiky se budeme zabývat dvěma základními úlohami • jak detekovat mezi daty neočekávané hodnoty • jak rozhodnout, zda mohou významně ovlivnit statistickou analýzu, případně jakým způsobem. Mějme regresní model Y = X/3 + e A Es = 0 A vare = a2In A h(X.) = k = p + 1 Vyrovnané hodnoty na základě odhadů metodou nejmenších čtverců dostaneme ze vztahu: Ý = X/3 = X(X'X)_1X'Y. H Projekční matice H X(X'X)_1X' je idempotentní H: Označme matici H symetrická h\n \ H' h. nn (Hi, • • • , : H H . Pak Y = HY. Díky symetrii platí: Yt = H^Y = ^2 hijYj a odtud vidíme, jak pozorování Y j ovlivňuje skrze híji-tou vyrovnanou hodnotu Yj. Definice 1: Sloupce matice H se nazývají vlivové vektory a ||Hj ||2 se nazývá vliv j-tého pozorování na odhad Y, stručně j-tý vliv. Věta 1: |H,||2 2 Stochastické modely časových řad (11. cvičení) Věta 2: Vľ=i 0 < hn < 1 Věta 3: írH = Věta 4: Průměrný vliv pozorování Y±, ■ ■ ■ , Yn je roven Definice 2: Definujme vektor reziduí: r = Y-Y Věta 5: r = (I - H)e Pokud jsou prvky ha ~ 1 (blízké k 1) =>• 1 — ha ~ 0. Pak neočekávaně velká chyba pozorování Yi (velká chyba a) se nemusí odrážet v r i. Na ostatní rezidua však vliv mít může. Věta 6: DY = cj2H Dr = a2(I - H) Pokud matice H není diagonální maticí, pak rezidua jsou korelovaná. Z předchozího je vidět, že rezidua r\ v některých situacích nemusí dobře identifikovat odlehlá pozorování. Proto se v literatuře zavádějí a používají další typy reziduí. Značení: Symbol (i) bude znamenat vynechání i-tého řádku, symbol [j] vynechání j-tého sloupce. Definice 3: Definujme: ■ standardizovaná rezidua (také normovaná rezidua) s^/l-hu kde s2 (Y-Y)'(Y-Y) n—k SSE n—k ' i-té predikované reziduum Yi-Y, (0 kde v lineárním modelu vynecháme i- té pozorování a značíme matici plánu X(j), vektor Y (j), odhad metodou nejmenších čtverců $^ = (X^-jX^-^^X^) Y(j) a i-té pozorování odhadneme pomoci $^ takto Y"(j) = XjjS^-j, kde Xj je i-tý řádek původní matice plánu. „■— /i _ h.. kde studentizovaná (jackknife) rezidua 2 _ (Y(i)-Y(i))/(Y(i)-Y(i)) ~ n-k-1 i-té DFFIT reziduum cL (0 r(i)Vl-hu yi-y(i) \fhiis(i) i-té parciální reziduum %1 Y-Y, kde v lineárním modelu vynecháme j- tý regresor, odpovídající matici plánu označíme Xyj, odhad metodou nejmenších čtverců 0yj = (Xjj.jX[J-])_1X[J-]Y a i-té pozorování odhadneme pomoci takto Yí\j\ = xi[j]0[j]> kde x^-j je i-tý řádek matice plánu Xyj. Stochastické modely časových řad (11. cvičení) 3 Věta 7: Nechť e ~ Nn(0, a2In) r ~ Nn(0, a2{I - H)) Protože s2(l — /ijj) je odhad rozptylu Dri = ^(l — hn), pak standardizovaná rezidua mají rozptyl přibližně roven 1. Pokud nastane situace, že chyba je příliš velká oproti modelu, pak pomocí příslušného standardizovaného rezidua je lze snadněji identifikovat. Věta 8: Věta 9: Věta 10: r(i) l-hii Dr « 1-hu (n - k)s2i} = (n- k)s2 - j^- W Vl-hus(i) Předchozí vzorce umožňují vypočítat statistiky r^, s2^ a pouze z hodnot známých z celého regresního modelu. Věta 11: Nechť e ~ Nn(0, a2Ir r*(i) ~ ř(n - k) Definice 4: Řekneme, že i-té pozorování Y;t je odlehlé, jestliže Eeí 7^ 0 Z věty 11 plyne, že pomocí i-tého studentizovaného rezidua lze testovat nulovou hypotézu Hq : Eeí = 0, že i-té pozorování není odlehlé proti alternativě H\ : Eeí 7^ 0, tj. že je odlehlé. Pokud \r\i)\ — ti-a/2(n ~ k)> Pak na hladině významnosti a zamítáme hypotézu Hq a tedy i-té pozorování na hladině významnosti a je odlehlé. Věta 12: DFFITS rezidua di pro i = 1, - • • ,n lze vyjádřit vztahem 1-hu Je-li n — k > 30, je na hladině významnosti a = 0.05 kvantil Studentova t rozdělení přibližně roven 2 a lze tedy v praktických situacích na základě věty 11 považovat i-té pózování za odlehlé na hladině významnosti a = 0.05, když | > 2 . Toto odpovídá také empirických zkušenostem (viz. Staudte, Sheather(1990)). Posuzujeme-li odlehlé pozorování pomocí i-tého DFFITS rezidiua, můžeme využít vztah z věty 12 a navíc uplatnit vliv i-tého pozorování ha a jeho průměrnou hodnotu. Pak platí \di 1 - hí |r«l>2lr^- Uvážíme-li, že průměrný vliv je - a dosadíme-li jej za ha, dostaneme 1 k \di 1 1 1. 2 /k\ 2 n-kl >2 U. Posledně uvedená nerovnost se v praxi užívá pro posouzení, zda i-té pozorování je odlehlé na základě DFFIT reziduí. 4 Stochastické modely časových řad (11. cvičení) Cookova vzdálenost. Pro měření vlivu i-tého pozorování na hodnotu odhadu vektoru f3 navrhl Cook použít statistiku = (Y-Y(i))'(Y-Y(i)) = r$ hjj Cookova vzdálenost souvisí s konfidenčním elipsoidem odhadů, což umožňuje její porovnání s kvantily F-rozdělení s k a n — k stupni volnosti. Jde zde však o posun odhadů, který vznikl vynecháním i-tého bodu. Orientačně platí, že pro Dí > 1 posun přesahuje 50%ní konfidenční oblast a daný bod je proto vlivný. Další možné vysvětlení Cookovy vzdálenosti vychází z toho, že jde o eukleidovskou vzdálenost mezi vektorem predikce Y z metody nejmenších čtverců a vektorem predikce Y(j), který odpovídá odhadům stanoveným metodou nejmenších čtverců při vynechání i-tého bodu. Cookova vzdálenost vyjadřuje vliv i-tého bodu pouze na odhady parametrů f3. Pokud proto i-tý bod neovlivní odhady regresních parametrů f3 výrazně, bude hodnota Cookovy vzdálenosti malá. Takový bod však může silně ovlivnit odhad reziduálního rozptylu a2, kde De = a2In. Welschova-Kuhova vzdálenost. Pro měření vlivu i-tého pozorování simultánně jak na odhad (3, tak na odhad a2, zvolili Welsch a Kuh statistiku DFFITSi = di V/Kíls(í) Parciální vliv. Pro měření vlivu i-tého pozorování na j-tou složku vektoru 0 navrhli Belsley, Kuh a Welsch statistiku DFBETAS,; -Pj(.i) Variační poměr. Pro stanovení míry vlivu i-tého pozorování na matici D0 je navržena statistika C OV RATIO,, l-hu Velké hodnoty statistiky signalizují vlivné body (vzhledem k D0). V literatuře se za vlivná pozorování doporučuje považovat ta, pro něž \C0VRATI0j - II > 3Í LITERATURA: Anděl, J.(1978): Matematická statistika, Praha, SNTL. Antoch, J., Vorlíčková, D. (1992): Vybrané metody statistické analýzy dat, Academia Praha Staudte, R.G.,Sheather, S.J.(1990): Robust Estimation and Testing, New York, Wiley Stochastické modely časových řad (11. cvičení) MÍRY INFLUENCE V PROSTŘEDÍ R: K dispozici je především funkce influence .measuresQ kterou lze použít pouze na objekt třídy lm. Výsledkem je objekt třídy infl, který je tvořen ze tří prvků: infmat, is.inf a call. Matice infmat : • každý řádek odpovídá jednomu pozorování (Yj,Xj) • prvních k sloupců tvoří matici statistik DFBETAS, sloupce jsou označeny dbf. a za tečkou následuje většinou přiměřeně zkrácený název příslušného regresoru • následuje sloupec statistik DFFITS označený dff it. • další sloupce nazvané cov.r, cook.d hat obsahují statistiky COVRATIO, Dí a diagonální prvky matice H. Samostatně lze jednotlivá rezidua a další statistiky získat z objektu typu lm pomocí funkcí rstandardO dffitsO dfbetasO covratioO rstudentO cooks. distance () hatvaluesO PŘÍKLAD 1: Koncentrace oxidu uhličitého v okolí sopky Mauna Loa V datovém souboru C02.dat jsou obsažena měsíční data, která udávají koncentraci oxidu uhličitého v okolí havajské sopky Mauna Loa od ledna 1965 do prosince 1980. Na prvním řádku souboru je uveden popis dat a data jsou pod ním ve dvanácti sloupcích. Nejjednodušší způsob, jak datový soubor načíst, je použít příkaz scan(). Při načítání bude nutné přeskočit první řádek, tudíž přidáme argument skip=l. > íileDat <- paste (data, library, "C02.dat", sep = "") > co2 <- scan(ÍileDat, skip = 1) Příkazem scan jsme data načetli do vektoru. Pro další práci bude výhodné data převést na časovou řadu (pomocí ts()). U nově vytvořené časové řady si prohlédneme její strukturu. > TS <- ts(co2, start = 1965, frequency = 12) > str(TS) Time-Series [1:192] from 1965 to 1981: 319 320 321 322 322 ... Časovou radu vykreslíme pomoci príkazu plot(). > par (mar = c (2, 2, 1, 0) + 0.5) > plot(TS, type = "o", pch = 20, cex = 3) 6 Stochastické modely časových řad (11. cvičení) 1965 1970 1975 1980 Obrázek 1: Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 Protože data hned na první pohled jasně vykazují sezónní chování, bude třeba k jejich dekompozici použít model zahrnující sezónnost. Pro ten účel zvolíme například (modifikovanou) metodu malého trendu Yjk = 11 + rrij + sk + ejk ejk ~ WN(0, a2), j = 1,..., r, k = l,...,d s dodatečnými podmínkami s± + • • • + Srf = 0 a ni\ + • • • + mr = 0. Dekompozici provedeme pomocí funkce SzSmallTrendModif (). > íileFun <- paste (data.library, "FunkceM5201.R", sep = "") > source(ÍileFun) > vysl <- SzSmallTrendModif(TS) Stochastické modely časových řad (11. cvičení) 7 i i i r 1965 1970 1975 1980 Obrázek 2: Metoda malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980. Pomocí funkcí summary () a anova() vypíšeme výsledný regresní dekompoziční model s příslušnými statistikami. > summary(vysl$model) Call: lm(formula = x grY + grS, data = data, contrasts = list(grY = contr.sum, grS = contr.sum)) Residuals: Min IQ Median 3Q Max -0.67042 -0.17260 -0.00146 0.18635 0.81125 Coefficients: ] (Intercept) 3'. grYl grY2 grY3 grY4 grY5 grY6 grY7 grY8 grY9 grYlO grYll grY12 grY13 grY14 imate Std. Error t value PrOltl) 46396 0 02054 15989 952 <2e- 16 *** 47229 0 07956 -106 491 <2e- 16 *** 80146 0 07956 -98 059 <2e- 16 *** 99063 0 07956 -87 868 <2e- 16 *** 16646 0 07956 -77 509 <2e- 16 *** 28229 0 07956 -53 826 <2e- 16 *** 95729 0 07956 -37 171 <2e- 16 *** 99312 0 07956 -25 052 <2e- 16 *** 86979 0 07956 -10 933 <2e- 16 *** 35437 0 07956 17 024 <2e- 16 *** 93604 0 07956 24 335 <2e- 16 *** 53354 0 07956 31 845 <2e- 16 *** 59854 0 07956 45 231 <2e- 16 *** 16354 0 07956 64 903 <2e- 16 *** 73187 0 07956 84 615 <2e- 16 *** 8 Stochastické modely časových řad (11. cvičení) ;rY15 8. 08437 0. 07956 101. 616 <2e-16 *** ;rSl -0. 63458 0. 06813 -9. 314 <2e-16 *** ;rS2 0. 08292 0. 06813 1. 217 0.225 ;rS3 0. 95104 0. 06813 13. 959 <2e-16 *** ;rS4 2. 10542 0. 06813 30. 903 <2e-16 *** ;rS5 2. 62667 0. 06813 38. 554 <2e-16 *** ;rS6 2. 22292 0. 06813 32. 628 <2e-16 *** ;rS7 0. 93667 0. 06813 13. 748 <2e-16 *** ;rS8 -1. 00458 0. 06813 -14. 745 <2e-16 *** ;rS9 -2. 63333 0. 06813 -38. 652 <2e-16 *** ;rS10 -2. 74208 0. 06813 -40. 248 <2e-16 *** ;rSll -1. 51458 0. 06813 -22. 231 <2e-16 *** Signif. codes: 0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.2846 on 165 degrees of freedom Multiple R-squared: 0.998, Adjusted R-squared: 0.9977 F-statistic: 3217 on 26 and 165 DF, p-value: < 2.2e-16 > anova(vysl$model) Analysis of Variance Table Response: x Df Sum Sq Mean Sq F value Pr(>F) grY 15 6195.3 413.02 5097.88 < 2.2e-16 *** grS 11 582.1 52.91 653.12 < 2.2e-16 *** Residuals 165 13.4 0.08 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Nyní, když máme pro data odhadnutý model, bychom chtěli ověřit, zda se pro naše data skutečně hodí, a také to, zda některá data nemají neobvykle velký vliv na odhady parametrů modelu. Proto postupně provedeme analýzu reziduí. Nejprve začneme tím, že se pomocí regresní diagnostiky pokusíme najít případná odlehlá a vybočující pozorování. Základní diagnostické grafy dostaneme v R tak, že na odhadnutý model zavoláme funkci plot. lm(). Stačí psát zkrácenou verzi plot (). Příkaz vykresluje šest různých diagnostických grafů, zajímá-li nás pouze některý z nich, použijeme argument which=" číslo grafu". Prohlédněme si nejprve první dva grafy. > par(mfrow = c(l, 2), mar = c(5, 5, 2, 0) + 0.05) > plot(vysl$model, which = 1:2) Stochastické modely časových řad (11. cvičení) 9 Residuals vs Fitted Normal Q-Q CÖ ■g C/3 CD ö o ö lo ö i 320 330 Fitted values 340 CC ■g C/í CD L_ T3 CD N T3 l_ CC T3 C 05 CO CO — C\J - 190^9 o-'' O - CNJ 1 _ 1 3 I I I I -10 1 2 2 Theoretical Quantiles Obrázek 3: Analýza reziduí pomocí funkce plot () - grafy 1 a 2 u metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980. Na prvním grafu jsou znázorněny na rr-ové ose vyrovnané hodnoty a na ose y proti nim příslušná rezidua. Rezidua by neměla být nijak závislá na velikosti odpovídajících vyrovnaných hodnot, tudíž by neměla vykazovat žádný trend a měla by kolísat kolem nuly. V grafu je zobrazena pomocná červená linie, která představuje neparametrický odhad střední hodnoty reziduí. První diagnostický graf pro naše rezidua vypadá přibližně tak, jak bychom očekávali. Druhým diagnostickým grafem je Q-Q plot, který slouží k vizuálnímu posouzení pravděpodobnostního rozložení zkoumaných dat. U reziduí nás obvykle zajímá, zda mají normální rozložení. V tom případě by body Q-Q plotu (konstruovaného pro normální rozdělení) měly ležet přibližně na přímce. Podíváme-li se na Q-Q plot pro naše rezidua, vidíme, že body na přímce zhruba leží, takže nic nenasvědčuje tomu, že by rezidua nemohla mít normální rozdělení. Prohlédněme si další dva diagnostické grafy. > par (mírow = c(l, 2), mar = c (5, 5, 2, 0) + 0.05) > plot(vysl$model, which = 3:4) 10 Stochastické modely časových řad (11. cvičení) Scale-Location 320 330 340 Fitted values Obrázek 4: Analýza reziduí pomocí funkce plc data Koncentrace CO2 u sopky Mauna Loa v Cook's distance o 1 0 50 100 150 Obs. number O — grafy 3 a 4 u metody malého trendu pro tech 1965-1980 Na třetím grafu vidíme odmocniny ze standardizovaných reziduích v závislosti na vyrovnaných hodnotách. Opět by mělo platit (podobně jako u prvního grafu), že by odmocniny ze standardizovaných reziduí neměly vykazovat závislost na vyrovnaných hodnotách. Čtvrtým grafem je graf Cookových vzdáleností Y od . Čím větší má některé pozorování Cookovu vzdálenost, tím větší je jeho vliv na odhady regresních koeficientů. V grafu jsou čísly označena pozorování 22, 25 a 190 s největšími Cookovými vzdálenostmi, která se nabízí jako vlivná pozorování. Z grafu však vidíme, že největší hodnota Cookovy vzdálenosti je 0,06, což je velmi malá hodnota v porovnání s hodnotou 1, která se orientačně používá jako dolní mez pro vlivná pozorování. Ještě se podívejme na pátý diagnostický graf, který znázorňuje závislost reziduí na faktorech vysvětlující proměnné. > par (mírow = c(l, 1), mar = c (5, 5, 2, 0) + 0.05) > plot(vysl$model, which = 5) Stochastické modely časových řad (11. cvičení) 11 Constant Leverage: Residuals vs Factor Levels 1965 1970 1975 198C Factor Level Combinations Obrázek 5: Analýza reziduí pomocí funkce plot () - graf 5 u metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 Další zajímavé grafy lze najít v knihovně car (příkazy inf luencePlot (), inf IndexPlot ()). První z nich je určen spíše pro klasickou regresi, druhý je již vhodnější pro naše účely, proto ho použijeme na náš model. Argumentem vars můžeme určovat, které grafy chceme vykreslit. Zvolíme graf Cookových vzdáleností, studentizovaných reziduí, Bon-ferroniho p-hodnot a vlivů ha. > library(car) > par (mírow = c(l, 1), mar = c (5, 5, 2, 0) + 0.05) > infIndexPlot(vysl$model, vars = c("Cook", "Studentized", "Bonf", "hat")) 12 Stochastické modely časových řad (11. cvičení) Diagnostic Plots Index Obrázek 6: Analýza reziduí pomocí funkce inf IndexPlot () - metoda malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980. První graf s Cookovými vzdálenostmi je prakticky totožný se čtvrtým diagnostickým grafem funkce plot. lm(). Na druhém grafu vidíme studentizovaná rezidua. Pokud je studentizované reziduum v absolutní hodnotě větší než 2, jemu odpovídající pozorování se považuje za odlehlé. V grafu najdeme takových reziduí hned několik (přibližně 10). Mezní hodnota 2 pro posouzení odlehlosti pozorování na základě studentizovaných reziduí je odvozena na základě testu s 5% hladinou významnosti, z jehož povahy vyplývá, že vždy 5% největších reziduí přesáhne tuto hranici. Bude se sice jednat o 5% „nejextrémnějších" pozorování, ale neznamená to, že všechna jsou skutečně výjimečná. Abychom zjistili, která z nich jsou skutečně odlehlá, použijeme Bonferroniho test. Jeho p-hodnoty jsou ve třetím grafu. Pokud je p-hodnota menší než zvolená hladina významnosti (obvykle 0,05), považuje se odpovídající pozorování za odlehlé. Z grafu vidíme, {z""!" ""!" ""!" ""!" ""!" ""!" ""!" ""!" ""!" "!!" ""!" ""! Stochastické modely časových řad (11. cvičení) 13 že p-hodnota není u žádného pozorování menší než 0, 05, tudíž žádné pozorování není skutečně odlehlé. V posledním grafu jsou vykresleny vlivy hu (diagonální prvky projekční matice H = X(X'X)_1X'). Ideálně by měly být vlivy všech pozorování přibližně stejné a pohybovat se kolem k/n. V grafu hledáme hlavně pozorování s neobvykle velkým vlivem. Pro náš model bychom takové pozorování v grafu našli těžko. 2 Testování normality V klasickém regresním modelu se předpoládá, že náhodné odchylky e i mají normální rozdělení. Existuje celá řada testů zabývající se normalitou. Jedna skupina testů je založena na empirických distribučních funkcích, jako zástupce můžeme uvést Kolmogorův-Smirnovův test, popř. Shapiro-Wilkův test. Další testy jsou založeny na momentových charakteristikách, především na šikmosti či špičatosti. Příkladem může být ďAgostinův test, popř. Jarque-Bera test. V základním balíku R base najdeme dva známé testy normality: Shapiro-Wilkův test a Kolmogorův-Smirnovův test. 2.1 Shapiro Wilkův test pro testování normality Shapiro-Wilkův test je založen na statistice kde jsou pořádkové statistiky a a;t jsou váhy, které jsou odvozeny ze středních hodnot a varianční matice pořádkových statistik prostého náhodného výběru z N(0,1) rozsahu n. Tyto hodnoty bývají tabelovány. Hlavní myšlenka stojící za testem spočívá v tom, že zjišťujeme, zda body v Q-Q plotu pro testovaná data jsou významně odlišné od regresní přímky proložené těmito body. Testová statistika dosahuje hodnoty 1 v případě, že data vykazují perfektní shodu s normálních rozdělením. Je-li W statisticky významně nižší než 1, zamítáme nulovou hypotézu o shodě s normálním rozdělením. V knihovnách tseries, popř. FitAR lze najít tzv. Jarque-Bera test normality, který je založen na výběrové šikmosti a špičatosti. 2.2 Jarque-Bera test pro testování normality Jedná se o test založený na porovnání výběrové šikmosti a špičatosti s teoretickou šikmostí a špičatosti pro normální rozdělení. W n _ i=i 14 Stochastické modely časových řad (11. cvičení) Označme výběrový průměr X = M\ = - X,t TI . výběrový k-tý centrální moment Mk = - J2 — X) í=i výběrová šikmost B\ m3 výběrová špičatost B2 = Pak pro Jarque-Bera statistiku platí PŘÍKLAD 1 (POKRAČOVÁNÍ): Koncentrace oxidu uhličitého v okolí sopky Mauna Loa Testy normality budeme provádět na standardizovaných reziduích modifikovaného modelu malého trendu. Ta získáme příkazem rstandard(). Pro Shapiro-Wilkův test normality standardizovaných reziduí použijeme funkci shapiro. test (), Jarque-Berův test provedeme příkazem jarque.bera.test(). > library(tseries) 'tseries' version: 0.10-25 'tseries' is a package for time series analysis and computational finance. See 'library(help="tseries")' for details. > res.standard <- rstandard(vysl$model) > shapiro.test(res.standard) Shapiro-Wilk normality test data: res.standard W = 0.9934, p-value = 0.5531 > jarque.bera.test(res.standard) Jarque Bera Test data: res.standard X-squared = 0.5919, df = 2, p-value = 0.7438 Stochastické modely časových řad (11. cvičení) 15 U prvního testu jsme dostali p-hodnotu 0,55, u druhého dokonce 0,74. Obě čísla jsou výrazně větší než hladina významnosti testu 0, 05, tudíž nezamítáme hypotézu o tom, že standardizovaná rezidua mají normální rozdělení. Rozložení standardizovaných reziduí si můžeme prohlédnout, když si vykreslíme jejich histogram. Když k vykreslení použijeme funkci HistFit (), vykreslí se nám do grafu kromě samotného histogramu také jádrový odhad hustoty reziduí a hustota normálního rozdělení se střední hodnotou a rozptylem, které dostaneme jako odhad střední hodnoty a rozptylu reziduí. > par (mírow = c(l, 1), mar = c (2, 2, 2, 0) + 0.05) > HistFit (res. standard) Histogram, Kernel Density Estimate, Normal Curve -3-2-10 1 2 3 x Obrázek 7: Grafické testování normality u standardizovaných reziduí metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 3 Testování homoskedasticity rozptylu K testování homoskedasticity rozptylu se často používá Breusch-Paganův test. Použijeme variantu vhodnou pro časové řady. Jestliže jsou rezidua homoskedastická, jejich rozptyl jev čase konstantní. Homoskedas-ticita je porušena například v případě, kdy rozptyl reziduí o\ s časem lineárně roste nebo klesá. Breusch-Paganův test je založen na myšlence popsat variabilitu reziduí pomocí regresního modelu ve tvaru log^t bt a následně testovat hypotézu 0 v s Pokud zamítneme nulovou hypotézu, prokázali jsme, že rozptyl reziduí není konstantní a tudíž jsou rezidua heteroskedastická. 16 Stochastické modely časových řad (11. cvičení) PŘÍKLAD 1 (POKRAČOVÁNÍ): Koncentrace oxidu uhličitého v okolí sopky Mauna Loa Testování homoskedascity budeme opět provádět na standardizovaných reziduích modifikovaného modelu malého trendu. Breush-Paganův test homoskedaticity lze v R provést příkazem bptestO, který najdeme v balíčku lmtest. > library(lmtest) > Time <- 1 :length(res. standard) > bptest(res. standard ~ Time) studentized Breusch-Pagan test data: res.standard ~ Time BP = 0.0191, df = 1, p-value = 0.8902 Homoskedasticita rozptylu nebyla zamítnuta, protože p-hodnota testu 0,89 je větší než hladina významnosti 0, 05. Podívejme se na problém homoskedasticity ještě s využitím grafu. Použitím příkazu boxplot Segment s () dostaneme krabicové grafy pro krátké úseky dat o čtyřech pozorováních. Z jejich rozměrů lze vidět, jak se variabilita reziduí mění v čase. > par (mírow = c(l, 1), mar = c (5, 2, 2, 0) + 0.05) > boxplotSegments(res.standard, seglen = 4, xlab = "Standardised Residuals") _l_ ) Úŕ n ^T ě w W u y x x ± J X X 1 ! I x ?! 1 4 7 10 14 18 22 26 30 34 38 42 46 Standardised Residuals Obrázek 8: Grafické testování homoskedasticity u standardizovaných reziduí metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 4 Testování nezávislosti reziduí Pokud je regresní model vyhovující, rezidua by měla být přibližně normálním bílým šumem. Stochastické modely časových řad (11. cvičení) 17 K testování bílého šumu se používají tzv. Portmanteauovy statistiky, nejčastěji Ljung-Boxova nebo Box-Piercova statistika. Mají přibližně \2 rozdělení a jsou založené na vztazích BP = n E r2 LB = n(n + 2) E ^- i=l PŘÍKLAD 1 (POKRAČOVÁNÍ): Koncentrace oxidu uhličitého v okolí sopky Mauna Loa Testování reziduí jako bílého šumu budeme provádět na standardizovaných reziduích modifikovaného modelu malého trendu. K provedení Box-Piercova testu je v R funkce Box.test(). Chceme-li provést Ljung-Boxův test, použijeme tutéž funkci a jako její argument uvedeme type="Ljung-Box". > Box.test(res.standard) Box-Pierce test data: res.standard X-squared = 13.8135, df = 1, p-value = 0.0002019 > Box. test (res. standard, type = "Ljung-Box") Box-Ljung test data: res.standard X-squared = 14.0305, df = 1, p-value = 0.0001799 Oba testy dávají jako výsledek p-hodnoty výrazně menší než hladina významnosti 0, 05, a proto v obou případech zamítáme nulovou hypotézu, že rezidua jsou bílým šumem. Tento výsledek by se měl rovněž projevit v grafech autokorelační a parciální autokorelační funkce. > par(mírow = c(2, 1), mar = c(2, 2, 1, 0) + 0.05) > acf(res.standard) > mtext ("ACF") > acf (res. standard, type = "partial") > mtext ("PACF") 18 Stochastické modely časových řad (11. cvičení) 0.2 0.6 1.0 ...... lil I . . CM o -I M' 1 I 1 | ---------------|---- I 0 1 1 1 1 5 10 15 20 PACF CM O O CM O - I Obrázek 9: Grafické testování bílého šumu u standardizovaných reziduí metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 Pokud by byla rezidua bílým šumem, měla by ACF významně nenulovou hodnotu pouze pro fc = 0a PACF by měla všechny hodnoty přibližně nulové. 5 Autokorelace reziduí V regresních modelech pro časové řady je třeba věnovat velkou pozornost problematice autokorelovaných reziduí. Ve většině případů se u časových řad s autokorelací reziduí setkáme, neboť hodnota pozorování v časovém okažiku t velmi pravděpodobně ovlivní následující hodnoty. Pro testování autokorelace reziduí prvního řadu je používán Durbin-Watsonův test 5.1 Durbin Watsonův test autokorelace reziduí 1. řádu Durbin-Watsonova statistika je definována vztahem n D = ^~n-• En2 1=1 Lze odvodit, že statistika D může nabývat pouze hodnot z intervalu (0,4). Pokud budou rezidua málo korelovaná, hodnota D se bude pohybovat kolem 2. Kladná korelace způsobí, že D E (0, 2) a záporná korelace způsobí, že D E (2,4). Přesné rozdělení statistiky D závisí na tvaru matice plánu X, proto jsou tabelovány intervaly a dy, ve kterých se nachází kritické hodnoty (pro různá n, k a a). Stochastické modely časových řad (11. cvičení) 19 Dolní a horní hranice Durbin-Watsonova testu na 5% hladině významnosti k= =1 k= =2 k= =3 =4 k= =5+ n dh du dL du dL du dL du dL du 50 1.50 1.59 1.46 1.63 1.42 1.67 1.38 1.72 1.34 1.77 60 1.55 1.62 1.51 1.65 1.48 1.69 1.44 1.73 1.41 1.77 70 1.58 1.64 1.55 1.67 1.52 1.70 1.49 1.74 1.46 1.77 80 1.61 1.66 1.59 1.69 1.56 1.72 1.53 1.74 1.51 1.77 90 1.63 1.68 1.61 1.70 1.59 1.73 1.57 1.75 1.54 1.78 100+ 1.65 1.69 1.63 1.72 1.61 1.74 1.59 1.76 1.57 1.78 kde k je počet nezávisle proměnných v regresní rovnici. Pro rychlé posouzení autokorelace prvního řadu vystačíme s následující tabulkou: Pokud hodnota Durbin-Watsonovy statistiky D bude v mezích 0 až dL dL až du du až (4 — du) (4 - du) až (4 - dL) (4 - dL) až 4 Zamítáme H0 KLADNÁ autokorelace Ani nezamítáme ani nepřijímáme Ho Nezamítáme nulovou hypotézu Ho Ani nezamítáme ani nepřijímáme Ho Zamítáme Ho NEGATIVNÍ autokorelace PŘÍKLAD 1 (POKRAČOVÁNÍ): Koncentrace oxidu uhličitého v okolí sopky Mauna Loa Testování autokorelace reziduí budeme provádět na standardizovaných reziduích modifikovaného modelu malého trendu. Durbin-Watsonův test provedeme příkazem dwtest () z baličku lmtest. Při použití funkce musíme zadat, že chceme testovat autokorelaci reziduí v závislosti na čase, oboustrannou alternativní hypotézu zadáme argumentem alternative="two.sided". > library(lmtest) > Time <- 1 :length(res. standard) > (DWtest <- dwtest(res.standard ~ Time, alternativě = "two.sided")) Durbin-Watson test data: res.standard ~ Time DW = 1.4615, p-value = 0.0001313 alternativě hypothesis: trne autocorelation is not 0 Vidíme, že p-hodnota je výrazně menší než hladina významnosti 0,05, takže se prokázala autokorelace reziduí prvního řádu. Abychom lépe pochopili, co představuje autokorelace 1. řádu, je třeba si uvědomit, že pokud pro rezidua (standardizovaná s nulovou střední hodnotou) platí autokorelace prvního řádu, pak můžeme psát r{ = ifn-x + rji, kde i]i jsou nekorelované náhodné veličiny s nulovou střední hodnotou. Proto použijeme k prozkoumání tohoto vztahu jednoduchý regresní model. I když předchozí regresní vztah neobsahuje konstantní člen, je lepší ho nevynechávat kvůli interpretaci koeficientu determinace. Odhadnutý absolutní člen (označovaný (Intercept)) by se ale v tomto případě neměl významně lišit od nuly. 20 Stochastické modely časových řad (11. cvičení) > n <- length(res.standard) > x <- res. standard[1: (n - 1)] > y <- res. standard[2:n] > m. res <- Im (y ~ x) > summary (m. res) Call: Im (formula = y ~ x) Residuals: Min 1Q Median 3Q Max -2.4131 -0.5804 -0.0001 0.6429 3.4238 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.001597 0.070248 0.023 0.981889 x 0.268755 0.070133 3.832 0.000173 *** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.9708 on 189 degrees of freedom Multiple R-squared: 0.0721, Adjusted R-squared: 0.06719 F-statistic: 14.68 on 1 and 189 DF, p-value: 0.0001729 Z výsledků je patrné, že podle očekávání se absolutní člen významně neliší od nuly. Ale směrnice již významná je (na 5% hladině), protože má p-hodnotu 0, 000173, která je menší než 0, 05. To je také důvod, proč Durbin-Watsonův test zamítl hypotézu, že rezidua jsou nekorelovaná. Výsledky znázorníme graficky > par(mfrow = c(l, 1), mar = c(2, 2, 1, 0) + 0.05) > txt <- pasteC'D =", round(DWtest$statistic, 2), " rhol =", roundtt - 0.5 * DWtest$statistic, 2), " p-value =", round(DWtest$p.value, 5)) > plot (x, y) > abline(h = 0, col = "gray") > äblineCm.res, col = "red", lwd = 2) > mtext(txt) Stochastické modely časových řad (11. cvičení) 21 D = 1.46 rho1=0.27 p-value = 0.00013 o ° ° o o o o _ -1-1-1-1-1- r2 -1 0 1 2 , 3 Obrázek 10: Grafické testování autokorelace u standardizovaných reziduí metody malého trendu pro data Koncentrace CO2 u sopky Mauna Loa v letech 1965-1980 6 Úkol: Načtěte si datový soubor zoo_jihlava.txt (měsíční údaje o počtu návštěvníků v ZOO Jihlava v letech 2006-2010), na načtená data aplikujte modifikovanou metodu malého trendu a získaná a rezidua analyzujte z hlediska • odlehlých a vybočujících pozorování • normality • heteroskedascity • nezávislosti • autokorelace reziduí