MASARYKOVA UNIVERZITA Přírodovědecká fakulta
Diplomová práce BAYESOVSKÉ DYNAMICKÉ MODELY
vypracovala: Bc. Tereza Stárková vedoucí práce: RNDr. Martin Kolář, Ph.D.
Brno 2008
Prohlašuji, že jsem tuto diplomovou práci vypracovala samostatně pod odborným vedením RNDr. Martina Koláře, Ph.D. a uvedla jsem všechny prameny a publikace, ze kterých jsem čerpala. V Brně dne 12. května 2008 Bc. Tereza Stárková
Chtěla bych poděkovat vedoucímu mé diplomové práce RNDr. Martinu Kolářovi, Ph.D. Děkuji za trpělivost, pečlivost a vstřícnost a za cenné rady, bez kterých by tato práce nemohla být napsána.
Obsah Úvod
5
1 Základní pojmy 1.1 Bayesovská teorie . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Základní značení . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Dynamické lineární modely . . . . . . . . . . . . . . . . . . . .
6 6 7 7
2 Polynomiální model prvního stupně 2.1 Dynamický lineární model a rekurentní vztahy 2.2 Konstantní model . . . . . . . . . . . . . . . . 2.3 Diskontní faktor a volba W . . . . . . . . . . 2.4 Limitní tvar modelu - ARIMA(0,1,1) . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
9 10 13 16 17
3 Dynamický regresní model
19
4 Dynamické lineární modely (DLMs)
27
5 Referenční analýza DLMs
32
6 Sezónní modely 34 6.1 Fourierovo vyjádření sezónnosti . . . . . . . . . . . . . . . . . 34 6.2 Modely plné sezónnosti ve Fourierově vyjádření . . . . . . . . 40 7 Rozšíření standardních DLMs 42 7.1 Polynomiální model druhého stupně . . . . . . . . . . . . . . . 42 7.2 Trendově-sezónní model . . . . . . . . . . . . . . . . . . . . . 43 Seznam použité literatury
48
Úvod Matematický model je jakési schéma, které pomocí matematického jazyka popisuje určitý systém. Matematické modely se využívají v přírodních vědách, v inženýrství, ekonomii, sociologii i politologii. Ve své diplomové práci jsem se zaměřila na třídu dynamických modelů, které pro své výpočty užívají Bayesovskou teorii. V úvodní kapitole si vysvětlíme, co znamená ”Bayesovský přístup”, jaké značení budeme používat, a uvedeme základní rysy Bayesovského dynamického modelování. Druhá kapitola představí polynomiální model prvního stupně. Zavedeme nejjednodušší definici dynamického lineárního modelu a ukážeme jeho speciální případ, tzv. konstantní model. Na závěr uvedeme vztah s Box-Jenkinsovými ARIMA procesy. Ve třetí kapitole se zaměříme na dynamický regresní model. Použijeme model vytvořený v programu Maple na různá data a budeme interpretovat získané výsledky. Ve čtvrté kapitole definujeme obecný dynamický lineární model a podrobněji se pak budeme věnovat jeho skalárnímu tvaru. Pátá kapitola nás v krátkosti seznámí s referenční analýzou, která je velmi užitečná pro určení počátečních hodnot v modelu. V šesté kapitole si podrobně představíme sezónní modely a ukážeme příklad takového modelu v programu Maple. Dále definujeme dynamické lineární modely vytvořené principem superpozice. V závěrečné kapitole zkombinujeme dva různé modely - model s trendovou složkou (polynomiální model druhého stupně) a model sezónní a použijeme tento model na reálná data. Na závěr pak porovnáme výsledky dvou různých modelů.
5
Kapitola 1 Základní pojmy 1.1
Bayesovská teorie
Bayesovská teorie je postavena na základním předpokladu, že všechny nejisté veličiny mohou být vyjadřovány a poměřovány pomocí pravděpodobností. Bayesovský přístup nabízí ucelený pohled na proces učení a nezávisí na žádných speciálních předpokladech. Klíčovým elementem takového přístupu je Bayesova věta. Uvedeme si ji pro spojité náhodné veličiny. Nechť p(Y, θ) označuje pravděpodobnostní hustotu náhodné veličiny pozorování Y a parametru θ, který také považujeme za náhodný. Potom p(Y, θ) = p(Y |θ)p(θ) = p(θ|Y )p(Y ) a tedy p(θ|Y ) =
p(θ)p(Y |θ) , p(Y )
kde p(Y ) 6= 0. Poslední výraz nám dává přímou úměrnost s konstantou úměrnosti p(Y1 ) : p(θ|Y ) ∝ p(θ)p(Y |θ), kde p(θ|Y ) je posteriorní pravděpodobnostní hustota θ (posterior ) podmíněná daty Y , p(θ) je apriorní pravděpodobnostní hustota θ (prior ) a p(Y |θ) je věrohodnostní funkce proměnné θ. Prior p(θ) je na datech nezávislý a shrnuje vše, co o θ víme ještě dříve, než obdržíme data. Funkce věrohodnosti vyjadřuje všechny informace obsažené v pozorování Y . Posterior p(θ|Y ) shrnuje vše, co o θ víme po prozkoumání dat, kombinuje datové i nedatové informace.
6
1.2
Základní značení
Modelujeme časovou řadu veličin vyjádřených reálnými čísly. Obecnou hodnotu značíme Y , hodnotu v čase t potom Yt . Budeme předpokládat, že pozorování začínají v čase t = 1 a řada se tedy rozvíjí jako Y1 , Y2, . . . . Neděláme rozdíl v pojmech náhodná veličina nebo realizace náhodné veličiny. Předtím, než pozorujeme hodnotu řady v čase t, Yt označuje neznámou náhodnou veličinu, kterou poznáme právě pozorováním. Ve všech případech uvažujeme normální rozdělení se spojitou hustotou p(.). Jestliže je náhodná veličina podmíněna jistou množinou informací nebo událostmi, zobrazíme to svislou čarou. Např. závisí-li veličina Y na množině informací D, pak její hustotu píšeme ve tvaru p(Y |D). Střední hodnotu značíme E[.], rozptyl V [.]a kovarianci C[., .]. Používáme římská písmena k označení veličin, které buď známe, nebo jsou neznámé, ale můžeme je pozorovat. Náhodné veličiny, které jsou neznámé a nepozorovatelné, popisujeme řeckými písmeny a představují pro nás parametry.
1.3
Dynamické lineární modely
Matematikové či statistici často modelují časové řady pomocí dynamických modelů. Pojem ”dynamický” znamená, že se procesy v čase mění. Nejznámějšími a nejpoužívanějšími modely jsou dynamické lineární modely s normálním rozdělením. Za předpokladu normality budeme jednoduše mluvit o dynamických lineárních modelech (DLMs). Bayesovský přístup k předpovídání a modelování zahrnuje i) flexibilitu modelu, ii) strukturaci pomocí parametrických modelů, iii) informace o parametrech dané pravděpodobnostmi, iv) předpovědi odvozené jako pravděpodobnostní hustoty. Flexibilita modelu a strukturace jsou pro modely časových řad typické. Jak se čas vyvíjí, dostáváme stále nové informace potřebné k předpovídání a k případnému poopravení modelu. Počáteční množina informací D0 představuje všechny dostupné informace, které na začátku modelování ovlivňují náš úsudek o modelu. Zahrnuje jeho historii a všechny veličiny, které jej určují. Obecně, v jakémkoliv čase t > 0 jsou naše rozhodnutí a předpovědi podmíněny množinou informací Dt . Ta obsahuje jak předchozí množinu informací
7
Dt−1 , tak nové pozorování Yt , Dt = {Yt , Dt−1 }. Jestliže obdržíme další relevantní informace, píšeme Dt = {It , Dt−1 }. Pro jednokrokovou předpověď je náš názor vyjádřen pomocí parametrického modelu, p(Yt |θ t , Dt−1 ), kde θ t je určující vektor parametrů v čase t. Tento matematicko-statistický zápis umožňuje komunikaci mezi matematikem, který model vytvořil, a klientem, který jej využívá.
8
Kapitola 2 Polynomiální model prvního stupně Polynomiální model prvního stupně je nejjednodušší netriviální model, na kterém se dají dobře vysvětlit důležité pojmy a obecně chování dynamických lineárních modelů. Pozorovaná veličina Yt je tvaru Yt = µt + νt , kde µt je aktuální hodnota veličiny v čase t a νt ∼ N[0, Vt ] je pozorovací chyba. Vývoj hodnoty veličiny v čase splňuje rovnici jednoduché náhodné procházky µt = µt−1 + ωt , kde ωt ∼ N[0, Wt ] nazýváme vývojovou chybou. Předpokládáme, že obě chyby mají pro každé t normální rozdělení. Dále předpokládáme, že pro každé i, j, i 6= j, jsou νi a νj nezávislé stejně jako ωi a ωj a že také νi a ωj jsou vzájemně nezávislé. Dále předpokládáme, že rozptyly Vt a Wt jsou známé pro každé t. Obrázky 2.1 a 2.2 ukazují příklady takových řad Yt (zeleně) společně s hodnotami µt (červeně). V obou případech je počáteční hodnota µ0 = 25 a rozptyly Vt = V = 1, Wt = W jsou v čase konstantní. Takový model zapíšeme ve tvaru (Yt |µt ) ∼ N[µt , Vt ] (µt |µt−1 ) ∼ N[µt−1 , Wt ] pro každé t = 1, 2, . . . . Tento model je v praxi dosti rozšířený. S úspěchem se využívá v mnoha aplikacích, zvláště u krátkodobých předpovědí např. pro plánování výroby nebo při regulaci zásob. 9
28
27
26
25
24
23
22 0
20
40
60
80
100
Obr. 2.1: Příklad řady pro W = 0.05
2.1
Dynamický lineární model a rekurentní vztahy
Připomeňme, že pro každé t (včetně t = 0) shrnuje Dt všechny dostupné informace o minulých hodnotách veličin, Dt = {Yt , Dt−1 }. Definice 2.1. Pro každé t definujeme dynamický lineární model (DLM) následovně pozorovací rovnice: rovnice systému: počáteční informace:
Yt = µt + νt , νt ∼ N[0, Vt ], µt = µt−1 + ωt , ωt ∼ N[0, Wt ], (µ0 |D0 ) ∼ N[m0 , C0 ],
kde posloupnosti chyb νt a ωt jsou nezávislé, vzájemně nezávislé a nezávislé na (µ0 |D0 ). Střední hodnota m0 je počáteční odhad hodnoty µ0 a C0 je její rozptyl. Klíčovými složkami modelu jsou rozdělení shrnutá v následující větě. Časový vývoj informací o modelu zachycují uvedené rovnice. 10
32
30
28
26
24
0
20
40
60
80
100
Obr. 2.2: Příklad řady pro W = 0.5 Věta 2.1. Pro každé t > 0 máme následující jednokrokovou předpověď a posteriorní rozdělení: a) posterior pro µt−1 : (µt−1 |Dt−1 ) ∼ N[mt−1 , Ct−1 ] pro nějakou střední hodnotu mt−1 a rozptyl Ct−1 ; b) prior pro µt : (µt |Dt−1 ) ∼ N[mt−1 , Rt ], kde Rt = Ct−1 + Wt ; c) jednokroková predikce: (Yt |Dt−1 ) ∼ N[ft , Qt ], kde ft = mt−1 a Qt = Rt + Vt ; d) posterior pro µt : (µt |Dt ) ∼ N[mt , Ct ], kde mt = mt−1 + At et a Ct = At Vt , kde At = Rt /Qt a et = Yt − ft . 11
Důkaz. Důkaz provedeme indukcí. Pro t = 1 máme (µ0 |D0 ) ∼ N[m0 , C0 ]. Podle definice µ1 = µ0 + ω1 , ω1 ∼ N[0, W1 ]. Odtud dostáváme (µ1 |D0 ) ∼ N[m0 , C0 + W1 ] a označíme R1 = C0 + W1 . Jednokrokovou předpověď Yt získáme jako Y1 = µ1 + ν1 , ν1 ∼ N[0, V1 ] a máme tedy (Y1 |D0 ) ∼ N[m0 , R1 + V1 ], kde označíme f1 = m0 a Q1 = R1 + V1 . Odvození posteriorního rozdělení µ1 si ukážeme obecně dále v důkazu. Nyní předpokládejme, že platí rozdělení uvedené v a) pro libovolné t > 1. Potom µt je součet dvou nezávislých normálně rozdělených náhodných veličin µt−1 a ωt . Má tedy normální rozdělení. Střední hodnotu a rozptyl získáme sečtením středních hodnot a rozptylů daných sčítanců, (µt |Dt−1 ) ∼ N[mt−1 , Rt ], kde Rt = Ct−1 + Vt . To nám dává b). Podobně Yt je součtem nezávislých normálně rozdělených veličin µt a νt , tedy má normální rozdělení (Yt |Dt−1 ) ∼ N[mt−1 , Qt ], kde Qt = Rt + Vt , což je právě c). Jakákoliv lineární funkce Yt a µt bude lineární kombinací nezávislých normálně rozdělených veličin νt , ωt a µt−1 . Bude tedy mít normální rozdělení a podle definice má (Yt , µt |Dt−1 ) dvojrozměrné normální rozdělení. Abychom určili vektor středních hodnot a varianční matici, uvědomme si, že v b) a v c) jsme již příslušné marginální hustoty určili. (Yt |Dt−1 ) ∼ N[ft , Qt ] a (µt |Dt−1 ) ∼ N[mt−1 , Rt ] a pro kovarianci platí C[Yt , µt |Dt−1 ] = C[µt + νt , µt |Dt−1 ] = V [µt |Dt−1 ] = Rt , což získáme aditivností a nezávislostí µt a νt . Dvojrozměrné rozdělení pak zapíšeme ve tvaru Qt Rt Yt Dt−1 µt−1 . , ∼N Rt Rt µt−1 µt Dt−1 12
Korelace ρt = Rt /(Rt Qt )1/2 je kladná, přičemž ρ2t = Rt /Qt = At . Rozdělení (µt |Yt , Dt−1 ) je normální se střední hodnotou mt−1 + ρ2t (Yt − mt−1 ) = mt a rozptylem (1 − ρ2t )Rt = Rt Vt /Qt = At Vt = Ct . Tím jsme dokázali d). et je chyba předpovědi v následujícím kroku, je to rozdíl mezi pozorovanou hodnotou Yt a očekávanou hodnotou ft . At má roli apriorního regresního koeficientu, v tomto případě je to čtverec jejich korelačního koeficientu. Díky tomu, že každá složka modelu má normální rozdělení, dají se výsledky spočíst jednoduše a elegantně. Poznamenejme, že mt lze vyjádřit i jinak, a to mt = At Yt + (1 − At )mt−1 , což nám ukazuje, že mt je vážený průměr apriorního odhadu hodnoty mt−1 a pozorování Yt . Váha, neboli tzv. adaptivní koeficient, At leží mezi 0 a 1, přičemž je blíže 0, jestliže více věříme apriorní hustotě nežli samotnému pozorování (Rt < Vt ), a blíže 1, jestliže je aprior méně informativní, jestliže dáváme větší váhu pozorování.
2.2
Konstantní model
Speciální případ modelu, ve kterém mají pozorovací a vývojová chyba konstantní rozptyl, se nazývá konstantní model. Poněvadž množina informací Dt = {Yt , Dt−1 } pro každé t > 0 neobsahuje o časové řadě žádné vnější informace, říkáme o modelu, že je uzavřený. Uzavřený konstantní DLM má V a W kladné, konečné a konstantní. Prior pro t = 0 je (µ0 |D0 ) ∼ N[m0 , C0 ] a pro t > 0 je Dt = {Yt , Dt−1 }. Důležitou konstantou je poměr r = W/V , který v inženýrské terminologii hraje roli poměru signál . Měří relativní rozptyl systémové a pozorovací rovšum nice. Příklad 1. Farmaceutická společnost dodává na trh lék Dobropyrin, kterého se průměrně za měsíc prodá 100 balení. Vedení se rozhodlo změnit obal výrobku a očekává nyní vyšší odbyt. ”Nový” Dobropyrin nahradí ten starý v lednu, t = 1, přičemž zůstane stejný jeho název a také jeho cena. Vedení požaduje krátkodobou předpověď budoucí poptávky. Pacienti užívají lék pravidelně, poptávka tedy bude v čase lokálně konstantní 13
a pro vyjádření celkového měsíčního prodeje Dobropyrinu použijeme polynomiální model prvního stupně. Očekáváme, že kolísání prodeje a pozorované odchylky značně převýší odchylku hodnoty poptávky, tedy ωt bude v porovnání s νt mnohem menší. U modelu pro starý Dobropyrin se s úspěchem předpovídalo pro V = 100 a W = 5 a stejné hodnoty si ponecháme i pro model nový. V prosinci (t = 0) předpokládáme, že v novém roce poptávka naroste přibližně o 30% na 130 balení za měsíc. Stejně tak předpokládáme, že se nesníží více než o 10 a nezvýší více než o 70 balení. Tedy rozmezí 80 balení představuje 95% odchylku pro µ0 , což jsou 4 směrodatné odchylky σ, a C0 = (80/4)2. Počáteční informace je tedy tvaru (µ0 |D0 ) ∼ N[130, 400]. Pracovní model prodeje Yt v měsíci t je Yt = µt + νt , µt = µt−1 + ωt , (µ0 |D0 ) ∼ N[130, 400].
νt ∼ N[0, 100], ωt ∼ N[0, 5],
Zde r = 0.05 je velmi malý, což je typické pro tento druh aplikace. Následující tabulka zobrazuje pozorování a hodnoty dalších složek modelu během prvních devíti měsíců. Měsíc t 0 1 2 3 4 5 6 7 8 9 10
Předpovídané rozdělení Qt ft 505 185 151 139 133 130 128 127 126 125
130.0 146.0 141.4 141.9 145.3 142.6 143.9 140.4 142.2 143.0
Adaptivní koeficient At
Pozorování
Chyba
Yt
et
0.80 0.46 0.34 0.28 0.25 0.23 0.22 0.21 0.21 0.20
150 136 143 154 135 148 128 149 146
20.0 -10.0 1.6 12.1 -10.3 5.3 -15.9 8.6 3.8
Posteriorní informace mt Ct 130.0 400 146.0 80 141.4 46 141.9 34 145.3 28 142.6 25 143.9 23 140.4 22 142.2 21 143.0 20
Z počátku vidíme, že rozptyl Ct je relativně velký, což zobrazuje naši počáteční nejistotu, o kolik se zvýší či sníží prodej. Jakmile ale dostáváme data 14
o prodeji, rychle tato hodnota klesá. Stejně tak adaptivní koeficient At je zpočátku roven 0.8, což znamená, že větší váhu dáváme hodnotám pozorování, později se však dostáváme až na limitní hodnotu 0.2, která reprezentuje jistotu našich předpovědí. Později si ukážeme, že At v tomto případě konverguje k hodnotě 0.2. Jednou z hlavních výhod bayesovského přístupu je, že můžeme jednoduše kombinovat modelové informace, data a informace, které dostaneme až během dalšího vývoje. V příkladu pro t = 9 máme (µ9 |D9 ) ∼ N[143, 20] s jednokrokovou predikcí (Y10 |D9 ) ∼ N[143, 125]. Předpokládejme nyní, že jsme v čase t = 9 obdrželi informace I9 o tom, že konkurenční firma stahuje z trhu svůj lék Zlodafen kvůli nežádoucím vedlejším účinkům. To znamená, že v čase t = 10 přejdou pacienti, kteří užívali Zlodafen, na Dobropyrin, a je známo, že uživatelé Zlodafenu tvoří asi 50% pacientů. Tedy poptávka po Dobropyrinu by měla vzrůst o 100% (m10 = 286). Naše nejistota zde je ale značná. Model se v čase t = 10 změní a (ω10 |D9 , I9 ) ∼ N[143, 900], což dává prior (µ10 |D9 , I9 ) ∼ N[286, 920]. Nová predikce bude (Y10 |D9 , I9 ) ∼ N[286, 1025]. Tím pádem opět vzroste adaptivní koeficient z 0.2 na 0.9 a celý proces se bude opakovat. Obr. 2.3 ukazuje průběh pozorování a jednokrokových predikcí včetně změny v čase t = 10, na obr. 2.4 jsou pak vykresleny hodnoty adaptivního koeficientu. Míra adaptace na nová data (vyjádřená adaptivním koeficientem At ) v uzavřeném modelu rychle konverguje ke konstantní hodnotě. Věta 2.2. Pro t → ∞ platí At → A a q 4 r 1+ r −1 . kde A = 2
Ct → C = AV,
Důkaz. Jelikož Ct = At V a 0 < At < 1, je 0 < Ct ≤ V pro každé t. Užitím rekurentních vztahů Ct−1 = Rt−1 + V −1 a Rt = Ct−1 + W dostáváme −1 −1 −1 −1 Ct−1 − Ct−1 = Rt−1 − Rt−1 = Kt Ct−1 , − Ct−2
kde Kt = Ct−1 Ct−2 /Rt Rt−1 > 0. Tedy Ct je monotonní a omezená posloupnost, která konverguje k C, a Rt = Ct−1 + W implikuje, že Rt konverguje k R = C + W. 15
350
300
250
200
150
2
4
6
8
10
12
14
Obr. 2.3: Predikce (zeleně) a pozorování (červeně) Dobropyrinu Dále, užitím Ct = Rt V / (Rt + V ) dostaneme C = RV / (R + V ), z čehož obdržíme kvadratickou rovnici C 2 + CW − V W = 0. p 1 + 4/r − 1 /2, z čehož Tato rovnice má jeden kladný kořen C = rV plyne, že jestliže At = Ct /V , potom At konverguje k A = C/V , což je požadovaná funkce proměnné r. Inverzí navíc získáme užitečný vztah r=
2.3
A2 . 1−A
Diskontní faktor a volba W
Definice 2.2. Diskontní faktor δ definujeme jako δ = 1 − A. 16
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2 2
4
6
8
10
12
14
Obr. 2.4: Adaptivní koeficient Z definice konstantního modelu víme, že Rt = Ct−1 + W a v limitě potom R = C + W = C/(1 − A). Odtud dostaneme W = AC/(1 − A) a W je tedy přímo úměrné C. Chyba ωt mezi pozorováními vede k růstu na W = 100A/(1 − A)% počáteční nejistoty C. Jelikož δ = 1 − A, je R = C/δ. Volba δ, které je spojeno s limitním přizpůsobením se datům, pak ovlivňuje také volbu W . Pro A = 0.1 máme δ = 0.9 a W je asi 11% z C. Pro δ = 0.8 je to už 25%. V konstantním modelu dosáhneme limitního chování velmi rychle, je proto vhodné a přirozené uvažovat konstantní poměr růstu nejistoty pro každé t, ne pouze v limitním případě. Tedy, pro diskontní faktor δ, jehož časté hodnoty jsou 0.8, 0.9, volíme ωt ∼ N[0, Wt ], kde Wt = Ct−1 (1 − δ)/δ, a dostáváme Rt = Ct−1 /δ.
2.4
Limitní tvar modelu - ARIMA(0,1,1)
Pravděpodobně nejznámějšími a široce používanými modely časových řad je třída Box-Jenkinsových ARIMA modelů. Proces ARIMA(0,1,1) si zde odvodíme jako limitní případ polynomiálního modelu prvního stupně. 17
Uvažujme následující vztahy mt−1 = Yt − et a mt = mt−1 + At et . Odtud dostaneme Yt − Yt−1 = et − (1 − At )et−1 . Z toho plyne, že lim [Yt − Yt−1 − et + (1 − At )et−1 ] = 0
t→∞
nebo Yt ≈ Yt−1 + et − δet−1 , kde δ = 1 − A, jak jsme definovali výše. Posloupnost náhodných veličin (et |Dt−1 ) má nulovou střední hodnotu, je nezávislá a pro t → ∞ rozptyl Qt konverguje ke Q = V /(1 − A). Tedy, alternativní limitní model časové řady je Yt = Yt−1 + at − δat−1 , kde at jsou nezávislé, at ∼ N[0, Q], což je právě předpis ARIMA(0,1,1) procesu.
18
Kapitola 3 Dynamický regresní model Regresní modelování zahrnuje konstrukci matematických či statistických modelů, které popisují účinky nezávislých, neboli regresních proměnných na časovou řadu Yt . Uvažujeme-li takovou proměnnou (zapsanou jako řadu pozorování Xt ), pak naším prvním úkolem je najít vztah mezi Xt a skutečnou hodnotou veličiny µt pomocí regresní funkce. Regresní modely se používají v mnoha úlohách predikce, interpolace, při odhadech a regulaci. Nejdříve se zaměříme na obecný dynamický regresní model. Mějme n regresních proměnných a označme je X1 , . . . , Xn . Hodnotu i-té proměnné Xi v čase t označíme Xti . Budeme předpokládat, že pokud je některý člen konstantní, je Xti = 1 pro každé t. Definice 3.1. Pro každé t definujeme obecný dynamický regresní model jako pozorovací rovnice: rovnice systému:
Yt = F′t θ t + νt , θ t = θ t−1 + ω t ,
νt ∼ N[0, Vt ], ω t ∼ N[0, Wt ],
kde F′t = (Xt1 , . . . , Xtn ) je regresní vektor v čase t, θ t je regresní vektor parametrů o rozměrech n × 1 a Wt nazýváme maticí rozptylu vývoje θ t . Standardní (statický) regresní model má stejný tvar jako dynamický, jen Wt = 0 pro každé t a θ t = θ je v čase konstantní. V dynamickém regresním modelu předpokládáme, že lineární regrese je pouze lokální a že regresní vektor parametrů splňuje rovnici náhodné procházky. Vývojová chyba ω t popisuje změny v prvcích vektoru parametrů na intervalu [t − 1, t]. Nulová střední hodnota tohoto vektoru vyjadřuje naše očekávání, že θ t je na tomto intervalu konstantní, zatímco matice rozptylu Wt určuje míru, rozsah pohybu v θ t , a tím rozsah časového intervalu, kde je rozumné očekávat lokální konstantnost. A samozřejmě předpokládáme, že posloupnosti chyb νt a ω t jsou vzájemně nezávislé. 19
Pro n = 1 se zahrnutím konstantního členu dostaneme regresní model určený vektory Ft = (1, Xt )′ a θ t = (αt , βt ). Máme tedy Yt = αt + βt Xt + νt , αt = αt−1 + ωt1 , βt = βt−1 + ωt2 ,
νt ∼ N[0, Wt ],
kde ω t = (ωt1 , ωt2 )′ ∼ N[0, Wt ]. Uvažujme nyní speciální případ regresního modelu, kde αt = 0 a θ t = θt = βt pro každé t. Předpokládáme, že přímka procházející počátkem modeluje lokální vztah, ale její sklon se v různých oblastech liší. Pro ilustraci si uveďme následující tři tabulky. Rok t Mléko: Yt (×109 lbs) Krávy: Ft (×106 ) Rok t Mléko: Yt (×109 lbs) Krávy: Ft (×106 )
1970 1971 117.0 118. 6
1972 1973 1974 1975 1976 1977 120.0 115.5 115.6 115.4 120.2 122.7
12.0
11.8
11.7
1978 121.5
1979 123.4
1980 1981 1982 128.5 130.0 135.8
10.8
10.7
10.8
11.4
10.9
11.2
11.1
11.0
11.0
11.0
V první tabulce reprezentuje řada Yt celkovou roční produkci mléka v USA během let 1970 až 1982 a Ft celkový počet krav, které dávají mléko. V tomto případě θt představuje průměrný roční odběr mléka od jedné krávy v čase t. Model není ani tak zaměřen na predikci Yt , jako spíše na stanovení měnící se produktivity. Rok t Yt Ft
1 2 12 11 4 4
3 4 5 9 5 3 3 2 1
6 7 8 9 10 0 -5 -7 -6 -3 -1 -3 -4 -3 -1
11 12 13 14 7 10 13 12 2 3 4 4
Ve druhém případě použijeme data Ft k predikci ročního prodeje určitého produktu. Yt je změna v ročním prodeji tohoto produktu mezi roky t − 1 a t a Ft je změna v indexu průmyslové výroby mezi roky t − 2 a t − 1.
20
Rok 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985
Tržby společnosti Y Čtvrtletí 1 2 3 4 71.2 52.7 44.0 64.5 70.2 52.3 45.2 66.8 72.4 55.1 48.9 64.8 73.3 56.5 50.0 66.8 80.2 58.8 51.1 67.9 73.8 55.9 49.8 66.6 70.0 54.8 48.7 67.7 70.4 52.7 49.1 64.8 70.0 55.3 50.1 65.6 72.7 55.2 51.5 66.2 75.5 58.5
Celkový obrat na trhu F Čtvrtletí 1 2 3 4 161.7 126.4 105.5 150.7 162.1 124.2 107.2 156.0 165.8 130.8 114.3 152.4 166.7 132.8 115.8 155.6 183.0 138.3 119.1 157.3 169.1 128.6 112.2 149.5 156.9 123.4 108.8 153.3 158.3 119.5 107.7 145.0 155.3 123.1 109.2 144.8 160.6 119.1 109.5 144.8 165.8 127.4
Třetí tabulka zobrazuje čtvrtletní tržby společnosti Yt v porovnání s celkovým prodejem na trhu Ft od roku 1975 až do poloviny roku 1985. Hlavní úkol je stanovit vztah mezi těmito dvěma řadami a určit krátkodobou předpověď. Na obr. 3.1 a 3.2, kde jsou řady vykresleny, jsou patrné sezónní vlivy. Jestliže ale data zakreslíme do jednoduchého grafu, kde závisí na sobě navzájem (obr. 3.3), je tato sezónnost odstraněna a vidíme, že užití regresního modelu je zcela vhodné. 80 180
78 76 74 72 70
160
68 66 64 62 60
140
58 56 54 52 120
50 48 46 44 10
20
30
10
40
20
30
Obr. 3.2: Obrat na trhu
Obr. 3.1: Tržby společnosti
21
40
80 78 76 74 72 70 68 66 64 62 60 58 56 54 52 50 48 46 44 120
140
160
180
Obr. 3.3: Závislost Y a F Definice 3.2. Pro každé t > 0 definujeme dynamický regresní model následovně: pozorovací rovnice: rovnice systému: počáteční informace:
Yt = Ft θt + νt , νt ∼ N[0, Vt ], θt = θt−1 + ωt , ωt ∼ N[0, Wt ], (θ0 |D0 ) ∼ N[m0 , C0 ]
pro nějakou střední hodnotu m0 a rozptyly C0 , Vt , Wt . Věta 3.1. Jednokroková predikce a posteriorní rozdělení jsou pro každé t > 0 dány následovně: a) posterior pro θt−1 : (θt−1 |Dt−1 ) ∼ N[mt−1 , Ct−1 ] pro nějakou střední hodnotu mt−1 a rozptyl Ct−1 ; b) prior pro θt : (θt |Dt−1 ) ∼ N[mt−1 , Rt ], kde Rt = Ct−1 + Wt ; 22
c) jednokroková predikce Yt : (Yt |Dt−1 ) ∼ N[ft , Qt ], kde ft = Ft mt−1 a Qt = Ft2 Rt + Vt ; d) posterior pro θt : (θt |Dt ) ∼ N[mt , Ct ], kde mt = mt−1 +At et a Ct = Rt Vt /Qt , kde At = Ft Rt /Qt a et = Yt −ft . Důkaz. Důkaz je částečně podobný důkazu věty 2.1. Předpokládejme, že věta platí pro t = 1. Dále předpokládejme, že platí a) v čase t. Jestliže zapíšeme Yt a θt jako lineární kombinace θt−1 , νt a ωt , pak okamžitě dostáváme, že Yt a θt mají dvourozměrné normální rozdělení (podmíněné množinou Dt−1 ) se střední hodnotou a rozptylem uvedenými výše v bodě b) a c). Pro kovarianci platí C[Yt , θt |Dt−1 ] = C[Ft θt + νt , θt |Dt−1 ] = C[Ft θt , θt |Dt−1 ] + C[θt , νt |Dt−1 ] = Ft V [θt |Dt−1 ] + 0 = Ft Rt , tedy Yt Dt−1 Qt Ft Rt mt−1 . , ∼N θt Dt−1 Ft Rt Rt ft
Odtud regresní koeficient θt na Yt je právě At = Ft Rt /Qt a (θt |Yt , Dt−1 ) ∼ N[mt , Ct ], kde mt = mt−1 + At (Yt − ft ) , a
Ct = Rt − (Rt Ft )2 /Qt .
Poslední rovnice dává Ct = Rt Vt /Qt a dostáváme d). Příklad 2. Použijte dynamický regresní model na data uvedená v prvních dvou tabulkách. Model jsme vytvořili v programu Maple.
23
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
with(plots): straightline:=proc(S::list,T::list,N,D::positive,V::positive, W::positive) local a,c,e,f,g,h,i,j,k,l,m,n,A,C,M,Q,R,plot1,plot2,plot3; c:=D; M:=N; A:=[]; C:=[]; m:=[]; for n from 1 to nops(S) do R:=c+W; Q:=(S[n])^2*R+V; f:=M*S[n]; e:=T[n]-f; a:=S[n]*R/Q; A:=[op(A),a]; c:=R*V/Q; C:=[op(C),c]; M:=M+a*e; m:=[op(m),M]; od; g:=[seq([j,m[j]],j=1..nops(m))]; h:=[seq([k,C[k]],k=1..nops(C))]; i:=[seq([l,A[l]],l=1..nops(A))]; plot1:=plot(g,color=red,title=‘mean‘); plot2:=plot(h,color=green,title=‘covariance‘); plot3:=plot(i,color=blue,title=‘adaptive coefficient‘); print(‘m=‘,m); print(plot1); print(‘C=‘,C); print(plot2); print(‘A=‘,A); print(plot3); end:
1. Produkce mléka v USA. Pro tento model jsme použili rozptyl V = 1 a W = 0.05. θt můžeme brát 24
jako množství mléka vyprodukovaného jednou krávou ročně. Jako počáteční informace uvažujeme m0 = 10 a C0 = 100. Takto velkým rozptylem dáváme najevo naši nejistotu ohledně θ0 . >straightline([12,11.8,11.7,11.4,11.2,11.1,11.0,11.0,10.8,10.7, 10.8,10.9,11.0],[117,118.6,120,115.5,115.6,115.4,120.2,122.7, 121.5,123.4,128.5,130,135.8],10,100,1,0.05); mean
adaptive coefficient
covariance 0.0076
12 0.082
0.0074
11.5 0.0072
0.08 0.007
11
0.0068
0.078
10.5
0.0066
10
0.076 0.0064
2
4
6
8
10
12
2
4
6
8
10
12
2
4
6
8
10
12
Řada mt je rostoucí kromě t = 4, to odráží rostoucí efektivitu produkce mléka. Hodnoty Ct a At se postupně liší pouze v malých odchylkách, to je dáno tím, že Ft se pohybuje pouze v malém intervalu v mezích od 10.7 do 12.0. Tento jednoduchý model ilustruje, jak dynamičnost θt dává racionální krátkodobou předpověď. Například v časech t = 11 a t = 12 je jednokroková predikce (Y12 |D11 ) ∼ N[129.2, 7.8], (Y13 |D12 ) ∼ N[131.1, 7.9]. Aktuální pozorování jsou Y12 = 130.0 a Y13 = 135.8, naše predikce jsou tedy vyhovující. 2. Roční prodej společnosti. V tomto modelu jsme počítali s V = 1 a W = 0.01 a s počátečními hodnotami m0 = 2 a C0 = 0.81. > straightline([4,4,3,2,1,-1,-3,-4,-3,-1,2,3,4,4],[12,11,9,5,3, 0,-5,-7,-6,-3,7,10,13,12],2,0.81,1,0.01);
25
mean covariance
adaptive coefficient
0.058 0.056
2.8
2.6
2.4
0.22
0.054
0.2
0.052
0.18
0.05
0.16
0.048
0.14
0.046
0.12
0.044
0.1
0.042
0.08
0.04
0.06
0.038
0.04
0.036
0.02
0.034
2.2
0
0.032
–0.02
0.03
–0.04
0.028
–0.06
0.026
–0.08
0.024
2
4
6
8
10
12
14
2
4
6
8
10
–0.1
2
4
6
8
10
12
14
–0.12
Kovariance ani adaptivní koeficient nejsou monotonní, ani se nepřibližují žádné konkrétní hodnotě. Ft totiž nabývá kladných i záporných hodnot a At přebírá znaménko stejné jako má aktuální hodnota Ft .
26
12
14
Kapitola 4 Dynamické lineární modely (DLMs) Polynomiální model prvního stupně a regresní model, které jsme si popsali v předchozích kapitolách, ilustrují základní pojmy a důležité rysy obecných dynamických lineárních modelů. Takové modely si nyní definujeme, následně se však zaměříme na jeden speciální případ - nejčastější skalární tvar. Definice 4.1. Nechť Yt je vektor pozorování o rozměrech r × 1 v časové řadě Y1 , Y2 , . . . . Obecný dynamický lineární model je určen čtveřicí {F, G, V, W}t = {Ft , Gt , Vt , Wt } pro každé t > 0, kde a) Ft je známá matice n × r, b) Gt je známá matice n × n, c) Vt je známá varianční matice r × r, d) Wt je známá varianční matice n × n. Tato čtveřice definuje model, který je v čase t tvaru: (Yt |θ t ) ∼ N[F′t θ t , Vt ], (θ t |θ t−1 ) ∼ N[Gt θ t−1 , Wt ].
(4.1a) (4.1b)
Rovnice 4.1 jsou podmíněny množinou apriorních informací Dt−1 . Ta zahrnuje hodnoty rozptylů Vt a Wt , minulá pozorování Yt−1 , Yt−2 , . . . i počáteční informační množinu D0 . Dt−1 není v rovnicích explicitně vyjádřena, 27
neměli bychom na ni ale zapomínat. Jiná možnost zápisu rovnic 4.1 je Yt = F′t θ t + ν t , θ t = Gt θ t−1 + ω t ,
ν t ∼ N[0, Vt ], ω t ∼ N[0, Wt ],
(4.2a) (4.2b)
kde posloupnosti chyb ν t a ω t jsou nezávislé a vzájemně nezávislé. Rovnice 4.2a je rovnice pozorování modelu, která určuje rozdělení Yt podmíněné veličinou θ t . Rovnice dává do vztahu pozorování a θ t pomocí dynamické lineární regrese s chybou ν t , která má normální rozdělení, a pozorovací varianční matici Vt . Matice Ft představuje regresní matici známých hodnot nezávislých proměnných a θ t je dynamický vektor regresních parametrů, tzv. stavový vektor modelu. Střední hodnota µt = F′t θ t je očekávaná hodnota Yt , která určuje hodnotu řady v čase t. O ν t mluvíme jako o pozorovací chybě. Rovnice 4.2b se nazývá vývojová rovnice systému a vyjadřuje časový vývoj stavového vektoru. Deterministickou složkou vývoje je přechod od stavu θ t−1 ke Gt θ t−1 . Jedná se o jednoduchou lineární transformaci pomocí transformační matice Gt . Na závěr přidáme vývojovou chybu ω t s vývojovou varianční maticí Wt . Poznámka. Jestliže jsou matice Ft a Gt konstantní v čase t, stejně jako pozorovací a vývojová varianční matice Vt a Wt , pak model nazýváme konstantním dynamickým lineárním modelem. Konstantní model je tedy určen na čase nezávislou (statickou) čtveřicí {F, G, V, W}. Tato množina DLMs zahrnuje všechny klasické lineární modely časových řad. Obecný skalární DLM je definován stejně jako v definici 4.1 pro r = 1 a je tedy charakterizován čtveřicí {Ft , Gt , Vt , Wt }, odkud dostáváme (Yt |θ t ) ∼ N[F′t θ t , Vt ], (θ t |θ t−1 ) ∼ N[Gt θ t−1 , Wt ]. Tyto rovnice spolu s počátečními informacemi dávají následující definici. Definice 4.2. Pro každé t > 0 definujeme skalární dynamický lineární model jako pozorovací rovnice: rovnice systému: počáteční informace:
νt ∼ N[0, Vt ], Yt = F′t θ t + νt , θ t = Gt θ t−1 + ω t ω t ∼ N[0, Wt ], (θ 0 |D0 ) ∼ N[m0 , C0 ], 28
pro nějaké apriorní hodnoty m0 a C0 . Pozorovací a vývojová posloupnost chyb jsou nezávislé, vzájemně nezávislé a nezávislé na (θ 0 |D0 ). Věta 4.1. Pro skalární DLM jsou pro každé t > 0 jednokroková předpověď a posteriorní rozdělení následující: a) posterior v čase t − 1: (θ t−1 |Dt−1 ) ∼ N[mt−1 , Ct−1 ] pro nějakou střední hodnotu mt−1 a varianční matici Ct−1 ; b) prior v čase t: (θ t |Dt−1 ) ∼ N[at , Rt], kde at = Gt mt−1 a Rt = Gt Ct−1 G′t + Wt ; c) jednokroková predikce: (Yt |Dt−1 ) ∼ N[ft , Qt ], kde ft = F′t at a Qt = F′t Rt Ft + Vt ; d) posterior v čase t: (θt |Dt ) ∼ N[mt , Ct ], kde mt = at + At et a Ct = Rt − At A′t Qt , kde At = Rt Ft Q−1 a t et = Yt − ft . Důkaz. Důkaz provedeme opět indukcí. Pro t = 1 máme (θ 0 |D0 ) ∼ N[m0 , C0 ]. Podle definice skalárního DLM je θ 1 = G1 θ 0 + ω 1 , kde ω 1 ∼ N[0, W1 ], a tedy (θ 1 |D0 ) ∼ N[G1 m0 , G1 C0 G′1 + W1 ]. Označíme a1 = G1 m0 a R1 = G1 C0 G′1 + W1 . Dále, z pozorovací rovnice dostaneme Y1 = F′1 θ 1 +ν1 , kde ν1 ∼ N[0, V1 ]. Čili pro jednokrokovou predikci platí (Y1 |D0 ) ∼ N[F′1 a1 , F′1 R1 F1 + V1 ] a opět zavedeme značení f1 = F′1 a1 a Qt = F′1 R1 F1 + V1 . lení pro θ 1 si odvodíme obecně dále v důkazu. Předpokládejme, že platí a). Lineární vývojová rovnice s lením v a) přímo dává prior v b). Jelikož je θ t lineární které jsou nezávislé a mají normální rozdělení, je také θ t 29
Posteriorní rozděnormálním rozděfunkcí θ t−1 a ω t , normální. Střední
hodnota at a rozptyl Rt jsou pak odvozeny triviálně. Pozorovací rovnice dává podmíněnou normální hustotu p (Yt |θ t ) = p (Yt |θt , Dt−1 ) . Odtud, užitím b) a vývojové rovnice plyne, že Yt společně s θ t mají normální rozdělení podmíněné Dt−1 s kovariancí C[Yt , θ t |Dt−1 ] = C[F′t θ t + νt , θ t |Dt−1 ] = F′t V [θ t |Dt−1 ] + 0′ = F′t Rt = A′t Qt , kde At = Rt Ft /Qt . Navíc střední hodnotu ft a rozptyl Qt dostaneme jednoduše z pozorovací rovnice a z b) a můžeme určit rozdělení pro jednokrokovou predikci jako Yt Dt−1 ft Qt Ft A′t . ∼N , At Qt Rt θ t Dt−1 at Podmíněné rozdělení θ t dáno Dt = {Yt , Dt−1 } lze získat přímo užitím standardní teorie normality (viz [4]). Regresní vektor θ t na Yt je právě At a dostáváme posterior pro θ t . V důkazu jsme odvodili následující identity. −1 a) At = Rt Ft Q−1 t = Ct Ft Vt ;
b) Ct = Rt − At A′t Qt = Rt (I − Ft A′t ); −1 ′ −1 c) C−1 t = Rt + Ft Ft Vt ;
d) Qt = (1 − F′t At )−1 Vt ; e) F′t At = 1 − Vt Q−1 t . et je chyba v jednom kroku a At je vektor adaptivního koeficientu v čase t. Definice 4.3. Předpovědní funkci ft (k) definujeme pro každé celé k ≥ 0 a čas t jako ft (k) = E[µt+k |Dt ] = E[F′t+k θ t+k |Dt ], kde µt+k = F′t+k θ t+k je funkce středních hodnot. 30
Pro k > 0 dává předpovědní funkce očekávanou hodnotu budoucího pozorování podmíněné současnými informacemi, ft (k) = E[Yt+k |Dt ],
k ≥ 1.
Definici raději uvádíme v očekávaných (středních) hodnotách µt+k namísto Yt+k = µt+k + νt+k . Definice tak zahrnuje i případ k = 0, který dává ft (0) = E[µt |Dt ], což je posteriorní bodový odhad aktuální hodnoty řady.
31
Kapitola 5 Referenční analýza DLMs Určit správné apriorní rozdělení v čase t = 0 k ”nastartování” modelu je možné pouze tehdy, máme-li relevantní, dostatečně vypovídající apriorní informace. To je sice pro Bayesovské předpovídání typické, mohou ale nastat situace, kdy máme prior nejasný či neinformativní. Pro takový případ uvedeme princip referenční analýzy, která model nastartuje, aniž by matematik musel určit informativní počáteční apriorní rozdělení. Obecně, Wt určíme pomocí diskontních faktorů. Tento postup však nemůžeme aplikovat pro t < n (n je počet parametrů, dimenze modelu), jelikož posteriorní kovariance zatím neexistují. Je tedy přirozené uvažovat počáteční kovariance nulové, Wt = 0 pro t = 1, 2, . . . , n + 1. V referenční analýze s n + 1 parametry (jako parametr uvažujeme také V ) potřebujeme alespoň n + 1 pozorování, abychom získali plně specifikované posteriorní rozdělení jedno pozorování pro každý parametr. Je nemožné odhadnout změny v parametrech během prvních n + 1 pozorování, kdy provádíme referenční analýzu. V čase t = n + 1 jsou posteriory plně určeny a v této chvíli se opět vracíme k dynamickému modelu s nenulovou varianční maticí. Referenční analýza DLMs je založena na referenčním neinformativním apriorním předpisu v čase t = 1 definovaným jako p(θ 1 |D0 ) ∝ konstanta. Věta 5.1. Nechť Gt je regulární matice a Wt = 0. Potom apriorní a posteriorní rozdělení stavového vektoru θ t pro t < n jsou tvaru 1 p(θ t |Dt−1 ) ∝ exp − (θ ′t Ht θ t − 2θ ′t ht ), 2 1 p(θ t |Dt ) ∝ exp − (θ ′t Kt θ t − 2θ ′t kt ), 2
32
kde −1 Ht = G′−1 t Kt−1 Gt ht = G′−1 t kt−1
a Kt = Ht + Ft F′t /Vt kt = ht + Ft Yt /Vt s počátečními hodnotami H1 = 0 a h1 = 0. Pro t ≥ n je posteriorní rozdělení (θ t |Dt ) definováno jako v kapitole 4 a platí Ct = K−1 t
a
mt = K−1 t kt .
Důkaz lze nalézt v [4]. Rovnice odvozené ve větě 5.1 jsou pro praktické použití klíčové. Předpoklad regularity Gt je pro výsledek rozhodující, avšak v modelech, které uvažujeme, je vždy splněn.
33
Kapitola 6 Sezónní modely U mnoha časových řad spojených s ekonomickými, obchodními, fyzikálními a biologickými systémy můžeme pozorovat cyklické (periodické) chování. Takové chování je v podstatě nekontrolovatelné a musíme jej akceptovat. Je důležité vyvolanou sezónnost rozpoznat a zahrnout ji jako jeden z faktorů do předpovědního modelu. Asi nejvýznamnějším cyklem je cyklus roční, mezi další důležité patří cykly den/noc, měsíční či čtvrtletní. Existují dvě různé třídy sezónních modelů, dva rozdílné přístupy. První se zabývá modely, u kterých nijak konkrétně tvar sezónní složky nespecifikujeme. Druhý přístup popisuje sezónní faktory pomocí trigonometrických funkcí. V praxi se s úspěchem využívá obou těchto přístupů.
6.1
Fourierovo vyjádření sezónnosti
Některé přístupy vyjádření sezónních vlivů používají lineární kombinace periodických funkcí. Speciální přístup, kterému se v této části budeme věnovat, užívá funkce trigonometrické. Sinové a kosinové vlny jsou totiž velmi jednoduché a výborně aproximují cyklické chování procesů. Sezónní vlivy můžeme pozorovat u elektrických a elektronických systémů, geofyzikálních jevů, včetně např. zemětřesení, a u mnoha procesů, které jsou důsledkem obíhání Země kolem Slunce. Jestliže sezónnost vyjadřujeme pomocí trigonometrických funkcí, mluvíme o Fourierově vyjádření sezónnosti. Poznámka. V této souvislosti je dobré připomenout základy teorie diskrétní Fourierovy tranasformace, kde funkce vyjadřujeme v komplexním tvaru. Nechť G je grupa celých čísel modulo N, G = {0, 1, . . . , N −1}. Nechť L2 (G) : {f : G → C} je konečně rozměrný vektorový prostor se skalárním součinem 1 X hf, gi = f (x)g(x) N x∈G 34
pro f, g ∈ L2 (G). Pro m ∈ G dále definujeme em (j) = (e
2πij N
)m .
Snadno lze ukázat, že em , m = 0, 1, . . . , N − 1, tvoří v L2 (G) ortonormální bázi, tj. ( 1, m = n, hem , en i = 0, m 6= n. Pro f ∈ L2 (G) definujeme její diskrétní Fourierovu transformaci (DFT) předpisem N −1 2πimj 1 X ˆ f (m) = hf, em i = f (j)e− N . N j=0
Každou funkci f ∈ L2 (G) lze tedy zapsat v komplexním tvaru pomocí exponenciální funkce N −1 X ˆ f (j) = f(m)e m. m=0
My vyjádříme periodické funkce ve tvaru reálném, pomocí funkcí sinus/kosinus, což v terminologii DFT nazýváme reálnou diskrétní Fourierovou transformací. Nechť g(t) je reálná funkce definovaná na množině nezáporných celých čísel t = 0, 1, . . . , kde t je časový index. Definice 6.1. Funkci g nazýváme cyklickou (periodickou) funkcí, jestliže pro nějaké celé číslo p ≥ 1 je g(t + np) = g(t) pro všechna celá čísla t ≥ 0 a n ≥ 0. Takové nejmenší p nazýváme periodou funkce g a interval (t, t + p − 1), kde t > 0, nazýváme cyklem. Definice 6.2. p hodnot vyjádřených v jednom cyklu ψj = g(j), kde j = 0, . . . , p − 1, nazýváme sezónními faktory. Uvažujme cyklickou funkci g definovanou pomocí sezónních faktorů ψ0 , . . . , ψp−1 . Tyto faktory můžeme zapsat jako lineární kombinace trigonometrických funkcí. Takové vyjádření závisí na periodě p. Označíme α = 2π/p.
35
Jestliže p = 2q − 1 pro nějaké celé číslo q > 0, pak existuje p reálných hodnot a0 , a1 , . . . , aq−1 , b1 , . . . , bq−1 takových, že pro j = 0, 1, . . . , p − 1 je ψj = a0 +
q−1 X
[ar cos(αrj) + br sin(αrj)] .
r=1
Jestliže p = 2q pro nějaké celé číslo q, pak existuje p reálných hodnot a0 , a1 , . . . , aq , b1 , . . . , bq−1 takových, že pro j = 0, 1, . . . , p − 1 je ψj = a0 +
q−1 X
[ar cos(αrj) + br sin(αrj)] + aq cos(πj).
r=1
Jestliže navíc definujeme bq = 0, můžeme pro libovolné p psát ψj = a0 +
q X
[ar cos(αrj) + br sin(αrj)] ,
r=1
kde j = 0, 1, . . . , p − 1, přičemž pro p liché je aq = 0. Z periodičnosti sinu a kosinu víme, že pro n ∈ Z a x ∈ R platí cos(x + 2πn) = cos x
a
sin(x + 2πn) = sin x.
Dále, pro každé r = 1, . . . , q je p−1 X j=0
cos(αrj) =
p−1 X
sin(αrj) = 0.
j=0
Pro h ∈ Z a k ∈ Z máme p−1 X
cos(αhj) sin(αkj) = 0;
j=0
0, X cos(αhj) cos(αkj) = p, p j=0 , 2 p−1 0, X sin(αhj) sin(αkj) = 0, p j=0 , 2 p−1
36
h 6= k, h = k = 2p , h = k 6= p2 ; h 6= k, h = k = p2 , h = k 6= 2p .
Koeficienty ar a br získáme vynásobením ψj příslušejícím sinem/kosinem a sumací přes j = 0, 1, . . . , p − 1. S použitím výše uvedených identit pak dostáváme a0 = p−1
p−1 X
ψj ,
j=0
p−1
aq = p−1
X
ψj cos(πj)
j=0
a pro r = 1, . . . , q − 1 ar = 2p−1
p−1 X
ψj cos(αrj),
j=0
p−1
br = 2p−1
X
ψj sin(αrj).
j=0
Všimněme si, že sezónní vlivy můžeme získat jako φj = ψj − a0 pro j = 0, 1, . . . , p − 1. Jestliže je součet faktorů roven nule, pak také a0 = 0. Veličiny ar a br se nazývají Fourierovy koeficienty Fourierova vyjádření sezónních faktorů. Funkci Srj definujeme pro j = 0, 1, . . . , p − 1 jako Srj = ar cos(αrj) + br sin(αrj) pro r = 1, . . . , q − 1 a ( aq cos(πj) = (−1)j , p = 2q, Sqj = 0, p = 2q − 1. Funkce Srj se nazývá r-tá harmonická. Lze ji také zapsat jako Srj = Ar cos(αrj + γr ), p kde Ar = a2r + b2r a γr = arctan(br /ar ). Harmonickou můžeme tedy popsat buď pomocí uspořádané dvojice (ar , br ), nebo pomocí dvojice (Ar , γr ). Ar se nazývá amplituda r-té harmonické. Je to maximální hodnota, které Srj nabývá pro j = 0, 1, . . . , p − 1. Fáze γr určuje umístění tohoto maxima. Jestliže p = 2q, pak q-tá harmonická je určena pomocí aq a bq = 0, tedy Aq = |aq | a γq = 0. Takovou funkci nazýváme Nyquistovou harmonickou. Frekvence r-té harmonické je definována jako αr = 2πr/p, kde απ je Nyquistova frekvence, a délka cyklu harmonické je p/r. Ještě poznamenejme, že 37
první harmonická se nazývá základní a má základní frekvenci α a základní délku cyklu p. Ve Fourierově vyjádření je funkce g součtem harmonických, a to tak, že pro každé t > 0 nebo j = p|t (j je zbytek po dělení čísla t číslem p) máme g(t) = ψj = a0 +
q X
Srj ,
r=1
( 0, p = 2q − 1, s posledním členem Sqj = j aq (−1) , p = 2q. Všimněme si, že pro j = p|t je Srt = Srj , r = 1, . . . , q. Funkci g pak můžeme psát v deterministickém tvaru DLM. Definice 6.3. Pro r = 1, . . . , q − 1 definujeme regulární matici typu 2 × 2 Gr a (2 × 1)-rozměrný vektor ar jako cos(ωr) sin(ωr) Gr = − sin(ωr) cos(ωr) a
ar = (ar , br )′ , kde ω = 2π/p. Lze ukázat, že každá Gr je p-periodická. Zvláště, pro každé celé číslo t > 0 a j = p|t máme cos(ωrj) sin(ωrj) cos(ωrt) sin(ωrt) t = Gjr . = Gr = − sin(ωrj) cos(ωrj) − sin(ωrt) cos(ωrt) Odtud dostáváme pro každé j ≥ 0 ar cos(ωrj) + br sin(ωrj) j Gr ar = −ar sin(ωrj) + br cos(ωrj)
pro r = 1, . . . , q −1. První složka vektoru je právě r-tá harmonická Srj funkce g v čase t = j. Užitím vektoru E2 = (1, 0)′ tedy můžeme psát Srj = E′2 Gjr ar pro r = 1, . . . , q − 1. Z tohoto vyjádření dostáváme, že funkce Srt proměnné t jsou právě předpovědní funkce DLMs tvaru {E2 , Gr , ., .} pro r < q. Jelikož je funkce g vyjádřena jako součet harmonických, můžeme ji získat jako předpovědní funkci DLM vytvořeného superpozicí jednotlivých DLMs. 38
Příklad 3. Použitím sezónního modelu proveďte na následujících datech jednokrokovou předpověď. Data představují čtvrtletní (p = 4) prodej kuřat v tisících starých jeden den v líhni ve městě Eiru během let 1974 až 1982. rok 1974 1975 1976 1977 1978 1979 1980 1981 1982
1 131.7 80.4 203.3 177.0 192.2 196.0 352.5 305.9 384.0
čtvrtletí 2 3 322.6 285.6 285.1 347.8 375.9 415.9 438.3 463.2 442.8 509.6 478.6 688.6 508.1 701.5 422.2 771.0 472.0 852.0
4 105.7 68.9 65.8 136.0 201.2 259.8 325.6 329.3
Model je opět vytvořen v programu Maple. > > > > > > > > > > > > > > > > > > > > > > > >
with(linalg): season:=proc(S::list,p) local A,C,F,G,H,K,Q,R,T,V,W,a,e,f,h,i,k,m,n,q,s,t; V:=1; W:=matrix(2,2,0); T:=[]; F:=matrix(2,1,[1,0]); G:=evalf(matrix(2,2,[cos(2*Pi/k),sin(2*Pi/k),-sin(2*Pi/k),cos (2*Pi/k)])); H:=matrix(2,2,0); h:=matrix(2,1,0); s:=1/nops(S)*sum(S[i],i=1..nops(S)); for n from 1 to nops(S) do t:=S[n]-s; T:=[op(T),t]; od; for n from 1 to 2 do K:=evalm(H + F &* transpose(F) /V); k:=evalm(h + F *T[n]/V); H:=evalm(transpose(inverse(G)) &* K &* inverse(G)); h:=evalm(transpose(inverse(G)) &* k); od; for n from 3 to nops(T) do a:=evalm(G &* m); 39
> > > > > > > > > > > > > >
R:=evalm(G &* C &* transpose(G)+W); f:=evalm(transpose(F) &* a); f:=f[1,1]; q:=evalm(transpose(F) &* R &* F); q:=q[1,1]; Q:=q+V; A:=evalm(R &* F *1/Q); e:=T[n]-f; m:=evalm(a+A*e); C:=evalm(R - A &* transpose(A)*Q); od; m[1,1]+s; end:
Jednokroková předpověď je 518.7.
6.2
Modely plné sezónnosti ve Fourierově vyjádření
Jak jsme popsali v předchozí části, jakékoliv p-sezónní faktory můžeme vyjádřit pomocí harmonických složek. Nyní si specifikujeme DLMs, které mají stejné periodické předpovědní funkce (vyjádřené Fourierovými koeficienty) zkombinované principem superpozice. Definice 6.4. Mějme dán (p − 1)-rozměrný vektor F a (p − 1) × (p − 1) matici G jako ( (E′2 , . . . , E′2 ), pro p = 2q − 1, F′ = (E′2 , . . . , E′2 , 1), pro p = 2q, a
( diag[G1 , . . . , Gq−1 ], pro p = 2q − 1, G= diag[G1 , . . . , Gq−1 , −1], pro p = 2q.
Potom model {F, G, ., .} definujeme jako model plné sezónnosti. Takový model je vytvořen superpozicí harmonických složek DLMs {E2 , Gr , ., .}
40
a pro předpovědní funkci a střední hodnotu platí následující. Předpovědní funkce je tvaru ft (k) =
q X r=1
Srk =
q−1 X
[atr cos(ωrk) + btr sin(ωrk)] + atq (−1)k ,
r=1
kde atq = 0 pro každé t, jestliže p = 2q − 1, a posteriorní střední hodnota stavového vektoru v čase t je ( (at1 , bt1 ; . . . ; at,q−1 , bt,q−1 )′ , pro p = 2q − 1, E[θ t |Dt ] = mt = ′ (at1 , bt1 ; . . . ; at,q−1 , bt,q−1 ; atq ) , pro p = 2q. θ t je vektor Fourierových koeficientů a a, b jsou hodnoty jejich odhadů.
41
Kapitola 7 Rozšíření standardních DLMs Standardní DLMs představují základní postupy pro analýzu časových řad. Při modelování řady je velmi užitečné pracovat s trendovou a sezónní složkou. Speciálně je tento přístup vhodný při retrospektivní analýze a následné krátkodobé predikci. Nejprve se v krátkosti zmíníme o polynomiálních modelech druhého stupně, abychom si později mohli definovat model, který kombinuje trendovou a sezónní složku. Pokud by čtenáře polynomiální modely zajímaly více, nechť nahlédne do [4].
7.1
Polynomiální model druhého stupně
Polynomiální modely nacházejí široké využití pro časové řady a pro predikci v mnoha odvětvích aplikované statistiky. Poskytují jednoduché flexibilní formy k popisu lokální trendové složky (trendem zde máme na mysli hladkou změnu v čase). Takové trendy jsou často dobře aproximovány polynomiálními funkcemi nízkého řádu. Polynomiální model prvního nebo druhého stupně je adekvátním modelem vhodným pro krátkodobou předpověď. Navíc má širokou škálu použitelnosti, pokud ho superpozicí zkombinujeme se sezónní či regresní složkou. Polynomiální DLMs je třída modelů definovaných v kapitole 4, jejichž regresní vektor F a matice systému G jsou v čase konstantní. Předpokládejme, že n ∈ N. Definice 7.1. DLM s předpovědní funkcí tvaru ft (k) = at0 + at1 k + · · · + at,n−1 k n−1 , kde t ≥ 0 a k ≥ 0, nazýváme polynomiálním modelem n-tého stupně. 42
Polynomiální model druhého stupně má předpovědní funkci ve tvaru ft (k) = at0 + at1 k,
(7.1)
kde t ≥ 0 a k ≥ 0, a zapíšeme ho jako 1 1 1 , , ., . . 0 0 1 θt1 µt Jestliže uvažujeme parametrizaci θ t = = , pak můžeme model θt2 βt zapsat pomocí obvyklých rovnic, a to Yt = µt + νt , µt = µt−1 + βt−1 + ωt1 , βt = βt−1 + ωt2 , kde νt ∼ N[0, Vt ] a ω t = (ωt1 , ωt2 )′ ∼ N[0, Wt ] a u kterých předpokládáme ”obvyklé” nezávislosti. µt je aktuální hodnota řady, βt−1 představuje změnu v této hodnotě. Teorie standardních DLMs nám dává informaci určující prior v čase t − 1 jako (θ t−1 |Dt−1 ) ∼ N[mt−1 , Ct−1 ], kde mt−1 a Ct−1
mt−1 = bt−1
Ct−1,1 Ct−1,3 = . Ct−1,3 Ct−1,4
Odtud vidíme, že v definici 7.1 je at0 = mt−1 a at1 = bt−1 .
7.2
Trendově-sezónní model
Nejzákladnější a pravděpodobně nejvhodnější deskriptivní model je takový, který v sobě shrnuje polynomiální trend a sezónní složku. V časové řadě předpokládáme čtvrtletní periodu a zároveň postupný nárůst trendu. Tento model je velmi často používán v praxi jako první krok v retrospektivní analýze časových řad. Jakýkoliv pohyb v trendu připisujeme změnám v hodnotách a růstu parametrů. První složkou modelu je lineární trend. Budeme uvažovat polynomiální model druhého stupně 1 1 1 , ., . . , 0 1 0 43
K tomuto přidáme sezónní složku, kterou můžememe vyjádřit v podobě 3 Fourierových koeficientů. Fourierovo vyjádření má právě dvě harmonické složky: první představuje základní roční cyklus, druhá přidává Nyquistovu frekvenci. Fourierův model je podle definice 6.4 tvaru ( ) 1 c s 0 0 , −s c 0 , ., . , 1 0 0 −1
kde pro ω = π/2 je c = cos(ω) = 0 a s = sin(ω) = 1. Pro 5-dimenzionální trendově-sezónní model je tedy F = (1, 0, 1, 0, 1)′ a
1 0 G= 0 0 0
1 0 0 0 1 0 0 0 0 0 1 0 0 −1 0 0 0 0 0 −1
Odpovídající vektor parametrů obsahuje 5 prvků, θ t = (θt1 , . . . , θt5 )′ , které mají následující význam. ⋄ θt1 je základní hodnota řady v čase t; ⋄ θt2 je nárůst mezi časy t − 1 a t; ⋄ θt3 a θt4 jsou Fourierovy koeficienty první harmonické v čase t, S1t = θt3 ; ⋄ podobně S2t = θt5 je koeficient Nyquistovy harmonické. Předpokládáme, že pozorovací rozptyl V je v čase konstantní a známý. Vývojové varianční matice Wt se v tomto případě v čase mění. Jelikož se tento model skládá ze dvou jednotlivých modelů, je vhodné uvažovat také dva diskontní faktory. Označíme δT diskontní faktor pro trendovou složku a δS diskontní faktor pro složku sezónní. Vývojové rovnice modelu jsou následující. Jestliže (θ t−1 |Dt−1 ) ∼ N[mt−1 , Ct−1 ], pak (θ t |Dt−1 ) ∼ N[at , Rt ], kde at = Gmt−1 a Rt = Pt + Wt , kde Pt = GCt−1 G′ 44
a vývojová varianční matice Wt je definována jako Wt = block diag{PtT (δT−1 − 1), PtS (δS−1 − 1)}, kde PtT je levý horní blok matice Pt o rozměrech 2 × 2 a PtS je pravý dolní blok Pt o rozměrech 3 × 3. Hodnoty trendového diskontního faktoru se většinou pohybují v rozmezí od 0.8 do 1.0, menší hodnoty naznačují větší změny. δS je obvykle větší než δT , což znamená, že v každé čtvrtině je o sezónních složkách získáváno menší množství informací než o složkách trendových. Navíc takový vztah vyjadřuje naše očekávání, že sezónní vlivy jsou stabilnější na rozdíl od relativně rychle se měnícího trendu. Nejsme schopni předvídat trvalý pohyb a změny v trendu a sezónnosti, můžeme je pouze identifikovat a následně je odhadnout. Příklad 4. Proveďte jednokrokovou predikci pomocí trendově-sezónního modelu. Data představují čtvrtletní produkci piva v megalitrech od března 1956 do ledna 1973 v Austrálii.
rok 1956 1957 1958 1959 1960 1961 1962 1963 1964
1 308.4 320.4 313.4 314.3 311.4 339.2 345.7 363.4
čtvrtletí 2 3 284.4 212.8 262.0 227.9 271.9 232.8 261.4 226.8 286.1 226.5 294.7 232.6 279.1 249.8 293.8 254.7 313.4 272.8
4 226.9 236.1 237.0 249.9 260.4 257.2 269.8 277.5 300.1
rok 1965 1966 1967 1968 1969 1970 1971 1972 1973
1 369.5 386.1 402.3 404.8 442.3 445.9 466.2 487.0 506.1
čtvrtletí 2 3 330.8 287.8 335.2 288.0 352.8 316.1 393.0 318.9 383.1 331.6 386.6 357.2 409.6 369.8 419.2 376.7
4 305.9 308.3 324.9 327.0 361.4 373.6 378.6 392.8
Model má následující tvar. > > > > > > > > >
with(linalg): trendseason:=proc(S::list,dT,dS) local A,C,F,G,H,K,M,N,Z,P,Q,R,T,V,W,a,e,f,h,i,k,m,n,q,s,t; V:=1; T:=[]; F:=matrix(5,1,[1,0,1,0,1]); G:=matrix(5,5,[1,1,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,-1,0,0,0,0,0, 0,-1]); H:=matrix(5,5,0); h:=matrix(5,1,0);
45
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
s:=1/nops(S)*sum(S[i],i=1..nops(S)); for n from 1 to nops(S) do t:=S[n]-s; T:=[op(T),t]; od; for n from 1 to 5 do K:=evalm(H + F &* transpose(F) /V); k:=evalm(h + F *T[n]/V); H:=evalm(transpose(inverse(G)) &* K &* inverse(G)); h:=evalm(transpose(inverse(G)) &* k); od; m:=evalm(inverse(K) &* k); C:=evalm(inverse(K)); for n from 6 to nops(T) do a:=evalm(G &* m); P:=evalm(G &* C &* transpose(G)); M:=matrix(2,2,[P[1,1],P[1,2],P[2,1],P[2,2]]); N:=matrix(3,3,[P[3,3],P[3,4],P[3,5],P[4,3],P[4,4],P[4,5],P[5,3], P[5,4],P[5,5]]); Z:=matrix(2,3,0); W:=blockmatrix(2,2,[M*(dT^(-1)-1),Z,transpose(Z),N*(dS^(-1)-1)]); R:=evalm(P + W); f:=evalm(transpose(F) &* a); f:=f[1,1]; q:=evalm(transpose(F) &* R &* F); q:=q[1,1]; Q:=q+V; A:=evalm(R &* F *1/Q); e:=T[n]-f; m:=evalm(a+A*e); C:=evalm(R- A &* transpose(A)*Q); od; m[1,1]+s; end:
Zvolili jsme δT = 0.85 a δS = 0.97. Jednokroková předpověď je 431.6. Na obr. 7.1 jsou vykreslena pozorování včetně predikce a jasně můžeme vidět sezónní i trendovou složku. 46
500
450
400
350
300
250
0
10
20
30
40
50
60
70
Obr. 7.1: Čtvrtletní produkce piva v Austrálii Porovnejme nyní model trendově-sezónní a sezónní. Data představují čtvrtletní prodej určitého výrobku a jsou vykreslena na obr. 7.2 a 7.3 společně s jednokrokovými predikcemi. 800
800
700
700
600
600
500
500
400
400
5
10
15
20
25
5
Obr. 7.2: Trendově-sezónní model
10
15
20
25
Obr. 7.3: Sezónní model
Je zřejmé, že model trendově-sezónní je pro tento případ vhodnější. Data se mění se čtvrtletní periodou a navíc postupně rostou. Sezónní model trendovou složku neuvažuje a v předpovědi se drží ”při zemi”.
47
Seznam použité literatury [1] Hyndman R. J. Time Series Data Library, Accessed on 26.4.2008. [2] Jaynes E. Probability Theory: The Logic of Science, Cambridge University Press, 2003. [3] Koop G. Bayesian Econometrics, Wiley, 2003. [4] West M., Harrison J. Bayesian Forecasting and Dynamic Models, Springer-Verlag, 1989. [5] Zellner A. An Introduction to Bayesian Inference in Econometrics, Wiley, 1996. [6] en.wikipedia.org/wiki/Mathematical_model
48
