Bc. Marie Turčičová. Modelování závislosti mezi hydrologickými a meteorologickými veličinami měřenými v několika stanicích

Univerzita Karlova v Praze Matematicko-fyzikální fakulta

DIPLOMOVÁ PRÁCE

Bc. Marie Turčičová Modelování závislosti mezi hydrologickými a meteorologickými veličinami měřenými v několika stanicích Katedra pravděpodobnosti a matematické statistiky

Vedoucí diplomové práce: Prof. RNDr. Daniela Jarušková CSc. Studijní program: Matematika Studijní obor: Pravděpodobnost, matematická statistika a ekonometrie Studijní plán: Matematická statistika

Praha 2012

Děkuji Prof. RNDr. Daniele Jaruškové, CSc. za vedení práce, za příjemnou spolupráci a především za užitečné rady a připomínky. Dále děkuji Prof. Rogeru Koenkerovi z University of Illinois za několik cenných rad a poskytnutí funkce rq.fit.hogg. Obrovský dík patří dále mé rodině a blízkým, kteří svou podporou přispěli k dokončení této práce.

Prohlašuji, že jsem tuto diplomovou práci vypracoval(a) samostatně a výhradně s použitím citovaných pramenů, literatury a dalších odborných zdrojů. Beru na vědomí, že se na moji práci vztahují práva a povinnosti vyplývající ze zákona č. 121/2000 Sb., autorského zákona v platném znění, zejména skutečnost, že Univerzita Karlova v Praze má právo na uzavření licenční smlouvy o užití této práce jako školního díla podle §60 odst. 1 autorského zákona. V Praze dne 13. 4. 2012

Podpis autora

Název práce: Modelování závislosti mezi hydrologickými a meteorologickými veličinami měřenými v několika stanicích Autor: Bc. Marie Turčičová Katedra: Katedra pravděpodobnosti a matematické statistiky Vedoucí diplomové práce: Prof. RNDr. Daniela Jarušková CSc., České vysoké učení technické v Praze, Stavební fakulta, katedra matematiky Abstrakt: Cílem této práce je průzkum závislosti denního průměrného průtoku řeky Opavy na vysokých denních srážkových úhrnech v jejím povodí. V práci jsou představeny tři metody, které lze použít při analýze závislosti vysokých hodnot veličin, a je předvedena jejich aplikace na studovaná data. V první řadě je to koeficient závislosti chvostů, který měří závislost vysokých hodnot dvou spojitých náhodných veličin. Konkrétní model pro vysoké kvantily průtoku při daných srážkách je určen nejprve neparametricky pomocí kvantilové regrese, ale dále také parametricky prostřednictvím metody špiček nad prahem (POT). Klíčová slova: závislost vysokých hodnot, koeficient závislosti chvostů, kvantilová regrese, metoda špiček nad prahem

Title: Modelling dependence between hydrological and meteorological variables measured on several stations Author: Bc. Marie Turčičová Department: Department of Probability and Mathematical Statistics Supervisor: Prof. RNDr. Daniela Jarušková CSc., Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics Abstract: The aim of the thesis is to explore the dependence of daily discharge averages of the Opava river on high daily precipitation values in its basin. Three methods are presented that can be used for analyzing the dependence between high values of random variables. Their application on the studied data is also given. First it is the tail-dependence coefficient that measures the dependence between high values of two continuous random variables. The model for the high quantiles of the discharge at a given precipitation value was first determined non-parametrically by quantile regression and then parametrically through the peaks-over-threshold (POT) method. Keywords: extremal dependence, tail-dependence coefficient, quantile regression, peaks over threshold method

Obsah 1 Úvod

2

2 Popis dat

4

3 Koeficienty závislosti 3.1 Korelační koeficienty . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Koeficient závislosti chvostů . . . . . . . . . . . . . . . . . . . . .

7 7 11

4 Kvantilová regrese 4.1 Kvantil a jeho odhad . . . . . . . . . . 4.2 Teoretické základy kvantilové regrese . 4.3 Hledání modelu pro zpracovávaná data 4.3.1 Výpočet v R . . . . . . . . . . . 4.4 Test rovnoběžnosti regresních kvantilů 4.5 Porovnání letního a zimního období . . 4.5.1 Výpočet v R . . . . . . . . . . .

. . . . . . .

15 15 16 18 23 23 27 28

5 Metoda špiček nad prahem 5.0.2 Výpočet v R . . . . . . . . . . . . . . . . . . . . . . . . . .

33 37

6 Diskuse výsledků

49

7 Závěr

53

Přílohy 7.1 Koeficient závislosti chvostů . 7.1.1 Letní období . . . . . . 7.1.2 Zimní období . . . . . 7.2 Rovnoběžné regresní kvantily

56 56 56 58 60

. . . .

1

. . . .

. . . .

. . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . . . . .

. . . .

. . . .

1. Úvod Modelování vztahu mezi srážkami a průtokem je jedním ze základních problémů v hydrologii. Průtok je přitom srážkami ovlivňován dvojím způsobem. Zřejmý je dopad aktuální hodnoty srážky, popřípadě součtu srážkových úhrnů za několik předchozích dní. Okamžitý průtok je ale srážkami ovlivňován i nepřímo, a to prostřednictvím veličin, jako je nasycenost půdy nebo výška spodní vody, které jsou ovlivňované množstvím srážek během dlouhé doby před okamžikem měření. Kromě srážek mají na velikost průtoku samozřejmě vliv i další veličiny, především teplota a vlhkost vzduchu, které ovlivňují vypařování vody do ovzduší. Data, která jsou předmětem zkoumání v této práci, obsahují denní průměrné průtoky řeky Opavy a denní srážkové úhrny v jejím povodí. Data pocházejí z let 1960-2003 a jejich podrobný popis lze nalézt v kapitole 2. K dispozici bohužel nebyly žádné informace o naměřených teplotách ani o jakýchkoli dalších meteorologických či hydrologických veličinách. Přes nesporný vliv těchto faktorů lze však říci, že při relativně vysokých hodnotách srážkových úhrnů je vliv srážek na okamžitý průtok poměrně značný, a má tedy cenu pokusit se závislost těchto dvou veličin hlouběji prozkoumat. Předkládaná diplomová práce se zabývá problémem, jak velké mohou být průtoky, je-li celkový srážkový úhrn za několik předchozích dní vysoký. Abychom přitom alespoň trochu zohlednili vliv ostatních veličin, rozhodli jsme se závislost zkoumat zvlášť v letním a zimním období. Z matematického hlediska se dá studovaný problém formálně zapsat tak, že nás zajímá odhad P (Yt > y|f (xt−1 , . . . , xt−d ) = s), kde Yt značí průměrný průtok ve dni t a xt−1 , . . . , xt−d jsou denní srážkové úhrny za předchozí dny. V této práci jsme pro jednoduchost předpokládali, že f je sumou těchto srážkových úhrnů. Jak již bylo naznačeno výše, hlavní pozornost byla věnována odhadu této podmíněné pravděpodobnosti pro hodnoty y a s velké, tj. snažili jsme se najít vhodný model jen pro situace, kdy součet srážkových úhrnů za předchozích d dní překročil určitou mez. Tuto mez jsme konkrétně zvolili jako 90% kvantil součtu srážkových úhrnů za předchozích d dní. První otázka, na kterou jsme hledali odpověď, se týkala počtu dní d, jejichž součet srážkových úhrnů by byl vhodným prediktorem do hledaného modelu. V tomto ohledu pro nás byly směrodatné hodnoty koeficientu závislosti chvostů, který měří závislost extrémně vysokých hodnot dvou náhodných veličin. O tomto koeficientu (a nejen o něm) je podrobně pojednáno v kapitole 3. Hlavní náplní práce je hledání vhodného modelu pro P (Yt > y|xt−1 + . . . + xt−d ). Jako hlavní nástroj jsme přitom použili kvantilovou regresi, která umožňuje nalezení regresního modelu pro libovolný kvantil odezvy. Důležitou otázkou zde bylo, zda hledat model přímo pro odezvu Yt a regresor v podobě součtu srážkových úhrnů, nebo zda na data aplikovat nějakou transformaci. Po vyřešení této otázky a nalezení uspokojivého modelu, jsme vypočetli parametry regresního mediánu a pokusili se ještě o další zjednodušení. Dále jsme provedli srovnání sklonu regresního mediánu v letním a zimním období, abychom mohli posoudit, zda se závislost průtoku a vysokého součtu srážkových úhrnů v těchto obdobích liší. Kvantilové regresi je věnována kapitola 4. Vzhledem k tomu, že náš zájem byl soustředěn hlavně na vysoké hodnoty y, použili jsme alternativně k odhadu chování vyšších kvantilů metodu špiček nad 2

prahem, kterou se zabýváme v kapitole 5. Za práh jsme zvolili 70% regresní kvantil vypočtený v kvantilové regresi. Jeho přesahy jsme v souladu s touto metodou modelovali zobecněným Paretovým rozdělením, jehož parametry jsme odhadli metodou maximální věrohodnosti. Předposlední kapitolou práce je diskuse výsledků, kde jsou shrnuty získané výsledky. Zajímavé jsou především odhady vysokých kvantilů průtoku vypočtených na základě kvantilové regrese a metody špiček nad prahem. Výsledky obou metod jsou porovnány s naivními odhady, které vznikly vypočtením výběrového kvantilu průtoku v okénku určité délky. Všechny výpočty byly provedeny ve statistickém softwaru R a příslušné skripty lze nalézt na přiloženém CD. Software R je v současné době hojně používaný, proto jsou v každé kapitole alespoň stručně zmíněny nejdůležitější příkazy, někdy doplněné podrobnějším komentářem. Tyto pasáže je samozřejmě možné přeskočit.

3

2. Popis dat Data, která jsme zpracovávali, poskytl Český hydrometeorologický ústav a obsahují údaje o průtoku severomoravské řeky Opavy a srážkových úhrnech v devíti různých místech jejího povodí. Průtok řeky byl měřen v městě Opavě jako průměrný denní průtok v m3 /s. Srážky byly sledovány na těchto srážkoměrných stanicích: 1. Heřmanovice (budeme značit zkratkou Herm), 2. Karlovice (Karl), 3. Krnov (Krno), 4. Lichnov (Lich), 5. Opava (Opav), 6. Praděd (Prad), 7. Rejvíz (Rejv), 8. Vidly (Vidl), 9. Albrechtice-Žáry (Zary). Tyto stanice pokrývají tok řeky Opavy od pramene až po místo měření průtoku. Celou situaci si lze prohlédnout na přiložené mapě. Srážkové úhrny byly měřeny standardně v milimetrech kapalné vody (1 mm = 1 l/m2 ) spadlé na zemský povrch za daný den. Všechny hodnoty byly zjišťovány denně v období 1. 1. 1960 až 31. 10. 2003, avšak ve většině záznamů o srážkových úhrnech se objevují delší časové úseky, v nichž pozorování chybí. Záznamy o průtoku jsou kompletní. Z hydrologie je známo (viz [7]), že momentální průtok není ovlivněn jen srážkami z několika posledních dní, ale velkou roli zde hrají i další meteorologické veličiny jako je vlhkost vzduchu a především teplota. Obzvlášť důležité je, zda teplota klesá nebo neklesá pod bod mrazu. Je-li teplota vzduchu větší nebo rovna 0◦ C, je třeba uvažovat jiný vztah mezi srážkami a průtoky, než v případě, že teplota klesne pod 0◦ C. Při nízkých teplotách totiž srážky zůstávají ležet na povrchu ve formě sněhu, kdežto při vyšších teplotách mají schopnost odtéci do blízkých toků. Leží-li při nezáporné teplotě vzduchu na zemi sněhová pokrývka, průtok je ovlivněn táním sněhu. Součástí našich dat bohužel nebyly žádné teplotní údaje, a proto jsme se alespoň rozhodli modelovat závislost zvlášť v letním a zimním období. Rozdělení na letní a zimní období jsme provedli na základě dlouhodobých průměrů měsíčních teplot naměřených v Ostravě: Leden Únor 0.0606 -1.7882 Červenec Srpen 18.9163 18.4628

Březen 3.3279 Září 13.7544

Duben Květen 9.0850 14.1177 Říjen Listopad 9.1921 4.6628

Červen 17.3115 Prosinec -0.6789

Je třeba poznamenat, že teplota vzduchu závisí na nadmořské výšce, která je samozřejmě u každé ze zmíněných srážkoměrných stanic jiná. Všechny jmenované stanice však leží ve vyšší nadmořské výšce než meteorologická stanice v Ostravě, a lze tedy předpokládat, že průměrné měsíční teploty v jednotlivých stanicích byly nižší než dlouhodobé průměry uvedené v tabulce výše. Na základě tohoto faktu jsme se rozhodli brát jako data odpovídající zimním podmínkám hodnoty pocházející z období od začátku listopadu do konce března. Letní podmínky jsme uvažovali pro hodnoty pocházející ze zbývajících měsíců. Jak již bylo řečeno dříve, srážková data nejsou kompletní a řada pozorování v nich chybí. Počty měření pro jednotlivé srážkoměrné stanice v létě a v zimě jsou uvedeny v tabulce níže, přičemž kompletní data by obsahovala 9416 pozorování v letním období a 6594 pozorování v zimním období. V záznamech o průtoku žádné pozorování nechybí a obsahují tedy 9416 a 6594 hodnot. 4

Herm Karl léto 8713 9386 zima 6081 6594

Krno Lich 9263 9416 6473 6594

Opav 9202 6442

Prad Rejv Vidl 8071 8682 8621 5687 6048 6050

Zary 9416 6504

Charakteristiky polohy pro srážkové úhrny z jednotlivých srážkoměrných stanic v zimním i letním období jsou uvedeny v tabulkách 2.1 a 2.2. Vidíme, že srážková data obsahují mnoho nulových hodnot. Průměr z denních srážkových úhrnů je nízký, zato letní maximální hodnoty jsou velmi vysoké. Z tabulek 2.1 a 2.2 lze tedy usuzovat, že data obsahují kromě nul velké množství malých hodnot a jen několik velmi vysokých hodnot. V zimním období jsou hodnoty maxim nižší než v létě, avšak situace je zcela analogická. Srážkové úhrny jsou tedy v letním i zimním období výrazně kladně zešikmeny. Charakteristiky polohy pro průtok jsou k nahlédnutí v tabulce 2.3. V průtokových datech se samozřejmě žádné nulové hodnoty nevyskytují, ale opět pozorujeme značné kladné zešikmení. Naším cílem je popis souvislosti vysokých hodnot srážkových úhrnů a odpovídajících průtoků, ovšem jak vidíme z příslušných tabulek, takovýchto hodnot je v datech relativně méně.

5

min Herm 0 Karl 0 Krno 0 Lich 0 Opav 0 0 Prad Rejv 0 Vidl 0 Zary 0

průměr 75% kv. 3.17 2.70 2.41 1.80 2.16 1.60 2.16 1.50 2.12 1.50 3.66 4.00 3.57 3.20 3.57 3.60 2.54 2.00

90% kv. 9.30 7.80 7.10 7.00 6.80 11.20 10.59 10.80 7.90

98% kv. 27.18 22.40 20.68 19.60 19.80 27.60 30.60 28.30 23.10

max 196.50 124.20 59.20 110.00 62.00 139.40 214.20 199.30 125.00

podíl kladných dat 41.9 % 39.5 % 39.1 % 40.2 % 40.3 % 44.5 % 41.1 % 47.2 % 47.4 %

Tabulka 2.1: Charakteristiky polohy pro srážkové úhrny z letního období.

min Herm 0 Karl 0 Krno 0 Lich 0 0 Opav Prad 0 Rejv 0 Vidl 0 Zary 0

průměr 75% kv. 1.57 1.60 1.43 1.30 0.95 0.70 1.02 0.80 0.89 0.50 2.38 3.10 1.77 1.90 2.51 2.90 1.38 1.40

90% kv. 5.10 4.70 3.00 3.10 2.60 7.20 5.60 7.70 4.40

98% kv. 12.68 12.30 9.26 9.51 8.80 15.73 14.20 18.90 10.99

max 39.2 35.7 28.0 28.0 39.4 51.8 42.1 51.8 34.6

podíl kladných dat 42.0 % 40.8 % 36.8 % 40.3 % 36.6 % 48.2 % 39.3 % 49.3 % 50.4 %

Tabulka 2.2: Charakteristiky polohy pro srážkové úhrny ze zimního období.

min léto 0.51 zima 0.90

25% kv. 2.90 2.68

medián 4.75 4.07

průměr 75% kv. 7.73 8.85 5.92 6.87

90% kv. 15.50 11.80

Tabulka 2.3: Charakteristiky polohy pro průtok.

6

98% kv. max 33.50 553.00 23.21 68.40

3. Koeficienty závislosti V této kapitole bychom si rádi učinili představu o síle závislosti průtoku a srážkových úhrnů. Je přirozené předpokládat, že vysoká hodnota průtoku není způsobena jen srážkovým úhrnem z daného dne, ale má na ni vliv i množství srážek z předchozích dní. Otázkou je, který srážkový úhrn, popř. který součet srážkových úhrnů za předchozí dny, ovlivňuje průtok nejvíce. Druhá otázka spočívá v tom, ze které srážkoměrné stanice pocházejí srážkové úhrny, které s průtokem v Opavě nejvíce souvisejí. K tomuto účelu nám poslouží různé korelační koeficienty a koeficient závislosti chvostů, jejichž výpočet provedeme v této kapitole. Zdůrazněme však, že půjde pouze o předběžnou analýzu, neboť žádný z těchto koeficientů nám neposkytuje informaci o tvaru zkoumané závislosti.

3.1

Korelační koeficienty

Nechť X a Y jsou náhodné veličiny s konečnými druhými momenty a s kladnými rozptyly. Pak Pearsonův korelační koeficient je definován vztahem cov(X, Y ) ρ= p . (var X)(var Y )

Jeho odhadem při daném náhodném výběru (X1 , Y1 )′ , . . . , (Xn , Yn )′ je výběrový korelační koeficient Pn ¯ i − Y¯ ) (Xi − X)(Y , r = pPn i=1 ¯ 2 Pn (Yi − Y¯ )2 (X − X) i i=1 i=1 ¯ a Y¯ značí výběrové průměry. Nevýhodou tohoto koeficientu je, že indikuje kde X pouze lineární závislost a je velmi citlivý na odlehlé hodnoty. Jako alternativu je proto možné použít Spearmanův korelační koeficient Pn ¯ ¯ (Ri − R)(Q i − Q) rS = pPn i=1 , P ¯ 2 n (Qi − Q) ¯ 2 (Ri − R) i=1

i=1

který pracuje s pořadími R1 , . . . , Rn , resp. Q1 , . . . , Qn , veličin X1 , . . . , Xn , resp. Y1 , . . . , Yn , díky čemuž je méně citlivý na odlehlá pozorování. Podrobnější informace týkající se především rozdělení těchto koeficientů za hypotézy nezávislosti lze nalézt v [1]. Hodnoty průtoku byly měřeny s přesností na dvě desetinná místa, srážkové úhrny s přesností na jedno desetinné místo. V důsledku toho se v datech vyskytuje mnoho shodných pozorování. Statistický software R při výpočtu Spearmanova korelačního koeficientu přiřazuje shodným pozorováním vždy průměr z odpovídajících pořadí, tedy shodným měřením je přiřazeno stejné pořadí. V takové situaci se v [1] doporučuje pracovat s tzv. korigovaným Spearmanovým korelačním koeficientem (více v [8]). Ten je definován vztahem P 6 ni=1 (Ri − Qi )2 rS,korig = 1 − 3 , n − n − Tx − Ty 7

kde

1X 3 1X 3 (tx − tx ), Ty = (ty − ty ). 2 2 Symbol tx označuje počty stejně velkých pozorování ve výběru X1 , . . . , Xn . Analogicky je definováno ty . Všechny výše uvedené korelační koeficienty jsme vypočetli pro zpracovávaná data, abychom zjistili, s jakou formou srážkových úhrnů průtok nejvíce koreluje. Nejprve jsme vypočítali korelaci průtoku a srážkového úhrnu z odpovídajícího dne, nebo některého předchozího dne, avšak největší hodnotu korelačních koeficientů jsme dostali pro průtok a součet srážkových úhrnů za několik předchozích dní. Pro každou stanici bylo přitom dosaženo nejvyšší korelace pro součet za jiný počet dní. Stejně tak byl velký rozdíl mezi zimním a letním obdobím. Výsledky jsou uvedeny v tabulkách 3.1, 3.2, 3.3, 3.4,3.5 a 3.6. Nejvyšší hodnota v řádku je vždy uvedena tučně. V souladu s naším očekáváním jsou hodnoty všech korelací kladné, tedy závislost velikosti průtoku na součtu srážkových úhrnů je pozitivní. Po porovnání hodnot každého koeficientu v letním a zimním období lze navíc prohlásit, že tato závislost je v zimě daleko slabší než v létě. Dále je patrné, že v letním období je pro většinu stanic nejvyšší Pearsonova korelace mezi průtokem a součtem srážkových úhrnů za poslední čtyři nebo pět dní. Výjimkou je pouze stanice Krnov. Nejvíce průtok koreloval se součty srážkových úhrnů ze stanic Heřmanovice, Rejvíz a Vidly. V zimním období hodnoty Pearsonova korelačního koeficientu u všech stanic neustále rostou. Největší hodnoty pak příslušejí stanicím Albrechtice-Žáry, Lichnov a Vidly. Hodnoty klasického i korigovaného Spearmanova korelačního koeficientu rostou v létě i v zimě pro všechny stanice. Co se týče klasického Spearmanova koeficientu, pak největší korelaci v letním období pozorujeme mezi průtokem a součty srážkových úhrnů ze stanic Rejvíz, Vidly a Albrechtice-Žáry. V zimě je závislost největší pro součty srážek naměřených ve stanicích Vidly, Lichnov a AlbrechticeŽáry. Pro korigovaný Spearmanův koeficient je největší korelace dosaženo pro součty srážkových úhrnů pocházejících ze stanic Praděd, Rejvíz a Vidly, a to v letním i zimním období. Příkaz pro výpočet výběrového i Spearmanova korelačního koeficientu je standardně obsažen v softwaru R. Funkce korigSpearman pro výpočet korigovaného Spearmanova korelačního koeficientu se nachází na přiloženém CD v souboru s názvem "korigSpearman.R". Tx =

8

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.6199 0.5675 0.4370 0.4560 0.4571 0.5568 0.6051 0.6020 0.5477

4 0.6329 0.5884 0.4810 0.4832 0.4788 0.5661 0.6192 0.6125 0.5740

5 0.6238 0.5866 0.5033 0.4896 0.4807 0.5565 0.6117 0.6034 0.5755

6 0.6083 0.5754 0.5136 0.4876 0.4752 0.5409 0.5971 0.5881 0.5675

7 0.5949 0.5640 0.5193 0.4851 0.4697 0.5281 0.5850 0.5753 0.5597

8 0.5835 0.5530 0.5217 0.4809 0.4637 0.5177 0.5738 0.5640 0.5520

9 0.5717 0.5418 0.5209 0.4734 0.4559 0.5065 0.5621 0.5521 0.5427

Tabulka 3.1: Pearsonovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v létě.

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.2692 0.2800 0.2857 0.2787 0.2734 0.2794 0.3065 0.2926 0.2999

4 0.2982 0.3095 0.3164 0.3064 0.2978 0.3041 0.3400 0.3237 0.3328

5 0.3243 0.3368 0.3420 0.3302 0.3185 0.3239 0.3685 0.3513 0.3609

6 0.3469 0.3600 0.3650 0.3515 0.3375 0.3409 0.3913 0.3750 0.3854

7 0.3661 0.3796 0.3841 0.3686 0.3528 0.3555 0.4110 0.3963 0.4071

8 0.3827 0.3964 0.3979 0.3818 0.3651 0.3672 0.4275 0.4148 0.4241

9 0.3963 0.4095 0.4080 0.3908 0.3748 0.3769 0.4406 0.4302 0.4379

Tabulka 3.2: Spearmanovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v létě.

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.4226 0.2874 0.3226 0.2786 0.3223 0.5468 0.4576 0.4583 0.3006

4 0.4451 0.3168 0.3523 0.3064 0.3448 0.5618 0.4833 0.4819 0.3336

5 0.4657 0.3440 0.3773 0.3302 0.3640 0.5742 0.5056 0.5031 0.3619

6 0.4836 0.3671 0.3998 0.3515 0.3816 0.5848 0.5236 0.5213 0.3866

7 0.4990 0.3866 0.4186 0.3686 0.3959 0.5940 0.5391 0.5378 0.4084

8 0.5123 0.4035 0.4323 0.3818 0.4074 0.6013 0.5521 0.5522 0.4256

9 0.5232 0.4166 0.4425 0.3908 0.4164 0.6075 0.5625 0.5640 0.4395

Tabulka 3.3: Korigované Spearmanovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v létě.

9

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.1796 0.1711 0.2017 0.2105 0.2173 0.1254 0.1898 0.2279 0.2105

4 0.1855 0.1773 0.2108 0.2210 0.2265 0.1279 0.1936 0.2395 0.2197

5 0.1863 0.1788 0.2140 0.2248 0.2290 0.1265 0.1920 0.2436 0.2239

6 0.1872 0.1781 0.2139 0.2243 0.2270 0.1238 0.1905 0.2437 0.2257

7 0.1898 0.1799 0.2150 0.2264 0.2255 0.1237 0.1923 0.2449 0.2299

8 0.1917 0.1827 0.2185 0.2301 0.2263 0.1262 0.1960 0.2470 0.2368

9 0.1942 0.1860 0.2214 0.2338 0.2262 0.1287 0.1981 0.2496 0.2416

Tabulka 3.4: Pearsonovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v zimě.

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.1489 0.1370 0.1396 0.1551 0.1496 0.1598 0.1528 0.1893 0.1638

4 0.1601 0.1508 0.1518 0.1715 0.1639 0.1671 0.1598 0.2078 0.1802

5 0.1707 0.1610 0.1624 0.1862 0.1728 0.1717 0.1654 0.2226 0.1928

6 0.1796 0.1683 0.1715 0.1996 0.1803 0.1769 0.1712 0.2345 0.2030

7 0.1880 0.1754 0.1794 0.2115 0.1863 0.1831 0.1789 0.2462 0.2136

8 0.1965 0.1825 0.1863 0.2226 0.1910 0.1895 0.1880 0.2568 0.2235

9 0.2051 0.1900 0.1939 0.2330 0.1943 0.1965 0.1971 0.2668 0.2337

Tabulka 3.5: Spearmanovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v zimě.

d Herm Karl Krno Lich Opav Prad Rejv Vidl Zary

3 0.3336 0.1369 0.1863 0.1551 0.2081 0.4619 0.3490 0.3737 0.1974

4 0.3415 0.1508 0.1973 0.1715 0.2208 0.4661 0.3535 0.3872 0.2130

5 0.3496 0.1610 0.2071 0.1862 0.2289 0.4689 0.3578 0.3983 0.2251

6 0.3565 0.1683 0.2156 0.1996 0.2358 0.4722 0.3624 0.4072 0.2349

7 0.3631 0.1754 0.2231 0.2115 0.2414 0.4762 0.3686 0.4159 0.2450

8 0.3697 0.1825 0.2296 0.2226 0.2457 0.4804 0.3758 0.4239 0.2545

9 0.3764 0.1900 0.2368 0.2330 0.2488 0.4849 0.3832 0.4314 0.2643

Tabulka 3.6: Korigované Spearmanovy korelační koeficienty pro průtok a součet srážek za d posledních dní v jednotlivých stanicích v zimě.

10

3.2

Koeficient závislosti chvostů

Náš zájem je soustředěn hlavně na vliv velkých srážkových úhrnů. Lze předpokládat, že větší množství srážek způsobuje vyšší průtok, a tedy ze statistického hlediska jde o problém závislosti velmi vysokých hodnot dvou náhodných veličin. Tuto závislost lze odhadnout pomocí tzv. koeficientu závislosti chvostů (anglicky "tail-dependence coefficient"), který měří závislost extrémně vysokých hodnot daných veličin X, Y . Podrobně je tento koeficient studován v [5]. Pokud by měly náhodné veličiny X a Y shodné marginální rozdělení, pak přirozenou mírou jejich extremální závislosti by byl například výraz χ = lim∗ P (Y > z|X > z), z→z

(3.1)

kde z ∗ je horní mez nosiče daného marginálního rozdělení. Tedy χ je pravděpodobnost, že jedna veličina nabývá velké hodnoty za podmínky, že druhá veličina nabývá velké hodnoty. Pro zobecnění výrazu (3.1) na nestejně rozdělené spojité náhodné veličiny X, Y stačí provést transformaci (U, V ) = {FX (X), FY (Y )}, kde FX a FY značí marginální distribuční funkce veličin X a Y . Veličiny U a V mají rovnoměrné rozdělení na intervalu [0, 1]. Podle (3.1) tedy dostáváme χ = lim P (V > u|U > u). u→1

(3.2)

Podle Sklarovy věty (viz [16]) pro dvojrozměrnou distribuční funkci F (X, Y ) s jednorozměrnými marginálami FX a FY existuje kopula C taková, že F (x, y) = C(FX (x), FY (y)) pro všechna reálná x a y. Funkce C je vlastně sdruženou distribuční funkcí veličin U a V . Jelikož u FX a FY předpokládáme spojitost, tak C je dána jednoznačně. Všimněme si dále, že platí P (U ≤ u nebo V ≤ u) = P (U ≤ u) + P (V ≤ u) − P (U ≤ u, V ≤ u) = 2u − C(u, u) a zároveň P (U ≤ u nebo V ≤ u) = 1 − P (U > u, V > u).

Odtud můžeme vyjádřit

P (U > u, V > u) = 1 − 2u + C(u, u). Výraz v (3.2) lze tedy ještě dále upravit P (U > u, V > u) 1 − 2u + C(u, u) = (3.3) P (U > u) 1−u log C(u, u) 1 − C(u, u) ∼2− , pro u → 1. (3.4) = 2− 1−u log u

P (V > u|U > u) =

Označíme-li

χ(u) = 2 −

log P (U ≤ u, V ≤ u) , pro 0 ≤ u ≤ 1, log P (U ≤ u)

(3.5)

pak χ = lim χ(u). u→1

11

(3.6)

Limita χ se nazývá koeficient závislosti chvostů veličin X a Y . Funkci χ(u) lze interpretovat jako míru závislosti v daném kvantilu u, přičemž znaménko hodnoty χ(u) udává, zda jsou veličiny v kvantilu u pozitivně či negativně závislé. Stejně jako klasický korelační koeficient je funkce χ(u) shora ohraničená hodnotou 1, neboť pro U a V zcela závislé je P (U ≤ u, V ≤ u) = P (U ≤ u). Pro spodní odhad χ(u) si stačí uvědomit, že C(u, u) = 2u − 1 + P (U > u, V > u) ≥ 2u − 1. Dolní mez χ(u) je tedy dána výrazem 2 − log(2u − 1)/ log(u) a interpretuje se jako −∞ pro u ≤ 1/2 a jako 0 pro u = 1. Pro nezávislé veličiny je χ(u) = 0 pro všechna u, při perfektní závislosti je χ(u) = 1 pro všechna u. Při výpočtu odhadu χ(u) ˆ se funkce FX a FY nahradí empirickými distribučními funkcemi FˆX a FˆY . Hodnoty (ui , vi ), i = 1, . . . , n, definované jako ui = FˆX (xi ) a vi = FˆY (yi ) jsou pak nezávislé realizace veličin U a V . Kopula C(u, u) se odhaduje pomocí n 1X ˆ C(u, u) = I[ui ≤ u & vi ≤ u], n i=1

kde I[ ] značí indikátor daného jevu. Odhad χˆ koeficientu závislosti chvostů se snažíme určit "od oka" z grafu funkce χ(u). ˆ Odhad této limity však nelze učinit na základě chování χ(u) ˆ těsně v blízkosti bodu 1, neboť zde jsou její funkční hodnoty ovlivněny značným úbytkem počtu pozorování. Číslo χˆ tedy odhadneme přibližně jako hodnotu, na níž se funkce χ(u) ˆ ustálila, než začalo v blízkosti pravého konce intervalu docházet ke kolísání. Ve statistickém softwaru R je funkce pro výpočet a vykreslení χ(u) ˆ obsažena v knihovně evd pod názvem chiplot. V případě studovaných dat jsme jako veličinu Y uvažovali průtok a jako veličinu X jsme zkoušeli brát součet srážkových úhrnů za předešlých d dní v jednotlivých stanicích. Poznamenejme, že distribuční funkce takto definované veličiny X nesplňuje předpoklad spojitosti, neboť nulový součet srážkových úhrnů má nenulovou pravděpodobnost výskytu, a tedy FX má v bodě x = 0 skok. To ale v tomto případě zjevně nevadí, protože koeficient závislosti chvostů měří limitní chování veličin U a V pro u → 1 (resp. před transformací X a Y pro jejich vysoké hodnoty). Pro formální ospravedlnění transformace U = FX (X) však lze veličinu X s libovolnou přesností v okolí nuly spojitě předefinovat, aniž by byl ovlivněn zbytek její distribuční funkce, a tedy i limitní chování pro vysoké hodnoty. Příslušné odhady funkce χ(u) jsou vykresleny v příloze na obrázcích 7.1 a 7.2. Z obrázku 7.1 je patrné, že v letním období se s rostoucím d celá funkce χ(u) ˆ posouvá směrem k vyšším funkčním hodnotám. Zdá se nám, že pro d = 3 a 4 je závislost slabší, než pro d ≥ 5, kde se jednotlivé χ(u) ˆ liší již relativně méně. Jelikož bylo třeba se subjektivně rozhodnout, zvolili jsme pro další použití součet srážkových úhrnů za posledních 5 dní. V zimě je situace poněkud odlišná. Z obrázku 7.2 vidíme, že funkce χ(u) ˆ se pro jednotlivá d příliš neliší a mají téměř konstantní průběh. Vzhledem k této podobnosti jednotlivých χ(u) ˆ jsme se rozhodli zvolit stejné d jako v letním období, tj. uvažovat součet srážkových úhrnů za předešlých 5 dní. Rádi bychom také posoudili, pro které srážkoměrné stanice je míra závislosti chvostů největší. Na obrázku 3.1 jsou za tímto účelem vykresleny odhady funkce 12

0.25

Herm Karl Krno Lich Opav

0.10 0.00

0.05

Chi(u)

0.15

0.20

Prad Rejv Vidl Zary

Herm Karl Krno

−0.05

Chi(u)

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

χ(u) příslušné součtu srážkových úhrnů za předchozích 5 dní pro jednotlivé stanice, a to pro letní i zimní období. Na obou grafech je funkce kladná, a tedy závislost vysokých hodnot průtoku a vysokých hodnot pětidenního součtu srážek je v létě i v zimě pozitivní. V letním období se zdá, že vysoký průtok je nejvíce ovlivněn pětidenními součty srážek ze stanic Heřmanovice, Karlovice a Albrechtice-Žáry. V zimním období je nejsilnější závislost chvostů průtoku a pětidenních součtů srážkových úhrnů ze stanic Opava, Krnov a Albrechtice-Žáry. Výrazná je dále křivka příslušející stanici Praděd, která leží nejníže. Příčinou by mohl být fakt, že meteorologická stanice Praděd, která se nacházela (před svým zrušením v roce 2000) na vrcholu stejnojmenné hory, byla se svojí nadmořskou výškou 1491 m n. m. nejvýše položenou stanicí patřící České republice. V zimním období se zde díky nízkým teplotám mohly držet srážky ve formě sněhu, který po většinu zimy ležel na zemském povrchu a nepřispíval nijak do toku řeky. Po porovnání obrázků 3.1a a 3.1b lze dále říci, že v zimě jsou hodnoty χ(u) ˆ daleko nižší než v létě, tedy závislost chvostů průtoku a pětidenního součtu srážek je v zimě slabší. Pro zpracovávaná data jsme vypočetli Pearsonův, Spearmanův i korigovaný Spearmanův korelační koeficient a graficky znázornili koeficient závislosti chvostů. Každá z těchto metod nám dala trochu jiné návrhy na formu srážkových úhrnů, která s průtokem nejvíce souvisí. Naším hlavním zájmem je prozkoumání závislosti vysokých hodnot srážek a příslušného průtoku, proto zvolíme pro další analýzu součet srážkových úhrnů, které nám doporučil koeficient závislosti chvostů, tj. součet za předchozích pět dní. Charakteristiky polohy pro data obsahující k danému dni součet srážkových úhrnů za předešlých pět dní se nacházejí v tabulkách 3.7 a 3.8. Kompletní data by obsahovala 9416 hodnot v letním období a 6589 hodnot v zimním období.

0.90

0.92

0.94

0.96

0.98

0.90

u

Lich Opav Prad 0.92

Rejv Vidl Zary 0.94

0.96

0.98

u

(a) Léto

(b) Zima

Obrázek 3.1: Koeficienty závislosti chvostů pro průtok a součet srážek za předchozích pět dní v jednotlivých stanicích.

13

min Herm 0 Karl 0 Krno 0 Lich 0 0 Opav Prad 0 Rejv 0 Vidl 0 Zary 0

průměr 75% kv. 15.84 20.7 12.05 16.9 10.81 15.3 10.80 15.7 10.62 15.3 18.24 25.0 17.89 23.8 17.86 24.8 12.72 17.7

90% kv. 39.40 30.90 29.20 28.20 28.30 42.40 44.30 43.45 32.80

98% kv. 80.13 56.34 48.81 49.80 51.70 78.00 89.70 76.07 57.70

max 427.50 320.40 106.20 173.40 203.30 455.00 511.00 501.20 290.90

počet dat 8708 9382 9245 9416 9202 8072 8679 8616 9411

podíl kladných dat 80.9 % 82.4 % 82.4 % 84.3 % 82.4 % 77.1 % 79.5 % 82.2 % 88.5 %

Tabulka 3.7: Počty pozorování a charakteristiky polohy pro součet srážkových úhrnů za předchozích pět dní v letním období.

min Herm 0 Karl 0 Krno 0 Lich 0 0 Opav Prad 0 Rejv 0 Vidl 0 Zary 0

průměr 75% kv. 7.85 11.2 7.16 10.6 4.77 6.3 5.10 6.8 4.41 5.6 11.93 17.8 8.80 12.6 12.55 18.2 6.89 10.0

90% kv. 19.60 19.40 12.60 13.80 12.10 28.40 22.42 31.50 17.70

98% kv. max 36.25 84.00 32.72 59.50 25.56 60.80 27.30 51.20 25.70 76.90 47.04 84.80 39.57 83.30 55.70 108.80 30.90 55.50

počet dat 6077 6589 6470 6589 6437 5682 6039 6050 6500

podíl kladných dat 79.3 % 80.9 % 79.1 % 81.9 % 79.3 % 77.6 % 76.7 % 82.2 % 87.4 %

Tabulka 3.8: Počty pozorování a charakteristiky polohy pro součet srážkových úhrnů za předchozích pět dní v zimním období.

14

4. Kvantilová regrese Nyní bude naší snahou odhad podmíněného rozdělení průtoku při daných srážkových úhrnech, především se zaměříme na vysoké kvantily tohoto rozdělení. K tomuto účelu použijeme kvantilovou regresi, která umožňuje nalezení regresního modelu pro libovolný kvantil odezvy. Vznik této metody sahá do roku 1978, kdy R. Koenker a G. Bassett v [14] definovali pojem regresního kvantilu, který byl rozšířením výběrového kvantilu na regresní model. Této metodě se věnuje zatím ne příliš rozsáhlá literatura, jejímž převážným autorem je R. Koenker.

4.1

Kvantil a jeho odhad

Nechť náhodná veličina X má distribuční funkci F (x) = P (X ≤ x). Pak pro libovolné τ ∈ (0, 1) se hodnota (4.1)

F −1 (τ ) = inf{x : F (x) ≥ τ }

nazývá τ -kvantilem veličiny X. Funkce F −1 proměnné τ definovaná tímto předpisem se nazývá kvantilová funkce. V [11] je ukázáno, že τ -kvantil je možné nalézt jako hodnotu minimalizující střední hodnotu ztrátové funkce ρτ (z) = |z|(τ I[z ≥ 0] + (1 − τ )I[z < 0]) = z(τ − I[z < 0]), z ∈ R, kde I[ ] značí indikátor daného jevu. Nahradíme-li F empirickou distribuční funkcí Fˆ (x) = n−1

n X i=1

I[Xi ≤ x],

kde n značí rozsah výběru, pak minimalizací výrazu Z

ρτ (x − ξ) dFˆ (x) = n−1

n X i=1

ρτ (Xi − ξ)

dostaneme odhad τ -kvantilu, tzv. výběrový τ -kvantil. Ten bývá spjat hlavně s pojmem uspořádaného náhodného výběru. Zde však vidíme, že výběrový τ -kvantil ξˆτ je řešením minimalizační úlohy min ξ∈R

n X i=1

(4.2)

ρτ (Xi − ξ),

to jest 

min  ξ∈R

X

i∈{i:Xi ≥ξ}

τ |Xi − ξ| +

X



(1 − τ )|Xi − ξ| .

i∈{i:Xi <ξ}

(4.3)

Podle [11] je možné tento princip dále zobecnit na odhadování kvantilů v regresním modelu. 15

4.2

Teoretické základy kvantilové regrese

Uvažujme lineární regresní model (4.4)

Yi = x′i β + Ui , i = 1, . . . , n,

kde Y1 , . . . , Yn jsou pozorování, β ∈ Rp je neznámý parametr a xi , i = 1, . . . , n, jsou sloupcové vektory mající délku p. Předpokládejme, že jde o model s absolutním členem, tedy xi1 = 1, i = 1, . . . , n. Veličiny U1 , . . . , Un představují náhodné chyby. Rozšířením pojmu výběrového τ -kvantilu na regresní model (4.4) je regresní ˆ ), τ ∈ (0, 1), který byl zaveden v [14] jako řešení minimalizační úlohy τ -kvantil β(τ minp

β∈R

se ztrátovou funkcí

n X i=1

(4.5)

ρτ (Yi − x′i β)

ρτ (z) = |z|(τ I[z ≥ 0] + (1 − τ )I[z < 0]), z ∈ R.

(4.6)

Tato úloha je zcela analogická úloze (4.2). Výraz (4.5) lze opět přepsat do tvaru   X X (4.7) τ |Yi − x′i β| + (1 − τ )|Yi − x′i β| . minp  β∈R

i∈{i:Yi ≥x′i β}

i∈{i:Yi <x′i β}

Zaveďme si dále nové proměnné {ri+ , ri− : i = 1, . . . , n} představující kladnou ˆ ) lze vypočítat jako a zápornou část vektoru residuí. Pak regresní kvantil β(τ + − složku β optimálního řešení (β, r , r ) úlohy min : τ

n X

ri+

i=1

p

za podmínek:

X j=1

+ (1 − τ )

n X

ri− ,

i=1

xij βj + ri+ − ri− = Yi , i = 1, . . . , n

(4.8)

βj ∈ R, j = 1, . . . , p, ri+ , ri− ≥ 0, i = 1, . . . , n 0 < τ < 1. Problém (4.8) je úlohou lineárního programování, kterou je možné řešit simpleˆ ), xovou metodou. Optimální řešení této úlohy odpovídá regresnímu kvantilu β(τ ˆ ) nulových. pro nějž je právě p z residuí Yi − x′i β(τ ′ Funkce QY (τ |x) = x β(τ ), která danému τ přiřazuje τ -kvantil podmíněného rozdělení odezvy Y při daném x, se nazývá podmíněná kvantilová funkce. Symbolem β(τ ) označujeme populační regresní kvantil, jehož odhadem je právě regresní ˆ ). Regresní kvantil definuje odhad podmíněné kvantilové funkce. kvantil β(τ Uvažujme, že v regresním modelu (4.4) jsou chyby Ui nezávislé a stejně rozdělené s distribuční funkcí F . V tomto jednoduchém případě má podmíněná kvantilová funkce tvar QY (τ |x) = x′ β + F −1 (τ ) (4.9) 16

a všechny populační regresní kvantily splňují β(τ ) = (β0 + F −1 (τ ), β1 , . . . , βp )′ . Nechť funkce F je spojitá a nechť jí odpovídající hustota f splňuje f (F −1 (τ )) > 0. Označme n 1X xi x′i = Q, lim n→∞ n i=1

kde Q je pozitivně definitní matice řádu p × p. Pak    √  τ (1 − τ ) d −1 ˆ n β(τ ) − β(τ ) −→ Np 0, . Q (f (F −1 (τ )))2

(4.10)

ˆ ) má tedy asymptoticky p-rozměrné normální rozdělení. Důkaz Regresní kvantil β(τ tohoto tvrzení pro regresní medián (tj. τ = 0.5) lze nalézt v [2]. Výraz (4.10) je speciálním případem poněkud obecnějšího tvrzení. Pro τi ∈ (0, 1), kde i = 1, . . . , m, označme ζ = (β(τ1 )′ , . . . , β(τm )′ )′ posloupnost populačních regresních kvantilů modelu (4.4) a nechť f (F −1 (τi )) > 0 pro i = 1, . . . , m. Dále označme Ω matici typu m × m s prvky ωij =

min(τi , τj ) − τi τj . f (F −1 (τi ))f (F −1 (τj ))

ˆ 1 )′ , . . . , β(τ ˆ m )′ )′ Pak sdružené asymptotické rozdělení regresních kvantilů ζˆ = (β(τ má tvar √ √ d −1 ˆ j ) − β(τj )))m n(ζˆ − ζ) = ( n(β(τ (4.11) j=1 −→ Nmp (0, Ω ⊗ Q ), kde Q je jako výše a ⊗ značí Kroneckerův součin. Důkaz tohoto tvrzení je poměrně zdlouhavý a čtenář ho může nalézt v [14]. Funkce s(t) = 1/f (F −1 (t)), jejíž hodnota v bodě t = τ vystupuje ve vzorcích pro varianční matici, se někdy nazývá kvantilová hustota. Zderivujeme-li totiž identitu F (F −1 (t)) = t, zjistíme, že 1 d −1 F (t) = , −1 dt f (F (t)) tedy s(t) je derivací kvantilové funkce. Pro s(t) se také v angličtině používá název "sparsity function" (sparsity = řídkost), neboť tato funkce odráží "řídkost" pozorování Yi v okolí požadovaného kvantilu. Je zcela přirozené, že přesnost odhadu regresního kvantilu závisí na hodnotě s(τ ). Vyskytují-li se totiž data v okolí τ kvantilu jen řídce, je těžké tento kvantil odhadnout. Na druhé straně, je-li hodnota s(τ ) nízká, tedy v okolí τ -kvantilu se vyskytuje mnoho pozorování, bude možné kvantil odhadnout mnohem přesněji. O tom, jak se varianční matice rozdělení (4.10) odhaduje v praxi, budeme hovořit v sekci 4.5.1. Viděli jsme, že pro model (4.4) se stejně rozdělenými chybami má limitní ˆ ) vcelku jednoduchý tvar. V případě nestejně rozdělených chyb je rozdělení β(τ asymptotické rozdělení poněkud složitější. Označme H(τ ) = lim n−1 n→∞

n X i=1

17

xi x′i fi (Fi−1 (τ )).

(4.12)

 √ ˆ ) − β(τ ) má tvar tzv. "Huberova sendPak limitní varianční matice pro n β(τ viče": 
√  d ˆ ) − β(τ ) −→ (4.13) Np 0, τ (1 − τ )H(τ )−1 QH(τ )−1 . n β(τ

Výraz fi (Fi−1 (τ )) značí hustotu odezvy Yi vyčíslenou v podmíněném τ -kvantilu. V případě stejně rozdělených chyb jsou všechny hustoty fi i funkce Fi−1 stejné a varianční matice v (4.13) se zjednoduší na tvar (4.10). ˆ 1 )′ , . . . , β(τ ˆ m )′ )′ regresních kvantilů v modelu Máme-li opět vektor ζˆ = (β(τ s nestejně rozdělenými chybami, pak podle [11] platí √ √ d ˆ j ) − β(τj )))m (4.14) n(ζˆ − ζ) = ( n(β(τ j=1 −→ Nmp (0, V ).

Asymptotická varianční matice V je matice s bloky tvaru √ √ ˆ i ) − β(τi )), n(β(τ ˆ j ) − β(τj ))) = acov( n(β(τ = [min(τi , τj ) − τi τj ]H(τi )−1 QH(τj )−1 , i, j = 1, . . . , m, (4.15) kde "acov" značí asymptotickou kovarianci.

4.3

Hledání modelu pro zpracovávaná data

V závěru třetí kapitoly jsme se rozhodli zabývat závislostí průtoku a součtu srážkových úhrnů za předchozích pět dní. Je nesmyslné předpokládat, že by mohl existovat stejný model v celém rozsahu srážkových úhrnů. Naším hlavním zájmem je prozkoumat průtok způsobený vysokými srážkovými úhrny, proto se zaměříme na hledání modelu jen pro 10% nejvyšších pětidenních součtů srážkových úhrnů z každého období. Abychom se vyhnuli problémům s odlehlými pozorováními (míněna pozorování s velkou hodnotou regresoru), která by ovlivňovala vysoké kvantily, budeme jakožto nezávisle proměnnou uvažovat pouze ty hodnoty součtu srážkových úhrnů za předchozích pět dní, které se nacházejí mezi 90% a 99.5% výběrovým kvantilem počítaným ze sduženého výběru obsahujícího srážkové souhrny za pět dní ze všech stanic. Z letního období tedy vybereme hodnoty součtu srážkových úhrnů za 5 předešlých dní, které se nacházejí mezi 35.0 a 105.3 mm, ze zimního období vezmeme hodnoty mezi 20.2 a 53.5 mm. Náhodnou veličinu odpovídající těmto vysokým pětidenním součtům srážkových úhrnů budeme značit X, jim příslušný průtok budeme označovat Y . Počet pozorování budeme značit n. Horní index L nebo Z u těchto veličin bude udávat, zda se veličina vztahuje k letnímu nebo zimnímu období, neboť model budeme samozřejmě hledat zvlášť pro zimní a letní období. Srážkoměrné stanice budeme v modelech značit pomocí dolního indexu k = 1, . . . , 9. Z pětidenních součtů srážkových úhrnů pro každou stanici vybíráme jen hodnoty nacházející se v určitém intervalu. To znamená, že každé stanici přísluší jiný soubor hodnot X. Když pak k těmto realizacím veličiny X přiřadíme odpovídající průtoky, dostaneme ke každé stanici i jiný soubor průtokových dat Y . V tabulkách 4.1 a 4.2 jsou proto pro každou stanici uvedeny charakteristiky polohy jak pro veličinu X, tak pro veličinu Y . Rozložení dat si lze prohlédnout na obrázku 4.1, kde je znázorněn 2D-histogram pro veličiny X a Y vztahující se k obci Heřmanovice. Z obrázku je patrné, že 18

Herm X Y Karl X Y Krno X Y Lich X Y Opav X Y Prad X Y Rejv X Y Vidl X Y Zary X Y

min 35.00 1.14 35.00 1.26 35.00 1.26 35.00 1.04 35.00 1.26 35.00 0.72 35.10 1.26 35.00 1.03 35.00 1.04

25% kv. 40.80 5.40 39.02 6.12 38.20 6.75 38.60 5.81 39.20 5.78 40.30 4.85 41.00 5.05 40.90 5.07 39.10 5.95

medián 47.50 9.48 45.90 13.10 42.40 12.85 44.45 12.50 45.00 10.90 47.60 8.40 49.00 8.87 48.70 8.90 44.50 11.90

průměr 75% kv. 53.67 62.30 14.25 17.90 50.24 56.00 20.79 26.98 46.83 51.10 21.53 27.55 48.39 53.98 21.58 27.00 48.27 54.50 21.26 24.60 51.11 57.00 13.28 16.10 53.88 61.60 13.36 16.98 52.52 59.40 13.64 17.00 49.57 55.20 19.95 25.90

90% kv. 78.44 31.27 71.90 45.66 63.90 50.20 66.40 52.17 64.10 54.70 70.82 30.20 79.99 29.78 73.22 30.20 69.25 45.01

max nL 105.2 972 194.0 972 103.8 690 194.0 690 104.9 596 171.0 596 102.3 562 379.0 562 104.9 581 194.0 581 104.7 1080 194.0 1080 103.2 1202 194.0 1202 105.0 1185 194.0 1185 103.9 810 194.0 810

Tabulka 4.1: Charakteristiky polohy veličiny X a příslušného průtoku Y pro letní období.

Herm X Y Karl X Y Krno X Y Lich X Y Opav X Y Prad X Y Rejv X Y Vidl X Y Zary X Y

min 20.20 1.35 20.20 1.31 20.20 1.30 20.20 1.30 20.20 1.30 20.20 1.26 20.20 1.30 20.20 1.31 20.20 1.30

25% kv. 22.83 3.23 22.60 3.24 22.30 3.39 22.60 3.41 23.00 3.45 22.90 2.86 23.50 3.28 23.60 3.27 22.50 3.39

medián 27.40 5.25 26.20 5.41 25.00 5.53 26.80 6.26 25.80 6.55 27.40 4.80 27.60 5.22 28.90 5.34 25.80 5.69

průměr 75% kv. 29.22 32.88 8.77 10.67 27.95 31.50 8.39 10.20 27.06 30.20 10.50 15.00 28.39 32.90 10.14 13.90 28.45 31.60 11.31 16.80 29.51 34.20 6.89 8.14 29.87 34.12 8.67 10.40 31.09 36.60 7.53 9.60 28.12 31.60 9.51 11.80

90% kv. 40.57 20.24 38.04 18.32 36.74 25.80 37.06 23.22 38.94 29.06 41.74 14.00 42.69 20.70 44.69 16.69 39.00 23.00

max 52.9 60.6 51.8 60.6 51.5 48.9 51.2 48.9 49.7 48.9 53.4 60.6 53.3 60.6 53.5 44.8 53.4 60.6

nZ 534 534 599 599 277 277 275 275 253 253 1099 1099 692 692 1172 1172 481 481

Tabulka 4.2: Charakteristiky polohy veličiny X a příslušného průtoku Y pro zimní období.

19

směrem k nižším hodnotám obou těchto veličin se hustota dat zvyšuje. Stejné pozorování by bylo možno učinit i pro ostatní srážkoměrné stanice, jejichž 2Dhistogramy by byly velmi podobné obrázku 4.1. Na obrázku 4.2 jsou pak vykresleny neparametrické odhady tří kvantilů průtoku Y v závislosti na X pro stanici Heřmanovice. Každý bod na tomto grafu představuje výběrový kvantil vypočtený z hodnot y odpovídajících hodnotám veličiny X z intervalu (x − 25, x + 25). Poznamenejme, že obrázky pro ostatní stanice by vypadaly podobně. Vidíme, že jednotlivé kvantily se hlavně v letním období při vyšších hodnotách X značně rozbíhají. To ilustruje přítomnost heteroskedasticity v závislosti veličiny Y na veličině X. Ve snaze modelovat tuto heteroskedasticitu jsme provedli logaritmickou transformaci závisle i nezávisle proměnné. Tato transformace se nakonec ukázala jako rozumná, a tudíž nadále budeme hledat model popisující závislost log(Y ) na log(X). Jako výchozí jsme zvolili následující tvar podmíněné kvantilové funkce Qlog YkL (τ |xLk ) = αkL (τ ) + βkL (τ ) log xLk

Qlog YkZ (τ |xZk ) = αkZ (τ ) + βkZ (τ ) log xZk , k = 1, . . . , 9.

(4.16) (4.17)

První model odpovídá letnímu období, druhý model zimnímu období. Koeficienty αkL , βkL , αkZ a βkZ jsou funkcí τ , ale nemohou být zcela libovolné. Aby totiž funkce Q(τ ) byla kvantilovou funkcí ke spojité distribuční funkci, musí být spojitá a rostoucí. Tomuto požadavku musejí obě dvojice koeficientů vyhovovat. V následujícím textu budeme regresním τ -kvantilem nazývat nejen vektor ˆ )), ale i funkci α ˆ ) log x proměnné x. Analogicky pro populační (ˆ α(τ ), β(τ ˆ (τ ) + β(τ ˆ )) budeme nazývat regresní τ -kvantil. Složky vektoru (α(τ ), β(τ )) nebo (ˆ α(τ ), β(τ po řadě absolutní člen a směrnice. Pro dané τ tedy předpis αkL (τ ) + βkL (τ ) log xLk , resp. αkZ (τ ) + βkZ (τ ) log xZk , k = 1, . . . , 9.

(4.18) (4.19)

představuje populační regresní τ -kvantil pro letní, resp. zimní období. Je patrné, že uvažované regresní kvantily mají pro různé srážkoměrné stanice různé koeficienty. Tyto koeficienty ale závisejí i na zvoleném τ , čímž povolujeme, aby v rámci každé stanice měly různé kvantily odlišné absolutní členy i směrnice. Tím vystihneme případnou další heteroskedasticitu závislosti log(Y ) na log(X). Odhady koeficientů v (4.18) a (4.19) pro τ = 0.5 (tj. složky regresního mediánu) a příslušné 95% intervaly spolehlivosti jsou uvedeny v tabulkách 4.3 a 4.4. Zajímavé jsou zejména údaje týkající se směrnice regresního mediánu. V letním období jsou všechna βˆL kladná a žádný z příslušných intervalů spolehlivosti neobsahuje nulu, tedy závislost log(Y L ) na log(X L ) je pozitivní. V zimním období jsou bodové odhady směrnic sice také kladné, ovšem téměř všechny intervaly spolehlivosti obsahují nulu. V tomto období tudíž nelze zamítnout hypotézu, že regresní medián logaritmu průtoku Y je konstantní funkcí a nezávisí na velikosti logaritmu X. Regresní medián spolu s 20%, 35%, 65% a 80% regresními kvantily pro jednotlivé stanice je vykreslen na obrázku 4.3. Poznamenejme, že intervaly spolehlivosti uvedené v tabulkách 4.3 a 4.4 byly pro směrnici i absolutní člen vypočteny stejnou metodou, a to inverzí pořadového 20

150 100 50

Y

40

50

60

70

80

90

100

X

30 10

20

Y

40

50

60

(a) Léto

25

30

35

40

45

50

X

(b) Zima

Obrázek 4.1: 2D-histogram veličin X a Y pro stanici Heřmanovice pro letní i zimní období. Zelená a žlutá barva odpovídají větší koncentraci dat. Podél obou os je znázorněno marginální rozdělení dat.

21

20 15

60

80% kvantil 50% kvantil 20% kvantil

10

Y

30

0

0

10

5

20

Y

40

50

80% kvantil 50% kvantil 20% kvantil

40

50

60

70

80

90

100

20

X

30

40

50

X

(a) Léto

(b) Zima

Obrázek 4.2: Neparametrické odhady 20%, 50% a 80% kvantilu veličiny Y v závislosti na X pro letní i zimní období. Data ze stanice Heřmanovice. α ˆL Herm -2.0489 Karl -3.7639 Krno -5.5364 Lich -4.1186 Opav -4.0470 Prad -1.3120 Rejv -1.1919 Vidl -1.3562 Zary -5.3776

95% int. spol. (-3.0085;-0.9836) (-5.5146;-2.3589) (-6.6308;-3.5236) (-5.9709;-2.9721) (-6.4820;-0.7124) (-2.3221;-0.2507) (-2.4946;-0.2246) (-2.2027;-0.5346) (-6.4803;-4.1828)

βˆL 1.1012 1.6414 2.1404 1.7291 1.6984 0.8769 0.8558 0.9024 2.0557

95% int. spol. (0.8134;1.3770) (1.2783;2.0882) (1.5766;2.4267) (1.4349;2.1972) (0.7966;2.3314) (0.6153;1.1492) (0.6060;1.1711) (0.7022;1.1135) (1.7405;2.3371)

Tabulka 4.3: Složky regresního mediánu a příslušné 95% intervaly spolehlivosti v letním období. α ˆZ Herm 0.4906 Karl 1.4238 Krno 1.1005 Lich 0.8172 Opav 1.5918 Prad 1.9015 Rejv 1.0630 Vidl 0.9094 Zary 1.3654

95% int. spol. (-0.9320;2.2602) (-0.9193;3.0964) (-0.2612;3.0807) (-0.3910;3.2921) (-0.3185;2.6459) (0.8422;2.4661) (-0.2721;2.2512) (0.1528;1.5145) (0.4915;2.9642)

βˆZ 0.3480 0.0832 0.1873 0.2984 0.0952 -0.0976 0.1773 0.2239 0.1131

95% int. spol. (-0.1804;0.9097) (-0.4072;0.8056) (-0.3551;0.5603) (-0.4821;0.6890) (-0.2400;0.5984) (-0.2472;0.1903) (-0.1615;0.5971) (0.0513;0.4490) (-0.3465;0.3745)

Tabulka 4.4: Složky regresního mediánu a příslušné 95% intervaly spolehlivosti v zimním období.

22

testu, který je detailně popsán v [10] a [11]. Tento typ intervalu spolehlivosti je kladně hodnocen v [10] a je také přednastaven v softwaru R.

4.3.1

Výpočet v R

Ve statistickém softwaru R je pro kvantilovou regresi určena knihovna quantreg, která obsahuje metody a testy popsané v [11]. Jejím autorem je prof. Roger Koenker. Složky regresního τ -kvantilu lze získat pomocí příkazu: rq(formula,tau= ,method= ,...) Parametr formula představuje požadovaný model, tau určuje regresní kvantil, který chceme vypočítat. Pomocí parametru method lze zvolit algoritmus, který bude pro výpočet použit. My jsme použili Barrodale-Robertsův simplexový algoritmus (method=br), který je nastaven jako výchozí a dle [12] se hodí pro rozsahy výběru menší než několik tisíc. Ostatní metody a další možné parametry jsou popsány v manuálu [12]. Označíme-li výsledek předchozího příkazu jako fit, pak 95% konfidenční intervaly pro jednotlivé složky regresního τ -kvantilu v případě nestejně rozdělených chyb dostaneme příkazem: summary(fit,se="rank",iid=FALSE,alpha=0.05) Tyto konfidenční intervaly jsou vypočteny inverzí pořadového testu, jak již bylo zmíněno dříve. Směrodatné odchylky, p-hodnoty a hodnoty testových statistik lze získat pomocí příkazu summary(fit,se=" ") Předpokládáme-li, že chyby v modelu jsou stejně rozdělené, volíme se="iid". Pokud tento předpoklad nemáme, použijeme se="nid". Detailně jsou všechny možnosti odhadu směrodatné odchylky popsány v [11]. Pro případ se="iid" se budeme praktickými stránkami výpočtu zabývat v sekci 4.5.1.

4.4

Test rovnoběžnosti regresních kvantilů

Model (4.16) a (4.17) je velmi obecný, a tak bychom jej rádi nějakým způsobem zjednodušili. Při pohledu na obrázek 4.3 by se nám vhodným zjednodušením zdál model, v němž by byly všechny regresní kvantily rovnoběžné. To by znamenalo, že chyby v tomto modelu by bylo možné považovat za stejně rozdělené. Jelikož je ale regresních kvantilů (nebo přesněji různých τ ) nekonečně mnoho, nelze testovat obecnou hypotézu o rovnoběžnosti všech regresních kvantilů, nýbrž je nutné se omezit jen na určitý konečný počet. My jsme pro test rovnoběžnosti zvolili tyto kvantily: τ = 0.2, 0.35, 0.5, 0.65 a 0.8. Podle [11] lze test obecné hypotézy H0 : Rζ = r

23

o vektoru ζ = (β(τ1 )′ , . . . , β(τm )′ )′ populačních regresních kvantilů založit na statistice T = n(Rζˆ − r)′ [RV R′ ]−1 (Rζˆ − r), (4.20)

kde V je matice typu mp × mp s ij-tým blokem tvaru

V (τi , τj ) = [min(τi , τj ) − τi τj ]H(τi )−1 QH(τj )−1 , jak již bylo definováno v (4.15). Matice R je typu s × mp. Z (4.14) víme, že √

d

n(ζˆ − ζ) −→ Nmp (0, V ).

Za nulové hypotézy tedy platí √

d

n(Rζˆ − r) −→ Ns (0, RV R′ ).

Z věty o rozdělení kvadratických forem (viz např. [1]) pak dostáváme, že statistika T má za H0 asymptoticky rozdělení χ2q , kde q je hodnost matice R. Náš záměr je testovat rovnost směrnic regresních kvantilů pro τ = 0.2, 0.35, 0.5, 0.65, 0.8, a tedy m = 5. Navíc regresní kvantily v našem případě jsou přímky, tudíž p = 2. Matice R odpovídající testu rovnoběžnosti je v našem případě tvaru   0 −1 0 1 0 0 0 0 0 0 0 0 0 −1 0 1 0 0 0 0   0 0 0 0 0 −1 0 1 0 0 0 0 0 0 0 0 0 −1 0 1 a vektor r = (0, 0, 0, 0)′ . Vektor

ζ = (α(0.2), β(0.2), α(0.35), β(0.35), α(0.5), β(0.5), α(0.65), β(0.65), α(0.8), β(0.8))′ , přičemž α a β značí absolutní člen a směrnici pro danou stanici v letním nebo zimním období. Hodnost matice R je q = (m − 1)(p − 1) = 4. V softwaru R se test rovnosti směrnic regresních kvantilů provádí pomocí funkce anova.rq. Tento příkaz ovšem nevrací přímo statistiku T definovanou v (4.20), ale dává ji podělenou příslušnými stupni volnosti (m − 1)(p − 1). Jde tedy o výraz připomínající F -statistiku. Udávaná p-hodnota je pak vypočtena jako 1 − G(T /[(m − 1)(p − 1)]), kde G je distribuční funkce F -rozdělení se stupni volnosti (m − 1)(p − 1) a nm − (p − 1)(m − 1). Pro takovéto zacházení s χ2 statistikou neexistuje žádné formální ospravedlnění. Podle [13] je to ale běžná praxe v literatuře zabývající se robustní analýzou, neboť se ukazuje, že tento přístup dává lepší výsledky. Hodnoty statistiky T /[(m − 1)(p − 1)] a příslušné p-hodnoty pro analyzovaná data jsou uvedeny v tabulce 4.5. Zvýrazněny jsou p-hodnoty nižší než 5 %. Z tabulky vidíme, že rovnoběžnost zvolených regresních kvantilů se zamítá jen pro modely odpovídající stanicím Heřmanovice, Krnov a Opava v letním období. V zimním období je nulová hypotéza zamítnuta pro modely náležící stanicím Karlovice a Albrechtice-Žáry. Jak jsme již nastínili v úvodu této podkapitoly, nezamítnutí rovnosti směrnic m regresních kvantilů nám nezaručuje platnost podmodelu, v němž by byly všechny regresní kvantily rovnoběžné. My jsme se však 24

Léto statistika T /q p-hodnota Herm 4.0381 0.0029 Karl 2.2772 0.0587 Krno 3.6949 0.0053 Lich 0.7307 0.5709 7.7824 3.0984 · 10−6 Opav Prad 0.9239 0.4488 Rejv 1.4175 0.2253 Vidl 2.1536 0.0716 Zary 2.1112 0.0768

Zima statistika T /q p-hodnota 1.4397 0.2182 4.2675 0.0019 0.6640 0.6171 1.1113 0.3496 0.4096 0.8018 2.0203 0.0888 1.7135 0.1441 1.3878 0.2354 2.4654 0.0431

Tabulka 4.5: Hodnoty statistiky T /q a příslušné p-hodnoty pro test rovnoběžnosti regresních kvantilů. snažili pro testování zvolit takové regresní kvantily, které by byly coby zástupci všech regresních kvantilů dostatečně "reprezentativní". Až na výjimky nebyla rovnoběžnost těchto kvantilů zamítnuta, proto budeme za účelem zjednodušení pro všechny srážkoměrné stanice předpokládat platnost modelu Qlog YkL (τ |xLk ) = αkL + βkL log xLk + [FkL ]−1 (τ )

Qlog YkZ (τ |xZk ) = αkZ + βkZ log xZk + [FkZ ]−1 (τ ),

(4.21) (4.22)

kde opět k označuje srážkoměrnou stanici, FkL značí distribuční funkci chyb v letním a FkZ distribuční funkci chyb v zimním období. Rovnoběžnost regresních kvantilů pro zlogaritmované veličiny znamená, že distribuční funkce Flog Y veličiny log(Y ) má bez ohledu na hodnotu log(X) stále stejný tvar F0 , jen se s rostoucí nezávisle proměnnou posouvá (pro β > 0) směrem vzhůru, tj. platí Flog Y (log y| log X = log x) = F0 (−α − β log x + log y). Pro původní veličiny X a Y to znamená, že  y    y  = F , FY (y|X = x) = FY (y| log X = log x) = F0 log 1 Axβ Axβ

kde α = log A a F1 = F0 ◦ log. Z důvodu názornosti jsme vynechali označení nebo Z i index k. Zvolme si hodnoty x1 a x2 takové, že x2 = cx1 . Pak   y 1 FY (y|X = x1 ) = F1 = F1 (by), kde b = β Ax1 Axβ1       y y 1 = F1 FY (y|X = x2 ) = F1 = F1 β by A(cx1 )β c Axβ2

L

Odtud plyne, že FY (y|X = x2 ) = FY (y/cβ |X = x1 ), kde cβ = (x2 /x1 )β . Podmíněná distribuční funkce veličiny Y se tedy s rostoucími srážkovými úhrny X mění ve smyslu změny měřítka, přičemž pro dané x1 a x2 tato změna odpovídá (x2 /x1 )β . 25

V tabulkách 4.6 a 4.7 uvádíme pro zajímavost koeficienty zvolených regresních kvantilů (τ =0.2, 0.35, 0.5, 0.65, 0.8) vypočtené za předpokladu shodné směrnice. Tyto hodnoty jsou výsledkem minimalizace min

(α,β)

5 X n X j=1 i=1

ρτj (Yi − αj − βxi ),

(4.23)

kde α = (α1 , . . . , α5 )′ značí vektor absolutních členů a β je společná směrnice. Funkce ρτj byla definována v (4.6). V článku [9] je dokázáno, že v modelu s nezávislými a stejně rozdělenými chybami má vypočtený odhad βˆ asymptoticky normální rozdělení se střední hodnotou odpovídající skutečné hodnotě směrnice. V příloze na obrázku 7.3 si lze prohlédnout grafické srovnání regresních kvantilů vypočtených se shodnou směrnicí s kvantily ze sekce 4.3. ˆ byly vypočteny v softwaru R pomocí funkce rq.fit.hogg, Odhady (α, ˆ β) kterou laskavě poskytl prof. Koenker v [13]. Funkce Hogg.se pro výpočet směrodatných odchylek se nachází na přiloženém CD.

τ1 = 0.2 Herm -2.6384 -4.4685 Karl Krno -5.8047 Lich -4.7369 Opav -4.3127 Prad -2.0427 -1.8162 Rejv Vidl -1.8589 Zary -5.7012

absolutní členy τ2 = 0.35 τ3 = 0.5 τ4 = 0.65 τ5 = 0.8 -2.3179 -2.0016 -1.6716 -1.2908 -4.0804 -3.6464 -3.2315 -2.8338 -5.4120 -5.0401 -4.6696 -4.2491 -4.3257 -3.9656 -3.5227 -3.0246 -3.8976 -3.5098 -3.0536 -2.5852 -1.7862 -1.4682 -1.1429 -0.6302 -1.5281 -1.1890 -0.8456 -0.4034 -1.5366 -1.2331 -0.8918 -0.4314 -5.3165 -4.9206 -4.6009 -4.1808

směrnice ˆ β se 1.0886 0.0847 1.6116 0.1426 2.0061 0.1668 1.6882 0.1846 1.5523 0.1958 0.9179 0.1077 0.8550 0.0860 0.8696 0.0888 1.9381 0.1255

Tabulka 4.6: Odhady absolutních členů a odhad společné směrnice se směrodatnou odchylkou (se) pro data z letního období. absolutní členy směrnice ˆ β se τ1 = 0.2 τ2 = 0.35 τ3 = 0.5 τ4 = 0.65 τ5 = 0.8 Herm -0.0024 0.2138 0.5234 0.8849 1.3872 0.3392 0.1451 Karl 0.2438 0.6487 0.9382 1.2652 1.7383 0.2317 0.1598 0.7875 1.1276 1.4260 1.9529 2.6118 0.0848 0.2669 Krno Lich 0.6602 0.9638 1.3298 1.7851 2.3607 0.1446 0.2587 Opav 1.2432 1.5680 2.0082 2.3821 3.0983 -0.0387 0.2737 Prad 1.0059 1.3855 1.6502 1.9296 2.3396 -0.0242 0.0984 Rejv 0.1194 0.3871 0.6643 1.0263 1.5187 0.2955 0.1386 0.2773 0.6292 0.8921 1.1686 1.6253 0.2290 0.0869 Vidl Zary 0.3235 0.6209 0.9170 1.3131 1.8732 0.2457 0.1607 Tabulka 4.7: Odhady absolutních členů a odhad společné směrnice se směrodatnou odchylkou (se) pro data ze zimního období.

26

4.5

Porovnání letního a zimního období

Na základě testu provedeného v předchozí podkapitole jsme se rozhodli v rámci každé stanice pokládat regresní kvantily za rovnoběžné. Nyní bude naším cílem pro každou stanici porovnat sklon regresních kvantilů v letním a zimním období. Pro každé z těchto období je v rámci dané stanice uvažován stejný model, avšak hodnoty odezvy i vysvětlující proměnné se liší. Nelze tedy použít klasické testy pro porovnání směrnic téhož regresního kvantilu ve dvou modelech, neboť tyto testy vyžadují, aby jeden z modelů byl podmodelem druhého. Místo toho využijeme asymptotické normální rozdělení regresních kvantilů. Uvažujme, že závislost velkých součtů srážkových úhrnů a odpovídajícího průtoku lze pro všechny stanice popsat modelem (4.21) a (4.22). Pak pro každou stanici k = 1, . . . , 9 podle (4.10) platí !   1 τ (1 − τ ) −1 βˆkL (τ ) − βkL (τ ) ≈ N2 0, L QLk nk (fkL ([FkL ]−1 (τ )))2 !   τ (1 − τ ) 1 −1 βˆkZ (τ ) − βkZ (τ ) ≈ N2 0, Z QZk , nk (fkZ ([FkZ ]−1 (τ )))2 kde βkL (τ ) = (αkL + [FkL ]−1 (τ ), βkL )′ a βkZ (τ ) = (αkZ + [FkZ ]−1 (τ ), βkZ )′ . Varianční matice těchto rozdělení jsou typu 2 × 2 a označme je po řadě VkL (τ ) a VkZ (τ ). Marginální rozdělení složek normálně rozděleného vektoru je opět normální. Proto dostáváme
βˆkL ≈ N βkL , VkL (τ )22 . . . letní období,
βˆZ ≈ N β Z , V Z (τ )22 . . . zimní období. k

k

k

Považujeme-li βˆkL a βˆkZ za nezávislé pro dané k, pak podle věty o součtu nezávislých veličin s normálním rozdělením platí
βˆkL − βˆkZ ≈ N βkL − βkZ , VkL (τ )22 + VkZ (τ )22 . Test hypotézy H0 : βkL = βkZ lze tedy založit na intervalu spolehlivosti



q  α∗ −1 L Z ˆ ˆ VkL (τ )22 + VkZ (τ )22 , 1− βk − βk − Φ 2  q  ∗ α L Z −1 Z L 1− βˆk − βˆk + Φ Vk (τ )22 + Vk (τ )22 . (4.24) 2

Symbol Φ značí distribuční funkci normovaného normálního rozdělení. Hladinu testu α∗ zvolíme jako 0.05. Na základě testu rovnoběžnosti pokládáme směrnice všech regresních kvantilů za shodné, avšak skutečnou hodnotu této společné směrnice není možné vypočítat. Odhady směrnic uvedené na konci předchozí podkapitoly pro jednotlivé stanice se týkají pouze pěti kvantilů. Proto zde společnou směrnici odhadneme prostřednictvím hodnoty vypočtené pro regresní medián a položíme βˆkL = βˆkL (0.5) a βˆkZ = βˆkZ (0.5). p p Odhady směrnic a jejich směrodatné odchylky VkL (τ )22 a VkZ (τ )22 jsou uvedeny v tabulce 4.8. V posledním sloupci této tabulky jsou vypočtené intervaly 27

Léto L ˆ β se Herm 1.1012 0.1197 Karl 1.6414 0.2422 Krno 2.1404 0.1739 Lich 1.7291 0.2216 Opav 1.6984 0.2413 Prad 0.8769 0.1303 Rejv 0.8558 0.1133 Vidl 0.9024 0.1093 Zary 2.0557 0.1557

Zima Z ˆ β se 0.3480 0.2018 0.0832 0.1925 0.1873 0.3413 0.2984 0.3687 0.0952 0.3356 -0.0976 0.0952 0.1773 0.1963 0.2239 0.1080 0.1131 0.2503

95% int. spol. pro βˆL − βˆZ (0.2932;1.2131) (0.9520;2.1645) (1.2024;2.7039) (0.5876;2.2738) (0.7931;2.4133) (0.6583;1.2907) (0.2342;1.1228) (0.3774;0.9796) (1.3648;2.5203)

Tabulka 4.8: Směrnice regresního mediánu, směrodatné odchylky a 95% intervaly spolehlivosti pro rozdíl směrnic regresního mediánu příslušného letnímu a zimnímu období. spolehlivosti pro (βˆkL − βˆkZ ). Vidíme, že žádný z intervalů neobsahuje nulu, tedy hypotézu o rovnosti směrnic βkL a βkZ zamítáme na hladině 5 %. Navíc z toho, že všechny intervaly jsou kladné, plyne, že regresní kvantily v letním období jsou daleko strmější než v zimním období. V zimním období tedy výše logaritmu srážek X ovlivňuje logaritmus průtoku Y mnohem méně.

4.5.1

Výpočet v R

Příkazy pro výpočet směrodatných odchylek v R byly zmíněny již na konci sekce 4.3.1. Nyní se zaměříme na praktické aspekty tohoto výpočtu. Uvažujme obecnou situaci, kdy máme model (4.9) a chceme odhadnout varianční matici P v (4.10). Matice Q−1 v této varianční matici se odhaduje přirozeně jako n( ni=1 xi x′i )−1 . Klíčový je však odhad hodnoty kvantilové hustoty s(τ ) = 1/f (F −1 (τ )) zmiňované v podkapitole 4.2. V knihovně quantreg je to provedeno následovně. Nejprve je odhadnuta kvantilová funkce F −1 pomocí residuí daného regresního modelu, tj.   i−1 i −1 ˆ , , F (t) = r(i) , t ∈ n n kde

ˆ ), i = 1, . . . , n. ri = Yi − x′i β(τ

Na základě Fˆ −1 se pak odhaduje kvantilová hustota s. Není však třeba odhadovat funkci s celou, neboť ve vzorci pro V (τ ) vystupuje pouze její hodnota v konkrétním bodě τ . V podkapitole 4.2 jsme ukázali, že funkce s je derivací kvantilové funkce F −1 . Pro dané τ ∈ (0, 1) tedy hodnotu s(τ ) odhadneme přirozeně jako směrnici tečny funkce Fˆ −1 v bodě τ . Kolem bodu (τ, Fˆ −1 (τ )) se vybere okolí h + 1 nejbližších bodů a těmi se proloží regresní medián. Jeho směrnice je pak odhadem s(τ ). Úskalím tohoto postupu je jen p nulových bodů (residuí), které se vyskytují kolem Fˆ −1 (τ ) a mohly by značně ovlivnit odhad směrnice. V R je toto ošetřeno jednoduše vynecháním těchto p nulových residuí z "výběru" r1 , . . . , rn . Odhad kvantilové funkce se tedy provádí na základě n − p hodnot. 28

Poslední otázkou je volba h. Software R bere h jako horní celou část z n · hn , kde hn značí Hallovu-Sheatherovu šířku. Ta je definována vztahem hn = n

−1/3

 −1 2/3 Φ (0.975)



1.5φ2 (Φ−1 (τ )) 2(Φ−1 (τ ))2 + 1

1/3

,

kde φ značí hustotu normovaného normálního rozdělení a Φ jeho distribuční funkci.

29

4

5

log(Y)

2

3

4 3

log(Y)

1

2 1 0 3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(b) Heřmanovice, zima

2

log(Y)

3 1

1

2

log(Y)

3

4

5

4

(a) Heřmanovice, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

log(X)

3.4

3.6

3.8

log(X)

(d) Karlovice, zima

3 2

log(Y)

3

1

2 1

log(Y)

4

5

4

(c) Karlovice, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6

3.8

log(X)

(e) Krnov, léto

(f) Krnov, zima

Obrázek 4.3: Regresní kvantily pro τ =0.2, 0.35, 0.5, 0.65 a 0.8. 30

4.0

4

6

3

5

2

log(Y)

4 3

log(Y)

1

2 1 0 3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(h) Lichnov, zima

2

log(Y)

3 1

1

2

log(Y)

4

3

5

4

(g) Lichnov, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(j) Opava, zima

2

log(Y) 0

1

1

2

log(Y)

3

3

4

5

4

(i) Opava, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6

3.8

log(X)

(k) Praděd, léto

(l) Praděd, zima

Obrázek 4.3: Regresní kvantily pro τ =0.2, 0.35, 0.5, 0.65 a 0.8. 31

4.0

4

5

2

log(Y)

3

4 3

log(Y)

1

2 1 3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

4.0

3.8

4.0

3.8

4.0

log(X)

(n) Rejvíz, zima

0

0.5

1

1.0

1.5

2.0

log(Y)

3 2

log(Y)

2.5

4

3.0

3.5

5

(m) Rejvíz, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6 log(X)

(p) Vidly, zima

2

log(Y)

3

1

2 1 0

log(Y)

3

4

5

4

(o) Vidly, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6 log(X)

(q) Albrechtice-Žáry, léto

(r) Albrechtice-Žáry, zima

Obrázek 4.3: Regresní kvantily pro τ =0.2, 0.35, 0.5, 0.65 a 0.8. 32

5. Metoda špiček nad prahem V předchozí kapitole jsme se snažili pomocí kvantilové regrese odhadnout kvantily logaritmu průtoku Y pro dané hodnoty log(X). Je však třeba zdůraznit, že velmi vysoké regresní kvantily v kvantilové regresi bývají zatíženy značnou nejistotou. Proto je vhodné zkusit k výpočtu vysokých kvantilů použít i některou z metod založených na modelování vysokých hodnot dané veličiny pomocí známého rozdělení s neznámými parametry. Velmi vysoké kvantily požadované veličiny se pak odhadnou prostřednictvím kvantilů tohoto rozdělení. Jedním z těchto postupů je metoda špiček nad prahem (anglicky "Peaks Over Threshold method"). Tato metoda pracuje s takovými hodnotami Yi z náhodného výběru Y1 , . . . , Yn , které přesahují určitou předem zvolenou hodnotu u zvanou "práh". Hodnoty Yi > u představují "špičky nad prahem". Nechť Y1 , . . . , Yn je náhodný výběr z rozdělení s distribuční funkcí F a označme Y libovolný člen této posloupnosti. Chování extrémně vysokých hodnot veličiny Y lze popsat pomocí podmíněné pravděpodobnosti P [Y > u + z|Y > u] =

1 − F (u + z) , z > 0. 1 − F (u)

Kdyby byla distribuční funkce F známa, bylo by známo i rozdělení přesahů prahu u. To se v praxi obvykle nestává, avšak podle teorie popsané v [4] lze pro dostatečně velké u distribuční funkci přesahů Z = Y − u za podmínky Y > u aproximovat funkcí  z −1/γ H(z) = 1 − 1 + γ (5.1) σ definované na množině {z : z > 0 & 1 + γ σz > 0}. Rozdělení s distribuční funkcí (5.1) se nazývá zobecněné Paretovo rozdělení s parametry γ ∈ R a σ > 0. Pro γ → 0 je  z H(z) = 1 − exp − , z > 0, (5.2) σ což odpovídá exponenciálnímu rozdělení. Odhady parametrů γ a σ založené na daném výběru Y1 , . . . , Yn a při zvoleném prahu u lze vypočítat například metodou maximální věrohodnosti. Hustota rozdělení (5.1) je tvaru h(z) =

1 z −1−1/γ 1+γ . σ σ

Nechť Z1 , . . . , Znu jsou přesahy prahu u pocházející z výběru Y1 , . . . , Yn . Pak věrohodnostní funkce pro případ γ 6= 0 je dána předpisem −1−1/γ nu  Zi 1 Y 1+γ . L(σ, γ) = h(Zi ) = nu σ i=1 σ i=1 nu Y

(5.3)

Odpovídající logaritmicko-věrohodnostní funkce je   nu   1 X Zi ℓ(σ, γ) = log L(σ, γ) = −nu log σ − 1 + log 1 + γ γ i=1 σ 33

(5.4)

za podmínky, že (1 + γ Zσi ) > 0 pro i = 1, . . . , nu . Maximálně věrohodné odhady σ ˆ a γˆ jsou prvky, v nichž věrohodnost (5.4) nabývá svého maxima. Provedeme-li reparametrizaci (σ, γ) → (ν, γ), kde ν =

γ , σ

doporučenou v [3], pak 

1 ℓ(ν, γ) = −nu log γ + nu log ν − 1 + γ

X nu

log(1 + νZi ).

i=1

∂ℓ ∂ℓ Položíme-li parciální derivace ∂ν a ∂γ rovny nule, zjistíme, že maximálně věrohodné odhady νˆ a γˆ parametrů ν a γ musejí vyhovovat rovnicím   nu 1 Zi 1 1 X − 1+ = 0, (5.5) νˆ γˆ nu i=1 1 + νˆZi nu 1 X log(1 + νˆZi ). γˆ = nu i=1

(5.6)

Analytický výpočet odhadů σ ˆ a γˆ bohužel není možný, a tak se při řešení využívají numerické metody. V případě γ = 0 se logaritmická věrohodnost odvodí na základě (5.2), přičemž dostaneme nu 1X Zi . ℓ(σ, 0) = −nu log σ − σ i=1

Maximálně věrohodný odhad σ pro exponenciální rozdělení P u lze vyjádřit explicitně a je roven výběrovému průměru přesahů, tj. σ ˆ = n1u ni=1 Zi . Z teorie maximální věrohodnosti je známo, že vypočtený maximálně věrohodný odhad má za platnosti určitých podmínek asymptoticky normální rozdělení. Podle [3] toto rozdělení sleduje i odhad (ˆ σ , γˆ )′ vektoru parametrů (σ0 , γ0 )′ , ale pouze v případě, že γ0 > −0.5. Při tomto omezení tedy platí     √ σ0 σ ˆ d − nu −→ N2 (0, IE (σ0 , γ0 )−1 ), (5.7) γ0 γˆ

kde IE (σ, γ) je Fisherova informační matice. Označíme-li θ = (θ1 , θ2 )′ = (σ, γ)′ , pak prvky této matice lze vyjádřit vztahem   ∂2 ℓ(θ) , i, j = 1, 2. ψij = E − ∂θi ∂θj Dle [3] je −1

IE (σ, γ)

 2σ 2 −σ . = (1 + γ) −σ 1 + γ 

Skutečná hodnota (σ0 , γ0 )′ obvykle není známa, a tak se někdy matice IE (σ0 , γ0 ) odhaduje jednoduše jako IE (ˆ σ , γˆ ). Dalším možným odhadem IE (σ0 , γ0 ) je výběrová Fisherova informační matice IO (ˆ σ , γˆ ), což je matice s prvky ∂2 ψ˜ij = − ℓ(θ), i, j = 1, 2, ∂θi ∂θj 34

vyčíslená v bodě maximálně věrohodného odhadu θˆ = (ˆ σ , γˆ )′ . Označíme-li zvolený odhad matice IE (σ0 , γ0 ) jako Iˆnu a prvky matice Iˆn−1 jako υ˜ij , pak z vlastností u mnohorozměrného normálního rozdělení plyne, že σ ˆ − σ0 γˆ − γ0

≈ ≈

N (0, υ˜11 /nu ) N (0, υ˜22 /nu ).

(5.8) (5.9)

Na základě (5.8) a (5.9) je pak možno vypočítat směrodatné odchylky obou odhadů nebo určit intervaly spolehlivosti. Důležitým prvkem metody špiček nad prahem je volba prahu u, která je do značné míry subjektivní. V řadě případů, především zajímáme-li se o modelování extrémů pouze jedné veličiny, je smysluplné uvažovat práh rovný pevně zvolenému číslu. Avšak pro data, v nichž existují dlouhodobé trendy, nebo trendy způsobené závislostí na jiné veličině, je vhodnější použít proměnlivý práh. Ten zajistí, že k jeho překonávání nebude docházet častěji (nebo téměř výhradně) pro určité hodnoty regresorů, ale bude přibližně rovnoměrné v celém rozsahu dat. Při volbě prahu je dále důležité si uvědomit, že vyšší u znamená, že se parametry zobecněného Paretova rozdělení budou odhadovat na základě nižšího počtu pozorování, což povede k vyššímu rozptylu. Na druhé straně příliš nízký práh nejspíš nebude naplňovat potřebnou asymptotiku, což bude mít za následek vychýlení odhadů. Obvykle se proto volí nejnižší možný práh, který však již vede k uspokojivé aproximaci přesahů zobecněným Paretovým rozdělením. Naším cílem je odhad velmi vysokých kvantilů průtoku při daném vysokém součtu srážkových úhrnů. Jak jsme viděli v předchozích kapitolách, je závislost průtoku a součtu srážkových úhrnů pozitivní. Nebudeme proto volit práh konstantní, ale určíme ho prostřednictvím kvantilové regrese jako regresní τ -kvantil pro vhodné τ =τtr . Po odhadnutí parametrů zobecněného Paretova rozdělení se vysoký τr -kvantil průtoku (τr > τtr ) při dané hodnotě regresoru odhadne tak, že se k hodnotě prahu přičte vhodný τp -kvantil odhadnutého rozdělení. Hodnotu τp určíme ze vztahu τr − τtr τp = . 1 − τtr Neboť pro všechny hodnoty regresoru přičítáme k prahu stejné číslo, je výsledný vysoký τr -kvantil odezvy pouhým posunem zvoleného prahu. Tedy v případě, že jako práh zvolíme lineární regresní kvantil, bude odhadnutý vysoký kvantil s ním rovnoběžný. Z tohoto důvodu budeme metodu špiček nad prahem aplikovat na veličiny log(X) a log(Y ) z předchozí kapitoly, pro něž předpokládáme platnost regresního modelu s rovnoběžnými kvantily. Vysoký kvantil odezvy určený metodou špiček nad prahem a chápaný jako funkce regresoru budeme rovněž označovat regresní kvantil, neboť stejně jako regresní kvantil určený kvantilovou regresí je odhadem populačního regresního kvantilu. Z kontextu bude přitom vždy zřejmé, kterou z metod byl v daném případě vypočten. V literatuře (např. [4]) je popsáno několik grafických metod, s jejichž pomocí lze vhodný práh odhadnout. My jsme použili dvě metody, jejichž podstatou je vizuální porovnání kvality aproximace pro několik různých prahů. První metodou je vykreslení grafu, na jehož horizontální ose jsou uspořádané přesahy Z(1) ≤ . . . ≤ Z(nu ) zvoleného prahu u a na vertikální ose jsou vyneseny 35

kvantily ˆ −1 H



i nu + 1



σ ˆ = γˆ



i 1− nu + 1

−ˆγ

!

− 1 , i = 1, . . . , nu

ˆ −1 (i/(nu + 1)), odhadu zobecněného Paretova rozdělení. Obě hodnoty, Z(i) i H ˆ dobrým odhadem H, jsou odhadem i/(nu + 1)100% kvantilu rozdělení H. Je-li H pak by se body grafu měly nacházet poblíž přímky y = x. Tento kvantil-kvantilový graf jsme vykreslili pro přesahy veličiny log(Y ) příslušné všem srážkoměrným stanicím, a to pro několik různých prahů v letním i v zimním období. Nejlepšího výsledku bylo dle našeho názoru dosaženo pro práh odpovídající regresnímu 70% kvantilu. Výsledné grafy si lze prohlédnout na obrázku 5.1 a 5.2. Druhou metodou je porovnání histogramu přesahů s hustotou odhadnutého zobecněného Paretova rozdělení. I zde se nám situace zdála dobrá pro práh představovaný regresním 70% kvantilem. Na obrázcích 5.3 a 5.4 vidíme, že pro většinu srážkoměrných stanic je shoda přijatelná, a tedy přesahy zvoleného prahu lze uspokojivě aproximovat zobecněným Paretovým rozdělením. V případě stanic Krnov, Lichnov a Opava v zimním období je bohužel shoda horší. Na tomto místě je třeba přiznat, že práh tvořený regresním 70% kvantilem neodpovídá zcela podstatě metody špiček nad prahem, která je formulovaná pro přesahy daleko vyšších prahů. Avšak při práci s vyššími prahy bychom se potýkali s velkým nedostatkem pozorování. Aproximace přesahů regresního 70% kvantilu se nám zdála dobrá, a proto jsme se rozhodli u tohoto nižšího prahu zůstat, i když použití zobecněného Paretova rozdělení zde nelze odůvodnit asymptotickým chováním špiček nad prahem. Odhady γˆ a σ ˆ parametrů zobecněného Paretova rozdělení pro jednotlivé srážkoměrné stanice v letním i zimním období jsou k dispozici v tabulkách 5.1 a 5.2, a to včetně směrodatných odchylek vypočtených z výběrové Fisherovy informační matice. U odhadů γˆ menších než -0.5 nebo blízkých této hranici je třeba počítat s tím, že vektor (ˆ σ , γˆ )′ nejspíše nemá asymptotické rozdělení (5.7). V posledních dvou sloupcích tabulek 5.1 a 5.2 jsou uvedeny hodnoty 83% a 97% kvantilu odhadu zobecněného Paretova rozdělení pro každou stanici. Tyto hodnoty po přičtení k absolutnímu členu zvoleného prahu dávají absolutní členy pro 95% a 99% regresní kvantil veličiny log(Y ). Parametry 95% a 99% regresního kvantilu určené metodou špiček nad prahem jsou uvedeny v tabulkách 5.3 a 5.4, kde je navíc možné i porovnání s parametry vypočtenými kvantilovou regresí. Graficky lze regresní kvantily z obou metod porovnat na obrázku 5.5. Je patrné, že v letním období jsou si odpovídající regresní kvantily vcelku podobné, ačkoli kvantily získané metodou špiček nad prahem jsou většinou položeny trochu výše než kvantily odhadnuté kvantilovou regresí. V zimním období je rozdíl především ve sklonu regresních kvantilů. Kvantily vzniklé metodou špiček nad prahem jsou rovnoběžné s prahem, tj. s regresním 70% kvantilem z kvantilové regrese, který je v zimním období téměř vodorovný. Oproti tomu vysoké regresní kvantily vypočtené kvantilovou regresí mají spíše rostoucí trend. Doposud byl trochu opomíjen předpoklad, že data lze modelovat pomocí posloupnosti nezávislých stejně rozdělených náhodných veličin. Z podstaty věci není přirozené předpokládat, že hodnoty denních průměrných průtoků (a tedy 36

ani jejich logaritmů) jsou nezávislé. Extrémně vysoký průtok v jednom dni bude pravděpodobně následován vysokou hodnotou průtoku i v následujícím dni. My jsme sice z hodnot logaritmu průtoku použili jen některé (a to ty, které odpovídaly vysokým pětidenním součtům srážek), ale i u těchto hodnot je velmi pravděpodobné, že budou pocházet z po sobě jdoucích dní. Z toho vyplývá, že přesahy daného prahu se budou nejspíše vyskytovat ve skupinách (shlucích). Nejrozšířenějším postupem, jak se vypořádat se závislostí přesahů, je tzv. metoda vedoucí k odstranění shluků (anglicky "declustering"). Principem této metody je výběr maximálního přesahu z každého shluku a předpoklad, že tato maxima lze již považovat za nezávislá, a lze je tedy modelovat zobecněným Paretovým rozdělením. Shluk je v našem případě definován jako skupina přesahů daného prahu pocházejících z po sobě jdoucích dní a je ukončen v momentě, kdy následující pozorování již nepochází z navazujícího dne, nebo leží pod prahem. To ale v principu umožňuje, aby dva sousední shluky měly mezi sebou mezeru pouze jediný den, což narušuje předpoklad nezávislosti shlukových maxim. Proto se často uvažuje, že shluk končí až v momentě, kdy prodleva mezi dvěma po sobě jdoucími přesahy je alespoň r dní. Volbu r lze založit například na hodnotách autokorelační funkce. Hodnoty výběrové autokorelační funkce pro logaritmus průtoku řeky Opavy jsou vypsány v tabulce 5.5 a je z nich patrné, že za nezávislé by bylo možné pokládat až hodnoty log(Y ) vzdálené od sebe alespoň o 80 dní. Dle [6] však není metoda vedoucí k odstranění shluků vhodná pro data s "dlouhou pamětí". Důvodem je fakt, že vysoká hodnota parametru r vede ke zřetězení shluků, jejichž maxima by mohla být považována za nezávislá. Dochází tak ke značné ztrátě informace. Na obrázcích 5.6 a 5.7 uvádíme pro ilustraci kvantil-kvantilový graf pro přesahy regresního 70% kvantilu po aplikaci metody vedoucí k odstranění shluků s parametrem r = 80. Tato metoda pro studovaná data není příliš vhodná, a proto pro naše závěry použijeme předchozí výsledky.

5.0.2

Výpočet v R

V softwaru R se metodě špiček nad prahem věnuje zejména knihovna POT. Maximálně věrohodný odhad parametrů zobecněného Paretova rozdělení pro přesahy daného prahu získáme příkazem fitgpd(data, threshold, est = "mle", std.err.type=" ") kde data představují původní pozorování a threshold je zvolený práh. Odhad ˆ Fisherovy informační matice získáme pomocí std.err.type="expected". IE (θ) Při std.err.type="observed" bude vypočtena výběrová Fisherova informační ˆ O dalších možných způsobech odhadu parametrů se lze dočíst matice IO (θ). v příslušném manuálu [15]. Funkce declust provádějící metodu vedoucí k odstranění shluků se nachází na přiloženém CD v souboru s názvem "MojeDeclustering.R". V tomto souboru je i stručný popis použití.

37

σ ˆ Herm 0.6125 Karl 0.5743 Krno 0.6992 Lich 0.7689 Opav 0.7770 Prad 0.7787 Rejv 0.6937 Vidl 0.6549 Zary 0.7061

se 0.0449 0.0440 0.0570 0.0695 0.0717 0.0468 0.0402 0.0393 0.0509

γˆ -0.2676 -0.2910 -0.3424 -0.4505 -0.3572 -0.2832 -0.3435 -0.2368 -0.4402

se 0.0468 0.0382 0.0408 0.0585 0.0589 0.0278 0.0301 0.0307 0.0436

83% kvantil 0.8643 0.7950 0.9289 0.9386 1.0201 1.0850 0.9207 0.9477 0.8687

97% kvantil 1.3933 1.2621 1.4274 1.3552 1.5536 1.7311 1.4139 1.5600 1.2613

Tabulka 5.1: Maximálně věrohodné odhady parametrů σ a γ, jejich směrodatné odchylky (se) a kvantily odhadu příslušného zobecněného Paretova rozdělení pro letní období.

σ ˆ Herm 0.9411 Karl 0.8455 Krno 1.3467 Lich 1.0505 Opav 1.2504 Prad 0.7598 Rejv 0.8528 Vidl 0.7896 Zary 0.9367

se 0.0798 0.0698 0.0607 0.1473 0.1658 0.0508 0.0667 0.0446 0.0831

γˆ -0.4419 -0.4217 -0.9177 -0.6997 -0.7282 -0.3124 -0.4381 -0.4354 -0.4753

se 0.0463 0.0471 0.0091 0.1156 0.1066 0.0417 0.0468 0.0297 0.0500

83% kvantil 1.1564 1.0553 1.1788 1.0668 1.2447 1.0338 1.0509 0.9750 1.1219

97% kvantil 1.6774 1.5480 1.4087 1.3723 1.5836 1.6187 1.5275 1.4194 1.5986

Tabulka 5.2: Maximálně věrohodné odhady parametrů σ a γ, jejich směrodatné odchylky (se) a kvantily odhadu příslušného zobecněného Paretova rozdělení pro zimní období.

38

95% regresní kvantil POT QR L L L α ˘ β˘ α ˆ βˆL Herm -1.0918 1.1877 0.2948 0.8250 Karl -1.6489 1.4421 -1.7884 1.4811 Krno -3.6601 2.0152 -1.7988 1.5062 Lich -2.9909 1.8379 -2.1082 1.5943 Opav -4.3117 2.1899 -0.2825 1.1592 Prad -0.6767 1.1141 -0.9133 1.1664 Rejv -0.6395 1.0718 0.0014 0.9084 Vidl 0.0143 0.9209 -0.4928 1.0438 Zary -3.8526 2.0018 -1.8290 1.4711

99% regresní kvantil POT QR L L L α ˘ β˘ α ˆ βˆL -0.5628 1.1877 -1.4529 1.4048 -1.1819 1.4421 -1.2473 1.4235 -3.1616 2.0152 0.6337 0.9856 -2.5744 1.8379 -2.9561 1.9217 -3.7783 2.1899 -1.6244 1.5837 -0.0306 1.1141 -1.6451 1.4563 -0.1463 1.0718 -0.5156 1.1347 0.6265 0.9209 -0.4464 1.1436 -3.4600 2.0018 -0.9476 1.3201

Tabulka 5.3: Absolutní členy a směrnice pro 95% a 99% regresní kvantil vypočtené kvantilovou regresní (QR) a metodou špiček nad prahem (POT) pro letní období.

95% regresní kvantil POT QR Z Z Z ˘ α ˘ β α ˆ βˆZ Herm 1.8827 0.4335 0.2295 0.9209 Karl 0.9016 0.7014 -0.2093 1.0288 Krno 3.6036 0.0000 0.8189 0.8162 Lich 3.4916 0.0000 0.1074 0.9904 Opav 5.7167 -0.6098 2.1505 0.4460 Prad 3.3819 -0.1180 3.4697 -0.1461 Rejv 1.8993 0.4121 0.6231 0.7973 Vidl 2.0853 0.2876 3.1174 -0.0114 Zary 1.8845 0.4614 0.1923 0.9675

99% regresní kvantil POT QR Z Z Z ˘ α ˘ β α ˆ βˆZ 2.4037 0.4335 2.7682 0.2808 1.3944 0.7014 2.5867 0.3304 3.8335 0.0000 3.0673 0.2176 3.7971 0.0000 0.3782 0.9988 6.0556 -0.6098 2.8450 0.2758 3.9668 -0.1180 1.7721 0.5093 2.3760 0.4121 1.8567 0.5353 2.5297 0.2876 2.2742 0.3447 2.3612 0.4614 1.8904 0.5390

Tabulka 5.4: Absolutní členy a směrnice pro 95% a 99% regresní kvantil vypočtené kvantilovou regresní (QR) a metodou špiček nad prahem (POT) pro zimní období.

k 0 R(k) 1.000 k ... R(k) . . .

1 2 3 4 5 6 0.974 0.932 0.893 0.858 0.826 0.797 77 78 79 80 81 82 0.120 0.114 0.108 0.103 0.099 0.094

... ... 83 0.090

Tabulka 5.5: Hodnoty výběrové autokorelační funkce pro logaritmus průtoku.

39

1.5

1.5 0.5

1.0 Empirical

1.5

0.0

0.0

0.5

0.5

Model

Model 1.0

1.0

1.5 Model 1.0 0.5 0.0 0.0

0.0

1.0 Empirical

1.5

0.0

0.5

(b) Karlovice

1.0 Empirical

1.5

(c) Krnov

0.0

Model 1.0

0.0

0.0

0.5

0.5

0.5

Model

Model 1.0

1.0

1.5

1.5

2.0

1.5

(a) Heřmanovice

0.5

0.0

0.5

1.0 Empirical

1.5

0.0

0.5

1.5

2.0

0.0

(e) Opava

0.5

1.0 1.5 Empirical

2.0

2.5

(f) Praděd

0.0

0.5

1.0 Empirical

(g) Rejvíz

1.5

1.0 Model 0.0

0.0

0.0

0.5

0.5

0.5

Model 1.0

Model 1.0

1.5

1.5

2.0

1.5

(d) Lichnov

1.0 Empirical

0.0

0.5

1.0 1.5 Empirical

(h) Vidly

2.0

2.5

0.0

0.5

1.0 Empirical

(i) Albrechtice-Žáry

Obrázek 5.1: Kvantil-kvantilové grafy pro přesahy 70% regresního kvantilu v letním období.

40

1.5

1.4 0.5

1.0 Empirical

1.5

2.0

0.0

0.0

0.2

0.4

0.5

Model 0.6 0.8

Model 1.0

1.0

1.2

1.5

1.5 Model 1.0 0.5 0.0 0.0

0.0

1.0 Empirical

1.5

0.0

0.2

0.4

1.2

1.4

(c) Krnov

1.5

Model 1.0 0.5 0.0

0.0

0.0

0.2

0.4

0.5

Model 0.6 0.8

Model 1.0

1.0

1.5

1.2

0.6 0.8 1.0 Empirical

2.0

(b) Karlovice

1.4

(a) Heřmanovice

0.5

0.0

0.2

0.4

0.6 0.8 1.0 Empirical

1.2

1.4

0.0

1.0 Empirical

1.5

1.0 Empirical

(g) Rejvíz

1.5

1.0 1.5 Empirical

2.0

1.5 Model 1.0 0.5

0.5

0.0

0.0 0.5

0.5

(f) Praděd

Model 1.0

1.5 Model 1.0 0.5 0.0 0.0

0.0

(e) Opava

1.5

(d) Lichnov

0.5

0.0

0.5

1.0 Empirical

(h) Vidly

1.5

0.0

0.5

1.0 Empirical

1.5

(i) Albrechtice-Žáry

Obrázek 5.2: Kvantil-kvantilové grafy pro přesahy 70% regresního kvantilu v zimním období.

41

0.5

1.0

1.5

2.0

0.0

0.0

0.4

0.5

0.8

1.0

1.2

1.5

1.5 1.0 0.5 0.0 0.0

0.0

1.0

1.5

2.0

1.0

1.5

2.0

1.0

1.5

2.0

0.0 0.2 0.4 0.6 0.8 1.0 1.2

1.2 0.4 0.0 0.5

0.5

(c) Krnov

0.8

0.8 0.4 0.0 0.0

0.0

(b) Karlovice

1.2

(a) Heřmanovice

0.5

0.0

0.5

1.5

2.0

0.0

(e) Opava

0.5

1.0

1.5

2.0

(f) Praděd

0.0

0.5

1.0

1.5

(g) Rejvíz

2.0

0.8 0.4 0.0

0.0

0.0

0.5

0.5

1.0

1.0

1.2

1.5

(d) Lichnov

1.0

0.0

0.5

1.0

1.5

(h) Vidly

2.0

0.0

0.5

1.0

1.5

2.0

(i) Albrechtice-Žáry

Obrázek 5.3: Histogramy přesahů 70% regresního kvantilu v letním období s hustotou příslušného odhadu zobecněného Paretova rozdělení.

42

0.0 0.2 0.4 0.6 0.8 1.0 1.2

1.2 1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.0

0.2

0.2 0.0 0.0

0.5

1.0

1.5

2.0

0.0

1.0

1.5

2.0

0.0

0.5

(b) Karlovice

1.0

1.5

2.0

(c) Krnov

0.0

0.5

1.0

1.5

2.0

0.0

0.0

0.2

0.4

0.4

0.6

0.8

0.8

1.2

1.0

0.0 0.2 0.4 0.6 0.8 1.0 1.2

(a) Heřmanovice

0.5

0.0

0.5

1.5

2.0

0.5

1.0

1.5

(g) Rejvíz

2.0

0.5

1.0

1.5

2.0

(f) Praděd

0.0

0.2

0.4

0.6

0.8

1.0

0.0 0.2 0.4 0.6 0.8 1.0 1.2

1.0 0.8 0.6 0.4 0.2 0.0 0.0

0.0

(e) Opava

1.2

(d) Lichnov

1.0

0.0

0.5

1.0

1.5

(h) Vidly

2.0

0.0

0.5

1.0

1.5

2.0

(i) Albrechtice-Žáry

Obrázek 5.4: Histogramy přesahů 70% regresního kvantilu v zimním období s hustotou příslušného odhadu zobecněného Paretova rozdělení.

43

4 2

log(Y)

3 1

1

2

log(Y)

3

4

5

POT − 99% QR − 99%

POT − 95% QR − 95%

POT − 99% QR − 99%

0

0

POT − 95% QR − 95% 3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

log(X)

3.4

3.6

3.8

4.0

log(X)

(b) Heřmanovice, zima

2

log(Y)

3

1

2

log(Y)

3

4

4

5

(a) Heřmanovice, léto

3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

1

POT − 95% QR − 95% POT − 99% QR − 99% 4.6

3.0

3.2

log(X)

3.4

3.6

POT − 99% QR − 99% 3.8

log(X)

(d) Karlovice, zima

2

log(Y)

3

1

2

log(Y)

3

4

4

5

(c) Karlovice, léto

3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

1

POT − 95% QR − 95% POT − 99% QR − 99% 4.6

3.0

log(X)

3.2

3.4

3.6

POT − 99% QR − 99% 3.8

log(X)

(e) Krnov, léto

(f) Krnov, zima

Obrázek 5.5: Porovnání 95% a 99% regresních kvantilů vypočtených kvantilovou regresní (QR) a metodou špiček nad prahem (POT).

44

4

6

3

5

2

log(Y)

4 3

log(Y)

1

2 0 3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

1

POT − 95% QR − 95% POT − 99% QR − 99% 4.6

3.0

3.2

3.4

log(X)

3.8

log(X)

(h) Lichnov, zima

2 1

2

3

log(Y)

3

4

5

4

(g) Lichnov, léto

log(Y)

3.6

POT − 99% QR − 99%

3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

1

POT − 95% QR − 95% POT − 99% QR − 99% 4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(i) Opava, léto

(j) Opava, zima

2

log(Y) 1

1

2

log(Y)

3

3

4

4

POT − 99% QR − 99%

5

POT − 99% QR − 99%

3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

0

POT − 95% QR − 95% 4.6

3.0

log(X)

3.2

3.4

3.6

POT − 99% QR − 99% 3.8

4.0

log(X)

(k) Praděd, léto

(l) Praděd, zima

Obrázek 5.5: Porovnání 95% a 99% regresních kvantilů vypočtených kvantilovou regresní (QR) a metodou špiček nad prahem (POT).

45

4

POT − 99% QR − 99%

2

log(Y)

3 1

1

2

log(Y)

3

4

5

POT − 95% QR − 95%

0

POT − 95% QR − 95%

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

4.0

(n) Rejvíz, zima

2

1

1

2

log(Y)

3

3

4

4

POT − 99% QR − 99%

5

3.8

log(X)

(m) Rejvíz, léto

log(Y)

3.6

POT − 99% QR − 99%

0 3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

POT − 95% QR − 95% 4.6

3.0

3.2

3.4

log(X)

3.6

POT − 99% QR − 99% 3.8

4.0

log(X)

(p) Vidly, zima

2

log(Y)

3

1

2

log(Y)

3

4

4

5

(o) Vidly, léto

0 3.6

3.8

4.0

4.2

4.4

POT − 95% QR − 95%

0

1

POT − 95% QR − 95% POT − 99% QR − 99% 4.6

3.0

log(X)

3.2

3.4

3.6

POT − 99% QR − 99% 3.8

4.0

log(X)

(q) Albrechtice-Žáry, léto

(r) Albrechtice-Žáry, zima

Obrázek 5.5: Porovnání 95% a 99% regresních kvantilů vypočtených kvantilovou regresní (QR) a metodou špiček nad prahem (POT).

46

1.5

1.5 0.5

1.0 Empirical

1.5

0.0

0.0

0.5

0.5

Model

Model 1.0

1.0

1.5 Model 1.0 0.5 0.0 0.0

0.0

1.0 Empirical

1.5

0.5

(b) Karlovice

1.0 Empirical

1.5

(c) Krnov

Model 1.0

0.0

0.5

0.5 0.0

0.0

0.5

Model

1.0

Model 1.0 1.5

1.5

2.0

1.5

(a) Heřmanovice

0.5

0.5

1.0 Empirical

1.5

0.0

0.5

1.5

2.0

0.0

(e) Opava

0.5

1.0 1.5 Empirical

2.0

2.5

(f) Praděd

0.5

1.0 Empirical

(g) Rejvíz

1.5

Model 0.6 0.8 1.0 0.4 0.0

0.5

1.0 1.5 Empirical

(h) Vidly

2.0

2.5

0.0

0.0

0.0

0.2

0.5

0.5

Model 1.0

Model 1.0

1.5

1.2

1.5

2.0

1.4

(d) Lichnov

1.0 Empirical

0.0

0.5

1.0 Empirical

(i) Albrechtice-Žáry

Obrázek 5.6: Kvantil-kvantilové grafy pro přesahy 70% regresního kvantilu v letním období při použití metody vedoucí k odstranění shluků s parametrem r = 80.

47

1.5

1.4 1.2

1.5 1.0 Empirical

1.5

2.0

0.5

1.0 Empirical

1.5

0.2

0.4

0.6

(b) Karlovice

1.2

1.4

Model 1.0 0.5 0.0

0.2

0.4

0.5

Model 0.6 0.8

Model 1.0

1.0

1.5

1.2

0.8 1.0 Empirical

(c) Krnov

1.5

1.4

(a) Heřmanovice

0.0

2.0

0.5

0.0

0.2

0.4

0.5

Model 0.6 0.8

Model 1.0

1.0

1.5 Model 1.0 0.5 0.0 0.0

0.0

0.2

0.4

0.6 0.8 1.0 Empirical

1.2

1.4

0.5

1.0 Empirical

(g) Rejvíz

1.5

1.0 1.5 Empirical

2.0

1.5 Model 1.0 0.5

0.5

0.0

0.0 0.5

0.5

(f) Praděd

Model 1.0

1.5 Model 1.0 0.0

0.0

(e) Opava

0.5 0.0

1.5

1.5

(d) Lichnov

1.0 Empirical

0.5

1.0 Empirical

(h) Vidly

1.5

0.0

0.5

1.0 Empirical

1.5

(i) Albrechtice-Žáry

Obrázek 5.7: Kvantil-kvantilové grafy pro přesahy 70% regresního kvantilu v zimním období při použití metody vedoucí k odstranění shluků s parametrem r = 80.

48

6. Diskuse výsledků Nyní se pokusíme zhodnotit získané výsledky a shrnout zkušenosti s jednotlivými metodami. Abychom získali první představu o síle zkoumané závislosti a o podobě vhodného prediktoru do zamýšleného regresního modelu, provedli jsme výpočet tří typů korelačních koeficientů. Hlavním kritériem pro nás ovšem byl koeficient závislosti chvostů, který je definován jako limita podílu pravděpodobností P (U > u, V > u) a P (U > u), kde U a V jsou veličiny s rovnoměrným rozdělením. Veličina V vznikla transformací průtoku, veličina U vznikla transformací zvolené funkce srážkových úhrnů (v našem případě součtu srážkových úhrnů za předchozí dny). Vysoká míra závislosti chvostů byla zjištěna mezi průtokem a součtem srážkových úhrnů za předchozích pět dní, přičemž její hodnoty se pohybovaly kolem 0.4 pro letní období a 0.15 pro zimní období. Tyto hodnoty nejsou příliš vysoké, což lze vysvětlit tím, že ve zkoumané závislosti hrají důležitou roli další faktory (teplota, vlhkost vzduchu, atp.), jejichž hodnoty jsme neměli k dispozici. Nízká míra závislosti v zimním období je v souladu s intuitivním očekáváním, neboť v tomto období se vyskytují zpravidla srážky ve formě sněhu, které dlouhodobě zůstávají na zemském povrchu a neodtékají skrz půdu do toku řeky. Na základě koeficientu závislosti chvostů jsme se rozhodli hledat model pro kvantil průtoku v závislosti na součtu srážkových úhrnů za předchozích pět dní. Přitom jsme se omezili jen na modelování průtoku způsobeného vysokými hodnotami tohoto součtu, což má z praktického hlediska největší význam. Veličinu představující hodnoty pětidenního součtu srážkových úhrnů, které se nacházejí mezi 90% a 99.5% výběrovým kvantilem, jsme značili X, jim příslušný průtok byl značen Y . Po zvážení různých možností jsme se rozhodli pracovat s logaritmy obou veličin a použít model kvantilové funkce, v němž jsou všechny regresní kvantily lineární a rovnoběžné: Qlog Y (τ |x) = α(τ ) + β log x, τ ∈ (0, 1).

Tento model se nám zdál být i přes svou jednoduchost uspokojivý, avšak jeho kvalitu pro zpracovávaná data je obtížné objektivně ohodnotit. Pro dané τ lze parametry modelu vypočíst jakožto prvky minimalizující určitou ztrátovou funkci, což je hlavní myšlenka kvantilové regrese. Lze tak získat předpis pro libovolný regresní kvantil závisle proměnné, což velmi přispívá k vytvoření širšího obrazu o zkoumané závislosti. To je výhoda kvantilové regrese oproti klasické regresi, která slouží jen k modelování střední hodnoty. Další výhodou pak je nepřítomnost požadavku normálního rozdělení dat a menší citlivost na odlehlá pozorování (míněna pozorování s velkou hodnotou odezvy). Pro srovnání jsme pro odhad vysokých populačních regresních kvantilů odezvy použili také metodu špiček nad prahem, která modeluje chování přesahů nad daný práh zobecněným Paretovým rozdělením. Práh jsme se rozhodli zvolit závislý na výši srážek, a to jako regresní 70% kvantil vypočtený v kvantilové regresi. Požadovaný vysoký regresní kvantil průtoku pak měl směrnici shodnou s prahem, zatímco jeho absolutní člen vznikl tak, že jsme k absolutnímu členu prahu přičetli vhodný kvantil odhadnutého zobecněného Paretova rozdělení. Odhady parametrů α a β vypočtené oběma metodami pro τ = 0.95 a 0.99 byly uvedeny v tabulkách 5.3 a 5.4. Odhady směrnic β v letním období vyšly 49

vždy kladné, což lze interpretovat tak, že v tomto období je závislost průtoku na srážkách pozitivní. Na základě regresního kvantilu vypočteného pro zvolené τ jsme se snažili odhadnout hodnotu průtoku y, která nebude s pravděpodobností τ při pětidenním součtu srážkových úhrnů x překročena. Odhad tohoto průtoku je dán vztaˆ hem exp(ˆ α(τ ))xβ . Z praktického hlediska je důležitá absolutní velikost průtoku, nikoli relativní, proto se budeme zabývat pouze odhady pro letní období, kdy jsou hodnoty průtoku daleko větší. Speciálně pro τ = 0.95 a za podmínky, že hodnota x je rovna 90%, 94%, 96% a 99% výběrovému kvantilu veličiny X v dané stanici, jsou vypočtené odhady průtoku uvedeny v tabulce 6.1. Těchto hodnot tedy může dosahovat průměrný denní průtok v letním dni, bude-li součet srážek naměřených za předchozích 5 dní v dané stanici roven předepsaným hodnotám x. Tabulka 6.2 pak udává odhady τ -kvantilu průtoku pro τ =0.8, 0.9, 0.95 a 0.99 při pětidenních srážkách odpovídajících 95% výběrovému kvantilu veličiny X. Pro přehlednost označíme výběrový q-kvantil veličiny X jako x˜q , kde q ∈ (0, 1). Pro srovnání je v obou tabulkách vždy uveden i příslušný výběrový kvantil y˜τ průtoku Y vypočtený z hodnot y, které odpovídají srážkám z intervalu šířky 10 se středem v právě uvažované hodnotě x. V tabulkách 6.1 a 6.2 pozorujeme, že odhady vypočtené kvantilovou regresí více odpovídají empirickým odhadům, což je dáno tím, že jde o neparametrickou metodu, která vychází přímo z použitých dat. Kvůli tomu ale mohou být její odhady vysokých kvantilů nepřesné, neboť jsou založeny na malém počtu pozorování. Metoda špiček nad prahem odhaduje vysoké kvantily parametricky, důsledkem čehož se její odhady zdají být na pohled méně důvěryhodné (neboť se značně liší od empirických). Tato metoda však umožňuje odhad i zcela extrémně vysokých kvantilů, jaké již kvantilová regrese pro nedostatek pozorování nedokáže s uspokojivou přesností zachytit. Na druhé straně je třeba připomenout, že u metody špiček nad prahem nebyl v našem případě naplněn předpoklad nezávislosti přesahů prahu, na jejichž základě byly odhadnuty koeficienty zobecněného Paretova rozdělení. Je také nutné si uvědomit, že získané regresní kvantily jsou rovnoběžné s prahem, tj. s regresním 70% kvantilem z kvantilové regrese, a jejich spolehlivost je tak podmíněna i platností modelu s rovnoběžnými kvantily. Z tabulky 4.5 je ovšem patrné, že tento předpoklad není splněn u dat ze stanic Heřmanovice, Krnov a Opava, a tudíž zejména v těchto stanicích nelze podmíněné odhady kvantilu veličiny Y považovat za důvěryhodné. V tabulce 6.1 jsou tyto hodnoty vyznačeny kurzívou. Posledním problematickým prvkem metody špiček nad prahem je volba prahu. Rozhodneme-li se pro jiný práh, získané výsledky se mohou trochu lišit. Výsledky použitých metod nejsou zcela konzistentní, ale nelze říci, která metoda je pro studovaná data vhodnější a na které z nich by tedy měl být odhad kvantilu průtoku založen. Určitou představu o možném průtoku si však na základě tabulek 6.1 a 6.2 učinit lze. Podíváme-li se v těchto tabulkách zvlášť na každý řádek s odhadem kvantilu průtoku, zjistíme, že podle výšky tohoto odhadu lze v každém z těchto řádků rozdělit srážkoměrné stanice do dvou stabilních skupin. První skupinu tvoří stanice dávající zpravidla nižší odhad podmíněného kvantilu průtoku. Patří do ní stanice Heřmanovice, Praděd, Rejvíz a Vidly, jejichž společným znakem je také vysoká nadmořská výška (viz přiložená mapa). Druhá 50

skupina je tvořena zbývajícími stanicemi. Stanice z první skupiny se nacházejí u pramenů zkoumané řeky, a tedy nejdále od místa měření průtoku v Opavě. Na obrázku 4.3 pak vidíme, že horní regresní kvantily pro tyto stanice jsou méně strmé, než u stanic ze druhé skupiny, což má za následek nižší předpovídaný průtok při daných srážkách. Možnou interpretací může být skutečnost, že se tyto stanice nacházejí na rozhraní povodí, a tudíž ne všechny zde spadlé srážky nutně přispívají k průtoku studované řeky. Dále v tabulce 4.1 pozorujeme, že stanice z první skupiny disponují větším počtem vysokých pětidenních součtů srážek než stanice ze druhé skupiny, a to přesto, že původní data obsahovala ve všech stanicích přibližně stejný počet pozorování (viz Popis dat). Jelikož byla dolní mez pro veličinu X volena stejně pro všechny stanice, znamená to, že ve stanicích s vyšší nadmořskou výškou jsou vysoké srážkové úhrny častější. Pozitivní závislost množství srážek na nadmořské výšce je skutečně na mnoha místech obecně pozorovaným jevem. V rámci každé ze dvou vytvořených skupin jsou již odhady kvantilu průtoku vcelku podobné. Na každé z těchto skupin můžeme založit přibližný odhad kvantilu průtoku například výpočtem průměru z příslušných hodnot z tabulky 6.2. Dostaneme tak, že při pětidenním součtu srážkových úhrnů x = x˜0.95 [mm] bude s pravděpodobností 1 − τ překročen přibližně tento průtok [m3 /s]: 1−τ skupina y˜τ QR POT

0.2 1 33.96 34.71 33.49

2 64.48 68.76 72.98

0.1 1 43.54 49.46 49.10

0.05

2 81.80 91.75 105.78

51

1 53.80 61.80 67.30

2 93.12 109.56 140.77

0.01 1 103.11 113.18 114.95

2 119.19 158.03 217.97

x˜0.9 y˜0.95 QR POT x˜0.94 y˜0.95 QR POT x˜0.96 y˜0.95 QR POT x˜0.99 y˜0.95 QR POT

Herm 78.44 45.08 49.10 59.70 88.10 54.94 54.03 68.52 94.80 81.10 57.40 74.76 104.00 130.18 61.96 83.45

Karl 71.90 94.10 94.04 91.51 78.47 72.82 107.04 103.80 82.13 81.71 114.53 110.87 96.86 110.00 146.21 140.63

Krno 63.90 83.67 86.75 111.92 70.07 93.82 99.67 134.76 74.94 84.00 110.29 154.31 93.01 168.30 152.70 238.47

Lich 66.40 93.96 97.60 112.20 71.53 108.50 109.91 128.66 78.53 92.72 127.53 152.72 90.60 109.75 160.18 198.63

Opav 64.10 86.17 93.70 121.42 70.24 119.80 104.18 148.35 77.32 107.05 116.45 183.07 91.26 156.00 141.11 263.18

Prad Rejv 70.82 79.99 66.25 51.96 57.74 53.62 58.53 57.81 78.13 88.86 50.35 59.20 64.75 58.99 65.29 64.70 82.80 92.19 47.16 62.40 69.29 60.99 69.66 67.30 97.04 99.60 56.83 105.50 83.38 65.43 83.13 73.12

Vidl 73.22 50.58 53.98 52.89 83.20 45.72 61.67 59.49 87.79 45.10 65.24 62.51 100.93 81.66 75.46 71.08

Zary 69.25 75.38 81.88 102.54 79.59 105.05 100.49 135.49 86.39 112.60 113.37 159.65 96.88 125.25 134.19 200.81

Tabulka 6.1: Odhady 95% kvantilu průtoku Y [m3 /s] při srážkách odpovídajících danému kvantilu x˜ veličiny X [mm]. K dispozici je výběrový kvantil (˜ y0.95 ) a odhad na základě kvantilové regrese (QR) a metody špiček nad prahem (POT).

x˜0.95 y˜0.8 QR POT y˜0.9 QR POT y˜0.95 QR POT y˜0.99 QR POT

Herm 91.40 51.52 36.88 38.17 56.98 47.81 54.03 77.32 55.70 72.12 164.54 132.96 118.43

Karl 80.05 55.64 59.44 60.10 64.80 81.44 82.79 72.82 110.26 107.57 95.38 147.15 166.72

Krno 72.28 60.58 65.61 73.83 90.36 93.97 107.48 94.66 104.43 144.53 103.27 128.08 230.89

Lich 74.29 71.90 72.23 71.74 86.60 101.33 105.04 93.20 116.73 138.86 109.04 204.94 205.64

Opav 73.50 50.86 70.62 79.21 62.82 90.92 119.68 92.32 109.80 165.18 173.66 177.91 272.74

Prad 80.80 17.50 33.17 30.87 25.02 51.42 47.78 31.46 67.34 68.42 75.42 115.69 125.42

Rejv Vidl 90.59 86.34 36.74 30.08 33.99 34.78 34.21 30.74 51.87 40.27 48.54 50.06 49.63 44.96 61.72 44.70 60.04 64.11 66.55 62.09 87.33 85.15 99.26 104.81 105.79 110.15

Zary 82.84 83.40 75.90 80.03 104.40 91.08 113.89 112.60 106.58 147.72 114.58 132.07 213.85

Tabulka 6.2: Odhady 80%, 90%, 95% a 99% kvantilu průtoku Y [m3 /s] při srážkách odpovídajících 95% kvantilu x˜0.95 veličiny X [mm]. K dispozici je výběrový kvantil (˜ y ) a odhad na základě kvantilové regrese (QR) a metody špiček nad prahem (POT).

52

7. Závěr Cílem této práce bylo modelovat podmíněné rozdělení průtoku při daných vysokých hodnotách srážkových úhrnů. Speciálně nás zajímal odhad velikosti průtoku, který při daných srážkách nebude překročen se zvolenou pravděpodobností. K prvnímu průzkumu závislosti vysokých hodnot obou veličin jsme použili koeficient závislosti chvostů. Jeho hodnoty vyšly kladné, ovšem nikterak velké. Domníváme se, že je to důsledek vlivu dalších hydrologických faktorů, jejichž hodnoty jsme neměli k dispozici. Při hledání vhodného regresního modelu pro kvantily průtoku jsme jako vysvětlující proměnnou použili přibližně 10 % nejvyšších součtů srážkových úhrnů za předešlých pět dní. Pro odhad vysokých kvantilů odezvy jsme použili kvantilovou regresi a také metodu špiček nad prahem. Tyto metody daly požadované odhady, které mohou být použity pro predikci extremálních hodnot průtoků. V několika případech však byly tyto odhady shledány jako nevěrohodné, zpravidla proto, že nebyly naplněny všechny předpoklady k aplikaci metod potřebné. Výsledky a vyhodnocení použití metod jsou podrobně diskutovány v předchozí kapitole. Zdůrazněme, že provedené analýzy jsou pouze prvotním průzkumem dané závislosti a nabízí se řada možností dalšího výzkumu. Mnoho variant vzniká kolem volby modelu. Jako prediktor by bylo možné zvolit jinou formu srážkových úhrnů, než jen prostý součet. Smysl by mělo třeba použití součtu váženého, neboť není důvod přisuzovat všem předchozím dnům stejnou důležitost. Podstatným prvkem v našem modelu byla také počáteční volba transformace, tj. logaritmu obou proměnných. Existuje však řada dalších použitelných transformací, například třída Coxových transformací. Kromě toho by mohlo být rozumné zaměřit se na hledání složitějšího tvaru modelu. Nabízí se třeba použití autoregresního modelu, v němž by ovšem díky silné závislosti průtoků z po sobě jdoucích dní mohlo dojít k vytěsnění vlivu srážkových úhrnů. Zajímavá je také myšlenka modelu, který by zahrnoval všechny srážkoměrné stanice současně, například seskupené podle nějakého kritéria (nadmořské výšky, apod.). V neposlední řadě je třeba opět připomenout, že velký vliv na průtok mají i další přírodní faktory, které nebyly obsaženy v datech, a nemohly být tudíž zahrnuty do modelu. Pokud by se v budoucnu podařilo získat záznamy alespoň o některých z těchto faktorů, velmi by to přispělo ke kvalitě sestavovaného modelu. Na závěr je třeba říci, že metody popsané v této práci mohou být užitečné v obecných situacích, kdy je zkoumána závislost extrémně vysokých hodnot několika veličin.

53

Seznam použité literatury [1] Anděl, J.: Základy matematické statistiky. Praha: Matfyzpress, 2007, ISBN 80-7378-001-1. [2] Bassett, G.; Koenker, R.: Asymptotic Theory of Least Absolute Error Regression. Journal of the American Statistical Association, ročník 73, č. 363, September 1978: s. 618–622. [3] Beirlant, J.; Goegebeur, Y.; Segers, J.; aj.: Statistics of Extremes, Theory and Applications. Wiley series in probability and statistics, John Wiley & Sons, 2004, ISBN 978-0-471-97647-9. [4] Coles, S.: An introduction to statistical modeling of extreme values. Springer series in statistics, Springer, 2001, ISBN 1852334592. [5] Coles, S.; Heffernan, J.; Tawn, J.: Dependence Measures for Extreme Value Analyses. Extremes, ročník 2, č. 4, 2000: s. 339–365. [6] Jarušková, D.; Hanek, M.: Peaks over threshold method in comparison with block-maxima method for estimating high return levels of several northern Moravia precipitation and discharges series. Journal of Hydrology and Hydromechanics, ročník 54, č. 4, 2006: s. 309–319. [7] Kašpárek, L.; Novický, O.: Background information BILAN, [online] [cit. 31.1.2012]. URL http://www.geo.uio.no/edc/software/BILAN/Background_ Information_BILAN.pdf [8] Kendall, M. G.: Rank correlation methods. Charles Griffin, 1962. [9] Koenker, R.: A note on L-estimates for linear models. Statistics & Probability Letters, ročník 2, č. 6, December 1984: s. 323–325. [10] Koenker, R.: Confidence Intervals for Regression Quantiles. In Asymptotic statistics: proceedings of the fifth Prague symposium, held from September 49, 1993, Contributions to statistics, New York: Physica-Verlag, 1994, ISBN 9783790807707, s. 349–359. [11] Koenker, R.: Quantile regression. Cambridge New York: Cambridge University Press, 2005, ISBN 0-521-60827-9. [12] Koenker, R.: Package "quantreg" reference manual. September 2011, [online] [cit. 28.12.2011]. URL http://www.r-project.org [13] Koenker, R.: personal communication, February 2012. [14] Koenker, R. W.; Bassett, J., Gilbert: Regression Quantiles. Econometrica, ročník 46, č. 1, January 1978: s. 33–50.

54

[15] Ribatet, M.: Package "POT" reference manual. October 2007, [online] [cit. 15.2.2012]. URL http://www.r-project.org [16] Rüschendorf, L.: On the distributional transform, Sklar’s Theorem, and the empirical copula process. Journal of Statistical Planning and Inference, , č. 139, 2009: s. 3921–3927.

55

Přílohy 7.1

Koeficient závislosti chvostů

Odhady funkce χ(u) pro průtok a součet srážkových úhrnů za posledních d dní vykreslené pro jednotlivé stanice. Koeficient χ závislosti chvostů je definován jako limita limu→1 χ(u), tedy nás zajímá hlavně chování funkce χ(u) ˆ v levostranném okolí bodu 1. V těsné blízkosti tohoto bodu je však výpočet χ(u) ˆ poznamenán značným úbytkem pozorování. Následující grafy se proto zaměřují na oblast před bodem 0.99.

Letní období

d=7 d=8 d=9

Chi(u)

d=3 d=4 d=5 d=6

0.90

0.92

0.94

0.96

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

Chi(u)

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

7.1.1

0.98

d=3 d=4 d=5 d=6

0.90

d=7 d=8 d=9

0.92

0.94

u

(b) Karlovice

Chi(u)

d=7 d=8 d=9

0.92

0.94

0.96

0.98

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

Chi(u)

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.90

0.98

u

(a) Heřmanovice

d=3 d=4 d=5 d=6

0.96

d=3 d=4 d=5 d=6

0.90

u

d=7 d=8 d=9

0.92

0.94

0.96 u

(c) Krnov

(d) Lichnov

Obrázek 7.1: Koeficient závislosti chvostů pro letní období.

56

0.98

Chi(u)

d=7 d=8 d=9

0.90

0.92

0.94

0.96

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

Chi(u)

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

d=3 d=4 d=5 d=6

0.98

d=3 d=4 d=5 d=6

0.90

d=7 d=8 d=9

0.92

0.94

u

0.94

0.96

0.98

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

Chi(u) 0.92

d=3 d=4 d=5 d=6

0.90

0.96

0.98

d=7 d=8 d=9

0.92

u

0.94 u

(g) Rejvíz

(h) Vidly

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

0.90

0.98

(f) Praděd

d=7 d=8 d=9

Chi(u)

Chi(u)

0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60

(e) Opava

d=3 d=4 d=5 d=6

0.96 u

d=3 d=4 d=5 d=6

0.90

d=7 d=8 d=9

0.92

0.94

0.96

0.98

u

(i) Albrechtice-Žáry

Obrázek 7.1: Koeficient závislosti chvostů pro letní období. 57

0.20 0.15 0.10

Chi(u)

0.05

0.10 0.05

Chi(u)

0.15

0.20

0.25

Zimní období

0.25

7.1.2

0.90

d=3 d=4 d=5 d=6

0.00

d=7 d=8 d=9

−0.05

−0.05

0.00

d=3 d=4 d=5 d=6

0.92

0.94

0.96

0.98

0.90

d=7 d=8 d=9

0.92

0.94

u

0.96

0.98

0.96

0.98

0.25 0.20 0.15 0.10 0.05

0.10

Chi(u)

0.15

0.20

0.25

(b) Karlovice

0.05

Chi(u)

0.98

u

(a) Heřmanovice

d=3 d=4 d=5 d=6

0.00

d=7 d=8 d=9

−0.05

−0.05

0.00

d=3 d=4 d=5 d=6 0.90

0.92

0.94

0.96

0.98

0.90

d=7 d=8 d=9

0.92

0.94

u

u

d=7 d=8 d=9

0.10

0.15

0.20

d=3 d=4 d=5 d=6

0.05

0.05

0.10

Chi(u)

0.15

0.20

0.25

(d) Lichnov

0.25

(c) Krnov

0.90

0.00

d=7 d=8 d=9

−0.05

0.00

d=3 d=4 d=5 d=6

−0.05

Chi(u)

0.96

0.92

0.94

0.96

0.98

0.90

u

0.92

0.94 u

(e) Opava

(f) Praděd

Obrázek 7.2: Koeficient závislosti chvostů pro zimní období. 58

0.25 0.05

0.10

Chi(u)

0.15

0.20

0.25 0.20 0.15

0.00

d=3 d=4 d=5 d=6

−0.05 0.92

0.94

0.96

0.98

0.90

d=7 d=8 d=9

0.92

0.94

u

0.96

u

(h) Vidly

0.05

0.10

0.15

0.20

0.25

(g) Rejvíz

0.00

d=3 d=4 d=5 d=6

−0.05

0.90

d=7 d=8 d=9

Chi(u)

Chi(u)

0.10 0.05 −0.05

0.00

d=3 d=4 d=5 d=6

0.90

d=7 d=8 d=9

0.92

0.94

0.96

0.98

u

(i) Albrechtice-Žáry

Obrázek 7.2: Koeficient závislosti chvostů pro zimní období.

59

0.98

Rovnoběžné regresní kvantily

log(Y)

2

3 0

1

1

2

log(Y)

3

4

5

4

7.2

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

4.0

log(X)

(b) Heřmanovice, zima

log(Y)

2

3 1

1

2

log(Y)

3

4

5

4

(a) Heřmanovice, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(d) Karlovice, zima

3 2

log(Y)

3 1

1

2

log(Y)

4

5

4

(c) Karlovice, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6

3.8

log(X)

(e) Krnov, léto

(f) Krnov, zima

Obrázek 7.3: Regresní kvantily (20%, 35%, 50%, 65% a 80%) počítané zvlášť (čárkovaně) a se společnou směrnicí (plnou čarou).

60

4

6

3

5

2

log(Y)

4 3

log(Y)

1

2 1 0

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(h) Lichnov, zima

2

log(Y)

3 1

1

2

log(Y)

4

3

5

4

(g) Lichnov, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

log(X)

(j) Opava, zima

log(Y)

2

0

1

1

2

log(Y)

3

3

4

5

4

(i) Opava, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6

3.8

4.0

log(X)

(k) Praděd, léto

(l) Praděd, zima

Obrázek 7.3: Regresní kvantily (20%, 35%, 50%, 65% a 80%) počítané zvlášť (čárkovaně) a se společnou směrnicí (plnou čarou).

61

4 5

log(Y)

2

3

4 3

log(Y)

1

2 1

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6

3.8

4.0

3.8

4.0

3.8

4.0

log(X)

(n) Rejvíz, zima

2.0

log(Y) 0

0.5

1

1.0

1.5

2

log(Y)

3

2.5

4

3.0

3.5

5

(m) Rejvíz, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

3.2

3.4

log(X)

3.6 log(X)

(p) Vidly, zima

log(Y)

2

0

1

1

2

log(Y)

3

3

4

5

4

(o) Vidly, léto

3.6

3.8

4.0

4.2

4.4

4.6

3.0

log(X)

3.2

3.4

3.6 log(X)

(q) Albrechtice-Žáry, léto

(r) Albrechtice-Žáry, zima

Obrázek 7.3: Regresní kvantily (20%, 35%, 50%, 65% a 80%) počítané zvlášť (čárkovaně) a se společnou směrnicí (plnou čarou).

62