Stochastické modely časových řad
1
M5201 - 10. CVIČENÍ: Modelování sezónnosti
1   Metoda malého trendu
Uvažujme regresní model ve tvaru:
Yt = Trt + Szt + et
et ~ WN(0,a2).
j = 1,..., r       r ... počet sezón k = 1,..., d      d... délka sezóny.
Přeindexujeme       Y±,... ,Yn   na Yjk,
Předpokládejme, že trend je konstantní pro j-tou sezónu,
a rovněž sezónní hodnota je konstantní pro fc-tou sezónní složku
Regresní model můžeme napsat ve tvaru
ejk~ WN(0,a2)
tj. Tr-j = m j tj. Szk = sk.
Mi
y3k
m j + sk + £jk
j = 1, • • • ,r
Maticově lze tento model rozepsat takto
fYn\
Y2i
	0			0	1	0		0
					0	1		
								0
1	0			0	0		0	1
0	1	0 ■■■		0	1	0		0
					0	1		
								0
0	1	0 ■■■		0	0		0	1
0	0	0 ■■■	0 1	0	1	0		0
					0	1		
								0
0	0	0 ■■■	... 1	0	0		0	1
0	0	0 ■■■	... 0	1	1	0		0
					0	1		
								0
0	0	0 ■■■	... 0	1	0		0	1
k = 1,
MA
( mi\
Sl
\SdJ
+
£21
£'2d
£rl
£rd \ )
Blokově lze psát
/ Yi \     / ld 0 0
\ 0
■    ■■■ 0	Id \
0	
0 ld	la )
I mi \
( £l \
si
+
\ sd J
\ £r J
označme XM/=X= (X(1j |X(2))
kde ld je sloupcový vektor délky d samých jedniček a 1^ je diagonální jednotková matice typu d x d.
2
Stochastické modely časových řad (10. cvičení)
Matice plánu Xjví7 = X = ^X^-JX^J však není plné hodnosti, neboť když sečteme prvních r sloupců, dostaneme vektor samých jedniček, což je rovno také součtu posledních d sloupců. Aby měl model plnou hodnost, přidejme ještě jednu podmínku, a to
Si-\-----h sd = 0
Potom
X'X :
/ d 0
0 d
0 ■■■
1 ■■■
'•• 0
0_d_
77. f
v \......
kde využíváme tzv. tečkové notace
1 ■■■
r Ô~
0 r
0 ■■■
1 \
■•• 0
0    r J
a    X'Y:
Y.i
i=i
Y,
ík-
i=l
, r, k = 1,..., d) do tvaru
Normálni rovnice X'X/3 = X'Y můžeme přepsat (pro j = 1,
dnij + EiLi «i = Yj. =0
ELi mi + rsk = y k
PŘÍKLAD 1
Typický příklad sezónního chování můžeme vysledovat například u návštěvnosti zoologických zahrad, kde je obvykle návštěvnost v zimních měsících poměrně malá, zatímco v létě velká. Proto zkusíme metodu malého trendu aplikovat na data v souboru zoo_jihlava.txt, v němž jsou uvedeny měsíční počty návštěvníků jihlavské ZOO v letech 2006-2010. Datový soubor načteme pomocí příkazu scan(). Příkazem str() vypíšeme strukturu načtených dat.
> íileDat <- paste(data.library,   "zoo_jihlava.txt", sep = "")
> zoo <- scan (ÍileDat)
> str(zoo)
num [1:60] 665 1309 2093 15624 33531 ...
Pomocí funkce ts() vytvoříme z vektoru časovou řadu. Příkazem str() se podíváme na její strukturu.
> zooTS <- ts(zoo, start = 2006, frequency = 12)
> str(zooTS)
Time-Series [1:60] from 2006 to 2011: 665 1309 2093 15624 33531
Stochastické modely časových řad
3
Časovou radu vykreslíme pomoci príkazu plot(). > plot(zooTS,  type = "o", pch = 20, cex = 3)
	O
	O
	O -
	O
	CO
	o
	o
	o -
	o
	LO
	o
	o
	o -
	o
	
<z>	o
1-	o
o	o -
o	o
N	CO
	o
	o
	o -
	o
	(M
	o
	o
	o -
	o
	
	o -
2006
2007
2008
2009
2010
2011
Time
Obrázek 1: Návštěvnost v ZOO Jihlava v letech 2006-2010
Pracujeme-li v prostředí R s regresním modelem pro kategóriami proměnné, je třeba specifikovat přesnou podobu matice plánu pomocí tzv. kontrastů.
Chceme-li v našem příkladu použít model
'jk = mj + sk + £jk      £jk ~ WN(0, a2)
i = i,
k = 1,..., d
s dodatečnou podmínkou
Si-\-----h sd = 0,
budeme muset odpovídajícím způsobem nastavit kontrasty.
Protože na roční úrovně (parametry mi,... ,mr) nejsou kladeny žádné dodatečné podmínky a matici plánu jsme (dost neobvykle) navrhli tak, že neobsahuje vektor jedniček, tak pro parametry mi,..., mr platí
/ 1    0   ■■■   0 \ / mi \     ( ml\
\ mr J
4
Stochastické modely časových řad (10. cvičení)
Vidíme, že matice kontrastů vztahující se k mi,..., mr bude jednotková diagonální matice řádu r x r.
Na rozdíl od ročních úrovní rrij jsou sezónní kolísání vázána podmínkou
Si-\-----h sd = 0.
Takže jednu složku můžeme vyjádřit pomocí ostatních, nejčastěji jde o poslední, tj.
Sd = -si-----Sd-i.
Díky tomu můžeme vypustit poslední složku. Celou reparametrizaci lze vyjádřit maticově takto
/	1     0 ••	.   ... o			
	o   '•. '•			í     51 \	í
					
		.   ■•. 0		V Sd-1 1	Sd-1
	o .....	■    0 1			V 51-----
\	-1 -1 ■■	■  -1 -1	)		
\
contr.sum
a matice typu d x (d—1), kterou jsme označili jako contr. sum představuje kontrasty pro měsíční úrovně. Tato matice je již v prostředí R v kontrastech předem nadefinovaná.
Vrátíme se k načteným datům. Vytvoříme novou proměnnou typu data-frame, kam uložíme načtená data spolu s dalšími dvěma kategoriálními proměnnými grY (rok pozorování) a grS (měsíc pozorování). Takto vytvořený datový rámec budeme používat při práci s regresním modelem.
> xts <- zooTS
> d <- frequency(xts)
> n <- length(xts)
> m <- ceiling(n/d)
> Shift <- start (xts) [2] - 1
K vytvoření vektoru, který by udával rok odpovídající pozorováním, použijeme funkci gl (n,k, length). Tato funkce vytvoří vektor kategoriálních proměnných (faktorů) s n úrovněmi, každou z nich použije fc-krát a celková délka vektoru je dána parametrem length.
> Season <- íactor(as.integer(gl (d, 1, n)) - shiít)
> Year <- íactor(íloor(time(xts)))
> data <- data.írame(x = as.vector(xts), grS = Season, grY = Year)
> str(data)
'data.frame': 60 obs. of    3 variables:
$ x    : num   665 1309 2093 15624 33531 ...
$ grS: Factor w/ 12 levels "1","2","3","4",..:  123456789 10 ... $ grY: Factor w/ 5 levels "2006","2007" ,..: 1111111111...
Pomocí funkce lm() odhadneme regresní model. Vysvětlujícími proměnnými budou ka-tegoriální proměnné rok (grY) a měsíc (grS). Při zadávání modelu nesmíme zapomenout potlačit v matici plánu první sloupec jedniček, který by se jinak nastavil automaticky, a taky správně nastavit kontrasty.
> Ml <- lm(x ~ grY + grS - 1, data, contrasts = list (grY = diag(l, length (levels (data$grY))),
grS = contr.sum))
> summary (Ml)
Stochastické modely časových řad
5
Call:
lm(formula = x     grY + grS - 1, data = data, contrasts = list(grY = diag(l, length(levels(data$grY))), grS = contr.sum))
Residuals:
Min IQ Median 3Q Max
-10486   -1835       397     1649 9116
Coefficients:
	Estimate Std.	Error	t value		PrOltl)	
;rY2006	19113	1225	15.	604	< 2e-16	***
;rY2007	22089	1225	18.	034	< 2e-16	***
;rY2008	20363	1225	16.	624	< 2e-16	***
;rY2009	22274	1225	18.	185	< 2e-16	***
;rY2010	18901	1225	15.	431	< 2e-16	***
;rSl	-19271	1817	-10.	607	1.05e-13	***
;rS2	-17776	1817	-9.	784	1.30e-12	***
;rS3	-14810	1817	-8.	152	2.44e-10	***
;rS4	4326	1817	2.	381	0.0216	*
;rS5	13649	1817	7.	513	2.04e-09	***
;rS6	13757	1817	7.	572	1.67e-09	***
;rS7	31617	1817	17.	403	< 2e-16	***
;rS8	36113	1817	19.	878	< 2e-16	***
;rS9	-1504	1817	-0.	828	0.4122	
;rS10	-10569	1817	-5.	817	6.26e-07	***
;rSll	-17500	1817	-9.	632	2.09e-12	***
Signif. codes:    0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '  ' 1
Residual standard error: 4243 on 44 degrees of freedom Multiple R-squared: 0.9835, Adjusted R-squared: 0.9775
F-statistic:  163.9 on 16 and 44 DF,    p-value: < 2.2e-16
Vidíme, že kromě jednoho jsou všechny koeficienty vysoce signifikantní, což znamená že se statisticky významně liší od nuly (testováno pomocí statistik t).
Také p-hodnota u statistiky F naznačuje významnost modelu jako celku.
Jestliže pracujeme s regresním modelem, kde matice plánu nemá první sloupec tvořený jednotkami (ve výpisu chybí proměnná (Intercept)), koeficient determinace (Multiple R-squared), popř. upravený koeficient determinace (Adjusted R-squared) nemá interpretaci. Proto je lepší takovému modelu se spíše vyhnout.
Chceme-li posoudit významnost faktorů grY a grS ne po jednotlivých složkách (jak to nabízí jednotlivé í-testy) ale jako celek, musíme použít F-testy, jejichž hodnotu získáme použitím funkce anova().
> anova(Ml)
Analysis of Variance Table
Response: x
Df Sum Sq       Mean Sq F value Pr(>F)
grY 5 2.5454e+10 5090860284   282.77 < 2.2e-16 ***
6
Stochastické modely časových řad (10. cvičení)
grS 11 2.1751e+10 1977338565   109.83 < 2.2e-16 ***
Residuals 44 7.9215e+08 18003488
Signif. codes:    0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '  ' 1
Výsledky F testů ukazují, že jak faktor roku (grY), tak faktor sezóny (grS) jsou statisticky významné (p-hodnoty v řádcích označených grY a grS jsou menší než hladina významnosti 0.05).
Abychom zjistili odhady koeficientů
rrij       (j = l,...,r)        a       sk       (k = l,...,d)
použijeme funkci predictO s volbou type="terms". Protože tentokrát neuvádíme parametr newdata (datový rámec s hodnotami vysvětlujících proměnných, pro které se má predikce počítat), počítají se predikce z dat, která vstupovala do modelu. Získanou matici Ml.terms si podrobně prohlédneme.
> Ml.terms <- predict(Ml, type = "terms")
> str(Ml.terms)
num [1:60, 1:2]  19113 19113 19113 19113 19113 ...
- attr(*,  "dimnames")=List of 2
..$ : chr [1:60]  "1" "2" "3" "4"  ... ..$ : chr [1:2]  "grY" "grS"
- attr(*,  "constant")= num 0
Ukážeme, že fitované hodnoty jsou součtem dvou sloupců matice Ml.terms označených grY a grS .
> pom <- cbind (Ml. terms, Ml. terms [, 11 + Ml.terms [, 2], fitted(Ml))
> colnames(pom) <- c(colnames(Ml .terms) ,   "grY+grS",   "fitted(Ml) ")
> print (pom [1:6, ])
grY grS
1 19112.75 -19271.267
2 19112.75 -17775.467
3 19112.75 -14810.067
4 19112.75 4326.333
5 19112.75 13649.133
6 19112.75 13757.333
grY+grS fitted(Ml) -158.5167 -158.5167 1337.2833 1337.2833 4302.6833 4302.6833 23439.0833 23439.0833 32761.8833 32761.8833 32870.0833 32870.0833
Vidíme, že ve 3. a 4. sloupci matice jsou stejná čísla, tudíž je zřejmé, že hodnoty vyrovnané na základě modelu Mj jsou dány součtem odpovídajících odhadnutých parametrů m j a s^.
Ze všech hodnot matice Ml .terms vypíšeme ty, keré se týkají koeficientů rrij a s^. Sloupec grY udává roční úroveň příslušející každému pozorování, tzn. že se vždy opakují 12-krát po sobě. Koeficienty rrij tudíž získáme tak, že vypíšeme pouze 1., 13., 25. .. .hodnotu z prvního sloupce matice.
> (Ml.mj <- Ml.terms[seq(l, n, d), 11)
Stochastické modely časových řad
7
1 13 25 37 49
19112.75 22088.58 20362.67 22273.58 18900.75
Koeficienty najdeme ve sloupci grS, kde jsou měsíční odchylky od průměrné roční úrovně, a opakují se pro pozorování v každém roce. Koeficienty s& proto dostaneme tak, že vypíšeme prvních 12 hodnot z druhého sloupce.
> (Ml.sk <- Ml.terms[l:d, 21)
12345678 -19271.267 -17775.467 -14810.067     4326.333   13649.133    13757.333   31616.933 36112.933
9 10 11 12
-1504.067 -10568.667 -17499.667 -18033.467
Sezónní kolísání graficky znázorníme.
> SK <- Ml.sk
> nf <- layout (matrix(c(l, 2), nrow = 1), width = c(2.75, 1.25))
> txtS <- paste("S", 1 -.length(SK),  " = ", sep = "")
> par (mar = c(2, 2, 0, 0) + 0.05)
> plot(l:length(SK), SK, type = "o")
> plot(c(0, 1), c(0, length (SK)), type = "n", axes = FALSE)
> text(rep(0, length (SK)) , rev(l: length(SK)) , txtS, adj = c(0, 1) , cex = 1.125)
> text(rep(l, length (SK)) , rev(l: length(SK)) , round (SK) , adj = c(l, 1) ,
cex = 1.125)
51 = -19271
52 = -17775 S3= -14810 S4 = 4326 S5= 13649 S6= 13757 S7= 31617 S8= 36113 S9= -1504 S10= -10569 S11 = -17500 S12= -18033
1 i        i        i        i i
i
2 4     ^  6 8 10 12
Obrázek 2: Sezónní složky získané metodou malého trendu pro data Návštěvnost v ZOO
Jihlava v letech 2006-2010
Výsledky celé dekompozice vykreslíme do jediného grafu.
> tt <- as. vector(time(xts))
> xx <- as.vector(xts)
> fit <- fitted(Ml)
> tr <- rep(Ml.mj, each = d, length.out = n)
> ylim <- range (range (xx) , range (fit) , tr)
> par(mfrow = c(l, 1), mar = c(2, 2, 2,  0) +0.5)
> plot(tt, xx, type = "o", pch = 20, col = "black", ylim = ylim)
> lines(tt, fit, type = "o", pch = 22, col = "red", bg = "yellow", cex = 0.85)
> lines(tt, tr, type = "s",  col = "dodgerblue", lwd = 2)
> ábline(v = as. numeric (as. character (levels (data$grY) )) , col = "gray",
lty = 2)
> legend(par("usr") [1], par("usr")[4], bty = "n", legend = c("puvodni",
"odhady"), lty = c(l, 1), pch = c(20, 22), bg = c("black", "yellow"), col = cC'black", "red"))
8
Stochastické modely časových řad (10. cvičení)
o o co
o o o
o o o
CM
i i i i i r
2006 2007 2008 2009 2010 2011
Obrázek 3: Metoda malého trendu pro data Návštěvnost v ZOO Jihlava v letech 2006-2010
Pokusme se výchozí model Mj modifikovat tak, aby matice plánu měla v prvním sloupci samé jedničky.
Toho lze dosáhnout, budeme-li uvažovat následující model
WN(0,a2),
M
modi f
Yjk = lí + rrij + Sk+ £jk
i = i,
, r, k = 1,
, d
Přidáním nového parametru /i (celkový průměr) dostaneme opět model s neúplnou hodnosti. Protože budeme chtít jediné řešení, přidáme dvě doplňující podmínky, a to
Si-\-----h sd = 0
mi + • • • + mr
0.
Odhad parametrů regresního modelu dostaneme obdobným způsobem jako u modelu Mj, pouze změníme kontrasty pro kategoriální proměnnou grY a nevynecháme konstantní člen. Protože doplňující podmínka pro m j je analogická podmínce pro s^, je i tvar matice kontrastů stejný (až na počet řádků a sloupců), v R se jedná o kontrasty contr.sum.
> Mim <- lm(x ~ grY + grS, data,  contrasts = UstCgrY = contr.sum, grS = contr.sum))
> summary (Mim)
Call:
lm(formula = x     grY + grS, data = data, contrasts = list(grY = contr.sum, grS = contr.sum))
Residuals:
Min 1Q Median 3Q Max
-10486   -1835       397     1649 9116
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept)    20547.7 547.8   37.511   < 2e-16 ***
grYl -1434.9 1095.6   -1.310 0.1971
grY2 1540.9 1095.6     1.407 0.1666
Stochastické modely časových rad
9
grY3	-185	,0	1095.	6 -0	,169	0.8667	
grY4	1725	,9	1095.	6 1	,575	0.1223	
grSl	-19271	,3	1816.	8 -10	,607	1.05e-13	***
grS2	-17775	,5	1816,	8 -9	,784	1.30e-12	***
grS3	-14810	,1	1816.	8 -8	,152	2.44e-10	***
grS4	4326	.3	1816.	8 2	,381	0.0216	*
grS5	13649	1	1816.	8 7	513	2.04e-09	***
grS6	13757	3	1816.	8 7	,572	1.67e-09	***
grS7	31616	9	1816.	8 17	403	< 2e-16	***
grS8	36112	9	1816.	8 19	878	< 2e-16	***
grS9	-1504.	,1	1816.	8 -0	828	0.4122	
grSlO	-10568.	7	1816.	8 -5	817	6.26e-07	***
grSU	-17499.	7	1816.	8 -9.	632	2.09e-12	***
Signif.	codes:    0 '	* * * }	0.001	'**' 0.01		'*' 0.05	C J
Residual standard error: 4243 on 44 degrees of freedom Multiple R-squared: 0.965, Adjusted R-squared: 0.9531
F-statistic: 80.99 on 15 and 44 DF,    p-value: < 2.2e-16
> anova(Mlm)
Analysis of Variance Table
Response: x
Df Sum Sq       Mean Sq   F value Pr(>F)
grY 4 1.2191e+08     30476274     1.6928 0.1687
grS 11 2.1751e+10 1977338565 109.8309 <2e-16 ***
Residuals 44 7.9215e+08 18003488
Signif. codes:    0 '***> 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '  ' 1
Parametry mi,..., mr v tomto případě mají jinou interpretaci. Nejsou roční hladinou, ale ukazují kolísání kolem průměrné hladiny dané parametrem (Intercept), který je odhadem parametru /x v modelu M™odif _
Pro získání koeficientů
m j       (j = l,...,r)        a       sk       (k = l,...,d)
opět použijeme funkci predict() s volbou parametru type="terms" a bez volby newdata. takže se predikce počítají z dat, která vstupovala do modelu.
> Mim.terms <- predict(Mim,  type = "terms")
> str(Mim.terms)
num [1:60, 1:2] -1435 -1435 -1435 -1435 -1435 ...
- attr(*,  "dimnames")=List of 2
..$ : chr [1:60]  "1" "2" "3" "4"  ... ..$ : chr [1:2]  "grY" "grS"
- attr(*,  "constant")= num 20548
10
Stochastické modely časových řad (10. cvičení)
Ukážeme, že tentokrát (protože model obsahuje (Intercept)) jsou fitované hodnoty součtem dvou sloupců matice Mim.terms označených grY a grS a konstanty, kterou získáme příkazem attr(Mim.terms, "constant").
> pom <- cbind (Mim. terms, Mim. terms [, 1] + Mim. terms í, 2] + attr (Mim. terms,
"constant"), fitted(Mlm))
> colnames(pom) <- c(colnames(Mlm. terms) ,   "grY+grS+constant",   "iitted(Mlm) ")
> print (pom [1:6, ])
		grY		grS	grY+grS+constant		fitted(Mlm)	
1	-1434.	917	-19271	267	-158	5167	-158	5167
2	-1434.	917	-17775	467	1337	2833	1337	2833
3	-1434.	917	-14810	067	4302	6833	4302	6833
4	-1434.	917	4326	333	23439	0833	23439	0833
5	-1434.	917	13649	133	32761	8833	32761	8833
6	-1434.	917	13757	333	32870	0833	32870	0833
Vybereme nyní příslušné koeficienty s& a /x + m j > (Mlm.sk <- Mim. terms[1:d, 21)
12345678 -19271.267 -17775.467 -14810.067     4326.333   13649.133    13757.333   31616.933 36112.933
9 10 11 12
-1504.067 -10568.667 -17499.667 -18033.467
> (Mlm.mj <- attr(Mim.terms,   "constant") + Mim.terms[seq(1, n, d), 11)
1 13 25 37 49
19112.75 22088.58 20362.67 22273.58 18900.75
Sezónní kolísání     vykresleme opět do grafu.
> SK <- Mlm.sk
> ní <- layout (matrix(c(1, 2), nrow = 1), width = c(2.75, 1.25))
> txtS <- paste("S", 1 -.length(SK),  " =", sep = "")
> par (mar = c(2, 2, 0, 0) + 0.05)
> plot(l:length(SK), SK, type = "o")
> plot(c(0, 1), c(0, length (SK)), type = "n", axes = FALSE)
> text(rep(0, length (SK)) , rev(l: length(SK)) , txtS, adj = c(0, 1) , cex = 1.125)
> text(rep(l, length (SK)) , rev(l: length(SK)) , round (SK) , adj = c(l, 1) ,
cex = 1.125)
51 = -19271
52 = -17775 S3= -14810 S4 = 4326 S5= 13649 S6= 13757 S7= 31617 S8= 36113 S9= -1504 S10= -10569 S11 = -17500 S12= -18033
2 4 6 8 10 12
Stochastické modely časových řad
11
Obrázek 4: Sezónní složky získané modifikovanou metodou malého trendu pro data Návštěvnost v ZOO Jihlava v letech 2006-2010
Výsledky celé dekompozice vykreslíme do jediného grafu.
> fit <- fitted(Mijnj)
> tr <- rep (Mim.mj, each = d, length.out = n)
> ylim <- range (range (xx) , range (fit) , tr)
> par (mírow = c(l, 1), mar = c(2, 2, 2,  0) +0.5)
> plot(tt, xx, type = "o", pch = 20, col = "black", ylim = ylim)
> lines(tt, fit, type = "o", pch = 22, col = "red", bg = "yellow", cex = 0.85)
> lines (tt, tr, type = "s",  col = "dodgerblue", lwd = 2)
> abline(v = as. numeric (as. character (levels (data$grY) )) , col = "gray",
lty = 2)
> ábline(h = attr(Mim. terms,   "constant"), lty = 2, col = "gray35")
> legend(par("usr") [1], par("usr")[4], bty = "n", legend = c("puvodni",
"odhady"), lty = c(l, 1), pch = c(20, 22), bg = c("black", "yellow"), col = cC'black", "red"))
o o o o co
o o in
o o o
o o o o
1 \ \ \ \ T
2006 2007 2008 2009 2010 2011
Obrázek 5: Modifikovaná metoda malého trendu pro data Návštěvnost v ZOO Jihlava v letech 2006-2010
Na závěr se můžeme porovnáním vyrovnaných hodnot pro model Mj a Mj že obě dvě metody dávají úplně stejné odhady.
modi f
přesvědčit,
> cbind (fitted (Ml), fitted (Mim)) [1:6, ]
[,1] [,2]
1 -158.5167 -158.5167
2 1337.2833 1337.2833
3 4302.6833 4302.6833
4 23439.0833 23439.0833
5 32761.8833 32761.8833
6 32870.0833 32870.0833
12
Stochastické modely časových řad (10. cvičení)
Pozn.: Obě dvě metody použité v Příkladu 1 jsou naprogramovány ve funkcích SzSmallTrendO a SzSmallTrendModif (), které jsou uloženy v souboru FunkceM5201 .R. Po načtení příkazem source šije můžete prohlédnout, pokud v R zadáte SzSmallTrend, resp. SzSmallTrendModif.
2   Modelování sezónnosti pomocí SARIMA modelů
Sezónnost je v Box-Jenkinsově metodologii stejně jako trend modelována stochasticky K tomu se zavádí sezónní diferenční operátor o délce L > 0:
AL = Yt - Yt-L = (1 - BL)Yt A| = AL(ALyt) = (1 - BL)2Yt
A£ = (1 - BL)DYt
Při konstrukci SARIMA modelů se postupuje způsobem, který si přiblížíme na příkladu dat, která vykazují sezónnost s periodou o délce L = 12.
1. Vezmeme nejprve pouze lednová měření, tj. řadu {Sf = B12Yt} a zkonstruujeme pro ně model ARIMA(P1, Di,Q{)
TniB^A^Yt = *i(i312H(1) ~ ARIMA(P1,D1,Q1), kde časový index t odpovídá lednovým obdobím. Přitom
tvi(B12) = 1 — tvi^B12 — ■ ■ ■ — tvipiB12'Pi (sezónní autoregresní operátor SAR(Pi)) ^i(B12) = 1 + vpi,iB12 + • • • + vpi,q1B12'®1 (sezónní operátor klouzavých součtů SMA(Qi)) = (1 — B12)Dl (sezónní diferenční operátor SI(Di))
2. Podobné modely zkonstruujeme i pro ostatní měsíce:
7r2(B12)Ag2Yt = y2(B12)42) ~ ARIMA(P2,D2,Q2)
7r12(B12)A^Yt = *12(B12)r,ľ2) ~ ARIMA(P12,D12,Q12)
3. Dále předpokládáme, že modely pro jednotlivé měsíce jsou přibližně stejné, tj.
Pl ps • • • ps P12 ps P Qi ps • • • ps Qi2 ps Q Dl ps ••• ps D12 ps D
th(B12) ps ... PS7T12CB12) ps7t(S12)
^!(B1Z) Ps ••• Ps ^12(B1Z) Ps V(B
Stochastické modely časových řad
13
4. Náhodné veličiny ryt ,j = 1,..., 12 by však v těchto modelech měly být pro různé měsíce mezi sebou korelované, neboť by měl existovat např. vztah mezi lednovými a únorovými hodnotami. Předpokládáme proto, že také řada ijt je popsána modelem ARIMA(p, d, q) tvaru:
^(B)A% = e(B)et ~ ARIMA(p, d, q), kde st ~ WN(0, a2) je bílý šum.
5. Oba předchozí model spojíme do jediného tzv. multiplikativního sezónního modelu řádu (p,d,q) X (P,D,Q)L
$(B)7r(BL)AdA%Yt = G(B)^(BL)et ~ SARIMA(p, d, q) x (P,D,Q)L,L = 12 PŘÍKLAD 2
Budeme zkoumat SARIMA proces zadaný předpisem
SARIMA(0,1,1) x (0,1,1)12 :  (1 - B)(l - B12)Yt = (1 - 0.5B)(1 - 0.7B12)et, kde
et ~ WN(0,a2),        a2 = 1.75. Pokud model rozepíšeme, dostaneme
Yt = yt_i + yt_i2 + yt-i3 + et - 0.5£í_i - 0.7£í_i2 + 0.35£í_i3
A nyní nasimulujeme časovou řadu s takto zadanými parametry. Použijeme k tomu další funkci ze souboru FunkceM5201 .R, která se jmenuje my.sarima.sim. Použijeme tyto argumenty:
• n=25 ... 25 sezón (let)
• period=12 .. .pro každou sezónu chceme 12 pozorování (odpovídá měsícům v roce)
• sd=sqrt (1.75) ... směrodatná odchylka bílého šumu a
(i)
• model=list(ma=-0.5,order=c(0,1,1)) ...parametry ARIMA procesu pro řadu ryt
• seasonal=list (ma=-0.7, order=c (0,1,1)) ... parametry ARIMA procesu pro časové řady odpovídající zvlášť jednotlivým měsícům
> fileFun <- paste(data.library,   "FunkceM5201.R", sep = "")
> source(fileFun)
> Sarima. sim <- my. sarima. sim(n = 25, period = 12, sd = sqrt(1.75), model = listOna = -0.5,
order = c(0, 1,  1)), seasonal = listOna = -0.7, order = c(0, 1, 1)))
Vygenerovanou časovou řadu vykreslíme do grafu.
> plot(Sarima.sim)
14
Stochastické modely časových řad (10. cvičení)
o
V dalším kroku by bylo dobré si prohlédnout výběrovou autokorelační funkci, výběrovou parciální autokorelační funkci a periodogram. Protože máme v souboru FunkceM5201 .R pro tento případ naprogramovanou funkci eda.ts, můžeme všechny grafy vykreslit naráz v jednom kroku.
> eda.tsCSarima.sim,  bands = TRUE)
Stochastické modely časových řad
15
5 10 15 20 25
1.0 -0.8 -0.6 -
LL
Obrázek 7: ACF, PACF a periodogram pro SARIMA(0,1,1) x (0,1,1)12 : (1 - B)(l -B12)Yt = (1 - 0.55)(1 - 0.7B12)et
U autokorelační funkce můžeme najít pro k = 12 lokální maximum a podobné chování by se mělo objevit v dalších násobcích čísla 12 (to už ale v grafu nevidíme). I u parciální autokorelační funkce jsou patrné vyšší hodnoty v bodech k = 12 a k = 24.
Předpokládejme, že máme k dispozici pouze časovou řadu Sarima.sim a nevíme, jakým procesem byla vygenerována. Provedeme pro ni odhad neznámých parametrů nejprve pomocí funkce arima().
> sarima.fit <- arijnaCSarijna.sijn, order = c(0,  1, 1), seasonal = list(order = c(0,
1, 1), period = 12))
> print(sarima..fit)
16
Stochastické modely časových řad (10. cvičení)
Series: Sarima.sim
ARIMA(0,1,1)(0,1,1)[12]
Call: arima(x
Sarima.sim, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1),
period = 12))
Coefficients: mal
smal
-0.6588 -0.8473 s.e.      0.0473 0.0450
sigma~2 estimated as 4.178: log likelihood = -620.3 AIC = 1246.6     AICc = 1246.68     BIC = 1257.57
Ve statistické knihovně máme ještě k dispozici funkci tsdiagO, kterou použijeme na odhadnutý SARIMA model, abychom ověřili, že rezidua již připomínají bílý šum.
> tsdiagCsarima.fit)
Stochastické modely časových řad
17
p values for Ljung-Box statistic
niii
2 4 6 8 10
Obrázek 8: Diagnostické grafy pro odhadnutý model SARIMA(0,1,1) x (0,1, l)i2 : (1 B)(l - B12)Yt = (1 - 0.55X1 - 0.7B12)et
Provedeme predikci na 2 roky dopředu, a to pomocí příkazu PlotPredictARIMAO. Výsledky zakreslíme do grafu.
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> PlotPredictARIMA(Sarima.sim, sarima.fit, n.ahead = 24)
18
Stochastické modely časových řad (10. cvičení)
10
5 -
0 -
-10
1-1-T
0 5 10 15 20
Obrázek 9: Predikce pro odhadnutý model SARIMA(0,1,1) x (0,1,1)12 B12)Yt = (1 - 0.5B)(1 - 0.7B12^
r
25
(1
B)(l
V knihovně forecast máme k dispozici funkci auto .arimaO, kterou také můžeme použít při odhadu modelu na naše simulovaná data (u této funkce nemusíme předem znát řád SARIMA procesu).
> library(forecastJ
> sarima.fitF <- auto.arima(Sarima.sim,  trace = TRUE)
ARIMA(2	1	2)	(1,0,1) [12]		with	drift		1304.535
ARIMACO	1	0)	with drift			: 1452	135	
ARIMA(1	1	0)	(1,0,o;	[12]	with	drift		1363.831
ARIMACO	1	1)	(0,0, i;	[12]	with	drift		1321.387
ARIMA(2	1	2)	(0,0, i;	[12]	with	drift		1320.695
ARIMA(2	1	2)	(2,0, i;	[12]	with	drift		1307.369
ARIMA(2	1	2)	(1,0,o;	[12]	with	drift		1324.994
ARIMA(2	1	2)	(1,0,2;	[12]	with	drift		1306.999
ARIMA(2	1	2)	with drift			: 1331	400	
ARIMA(2	1	2)	(2,0,2;	[12]	with	drift		le+20
ARIMA(1	1	2)	(l,0,i;	[12]	with	drift		1295.781
ARIMA(1	1	1)	(l,0,i;	[12]	with	drift		1295.731
ARIMACO	1	0)	(l,0,i;	[12]	with	drift		le+20
ARIMAd	1	1)	(l,0,i;	[12]				1293.779
ARIMAd	1	1)	(0,0,i;	[12]				1320.990
ARIMAd	1	1)	(2,0,i;	[12]				1301.367
ARIMAd	1	1)	(1,0,0;	[12]				1320.275
ARIMAd	1	1)	(1,0,2;	[12]				1295.749
ARIMAd	1	1)				: 1336	943	
ARIMAd	1	1)	(2,0,2;	[12]				le+20
ARIMACO	1	1)	(l,0,i;	[12]				1297.169
ARIMAC2	1	1)	(l,0,i;	[12]				1303.838
ARIMACl	1	0)	(l,0,i;	[12]				le+20
ARIMACl	1	2)	(l,0,i;	[12]				1293.807
ARIMACO	1	0)	(l,0,i;	[12]				le+20
ARIMAC2	1	2)	(l,0,i;	[12]				1303.339
Stochastické modely časových řad
19
Best model: ARIMA(1,1,1)(1,0,1)[12]
> print(sarima..fitF)
Series: Sarima.sim ARIMA(1,1,1)(1,0,1)[12]
Call: auto.arima(x = Sarima.sim, trace = TRUE)
Coefficients:
arl mal sari smal
0.0028 -0.6554 0.9666 -0.8050
s.e.    0.0798 0.0651 0.0356 0.0627
sigma~2 estimated as 4.31:    log likelihood = -642.69 AIC = 1293.78     AICc = 1293.98     BIC = 1312.28
Pro predikci o 2 roky dopředu použijeme funkci forecast ().
> par(mfrow = c(l, 1), mar = c(2, 2, 1,  0) + 0.05)
> plot(íorecast(sarima.íitF), n = 24)
Forecasts from ARIMA(1,1,1)(1,0,1)M21
i-1-1-1-1-1-
0 5 10 15 20 25
Obrázek 10: Predikce na základě odhadnutého modelu SARIMA{0,1,1) x (0,1, l)i2 : (1 — B)(l - B12)Yt = (1 - 0.5B)(1 - 0.7B12)et
3 Úkoly:
Načtěte do R datový soubor rusove_prenocovani .txt. V něm jsou uvedeny čtvrtletní údaje o počtu přenocování návštěvníků z Ruska v lázeňských zařízeních v CR (vyjádřeno v počtu nocí strávených v ubytovacích zařízeních). Data se vztahují k období od 1. čtvrtletí roku 2000 do 2. čtvrtletí 2011. Na tato data aplikujte metodu malého trendu.