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("priceR","astsa","ggplot2","forecast","datasets")

# 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 Numerická kontrola - Cvičenie 3, Príklad 5

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.

3 Numerické simulácie - Cvičenie 3, Príkad 6

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

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

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

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

4 Balík priceR, modelovanie inflácie

4.1 Príprava dát a hľadanie AR modelu

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

Príklad 1. Diferencie sú “ar(0)” s danou hodnotou \(x_0 = u_0\):

\[x_t - x_{t-1} = c + u_t \Rightarrow x_t = ct + \sum_{i=0}^t u_i \] \[x_t - x_{t-1} = u_t \Rightarrow x_t = \sum_{i=0}^t u_i\]

Príklad 2. Simulácie priebehu, ak sú diferencie AR(1)

set.seed(123)
x1 <- arima.sim(model = list(ar = c(0.8)), n = 200)
y1 <- ts(cumsum(x1)) # diferencie su ar(1) bez konst.
plot(y1)

set.seed(123)
x2 <- 0.5 + arima.sim(model = list(ar = c(0.8)), n = 200)
y2 <- ts(cumsum(x2)) # diferencie su ar(1) s konst.
plot(y2)

  • 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)

4.1.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.2 Woldova reprezentácia pre autoregresné procesy rôznych rádov

Najskôr si zapíšeme vektory AR koeficientov pre procesy dané modelmi rádov 1 až 5 (budeme ich ešte potrebovať a nechceme ich odhadovať znovu).

Použijeme podobnú funkciu ako už známu sapply. Rozdiel bude v tom, že teraz sa nedajú zapísať výsledky do matice, pretože vektory koeficientov budú mať rôzne dĺžky. Riešením je použitie funkcie lapply, ktorej výstupom je list (zoznam).

Doplňte nasledovný úryvok kódu:

koefAR <- lapply(1:5, function(p) ... )
  # function(p) ... 
  # vstup: rad procesu p pre pre diferencie bez konst.
  # vystup: vektor ar koeficientov

Mali by ste dostať nasledovné hodnoty koeficientov:

## [[1]]
##      ar1 
## 0.119085 
## 
## [[2]]
##        ar1        ar2 
##  0.2090878 -0.4685810 
## 
## [[3]]
##        ar1        ar2        ar3 
##  0.1646598 -0.4443555 -0.1102319 
## 
## [[4]]
##         ar1         ar2         ar3         ar4 
##  0.14465097 -0.51337332 -0.07845785 -0.18312873 
## 
## [[5]]
##        ar1        ar2        ar3        ar4        ar5 
##  0.1640901 -0.5019967 -0.0230789 -0.2018162  0.1202135

Teraz použijeme jednotlivé zložky zoznamu ako vstup do funkcie, ktorá pre zadaný vektor AR koeficientov vypočíta koeficienty Woldovej reprezentácie (vypočítajme prvých 15). Tu už bude výstupom matica (rozmery sú dané počtom modelov a počtom Woldových koeficientov), preto môžeme použiť funkciu sapply.

koefWold <- sapply(koefAR, function(ARvektor) ... ) 
  # vstup funkcie: vektor AR koef. modelu
  # vystup: prvych 15 hodnot Woldovej reprezetacie

Dostávame maticu Woldových koeficientov:

##               [,1]         [,2]         [,3]         [,4]         [,5]
##  [1,] 1.190850e-01  0.209087823  0.164659832  0.144650972  0.164090076
##  [2,] 1.418125e-02 -0.424863291 -0.417242667 -0.492449415 -0.475071156
##  [3,] 1.688774e-03 -0.186808324 -0.252102560 -0.223951089 -0.183406041
##  [4,] 2.011078e-04  0.160023524  0.125742146  0.025937911  0.002785868
##  [5,] 2.394893e-05  0.120993803  0.178721318  0.130869233  0.190587979
##  [6,] 2.851959e-06 -0.049685654  0.001343761  0.113366888  0.149710792
##  [7,] 3.396257e-07 -0.067084063 -0.093055344 -0.011809296 -0.091268454
##  [8,] 4.044434e-08  0.009255293 -0.035620383 -0.074925458 -0.117139253
##  [9,] 4.816315e-09  0.033369487  0.035336285 -0.037635902 -0.014988921
## [10,] 5.735511e-10  0.002640299  0.031904253  0.013186458  0.051147558
## [11,] 6.830136e-11 -0.015084254 -0.006522020  0.029269814  0.055037352
## [12,] 8.133670e-12 -0.004391128 -0.019145933  0.014138167 -0.003630003
## [13,] 9.685984e-13  0.006150063 -0.003771338 -0.007123608 -0.040461362
## [14,] 1.153456e-13  0.003343502  0.008605548 -0.012999861 -0.018211535
## [15,] 1.373593e-14 -0.002182717  0.005203297 -0.004692766  0.012448117

Teraz môžeme koeficienty Woldovej reprezentácie pre všetých 5 modelov vykresliť na jednom obrázku. Pri kreslení sa môže zísť funkcia matplot (spoločná x-ová os 1:15, na y-ovej stĺpce zadanej matice).

Dostávame:

4.2.1 Úlohy

  1. Porovnajte Woldove reprezentácie.
  2. Ktoré z týchto modelov mali dobré rezíduá?

4.3 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 5 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 piatich 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:

Iný postup (len doplnenie dát do obrázku zo sarima.for): Po vykreslení grafu do neho štandardným spôsobom môžeme pridávať lines, points a pod.

Krajšie grafy sa tiež dajú spraviť pomocou funkcií autoplot a autolayer

pallete <-  c("#D55E00","#009E73", "gray", "gray")

UB <- predikcie$pred + 2 * predikcie$se
LB <- predikcie$pred - 2 * predikcie$se

autoplot(window(x, start = 1960, end = 2013)) +
  autolayer(window(x, start = 2013, end = 2023), series = "Data") +
  autolayer(predikcie$pred, series = "Forecasts") +
  autolayer(UB, series = "UB") +
  autolayer(LB, series = "LB") +
  scale_colour_manual(values = pallete) +
  theme(axis.title.y = element_blank())

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

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

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

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

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

library(astsa)
x <- lead
plot(x)

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

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

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

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

6 Dodatky

6.1 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)  355.4  380.6  440.302  391.3  408.6 2826.0   100
##  verzia2(10^2) 1597.2 1658.3 1818.211 1691.0 1742.3 6220.7   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.1782  24.31025  25.64945  25.1309  27.01865  30.5049   100
##  verzia2(10^4) 157.5383 162.72055 170.02484 169.8969 171.60480 238.5488   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

6.2 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