Univerzita Karlova v Praze Matematicko-fyzikální fakulta
DIPLOMOVÁ PRÁCE
Eduard Hybler Moderní přístupy k analýze finančních časových řad
Katedra pravděpodobnosti a matematické statistiky Vedoucí diplomové práce: RNDr. Jitka Zichová, Dr. Studijní program: Matematika Studijní obor: Finanční a pojistná matematika
Rád bych na tomto místě poděkoval vedoucí mé diplomové práce RNDr. Jitce Zichové, Dr., za cenné připomínky a podněty k této práci, poskytnutí potřebné literatury, za umožnění práce s potřebným softwarem a za její čas strávený na konzultacích.
Prohlašuji, že jsem svou diplomovou práci napsal samostatně a výhradně s použitím citovaných pramenů. Souhlasím se zapůjčováním práce.
V Praze dne 10. 12. 2007
Eduard Hybler
2
Obsah Úvod
6
1 Korelace a regrese 1.1 Podmíněná nezávislost . . . . . . . . . . . . . . . 1.2 Korelační analýza . . . . . . . . . . . . . . . . . . 1.3 Odhady parametrů u mnohorozměrného rozdělení 1.4 Lineární regrese . . . . . . . . . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
7 7 8 11 12
2 Grafické modely 2.1 Základní pojmy z teorie grafů . . . . . 2.2 Graf podmíněných nezávislostí . . . . . 2.3 Markovské vlastnosti . . . . . . . . . . 2.4 Grafické modely a modelování . . . . . 2.5 Orientované acyklické grafy nezávislosti 2.6 Wermuthova podmínka a morální graf 2.7 Deviance . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
15 15 16 18 18 19 20 22
. . . . . . . . .
23 23 25 26 27 27 29 31 32 32
4 Grafické modely a vektorová autoregrese 4.1 Vektorová autoregrese . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Identifikace strukturálních autoregresí . . . . . . . . . . . . . . . . . . . .
34 34 35
5 Zpracování finančních dat 5.1 Identifikace kanonického ARMA modelu . . . . . . . . . . . . . . . . . . 5.2 Identifikace strukturální autoregrese pomocí grafických modelů . . . . . .
37 38 43
Závěr
50
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
3 Mnohorozměrné časové řady 3.1 Základní charakteristiky mnohorozměrného náhodného procesu . 3.2 Mnohorozměrný ARMA model . . . . . . . . . . . . . . . . . . . 3.2.1 Stacionarita a invertibilita posloupností ARMA . . . . . 3.3 Identifikace ARMA modelu . . . . . . . . . . . . . . . . . . . . 3.3.1 Výběrová autokovarianční a autokorelační matice . . . . 3.3.2 Výběrová parciální autokorelační matice . . . . . . . . . 3.3.3 Shrnutí volby řádu modelu . . . . . . . . . . . . . . . . . 3.3.4 Odhad parametrů modelu . . . . . . . . . . . . . . . . . 3.3.5 Verifikace modelu . . . . . . . . . . . . . . . . . . . . . .
3
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
Literatura
51
Přílohy Příloha č.1: Soubor vstupních dat . . . . . . . . . . . . . . . . . . . . . . . . . Příloha č.2: Obsah přiloženého CD . . . . . . . . . . . . . . . . . . . . . . . .
52 53 55
4
Název práce: Moderní přístupy k analýze finančních časových řad Autor: Eduard Hybler Katedra: Katedra pravděpodobnosti a matematické statistiky Vedoucí diplomové práce: RNDr. Jitka Zichová, Dr. e-mail vedoucího:
[email protected] Abstrakt: Práce se zabývá identifikací mnohorozměrného kanonického ARM A modelu a dále pak určením grafického modelu, který obsahuje i vztahy mezi proměnnými v současnosti, nikoli jenom závislost na zpožděných proměnných. V praktické části je pomocí výpočtů v softwaru Mathematica identifikovaný kanonický AR model, u kterého se pomocí metodologie grafických modelů určí strukturální autoregrese. Zpracována byla časová řada měnových kurzů, která je spolu s programem obsažena na přiloženém CD. Klíčová slova: Mnohorozměrný ARM A model, graf podmíněných nezávislostí, grafický model.
Title: Modern Approach To Financial Time Series Analysis Author: Eduard Hybler Department: Department of probability and mathematical statistics Supervisor: RNDr. Jitka Zichová, Dr. Supervisor’s e–mail address:
[email protected] Abstract: This thesis is devoted to the multivariate canonical ARM A model and then continues with determination of a graphical model. Graphical model includes also relations between contemporaneous variables, not only dependence on past variables. In the practical part of this work canonical AR model is identified using software Mathematica. In this model we specify structural autoregresion according to the graphical models methodology. A time series of exchange rates was processed. It is together with program included on the enclosed CD. Keywords: Multivariate ARM A model, conditional independence graph, graphical model.
5
Úvod Zpracování časových řad je jedním z úkolů řešených ve finanční analýze. Jedná se např. o studium vývoje cen akcií, měnových kurzů či hodnoty burzovních indexů. Finanční časové řady vykazují značnou nestabilitu, a proto u nich většinou dobře nefungují klasické přístupy.
V posledních letech byly vyvinuty nové modely vhodné právě pro práci s finančními časovými řadami, a to jak v jednorozměrném, tak v mnohorozměrném případě. Navíc v mnohorozměrném případě nás zajímají nejen vlivy minulých pozorování, ale i vztahy mezi jednotlivými složkami pozorovaného vektoru ve stejném čase.
Práce je rozdělena do pěti kapitol. V prvních čtyřech je vyložena teorie, která je pak v poslední kapitole aplikována na konkrétní data z finanční praxe. První kapitola se zabývá teorií korelace a regrese. Připomíná pojmy a metody potřebné v nasledujícím textu. Druhá kapitola obsahuje základy z teorie grafů, které jsou pak využité ve spojitosti s grafem podmíněných nezávislostí a jeho dalšími vlastnostmi. Třetí kapitola pojednává o mnohorozměrných časových řadách s důrazem na mnohorozměrný ARM A proces - zavedení základních pojmů, vlastností, způsoby identifikace, odhad parametrů a verifikace modelu. Čtvrtá kapitola uvádí postup kontrukce strukturálního AR modelu, jak je uvedeno v [9]. Poslední pátá kapitola aplikuje tento postup na časovou řadu - napozorovaných hodnot dvourozměrného vektoru kurzu eura a dolaru vůči české koruně - a popisuje získané výsledky.
6
Kapitola 1 Korelace a regrese V první kapitole připomeneme základní statistické pojmy jako jsou nezávislost, korelace a regrese - všechno se zřetelem na vícerozměrný případ - tj. uvažováním náhodných vektorů. Bez důkazů uvedeme některé vztahy potřebné pro naše další účely. Důkazy či podrobnější odvození vět a lemmat jsou uvedena např. v knize [2], kapitoly 2, 5 a 6 nebo v [4], zejména kapitoly 3 a 7. Důkazy vět a lemmat o parciálních kovariancích jsou v knize [11]. Poznamenejme ještě, že v celé práci budeme uvažovat sloupcové vektory. Transpozici vektoru X budeme značit X T .
1.1
Podmíněná nezávislost
Definice 1.1: Nezávislost náhodných veličin Náhodné veličiny X1 , . . . Xn , jsou nezávislé právě když FX (x1 , . . . xn ) = FX1 (x1 ) . . . FXn (xn ) ∀(x1 , . . . , xn )T ∈ Rn (tj. sdružená distribuční funkce se rovná součinu marginálních distribučních funkcí). ¤ Věta 1.2: Nechť náhodné veličiny X1 , . . . , Xn mají sdruženou hustotu fX a marginální hustoty fX1 , . . . , fXn . Pak X1 , . . . , Xn jsou nezávislé tehdy a jen tehdy, platí-li n Y fX (x1 , . . . xn ) = fXi (xi ). i=1
¤ Definice 1.3: Nezávislost náhodných vektorů Náhodné vektory X = (X1 , . . . Xn )T a Y = (Y1 , . . . Yk )T jsou nezávislé, právě když Xi a Yj jsou nezávislé náhodné veličiny ∀i = 1, . . . n a ∀j = 1, . . . k. ¤ Nezávislost budeme dále značit X ⊥ Y . Definici 1.3 můžeme ekvivalentně zapsat pomocí podmíněné hustoty: X ⊥ Y ⇔ fX|Y (x; y) = fX (x) ∀x. 7
Definice 1.4: Podmíněná nezávislost Náhodné vektory Y a Z jsou podmíněně nezávislé při pevné hodnotě vektoru X, právě tehdy když fY Z|X (y, z; x) = fY |X (y; x)fZ|X (z; x) pro všechny hodnoty y, z a pro všechna x taková, že fX (x) > 0.
¤
Podmíněnou nezávislost budeme dále značit Y ⊥ Z | X. Lemma 1.5: Nechť Y ⊥ Z | X. Potom jsou následující tvrzení ekvivalentní: (i) fY,Z|X (y, z; x) = fY |X (y; x)fZ|X (z; x) (ii) fY |X,Z (y; x, z) = fY |X (y; x) a analogicky fZ|X,Y (z; x, y) = fZ|X (z; x) f (x,y)fX,Z (x,z) (iii) fX,Y,Z (x, y, z) = X,Y fX (x) (iv) fX,Y,Z (x, y, z) = fX,Y (x, y)fZ|X (z; x) = fX,Z (x, z)fY |X (y; x) (v) fX,Y,Z (x, y, z) = fY |X (y; x)fZ|X (z; x)fX (x).
1.2
¤
Korelační analýza
Definice 1.6: Kovariance Kovarianci dvou náhodných vektorů X, Y s konečnými druhými momenty definujeme předpisem cov(X, Y ) = E(X − EX)(Y − EY )T . ¤ Pro praktické výpočty se používá vzorec cov(X, Y ) = EXY T − (EX)(EY )T . Definice 1.7: Varianční matice Mějme náhodný vektor X = (X1 , . . . , Xn )T . Matici V = varX = (cov(Xi , Xj )), která je typu n × n, nazýváme varianční maticí vektoru X. Platí pro ni varX = cov(X, X) = E(X − EX)(X − EX)T = EXX T − (EX)(EX)T . ¤ Definice 1.8: Korelační koeficient Nechť X, Y jsou náhodné veličiny s konečnými druhými momenty a kladnými rozptyly. Pak definujeme korelační koeficient cov(X, Y ) . cor(X, Y ) = √ varXvarY ¤ 8
Definice 1.9: Korelační matice Mějme náhodné vektory X = (X1 , . . . , Xn )T a Y = (Y1 , . . . , Ym )T s konečnými druhými momenty. Nechť všechny složky těchto dvou vektorů mají kladné rozptyly. Matici typu n × m, jejíž (i, j)-tý prvek je roven korelačnímu koeficientu veličin Xi a Yj , nazveme korelační maticí vektorů X a Y a označíme ji cor(X, Y ). Matice cor(X, X) se nazývá korelační matice vektoru X a značí se corX. ¤ Definice 1.10: Parciální kovariance Nechť X, Y , Z jsou náhodné vektory s konečnými druhými momenty a varianční matice varX je regulární matice. Pak parciální kovarianci Y a Z při pevném X počítáme jako ˆ cov(Y, Z | X) = cov(Y − Yˆ , Z − Z), kde Yˆ , Zˆ jsou nejlepší lineární odhady Y, Z získané z X. Je tedy Yˆ = α1 + β1T X, kde α1 = EY − β1T EX, β1 = (varX)−1 cov(X, Y ) Zˆ = α2 + β2T X, kde α2 = EZ − β2T EX, β2 = (varX)−1 cov(X, Z).
¤
Vektory Y − Yˆ , Z − Zˆ (tj. ty části Y a Z, které vektor X nevysvětlí) nazýváme vektory reziduí. Definice 1.11: Parciální rozptyl Nechť X a Y jsou náhodné vektory. Pak parciální rozptyl Y při pevném X definujeme jako var(Y | X) = cov(Y, Y | X). ¤ Definice 1.12: Parciální korelační koeficient Nechť Y , Z jsou náhodné veličiny s konečnými druhými momenty a X nሠ> 0 definujeme hodný vektor. Za předpokladu var(Y − Yˆ ) > 0 a var(Z − Z) parciální korelační koeficient veličin Y a Z při pevném X jako cov(Y, Z | X)
cor(Y, Z | X) = p
var(Y | X)var(Z | X)
. ¤
Definice 1.13: Parciální korelační matice Nechť X, Y = (Y1 , . . . , Ym )T a Z = (Z1 , . . . , Zk )T jsou náhodné vektory s konečnými druhými momenty. Matici typu m × k, jejíž (i, j)-tý prvek je roven koeficientu parciální korelace veličin Yi a Zj při pevném X, nazveme parciální korelační maticí vektorů Y a Z při pevném X a označíme ji cor(Y, Z | X). Platí, že ˆ cor(Y, Z | X) = cor(Y − Yˆ , Z − Z). ¤ 9
Protože na hodnotách α1 a α2 , jak byly zavedeny v definici 1.10, parciální korelační matice ˆ nezávisí, platí rovněž cor(Y − Yˆ , Z − Z) cor(Y, Z | X) = cor(Y − β1T X, Z − β2T X). Lemma 1.14: Platí cov(Y, Z | X) = cov(Y, Z) − cov(Y, X)var(X)−1 cov(X, Z). ¤ Lemma 1.15: Platí var(Y | X) = var(Y ) − var(Yˆ ) = cov(Y, Y ) − cov(Y, X)var(X)−1 cov(X, Y ) = var(Y − Yˆ ) a číslo var(Y − Yˆ ) se nazývá reziduální rozptyl.
¤
Věta 1.16: Nechť X, Y, Z jsou náhodné vektory s normálním rozdělením. Pak platí X⊥Y
⇔
cov(X, Y ) = 0
X⊥Y |Z
⇔
cov(X, Y | Z) = 0. ¤
Věta 1.17: Nechť X = (X1 , . . . , Xn )T je náhodný vektor s konečnými druhými momenty, jehož√všechny složky mají kladné rozptyly. Označme V = varX, P = corX, σi = varXi pro i = 1, . . . , n a σ1 0 . . . 0 D = Diag{σ1 , . . . , σn } = . . . . . . . . . . . . . 0 Pak P = D−1 V D−1 .
0
. . . σn ¤
Věta 1.18: Lemma o inverzní varianci Mějme náhodné vektory X a Y s konečnými µ ¶ druhými momenty. Inverzní X varianční matici náhodného vektoru lze za předpokladu regularity Y matic varX a var(Y | X) psát ve tvaru ¶ µ var(X)−1 + B T var(Y | X)−1 B −B T var(Y | X)−1 T T T −1 , (var(X , Y ) ) = −var(Y | X)−1 B var(Y | X)−1 kde B = cov(X, Y )(varX)−1 .
¤ 10
Důsledek 1.19: Mějme k-rozměrný náhodný vektor X, K = {1, 2, . . . , k}. Pak každý diagonální prvek inverzní varianční matice V −1 = (wij ) je převrácenou hodnotou parciálního rozptylu. Tedy wii =
1 var(Xi | XK\{i} )
∀i ∈ K.
A dále každý mimodiagonální prvek inverzní varianční matice škálované tak, aby na diagonále byly jedničky, je záporně vzatá parciální korelace odpovídajících složek vektoru X při pevně daných ostatních složkách. Tj. C = DV −1 D = (cij ), kde
½
¾ 1 1 D = Diag √ ,..., √ . w11 wnn cij = √
wij = −cor(Xi , Xj | XK\{i,j} ) ∀i ∈ K, ∀j ∈ K, i 6= j. wii wjj
Odhady parciálních korelací získáme škálováním výběrové varianční matice. ¤ Důsledek 1.20: Mějme k-rozměrný náhodný vektor X s normálním rozdělením, K = {1, 2, . . . , k}. Pak Xi ⊥ Xj | XK\{i,j} ⇔ cij = 0 ∀i ∈ K, ∀j ∈ K, i 6= j, kde cij je jako v důsledku 1.19.
¤
Jestliže pro parciální korelační koeficient platí ρ = 0 a pro počet pozorování N > p + 2, kde p je rozměr podmiňujíciho vektoru, pak p r T =√ N − p − 2 ∼ tn−p−2 . (1.1) 1 − r2 Vztah (1.1) nám umožňuje testovat nulovost parciálních korelačních koeficientů.
1.3
Odhady parametrů u mnohorozměrného rozdělení
Mějme nějaké m-rozměrné rozdělení se střední hodnotou µ = (µ1 , . . . , µm )T a varianční maticí V typu m × m. Nechť xN,1 x1,1 x1 = ... , . . . , xN = ... xN,m
x1,m 11
je náhodný výběr z tohoto rozdělení. Omezíme se na případ N > m. Zaveďme výběrový průměr x¯ a výběrovou varianční maticí Vˆ = (ˆ vij ) pomocí vzorců: x¯ =
N 1 X xi N i=1
N
1 X 1 Vˆ = (xi − x¯)(xi − x¯)T = N − 1 i=1 N −1
µX N
xi xTi
¶ − N x¯x¯ . T
i=1
V případě nulové střední hodnoty můžeme uvedený vztah přepsat na Vˆ = kde
1 XX T , N −1
(1.2)
x1,1 . . . xN,1 .. . X = ... . x1,m . . . xN,m
Jsou-li všechny diagonální prvky matice Vˆ kladné, definujeme výběrovou korelační matici µ ¶m v ˆ ij ˆ XX = p R . vˆii vˆjj i,j=1 Dále zavedeme výběrovou kovarianční matici vektorů X a Y : N
VˆXY =
1.4
1 X (xi − x¯)(yi − y¯)T . N − 1 i=1
Lineární regrese
Mějme náhodné veličiny Y1 , . . . , Yn a matici daných čísel X = (xi,j ) typu n × k, kde k < n. Požaduje se, aby matice X měla lineárně nezávislé sloupce. Předpokládejme, že pro náhodný vektor Y = (Y1 , . . . , Yn )T platí Y = Xβ + e
(1.3)
kde β = (β1 , . . . , βk )T je vektor neznámých parametrů a e = (e1 , . . . , en )T je vektor náhodných chyb s normálním rozdělením s parametry (0, σ 2 I), kde 0 je nulový vektor, σ 2 > 0 je neznámý parametr a I je jednotková matice typu n × n. Tento model budeme nazývat regresní. Protože Y závisí na β lineárně, mluvíme o lineárním regresním modelu. Parametry (β1 , . . . , βk )T se odhadují metodou nejmenších čtverců, tj. minimalizací výrazu (Y − Xβ)T (Y − Xβ) jakožto funkce β. Tento odhad označme b = (b1 , . . . , bk )T . Velmi často se stává, že první sloupec matice X je tvořen samými jedničkami. Položme r = k − 1. Pak můžeme psát X = (1, X1 ), kde X1 je matice typu n × r a 1 je sloupcový vektor jedniček. Vektory β a b se pak zapisují ve tvaru β = (β0 , β1 , . . . , βr )T a b = (b0 , b1 , . . . , br )T . 12
Věta 1.21: Odhad vektoru parametrů β metodou nejmenších čtvrců je b = (X T X)−1 X T Y. Pro tento odhad platí E(b) = β, var(b) = σ 2 (X T X)−1 .
¤
Vektor Yˆ = Xb = X(X T X)−1 X T Y může být považován za nejlepší aproximaci vektoru Y , jaká se dá vytvořit lineární kombinací sloupců matice X. Ještě uveďme, že se veličině Se = (Y − Yˆ )T (Y − Yˆ ) říká reziduální součet čtverců. Věta 1.22: Pro reziduální součet čtverců Se platí Se = Y T (I − X(X T X)−1 X T )Y = Y T Y − bT X T Y. ¤ Věta 1.23: Náhodná veličina s2 = Se /(n − k) je nestranným odhadem parametru σ 2 . ¤ Věta 1.24: Rozdělení a vlastnosti odhadnutých parametrů Platí: (i) b ∼ N (β, σ 2 (X T X)−1 ) (ii) Se /σ 2 ∼ χ2n−k (iii) vektor b a veličina s2 jsou nezávislé. Označme ST =
n X
(Yi − Y¯ )2 =
i=1
¤ n X
Yi2 − nY¯ 2 .
i=1
Kvalitu odhadu Yˆ náhodné veličiny Y posuzujeme například koeficientem determinace ¯ 2 , které jsou dány vzorci R2 a korigovaným koeficientem determinace R R2 = 1 −
Se , ST
¯2 = 1 − R
n − 1 Se . n − r − 1 ST
Čím jsou tyto koeficienty blíže jedné, tím těsnější je regresní závislost.
13
T-test pro jeden regresní parametr Věta 1.25: Označme si uij prvky matice (X T X)−1 . Potom bi − βi Ti = √ 2 ∼ tn−k , s uii
i = 1, . . . , k. ¤
Předchozí věta umožňuje testovat hypotézy tvaru H0 : βi = βi0 proti alternativám H1 : βi 6= βi0 , kde i ∈ {1, 2, . . . , k} a βi0 je dané reálné číslo. Jako testovou statistiku použijeme bi − β 0 Ti = √ 2 i . s uii V případě, že | Ti |≥ tn−k (α), zamítáme hypotézu H0 na hladině α. Nejčastěji volíme βi0 = 0 a ověřujeme, zda Y vůbec závisí na Xi . tn−k (α) je kritická hodnota příslušného Studentova rozdělení.
14
Kapitola 2 Grafické modely Druhá kapitola je věnovaná připomenutí základních pojmů z teorie grafů, o kterých najdeme více např. v knize [7]. Tyto poznatky pak tvoří základ pro další teorii týkající se např. grafu podmíněných nezávislostí, markovských vlastností či grafických modelů. Věty jsou opět uváděny bez důkazu, s odvoláním se na [11], kapitoly 3 a 5. Poslední část se zabývá deviancí, jak je vyložena v článku [9].
2.1
Základní pojmy z teorie grafů
Definice 2.1: Graf Graf G je uspořádaná dvojice (K, E), kde K je nějaká neprázdná množina a E je podmnožina kartézského součinu K × K. Prvky množiny K se jmenují vrcholy grafu G a prvky množiny E hrany grafu G. Říkáme, že vrcholy i, j jsou spojeny hranou právě tehdy, když (i, j) ∈ E. ¤ Množina K je obvykle podmnožina přirozených čísel, K = {1, 2, . . . , k}. Definice 2.2: Úplný graf Graf označujeme jako úplný, jsou-li všechny dvojice vrcholů spojené hranou, tj. ∀ i, j ∈ K; (i, j) ∈ E. ¤ Definice 2.3: Orientovaný graf Orientovaný graf G je dvojice (K, E), kde E je podmnožina kartézského součinu K × K. Prvky E nazýváme šipky (nebo orientované hrany). Šipka e má tvar (i, j). Říkáme, že tato šipka vychází z i a končí v j. O vrcholu i pak mluvíme jako o rodiči vrcholu j. Je-li alespoň jedna hrana grafu orientovaná, mluvíme o orientovaném grafu. ¤ Definice 2.4: Podgraf, faktor, klika Řekneme, že graf G1 = (K1 , E1 ) je podgrafem grafu G = (K, E), jestliže K1 ⊆ K a E1 ⊆ E. Pokud platí K1 = K říkáme, že graf G1 je faktorem grafu G. Klika je maximální úplný podgraf. ¤ 15
Definice 2.5: Hranice vrcholu Hranice vrcholu i je podmnožina množiny vrcholů K obsahující všechny vrcholy, které jsou s vrcholem i spojeny hranou. Tuto hranici vrcholu budeme značit bd(i). ¤ Definice 2.6: Cesta, dostupnost, souvislost Cesta je posloupnost vrcholů i1 , i2 , . . . , im takových, že (il , il+1 ) ∈ E
∀ l = 1, 2, . . . , m − 1.
Nejkratší cesta je cesta s nejmenším počtem hran. Vrcholy označujeme jako dostupné právě tehdy, když existuje cesta z i do j a z j do i. Graf označujeme jako souvislý právě tehdy, když všechny dvojice vrcholů jsou navzájem dostupné. ¤ Definice 2.7: Cyklus Orientovaný graf G = (K, E) s množinou vrcholů K = k1 , . . . , kn a množinou orientovaných hran E = (e1 , . . . , en ) nazýváme cyklus, jestliže platí ei = (ki , ki+1 ) pro i = 1, . . . , n − 1 a en = (kn , k1 ). ¤ Definice 2.8: Acyklický graf Orientovaný graf neobsahující cykly nazveme acyklický.
¤
Definice 2.9: Separace Říkáme, že A ⊆ K separuje dva vrcholy i, j právě tehdy, když každá cesta z i do j (z j do i) obsahuje alespoň jeden vrchol ze separující množiny A. Množina A separuje dvě podmnožiny vrcholů B a C právě tehdy, když separuje každou dvojici vrcholů i ∈ B, j ∈ C. ¤
2.2
Graf podmíněných nezávislostí
Definice 2.10: Graf podmíněných nezávislostí Grafem podmíněných nezávislostí náhodného vektoru X = (X1 , . . . , Xk )T nazýváme takový graf G = (K, E), v němž K = {1, 2, . . . , k} a (i, j) ∈ / E ⇔ Xi ⊥ Xj | XK\{i,j}
∀i, j ∈ K. ¤
Graf podmíněných nezávislostí se označuje CIG (z ang. conditional independence graph).
16
Lemma 2.11: Předpokládejme, že množina vrcholů K = {1, 2, . . . , k} může být rozdělena do dvou podmnožin B, C takových, že neexistuje žádná cesta spojující libovolný vrchol z B s libovolným vrcholem z C. Potom Xi ⊥ Xj
∀i ∈ B, ∀j ∈ C. ¤
Lemma 2.12: Předpokládejme, že množina vrcholů K = {1, 2, . . . , k} může být rozdělena do dvou podmnožin B, C takových, že neexistuje žádná cesta spojující libovolný vrchol z B s libovolným vrcholem z C. Potom Xi ⊥ Xj | XA
∀i ∈ B, ∀j ∈ C
a A je libovolná podmnožina K neobsahující i ani j.
¤
Lemma 2.13: Je-li A podmnožina vrcholů K, která separuje dva vrcholy i a j, pak Xi ⊥ Xj | XA . ¤ Věta 2.14: Věta o separaci Nechť XA , XB a XC jsou vektory obsahující disjunktní podmnožiny složek náhodného vektoru X a nechť v grafu podmíněných nezávislostí vektoru X každý vrchol v B je separován od všech vrcholů z C podmnožinou A. Potom XB ⊥ XC | XA . ¤ Definice 2.15: Minimální podmíněná nezávislost, minimální separační množina Říkáme, že podmíněná nezávislost mezi dvěma proměnnými je minimální, pokud není možné eliminovat jakoukoliv proměnnou z podmiňujícího vektoru XA . Množinu A jakožto separační množinu nazýváme minimální. ¤
17
2.3
Markovské vlastnosti
V grafických modelech zavádíme 3 typy markovských vlastností: Definice 2.16: Lokální markovská vlastnost Říkáme, že náhodný vektor X s grafem podmíněných nezávislostí G = (K, E) má lokální markovskou vlastnost pokud platí Xi ⊥ XB | XA , kde i ∈ K je libovolný vrchol, A = bd(i) je hranice vrcholu i a B = K \ ({i} ∪ A). Zjednodušeně můžeme psát i ⊥ zbytek | hranice. ¤ Definice 2.17: Párová markovská vlastnost Pro všechny vrcholy i, j, které nejsou spojeny hranou, platí Xi ⊥ Xj | XA , kde A = K \ {i, j}.
¤
Definice 2.18: Globální markovská vlastnost Nechť A, B, C jsou po dvou disjunktní podmnožiny K. Pokud A separuje B a C, pak platí XB ⊥ XC | XA . ¤ Věta 2.19: Ekvivalence markovských vlastností Uvedené tři markovské vlastnosti jsou ekvivalentní.
2.4
¤
Grafické modely a modelování
Definice 2.20: Obecný grafický model Uvažujme k-rozměrný náhodný vektor X = (X1 , X2 , . . . , Xk )T a graf podmíněných nezávislostí G = (K, E). Grafický model pro náhodný vektor X je rodina pravděpodobnostních rozdělení vektoru X, která splňuje podmíněné nezávislosti dané grafem G. ¤
18
Definice 2.21: Grafický Gausovský model V případě, že vektor X z definice 2.20 má mnohorozměrné normální rozdělení, mluvíme o grafickém Gausovském modelu. V tomto případě platí, že podmínky podmíněné nezávislosti jsou ekvivalentní s určením nulových prvků v inverzní varianční matici. ¤ Jednoduchá procedura pro nalezení grafu podmíněných nezávislostí má nasledující kroky: 1. Odhadne se varianční matice V = var(X) pomocí výběrové varianční matice Vˆ . 2. Spočte se inverzní matice Vˆ −1 . 3. Škálováním inverzní matice, jak je uvedeno v důsledku 1.19, získáme diagonální matici, ze které určíme výběrové parciální korelace cor(Xi , Xj | zbytek). 4. Testujeme nulovost prvků výběrové parciální korelační matice. Nakreslí se výsledný graf podmíněných nezávislostí. Definice 2.22: Saturovaný model Saturovaný model je grafický model určený úplným grafem.
2.5
¤
Orientované acyklické grafy nezávislosti
Dostáváme se k situaci, kdy X ovlivňuje Y , ale Y neovlivňuje X. To můžeme vyjádřit orientovaným grafem na obrázku 2.23. ²¯
X ±°
²¯ - Y ±°
Obrázek 2.23: Jednoduchý orientovaný graf Pokud se budeme těmito modely dále zabývat, můžeme se dostat i k situacím znázorněnými na obrázku 2.24. ²¯
²¯
W ±°
Y ±°
* ©© © ²¯ © X ±° YH H ? HH²¯
6
²¯ - X ±°
²¯ Z ¾ ±°
Z ±°
Obrázek 2.24: Cyklické grafy
19
? ²¯
Y ±°
Tyto cykly nebudeme v dalším textu uvažovat, neboť neumožňují dobře modelovat sdružené rozdělení příslušných proměnných pomocí podmíněných rozdělení. Potřeba vyloučit orientované cykly nás vede k následujícímu: Na prvky množiny K zavedeme relaci úplného uspořádání  splňující (pro ∀i ∈ K, ∀j ∈ K): (i) buď i ≺ j nebo i  j (ii) relace  je ireflexivní, tj. non(i  i), (iii) relace  je tranzitivní, tj. jestliže (i ≺ j) ∧ (j ≺ k) ⇒ (i ≺ k). Pokud uvedené aplikujeme na orientovaný graf, dostáváme, že každá strana grafu může mít jen jednu možnou orientaci. Lemma 2.25: V orientovaném grafu jsou následující podmínky ekvivalentní: (i) graf neobsahuje orientované cykly (ii) existuje úplné uspořádání vrcholů grafu.
¤
Definice 2.26: Orientovaný graf podmíněných nezávislostí Orientovaným grafem podmíněných nezávislostí náhodného vektoru X nazýváme orientovaný graf GÂ = (K, E Â ), kde K = {1, 2, . . . , k}, K(j) = {1, 2, . . . , j} a hrana (i, j), i ≺ j, není obsažena v množině hran právě tehdy když Xj ⊥ Xi | XK(j)\{i,j} . ¤ Orientovaný graf podmíněných nezávislosti se označuje DAG (z ang. directed acyclic graph). Poznámka 2.27: DAG implikuje jediný CIG pro složky vektor X. Naopak k danému CIG může existovat více DAG, nebo nemusí existovat žádný.
2.6
Wermuthova podmínka a morální graf
Definice 2.28: Asociovaný neorientovaný graf Je-li GÂ = (K, E Â ) jako v definici 2.26, potom asociovaný neorientovaný graf definujeme jako Gu = (K, E u ) se stejnou množinou vrcholů, v němž orientované hrany jsou nahrazeny neorientovanými. ¤ Definice 2.29: Wermuthova podmínka Říkáme, že orientovaný graf splňuje Wermuthovu podmínku jestliže neobsahuje podgraf znázorněný na obrázku 2.30.
20
²¯ HH ²¯ ±° j H *±° © ²¯ ©© ±°
Obrázek 2.30: Podgraf porušující Wermuthovu podmínku ¤ Definice 2.31: Morální graf Morální graf k orientovanému grafu G = (K, E  ) je takový neorientovaný graf Gm = (K, E m ), který má stejnou množinu vrcholů K a obsahuje množinu hran E  a navíc obsahuje všechny hrany, tzv. morální, nezbytné k tomu, aby byla v G splněna Wermuthova podmínka. ¤ Název morální vystihuje skutečnost, že tento graf ”sezdává rodiče vrcholu”. Příklad 2.32: Uvažujme orientovaný graf G1 = (K1 , E1 ), kde K1 = {1, 2, 3, 4} a E1 = {(1, 4)(2, 4)(2, 3)} (obrázek 2.33a). ²¯
²¯ - 4 ±° ¡ µ ¡
1 ±° ¡
¡ ²¯ ¡ 2 ±°
(a)
²¯
- 3 ±°
²¯
²¯
1 ±°
4 ±°
¡ ¡
¡
¡ ²¯ ¡ 2 ±°
(b)
²¯
3 ±°
²¯
1 ±°
²¯
¡ ²¯ ¡ 2 ±°
4 ±° ¡ ¡ ¡
(c)
²¯
3 ±°
Obrázek 2.33: (a) graf G1 , (b) asociovaný graf Gu1 k (a), (c) morální graf Gm 1 odvozený od (b) Sestrojme si jeho asociovaný graf Gu1 (obrázek 2.33b). V grafu G1 nám vrcholy 1, 2, 4 porušují Wermuthovu podmínku. Při konstrukci morálního grafu tedy musíme přidat hranu (1, 2). Morální graf je na obrázku 2.33(c). ¤
Markovské vlastnosti Věta 2.34: Markovský teorém pro orientované grafy podmíněných nezávislostí Orientovaný graf GÂ má markovské vlastnosti svého morálního grafu Gm . ¤ 21
Důsledek 2.35: Pro Gm = Gu má orientovaný graf GÂ právě jenom markovské vlastnosti svého morálního grafu Gm . ¤ Pokud má totiž morální graf více hran než asociovaný, obsahuje graf G i některé nezávislosti, které nelze z markovských vlastností morálního grafu odvodit.
Vraťme se k příkladu 2.32: Z markovských vlastností morálního grafu dostáváme X1 ⊥ X3 | X2 , X4 X3 ⊥ X4 | X1 , X2 . V grafu G1 nejsou vrcholy 1, 2 spojeny hranou a z definice 2.10 tedy pro Gu1 vyplývá, že X1 ⊥ X2 | X3 , X4 a odtud dle definice 2.26 pak X1 ⊥ X2 pro graf G1 navíc k podmíněným nezávislostem odvozeným z morálního grafu.
2.7
Deviance
Deviance představuje základní způsob, jak měřit shodu grafického modelu s daty, případně způsob, jak mezi sebou porovnat dva grafické modely, pokud je jeden faktorem druhého. Definice 2.36: Deviance Mějme náhodný vektor X = (X1 , . . . , Xk )T s rozdělením N (0, V ) a pozorované realizace x1 , . . . , xN . Devianci modelu DAG definujeme vztahem X Dev = N log s2r , r
kde N je rozsah výběru, r je počet regresních rovnic (pro popis závislostí každé ze složek mnohorozměrné časové řady v čase t na ostatních složkách v čase t a na zpožděných hodnotách. Tyto rovnice jsou dány modelem DAG, jak bude popsáno v kapitolách 4 a 5) a s2r je odhad reziduálního rozptylu jednotlivých rovnic. ¤ Saturovaný model má pak devianci Dev = N log det Vˆ , kde Vˆ je výběrová varianční matice vektoru X. Model s méně hranami má větší devianci. Pokud dojde po vynechání hrany k výraznímu nárůstu deviance, neměli bychom ji z modelu vynechávat. 22
Kapitola 3 Mnohorozměrné časové řady V kapitole o mnohorozměrných časových řadách jsou nejprve uvedeny základní definice jako střední hodnota, autokorelační a autokovarianční maticové funkce. Následuje podrobnější popis mnohorozměrných ARM A posloupností. Věty jsou opět uváděné bez důkazů. Ty je možné najít v [10], kapitoly 2 a 3, nebo [3], 11. a 12. kapitola, a také v [8].
3.1
Základní charakteristiky mnohorozměrného náhodného procesu
Definice 3.1: Náhodný proces Nechť Xt = (Xt,1 , . . . , Xt,m )T je m-rozměrný náhodný vektor. Systém {Xt , t ∈ T } se nazývá m-rozměrný náhodný proces. Je-li speciálně T podmnožina celých čísel, nazývá se {Xt , t ∈ T } proces s diskrétním časem nebo také časová řada. ¤ Pokud je z kontextu zřejmé, o jakou množinu T se jedná, budeme psát pouze {Xt }. Definice 3.2: Střední hodnota Nechť {Xt , t ∈ T } je m-rozměrný náhodný proces takový, že pro každé t ∈ T existuje střední hodnota EXt . Potom vektorová funkce µt = EXt definovaná na T nazývá střední hodnota procesu {Xt }. Proces, jehož střední hodnota je identicky rovna nule, se nazývá centrovaný. ¤ Platí-li E|Xt,k |2 < ∞ pro všechna t ∈ T a k = 1, 2, . . . , m, říká se, že proces {Xt , t ∈ T } má konečné druhé momenty. Je-li tato podmínka splněna, existuje i střední hodnota procesu. Definice 3.3: Autokovarianční maticová funkce procesu Jestliže {Xt , t ∈ T } je m-rozměrný náhodný proces s konečnými druhými momenty, potom maticová funkce dvou proměnných definována na T × T předpisem R(s, t) = E(Xs − µs )(Xt − µt )T 23
se nazývá autokovarianční funkce procesu {Xt }. Funkce R(t, t) se nazývá rozptylová maticová funkce náhodného procesu {Xt } v čase t. ¤ V některých případech se stává, že funkce R(s, t) závisí na svých argumentech pouze prostřednictvím jejich rozdílu. V tomto případě zavádíme funkci jedné proměnné vztahem R(s − t) = R(s, t). Definice 3.4: Slabá stacionarita Náhodný proces {Xt , t ∈ T } s konečnými druhými momenty se nazývá slabě stacionární, má-li konstantí střední hodnotu µt = µ pro všechna t ∈ T a je-li jeho autokovarianční funkce R(s, t) funkcí pouze rozdílu s − t. Pak zkráceně píšeme R(k) = R(t, t + k) a R(0) nazýváme rozptyl procesu. Je-li splněna pouze podmínka na autokovarianční funkci, mluvíme o kovarianční stacionaritě. ¤ Autokovarianční funkci slabě stacionárních procesů lze definovat jako funkci jedné proměnné vztahem R(t) := R(t, 0), t ∈ T ; autokorelační funkce je potom definována předpisem r(t) = R(t)/R(0). Pro t ∈ Z zřejmě platí R(−k) = R(k)T , proto často matici R(k) počítáme pouze pro t ∈ N0 . Definice 3.5: Autokorelační maticová funkce procesu Nechť {Xt , t ∈ T } je slabě stacionární náhodná posloupnost s autokovarianční maticovou funkcí R(k) = [Rij (k)]i,j=1,2,...,m . Označme D diagonální matici n o 1 1 1 ∆ = Diag p ,p ,..., p . R11 (0) R22 (0) Rmm (0) Pak pro každé k ∈ Z definujeme autokorelační maticovou funkci r(k) jako r(k) = ∆R(k)∆. ¤ Definice 3.6: Parciální autokorelační maticová funkce Nechť {Xt } je slabě stacionární náhodná posloupnost, mějme t ∈ Z a k ∈ N a zaveďme vektory uk−1,t+k = Xt+k − αk−1,1 Xt+k−1 − αk−1,2 Xt+k−2 − . . . − αk−1,k−1 Xt+1 vk−1,t = Xt − βk−1,1 Xt+1 − βk−1,2 Xt+2 − . . . − βk−1,k−1 Xt+k−1 , kde αk−1,s , βk−1,s pro s = 1, 2, . . . , k−1 jsou matice minimalizující E(uk−1,t+k )2 , resp. E(vk−1,t )2 . Dále označíme Vu (k) = var(uk−1,t+k ) Vv (k) = var(vk−1,t ) Vvu (k) = cov(uk−1,t+k , vk−1,t ). 24
Označme Dv (k) diagonální matici, která má na diagonále odmocniny příslušných diagonálních prvků matice Vv (k), a Du (k) diagonální matici analogicky odvozenou z Vu (k). Pak se funkce ξ(k) = (Dv (k))−1 Vvu (k)(Du (k))−1 definovaná pro k = 1, 2, . . . nazývá parciální autokorelační maticová funkce náhodné posloupnosti {Xt }. ¤ Parciální autokorelační maticová funkce ξ(k) posloupnosti {Xt , t ∈ Z} je mnohorozměrnou analogií parciálního autokorelačního koeficientu, a to v tom smyslu, že parciální autokorelační maticová funkce v bodě k je korelační maticí mezi Xt a Xt+k při pevných hodnotách Xt+1 , . . . , Xt+k−1 .
3.2
Mnohorozměrný ARMA model
Definice 3.7: Bílý šum Bílý šum {Yt , t ∈ Z} je náhodná posloupnost, EYt = 0 pro všechna t ∈ Z, autokovarianční maticová funkce R(s, t) je nulová matice pro všechna s 6= t, s, t ∈ Z a R(t, t) je konečná, symetrická, pozitivně definitní matice a je stejná pro všechna t ∈ Z. ¤ Úmluva 3.8: Pro bílý šum se někdy zavádí označení W N (0, Σ) (z anglického white noise), kde Σ značí varianční matici bílého šumu R(0). ¤ Definice 3.9: Operátor zpětného posunutí Nechť {Xt , t ∈ Z} je systém m-rozměrných náhodných vektorů. Pak se operátor B, který pro všechna t ∈ Z vektoru Xt přirazuje vektor Xt−1 , nazývá operátor zpětného posunutí. Pro n = 1, 2, . . . definujeme rekurentně operátor n-násobného zpětného posunutí: Nejprve položíme B 1 Xt = BXt = Xt−1 , dále pak B n Xt = B(B n−1 Xt ) = . . . = Xt−n . ¤ Definice 3.10: Klouzavé součty Náhodná centrovaná m-rozměrná posloupnost {Xt , t ∈ Z} definovaná předpisem Xt = Yt + Θ1 Yt−1 + Θ2 Yt−2 + . . . + Θq Yt−q , t ∈ Z kde {Yt , t ∈ Z} je m-rozměrný bílý šum a Θ1 , Θ2 , . . . , Θq jsou matice typu m × m, přičemž Θq je nenulová, se nazývá posloupnost klouzavých součtů řádu q, značíme M A(q). ¤ P Označíme-li Θ(B) = qi=0 Θi B i , kde B je operátor zpětného posunutí a Θ0 = I je jednotková matice, lze ekvivalentně psát Xt = Θ(B)Yt . 25
Definice 3.11: Autoregresní posloupnost Náhodná centrovaná m-rozměrná posloupnost {Xt , t ∈ Z} se nazývá autoregresní posloupnost řádu p (značíme AR(p)), jestliže vyhovuje rovnici Xt = Yt − Φ1 Xt−1 − Φ2 Xt−2 − . . . − Φp Xt−p ,
t ∈ Z,
kde Φ1 , Φ2 , . . . , Φp jsou matice typu m × m, přičemž Φp je nenulová a {Yt , t ∈ Z} je m-rozměrný bílý šum. ¤ Pp Označíme-li Φ(B) = i=0 Φi B i , kde Φ0 = I je jednotková matice, lze ekvivaletně psát Φ(B)Xt = Yt . Definice 3.12: Posloupnost ARMA Náhodná centrovaná posloupnost {Xt , t ∈ Z} se řídí modelem ARM A(p, q), jestliže Φ(B)Xt = Θ(B)Yt , kde p, q, Θ0 , Θ1 , . . . , Θq , Φ0 , Φ1 , . . . , Φp , Θ(B) a Φ(B) jsou jako v předchozích definicích. ¤ Model ARM A(p, q) je tedy smíšený model autoregrese a klouzavých součtů. Odtud plyne i označení a je zřejmé, že posloupnosti AR(p) a M A(q) lze považovat za jeho speciální případy.
3.2.1
Stacionarita a invertibilita posloupností ARMA
Definice 3.13: Sčitatelnost v absolutní hodnotě Posloupnost matic {Ψs , s ∈ NP 0 }, kde Ψs = [ψij,s ]i,j=1,2,...,m se nazývá sčitatelná v absolutní hodnotě, jestliže ∞ ¤ s=0 |ψij,s | < ∞ pro všechna i a j. Definice 3.14: Stacionarita Nechť {Xt } je posloupnost ARM A(p, q). Posloupnost {Xt } se nazývá stacionární, jestliže existuje posloupnost matic {Ψs , s ∈ N0 } sčitatelná v absolutní P∞ hodnotě a taková, že Xt = s=0 Ψs Yt−s . ¤ Nejjednodušším příkladem stacionární posloupnosti je libovolná posloupnost M A(q). Věta 3.15: Stacionární posloupnost ARM A(p, q) je centrovaná slabě stacionární náhodná posloupnost. ¤ Věta 3.16: Podmínka stacionarity Posloupnost ARM A(p, q) daná rovnicí Φ(B)Xt = Θ(B)Yt je stacionární, pokud kořeny polynomu det(Φ(B)) leží vně jednotkového kruhu. ¤ 26
Definice 3.17: Invertibilita Nechť {Xt } je posloupnost ARM A(p, q). Posloupnost {Xt } se nazývá invertibilní, jestliže existuje posloupnost matic {Πs , s ∈ N0 } sčitatelná v absolutní P∞ hodnotě a taková, že Yt = s=0 Πs Xt−s . ¤ Věta 3.18: Podmínka invertibility Posloupnost ARM A(p, q) s konečnou kovarianční maticovou funkcí daná rovnicí Φ(B)Xt = Θ(B)Yt je invertibilní, pokud kořeny polynomu det(Θ(B)) leží vně jednotkového kruhu. ¤
3.3
Identifikace ARMA modelu
Identifikace modelu je první fází výstavby modelu a jejím úkolem je rozhodnout, jaký typ modelu vybrat (tj. zda AR, M A nebo ARM A) a explicitně určit řád modelu. Výsledkem identifikační fáze je tedy např. stanovení modelu AR(2) jako vhodného modelu pro modelování předložené časové řady. Určení řádů p a q modelu je založeno na zkoumání výběrové autokorelační maticové funkce a výběrové parciální autokorelační maticové funkce. Obvykle stačí použít jen několik prvních hodnot těchto funkcí, například 10 až 20, přitom by se jich nemělo používat víc než N4 . V případě nejasné volby mezi několika modely je vhodné pokračovat se všemi modely a pro ten správný se rozhodnout až ve fázi verifikace modelu. Dále budeme pojmem časová řada rozumět řadu pozorovaných hodnot, přičemž hodnotu v čase t budeme značit xt . Pro teoretický popis budeme i nadále používat pojem náhodná posloupnost a příslušné náhodné vektory značit Xt . Složky těchto jednotlivých vektorů budeme značit xt,i , resp. Xt,i pro i = 1, 2, . . . , m. Rozsah výběru xt budeme označovat N , tzn. předpokládáme, že známe vektory xt pro t = 1, 2, . . . , N .
3.3.1
Výběrová autokovarianční a autokorelační matice
Úmluva 3.19: V dalším textu bude x¯i =
N 1 X xt,i , N t=1
i = 1, 2, . . . , m
značit odhad i-té složky vektoru středních hodnot časové řady.
27
¤
Definice 3.20: Výběrová autokovarianční matice Nechť časová řada {xt , t = 1, 2, . . . , N }, N ∈ N je realizací náhodné posloupˆ nosti {Xt , t ∈ N}. Pak se matice R(k) s prvky N −k 1 X ˆ Rij (k) = (xt,i − x¯i )(xt+k,j − x¯j ) N t=1
nazývá výběrová autokovarianční matice.
¤
Definice 3.21: Výběrová autokorelační matice Nechť časová řada {xt , t = 1, 2, . . . , N }, N ∈ N je realizací náhodné posloupnosti {Xt , t ∈ N}. Pak se matice rˆ(k) s prvky PN −k
(xt,i − x¯i )(xt+k,j − x¯j ) rˆij (k) = qP t=1 P N ¯ i )2 N ¯j )2 t=1 (xt,i − x t=1 (xt,j − x nazývá výběrová korelační matice.
¤
Autokovarianční maticová funkce je vhodným nástrojem pro určení řádu q, chceme-li popsat pozorovanou časovou řadou modelem M A(q). Věta 3.22: Autokovarianční maticová funkce R(k) posloupnosti M A(q) je nulová pro všechna |k| > q. ¤ Řád q se tedy volí takový, jaké je poslední k, pro které je ještě výběrová autokovarianční ˆ matice R(k) statisticky významně nenulová. Věta 3.23: Asymptotické vlastnosti výběrové autokorelační matice Pokud pro prvky autokorelační maticové fuknce platí rij (k) = 0 pro |k| > q pro nějaké q, pak pro k > q a pro velká N má rˆij (k) asymptoticky normální rozdělení s nulovou střední hodnotou a rozptylem k ´ X 1 ³ σij (k) = 1+2 rii (s)rjj (s) . N −k s=1 2
Navíc pro k > q má statistika m X m ³ X rˆij (k) ´2 i=1 j=1
σij (k)
asymptoticky χ2 rozdělení o m2 stupních volnosti.
28
¤
Test nulovosti prvků autokorelační matice S využitím věty 3.23 můžeme testovat hypotézu o nulovosti jednotlivých autokorelačních koeficientů. Testujeme nulovou hypotézu, že rij (k) = 0 pro |k| > q pro zvolené q, proti alternativě rij (k) 6= 0 pro nějaké |k| > q. Pak vzhledem k asymptoticky normálnímu rozdělení rˆij (k) zamítáme na hladině α nulovou hypotézu ve prospěch alternativy, pokud |ˆ rij (k)| > u1− α2 σ ˆij (k) pro některé k, kde uβ je β-kvantil standardního normálního rozdělení a k ³ ´1/2 X 1 σ ˆij (k) = √ 1+2 rˆii (s)ˆ rjj (s) . N −k s=1
Vzhledem k tomu, že rij (k) = rij (−k) pro k 6= 0, stačí rˆij (k) uvažovat pro kladná k.
Test nulovosti autokorelační matice S využitím věty 3.23 lze také testovat nulovost celé výběrové autokorelační matice. Její nulovost zamítáme, pokud m X m ³ X rˆij (k) ´2 i=1 j=1
σ ˆij (k)
> χ2m2 (1 − α),
kde χ2m2 (β) je β-kvantil χ2 rozdělení o m2 stupních volnosti.
3.3.2
Výběrová parciální autokorelační matice
Definice 3.24: Výběrová parciální autokorelační matice Nechť časová řada {xt , t = 1, 2, . . . , N }, N ∈ N je realizací náhodné poˆ sloupnosti {Xt , t ∈ N}. Pak se matice ξ(k) vypočtená podle následujícího rekurentního algoritmu nazývá výběrová parciální autokorelační matice. ¤ V [10] na str. 359-361 je uvedený algoritmus výpočtu výběrové parciání autokorelační matice. Pro k = 1 je ˆ Vu (1) = Vv (1) = R(0), ˆ Vvu (1) = R(1), −1 ˆ T (R(0)) ˆ α1,1 = R(1) , −1 ˆ ˆ β1,1 = R(1)( R(0)) . Dále pro k ≥ 2 a s = 1, 2, . . . , k − 1 je ˆ Vu (k) = R(0) − ˆ Vv (k) = R(0) −
k−1 X s=1 k−1 X s=1
29
ˆ αk−1,s R(s), ˆ T, βk−1,s R(s)
ˆ Vvu (k) = R(k) −
k−1 X
ˆ − s)αT , R(k k−1,s
s=1
αk,k αk,s βk,k βk,s
= = = =
Vvu (k)T (Vv (k))−1 , αk−1,s − αk,k βk−1,k−s , Vvu (k)(Vu (k))−1 , βk−1,s − βk,k βk−1,k−s .
Nechť Dv (k) a Du (k) jsou jako v definici 3.6. Pak ˆ = (Dv (k))−1 Vvu (k)(Du (k))−1 . ξ(k) Parciální autokorelační funkce je vhodným nástrojem pro určení řádu p, chceme-li popsat pozorovanou časovou řadou modelem AR(p). Věta 3.25: Parciální autokorelační maticová funkce ξ(k) posloupnosti AR(p) je nulová pro k > p. ¤ Řád p se tedy volí takový, jaké je poslední k, pro které je ještě výběrová parciální autoˆ kovarianční matice ξ(k) statisticky významně nenulová. Věta 3.26: Asymptotické vlastnosti výběrové parciální autokorelační matice Nechť časová řada {xt , t = 1, 2, . . . , N }, N ∈ N je realizací posloupnosti AR(p). Pak pro k > p a pro velká N má ξˆij (k) asymptoticky normální rozdělení s nulovou střední hodnotou a rozptylem N 1−k . Navíc pro k > p má statistika m X m ³ ´2 X ξˆij (k) (N − k) i=1 j=1
asymptoticky χ2 rozdělení o m2 stupních volnosti.
¤
Test nulovosti prvků parciální autokorelační matice Test nulovosti jednotlivých prvků parciální autokorelační matice využívá jejich asymptoticky normálního rozdělení. Nulovost zamítáme na hladině α, pokud |ξˆij (k)| > u1− α2 √
1 . n−k
Test nulovosti parciální autokorelační matice Test nulovosti celé parciální autokorelační matice je také založen na větě 3.26. Nulovost zamítáme na hladině α, pokud m X m ³ ´2 X (N − k) ξˆij (k) > χ2m2 (1 − α). i=1 j=1
30
3.3.3
Shrnutí volby řádu modelu
Přestože výpočet výběrové autokorelační matice a výběrové parciální autokorelační matice je k odhadu řádu příslušné posloupnosti nezbytný, univerzální návod k volbě řádů p a q obvykle nedává. Zřejmým nedostatkem výše popsaných metod je fakt, že nedávají návod, jak volit řády smíšeného modelu ARM A(p, q). Navíc i v případě, že se podaří zvolit konkrétní model AR, M A nebo ARM A, jedná se do značné míry o subjektivní odhad. Volba řádů smíšeného modelu ARMA Jednou z možností vhodnou zejména pro mnohorozměrné časové řady s méně složkami je zkoumat korelační strukturu po složkách a přistupovat k ní jako v případě jednorozměrných časových řad. Proto tu zapíšeme shrnutí poznatků o tvaru autokorelační a parciální autokorelační funkce, jak je uvedena v [3] na straně 120.
r(k)
AR(p) neexistuje k0 , r(k) ve tvaru U
M A(q) k0 = q
ξ(k)
k0 = p
neexistuje k0 , ξ(k) omezená křivkou U
ARM A(p, q) neexistuje k0 , r(k) ve tvaru U po prvních q − p krocích neexistuje k0 , ξ(k) omezená křivkou U po prvních p − q krocích
Tabulka 3.27: Tvar autokorelační a parciální autokorelační funkce procesů AR(p), M A(q) a ARM A(p, q) V tabulce uvedený symbol U označuje křivku ve tvaru lineární kombinace klesajících geometrických posloupností a sinusoid s geometricky klesající amplitudou. Dále k0 je takové, že pro |k| > k0 je r(k) = 0, resp. ξ(k) = 0. Za p a q pro zkoumanou mnohorozměrnou řadu pak volíme maximální hodnoty z p a q nalezených pro jednotlivé složky. Další možností v případě, že se nepodařilo zvolit žádný AR ani M A model, je zvolit například model ARM A(1, 1), nebo ARM A(1, 2) či ARM A(2, 1). Pro tento model by se pak provedl odhad parametrů a verifikace modelu. Pokud by byly zvolené řády příliš vysoké, projeví se to při odhadech parametrů. Volba příliš nízkých řádů se odhalí ve fázi verifikace. Jiným možným postupem při identifikaci modelu jsou statistická rozhodovací kritéria. Statistická rozhodovací kritéria Nejpoužívanějším kritériem je Akaikeho informační kritérium, které má tvar ³ ´ 2 ˆ p,q ) + 2m (p + q) , AIC(p, q) = log det(Σ N
(3.1)
ˆ p,q příslušný odhad varianční matice bílého šumu. Za kde p a q jsou řády modelu a Σ optimální model by pak měl být považovám ten, který AIC minimalizuje.
31
Z dalších kritérií připomeňme Baysovské informační kritérium, které má tvar ³ ´ ˆ p,q ) + 2m2 (p + q) log N BIC(p, q) = log det(Σ N a se kterým se zachází stejně jako s AIC.
3.3.4
Odhad parametrů modelu
V případě, že jsou již určeny řády modelu, je možno přistoupit k odhadu parametrů Θ1 , Θ2 , . . . , Θq , Φ1 , Φ2 , . . . , Φp a Σ modelu ARM A(p, q). V praxi se nejčastěji používá odhad momentovou metodou a metodou maximální věrohodnosti. Podrobněji se těmito metodami nebudeme zabývat, více o tom např. v [10], kapitola 7 nebo [6] strany 29-30.
3.3.5
Verifikace modelu
Verifikace modelu spočívá v analýze reziduí ˆ 1 xt−1 + . . . + Φ ˆ p xt−p − Θ ˆ 1 yˆt−1 − . . . − Θ ˆ q yˆt−p . yˆt = xt + Φ Pokud je model zvolen správně, rezidua yˆt budou realizacemi mnohorozměrného bílého šumu, tudíž by měla mít nulovou střední hodnotu a nebýt navzájem korelovaná. Jednotlivé složky vektoru yˆt však korelované být mohou. K dispozici máme pozorování xt pro t = 1, 2, . . . , N , proto je možné vypočítat rezidua yˆt pouze pro t = r + 1, r + 2, . . . , N , kde r = max(p, q).
Test nulovosti střední hodnoty reziduí Test nulovosti reziduí je založen na skutečnosti, že za nulové hypotézy jsou rezidua yˆt realizacemi bílého šumu W N (0, Σ). Předpokládáme také nekorelovanost pro různá t. ˆ rozptylové matice Σ (jejich prvky označme σ Máme odhad Σ ˆij , resp. σij ). Za předpokladu normality a za nulové hypotézy mají výběrové průměry N X 1 y¯i = yˆt,i N − r i=r+1
i = 1, 2, . . . , m
asymptoticky normální rozdělení s nulovou střední hodnotou a rozptylem střední hodnoty reziduí zamítáme na hladině α, když pro nějaké i je r r σ ˆii σ ˆii , resp. |¯ yi | > 2 pro α = 0, 05. |¯ yi | > u1− α2 N −r N −r
σii . N −r
Nulovost
Test normality reziduí Testování normality reziduí je do značné míry subjektivním procesem. Je možné například posuzovat histogramy jednotlivých složek reziduí, porovnávat výběrovou šikmost 32
a špičatost s příslušnými ukazateli normálního rozdělení. Test nekorelovanosti reziduí Testování nekorelovanosti reziduí se provádí pomocí výběrové autokorelační matice řady {ˆ yt }. Za nulové hypotézy jsou rezidua realizací bílého šumu, jehož autokorelační maticová funkce r(k) je nulová pro k 6= 0. Z věty 3.23 pak plyne, že rˆij (k) mají asymptoticky normální rozdělení s nulovou střední hodnotou a rozptylem N 1−k . Nulovost rij (k) zamítáme na hladině α, pokud |ˆ rij (k)| > u1− α2 √
1 2 , resp. |ˆ rij (k)| > √ pro α = 0, 05. N −k N −k
33
Kapitola 4 Grafické modely a vektorová autoregrese Tato kapitola popisuje identifikaci a odhad strukturální vektorové autoregrese s použitím metodologie grafických modelů. Citujeme zde pouze z článku [9].
4.1
Vektorová autoregrese
Definice 4.1: Vektorová autoregrese Kanonický vektorový autoregresní model p-tého řádu, AR(p), stacionární mrozměrné časové řady Xt = (Xt,1 , Xt,2 , . . . , Xt,m )T má tvar: Xt = c − Φ1 Xt−1 − Φ2 Xt−2 − . . . − Φp Xt−p + et ,
(4.1)
kde et je mnohorozměrný bílý šum s varianční maticí V . Kostanta c je pro centrovaný proces nulová. ¤ V dalším budeme předpokládat, že jde o centrovanou Gaussovskou časovou řadu a že et jsou nezávislé, stejně rozdělené náhodné vektory. Řád autoregrese p může být určen různými metodami, včetně zkoumání výběrové parciální autokorelační matice, nebo minimalizací některého rozhodovacího kritéria, např. AIC, jak bylo zmíněno v kapitole 3. Kanonický model však neodhalí vztahy mezi složkami vektoru Xt . Proto se uvažuje strukturální model AR(p) ve tvaru, jak je uveden v (4.1), ale s přidáním závislostí mezi promněnnými v čase t pomocí matice Φ0 : Φ0 Xt = d − Φ∗1 Xt−1 − Φ∗2 Xt−2 − . . . − Φ∗p Xt−p + at .
(4.2)
Vztah mezi tvarem (4.1) a (4.2) je dán rovnicemi Φ∗i = Φ0 Φi a Φ0 et = at . V centrované řadě je d = 0. Ve vztahu (4.2) je požadováno, aby varianční matice D vektoru at byla diagonální. Další podmínka na Φ0 je, že reprezentuje závislost každé složky 34
vektoru Xt na ostatních současných složkách. To je ekvivalentní s existencí permutace složek vektoru Xt tak, aby matice Φ0 byla trojúhelníková s jedničkami na diagonále. Každé možné přeuspořádání složek Xt proto může dávat odlišný tvar (4.2), ale všechny jsou statisticky ekvivalentní, odpovídají vztahu: V −1 = ΦT0 D−1 Φ0 . Význam (4.2) spočívá v možnosti, že existuje vztah, který obsahuje méně parametrů než buď (4.1) nebo jiný tvar (4.2). To je zohledněno v možnosti vyloučit mnohé prvky Φ0 a Φ∗i . V dalším textu se budeme zabývat identifikací strukturálního modelu s využitím teorie grafických modelů
4.2
Identifikace strukturálních autoregresí
Předpokládejme, že k dané časové řadě máme určený řád p kanonického autoregresního modelu AR(p). Identifikace strukturální vektorové autoregrese začíná sestrojením CIG. Využijeme při tom matici dat X typu (N − p) × (p + 1)m tvořenou sloupci (Xp+1−k,i , . . . , XN −k,i )T pro každé i = 1, . . . , m a k = 0, . . . , p. Předpokládáme, že časová řada je centrovaná, proto pro výpočet výběrové varianční matice můžeme použít vztah (1.2), kde N nahradíme N − p. Podle postupu uvedeného v kapitole 2.4 zkonstruujeme CIG s tím rozdílem, že nás zajímají pouze hrany mezi proměnnými v čase t a hrany směřující ze zpožděných proměnných do přítomnosti. Pro stacionární AR(p) model CIG sestávající pouze z výše specifikovaných hran se nezmění, pokud maximální zpoždění použité při konstrukci grafu je větší než řád p. Ovšem může zde dojít ke snížení efektivnosti v statistickém usuzování. K nalezenému CIG určíme DAG. Orientaci hran volíme z minulosti do přítomnosti. Problémem je volba orientace hran mezi současnými proměnnými v DAG. Vycházet by se mělo z logických souvislostí, přičemž se snažíme zabránit nutnosti moralizace. DAG určí regresní rovnice, u nichž odhadneme parametry. Testování nulovosti regresních koeficientů naznačí možné zjednodušení modelu. Kvalitu modelu pak posuzujeme pomocí deviance X Dev = N 0 log σ ˆr2 , kde N 0 = N −p a σ ˆr2 jsou odhady reziduálních rozptylů v jednotlivých regresních rovnicích. Jako příklad uvažujme strukturální AR(1) model pro řadu Xt = (xt , yt , zt )T . Φ0 (xt , yt , zt )T = Φ∗1 (xt−1 , yt−1 , zt−1 )T , kde matice Φ0 = [φij0 ]i,j=1,2,3 1 0 −φ210 1 0 −φ320
a Φ∗1 = [φ∗ij1 ]i,j=1,2,3 jsou uvedené ve vztahu (4.3). ∗ xt−1 0 φ111 0 φ∗131 xt 0 yt = 0 φ∗221 0 yt−1 zt−1 0 0 φ∗331 zt 1 35
(4.3)
Rozepsáním po složkách pak získáváme rovnice: xt = φ∗111 xt−1 + φ∗131 zt−1 yt = φ210 xt + φ∗221 yt−1 zt = φ320 yt + φ∗331 zt−1 Příslušný DAG je na obrázku 4.2(a). Orientace hran v něm byla zavedena následujícím způsobem. Situace dovoluje 3 možné orientace hran v DAG, xt → yt → zt , xt ← yt → zt a xt ← yt ← zt , se 4 možností xt → yt ← zt vyloučenou kvůli nutnosti přidat morální hranu mezi xt a zt v CIG. Hrana xt−1 → xt určitě není morální hrana. To implikuje orientaci xt → yt , jinak by musela být přidána morální hrana mezi xt−1 a yt . Pak xt → yt → zt je jediná možná orientace hran mezi současnými proměnnými. CIG, který jsme obdrželi moralizací DAG, je na obrázku 4.2(b). Poznamenejme ještě, že hranu mezi xt−1 a zt−1 , která by měla být při moralizaci přidána, nebudeme uvažovat. Finální DAG je na obrázku 4.2(c).
º· xt ¾ ¹¸ AK A º· ? A yt ¾ AA ¹¸A
º· º· º· º· º· xt−1 xt−1 xt−1 xt xt ¾ ¹¸ ¹¸ ¹¸ ¹¸ ¹¸ I AA@ AK@ @ @ A @ A @ A º· º· ? A @º· º· @º· A A ¾ yt−1 yt−1 yt−1 yt yt A A ¹¸ ¹¸A ¹¸ ¹¸A ¹¸ @ I @ A @ A @ A A @A @A º· º· º· @º· º· ? ? A A A @º· zt−1 zt−1 zt−1 zt ¾ zt zt ¾ ¹¸ ¹¸ ¹¸ ¹¸ ¹¸ ¹¸ (a) (c) (b)
Obrázek 4.2: (a) DAG reprezentující model (4.3), (b) morální graf odvozený od (a), (c) DAG reprezentující model (b) Orientované hrany xt−1 → xt , yt−1 → yt a zt−1 → zt nemohou být interpretovány jako morální hrany. Každá ze zbývajících hran yt−1 → xt , zt−1 → xt a zt−1 → yt může být takto interpretována.
Odhady parametrů ve vybraném modelu jsou konstruovány na základě jednotlivých regresních rovnic vytvořených podle DAG z obrázku 4.2(c). V našem případě jsou to rovnice: (1)
(1)
(1)
xt = β1X xt−1 + β1Y yt−1 + β1Z zt−1 (2)
(2)
(3)
(3)
(2)
yt = β0X xt + β1Y yt−1 + β1Z zt−1 zt = β0Y yt + β1Z zt−1 Konkrétně ukážeme popsaný přístup včetně odhadů a testování nulovosti v regresních rovnicích na zpracování finančních dat v následující kapitole. 36
Kapitola 5 Zpracování finančních dat V této kapitole aplikujeme dosud popsanou teorii na konkrétní mnohorozměrnou časovou řadu. Nejprve se budeme zabývat identifikováním mnohorozměrného kanonického ARM A modelu. Pro další účely nám postačí jenom určení řádu modelu, explicitní odhad parametrů nebudeme potřebovat (jen pro výpočty ve fázi verifikace). K samotné identifikaci kanonického modelu využijeme část programu uvedeného v [6], který jsme doplnili o identifikaci strukturálních modelů. A v dalším budeme pokračovat s grafickými modely. Všechny potřebné výpočty budou prováděné v systému Mathematica 5.0, část programu na identifikaci ARM A modelu vyžaduje speciální knihovnu Time Series Pack. Pro práci v této knihovně je užitečná příručka [5]. Pracovat budeme s dvourozměrnou časovou řadou tvořenou měsíčními průměrnými kurzy CZK/EUR a CZK/USD od ledna 1999 (zavedení eura) do srpna 2007. Data pocházejí z internetových stránek ČNB [12]. Celkem máme k dispozici 104 pozorování vektoru (CZK/EUR, CZK/USD)T . Složky vektoru jsou znázorněné na obrázku 5.1.
Obrázek 5.1: Měsíční průměrné kurzy CZK/EUR a CZK/USD Značení 5.2: Měsíční průměrné kurzy CZK/EUR (později pak transformované) budeme značit yt , t = 1, ..., 104 (po transformaci t = 1, . . . , 103). Měsíční průměrné kurzy CZK/USD (později pak transformované) budeme značit zt , t = 1, ..., 104 (po transformaci t = 1, . . . , 103) Dvourozměrný vektor (yt , zt )T budeme značit Xt pro t = 1, . . . , 103. ¤
37
5.1
Identifikace kanonického ARMA modelu
Na základě věty 3.15 musíme zabezpečit, aby daná časová řada byla slabě stacionární a centrovaná. Z obrázku 5.1 je vidět, že rozptyl řady je větší, když kurz nabývá vyšších hodnot. Proto jako transformaci použijeme zlogaritmování řady. Výsledek je graficky znázorněn na obrázku 5.3.
Obrázek 5.3: Grafy zlogaritmovaných řad Dále je potřeba zabezpečit, aby řada byla centrovaná. Z tohoto důvodu přejdeme k řadě prvních diferencí. Graf prvních diferencí je znázorněn na obrázku 5.4. Poznamenejme ještě, že uvedenou transformací se zkrátila časová řada o jednu hodnotu.
Obrázek 5.4: Grafy prvních diferencí logaritmovaných řad Jak bylo uvedeno v kapitole 3.3, vhodným nástrojem pro určení řádu p a q modelu ARM A(p, q) je zkoumání výběrové autokorelační matice rˆij (k) a výběrové parciální autokorelační matice ξˆij (k), postačí pro k = 1, . . . 20. Na obrázku 5.5 jsou uvedené prvky výběrové autokorelační matice. Čárkovanými horizontálními čarami jsou znázorněné odhadnuté směrodatné odchylky. Test nulovosti jednotlivých prvků na hladině α = 0.05 je založen na porovnání právě se směrodatnou odchylkou. Test je tedy přísnější než postup uvedený v kapitole 3.3 (tam se absolutní hodnota autokorelací porovnává se směrodatnou odchylkou vynásobenou kvantilem normálního rozdělení u1− α2 , který je při volbě α = 0.05 roven 1.96. Někdy se 1.96 může zaokrouhlit na 2). Statisticky nevýznamná korelace podle našeho testu tedy tím spíše leží v pásu ±u1− α2 × směrodatná odchylka.
38
Obrázek 5.5: Grafy prvků autokorelační matice Z obrázku 5.5 je vidět, že žádný prvek se nejeví jako statisticky významně nenulový. Pouze prvek rˆ22 (1) nepatrně překročí svou směrodatnou odchylku. Uvedené podle věty 3.22 poukazuje spíš na AR model, ale vzhledem k prvku rˆ22 (1) můžeme připustit i model M A(1). Dále se budeme zabývat výběrovou parciální autokorelační maticí. Grafické znázornění je na obrázku 5.6.
Obrázek 5.6: Grafy prvků parciální autokorelační matice Tady je situace komplikovanější. Graf ξˆ22 (k) indikuje AR(1) nebo AR(2) model. Zbylé grafy poukazují spíš na jiný než AR model. Zkoumané výběrové matice nevedou k jedinému modelu, proto dál budeme uvažovat více možností. A sice modely M A(1), AR(1) a AR(2). Také vezmeme v úvahu smíšené modely ARM A(1, 1) a ARM A(2, 1). 39
Pro první porovnání zmíněných modelů použijeme Akaikeho informační kritérium, jak bylo zavedené v (3.1). Výsledky jsou v tabulce 5.7. AIC
AR(1) -16.5707
AR(2) -16.5306
M A(1) -16.3461
ARM A(1, 1) ARM A(2, 1) -16.551 -16.6061
Tabulka 5.7: Akaikeho informační kritérium Na základě tabulky 5.7 se jako nejlepší jeví model ARM A(2, 1) s nejmenší hodnotou AIC, následuje model AR(1) a ARM A(1, 1). Podíváme se i na model AR(2), který připouští hodnoty parciální autokorelační matice. Modely pak porovnáme ve fázi verifikace, která je založená na analýze reziduí. K tomuto účelu musíme odhadnout i parametry jednotlivých modelů. Výsledky jsou uvedené v tabulce 5.8 a 5.9. Koeficienty, které nepřesahovaly odhad své směrodatné odchylky, byly nahrazené nulami.
ˆ1 Φ
AR(1) -0.215021 0.0553711 0 -0.234298
ˆ2 Φ ˆ Σ
0.000119 0.000144 0.000144 0.000668
AR(2) -0.151075 0 0 -0.295504 -0.155152 0.0911654 -0.338896 0.229181 0.000116 0.000135 0.000135 0.000646
Tabulka 5.8: Odhady parametrů v modelech AR(1) a AR(2) ˆ1 Φ
ARM A(1, 1) -0.76746 0.329646 -1.0784 0
ˆ2 Φ ˆ1 Θ ˆ Σ
0.650639 1.3924 0.000114 0.000132
-0.301586 -0.485695 0.000132 0.000638
ARM A(2, 1) -1.48564 0.45628 0 0 0 0 0 0 1.44629 -0.426014 1.27638 0 0.000994 0.000118 0.000118 0.000629
Tabulka 5.9: Odhady parametrů v modelech ARMA(1,1) a ARMA(2,1) Nyní můžeme přejít k verifikaci modelu a ověřit tak správnost volby řádu modelu. V tabulce 5.10 jsou uvedené průměry jednotlivých složek reziduí pro uvažované modely a zároveň je tam kritická hodnota pro test nulovosti na hladině α = 0.05. Je to dvojnásobek odhadnutých směrodatných odchylek.
40
e¯1 Kritická hodnota e¯2 Kritická hodnota
AR(1) 0.00212475 0.00215145 0.00303077 0.0050971
AR(2) 0.00204373 0.00212299 0.00287863 0.00501001
ARM A(1, 1) 0.00248837 0.00210813 0.00320546 0.00497948
ARM A(2, 1) -0.0102512 0.00196513 -0.0095084 0.00494418
Tabulka 5.10: Rezidua Test nulovosti splňují modely AR(1) a AR(2) v podstatě stejně. Nulovost střední hodnoty první ani druhé složky vektoru reziduí v těchto modelech nelze zamítnout. U modelů ARM A(1, 1) a ARM A(2, 1) je první složka vektoru průměru reziduí statisticky významná. Proto se těmito modely nebudeme dál zabývat. Obrázek 5.11 a 5.12 obsahují grafické znázornění jednotlivých složek reziduí u dvou uvažovaných modelů AR(1) a AR(2).
Obrázek 5.11: Rezidua modelu AR(1) Vodorovnou čárkovanou čarou jsou vyznačené kritické hodnoty pro testy nulovosti na hladině 0.05.
Obrázek 5.12: Rezidua modelu AR(2) Normalitu reziduí lze posuzovat například podle histogramu. Histogramy uvedené na obrázku 5.13 a 5.14 nejsou v zásadním rozporu s tvarem hustoty normálního rozdělení, zanedbáme-li ojedinělé odlehlé hodnoty.
41
Obrázek 5.13: Histogramy reziduí modelu AR(1)
Obrázek 5.14: Histogramy reziduí modelu AR(2) Následuje posouzení nekorelovanosti reziduí. To se provede podobně jako u samotné časové řady pomocí autokorelační matice.
Obrázek 5.15: Autokorelační matice reziduí modelu AR(1) Uveďme, že vodorovnou čárkovanou čarou jsou opět vyznačené kritické hodnoty pro testy nulovosti na hladině 0.05.
42
Obrázek 5.16: Autokorelační matice reziduí modelu AR(2) Co se týče identifikace ARM A modelu můžeme učinit závěr, že zkoumané časové řadě nejvíce vyhovují modely AR(1) a AR(2).
5.2
Identifikace strukturální autoregrese pomocí grafických modelů
Dále budeme pokračovat s oběma navrženými modely, AR(1) i AR(2). AR(2) model Začneme s modelem AR(2). Vytvoříme matici dat X, která má tvar y3 z3 y2 z2 y1 z1 X3T X2T X1T .. .. .. .. .. = .. .. .. X = ... . . . . . . . . T T T y103 z103 y102 z102 y101 z101 X103 X102 X101 a je typu (N − p) × (p + 1)m, v tomto případě 101 × 6, kde N = 103 je počet pozorování, p = 2 je řád autoregrese, m = 2 je počet složet vektoru Xt . T T Dál spočteme výběrovou varianční matici Vˆ vektoru (XtT , Xt−1 , Xt−2 ), od kterého máme N − p pozorování v řádcích matice X Vˆ =
1 XT X N −p
(5.1)
Podle důsledku 1.19 škálováním Vˆ −1 vytvoříme matici výběrových parciálních korelací. Označme prvky inverzní varianční matice Vˆ −1 = (wij )6i,j=1 . Vytvoříme si pomocnou matici D. ½ ¾ 1 1 D = Diag √ ,..., √ . w11 w66
43
Matice výběrových parciálních korelací C má pak tvar µ ¶6 wij −1 ˆ C = DV D = √ . wii wjj i,j=1 Konkrétně
C=
1 −0.494441 −0.193876 0.204407 −0.0891013 0.0994488 −0.494441 1 0.168435 −0.314146 −0.070832 0.119943 −0.193876 0.168435 1 −0.538676 −0.221704 0.25135 0.204407 −0.314146 −0.538676 1 0.133603 −0.31206 −0.0891013 −0.070832 −0.221704 0.133603 1 −0.584042 0.0994488 0.119943 0.25135 −0.31206 −0.584042 1
Rozepsáno po prvcích: −c12 −c13 −c14 −c15 −c16 −c23 −c24 −c25
= cor(yt , zt | zbytek) = cor(yt , yt−1 | zbytek) = cor(yt , zt−1 | zbytek) = cor(yt , yt−2 | zbytek) = cor(yt , zt−2 | zbytek) = cor(zt , yt−1 | zbytek) = cor(zt , zt−1 | zbytek) = cor(zt , yt−2 | zbytek)
−c26 −c34 −c35 −c36 −c45 −c46 −c56
= cor(zt , zt−2 | zbytek) = cor(yt−1 , zt−1 | zbytek) = cor(yt−1 , yt−2 | zbytek) = cor(yt−1 , zt−2 | zbytek) = cor(zt−1 , yt−2 | zbytek) = cor(zt−1 , zt−2 | zbytek) = cor(yt−2 , zt−2 | zbytek)
Nás zajímají pouze korelace veličin yt a zt s veličinami v čase t a zpožděnými, nikoli korelace zpožděných proměnných mezi sebou. Tzn. zajímají nás jenom prvky uvedené v prvních dvou řádcích nasledující tabulky. Pro ně budeme testovat nulovost. yt yt zt yt−1 zt−1 yt−2 zt−2
zt −c12
yt−1 −c13 −c23
zt−1 −c14 −c24 −c34
yt−2 −c15 −c25 −c35 −c45
zt−2 −c16 −c26 −c36 −c46 −c56
Tabulka 5.17: Koeficienty parciální korelace mezi jednotlivými složkami vektoru T T ) , Xt−2 (XtT , Xt−1 V kapitole 1.2 je uvedený test nulovosti parciálních korelačních koeficientů. Do vztahu (1.1) dosadíme N = 101 je počet pozorování vektoru (yt , zt , yt−1 , zt−1 , yt−2 , zt−2 )T , p = 4, p + 2 = 6 je počet složek vektoru (yt , zt , yt−1 , zt−1 , yt−2 , zt−2 )T . Test provádíme na hladině α = 0.05. Za r postupně volíme −c12 , −c13 , . . . , −c26 . V případě, že |T | > t1− α2 ,95 = 1.9852, zamítáme nulovost parciálního korelačního koeficientu. Nulovost jsme zamítli u −c12 , −c14 a −c24 . To znamená, že graf podmíněných nezávislostí CIG obsahuje hrany (yt , zt ), (yt , zt−1 ) a (zt , zt−1 ). Příslušný graf je na obrázku 5.18. 44
º·
º·
yt
yt−1
zt
zt−1
º·
yt−2
¹¸ ¹¸ Z Z Z º· ZZº·
¹¸
¹¸
¹¸
¹¸
º·
zt−2
Obrázek 5.18: Graf podmíněných nezávislostí pro model AR(2) Nalezení orientovaného acyklického grafu DAG k danému CIG spočívá v určení orientace hran. Mezi vrcholy reprezentujícími proměnné ve stejném čase volíme takovou orientaci, aby pokud možno nemusela být prováděná moralizace, tzn. aby do vrcholu yt nebo zt nesměrovaly dvě orientované hrany tvořící ”neoddané rodiče”, zejména, je-li jeden z rodičů zt nebo yt . Případně volíme logicky interpretovatelnou orientaci. Mezi vrcholy reprezentujícími proměnné v různých časech má smysl pouze orientace od minulosti k budoucnosti. Hrany mezi vrcholy reprezentujícími dvojice zpožděných proměnných nás zpravidla nezajímají. º·
yt
º·
yt−1
¹¸ ¹¸ Z } 6 Z Z Z º· º· ? Z zt−1 zt ¾ ¹¸ ¹¸
º·
yt−2
¹¸ º·
zt−2
¹¸
Obrázek 5.19: Orientovaný acyklický graf pro model AR(2) Z obrázku 5.19 vidíme, že v grafu DAG není potřeba moralizace ani při orientaci yt → zt ani při orientaci zt → yt . Zároveň je vidět, že se neprokázala závislost proměnných yt−2 a zt−2 na yt a zt . Proto se v dalším budeme zabývat případem AR(1). AR(1) model Matice dat v případě modelu AR(1) má tvar: y2 z2 y1 z1 X2T X1T .. .. .. = .. .. X = ... . . . . . T T y103 z103 y102 z102 X103 X102 a je typu (N − p) × (p + 1)m, teda 102 × 4, kde N = 103 je počet pozorování, p = 1 je řád autoregrese, m = 2 je počet složet vektoru Xt . T ), Dál s využitím vztahu (5.1) spočteme výběrovou varianční matici Vˆ vektoru (XtT , Xt−1 od kterého máme N − p pozorování v řádcích matice X. Podle důsledku 1.19 škálováním Vˆ −1 vytvoříme matici výběrových parciálních korelací C.
45
1 −0.511881 −0.233198 0.250707 −0.511881 1 0.124975 −0.300958 C= −0.233198 0.124975 1 −0.570427 0.250707 −0.300958 −0.570427 1 Rozepsáno po prvcích: −c12 = cor(yt , zt | zbytek) −c23 = cor(zt , yt−1 | zbytek) −c13 = cor(yt , yt−1 | zbytek) −c24 = cor(zt , zt−1 | zbytek) −c14 = cor(yt , zt−1 | zbytek) −c34 = cor(yt−1 , zt−1 | zbytek) Opět nás zajímají pouze první dva řádky tabulky 5.20. yt yt zt yt−1 zt−1
zt −c12
yt−1 −c13 −c23
zt−1 −c14 −c24 −c34
Tabulka 5.20: Koeficienty parciální korelace mezi jednotlivými složkami vektoru T (XtT , Xt−1 ) Test nulovosti parciálních korelačních koeficientů provedeme stejně jako u AR(2) modelu s využitím vztahu (1.1), ve kterém N = 102 je počet pozorování vektoru (yt , zt , yt−1 , zt−1 )T , p = 2, p + 2 = 4 je počet složek vektoru (yt , zt , yt−1 , zt−1 )T . Test provedeme opět na hladině α = 0.05. Za r postupně volíme −c12 , −c13 , . . . , −c24 . V případě, že |T | > t1− α2 ,97 = 1.9847, zamítáme nulovost parciálního korelačního koeficientu. Nulovost jsme zamítli u −c12 , −c13 , −c14 a −c24 . Graf podmíněných nezávislostí CIG tedy obsahuje hrany (yt , zt ), (yt , yt−1 ), (yt , zt−1 ) a (zt , zt−1 ), znázorněn je na obrázku 5.21. º·
º·
yt−1
yt
¹¸ ¹¸ Z Z Z º· ZZº·
zt−1
zt
¹¸
¹¸
Obrázek 5.21: Graf podmíněných nezávislostí pro model AR(1) K hraně (yt , yt−1 ), která u AR(2) modelu nebyla, poznamenejme ještě, že při testu nulovosti parciálních korelací u AR(2) modelu vyšla cor(yt , yt−1 | zbytek) = −0.193876, nulovost jsme nezamítli, ale již u cor(yt , zt−1 | zbytek) = 0.204407 jsme nulovost zamítali.
46
Uplatíme-li při konstrukci DAG stejné zásady, jaké jsme uvedli u AR(2) modelu, získáme graf na obrázku 5.22. Označme si tento model M . º· º· yt−1 yt ¾ ¹¸ ¹¸ Z } Z Z Z º· º· ? Z ¾ zt−1 zt ¹¸ ¹¸
Obrázek 5.22: Orientovaný acyklický graf pro model M Z obrázku 5.22 je vidět, že vrcholy yt−1 , zt−1 a yt nesplňují Wermuthovu podmínku z definice 2.28. Správně by měla být provedena moralizace, tj. přidání hrany (yt−1 , zt−1 ). Jak však už bylo zmíněno, hrany mezi zpožděnými vrcholy nás nezajímají, tudíž tuto hranu nemusíme uvažovat. Co se týká orientace hrany (yt , zt ), je teď situace jiná než u AR(2) modelu. Orientace zt → yt by si při moralizaci vyžádala přidání morální hrany (zt , yt−1 ). Naproti tomu orientace yt → zt žádné přidání morální hrany nevyžaduje. Na základě věty o separaci 2.14 a věty 2.34 můžeme v grafu na obrázku 5.22 identifikovat nasledující podmíněnou nezávislost: zt ⊥ yt−1 | yt , zt−1 . Nyní můžeme setavit regresní rovnice pro yt a zt . Obecně mají tvar: (1)
(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
(2)
yt = β0 + β0Z zt + β1Y yt−1 + β1Z zt−1 + et
zt = β0 + β0Y yt + β1Y yt−1 + β1Z zt−1 + et . (1)
(1)
V regresních rovnicích jsou však absolutní členy β0 a β0 nulové, neboť předpokládáme centrovanou časovou řadu. S nenulovým koeficientem β jsou v rovnicích pak pouze ty proměnné, z jejichž vrcholů vede hrana do yt a zt (jak je znázorněno na obrázku 5.22). Tj.: (1)
(1)
(1)
yt = β1Y yt−1 + β1Z zt−1 + et (2)
(2)
(2)
zt = β0Y yt + β1Z zt−1 + et . Odhad parametrů β provedeme metodou nejmenších čtverců, jak bylo uvedené ve větě 1.21. Dostáváme: yˆt = 0.215021yt−1 − 0.0553711zt−1 zˆt = 1.16182yt + 0.233692zt−1 .
47
º· º· yt ¾ 0.215 yt−1 ¹¸ ¹¸ Z } Z 1.162 Z−0.055 Z º· º· ? Z 0.234 zt−1 zt ¾ ¹¸ ¹¸
Obrázek 5.23: Orientovaný acyklický graf pro model M s vyznačenými regresními koeficienty Odhady reziduálních rozptylů, provedené podle věty 1.23 jsou s21 = 0.000121574 a s22 = 0.000511567. A deviance navrženého modelu M je Dev(M ) = (N − p)
2 X
log s2r = −1692.49.
r=1
Identifikace strukturálního modelu pomocí DAG prokázala silnou provázanost proměnných v čase t a t − 1. DAG obsahuje kromě hrany (zt , yt−1 ) všechny možné. Otestujme dále nulovost regresních koeficientů. Využijeme přitom větu 1.25. Testujeme (1) (1) hypotézu H0 : βi = 0 proti alternativě H1 : βi 6= 0. Za βi postupne volíme β1Y , β1Z , (2) (2) β0Y a β1Z . Test prováděný na hladině α = 0.05 nezamítá nulovost regresního koeficientu (1) β1Z , u ostatních koeficientů jsme nulovou hypotézu zamítli. To znamená, že hranu zt−1 → yt lze vynechat. Poznamenejme, že tuto hranu lze interpretovat jako morální. Model s vynechanou hranou si označme M 0 . Nové regresní rovnice pak mají tvar: yˆt = 0.14509yt−1 zˆt = 1.16182yt + 0.233692zt−1 . º· yt ¾ 0.145 ¹¸
º·
yt−1
¹¸
1.162
º· ? zt ¾ 0.234 ¹¸
º·
zt−1
¹¸
Obrázek 5.24: Orientovaný acyklický graf pro model M 0 s vyznačenými regresními koeficienty Odhady reziduálních rozptylů jsou s21 = 0.00012329 a s22 = 0.000511567. 48
Deviance modelu M 0 je Dev(M 0 ) = −1691.06. Minimální nárůst deviance modelu M 0 s vynechanou hranou oproti původnímu modelu M, Dev(M 0 ) − Dev(M ) = 1.43, není v rozporu s odstraněním dané hrany. Při provedené analýze 2-rozměrné časové řady kurzů CZK/EUR a CZK/USD se ukázalo ovlivnění obou kurzů v současnosti hodnotou zpožděnou o jeden měsíc. Dále se prokázala souvislost kurzů obou měn vůči české koruně ve stejném čase, která je modelem interpretována jako ovlivnění hodnoty CZK/USD hodnotou CZK/EUR. Vliv zpožděných hodnot jedné na současný kurz druhé měny se výrazně neprokázal.
49
Závěr Cílem práce bylo vyložit základní teorii mnohorozměrných časových řad, zejména tzv. ARM A modelů, ukázat jejich propojení s grafickými modely a aplikovat ji na časovou řadu měsíčních průměrů kurzů eura a dolaru vůči české koruně. V teoretické části jsme bez důkazů uváděli základy korelace, regrese, definice a vlastností ARM A modelu a teorie grafů. Propojenost jednotlivých teoretických kapitol jsme ukázali v souvislosti s grafem podmíněných nezávislostí. Na základě tohoto grafu je pak odvozený strukturální AR model. Význam strukturální autoregrese spočívá v identifikování závislostí i mezi jednotlivými složkami pozorovaného vektoru, nejen závislosti na minulých hodnotách, jak je to u kanonických modelů. Uvedený postup byl pak prakticky aplikován na napozorovanou časovou řadu. Výsledky a spůsob výpočtů je obsažen i na přiloženém CD v souboru zpracovani.nb. Výsledky získáné použitím jednotlivých metod jsou si velice blízké a navzájem se doplňují.
50
Literatura [1] Anděl, J. (1976): Statistická analýza časových řad. SNTL, Praha. [2] Anděl, J. (2007): Základy matematické statistiky. Matfyzpress, Praha. [3] Cipra, T. (1986): Analýza časových řad s aplikacemi v ekonomii. SNTL, Praha. [4] Dupač, V., Hušková M. (2005): Pravděpodobnost a matematická statistika. Karolinum, Praha. [5] He, Y. (1995): Time Series Pack - Reference and User’s Guide. Wolfram Research, Illinois. [6] Hrba, M. (2006): Aplikace modelů mnohorozměrných časových řad ve finanční analýze. Diplomová práce, MFF UK. [7] Matoušek, J., Nešetřil, J. (2002): Kapitoly z diskrétní matematiky. Karolinum, Praha. [8] Prášková, Z. (2004): Základy náhodných procesů II. Karolinum, Praha. [9] Reale, M., Wilson, G. T. (2000): Identification of vector AR models with recursive structural errors using conditional independence graphs. Statistical Methods and Applications 10(2001), 49-56. [10] Wei, W. W. S. (1990): Time Series Analysis. Addison-Wesley Publishing Company, California. [11] Whittaker, J. (1990): Graphical models in Applied multivariate Statistics. John Wiley, New York. [12] www.cnb.cz/cs/financni trhy/devizovy trh/kurzy devizoveho trhu/prumerne form.jsp
51
Přílohy
52
Příloha č.1: Soubor vstupních dat Měsíc I.99 II.99 III.99 IV.99 V.99 VI.99 VII.99 VIII.99 IX.99 X.99 XI.99 XII.99 I.00 II.00 III.00 IV.00 V.00 VI.00 VII.00 VIII.00 IX.00 X.00 XI.00 XII.00 I.01 II.01 III.01 IV.01 V.01 VI.01 VII.01 VIII.01 IX.01 X.01 XI.01 XII.01 I.02 II.02 III.02 IV.02 V.02
EUR/CZK 35.638 37.715 37.989 37.997 37.692 37.152 36.521 36.415 36.356 36.587 36.403 36.054 36.025 35.709 35.595 36.31 36.555 36.017 35.619 35.356 35.425 35.275 34.617 34.817 35.139 34.64 34.601 34.55 34.382 33.975 33.855 34.034 34.188 33.562 33.325 32.592 32.078 31.789 31.388 30.356 30.558
USD/CZK 30.656 33.597 34.845 35.465 35.434 35.762 35.305 34.305 34.635 34.077 35.12 35.63 35.45 36.253 36.844 38.226 40.319 37.979 37.871 39.004 40.658 41.125 40.475 38.942 37.425 37.551 37.955 38.727 39.27 39.777 39.335 37.87 37.559 37.009 37.475 36.475 36.325 36.539 35.844 34.269 33.313
Měsíc VI.02 VII.02 VIII.02 IX.02 X.02 XI.02 XII.02 I.03 II.03 III.03 IV.03 V.03 VI.03 VII.03 VIII.03 IX.03 X.03 XI.03 XII.03 I.04 II.04 III.04 IV.04 V.04 VI.04 VII.04 VIII.04 IX.04 X.04 XI.04 XII.04 I.05 II.05 III.05 IV.05 V.05 VI.05 VII.05 VIII.05 IX.05 X.05
53
EUR/CZK 30.295 29.749 30.796 30.193 30.653 30.756 31.192 31.49 31.645 31.758 31.625 31.391 31.41 31.877 32.289 32.354 31.985 31.974 32.313 32.723 32.857 32.984 32.514 31.974 31.614 31.521 31.634 31.6 31.484 31.287 30.647 30.31 29.961 29.782 30.13 30.216 30.032 30.191 29.592 29.305 29.677
USD/CZK 31.726 29.959 31.495 30.787 31.236 30.716 30.654 29.653 29.374 29.393 29.157 27.095 26.936 28.035 28.999 28.849 27.354 27.343 26.32 25.949 25.985 26.902 27.117 26.633 26.048 25.709 25.984 25.876 25.233 24.091 22.87 23.102 23.024 22.585 23.288 23.809 24.688 25.05 24.073 23.896 24.715
Měsíc XI.05 XII.05 I.06 II.06 III.06 IV.06 V.06 VI.06 VII.06 VIII.06 IX.06 X.06 XI.06 XII.06 I.07 II.07 III.07 IV.07 V.07 VI.07 VII.07 VIII.07
EUR/CZK 29.261 28.975 28.721 28.409 28.65 28.508 28.271 28.385 28.445 28.193 28.38 28.29 28.03 27.777 27.841 28.231 28.055 28.01 28.231 28.545 28.33 27.858
USD/CZK 24.82 24.439 23.733 23.796 23.834 23.251 22.138 22.441 22.441 22.007 22.297 22.433 21.754 21.022 21.419 21.593 21.189 20.731 20.899 21.272 20.641 20.45
54
Příloha č.2: Obsah přiloženého CD Přiložené CD obsahuje tyto soubory: 1. diplomova prace.pdf - elektronická verze této diplomové práce. 2. zpracovani.nb - notebook programu Mathematica, částečně převzatý z [6] a upravený pro naše potřeby. Je členěný na následující části: Načtení knihoven a dat ze souboru - načte potřebné knihovny a data určené ke zpracování. Transformace dat - provede se transformace dat - zlogaritmování a vytvoření řady prvních diferencí. Určení řádů modelu - provádí se výpočet výběrové autokorelační matice a výběrové parciální autokorelační matice. Na základě toho pokračujeme se 4 modely: AR(1), AR(2), ARM A(1, 1) a ARM A(2, 1). Odhad parametrů modelu - pro všechny modely jsou odhadnuty parametry. Verifikace modelu - verifikace všech modelů. Strukturální autoregrese v modelu - výpočty potřebné k stanovení struturálního AR(1) a AR(2) modelu. Uvedený soubor je uložený i se všemi výstupy, protože k práci je potřebná knihovna Time Series Pack, která není standardní součásti programu Mathematica. 3. data.txt - textový soubor zdrojových dat, udenený též v příloze 1.
55