ARIMA modely, časť 1
Beáta Stehlíková, FMFI UK Bratislava, Časové rady
- 1 Príklad: Cestovanie lietadlom
- 1.1 Zobrazenie dát, transformácia
- 1.2 Testovanie autokorelácií I.
- 1.3 Testovanie autokorelácií II.
- 1.4 MA(1) proces - definícia a vlastnosti
- 1.5 Knižnica astsa, funkcia sarima
- 1.6 MA(1) proces - odhadovanie v R-ku
- 1.7 Predikcie
- 1.8 MA(q) proces vyššieho rádu
- 1.9 Výpočet ACF zadaného procesu v R-ku
- 1.10 Autoregresné AR(p) procesy
- 1.11 AR(p) proces - odhadovanie v R-ku
- 2 Stacionarita AR procesov
- 3 Stredná hodnota AR procesu
- 4 Cvičenia
- 4.1 Dáta Svetovej banky a grafy pomocou balíka ggplot2
- 4.2 Stacionarita a parameter procesu I.
- 4.3 Simulačné zistenie stacionarity pre AR(2) proces
- 4.4 Stacionarita a parameter procesu II.
- 4.5 Hľadanie modelu pre dáta - úvod
- 4.6 Hľadanie modelu pre dáta I. - ceny benzínu
- 4.7 Hľadanie modelu pre dáta II. - indikátor od Boxa a Jenkinsa
- 4.8 Nájdenie príkladu procesu s danou vlastnosťou
- 5 Nepovinné dodatky
1 Príklad: Cestovanie lietadlom
1.1 Zobrazenie dát, transformácia
Budeme pracovať s dátami z knižnice datasets, ktoré sú uložené v premennej airmiles. Podľa popisu v helpe: The revenue passenger miles flown by commercial airlines in the United States for each year from 1937 to 1960.
library(datasets)
airmiles## Time Series:
## Start = 1937
## End = 1960
## Frequency = 1
## [1] 412 480 683 1052 1385 1418 1634 2178 3362 5948 6109
## [12] 5981 6753 8003 10566 12528 14760 16769 19819 22362 25340 25343
## [23] 29269 30514plot(airmiles)
Podobne ako v prípade počtu cestujúcich z minulého týždňa spravíme logaritmickú transformáciu:
plot(log(airmiles))
Budeme robiť aj predikcie modelov, preto z našich dát vynecháme posledné štyri pozorovania, aby sme mohli predikcie porovnať s reálne pozorovanými hodnotami. Použijeme na to funkciu window:
# DOPLNTE:
y <- window(log(airmiles), end = ...)Minulý týždeň sme hovorili o stacionarite. Modely, ktorými sa budeme zaoberať, budú modely pre stacionárne dáta, preto ešte tieto dáta musíme upravi - zoberieme diferencie:
plot(diff(y))
1.2 Testovanie autokorelácií I.
Cvičenie: Zistite, či diferencie premennej y môžeme považovať za biely šum. Zobrazte odhadnutú autokorelačnú funkciu a P hodnoty Ljung-Boxovho testu.
1.3 Testovanie autokorelácií II.
Môže byť prirodzené čakať, že ak nám teda Ljung-Boxov test biely šum zamietol, významná korelácia bude zodpovedať lagu 1 (teda korelácia medzi rastom v tomto roku a v predchádzajúcom.)
Pripomeňme si zo slajdov (str. 25):
Cvičenie: Chceme zistiť, či pre náš proces platí \(\rho(k)=0\) pre \(k \geq 2\), pričom prvá autokorelácia môže byť nenulová. Vypočítajte intervaly spoľahlivosti a testujte tieto hypotézy.
1.4 MA(1) proces - definícia a vlastnosti
Proces s vlastnosťou, že len prvá autokorelácia je nenulová, sme už videli (str.17):
Táto vlastnosť sa zachová aj vtedy, keď pridáme konštantný člen a hodnotu \(u_{t-1}\) zoberieme s iným (nie jednotkovým) koeficientom. Proces \[x_t = \mu + u_t - \beta u_{t-1}\] sa nazýva MA(1) proces, t. j. moving average proces (proces kĺzavých priemerov) rádu 1.
Cvičenie: Odvoďte jeho autokorelačnú funkciu.
1.5 Knižnica astsa, funkcia sarima
Nainštalujte a načítajte si knižnicu astsa - applied statistical time series analysis
library(astsa)Na odhadovanie našich modelov budeme používať funkciu sarima:
szodpovedá tomu, že budeme odhadovať aj všeobecnejšie modely ako ARIMA, a to také, v ktorých je aj sezónnosť (preto SARIMA)arvyjadruje počet autoregresných členov - budeme sa nimi zaoberať o chvíľuivyjadruje to, koľkokrát sme proces diferencovali, aby sme získali stacionárny časový rad (ide vtedy o tzv. integrovaný proces)mavyjadruje počet moving average členov
Terminológia ohľadom ARIMA modelov: ARIMA(p,k,q) model je taký, pri ktorom sa dáta najskôr k-krát diferencujú a tieto diferencie sa modelujú použitím p autoregresných a q moving average členov.
Syntax funkcie sarima (v prípade modelov bez sezónnosti):
sarima(data, p, k, q) # ARIMA(p,k,q) model pre data v premennej `data`resp. sarima(data, p, k, q, details = FALSE), ak nechceme vypisovať priebeh konvergencie odhadov a kresliť grafy s analýzou rezíduí. To sa môže hodiť napríklad vtedy, ak chceme uložiť model do nejakej premmenej a neskôr s ním pracovať.
1.6 MA(1) proces - odhadovanie v R-ku
Pripomeňme si, že pre prvé diferencie premennej y chceme odhadnúť MA(1) model. To znamená, že funkciu sarima budeme volať nasledovne:
sarima(y, 0, 1, 1) ## initial value -1.901068
## iter 2 value -1.996142
## iter 3 value -2.024002
## iter 4 value -2.031118
## iter 5 value -2.042362
## iter 6 value -2.063857
## iter 7 value -2.072419
## iter 8 value -2.074568
## iter 9 value -2.075498
## iter 10 value -2.075586
## iter 11 value -2.075595
## iter 12 value -2.075595
## iter 13 value -2.075596
## iter 14 value -2.075596
## iter 14 value -2.075596
## iter 14 value -2.075596
## final value -2.075596
## converged
## initial value -1.981058
## iter 2 value -2.011405
## iter 3 value -2.028875
## iter 4 value -2.029789
## iter 5 value -2.030228
## iter 6 value -2.031404
## iter 7 value -2.031543
## iter 8 value -2.031773
## iter 9 value -2.032016
## iter 10 value -2.032103
## iter 11 value -2.032127
## iter 11 value -2.032127
## final value -2.032127
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ma1 constant
## 0.6984 0.2034
## s.e. 0.4263 0.0494
##
## sigma^2 estimated as 0.01658: log likelihood = 11.65, aic = -17.3
##
## $degrees_of_freedom
## [1] 17
##
## $ttable
## Estimate SE t.value p.value
## ma1 0.6984 0.4263 1.6383 0.1197
## constant 0.2034 0.0494 4.1212 0.0007
##
## $AIC
## [1] -2.899462
##
## $AICc
## [1] -2.724462
##
## $BIC
## [1] -3.799889Odhadnutý model zapíšeme pomocou hodnôt constant a ma1 \[\Delta y_t = 0.2034 + u_t + 0.6984 u_{t-1} \] Pozrime sa ešte na výstup. Okrem odhadnutých koeficientov máme informácie o rezíduách na obrázku:
- priebeh
- ACF
- testovanie normality
- výsledok Ljung-Boxovho testu, spolu s hranicou 0,05 pre P hodnoty
Všimnime si, že Ljung-Boxov test začína od lagu 2. To súvisí s poznámkou zo sladov z predchádzajúceho týždňa:
Konkrétne, počet stupňov voľnosti sa znáži o toľko, koľko je v modeli spolu AR a MA koeficientov.
Cvičenie: Vypočítajte (priamo dosadením do vzorca zo slajdov) P hodnotu, ktorá je zobrazená na obrázku, pre lagy 2 (prvý, ktorý sa dá uvažovať) a 5 (zdá sa, že je blízko hranice 0,05). Ak si do premennej model.1 uložíme odhadovaný model model.1 <- sarima(y, 2, 1, 1, details = FALSE), tak pomocou str(model.1) si môžeme pozrieť, aká je štruktúra tohto objektu a ako pristupovať k jeho zložkám. Napríklad rezíduá, ktoré tu potrebujeme, dostaneme ako model.1$fit$residuals.
1.7 Predikcie
Funkcia sarima.for (for ako forecast):
sarima.for(data, n, p, k, q) # n predikovanych hodnot pre `data` podla ARIMA(p,k,q) modeluCvičenie: Predikujte ďalšie pozorovania pre náš časový rad a porovnajte ich so skutočnými hodnotami
1.8 MA(q) proces vyššieho rádu
MA(1) proces sa dá zovšeobecniť prirodzeným spôsobom: Proces \[x_t = \mu + u_t - \beta_1 u_{t-1} - \beta_2 u_{t-2} + \dots + \beta_q u_{t-q}\] sa nazýva MA(q) proces, t. j. moving average proces (proces kĺzavých priemerov) rádu \(q\).
Cvičenie:
- Dokážte, že pre ľubovoľné hodnoty parametrov je proces MA(q) stacionárny.
- Ktoré autokorelácie MA(q) procesu sú nulové?
1.9 Výpočet ACF zadaného procesu v R-ku
ACF procesu sa dá vypočítať pomocou funkcie ARMAacf.
Príklad: Vypočítajte prvých päť hodnôt autokorelačnej funkcie MA(2) procesu \[x_t = 0.2 + u_t - 0.3 u_{t-1} + 0.4 u_{t-2}.\]
Riešenie:
ARMAacf(ma = c(-0.3, 0.4), lag.max = 5)## 0 1 2 3 4 5
## 1.000 -0.336 0.320 0.000 0.000 0.000# ARMAacf(ma = c(-0.3, 0.4), lag.max = 5)[-1] - ak chceme vynechat lag 0Užitočným cvičením je získať tieto hodnoty ručne a vypočítať aj strednú hodnotu a disperziu tohto procesu.
1.10 Autoregresné AR(p) procesy
Ak sa vrátime k dátam o raste precestovaných míľ, mohli by sme chcieť skúsiť tento rast modelovať nie ako kombináciu hodnôt bieleho šumu, ale pomocou jeho starších hodnôt. Teda zistiť, či by sa nedala hodnota v čase \(t\) vysvetliť pomocou hodnoty v čase \(t-1\), od ktorej by závisela. Prípadne pridáme hodnotu v čase \(t-2\), alebo aj staršie.
AR(p) proces je definovaný ako \[x_t = \delta + \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \dots + \alpha_p x_{t-p} + u_t\]
1.11 AR(p) proces - odhadovanie v R-ku
Pripomeňme si, že modelujeme prvé diferencie premennej y. Sksme odhadnúť AR(1) model. To znamená, že funkciu sarima budeme volať nasledovne:
sarima(y, 1, 1, 0) ## initial value -1.877941
## iter 2 value -1.934387
## iter 3 value -1.936101
## iter 4 value -1.936681
## iter 5 value -1.936689
## iter 5 value -1.936689
## iter 5 value -1.936689
## final value -1.936689
## converged
## initial value -1.956384
## iter 2 value -1.956770
## iter 3 value -1.956909
## iter 4 value -1.956909
## iter 4 value -1.956909
## iter 4 value -1.956909
## final value -1.956909
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ar1 constant
## 0.3210 0.2067
## s.e. 0.2126 0.0466
##
## sigma^2 estimated as 0.01985: log likelihood = 10.22, aic = -14.44
##
## $degrees_of_freedom
## [1] 17
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.3210 0.2126 1.5097 0.1495
## constant 0.2067 0.0466 4.4390 0.0004
##
## $AIC
## [1] -2.719541
##
## $AICc
## [1] -2.544541
##
## $BIC
## [1] -3.619967Model síce nie je dobrý - vysvetlite prečo - ale použijeme ho na to, aby sme si ukázali, ako zapísať odhadnutý model pomocou hodnôt constant a ar1:
\[\Delta y_t = 0.2067 + 0.3210 \Delta y_{t-1} + u_t\]
Cvičenie: Vyskúšajte, či sa pridaním ďalších AR členov dá získať model s dobrými rezíduami.
2 Stacionarita AR procesov
2.1 Poznámky
V prípade MA procesov nebol so stacionaritou žiadny problém, MA proces je vždy stacionárny. V prípade AR procesov to tak nie je. Presne to spravíme na prednáške, tu uvedieme niekoľko poznámok:
- Pripomeňme si príklad procesu z minulého týždňa, ktorý nebol stacionárny kvôli rastúcej disperzii:
Ďalšia úvaha tiež súvisí s disperziou. Predpokladajme, že AR(1) proces \[x_t = \delta + \alpha x_{t-1} + u_t\] je stacionárny. Potom z nezávislosti \(x_{t-1}\) a \(u_t\) (proces v minulosti a biely šum dnes) dostaneme \[D[x_t] = \alpha^2 D[x_{t-1} + \sigma^2 \Rightarrow D[x_t] = \frac{\sigma^2}{1-\alpha^2}\]. Aby bola táto dispezia definovaná a aby nebola záporná, musí platiť \(|\alpha|<1\).
Vieme, že korelácia nadobúda hodnoty medzi -1 a 1, ale niekedy dostaneme takéto zvláštne výsledky
ARMAacf(ar = c(0.8, 0.4, 0.3), lag.max = 10)## 0 1 2 3 4 5 6
## 1.000000 3.407407 4.148148 4.981481 6.666667 8.570370 11.017407
## 7 8 9 10
## 14.242074 18.371733 23.699439 30.580866- Funkciou
arima.simmôžeme simulovať ARIMA procesy. Napríklad:
set.seed(1234)
x1 <- arima.sim(model = list(ar = (c(0.2, 0.4))), n = 200)
plot(x1)
Niekedy ale vznikne problém:
set.seed(1234)
x2 <- arima.sim(model = list(ar = (c(0.8, 0.4))), n = 200)Prestavme si, že nestacionárny model z predchádzajúceho bodu \[x_t = 0.8 x_{t-1} + 0.4 x_{t-2} + u_t\] budeme uvažovať bez náhodnej zložky, teda ako deterministickú postupnosť s nejakými štartovacími hodnotami. Aký bude jej priebeh? Kedy deterministický proces \[x_t = \delta + \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \dots + \alpha_p x_{t-p} \] utečie do nekonečne alebo bude oscilovať, a kedy konverguje ku konečnej limite? Stacionárny proces by sa mal ustáliť okolo svojej strednej hodnoty, príslušný deterministický proces by mal mať konečnú limitu.
Dokážte, že na zistenie stacionarity môžeme ekvivalentne počítať korene polynómu \[1 - \alpha_1 z - \alpha_2 z^2 - \dots \alpha_p z^p = 0.\] Tento zápis podmienky stacionarity (prípadne v premennej
Lnamiestoz) je štandardný v literatúre týkajúcej sa časových radov.
2.2 Výpočet koreňov polynómu v R-ku
Funkcia polyroot. Napríklad \[(z-1)(z-2)(z^2+2) = z^4 - 3 z^3 + 4 z^2 - 6 z + 4\] má korene \(1, 2, \sqrt{2}, -\sqrt{2}\).
V R-ku:
polyroot(c(4, -6, 4, -3, 1))## [1] 1+0.000000i 0+1.414214i 0-1.414214i 2+0.000000iCvičenie: Overte stacionaritu niekoľkých zvolených AR procesov a overte svoj záver simuláciou pomocou funkcie arima.sim, ktorá nestacionárny proces odhalí.
3 Stredná hodnota AR procesu
3.1 Vypočet strednej hodnoty
Predpokladajme, že daný AR proces \[x_t = \delta + \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \dots + \alpha_p x_{t-p} + u_t\] je stacionárny. Je dôležité uvedomiť si, že parameter \(\delta\) nie je stredná hodnota procesu. Odvoďte, čomu sa stredná hodnota rovná.
Cvičenie: Funkcia arima.sim generuje proces s nulovou strednou hodnotou. Predpokadajme však, že chceme simulovať proces \[x_t = 0.3 + 0.4 x_{t-1} + 0.2 x_{t-2} - 0.1 x_{t-3} + u_t.\] Doplňte hodnotu konštatnty K v nasledovnom kóde, pri ktorej dostaneme požadovanú simuláciu:
K <- ...
x <- K + arima.sim(model = list(ar = c(0.4, 0.2, -0.1)), n = 250)3.2 Odhadovanie AR modelu v R-ku
V našom príklade s počtom prelietaných míľ sme odhadovali model pre diferencie. Pripomeňme si, že vo výstupe bol odhadnutý parameter constant.
Teraz však odhadneme model pre pôvodnú premennú, nie pre diferencie:
# simulovane data
set.seed(1234)
x <- 10 + arima.sim(model = list(ar = c(0.5, -0.4)), n = 250)
plot(x)
# odhadnuty AR(2) model
sarima(x, 2, 0, 0)## initial value 0.212996
## iter 2 value 0.051310
## iter 3 value 0.023018
## iter 4 value 0.006656
## iter 5 value 0.006166
## iter 6 value 0.006162
## iter 7 value 0.006162
## iter 8 value 0.006162
## iter 9 value 0.006162
## iter 10 value 0.006162
## iter 11 value 0.006162
## iter 11 value 0.006162
## iter 11 value 0.006162
## final value 0.006162
## converged
## initial value 0.006538
## iter 2 value 0.006529
## iter 3 value 0.006524
## iter 4 value 0.006524
## iter 5 value 0.006524
## iter 6 value 0.006524
## iter 6 value 0.006524
## iter 6 value 0.006524
## final value 0.006524
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 xmean
## 0.6169 -0.4317 10.0222
## s.e. 0.0569 0.0569 0.0781
##
## sigma^2 estimated as 1.011: log likelihood = -356.37, aic = 720.73
##
## $degrees_of_freedom
## [1] 247
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.6169 0.0569 10.8460 0
## ar2 -0.4317 0.0569 -7.5894 0
## xmean 10.0222 0.0781 128.2821 0
##
## $AIC
## [1] 1.034576
##
## $AICc
## [1] 1.043229
##
## $BIC
## [1] 0.07683338Vo výstupe je stredná hodnota procesu xmean. Nie je tam konštanta constant, zodpovedajúca parametru \(\delta\) zo zápisu \[x_t = \delta + \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + u_t\]
Cvičenie: Dopočítajte na základe výstupu odhadnutý parameter \(\delta\). Čomu sa rovná jeho presná hodnota (na základe zadaného simulovaného procesu, z ktorého realizácie sme získali tento odhad)?
3.3 Príklad: Spread úrokových mier
- Zdroj: Mills, Markellos: The Econometric Modelling of Financial Time Series. Cambridge University Press, 2008
- Štvrťročné dáta, 1952Q1 - 2005Q4
- Premenné:
- krátkodobá úroková miera (3 mesiace)
- dlhodobá úroková miera (20 rokov)
- Budeme modelovať rozdiel dlhodobej a krátkodobej úrokovej miery, tzy. spread
Načítajte dáta a vytvorte premennú spread:
rs <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr17/data/RSQ.txt") # short term rate
r20 <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr17/data/R20Q.txt") # 20Y rate
spread <- r20 - rsZ premennej spread spravíme objekt typu time-series, ktorý bude obsahovať informáciu o časovej štruktúre dát.
Použijeme funkciu ts:
spread <- ts(spread, # data su premennej spread
frequency=4, # su kvartalne
start=c(1952,1)) # zacinaju v 1. kvartali roku 1952
class(spread) # budeme vidiet, ze uz je to casovy rad## [1] "ts"Vykreslíme teraz priebeh - na x-ovej osi budú správne časy:
plot(spread)
Cvičenie: Odhadnite AR(1), AR(2) a AR(3) model a zhodnoťte, či sú vyhovujúce (stacionárne a rezíduá sa dajú považovať za biely šum)
4 Cvičenia
4.1 Dáta Svetovej banky a grafy pomocou balíka ggplot2
V tomto príklade budeme vidieť:
- prístup k dátam Svetovej banky priamo v R-ku
- príklad použitia knižnice ggplot2 na zostrojenie pekných grafov
Pripravíme si tak dáta, ktoré potom použijeme na modelovanie AR procesom.
Nainštalujte a načítajte knižnice
WDI- World Development Indicators, prístup k dátam budeme mať priamo v Rggplot2- pekné grafy :)
library(WDI)
library(ggplot2)4.1.1 Vyhľadávanie dát vo WDI
Hľadáme indikátory, ktore maju v nazve gdp:
WDIsearch('gdp') Kto pozná regulárne výrazy (regular expressions) - dajú sa používať aj tie. Ma ukážku - hľadajme indikátory, ktoré hovoria o HDP na obyvateľa - vypíšme si tie, ktoré obsahujú gdp a capita
WDIsearch('gdp.*capita')## indicator
## [1,] "6.0.GDPpc_constant"
## [2,] "SE.XPD.TERT.PC.ZS"
## [3,] "SE.XPD.SECO.PC.ZS"
## [4,] "SE.XPD.PRIM.PC.ZS"
## [5,] "NY.GDP.PCAP.PP.KD.ZG"
## [6,] "NY.GDP.PCAP.PP.KD.87"
## [7,] "NY.GDP.PCAP.PP.KD"
## [8,] "NY.GDP.PCAP.PP.CD"
## [9,] "NY.GDP.PCAP.KN"
## [10,] "NY.GDP.PCAP.KD.ZG"
## [11,] "NY.GDP.PCAP.KD"
## [12,] "NY.GDP.PCAP.CN"
## [13,] "NY.GDP.PCAP.CD"
## [14,] "UIS.XUNIT.GDPCAP.5T8.FSHH"
## [15,] "UIS.XUNIT.GDPCAP.5T8.FSGOV"
## [16,] "UIS.XUNIT.GDPCAP.3.FSGOV"
## [17,] "UIS.XUNIT.GDPCAP.23.FSHH"
## [18,] "UIS.XUNIT.GDPCAP.23.FSGOV"
## [19,] "UIS.XUNIT.GDPCAP.2.FSGOV"
## [20,] "UIS.XUNIT.GDPCAP.1.FSHH"
## [21,] "UIS.XUNIT.GDPCAP.1.FSGOV"
## [22,] "UIS.XUNIT.GDPCAP.02.FSGOV"
## [23,] "FB.DPT.INSU.PC.ZS"
## [24,] "NV.AGR.PCAP.KD.ZG"
## [25,] "NE.GDI.FTOT.CR"
## name
## [1,] "GDP per capita, PPP (constant 2011 international $) "
## [2,] "Government expenditure per student, tertiary (% of GDP per capita)"
## [3,] "Government expenditure per student, secondary (% of GDP per capita)"
## [4,] "Government expenditure per student, primary (% of GDP per capita)"
## [5,] "GDP per capita, PPP annual growth (%)"
## [6,] "GDP per capita, PPP (constant 1987 international $)"
## [7,] "GDP per capita, PPP (constant 2011 international $)"
## [8,] "GDP per capita, PPP (current international $)"
## [9,] "GDP per capita (constant LCU)"
## [10,] "GDP per capita growth (annual %)"
## [11,] "GDP per capita (constant 2010 US$)"
## [12,] "GDP per capita (current LCU)"
## [13,] "GDP per capita (current US$)"
## [14,] "Initial household funding per tertiary student as a percentage of GDP per capita"
## [15,] "Initial government funding per tertiary student as a percentage of GDP per capita"
## [16,] "Initial government funding per upper secondary student as a percentage of GDP per capita"
## [17,] "Initial household funding per secondary student as a percentage of GDP per capita"
## [18,] "Initial government funding per secondary student as a percentage of GDP per capita"
## [19,] "Initial government funding per lower secondary student as a percentage of GDP per capita"
## [20,] "Initial household funding per primary student as a percentage of GDP per capita"
## [21,] "Initial government funding per primary student as a percentage of GDP per capita"
## [22,] "Initial government funding per pre-primary student as a percentage of GDP per capita"
## [23,] "Deposit insurance coverage (% of GDP per capita)"
## [24,] "Real agricultural GDP per capita growth rate (%)"
## [25,] "GDP expenditure on gross fixed capital formation (in IDR Million)"Ak je veľa výsledkov, môžeme chciet vypísať niekoľko prvých:
WDIsearch('gdp.*capita')[1:5,]## indicator
## [1,] "6.0.GDPpc_constant"
## [2,] "SE.XPD.TERT.PC.ZS"
## [3,] "SE.XPD.SECO.PC.ZS"
## [4,] "SE.XPD.PRIM.PC.ZS"
## [5,] "NY.GDP.PCAP.PP.KD.ZG"
## name
## [1,] "GDP per capita, PPP (constant 2011 international $) "
## [2,] "Government expenditure per student, tertiary (% of GDP per capita)"
## [3,] "Government expenditure per student, secondary (% of GDP per capita)"
## [4,] "Government expenditure per student, primary (% of GDP per capita)"
## [5,] "GDP per capita, PPP annual growth (%)"4.1.2 Pomocou funkcie WDI, ukážka použitia:
data.hdp <- WDI(indicator='NY.GNP.PCAP.CD',
country=c('CA','US','FR'),
start=1975)indicatorsme našli pomocouWDIsearchcountryvo formáte iso2cstart- ak zvolíme príliš skorý začiatok a dáta nebudú dostupné, upravímestartalebo v načítaných dátach odstránime tie sNAiso2ckódy štátov sa dajú nájsť napríklad tu: https://datahub.io/core/country-list alebo pomocou knižnicecountrycodev R-ku (https://cran.r-project.org/web/packages/countrycode/)
head(data.hdp)## iso2c country NY.GNP.PCAP.CD year
## 1 CA Canada 44860 2018
## 2 CA Canada 43000 2017
## 3 CA Canada 43950 2016
## 4 CA Canada 47590 2015
## 5 CA Canada 52190 2014
## 6 CA Canada 52810 20134.1.3 Graf pomocou ggplot
ggplot(data.hdp, aes(year, NY.GNP.PCAP.CD, color = country)) +
geom_line() +
xlab('Year') + ylab('GDP per capita') +
labs(title = 'GDP per capita (current USD)')
dataje data frame, v ktorom sú naše premennéyear,NY.GNP.PCAP.CD- z dát uložených vdatabude na x-ovej osi premennáyeara na y-ovejNY.GNP.PCAP.CDcolor = country- grafy budú odlíšené farebne (pretocolor) podľa premennejcountry
4.1.4 AR(1) model pre rýchlosť rastu HDP
Jednoduchý AR(1) model je často dobrým modelom pre rýchlosť rastu HDP. Vyskúšame ho pre nasledovné dáta.
# GDP per capita (constant 2000 US$)
data <- WDI(indicator='NY.GNP.PCAP.KD',
country=c('US'),
start=1965)
# zoradime rastuco podla rokov
data <- data[order(data$year),]
# nasa premenna na modelovanie
log.y <- log(data$NY.GNP.PCAP.KD)
# doplnime casovu strukturu
log.y <- ts(log.y, start=1965, frequency=1)Zadanie:
- Z priebehu dát vidíme, ze nie sú stacionárne, pretovbudeme pracovat s diferenciami
- Zoberieme diferencie logaritmov - to je rýchlost rastu HDP.
- Vyskúšajte, či je AR(1) model dobrým modelom pre rýchlosť rastu HDP.
- Zopakujte pre iný štát.
4.2 Stacionarita a parameter procesu I.
Zadanie: Určte všetky hodnoty parametra \(k\), pre ktoré je stacionárny AR(2) proces \(x_t=x_{t−1}+ k x_{t−2} +u_t\).
Skúsime najskôr numericky: Pre nejaký rozsah parametra k nájdeme absolútnu hodnotu koreňov.
Otázky k nasledujúcemu obrázku:
Preco pre niektoré k je len jedna absolútna hodnota (napríklad k = -1.5 - modrá), kým pre iné sú dve (k = -0.2 - zelená, k = 0.5 - hnedá)?
Ktoré k vyhovovujú podmienke stacionarity a ktoré nie?
Analytický výpočet:
- Odvoďte korene a ich absolútne hodnoty.
- Zistite, kedy je splnená podmienka stacionarity.
- Skontrolujte si, či sa vaše výpočty zhodujú s numerickými na predchádzajúcom obrázku
Precvičenie programovania v R-ku: Spravte samostatne graf z predchádzajúceho obrázku
4.3 Simulačné zistenie stacionarity pre AR(2) proces
Uvažujme proces \[x_t=\delta +\alpha_1 x_{t-1} +\alpha_2 x_{t-2} + u_t\] Má dva paramametre, takže ich môžeme znázorniť do roviny a odlíšiť tie, ktoré zodpovedajú stacionárnemu a tie, ktoré zodpovedajú nestacionárnemu AR(2) procesu.
df <- expand.grid(alpha1 = seq(from = -2.5, to = 2.5, by = 0.1),
alpha2 = seq(from = -2.5, to = 2.5, by = 0.1))
df$stacionarita <- NA
head(df)## alpha1 alpha2 stacionarita
## 1 -2.5 -2.5 NA
## 2 -2.4 -2.5 NA
## 3 -2.3 -2.5 NA
## 4 -2.2 -2.5 NA
## 5 -2.1 -2.5 NA
## 6 -2.0 -2.5 NANapíšte funkciu, ktorá rozhodne, či je AR(2) proces so zadanými parametrami stacionárny. Napríklad stac.ar.2(c(1,1)) vráti FALSE.
stac.ar.2 <- function(koef) # doplnte definiciu funkcie
for (i in 1:nrow(df)) {
# doplnte cyklus
# TRUE, FALSE podla toho, ci je prislusny proces stacionarny
# bez cyklu v nepovinnych dodatkoch
df$stacionarita[i] <-
}
plot(...) # body (alpha1, alpha2) odlisene farebne podla stlpca `stacionarita`4.4 Stacionarita a parameter procesu II.
V jednej z domácich úloh bude úlohou nájsť príklady z kurzov časových radov na zahraničných univerzitách a samostatne ich vyriešiť. Takto sa našlo aj nasledovné zadanie:
Zistite, pre ktoré hodnoty parametra \(k\) je stacionárny proces \[x_t=x_{t−1}+k x_{t−2}+ k x_{t−3}+u_t.\]
Zo vzorového riešenia vyplýva, že malo ísť o proces \[x_t=x_{t−1}+k x_{t−2} - k x_{t−3}+u_t.\]
Vyriešte obidve verzie zadania (preklepom sa úloha stala o dosť zložitejšou, ale podmienka stacionarity sa dá odvodiť analyticky aj pre prvý proces).
4.5 Hľadanie modelu pre dáta - úvod
Nasledujúce dve zadania sú zo starých skúšok. Úlohou je
nájsť vhodný AR alebo MA model pre zadané dáta
vysvetliť, prečo sú rezíduá vyhovujúce - skomentovať ACF rezíduí aj Ljung-Boxov test, pričom treba presne povedať, aká hypotéza sa testuje a či sa v tomto prípade zamieta alebo nie (a prečo sme s tým výsledkom spokojní)
overiť stacionaritu získaného modelu (pre AR model napísať polynóm, ktorého korene overujete, aké absolútne hodnoty koreňov vyšli a prečo sme s tým spokojní)
Na skúške je súčasťou kostry zistenie, či môžeme pracovať priamo so zadanými dátami alebo s diferenciami (trend nie je jediným dôvodom diferencovania - budeme sa tomu venovať neskôr). V týchto zadaniach je povedané, s akými dátami máte pri hľadaní AR alebo MA modelu pracovať.
Takisto je na skúške na výber širšia trieda modelov, ako sú AR a MA modely (AR a MA členy sa dajú spojiť - o tom neskôr). V niektorých prípadoch - vrátane týchto - to stačí.
4.6 Hľadanie modelu pre dáta I. - ceny benzínu
Dáta z knižnice astsa, z popisu v helpe: New York Harbor conventional regular gasoline weekly spot price FOB (in cents per gallon) from 2000 to mid-2010. Zoberieme dáta od roku 2006.
library(astsa)
data(gas)
x <- window(gas, start=c(2006,1))Nájdite model pre premennú x tak, že budete modelovať jej diferencie.
4.7 Hľadanie modelu pre dáta II. - indikátor od Boxa a Jenkinsa
Dáta z knižnice astsa, z popisu v helpe: Leading indicator, 150 months; taken from Box and Jenkins (1970).
library(astsa)
data(lead)
x <- leadNájdite model pre premennú x tak, že budete modelovať jej diferencie.
4.8 Nájdenie príkladu procesu s danou vlastnosťou
Nájdite príklad
- stacionárneho AR(3) procesu,
- nestacionárneho AR(4) procesu,
- stacionárneho AR(2) procesu, ktorého stredná hodnota je rovná 10.
V každom z týchto prípadov dokážte, že váš proces naozaj má požadovanú vlastnosť.
Príklady tohto typu “nájdite príklad procesu, ktorý … a ukážte, že má požadovanú vlastnosť” budú aj na skúške.
5 Nepovinné dodatky
5.1 Viac o grafoch pomocou ggplot
Z fakultnej siete je na stránke vydavateľstva Springer dostupná kniha ggplot2: Elegant Graphics for Data Analysis od Hadleya Wickhama v pdf formáte: https://link.springer.com/book/10.1007/978-3-319-24277-4
5.2 Iné pekné grafy pomocou dygraphs
Ukážka:
x <- Nile
library(dygraphs)
dygraph(x) %>% dySeries(label="Annual flow", color="black") %>% dyRangeSelector()O knižnici dygraphs: https://rstudio.github.io/dygraphs/index.html
5.3 Regulárne výrazy pre pristupovanie k parametrov odhadnutého modelu
Ak odhadneme model, hodnoty koeficientov nemusíme odpisovať alebo kopírovať z výstupu, tak by sme navyše dostali iba ich približné hodnoty. Dajú sa získať presné. Okrem toho, budeme mať aj zložitejšie modely a môže byť praktické mať automatický prístup ku koeficientom jednotlivých typov.
Vygenerujme si dáta a zložitejší model. Nebude nás teraz zaujímať, čo znamenajú jednotlivé vstupné parametre, iba vektor s odhadnutými koeficientami a názvy jeho zložiek.
y <- arima.sim(model = list(ar=c(0.5,-0.4)), n = 250)
model <- sarima(y, 3, 0, 2, 2, 0, 1, 12, details = FALSE)
model## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 sar1 sar2 sma1
## 1.1888 -0.6390 0.2847 -0.6595 -0.2477 0.1612 0.0233 -0.1945
## s.e. 0.1518 0.1738 0.0993 0.1489 0.1119 5.3793 0.1777 5.3805
## xmean
## -0.0169
## s.e. 0.0378
##
## sigma^2 estimated as 1.069: log likelihood = -363.5, aic = 746.99
##
## $degrees_of_freedom
## [1] 241
##
## $ttable
## Estimate SE t.value p.value
## ar1 1.1888 0.1518 7.8301 0.0000
## ar2 -0.6390 0.1738 -3.6766 0.0003
## ar3 0.2847 0.0993 2.8661 0.0045
## ma1 -0.6595 0.1489 -4.4292 0.0000
## ma2 -0.2477 0.1119 -2.2146 0.0277
## sar1 0.1612 5.3793 0.0300 0.9761
## sar2 0.0233 0.1777 0.1311 0.8958
## sma1 -0.1945 5.3805 -0.0361 0.9712
## xmean -0.0169 0.0378 -0.4480 0.6546
##
## $AIC
## [1] 1.138679
##
## $AICc
## [1] 1.150361
##
## $BIC
## [1] 0.2654513Koeficienty:
model$fit$coef## ar1 ar2 ar3 ma1 ma2 sar1
## 1.18876175 -0.63901327 0.28472370 -0.65951300 -0.24771334 0.16122962
## sar2 sma1 xmean
## 0.02330651 -0.19447400 -0.01693443Takže ak definujeme
koef <- model$fit$coeftak máme napríklad
koef[1:3]## ar1 ar2 ar3
## 1.1887618 -0.6390133 0.2847237Povedzme, že z koeficientov chceme vybrať tie, ktorých mená sa začínajú na ma.
Názvy zložiek vektora sú
names(koef)## [1] "ar1" "ar2" "ar3" "ma1" "ma2" "sar1" "sar2" "sma1" "xmean"Ak by sme chceli indexy tých, ktoré obsahujú ma, spravili by sme:
grep("ma", # co sa ma hladat
names(koef)) # kde sa to ma hladat## [1] 4 5 8My ale checeme iba tie, ktoré začínajú ma (lebo je tam aj sma1 koeficient, ktorý necheme). Tu prichádzajú na rad regulárne výrazy:
grep("^ma", # co sa ma hladat
names(koef)) # kde sa to ma hladat## [1] 4 5Príslušné koeficienty potom sú
koef[grep("^ma", names(koef))]## ma1 ma2
## -0.6595130 -0.2477133Tutoriál o regulárnych výrazov v R-ku: https://github.com/daattali/UBC-STAT545/blob/master/reference/regex/regularExpressions.md
5.4 Trieda funkcií apply
Vráťme sa k príkladu, kde sme pre zadané kombinácie parametrov zisťovali, či je AR(2) proces s týmito koeficientami stacionárny.
df <- expand.grid(alpha1 = seq(from = -2.5, to = 2.5, by = 0.1),
alpha2 = seq(from = -2.5, to = 2.5, by = 0.1))
df$stacionarita <- NA
head(df)## alpha1 alpha2 stacionarita
## 1 -2.5 -2.5 NA
## 2 -2.4 -2.5 NA
## 3 -2.3 -2.5 NA
## 4 -2.2 -2.5 NA
## 5 -2.1 -2.5 NA
## 6 -2.0 -2.5 NANamiesto toho, aby sme hodnoty stĺpca stacionarita vkladali v cykle, môžeme použiť funkciu apply.
# Predpokladame fungujucu funkciu stac.ar.2, ktora ma ako vstup dvojzlozkovy vektor a vystup TRUE/FALSE
apply(df[, 1:2], # vstupne data = 1. a 2. stlpec df
MARGIN = 1, # znamena, ze sa budu brat po riadkoch
FUN = stac.ar.2) # aplikuje sa funkcia stac.ar.2Výsledok teda môžeme vložiť do df:
df$stacionarita <- apply(...)