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 10
packages <- c("astsa", "urca", "forecast", "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 Zhlukovanie

Zhlukovanie je metóda strojového učenia, ktorá sa používa na zoskupenie podobných objektov na základe ich vlastností. Cieľom je, aby objekty v rovnakom zhluku boli čo najviac podobné, zatiaľ čo objekty z rôznych zhlukov sa od seba odlišovali. Táto technika je často využívaná pri analýze dát, keď chceme identifikovať vzory alebo prirodzené skupiny bez predchádzajúceho zadania kategórií. Jednotlivým technikám zhlukovania sa budete podrobne venovať na predmete Analýza zhlukov a klasifikácia dát v letnom semestri.

V prípade časových radov zhlukovanie pomáha identifikovať rady s podobnou dynamikou, sezónnosťou, trendmi alebo inými charakteristikami. V tomto cvičení si pomocou funkcie `hclust˙ ukážeme tzv. hierarchické zhlukovanie. Cieľom hierarchického zhlukovania je vytvoriť tzv. “strom podobností”, nazývaný dendrogram. Tento výstup je často informatívnejšíako výstupy iných metód zhlukovania a keďže zhluky dokáže identifikovať na rôznej granularite, netreba algoritmus zbiehať opakovane.

2.1 Modelový príklad

Princíp zhlukovania si najskôr ukážeme na modelovom príklade. Uvažujme nasledovných 5 dátových bodov v rovine uložených v premennej data.

data <- cbind(x1 = c(1,1,2,7,8), x2 = c(1,2,5,5,3))
data
##      x1 x2
## [1,]  1  1
## [2,]  1  2
## [3,]  2  5
## [4,]  7  5
## [5,]  8  3

Dátové body vykreslíme.

plot(data, pch = 21, bg = "red")
text(data[, 1] + 0.2, data[, 2], 1:5)

Naším cieľom je aplikovať postup hierarchického zhlukovania. V R-ku na to existuje napríklad funkcia hclust. V popise v dokumentácii je uvedené, že povinným argumentom funkcie je argument d, ktorý je v helpe opísaný ako “a dissimilarity structure as produced by dist”.

Aplikujme teda funkciu dist na naše dáta, čím získame maticu euklidovských vzdialeností medzi bodmi v rovine.

dist(data)
##          1        2        3        4
## 2 1.000000                           
## 3 4.123106 3.162278                  
## 4 7.211103 6.708204 5.000000         
## 5 7.280110 7.071068 6.324555 2.236068

Ako bolo uvedené, táto matica bude vstupom do funkcie hclust. Pomocou defaultne zvolenej metódy hierarchického zhlukovania (method="complete") dostávame nasledový dendrogram.

zhlukovanie <- hclust(dist(data))
plot(zhlukovanie, hang = -1)

Teda ak chceme dátové body rozdeliť do dvoch zhlukov, prvý zhluk budú tvoriť body 1, 2, 3 a druhý zhluk body 4 a 5. Farebne rozlíšené zhluky môžeme vykresliť pričom dostávame nasledovný výsledok.

plot(data, pch = 21, bg = c("red", "red", "red", "green", "green"))
text(data[, 1]+0.2, data[, 2], 1:5)

2.2 Reálne dáta

Postup uvedený vyššie pre modelový príklad zopakujeme na reálnych dátach. Budeme pracovať s dátami o výmenných kurzoch. Načítajte dáta kurzy.Rdata do R-ka pomocoou príkazu:

load(url("https://anna-hlubinova.github.io/kurzy.Rdata"))

Dáta obsahujú 8 časových radov - CNY, CZK, GBP, HUF, CHF, JPY, KRW, PLN. Každý časový rad zodpovedá jednému výmennému kurzu s eurom (napr. USD je výmenný kurz EUR/USD).

Naším cieľom je pre jednotlivé časové rady nájsť vhodné ARMA modely a následne pre ne nájsť Woldovu reprezentáciu. Potom budeme môcť jednotlivé časové rady zhlukovať.

Inšpiráciou je článok: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1467-9892.1990.tb00048.x , kde síce autor nepoužíva Woldovu reprezentáciu MA(\(\infty\)), ale používa AR(\(\infty\)) repzerentáciu.

Zobrazte si priebeh jednotlivých výmenných kurzov:

# Nastavenie mriežky na 2 riadky a 4 stĺpce pre zobrazenie 8 grafov
par(mfrow = c(2, 4))

# Vykreslenie grafov jednotlivých kurzov
plot(CHF, type = "l", main = "CHF", xlab = "Čas", ylab = "Kurz")
plot(CZK, type = "l", main = "CZK", xlab = "Čas", ylab = "Kurz")
plot(CNY, type = "l", main = "CNY", xlab = "Čas", ylab = "Kurz")
plot(GBP, type = "l", main = "GBP", xlab = "Čas", ylab = "Kurz")
plot(HUF, type = "l", main = "HUF", xlab = "Čas", ylab = "Kurz")
plot(PLN, type = "l", main = "PLN", xlab = "Čas", ylab = "Kurz")
plot(JPY, type = "l", main = "JPY", xlab = "Čas", ylab = "Kurz")
plot(KRW, type = "l", main = "KRW", xlab = "Čas", ylab = "Kurz")

Vašou prvou úlohou bude pre každý výmenný kurz nájsť ARIMA model vhodný na modelovanie jednotlivých výmenných kurzov.

Urobiť to môžete napríklad pomocou funkcie auto.arima z balíka forecast.

Vsuvka: Funkcia auto.arima

Popis funkcie z helpu: Returns best ARIMA model according to either AIC, AICc or BIC value. The function conducts a search over possible model within the order constraints provided.

Táto funkcia teda pri výbere nekontroluje rezíduá, ale model vyberá na základe informačného kritéria.

Môže sa hodiť ako inšpirácia pre vyskúšanie modelu pre dáta (napr. ak na základe ACF a PACF nemáme konkrétny tip). Výstup však je potrebné skontrolovať - najmä rezíduá.

Funkcia auto.arima umožňuje zadať viaceré vstupné parametre na doladenie vyhľadávania. Detailný popis týchto parametrov nájdete v dokumentácii. Príklad použitia funkcie:

auto.arima(
  y,
  d = NA,
  D = NA,
  max.p = 5,
  max.q = 5,
  max.P = 2,
  max.Q = 2,
  max.order = 5,
  max.d = 2,
  max.D = 1,
  start.p = 2,
  start.q = 2,
  start.P = 1,
  start.Q = 1,
  stationary = FALSE,
  seasonal = TRUE,
  ic = c("aicc", "aic", "bic"),
  stepwise = TRUE,
  nmodels = 94,
  trace = FALSE,
  approximation = (length(x) > 150 | frequency(x) > 12),
  method = NULL,
  truncate = NULL,
  xreg = NULL,
  test = c("kpss", "adf", "pp"),
  test.args = list(),
  seasonal.test = c("seas", "ocsb", "hegy", "ch"),
  seasonal.test.args = list(),
  allowdrift = TRUE,
  allowmean = TRUE,
  lambda = NULL,
  biasadj = FALSE,
  parallel = FALSE,
  num.cores = 2,
  x = y,
  ...
)

Vytvorte dataframe kurzy, stĺpce ktorého bude tvoriť 8 výmenných kurzov.

kurzy <- data.frame(CHF, CZK, CNY, GBP, HUF, PLN, JPY, KRW)

Pomocou funkcií apply a auto.arima nájdite SARIMA modely pre kurzy v stĺpcoch dataframu kurzy. Pritom

  • nemá sa použiť aproximácia

  • majú sa prejsť všetky možné modely s prípustným počtom parametrov

  • uvažované kombinácie počtu AR a MA členov nechajte defaultné

  • ako kritérium výberu modelu použite Bayesovo informačné kritérium

Výstup pre kontrolu:

## $EUR.CHF
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.4518
## s.e.  0.0654
## 
## sigma^2 = 8.014e-06:  log likelihood = 792.24
## AIC=-1580.47   AICc=-1580.41   BIC=-1574.11
## 
## $EUR.CZK
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.6042
## s.e.  0.0605
## 
## sigma^2 = 0.0006962:  log likelihood = 395.17
## AIC=-786.34   AICc=-786.27   BIC=-779.97
## 
## $EUR.CNY
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.4259
## s.e.  0.0688
## 
## sigma^2 = 0.0005215:  log likelihood = 420.62
## AIC=-837.23   AICc=-837.16   BIC=-830.87
## 
## $EUR.GBP
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.2679
## s.e.  0.0713
## 
## sigma^2 = 8.142e-06:  log likelihood = 790.89
## AIC=-1577.77   AICc=-1577.7   BIC=-1571.41
## 
## $EUR.HUF
## Series: kurz 
## ARIMA(1,1,3) 
## 
## Coefficients:
##           ar1     ma1     ma2      ma3
##       -0.9135  1.6436  0.5034  -0.2140
## s.e.   0.0404  0.0779  0.1407   0.0764
## 
## sigma^2 = 3.341:  log likelihood = -359.36
## AIC=728.72   AICc=729.07   BIC=744.63
## 
## $EUR.PLN
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.6797
## s.e.  0.0510
## 
## sigma^2 = 0.0001172:  log likelihood = 553.34
## AIC=-1102.69   AICc=-1102.62   BIC=-1096.32
## 
## $EUR.JPY
## Series: kurz 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       0.4525
## s.e.  0.0651
## 
## sigma^2 = 0.3939:  log likelihood = -169.24
## AIC=342.48   AICc=342.55   BIC=348.85
## 
## $EUR.KRW
## Series: kurz 
## ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.2736
## s.e.  0.0719
## 
## sigma^2 = 23.2:  log likelihood = -531.91
## AIC=1067.83   AICc=1067.89   BIC=1074.19

Aké modely vybrala funkcia auto.arima za najlepšie? Majú tieto modely vyhovujúce rezíduá? Sú stacionárne a invertovateľné?

Nájdené modely uložte do zoznamu, napr. takto:

chf <- sarima(CHF, p, d, q)
czk <- sarima(CZK, p, d, q)
cny <- sarima(CNY, p, d, q)
gbp <- sarima(GBP, p, d, q)
huf <- sarima(HUF, p, d, q)
pln <- sarima(PLN, p, d, q)
jpy <- sarima(JPY, p, d, q)
krw <- sarima(KRW, p, d, q)

models <- list(chf = chf, czk = czk, cny = cny, gbp = gbp, huf = huf, pln = pln, jpy = jpy, krw = krw)

Následne vytvorte funkciu calculate_wold na výpočet Woldovej reprezentácie. Ako by ste zvolili parameter lag.max?

calculate_wold <- function(...) {
  ....
  return(wold)
}

Pomocou funkcie laplly vypočítajte koeficienty Woldovej reprezentácie pre každý model. Výsledok uložte do objektu wold_vectors.

wold_vectors <- lapply(...)

Nasledovným príkazom dostaneme maticu Woldových koefientov.

wold_matrix <- do.call(rbind, wold_vectors)

Na vytvorenie dendrogramu chceme použiť funkciu hclust.

clusters <- hclust(...)

Dendrogram vykreslite. Mali by ste dostať:

Na určenie počtu zhlukov možno použiť napríklad tzv. ľakťový diagram (viac o ňom sa dozviete napr. na predmete Analýza zhlukov a klasifikácia dát). Na osi x sú v diagrame zobrazené počty zhlukov. Na osi y tohto diagramu sa zobrazuje výška spojení, teda vzdialenosť medzi zhlukmi, ktoré sa spájajú na každej úrovni hierarchie.

Pre naše dáta vykreslíme diagram pomocou nasledujúcich riadkov kódu.

# Extrahujte výšky spojení (vzdialenosti medzi zhlukmi)
heights <- rev(clusters$height)  # Vzdialenosti medzi spojenými zhlukmi, v obrátenom poradí
# Vykreslenie výšok spojení, aby sme identifikovali zlom
plot(1:length(heights), heights, type = "b", pch = 19,
     xlab = "Počet zhlukov",
     ylab = "Výška spojenia (Vzdialenosť medzi zhlukmi)",
     main = "Elbow method pre stanovenie počtu zhlukov")

Princíp tejto metódy na určovanie počtu zhlukov spočíva v tom, že hľadáme tzv. ľakeť. Teda bod, kde dochádza k výraznej zmene sklonu alebo k veľkému skoku vo výške spojenia. V našom prípade, by sme za ľakeť pravdepodobne zvolili bod 2 alebo 4.

Výsledné zaradenie do zhlukov si potom môžeme vypísať napríklad nasledovne.

optimal_clusters <- cutree(clusters, k = 4)
print(optimal_clusters)
## chf czk cny gbp huf pln jpy krw 
##   1   2   1   3   4   2   1   3

Ako by vyzerali zhluky, ak by sme namiesto koeficientov Woldovej reprezentácie počítali s koeficientami AR(\(\infty\)) reprezentáciou procesu ako je uvedené v článku?

3 Backcasting

Rovnako ako môžeme pri analýze náhodných procesov využívať forecasting na predpovede budúceho vývoja, môžeme aplikovať aj backcasting, teda odhady vývoja smerom do minulosti. Táto technika je užitočná, ak chceme zrekonštruovať priebeh udalostí alebo hodnôt v období, pre ktoré nemáme priamo dostupné údaje. Ako sa dajú takéto predikcie do minulosti zostrojiť?

Ilustrujte na dátach BJsales z knižnice datasets, ktoré sme modelovali na Cvičení 8. Ich popis v helpe: The sales time series BJsales and leading indicator BJsales.lead each contain 150 observations.

x <- BJsales
plot(x)

Vynechajte z dát prvých 15 pozorovaní. Tieto sa pokúste zrekonšturovať. Zrekonštuované dáta porovnajte so skutočnosťou. Výstup pre kontrolu: