Несезонски 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(p,d,q)\) модели

\(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,1,1)\) без константе је експоненцијално равнање.

Из једнаине \(ARIMA\) модела видимо да када је \(d\ge 1\) поред ауторегресионог полинома имамо полином \((1-B)^d\), а нула овог полинома је у \(B=1\). Ако је ово случај кажемо да серија има јединични корен, па се ове серије са стохастичким трендом зову и серије са јединичним кореном.

Симулација и одабир \(ARIMA(p,d,q)\) модела

\(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\) \(АR\) ред \(p\) претходног модела.

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

  • Ако ни након поправљања модела серија резидуала не изгледа стационарно можда је потребно повећати ред диференцирања, ког смо у првом кораку, због опасности од непотребног диференцирања фиксирали на \(0\) или \(1\).

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

Предвиђање зависи доста од вредности \(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\), онда је \(\mu\) баш средња вредност серије \(x_t\). У зависности од тога да ли је укључена константа \(c\) у \(ARIMA\) модел и који је ред диференцирања \(d\) имамо следеће случајеве

  • \(c=0, d=0\) - дугорочна предвидања теже ка 0
  • \(c=0, d=1\) - дугорочна предвидања теже ка не-нула константи
  • \(c=0, d=2\) - дугорочна предидања прате праву
  • \(c>0, d=0\) - дугорочна предвидања теже средњој вредности
  • \(c>0, d=1\) - дугорочна предвидања прате праву
  • \(c>0, d=2\) - дугорочна предвидања прате квадратни тренд

Одабрана вредност \(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.

Симулирајмо \(ARIMA(1,1,1)\) модел \(x_t=0.5x_{t-1}+x_{t-1} -0.5x_{t-2}+w_t+0.3w_{t-1}\).

set.seed(1)
x <- w <- rnorm(1000)
for (i in 3:1000) x[i] <- 0.5*x[i-1] + x[i-1] - 0.5*x[i-2] + w[i] + 0.3*w[i-1]
arima(x, order = c(1,1,1))
## 
## Call:
## arima(x = x, order = c(1, 1, 1))
## 
## Coefficients:
##          ar1     ma1
##       0.4235  0.3308
## s.e.  0.0433  0.0450
## 
## sigma^2 estimated as 1.067:  log likelihood = -1450.13,  aic = 2906.26

Пример 1

library(fpp)
## Warning: package 'fpp' was built under R version 4.1.3
## Loading required package: forecast
## Warning: package 'forecast' was built under R version 4.1.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: fma
## Warning: package 'fma' was built under R version 4.1.3
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 4.1.3
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.1.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.1.3
## 
## 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 4.1.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 = 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

Пример 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)
## 
##  KPSS Test for Level Stationarity
## 
## data:  oil
## KPSS Level = 0.64075, Truncation lag parameter = 3, p-value = 0.01893
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 = 2427:  log likelihood = -239.23
## AIC=480.45   AICc=480.55   BIC=482.26
plot(forecast(fit1,h=6))

Пример 3

library(astsa)
## 
## Attaching package: 'astsa'
## The following object is masked from 'package:fpp':
## 
##     oil
## The following objects are masked from 'package:fma':
## 
##     chicken, sales
## The following object is masked from 'package:forecast':
## 
##     gas
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] 52
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 = -2.775, Lag order = 8, p-value = 0.2503
## 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 = 6.3174, Truncation lag parameter = 6, 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

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(1,1,1) 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##           ar1     ma1
##       -0.5253  0.7142
## s.e.   0.0872  0.0683
## 
## sigma^2 = 0.002112:  log likelihood = 904.58
## AIC=-1803.15   AICc=-1803.11   BIC=-1790.25
# predvidjanje
plot(forecast(fit2,h=52))

Задатак

  • Временску серију austa(пакет: fpp2) уклопити у \(ARIMA\) модел, а затим извршити предвиђање за наредних 10 година.