ZÁPADOČESKÁ UNIVERZITA V PLZNI Fakulta aplikovaných věd Katedra matematiky
DIPLOMOVÁ PRÁCE Předpovědi vysokofrekvenčních časových řad pomocí modelů typu ARCH a metody Monte Carlo
Plzeň 2013
Vypracoval: Bc. Tomáš Pakosta Vedoucí práce: RNDr. Blanka Šedivá, Ph.D.
Čestné prohlášení Prohlašuji, že jsem zadanou diplomovou práci vypracoval samostatně s přispěním vedoucího práce. Uvedl jsem všechny literární prameny a publikace, ze kterých jsem čerpal.
Datum: 20. května 2013
....................... podpis
Poděkování Rád bych na tomto místě poděkoval všem, kteří mi s prací přímo i nepřímo pomohli, protože bez nich by tato práce nevznikla. Především bych chtěl poděkovat vedoucí práce RNDr. Blance Šedivé, Ph.D. za vedení, pomoc a odborné rady při zpracování této práce.
!!! Originál zadání !!!
Abstrakt Hlavním cílem této práce je získat model typu ARCH pro zvolená a popsaná vysokofrekvenční data, konkrétně data časové řady ceny zlata. Dalším cílem je získat za pomoci vybraného a odhadnutého modelu typu ARCH předpověď pro danou časovou řadu zkonstruovanou pomocí metody Monte Carlo. Metoda Monte Carlo je v této práci také popsána a hlavně její aplikace pro modely typu ARCH. Obsahem práce je také zpracovaný popis zvoleného software Project R, který poslouží jako dobrý nástroj pro zpracování časové řady z pohledu modelování a předpovídání. Project R lze také využít k velmi dobré interpretaci výsledků pomocí programovatelných grafů. V závěrečné části této práce je zpracováno porovnání předpovědí se skutečnými hodnotami v odhadovaném období. Klíčová slova: časová řada, volatilita, ARCH, GARCH, Project R, Monte Carlo
Abstract The main goal of this thesis is to obtain a model ARCH type for chosen and described high-frequency data, namely data time series prices of gold. The next target is to obtain forecast of the time series constructed using Monte Carlo methods with the help of selected estimated ARCH type model. Monte Carlo method is also described in this thesis and especially application to ARCH type models. The thesis also include description of the selected software Project R which will serve as a good tool for processing time series from the perspective of modeling and forecasting. Project R would be used for a very good interpretation of the results using programmable graphs. In the final section of this thesis is processed comparing between predictions with the actual values in the estimated period. Keywords: time series, volatility, ARCH, GARCH, Project R, Monte Carlo
Obsah 1 Úvod
1
2 Časové řady typu ARCH 2.1 Finanční časové řady . . . . . . . . . . . . . . . . . . . 2.2 Charakteristické rysy chování finančních časových řad 2.2.1 Předpoklad normality . . . . . . . . . . . . . . 2.2.2 Předpoklad lineární závislosti . . . . . . . . . .
. . . .
2 6 6 8 10
3 Modely volatility 3.1 Základní reprezentace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Lineární modely volatility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Modely ARCH (Autoregressive Conditional Heteroscedasticity) . . . . . . . . . . . . . . .
11 11 12 12
3.3 3.4 3.5 3.6 3.7 3.8
3.2.2 Modely GARCH (Generalized ARCH) Výstavba modelů volatility . . . . . . . . . . Testování podmíněné heteroskedasticity . . . Volba řádu modelu . . . . . . . . . . . . . . . Odhad parametrů . . . . . . . . . . . . . . . . Verifikace modelu . . . . . . . . . . . . . . . . 3.7.1 Kritéria hodnocení modelu . . . . . . Konstrukce předpovědí . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . .
. . . . . . . .
. . . . . . . .
14 15 15 16 17 18 18 19
4 Monte Carlo 21 4.1 Předpovědi modelů GARCH pomocí MC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5 Project R 23 5.1 Počáteční nastavení a načtení dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.2 Použité nástroje a funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6 Aplikační část 6.1 Základní popis a předzpracování dat 6.2 Předpoklad normality . . . . . . . . 6.3 Fáze 1: Odhad úrovňového modelu . 6.3.1 Test residuí . . . . . . . . . . 6.4 Fáze 2: Odhad modelu volatility . . 6.4.1 Test residuí a úprava modelu 6.5 Celkový konečný model . . . . . . . 6.6 Monte Carlo předpovědi . . . . . . . 6.6.1 Porovnání se skutečností . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
28 29 35 37 40 42 44 47 49 53
7 Závěr
58
Literatura
59
A Přílohy obsažené na CD
60
B Přílohy tištěné
61
Seznam obrázků 1 2 3 4 5 6 7 8 9 10 11 12
Vývoj časové řady logaritmů výnosů ceny zlata . . . . . Simulace modelu ARCH(2) . . . . . . . . . . . . . . . . Vývoj časové řady zlata za období 9.2.2005-31.10.2012 . Vývoj časové řady zlata a odhad exponenciálního trendu Graf ACF pro řadu cen zlata . . . . . . . . . . . . . . . Graf PACF pro řadu cen zlata . . . . . . . . . . . . . . Vývoj logaritmů cen zlata a odhad lineárního trendu . . Vývoj časové řady logaritmů výnosů ceny zlata . . . . . Graf ACF logaritmů výnosů časové řady ceny zlata . . . Graf PACF logaritmů výnosů časové řady ceny zlata . . caption . . . . . . . . . . . . . . . . . . . . . . . . . . . Q-Q graf výnosů . . . . . . . . . . . . . . . . . . . . . .
13 14 15
Řada výnosů a odhad modelu ARMA(1,0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Histogram residuí, gaussovský jádrový odhad hustoty rozdělení a teoretické normální rozdělení . 40 ACF residuí ARMA(1,0) modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
16
PACF residuí ARMA(1,0) modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
17
ACF kvadrátů residuí ARMA(1,0) modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
18
PACF kvadrátů residuí ARMA(1,0) modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
19
21
Histogram residuí modelu GARCH(1,1), gaussovský jádrový odhad hustoty rozdělení a teoretické normální rozdělení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Histogram residuí {ˆ } modelu GARCH(1,1), gaussovský jádrový odhad hustoty rozdělení a teoretické Studentovo t-rozdělení . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Podmíněný rozptyl odhadnutého modelu GARCH(1,1) a nepodmíněný rozptyl . . . . . . . . . . 48
22 23 24 25 26 27 28 29 30 31
ˆ t na 20 dní dopředu pomocí metody MC . Simulované předpovědi podmíněného rozptylu h Simulované předpovědi na 20 dní pro residua εˆt modelu ARMA(1,0) pomocí metody MC Simulované hodnoty residuí ARMA modelu a intervaly empirické kvantilové funkce. . . . Residua ARMA modelu a empirické intervaly spolehlivosti metody Monte Carlo. . . . . . Residua ARMA modelu a 95% intervaly spolehlivosti pomocí příkazu predict. . . . . . . . Graf odhadů intervalů metody MC (červeně) a pomocí příkazu predict (modře) . . . . . . Graf odchylek mezi hornímy hranicemi intervalových pásů daných metod . . . . . . . . . Graf odchylek mezi dolními hranicemi intervalových pásů daných metod . . . . . . . . . . Graf historické střední hodnoty logaritmů výnosů časové řady ceny zlata . . . . . . . . . . Graf historického rozptylu hodnoty logaritmů výnosů časové řady ceny zlata . . . . . . . .
20
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
5 13 28 29 30 31 32 33 34 34 35 36
50 50 52 53 53 55 56 56 61 61
Seznam tabulek 1 2 3 4 5
Základní údaje časové řady cen zlata . . . . . . . . Statistické charakteristiky řady výnosů . . . . . . . Výsledky testů normality pro řadu výnosů . . . . . Tabulka kritérií pro různé modely typu ARMA . . Tabulka odhadů parametrů pro ARMA(1,0) model
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
29 35 35 38 38
6 7 8
Výsledky testů normality pro řadu {εt } . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Hodnoty maximalizované věrohodnostní funkce, AIC a BIC modelů GARCH . . . . . . . . . . . 43 Odhadnuté parametry modelu GARCH(2,1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
9
Odhadnuté parametry modelu GARCH(1,1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
10
Testy standardizovaných residuí {ˆ } modelu GARCH(1,1) . . . . . . . . . . . . . . . . . . . . . . 44
11
Odhadnuté parametry modelu GARCH(1,1) s předpokladem Studentova t-rozdělení pro residua {ˆ } . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
12
Testy standardizovaných residuí {ˆ } modelu GARCH(1,1) s předpokladem Studentova t-rozdělení 46
1: Úvod
1
Úvod
S pojmem „časová řadaÿ je svázána spousta teorie, dat, modelů či aplikací. Z tohoto také vyplývá mnoho různých názorů či přístupů, jak se postavit ke zpracování časových řad. Velký zájem, zejména v ekonomii, je kladen na zpracování krátkodobých časových řad, tj. řad, jejichž frekvence hodnot je v obdobích kratších než 1 rok. Zvláštní postavení, v této skupině časových řad, pak mají denní časové řady, u kterých bylo empirickým výzkumem zjištěno, že právě tyto časové řady se vyznačují v čase proměnlivou variabilitou (tzv. volatilita). Tato vlastnost však porušuje předpoklady pro použití klasických lineárních modelů. Počátkem 80. let přišel Robert Engle s novou koncepcí modelování volatility, kde zavádí modely typu ARCH neboli autoregresní modely s podmíněnou heteroskedasticitou (Autoregressive Conditional Heteroscedasticity). Koncepce modelů typu ARCH se postupně rozvíjela v další modely tohoto typu tak, aby co nejvíce reflektovaly různá typická chování určitých druhů vysokofrekvenčních časových řad. Další vývoj modelů typu ARCH pak vychází z myšlenky modifikace předpokladu normality, na kterém je založena většina modelů časových řad. Cílem této diplomové práce je zpracovat finanční časovou řadu ceny zlata za období 9.2.2005 - 31.10.2012 a určit její předpověď. Dosažením tohoto cíle se předpokládá užití modelů podmíněné heteroskedasticity a užití softwaru Project R, pomocí kterého budou získány všechny výsledky a grafické výstupy. Pro určení charakteru předpovědi časové řady cen zlata bude využita metoda Monte Carlo a na konci práce bude popsáno zhodnocení použití metody Monte Carlo pro předpovídání dalšího vývoje časové řady. V první části diplomové práce je popsána potřebná teorie, která je využita v aplikační části. Jedná se o popis finančních časových řad, jejich vlastností, předpokladů a charekateristických rysů chování těchto řad. V druhé části jsou popsány samotné modely volatility použité pro zpracování podmíněné heteroskedasticity v časové řadě. K těmto modelům je také popsána jejich výstavba, volba řádu, odhad parametrů, verifikace modelu a předpovědi. Další částí je pak popis metody Monte Carlo, která je dále využita k odhadu předpovědi časové řady ceny zlata. Před zpracováním časové řady ceny zlata je nezbytný popis programu, ve kterém je aplikační část zpracována a tím je Project R. Samotná aplikační část je poslední částí této práce, ve které je vyšetřen průběh časové řady za pomoci popsané teorie v předcházejících částech. Odhadnutý model časové řady cen zlata je ověřen a použit k předpovědím. Závěrem diplomové práce bude zhodnocení předpovědí pro danou časovou řadu pomocí metody Monte Carlo. Také bude stručně popsáno možné rozšíření této práce.
—1—
2: Časové řady typu ARCH
2
Časové řady typu ARCH
Ekonomické časové řady můžeme popsat jako empirická pozorování v oblasti ekonomiky, která jsou uspořádána v čase, do řady hodnot, směrem od minulosti do přítomonosti. Pokud ekonomické časové řady budeme rozlišovat z pohledu délky intervalu sledovaných hodnot, zjistíme, že spadají do dvou skupin a těmi jsou dlouhodobé časové řady, které mají sledované hodnoty v ročních nebo delších časových intervalech. Naopak krátkodobé časové řady mají interval sledovanosti kratší než jeden rok. S tímto rozdělením souvisí také tvar ekonomických časových řad. Tvar časové řady bude vyhlazenější, čím delší interval sledovaných hodnot budeme mít. To vyplývá z poznatku, že časové řady mají časově svázané jednotlivé hodnoty. Kdybychom časovou návaznost neuvažovali, pojem časové řady by zcela ztratil smysl a neurčili bychom jejich tvar a charakteristické vlastnosti. Charakteristickými vlastnostmi ekonomických časových řad jsou: trend, sezónnost, nelinearita, podmíněná heteroskedasticita a společné vlastnosti více časových řad jako je například tzv. společný trend. Tyto vlastnosti se zpravidla nevyskytují ve všech časových řadách současně. Výskyt některé z vlastností závisí na typu časové řady, např. heteroskedasticita, která zde bude později zkoumána se vyskytuje u vysokofrekvenčních řad (spadají pod krátkodobé časové řady), kterými jsou právě například finanční časové řady. Pro potřeby teoretických modelů lze časovou řadu chápat jako speciální typ náhodného procesu. Stochastický proces Mějme v čase uspořádanou řadu náhodných veličin {X(s, t), s ∈ S, t ∈ T}, kde S je výběrový prostor a T je indexní řada, pak ji nazmeme stochastickým procesem, kde pro každé t ∈ T je X(., t) náhodná veličina definovaná na výběrovém prostoru S a pro každé s ∈ S je X(s, .) realizace stochastického procesu definovaná na indexní řadě T (tj. uspořádaná řada čísel, z níž každé odpovídá jedné hodnotě indexní řady). Budeme předpokládat, že indexní řada je řadou celých čísel, tj. T = {0, ±1, ±2, ±3, . . .} a že uvažované stochastické procesy mají nespojitý charakter. Jelikož časovou řadu lze chápat jako realizaci stochastického procesu, budou se dále stochastické procesy označovat jako {Xt , t = 0, ±1, ±2, ±3, . . .}, resp. {Xt } a jejich realizace, tj. časové řady jako Xt [1]. Časová řada je tedy systémem náhodných veličin {X1 , X2 , ...} (resp. {X1 , X2 , ..., XT }, t = 1, 2, ..., T ) a {x1 , x2 , ...} (resp. {x1 , x2 , ..., xT }, t = 1, 2, ..., T ) jsou jednotlivé realizace časové řady. Časovou řadu Xt nazveme bílým šumem, pokud Xt je posloupnost nekorelovaných stejně rozdělených náhodných veličin s nulovou střední hodnotou a s konečným rozptylem σ 2 . Bílý šum1 označíme WN(0, σ 2 ). Má-li Xt navíc normální rozdělení N(0, σ 2 ), nazývá se Xt Gaussovský bílý šum. Stacionarita Striktně stacionárním procesem nazveme takový stochastický proces, jestliže pro jakoukoliv konečnou indexní část t1 , t2 , ..., tn z T = 0, ±1, ±2, ±3, . . . a jakékoliv reálné číslo k, pro které ti + k ∈ T, i = 1, 2, ..., n, platí: F (Xt1 , Xt2 , ..., Xtn ) = F (Xt1 +k , Xt2 +k , ..., Xtn +k ), kde F (.) je sdružená distribuční funkce. Striktní (silná, kompletní) stacionarita tedy požaduje, aby libovolná skupina náhodných veličin z časové řady měla v případě libovolného časového posunu stejné rozdělení jako neposunutá [1]. 1 White
Noise.
—2—
2: Časové řady typu ARCH
Časovou řadu {Xt : t ∈ T} pak nazveme slabě stacionární, jestliže platí: (i) EXt2 < ∞, (ii) EXt = µ pro všechna t ∈ T a nějaké µ ∈ R, (iii) cov(Xp , Xr ) = γ(|p − r|) pro všechna p, r ∈ T. Bod (iii) znamená, že autokovarianční funkce Xt závisí pouze na vzdálenosti indexů p a r. U stacionárních časových řad proto autokovarianční funkci píšeme jen s jedním argumentem, který udává tuto vzdálenost: ozn.
γ(p, r) = γ(l) pro l = |p − r|. Striktně stacionární proces, který má první dva obecné momenty konečné, je také slabě stacionární proces. Stochastický proces se označuje jako normální nebo gaussovský, jestliže jeho sdružené pravděpodobnostní rozdělení je normální. Striktní a slabá stacionarita jsou pak ekvivalentními pojmy, neboť normální rozdělení je jednoznačně charakterizováno prvními dvěma obecnými momenty. V dalších částech se zabýváme gaussovskými stochastickými procesy, protože většina zde uvedených výsledků z těchto procesů vychází [1]. Autokorelační funkce (ACF) Pro stacionární stochastický proces Xt lze vyjádřit autokovarianční funkci mezi veličinami Xt a Xt−k takto: γk = cov(Xt , Xt−k ) = E[(Xt − µ)(Xt−k − µ)],
k = 0, 1, 2, ..., T − 1.
Autokrelační funkci (ACF) pak jako: ρk = p
cov(Xt , Xt−k ) γk p = , γ0 D(Xt ) D(Xt−k )
kde vzhledem ke stacionaritě procesu var(Xt ) = Xt−k = γ0 . V případě stacionárního stochastického procesu má autokorelační funkce následující vlastnosti: (i) ρ0 = 1, (ii) |γk | ≤ γ0 ; |ρk | ≤ 1 pro všechna k, (iii) (iii) γk = γ−k a ρk = ρ−k pro všechna k. Autokovarianční a autokorelační funkce je tedy symetrická kolem k = 0. Tato vlastnost vyplývá ze skutečnosti, že časová vzdálenost dvojic veličin Xt a Xt+k , Xt−k a Xt je stejná. Z tohoto důvodu a z důvodu (i) je autokorelační funkce vyjadřována pouze pro k > 0. Graf autokorelační funkce se nazývá korelogram. Výběrová autokorelační funkce Ve většině případů jsou parametry µ, γ0 a ρk neznámé. Za předpokladu stacionarity může být střední hodnota procesu µ odhadována prostředictvím výběrového průměru a rozptyl procesu γ0 může být odhadován pomocí výběrového rozptylu: T X ¯= 1 X Xt , T t=1
s2 =
T 1X ¯ 2, (Xt − X) T t=1
kde T je počet hodnot časové řady. Odhad ρk je dán výběrovou autokorelací se zpožděním k: PT ¯ ¯ (Xt − X)(X t−k − X) ρˆk = t=k+1 , k = 1, 2, ..., T − 1. PT ¯ 2 (Xt − X) t=1
—3—
2: Časové řady typu ARCH
Parciální autokorelační funkce (PACF) Korelace mezi veličinami Xt a Xt−k může být způsobena jejich korelací s veličinami Xt−1 , Xt−2 , ..., Xt−k+1 . Parciální autokorelace podávají informaci o korelaci veličin Xt a Xt−k očištěnou o vliv veličin ležících mezi nimi. Parciální korelaci se zpožděním k vyjadřuje parciální regresní koeficient φkk v autoregresi k-tého řádu: Xt = φk1 Xt−1 + φk2 Xt−2 + ... + φkk Xt−k + et , kde veličina et je nekorelovaná s veličinami Xt−j , j ≥ 1. Vynásobíme-li obě strany předchozí rovnice Xt−j , potom střední hodnota této rovnice má tvar: γj = φk1 γj−1 + φk2 γj−2 + ... + φkk γj−k takže platí: ρj = φk1 ρj−1 + φk2 ρj−2 + ... + φkk ρj−k a pokud dosadíme za j = 1, 2, ..., k, získáme systém rovnic. Dále pak postupně dosazujeme za k = 1, 2, ... a použitím Cramerova pravidla získáme φ11 , φ22 , ...φkk . Například φkk vypadá takto: 1 ρ1 ρ2 · · · ρk−2 ρ1 ρ1 1 ρ1 · · · ρk−3 ρ2 .. .. .. .. .. .. . . . . . . ρk−1 ρk−2 ρk−3 · · · ρ1 ρk φkk = 1 ρ1 ρ2 · · · ρk−2 ρk−1 ρ1 1 ρ1 · · · ρk−3 ρk−2 .. .. .. .. .. .. . . . . . . ρk−1 ρk−2 ρk−3 · · · ρ1 1 Výběrová parciální autokorelační funkce Nahrazením ρi odhadem ρˆi v předcházejícím vztahu pro φkk získáme výběrovou parciální funkci φˆkk . Namísto komplikovaných počítání determinantů navrhl Durbin (1960) rekurzivní vztah [1]: φˆkk =
Pk−1 ˆ ˆk−j j=1 φk−1j ρ Pk−1 ˆ 1 − j=1 φk−1j ρˆj
ρˆk −
φˆkj = φˆk−1j − φˆkk φˆk−1k−j ,
j = 1, 2, ..., k − 1.
Výběrová ACF a výběrová PACF tedy slouží jako odhad pro ACF a PACF. Podle [3] lze ukázat, že za určitých předpokladů má rozdíl ρˆk − ρk asymptoticky normální rozdělení. Toho je samozřejmě možné využít k testování hypotézy H0 : ρk = 0 pro konkrétní k ≥ 1. Pokud data odpovídají realizaci Gaussovského bílého šumu, √ √ pak ρˆk leží s pravděpodobností 95% (α = 5%) v intervalu (u α2 / n, u1− α2 / n), kde up je kvantil normálního normovaného rozdělení. Tento interval lze nazvat intervalem spolehlivosti pro výběrové ACF a PACF. Leží-li koeficienty ACF a PACF v tomto intervalu, pak je lze považovat za nulové. Podmíněná heteroskedasticita V analýzách zabývajících se finančními časovými řadami, se často vychází z jednoho ze základních předpokladů, tj. že logaritmy tzv. výnosů (koeficienty růstu) mají normální rozdělení s konstantní střední hodnotou a konstantním rozptylem v čase. Logaritmování se provádí z toho důvodu, že ceny nemohou být záporné a předpokládá se tedy jejich logaritmicko-normální rozdělení. —4—
2: Časové řady typu ARCH
Ze skutečnosti vyplývá, že časové řady logarimů výnosů mají rozdělení špičatější s „tlustšímiÿ konci než je normální rozdělení. Tato vlastnost může být způsobena charakteristickým rysem jejich chování, resp. chováním logaritmů jejich koeficientů růstu. Variabilita (volatilita) těchto časových řad se v průběhu času mění, období s vysokou variabilitou střídají období s variabilitou nižší. To vyplývá z rostoucí a klesající nejistoty na trhu. Můžeme si představit, že logarimus výnosu má normální rozdělení s rozptylem, který se v čase mění. Tuto vlastnost nazveme podmíněnou heteroskedasticitou. Nepodmíněné rozdělení logaritmů výnosů je pak směsicí normálních rozdělení, kde ta s malým podmíněným rozptylem koncentrují výnosy v blízkosti střední hodnoty a naopak ta s velkým podmíněným rozptylem posouvají výnosy do konců rozdělení. Výstupem je nepodmíněné špičaté rozdělení s „tlustýmiÿ konci. Například z obrázku logaritmů výnosů (1), který je uveden také v kapitole 6.1, je patrné, že vysoká volatilita přechází do nízké volatility postupně a stejně tak naopak. To znamená, že volatilita určité úrovně se udržuje po určitý čas, volatilita v čase t tedy závisí na volatilitě v čase t − 1.
Obr. 1: Vývoj časové řady logaritmů výnosů ceny zlata Podmíněnou heteroskedasticitu můžeme vyjádřit jako: (ln Xt − ln Xt−1 )2 = α + ρ(ln Xt−1 − ln Xt−2 )2 + εt .
(1)
Pokud je parametr ρ = 0, pak podmíněná heteroskedasticita v časové řadě není. Naopak pokud ρ 6= 0 a zároveň |ρ| ≤ 1, pak podmíněná heteroskedasticita v časové řadě existuje. Čím je parametr |ρ| blíže k jednotce, tím významnější je autokorelace kvadrátů logaritmů [1].
—5—
2: Časové řady typu ARCH
2.1
Finanční časové řady
Finanční časové řady patří mezi vysokofrekvenční časové řady, které spadají pod krátkodobé časové řady, jejichž frekvence sledování je kratší než jeden rok. Sledovaným prvkem finančních časových řad je základní informace finančních trhů, kterou je cena. Může to být například cena akcie, cena měny, cena dluhopisu. Základní finanční časové řady vycházejí tedy z cen dosahovaných na veřejných trzích nebo charakterizují ceny a jejich vývoj. Základním rysem finančních časových řad je vysoká frekvence zaznamenávaných hodnot. Tyto hodnoty jsou zaznamenávány nejčastěji v měsíční, týdenní nebo denní frekvenci, ale mohou být zaznamenávány například na burzách i v hodinnových nebo minutových intervalech. Vliv na dynamiku takto rychlé frekvence záznamu dat mají jak systematické faktory (výrazně se projevuje trendová a cyklická složka), tak nesystematické faktory, což má za následek jejich vysokou a proměnlivou variabilitu. Z následující kapitoly vyplyne, že pro investora poskytuje výnos stejnou informaci jako samotná cena aktiva a navíc jej lze interpretovat v procentech, což může být jasnějším ukazatelem výhodnosti dané investice. Proto se při zpracování finančních časových řad většinou namísto cen aktiv pracuje s jejich výnosy.
2.2
Charakteristické rysy chování finančních časových řad
Základní a primární hypotézou o chování finančního trhu je hypotéza efektivního trhu2 . V literatuře lze nalézt například formulaci efektivního trhu takovou, že pokud ceny plně zahrnují očekávání a informace všech účastníků trhu, jsou jejich změny nepredikovatelné. S hypotézou efektivního trhu souvisí zřejmě nejstarší model „chováníÿ cen aktiv. Tento model je známý pod jménem martingál a lze ho popsat například takto: Jestliže Pt je cena aktiva v čase t, potom očekávaná cena v čase t + 1 je cena v čase t, za podmínky znalosti všech cen aktiva v minulosti, tj. v čase t − 1, t − 2, . . .. Z hlediska tvorby předpovědí model martingálu implikuje, že nejlepší (podle minimální střední čtvercové chyby) předpovědí zítřejší ceny je cena dnešní [1]. Představme si časovou řadu cen aktiva jako realizaci stochastického procesu {Pt }, tzn. řadu náhodných veličin uspořádaných v čase, pak martingál vyjádříme jako: E[Pt+1 |Pt , Pt−1 , ...] = Pt ,
(2)
E[Pt+1 − Pt |Pt , Pt−1 , ...] = 0,
(3)
nebo také jako:
kde E[.|.] je podmíněnou střední hodnotou. Dále model martingálu předpokládá, že nepřekrývající se cenové změny ve všech časových posunech (vpřed i vzad) nejsou korelované, tj. jsou lineárně nezávislé. Avšak mohou být stále nelineárně závislé. Tato závislost se může projevit tak, že cenové změny nejsou v čase konstantní. Dříve se předpokládalo, že martingál je nutnou podmínkou efektivního trhu. V současnosti se v souvislosti s podmínkou efektivního trhu hovoří nejen o nepredikovatelnosti, ale také o vztahu očekávané cenové změny a rizika (rozptylu). Z hlediska investora může jistě být zajímavé, že očekávaná cenová změna je kladná, i když je zatížená vysokým rizikem. Ale naopak z hlediska trhu je budoucí cena stále prakticky nepredikovatelná, a jedná se tedy o efektivní trh. Martingál klade restrikci 2 První
formulace Bachelier (1900), Cowles (1933).
—6—
2: Časové řady typu ARCH
pouze na očekávanou cenovou změnu nikoliv na riziko, nemůže proto být nutnou podmínkou efektivního trhu, jelikož mohou nastat situace, kdy lze trh stále považovat za efektivní, ale nejedná se již o model martingálu. Martingál můžeme také vyjádřit vztahem: Pt = Pt−1 + at ,
(4)
kde at se označuje jako přírůstek martingálu nebo také jako diference martingálu. Tato forma zápisu je podobná modelu náhodné procházky. Zde se na rozdíl od martingálu předpokládá, že {at } je proces bílého šumu3 . V tomto procesu jsou náhodné veličiny nejen nekorelované, ale také stejně rozdělené s nulovou střední hodnotou a konstantním rozptylem. Můžeme vycházet z představy, že rozdělení, které mají tyto náhodné veličiny, je normální, tj. at ∼ N (0, σa2 ). Toto je sice zřejmé, ale je zde zásadní nedostatek. Cena aktiva nemůže být menší než nula, minimální dosažitelný relativní přírůstek ceny neboli minimální dosažitelný jednoduchý výnos aktiva je tedy: Rt =
Pt − Pt−1 = −1. Pt−1
(5)
Protože náhodná veličina mající normální rozdělení může nabývat jakéhokoliv reálného čísla a ze vztahu (4) vyplývá, že cena aktiva Pt má normální rozdělení, není jeho dolní mez, a tedy ani dolní mez jednoduchého čistého výnosu zaručena. Těmto problémům lze předejít úvahou, že jednoduché výnosy aktiva definované jako: Rt + 1 =
Pt Pt−1
(6)
a nazývané také jako koeficienty růstu ceny aktiva, by měly mít rozdělení nezáporné náhodné veličiny. Pro tento případ se logicky nabízí rozdělení logaritmicko-normální4 . Tedy za předpokladu, že jednoduchý výnos Rt + 1 má logaritmicko-normální rozdělení, pak jeho logaritmus: Pt rt = ln (Rt + 1) = ln = ln Pt − ln Pt−1 = pt − pt−1 (7) Pt−1 má normální rozdělení. A výnos aktiva za k období od času t−k do času t lze vyjádřit jako součin k jednoduchých výnosů aktiva: Rt (k) + 1 = (Rt + 1) · (Rt−1 + 1) · ... · (Rt−k+1 + 1) = =
Pt Pt−1 Pt−2 Pt−k+1 Pt · · · ... · = . Pt−1 Pt−2 Pt−3 Pt−k Pt−k
(8)
Za předpokladu logaritmicko-normálního rozdělení jednoduchých výnosů má také celý výnos stejné rozdělení. Jeho logaritmická transformace má normální rozdělení a je rovna součtu k logaritmovaných jednoduchých výnosů [1], tj.: rt (k) = rt + rt−1 + rt−2 + ... + rt−k+1 .
(9)
Logaritmické výnosy jsou v aplikacích finančních časových řad často používané a tak si lze pro upřesnění namísto Xt uvažovat rt . 3 Obvykle se také vychází ze silnějšího předpokladu, že {a } je proces striktního bílého šumu, ve kterém jsou náhodné veličiny t nezávislé a stejně rozdělené s nulovou střední hodnotou a konstantním rozptylem. 4 Logaritmus náhodné veličiny s logaritmicko-normálním rozdělením má normální rozdělení.
—7—
2: Časové řady typu ARCH
2.2.1
Předpoklad normality
Jak bylo uvedeno v předcházející kapitole základním předpokladem, ze kterého se často vychází je, že logaritmy výnosů mají normální rozdělení s konstantní střední hodnotou µ a konstantním rozptylem σr2 , tj. předpoklad rt ∼ N (µ, σr2 ). O tomto rozdělení víme, že je symetrické a že jeho šikmost je rovna nule a jeho špičatost je rovna číslu tři5 . Avšak ve většině případů je odhad špičatosti logaritmů výnosů finančních časových řad větší než tři, to znamená, že skutečné rozdělení logaritmů výnosů je špičatější než normální rozdělení, tzn. že rozdělení výnosů má tedy sklon k leptokurticitě, neboli výnosy blízko střední hodnoty a výnosy velmi vysokých záporných a kladných hodnot se objevují častěji než by předpokládala normalita. Tato vlastnost bude demonstrována na reálných datech v aplikační části. Další odchylkou (i když ne tak výraznou) může být nesymetričnost rozdělení výnosů, tj. koeficient šikmosti není nulový. Z toho vyplývá, že kladné a záporné výnosy se nevyskytují se stejnou četností. Tento jev byl pozorován u většiny aktiv kromě směnných kurzů [1]. Není žádným překvapením, že i přes uvedené nedostatky je předpoklad normálního rozdělení výnosů poměrně běžný, protože ve prospěch předpokladu normálního rozdělení jednoznačně hovoří příznivé statistické vlastnosti, například v podobě snadného testování hypotéz nebo při odhadu parametrů pomocí metody maximální věrohodnosti. K těmto faktům pak někteří autoři poukazují na to, že i kdyby např. denní výnosy nebyly normálně rozdělené, pak tvrzení centrální limitní věty (rozdělení součtu velkého množství stejně rozdělených a nezávislých náhodných veličin s konečným rozptylem je přibližně normální) a z logaritmu vztahu (8) vyplývá, že řady výnosů s nižší frekvencí záznamu (měsíční, roční) konvergují k normálnímu rozdělení6 . Rychlost konvergence ale není příliš vysoká, i když jsou splněny všechny předpoklady centrální limitní věty. Přesto se však i ekonomové shodují, že pro delší frekvence záznamu cen aktiv se rozdělení výnosů blíží k normálnímu rozdělení (Například v [13], [12]). Naopak pro denní a menší frekvence výnosů už podobná shoda neexistuje a nejčastěji se uvádí srovnání například se Studentovo t-rozdělením7 . K-S test a Lilliefors normality test Oba tyto testy patří do kategorie testů dobré shody, pomocí kterých se zkoumá průběh celé distribuční funkce. Test Kolomogorova-Smirnova8 je obecný pro jakýkoliv typ rozdělení a Lillieforsova modifikace se používá pro test normálního rozdělení dat, kdy parametry neznáme a musíme je odhadnout přímo z dat pomocí výběrového průměru a rozptylu (využívá tedy jiných kritických hodnot než K-S test). K-S test je založen na výpočtu suprema vzdálenosti empirické a teoretické distribuční funkce a jeho následném porovnání s kritickou hodnotou K-S testu Dn (α), např. v [4]. Testujeme hypotézu takovou, že náhodný výběr X1 , ..., Xn pochází z rozložení s distribuční funkcí φ(x). Nechť Fn (x) je výběrová distribuční funkce. Testová statistika: Dn =
sup −∞<x<∞
|Fn (x) − φ(x)|.
5 Vzorce
pro šikmost a špičatost a jejich bodové odhady jsou uvedeny např. v [1, str. 20]. fakt je ověřen např. v práci [14]. 7 Také například Lévyho stabilní rozdělení, rozdělení s proměnlivým rozptylem. 8 Například zde: http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test. 6 Tento
—8—
2: Časové řady typu ARCH
Nulovou hypotézu zamítáme na hladině významnosti α, když Dn ≥ Dn (α), kde Dn (α) je tabelovaná kritická hodnota9 . Pro n ≥ 30 lze Dn (α) aproximovat výrazem: r 1 2 ln . 2n α Nulová hypotéza musí specifikovat distribuční funkci zcela přesně, tj. včetně všech jejich parametrů. Pokud je neznáme a odhadujeme je (Lillieforsův test), pak se nám mění i rozložení testové statistiky Dn . √ Pro hladinu významnosti 0.05 (resp. 0.01) má příslušná kritická mez asymptoticky platnou hodnotu 1, 358/ n √ (resp. 1, 628/ n). Lilliefors pak dále tuto hodnotu upravil pro případ testu normality rozdělení dat, kdy se nejdříve odhadnou teoretické parametry normálního rozdělení pomocí průměru a výběrové směrodatné odchylky. √ √ Asymptoticky platná kritická hodnota je pak 0, 89/ n (1, 04/ n) pro hladinu významnosti 0.05 (0.01). Shapiro-Wilk normality test Ve statistice se Shapiro Wilkův test používá pro testování hypotézy takové, že náhodný výběr X1 , ..., Xn pochází z normálního rozdělení s blíže nespecifikovanými parametry µ a σ 2 , N (µ, σ 2 ). Test je založen na porovnání bodů sestrojeného Q-Q grafu a regresní přímky proložené těmito body. Testovací statistika: W =
b2 = S2
k P
2 an−i+1 (yn−i+1 − yi )
i=1 n P
, (yi − y¯)
2
i=1
kde yi jsou uspořádané hodnoty náhodného výběru X1 , ..., Xn , an−i+1 jsou tabelizované váhy, y¯ je výběrový průměr a k =
n 2
je-li n sudé, resp. k =
n−1 2 ,
je-li n liché. Pokud je výsledek testu (W ) blízký jednotce
(pokud W = 1, pak jde o dokonalou shodu), p-value není dostatečně nízké a tedy hodnota testové statistiky nepřekročí tabelovanou kritickou hodnotu Shapiro-Wilkova testu, pak není možno zamítnout H0 předpokládající normální rozdělení, resp. čím je hodnota testovací statistiky blíže k jednotce, tím je lepší shoda mezi teoretickým a empirickým rozložením. Podle [10] postupujeme tímto způsobem: 1. Hodnoty jednotlivých pozorování seřadíme od nejmenší po největší tak, abychom získali uspořádaný výběr y = (y1 , ..., yn ), kde y1 ≤ · · · ≤ yn 2
2. Vypočteme S =
n P i=1
2
(yi − y¯) =
n P
(yi2 ) i=1
−
1 n
n P
2 yi
.
i=1
3. Vypočteme b podle vzorce b = an (yn − y1 ) + · · · + ak+2 (yk+2 − yk ) 4. Vypočteme testovou statistiku W =
b2 S2
5. Hodnotu testové statistiky porovnáme s tabelovanou kritickou hodnotou Shapiro-Wilkova testu a učiníme závěr o zamítnutí, resp. nezamítnutí nulové hypotézy na hladině významnosti α.
9 Například
zde: http://www.kmt.zcu.cz/person/Kohout/info_soubory/letnisem/tabulky.htm, odkaz platí i pro tabelizované váhy a kritické hodnoty Shapiro-Wilk testu.
—9—
2: Časové řady typu ARCH
Jarque-Bera test Tento test je testem dobré shody založený na porovnávání šikmosti a špičatosti s normálním rozdělením. Testové kritérium:
n 1 2 SK 2 + (KU − 3) , 6 4 kde SK je šikmost, KU je špičatost a n je počet pozorování. JB =
Kritická hodnota: χ21−α (ν), kde ν = 2 je počet stupňů volnosti. Pokud je JB < χ21−α , pak je přijímána hypotéza o tom, že data se řídí normálním rozdělením10 . 2.2.2
Předpoklad lineární závislosti
V klasických analýzách finančních časových řad předpokládáme, že logaritmy výnosů jsou nekorelované stejně rozdělené náhodné veličiny s nulovou střední hodnotou a konstantním rozptylem (proces bílého šumu), nebo nezávislé stejně rozdělené náhodné veličiny s nulovou střední hodnotou a konstantním rozptylem (proces striktního bílého šumu). Jednou ze společných vlastností logaritmů výnosů časových řad je proměnlivá variabilita neboli volatilita. Variabilita se mění ve velmi krátkých časových úsecích (pak se mezi nimi mohou objevit výnosy, které můžeme považovat za extrémní hodnoty) nebo může zůstat po určitou dobu téměř stabilní a pak se mění, tzn. variabilita se mění ve shlucích. Z proměnlivé variability logarimů výnosů vyplývají další zajímavé poznatky. Zjistilo se například, že variabilita může mít vliv na úroveň a sílu autokorelace v časových řadách výnosů. Zajímavé je, že bylo zjištěno, že vysoká variabilita následuje po záporném výnosu11 . Uvedené vlastnosti finančních časových řad a samozřejmě i jiné vlastnosti vyplývají z ekonomické podstaty trhu. Zároveň také z chování a jednání ekonomických subjektů (investorů). Klasické lineární modely však předpokládají pouze jeden typ závislosti (korelační závislost) a proto nemohou zachytit charakteristické rysy analyzovaných časových řad12 . K odhalení lineárních závislostí lze využít následující test, který se také používá k verifikaci modelů AR, MA, ARMA. Portmanteaův test (Ljung-Boxova statistika) Nechť X1 , ..., Xn jsou pozorované hodnoty časové řady {Xt } s autokorelační funkcí ρ(·). Za platnosti hypotézy: H0 : ρ1 = ρ2 = · · · = ρk = 0 platí, že: Qk = n(n + 2)
k X ρˆi ζ 2 − → χk , n −i i=1
k ∈ N,
ζ
kde − → značí konvergenci v distribuci pro n −→ ∞. Číslo k udává počet korelačních koeficientů, které testujeme. Výsledky ze simulačních studií ukázaly, že optimální volbou k je přibližně ln(n) [2]. 10 Například
zde: http://en.wikipedia.org/wiki/Jarque-Bera_test. informací je uvedeno v knize [1, str. 25], kde je tato skutečnost otestována. 12 Jiné typy závislostí umožnuje uvažovat třída nelineárních modelů. 11 Více
— 10 —
3: Modely volatility
3
Modely volatility
Ještě před popisem samotných modelů volatility by bylo dobré se seznámit s některými důležitými pojmy, které jsou obecně platné pro modelování různých typů ekonomických a finančních časových řad. Základem pro podávání informace o charakteru stochastického procesu je autokorelační funkce a parciální autokorelační funkce (např. modely stacionárních časových řad třídy AR, MA a ARMA jsou charakteristické specifickou formou autokorelační funkce a parciální autokorelační funkce, takže odhady těchto funkcí lze použít pro identifikaci modelu) [1].
3.1
Základní reprezentace
Základními charakteristikami náhodných veličin stochastického procesu jsou nepodmíněná a podmíněná střední hodnota a nepodmíněný a podmíněný rozptyl. Uvažujme například proces AR(1) ve formě: Xt = φ1 Xt−1 + at ,
(10)
kde |φ1 | < 1 a {at } je proces bílého šumu s nulovou střední hodnotou a rozptylem σa2 . Tento proces lze podle [1] vyjádřit ve formě: Xt = at + φ1 at−1 + φ21 at−2 + φ31 at−3 + . . . .
(11)
Z tohoto vyjádření je zřejmé, že E(Xt ) = 0. Podmíněná střední hodnota je střední hodnota veličiny Xt za předpokladu, že náhodné veličiny nabyly konkrétních hodnot v časech t − 1, t − 2, . . .. Podmíněná střední hodnota závisí na volbě podmínky, a je tedy její funkcí, tj. funkce regresní. V případě procesu AR(1) je podmínkou určitá hodnota náhodné veličiny v čase t − 1, E(Xt |Xt−1 ) = φ1 Xt−1 . Podmíněná střední hodnota veličiny Xt je tedy závislá na čase (proměnlivá v čase). Ze vztahu (11) vyplývá, že nepodmíněný rozptyl veličiny Xt je neměnný v čase, tj. var(Xt ) =
2 σa . 1−φ21
Pod-
míněný rozptyl se označuje jako funkce skedastická a její průběh charakterizuje měnlivost rozptylu veličiny Xt v závislosti na hodnotách veličin Xt−1 , Xt−2 , . . .. Ze vztahu (10) vyplývá, že podmíněný rozptyl veličiny Xt je var(Xt |Xt−1 ) = σa2 a je také v čase neměnný (označuje se jako funkce homoskedastická, v případě proměnlivého podmíněného rozptylu se označuje jako funkce heteroskedastická). Všechny tyto vlastnosti jsou charakteristické pro modely typu AR, MA a ARMA. Modely volatility vycházejí z představy, že např. model (10) lze modifikovat do tvaru: Xt = φ1 Xt−1 + εt ,
(12)
kde |φ1 | < 1 a {εt } je podmíněně heteroskedastický proces s podmíněnou střední hodnotou E(εt |Ωt−1 ) = 0 a podmíněným rozptylem var(εt |Ωt−1 ) = E(ε2t |Ωt−1 ) = ht , kde Ωt−1 je relevantní minulá informace až do času t − 1. Tyto požadavky splňuje například model procesu εt ve tvaru: 1/2
εt = et + ht ,
(13)
kde veličiny procesu {et } jsou nezávislé náhodné veličiny s nulovou střední hodnotou a jednotkovým rozptylem a ht je podmíněný rozptyl. Pokud je rozdělení náhodné veličiny et normované normální, za podmínky informace, která je k dispozici v čase t − 1, tj. et ∼ Nt−1 (0, 1), potom je rozdělení náhodné veličiny Xt také normální, za podmínky informace, která je k dispozici v čase t − 1, ale s podmíněným rozptylem, který se mění v závislosti na čase, tj. Xt ∼ Nt−1 (0, ht ). — 11 —
3: Modely volatility
3.2
Lineární modely volatility
Modely volatility poprvé popsal Engle(1982)13 . Tyto modely a modely z nich odvozené se označují jako lineární, protože jsou charakteristické tím, že podmíněný rozptyl je lineární funkcí veličin ε2t−1 , ε2t−2 , . . . , ε2t−q . V dalších částech budou popsány jednotlivé modely volatility, které formulují vývoj podmíněného rozptylu ht v čase a jejich úplný tvar se bude vždy řídit podle modifikace zvoleného úrovňového modelu, tj. například (12), kde εt ∼ N (0, ht ), jak bylo uvedeno v předchozí kapitole, proto jsou formulovány pouze vývoje podmíněného rozptylu. 3.2.1
Modely ARCH (Autoregressive Conditional Heteroscedasticity)
ARCH(1) Podmíněný rozptyl tohoto modelu má tvar: ht = ω + α1 ε2t−1 ,
(14)
kde ω > 0 a α1 ≥ 0, aby bylo zaručeno, že podmíněný rozptyl je kladné číslo. Dále platí, že když α1 = 0, pak je podmíněný rozptyl v čase konstantní a proces {εt } se označuje za podmíněně homoskedastický. Pokud k oběma stranám rovnice (14) přičteme ε2t a odečteme ht , dostaneme autoregresivní tvar modelu: ε2t = ω + α1 ε2t−1 + νt ,
(15)
kde νt = ε2t − ht = ht (e2t − 1) a εt ∼ N (0, 1). Podmíněná a nepodmíněná střední hodnota procesu {νt } je nulová (tj. proces není autokorelovaný). Z teorie autoregresních procesů vyplývá, že proces ARCH(1) je stacionární v kovariancích, je-li α1 < 1, pak nepodmíněná střední hodnota procesu ε2t , to je nepodmíněný rozptyl procesu εt , má tvar: var(εt ) =
ω . 1 − α1
(16)
Je tedy konstantní v čase a proces εt je nepodmíněně homoskedastický, autokorelační funkce je ve zpoždění k rovna α1k . Prostřednictvím modelu ARCH lze zachytit shluky volatility v časové řadě. To vyplývá z autoregresivního tvaru modelu (15). Je-li εt−1 v absolutní hodnotě vysoké, potom lze očekávat, že εt v absolutní hodnotě bude vysoké. Tedy vysoké (resp. nízké) hodnoty časové řady budou následovány vysokými (resp. nízkými) hodnotami jakéhokoliv znaménka. Tento model umožňuje také zachytit vyšší špičatost pravděpodobnostního rozdělení. ARCH(q) Podmíněný rozptyl tohoto modelu má tvar: ht = ω + α1 ε2t−1 + α2 ε2t−2 + . . . + αq ε2t−q .
(17)
Podmínky ω > 0 a αi ≥ 0 pro i = 1, 2, . . . , q zaručují kladný rozptyl. Opět lze převést vztah (17) na autoregresní tvar (konkrétně je to tvar modelu AR(q) procesu {ε2t }): ε2t = ω + α1 ε2t−1 + α2 ε2t−2 + . . . + αq ε2t−q + νt , 13 http://finance.martinsewell.com/arch-garch/Engle2001.pdf
— 12 —
(18)
3: Modely volatility
kde νt = ε2t − ht . Nepodmíněný rozptyl procesu {εt } má tvar: var(εt ) =
ω , 1 − α1 − . . . − αq
(19)
a to znamená, že je konstantní v čase a proces {εt } je nepodmíněně homoskedastický. Následující obrázek ukazuje příklad (jedné realizace) nasimulovaného modelu ARCH(2) s nastavenými parametry α1 = 0.15, α2 = 0.34.
Obr. 2: Simulace modelu ARCH(2)
R code: 1 # Arch simulace 2 k = 1000 3 4 # ARCH(2) − use default omega and specify alpha, set beta=0! 5 spec = garchSpec(model = list(alpha = c(0.15, 0.34), beta = 0)) 6 sim = garchSim(spec, n = k) 7 8 # Graf sumulace ARCH(1) 9 ggplot(data.frame(Time=seq(1:k), Innovations = sim$garch), aes(Time, Innovations)) 10 + geom line(colour=”blue”, size = 1) 11 + geom point(colour = ”black”, size = 1) 12 + labs(title = ”ARCH(2) simulation (k = 1000)”)
— 13 —
3: Modely volatility
3.2.2
Modely GARCH (Generalized ARCH)
GARCH(1,1) Při modelování časových řad, u kterých je třeba použít modely ARCH(q) s vysokým číslem q, vzniká problém s odhadováním velkého množství parametrů, umocněný existencí řady omezujících podmínek14 . Tomuto problému se můžeme vyhnout rozšířením modelu ARCH o zpožděný podmíněný rozptyl15 . Rozšíření modelu ARCH(1) o podmíněný rozptyl v prvním zpoždění má model podmíněného rozptylu tvar: ht = ω + α1 ε2t−1 + β1 ht−1 .
(20)
Podmínky ω > 0, α1 > 0 a β1 ≥ 0 zaručují kladný podmíněný rozptyl. Model (20) se označuje jako GARCH(1,1) a lze ho použít tam, kde by bylo vhodné volit model ARCH s mnoha zpožděními. Abychom blíže nastínili myšlenku modelů GARCH, přepišme vztah (20) takto: ht = ω + α1 ε2t−1 + β1 ht−1 = ω + (α + β)ht−1 + α(ε2t−1 − ht−1 ).
(21)
Podmíněný rozptyl v této podobě se pak rovná váženému součtu rozptylu ht−1 předpovězeného v minulém období a neočekávaného minulého šoku ε2t−1 − ht−1 . Parametr α měří vliv tohoto šoku na předpověď pro další období, (α + β) pak představuje rychlost, s jakou vliv šoku v následujících obdobích vymizí. Čím více se (α + β) blíží jedné, tím je čas doznívání šoku delší. Pokud přičteme ε2t k oběma stranám modelu (20) a odečteme ht , lze tento model přepsat do tvaru modelu ARMA(1,1): ε2t = ω + (α1 + β1 )ε2t−1 + νt − β1 νt−1 ,
(22)
kde νt = ε2t − ht . Jestliže α1 + β1 < 1, pak z tohoto zápisu vyplývá, že model GARCH(1,1) je stacionární v kovariancích. Nepodmíněný rozptyl procesu {εt } má tvar: var(εt ) =
ω 1 − α1 − β1
(23)
a je tedy konstantní v čase a proces {εt } je nepodmíněně homoskedastický. Podmínka α1 > 0 zaručuje, aby nenastala situace, že α1 = 0. Důvod této podmínky vyplývá z modelu (22). Nebylo by pak možné identifikovat parametr β1 16 . Bollerslev (1988) také odvodil rekurentní vztah pro autokorelační funkci procesu ε2t : ρ1 = α1 +
α12 β1 , 1 − 2α1 β1 − β12
ρk = (α1 + β1 )k−1 ρ1 ,
k = 2, 3, . . . .
(24) (25)
Hodnoty autokorelační funkce s rostoucím zpožděním k exponenciálně klesají. Rychlost klesání těchto hodnot záleží na hodnotě součtu α1 + β1 . Pokles autokorelační funkce je velmi pozvolný pokud se tento součet blíží k hodnotě jedna. Hodnoty parciální autokorelační funkce se chovají obdobně. Obecně tedy lze konstatovat, že tvar ACF a PACF procesu ε2t odpovídá tvaru těchto funkcí modelu ARMA(p,q). 14 Např.
nenulovost podmíněného rozptylu, stacionarita. rozšíření navrhl Bollerslev (1986). 16 Podrobněji viz.[1, str. 167]. 15 Toto
— 14 —
3: Modely volatility
GARCH(p,q) Podmíněný rozptyl modelu GARCH má obecný tvar: ht = ω +
q X
αi ε2t−i +
i=1
p X
βi ht−i .
(26)
i=1
Podmínky ω > 0, αi > 0 pro i = 1, 2, . . . , q a βi ≥ 0 pro i = 1, 2, . . . , p zaručují kladný podmíněný rozptyl. Model (26) se označuje jako GARCH(p,q) a lze ho přepsat do tvaru: ε2t = ω +
m X
(αi + βi )ε2t−i −
i=1
p X
βi νt−i + νt ,
(27)
i=1
kde νt = ε2t − ht . Nepodmíněný rozptyl procesu {εt } je: var(εt ) = kde α(1) =
q X
αi
ω , 1 − α(1) − β(1)
a
β(1) =
i=1
p X
(28)
βj .
j=1
Další typy a modifikace V literatuře lze nalézt další typy těchto modelů. Zde jsou uvedeny jen ty nejdůležitější pro aplikační část této práce. Mimojiné byly vyvinuty další modifikace těchto modelů pro zachycení nejrůznějších typů asymetrických efektů (QGARCH, VS-GARCH nebo ANST-GARCH), které lze nalézt např. v [1].
3.3
Výstavba modelů volatility
Výstavbou modelů volatility rozumíme určitý postup pro správnou aplikaci modelů volatility. Tento postup lze obecně shrnout do těchto kroků: 1. Pro časovou řadu se určí vhodný lineární nebo nelineární úrovňový model. 2. Analyzuje se kvalita úrovňového modelu z hlediska výsledku testu podmíněné heteroskedasticity a normality, případně se otestuje hypotéza podmíněné heteroskedasticity nelineárního typu. 3. Odhadnou se parametry lineárního nebo nelineárního modelu podmíněné heteroskedasticity. 4. Ověří se vhodnost daného modelu diagnostickými testy. 5. Pokud je třeba, model se dále modifikuje. 6. Konečný model se použije k popisným nebo predikčním účelům.
3.4
Testování podmíněné heteroskedasticity
Test podmíněné heteroskedasticity lineárního typu vychází z formulace modelu ARCH. Podmíněný rozptyl ht modelu ARCH(q) ve tvaru (17) je konstantní, jsou-li parametry odpovídající veličinám ε2t−1 , . . . , ε2t−q rovny nule. Nulová hypotéza, tj. hypotéza podmíněné heteroskedasticity, je H0 : α1 = α2 = . . . = αq = 0, alternativní hypotézou je, že alespoň jeden parametr je různý od nuly, tj. H1 : nonH0 . — 15 —
3: Modely volatility
Test je prováděn v následujících krocích: 1. Odhadnou se parametry lineárního nebo nelineárního úrovňového modelu a získají se rezidua εbt a reziduální P 2 součet čtverců RSE0 = εbt . 2. Konstruuje se regresní model εb2t = ω + α1 εb2t−1 + α2 εb2t−2 + . . . + αq εb2t−q + ut ,
(29)
z kterého se získá reziduální součet čtverců RSE1 a index determinace R2 . 3. Testové kritérium LM ve tvaru T ·R2 má za předpokladu platnosti nulové hypotézy asymptoticky rozdělení χ2 (q). 4. F -verze tohoto testového kritéria pro malé výběry má tvar FLM =
(RSE0 − RSE1 )/q , RSE1 /(T − q − 1)
(30)
její rozdělení lze za předpokladu platnosti nulové hypotézy aproximovat rozdělením F (q, T − q − 1). Tento test se často označuje jako ARCH test. Lze jej také interpretovat jako test autokorelace čtverce nesystematické složky. ARCH-LM test Pro řadu {εt } se zkonstruuje lineární regresní model, kde jako vysvětlující proměnné vystupuje q zpožděných hodnot εt−1 , εt−2 , ..., εt−q , tj.: εt = ω + α1 εt−1 + α2 εt−2 + · · · + αq εt−q + νt . Dále pak odhadneme parametry modelu a vypočítáme koeficient determinace R2 . Za platnosti nulové hypotézy: H0 : αi = 0 pro i = 1, ..., q má statistika LM = T · R2 asymptoticky χ2 (q) rozdělení s q stupni volnosti. Tento test lze interpretovat tak, že pokud by pro model ARCH(q) ve tvaru (17) platilo α1 = ... = αq = 0, byl by jeho podmíněný rozptyl konstantní.
3.5
Volba řádu modelu
Při výstavbě modelu typu ARCH narazíme na zásadní problém a tím je volba řádu modelu. Tento řád lze přibližně odhadnout z průběhu autokorelační funkce nebo parciální autokorelační funkce řady {ε2t }. Avšak žádný formální test, podle kterého by bylo možné určit řád modelu není k dispozici. Obvyklý postup při určování řádu modelu typu ARCH je takový, že se nejprve odhadne model nízkého řádu a poté se tento řád upraví např. podle výsledků statistické významnosti parametrů nebo analýzy standardizovaných reziduí. Ve velké většině případů jsou postačující modely nízkých řádů např. ARCH(1), ARCH(2) nebo GARCH(1,1), GARCH(2,1) atd.
— 16 —
3: Modely volatility
3.6
Odhad parametrů
Odhady parametrů modelů volatility se obvykle provádějí pomocí metody maximální věrohodnosti17 . V aplikacích modelu ARCH se běžně uvažuje normované normální rozdělení ∼ N (0, 1) nebo standardizované Studentovo rozdělení ∼
√ zt
varzt , kde zt ∼ tν . Je samozřejmě možné využít i jiná rozdělení, například v [5].
Uvažujme tedy model ARCH(m) a ∼ N (0, 1). Dále pak označme pozorované hodnoty y1 , ..., yT náhodných veličin Y1 , ..., YT . Podle vlastností podmíněné hustoty využijeme: f (y1 , ...yT |α) = f (yT , ...ym+1 |α, Fm ) · f (ym , ...y1 |α) = = f (yT , ...ym+2 |α, Fm+1 ) · f (ym+1 ||α, Fm ) · f (ym , ...y1 |α) = .. . = f (yT |α, FT −1 ) · · · f (ym+1 |α, Fm ) · f (ym , ...y1 |α), kde f (y1 , ...yT |α) je sdružená hustota veličin Y1 , ..., YT v bodě y1 , ..., yT , která je podmíněná vektorem parametrů α = (α1 , ..., αm ) a f (yt |α, Ft−1 ) je hustota náhodné veličiny Yt podmíněná vektorem parametrů α a σ-algebrou Ft−1 = σ(Yt−1 , ..., Y1 ) generovanou veličinami Yt−1 , ..., Y1 . Protože ∼ N (0, 1) a Yt = σt t , jsou podmíněná rozdělení Yt |Ft−1 ∼ N (0, σt2 ) a tedy: f (yt |Ft−1 ) = p
1 2πσt2
e
−
2 yt 2 2σt
.
2 2 , nebo-li σt2 = σt2 (α). Pak můžeme psát: + · · · + αm Yt−m Stále uvažujeme model, kde σt2 = α0 + α1 Yt−1
f (y1 , ...yT |α) = p
1 2πσt2
e
−
2 yt 2 2σt
f (ym , ...y1 |α).
A protože vyjádření hustoty f (ym , ...y1 |α) by bylo složité, podmíníme sdruženou hustotu f (y1 , ...yT |α) hodnotami y1 , ...ym a budeme je odhadovat společně s α. Dostaneme potom: T Y
f (y1 , ...yT |α, y1 , ...ym ) = f (y1 , ...yT |α, Fm ) =
1 p
t=m+1
2πσt2
e
−
2 yt 2 2σt
.
Metoda maximální věrohodnosti tedy staví na myšlence, že náš výběr z rozdělení s hustotou f je „nejvíce pravděpodobnýÿ, a tedy hustota f je v bodě určeném naším výběrem maximální. Tato hustota závisí na neznámém vektoru α. Hledáme tedy takové jejich hodnoty, aby hustota byla maximální, které jsou pak odhadem metodou maximální věrohodnosti. Druhým případem je, že t má standardizované Studentovo t-rozdělení s ν stupni volnosti. V tomto případě √ zt
ν varzt , kde zt ∼ tν . Z vlastností Studentova rozdělení víme, že varzt = ν−2 pro ν > 2. Hustotu t lze získat jako hustotu lineární transformace zt s využitím věty o transformaci náhodných veličin (viz. [6]).
je t ve tvaru
Z tohoto vyplývá: − ν+1 2 Γ ν+1 x2 2 ft (x) = 1+ , p ν ν−2 Γ 2 (ν − 2) π kde ν > 2, x ∈ R a Γ(x) =
R∞ 0
y x−1 e−y dy pro x > 0. Podobně jako v případě t ∼ N (0, 1) dostáváme:
f (y1 , ...yT |α, Fm ) =
T Y t=m+1
17 Například
Γ
Γ p ν 2
ν+1 2
(ν − 2) πσt2
v knize [1, str. 190] je tato metoda podrobně popsána.
— 17 —
1+
x2 (ν − 2) σt2
− ν+1 2 .
3: Modely volatility
Tuto funkci maximalizujeme, pokud je počet stupňů ν zadán a jedinými neznámými parametry jsou α. Pokud není parametr ν zadán, pak je přidán do podmínky podmíněné hustoty a odhadne se spolu s α. V praxi se ukázalo, že odhadnutý počet stupňů volnosti je nejčastěji mezi hodnotami 3 a 6 (uvedeno například v [2]). Pokud bychom potřebovali počáteční nastavení tohoto parametru např. pro nějaký iterační proces, volba tohoto parametru by byla mezi uvedenými čísly. Výsledný odhad má tvar: 2 2 σ ˆt = α ˆ0 + α ˆ 1 Yt−1 + ··· + α ˆ m Yt−m .
Je dobré si také uvědomit, že odhady parametrů dané časové řady pro model volatility mohou být zkreslené. To může být způsobeno např. odlehlými pozorováními, které také mohou zkreslit výsledky testů podmíněné heteroskedasticity. Těmto problémům lze předejít pomocí metod robustních odhadů parametrů, které jsou vzhledem k odlehlým pozorováním robustní18 . Tento přístup však není cílem práce.
3.7
Verifikace modelu
Ověřit platnost a vhodnost modelu pro konkrétní časovou řadu můžeme několika způsoby: 1. Testování nesystematické složky - konkrétně se jedná o testy autokorelace (např. portmanteau uvedený v kapitole 2.2.2), normality (např. Jarque-Bera test uvedený v kapitole 2.2) a podmíněné heteroskedasticity (Ljung-Boxův Q-test19 nebo ARCH-LM test uvedený v kapitole 3.4) 2. Testování linearity - pro testování lineárního modelu GARCH proti nějakému nelineárnímu modelu volatility lze použít testy jako jsou SB, NSB nebo PSB20 . Tyto testy souvisejí s testováním podmíněné heteroskedasticity nelineárního typu. Velkou vahou ve správném specifikování modelu je to, že pro jeho standardizovaná residua by mělo platit ebt ∼ IID(0, 1)21 . To lze ověřit právě testem ARCH. 3.7.1
Kritéria hodnocení modelu
Další posouzení vhodnosti určitého modelu lze provést několika způsoby. Jedním z nich je porovnáním kritérií, které určitým způsobem modely hodnotí. V této části uvedeme kritéria AIC a BIC, která jsou definována jako22 : AIC = −2lnL(ϕ) ˆ + 2M BIC = −2lnL(ϕ) ˆ + M ln(T ), kde lnL(ϕ) ˆ je hodnota maximalizované logaritmické věrohodnostní funkce, M = p + 1 a p je počet parametrů modelu, T je počet pozorování23 . Bylo dokázáno, že AIC vede k nadhodnocení řádu autoregrese ([1, str. 105], Akaike (1979)). Vybere se takový model, který vede k minimální hodnotě těchto kritérií. 18 Jedna
z metod je uvedena v knize [1, str. 193]. například: http://en.wikipedia.org/wiki/LjungBox_test. 20 Tyto testy jsou popsány například v knize [1, str. 187]. 21 IID = independent and identically distributed (nezávislé a stejně rozdělené). 22 Akaikeho informační kritérium (Akaike Information Criterion), Bayesovo informační kritérium (Bayesian Information Criterion) 23 Podle manuálu pro R: http://tolstoy.newcastle.edu.au/R/e2/help/06/09/1750.html.
19 Viz.
— 18 —
3: Modely volatility
3.8
Konstrukce předpovědí
Pokud bychom uvažovali proces ARMA(P,Q)-GARCH(p,q) ve tvaru: εt ∼ GARCH(p, q),
φ(B)Xt = Θ(B)εt ,
(31)
ˆ t+k nejlepší (ve smyslu minimální kde φ(B) = 1−φ1 B −· · ·−φP B P , Θ(B) = 1+Θ1 B +· · ·+ΘQ B Q . Označme X ˆ t+k = Et (Xt+k ). střední kvadratické chyby) k-krokovoupredikci tohoto procesu konstruovanou v čase t, tj. X ˆ t+k vliv. Naopak se ale projeví Lze ukázat, že podmíněná heteroskedasticita nemá na bodové předpovědi X na velkosti chyb předpovědí a předpovědních intervalů. Pro výpočet kvadratické chyby předpovědi budeme předpokládat, že kořeny polynomu 1 − φ1 z − · · · − φP z P , z ∈ C leží vně jednotkového kruhu. Pak můžeme vyjádřit proces Xt v kauzálním tvaru: ∞
Xt =
X Θ(B) εt = ψj εt−j . φ(B) j=0
Tento proces v čase t + k je možné zapsat ve formě: Xt+k =
∞ X
ψj εt+k−j ,
j=0
ˆ t+k je tedy: protože Et (εt+k−j ) = 0 pro j < k, nejlepší předpověď X ˆ t+k = X
∞ X
ψj εt+k−j .
j=k
Chybu předpovědi můžeme vypočítat jako: ˆ t+k = rt,k = Xt+k − X
k−1 X
ψj εt+k−j
j=0
a minimální nepodmíněná střední kvadratická chyba (Mean Square Error, MSE) predikce je potom: ˆ t+k ) = E(r2 ) = var(rt,k ) = σ 2 M SE(X t,k ε
k−1 X
ψj2 ,
(32)
j=0
kde σε2 značí nepodmíněný rozptyl procesu εt . 2 ˆ t+k ) = Et (rt,k M SEt (X )=
k−1 X
ψi2 Et (ht+k−j ),
(33)
j=0
kde Et (ht+k−j ) je nejlepší předpověď podmíněného rozptylu ht+k−j . Při výpočtu této předpovědi můžeme postupovat následovně: nechť ht je podmíněný rozptyl procesu GARCH(p,q) ve tvaru: ht = ω + α1 ε2t−1 + · · · + αq ε2t−q + β1 ht−1 + · · · + βp ht−p . Hodnotu tohoto podmíněného rozptylu v čase t + k můžeme zapsat jako: ht+k = ω + α1 ε2t+k−1 + · · · + αq ε2t+k−q + β1 ht+k−1 + · · · + βp ht+k−p , — 19 —
3: Modely volatility
ˆ t+k = Et (ht+k ) má tvar: potom nejlepší předpověď h ˆ t+k = ω + α1 εˆ2 h t+k−1 + · · · + ˆ t+k−1 + · · · + βp h ˆ t+k−p , αq εˆ2t+k−q + β1 h kde ˆ t+j = Et (ht+j ) = εˆ2 = Et (ε2 ) h t+j t+j ˆ t+j = Et (ht+j ) = ht+j h εˆt+j = Et (ε2t+j ) = ε2t+j
j ≥ 1, j ≤ 0, j ≤ 0.
Dále budeme uvažovat předpověď vypočítanou na základě modelu GARCH(1,1). Podmíněný rozptyl ht je: ht = ω + αε2t−1 + βht−1 , ˆ t+k má potom tvar: předpověď h ˆ t+k = ω + αˆ ˆ t+k−1 , h ε2t+k−1 + β h a opakovanou substitucí lze získat: ˆ t+k = ω h
k−1 X
(α + β)j + (α + β)k−1 αε2t + (α + β)k−1 βht
j=0
=ω
k−2 X
(α + β)j + (α + β)k−1 ht+1 .
j=0
Jestliže je proces GARCH(1,1) stacionární, tzn. že platí α + β < 1 a σε2 = ω/(1 − α − β), lze předchozí tvar vyjádřit jako: ˆ t+k = σ 2 + (α + β)k−1 (ht+1 − σ 2 ) h ε ε a pokud k −→ ∞, předpověď podmíněného rozptylu konverguje k nepodmíněnému rozptylu σε2 (více například v [1] nebo [9]).
— 20 —
4: Monte Carlo
4
Monte Carlo
Metoda Monte Carlo je poměrně široká skupina numerických výpočetních metod, která je založena na využití náhodných veličin a teorie pravděpodobnosti. Můžeme říci, že jde o simulaci pomocí stochastických metod, které využívají pseudonáhodná čísla. Tato metoda je široce využívaná např. pro simulace experimetnů, počítání integrálů, výpočtů čísla π, řešení diferenciálních rovnic atd. Simulování je často jednoduchou cestou k řešení velmi složitých, ale praktických úloh. Jinak řečeno základním principem metody Monte Carlo je simulace velkého množství možných budoucích stavů, jejichž počet může jít až do stovek tisíců. Následně se pak provádí analýza hlavních charakteristik výsledků. Tato metoda spadá do skupiny metod, které jsou využívané v případech, kdy matematické řešení problému je buď příliš složité nebo jej dokonce nelze nalézt. Do skupiny úloh řešitelných pomocí metody Monte Carlo můžeme zařadit téměř libovolnou funkci, která převádí vstupní veličiny X na výstupní Y . Co je potřeba znát je distribuční funkce vstupních veličin X, pomocí které se generují náhodné vstupní veličny a tak zaznamenáváme výstupní veličiny. Po dostatečném24 počtu opakování, lze pomocí statistické analýzy výstupů odhadnout hledaný parametr. Pro Monte Carlo existují dvě varinaty: • Geometrický model - v tomto případě se při výpočtu nepoužívá model reálného děje (př. výpočet integrálu), • Analogový model - musíme znát pravděpodobnostní rozdělení zkoumaných jevů. Opakovanou simulací pak získáváme výsledky (náhodné veličiny). Pro generátory náhodných čísel, které se využívají při této metodě, je potřeba využít nějaký hardwarový generátor. Ty jsou většinou založeny na nějakých náhodných procesech (např. se využívá aktuální datum a čas, měření nějakého šumu, např. polovodičů či kolísání napětí, rozpad atomového jádra izotopu). Metodu Monte Carlo lze využít v mnoha aplikacích a modelech. O této metodě je sepsáno mnoho článků a knih v různých aplikacích či verzích. Další informace o této metodě lze nalézt např. v [7] nebo [8].
4.1
Předpovědi modelů GARCH pomocí MC
V této části bude stručně popsáno, jakým způsobem lze pomocí metody Monte Carlo předpovědět budoucí hodnotu časové řady modelované pomocí GARCH modelu. Uveďmě nejdříve, co v tomto případě předchází použití metody MC: • Pro stochastickou časovou řadu máme již odhadnutý a verifikovaný model podmíněné střední hodnoty (ARMA) a odhad εˆt , • Pro řadu residuí εˆt máme odhadnutý a verifikovaný model podmíněného rozptylu (GARCH), • Pro standardizovaná residua ˆt modelu GARCH jsme určili jakým rozdělením se řídí (nejčastěji normální nebo studentovo rozdělení). V tuto chvíli je vše připraveno pro aplikaci metody Monte Carlo. 24 Zde
se objevuje problém, co znamená dostatečný. Obecně a logicky však platí, že přesnost odhadu pomocí metody Monte Carlo roste s počtem opakování.
— 21 —
4: Monte Carlo
Jak už bylo uvedeno výše, metoda je založena na generování náhodných čísel a budeme tedy postupovat v následujících krocích: • Budeme tedy generovat náhodná čísla R, resp. pravděpodobnosti, v intervalu h0, 1i. • Dále budeme potřebovat modelovou distribuční funkci p = F (ˆ t ), resp. inverzní distribuční funkci neboli kvantilovou funkci ˆt = F −1 (p) rozdělení residuí ˆt (většinou kvantilová funkce normálního nebo studentova rozdělení). • Pro generovanou náhodnou pravděpodobnost R určíme pomocí kvantilové funkce první simulovanou hodnotu ˆ1t+1 = F −1 (R). • Tím získáme jednu simulovanou přepověď ˆ1t+1 = F −1 (R) pro standardizovaná residua GARCH modelu. Pokud uvažujeme např. model GARCH(1,1), pak simulovanou předpověď pro t + 1 podmíněného rozptylu ˆ t+1 = ω ˆ t. bude h ˆ +α ˆ 1 εˆt + βˆ1 h q ˆ t+1 , tj. simulované residuum ARMA • Poté lze určit jednu simulovanou předpověď pro εˆ1t+1 = ˆt+1 · h modelu v čase t + 1. • Tento postup pak opakujeme k-krát, abychom dostali k náhodných simulací předpovědi v čase t + 1, tj. ˆ1t+1 , ..., ˆkt+1 a dále pak εˆ1t+1 , ..., εˆkt+1 . Pro delší předpověď až do času t + h lze postupovat krok po kroku, kdy popsaný postup opakujeme, ale až poté, co dokončíme simulaci do času t + h. Tedy v momentě, kdy pomocí náhodného čísla získáme odhad εˆ1t+1 , vygeˆ t+2 za pomoci nerujeme další náhodnou pravděpodobnost R, pro kterou určíme ˆ1 = F −1 (R) a také určíme h t+2
εˆ1t+1 , který jsme odhadli v přechozím kroku. Pomocí těchto simulovaných odhadů můžeme dále vypočítat εˆ1t+2 . Takto bychom pokračovali až do času t + h a poté celý tento postup opakovali a generovali řadu ˆ2t+1 , ..., ˆ2t+h , resp. εˆ2t+1 , ..., εˆ2t+h . Podle zvolené hodnoty k počtu simulací pak postup opakujeme k-krát a tedy poslední získanou řadou simulovaných hodnot bude ˆkt+1 , ..., ˆkt+h , resp. εˆkt+1 , ..., εˆkt+h .
— 22 —
5: Project R
5
Project R
R je jazyk a prostředí pro statistické výpočty a grafiku. Je to GNU projekt, který je podobný jazyku S a prostředí, které bylo vyvinuto v Bell laboratories (dříve AT & T, nyní Lucent Technologies) Johnem Chambersem a kolegy. Projekt R nabízí širokou řadu statistických funkcí (lineární a nelineární modelování, klasické statistické testy, analýzy časových řad, klasifikace, clustering, . . .) a grafické techniky, a je vysoce rozšiřitelný. Jazyk S je často preferovaným prostředkem pro výzkum statistické metodiky a R poskytuje Open Source cestu k účasti na této činnosti. Jednou ze silných stránek R je snadnost, s níž mohou být produkovány dobře navržené a publikovatelně kvalitní grafy, včetně matematických symbolů a vzorců. R je k dispozici jako Free software v souladu s podmínkami GNU General Public License (Free Software Foundation) ve formě zdrojového kódu. Běží na nejrůznějších platformách UNIX a podobných systémech (včetně FreeBSD a Linux), Windows a MacOS. R je integrovaná sada softwarových zařízení pro manipulaci s daty, výpočty a grafické zobrazení. Zahrnuje tedy: • efektivní zpracování dat a skladovací prostory, • sadu operátorů pro výpočty polí, zejména matic, • velkou, soudržnou, integrovanou kolekci dílčích nástrojů pro analýzu dat, • grafické zařízení pro analýzu dat a zobrazení buď na obrazovce nebo pro tisk, • dobře vyvinutý, jednoduchý a efektivní programovací jazyk, který zahrnuje podmínky, smyčky, uživatelem definované rekurzivní funkce a vstupní a výstupní zařízení. Termín „prostředíÿ je určen k charakterizování plně plánovaného a soudržného systému, než přídavného narůstání velmi specifických a neflexibilních nástrojů, jak je tomu s jinými softwary pro analýzy dat. Mnoho uživatelů smýšlí o R jako o statistickém systému. Vývojáři dávají přednost spíše pojmu prostředí, ve kterém jsou prováděny statistické metody. R může být rozšířen (snadno) pomocí balíčků. Existuje osm balíků dodávaných s distribucí R a mnoho dalších je k dispozici prostřednictvím internetových stránek skupiny CRAN obsahujicích velmi širokou řadu moderních metod statistiky. R má svou vlastní dokumentaci ve formátu LaTeX-like, který se používá k zásobování komplexní dokumentace, a to jak on-line formátů tak tištěných podobách [18].
5.1
Počáteční nastavení a načtení dat
Ze stránek http://www.r-project.org/ lze stáhnout instalační soubor, který nám umožní nainstalovat Project R25 . Po instalaci spustíme Project R a do hlavní konzole napíšeme příkazy: R code: 1 install.packages(”Rcmdr”) 2 library(Rcmdr) 25 Podrobnější
informace a návody lze nalézt například na stránkách http://www.karlin.mff.cuni.cz/~kulich/vyuka/Rdoc/
index.html.
— 23 —
5: Project R
Prvním příkazem nainstalujeme balíček, který umožní použití uživatelsky příjemnější konzole. Druhým příkazem tuto konzoli spustíme. V této konzoli načteme vytvořený soubor v projektu R, který umožní zpracování časové řady. Pomocí sekvence příkazů ”File” - ”Open script file” - ”. . \GARCH” - MainGARCH.R otevřeme zdrojový soubor, který je na CD přiloženém k této diplomové práci. V první části tohoto skriptu je načtení balíčků potřebných v dalších částech skriptu a také nastavení pracovního adresáře a načtení dat. R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
# Instalace použitých balíčků install.packages(”tseries”) install.packages(”FinTS”) install.packages(”nortest”) install.packages(”ggplot2”) install.packages(”xtable”) install.packages(”lmtest”) install.packages(”fGarch”) install.packages(”reshape”) install.packages(”quantreg”) install.packages(”scales”)
# balíček pro časové řady # balíček pro časové řady # balíček pro testy normality # balíček pro lepší grafy # balíček pro sázení tabulek do LaTEXu # balíček pro Likelihood ratio test # balíček pro modely typu ARCH # balíček pro funkci melt na zpracovaní dat # balíček pro funkci quantile # balíček pro vykreslovaní v ggplot2 funkce geom quantile
# Načtení použitých knihoven library(xtable) # knihovna pro sázení tabulek do Latexu library(fGarch) # knihovna pro modelovaní volatility library(lmtest) # knihovna pro Likelihood ratio test library(FinTS) # knihovna s nástroji pro finanční časové řady library(nortest) # knihovna pro testy normality library(tseries) # knihovna pro časové řady library(ggplot2) # knihovna pro grafický nástroj library(reshape) library(quantreg) library(scales) library(relimp, pos=4) require(graphics) require(grid)
# Nastavení adresáře s daty getwd() # zjištění aktualního pracovního adresáře setwd(”C:/Users/pakosta/Documents/R”) # nastavení pracovního adresáře # Nastavení cest k použitým funkcím filepath <− ”.//functions//” # nastavení cesty k funkcím, # složka s funkcemi musí být v aktualním pracovním adresáři source(paste(filepath, source(paste(filepath, source(paste(filepath, source(paste(filepath, source(paste(filepath,
”QQplot.r”, sep = ””)) ”qpacf.r”, sep = ””)) ”qacf.r”, sep = ””)) ”NKH.r”, sep = ””)) ”NKHst.r”, sep = ””))
# Příprava dat # Načtení dat data <− read.csv(file = ”gold.csv”, head = FALSE, sep = ”;”)
# Rozdělení dat do skupin # Převod formatu datumů na format datumů v R dates = rev(as.Date(data$V1, format = ”%d.%m.%Y”)) serie = rev(data$V2) # uložení dat do pomocné proměnné
Jak je vidět z příkazů první části skriptu, data jsou načtena ze souboru .csv, který je také na CD přiloženém k této práci. Soubor obsahuje 2 sloupce s datumy (formát DD.MM.RRRR) a cenami zlata v daných časech (s desetinnou tečkou). Dále je pak nastavena cesta do vlastního pracovního adresáře, který je zvolen pomocí příkazu setwd a také obsahuje složku s dalšími funkcemi, které hlavní skript využívá.
— 24 —
5: Project R
5.2
Použité nástroje a funkce
V této části budou uvedeny některé důležité balíčky a nejdůležitější použité funkce. K těmto vybraným funkcím budou vždy stručně napsané základní informace. Pokud požadujeme detailnější informace o balíčku nebo funkci, v příkazovém řádku programu napíšeme před název balíčku nebo funkce otazník (např. ?F inT S nebo ?shapiro.test). Chybějící balíčky lze doinstalovat pomocí již uvedené funkce install.packages(”balicek”). Balíček „tseriesÿ - pro analýzu časových řad a výpočetní finančnictví. Podrobné informace dostupné na: http://cran.r-project.org/web/packages/tseries/tseries.pdf Balíček „FinTSÿ - R-kový protějšek knihy Analysis of Financial Time Series, 2nd ed. (Wiley), od autora Tsay (2005). Zahrnuje soubory dat, funkce a soubory skriptů potřebné k práci s některými příklady. Podrobné informace dostupné na: http://math.uncc.edu/~zcai/FinTS.pdf Balíček „nortestÿ - soubor 5 testů pro složené hypotézy normality. Podrobné informace dostupné na: http://cran.r-project.org/web/packages/nortest/nortest.pdf Balíček „ggplot2ÿ - je kreslící systém pro R, který používá jen dobré vlastnosti ze základní R-kové grafiky a balíčku Lattice grafiky. Tento balíček poskytuje silný model pro grafiku, který umožňuje lehce vytvářet komplexní vícevrstvou grafiku. Podrobné informace dostupné na: http://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf Balíček „fGarchÿ - Prostředí pro výuku ”Finanční inženýrství a výpočetní finance”. Podrobné informace dostupné na: http://cran.r-project.org/web/packages/fGarch/fGarch.pdf Další použité balíčky lze dohledat na stránkách: http://cran.r-project.org/. V následující části budou uvedeny funkce vytvořené v programu R, které jsou použity v hlavním výpočtovém skriptu a které již v dalších částech budou jen volány. V aplikační části této práce budou uváděny další R kódy s příkazy a funkcemi, které byly využity k řešení, výpočtu nebo sestavení grafu. Pokud budou pro další výsledky použity stejné příkazy nebo funkce, které již byly uvedeny a jsou tedy použity pouze pro jiná data, pak už R kód nebude uváděn. Celý zdrojový kód pro získání všech výsledků lze nalézt na CD, které je přílohou této práce. Funkce pro vytvoření Q-Q grafu: R code: 1 # Q−Q graf 2 QQplot <− function(y){ 3 x <− qnorm(ppoints(y)) 4 d <− data.frame(x = x, y = sort(y)) 5 xx <− qnorm(c(0.25, 0.75)) 6 yy <− quantile(y, c(0.25, 0.75)) 7 slope <− diff(yy) / diff(xx) 8 int <− yy[1L] − slope ∗ xx[1L] 9 p <− ggplot(d, aes(x = x, y = y)) + geom point() + geom hline(yintercept = 0) 10 + xlab(’Theoretical Quantiles’) + ylab(’Sample Quantiles’) 11 + geom abline(intercept = int, slope = slope, color = ’blue’) 12 return(p) 13 }
— 25 —
5: Project R
Funkce pro vytvoření ACF grafu: R code: 1 # ACF 2 qacf <− function(x, conf.level = 0.95, max.lag = NULL, min.lag = 0, title = ””) { 3 ciline <− qnorm((1 − conf.level) / 2) / sqrt(length(x)) 4 bacf <− acf(x, plot = FALSE, lag.max = max.lag, type = ”correlation”) 5 bacfdf <− with(bacf, data.frame(lag, acf)) 6 if (min.lag > 0) { bacfdf <− bacfdf[−seq(1, min.lag), ] } 7 significant <− (abs(bacfdf[, 2]) > abs(ciline))ˆ2 8 bacfdf <− cbind(bacfdf, significant) 9 q <− qplot(lag, acf, data = bacfdf, geom = ”bar”, stat = ”identity”, 10 position = ”identity”, ylab = ”Autocorrelation”, main = title, fill = factor(significant)) 11 + geom hline(yintercept = −ciline, color = ”blue”, size = 0.2) 12 + geom hline(yintercept = ciline, color = ”blue”, size = 0.2) 13 + geom hline(yintercept = 0, color = ”red”, size = 0.3) 14 + scale fill hue(name = paste(”Significant at the\n”, conf.level, ”level”), 15 breaks = 0:1, labels = c(”False”, ”True”)) 16 + theme(legend.justification = c(1,1), legend.position = c(1,1)) 17 return(q) 18 }
Funkce pro vytvoření PACF grafu: R code: 1 # PACF 2 qacf <− function(x, conf.level = 0.95, max.lag = NULL, min.lag = 0, title = ””) { 3 ciline <− qnorm((1 − conf.level) / 2) / sqrt(length(x)) 4 bacf <− acf(x, plot = FALSE, lag.max = max.lag, type = ”correlation”) 5 bacfdf <− with(bacf, data.frame(lag, acf)) 6 if (min.lag > 0) { bacfdf <− bacfdf[−seq(1, min.lag), ] } 7 significant <− (abs(bacfdf[, 2]) > abs(ciline))ˆ2 8 bacfdf <− cbind(bacfdf, significant) 9 q <− qplot(lag, acf, data = bacfdf, geom = ”bar”, stat = ”identity”, 10 position = ”identity”, ylab = ”Autocorrelation”, main = title, fill = factor(significant)) 11 q <− q + geom hline(yintercept = −ciline, color = ”blue”, size = 0.2) 12 + geom hline(yintercept = ciline, color = ”blue”, size = 0.2) 13 + geom hline(yintercept = 0, color = ”red”, size = 0.3) 14 + scale fill hue(name = paste(”Significant at the\n”, conf.level, ”level”), breaks = 0:1, labels = c(”False”, ”True”)) 15 + theme(legend.justification = c(1,1), legend.position = c(1,1)) 16 return(q) 17 }
Funkce pro vytvoření grafu histogramu, jádrového odhadu funkce hustoty a hustoty odhadnutého teoretického normálního rozdělení: R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
NKH <− function(x, title = ””) {
# Jádrový odhad distribuční funkce kernel = density(x, bw = ”nrd0”, adjust = 1, kernel = c(”gaussian”), weights = NULL, window = kernel, width, give.Rkern = FALSE, cut = 3, na.rm = FALSE) xfit<−seq(min(x),max(x),length=length(kernel$y)) yfit<−dnorm(xfit, mean = mean(x), sd = sd(x)) df <− data.frame(type = factor( rep(c(”kernel (gaussian)”, ”Teoretical normal”), each=length(kernel$x))), LogP = c(kernel$y, yfit), axisX = c(kernel$x, kernel$x))
# vykresleni jadroveho odhadu pomoci ggplot2 H = ggplot(data.frame(x = x), aes(x = x)) H = H + geom histogram(aes(y = . .density. .,), colour=I(’#FFFFFF’), binwidth = density(x)$bw) + xlab(’Data(LogProfit)’) + ylab(’Density’) + geom line(data = df, aes(x=axisX, y=LogP, colour=type), size = 1) + geom polygon(data = df, aes(x=axisX, y=LogP, fill=type ,group=type), alpha = 0.2) + theme(legend.justification = c(1,1), legend.position = c(1,1)) return(H) }
V práci je také použita vytvořená funkce pro vykreslení histogramu, jádrového odhadu funkce hustoty a hus— 26 —
5: Project R
toty odhadnutého teoretického studentovo rozdělený. Tato funkce je podobná jako předcházející, ale je modifikovaná právě pro studentovo rozdělení: R code: 1 NKHst <− function(x, t1, t2, t3, title = ””) { 2 # Jádrový odhad distribuční funkce 3 kernel = density(x, bw = ”nrd0”, adjust = 1, kernel = c(”gaussian”), weights = NULL, window = kernel, 4 give.Rkern = FALSE, cut = 3, na.rm = FALSE) 5 6 xfit<−seq(min(x),max(x),length=length(kernel$y)) 7 yfit<−dstd(xfit,t1,t2,t3) 8 9 df <− data.frame(type = factor( rep(c(”kernel (gaussian)”, ”Teoretical student”), each=length(kernel$x))), 10 LogP = c(kernel$y,yfit), axisX = c(kernel$x, kernel$x)) 11 12 H = ggplot(data.frame(x = x), aes(x = x)) 13 H = H + geom histogram(aes(y = . .density. .,),colour=I(’#FFFFFF’), binwidth = density(x)$bw) 14 H = H + xlab(’Data(LogProfit)’) + ylab(’Density’) 15 H = H + geom line(data = df, aes(x = axisX, y = LogP, colour = type), size = 1) 16 H = H + geom polygon(data = df, aes(x = axisX, y = LogP, fill = type, group = type), alpha = 0.2) 17 H = H + theme(legend.justification=c(1,1), legend.position=c(1,1)) # pozice legendy 18 19 return(H) 20 }
— 27 —
6: Aplikační část
6
Aplikační část
Doposud probranou teorii lze aplikovat na vysokofrekvenční časové řady, kterými jsou např. devizové kurzy, komodity, akcie a jiné. Aplikační část této práce bude zaměřena na zpracování vývoje komodity zlata26 . Investiční zlato je kupované jako výhodná investice na ochranu proti znehodnocení úspor. Obchoduje se zejména jako zlaté slitky, cihly nebo mince s vysokou ryzostí kovu (až 999/1000), přičemž v ČR je dodávka investičního zlata osvobozena od DPH27 . Investiční zlato se obchoduje ve standardizovaných hmotnostech, jednou z nich je například 1 oz (trojská unce = 31,103 g). Zlato patří k jedné z nejatraktivnějších investic, protože pro mnohé investory má nemalý význam, jako diverzifikační nástroj či jako uchovatel hodnoty v dobách nejistoty nebo při očekávání vysoké inflace, kterou by zlato mělo zahrnout ve své ceně. V [15] je zpracována analýza této řady a možné vlivy na vývoj ceny zlata, které se však ukázaly jako nevýznamné. Data ceny zlata jsou stažena z [11] a obsahují údaje od 9.2.2005 do 31.10.2012. Zahrnují tak celkem 1858 pozorování s denní frekvencí, ve smyslu čistě obchodních dnů. Kvůli existenci neobchodních dnů28 nemohou být časové okamžiky pořízení jednotlivých pozorování stejně vzdáleny a proto bude vždy existovat nakumulovaná informace po průběhu neobchodních dnů, která by teoreticky měla ovlivnit počátek obchodování zvýšenou volatilitou. Ve skutečnosti ovšem není zvykem zveřejňovat kurzotvorné informace v neobchodních dnech a tento efekt by tak měl být potlačen [16].
Obr. 3: Vývoj časové řady zlata za období 9.2.2005-31.10.2012 26 Informace
o investičním zlatě: http://wiki.kurzy.cz/Investiční_zlato/. v §92 „Zvláštní režim pro investniční zlatoÿ zákona o DPH č. 235/2004 Sb. http://zakony.kurzy.cz/ 235-2004-zakon-o-dani-z-pridane-hodnoty-dph/paragraf-92/. 28 Svátky či víkendy. 27 Podrobnosti
— 28 —
6: Aplikační část
R code: 1 # Vykreslení časové řady 2 qplot(dates, serie, geom = ”line”, xlab = ”Time”, ylab = ”Price of gold”, 3 main = ”Time serie of GOLD (1 USD/troy once)”, ylim = c(0,2000)) 4 + geom point(colour = ”black”, size = 2, alpha = 0.3)
6.1
Základní popis a předzpracování dat
V této části budou popsány základní charakteristiky zpracovávaných dat ceny zlata. Rok 2005 2006 2007 2008 2009 2010 2011 2012 2005-2012
Min. 412.9 528.0 607.2 712.3 812.7 1063.0 1312.0 1539.0 412.9
25% Kvantil 426.8 571.4 658.8 828.0 921.3 1139.0 1436.0 1613.0 657.8
Medián 438.2 612.0 677.1 884.9 949.4 1217.0 1546.0 1655.0 933.2
Průměr 446.5 606.5 700.6 872.3 975.1 1229.0 1574.0 1663.0 1025.0
75% Kvantil 463.0 634.6 743.2 921.8 1016.0 1319.0 1699.0 1721.0 1382.0
Max. 527.9 721.4 837.2 1008.0 1216.0 1426.0 1903.0 1792.0 1903.0
Počet 210 220 208 237 258 254 257 214 1858
Sm. odchylka 25.921 40.336 59.356 68.617 86.073 99.088 149.915 65.962 418.524
Tab. 1: Základní údaje časové řady cen zlata Časová řada ceny zlata je nestacionární, má totiž určitý trendový vývoj v čase a jak je z obrázku 4 vidět, mohlo by se jednat o exponenciální trend. Pomocí autokorelační funkce (ACF), resp. parciální autokorelační funkce (PACF)29 , lze potvrdit, že časová řada je nestacionární.
Obr. 4: Vývoj časové řady zlata a odhad exponenciálního trendu 29 Podle
předchozí kapitoly 2.
— 29 —
6: Aplikační část
R code: 1 # Získání základních statistických charakteristik 2 summary(serie) # informace o datech 3 sqrt(var(serie)) # informace o směrodatné odchylce 4 length(serie) # informace o délce řady 5 6 # Řada cen zlata + exponenciální trend 7 ts = serie 8 tr <− seq(from = 1, to = length(ts), by = 1) #cas 9 yt <− ts ∗ log(ts) 10 xt1 <− ts 11 xt2 <− tr ∗ ts 12 summary(a<−lm(formula = yt˜xt1 + xt2 −1)) 13 exp(a$coef[1]) 14 exp(a$coef[2]) 15 expT = exp(a$coef[1])∗(exp(a$coef[2])ˆtr) 16 17 df <− data.frame(x = dates, trend = expT) 18 G = ggplot(data.frame(Time = dates, Gold = serie), aes(Time,Gold), ylim = c(0,2000)) 19 + geom line(colour = ”black”, size = 1) 20 + geom point(colour = ”black”, size = 2, alpha = 0.3) 21 + geom line(aes(dates, trend), data = df, colour = ”red”, size = 1) 22 print(G)
Charakter ACF na obrázku 5 podporuje tvrzení, že řada cen zlata je nestacionární. Graf ACF ukazuje velmi pozvolně klesající průběh autokorelační funkce, resp. koeficentů korelace, které přesahují hranici statistické významnosti30 .
Obr. 5: Graf ACF pro řadu cen zlata
30 Pravděpodobností
√ √ 95% interval (−1.96/ n, 1.96/ n) uvedený v kapitole 2.
— 30 —
6: Aplikační část
Obr. 6: Graf PACF pro řadu cen zlata R code: 1 #ACF a PACF původních dat 2 ACF1 = acf(serie, lag.max = NULL, type = ”correlation”, 3 plot = FALSE, na.action = na.fail) 4 PACF1 = acf(serie, lag.max = NULL, type = ”partial”, 5 plot = FALSE, na.action = na.fail) 6 LAGACF1 = length(ACF1$lag) # nastavení délky zpoždění 7 LAGPACF1 = length(PACF1$lag) # nastavení délky zpoždění 8 9 # vykreslení ACF pomocí ggplot2 10 source(’.//Functions//qacf.r’) 11 G = qacf(serie, conf.level = 0.95, max.lag = LAGACF1, 12 title = ”Graph of ACF for price of gold”) 13 print(G) 14 15 # vykreslení PACF pomocí ggplot2 16 source(’.//Functions//qpacf.r’) 17 G = qpacf(serie, conf.level = 0.95, max.lag = LAGPACF1, 18 title = ”Graph of PACF for price of gold”) 19 print(G)
— 31 —
6: Aplikační část
Nestacionární finanční časovou řadu ceny zlata lze podle vzorce31 rt = ln Pt − ln Pt−1 = pt − pt−1 transformovat na stacionární. Získáme tak řadu logaritmických výnosů, která má v tomto případě intuitivnější význam než jednoduchá diferencovaná řada. Nejprve tedy logaritmováním linearizujeme předpokládaný exponenciální trend a dostáváme tak řadu logaritmů cen zlata, která v čase vypadá následovně.
Obr. 7: Vývoj logaritmů cen zlata a odhad lineárního trendu
R code: 1 # Logaritmovaní původní časové řady 2 logSerie = log(serie) # logaritmy dat 3 4 # Graf logaritmů cen zlata + linearní trend 5 G = ggplot(data.frame(Time = dates, serie = logSerie), aes(Time, serie)) 6 + geom point(colour = ”black”, size = 2, alpha = 0.3) 7 + geom line() 8 + geom smooth(method = ”lm”, se = FALSE, colour = ”red”, size = 1) 9 + scale fill discrete(””) + ylim(c(5,8)) + ylab(”Serie log(Price of gold)”) 10 print(G)
31 Uvedený
v kapitole 2.2.
— 32 —
6: Aplikační část
Po logaritmické transformaci lze diferencí 1. řádu získat logaritmy výnosů původní časové řady, které v čase vypadají takto:
Obr. 8: Vývoj časové řady logaritmů výnosů ceny zlata
R code: 1 # Dokončení převodu na stacionarní časovou řadu (diference 1. řádu) 2 lags = 1 # nastavení řádu diference 3 LogProfit = diff(logSerie, lag=lags) # výnosy logaritmů dat 4 5 # Graf logaritmů výnosů 6 G = ggplot(data.frame(Time=dates[2:length(dates)],LogProfit), aes(Time,LogProfit)) 7 + geom line(colour = ”red”, size = 1) 8 + geom point(colour = ”black”, size = 1) 9 print(G)
— 33 —
6: Aplikační část
Obr. 9: Graf ACF logaritmů výnosů časové řady ceny zlata
Obr. 10: Graf PACF logaritmů výnosů časové řady ceny zlata Řada logaritmů výnosů zlata (dále už jen výnosů) můžeme již považovat za stacionární, což je patrné z přecházejících obrázků 8, 9, 10. Autokorelaci již neprokazuje ani sestavená ACF a PACF, kde je řada výnosů charakteristická nesystematickými výkyvy kolem nulové hodnoty. Zároveň můžeme na obrázku 8 pozorovat proměnlivou volatilitu (variabilitu), která naznačuje měnící se rozptyl. — 34 —
6: Aplikační část
6.2
Předpoklad normality
Modelování časových řad je často svazováno s předpokladem, že data (resp. residua dat) pocházejí z normálního rozdělení (nebo alespoň se k němu blíží). Rok 2005 2006 2007 2008 2009 2010 2011 2012 2005-2012
Počet 209 220 208 237 258 254 257 214 1857
Průměr 0.00109 0.00095 0.00128 0.00024 0.00085 0.00102 0.00037 0.00045 0.00077
Sm. odchylka 0.00836 0.01624 0.01201 0.02001 0.01266 0.01023 0.01264 0.00967 0.01324
Šikmost -0.15630 -0.56442 -1.10897 0.04469 -0.13740 -0.56905 -0.52707 -0.41322 -0.35826
Špičatost 1.54873 2.38897 4.68733 1.83305 0.51992 2.11572 1.67968 4.83308 3.76726
Min. -0.03036 -0.07450 -0.06414 -0.06683 -0.04202 -0.04424 -0.05448 -0.05293 -0.07450
Max. 0.02695 0.05535 0.03391 0.07355 0.04166 0.03167 0.03555 0.03802 0.07360
Tab. 2: Statistické charakteristiky řady výnosů Liliefors test Shapiro-Wilk test Jarqu-Bera test
D = 0.076 W = 0.9528 χ2 = 1142.172
p-value< 2.2 · 10−16 p-value< 2.2 · 10−16 p-value< 2.2 · 10−16
Tab. 3: Výsledky testů normality pro řadu výnosů Tabulka 3 ukazuje testy normality popsané v kapitole 2.2, které neprokazují normalitu výnosů. Odlišnost od normálního rozdělení ukazuje také tabulka 2, kde šikmost a špičatost vychází odlišně od normálního rozdělení (a to i v dílčích rocíh), které má šikmost rovnu nule a špičatost rovnu třem.
Obr. 11: Histogram výnosů, gaussovský jádrový odhad hustoty rozdělení a teoretické normální rozdělení32 32 Základní informace o odhadech jádrových hustot lze nalézt např. v disertační práci Mgr. Martina Řezáče, Ph.D., dostupné z: http://is.muni.cz/th/20411/prif_d/disertace.pdf nebo také v diplomové práci Karolíny Mladé, dostupné z: http://is.muni. cz/th/270467/prif_m/diplomkaMlada.pdf.
— 35 —
6: Aplikační část
Na předchozím obrázku je pak histogram výnosů a odhadnuté rozdělení výnosů pomocí gaussovského jádrového odhadu a také teoretické normální rozdělení s parametry odhadnutými z dat výnosů. Z tohoto porovnání jsou patrné vlastnosti finančních časových řad popsané v kapitole 2.2. Histogram33 na obrázku 11 ukazuje větší špičatost oproti normálnímu rozdělení a těžké konce34 , které ukazuje také Q-Q graf na následujícím obrázku 12. Tyto obrázky a tabulky demonstrují, jak se data výnosů liší od normálního rozdělení s odhadnutými parametry střední hodnoty a rozptylu z dat výnosů, ale také indikují použití modelů volatility.
Obr. 12: Q-Q graf výnosů
R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
# Testy normality x = LogProfit TSstat = FinTS.stats(x); TSstat LLT = lillie.test(x); LLT SWT = shapiro.test(x); SWT JBT = jarque.bera.test(x); JBT
# Grafické testy normality logaritmů výnosů (LogProfit) # Jadrový odhad distribuční funkce source(’.//functions//NKH.r’) G = NKH(LogProfit) print(G)
# QQ graf logaritmů výnosů source(’.//functions//QQplot.r’) G = QQplot(LogProfit) print(G)
33 Četnosti histogramu vyplývají z funkce v R. Zde je implementováno rule-of-thumb (pravidlo palce) pro volbu šířky pásma gaussovského jádrového odhadu hustoty. (Silvermanovo ”pravidlo”, Silverman (1986, strana 48, kapitola (3.31))). 34 Nebo také tlusté konce, anglicky: heavy tails, fat tails.
— 36 —
6: Aplikační část
6.3
Fáze 1: Odhad úrovňového modelu
V dalších částech již začíná samotné modelování časové řady, které se bude více méně řídit podle kapitoly 3.3. Nejprve se pokusíme z dat výnosů odstranit lineární závislosti pomocí tzv. „úrovňového modeluÿ. Z předcházejících obrázků lze vidět, že úroveň časové řady se v čase mění a proto lze použít jeden z lineárních modelů ARMA(p,q)35 , které jsou založeny na podmínce, že se sice podmíněná střední hodnota mění v čase, podmíněný rozptyl je však konstantní, což neodpovídá realitě. Tento problém lze vyřešit právě uvažovanými modely volatility v další části. Podstatným rysem této koncepce je, že se nemění původní požadavek normality. Opodstatnění použití ARMA modelu vyplývá již z obrázku 10 na straně 34, kde PACF výnosů ukazuje porušení nekorelovanosti už při k = 1 (dále pak k = 21, 23, 27). Tento fakt lze ještě ověřit například pomocí Portmentau testu. Řada výnosů obsahuje 1857 hodnot a přirozený logaritmus tohoto počtu je 7.5267. Z toho plyne použití Portmentau testu pro k = 8 korelačních koeficientů. Výsledkem testu je p-hodnota rovna 0.04551 a to je na hranici zamítnutí (na hladině významnosti α = 5%). Hypotézu o nekorelovanosti tedy zamítáme a pomocí ARMA modelu se pokusíme závislosti odstranit. R code: 1 # Test autokorelace výnosů (LogProfit) 2 x = LogProfit 3 # Nastavení hodnoty zpožděni podle vzorce log(n) 4 k = round(log(length(x))) 5 # Portmentau test − hypotéza o nekorelovanosti 6 Q = Box.test(x, lag = k, type = ’Ljung’); Q
Řád modelu ARMA lze určit několika způsoby. Jedním z nich je určení řádu podle tvaru ACF a PACF dat výnosů. Možné tvary ACF a PACF a k nim určené řády modelu jsou36 : • ARMA(1,0): ACF exponenciálně klesá k nule, PACF má pouze jeden vrchol na pozici 1. • ARMA(2,0): ACF má sinusový průběh nebo několik exponenciálních poklesů, PACF má vrcholy na pozicích 1 a 2. Na dalších pozicích nejsou již korelace. • ARMA(0,1): ACF má pouze jeden vrchol na pozici 1, PACF jde exponenciálně k nule. • ARMA(0,2): ACF má dva vrcholy na pozici 1 a 2, bez dalších korelací, PACF má sinusový průběh, nebo exponenciálně tlumený průběh. • ARMA(1,1): ACF má exponenciální pokles začínající na pozici 1, PACF má oscilující pokles začínající na pozici 1. • ARMA(1,1): ACF má oscilující pokles začínající na pozici 1, PACF má exponenciální pokles začínající na pozici 1. • ARMA(p,q): ACF přímý nebo oscilující pokles k nule začínající na pozici q. PACF přímý nebo oscilující pokles k nule začínající na pozici p. Pokud se podíváme na ACF a PACF na obrázku 9 a 10 (na straně 34), lze vidět, že tvary těchto funkcí odpovídají spíše modelu ARMA(1,1). 35 Často
je využíván pouze jednoduchý model Xt = µ + εt , zvláště pak u finančních časových řad. zde: www.statsoft.cz/file1/PDF/newsletter/2012_07_23_ARMA.pdf.
36 Například
— 37 —
6: Aplikační část
V následující tabulce budou shrnuty výsledky odhadnutých ARMA(p,q) modelů. Požadujeme co nejjednodušší model, proto bude využit model maximálně do řádu ARMA(1,1). Kritéria lnL(ϕ), ˆ AIC a BIC jsou spočtena podle kapitoly 3.7.1. model ARMA(1,0) ARMA(0,1) ARMA(1,1)
lnL(ϕ) ˆ 5400.62 5400.78 5400.87
AIC -10795.23 -10795.97 -10793.74
BIC -10778.65 -10778.97 -10771.64
Tab. 4: Tabulka kritérií pro různé modely typu ARMA Podle AIC kritéria lze volit model ARMA(0,1), který vede k nejnižší hodnotě AIC a stejně tak podle BIC kritéria. Volba konečného modelu bude však z jiného pohledu. Protože se kritéria jednotlivých modelů nijak významně od sebe neliší, můžeme využít model ARMA(1,0), kvůli jeho snadné interpretovatelnosti. Obecný model ARMA(1,0) tedy vypadá takto: Xt = µ + φXt−1 + εt , kde Xt jsou výnosy, |φ| < 1 a εt ∼ W N (0, σε2 ) je proces bílého šumu se střední hodnotou 0 a rozptylem σε2 . Odhady parametrů modelu ARMA(1,0) jsou37 :
µ φ
Estimate -0.0007192 0.0608666
Std. Error 0.0003071 0.0231723
t-value -2.342 2.627
P r(> |t|) 0.01920 * 0.00862 **
Tab. 5: Tabulka odhadů parametrů pro ARMA(1,0) model V předcházející tabulce si můžeme všimnou, že parametry sice jsou statisticky významné, ale dalo by se zde uvažovat o přímém použití modelů typu ARCH, díky nepříliš vysoké významnosti parametrů u tohoto konkrétního modelu ARMA(1,0). V tomto případě bylo však rozhodnuto o ponechání modelu ARMA(1,0), který by mohl mít svá určitá ekonomická opodstatnění, která však v této práci nebudeme rozebírat ani podrobně zkoumat. 37 Významnost parametrů: ’***’ pro interval h0, 0.001i, ’**’ pro interval (0.001, 0.01i, ’*’ pro inteval (0.01, 0.05i, ’.’ pro interval (0.05, 0.1i, ’ ’ pro interval (0.1, 1i.
— 38 —
6: Aplikační část
Obr. 13: Řada výnosů a odhad modelu ARMA(1,0) R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
#ARMA model x = LogProfit p=1 q=0 ARMAm = arma(x, order = c(p, q), include.intercept = TRUE) summary(ARMAm)
# M = počet parametrů+1, tj. konstanta+p+q+1 M=1+p+q+1
# AIC kriterium ARMAmAIC = summary(ARMAm10)$aic
# Max. log−likelihood kriterium ARMAmLLH = (summary(ARMAm10)$aic − 2 ∗ M) / (−2)
# BIC kriterium ARMAmBIC = ((summary(ARMAm10)$aic − 2∗M)) + M ∗ log(length(x))
# Příprava dat pro graf df <− data.frame(dates = dates[(max(p,q) + 2):length(dates)], LogP = LogProfit[(max(p,q) + 1 + lags):length(dates) − lags], Res = ARMAm$residuals[(max(p,q) + 1 + lags):length(dates) − lags], arma = ARMAm$fitted.value[(max(p,q) + 1 + lags):length(dates) − lags])
# Graf ARMA modelu G = ggplot(df, aes(x = dates)) + xlab(’Time’) + ylab(’LogProfit − ARMA(1,0)’) + geom point(aes(y = LogP), size = 2, colour = ”#999999”) + geom line(aes(y = arma), size = 1, colour = ”#E69F00”) print(G)
— 39 —
6: Aplikační část
6.3.1
Test residuí
Po odhadu úrovňového modelu následuje test residuí εˆt , které by v případě, že v modelu není objevena podmíněná heteroskedasticita, měly vést k procesu gaussovského bílého šumu. Následující údaje výsledků testů normality, autokorelace a homoskedasticity pro řadu {ˆ εt } potvrzují fakt zmíněný na začátku kapitoly 6.3, že tyto modely neuvažují nekonstantnost rozptylu. Díky tomu je vidět porušení předpokladu gaussovského bílého šumu residuí ARMA modelu. Liliefors test Shapiro-Wilk test Jarqu-Bera test
D = 0.0734 W = 0.9545 χ2 (2) = 1086.963
p-value< 2.2 · 10−16 p-value< 2.2 · 10−16 p-value< 2.2 · 10−16
Tab. 6: Výsledky testů normality pro řadu {εt } Testy normality residuí ARMA(1,0) modelu zamítají hypotézu o normalitě dat (tabulka 6). Na grafu 14 je vidět, že residua mají rozdělení špičatější a také s těžšími konci. Při porovnání následujících obrázků ACF a PACF residuí modelu ARMA(1,0) a obrázků ACF a PACF logaritmů výnosů z kapitoly 6.1 (obrázek 9 a 10), můžeme vidět, že některé závislosti se podařilo odstranit. Pro ověření můžeme opět použít Portmanteau test z kapitoly 2.2.2. V tomto případě však s modifikací - testová statistika Q(k) v případě modelu ARMA(p,q) konverguje v distribuci k χ2 -rozdělení o k − (p + q) stupních volnost, tj. χ2m−(p+q) -rozdělení. Tím se zahrne odhadovaný počet parametrů [2]. Podle Portmanteauova testu (k = 8 − (1 + 0) = 7), který dává p-hodnotu 0.2295, tak lze přijmout hypotézu o nekorelovanosti residuí ARMA(1,0) modelu.
Obr. 14: Histogram residuí, gaussovský jádrový odhad hustoty rozdělení a teoretické normální rozdělení
— 40 —
6: Aplikační část
Obr. 15: ACF residuí ARMA(1,0) modelu
Obr. 16: PACF residuí ARMA(1,0) modelu Posledním krokem je pak test podmíněné heteroskedasticity38 pomocí ARCH-LM testu (popsaný v kapitole 3.4), který vrátí téměř nulovou p-hodnotu < 2.2 · 10−16 (χ2 (8) = 148.54). Zamítáme tedy hypotézu o nepřítomnosti podmíněné heteroskedasticity. 38 Testování
tzv. ARCH efektu.
— 41 —
6: Aplikační část
R code: 1 # Test podmíněné heteroskedasticity výnosů (LogProfit) 2 x = LogProfit 3 # Nastavení hodnoty zpoždění podle vzorce log(n) 4 k = round(log(length(x))) 5 # Arch−LM test 6 AT = ArchTest(x, lags = k); AT
6.4
Fáze 2: Odhad modelu volatility
Stejně jako v předcházející části, kde byla odstraněna měnící se úroveň časové řady výnosů, zde se budeme zabývat odstraněním variability z residuí ARMA modelu, která se také v čase mění. V této fázi je tedy potřeba zvolit model volatility. Jako první se nabízí využítí modelu ARCH(q). Použití tohoto modelu však často vede k výběru modelu s vysokým řádem q. Proto lze využít modelu volatility typu GARCH(p,q). V praxi bývají často využity tyto modely s nízkým řádem. Z tohoto důvodu bude použit model GARCH nejvýše do řádu GARCH(2,2). K určení řádu modelu GARCH můžeme využít ACF a PACF kvadrátů residuí ARMA(1,0) modelu, které jsou na následujících obrázcích.
Obr. 17: ACF kvadrátů residuí ARMA(1,0) modelu
— 42 —
6: Aplikační část
Obr. 18: PACF kvadrátů residuí ARMA(1,0) modelu Předcházející obrázky naznačují řád modelu GARCH(1,1), avšak řád můžeme určit přesněji pomocí hodnoty Akaikeho (AIC) a Bayesova informačního kritéria (BIC) a hodnoty maximalizované logaritmické věrohodnostní funkce. Tyto hodnoty jsou uvedeny v následující tabulce pro různá p a q modelu GARCH.
Ln L(ϕ) ˆ AIC BIC
GARCH(1,1) 5558.165 −5.98617 −5.97724
GARCH(2,1) 5559.668 −5.98671 −5.97480
GARCH(1,2) 5558.052 −5.98497 −5.97306
GARCH(2,2) 5559.992 −5.98598 −5.97110
Tab. 7: Hodnoty maximalizované věrohodnostní funkce, AIC a BIC modelů GARCH Podle kritéria BIC v předcházející tabulce volíme model GARCH(1,1). Podle AIC bychom volili model GARCH(2,1) a pokud budeme dále zkoumat tento model, zjistíme, že v tomto případě má model parametry α1 a α2 statisticky nevýznamné (na hladině významnosti α = 5%, resp. α = 1%).
ω α1 α2 β1
Odhad 1.960 · 10−06 1.883 · 10−02 3.583 · 10−02 9.344 · 10−01
Sm. odchylka 7.076 · 10−07 1.484 · 10−02 1.799 · 10−02 1.145 · 10−02
t-stat. 2.770 1.269 1.992 81.579
p-hodnota 0.00561 ** 0.20443 0.04637 * < 2 · 10−16 ***
Tab. 8: Odhadnuté parametry modelu GARCH(2,1) To nám napovídá k použití modelu GARCH(1,1), jak už ukazuje kritérium BIC, které má větší vypovídací schopnost než AIC podle kapitoly 3.7.1. Volíme tedy model GARCH(1,1) jako finální model volatility. Kritéria však neříkají nic o tom, jaké předpoklady splňují residua modelu nebo jestli jsou parametry modelu statisticky významné.
— 43 —
6: Aplikační část
odhad 1.684 · 10−06 4.792 · 10−02 9.426 · 10−01
ω α1 β1
sm. odch. 5.811 · 10−07 7.197 · 10−03 8.645 · 10−03
t-stat. 2.897 6.658 109.027
p-hodnota 0.00376 ** 2.77 · 10−11 *** < 2 · 10−16 ***
Tab. 9: Odhadnuté parametry modelu GARCH(1,1) Podle tabulky 9 jsou parametry v pořádku, všechny jsou statisticky významné a pro model mají smysl. R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14
# Odhad parametrů GARCH GARCHm = garchFit(˜ garch(1,1), data = ARMAe, include.mean = FALSE, cond.dist = ”norm”)
# Hodnota maximalizované logaritmické věrohodnostní funkce GARCHm@fit$llh
# Hodnota Akaikeho (AIC) informačního kriteria GARCHm@fit$ics[1]
# Hodnota Bayesova (BIC) informačního kriteria GARCHm@fit$ics[2] summary(GARCHm)
# Všechny informace (testy residuí, odhady parametrů,. . .) jsou obsaženy v: GARCHm@fit summary(GARCHm)
6.4.1
Test residuí a úprava modelu
p ˆ t , pro které by Pro ověření vhodnosti modelu GARCH(1,1) lze využít standardizovaných reziduí ˆt = εˆt / h mělo platit ˆt ∼ N (0, 1). Následující tabulka ukazuje, že model GARCH(1,1) uspěl při modelování volatility, resp. testy ukázaly, že v řadě residuí ˆt již není podmíněná heteroskedasticita prokázána.
Jarque-Bera Test Shapiro-Wilk Test Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test LM Arch Test
ˆ ˆ ˆ ˆ ˆ ˆ2 ˆ2 ˆ2 ˆ
2
χ W Q(10) Q(15) Q(20) Q(10) Q(15) Q(20) T R2
Statistic 448.17190 0.97587 10.79942 12.24231 15.59231 4.82837 8.87266 9.45254 6.09475
p-Value 0 0 0.37336 0.66060 0.74157 0.90234 0.88408 0.97705 0.91123
Tab. 10: Testy standardizovaných residuí {ˆ } modelu GARCH(1,1) Testy normality residuí (Jarque-Bera test, Shapiro-Wilk test popsaný v kapitole 2.2) však ukazují nesplnění podmínky normality. Tento fakt můžeme ilustrovat také na následujícím grafu histogramu residuí ˆt .
— 44 —
6: Aplikační část
Obr. 19: Histogram residuí modelu GARCH(1,1), gaussovský jádrový odhad hustoty rozdělení a teoretické normální rozdělení Jak je vidět rozdělení hodnot residuí sestaveného modelu GARCH(1,1) je špičatější s těžšími konci. V tomto případě by model nebyl vyhovující a v hodně případech se tento fakt přehlíží nebo je znám, ale přesto se model použije pro predikci. Nyní máme tedy dvě možnosti dalšího postupu. Buď model zamítneme a pokusíme se upravit parametry, tak aby model splnil tyto podmínky39 (nebo lze použít úplně jiný model) a nebo se pokusíme najít jiné rozdělení, kterým by se residua tohoto modelu mohla řídit. Provedeme tedy odhad modelu GARCH(1,1) znovu s předpokladem Studentova t-rozdělení (podle kapitoly 3.6).
ω α1 β1 ν
Odhad 1.240 · 10−06 4.944 · 10−02 9.467 · 10−01 4.719
Sm. odchylka 6.207 · 10−07 8.985 · 10−03 9.239 · 10−03 0.572 · 10−01
t-stat. 1.998 5.502 102.466 8.255
p-hodnota 0.0457 * 3.76 · 10−8 *** < 2 · 10−16 *** 2.22 · 10−16 ***
Tab. 11: Odhadnuté parametry modelu GARCH(1,1) s předpokladem Studentova t-rozdělení pro residua {ˆ }
R code: 1 # Odhad modelu GARCH s předpokladem studentova t−rozdělení 2 GARCHm = garchFit(˜ garch(1,1), data = ARMAe, include.mean = FALSE, cond.dist = ”std”) 3 # Residua modelu 4 GARCHe = residuals(GARCHm, standardize = TRUE)
39 V tomto případě je tato možnost vyloučena. Žádná ze základních úprav parametrů (přidávání, ubírání) ať už ARMA modelu nebo GARCH modelu nevede ke zlepšení či splnění podmínek normality standardizovaných residuí. Změnou počtu parametrů (nehledě na jejich významnost) dosáhneme minimálních změn.
— 45 —
6: Aplikační část
V tomto případě byl společně s parametry modelu GARCH(1,1) odhadován také parametr ν, což je tzv. shape parametr, který je odhadován jako počet stupňů volnosti t-rozdělení. Všimněme si také, že tento parametr skutečně vyšel v mezích (3, 6) uvedených v kapitole 3.6.
Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test Ljung-Box Test LM Arch Test
ˆ ˆ ˆ ˆ2 ˆ2 ˆ2 ˆ
Q(10) Q(15) Q(20) Q(10) Q(15) Q(20) T R2
Statistic 10.81694 12.29773 15.58943 4.80247 8.89128 9.47362 5.97691
p-Value 0.37196 0.65637 0.74175 0.90398 0.88313 0.97674 0.91724
Tab. 12: Testy standardizovaných residuí {ˆ } modelu GARCH(1,1) s předpokladem Studentova t-rozdělení Testy standardizovaných residuí prokázaly nekorelovanost a nepřítomnost podmíněné heteroskedasticity40 . Ovšem je opět nutné ověřit shodu s předpokládaným rozdělením daných residuí. Vhodným testem je například Kolmogorovův-Smirnovův test pro rovnost distribučních funkcí (popsaný v kapitole 2.2). Výsledkem K-S testu je p-hodnota rovna 0.9133. Tedy hypotéza o tom, že standardizovaná residua pocházejí ze studentova t-rozdělení, s počtem stupňů volnosti ν = 4.87399, není zamítnuta a model lze považovat za adekvátní. R code: 1 # Test t−rozdělení standardizovaných residuí 2 # Parametry 3 t1 = 0 4 t2 = 1 5 t3 = GARCHm@fit$coef[4] 6 # K−S test 7 ks.test(GARCHe, dstd(t1, t2, t3))
40 Také
bychom tento fakt mohli opět ověřit graficky pomocí ACF a PACF standardizovaných residuí a kvadrátů těchto residuí.
— 46 —
6: Aplikační část
Výsledek Kolmogorova-Smirnova testu ještě ilustrujeme následujícím obrázkem grafu odhadu hustoty standardizovaných residuí a t-rozdělení s počtem stupňů volnosti 4.87399.
Obr. 20: Histogram residuí {ˆ } modelu GARCH(1,1), gaussovský jádrový odhad hustoty rozdělení a teoretické Studentovo t-rozdělení
R code: 1 # Jádrový odhad distribuční funkce (student) 2 source(’.//functions//NKHst.r’) 3 G = NKHst(GARCHe, t1, t2, t3) 4 print(G)
6.5
Celkový konečný model
ARMA-GARCH je konečným modelem, kterým se rozumí klasický lineární model ARMA proces s podmíněně heteroskedstickými inovacemi. Při konstrukci tohoto modelu se nejprve odhadnul daný ARMA model, pro odstranění lineárních závislostí, a pro jeho residua se použil jeden z modelů podmíněné heteroskedasticity (GARCH). Model ARMA(1,0)-GARCH(1,1), resp. AR(1)-GARCH(1,1) má tvar: Xt =√µ + φXt−1 + εt , εt = ht t , ht = ω + α1 ε2t−1 + β1 ht−1 , kde Xt jsou logaritmy výnosů ceny zlata, |φ| < 1, t ∼ IID(0, 1), ω > 0, α1 ≥ 0, β1 ≥ 0. Pod označením IID41 náhodných veličin t se v tomto případě skrývá standardizované Studentovo t-rozdělení s odhadnutým parametrem ν = 4.719 stupňů volnosti [17]. 41 Independent
and Identically Distributed, tj. nezávislé a stejně rozdělené.
— 47 —
6: Aplikační část
Model s odhadnutými parametry: Xt =p 0.0007122 + 0.0608662Xt−1 + εˆt , ˆ t ˆt , εˆt = h ˆ t = 0.00000124 + 0.04944ˆ ˆ t−1 . h ε2t−1 + 0.9467h
Obr. 21: Podmíněný rozptyl odhadnutého modelu GARCH(1,1) a nepodmíněný rozptyl
R code: 1 # Vykreslení podmíněného rozptylu 2 3 df <− data.frame(dates = dates[(length(dates)−length(
[email protected])+1):length(dates)], 4 CondVar =
[email protected]) 5 { 6 G = ggplot(df, aes(x = dates)) 7 G = G + xlab(’Time’) + ylab(’Conditional variance’) 8 + geom line(aes(y=CondVar),size = 0.5, colour = ”grey”) 9 + geom point(aes(y=CondVar, colour = CondVar), size = 2) 10 + scale colour gradient(low = ”green”, high = ”red”) 11 + theme(legend.justification=c(1,1), legend.position=c(1,1)) 12 + geom hline(yintercept=0, color=”black”) 13 + geom hline(yintercept=var(ARMAe), linetype=’dashed’, color=”red”) 14 } 15 print(G)
— 48 —
6: Aplikační část
6.6
Monte Carlo předpovědi
V této části se pokusíme blíže analyzovat charakter předpovědi časové řady pomocí odhadnutého modelu z předchozích částí a k tomu využít metodu Monte Carlo. Postup aplikace této metody byl již popsán v kapitole 4.1. Tento postup tedy převedeme do jazyka R: R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
# MONTE CARLO # Parametry n.days = 20 # počet odhadovaných časových jednotek dopředu t3 # odhadnutý stupeň volnosti pro t−rozdělení k = 1000 # počet simulací # Inicializace k metodě Monte Carlo h.t = 0 h.t =
[email protected][length(
[email protected])] # hodnota podm. rozptylu v čase t e=0 e = ARMAe[length(ARMAe)] # hodnota residua ARMA modelu v čase t logProfit = 0 logProfit = LogProfit[length(LogProfit)] # hodnota residua logaritmu vynosů v čase t i=0 MCsim = seq(1:(n.days+1)) HT = seq(1:(n.days+1)) MClogProfit = seq(1:(n.days+1))
# Iterace metody Monte Carlo for (i in 1:k) { j=0 for (j in 1:n.days) { R = runif(k); # generovani nahodneho cisla MC = qt(R, df=t3);# hodnoty kvantilové funkce t−rozdělení pro simulované náhodné pravděpodobnosti
# simulovaný podm. rozptyl h.t[j+1] = GARCHm@fit$coef[1] + GARCHm@fit$coef[2]∗(e[j]ˆ2) + GARCHm@fit$coef[3]∗h.t[j] e[j+1] = MC∗sqrt(h.t[j+1]) # simulované residuum ARMA modelu logProfit[j+1] = ARMAm$coef[2] + ARMAm$coef[1]∗logProfit[j] + e[j+1] # simulované logaritmy výnosů }
# Uložení i−té iterace do pom. proměnné MCsim = data.frame(MCsim, e) HT = data.frame(HT, h.t) MClogProfit = data.frame(MClogProfit, logProfit) }
Nasimulovaná předpověď pro podmíněný rozptyl s odhadnutými parametry z předchozí kapitoly a simulovaná residua ARMA modelu lze vidět na následujících obrázcích, kde první hodnota je vždy hodnota odhadovaného modelu časové řady zlata z času t, což je čas poslední známé hodnoty časové řady cen zlata. Dále byl nastaven počet simulovaných hodnot na k = 1000 a počet předpovídaných dnů na h = 20.
— 49 —
6: Aplikační část
ˆ t na 20 dní dopředu pomocí metody MC Obr. 22: Simulované předpovědi podmíněného rozptylu h
Obr. 23: Simulované předpovědi na 20 dní pro residua εˆt modelu ARMA(1,0) pomocí metody MC R code: 1 # Graf MC předpovědí pro podmíněný rozptyl 2 A = melt(HT, id = ’HT’, variable name = ’series’) 3 G = ggplot(A, aes(HT, value)) + geom line(aes(colour = series)) 4 + theme(legend.position = ”none”) 5 + xlab(”Predicted days”) 6 + ylab(”Condtitional variance”) 7 + labs(title = paste(”Predicted conditional variance by Monte Carlo ”, k,” simulations for ”, n.days,” days”)) 8 print(G) 9 10 # Graf MC předpovědí pro residua ARMA modelu analogicky
— 50 —
6: Aplikační část
Z předcházejících obrázků je patrné, že simulace metody Monte Carlo vlastně vymezují oblast, ve které by se mohla nová budoucí data pohybovat. Toho lze využít několika způsoby. Zde bude uveden způsob určení intervalů pomocí empirických kvantilových funkcí. Obecně lze kvantilovou funkci (resp. inverzní distribuční funkci) napsat jako: Q(p) = F −1 (p)
nebo také
xp = F −1 (p),
kde pro danou pravděpodobnost p lze pomocí kvantilové funkce nalézt odpovídající hodnotu p-kvantilu xp náhodné veličiny X (tj. platí xp = F −1 (p) a F (xp ) = p). Inverzní funkci lze vytvořit pouze ke spojité rostoucí distribuční funkci, resp. pro spojitá rozdělení42 . Empirické kvantilové funkce lze získat z n dat vzestupným tříděním jejich hodnot xi a vytvořením tzv. schodovité funkce, kde osa y pak udává hodnoty xi a osa x pravděpodobnost 1/n jednotlivých xi , i = 1, ..., n. R code: 1 # Empirické intervalové odhady pro simulaci předpovědí MC 2 empInt = 0 3 i=0 4 j=0 5 Means = 0 # pro ukládaní průměrů simulací v daných dnech 6 MCsim1 = data.frame(days = seq(1:(n.days))) # pomocná proměnná pro graf 7 8 for (i in 1:n.days){ 9 P=0 10 for (j in 1:k) { 11 P[j] = MCsim[i+1, j+1] 12 MCsim1[i,j+1] = MCsim[i+1, j+1] 13 } 14 Means[i] = mean(P) 15 P = ecdf(P) # empirical cumulative distribution function 16 if (i==1) { 17 empInt = quantile(P, probs = c(seq(alpha/2, 0.5, 0.05), rev(1−seq(alpha/2, 0.5, 0.05))), 18 na.rm = FALSE, names = TRUE, type = 1) 19 } 20 else { 21 empInt = data.frame(empInt, empInt = quantile(P, probs = c(seq(alpha/2, 0.5, 0.05), 22 rev(1−seq(alpha/2, 0.5, 0.05))), na.rm = FALSE, names = TRUE, type = 1)) 23 } 24 }
42 Pro
případy diskrétních a smíšených rozdělení je třeba definici kvantilové funkce zobecnit: Q(p) = inf x|F (x) ≥ p.
— 51 —
6: Aplikační část
Následující graf ukazuje simulované hodnoty v jednotlivých dnech a 97.5%, 92.5%,. . ., 7.5%, 2.5% kvantily zkonstruované pomocí empirických kvantilových funkcí, vytvořených výše uvedeným postupem. Pro simulovaná data tedy byla vytvořena empirická kvantilová funkce v jednotlivých budoucích dnech a z těchto funkcí pak bylo možné získat jednotlivé procentní meze, které tvoří předpovědní intervaly.
Obr. 24: Simulované hodnoty residuí ARMA modelu a intervaly empirické kvantilové funkce.
R code: 1 # Příprava pro graf 2 P = names(quantile(P, probs = c(seq(alpha/2, 0.5, 0.05), rev(1−seq(alpha/2, 0.5, 0.05))), 3 na.rm = FALSE, names = TRUE, type = 1)) 4 5 empInt = melt(empInt); 6 empInt = data.frame(day = rep(seq(1:n.days), each = length(P)), quantiles = rep(P, n.days), empInt) 7 names(empInt) = c(”days”, ”percentage”, ”variable”, ”Quantiles”) 8 9 # Graf empirických intervalů vytvořených simulacemi 10 A = melt(MCsim1, id = ’days’, variable name=’series’) 11 G = ggplot(data = A, aes(days, value, group = series)) 12 + geom point(col = ”black”, alpha = 0.3) 13 + geom line(data = empInt, aes(days, Quantiles, group = percentage, colour = Quantiles)) 14 + scale colour gradient2(name = ”Values of\n quantiles ”, 15 midpoint = 0, low = ”red”, 16 mid = ”green”, high = ”red”, 17 guide = ”colourbar”) 18 + theme(legend.key.height = unit(2, ”cm”)) 19 print(G)
— 52 —
6: Aplikační část
6.6.1
Porovnání se skutečností
V této fázi porovnáme předpovědi sestavené pomocí metody Monte Carlo se skutečností. V projektu R načteme a přidáme další data do časové řady a opět provedeme několik prvních kroků, které jsou stejné jako na počátku zpracování původní časové řady. Jedná se tedy o zlogaritmování a diferencování nově přidaných dat, odstranění lineárních závislostí pomocí již odhadnutého ARMA(1,0) modelu. Tím získáme nová residua εt+1 , εt+2 , ..., εt+h (resp. jejich odhady), kde h = 20.
Obr. 25: Residua ARMA modelu a empirické intervaly spolehlivosti metody Monte Carlo.
Obr. 26: Residua ARMA modelu a 95% intervaly spolehlivosti pomocí příkazu predict.
— 53 —
6: Aplikační část
Předcházející graf na obrázku 25 ukazuje, jak dobře metoda Monte Carlo vymezila předpovědní oblast pro skutečné hodnoty. Pro porovnání pak graf na obrázku 26 ukazuje předpovědní interval sestavený pomocí příkazu predict, který užívá metodu předpovědi popsanou v kapitole 3.8. R code: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
# Příprava nových dat pro srovnaní predikce # Načtení dat s přidanými hodnotami (+20 dní) data2 <− read.csv(file = ”goldPred.csv”, head = FALSE, sep = ”;”)
# Rozdělení dat do skupin dates2 = rev(as.Date(data2$V1, format = ”%d.%m.%Y”)) serie2 = rev(data2$V2)
# Převod na stacionarní časovou řadu logSerie2 = log(serie2) # logaritmy dat lags2 = 1 # nastavení řádu diference LogProfit2 = diff(logSerie2, lag = lags2) # logaritmy výnosů
# Určení residuí ARMAe pomocí již odhadnutého modelu ARMA pro další doplněná data model = 0 i=0 a = 100 # počet posledních znamých dat
# Příprava pro Graf Predict for (i in 2:length(LogProfit2)) { model[i] = ARMAm$coef[2] + ARMAm$coef[1] ∗ LogProfit2[i] } ARMAePred = LogProfit2[(max(p,q)+1):length(LogProfit2)] − model[1:(length(model)−1)] K = ARMAePred[(length(ARMAePred)−(a+n.days)):length(ARMAePred)] # výběr pouze posledních a+n.days dat
# Graf předpovědí pomocí ”predict” + porovnaní se skutečností m = seq(1:(a+n.days)) predict(GARCHm, n.ahead = n.days, plot = TRUE, mse = ”cond”, nx = a, conf = (1−alpha)) lines(K[2:(a+n.days+1)]) points(K[2:(a+1)], col = ”black”, pch = 20, cex = 1.8) points(m[(a+1):(a+n.days)],K[(a+2):(a+n.days+1)], col = ”red”, pch = 20, cex = 2)
# Příprava pro Graf Monte Carlo K = ARMAePred[(length(ARMAePred)−(n.days)):length(ARMAePred)] # výběr nových (n.days) hodnot residuí K = data.frame(days = seq(1:n.days), real = K[2:length(K)]) # uložení do datového rámce s příslušnými dny M = ARMAe[(length(ARMAe)−a):length(ARMAe)] # výběr pouze (a) posledních hodnot residuí A3 = data.frame( days = rev(A$days[1:(length(A2)∗n.days)]), gr = rep(0, length(A2)∗n.days), residuals = rev(rep(Means, length(A2))), real = rep(rev(K$real), length(A2)), Intervals = rev(rep(rev(A2), each = n.days)), Intervals2 = (rep(rev(A2), each = n.days)), Upper = rev(melt(A1[length(A1):((length(A1)−1)/2+2)])$value), Lower = rev(melt(A1[2:((length(A1)−1)/2+1)])$value)) A4 = with(A3, A3[order(days), ]) # seřazení podle dnů i=0 for (i in 1:a) { A4 = data.frame( days = c(rep(−i+1, length(A2)), A4$days), gr = c(rep(1, length(A2)), A4$gr), residuals = c(rep(M[length(M)−i+1], length(A2)), A4$residuals), real = c(rep(M[length(M)−i+1], length(A2)),A4$real), Intervals = rep(levels(A4$Intervals),(length(P)+i)), Intervals2 = rev(rep(levels(A4$Intervals2), (length(P)+i))), Upper = c(rep(M[length(M)−i+1], length(A2)), A4$Upper), Lower = c(rep(M[length(M)−i+1], length(A2)), A4$Lower)) }
# Graf (A4) G = ggplot(A4, aes(x = days, y = residuals)) + geom ribbon(aes(ymin = Lower, ymax = Upper, fill = Intervals, group = Intervals2), alpha = 0.25) + geom line(aes(days, real), col = ”black”) + geom point(aes(days, real, colour = abs(real)), size = 3) + scale colour gradient(name = ”Residuals”, low = ”black”, high = ”red” ) print(G)
— 54 —
6: Aplikační část
Při porovnání předešlých obrázků 25 a 26, které ukazují posledních 100 + 20 hodnot residuí ARMA modelu zpracované časové řady, lze vidět, že 95% intervaly se u těchto dvou metod sestavení předpovědi liší. A protože příkaz predict neumí vykreslit graf s více intervaly, podívejme se na porovnání všech předpovědních intervalů daných metod v jiném grafu.
Obr. 27: Graf odhadů intervalů metody MC (červeně) a pomocí příkazu predict (modře)
R code: 1 # Porovnání predikčních metod pro danou simulaci 2 # Načtení dat (data simulace byla upravena v excelu a znovu načtena − v příloze na CD) 3 intervals <− read.csv(file = ”intervals.csv”, head = TRUE, sep = ”;”) 4 5 # Graf empirických intervalů 6 GA = ggplot(intervals, aes(x = days, y = lowerInterval2, group = Intervals)) 7 + scale shape identity() 8 + geom line(col = ”red”) 9 + geom point(aes(shape = 25), fill = ”red”, size = 4) 10 + geom line(aes(y = upperInterval2), col = ”red”) 11 + geom point(aes(y = upperInterval2, shape = 24), fill = ”red”, size = 4) 12 + geom line(aes(y = lowerInterval), col = ”blue”) 13 + geom point(aes(y = lowerInterval, shape = 25), fill = ”blue”, size = 3) 14 + geom line(aes(y = upperInterval), col = ”blue”) 15 + geom point(aes(y = upperInterval, shape = 24), fill = ”blue”, size = 3) 16 print(GA)
Z předchozího grafu vidíme jasné rozdíly mezi výsledky odhadu předpovědních intervalových mezí obou metod. Graf tedy ukazuje jak velké oblasti pokrývají jednotlivé pásy obou metod při dané α = 5%, 10%, ..., 50%. Z tohoto grafu však poměrně těžko posoudíme konkrétní rozdíly při dané α. Lepší interpretaci umožní následující grafy.
— 55 —
6: Aplikační část
Obr. 28: Graf odchylek mezi hornímy hranicemi intervalových pásů daných metod
Obr. 29: Graf odchylek mezi dolními hranicemi intervalových pásů daných metod R code: 1 # Graf odchylek mezi danými metodamy intervalových odhadů 2 G = ggplot(intervals, aes(x = days, y = varUpper, group = Intervals)) 3 + scale shape identity() 4 + geom line(col = ”red”) 5 + geom point(size = 3) 6 + geom text(aes(x = days+0.4, label = Intervals), size = 3) 7 print(G)
— 56 —
6: Aplikační část
Předchozí graf 28 ukazuje odchylky mezi jednotlivými horními hranicemi pásů intervalových odhadů. Naopak graf 29 ukazuje odchylky mezi dolními hranicemi pásů intervalových odhadů. Pokud se zaměříme na odchylky, které jsou v grafu nulové, resp. blízké nule, můžeme určit pro jaké pravděpodobnostní pásy dávají obě metody podobné výsledky. Pokud je odchylka kladná, pak metoda Metoda Carlo dává větší horní hranici a tedy širší intervalový pás. Pro zápornou odchylku platí, že metoda Monte Carlo dává nižší horní hranici a tedy užší intervalový pás pro predikci. Pro dolní intervalovou hranici to platí zdrcadlově, záporná odchylka říká, že metoda Monte Carlo dává širší pás a kladná odchylka užší pás pro dané pravděpodobnosti. Pokud se podíváme na pásy v čase t + 1, pak je z grafů odchylek patrné, že obě metody dávají podobné výsledky pro 85% pás v horní hranici intervalového odhadu a pro 90% pás v dolní hranici intevalového odhadu. Čím delší je předpověď, tím se hodnoty intervalových pásů obou metod potkávají v užších pravděpodobnostních pásech. V čase t + 20 jsou hodnoty intervalových odhadů podobné při 75% pravděpodobnosti horní i dolní hranice předpovědních pásů. V tomto případě si však uvědomme, že tyto výsledky platí pro danou simulaci metody Monte Carlo.
— 57 —
7: Závěr
7
Závěr
Cílem mé diplomové práce bylo stručně popsat základní charakteristiku finančních časových řad, modelů typu ARCH a popis aplikace těchto modelů, které jsou využity v aplikační části. Před zpracováním aplikační části jsem si vybral konkrétní časovou řadu ceny zlata, na kterou jsem aplikoval uvedené postupy a metody k získání výsledků ke všem cílům této práce. Dále jsem si vybral software Project R, který jsem s jeho vlastnostmi také popsal v teoretické části. V tomto softwaru jsem zpracoval všechny postupy a metody aplikační části a zároveň jsem ho také použil pro zpracování a interpretaci výsledků. V aplikační části jsem stručně popsal charakteristiku časové řady ceny zlata za období 9.2.2005 - 31.10.2012 a provedl předzpracování dat jako přípravu pro modelování této časové řady, ze kterého vyšly najevo některé typické vlastnosti finančních časových řad. Prvním krokem při modelování této řady bylo stanovení úrovňového modelu ARMA, druhým krokem stanovení modelu GARCH pro residua modelu ARMA. Při aplikaci těchto modelů jsem postupoval podle známých metod popsaných v textu práce, kde jsou popsány i postupy pro výběr a ověření modelů. Výsledkem zpracované finanční řady ceny zlata je, že tuto řadu lze modelovat pomocí modelu ARMA(1,0)-GARCH(1,1). Dalším cílem této práce bylo určení předpovědí pro zvolenou časovou řadu a k tomuto cíli využít metodu Monte Carlo. Základní princip této metody jsem také popsal v teoretické části a zároveň nastínil jakým způsobem je tuto metodu možné využít k předpovědím při modelování časové řady pomocí modelů typu ARCH. Předpovědní intervaly sestavené pomocí metody Monte Carlo pro časovou řadu ceny zlata jsou uvedeny v poslední části této práce a jsou porovnány se skutečnými hodnotami, pro které byly předpovědní intervaly odhadovány. Toto porovnání ukázalo, jak dobře dokáže tato metoda odhadnout předpovědi, ale také jak dobře byl celkový model pro časovou řadu odhadnut. Nové skutečné hodnoty nejsou nijak významně vychýleny z předpovědních intervalů sestavených pomocí metody Monte Carlo a můžeme tak tvrdit, že sestavený model a metodu odhadu přepovědi lze v tomto případě využít k dalšímu předpovězení vývoje časové řady ceny zlata. V této části bych chtěl také ohodnotit software Project R, který jsem si vybral a tímto způsobem také otestoval jeho využití k modelovaní finančních časových řad. Pomocí tohoto software je možné poměrně snadno připravit model pro určitou časovou řadu. Výsledky lze jednoduchým způsobem získat pomocí řady příkazů a funkcí, ale také mnoha manuálů a rad, které lze snadno získat z již vytvořené komunity okolo tohoto softwaru. Project R také slouží jako dobrý a pěkný nástroj pro interpretaci výsledků ve formě grafů, které lze naprosto jakkoliv naprogramovat a nadefinovat podle představ každého uživatele. Způsobů, jak řešit modelování časových řad je obecně mnoho a proto tato práce obsahuje mnoho potenciálu k dalšímu rozšíření či zkoumání. Například by bylo možné zkoumat a využít některé nelineární modely volatility a porovnat výsledky s lineárními modely volatility. Poměrně užitečným rozšířením by bylo zkoumání stability parametrů daných modelů pro časovou řadu ceny zlata nebo také citlivosti výsledků na změnu počtu simulací u metody Monte Carlo. Dalším rozšířením práce by také mohlo být zpracování dalších jiných typů předpovědí a porovnání výsledků s předpověďmi pomocí metody Monte Carlo.
— 58 —
LITERATURA A ZDROJE
Literatura a zdroje [1] J. Arlt, M. Arltová:Finanční časové řady. Grada Publishing a.s., Praha, 2003. ISBN 80-247-0330-0. [2] R. S. Tsay: Analysis of financial time series. Wiley, New York, 2002. [3] Z. Prášková: Základy náhodných procesů II. Karolinum, Praha, 2007. [4] J. Hendl: Přehled statistických metod. Portál, s.r.o. Praha, 2009. [5] J. D. Hamilton: Time Series Analysis. Princeton univerzity Press, Princeton, 1994. [6] J. Anděl: Matematická statistika. Státní nakladatelství technické literatury, Praha, 1985. [7] M. H. Kalos, P. A. Whitlock: Monte Carlo Methods. Second, Revised and Enlarged Edition, WILEY-VCH, Weinheim, 2008. [8] J. Mun: Modelling Risk: Applying Monte Carlo Simulation, Real Options Analysis, Forecasting, and Optimization Techniques. John Wiley & Sons, Inc., Hoboken, 2006. [9] R. T. Baillie, T. Bollerslev: Prediction in Dynamic Models with Time-Dependent Conditional Variances, Journal of GARCH, 1992, 52:91-113. Dostupné z: https://www.msu.edu/user/baillie/J_ Econometrics.1992.pdf [10] Shapiro, S.S., Wilk, M.B.: Biometrika. An Analysis of Variance Test for Normality (Complete Samples). 1965, roč. 52, č. 3/4, s. 591-611. Dostupné z: http://www.jstor.org/stable/2333709. [11] KURZYCZ [online]. [cit. 2012-12-08]. Dostupné z: http://www.kurzy.cz/komodity/index.asp?A=5&idk= 87&od=29.9.2003&curr=USD&unit=&lg=1&MAXROWS=50 [12] E. F. Fama: The Behaviour of Stock-Market Prices. Journal of Business, 1965 Volume 38. Dostupné z: http://www.ifa.com/Media/Images/PDFfiles/FamaThe_Behavior_of_Stock_Market_Prices.pdf [13] R. Cont: Empirical properties of asset returns: stylized facts and statistical issues. Quantitative Finance, 2001 Volume 1. Dostupné z: http://www.cmap.polytechnique.fr/~rama/papers/empirical.pdf [14] Analýza finančních časových řad typu ARCH. Plzeň, 2010. Dostupné z: https://stag-ws.zcu.cz/ws/ services/rest/kvalifikacniprace/downloadPraceContent?adipIdno=36997. Bakalářská práce. Západočeská univerzita. Vedoucí práce RNDr. Blanka Šedivá, Ph.D. [15] Analýza vývoje ceny ropy, zlata, vybraných blue chips a akciových indexů. Brno, 2011. Dostupné z: http:// is.muni.cz/th/206665/esf_m/Diplomka4.pdf. Diplomová práce. Masarykova univerzita. Vedoucí práce Ing. Gabriela Oškrdalová. [16] Struktura a vlastnosti modelu GARCH(1,1). Praha, 2009. Dostupné z: http://is.muni.cz/th/206665/ esf_m/Diplomka4.pdf. Bakalářská práce. Vysoká škola ekonomická v Praze. Vedoucí práce Ing. Jan Pígl. [17] Studentovo t-rozdělení a jeho aplikace. Olomouc, 2008. Dostupné z: http://mant.upol.cz/soubory/ OdevzdanePrace/B09/b09-13-aw.pdf. Bakalářská práce. Univerzita palackého. Vedoucí práce RNDr. Karel Hron, Ph.D. [18] The R Project: for Statistical Computing. [online]. [cit. 2013-01-05]. Dostupné z: http://www.r-project. org/.
— 59 —
A: Přílohy obsažené na CD
A
Přílohy obsažené na CD
P1 DPFIS PAKOSTA A10N0162P.pdf - soubor s kompletním textem diplomové práce v elektronické podobě. P2 DATA - složka obsahující všechna použitá data. P3 R - složka obsahující všechny zdrojové soubory, které obsahují celé zpracování časové řady v projektu R. P4 LATEX - složka obsahující zdrojový kód textu práce v LATEXu.
— 60 —
B: Přílohy tištěné
B
Přílohy tištěné
Obr. 30: Graf historické střední hodnoty logaritmů výnosů časové řady ceny zlata
Obr. 31: Graf historického rozptylu hodnoty logaritmů výnosů časové řady ceny zlata
— 61 —