Линеарни процес са средњом вредношћу \(\mu\) је процес који се може записати у облику: \[X_t=\mu + \sum_{j=0}^{\infty}\psi_jZ_{t-j},\] где је \(Z_t\) бели шум и \(\sum_{j=0}^{\infty}|\psi_j|<\infty.\)

Напомене:

  1. При наведеним претпоставкама ови процеси су стационарни.

  2. За линеарне процесе може се извести коваријациона функција као функција коефицијената \(\psi_j\).

ВОЛДОВА (Wold) ТЕОРЕМА О ДЕКОМПОЗИЦИЈИ Сваки слабо стационарни процес са средњом вредношћу 0 може се записати као \[X_t=\sum_{j=0}^{\infty}\psi_jZ_{t-j}+V_t,\] где је \(Z_t\) iid са очекивањем 0 и дисперзијом \(\sigma^2\), \(\sum_{j=0}^{\infty}\psi_j^2<\infty\). \(V_t\) је детерминистички процес, тј. предвидив на основу своје прошлости, функција је прошлих тренутака. Процеси \(V_t\) и \(Z_t\) су независни.

Напомене:

  1. Ако је процес Гаусов, онда случајне величине \(Z_t\) имају нормалну расподелу.

  2. Први члан збира је управо линеарни процес са средњом вредношћу 0.

\(\color{deepskyblue}{\text{$AR(1)$ модел}}\)

\[x_{t}=\varphi _0+\varphi _{1}x_{t-1}+w_{t}\] \[\left (1-\varphi _{1}L\right )x_{t}=w_{t}\]

Теоријски резултати:

Напомена:

Симулација једног \(AR(1)\) процеса

par(mfrow=c(1,1))
set.seed(3)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
plot(x, type = "l")

Како изгледа корелограм \(AR(1)\) процеса?

acf(x)

На корелограму се види експоненцијални пад. Ово је график узорачке аутокорелационе функције. Њега можемо добити када имамо задат узорак, а рачуна се по формулама за оцену аутокорелационе функције (формуле са другог часа).

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

Наредна два графика нису везана ни за једну временску серију (узорак), само за процес \(AR(1)\).

rho <- function(k, alpha) alpha^k
par(mfrow=c(1,2))
plot(0:10, rho(0:10, 0.7), type = "h")
abline(h=0)
plot(0:10, rho(0:10, -0.7), type = "h")
abline(h=0)

Наредна два графика су везана за узорак. Користимо функцију за симулирање \(AR(1)\) процеса.

plot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x",
     main=(expression(AR(1)~~~phi==+.9)))

plot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100), ylab="x",
     main=(expression(AR(1)~~~phi==-.9)))

\(\color{deepskyblue}{\text{$AR(2)$ процеси}}\)

\[x_t=\varphi_0+\varphi_1x_{t-1}+\varphi_2x_{t-2}+w_t\]

\(\color{deepskyblue}{\text{$AR(p)$ процеси}}\)

Ауторегресиони модел реда \(p\) \(AR(p)\) је: \[x_{t}=\varphi _{1}x_{t-1}+\varphi _{2}x_{t-2}+\varphi _{3}x_{t-3}+...+\varphi _{p}x_{t-p}+w_{t}\] где је \(w_{t}\) бели шум, \(\varphi _{i}\) параметри модела и \(\varphi _{p}\neq 0\) ако је модел реда \(p\).

Еквивалентна дефиниција преко оператора \(L\): \[\left ( 1-\varphi _{1}L+\varphi _{2}L^{2}+..+\varphi _{p}L^{p} \right )x_{t}=w_{t},\]

где је \(Lx_t=x_{t-1}\) backshift оператор (Lag оператор), који враћа уназад.

Сада \(AR(p)\) процесе можемо да запишемо као (гледамо оне са средњом вредношћу 0) \(\varphi(L)x_t=w_t\), где је \(\varphi(L)\) ауторегресиони оператор који дефинишемо као \[\varphi(L)=1-\varphi_1L-\varphi_2L^2-...-\varphi_pL^p.\]

Теоријски резултати:

  1. \(AR(p)\) процес је каузалан ако полином \(\varphi(L)\) има корене који су сви по норми строго већи од \(1\). (Погледати предавања!)

  2. Аутокорелациона функција \(AR(2)\) процеса у зависности од тога да ли су нуле полинома реалне или коњуговано-комплексне изгледа другачије:

Илустровање теоријског резултата 2:

z=c(1,-1.5,0.75) 
(a=polyroot(z)) #vraća nule polinoma koji je zadat vektorom koeficijenata, na i-tom mestu je u vektoru koeficijent uz x^{i-1}
## [1] 1+0.57735i 1-0.57735i
par(mfrow=c(1,1))
r<-arima.sim(list(order=c(2,0,0), ar=c(1.5,-0.75)), n=100)
acf(r)

Компликовано је рачунање аутокорелационе функције за \(AR(p)\) моделе већег реда.

У \(R\)-у постоји функција која рачуна за прослеђене параметре теоријску аутокорелациону функцију.

ACF=ARMAacf(ar=c(1.5,-0.75), ma=0, 50)
plot(ACF,type="h", xlab="lag")
abline(h=0)

\(\color{deepskyblue}{\text{Идентификовање реда $AR(p)$ процеса}}\)

У пракси ред \(p\) \(AR\) модела није познат, требало би га емпиријски одредити и онда неком методом одредити непознате параметре модела.

\(\color{deepskyblue}{\text{PACF - Парцијална аутокорелациона функција}}\)

Иако код \(AR(1)\) модела \(x_t\) може да се изрази само преко \(x_{t-1}\) корелација између \(x_t\) и \(x_{t-2}\) није једнака нули, већ \(x_t\) зависи од \(x_{t-2}\) преко \(x_{t-1}\).

Желимо неку функцију која мери зависност између \(x_t\) и \(x_s\) при уклођеном утицају свега између.

Парцијална корелација - мери линеарну везу између \(x_t\) и \(x_s\) при уклоњеном линеарном ефекту свега између, тј. сматра да је њихов утицај познат. Ово је било неформално објашњење. Формаалну дефиницију парцијалну аутокорелациону функцију за слабо стационарне процесе нећемо наводити.

  • Код \(AR\) модела ред модела је управо последњи члан парцијалне аутокорелационе функције који није нула.

По овој дефиницији могуће је једноставно израчунати на пример PACF за \(AR(1)\) модел. За моделе већег реда иде мало теже, а ово нам свакако не сугерише метод за оцењивање ове функције.

Једна добра оцена за PACF(h) може се добити на следећи начин:

  • Направити линеарни модел, тј. оценити параметре \(\varphi_0,...,\varphi_h\) методом најмањих квадрата \(x_t=\varphi_0+\varphi_1x_{t-1}+...+\varphi_hx_{t-h}+e_t\), па је \(\varphi_h\) оцена за PACF(h).

Други начин за рачунање исте оцене је преко Yule-Walker система једначина везан за \(AR(p)\) моделе.

\[ \begin{bmatrix} 1&\rho_1&\rho_2&\dots &\rho_{p-1}\\ \rho_1&1&\rho_1& \dots &\rho_{p-2} \\ \rho_2&\rho_1&1& \dots &\rho_{p-3}\\ \dots&\dots&\dots& \dots &\dots\\ \rho_{p-1}&\rho_{p-2}&\rho_{p-3}&\dots &1\\ \end{bmatrix} \begin{bmatrix} \varphi_1\\ \varphi_2\\ \varphi_3\\ \dots\\ \varphi_{p}\\ \end{bmatrix} = \begin{bmatrix} \rho_1\\ \rho_2\\ \rho_3\\ \dots\\ \rho_{p}\\ \end{bmatrix} \]

Парцијална аутокорелациона функција:

our_pacf<-function(x)
{
  lag_max=length(x)/5
  pacf_array<-c() 
  r=acf(x, lag.max = lag_max, plot = FALSE)$acf
  for(i in 1:lag_max)
  {
    ri=r[2:(i+1)] #niz ocenjenih autokorelacija
    Ri=toeplitz(r[1:i]) #Teplicova matrica u R-u, dovoljno je zadati prvu kolonu
    Phi<-solve(Ri, ri)
    pacf_array<-c(pacf_array, Phi[i])
  }
  return(pacf_array)
}

Одредимо парцијалну аутокорелациону функцију са кораком \(5\).

sacf <- acf(x, lag.max = 10, plot = FALSE)$acf #samo niz autokorelacija 
#lag.max - poredicemo samo prvih 5 parcijalnih autokorelacija
res1 <- solve(toeplitz(sacf[1:5]), sacf[2:6])
res1
## [1]  0.709330142 -0.019633857  0.072123282 -0.021849294 -0.004276498
res2 <- pacf(x, lag.max = 5, plot = FALSE)$acf
res2
## , , 1
## 
##              [,1]
## [1,]  0.723049621
## [2,]  0.019252410
## [3,]  0.054589317
## [4,] -0.024883198
## [5,] -0.004276498

Парцијална акутокорелација са корак \(5\) је последња у овом низу. Упоредимо да ли се исто добиија и применом функције \(pacf\).

pacf(x)

\(\color{deepskyblue}{\text{YULE-WALKER}}\)

Користи се ако нам је претходно познато да је у питању ауторегресиони модел.

Yule_Walker<-function(x,p)
{
  y<-acf(x,lag.max=p,plot="false")$acf
  vector<-solve(toeplitz(y[1:p]),y[2:(p+1)])
  v<-c(1,-vector)
  sigma2<-crossprod(acf(x,lag.max=p,plot="false",type="cov")$acf, v)
  return (c(vector,sigma2))
}
Yule_Walker(x,4)
## [1]  0.70943655 -0.01994266  0.07220857 -0.02488320  0.72302410
ar(x,method="yule-walker") 
## 
## Call:
## ar(x = x, method = "yule-walker")
## 
## Coefficients:
##     1  
## 0.723  
## 
## Order selected 1  sigma^2 estimated as  0.7407

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

Временску серију коју смо симулирали раније, тражимо најбољи AR модел “mle” методом.

set.seed(3)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
x.ar <- ar(x,method = "mle")
x.ar$order
## [1] 1
x.ar$ar
## [1] 0.7213527
x.ar$ar+c(-2,2)*sqrt(x.ar$asy.var)
## [1] 0.5837605 0.8589450

Расподела оцена AR модела добијених методама YW и MLE је асимптотски нормална, па се може одредити 95% интервал поверења за параметре.

Још један начин оцењивања реда модела је информациони критеријум AIC - Akaike Information Criterion. AIC је број који се рачуна на основу узорка (једне серије) и предвиђеног модела. Дакле, зависи од узорка и од модела. То је мера колико је неки модел добар за добијену серију. \[AIC=\frac{-2}{n}ln(likelihood)+\frac{2}{n}(broj \quad parametara)\]

Напомена: функција веродостојности (likelihood) је n-димензиона густина узорка. У статистици, због претпоставке о ПСУ функција веродостојности једнака је производу густина. У теорији временских серија, немамо ту претпоставку о независности, па је то само \(n\)-димензиона густина. Ред \(AR\) модела бира се рачунањем AIC, пошто се израз за AIC у случају \(AR\) модела поједностављује \[AIC(l)=ln(\hat{\sigma}_l^2)+\frac{2l}{n},\]

где је \(\hat{\sigma}_l^2\) оцена за \(\sigma_w^2\) методом максималне веродостојности.

set.seed(4)
polyroot(c(1,0.4,0.6,0.6))
## [1]  0.193159+1.079312i -1.386318-0.000000i  0.193159-1.079312i
x<-arima.sim(list(ar=c(-0.4,-0.6,-0.8)),n=200)
plot(x)

acf(x,lag.max=50) #ne vidi se nista, neka kombinacija opadanja u sinusoidnom stilu

pacf(x) #vidi se ono sto zelimo! Korelacija reda 3 je poslednja znacajna

set.seed(2)
x<-arima.sim(list( ar=c(0.5,0.1,0.3,-0.3)), n=300)
pacf(x)$acf #rekli bismo da je red 4

## , , 1
## 
##               [,1]
##  [1,]  0.552149438
##  [2,]  0.130807401
##  [3,]  0.088130014
##  [4,] -0.356604564
##  [5,]  0.001010459
##  [6,]  0.013748379
##  [7,] -0.114343815
##  [8,] -0.023198900
##  [9,] -0.046257055
## [10,]  0.016273162
## [11,]  0.038294132
## [12,] -0.027601487
## [13,]  0.144193814
## [14,]  0.066313896
## [15,] -0.082340112
## [16,]  0.034246833
## [17,] -0.055189599
## [18,] -0.002776617
## [19,] -0.044540505
## [20,] -0.047537435
## [21,] -0.045679215
## [22,] -0.022083239
## [23,] -0.078656569
## [24,]  0.086631745

\(\color{deepskyblue}{\text{Пример 1}}\)

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
#srednju vrednost i ne ocenjuje taj parametar nijednom metodom
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 za koeficijente
## [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\) прво одузима средњу вредност, а уклапа у модел, тј. оцењује параметре. Како изгледа једначина одабраног модела?
mean(Z.ts)
## [1] 2.823251
plot(Z.ts)
lines(Z.ar$resid+Z.ts, col="red") #na grafik dodajemo predvidjene vrednosti

\(\color{deepskyblue}{\text{Пример 2}}\)

www<-"global.txt"
Global<-scan(www) #ucita u niz
Global.ts<-ts(Global,frequency=12,start=c(1856,1),end=c(2005,12))
plot(Global.ts)

Global.ar <- ar(aggregate(Global.ts, FUN = mean),method="mle")
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1
mean(aggregate(Global.ts, FUN = mean))
## [1] -0.1382628
Global.ar$order
## [1] 4
Global.ar$ar
## [1] 0.58762026 0.01260254 0.11116731 0.26763656
acf(Global.ar$res[-(1:Global.ar$order)], lag = 50)

  • Како изгледа једначина одабраног модела?

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

  • Симулирати 1000 вредности временске серије која је дата моделом \(x_t = \frac{5}{6}x_{t-1}-\frac{1}{6}x_{t-2}+w_t\), где је \(w_t\) Гаусов бели шум.
  • Нацртати корелограм и парцијални корелограм симулиране серије.
  • Уклопити серију \(\{x_t\}\) у одговарајући \(AR\) модел и написати једначину одговарајућег \(AR\) модела.
  • Одредити 95\(\%\) интервал поверења за оцењене параметре.
  • Испитати стационарност (каузалност) временске серије.
  • Нацртати корелограм резидуала фитованог модела.

Све резултате и графике прокоментарисати.