Stochastické modely časových řad (5. cvičení) 1 M5201 - 5. CVIČENÍ: Transformace časových řad stabilizující rozptyl 1 Transformace stabilizující rozptyl Nechť náhodná veličina X má rozdělení, které závisí na nějakém parametru 9. Předpokládejme, že tento parametr je zvolen tak, aby platilo EqX = 9. Ve většině případů (ne však u normálního rozdělení) na 9 závisí i rozptyl veličiny X, takže můžeme psát DeX = a2{9). Přitom cr{0) bývá obvykle hladká funkce proměnné 9. Vzniká otázka, zda lze najít netriviální funkci [~g~\ tak, aby náhodná veličina Y = g(X) měla rozptyl nezávisející na 9. (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 funkce g v bodě 9 dostaneme aproximaci g(X)^g(9) + g'(9)(X-9). Na základě rozvoje lze střední hodnotu aproximovat takto Egg{X) « E [g{9) + g'(9)(X - 9)] = g(9) a rozptyl De [g(X)] « [g'(9)]2DeX = [g'{9)}2a2{9). Chceme, aby po transformaci byl rozptyl konstantní a nezávisel na střední hodnotě, tj. „2 n r \Jia\~}2 2/- c c2 = De[g{Yt)} = [g'{9)Ya2{9) g'{9) kde c je nějaká konstanta. Odtud snadno dostaneme tvar transformace stabilizující rozptyl 9(0) = c J -Ljd6 + K. Konstanty c a K se volí tak, aby funkce fg] vypočtená podle předchozího vzorce měla výhodný tvar. Ukázalo se, že takto vypočtená funkce \g nejen výrazně stabilizuje rozptyl, takže rozptyl Dgg(X) závisí na 9 jen velmi málo, ale zároveň také rozdělení náhodné veličiny Y = g(X) bývá již velmi blízké normálnímu, i když třeba samotné rozdělení veličiny X je výrazně nenormální. 2 Stochastické modely časových řad (5. cvičení) 2 Příklad transformace stabilizující rozptyl — Poissonovo rozdělení Nechť náhodná veličina X má Poissonovo rozdělení s parametrem A > 0, tj. X ~ -Po(A) s pravděpodobnostní funkcí px{x) = P(X = x) = '-Lre~X pro x = 0,1, 2,.... x\ a2(X) = A . Pak Lze spočítat, že EX = DX = A, tj. g(\) = c í -4rrdX + K = c [ -^=d\ + K = 2cV\ + K. J o-(A) J VA Obvykle se volí c = h, K = 0 & pracuje se s velmi známou odmocninovou transformací Y = g(X) = VX Ověřme střední hodnotu a rozptyl náhodné veličiny Y: EY = Eg{X) « g{\) = VX DY = Dg{X)^[g\\)]2a\\) 1 1 2 Ta, 3 Mocninné transformace Mějme kladnou náhodnou veličinu X z rozdělení, které závisí na parametru 0 se střední hodnotou a rozptylem A* 2, , , m ^K- tj. I-^a2^). Podle obecného vzorce se transformace stabilizující rozptyl vypočítá takto: cr(/i) a J /i Dále budeme pracovat s parametrem cd^ , K _ c f j §ln\fi\ + K ů = l, A = 1 -ů a budeme ho nazývat transformační parametr pro mocninnou transformaci. 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ě a odtud _ c = o a K Stochastické modely časových řad (5. cvičení) 3 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: N , sn. íln(X + a) X = 0(ů = l), A^O(^l). Mocninná transformace se znaménkem lze opět použít v případě, že náhodné veličiny nejsou kladné: (Y, . ,V,|V|(A) ,sign(X)ln\X\ A = 0 (ů = 1), 3.1 Odhad transformačního parametru mocninné transformace 3.1.1 Parametrický přístup pomocí metody maximální věrohodnosti Mějme nezávislé realizace náhodné veličiny X ~£(ft(7VM). Předpokládejme, že existuje takové A = 1 -ů, že transformovaný náhodný vektor Y = (Y1=g(X1),...,Yn=g(Xn)y je výběr z normálního rozdělení se střední hodnotou fi a rozptylem a2. Označme y = (yi, • • • ,yn)' realizaci náhodného výběru. Hledejme maximum věrohodnostní funkce pro 6 = (fi,a2)', tj. pro funkci i ex 11 / y i - v (27ra2)-f exp í=i y i - v což je stejná úloha jako hledat maximum logaritmu věrohodnostní funkce .In^-^hV)-^ l(fi,a í=i Maxima nalezneme, položíme-li l- = 0a = 0. Lze dokázat, že funkce l(fi, a2) nabývá v bodě (jl, a2) = (y, s2) svého maxima. Celkově jsme tedy dostali, že maxl(fi,az) = l(y,s2) -ln(27r)--ln(S2 n 4 Stochastické modely časových řad (5. cvičení) max L(fi, a2) = L(y, s2) = (2irs n2\ 2e-f_ Nyní toto maximum vyjádřeme v původních proměnných xí, kdy [ ln Xi A = 0, Nejprve vypočtěme jakobián této transformace: dy. Pak ' H- II n Ax^1 n n „A-l i=l i=l max L(/x, u (2tts2(A)) 2 e -§ + (A-l) ^lnii max /(/x, a2, A) ■§ ln(27r) - f ln(s2(A)) - § + (A - 1) E lnxť. í=i Nyní hledejme maximum funkce a2, A) = Z(y, s2, A) pro parametr A. Protože maximum vzhledem k A nezávisí na konstantách, budeme maximalizovat funkci r(A) = -^ln(S2(A)) + (A-l)^lnxí. í=i Teoretickým odvozením maximálně věrohodného odhadu parametru A, se zde dále nebudeme zabývat, ale vysvětlíme si jednodušší přístup: pro různé hodnoty A £ (Ai,A2) (Ai,A2 G M, Ai < A2) se vykreslí do grafu hodnoty Z*(A) a hledá se maximum A v daném intervalu. Pro tento případ odvodili Box a Cox (1964) asymptotické rozdělení statistiky K Na jejím základě lze odvodit interval spolehlivosti pro parametr A: ( l-a = P (K < xUW)=P (-2 - Z*(A)1 < xUW) =P i*W-^(i) Da leží v intervalu spolehlivosti a jsou tedy přijatelná. Testování hypotéz typu Hq : A = Ao proti alternativě H\ : A > Ao: . Pokud hypotézu nezamítneme, tj. Hl : A = 1 1. Nejprve budeme testovat hypotézu > Da, nemusíme data transformovat. 2. Pokud předchozí hypotézu zamítneme, můžeme testovat další hypotézu Hq : x = 0 . Pokud tuto hypotézu nezamítneme, tj. Z*(0) > Da A < Da, transformace bude tvaru Ví = Inxí. Stochastické modely časových řad (5. cvičení) 5 Pokud však se Z*(0) < Da A /*(!) < Da, provedeme transformaci Vi = A 3.1.2 Jednoduchý algoritmus v praktických úlohách — funkce powtr() Další možností, jak odhadnout transformační parametr A je následující algoritmus: 1. Algoritmus nejprve zkontroluje vstupní data tak, aby byla nezáporná, tj. případně přičte kladnou konstantu. Upravený vektor dat rozdělí (podle nějakého dalšího kritéria, pokud nejsou opakovaná pozorování; např. u časových řad jsou data uspořádána podle časového kritéria) na krátké úseky o délce 4 až 12 údajů. V každém úseku dat se provede pokud možno robustní odhad polohy p, (průměr, medián) a robustní odhad variability fileDat<-paste(data.library, "UKgas. txt", sep=" ") > con<-file(fileDat) > (popis<-readLines(con,n=l)) [1] "UK Quarterly Gas Consumption" > close(con) > UKgas <- scan(íileDat,skip=l) Z vektoru UKgas vytvoříme objekt typu časové řady a data vykreslíme. > UKgasTS<-ts(UKgas,start=1963,frequency=4) > par(mar=c(2,2,l,0)+0.5) > plot(UKgasTS) O O CM O O O O O oo o o CD O O o o CM MM I I -1-1-1-1-1-r 1965 1970 1975 1980 1985 1990 Obrázek 1: Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje Z grafu je na první pohled patrné, že s rostoucí střední hodnotou dat roste také jejich variabilita, a tudíž bude nutné data transformovat. K tomu použijeme mocninnou transformaci. Nejprve si můžeme udělat představu o transformačním parametru A, který bychom získali při odhadu metodou maximální věrohodnosti. Obrázek s věrohodnostní funkcí, vyznačeným odhadem A metodou maximální věrohodnosti a 95% intervalem spolehlivosti získáme pomocí funkce boxCoxQ z balíčku car. Stochastické modely časových řad (5. cvičení) 7 > libraxy (car) > boxCox(UKgELsTS ~1) Prostřední svislá přerušovaná čára v obrázku vyznačuje maximálně věrohodný odhad A, který odpovídá maximu věrohodnostní funkce. Dvě krajní přerušované čáry vymezují 95% interval spolehlivosti pro odhad A, který má přibližně meze —0.5 a 0. Vidíme, že interval ještě těsně obsahuje nulu, tudíž není vyloučena možnost použití logaritmické transformace dat. Dále při odhadu parametru A = 1 — ů vyzkoušíme jednoduchý přístup založený na regresním modelu (v R je implementován ve funkci powtrO). Funkce powtrO rozdělí nezávisle proměnnou, tj. čas, na subintervaly obsahující takový počet pozorování, jaký určíme parametrem seglen (doporučuje se volit číslo mezi 4 až 12). V každém subintervalu (tj. na základě 4 až 12 pozorování v závislosti na hodnotě parametru seglen) se provede odhad polohy a variability (tj. odhad /x a a). Pomocí parametrů location a variability specifikujeme, jaký typ statistiky se při tom má použít. Možné volby jsou: location="mean" odhad polohy pomocí průměru location="median" odhad polohy pomocí mediánu variability="sd" odhad variability pomocí směrodatné odchylky variability="iqr" odhad variability pomocí interkvartilového rozpětí variability="range" odhad variability pomocí rozdílu mezi maximem a minimem 8 Stochastické modely časových řad (5. cvičení) Funkce powtr() využívá vztah = 07/ ln( fileSkript<-paste(data.library,"FunkceM5201.R",sep="") > source(íileSkript) > outp<-powtr (UKgas, seglen=8,íigure=TRUE,location="median", variabilitj="iqr ") > str(outp) List of 3 $ lambda : Named mm -0.435 ..- attr(*, "names")= chr "Estimate" $ transfx: imm [1:108] 0.11 0.121 0.145 0.125 0.11 ... $ txt : chr "POWER GROWTH OF VARIANCE (power transform): transx=x.~-0.434584" Stochastické modely časových řad (5. cvičení) 9 POWER GROWTH OF VARIANCE (power transform): transx=x.A-0.434584 o CD LO CO ° , ''b = 1.434584. lambda = -0.434584, Cl = (-0.573656 , -0.295512) 5.0 5.5 6.0 6.5 logLocation Obrázek 2: Mocninná transformace dat pomocí funkce powtr - regresní přímka pro logaritmy polohy a variability (volba f igure=TRUE) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje " Na obrázku vidíme body, jejichž x-ové souřadnice jsou dány jako logaritmy odhadů polohy a y-ové souřadnice jako logaritmy odhadů variability pro úseky dat o délce 8 pozorování. Body je proložena regresní přímka. Výsledky regrese jsou uvedeny v obrázku dole - odhad směrnice regresní přímky b, odhad transformačního parametru A = 1 — b a 95%ní interval spolehlivosti pro A. Protože interval spolehlivosti pro A neobsahuje ani jedničku, ani nulu, byla provedena mocninná transformace s odhadnutým parametrem A = —0.434584. Funkce powtr nabízí další zajímavý graf, který názorně ukazuje, jak vypadá variabilita dat v jednotlivých segmentech před a po transformaci (stačí zadat f igure2=TRUE). Tento graf by nám měl také ukázat, zda je vůbec mocninná transformace vhodná. Pokud by nedošlo ke stabilizaci rozptylu ani po transformaci, pak by bylo třeba hledat jiný typ transformace než je mocninná. > outp<-powtr(UKgas, seglen=8,íigure2=TRUE,location="median", variability="iqr") 10 Stochastické modely časových řad (5. cvičení) c/í CD O O CM O O O O O 00 CO > co o E O) O O o o CM o c/í CD > o T3 CD O E i_ o M— C/3 H—' O CD O seglen = 8, location = median, variability = iqr o o ■H o o ° w ■ o o 4 ~~r 6 8 10 12 POWER GROWTH OF VARIANCE (power transform): transx=x.A-0.434584 r~ 14 lambda = -0.434584, Cl = (-0.573656 , -0.295512)_ -1-1-1- -1-1-1- 2 4 6 8 10 12 14 Obrázek 3: Mocninná transformace pomocí funkce powtr - odhady polohy a variability v jednotlivých segmentech (volba f igure2=TRUE) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje " Stochastické modely časových řad (5. cvičení) 11 Obrázek naznačuje, že variabilita transformovaných dat je nezávislá na střední hodnotě v jednotlivých segmentech a je víceméně konstantní. Obdobný výstup, pouze vyjádřený pomocí krabicových grafů za jednotlivé segmenty vstupních dat, získáme volbou f igure3=TRUE. > outp<-powtr(UKgas, seglen=8,íigure3=TRUE,location="median",variability="iqr") o o o o o o o oo o o CD o o o o seglen = 8, location = median, variability = iqr POWER GROWTH OF VARIANCE (power transform): transx=x.A-0.434584 1 2 3 4 5 6 7 8 9 10 11 12 13 Obrázek 4: Mocninná transformace pomocí funkce powtr - krabicové grafy v jednotlivých segmentech (volba f igure3=TRUE) pro data „Spotřeba zemního plynu ve Velké Británii -čtvrtletní údaje " 12 Stochastické modely časových řad (5. cvičení) Abychom vyzkoušeli robustnost funkce powtr() na našich datech (neboli zda je odhad A citlivý na to, jakou hodnotu seglen zvolíme), provedeme odhad mocninné transformace postupně pro seglen=4,6,10,12 . Ostatní parametry při tom necháme beze změny K znázornění výsledku vybereme třetí typ grafu. > outp<-powtr(UKgas, seglen=4, í igure3=ŤRUE, location="median ", variability="iqr") o o o o o o o oo o o CD o o o o seglen = 4, location = median, variability = iqr i—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I 1 3 5 7 9 11 13 15 17 19 21 23 25 27 co CÖ CM Ö O Ö 00 O CD O o POWER GROWTH OF VARIANCE (power transform): transx=x.A-0.470665 i i i i i _i i i i lambda = -0.470665, Cl = (-0.584881 , -0.356448) ~i—i—i—i—i—i—i—i—i—i—i—i—i—rn—i—i—i—i—i—i—rn—i—i—i—r~ 1 3 5 7 9 11 13 15 17 19 21 23 25 27 Obrázek 5: powtr - krabicové grafy (volba seglen=4 a f igure3=TRUE) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Totéž znovu zopakujeme pro volbu seglen=6 Stochastické modely časových řad (5. cvičení) 13 > outp<-powtr(UKgas, seglen=6,íigure3=TRUE,location="median", variability="iqr") o o o o o o o oo o o CD o o o o seglen = 6, location = median, variability = iqr n—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 POWER GROWTH OF VARIANCE (power transform): transx=x.A-0.426178 ö ö o ö 00 o CD O 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Obrázek 6: powtr - krabicové grafy (volba seglen=6 a f igure3=TRUE) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Další volba je seglen=10 . > outp<-powtr(UKgas, seglen=10, íigure3=TRUE,location="median", variabilitj= "iqr ") 14 Stochastické modely časových řad (5. cvičení) o o o O O O O 00 O o CD O O o o seglen = 10, location = median, variability = iqr o 00 o CD O ^-O outp<-powtr(UKgas, seglen=12, íigure3=TRUE,location="median", variabilitj= "iqr ") Stochastické modely časových řad (5. cvičení) 15 seglen = 12, location = median, variability = iqr o HistFit(UKgELs) > mtext("originál values",side=3,line=-0.5,cex=0.95) Histogram, Kernel Density Estimate, Normal Curve original values - Histogram ---- Kernel Density Estimate - Normal Density 0 200 400 600 800 1000 1200 x Obrázek 9: Testování normality pro netrasformovaná data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje " Nejprve výsledky pro mocninnou transformaci s A = —0.5. > boxplotSegments (UKgas ~ (-0. 5) , seglen=8) Stochastické modely časových řad (5. cvičení) 17 X n-1-1-1-1-1-1-1-1-1-1-1-1 1 2 3 4 5 6 7 8 9 10 11 12 13 segments Obrázek 10: Porovnání variability pro trasformovaná data (Y = Pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" > HistFit (UKgas~ (-0. 5) ,xlab=expression(y == x~-0.5)) Histogram, Kernel Density Estimate, Normal Curve Histogram Kernel Density Estimate Normal Density 0.02 0.04 0.06 0.08 0.10 y = x Obrázek 11: Testování normality pro trasformovaná data (Y = pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Obdobně dostaneme výsledky pro transformaci Y = X > boxplotSegments (UKgas ~ (-1/3) , seglen=8) 18 Stochastické modely časových řad (5. cvičení) n-1-1-1-1-1-1-1-1-1-1-1-1 1 2 3 4 5 6 7 8 9 10 11 12 13 segments Obrázek 12: Porovnání variability pro trasformovaná data (Y = -37^) Pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" > HistFit(UKgas~(-l/3),xlab=expression(y == x~-0.33)) Histogram, Kernel Density Estimate, Normal Curve - Histogram --- Kernel Density Estimate — Normal Density ✓ / / N W \ \ \ \ \ \ V \ \ / / / y > / / \ \ t / / j / / \ \ \ \ \ \ y \ 1 1 0.10 0.15 0.20 Obrázek 13: Testování normality pro trasformovaná data (Y = -37^) „Spotřeba plynu ve Velké Británii - čtvrtletní údaje " A konečně výsledky pro transformaci Y = X-1/4. > boxplotSegments(UKgas~(-0.25),seglen=8) Stochastické modely časových řad (5. cvičení) 19 ~~i-1-1-1-1-1-1-1-1-1-1-1-r 1 2 3 4 5 6 7 8 9 10 11 12 13 segments Obrázek 14: Porovnání variability pro trasformovaná data (Y = -37^) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" > HistFit (UKgcLS~ (—0.25) ,xla.b=expression(y == x~-0.25)) Histogram, Kernel Density Estimate, Normal Curve - Histogram ---- Kernel Density Estimate - Normal Density Obrázek 15: Testování normality pro trasformovaná data (Y = -37^) pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Z grafů je ihned vidět, že ačkoli je rozdělení transformovaných dat méně sešikmené, stále má k normalitě daleko. Je třeba si však uvědomit, že to nemusí být ani tak volbou transformace, jako spíše faktem, že časová řada má výrazný deterministický trend, který pak přehluší stochastické vlastnosti kolísání kolem trendu. Ihned nás napadne myšlenka nejprve odstranit trend a teprve pro rezidua hledat vhodnou mocninnou transformaci. Tento postup nyní vyzkoušíme. Pokud si pozorně prohlédneme graf dat, napadne nás, že trend by mohl být kvadratický. Zcela korektní postup by vyžadoval, abychom stupeň prokládaného polynomu určili exaktnějším způsobem než je pohled „od oka" (viz např. 3. cvičení), ale tím se zde nyní nebudeme zabývat a rovnou odhadneme model s kvadratickým trendem. 20 Stochastické modely časových řad (5. cvičení) > y <- as.vector(UKgasTS) > x <- as.vector(time(UKgasTS)) > n<-length(y) > nn<-300 > data<-data.frame(x,y) > LinTrend <- lm(y ' x+I(x~2),data=data) > print (summary (LinTrend) ) Call: lm(formula = y " x + I(x~2) , data = data) Residuals: Min 1Q Median 3Q Max -412.12 -65.52 2.89 58.20 457.79 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.652e+06 1.149e+06 2.309 0.0229 * x -2.707e+03 1.162e+03 -2.329 0.0218 * I(x~2) 6.910e-01 2.941e-01 2.349 0.0207 * Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 166 on 105 degrees of freedom Multiple R-squared: 0.5717, Adjusted R-squared: 0.5636 F-statistic: 70.09 on 2 and 105 DF, p-value: < 2.2e-16 > new <- data.frame(x = seq(x[l],x[n],length.out=nn)) > pred.w.plim <- predict (LinTrend, new, interval="prediction") > pred.w.clim <- predict (LinTrend, new, interval="confidence") > par(mar=c(2,2,l,0)+0.5) > matplot (new$x, cbind(pred. w. dim, pred.w.plim[, -1]) , col=c (2,3,3,4,4) , lty=c(l,2,2,3,3), type="l", ylab="predicted y") > lines (x,y) -1-1-1-1-1-r1 1965 1970 1975 1980 1985 1990 Obrázek 16: Lineární trend pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje " Dříve než použijeme mocninnou transformaci, podívejme se pomocí funkce boxplotSegments (), zda to pro rezidua má vůbec smysl. Stochastické modely časových řad (5. cvičení) 21 > x<-resid(LinTrend) > boxplotSegments(x, seglen=8) segments Obrázek 17: boxplotSegments (seglen=8) pro rezidua po lineárním trendu u dat „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Vyzkoušíme mocninnou transformaci (i když to nejspíše nebude mít smysl) a pomocí funkce boxplotSegments() si prohlédneme variabilitu pro transformovaná data. Ještě předtím, než provedeme transformaci, přičteme k reziduím absolutní člen z odhadnuté regresní přímky, abychom prováděli transformaci pro nezáporná data. > x<-resid(LinTrend)+coef(LinTrend) [1] > outp<-powtr(x, seglen=8, íigure=TRUE, location="median ", variability="iqr") CONSTANT VARIANCE (data not transformed): transfx = x CD O seglen = 8, location = median, o variability = iqr LO LO O LÓ o o ——— LO o o O o o o LO b = 20395.496984, lambda = -20394.496984, Cl = (-971066.808085 , 30277.814116) 14.79092 14.79093 14.79093 14.79094 14.79094 14.7909 logLocation Obrázek 18: powtr - regresní přímka pro logaritmy polohy a variability (volba f igure=TRUE) pro rezidua po lineárním trendu u dat „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje " Protože interval spolehlivosti obsahuje jedničku, nemělo by být třeba provádět transformaci. Všimněme si ale dole v obrázku, jak velkých hodnot nabývá odhad A a meze intervalu 22 Stochastické modely časových řad (5. cvičení) spolehlivosti. Nyní se podívejme na normalitu pomocí příkazu HistFit (). > HistFit (outp$transíx,xlab=x) Histogram, Kernel Density Estimate, Normal Curve ,"\ - Histogram ' 1 ---- Kernel Density Estimate I - Normal Density ....._I_J_I_U_........,................................■............,_LJ_I_U_,_J_O_I_L 2651800 2652000 2652200 2652400 2652600 2652800 2652350.1676275 , , , 2652318.39842329 Obrázek 19: lestovam normality pro netrasíormovana data po odstraněni lineárního trendu (Y = X) „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Ačkoli by se třeba pouhým pohledem mohlo zdát, že rezidua mají k normálnímu rozdělení trochu blíže než data po předchozích transformacích, je třeba si uvědomit, že hlavním cílem bylo stabilizovat rozptyl, což se nám, jak je patrné z obrázku s krabicovými grafy pro rezidua, nepovedlo. Tento postup se tedy rozhodně neosvědčil. Byl totiž od samotného počátku nekorektní. Klasický regresní model předpokládá homoskedastická rezidua, což evidentně nebylo splněno. Naštěstí existují postupy, které v jednom kroku hledají v regresním modelu všechny neznámé parametry. V prostředí R balíček car nabízí funkci powerTransf orm(), která hledá parametr A pro Box-Coxovu transformaci. > library(car) > TT<-time(UKgasTS) > data<-data.fráme(TIME=TT-mean(TT),X=UKgas) > transfl <-powerTransform (X ~ TIME+I(TIME~2) , data=dataj > summary(transfl) bcPower Transformation to Normality Est.Power Std.Err. Wald Lower Bound Wald Upper Bound Yl -0.362 0.1278 -0.6126 -0.1115 Likelihood ratio tests about transformation parameters LRT df pval LR test, lambda = (0) 7.947835 1 0.004814494 LR test, lambda = (1) 102.969916 1 0.000000000 > str(transfl) Stochastické modely časových řad (5. cvičení) 23 List of 13 $ value : imm 499 $ counts : Named int [1:2] 3 ..- attr(*, "names")= chr [1:2] $ convergence: int 0 $ message : chr "CONVERGENCE: $ hessian : num [1, 1] 61.2 $ start : num -0.362 $ lambda : Named num -0.362 ..- attrO, "names")= chr "Yl" $ roundlam : Named num -0.5 ..- attrO, "names")= chr "Yl" $ family : chr "bcPower" $ xqr :List of 4 REL_REDUCTI0N_0F_F <= FACTR*EPSMCH 3 "function gradient" ii ..$ qr : num [1:108, 1:3] -10.3923 0.0962 0.0962 0.0962 0.0962 ... .. ..- attrO, "assign")= int [1:3] 0 12 .. ..- attr(*, "dimnames")=List of 2 ......$ : chr [1:108] "1" "2" "3" "4" ... ......$ : chr [1:3] "(Intercept)" "TIME" "I(TIME~2)" ..$ rank : int 3 ..$ qraux: num [1:3] 1.1 1.15 1.15 ..$ pivot: int [1:3] 12 3 ..- attr(*, "class")= chr "qr" $ y : num [1:108, 1] 160.1 129.7 84.8 120.1 160.1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:108] "1" "2" "3" "4" ... .. ..$ : NULL $ x : num [1:108, 1:3] 1111111111... ..- attr(*, "dimnames")=List of 2 ....$: chr [1:108] "1" "2" "3" "4" ... .. ..$ : chr [1:3] "(Intercept)" "TIME" "I(TTME~2)" ..- attr(*, "assign")= int [1:3] 0 12 $ weights : num [1:108] 1111111111... - attr(*, "class")= chr "powerTransform" Všimněme si, že transf $roundlam nabízí vhodnou volbu parametru A. Nejedná se přímo o odhad A, ale o takovou hodnotu, která leží v intervalu spolehlivosti a má rozumnou interpretaci. > print (transfl$roundlam) Nyní ukážeme ještě trochu jiný, ale zcela ekvivalentní postup. > ml <- lm(X ~ TIME+I(TIME~2), data=data.) > summary(ml) Call: lm(formula = X ~ TIME + I(TIME~2), data = data) Yl -0.5 Residuals: 24 Stochastické modely časových řad (5. cvičení) Min IQ Median 3Q Max -412.12 -65.52 2.89 58.20 457.79 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 295.6571 23.9670 12.336 <2e-16 *** TIME 23.7878 2.0499 11.604 <2e-16 *** I(TIME~2) 0.6910 0.2941 2.349 0.0207 * Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 166 on 105 degrees of freedom Multiple R-squared: 0.5717, Adjusted R-squared: 0.5636 F-statistic: 70.09 on 2 and 105 DF, p-value: < 2.2e-16 > outT <- powerTrELnsform(ml) > summary(outT) bePower Transformation to Normality Est.Power Std.Err. Wald Lower Bound Wald Upper Bound Yl -0.362 0.1278 -0.6126 -0.1115 Likelihood ratio tests about transformation parameters LRT df pval LR test, lambda = (0) 7.947835 1 0.004814494 LR test, lambda = (1) 102.969916 1 0.000000000 > print(outT$roundlam) Yl -0.5 > m2<-update(ml, basicPower(outT$y,outT$roundlam) ~.) > summary (m2) Call: lm(formula = basicPower(outT$y, outT$roundlam) ~ TIME + I(TIME~2), data = data) Residuals: Min 1Q Median 3Q Max -0.018379 -0.009724 -0.003079 0.010398 0.023211 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.341e-02 1.736e-03 36.53 <2e-16 *** TIME -2.195e-03 1.485e-04 -14.78 <2e-16 *** I(TIME~2) 2.599e-05 2.130e-05 1.22 0.225 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.01203 on 105 degrees of freedom Multiple R-squared: 0.677, Adjusted R-squared: 0.6708 F-statistic: 110 on 2 and 105 DF, p-value: < 2.2e-16 Stochastické modely časových řad (5. cvičení) 25 Na závěr se ještě podívejme, jak dopadla rezidua u transformovaného modelu. > res<-resid(m2) > boxplotSegments(res,seglen=8,xlab="Residuals oí Box-Cox model") | ~l I I I I I I I I I I I I 1 2 3 4 5 6 7 8 9 10 11 12 13 Residuals of Box-Cox model Obrázek 20: boxplotSegments (volba seglen=8) rezidua v Box-Coxově modelu pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" > HistFit(res,xlab="Residuals oí Box-Cox model") Histogram, Kernel Density Estimate, Normal Curve - Histogram ---- Kernel Density Estimate - Normal Density -0.02 -0.01 0.00 0.01 0.02 Residuals of Box-Cox model Obrázek 21: Testování normality reziduí v Box-Coxově modelu pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Na závěr můžeme říci, že se nám podařilo najít transformaci, po níž mají data konstantní rozptyl. Normalita se ale ani v tomto případě příliš nezlepšila. Protože z výpisu informací o modelu m2 je vidět, že koeficient u druhé mocniny času se neliší od nuly na hladině významnosti 5 %, vyzkoušíme ještě jako poslední možnost odstranit z dat pouze lineární trend a zároveň najít vhodný transformační parametr. 26 Stochastické modely časových řad (5. cvičení) > m3 <- lm(X ~ TIME, data=dataj > summary (m3) Call: lm(formula = X ~ TIME, data = data) Residuals: Min 1Q Median 3Q Max -368.57 -85.50 1.20 69.86 525.95 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 337.631 16.314 20.70 <2e-16 *** TIME 23.788 2.093 11.36 <2e-16 *** Signif. codes: 0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 169.5 on 106 degrees of freedom Multiple R-squared: 0.5492, Adjusted R-squared: 0.545 F-statistic: 129.2 on 1 and 106 DF, p-value: < 2.2e-16 > outT <- powerTransform(m3) > summary (outT) bePower Transformation to Normality Est.Power Std.Err. Wald Lower Bound Wald Upper Bound Yl -0.3275 0.121 -0.5647 -0.0903 Likelihood ratio tests about transformation parameters LRT df pval LR test, lambda = (0) 7.299312 1 0.006898103 LR test, lambda = (1) 107.836554 1 0.000000000 > print(outT$roundlam) Yl -0.5 > m4<-update(m3,basicPower(outT$y,outT$roundlam)~.) > summary (m4) Call: lm(formula = basicPower(outT$y, outT$roundlam) ~ TIME, data = data) Residuals: Min 1Q Median 3Q Max -0.018188 -0.010034 -0.002811 0.009518 0.024849 Coefficients: Stochastické modely časových řad (5. cvičení) 27 Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0649838 0.0011598 56.03 <2e-16 *** TIME -0.0021950 0.0001488 -14.75 <2e-16 *** Signif. codes: 0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.01205 on 106 degrees of freedom Multiple R-squared: 0.6724, Adjusted R-squared: 0.6693 F-statistic: 217.6 on 1 and 106 DF, p-value: < 2.2e-16 I v případě jednoduššího lineárního trendu je navržena stejná transformace jako v případě kvadratického trendu. Na závěr si prohlédněme transformovaná rezidua pro lineární trend. > res<-resid(m4) > boxplotSegments(res,seglen=8,xlab="Residuals oí Box-Cox model") 10 "T" 11 12 13 Residuals of Box-Cox model Obrázek 22: boxplotSegments (volba seglen=8) rezidua v Box-Coxově modelu pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" > HistFit(res,xlab="Residuals oí Box-Cox model") 28 Stochastické modely časových řad (5. cvičení) Histogram, Kernel Density Estimate, Normal Curve Histogram Kernel Density Estimate Normal Density -0.02 0.02 -0.01 0.00 0.01 Residuals of Box-Cox model Obrázek 23: Testování normality reziduí v Box-Coxově modelu pro data „Spotřeba zemního plynu ve Velké Británii - čtvrtletní údaje" Oproti modelu s kvadratickým trendem není na první pohled patrný žádný výrazný rozdíl. 4 Úkol: Načtěte data uložená v souboru j j .dat. • Data vykreslete do grafu a posuďte, zda bude třeba je transformovat. • Použijte funkci boxCoxO ke zjištění odhadu parametru mocninné transformace A metodou maximální věrohodnosti. Obsahuje interval spolehlivosti některou z hodnot 1 nebo 0, což by znamenalo, že data není třeba transformovat, resp. lze použít logaritmickou transformaci? • Dále použijte funkci powtr s různými volbami parametru seglen a na základě výsledků zvolte vhodný transformační parametr A. • Na závěr pomocí funkce powerTransf orm zároveň odhadněte lineární trend i transformační parametr A. • U všech možných transformací vykreslete histogram a porovnávejte, zda se rozdělení dat blíží normálnímu.