1 Inštalácia knižníc v R-ku

Začneme s inštaláciou knižníc potrebných na toto cvičenie. Spustením nasledovného príkazu sa nainštalujú tie z knižníc, ktoré nainštalované nie sú, a zároveň sa aj načítajú.

# názvy balíkov, ktoré budú potrebné pre cvičenie 1
packages <- c("astsa","priceR","datasets","forecast","WDI","ggplot2")

# inštalácia balíkov, ak nie sú nainštalované
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages],dependencies=TRUE)
}

# načítanie balíkov
invisible(lapply(packages, library, character.only = TRUE))

2 Overovanie stacionarity AR procesu v R-ku

V prípade AR(1) procesu, je overovanie stacionarity zrejmé - stačí sa pozrieť na absolútnu hodnotu autoregresného koeficientu, ktorý musí byť v absolútnej hodnote menší ako 1.

Pre všeobecný AR(p) proces je v prednáškach odvodená podmienka stacionarity: proces \((1-\alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p)x_t = \delta + u_t\) je stacionárny práve vtedy, keď sú všetky korene polynómu \(1-\alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p\) v absolútnej hodnote väčšie ako 1 (geometricky: mimo jednotkového kruhu).

Funkcia v R-ku: Korene polynómu \(1-\alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p\) nájdeme pomocou funkcie polyroot, ktorá je vysvetlená v slajdoch z prenášky.

2.1 Príklad: Overenie stacionarity AR(2) procesu

Zadanie: Zistíme, či je proces \((1−0.9L+0.6L2^)x_t=u_t\) stacionárny.

Výpočet v R-ku: Hľadáme korene (a ich absolútne hodnoty) polynómu \(1−0.9L+0.6L^2\), teda:

polyroot(c(...)) # korene polynómu 
abs(polyroot(c(...))) # absolútne hodnoty koreňov polynómu

## [1] 0.75+1.050793i 0.75-1.050793i

## [1] 1.290994 1.290994

Záver: Korene sú mimo jednotkového kruhu, lebo absoútne hodnota je pre každý koreň väčšia ako 1. Teda náš proces JE STACIONÁRNY.

2.2 Cvičenie: Overenie stacionarity AR procesu (súčasť kostry na skúške)

Zistite, či sú nasledovné AR procesy stacionárne:

\((1 + 0.3L + 0.2L^2)x_t = 5+u_t\)
\((1 - 0.25L + 0.6L^2 - 0.55L^3)x_t = u_t\)
\(x_t = 0.8x_{t-1} + 0.3x_{t-2} + 0.2x_{t-3} + u_t\)
\(x_t = 0.25 + 0.1x_{t-1} + 0.2x_{t-3} + u_t\)

Vypočítajte strednú hodnotu stacionárnych procesov.

2.3 Cvičenie: Grafické znázornenie koreňov

Znázornite korene polynómov z predošlého cvičenia v komplexnej rovine spolu s jednotkovým kruhom.

Užitočné funkcie a parametre:

ak dostne plot ako vstup komplexné čísla, nakreslí na x-ovú os ich reálnu a na y-ovú os imaginárnu časť
parameter asp = 1 zabezpečí rovnakú mierku na obidvoch osiach
funkcia curve nakreslí zadanú krivku pre zadaný rozsah nezávislej premennej, s parameterom add = TRUE ju dokreslí do pôvodného obrázku
osi môžeme kresliť pomocou abline (horizontálnu s parametrom h, vertikálnu s parametrom v)

2.4 Cvičenie: 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á 15

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.

3 Hľadanie AR modelu

V tejto časti budeme modelovať infláciu. Z balíka priceR použijeme príkaz retrieve_inflation_data na získanie dát z World Bank API.

retrieve_inflation_data(country, countries_dataframe)

Budeme skúmať infláciu v USA. Zvoľme teda `country=“US”’.

data <- retrieve_inflation_data("US")

## Validating iso2Code for US 
## Generating URL to request all 296 results
## Retrieving inflation data for US 
## Generating URL to request all 64 results

Vašou prvou úlohou je vytvoriť z týchto dát pomocou príkazu ts(...) časový rad \(x\) s hodnotami inflácie od roku 1960 po rok 2023. Užitočné bude začať so zobrazením štruktúry objektu data.

str(data)

## 'data.frame':    64 obs. of  8 variables:
##  $ indicator      :'data.frame': 64 obs. of  2 variables:
##   ..$ id   : chr  "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" ...
##   ..$ value: chr  "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" ...
##  $ country        :'data.frame': 64 obs. of  2 variables:
##   ..$ id   : chr  "US" "US" "US" "US" ...
##   ..$ value: chr  "United States" "United States" "United States" "United States" ...
##  $ countryiso3code: chr  "USA" "USA" "USA" "USA" ...
##  $ date           : chr  "2023" "2022" "2021" "2020" ...
##  $ value          : num  4.12 8 4.7 1.23 1.81 ...
##  $ unit           : chr  "" "" "" "" ...
##  $ obs_status     : chr  "" "" "" "" ...
##  $ decimal        : int  1 1 1 1 1 1 1 1 1 1 ...

Po vykreslení časového radu \(x\) by ste mali dostať nasledovný priebeh:

plot(x)

Neskôr uvidíme, že takéto dáta treba diferencovať, aj keď v nich nie je trend. Preto teraz budeme pracovať s diferenciami . Pri modelovaní si treba uvedomiť, že konštantný člen v modeli pre diferencie znamená lineárny trend v pôvodných dátach a tento trend sa dostane aj do predikcií.

V dátach o inflácii nie je trend \(\Rightarrow\) chceme odhadnúť model pre diferencie bez konštanty
Vo funkcii sarima pridáme parameter no.constant na TRUE
Ak sa diferencuje viac ako raz, spraví sa to automaticky

# ar(p) model bez konstanty pre 1. diferencie
sarima(data, p, 1, 0, no.constant = TRUE)

3.0.1 Úlohy

Pre diferencie zistite či sú pre ne na základe rezíduí dobrými modelmi AR(p) pre \(p \in \{0, 1, 2,3,4,5 \}\) bez konštanty
Z vyhovujúcich modelov vyberte ten, ktorý má najnižšie BIC.

Zopakujte pre infláciu v štátoch (všade treba diferencovať):

Francúzsko (FR)
Japonsko (JP)
Kanada (CA)

4 Konštrukcia predikcií

Cieľom modelovania je často konštrukcia predikcií do budúcnosti. V R-ku na to existuje funkcia sarima.for.

sarima.for(data, n, p, 0, 0)  # predikcie pre n pozorovani z AR(p) modelu
sarima.for(data, n, p, k, 0)  # predikcie pre n pozorovani premennej data, ak jej k-te diferencie modelujeme ako AR(p)

Vynechajte z dát posledných 10 pozorovaní.

x_skratene <- window(x, start=1960, end=2013)

Zo zostávajúcich dát odhadnite AR model pre diferencie bez konštanty, ktorý má najmenšie BIC. Pomocou funkcie sarima.for spravte predikcie pre nasledujúcich 10 rokov. Mali by ste dostať výstup:

Nakreslite do jedného grafu farebne odlíšené: dáta použité na odhadovanie modelu, predikcie, predikcie +/- 2 ′× štandardná odchýlka, skutočný vývoj v posledných desiatich rokoch. Na to sa dá dobre využiť funkcia ts.plot. Príklad použitia:

 ts.plot(data1, data2, data3,
        gpars=list(xlab="time", ylab="popis dat",
                  col=c("red","blue","black"))
)

Výsledok pre porovnanie:

5 Zautomatizovanie testovania viacerých AR modelov

Uvažujme dáta o prietokoch Nílu z Cvičenia 2.

x <- Nile
plot(x)

Prípravná úloha: Odhadujeme pre dáta AR(2) model. Čomu sa rovná p-hodnota z Ljung-Boxovho testu pre rezíduá, ak testujeme nulovosť prvých piatich autokorelácií?

# POSTUP
model <- sarima(...)    # odhadneme model
model$...               # chceme pristupovat k reziduam
Box.test(..., fitdf = ) # Ljung-Boxov test pre rezidua

Napíšte funkciu, ktorá pre zadané dáta, rád autoregresného procesu a maximálny uvažovaný počet lagov vráti TRUE/FALSE podľa toho, či sú všetky p-hodnoty Ljung-Boxovho testu pre rezíduá väčšie ako 0.05.

ARmodel <- function(data, p, lag.max = 20){ # nastavenie defaultnej hodnoty vstupneho parametra
  
  # DOPLNTE
  
}

# testujeme prvych k autokorelacii (k<=15) rezidui z AR(5) modelu pre data x
ARmodel(data = x, p = 5, lag.max = 15)

Pomocou funkcie sapply zistite, ktoré AR(p) modely s rádom nanajvýš 5 majú dobré rezíduá. Upravte funkciu ARmodel tak, aby vrátila aj hodnotu Bayesovho informačného kritéria. Pomocou tejto funkcie potom spomedzi modelov s dobrými rezíduami vyberte ten, ktorý má najnižšiu hodnotu informačného kritéria.

##                [,1]     [,2]    [,3]     [,4]     [,5]     [,6]    
## good_residuals FALSE    TRUE    TRUE     TRUE     TRUE     TRUE    
## bic            13.18242 12.9372 12.94383 12.97586 13.02168 13.06053

6 Woldova reprezentácia AR procesu

Uvažujte AR(3) proces \[x_t=1.5x_{t−1}−0.8x_{t−2}+0.2x_{t−3}+u_t\].

Ukážte, že je tento proces stacionárny.
Pomocou funkcie ARMAtoMA nájdite prvých 15 koeficinetov Woldovej reprezentácie procesu.

Výstup pre kontrolu:

##  [1] 1.5000000 1.4500000 1.1750000 0.9025000 0.7037500 0.5686250 0.4704375
##  [8] 0.3915062 0.3246344 0.2678341 0.2203448 0.1811769 0.1490563 0.1227119
## [15] 0.1010582

7 Autokorelačná funkcia AR procesu

Nájdeme autokorelačnú funkciu AR(3) procesu \[x_t=1.5x_{t−1}−0.8x_{t−2}+0.2x_{t−3}+u_t\]

Ukážeme si dva spôsoby riešenia:

priamo z Yule-Wolkerových rovníc a diferenčnej rovnice z prednášky,
pomocou funkcie v R-ku.

7.1 Z Yule Wolkerových rovníc

Pri tomto spôsobe riešenia je prvým krokom zostavenie Yule-Wolkerových rovníc. Tie následne vyriešime v R-ku. Dostaneme ACF(k) pre \(k=1,2,3\) (pozrite ?matrix, ?solve). Následne v cykle z diferenčnej rovnice vypočítame ďalšie hodnoty ACF.

rho <- rep(0, times=10)
rho <- ....                   # prve tri zlozky ako riesenie sustavy
for (i in 4:19) rho[i] <- ... # ostatne zlozky z diferencnej rovnice

Výstup pre kontrolu:

##  [1] 0.9178082 0.7602740 0.6061644 0.4845890 0.3940068 0.3245719 0.2685702
##  [8] 0.2219991 0.1830569 0.1507001

7.2 Funkcia v R-ku

Použijeme funkciu ARMAacf pre náš proces.

ARMAacf(ar=c(...), lag.max=10)

Výstup pre kontrolu:

##         0         1         2         3         4         5         6         7 
## 1.0000000 0.9178082 0.7602740 0.6061644 0.4845890 0.3940068 0.3245719 0.2685702 
##         8         9        10 
## 0.2219991 0.1830569 0.1507001

Porovnajte výsledky z oboch postupov. Funkcia ARMAacf by už nemala byť “čiernou skrinkou” :).

8 Parciálna autokorelačná funkcia AR procesu

8.1 Opakovanie z prednášky

Čo je parciálna autokorelačná funkcia procesu?
Čím je charakteristická PACF autoregresného procesu rádu p? Prečo? Ako to vyplýva z definície PACF?

8.2 Výpočet v R-ku

Nájdeme parciálnu autokorelačnú funkciu AR(3) procesu \[x_t=1.5x_{t−1}−0.8x_{t−2}+0.2x_{t−3}+u_t\]

Podobne ako v prípade ACF, aj teraz použijeme funkciu ARMAacf, ale pridáme parameter pacf=TRUE.

ARMAacf(ar=c(...), lag.max=10, pacf=TRUE)

##  [1]  9.178082e-01 -5.208333e-01  2.000000e-01  0.000000e+00  5.033930e-16
##  [6] -1.123254e-15  8.351739e-16  9.367857e-17 -6.465462e-16  5.159779e-16

8.3 Výberová PACF z dát

Môžeme

použiť funkciu pacf,
pomocou funkcie acf2 z knižnice astsa vykresliť naraz výberovú ACF aj PACF (tu bude mať navyše ACF odstránenú hodnotu pre lag 0, ktorá sa vždy rovná 1).

Napríklad:

x <- rnorm(100)
pacf(x)    # moznost 1

acf2(x)    # moznost 2

##       [,1]  [,2] [,3]  [,4]  [,5] [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.04 -0.11 0.09 -0.01 -0.04 0.06 -0.07 0.02  0.00  0.13  0.01  0.04 -0.05
## PACF -0.04 -0.11 0.08 -0.01 -0.02 0.05 -0.07 0.03 -0.02  0.15  0.01  0.07 -0.07
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## ACF  -0.02 -0.02  0.01  0.01  0.03 -0.05  0.12
## PACF -0.02 -0.03  0.00  0.03  0.02 -0.03  0.09

8.4 ACF, PACF: Kontrolná otázka

ACF aj PACF sa v R-ku počítajú numericky. Ktoré z nasledujúcich postupností sa po určitom počte členov vynulujú (a teda tie malé čísla sú v skutočnosti nuly) a ktoré majú hodnoty, ktoré konvergujú k nule, ale neexistuje index, od ktorého by boli nulové?

Príklad 1

ARMAacf(ar=c(0.05,-0.2), lag.max=10)

##             0             1             2             3             4 
##  1.0000000000  0.0416666667 -0.1979166667 -0.0182291667  0.0386718750 
##             5             6             7             8             9 
##  0.0055794271 -0.0074554036 -0.0014886556  0.0014166479  0.0003685635 
##            10 
## -0.0002649014

barplot(ARMAacf(ar=c(0.05,-0.2), lag.max=10))

Príklad 2

ARMAacf(ar=c(0.3,0.05), lag.max=10, pacf=TRUE)

##  [1]  3.157895e-01  5.000000e-02 -3.863413e-18  3.927096e-19 -4.622862e-19
##  [6] -1.450362e-19 -4.829267e-20  6.640242e-20  5.131096e-20 -2.263719e-20

barplot(ARMAacf(ar=c(0.3,0.05), lag.max=10,pacf=TRUE))

9 Stacionarita AR(2) procesu - simulácie

Vytvoríme funkciu, ktorá pre zadaný vektor parametrov \((\alpha_1, \alpha_2)\) vráti TRUE, resp. FALSE podľa toho, či je AR(2) proces s týmito koeficientami stacionárny alebo nie.

stacAR2 <- function(alpha12) {
# vstup: dvojzložkový vektor koeficientov
# výstup: TRUE/FALSE podľa stacionarity
....
}

Otestujte:

stacAR2(c(1.2, -0.5))

## [1] TRUE

stacAR2(c(1.2, 0.5))

## [1] FALSE

stacAR2(c(0.5, 0)) # AR(1), ale stacionárny

## [1] TRUE

stacAR2(c(0, 0))   # konštanta + biely šum -> stacionárny

## [1] TRUE

Spravíme sieť bodov \((\alpha_1, \alpha_2)\) a zakreslíme, ktoré sú stacionárne a ktoré nie sú. Na vytvorenie dvojíc bodov použijeme náhodne generované čísla z rovnomerného rozdelenia na štvorci \((-2,2)\times(-2,2)\).

set.seed(123)
alfa1 <- runif(1000, min = -2, max = 2)
alfa2 <- runif(1000, min = -2, max = 2)
df <- data.frame(alfa1, alfa2)
head(df)

##        alfa1       alfa2
## 1 -0.8496899 -0.90550906
## 2  1.1532205  0.37546773
## 3 -0.3640923 -1.35926075
## 4  1.5320696  1.41372096
## 5  1.7618691  1.39095663
## 6 -1.8177740 -0.08845275

Teraz pridáme do data framu tretí stĺpec, v ktorom bude výsledok overovania stacionarity. Namiesto for-cyklu použijeme funkciu apply. Tá umožňuje aplikovať zadanú funkciu na riadky alebo stĺpce objektu (matica, data frame).

Vsuvka: Ukážka použitia funkcie apply

set.seed(123)
dfUkazka <- data.frame(a = sample.int(5),
                       b = sample.int(5),
                       c = sample.int(5))
dfUkazka # nejaké čísla

##   a b c
## 1 3 3 2
## 2 2 1 3
## 3 5 2 1
## 4 4 5 4
## 5 1 4 5

f <- function(x) (x[1] + x[2]) * x[3] # nejaká funkcia

# POUŽITIE APPLY
apply(dfUkazka,   # vstupné dáta pre funkciu 
      MARGIN = 1, # po riadkoch
      FUN = f)    # použitá funkcia

## [1] 12  9  7 36 25

# priradenie ako nový stĺpec
dfUkazka$vysledok <- apply(dfUkazka, 1, f)
dfUkazka

##   a b c vysledok
## 1 3 3 2       12
## 2 2 1 3        9
## 3 5 2 1        7
## 4 4 5 4       36
## 5 1 4 5       25

# ĎALŠIE POUŽITIE APPLY: teraz nie je vstupom celý riadok
g <- function(x) x[2]/x[1] + x[1] # zase nejaká funkcia
dfUkazka$vysledok2 <- apply(dfUkazka[ , 3:4], 1, g)
dfUkazka

##   a b c vysledok vysledok2
## 1 3 3 2       12         8
## 2 2 1 3        9         6
## 3 5 2 1        7         8
## 4 4 5 4       36        13
## 5 1 4 5       25        10

Vrátime sa k overovaniu stacionarity.

Pomocou apply a vytvorenej funkcie stacAR2 pridajte do data framu df stĺpec stacionarita s informáciou o stacionarite procesu s danymi parametrami.

Mali by ste dostať

head(df)

##        alfa1       alfa2 stacionarita
## 1 -0.8496899 -0.90550906         TRUE
## 2  1.1532205  0.37546773        FALSE
## 3 -0.3640923 -1.35926075        FALSE
## 4  1.5320696  1.41372096        FALSE
## 5  1.7618691  1.39095663        FALSE
## 6 -1.8177740 -0.08845275        FALSE

Vykreslite body farebne odlíšené podľa toho, či je proces stacionárny (modrá) alebo nie (červená). Do obrázka pridajte čiernou nerovnosti odvodené teoreticky.

10 Ďalšie príklady

10.1 Nájdenie procesu s danou vlastnosťou

Nájdite príklad procesov, ktoré majú nasledovné vlastnosti. Pre každý proces dokážte, že má naozaj požadovanú vlastnosť.

Autoregrený proces, ktorého PACF pre lag 3 je nulová.
Autoregrený proces prvého rádu, ktorého PACF pre lag 1 je rovná 0.5.
Autoregresný proces, ktorého ACF je vždy kladná.
Autoregresný proces, ktorého ACF nie je monotónna.

10.2 Overovanie stacionarity AR procesov II.

10.2.1 Analyticky

Odvoďte podmienku stacionarity AR(2) procesu v tvare lineárnych nerovníc pre autoregresné koeficienty a zakreslite do roviny (na jednotlivých osiach sú hodnoty autoregresných koeficientov) tie body, ktoré zodpovedajú stacionárnym procesom. Porovnajte so simuláciami, ktoré sa robili v cvičení 9.

10.2.2 Simulačne

Ak riešime úlohy typu “pre aké hodnoty parametra je daný proces stacionárny” analyticky, je dobré mať numerickú kontrolu.

Príklad: Pre aké hodnoty parametra \(k\) je proces \(x_t = x_{t−1} + kx_{t−2} + u_t\) stacionárny?

Vytvorte funkciu armaRoots, ktorá vypočíta korene charakteristickej rovnice \(1-L-kL^2=0\) a vráti ich absolútne hodnoty. (Hint: funkcia Mod)

armaRoots <- function(k) {
  ...
  return(...)
}

Následne vytvorte vektor hodnôt kVektor (od -2 po 2) a pre každú jeho zložku nájdite absolútne hodnoty koreňov. Výsledok zapíšte do objektu root_abs_vals. Dostanete

## [[1]]
## [1] 0.7071068 0.7071068
## 
## [[2]]
## [1] 0.7254763 0.7254763
## 
## [[3]]
## [1] 0.745356 0.745356
## 
## [[4]]
## [1] 0.766965 0.766965
## 
## [[5]]
## [1] 0.7905694 0.7905694
## 
## [[6]]
## [1] 0.8164966 0.8164966

Keďže objekt root_abs_vals je zoznam, prvý a druhý koreň z neho extrahujeme nasledovným príkazom, čím dostávame dva vektory rovnakej dĺžky.

root1 <- sapply(root_abs_vals, function(x) x[1])  # Prvý koreň
root2 <- sapply(root_abs_vals, function(x) ifelse(length(x) == 2, x[2], x[1]))  # Druhý koreň (osetreny pripad dvojnasobneho korena)

Vektory root a root2 teraz môžeme vykresliť vzhľadom k parametru \(k\). Doplníme tiež horizontálnu čiaru znázorňujúcu hranicu stability a vertikálne čiary pre tri konkrétne hodnoty parametra \(k\).

Skontrolujte, či sa vaše výpočty zhodujú s numerickými na nasledujúcom obrázku. Prečo pre niektoré 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á)?

plot(kVektor, root1, type="p", col='black', ylim=c(0, 12), 
     xlab="parameter k", ylab="absolútna hodnota koreňov")
lines(kVektor, root2, type="p", col='black')  # Pridaj druhý koreň

abline(h=1, col="red")              # horizontálna čiara znázorňujúca hranicu stability (absolútna hodnota koreňov > 1)
abline(v=-1.5, col="blue", lty=2)   # Vertikálna čiara pre k = -1.5
abline(v=-0.2, col="green", lty=2)  # Vertikálna čiara pre k = -0.2
abline(v=0.5, col="brown", lty=2)   # Vertikálna čiara pre k = 0.5

10.3 Hľadanie AR modelu pre zadané dáta

Všetky nasledujúce dáta majú spoločné zadanie:

nájsť vhodný autoregresný model pre zadané dáta
vysvetliť, prečo sú rezíduá vyhovujúce - okomentovať 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 (napísať polynóm, ktorého korene overujete, aké absolútne hodnoty koreňov vyšli a prečo sme s tým spokojní)
spraviť predikcie

V prípade modelovania diferencií sa treba na základe prítomnosti, resp. neprítomnosti trendu treba rozhodnúť, či do modelu zahrnúť konštantu alebo nie.

10.3.1 Hladina Hurónskeho jazera

Dáta z knižnice datasets. Z popisu v helpe: Annual measurements of the level, in feet, of Lake Huron 1875–1972. Použijeme dáta od roku 1900.

x <- window(LakeHuron, start = 1900)
plot(x)

Dáta nediferencujeme, AR model hľadáme priamo pre tieto dáta.

10.3.2 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).

x <- lead
plot(x)

Nájdite model pre premennú x tak, že budete jej diferencie modelovať AR procesom.

10.3.3 Ceny ropy

Dáta z knižnice astsa, z popisu v helpe: Crude oil, WTI spot price FOB (in dollars per barrel), weekly data from 2000 to mid-2010.

Zoberieme dáta od roku 2006.

x <- window(oil, start = c(2006, 1))
plot(x)

Nájdite model pre premennú x tak, že budete jej diferencie modelovať AR procesom.

11 Dodatky

11.1 Pristupovanie ku koeficientom modelu

Pri AR(1) modeli sme hneď videli, či je stacionárny alebo nie. Pri modeloch vyššieho rádu použijeme funkciu polyroot. Hodnoty koeficientov však nemusíme odpisovať, tak by sme navyše dostali iba ich približné hodnoty.

model.ar3 <- sarima( x, 3,1,0, details=FALSE)

Pozrite si štruktúru str(model.ar3) objektu model.ar3. Následne vyberte z výstupu model.ar3 vektor odhadnutých parametrov modelu. Z tohto chceme tie zložky vektora, ktorých názvy obsahujú ar

##        ar1        ar2        ar3   constant 
## 0.16557547 0.05630823 0.08710237 0.02120166

Pre AR(p) model sa dá napísať koef[1:p]. Budeme mať aj všeobecnejšie modely (koeficienty budú s názvami ma1, ma2, … neskôr aj napríklad sar1, sma1). Napríklad

model.zlozitejsi <- sarima( x, 3,1,2,1,0,1, 1,details=FALSE)
koef.zlozitejsie <- model.zlozitejsi$fit$coef
koef.zlozitejsie

##          ar1          ar2          ar3          ma1          ma2         sar1 
## -0.114855736  0.884714786  0.028837190 -0.063466713 -0.585713976 -0.624833707 
##         sma1     constant 
##  0.960964358 -0.002350575

Počet členov jednotlivých typov je daný našou špecifikáciou modelu, takže sa znovu dá spočítať, ktoré zložky zodpovedajú členom ar, ktor ma atď. Alternatívou je použiť regulárne výrazy a automaticky vybrať napríklad zložky, ktorých názov začína ar, zložky, ktorých názov začína ma atď. Ak by sme chceli indexy tých, ktoré obsahujú ma, spravili by sme

grep("ma",                     # co sa ma hladat
    names(koef.zlozitejsie))   # kde sa to ma hladat

## [1] 4 5 7

My ale chceme iba tie, ktoré začínajú ma (lebo je tam aj sma1 koeficient, ktorý nechceme). Tu prichádzajú na rad regulárne výrazy

grep("^ma",                     # co sa ma hladat
    names(koef.zlozitejsie))   # kde sa to ma hladat

## [1] 4 5

Príslušné koeficienty potom sú

koef.zlozitejsie[grep("^ma",names(koef.zlozitejsie))]

##         ma1         ma2 
## -0.06346671 -0.58571398

Viac o regulárnych výrazov v R-ku: https://r4ds.had.co.nz/strings.html#matching-patterns-with-regular-expressions

11.2 Dáta Svetovej banky a grafy pomocou balíka `ggplot2`

Knižnica WDI v R-ku slúži na priamy prístup k dátam zo Svetovej banky (World Bank). Knižnica obsahuje informácie o viac ako 200 krajinách a viacerých regiónoch, pokrýva obdobie od roku 1960 po súčasnosť a obsahuje rôzne ukazovatele, ako sú HDP, inflácia, nezamestnanosť, úroveň vzdelania, prístup k vode, emisie CO2, atď.

11.2.1 Vyhľadávanie dát

Na vyhľadávanie názvov a popisov dostupných databáz v rámci knižnice WDI slúži funkcia WDI.search.

Ak napríklad hľadáme indikátory, ktoré majú v názve “gdp” použijeme

WDIsearch('gdp')

Dajú sa použiť aj regulárne výrazy (viac o regulárnych výrazov v R-ku: https://r4ds.had.co.nz/strings.html#matching-patterns-with-regular-expressions). Napríklad, ak je úlohou nájsť indikátory, ktoré hovoria o HDP na obyvateľa, vypíšme si tie, ktoré obsahujú reťazce “gdp” a “capita”.

WDIsearch('gdp.*capita')

Ak je veľa výsledkov, môžeme vypísať niekoľko prvých

WDIsearch('gdp.*capita')[1:6,]

11.2.2 Načítanie dát

Dáta z knižnice WDI načítame príkazom WDI. Ukážka:

data.hdp <- WDI( country = c("FI","US","FR", "DE"), 
                 indicator = "NY.GDP.PCAP.KD", 
                 start = 1975)

Argumenty funkcie:

country - vektor kódov štátov vo formáte iso2c. Kódy sa dajú nájsť napríklad tu: https://datahub.io/core/country-list alebo pomocou knižnice countrycode v R-ku (https://vincentarelbundock.github.io/countrycode/#/)
indicator - vektor indikátorov, ktoré sme našli pomocou funkcie WDIsearch
start - počiatočný rok (default je 1960)
ďalšie argumenty (pozri ?WDI)

head(data.hdp)

##   country iso2c iso3c year NY.GDP.PCAP.KD
## 1 Germany    DE   DEU 2023       42878.81
## 2 Germany    DE   DEU 2022       43361.18
## 3 Germany    DE   DEU 2021       42900.02
## 4 Germany    DE   DEU 2020       41601.97
## 5 Germany    DE   DEU 2019       43292.68
## 6 Germany    DE   DEU 2018       42928.74

11.2.3 Grafy pomocou balíka `ggplot`

ggplot(data.hdp, aes(year, data.hdp$NY.GDP.PCAP.KD, color = country)) +
  geom_line() +
  xlab('Year') + ylab('GDP per capita') +
  labs(title = 'GDP per capita (current USD)')

Viac o balíku ggplot2 napríklad tu: https://www.rdocumentation.org/packages/ggplot2/versions/3.5.0

11.3 Porovnanie rýchlosti for-cyklu a funkcie apply

Nasledovný kód porovnáva dva spôsoby určovania stability AR(2) modelu – pomocou vektorovej operácie apply a pomocou explicitného for cyklu.

Funkcia stacAR2 testuje, či je AR(2) proces stabilný na základe jeho parametrov. Vstupom tejto funkcie je vektor parametrov, výstupom je logický operátor TRUE/FALSE.

stacAR2 <- function(par) all(abs(polyroot(c(1, -par))) > 1)

Vo funkciách nazvaných verzia1 a verzia2 sa následne generujú dva náhodné vektory \(x\) a \(y\) s normálnym rozdelením, z ktorých sa vytvorí dataframe df. Rozdiel nastáva v tom, že zatiaľ čo verzia1 aplikuje funkciu stacAR2 po riadkoch dataframu df pomocou príkazu apply, verzia2 aplikuje na každý riadok explicitný for cyklus. Výsledkom je v oboch prípadoch logický vektor, ktorý pre každý riadok udáva, či sú príslušné parametre stabilné.

verzia1 <- function(N){
  set.seed(12345)
  x <- rnorm(N)
  y <- rnorm(N)
  df <- data.frame(x,y)
  stac <- apply(df, 1, stacAR2)
  return(stac)
}

verzia2 <- function(N){
  set.seed(12345)
  x <- rnorm(N)
  y <- rnorm(N)
  df <- data.frame(x,y)
  stac <- rep(0, N)
  for(i in 1:N) stac[i] <- stacAR2(as.numeric(df[i,]))
  return(stac)
}

Použijeme balík microbenchmark na meranie času, ktorý zaberie vykonanie funkcií verzia1 a verzia2 s veľkosťou vstupu N = 100 (10^2). Výsledky vizualizujeme pomocou boxplotu.

library(microbenchmark)
cas2 <- microbenchmark(verzia1(10^2), verzia2(10^2))
cas2

## Unit: microseconds
##           expr    min      lq     mean  median      uq    max neval
##  verzia1(10^2)  347.8  374.65  461.944  386.25  404.25 3858.8   100
##  verzia2(10^2) 1550.3 1593.75 1785.479 1633.25 1686.30 5537.2   100

boxplot(cas2)

Zo štatistík vidíme, že verzia1 je podstatne efektívnejšia ako verzia2 – napríklad priemerný čas vykonania je pri verzii1 asi štvornásobne kratší.

cas4 <- microbenchmark(verzia1(10^4), verzia2(10^4))
cas4

## Unit: milliseconds
##           expr      min        lq      mean   median        uq      max neval
##  verzia1(10^4)  22.0418  23.78085  25.86447  25.4669  27.71055  31.5694   100
##  verzia2(10^4) 157.9895 165.25480 170.35826 169.0584 171.08775 249.6299   100

boxplot(cas4)

Pre väčšie vstupy je tento rozdiel ešte výraznejší, čo ukazuje, že vektorové operácie sú v R efektívnejšie ako explicitné cykly.

Viac o funkcii apply a príbuzných funkciách napríklad tu: https://www.datacamp.com/community/tutorials/r-tutorial-apply-family

Cvičenie 2: Autoregresné procesy

1 Inštalácia knižníc v R-ku

2 Overovanie stacionarity AR procesu v R-ku

2.1 Príklad: Overenie stacionarity AR(2) procesu

2.2 Cvičenie: Overenie stacionarity AR procesu (súčasť kostry na skúške)

2.3 Cvičenie: Grafické znázornenie koreňov

2.4 Cvičenie: Nájdenie príkladu procesu s danou vlastnosťou

3 Hľadanie AR modelu

3.0.1 Úlohy

4 Konštrukcia predikcií

5 Zautomatizovanie testovania viacerých AR modelov

6 Woldova reprezentácia AR procesu

7 Autokorelačná funkcia AR procesu

7.1 Z Yule Wolkerových rovníc

7.2 Funkcia v R-ku

8 Parciálna autokorelačná funkcia AR procesu

8.1 Opakovanie z prednášky

8.2 Výpočet v R-ku

8.3 Výberová PACF z dát

8.4 ACF, PACF: Kontrolná otázka

9 Stacionarita AR(2) procesu - simulácie

10 Ďalšie príklady

10.1 Nájdenie procesu s danou vlastnosťou

10.2 Overovanie stacionarity AR procesov II.

10.2.1 Analyticky

10.2.2 Simulačne

10.3 Hľadanie AR modelu pre zadané dáta

10.3.1 Hladina Hurónskeho jazera

10.3.2 Indikátor od Boxa a Jenkinsa

10.3.3 Ceny ropy

11 Dodatky

11.1 Pristupovanie ku koeficientom modelu

11.2 Dáta Svetovej banky a grafy pomocou balíka ggplot2

11.2.1 Vyhľadávanie dát

11.2.2 Načítanie dát

11.2.3 Grafy pomocou balíka ggplot

11.3 Porovnanie rýchlosti for-cyklu a funkcie apply

11.2 Dáta Svetovej banky a grafy pomocou balíka `ggplot2`

11.2.3 Grafy pomocou balíka `ggplot`