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

Напомене:

  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.

Ауторегресиони модели

\(AR(1)\) модел

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

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

  • Ако је \(|\varphi_1|<1\), овај процес се може приказати као линеарaн процес, тј. тада је \(AR(1)\) процес стационаран.

  • Записивање \(AR(1)\) процеса као линеарног процеса, итерацијама уназад \(x_t=\sum_{i=0}^{\infty}\varphi_1^iw_{t-i}\)

  • Члан \(\varphi_0\) има везе са средњом вредношћу, али није средња вредност. Добија се да је средња вредност \(\mu=\frac{\varphi_0}{1-\varphi_1}\).

  • Аутоковаријациона функција \(\gamma_k= \frac{\sigma_w^2 \varphi_1^k}{(1-\varphi_1^2)}, k\ge 0\)

  • Аутокорелациона функција \(\rho_k=\varphi_1^k, k\ge 0\)

Напомена:

  • Услов \(|\varphi_1|<1\) се у литератури некад назива услов стационарности, некад услов каузалности. Вероватно је мање правилно израз ‘услов стационарности’, јер и када је \(|\varphi_1|>1\), процес \(x_t\) је стационаран, може се записати у облику \(x_t=-\sum_{j=1}^{+\infty}\varphi_1^{-j}w_{t+j}\), само што је овакав процес неупотребљив, јер зависи од будућности. Дакле, не може се користити за предвиђање. Израз ‘услов каузалности’ ће односи на то да претходне вредности узрокују (causallity (енг.) = узрочност) садашњу, док у случају \(|\varphi_1|>1\) процес зависи од будућности. Ми ћемо користити израз услов каузалности.

Симулација једног \(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)))

\(AR(2)\) модел

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

\(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)\) процеса у зависности од тога да ли су нуле полинома реалне или коњуговано-комплексне изгледа другачије:

    • ако су реалне, ACF је линеарна комбинација две функције које опадају експоненцијалном брзином \(\rho(h)=c_1z_1^{-h}+c_2z_2^{-h}\)
    • ако нуле нису реалне, онда је \(\rho(h)=a|z_1|^{-h}cos(h\theta +b)\), где су a и b неке константе и овде функција експоненцијално опада, али има и облик синусоиде

Илустровање теоријског резултата 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)

Идентификовање реда \(AR(p)\) процеса

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

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)\) моделе.

YULE-WALKER

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

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

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

MLE

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

set.seed(1)
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.6009459

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

x.ar$ar+c(-2,2)*sqrt(x.ar$asy.var)
## Warning in c(-2, 2) * sqrt(x.ar$asy.var): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## [1] 0.4404031 0.7614886

Још један начин оцењивања реда модела је информациони критеријум 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)
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

Пример 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
## Warning in c(-2, 2) * sqrt(Z.ar$asy.var): Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## [1] 0.7405097 1.0400123
acf(Z.ar$res[-1]) #prvi rezidual (greska predvidjanja jedna korak unapred) ne postoji, jer ne postoji predvidjanje prvog elementa

  • Функција \(ar\) прво одузима средњу вредност, а уклапа у модел, тј. оцењује параметре. Како изгледа једначина одабраног модела?

\[\bar{z_t}=2.8+0.89(z_{t-1}-2.8)\]

Пример 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)

  • Како изгледа једначина одабраног модела? \[\bar{x}_t = -0.14+0.59(x_{t-1}+0.14)+0.013(x_{t-2}+0.14)+0.11(x_{t-3}+0.14)+0.27(x_{t-4}+0.14)\]

Задатак

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

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