Autoregresné procesy

1 Opakovanie z prednášky

Zopakujeme otázky to slajdov k príkladu o modelovaní volebných preferencií z prednášky.

2 Odhadovanie parametrov AR procesu

Funkcia v R-ku: Budeme používať funkciu sarima (AR proces je špeciálnym prípadom SARIMA procesu) z knižnice astsa (astsa = applied statistical time series analysis).

library(astsa)  # ak treba, najskor nainstalujte

2.1 Príklad z prednášky: modelovanie spreadu

Zopakujeme najskôr príklad z prednášky - modelovanie spreadu, teda rozdielu medzi dlhodobou a krátkodobou úrokovou mierou

Načítajte dáta a vytvorte premennú spread:

rs <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr20/data/RSQ.txt")   # short term rate
r20 <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr20/data/R20Q.txt") # 20Y rate
spread <- r20 - rs

Z premennej spread spravíme objekt typu time-series, ktorý bude obsahovať informáciu o časovej štruktúre dát.

Použijeme funkciu ts:

spread <- ts(spread,             # data su premennej spread
             frequency=4,        # su kvartalne
             start=c(1952,1))    # zacinaju v 1. kvartali roku 1952

class(spread)  # budeme vidiet, ze uz je to casovy rad

## [1] "ts"

Vykreslíme teraz priebeh - na x-ovej osi budú správne časy:

plot(spread)

Úloha: Zobrazte výberovú autokorelačnú funkciu.

Na prednáške sme videli, že hoci sa podobá na ACF procesu AR(1), nie je to dobrý model pre dáta. Zopakujeme to - odhadneme model a skontrolujeme rezíduá.

Použitie funkcie sarima:

sarima(data, p, k, 0) # AR(p) model pre k-te diferencie premennej data

sarima(data, p, 0, 0)  # AR(p) model pre premennu data
model <- sarima(data,p, 0, 0, details=FALSE)

Výstup pri spustení funkcie sarima:

• konvergencia optimalizačného algoritmu z odhadovania parametrov
• odhady parametrov, štandardné odchýlky, informačné kritériá
• graficky rezíduá, ich ACF, p-hodnoty Ljung-Boxovej štatistiky

Výstup vypísaní objektu model (vytvorený s parametrom details = FALSE):

• odhady parametrov, štandardné odchýlky, informačné kritériá

Úloha:

• Na prednáške sme si ukázali, že AR(1) proces nie je dobrým modelom pre tieto dáta. Odhadnite ho teraz samostatne pomocou uvedenej funkcie a vysvetlite, na základe čoho tento model zamietneme.
• Potom odhadnite a skomentujte AR(2) model.

2.2 Cvičenie: Modelovanie prietoku Nílu

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

library(datasets)
x <- Nile
# 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=1910, end=1960)

Úloha:

• Zistite, či je AR(1) proces dobrým modelom pre tieto dáta
• Odhadnite aj AR(2) a AR(3) modely, zhodnoťte rezíduá. Porovnajte ich informačné kritériá.
• Zistite, či by nestačilo dáta modelovať bez akýchkoľvek autoregresných členov, t. j. ako konštantu plus biely šum.

2.3 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.8 x_{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:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xmean, include.mean = FALSE, 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] 1.457637
## 
## $AICc
## [1] 1.480137
## 
## $BIC
## [1] 0.5097405

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

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

3 Dáta Svetovej banky a grafy pomocou balíka ggplot2

Obsah tejto časti:

• prístup k dátam Svetovej banky priamo v R-ku
• príklad použitia knižnice ggplot2 na zostrojenie pekných grafov
• pripravíme si dáta, ktoré potom použijeme na modelovanie AR procesom

3.1 Dáta Svetovej banky

Nainštalujte a načítajte knižnice

• WDI - World Development Indicators, prístup k dátam budeme mať priamo v R
• ggplot2 - pekné grafy :)

library(WDI)
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.4.2

3.2 Vyhľadávanie dát vo WDI

# hladame indikatory, ktore maju v nazve gdp
WDIsearch('gdp') 

Kto pozná regulárne výrazy (regular expressions)

• dajú sa používať aj tie
• na ukážku - hľadajme indikátory, ktoré hovoria o HDP na obyvateľa - vypíšme si tie, ktoré obsahujú gdp a capita

WDIsearch('gdp.*capita')

Ak je veľa výsledkov, môžeme chciet vypísať niekoľko prvých:

WDIsearch('gdp.*capita')[1:3,]

##      indicator          
## [1,] "6.0.GDPpc_constant"          
## [2,] "FB.DPT.INSU.PC.ZS"          
## [3,] "NE.GDI.FTOT.CR"
   
##      name                                                                 
## [1,] "GDP per capita, PPP (constant 2011 international $) "                             
## [2,] "Deposit insurance coverage (% of GDP per capita)"
## [3,] "GDP expenditure on gross fixed capital formation (in IDR Million)" 

3.3 Načítanie dát

Pomocou funkcie WDI, ukážka použitia:

data.hdp <- WDI(indicator = 'NY.GNP.PCAP.CD',
                country = c('CA','US','FR', 'DE'),
                start = 1975)

• indicator sme našli pomocou WDIsearch
• country vo formáte iso2c
• defaultný start je 1960 a defaultný end je 2020, môžeme ich takto zmeniť
• iso2c kódy štátov sa dajú nájsť napríklad tu: https://datahub.io/core/country-list alebo pomocou knižnice countrycode v R-ku (https://cran.r-project.org/web/packages/countrycode/)

head(data.hdp)

##   iso2c country NY.GNP.PCAP.CD year
## 1    CA  Canada             NA 2020
## 2    CA  Canada          46370 2019
## 3    CA  Canada          43500 2018
## 4    CA  Canada          42990 2017
## 5    CA  Canada          43870 2016
## 6    CA  Canada          47570 2015

3.4 Graf pomocou ggplot

ggplot(data.hdp, aes(year, NY.GNP.PCAP.CD, color = country)) +
    geom_line() +
    xlab('Year') + ylab('GDP per capita') +
    labs(title = 'GDP per capita (current USD)')

• data.hdp je data frame, v ktorom sú naše premenné
• year, NY.GNP.PCAP.CD - z dát uložených v data.hdp bude na x-ovej osi premenná year a na y-ovej NY.GNP.PCAP.CD
• color=country - grafy budú odlíšené farebne (preto color) podľa premennej country

3.5 AR(1) model pre rýchlosť rastu HDP

# GDP per capita (constant 2010 US$) 
data <- WDI(indicator = 'NY.GNP.PCAP.KD',
            country= c ('FR'),
            start = 1975, end = 2019)

# zoradime rastuco podla rokov
data <- data[order(data$year),]

# nasa premenna na modelovanie
log.y <- log(data$NY.GNP.PCAP.KD)

# doplnime casovu strukturu
log.y <- ts(log.y, start = 1975, frequency = 1)

Zadanie:

• Z priebehu dát vidíme, ze nie sú stacionárne, pretovbudeme pracovat s diferenciami
• Zoberieme diferencie logaritmov - to je rýchlost rastu HDP.
• Vyskúšajte, či je AR(1) model dobrým modelom pre rýchlosť rastu HDP.
• Ak je AR(1) model dobrý, zistite, či je autoregresný člen nutný a či sa rýchlosť rastu nedá modelovať ako posunutý biely šum. Ak AR(1) nevyhovuje, vyskúšajte AR modely vyššieho rádu.
• Zopakujte pre iný štát.

4 Overovanie stacionarity AR procesu v R-ku

V prípade AR(1) procesu, je overovanie stacionarity zrejmé - stačí sa pozrieť na absolútnu hodnotu autoregresného koeficientu, musí byť v absolútnej hodnote menší ako 1.

Pre všeobecný AR(p) proces je v prednáškach odvodená podmienka stacionarity: proces \((1 - \alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p)x_t = \delta + u_t\) je stacionárny práve vtedy, keď sú všetky korene polynómu \(1 - \alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p\) v absolútnej hodnote väčšie ako 1 (geometricky: mimo jednotkového kruhu).

Otázka: Ukážte, že toto kritériu je v súlade s tým, čo je povedané pre AR(1) proces o jeho autoregresnom koeficiente.

Funkcia v R-ku: Korene polynómu \(1 - \alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p\) nájdeme pomocou funkcie polyroot, ktorá je vysvetlená v slajdoch z prenášky.

4.1 Príklad: Overenie stacionarity AR(2) procesu

Zadanie: Zistíme, či je proces \((1 - 0.9 L + 0.6 L^2)x_t = u_t\) stacionárny.

Výpočet v R-ku: Hľadáme korene (a ich absolútne hodnoty) polynómu \(1- 0.9 L + 0.6 L^2\), teda:

polyroot(c(1, -0.9, 0.6))
abs(polyroot(c(1, -0.9, 0.6)))

Záver: Korene sú mimo jednotkového kruhu, lebo absoútne hodnota je pre každý koreň väčšia ako 1 \(\Rightarrow\) náš proces JE STACIONÁRNY

Cvičenie: Zistite, či sú stacionárne nasledovné AR(2) procesy:

• \((1 - 0.25 L + 0.8 L^2)x_t = -3 + u_t\)
• \(x_t = 2 + 0.8 x_{t-1} - 0.1 x_{t-2} + u_t\)

4.2 Cvičenie: Overenie stacionarity AR procesu (súčasť kostry na skúške)

Zistite, či sú stacionárne nasledovné AR procesy:

• \((1 - 0.25 L + 0.6 L^2 - 0.55 L^3)x_t = u_t\)
• \((1 + 0.3 L + 0.2 L^2)x_t = 5 + u_t\)
• \(x_t = 2 + 0.3 x_{t-1} - 0.3 x_{t-2} + u_t\)
• \(x_t = 0.25 + 0.1 x_{t-1} + 0.2 x_{t-3} + u_t\)

Vypočítajte strednú hodnotu stacionárnych procesov.

4.3 Cvičenie: Nájdenie príkladu procesu s danou vlastnosťou

Nájdite príklad

• stacionárneho AR(3) procesu,
• nestacionárneho AR(4) procesu,
• stacionárneho AR(2) procesu, ktorého stredná hodnota je rovná 10.

V každom z týchto prípadov dokážte, že váš proces naozaj má požadovanú vlastnosť.

Príklady tohto typu “nájdite príklad procesu, ktorý … a ukážte, že má požadovanú vlastnosť” budú aj na skúške.

5 Overenie stacionarity odhadnutého AR modelu

5.1 Pristupovanie ku koeficientom modelu

Samozrejme pri AR(1) modeli sme hneď videli, či je stacionárny alebo nie, pri modeloch vyššieho rádu použijeme funkciupolyroot

Hodnoty koeficientov však nemusíme odpisovať, tak by sme navyše dostali iba ich približné hodnoty

Napríklad ak odhadneme model:

model.ar3 <- sarima(spread, 3,0,0, details=FALSE)

Pozrite si str(model.ar3) - štruktúru objektu model.ar3 - a potom vyberte z model.ar3 vektor odhadnutých parametrov modelu.

Z tohto chceme tie zložky vektora, ktorých názvy obsahujú ar:

##         ar1         ar2         ar3       xmean 
##  1.18493093 -0.30509964  0.01391889  1.04510079

Dá sa spraviť koef[1:p], ak je to AR(p) model

Budeme mať aj všeobecnejšie modely (koeficienty budú s názvami ma1, ma2, … neskôr aj napríklad sar1, sma1), napríklad:

koef.zlozitejsie

##        ar1        ar2        ar3        ma1        ma2       sar1 
##  1.3167189 -0.7961300  0.3590080 -0.1168011  0.3801231  0.9967049 
##       sma1      xmean 
## -0.9797410  1.0584042

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ď. Postup je v nepovinných dodatkoch.

5.2 Cvičenie: AR(2) model pre spread, overenie stacionarity

Na prednáške sme videli, že ak pre dáta v premennej spread odhadneme AR(2) model, bude mať dobré rezíduá.

Úloha: Overte stacionaritu získaného modelu.

6 Konštrukcia predikcií

Cieľom modelovania je často konštrukcia predikcií do budúcnosti.

Funkcia v R-ku: 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) 

Úloha: Pomocou AR(2) modelu pre dáta spreadspravte predikcie na nasledujúce dva roky (máme kvartálne dáta, takže 2 roky predstavujú 8 pozorovaní)

## $pred
##             Qtr1        Qtr2        Qtr3        Qtr4
## 2006 -0.11864092  0.03883299  0.19265253  0.32884850
## 2007  0.44528788  0.54348324  0.62583626  0.69474654
## 
## $se
##           Qtr1      Qtr2      Qtr3      Qtr4
## 2006 0.6838752 1.0582313 1.3006955 1.4585690
## 2007 1.5627848 1.6325528 1.6798046 1.7120892

Ak uložíme výstup funkcie sarima.for napr. do premennej predikcie, tak máme prístup k predikciám a štandardným odchýlkam: predikcie$pred, predikcie$se. To je užitočné napríklad vtedy, keď chceme vynechať nejaké posledné dáta, odhadnúť model pomocou zostávajúcich a potom porovnať predikcie s tým, čo v skutočnosti nastalo.

Úloha:

• Vynechajte z dát v premennej spread posledné dva roky.

• Zo zostávajúcich dát odhadnite AR(2) model. Overte jeho stacionaritu a či môžeme považovať rezíduá za biely šum.

• Spravte z tohto modelu predikcie na nasledujúce dva roky.

• Nakreslite do jedného grafu farebne odlíšené: dáta použité na odhadovanie modelu, predikcie, predikcie +/- 2 \('\times\) štandardná odchýlka, skutočný vývoj v posledných dvoch rokoch. Na to sa dá dobre využiť funkcia ts.plot, môžeme spraviť napríklad:

 ts.plot(data1, data2, data3,
        gpars=list(xlab="time", ylab="popis dat",
                  col=c("red","blue","black"))
)

Pritom jednotlivé dáta nemusia mať rovnaký rozsah časov, pre ktoré sú dostupné, R-ko sa postará aj o rozsah y-ovej osi (čo nie je zaručené pri kreslení pomocou plot a potom lines).

Výstup (defaultne sa kreslí posledných 100 pozorovaní + predikcie, parametrom plot.all = TRUE funkcie sarima.for môžeme dosiahnuť vykreslenie všetkých dát):

Doplňte legendu.

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.

Nepovinné rozšírenie: Zvoľte si vlastný rozsah použitých dát, ktoré chcete zobraziť (teda nie nutne 100, resp. všetky). Namiesto čiar, ktoré ohraničujú predikčné intervaly spravte vyfarbenú oblasť:

Návod: https://stats.stackexchange.com/questions/154346/fitted-confidence-intervals-forecast-function-r

Samozrejme, graf môžete ďalej vylepšovať.

Cvičenie: Modelovali sme rýchlosť rastu HDP autoregresným procesom. Vynechajte teraz z dát posledných 5 rokov. Odhadnite vhodný model a použitím nájdeného modelu spravte predikcie pre:

• logaritmus HDP
• rýchlosť rastu HDP

Spravte porovnanie s reálnymi hodnotami týchto premenných.

7 Autokorelačná funkcia AR procesu

Nájdeme autokorelačnú funkciu AR(3) procesu \[x_t = 1.5 x_{t-1} - 0.8 x_{t-2} + 0.2 x_{t-3} + u_t,\] a to dvoma spôsobmi:

• priamo z Yule-Wolkerových rovníc a diferenčnej rovnice z prednášky

• pomocou funkcie v R-ku

7.1 Z Yule Wolkerových rovníc

Postup:

• Vyriešime v R-ku Yule-Wolkerove rovnice - dostaneme ACF(k) pre k=1,2,3 (pozrite ?matrix, ?solve)

• V cykle vypočítame nasledjúce hodnoty ACF

rho <- rep(0, times=10)
rho <- ....                   # prve tri zlozky ako riesenie sustavy
for (i in 4:10) rho[i] <- ... # z diferencnej rovnice

7.2 Funkcia v R-ku

Použite funkciu ARMAacf pre náš proces

ARMAacf(ar=c(...), lag.max=10)

a porovnajte výsledky z oboch postupov. Funkcia ARMAacf by už nemala byť “čiernou skrinkou” ;-)

8 Parciálna autokorelačná funkcia AR procesu

8.1 Opakovanie z prednášky

• Čo je parciálna autokorelačná funkcia procesu?

• Čím je charakteristická PACF autoregresného procesu rádu p? Prečo? Ako to vyplýva z definície PACF?

8.2 Výpočet v R-ku

Nájdeme autokorelačnú funkciu AR(3) procesu \[x_t = 1.5 x_{t-1} - 0.8 x_{t-2} + 0.2 x_{t-3} + u_t,\] Podobne ako v prípade ACF, aj teraz použijeme funkciu ARMAacf, ale pridáme parameter pacf=TRUE:

ARMAacf(ar=c(...), lag.max=10, pacf=TRUE)

8.3 Kontrolná otázka

ACF aj PACF sa v R-ku počítajú numericky.

Ktoré z nasledujúcich postupností

• sa po určitom počte členov vynulujú (a teda tie malé čísla sú v skutočnosti nuly)

a ktoré

• majú hodnoty, ktoré konvergujú k nule, ale neexistuje index, od ktorého by boli nulové?

Príklad 1

ARMAacf(ar=c(0.05,-0.2), lag.max=10)

##             0             1             2             3             4 
##  1.0000000000  0.0416666667 -0.1979166667 -0.0182291667  0.0386718750 
##             5             6             7             8             9 
##  0.0055794271 -0.0074554036 -0.0014886556  0.0014166479  0.0003685635 
##            10 
## -0.0002649014

barplot(ARMAacf(ar=c(0.05,-0.2), lag.max=10))

Príklad 2

ARMAacf(ar=c(0.3,0.05), lag.max=10, pacf=TRUE)

##  [1]  3.157895e-01  5.000000e-02 -3.863413e-18  3.927096e-19 -4.622862e-19
##  [6] -1.450362e-19 -4.829267e-20  6.640242e-20  5.131096e-20 -2.263719e-20

barplot(ARMAacf(ar=c(0.3,0.05), lag.max=10,pacf=TRUE))

8.4 Výberová PACF z dát

Môžeme

• použiť funkciu pacf

• pomocou funkcie acf2 z knižnice astsa vykresliť naraz výberovú ACF aj PACF (tu bude mať navyše ACF odstránenú hodnotu pre lag 0, ktorá sa vždy rovná 1)

x <- rnorm(100)

# moznost 1
pacf(x)

# moznost 2
library(astsa)
acf2(x)

##         ACF  PACF
##  [1,] -0.09 -0.09
##  [2,] -0.13 -0.14
##  [3,]  0.04  0.01
##  [4,] -0.20 -0.22
##  [5,]  0.08  0.05
##  [6,] -0.18 -0.25
##  [7,]  0.02  0.01
##  [8,]  0.15  0.03
##  [9,] -0.06 -0.01
## [10,]  0.17  0.12
## [11,] -0.13 -0.11
## [12,] -0.03  0.02
## [13,]  0.14  0.08
## [14,] -0.10  0.03
## [15,]  0.15  0.13
## [16,] -0.12 -0.09
## [17,] -0.08 -0.04
## [18,]  0.02 -0.10
## [19,] -0.06  0.04
## [20,] -0.01 -0.13

9 Ďalšie príklady

9.1 Stacionarita a parameter procesu

Zadanie: Urcte všetky hodnoty parametra k, pre ktoré je stacionárny AR(2) proces \(x_t = x_{t-1} + k x_{t-2} + u_t\).

Skúsime najskôr numericky: Pre nejaký rozsah parametra k nájdeme absolútnu hodnotu koreňov (pomocou armaRoots).

Otázky k nasledujúcemu obrázku:

• Preco pre niektoré k je len jedna absolútna hodnota (napríklad k = -1.5 - modrá), kým pre iné sú dve (k = -0.2 - zelená, k = 0.5 - hnedá)?

• Ktoré k vyhovovujú podmienke stacionarity a ktoré nie?

Analytický výpočet:

• Odvoďte korene a ich absolútne hodnoty.
• Zistite, kedy je splnená podmienka stacionarity.
• Skontrolujte si, či sa vaše výpočty zhodujú s numerickými na predchádzajúcom obrázku

Precvičenie programovania v R-ku: Spravte samostatne graf z predchádzajúceho obrázku.

9.2 Príklad: Zmena príkladu kvôli preklepu

V jednej z domácich úloh je úlohou nájsť príklady z kurzov časových radov na zahraničných univerzitách a samostatne ich vyriešiť. Takto sa našlo aj nasledovné zadanie:

Zistite, pre ktoré hodnoty parametra \(k\) je stacionárny proces \[x_t = x_{t-1} + k x_{t-2} + k x_{t-3} + u_{t}.\]

Zo vzorového riešenia vyplýva, že malo ísť o proces \[x_t = x_{t-1} + k x_{t-2} - k x_{t-3} + u_{t}.\]

Vyriešte obidve verzie zadania (preklepom sa úloha stala zložitejšou, ale podmienka stacionarity sa dá odvodiť analyticky aj pre prvý proces).

9.3 Simulačné a analytické zistenie podmienky stacionarity pre AR(2) proces

Uvažujme proces \[x_t = \delta + \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + u_t.\] Má dva paramametre, takže ich môžeme znázorniť do roviny a odlíšiť tie, ktoré zodpovedajú stacionárnemu a tie, ktoré zodpovedajú nestacionárnemu AR(2) procesu.

n <- 2000          # najskor menej, po odladeni napr. týchto 2000
set.seed(123)
alpha1 <- runif(n, min = -2.5, max = 2.5) # rovnomerne na (min, max)
alpha2 <- runif(n, min = -2.5, max = 2.5)
stationary <- rep(NA,n)
df <- data.frame(alpha1, alpha2, stationary)
head(df)

##       alpha1      alpha2 stationary
## 1 -1.0621124 -1.70163010         NA
## 2  1.4415257 -1.77742074         NA
## 3 -0.4551154 -1.75409804         NA
## 4  1.9150870  0.07217130         NA
## 5  2.2023364 -0.03586347         NA
## 6 -2.2722175  0.58171384         NA

Napíšte funkciu, ktorá rozhodne, či je AR(2) proces so zadanými parametrami stacionárny. Napríklad stat.ar.2(c(1,1)) vráti FALSE.

stat.ar.2 <- function(koef) # doplnte definiciu funkcie

for (i in 1:n) {
  # doplnte cyklus
  # TRUE, FALSE podla toho, ci je prislusny proces stacionarny
  df$stationary[i] <-
}

Potom môžeme pomocou knižnice `gplot2 kresliť:

qplot(alpha1, alpha2, colour=stationary, data=df)

Aby bol obrázok výstižný, treba viac vygenerovaných bodov (napríklad tých 2000 spomínaných hore), ale na nasledujúcom grafe pre 100 bodov vidíme, aký typ grafu dostaneme:

Sformulujte tvrdenie o podmienke stacionarity AR(2) procesu (dá sa vyjadriť pomocou troch nerovností pre parametre, ktoré musia súčasne platiť) a dokážte ho. Z grafu môžeme získať hypotézu, ako vyzerá podmienka stacionarity pre AR(2) proces a následne sa ju snažiť dokázať. Alebo ho môžeme použiť ako skúšku správnosti našich výpočtov.

9.4 Hľadanie vhodného autoregresného modelu pre zadané dáta

Všetky nasledujúce zadania sú zo starých skúšok. Úlohou je

• nájsť vhodný autoregresný model pre zadané dáta

• vysvetliť, prečo sú rezíduá vyhovujúce - skomentovať 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í)

Na skúške je súčasťou kostry zistenie, či môžeme pracovať priamo so zadanými dátami alebo s diferenciami (trend nie je jediným dôvodom diferencovannia - budeme sa tomu venovať neskôr). V týchto zadaniach je povedané, s akými dátami máte pri hľadaní AR modelu pracovať.

Takisto je na skúške na výber širšia trieda modelov, ako sú AR modely. V niektorých prípadoch - vrátane týchto - však AR model stačí.

9.4.1 Ceny benzínu

Dáta z knižnice astsa, z popisu v helpe: New York Harbor conventional regular gasoline weekly spot price FOB (in cents per gallon) from 2000 to mid-2010. Zoberieme dáta od roku 2006.

library(astsa)
data(gas)
x <- window(gas, start=c(2006,1))

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

9.4.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)
data(lead)
x <- lead

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

9.5 Nájdenie procesu s danou vlastnosťou

Nájdite príklad procesov, ktoré majú nasledovné vlastnosti. Pre každý proces dokážte, že má naozaj požadovanú vlastnosť.

1. Autoregrený proces, ktorého PACF pre lag 3 je nulová.

2. Autoregrený proces prvého rádu, ktorého PACF pre lag 1 je rovná 0.5.

3. Autoregresný proces, ktorého ACF je vždy kladná.

4. Autoregresný proces, ktorého ACF nie je monotónna.

10 Nepovinné dodatky

10.1 Iné pekné grafy pomocou dygraphs

Napríklad zvýraznená časť grafu s priebehom časového radu, keď sme modelovali prietok Nílu, alebo takýto interaktívny graf (s tými istými dátami):

x <- Nile
library(dygraphs)
dygraph(x)  %>% dySeries(label="Annual flow", color="black") %>% dyRangeSelector()

O knižnici dygraphs: https://rstudio.github.io/dygraphs/index.html

10.2 Regulárne výrazy pre pristupovanie k parametrov odhadnutého modelu

Vráťme sa k odhadnutému AR modelu z časti 4 (sarima(spread,3,0,2,2,0,1))s nasledovnými koeficientami:

koef.zlozitejsie

##        ar1        ar2        ar3        ma1        ma2       sar1 
##  1.3167189 -0.7961300  0.3590080 -0.1168011  0.3801231  0.9967049 
##       sma1      xmean 
## -0.9797410  1.0584042

Z nich chceme vybrať tie, ktoré sa začínajú na ma.

Názvy zložiek vektora sú

names(koef.zlozitejsie)

## [1] "ar1"   "ar2"   "ar3"   "ma1"   "ma2"   "sar1"  "sma1"  "xmean"

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 checeme iba tie, ktoré začínajú ma (lebo je tam aj sma1 koeficient, ktorý necheme). 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.1168011  0.3801231

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

10.3 Trieda funkcií apply

Vráťme sa k príkladu, kde sme pre zadané kombinácie parametrov \(\alpha_1, \alpha_2\) zisťovali, či je AR(2) proces s týmito koeficientami stacionárny.

n <- 10          # najskor menej, po odladeni napr. 2000
set.seed(12345)
alpha1 <- runif(n, min = -2.5, max = 2.5) # rovnomerne na (min, max)
alpha2 <- runif(n, min = -2.5, max = 2.5)
stationary <- rep(NA, n)
df <- data.frame(alpha1, alpha2, stationary)
df

##         alpha1     alpha2 stationary
## 1   1.10451948 -2.3273228         NA
## 2   1.87886597 -1.7381325         NA
## 3   1.30491164  1.1784248         NA
## 4   1.93062283 -2.4943171         NA
## 5  -0.21759520 -0.5439833         NA
## 6  -1.66814107 -0.1875267         NA
## 7  -0.87452307 -0.5592801         NA
## 8   0.04612168 -0.4875743         NA
## 9   1.13852627 -1.6051821         NA
## 10  2.44868469  2.2582938         NA

Namiesto toho, aby sme hodnoty stĺpca stationaryvkladali v cykle, môžeme použiť funkciu apply.

# Predpokladame fungujucu funkciu stac.ar.2, ktora ma ako vstup dvojzlozkovy vektor a vystup TRUE/FALSE

apply(df[, 1:2],        # vstupne data = 1. a 2. stlpec df
      MARGIN = 1,        # znamena, ze sa budu brat po riadkoch
      FUN = stac.ar.2)   # aplikuje sa funkcia stac.ar.2

##  [1] FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE

Takže sa dá spraviť

df$stationary <- apply(df[, 1:2], MARGIN = 1,FUN = stac.ar.2)
df

##         alpha1     alpha2 stationary
## 1   1.10451948 -2.3273228      FALSE
## 2   1.87886597 -1.7381325      FALSE
## 3   1.30491164  1.1784248      FALSE
## 4   1.93062283 -2.4943171      FALSE
## 5  -0.21759520 -0.5439833       TRUE
## 6  -1.66814107 -0.1875267      FALSE
## 7  -0.87452307 -0.5592801       TRUE
## 8   0.04612168 -0.4875743       TRUE
## 9   1.13852627 -1.6051821      FALSE
## 10  2.44868469  2.2582938      FALSE

Viac o funkcii apply a príbuzných funkciách napríklad tu: https://www.datacamp.com/community/tutorials/r-tutorial-apply-family