Autoregresné procesy
Beáta Stehlíková, 2-EMM-105 Časové rady, FMFI UK Bratislava
- 1 Opakovanie z prednášky
- 2 Odhadovanie parametrov AR procesu
- 3 Dáta Svetovej banky a grafy pomocou balíka ggplot2
- 4 Overovanie stacionarity AR procesu v R-ku
- 5 Overenie stacionarity odhadnutého AR modelu
- 6 Konštrukcia predikcií
- 7 Autokorelačná funkcia AR procesu
- 8 Parciálna autokorelačná funkcia AR procesu
- 9 Ďalšie príklady
- 10 Nepovinné dodatky
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 nainstalujte2.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/cr17/data/RSQ.txt") # short term rate
r20 <- read.table("http://www.iam.fmph.uniba.sk/institute/stehlikova/cr17/data/R20Q.txt") # 20Y rate
spread <- r20 - rsZ 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
sarima(data,p,0,0, details=FALSE)Výstup:
- odhady parametrov, štandardné odchýlky, informacné kritériá
- graficky rezíduá, ich ACF, p-hodnoty Ljung-Boxovej štatistiky
Úloha: Na prednáške sme si ukázali, že AR(1) proces nie je dobrým modelom pre tieto dáta. Odhadnite ho teraz samostante pomocou uvedenej funkcie a vysvetlite, na základe čoho tento model zamietneme.
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: ?NileKonkré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
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.5097405Odhadnutý 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\).
Cvičenie: Chceme simulovať proces \[x_t = 10 + 0.5 x_t + u_t\] pričom disperzia bieleho šumu je 1. Doplňte hodnoty parametrov k,alpha, sigma.
set.seed(42)
k <-
alpha <-
sigma <-
poc <- 200 # chceme 200 hodnot procesu
x <- arima.sim(model=list(ar=c(alpha)), n=poc, sd=sigma) + kVypočítajte strednú hodnotu tohto procesu a porovnajte ju s odhadnutou hodnotou xmean, ak pre vygenerované dáta odhadneme AR(1) model.
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 Rggplot2- pekné grafy :)
library(WDI)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.23.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')## indicator
## [1,] "GDPPCKD"
## [2,] "GDPPCKN"
## [3,] "NV.AGR.PCAP.KD.ZG"
## [4,] "NY.GDP.PCAP.CD"
## [5,] "NY.GDP.PCAP.KD"
## [6,] "NY.GDP.PCAP.KD.ZG"
## [7,] "NY.GDP.PCAP.KN"
## [8,] "NY.GDP.PCAP.PP.CD"
## [9,] "NY.GDP.PCAP.PP.KD"
## [10,] "NY.GDP.PCAP.PP.KD.ZG"
## [11,] "SE.XPD.PRIM.PC.ZS"
## [12,] "SE.XPD.SECO.PC.ZS"
## [13,] "SE.XPD.TERT.PC.ZS"
## name
## [1,] "GDP per Capita, constant US$, millions"
## [2,] "Real GDP per Capita (real local currency units, various base years)"
## [3,] "Real agricultural GDP per capita growth rate (%)"
## [4,] "GDP per capita (current US$)"
## [5,] "GDP per capita (constant 2000 US$)"
## [6,] "GDP per capita growth (annual %)"
## [7,] "GDP per capita (constant LCU)"
## [8,] "GDP per capita, PPP (current international $)"
## [9,] "GDP per capita, PPP (constant 2005 international $)"
## [10,] "GDP per capita, PPP annual growth (%)"
## [11,] "Expenditure per student, primary (% of GDP per capita)"
## [12,] "Expenditure per student, secondary (% of GDP per capita)"
## [13,] "Expenditure per student, tertiary (% of GDP per capita)"Ak je veľa výsledkov, môžeme chciet vypísať niekoľko prvých:
WDIsearch('gdp.*capita')[1:5,]## indicator
## [1,] "GDPPCKD"
## [2,] "GDPPCKN"
## [3,] "NV.AGR.PCAP.KD.ZG"
## [4,] "NY.GDP.PCAP.CD"
## [5,] "NY.GDP.PCAP.KD"
## name
## [1,] "GDP per Capita, constant US$, millions"
## [2,] "Real GDP per Capita (real local currency units, various base years)"
## [3,] "Real agricultural GDP per capita growth rate (%)"
## [4,] "GDP per capita (current US$)"
## [5,] "GDP per capita (constant 2000 US$)"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)indicatorsme našli pomocouWDIsearchcountryvo formáte iso2c- defaultný
startje 2005, ale mnohé dáta sú dostupné aj skôr (pozrieme si začiatok dát) iso2ckódy štátov sa dajú nájsť napríklad tu: https://datahub.io/core/country-list alebo pomocou knižnicecountrycodev R-ku (https://cran.r-project.org/web/packages/countrycode/)
head(data.hdp)## iso2c country NY.GNP.PCAP.CD year
## 1 CA Canada 47060 2011
## 2 CA Canada 44370 2010
## 3 CA Canada 43110 2009
## 4 CA Canada 44810 2008
## 5 CA Canada 41420 2007
## 6 CA Canada 37780 20063.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)')
datajedata frame, v ktorom sú naše premennéyear, NY.GNP.PCAP.CD- z dát uložených vdatabude na x-ovej osi premennáyeara na y-ovejNY.GNP.PCAP.CDcolor=country- grafy budú odlíšené farebne (pretocolor) podľa premennejcountry
3.5 AR(1) model pre rýchlosť rastu HDP
# GDP per capita (constant 2000 US$)
data <- WDI(indicator='NY.GNP.PCAP.KD',
country=c('US'),
start=1965)
# 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=1965, 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.
- 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 armaRoots. Funkcia armaRoots bola pôvodne v knižnici fArma, táto knižnica však od septembra 2017 nie je dostupná, preto si z jej pôvodnej verzie zoberieme túto funkciu (grafický výstup je vynechaný, kto chce, môže použiť funkciu z armaRoots.R na stránke cvičení, ktorá má parameter doplot s defaultnou hodnotou TRUE):
armaRoots <- function(coefficients, digits = 4)
{
# z povodnej funkcie armaRoots v kniznici fArma
root = polyroot(c(1, -coefficients))
real.root = Re(root)
im.root = Im(root)
distance = sqrt(real.root^2 + im.root^2)
root.mat = cbind(round(real.root, digits = digits), round(im.root,
digits = digits), round(distance, digits = digits))
dimnames(root.mat) = list(1:nrow(root.mat), c("re", "im",
"dist"))
result = root.mat
data.frame(result)
}Ak si funkciu uložíte do R-súboru, tak ju nemusíte kopírovať do každého skriptu, stačí dať tento súbor do pracovaného adresára a potom v skrite načítať pomocou source.
Korene polynómu \(1 - \alpha_1 L - \alpha_2 L^2 - \dots \alpha_p L^p\) teraz nájdeme ako:
armaRoots(c(alfa1, alfa2, ..., alfap))Výstup: tabuľka s koreňmi a ich absolútnymi hodnotami
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: V zápise polynómu \(1-\alpha_1 L - \alpha_2 L^2\) je \(\alpha_1 = 0.9\), \(\alpha_2 = -0.6\), teda:
armaRoots(c(0.9, -0.6))Záver: Korene sú mimo jednotkového kruhu, lebo dist je vždy viac ako 1 \(\Rightarrow\) náš proces JE STACIONÁRNY
## re im dist
## 1 0.75 1.0508 1.291
## 2 0.75 -1.0508 1.2914.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.
Upozornenie: Pri používaní funkcie armaRoots si dajte pozor na znamienka vektora, ktorý do nej vstupuje.
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 funkciuarmaRoots
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.04510079Dá 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.0584042Poč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:
Odhadnite AR(2) model pre premennú
spread. Vysvetlite, prečo rezíduá považujeme za dobré.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Ú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.7120892Ak 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
spreadposledné 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:
Doplňte legendu.
Nepovinné rozšírenie: 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ť.
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 rovnice7.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.0002649014barplot(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-20barplot(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
pacfpomocou funkcie
acf2z knižniceastsavykresliť 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.139 Ď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é
kje len jedna absolútna hodnota (napríkladk = -1.5- modrá), kým pre iné sú dve (k = -0.2- zelená,k = 0.5- hnedá)?Ktoré
kvyhovovujú 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é zistenie 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 NANapíš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:
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 <- leadNá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ť.
Autoregrený proces, ktorého PACF pre lag 3 je nulová.
Autoregrený proces prvého rádu, ktorého PACF pre lag 1 je rovná 0.5.
Autoregresný proces, ktorého ACF je vždy kladná.
Autoregresný proces, ktorého ACF nie je monotónna.
10 Nepovinné dodatky
10.1 Viac o grafoch pomocou ggplot
Z fakultnej siete je na stránke vydavateľstva Springer dostupná kniha ggplot2: Elegant Graphics for Data Analysis od Hadleya Wickhama v pdf formáte: https://link.springer.com/book/10.1007/978-0-387-98141-3
10.2 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.3 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.0584042Z 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 7My 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 5Príslušné koeficienty potom sú
koef.zlozitejsie[grep("^ma",names(koef.zlozitejsie))]## ma1 ma2
## -0.1168011 0.3801231Tutoriál o regulárnych výrazov v R-ku: http://stat545.com/block022_regular-expression.html
10.4 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. týchto 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 NANamiesto 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 FALSETakž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 FALSEViac o funkcii apply a príbuzných funkciách napríklad tu: https://www.datacamp.com/community/tutorials/r-tutorial-apply-family