ARMA procesy
Beáta Stehlíková, 2-EMM-105 Časové rady, FMFI UK Bratislava
1 AR procesy - opakovanie
- Zistite, či sú stacionárne nasledovné AR procesy:
- \(x_t = 2 + 0.3 x_{t-1} + u_t\)
- \(x_t = 1 + 0.3 x_{t-1} - 0.8 x_{t-2} + u_t\)
Vypočítajte prvých 10 hodnôt ACF a PACF autoregresného procesu \(x_t = 1 + 0.3 x_{t-1} + 0.4 x_{t-2} + u_t\)
Zo stránky predmetu si stiahnite dáta
ovce.txt. Udávajú počet oviec v tisícoch v Anglicku a Walese, sú to ročné dáta z obdobia 1867-1939. Majú klesajúci trend, preto ich na dosiahnutie stacionarity zdiferencujeme. Nájdite pre tieto diferencie vhodný AR model - stacionárny a s dobrými rezíduami (na základe ACF a Ljung-Boxovho testu).
2 MA (moving average) procesy
2.1 Opakovanie z prednášky
Ako vyzerá vo všeobecnosti MA(q) proces?
Čo je charakteristické pre ACF moving average procesu?
Čo je invertovateľnosť MA procesu a ako vyzerá podmienka invertovateľnosti?
2.2 Overovanie invertovateľnosti MA procesu
Keďže ide o výpočet koreňov polynómu, rovnako ako v prípade overovania stacionarity AR procesu, môžeme použiť funkciu armaRoots z knižnice fArma. Treba si dať pozor na znamienka koeficientov, ktoré do nej vkladáme.
Príklad: Zistite, či je invertovateľný proces
\(x_t = 10 + u_t + 0.2 u_{t-1} - 0.3 u_{t-2}\)
\(x_t = 10 + u_t + 0.2 u_{t-1} + 0.3 u_{t-2} - 0.8 u_{t-3}\)
2.3 ACF a PACF moving average procesu
Analogicky ako v prípade AR procesov, napríklad pre proces \[x_t = 10 + u_t + 0.2 u_{t-1} - 0.3 u_{t-2}:\]
ARMAacf(ma=c(0.2,-0.3), lag.max=10)## 0 1 2 3 4 5
## 1.0000000 0.1238938 -0.2654867 0.0000000 0.0000000 0.0000000
## 6 7 8 9 10
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000ARMAacf(ma=c(0.2,-0.3), lag.max=10, pacf=TRUE)## [1] 0.123893805 -0.285214348 0.085799146 -0.100209762 0.045091942
## [6] -0.038884292 0.021248785 -0.015899253 0.009550012 -0.0066785103 Zmiešané ARMA procesy
3.1 Stacionarita a invertovateľnosť
Napíšte podmienky stacionarity a invertovateľnosti pre ARMA proces.
Pomocou funkcie
armaRootszistite, či je nasledujúci ARMA proces stacionárny a či je invertovateľný: \[x_t = 10 + 0.2 x_{t-1} - 0.3 x_{t-2} + u_t - 0.8 u_{t-1} + 0.1 u_{t-2}\]
3.2 ACF a PACF
- Nasledovný výpočet je výpočtom ACF a PACF procesu \[x_t = 5 + 0.3 x_{t-1} - 0.2 x_{t-2} + u_t + 0.8 u_{t-1} + 0.2 u_{t-2}\]
ARMAacf(ar=c(0.3,-0.2), ma=c(0.8,0.2), lag.max=10)## 0 1 2 3 4
## 1.000000e+00 6.125592e-01 6.907583e-02 -1.017891e-01 -4.435190e-02
## 5 6 7 8 9
## 7.052251e-03 1.098605e-02 1.885366e-03 -1.631601e-03 -8.665535e-04
## 10
## 6.635415e-05ARMAacf(ar=c(0.3,-0.2), ma=c(0.8,0.2), lag.max=10, pacf=TRUE)## [1] 0.612559242 -0.490024200 0.285060863 -0.137613544 0.054214131
## [6] -0.015934991 0.001908003 0.001660572 -0.001710051 0.001035927Úloha: Znázornite graficky túto ACF a PACF a všimnite si, že pri zmiešaných ARMA procesoch sa ani ACF ani PACF nevynuluje.
3.3 Odhadovanie ARMA modelu v R-ku
library(astsa)
sarima(data, p,k,q) # odhadujeme ARMA(p,q) model pre k-te diferencie premennej dataKeďže ACF a PACF nám neidentifikujú rád procesu, v prípade takejto ACF/PACF (keď sa ani jedna nezdá byť od istého indexu nulová), treba skúsiť zmiešaný ARMA proces najskôr nižšieho rádu, potom si ohraničiť parametre p a q a spomedzi modelov s vyhovujúcimi rezíduami sa rozhodnúť napríklad pomocou Akaikeho alebo Bayesovho informačného kritéria. Pripomeňme si, že lepšia je menšia hodnota AIC a BIC, cieľom je vybrať jednoduchší model, s malým počtom parametrov.
Cvičenie:
Načítajme dáta:
Uvažujme ARMA(p,q) modely, kde p,q sú menšie alebo rovné 3.
Ktoré z týchto modelov majú dobré rezíduá?
Ktorý z modelov s dobrými rezíduami má najnižšiu hodnotu ACF?
Ktorý má najnižšiu hodnotu BIC?
3.4 Postup pri ARMA modelovaní - zhrnutie
Zobrazíme výberovú ACF a PACF, naraz sa to dá spraviť funkciou
acf2(data)z balíkaastsa.“Podozrenie” na AR proces máme vtedy, ak sa PACF vynuluje, počet signifikantných hodnôt je rád procesu.
“Podozrenie” na MA proces máme vtedy, ak sa ACF vynuluje, počet signifikantných hodnôt je rád procesu.
Môžeme odhadovať aj zmiešané ARMA procesy s AR aj MA členmi.
ARMA(p,q) model pre k-te diferencie sa odhadne funkciou
sarima(data,p,k,q), analogicky predikciesarima.forTreba overiť nekorelovanosť rezíduí (ACF a Ljung-Boxovu štatistiku)
Ak existujú blízke korene z AR a MA časti, treba skúsiť ARMA(p-1,q-1) model.
Príklad zo skúšky 2017: Nájdite model pre dáta (bez diferencovania):
library(datasets)
y <- LakeHuronMusí mať dobré rezíduá (vysvetlite, ako ste ich testovali), byť stacionárny a invertovateľný (ukážte, ako ste overovali stacionaritu a inverovateľnosť).
4 Ďalšie príklady
4.1 Overovanie stacionarity a invertovateľnosti (kostra na skúške)
Príklad 3 z kapitoly 6.8 knihy Introductory Time Series with R od Paula S. P. Cowperweita a Andrewa V. Metcalfa (dostupná z fakultnej siete, odkaz na stránke cvičenia).
4.2 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, 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 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.1307446Graficky 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?
4.3 Porovnanie Woldových reprezentácií
Ukážte, že aj AR(4) aj MA(3) sú dobré modely pre diferencie premennej gas:
library(astsa)
data(gas)Porovnáme teraz ich Woldove reprezentácie a uvidíme, že tieto dva modely nie sú v skutočnosti až také odlišné.
Postup:
- Definujme
wold <- expand.grid(k=as.factor(1:10), model=c("ar4","ma3"))Pozrite si data frame wold, aby ste videli, akú má štruktúru.
- Potom pridajte hodnoty premennej
psis koeficientami \(\psi_k\) Woldovej reprezentácie pre príslušný model.
Keďže Woldova reprezentácia je vlastne MA(\(\infty\)) model, na prevod ARMA procesu do tohto tvaru slúži funkcia ARMAtoMA - pozrite si v helpe jej použitie.
wold$psi <- ...- Porovnanie pomocou ggplot2:
library(ggplot2)
ggplot(wold, aes(x=k, y=psi, fill=model)) +
geom_bar(stat="identity", position="dodge")4.4 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+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 generovaliAle 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.
4.5 Kniha “tsaEZ”, Problem 3.3
Odkaz na knihu je na stránke predmetu.
Poznámka: \(\psi\)-weights v zadaní sú parametre Woldovej MA(\(\infty\)) reprezentácie stacionárneho (causal) procesu, \(\pi\)-weights sú parametre AR(\(\infty\)) reprezentácie invertovateľného procesu.
4.6 Kniha “tsaEZ”, Problem 3.13, 3.14
Ďalšie príklady na hľadanie ARMA modelov pre zadané dáta.
Hľadanie ARMA modelu pre zadané dáta je súčasťou kostry na skúške.
4.7 Nájdenie procesov s danými vlastnosťami
Nájdite príkladov procesov s danými vlastnosťami. Pre každý z nich dokážte, že má požadovanú vlastnosť. Ak nie je povedané inak, požadujeme aj stacionaritu daného procesu.
ARMA(2,2) proces, ktorý je stacionárny, ale nie je invertovateľný.
ARMA(2,2) proces, ktorého stredná hodnota je 25.
Proces, ktorého Woldova reprezentácia obsahuje len konečne veľa nenulových členov.
Proces, ktorého Woldova reprezentácia obsahuje nekonečne veľa nenulových členov.
Nestacionárny AR proces.
MA(2) proces so strednou hodnotou 17.
AR(2) proces, ktorého autokorelačná funkcia nadobúda aj kladné, aj záporné hodnoty.
Proces, ktorého ACF má len konečne veľa nenulových členov.
Proces, ktorého ACF aj PACF má nekonečne veľa nenulových členov.