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("quantmod","astsa","plotly","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 Opakovanie - základné pojmy

  • Čo je biely šum? (stredná hodnota, disperzia, autokovariancia)
  • Kedy považujeme proces za stacionárny? (slabá a silná stacionarita)
  • Ako je definovaná autokorelačná funkcia?

3 Generovanie nezávislých dát

Vygenerujeme dáta - ide o biely šum posunutý o konštantu.

N <- 250        # pocet dat
set.seed(12345) # aby sme mali rovnake vysledky
x <- rnorm(N, mean = 10, sd = 1) # nezavisle N(10, 1)

Zobrazte priebeh vygenerovaných dát

4 Testovanie nulovosti autokorelácií - každá samostatne

Časový rad vieme veľmi dobre reprezentovať pomocou autokorelačnej funkcie (ACF).

V jazyku R máme niekoľko možností ako autokorelačnú funkciu vypočítať. Súčasťou základnej syntaxe jazyka R je funkcia acf. Zjednodušením tejto funkcie je potom napríklad funkcia acf1 z balíka astsa, ktorá vynecháva lag 1. Ďalšou možnosťou je použiť funkciu acf2. Tá bude užitočná neskôr, zatiaľ ju však používať nebudeme.

acf(x = časový rad,
    lag.max = maximálny počet lagov ktorý nás zaujíma,
    type = jedno z c("correlation", "covariance", "partial"),
    plot = TRUE/FALSE podľa toho či chceme vykresliť aj graf, 
    a ďaľšie argumenty)

acf1(series = časový rad, 
     max.lag = maximálny počet lagov ktorý nás zaujíma,
     plot = TRUE/FALSE podľa toho či chceme vykresliť aj graf,
     pacf = TRUE/FALSE podľa toho či chceme autokorelácie (FALSE) alebo parciálne autokorelácie (TRUE),
     a ďaľšie argumenty)

Zobrazte odhadnutú autokorelačnú funkciu.

acf(x)

Hodnoty autokorelácií vieme dostať pomocou príkazu

acf(x, plot = FALSE)$acf
## , , 1
## 
##                [,1]
##  [1,]  1.000000e+00
##  [2,] -4.802202e-03
##  [3,] -5.694779e-02
##  [4,] -2.009150e-02
##  [5,] -3.201563e-02
##  [6,]  6.516266e-02
##  [7,] -1.014845e-01
##  [8,]  5.320290e-02
##  [9,] -8.837230e-02
## [10,] -6.278643e-02
## [11,]  8.468051e-02
## [12,] -3.144206e-04
## [13,]  4.505565e-02
## [14,]  1.599475e-02
## [15,] -4.271585e-02
## [16,] -3.798188e-05
## [17,] -1.013824e-01
## [18,]  5.826572e-02
## [19,]  7.121903e-02
## [20,] -1.792323e-02
## [21,]  2.131789e-03
## [22,]  1.496130e-02
## [23,] -1.538814e-02
## [24,] -3.465120e-02

4.1 Otázky

  • Akú hypotézu vieme pomocou tohto výstupu testovať pre každú autokoreláciu?
  • Kedy sa táto nulová hypotéza zamieta?
  • Čo dostávame v našom prípade?

5 Testovanie nulovosti autokorelácií - Ljung-Boxov test

Ak nám vyjdú hodnoty autokorelácií na hranici kritickej oblasti pre zamietanie nulovej hypotézy, je dobré pozrieť sa na autokorelácie aj skupinovo, nielen jednotlivo. Vždy pracujeme len s odhadmi skutočných autokorelácií, preto je užitočné testovať nulovosť skupiny prvých autokorelácií. Napríklad, môžeme testovať nulovosť prvých štyroch autokorelácií, teda \(\rho_1=\rho_2=\rho_3=\rho_4=0\).

Na to nám poslúži Ljung-Boxov test, v Rku príkaz Box.test.

Box.test(x = testovaný časový rad,
         lag = počet prvých autokorelácií, ktoré chceme testovať, 
        type = jedno z c("Box-Pierce", "Ljung-Box") pričom budeme používať druhý menovaný,
         fitdf = počet parametrov modelu, ktoré sa majú odpočítať (ak x sú reziduá))

Argument fitdf bude pre nás mať zmysel neskôr, pri modelovaní ARIMA modelov. Slúži na to, že ak testujeme namiesto pôvodného časového radu reziduá z modelu, stupne voľnosti pre chi-kvadrát štatistiku sa zmenia. To zohľadňuje počet parametrov modelu, ktoré sú už v reziduách vysvetlené.

Otestujte nulovosť prvých troch autokorelácií časového radu x.

## 
##  Box-Ljung test
## 
## data:  x
## X-squared = 0.93263, df = 3, p-value = 0.8175

5.1 Otázky

  • Ako sa vypočíta testovacia štatistika?
  • Aké je pravdepodobnostné rozdelenie štatistiky za platnosti nulovej hypotézy \(H_0\)?
  • Pre aké hodnoty testovej štatistiky zamietame hypotézu \(H_0\)?
  • Ako sa určí p-hodnota?
  • Čo dostávame v našom prípade?

6 Ljung-Boxov test pre viac lagov a grafické znázornenie I. (simulované dáta)

Cieľom je zrekonštruovať grafický výstup na strane 44 (dáta sú generované na strane 43) slajdov z prednášky: (http://www.iam.fmph.uniba.sk/institute/stehlikova/cr23/prednasky/cr01_uvod.pdf)

Keďže cykly sú vo všeobecnosti v R pomalé, zvolíme iný prístup. Využijeme funkciu sapply z triedy funkcií apply(s v názve označuje zjednodušenie výstupu, v našom prípade na vektor)

6.1 Príklad použitia funkcie sapply

# sucet 1^2 + 2^2 + ... + k^2
k <- 1:10
sapply(k, FUN = function(n) sum((1:n)^2))
##  [1]   1   5  14  30  55  91 140 204 285 385
# sucet 1^p + 2^p + ... + k^p
k <- 1:10
sapply(k, FUN = function(n, p) sum((1:n)^p), p = 2)
##  [1]   1   5  14  30  55  91 140 204 285 385

6.2 V našom prípade

Použité dáta môžu byť:

  • pevne zvolené (analógia prvého vzorového príkladu)
  • alebo môžu byť parametrom funkcie FUN (analógia druhého vzorového príkladu)

Oplatí sa nastaviť y-ovú os na interval \((0,1)\), aby bol skript univerzálny pre všetky dáta, aby sa vždy (aj pri vysokých p hodnotách) dala vidieť vyznačená hodnota 0.05, s ktorou p hodnoty porovnávame. Tú vieme do existujúceho grafu dokresliť použítím jedného z príkazov

# prvý spôsob
abline( h = 0.05, # horizotálna čiara na hodnote h
        lty = "dashed", # prerušovaná čiara ako typ
        col = "red") # farba

# druhý spôsob
abline( a = 0.05, # intercept
        b = 0, # sklon priamky
        lty = "dashed", # prerušovaná čiara
        col = "red") # farba

Dostávame

7 Analýza výnosov akcií

Pomocou knižnice quantmod načítame priamo do R-ka ceny akcií a zistíme, či sú korelované alebo nie.

7.1 Načítanie dát o cenách akcií pomocou quantmod

Načítajte (ak treba, aj nainštalujte) knižnicu quantmod.

library(quantmod)

Na získanie cien akcií sa použije funkcia getSymbols.

getSymbols(Symbols = skratka pre danú akciu - dá sa nájsť na finance.yahoo, 
           from = dátum odkedy chceme dáta,
           to = dátum dokedy chceme dáta, 
           auto.assign = TRUE/FALSE podľa toho či chceme vytvoriť dataset s daným názvom rovno)

Napríklad:

getSymbols("EBAY", from = "2023-01-01", to = "2023-12-31")
## [1] "EBAY"

Pozrime sa, ako vyzerajú naše dáta, ktoré sú uložené v premennej EBAY, zobrazíme ich začiatok:

head(EBAY)
##            EBAY.Open EBAY.High EBAY.Low EBAY.Close EBAY.Volume EBAY.Adjusted
## 2023-01-03     42.08     42.66    41.54      42.15     4494900      40.56991
## 2023-01-04     42.79     43.34    42.28      43.09     3743500      41.47466
## 2023-01-05     43.01     43.36    42.53      43.10     3748100      41.48429
## 2023-01-06     43.62     45.45    43.06      45.11     5293700      43.41894
## 2023-01-09     44.76     45.44    43.44      43.52     7096400      41.88854
## 2023-01-10     43.32     45.13    43.32      45.09     4432600      43.39969

7.2 Jednoduché ukážky funkcií z balíka

chartSeries(EBAY)

chartSeries(EBAY, theme="white")

Skúste aj:

chartSeries(EBAY, subset="2023-01::2023-06", up.col = "#009E73", dn.col = "#D55E00") # od januara do juna

Funkciou to.monthly transformujeme dáta na mesačné a zobrazíme na grafe.

EBAY.mesacne <- to.monthly(EBAY)  # mesacne data
EBAY.mesacne                      # vypiseme
##          EBAY.Open EBAY.High EBAY.Low EBAY.Close EBAY.Volume EBAY.Adjusted
## jan 2023     42.08     49.89    41.54      49.50    89423000      47.64437
## feb 2023     49.26     52.23    43.71      45.90    95454900      44.17933
## mar 2023     45.68     46.27    40.13      44.37   118451300      42.94928
## apr 2023     44.15     46.67    42.49      46.43    88021200      44.94331
## máj 2023     46.02     46.38    42.01      42.54    99430900      41.41524
## jún 2023     42.80     46.54    42.53      44.69   103529300      43.50839
## júl 2023     44.66     49.48    43.30      44.51   125855200      43.33315
## aug 2023     44.15     45.34    42.23      44.78   104401600      43.84093
## sep 2023     45.21     45.54    42.41      44.09   106480900      43.16541
## okt 2023     43.91     44.04    37.93      39.23   107359500      38.40732
## nov 2023     39.14     41.86    37.17      41.01   164004100      40.39850
## dec 2023     40.99     44.27    40.60      43.62   143449600      42.96959
chartSeries(EBAY.mesacne)         # graf

7.3 Výnosy akcií

Budeme pracovať s výnosmi akcií. Z našich dát EBAY budeme potrebovať posledný stĺpec EBAY.Adjusted, z ktorého vypočítame výnosy. Budeme pracovať zo spojitými výnosmi, teda denné výnosy sa budú počítať ako logaritmus podielu cien v dvoch po sebe idúcich dňoch. Ekvivalentne - ako rozdiel logaritmov (takto sa dá použiť funkcia diff na výpočet diferencií).

ceny <- EBAY$EBAY.Adjusted
vynosy <- diff(log(ceny))
head(vynosy) 
##            EBAY.Adjusted
## 2023-01-03            NA
## 2023-01-04  0.0220560923
## 2023-01-05  0.0002321219
## 2023-01-06  0.0455809586
## 2023-01-09 -0.0358834528
## 2023-01-10  0.0354400235
# vidime, ze prva hodnota je NA, vymazeme ju
vynosy <- vynosy[-1]

Priebeh výnosov:

chartSeries(vynosy, theme=chartTheme('white',up.col="#009E73"))

7.4 Cvičenie: Testovanie autokorelácií výnosov akcií

Zobrazte odhad ACF. Sú autokorelácie signifikantné? Pomocou Ljung-Boxovho testu testujte hypotézu, že výnosy sú nekorelované (teda že sú súčtom konštanty - priemerný výnos - a bieleho šumu).

Výsledok pre porovnanie:

8 AR(1) procesy

V knižnici datasets je viacero zaujímavých časových radov, my sa teraz pozrieme na prietok Nílu:

x <- Nile
x
## Time Series:
## Start = 1871 
## End = 1970 
## Frequency = 1 
##   [1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935 1110  994 1020
##  [16]  960 1180  799  958 1140 1100 1210 1150 1250 1260 1220 1030 1100  774  840
##  [31]  874  694  940  833  701  916  692 1020 1050  969  831  726  456  824  702
##  [46] 1120 1100  832  764  821  768  845  864  862  698  845  744  796 1040  759
##  [61]  781  865  845  944  984  897  822 1010  771  676  649  846  812  742  801
##  [76] 1040  860  874  848  890  744  749  838 1050  918  986  797  923  975  815
##  [91] 1020  906  901 1170  912  746  919  718  714  740
# podrobnejsi popis v helpe: ?Nile

Konkrétne sa zameriame na vyznačenú časť tohto časového radu - od roku 1910 do roku 1960. Aby sme získali časť časového radu, použijeme funkciu window:

x <- window(x, start=1880, end=1920)

Našou úlohou bude zistiť, či je AR(1) proces dobrým modelom pre tieto dáta. Budeme modelovať dáta pomocou funkcie sarima, zatiaľ iba so 4 argumentami.

sarima(xdata = dáta ktoré chceme modelovať,
       p = autoregresný rád,
       d = rád diferencovania,
       q = moving average rád)

8.1 Otázky

  • Je AR(1) proces dobrým modelom pre tieto dáta?
  • Nestačilo by dáta modelovať bez akýchkoľvek autoregresných členov, t. j. ako konštantu plus biely šum?

9 Zapísanie odhadnutého modelu na základe výstupu z R-ka

Vygenerujme si dáta - 100 hodnôt autoregresného procesu \(x_t=0.8x_{t−1}+u_t\), kde \(D[u_t]=2\):

set.seed(123)
x <- arima.sim(model=list(ar=c(0.8)), n=100, sd=sqrt(2))
plot(x)

Odhadneme AR(1) model:

sarima(x,1,0,0, details=FALSE)
## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = xmean, include.mean = FALSE, transform.pars = trans, fixed = fixed, 
##     optim.control = list(trace = trc, REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1   xmean
##       0.7930  0.0204
## s.e.  0.0609  0.5756
## 
## sigma^2 estimated as 1.518:  log likelihood = -163.27,  aic = 332.54
## 
## $degrees_of_freedom
## [1] 98
## 
## $ttable
##       Estimate     SE t.value p.value
## ar1     0.7930 0.0609 13.0132  0.0000
## xmean   0.0204 0.5756  0.0355  0.9717
## 
## $AIC
## [1] 3.325426
## 
## $AICc
## [1] 3.326664
## 
## $BIC
## [1] 3.403582

Odhadnutý model je \[ x_t = \delta + 0.7930 x_{t-1} + u_t \]

kde \(\delta\) je taká konštanta, aby platilo, že \(E[x_t]=0.0204\) (xmean z výstupu).

9.1 Úloha

Dopočítajte konštantu \(\delta\).

10 Teoretický príklad (str. 36 prednáškových slajdov)

(Jeden teoretický príklad tohto typu je v hodnote 10 bodov súčasťou skúšky)

Nech \(z_t\) je proces, ktorého hodnoty sú nezávislé náhodné premenné s rozdelením \(N (0,1)\). Ukážte, že nasledujúci proces je stacionárny a vypočítajte jeho ACF: \[ y_t = \begin{cases} z_t, & \text{pre } t \text{ nepárne}, \\ \frac{1}{\sqrt{2}} \left(z_t^2-1\right), & \text{pre } t \text{ párne}. \end{cases} \]

11 Poznámka: R markdown

R markdown bol použitý aj na vytvorenie tejto stránky priamo v R-ku, takisto na vytvorenie slajdov k prednáške. Užitočné odkazy na tvorbu R markdownu/syntax: