Stochastické modely časových řad (4. cvičení) 1 M5201 - 4. CVIČENÍ: Identifikace periodických komponent 1 Regresní modely cyklických trendů 1.1 Nelineární model Uvažujme následující model s nekorelovanými a homoskedastickými rezidui (tj. majícími konstantní rozptyl): Yí=f(tí)+Eí Ea=0; Da = a ; C(eí, Ej) = 0; i ^ j; i, j = 1,..., n Je-li f (t) periodická funkce s periodou T, pak frekvencí rozumíme veličinu \ _ 2tt Dále uvažujme takovou f(ti), kterou lze zapsat ve tvaru nebo ekvivalentně (a) f (ti) = p + Yľj=i(aj cos ujjti + Pj sinujjti (b) f (U) = /x + Ej=i P j cos(lújU + 6j) pj = yaj + j3j; Oj = arctan Jde o nelineární regresní model vzhledem k (3p + 1) neznámých parametrů: (a) ai,...,ap Pi,...,f3p p ui,...,up (b) pi,...,pp p ui,...,up 6i,...,6p Odhad vektoru neznámých parametrů pomocí metody nejmenších čtverců minimalizuje výraz (a) S(fi,a1,... ,ap,/3i,... ,/3p,tJi, ...,up) = Ya=i (Yí ~ f(tí)Ý (b) 5(p,pi,... ,pp,6i, ■ ■ .,6p,u}i, ...,u>p) = Eľ=i (Y - f (U))2 Numericky lze systém nelineárních rovnic řešit např. pomocí Gaussovy-Newtonovy iterační metody. V prostředí R lze v základním statistickém balíčku stats najít funkci nls(), kterou je možno použít k odhadu parametrů nelineárních regresních modelů jako je tento. 1.2 Lineární model pro známé frekvence Pokud jsou frekvence ui,... ,up známé, situace je podstatně jednodušší. Z modelu (a) se v tom případě stane lineární regresní model, jehož koeficienty můžeme odhadnout metodou nejmenších čtverců. Vektor neznámých parametrů lze psát takto /3 = (p, ct\,..., ap, [3\,..., [3P)' a matice plánu je tvaru (1 coscjiíi sincjiíi ••• costjpíi sintJpíiN : : : ' •. : : 1 costJiín sintJiín ••• cos loptn sin tjptn, 2 Stochastické modely časových řad (4. cvičení) Platí X'X /n O O - u 2 \o o ze vztahu j3 = (X'X) 1 X'Y dostaneme 1 v^n v = I Eľ=i yí cos Pro j = 1, Neznámé parametry modelu (b) získáme na základě odhadů pro model (a) ze vztahů 7i pro j = 1,... ,p čí;,- = arctan t^- 1.3 Metoda skrytých period Metoda skrytých period slouží k nalezení významných period lji, ..., ljp, pokud nejsou předem známé. Uvažujme takové časové řady, které lze rozložit na součet harmonických frekvencí, jejichž délky period můžeme vyjádřit jako podíl a 0 < k < n. , kde n je počet naměřených hodnot Označme dále f =- = - Jk Tk n Lúk = 'žit f k = 2lľ k n k-tá frekvence k-tá uhlová frekvence Maximální délka periody je rovna počtu pozorování, tj. Tmax = n , tedy k = 1 a minimální frekvence má velikost túr, 2tt Nejkratší zjistitelná perioda je Nyquistova frekvence). ± mm — ^ jí odpovídající frekvence je (tz může nabývat hodnot: k 1 2 ^ 1 2 pro n sudé pro n liché. Model časové řady pak můžeme zapsat ve tvaru Yt = st + £t kde t označuje ekvidistatní časové okamžiky měření, pro jednoduchost předpokládáme, že intervaly mají jednotkovou velikost; Stochastické modely časových řad (4. cvičení) • n je počet naměřených hodnot; • et je bílý šum: et ~ WN(0,a2), tj. Eet = 0, Det = a2, C(et,es)=Eeteh = 0; t ^ h, • st je periodická funkce tvaru st Ylj=i (aj cos t^j + P j sin t^j) Pro n liché Y^j=i(aj costujj + [3j sin. tu j) + a »( —1)* pro n sudé kde oij,f3j é i, j) é N jsou daná čísla a nazýváme je parametry modelu. Ekvivalentně můžeme st napsat ve tvaru: st Yfj=i P j cos(tu>j — Oj) pro n liché YTj=\ P j cos(tujj — Oj) + a» ( — 1)* pro n sudé kde Pk = \/a2k + Pl je amplituda k-té harmonické složky, I arctan — au > 0 0k = l k y arctan ^ + tv ak < 0 je fázový posun k-té harmonické složky. Pokud zahrneme do počtu frekvencí i frekvenci nulovou (kdy k = 0), pak přibude člen, který označíme přitom [3q = 0. 1.4 Periodogram Ukazatelem významných periodicit náhodné posloupnosti {Yt,t £ {t±,... ,tn}} je periodogram: Za platnosti uvažovaného modelu Yt = st + £t = -^- + cos tuJi + Po sin + £t £t ~ ^-^(0) °"2 bude mít periodogram pro velká n v bodech cj = UJl, výrazně velké hodnoty, jinde jeho hodnoty budou relativně malé, budou kolísat kolem hodnoty 2tt Odtud plyne, že významná lokální maxima v průběhu periodogramu by měla identifikovat periodickou strukturu uvažovaného modelu tak, že vyznačí neznámé frekvence uji, ... ,ujp. Vhodným statistickým testem pak rozhodneme, které hodnoty periodogramu můžeme opravdu považovat za významně velké ve srovnání s hodnotami ostatními. 4 Stochastické modely časových řad (4. cvičení) 1.5 Test R. A. Fishera R. A. Fisher odvodil test, kterým se dá zjistit významnost nejvyšších hodnot periodogramu. Uvažujme posloupnost nezávislých náhodných veličin Y±,..., Yn. Nulová hypotéza: H0: Yt = et st^WN(0,a2). Předpokládejme, že n = 2m + 1 je liché číslo (je-li n sudé, obvykle se vynechá první člen) a uvažujme hodnoty periodog ramu In(Lúfc) v bodech lo^ — — 1,..., m). Srovnejme je podle velikosti a označme Vi,..., Vm. Test se provádí na základě tzv. Fisherovy statistiky W Vi Vi + • • • + v„ která nabývá hodnot mezi nulou a jedničkou. Velké hodnoty (blízké 1) budou tvořit kritický obor nulové hypotézy proti alternativě Hi : Yt = £?=1 Pj cos(í^- - Oj) +st et ~ WN(0, a2) R. A. Fisher odvodil distribuční funkci statistiky W za platnosti hypotézy Hq : (za předpokladu, že uvažujeme gaussovský bílý šum) W\H0\ P{W > x\Hq) = m(l - x iín—1 m , 2 , (1 - 2x)m-1 + kde 0 < x < 1 a na pravé straně sčítáme tak dlouho, dokud jsou členy (1 — kx) kladné, což lze také zapsat takto P(W > x\H0) = £ Q max(0,1 — kx)] m—l E fc=l,2,... m , k , [(1 - kx). im-1 Hypotézu Hq : Yt = et £t ~ WN(0, a2) zamítáme na hladině významnosti a, pokud 1 - Fw\Ho{w) = P(W > w\H0) = aw{H(j < a, kde w je skutečná hodnota Fisherovy statistiky při daných hodnotách konkrétní časové řady a ctw\Ho Je ^zv p-value. 1.5.1 Modifikace Fisherova testu pro testování dalších periodicit V případě, že pomocí Fisherova testu zjistíme signifikantní periodicitu určité frekvence, postupujeme při testování dalších periodicit následujícím způsobem. Vynecháme V\, na základě veličin V2,..., Vm položíme a stanovíme P(W^ > w^) podle stejného vzorce, kde místo m dosadíme m — 1. Vyjde-li i tato druhá největší hodnota významná, opět se vynechá a m se zmenší o další jedničku. Když takto stanovíme frekvence ui,... ,up, získáme model se známými frekvencemi a zbylé neznámé koeficienty odhadneme metodou nejmenších čtverců (viz odstavec Lineární model pro známé frekvence). Stochastické modely časových řad (4. cvičení) 5 1.6 Siegelův test V případě více významných period není Fisherův test dostatečně silný. Vhodnější by měl být Siegelův test. Zde se používá statistika T\ tvaru kde (#)+ = max(x,0), gp je 100(1 — a)% kritická hodnota Fisherova testu (tj. P(W > gp\Ho) = a)aO t\ , kde t\ je kritická hodnota tohoto testu. Kritické hodnoty bývají tabelovány pro různá A a m. Jako významné jsou uznány pouze ty periodicity, jejichž odpovídající sčítance přispěly do celkové hodnoty testovací statistiky Tx._ Pokud počet dat odpovídá násobku délky periody, najde Fisherův (či Siegelův) test skutečnou periodu přesně. V případě, že však počet dat není násobkem délky periody nebo snad je dat méně než je délka periody, uvedené testy skutečnou periodu neobjeví. 1.7 Iterativní metoda (Damsleth a Spj0tvoll, 1982) Následující iterativní metoda představuje alternativní způsob určení významných frekvencí a odhadu regresních koeficientů v modelech s trigonometrickým trendem. Postupujeme takto: (1) Pomocí Fisherova (nebo jiného) testu zjistíme významnou frekvenci ujq. Pokud taková neexistuje, končíme. (2) Minimalizujeme nelineární regresní model Yi = /x + a cos tiu + [3 sin Uu + e% vzhledem k neznámým parametrům (/i, a, 0, uj), přičemž jako počáteční hodnoty bereme • ujq frekvenci nabídnutou v kroku (1) Fisherovým či jiným testem, • fiQ, ao, [3q získáme jako řešení lineárního regresního modelu Yi = fio + a0 cos tiUQ + /30 sin tu q + d Získané hodnoty označme fi, a, 0, u). (3) Z řady odstraníme vliv ú, tj. Zi = Yf — fi — á cos tíu) — [3 sin ijú; + r)i a s řadou Z{ jdeme na bod (1) 6 Stochastické modely časových řad (4. cvičení) Algoritmus při odhadu frekvencí vychází z významných frekvencí, které určil Fisherův či jiný test. Současně s odhadem dalších parametrů odhadnuté frekvence ještě zpřesňuje. Jeho nevýhodou je, že odhaluje frekvence postupně, tj. v každém kroku hledá pouze jedinou významnou frekvenci. Datový soubor, kterým se budeme v příkladu zabývat, se týká vzdáleností letokruhů jednoho stromu, který rostl v Kašmíru v letech 1753-1980. Vzdálenosti byly transformovány tak, že jsou vyjádřeny v bezrozměrných jednotkách. Data jsou uložena v souboru treerings.txt. První řádek obsahuje nadpis, druhý řádek je volný a od třetího řádku dál jsou v jednom sloupci uvedeny vzdálenosti letokruhů. Nejprve načteme nadpis. > íileDat <- paste (data.library, "treerings. txt", sep=" ") > (TXT <- paste(scan(fileDat,what="", nlines=l) ,collapse=" ")) [1] "Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980" Pak načteme vlastní data do datového rámce a vypíšeme jeho strukturu. > rings<-read.table(ÍileDat,header=F,skip=2) > str(rings) 'data, frame': 227 obs. of 1 variable: $ VI: num 0.912 0.734 1.024 0.994 0.98 ... Z načtených dat vytvoříme časovou řadu a vykreslíme ji. > ringsTS<-ts(rings, start=l753,írequency=l) > par(mar=c(2,2,l,0)+0.5) > plot(ringsTS,main=TXT,cex.main=0.85) 2 Příklad Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980 o 1750 1800 1850 1900 1950 Obrázek 1: Šířky letokruhů Stochastické modely časových řad (4. cvičení) 7 Abychom zjistili významné periody, budeme chtít vypočítat periodogram. Než tak budeme moci učinit, bude nutné z časové řady odstranit lineární trend. Proto nejprve odhadneme model s lineárním trendem a rezidua, která na základě tohoto modelu získáme, použijeme k výpočtu periodogramu. > 7 <- as.vector(ringsTS) > x <- as.vector(time(ringsTS) ) > n<-length(y) > nn<-300 > data<-data. frame (x,y) > LinTrend <- lm( y ~ x,data=dataj > print (summary (LinTrend) ) Call: lm(formula = y ~ x, data = data) Residuals: Min 1Q Median 3Q Max -0.35750 -0.13113 0.00535 0.10955 0.65982 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.3470184 0.3420518 3.938 0.00011 *** x -0.0001848 0.0001832 -1.009 0.31415 Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1809 on 225 degrees of freedom Multiple R-squared: 0.004503, Adjusted R-squared: 7.826e-05 F-statistic: 1.018 on 1 and 225 DF, p-value: 0.3142 Odhadnutý lineární trend zakreslíme do grafu a spolu s ním vyneseme také intervaly spolehlivosti kolem trendu a predikční interval. > 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) > plot(x,y,type="l") > matlines(new$x, cbind(pred. w. dim, pred.w.pliml, -11) , col=c(2,3,3,4,4) , lty=c(l,2,2,3,3), type="l", ylab="predicted y") > title (main=TXT, cex.main=0. 85) 8 Stochastické modely časových řad (4. cvičení) Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980 CO CD O CM O 00 O 1750 1800 1850 1900 1950 Obrázek 2: Lineární trend pro šířky letokruhů. Přestože není lineární trend statisticky významný (p — value = 0.31415 > 0.05), odstraníme jej z časové řady a vstupem do funkce, která vypočítává periodogram, budou rezidua. Ta získáme pomocí funkce resid. > res<-resid(LinTrend) K výpočtu periodogramu slouží funkce tperiodO. Tato funkce není součástí R, ale je obsažena v souboru FunkceM5201 .R. Abychom ji mohli použít, musíme zmíněný soubor nejprve načíst. Načtení provedeme příkazem source, kde jako argument zadáme cestu k souboru. > fileFun<-paste(data.library, "FunkceM5201 .R", sep=" ") > source(íileFun) Nyní již můžeme funkci tperiodO použít. Tato funkce na požádání (pokud zadáme argument f ig=TRUE) periodogram také vykreslí do grafu. > pdgV<-tperiod(res,fig=TRUE) > str(pdgV) List of 3 $ I : Named imm [1:113] 0.382 0.356 0.222 0.726 0.657 ... ..- attrO, "names")= chr [1:113] "2" "3" "4" "5" ... $ kF: Named imm [1:8] 12 4 5 11 1 16 2 17 ..- attrO, "names")= chr [1:8] "13" "5" "6" "12" ... $ kS: Named imm [1:9] 12 4 5 11 1 16 2 17 3 ..- attrO, "names")= chr [1:9] "13" "5" "6" "12" ... Stochastické modely časových řad (4. cvičení) 9 PERIODOGRAM (Siegeľs test) CO Ö f 18.9 5618 # 45.4 CD Ö ö ti Ö o o Ö 0 20 40 60 80 100 Index Obrázek 3: Periodogram pro šířky letokruhů Na obrázku jsou vyneseny hodnoty periodogramu pro periody ^, na ose x jsou uspořádány podle k. Ty periody, které jsou významné na základě Fisherova testu, jsou v grafu označeny hvězdičkou a popsány (číslo udává periodu). Výstupem funkce tperiod je (kromě grafu) objekt typu seznam obsahující tři položky První položka I obsahuje hodnoty periodogramu pro všechny periody ^. Ve druhé položce kF jsou uložena ta k, jimž odpovídající periody označil Fisherův test jako významné. A třetí položka kS obsahuje k, jimž přísluší periody, které jsou významné podle Siegelova testu. Vypišme hodnoty významných period, jak nám je nabízí Fisherův a Siegelův test. > (TkF<-n/pdgV$kF) 13 5 6 12 2 17 3 18 18.91667 56.75000 45.40000 20.63636 227.00000 14.18750 113.50000 13.35294 > (TkS<-n/pdgV$kS) 13 5 6 12 2 17 3 18 18.91667 56.75000 45.40000 20.63636 227.00000 14.18750 113.50000 13.35294 4 75.66667 Chceme-li explicitně zjistit a do grafu vykreslit trigonometrický trend, použijeme k tomu dvě funkce, a to trigof it() a trigoval(). Opět se jedná o funkce, které nejsou součástí R. Protože jsou ale uloženy ve stejném souboru jako tperiod(), máme je již načtené, a tudíž je můžeme rovnou používat. První z těchto dvou funkcí - trigof it () - opět hledá pro zadanou časovou řadu významné periody, ale navíc odhadne koeficienty v lineárním regresním modelu s trigonometrickým trendem, v němž byly použity nalezené významné periody Jejím výstupem je seznam s položkami: 10 Stochastické modely časových řad (4. cvičení) • periodogram (hodnoty periodogramu), • freq (odpovídající úhlové frekvence), • signA (odhady koeficientů a j pro významné periody v modelu dle zápisu (a)), • signB (odhady koeficientů [3j pro významné periody v modelu dle zápisu (a)), • signFreq (významné frekvence), • amplituda (odhady amplitud p j pro v modelu ve tvaru (b) se všemi periodami), • faze (odhady fází Oj v modelu (b) se všemi periodami), • signL (logický vektor s hodnotami TRUE na pozicích k odpovídajících významným frekvencím) , • A (odhady koeficientů a j v modelu (a) se všemi periodami) • B (odhady koeficientů [3j v modelu (a) se všemi periodami) > TrigoTrend<-trigofit(res) > str(TrigoTrend) List of 10 $ periodogram $ freq $ signA $ signB $ signFreq $ amplituda $ faze $ signL $ A $ B imm [1:113] 0.0304 0.0284 0.0177 0.0578 0.0523 imm [1:113] 0.0277 0.0554 0.083 0.1107 0.1384 . imm [1:8] -0.0704 -0.0497 0.0241 0.0566 0.0261 imm [1:8] 0.0428 0.0627 0.0722 0.0245 -0.0518 . imm [1:8] 0.3322 0.1107 0.1384 0.3045 0.0277 .. imm [1:113] 0.00545 0.00526 0.00415 0.00751 0.00714 imm [1:113] 2.26 -5.6 6.58 -3.13 -0.15 ... logi [1:113] TRUE TRUE FALSE TRUE TRUE FALSE ... imm [1:113] 0.0261 -0.0327 0.0255 -0.0497 0.0241 ... imm [1:113] -0.0518 0.0455 0.0362 0.0627 0.0722 ... > (FisherPeriods<-2*pi/TrigoTrend$signFreq) [1] 18.91667 56.75000 45.40000 20.63636 227.00000 14.18750 113.50000 [8] 13.35294 Významné frekvence, které vypočítá funkce trigofitO, jsou úhlové frekvence, zatímco funkce tperiodO dává jako výsledek indexy k odpovídající významným periodám. Abychom mohli výsledky navzájem srovnat, musíme je přepočítat. Indexy k odpovídající významným periodám: > pdgV$kF # v prvním radku jsou jména složek vektoru, výsledky jsou az ve druhem radku 13 5 6 12 2 17 3 18 12 4 5 11 1 16 2 17 Stochastické modely časových řad (4. cvičení) > TrigoTrend$signFreq*n/(2*pi) [1] 12 4 5 11 1 16 2 17 Významné frekvence ve tvaru ^: > pdgV$kF/n 13 5 6 12 2 17 0.052863436 0.017621145 0.022026432 0.048458150 0.004405286 0.070484581 3 18 0.008810573 0.074889868 > TrigoTrend$signFreq/(2*pi) [1] 0.052863436 0.017621145 0.022026432 0.048458150 0.004405286 0.070484581 [7] 0.008810573 0.074889868 Významné úhlové frekvence (tvar 27r^): > 2*pi*pdgV$kF/n 13 5 6 12 2 17 3 0.33215077 0.11071692 0.13839615 0.30447153 0.02767923 0.44286769 0.05535846 18 0.47054692 > TrigoTrend$signFreq [1] 0.33215077 0.11071692 0.13839615 0.30447153 0.02767923 0.44286769 0.05535846 [8] 0.47054692 Významné periody ve tvaru ^: > n/pdgV$kF 13 5 6 12 2 17 3 18 18.91667 56.75000 45.40000 20.63636 227.00000 14.18750 113.50000 13.35294 > 2*pi/TrigoTrend$signFreq [1] 18.91667 56.75000 45.40000 20.63636 227.00000 14.18750 113.50000 [8] 13.35294 12 Stochastické modely časových řad (4. cvičení) Nyní s využitím funkce trigoval () vykreslíme pro rezidua do grafu trigonometrický trend se zjištěnými významnými periodami. Za účelem vytvoření hladké křivky nejprve vytvoříme vektor TTIME se 300 ekvidistantními časovými okamžiky. Pak na něj aplikujeme již zmíněnou funkci trigoval, která vyžaduje ještě další tři argumenty: A - vektor odhadů ctj, B - odhady koeficientů 0j, freq - vektor s významnými úhlovými frekvencemi. Dosadíme do nich výsledky funkce trigofit(). Výsledkem funkce trigoval() pak bude vektor s hodnotami odhadnutého modelu s trigonometrickým trendem pro zadané časové okamžiky. Výsledky znázorníme graficky obvyklým způsobem. > TTIME<-seq(1, length (res) , length=nn) > predTrigoTrend<-trigoval(TTIME,A=TrigoTrend$signA,B=TrigoTrend$signB,freq=TrigoTrend$signFreq) > par(mar=c(2,2,l,0)+0.5) > plot(res~x, type="l",mELÍn=TXT,cex.mELÍn=0.85) > Unes (new$x,predTrigoTrend, col="dodgerblue ",lwd=2) Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980 co o 1750 1800 1850 1900 1950 Obrázek 4: Trigonometrický trend na reziduích pro šířky letokruhů Na tatáž data zkusíme aplikovat iterativní metodu Damsletha a Spj0tvolla a porovnat výsledky obou metod. K tomu použijeme funkci damslethO (další funkce ze souboru FunkceM5201 .R). Jako vstup musíme této funkci zadat data, pro která má být model odhadován (v našem případě rezidua z modelu s lineárním trendem), a hladinu významnosti pro Fisherův test významných period (argument pstlevel). > betaD<-damsleth(res,pstlevel=0. 05) > str(betaD) List of 3 $ A : Named num [1:7] -0.0818 -0.0478 -0.0182 0.0454 0.0715 ... ..- attr(*, "names")= chr [1:7] "A" "A" "A" "A" ... $ B : Named num [1:7] 0.01663 0.06242 0.07955 0.05614 -0.00469 ... ..- attr(*, "names")= chr [1:7] "B" "B" "B" "B" ... $ freq: Named num [1:7] 0.3367 0.1107 0.1432 0.3094 0.0373 ... ..- attr(*, "names")= chr [1:7] "freq" "freq" "freq" "freq" ... Stochastické modely časových řad (4. cvičení) 13 Funkce vrací jako výsledek seznam se třemi vektory: f req - významné úhlové frekvence, A a B - odhady koeficientů a j a [3j pro trigonometrický trend s těmito periodami. Podívejme se, k jakým významným periodám dospěl iterační algoritmus. > (DamslethPeriods<-2*pi/betaD$freq) freq freq freq freq freq freq freq 18.66019 56.77548 43.87882 20.30496 168.56434 14.27805 13.41903 Nyní pro rezidua graficky znázorníme výsledky klasické i iterativní metody do jednoho grafu. > predTrigD<-trigovELl (TTIME, A=betaĽ$A,B=betaĽ$B,íreq=betaD$íreq) > par(mar=c(2,2,l,0)+0.5) > plot(res~x, type="l",mELÍn=TXT,cex.mELÍn=0.85) > Unes (new$x,predTrigoTrend, col="dodgerblue ",lwd=2) > Hnes(new$x,predTrigD,lwd=2, col="red",lty=2) > TFP<-paste("Fisher's periods: ", paste (round(FisherPeriods, 1) , collapse=", ")) > TDS<-paste("Damsleth-Spjotvoll: ", paste (round(DamslethPeriods, 1) , collapse=", ")) > legend("topright", legend = c(TFP,TDS) ,bty="n", lty = c(l,2), xjust = 1, yjust = 1, col=c("dodgerblue", "red")) Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980 CO o o cg o o o o l Fisher's periods: 18.9, 56.8, 45.4, 20 Damsleth-Spjotvoll: 18.7, 56.8, 43.9, 20 6, 227, 14.2, 113.5, 13.4 |3, 168.6, 4.3, 13.4 1-1-1-1-1- 1750 1800 1850 1900 1950 Obrázek 5: Trigonometrický trend (metoda skrytých period i iterativní metoda) na reziduích pro šířky letokruhů. Na závěr přidáme i lineární trend a vše vykreslíme pro původní hodnoty časové řady. > pred.lintrend<-predict(LinTrend, new) > par(mar=c(2,2,l,0)+0.5) > plot(y~x, type="l",main=TXT, cex.main=0.85) > Unes (new$x,predTrigoTrend+pred. lintrend, col="'dodgerblue ",lwd=2) > Unes(new$x,predTrigD+pred.lintrend,lwd=2, col="red",lty=2) > Unes (new$x,pred.lintrend, lwd=l .5, col="darkgreen",lty=2) > legend("topright", legend = c(TFP,TDS, "linear trend") , bty="n", cex=0. 8, lty = c(l,2,2), xjust = 1, yjust = 1, col=c("dodgerblue", "red", "darkgreen")) 14 Stochastické modely časových řad (4. cvičení) CD CM 00 O CD Tree rings of Gulmarg, Kashmir (Himalayan Fir), Years: 1753-1980 • Fisher's periods: 18.9, 56.8, 45 Damsleth-Spjotvoll: 18.7, 56.8, 43. — linear trend T X, 20.6, 227, 14.2, 113.5, 13.4 ), 20.3, 168.6, 14.3, 13.4 T T T T 1750 1800 1850 1900 1950 Obrázek 6: Lineární i trigonometrický trend (metoda skrytých period i iterativní metoda) pro šířky letokruhů odpovídajících rokům 1753-1980 3 Ukol Datový soubor andrews2. dat obsahuje údaje o ročním počtu získaných rysích kožešin týkající se kanadské společnosti Hudson's Bay company. • Datový soubor načtěte. • Časovou řadu vykreslete do grafu. • Vypočítejte periodogram. (Předtím nezapomeňte z řady odstranit lineární trend.) • Pokud naleznete nějaké významné periody: — Vykreslete do grafu časovou řadu (po odečtení lineárního trendu) spolu s odhadnutým trigonometrickým trendem. — Vyzkoušejte také odhadnout parametry v nelinární model s využitím algoritmu, který navrhli Damsleth a Spj0tvoll a výsledky znázorněte graficky — Regresní křivky získané oběma metodami vyneste do jednoho grafu, pracujte při tom jak s časovou řadou po odstranění lineárního trendu, tak s časovou řadou bez odstraněného lineárního trendu.