ARMA modeli
Ivana Bogicevic
April 15, 2019
Podsecanje
Pre nego sto teorijski i prakticno objasnimo i prikazemo ARMA modele prisetimo se nekih stvari sa prethodnih casova. Pogledajmo neke primere vremenskih serija(sve vremenske serije sa nalaze u paketu astsa):
Johnson & Johnson vremenska serija (jj)
Ova vremenska serija predstavlja zaradu po akciji kvartalno gledano. Vremenska serija ima 84 kvartala, odnosno podaci su sakupljeni kroz 21 godinu. Serija pocinje od prvog kvartala 1960. godine do poslednjeg kvartala 1980. godine.
plot(jj, main = "Johnson & Johnson", type = "c")
text(jj, labels = 1:4, col = 1:4)
Vremenska serija ima neke karakteristicne i vrlo ceste osobine za ovakvu vrstu podataka. Serija ima trend (rastuci), sezonsku komponentu i drugi i treci kvartal su obicno visoki dok je 4 kvartal nizak. Ova serija ima i osobinu heteroskedasticnosti sto znaci da kako raste vrednost akcija tako procentualno male razlike postaju velike razlike po apsolutnoj vrednosti.
Globtemp vremenska serija
Ova serija predstavlja odstupanja od godisnjih temperatura (iz perioda 1951-1980) u periodu 1880-2015 godine.
plot(globtemp, main = "Global temperature deviations", type = "o")
Primeticemo da serija ima relativno pozitivan trend, ali vidimo da trend nije uvek pozitivan. Za raziku od Johnson & Johnson serije ova serija nema sezonsku komponentu.
S&P 500 vremenska serija
plot(sp500w, main = "S&P 500", type = "c")
Za razliku od prethodne dve serije ova serija nema ni trend ni sezonsku komponentu. Zapravo, uocavamo da ova serija nema nikakav obrazac “ponasanja”, osim da je s vremena na vreme disperzija velika. Ova je primer posebne vrste procesa koji se naziva sum ( noise ).
Regresioni modeli vremenskih serija
Regresija: \(Y_{i}=\beta X_{i} + \epsilon_{i}\), pri cemu je \(\epsilon_{i}\) beli sum
AutoRegresija (AR): \(X_{t}=\phi X_{t-1}+\epsilon_{t}\) (\(\epsilon_{t}\) je beli sum)
Obicno su podaci iz vremenskih serija korelisani i pretpostvljati da greske nisu korelisane moglo bi da dovede do gresaka u predvidjanju. Jedan od nacina da se ovaj problem resi je da se koristi metoda pokretnih proseka za greske modela.
Moving average (MA): \(\epsilon_{t}=W_{t}+\theta W_{t-1}\), (\(W_{t}\) je beli sum).
Sklapanje i razumevanje svega pomenutog nas dovodi da ARMA modela. Drugim recima, model ARMA je autoregresija sa autokorelisanim greskama.
ARMA: \(X_{t}=\phi X_{t-1}+W_{t}+\theta W_{t-1}\)
Stacionarnost i nestacionarnost
U kontekstu vremenskih serija stacionarnost predstavlja konstantnost srednje vrednosti. Kazemo da je vremenska serija stabilna ako je srednja vrednost konstanta (odnosno da nema trenda) i da je korelaciona struktura podataka iz serije konstantna tokom vremena.
Vremenska serija globtemp nije stacionarna kao sto smo videli. Ali diferenciranje ove serije dovodi do stacionarnosti.
par(mfrow=c(2,1))
plot(globtemp, main = "globtemp")
plot(diff(globtemp), main = "diff(globtemp)")
U slucaju da postoji i trend i heteroskedasticnost, kao kod Johnson&Johnson serije, logaritmovanje pa diferenciranje dace nam stacionarnu vremensku seriju. Prvo logaritmovanje stabilizuje disperziju, a zatim ce diferenciranje ukloniti trend.
par(mfrow=c(3,1))
plot(jj, main="J&J")
plot(log(jj), main="log(J&J)")
plot(diff(log(jj)), main="diff(log(J&J))")
Stacionarne vremenske serije: ARMA
ARMA(p, q) procesi su stacionarni procesi koji se mogu zapisati u obliku
\(x_{t}-\mu = \phi_{1}(x_{t-1}-\mu)+ ... +\phi_{p}(x_{t-p}-\mu)+ \omega_{t}+ \theta_{1} \omega_{t-1}+ ... + \theta_{q} \omega_{t-q}\)
gde su \(\phi_{p}\neq 0, \theta_{q}\neq0\), a \(\omega_{t}\) je beli sum. Vrednost \(\mu\) predstavlja srednju vrednost procesa i ako se uvede oznaka \(\alpha=\mu(1-\phi_{1}-...-\phi_{p})\) prethodna jednakost se svodi na
\(x_{t}= \alpha+ \phi_{1}x_{t-1}+ ... +\phi_{p}x_{t-p}+ \omega_{t}+ \theta_{1} \omega_{t-1}+ ... + \theta_{q} \omega_{t-q}\)
AR i MA modeli su specijalan slucaj ARMA modela. S obzirom da smo na prethodnim casovima definisali AR i MA polinome, ARMA model se moze zapisati kao
\(\Phi(B)x_{t}=\Theta(B)\omega_{t}\)
ARMA model je kauzalan ako se moze zapisati kao bekonacni MA model, dok je ARMA model invertibilan ako se moze zapisati kao beskonacni AR model. Pre nego sto se zapocne ispitivanje invertibilnosti i kauzalnosti ARMA modela bitno je proveriti redudantnost parametara. Moze da se desi da AR i MA polinomi imaju iste faktore. Kada se dobije jednacina nekog ARMA modela i treba da se ispita koji je model u pitanju, najpre treba proveriti AR i MA polinome i ako oni imaju zajednican faktor-skratiti ga! Ukoliko to ne uradimo desava se da pogresno prepoznamo model tj. da neki jednostavni model prepoznamo kao neki znatno komplikovaniji.
Kazemo da je ARMA(p,q) model kauzalan ako se moze zapisati kao linearan proces
\(x_{t}=\Psi(B)\omega_{t}\), gde je \(\Psi(B)=\sum_{j=0}^{\infty}\psi_{j}B^{j}\) i \(\sum_{j=o}^{\infty}\mid\psi_{j}\mid<\infty\), \(\psi_{0}=1\). Ovo vazi kada su svi koreni AR polinoma \(\Psi(B)\) veci od 1 po normi.
Kazemo da je ARMA(p,q) model invertibilan ako se moze zapisati kao
\(\omega_{t}=\pi(B)x_{t}\), gde je \(\pi(B)=\sum_{j=0}^{\infty}\pi_{j}B^{j}\) i \(\sum_{j=o}^{\infty}\mid\pi_{j}\mid<\infty\), \(\pi_{0}=1\). Ovo vazi kada su svi koreni AR polinoma \(\Phi(B)\) veci od 1 po normi.
Odgovor na to zasto je validno koristiti ARMA model za stacionarne vremenske serije dao je Herman Wold koji je pokazao da se svaki stacionarni model moze prestaviti kao linearna kombinacija belih sumova. Takodje se moze pokazati da je ARMA model takodje linearna kombinacija belih sumova. R kao sto vec znamo da pruza vrlo lak nacin da simuliramo ARMA modele arima.sim().
arima.sim(model, n,..)
arima.sim(list(order=c(p ,d, q)))
Generisemo MA(1) model \(X_{t}=W_{t}-0.9W_{t-1}\)
x=arima.sim(list(order=c(0, 0, 1), ma=0.9), n=100)
plot(x)
Generisemo AR(2) model \(X_{t}=-0.9X_{t-2}+W_{t}\)
x=arima.sim(list(order=c(2, 0, 0), ar=c(0, -0.9)), n=100)
plot(x)
Generisemo ARMA(2,1) model
x=arima.sim(list(order=c(2, 0, 1), ar=c(0, -0.9), ma=0.9), n=100)
plot(x)
Kako da prepoznamo u koji ARMA model da uklopimo seriju?
Kao sto vidimo AR(2) i MA(1) na grafiku izgledaju vrlo slicno. Jasno je da se samo na grafika vremenske serije ne moze prepoznati koji je model u pitanju. Moraju se posmatrati ACF i PACF funkcije. Pogledajmo sta se desava u slucaju AR(2) modela 
## ACF PACF
## [1,] 0.01 0.01
## [2,] -0.89 -0.89
## [3,] -0.01 0.03
## [4,] 0.78 -0.02
## [5,] 0.01 -0.02
## [6,] -0.69 -0.03
## [7,] -0.01 0.00
## [8,] 0.63 0.04
## [9,] 0.01 0.04
## [10,] -0.57 -0.02
## [11,] -0.02 0.02
## [12,] 0.51 0.00
## [13,] 0.02 0.00
## [14,] -0.46 0.03
## [15,] 0.00 0.06
## [16,] 0.41 0.01
## [17,] 0.00 0.02
## [18,] -0.38 -0.04
## [19,] 0.01 0.03
## [20,] 0.34 -0.04
## [21,] -0.03 -0.07
## [22,] -0.30 -0.01
## [23,] 0.04 -0.04
## [24,] 0.27 -0.03
## [25,] -0.04 0.03
## [26,] -0.23 0.05
## [27,] 0.05 0.05
## [28,] 0.19 0.01
## [29,] -0.05 0.03
## [30,] -0.16 0.02
## [31,] 0.05 0.03
## [32,] 0.14 0.00
## [33,] -0.06 -0.01MA(1) model
## ACF PACF
## [1,] 0.46 0.46
## [2,] -0.06 -0.36
## [3,] -0.05 0.23
## [4,] -0.05 -0.23
## [5,] -0.06 0.13
## [6,] -0.05 -0.15
## [7,] -0.07 0.03
## [8,] -0.08 -0.10
## [9,] -0.02 0.07
## [10,] 0.05 0.00
## [11,] 0.13 0.15
## [12,] 0.10 -0.10
## [13,] -0.01 0.06
## [14,] 0.04 0.06
## [15,] 0.08 0.01
## [16,] 0.01 -0.01
## [17,] -0.05 -0.03
## [18,] -0.08 -0.05
## [19,] -0.09 -0.01
## [20,] -0.04 -0.01
## [21,] -0.04 -0.07
## [22,] -0.07 -0.04
## [23,] 0.00 0.07
## [24,] 0.04 -0.06
## [25,] 0.02 0.04
## [26,] 0.00 -0.10
## [27,] -0.04 0.02
## [28,] -0.03 -0.03
## [29,] -0.07 -0.09
## [30,] -0.10 -0.01
## [31,] -0.05 -0.01
## [32,] 0.01 0.04
## [33,] 0.01 -0.02Primer: pounds_nz podaci
U dole prikazanom kodu uporedicemo MA(1), AR(1) i ARMA(1,1) modele. Ovi modeli su poredjeni na osnovi AIC. Ispostavlja se se da je najbolji model ARMA(1,1), pa zatim model AR(1) i najlosiji model je MA(1).
MA(1), AR(1) I ARMA(1,1) modeli
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.arma <- arima(x.ts, order = c(1, 0, 1))
AIC(x.ma)## [1] -3.526895AIC(x.ar)## [1] -37.40417AIC(x.arma)## [1] -42.27357x.arma##
## Call:
## arima(x = x.ts, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.8925 0.5319 2.9597
## s.e. 0.0759 0.2021 0.2435
##
## sigma^2 estimated as 0.01505: log likelihood = 25.14, aic = -42.27acf(resid(x.arma))#reziduali imaju male autokorelacije
Kao sto vidimo na korelogramu reziduali ARMA(1,1) modela imaju male autokorelacije, sto se poklapa sa sa realizacijom belog suma i time se opravdava koriscenje ovog modela.
primer 2: cbs podaci
Razmotrimo sad seriju korisecnu na drugom casu vezbi, seriju o proizvodnji elektricne energije u Australiji. Vidimo da serija ima pozitivan trend i sezonsku komponentu. Dodatno i disprezija vremenom raste, sto nam sugerise da je potrebno izvrsiti log-transformaciju.
www <- "http://www.math.rs/p/files/1550577090-70-cbe.txt"
CBE <- read.table(www, header = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
plot(Elec.ts)
Time <- 1:length(Elec.ts)
Imth <- cycle(Elec.ts)
Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth))
acf(resid(Elec.lm))
Cini se da korelogram reziduala ima ciklus sa periodom od 12 meseci. Ovakve situacije mogu da se dogode kada se koristi nestacionarni model kada je pristuna sezonska komponenta. O ovome ce se diskutovati na nekim sledecim vezbama.
Da bismo pronasli najbolji ARMA model mozemo da isprobavamo razlicite vrednosti za p i q pa zatim da poredimo modele na osnovu AIC(uzimamo model koji ima najmanji AIC). Ili mozemo da koristimo for petlju gde cemo tako da isprobavamo vrednosti. Postavicemo da p i q idu najdalje do vrednosti 2, poredimo modele na osnovu AIC i cuvamo onaj koji ima najmanji AIC.
best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:2) for (j in 0:2) {
fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0, j)))
if (fit.aic < best.aic) {
best.order <- c(i, 0, j)
best.arma <- arima(resid(Elec.lm), order = best.order)
best.aic <- fit.aic }
}primer 3: Wave tank podaci
Date su razlike u vrednosti nivoa vode u sredini bunara koje su posledica simulatora talasa. Merenja su vrsena na svakih 0.1s i poslednji trenutak je 39.7s. Nas cilj je da podatke uklopimo u najbolji ARMA model. Pogedajmo seriju, acf i pacf.
www <- "http://www.math.rs/p/files/1551218509-70-wave.txt"
wave.dat <- read.table(www, header = T)
attach (wave.dat)
layout(1:3)
plot (as.ts(waveht), ylab = 'Wave height')
acf (waveht)
pacf (waveht) 
Grafik pacf nam sugerise da p treba da bude najmanje 2. Isprobavanjem razlicitih vrednosti za p i q ispostavlja se da najbolji model ARMA(p,q) se dobija kada se za p i q uzme vrednost 4. ACF i PACF reziduala modela se poklapaju sa time da su reziduali realizacija belog suma.
wave.arma <- arima(waveht, order = c(4,0,4))
acf (wave.arma$res[-(1:4)]) 
pacf (wave.arma$res[-(1:4)]) 
hist(wave.arma$res[-(1:4)], xlab='mm')