\(\color{deepskyblue}{\text{Временске серије и примене у финансијама}}\)
\(\color{deepskyblue}{\text{Бојана Тодић}}\)
\(\color{deepskyblue}{\text{Сезонски ARIMA модели}}\)
Често се на корелограму зависност од прошлости може видети на задршкама које су умношци сезонског периода. На пример, код месечних временских серија, значајне корелације на задршкама l=12,24,… упућују на постојање сезонске компоненте. Слично, код кварталних серија, на задршкама l=4,8,… На пример, месечна серија коју чине просечне температуре ће имати јаку позитивну корелациију са задршком 12. Ово није нови појам, и већ смо објаснили један метод уклањања сезонске компоненте - сезонског прилагодавања. Увешћемо још један метод. До сада смо се бавили ARIMA моделима којима се моделирају нестационарне серије без сезонске компоненте. Нестационарност потиче од тренда који се диференцирањем убија. Постојање сезонске компоненте такође чини серију нестационарном, пошто њена функција средње вредности није константна. ARIMA модели се могу генерализовати тако да се њима описује и сезонска компонента. Њихова ознака је \[ARIMA(p,d,q)(P,D,Q)_s.\] Kао и до сада, мала слова означавају редом, ред AR полинома, број (обичних) диференцирања до стационарности, и ред МА полинома. Великим словима означени су одговарајући сезонски параметри. У полиномијалном запису то је \[\Theta_P(B^s)\theta_p(B)(1-B^s)^Dx_t=c+\Phi_Q(B^s)\phi_q(B)w_t,\] где су \(\theta_p\) и \(\phi_q\) редом АR и МА полиноми које смо већ дефинисали, а \(\Theta_P\) и \(\Phi_Q\) су сезонски АR и сезонски МА полином, \(s\) је сезонски период. \[\Phi_P(B^s)=1-\Phi_1B^s-\Phi_2B^{2s}-...-\Phi_PB^{Ps}\] \[\Theta_Q(B^s)1+\Theta_1B^s+\Theta_2B^{2s}+...+\Theta_QB^{Qs}\] Уопште ови модели су нестационарни, осим када је \(d=D=0\) и корени полинома са леве стране једначине ARIMA модела су по апсолутној вредности већи од 1. Инвертибилост важи ако су корени полинома са десне стране по апсолутној вредности већи од 1. Навешћемо неке једноставне примере сезонских ARIMA модела
а) \(ARIMA(0,0,0)(1,0,0)_{12}\) \[x_t=\Phi x_{t-12} +w_t\] Овај модел је одговарајући за месечну временску серију код које само вредност у одговарајућем месецу прошле године утиче на садашњост. Овај модел је стациоаран ако је \(|\Phi^{-1/12}|>1\).
б) \(ARIMA(0,1,0)(1,0,0)_{12}\) \[x_t=x_{t-1}+\Phi x_{t-12}-\alpha x_{t-13}+w_t\] Ред модела идентификујемо записом преко полинома. \[(1-\Phi B^{12})(1-B)x_t=w-t \] Овај модел није стационаран јер постоји јединични корен.
в) \(ARIMA(0,0,0)(1,0,1)_{12}\) \[x_t=\Phi x_{13} +w_t+\theta w_{t-12}\] Ред модела идентификујемо записом преко полинома: \[(1-\Phi B^{12})x_t=(1+\Theta B^{12}) w_{t-12}. \] Овај модел је каузалан ако је \(|\Phi|<1\) и инвертибилан ако \(|\Theta|<1\). У случају када су модели чисто сезонски (p, d, q = 0) рачунање аутокорелационе функције (као и оцењивање параметара и предвидање о коме нисмо детаљно причали код ARIMA модела) је аналогно несезонском случају.
г) \(ARIMA(0, 1, 1)(0, 1, 1)_{12}\) \[x_t=x_{t-1}+x_{t-12}-x_{t-13}+w_t+\theta w_{t-1}+\Theta w_{t-12} +\Theta \theta w_{t-13}\] Због мултипликативне природе модела, коефицијент уз \(w_{t-13}\) није самосталан параметар већ је производ коефицијената уз \(w_{t-1}\) и \(w_{t-12}\).
Сезонски ARIMA модели су још комплекснији, питање је како уопште одредити ред модела. Неки једноставни сезонски ARIMA модели (са много нула) ипак имају карактеристичан облик ACF, PACF функција. На пример, за \(ARIMA(0, 0, 0)(0, 0, 1)_{12}\) важи
Велика позитивна вредност ACF(12), остале вредности ACF блиске 0.
Експоненцијално опадање сезонских задршки на графику PACF(12,24,36…).
Слично, за \(ARIMA(0, 0, 0)(1, 0, 0)_{12}\) важи
Велика позитивна вредност PACF(12), остале вредности PACF блиске 0.
Експоненцијално опадање сезонских задршки на графику ACF (12,24,36…).
Процес уклапања у модел је сличан као и код несезонских ARIMA модела. Прво је потребно одредити број потребних диференцирања, \(d\) - несезонских и \(D\) - сезонских.
У пракси \(d+D\) не би требало да буде веће од 2, а ако је \(d+D=2\) параметар \(c\) (intercept) не оцењујемо. Дакле, ако је \(d=0, D=1\), онда је серија \(x_t-x_{t-s}\) стационарна, а ако је \(d=1,D=1\), онда је серија \((x_t − x_{t-1})−(x_{t-s} − x_t − s − 1)\) стационарна. Значење сезонских редова је, ако је \(P=1\), онда се за предвидање у неком тренутку \(t\), \(x_t\), користи вредност \(x_t-s\), а ако је \(Q = 1\), онда се за предвидање \(x_t\) користи \(w_{t-s}\). Сада ћемо навести кораке у идентификацији сезонских модела.
Ако се чини да је сезонска компонента мултипликативна логаритмовати почетну серију.
Покушати разне комбинације диференцирања \(d=0,1, \; D=0,1\) да би се дошло до стационарности. Уколико се одлучимо за \(d=1,D=1\) није битно којим редоследом смо дошли до тога (да ли смо прво обично диференцирали, па нисмо били задовољни, па онда и сезонски или обрнуто). Ипак, уколико на серији видимо јаку сезонску компоненту најбоље је прво сезонски диференцирати, јер се некад цела нестационарност уклања само тим диференцирањем, и нема потребе за обичним, док ако прво обично диференцирамо, тиме нећемо уклонити сезонски ефекат. Тестовима јединичног корена које смо поменули (ADF (adf.test), KPSS (kpss.test)) тестира се постојање “обичног” јединичног корена, тј полинома \((1-B)\). Постоје и тестови сезонског јединичног корена о којима нећемо причати. У пакету fpp налази се функција nsdis која користећи тестове сезонског јединичног корена утврдује број потребних сезонских диференцирања.
Након диференцирања посматрамо ACF и PACF на умношцима сезонског периода. Шпицеви (велике позитивне вредности) у ACF на задршкама s,2s,3s,… и позитиван шпиц у PACF(s) сугерише на један сезонски AR параметар (SAR=1). Негативни шпиц у ACF на задршци s и негативни шпицеви у PACF на задршкама s,2s,3s,… сугерише на један сезонски МА параметар (SМА=1). Често SМА=1 иде у комбинацији са сезонским диференцирањем Дакле, принцип је сличан као у несезонском случају, само што се сада фокусирамо на понашање функција ACF, PACF у умношцима сезонског периода.
Често се добије да је \(d=1, D=1\). Због оба диференцирања често ћемо имати негативне аутокорелације и на lag=1 и lag=s. То сугерише додавање МА параметара, и један уобичајени сезонски модел \(ARIMA(0, 1, q)(0, 1, 1)_s\).
Често, ако је \(d=0, D=1\), имамо позитивне аутокорелације, тада пробати додавањем AR чланова (p=1,2,3) ради њиховог убијања. Тада се може наћи негативна аутокорелација за lag=s. У том случају додати један сезонски МА параметар SМА=1. Добијени модел је још један уобичајени сезонски модел \(ARIMA(p, 0, 0)(0, 1, 1)_s + c\).
Пример ACF и PACF функција неког сезонског ARIMA модела.
phi=c(rep(0,11),.8)
ACF=ARMAacf(ar=phi, ma=-0.5, 50)[-1]
PACF=ARMAacf(ar=phi, ma=-0.5, 50, pacf=TRUE)
plot(ACF, type="h")
abline(h=0)
plot(PACF, type="h")
abline(h=0)
\(\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.3class(euretail)## [1] "ts"frequency(euretail)## [1] 4plot(euretail)
abline(v=1996:2011,lty="dotted") #Vidi se i sezonska komponenta
Диференцирање временске серије
nsdiffs(euretail) #prvo sezonski diferenciramo## [1] 1tsdisplay(diff(euretail,4))
ndiffs(diff(euretail,4)) # ni ovo ne izgleda stacionarno, pa cemo i obicno diferencirati## [1] 1tsdisplay(diff(diff(euretail,4))) #ovo izgleda malo stacionarnije, srednja vrednost prblizno nula
На ACF и PACF се види значајност на сезонским задршкама, што нам више одговара. Одабрали смо \(d=1,D=1\). Остаје да одаберемо и остале параметре модела.
Имамо негативну вредност ACF(4), али позитивну ACF(1). Покушаћемо са \(АRIMA(0,1,1)(0,1,1)_4\) моделом.
Функција Arima, слично као arima, има још неке могућности које нећемо користити.(arima је из stats пакета, а Arima из forecast).
fit <- Arima(euretail, order=c(0,1,1), seasonal=c(0,1,1)) #probati i 1,1,0 0,1,1
tsdisplay(residuals(fit))
По резидуалима видимо да смо одабраним SMA=1 решили параметре, али значајне корелације на мањим задршкама сугеришу да треба додати још неки несезонски параметар. Питање је који параметар додати. Код избора модела несезонских arima модела речено је да ако је негативни шиљак на ACF за мале lag, онда се додаје један MA параметар, а ако је позитивни шиљак на PACF за мале lag, онда се додаје један AR параметар. Према томе, овде бисмо могли да додамо и један AR несезонски или MA несезонски параметар. Овде нема јединственог правила. Ми ћемо додати MA параметар.
fit4 <- Arima(euretail, order=c(1,1,1), seasonal=c(0,1,1))
res <- residuals(fit4)
tsdisplay(res)
fit4$aic## [1] 68.37079fit2 <- Arima(euretail, order=c(0,1,2), seasonal=c(0,1,1))
res <- residuals(fit2)
tsdisplay(res)
fit2$aic## [1] 73.61446Box.test(res, lag=16, fitdf=3, type="Ljung") #lag=16, znaci posmatramo ukupno 16 autokor, #fitdf=3 jer smo ocenili 3 parmetara 2xAR, 1XSMA, #nulta hipoteza je da su autokorelacije 0##
## Box-Ljung test
##
## data: res
## X-squared = 18.665, df = 13, p-value = 0.1339qqnorm(res)
Гледајући p вредност Љунг-Бокс теста, рекли бисмо да је модел добар, али ипак видимо да су још значајне ACF 3 и PACF 3, па ћемо покушати додати још један МА параметар.
fit3 <- Arima(euretail, order=c(0,1,3), seasonal=c(0,1,1))
res <- residuals(fit3)
tsdisplay(res)
Box.test(res, lag=16, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: res
## X-squared = 7.0105, df = 12, p-value = 0.8569fit3$aic #Sada i reziduali bolje izgledaju, i veca je p vrednost testa, a mozemo da uporedimo AIC vrednosti, ili AICc## [1] 67.39644fit2$aicc## [1] 74.3552fit3$aic #ovo je bolje jer je manje## [1] 67.39644Извршимо предвиђање на основу изабраног модела. Наставићемо да радимо са пакетом forecast, у ком се налази функција за предвиђање forecast. Иначе, постоји и функција predict за предвиђање.
plot(forecast(fit3,h=12)) #predvidjanje za naredne 3 godina
Предвиђање прате скорашњи опадајући тренд, али интервали поверења су велики и дозвољавају у неком тренутку да постане и растући.
Функција аутоматски одређује модел (може одредити и сезонске компоненте)
auto.arima(euretail)## Series: euretail
## ARIMA(1,1,2)(0,1,1)[4]
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.7345 -0.4655 0.2162 -0.8413
## s.e. 0.2239 0.1995 0.2096 0.1869
##
## sigma^2 estimated as 0.1592: log likelihood=-29.69
## AIC=69.37 AICc=70.51 BIC=79.76Након мало чекања, добије се да је најбољи модел \(ARIMA(1,1,2)(0,1,1)\), а његова вредност \(AICc\) је чак већа од нашег одабраног модела. То се дешава јер ова функција не проверава баш све моделе, него прави неке пречице у налажењу модела, па некад и не да стварно најбољи модел. Те пречице се могу искључити, па следећим позивом добијамо најбољи модел који смо и сами идентификовали.
auto.arima(euretail,stepwise=FALSE, approximation=FALSE)## Series: euretail
## ARIMA(0,1,3)(0,1,1)[4]
##
## Coefficients:
## ma1 ma2 ma3 sma1
## 0.2625 0.3697 0.4194 -0.6615
## s.e. 0.1239 0.1260 0.1296 0.1555
##
## sigma^2 estimated as 0.1564: log likelihood=-28.7
## AIC=67.4 AICc=68.53 BIC=77.78\(\color{deepskyblue}{\text{Пример 2}}\)
plot(h02) #Mesecna serija, prodaja lekova sa kortikosteriodima u Australiji
Сезонска компонента је јасно уочљива и изгледа као да расте са порастом тренда. Серија је позитивна, па је можемо логаритмовати.
lh02 <- log(h02)
par(mfrow=c(2,1))
plot(h02, ylab="H02 sales (million scripts)", xlab="Year")
plot(lh02, ylab="Log H02 sales", xlab="Year")
Са графика се може видети како смо логаритамском трансформацијом мало смирили дисперзију, али идаље имамо благо растући тренд и очигледну сезонску компоненту.
Диференцирање временске серије
tsdisplay(diff(lh02,12),
main="Seasonally differenced H02 scripts", xlab="Year") #prvo sezonsko
nsdiffs(lh02) ## [1] 1#da li je potrebno i obicno diferencirati?
kpss.test(diff(lh02,12)) #kpss test kaze da, a kod adf je p vrednost 0.06## Warning in kpss.test(diff(lh02, 12)): p-value smaller than printed p-value##
## KPSS Test for Level Stationarity
##
## data: diff(lh02, 12)
## KPSS Level = 1.012, Truncation lag parameter = 3, p-value = 0.01#to takodje kaze i funkcija ndiffs
ndiffs(diff(lh02,12)) ## [1] 1Погледајмо график ACF и PACF када и обично диференцирамо.
tsdisplay(diff(diff(lh02,12),1)) # ako su ACF vrednosti manje od -0.5 onda je previse puta diferencirano, ipak necemo jos jednom diferencirati
tsdisplay(diff(lh02,12)) 
Ово одговара другом описаном уобичајеном сезонском моделу ARIMA(p,0,0)(0,1,1)+c. Имамо доста позитивних аутокорелација које убијамо додавањем AR параметара.
fit1<-Arima(h02, order=c(1,0,0),seasonal=c(0,1,1),include.mean=TRUE,lambda=0)
res<-residuals(fit1)
tsdisplay(res) #od mnogo pozitivnih, dodavanjem AR parametara, dobili jako male ACF i PACF 1
Мораћемо да поправимо ову ситуацију, али не додавањем AR параметара, јер ћемо тако направити само још негативније ACF (покушајте). Дакле, одустајемо од уобичајеног модела.
Покушајмо додавањем MA несезонског параметра.
fit2<-Arima(h02, order=c(1,0,1),seasonal=c(0,1,1),lambda=0)
res<-residuals(fit2)
tsdisplay(res)
fit3<-Arima(h02, order=c(1,0,2),seasonal=c(0,1,1),lambda=0)
res<-residuals(fit3)
tsdisplay(res) 
Никако се не поправља ситуација, на великим lag значајне аутокорелације, па ћемо променити и сезонске компоненте. Крећемо испочетка. Пробаћемо са додавањем још једног SMA параметра. Ставићемо један AR параметар због позитивних аутокорелација на малим lag.
tsdisplay(diff(lh02,12))
fit4<-Arima(h02, order=c(1,0,0),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit4)
tsdisplay(res)
fit4$aicc## [1] -398.8206Додајемо још један MA због негативних аутокорелација на малим lag.
fit5<-Arima(h02, order=c(1,0,1),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit5)
tsdisplay(res)
fit5$aicc## [1] -474.8237Додајемо још два AR чиме се поправља ситуација на малим lag. Иако остају неке значајне ACF и PACF за велике lag, али то није на сезонским задршкама, па их игноришемо.
fit6<-Arima(h02, order=c(3,0,1),seasonal=c(0,1,2),lambda=0)
res<-residuals(fit6)
tsdisplay(res)
fit6$aicc## [1] -485.4754#uzimamo za lag (broj prvih autokorelacije koje se zajedno posmatraju) da je neki
#umnozak sezonskog perioda, fitdf je broj ocenjenih parametara, a to je 7 zbog mean
Box.test(res, lag=36, fitdf=7, type="Ljung")##
## Box-Ljung test
##
## data: res
## X-squared = 50.712, df = 29, p-value = 0.007561Резидуали нашег најбољег модела нису прошли тест. Некада није могуће наћи модел који ће проћи све тестове. Покушајмо са функцијом auto.arima.
#nametnucemo nase diferenciranje i transformaciju
# fit<-auto.arima(h02,lambda=0,d=0,D=1,max.order=9,stepwise=FALSE, approximation = FALSE)
#max.order je najveca vrednost za p+q+P+QДобијамо да је најбољи ARIMA(4,0,3)(0,1,1) модел и он пролази Љунг-Боксов тест (проверити).
Одабрали смо ARIMA(3,0,1)(0,1,2) модел и предвиђамо вредности.
plot(forecast(fit6)) #po defaultu predvidjanje je 2*sezonski period, dakle 2 godine
\(\color{deepskyblue}{\text{Задатак }}\)
- Уклопити временску серију cbe у одговарајући сезонски АRIMA модел.