Cvičenie 7: Testovanie jednotkového koreňa, ARIMA modelovanie
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 7
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 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"))3 Simulované dáta
Predtým ako prejdeme k práci s reálnymi dátami si postup ukážeme na simulovaných dátach, s ktorými ste sa stretli na prednáške. Dáta do R-ka načítajte príkazom
load(url("https://anna-hlubinova.github.io/ARIMAmodely.Rdata"))V premenných y1, y2, y3, y4 sú uložené 4 simulované časové rady dĺžky 100. Každý časový rad si vykreslite a otestujte na prítomnosť jednotkového koreňa. Napokon nájdite vhodný ARIMA model pre dáta.
- Časový rad y1:
plot(y1)
- Časový rad y2:
plot(y2)
- Časový rad y3:
plot(y3)
- Časový rad y4:
plot(y4)
4 Modelový príklad
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.
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:
- 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(model.y$fit$residuals, fitdf = 2, lag = 4, type = "Ljung-Box")$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 Dáta na samostatnú prácu
Zadanie a otázky sú vždy rovnaké ako v modelovom príklade.
6.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)
6.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)
6.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)
6.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)
6.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)