Подсетимо се да постоје строго и слабо стационарни модели. Овде се бавимо слабо стационарним моделима, тј. моделима који имају константно очекивање и дисперзију у времену и \(cov(x_{t},x_{t+k})\) зависи само од кашњења \(|k|\). Стационарност је идеализација модела и ми увек настојимо да “очистимо” податке од трендова и сезонских компоненти да бисмо могли да им доделимо стационарни модел.

\(MA(q)\) модел

\(MA\) модел реда \(q\) је модел који се може записати у облику:

\[x_{t}= e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}-\dots-\theta_{q}e_{t-q},\] где су \(e_{t}\) iid грешке модела очекивања 0 и константне дисперзије тј. бели шум. Најчешће је то Гаусов бели шум. У \(MA\) моделима тренутна вредност серије зависи од претходних грешака. У некој литератури може се пронаћи и овај облик модела: \[x_{t}= e_{t}+\theta_{1}e_{t-1}+\theta_{2}e_{t-2}+\dots+\theta_{q}e_{t-q},\] али видимо да је то исто само шти коефицијенти тада мењају знак.
Неке особине:
1. Средња вредност: \(E(x_{t})=0\)
2. Дисперзија: \(D(x_{t})=(1+\theta_{1}^2+\dots+\theta_{q}^2)\sigma^2\)
3. Аутокорелациона функција: \[\rho(k) = \begin{cases}1, \; k=0\\ \frac{\sum_{i=0}^{q-k}\theta_i\theta_{i+k}}{\sum_{i=0}^q \theta_i^2}, \; k=1, \dots,q\\ 0, \; k>q\end{cases}\]

Приметимо да је \(MA(q)\) специјалан случај линеарног модела када узмемо да су коефицијенти у том моделу \(\theta_{j}=0\) за \(j>q\).

\(MA(q)\) модел је увек стационаран (за свако \(q\) и све коефицијенте \(\theta\)). Објашњење за то је некорелисаност \(x_{t}\) и свих пре \(x_{t-q}\), јер се не појављују исти шумови.

Линеарни модел је заправо MA(\(\infty\)): \[x_{t}=\displaystyle\sum_{i=0}^{\infty} \theta_{i}e_{t-i},\] где ред \(\displaystyle\sum_{i=0}^{\infty} \theta_{i}\) апсолутно конвергира.

Подсетимо се \(AR\) модела. Ауторегресиона серија реда \(q\) (ознакa \(AR(q)\)) је облика: \[x_{t}= e_{t}+\theta_{1}x_{t-1}-\theta_{2}x_{t-2}-\dots-\theta_{q}x_{t-q}\]

Видимо да њена тренутна вреност зависи од њених претходних вредности. Сада ћемо видети да су те зависности (од грешака и вредности модела) на неки начин исте. Следећа теорема нам даје веома важно својство стационарних процеса.

Волдова теорема о декомпозицији Сваки стационарни модел се може приказати у облику \(MA(\infty)\) модела.

Закључујемо да се сваки стационарни \(AR(q)\) модел може конвертовати у \(MA(\infty)\) и обрнуто, сваки \(MA(q)\) у \(AR(\infty)\) модел ако важи услов инвертибилности. \(MA(q)\) је инвертибилан и има јединствен инвертибилнии модел ако су нуле полинома са коефицијентима \(\theta\) ван јединичног корена (Више о овоме теоријски на предавањима).

Не постоји правило, а одговор на ово питање налазимо у искуству. Свакако што се тиче стационарних процеса ако нам је један модел одговарајући, онда ће и други бити скоро сигурно.

Посматрајмо сада најједноставнији \(MA(1)\) модел: \[x_{t}= e_{t}-\theta_{1}e_{t-1}\] ACF овог модела дата са: \[\rho_{0}=1,\quad \rho_{1}=\frac{\theta_{1}}{1+\theta_{1}^2}, \quad \rho_{j}=0, \; j\geq2. \] ACF узима исте вредности за \(\theta\) и \(\frac{1}{\theta}\) (види услов инвертибилности који нам каже који модел користимо).

Теоријски ACF ћемо израчунати у R-у:

rho <- function(lag, theta){#lag je kašnjenje u ACF, theta koeficijenti modela sa dodatom jedinicom
  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)
  }

Нацртајмо сада добијене теоријске вредности ACF за један \(MA(3)\) модел са параметрима \(\beta_1=0.7\), \(\beta_2=0.5\) и \(\beta_3=0.2\) и за један \(MA(3)\) модел са параметрима \(\beta_1=-0.7\), \(\beta_2=0.5\) и \(\beta_3=-0.2\):

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)}
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()

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)}
rho_k
##  [1] -0.6460674  0.3595506 -0.1123596  0.0000000  0.0000000  0.0000000
##  [7]  0.0000000  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()

Примери

Симулирајмо један \(MA(3)\) модел.

set.seed(1)
b <- c(0.8, 0.6, 0.4)
x <- w <- rnorm(1000)
for (t in 4:1000){
  for(j in 1:3) x[t] <- x[t] + b[j] * w[t - j]
}
plot(x, type = 'l')

acf(x)

Функција у R-у која симулира \(MA\) модел је arima.sim. Њени обавезни елементи су листа са елементима:
1. order = (0,0,q), где је q ред MA модела
2. maкоефицијенти \(\theta\) у моделу
3. n представља број серија који желимо да симулирамо.

Генеришимо сада два \(MA(1)\) модела:

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

Нацртајмо сада графике узорачких (оцењених) функција аутокорелације.

acf(ma1,lag.max = 50)

acf(ma2,lag.max = 50)

Плава испрекидана линија нам говори да су већ од \(lag=2\) вредности ван критичне области, па за њих прихватамо хипотезу да је аутокорелацја нула што је очекивано, јер су ово MA(1) модели. Видимо да има неких вредности за \(lag\geq2\) које су ван плавог испрекиданог појаса, али ако оне само мало излазе то занемарујемо, јер тада би опет прихватили хипотезу ако само мало повећамо праг значајности.

Упоредимо сада са теоријским вредностима за 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

Можемо приметити да је са повећањем обима узорка узорачка ACF ближа теоријској.

Фитовање МА модела

Функција arima се користи за фитовање \(MA(q)\) модела за неке податке тако што је параметар order=(0,0,q) где је \(q\) ред модела. За разлику од функције ar функција arima иницијално не одузима средњу вредности, па ћемо имати и оцењени intercept. Једна од повратних вредности су коефицијенти који се оцењују методом најмањих квадрата. У R-у се нумеричким методама тражи вредност коефицијената које минимизују:

\[ 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},\]

где су оцењене грешке заправо резидуали који се добијали итеративном методом, а за њихову почетну вредност узете су нуле.

x.ma <- arima(x, order = c(0, 0, 3))
x.ma
## 
## Call:
## arima(x = x, order = c(0, 0, 3))
## 
## Coefficients:
##          ma1     ma2     ma3  intercept
##       0.7898  0.5665  0.3959    -0.0322
## s.e.  0.0307  0.0351  0.0320     0.0898
## 
## sigma^2 estimated as 1.068:  log likelihood = -1452.41,  aic = 2914.83

За \(MA(q)\) модел теоријске вредности ACF су једнаке 0 за задршке веће од \(q\). Зато би узорачке требало да буду блиске 0. Ако серију уклапамо у \(MA\) модел, има смисла за ред модела изабрати ред оне узорачке аутокорелације која је последња значајна. Некад и код серија тачно симулираних из неког \(MA(q)\) модела, имаћемо значајне корелације и реда већег од \(q\). Битно је и колико те вредности испадају из граница. Ако је само мало, онда са мало већим нивоом значајности вредности узорачке ACF упадају у област прихватања хипотезе.

Посматрајмо сада базу са прошлог часа - индекси потрошачких цена.

data<-read.table("CPIAUCSL.txt")
data.ts<-ts(data$V2,start = c(1947,1), end = c(2015,3), frequency = 12 )
plot(data.ts)

acf(data.ts,lag.max = 30)

Када нацртамо корелограм резидуала видимо да модел није баш задовољавајући, јер не добијамо бели шум, али није ни тако лош што је последица тога да смо га ми на прошлом часу фитовали линеарно тј. за MA(\(\infty\)).

x.ts <- arima(data.ts, order = c(0, 0, 10)) # isprobati koji red najvise odgovara
## Warning in arima(data.ts, order = c(0, 0, 10)): possible convergence problem:
## optim gave code = 1
x.ts
## 
## Call:
## arima(x = data.ts, order = c(0, 0, 10))
## 
## Coefficients:
##          ma1     ma2      ma3      ma4      ma5      ma6     ma7     ma8
##       3.3654  6.8150  10.4034  12.9768  13.5686  11.9834  8.8217  5.1968
## s.e.  0.0320  0.0972   0.1847   0.2643   0.3071   0.2963  0.2377  0.1544
##          ma9    ma10  intercept
##       2.2301  0.5169   100.6168
## s.e.  0.0767  0.0256     3.2114
## 
## sigma^2 estimated as 1.446:  log likelihood = -1322.62,  aic = 2669.25

Напомена: У R-у постоји функција ma која изравњава временску серију методом покретних просека. Не треба мешати са моделовањем и функцијом arima!

Задатак