\[x_t = \delta + \alpha_1 x_{t-1} + \dots \alpha_p x_{t-p} + u_t - \beta_1 u_{t-1} - \dots - \beta_q u_{t-q}\]
Stacionarita závisí iba od AR časti (stacionarita, resp. nestacionarita je rovnaká ako pre proces \(x_t = \delta + \alpha_1 x_{t-1} + \dots \alpha_p x_{t-p} + u_t\))
Uvažujme procesy \[x_t = 10 + 0.2 x_{t-1} - 0.1 x_{t-2} + u_t - 0.3 u_{t-1} + 0.1 u_{t-2}\] \[x_t = 5 - 0.8 x_{t-1} + 0.1 x_{t-2} + 0.7 x_{t-3} + u_t - 0.2 u_{t-1} - 0.5 u_{t-2}\] \[x_t = 1.5 + 1.8 x_{t-1} + u_t - 0.3 u_{t-1} + 0.1 u_{t-2}\]
Zistite, či sú stacionárne. Pre stacionárne procesy vypočítajte ich stednú hodnotu a pomocou funkcie ARMAacf
vypočítajte prvých 10 hodnôt ich autokorelačnej funkcie.
Uvažujme nasledovné dáta (podľa popisu v helpe: Annual measurements of the level, in feet, of Lake Huron 1875–1972.)
library(datasets)
x <- window(LakeHuron, start=1900, end=1967)
plot(x)
LakeHuron
Pripomeňme si Woldovu reprezentáciu procesu z prvého týždňa - každý stacionárny proces sa dá zapísať v tvare nekonečnej sumy \[y_t = \mu + u_t + \psi_1 u_{t-1} + \psi_2 u_{t-2} + \psi_3 u_{t-3} + \dots\]
Na výpočet koeficientov v R-ku sa dá použiť funkcia ARMAtoMA
(keďže ide o nekonečný MA proces).
Príklad: Woldova reprezentácia procesu \[y_t = 10 + 0.5 y_{t-1} + 0.1 y_{t-2} + u_t + 0.1 u_{t-1} - 0.15 u_{t-2}\] má prvých 10 koeficientov Woldovej reprezentácie nasledovných:
ARMAtoMA(ar = c(0.5, 0.1), ma = c(0.1, -0.15), lag.max = 10)
## [1] 0.600000000 0.250000000 0.185000000 0.117500000 0.077250000
## [6] 0.050375000 0.032912500 0.021493750 0.014038125 0.009168438
Vráťme sa k dátam o Hurónskom jazere:
library(datasets)
x <- window(LakeHuron, start=1900, end=1967)
wold1 <- ...
wold2 <- ...
barplot(rbind(wold1, wold2), beside = TRUE,
legend.text = c("ARMA(1,1)","AR(2)"),
col=c("red","blue"))
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:
library(astsa)
sarima(y,3,0,2)
## initial value 0.575413
## iter 2 value 0.240783
## iter 3 value 0.156571
## iter 4 value 0.078853
## iter 5 value 0.037933
## iter 6 value 0.019129
## iter 7 value 0.009254
## iter 8 value 0.008716
## iter 9 value 0.008098
## iter 10 value 0.007761
## iter 11 value 0.007584
## iter 12 value 0.007302
## iter 13 value 0.006680
## iter 14 value 0.005467
## iter 15 value 0.005314
## iter 16 value 0.004718
## iter 17 value 0.004672
## iter 18 value 0.004523
## iter 19 value 0.004447
## iter 20 value 0.004386
## iter 21 value 0.004135
## iter 22 value 0.003971
## iter 23 value 0.003749
## iter 24 value 0.003514
## iter 25 value 0.003035
## iter 26 value 0.002584
## iter 27 value 0.002240
## iter 28 value 0.002168
## iter 29 value 0.002144
## iter 30 value 0.002061
## iter 31 value 0.002055
## iter 32 value 0.002037
## iter 33 value 0.002036
## iter 34 value 0.002036
## iter 35 value 0.002036
## iter 36 value 0.002036
## iter 37 value 0.002036
## iter 38 value 0.002036
## iter 39 value 0.002036
## iter 40 value 0.002036
## iter 41 value 0.002036
## iter 42 value 0.002036
## iter 43 value 0.002036
## iter 44 value 0.002036
## iter 45 value 0.002036
## iter 46 value 0.002036
## iter 47 value 0.002036
## iter 48 value 0.002036
## iter 49 value 0.002036
## iter 50 value 0.002036
## iter 51 value 0.002036
## iter 52 value 0.002036
## iter 53 value 0.002036
## iter 53 value 0.002036
## final value 0.002036
## converged
## initial value 0.003822
## iter 2 value 0.003766
## iter 3 value 0.003724
## iter 4 value 0.003674
## iter 5 value 0.003654
## iter 6 value 0.003643
## iter 7 value 0.003641
## iter 8 value 0.003641
## iter 9 value 0.003641
## iter 10 value 0.003640
## iter 11 value 0.003639
## iter 12 value 0.003637
## iter 13 value 0.003634
## iter 14 value 0.003629
## iter 15 value 0.003627
## iter 16 value 0.003626
## iter 17 value 0.003626
## iter 18 value 0.003626
## iter 19 value 0.003626
## iter 20 value 0.003626
## iter 21 value 0.003626
## iter 22 value 0.003626
## iter 23 value 0.003626
## iter 24 value 0.003626
## iter 25 value 0.003625
## iter 26 value 0.003625
## iter 27 value 0.003625
## iter 28 value 0.003625
## iter 29 value 0.003625
## iter 30 value 0.003625
## iter 31 value 0.003624
## iter 32 value 0.003624
## iter 33 value 0.003624
## iter 34 value 0.003624
## iter 35 value 0.003624
## iter 36 value 0.003624
## iter 37 value 0.003624
## iter 38 value 0.003624
## iter 39 value 0.003624
## iter 40 value 0.003624
## iter 41 value 0.003624
## iter 42 value 0.003624
## iter 43 value 0.003624
## iter 44 value 0.003624
## iter 44 value 0.003624
## iter 44 value 0.003624
## final value 0.003624
## converged
## $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 ar2 ar3 ma1 ma2 xmean
## 1.2508 -0.5125 0.1126 0.0941 -0.6013 0.0233
## s.e. 0.4442 0.2458 0.0821 0.4424 0.3753 0.2077
##
## sigma^2 estimated as 0.9982: log likelihood = -355.64, aic = 725.28
##
## $degrees_of_freedom
## [1] 244
##
## $ttable
## 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
##
## $AIC
## [1] 1.046229
##
## $AICc
## [1] 1.056081
##
## $BIC
## [1] 0.1307446
Cvičenie: 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íšite, 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?
Vráťme sa k dátam o Hurónskom jazere:
library(datasets)
x <- window(LakeHuron, start=1900, end=1967)
pri ktorých sme videli, že aj ARMA(1,1) aj AR(2) sú vhodné modely.
Porovnajte autokorelačnú funkciu týchto dvoch modelov (podobne, ako sme porovnávali ich Woldove reprezentácie) a výberovú autokorelačnú funkciu odhadnutú z modelovaných dát. (Teda na stĺpcovom grafe budú pri každom lagu tri hodnoty - odhad ACF, ACF ARMA(1,1) modelu, ACF AR(2) modelu.)
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 + \varepsilon\), kde \(\varepsilon\) je proces definovaný tak, ako v simuláciách hore, tak \(x\) sa dá zapísať ako ARMA(1,1) proces.