Cvičenie 4: ARMA procesy, testovanie jednotkového koreňa
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 4
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 Autokorelačná funkcia AR procesu
Nájdeme autokorelačnú funkciu AR(3) procesu \[x_t=1.5x_{t−1}−0.8x_{t−2}+0.2x_{t−3}+u_t\]
Ukážeme si dva spôsoby riešenia:
- priamo z Yule-Wolkerových rovníc a diferenčnej rovnice z prednášky,
- pomocou funkcie v R-ku.
2.1 Z Yule Wolkerových rovníc
Pri tomto spôsobe riešenia je prvým krokom zostavenie
Yule-Wolkerových rovníc. Tie následne vyriešime v R-ku. Dostaneme ACF(k)
pre \(k=1,2,3\) (pozrite
?matrix, ?solve). Následne v cykle z
diferenčnej rovnice vypočítame ďalšie hodnoty ACF.
rho <- rep(0, times=10)
rho <- .... # prve tri zlozky ako riesenie sustavy
for (i in 4:19) rho[i] <- ... # ostatne zlozky z diferencnej rovniceVýstup pre kontrolu:
## [1] 0.9178082 0.7602740 0.6061644 0.4845890 0.3940068 0.3245719 0.2685702
## [8] 0.2219991 0.1830569 0.15070012.2 Funkcia v R-ku
Použijeme funkciu ARMAacf pre náš proces.
ARMAacf(ar=c(...), lag.max=10)Výstup pre kontrolu:
## 0 1 2 3 4 5 6 7
## 1.0000000 0.9178082 0.7602740 0.6061644 0.4845890 0.3940068 0.3245719 0.2685702
## 8 9 10
## 0.2219991 0.1830569 0.1507001Porovnajte výsledky z oboch postupov. Funkcia ARMAacf by
už nemala byť “čiernou skrinkou” :).
3 Parciálna autokorelačná funkcia AR procesu
3.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?
3.2 Výpočet v R-ku
Nájdeme parciálnu autokorelačnú funkciu AR(3) procesu \[x_t=1.5x_{t−1}−0.8x_{t−2}+0.2x_{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)## [1] 9.178082e-01 -5.208333e-01 2.000000e-01 0.000000e+00 5.033930e-16
## [6] -1.123254e-15 8.351739e-16 9.367857e-17 -6.465462e-16 5.159779e-163.3 Výberová PACF z dát
Môžeme
- použiť funkciu
pacf, - pomocou 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).
Napríklad:
x <- rnorm(100)
pacf(x) # moznost 1
acf2(x) # moznost 2
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## ACF -0.13 -0.14 0.03 0.01 0.12 -0.14 0.12 -0.02 -0.01 0.07 -0.02 0.01
## PACF -0.13 -0.16 -0.01 -0.01 0.13 -0.11 0.13 -0.03 0.03 0.05 0.03 -0.01
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## ACF -0.01 0.01 0.04 0.14 -0.05 -0.07 0.03 -0.07
## PACF 0.03 -0.02 0.04 0.17 -0.02 -0.03 0.00 -0.113.4 ACF, PACF: 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))
4 Jednotkový koreň
Cieľom tejto časti bude zisťovanie rádu diferencovania časového radu.
Na predošlých cvičeniach sme ukázali, že časový rad môže obsahovať deterministický trend (napríklad lineárny v podobe \(\mu t\)), ktorý spôsobuje nestacionaritu. Vo väčšine prípadov pozorovaný trend vieme aproximovať polýnomom v premennej \(t\) rádu \(p\). Následne pri ARMA modelovaní pracujeme s \(p\)-tými diferenciami časového radu.
Časový rad však môže vykazovať aj iný typ nestacionarity, spôsobenej takzvaným jednotkovým koreňom. Procesy s jednotkovým koreňom nazývame integrované procesy a označujeme ich ako \(I(d)\), kde \(d\) predstavuje počet diferenciácií potrebných na dosiahnutie stacionarity. To znamená, že pomocou \(d\)-tych diferencií vieme dostať stacionárny a invertovateľný proces z pôvodného procesu. Vo všeobecnosti, hovoríme, že \(y_t\) je \(I(d)\), ak ho vieme zapísať ako \[(1−L)^dy_t=\delta+x_t,\]
pričom \(x_t\) je stacionárny a invertovateľný ARMA(p,q) proces.
Na odhalenie prítomnosti jednotkového koreňa v časovom rade budeme
používať štatistické testy z balíka urca (Unit Root and
Cointegration Tests for Time Series Data).
Konkrétne funckiu ur.df (ur = unit root, df =
Dickey-Fuller):
ur.df(y = časový rad,
type = jedno z c(
"trend" (konštanta aj lineárny trend),
"drift" (konštanta bez lineárneho trendu) ,
"none" (nič)),
lags = maximálny počet lagov, ktorý uvažujeme,
selectlags = kritérium, podľa ktorého sa vyberá počet lagov - jedno z c("Fixed", "AIC", "BIC"))4.1 Modelový príklad z minimálnej kostry na skúške
Budeme pracovať s dátami o nezamestnanosti, ktoré sú dostupné v
balíku astsa:
y <- econ5[,"unemp"]
plot(y)
Cieľom je nájsť vhodný ARIMA modelu pre uvedené dáta. (Na skúške bývajú dáta podobné týmto, ktoré sú vždy vybrané tak, aby sa pre ne dal nájsť nesezónny ARIMA(p,d,q) model, pričom p ani q nie sú väčšie ako 5).
Zadanie úlohy bude v nasledovnom tvare:
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:
- Napíšte, s akými parametrami ADF testu ste ich tieto dáta testovali a aká regresia sa odhadla.
- Aká hypotéza o koeficientoch tejto regresie sa testuje?
- Odvoďte, že táto hypotéza predstavuje hypotézu o jedotkovom koreni daného časového radu.
- Kedy túto hypotézu zamietame (ako vyzerá kritérium založené na testovacej štatistike a kritickej hodnote)?
- Č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.
Vzorové riešenie:
- Nami skúmané dáta podľa grafu nevykazujú prítomnosť trendu. To
znamená, že dáta
yotestujeme už iba na prítomnosť jednotkového koreňa.
- Stredná hodnota dát vyzerá byť niekde v okolí hodnoty 5, preto
zvolíme
type = "drift". Je dôležité zvoliť primeraný počet lagov, aby sa nám nestalo, že všetky lagy sa nám dostanú do výsledného modelu. V takom prípade nevieme, či je to najlepší možný model pre naše dáta podľa informačného kritéria, alebo sme ho dostali iba kvôli tomu, že sme väčší počet lagov nepovolili. V našom prípade zvolímelags = 4. Na výber počtu lagov zvolíme Bayesovo informačné kritérium, t. j.selectlags = "BIC".
summary(ur.df(y, lags = 5, selectlags = "BIC", type = "drift"))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82251 -0.22532 -0.02061 0.24447 1.89015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.34661 0.12149 2.853 0.00493 **
## z.lag.1 -0.06217 0.02033 -3.059 0.00263 **
## z.diff.lag 0.47257 0.06930 6.819 2.03e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4275 on 152 degrees of freedom
## Multiple R-squared: 0.2507, Adjusted R-squared: 0.2408
## F-statistic: 25.42 on 2 and 152 DF, p-value: 2.987e-10
##
##
## Value of test-statistic is: -3.0585 4.7181
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Regresia, ktorá sa v teste odhadla je: \[ \Delta z_t = \alpha + c_1 z_{t-1} + c_2 \Delta z_{t-1} + \epsilon_t\]
- Testuje sa nasledovná hypotéza o koeficientoch tejto regresie:
\[ H_0: c_1=0 \qquad \text{vs.} \qquad H_1: c_1<0\]
- Odvodenie, že táto hypotéza predstavuje hypotézu o jedotkovom koreni daného časového radu:
\[ \begin{align} \Delta z_t &=\alpha + c_1 z_{t-1} + c_2 \Delta z_{t-1} + \epsilon_t, \\ (z_t - z_{t-1}) &=\alpha + c_1 z_{t-1} + c_2 (z_{t-1} - z_{t-2}) + \epsilon_t, \\ z_t &=\alpha + (1+c_1+c_2) z_{t-1} + (-c_2) z_{t-2} + \epsilon_t, \\ [1 - (1+c_1+c_2)L + c_2 L^2] z_t &= \alpha + \epsilon_t. \end{align} \]
Ak by polynóm \(1 - (1+c_1+c_2)L + c_2 L^2\) mal jednotkový koreň, platilo by \[ \begin{align} 1 - (1+c_1+c_2) + c_2 &= 0, \\ - (1+c_1+c_2) + c_2 &= -1, \\ c_1 &= 0. \end{align}\]
Teda vidíme, že hypotéza o testovaní koeficientu \(c_1\) z regresie, je totožná s hypotézou o testovaní jednotkového koreňa.
Hypotézu zamietame, ak je testová štatistika menšia ako kritická hodnota.
V našom prípade vyšlo, že testová štatistika (
Value of test-statistic is: -3.4904) je na hladine významnosti 5% menšia ako kritická hodnota (Critical values for test statistics: (5pct, tau2) -2.88). Preto hypotézu \(H_0\) zamietame. Znamená to, že v dátach nie je prítomný jednotkový koreň, teda dáta nediferencujeme.
- Teraz je úlohou nájsť vhodný ARIMA model. Na to sa nám zíde zobraziť výberovú ACF a PACF.
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF 0.95 0.86 0.76 0.66 0.59 0.54 0.51 0.48 0.46 0.44 0.41 0.38 0.36
## PACF 0.95 -0.48 -0.01 0.12 0.16 -0.03 0.05 -0.10 0.16 -0.09 0.01 0.01 0.10
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## ACF 0.34 0.34 0.33 0.33 0.31 0.29 0.25 0.21 0.18 0.14
## PACF -0.03 0.11 -0.08 -0.01 -0.09 -0.01 0.00 -0.03 -0.03 -0.02Vyzerá to tak, že vhodný model pre naše dáta by mohol byť AR(2). Dôvod je ten, že PACF(k), pre k > 2 je takmer nulová, zatiaľ čo ACF je nenulová a postupne klesá.
model.y <- sarima(y,2,0,0)## initial value 0.529468
## iter 2 value 0.509932
## iter 3 value -0.465627
## iter 4 value -0.609804
## iter 5 value -0.824525
## iter 6 value -0.824766
## iter 7 value -0.848374
## iter 8 value -0.848386
## iter 9 value -0.848405
## iter 10 value -0.848418
## iter 11 value -0.848425
## iter 12 value -0.848433
## iter 13 value -0.848448
## iter 14 value -0.848457
## iter 15 value -0.848460
## iter 16 value -0.848460
## iter 17 value -0.848460
## iter 18 value -0.848460
## iter 18 value -0.848460
## iter 18 value -0.848460
## final value -0.848460
## converged
## initial value -0.839506
## iter 2 value -0.839553
## iter 3 value -0.839778
## iter 4 value -0.839822
## iter 5 value -0.839886
## iter 6 value -0.839907
## iter 7 value -0.839912
## iter 8 value -0.839912
## iter 8 value -0.839912
## final value -0.839912
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 1.4389 0.0674 21.3471 0
## ar2 -0.5064 0.0677 -7.4840 0
## xmean 5.6034 0.4816 11.6340 0
##
## sigma^2 estimated as 0.1829345 on 158 degrees of freedom
##
## AIC = 1.207742 AICc = 1.208692 BIC = 1.284299
## 
Pre lag=4 z grafu nie je jasné, či je p-hodnota väčšia
ako 5% alebo nie. Overme si to priamo výpočtom
Box.test(...)$p.value## [1] 0.05157508Dostávame, že uvedená p-hodnota Ljung-Boxovho je väčšia ako 5%. Z grafu vidíme, že všetky ostatné p-hodnoty sú tiež väčšie ako 5%. Invertovateľnosť overovať netreba, keďže nemáme MA členy. Overíme ešte stacionaritu. Použitím regulárnym výrazov vieme k odhadnutým AR koeficientom pristupovať nasledovne:
koef <- model.y$fit$coef
koef_ar <- koef[grep("^ar",names(koef))]Absolútne hodnoty koreňom polynómu potom sú
abs(polyroot(c(1,-koef_ar)))## [1] 1.211715 1.629720Oba korene ležia mimo jednotkového kruhu, preto proces
sarima(y,2,0,0) je stacionárny.
5 Blízke koreňe AR a MA časti
Vygenerujeme dáta z ARMA(2,1) procesu:
set.seed(1234)
y <- arima.sim(model=list(ar=c(0.4,-0.1), ma=c(0.8)), n=250)
plot(y)
Odhadneme ARMA(3,2) model:
model32 <- sarima(y,3,0,2, details = FALSE)## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 1.2508 0.4442 2.8161 0.0053
## ar2 -0.5125 0.2458 -2.0851 0.0381
## ar3 0.1126 0.0821 1.3714 0.1715
## ma1 0.0941 0.4424 0.2126 0.8318
## ma2 -0.6013 0.3753 -1.6020 0.1105
## xmean 0.0233 0.2077 0.1122 0.9107
##
## sigma^2 estimated as 0.9982311 on 244 degrees of freedom
##
## AIC = 2.901126 AICc = 2.902509 BIC = 2.999727
## Graficky znázornite korene polynómov pre AR a pre MA časť procesu (ako body v komplexnej rovine, na x-ovej osi bude reálna časť čísla, na y-ovej imaginárna), pričom vhodným spôsobom odlíšte, ktorý koreň zodpovedá ktorému polynómu. Na grafe by ste mali vidieť, že jeden z koreňov AR časti je veľmi blízko jednému z koreňov MA časti.
Čomu tento výsledok nasvedčuje? Ako by malo ďalej pokračovať modelovanie týchto dát?
6 AR(1) proces pozorovaný s chybou
Vygenerjeme AR(1) proces \(y\):
set.seed(12345)
y <- arima.sim(n=200, list(ar = c(0.9)), sd=0.1)Pozrite si data frame wold, aby ste zistili, akú má
štruktúru.
Pozorovať ho ale budeme s chybou eps, vidíme teda proces
\(x\):
eps <- rnorm(200)/15
x <- y+epsZobrazte 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.
7 Dáta na samostatnú prácu
Zadanie a otázky sú vždy rovnaké ako v modelovom príklade.
7.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 <- BJsales
plot(y)
7.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)
7.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)
7.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)
7.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)