Autoregresioni modeli
Vladana Đakovic
15 March 2019
definicija: Serija je autoregresiona reda p (AR(p)) ako važi \[x_{t}=\varphi _{1}x_{t-1}+\varphi _{2}x_{t-2}+\varphi _{3}x_{t-3}+...+\varphi _{p}x_{t-p}+e_{t}\] gdje je \(e_{t}\) bijeli šum i \(\varphi _{i}\) parametri modela i \(\varphi _{p}\neq 0\) ako je reda p. Ekvivalentna definicja preko operatora L: \[\left ( 1-\varphi _{1}L+\varphi _{2}L^{2}+..+\varphi _{p}L^{p} \right )x_{t}=e_{t}\] Ako posmatramo slučaj \(p=1\) tj AR(1)model: \[x_{t}=\varphi _{1}x_{t-1}+e_{t}\] \[\left ( 1-\varphi _{1}L\right )x_{t}=e_{t}\] I slučaj \(p=2\) tj AR(2)model: \[x_{t}=\varphi _{1}x_{t-1}+\varphi _{2}x_{t-2}+e_{t}\] \[\left ( 1-\varphi _{1}L+\varphi _{2}L^{2}\right )x_{t}=e_{t}\]
Uslov stacionarnosti
To je uslov koji treba da zadovoljavaju koeficijenti AR modela da bi taj model mogao da se zapiše kao linearni proces. \[1-\varphi _{1}L+\varphi _{2}L^{2}+..+\varphi _{p}L^{p} =\left ( 1-g_{1}L \right ) \left ( 1-g_{2}L \right )...\left ( 1-g_{p}L \right )\] gdje su \(g_{1}...g_{p}\) rjesenja karakteristicne jednacine \[g^{p}-\varphi _{1}g^{p-1}-\varphi _{2}g^{p-2}-....-\varphi _{p}=0\] pa je uslov da su \(\left |g_{i} \right |< 1\) Ako vazi uslov stacionarnosti tj AR(p) model se može predstaviti kao linearni onda je $EX_{t} Ako posmatramo slučaj \(p=1\):
Autokorelacioni i autokovarijacioni koeficijent
Autokovarijacioni koeficijent AR(p)modela je \[ \gamma_{k}=\left\{\begin{matrix} \varphi _{1}\gamma_{k-1}+\varphi _{2}\gamma_{k-2}+...+\varphi _{p}\gamma_{k-p}+\sigma^{2} &k=0 \\ \varphi _{1}\gamma_{k-1}+\varphi _{2}\gamma_{k-2}+...+\varphi _{p}\gamma_{k-p} &k>0 \end{matrix}\right. \] Pa je autokorelacioni koeficijet AR(p) modela \[\rho_{k}=\rho _{1}\gamma_{k-1}+\rho _{2}\gamma_{k-2}+...+\rho _{p}\gamma_{k-p} k>0\]
Autokorelacioni i autokovarijacioni koeficijent AR(1)
\[\gamma_{k}=\sigma^{2}* \frac{\varphi _{1}^{k}} {1-\varphi _{1}^{2}}\] i \[\rho_{k}=\varphi _{1}^{k}\] ###### Autokovarijacioni koeficijent AR(2)
\[k=0: \gamma_{0}= \varphi _{1}\gamma_{1}+\varphi _{2}\gamma_{2}+\sigma^{2}\] \[ k=1: \gamma_{1}= \varphi _{1}\gamma_{0}+\varphi _{2}\gamma_{1}\] \[ k=2: \gamma_{2}= \varphi _{1}\gamma_{1}+\varphi _{2}\gamma_{2}\] odakle se dobijaju \(\gamma_{0},\gamma_{1},\gamma_{2}\)
Primjeri
AR(1) PRIMJERI
simulacija jednog AR(1) procesa
par(mfrow=c(1,1))
set.seed(3)
x <- e <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + e[t]
plot(x, type = "l")
acf(x)
vidi se eksponencijalni pad u ACF
Pažnja! Ovo je grafik UZORAČKE ACF. Njega možemo da dobijemo kada imamo uzorak. Računa se po formulama za ocenu ACF. Kada imamo zadat stacionaran proces za njega možemo računati TEORIJSKU ACF
Sledeća dva grafika nisu vezana ni za jednu seriju-uzorak, samo za proces 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)
Sledeća dva grafika vezana su za uzorak.Koristimo funkciju za simuliranje AR(1) procesa, mada smo već napravili jednu svoju
par(mfrow=c(2,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) PROCESI
Teorijski rezultati
AR(p) proces je kauzalan ako polinom (B) ima korene koji su svi po normi strogo veći od 1
ACF AR(2) procesa u zavisnosti od toga da li su nule polinoma realne ili konj kompleksne
izgleda drugačije - ako su realne, ACF je linearna kombinacija dve funkcije koje opadaju eksponencijalnom brzinom \(\rho(h)=c_1z_1^{-h}+c_2z_2^{-h}\) - ako nisu realne nule onda je \(\rho(h)=a|z_1|^{-h}cos(h\theta +b)\), gde su a i b neke konstante i ovde funkcija eksponencijalno opada ali ima i oblik sinusoide
z=c(1,-1.5,0.75)
(a=polyroot(z)) #vraća nule polinoma koji je zadat vektorom koeficijenata, na i-tom mestu je## [1] 1+0.57735i 1-0.57735i# u vektoru koeficijent uz x^{i-1}
par(mfrow=c(1,1))
r<-arima.sim(list(order=c(2,0,0), ar=c(1.5,-0.75)), n=100)
acf(r)
Komplikovano je računanje ACF za AR(p) modele većeg reda. U R-u postoji funkcija koja računa za prosleđene parametre teorijsku ACF
ACF=ARMAacf(ar=c(1.5,-0.75), ma=0, 50)
plot(ACF,type="h", xlab="lag")
abline(h=0)
Identifikovanje reda AR(p) procesa
U praksi red p AR modela nije poznat, trebalo bi ga empirijski odrediti i onda nekom metodom oceniti nepoznate parametre modela.
PACF - parcijalna autokorelaciona funkcija
Iako kod AR(1) modela \(x_t\) može da se izrazi samo preko \(x_{t-1}\) korelacija između \(x_t\) i \(x_{t-2}\) nije jednaka 0.\(x_t\) zavisi od \(x_{t-2}\) preko \(x_{t-1}\). Želimo neku funkciju koja meri zavisnost između \(x_t\) i \(x_s\) pri uklonjenom uticaju svega između. Parcijalna korelacija radi to - meri linearnu vezu između \(x_t\) i \(x_s\) pri uklonjenom linearnom efektu svega između, tj smatra da je njihov uticaj poznat.Ovo je bilo neformalno objašnjenje Formalnu definiciju PACF za slabo stacionarne procese nećemo navoditi
U specijalnom slučaju, kada se radi o Gausovim slabo stacionarnim procesima,formalno to je uslovna korelacija, korelacija uslovne raspodele vektora (\(x_t, x_{t+l}\)) pri uslovu poznatih vrednosti (\(x_{t+1},...,x_{t+l-1}\)). Za l=1 definišemo PACF(1)=ACF(1), a za l>1 kao što je navedeno. Kod AR modela red modela je upravo poslednji član PACF koji nije nula
Po ovoj definiciji moguće je jednostavno izračunati na primer PACF za AR(1) model za modele većeg reda vec ide malo teže, a ovo nam svakako ne sugerise metod za ocenjivanje ove funkcije.
Jedna dobra ocena za PACF(h) moze se dobiti na sledeći nacin 1. Napraviti linearni model, tj oceniti parametre \(\beta_0,...,\beta_h\)metodom najmanjih kvadrata \(r_t=\beta_0+\beta_1r_{t-1}+...+\beta_hr_{t-h}+e_t\) 2. \(\beta_h\) je ocena za PACF(h)
Drugi način za računanje iste ocene je preko Yule-Walker sistema jednačina Vezan je za AR(p) modele. Na kraju je naveden algoritam koji ćemo implementirati
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)
{
#niz ocenjenih autokorelacija
ri=r[2:(i+1)]
#R je tzv Teplicova matrica, dovoljno je zadati one koji ce biti u prvoj koloni
#sve ostalo se onda zna
Ri=toeplitz(r[1:i])
Phi<-solve(Ri, ri)
pacf_array<-c(pacf_array, Phi[i])
}
return(pacf_array)
}
set.seed(1)
x<-arima.sim(list( ar=c(0.5,0.1,0.3,-0.1)), n=10000)
#R javlja ako smo izabrali koeficijente td AR nije eksplozivan
#mozemo da proverimo eksplozivan nalazenjem nula polinoma
polyroot(c(1,0.5,0.1,0.3,-0.1))## [1] 0.159563+1.53047i -1.113153+0.00000i 0.159563-1.53047i
## [4] 3.794028-0.00000i#sve nule su vece od 1 po normi
#all.equal(our_pacf(x),pacf(x,plot=FALSE)$acf[1:(length(x)/5)])
pacf(x)
set.seed(1)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
plot(x, type = "l")
acf(x)
pacf(x)
prvi primjer sa pacf funkcijom
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")
acf(x)
#kako izgleda autokorelaciona funkcija za AR(1) model
rho <- function(k, alpha) alpha^k
par(mfrow=c(1,2))
plot(0:10, rho(0:10, 0.7), type = "b")
plot(0:10, rho(0:10, -0.7), type = "b")
#kako na grafik utice znak koeficijenta
par(mfrow=c(2,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)))
sacf <- acf(x, lag.max = 10, plot = FALSE)$acf #samo niz autokorelacija
#ako npr zelimo pacf sa korakom 5, resavamo sistem
#pravljenje M=Teplicove matrice, i vektora y, resava se sistem Mx=y
res1 <- solve(toeplitz(sacf[1:5]), sacf[2:6])
res1## [1] 0.709330142 -0.019633857 0.072123282 -0.021849294 -0.004276498#znaci pacf sa korakom 5 je poslednji u ovom nizu
#uporedimo da li se isto dobija i primenom funkcije pacf
#lag.max - poredicemo samo prvih 5 parcijalnih autokorelacija
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#jednaki!
pacf(x)
YULE-WALKER ako prethodno znamo da je AR(p) model
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.72302410ar(x,method="yule-walker") #isti koeficijenti, ali razlicito sigma (?)##
## Call:
## ar(x = x, method = "yule-walker")
##
## Coefficients:
## 1
## 0.723
##
## Order selected 1 sigma^2 estimated as 0.7407CSS - ZA AR(1)
set.seed(199)
x <- w <- rnorm(100)
for (t in 2:100) x[t] <- 0.7 * x[t - 1] + w[t]
#vise o tome na: http://econ.nsysu.edu.tw/ezfiles/124/1124/img/Chapter17_MaximumLikelihoodEstimation.pdf
CSS_function<-function(x)
{
n<-length(x)
x1<-sum(x[1:(n-1)])/(n-1)
x2<-sum(x[2:n])/(n-1)
A<-0
#brojilac
B<-0
#imenilac
for(i in 2:n) A<-A+(x[i-1]-x1)*(x[i]-x2)
for(i in 2:n) B<-B+(x[i-1]-x1)*(x[i-1]-x1)
phi<-A/B
mi<-(x2-phi*x1)/(1-phi)
S<-0
for(i in 2:n)
S<-S+((x[i]-mi)-phi*(x[(i-1)]-mi))^2
sigma_2 <- S/(n-1)
return (c(phi,mi,sigma_2))
}
CSS_function(x)## [1] 0.6086996 0.1824780 0.8287996arima(x, order=c(1,0,0),method="CSS")##
## Call:
## arima(x = x, order = c(1, 0, 0), method = "CSS")
##
## Coefficients:
## ar1 intercept
## 0.6087 0.1825
## s.e. 0.0774 0.2329
##
## sigma^2 estimated as 0.8288: part log likelihood = -132.51#nema metod CSS u ar funkciji, ali ima u arima
#(kod coefficients treba da bude mean a ne intercept)
#nema metod CSS u ar funkciji, ali ima u arima
#jos jedan nacin za ocenjivanje reda AR(1) modela
#preko linearne regresije, mada neki uslovi za primenu tog modela nisu ispunjeni
lm(x~c(NA,x[1:(length(x)-1)]))##
## Call:
## lm(formula = x ~ c(NA, x[1:(length(x) - 1)]))
##
## Coefficients:
## (Intercept) c(NA, x[1:(length(x) - 1)])
## 0.0714 0.6087ar(x,method="ols")##
## Call:
## ar(x = x, method = "ols")
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.5845 -0.1363 0.2827 -0.1708 0.1218 -0.3677 0.1373 -0.1020
## 9 10 11 12
## 0.2990 -0.0948 0.0741 -0.2642
##
## Intercept: 0.07347 (0.08502)
##
## Order selected 12 sigma^2 estimated as 0.6194MLE - ZA AR(1)
Pažnja, kod korišćenja arima.sim funkcije: dva oblika zapisa AR(1) modela (može i za više redove) \(x_t= \mu + \phi*(x_{t-1} - \mu) + e_t\) \(x_t=\alpha + \phi*x_{t-1} + e_t\), gde je \(\alpha=\mu*(1-phi)\) \(\mu\) je srednja vrednost serije \(x_t\) - mean \(\alpha\) - zovemo intercept (kao pomeraj, kod jednačine prave) ipak funkcija arima.sim parametar mean posmatra kao da se zadaje intercept
set.seed(101)
data <- arima.sim(n=1000,list(ar=0.1),mean=10)
#to mozete da proverite sa:
mean(data) ## [1] 11.07179#ova funkcija u stvari vraca negtivni loglikelihood (funkciju verodostojnosti)
#zbog log=TRUE se logaritmuje
logl <- function(sigma,alpha,phi) {
-sum(dnorm(data[-1],alpha+phi*data[1:length(data)-1],sigma,log=TRUE))-dnorm(data[1],alpha/(1-phi),sigma/sqrt(1-phi^2),log=TRUE)
}
library(stats4)
ocene<-mle(logl,start=list(sigma=1,alpha=10,phi=0.05),method="L-BFGS-B",lower=c(0, -Inf, -1), upper = c(Inf, Inf, 1))
#logl je u stvari negativni log likelihood, pa trazimo minimum
#trazi minimum funkcije nekim numerickim metodama, iterativno, potrebno je
#postaviti pocetne vrednosti
#ovako trazimo minimum funkcije 3 promenljive
#ako hocemo da prvo ocenimo srednju vrednost pa onda za pomerenu seriju da ocenjujemo
#preostale parametre sigma i phi
data1<-data-mean(data)
mu=mean(data)
logl1<-function(sigma,phi)
{
-sum(dnorm(data1[-1], phi*data1[1:length(data1)-1],sigma,log=TRUE))-dnorm(data1[1],0,sigma/sqrt(1-phi^2), log=TRUE)
}
ocene1<-mle(logl1, start=list(sigma=1,phi=0.5),method="L-BFGS-B",lower=c(0,-1),upper=c(Inf,1))
#ovaj drugi nacin ocenjivanja je kako radi funkcija ar, ako ima arguemnt demean=TRUE - prvo oceni parametar mu sa srednjom vrednoscu
#a onda ocenjuje ostale
#ocene1 bi trebalo da se poklapaju sa x.ar1
x.ar1<-ar(data, order=1, method = "mle",demean=TRUE)
x.ar1$x.mean*(1-x.ar1$ar) #ovo e mean u stvari## [1] 10.02769#i kod funkcije ar je mean u stvari intercept
#smatra se da je srednja vrednost 0
x.ar=ar(data, order=1, method = "mle",demean=FALSE)
x.ar$x.mean #ne ocenjuje se srednja vrednost## [1] NA#koeficijenti dobijeni su cudni - phi=0.9922 jako je veliko, skoro 1
#sto je kao Braunovo kretanje
#funkciju ar sa metodom mle mozemo koristiti i kad ne znamo unapred broj parametara
#u tom slucaju trazi se model sa najmanjom vrednosti AIC (Akaike Information Criterion)
x.ar$aic## 0 1
## 4289.381 0.000#ispisuje se niz brojeva: na k-tom mestu je razlika u aic izmedju modela sa k-1 parametara
#i najboljeg modela, dakle x.ar$aic[x.ar$order+1] je UVEK 0
#funkcijom ar se ne moze dobiti vrednost aic, to moze funkcijom arima
#arima je druga funkcija, uopstenija od ar koja se koristi za nalazenje najboljeg modela
#ona ne oduzima prvo srednju vrednost nego i nju ocenjuje
#po definiciji AIC=-2*loglikelihood+2*#parameters
arima(data,order=c(1,0,0),method="ML") #mozete proveriti formulu za AIC i odavde##
## Call:
## arima(x = data, order = c(1, 0, 0), method = "ML")
##
## Coefficients:
## ar1 intercept
## 0.0943 11.0719
## s.e. 0.0315 0.0336
##
## sigma^2 estimated as 0.9237: log likelihood = -1379.26, aic = 2764.52#ako se ne navede method u funkciji arima
#podrazumevano je da se prvo primeni CSS za nalazenje pocetnih vrednosti parametara
#pa onda ML metoda
#vrednost AIC se moze naci i preko rucno odredjenih ocena mle
#poklapaju se!
2*logl1(0.99948791,0.04660395)+2*3## [1] 2769.634#raspodela ocena AR modela dobijenih metodama YW, MLE, CSS je asimptotski normalna
#srednja vrednost -stvarna vrednost parametra, disperzija...
#asy.var - asymptotic variance-covariance matrix
#na osnovu ovog moze se napraviti 95%-ni interval poverenja za parametre
x.ar$order## [1] 1x.ar$ar## [1] 0.9932656x.ar$asy.var## [,1]
## [1,] 0.001805989#95%-ni interval poverenja
x.ar$ar + c(-2, 2) * sqrt(x.ar$asy.var)## [1] 0.9082717 1.0782595PRIMER1
www <- "http://www.math.rs/p/files/1550577126-70-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.823251Z.ar$order #nismo naveli red modela, ar funkcija je odredila najbolji na osnovu AIC## [1] 1Z.ar$ar## [1] 0.890261Z.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])
#posto funkcija ar prvo oduzima srednju vrednost pa onda uklapa u model
#tj onda ocenjuje parametre
#vrednosti predvidjene modelom dobijaju se tako sto se doda ta srednja vrednost
#\hat{z}_t= 2.82+ 0.89*(z_{t-1}-2.82)
mean(Z.ts)## [1] 2.823251plot(Z.ts)
#na grafik dodajemo predvidjene vrednosti
lines( Z.ar$resid+Z.ts, col="red") 
#PRIMER2
www<-"http://www.math.rs/p/files/1550577150-70-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))
mean(aggregate(Global.ts, FUN = mean))## [1] -0.1382628Global.ar$order## [1] 4Global.ar$ar## [1] 0.68245556 0.00319301 0.06716095 0.17296979acf(Global.ar$res[-(1:Global.ar$order)], lag = 50)
#ovde smo uklapali Global u AR model, za koji se smatra da nema deterministicku komponetu #iako smo ranije pomenuli da izgleda kao da postoji rastuci trend #ova serija moze se uklopiti u razlicite modele, i nijedan nije pogresan #ali ukoliko imamo neko objasnjenje, dokaz da postoji trend #onda ne bismo uklapali u ovaj model