\(\color{deepskyblue}{\text{Временске серије и примене у финансијама}}\)
\(\color{deepskyblue}{\text{Бојана Тодић}}\)
\(\color{deepskyblue}{\text{Holt-Winters метод}}\)
Holt-Winters метод је метод који се користи када временска серија има и тренд и сезонску компоненту који се мењају током времена. Овај модел може да се користи и за предвиђање.
\(\color{deepskyblue}{\text{Једначине адитивног Holt-Winters метода}}\)
\[a_t=\alpha(x_t-s_{t-p})+(1-\alpha)(a_{t-1}+b_{t-1})\] \[b_t=\beta(a_t-a_{t-1})+(1-\beta)b_{t-1}\] \[s_t=\gamma(x_t-a_t)+(1-\gamma)s_{t-p}\] где су \(a_t\), \(b_t\) и \(s_t\) оцењени ниво (тренд), нагиб и сезонски ефекат у тренутку \(t\), а \(\alpha\), \(\beta\) и \(\gamma\) су параметри изравнања.
До сада смо имали само појам тренда и сезонске компоненте. У овом моделу посматрамо тренд и његов нагиб ради бољег описивања временске серије. Слободно речено, нагиб представља извод тренда тј. брзину промене функције тренда, а ниво саму вредност тренда у тренутку \(t\).
У наведеним рекурентним једначинама по којима се одређују \(a_t\), \(b_t\) и \(s_t\) користи се линеарна комбинација две различите могуће оцене за те компоненте на основу вредности претходно оцењених компоненти.
У једначини за ниво у тренутку \(t\) прва могућа оцена је разлика вредности опсервације у тренутку \(t\) и одговарајуће сезонске компоненте оцењене за претходни период (ово је добра оцена за ниво ако се сезонска компонента не мења пуно по годинама). Друга могућа оцена за ниво добијена је на основу претпоставке да се тренд оцењен претходног тренутка линеарно наставља брзином која је такође оцењена за претходни тренутак. То је управо вредност \(a_{t−1} + b_{t−1}\).
У једначини за нагиб у тренутку \(t\) прва могућа оцена је разлика нивоа оцењених у тренутку \(t\) и \(t−1\). Ово је логична оцена као извод, али у дискретном времену, мера промене тренда. Друга оцена је само оцена нагиба у претходном тренутку, и то је добра оцена уколико је тренд линеаран, јер је тада извод (нагиб) константан. Уколико се мења непредвидиво, за оцену нагиба је боља прва оцена.
У једначини за сезонски ефекат у тренутку \(t\) прва оцена је разлика опсервације у тренутку \(t\) и нивоа у том тренутку (који је оцењен у првој једначини). Друга оцена је само сезонска компонента из истог месеца претходне године, што смо већ имали оцењено. Прва оцена је боља ако се сезонска компонента значајно променила у односу на ону претходну, а друга оцена је боља ако она остаје приближно иста.
Параметре равнања \(\alpha\), \(\beta\) и \(\gamma\) можемо експлицитно задати у функцији HoltWinters() или препустити R-у да их сам оцени.
Прво је потребно прецизирати почетне вредности за \(a_t\), \(b_t\) и \(s_t\). Најчешће се узима \(a_1 = x_1\), а вредности \(b_1\), \(s_1\), …, \(s_p\) се постављају на 0 или се оцењују функцијом decompose, што је по default-у у овој функцији у R-у. Осим ове функције постоји и функција hw у пакету forecast где се почетне вредности бирају мало другачије.
Сада кад смо видели смисао и знамо да израчунамо ова три низа, питање је како се овај модел користи за прогнозирање. За адитивни Holt-Winters модел, једначина изравнате серије је заправо њена оцењена вредност један корак унапред \[\bar{x_t} = a_{t−1} + b_{t−1} + s_{t−p}\] Вратимо се на оцене параметара равнања. Уколико се они не наведу, рачунају се тако да грешка између оцењене и изравнате серије буде минимизирана. Дакле, тражи се минимум по \((\alpha,\beta,\gamma)\) функције \[\sum_{t=1}^{n}(x_t-\bar{x_t})^2.\]
Вредност те минимизиране суме може се добити позивом model.hw$SSE, где је SSE скраћено од “Sum of Squared Errors”, а са model.hw обележен је објекат који враћа функција HoltWinters.
wine <- read.delim("wine.txt")
attach (wine)
sweetw.ts <- ts(sweetw, start = c(1980,1), freq = 12)
plot(sweetw.ts, xlab= "Time (months)", ylab = "sales (1000 litres)") 
Временска серија “wine.txt” садржи податке о продаји различитих врста аустралијског вина (multiple ts) у периоду од 1.1980. до 7.1995. године. Са графика се уочава да продаја слатког вина бележи велики раст након 1985. године.
sweetw.hw <- HoltWinters(sweetw.ts)
sweetw.hw # ispisuju se parametri modela i koeficijenti (vrednosti koje sluze za predvidjanje)## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = sweetw.ts)
##
## Smoothing parameters:
## alpha: 0.3633865
## beta : 0
## gamma: 0.4975705
##
## Coefficients:
## [,1]
## a 276.000769
## b 1.350962
## s1 -7.738580
## s2 1.385719
## s3 16.493531
## s4 65.887912
## s5 164.900275
## s6 -66.482340
## s7 -12.045735
## s8 -53.525533
## s9 -36.901263
## s10 -69.418391
## s11 -62.164586
## s12 -8.661942sweetw.hw$SSE # minimizirana greska predvidjanja jedan korak unapred## [1] 548748.9head(sweetw.hw$fitted) # pocinje sa godinom zakasnjenja, moguce je ocenjivat nakon p godina (ovde 1 godina)## xhat level trend season
## [1,] 140.22943 117.9028 1.350962 20.9756944
## [2,] 125.46422 122.8042 1.350962 1.3090278
## [3,] 151.10686 126.5302 1.350962 23.2256944
## [4,] 130.12254 129.2959 1.350962 -0.5243056
## [5,] 91.53074 131.3291 1.350962 -41.1493056
## [6,] 95.89717 132.4872 1.350962 -37.9409722plot(sweetw.hw) 
Важи једнакост \(\bar{x}= level+trend+season\), у шта се можемо уверити и позивом sweetw.hw$fitted и поређењем вредности.
На графику имамо две линије, црна је оригинлна временска серија, а црвена је изравната временска серија предвиђена моделом. Видимо да се црвена почиње мало после црне, а то је због индекса \(p\), јер је могуће оцењивати после прве године, тј. у општем случају за \(t ≥ p\), где је \(p\) период.
plot (sweetw.hw$fitted) 
На графику добијамо приказане све компоненте модела.
sweetw.hw$coef # niz ocenjenih koeficijenata modela u poslednjim trenucima## a b s1 s2 s3 s4
## 276.000769 1.350962 -7.738580 1.385719 16.493531 65.887912
## s5 s6 s7 s8 s9 s10
## 164.900275 -66.482340 -12.045735 -53.525533 -36.901263 -69.418391
## s11 s12
## -62.164586 -8.661942Позивом \(sweetw.hw\$coef\) добија се низ оцењених коефицијената модела у последњим тренуцима, који се такође рачунају по претходној формули. Дакле, ако имамо n опсервација добијени коефицијенти су \(a = a_n\), \(b=b_n\) и \(s_{n−p+1}\),…,\(s_n\) оцене за 12 сезонских компоненти (код годишње серије). Они су посебно издвојени, јер се на основу њих врши предвиђење.
Једначина предвиђања k корака унапред (k корака након завршетка серије) je \[\bar{x}_{n+k|n} = a_n + k·b_n + s_{n+k−p}, \quad k \leq p.\] Као што можемо видети из ове једначине предвиђања, претпоставља се да се тренд након завршетка серије линеарно наставља следећих k корака у будућности почев од последњег нивоа, нагибом оцењеним за последњи тренутак. За сезонску компоненту користи се одговарајућа последња оцењена вредност. Имајући то у виду, малом променом индекса код сезонске компоненте добија се формула предвиђања и у случају \(k > p\).
\(\color{deepskyblue}{\text{Једначине мултипликативног Holt-Winters метода}}\)
\[a_t=\alpha(\frac{x_t}{s_{t-p}})+(1-\alpha)(a_{t-1}+b_{t-1})\] \[b_t=\beta(a_t-a_{t-1})+(1-\beta)b_{t-1}\] \[s_t=\gamma(\frac{x_t}{a_t})+(1-\gamma)s_{t-p}\] где су \(a_t\), \(b_t\) и \(s_t\) оцењени ниво, нагиб и сезонски ефекат у тренутку t, а \(\alpha\), \(\beta\) и \(\gamma\) су параметри изравнања.
Мултипликативни модел користимо када на основу графика серије сумњамо да сезонска компонента расте са порастом тренда. Рекурентне једначине овог модела су аналогне адитивном моделу и њихов смисао се лако може описати као што смо то већ урадили и случају адитивног модела.
Формула по којој се рачуна изравната серија је \[x_t = (a_{t−1}+b_{t−1})·s_{t−p}\] Позивима функција у R-у за HoltWinters модел (само треба навести \(seasonal="mult"\) као аргумент функције) добијамо одговарајуће грешке и коефицијенте за мултипликативни модел. Код мултипликативног модела је \(\bar{x}=(level+trend)season\). Ако је дато n опсервација укупно, једначина предвиђања k корака унапред је \[\bar{x}_{n+k|n} = (a_n + k · b_n) · s_{n+k−p}, \quad k \leq p.\]
Уколико је временска серија добро оцењена Holt-Winters моделом, резидуали (грешке предвижања за један корак унапред) би требало да буду некорелисане (што се проверава корелограмом), или ако није баш очигледно Ljung-Box тестом. Ljung-Box тест за разлику од плавих линија на корелограму (које исто представњају тестирање) не проверава појединачно аутокорелације, већ их посматра заједно. Тест статистика је збир квадрата (скалираних) оцењених аутокорелација што при нултој хипотезио некорелисаности има \(\chi^2\) расподелу. Број степени слободе је број оцењених аутокорелација у суми, што се у R-у задаје параметром lag у функцији Box.test, а критичне су велике вредности тест статистике. Такође, резидуали би требало да буду приближно нормално расподељени са константном дисперзијом. Константност дисперзије у свим временским тренуцима мође се закњучити на основу графика резидуала, а нормалност се може испитати Q-Q графиком.
sweetw.hw <- HoltWinters(sweetw.ts, seasonal="mult")
sweetw.hw## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = sweetw.ts, seasonal = "mult")
##
## Smoothing parameters:
## alpha: 0.4086698
## beta : 0
## gamma: 0.4929402
##
## Coefficients:
## [,1]
## a 285.6890314
## b 1.3509615
## s1 0.9498541
## s2 0.9767623
## s3 1.0275900
## s4 1.1991924
## s5 1.5463100
## s6 0.6730235
## s7 0.8925981
## s8 0.7557814
## s9 0.8227500
## s10 0.7241711
## s11 0.7434861
## s12 0.9472648plot(sweetw.hw)
sweetw.hw$SSE## [1] 477693.9Видимо да је \(\beta=0\), што значи да се вредност нагиба узима последња оцењена вредност. Значи, да се нагиб тренда споро мења за разлику од нивоа и сезонске компоненте. Мултипликативнии модел је бољи од адитивног - ПРОВЕРИТИ! Мултипликативни модел ћемо користити за предвиђање.
sweetw.hw.predict <- predict(sweetw.hw, n.ahead=4*12)
ts.plot(sweetw.ts,sweetw.hw.predict,lty=1:2)
library(forecast) # Drugi nacin za prognoziranje - paket forecast## Warning: package 'forecast' was built under R version 3.3.3
sweetw.hw2 <- forecast(sweetw.hw, h=4*12) # argument h - za koliko vrednosti unapred predvidjamo
plot(sweetw.hw2)
Свеједно је да ли се за предвиђање користи функција predict или forecast, добијају се исте вредности. Само код функције forecast добијамо и интервале поверења.
У овом примеру нисмо се бавили резидуалима, тј. испитивањем да ли је модел добар. Видети у наредним примерима шта се проверава за резидуале, па тестирајте то и на овом примеру.
\(\color{deepskyblue}{\text{Експоненцијално изравнање}}\)
Експоненцијално изравнање је најједноставнији специјални случај Holt-Winters метода. Користимо га када имамо серију без систематског тренда (очигледан дугорочни растући или опадајући) и сезонске компоненте. Средња вредност те серије би требало у том случају да буде приближно константна, јер да је средња вредност на пример растућа, онда би постојао тренд, што је супротно претпоставци. Једначина овог модела је \[x_t = \mu_t + w_t, \quad k=1,2,...\] где је \(\mu_t\) функција средње вредности, а \(w_t\) процес грешке са независним вредностима, очекивањем 0 и дисперзијом \(\sigma^2\). Без знања Holt-Winters метода, ту функцију средње вредности бисмо могли да оценимо методом покретних просека (MA), али то не би било добро, јер претпостављамо да се та функција средње вресности мења непредвидиво, па не би било згодно користити једнаке тежине свих опсервација у околини t за оцену \(\mu_t\). Једина компонента Holt-Winters метода је ниво, пошто нема сезонске компоненте ни тренда који би имао свој нагиб. Ниво се рачуна по формули: \[a_t=\alpha x_t+(1-\alpha)a_{t-1},\] где је \(\alpha \in (0,1)\) параметар равнања. Што је \(\alpha\) мање, мање вредност \(x_t\) утиче на \(a_t\). Велико \(\alpha\) је добро користити када су промене у \(\mu_t\) веће у односу на девијацију \(\sigma\). Тада се почетна серија веома мало изравна (што нам одговара, јер већи део варијабилности потиче од \(\mu_t\)). Мале вредности \(\alpha\) добро користити када се \(\mu_t\) споро мења, а \(\sigma\) је велико. Тада се серија више изравна - вредности у тренутку t су блиске оним у тренутку t−1, а само \(x_t\) мало утиче на \(a_t\). Питање може бити како знамо одакле потиче већи део варијабилности. Ако се бавимо серијом чија је функција средње вреднодти приближно константна, више нам одговарају мање вредности за \(\alpha\), често се узима на пример 0,2. Прогнозирање на основу HW метода је дато са \[\bar{x}_{n+k|n} = a_n.\]
Ово је логична оцена, јер претпоставњамо да се функција средње вредности споро мења (приближно је константна), па је најбоље будућност оценити последњом познатом вредношћу нивоа. Додатно, пошто процес грешке има независне вредности са средњом вредности 0, ту грешку никако не мођемо оценити и укључити у прогнозу. Ради бољег разумевања експоненцијалног изравнања једначину за ниво можемо записати у другачијем облику: \[a_t=\alpha x_t+\alpha(1-\alpha)x_{t-1}+\alpha(1-\alpha)^2x_{t-2}+...\] Одавде се види да вредности у најближим тренуцима највиже утичу на вредности нивоа, а оне у дањим све мање и мање. Тежине експоненцијано опадају, пошто је геометријски низ експоненцијална функција у природним тачкама, одакле и потиче име овог метода. Ово је добро ако се процес мења у непредвидивим тренуцима. За разлику од MA метода овде утичу све вредности из прошлости (код MA метода утиче само претходних и будућих 6 у случају годишње серије).
\(\color{deepskyblue}{\text{Пример 1}}\)
Motor <- read.csv("motororg.txt", sep="")
attach(Motor) # Broj zalbi upucenih AMS-u u periodi od 1996. do 1999. godine
Comp.ts <- ts(complaints, start = c(1996, 1), freq = 12)
plot(Comp.ts, xlab = "Time / months", ylab = "Complaints")
На основу графика видимо да је средња вредност приближно 20, али не видимо јасно ни тренд ни сезонску компоненту. Зато користимо експоненцијално изравнање. Потребно је подесити параметре \(beta\) и \(gamma\) на FALSE, а пошто се не наводи вредност аргумента \(alpha\), она се оцењује тако да сума квадрата резидуала предвиђања за један корак буде минимална. Желимо да предвидимо број жалби у 2000. години и оценимо тренд и сезонску компоненту.
Comp.hw1 <- HoltWinters(complaints, beta = FALSE, gamma = FALSE)
Comp.hw1## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = complaints, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.1429622
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 17.70343plot(Comp.hw1)
Ако сами задамо alpha: За \(\alpha=1\) - иста серија, померена један корак унапред, јер немамо \(a_{1-1} (a_0)\) за оцењивање \(a_{1}\), па изравњена серија код експоненцијалног изравнања креће од \(a_2\). Када имамо и сезонску компонент, на пример годишњу, серија је померена целу годину унапред.
Comp.hw2 <- HoltWinters(complaints, alpha = 0.3, beta=FALSE, gamma=FALSE)
plot(Comp.hw2)
Comp.hw1$fitted # Fitovane vrednosti, tj. vrednosti izravnate serije su bas jednake nivou## Time Series:
## Start = 2
## End = 48
## Frequency = 1
## xhat level
## 2 27.00000 27.00000
## 3 28.00074 28.00074
## 4 28.42952 28.42952
## 5 27.79626 27.79626
## 6 26.39577 26.39577
## 7 25.33845 25.33845
## 8 24.14637 24.14637
## 9 22.40990 22.40990
## 10 22.92315 22.92315
## 11 21.64747 21.64747
## 12 21.12602 21.12602
## 13 22.82355 22.82355
## 14 23.99248 23.99248
## 15 24.99429 24.99429
## 16 24.13733 24.13733
## 17 23.11696 23.11696
## 18 21.24173 21.24173
## 19 21.77902 21.77902
## 20 22.38246 22.38246
## 21 21.75594 21.75594
## 22 21.21898 21.21898
## 23 19.61509 19.61509
## 24 17.38272 17.38272
## 25 17.75689 17.75689
## 26 19.22127 19.22127
## 27 19.47556 19.47556
## 28 19.26461 19.26461
## 29 19.79863 19.79863
## 30 19.68446 19.68446
## 31 19.15772 19.15772
## 32 17.84851 17.84851
## 33 18.72794 18.72794
## 34 17.62314 17.62314
## 35 16.67628 16.67628
## 36 16.15071 16.15071
## 37 17.70175 17.70175
## 38 17.74439 17.74439
## 39 18.06685 18.06685
## 40 18.48618 18.48618
## 41 16.70113 16.70113
## 42 15.60016 15.60016
## 43 17.51583 17.51583
## 44 16.72728 16.72728
## 45 17.05219 17.05219
## 46 16.61584 16.61584
## 47 16.95669 16.95669
## 48 16.81992 16.81992Comp.hw1.predict <- predict(Comp.hw1, n.ahead=12)
Comp.hw1.predict## Time Series:
## Start = 49
## End = 60
## Frequency = 1
## fit
## [1,] 17.70343
## [2,] 17.70343
## [3,] 17.70343
## [4,] 17.70343
## [5,] 17.70343
## [6,] 17.70343
## [7,] 17.70343
## [8,] 17.70343
## [9,] 17.70343
## [10,] 17.70343
## [11,] 17.70343
## [12,] 17.70343Добијамо да је 17.7 Оцена на крају 1999. године (за било који корак унапред).
\(\color{deepskyblue}{\text{Пример 2}}\)
rain <- read.delim("rain.txt")
attach(rain) # U seriji su date godisnje vrednosti kolicine kise (u incima) u Londonu od 1813
rainseries <- ts(Rainfall,start=c(1813))
plot(rainseries)
Са графика не видимо јасан тренд ни неки циклус, а ниво је приближно константан.
rainseriesforecasts <- HoltWinters(rainseries, beta=FALSE, gamma=FALSE)
rainseriesforecasts # Dobija se jako mala vrednost alpha## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = rainseries, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.02412151
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 24.67819rainseriesforecasts$fitted## Time Series:
## Start = 1814
## End = 1912
## Frequency = 1
## xhat level
## 1814 23.56000 23.56000
## 1815 23.62054 23.62054
## 1816 23.57808 23.57808
## 1817 23.76290 23.76290
## 1818 23.76017 23.76017
## 1819 23.76306 23.76306
## 1820 23.82691 23.82691
## 1821 23.79900 23.79900
## 1822 23.98935 23.98935
## 1823 23.98623 23.98623
## 1824 23.98921 23.98921
## 1825 24.19282 24.19282
## 1826 24.17032 24.17032
## 1827 24.13171 24.13171
## 1828 24.10442 24.10442
## 1829 24.19549 24.19549
## 1830 24.22261 24.22261
## 1831 24.24329 24.24329
## 1832 24.32812 24.32812
## 1833 24.21938 24.21938
## 1834 24.23290 24.23290
## 1835 24.13369 24.13369
## 1836 24.13867 24.13867
## 1837 24.21782 24.21782
## 1838 24.10257 24.10257
## 1839 24.04293 24.04293
## 1840 24.12608 24.12608
## 1841 24.01280 24.01280
## 1842 24.18448 24.18448
## 1843 24.15808 24.15808
## 1844 24.19889 24.19889
## 1845 24.16153 24.16153
## 1846 24.12748 24.12748
## 1847 24.18133 24.18133
## 1848 24.02499 24.02499
## 1849 24.16454 24.16454
## 1850 24.13476 24.13476
## 1851 24.01621 24.01621
## 1852 23.93453 23.93453
## 1853 24.20964 24.20964
## 1854 24.25018 24.25018
## 1855 24.11509 24.11509
## 1856 24.08964 24.08964
## 1857 24.04430 24.04430
## 1858 23.99933 23.99933
## 1859 23.87319 23.87319
## 1860 23.97780 23.97780
## 1861 24.17710 24.17710
## 1862 24.13110 24.13110
## 1863 24.21405 24.21405
## 1864 24.15075 24.15075
## 1865 23.97658 23.97658
## 1866 24.10933 24.10933
## 1867 24.29001 24.29001
## 1868 24.33729 24.33729
## 1869 24.31468 24.31468
## 1870 24.34134 24.34134
## 1871 24.26847 24.26847
## 1872 24.28659 24.28659
## 1873 24.51752 24.51752
## 1874 24.47295 24.47295
## 1875 24.33660 24.33660
## 1876 24.43558 24.43558
## 1877 24.47717 24.47717
## 1878 24.56625 24.56625
## 1879 24.79573 24.79573
## 1880 25.01341 25.01341
## 1881 25.14045 25.14045
## 1882 25.20750 25.20750
## 1883 25.25411 25.25411
## 1884 25.23351 25.23351
## 1885 25.11571 25.11571
## 1886 25.15248 25.15248
## 1887 25.19729 25.19729
## 1888 25.05286 25.05286
## 1889 25.11768 25.11768
## 1890 25.08710 25.08710
## 1891 24.99407 24.99407
## 1892 25.07019 25.07019
## 1893 25.01085 25.01085
## 1894 24.88515 24.88515
## 1895 24.95884 24.95884
## 1896 24.87469 24.87469
## 1897 24.84201 24.84201
## 1898 24.79420 24.79420
## 1899 24.62284 24.62284
## 1900 24.57259 24.57259
## 1901 24.54141 24.54141
## 1902 24.48421 24.48421
## 1903 24.39631 24.39631
## 1904 24.72686 24.72686
## 1905 24.62852 24.62852
## 1906 24.58852 24.58852
## 1907 24.58059 24.58059
## 1908 24.54271 24.54271
## 1909 24.52166 24.52166
## 1910 24.57541 24.57541
## 1911 24.59433 24.59433
## 1912 24.59905 24.59905Утицај различитих вредности за \(\alpha\) код експоненцијалног изравнања.
plot(rainseriesforecasts) 
plot(HoltWinters(rainseries, alpha=0.25, beta=FALSE, gamma=FALSE), main="Alpha=0.25")
plot(HoltWinters(rainseries, alpha=0.5, beta=FALSE, gamma=FALSE), main="Alpha=0.5")
plot(HoltWinters(rainseries, alpha=0.75, beta=FALSE, gamma=FALSE), main="Alpha=0.75")
plot(HoltWinters(rainseries, alpha=1, beta=FALSE, gamma=FALSE), main="Alpha=1")
library(forecast)
(rainseriesforecasts2 <- forecast(rainseriesforecasts, h=50))## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1913 24.67819 19.17493 30.18145 16.26169 33.09470
## 1914 24.67819 19.17333 30.18305 16.25924 33.09715
## 1915 24.67819 19.17173 30.18465 16.25679 33.09960
## 1916 24.67819 19.17013 30.18625 16.25434 33.10204
## 1917 24.67819 19.16853 30.18785 16.25190 33.10449
## 1918 24.67819 19.16694 30.18945 16.24945 33.10694
## 1919 24.67819 19.16534 30.19105 16.24701 33.10938
## 1920 24.67819 19.16374 30.19265 16.24456 33.11182
## 1921 24.67819 19.16214 30.19425 16.24212 33.11427
## 1922 24.67819 19.16054 30.19584 16.23968 33.11671
## 1923 24.67819 19.15895 30.19744 16.23724 33.11915
## 1924 24.67819 19.15735 30.19904 16.23479 33.12159
## 1925 24.67819 19.15576 30.20063 16.23235 33.12403
## 1926 24.67819 19.15416 30.20223 16.22991 33.12647
## 1927 24.67819 19.15257 30.20382 16.22748 33.12891
## 1928 24.67819 19.15097 30.20542 16.22504 33.13135
## 1929 24.67819 19.14938 30.20701 16.22260 33.13379
## 1930 24.67819 19.14778 30.20860 16.22016 33.13623
## 1931 24.67819 19.14619 30.21020 16.21773 33.13866
## 1932 24.67819 19.14460 30.21179 16.21529 33.14110
## 1933 24.67819 19.14301 30.21338 16.21286 33.14353
## 1934 24.67819 19.14142 30.21497 16.21042 33.14597
## 1935 24.67819 19.13982 30.21656 16.20799 33.14840
## 1936 24.67819 19.13823 30.21815 16.20556 33.15083
## 1937 24.67819 19.13664 30.21974 16.20312 33.15326
## 1938 24.67819 19.13505 30.22133 16.20069 33.15570
## 1939 24.67819 19.13346 30.22292 16.19826 33.15813
## 1940 24.67819 19.13188 30.22451 16.19583 33.16056
## 1941 24.67819 19.13029 30.22610 16.19340 33.16298
## 1942 24.67819 19.12870 30.22769 16.19097 33.16541
## 1943 24.67819 19.12711 30.22928 16.18855 33.16784
## 1944 24.67819 19.12552 30.23086 16.18612 33.17027
## 1945 24.67819 19.12394 30.23245 16.18369 33.17269
## 1946 24.67819 19.12235 30.23404 16.18127 33.17512
## 1947 24.67819 19.12077 30.23562 16.17884 33.17755
## 1948 24.67819 19.11918 30.23721 16.17642 33.17997
## 1949 24.67819 19.11760 30.23879 16.17399 33.18239
## 1950 24.67819 19.11601 30.24038 16.17157 33.18482
## 1951 24.67819 19.11443 30.24196 16.16915 33.18724
## 1952 24.67819 19.11284 30.24354 16.16673 33.18966
## 1953 24.67819 19.11126 30.24513 16.16431 33.19208
## 1954 24.67819 19.10968 30.24671 16.16189 33.19450
## 1955 24.67819 19.10810 30.24829 16.15947 33.19692
## 1956 24.67819 19.10652 30.24987 16.15705 33.19934
## 1957 24.67819 19.10493 30.25145 16.15463 33.20176
## 1958 24.67819 19.10335 30.25303 16.15221 33.20418
## 1959 24.67819 19.10177 30.25461 16.14980 33.20659
## 1960 24.67819 19.10019 30.25619 16.14738 33.20901
## 1961 24.67819 19.09861 30.25777 16.14496 33.21142
## 1962 24.67819 19.09704 30.25935 16.14255 33.21384plot(rainseriesforecasts2)
Видимо да су све прогнозе исте, једнаке последњој вредности нивоа. Овде су дате и границе интервала за ове оцењене вредности.
Границе интервала поверења се рачунају под претпоставком да су грешке независне, нормално расподењене са константном дисперзијом.
Позабавимо се сада серијом резидуала у циљу провере да ли је HW одговарајући моделна основу кога предвиђамо.
Потребно је проверити следеће:
- Да ли постоји серијска корелација у серији резидуала? Корелограм
acf(rainseriesforecasts2$residuals[-1], lag.max=20)
Осим тестом датим на графику корелограма, тестирање значајности корелација може се спровести још једним тестом. То је Ljung-Box тест који такође тестира да ли серија има статистички значајне корелације.
Box.test(rainseriesforecasts2$residuals, lag=20, type="Ljung-Box")##
## Box-Ljung test
##
## data: rainseriesforecasts2$residuals
## X-squared = 17.401, df = 20, p-value = 0.6268Нулта хипотеза је да су све корелације једнаке 0 и пошто је овде p-вредност теста велика прихватамо нулту хипотезу.
- Да ли су резидуали нормално расподељени?
Нормалност резидуала можемо проверити на више начина. На основу qqplot-a, где поредимо квантиле серије резидуала и теоријске квантиле нормалне расподеле.
qqnorm(rainseriesforecasts2$residuals[-1])
- Да ли је дисперзија серије резидуала константна током времена?
plot(rainseriesforecasts2$residuals[-1])
Може се закључити на основу графика. Са графика бисмо рекли да је можда ипак мања дисперзија на почетку, а већа касније.
\(\color{deepskyblue}{\text{Пример 3}}\)
skirts <- read.delim("skirts.txt")
attach(skirts) # Godisnja serija precnik suknji, od 1866 god## The following object is masked from rain:
##
## Yearskirts.ts <- ts(Diameter,start=c(1866))
plot(skirts.ts)
Са графика не видимо сезонску компоненту, али постоје неки краткорочни трендови, на почетку растућу, а касније опадајући. Применићемо на податке HW модел, само без параметара \(\gamma\) и сезонске конпоненте.
skirts.ts.hw <- HoltWinters(skirts.ts, gamma=FALSE)
skirts.ts.hw## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = skirts.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.8383481
## beta : 1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 529.308585
## b 5.690464Видимо да је \(\beta=1\), дакле нагиб се оцењује као разлика вредности серије у узастопним тренуцима. Рекли смо да је то добра оцена када је тренд непредвидив.
library(forecast)
skirts.ts.hw.forecast <- forecast(skirts.ts.hw, h=19)
skirts.ts.hw.forecast## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1912 534.9990 509.55210 560.4460 496.08130 573.9168
## 1913 540.6895 491.01052 590.3685 464.71204 616.6670
## 1914 546.3800 465.36129 627.3987 422.47258 670.2874
## 1915 552.0704 434.40205 669.7388 372.11216 732.0287
## 1916 557.7609 398.94120 716.5806 314.86713 800.6547
## 1917 563.4514 359.47147 767.4313 251.49103 875.4117
## 1918 569.1418 316.34076 821.9429 182.51596 955.7677
## 1919 574.8323 269.81480 879.8498 108.34829 1041.3163
## 1920 580.5228 220.10648 940.9391 29.31362 1131.7319
## 1921 586.2132 167.39191 1005.0345 -54.31870 1226.7452
## 1922 591.9037 111.82029 1071.9871 -142.32052 1326.1279
## 1923 597.5942 53.52019 1141.6681 -234.49517 1429.6835
## 1924 603.2846 -7.39593 1213.9652 -330.67069 1537.2399
## 1925 608.9751 -70.82861 1288.7788 -430.69495 1648.6451
## 1926 614.6655 -136.68903 1366.0201 -534.43211 1763.7632
## 1927 620.3560 -204.89720 1445.6092 -641.75986 1882.4719
## 1928 626.0465 -275.38060 1527.4736 -752.56728 2004.6602
## 1929 631.7369 -348.07309 1611.5470 -866.75319 2130.2271
## 1930 637.4274 -422.91396 1697.7688 -984.22478 2259.0796plot(skirts.ts.hw.forecast)
skirts.ts.hw.predict <- predict(skirts.ts.hw, n.ahead= 19)
ts.plot(skirts.ts, skirts.ts.hw.predict, lty = 1:2) 
acf(skirts.ts.hw.forecast$residuals[-(1:2)],lag.max=20) # kao argument acf mozemo zadati koliko vrednosti autokorelacione funkcije zelimo sa lag.max
Провера да ли је модел добар:
имамо једну статистички значајну корелацију, што је добро, јер једозвољено да \(5\%\) корелација упадне у критичну област, (тј. 1 корелација ако гледамо њих 20)
провера Ljung-Box тестом потврђује да је модел добар
Box.test(skirts.ts.hw.forecast$residuals, lag=20, type="Ljung-Box")##
## Box-Ljung test
##
## data: skirts.ts.hw.forecast$residuals
## X-squared = 19.731, df = 20, p-value = 0.4749plot(skirts.ts.hw.forecast$residuals[-(1:2)])
- Провера да ли је дисперзија константна током времена. На основу графика рекли бисмо да јесте.
\(\color{deepskyblue}{\text{Пример 4}}\)
data(AirPassengers) # serija analizirana na prvom casu
# Posto se cinilo da se sezonski efekat povecava kako trend raste, koristicemo multiplikativni HW model
AP <- AirPassengers
AP.hw <- HoltWinters(AP, seasonal = "mult")
AP.hw## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = AP, seasonal = "mult")
##
## Smoothing parameters:
## alpha: 0.2755925
## beta : 0.03269295
## gamma: 0.8707292
##
## Coefficients:
## [,1]
## a 469.3232206
## b 3.0215391
## s1 0.9464611
## s2 0.8829239
## s3 0.9717369
## s4 1.0304825
## s5 1.0476884
## s6 1.1805272
## s7 1.3590778
## s8 1.3331706
## s9 1.1083381
## s10 0.9868813
## s11 0.8361333
## s12 0.9209877plot(AP.hw)
AP.predict <- predict(AP.hw, n.ahead= 4 * 12)
ts.plot(AP, AP.predict, lty = 1:2) 
plot(AP.predict)
\(\color{deepskyblue}{\text{Задатак}}\)
- Учитати временску серију global (1856-2005) и нацртати график временске серије.
- Издвојити тренд, сезонску и случајну компоненту и нацратати корелограм.
- Уклопити временску серију у одговарајући HW модел.
- Користећи фитовани модел, предвидети вредности за период 2005-2010 и те вредности додати на график.