\(\color{deepskyblue}{\text{Сврха регресије}}\)

Трендови у временским серијама могу бити стохастички или детерминистички. Можемо сматрати да је тренд стохастички ако на њему видимо промене у кретању које се не могу објаснити неком погодном функцијом времена, а уочљиве пролазне трендове приписујемо јакој серијској корелацији са случајном грешком. Трендови овог типа су чести у серијама везаним за финансије и могу бити симулирани у R-у користећи моделе као што су случајно лутање или ауторегресивни процес. С друге стране, када имамо неко разумно физичко објашњење тренда, то желимо и да искористимо и објаснимо тренд детерминистички. На пример, детерминистички растући тренд може бити последица пораста популације, а неке промене које се циклично понављају могу настати због одређене сезонске компоненте. Детерминистички трендови и сезонске варијације могу се моделовати помоћу регресије.

Регресија код временских серија се најчешће разликује од стандардне регресионе анализе, јер су резидуали који такође чине једну временску серију, најчешће корелисани. Kада је серијска корелација позитивна, оцене стандардних грешака оцењених параметара, које у R-у добијамо примењујући стандардни поступак за регресиону анализу, теже да буду мање од њихових стварних вредности. Због тога се статистичким тестовима додељује већа значајност него што би требало (p-вредности тестова су мање него што треба да буду) у стандардном излазу који се добија у R-у.

\(\color{deepskyblue}{\text{Линеарни модели}}\)

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

Напомене:

  1. При наведеним претпоставкама ови процеси су стационарни

  2. За линеарне процесе може се извести коваријациона функција као функција од коефицијената \(\psi_j\)

ВОЛДОВА ТЕОРЕМА О ДЕKОМПОЗИЦИЈИ (Wold decompositon):

Сваки слабо стационарни процес са средњом вредношћу 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.

Недетерминистички стационарни процеси су они који немају детерминистичку компоненту и то су онда линеарни процеси.

Kажемо да је модел временске серије \(\{x_t : t = 1, ... , n\}\) линеаран ако се може написати у облику: \[x_t = \alpha_0 + \alpha_1 u_{1,t} + \alpha_2 u_{2,t} + ... + \alpha_m u_{m,t} + z_t,\] где је \(u_{i,t}\) вредност \(i\)-те променљиве објашњења у тренутку \(t\), \((i=1,\dots, m; t=1,\dots,n)\), \(\{z_t\}\) грешка у тренутку \(t\), а \(\alpha_0, \alpha_1,\dots, \alpha_m\) су параметри модела, који могу бити оцењени методом најмањих квадрата. Грешке формирају временску серију \(\{z_t\}\), са функцијом средње вредности која је идентички једнака нули, која не мора да буде Гаусова, нити бели шум.

Пример линеарног модела је полином \(p\)-тог степена у функцији од \(t\):

\[x_t = \alpha_0 + \alpha_1 t + \alpha_2 t^2 + ... + \alpha_p t^p + z_t \]

Променљиве објашњења у овом случају су \(u_{i,t} = t^i (i = 1, ... , p)\). Израз “линеаран” овде се односи на линеарност у односу на параметре, а не предикторе.

Једноставан пример линеарног модела је права, чија се једначина добија заменом \(p=1\) у претходној једначини:

\[x_t = \alpha_0 + \alpha_1 t + z_t.\]

У овом случају, вредност тренда у тренутку \(t\), је тачка са праве, тј. \(m_t = \alpha_0 + \alpha_1 t\). За полиноме вишег реда, вредност тренда у тренутку \(t\) је вредност полинома из једначине у тачки \(t\), тј. \(m_t = \alpha_0 + \alpha_1 t + \alpha_2 t^2 + ... + \alpha_p t^p\)

Многи нелинеарни модели могу бити трансформисани у линеарне. На пример, модел \(x_t = e^{\alpha_0 + \alpha_1 t + z_t}\) временске серије \(\{x_t\}\) може бити трансформисан у линеарни логаритмовањем.

\[y_t = \log{x_t} = \alpha_0 + \alpha_1 t + z_t \quad (3)\]

\(\color{deepskyblue}{\text{Стационарност}}\)

Линеарни модел временске серије биће нестационаран ако садржи неке функције које зависе од времена. Диференцирање често може трансформисати нестационарну серију са детерминистичким трендом у стационарну. На пример, ако је \(\{x_t\}\) временска серија које се моделује линијским трендом на који је додат бели шум \(x_t = \alpha_0 + \alpha_1 t + z_t\),𝑥𝑡 онда је једном диференцирана серија задата са:

\[∇x_t = x_t - x_{t-1} = z_t - z_{t-1} + \alpha_1 \quad (4)\]

Ако важи да је серија грешака \(\{z_t\}\) стационарна, онда је и серија \(\{∇x_t \}\), јер је функција средње вредности константна и корелациона функција зависи само од дужине временског интервала. Диференцирање може да трансформише нестационарну серију са стохастичким трендом у стационарну (пример је случајно лутање), што значи да диференцирање може уклонити и стохастички и детерминистички тренд из серије. Ако је тренд полином степена \(m\)𝑚, почетну серију треба диференцирати \(m\) 𝑚пута да би се уклонио тренд.

\(\color{deepskyblue}{\text{Симулације}}\)

Уобичајено је да временскa серијa грешака \(\{z_t\}\) буде аутокорелисана. Наводимо код којим се симулира и графички представља временска серија са растућим линеарним трендом \[m_t = 50 + 3t\] и аутокорелисаним резидуалима.

set.seed(1)
z <- w <- rnorm(100, sd = 20)
for (t in 2 : 100) {
  z[t] <- 0.8 * z[t-1] + w[t]
}
t <- 1 : 100
x <- 50 + 3 * t + z
plot(x, xlab = "vreme", type = "l")

Модел који одговара наведеном коду се може изразити као \[x_t = 50 + 3t + z_t,\] где је \(\{z_t\}\) један AR(1) процес задат са \(z_t = 0.8z_{t-1} + w_t\), при чему је \(\{w_t\}\) Гаусовски бели шум са стандардним одступањем \(\sigma = 20\). График добијене временске серије (у зависности од времена) приказан је на претходној слици.

\(\color{deepskyblue}{\text{Фитовани модели}}\)

Линеарни модели се најчешће формирају минимизирањем суме квадрата резидуала

\[\sum_{i=1}^{n} z_t^2 = \sum_{i=1}^{n} (x_t - \alpha_0 - \alpha_1 u_{1,t} - \alpha_2 u_{2,t} - ... - \alpha_m u_{m,t})^2, \]

што се у R-у постиже применом функције lm, на следећи начин:

x.lm <- lm(x~t)
coef(x.lm)
## (Intercept)           t 
##   58.551218    3.063275
vcov(x.lm)
##             (Intercept)            t
## (Intercept)   23.815013 -0.355447951
## t             -0.355448  0.007038573
sqrt(diag(vcov(x.lm)))
## (Intercept)           t 
##  4.88006278  0.08389621

Оцене параметара линеарног модела се могу издвојити функцијом coef. Приметимо да су, као што је и очекивано, оцене блиске вредностима параметара које смо задали приликом симулације: оцењена вредност за intercept (пресек са \(x\)𝑥-осом) је 58.55, што је релативно блиско броју 50, а нагиб праве (коефицијент уз \(t\)𝑡) је 3.06 ≈ 3. Стандардне грешке се могу добити као квадратни корени елемената на главној дијагонали матрице коју чине оцењене коваријације између оцењених параметара линеарног модела. Уграђена функција у R-у која даје ту матрицу је vcov, а као аргумент јој се прослеђује посматрани линеарни или нелинеарни модел. Међутим, овако добијене вредности стандардних грешака су често потцењене због аутокорелације резидуала.

Након фитовања одговарајућег модела, можемо размотрити разне дијагностичке графике. У случају регресије временских серија, један од веома битних дијагностичких графика је корелограм резидуала.

Аутокорелација (серијска корелација) је корелација између вредности временске серије у различитим временским тренуцима. Парцијална аутокорелација са кашњењем k је корелација између \(x_t\) и \(x_{t-k}\) након што је елиминисан утицај свих елемената између. График аутокорелационе функције назива се корелограм, а график парцијалне аутокорелационе функције – парцијални корелограм.

acf(resid(x.lm))

pacf(resid(x.lm))

Kао што је очекивано, са првог графика видимо да је серија резидуала аутокорелисана, зато што више од 5% вредности излази ван сегмента ограниченог плавим испрекиданим линијама. Други график указује на то да је само парцијална аутокорелација са кашњењем 1 значајна, што имплицира да серија резидуала представља AR(1) процес. Овакав резултат је очекиван, с обзиром да смо приликом симулације резидуала користили управо AR(1) процес. Узгред, приметимо да acf рачуна корелације за кашњења почевши од нуле, док pacf рачуна почевши од јединице.

\(\color{deepskyblue}{\text{Пример - индекси потрошачких цена}}\)

Индекси потрошачких цена дефинишу се као мера просечне промене малопродајних цена робе и услуга које се користе за личну потрошњу типичног потрошача. Kористи се за праћење промена трошкова живота у времену.

baza <- read.table("CPIAUCSL.txt", header = F)
serija <- ts(baza$V2, start = c(1947,1), end = c(2015,3), frequency = 12)
plot(serija, xlab = "vreme", ylab = "CPI", main = "Consumer Price Index for All Urban Consumers")



Са графика се јасно види да постоји растући детерминистички тренд, који је у периоду од 1947. до око 1975. године растао споријим темпом у односу на наредне године у којима је раст био бржи, све до данашњег дана. Такође, можемо уочити скок и одмах затим пад нивоа у 2009. години, што је вероватно била последица светске економске кризе.

Издвојићемо податке од 1980. до 2007. године како бисмо им придружили одговарајући линеарни модел, с обзиром да су индекси потрошачких цена у том периоду пратили приближно линеарни тренд. Време је једини предиктор модела, па ћемо функцијом time издвојити временске тренутке из ts објекта.

cpi <- window(serija, start = c(1980,1), end = c(2007,12))
cpi.lm <- lm(cpi ~ time(cpi))
plot(cpi, xlab = "vreme", ylab = "CPI", main = "Consumer Price Index (1980 - 2007)")
abline(cpi.lm, lty = 2, lwd = 2, col = "goldenrod1")

summary(cpi.lm)
## 
## Call:
## lm(formula = cpi ~ time(cpi))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4168 -1.9114  0.2692  1.4515  5.4432 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.611e+03  2.617e+01  -329.0   <2e-16 ***
## time(cpi)    4.391e+00  1.313e-02   334.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.945 on 334 degrees of freedom
## Multiple R-squared:  0.997,  Adjusted R-squared:  0.997 
## F-statistic: 1.119e+05 on 1 and 334 DF,  p-value: < 2.2e-16

Како је Adjusted R-squared: 0.997, можемо закључити да модел готово савршено одговара подацима.

confint(cpi.lm)
##                    2.5 %       97.5 %
## (Intercept) -8662.458216 -8559.487895
## time(cpi)       4.365286     4.416926
acf(resid(cpi.lm))

Функцијом confint одређујемо 95%-тни интервал поверења за параметре модела. Интервал поверења за нагиб праве не садржи нулу, што представља статистичко објашњење за постојање тренда, у случају да је аутокорелација резидуала занемарљива. Међутим, серија резидуала је позитивно аутокорелисана, што доводи до тога су вредности стандардних грешака оцењених параметара мање него што би требало да буду, а самим тим су и интервали поверења ужи. Уколико у временској серији постоји изражени тренд, онда аутокорелациона функција опада од јединице скоро линеарно, као што је случај у овом примеру.

\(\color{deepskyblue}{\text{Уопштени метод најмањих квадрата}}\)

У регресији временских серија уобичајено је и очекивано да серија резидуала буде аутокорелисана. Kао што је већ речено, у случају позитивне серијске корелације у серији резидуала, стандардне грешке оцењених параметара имају тенденцију да буду потцењене, па би стога требало извршити корекцију.

Уопштени метод најмањих квадрата (generalised least square (GLS)) се може користити за добијање бољих оцена стандардних грешака параметара регресије, тако да те оцене објашњавају аутокорелацију резидуала. Овај метод се базира на максимизирању веродостојности уколико је дата аутокорелација података. Одговарајућа функција у R-у се назива gls и налази се у пакету nlme.

\(\color{deepskyblue}{\text{Фитовање модела симулираној серији помоћу функције GLS}}\)

Наредни пример илуструје употребу функције gls којим се фитује модел временској серији са којом смо се већ сусрели.

library(nlme)
x.gls <- gls(x ~ t, cor = corAR1(0.8))
summary(x.gls)
## Generalized least squares fit by REML
##   Model: x ~ t 
##   Data: NULL 
##        AIC      BIC    logLik
##   862.8866 873.2265 -427.4433
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  Parameter estimate(s):
##       Phi 
## 0.7161368 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 58.23302 11.924568  4.883449       0
## t            3.04225  0.202445 15.027538       0
## 
##  Correlation: 
##   (Intr)
## t -0.857
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6171329 -0.6195428  0.0353972  0.5836326  2.3184155 
## 
## Residual standard error: 25.58595 
## Degrees of freedom: 100 total; 98 residual
sqrt(diag(vcov(x.gls)))
## (Intercept)           t 
##  11.9245679   0.2024447

Аргументом \(cor = corAR1(0.8)\) прецизира се да је аутокорелација са кашњењем 1 једнака 0.8 (зато што је та вредност коришћена приликом раније симулације података). У случају реалних података, тј. неке временске серије која је настала на основу података из прошлости, аутокорелација са кашњењем 1 се процењује на основу корелограма резидуала фитованог линеарног модела. То значи да прво мора да се направи линеарни модел обичном методом најмањих квадрата и да се, затим, из корелограма резидуала тог модела очита вредност аутокорелације са кашњењем 1.

У наведеном примеру, стандардне грешке параметара модела су 11.925 за пресек са \(x\)𝑥-осом и 0.202 за нагиб праве, тј. знатно су веће од стандардних грешака добијених применом обичне методе најмањих квадрата (у поглављу 3), које су износиле 4.88 за пресек са \(x\)𝑥-осом и 0.084 за нагиб. Оцене стандардниx грешака које су добијене помоду gls су тачније, зато што узимају у обзир серијску корелацију резидуала.

Могу се прочитати из излаза који даје функција суммарy или се могу “ручно” одредити наредбом sqrt(diag(vcov(x.gls))).

Такође, можемо приметити да се оцене самих параметара два посматрана модела мало разликују, због пондерисања. На пример, нагиб праве износи 3.06 према lm моделу, односно 3.04 према gls моделу.

\(\color{deepskyblue}{\text{Интервали поверења за индексе потрошачких цена}}\)

Да бисмо одредили 95%-тни интервал поверења за тренд у серији са индексима потрошачких цена, тако да серијска корелација у серији резидуала буде узета у обзир, користимо функцију gls. Kао аргумент функције gls наводимо да серија резидуала отприлике прати АР(1) процес са аутокорелацијом са кашњењем 1 која износи 0.959, што се може очитати са корелограма.

acf(resid(cpi.lm))[1]

## 
## Autocorrelations of series 'resid(cpi.lm)', by lag
## 
##     1 
## 0.959
cpi.gls <- gls(cpi ~ time(cpi), cor = corAR1(0.959))
confint(cpi.gls)
##                     2.5 %       97.5 %
## (Intercept) -10352.056812 -8421.215908
## time(cpi)        4.318123     5.242116

Иако су интервали поверења сада шири него приликом претходне, нула и даље није садржана у интервалима, што значи да су оцене параметара статистички значајне. Стога, постоји статистички доказ растућег тренда у периоду од 1980. до 2007. године, и очекује се да се такав тренд настави у наредним годинама.

\(\color{deepskyblue}{\text{Линеарни модели са сезонским компонентама}}\)

Како временске серије представљају опсервације мерене у току времена, сезонске компоненте се често јављају у подацима. Размотримо линеарне регресионе моделе са предикторима који описују ефекат сезонских промена.

\(\color{deepskyblue}{\text{Адитивне сезонске индикаторске случајне величине}}\)

Нека временска серија садржи \(s\) сезона. Сезонски индикаторски модел за серију \(\{x_t : t = 1, ... , n\}\) која има \(s\) сезона и тренд \(m_t\) задат је са: \[x_t = m_t + s_t + z_t\] где је \(s_t = \beta_i\) уколико се тренутак \(t\) јавња у \(i\)-тој сезони \((t = 1,\dots,n;i = 1,\dots, s\), \(\{z_t\}\) серија резидуала која може бити аутокорелисана. Овај модел има исти облик као адитивна декомпозиција модела, али се разликује у параметрима који дефинишу тренд. У претходној једнакости \(m_t\) нема слободан члан (intercept), нпр. \(m_t\) би могао да буде полином степена \(p\) са коефицијентима \(\alpha_1,\dots,\alpha_p\). Дакле, ова једнакост је еквивалентна полиномијалном тренду у коме слободан члан зависи од сезоне, тако да сезонским параметрима \(\beta_1, ... , \beta_s\) одговара \(s\) могућих слободних чланова једнакости (2). Претходна једнакост се другачије може записати: \[x_t = m_t + \beta_{1 + ((t - 1) \mod s)} + z_t\]

На пример, за временску серију \(\{x_t\}\) чије су вредности бележене сваког месеца почевши са \(t\) = 1 у јануару, сезонски индикаторски модели са линеарним трендом задат је на следећи начин:

\[x_t = \alpha_1t + s_t + z_t =\left\{\begin{matrix} \alpha_1t + \beta_1 + z_t &t=1, 13, ... \\ \alpha_1t + \beta_2 + z_t &t=2, 14, ... \\ ... \\ \alpha_1t + \beta_{12} + z_t &t=12, 24, ... \end{matrix}\right. \]

Параметри претходног модела могу се оценити обичном или уопштеном методом најмањих квадрата, третирајући сезонску компоненту као фактор. У R-у се функција factor може применити на сезонске индексе који из временске серије могу издвојити функцијом cycle.

\(\color{deepskyblue}{\text{Пример - незапосленост жена}}\)

Илустроваћемо описани модел на примеру временске серије која бележи број незапослених жена у Сједињеним Америчким Државама у периоду од јануара 1948. до априла 2015. године.

women <- read.table("LNU03000002.txt", header = T)
wmn <- ts(women$X562, start = c(1948, 1), end = c(2015, 4), frequency = 12)
plot(wmn, xlab = "vreme", ylab = "hiljada osoba", main = "Nezaposlenost žena")

На \(y\)𝑦-оси је назначено да ова временска серија прати колико хиљада жена је незапослено у одређеном временском тренутку. Kао што можемо видети са графика, гледано на дуже стазе, делује као да постоји известан растући тренд, што је вероватно последица пораста броја становника у САД-у у последњих шездесетак година (с обзиром да се посматра број, а не проценат, незапослених). Такође, уочава се и нагли скок у броју незаполених у 2009. и 2010. години, због светске економске кризе. Уколико се нацрта график за краћи временски интервал, постојање сезонске компоненте је сасвим уочљиво. У оквиру једне године, највиши је ниво незаполености у летњим месецима (јул и август), а најмањи у периоду пре Нове године (новембар и децембар) и у априлу.

sezone <- cycle(wmn)
vreme <- time(wmn)
wmn.lm <- lm(wmn ~ 0 + vreme + factor(sezone), data = wmn)

\(\color{deepskyblue}{\text{Прогнозирање помоћу регресије}}\)

Битна ставка у регресионој анализи временских серија је предвиђање будућности на основу података из прошлости. Притом се користи екстраполација одабраног модела, тј. налазе се вредности у будућим тренуцима. Главни проблем оваквог приступа је у томе што се тренд може променити у будућности. Због тога је боље гледати на предвиђање, које се базира на регресионом моделу, као на очекивану вредност под условом да се тренд из прошлости наставља и у будућности.

Сада ћемо одредити двадесетомесечну прогнозу будућег броја незпослених жена, помоћу већ направљеног модела. У наредном коду, прво креирамо временске тренутке који одговарају месецима почев од јануара 2015. до децембра 2016. а затим одстрањујемо прва 4 месеца, зато што већ имамо реализоване вредности закључно са априлом 2015. и тиме добијамо вектор new.t са будућим временским тренуцима.

Прво наводимо “ручно израчунату” прогнозу: нагиб праве која представља тренд је \(\alpha\), док \(\beta\) чине сезонски параметри, па су тада прогнозиране вредности \((\alpha * new.t + \beta)\).

new.t <- seq(2015, len = 2 * 12, by = 1/12)[5 : 24]
new.t
##  [1] 2015.333 2015.417 2015.500 2015.583 2015.667 2015.750 2015.833
##  [8] 2015.917 2016.000 2016.083 2016.167 2016.250 2016.333 2016.417
## [15] 2016.500 2016.583 2016.667 2016.750 2016.833 2016.917
alpha <- coef(wmn.lm)[1]
beta <- rep(coef(wmn.lm)[2 : 13], 2)[5 : 24]
(alpha * new.t + beta)
##  factor(sezone)5  factor(sezone)6  factor(sezone)7  factor(sezone)8 
##         5219.573         5212.543         5106.857         5005.290 
##  factor(sezone)9 factor(sezone)10 factor(sezone)11 factor(sezone)12 
##         4885.021         4859.633         4663.305         5021.170 
##  factor(sezone)1  factor(sezone)2  factor(sezone)3  factor(sezone)4 
##         4976.434         4890.140         4764.346         4842.949 
##  factor(sezone)5  factor(sezone)6  factor(sezone)7  factor(sezone)8 
##         5281.634         5274.604         5168.918         5067.351 
##  factor(sezone)9 factor(sezone)10 factor(sezone)11 factor(sezone)12 
##         4947.082         4921.694         4725.365         5083.231

Алтернативно, може се користити генеричка функција за предвиђање у R-у, која се назива predict. Kао параметри јој се прослеђују фитовани модел и нови подаци (нови временски тренуци, у случају анализе временских серија) у којима желимо да одредимо прогнозиране вредности. Посебно треба обратити пажњу на то да нови подаци буду адекватно дефинисани и именовани у оквиру објекта класе data.frame.

Неопходно је да колоне буду именоване исто као предиктори модела (време и сезоне, у овом примеру). У супротном, R ће јавити грешку.

new.dat <- data.frame(vreme = new.t, sezone = c(5 : 12, 1 : 12))
predict(wmn.lm, new.dat)
##        1        2        3        4        5        6        7        8 
## 5219.573 5212.543 5106.857 5005.290 4885.021 4859.633 4663.305 5021.170 
##        9       10       11       12       13       14       15       16 
## 4976.434 4890.140 4764.346 4842.949 5281.634 5274.604 5168.918 5067.351 
##       17       18       19       20 
## 4947.082 4921.694 4725.365 5083.231

На оба начина добијају се исте вредности.

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

  • Уклопити временску серију global у линеарни модел. Одредити 95% интервале поверења за оцењене параметре модела. Нацртати корелограм резидуала добијеног модела.
  • Уклопити временску серију global (од 1986-2005) у линеарни модел са трендом који је линеарна функција времена и са адитивном сезонском индикаторском случајном променљивом. На основу добијеног модела извршити предвиђање за наредне две године.