\(\color{deepskyblue}{\text{Несезонски ARIMA модели }}\)

Формалну дефиницију слабе стационарности је битно знати, али у анализи временских серија битно је знати препознати је (колико је то могуће) са графика. Из дефиниције слабе стационарности следи да серија мора имати константну средњу вредност и дисперзију, и то се може видети са графика, док се то да корелациона функција зависи само од дужине временске задршке и не може видети. Предвиђање стационарне временске серије је једноставно, јер се претпоставља да ће се обрасци из прошлости наставити и у будућности. Зато нам је битно да трансформишемо серију, диференцирањем, тако да од нестационарне дођемо до стационарне серије. Ако радимо са нестационарном серијом, стандардне статистике, узорачка средња вредност, дисперзија, немају смисла. На пример, ако средњу вредност нестационарне серије са растућим трендом оцењујемо узорачком средњом вредности, оцене ће бити мање него што би требало да буду.

Од почетка увођења AR модела бавили смо се стационарним временским серијама. Ако важи услов каузалности AR процеси су сигурно стационарни. Видели смо на примеру AR(1) процеса да ако је \(\phi = 1\) ради се о дискретном Брауновом кретању и да то није стационаран процес. Може се доказати итерацијама уназад да ако је \(|\phi| < 1\), такав AR(1) процес јесте стационаран, јер се може записати као линеарни процес. Такође, на сличан начин AR(1) процес код кога је \(|\phi| > 1\) се може записати као бесконачна линеарна комбинација белих шумова из будућности, која јесте стационарни процес, али неупотребљив за предвидање.

Погледајмо још једном дискретно Брауново кретање \[x_t=x_{t-1}+w_t.\] Свака вредност је од претходне удаљена за случајан “корак”. Може се написати и као \[x_t-x_{t-1}=w_t.\] Разлике узастопних вредности чине нову временску серију, означимо је са \(\bigtriangledown x_t=x_t-x_{t-1}\). Ову временску серију зовемо диференцирана временска серија, њом се мери промена између узастопних опсервација. Видимо да је диференцирана серија случајног лутања стационарни процес, јер је баш бели шум. Дакле, почетна серија није била стационарна али диференцирана јесте. Још један пример нестационарности је када серија има детерминистичку компоненту и стационарну случајну компоненту. На пример то је серија облика \[x_t=\mu_t+y_t,\] где је \(\mu_t= \beta_0+\beta_1t\), а \(y_t\) је стационарни процес. Диференцирањем добијамо \[\bigtriangledown x_t=\beta+y_t-y_{t-1}=\beta+\bigtriangledown y_t,\] а ово јесте стационарни процес (проверите). У једначини \(\mu_t\) може да буде и стохастичко, на пример случајно лутање \[\mu_t=\mu_{t-1}-v_t,\] где је \(v_t\) стационарно. У том случају диференцирани процес је \[\bigtriangledown x_t=v_t+\bigtriangledown y_t,\] и он је такође стационаран. Ако је \(\mu_t\) детерминистичко, на пример полином \(k\)-тог степена тада се може доказати да је \(k\) пута диференцирана серија \(\bigtriangledown ^kx_t\) стационарна. Kао и овај детерминистички модел и стохастички модели могу да воде до диференцирања вишег реда. Овде можемо направити разлику између серија које су нестационарне, јер имају детерминистички тренд (trend stationary) и оних које су нестационарне због стохастичког тренда (difference stationary). Многе економске серије имају тренд, битно је раздвојити та два случаја. У првом случају шокови, изненадне промене понашања серије, имају пролазан краткорочан утицај, јер те промене потичу од стационарног процеса који има очекивану вредност нула. У другом случају шокови имају дугорочан утицај чији се ефекат не смањује са временом као што је о код стационарних процеса. Детерминистички тренд наравно не мора да буде ограничен на линеарне функције времена, али бавићемо се само линеарним моделима у временским серијама. Сада ћемо се бавити нестационарним серијама са стохастичким трендом. ARIMA процесе добијамо када укомбинујемо AR, MA и диференцирање.

Процес \(x_t\) је ARIMA(p,d,q) са средњом вредности 0 ако је \[\bigtriangledown^d x_t- (1-B)^d x_t \] ARMA(p,q) процес. Једначину модела можемо записати као \[\phi(B)(1-B)^dx_t=\theta(B)w_t.\] Ако је средња вредност диференцираног процеса \(E(\bigtriangledown^d x_t)=\mu\), онда је једначина модела \[\phi(B)(1-B)^dx_t=c+\theta(B)w_t,\] где је \(c=\mu(1-\phi_1-\dots -\phi_p)\).

Неке већ поменуте моделе сада можемо идентификовати: ARIMA(0,0,0) - бели шум, ARIMA(0,1,0) - случајно лутање, са или без дрифта зависи да ли има константе, ARIMA(п,0,q) - ARMA(п,q), ARIMA(0,1,1) без константе је експоненцијално равнање.

Из једнаине ARIMA модела видимо да када је \(d\ge 1\) поред ауторегресионог полинома имамо полином \((1-B)^d\), а нула овог полинома је у \(B+1\). Ако је ово случај кажемо да серија има јединични корен, па се ове серије са стохастичким трендом зову и серије са јединичним кореном. ARIMA модели је најсложенија класа модела. У односу на до сад посматране моделе, почевши од AR модела, ARIMA су најсложенији и уопштење претходно посматраних. Kао и код претходних модела поставља се питање оцењивање реда модела, оцењивање параметара, како мерити колико је изабрани модел добар и на крају предвидање. На ова питања није било лако одговорити ни само посматрајући AR и MA моделе нижих редова. Описаћемо уобичајени поступак налажења најбољег ARIMA модела за добијену временску серију, притом се ослањајући на функције доступне у R-у.

Да би одредили о ком ARIMA моделу је ради потребно је одредити ред ARIMA модела, тј. бројеве p, d, q. Први корак и најбитнији је одредивање реда диференцирања. Ако са графика временске серије \(x_t\) видимо очигледну нестационарнот, на пример неконстантан тренд, који може бити детерминистички или случајан, можемо пробати да диференцирањем дођемо до стационарне серије. Ако и график диференциране серије \(\bigtriangledown x_t\) не изгледа стационарно можемо диференцирати још једанпут, ово може бити случај ако је почетна серија споро променљива, глатка. Kолико је диференцирања потребно можемо видети и на графику ACF. Ако серија има јединични корен, график ACF ће јако споро опадати. Циљ је диференцирањем доћи до серије чији ће ACF брзо опадати и која има константну средњу вредности. Ипак, најбоље је првo узети да је \(d=0\) или \(d=1\), пошто као што знамо и ACF AR модела може споро опадати, на пример за неке веће вредности параметра AR(1) модела. Веома битно је пазити да се не диференцира превише пута, јер тако можемо увести корелације у диференцирану серију тамо где их испрва није ни било. На пример \(x_t=w_t\), диференцирањем добијамо \(\bigtriangledown x_t =w_t-w_{t-1}\), и ово је МА(1) модел.

Диференцирањем убијамо позитивне аутокорелације, а тако се може јавити и негативна аутокорелација реда 1, додатним диференцирањем та прва аутокорелација још више би се смањила. Ако добијете да је ACF(1)=-0.5 или мање, превише пута сте диференцирали. Оптимални број диференцирања је онај за који деференцирана серија има најмању стандардну девијацију. У вези са диференцирањем је константа \(c\) у једначини модела, да ли ћемо је укључити или не. Ако нисмо диференцирали, константом дозвољавамо да серија има не-нула средњу вредност. Ако смо једном диференцирали значи да у почетној серији имамо тренд не-нула средње вредности. Ако смо два пута диференцирали ствар се компликује и у том случају искључујемо константу \(c\). У случају диференцирања мањег реда, сами одлучујемо о укључивању константе.

Нека је одабран ред диференцирања d. Следећи корак је посматрати узорачке ACF и PACF серије \(\bigtriangledown^d x_t\). Ако нема значајних корелација ради се о константном моделу (само средња вредност) или о случајном лутању у зависности да ли смо диференцирали или не. Ако имамо значајне корелације само реда 1 и на ACF и на PACF узима се да је модел AR(1) ако је она позитивна, односно MA(1) ако је она негативна. Ако ACF опадa, а PACF нестаје после задршке p ради се о AR(p) моделу, а ако PACF опада, а ACF нестаје после задршке q ради се о МА(q) моделу. Такође, пошто радимо са оценама, некада на основу ових графика неће бити јасно ни да ли су корелације нула, или само опадају на графицима ACF или PACF. Често, што не мора увек бити случај, знак да се ради о МА моделу су негативне вредности ACF, а знак AR модела су позитивне вредности ACF.

Након што одаберемо ред модела, следи дијагностика. Ово подразумева посматрање резидуала и упоредивање различитих модела. Резидуали су дефинисани као грешке предвидања један корак унапред \[e_t=x_t-\hat{x}_{t|t-1},\] где је \(\hat{x}_{t|t-1}\) предвидање један корак унапред. Предвидамо вредност серије у тренутку t, ако знамо претходних \(t-1\) вредности \(x_1, \dots x_{t-1}\). Осим резидуала могу се посматрати и стандардизовани резидуали које добијамо када резидуале поделимо са оцењеном стандардном девијацијом предвидања један корак унапред. Рачунањем ових вредности се нећемо бавити, оне се добијају позивима функција у R-у. Ако је модел добар стандардизовани резидуали би требало да буду бели шум са дисперзијом 1. Да ли резидуали ово задовољавају може се видети са самог графика резидуала. За испитивање корелисаности можемо користити узорачки ACF. Осим ако су резидуали нормално расподељени, није довољно проверити да ли су некорелисани (на пример само испитивањем ACF). На пример могуће је да ако процес није Гаусов и није корелисан да ипак различите вредности у времену буду јако зависне. За испитивање нормалности можемо да користимо Q-Q plot, испитујући тако како се квантили резидуала разликују од квантила нормалне расподеле, могуће је користити и хистограм. Веома је могуће да су резидуали некорелисани, али ипак зависни. Тим проблемом бавићемо се касније, уводењем ARCH и GARCH модела. За сада покрићемо само случај када корелограм не изгледа добро.

За почетак довољно је одабрати неки једноставнији модел који је можда и довољно добар, и прећи на следећи корак, дијагностику тог модела. Међутим, ако се први једноставни модел не покаже као добар морамо га надоградити. Тај модел није покупио све корелације у подацима, па је потребн додати неке параметре. Питање је како поправити модел, повећавањем реда AR или МА. - Ако постоји шиљак (spike) на мањим временским задршкама (lag=1,2,3) графика ACF и/или аутокорелација је негативна, повећати за 1 МА ред q претходног модела. - Ако постоји шиљак на мањим временским задршкама (lag=1,2,3) графика PACF и/или парцијална аутокорелација је позитивна, повећати за 1 АР ред p претходног модела.

Уколико се ради о серији која нема сезонску компоненту не морамо обраћати пажњу о задршкама већим од 3. Овај поступак не треба пуно пута понављати, јер углавном би се требало држати модела за које је \(p+q\le3\). Проблем са компликованијим моделима је преприлагоденост (overfitting) и могућа појава истог фактора у AR и MA полиномима. Kао што је поменуто, негативне аутокорелације су често последица диференцирања. Диференцирањем убијамо позитивне аутокорелације које онда могу прећи у негативне. Зато је модел ARIMA(0,1,1) у пракси чешће коришћен него ARIMA(1,1,0). У вези са повећавањем реда диференцирања су и следећа два правила.

Предвидање зависи доста од вредности c и d. Запис ARIMA модела еквивалентан оном који је дат на почетку је \[(1-\phi_1B-\dots -\phi_p B^p)(1-B)^d(x_t-\mu t^d/d!) = (1+\theta B+\dots +\theta_qB^q)w_t,\] где је \(c=\mu(1-\phi_1-\dots-\phi_p)\) и \(\mu\) је средња вредност процеса \((1-B )^d\). Укључивањем константе функција предвиђања укључује полином степена d, иначе полином степена d-1. Ако је d=0, онда је му баш средња вредност серије \(y_t\). У зависности од тога да ли је укључена константа c у ARIMA модел и који је ред диференцирања d имамо следеће случајеве

Одабрана вредност d утиче на интервале поверења предвидања. Што је d веће интервали су шири са порастом времена. Ако је d=0 стандардна девијација интервала поверења се не повећава. Подразумевано у функцији Arima је да је \(c=\mu=0\) када је d>0, а ако је d=0, онда се оцењује \(\mu\) средња вредност серије. Она се оцењује методом МВ, дакле оцена неће бити једнака узорачкој средини (осим за p=q=0). Функцију Arima можемо и да натерамо да не укључује константу у модел аргументом include.mean=FALSE, али није могуће натерати је да укључи консанту ако је d>0. С друге стране ту је аргумент include.drift који дозвољава \(\mu\neq0\) када је d=1. За d>1 константа није дозвољена јер би то значило да предвидања прате квадратни тренд, који се не користи иначе за предвидање. Дакле, када је d=1 параметар се зове дрифт. Један аргумент који покрива оба случаја је include.constant који ако је TRUE поставља include.mean=ТRUE ако је d=0, односно include.drift=TRUE ако је d=1.

Приметимо како смо причали о предвидању, а не и о инверзној трансформацији података, што је битно ако су подаци логаритмовани или диференцирани. Разлог за то је што ће то све софтвер да уради за вас! Препоручено је коришћење функције Arima, јер је у њој исправљена већ поменута грешка око термина mean/intercept, а и има могућност да се и када је d=1 укључи константа, што није могуће код функције arima, она има само параметар include.mean који се игнорише ако је d>0. Добијање предвидања и њихово додавање на график је јако једноставно, користи се функција forcast.

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

library(fpp)
## Warning: package 'fpp' was built under R version 3.3.3
## Loading required package: forecast
## Warning: package 'forecast' was built under R version 3.3.3
## Loading required package: fma
## Warning: package 'fma' was built under R version 3.3.3
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 3.3.3
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.3.2
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
## Warning: package 'tseries' was built under R version 3.3.3
plot(usconsumption[,1]) # Promene u procentima, potrosnje u SAD, kvartalna serija
abline(v=1970:2010,lty="dotted")

Са графика закључујемо са серија нема сезонску компоненту, као и да није потребно диференцирати серију, нема очигледног тренда.

usc <- usconsumption[,1]
tsdisplay(usc)

fit1 <- Arima(usc, order=c(3,0,0))
fit2 <- Arima(usc, order=c(0,0,3))
fit1$aic
## [1] 318.1607
fit2$aic # malo je bolji AR(3) model
## [1] 319.4611
#***automatsko odredjivanje reda***
fit <- auto.arima(usconsumption[,1],seasonal=FALSE)
fit
## Series: usconsumption[, 1] 
## ARIMA(0,0,3) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3    mean
##       0.2542  0.2260  0.2695  0.7562
## s.e.  0.0767  0.0779  0.0692  0.0844
## 
## sigma^2 estimated as 0.3953:  log likelihood=-154.73
## AIC=319.46   AICc=319.84   BIC=334.96
#***popravljanje***
tsdisplay(residuals(fit1))

tsdisplay(residuals(fit2)) # i reziduali slicno izgledaju

#***predvidjanje***
# odaberemo ARIMA(3,0,0)
plot(forecast(fit1,h=10)) # prikazani intervali poverenja su 95%-ni i 80%-ni

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

class(oil) # godisnja proizvodnja nafte u Saudijskoj Arabiji
## [1] "ts"
plot(oil)

adf.test(oil) # nije stacionaran
## 
##  Augmented Dickey-Fuller Test
## 
## data:  oil
## Dickey-Fuller = -3.0058, Lag order = 3, p-value = 0.174
## alternative hypothesis: stationary
kpss.test(oil)
## Warning in kpss.test(oil): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  oil
## KPSS Level = 1.1292, Truncation lag parameter = 1, p-value = 0.01
ndiffs(oil) # ocekivano vraca 1
## [1] 1
doil <- diff(oil) # jednom diferenciramo
tsdisplay(doil)

fit1 <- Arima(oil, order=c(0,1,0))
fit1$aic
## [1] 480.4533
res <- residuals(fit1)
tsdisplay(res) # nema znacajnih autokorelacija

auto.arima(oil) # i ovim smo dosli do istog modela
## Series: oil 
## ARIMA(0,1,0) 
## 
## sigma^2 estimated as 2427:  log likelihood=-239.23
## AIC=480.45   AICc=480.55   BIC=482.26
plot(forecast(fit1,h=6))

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

# paket astsa
class(oil) # cene WTI - West Texas Intermediate (tip sirove nafte) cija cena dalje utice i na cene preradjevina
## [1] "ts"
#cena je u dolarima po barelu serija je nedeljna od 2000. do sredine 2010. god
frequency(oil) # 52 jer je nedeljna vs
## [1] 1
plot(oil)
abline(v=2000:2010,lty="dotted") # ne vidimo sezonsku komponentu

# serija ne izgleda stacionarno, i izgleda da disperzija raste sa porastom trenda
# prvo logaritmujemo
logoil <- log(oil)
plot(logoil) # smirili smo disperziju

adf.test(logoil)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  logoil
## Dickey-Fuller = -3.1828, Lag order = 3, p-value = 0.1037
## alternative hypothesis: stationary
kpss.test(logoil) # napomena: ne moze se uvek verovati testovima jedinicnog korena, ovde ipak mozemo
## Warning in kpss.test(logoil): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  logoil
## KPSS Level = 1.12, Truncation lag parameter = 1, p-value = 0.01
ndiffs(oil) # ocigledno je da serija nije stacionarna, pa cemo diferencirati jednom
## [1] 1
dlogoil <- diff(logoil)
tsdisplay(dlogoil)

# zbog velikog broja autokorelacija, nije jasno sta se desava tacno
# acf i pacf izgledaju slicno
# imamo pozitivnu acf(1) pokusacemo dodavanjem jednog ar parametra da ubije
# pozitivnu autokorelaciju
fit1 <- Arima(oil,order=c(1,1,0),lambda=0)
res <- residuals(fit1)
tsdisplay(res) 

# sada cemo dodati jedan MA parametar jer imamo negativnu prvu autokorelaciju
fit2 <- Arima(oil,order=c(1,1,1),lambda=0)
res <- residuals(fit2)
tsdisplay(res)

#ovim smo prvih nekoliko acf i pacf vrednosti smirili, tako da su u granicama
Box.test(res,type="Ljung-Box",lag=104, fitdf=2)
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = NA, df = 102, p-value = NA
# prolazi test! iako ima i znacajnih autokorelacija
# sveukuono ih ima puno, pa je ok da neke ne budu u intervalu
qqnorm(res) # i ne izgleda bas normalno ali nema veze

auto.arima(oil,d=1,lambda = 0,seasonal=FALSE) #dobijamo isti model
## Series: oil 
## ARIMA(0,1,0) 
## Box Cox transformation: lambda= 0 
## 
## sigma^2 estimated as 0.02371:  log likelihood=20.34
## AIC=-38.68   AICc=-38.59   BIC=-36.87
# predvidjanje
plot(forecast(fit2,h=52))

# skoro da je ravna linija, predvidjanje za narednu godinu
# nema puno smisla jer vidimo kako se menja puno i tokom nedelja
# ipak znace nam intervali poverenja, kao donja i gornja granica
# intervala u kom ce se kretati serija

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

  • Симулирати временску серију \(x_t=0.5x_{t-1} -0.5x_{t-2}+w_t+0.3w_{t-1}\) и уклопити је у одговарајући ARIMA модел.