Modeli pokretnih proseka (moving average models)

Podsetimo se da postoje strogo i slabo stacionarni modeli. Mi se ovde bavimo slabo stacionarnim tj. modelima koji imaju konstantno očekivanje i disperziju u vremenu i \(cov(x_{t},x_{t+k})\) zavisi samo od kašnjenja \(|k|\). Stacionarnost je idealizacija modela i mi uvek nastojimo da “očistimo” podatke od trendova i sezonskih komponenti da bi mogli da im dodelimo stacionaran model.

MA model reda q (oznaka MA(q)) se može napisati u obliku:

\[x_{t}= e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}-\dots-\theta_{q}e_{t-q},\] gde su \(e_{t}\) iid greške modela tj. beli šum. Najčešće je to Gausov šum. U MA modelima trenutna vrednost serije zavisi od prethodnih grešaka. U nekoj literaturi možete pronaći i ovaj oblik modela: \[x_{t}= e_{t}+\theta_{1}e_{t-1}+\theta_{2}e_{t-2}+\dots+\theta_{q}e_{t-q},\] ali vidimo da je to isto samo što koeficijenti tada menjaju znak.
MA(q) model je uvek stacionaran (za svako q i sve koeficijente \(\theta\)). Objašnjenje za to je nekorelisanost \(x_{t}\) i svih pre \(x_{t-q}\) jer se ne pojavljuju isti šumovi.
Neke osobine:
1. E\(x_{t}\)=0
2. D\(x_{t}=(1+\theta_{1}^2+\dots+\theta_{q}^2)\sigma^2\)
Primetimo da je ovo specijalan slučaj linearnog modela kada uzmemo da su koeficijenti u tom modelu \(\theta_{j}=0\) za \(j>q\). Linearni model je zapravo MA(\(\infty\)): \[x_{t}=\displaystyle\sum_{i=0}^{\infty} \theta_{i}e_{t-i},\] gde red \(\displaystyle\sum_{i=0}^{\infty} \theta_{i}\) apsolutno konvergira.
Podsetimo se AR modela. Autoregresiona serija reda q (oznaka AR(q)) je oblika: \[x_{t}= e_{t}+\theta_{1}x_{t-1}-\theta_{2}x_{t-2}-\dots-\theta_{q}x_{t-q}\] Vidimo da njena trenutna vrednost zavisi od njenih prethodnih vrednosti. Sada ćemo videti da su te zavisnosti (od grešaka i vrednosti modela) na neki način iste. Sledeća teorema nam daje veoma važno svojstvo stacionarnih procesa.
Voldova teorema o dekompoziciji: Svaki stacionarni model se može prikazati u obliku MA(\(\infty\)) modela.
Zaključujemo da se svaki stacionarni AR(q) model može konvertovati u MA(\(\infty\)) i obrnuto, svaki MA(q) u AR(\(\infty\)) model ako važi uslov invertibilnosti. MA(q) je invertibilan i ima jedinstven invertibilni model ako su nule polinoma sa koeficijentima \(\theta\) van jediničnog kruga (više o ovome teorijski na predavanjima).
Sada se pitamo kada onda koristimo MA modele? Ne postoji pravilo, a odgovor na ovo pitanje nalazimo u iskustvu. Svakako što se tiče stacionarnih procesa ako nam je jedan model odgovarajući onda će i drugi biti skoro sigurno.
Posmatrajmo sada najjednostavniji MA(1) model: \[x_{t}= e_{t}-\theta_{1}e_{t-1}\] ACF ovog modela je data sa: \[\rho_{0}=1,\rho_{1}=\frac{\theta_{1}}{1+\theta_{1}^2}, \rho_{j}=0, |j|\geq2 \] Fun fact: ACF uzima iste vrednosti za \(\theta\) i \(\frac{1}{\theta}\) (vidi uslov invertibilnosti koji nam kaže koji model koristimo).
Formula za opšti MA(q) model je malo komplikovanija pa je ovde ne navodimo, ali je jasno da se analogno računa. Teorijski ACF ćemo izračunati u R-u:

rho <- function(lag, theta){#lag je kašnjenje u ACF, theta koeficijenti modela sa dodatim kecom
  q <- length(theta) - 1
  if (lag > q){
    ACF <- 0
    }else {
    s1 <- 0
    s2 <- 0
    for (i in 1:(q-lag+1)) {
      s1 <- s1 + theta[i] * theta[i+lag]}
    for (i in 1:(q+1)) {
      s2 <- s2 + theta[i]^2}
    ACF <- s1 / s2}
  return(ACF)}

Nacrtajmo sada dobijene teorijske vrednosti ACF za jedan MA(3) model:

theta <- c(1, 0.7, 0.5, 0.2) 
polyroot(c(1,0.7,0.5,0.2))     #sve nule su veće od 1
## [1] -0.25+1.561249i -2.00+0.000000i -0.25-1.561249i
rho_k <- rep(0, 10)            #ovde ćemo smestiti vrednosti funkcije
for (k in 1:10){
  rho_k[k] <- rho(k, theta)}
print(rho_k)
##  [1] 0.6460674 0.3595506 0.1123596 0.0000000 0.0000000 0.0000000 0.0000000
##  [8] 0.0000000 0.0000000 0.0000000
library(ggplot2)
data<-as.data.frame(cbind(lag=0:10, rho_k=c(1,rho_k)))
ggplot(data,aes(x=lag,y=rho_k))+
  geom_segment(aes(y=0, x=lag,xend=lag,yend=rho_k),col="red")+
  geom_point()

Primeri:
Funkcija u R-u koja simulira MA model je arima.sim. Njeni obavezni parametri su lista sa elementima:
1. order = (0,0,q), gde je q red MA modela
2. ma = koeficijenti \(\theta\) u modelu
i n koji predstavlja broj serija koji želimo da simuliramo.
Generišimo sada dva MA(1) modela:

set.seed(5)
ma1<-arima.sim(list(order=c(0,0,1), ma=0.5), n=100)
ma2<-arima.sim(list(order=c(0,0,1), ma=-0.5), n=100)
par(mfrow = c(1,2))
plot.ts(ma1,main="MA(1) sa koef. 0.5")
plot.ts(ma2,main="MA(1) sa koef. -0.5")

Nacrtajmo sada grafike uzoračkih (ocenjenih) funkcija autokorelacije.

acf(ma1,lag.max = 50)

acf(ma2,lag.max = 50)

Plava isprekidana linija nam govori da su već od \(lag=2\) vrednosti van kritične oblasti pa za njih prihvatamo hipotezu da je autokorelacija nula što je i očekivano jer su ovo MA(1) modeli. Vidimo i da ima nekih vrednosti za \(lag\geq2\) koje su van plavog isprekidanog pojasa, ali ako one samo malo izlaze to zanemarujemo jer tada bi opet prihvatili hipotezu ako samo malo povećamo prag značajnosti.
Uporedimo to sada sa teorijskim vrenostima za ACF:

theta1 <- c(1,0.5)   #ponavljamo kod od gore
rho_k1 <- rep(0, 50)            
for (k in 1:50){
  rho_k1[k] <- rho(k, theta1)}
theta2 <- c(1,-0.5) 
rho_k2 <- rep(0, 50)            
for (k in 1:50){
  rho_k2[k] <- rho(k, theta2)}
par(mfrow = c(1,2))
plot(0:50,c(1,rho_k1), type = "h", xlab = "lag", ylab = "ma1")
plot(0:50,c(1,rho_k2), type = "h",xlab = "lag", ylab = "ma2")

ma1<-arima.sim(list(order=c(0,0,1), ma=0.5), n=100)
dist(rbind(rho_k1, acf(ma1,lag.max = 50,plot = F)$acf[1:50]), method = "euclidian")
##    rho_k1
##  1.068124
ma1<-arima.sim(list(order=c(0,0,1), ma=0.5), n=1000)
dist(rbind(rho_k1, acf(ma1,lag.max = 50,plot = F)$acf[1:50]), method = "euclidian")
##     rho_k1
##  0.7399258

Možemo primetiti očekivano da je sa povećanjem obima uzorka uzoračka ACF bliža teorijskoj.
Fitovanje MA modela: Funkcija arima se koristi za fitovanje MA(q) modela za neke podatke tako što je, isto kao gore, parametar order=(0,0,q) gde je q red modela. Za razliku od funkcije ar funkcija arima inicijalno ne oduzima srednju vrednost pa ćemo imati i ocenjeni intercept. Jedna od povratnih vrednosti su koeficijenti koji se ocenjuju metodom najmanjih kvadrata. R numeričkim metodama traži vrednosti koeficijenata koje minimizuju:

\[ S(\hat\theta_{1},\dots,\hat\theta_{q})=\displaystyle\sum_{t=1}^{n} {\big(x_{t}-(\hat\theta_{1}\hat{e_{t-1}}+\dots+\hat\theta_{q}\hat{e_{t-q}})\big)^2},\]

gde su ocenjene greške zapravo reziduali koji se dobijaju iterativnom metodom, a za njihovu početnu vrednost uzete su nule.
Kako biramo red q?
Komentarisali smo već kada vrednosti autokorelacione funkcije upadaju u kritičnu oblast. Zato ima smisla da za red MA modela uzmemo indeks poslednje značajne autokorelacije (jer za ostale smatramo da su nula).
Posmatrajmo sada bazu sa prošlog časa-indeksi potrošačkih cena.

library(DT)
## Warning: package 'DT' was built under R version 3.5.3
data<-read.table("http://www.math.rs/p/files/70-CPIAUCSL.txt")
datatable(data)
data.ts<-ts(data$V2,start = c(1947,1), end = c(2015,3), frequency = 1 ) #uzećemo da je frekvencija 1
plot(data.ts)

acf(data.ts,lag.max = 30) #odavde čitamo q približno 10

#sada fitujemo 
data.ma<-arima(data.ts, order = c(0,0,10))
#da vidimo rezultate
data.ma
## 
## Call:
## arima(x = data.ts, order = c(0, 0, 10))
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       1.8471  2.3491  3.3165  3.5919  3.4901  3.3034  2.5159  1.3414
## s.e.  0.1551  0.2905  0.3690  0.4361  0.4892  0.5404  0.4862  0.3584
##          ma9    ma10  intercept
##       0.7633  0.4620    24.4165
## s.e.  0.2837  0.1531     0.3524
## 
## sigma^2 estimated as 0.01725:  log likelihood = 33.37,  aic = -42.74

Kada nacrtamo korelogram reziduala vidimo da model nije baš zadovoljavajući jer ne dobijamo beli šum, ali nije ni tako loš što je posledica toga da smo ga mi na prošlom času fitovali linearno tj.za MA(\(\infty\)):

acf(data.ma$res,lag.max = 30)

Napomena: U R-u postoji funkcija ma koja “izravnjava” vremensku seriju metodom pokretnih proseka. Ne treba mešati sa modelovanjem i funkcijom arima!