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 65 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':    65 obs. of  8 variables:
##  $ indicator      :'data.frame': 65 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': 65 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  "2024" "2023" "2022" "2021" ...
##  $ value          : num  2.95 4.12 8 4.7 1.23 ...
##  $ 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 diferencie bez konst.
  # vystup: vektor ar koeficientov

Mali by ste dostať nasledovné hodnoty koeficientov:

## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##     Estimate     SE t.value p.value
## ar1   0.1325 0.1233  1.0746  0.2867
## 
## sigma^2 estimated as 3.104685 on 63 degrees of freedom 
##  
## AIC = 4.033566  AICc = 4.034574  BIC = 4.101031 
##  
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##     Estimate     SE t.value p.value
## ar1   0.1835 0.1119  1.6397  0.1061
## ar2  -0.4455 0.1140 -3.9081  0.0002
## 
## sigma^2 estimated as 2.498012 on 62 degrees of freedom 
##  
## AIC = 3.85429  AICc = 3.857364  BIC = 3.955488 
##  
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##     Estimate     SE t.value p.value
## ar1   0.1404 0.1248  1.1256  0.2648
## ar2  -0.4212 0.1179 -3.5725  0.0007
## ar3  -0.1010 0.1320 -0.7652  0.4471
## 
## sigma^2 estimated as 2.47466 on 61 degrees of freedom 
##  
## AIC = 3.87644  AICc = 3.88269  BIC = 4.01137 
##  
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##     Estimate     SE t.value p.value
## ar1   0.1215 0.1235  0.9837  0.3292
## ar2  -0.4921 0.1267 -3.8850  0.0003
## ar3  -0.0695 0.1318 -0.5276  0.5997
## ar4  -0.1861 0.1298 -1.4339  0.1568
## 
## sigma^2 estimated as 2.393485 on 60 degrees of freedom 
##  
## AIC = 3.876218  AICc = 3.886811  BIC = 4.04488 
##  
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##     Estimate     SE t.value p.value
## ar1   0.1436 0.1246  1.1527  0.2537
## ar2  -0.4820 0.1262 -3.8181  0.0003
## ar3  -0.0127 0.1435 -0.0883  0.9300
## ar4  -0.2051 0.1301 -1.5770  0.1201
## ar5   0.1256 0.1299  0.9670  0.3375
## 
## sigma^2 estimated as 2.355923 on 59 degrees of freedom 
##  
## AIC = 3.893015  AICc = 3.909179  BIC = 4.09541 
## 
## [[1]]
##       ar1 
## 0.1325166 
## 
## [[2]]
##        ar1        ar2 
##  0.1835071 -0.4455040 
## 
## [[3]]
##        ar1        ar2        ar3 
##  0.1404300 -0.4211714 -0.1010244 
## 
## [[4]]
##         ar1         ar2         ar3         ar4 
##  0.12148906 -0.49205391 -0.06953632 -0.18610615 
## 
## [[5]]
##         ar1         ar2         ar3         ar4         ar5 
##  0.14357752 -0.48196095 -0.01266644 -0.20514028  0.12560244

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.325166e-01  0.1835070966  0.1404299897  0.121489056  0.143577518
##  [2,] 1.756066e-02 -0.4118291453 -0.4014508074 -0.477294324 -0.461346446
##  [3,] 2.327080e-03 -0.1573267163 -0.2165452641 -0.187301522 -0.148104176
##  [4,] 3.083768e-04  0.1546009625  0.1244832844  0.017545405 -0.005872356
##  [5,] 4.086506e-05  0.0984600551  0.1492401981  0.104873452  0.172529815
##  [6,] 5.415301e-06 -0.0508072283 -0.0095946349  0.105959358  0.142152027
##  [7,] 7.176175e-07 -0.0531878353 -0.0767789298 -0.005092568 -0.090232523
##  [8,] 9.509627e-08  0.0128744782 -0.0218179857 -0.063314230 -0.101050015
##  [9,] 1.260184e-08  0.0260579515  0.0302424816 -0.032071786 -0.008950915
## [10,] 1.669953e-09 -0.0009538125  0.0211926109  0.007892074  0.041068995
## [11,] 2.212966e-10 -0.0117839530 -0.0075570401  0.022090245  0.047855486
## [12,] 2.932548e-11 -0.0017375117 -0.0130421862  0.012813719 -0.003413300
## [13,] 3.886115e-12  0.0049309525 -0.0007896766 -0.003892894 -0.034930682
## [14,] 5.149749e-13  0.0016789332  0.0061455471 -0.009782822 -0.013525504
## [15,] 6.824275e-14 -0.0018886629  0.0025131878 -0.004275141  0.010277779

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=2014)

Zo skrátených 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 +/- dvakrát š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 = 2014)) +
  autolayer(window(x, start = 2014, end = 2024), 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 Dodatky

5.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)  327.2  363.15  422.255  404.20  428.55 2462.0   100
##  verzia2(10^2) 1531.2 1641.05 1904.292 1733.35 1831.30 7266.1   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)  21.5066  23.70885  25.09232  24.9673  26.7382  31.0364   100
##  verzia2(10^4) 154.9174 160.85120 171.04452 171.4925 174.9166 251.9981   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

5.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)
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ar1        0.1656 0.0655  2.5283  0.0121
## ar2        0.0563 0.0665  0.8473  0.3977
## ar3        0.0871 0.0658  1.3241  0.1868
## constant   0.0212 0.3308  0.0641  0.9490
## 
## sigma^2 estimated as 12.19486 on 228 degrees of freedom 
##  
## AIC = 5.382287  AICc = 5.383047  BIC = 5.456571 
##

Pozrite si štruktúru str(model.ar3) objektu model.ar3. Následne vyberte z výstupu model.ar3 vektor odhadnutých parametrov modelu.

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

Z tohto chceme tie zložky vektora, ktorých názvy obsahujú ar. 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)
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ar1       -0.1149 0.2058 -0.5580  0.5774
## ar2        0.8847 0.0668 13.2472  0.0000
## ar3        0.0288 0.1648  0.1750  0.8612
## ma1       -0.0635 0.3945 -0.1609  0.8723
## ma2       -0.5857 0.3336 -1.7556  0.0805
## sar1      -0.6248 0.5016 -1.2457  0.2142
## sma1       0.9610 0.0542 17.7252  0.0000
## constant  -0.0024 0.4656 -0.0050  0.9960
## 
## sigma^2 estimated as 11.85058 on 224 degrees of freedom 
##  
## AIC = 5.38931  AICc = 5.392093  BIC = 5.523019 
## 
koef.zlozitejsie <- model.zlozitejsi$fit$coef
koef.zlozitejsie
##          ar1          ar2          ar3          ma1          ma2         sar1 
## -0.114855734  0.884714782  0.028837186 -0.063466694 -0.585713988 -0.624833726 
##         sma1     constant 
##  0.960964357 -0.002350561

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.06346669 -0.58571399

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