Autoregresné procesy
Časové rady, Beáta Stehlíková, FMFI UK Bratislava
1 Opakovanie z prednášky: odhadnutie AR(1) modelu, interpretácia výstupu pre rezíduá
x <- window(Nile, start = 1905)
plot(x)
- Pre dáta o prietokoch Nílu sme Ljung-Boxovým testom zistili, že majú signifikantné autokorelácie. Nejde preto iba o biely šum posunutý o konštantu, ale treba modelovať závislosť pozorovaných hodnôt.
- Odhadnite pre tieto dáta AR(1) model a na základe Ljung-Boxovho testu pre rezíduá (ktorý je súčasťou výstupu) rozhodnite, či ide o dobrý model pre tieto dáta.
- Prečo sa Ljung-Boxov test začína až od testovania prvých dvoch autokorelácií (pri testovaní pôvodných dát sme testovali aj prvú koreláciu)?
- Zopakujte pre modely AR(2), AR(3) a AR(4) - majú dobré rezíduá?
- Ktorý z modelov s dobrými rezíduami má najlepšiu hodnotu BIC?
- Pomocou modelu z predchádzajúceho bodu spravte predikcie pre niekoľko nasledujúcich rokov.
2 Zapísanie odhadnutého modelu
Uvažujme AR(2) model z predchádzajúceho príkladu. Zapíšeme rovnicu, ktorá je výsledkom jeho odhadovania, pričom konštantný člen vypočítame pomocou pristupovania k presným hodnotám parametrov (nie použitím zaokrúhlených hodnôt z výstupu).
3 Hľadanie modelu, porovnanie predikcií s dátami
Úloha typu "Nájdite dobrý model pre zadané dáta" je súčasťou kostry na skúške.
Uvažujme nasledovné dáta (Passenger Miles on Commercial US Airlines, 1937–1960), pričom kvôli rastúcej disperzii budeme modelovať a predikovať logaritmy:
par(mfrow = c(1, 2)) # kreslenie obrazkov: 1 riadok, 2 stlpce
plot(airmiles)
plot(log(airmiles))
par(mfrow = c(1, 1)) # vratime naspat
- Vynechajme z dát posledné štyri roky, použijeme ich na porovnanie predikcií s realizovanými hodnotami.
- Čo pre naše modelovanie vyplýva z toho, že dáta majú rastúci trend?
- Nájdite model s dobrými rezíduami a zapíšte ho.
- Spravte pomocou modelu z predchádzajúceho bodu predikcie
log(airmiles)na štyri roky. - Do grafu s predikciami pridajte skutočné hodnoty.
Ukážka možného výstupu:
4 Zautomatizovanie testovania viacerých AR modelov
- Prípravná úloha: Odhadujeme pre dáta o Níle AR(2) model. Čomu sa rovná p-hodnota z Ljung-Boxovho testu pre rezíduá, ak testujeme nulovosť prvých piatich autokorelácií?
# POSTUP
model <- sarima(...) # odhadneme model
model$... # chceme pristupovat k reziduam
Box.test(..., fitdf = ) # Ljung-Boxov test pre rezidua- Napíšte funkciu, ktorá pre zadané dáta, rád autoregresného procesu a maximálny uvažovaný počet lagov vráti TRUE/FALSE podľa toho, či sú všetky p hodnoty Ljung-Boxovho testu pre rezíduá väčšie ako 0.05.
ARmodel <- function(data, p, lag.max = 20){ # nastavenie defaultnej hodnoty vstupneho parametra
# DOPLNTE
}
# testujeme prvych k autokorelacii (k<=15) rezidui z AR(5) modelu pre data x
ARmodel(data = x, p = 5, lag.max = 15)
ARmodel(x, 5, 15) - Pomocou funkcie
sapplyzistite, ktoré AR(p) modely s rádom nanajvýš 5 majú dobré rezíduá. - Upravte funkciu
ARmodeltak, aby vrátila aj hodnotu Bayesovho informačného kritéria. Pomocou tejto funkcie potom spomedzi modelov s dobrými rezíduami vyberte ten, ktorý má najnižšiu hodnotu informačného kritéria.
5 Overovanie stacionarity AR procesov I.
Overovanie stacionarity procesu je súčasťou kostry na skúške.
- Zistite, či je stacionárny proces \(x_t = 10 + 0.5 x_{t-1} + 0.1 x_{t-2} + u_t\).
- Zistite, či je stacionárny proces \(x_t = 10 + 0.8 x_{t-1} + 0.1 x_{t-2} - 0.2 + u_t\).
- Zistite, či je stacionárny proces \((1 - 0.5 L - 0.4 L^2 + 0.2 L^3) x_t = 10 + u_t\).
6 Grafické znázornenie koreňov
Výstup z funkcie polyroot (komplexné čísla) sa dá priamo kresliť funkciou plot. Na x-ovej osi bude reálna časť, na y-ovej osi bude imaginárna časť, tomuto budú zodpovedať aj popisy osí (tie sa dajú zmeniť štandardným spôsobom nastavením parametrov xlab, ylab). Keďže budeme kvôli stacionarite kresliť jednotkovú kružnicu, oplatí sa nastaviť pri kreslení koreňov parameter asp = 1, aby potom jednotková kružnica vyzerala naozaj ako kružnica, nie ako elipsa.
z <- polyroot(c(1, 0.8, -0.4, 0.2)) # korene 1 + 0.8*x - 0.4*x^2 + 0.1*x^3 = 0
z## [1] -0.8008768-0.000000i 1.4004384+2.069282i 1.4004384-2.069282iplot(z, asp = 1)
Na kreslenie jednotkovej kružnice sa dá výhodne použiť funkcia curve s parametrom add = TRUE (pridá sa do predchádzajúceho grafu).
- Pridajte do obrázku jednotkovú kružnicu.
- Pridajte do obrázku osi.
- Prípadne ďalšie formátovanie.
7 Overovanie stacionarity AR procesov II.
Ak riešime úlohy typu "pre aké hodnoty parametra je daný proces stacionárny" analyticky, je dobré mať numerickú kontrolu.
Príklad 1. Pre aké hodnoty parametra \(k\) je stacionárny proces \(x_t= x_{t−1} + k x_{t−2} + u_t\)?
Numerická simulácia. Vytvoríme vektor hodnôt kVektor a pre každú jeho zložku zistíme, či je získaný proces stacionárny.
kVektor <- seq(from = -3, to = 3, by = 0.1)
stacionarita <- sapply(kVektor, function(k) .... )Dopíšte funkciu, ktorá pre zadanú hodnotu k vráti hodnotu TRUE alebo FALSE podľa toho, či je proces stacionárny. Tieto hodnoty sa potom dajú napríklad kresliť do grafu ako hodnoty 1 a 0:
plot(kVektor, stacionarita)
Alebo si vypíšeme:
kVektor[stacionarita] # znamena: kVektor[stacionarita == TRUE]## [1] -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1kVektor[!stacionarita] # znamena: kVektor[stacionarita == FALSE]## [1] -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 0.0 0.1 0.2 0.3
## [16] 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## [31] 1.9 2.0Príklad 2. Pre aké hodnoty parametra \(k\) je stacionárny proces \(x_t = x_{t−1} + k x_{t−2} − k x_{t−3} + u_t\)?
Spravte numerickú simuláciu pre tento príklad.
8 Stacionarita AR(2) procesu
Vytvorte funkciu, ktorá pre zadaný vektor parametrov \((\alpha_1, \alpha_2)\) vráti TRUE, resp. FALSE podľa toho, či je AR(2) proces s týmito koeficientami stacionárny alebo nie.
stacAR2 <- function(alpha12) { # vstup: dvojzlozkovy vektor
....
}Otestujte:
stacAR2(c(1.2, -0.5))## [1] TRUEstacAR2(c(1.2, 0.5))## [1] FALSEstacAR2(c(0.5, 0)) # ar(1), ale stacionarny## [1] TRUEstacAR2(c(0, 0)) # konstanta + biely sum -> stacionarny## [1] TRUESpravíme sieť bodov \((\alpha_1, \alpha_2)\) a zakreslíme, ktoré sú stacionárne a ktoré nie sú. Na vytvorenie dvojíc bodov použijeme náhodne generované čísla z rovnomerného rozdelenia na štvorci \((-2,2)\times(-2,2)\).
set.seed(123)
alfa1 <- runif(1000, min = -2, max = 2)
alfa2 <- runif(1000, min = -2, max = 2)
df <- data.frame(alfa1, alfa2)
head(df)## alfa1 alfa2
## 1 -0.8496899 -0.90550906
## 2 1.1532205 0.37546773
## 3 -0.3640923 -1.35926075
## 4 1.5320696 1.41372096
## 5 1.7618691 1.39095663
## 6 -1.8177740 -0.08845275Teraz pridáme do data framu tretí stĺpec, v ktorom bude výsledok overovania stacionarity. Namiesto for-cyklu použijeme funkciu apply. Tá umožňuje aplikovať zadanú funkciu na riadky alebo stĺpce objektu (matica, data frame).
Ukážka použitia apply
set.seed(123)
dfUkazka <- data.frame(a = sample.int(5),
b = sample.int(5),
c = sample.int(5))
dfUkazka # nejake cisla## a b c
## 1 3 3 2
## 2 2 1 3
## 3 5 2 1
## 4 4 5 4
## 5 1 4 5f <- function(x) (x[1] + x[2]) * x[3] # nejaka funkcia
# POUZITIE APPLY
apply(dfUkazka, # vstupne data pre funkciu
MARGIN = 1, # po riadkoch
FUN = f) # pouzita funkcia## [1] 12 9 7 36 25# priradenie ako novy stlpec
dfUkazka$vysledok <- apply(dfUkazka, 1, f)
dfUkazka## a b c vysledok
## 1 3 3 2 12
## 2 2 1 3 9
## 3 5 2 1 7
## 4 4 5 4 36
## 5 1 4 5 25# DALSIE POUZITIE APPLY: teraz nie je vstupom cely riadok
g <- function(x) x[2]/x[1] + x[1] # zase nejaka funkcia
dfUkazka$vysledok2 <- apply(dfUkazka[ , 3:4], # vstupom do g bude 3. a 4. stlpec
1, g)
dfUkazka## a b c vysledok vysledok2
## 1 3 3 2 12 8
## 2 2 1 3 9 6
## 3 5 2 1 7 8
## 4 4 5 4 36 13
## 5 1 4 5 25 10Vrátime sa k overovaniu stacionarity: Pomocou apply a vytvorenej funkcie stacAR2 pridajte do data framu df stĺpec s informáciou o stacionarite procesu s danymi parametrami.
# na kontrolu
head(df)## alfa1 alfa2 stacionarita
## 1 -0.8496899 -0.90550906 TRUE
## 2 1.1532205 0.37546773 FALSE
## 3 -0.3640923 -1.35926075 FALSE
## 4 1.5320696 1.41372096 FALSE
## 5 1.7618691 1.39095663 FALSE
## 6 -1.8177740 -0.08845275 FALSEVykreslite body farebne odlíšené podľa toho, či je proces stacionárny alebo nie.
9 Ďalšie príklady na precvičenie
9.1 Hľadanie autoregresného modelu pre zadané dáta
Všetky nasledujúce zadania sú zo starých skúšok. Úlohou je
- nájsť vhodný autoregresný model pre zadané dáta
- vysvetliť, prečo sú rezíduá vyhovujúce - skomentovať ACF rezíduí aj Ljung-Boxov test, pričom treba presne povedať, aká hypotéza sa testuje a či sa v tomto prípade zamieta alebo nie (a prečo sme s tým výsledkom spokojní)
- overiť stacionaritu získaného modelu (napísať polynóm, ktorého korene overujete, aké absolútne hodnoty koreňov vyšli a prečo sme s tým spokojní)
Na skúške je súčasťou kostry zistenie, či môžeme pracovať priamo so zadanými dátami alebo s diferenciami (trend nie je jediným dôvodom diferencovannia - budeme sa tomu venovať neskôr). V týchto zadaniach je povedané, s akými dátami máte pri hľadaní AR modelu pracovať. V prípade modelovania diferencií sa treba na základe prítomnosti, resp. neprítomnosti trendu treba rozhodnúť, či do modelu zahrnúť konštantu alebo nie.
Takisto je na skúške na výber širšia trieda modelov, ako sú AR modely. V niektorých prípadoch - vrátane týchto - však AR model stačí.
9.1.1 Ceny benzínu
Dáta z knižnice astsa, z popisu v helpe: New York Harbor conventional regular gasoline weekly spot price FOB (in cents per gallon) from 2000 to mid-2010. Zoberieme dáta od roku 2006.
library(astsa)
x <- window(gas, start = c(2006, 1))
plot(x)
Nájdite model pre premennú x tak, že budete jej diferencie modelovať AR procesom.
9.1.2 Indikátor od Boxa a Jenkinsa
Dáta z knižnice astsa, z popisu v helpe: Leading indicator, 150 months; taken from Box and Jenkins (1970).
library(astsa)
x <- lead
plot(x)
Nájdite model pre premennú x tak, že budete jej diferencie modelovať AR procesom.
9.1.3 Hladina Hurónskeho jazera
Dáta z knižnice datasets, z popisu v helpe: Annual measurements of the level, in feet, of Lake Huron 1875–1972. Použijeme dáta od roku 1900.
library(datasets)
x <- window(LakeHuron, start = 1900)
plot(x)
9.2 Stacionarita AR(2) procesu
Odvoďte podmienku stacionarity AR(2) procesu v tvare lineárnych nerovníc pre autoregresné koeficienty a zakreslite do roviny (na jednotlivých osiach sú hodnoty autoregresných koeficientov) tie body, ktoré zodpovedajú stacionárnym procesom. Porovnajte so simuláciami, ktoré sa robili v cvičení 8.
9.3 Teoretické cvičenie
Nech \(x_t\) je stacionárny proces. Odvoďte strednú hodnotu, disperziu a autokorelačnú funkciu procesu \(y_t = x_t - x_{t-1}\).
9.4 Príprava dát a hľadanie AR modelu
Pomocou balíka priceR a funkcie retrieve_inflation_data získame dáta o inflácii z World Bank API.
library(priceR) # ak treba, nainstalujte
data <- retrieve_inflation_data("US")## Validating iso2Code for US
## Generating URL to request all 297 results
## Retrieving inflation data for US
## Generating URL to request all 63 resultsstr(data)## 'data.frame': 63 obs. of 8 variables:
## $ indicator :'data.frame': 63 obs. of 2 variables:
## ..$ id : chr "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" "FP.CPI.TOTL.ZG" ...
## ..$ value: chr "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" "Inflation, consumer prices (annual %)" ...
## $ country :'data.frame': 63 obs. of 2 variables:
## ..$ id : chr "US" "US" "US" "US" ...
## ..$ value: chr "United States" "United States" "United States" "United States" ...
## $ countryiso3code: chr "USA" "USA" "USA" "USA" ...
## $ date : chr "2022" "2021" "2020" "2019" ...
## $ value : num 8 4.7 1.23 1.81 2.44 ...
## $ unit : chr "" "" "" "" ...
## $ obs_status : chr "" "" "" "" ...
## $ decimal : int 1 1 1 1 1 1 1 1 1 1 ...Príprava dát: Vytvorte z týchto dát časový rad hodnôt inflácie.
Hľadanie modelu.
- Neskôr uvidíme, že takéto dáta treba diferencovať, aj keď v nich nie je trend.
- Preto teraz budeme pracovať s diferenciami. Má ísť o model s konštantou alebo bez konštanty?
- Pre diferencie zistite či sú pre ne na základe rezíduí dobrými modelmi AR(p) pre \(p \in \{0, 1, 2,3,4,5 \}\) bez konštanty
- Z vyhovujúcich modelov vyberte ten, ktorý má najnižšie BIC.
Ďalšie dáta. Zopakujte pre infláciu v štátoch (všade treba diferencovať):
- Francúzsko (FR)
- Japonsko (JP)
- Kanada (CA)