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","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 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

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.

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.

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.

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.

4 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-Walkerových rovníc a diferenčnej rovnice z prednášky,
pomocou funkcie v R-ku.

4.1 Z Yule-Walkerových rovníc

Riešenie: https://anna-hlubinova.github.io/ACF_Yule-Walker.html

Pri tomto spôsobe riešenia je prvým krokom zostavenie Yule-Walkerový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

4.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” :).

5 Parciálna autokorelačná funkcia AR procesu

5.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?

5.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

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

5.3.1 Príklad

Pozrieme sa na mesačné dáta o exporte malých ostrovných ekonomík. Link na dataset: https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/UA2LRO.

Dáta do R-ka načítajte príkazom

load(url("https://anna-hlubinova.github.io/NLucic_SIHTD_Continuous_Foreign_Trade_series_1900-2021-14032024.RData"))

Vytvorí sa objekt x (skontrolujte jeho štruktúru).

Na ukážku si zoberme ostrov Seychelles. Načítajte dáta pre Seychelles do premennej data. Skontrolujte, či sú v dátach nejaké NA hodnoty na konci - ak áno, treba ich vynechať. Z vektora hodnôt spravte časový rad pomocou funkcie ts. Vykreslite ich priebeh.

S rastúcimi hodnotami rastie disperzia dát, preto dáta zlogaritmujte

plot(log(data))

V týchto dátach je rastúci trend, preto nie sú stacionárne a kvôli tomu ich zdiferencujte

plot(diff(log(data)))

Zobrazte výberovú parciálnu autokorelačnú funkciu. Vyskúšajte oba spôsoby - funkcie pacf aj acf2. Na porovnanie

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

Úloha: Nájdite vhodný AR(p) model pre naše dáta a ukážte, že daný proces je stacionárny.

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

6 Dodatok: 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ď.

6.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,]

6.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

6.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

Cvičenie 5: AR(p) procesy