Погледајмо прво неке примере временских серија (све у пакету astsa):

plot(jj, main = "Johnson & Johnson", type = "c")
text(jj, labels = 1:4, col = 1:4)

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

plot(globtemp, main = "Global temperature deviations", type = "o")

Приметимо да серија има углавном позитиван тренд, видимо да тренд није увек позитиван. За разлику од Johnson & Johnson серије, ова серија нема сезонску комппоненту.

plot(sp500w, main = "S&P 500")

За разлику од претходне две серије, ова серија нема ни тренд ни сезонску компоненту. Заправо, уочавамо да ова серија нема никакво правило у понашању. Видимо да је дисперзија повремено велика. Ово је пример серије која се назива шум.

Временска серија globtemp није стационарна, као сто смо видели, али диференцирање ове серије доводи до стационарности.

plot(globtemp, main = "globtemp")

plot(diff(globtemp), main = "diff(globtemp)")

У случају да постоји и тренд и хетероскедастичност, као код Johnson&Johnson серије, логаритмовање, па диференцирање даће нам стационарну временску серију. Прво логаритмовање стабилизује дисперзију, а затим диференцирање уклања тренд.

plot(jj, main="J&J")

plot(log(jj), main="log(J&J)")

plot(diff(log(jj)), main="diff(log(J&J))")

\(ARMA\) модели

ARMA(p,q) процеси су стационарни процеси који се могу записати у облику \[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}\]

где су \(\phi_{p}\neq 0, \theta_{q}\neq0\), а \(\omega_{t}\) је бели шум.

Вреднодт \(\mu\) представља средњу вредност процеса и ако се уведе ознака \(\alpha=\mu(1-\phi_{1}-...-\phi_{p})\) претходна једнакост се своди на: \[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\) и \(MA\) модели су специјани случајеви ARMA модела, па се ARMA модел може записати као: \[\Phi(B)x_{t}=\Theta(B)\omega_{t}.\]

Стационарност и инвертибилност

Код \(ARMA\) модела услов стационарности зависи од \(AR\) компоненте, а услов инвертибилности од \(MA\) компоненте.

Кажемо да је \(ARMA(p,q)\) модел стационаран ако се може записати као линеарни процес: \[x_{t}=\Psi(B)\omega_{t},\] где је \(\Psi(B)=\sum_{j=0}^{\infty}\psi_{j}B^{j}\) и \(\sum_{j=o}^{\infty}\mid\psi_{j}\mid<\infty\), \(\psi_{0}=1\).

Кажемо да је \(ARMA(p,q)\) модел инвертибилан ако се може записати као \[\omega_{t}=\pi(B)x_{t},\] где је \(\pi(B)=\sum_{j=0}^{\infty}\pi_{j}B^{j}\) и \(\sum_{j=o}^{\infty}\mid\pi_{j}\mid<\infty\), \(\pi_{0}=1\).

ARMA(1,1) \[X_t= \alpha X_{t-1}+\omega_t+\beta\omega_{t-1},\] где је \(\omega_t\) белли шум са очекивањем 0 и дисперзијом \(\sigma^2_\omega\).

Процес \(X_t\) може се записати преко белог шума, као: \[X_t= (1-\alpha B)^{-1}(1+\beta B) \omega_t = \dots = \omega_t+(\alpha+\beta)\sum_{i=1}^{\infty}\alpha^{i-1}\omega_{t-i}.\]

Јасно је да је очекивање процеса \(X_t\) једнако 0, a дисперзија процеса је \(\sigma_\omega^2(1+(\alpha+\beta)^2(1-\alpha^2)^{-1})\).

Аутоковаријациона функција је \[cov(X_t,X_{t+k}) = (\alpha+\beta)\alpha^{k-1}\sigma_\omega^2+(\beta+\alpha)^2\sigma_\omega^2 \alpha^k (1-\alpha^2)^{-1}.\] Аутокорелациона функција је \[\gamma_k = \frac{\alpha^{k-1}(\alpha+\beta)(1+\alpha\beta)}{1+\alpha\beta+\beta^2}.\]

Симулације \(ARMA\) модела

x=arima.sim(list(order=c(0, 0, 1), ma=0.9), n=100)
plot(x)

x=arima.sim(list(order=c(2, 0, 0), ar=c(0, -0.9)), n=100)
plot(x)

x=arima.sim(list(order=c(2, 0, 1), ar=c(0, -0.9), ma=0.9), n=100)
plot(x)

set.seed(1)
x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))
coef(arima(x, order = c(1, 0, 1)))
##          ar1          ma1    intercept 
## -0.596966371  0.502703368 -0.006571345

Пример 1

Упоредимо \(MA(1)\), \(AR(1)\) и \(ARMA(1,1)\) моделе на основу \(AIC\).

www <- "pounds_nz.txt"
x <- read.table(www, header = TRUE)
x.ts <- ts(x, st = 1991, fr = 4)
x.ma <- arima(x.ts, order = c(0, 0, 1))
x.ma
## 
## Call:
## arima(x = x.ts, order = c(0, 0, 1))
## 
## Coefficients:
##         ma1  intercept
##       1.000     2.8329
## s.e.  0.072     0.0646
## 
## sigma^2 estimated as 0.04172:  log likelihood = 4.76,  aic = -3.53
AIC(x.ma)
## [1] -3.526895
x.ar <- arima(x.ts, order = c(1, 0, 0))
x.ar
## 
## Call:
## arima(x = x.ts, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.9439     3.0111
## s.e.  0.0458     0.2953
## 
## sigma^2 estimated as 0.01818:  log likelihood = 21.7,  aic = -37.4
AIC(x.ar)
## [1] -37.40417
x.arma <- arima(x.ts, order = c(1, 0, 1))
x.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.27
AIC(x.arma)
## [1] -42.27357
acf(resid(x.arma))

На основу вредности \(AIC\) видимо да је \(ARMA(1,1)\) најбољи од понуђених модела.

Пример 2

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

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

Корелограм резидуала има циклус са периодом од 12 месеци. Овакве ситуације се могу јавити када се користе нестационарни модели (наредни часови).

Да бисмо пронашли најбољи \(ARMA\) модел можемо да испрбамо различите вредности за \(p\) и \(q\), па затим да поредимо моделе на основу \(AIC\) (узимамо модел који има најмањи \(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 }
}
best.order
## [1] 2 0 0
acf(resid(best.arma))

Пример 3

Дате су разлике у вредностима нивоа воде у средини бунара који су последица симулатора таласа. Мерења су вршена на сваких 0.1s и последњи тренутак је 39.7s. Наш циљ је серију уклопити у најбољи \(ARMA\) модел.

www <- "wave.txt" 
wave.dat <- read.table(www, header = T) 
attach (wave.dat) 
plot (as.ts(waveht), ylab = 'Wave height') 

acf (waveht) 

pacf (waveht)

График \(pacf\) нам сугерише да \(p\) треба да буде најмање 2. Испробавањем различитих вредности за p и q испоставља се да најбољи модел \(ARMA(p,q)\) се добија када се за \(p\) и \(q\) узме вредност 4. Са графика \(ACF\) и \(PACF\) резидуала модела закључујемо да су резидуали реализација белог шума.

best.order <- c(0, 0, 0)
best.aic <- Inf
for (i in 0:5) for (j in 0:5) {
  fit.aic <- arima(waveht, order = c(i,0,j))$aic
  if (fit.aic < best.aic) {
    best.order <- c(i, 0, j)
    best.arma <- arima(waveht, order = best.order)
    best.aic <- fit.aic
  }
}
## Warning in log(s2): NaNs produced
best.order
## [1] 4 0 4
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')

Задатак