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

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

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

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

Globtemp временска серија представља одступања од годишњих температура (из периода 1951-1980) у периоду 1880-2015. године.

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

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

S&P 500 временска серија

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

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

\(\color{deepskyblue}{\text{Регресиони модели временских серија}}\)

Регресија \(Y_{i}=\beta X_{i} + \epsilon_{i}\), при чему је \(\epsilon_{i}\) бели шум

Ауторегресија (AR): \(X_{t}=\phi X_{t-1}+\epsilon_{t}\) (\(\epsilon_{t}\) је бели шум)

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

Moving average (MA): \(\epsilon_{t}=W_{t}+\theta W_{t-1}\), (\(W_{t}\) је бели шум). Све претходно поменуто доводи нас до ARMA модела. Другим речима, модел ARMA је ауторегресија са аутокорелисаним грешкама.

ARMA \(X_{t}=\phi X_{t-1}+W_{t}+\theta W_{t-1}\)

\(\color{deepskyblue}{\text{Стационарност}}\)

У контексту временских серија стационарност представља константност средње вредности. Kажемо да је временска серија стабилна ако је средња вредност константа (односно да нема тренда) и да је корелациона структура података из серије константна током времена.

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

par(mfrow=c(2,1))
plot(globtemp, main = "globtemp")
plot(diff(globtemp), main = "diff(globtemp)")

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

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

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(p,q) модел каузалан ако се може записати као линеарни процес: \(x_{t}=\Psi(B)\omega_{t}\), gde je \(\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}\), gde je \(\pi(B)=\sum_{j=0}^{\infty}\pi_{j}B^{j}\) и \(\sum_{j=o}^{\infty}\mid\pi_{j}\mid<\infty\), \(\pi_{0}=1\).

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)

\(\color{deepskyblue}{\text{Како препознати у који ARMA модел да уклопимо серију?}}\)

Као што видимо на графицима, AR(2) и MA(1) изгледају врло слично. Јасно је да се само са графика временске серије не може препознати који је модел у питању. Морају се посматрати ACF и PACF функције.

Погледајмо шта се дешава у случају AR(2) и MA(1) модела.

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

##         ACF  PACF
##  [1,]  0.01  0.01
##  [2,] -0.92 -0.92
##  [3,] -0.04  0.03
##  [4,]  0.85 -0.01
##  [5,]  0.05  0.00
##  [6,] -0.78  0.02
##  [7,] -0.06  0.03
##  [8,]  0.71 -0.08
##  [9,]  0.07  0.04
## [10,] -0.63  0.02
## [11,] -0.09 -0.04
## [12,]  0.56 -0.01
## [13,]  0.09 -0.02
## [14,] -0.50  0.04
## [15,] -0.12 -0.13
## [16,]  0.44  0.03
## [17,]  0.13 -0.05
## [18,] -0.38  0.02
## [19,] -0.14  0.02
## [20,]  0.33  0.04
## [21,]  0.14 -0.01
## [22,] -0.29  0.01
## [23,] -0.15 -0.04
## [24,]  0.25 -0.03
## [25,]  0.15 -0.06
## [26,] -0.21  0.00
## [27,] -0.15 -0.02
## [28,]  0.17 -0.07
## [29,]  0.14  0.01
## [30,] -0.14 -0.01
## [31,] -0.14 -0.05
## [32,]  0.10 -0.06
## [33,]  0.14 -0.01
x=arima.sim(list(order=c(0, 0, 1), ma=0.9), n=500)
acf2(x)

##         ACF  PACF
##  [1,]  0.53  0.53
##  [2,]  0.07 -0.28
##  [3,]  0.11  0.31
##  [4,]  0.07 -0.24
##  [5,]  0.00  0.16
##  [6,] -0.02 -0.16
##  [7,] -0.07  0.02
##  [8,] -0.11 -0.09
##  [9,] -0.04  0.09
## [10,] -0.05 -0.14
## [11,] -0.11  0.03
## [12,] -0.14 -0.16
## [13,] -0.15  0.00
## [14,] -0.09 -0.01
## [15,]  0.00  0.05
## [16,]  0.03  0.01
## [17,]  0.06  0.08
## [18,]  0.03 -0.13
## [19,] -0.03  0.03
## [20,]  0.00 -0.04
## [21,]  0.01  0.00
## [22,]  0.00  0.01
## [23,]  0.00 -0.02
## [24,] -0.02 -0.04
## [25,] -0.04 -0.01
## [26,] -0.01 -0.01
## [27,]  0.03  0.06
## [28,]  0.08  0.10
## [29,]  0.11  0.03
## [30,]  0.05 -0.04
## [31,]  0.04  0.08
## [32,]  0.10 -0.03
## [33,]  0.04 -0.03

Симулација ARMA(1,1) процеса и уклапање у ARMA(1,1) модел

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

Упоредимо 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) најбољи од понуђених модела.

\(\color{deepskyblue}{\text{pounds_nz}}\)

www <- "pounds_nz.txt"
Z <- read.table(www, header = TRUE)
Z.ts <- ts(Z, st = 1991, fr = 4)
plot(Z.ts)

Z.ar <- ar(Z.ts)
Z.ar1<-ar(Z.ts-mean(Z.ts)) #vidi se da su ovo isti modeli, dakle funkcija ar prvo oduzima
#tako ocenjuje srednju vrednost
Z.ar$x.mean; mean(Z.ts)  # ove dve su jednake
##    xrate 
## 2.823251
## [1] 2.823251
Z.ar$order #nismo naveli red modela, ar funkcija je odredila najbolji na osnovu AIC
## [1] 1
Z.ar$ar
## [1] 0.890261
Z.ar$ar + c(-2, 2) * sqrt(Z.ar$asy.var) #95%-ni IP, asimptotski normalna raspodela AR parametara.
## [1] 0.7405097 1.0400123
#Prvi rezidual (greska predvidjanja jedna korak unapred) ne postoji
#jer ne postoji predvidjanje prvog elementa,
#zato ako ispitujemo korelogram reziduala
#radi provere da li je model dobijen fjom ar dobar moramo iskljuciti prvi clan
acf(Z.ar$res[-1])

Пошто функција ar прво одузима среднју вредност, па онда уклапа у модел, тј. оцењује параметре модела, вредности предвиђене моделом добијају се тако што се дода средња вредност. \[\hat{z}_t= 2.82+ 0.89*(z_{t-1}-2.82)\]

\(\color{deepskyblue}{\text{cbе}}\)

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

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

\(\color{deepskyblue}{\text{Wave tank}}\)

Дате су разлике у вредностима нвоа воде у средини бунара који су последица симулатора таласа. Мерења су вршена на сваких 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 резидуала модела се поклапају са тиме да су резидуали реализација белог шума.

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

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

\(\color{deepskyblue}{\text{Задатак}}\)