Stochastické modely časových řad (8. cvičení)
1
M5201 - 8. CVIČENÍ: Procesy diferenčně stacionární - ARIMA modely. Jejich generování a následné odhady
1   Procesy nestacionární ve střední hodnotě
U procesů nestacionárních ve střední hodnotě je třeba odlišovat dva různé pojmy:
• deterministický trend
• stochastický trend
Připomeňme pojem stochastického trendu a uveďme jeho souvislost s ARIMA procesy. 1.1    Stochastický trend
U ARMA procesů jsme jako podmínku stacionarity požadovali, aby všechny kořeny polynomu
$(z) = 1 - (fiZ - (f2Z - ■ ■ ■ - (fpZP
ležely vně jednotkové kružnice, tj. aby proces byl kauzální. Pokud ale nějaký kořen leží
• na jednotkové kružnici, mluvíme o procesu nestacionárním se stochastickým trendem,
• uvnitř jednotkové kružnice, mluvíme o procesu nestacionárním explozivního typu.
První případ, tedy nestacionární proces se stochastickým trendem, lze převést na stacionární proces pomocí diferencování. Pomocí operátoru zpětného chodu B, který je definován vztahem
BYt = yt_i
můžeme zavést diferenční operátor a vyjádřit 1. diferenci, 2. diferenci, ... d-tou diferenci: AYt = Yt- yt_i = (1 - B)Yt
A2Yt = A(Yt - yt_i) = Yt- 2Yt-! + Yt-2 = (1 - B)2Yt AdYt = (1 - B)dYt
Nestacionární proces se stochastickým trendem nazýváme integrovaný smíšený model
a značíme jej
ARIMA(p,d, q). Formální zápis pomocí operátoru zpětného chodu je
ARIMA(p, d, q) : $(5)(1 - B)dYt = G(B)et.
2
Stochastické modely časových řad (8. cvičení)
Provedeme-li d-krát diferencování a položíme-li
Wt = (l- B)dYt,
pak Wt je stacionární ARMA(p, q) proces.
Některé náhodné procesy, se kterými jsme se již setkali, jsou speciálními případy ARIMA(p, d, q) procesu. Přehled je uveden v následující tabulce.
p	d	q	Zkratka	Název	Formální zápis
0			IMA(d, q)	Integrovaný proces klouzavých součtů	(1 - B)dYt = e(B)et
0	0		MA{q)	Proces klouzavých součtů	Yt = e{B)et
		0	ARI(p, d)	Integrovaný autoregresní proces	- B)dYt = et
	0	0	AR{p)	Autoregresní proces	HB)Yt = et
0		0	I(d)	Integrovaný proces	(1 - B)dYt = et
0	1	0		Náhodná procházka	(1 - B)Yt = et
Operátor
v{B) = - B)d
se někdy nazývá zobecněný autoregresní operátor.
Pokud v{B) chápeme jako polynom proměnné B, pak vzhledem ke kauzalitě procesu $>(B)Wt = @(B)et má v(B) právě p kořenů ležících vně jednotkového kruhu a d kořenů rovných 1.
V praxi se nejprve <i-násobným diferencováním časové řady získá stacionární řada Wt a pro ni se pak vybuduje proces ARM A(p, q). Pokud jsme původně měli n pozorování Y±, Y2,..., Yn, po diferencování zůstane n — d pozorování Wd+i, • • •, Wn.
PŘÍKLAD 1 - ARI proces
Zaměřme se nejprve na konkrétní ARI proces řádu (1,1,0) (tzn. po jednom diferencování dostaneme AR(1) proces). Zvolme konkrétní hodnotu tp = 0.9, tj. proces lze zapsat ve tvaru
(1 - B)(l - 0.9B)Yt = et.
AR(1) proces s tímto koeficientem by měl být kauzální, protože u AR{1) procesů je podmínka pro kauzalitu splněna, jestliže \ip\ < 1. Můžeme ale ověřit i obecnou podmínku kauzality pro AR proces, tj. zda kořen polynomu 1 — tpz = 0 leží vně jednotkové kružnice.
1 - 0.9z = 0
Tento výpočet můžeme provést také pomocí R:
> ARphi <- c(0.9)
> ARroots <- polyroot(c(l, -ARphi))
> Mod(ARroots)
[1] 1.111111
Stochastické modely časových řad (8. cvičení)
3
Z předchozích cvičení víme, že kauzální proces je stacionární. ARI proces se od AR procesu liší v tom, že polynom na levé straně rovnice popisující AR proces je navíc vynásoben výrazem (1 — B), a proto má polynom na levé straně rovnice navíc ještě jeden kořen rovný 1.
ARI proces můžeme postupně rozepsat následujícím způsobem:
(l-B)(l-0.9B)Yt = et (1 - 1.95 + 0.9B2)Yt = et Yt-1.9Yt-i + 0.9Yt-2=£t
Yt = 1.9^-0.9^-2+^
Nyní vygenerujeme konkrétní realizaci ARI procesu.
> ARphi <- c(0.9)
> ari.sim <- arima.sim(model = list(order = c(l, 1, 0),
ar = ARphi), n = 300)
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> plot(ari.sim)
0 50 100 150 200 250
Obrázek 1: Simulace ARI procesu (1 — B)(l — 0.95)1^ = Et Prohlédneme si autokorelační a parciální autokorelační funkci simulovaného ARI procesu.
> vARIacf <- acf(ari.sim, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c(2, 2, 3,  0) + 0.05)
> plot(vARlací, main = "Výběrová ACF")
4
Stochastické modely časových řad (8. cvičení)
Výberová ACF
~i
10
~i
15
I
20
Obrázek 2: ACF pro simulovanou řadu ARI procesu (1 — B)(l — 0.9B)Yf = Et
> vARIpací <- pací(ari.sim, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vARIpací, main = "Výberová PACF")
Výberová PACF
1     I T
~1
10
"T-1-r-
I
15
I
20
Obrázek 3: ACF pro simulovanou řadu ARI procesu (1 — B)(l — 0.9B)Yt = et
Vidíme, že autokorelační funkce klesá velmi pomalu, což odpovídá procesům s jednotkovým kořenem.
Parametry ARI procesu odhadneme pomocí funkce arima, kde jako parametr order zadáme řád ARIMA procesu c (1,1,0).
> ari.fit <- arima (ari. sim, order = c(l, 1, 0))
> priu t (ari . f i t)
Call:
arima(x = ari.sim, order = c(l, 1, 0))
Coefficients: ari
Stochastické modely časových řad (8. cvičení)
5
0.9219 s.e. 0.0218
sigma~2 estimated as 0.8711:    log likelihood = -405.93,    aic = 815.86
Nyní, když máme odhadnuty parametry ARI procesu, využijeme jejich znalosti k predikování budoucích hodnot procesu. K tomu použijeme funkci PlotPredictARIMA, jejíž zdrojový kód musíme nejprve načíst ze souboru FunkceM5201 .R.
> fileFun <- paste (data.library,   "FunkceM5201.R", sep = "")
> source(fileFun)
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> PlotPredictARIMA (ari.sim, ari.íit, n.ahead = 10)
0 50 100 150 200 250
Obrázek 4: Predikované hodnoty pro simulovanou řadu ARI procesu (1 — B)(l — 0.9B)Yf = Et
Pokud bychom neznali přesný řád odhadovaného ARI procesu, museli provést diferencování tolikrát, dokud bychom nezískali stacionární řadu, a poté odhadnout parametry ARMA procesu pro diferencovanou řadu.
Tento postup nyní vyzkoušíme. Nejprve spočítáme první diference simulované řady.
> dif.ari <- diff(ari.sim)
> par(mírow = c(l, 1), mar = c(2, 2, 1,  0) + 0.05)
> plot(dií. ari)
6
Stochastické modely časových řad (8. cvičení)
Obrázek 5: 1. diference simulované řady ARI procesu (1 — B)(l — 0.9B)Yf = Et Prohlédneme si autokorelační a parciální autokorelační funkci diferencované řady.
> difARlELcf <- acf(dif.ari, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mař = c(2, 2, 3,  0) + 0.05)
> plot (difARlELcf, main = "Výběrová ACF")
Výběrová ACF
~i
10
~i
15
I
20
Obrázek 6: ACF pro simulovanou řadu ARI procesu (1 — B)(l — 0.9B)Yt = et
> vARIpací <- pací(dií.ari, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vARIpací, main = "Výberová PACF")
Stochastické modely časových řad (8. cvičení)
7
Výberová PACF
J_L
I    . I
10
15
I
20
Obrázek 7: PACF pro simulovanou řadu ARI procesu (1 — B){1 — 0.9B)Yf = et
Vidíme, že autokorelační funkce klesá exponenciálně k nule, zatímco parciální autokore-lační funkce je významně nenulová pouze v bodě k = 1, což by odpovídalo AR{1) procesu. Proto odhadneme parametry takového modelu pro diferencovanou řadu.
> dif.ari.fit <- arimaCdif. ari, order = c(l, 0, 0))
> print (dif. ari . fit)
Call:
arima(x = dif.ari, order = c(l, 0, 0))
Coefficients:
ari intercept
0.9201 0.3020
s.e.    0.0222 0.6509
sigma~2 estimated as 0.8706:    log likelihood = -405.83,    aic = 817.65
Z výsledku vidíme, že jsme opět dostali podobný odhad parametru tp.
příklad 2 - IMA proces
Nyní se podívejme na konkrétní případ IMA procesu řádu (0,1,2) (tzn. proces, u nějž po jednom diferencování dostaneme MA(2) proces). Zvolíme konkrérní hodnoty parametrů 9\ = 0.75 a 62 = 0.2, tj. proces lze zapsat ve tvaru
(1 - B)Yt = st + 0.75£í_i + 0.2£í_2 Yt - yt_i = et + 0.75£í_i + 0.2£í_2
Yt = Yt-! +et + 0.75£í_i + 0.2£í_2
Ještě předtím, než vygenerujeme realizaci tohoto procesu, bylo by vhodné ověřit, že polynom definující MA část procesu má kořeny ležící vně jednotkové kružnice a že MA proces, který bychom získali diferencováním, je tím pádem invertibilní.
> MAtheta <- c(0.75, 0.2)
> Mod(polyroot(c(l, MAtheta)))
8
Stochastické modely časových řad (8. cvičení)
[1] 2.236068 2.236068
Absolutní hodnota obou kořenů je větší než 1, tudíž leží vně jednotkové kružnice a podařilo se nám ověřit invertibilitu.
Analogicky jako v předchozím příkladu vygenerujeme 300 hodnot IMA(1,2) procesu a takto získanou řadu vykreslíme do grafu.
> ima.sim <- arima.sim(model = list(order = c(0, 1, 2),
ma = MAtheta), n = 300)
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> plot (ima. sim)
Obrázek 8: Simulovaná řada IMA procesu (1 - B)Yt = (1 + 0.755 + Q.2B2)et Opět si můžeme prohlédnout grafy výběrové autokorelační a parciální autokorelační funkce.
> vIMAací <- acf Cima. sim, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vIMAací, main = "Výberová ACF")
Stochastické modely časových řad (8. cvičení)
9
Výberová ACF
~i
10
~i
15
I
20
Obrázek 9: ACF simulované řady IMA procesu (1 - B)Yt = (1 + 0.755 + 0.2B2)et
> vIMApacf <- pacfdma.sim, lag.max = 20, plot = FALSE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vIMApacf, main = "Výberová PACF")
Výberová PACF
I    I 1
—1
10
i
15
I
20
Obrázek 10: PACF simulované řady IMA procesu (1 - B)Yt = (1 + 0.75S + 0.2B2)et
Opět vidíme, že autokorelační funkce klesá k nule velmi pomalu, což obvykle indikuje přítomnost jednotkového kořenu.
Podíváme se, jak by dopadl odhad parametrů IMA procesu, a na základě odhadnutého modelu budeme predikovat 10 budoucích hodnot procesu.
> ima.fit <- arimaCima. sim, order = c(0, 1, 2))
> print(ima.fit)
Call:
arima(x = ima.sim, order = c(0, 1, 2))
Coefficients:
10
Stochastické modely časových řad (8. cvičení)
mal ma2 0.7251 0.2331 s.e.    0.0534 0.0590
sigma~2 estimated as 1.078:    log likelihood = -437.19,    aic = 880.37
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> PlotPredictARIMAdma. sim, ima.fit, n. ahead = 10)
0 50 100 150 200 250
Obrázek 11: Predikce pro simulovanou řadu IMA procesu (1 — B)Yt = (1 + 0.75-B + 0.2B2)et
Opět se přesvědčíme, že diferencování získáme z původního IMA(1, 2) procesu MA(2) proces. Nejprve si prohlédneme graf diferencované časové řady.
> dif.ima <- diff(ima.siiiO
> parGnfrow = c(l, 1), mar = c(2, 2, 1,  0) + 0.05)
> plot (dif. ima)
Stochastické modely časových řad (8. cvičení)
11
~~i—
100
50
150
200
250
300
Obrázek 12: 1. diference simulované řady IMA procesu (1 — B)Yt = (1 + 0.75-B + 0.2B2)et Dále se jako obvykle podíváme na výběrovou autokorelační a parciální autokorelační funkci
> diíIMAací <- ací(dif.ima, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (diíIMAací, main = "Výberová ACF")
Výberová ACF
n
J_L
J_L
10
15
I
20
Obrázek 13: ACF diferencovaného IMA procesu (1 - B)Yt = (1 + 0.75S + 0.2B2)et
> vlMApací <- pací(dií.ima, lag.max = 20, plot = FALŠE)
> par(mírow = c(l, 1), mar = c(2, 2, 3,  0) + 0.05)
> plot (vlMApací, main = "Výberová PACF")
12
Stochastické modely časových řad (8. cvičení)
Výberová PACF
CD O
O
CM O
CM
O -
5 10 15 20
Obrázek 14: PACF diferencovaného IMA procesu (1 - B)Yt = (1 + 0.755 + 0.2B2)et
Z vlastností MA(2) procesu vyplývá, že autokorelační funkce by měla být významně nenulová pouze v bodech k = 0,1,2 a parciální autokorelační funkce by měla exponenciálně klesat k nule. Ani jeden z předchozích dvou grafů není v rozporu s očekáváním.
Na závěr ještě ověříme, že když odhadneme MA(2) model pro diferencovanou řadu, dostaneme podobné odhady parametrů 0\ a 62, jako když jsme odhadovali parametry IMA(1, 2) modelu pro původní nediferencovaná data.
> dif. ima.fit <- arimaCdif. ima., order = c(0, 0, 2))
> print(dif.ima. fit)
Call:
arima(x = dif.ima, order = c(0, 0, 2))
Coefficients:
mal        ma2 intercept 0.7243   0.2324 0.0756 s.e.    0.0534   0.0591 0.1170
sigma~2 estimated as 1.076:    log likelihood = -436.98,    aic = 881.96
příklad 3 - ARIMA proces
Ještě se podíváme na konkrétní případ ARIMA procesu řádu (2,1,2) (tzn. proces, u nějž po jednom diferencování dostaneme ARMA(2,2) proces). Zvolíme konkrérní hodnoty parametrů tpi = 0.9, (f2 = —0.25, 61 = 0.75 a 62 = 0.2, tj. proces lze zapsat ve tvaru
(1 -B)(l- 0.9B + 0.25B2)Yt = (1 + 0.75S + 0.2B2)et (1 - 0.9B + 0.25B2 -B + 0.9B2 - 0.25B3)Yt = (1 + 0.75B + 0.2B2)et (1 - 1.95 + 1.1552 - 0.25B3)yt = (1 + 0.755 + 0.252)£í Yt - l.9Yt-i + 1.15yt_2 - 0.25yt_3 = et + 0.75£í_i + 0.2£í_2
Yt = l.9Yt-i ~ 1.15yt_2 + 0.25yt_3 + et + 0.75£í_i + 0.2£í_2
Stochastické modely časových řad (8. cvičení)
13
Nejprve ověříme kauzalitu a invertibilitu AR a MA části procesu, tj. zjistíme, zda kořeny polynomu $>(z) a Q (z) leží vně jednotkové kružnice.
> ARphi <- c(0.9, -0.25)
> Mod(polyroot(c(l, -ARphi)))
[1] 2 2
> MAthetEL <- c(0.75, 0.2)
> Mod(polyroot(c(l, MAtheta)))
[1] 2.236068 2.236068
Všechny kořeny leží vně jednotkové kružnice, takže se nám podařilo ověřit kauzalitu resp. invertibilitu AR a MA části ARIMA procesu. Teď už můžeme vygenerovat 300 hodnot procesu a získanou časovou řadu znázornit graficky.
> arimaSim <- arima. sim (model = list(order = c (2, 1, 2),
ar = ARphi, ma = MAtheta), n = 300)
> par (mírow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> plot (arimaSim)
o
CD
0 50 100 150 200 250 300
Obrázek 15: Simulovaný ARIMA proces (1 -B)(l -0.9B+ 0.25B2)Yt = (l+0.75B + 0.2B2)et Opět si prohlédneme výběrovou autokorelační a parciální autokorelační funkci.
> vARIMAacf <- acf(arimaSim, lag.max = 20, plot = FALŠE)
> par(mírow = c(l, 1), mar = c(2, 2, 3,  0) + 0.05)
> plot (vARIMAacf, main = "Výběrová ACF")
14
Stochastické modely časových řad (8. cvičení)
Výberová ACF
~i
10
~i
15
I
20
Obrázek 16: ACF simulovaného ARIMA procesu (1 - B)(l - 0.95 + 0.25B2)Yt 0.2B2)st
> vARIMApacf <- pacf(arimaSim, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vARIMApacf, main = "Výberová PACF")
Výberová PACF
[1 + 0.75B +
J_L
—1—l—r
10
i
15
I
20
Obrázek 17: PACF simulovaného ARIMA procesu (l-JB)(l-0.9JB + 0.25JB2)yí = (1+0.755 + 0.2B2)st
Stejně jako v případě předchozích procesů s jednotkovým kořenem vidíme, že autokorelační funkce klesá k nule velmi pozvolna.
Pomocí funkce arima opět odhadneme parametry ARIMA(2,1,2) procesu, abychom je mohli porovnat s teoretickými hodnotami zadanými při simulaci časové řady. Potom na základě odhadnutého modelu provedeme predikci 10 budoucích hodnot procesu.
> arima.fit <- arimaCarimaSim, order = c(2, 1, 2))
> priut(arima.fit)
Call:
arima(x = arimaSim, order = c(2,  1, 2))
Stochastické modely časových řad (8. cvičení)
15
Coefficients:
arl ar 2        mal ma2
0.7284   -0.1648   0.8108 0.2638 s.e.    0.2767     0.2095   0.2814 0.1973
sigma~2 estimated as 0.9598:    log likelihood = -420.72,    aic = 851.44
> par(mfrow = c(l, 1), mar = c(2, 2, 1,  0) + 0.05)
> PlotPredictARIMA(arimaSim, arima.fit, n.ahead = 10)
60 -
0 50 100 150 200 250
Obrázek 18: Predikované hodnoty ARIMA procesu (1 - B)(l - 0.9B + 0.25B2)Yt = (1 + 0.755 + 0.2B2)st
V závěrečné části příkladu ověříme, že diferencováním ARIMA(2,1, 2) procesu dostaneme ARMA(2,2) proces s týmiž parametry. Nejprve proto spočítáme 1. diference simulované časové řady a vykreslíme je do grafu.
> dif.arima <- diff(arimaSiiiO
> par (mfrow = c(l, 1), mar = c (2, 2, 1,  0) + 0.05)
> plot (dif.arima)
16
Stochastické modely časových řad (8. cvičení)
CD I
O 50 100 150 200 250
Obrázek 19: 1. diference simulovaného ARIMA procesu (1 - 5)(1 - 0.95 + 0.25B2)Yt = (1 + 0.755 + 0.2B2)et
Podíváme se na výběrovou autokorelační a parciální autokorelační funkci diferencované řady. Jestliže se jedná o ARMA proces, jak předpokládáme, měly by obě funkce exponenciálně klesat k nule.
> diíARIMAací <- ací (dií.arima, lag.max = 20, plot = FALŠE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (diíARIMAací, main = "Výberová ACF")
Výberová ACF
I    I    I I
~i
10
~i
15
I     I I
I
20
Obrázek 20: ACF simulovaného ARIMA procesu (1 - B)(l - 0.95 + 0.2552)yt = (1 + 0.755 + 0.252)£í
Stochastické modely časových řad (8. cvičení)
17
> vARIMApací <- pací(dií.arima, lag.max = 20, plot = FALSE)
> par (mírow = c(l, 1), mar = c (2, 2, 3,  0) + 0.05)
> plot (vARIMApací, main = "Výberová PACF")
Výberová PACF
_l_L
10
I
15
I
20
Obrázek 21: PACF simulovaného ARIMA procesu (1- B)(l-0.9B + 0.25B2)Yt = (1+0.755 + 0.2B2)ct
Vidíme, že průběhy obou funkcí odpovídají ARMA procesu. Na samotný závěr odhadneme parametry ARMA(2,2) procesu pro diferencovanou řadu, abychom se přesvědčili, že při odhadu parametrů pro původní časovou řadu (ARIMA proces) před diferencováním i při odhadu parametrů ARMA procesu pro časovou řadu po diferencování dostaneme podobné výsledky, které odpovídají teoretickým hodnotám použitým při simulaci procesu.
> dií.arima.fit <- arima(dif.arima, order = c(2, 0, 2))
> print(dif.arima.fitJ
Call:
arima(x = dif.arima, order = c(2, 0, 2))
Coefficients:
arl ar2        mal        ma2 intercept
0.7292   -0.1656   0.8098   0.2632 -0.0433 s.e.    0.2768     0.2096   0.2815   0.1974 0.2673
sigma~2 estimated as 0.9597: log likelihood = -420.71,    aic = 853.41
Příklad 4 - Produkce piva v Austrálii
Nyní vyzkoušíma aplikovat ARIMA modely na reálná data. K tomu použijeme měsíční údaje o množství piva vyrobeného v Austrálii, které byly zjišťovány od roku 1958. Data jsou uložena v souboru beer.txt. Data načteme pomocí funkce read.table jako obvykle a pak z nich vytvoříme časovou řadu, kterou vykreslíme do grafu.
> íileDat <- paste(data.library,   "beer.txt", sep = "")
> beer <- read.table(ÍileDat, header = TRUE)
> beerTS <- ts(beer, start = 1958, frequency = 12)
> par(mar = c(2, 2, 3, 0) + 0.05)
> plot (beerTS, main = "Produkce piva v Austrálii")
18
Stochastické modely časových řad (8. cvičení) Produkce piva v Austrálii
1960       1965       1970       1975       1980       1985 1990
Obrázek 22: Měsíční údaje o množství vyrobeného piva v Austrálii (1958-1991)
V časové řadě lze vysledovat rostoucí trend ve výrobě piva ve zkoumaném období. Z toho důvodu by mohl být vhodným modelem jednoduchý integrovaný proces klouzavých součtů IMA(1,1) (tj. ARIMA{0,1,1)), protože umožňuje modelovat tento trend tak, že přenáší pozorování v čase t—1 do pozorování v čase t a přidává k němu náhodnou odchylku. IMA(1,1) proces bývá pro podobná data často vhodný z toho důvodu, že představuje lineární trend s přidaným bílým šumem.
Parametry IMA procesu odhadneme pomocí funkce arima.
> ima.beer <- arima(beerTS, order = c(0, 1, 1))
> ima.beer
Call:
arima(x = beerTS, order = c(0, 1, 1))
Coefficients: mal -0.3334 s.e. 0.0558
sigma~2 estimated as 360.4:    log likelihood = -1723.27,    aic = 3450.53
Odhadnutý model má podobu
yt = yt_i+et- 0.333et_i.
Odhadnutý model můžeme použít k predikci budoucího vývoje. K tomu v R použijeme jako obvykle funkci predict() s parametrem n.ahead udávajícím počet predikovaných hodnot.
Stochastické modely časových řad (8. cvičení)
19
Pokud by nás například na konci roku 1991 zajímalo, jaké je predikované množství vyrobeného piva za celý následující rok 1992, a to dohromady za všech dvanáct měsíců, spočítáme predikci 12 budoucích hodnot a sečteme je.
> beerl992 <- predictdma.beer, n.ahead = 12)
> sum(beerl992$pred)
[1] 2365.412
Poté, co jsme odhadli model pro časovou řadu, bychom měli správně ověřit vhodnost tohoto modelu. Diagnostickým grafem může být například graf výběrové autokorelační funkce pro rezidua.
> acf(resid(ima.beer))
Series resid(ima.beer)
ACF .0       0.4 0.8 ......					__i__
		■ i	r r		
o o -					
	................................................ r      f í			...................	
0.0 0.5 1.0 1.5 2.0
Lag
Obrázek 23: ACF reziduí IMA(1,1) modelu pro data Měsíční údaje o množství vyrobeného piva v Austrálii (1958-1991)
Graf autokorelační funkce neklesá exponenciálně k nule, ale má vrcholy s periodou 12, což naznačuje, že do modelu by bylo třeba zahrnout ještě sezónní složku (k tomu jsou vhodné například SARIMA modely).
Příklad 5 - Analýza amerického HNP
V posledním příkladu použijeme ARIMA modely k analýze hrubého národního produktu USA. Data uložená v souboru gnp96.dat jsou čtvrtletní údaje o velikosti HNP USA od 1. čtvrtletí roku 1947 do 3. čtvrtletí roku 2002. Pozorování je celkem 223. Konkrétně se jedná o hodnoty reálného HNP v bilionech dolarů, v cenách odpovídajících roku 1996, údaje jsou sezónně očištěné.
Data musíme nejdřív načíst. Soubor je bez hlavičky a obsahuje dva sloupce, v prvním sloupci je uvedeno období a ve druhém odpovídající hodnota HNP. Data načteme funkcí read.table, vytvoříme z nich časovou řadu a vykreslíme je do grafu.
> íileDat <- paste (data. library,   "gnp96 .dat", sep = "")
> gnp96 <- read. table (ÍileDat)
> gnpTS <- ts(gnp96[, 2], start = 1947, frequency = 4)
> plot(gnpTS)
20
Stochastické modely časových řad (8. cvičení)
Obrázek 24: HNP USA v letech 1947-2002
Na první pohled je zřejmé, že v datech je obsažen výrazný trend, který ale přehluší většinu dalších charakteristik časové řady. Z obrázku například není vůbec patrné, že s časem roste také variabilita dat. Podívejme se na autokorelační funkci.
> ELCÍ(gnpTS, lELg.max = 50)
Series gnpTS
o <
n—
6
Lag
I
10
12
Obrázek 25: Autokorelační funkce HNP USA v letech 1947-2002 Pokusme se odstranit trend pomocí diferencování.
> gnp.dif <- diff(gnpTS)
> plot (gnp.dif)
Stochastické modely časových řad (8. cvičení)
21
O LO
O O
O
ľ!= LO "D
Q.
C
O)
o
LO
I
O
o
T
1950       1960       1970       1980       1990 2000 Time
Obrázek 26: 1. diference HNP USA v letech 1947-2002
Z tohoto obrázku je již patrné, že variabilita druhé poloviny dat je vyšší. Navíc se zdá, že ani diferencováním se nám nepodařilo z řady úplně eliminovat trend.
Proto zvolíme poněkud odlišný postup, který je v ekonometrii poměrně obvyklý. Zkoumanou časovou řadu, HNP USA, zlogaritmujeme a pro zlogaritmovaná data spočítáme 1. diference. Budeme tedy pracovat s časovou řadou
Zt = A lny.
Takto vzniklá časová řada má zároveň dobrou interpretaci, protože jednotlivé údaje vyjadřují tempo růstu HNP vztahující se k danému období.
Podívejme se na obrázek nově vytvořené časové řady.
> gnpgr <- diíí(log(gnpTS))
> plot (gnpgr)
22
Stochastické modely časových řad (8. cvičení)
"T
1950
1960
1970 1980 Time
1990
2000
Obrázek 27: Tempa růstu HNP USA v letech 1947-2002
Časová řada se od pohledu jeví jako stacionární. Tomu by měla odpovídat také autokore-lační a parciální autokorelační funkce, kdy by obě měly poměrně rychle klesat k nule.
> par (mírow = c(2, 1), mar = c(2, 2, 3,  0) + 0.05)
> ací(gnpgr, lag.max = 24)
> pací(gnpgr, lag.max = 24, main = "")
Series gnpgr
II
i 1 r
Stochastické modely časových řad (8. cvičení)
23
Obrázek 28: ACF a PACF pro tempa růstu HNP USA v letech 1947-2002
Z výše uvedených grafů ACF a PACF bychom mohli usoudit, že ACF se pro k > 2 významně neliší od nuly, zatímco PACF klesá exponenciálně k nule. To by nás vedlo k názoru, že tempa růstu HNP se dají modelovat pomocí MA{2) procesu, a tudíž zlogaritmované hodnoty HNP můžeme modelovat jako IMA(1,2) proces.
Grafy ACF a PACF ale můžeme zároveň interpretovat tak, že ACF klesá exponenciálně k nule a PACF se neliší významně od nuly pro k > 1. To by vedlo k modelování temp růstu pomocí AR{1) procesu a zlogaritmovaného HNP pomocí ARIMA(1,1,0) procesu.
Postupně odhadneme oba dva modely.
> mEL.gnpgr <- ar ima.(gnpgr, order = c(0,  0, 2))
> mEL.gnpgr
Call:
arima(x = gnpgr, order = c(0, 0, 2))
Coefficients:
mal        ma2 intercept 0.3028   0.2035 0.0083 s.e.    0.0654   0.0644 0.0010
sigma~2 estimated as 8.919e-05:    log likelihood = 719.96,    aic = -1431.93
> ar. gnpgr <- arima(gnpgr, order = c(l,  0, 0))
> &r-gnpgr
Call:
arima(x = gnpgr, order = c(l, 0, 0))
Coefficients:
arl intercept 0.3467 0.0083 s.e.    0.0627 0.0010
sigma~2 estimated as 9.03e-05:    log likelihood = 718.61,    aic = -1431.22
Odhadnuté modely jsou:
MA{2) :    Zt = 0.0083 + et + 0.3028et_i + 0.2035^-2 AR(1) :    Zt = 0.0083 + 0.3467Zt + et
Funkce arima automaticky odhaduje ARIMA model se střední hodnotou, což se projeví v odhadnuté konstantě označené intercept. V tomto konkrétním případě je takový postup zcela na místě, protože průměrný růst HNP USA je kladný. Pokud bychom konstantu do modelu nezahrnuli, znamenalo by to, že předpokládáme průměrný růst HNP rovný 0.
Odhadli jsme dva modely. Pomocí analýzy reziduí, kterou se budeme podrobně zabývat v jednom z pozdějších cvičení, můžeme ověřit, zda jsou oba modely vyhovující. Nyní však
24
Stochastické modely časových řad (8. cvičení)
předpokládejme, že oba modely dobře popisují data. Naskýtá se otázka, který z modelů je lepší, vhodnější?
Ve skutečnosti jsou oba modely velmi podobné. Dá se to snadno ověřit, pokud vyjádříme odhadnutý AR(1) proces pomocí ekvivalentní MA reprezentace (což je, jak víme, u kauzálních AR procesů možné). Budeme při tom uvažovat odhadnutý AR(1) model, v němž budeme prozatím ignorovat odhadnutou konstantu, a vyjádříme ho díky vlastnosti kauzality jako proces MA{<x>).
oo j=0
Prvních deset koeficientů ipj MA reprezentace kauzálního AR{1) procesu zjistíme příkazem ARMAtoMAO.
> ARMAtoMA(ar = 0.35, ma = 0, 10)
[1] 3.500000e-01 1.225000e-01 4.287500e-02 1.500625e-02 5.252187e-03 [6]  1.838266e-03 6.433930e-04 2.251875e-04 7.881564e-05 2.758547e-05
Odtud vidíme, že
yt «et + 0.35et_i +0.12et_2,
což je velmi blízké odhadnutému MA(2) procesu.
Na závěr na základě odhadnutého MA(2) modelu spočítáme predikci pro tempa růstu amerického HNP do roku 2010.
> PlotPredictARIMA(gnpgr, ma.gnpgr, n.ahead = 30)
.04 -
.02 -
1950       1960       1970       1980       1990       2000 2010
Obrázek 29: Predikce tempa růstu amerického HNP
Stochastické modely časových řad (8. cvičení)
25
2 Ukol:
1. Simulujte 300 hodnot IMA(2,2) procesu s hodnotami parametrů 0\ = 0.2, 62 = —0.75
(a) Vykreslete simulovanou časovou řadu.
(b) Vykreslete ACF a PACF.
(c) Odhadněte parametry IMA(2, 2) procesu ze simulovaných dat.
(d) Vykreslete predikci pro 10 budoucích hodnot.
(e) Spočítejte 1. diference simulované řady a prohlédněte si autokorelační funkci.
(f) Spočítejte 2. diference a opět se podívejte na autokorelační funkci.
2. Časová řada v datovém souboru robot.dat udává x-ové souřadnice odpovídající pozicím, v nichž se nacházel průmyslový robot po ukončení série operací. Data udávají odchylky od cílové pozice.
(a) Data načtěte a vykreslete.
(b) Prohlédněte si graf 1. diferencí.
(c) Pro původní (nediferencovaná) data odhadněte parametry modelu IMA(1,1) a vykreslete graf s predikcí pro 10 budoucích hodnot.