Analýza časových řad

Tyto příklady volně navazují na teoretickou prezentaci a demonstrují vysvětlovanou problematiku. Prakticky budeme pracovat jen v programu IBM SPSS Statistics. Potřebné podklady jsou ke stažení na disku Adelka/Public/Casoprostorove_analyzy.

Otevřete si soubor hrebiky.sav

Zobrazení časové řady

Data je potřeba si vždy nejdříve zobrazit v Sekvenčním grafu a prohlédnout si průběh. Zaměřte se na existenci trendu, heteroskedasticitu, přítomnost extrémů, sezónnost apod.

V SPSS naleznetu tuto volnu v nabídce Analyze/Forecasting/Sequence Charts.

Chybějící hodnoty

Časová řada musí být úplná, tedy nesmín obsahovat žádné chybějící hodnoty. Pokud obsahuje, je možné tyto údaje imputovat. Zkopírujte si sloupec s cenou hřebíků a vymažte si několik hodnot, které se pokusíme dopočítat.

Spusťte si nástroj Transform/Replace Missing Values

Z nabízených nástrojů není vhodný Series Mean (využívat jen v nouzi a pokud v datech není trend), lepší volbou je Nearby Mean. Linear Trend proloží časovou řadou přímku a a doplní podle ní chybějící hodnoty. Funguje jen tehdy, pokud je lineární trend, což většinou není. Navíc pracuje se všemi hodnotami. Linear interpolation spojuje jen sousední body. Vyzkoušejte všechny metody a posuďte kvalitu doplněných hodnot s originálními hodnotami. Zobrazte si všechny řady v jednom sekvenčním grafu a vyberte nejvhodnější metodu.

Testování heteroskedasticity

Pro hodnocení heteroskedasticity je potřeba si data rozdělit do skupin. Konkrétní pravidlo neexistuje, pokud řada není področní (náš případ) a ze sekvenčního grafu není jasná sezónnost, tak využijeme např. rozdělení do 10 skupin. V tomto případě stačí vydělit rok hodnotou 10 a hodnoty zaokrouhlit na celá čísla.

V nabídce Analyze/Explore pak zadejte cenu hřebíků a příslušenství do skupiny a v nabídce Plots vyberte Spread vs. Level with Levine’s Test. Pokud není trend v datech, tak budeme pracovat s untransformed, jinak je nutné data logaritmovat. V SPSS je veriabilita hodnocena IQR a centralita pak mediánem.

Vytvořte graf pro transformovaná i surová data a posuďte rozdíly mezi grafy. Levenův test má jako nulovou hypotézu skutečnost, že rozptyly jsou stejné (homoskedasticita). V našem případě však variabilita roste. Dle testu je ale v obou případech problém. Levenova statistika se ale snížila, což je dobře.

Vyhlazení časových řad

Vyhlazení časové řady se v SPSS realizuje v nabídce Transformation/Create Time Series a vybrat z nabízených metod Centered Moving Average. Nastavte vyhlazení s využitím 5, 10 a 30 kroky a vykreslete tyto tři vyhlazené řady spolu s originální časovou řadou. Porovnejte rozdíly a posuďte, zda se vyskytují části, kde by byl vývoj originální a vyhlazené řady v rozporu. Všimněte si, že vyhlazená řada je zkrácená.

Dekompozice časové řady a predikce s využitím regresního modelu

Pro predikci s využitím regresního modelu si otevřete data auta_benzin.sav. Časovou řadu (prodaná auta) si zobrazte a posuďte přítomnost trendu a sezónnosti.

Nadefinujte časovou složku v nabídce Data/Define Dates. Časová řada začíná v lednu/1960 a je patrný trend co 12 měsíců.

Dekompozice časové řady se provádí v nabídce Analyze/Forecasting/Seasional Decomposition. Je potřeba zvolit vhodný typ modelu. Posuďte, zda se jedná o aditivní či multiplikativní model. Jelikož aplituda se zdá být stejná, tak aditivní typ modelu. Pokud si nejsme jistí, je lepší otestovat heteroskedasticitu nebo vyzkoušet oba typy modelu a určit, který je vhodnější. Pro vyhlazení zvolte variantu Endpoints Weighted by 0.5.

Prohledněte si výsledek, vykreslete si sezónní faktory v grafu (odchylky v % od průměru). Rovněž si v grafu vykreslete všechny čtyři komponenty časové řady a zaměřte se na sezónní složku a trend.

V dalším kroku přistoupíme k modelování trendu, o který nám jde. Využijeme klasickou lineární regresi, tedy v nabídce Analyze/Regression/Curve Estimation a jako závislou proměnnou vybereme trendovou komponentu a jako nezávislou proměnnou pak zvolte položku Time. Můžeme vybrat různé křivky pro hledání závislosti, vyberte lineární, kvadratickou a exponenciální a nechte vykreslit výsledky. Ty zhodnoťme a vyberme nejvhodnější model.

V nabídce Save je možné uložit predikci, nastavte predikci na dva roky. Rozhodně není vhodné predikovat na dlouhá období. Následně zkopírujeme sezónní složku do prázdných buněk a sečteme predikovanou řadu a sezónní komponentu (aditivní model). V sekvenčním grafu následně zobrazíme originální řadu a predikci.

Všimněte si odchylky v 10/1962, jedná se o velkou odchylku, pro kterou je potřeba dělat intervenci? To je na zváženou, je potřeba každý zásah odůvodnit. Tato odchylka však ovlivnila všechny říjny v modelu.

Otevřete si data Trends Chapter 12 (od 1/1951 v měsících). Zobrazte si řadu a opět popište. Evidentně se jedná o multiplikativní model. Proveďte dekompozici s výběrem multiplikativního modelu. Vykreslete následně sezónní faktory a originální řadu. Všimněte si, že sezónnost hodně kopíruje originální data, tedy sezónnost není pravidelná. Odchylky v měsících jsou někdy malé, jindy velké, což je problém. Tato řada je již komplikovanější a není možné na ni použít tento jednoduchý přístup (podobně jako na téměř žádnou dnešní časovou řadu). Je tedy potřeba využít jiné metody (viz dále).

Vraťme se k datům auta_benzin.sav. Pokusíme se predikovat data čistě regresním modelem. Již víme, že je v datech trend, tedy musíme přistoupit k sezónnímu očišťování. Vytvoříme si časové řady sezónních indexů. V prvním kroku si vytvoříme novou proměnnou t. Ta bude nabývat hodnot 1 – N, tedy v nabídce Transform/Compute Variable vyberte funkci z All s názvem $CASENUM a vytvořte tuto proměnnou. Ta určuje jen posloupnost záznamů. Zkopírujte si rovněž proměnnou měsíc do nové proměnné. Pro vytvoření sezónních indexů si vybereme nástroj Data/Restructure, zvolíme druhou možnost, jako Identifier vybereme proměnnou t a jako Index nově zkopírovanou proměnnou s měsícem, další krok přeskočíme a na další kartě vybereme Create Indicator Variables a pojmenujeme ji M. Vytvoříme 12 nových sezónních proměnných s 0 a 1.

Následně přistoupíme k regresnímu modelování, v nabídce Analyze/Regression vybereme Linear Regression. Jako závislou proměnnou vybereme počet prodaných aut a jako nezávislou pak proměnnou t a všechny sezónní proměnné s výjimkou ledna, který si vybereme jako referenční, tedy M2 – M12. Posuďte kvalitu modelu a jeho parametry rovnice. Konstanta u t ozančuje měsíční nárůst prodeje a u jednotlivých měsíců pak změnu oproti lednu. Posuďte statistickou význmanost jednotlivých sezónních složek. Je patrné, že únor, srpen a prosinec nejsou statisticky významné a tedy ne různé od ledna. Nicméně tyto proměnné se neodstraňují (je to možné, ale nedělá se to). Uložte si Unstandardized Predicted Values.

Všimněte si, že v 9/1962 je evidentní nějaký problém. Podobné problémy mohou být identifikovány také analýzou reziduí (podstatná část regresního modelování obecně). Tedy uložíme si rovněž Unstandardized Residual Values a vykreslíme si sekvenční graf těchto reziduí. Není to největší odchylka v modelu, ale je velká. Pokud se objeví podobná odchylka, je vhodné (pokud je odůvodněná) přistoupit k intervenci. Tedy vytvoříme si novou proměnnou, která bude mít všude hodnoty 0 a v řádku 9/1962 pah hodnotu 1. Tuto proměnnou přidáme do modelu. Pokud bychom chtěli přidat další intervence, tak musíme přidat další proměnné, ale pozor na to. Ve zhodnocení modelu je možné posoudit, zda je daná intervence významná, je patrné, že ne, i když jen hraničně. Pozor na interpretaci hraniční hladiny významnosti.

Tento model lze ještě zpřesnit a to přidáním jiného než lineátrního trendu. Můžeme posoudit kvadratický trend, tedy vypočítáme novou proměnnou t2, tedy umocněná proměnná t a obě tyto proměnné spolu se sezónními faktory přidáme do modelu. Vidíme zlepšení, ale podobné komplikace modelu musí být odůvodněné (viz výše diskuze nad tvarem kvadratické křivky).

Pokud bych chtěl nyní predikovat do budoucna, musíme si protáhnout proměnnou t i všechny komponenty M1 – M12 a následně model spustit a predikce se vytvoří. Predikujte hodnoty o 2 roky dopředu s využitím lineárního modelu.

Vraťme se nyní zpět k datům Trends Chapter 12 a pokusme se vytvořit lepší model. Ze sekvenčního grafu je patrná změna sezónnosti mezi obdobími 1951-1963 a 1964-1968. Něco se zde stalo, nevíme co. Vytoříme tedy dvě sady sezónních faktorů pro tato dvě období (MA1 – MA12 a MB1 – MB12). První sada bude mít samé 0 v letech druhého období a naopak. Uvnitř daného období budou klasicky 1 a 0. Do modelu pak nezahrňte hodnotu MA1 a porovnejte výsledky. Tento model vychází lépe.

Rezidua je možné posoudit ještě uložením tzv. Studentizied Residuals a u těch která jsou mimo interval <-2;2> nebo <-3;3> je vhodné přistoupit k intervenci.

Hodnocení autokorelace reziduí

Spusťte si znova regresní model pro soubor benzin_auta.sav a uložte si unstandardized residuals. Zaměříme se na kontrolu reziduí v modelu. V nabídce Analyse/Forecasting/Autocorrelation vyberta proměnnou s rezidui a posuďte výsledek. Je patrné, že stále existují významné závislosti o 3 období nejvýznamnější, tedy v reziduích modelu jsou stále informace, které je ptřeba vytěžit. Zatím nemáme znalosti na nic lepšího, než přidat nějaké intervence.

Hodnocení autokorelace lze použít také pro originální časovou řadu, vyzkoušejme a zjistíme, že nejvyšší hodnota je ve 12. měsíci, tedy odpovídá sezónní složce. V tomto případě je to patrné, ale často jsou tyto závislosti skryty. Pak je potřeba odstranit trend, který je velmi silný a vyzkoušet na takto očištěných datech (v nabídce zapnout differences).

Je zároveň možné zvětšit také počet pozorování. V nabídce Options nastavíme 36. Jsou tak vidět opakované sezónnosti také v dalších sezónách, někdy je totiž patrná sezónnost i ob rok.

Exponenciální vyhlazování

Otevřete si soubor hrebiky.sav a zobrazte v sekvenčním grafu. Řada má patrný trend se změnou kolem roku 1900, který je dále následován mírným nárůstem. V prvním kroku se zbavíme/potlačíme heteroskedasticitu, tedy řadu zlogaritmujeme. Vyzkoušíme klasickou regresi na celou řadu, tedy Analyze/Regression/Curve Estimation a vybereme vhodné křivky (linear, quadratic, cubic, log). Posuďte kvalitu výsledných modelů, neřiďte se jen hodnotami v tabulkách, ale posuďte také vývoj křivek v grafu. Žádná z křivek není vhodná, navíc příliš historické záznamy nemají již pro vývoj příliš význam.

Přistoupíme tedy k exponenciálnímu vyhlazování. Klikněte na Analyze/Forecasting/Create Traditional Models. Jako závislou proměnnou vybereme cenu hřebíků a jako nezávislou zatím nic (v klasickém modelu to byl čas). V nabídce Method vybereme Exponential Smoothing a jako Criteria pak Holt’s Model. Kromě Brownova modelu se zde vyskystuje ještě možnost Damped Trend, což je model s potlačeným trendem. V praxi se moc nepoužívá. Dále v nabídce Statistics vyberte Parameter Estimates, v nabídce Plots: Fit Values, ACF a PACF (For Individual Models). Model spusťte. Ve výsledku se zaměříme jen na tabulku Exponential Smoothing Model Parameters, tedy konkrétně na hodnoty alfa (vliv předchozí hodnoty) a gama (trend za celé období). Výsledek je extrémní, prakticky vliv má jen předchozí hodnota, což potvrzuje také sekvenční graf, kdy je model posunut o jedno období dopředu = primitivní model. Rovněž u autokorelace je vidět výrazný vliv pro 2. období. V modelu zvolte na záložce Save uložit Predicted Values a na kartě Options vyberte …though specific date (dáme rok 2010). Predikován je meziroční pokles o 0,00387 (Transformation/Create Time Series/Difference), stále stejný pokles. Model stále není ideální, později se k němu ještě vrátíme.

Pozn.

Alpha Smoothing Constant – When the Search Method is set to Specified Value, this option specifies the value of alpha. Alpha is the smoothing constant for the level of the series. The limits of this value are zero and one. Usually, a value between 0.1 and 0.3 are used. As the value gets closer to one, more and more weight is given to recent observations. 

Gamma Smoothing Constant – When the Search Method is set to Specified Value, this option specifies the value of beta. Gamma is the smoothing constant for the trend. The limits of this value are zero and one. Usually, a value between 0.1 and 0.3 are used. As
the value gets closer to one, more and more weight is given to recent observations.

Delta Smoothing Constant – When the Search Method is set to Specified Value, this option specifies the value of gamma. Delta is the smoothing constant for the seasonal factors. The limits of this value are zero and one. Usually, a value between 0.1 and 0.3 are used. As the value gets closer to one, more and more weight is given to recent observations.

Další soubor, který použijeme je Trends Chapter 12 a vyzkoušíme exponenciální vyhlazování s trendem i sezónností. Nadefinujeme časovou složku a vytvoříme nový model. Pod tlačítkem Criteria máme nové možnosti – Simple Seasonal (bez trendu, jen sezónnost, málo časté) a dva typy Winter’s models, vyzkoušet můžeme oba a porovnat výsledky. Ostatní nastavení zůstává stejné. Posuďte opět vyrovnávací konstanty, kdy hodnota delta je vysoká, tedy významný vliv nedávných sezón (v roce 1964 je vidět zajímavý rozdíl). Predikce se doporučuje dělat na kratší období, v tomto případě maximálně na 1 rok, aby se chyba odhadu nerozšiřovala do budoucna, a dělat predikce častěji např. každý další měsíc. Prozkoumejte autokorelaci v reziduích, jak je patrné, tak stále v reziduích existuje nevysvětlená informace.

Model ARIMA

Vytvoření modelu ARIMA je často otázkou zkušeností a umění, my si to vyzkoušíme na datech s názvem Trends Chapter 6. V těchto datech nevíme časové měřítko, ale to nevadí, jde nám o posloupnost dat. Náš model si otestujeme tzv. pseudopředpovědí, tedy vytvoříme model pro většinu dat a poslední záznamy si “odtrhneme” a ověříme na nich kvalitu předpovědi. Pokud pseudopředpověď sedí, tak aplikujeme model na všechna data. Vybereme nástrojem Select Cases prvních 100 záznamů a ostatní “zapomeneme”. Data si následně zobrazíme v sekvenčním grafu, posoudíme vývoj řady. Zdá se, že v datech není trend a hodnoty kolísají náhodně.

V nabídce Analyze/Forecasting/Create Traditional Models vybereme možnost ARIMA, začínáme s prázdným modelem, tedy všechny parametry modelu budou nulová, další nastavení zvolíme stejné jako v případě exponenciálního vyhlazování. Výstupem je konstanta (není integrovaný proces), proto linie. Podle korelogramů je potřeba vybrat vhodný typ modelu, postupuje se metodou pokus-omyl a obecně dle instrukcí na přednášce. Podle tvaru korelogramu volíme v prvním kroku model MA, tedy parametr MA = 1. Opět spustíme model a posoudíme korelogramy, je vidět, že 10. zpoždění je nad hranicí. Interpretujeme také Box-Ljungovu statistiku (v SPSS se měří v 18. zpoždění) a ta testuje nulovou hypotézu, že všechny autokorelace pro předchozí zpoždění jsou nulové. V našem případě je významnost <0,05 a tedy přijímáme alternativní hypotézu. Posoudíme kvalitu modelu podle stacionárního koeficientu determinace (0,37), je možné použít i další ukazatele. V grafu je vhodné zobrazit i intervaly spolehlivosti naší předpovědi a to v nabídce Plots zapnout Confidence Intervals for Forecast. V grafu si zobrazíme celou řadu, předpověď i intervaly spolehllivosti a posoudíme výsledek. Jak je vidět z dat, je zde patrný výkyv u záznamů 160–166, zde tedy zavedeme intervenci (hodnota 1 u těchto hodnot, jinak 0).

Nyní zvolíme kompletní data a vytvoříme prázdný model, zatím bez intervence a posuďte korelogramy. Je vidět, že korelogramy jsou komplikované a to z důvodu chyb v záznamech 160–166. Vždy je potřeba data kontrolovat pro přítomnost anomálních hodnot. Spustíme tedy prázdný model ARIMA (0,0,0) s itervencí a již je vidět opět MA proces, tedy další model ARIMA (0,0,1). Korelogramy jsou již v pořádku. V Options zadáme predikci do řádku 220, u stacionárních řad je možné volit delší predikce. V proměnné intervence musíme prodloužit hodnotu 0 do řádku 220 a následně spustíme model. Do datové matice jsou vloženy predikované hodnoty. Jedná se o jednoduchou lineární predikci. Zobrazíme opět model a intervaly spolehlivosti, je vidět, že predikce vypadá již lépe.

Další datovým souborem, který vyzkoušíme je soubor Trends Chapter 10 – tržní podíly zubních past, kdy od 136. týdne začala reklamní kampaň. Jako obvykle si v prvním kroku vykreslíme časovou řadu, kde patrný vliv kampaně (pokud bychom o ní nevěděli, je potřeba zjistit o datech maximum informací) a tedy je potřeba zavést intervenci od řádku 136+ (Compute Variable: $Casenum>=136). Spustíme klasicky model ARIMA (0,0,0) a zhodnotíme korelogramy. Je zde patrný trend, tedy nastavíme integrovaný proces a tedy ARIMA (0,1,0), pozor, musíme diferencovat také všechny nezávislé proměnné, v tomto případě intervenci (záložka Transfer Function pod tlačítkem Criteria). Korelogramy jsou lepší, dalším krokem by bylo vhodné zkusit MA, tedy ARIMA (0,1,1). Posuďte výsledné korelogramy, Ljung-Box statistika je větší než 0,05 a tedy autokorelace všech kroků jsou nulové. Testovány jsou i přínosy proměnných (ARIMA parameters), zde je patrné, že konstanta je zbytečná a je vhodné ji tedy vypnout (Criteria/Model a odškrtnout Include Constant in Model). Numerator v Lag 0 ukazuje 0,162 nárůst mezi 2 řadami (před a po intervenci), což odpovídá nárůstu o 16,2 % (vliv kampaně). Pro predikci musím opět protáhnout intervenci a to hodnotou 1, pokud budu v kampani pokračovat nebo hodnotou 0, pokud nebudu, protáhnout musím i weeks do hodnoty 300. Zobrazíme si predikci a to včetně intervalů spolehlivosti, které se rozpínají vzhledem k zakomponovanému integrovanému procesu.

Zkusíme vytvořit ještě jiný model a to bez integrovaného procesu, tedy zkusíme ARIMA (0,0,1) a pak (1,0,1). Porovnáme oba modely, z hlediska statistiky jsou oba modely v pořádku, ale logicky je vhodnější predikce modelu s integrovaným procesem. Je možné vyzkoušet také pseudopredikci – samostatný úkol.

Model SARIMA

Model SARIMA si vyzkoušíme na souboru registrace.sav. Postupovat budeme k již známých krocích:

1. data si zobrazíme v sekvenčním grafu, všimněte si, že sezónní složka není patrná, ale v případě měsíčních dat je téměř vždy
2. vypočítáme si autokorelaci a zobrazíme korelogramy, kdy zvětšíme Number of Lags na 65 (studuji sezónnost, tedy v tomto případě 5 let a kousek)
3. ACF – vykazuje pozvolný pokles a jasně viditelné sezóny, aby byla sezónnost jasně patrná, zaškrtneme differences a sezónnost je již jasně patrná
4. v options je možné zaškrtnout Display Autocorrelation at priodic lags
5. bude nutná diference (pozvolný pokles ACF), tedy v prvním kroku ARIMA (0,1,0), jinak nastavení klasické, ale v Options nastavíme Max. N of lags 65
6. zaměříme se nejdříve na sezónní část, kterou je vhodné vymodelovat nejdřív a následně se zaměřit na nesezónní část modelu. ACF postupně klesá, ale ne od 1 pro integrální, možná AR (exponenciální pokles), ale nejdříve zkusíme pro sezónní složku diferenci, tedy rovněž (0,1,0)
7. v dalším kroku vypadá vhodné modelovat sezónnost procesem MA – ACF jedna významní a PACF postupně klesá, tedy (0,1,1), toto již vyčistilo sezónnost
8. teď se již zaměříme na nesezónní část modelu, vidíme první dvě významná pozorování, takže asi budeme muset volit 2 do modelu, ale půjdeme postupně metodou pokus-omyl. Začneme AR (lépe se interpretuje), tedy (1,1,0)(0,1,1). Odstranil se první statisticky významný prvek, ale druhý zůstal, proto zkusíme (2,1,0)(0,1,1). Tento model již vypadá slušně, opíšeme si hodnoty kvality modelu.
9. Budeme ale testovat dále a místo procesu AR vybereme proces MA a vyzkoušíme (0,1,2)(0,1,1), ale výsledky jsou horší.
10. Zkusíme další model úplně od začátku, kdy v sezónní složce zkusíme nepoužít integrovaný proces (na začátku jsme váhali) a vyzkoušíme (0,1,0)(1,0,0) a následně pak (0,1,0)(2,0,0), sezónní je již vyřešena a vypadá lépe než v předchozím případě.
11. Dalším zlepšováním modelu dojdeme k výsledku (2,1,0)(2,0,0) a porovnáme s modelem (0,1,2)(2,0,0)
12. Jak vybrat ten nejvhodnější ze čtyř výsledných modelů.
    1. Model (0,1,2)(0,1,1) – zavržen z důvodu nízké Box-Ljungovy statistiky
    2. Model (0,1,2)(2,0,0) – rovněž nevhodný, jelikož proces MA se špatně interpretuje a BIC je vyšší
    3. Modely (2,1,0)(0,1,1) a (2,1,0)(2,0,0) budeme dále analyzovat
13. v prvním kroku otestujeme normalitu reziduí – v Options vybereme možnost Noise Residuals pro oba modely a v nabídce Explore, kde je patrný vyšší průměr, ale také extrémní hodnota 96 (kandidát na intervenci), další odlehlé hodnoty již nejsou extrémní. Podobně také testy nenaznačují normální rozdělení.
14. Zavedeme danou intervenci pro řádek 96 a spustíme oba modely. Pozor, intervenci musíme sezónně i nesezónně diferencovat! V modelu je vidět, že intervence je statisticky významná. Porovnáme oba modely, ale oba jsou stále dobré.
15. Vyzkoušíme pseudopredikci pro jeden rok, ale i zde je vidět, že oba modely predikují podobný vývoj.
16. poslední rok v datech se chová jinak, po poklesu 4 roky se 5. rok zvedl, což model neočekával a nepracuje s tím dobře.

Vyzkoušíme ještě soubor hrebiky.sav, který se nám nepodařilo dobře predikovat. Zobrazíme a vidíme, že je nutné diferencovat, tedy (0,1,0), dále model AR(2,1,0), lag 1 je ale nevýznamný a korelace 4 je vysoký. Vyzkoušíme MA (0,1,2), který je lepší, ale lag 1 je opět nevýznamný. Čistota autokorelační funkce je primární. Mohu využít neúplné procesy, tedy vynechám nějaký krok (zpoždění), bohužel toto se musí provést v syntaxy, kde smažu 1 v MA=[1,2]. Tím, nastavím korektně model a vyhodnotím.