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 8
packages <- c("astsa", "urca", "datasets", "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 Dáta na samostatnú prácu z Cvičenia 7

Cieľom je nájsť vhodný ARIMA model pre uvedené dáta.

Zadanie úlohy:

  • Vysvetlite, koľkokrát ste dáta diferencovali a prečo. Teda pre každý časový rad (pôvodné dáta, prípadné prvé diferencie, druhé diferencie atď.) napíšte, či ste ich diferencovali a prečo. Skončíme teda tým, že určitý časový rad už diferencovať netreba.

  • Súčasťou predchádzajúceho bodu bolo testovanie jednotkového koreňa. V poslednom kroku (po prípadnom predchádzajúcom diferencovaní) nastala situácia, že v dátach nebol ani trend, ani jednotkový koreň, a preto ich nebolo potrebné diferencovať. Podrobne vysvetlite, čo sa tam dialo:

  1. Napíšte, s akými parametrami ADF testu ste ich tieto dáta testovali a aká regresia sa odhadla.
  2. Aká hypotéza o koeficientoch tejto regresie sa testuje?
  3. Odvoďte, že táto hypotéza predstavuje hypotézu o jedotkovom koreni daného časového radu.
  4. Kedy túto hypotézu zamietame (ako vyzerá kritérium založené na testovacej štatistike a kritickej hodnote)?
  5. Čo vyšlo v našom prípade (zamietame alebo nie) a čo to znamená pre diferencovanie nášho časového radu (diferencujeme alebo nie)?
  • Nájdite vhodný ARIMA model pre dáta y. Požiadavky sú: stacionarita, invertovateľnosť a p-hodnoty Ljung-Boxovho testu nad 5 percent. Odpoveď napíšte zapísaním parametrov funkcie sarima, ktorou model odhadnete.

2.1 Dáta BJsales z knižnice datasets

Popis dát z helpu: The sales time series BJsales and leading indicator BJsales.lead each contain 150 observations.

y <- window(BJsales, end=130)
plot(y)

2.2 Dáta WWWusage z knižnice datasets

Popis dát z helpu: A time series of the numbers of users connected to the Internet through a server every minute.

y <- window(WWWusage, start=1, end=90)
plot(y)

2.3 Dáta faithful z knižnice datasets

Popis dát z helpu: Waiting time between eruptions and the duration of the eruption for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA.

y <- ts(faithful[, 1])
plot(y)

2.4 Dáta blood z knižnice astsa

Popis dát z helpu: Multiple time series of measurements made for 91 days on the three variables, log(white blood count) [WBC], log(platelet) [PLT] and hematocrit [HCT]. Missing data code is NA.

y <- ts(blood[1:36, 3])
plot(y)

2.5 Dáta gnp z knižnice astsa

Popis dát z helpu: Seasonally adjusted quarterly U.S. GNP from 1947(1) to 2002(3).

y <- log(gnp)
plot(y)

3 AR(1) proces pozorovaný s chybou

Vygenerujeme AR(1) proces \(y\):

set.seed(12345)
y <- arima.sim(n=200, list(ar = c(0.9)), sd=0.1)

Pozorovať ho ale budeme s chybou eps, vidíme teda proces \(x\):

eps <- rnorm(200)/15
x <- y+eps

Zobrazte vygenerované procesy \(y\) (pôvodný - čiernou) a \(x\) (pozorovaný - modrou).

Ukážte, že:

  • AR(1) je dobrý model pre dáta \(y\) - to je prirodzené, tak sme ich generovali

  • Ale AR(1) nie je dobrý model pre dáta \(x\) (chyba pri pozorovaní procesu, aj keď nie je veľká, teda spôsobuje problém)

  • ARMA(1,1) je dobrý model pre dáta \(x\).

Poslednú vlastnosť dokážte aj analyticky (nielen pre vygenerované dáta) - teda ak \(y\) je AR(1) a \(x=y+\epsilon\), kde \(\epsilon\) je proces definovaný tak, ako v simuláciách hore, tak \(x\) sa dá zapísať ako ARMA(1,1) proces.

4 Porovnanie Woldových reprezentácií

Budeme pracovať s dátami gas z knižnice astsa. Popis dát z helpu: New York Harbor conventional regular gasoline weekly spot price FOB (in cents per gallon) from 2000 to mid-2010.

y <- gas
plot(y)

Ako prvé ukážte, že aj AR(4) aj MA(3) sú dobré modely pre diferencie premennej gas.

Porovnáme teraz ich Woldove reprezentácie a uvidíme, že tieto dva modely nie sú v skutočnosti až také odlišné.

Definujte

wold <- expand.grid(k=as.factor(1:10), model=c("ar4","ma3"))

Pozrite si data frame wold, aby ste zistili, akú má štruktúru.

Do dataframu wold pridajte stĺpec psi s koeficientami \(\psi_k\) Woldovej reprezentácie pre príslušné modely.

wold$psi <- ...

Porovnanie vykreslite pomocu balíka ggplot2:

ggplot(wold, aes(x=k, y=psi, fill=model)) + 
  geom_bar(stat="identity", position="dodge")

Výstup pre kontrolu:

5 Otázky k DÚ 3?