1 Zmiešané ARMA procesy

\[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\))

1.1 Príklad: stacionarita a momenty

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.

1.2 Odhadovanie v R-ku

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)

  • Odhadnite ARMA(1,1) model pre tieto dáta.
  • Overte jeho stacionaritu a zhodnoťte rezíduá.
  • Spravte predikciu pre nasledujúcich 5 rokov a graficky porovnajte s pozorovanými hodnotami v premennej LakeHuron

2 Woldova reprezentácia procesu

2.1 Definícia

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\]

2.2 Výpočet v R-ku

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

2.3 Porovnanie Woldových reprezentácií odhadnutých modelov

Vráťme sa k dátam o Hurónskom jazere:

library(datasets)
x <- window(LakeHuron, start=1900, end=1967)
  • Ukážte, že AR(1) model nie je pre tieto dáta vyhovujúci.
  • Ukážte, že aj ARMA(1,1) aj AR(2) modely majú dobré rezíduá.
  • Porovnajte ich Woldove reprezentácie
wold1 <- ...
wold2 <- ...

barplot(rbind(wold1, wold2), beside = TRUE,
        legend.text = c("ARMA(1,1)","AR(2)"),
        col=c("red","blue"))

3 Cvičenia

3.1 Blízke korene 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:

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?

3.2 Porovnanie autokorelačných funkcií

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

3.3 AR(1) proces pozorovaný s chybou

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.